Month: September 2025

Part 3/3. 2D Green’s functions for the Helmholtz (wave equation) operator.

September 20, 2025 math and physics play , , , , , , , , ,

[Click here for a PDF version of this post]

Having found the 1D and 3D Green’s function for the wave equation (Helmholtz) operator, we are now ready to attempt the harder 2D case again.

2D Green’s function.

Our starting place is
\begin{equation}\label{eqn:helmholtzGreens:680}
G(\Br) = -\inv{(2 \pi)^2} \int \frac{e^{j \Bp \cdot \Br}}{\Bp^2 – k^2} d^2 p.
\end{equation}
With a change of variables to polar coordinates, letting
\begin{equation}\label{eqn:helmholtzGreens:700}
\begin{aligned}
\Bp &= p \lr{ \cos\phi, \sin\phi } \\
\Br &= \Abs{\Br} \Be_2,
\end{aligned}
\end{equation}
we can make the integral explicit
\begin{equation}\label{eqn:helmholtzGreens:720}
G(\Br) = -\inv{(2 \pi)^2} \int_0^\infty \frac{p dp}{p^2 – k^2} \int_0^{2 \pi} d\phi e^{j p \Abs{\Br} \sin\phi}.
\end{equation}
Unlike the 3D case, where the angular dependence could be trivially evaluated, we are no longer so lucky. What on earth can we do with the \( \phi \) integral? Just like Hilter’s lament about “undoable integrals in Jackson”, we are faced with the same enemy. As it turns out, due to the cylindrical symmetry of the problem, we are also staring down the gun of Bessel functions. Both Mathematica and Grok point out that we can evaluate integrals of this form, like so:
\begin{equation}\label{eqn:helmholtzGreens:740}
\int_0^{2 \pi} d\phi e^{j a \sin\phi} = 2 \pi J_0(a).
\end{equation}

From [2] we have two representations of \( J_n \), a series representation and integral representation
\begin{equation}\label{eqn:helmholtzGreens:760}
J_n(z) = \sum_{m=0}^\infty \frac{(-1)^m (z/2)^{2m + n}}{(n + m)!m!} = \inv{\pi} \int_0^\pi \cos(n \theta – z \sin\theta) d\theta.
\end{equation}

In particular, this means that
\begin{equation}\label{eqn:helmholtzGreens:800}
J_0(z) = \sum_{m=0}^\infty \frac{(-1)^m (z/2)^{2m}}{(m!)^2} = J_0(z) = \inv{\pi} \int_0^\pi \cos(z \sin\theta) d\theta.
\end{equation}
This is a damped sine like function, as illustrated in fig. 4.

fig. 4. Bessel function of zeroth order.

 

In section 6.9, both of these are derived from a generating function representation of the Bessel functions, and one of the intermediate steps in that construction has
\begin{equation}\label{eqn:helmholtzGreens:840}
J_n(z) = \inv{2\pi} \int_{-\pi}^\pi e^{-j(n\theta – z \sin\theta)} d\theta,
\end{equation}
where the \( [-\pi, \pi] \) range was the result of a contour integration using a unit circle parameterization, which could have also used \( [0, 2 \pi] \). That means, sure enough, that we have
\begin{equation}\label{eqn:helmholtzGreens:860}
J_0(z) = \inv{2\pi} \int_{0}^{2\pi} e^{j z \sin\theta} d\theta,
\end{equation}
as claimed by both Grok and Mathematica.

