Note by the way that the result extends to any number of dimensions. Say we're working in an $n$-dimensional (oriented, Riemannian) manifold, and (in terms of positively oriented coordinates $(x^1,\dots, x^n)$) we have a vector field $A=A^i\frac{\partial}{\partial x^i}$ and the corresponding $1$-form $\alpha=A_i\,dx^i$, and we wish to calculate $\text{div}(A):=\star d\star\alpha$. We have in general that for any multindex $I=(i_1,\dots, i_k)$,
\begin{align}
\star(dx^I)&=\sqrt{|g|}\sum_{|J|=n-k}\text{sgn}(J',J)\det\left(g^{ij'}\right)_{\substack{i\in I\\j'\in J'}}\,dx^J,
\end{align}
where the sum is over all increasing tuples $J=(j_1,\dots, j_{n-k})$ and $J'=\{1,\dots, n\}\setminus J$ is the remaining tuples arranged in increasing order. Also, $\text{sgn}$ denotes the sign of the permutation, so you may write it as $\epsilon_{J'J}$ (the Levi-Civita symbol, which is just another name for the sign of a permutation). In the special case that $k=1$ and that $I=(i)$ consists of a single index, the formula above can be specialized because the sum over all indices $|J|=n-1$ can be replaced by a sum over all the complementary indices $j=1,\dots, n$:
\begin{align}
\star(dx^i)&=\sqrt{|g|}\sum_{j=1}^n(-1)^{j-1}\cdot g^{ij} \,dx^1\wedge\cdots\wedge
\widehat{dx^j}\wedge \cdots \wedge dx^n,
\end{align}
where the $\widehat{dx^j}$ indicates that $dx^j$ is removed. Therefore,
\begin{align}
\star\alpha&=\sum_{j=1}^n(-1)^{j-1}\left(A^j\sqrt{|g|}\right)\,dx^1\wedge\cdots\wedge
\widehat{dx^j}\wedge \cdots \wedge dx^n,
\end{align}
where we have used $A_ig^{ij}=A^j$. Therefore,
\begin{align}
d\star\alpha&=\sum_{j=1}^n(-1)^{j-1}d\left(A^j\sqrt{|g|}\right)\wedge \,dx^1\wedge\cdots\wedge
\widehat{dx^j}\wedge \cdots \wedge dx^n\\
&=\sum_{j=1}^n(-1)^{j-1}\frac{\partial(A^j\sqrt{|g|})}{\partial x^j}\,dx^j
\wedge dx^1\wedge\cdots\wedge
\widehat{dx^j}\wedge \cdots \wedge dx^n\\
&=\sum_{j=1}^n\frac{\partial(A^j\sqrt{|g|})}{\partial x^j}\,dx^1\wedge \cdots\wedge dx^n,
\end{align}
of course at this stage we can drop the summation signs. Finally, we can calculate the Hodge star of this expression to finally arrive at
\begin{align}
\text{div}(A)&:=\star d\star \alpha = \frac{\partial(A^k\sqrt{|g|})}{\partial x^k}\cdot\star{dx^1\wedge \cdots \wedge dx^n}=
\frac{\partial(A^k\sqrt{|g|})}{\partial x^k}\cdot \frac{1}{\sqrt{|g|}},
\end{align}
which is the desired expression. Note that in the pseudo-Riemannian case, this will all go through, except for some annoying minus signs that (depending on your convention for the definition of the Hodge star) we have $\star \mu = (-1)^{\#}$, where $\mu$ is the volume form (in local positively oriented coordinates $\sqrt{|g|}\,dx^1\wedge \cdots\wedge dx^n$) and $\#$ denotes the number of minus signs in the diagonalized form of the metric (so $(-1)^{\#}=-1$ if you're working with a metric of signature $(-,+,\cdots ,+)$), so one should probably define $\text{div}(A)= (-1)^{\#}\star d\star \alpha$.
