Consistency of Coulomb gauge

Here’s a subtlety that I only noticed embarrassingly late: Let’s say we’re doing magnetostatics, and trying to find $\mathbf B$ given a current density $\mathbf J$, with some BCs involving $\mathbf B$. We can write $\mathbf B = \nabla \times \mathbf A$ for some vector potential $\mathbf A$ (guaranteed by the Helmholtz decomposition). Then the 4th Maxwell equation becomes $$\nabla \times \mathbf{B} = \nabla \times (\nabla \times \mathbf A) = \nabla (\nabla \cdot \mathbf{A}) - \nabla^2 \mathbf{A} = \mu_0 \mathbf J. \label{eq:curlcurlA}\tag{1}$$ Now we note that if we can find some solution $\mathbf{A}_1$ to this equation that satisfies the BCs¹, then we could find another one that has zero divergence $\mathbf{A}_2 = \mathbf{A}_1 + \nabla f$ if we just solve $\nabla^2 f = -\nabla \cdot \mathbf{A}_1$. That latter equation doesn’t even come with any BCs so it can definitely be satisfied. Thus, we conclude that if there’s a solution at all to our problem, then there’s a solution that has zero divergence. Thus we can impose $\nabla \cdot \mathbf A = 0$ from the start, which is called Coulomb gauge.

When we make that gauge choice, eq. \ref{eq:curlcurlA} becomes $\nabla^2 \mathbf{A} = - \mu_0 \mathbf J$. That’s just a Poisson equation for each component of $\mathbf A$, so we can solve it with a Green’s function. If we have all of 3D space as our domain, and fields die off quickly enough, this gives us $$\mathbf{A}(\mathbf{r}) = \frac{\mu_0}{4\pi} \int \frac{\mathbf{J}(\mathbf{r}’)}{|\mathbf{r}-\mathbf{r}’|} \, \mathrm{d}V’ . \label{eq:green}\tag{2}$$

All seems well so far, but we should be careful: it’s absolutely not trivial or obvious that the result eq. \ref{eq:green} has zero divergence! Just because we assumed $\nabla \cdot \mathbf{A} = 0$ to derive the Poisson equation from eq. \ref{eq:curlcurlA}, does not mean that the resulting Poisson equation will spit out a solution satisfying $\nabla \cdot \mathbf{A} = 0$. How can we be sure that in fact we are ok, and are not going to get inconsistent results?

Well, first off we can calculate the divergence of eq. \ref{eq:green} directly. Using $\nabla$ and $\nabla’$ to represent derivatives with respect to components of $\mathbf r$ and $\mathbf{r}’$ respectively, we find $$\begin{aligned} \nabla \cdot \mathbf{A} &= \frac{\mu_0}{4\pi} \int \nabla \cdot \left( \frac{\mathbf{J}(\mathbf{r}’)}{|\mathbf{r}-\mathbf{r}’|} \right) \, \mathrm{d}V’ \\ &= \frac{\mu_0}{4\pi} \int \mathbf{J}(\mathbf{r}’) \cdot \, \nabla \left( \frac{1}{|\mathbf{r}-\mathbf{r}’|} \right) \, \mathrm{d}V’ \\ &= - \frac{\mu_0}{4\pi} \int \mathbf{J}(\mathbf{r}’)\cdot \, \nabla’ \left( \frac{1}{|\mathbf{r}-\mathbf{r}’|} \right) \, \mathrm{d}V’ \\ &= - \frac{\mu_0}{4\pi} \int \nabla’ \cdot \left( \frac{\mathbf{J}(\mathbf{r}’) }{|\mathbf{r}-\mathbf{r}’|} \right) \, \mathrm{d}V’ + \frac{\mu_0}{4\pi} \int \frac{\cancelto{0}{ \nabla’ \cdot \mathbf{J}(\mathbf{r}’) } }{|\mathbf{r}-\mathbf{r}’|} \, \mathrm{d}V’ \\ &= - \frac{\mu_0}{4\pi} \oint \left( \frac{\mathbf{J}(\mathbf{r}’) }{|\mathbf{r}-\mathbf{r}’|} \right) \cdot \mathrm{d}\mathbf{S}’ , \end{aligned}$$ where crucially we’ve used the fact that we’re doing magnetostatics so $\nabla \cdot \mathbf{J} = 0$. Thus, as long as $\mathbf{J}$ decays fast enough as we send the boundary surface off to infinity (which I think it has to for eq. \ref{eq:green} to hold), we have $\nabla \cdot \mathbf A = 0$. Phew!