This means that the evaluation of the Green’s function is now reduced to the limit of one final integral
\begin{equation}\label{eqn:helmholtzGreens:880}
G(\Br) = -\inv{2 \pi} \int_0^\infty \frac{p J_0(p \Abs{\Br} ) dp}{p^2 – \lr{k + j \epsilon}^2},
\end{equation}
where we’ve also displaced the problematic pole by a small imaginary amount as before. Grok incorrectly claimed that this was an even integral, and then argued that the end result is a Hankel function (that may be the case, but it’s reasoning to get there was clearly wrong.) Mathematica, on the other hand, can evaluate this integral
\begin{equation}\label{eqn:helmholtzGreens:900}
G(\Br) = -\inv{2 \pi} K_0\lr{\frac{\Abs{\Br}}{\sqrt{\frac{1}{(\epsilon – j k)^2}}}}, \epsilon \neq 0.
\end{equation}
It’s not clear to me why Mathematica writes the argument as 1 over a reciprocal root. Perhaps that has something to do with the branch cut that Mathematica uses for it’s square root function? If I plug in representitive numeric values, it simplifies in the expected way, as illustrated in fig. 5.

fig. 5. Mathematica weird Bessel argument.

The take away appears to be that the limiting form of the 2D Green’s function, for \( k > 0 \), is
\begin{equation}\label{eqn:helmholtzGreens:920}
G(\Bx, \Bx’) = -\inv{2 \pi} K_0\lr{-j k \Abs{\Bx – \Bx’} }.
\end{equation}
A peek at [1] shows that \( K_0 \) can be expressed in terms of a Hankel function of the first kind (order 0) \( H_0^{(1)}(z) = J_0(z) + j Y_0(z) \), plotted in fig. 6.

fig. 6. Hankel function of the first kind (order 0).

For real positive \( \alpha \), we have
\begin{equation}\label{eqn:helmholtzGreens:940}
K_0(-j \alpha) = \frac{j\pi}{2} H_0^{(1)}(\alpha),
\end{equation}
so
\begin{equation}\label{eqn:helmholtzGreens:960}
\boxed{
G(\Bx, \Bx’) = -\frac{j}{4} H_0^{(1)}(k \Abs{\Bx – \Bx’}).
}
\end{equation}

References

[1] M. Abramowitz and I.A. Stegun. Handbook of mathematical functions with formulas, graphs, and mathematical tables, volume 55. Dover publications, 1964.

[2] F.W. Byron and R.W. Fuller. Mathematics of Classical and Quantum Physics. Dover Publications, 1992.

Part 2. 3D Green’s function for the Helmholtz (wave equation) operator

September 20, 2025 math and physics play , , , , , , , ,

[Click here for a PDF version of this post]

This is a continuation of yesterday’s post on wave function Green’s functions. We derived the 1D Green’s function, now it’s time for the 3D.

Next up after this will be the 2D Green’s function, which looks harder to evaluate than the 1D or 3D versions.

3D Green’s function.

The 3D Green’s function that we wish to try to evaluate is
\begin{equation}\label{eqn:helmholtzGreens:540}
G(\Br) = -\inv{(2 \pi)^3} \int \frac{e^{j \Bp \cdot \Br}}{\Bp^2 – k^2} d^3 p.
\end{equation}
We will have to displace the pole again, but we will get to that in a bit. First let’s make a spherical change of variables to evaluate the integral, with
\begin{equation}\label{eqn:helmholtzGreens:560}
\begin{aligned}
\Bp &= p \lr{ \sin\alpha \cos\phi, \sin\alpha \sin\phi, \cos\alpha } \\
\Br &= \Abs{\Br} \Be_3.
\end{aligned}
\end{equation}
This gives
\begin{equation}\label{eqn:helmholtzGreens:580}
G(\Br)
= -\inv{(2 \pi)^3} \int_0^\infty p^2 dp \int_0^\pi \sin\alpha d\alpha \int_0^{2 \pi} d\phi \frac{e^{j p \Abs{\Br} \cos\alpha}}{p^2 – k^2}.
\end{equation}
Let \( t = \cos\alpha \), to find
\begin{equation}\label{eqn:helmholtzGreens:600}
\begin{aligned}
G(\Br)
&= -\inv{(2 \pi)^2} \int_0^\infty p^2 dp \int_1^{-1} (-dt) \frac{e^{j p \Abs{\Br} t}}{p^2 – k^2} \\
&= \inv{(2 \pi)^2} \int_0^\infty p^2 dp \evalrange{\frac{e^{j p \Abs{\Br} t}}{\lr{p^2 – k^2} j p \Abs{\Br}}}{1}{-1} \\
&= \inv{j (2 \pi)^2 \Abs{\Br}} \int_0^\infty p dp \frac{e^{-j p \Abs{\Br}} – e^{j p \Abs{\Br}} }{p^2 – k^2} \\
&= -\inv{j (2 \pi)^2 \Abs{\Br}} \int_{-\infty}^\infty p dp \frac{e^{j p \Abs{\Br}} }{p^2 – k^2} \\
&\sim -\inv{j (2 \pi)^2 \Abs{\Br}} \int_{-\infty}^\infty p dp \frac{e^{j p \Abs{\Br}} }{p^2 – \lr{k + j \epsilon}^2}.
\end{aligned}
\end{equation}
In the last step, we’ve displaced the pole so that we can evaluate it using an infinite upper half plane semicircular contour, as illustrated in fig 3.

