Heteroscedasticity that depends on the regression parameters Consider a vector of observations $\mathbf{Y}$ that can be modeled as
\begin{equation}
  \mathbf{Y} \sim \mathcal{N}( \mathbf{H}\boldsymbol{\beta} , \boldsymbol{\Sigma} )
\end{equation}
with $\mathbf{H}$ a linear parametric model, and $\boldsymbol{\Sigma}$ a diagonal matrix whose elements are given by $a_1 \mathbf{H}\boldsymbol{\beta} + a_2$, and where the constants $a_1$ and $a_2$ are known.
To properly estimate $\boldsymbol{\beta}$, and since the variance is heteroscedastic but depends on the regression parameters, a first estimation can be obtained applying OLS:
$$
  \widetilde{\boldsymbol{\beta}}^{\mathrm{ols}} = \left( \mathbf{H}^\mathrm{T} \mathbf{H} \right)^{-1} \mathbf{H}^\mathrm{T} \mathbf{Y} .
$$
Then, a rough estimation of the noise variance is obtained from $\widetilde{\boldsymbol{\beta}}^{\mathrm{ols}}$ by:
$$
  \widetilde{\boldsymbol{\Sigma}}^{\mathrm{ols}} = \mathbf{I}\times (a_1 \mathbf{H}\, \widetilde{\boldsymbol{\beta}}^{\mathrm{ols}} + a_2 ),
$$
where $\mathbf{I}$ denotes the identity matrix. This rough estimation of the covariance is thus reused to update the estimation of the expectation using WLS:
\begin{equation}
\label{eq:estimated_expectation_detect}
\begin{cases}
  \widetilde{\boldsymbol{\beta}} = \left( \mathbf{H}^\mathrm{T} {\widetilde{\boldsymbol{\Sigma}}^{{\mathrm{ols}}^{-1}}} \mathbf{H} \right)^{-1} \mathbf{H}^\mathrm{T} {\widetilde{\boldsymbol{\Sigma}}^{{\mathrm{ols}}^{-1}}} \mathbf{Y} , \\
  \widetilde{\boldsymbol{\Sigma}} = \mathbf{I}\times (a_1 \mathbf{H}\, \widetilde{\boldsymbol{\beta}} + a_2 ).
\end{cases}
\end{equation}
We can prove that the estimation $\widetilde{\boldsymbol{\beta}}^{\mathrm{ols}}$ is unbiased. However, is the estimation $\widetilde{\boldsymbol{\beta}}$ unbiased?
Another question is, is there a better way to estimate $\boldsymbol{\beta}$?
 A: What you have here is a heteroscedastic linear model, which gives a non-linear function of the coefficients in the log-likelihood.  This occurs because your coefficients affect both the slope of the regression line and also the variance matrix, giving a complicated form in the  log-likelihood.  If you were to keep iterating your estimators in the way you have specified (though you have only done one iteration here), you would be doing iteratively reweighted least-squares (IRLS) estimation, which is a standard method of parameter estimation for this kind of model.
The use of IRLS estimation in heteroscedastic linear models is analysed in Chen and Shao (1993).  In that paper, the authors show that, under certain broad conditions, the IRLS estimator converges in distribution to normality, around the true coefficients (assuming the model is correctly specified).  The estimator is not generally unbiased, which is standard whenever you have non-linearity in a regression.  However, under broad conditions, it is asymptotically unbiased.  If you would like to know the details, I recommend reading up on IRLS estimation, and the linked paper.
