Why doesn't homoskedacticity bias an estimator? I keep reading that homoskedasticity biases the SE, but not the estimator. Why? I'm imagining a plot, where a bunch of errors are clustered on the top left. That would "pull" the OLS line up  towards there, which should affect the estimator right?
 A: What heteroskedasticity describes is that the variation of the errors may depend on the values of the regressors. That is, that for certain values of $x$ we expect that, while we still expect zero errors on average, any given error tends to be further away from the true regression line in either direction.
The situation you describe rather concerns the situation in which the errors systematically deviate from the regression line in one direction (or in one direction for some range of $x$, and in another for another range of $x$), so that errors would no longer have mean zero for such predictor values, for example due to omitted nonlinearities or omitted variables.
Here is an example in which the error term $u$ of the model is generated such that it correlates with the regressor $X$ (see code below). This causes the scatter plot not to scatter around the true (red) regression line, such that, despite the huge sample size of $n=10,000$, the (blue) estimated OLS line is pretty far away from the true value $\beta_1=0.5$.

library(mvtnorm)
# truth 
beta0 <- 1
beta1 <- 0.5

# generate some data with correlation between X and u
n <- 10000 
errors <- rmvnorm(n, mean = rep(0, 2), sigma = matrix(c(1,-0.5,-0.5,1),2,2))
u <- errors[,1]
X <- errors[,2]
y <- beta0 + beta1*X + u 

plot(X,y,xlab="x",ylab="y")
abline(a = beta0, b = beta1, col="red", lwd=4) # the truth
regr <- lm(y~X) 
abline(regr, col="blue", lwd=3)

A: Consider a heteroskedastic linear regression with model form:
$$\boldsymbol{Y} = \boldsymbol{x} \boldsymbol{\beta} + \boldsymbol{\varepsilon} \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \boldsymbol{\varepsilon}  \sim \text{N}(\boldsymbol{0}, \sigma^2 \text{diag}(\boldsymbol{\tau} ) ),$$
where $\boldsymbol{\tau} = (\tau_1, ..., \tau_n) \in \mathbb{R}^n$ gives the underlying variance scales for the variables.  (To simplify our analysis we will also assume that none of these weightings is zero, which means that without loss of generality, we can let $|\text{diag}(\boldsymbol{\tau})| = 1$.)  The ordinary least squares (OLS) estimator is:
$$\hat{\boldsymbol{\beta}} = (\boldsymbol{x^\text{T}} \boldsymbol{x})^{-1} (\boldsymbol{x^\text{T}} \boldsymbol{Y}).$$
Since $\mathbb{E}(\boldsymbol{\varepsilon}) = \boldsymbol{0}$, this has expected value:
$$\begin{equation} \begin{aligned}
\mathbb{E}( \hat{\boldsymbol{\beta}} ) = \mathbb{E}( (\boldsymbol{x^\text{T}} \boldsymbol{x})^{-1} (\boldsymbol{x^\text{T}} \boldsymbol{Y}) ) &= (\boldsymbol{x^\text{T}} \boldsymbol{x})^{-1} \boldsymbol{x^\text{T}} \mathbb{E}(\boldsymbol{Y})\\[6pt]
&= (\boldsymbol{x^\text{T}} \boldsymbol{x})^{-1} \boldsymbol{x^\text{T}} \mathbb{E}(\boldsymbol{x} \boldsymbol{\beta} + \boldsymbol{\varepsilon} )\\[6pt]
&= \boldsymbol{\beta} + (\boldsymbol{x^\text{T}} \boldsymbol{x})^{-1} \boldsymbol{x^\text{T}} \mathbb{E} (\boldsymbol{\varepsilon}) \\[6pt]
&= \boldsymbol{\beta}.
\end{aligned} \end{equation}$$
It also has variance:
$$\begin{equation} \begin{aligned}
\mathbb{V}( \hat{\boldsymbol{\beta}} ) = \mathbb{V}( (\boldsymbol{x^\text{T}} \boldsymbol{x})^{-1} (\boldsymbol{x^\text{T}} \boldsymbol{Y}) ) &= (\boldsymbol{x^\text{T}} \boldsymbol{x})^{-1} \boldsymbol{x^\text{T}} \mathbb{V}(\boldsymbol{Y}) \boldsymbol{x}  (\boldsymbol{x^\text{T}} \boldsymbol{x})^{-1} \\[6pt]
&= (\boldsymbol{x^\text{T}} \boldsymbol{x})^{-1} \boldsymbol{x^\text{T}} \mathbb{V}(\boldsymbol{\varepsilon}) \boldsymbol{x}  (\boldsymbol{x^\text{T}} \boldsymbol{x})^{-1} \\[6pt]
&= \sigma^2 (\boldsymbol{x^\text{T}} \boldsymbol{x})^{-1} (\boldsymbol{x^\text{T}} \text{diag}(\boldsymbol{\tau}) \boldsymbol{x})  (\boldsymbol{x^\text{T}} \boldsymbol{x})^{-1}. \\[6pt]
\end{aligned} \end{equation}$$
As you can see, nothing in the expected value derivation depends on the variance of the error vector in the model.  Heteroscedasticity affects the variance matrix of the OLS coefficient estimator, but it does not affect the expected value.  Intuitively, this is because positive and negative variation in the (symmetric) distribution of the error terms "balance out" in the expected value.
A: I would like to go a bit theoretical on this. One of the most important assumptions for the unbiasedness of estimators is:
\begin{equation}
E(u|x) = 0
\end{equation}
This assumption implies that we expect no variance between a parameter's true value and its estimator. This can also be expressed as:
\begin{equation}
E(\hat{\mu}) = \mu
\end{equation}
with $\hat{\mu}$ is the estimator of parameter $\mu$. This is called the zero conditional mean assumption, and it is distinct from the homoskedasticity assumption, which is:
\begin{equation}
Var(u|x) = {\sigma}^2
\end{equation}
Woolridge (2002) clearly distinguishes between the assumptions of unbiasedness and homoskedasticity. The first concerns the expected value of $u$, while the second concerns the variance of $u$. @Christopher Hanck illustrates a good example of heteroskedasticity: the variation of $\hat{\beta}$ should be spread out unequally along the true $\beta line in every direction.
In your case, when you imagine a bunch of errors is clustered on the top left, such data will violate the unbiasedness assumption as the $E(u|x) > 0$. Yet, it does not necessarily means that $Var(u|x)$ is non-constant.