fig 3. Contours for 3D Green’s function evaluation

Which pole we choose depends on the sign we pick for the “small” pole displacement \( \epsilon \). For the \( \epsilon > 0 \) case, we find
\begin{equation}\label{eqn:helmholtzGreens:620}
\begin{aligned}
G(\Br)
&= -\inv{j (2 \pi)^2 \Abs{\Br}} \int_{-\infty}^\infty p dp \frac{e^{j p \Abs{\Br}} }{\lr{p – (k + j \epsilon)}\lr{p – (-k – j \epsilon)}} \\
&= -\frac{2 \pi j}{j (2 \pi)^2 \Abs{\Br}} \evalbar{\frac{p e^{j p \Abs{\Br}} }{p + k + j \epsilon}}{p = k + j \epsilon} \\
&= -\frac{1}{2 \pi \Abs{\Br}} (k + j \epsilon) \frac{e^{j (k + j \epsilon) \Abs{\Br}} }{2 (k + j \epsilon)} \\
&= -\frac{1}{4 \pi \Abs{\Br}} e^{j k \Abs{\Br}} e^{-\epsilon \Abs{\Br}} \\
&\rightarrow -\frac{e^{j k \Abs{\Br}} }{4 \pi \Abs{\Br}}.
\end{aligned}
\end{equation}
whereas for \( \epsilon < 0 \), we have
\begin{equation}\label{eqn:helmholtzGreens:640}
\begin{aligned}
G(\Br)
&= -\inv{j (2 \pi)^2 \Abs{\Br}} \int_{-\infty}^\infty p dp \frac{e^{j p \Abs{\Br}} }{\lr{p – (k + j \epsilon)}\lr{p – (-k – j \epsilon)}} \\
&= -\frac{2 \pi j}{j (2 \pi)^2 \Abs{\Br}} \evalbar{\frac{p e^{j p \Abs{\Br}} }{p – k – j \epsilon}}{p = -k – j \epsilon} \\
&= -\frac{1}{2 \pi \Abs{\Br}} (-k – j \epsilon) \frac{e^{j (-k – j \epsilon) \Abs{\Br}} }{2 (-k – j \epsilon)} \\
&= -\frac{1}{4 \pi \Abs{\Br}} e^{-j k \Abs{\Br}} e^{\epsilon \Abs{\Br}} \\
&\rightarrow -\frac{e^{-j k \Abs{\Br}} }{4 \pi \Abs{\Br}}.
\end{aligned}
\end{equation}

The Green’s function has the structure of either an outgoing or incoming spherical wave, with inverse radial amplitude:
\begin{equation}\label{eqn:helmholtzGreens:660}
\boxed{
G(\Bx, \Bx’) = -\frac{e^{\pm j k \Abs{\Br}} }{4 \pi \Abs{\Br}}.
}
\end{equation}

Part 1. Green’s functions for the Helmholtz (wave equation) operator in various dimensions.

