Bessel function

Derivation of the 2D Green’s function for the wave equation operator.

October 22, 2025 math and physics play , , , , , , ,

[Click here for a PDF version of this post]

While it was difficult to attempt to verify the 2D Green’s function, it actually turns out to be fairly easy to derive it, provided we pick an alternate pole displacement from the 1D evaluation to make our lives easier.

With \( \Br = \Bx – \Bx’ \), and \( \tau = t – t’ \), and \( \epsilon > 0 \), we can form
\begin{equation}\label{eqn:waveEquationGreens:1340}
G_\epsilon(\Br, \tau) = \frac{c^2}{\lr{2 \pi}^3} \int d^2 \Bk d\omega \frac{ e^{j \Bk \cdot \Br + j \omega \tau}}{\lr{\omega -j \epsilon}^2 – \Bk^2 c^2 }
\end{equation}
This pole displacement has the nice property that both poles live in the upper half plane, so for \( \tau > 0 \), we have
\begin{equation}\label{eqn:waveEquationGreens:1360}
\begin{aligned}
G_\epsilon(\Br, \tau)
&= \frac{c^2}{\lr{2 \pi}^3} \int d^2 \Bk d\omega \frac{ e^{j \Bk \cdot \Br + j \omega \tau}}{
\lr{\omega -\lr{ \Abs{\Bk} c – j \epsilon}}
\lr{\omega -\lr{ -\Abs{\Bk} c – j \epsilon}}
} \\
&=
\Theta(\tau) \frac{c^2 j}{\lr{2 \pi}^2} \int d^2 \Bk e^{j \Bk \cdot \Br}
\lr{
\evalbar{
\frac{e^{ j \omega \tau}}{ \lr{\omega -\lr{ -\Abs{\Bk} c – j \epsilon}} }
}
{\omega = \Abs{\Bk} c – j \epsilon}
+
\evalbar{
\frac{e^{ j \omega \tau}}{ \lr{\omega -\lr{ \Abs{\Bk} c – j \epsilon}} }
}
{\omega = -\Abs{\Bk} c – j \epsilon}
}
\\
&=
\Theta(\tau) \frac{c^2 j}{\lr{2 \pi}^2} \int d^2 \Bk e^{j \Bk \cdot \Br} e^{ -j \epsilon \tau }
\lr{
\frac{e^{ j \Abs{\Bk} c \tau}}{ 2 \Abs{\Bk} c }
+
\frac{e^{ -j \Abs{\Bk} c \tau}}{ -2 \Abs{\Bk} c }
} \\
&=
\Theta(\tau) \frac{j^2 c}{\lr{2 \pi}^2} \int_{k=0}^\infty k dk \int_{\phi=0}^{2 \pi} d\phi e^{j k\Abs{\Br} \cos\phi } e^{ -j \epsilon \tau } \frac{\sin\lr{ k c \tau }}{k}.
\end{aligned}
\end{equation}
We’ve now successfully removed the singularity, and can evaluate the \(\epsilon \rightarrow 0 \) limit. We may also evaluate the \( \phi \) integral, remembering that
\begin{equation}\label{eqn:waveEquationGreens:1380}
\int_0^{2 \pi} e^{j \Abs{a} \cos\phi} d\phi = 2 \pi J_0(\Abs{a}),
\end{equation}
to find
\begin{equation}\label{eqn:waveEquationGreens:1400}
G(\Br, \tau) = -\Theta(\tau) \frac{c}{2 \pi} \int_{k=0}^\infty dk J_0(k\Abs{\Br}) \sin\lr{ k c \tau }.
\end{equation}
This integral yields easily to Mathematica, and we find
\begin{equation}\label{eqn:waveEquationGreens:1420}
G(\Br, \tau) = -\Theta(\tau) \frac{c}{2 \pi} \frac{\Theta(c \tau – \Abs{\Br})}{\sqrt{(c\tau)^2 – \Br^2}}.
\end{equation}
However, since \( \Theta(c \tau – \Abs{\Br}) = 1 \) only for \( \tau > \Abs{\Br}/c \), the \( \Theta(\tau) \) factor is redundant, and we find
\begin{equation}\label{eqn:waveEquationGreens:1440}
\boxed{
G(\Br, \tau) = – \frac{1}{2 \pi} \frac{\Theta(c \tau – \Abs{\Br})}{\sqrt{\tau^2 – \Br^2/c^2}},
}
\end{equation}
which matches the retarded Green’s function claimed by Grok.

