Why does the covariant derivative not affect the vacuum expectation value of the Higgs field? Consider the following Lagrangian density:
$$
\mathcal{L}_{1} = (\partial_{\mu}\phi^{*})(\partial^{\mu}\phi)-m^{2}\phi^{*}\phi-h(\phi^{*}\phi)^{2}
$$
for the complex scalar field $\phi$ with $h>0$. I obtain the vacuum expectation value as the minimum of the Hamiltonian density
\begin{align}
\mathcal{H}_{1} &= \pi_{\phi}\dot{\phi}+\pi_{\phi^{*}}\dot{\phi^{*}} - \mathcal{L}_{1} \\ &=
\frac{1}{c^{2}}|\dot{\phi}|^{2}+|\nabla\phi|^{2}+m^{2}|\phi|^{2}+h|\phi|^{4}
\end{align}
with respect to $|\phi|$. Here, the canonical momenta are
\begin{align}
\pi_{\phi} &= \frac{\partial\mathcal{L}_{1}}{\partial\dot{\phi}}=\frac{1}{c^{2}}\dot{\phi^{*}},\\
\pi_{\phi^{*}} &= \frac{\partial\mathcal{L}_{1}}{\partial\dot{\phi^{*}}}=\frac{1}{c^{2}}\dot{\phi}.
\end{align} For the symmetry-broken case $m^{2}<0$, I find
$$
|\phi_{0}| = \sqrt{\frac{-m^{2}}{2h}}.
$$
Replacing the partial derivative by the covariant derivative
$$
\partial_{\mu}\phi\mapsto D_{\mu}\phi = (\partial_{\mu}-igA_{\mu})\phi.
$$
and adding the free part of the vector field $A_{\mu}$,
$$
-\frac{1}{4}F_{\mu\nu}F^{\mu\nu},
$$
the Lagrangian density is
$$
\mathcal{L}_{2} = [(\partial_{\mu}+igA_{\mu})\phi^{*}][(\partial^{\mu}-igA^{\mu})\phi]-m^{2}\phi^{*}\phi-h(\phi^{*}\phi)^{2}-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}
$$
and I obtain a similar Hamiltonian density
\begin{align}
\mathcal{H}_{2} = &\pi_{\phi}\dot{\phi}+\pi_{\phi^{*}}\dot{\phi^{*}}+(\pi_{A})_{\gamma}\dot{A^{\gamma}} - \mathcal{L}_{2} \\
 = &\frac{1}{c^{2}}|\dot{\phi}|^{2}+|\nabla\phi|^{2}+m^{2}|\phi|^{2}+h|\phi|^{4}\\ &+ 2gA^{n}\text{Im}[\phi^{*}(\partial_{n}\phi)]-g^{2}A_{\mu}A^{\mu}|\phi|^{2}-(\partial_{n}A_{\mu})F^{n\mu},
\end{align}
where $\mu\in\{0,1,2,3\}$ and $n\in\{1,2,3\}$. Here, the canonical momenta are
\begin{align}
\pi_{\phi} &= \frac{\partial\mathcal{L}_{2}}{\partial\dot{\phi}}=\frac{1}{c}\left[\frac{1}{c}\dot{\phi^{*}}+ig\phi^{*}A^{0}\right],\\
\pi_{\phi^{*}} &= \frac{\partial\mathcal{L}_{2}}{\partial\dot{\phi^{*}}}=\frac{1}{c}\left[\frac{1}{c}\dot{\phi}-ig\phi A^{0}\right],\\
(\pi_{A})_{\gamma} &=\frac{\partial \mathcal{L}_{2}}{\partial (\dot{A^{\gamma}})} =-\frac{1}{c}g_{\mu\gamma}F^{0\mu}.
\end{align}
If I require the vacuum expectation value to be real (similar to the unitary gauge) and neglect $(\partial_{n}A_{\mu})F^{n\mu}$ as an offset, there is still the contribution $g^{2}A_{\mu}A^{\mu}|\phi|^{2}$ modifying $m^{2}$. This changes the vacuum expectation value in the symmetry-broken case $m^{2}<0$ and makes it dependent on the gauge field. Several books ("Fields, symmetries and quarks" by Mosel, chap. 9.1; "Symmetries in Physics" by Ludwig, Falter, chap. 14.4.3; "Group theory in physics, Vol. II" by Cornwell, chap. 19.4) claim, that the vacuum expectation value does not change compared to the unmodified Lagrange density. It also should not depend on $A_{\mu}$, since otherwise the gauge transformation to remove the massless Goldstone boson would affect the vacuum expectation value.
Question: What is the argument for the vacuum expectation value not to depend on the gauge field $A_{\mu}$, starting from the Hamiltonian density $\mathcal{H}_{2}$? I am interested in a more technical solution.
Disclaimer: I already found:


