# 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^{1}, 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.

If you liked this post, please subscribe to this blog to be notified via email about future posts, and/or follow me on twitter.

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

Mathematicans spend ages on this stuff and it’s mostly beyond me, but regarding existence of solutions I think a very helpful thing must be to realize that in physics we often have a variational principle, e.g. minimizing energy, extremising action…etc. I think that means there’s usually a nice ‘weak’ form of the problem, and then powerful things like the Lax-Milgram theorem can guarantee existence and/or uniqueness. ↩︎