September 18, 2025 math and physics play , , , , , , , , ,

[Click here for a PDF version of this post]

My favorite book on mathematical physics derives the Green’s function for the 3D Helmholtz (wave equation) operator. I tried to derive the 2D Green’s function the same way and had trouble. In this series of blog posts, I’ll attempt that again, but will start with the easier 1D and 3D cases. Presuming that I don’t hit any conceptual troubles trying both of those from first principles, I’ll attempt the seemingly trickier 2D case again.

Motivation and background.

We seek a solution to non-homogeneous Helmholtz equation
\begin{equation}\label{eqn:helmholtzGreens:20}
0 = \lr{ \spacegrad^2 + k^2 } U(\Bx) – V(\Bx).
\end{equation}

This is a problem that can be solved using Fourier transform techniques. Following [1], let’s write our transform pair in the symmetrical form:
\begin{equation}\label{eqn:helmholtzGreens:40}
\begin{aligned}
F(\Bx) &= \lr{\inv{\sqrt{2 \pi}}}^N \int \hat{F}(\Bp) e^{j \Bp \cdot \Bx} d\Bp \\
\hat{F}(\Bp) &= \lr{\inv{\sqrt{2 \pi}}}^N \int F(\Bx) e^{-j \Bp \cdot \Bx} d\Bx.
\end{aligned}
\end{equation}

Expressing \(U(\Bx), V(\Bx)\), in terms of their Fourier transforms, \ref{eqn:helmholtzGreens:20} becomes
\begin{equation}\label{eqn:helmholtzGreens:80}
\begin{aligned}
0 &=
\lr{ \spacegrad^2 + k^2 }
\lr{\inv{\sqrt{2 \pi}}}^N
\int \hat{U}(\Bp) e^{j \Bp \cdot \Bx} d\Bp
–
\lr{\inv{\sqrt{2 \pi}}}^N
\int \hat{V}(\Bp) e^{j \Bp \cdot \Bx} d\Bp \\
&=
\lr{\inv{\sqrt{2 \pi}}}^N
\int \lr{ \hat{U}(\Bp) \lr{ -\Bp^2 + k^2 } – \hat{V}(\Bp) } e^{j \Bp \cdot \Bx} d\Bp,
\end{aligned}
\end{equation}
which requires
\begin{equation}\label{eqn:helmholtzGreens:100}
\hat{U}(\Bp) = \frac{\hat{V}(\Bp) }{k^2 – \Bp^2}.
\end{equation}

We can now inverse transform to find \( U(\Bx) \), which gives
\begin{equation}\label{eqn:helmholtzGreens:120}
\begin{aligned}
U(\Bx) &=
\lr{\inv{\sqrt{2 \pi}}}^N \int \hat{U}(\Bp) e^{j \Bp \cdot \Bx} d\Bp \\
&=
\lr{\inv{\sqrt{2 \pi}}}^N \int
\frac{\hat{V}(\Bp) }{k^2 – \Bp^2} e^{j \Bp \cdot \Bx} d\Bp \\
&=
\lr{\inv{2 \pi}}^N \int
\inv{k^2 – \Bp^2} e^{j \Bp \cdot \Bx} d\Bp \int V(\Bx’) e^{-j \Bp \cdot \Bx’} d\Bx’ \\
&=
\int V(\Bx’) d\Bx’
\lr{\inv{2 \pi}}^N
\int
\inv{k^2 – \Bp^2} e^{j \Bp \cdot \lr{\Bx – \Bx’}}
d\Bp
.
\end{aligned}
\end{equation}

The general solution is given by
\begin{equation}\label{eqn:helmholtzGreens:140}
U(\Bx) = \int G(\Bx, \Bx’) V(\Bx’) d\Bx’,
\end{equation}
where \( G(\Bx, \Bx’) \) is called the Green’s function, and has the form
\begin{equation}\label{eqn:helmholtzGreens:160}
G(\Bx, \Bx’) = \lr{\inv{2 \pi}}^N
\int
\inv{k^2 – \Bp^2} e^{j \Bp \cdot \lr{\Bx – \Bx’}}
d\Bp.
\end{equation}