Here’s a slicker way to see that we’re ok: take the divergence of both sides of $\nabla^2 \mathbf{A} = - \mu_0 \mathbf J$ to find $\nabla^2 (\nabla \cdot \mathbf{A}) = 0$, where again we’ve used that fact that $\nabla \cdot \mathbf J = 0$ in magnetostatics. In the far field we have a BC that $\nabla \cdot \mathbf A = 0$ because we’ve assumed the fields die off at infinity, so one solution to this Laplace equation is just $\nabla \cdot \mathbf A = 0$ everywhere …but it’s a Laplace equation so the uniqueness theorem holds, so this is the only solution satisfying that BC. Thus, we again find that eq. \ref{eq:green} is guaranteed to consistently give $\nabla \cdot \mathbf{A} = 0$.

What if we want to make the same argument on a bounded domain? We need to enforce $\nabla \cdot \mathbf{A} = 0$ on the boundary, which will ensure that the solution to the Poisson equation has $\nabla \cdot \mathbf{A} = 0$ everywhere. Let’s think about a concrete example: a localized current distribution (e.g. loop of wire) inside a spherical cavity with superconducting walls. Superconductors have $\mathbf B = 0$ inside (Meissner effect), and $B_\perp$ is continuous across any interface, so we have a first physical boundary condition $\mathbf B \cdot \mathbf{n} = (\nabla \times \mathbf A) \cdot \mathbf n = \mathbf{0}$ at the walls. Let’s assume there isn’t current flowing into the walls, so $\mu_0 \mathbf J \cdot \mathbf n = \nabla \times \mathbf B \cdot \mathbf n = \nabla \times (\nabla \times \mathbf {A}) \cdot \mathbf n = 0$ at the walls. Using $\nabla^2 \mathbf{A} = - \mu_0 \mathbf J$ this gives us a second physical BC: $\nabla (\nabla \cdot \mathbf A ) \cdot \mathbf n = 0$. But via the same slick argument as before, after imposing $\nabla \cdot \mathbf A = 0$ on the boundary, the solution we find will have zero divergence everywhere, so this second physical BC will be satisfied automatically. Thus we actually still have some freedom in choosing the final BC, as long as it’s consistent with everything else.

A natural choice reveals itself by rewriting our first physical BC. Applying Stokes’s theorem, that BC is the same as requiring $\oint \mathbf{A} \cdot \mathrm{d}\mathbf{l} = 0$ for any loop on the boundary. This means that $\mathbf{A}_{||} = \nabla_S \zeta$ for some $\zeta$, where $\mathbf{A}_{||}$ is the projection of $\mathbf{A}$ into the tangent plane at the boundary, and $\nabla_S$ is the 2D gradient operator on the curved boundary surface (note this is just like an electrostatic field in a curved 2D universe). One neat way to satisfy this while adding one more BC is to take $\zeta$ to be constant, so $\mathbf{A}_{||} = \mathbf{0}$, i.e. require $\mathbf{A}$ to be perpendicular to the boundary. Thus we have a natural choice of three scalar BCs: $\nabla \cdot \mathbf A = 0$ and $\mathbf{A}_{||} = 0$. These will ensure that our physical BCs are satisfied, while also making the solutions to our three scalar Poisson equations fully determined.