Repeating this analysis for \( \tau < 0, \epsilon < 0 \), we find
\begin{equation}\label{eqn:waveEquationGreens:1460}
G(\Br, \tau) = -\Theta(-\tau) \frac{c}{2 \pi} \frac{\Theta(-c \tau – \Abs{\Br})}{\sqrt{(c\tau)^2 – \Br^2}},
\end{equation}
which we also see matches the Grok result for the advanced Green’s function. Both of these computations can be trivially performed in Mathematica following the same steps (taking all the fun from the story.) The advanced integral evaluation is shown in fig. 1 as an example.

fig. 1. Advanced 2D Green’s function for wave equation operator.

A trilogy in four+ parts: The 2D Laplacian Green’s function.

September 22, 2025 math and physics play , , , , , , ,

[Click here for a PDF version of this post]

I was questioning the correctness of the 1D and 2D Helmholtz Green’s functions derived above, since they are both seemingly malformed for \( k \rightarrow 0 \).

Let’s try to carefully expand the 2D Green’s function in the neighbourhood of \( k = 0 \) to validate that result, and as a side effect, obtain the Green’s function for the 2D Laplacian.

In [1], section 9.1.7, 9.1.8, we have for \( z \rightarrow 0 \)
\begin{equation}\label{eqn:helmholtzGreens:980}
\begin{aligned}
J_\nu(z) &\sim \lr{\frac{z}{2}}^\nu/\Gamma(\nu+1) \\
Y_0(z) &\sim \frac{2}{\pi} \ln z,
\end{aligned}
\end{equation}
so for \( k \ll r \)
\begin{equation}\label{eqn:helmholtzGreens:1000}
H_0^{(1)}(k r) \sim 1 + \frac{2j}{\pi} \ln\lr{k r},
\end{equation}
and
\begin{equation}\label{eqn:helmholtzGreens:1020}
\begin{aligned}
G(\Br)
&\sim -\frac{j}{4} \lr{ 1 + \frac{j}{2 \pi} \ln\lr{k r} } \\
&= -\frac{j}{4} + \frac{1}{2 \pi} \ln k + \frac{1}{2 \pi} \ln r.
\end{aligned}
\end{equation}
Here is where we have to get sneaky. Since we seek a Green’s function for the Laplacian operator, we are free to add any solution \( f(x,y) \) that satisfies \( \spacegrad^2 f = 0 \). Constants are clearly in that homogeneous solution space, so we may adjust this expansion of the Green’s function, throwing away the leading constant imaginary term, and treating \( k \) as a small constant, the \( \ln k \) term. That leaves us with
\begin{equation}\label{eqn:helmholtzGreens:1040}\boxed{
G(\Bx, \Bx’) = \frac{1}{2 \pi} \ln \Abs{\Bx – \Bx’}.
}
\end{equation}

Verifying the Laplacian Green’s function.

Let’s try to verify that this Green’s function is correct, since we’ve had lots of opportunities to screw up signs. We want to evaluate the Laplacian of the convolution and if all goes well, it should be \( V(\Bx) \). That is
\begin{equation}\label{eqn:helmholtzGreens:1060}
\begin{aligned}
\spacegrad^2 \int G(\Bx, \Bx’) V(\Bx’) d^2 \Bx’
&= \inv{2 \pi} \int \spacegrad^2 \ln \Abs{\Bx – \Bx’} V(\Bx’) d^2 \Bx’ \\
&= \inv{2 \pi}
\int V(\Bx’) \lr{\spacegrad’}^2 \ln \Abs{\Bx – \Bx’} d^2 \Bx’ \\
\end{aligned}
\end{equation}
We can verify that \( \lr{\spacegrad’}^2 \ln \Abs{\Bx – \Bx’} \) is zero whenever \( \Abs{\Bx – \Bx’} \ne 0 \). A nice way of doing that is in polar coordinates. Write
\begin{equation}\label{eqn:helmholtzGreens:1160}
\begin{aligned}
\Br &= \Bx’ – \Bx \\
r &= \Abs{r} \\
\end{aligned}
\end{equation}
and recall that
\begin{equation}\label{eqn:helmholtzGreens:1180}
\spacegrad^2 f = \inv{r} \PD{r}{} \lr{ r \PD{r}{f} } + \frac{\partial^2 f}{\partial \theta^2},
\end{equation}
but \( r \PD{r}{\ln r} = 1 \), and \( \ln r \) has no angular dependence. That means that
\begin{equation}\label{eqn:helmholtzGreens:1200}
\spacegrad^2 \int G(\Bx, \Bx’) V(\Bx’) d^2 \Bx’
= \inv{2 \pi} \int_{\Abs{\Bx – \Bx’} < \epsilon} V(\Bx’) \lr{\spacegrad’}^2 \ln \Abs{\Bx – \Bx’} d^2 \Bx’,
\end{equation}
where we let \( \epsilon \rightarrow 0 \). Such a region is illustrated in fig. 7.

