Green's functions for partial differential equations

Green's functions for the Laplace, Helmholtz, wave and heat equations.

In this chapter 5

Tick a unit to add it to your PDF basket.

Green's functions in several dimensions

Green's functions in several dimensions

Theorem

Green's functions for Laplace and Poisson problems

Adapted

Add to PDF

Defines the Dirichlet and Neumann Green's functions of a bounded region for $\nabla^2$ , derives the representation formula that expresses the solution of $\nabla^2 u = \rho$ from the data alone, and proves that the Dirichlet Green's function is symmetric in its two arguments.

Prerequisites: Fundamental solution of Laplace's equation, Green's identities, The Dirac delta in calculations

Note(Sign and naming conventions)Original

Books disagree about the sign of the Green's function of the Laplacian, and the disagreement is pure bookkeeping: a change of sign in $G$ is compensated by a change of sign in the integral that uses it. Throughout this unit we take

\nabla^2 G(\mathbf{r}, \mathbf{r}_0) = \delta(\mathbf{r} - \mathbf{r}_0) ,

which is the convention of Riley, Hobson and Bence, of Arfken, Weber and Harris, and of Hassani. In three dimensions $G$ is then negative near its pole, behaving like $-1/(4\pi|\mathbf{r}-\mathbf{r}_0|)$ . Evans instead defines $G$ by $-\Delta G = \delta_x$ , so his Green's function is the negative of ours; Strauss writes down the same function as ours but calls it "the Green's function for the operator $-\Delta$ ", because he thinks of the Dirichlet problem in the form $-\Delta u =$ source. Arfken comments explicitly that previous editions of that book used the other sign. The remark at the end of the unit lists the resulting differences in the formulas.