Equivalently, if we presume that a solution of the form \ref{eqn:helmholtzGreens:140} can be found, and operate on that with the Helmholtz operator \( \spacegrad^2 + k^2 \), we find
\begin{equation}\label{eqn:helmholtzGreens:60}
\lr{ \spacegrad^2 + k^2 } U(\Bx) = \int \lr{ \spacegrad^2 + k^2 } G(\Bx, \Bx’) V(\Bx’) d\Bx’ = V(\Bx),
\end{equation}
which requires that our Green’s function \( G(\Bx, \Bx’) \) has the functional form
\begin{equation}\label{eqn:helmholtzGreens:180}
\lr{ \spacegrad^2 + k^2 } G(\Bx, \Bx’) = \delta(\Bx – \Bx’).
\end{equation}

Evaluating the Green’s function in 1D.

For the one dimensional case, we want to evaluate
\begin{equation}\label{eqn:helmholtzGreens:200}
G(u) = -\inv{2 \pi} \int \inv{p^2 – k^2} e^{j p u} dp,
\end{equation}
an integral which is unfortunately non-convergent. Since we are dealing with delta functions, it is not surprising that we have convergence problems. The technique used in the book is to displace the pole slightly by a small imaginary amount, and then take the limit.

That is
\begin{equation}\label{eqn:helmholtzGreens:220}
G(u) = \lim_{\epsilon \rightarrow 0} G_\epsilon(u),
\end{equation}
where
\begin{equation}\label{eqn:helmholtzGreens:240}
\begin{aligned}
G_\epsilon(u)
&= -\inv{2 \pi} \int_{-\infty}^\infty \inv{p^2 – \lr{k + j \epsilon}^2} e^{j p u} dp \\
&= -\inv{2 \pi} \int_{-\infty}^\infty \frac{e^{j p u}}{\lr{ p – k -j \epsilon}\lr{p + k + j \epsilon}} dp.
\end{aligned}
\end{equation}
For \( u > 0 \) we can use an upper half plane infinite semicircular contour integral, as illustrated in fig. 1, where we
let \( R \rightarrow \infty \).

fig.1 Contour for u > 0.

The residue calculation for that contour gives
\begin{equation}\label{eqn:helmholtzGreens:260}
\begin{aligned}
G_\epsilon(u)
&= -\frac{2 \pi j}{2 \pi} \evalbar{\frac{e^{j p u}}{p + k + j \epsilon} }{p = k + j \epsilon} \\
&= -j \frac{e^{j \lr{k + j \epsilon} u}}{2\lr{k + j \epsilon}} \\
&= -j \frac{e^{j k u} e^{-\epsilon u}}{2\lr{k + j \epsilon}} \\
&\rightarrow -\frac{j}{2k} e^{j k u}.
\end{aligned}
\end{equation}

For \( u < 0 \) we can use a lower half plane infinite semicircular contour, as illustrated in fig. 2.

fig 2. Contour for u < 0.

For this contour, we find
\begin{equation}\label{eqn:helmholtzGreens:280}
\begin{aligned}
G_\epsilon(u)
&= -\frac{2 \pi j}{2 \pi} \evalbar{\frac{e^{j p u}}{p – k – j \epsilon} }{p = -k – j \epsilon} \\
&= j \frac{e^{-j \lr{k + j \epsilon} u}}{2\lr{k + j \epsilon}} \\
&= j \frac{e^{-j k u} e^{\epsilon u}}{2\lr{k + j \epsilon}} \\
&\rightarrow \frac{j}{2k} e^{-j k u} \\
&= \frac{j}{2k} e^{j k \Abs{u}}.
\end{aligned}
\end{equation}
We find that our Green’s function is
\begin{equation}\label{eqn:helmholtzGreens:300}
\boxed{
G(u) = \frac{-j \mathrm{sgn}(u)}{2k} e^{j k \Abs{u}}.
}
\end{equation}