fig. 7. Neighborhood around x

 

We can now apply Green’s theorem, which for 2D is
\begin{equation}\label{eqn:helmholtzGreens:1080}
\int_A \lr{ u \spacegrad^2 v – v \spacegrad^2 u } dA = \int_{\partial A} \lr{ u \spacegrad v – v \spacegrad u } \cdot \mathbf{\hat{n}} \, dS,
\end{equation}
or
\begin{equation}\label{eqn:helmholtzGreens:1090}
\int_A v \spacegrad^2 u \, dA = \int_A u \spacegrad^2 v \, dA -\int_{\partial A} \lr{ u \spacegrad v – v \spacegrad u } \cdot \mathbf{\hat{n}} \, dS,
\end{equation}
With
\begin{equation}\label{eqn:helmholtzGreens:1161}
\begin{aligned}
\mathbf{\hat{r}} &= \Br/r = \mathbf{\hat{n}} \\
u &= \ln r \\
v &= V(\Bx’) \\
dA &= r dr d\theta \\
dS &= r d\theta,
\end{aligned}
\end{equation}
we have
\begin{equation}\label{eqn:helmholtzGreens:1100}
\begin{aligned}
\int V(\Bx’) \lr{ \lr{\spacegrad’}^2 \ln \Abs{\Bx – \Bx’} } d^2 \Bx’
&=
\int_{r=0}^\epsilon r dr d\theta \ln r \lr{ \spacegrad’}^2 V(\Bx’) \\
&\quad-
\int_{\theta = 0}^{2 \pi}
\epsilon d\theta
\evalbar{
\lr{
\ln \epsilon \spacegrad’ V(\Bx’) – V(\Bx’) \spacegrad’ \ln r
}
\cdot \mathbf{\hat{r}}
}
{\,r = \epsilon}
\end{aligned}
\end{equation}
We have \( r \ln r \), or \( \epsilon \ln \epsilon \) dependence in two of the integrand terms, and with \( r < \epsilon \), and \( \epsilon \ln \epsilon \rightarrow 0 \), in the limit, we are left with
\begin{equation}\label{eqn:helmholtzGreens:1120}
\begin{aligned}
\int V(\Bx’) \lr{\spacegrad’}^2 \ln \Abs{\Bx – \Bx’} d^2 \Bx’
&=
\int_{\theta = 0}^{2 \pi} \epsilon \evalbar{ d\theta V(\Bx’) \lr{ \spacegrad’ \ln r } \cdot \mathbf{\hat{r}} }{\,r = \epsilon} \\
&=
\int_{\theta = 0}^{2 \pi}
\epsilon
\evalbar{
d\theta V(\Bx’) \lr{ \lr{ \mathbf{\hat{r}} \partial_r + \frac{\thetacap}{r^2} \partial_\theta } \ln r } \cdot \mathbf{\hat{r}}
}{\,r = \epsilon}
\\
&=
\int_{\theta = 0}^{2 \pi}
\epsilon
\evalbar{
d\theta V(\Bx’) \frac{\mathbf{\hat{r}}}{r} \cdot \mathbf{\hat{r}}
}{\,r = \epsilon}
\\
&=
\int_{\theta = 0}^{2 \pi} d\theta V(\Bx + \epsilon \mathbf{\hat{r}}).
\end{aligned}
\end{equation}

In the limit we are left with just \( 2 \pi V(\Bx) \), so
\begin{equation}\label{eqn:helmholtzGreens:1140}
\spacegrad^2 \int G(\Bx, \Bx’) V(\Bx’) d^2 \Bx’ = V(\Bx),
\end{equation}
as desired.

References

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

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.