Definition(Fundamental solution, Dirichlet and Neumann Green's functions)AdaptedAdapted from Riley, Hobson & Bence (2006), §21.5.2-21.5.4, Eqs. (21.82), (21.85) and p. 766; Strauss (2008), §7.3, the three properties, p. 188; Evans (2010), §2.2.4, the corrector problem (26) and the definition on p. 33

Let $V \subseteq \mathbb{R}^3$ be a bounded domain whose boundary $S = \partial V$ is of class $C^1$ , with unit outward normal $\hat{\mathbf{n}}$ and total area $A = \oint_S dS$ . Write

F(\mathbf{r}, \mathbf{r}_0) = -\frac{1}{4\pi\lvert\mathbf{r}-\mathbf{r}_0\rvert}

for the fundamental solution of $\nabla^2$ in three dimensions, which satisfies $\nabla^2 F = \delta(\mathbf{r}-\mathbf{r}_0)$ and decays at infinity.

Fix $\mathbf{r}_0 \in V$ . A function $G(\cdot\,, \mathbf{r}_0)$ defined on $\overline{V}\setminus\{\mathbf{r}_0\}$ is a Green's function for $\nabla^2$ on $V$ with pole at $\mathbf{r}_0$ if

$G(\cdot\,,\mathbf{r}_0)$ is $C^2$ and harmonic on $V \setminus \{\mathbf{r}_0\}$ ; and
the corrector $H(\mathbf{r},\mathbf{r}_0) = G(\mathbf{r},\mathbf{r}_0) - F(\mathbf{r},\mathbf{r}_0)$ extends to a function of class $C^2(\overline{V})$ that is harmonic on all of $V$ , including at $\mathbf{r}_0$ .

Condition 2 says that $G$ has exactly the singularity of $F$ at its pole and nothing else; with it, $\nabla^2 G = \delta(\mathbf{r}-\mathbf{r}_0)$ holds in the sense of distributions. A Green's function is pinned down by a boundary condition, and two are standard:

the Dirichlet Green's function $G_D$ , required to satisfy $G_D(\mathbf{r},\mathbf{r}_0) = 0$ for $\mathbf{r} \in S$ ;
the Neumann Green's function $G_N$ , required to satisfy $\frac{\partial G_N}{\partial n}(\mathbf{r},\mathbf{r}_0) = \frac{1}{A} \qquad\text{for } \mathbf{r} \in S .$

Strauss states the Dirichlet case as three conditions: $G$ harmonic in $D$ except at $x_0$ and twice continuously differentiable there, $G = 0$ on the boundary, and $G(x) + 1/(4\pi\lvert x - x_0\rvert)$ finite at $x_0$ , with continuous second derivatives everywhere and harmonic at $x_0$ . Our conditions 1 and 2 are that third condition restated through the corrector, which is how Evans and Riley, Hobson and Bence organise it: Evans builds $G$ by solving $\Delta\phi^x = 0$ in $U$ with $\phi^x = \Phi(y-x)$ on $\partial U$ and then subtracting, and Riley, Hobson and Bence write $G = F + H$ with $H$ a solution of Laplace's equation in $V$ adjusted so that the sum vanishes on $S$ .

Remark(Why the Neumann condition cannot be

\partial G/\partial n = 0

)AdaptedAdapted from Riley, Hobson & Bence (2006), §21.5.4, the consistency condition and the choice $1/A$, p. 766

The obvious analogue of the Dirichlet condition would be $\partial G/\partial n = 0$ on $S$ , and it is not available. Integrating $\nabla^2 G = \delta(\mathbf{r}-\mathbf{r}_0)$ over $V$ and converting the left side with the divergence theorem gives

\oint_S \frac{\partial G}{\partial n}\, dS = \int_V \nabla^2 G\; dV = 1 ,

so the normal derivative of any Green's function has total flux $1$ through $S$ and cannot vanish identically. The cheapest way to meet the constraint is to spread the flux uniformly, $\partial G/\partial n = 1/A$ ; this is the choice Riley, Hobson and Bence call the Neumann Green's function, and the one used below.

Theorem(Representation formula and the solution of Poisson's equation)AdaptedAdapted from Riley, Hobson & Bence (2006), §21.5.2 Eq. (21.84) p. 756, §21.5.3 Eq. (21.86) p. 756, §21.5.4 Eq. (21.106) p. 766; Strauss (2008), §7.3, Theorems 1 and 2, pp. 188-190; Evans (2010), §2.2.4, Theorem 12, p. 34; Hassani (2013), §22.1, Eq. (22.2), p. 666

Let $V$ be a bounded domain with $C^1$ boundary $S$ , let $\rho$ be continuous on $\overline{V}$ , and let $u \in C^2(\overline{V})$ solve

\nabla^2 u = \rho \quad\text{in } V .

Let $G(\cdot\,,\mathbf{r}_0)$ be any Green's function for $\nabla^2$ on $V$ with pole at $\mathbf{r}_0 \in V$ , in the sense of the definition above. Then

u(\mathbf{r}_0) \;=\; \int_V G(\mathbf{r},\mathbf{r}_0)\,\rho(\mathbf{r})\; dV(\mathbf{r}) \;+\; \oint_S \left[ u(\mathbf{r})\,\frac{\partial G(\mathbf{r},\mathbf{r}_0)}{\partial n} - G(\mathbf{r},\mathbf{r}_0)\,\frac{\partial u(\mathbf{r})}{\partial n} \right] dS(\mathbf{r}) , \tag{R}

all derivatives and integrations being with respect to the first argument $\mathbf{r}$ . Two specialisations remove the term that is not data:

Dirichlet problem. If $u = f$ on $S$ and $G = G_D$ vanishes on $S$ , the second boundary term drops out and $u(\mathbf{r}_0) = \int_V G_D(\mathbf{r},\mathbf{r}_0)\,\rho(\mathbf{r})\; dV + \oint_S f(\mathbf{r})\,\frac{\partial G_D(\mathbf{r},\mathbf{r}_0)}{\partial n}\; dS . \tag{D}$
Neumann problem. If $\partial u/\partial n = f$ on $S$ and $G = G_N$ has $\partial G_N/\partial n = 1/A$ on $S$ , the first boundary term becomes $A^{-1}\oint_S u\, dS = \langle u\rangle_S$ , the average of $u$ over $S$ , and $u(\mathbf{r}_0) = \int_V G_N(\mathbf{r},\mathbf{r}_0)\,\rho(\mathbf{r})\; dV + \langle u\rangle_S - \oint_S G_N(\mathbf{r},\mathbf{r}_0)\, f(\mathbf{r})\; dS . \tag{N}$ Here $\langle u \rangle_S$ is not data: it is the one constant that Neumann data leave undetermined (Green's identities), and may be assigned freely.

In (D) the solution has been expressed entirely through the data $\rho$ and $f$ ; the two terms are, in Hassani's reading, the response to the source inside and the response to the boundary values. Both reduce to $u(\mathbf{r}_0) = \int_V G\rho\, dV$ when the data on $S$ vanish, which is the elementary form of the Green's function method quoted by Riley, Hobson and Bence as Eq. (21.78).

Idea of the proofOriginal

Everything comes from Green's second identity, which pairs two functions and trades their Laplacians inside the region against their boundary values and normal derivatives. Put $u$ , the unknown, in one slot. In the other slot we are free to put whatever we like, and the whole art is in the choice.

Choosing the fundamental solution $F$ is the first step: $F$ is harmonic, so its Laplacian contributes nothing to the volume integral, except that its singularity at $\mathbf{r}_0$ is exactly strong enough to leave behind the single number $u(\mathbf{r}_0)$ . (Technically we cut a small ball around the pole out of $V$ , apply the identity on what is left, and watch the sphere of radius $\varepsilon$ contribute the spherical average of $u$ , which tends to $u(\mathbf{r}_0)$ .) That already yields a formula for $u(\mathbf{r}_0)$ , but it involves both $u$ and $\partial u/\partial n$ on $S$ , and in a boundary-value problem only one of them is known.

The second step repairs this. Adding to $F$ any function $H$ that is harmonic throughout $V$ changes nothing in the volume term, because Green's second identity applied to the pair $u, H$ has no singularity to produce; it only changes the boundary terms. So we may shop among correctors $H$ for the one that makes the unwanted boundary term vanish — $G = 0$ on $S$ kills the $\partial u/\partial n$ term, and $\partial G/\partial n$ constant converts the $u$ term into the harmless surface average. The Green's function is $F$ plus the corrector that the boundary condition selects.

Proof(Proof of (R), (D) and (N))AdaptedAdapted from Strauss (2008), §7.2, pp. 186-187 and §7.3, pp. 188-189; Riley, Hobson & Bence (2006), §21.5.2, pp. 754-756; Evans (2010), §2.2.4, pp. 33-34

Write $r = \lvert\mathbf{r}-\mathbf{r}_0\rvert$ and $F = F(\mathbf{r},\mathbf{r}_0) = -1/(4\pi r)$ , and split $G = F + H$ with $H$ as in condition 2 of the definition.

Step 1 (the corrector contributes only boundary terms). $H(\cdot\,,\mathbf{r}_0)$ and $u$ both lie in $C^2(\overline{V})$ and $\nabla^2 H = 0$ in $V$ , so Green's second identity applies to the pair on $V$ with no excision:

\int_V \bigl( u\,\nabla^2 H - H\,\nabla^2 u \bigr)\, dV = \oint_S \left( u\,\frac{\partial H}{\partial n} - H\,\frac{\partial u}{\partial n}\right) dS .

Since $\nabla^2 H = 0$ and $\nabla^2 u = \rho$ , this reads

- \int_V H\rho\; dV = \oint_S \left( u\,\frac{\partial H}{\partial n} - H\,\frac{\partial u}{\partial n}\right) dS . \tag{1}

Step 2 (the fundamental solution, with the pole cut out). Let $B_\varepsilon$ denote the open ball $\lvert\mathbf{r}-\mathbf{r}_0\rvert < \varepsilon$ , with $\varepsilon$ small enough that its closure lies inside $V$ , and excise it: $V_\varepsilon = V\setminus\overline{B_\varepsilon}$ . On $V_\varepsilon$ the function $F$ is smooth and harmonic, so Green's second identity gives

- \int_{V_\varepsilon} F\rho\; dV = \oint_S \left( u\,\frac{\partial F}{\partial n} - F\,\frac{\partial u}{\partial n}\right) dS \;+\; \oint_{\partial B_\varepsilon} \left( u\,\frac{\partial F}{\partial n} - F\,\frac{\partial u}{\partial n}\right) dS , \tag{2}

the boundary of $V_\varepsilon$ consisting of $S$ and the sphere $\partial B_\varepsilon$ . The hypothesis $u\in C^2(\overline V)$ is used here and nowhere else: it is what lets the identity be applied at all and what makes $u$ and $\nabla u$ bounded near $\mathbf{r}_0$ .

Step 3 (the small sphere delivers $-u(\mathbf{r}_0)$ ). On $\partial B_\varepsilon$ the outward normal of $V_\varepsilon$ points towards $\mathbf{r}_0$ , so $\partial/\partial n = -\partial/\partial r$ there. With $F = -1/(4\pi r)$ ,

F\big|_{r=\varepsilon} = -\frac{1}{4\pi\varepsilon}, \qquad \frac{\partial F}{\partial n}\bigg|_{r=\varepsilon} = -\frac{\partial}{\partial r}\left(-\frac{1}{4\pi r}\right)\bigg|_{r=\varepsilon} = -\frac{1}{4\pi\varepsilon^2} .

Both are constant on the sphere, whose area is $4\pi\varepsilon^2$ , so

\oint_{\partial B_\varepsilon}\left( u\,\frac{\partial F}{\partial n} - F\,\frac{\partial u}{\partial n}\right) dS = -\frac{1}{4\pi\varepsilon^2}\oint_{\partial B_\varepsilon} u\, dS + \frac{1}{4\pi\varepsilon}\oint_{\partial B_\varepsilon}\frac{\partial u}{\partial n}\, dS = -\overline{u}_\varepsilon + \varepsilon\,\overline{\left(\frac{\partial u}{\partial n}\right)}_{\!\varepsilon} ,

where a bar denotes the average over $\partial B_\varepsilon$ . As $\varepsilon \to 0$ the first term tends to $-u(\mathbf{r}_0)$ , because $u$ is continuous at $\mathbf{r}_0$ ; the second tends to $0$ , because $\lvert\nabla u\rvert$ is bounded near $\mathbf{r}_0$ and the prefactor is $\varepsilon$ .

Step 4 (the volume integral converges). $\lvert F\rho\rvert \le C/r$ near $\mathbf{r}_0$ with $C = \max\lvert\rho\rvert/4\pi$ , and $\int_{B_\delta} r^{-1} dV = 4\pi\int_0^\delta r\, dr = 2\pi\delta^2$ is finite, so $\int_{V_\varepsilon} F\rho\, dV \to \int_V F\rho\, dV$ . Letting $\varepsilon\to0$ in (2) therefore gives

u(\mathbf{r}_0) = \int_V F\rho\; dV + \oint_S \left( u\,\frac{\partial F}{\partial n} - F\,\frac{\partial u}{\partial n}\right) dS . \tag{3}

This is the representation formula with the bare fundamental solution (Strauss's Eq. (7.2.1), there written for a harmonic $u$ ).

Step 5 (add the two). Adding (1) to (3) and using $G = F + H$ in each of the three places where the pair $F, H$ occurs gives exactly (R).

Step 6 (the two specialisations). If $G = G_D$ vanishes at every point of $S$ , the term $-\oint_S G_D\,\partial u/\partial n\, dS$ in (R) vanishes identically, whatever $\partial u/\partial n$ may be; substituting the Dirichlet datum $u = f$ in the surviving term gives (D). If instead $\partial G_N/\partial n = 1/A$ on $S$ , then

\oint_S u\,\frac{\partial G_N}{\partial n}\; dS = \frac{1}{A}\oint_S u\; dS = \langle u\rangle_S ,

and substituting the Neumann datum $\partial u/\partial n = f$ in the other term gives (N). $\square$

Theorem(Symmetry of the Dirichlet Green's function)AdaptedAdapted from Strauss (2008), §7.3, Eq. (4) and its proof, pp. 189-190; Evans (2010), §2.2.4, Theorem 13, pp. 34-35; Arfken, Weber & Harris (2013), §10.2, Eqs. (10.34)-(10.36), p. 460

Let $V$ be a bounded domain with $C^1$ boundary and suppose a Dirichlet Green's function $G_D(\cdot\,,\mathbf{r}_0)$ exists for every $\mathbf{r}_0 \in V$ . Then

G_D(\mathbf{a},\mathbf{b}) = G_D(\mathbf{b},\mathbf{a}) \qquad\text{for all } \mathbf{a},\mathbf{b}\in V,\ \mathbf{a}\neq\mathbf{b} .

In electrostatics $G_D(\mathbf{r},\mathbf{r}_0)$ is (up to the factor $-1/\varepsilon_0$ ) the potential inside a grounded conductor $S$ due to a unit point charge at $\mathbf{r}_0$ , and the symmetry is the principle of reciprocity: a source at $\mathbf{a}$ produces at $\mathbf{b}$ what a source at $\mathbf{b}$ would produce at $\mathbf{a}$ . Arfken, Weber and Harris obtain the same conclusion, in the form $G(\mathbf{r}_1,\mathbf{r}_2) = G^*(\mathbf{r}_2,\mathbf{r}_1)$ , from the Hermiticity of the operator together with its boundary conditions, and note that it holds for any problem whose definition makes the operator Hermitian. For real $G$ the two statements agree.

Idea of the proofOriginal

Green's second identity is antisymmetric in its two slots, so if we feed it two Green's functions with poles at $\mathbf{a}$ and $\mathbf{b}$ , the left-hand side is the difference of two numbers and the right-hand side vanishes, because both functions are zero on $S$ . The two numbers turn out to be $G_D(\mathbf{a},\mathbf{b})$ and $G_D(\mathbf{b},\mathbf{a})$ . Formally: $\int_V(v\,\nabla^2 w - w\,\nabla^2 v)\,dV$ with $\nabla^2 v = \delta_{\mathbf{a}}$ and $\nabla^2 w = \delta_{\mathbf{b}}$ picks out $v(\mathbf{b}) - w(\mathbf{a})$ , which must equal zero.

The work lies in making this legitimate: both functions are singular, so two balls have to be cut out, and the limit on each small sphere has to be computed. Each sphere reproduces the calculation of Step 3 of the previous proof, with the other Green's function playing the role that $u$ played there.

ProofAdaptedAdapted from Strauss (2008), §7.3, pp. 189-190; Evans (2010), §2.2.4, proof of Theorem 13, pp. 34-35

Fix distinct $\mathbf{a},\mathbf{b}\in V$ and write

v(\mathbf{r}) = G_D(\mathbf{r},\mathbf{a}), \qquad w(\mathbf{r}) = G_D(\mathbf{r},\mathbf{b}) .

Let $\varepsilon$ be small enough that the closed balls $\overline{B_\varepsilon(\mathbf{a})}$ and $\overline{B_\varepsilon(\mathbf{b})}$ are disjoint and contained in $V$ , and put $V_\varepsilon = V\setminus\bigl(\overline{B_\varepsilon(\mathbf{a})}\cup \overline{B_\varepsilon(\mathbf{b})}\bigr)$ .

Step 1 (the identity on the punctured region). On $V_\varepsilon$ both $v$ and $w$ are $C^2$ and harmonic, so Green's second identity has vanishing left-hand side:

0 = \oint_{S}\left(v\,\frac{\partial w}{\partial n} - w\,\frac{\partial v}{\partial n}\right) dS + A_\varepsilon + B_\varepsilon ,

where $A_\varepsilon$ and $B_\varepsilon$ are the same integrand over $\partial B_\varepsilon(\mathbf{a})$ and $\partial B_\varepsilon(\mathbf{b})$ , with the normal of $V_\varepsilon$ , which points into each small ball.

Step 2 (the outer boundary drops out). On $S$ both $v$ and $w$ vanish, by the Dirichlet condition. Each of the two products $v\,\partial w/\partial n$ and $w\,\partial v/\partial n$ therefore vanishes pointwise on $S$ , and $A_\varepsilon + B_\varepsilon = 0$ for every admissible $\varepsilon$ .

Step 3 (the sphere about $\mathbf{a}$ ). Near $\mathbf{a}$ write $v = F_{\mathbf{a}} + H_{\mathbf{a}}$ with $F_{\mathbf{a}}(\mathbf{r}) = -1/(4\pi\lvert\mathbf{r}-\mathbf{a}\rvert)$ and $H_{\mathbf{a}}$ of class $C^2$ near $\mathbf{a}$ ; $w$ is also $C^2$ near $\mathbf{a}$ , since its only singularity is at $\mathbf{b}$ . Of the four pieces of $A_\varepsilon$ :

$\oint v\,\partial w/\partial n\, dS$ is $O(\varepsilon^{-1})\cdot O(1)\cdot O(\varepsilon^{2}) = O(\varepsilon)$ , because $\lvert v\rvert \le 1/(4\pi\varepsilon) + O(1)$ while $\partial w/\partial n$ is bounded and the area is $4\pi\varepsilon^2$ ;
$-\oint w\,\partial H_{\mathbf{a}}/\partial n\, dS$ is $O(\varepsilon^{-1})\cdot O(1) \cdot O(\varepsilon^2) = O(\varepsilon)$ for the same reason, with the roles exchanged;
$-\oint w\,\partial F_{\mathbf{a}}/\partial n\, dS$ is the surviving term. As in Step 3 of the previous proof, $\partial F_{\mathbf{a}}/\partial n = -1/(4\pi\varepsilon^2)$ on the sphere, so this equals $+\overline{w}_\varepsilon$ , the average of $w$ over $\partial B_\varepsilon(\mathbf{a})$ , which tends to $w(\mathbf{a})$ by continuity of $w$ at $\mathbf{a}$ .

Hence $A_\varepsilon \to w(\mathbf{a}) = G_D(\mathbf{a},\mathbf{b})$ .

Step 4 (the sphere about $\mathbf{b}$ ). The same count with $v$ and $w$ exchanged applies at $\mathbf{b}$ , but the integrand $v\,\partial w/\partial n - w\,\partial v/\partial n$ is antisymmetric, so the surviving term is now $+\oint v\,\partial F_{\mathbf{b}}/\partial n\, dS \to -v(\mathbf{b}) = -G_D(\mathbf{b},\mathbf{a})$ . Hence $B_\varepsilon \to -G_D(\mathbf{b},\mathbf{a})$ .

Step 5. Letting $\varepsilon\to0$ in $A_\varepsilon + B_\varepsilon = 0$ gives $G_D(\mathbf{a},\mathbf{b}) - G_D(\mathbf{b},\mathbf{a}) = 0$ . $\square$

Theorem(Uniqueness of the Dirichlet Green's function)AdaptedAdapted from Strauss (2008), §7.3, p. 188, where uniqueness is Exercise 1

For each $\mathbf{r}_0\in V$ there is at most one Dirichlet Green's function with pole at $\mathbf{r}_0$ .

ProofOriginal

Let $G_1, G_2$ be two of them and put $D = G_1 - G_2$ . By condition 2 of the definition, $G_i = F + H_i$ with $H_i \in C^2(\overline{V})$ harmonic in $V$ , so $D = H_1 - H_2$ is harmonic on the whole of $V$ — the singularities, both equal to $F$ , have cancelled — and $D \in C^2(\overline{V})$ . On $S$ , $D = 0 - 0 = 0$ . By the uniqueness theorem for the Dirichlet problem (Green's identities, part 1, applied with $f = 0$ and $h = 0$ ), $D \equiv 0$ . $\square$

Note(Existence)AdaptedAdapted from Strauss (2008), §7.3, p. 188; Evans (2010), §2.2.4, p. 34

That a Dirichlet Green's function exists for a bounded domain is asserted but not proved here; Strauss records that it "can be shown", and Evans points out that constructing $G$ amounts to solving the corrector problem $\Delta\phi^x = 0$ in $U$ with $\phi^x = \Phi(y-x)$ on $\partial U$ , which is itself a Dirichlet problem and "in general a difficult matter", possible explicitly only for domains of simple geometry. The theorems above are therefore conditional on existence, and the next unit, the method of images, exhibits $G_D$ for the two geometries where the construction is elementary. For bounded domains the existence proof belongs to the general solvability theory of elliptic equations; see The Lax-Milgram theorem and the Dirichlet problem.

Remark(Where the hypotheses are used, and what the conventions cost)Original

$u \in C^2(\overline{V})$ . Used in Step 2 of the first proof, to apply Green's second identity on $V_\varepsilon$ , and in Step 3, where boundedness of $u$ and $\nabla u$ near $\mathbf{r}_0$ kills the second small-sphere term. Without continuity of $u$ up to $S$ the boundary integrals in (R) need not even be defined.
$\rho$ continuous (or merely bounded) on $\overline{V}$ . Used only in Step 4, to make $\int_V F\rho\, dV$ converge against the $1/r$ singularity of $F$ . Boundedness is more than needed; $\rho$ integrable against $1/r$ suffices for that step.
Boundedness of $V$ . The formulas extend to unbounded domains, but only with an added condition at infinity. Strauss handles the half-space by demanding of every function in play that it, along with its derivatives, decay to $0$ at large $\lvert x\rvert$ ; Riley, Hobson and Bence treat an exterior region formally, taking $S$ to consist of the inner surface together with a sphere of radius $R$ and letting $R\to\infty$ . The content of the condition is that the contribution of the sphere at infinity vanishes.
The Dirichlet condition in Step 2 of the symmetry proof. It is used as $G_D = 0$ on $S$ for both poles. For the Neumann Green's function the same computation leaves the outer boundary term $A^{-1}\oint_S(v - w)\, dS$ , so it gives instead $G_N(\mathbf{a},\mathbf{b}) - G_N(\mathbf{b},\mathbf{a}) = \langle G_N(\cdot\,,\mathbf{b})\rangle_S - \langle G_N(\cdot\,,\mathbf{a})\rangle_S .$ For each pole, a Neumann Green's function is fixed only to within one added constant (adding a constant changes neither $\nabla^2 G$ nor $\partial G/\partial n$ ), so the normalisation $\langle G_N(\cdot\,,\mathbf{r}_0)\rangle_S = 0$ may be imposed for every $\mathbf{r}_0$ , and with it $G_N$ is symmetric too. Without that normalisation it need not be.
The Neumann formula (N) does not depend on the normalisation. Replacing $G_N(\mathbf{r},\mathbf{r}_0)$ by $G_N(\mathbf{r},\mathbf{r}_0) + c(\mathbf{r}_0)$ changes the right-hand side of (N) by $c(\mathbf{r}_0)\left[\int_V \rho\, dV - \oint_S f\, dS\right]$ , which is zero by the compatibility condition that Neumann data must satisfy for a solution to exist (Green's identities). So the ambiguity in $G_N$ is invisible in the answer — as it must be, since the answer is $u(\mathbf{r}_0)$ .
Signs under the other convention. With Evans's $G^{\mathrm E} = -G$ , satisfying $-\Delta G^{\mathrm E} = \delta$ , and his form $-\Delta u = f$ of the equation (so $f = -\rho$ ), formula (D) becomes $u(x) = -\oint g\,\partial G^{\mathrm E}/\partial\nu\, dS + \int f G^{\mathrm E}\, dy$ , which is his Eq. (30): both terms flip sign twice and one overall sign remains in the boundary term. Strauss's Theorem 2 reads $u(x_0) = \oint h\, \partial G/\partial n\, dS + \int f G\, dx$ for $\Delta u = f$ , which is our (D) with his $G$ equal to ours. The three statements are the same statement.
Dimension. Nothing above is special to $\mathbb{R}^3$ except the explicit $F$ . In $\mathbb{R}^2$ the fundamental solution is $F = (2\pi)^{-1}\ln\lvert\mathbf{r}-\mathbf{r}_0\rvert$ (plus an arbitrary constant, which here matters: Riley, Hobson and Bence note that in two dimensions one cannot impose $F\to0$ at infinity), and the small-circle computation of Step 3 produces $u(\mathbf{r}_0)$ in just the same way. Hassani states (D) for a general dimension $m$ , and Evans works throughout in $\mathbb{R}^n$ with $\Phi$ the $n$ -dimensional fundamental solution.

Note(Radial integrals used in the examples)Original

On the ball $B_a$ of radius $a$ about the origin and its boundary sphere $S_a$ , for a function of $r = \lvert\mathbf{r}\rvert$ alone,

\int_{B_a} g(r)\, dV = 4\pi\int_0^a g(r)\, r^2\, dr , \qquad \oint_{S_a} g(a)\, dS = 4\pi a^2 g(a) .

In particular $\int_{B_a} r^{-1}\, dV = 4\pi\int_0^a r\, dr = 2\pi a^2$ and $\int_{B_a} dV = \tfrac{4}{3}\pi a^3$ .

Example(The Green's functions of a ball with pole at the centre, and a test of (D))Original

Take $V = B_a$ , the ball of radius $a$ about the origin, and put the pole at $\mathbf{r}_0 = \mathbf{0}$ . Then $F = -1/(4\pi r)$ is already radial, and the corrector need only be a constant, so both Green's functions can be written down at once.

Dirichlet. Try $G_D(\mathbf{r},\mathbf{0}) = -\dfrac{1}{4\pi r} + c$ . A constant is harmonic and $C^2$ on $\overline{B_a}$ , so condition 2 holds; $G_D = 0$ at $r = a$ forces $c = 1/(4\pi a)$ , giving $G_D(\mathbf{r},\mathbf{0}) = \frac{1}{4\pi}\left(\frac{1}{a} - \frac{1}{r}\right), \qquad \frac{\partial G_D}{\partial n}\bigg|_{r=a} = \frac{\partial}{\partial r}\left(-\frac{1}{4\pi r}\right)\bigg|_{r=a} = \frac{1}{4\pi a^2} .$
Neumann. For $G_N = -1/(4\pi r) + c$ the normal derivative at $r = a$ is the same $1/(4\pi a^2)$ , which is precisely $1/A$ with $A = 4\pi a^2$ . So the bare fundamental solution is a Neumann Green's function of the ball when the pole sits at the centre, for any $c$ ; take $c = 0$ .

Test of (D). Let $u = r^2 = x^2+y^2+z^2$ on $\overline{B_a}$ . It is a polynomial, hence in $C^2(\overline{B_a})$ , with $\rho = \nabla^2 u = 6$ and boundary value $f = a^2$ . Formula (D) should return $u(\mathbf{0}) = 0$ .

Boundary term. $\displaystyle \oint_{S_a} f\,\frac{\partial G_D}{\partial n}\, dS = a^2\cdot\frac{1}{4\pi a^2}\cdot 4\pi a^2 = a^2 .$
Volume term. $\displaystyle \int_{B_a} G_D\,\rho\; dV = 6\cdot\frac{1}{4\pi}\left[\frac{1}{a}\int_{B_a} dV - \int_{B_a}\frac{dV}{r}\right] = \frac{6}{4\pi}\left[\frac{1}{a}\cdot\frac{4\pi a^3}{3} - 2\pi a^2\right] = \frac{6}{4\pi}\left[\frac{4\pi a^2}{3} - 2\pi a^2\right] = -a^2 .$
Sum. $a^2 - a^2 = 0 = u(\mathbf{0})$ , as required. Note that both terms are of size $a^2$ and cancel exactly; with $a = 2$ they are $+4$ and $-4$ .

Test of (N). The same $u = r^2$ has $f = \partial u/\partial n = 2a$ on $S_a$ and $\langle u\rangle_{S_a} = a^2$ . Formula (N) with $G_N = -1/(4\pi r)$ gives

\underbrace{6\cdot\left(-\frac{1}{4\pi}\right)\int_{B_a}\frac{dV}{r}}_{=\;-3a^2} \;+\;\underbrace{a^2}_{\langle u\rangle_S} \;-\;\underbrace{\left(-\frac{1}{4\pi a}\right)(2a)\,4\pi a^2}_{=\;-2a^2} \;=\; -3a^2 + a^2 + 2a^2 = 0 = u(\mathbf{0}) .

The data pass the compatibility test, $\oint_{S_a} f\, dS = 2a\cdot4\pi a^2 = 8\pi a^3$ and $\int_{B_a}\rho\, dV = 6\cdot\tfrac43\pi a^3 = 8\pi a^3$ , as they must for a solution to exist.

Example(Symmetry made explicit, and a non-zero value)Original

Stay with the unit ball, $a = 1$ . The Dirichlet Green's function with pole at a general interior point is derived in the next unit; for the pole at the centre we have $G_D(\mathbf{r},\mathbf{0}) = (4\pi)^{-1}(1 - 1/r)$ from the previous example, and the symmetry theorem then predicts

G_D(\mathbf{0},\mathbf{b}) = G_D(\mathbf{b},\mathbf{0}) = \frac{1}{4\pi}\left(1 - \frac{1}{\lvert\mathbf{b}\rvert}\right) .

So the potential at the centre of a grounded unit sphere due to a unit source at distance $\lvert\mathbf{b}\rvert = \tfrac12$ is $(4\pi)^{-1}(1-2) = -1/(4\pi) = -0.0795775$ , and this equals the potential at distance $\tfrac12$ from the centre due to a unit source at the centre — which is immediate from the explicit formula. The content of the theorem is that the first of those two numbers, which would otherwise require the Green's function with an off-centre pole, is obtainable from the second.

For a case in which (D) returns a non-zero value, take $u = r^2 + 1$ on the unit ball: then $\rho = 6$ , $f = 2$ , and (D) gives $2\cdot\tfrac{1}{4\pi}\cdot 4\pi + (-1) = 2 - 1 = 1 = u(\mathbf{0})$ .

ApplicationOriginal

Electrostatics with conductors. The reason the Green's function method is the standard tool of potential theory is that it separates the geometry from the data. The potential inside a grounded conducting surface $S$ obeys $\nabla^2\phi = -\rho_{\mathrm e}/\varepsilon_0$ with $\phi = 0$ on $S$ , so by (D),

\phi(\mathbf{r}_0) = -\frac{1}{\varepsilon_0}\int_V G_D(\mathbf{r},\mathbf{r}_0)\,\rho_{\mathrm e}(\mathbf{r})\; dV .

One solves the geometry once, for $G_D$ , and then any charge distribution inside is handled by an integration. If the conductor is instead held at a prescribed potential $f$ that varies over the surface, the boundary term of (D) supplies the rest. The reciprocity theorem $G_D(\mathbf{a},\mathbf{b}) = G_D(\mathbf{b},\mathbf{a})$ is used in reverse in measurement: a four-point resistivity probe gives the same reading when source and sense contacts are interchanged, and in capacitance and antenna work the same symmetry underlies the reciprocity of transmitting and receiving patterns.

Steady heat flow and groundwater. $\nabla^2 T = -q/k$ with $T$ prescribed on the boundary is the steady heat equation with volumetric heating $q$ and conductivity $k$ ; $G_D$ of the region is the temperature response to a point heater with the boundary held at zero, and (D) is the superposition of such responses. The Neumann form (N) is the one needed when it is the flux that is prescribed — an insulated or uniformly heated surface, or in hydrology a confined aquifer with prescribed recharge across its boundary — and the fact that only $\langle u \rangle_S$ is left free is the statement that steady flux data fix the temperature field but not its datum level.

Lattice sums and image systems. Because $G_D$ for a region and $F$ for free space differ by a harmonic corrector, the practical content of the theory is the search for correctors. For half-spaces, slabs, wedges of rational angle and spheres they are finite or infinite sets of image sources (The method of images); for a general region one falls back on an eigenfunction expansion or on numerical solution, and in boundary-element codes the discretised form of (R) — unknowns on the surface only — is the whole method.

Why the formula is more than a solution formula. Arfken, Weber and Harris stress that a Green's function turns a differential operator into an integral operator, and integral operators are bounded while differential operators are not. That trade is what makes the Green's function the starting point for the perturbative treatment of the Helmholtz and Schrödinger equations (Green's functions for the Helmholtz equation, Scattering and the Born approximation) and for the Fredholm theory of integral equations.

Method

The method of images

Adapted

Add to PDF

Builds the corrector of a Dirichlet Green's function out of point sources placed outside the region, and carries the construction through for a half-space (one image of strength $-1$ ) and for a sphere of radius $a$ (one image at the inverse point $a^2\mathbf{r}_0/r_0^2$ of strength $-a/r_0$ ), recovering the Poisson kernels of both.

Prerequisites: Green's functions for Laplace and Poisson problems, Fundamental solution of Laplace's equation

Note(What the method does)Original

By the previous unit, finding a Dirichlet Green's function means finding a corrector: a function $H$ , harmonic throughout $V$ and continuous up to $S$ , such that $F + H$ vanishes on $S$ , where $F$ is the fundamental solution with pole at $\mathbf{r}_0$ . The method of images produces $H$ from a stock of functions that are automatically harmonic inside $V$ : fundamental solutions with poles outside $V$ . Nothing else about them is needed. The only question is whether finitely many of them, with suitable positions and strengths, can be made to cancel $F$ on $S$ — and for a plane boundary and a spherical boundary the answer is yes, with a single image in each case.

Method(Constructing a Green's function by images)AdaptedAdapted from Riley, Hobson & Bence (2006), §21.5.3, steps (i)-(iv), pp. 758-759

To solve $\nabla^2 u = \rho$ in a region $V$ with $u = f$ given on $S = \partial V$ :

Place images. Choose points $\mathbf{r}_1,\dots,\mathbf{r}_N$ outside $V$ and strengths $q_1,\dots,q_N$ , to be determined, and think of them as extra point sources $\sum_n q_n\delta(\mathbf{r}-\mathbf{r}_n)$ added to the unit source at $\mathbf{r}_0\in V$ .
Form the candidate. Set $G(\mathbf{r},\mathbf{r}_0) = F(\mathbf{r},\mathbf{r}_0) + \sum_{n=1}^N q_n\, F(\mathbf{r},\mathbf{r}_n), \qquad F(\mathbf{r},\mathbf{r}') = -\frac{1}{4\pi\lvert\mathbf{r}-\mathbf{r}'\rvert}$ in three dimensions, or $F = (2\pi)^{-1}\ln\lvert\mathbf{r}-\mathbf{r}'\rvert$ in two, where an additive constant is also at one's disposal. Each image term is harmonic on $V$ , because its pole is outside $V$ ; so the sum satisfies $\nabla^2 G = \delta(\mathbf{r}-\mathbf{r}_0)$ in $V$ and has the required singularity there, whatever the $q_n$ and $\mathbf{r}_n$ turn out to be.
Impose the boundary condition. Fix the positions and strengths by requiring $G(\mathbf{r},\mathbf{r}_0) = 0$ at every $\mathbf{r}\in S$ (for a Neumann Green's function, require $\partial G/\partial n = 1/A$ instead). If $V$ is unbounded, require in addition that $G\to0$ in the directions in which $V$ extends to infinity.
Integrate. Substitute into the representation formula (D) or (N) of Green's functions for Laplace and Poisson problems.

For a boundary made of planes (in three dimensions) or straight lines (in two), step 3 needs no calculation: the image positions are the geometrical mirror images of $\mathbf{r}_0$ in those planes, and the strengths alternate in sign, one image for a single plane, three for a right-angled corner, and an infinite row for a slab. Riley, Hobson and Bence stress that for a general curved boundary the positions and strengths cannot be found at all, and that the method is confined to a short list of geometries.

Method(Green's function of a half-space)AdaptedAdapted from Strauss (2008), §7.4, Eq. (1) and the three verifications, pp. 191-192; Riley, Hobson & Bence (2006), §21.5.3, Eq. (21.94), p. 759; Hassani (2013), §22.1, Example 22.1.1, pp. 667-668; Evans (2010), §2.2.4b, the definition on p. 35

Let $V = \{\mathbf{r} = (x,y,z) : z > 0\}$ and for $\mathbf{r}_0 = (x_0,y_0,z_0)\in V$ let

\mathbf{r}_0^* = (x_0,y_0,-z_0)

be its mirror image in the plane $z = 0$ . Then

G_D(\mathbf{r},\mathbf{r}_0) = -\frac{1}{4\pi\lvert\mathbf{r}-\mathbf{r}_0\rvert} + \frac{1}{4\pi\lvert\mathbf{r}-\mathbf{r}_0^*\rvert}

is a Dirichlet Green's function for $\nabla^2$ on $V$ : it is harmonic on $V\setminus\{\mathbf{r}_0\}$ , differs from $F(\mathbf{r},\mathbf{r}_0)$ by a function harmonic on $V$ , vanishes at every point of the plane $z=0$ , and tends to $0$ as $\lvert\mathbf{r}\rvert\to\infty$ . Its normal derivative on the boundary is

\frac{\partial G_D}{\partial n}\bigg|_{z=0} = \frac{z_0}{2\pi\bigl[(x-x_0)^2+(y-y_0)^2+z_0^2\bigr]^{3/2}} ,

so the solution of $\nabla^2 u = 0$ in $z>0$ with $u(x,y,0) = f(x,y)$ is

u(x_0,y_0,z_0) = \frac{z_0}{2\pi}\int_{-\infty}^{\infty}\!\!\int_{-\infty}^{\infty} \frac{f(x,y)\; dx\, dy}{\bigl[(x-x_0)^2+(y-y_0)^2+z_0^2\bigr]^{3/2}} . \tag{H}

The kernel in (H) is the Poisson kernel of the half-space. Strauss derives (H) in this way; Evans, who writes it in $\mathbb{R}^n$ and calls the construction "reflecting the singularity", proves separately (his Theorem 14) that for $f$ continuous and bounded the formula really does define a bounded harmonic function in $z>0$ that attains the boundary values $f$ in the limit — the Green's function derivation presupposes a solution and so cannot establish existence by itself. Hassani notes that the whole construction is the electrostatics of two opposite point charges placed symmetrically about the plane, which is why their potential vanishes on it.

Idea of the proofOriginal

A plane is a mirror, and the potential of a charge plus an equal and opposite charge at the mirror point vanishes on the mirror, because every point of the plane is equidistant from the two. That one sentence is the whole construction; everything else is checking the bookkeeping, namely that the second charge sits outside the region (so that it adds nothing to the Laplacian there) and that its singularity therefore does not disturb the required behaviour at $\mathbf{r}_0$ .

Hassani makes the same point from the other direction, without guessing: the corrector must be harmonic and, by the symmetry of the problem in $x$ and $y$ and the evenness of $\nabla^2$ in $z$ , must be a fundamental solution whose pole has the same $x_0,y_0$ and $z$ -dependence through either $(z-z_0)^2$ or $(z+z_0)^2$ . The first choice gives $G\equiv0$ , so the second is forced.

ProofAdaptedAdapted from Strauss (2008), §7.4, the verification of properties (i)-(iii) and the computation of $-\partial G/\partial z$, pp. 191-192; Riley, Hobson & Bence (2006), §21.5.3, p. 760

Step 1 (the corrector is harmonic in $V$ ). The function $\mathbf{r}\mapsto (4\pi\lvert\mathbf{r}-\mathbf{r}_0^*\rvert)^{-1}$ is smooth and harmonic except at $\mathbf{r}_0^*$ , which has $z = -z_0 < 0$ and so lies outside $V$ . Hence $H = G_D - F$ is $C^\infty$ and harmonic on a neighbourhood of $\overline V$ , which is condition 2 of the definition of a Green's function; in particular $G_D$ keeps exactly the singularity of $F$ at $\mathbf{r}_0$ .

Step 2 (it vanishes on the boundary). Let $\mathbf{r} = (x,y,0)$ . Then

\lvert\mathbf{r}-\mathbf{r}_0\rvert^2 = (x-x_0)^2+(y-y_0)^2+z_0^2 = \lvert\mathbf{r}-\mathbf{r}_0^*\rvert^2 ,

the two differing only in the sign of $z_0$ , which is squared. So the two terms of $G_D$ cancel and $G_D = 0$ on $z=0$ . The same computation shows that the cancellation fails at every other point, as it must.

Step 3 (behaviour at infinity). Both terms are $O(\lvert\mathbf{r}\rvert^{-1})$ , so $G_D\to0$ ; this is the extra condition that an unbounded region needs, and it is what makes the sphere at infinity contribute nothing to the representation formula.

Step 4 (the normal derivative). Using $\nabla\lvert\mathbf{r}-\mathbf{r}'\rvert = (\mathbf{r}-\mathbf{r}')/\lvert\mathbf{r}-\mathbf{r}'\rvert$ ,

\nabla G_D = \frac{\mathbf{r}-\mathbf{r}_0}{4\pi\lvert\mathbf{r}-\mathbf{r}_0\rvert^3} - \frac{\mathbf{r}-\mathbf{r}_0^*}{4\pi\lvert\mathbf{r}-\mathbf{r}_0^*\rvert^3} .

On $z=0$ the two denominators are equal, call the common value $4\pi\rho^3$ , while the $z$ -components of the numerators are $z-z_0 = -z_0$ and $z+z_0 = +z_0$ . Hence $\partial G_D/\partial z = -2z_0/(4\pi\rho^3)$ there. The outward normal of $\{z>0\}$ is $-\hat{\mathbf{z}}$ , so $\partial G_D/\partial n = -\partial G_D/\partial z = z_0/(2\pi\rho^3)$ , which is the formula stated.

Step 5. Substituting into formula (D) of Green's functions for Laplace and Poisson problems with $\rho \equiv 0$ , and noting that the part of $S$ at infinity contributes nothing by Step 3, gives (H). $\square$

Method(Green's function of a ball or of the region outside a sphere)AdaptedAdapted from Strauss (2008), §7.4, Eqs. (5)-(8), (11) and (13)-(16), pp. 193-195; Riley, Hobson & Bence (2006), §21.5.3, Eqs. (21.97)-(21.99) and p. 764, pp. 762-764; Hassani (2013), §22.1, Example 22.1.2 and Eq. (22.6), pp. 668-669

Fix $a > 0$ . For $\mathbf{r}_0\neq\mathbf{0}$ with $r_0 = \lvert\mathbf{r}_0\rvert$ , the inverse point of $\mathbf{r}_0$ in the sphere $\lvert\mathbf{r}\rvert = a$ is

\mathbf{r}_1 = \frac{a^2}{r_0^2}\,\mathbf{r}_0 ,

the point on the same ray from the origin with $r_0\lvert\mathbf{r}_1\rvert = a^2$ . Write $\rho = \lvert\mathbf{r}-\mathbf{r}_0\rvert$ and $\rho^* = \lvert\mathbf{r}-\mathbf{r}_1\rvert$ . Then:

The geometrical identity. For every $\mathbf{r}$ with $\lvert\mathbf{r}\rvert = a$ , $\rho^* = \frac{a}{r_0}\,\rho . \tag{G}$
The Green's function. For $0 < r_0 < a$ , $G_D(\mathbf{r},\mathbf{r}_0) = -\frac{1}{4\pi\rho} + \frac{a/r_0}{4\pi\rho^*}$ is a Dirichlet Green's function for the ball $V = \{\lvert\mathbf{r}\rvert<a\}$ . The same formula, read with $r_0 > a$ (so that the image lies inside the sphere), is a Dirichlet Green's function for the exterior region $\{\lvert\mathbf{r}\rvert>a\}$ , which then also satisfies $G_D\to0$ at infinity. For the centred case $\mathbf{r}_0 = \mathbf{0}$ , where the inverse point is undefined, the Green's function of the ball is $G_D(\mathbf{r},\mathbf{0}) = (4\pi)^{-1}(a^{-1}-r^{-1})$ .
The Poisson kernel. On $\lvert\mathbf{r}\rvert = a$ , $\frac{\partial G_D}{\partial n} = \frac{a^2-r_0^2}{4\pi a\,\rho^3} ,$ so the solution of $\nabla^2 u = 0$ in $\lvert\mathbf{r}\rvert<a$ with $u = f$ on $\lvert\mathbf{r}\rvert = a$ is $u(\mathbf{r}_0) = \frac{a^2-r_0^2}{4\pi a}\oint_{\lvert\mathbf{r}\rvert=a}\frac{f(\mathbf{r})}{\lvert\mathbf{r}-\mathbf{r}_0\rvert^3}\; dS(\mathbf{r}) . \tag{P}$

Formula (P) is Poisson's integral formula in three dimensions, obtained in Poisson's integral formula by a different route; Evans asserts (his Theorem 15) that for continuous $f$ it does define a smooth harmonic function in the ball attaining the boundary values $f$ in the limit, remarking that the proof is like that of his Theorem 14 and leaving it as an exercise. Strauss notes that for $\mathbf{r}_0 = \mathbf{0}$ formula (P) reduces to $u(\mathbf{0}) = \langle f\rangle$ , the mean value property of harmonic functions (Harmonic functions).

Idea of the proofOriginal

A sphere is not a mirror, so the strength of the image cannot be $-1$ : no pair of equal and opposite point sources has a sphere as an equipotential. What saves the construction is that the sphere is an equipotential of a pair of unequal opposite sources, and the required ratio of distances is supplied by elementary geometry: inversion in the sphere, $\mathbf{r}_0\mapsto a^2\mathbf{r}_0/r_0^2$ , is exactly the map for which the ratio $\rho^*/\rho$ of distances from the two points is the same constant $a/r_0$ at every point of the sphere. Dividing the image strength by that constant then makes the two terms cancel there.

One can also run the argument backwards and derive the position and strength, as Hassani does: write the image as an unknown strength $k$ at an unknown point on the ray through $\mathbf{r}_0$ , demand that $k\rho = \rho^*$ hold identically in the angle between $\mathbf{r}$ and $\mathbf{r}_0$ , and compare coefficients. Two equations emerge, with two roots; the root $k=1$ puts the image on top of the source and is discarded, leaving $k = a/r_0$ . Riley, Hobson and Bence shortcut the comparison by imposing the ratio condition at just the two points where the sphere meets the line through the centre and $\mathbf{r}_0$ , which already fixes $\lvert\mathbf{r}_1\rvert = a^2/r_0$ , and then verifying that it holds everywhere.

ProofAdaptedAdapted from Riley, Hobson & Bence (2006), §21.5.3, Eqs. (21.98)-(21.99) and pp. 762-764; Strauss (2008), §7.4, Eqs. (7)-(8) and (13)-(15), pp. 193-195; Evans (2010), §2.2.4c, the identity on p. 38

Step 1 (proof of (G)). Expand, using $\lvert\mathbf{r}\rvert^2 = a^2$ and $\lvert\mathbf{r}_1\rvert^2 = a^4/r_0^2$ :

\rho^{*2} = \lvert\mathbf{r}\rvert^2 - 2\,\mathbf{r}\cdot\mathbf{r}_1 + \lvert\mathbf{r}_1\rvert^2 = a^2 - \frac{2a^2}{r_0^2}\,\mathbf{r}\cdot\mathbf{r}_0 + \frac{a^4}{r_0^2} = \frac{a^2}{r_0^2}\Bigl( r_0^2 - 2\,\mathbf{r}\cdot\mathbf{r}_0 + a^2\Bigr) .

Now use $a^2 = \lvert\mathbf{r}\rvert^2$ once more in the last bracket: it becomes $\lvert\mathbf{r}\rvert^2 - 2\mathbf{r}\cdot\mathbf{r}_0+r_0^2 = \rho^2$ . Taking square roots gives (G). The step is reversible only on the sphere: it is the single place where $\lvert\mathbf{r}\rvert = a$ is used, and off the sphere the ratio $\rho^*/\rho$ is not constant.

Step 2 (the corrector is harmonic in $V$ ). For the ball, $\lvert\mathbf{r}_1\rvert = a^2/r_0 > a$ because $r_0 < a$ , so the image point lies outside $\overline V$ and the image term is smooth and harmonic on a neighbourhood of $\overline V$ . For the exterior region, $r_0 > a$ gives $\lvert\mathbf{r}_1\rvert < a$ , again outside the region in question. In both cases condition 2 of the definition holds, and $\nabla^2 G_D = \delta(\mathbf{r}-\mathbf{r}_0)$ in $V$ .

Step 3 (it vanishes on the sphere). At $\lvert\mathbf{r}\rvert = a$ , by (G),

G_D = -\frac{1}{4\pi\rho} + \frac{a/r_0}{4\pi\,(a/r_0)\rho} = -\frac{1}{4\pi\rho} + \frac{1}{4\pi\rho} = 0 .

This is where the strength $a/r_0$ is used, and it is the only reason for that particular value.

Step 4 (the centred case). If $\mathbf{r}_0 = \mathbf{0}$ the fundamental solution $-1/(4\pi r)$ is already constant on the sphere, so a constant corrector suffices; $G_D = (4\pi)^{-1}(a^{-1}-r^{-1})$ vanishes at $r=a$ and is the Green's function, as computed in Green's functions for Laplace and Poisson problems.

Step 5 (the normal derivative on the sphere). Differentiating $G_D$ with respect to $\mathbf{r}$ ,

\nabla G_D = \frac{\mathbf{r}-\mathbf{r}_0}{4\pi\rho^3} - \frac{a}{r_0}\,\frac{\mathbf{r}-\mathbf{r}_1}{4\pi\rho^{*3}} .

Restrict to $\lvert\mathbf{r}\rvert = a$ and replace $\rho^{*3}$ by $(a/r_0)^3\rho^3$ using (G); the coefficient of the second term becomes $\frac{a}{r_0}\bigl(\frac{r_0}{a}\bigr)^3 = \frac{r_0^2}{a^2}$ , so

\nabla G_D = \frac{1}{4\pi\rho^3}\left[ (\mathbf{r}-\mathbf{r}_0) - \frac{r_0^2}{a^2}\left(\mathbf{r} - \frac{a^2}{r_0^2}\mathbf{r}_0\right)\right] = \frac{1}{4\pi\rho^3}\left(1 - \frac{r_0^2}{a^2}\right)\mathbf{r} .

The $\mathbf{r}_0$ terms have cancelled, which is the small miracle that makes the kernel so simple: on the sphere, $\nabla G_D$ is radial. The outward normal there is $\hat{\mathbf{n}} = \mathbf{r}/a$ , so

\frac{\partial G_D}{\partial n} = \frac{\mathbf{r}}{a}\cdot\nabla G_D = \frac{1}{4\pi\rho^3}\left(1-\frac{r_0^2}{a^2}\right)\frac{\lvert\mathbf{r}\rvert^2}{a} = \frac{a^2-r_0^2}{4\pi a\rho^3} .

Step 6. Substituting into formula (D) of Green's functions for Laplace and Poisson problems with $\rho\equiv0$ gives (P). $\square$

Remark(Scope and variants)Original

Where the hypotheses bite. Step 1 uses $\lvert\mathbf{r}\rvert=a$ and nothing else; Step 2 uses $r_0\neq a$ , so that source and image are on opposite sides of the sphere; Step 3 uses both. If $\mathbf{r}_0$ is placed on the sphere the construction degenerates, which is as it should be, since a Green's function has its pole in the interior.
The same formula twice. As Riley, Hobson and Bence observe, the interior and exterior problems for a sphere are covered by one formula, read with $r_0<a$ or $r_0>a$ . Arfken, Weber and Harris put the interior case as an exercise in electrostatics with the image charge $q' = -qr_0/a$ at $r_0^2/a$ from the centre, where their $r_0$ is the sphere radius and their $a$ the distance of the charge; that is our $-a/r_0$ at $a^2/r_0$ with the letters exchanged, and the two agree.
Two dimensions. The image sits at the same inverse point, but the strength is $-1$ , not $-a/r_0$ , because the logarithm turns the ratio (G) into an additive constant rather than a factor. For the disc $\lvert\mathbf{r}\rvert<a$ Riley, Hobson and Bence obtain $G_D(\mathbf{r},\mathbf{r}_0) = \frac{1}{2\pi}\left[\ln\lvert\mathbf{r}-\mathbf{r}_0\rvert - \ln\left\lvert\mathbf{r}-\frac{a^2}{r_0^2}\mathbf{r}_0\right\rvert - \ln\frac{r_0}{a}\right],$ where the last term is the additive constant that two-dimensional fundamental solutions carry, chosen so that (G) makes the bracket vanish on $\lvert\mathbf{r}\rvert = a$ ; the normal derivative is $(a^2-r_0^2)/(2\pi a\lvert\mathbf{r}-\mathbf{r}_0\rvert^2)$ and (D) then gives the classical two-dimensional Poisson formula $u(\rho_0,\phi_0) = \frac{1}{2\pi}\int_0^{2\pi}\frac{(a^2-\rho_0^2)\,f(\phi)\; d\phi}{a^2+\rho_0^2-2a\rho_0\cos(\phi-\phi_0)} .$ Strauss reaches the same Green's function and remarks that it reproduces the Poisson formula he had found earlier by separation of variables.
Neumann by images. For the disc, Riley, Hobson and Bence leave the image strength as a parameter $q$ , compute $\partial G/\partial\rho$ at $\rho=a$ , and find that $q = +1$ — the image with the same sign as the source — makes it equal to $1/(2\pi a) = 1/L$ , which is the Neumann condition. So the Neumann and Dirichlet Green's functions of a disc use the same image point and differ only in the sign of the image.
More images, and when they run out. A right-angled corner in two dimensions needs three images: the reflections of $\mathbf{r}_0$ in each of the two lines, both of strength $-1$ , and the reflection in the corner (the image of either of those in the other line), of strength $+1$ . This is the arrangement in Riley, Hobson and Bence's figure (Riley, Hobson & Bence (2006), Fig. 21.13, p. 762) for the quarter-plane $x>0$ , $y>0$ with both axes earthed. A wedge of angle $\pi/m$ needs $2m-1$ ; a slab needs an infinite row, and a half-ball is handled by reflecting the ball's Green's function in the flat face. A region with no such symmetry has no image system, and one must turn to eigenfunction expansions or to numerical solution of the corrector problem.
What the method does not prove. All of the above presupposes that the boundary-value problem has a solution of class $C^2(\overline V)$ , since that is the hypothesis of the representation formula. Evans's Theorems 14 and 15 supply the missing half, by asserting directly that (H) and (P) define harmonic functions with the right boundary limits; he proves Theorem 14, the half-space case, and leaves Theorem 15, the ball, as an exercise on the ground that its proof is similar.

Example(A charged patch on a grounded plane)Original

Solve $\nabla^2 u = 0$ in the half-space $z>0$ with boundary data

f(x,y) = \begin{cases} 1, & x^2+y^2 < 1,\\ 0, & x^2+y^2 > 1,\end{cases}

and evaluate $u$ on the axis, at $\mathbf{r}_0 = (0,0,z_0)$ .

Set up. With $x_0 = y_0 = 0$ , formula (H) becomes an integral over the unit disc in polar coordinates $x = \varrho\cos\phi$ , $y = \varrho\sin\phi$ , with $dx\,dy = \varrho\,d\varrho\,d\phi$ and $(x-x_0)^2+(y-y_0)^2 = \varrho^2$ : $u(0,0,z_0) = \frac{z_0}{2\pi}\int_0^{2\pi}\!\!\int_0^1 \frac{\varrho\; d\varrho\, d\phi}{(\varrho^2+z_0^2)^{3/2}} = z_0\int_0^1\frac{\varrho\; d\varrho}{(\varrho^2+z_0^2)^{3/2}} .$
Integrate. The substitution $s = \varrho^2+z_0^2$ , $ds = 2\varrho\,d\varrho$ , turns this into $\tfrac12 z_0\int_{z_0^2}^{1+z_0^2}s^{-3/2}\,ds = z_0\bigl[-s^{-1/2}\bigr]_{z_0^2}^{1+z_0^2}$ , so $u(0,0,z_0) = 1 - \frac{z_0}{\sqrt{1+z_0^2}} .$
The number. At $z_0 = 1$ , $u = 1 - 1/\sqrt{2} = 0.2928932$ to seven figures.
Three checks. As $z_0\to0^+$ , $u\to1$ , which is the boundary value at the centre of the patch. As $z_0\to\infty$ , $u = 1-(1+z_0^{-2})^{-1/2} = \tfrac12 z_0^{-2}+O(z_0^{-4})$ , the far field of a dipole-like source, decaying faster than a point charge because the total "charge" of a Dirichlet problem is not free. And the answer equals $\Omega/2\pi$ , where $\Omega = 2\pi(1 - z_0/\sqrt{1+z_0^2})$ is the solid angle that the unit disc subtends at $\mathbf{r}_0$ — a general feature of (H) for data that are the indicator function of a region, since the Poisson kernel integrated over a patch is $1/2\pi$ times the solid angle it subtends.
A range check. The one-dimensional restriction $z_0\mapsto 1 - z_0(1+z_0^2)^{-1/2}$ is of course not harmonic by itself; the transverse curvature that balances its axial curvature is invisible in a formula evaluated only on the axis. A cheap check of the closed form instead: the Poisson kernel of (H) is normalised, $\int_{\mathbb{R}^2}\frac{z_0\; dx\, dy}{2\pi\bigl[(x-x_0)^2+(y-y_0)^2+z_0^2\bigr]^{3/2}} = 1 \qquad\text{for every } z_0 > 0 ,$ which is the integral of step 2 with the upper limit $\infty$ instead of $1$ , and indeed gives $z_0\cdot(1/z_0) = 1$ . Equivalently $\int_{z=0}\partial G_D/\partial n\, dS = 1$ , by the formula for $\partial G_D/\partial n$ in the statement. So (H) applied to $f\equiv1$ returns $u\equiv1$ , as it must, since the constant function $1$ is the solution of that Dirichlet problem. The kernel is also positive, so (H) carries data with $0\le f\le1$ to solutions with $0\le u\le1$ , and the value $0.2928932$ is in range.

Example(The Poisson kernel of the unit ball, with a value one can check)Original

Take $a = 1$ and the boundary datum $f = \cos^2\theta$ on the unit sphere, where $\theta$ is the polar angle from the $z$ -axis. Evaluate (P) at the interior point $\mathbf{r}_0 = (0,0,\tfrac12)$ , so $r_0 = \tfrac12$ .

Reduce to one integral. With $\mathbf{r}_0$ on the axis, $\rho^2 = \lvert\mathbf{r}-\mathbf{r}_0\rvert^2 = 1 + r_0^2 - 2r_0\cos\theta$ , which depends on $\theta$ only; using $dS = \sin\theta\, d\theta\, d\phi$ on the unit sphere and doing the $\phi$ integral, $u(\mathbf{r}_0) = \frac{1-r_0^2}{4\pi}\cdot 2\pi\int_0^\pi\frac{\cos^2\theta\,\sin\theta\; d\theta}{(1+r_0^2-2r_0\cos\theta)^{3/2}} = \frac{1-r_0^2}{2}\int_{-1}^{1}\frac{t^2\; dt}{(1+r_0^2-2r_0 t)^{3/2}} ,$ substituting $t = \cos\theta$ .
Put in the numbers. With $r_0 = \tfrac12$ : $1+r_0^2 = \tfrac54$ , $2r_0 = 1$ , $(1-r_0^2)/2 = \tfrac38$ , so $u = \frac38\int_{-1}^{1}\frac{t^2\; dt}{\bigl(\tfrac54 - t\bigr)^{3/2}} .$
Evaluate. Put $s = \tfrac54 - t$ , so that $t = \tfrac54 - s$ and $dt = -ds$ , and expand the numerator: $\frac{t^2}{s^{3/2}} = \frac{\bigl(\tfrac54-s\bigr)^2}{s^{3/2}} = \frac{25}{16}s^{-3/2} - \frac{5}{2}s^{-1/2} + s^{1/2} .$ An antiderivative with respect to $s$ is $-\tfrac{25}{8}s^{-1/2} - 5s^{1/2} + \tfrac23 s^{3/2}$ , so an antiderivative with respect to $t$ is its negative, $\Phi(s) = \frac{25}{8}\,s^{-1/2} + 5\,s^{1/2} - \frac23\, s^{3/2} .$ As $t$ runs from $-1$ to $1$ , $s$ runs from $\tfrac94$ to $\tfrac14$ , so the integral equals $\Phi(\tfrac14) - \Phi(\tfrac94)$ : $\Phi\!\left(\tfrac14\right) = \frac{25}{4} + \frac52 - \frac{1}{12} = \frac{104}{12}, \qquad \Phi\!\left(\tfrac94\right) = \frac{25}{12} + \frac{15}{2} - \frac94 = \frac{88}{12} .$ Hence the integral is $\tfrac{104-88}{12} = \tfrac43$ , and $u = \tfrac38\cdot\tfrac43 = \tfrac12$ .
Independent check. Since $P_2(t) = \tfrac12(3t^2-1)$ , the datum is $f = \tfrac13 + \tfrac23 P_2(\cos\theta)$ , and $1$ and $r^2P_2(\cos\theta)$ are both harmonic (they are solid harmonics of degree $0$ and $2$ ). So the solution of this Dirichlet problem is $u = \tfrac13 + \tfrac23 r^2P_2(\cos\theta)$ , which at $r = \tfrac12$ , $\theta = 0$ gives $\tfrac13 + \tfrac23\cdot\tfrac14\cdot1 = \tfrac13+\tfrac16 = \tfrac12$ . The two computations agree exactly.
A second value. At the centre, (P) gives the average of $f$ over the sphere, $\tfrac{1}{4\pi}\oint\cos^2\theta\, dS = \tfrac12\int_{-1}^1 t^2\, dt = \tfrac13$ , matching $u(\mathbf{0}) = \tfrac13$ from the closed form.

Example(Image charge, image force, and induced charge)Original

Electrostatics is where the images are literally charges. Two standard numbers follow at once from the two Green's functions above. Throughout, $\phi$ is the potential, $\nabla^2\phi = -\rho_{\mathrm e}/\varepsilon_0$ , and a grounded conductor is the surface $\phi = 0$ .

Charge above a grounded plane. A charge $q$ at height $h$ above the grounded plane $z=0$ has potential, for $z>0$ , $\phi(\mathbf{r}) = \frac{q}{4\pi\varepsilon_0}\left(\frac{1}{\lvert\mathbf{r}-h\hat{\mathbf{z}}\rvert} - \frac{1}{\lvert\mathbf{r}+h\hat{\mathbf{z}}\rvert}\right),$ which is $-q/\varepsilon_0$ times the half-space Green's function. The force on the charge is that of the image alone, attractive, of magnitude $q^2/[4\pi\varepsilon_0(2h)^2] = q^2/(16\pi\varepsilon_0 h^2)$ . The induced surface charge is $\sigma = \varepsilon_0 E_z|_{z=0}$ with $E_z = -\partial\phi/\partial z$ ; differentiating and setting $z=0$ , $\sigma(\varrho) = -\frac{q h}{2\pi(\varrho^2+h^2)^{3/2}} , \qquad \int_0^\infty\sigma\,2\pi\varrho\; d\varrho = -qh\Bigl[-(\varrho^2+h^2)^{-1/2}\Bigr]_0^\infty = -q .$ The plane carries exactly the image charge, spread out; and $\sigma = -q\,\partial G_D/\partial n$ , so the Poisson kernel of (H) is the induced charge density, up to the factor $-q$ .
Charge outside a grounded sphere. A charge $q$ at distance $d>a$ from the centre of a grounded sphere of radius $a$ has image $q' = -qa/d$ at distance $a^2/d$ from the centre, on the same ray. The separation of the two is $d - a^2/d = (d^2-a^2)/d$ , so the force on $q$ is attractive with magnitude $F = \frac{1}{4\pi\varepsilon_0}\frac{q\lvert q'\rvert}{\bigl[(d^2-a^2)/d\bigr]^2} = \frac{q^2 a d}{4\pi\varepsilon_0\,(d^2-a^2)^2} .$ For $a = 1$ and $d = 2$ : $q' = -q/2$ (so the sphere carries half the charge's worth, with the opposite sign), the separation is $\tfrac32$ , and $F = 2q^2/(4\pi\varepsilon_0\cdot 9) = q^2/(18\pi\varepsilon_0)$ .
Limits as checks. As $d\to\infty$ , $F \approx q^2a/(4\pi\varepsilon_0 d^3)$ : the sphere looks like an induced dipole, falling off as $d^{-3}$ . As $d\to a^+$ , write $d = a+\epsilon$ ; then $d^2-a^2 \approx 2a\epsilon$ and $F\approx q^2 a\cdot a/(4\pi\varepsilon_0\cdot4a^2\epsilon^2) = q^2/(16\pi\varepsilon_0\epsilon^2)$ , which is exactly the plane result of part 1 with $h = \epsilon$ — as it must be, since close up a sphere is a plane.

ApplicationOriginal

Conductors and shields. Every undergraduate electrostatics problem with a grounded plane, a grounded sphere, two intersecting planes or a sphere-and-plane pair is an image problem, and the "image charge" of those problems is the corrector of a Green's function. The practical pay-off is the force and the induced charge density computed above: the attraction of an ion to a metal surface (and hence the image-potential contribution to the work function and to physisorption energies), the capacitance of a sphere above a plane, and the field enhancement at a grounded electrode.

Method of images outside electrostatics. The same construction solves any problem with the same operator and the same geometry. In steady heat conduction, an image sink across an isothermal wall enforces fixed wall temperature, and an image source of the same sign across an insulated wall enforces zero flux — the Neumann variant noted in the remark. In groundwater hydraulics the Theis solution for a pumping well near a constant-head river (image recharge well) or near an impermeable barrier (image pumping well) is the textbook application, and the whole industry of well-field design rests on it. In fluid mechanics the image vortex or image source enforces a no-through-flow boundary at a wall, and the image of a vortex in a cylinder uses the same inversion point $a^2\mathbf{r}_0/r_0^2$ .

Why the sphere result keeps reappearing. The map $\mathbf{r}\mapsto a^2\mathbf{r}/\lvert\mathbf{r}\rvert^2$ is inversion in the sphere, a conformal map in two dimensions and a map preserving the harmonicity of $\lvert\mathbf{x}\rvert^{2-n}u(a^2\mathbf{x}/\lvert\mathbf{x}\rvert^2)$ in $n$ . Evans builds the ball's corrector in exactly that form. That is why the Kelvin inversion appears again in the Poisson integral formula and in the theory of harmonic functions on balls generally, and why the image system of a sphere is a single point rather than a continuum.

Limits of the method, and what replaces it. The list of regions with finite image systems is short: half-spaces, slabs, wedges of angle $\pi/m$ , balls and the regions assembled from them by reflection. For anything else the corrector must be computed, and the Green's function is then obtained from an eigenfunction expansion, from a conformal map in two dimensions, or numerically; the discretised representation formula with unknowns on the surface alone is the boundary-element method, and the image Green's functions above are what such codes use for layered geometries.

Method

Green's functions for the Helmholtz equation

Adapted

Add to PDF

Computes the free-space outgoing Green's function of $\nabla^2+k^2$ by Fourier transforming and evaluating the remaining integral by residues with the poles displaced off the real axis, giving $-\tfrac{i}{2k}e^{ik|x|}$ in one dimension, $-\tfrac{i}{4}H_0^{(1)}(k\varrho)$ in two and $-e^{ikr}/(4\pi r)$ in three, and reads off the modified Helmholtz case by $k\mapsto i\kappa$ .

Prerequisites: The Helmholtz equation, Infinite real integrals and Jordan's lemma, Green's functions for Laplace and Poisson problems, Fourier transform in several dimensions

Note(Which Green's function, and the sign convention)Original

The operator $\nabla^2 + k^2$ with $k > 0$ real has non-trivial solutions of the homogeneous equation that decay in no direction — plane waves — so a Green's function is not fixed by asking it to vanish at infinity. One must say which waves are allowed, and the physical choice is outgoing: the Green's function should behave like a wave travelling away from its source, so that when multiplied by the time factor $e^{-i\omega t}$ with $\omega = ck$ it describes radiation leaving the source rather than converging on it. That is the convention of Arfken, Weber and Harris and of Hassani, and the one adopted here. The opposite choice, $e^{-ikr}$ in place of $e^{ikr}$ , gives the incoming Green's function, which is the complex conjugate; their average is the standing-wave Green's function.

As in Green's functions for Laplace and Poisson problems we write the defining equation as

\bigl(\nabla^2 + k^2\bigr) G(\mathbf{r}) = \delta(\mathbf{r}) ,

which is the convention of Riley, Hobson and Bence, of Arfken and of Hassani, and which makes $G$ negative-definite in sign at short range in three dimensions. Because the free-space problem is translation invariant, $G(\mathbf{r},\mathbf{r}') = G(\mathbf{r}-\mathbf{r}')$ and we may put the source at the origin throughout; and because it is also rotation invariant, $G$ depends only on $\lvert\mathbf{r}\rvert$ .

Definition(Free-space Green's function and the outgoing condition)AdaptedAdapted from Arfken, Weber & Harris (2013), §10.2, Table 10.1 and its footnotes, p. 463; Hassani (2013), §22.4.2, p. 683

Let $n \in\{1,2,3\}$ and $k>0$ . A distribution $G$ on $\mathbb{R}^n$ is the outgoing (or retarded) free-space Green's function of the Helmholtz operator if

$\bigl(\nabla^2+k^2\bigr)G = \delta$ on $\mathbb{R}^n$ in the sense of distributions, so that $G$ is smooth and solves $(\nabla^2+k^2)G=0$ away from the origin and carries at the origin the same singularity as the fundamental solution of $\nabla^2$ ; and
$G$ satisfies the outgoing radiation condition: as $r = \lvert\mathbf{r}\rvert\to\infty$ its asymptotic dependence on $r$ is through $e^{+ikr}$ and not $e^{-ikr}$ (equivalently, $\partial G/\partial r - ikG = o(r^{-(n-1)/2})$ , the Sommerfeld condition of The Helmholtz equation). In one dimension, where this reads $o(1)$ , the symbols are to be understood in the same radial sense: $r = \lvert x\rvert$ and $\partial/\partial r$ is the outward derivative, $+\partial_x$ on $x>0$ and $-\partial_x$ on $x<0$ . With that reading the condition does its work even at $n=1$ : for $G = -\tfrac{i}{2k}e^{ik\lvert x\rvert}$ the combination $\partial G/\partial r - ikG$ vanishes identically, while for the incoming function $-\tfrac{i}{2k}e^{-ik\lvert x\rvert}$ it equals $-e^{-ik\lvert x\rvert}$ , of modulus $1$ , which is not $o(1)$ .

For the modified Helmholtz operator $\nabla^2-\kappa^2$ with $\kappa>0$ , condition 2 is replaced by the requirement that $G\to0$ in every direction, which is now available because the homogeneous solutions grow or decay exponentially.

Condition 1 has the useful corollary, which Arfken, Weber and Harris set as an exercise, that the term $\pm k^2 G$ contributes nothing to the delta function: its integral over a small ball about the origin tends to $0$ , so the local structure of $G$ at its pole is exactly that of the Laplace Green's function, namely $\tfrac12\lvert x\rvert$ , $(2\pi)^{-1}\ln\varrho$ and $-1/(4\pi r)$ in one, two and three dimensions. This is what fixes the normalisation in every derivation below.

Method(Fourier transform and residues)AdaptedAdapted from Hassani (2013), §22.4.1-22.4.2, pp. 682-683 (the Fourier representation and its angular reduction); Arfken, Weber & Harris (2013), §10.1, Example 10.1.3, pp. 454-455 (the identification of the outgoing solution)

Transform. Write $G(\mathbf{r}) = \displaystyle\int\frac{d^n q}{(2\pi)^n}\,\tilde G(\mathbf{q})\,e^{i\mathbf{q}\cdot\mathbf{r}}$ and use $\delta(\mathbf{r}) = \int\frac{d^nq}{(2\pi)^n}e^{i\mathbf{q}\cdot\mathbf{r}}$ (Fourier transform in several dimensions). Since $\nabla^2 \mapsto -q^2$ , $\bigl(k^2 - q^2\bigr)\tilde G(\mathbf{q}) = 1, \qquad \tilde G(\mathbf{q}) = \frac{1}{k^2-q^2} .$
Notice the obstruction. $\tilde G$ is singular on the whole sphere $q = k$ , so the inverse transform is not an ordinary integral. This is the analytic face of the fact that the Green's function is not unique: the homogeneous solutions live exactly on that sphere.
Displace the poles. Impose the outgoing condition by giving $k$ a small positive imaginary part, $k \to k+i\delta$ with $\delta\downarrow0$ (equivalently $k^2\to k^2+i\epsilon$ ). The poles then sit at $q = \pm(k+i\delta)$ , one in each half-plane, and the integral is well defined.
Reduce to one radial integral. Do the angular integration. In three dimensions $\int e^{i\mathbf{q}\cdot\mathbf{r}}\, d\Omega_{\mathbf q} = 4\pi\dfrac{\sin qr}{qr}$ , in two it gives $2\pi J_0(q\varrho)$ , and in one there is nothing to do.
Close the contour. Extend the radial integral to the whole real line using the parity of the integrand, write the trigonometric factor in exponentials, and apply Jordan's lemma and the residue theorem, closing in the half-plane in which the exponential decays. Only one pole is enclosed, and it is the one that produces $e^{+ikr}$ .
Let $\delta\to0$ and check the normalisation against the short-range behaviour of step 1 of the definition.
For $\nabla^2-\kappa^2$ , substitute $k = i\kappa$ in the answer. The outgoing condition becomes decay, and no pole displacement is needed, since $\tilde G = -1/(\kappa^2+q^2)$ has no real poles at all.

Method(The nine free-space Green's functions)AdaptedAdapted from Arfken, Weber & Harris (2013), §10.2, Table 10.1, p. 463, and §10.1, Eq. (10.28), p. 454 and Exercise 10.1.9, p. 457; Hassani (2013), §22.4.2, p. 683

With the conventions above — outgoing waves for $\nabla^2+k^2$ , decay at infinity for $\nabla^2-\kappa^2$ and for $\nabla^2$ in three dimensions, and an arbitrary additive constant for $\nabla^2$ in one and two dimensions — the free-space Green's functions are

	$\nabla^2$	$\nabla^2+k^2$	$\nabla^2-\kappa^2$
one dimension	$\tfrac12\lvert x\rvert$	$-\dfrac{i}{2k}e^{ik\lvert x\rvert}$	$-\dfrac{1}{2\kappa}e^{-\kappa\lvert x\rvert}$
two dimensions	$\dfrac{1}{2\pi}\ln\varrho$	$-\dfrac{i}{4}H_0^{(1)}(k\varrho)$	$-\dfrac{1}{2\pi}K_0(\kappa\varrho)$
three dimensions	$-\dfrac{1}{4\pi r}$	$-\dfrac{e^{ikr}}{4\pi r}$	$-\dfrac{e^{-\kappa r}}{4\pi r}$

Here $H_0^{(1)} = J_0+iY_0$ is the Hankel function of the first kind and $K_0$ the modified Bessel function. This is Table 10.1 of Arfken, Weber and Harris. Each column is obtained from the next by a substitution: setting $k=i\kappa$ in the middle column gives the right-hand one, and letting $k\to0$ in the middle column recovers the left-hand one up to the additive constants that the one- and two-dimensional Laplace Green's functions carry anyway. In three dimensions Hassani writes the middle entry as $-e^{\pm i\mu r}/(4\pi r)$ and identifies the two signs with incoming and outgoing waves; with the time convention $e^{-i\omega t}$ it is the upper sign that is outgoing.

Once $G$ is known, the solution of $(\nabla^2+k^2)\psi = f$ on all of $\mathbb{R}^n$ with outgoing behaviour is

\psi(\mathbf{r}) = \int G(\mathbf{r}-\mathbf{r}')\, f(\mathbf{r}')\; d^n r' ,

for $f$ of compact support, by the superposition principle of Riley, Hobson and Bence; any solution of the homogeneous equation may be added, and is added in scattering problems to supply the incident wave.

Idea of the proofOriginal

In free space the operator is diagonal in Fourier variables, so inverting it is a division — except that $k^2-q^2$ vanishes on a sphere, which is precisely the statement that the homogeneous equation has solutions and the Green's function is therefore not unique. Everything interesting is in how one steps around that sphere of poles, and the rule is: the choice of contour is the choice of boundary condition.

Pushing the pole at $q=+k$ up into the upper half-plane and the one at $q=-k$ down means that, when the remaining radial integral is closed in the upper half-plane (which is where $e^{iqr}$ decays, for $r>0$ ), only the pole at $+k$ is caught, and it contributes $e^{+ikr}$ — an outgoing wave. Pushing them the other way catches $-k$ and gives $e^{-ikr}$ , the incoming solution. Taking the principal value catches half of each and gives the standing wave, in three dimensions $-\cos(kr)/(4\pi r)$ . Being the half-sum of two Green's functions, that is itself a Green's function: it carries the same $-1/(4\pi r)$ singularity and solves the homogeneous equation away from the origin. So nothing is wrong with it as a solution of $(\nabla^2+k^2)G=\delta$ ; what it is not is outgoing. It is as much incoming as outgoing, so it fails condition 2 of the definition and describes a cavity, where waves return, rather than radiation into free space. The contour is thus not chosen to make $G$ a Green's function at all, but to say which of the many Green's functions one wants.

In two dimensions the contour integral runs into Bessel functions, and it is quicker to work in position space: away from the source $G$ obeys Bessel's equation of order zero, so it is a combination of the two Hankel functions; the outgoing condition picks $H_0^{(1)}$ by its asymptotics, and the strength is fixed by matching the logarithmic singularity at the origin to the Laplace Green's function, which is legitimate because the $k^2G$ term is too weak to affect the singularity.

Proof(One dimension)AdaptedAdapted from Hassani (2013), §22.4.2, the Fourier representation and the pole displacement, pp. 682-683; Arfken, Weber & Harris (2013), §10.1, Example 10.1.3, Eq. (10.28), pp. 454-455

Step 1 (transform). With $n=1$ , step 1 of the method gives

G(x) = \int_{-\infty}^{\infty}\frac{dq}{2\pi}\;\frac{e^{iqx}}{k^2-q^2} ,

and after the displacement $k\to k+i\delta$ the integrand acquires simple poles at $q = k+i\delta$ , above the real axis, and at $q = -k-i\delta$ , below it.

Step 2 (close the contour). The factor $1/(k^2-q^2)$ tends to $0$ as $\lvert q\rvert\to\infty$ , so for $x>0$ Jordan's lemma applies with the factor $e^{iqx}$ and the semicircle in the upper half-plane contributes nothing; the contour may be closed there. For $x<0$ the exponential decays in the lower half-plane and the contour closes there instead, with the opposite orientation.

Step 3 (residues). Writing $k^2-q^2 = -(q-q_0)(q+q_0)$ with $q_0 = k+i\delta$ , the residue at $q=q_0$ is

\lim_{q\to q_0}(q-q_0)\,\frac{e^{iqx}}{-(q-q_0)(q+q_0)} = -\frac{e^{iq_0x}}{2q_0} ,

and at $q=-q_0$ it is $+e^{-iq_0x}/(2q_0)$ by the same computation. Hence for $x>0$ ,

G(x) = \frac{1}{2\pi}\cdot 2\pi i\left(-\frac{e^{iq_0x}}{2q_0}\right) \xrightarrow[\delta\to0]{} -\frac{i}{2k}\,e^{ikx} ,

while for $x<0$ the clockwise contour supplies $-2\pi i$ and

G(x) = \frac{1}{2\pi}(-2\pi i)\,\frac{e^{-iq_0x}}{2q_0}\xrightarrow[\delta\to0]{} -\frac{i}{2k}\,e^{-ikx} .

The two cases combine into $G(x) = -\dfrac{i}{2k}e^{ik\lvert x\rvert}$ , which is outgoing on both sides, as the displacement was chosen to arrange.

Step 4 (check). $G$ is continuous at $0$ with $G(0) = -i/2k$ ; differentiating the two branches gives $G'(x) = \tfrac12 e^{ikx}$ for $x>0$ and $G'(x) = -\tfrac12 e^{-ikx}$ for $x<0$ , so $G'(0^+)-G'(0^-) = \tfrac12+\tfrac12 = 1$ , which is the jump that $G''=\delta$ demands; away from $0$ , $G'' = -k^2G$ . Arfken, Weber and Harris reach the same function by the one-dimensional construction from the two homogeneous solutions $e^{\mp ikx}$ with the Wronskian factor $A = -i/2k$ , a route that does not use a contour at all. $\square$

Proof(Three dimensions)AdaptedAdapted from Hassani (2013), §22.4.2, pp. 682-683, the angular reduction and the three-dimensional result; Arfken, Weber & Harris (2013), §10.2, Table 10.1 and Exercise 10.2.4, pp. 463, 466

Step 1 (angular integration). With the source at the origin and $r = \lvert\mathbf{r}\rvert$ , taking the polar axis of $\mathbf{q}$ along $\mathbf{r}$ gives $\int e^{i\mathbf{q}\cdot\mathbf{r}}d\Omega_{\mathbf{q}} = 2\pi\int_{-1}^{1}e^{iqr\mu}d\mu = 4\pi\,\dfrac{\sin qr}{qr}$ . Hence

G(r) = \frac{1}{(2\pi)^3}\int_0^\infty dq\; q^2\,\frac{1}{k^2-q^2}\cdot\frac{4\pi\sin qr}{qr} = \frac{1}{2\pi^2 r}\int_0^\infty \frac{q\,\sin qr}{k^2-q^2}\; dq ,

using $4\pi/(2\pi)^3 = 1/(2\pi^2)$ .

Step 2 (extend to the whole line). The integrand $q\sin(qr)/(k^2-q^2)$ is even in $q$ , so

G(r) = \frac{1}{4\pi^2 r}\int_{-\infty}^{\infty}\frac{q\,\sin qr}{k^2-q^2}\; dq .

Step 3 (to a single exponential). Put $\sin qr = (e^{iqr}-e^{-iqr})/2i$ and let $J = \int_{-\infty}^\infty q\,e^{iqr}(k^2-q^2)^{-1}dq$ . Substituting $q\to-q$ in the second piece turns $\int q e^{-iqr}(k^2-q^2)^{-1}dq$ into $-J$ , so the whole integral is $\tfrac{1}{2i}\bigl(J-(-J)\bigr) = -iJ$ .

Step 4 (residues). After the displacement $k\to k+i\delta$ , the integrand of $J$ has poles at $\pm(k+i\delta)$ , and $q/(k^2-q^2) \to 0$ at large $\lvert q\rvert$ , so Jordan's lemma permits closing in the upper half-plane for $r>0$ . Only $q_0 = k+i\delta$ is enclosed, with residue

\lim_{q\to q_0}(q-q_0)\frac{q\,e^{iqr}}{-(q-q_0)(q+q_0)} = -\frac{q_0 e^{iq_0 r}}{2q_0} = -\frac{e^{iq_0r}}{2} ,

so $J = 2\pi i\cdot(-\tfrac12 e^{iq_0r}) \to -i\pi e^{ikr}$ as $\delta\to0$ .

Step 5 (assemble). The radial integral is $-iJ = -i(-i\pi e^{ikr}) = -\pi e^{ikr}$ , whence

G(r) = \frac{-\pi e^{ikr}}{4\pi^2 r} = -\frac{e^{ikr}}{4\pi r} .

Step 6 (checks). As $k\to0$ this is the Laplace Green's function $-1/(4\pi r)$ , and as $r\to0$ it behaves like $-1/(4\pi r)$ for any $k$ , confirming the normalisation demanded by condition 1 of the definition. Away from the origin, $\nabla^2(e^{ikr}/r) = r^{-1}(d^2/dr^2)(e^{ikr}) = -k^2 e^{ikr}/r$ , so $(\nabla^2+k^2)G=0$ there. Finally $\partial G/\partial r - ikG = e^{ikr}/(4\pi r^2) = O(r^{-2})$ , so the Sommerfeld condition holds. Had we displaced the poles the other way, the enclosed pole would have been $-k$ and the answer $-e^{-ikr}/(4\pi r)$ : the incoming Green's function, which is the complex conjugate, as Hassani's two signs record. $\square$

Proof(Two dimensions)AdaptedAdapted from Arfken, Weber & Harris (2013), §10.2, Table 10.1 and footnote b, p. 463, and Exercise 10.2.3, p. 466; Hassani (2013), §22.4.2, the general-$m$ Hankel form, p. 682

Here the contour route runs into Bessel integrals — the method of step 4 leaves $G = (2\pi)^{-1}\int_0^\infty q J_0(q\varrho)(k^2-q^2)^{-1}dq$ , and Hassani's general-dimension computation of that integral produces a Hankel function $H^{(1)}_{m/2-1}$ , whose $m=2$ case is $H_0^{(1)}$ . We instead fix the answer in position space, which needs only the radial equation and the two facts about $H_0^{(1)}$ that are proved in the Bessel chapter.

Step 1 (the radial equation). Away from the source $G$ depends only on $\varrho$ and satisfies

\frac{1}{\varrho}\frac{d}{d\varrho}\left(\varrho\,\frac{dG}{d\varrho}\right) + k^2 G = 0 ,

which is Bessel's equation of order zero in the variable $k\varrho$ . Its solutions are spanned by $H_0^{(1)}(k\varrho)$ and $H_0^{(2)}(k\varrho)$ , so $G = A\,H_0^{(1)}(k\varrho) + B\,H_0^{(2)}(k\varrho)$ for $\varrho>0$ .

Step 2 (the outgoing condition fixes $B=0$ ). By the large-argument asymptotics, $H_0^{(1)}(z)\sim\sqrt{2/\pi z}\;e^{i(z-\pi/4)}$ and $H_0^{(2)}(z)\sim\sqrt{2/\pi z}\;e^{-i(z-\pi/4)}$ . The first depends on $\varrho$ through $e^{+ik\varrho}$ and the second through $e^{-ik\varrho}$ , so condition 2 of the definition forces $B=0$ .

Step 3 (the singularity fixes $A$ ). By the corollary in the definition, $G$ must have the same singularity at $\varrho=0$ as the two-dimensional Laplace Green's function, that is $G = (2\pi)^{-1}\ln\varrho + O(1)$ . The small-argument behaviour of the Hankel function is $H_0^{(1)}(z) = 1 + \tfrac{2i}{\pi}\bigl(\ln(z/2)+\gamma\bigr) + O(z^2\ln z)$ , since $J_0(z)\to1$ and $Y_0(z) = \tfrac{2}{\pi}(\ln(z/2)+\gamma) + O(z^2)$ . So

A\,H_0^{(1)}(k\varrho) = \frac{2iA}{\pi}\ln\varrho + \text{(constant)} + o(1) ,

and matching the coefficient of $\ln\varrho$ gives $2iA/\pi = 1/(2\pi)$ , that is

A = \frac{1}{4i} = -\frac{i}{4}, \qquad G(\varrho) = -\frac{i}{4}\,H_0^{(1)}(k\varrho) ,

which is the entry in Arfken's Table 10.1. $\square$

Remark(Where the choices enter, and what else they could have been)Original

The pole displacement is the boundary condition, nothing more. $\tilde G = (k^2-q^2)^{-1}$ is the same in all cases; outgoing, incoming and standing-wave Green's functions differ only in how the contour passes the two poles. In three dimensions they are $-e^{ikr}/(4\pi r)$ , its conjugate, and the average $-\cos(kr)/(4\pi r)$ ; the last is what Arfken, Weber and Harris ask for in the exercise on standing waves, and it is the one relevant to a cavity rather than to radiation into free space.
A principal value is not a radiation Green's function — and for the wave operator it is not a Green's function at all. For $\nabla^2+k^2$ the symmetric contour gives the standing wave $-\cos(kr)/(4\pi r)$ of the previous bullet, which does solve $(\nabla^2+k^2)G=\delta$ ; it is only the radiation condition that it fails, so here one must choose a side, and the physics chooses it. The analogous prescription for the wave operator is worse, and it is there that Hassani places his warning: the principal value of his $\omega$ integral produces a function with no delta source whatever, so it solves the homogeneous wave equation and is no Green's function at all (his Problem 22.25); the useful ones, he concludes, come from displacing the poles instead. The difference is that for the wave operator the contour must be closed in opposite half-planes according to the sign of $t$ , so the symmetric prescription combines the retarded and advanced functions in a way that cancels their delta functions; in the spatial integral above the contour closes in the upper half-plane for every $r>0$ , and no such cancellation can occur.
Why $\pm k^2$ cannot change the singularity. Near the pole, $G$ blows up like the Laplace Green's function, which in $n$ dimensions is integrable; multiplying by $k^2$ and integrating over a ball of radius $\varepsilon$ gives $O(\varepsilon^2)$ in three dimensions, $O(\varepsilon^2\ln\varepsilon)$ in two and $O(\varepsilon^2)$ in one. So the $k^2G$ term contributes nothing to the delta function, which is why the normalisation can always be read off from the Laplace case. The qualitative statement is Exercise 10.2.3 of Arfken, Weber and Harris, which asks one to show that neither $+k^2$ nor $-k^2$ alters how $G$ behaves in "the immediate vicinity of the singular point"; the $\varepsilon$ -ball estimates just given are our own, and are not the integral that the exercise goes on to display.
Dimension matters more than one expects. The three answers in the Helmholtz column are of three different kinds: a bounded oscillation of constant amplitude in one dimension, a logarithmically singular function with $\varrho^{-1/2}$ decay in two, and a $1/r$ decay in three. This is the same difference that Arfken, Weber and Harris stress for the Laplace column: "flatland" physics is genuinely different, even for problems with translational symmetry in one direction. Indeed the two-dimensional Green's function is the three-dimensional one integrated along a line, which is how $H_0^{(1)}$ arises from $e^{ikr}/r$ .
Why $-i/4$ and not $-1/4$ . The factor of $i$ is an artefact of writing the outgoing solution as a Hankel function rather than as an exponential: $H_0^{(1)}$ is complex at small argument, and its real part $J_0$ is regular there while the singular part sits in $iY_0$ . Matching a real logarithm therefore demands an imaginary coefficient.
Bounded regions. All of the above is free space. On a bounded region one adds a solution of the homogeneous Helmholtz equation to meet the boundary condition, exactly as for Laplace's equation; but unlike the Laplace case this fails when $k^2$ is an eigenvalue of $-\nabla^2$ on the region, for then the corrector problem is unsolvable and no Green's function exists. Arfken, Weber and Harris note the same phenomenon in the form that the Dirichlet problem for Helmholtz's equation on a closed boundary has solutions "only for certain values of $k$ ", and the resolution is the Fredholm alternative; the eigenfunction expansion $G = \sum_n \varphi_n(\mathbf{r})\varphi_n^*(\mathbf{r}')/(\lambda_n-\lambda)$ of Eigenfunction expansion of Green's functions shows the breakdown as a vanishing denominator.
Modified Helmholtz needs no prescription. For $\nabla^2-\kappa^2$ the transform is $-(\kappa^2+q^2)^{-1}$ , which has no real poles; the contour is unambiguous and the answer decays. Substituting $k=i\kappa$ in the middle column reproduces the right-hand one: $-\tfrac{i}{2k}e^{ik\lvert x\rvert}\mapsto -\tfrac{1}{2\kappa}e^{-\kappa\lvert x\rvert}$ , $-e^{ikr}/(4\pi r)\mapsto -e^{-\kappa r}/(4\pi r)$ , and, using $K_0(z) = \tfrac{i\pi}{2}H_0^{(1)}(iz)$ , $-\tfrac i4 H_0^{(1)}(i\kappa\varrho)\mapsto -\tfrac{1}{2\pi}K_0(\kappa\varrho)$ . Hassani obtains the three-dimensional case directly as the Yukawa potential of a unit charge.

Example(A uniform source in a slab, in one dimension)Original

Solve $\psi'' + k^2\psi = f$ on $\mathbb{R}$ with outgoing behaviour, where $f(x) = 1$ for $\lvert x\rvert<L$ and $f=0$ outside.

Superpose. $\displaystyle\psi(x) = -\frac{i}{2k}\int_{-L}^{L}e^{ik\lvert x-y\rvert}\, dy .$
Outside, $x>L$ . Then $\lvert x-y\rvert = x-y$ throughout, so $\psi(x) = -\frac{i}{2k}e^{ikx}\int_{-L}^{L}e^{-iky}\, dy = -\frac{i}{2k}e^{ikx}\cdot\frac{2\sin kL}{k} = -\frac{i\sin kL}{k^2}\,e^{ikx} ,$ and by symmetry $\psi(x) = -i\sin(kL)k^{-2}e^{-ikx}$ for $x<-L$ ; both are outgoing, with the amplitude $\lvert\sin kL\rvert/k^2$ .
Inside, $\lvert x\rvert<L$ . Split the integral at $y=x$ : $\int_{-L}^{x}e^{ik(x-y)}dy = \frac{e^{ik(x+L)}-1}{ik}, \qquad \int_{x}^{L}e^{ik(y-x)}dy = \frac{e^{ik(L-x)}-1}{ik} ,$ whose sum is $\bigl(2e^{ikL}\cos kx - 2\bigr)/ik$ . Multiplying by $-i/2k$ , $\psi(x) = \frac{1 - e^{ikL}\cos kx}{k^2} .$
Check the equation. Inside, $\psi'' = e^{ikL}\cos kx$ and $k^2\psi = 1-e^{ikL}\cos kx$ , so $\psi''+k^2\psi = 1 = f$ . Outside, $\psi''+k^2\psi=0$ . Both hold, so $\psi$ solves the problem.
Check the matching. At $x=L$ the inside formula gives $(1-e^{ikL}\cos kL)/k^2$ and the outside one $-i\sin(kL)e^{ikL}/k^2$ ; these agree because $e^{ikL}\cos kL - ie^{ikL}\sin kL = e^{ikL}e^{-ikL} = 1$ . The first derivative matches too, as it must: $f$ is bounded, so $\psi''$ is bounded and $\psi'$ continuous.
Numbers, and the silent slab. Take $k=1$ and $L = \pi/2$ . Then $\psi(0) = 1 - e^{i\pi/2} = 1-i$ , of modulus $\sqrt2 = 1.414214$ , and for $x>\pi/2$ , $\psi(x) = -ie^{ix}$ , of modulus exactly $1$ — the largest value the outgoing amplitude $\lvert\sin kL\rvert/k^2$ can take at $k=1$ . For $k=1$ and $L=\pi$ , by contrast, $\sin kL$ vanishes and the slab radiates nothing at all: $\psi = (1-e^{i\pi}\cos x)/1 = 1+\cos x$ inside, identically $0$ outside. The silent case is $kL = n\pi$ , that is $L = n\lambda/2$ with $\lambda = 2\pi/k$ : the half-width $L$ must be a whole number of half-wavelengths, so that the slab's full width $2L$ is a whole number of full wavelengths. (The $L=\pi/2$ case above has width $2L = \pi = \lambda/2$ , one half-wavelength, and radiates the maximum, which shows that the condition is on $L$ and not on the width.) The two edges then interfere destructively.

Example(A spherical shell source in three dimensions)Original

Solve $(\nabla^2+k^2)\psi = \sigma\,\delta(r-b)$ on $\mathbb{R}^3$ with outgoing behaviour: a uniform source of strength $\sigma$ per unit area spread over the sphere of radius $b$ about the origin. (Total source strength $4\pi b^2\sigma$ .)

The one integral needed. By symmetry $\psi$ depends only on $r$ , and $\psi(r) = -\frac{\sigma}{4\pi}\oint_{\lvert\mathbf{r}'\rvert=b}\frac{e^{ik s}}{s}\; dS', \qquad s = \lvert\mathbf{r}-\mathbf{r}'\rvert .$ Write $dS' = 2\pi b^2\sin\theta\, d\theta$ with $\theta$ the angle at the origin, so $s^2 = r^2+b^2-2rb\cos\theta$ and $s\, ds = rb\sin\theta\, d\theta$ . The measure converts exactly: $\oint\frac{e^{iks}}{s}\, dS' = 2\pi b^2\int\frac{e^{iks}}{s}\cdot\frac{s\, ds}{rb} = \frac{2\pi b}{r}\int_{\lvert r-b\rvert}^{r+b}e^{iks}\, ds = \frac{2\pi b}{ikr}\Bigl(e^{ik(r+b)} - e^{ik\lvert r-b\rvert}\Bigr) .$ This is the whole computation; no Bessel functions are needed.
Outside, $r>b$ . Then $\lvert r-b\rvert = r-b$ and the bracket is $e^{ikr}(e^{ikb}-e^{-ikb}) = 2i\sin(kb)\,e^{ikr}$ , so $\psi(r) = -\frac{\sigma}{4\pi}\cdot\frac{4\pi b\sin kb}{kr}e^{ikr} = -\frac{\sigma b\sin kb}{k}\,\frac{e^{ikr}}{r} .$
Inside, $r<b$ . Then $\lvert r-b\rvert = b-r$ and the bracket is $e^{ikb}(e^{ikr}-e^{-ikr}) = 2i\sin(kr)e^{ikb}$ , so $\psi(r) = -\frac{\sigma b\, e^{ikb}}{k}\,\frac{\sin kr}{r} ,$ which is finite at the origin, with $\psi(0) = -\sigma b\,e^{ikb}$ .
Check the equation away from the shell. For any $g$ , $\nabla^2\bigl(g(r)/r\bigr) = r^{-1}g''(r)$ . With $g = e^{ikr}$ outside and $g=\sin kr$ inside, $g'' = -k^2g$ in both cases, so $(\nabla^2+k^2)\psi=0$ on $r\neq b$ . The outside form is outgoing and the inside form regular, as required.
Check the jump. Integrating $(\nabla^2+k^2)\psi = \sigma\delta(r-b)$ over a thin shell and using the divergence theorem gives $4\pi b^2[\psi'(b^+)-\psi'(b^-)] = 4\pi b^2\sigma$ , so the radial derivative must jump by $\sigma$ . From steps 2 and 3, $\psi'(b^+) = -\frac{\sigma\sin kb\;e^{ikb}(ikb-1)}{kb}, \qquad \psi'(b^-) = -\frac{\sigma e^{ikb}(kb\cos kb - \sin kb)}{kb} ,$ and the difference is $-\dfrac{\sigma e^{ikb}}{kb}\bigl[ikb\sin kb - kb\cos kb\bigr] = \sigma e^{ikb}\bigl(\cos kb - i\sin kb\bigr) = \sigma$ , using $e^{ikb}e^{-ikb}=1$ . The two formulas also agree at $r=b$ , both giving $-\sigma b\sin(kb)e^{ikb}/(kb)$ , so $\psi$ is continuous there.
Numbers. Take $\sigma=1$ and $k=1$ . For $b = 0.7$ : $\psi(0) = -0.7\,e^{0.7i} = -0.5353895 - 0.4509524\,i$ , and the far field has amplitude $b\sin(kb)/k = 0.7\sin 0.7 = 0.4509524$ . For $b = \pi$ , $\sin kb = 0$ and the shell is silent: $\psi\equiv0$ outside, while inside $\psi = -\pi e^{i\pi}\sin(r)/r = \pi\sin(r)/r$ , which rises to $\psi(0) = \pi$ . The silent radii are $kb = n\pi$ , that is $b = n\lambda/2$ : a breathing sphere whose radius is a whole number of half-wavelengths radiates nothing, and $b=\pi$ at $k=1$ is the first of them ( $\lambda = 2\pi$ , $b = \lambda/2$ ), with $b = \lambda, \tfrac32\lambda,\dots$ equally silent. This is the direct analogue of the one-dimensional slab of the previous example, and the reason a monopole radiator has a comb of blind sizes.

ApplicationOriginal

Radiation and antennas. The three-dimensional outgoing Green's function is the radiation field of a point monopole: a time-harmonic source $f(\mathbf{r})e^{-i\omega t}$ produces $\psi(\mathbf{r})e^{-i\omega t}$ with $\psi = \int G f$ , and at large $r$ , $e^{ikr}/r$ times an angular factor $\int f(\mathbf{r}')e^{-i\mathbf{k}\cdot\mathbf{r}'}d^3r'$ — the Fourier transform of the source evaluated on the sphere $\lvert\mathbf{k}\rvert = k$ . That one observation is the whole far-field theory of acoustic and electromagnetic radiators: the radiation pattern is the source's spatial spectrum restricted to the shell of radiating wavenumbers, the rest of the spectrum being evanescent. The silent shell of the second example is the extreme case of destructive interference between parts of a source, and the same arithmetic sets the directivity of a phased array.

Scattering. Writing the Schrödinger equation as $(\nabla^2+k^2)\psi = (2m/\hbar^2)V\psi$ and applying the outgoing Green's function turns it into the integral equation

\psi(\mathbf{r}) = e^{i\mathbf{k}_0\cdot\mathbf{r}} - \frac{2m}{\hbar^2}\int\frac{e^{ik\lvert\mathbf{r}-\mathbf{r}'\rvert}}{4\pi\lvert\mathbf{r}-\mathbf{r}'\rvert}V(\mathbf{r}')\psi(\mathbf{r}')\; d^3r' ,

the plane wave being the added homogeneous solution that supplies the incident beam. This is the Lippmann-Schwinger equation; iterating it once gives the Born approximation (Scattering and the Born approximation), and iterating it indefinitely gives the Neumann series. The outgoing condition is what makes the scattered part a wave leaving the target, and hence what makes the cross-section come out positive.

Screening, and why the modified operator appears. $\nabla^2-\kappa^2$ with the Green's function $-e^{-\kappa r}/(4\pi r)$ is the screened Coulomb or Yukawa kernel: $\kappa^{-1}$ is the Debye length in a plasma or electrolyte, the Thomas-Fermi screening length in a metal, and $\hbar/m_\pi c$ in Yukawa's original nuclear force. In two dimensions the same Green's function is $-K_0(\kappa\varrho)/2\pi$ , which is the response of a membrane on an elastic foundation and the kernel of linearised groundwater flow with leakage. The replacement $k\to i\kappa$ that connects the two columns of the table is, physically, the passage from a propagating wave to an evanescent one, and in a waveguide it happens exactly at cut-off.

Seismology, acoustics and optics in two dimensions. $-\tfrac i4 H_0^{(1)}(k\varrho)$ is the two-dimensional point-source field used in full-waveform inversion, in the boundary-element modelling of noise barriers, and in the analysis of surface waves; its $\varrho^{-1/2}$ decay and its $-\pi/4$ phase lag relative to a three-dimensional source are the signatures of a line source, and getting that phase right is what distinguishes a two-dimensional from a three-dimensional forward model. The same function, continued to imaginary argument, is the $K_0$ of the screened problem.

Time-dependent problems. Because a time-harmonic factor $e^{-i\omega t}$ converts the wave equation into Helmholtz's equation with $k=\omega/c$ , the outgoing Green's function here is the Fourier transform in time of the retarded Green's function of the wave equation; that correspondence, and the delta-function on the light cone that it implies, is the subject of Green's functions for the wave and heat equations.

Method

Green's functions for the wave and heat equations

Adapted

Add to PDF

Shows that the causal Green's function of the heat operator is the heat kernel and that of the three-dimensional wave operator is $\delta(t-r/c)/(4\pi c^2 r)$ , a delta function on the light cone, derives the retarded potential from it, and identifies Duhamel's principle as convolution in time with the causal Green's function.

Prerequisites: Green's functions for the Helmholtz equation, The heat kernel, Delta-function identities, Gaussian integrals

Note(Conventions: causal, retarded, advanced)Original

We write the two evolution operators as

L_{\mathrm{heat}} = \frac{\partial}{\partial t} - \kappa\nabla^2 , \qquad L_{\mathrm{wave}} = \frac{\partial^2}{\partial t^2} - c^2\nabla^2 ,

with $\kappa>0$ the diffusivity and $c>0$ the wave speed, and keep the convention of the previous units, $LG = \delta$ :

L\,G(\mathbf{r},t) = \delta^{(n)}(\mathbf{r})\,\delta(t) .

Both operators have translation-invariant coefficients, so the Green's function depends only on the separations $\mathbf{r}-\mathbf{r}'$ and $t-t'$ and we may put the source at the origin of space and time. Neither operator has a unique Green's function, and the condition that selects one is causality: we require

G(\mathbf{r},t) = 0 \qquad\text{for } t<0 ,

so that the response never precedes the impulse. The resulting $G$ is the causal or retarded Green's function. Reversing the requirement to $G=0$ for $t>0$ gives the advanced Green's function; Hassani obtains the two by displacing the poles of the $\omega$ integral downwards and upwards respectively, and notes that the advanced solution is not merely a curiosity, since relativistic quantum field theory builds its propagator from a combination of the two.

Unlike the Green's functions of the previous units, these are genuine distributions: Hassani stresses that in more than one dimension a Green's function need not be a function at all, and the wave case below contains a delta function outright.

Definition(Causal Green's function of an evolution operator)AdaptedAdapted from Hassani (2013), §22.4.3-22.4.4, the step function in Eqs. (22.45) and (22.48), pp. 685-686; Evans (2010), §2.3.1, the symbolic statement of the initial condition for $\\Phi$, p. 46

Let $L$ be $L_{\mathrm{heat}}$ or $L_{\mathrm{wave}}$ on $\mathbb{R}^n\times\mathbb{R}$ . A distribution $G$ is the causal Green's function of $L$ if $LG = \delta^{(n)}(\mathbf{r})\delta(t)$ and $G = 0$ for $t<0$ . Equivalently, and more usably:

for $L_{\mathrm{heat}}$ : $G$ vanishes for $t<0$ , solves $G_t = \kappa\nabla^2 G$ for $t>0$ , and $G(\cdot\,,t)\to\delta^{(n)}$ as $t\downarrow0$ ;
for $L_{\mathrm{wave}}$ : $G$ vanishes for $t<0$ , solves $G_{tt} = c^2\nabla^2 G$ for $t>0$ , and $G(\cdot\,,t)\to0$ , $G_t(\cdot\,,t)\to\delta^{(n)}$ as $t\downarrow0$ .

The second description is the one Evans uses when he writes the heat kernel's initial condition as $\Phi = \delta_0$ on $\mathbb{R}^n\times\{t=0\}$ ; it says that a unit impulse in a first-order equation sets the field itself, while in a second-order equation it sets the field's rate of change — exactly as for the one-dimensional causal Green's functions of ordinary differential operators, where the $p$ th-order operator has a jump in the $(p-1)$ st derivative.

Given $G$ , the solution of $Lu = f$ with zero data at $t=0$ is the space-time convolution

u(\mathbf{r},t) = \int_{-\infty}^{\infty}\!\!\int_{\mathbb{R}^n} G(\mathbf{r}-\mathbf{r}',\,t-t')\, f(\mathbf{r}',t')\; d^nr'\, dt' ,

in which causality truncates the $t'$ integral at $t' = t$ .

Method(The space-time Fourier transform)AdaptedAdapted from Hassani (2013), §22.4.3, Eqs. (22.42)-(22.43), p. 684, and §22.4.4, Eq. (22.46) and the pole displacement, pp. 685-686

Transform space, keep time. Set $\tilde G(\mathbf{q},t) = \int G(\mathbf{r},t)e^{-i\mathbf{q}\cdot\mathbf{r}}d^nr$ . The spatial Laplacian becomes $-q^2$ , so each Fourier mode obeys an ordinary differential equation in $t$ with a delta forcing: $\tilde G_t + \kappa q^2\tilde G = \delta(t) \qquad\text{or}\qquad \tilde G_{tt} + c^2q^2\tilde G = \delta(t) ,$ to be solved with $\tilde G = 0$ for $t<0$ . These are the causal Green's functions of a relaxation and of an undamped oscillator, and are written down in Causal Green's functions for initial-value problems.
Or transform time as well. Then $\partial_t\mapsto -i\omega$ and the equation becomes algebraic: $\tilde G = (\kappa q^2 - i\omega)^{-1}$ or $\tilde G = -(\omega^2-c^2q^2)^{-1}$ . The first has its single pole at $\omega = -i\kappa q^2$ , already in the lower half-plane, so no prescription is needed and causality is automatic. The second has two poles on the real axis, and causality must be imposed by pushing both into the lower half-plane; closing the $\omega$ contour in the lower half-plane for $t>0$ then catches them, while for $t<0$ the contour closes above and encircles nothing, which is where $G=0$ for $t<0$ comes from.
Invert the spatial transform, reducing the $\mathbf{q}$ integral to a radial one with the angular formula $\int e^{i\mathbf{q}\cdot\mathbf{r}}d\Omega_{\mathbf{q}} = 4\pi\sin(qr)/qr$ in three dimensions.
Convolve with the source, as in the definition.

Step 2 is Hassani's route, and he warns at it: taking the principal value of the $\omega$ integral for the wave operator — the apparently neutral choice — produces something that solves the wave equation with zero on the right-hand side, so that it is no Green's function at all. The contour is not a matter of taste.

Method(The heat kernel is the causal Green's function of the heat operator)AdaptedAdapted from Hassani (2013), §22.4.3, Eq. (22.45), p. 685; Evans (2010), §2.3.1, the definition of $\\Phi$ and Theorem 1, pp. 45-46

On $\mathbb{R}^n$ , the causal Green's function of $\partial_t - \kappa\nabla^2$ is

G(\mathbf{r},t) = \frac{\theta(t)}{(4\pi\kappa t)^{n/2}}\,\exp\!\left(-\frac{\lvert\mathbf{r}\rvert^2}{4\kappa t}\right) ,

$\theta$ being the Heaviside step. For $\kappa = 1$ this is Hassani's $G_s(r,t) = (4\pi t)^{-m/2}e^{-r^2/4t}\theta(t)$ in $m$ dimensions and Evans's fundamental solution $\Phi$ ; for $n=1$ it is the heat kernel. Consequently the solution of

u_t - \kappa\nabla^2 u = f \quad\text{in } \mathbb{R}^n\times(0,\infty), \qquad u = 0 \quad\text{on } \mathbb{R}^n\times\{t=0\}

u(\mathbf{r},t) = \int_0^t\!\!\int_{\mathbb{R}^n}\frac{\exp\bigl(-\lvert\mathbf{r}-\mathbf{r}'\rvert^2/4\kappa(t-s)\bigr)}{\bigl(4\pi\kappa(t-s)\bigr)^{n/2}}\; f(\mathbf{r}',s)\; d^nr'\, ds , \tag{H}

which is Evans's formula (13). He proves (his Theorem 2) that if $f$ has two continuous spatial derivatives and one continuous time derivative on $\mathbb{R}^n\times[0,\infty)$ and compact support in the space variable, then (H) is of class $C_1^2$ , satisfies the equation, and tends to $0$ as $t\downarrow0$ .

Idea of the proofOriginal

Each spatial Fourier mode of the heat equation decays independently, at the rate $\kappa q^2$ ; a unit impulse at $t=0$ excites every mode with amplitude $1$ , so the Green's function is the inverse transform of $e^{-\kappa q^2 t}$ — a Gaussian whose width in $\mathbf{q}$ shrinks like $t^{-1/2}$ and whose inverse transform is therefore a Gaussian whose width in $\mathbf{r}$ grows like $t^{1/2}$ . That is the whole content of the heat kernel, and the step function appears because the mode equation is first order, so the impulse creates the amplitude rather than changing it.

ProofAdaptedAdapted from Hassani (2013), §22.4.3, pp. 684-685; Evans (2010), §2.3.1, p. 45

Step 1 (the mode equation). By step 1 of the method, $\tilde G_t + \kappa q^2\tilde G = \delta(t)$ with $\tilde G=0$ for $t<0$ . Away from $t=0$ the solution is $Ce^{-\kappa q^2t}$ ; integrating the equation across $t=0$ shows that $\tilde G$ jumps by $1$ there, since $\int\kappa q^2\tilde G\, dt\to0$ over a shrinking interval. Hence

\tilde G(\mathbf{q},t) = \theta(t)\,e^{-\kappa q^2 t} .

Step 2 (invert). For $t>0$ , by the Gaussian integral in $n$ dimensions,

G(\mathbf{r},t) = \int\frac{d^nq}{(2\pi)^n}\,e^{i\mathbf{q}\cdot\mathbf{r}}e^{-\kappa t q^2} = \frac{1}{(2\pi)^n}\left(\frac{\pi}{\kappa t}\right)^{n/2}e^{-\lvert\mathbf{r}\rvert^2/4\kappa t} = \frac{e^{-\lvert\mathbf{r}\rvert^2/4\kappa t}}{(4\pi\kappa t)^{n/2}} ,

the middle step being the completion of the square $-\kappa tq^2+i\mathbf{q}\cdot\mathbf{r} = -\kappa t\lvert\mathbf{q}-i\mathbf{r}/2\kappa t\rvert^2 - \lvert\mathbf{r}\rvert^2/4\kappa t$ applied coordinate by coordinate, and the arithmetic of the constants being $(2\pi)^{-n}\pi^{n/2} = (2^n\pi^{n/2})^{-1}$ , so that the denominator is $(4\pi\kappa t)^{n/2}$ .

Step 3 (checks). $\int_{\mathbb{R}^n}G(\mathbf{r},t)\, d^nr = \tilde G(\mathbf{0},t) = 1$ for every $t>0$ : heat is conserved. And $G(\cdot\,,t)\to\delta^{(n)}$ as $t\downarrow0$ , since the mass is $1$ and concentrates in a ball of radius $O(\sqrt{\kappa t})$ . Direct differentiation confirms $G_t = \kappa\nabla^2G$ for $t>0$ .

Step 4 (the solution formula). By the definition's convolution formula, with the $t'$ integral cut at $t'=t$ by the step function and at $t'=0$ because $f$ is given only for $t\ge0$ , we get (H). $\square$

Method(The retarded Green's function of the three-dimensional wave equation)AdaptedAdapted from Hassani (2013), §22.4.4, Eq. (22.53) with $m=3$, p. 687; Evans (2010), §2.4.2, Eq. (44), p. 79

On $\mathbb{R}^3$ , the causal Green's function of $\partial_t^2 - c^2\nabla^2$ is

G(\mathbf{r},t) = \frac{\delta\bigl(t - \lvert\mathbf{r}\rvert/c\bigr)}{4\pi c^2\,\lvert\mathbf{r}\rvert} \;=\;\frac{\delta\bigl(\lvert\mathbf{r}\rvert - ct\bigr)}{4\pi c\,\lvert\mathbf{r}\rvert}\;\theta(t) .

It is supported entirely on the forward light cone $\lvert\mathbf{r}\rvert = ct$ . For $c=1$ this is Hassani's $G^{(\mathrm{ret})}_s(r,t) = \delta(t-r)/(4\pi r)$ . The advanced Green's function is the same expression with $t\mapsto -t$ , supported on the backward cone.

Consequently the solution of $u_{tt} - c^2\nabla^2 u = f$ with $u = u_t = 0$ at $t=0$ is the retarded potential

u(\mathbf{r},t) = \frac{1}{4\pi c^2}\int_{\lvert\mathbf{r}-\mathbf{r}'\rvert < ct}\frac{f\bigl(\mathbf{r}',\, t - \lvert\mathbf{r}-\mathbf{r}'\rvert/c\bigr)}{\lvert\mathbf{r}-\mathbf{r}'\rvert}\; d^3r' , \tag{W}

which for $c=1$ is Evans's formula (44). Each volume element of the source contributes as a static source would, but with its strength evaluated at the earlier time $t - \lvert\mathbf{r}-\mathbf{r}'\rvert/c$ at which a signal leaving it would have had to depart to arrive now.

Idea of the proofOriginal

Fourier transform in space again, and each mode is now an undamped oscillator of frequency $cq$ driven by an impulse. The causal response of such an oscillator is $\theta(t)\sin(cqt)/cq$ , so the Green's function is the inverse transform of that — and the inverse transform of $\sin(cqt)/cq$ in three dimensions is not a smooth function but a shell: the radial integral turns into $\int_0^\infty\sin(qr)\sin(cqt)\, dq$ , which is a delta function in $r-ct$ .

Why a shell, and not a blob as in the heat case? Because the modes all travel at the same speed $c$ and therefore stay in step: a wave packet assembled at $t=0$ from all wavenumbers with equal amplitude and the same speed remains a sharp front. In the heat equation the decay rate $\kappa q^2$ differs between modes and the packet smears. The sharpness of the front in three dimensions is Huygens' principle, and Evans notes that it holds in odd dimensions $n\ge3$ but fails in even ones, where a disturbance leaves a trailing wake behind the front.

ProofAdaptedAdapted from Hassani (2013), §22.4.4, Eqs. (22.46)-(22.48) and (22.53), pp. 685-687; Evans (2010), §2.4.2, p. 79

Step 1 (the mode equation). By step 1 of the method, $\tilde G_{tt} + c^2q^2\tilde G = \delta(t)$ with $\tilde G=0$ for $t<0$ . The homogeneous solutions are $\cos(cqt)$ and $\sin(cqt)$ ; causality forces $\tilde G$ and $\tilde G_t$ to vanish just before $0$ , and integrating the equation across $t=0$ shows that $\tilde G$ stays continuous while $\tilde G_t$ jumps by $1$ . The combination with $\tilde G(0)=0$ and $\tilde G_t(0^+)=1$ is

\tilde G(\mathbf{q},t) = \theta(t)\,\frac{\sin(cqt)}{cq} ,

which is the causal Green's function of an oscillator; Arfken, Weber and Harris exhibit the same function in one variable as $\theta(x-t)\sin(x-t)$ for $d^2/dx^2+1$ and point out that it is causality.

Step 2 (angular integration). For $t>0$ , taking the polar axis of $\mathbf{q}$ along $\mathbf{r}$ and writing $r=\lvert\mathbf{r}\rvert$ ,

G(\mathbf{r},t) = \frac{1}{(2\pi)^3}\int_0^\infty\!\! dq\; q^2\;\frac{\sin(cqt)}{cq}\cdot\frac{4\pi\sin(qr)}{qr} = \frac{1}{2\pi^2 c\,r}\int_0^\infty \sin(qr)\,\sin(cqt)\; dq .

Step 3 (the radial integral is a delta function). By the product formula $2\sin\alpha\sin\beta = \cos(\alpha-\beta)-\cos(\alpha+\beta)$ and the identity $\int_0^\infty\cos(qa)\, dq = \pi\delta(a)$ (Delta-function identities),

\int_0^\infty \sin(qr)\sin(cqt)\; dq = \frac12\int_0^\infty\Bigl[\cos\bigl(q(r-ct)\bigr) - \cos\bigl(q(r+ct)\bigr)\Bigr] dq = \frac{\pi}{2}\Bigl[\delta(r-ct) - \delta(r+ct)\Bigr] .

For $t>0$ and $r>0$ the argument $r+ct$ is strictly positive, so the second delta function contributes nothing and

G(\mathbf{r},t) = \frac{1}{2\pi^2 c r}\cdot\frac{\pi}{2}\,\delta(r-ct) = \frac{\delta(r-ct)}{4\pi c\,r} \qquad (t>0) .

Step 4 (rewrite on the cone). $\delta(r-ct) = c^{-1}\delta(t-r/c)$ , by the scaling rule for a delta function composed with $t\mapsto r-ct$ ; hence $G = \delta(t-r/c)/(4\pi c^2 r)$ , and this expression already vanishes for $t<0$ when $r>0$ , so the step function may be dropped. For $c=1$ and $m=3$ Hassani's general odd-dimensional formula reduces to the same thing.

Step 5 (the retarded potential). Convolving as in the definition,

u(\mathbf{r},t) = \int_0^\infty\!\!\int \frac{\delta\bigl(t-t'-\lvert\mathbf{r}-\mathbf{r}'\rvert/c\bigr)}{4\pi c^2\lvert\mathbf{r}-\mathbf{r}'\rvert} f(\mathbf{r}',t')\; d^3r'\, dt' ,

and the $t'$ integral is done by the delta function, which fires at $t' = t-\lvert\mathbf{r}-\mathbf{r}'\rvert/c$ . That value of $t'$ is non-negative only for $\lvert\mathbf{r}-\mathbf{r}'\rvert\le ct$ , which is the restriction on the spatial domain in (W). $\square$

Method(Duhamel's principle)AdaptedAdapted from Evans (2010), §2.3.1c, the family $(12_s)$ and Theorem 2, pp. 47-48, and §2.4.2, the family $(40_s)$, Eq. (41) and Theorem 4, pp. 77-78; Strauss (2008), §3.4, Eqs. (13)-(14) and the operator method, pp. 77-78

Duhamel's principle builds the solution of a forced problem out of solutions of the unforced one.

Heat equation. For $s\ge0$ let $u(\cdot\,,\cdot\,;s)$ solve the homogeneous initial-value problem started at time $s$ with initial datum $f(\cdot\,,s)$ : $u_t(\cdot\,;s) - \kappa\nabla^2u(\cdot\,;s) = 0 \text{ in } \mathbb{R}^n\times(s,\infty), \qquad u(\cdot\,;s) = f(\cdot\,,s) \text{ at } t=s .$ Then $u(\mathbf{r},t) = \int_0^t u(\mathbf{r},t;s)\, ds$ solves $u_t-\kappa\nabla^2u = f$ with $u=0$ at $t=0$ . Evans proves this (Theorem 2) for $f\in C_1^2(\mathbb{R}^n\times[0,\infty))$ with compact support in $\mathbf{r}$ , and it is the same thing as formula (H).
Wave equation. For $s\ge0$ let $u(\cdot\,,\cdot\,;s)$ solve $u_{tt}(\cdot\,;s) - c^2\nabla^2u(\cdot\,;s) = 0 \text{ in } \mathbb{R}^n\times(s,\infty), \qquad u(\cdot\,;s) = 0,\quad u_t(\cdot\,;s) = f(\cdot\,,s) \text{ at } t=s .$ Then $u(\mathbf{r},t) = \int_0^t u(\mathbf{r},t;s)\, ds$ solves $u_{tt}-c^2\nabla^2u = f$ with $u = u_t = 0$ at $t=0$ . Evans proves this (Theorem 4) for $n\ge2$ and $f\in C^{[n/2]+1}(\mathbb{R}^n\times[0,\infty))$ .

Note the asymmetry: in the first-order problem the source is fed in as a value, in the second-order problem as a velocity. That is the same asymmetry as in the definition of the causal Green's function, and for the same reason.

Strauss presents the principle as an operator identity. The ordinary differential equation $\ddot u + A^2u = f$ with $u(0)=\varphi$ , $\dot u(0) = \psi$ has solution $u(t) = S'(t)\varphi + S(t)\psi + \int_0^t S(t-s)f(s)\, ds$ with $S(t) = A^{-1}\sin tA$ ; the partial differential equation has the same form with $S(t)$ the solution operator that maps an initial velocity to the solution, which for $u_{tt} = c^2u_{xx}$ on the line is $S(t)\psi = (2c)^{-1}\int_{x-ct}^{x+ct}\psi$ . His summary is the one to remember: a solver for the free problem is already a solver for the forced one.

Idea of the proofOriginal

A source acting over a time interval can be chopped into a succession of brief impulses. Each impulse, once it has been delivered, is no longer a source: from then on the system evolves freely from whatever state the impulse left it in. So the forced solution is the superposition over $s$ of free evolutions started at time $s$ from the state that the source deposited there — and "the state the source deposited" is $f(\cdot,s)\,ds$ of field for a first-order equation and of velocity for a second-order one.

Stated in terms of Green's functions there is nothing left to prove: Duhamel's formula is convolution with the causal Green's function, because the causal Green's function is by definition the free evolution from a point impulse. The two statements of the principle are what the formulas (H) and (W) look like when the spatial convolution is left as "the solution operator" instead of being written out.

Proof(Wave equation)AdaptedAdapted from Evans (2010), §2.4.2, proof of Theorem 4, step 2, p. 78

Assume the regularity Evans assumes, which makes $u(\cdot\,,\cdot\,;s)$ of class $C^2$ on $\mathbb{R}^n\times[s,\infty)$ and allows the differentiations below under the integral sign. Put $u(\mathbf{r},t) = \int_0^t u(\mathbf{r},t;s)\, ds$ .

Step 1 (first time derivative). Differentiating with respect to $t$ picks up a boundary term from the upper limit plus the integral of the derivative:

u_t(\mathbf{r},t) = u(\mathbf{r},t;t) + \int_0^t u_t(\mathbf{r},t;s)\, ds = \int_0^t u_t(\mathbf{r},t;s)\, ds ,

the boundary term vanishing because $u(\cdot\,;s) = 0$ at $t=s$ — this is the first place the initial condition of the family is used.

Step 2 (second time derivative). Differentiating again,

u_{tt}(\mathbf{r},t) = u_t(\mathbf{r},t;t) + \int_0^t u_{tt}(\mathbf{r},t;s)\, ds = f(\mathbf{r},t) + \int_0^t u_{tt}(\mathbf{r},t;s)\, ds ,

the boundary term now being $u_t(\cdot\,;s)$ at $t=s$ , which is $f(\cdot\,,s)$ by the second initial condition of the family. This is where the source enters, and it enters once, undifferentiated.

Step 3 (the Laplacian). The Laplacian acts in $\mathbf{r}$ only, so it passes through the $s$ integral, and each member of the family satisfies the homogeneous equation:

c^2\nabla^2 u(\mathbf{r},t) = \int_0^t c^2\nabla^2u(\mathbf{r},t;s)\, ds = \int_0^t u_{tt}(\mathbf{r},t;s)\, ds .

Step 4 (subtract). The two integrals are identical, so $u_{tt} - c^2\nabla^2u = f$ . Finally $u(\mathbf{r},0)=0$ because the integral is empty, and $u_t(\mathbf{r},0)=0$ by Step 1 for the same reason. $\square$

Proof(Heat equation, and its identification with (H))AdaptedAdapted from Evans (2010), §2.3.1c, pp. 47-48; Hassani (2013), §22.4.3, p. 685

Step 1 (what the family is). By Step 2 of the heat proof above, the solution of the homogeneous heat equation started at time $s$ with datum $g$ is the convolution $\bigl(G(\cdot\,,t-s)*g\bigr)(\mathbf{r})$ . With $g = f(\cdot\,,s)$ ,

u(\mathbf{r},t;s) = \int_{\mathbb{R}^n} G(\mathbf{r}-\mathbf{r}',\,t-s)\, f(\mathbf{r}',s)\; d^nr' ,

so integrating over $s\in(0,t)$ reproduces exactly formula (H). Duhamel's principle for the heat equation and the convolution formula (H) are therefore the same statement, written in two ways.

Step 2 (why it solves the problem). Formally the computation is Steps 1-4 of the wave proof with one derivative fewer: $u_t = u(\mathbf{r},t;t) + \int_0^t u_t\, ds = f(\mathbf{r},t) + \int_0^t u_t(\mathbf{r},t;s)\, ds$ , now with the boundary term supplying $f$ because the family carries the source as its initial value; and $\kappa\nabla^2u = \int_0^t\kappa\nabla^2u\, ds = \int_0^t u_t\, ds$ . Subtracting gives $u_t-\kappa\nabla^2u=f$ .

Step 3 (why this needs Evans's theorem). The formal computation differentiates under an integral whose integrand is singular at the endpoint $s=t$ , where $G(\cdot\,,t-s)$ blows up. Evans handles this by first moving the singularity onto $f$ — changing variables so that the formula reads $\int_0^t\int G(\mathbf{r}',s)f(\mathbf{r}-\mathbf{r}',t-s)\, d^nr'\, ds$ , after which the derivatives fall on the smooth factor $f$ — and then integrating by parts on $(\varepsilon,t)$ and letting $\varepsilon\to0$ . This is where his hypotheses on $f$ (two spatial and one time derivative, continuous, with compact support in $\mathbf{r}$ ) are spent. The wave proof above needs no such care, because there the regularity of the family members is supplied by the explicit solution formulas for the homogeneous problem. $\square$

Remark(What the two Green's functions say about the two equations)Original

Support, and propagation speed. The wave Green's function is supported on the single sphere $\lvert\mathbf{r}\rvert = ct$ : a disturbance reaches a point exactly once and then is gone. The heat Green's function is positive everywhere for every $t>0$ , however large $\lvert\mathbf{r}\rvert$ : Evans's "infinite propagation speed". These are the opposite extremes, and the contrast is made concrete in the second example below.
Dimension. The sharp front is special to odd dimensions $n\ge3$ . In two dimensions the Green's function is supported on the whole solid cone $\lvert\mathbf{r}\rvert\le ct$ rather than on its boundary, which is what Evans draws from the representation formulas: in even dimensions the data at a point affect the solution throughout the interior of the cone, so a disturbance has a wake. This is Huygens' principle and its failure. Hassani's reduction to delta functions, (22.52)-(22.53), is derived only for odd $m$ , the even case being left as one of his problems, so the even-dimensional statement here is Evans's and not his.
The one-dimensional wave case. There $G = \tfrac{1}{2c}\theta(ct-\lvert x\rvert)$ — a plateau, not a delta function — and (W) becomes the integral of $f$ over the backward characteristic triangle, which is Strauss's formula $u = (2c)^{-1}\iint_\Delta f$ . So even the "odd-dimensional" sharpness is a statement about $n\ge3$ : in one dimension the front is sharp but the Green's function is a step, and the wake is constant rather than absent.
Relation to the Helmholtz Green's function. Fourier transforming the retarded Green's function in time at frequency $\omega$ gives $\int G(\mathbf{r},t)\,e^{i\omega t}\, dt = \frac{e^{i\omega r/c}}{4\pi c^2 r} = -\frac{1}{c^2}\left(-\frac{e^{ikr}}{4\pi r}\right),\qquad k=\frac{\omega}{c},$ which is $-c^{-2}$ times the outgoing Helmholtz Green's function of Green's functions for the Helmholtz equation. The factor $-c^{-2}$ is exactly what the algebra of substituting $u = \psi e^{-i\omega t}$ into $u_{tt}-c^2\nabla^2u = f$ demands, as the first example below checks. So retarded and outgoing are the same condition seen in the two domains, and the contour prescription $k\to k+i\delta$ of the previous unit is the Fourier transform of the prescription "push both $\omega$ poles down" used here.
Where the delta function went. In the heat case the Green's function is a smooth function for $t>0$ ; in the wave case it is not a function at all. Hassani draws the general moral: in more than one dimension Green's functions must be handled as distributions. The practical consequence is that (W) must never be differentiated naively under the integral sign — one differentiates the retarded argument, which produces the extra terms responsible for radiation fields in electrodynamics.
What is and is not proved here. That the stated $G$ are the causal Green's functions is established above. That formulas (H) and (W) actually define solutions of the stated regularity is not; that is Evans's Theorems 2 and 4, whose hypotheses on $f$ are quoted in the statement, and we use them rather than reprove them.

Example(A steady spatial source in the heat equation, and the approach to the steady state)Original

Solve $u_t - u_{xx} = \cos x$ on $\mathbb{R}$ with $u(x,0)=0$ ( $n=1$ , $\kappa=1$ ).

The one integral needed. For $\tau>0$ , substituting $z = x-y$ and using that $G$ is even, $\int_{\mathbb{R}}G(x-y,\tau)\cos y\; dy = \cos x\int_{\mathbb{R}}G(z,\tau)\cos z\; dz = \cos x\,\operatorname{Re}\int\frac{e^{-z^2/4\tau}}{\sqrt{4\pi\tau}}\,e^{iz}\, dz = e^{-\tau}\cos x ,$ the odd part dropping out and the Gaussian integral giving $e^{-\tau}$ (the Fourier transform of $G$ at wavenumber $1$ , which is $e^{-\kappa q^2\tau}$ with $\kappa=q=1$ — exactly Step 1 of the proof).
Duhamel. By (H), $u(x,t) = \int_0^t e^{-(t-s)}\cos x\; ds = \cos x\int_0^t e^{-\tau}\, d\tau = \bigl(1-e^{-t}\bigr)\cos x .$
Check the equation. $u_t = e^{-t}\cos x$ and $u_{xx} = -(1-e^{-t})\cos x$ , so $u_t - u_{xx} = e^{-t}\cos x + \cos x - e^{-t}\cos x = \cos x ,$ and $u(x,0)=0$ . Both conditions hold.
Numbers and limits. At $x=0$ , $t=\ln 2$ : $u = 1-\tfrac12 = 0.5$ exactly. As $t\to\infty$ , $u\to\cos x$ , which is the steady state: it satisfies $-u_{xx} = \cos x$ . The approach is exponential with rate $1 = \kappa q^2$ for $q=1$ , so a source of spatial wavelength $2\pi/q$ equilibrates in a time $\sim 1/\kappa q^2$ — the diffusive scaling, and the reason fine spatial structure in a temperature field relaxes first.

Example(A point source switched on: light cone versus diffusion)Original

Compare the responses of the two equations in three dimensions to the same source, a point source of constant strength switched on at $t=0$ :

f(\mathbf{r},t) = Q\,\delta^{(3)}(\mathbf{r})\,\theta(t) .

Wave equation. Formula (W) with the delta function in space gives $u(\mathbf{r},t) = \frac{Q}{4\pi c^2 r}\,\theta\!\left(t-\frac{r}{c}\right) .$ Nothing at all happens at distance $r$ until the moment $t=r/c$ ; from then on the field is frozen at the static value $Q/(4\pi c^2 r)$ , which is the solution of $-c^2\nabla^2u = Q\delta^{(3)}$ . With $Q = 4\pi$ , $c=1$ : at $r=2$ the field is $0$ for $t<2$ and $Q/(4\pi c^2 r) = 4\pi/(4\pi\cdot 2) = 0.5$ for $t>2$ .
Heat equation. Formula (H) with the same source gives, with $\tau = t-s$ , $u(\mathbf{r},t) = Q\int_0^t\frac{e^{-r^2/4\kappa\tau}}{(4\pi\kappa\tau)^{3/2}}\; d\tau .$ Substitute $\xi = r/(2\sqrt{\kappa\tau})$ , so that $\tau = r^2/(4\kappa\xi^2)$ and $d\tau = -r^2\,d\xi/(2\kappa\xi^3)$ , while $(4\pi\kappa\tau)^{3/2} = \pi^{3/2}r^3/\xi^3$ . The powers of $\xi$ cancel exactly and $u(\mathbf{r},t) = \frac{Q}{2\pi^{3/2}\kappa r}\int_{r/2\sqrt{\kappa t}}^{\infty}e^{-\xi^2}\, d\xi = \frac{Q}{4\pi\kappa r}\operatorname{erfc}\!\left(\frac{r}{2\sqrt{\kappa t}}\right) ,$ using $\int_a^\infty e^{-\xi^2}d\xi = \tfrac{\sqrt\pi}{2}\operatorname{erfc}(a)$ .
Checks. As $t\to\infty$ the complementary error function tends to $1$ and $u\to Q/(4\pi\kappa r)$ , the solution of $-\kappa\nabla^2u = Q\delta^{(3)}$ : the same static limit as in part 1 with $\kappa$ in place of $c^2$ . At short times $\operatorname{erfc}(a) \sim e^{-a^2}/(a\sqrt\pi)$ , so $u$ is exponentially small but not zero at any $r$ — the contrast with part 1, where it is exactly zero outside the cone.
A number. Take $Q = 4\pi$ , $\kappa = 1$ , $r = 1$ , $t = \tfrac14$ . Then the argument is $r/2\sqrt{\kappa t} = 1$ and $u = \operatorname{erfc}(1) = 0.1572992 ,$ against the eventual steady value $1$ . The front of the wave problem, by contrast, would not have arrived at $r=1$ until $t=1/c$ .
Diffusive versus ballistic. The error-function argument shows that the heat response at distance $r$ reaches half its final value when $r/2\sqrt{\kappa t} = \operatorname{erfc}^{-1}(1/2) = 0.4769$ , that is at $t \approx 1.10\,r^2/\kappa$ — a time quadratic in $r$ , not linear. That is the whole difference between diffusion and propagation, and it is readable straight off the two Green's functions.

Example(Time-harmonic source: consistency with the Helmholtz Green's function)Original

Take $f(\mathbf{r},t) = q\,\delta^{(3)}(\mathbf{r})\cos\omega t$ , a point source oscillating for all time (so we drop the initial conditions and use the Green's function directly). Convolving,

u(\mathbf{r},t) = \frac{q}{4\pi c^2 r}\cos\bigl(\omega t - k r\bigr), \qquad k = \frac{\omega}{c} ,

an outgoing spherical wave whose phase lags the source by $kr$ : the retardation made visible.

The same answer must follow from the previous unit. Writing $u = \operatorname{Re}\bigl[\psi(\mathbf{r})e^{-i\omega t}\bigr]$ and $f = \operatorname{Re}\bigl[q\delta^{(3)}e^{-i\omega t}\bigr]$ in $u_{tt}-c^2\nabla^2u = f$ gives $-\omega^2\psi - c^2\nabla^2\psi = q\delta^{(3)}$ , that is

\bigl(\nabla^2 + k^2\bigr)\psi = -\frac{q}{c^2}\,\delta^{(3)} ,

so $\psi = -\dfrac{q}{c^2}\,G_{\mathrm{Helmholtz}} = -\dfrac{q}{c^2}\left(-\dfrac{e^{ikr}}{4\pi r}\right) = \dfrac{q\,e^{ikr}}{4\pi c^2 r}$ , whose real part against $e^{-i\omega t}$ is the formula above. The two units agree, including the sign, and the agreement is a check on both conventions at once. With $q = 4\pi$ , $c=1$ , $\omega=1$ and $r=\pi$ : $u(\mathbf{r},0) = \cos(-\pi)/\pi = -1/\pi = -0.3183099$ .

ApplicationOriginal

Electrodynamics. In the Lorenz gauge each component of the four-potential satisfies $\Box A^\mu = \mu_0 J^\mu$ with $\Box = c^{-2}\partial_t^2-\nabla^2$ , so (W) gives the retarded potentials

\phi(\mathbf{r},t) = \frac{1}{4\pi\varepsilon_0}\int\frac{\rho_{\mathrm e}(\mathbf{r}',t-\lvert\mathbf{r}-\mathbf{r}'\rvert/c)}{\lvert\mathbf{r}-\mathbf{r}'\rvert}\, d^3r', \qquad \mathbf{A}(\mathbf{r},t) = \frac{\mu_0}{4\pi}\int\frac{\mathbf{J}(\mathbf{r}',t-\lvert\mathbf{r}-\mathbf{r}'\rvert/c)}{\lvert\mathbf{r}-\mathbf{r}'\rvert}\, d^3r' ,

from which the Liénard-Wiechert fields of a moving charge follow by differentiating the retarded argument — and the terms generated by that differentiation, proportional to the acceleration and falling off as $1/r$ , are the radiation fields. The choice of the retarded rather than the advanced Green's function is the arrow of time in classical electrodynamics; the symmetric combination is what Feynman-Wheeler absorber theory and the Feynman propagator of field theory are built on, which is Hassani's point about the advanced solution not being discardable.

Heat sources and thermal design. The heat Green's function is the workhorse of transient thermal analysis: the temperature rise under a laser spot, a buried power cable, a frictional contact or a weld pool is computed as a convolution of the source history with the Gaussian kernel, and the error-function profile of the second example is the standard "continuous point source" solution used in geothermal borehole design and in machining-temperature estimates. In geophysics the same kernel with $\kappa\sim10^{-6}\,\mathrm{m^2\,s^{-1}}$ gives the diurnal and annual skin depths of the ground temperature, and in semiconductor processing it sets the thermal budget of an anneal.

Duhamel in engineering practice. Structural and mechanical engineering uses Duhamel's integral under that name: the response of a linear oscillator to an arbitrary load is $\int_0^t h(t-s)F(s)\, ds$ with $h$ the unit-impulse response, which is exactly the one-mode case of the statement above and is what seismic response spectra are computed from. In control theory the same formula is the convolution of the input with the impulse response, and the step-response form is the $\theta(t-r/c)$ of the second example. The content in every case is Strauss's: solving the free problem once suffices.

Acoustics and seismology. Convolution with the retarded Green's function is forward modelling: a seismogram is the Earth's Green's function convolved with the source time function, and the central task of exploration seismology is to deconvolve the one from the other. The sharpness of the three-dimensional kernel — Huygens' principle — is what makes travel-time tomography possible at all, since it is what lets an arrival be assigned a single ray path; in two-dimensional modelling the wake of the even-dimensional Green's function has to be corrected for explicitly.

Diffusion-limited processes and probability. The heat Green's function is the transition density of Brownian motion, and (H) reads as the expected accumulated source along Brownian paths; this is the Feynman-Kac link that makes Brownian motion and the heat equation two views of one object, and it is why Monte Carlo path sampling is a legitimate solver for the heat equation with sources.

Application · advanced

Scattering and the Born approximation

Adapted

Add to PDF

Turns the stationary scattering problem into the Lippmann-Schwinger integral equation by inverting $\nabla^2+k^2$ with the outgoing Green's function, reads the scattering amplitude off the far field, identifies its Neumann series as the Born series, and evaluates the first Born amplitude $-\tfrac{1}{4\pi}\mathcal{F}U(\mathbf q)$ at momentum transfer $q=2k\sin(\theta/2)$ on a Gaussian potential and on a one-dimensional delta potential where the whole series can be summed.

Prerequisites: Green's functions for the Helmholtz equation, Neumann series, The Helmholtz equation, Fourier transform in several dimensions

Note(Conventions)AdaptedAdapted from Arfken, Weber & Harris (2013), §10.2, Example 10.2.2, Eqs. (10.45)-(10.47) and Table 10.1, pp. 463, 465; Strauss (2008), §13.3, Eqs. (7)-(10) and the paragraph on the radiation condition, pp. 367-368

Four conventions are fixed once and used throughout.

Time factor and the meaning of "outgoing". We suppress a time dependence $e^{-i\omega t}$ with $\omega = \hbar k^2/2m$ , so that an outgoing spherical wave carries $e^{+ikr}$ . This is the convention of Arfken, Weber and Harris, and it is the one Green's functions for the Helmholtz equation uses. Strauss suppresses $e^{+i\omega t}$ instead; in his book an outgoing spherical wave carries $e^{-ikr}$ and his radiation condition reads $r(\partial v/\partial r+ikv)\to0$ , with a plus sign where ours has a minus. The two conventions are complex conjugates of each other and describe the same physics; mixing them silently is the one thing that must not be done, so every formula below is in the $e^{-i\omega t}$ convention, and Strauss's results are conjugated before being used.

Reduced potential. The time-independent Schrodinger equation $-\tfrac{\hbar^2}{2m}\nabla^2\psi+V\psi = E\psi$ is rearranged, following Arfken, Weber and Harris, into a Helmholtz equation with a source on the right:

\bigl(\nabla^2+k^2\bigr)\psi = U\psi , \qquad U = \frac{2m}{\hbar^2}\,V , \qquad k^2 = \frac{2mE}{\hbar^2} \;>\;0 .

The reduced potential $U$ has the dimensions of an inverse squared length; nothing below needs $\hbar$ or $m$ again.

Green's function. We use $\bigl(\nabla^2+k^2\bigr)G = \delta$ , so that in three dimensions

G(\mathbf{r}) = -\frac{e^{ikr}}{4\pi r}, \qquad r = \lvert\mathbf{r}\rvert ,

as in Green's functions for the Helmholtz equation. Riley, Hobson and Bence write the same object with the operator $-\nabla^2-K^2$ and therefore quote it with the opposite sign, $+e^{iKr}/(4\pi r)$ ; the two agree.

Fourier transform. Scattering formulas are cleanest with no prefactor on the transform, so throughout this unit

\mathcal{F}U(\mathbf{q}) \;=\; \int e^{-i\mathbf{q}\cdot\mathbf{r}}\,U(\mathbf{r})\,\mathrm{d}^{3}r .

This is not the transform of Fourier transform in several dimensions, which is the symmetric one, $\tilde U(\mathbf{q}) = (2\pi)^{-3/2}\int e^{-i\mathbf{q}\cdot\mathbf{r}}U\,\mathrm{d}^{3}r$ , and for which this handbook reserves the tilde. The two differ by a fixed factor,

\mathcal{F}U = (2\pi)^{3/2}\,\tilde U , \qquad (2\pi)^{3/2} = 15.7496\ldots ,

so any formula taken from that unit must be multiplied by $(2\pi)^{3/2}$ before it is used here, and any value computed here must be divided by it before it is compared with one quoted there. The two places below where that unit is cited for a transform say so explicitly. The unnormalised choice is what makes the first Born amplitude exactly $-\mathcal{F}U/4\pi$ , with no stray powers of $2\pi$ .

Definition(The scattering set-up, the amplitude and the cross-sections)AdaptedAdapted from Arfken, Weber & Harris (2013), §10.2, Example 10.2.2, Eqs. (10.46)-(10.47) and the paragraph after them, pp. 465; Riley, Hobson & Bence (2006), Problem 4.36 and its closing note, p. 149; Strauss (2008), §13.3, opening paragraph, p. 366

A scatterer is a reduced potential $U$ that is measurable, bounded, and supported in the ball $\lvert\mathbf{r}\rvert\le a$ . Fix an energy through $k>0$ and a beam direction through a unit vector $\hat{\mathbf{e}}_{0}$ , and put $\mathbf{k}_{0} = k\hat{\mathbf{e}}_{0}$ . The incident wave is the plane wave

\psi_{\mathrm{inc}}(\mathbf{r}) = e^{i\mathbf{k}_{0}\cdot\mathbf{r}},

which solves $(\nabla^2+k^2)\psi_{\mathrm{inc}} = 0$ precisely because $\lvert\mathbf{k}_{0}\rvert = k$ . The same $k$ appears in the incident and the scattered wave because the equation contains nothing that could absorb energy: the scattering is elastic, which is the reading Arfken, Weber and Harris give of their Eq. (10.45).

A scattering solution is a bounded continuous $\psi$ with

\psi(\mathbf{r}) = e^{i\mathbf{k}_{0}\cdot\mathbf{r}} + f(\hat{\mathbf{n}})\,\frac{e^{ikr}}{r} + O\!\left(r^{-2}\right), \qquad r\to\infty,\quad \hat{\mathbf{n}} = \mathbf{r}/r ,

the error being uniform in $\hat{\mathbf{n}}$ . The coefficient $f$ is the scattering amplitude; it has the dimensions of a length and depends on $k$ and on $\hat{\mathbf{n}}$ , hence on the scattering angle $\theta$ between $\hat{\mathbf{n}}$ and $\hat{\mathbf{e}}_{0}$ and on the azimuth $\varphi$ . The factor $1/r$ is forced by conservation: it makes the flux carried by the scattered wave through a large sphere independent of its radius.

Following Riley, Hobson and Bence, the differential cross-section and the total cross-section are

\frac{\mathrm{d}\sigma}{\mathrm{d}\Omega} = \lvert f(\theta,\varphi)\rvert^{2}, \qquad \sigma = \int \lvert f\rvert^{2}\,\mathrm{d}\Omega = 2\pi\int_{0}^{\pi}\lvert f(\theta)\rvert^{2}\sin\theta\,\mathrm{d}\theta \quad(\text{central } U).

Arfken, Weber and Harris call $\lvert f\rvert^{2}$ itself the scattering cross section, meaning the differential one.

Finally, the momentum transfer is the vector

\mathbf{q} = k\hat{\mathbf{n}} - \mathbf{k}_{0}, \qquad q = \lvert\mathbf{q}\rvert = 2k\sin\frac{\theta}{2},

the second equality because $\lvert k\hat{\mathbf{n}}-\mathbf{k}_{0}\rvert^{2} = 2k^{2}(1-\cos\theta)$ . Riley, Hobson and Bence call this quantity $K$ . Note the range: $0\le q\le 2k$ . An experiment at wavenumber $k$ therefore samples the scatterer only at momentum transfers of magnitude at most $2k$ , that is only at spatial periods $2\pi/q$ down to $\pi/k$ and no finer; the application block reads this as the resolution limit.

Application(Lippmann-Schwinger equation, Born series, first Born amplitude)AdaptedAdapted from Arfken, Weber & Harris (2013), §10.2, Example 10.2.2, Eqs. (10.48)-(10.51), pp. 465-466; Arfken, Weber & Harris (2013), §21.3, p. 1067; Riley, Hobson & Bence (2006), Problem 4.36, p. 149

Let $U$ be a scatterer and $\mathbf{k}_{0},k$ as in the definition. Write

(\mathcal{T}\phi)(\mathbf{r}) \;=\; -\frac{1}{4\pi}\int_{\lvert\mathbf{r}'\rvert\le a} \frac{e^{ik\lvert\mathbf{r}-\mathbf{r}'\rvert}}{\lvert\mathbf{r}-\mathbf{r}'\rvert}\,U(\mathbf{r}')\,\phi(\mathbf{r}')\,\mathrm{d}^{3}r' ,

the convolution with the outgoing Green's function of the source $U\phi$ .

The integral equation. Every bounded continuous solution of the Lippmann-Schwinger equation
$\psi = \psi_{\mathrm{inc}} + \mathcal{T}\psi \tag{LS}$
satisfies $(\nabla^{2}+k^{2})\psi = U\psi$ in the sense of distributions on $\mathbb{R}^{3}$ , and is a scattering solution in the sense of the definition. In the other direction, (LS) is what one gets by inverting $\nabla^{2}+k^{2}$ with the outgoing Green's function and adding to the result the homogeneous solution $\psi_{\mathrm{inc}}$ ; Arfken, Weber and Harris obtain it exactly this way and call the resulting equation exact.
The amplitude. For a solution of (LS) the amplitude of the definition exists and is
$f(\hat{\mathbf{n}}) = -\frac{1}{4\pi}\int e^{-ik\hat{\mathbf{n}}\cdot\mathbf{r}'}\,U(\mathbf{r}')\,\psi(\mathbf{r}')\,\mathrm{d}^{3}r' . \tag{A}$
This is exact, but it is not a formula for $f$ : the unknown $\psi$ stands inside it.
The Born series. Equation (LS) is a Fredholm equation of the second kind with parameter $1$ and kernel $G(\mathbf{r}-\mathbf{r}')U(\mathbf{r}')$ . Its Neumann series (Neumann series) is the Born series
$\psi = \sum_{n\ge0}\mathcal{T}^{\,n}\psi_{\mathrm{inc}} = \psi_{\mathrm{inc}} + \mathcal{T}\psi_{\mathrm{inc}} + \mathcal{T}^{2}\psi_{\mathrm{inc}} + \cdots ,$
and whenever it converges uniformly on the support of $U$ it may be inserted in (A), giving $f = \sum_{n\ge1}f^{(n)}$ with $f^{(n)}$ built from $\mathcal{T}^{\,n-1}\psi_{\mathrm{inc}}$ . Arfken, Weber and Harris record the general rule for such a series: the condition $\lvert\lambda\rvert\lvert K\rvert_{\max}\lvert b-a\rvert<1$ is sufficient but crude, and the sharp condition is that $\lvert\lambda\rvert$ be smaller than the smallest characteristic value of the homogeneous equation.
The first Born approximation. Truncating after the first correction, that is replacing $\psi$ by $\psi_{\mathrm{inc}}$ inside (A), gives
$f_{\mathrm{B}}(\mathbf{q}) = -\frac{1}{4\pi}\int e^{-i\mathbf{q}\cdot\mathbf{r}'}\,U(\mathbf{r}')\,\mathrm{d}^{3}r' = -\frac{1}{4\pi}\,\mathcal{F}U(\mathbf{q}), \qquad \mathbf{q} = k\hat{\mathbf{n}}-\mathbf{k}_{0} . \tag{B}$
The amplitude in this approximation depends on $\hat{\mathbf{n}}$ and $k$ only through the single variable $\mathbf{q}$ . If moreover $U$ is central, $U = U(r)$ , then $f_{\mathrm{B}}$ depends only on $q = 2k\sin(\theta/2)$ and
$f_{\mathrm{B}}(q) = -\frac{1}{q}\int_{0}^{\infty} r\,U(r)\,\sin(qr)\,\mathrm{d}r , \tag{B$'$}$
which is real whenever $U$ is real, as step 5 of the proof records; the definition of a scatterer does not require reality, and for a complex $U$ — an absorbing medium, say — the right-hand side of (B $'$ ) is complex. Riley, Hobson and Bence quote the special case of a uniform well of depth $V_{0}$ and radius $b$ ; our (B $'$ ) with $U = -2mV_{0}/\hbar^{2}$ on $r<b$ reproduces their $f(\theta) = (2mV_{0}/\hbar^{2}K^{3})(\sin Kb-Kb\cos Kb)$ exactly, including the sign, which is the check that fixes all of our conventions together.

Idea of the proofOriginal

A scattering experiment is an inhomogeneous Helmholtz problem in disguise. Move the potential term to the right-hand side and it becomes $(\nabla^{2}+k^{2})\psi = U\psi$ : a known operator acting on $\psi$ , driven by a source that happens to contain $\psi$ itself. Inverting the operator with the outgoing Green's function therefore does not solve the problem, it converts it into an integral equation; the incident beam enters not as data but as the homogeneous solution one is free to add. Two things then follow almost for free. First, because the source is confined to a ball, the far field of the scattered part is a single outgoing spherical wave whose coefficient is the Fourier transform of the source evaluated at the outgoing wavevector, so the amplitude is read off by one Taylor expansion of $\lvert\mathbf{r}-\mathbf{r}'\rvert$ . Second, the integral equation is of the second kind, so it may be iterated: guess that the particle has not been scattered, feed that guess back in, repeat. The first round is the Born approximation, and the amplitude it produces is just the Fourier transform of the potential at the momentum the particle lost.

ProofAdaptedAdapted from Arfken, Weber & Harris (2013), §10.2, Example 10.2.2, the construction of Eqs. (10.48)-(10.51), pp. 465-466

Throughout, $U$ is bounded and supported in $\lvert\mathbf{r}\rvert\le a$ , and $\psi$ is bounded and continuous, so every integral written converges absolutely: the only singularity of the kernel is the integrable $\lvert\mathbf{r}-\mathbf{r}'\rvert^{-1}$ , and the domain of integration is a ball.

Step 1: a solution of (LS) solves the differential equation. By definition $\mathcal{T}\psi = G*(U\psi)$ with $G$ the outgoing Green's function of the note. Applying $\nabla^{2}+k^{2}$ and using $(\nabla^{2}+k^{2})G = \delta$ (that is Green's functions for the Helmholtz equation), the convolution theorem for distributions gives

\bigl(\nabla^{2}+k^{2}\bigr)\mathcal{T}\psi = \delta * (U\psi) = U\psi .

The other term contributes nothing, since $(\nabla^{2}+k^{2})e^{i\mathbf{k}_{0}\cdot\mathbf{r}} = (k^{2}-\lvert\mathbf{k}_{0}\rvert^{2})e^{i\mathbf{k}_{0}\cdot\mathbf{r}} = 0$ — here, and only here, the hypothesis $\lvert\mathbf{k}_{0}\rvert = k$ is used. Hence $(\nabla^{2}+k^{2})\psi = U\psi$ .

Step 2: the far field, and part 2. Let $r = \lvert\mathbf{r}\rvert > 2a$ and $\lvert\mathbf{r}'\rvert\le a$ . Then

\lvert\mathbf{r}-\mathbf{r}'\rvert = r\sqrt{1-\frac{2\hat{\mathbf{n}}\cdot\mathbf{r}'}{r}+\frac{\lvert\mathbf{r}'\rvert^{2}}{r^{2}}} = r - \hat{\mathbf{n}}\cdot\mathbf{r}' + R(\mathbf{r},\mathbf{r}'), \qquad \lvert R\rvert \le \frac{2a^{2}}{r},

the bound coming from expanding the square root and estimating the remainder, which is legitimate because $\lvert\mathbf{r}'\rvert/r<\tfrac12$ . Consequently

e^{ik\lvert\mathbf{r}-\mathbf{r}'\rvert} = e^{ikr}e^{-ik\hat{\mathbf{n}}\cdot\mathbf{r}'}\bigl(1+O(ka^{2}/r)\bigr), \qquad \frac{1}{\lvert\mathbf{r}-\mathbf{r}'\rvert} = \frac{1}{r}\bigl(1+O(a/r)\bigr),

with constants depending only on $k$ and $a$ , uniformly in $\hat{\mathbf{n}}$ and in $\mathbf{r}'$ in the support. Substituting into $\mathcal{T}\psi$ and using that $\lvert U\psi\rvert$ is bounded on a set of finite volume,

(\mathcal{T}\psi)(\mathbf{r}) = \frac{e^{ikr}}{r}\left[-\frac{1}{4\pi}\int e^{-ik\hat{\mathbf{n}}\cdot\mathbf{r}'}U(\mathbf{r}')\psi(\mathbf{r}')\,\mathrm{d}^{3}r'\right] + O\!\left(r^{-2}\right).

Adding $\psi_{\mathrm{inc}}$ gives exactly the asymptotic form required of a scattering solution, with the bracket as $f(\hat{\mathbf{n}})$ . That is part 2, and it also finishes part 1: the scattered part is a multiple of $e^{ikr}/r$ to leading order, so $\partial_{r}(\mathcal{T}\psi)-ik\,\mathcal{T}\psi = O(r^{-2}) = o(r^{-1})$ , which is the Sommerfeld outgoing condition of The Helmholtz equation.

Step 3: the series, part 3. Iterating (LS) in the manner of Neumann series, with $\lambda = 1$ and the kernel $K(\mathbf{r},\mathbf{r}') = G(\mathbf{r}-\mathbf{r}')U(\mathbf{r}')$ , produces the partial sums $\sum_{n=0}^{N}\mathcal{T}^{\,n}\psi_{\mathrm{inc}}$ ; these are the Born approximations of successive orders. If they converge uniformly on the support of $U$ then the limit may be taken inside the integral in (A), because the support has finite volume and the remaining factor $e^{-ik\hat{\mathbf{n}}\cdot\mathbf{r}'}U(\mathbf{r}')$ is bounded there; this splits $f$ into the stated sum. The convergence itself is not supplied by the theorem of Neumann series, whose hypotheses are not met here; see the remark.

Step 4: the first term, part 4. Putting $\psi = \psi_{\mathrm{inc}}$ in (A) and combining the two exponentials,

f_{\mathrm{B}} = -\frac{1}{4\pi}\int e^{-ik\hat{\mathbf{n}}\cdot\mathbf{r}'}U(\mathbf{r}')e^{i\mathbf{k}_{0}\cdot\mathbf{r}'}\,\mathrm{d}^{3}r' = -\frac{1}{4\pi}\int e^{-i(k\hat{\mathbf{n}}-\mathbf{k}_{0})\cdot\mathbf{r}'}U(\mathbf{r}')\,\mathrm{d}^{3}r' ,

which is (B).

Step 5: central potentials. For $U = U(r)$ take the polar axis of the $\mathbf{r}'$ integration along $\mathbf{q}$ and write $\mu = \cos\angle(\mathbf{q},\mathbf{r}')$ . Then

\int e^{-i\mathbf{q}\cdot\mathbf{r}'}U(r')\,\mathrm{d}^{3}r' = 2\pi\int_{0}^{\infty}\!\!\int_{-1}^{1} e^{-iqr'\mu}\,U(r')\,r'^{2}\,\mathrm{d}\mu\,\mathrm{d}r' = 4\pi\int_{0}^{\infty} r'^{2}U(r')\,\frac{\sin qr'}{qr'}\,\mathrm{d}r' ,

using $\int_{-1}^{1}e^{-is\mu}\mathrm{d}\mu = 2\sin s/s$ . Dividing by $-4\pi$ gives (B $'$ ). The result is real whenever $U$ is real, and depends on the direction only through $q = 2k\sin(\theta/2)$ .

Step 6: the check against the uniform well. With $V = -V_{0}$ for $r<b$ and $0$ outside, so $U = -2mV_{0}/\hbar^{2}$ on $r<b$ , formula (B $'$ ) gives

f_{\mathrm{B}}(q) = \frac{2mV_{0}}{\hbar^{2}q}\int_{0}^{b} r\sin(qr)\,\mathrm{d}r = \frac{2mV_{0}}{\hbar^{2}q}\cdot\frac{\sin qb - qb\cos qb}{q^{2}} ,

because $\int_{0}^{b}r\sin(qr)\,\mathrm{d}r = [\,q^{-2}\sin qr - q^{-1}r\cos qr\,]_{0}^{b}$ . With $q$ written as $K$ this is the expression Riley, Hobson and Bence quote. $\square$

Remark(Where the hypotheses are used, and where the method breaks)Original

Elasticity. $\lvert\mathbf{k}_{0}\rvert = k$ is used in step 1, to make the incident plane wave a homogeneous solution, and already in the definition, where it is what reduces $q = \lvert k\hat{\mathbf{n}}-\mathbf{k}_{0}\rvert$ to $2k\sin(\theta/2)$ . Without it the plane wave would contribute a spurious source and (LS) would not solve the right equation.
Compact support. Used only in step 2, but used twice there: to make the error terms uniform, and to keep the far-field integral finite. It can be weakened to rapid decay, but not to arbitrary decay. The Coulomb potential $U\propto1/r$ fails it badly: the expansion of $\lvert\mathbf{r}-\mathbf{r}'\rvert$ is then applied where $\lvert\mathbf{r}'\rvert$ is comparable with $r$ , the scattered wave acquires a logarithmic phase, and no amplitude of the assumed form exists. It is a well-known accident that formula (B $'$ ), applied to a screened Coulomb potential and then unscreened, still produces the correct Rutherford cross-section; the application block does this calculation. The accident does not license the derivation.
Outgoing, not incoming. Replacing $G$ by its conjugate would give a perfectly good solution of the same differential equation with an incoming scattered wave. Nothing in the mathematics prefers one; the physical picture of a beam that arrives and is scattered away does. This is the same non-uniqueness that Green's functions for the Helmholtz equation resolves by the radiation condition, and the reason Strauss states his problem with a radiation condition attached.
Convergence of the Born series is not covered by our Neumann theorem. The theorem of Neumann series assumes a finite interval and a kernel bounded and continuous on the square. Here the domain is $\mathbb{R}^{3}$ , the kernel $G(\mathbf{r}-\mathbf{r}')U(\mathbf{r}')$ has an integrable but unbounded singularity on the diagonal, and on $L^{2}(\mathbb{R}^{3})$ it is not even Hilbert-Schmidt, since $\int\lvert G(\mathbf{r}-\mathbf{r}')\rvert^{2}\mathrm{d}^{3}r = (16\pi^{2})^{-1}\int r^{-2}\mathrm{d}^{3}r$ diverges. The divergence is at large $r$ only, and it is removed by restricting to the ball $\lvert\mathbf{r}\rvert\le a$ , which is the space the Born series is used on: part 3 of the statement and step 3 of the proof both ask only for convergence on the support of $U$ , and nothing more is needed, since (A) integrates over that support alone. On $L^{2}(\{\lvert\mathbf{r}\rvert\le a\})$ the kernel is Hilbert-Schmidt, because $\lvert G(\mathbf{r}-\mathbf{r}')\rvert^{2} = (16\pi^{2})^{-1}\lvert\mathbf{r}-\mathbf{r}'\rvert^{-2}$ has an integrable singularity and $\int_{\lvert\mathbf{r}'\rvert\le a}\int_{\lvert\mathbf{r}\rvert\le a}\lvert\mathbf{r}-\mathbf{r}'\rvert^{-2}\lvert U(\mathbf{r}')\rvert^{2}\,\mathrm{d}^{3}r\,\mathrm{d}^{3}r' \le \lVert U\rVert_{\infty}^{2}\,\tfrac{4}{3}\pi a^{3}\cdot 8\pi a < \infty$ . So $\mathcal{T}$ is compact there and does have a discrete set of characteristic values, which is what makes the reading below legitimate; what is lost is the uniform bound of our Neumann theorem, so the convergence must be justified case by case. What does carry over is the mechanism that Arfken, Weber and Harris describe on p. 1067: the series converges while the parameter stays inside the smallest characteristic value of the kernel. Here the parameter is fixed at $\lambda = 1$ and the kernel moves with the energy instead, so the statement becomes one about $k$ : the wavenumbers at which the homogeneous equation $\chi = \mathcal{T}\chi$ has a nonzero solution are exactly the bound states and resonances of $U$ , and these are what bound the convergence. The one-dimensional example below exhibits this with everything in closed form: the series diverges precisely at energies below the binding energy of the potential's single bound state.
The first Born approximation is not unitary. By (B $'$ ) a real central potential gives a real $f_{\mathrm{B}}$ . A real amplitude cannot satisfy the exact conservation identities, which couple the imaginary part of the forward amplitude to the total scattered intensity; the one-dimensional analogue is Strauss's $\lvert R\rvert^{2}+\lvert T\rvert^{2} = 1$ , and the second example below shows the first Born approximation overshooting it by a computable amount. The defect is of second order in $U$ , which is consistent: the approximation is first order.
When to expect accuracy. Arfken, Weber and Harris state the practical rule: the Born approximation is expected to be best for weak potentials and high incident energy. Both halves are visible in (B): the amplitude is first order in $U$ , so weak potentials are favoured; and at large $k$ the factor $e^{-i\mathbf{q}\cdot\mathbf{r}'}$ oscillates rapidly over the scatterer except in a narrow forward cone, so the scattered amplitude is small compared with the incident wave. A usable quantitative test is to measure the first neglected term against the incident wave at the scattering centre, $\lvert(\mathcal{T}\psi_{\mathrm{inc}})(\mathbf{0})\rvert\ll1$ ; the first example evaluates it.
The alternative. When the series is useless — strong potential, low energy — one solves exactly instead, by separating variables and matching partial waves. Strauss carries this out for a hard sphere and reaches a closed-form series in spherical Bessel and Hankel functions (Spherical Bessel functions). That route needs a geometry that separates; Born needs only that the potential be weak.

Method(A first Born cross-section)Original

Given a potential energy $V$ and a beam energy $E$ :

Reduce. Form $U = 2mV/\hbar^{2}$ and $k = \sqrt{2mE}/\hbar$ . Confirm that $U$ decays fast enough to be treated as compactly supported; if $U$ is Coulombic, screen it and unscreen at the end.
Momentum transfer. Set $q = 2k\sin(\theta/2)$ , noting that $q$ runs over $[0,2k]$ as $\theta$ runs over $[0,\pi]$ .
Transform. Compute the unnormalised transform $\mathcal{F}U(\mathbf{q}) = \int e^{-i\mathbf{q}\cdot\mathbf{r}}U(\mathbf{r})\,\mathrm{d}^{3}r$ of the note. For a central $U$ use the radial form of step 5 of the proof, $\mathcal{F}U(q) = \frac{4\pi}{q}\int_{0}^{\infty} r\,U(r)\sin(qr)\,\mathrm{d}r ,$ the one-dimensional sine transform of $rU(r)$ . The same radial reduction is carried out in Fourier transform in several dimensions, but for the symmetric transform, where it reads $\tilde U(q) = (2\pi)^{-3/2}(4\pi/q)\int_{0}^{\infty} rU(r)\sin(qr)\,\mathrm{d}r$ ; a transform taken from that unit, or from its table of examples, must therefore be multiplied by $(2\pi)^{3/2} = 15.7496\ldots$ before it enters step 4.
Amplitude. $f_{\mathrm{B}} = -\mathcal{F}U(q)/4\pi$ , and $\mathrm{d}\sigma/\mathrm{d}\Omega = \lvert f_{\mathrm{B}}\rvert^{2}$ .
Total cross-section. Change variable from $\theta$ to $q$ . Since $q^{2} = 2k^{2}(1-\cos\theta)$ gives $q\,\mathrm{d}q = k^{2}\sin\theta\,\mathrm{d}\theta$ , $\sigma = 2\pi\int_{0}^{\pi}\lvert f_{\mathrm{B}}\rvert^{2}\sin\theta\,\mathrm{d}\theta = \frac{2\pi}{k^{2}}\int_{0}^{2k}\lvert f_{\mathrm{B}}(q)\rvert^{2}\,q\,\mathrm{d}q .$ This form is almost always the easier integral, and it displays the high-energy behaviour at a glance.
Test the approximation. Evaluate $\varepsilon = \lvert(\mathcal{T}\psi_{\mathrm{inc}})(\mathbf{0})\rvert = \frac{1}{k}\left\lvert\int_{0}^{\infty} e^{ikr}\,U(r)\,\sin(kr)\,\mathrm{d}r\right\rvert \quad(\text{central } U),$ the size of the first correction to the wave at the scattering centre relative to the incident wave there. The radial form follows from the angular integral $\int e^{i\mathbf{k}_{0}\cdot\mathbf{r}'}\mathrm{d}\Omega' = 4\pi\sin(kr')/(kr')$ . Accept the answer only if $\varepsilon\ll1$ .

Example(A Gaussian potential)Original

Take the repulsive Gaussian

U(r) = U_{0}\,e^{-r^{2}/a^{2}}, \qquad U_{0}>0 .

It is not compactly supported, so the proof does not apply to it verbatim: step 2 bounds the remainder by $\lvert R\rvert\le2a^{2}/r$ using $\lvert\mathbf{r}'\rvert\le a$ exactly, and here $U$ is nowhere zero. The repair is to cut the potential at radius $Ca$ . The truncated potential is a scatterer in the sense of the definition, so every estimate of the proof applies to it with $a$ replaced by $Ca$ ; and the discarded tail is small, since $\lvert U\rvert\le U_{0}e^{-C^{2}}$ there and

\int_{\lvert\mathbf{r}'\rvert>Ca}\lvert U\rvert\,\mathrm{d}^{3}r' = 4\pi U_{0}a^{3}\int_{C}^{\infty}u^{2}e^{-u^{2}}\,\mathrm{d}u = O\bigl(U_{0}a^{3}C\,e^{-C^{2}}\bigr),

so cutting it changes the amplitude — and every other integral below — by $O(e^{-C^{2}})$ , uniformly in $q$ . Since $C$ may be taken as large as one pleases, the formulas computed below from the full Gaussian are exact.

The amplitude. Step 3 of the method needs $\int_{0}^{\infty}re^{-r^{2}/a^{2}}\sin(qr)\,\mathrm{d}r$ . The integrand is even in $r$ , so the integral is half the integral over $\mathbb{R}$ , which is the imaginary part of

\tfrac12\int_{-\infty}^{\infty} r\,e^{-r^{2}/a^{2}}e^{iqr}\,\mathrm{d}r = \tfrac12\,\frac{1}{i}\frac{\mathrm{d}}{\mathrm{d}q}\int_{-\infty}^{\infty} e^{-r^{2}/a^{2}}e^{iqr}\,\mathrm{d}r = \tfrac12\,\frac{1}{i}\frac{\mathrm{d}}{\mathrm{d}q}\Bigl(a\sqrt{\pi}\,e^{-q^{2}a^{2}/4}\Bigr) = \frac{i\sqrt{\pi}a^{3}q}{4}e^{-q^{2}a^{2}/4}.

Hence $\int_{0}^{\infty}re^{-r^{2}/a^{2}}\sin(qr)\,\mathrm{d}r = \tfrac14\sqrt{\pi}a^{3}q\,e^{-q^{2}a^{2}/4}$ and, by (B $'$ ),

f_{\mathrm{B}}(q) = -\frac{\sqrt{\pi}}{4}\,U_{0}a^{3}\,e^{-q^{2}a^{2}/4} .

The amplitude is a Gaussian in the momentum transfer, of width $2/a$ : a scatterer of size $a$ deflects through angles up to about $1/(ka)$ and little beyond. The sign is negative, as it must be for a repulsive potential in this convention.

The cross-sections. $\mathrm{d}\sigma/\mathrm{d}\Omega = \tfrac{1}{16}\pi U_{0}^{2}a^{6}e^{-q^{2}a^{2}/2}$ , and step 5 gives

\sigma = \frac{2\pi}{k^{2}}\cdot\frac{\pi U_{0}^{2}a^{6}}{16}\int_{0}^{2k} e^{-q^{2}a^{2}/2}q\,\mathrm{d}q = \frac{\pi^{2}U_{0}^{2}a^{4}}{8k^{2}}\Bigl(1-e^{-2k^{2}a^{2}}\Bigr).

Two limits are worth noting. At low energy, $ka\ll1$ , the bracket is $2k^{2}a^{2}+O(k^{4}a^{4})$ and $\sigma\to\tfrac14\pi^{2}U_{0}^{2}a^{6}$ , a constant: the scattering becomes isotropic and energy-independent, which is what $f_{\mathrm{B}}(q)\to f_{\mathrm{B}}(0)$ says. At high energy the bracket saturates at $1$ and $\sigma\sim\pi^{2}U_{0}^{2}a^{4}/8k^{2}$ , falling off as $1/E$ .

Numbers. Measure lengths in units of $a$ and take the dimensionless strength $U_{0}a^{2} = 1$ at $ka = 1$ . Then, in units of $a$ and $a^{2}$ ,

f_{\mathrm{B}}(0) = -\frac{\sqrt{\pi}}{4} = -0.4431135, \qquad \left.\frac{\mathrm{d}\sigma}{\mathrm{d}\Omega}\right|_{\theta=0} = \frac{\pi}{16} = 0.1963495 ,

and at $\theta = \pi$ , where $q = 2k$ , so that the amplitude's exponent $-q^{2}a^{2}/4$ is $-1$ and the cross-section's is $-2$ ,

\left.\frac{\mathrm{d}\sigma}{\mathrm{d}\Omega}\right|_{\theta=\pi} = \frac{\pi}{16}e^{-2} = 0.02657302 ,

so forward scattering beats backward scattering by the factor $e^{2} = 7.389$ . The total cross-section is

\sigma = \frac{\pi^{2}}{8}\bigl(1-e^{-2}\bigr) = 1.2337006\times0.8646647 = 1.0667373 \;a^{2} ,

the number to check the whole calculation against. It is about a third of the geometric cross-section $\pi a^{2} = 3.1415927\,a^{2}$ of a disc of radius $a$ , which is the right order of magnitude for a soft scatterer of that size.

Is it any good? Step 6 gives

\varepsilon = \frac{1}{k}\left\lvert\int_{0}^{\infty} e^{ikr}U_{0}e^{-r^{2}/a^{2}}\sin(kr)\,\mathrm{d}r\right\rvert = \left\lvert\int_{0}^{\infty} e^{ir-r^{2}}\sin r\,\mathrm{d}r\right\rvert = \lvert 0.2690398+0.2801011\,i\rvert = 0.3884 ,

in the same units. That is not small: the wave at the centre of the scatterer differs from the incident wave by nearly $40\%$ , and the computed cross-section should be read as an order of magnitude, not a prediction. Weakening the potential to $U_{0}a^{2} = 0.2$ at the same energy scales $\varepsilon$ to $0.0777$ and $\sigma$ to $0.04\times1.0667373 = 0.04266949\,a^{2}$ , and there the first Born approximation is worth quoting.

Example(A delta potential in one dimension: summing the whole Born series)Original

In one dimension everything above can be carried out exactly, which makes the example the honest test of the method. Two warnings first: the potential below is one-dimensional and it is a measure, not a bounded function, so it is not a scatterer in the sense of the definition and neither the definition nor the statement applies to it. Nothing is inherited from them here. The Lippmann-Schwinger equation, the Born series and the amplitudes are all rederived from scratch in closed form, and the one general fact invoked at the end — that the series converges up to the smallest characteristic value — is verified directly against the exact answer rather than quoted. The equation is $\psi''+k^{2}\psi = U\psi$ with

U(x) = u_{0}\,\delta(x),

$u_{0}$ real with the dimensions of an inverse length; in the Schrodinger reading, $u_{0} = 2m\alpha/\hbar^{2}$ for a potential energy $\alpha\delta(x)$ . The outgoing Green's function of $\mathrm{d}^{2}/\mathrm{d}x^{2}+k^{2}$ is $G_{1}(x) = -\tfrac{i}{2k}e^{ik\lvert x\rvert}$ (Green's functions for the Helmholtz equation), the incident wave is $e^{ikx}$ travelling to the right, and (LS) reads

\psi(x) = e^{ikx} - \frac{iu_{0}}{2k}\,e^{ik\lvert x\rvert}\,\psi(0).

The delta has collapsed the integral entirely: the only unknown left is the single number $\psi(0)$ .

Exact solution. Setting $x = 0$ ,

\psi(0) = 1 - \frac{iu_{0}}{2k}\psi(0) \quad\Longrightarrow\quad \psi(0) = \frac{1}{1+iu_{0}/2k} = \frac{2k}{2k+iu_{0}} .

Reading off the two ends — with the beam coming in from the left, so that $e^{ikx}$ is the incident wave and the reflection and transmission coefficients $R$ and $T$ are defined by the far-field behaviour exactly as in Strauss's §13.4 (his beam arrives from the other side, which only swaps the roles of the two ends) —

\psi(x)\to T e^{ikx}\ (x\to+\infty), \qquad \psi(x)\to e^{ikx}+Re^{-ikx}\ (x\to-\infty),

T = \frac{2k}{2k+iu_{0}}, \qquad R = -\frac{iu_{0}}{2k}\,\psi(0) = \frac{-iu_{0}}{2k+iu_{0}} .

These satisfy $\lvert R\rvert^{2}+\lvert T\rvert^{2} = (u_{0}^{2}+4k^{2})/(4k^{2}+u_{0}^{2}) = 1$ , which is Strauss's Theorem 1 — a check on the exact solution, not an assumption in it.

The Born series. Iterating instead of solving gives $\psi(0) = \sum_{n\ge0}(-iu_{0}/2k)^{n}$ , a geometric series with ratio $-iu_{0}/2k$ . So the Born series converges if and only if

\frac{\lvert u_{0}\rvert}{2k} < 1,

and when it does it sums to the exact $\psi(0)$ above. The $n$ -th Born approximation keeps $n$ terms.

Where the divergence comes from. For an attractive delta, $u_{0}<0$ , the transmission coefficient $T(k) = 2k/(2k+iu_{0})$ has a pole at $k = i\lvert u_{0}\rvert/2$ , that is at energy $k^{2} = -u_{0}^{2}/4$ : the single bound state of the delta well, whose wavefunction $e^{-\lvert u_{0}\rvert\lvert x\rvert/2}$ one checks directly from the jump condition $\psi'(0^{+})-\psi'(0^{-}) = u_{0}\psi(0)$ . The Born series converges exactly when $k>\lvert u_{0}\rvert/2$ , that is when the incident energy $k^{2}$ exceeds the binding energy $u_{0}^{2}/4$ . This is the general rule of Neumann series, made concrete: the radius of convergence is set by the characteristic value of the kernel, and the characteristic value is the bound state.

Numbers. Take $u_{0} = 1$ and $k = 2$ , so that the ratio is $u_{0}/2k = 1/4$ and the series converges comfortably. The exact coefficients are

T = \frac{4}{4+i} = \frac{16-4i}{17}, \qquad R = \frac{-i}{4+i} = \frac{-1-4i}{17},

\lvert T\rvert^{2} = \frac{16}{17} = 0.9411765, \qquad \lvert R\rvert^{2} = \frac{1}{17} = 0.0588235, \qquad \lvert R\rvert^{2}+\lvert T\rvert^{2} = 1 .

The first Born approximation replaces $\psi(0)$ by $1$ :

R_{1} = -\frac{i}{4}, \qquad T_{1} = 1-\frac{i}{4}, \qquad \lvert R_{1}\rvert^{2} = \frac{1}{16} = 0.0625, \qquad \lvert T_{1}\rvert^{2} = \frac{17}{16} = 1.0625 .

Three things are visible at once. The reflected intensity is too large by $1/16-1/17 = 1/272$ , a relative error of $6.25\%$ . That is the square of the series ratio $u_{0}/2k = \tfrac14$ , and exactly so: $R = R_{1}/(1+iu_{0}/2k)$ , hence $\lvert R\rvert^{2} = \lvert R_{1}\rvert^{2}/(1+u_{0}^{2}/4k^{2})$ , and the intensity is quadratic in an amplitude whose next correction is one ratio down. The transmitted intensity exceeds $1$ , which is impossible. And the total

\lvert R_{1}\rvert^{2}+\lvert T_{1}\rvert^{2} = \frac{18}{16} = 1.125

breaks Strauss's identity by $12.5\%$ , which is $2(u_{0}/2k)^{2}$ : the unitarity defect is second order in the potential, exactly the order the first Born approximation drops.

Going to second order, $\psi(0)\approx1-i/4$ , gives $R_{2} = -\tfrac{1}{16}-\tfrac{i}{4}$ and $T_{2} = \tfrac{15}{16}-\tfrac{i}{4}$ , with

\lvert R_{2}\rvert^{2}+\lvert T_{2}\rvert^{2} = \frac{17}{256}+\frac{241}{256} = \frac{258}{256} = 1.0078125 ,

a defect of $0.78\%$ instead of $12.5\%$ . The series repairs unitarity order by order, and no truncation of it ever has unitarity exactly.

ApplicationOriginal

Diffraction measures a Fourier transform. Formula (B) is the reason that scattering is the standard tool for looking at structure. In the first Born approximation the measured intensity is

\frac{\mathrm{d}\sigma}{\mathrm{d}\Omega} = \frac{1}{16\pi^{2}}\bigl\lvert\mathcal{F}U(\mathbf{q})\bigr\rvert^{2},

the squared modulus of the Fourier transform of the scatterer at the momentum transfer $\mathbf{q}$ . Rotating the sample and scanning the angle scans $\mathbf{q}$ over a ball of radius $2k$ , so a diffraction experiment is a measurement of $\mathcal{F}U$ on that ball. Three consequences are generic:

Resolution is set by the wavelength. Since $q\le2k$ , the finest spatial period recorded is the $\pi/k$ of the definition block — half a wavelength. This is why structural work uses X-rays and neutrons rather than light.
Only the modulus is measured. The phase of $\mathcal{F}U(\mathbf{q})$ is lost, and recovering it is the phase problem of crystallography. For a central real potential $f_{\mathrm{B}}$ is real by (B $'$ ), so the lost information is only a sign — one reason that spherically averaged scattering is easier to invert than single-crystal diffraction.
Weak scattering is a feature. The approximation is linear in $U$ , so the inverse problem — infer $U$ from the data — is linear too, and can be solved by an inverse Fourier transform. Multiple scattering, that is the later terms of the Born series, makes the inverse problem nonlinear. The same linearisation under the same name is used in diffraction tomography and in seismic imaging, where it plays the role that the projection-slice theorem plays for straight-ray tomography (The Radon transform).

Rutherford's cross-section as a screened limit. Take the screened Coulomb (Yukawa) potential

V(r) = \frac{Z_{1}Z_{2}e^{2}}{4\pi\epsilon_{0}}\frac{e^{-\mu r}}{r}, \qquad\text{so}\qquad U(r) = \beta\frac{e^{-\mu r}}{r}, \quad \beta = \frac{2m}{\hbar^{2}}\frac{Z_{1}Z_{2}e^{2}}{4\pi\epsilon_{0}} .

Then (B $'$ ) needs only $\int_{0}^{\infty}e^{-\mu r}\sin(qr)\,\mathrm{d}r = q/(\mu^{2}+q^{2})$ , giving

f_{\mathrm{B}}(q) = -\frac{\beta}{\mu^{2}+q^{2}} .

Equivalently $\mathcal{F}\bigl(e^{-\mu r}/r\bigr)(q) = 4\pi/(\mu^{2}+q^{2})$ . The same Yukawa transform is listed in Fourier transform in several dimensions in that unit's symmetric normalisation, as $(2\pi)^{-3/2}\cdot4\pi/(q^{2}+\mu^{2})$ ; it is our value divided by $(2\pi)^{3/2}$ , as the note says it must be. Letting the screening length $1/\mu$ grow without bound and using $q = 2k\sin(\theta/2)$ together with $\hbar^{2}k^{2}/2m = E$ ,

\frac{\mathrm{d}\sigma}{\mathrm{d}\Omega} = \frac{\beta^{2}}{q^{4}} = \left(\frac{\beta}{4k^{2}}\right)^{2}\frac{1}{\sin^{4}(\theta/2)} = \left(\frac{Z_{1}Z_{2}e^{2}}{16\pi\epsilon_{0}E}\right)^{2}\frac{1}{\sin^{4}(\theta/2)} ,

which is Rutherford's formula, including the constant. The result is famously better than it has any right to be: the remark above explains that the unscreened Coulomb potential does not satisfy the hypotheses under which $f$ was defined, and in fact the exact Coulomb amplitude differs from $f_{\mathrm{B}}$ by a phase — which the cross-section then discards.

Beyond quantum mechanics. Nothing in the derivation used quantum mechanics; it used $(\nabla^{2}+k^{2})\psi = U\psi$ . The same equation governs a time-harmonic acoustic field in a medium whose sound speed varies, with $U = k^{2}(1-c_{0}^{2}/c(\mathbf{r})^{2})$ , and a scalar electromagnetic field in a medium of varying refractive index, with $U = -k^{2}(n(\mathbf{r})^{2}-1)$ . In both cases the first Born approximation is the statement that the field scattered by a weak inhomogeneity is the Fourier transform of the contrast, and the convergence condition of the Born series is the statement that multiple scattering inside the inhomogeneity may be neglected. Ultrasound imaging, optical coherence tomography and the interpretation of radar returns from weak targets all rest on it; when the contrast is large — a metal object, a gas bubble — the series fails and one is back to solving (LS), typically by discretising it into a linear system and applying the iterative methods of Neumann series.