Let’s plug this into the convolution integral to see the form of the general solution
\begin{equation}\label{eqn:helmholtzGreens:320}
U(x) = -\frac{j}{2k} \int_{-\infty}^\infty \mathrm{sgn}(x – x’) e^{j k \Abs{x – x’}} V(x’) dx’.
\end{equation}
We want to break this integral into two regions
\begin{equation}\label{eqn:helmholtzGreens:340}
\int_{-\infty}^\infty = \int_{-\infty}^x + \int_x^\infty,
\end{equation}
separating the integral into regions where \( x > x’ \) and \( x < x’ \) respectively. That is
\begin{equation}\label{eqn:helmholtzGreens:360}
U(x) =
-\frac{j}{2k} \int_{-\infty}^x e^{j k \lr{x – x’}} V(x’) dx’
+\frac{j}{2k} \int_x^\infty e^{-j k \lr{x – x’}} V(x’) dx’.
\end{equation}
This isn’t the most general solution, as we can also add any solution to the homogeneous Helmholtz equation. That is
\begin{equation}\label{eqn:helmholtzGreens:400}
U(x) = A e^{j k x} + B e^{-j k x} – \frac{j}{2k} \int_{-\infty}^x e^{j k \lr{x – x’}} V(x’) dx’
+\frac{j}{2k} \int_x^\infty e^{-j k \lr{x – x’}} V(x’) dx’.
\end{equation}

The real and imaginary parts of this equation must also be independent solutions. For example, taking the real parts, we find the following general solution
\begin{equation}\label{eqn:helmholtzGreens:380}
U(x) =
A’ \cos\lr{ k x } +
B’ \sin\lr{ k x }
+\inv{2} \int_{-\infty}^x \frac{\sin\lr{k \lr{x – x’}}}{k} V(x’) dx’
-\inv{2} \int_x^\infty \frac{\sin\lr{k \lr{x – x’}}}{k} V(x’) dx’.
\end{equation}

A strictly causal solution.

It is interesting that the specific solution above has equal causal and acausal contributions. Such a solution (outside of QFT) is generally undesirable. We can construct a specific solution that is either causal or acausal by picking just one of the integrals above, instead of averaging. For example, let
\begin{equation}\label{eqn:helmholtzGreens:440}
f(x) = \int_{-\infty}^x \frac{\sin\lr{k \lr{x – x’}}}{k} V(x’) dx’.
\end{equation}
We can verify that this is a specific solution to our equation using the identity
\begin{equation}\label{eqn:helmholtzGreens:460}
\frac{d}{dx} \int_a^x g(x, x’) dx’
=
\evalrange{g(x, x’) }{a}{x} + \int_a^x \frac{\partial g(x,x’)}{dx} dx’.
\end{equation}
Taking the first derivative of \( f(x) \), we find
\begin{equation}\label{eqn:helmholtzGreens:480}
\begin{aligned}
\frac{df}{dx}
&= \evalrange{ \frac{\sin\lr{k \lr{x – x’}}}{k} V(x’) }{-\infty}{x} + k \int_{-\infty}^x \frac{\cos\lr{k \lr{x – x’}}}{k} V(x’) dx’ \\
&= \frac{\sin\lr{k \lr{x + \infty}}}{k} V(-\infty) + k \int_{-\infty}^x \frac{\cos\lr{k \lr{x – x’}}}{k} V(x’) dx’ \\
&= k \int_{-\infty}^x \frac{\cos\lr{k \lr{x – x’}}}{k} V(x’) dx’,
\end{aligned}
\end{equation}
where we have assumed that our forcing function \( V(x) \) is zero at \( -\infty \). Taking the second derivative, we have
\begin{equation}\label{eqn:helmholtzGreens:500}
\begin{aligned}
\frac{d^2f}{dx^2}
&=
\evalrange{ k \frac{\cos\lr{k \lr{x – x’}}}{k} V(x’) }{-\infty}{x} – k^2 \int_{-\infty}^x \frac{\sin\lr{k \lr{x – x’}}}{k} V(x’) dx’ \\
&= V(x) – k^2 f(x),
\end{aligned}
\end{equation}
or
\begin{equation}\label{eqn:helmholtzGreens:520}
\frac{d^2}{dx^2} f(x) + k^2 f(x) = V(x).
\end{equation}
This verifies that \ref{eqn:helmholtzGreens:440} is also a specific solution to the wave equation, as expected and desired.