The annoying thing is that although those three equations are uncoupled, the divergence BC couples together all three components of $\mathbf A$, so we actually do not straightforwardly have a standard Poisson equation problem for each component! If we’d chosen BCs $\mathbf A = \mathbf 0$ instead, everything would be uncoupled, but that would be wrong because it means we lose our consistency guarantee that the solution will have zero divergence. People have thought about this kind of thing, and actually found complex ways to decouple the problem by introducing some new fields (first paper below); but the point is that actually using the Coulomb gauge in a finite domain is tricky, and one has to be very careful! I think similar things could be said about other gauges like Lorenz.

For the non-believers, let’s see explicitly that you can get the wrong physical answer if you don’t enforce $\nabla \cdot \mathbf{A} = 0$ in a simple example: Suppose we want to find the $\mathbf{B}$ field produced by a current density that happens to be \begin{equation} \mathbf{J} = -\nabla^2 \left( \begin{pmatrix} x^2 \\ y^2 \\ z^2 \end{pmatrix} (1 - r)^3 / \mu_0 \right) \end{equation} inside a spherical cavity of radius 1 with superconducting ($B_\perp = 0$) walls, where $r \equiv \sqrt{x^2 +y^2 +z^2}$. Imagine someone’s evaluated the above expression and given $\mathbf J$ to us as a disgusting mess (which it is)! We want to solve $\nabla \times \nabla \times \mathbf{A} = \mu_0 \mathbf{J}$. We proclaim that we’ll use Coulomb gauge, so that we only have to solve $\nabla^2 \mathbf{A} = -\mu_0 \mathbf{J}$. Yay. In our carelessness, we then decide to use BCs $\mathbf{A} = \mathbf{0}$ on the boundary, feeling good because, via our previous loop-on-the-boundary argument, this BC ensures that $(\nabla \times \mathbf{A}) \cdot \mathbf{n} = 0$ on the boundary. Our vector Poisson problem has the unique solution \begin{equation} \mathbf{A} = \begin{pmatrix} x^2 \\ y^2 \\ z^2 \end{pmatrix} (1 - r)^3 . \end{equation} This $\mathbf{A}$ has $\nabla(\nabla \cdot \mathbf{A}) \neq \mathbf{0}$, as you can check. Thus this $\mathbf{A}$ does not satisfy $\nabla \times \nabla \times \mathbf{A} = \mu_0 \mathbf{J}$, so its curl is not the correct $\mathbf{B}$ field for this physical problem!

For more on this kind of issue, see

• Zhu et al. - ‘Finite element solution of vector Poisson equation with a coupling boundary’, which develops an approach to uncoupling the annoying finite domain problem by splitting into fields with different properties.

• Jiang et al - ‘The Origin of Spurious Solutions in Computational Electromagnetics’, which is a deep dive, and makes clear that lots of people have made real mistakes in this area when doing simulations.

• Kangro et al. - ‘Divergence boundary conditions for vector Helmholtz equations with divergence constraints’, which finds that actually the slick thing of enforcing the gauge condition via a BC rather than a bulk equation can fail for domains with sharp concave corners, so that you have to be even more careful apparently!

• Mayergoyz et al. - ‘A New Point of View on Mathematical Structure of Maxwell’s Equations’, which discusses essentially the same feature of the vector Helmholtz equation for $\mathbf E$ used in scattering theory.

• Dong et al. ‘Meshfree Finite Differences for Vector Poisson and Pressure Poisson Equations with Electric Boundary Conditions’ has some discussion of issues in this vein.

• Zhu et al. ‘3D vector Poisson-like problem with a triplet of intrinsic scalar boundary conditions’, which looks at potential inconsistency of BCs.

• Greengard et al - ‘A Consistency Condition for the Vector Potential in Multiply-Connected Domains’ is in a similar vein.

• Snider - ‘Note on the Coulomb gauge condition in magnetostatics’, which generally hammers home the need to actually make sure your solution satisfies your gauge condition.

Using the Coulomb gauge in a finite domain is a bit thorny --- if you start assuming div(A)=0 you'd better actually keep imposing it, otherwise you'll get the wrong answer: https://t.co/hhUeS2ykDr pic.twitter.com/Z7ZYqtApjH

— Daniel Duffy (@DDuffeh) September 17, 2022