The expression of the supercurrent in a superconductor is $$\vec{j}_s=-\frac{n_se^2}{m}\vec{A}$$ where $\vec{A}$ is the vector potential, $n_s$ is the number density of superconducting carriers and $e,m$ are the charge and mass of the electron. Wikipedia article of London equations notes that this equation suffers from the disadvantage that in this form $\vec{j}_s$ does not seem to be gauge invariant. However, it asserts that this expression is true only in the Coulomb gauge (${\rm div}~\vec{A}=0$). I want to show that this is true only in the Coulomb gauge.
I started from the general expression of the supercurrent $$\vec{j}_s=\frac{-e}{2m}\Big\{\psi^*\Big(-i\hbar\vec{\nabla}-q\vec{A}\Big)\psi+\psi\Big(-i\hbar\vec{\nabla}-q\vec{A}\Big)^*\psi^*\Big\}\\=\frac{ie\hbar}{2m}(\psi^*\vec{\nabla}\psi-\psi\vec{\nabla}\psi^*)-\frac{2e^2}{m}\vec{A}|\psi|^2$$ where $q=-2e$ has been used. Now, assuming that the macroscopic wavefunction has the form $$\psi(\vec{r})=\rho^{1/2}\exp[i\theta(\vec{r})]$$ with a spatially uniform modulus $\sqrt{\rho}$. With direct subbstitution, $\vec{j}_s$ simplifies to [Ref. Aschroft & Mermin, Eqn. $(34.29)$] $$\vec{j}_s=-\Big[\frac{e\hbar}{m}\vec{\nabla}\theta+\frac{2e^2}{m}\vec{A}\Big]\rho.$$
- The last expression for $\vec{j}_s$ is gauge invariant but the first expression for $\vec{j}_s$ is not. Please explain how the first expression for $\vec{j}_s$ arises when one chooses the gauge ${\rm div}~\vec{A}=0$.