It appears that the general solution is likely of the following form
\begin{equation}\label{eqn:helmholtzGreens:380b}
U(x) =
A’ \cos\lr{ k x } +
B’ \sin\lr{ k x }
+\alpha \int_{-\infty}^x \frac{\sin\lr{k \lr{x – x’}}}{k} V(x’) dx’
-(1-\alpha)\int_x^\infty \frac{\sin\lr{k \lr{x – x’}}}{k} V(x’) dx’,
\end{equation}
with \( \alpha \in [0,1] \).

It’s pretty cool that we can completely solve the 1D forced wave equation, for any forcing function, from first principles. Yes, I took liberties that would make a mathematician cringe, but we are telling a story, and leaving the footnotes to somebody else.

More specific boundary constraints.

Just as we have the freedom to add any homogeneous solution to our specific convolution based solution, we may do so for the Green’s function itself. Our process above, implicitly assumes that we are interested in infinite boundary value constraints. Should we wish to impose different boundary constraints, we can form
\begin{equation}\label{eqn:helmholtzGreens:420}
G(u) = A e^{ j k u} + B e^{-j k u} – \frac{j \mathrm{sgn}(u)}{2k} e^{j k \Abs{u}},
\end{equation}
but must then use the boundary value constraints to determine the desired form of the Green’s function, using the two degrees of freedom to do so. That’s also an interesting topic, and would be good to also visit in a followup post.

References

[1] F.W. Byron and R.W. Fuller. Mathematics of Classical and Quantum Physics. Dover Publications, 1992.

An “easy” integral from Jackson’s electrodynamics

September 6, 2025 math and physics play , , ,

Screenshot

[Click here for a PDF version of this post]

Once again, I was reading my Jackson [1], which characteristically had the statement “the […] integral can easily be shown to have the value \( 4 \pi \)”, in a discussion of electrostatic energy and self energy.

The integral is
\begin{equation}\label{eqn:selfEnergyIntegral:20}
I = \int \frac{\Brho}{\rho^3} \cdot \frac{\Brho + \Bn}{\Norm{\Brho + \Bn}^3} d^3 \rho.
\end{equation}

This is something that I once figured out once before (see [2] appendix C). However, trying to do it a second time around, I think that I found the “easy” way.

As Jackson hints, the starting point is
\begin{equation}\label{eqn:selfEnergyIntegral:40}
\frac{\Bx}{\Norm{\Bx}^3}
=
-\spacegrad \inv{\Norm{\Bx}},
\end{equation}
but we don’t have to apply it to both the vector terms, as I did in my initial attempt (which results in a Laplacian to reduce.) Inserting this and applying chain rule, we find
\begin{equation}\label{eqn:selfEnergyIntegral:60}
\begin{aligned}
I
&= -\int \frac{\Brho}{\rho^3} \cdot \spacegrad_\Brho \inv{\Norm{\Brho + \Bn}} d^3 \rho \\
&=
-\int
\spacegrad_\Brho \cdot \lr{
\frac{\Brho}{\rho^3} \cdot
\inv{\Norm{\Brho + \Bn}}
}
d^3 \rho
+
\int
\lr{
\spacegrad_\Brho \cdot
\frac{\Brho}{\rho^3}
}
\inv{\Norm{\Brho + \Bn}}
d^3 \rho
\\
\end{aligned}
\end{equation}