*

*Determination of the ground state of a field theory

*What is the role of the vacuum expectation value in symmetry breaking and the generation of mass?
but both address $\mathcal{L}_{1}$.
 A: Good question; first let's recap the setup. Finding the vevs of fields in the ground state is really a quantum field theory question. However, in introductory textbooks we ignore quantum effects, hoping that they'll be small, and treat every field as classical. This is reasonable because if you do have a nonzero classical vev, the quantum effects should come in a power series in $\hbar$ and hence be "small". The power series is really in $\hbar g^2$, so this is actually a weak coupling approximation.
In any case, this means that in a theory with $n$ fields, we can think of each field as just an ordinary function on space. Next, gradient terms in the Hamiltonian mean that to minimize the energy, the fields must be constant. So the problem boils down to minimizing a simple function of $n$ variables.
However, the case of gauge theories is more subtle. First, there are different classical configurations of the fields that correspond to the exact same physical state. For example, a configuration where $A_\mu = 0$ is gauge equivalent to a configuration $A_\mu = \partial_\mu \Lambda$ for any function $\Lambda$, so they must have the same energy, even though the latter can have rapid spatial variation. Second, we see there are terms in $\mathcal{H}$ that are linear in $\partial_\mu \phi$, so the minimum energy configuration need not have constant $\phi$ either.
These issues are related, and make it very hard to find the classical minimum of the potential without fixing the gauge. There is a classical minimum where $\phi$ is just constant and equal to what it would be without the gauge field, but it's gauge equivalent to all kinds of complicated configurations where both $\phi$ and $A_\mu$ vary.
To get out of this mess, we need to gauge fix. We have the gauge transformations
$$A_\mu \to A_\mu + \partial_\mu \Lambda, \quad \phi \to e^{i g \Lambda} \phi.$$
Therefore, we may choose $\Lambda$ so that $\phi$ is real at time $t = 0$. This makes the problematic $A^i \text{Im}(\phi^* \partial_i \phi)$ term vanish. Since all gradient terms are now quadratic, the vevs must be constant, so we have
$$\mathcal{H} = m^2 |\phi|^2 + h |\phi|^4 - g^2 A_\mu A^\mu |\phi|^2.$$
Now we still have a problem because $A_\mu A^\mu$ has indefinite sign, so this potential doesn't seem bounded below. But we still haven't used up all the gauge freedom, since we haven't specified the temporal dependence of $\Lambda$. So we still have the freedom to transform
$$A_0 \to A_0 + \partial_0 \Lambda$$
which we may use to set $A_0 = 0$ at time $t = 0$. Then we have
$$\mathcal{H} = m^2 |\phi|^2 + h |\phi|^4 - g^2 A_i A^i |\phi|^2$$
which is indeed positive definite thanks to the $(+---)$ signature. Minimizing this function is now straightforward; we must have $A_i = 0$, and finding the value of $|\phi|$ proceeds as in the uncharged case. Of course the values of the vevs will be gauge-dependent, but all values of gauge-invariant observables we derive will not be.
A: VEV actually represents a solution to the equations of motion. Just find the equations of motion and set all partial derivatives to zero, because we want to study the vacuum behaviors. 
For the 4-potential, its equation of motion is $\propto A^\mu|\phi|^2=0$, so its VEV is zero. That is why Higgs' VEV does not change. 