The first integral can be evaluated using an infinite spherical shell
\begin{equation}\label{eqn:selfEnergyIntegral:80}
\begin{aligned}
I_1
&= -\int
\spacegrad_\Brho \cdot \lr{
\frac{\Brho}{\rho^3} \cdot
\inv{\Norm{\Brho + \Bn}}
}
d^3 \rho \\
&=
\lim_{\rho \rightarrow \infty}
-\frac{\Brho}{\rho} \cdot \lr{
\frac{\Brho}{\rho^3} \cdot
\inv{\Norm{\Brho + \Bn}}
} 4 \pi \rho^2 \\
&=
\lim_{\rho \rightarrow \infty}
\frac{-4 \pi}{\Norm{\Brho + \Bn}} \\
&=
0.
\end{aligned}
\end{equation}

The divergence term in the second integral, provided \( \Bx \ne 0 \), has the form
\begin{equation}\label{eqn:selfEnergyIntegral:100}
\begin{aligned}
\spacegrad \cdot \frac{\Bx}{\Norm{\Bx}^3}
&=
\inv{\Norm{\Bx}^3} \spacegrad \cdot \Bx
+
\lr{ \Bx \cdot \spacegrad } \inv{\Norm{\Bx}^3} \\
&=
\frac{3}{\Norm{\Bx}^3}
+
2 \frac{x_k x_j}{\Norm{\Bx}^5} \lr{-\frac{3}{2}} \partial_k x_j \\
&=
\frac{3}{\Norm{\Bx}^3}
– \frac{3}{\Norm{\Bx}^3}
\end{aligned}
\end{equation}
However, in a neighbourhood of the origin, this actually has a delta function structure. We can see that from Gauss’s law, where we have
\begin{equation}\label{eqn:selfEnergyIntegral:120}
\spacegrad \cdot \BE = \frac{\rho}{\epsilon_0}.
\end{equation}
If we plug in the integral representation of \( \BE \) on the LHS, we have
\begin{equation}\label{eqn:selfEnergyIntegral:140}
\begin{aligned}
\spacegrad \cdot \BE
&=
\spacegrad \cdot \int \frac{\rho(\Bx’)}{4 \pi \epsilon_0} \frac{\Bx – \Bx’}{\Norm{\Bx – \Bx’}^3} d^3 x \\
&=
\int \frac{\rho(\Bx’)}{4 \pi \epsilon_0} \spacegrad \cdot \frac{\Bx – \Bx’}{\Norm{\Bx – \Bx’}^3} d^3 x \\
&=
-\int \frac{\rho(\Bx’)}{4 \pi \epsilon_0} \spacegrad’ \cdot \frac{\Bx – \Bx’}{\Norm{\Bx – \Bx’}^3} d^3 x.
\end{aligned}
\end{equation}
Comparing the LHS and RHS, we must have
\begin{equation}\label{eqn:selfEnergyIntegral:160}
\spacegrad’ \cdot \frac{\Bx’ – \Bx}{\Norm{\Bx’ – \Bx}^3} = 4 \pi \delta^3\lr{\Bx’ – \Bx}.
\end{equation}

We can now substitute that into the second integral to find
\begin{equation}\label{eqn:selfEnergyIntegral:n}
\begin{aligned}
I_2 &=
\int
\lr{
\spacegrad_\Brho \cdot
\frac{\Brho}{\rho^3}
}
\inv{\Norm{\Brho + \Bn}}
d^3 \rho \\
&=
\frac{4 \pi}{\Norm{\Bn}} \\
&=
4 \pi.
\end{aligned}
\end{equation}

Sure enough, the integral has a \( 4 \pi \) value. But was that easy?  I think Hitler would disagree.

References

[1] JD Jackson. Classical Electrodynamics. John Wiley and Sons, 2nd edition, 1975.

[2] Peeter Joot. Electromagnetic Theory. Kindle Direct Publishing, Toronto, 2016.