Shows that a harmonic function equals its average over every sphere and every ball lying in its domain, that the converse holds for C2 functions, that consequently it attains its maximum and minimum only on the boundary unless it is constant, and that the Dirichlet problem for Poisson's equation therefore has at most one solution, which depends continuously on its data.
Let U⊆Rn be open. A function u:U→R is harmonic in U if u∈C2(U) and
∇2u=i=1∑n∂xi2∂2u=0at every point of U.
The equation ∇2u=0 is Laplace's equation; the inhomogeneous equation ∇2u=f with f given is Poisson's equation. Ahlfors, working in the plane, builds the continuity of u and of its first and second partial derivatives into the definition and remarks that the regularity demanded can be relaxed; Strauss calls a solution of Laplace's equation a harmonic function without naming a smoothness class.
Two pieces of notation are used throughout. B(x,r) is the open ball of radius r about x and ∂B(x,r) its bounding sphere. If α(n) denotes the volume of the unit ball of Rn, then
∣B(x,r)∣=α(n)rn,∣∂B(x,r)∣=nα(n)rn−1,
so that α(2)=π, α(3)=34π, and the two sphere areas are 2πr and 4πr2. We write ∂u/∂n=n^⋅∇u for the outward normal derivative, as in Green's identities.
Signs. Evans writes Poisson's equation as −Δu=f, so that the operator −Δ is the non-negative one; Strauss and Ahlfors write Δu=f and Δu=0. We use ∇2 for the Laplacian and put the sign where the source we are following puts it, saying each time which it is.
Theorem(Mean value property)AdaptedAdapted from Evans (2010), §2.2.2, Thm 2, p. 25; Strauss (2008), §7.1, Eq. (6), p. 181
A harmonic function carries no information at the centre of a ball that is not already carried, on average, by its values on the surrounding sphere. Precisely: let U⊆Rn be open, let u be harmonic in U, and let x and r>0 be such that the ball B(x,r) and its boundary sphere lie in U. Then
In particular, the right-hand sides do not depend on r. Evans states the two equalities together in Rn for u∈C2(U) and "each ball B(x,r)⊆U"; we read that as including the boundary sphere, which the surface integral needs. Strauss derives the spherical form in three dimensions only, and remarks that the same idea works in n dimensions.
Theorem(Converse)AdaptedAdapted from Evans (2010), §2.2.2, Thm 3, p. 26
Let U⊆Rn be open and let u∈C2(U) satisfy the first equality of (MV) for every ball B(x,r)⊆U. Then u is harmonic in U.
Idea of the proofOriginal
Write ϕ(r) for the average of u over the sphere ∂B(x,r) and watch it as r grows. Rescaling the integral to the fixed unit sphere lets us differentiate in r under the integral sign, and the derivative that comes out is exactly the average of the normal derivative ∂u/∂n over the sphere. By Green's first identity with v≡1, that integral equals the integral of ∇2u over the ball, which is zero. So ϕ is constant in r; and as r→0+ the average over a shrinking sphere tends to u(x) by continuity. That gives the sphere form of (MV), and the ball form follows by integrating the sphere form in r.
The converse is the same computation read backwards. The identity ϕ′(r)=nr×(average of ∇2u over B(x,r)) holds for every C2 function, harmonic or not. If the averages are constant in r then ϕ′≡0, so the average of ∇2u over every ball vanishes; a continuous function whose average over every ball vanishes is identically zero.
Step 1 (the spherical average, rescaled). Fix x with B(x,R)⊆U and for 0<r<R put
ϕ(r)=∣∂B(x,r)∣1∫∂B(x,r)u(y)dS(y).
Substituting y=x+rz maps ∂B(x,r) onto the unit sphere ∂B(0,1) and scales the surface element by rn−1, so, since ∣∂B(x,r)∣=nα(n)rn−1,
ϕ(r)=nα(n)1∫∂B(0,1)u(x+rz)dS(z).(1)
The domain of integration in (1) no longer depends on r.
Step 2 (differentiate). The integrand of (1) is C1 in (r,z) on [0,R)×∂B(0,1), because u∈C2, and the unit sphere is compact; so we may differentiate under the integral sign:
ϕ′(r)=nα(n)1∫∂B(0,1)∇u(x+rz)⋅zdS(z).
Undoing the substitution, z=(y−x)/r is the outward unit normal to ∂B(x,r) at y, and dS(z)=r−(n−1)dS(y), so
The derivative of the spherical average is the average of the normal derivative.
Step 3 (use Green's identity). The ball B(x,r) has a smooth boundary and u∈C2 on its closure, so part 2 of Green's identities — the case v≡1 — applies. That unit displays its identities for a bounded region of R3; what we need here is the n-dimensional version, which is the form recorded in the closing paragraph of the same statement block (Evans's hypotheses: U⊂Rn bounded and open with C1 boundary, u∈C2(U)), and the ball B(x,r) with B(x,r)⊆U meets them. It turns (2) into
using ∣B(x,r)∣=α(n)rn in the second form. Note that (3) has used nothing about u beyond u∈C2.
Step 4 (harmonicity). Now let u be harmonic. Then the integrand in (3) vanishes, so ϕ′(r)=0 for all 0<r<R and ϕ is constant there.
Step 5 (the limit at the centre). Let ε>0. Since u is continuous at x, there is δ>0 with ∣u(y)−u(x)∣<ε for ∣y−x∣<δ. For r<δ, averaging this inequality over ∂B(x,r) gives ∣ϕ(r)−u(x)∣≤ε, because an average of u(x) is u(x). Hence ϕ(r)→u(x) as r→0+, and by Step 4 the constant value of ϕ is u(x). This is the first equality of (MV).
Step 6 (from spheres to the ball). Integrating in polar coordinates about x and inserting the sphere result,
Step 7 (the converse). Let u∈C2(U) satisfy the sphere equality for every ball in U. Then ϕ is constant on (0,R) for every admissible R, so ϕ′≡0, and (3) — which needed only u∈C2 — says that
∫B(x,r)∇2udV=0for every B(x,r)⊆U.
Suppose ∇2u(x0)=0 for some x0∈U; replacing u by −u if necessary, assume ∇2u(x0)>0. Since ∇2u is continuous, it stays above 21∇2u(x0)>0 on some ball B(x0,ρ)⊆U, whose integral is then strictly positive. That contradicts the display. Hence ∇2u≡0 on U. □
RemarkAdaptedAdapted from Ahlfors (1979), §4.6.1, Thm 19 and Eq. (60), p. 164, and §4.6.2, pp. 164–165
Ahlfors reaches the plane case of (MV) by a route with no divergence theorem in it. Write ∗du=−uydx+uxdy for the conjugate differential. His Theorem 19 (Eq. (60), p. 164) is that ∫γ(u1∗du2−u2∗du1)=0 whenever u1 and u2 are harmonic in a region Ω and the cycle γ is homologous to zero in Ω. He proves it by reducing the claim to the boundary of a rectangle, where single-valued conjugate functions exist and the integrand is the imaginary part of an analytic differential, so Cauchy's theorem applies.
Two specialisations give the plane mean value property for a u harmonic on the disc ∣z∣<ρ. On a circle ∣z∣=r one has ∗du=r(∂u/∂r)dθ.
With Ω the punctured disc 0<∣z∣<ρ, u1=logr, u2=u, and γ the difference of two concentric circles — a cycle homologous to zero there — Theorem 19 gives the same value on both circles for
∫∣z∣=rudθ−logr∫∣z∣=rr∂r∂udθ,
so this combination is independent of r. Harmonicity in the annulus between the two circles is all that is used.
With u1=1 and u2=u, Theorem 19 collapses to his earlier Eq. (58), ∫γ∗du=0. On the punctured disc that only repeats that ∫∣z∣=rr(∂u/∂r)dθ is independent of r; a single circle is not homologous to zero there. To get the stronger conclusion that this integral vanishes, one needs u harmonic on the whole disc, for then the circle ∣z∣=r is by itself homologous to zero and (58) applies to it alone.
Subtracting the second from the first leaves ∫∣z∣=rudθ independent of r, which is the mean value property; in the annulus, where only the first bullet is available, Ahlfors gets instead a circular mean of the form alogr+β (his Theorem 20, Eq. (61), p. 165). The conclusion is the same as Step 5's; the hypotheses differ in that Ahlfors's argument is about a region in the plane and is sensitive to the topology of γ, whereas the argument above is dimension-free.
Let U⊆Rn be open and bounded, and let u∈C2(U)∩C(U) be harmonic in U. Then:
Weak maximum principle.Umaxu=∂Umaxu,Uminu=∂Uminu.
Strong maximum principle. If in addition U is connected and u(x0)=maxUu for some x0∈U, then u is constant on U.
Part 2 applied to −u gives the same statement for an interior minimum. Equivalently: a nonconstant harmonic function on a connected bounded U attains its maximum and its minimum on ∂U and at no interior point.
The sources state this with different scaffolding. Evans assumes U open and bounded and proves part 2 first, deducing part 1 from it. Strauss states it for a connected bounded open set in two or three dimensions, with the interior conclusion included, and proves the boundary bound of part 1 by the device below. Stakgold and Holst give the weak version for a bounded domain Ω with boundary Γ and a continuous u, in the form: if −∇2u≤0 in Ω then maxu is on Γ, if −∇2u≥0 then minu is on Γ, and if ∇2u=0 then both are — so for them the harmonic case is the overlap of two one-sided principles for subsolutions and supersolutions. Ahlfors states it as Theorem 21 without a boundedness hypothesis, for a harmonic function on a plane region: if it is not constant it attains no maximum value and no minimum value at any point where it is harmonic, and the extrema over a closed bounded subset are therefore attained on that subset's boundary.
Idea of the proofOriginal
There are two proofs, and they are worth having both because they need different things.
The first proves only part 1, and needs no mean value property. At an interior maximum of a C2 function each second pure derivative is ≤0, so the Laplacian is ≤0 there. That nearly contradicts ∇2u=0, but not quite, since ≤0 allows =0. The repair is to tilt u by a small multiple of ∣x∣2, whose Laplacian is the positive constant 2n: the tilted function has strictly positive Laplacian, so it has no interior maximum at all, and letting the tilt go to zero transfers the boundary bound to u.
The second proves part 2 from the ball form of (MV). If u reaches its global maximum M at an interior point, then M is also the average of u over a small ball around that point; but u≤M there, and a function ≤M whose average is M must equal M throughout. So the set where u=M swallows a whole ball around each of its points: it is open. It is also closed in U, being a level set of a continuous function, and connectedness then forces it to be all of U.
Step 1 (no interior maximum for a strict supersolution). Let w∈C2(U) satisfy ∇2w>0 at every point of U. If w had a local maximum at x1∈U, then each function t↦w(x1+tei) would have a local maximum at t=0, so ∂2w/∂xi2(x1)≤0 by the second-derivative test, for each i. Summing over i gives ∇2w(x1)≤0, contradicting the hypothesis. So w has no local maximum in U.
Step 2 (the tilted function). For ε>0 set vε(x)=u(x)+ε∣x∣2. Since ∇2∣x∣2=2n and u is harmonic,
∇2vε=0+2nε>0in U,
so Step 1 applies to vε: it has no local maximum in the open set U.
Step 3 (where the maximum is). Here boundedness of U and continuity of u on U are used: U is compact and vε is continuous on it, so vε attains a maximum on U. By Step 2 that maximum is not attained in U, hence it is attained at some point of ∂U.
Step 4 (remove the tilt). Let L=maxy∈∂U∣y∣, finite because U is bounded. For every x∈U,
the first inequality because ε∣x∣2≥0. The left and right ends do not depend on ε except through the last term, so letting ε→0+ gives u(x)≤max∂Uu for all x∈U. Since ∂U⊆U, the reverse inequality is automatic, and the two maxima are equal.
Step 5 (the minimum). Apply Steps 1–4 to −u, which is harmonic as well, and use min∂Uu=−max∂U(−u). □
Proof(Part 2, from the ball mean value property)AdaptedAdapted from Evans (2010), §2.2.3a, p. 26; Stakgold & Holst (2011), §8.3, pp. 495–496
Step 1 (set-up). Let M=maxUu, which exists because U is compact and u∈C(U), and let
E={x∈U:u(x)=M},
assumed nonempty: x0∈E.
Step 2 (E is open). Let x∈E and let r be any radius with 0<r<dist(x,∂U), so that B(x,r)⊆U. The ball form of (MV) gives
M=u(x)=∣B(x,r)∣1∫B(x,r)udV.
Hence ∫B(x,r)(M−u)dV=0. The integrand is continuous and non-negative, since M is the maximum of u over all of U. A continuous non-negative function with vanishing integral vanishes identically: if M−u(y)>0 at some y∈B(x,r), it would exceed half that value on a small ball about y and the integral would be positive. So u≡M on B(x,r), and therefore B(x,r)⊆E. Thus E is open.
Step 3 (E is closed in U).E=u−1({M})∩U and u is continuous, so E is the intersection of U with a closed set: it is relatively closed in U.
Step 4 (connectedness).U is connected, and E is a nonempty subset of U that is both open and relatively closed. The only such subset is U itself, for otherwise U=E∪(U∖E) would split U into two disjoint nonempty relatively open pieces. Hence E=U and u≡M on U. □
RemarkOriginal
Where the hypotheses are used, and what fails without them.
Boundedness of U, in Steps 3 and 4 of the first proof. It is what makes U compact and L finite. Both conclusions genuinely fail on an unbounded domain: on the half-plane U={(x,y):y>0} the function u=y is harmonic and vanishes on ∂U, so max∂Uu=0 while u is unbounded above on U. Ahlfors's Theorem 21 avoids the issue by asserting only the local statement — no interior extremum — and then drawing the boundary conclusion on a closed bounded set.
Continuity on U. Without it there is nothing to maximize on the boundary. On the punctured disc U={0<∣x∣<1}⊂R2 the function u=log(1/∣x∣) is harmonic in U and zero on the outer boundary circle, but it has no continuous extension to the boundary point 0, and its supremum over U is +∞.
Connectedness, in Step 4 of the second proof only. If U is the union of two disjoint open balls, u=0 on one and u=1 on the other is harmonic on U, continuous on U, and attains its maximum at every interior point of the second ball without being constant on U. Part 1 survives: each component has its own boundary. This is why Evans attaches connectedness to part 2 and not to part 1.
u∈C2, in Step 1 of the first proof. The second-derivative test needs the pure second derivatives to exist. (By Evans's smoothness theorem, a continuous function with the mean value property is automatically C∞ inside U, so for harmonic functions this costs nothing in the interior; see Poisson's integral formula.)
Which proof proves what. The first proof never uses the mean value property and gives only part 1; it is the one Strauss and Stakgold–Holst give, and it extends at once to subsolutions, ∇2u≥0, since Step 2 only needs ∇2vε>0. The second proof uses (MV) and gives the sharper part 2 but says nothing about subsolutions directly. Strauss states the interior conclusion with the weak proof and defers the interior argument to his §6.3; Stakgold and Holst likewise prove the weak version in their §1.3 and the strong version only in §8.3, after the mean value theorem is available.
Theorem(Uniqueness, comparison and stability for the Dirichlet problem)AdaptedAdapted from Evans (2010), §2.2.3a, Thm 5, p. 27; Strauss (2008), §6.1, pp. 155–156; Stakgold & Holst (2011), §1.3, pp. 73–74
Let Ω⊂Rn be a bounded domain with boundary Γ, let f∈C(Ω) and h∈C(Γ), and consider
−∇2u=f in Ω,u=h on Γ.(D)
Solutions are required to lie in C2(Ω)∩C(Ω).
Uniqueness. Problem (D) has at most one solution.
Comparison. If u1 solves (D) with data (f1,h1), u2 solves it with data (f2,h2), and f1≥f2 in Ω together with h1≥h2 on Γ, then u1≥u2 in Ω.
Continuous dependence on the data. Let v solve −∇2v=1 in Ω, v=0 on Γ. Then every solution of (D) obeys
∥u∥∞≤∥v∥∞∥f∥∞+∥h∥∞,Γ,
where ∥⋅∥∞ is the supremum over Ω and ∥h∥∞,Γ=supΓ∣h∣. In particular the solution map is continuous from the supremum norm on the data to the supremum norm on the solution.
Stakgold and Holst prove 2 and 3 from the weak maximum principle alone, state 3 first for h=0 (with the general case as their Exercise 1.3.2(a)), and observe that no constant smaller than ∥v∥∞ works in 3, because taking f≡1 makes it an equality. Evans states 1 for f∈C(U), g∈C(∂U) and solutions in C2(U)∩C(U); Strauss states and proves 1 for the form ∇2u=f.
Idea of the proofOriginal
All three parts are the weak maximum principle applied to a difference. Subtracting two solutions of (D) kills the source and the boundary data, leaving a harmonic function with zero boundary values; the maximum principle squeezes it between 0 and 0. For comparison, the difference is not harmonic but a subsolution, −∇2w≤0, with w≥0 on Γ, and the one-sided principle still puts its maximum on the boundary. For the bound, the function v is used as a barrier: ∥f∥∞v± (the boundary bound) dominates u by comparison, and evaluating the barrier gives the estimate.
Step 1 (uniqueness). Let u1,u2 both solve (D) and put w=u1−u2. By linearity −∇2w=f−f=0 in Ω and w=h−h=0 on Γ, and w∈C2(Ω)∩C(Ω). Part 1 of the maximum principle, applied to the bounded domain Ω, gives
Ωmaxw=Γmaxw=0,Ωminw=Γminw=0,
so 0≤w≤0: w≡0 and u1=u2.
Step 2 (a one-sided principle for subsolutions). Let w∈C2(Ω)∩C(Ω) satisfy −∇2w≤0 in Ω, that is ∇2w≥0. The proof of part 1 of the maximum principle goes through unchanged: vε=w+ε∣x∣2 has ∇2vε≥2nε>0, so Steps 1–4 there give maxΩw=maxΓw. (Only the upper bound survives; a subsolution may well take its minimum inside.)
Step 3 (comparison). Put w=u1−u2. Then −∇2w=f1−f2≥0 in Ω. Careful with the sign: this says ∇2w≤0, so it is −w that is a subsolution in the sense of Step 2. Applying Step 2 to −w,
Ωmax(−w)=Γmax(−w)=Γmax(h2−h1)≤0,
because h1≥h2 on Γ. Hence −w≤0, i.e. u1≥u2 throughout Ω.
Step 4 (the barrier). Let u solve (D), set F=∥f∥∞ and H=∥h∥∞,Γ, and let v be as in part 3. Consider
w±=Fv±u+H.
Then −∇2w±=F±f in Ω, which is ≥0 since ∣f∣≤F; and on Γ, w±=0±h+H≥0 since ∣h∣≤H. So w± is a nonnegative-data supersolution, and comparison against the zero solution — Step 3 with (f2,h2)=(0,0), u2≡0 — gives w±≥0 in Ω. Reading the two signs,
−Fv−H≤u≤Fv+Hin Ω.
Step 5 (evaluate). Note that v≥0 by Step 3 (its data 1≥0 and 0≥0 dominate the zero data), so ∣Fv+H∣≤F∥v∥∞+H. Taking the supremum of Step 4 over Ω,
∥u∥∞≤∥v∥∞∥f∥∞+∥h∥∞,Γ,
which is part 3. If two data sets (f1,h1) and (f2,h2) are given, applying this to the difference of the solutions — which solves (D) with data (f1−f2,h1−h2) — bounds ∥u1−u2∥∞ by ∥v∥∞∥f1−f2∥∞+∥h1−h2∥∞,Γ: the solution depends continuously on the data. □
RemarkOriginal
Uniqueness has two independent proofs in this handbook, and they assume different things. The one above uses the maximum principle and needs only u∈C2(Ω)∩C(Ω), with no regularity of Γ and no derivatives of u on the boundary. The energy proof in Green's identities puts v=u in Green's first identity and needs u∈C2(Ω) and a piecewise C1 boundary, because it integrates by parts; in exchange it handles the Neumann and Robin conditions, which the maximum principle does not. Strauss gives both, the maximum-principle proof in his §6.1 and the energy proof in his §7.1; Stakgold and Holst give both as well, in their §1.3 and §8.3.
Example(The mean value property tested in two and three dimensions)Original
(a) A plane polynomial. The function
u(x,y)=x3−3xy2
is harmonic on R2: uxx=6x and uyy=−6x. (It is Re(x+iy)3, so harmonicity is also Strauss's third example.) Take the centre (x0,y0)=(1,2), where u=1−3⋅1⋅4=−11. Parametrize the circle of radius r about it by (1+rcosθ,2+rsinθ) and expand:
independent of r, as (MV) requires. The two r2 terms cancel because uxx+uyy=0: indeed for any C2 function the O(r2) term of the circular average is 41r2∇2u(x0,y0), which is the two-dimensional case of (3). Numerically the average is −11.000000 for r=0.5, 1 and 3, and the average over the solid disc of radius 2 is also −11.000000.
(b) Newton's theorem in three dimensions. Let a=(0,0,2) and
u(x)=∣x−a∣1,
the Coulomb potential of a unit point source at a. It is harmonic on R3∖{a} (see the fundamental solution), hence on a neighbourhood of the closed unit ball, which does not contain a. So (MV) applies with x=0, r=1 and predicts that the average of u over the unit sphere is u(0)=1/∣a∣=21. Checking directly, with θ the angle from the z-axis and dS=sinθdθdφ on the unit sphere,
using t=cosθ. This is Newton's theorem: a sphere of uniform surface density attracts an external point exactly as if its mass were concentrated at its centre, read here in the reciprocal direction. Note that the mean value property says nothing when the singularity is inside the sphere, and indeed for ∣a∣<1 the same integral gives 1, not 1/∣a∣.
Example(The stability constant for the unit disc is sharp)Original
Take Ω to be the unit disc in R2 and compute the barrier v of part 3 explicitly. We need −∇2v=1 with v=0 on r=1. Looking for a radial solution and using ∇2v=v′′+v′/r for radial functions in the plane,
v(r)=41−r2⟹v′′+rv′=−21−21=−1,
and v(1)=0. So ∥v∥∞=v(0)=41, and part 3 reads
∥u∥∞≤41∥f∥∞+∥h∥∞,Γ
for every solution of −∇2u=f in the unit disc with u=h on its boundary.
The constant 41 cannot be reduced. Take f≡8 and h≡0. Then u=2(1−r2) solves the problem, since ∇2u=−8 and u(1)=0, and ∥u∥∞=u(0)=2, while the bound gives 41⋅8+0=2 exactly. This is the sharpness remark of Stakgold and Holst, which they make for f≡1.
A case with both terms active. Take f≡8 and h=3cosθ on the boundary circle, i.e. h=3x. Since 3x is harmonic and equals 3cosθ on r=1, the solution is
u=2(1−r2)+3x,
whose supremum over the closed disc is attained where 2(1−x2−y2)+3x is largest: setting the gradient to zero gives y=0 and −4x+3=0, so x=43 and ∥u∥∞=2(1−169)+49=87+49=825=3.125. The bound gives 41⋅8+3=5. The estimate holds with room to spare, as it must: the two extreme cases that make each term sharp cannot be arranged at the same place.
ApplicationOriginal
Electrostatics and the shielding of fields. Strauss's first example is Maxwell's equations in a charge-free region: ∇×E=0 gives E=−∇ϕ, and ∇⋅E=4πρ (Gaussian units) then gives ∇2ϕ=−4πρ, so the potential is harmonic wherever the charge density vanishes. The maximum principle is then a physical statement with teeth: inside a charge-free cavity the potential can have no interior maximum or minimum, so a charged particle cannot sit in stable electrostatic equilibrium in free space. That is Earnshaw's theorem, and it is why electrostatic traps for ions must be time-dependent. The same argument shows that a conductor held at a single potential encloses a cavity at exactly that potential — a Faraday cage — since the boundary data are constant and uniqueness forces the constant solution.
Steady flow and steady temperature. For an incompressible irrotational flow, Strauss's second example, the velocity potential is harmonic, and the maximum principle says the extremes of the potential sit on the boundary. For steady heat conduction with no sources the temperature is harmonic, and the weak maximum principle says that the hottest and coldest points of a body in steady state are on its surface: a result that is obvious physically and would be hard to prove without it. The comparison principle is the formal version of the statement that raising the temperature anywhere on the boundary, or adding heat anywhere inside, raises the temperature everywhere.
Error control for numerical solvers. Part 3 of the uniqueness statement is what makes a discretized Poisson solver trustworthy. Suitable discretizations of −∇2 satisfy a discrete maximum principle — the five-point finite-difference stencil on a uniform grid, and piecewise-linear finite elements on a triangulation with no obtuse angle, the condition that keeps the off-diagonal stiffness entries non-positive — and the discrete analogue of the barrier estimate — with ∥v∥∞=41 on the unit disc, as computed above — converts a bound on the residual f into a bound on the error in the computed solution, in the supremum norm and with no reference to derivatives of the exact solution. The mean value property is the statement that the exact solution is the fixed point of the averaging operator, which is precisely what the five-point stencil on a uniform grid discretizes, each value being replaced by the average of its four neighbours; Strauss's fourth example makes the probabilistic version of this explicit, with u the probability that a Brownian particle started at a point leaves the region through a chosen piece of its boundary, and that interpretation is the basis of Monte Carlo (walk-on-spheres) solvers for the Dirichlet problem.
Harmonic functions as the real parts of analytic functions. In the plane, Strauss's and Ahlfors's third setting, an analytic function splits into a harmonic real part and a harmonic imaginary part, so the maximum principle for harmonic functions and the maximum modulus principle for analytic functions are the same theorem seen twice. Ahlfors notes the direction that is less often used: the maximum principle for an analytic f with no zeros follows from the harmonic one applied to log∣f∣.
Finds all radial harmonic functions, normalizes them into the fundamental solution Φ with −∇2Φ=δ0 — the logarithm in the plane and 1/(4πr) in space — and proves that for f∈Cc2(Rn) the Newtonian potential u=Φ∗f is C2 and solves Poisson's equation −∇2u=f in all of Rn.
Theorem(The radial harmonic functions)AdaptedAdapted from Evans (2010), §2.2.1a, Eq. (5), p. 21; Strauss (2008), §6.1, pp. 158, 160
Laplace's equation is unchanged by rotations about the origin, so it is natural to look for solutions that depend only on r=∣x∣. Let v be twice differentiable on (0,∞) and put u(x)=v(∣x∣) on Rn∖{0}. Then
∇2u=v′′(r)+rn−1v′(r),(1)
so u is harmonic away from the origin if and only if
Strauss obtains the two cases n=2 and n=3 separately from the polar and spherical forms of the Laplacian, writing (2) as (rur)r=0 and (r2ur)r=0 and integrating; Evans does the general n at once, as above. For n=1 the equation reads v′′=0, and the harmonic functions are the affine ones; no singular solution appears, which is why the one-dimensional case gives no hint of what follows.
Idea of the proofOriginal
Writing out ∇2 for a radial function is a chain-rule computation, and the only thing worth noticing is where the n−1 comes from: each of the n coordinate directions contributes a term v′′xi2/r2 plus a term v′(1/r−xi2/r3), and summing xi2=r2 collapses the first group to v′′ and the second to (n−1)v′/r. Equation (2) is then first order in v′: it says that rn−1v′ is constant, which is the statement that the flux of ∇u through the sphere of radius r does not depend on r — exactly what Green's identity with v≡1 predicts for a harmonic function. Integrating once more gives (3), and the split at n=2 is the familiar split between ∫r−1dr and ∫r−pdr.
ProofAdaptedAdapted from Evans (2010), §2.2.1a, p. 21
Step 1 (derivatives of r). For x=0 and each i, differentiating r=(x12+⋯+xn2)1/2 gives ∂r/∂xi=xi/r.
Step 2 (the Laplacian of a radial function). By the chain rule, uxi=v′(r)xi/r, and differentiating again, using the product rule and Step 1,
uxixi=v′′(r)r2xi2+v′(r)(r1−r3xi2).
Summing over i=1,…,n and using ∑ixi2=r2,
∇2u=v′′(r)r2r2+v′(r)(rn−r3r2)=v′′+rn−1v′,
which is (1). Hence (2) is equivalent to harmonicity off the origin.
Step 3 (solve (2)). Suppose first that v′ never vanishes on an interval. Then (2) may be divided by v′ and reads (log∣v′∣)′=(1−n)/r, whose solutions are log∣v′∣=(1−n)logr+const, that is
v′(r)=rn−1a(4)
for a constant a=0. Equation (4) also holds, with a=0, when v′≡0; and these are the only possibilities, since (2) is a linear first-order equation for v′ and (4) is its general solution. (Equivalently, multiply (2) by rn−1 to get (rn−1v′)′=0 directly, which is Strauss's form.)
Step 4 (integrate (4)). For n=2, v′=a/r integrates to alogr+c. For n≥3, v′=ar1−n integrates to 2−nar2−n+c. Renaming the constant b in each case gives (3). □
Definition(The fundamental solution)AdaptedAdapted from Evans (2010), §2.2.1a, Eq. (6), p. 22; Stakgold & Holst (2011), §2.5, pp. 175–176 and §8.3, Eqs. (8.3.16)–(8.3.17), p. 496
Among the radial harmonic functions (3) we pick out one, normalized so that it is the potential of a unit source. Following Stakgold and Holst, a fundamental solution for a linear operator L with pole at ξ is a distributional solution of
LE=δ(x−ξ),
which by their criterion means (E,L∗ϕ)=ϕ(ξ) for every test function ϕ∈C0∞(Rn). Such an E is never unique — any solution of the homogeneous equation may be added — and for an isotropic problem one singles out a particular one by requiring spherical symmetry about the pole and prescribed behaviour at infinity. For a constant-coefficient L it is enough to find the solution with pole at the origin and translate: E(x,ξ)=E(x−ξ,0).
For L=−∇2 this selection gives Evans's fundamental solution of Laplace's equation
defined for x=0. Here, as in harmonic functions, α(n) denotes the volume of the n-dimensional unit ball, so that a ball and a sphere of radius r measure ∣B(x,r)∣=α(n)rn and ∣∂B(x,r)∣=nα(n)rn−1. In three dimensions n(n−2)α(n)=3⋅1⋅34π=4π, so
Φ(x)=4π∣x∣1(n=3),
which is Stakgold and Holst's E(x,ξ)=1/(4π∣x−ξ∣). By construction Φ is harmonic on Rn∖{0}, and so is x↦Φ(x−y) away from x=y, since Laplace's equation is translation invariant. Evans records the size of its derivatives: there is a constant C with
∣DΦ(x)∣≤∣x∣n−1C,∣D2Φ(x)∣≤∣x∣nC(x=0).(6)
Sign conventions, which differ between our sources. Evans and Stakgold–Holst normalize against −∇2, so that −∇2Φ=δ0 and Φ>0 for n≥3: the potential of a positive point source is positive. Riley, Hobson and Bence normalize against ∇2 itself, defining F by ∇2F(r,r0)=δ(r−r0), and therefore get the opposite sign,
So F=−Φ. We follow Evans and Stakgold–Holst: in this unit −∇2Φ=δ0. Arfken writes the same content as the identity ∇2(1/r)=−4πδ(r), proved in the Dirac delta in calculations, which is −∇2Φ=δ0 divided by 4π.
Theorem(The normalizing constant is a flux)AdaptedAdapted from Riley, Hobson & Bence (2006), §21.5.3, Eqs. (21.87)–(21.90), p. 757
Let n≥2 and let Φ be as in (5). Then for every r>0
∮∂B(0,r)∂n∂ΦdS=−1,(7)
where ∂/∂n is the derivative along the outward normal of the ball. This single equation is what fixes the constants in (5), and it is how Riley, Hobson and Bence derive them: integrate ∇2F=δ over a ball centred on the pole, get 1 on the right, and convert the left side by the divergence theorem into the flux of ∇F through the sphere; radial symmetry then turns that flux into (area)×F′(r) and determines F.
Step 1 (the normal derivative).Φ is radial, Φ(x)=v(r), so on ∂B(0,r) the outward normal is radial and ∂Φ/∂n=v′(r), a constant on the sphere.
Step 2 (evaluate v′). For n≥3, v(r)=n(n−2)α(n)1r2−n, so
v′(r)=n(n−2)α(n)2−nr1−n=−nα(n)1r1−n.
For n=2, v(r)=−2π1logr gives v′(r)=−2πr1, which is the same formula since nα(n)=2π when n=2.
Step 3 (multiply by the area). The sphere ∂B(0,r) has area nα(n)rn−1, so
∮∂B(0,r)∂n∂ΦdS=nα(n)rn−1⋅(−nα(n)1r1−n)=−1,
independently of r, which is (7). □
RemarkOriginal
Equation (7) is the whole content of the normalization, and it also explains why the flux is independent of r: by part 2 of Green's identities, the flux of ∇Φ through ∂B(0,r) equals ∫B(0,r)∇2Φ, which gets no contribution from the region where Φ is harmonic. All of the −1 sits at the origin. Put the other way round: no classical function can have ∇2Φ=0 everywhere and a nonzero flux through a closed surface, so the singularity is not a blemish on Φ but the source itself.
Theorem(The Newtonian potential solves Poisson's equation)AdaptedAdapted from Evans (2010), §2.2.1b, Thm 1, pp. 22–23; Stakgold & Holst (2011), §8.3, Eq. (8.3.18), p. 496
Let f∈Cc2(Rn): twice continuously differentiable with compact support. Define the Newtonian potential of f as the convolution
In three dimensions (8) reads u(x)=∫f(y)dy/(4π∣x−y∣), which is Stakgold and Holst's superposition of point-source potentials: the volume element at y carries source strength f(y)dy and contributes f(y)dy/(4π∣x−y∣) to the potential at x. Evans notes that the smoothness demanded of f is far more than is needed, and refers to Gilbarg and Trudinger for weaker hypotheses.
Idea of the proofOriginal
The formula is the superposition one would write down without thinking, but the obvious proof of it is wrong, and seeing why is the point. One cannot simply move two derivatives inside the integral in (8): by the Hessian estimate (6), ∣D2Φ∣ is of size ∣x−y∣−n near the singularity, which is not integrable in n dimensions, so the differentiated integral does not even converge.
The repair has two moves. First, write the convolution the other way round, u(x)=∫Φ(y)f(x−y)dy, so that the derivatives fall on f, which is smooth and compactly supported; difference quotients then converge uniformly and give u∈C2 with ∇2u(x)=∫Φ(y)∇x2f(x−y)dy. Second, cut a ball B(0,ε) out around the singularity. On the small ball, Φ is integrable — that is the whole difference between Φ and D2Φ — so that piece is O(ε2) or O(ε2∣logε∣) and disappears. Outside the ball, Φ is smooth and harmonic, so Green's second identity moves both derivatives from f onto Φ, killing the volume term and leaving only boundary terms on the sphere ∂B(0,ε). One of those terms carries a factor Φ and vanishes with ε; the other carries ∂Φ/∂ν, which by the normalization (7) is exactly 1/∣∂B(0,ε)∣, so that term is an average of f over a shrinking sphere, and it tends to −f(x). That limit is the theorem.
ProofAdaptedAdapted from Evans (2010), §2.2.1b, pp. 23–24
Throughout, C denotes a constant that may change from line to line but does not depend on ε.
Step 1 (Φ is locally integrable). Near the origin ∣Φ(y)∣≤C∣y∣2−n for n≥3 and C∣log∣y∣∣ for n=2. In polar coordinates dy=nα(n)ρn−1dρ, so
In particular the integral (8) converges absolutely for every x: this is Stakgold and Holst's remark that the singularity of the kernel is mild enough to be cancelled by the volume element.
Step 2 (move the derivatives onto f). Substituting y↦x−y in (8),
u(x)=∫RnΦ(y)f(x−y)dy.(11)
Fix i and h=0. Subtracting two copies of (11),
hu(x+hei)−u(x)=∫RnΦ(y)hf(x+hei−y)−f(x−y)dy.
Because f∈Cc2, its first derivatives are uniformly continuous, so the difference quotient of f converges to fxi(x−y) uniformly in y as h→0; and Φ is integrable on the (bounded) region where the quotients are not zero, by Step 1. Hence we may pass to the limit under the integral:
uxi(x)=∫RnΦ(y)fxi(x−y)dy.
Repeating the argument with fxi in place of f, which is legitimate because fxi∈Cc1,
The right side of (12) is continuous in x, again by uniform continuity of fxixj and integrability of Φ, so u∈C2(Rn). This is part 1. Note that no derivative has been taken of Φ so far.
Step 3 (split off the singularity). Summing (12) over i=j and splitting the domain at radius ε,
By (9), ∣Iε∣≤∥D2f∥∞∫B(0,ε)∣Φ∣≤Cε2 for n≥3 and Cε2∣logε∣ for n=2; either way Iε→0.
Step 4 (integrate by parts once). In Jε the x-Laplacian of f(x−y) equals its y-Laplacian, because each differentiation in y brings one factor −1 and they appear in pairs. On the region Ωε=Rn∖B(0,ε) — intersected with a large ball containing the support of y↦f(x−y), outside which everything vanishes — the function Φ is smooth and f(x−⋅) is C2, so Green's first identity applies:
where ν is the unit normal on ∂B(0,ε) pointing into the ball, i.e. outward from Ωε (the sign in Green's identity is fixed by that choice, and there is no contribution from the large sphere, where f vanishes). Call the two terms Kε and Lε. By (10), ∣Lε∣≤∥Df∥∞∮∂B(0,ε)∣Φ∣dS≤Cε for n≥3 and Cε∣logε∣ for n=2, so Lε→0.
Step 5 (integrate by parts again; the volume term dies). Applying Green's first identity to Kε in the other direction,
substituting z=x−y, which maps the sphere of radius ε about 0 onto the sphere of radius ε about x. The right side is minus the average of f over a sphere shrinking to x, and f is continuous, so by the argument of Step 5 of harmonic functions it tends to −f(x).
Step 7 (collect).∇2u(x)=Iε+Kε+Lε for every ε>0, and the left side does not depend on ε. Letting ε→0+ and using Steps 3, 4 and 6,
∇2u(x)=0+(−f(x))+0=−f(x),
that is −∇2u=f. □
RemarkOriginal
Where the hypotheses are used, and what fails without them.
Compact support of f, in Steps 2 and 4. It makes (8) converge at infinity, makes the difference quotients vanish outside a bounded set, and kills the outer boundary term in Step 4. Without decay the convolution need not exist at all: in three dimensions f≡1 gives a divergent integral.
f∈C2, in Steps 2, 3 and 4. Two derivatives are moved onto f, and ∥D2f∥∞ is used to bound Iε. Evans states that much less suffices, and the formula remains correct for f merely Hölder continuous with compact support, but that is not what is proved here.
The excision of B(0,ε), in Step 5. This is the hinge. The term ∫∇2Φf vanishes only because the origin has been removed; if one keeps it and writes ∇2Φ=−δ0 in the sense of distributions, the same term produces −f(x) directly, which is Evans's "formal computation" and Arfken's way of using ∇2(1/r)=−4πδ(r). The proof above is that formal computation made legitimate.
n≥2. The constants in (5) are singular at n=2 (where r2−n is constant) and the logarithm replaces the power. In n=1 there is no singular radial solution at all, and the fundamental solution of −d2/dx2 is −21∣x∣, which grows rather than decays.
Behaviour at infinity, and why two dimensions is different. For n≥3, Φ→0 at infinity, and this is the extra condition that picks Φ out of the family (3): Riley, Hobson and Bence impose F→0 as ∣r∣→∞ and conclude that the additive constant is zero. For n=2 they point out that this condition cannot be imposed, because log∣x∣ is unbounded in both directions, and the additive constant stays undetermined. The non-uniqueness is exactly Stakgold and Holst's remark that a fundamental solution is determined only up to a solution of the homogeneous equation; adding a constant, or any entire harmonic function, gives another fundamental solution.
What the theorem does not say. It gives a solution of −∇2u=f on all of Rn, not the only one: u+h solves the same equation for any harmonic h. Uniqueness needs a condition at infinity (for n≥3, decay, by the maximum principle applied on large balls), or boundary data on a bounded region, which is the subject of Green's functions for Laplace and Poisson problems.
Example(A smooth source in three dimensions, solved twice)Original
Take n=3 and the radial source
f(x)={(1−r2)3,0,r=∣x∣≤1,r>1,
which is in Cc2(R3) but not C3: with s=r2 the profile (1−s)3 and its first two derivatives in s vanish at s=1, while the third does not. So Evans's hypotheses hold exactly and no more.
Route 1: solve the ODE. For a radial u, −∇2u=−r21(r2u′)′, so −∇2u=f says (r2u′)′=−r2f and hence r2u′(r)=−M(r) with
the last value from 31−53+73−91=315105−189+135−35. Outside the support M is constant, so u′=−M(1)/r2 and the decaying solution is
u(r)=315r16(r≥1).
Integrating u′=−M(r)/r2 inwards from r=1 gives, for r≤1,
u(r)=81−6r2+203r4−14r6+72r8.
(Check the match at r=1: 81−61+203−141+721=31516=0.0507937, and the two expressions for u′ both give −31516 there. Check the equation: u′′+2u′/r=−1+3r2−3r4+r6=−(1−r2)3, so −∇2u=f inside, as required.)
Route 2: evaluate the convolution at the centre. Formula (8) at x=0 with Φ=1/(4π∣y∣) and dy=4πs2ds gives
The two routes agree: u(0)=81 exactly. Notice that the integrand s(1−s2)3 is bounded even though Φ blows up at the origin — the factor s2 from the volume element beats the factor s−1 from the kernel, which is Step 1 of the proof in one line.
The far field. The total source strength is Q=∫fdy=4πM(1)=64π/315, and outside the support u=16/(315r)=Q/(4πr)=QΦ(x). A localized source looks like a point source of the same strength from far away; the systematic version of this statement is the multipole expansion.
Example(The same problem in the plane: the potential does not decay)Original
Take n=2 and a uniform source on the disc of radius a: f=ρ0 for r≤a and f=0 for r>a. (This f is discontinuous, so it falls outside Evans's hypotheses; we solve the ODE directly and verify the result, and the point of the example is the behaviour at infinity, not the regularity.) With −∇2u=−r1(ru′)′,
u(r)=4ρ0(a2−r2)(r≤a),u(r)=−2ρ0a2logar(r≥a).
Check. Inside, u′=−ρ0r/2 and u′′=−ρ0/2, so u′′+u′/r=−ρ0 and −∇2u=ρ0. Outside, u is a multiple of logr, which is harmonic by (3), so −∇2u=0. At r=a both expressions give u=0, and both give u′=−ρ0a/2: the solution is C1 across the edge, as it must be, and its second derivative jumps there by exactly the jump in f.
What this shows. The total source strength is Q=πa2ρ0, and the exterior potential is −2πQlog(r/a), which tends to −∞ as r→∞ rather than to 0. In two dimensions the Newtonian potential of a nonzero net source cannot be made to vanish at infinity, which is why Riley, Hobson and Bence keep an undetermined constant in the planar fundamental solution and impose no condition at infinity. The field −u′(r)=Q/(2πr) does decay, and it is the field, not the potential, that is physical here: this is the potential of a uniformly charged infinite cylinder of radius a, read in the plane transverse to its axis.
ApplicationOriginal
Electrostatics and gravitation. Arfken's route to the fundamental solution is the physicist's: Gauss's law ∇⋅E=ρ/ε0 with E=−∇φ gives Poisson's equation ∇2φ=−ρ/ε0, and a point charge q at the origin means ρ=qδ(r). Comparing with the Coulomb potential φ=q/(4πε0r) produces ∇2(1/r)=−4πδ(r), which Arfken emphasises is meaningful only under an integral sign. In the normalization of this unit, φ=Φ∗(ρ/ε0), which is Riley, Hobson and Bence's Eq. (21.91),
φ(x)=∫4πε0∣x−y∣ρ(y)dV(y),
the integral over all space. The same formula with −4πG in place of 1/ε0 is the Newtonian gravitational potential of a mass distribution, which is where the name "Newtonian potential" comes from. The theorem above is the statement that this universally used formula really does solve the field equation, and the proof is where one sees that the 1/r singularity of each point contribution is harmless after integration but fatal if one differentiates first.
Boundary-value problems. Riley, Hobson and Bence introduce the fundamental solution not for its own sake but as the singular part of a Green's function: for a region V with boundary S one writes G=F+H with F the fundamental solution, carrying the delta, and H harmonic in V, chosen to make G satisfy the boundary condition. Every construction in Green's functions for Laplace and Poisson problems and the method of images begins this way, and the representation formula they derive from Green's second identity reduces to (8) when the boundary is pushed to infinity and the data vanish there.
Screened potentials and other operators. The pattern of the definition — solve the homogeneous equation away from the pole, then fix the strength of the singularity by the integrated form of the equation — is Stakgold and Holst's general recipe, and it is not special to the Laplacian. For −∇2+q2, which governs steady diffusion in an absorbing medium (and, in particle physics, a massive scalar field), the same recipe gives e−qr/(4πr) in three dimensions — the Yukawa or Debye-screened potential — and K0(q∣x∣)/(2π) in two, as Stakgold and Holst set as their Exercise 2.5.3. The Laplacian is the q→0 limit, and the comparison shows what the long range of the Coulomb and Newtonian potentials is really due to: the absence of a mass or screening term.
Numerics: boundary integral and fast multipole methods. Because (8) is a convolution with an explicit kernel, a charge distribution given on a mesh of N points produces a potential computable at all N points in O(NlogN) or O(N) operations instead of O(N2), by grouping distant sources — the fast multipole method, and more broadly the boundary-element method, in which an elliptic problem on a region is converted into an integral equation on its boundary using Φ as the kernel. The local integrability established in Step 1 is what makes those singular boundary integrals convergent, and the inverse-power bounds (6) are what make the far-field grouping accurate.
Solves the Dirichlet problem for Laplace's equation in a disk by separating variables and summing the resulting Fourier series into the Poisson kernel, proves from the kernel's three properties that the integral has the prescribed continuous boundary values, records the corresponding formula for a ball in three dimensions, and deduces that every harmonic function is infinitely differentiable inside its domain.
Definition(The Dirichlet problem for a disk, and the Poisson kernel)AdaptedAdapted from Strauss (2008), §6.3, Eqs. (1)–(2), p. 165 and Eq. (17), p. 170; Stein & Shakarchi (2003), §2.5.4, Lemma 5.5, p. 55; Stakgold & Holst (2011), §8.3, Eq. (8.3.2), p. 491
Let D={x∈R2:∣x∣<a} be the open disk of radius a and C=∂D its boundary circle, with polar coordinates (r,θ) on D. Given a function h on C, the Dirichlet problem for the disk is
For a=1 and with the angle written as a single variable, Stein and Shakarchi's normalization is
Pr(θ)=n=−∞∑∞r∣n∣einθ=1−2rcosθ+r21−r2,
the same function: they arrive at it as the kernel whose convolution with f produces the Abel meansAr(f)(θ)=∑nr∣n∣aneinθ of the Fourier series of f. So the Poisson kernel of the disk and the Abel-summation kernel of Fourier series are one object, and Poisson's formula below is Abel summation of the boundary data's Fourier series.
Stakgold and Holst write the interior Dirichlet problem for −∇2 and warn about an artefact of polar coordinates: the lines θ=π and θ=−π are the same ray and r=0 is an interior point, so a candidate solution must be checked for good behaviour at both — which is where the periodic conditions and the boundedness requirement at the origin come from, rather than from any physical boundary.
Theorem(Poisson's integral formula for a disk)AdaptedAdapted from Strauss (2008), §6.3, Thm 1 with Eqs. (13)–(15), pp. 168–169; Ahlfors (1979), §4.6.3, Thm 22, Eq. (64), p. 168 and §4.6.4, Thm 23, p. 169
Let h be continuous on the circle C={r=a} and define, for r<a,
Then u is harmonic in D, it is infinitely differentiable there, and
x→x0limu(x)=h(x0)for every x0∈C,(4)
so that u extends to a continuous function on D with boundary values h. It is the only harmonic function in D with this property.
Writing x for the point with polar coordinates (r,θ) and x′ for the boundary point with polar coordinates (a,φ), the law of cosines gives ∣x−x′∣2=a2+r2−2arcos(θ−φ), and since arc length on C is ds=adφ, formula (3) takes the geometric form
u(x)=2πaa2−∣x∣2∮∣x′∣=a∣x−x′∣2u(x′)ds.(5)
Ahlfors states (5) the other way round — as a representation of a given harmonic function through its boundary values, for u harmonic in ∣z∣<R and continuous on ∣z∣≤R — and then proves separately (his Theorem 23, due to Schwarz) that the integral on the right has the right boundary limits whenever the data are merely piecewise continuous and continuous at the point in question. Stein and Shakarchi prove the same for any integrable f on the circle, with the boundary limit holding at every point of continuity of f and uniformly if f is continuous.
Idea of the proofOriginal
Two different things have to happen, and it is worth keeping them apart.
Finding the formula. Separating variables in polar coordinates gives Θ′′+λΘ=0 with periodic conditions, hence λ=n2, and an Euler equation for the radial factor, hence r±n (and 1,logr for n=0). Half of these blow up at the centre and are discarded — the regularity requirement at r=0 is doing the work that a boundary condition would do. What survives is rncosnθ and rnsinnθ, and matching at r=a identifies the coefficients as the Fourier coefficients of h divided by an. Substituting them back and exchanging the sum with the integral produces an explicit series in (r/a)ncosn(θ−φ), which is the real part of two geometric series and therefore sums in closed form. The result is the kernel (2).
Proving that the formula works. Harmonicity is easy, because each term of the series is harmonic and one may differentiate the integral (3) under the integral sign: the integrand is smooth in x for x inside and x′ on the circle, so they never collide. The boundary limit is the real content, and it follows from three properties of the kernel: it is positive, it has total mass one, and for angles bounded away from zero it becomes uniformly small as r→a. A kernel with those three properties concentrates all of its mass near θ=0, so the integral (3) is a weighted average of values of h taken almost entirely from a small arc around the point being approached; continuity of h then pins the average to h at that point. Uniqueness is immediate from the maximum principle.
Step 1 (separate). In polar coordinates ∇2u=urr+r−1ur+r−2uθθ (see cylindrical and spherical coordinates). Substituting u=R(r)Θ(θ), dividing by RΘ and multiplying by r2 gives
Rr2R′′+rR′=−ΘΘ′′.
The two sides depend on different variables, so each is a constant λ:
Θ′′+λΘ=0,r2R′′+rR′−λR=0.(6)
Step 2 (the angular factor). The coordinate θ is an angle, so u must be unchanged by θ↦θ+2π; the natural conditions on Θ are therefore the periodic ones, Θ(θ+2π)=Θ(θ) for all θ (Stakgold and Holst impose Θ(−π)=Θ(π) and Θ′(−π)=Θ′(π), which is the same requirement). The solutions of Θ′′+λΘ=0 are periodic with period 2π exactly when λ=n2 with n=0,1,2,…, and then
Θ(θ)=Acosnθ+Bsinnθ(n≥1),Θ=A(n=0).
Step 3 (the radial factor). With λ=n2 the radial equation in (6) is of Euler type, so we try R=rα: substituting gives α(α−1)+α−n2=0, that is α2=n2 and α=±n. For n≥1 this gives R=Crn+Dr−n; for n=0 the two roots coincide and the second solution is logr, so R=C+Dlogr.
Step 4 (reject the singular solutions). We want u harmonic on all of D, including the centre. Both r−n and logr are unbounded as r→0, so D=0 in each case. Strauss calls this the "boundary condition" at r=0; it is not a boundary condition at all, but the demand that the solution be a genuine function on the disk. Superposing what is left,
u(r,θ)=21A0+n=1∑∞rn(Ancosnθ+Bnsinnθ).(7)
Step 5 (match the boundary data). Setting r=a in (7), the requirement u=h on C reads
h(θ)=21A0+n=1∑∞an(Ancosnθ+Bnsinnθ),
which is the Fourier series of h with coefficients anAn and anBn. Hence (Fourier series)
Step 6 (substitute and sum). Putting (8) into (7) and using cosnφcosnθ+sinnφsinnθ=cosn(θ−φ),
u(r,θ)=∫02πh(φ)[1+2n=1∑∞(ar)ncosn(θ−φ)]2πdφ.
The interchange of sum and integral is legitimate because r/a<1 is fixed and the series converges uniformly in φ (h is bounded, being continuous on a circle). The bracket is the second expression in (2); to sum it, write cosnψ=21(einψ+e−inψ) with ψ=θ−φ and use the geometric series with ratio (r/a)e±iψ, of modulus r/a<1:
since the two added terms contribute 2arcosψ−2r2. So the bracket equals (a2−r2)/(a2−2arcosψ+r2)=P(r,ψ), which turns the display into (3). □
Proof(The three kernel properties, and the boundary limit)AdaptedAdapted from Strauss (2008), §6.3, properties (i)–(iii) and Eqs. (18)–(20), pp. 170–171; Stein & Shakarchi (2003), §2.5.4, Lemma 5.5, p. 56
so the denominator of P is positive and bounded below, and since a2−r2>0,
P(r,θ)>0(r<a).
The middle expression in (9) is Stein and Shakarchi's split; the right-hand one is Strauss's, and it is the sharper of the two for what follows.
Step 2 (unit mass). Integrating the series in (2) term by term — legitimate for fixed r<a, where it converges uniformly — and using ∫02πcosnθdθ=0 for n≥1,
∫02πP(r,θ)2πdθ=1(0≤r<a).(10)
(Ahlfors obtains the same identity as the special case u≡1 of the representation formula.)
Step 3 (concentration). Let δ>0. For δ≤∣θ∣≤π we have sin2(θ/2)≥sin2(δ/2), so by (9)
0<P(r,θ)≤4arsin2(δ/2)a2−r2(δ≤∣θ∣≤π).(11)
The right-hand side tends to 0 as r→a−, uniformly in such θ. So: given δ>0 and η>0, there is r0<a with P(r,θ)<η for all r0<r<a and all δ≤∣θ∣≤π.
Step 4 (harmonicity). Each term rncosnθ, rnsinnθ of (7) is harmonic — they are the real and imaginary parts of zn, or directly: ∇2(rncosnθ)=n(n−1)rn−2cosnθ+nrn−2cosnθ−n2rn−2cosnθ=0. On any disk r≤ρ<a the series (7) and all of its term-by-term derivatives converge uniformly, because the coefficients anAn,anBn are bounded (they are Fourier coefficients of a bounded function) and (ρ/a)n decays geometrically. Hence u may be differentiated term by term any number of times, and ∇2u=0 on r<ρ for every ρ<a, so on D. The same argument applied to the series in (2) shows that P(r,θ) is itself harmonic in the disk, which is Strauss's third property of the kernel: each term (r/a)ncosnθ is harmonic, and so is the sum.
Step 5 (the boundary limit: set-up). Fix θ0. By (10),
Let ε>0 and let H=maxC∣h∣, finite since h is continuous on a circle. By continuity of h at θ0 choose δ>0 with ∣h(φ)−h(θ0)∣<ε whenever ∣φ−θ0∣<δ.
Step 6 (split the integral). Split (12) at ∣φ−θ0∣=δ and bound each piece.
Near part. On ∣φ−θ0∣<δ the bracket is smaller than ε in absolute value, and P>0, so this piece is at most ε∫02πP(r,θ0−φ)2πdφ=ε by (10). (Positivity, Step 1, is what allows the bracket to be bounded and the kernel left alone.)
Far part. On ∣φ−θ0∣≥δ the bracket is at most 2H, and by Step 3 the kernel is less than η once r is close enough to a; so this piece is at most 2Hη.
Choosing η=ε gives ∣u(r,θ0)−h(θ0)∣≤ε(1+2H) for all r sufficiently close to a. Since ε was arbitrary, u(r,θ0)→h(θ0) as r→a−.
Step 7 (from radial to unrestricted limits). The estimate in Step 6 is uniform in θ0 if δ is chosen by uniform continuity of h, which holds since C is compact. Hence u(r,θ)→h(θ) uniformly in θ, and therefore for any sequence xk→x0∈C with xk=(rk,θk),
the first term by uniformity and the second by continuity of h. This is (4), and it also shows u extends continuously to D.
Step 8 (uniqueness). If u1,u2 are harmonic in D and both satisfy (4) with the same h, then each extends continuously to D with the same boundary values, so w=u1−u2 is harmonic in D, continuous on D, and zero on C. Since D is bounded, the weak maximum principle of harmonic functions gives maxDw=minDw=0, hence w≡0. □
RemarkOriginal
Ahlfors's derivation of (5) is shorter than the one above and uses no Fourier series at all, at the price of assuming that a solution exists. The mean value property gives the value of a harmonic function at the centre of the disk; a Möbius transformation of the disk onto itself carries any interior point to the centre, and transplanting the mean value identity through it produces (5), with the kernel appearing as the Jacobian dθ∗/dθ of the induced map of the boundary circle. Ahlfors even gives the resulting picture: to compute PU(z), replace each boundary value U(θ) by the value at the point diametrically opposite with respect to z — the second intersection with the circle of the line through eiθ and z — and take the ordinary average. Our Step 6 is instead his and Schwarz's arc-splitting argument in the language of kernels.
Theorem(Poisson's formula for a ball in three dimensions)AdaptedAdapted from Strauss (2008), §7.4, Eqs. (16)–(17), p. 195
Let D={x∈R3:∣x∣<a} and let h be continuous on the sphere ∣x∣=a. Then the solution of ∇2u=0 in D, u=h on ∂D, is
Compare the disk, formula (5): the kernel is ∣x−x0∣−2 over 2πa in two dimensions and ∣x−x0∣−3 over 4πa in three, in each case multiplied by a2−∣x0∣2. At the centre, x0=0, formula (13) collapses to u(0)=4πa21∮hdS, the mean value property.
NoteOriginal
Strauss obtains (13) not by separating variables but from the Green's function of the ball, which he constructs by reflecting the source point x0 across the sphere to x0∗=a2x0/∣x0∣2 and taking G=−4πρ1+∣x0∣a4πρ∗1 with ρ=∣x−x0∣ and ρ∗=∣x−x0∗∣; formula (13) is then ∮h∂G/∂ndS. That construction belongs to the method of images, and the general representation formula it uses to Green's functions for Laplace and Poisson problems, so we state (13) here and do not reprove it. Strauss notes that the same reflection carried out in two dimensions, starting from 2π1logr, reproduces formula (5) — the disk formula obtained in this unit by a completely different route.
Theorem(Harmonic functions are infinitely differentiable)AdaptedAdapted from Strauss (2008), §6.3, 'Differentiability', p. 170; Evans (2010), §2.2.3b, Thm 6, p. 27
Let U be any open subset of the plane and let u be harmonic in U. Then u has partial derivatives of every order at every point of U; that is, u∈C∞(U).
Evans proves the corresponding statement in Rn under a weaker hypothesis — u∈C(U) satisfying the mean value property for every ball in U — and warns that nothing follows about the boundary: u may fail to extend smoothly to ∂U, and may fail even to extend continuously. The conclusion is interior and only interior.
Idea of the proofOriginal
Smoothness is not visible in Laplace's equation, which mentions only two derivatives; it has to be imported from a formula. Poisson's formula is that formula. Its integrand, as a function of the interior point x, is a smooth function of x — the only thing that could spoil it is the vanishing of ∣x−x′∣, and that cannot happen because x is inside and x′ is on the circle. So one may differentiate under the integral sign as often as one likes. Since every point of an open set is the centre of some small disk contained in it, and a function differentiable on each of a family of disks covering U is differentiable on U, the local statement gives the global one.
Step 1 (a disk about the point). Let x0∈U. Since U is open, choose ρ>0 with the closed disk B(x0,ρ)⊂U, and let B=B(x0,ρ).
Step 2 (u on B is a Poisson integral).u is harmonic on B and continuous on B (it is C2 on U⊇B). Let v be the Poisson integral (5) on B built from the boundary values u∣∂B. By Step 8 of the earlier proof — uniqueness — v=u on B. So on B,
u(x)=2πρρ2−∣x−x0∣2∮∂B∣x−x′∣2u(x′)ds(x′).(15)
Step 3 (differentiate under the integral sign). Fix a closed disk B(x0,σ) with σ<ρ. For x in it and x′∈∂B we have ∣x−x′∣≥ρ−σ>0, so the integrand of (15) is a C∞ function of x there, and every x-derivative of it is bounded on B(x0,σ)×∂B (a product of compact sets) and continuous. The hypotheses of the theorem on differentiating an integral with respect to a parameter are therefore met, for derivatives of every order, and
for every multi-index α, each such derivative being continuous in x. Hence u∈C∞(B(x0,σ)), and in particular all derivatives of u exist at x0.
Step 4 (globalize).x0 was an arbitrary point of U, and differentiability of all orders is a pointwise property, so u∈C∞(U). □
RemarkOriginal
Where the hypotheses are used, and what fails without them.
Continuity of h, in Steps 5–6 of the kernel proof. It is used only at the point θ0 being approached, which is why Ahlfors's Theorem 23 and Stein and Shakarchi's Theorem 5.7 can take h merely piecewise continuous or integrable and still conclude the limit at every point of continuity. The example below has a two-point jump and illustrates exactly this.
Unrestricted versus radial limits, in Step 7. Condition (4) demands u(x)→h(x0) as x→x0 in any manner, and that is what makes the solution unique. If one asks only for the radial limit limr→1u(r,θ)=h(θ) for each fixed θ, uniqueness fails: Stein and Shakarchi's Exercise 18 points out that
u(r,θ)=∂θ∂Pr(θ)=(1−2rcosθ+r2)2−2r(1−r2)sinθ
is harmonic in the unit disc (a derivative of a harmonic function with respect to the rotation parameter), satisfies limr→1u(r,θ)=0 for every θ, and is not identically zero. It is unbounded near z=1, so it has no continuous extension to D and (4) fails for it at that one boundary point. One boundary point is enough.
Boundedness of the domain, in Step 8 only. Uniqueness is the maximum principle's contribution and inherits its need for a bounded domain.
Existence can fail for other domains. Poisson's formula settles the disk and the ball completely, and it is tempting to assume a continuous h on a nice-looking boundary always has a solution. Stakgold and Holst give three warnings. For the punctured unit disk 0<∣x∣<1, whose boundary consists of the circle ∣x∣=1 plus the one extra point 0, the data h=0 on the circle and h=1 at the puncture admit no solution at all. For the unit disk with suitable continuous h, ∣∇u∣ can blow up at the boundary, so the solution lies in C2(Ω)∩C(Ω) but not in C2(Ω). And for a three-dimensional domain with a sharp inward spike — in the extreme case, the unit ball with a radial segment removed — there are continuous data with no solution; they interpret this physically as charge leaking out along the spike. Their sufficient condition is the exterior sphere condition: if at every boundary point there is an open ball outside Ω whose boundary touches Γ at that point and nowhere else, then continuous data on Γ always admit a harmonic solution lying in C2(Ω)∩C(Ω).
The kernel's two-sided bound gives more than the boundary limit. Ahlfors notes the elementary inequalities, in our notation with ∣x∣=r and radius a,
a+ra−r≤∣x−x′∣2a2−r2≤a−ra+r,
which follow from a−r≤∣x−x′∣≤a+r. For u≥0 they turn (5) into Harnack's inequalitya+ra−ru(0)≤u(x)≤a−ra+ru(0),
since the average of u over the circle is u(0). Ahlfors emphasises that this holds only for positive harmonic functions, and uses it to prove Harnack's principle: a monotone sequence of harmonic functions either tends to +∞ uniformly on compact sets or converges to a harmonic limit, uniformly on compact sets.
Conventions. Strauss writes the kernel as P(r,θ) for a disk of radius a and divides by 2π inside the integral; Stein and Shakarchi write Pr(θ) for the unit disc and define convolution with the same 2π1; Ahlfors writes the kernel as (R2−∣a∣2)/∣z−a∣2 with the interior point called a and the radius R. Only names differ.
Example(A jump on the boundary circle: harmonic measure of a semicircle)Original
Take a=1 and the discontinuous data
h(φ)={1,0,0<φ<π,−π<φ<0,
so h is continuous except at the two points φ=0 and φ=π. Poisson's formula still makes sense — h is bounded and integrable — and gives
u(r,θ)=2π1∫0πP(r,θ−φ)dφ.
Evaluation. Write z=reiθ. Then
1−z1+z=∣1−z∣2(1+z)(1−zˉ)=∣1−z∣2(1−r2)+2irsinθ,
whose real part is positive for ∣z∣<1; so w=log1−z1+z is analytic on the unit disc with the principal branch, and
Imw=arg1−z1+z=arctan1−r22rsinθ.
Set
U(r,θ)=21+π1arctan1−r22rsinθ.
U is harmonic: it is an affine function of the imaginary part of an analytic function. U has the right boundary values: at z=eiφ with sinφ=0, the real part of (1+z)/(1−z) is (1−r2)/∣1−z∣2=0 while the imaginary part has the sign of sinφ, so arg=±π/2 and U→21±21, which is 1 on the upper semicircle and 0 on the lower. U is bounded, by 0≤U≤1. By the uniqueness argument of Step 8, applied on the disc minus two boundary points (or simply by checking against the integral), u=U.
Numerical check of the identification. Evaluating the integral 2π1∫0πP(r,θ−φ)dφ numerically against U:
r
θ
integral
U(r,θ)
0
0.3
0.500000
0.500000
0.5
π/2
0.7951672
0.7951672
0.5
0.3
0.6194765
0.6194765
0.9
−1.0
0.0397219
0.0397219
At the centre, u(0,θ)=21, which is the mean value property: half the boundary circle carries the value 1. At r=21, θ=2π the exact value is
What the example illustrates. The function u is C∞ everywhere inside the disc, as the smoothness theorem demands, even though its boundary data have two jumps: smoothness is an interior statement and knows nothing about the boundary. The boundary limit (4) holds at every point of the circle except z=±1, in agreement with Ahlfors's and Stein–Shakarchi's versions of the theorem, which ask for continuity of the data only at the point approached. At z=1 the limit genuinely fails: approaching along the upper semicircle gives 1, along the lower gives 0, and approaching along the real axis gives 21 — the average of the two one-sided values, which is the general rule for a jump (Stein and Shakarchi's Exercise 17). The level sets u=c are the circular arcs through ±1, so u is the harmonic measure of the upper semicircle.
Example(A ball in three dimensions, checked against an explicit harmonic polynomial)Original
Take a=1 in (13) and the axially symmetric data h=cos2θ on the unit sphere, i.e. h=z2 in Cartesian coordinates.
The answer, found by inspection. Write cos2θ=31(2P2(cosθ)+1), using P2(t)=21(3t2−1). The function
substituting t=cosθ. The integral evaluates to exactly 34 (numerically 1.3333333), so u(x0)=83⋅34=21; and the polynomial gives 31+41(3−1)/3=31+61=21. The two agree.
The centre, and the mean value property. Setting x0=0 in (13) gives u(0)=4π1∮cos2θdS=21∫−11t2dt=31, which is the constant term of the polynomial, as it must be: the degree-2 solid harmonic r2P2(cosθ) vanishes at the origin and averages to zero on every sphere about it.
ApplicationOriginal
Heat, and why the Dirichlet problem is the natural one. Stein and Shakarchi reach Laplace's equation in the disc as the long-time limit of the heat equation: hold a prescribed temperature distribution f on the rim of a disc, wait, and the interior settles into the steady state u with ∇2u=0 and u=f on the rim. Poisson's formula says what the equilibrium temperature at an interior point is — a weighted average of the rim temperatures, weighted by (a2−r2)/∣x−x′∣2, which favours the part of the rim nearest the point and does so more and more sharply as the point approaches the rim. The same formula with the same reading describes the electrostatic potential inside a cylindrical capacitor whose surface potential is prescribed and which is long enough to be treated as two-dimensional, which is Strauss's reading of his formula (14).
Probability: harmonic measure and exit distributions. Because the kernel is positive and integrates to one, Poisson's formula exhibits u(x) as the expectation of h against a probability measure on the boundary — the harmonic measure seen from x. That measure is exactly the exit distribution of a Brownian particle started at x: the probability that it first touches the circle in the arc A is 2πa1∫A∣x−x′∣2a2−∣x∣2ds. The worked example computes exactly this for A a semicircle, and the answer 21+π1arctan1−r22rsinθ is the probability that a Brownian particle started at reiθ leaves the disc through the upper half of the circle. This is the probabilistic reading of Strauss's fourth motivating example for Laplace's equation, and it underlies walk-on-spheres Monte Carlo solvers, which sample exit points of successive spheres instead of discretizing the equation.
Summability of Fourier series. Read with the roles reversed, (3) is a statement about Fourier series, and this is Stein and Shakarchi's reason for introducing it. A continuous function's Fourier series may diverge at some points; but u(r,θ), which is its Abel mean ∑nr∣n∣aneinθ, converges to f at every point of continuity and uniformly for continuous f. So "evaluate the Fourier series inside the disc and then let r→1" always works where termwise summation on the circle may not. Damping the nth coefficient by r∣n∣ is the oldest and simplest of the summability methods, and solving a Dirichlet problem is what it amounts to geometrically.
Interior regularity as a tool. The smoothness theorem is used constantly without being mentioned: it is what licenses differentiating a harmonic function as many times as an argument needs, and what makes the term-by-term manipulations of multipole and spherical-harmonic expansions legitimate inside the region of harmonicity (multipole expansions). It also marks the difference between elliptic and hyperbolic equations: Strauss points out that the analogous smoothing holds for the diffusion equation but is simply false for the wave equation, where a kink in the initial data travels along a characteristic for ever (the one-dimensional wave equation). In numerical analysis, interior smoothness is why spectral methods converge so fast for elliptic problems away from the boundary, and why the loss of accuracy in such a computation is almost always a boundary effect.
Gives the recipe for reducing Laplace's and Helmholtz's equations to ordinary differential equations by seeking product solutions, proves the lemma that makes the separation constants constant, carries the spherical case through in full, and catalogues which special functions - trigonometric, hyperbolic, Bessel, Legendre and spherical harmonic - each of the Cartesian, cylindrical and spherical systems produces.
Definition(Separable solutions and separation constants)AdaptedAdapted from Riley, Hobson & Bence (2006), §21.1, Eq. (21.1), pp. 713–714; Hassani (2013), §13.1, pp. 396–397
Let a linear partial differential equation be given in independent variables x1,…,xN. A solution of the product form
u(x1,…,xN)=X1(x1)X2(x2)⋯XN(xN)(1)
is called separable in those variables, and looking for solutions of the form (1) is the method of separation of variables. The requirement is only that each factor depend on one variable alone; the factors may, and usually do, share numerical parameters. Riley, Hobson and Bence illustrate the distinction with three functions of (x,y,z,t): xyz2sinbt separates completely; (x2+y2)zcosωt separates in the last two variables only, since the first two occur in a sum that cannot be factored out; and xy+zt separates in none of them.
The parameters that appear when (1) is substituted into the equation are the separation constants. In the examples below they turn out to be the eigenvalues of the ordinary differential operators that the method produces, so the special functions of mathematical physics enter as the eigenfunctions.
Throughout we use the Laplacian in the three standard systems (cylindrical and spherical coordinates): with ρ the distance from the z-axis, φ the azimuth and r the distance from the origin,
The two equations to be separated are Laplace's equation∇2u=0 and the Helmholtz equation∇2u+k2u=0, which reduces to Laplace's when k=0.
Method(The separation lemma)Original
Let A⊆Rp and B⊆Rq be nonempty sets, let F:A→R and G:B→R, and suppose
F(s)+G(t)=0for all s∈A,t∈B.(2)
Then F and G are each constant, with G≡−λ and F≡λ for a single number λ.
Idea of the proofOriginal
This is the only genuine piece of logic in the method, and it is a one-line argument: in (2) the variables s and t may be moved independently, so freezing one of them exhibits the other's function as a constant. Everything else in separation of variables is calculus and bookkeeping. It is worth stating the lemma separately because the textbook phrase "a function of x alone equals a function of y alone, so both are constant" is used dozens of times in what follows, and because it makes visible what the method actually needs: that the variables range over a product set, which is exactly the condition that fails on a domain not bounded by coordinate surfaces.
ProofOriginal
Step 1. Fix any t0∈B and set λ=−G(t0). Putting t=t0 in (2) gives F(s)=−G(t0)=λ for every s∈A. So F is the constant λ.
Step 2. Substituting F≡λ back into (2) gives G(t)=−λ for every t∈B. So G is the constant −λ. □
Choose coordinates that fit the boundary. Pick a coordinate system in which every piece of the boundary is a level surface of one coordinate: a box in Cartesians, a disk or annulus or wedge in plane polars, a cylinder in cylindrical polars, a ball or spherical shell in spherical polars. This is the step that decides whether the method will work at all.
Separate the time, if there is any. Write Ψ(r,t)=u(r)T(t) and let S denote the whole spatial operator, so that the equation reads SΨ=∂Ψ/∂t (heat, Schrödinger) or SΨ=∂2Ψ/∂t2 (wave, Klein–Gordon). Dividing by uT puts a function of r on one side and a function of t on the other; by the separation lemma both equal a constant, Hassani's α, and one is left with an ODE in t and a homogeneous equation in space,
∇2u(r)+f(r)u(r)=0,(3)
which is Hassani's master equation. (Hassani writes the product as Ψ=RT and then lets the letter Ψ stand for the spatial factor as well in his master equation; we reserve Ψ for the full solution and call the spatial factor u, which is the letter steps 3–7 and the rest of this unit use, so that u in (3) and u in step 3 are the same object.) For the heat and wave equations f is the constant k2 and (3) is Helmholtz's equation; for the Schrödinger equation f(r)=2m(E−V(r))/ℏ2. Hassani notes that α is conventionally taken real except for the Schrödinger equation, where a purely imaginary α is more convenient because it cancels the i in the equation, and that in every case the admissible α are fixed by the boundary conditions, not by the separation. This time-separation constant is used nowhere else in the unit; the letter α is reused below, in the cylindrical catalogue and the cylinder example, for the unrelated axial separation constant.
Posit a product and divide by it. Substitute u=X1X2⋯XN into (3) and divide the result by u. Then multiply by whatever factor (typically r2 or ρ2) is needed to make one variable appear in one group of terms only.
Introduce separation constants. Apply the separation lemma to each split, naming the constants in whatever algebraic form will be convenient later — −n2, k2, ℓ(ℓ+1). With N variables there are N−1 independent constants, and in a Cartesian separation of the wave equation they satisfy one algebraic relation among themselves, as in Riley, Hobson and Bence's example where μ2=l2+m2+n2.
Impose the conditions the coordinates hide. These are not boundary conditions in the physical sense, and forgetting them is the commonest error:
Single-valuedness in an angle.Φ(φ+2π)=Φ(φ) forces the azimuthal separation constant to be m2 with m an integer. (The exception is a genuinely multivalued potential, such as the magnetic scalar potential of a current, where Riley, Hobson and Bence keep the term Φ=Cφ+D with C=I/2π.)
Regularity at the origin or on the axis. If the region contains ρ=0 or r=0, discard ρ−m, lnρ, r−(ℓ+1), Ym(kρ), yℓ(kr). If it contains the polar axis θ=0,π, discard Qℓm(cosθ).
Behaviour at infinity. On an unbounded exterior region, discard the solutions that grow.
What is left after these deletions is an eigenvalue problem of Sturm–Liouville type in each variable, and the admissible separation constants are its eigenvalues.
Read off the special functions. Each separated ODE is one of a short list, and its solutions are the classical special functions; the catalogue is in the note below.
Superpose, and fit the inhomogeneous data by orthogonality. A single product rarely satisfies the given data. Sum over all admissible separation constants and determine the coefficients from the one inhomogeneous condition by the orthogonality of the eigenfunctions: Fourier coefficients for cosmφ,sinmφ, Fourier–Bessel coefficients for Jm(knρ) (orthogonality and Fourier–Bessel series), Legendre coefficients for Pℓ(cosθ), and the full double sum over ℓ and m in spherical harmonics (spherical harmonics) when the data are not axially symmetric.
Note(Catalogue of separated equations)AdaptedAdapted from Riley, Hobson & Bence (2006), §21.1, Eqs. (21.6)–(21.7), p. 715, §21.3.1, Eqs. (21.27)–(21.49), pp. 726–734, and §21.3.2, Eqs. (21.52)–(21.57), pp. 737–741; Arfken, Weber & Harris (2013), §9.4, Eqs. (9.29)–(9.41), pp. 415–417, Eqs. (9.58)–(9.63), pp. 421–422, and Eqs. (9.74)–(9.84), pp. 424–427
For the Helmholtz equation ∇2u+k2u=0; Laplace's equation is the case k=0 throughout. In each system the product is substituted and the equation divided by it, which is legitimate only where no factor vanishes; as in the proof below, the separated ODEs listed here then hold on the whole of each coordinate interval because both sides of each are continuous.
Cartesian, u=X(x)Y(y)Z(z). With constants −l2,−m2,−n2,
X′′=−l2X,Y′′=−m2Y,Z′′=−n2Z,l2+m2+n2=k2.
Each factor is Aeilx+Be−ilx: a sine and cosine pair if l2>0; a pair of exponentials e±βx, equivalently sinhβx and coshβx, if l2<0, where β=−l2 is real and l=iβ is imaginary (the exponent is the real number −l2, not l2=±l, which is imaginary here); and the linear function A+Bx if l2=0. For Laplace's equation (k=0) the three constants sum to zero, so they cannot all share a sign: either one of them vanishes, or one has the opposite sign to another. So Laplace's equation in a box never produces an oscillating factor in every direction; some direction carries hyperbolic or linear behaviour instead, which is why a Dirichlet datum on one face decays into the box.
Plane polars, u=P(ρ)Φ(φ). With constant n2,
Φ′′+n2Φ=0,ρ2P′′+ρP′+(k2ρ2−n2)P=0.
Single-valuedness makes n an integer. For k=0 the radial equation is the Euler equation ρ2P′′+ρP′−n2P=0 with solutions ρn,ρ−n (and 1,lnρ when n=0); for k=0 the substitution μ=kρ turns it into Bessel's equation of order n, with solutions Jn(kρ) and Yn(kρ). Riley, Hobson and Bence warn against confusing Yn, the Bessel function of the second kind, with a spherical harmonic Yℓm; the latter carries two indices.
Cylindrical polars, u=P(ρ)Φ(φ)Z(z). With constants m2 (azimuthal) and α2 (axial),
The sign convention in the axial constant. We write the axial equation as Z′′=α2Z, so that a real α gives exponential behaviour in z and an imaginary α oscillatory behaviour; substituting back into ∇2u+k2u=0 then puts k2+α2, not k2−α2, in the radial equation. Riley, Hobson and Bence write the axial factor as Eeiαz+Fe−iαz (p. 740), that is Z′′=−α2Z, and their radial argument is correspondingly k2−α2. The two conventions differ only by α↦iα; everything here and in the cylinder example below is in ours.
For Laplace's equation (k=0) the radial argument is α and the solution is [AJm(αρ)+BYm(αρ)][Ccosmφ+Dsinmφ][Eeαz+Fe−αz]: Bessel functions radially, trigonometric functions azimuthally, exponentials axially. (Riley, Hobson and Bence call the axial constant k in this case and write the same solution as [AJm(kρ)+BYm(kρ)][Ccosmφ+Dsinmφ][Ee−kz+Fekz], Eq. (21.35).) When k2+α2<0 the radial argument becomes imaginary and the Bessel functions are replaced by the modified Bessel functionsIm,Km.
Spherical polars, u=R(r)Θ(θ)Φ(φ). With constants ℓ(ℓ+1) and m2,
The polar equation becomes the associated Legendre equation under μ=cosθ, with solutions Pℓm(μ) and Qℓm(μ) (associated Legendre functions), finite polynomial solutions requiring ℓ a non-negative integer and 0≤∣m∣≤ℓ; the product ΘΦ regular on the axis is the spherical harmonic Yℓm(θ,φ). The radial equation for k=0 is the Euler equation with solutions rℓ and r−(ℓ+1); for k=0 the substitution R=r−1/2S with μ=kr reduces it to Bessel's equation of order ℓ+21, so that R=r−1/2Jℓ+1/2(kr) or r−1/2Yℓ+1/2(kr), which after normalization are the spherical Bessel functionsjℓ(kr) and yℓ(kr).
Method(The separated solutions in spherical polars)AdaptedAdapted from Riley, Hobson & Bence (2006), §21.3.1, Eqs. (21.42)–(21.49), pp. 732–734 and §21.3.2, Eqs. (21.56)–(21.57), pp. 740–741
Let u=R(r)Θ(θ)Φ(φ) be a solution of Laplace's equation ∇2u=0 of product form on a region that is a product of coordinate intervals in r∈(0,∞), θ∈(0,π) and φ, where the spherical coordinates are regular. Let the three factors be twice continuously differentiable, and let none of them vanish identically. (Note what this last hypothesis does not say: it does not make u itself nonvanishing, and indeed u vanishes somewhere inside the region in all but the simplest cases. Step 1 of the proof shows that this weaker hypothesis is nevertheless all the separation needs.)
1. What the separation alone gives. There are constants m2 and λ=ℓ(ℓ+1) — complex numbers, subject to no further restriction — such that Φ′′=−m2Φ, such that M(μ)=Θ(θ) with μ=cosθ solves the associated Legendre equation with those parameters, and such that R solves the Euler equation r2R′′+2rR′−ℓ(ℓ+1)R=0. Those three ODEs are the entire content of the separation.
Naming their solutions takes a little care, because the parameters are so far unrestricted. Here Pℓm and Qℓm denote that equation’s two standard solutions — Legendre’s function of first kind, and the one of second kind — taken at complex degree ℓ and order m; associated Legendre functions defines them only for an integer ℓ≥0 with ∣m∣≤ℓ, which is the only case used from part 2 onwards. Generically — that is, for parameters with ℓ=−21 and m=0 at which Pℓm and Qℓm are moreover linearly independent —
Each of the three provisos is a real restriction, and the third is not claimed to be a complete list. When ℓ=−21 the two radial exponents ℓ and −(ℓ+1) coincide and the radial factor is r−1/2(A+Blnr) instead. When m=0 the azimuthal factor is Cφ+D. And Pℓm,Qℓm really can fail to be an independent pair: if ℓ is a non-negative integer and m an integer with m>ℓ, then Pℓm≡0 and the second solution has to be found some other way. No such parameter pair arises in part 2, where ℓ and m are integers with ∣m∣≤ℓ.
2. With the two regularity conditions. Suppose in addition that
u is single-valued in the azimuth, Φ(φ+2π)=Φ(φ), and
u is finite on the polar axis θ=0,π.
Then m is an integer, ℓ is a non-negative integer with 0≤∣m∣≤ℓ, and F=0, so that
The three factors are then linked only through ℓ and m. Finiteness at the origin sets B=0 as well, and decay at infinity sets A=0.
Both extra conditions are needed. Without single-valuedness, u=φ is harmonic and of product form but is not of the form (4a) — it is the m=0, C=0 case, and it is exactly the magnetic scalar potential kept in step 5 of the method. Without finiteness on the axis, u=r1/2P1/2(cosθ) is a product solution with ℓ=21, so the indices need not be integers at all. This does not contradict part 2, because the second condition fails: P1/2 is perfectly finite at the north pole, P1/2(1)=1, but diverges logarithmically at the south pole θ→π, where cosθ→−1. Numerically P1/2(−0.9)=−0.50127, P1/2(−0.99999)=−3.4946 and P1/2(−1+10−7)=−4.9605 — slow divergence, but divergence.
The Helmholtz equation. For ∇2u+k2u=0 the angular factors and the two conditions on them are unchanged, and only the radial factor differs:
Under the same two regularity conditions, ℓ and m are integers with 0≤∣m∣≤ℓ and the coefficient of Qℓm vanishes; finiteness at the origin sets B=0 in addition. (Riley, Hobson and Bence call the solution F in their Eq. (21.57), where it collides with their coefficient F of Qℓm; we write the solution as u, as in (4), so that F here is the coefficient only.)
Riley, Hobson and Bence are careful to point out that nothing in the derivation of (4) forces ℓ to be an integer; integrality is imposed afterwards, by the requirement that the polar factor be finite on the axis, and is then inherited by the radial factor.
Idea of the proofOriginal
The separation is done in two cuts, and the order matters. Multiplying the equation by r2 after dividing by u isolates all the r-dependence in one group: that is the first cut, and its separation constant is written ℓ(ℓ+1) because the radial equation it leaves is an Euler equation whose indices are then the tidy pair ℓ,−(ℓ+1). The remaining equation involves θ and φ only; multiplying it by sin2θ isolates the φ-dependence, which is the second cut, with constant m2. Single-valuedness in φ makes m an integer, and what is left for θ is the associated Legendre equation once cosθ is used as the variable in place of θ.
The one subtlety is that the two conditions fixing the integers come from different places. The integrality of m comes from the geometry of the azimuth — φ and φ+2π are the same direction. The integrality of ℓ comes from the behaviour of the Legendre equation at its singular points cosθ=±1, which are the poles of the sphere: only for integer ℓ does that equation have a solution bounded there.
Step 1 (first cut, and where the division is legal). Substitute u=R(r)Θ(θ)Φ(φ) into ∇2u=0 in spherical polars. One wants to divide by u, and u is not nonvanishing: in all but the simplest cases it is zero on a surface inside the region, because Pℓm(cosθ) and cosmφ have interior zeros and Arℓ+Br−(ℓ+1) has one whenever AB<0. So divide only where the division is legal. Put
IR={r:R(r)=0},IΘ={θ:Θ(θ)=0},IΦ={φ:Φ(φ)=0}.
Each is open in its own coordinate interval, the three factors being continuous, and each is nonempty — this, and only this, is what the hypothesis that no factor vanishes identically gives, and it is all that is needed. Their product IR×IΘ×IΦ is exactly the subset of the region where u=0, and it is again a product set, which is what the separation lemma requires. On it, dividing by u and multiplying by r2 gives
depends on r onlyR1drd(r2drdR)+depends on (θ,φ) onlyΘsinθ1dθd(sinθdθdΘ)+Φsin2θ1dφ2d2Φ=0,
the factor sin2θ being nonzero because 0<θ<π on the region. By the separation lemma, with s=r ranging over IR and t=(θ,φ) over IΘ×IΦ, each group is a constant; call the first λ, so the second is −λ.
Step 1a (an identity true where the division was legal is true everywhere). Each separated equation below is first obtained only on the set where its own factor is nonzero, and then holds on the whole coordinate interval, by the following argument, which we use without further comment. Take the radial equation r2R′′+2rR′−λR=0, which multiplying the first group of Step 1 by R gives on IR. Let Z be the complement of IR in the r-interval, a closed set. A point of ∂Z is a limit of points of IR, and both sides of the equation are continuous in r (R being C2), so the equation passes to the limit and holds there. On the interior of Z the function R vanishes on a neighbourhood, so R=R′=R′′=0 and both sides are 0. Since the interval is IR∪∂Z∪intZ, the equation holds throughout it. The same argument applies verbatim to Θ on its interval and to Φ on its own.
Step 2 (the radial equation). By Steps 1 and 1a, on the whole r-interval
r2R′′+2rR′−λR=0,
an Euler equation. Substituting r=et and R(r)=S(t), so that rd/dr=d/dt and r2d2/dr2=d2/dt2−d/dt, it becomes the constant-coefficient equation S′′+S′−λS=0, with solutions eλ1t,eλ2t where λ1+λ2=−1 and λ1λ2=−λ. Hence
R(r)=Arλ1+Brλ2,λ1+λ2=−1.
Writing λ1=ℓ, the relation λ1+λ2=−1 gives λ2=−(ℓ+1) and λ=−λ1λ2=ℓ(ℓ+1). Nothing so far requires ℓ to be an integer, or even real: it is a name for the first root.
The two roots coincide exactly when λ1=λ2=−21, that is when ℓ=−21 and λ=−41. There S′′+S′−λS=0 has the repeated root −21, so its general solution is S(t)=(A+Bt)e−t/2 and
R(r)=r−1/2(A+Blnr),
which is the exceptional radial factor recorded in part 1 of the statement. Every other value of ℓ gives two distinct powers.
Step 3 (second cut). The second group of Step 1, with λ=ℓ(ℓ+1), reads, on IΘ×IΦ,
By the separation lemma again, applied on the product set IΘ×IΦ, each of the two groups is a constant: Φ′′/Φ≡−m2 on IΦ, and the θ-group ≡m2 on IΘ.
Step 4 (the azimuthal factor). Multiplying by Φ gives Φ′′+m2Φ=0 on IΦ, hence on the whole φ-interval by Step 1a. It gives Φ=Ccosmφ+Dsinmφ for m=0 and Φ=Cφ+D for m=0. Since φ and φ+2π label the same half-plane, a single-valued u requires Φ(φ+2π)=Φ(φ): hence C=0 in the m=0 case, and m must be an integer in the other.
Step 5 (the polar factor). With the constant m2 inserted, the θ-group of Step 3 reads, on IΘ,
Θsinθdθd(sinθdθdΘ)+ℓ(ℓ+1)sin2θ=m2.
Change the variable to μ=cosθ. Then dμ/dθ=−sinθ, so d/dθ=−sinθd/dμ=−(1−μ2)1/2d/dμ, and sinθdΘ/dθ=−(1−μ2)dM/dμ where M(μ)=Θ(θ). Substituting, multiplying through by M and dividing by sin2θ=1−μ2,
dμd[(1−μ2)dμdM]+[ℓ(ℓ+1)−1−μ2m2]M=0,(6)
which is the associated Legendre equation — on the image of IΘ under θ↦cosθ to begin with, and then, by Step 1a, on the whole interval. Its solutions are M=EPℓm(μ)+FQℓm(μ).
Step 6 (integrality of ℓ, and the axis). Equation (6) is singular at μ=±1, which are the two poles θ=0,π of the sphere. The function Qℓm diverges there, so a solution finite on the polar axis has F=0; and Pℓm is finite there only when ℓ is a non-negative integer with 0≤∣m∣≤ℓ, which is the condition established in associated Legendre functions. This is where integrality of ℓ enters, and through λ=ℓ(ℓ+1) it feeds back into the radial exponents of Step 2.
Step 7 (assemble). Multiplying the three factors gives (4), which is part 1 of the statement; imposing the single-valuedness of Step 4 and the finiteness on the axis of Step 6 gives part 2 and the form (4a).
Step 8 (the Helmholtz case). Repeating Steps 1–6 for ∇2u+k2u=0 changes only Step 1, where the extra term contributes k2r2 to the r-group. So Steps 3–6 are untouched — the angular factors are identical — and the radial equation becomes
r2R′′+2rR′+[k2r2−ℓ(ℓ+1)]R=0.(7)
Put R(r)=r−1/2S(r). Then R′=r−1/2S′−21r−3/2S and R′′=r−1/2S′′−r−3/2S′+43r−5/2S, so
r2R′′+2rR′=r3/2S′′+r1/2S′−41r−1/2S,
and multiplying (7) by r1/2 gives
r2S′′+rS′+[k2r2−(ℓ+21)2]S=0,
using ℓ(ℓ+1)+41=(ℓ+21)2. With μ=kr this is Bessel's equation of order ℓ+21, whose solutions are Jℓ+1/2(kr) and Yℓ+1/2(kr). Hence (5). □
Remark(Where the hypotheses are used: nonvanishing factors, and a product domain for the separation lemma)Original
Two hypotheses of the statement are consumed in Step 1, and nowhere else. The first is that no factor vanishes identically. It is far too weak to make u=0 on the region — no hypothesis short of assuming that outright would — but it is exactly enough to make the set IR×IΘ×IΦ on which one may divide by u nonempty, and continuity then carries each separated equation across the zeros of its own factor. The second is the product form of the region: Step 1 applies the lemma with t=(θ,φ) ranging over its angular factor independently of r. If the region were, say, a ball dented on one side, the admissible (θ,φ) would depend on r, the two groups would not be free to vary independently, and the conclusion that each is constant would be unjustified. (The restriction to IR×IΘ×IΦ does no damage here, because cutting a product set down factor by factor leaves a product set.) This is the precise reason the method demands a region bounded by coordinate surfaces, and why step 1 of the recipe is the decisive one.
Remark(Separation finds solutions; it does not prove they are all of them)AdaptedAdapted from Evans (2010), §4.1.1, pp. 159–161
A single product (4) is one solution; a sum of products is another. That a superposition can meet arbitrary data is a completeness statement about the eigenfunctions, proved separately — Fourier series for the trigonometric factors, Fourier–Bessel series for the radial ones, the completeness of spherical harmonics for the angular ones. Uniqueness then comes from the maximum principle or an energy argument (harmonic functions), not from the separation. Evans makes the same point about his separated heat equation: only the single product u=e−λtw is produced by the separation itself, and the finite sums and infinite series built from it rest on the linearity of the equation and on a completeness statement he defers to his later chapters.
Remark(What the sources add: angular factors, time factors, and the Schrödinger exception)AdaptedAdapted from Riley, Hobson & Bence (2006), §21.3.2, pp. 737–741; Hassani (2013), §13.2, Eqs. (13.11)–(13.14), pp. 398–400
Why every equation built from ∇2 has the same angular factors. Riley, Hobson and Bence make this the general point of their section: provided θ and φ enter the equation only through ∇2, the second cut of the proof is untouched, so the angular factors are spherical harmonics whatever the radial term. Hassani's operator version makes it structural: ∇2=−r−2L2+r−2(r⋅∇)(r⋅∇)+r−2r⋅∇ with L=r×(−i∇) depending only on the angles, so that for any f=f(r) the master equation (3) splits into the eigenvalue problem L2Y=λY on the sphere and a radial ODE containing f and λ. (Hassani writes this angular eigenvalue α; we write λ, which is the λ=ℓ(ℓ+1) of the proof above, to keep α for the axial constant.) In quantum mechanics this is the statement that angular momentum is conserved in any central potential, and it is why the hydrogen atom, the three-dimensional oscillator and a free particle all have the same angular wavefunctions and differ only radially.
Which equations separate the time cleanly. Step 2 works because the equation is linear with coefficients independent of t. Riley, Hobson and Bence note that for the wave equation the resulting time factor e±iωt with ω=kc accompanies the spatial Helmholtz solution of the samek, so that a separated solution has a definite frequency; but in spherical polars the zeros of jℓ(kr) are not equally spaced in r except when ℓ=0, so such a solution has no definite wavelength. For the diffusion equation the time factor is the decaying e−k2κt instead.
The Schrödinger equation is close to Helmholtz but not equal to it. Riley, Hobson and Bence point this out for the hydrogen atom: the stationary-state equation is ∇2u+ℏ22m(E+4πε0re2)u=0, which has an r-dependent coefficient rather than the constant k2, so the radial equation is not Bessel's. The angular factors are nevertheless spherical harmonics, by the previous point. See the Schrödinger equation.
Remark(Not every separable system is one of these three)Original
The method also separates in elliptic, parabolic, prolate and oblate spheroidal, and other coordinate systems, each producing its own family of special functions (Mathieu, Lamé, spheroidal wave functions). The three systems treated here cover the geometries that occur most often, and they are the ones whose special functions this handbook develops.
Example(Cartesian: Laplace's equation in a cube)Original
Solve ∇2u=0 in the unit cube 0<x,y,z<1 with u=0 on the five faces x=0, x=1, y=0, y=1, z=0, and
u(x,y,1)=sinπxsinπy
on the top face.
Separate. With u=X(x)Y(y)Z(z) and the constants of the catalogue, X′′=−l2X, Y′′=−m2Y, Z′′=−n2Z and l2+m2+n2=0 (Laplace, so k=0). The four homogeneous side faces force X(0)=X(1)=0 and Y(0)=Y(1)=0, so X=sinpπx and Y=sinqπy with p,q positive integers, and l2=p2π2, m2=q2π2. Then n2=−(p2+q2)π2<0, so the z-factor is hyperbolic, as the catalogue predicts:
Z(z)=Esinhγpqz+Fcoshγpqz,γpq=πp2+q2.
The face z=0 forces F=0.
Fit the data. The datum is already a single eigenfunction, p=q=1, so only one term survives and
The datum has decayed by a factor of about nine over half the height of the cube — the exponential decay of a Dirichlet datum into the interior at the rate set by γ, which is the Cartesian signature of harmonicity. Had the datum been sin3πxsin4πy instead, γ would be 5π and the centre value 1/(2cosh(5π/2))=3.9×10−4: finer structure on the boundary dies faster, which is the mechanism behind interior smoothing.
Example(Cylindrical: Dirichlet eigenvalues of a cylinder)Original
Find the lowest frequencies at which a scalar field can resonate inside a cylinder of radius a=1 and height L=2 while vanishing on the whole boundary; that is, find the smallest k with a nonzero solution of the Helmholtz equation ∇2u+k2u=0 in the cylinder and u=0 on its curved surface, bottom and top.
Separate. From the catalogue, with the region containing the axis ρ=0 we drop Ym, and with u vanishing at z=0 and z=L the axial factor must be sin(pπz/L), so α2=−(pπ/L)2 with p a positive integer. The radial factor is then Jm(κρ) with the catalogue's radial argument κ=k2+α2, that is
κ2=k2+α2=k2−L2p2π2,
and the condition u=0 on ρ=a requires Jm(κa)=0, i.e. κa=jm,n, the nth positive zero of Jm. Hence the admissible values are
Numbers. With a=1 and L=2, the relevant Bessel zeros are j0,1=2.4048256, j1,1=3.8317060, j2,1=5.1356223, j0,2=5.5200781, and (pπ/L)2=2.4674011p2. The seven lowest eigenvalues are
(m,n,p)
jm,n2
(pπ/L)2
k2
k
multiplicity
(0,1,1)
5.783186
2.467401
8.250587
2.872384
1
(0,1,2)
5.783186
9.869604
15.652790
3.956361
1
(1,1,1)
14.681971
2.467401
17.149372
4.141180
2
(1,1,2)
14.681971
9.869604
24.551575
4.954955
2
(0,1,3)
5.783186
22.206610
27.989796
5.290538
1
(2,1,1)
26.374616
2.467401
28.842018
5.370476
2
(0,2,1)
30.471262
2.467401
32.938663
5.739222
1
The modes with m=0 come in pairs, cosmφ and sinmφ, so they are doubly degenerate; those with m=0 are simple. The lowest mode is J0(2.4048256ρ)sin(πz/2), nodeless inside, with k=2.8723835.
A check that the separation constants are doing what they should. Each k2 is a sum of a radial and an axial piece, and the two pieces are the eigenvalues of the two separated operators; the algebraic relation among the separation constants in step 4 of the recipe is here simply κ2+(pπ/L)2=k2. Lengthening the cylinder, L→∞, drives the axial piece to zero and the spectrum down onto the two-dimensional values (jm,n/a)2 — the drumskin eigenvalues — as it must.
Example(Spherical: a sphere held at potential V0cos3θ)Original
A sphere of radius a carries the surface potential h(θ)=V0cos3θ. Find the potential inside and outside, with u bounded at the origin and u→0 at infinity.
Expand the datum in Legendre polynomials. From P1(μ)=μ and P3(μ)=21(5μ3−3μ) we get μ3=52P3(μ)+53P1(μ). So with μ=cosθ,
h=V0[53P1(μ)+52P3(μ)].
Choose the surviving terms of (4). The datum is independent of φ, so m=0 and the polar factors are the Legendre polynomials Pℓ; the solution must be finite on the axis, so F=0. Inside the sphere finiteness at the origin kills r−(ℓ+1); outside, decay at infinity kills rℓ. Hence
Checks. Both are harmonic, being superpositions of the solid harmonics rℓPℓ and r−(ℓ+1)Pℓ of (4). At r=a both reduce to V0[53P1+52P3]=V0cos3θ, so the potential is continuous across the surface. Numerically, on the polar axis θ=0 (where Pℓ(1)=1 for all ℓ):
At the centre, uin(0)=0, which is the mean value property: cos3θ is odd under θ↦π−θ and so averages to zero over the sphere. Far away, uout≈53V0a2cosθ/r2: the leading behaviour is a dipole, because the datum has no ℓ=0 component and hence no net charge — the systematic version of this reading is multipole expansions.
The discontinuity in the normal derivative. Differentiating the two expressions and setting r=a,
The radial derivative therefore jumps by −aV0[59P1+514P3] across the surface, which in electrostatics equals −σ/ε0 for the surface charge density σ that holds the sphere at the prescribed potential. The potential itself is continuous and its tangential derivative is continuous; only the normal derivative jumps, and the jump measures the source sitting on the surface.
Application(Electrostatics of conductors)AdaptedAdapted from Riley, Hobson & Bence (2006), §21.3.1, pp. 734, 735–737
Riley, Hobson and Bence's two spherical examples are the standard pair. An uncharged conducting sphere in a uniform field E must have a potential that approaches −Ercosθ at large r and is constant on r=a; since cosθ=P1(cosθ), only the ℓ=1 terms of (4) can appear, and the two conditions fix them completely, giving v=−Er(1−a3/r3)cosθ — the sphere behaves as a dipole of moment 4πε0a3E. A hollow split sphere with its two hemispheres held at v0 and 0 needs the full Legendre series instead, and the coefficients come out as v0/2,3v0/4a,0,−7v0/16a3,…; on the equator the series sums to v0/2, the average of the two values, because a Legendre series behaves at a jump exactly as a Fourier series does.
The cylindrical example above is a cavity calculation for a single scalar field required to vanish on the whole boundary: an acoustic cavity with pressure-release walls, or — Arfken's reading of the same eigenvalue problem — a particle in a cylindrical box, with k2=2mE/ℏ2. With ω=ck and c=2.998×1010cms−1, a cylinder of radius a=1 cm and height L=2 cm has its lowest such mode at k=2.8723835cm−1, that is ω/2π=ck/2π=13.71 GHz.
An electromagnetic cavity of the same dimensions is not this problem, and resonates lower. A conducting wall forces only the tangential electric field to vanish on it, not every field component, so for a TM mode the axial field Ez∝Jm(jm,nρ/a)cos(pπz/L) need have no axial node: k2=(jm,n/a)2+(pπ/L)2 with p=0,1,2,…, whereas the TE modes are built from the zeros jm,n′ of Jm′ and need p≥1. The lowest mode of this cylinder is therefore TM010, with k=j0,1/a=2.4048256cm−1 and ck/2π=11.47 GHz, independent of the height L; the lowest TE mode, TE111, has k=[(j1,1′/a)2+(π/L)2]1/2=2.4202cm−1 with j1,1′=1.8411838, giving 11.55 GHz. The separation itself is word for word the one above; only the admissible axial index changes.
The same separation with the z-dependence dropped gives the modes of a circular drumskin, whose frequencies are ω=jm,nc/a with c2=T/σ; Riley, Hobson and Bence compute the four lowest and note that the second and third are two-fold degenerate. The degeneracy is not an accident of the numbers but a consequence of the rotational symmetry of the disk, which pairs cosmφ with sinmφ; breaking the symmetry — an elliptical drum — splits each pair. See normal modes and standing waves and normal modes.
Application(Steady heat conduction with mixed data)AdaptedAdapted from Riley, Hobson & Bence (2006), §21.3.1, pp. 729–731
Riley, Hobson and Bence's semi-infinite cylinder, curved surface at 0∘ and base at T0, is the model problem for which a Fourier–Bessel series rather than a Fourier series is required: the base datum is expanded in J0(knρ) with kna the zeros of J0, and the coefficients come out as 2T0/(knaJ1(kna)) using the recurrence dzd[zJ1(z)]=zJ0(z). The same calculation with a non-axisymmetric base temperature needs the double series over m and n, with coefficients extracted by the orthogonality of both the Bessel and the trigonometric factors — step 7 of the recipe in its general form.
ApplicationOriginal
Quantum mechanics of central potentials. The reason the periodic table has the shape it has is step 5 and step 6 of this recipe. For any potential V(r) depending on r alone, the angular factors of the stationary Schrödinger equation are the spherical harmonics Yℓm, with ℓ a non-negative integer and 2ℓ+1 values of m; the potential enters only the radial equation and so fixes the energies but not the angular degeneracy. The s,p,d,f labels are the values ℓ=0,1,2,3, and the 2(2ℓ+1) electrons per subshell are the 2ℓ+1 values of m times two spin states. The operator derivation of Hassani (Hassani (2013), §13.2, pp. 398–400) shows why: L2 commutes with any ∇2+f(r), so angular momentum and energy can be diagonalized together. See the Schrödinger equation and the Helmholtz equation.
Numerical methods. Separation is also the basis of the fast direct solvers. On a box, a rectangular grid, or a disk, the discrete Laplacian is diagonalized by the discrete analogues of the separated eigenfunctions, so a Poisson solve becomes a forward transform, a division by the eigenvalues l2+m2+n2, and an inverse transform — O(NlogN) work instead of the O(N3/2) or worse of an iterative solve (the discrete Fourier transform). The same spectral idea, applied with Chebyshev or Legendre expansions in a non-periodic direction, is what makes spectral methods for elliptic problems on tensor-product domains converge faster than any fixed power of the mesh size.
Expands the potential of a localized source outside a sphere containing it, first in Legendre polynomials from the generating function and then in spherical harmonics by the addition theorem, defines the monopole, dipole and quadrupole moments that the terms carry, and shows that the lowest nonvanishing moment is independent of the choice of origin while all higher ones are not.
A point dipole at ξ with axis the unit vector l is the limit, as ε→0, of a source of strength +1/ε at ξ+2εl together with one of strength −1/ε at ξ−2εl. Stakgold and Holst take the limit in the sense of distributions: for a test function ϕ,
ε1[ϕ(ξ+2εl)−ϕ(ξ−2εl)]⟶∂l∂ϕ(ξ),
so the point dipole is the distribution ϕ↦∂ϕ/∂l(ξ). Equivalently, its volume source density is a derivative of the delta: (∂/∂ξi)δ(x−ξ) is the density of a unit dipole at ξ with axis along xi. Note that the charges must be made stronger as they are brought together, or nothing survives.
Arfken takes the same limit at the level of the potential: the two charges ±q at z=±a produce a potential whose a→0 limit at fixed μ=2qa is
ψ=4πε0μr2P1(cosθ),
the potential of a point dipole of moment μ at the origin, oriented along the polar axis. A linear 2n-pole is obtained by iterating the construction n times, and contributes Pn(cosθ)/rn+1; n=0,1,2 are the monopole, dipole and quadrupole.
Throughout, r> and r< denote the larger and the smaller of two radii, and Ω=(θ,φ) a direction. We use the physicists' spherical harmonics Yℓm of spherical harmonics, in which the addition theorem reads as in (3) below.
Application(The Laplace expansion of the inverse distance)AdaptedAdapted from Arfken, Weber & Harris (2013), §15.3, Eqs. (15.60)–(15.63), pp. 736–737 and §16.3, Eq. (16.66), p. 800
Let r1,r2 be measured from a common origin, with magnitudes r1,r2 and angular coordinates Ω1,Ω2, and let χ be the angle between them. Then
which holds in a coordinate system of arbitrary orientation, the angle χ having been traded for the two directions separately. The addition theorem used here is
Pℓ(cosχ)=2ℓ+14πm=−ℓ∑ℓYℓm(Ω1)∗Yℓm(Ω2).(3)
Arfken notes that the rate of convergence of (1) improves as r>/r< grows, and that the expansion of a point charge's potential inside the sphere through the charge — the case r<=r — is obtained from the same generating function by rearranging the square root the other way.
Idea of the proofOriginal
Everything follows from the Legendre generating function
(1−2ht+h2)−1/2=ℓ=0∑∞Pℓ(t)hℓ(∣h∣<1),
because the law of cosines makes ∣r1−r2∣ exactly of that shape. Writing ∣r1−r2∣2=r12+r22−2r1r2cosχ and pulling out whichever of r1,r2 is larger leaves (1−2ht+h2)−1/2 with t=cosχ and h=r</r>, which is less than one. Pulling out the smaller radius would give h>1 and a divergent series, which is why the expansion must be organized around r> and why a multipole expansion is a statement about the field outside the source, never inside it.
The passage from (1) to (2) is a change of bookkeeping with one real gain. In (1) the single function Pℓ(cosχ) couples the two directions, so the expansion is only useful when one of them is fixed — Arfken's charge on the polar axis. The addition theorem (3) factorizes that coupling into a sum of products, each factor depending on one direction alone. Once that is done, the source variables and the field variables appear in separate factors, and the sum over the source can be carried out once and for all. The numbers it produces are the multipole moments.
Step 1 (the law of cosines). With χ the angle between r1 and r2,
∣r1−r2∣=(r12+r22−2r1r2cosχ)1/2.
Step 2 (pull out the larger radius). Let r>=max(r1,r2), r<=min(r1,r2) and h=r</r>∈[0,1). The expression inside the square root is symmetric in r1,r2, so
∣r1−r2∣=r>(1−2hcosχ+h2)1/2,
and therefore
∣r1−r2∣1=r>1(1−2hcosχ+h2)−1/2.
Step 3 (insert the generating function). Since ∣h∣<1 and ∣cosχ∣≤1, the generating function of the Legendre polynomials applies with t=cosχ and gives
∣r1−r2∣1=r>1ℓ=0∑∞Pℓ(cosχ)hℓ,
which is (1). The generating-function expansion is guaranteed only for ∣h∣<1; at h=1, that is at r1=r2, absolute convergence is lost, and Arfken excludes that case. Note that ∣r1−r2∣ itself is perfectly finite there unless also χ=0, so what the exclusion limits is the expansion rather than the function.
Step 4 (factorize the angle). Substituting the addition theorem (3) into (1) and exchanging the two finite inner sums with the outer one gives (2) term by term. □
Application(The multipole expansion)AdaptedAdapted from Arfken, Weber & Harris (2013), §16.3, Eqs. (16.72)–(16.74), p. 801 and §15.3, Eq. (15.67), p. 739
Let charges qi sit at points ri, all inside the sphere ∣r∣=a centred at the origin. For a field point with r>a,
The numbers Mℓm are the multipole moments of the distribution. For a continuous density ρ the sums over i become integrals, Mℓm=2ℓ+14π∫ρ(r′)r′ℓYℓm(Ω′)∗dV′.
If the charges lie on the polar axis, at z=ai, the sum collapses to m=0 and (4) becomes Arfken's linear multipole expansion
with μ0 the monopole moment (the total charge), μ1 the dipole moment and μ2 the quadrupole moment. Both expansions converge for r larger than every ri, that is at field points farther from the origin than any part of the source.
Normalizations differ between authors. Arfken warns of this explicitly: his Mℓm of (4) is the choice that makes the expansion formula simplest, and it does not reproduce the traditional scalings of the low-order moments. With his convention
M00=4πi∑qi,M10=34πi∑qizi,
so the monopole moment is 4π times the total charge and M10 is 4π/3 times the z-component of the usual dipole vector. We follow Arfken's convention and state the factors wherever a traditional moment is wanted.
Idea of the proofOriginal
Take the Laplace expansion (2), put r1=ri (a source point) and r2=r (the field point), and note that the hypothesis r>a decides once and for all which radius is the larger: r>=r and r<=ri for every i. The ℓth term then factors into riℓYℓm(Ωi)∗, which involves the source only, times Yℓm(Ω)/rℓ+1, which involves the field point only. Summing over the charges touches only the first factor, and what comes out is a finite list of numbers that is all the far field ever needs to know about the source. That is the whole content of the expansion: the exterior field of a localized source is determined by countably many numbers, and the first few of them dominate.
Step 1 (superpose). By linearity, the potential of the charges is the sum of their individual Coulomb potentials,
ψ(r)=4πε01i∑∣r−ri∣qi,
each term being qi/ε0 times the fundamental solution Φ(x)=1/(4π∣x∣) of Laplace's equation centred at ri, since (qi/ε0)Φ(r−ri)=qi/(4πε0∣r−ri∣).
Step 2 (expand each term). For each i we have ri≤a<r, so in (2) the roles are fixed: r<=ri, r>=r. Hence
Step 3 (exchange the sums). The sum over i is finite, and each expansion converges absolutely for r>a. To see that, bound the inner sum by the Cauchy–Schwarz inequality and then use ∑m=−ℓℓ∣Yℓm(Ω)∣2=(2ℓ+1)/4π, which is the addition theorem (3) at χ=0:
The ℓth term of the expansion of 1/∣r−ri∣ is therefore at most 2ℓ+14π⋅rℓ+1riℓ⋅4π2ℓ+1=riℓ/rℓ+1≤r−1(a/r)ℓ in modulus, a convergent geometric series since a<r. (Applying the pointwise bound ∣Yℓm∣≤(2ℓ+1)/4π to each of the 2ℓ+1 terms separately would leave a factor 2ℓ+1, which Cauchy–Schwarz removes.) So the order of summation may be changed:
Step 4 (the axial case, for ai of either sign). Let every charge sit on the polar axis, at z=ai with ai positive, negative or zero. A charge with ai>0 has ri=ai and θi=0; one with ai<0 has ri=∣ai∣=−ai and θi=π. By part 3 of spherical harmonics the two axis values are
because (−1)ℓ=(sgnai)ℓ when ai<0. (A charge at the origin, ai=0, contributes 0 for ℓ≥1 and qi for ℓ=0, which is what aiℓ gives with 00=1.) At the field point Yℓ0(θ,φ)=(2ℓ+1)/4πPℓ(cosθ). Substituting both into (4),
Application(Only the lowest nonvanishing moment is independent of the origin)AdaptedAdapted from Arfken, Weber & Harris (2013), §15.3, p. 739
Move the expansion centre from 0 to rp, so that each ri is replaced by ri−rp and the moments change accordingly. Then the lowest nonvanishing moment of the distribution is unchanged, while every higher moment is in general changed. In particular:
the monopole moment — the total charge — is always independent of the choice of centre;
the dipole moment is independent of it only when the total charge is zero;
the quadrupole moment is independent of it only when the total charge and the dipole moment both vanish.
Idea of the proofOriginal
A moment is a weighted sum of a solid harmonic over the charges, and a solid harmonic of degree ℓ is a homogeneous harmonic polynomial of that degree. Translating such a polynomial keeps it harmonic and keeps its degree-ℓ part exactly where it was; everything the translation creates is of lower degree. So moving the centre adds to each moment a combination of moments of strictly lower order and nothing else. If the lower ones are all zero there is nothing to add, which is the theorem; if they are not, the added terms are visible, which gives the converse in the cases named.
ProofOriginal
Step 1 (moments as a linear functional on harmonic polynomials). For ℓ≥0 and ∣m∣≤ℓ put
Sℓm(x)=rℓYℓm(Ω)∗,x=rΩ.
By part 7 of spherical harmonics, rℓYℓm is a homogeneous polynomial of degree ℓ in x,y,z which is harmonic on R3, and by the conjugation relation of part 2 there, Yℓm∗=(−1)mYℓ−m, so Sℓm=(−1)mrℓYℓ−m is one too. Write Hℓ for the complex span of Sℓ−ℓ,…,Sℓℓ; relabelling m↦−m shows this is also the span of the solid harmonics rℓYℓm.
Hℓ is in fact all of the homogeneous harmonic polynomials of degree ℓ. For if P is one, then P is harmonic on the whole of R3 and regular at the origin, so the second sentence of part 7 writes it as P=∑l,mAlmrlYlm; each rlYlm is homogeneous of degree l, and homogeneous polynomials of distinct degrees are linearly independent, so matching degrees against the homogeneous P leaves only l=ℓ and puts P in Hℓ.
Define the linear functional Λ(P)=∑iqiP(ri) on polynomials. Then (4) reads Mℓm=2ℓ+14πΛ(Sℓm), and the moments about the new centre rp — got by replacing each ri by ri−rp — are
In this language, all the moments of order ℓ vanish precisely when Λ vanishes on Hℓ.
Step 2 (what a translation does to a harmonic polynomial). Let P be a homogeneous harmonic polynomial of degree N and set Q(x)=P(x−rp). Then Q is a polynomial of degree N, and it is harmonic, because ∇2 is translation invariant. Split Q into its homogeneous parts, Q=∑d=0NQd. Each ∇2Qd is homogeneous of degree d−2, and homogeneous polynomials of distinct degrees are linearly independent, so ∇2Q=0 forces ∇2Qd=0 for every d separately: each homogeneous part of a harmonic polynomial is itself harmonic. Finally, expanding P(x−rp) and collecting terms of top degree gives QN=P, since P is homogeneous of degree N. Hence
P(x−rp)=P(x)+d<N∑Qd(x),Qd∈Hd.(7)
This is the step that the axial moments cannot supply on their own: a shift off the axis destroys the axial symmetry, so the scalar binomial expansion of (ai−zp)n in Step 5 below settles only the axial case.
Step 3 (the lowest nonvanishing moment is unchanged). Let N be the smallest order at which some moment is nonzero, so by Step 1 Λ vanishes on Hd for every d<N. Apply (7) with P=SNm, whose lower homogeneous parts Qd lie in Hd with d<N:
for every m and every rp. So the whole set of order-N moments survives the shift unchanged; in particular it is still not all zero, so the order N itself does not depend on the centre either.
Step 4 (the three named cases). The converses in the three bullets of the statement come out of (7) at low order, and need no choice of phase convention.
Monopole.H0 consists of the constants, so P(x−rp)=P(x) and M~00=M00 for every rp. This is Step 3 with N=0: the total charge is always centre-independent.
Dipole. Every linear form is harmonic, so H1 is the space of all P(x)=αx+βy+γz, and for such a P, P(x−rp)=P(x)−P(rp). Hence Λ(P(⋅−rp))=Λ(P)−P(rp)∑iqi. The three order-1 moments are the values of Λ on the basis S1−1,S10,S11 of H1, so they are all unchanged exactly when Λ(P(⋅−rp))=Λ(P) for every P∈H1; by the displayed identity that happens for every rp if and only if the total charge ∑iqi vanishes, since otherwise a P and an rp with P(rp)=0 move one of them.
Quadrupole. Take the real element S20(x)=5/16π(2z2−x2−y2) of H2 and rp=(xp,yp,zp). Since
2(z−zp)2−(x−xp)2−(y−yp)2=(2z2−x2−y2)−4zpz+2xpx+2ypy+(2zp2−xp2−yp2),
we get, with c=54π5/16π=214π/5,
M~20−M20=c[−4zpi∑qizi+2xpi∑qixi+2ypi∑qiyi+(2zp2−xp2−yp2)i∑qi].
The right-hand side is a polynomial in rp, so it vanishes for every shift only if its linear part does, which kills all three components of the dipole vector ∑iqiri, and its quadratic part does, which kills the total charge. Together with Step 3 this is the third bullet.
Step 5 (the axial moments, and how they change). On the polar axis the general moments collapse to the μn of (5), and the shift is a scalar computation. With the centre moved to z=zp and all charges at z=ai, the binomial theorem gives
which is (7) read in the one-dimensional case: μ~n is μn plus a combination of the strictly lower moments. For n>N it gives μ~n=μn+(Nn)(−zp)n−NμN+(terms with μj,N<j<n), whose second term is nonzero for every zp=0. So moments above the lowest really do move, and the "in general" in the statement cannot be dropped. Arfken reaches the same conclusion from the shape of the binomial expansion of ∣r−rp∣−n, which feeds each term of the expansion into the later ones and never into the earlier ones. □
Remark(How many components a moment has, and why)AdaptedAdapted from Arfken, Weber & Harris (2013), §16.3, pp. 801–802
For a unit charge at (x,y,z) the moments Mℓm of (4) are, in Cartesian form, homogeneous polynomials of degree ℓ. Arfken tabulates them for ℓ≤2; for example M00=(4π)1/2, M10=(4π/3)1/2z, M1±1=∓(2π/3)1/2(x±iy), and M20=(4π/5)1/221(2z2−x2−y2). There are 2ℓ+1 values of m, but the space of homogeneous polynomials of degree ℓ in three variables has dimension (2ℓ+2), which is larger from ℓ=2 on: six monomials x2,y2,z2,xy,xz,yz against five values of m. Arfken resolves the discrepancy by identifying the missing direction: the combination r2=x2+y2+z2 is rotationally invariant and so belongs with M00, but it carries the wrong power of r to contribute to a solution of Laplace's equation, and therefore cannot appear in a potential. The same happens at ℓ=3, where ten monomials split into seven that span the M3m and three with the symmetry of M1m but again the wrong radial behaviour. In Cartesian language this is the statement that the quadrupole tensor may be taken traceless.
Remark(Symmetry makes moments vanish, and it is worth checking before computing)AdaptedAdapted from Riley, Hobson & Bence (2006), §29.1, pp. 1077–1078
Riley, Hobson and Bence open their chapter on representation theory with this, using three molecules. Hydrogen chloride has no symmetry operation to forbid a dipole moment, and has one. Carbon dioxide is linear and symmetric: a rotation about the long axis kills any perpendicular component, and a rotation by π about a perpendicular axis through the carbon, or a reflection in the plane through the carbon, reverses the parallel component while carrying the molecule to an indistinguishable copy — so that component is zero too, and CO2 has no permanent dipole moment. Ozone is bent, and only the component perpendicular to its two-fold axis is forbidden. The moral to carry away is one-sided: a symmetry operation can rule a component out, but no symmetry argument ever produces a nonzero moment. What survives the symmetry test is merely permitted, and whether it is actually there has to be settled by computation or by experiment.
Remark(What the expansion is not)AdaptedAdapted from Arfken, Weber & Harris (2013), §15.3, Eq. (15.62), p. 736 and p. 739
It is an expansion of the field outside a sphere containing the whole source. Arfken's Eq. (15.62) is the companion statement for a field point closer to the origin than the charge: one single charge q at z=a, seen from r<a, with the roles of r and a in (1) exchanged and the series running in powers of r/a. No shell appears in that statement, and none could: inside a uniform spherical shell of radius a carrying total charge q the potential is the constant q/4πε0a, with no Pn series at all. A field point in the middle of a distribution needs both expansions, split at r, as in the Green's-function expansion of Green's functions for Laplace and Poisson problems. Nor is the expansion restricted to electrostatics: Arfken's third closing observation is that it applies to any inverse-square force, so that the shape of the Earth is described by mass multipoles, and the leading term of gravitational radiation is set by the time dependence of the mass quadrupole.
Example(A planar quadrupole: four charges, no monopole, no dipole)Original
Take a>0 and place charges
+q at (±a,0,0),−q at (0,±a,0),
four equal magnitudes at the corners of a square in the equatorial plane, alternating in sign. Work in units with 4πε0=1.
The low moments vanish. The total charge is q+q−q−q=0, so M00=0. The dipole vector is q(a,0,0)+q(−a,0,0)−q(0,a,0)−q(0,−a,0)=0, so M1m=0 for all m. Both vanishings are forced by symmetry: the configuration is carried to itself by r↦−r with the charges unchanged, so every odd-ℓ moment is zero, and by the rotation through π/2 about the z-axis combined with q↦−q, which kills the m=0 terms.
The quadrupole moments. All four charges have ri=a and θi=π/2, with φi=0,π (charge +q) and π/2,3π/2 (charge −q). Since Y2±1∝sinθcosθe±iφ vanishes at θ=π/2, M2±1=0. For m=0, Y20(π/2,φ) is the same constant for all four charges, and the charges sum to zero, so M20=0. Only m=±2 survives. With Y2±2(θ,φ)=15/32πsin2θe±2iφ,
Cross-check against the Cartesian form. The traceless quadrupole tensor Qjk=∑iqi(3xjxk−ri2δjk) of this configuration is diagonal with Qxx=6qa2, Qyy=−6qa2, Qzz=0 (each +q contributes (2,−1,−1)qa2 and each −q contributes (1,−2,1)qa2 to the diagonal). The standard Cartesian quadrupole potential 21r−3Qjkn^jn^k is then 21r−3⋅6qa2(n^x2−n^y2)=3qa2sin2θcos2φ/r3, which is (6).
Numerical check and the rate of convergence. Set a=1, q=1. At the point r=5, θ=π/2, φ=0 — that is, (5,0,0) — the exact potential is
41+61−262=0.02443440,
while (6) gives 3/125=0.02400000. Comparing at several points:
r
θ
φ
exact
from (6)
error
5
π/2
0
0.02443440
0.02400000
4.34×10−4
10
π/2
0
0.00301276
0.00300000
1.28×10−5
20
π/2
0
0.00037539
0.00037500
3.93×10−7
5
π/3
0.3
0.01465879
0.01485604
−1.97×10−4
10
π/3
0.3
0.00185111
0.00185701
−5.89×10−6
Each doubling of r divides the error by about 32=25, so the error decays like r−5 — exactly the size of the next surviving term, ℓ=4. The odd-ℓ terms are absent by the inversion symmetry noted above, so the first correction to a quadrupole here is a hexadecapole, not an octupole.
The moments are origin-independent here. Because the monopole and dipole moments both vanish, the previous statement says the quadrupole moments cannot change if the expansion centre is moved. Checking directly: shifting all four charges by a vector d changes ∑qi(3xjxk−r2δjk) by terms proportional to ∑qixj (the dipole, zero) and ∑qi (the total charge, zero). So Qjk is unchanged, as claimed.
Example(When the dipole moment depends on where you put the origin)Original
Place +2q at z=a and +q at z=−a — two charges of the same sign, so the monopole moment does not vanish. Using the axial moments (5) with the origin at z=0:
μ0=2q+q=3q,μ1=2q⋅a+q⋅(−a)=qa,μ2=2qa2+qa2=3qa2.
Now move the centre to z=zp. By Step 1 of the proof,
The monopole moment is unmoved, as it must be. The dipole moment is not: it is qa about the origin, zero about zp=a/3 — the charge-weighted centre — and −2qa about zp=a. So the sentence "this configuration has dipole moment qa" is meaningless without saying about which point, whereas "this configuration has total charge 3q" is not.
The field does not care. At r=10a on the positive axis, the exact potential (units 4πε0=1, q=a=1) is 2/9+1/11=0.31313131. From the origin, (5) gives μ0/r+μ1/r2+μ2/r3=0.3+0.01+0.003=0.31300000. From the centre zp=1/3, the field point is at distance r′=10−1/3=9.6666667 and the moments are μ~0=3, μ~1=0, μ~2=3−2/3+1/3=2.6666667, giving 3/9.6666667+0+2.6666667/9.66666673=0.31034483+0.00295215=0.31329698. The two truncations miss the exact value by −1.31×10−4 and +1.66×10−4: both are accurate to four digits, and they bracket the answer. Including the next (ℓ=3) term improves each to 0.3131000 and 0.3130934 respectively, against the exact 0.3131313. Moving the origin reshuffles the moments but is bound to leave the sum of the series alone; the choice zp=a/3 buys a vanishing dipole term at the price of a changed radial variable and a changed quadrupole.
Application(Induced dipoles and polarizability)AdaptedAdapted from Riley, Hobson & Bence (2006), §21.3.1, p. 734
Riley, Hobson and Bence's conducting sphere is the cleanest instance of a multipole moment being produced rather than measured. An uncharged conducting sphere of radius a in a uniform field E must have a potential that tends to −Ercosθ far away and is constant on the surface; since cosθ=P1(cosθ), only the ℓ=1 terms of the separated solution can appear (separation of variables), and the two conditions determine them, giving v=−Er(1−a3/r3)cosθ. Comparing the second term with the point-dipole potential μP1(cosθ)/4πε0r2 identifies an induced dipole moment 4πε0a3E: the polarizability of a conducting sphere is 4πε0a3, three times its volume times ε0. The same calculation with a dielectric instead of a conductor gives the Clausius–Mossotti factor, and that is how the refractive index of a dilute suspension is computed.
ApplicationOriginal
Nuclear and molecular moments as experimental quantities. The ℓ=2 moment of a nucleus is measurable through its interaction with the electric field gradient of the surrounding electrons, and a nonzero value is direct evidence that the nucleus is not spherical — prolate for a positive quadrupole moment, oblate for a negative one. The origin-independence statement is what makes such a number well defined: a nucleus has zero dipole moment (by parity), so its quadrupole moment is its lowest nonvanishing one and does not depend on where the coordinate origin is put. By the same token a quoted molecular dipole moment is meaningful only for a neutral molecule, which is why one reads of the dipole moment of water but of the quadrupole moment of an ion.
Gravity: geodesy and gravitational waves. The expansion applies unchanged to Newtonian gravity, with ρ the mass density and the sign reversed. Because mass density is positive, the monopole moment never vanishes, so the gravitational dipole moment of a body is always origin-dependent — and the origin that makes it vanish is the centre of mass. This is why the gravitational field of the Earth is quoted as a series of Jℓ coefficients about the centre of mass, with the flattening appearing as J2≈1.08×10−3 and everything beyond it some three orders of magnitude smaller (J3≈−2.5×10−6, J4≈−1.6×10−6, each about 1/430 of J2); satellite orbits are fitted to those coefficients, and the measured drift of the coefficients tracks the melting of ice sheets. In general relativity the leading radiative term is the second time derivative of the mass quadrupole, as Arfken notes; the mass monopole and dipole cannot radiate because they are conserved.
Fast summation. The expansion is the mathematical content of the fast multipole method. To evaluate the potential of N sources at N targets, group the sources into cells; for a target far from a cell, replace the sources in the cell by the first few of its multipole moments, computed once for the whole cell. The error is controlled by the convergence rate in (1), which is governed by r</r> — the ratio of the cell radius to the distance — and the example above shows how fast that is in practice: truncating at the quadrupole gave three-digit accuracy at five cell radii. This turns an O(N2) computation into O(N), and it is the standard tool in astrophysical N-body simulation, molecular dynamics with Coulomb interactions, and boundary-element solvers built on the fundamental solution of Laplace's equation.
Application(Partial waves in scattering)AdaptedAdapted from Arfken, Weber & Harris (2013), §16.3, Eq. (16.61), pp. 798–799
The same factorization of a two-point function into products of one-point functions is what makes partial-wave analysis possible. Arfken's Eq. (16.61) is the companion expansion eik⋅r=4π∑ℓmiℓjℓ(kr)Yℓm(Ωk)∗Yℓm(Ωr), reached by expanding eikrcosχ in Legendre polynomials, recognising the coefficient integral as the integral representation of jℓ, and then applying the same addition theorem. Spherical Bessel functions take the place of the powers of r, and the identity converts an incident plane wave into a superposition of spherical waves, each of which is scattered independently by a central potential. See the Helmholtz equation and spherical Bessel functions.
Reduces the wave and diffusion equations to the Helmholtz equation by separating a time-harmonic factor, collects its elementary solutions - plane waves, cylindrical waves and the spherical waves built from spherical Bessel and Hankel functions - and states and proves the Sommerfeld radiation condition that picks out the outgoing solution and forces a field that is source-free on all of space to vanish.
Definition(Time-harmonic reduction and the Helmholtz equation)AdaptedAdapted from Riley, Hobson & Bence (2006), §21.3.2, Eqs. (21.52)–(21.54), p. 737; Hassani (2013), §19.7.2, Eq. (19.38), p. 593; Evans (2010), §6.6, Problem 16, p. 347
A field oscillating at a single angular frequency ω is time-harmonic: it has the form
u(x,t)=e−iωtw(x),(1)
the physical field being Reu. Substituting (1) into the wave equation ∇2u=c−2utt gives
∇2w+k2w=0,k=cω,(2)
the Helmholtz equation, with k the wavenumber. Equation (2) reduces to Laplace's equation when k=0. We also meet the modified Helmholtz equation∇2w−k2w=0, which is (2) with k2 replaced by −k2: its radial solutions grow or decay exponentially instead of oscillating, and its only bounded solution on all of Rn is w≡0. (It does have solutions that oscillate in some directions while decaying in others — cos(ax)e−k2+a2z is one — and those evanescent fields are what the waveguide and tunnelling paragraphs below are about; what it has not got is a bounded oscillation filling space.)
Sign and time conventions, which differ between our sources. Riley, Hobson and Bence separate time without choosing a sign, writing T(t)=Aeiωt+Be−iωt with ω=kc, and get (2) either way, since only ω2 enters. Evans fixes u=e−iσtw and writes the equation as −Δw=λw with λ=σ2, so his σ is our k in units with c=1. Arfken fixes the convention at the other end, by declaring an outgoing wave to be one whose asymptotic radial dependence is exp(+ikr). We use e−iωt throughout, which makes ei(kr−ωt) travel outwards and so agrees with Arfken: with the opposite time convention e+iωt, the outgoing solution would be e−ikr instead, and every sign in the radiation condition below would flip.
Hassani records why (2) is the equation to study: separating time out of the diffusion equation, the wave equation, or the Schrödinger equation of a free particle leaves (2) in each case, the only difference being the time factor one has discarded. For the wave equation it is e±iωt with ω=kc; for the diffusion equation κ∇2u=ut it is the decaying e−k2κt, so that there k2 is a decay rate rather than a frequency.
Strip the time. Write the field as e−iωtw(x) and reduce the equation to (2) with k=ω/c. Everything from here on is a problem for a complex-valued function of position; the physical field is recovered at the end as Re(e−iωtw), and amplitudes and phases are the modulus and argument of w.
Choose coordinates and separate. Apply the recipe of separation of variables. The special functions are those of the catalogue there, now with k=0: Jm(kρ),Ym(kρ) in plane polars, Jm(k2−α2ρ) in cylindrical polars, and the spherical Bessel functions jℓ(kr),yℓ(kr) in spherical polars.
Interior problem: discard what is singular, then quantize. If the region contains the origin, drop Ym and yℓ. A homogeneous boundary condition on a closed surface then admits a solution only for special values of k — the eigenvalues — and the frequencies ω=ck are the resonances of the cavity.
Exterior problem: impose a condition at infinity. The Helmholtz equation in an unbounded region does not have enough boundary data to be determined, because its solutions oscillate rather than decay and the usual "vanishes at infinity" is too weak to rule anything out. The missing condition is the Sommerfeld radiation condition, which in Evans's form reads
r→∞limr(∂r∂w−ikw)=0(n=3),(3)
and which Arfken states in the equivalent working form that the asymptotic radial dependence must be e+ikr: energy flows outwards from the source and none comes in from infinity. Statement (3) is what makes the exterior problem well posed; see the uniqueness theorem below.
For the modified Helmholtz equation, use decay instead. With ∇2−k2 the radial solutions are exponentially growing or decaying rather than oscillatory, and the right condition at infinity is simply that w→0 — Arfken's choice in his Table 10.1, where the three-dimensional modified-Helmholtz Green's function is −e−k∣r1−r2∣/4π∣r1−r2∣, the screened potential.
Build the solution from the elementary waves. Plane waves eik⋅x are the natural basis for a problem with planar symmetry, outgoing spherical waves hℓ(1)(kr)Yℓm for a radiating or scattering problem, and jℓ(kr)Yℓm for the field inside a sphere. The plane-wave expansion of spherical Bessel functions converts between the first two.
Plane waves. For a constant vector k, the function w(x)=eik⋅x solves (2) if and only if k⋅k=k2. Writing k=kω with ∣ω∣=1, the corresponding time-harmonic field is the travelling wave u=eik(ω⋅x−ct), moving with speed c in the direction ω.
Spherical waves. In R3 the functions
w±(x)=re±ikr,r=∣x∣,
solve (2) at every point with r>0; w+ is outgoing and w− incoming. More precisely, Evans records that
Φ(x)=4πreikrsatisfies−∇2Φ=k2Φ+δ0 in R3,
so Φ is the outgoing fundamental solution of the Helmholtz operator: it reduces to the fundamental solution 1/4πr of Laplace's equation as k→0. Arfken's Table 10.1 lists the same object as −eik∣r1−r2∣/4π∣r1−r2∣, the sign difference being his convention ∇2G+k2G=δ against Evans's −∇2Φ−k2Φ=δ.
The separated solutions in spherical polars. The radial equation is
R′′+r2R′+[k2−r2ℓ(ℓ+1)]R=0,
whose solutions are the spherical Bessel functions R=Ajℓ(kr)+Byℓ(kr), which Hassani builds out of the half-integer-order Bessel functions through the normalization zℓ(x)=π/2Zℓ+1/2(x)/x. So
with every Bℓm=0 if the region contains the origin, since yℓ is singular there. The outgoing combination is the spherical Hankel function hℓ(1)=jℓ+iyℓ, and h0(1)(kr)=−ieikr/kr recovers part 2.
Cylindrical waves. In cylindrical polars (ρ,φ,z) the separated solutions are
w=Zm(κρ)eimφeiαz,κ=k2−α2,m∈Z,
where Zm stands, as in Hassani's normalization above, for either of Jm and Ym — only Jm if the region contains the axis, since Ym is singular there — and α is the separation constant of the z direction. The radial argument κρ=k2−α2ρ is Riley, Hobson and Bence's. The three wavenumbers are tied by κ2+α2=k2, which is the dispersion relation of a waveguide. For α2<k2 the factor Zm(κρ) oscillates in ρ and the wave propagates along the axis; for α2>k2 the argument is imaginary and the radial factor is no longer oscillatory. In the plane, α=0 and the solutions are Jm(kρ)eimφ and Ym(kρ)eimφ.
Lower dimensions. In two dimensions the role of eikr/4πr is played by −4iH0(1)(kρ), with a Hankel function in place of the exponential; in one dimension by −2kieik∣x∣. Arfken tabulates all three. He also warns that the step from two dimensions to three is a genuine change and not a change of notation — a remark he makes on the previous page about the Laplace Green's functions, where 2π1ln∣ρ1−ρ2∣ grows at infinity while −1/4π∣r1−r2∣ decays, and which applies just as well to the entries here.
Idea of the proofOriginal
Parts 1 and 2 are two-line verifications, and the only thing worth watching is where each one fails. A plane wave is bounded but does not decay in any direction, so it is a legitimate solution of (2) on all of R3 with no source anywhere: the Helmholtz equation on all space has a huge supply of global solutions, which is exactly why a condition at infinity is needed to make an exterior problem well posed. A spherical wave does decay, like 1/r, but it has a singularity at the origin, and the singularity is the source. Part 2 says that the singularity is precisely a unit point source: near r=0 the factor eikr is 1+O(r), so Φ looks like the Laplace fundamental solution 1/4πr, whose Laplacian is −δ0, and the correction accounts for the term k2Φ.
Part 3 is the separation already carried out in spherical polars, read with k=0: nothing changes in the angular factors, and the radial factor becomes a Bessel function of half-integer order.
Step 1 (plane waves). For w=eik⋅x we have ∂w/∂xj=ikjw and hence ∂2w/∂xj2=−kj2w. Summing over j,
∇2w=−∣k∣2w,
so ∇2w+k2w=(k2−∣k∣2)w, which vanishes identically if and only if ∣k∣=k (as w is nowhere zero). Writing k=kω and restoring the time factor, u=e−iωteikω⋅x=eik(ω⋅x−ct) using ω=kc, which is constant on the planes ω⋅x−ct=const moving with speed c.
Step 2 (spherical waves, away from the origin). For a radial w=v(r), the Laplacian is v′′+2v′/r (the case n=3 of the formula in the fundamental solution), which can be written r1(rv)′′. With v=e±ikr/r we get rv=e±ikr and (rv)′′=−k2e±ikr, so
∇2v=r1(−k2e±ikr)=−k2v,
and ∇2v+k2v=0 for r>0.
Step 3 (the source at the origin). Put Φ=eikr/4πr and let ϕ be a test function. Write Φ=4πr1+Ψ with
Ψ=4πreikr−1,
which is bounded near the origin, since eikr−1=O(r); indeed Ψ→ik/4π as r→0. Then:
4πr1 is the Laplace fundamental solution, so −∇2(4πr1)=δ0 in the sense of distributions (the Dirac delta in calculations).
By Step 2, −∇2Φ=k2Φ pointwise for r>0; subtracting the previous item, −∇2Ψ=k2Φ for r>0. This also holds in the sense of distributions on all of R3, with no term supported at 0, and the reason is that Ψ is too mild a function to carry one. Indeed, work outside a small ball, where Ψ is smooth: Green's second identity on {r>ε} turns ∫Ψ∇2ϕ into ∫ϕ∇2Ψ plus a boundary term on r=ε, giving
−∫r>εΨ∇2ϕdV=k2∫r>εΦϕdV+∮r=ε(Ψϕr−ϕΨr)dS.
Expanding the exponential, Ψ=4πik−8πk2r+O(r2), so Ψ and Ψr are both bounded near the origin while the sphere r=ε has area 4πε2: the boundary term is O(ε2) and drops out as ε→0. Both volume integrals converge, the first because Ψ is bounded and the second because ∣Φ∣≤1/4πr is locally integrable. So −∇2Ψ=k2Φ as distributions, and the only singular contribution in Φ sits in 4πr1.
Adding the two statements gives −∇2Φ=k2Φ+δ0 in R3, which is Evans's part (b).
Step 4 (the radial equation and the spherical Bessel functions). This is Step 8 of the proof in separation of variables: the substitution R=r−1/2S turns the radial equation into Bessel's equation of order ℓ+21 in the variable kr. Hassani's normalization zℓ(x)=π/2Zℓ+1/2(x)/x is exactly the factor r−1/2 together with a constant chosen so that jℓ(x)∼xℓ/(2ℓ+1)!! at small x; with it, R=Ajℓ(kr)+Byℓ(kr). □
Method(The Sommerfeld radiation condition and uniqueness)AdaptedAdapted from Evans (2010), §6.6, Problems 16(c) and 17, pp. 347–348
Throughout, k is real and nonzero, w is complex-valued and defined at least outside some ball, and wr=∇w⋅x/∣x∣ is its radial derivative. The Sommerfeld radiation condition on w is
r→∞limr(wr−ikw)=0,(3)
the limit being uniform in the direction x/r: for every ε>0 there is an R0 with ∣r(wr−ikw)∣≤ε at every point with r≥R0. (This is Evans's reading, and it is what makes the surface estimate of Step 6 below work.) Then:
Of the three solutions of (2) found above — all of them solutions away from the origin, where the first two are singular — the outgoing fundamental solution Φ=eikr/4πr satisfies (3), while the incoming one e−ikr/4πr and the plane wave eikω⋅x both fail it.
Uniqueness. If w solves −∇2w=k2w on all of R3 and satisfies (3), then w≡0.
Part 2 is the statement that makes the condition useful: it says that the exterior Helmholtz problem — the equation outside an obstacle, data on the obstacle, and (3) at infinity — has at most one solution, because the difference of two solutions would solve the source-free problem everywhere and hence vanish. Evans sets both parts as problems, with the hints reproduced in the proof below, and refers to Colton for a fuller treatment.
Idea of the proofOriginal
Condition (3) is a statement about which way the energy goes. For a field behaving like A(θ,φ)eikr/r at large r, the radial derivative is ikAeikr/r plus a term of order r−2, so wr−ikw is O(r−2) and (3) holds; for Ae−ikr/r the two terms add instead of cancelling, leaving −2ikw, which is O(r−1), and r(wr−ikw) does not go to zero. So (3) selects e+ikr over e−ikr, which with the time convention e−iωt means outgoing over incoming. A plane wave fails for a different reason: its radial derivative is ik(ω⋅x^)w, so wr−ikw is ik(ω⋅x^−1)w, which does not decay at all in any direction other than ω — the plane wave is coming in from infinity on one side as much as it is going out on the other.
The uniqueness proof has two halves. First, Green's second identity applied to w and wˉ on a large ball gives a conserved flux: because wˉΔw−wΔwˉ=0 for a solution with real k2, the boundary integral ∫(wˉwr−wwˉr)dS vanishes for every R. That integral is 2iIm∫wˉwr, so the radiation condition — which controls wr−ikw — combines with it to force ∫∂B(0,R)(∣wr∣2+k2∣w∣2)dS→0: the field carries no energy at all through large spheres, and both w and wr die faster than 1/r in the mean-square sense. Second, Green's identity against the fundamental solution represents w at any point as an integral over ∂B(0,R); the first half shows that integral tends to zero, so w vanishes at that point.
ProofAdaptedAdapted from Evans (2010), §6.6, Problems 16(c) and 17, pp. 347–348
Step 1 (Φ satisfies (3), exactly). With Φ=eikr/4πr,
Φr=drd4πreikr=4πrikeikr−4πr2eikr=ikΦ−rΦ,
so r(Φr−ikΦ)=−Φ=−eikr/4πr, whose modulus is 1/4πr→0. Condition (3) holds, with the quantity in it of size exactly 1/4πr.
Step 2 (the incoming wave fails). With Ψ=e−ikr/4πr the same computation gives Ψr=−ikΨ−Ψ/r, so
r(Ψr−ikΨ)=r(−2ikΨ−rΨ)=−2ikrΨ−Ψ,
of modulus −2ikrΨ−Ψ→k/2π as r→∞, since ∣rΨ∣=1/4π and Ψ→0. The limit is nonzero because k=0.
Step 3 (the plane wave fails). For w=eikω⋅x, ∇w=ikωw and so wr=ik(ω⋅x^)w with x^=x/r. Hence
r(wr−ikw)=ikr(ω⋅x^−1)w,
and ∣w∣=1. Along any ray with ω⋅x^=1 the modulus is kr∣ω⋅x^−1∣, which grows without bound. So (3) fails, and badly.
Step 4 (a vanishing flux). Now let w solve −∇2w=k2w on R3 and satisfy (3). Then Δw=−k2w and, taking complex conjugates and using that k2 is real, Δwˉ=−k2wˉ. Therefore
Step 6 (both sides tend to zero). By (3), whose limit is uniform in direction, there is for each ε>0 an R0 with ∣wr−ikw∣≤ε/r at every point with r≥R0. The sphere ∂B(0,R) has area 4πR2, so the left side of (5) is at most 4πR2⋅ε2/R2=4πε2 for R≥R0. Hence
∮∂B(0,R)(∣wr∣2+k2∣w∣2)dS⟶0(R→∞).(6)
Both terms are non-negative, so each tends to zero separately; and since k=0, (6) controls ∣w∣ as well as ∣wr∣.
Step 7 (represent w on a large sphere). Fix x0 and let R>∣x0∣. Apply Green's second identity on B(0,R) to the pair w and Φ(⋅−x0), excising a small ball about x0 exactly as in the proof of the representation formula for Laplace's equation. The volume terms cancel, because both functions satisfy the same equation Δv=−k2v away from x0; the excised sphere contributes w(x0) in the limit, because Φ has the same singularity as the Laplace fundamental solution (Step 3 of the previous proof). What remains is Evans's representation
w(x0)=∮∂B(0,R)(Φwr−wΦr)dS,Φ=Φ(x−x0).(7)
Step 8 (let R→∞). On ∂B(0,R) with R>2∣x0∣ we have ∣x−x0∣≥R/2, so ∣Φ∣≤1/(2πR) and, from the identity ∇Φ=(ik−1/ρ)Φρ^ of Step 1 with ρ=∣x−x0∣, also ∣Φr∣≤∣∇Φ∣≤k∣Φ∣+∣Φ∣/ρ≤k/(2πR)+1/(πR2)≤C/R for large R. By the Cauchy–Schwarz inequality on the sphere,
whose first factor is bounded and whose second tends to 0 by (6). The same estimate applies to ∮wΦrdS, using the ∣w∣2 half of (6). Hence the whole right-hand side of (7) vanishes in the limit, while its left-hand side is a fixed number independent of R. Hence w(x0)=0, and since x0 was arbitrary, w≡0. □
RemarkOriginal
Where the hypotheses are used, and what fails without them.
k real and nonzero, in Steps 2, 4 and 6. Step 4 needs k2=k2, which is where "real" is used; a complex k (an absorbing medium) makes the flux identity (4) false, and then the right condition at infinity is plain decay, as for the modified Helmholtz equation. Step 2 needs k=0 to make the incoming wave fail, and Step 6 needs it to extract ∣w∣ from (6). At k=0 the statement collapses: the outgoing and incoming fundamental solutions e+ikr/4πr and e−ikr/4πr both become the single function 1/4πr, so there is nothing left to choose between, and the condition at infinity reverts to w→0.
Three dimensions, in Steps 6 and 8. The factor r in (3) and the area 4πR2 are matched to n=3. In n dimensions the condition carries r(n−1)/2 instead, matching the r−(n−1)/2 decay of an outgoing wave; in two dimensions the outgoing solution is H0(1)(kρ), which decays like ρ−1/2, and the condition reads limρ1/2(wρ−ikw)=0.
The condition is not implied by decay. This is the whole point. Both eikr/r and e−ikr/r tend to zero at infinity, and so does any combination of them; "w→0" therefore fails to distinguish a radiating source from one absorbing energy from infinity, and fails to make the exterior problem uniquely solvable. Nor is boundedness enough, as the plane wave shows: it is bounded, source-free and not identically zero, so part 2 of the statement would be false with (3) weakened to boundedness.
What part 2 does and does not say. It is a statement about solutions on all of R3 with no source. It does not say that an exterior problem has no solution — only that it has at most one, which is what one wants. The existence of that one solution is proved by constructing it from the outgoing Green's function, which is the subject of Green's functions for the Helmholtz equation.
Interior resonances are the other side of the same coin. For a bounded region with homogeneous boundary data, (2) has nonzero solutions for a discrete set of k and only for those; for an exterior region with the radiation condition, it has none. That contrast is why cavity problems are eigenvalue problems and radiation problems are not.
A separated solution has a frequency but not a wavelength. Riley, Hobson and Bence make this observation about (2) in spherical polars: the time factor e−iωt with ω=kc is a pure frequency, but the zeros of jℓ(kr) are not equally spaced in r except for ℓ=0, so the radial profile has no single wavelength. It acquires one only asymptotically, where jℓ(kr)≈sin(kr−ℓπ/2)/kr.
Take k=2 and compare three solutions of ∇2w+4w=0 in R3 (away from the origin where necessary), computing Q(r)=sup∣x∣=r∣r(wr−2iw)∣ exactly.
w
wr
r(wr−ikw)
∣⋅∣ at r=10
r→∞
4πre2ir
2iw−rw
−w
40π1=0.0079577
→0
4πre−2ir
−2iw−rw
−4irw−w
4πr16r2+1r=10=0.3184093
→2πk=0.3183099
e2iz (i.e. ω=z^)
2icosθw
2ir(cosθ−1)w
0 at θ=0, 40 at θ=π
unbounded
The first row is Step 1 of the proof with the numbers put in: the quantity in (3) has modulus exactly 1/4πr, so it decays like 1/r and the limit is zero. The second row is Step 2: the modulus is 16r2+1/4πr, which tends to k/2π=1/π=0.3183099 and at r=10 is already 0.3184093, within 0.04% of its limit. The third row is Step 3: along the −z axis, cosθ=−1 and the modulus is 4r, which at r=10 is 40 and grows without bound; only along the single ray θ=0 does it vanish.
A superposition makes the selectivity vivid. The standing wave
w=2rsin2r=4i1(re2ir−re−2ir)=j0(2r)
is a solution of ∇2w+4w=0 on all of R3, with no singularity at the origin (where w→1), and it decays like 1/r. It is a counterexample to any uniqueness statement based on decay alone; and indeed it is half incoming and half outgoing, so it fails (3). Explicitly, combining the first two rows of the table linearly,
r(wr−ikw)=e−ikr−krsinkr,
whose modulus tends to 1 rather than to 0 — numerically 0.98226 at r=10, 1.00437 at r=50 and 1.00213 at r=100. Part 2 of the statement is consistent with the existence of this solution precisely because the solution fails (3).
Example(Resonances of a spherical shell cavity)Original
Find the values of k for which ∇2w+k2w=0 has a nonzero solution in the shell 1≤r≤2 with w=0 on both bounding spheres.
Separate. The region excludes the origin, so both spherical Bessel functions survive and the radial factor for each ℓ is R(r)=Ajℓ(kr)+Byℓ(kr). The two boundary conditions give a homogeneous 2×2 system for (A,B), which has a nonzero solution exactly when
Dℓ(k)=jℓ(k)yℓ(2k)−yℓ(k)jℓ(2k)=0.(8)
The case ℓ=0 in closed form. Here j0(x)=sinx/x and y0(x)=−cosx/x, so
So D0(k)=0 exactly when k=nπ, n=1,2,…. This is no accident: for ℓ=0 the substitution R=χ/r turns the radial equation into χ′′+k2χ=0 with χ(1)=χ(2)=0, which is a string of length 1 — the shell behaves, in its spherically symmetric modes, exactly like a one-dimensional box of width 1, the thickness of the shell.
Higher ℓ, numerically. Solving (8) for the first few roots:
ℓ
kℓ,1
kℓ,2
kℓ,3
degeneracy 2ℓ+1
0
3.141593
6.283185
9.424778
1
1
3.286007
6.360678
9.477196
3
2
3.555788
6.513068
9.581270
5
3
3.922520
6.735565
9.735534
7
The ℓ=0 roots are π,2π,3π to six figures, confirming the closed form. The lowest resonance of the cavity is k=π, simple; the next is k=3.286007, threefold degenerate. Reading down any column, the roots rise with ℓ, as they must: a larger ℓ adds the repulsive term ℓ(ℓ+1)/r2 to the radial equation, which pushes the eigenvalue up. Reading along a row instead — successive roots at fixed ℓ — the spacing tends to π, the ℓ=0 value, because at large k the centrifugal term becomes negligible beside k2: at ℓ=3 the two gaps printed are 2.813045 and 2.999969.
Example(A pulsating sphere: the radiation condition doing work)Original
A sphere of radius a in air oscillates radially at angular frequency ω, its surface moving with radial velocity amplitude U0. Find the radiated sound field and the acoustic power.
Set up. The acoustic pressure satisfies the wave equation, so p=Re(e−iωtw) with ∇2w+k2w=0, k=ω/c, outside the sphere. The problem is spherically symmetric, so only ℓ=0 appears; of the two ℓ=0 solutions e±ikr/r, the radiation condition keeps only the outgoing one. This is the step that makes the problem solvable at all: without (3), any combination of the two would do, and the answer would not be determined by the data on r=a. So
w(r)=Areikr.
Apply the surface condition. Linearized momentum balance, ρ0∂v/∂t=−∇p, gives for the time-harmonic parts −iωρ0vr=−w′(r), so vr=w′(r)/(iωρ0). Using w′=ikw−w/r and imposing vr=U0 at r=a,
So the pressure amplitude at 1 m is 0.0188 Pa — root-mean-square 0.01327 Pa, which against the reference 2×10−5 Pa is 56.4 dB SPL — and the sphere radiates about 5.4 microwatts.
The two limits, and why small loudspeakers are inefficient. For ka≫1 the formula gives P→2πρ0ca2U02=(4πa2)⋅21ρ0cU02: the whole surface radiates as if it were a piece of plane wave, at the best possible rate. For ka≪1,
P≈2πρ0ck2a4U02,
smaller by the factor k2a2/(1+k2a2), which here is 0.00832. A source much smaller than a wavelength radiates with efficiency proportional to (ka)2, and since k∝f the power at fixed surface velocity falls off as f2 — which is why a small driver needs a large enclosure, a horn, or a great deal of excursion to produce bass.
ApplicationOriginal
Acoustics, antennas and scattering. Every steady-state radiation calculation is an exterior Helmholtz problem closed by the radiation condition: a loudspeaker in a baffle, a dipole antenna, an ultrasound transducer. For scattering, one splits the field into a prescribed incident part and an unknown scattered part, w=winc+wsc, and imposes (3) on wsc alone — not on w, which contains the incident plane wave and therefore cannot satisfy it, as Step 3 of the proof shows. That split is exactly what the uniqueness theorem makes legitimate: the scattered field is the unique outgoing solution with the boundary data that the incident field forces on the obstacle. The far-field pattern A(θ,φ) in wsc∼A(θ,φ)eikr/r is the scattering amplitude, and ∣A∣2 integrated over directions is the scattering cross-section; see scattering and the Born approximation.
Waveguides and optical fibres. In a cylinder, step 2 of the recipe with the catalogue of separation of variables gives modes Jm(κρ)eimφeiβz with κ2+β2=k2. A mode propagates only when β2>0, that is when k exceeds the cutoff κ fixed by the boundary condition; below cutoff, β is imaginary and the mode is evanescent. The relation κ2+β2=k2 is the algebraic relation among separation constants of step 4 of that recipe, and here it is the dispersion relation of the guide — which is where waveguide dispersion, and hence pulse spreading in a fibre, comes from.
Cavity resonances and the density of states. The interior problem of step 3 gives the resonant frequencies of a cavity: a microwave oven, an organ pipe, a laser cavity. Counting eigenvalues below k in a box of volume V gives Weyl's asymptotic N(k)≈Vk3/6π2, which multiplied by two polarizations is the mode count behind the Rayleigh–Jeans and then the Planck radiation law. The shell example above shows the structure in miniature: a discrete set of k, with multiplicity 2ℓ+1 forced by the rotational symmetry.
Quantum mechanics of a free particle, and the modified equation in a barrier. The time-independent Schrödinger equation −2mℏ2∇2ψ=Eψ is the Helmholtz equation with k2=2mE/ℏ2; Hassani reaches it exactly this way, and the solutions regular at the origin are jℓ(kr)Yℓm. Where the potential exceeds the energy, k2 becomes negative and the equation turns into the modified Helmholtz equation, whose solutions decay exponentially — quantum tunnelling, in the form in which it is actually computed. The same substitution turns the Coulomb-screening problem into the modified equation, whose three-dimensional Green's function Arfken tabulates as the Yukawa potential e−kr/4πr. See the Schrödinger equation.
Numerical methods, and why the radiation condition is a nuisance. A computational domain is finite, so the condition at infinity has to be faked: by a non-reflecting or absorbing boundary condition, which enforces (3) approximately on an artificial outer boundary; by a perfectly matched layer, which continues the solution into complex coordinates so that outgoing waves decay and incoming ones would grow; or by reformulating the problem as a boundary integral equation on the obstacle using the outgoing Green's function, which satisfies (3) exactly and so needs no outer boundary at all. The last route is the fully faithful one, and it is why the boundary-element method dominates exterior acoustics and electromagnetic scattering. It also inherits the one difficulty of the Helmholtz equation that the Laplace equation does not have: the interior resonances of step 3 appear as frequencies at which the integral operator becomes singular, even though the exterior problem itself remains uniquely solvable.
States the time-dependent and time-independent Schrodinger equations, shows that separating the time turns the first into an eigenvalue problem whose solutions are the stationary states, works out the free particle, the infinite well and the harmonic oscillator with their normalizations and energy levels, and derives the probability current and the continuity equation that makes the normalization persist.
A non-relativistic particle of mass m moving in a potential V(r) is described by a complex-valued wave functionΨ(r,t) satisfying the time-dependent Schrodinger equation
−2mℏ2∇2Ψ+V(r)Ψ=iℏ∂t∂Ψ,(1)
where ℏ=h/2π and h is Planck's constant. The quantity ∣Ψ(r,t)∣2 is the probability density for finding the particle at r at time t, so a physically admissible Ψ must satisfy the normalization
∫∣Ψ(r,t)∣2dV=1.(2)
The time-independent Schrodinger equation is the eigenvalue problem
−2mℏ2∇2ψ+V(r)ψ=Eψ,that isH^ψ=Eψ,(3)
with H^=−2mℏ2∇2+V the Hamiltonian operator. Kreyszig notes at once what (3) is for, and he puts the claim carefully: the operator form suggests that which energies are possible depends on the spectrum of H^, and a solution counts as physically meaningful only if it stays finite and tends to zero at infinity — it is that restriction, not the differential equation alone, that quantizes the energy. The identification of the admissible energies with the spectrum of H^ is the stronger statement, and it is Reed and Simon's rather than his: it presupposes a domain on which H^ is self-adjoint, and the remark below returns to it. Kreyszig's §11.4 records how H^ is written down: take the classical Hamiltonian H=Ekin+V in terms of positions and momenta, and replace each momentum pj by the operator 2πih∂/∂qj=−iℏ∂/∂qj (his §11.2, Eq. (2)) and each coordinate qj by multiplication by qj. Note the sign: 1/i=−i, so Kreyszig's 2πih∂/∂qj, written without a leading minus, is the physicists' −iℏ∂/∂qj.
Units and conventions. Kreyszig writes h where we write 2πℏ, so his equation (5) reads Δψ+h28π2m(E−V)ψ=0, which is (3). Hassani and Riley, Hobson and Bence use ℏ as above. Reed and Simon work in units with ℏ=1 and write the dynamics as dtd[U(t)φ]=iH[U(t)φ], which is (1) with the opposite sign of i — their U(t)=eitH against the physicists' e−iH^t/ℏ. We follow the physicists' convention (1), so that a stationary state carries the time factor e−iEt/ℏ, as in Riley, Hobson and Bence. Nothing physical depends on the choice, since the two are exchanged by complex conjugation, but the sign of the probability current below does.
Separable solutions of (1) of the form Ψ(r,t)=ψ(r)T(t) exist exactly when ψ solves (3) for some constant E∈C, and then
Ψ(r,t)=ψ(r)e−iEt/ℏ.(4)
If in addition E is real, such a Ψ is a stationary state: its probability density ∣Ψ∣2=∣ψ∣2 does not depend on t, because the time factor then has modulus one.
That hypothesis has to be carried, because (3) as written is a differential equation with no normalization and no boundary condition attached, and so does not by itself force E to be real: for complex E the time factor has modulus e(ImE)t/ℏ and the density grows or decays, and the remark below exhibits such a separable solution for the free particle. Reality is automatic at the one place it matters. If ψ is an eigenfunction of H^as a self-adjoint operator — in particular if V is real and ψ is normalizable — then E is real, as Step 2 of the proof shows in a line.
If the eigenfunctions ψn of (3) with eigenvalues En form a complete orthonormal set, the general solution of (1) with initial datum Ψ(r,0)=f(r) is
Ψ(r,t)=n∑cnψn(r)e−iEnt/ℏ,cn=∫ψnfdV,(5)
and ∑n∣cn∣2=1 when f is normalized. A superposition of two or more stationary states with different energies is not stationary: its density oscillates at the beat frequencies (En−En′)/ℏ.
Kreyszig arrives at (3) along a different road, which is worth recording because it explains the equation's shape. He starts from the wave equation for a scalar field, imposes the time-harmonic factor to get the Helmholtz equation Δψ+k2ψ=0 with k=2π/λ, and then substitutes de Broglie's λ=h/mv for the wavelength of matter waves; since mv2/2=E−V, the wavenumber becomes k2=8π2m(E−V)/h2, and the Helmholtz equation becomes (3). So the time-independent Schrodinger equation is the Helmholtz equation with a position-dependent wavenumber. Riley, Hobson and Bence make the same comparison for hydrogen, where the coefficient involves 1/r, and record that what one gets there resembles the Helmholtz equation without being it.
Idea of the proofOriginal
Equation (1) is first order in t and second order in space, like the heat equation, but with an i in front of the time derivative — and that single factor of i changes everything. Separating variables produces T′/T=−iE/ℏ instead of the heat equation's real decay rate, so for real E the time factor has modulus one and oscillates instead of decaying. The result is that a Schrodinger eigenmode persists for ever rather than dying out, and an initial state does not relax to anything; it merely has its modes rotate at different rates, which is why quantum evolution is reversible while diffusion is not.
The separation itself is the ordinary argument: divide (1) by ψT, note that one side depends on r alone and the other on t alone, and apply the separation lemma. The content is in which values of the separation constant are admissible, and that is decided by the requirement that ψ be normalizable — the three examples below show it deciding differently in three cases.
Step 1 (separate). Put Ψ=ψ(r)T(t) into (1) and divide by ψT at the points where neither factor vanishes — ψ may well have interior nodes, as the well's own ψ2 does at x=L/2:
ψ1[−2mℏ2∇2ψ+Vψ]=iℏT(t)T′(t).
The left side depends only on r and the right only on t, so by the separation lemma of separation of variables both equal a constant, which we call E. This gives −2mℏ2∇2ψ+Vψ=Eψ and iℏT′=ET wherever the division was legitimate; both sides of each identity are continuous, so each extends to the whole domain, and (3) holds there.
Step 2 (the time factor). The equation T′=−iET/ℏ has the general solution T=Ae−iEt/ℏ, for any E∈C; absorbing A into ψ gives (4). Now impose the hypothesis that E is real. Then ∣T∣=1, so ∣Ψ∣2=∣ψ∣2 does not depend on t and Ψ is stationary.
Reality is a genuine extra hypothesis and not a consequence of (3) — the remark below exhibits a separable solution with E purely imaginary — but it is automatic in the case the rest of the unit works in. If V is real and ψ is a normalizable eigenfunction of H^, then, with the inner product linear in its first slot and H^ symmetric on its domain,
E∥ψ∥2=(H^ψ,ψ)=(ψ,H^ψ)=E∥ψ∥2,
and ∥ψ∥=0, so E=E. Finding a domain on which H^ really is self-adjoint, rather than merely symmetric, is Reed and Simon's first mathematical problem.
Step 3 (superposition). Each term of (5) solves (1), by Step 2, and (1) is linear, so any finite sum does; for an infinite sum, the convergence of ∑∣cn∣2 and the completeness of {ψn} give convergence in the mean-square sense (eigenfunction expansions). At t=0 the sum is ∑cnψn, which is the expansion of f, so the initial condition holds; the coefficients are cn=(f,ψn) by orthonormality, which is the integral in (5) with our convention that an inner product is linear in its first slot.
Step 4 (normalization is preserved, and a superposition is not stationary). By orthonormality,
in which every off-diagonal pair with En=En′ contributes a term oscillating at the frequency ∣En−En′∣/ℏ. Only if a single energy is present does the density stand still. □
Application(The probability current and the continuity equation)AdaptedAdapted from Hassani (2013), §14.8.1, Eqs. (14.49)–(14.52), pp. 452–453
For a solution Ψ of (1) with real V, define the probability density and probability current density
ρ=∣Ψ∣2,J=mℏIm(Ψ∇Ψ)=2miℏ(Ψ∇Ψ−Ψ∇Ψ).(6)
Then
∂t∂ρ+∇⋅J=0.(7)
Consequently dtd∫ΩρdV=−∮∂ΩJ⋅n^dS for any fixed region Ω, and if Ψ and ∇Ψ decay suitably at infinity then ∫ρdV is constant in time: normalization, once imposed, is permanent.
Hassani obtains (6) and (7) in polar form. Writing Ψ=Aexp(iS/ℏ) with A and S real, equation (1) is exactly equivalent to the pair
and with ρ=A2 and J=A2∇S/m=ρv, where v=∇S/m, the second of these is (7). (The italicized claim is Hassani's (Hassani (2013), §14.8.1, after Eq. (14.50), p. 452), who calls the pair "completely equivalent" to the original equation, and we report it as he makes it. It should be read on a region where Ψ has no zero: at a node A=0, so the phase S is undefined there and the first equation divides by A; and in any case the polar form fixes S only modulo 2πℏ.) His Proposition 14.8.2 then reads the first equation in the limit ℏ→0: the right side disappears, S satisfies the classical Hamilton–Jacobi equation, and the pair describes a fluid of non-interacting classical particles of mass m in the potential V, with density ρ and current J.
Idea of the proofOriginal
Equation (7) is Green's second identity in disguise. Multiply (1) by Ψ, multiply the conjugate equation by Ψ, and subtract: the potential term cancels, because V is real, and what is left is iℏ∂t∣Ψ∣2 on one side and −2mℏ2(Ψ∇2Ψ−Ψ∇2Ψ) on the other. The bracket is exactly the integrand of Green's second identity, and the product rule turns it into the divergence of Ψ∇Ψ−Ψ∇Ψ. So the cancellation of V is where reality of the potential is used, and the identity Ψ∇2Ψ−Ψ∇2Ψ=∇⋅(Ψ∇Ψ−Ψ∇Ψ) is where the second-order operator becomes a divergence. A complex potential — an absorbing medium, or a decaying state — breaks the first step and probability is then not conserved, which is exactly what one wants it to mean.
ProofAdaptedAdapted from Hassani (2013), §14.8.1, pp. 452–453
Step 1 (two equations). Write (1) as
iℏ∂tΨ=−2mℏ2∇2Ψ+VΨ,
and take complex conjugates, using V=V:
−iℏ∂tΨ=−2mℏ2∇2Ψ+VΨ.
Step 2 (combine). Multiply the first by Ψ and the second by Ψ and subtract the second from the first. The terms V∣Ψ∣2 cancel — this is where the reality of V enters — leaving
iℏ(Ψ∂tΨ+Ψ∂tΨ)=−2mℏ2(Ψ∇2Ψ−Ψ∇2Ψ).
The left bracket is ∂t(ΨΨ)=∂tρ.
Step 3 (make the right side a divergence). By the product rule for the divergence,
∇⋅(Ψ∇Ψ)=∇Ψ⋅∇Ψ+Ψ∇2Ψ,
and the same with Ψ and Ψ exchanged. Subtracting, the symmetric first terms cancel:
Ψ∇2Ψ−Ψ∇2Ψ=∇⋅(Ψ∇Ψ−Ψ∇Ψ).
Step 4 (assemble). Substituting into Step 2 and dividing by iℏ,
∂tρ=−2miℏ∇⋅(Ψ∇Ψ−Ψ∇Ψ)=−∇⋅J,
which is (7), with J as in (6). The two expressions for J in (6) agree because Ψ∇Ψ=Ψ∇Ψ, so their difference is 2iIm(Ψ∇Ψ).
Step 5 (global conservation). Integrating (7) over a fixed region Ω with piecewise smooth boundary and applying the divergence theorem,
dtd∫ΩρdV=−∫Ω∇⋅JdV=−∮∂ΩJ⋅n^dS.
Taking Ω to be a large ball and letting its radius grow, the surface term vanishes provided ∣Ψ∇Ψ∣=o(r−2), and ∫ρdV is constant. □
RemarkOriginal
Where the hypotheses are used, and what fails without them.
Reality of V, in Step 2. Without it the potential terms do not cancel and (7) acquires a source ℏ2Im(V)ρ: probability leaks away where ImV<0. Optical-potential models of nuclear reactions use precisely this to represent absorption.
Reality of the separation constant E, in Step 2 of the separation proof. The constancy of ∣Ψ∣2 for a separable solution is a consequence of self-adjointness and not of the differential equation, and (3) on its own does not deliver it. Take V≡0 in one dimension and E=iε with ε>0 real. Then ψ=eαx with α2=−2miε/ℏ2, that is α=e−iπ/42mε/ℏ, solves (3), and
Ψ=eαxe−iEt/ℏ=eαxeεt/ℏ
solves (1) with density ∣Ψ∣2=e2(Reα)xe2εt/ℏ, which grows without bound in t. It is a separable solution and it is not stationary. What excludes it is not the equation but the normalization (2): Reα=mε/ℏ2>0, so ψ∈/L2(R) and ψ is not an eigenfunction of the self-adjoint H^ at all. Where a suitable domain is to be found is the first of the three problems recorded in the next remark.
Normalizability as the quantization condition. Equation (3) is a differential equation with a solution for every E; what rules out most E is the demand that the solution be in L2, or satisfy a boundary condition. Kreyszig makes the criterion explicit — a solution counts as physically meaningful only if it stays bounded everywhere and dies away at large distance — and observes that, once a potential field is fixed, (3) possesses such solutions for certain values of E and no others. The three examples below show the three outcomes: a continuous spectrum with no normalizable eigenfunctions at all (free particle), a discrete spectrum from a boundary condition (infinite well), and a discrete spectrum from decay at infinity (oscillator).
Completeness, in Step 3 of the separation proof. Formula (5) presumes that the eigenfunctions span the space of initial data. For the well and the oscillator this is the completeness of a Fourier sine series and of the Hermite functions respectively; for a potential with both bound and scattering states, the sum must be supplemented by an integral over the continuous spectrum, and (5) as written is incomplete.
The separation of angular variables survives any central potential. Riley, Hobson and Bence's general remark applies verbatim to (3): if V=V(r), then θ and φ enter only through ∇2, so the angular factors are spherical harmonics Yℓm and only the radial equation feels the potential. That is why hydrogen, the three-dimensional oscillator and a free particle share their angular wavefunctions, and why the orbital labels s,p,d,f mean the same thing in every atom.
What the time-dependent equation is not. It is not a wave equation: it is first order in time, so one initial condition determines the evolution, not two, and there is no finite propagation speed. It is not the heat equation either, despite the shape, because the factor i makes the evolution unitary rather than dissipative: ∣cn∣ is constant in (5), where a heat equation would have e−Ent and lose the high modes.
Remark(What the operator theory adds: three problems the formal equation leaves open)AdaptedAdapted from Reed & Simon (1980), §VIII.11, pp. 302–304
Reed and Simon set the model up in Hilbert-space terms: a separable Hilbert space, physical states as unit vectors determined up to a phase, observables as self-adjoint operators whose measurement probabilities come from a projection-valued measure, and the dynamics as a continuous one-parameter group of unitary operators U(t) whose self-adjoint generator H is the Hamiltonian. Written out for φ in the domain of H, their dynamical law is dtd[U(t)φ]=iH[U(t)φ], which is (1) in their units ℏ=1 and with their sign of i. Three mathematical problems are then left open, and they are the three this unit keeps running into.
Self-adjointness. The physics supplies a formal differential expression for H^ and no domain. One must find a domain on which it is essentially self-adjoint; where there is none, its self-adjoint extensions have to be surveyed and one of them chosen as the observable. The free-particle solution with imaginary E above is the elementary face of this: the differential equation admits a solution that the operator does not.
Spectral analysis. The point spectrum of H is where the stationary states sit, and a system disturbed from one of them moves to another, radiating at a frequency set by the gap between the two energies. This is the sense in which (3) is an experiment.
Scattering theory. The behaviour of U(t)φ for large ∣t∣.
Their formulation also makes the reversibility of (1) structural rather than observed: a one-parameter group of unitary operators has an inverse at every t, and the norm it preserves is exactly the normalization (2).
Example(The free particle: a continuous spectrum and no normalizable eigenstates)Original
Take V≡0 in one dimension. Then (3) is −2mℏ2ψ′′=Eψ, whose solutions for E>0 are ψ=e±ikx with
E=2mℏ2k2,
and for E<0 are the real exponentials e±κx with E=−ℏ2κ2/2m.
No eigenfunction is normalizable.∫−∞∞∣eikx∣2dx=∫1dx=∞, and eκx blows up at one end. So (3) has no solution satisfying (2) for any E: the free particle has no bound states, and its spectrum is the continuous half-line E≥0. Kreyszig's criterion — finite and tending to zero at infinity — excludes every E, and the resolution is that the physical states are superpositions,
Ψ(x,t)=2π1∫−∞∞f^(k)ei(kx−ℏk2t/2m)dk,
a wave packet, normalizable when f^∈L2. The exponent shows the dispersion relation ω=ℏk2/2m. Unlike the wave equation, where ω=ck and every component of the packet travels at the same speed, here the group velocity
dkdω=mℏk=mp
— the classical velocity of a particle of momentum p=ℏk — depends on k, and the packet therefore spreads: its envelope is broadened at the rate set by ω′′=ℏ/m=0, and a packet of initial width Δx0 doubles in width after a time of order m(Δx0)2/ℏ. It is ω′′=0 that does this and not the k-dependence of the phase velocity ω/k=ℏk/2m, which is exactly half the group velocity here: a medium with ω=ck+a has phase velocity c+a/k, varying with k, and spreads nothing at all, because every component of its envelope still moves at c.
The current of a plane wave. For Ψ=Aei(kx−ωt), formula (6) gives
J=mℏIm(Ae−ikx⋅ikAeikx)=mℏk∣A∣2=vρ,
with v=ℏk/m=p/m the classical velocity. So the current is density times velocity, exactly as in fluid mechanics, and the sign of k decides the direction of flow. Here ρ and J are both constant, so (7) holds trivially, 0=0; the example is still worth doing because it fixes the interpretation of J and the normalization of (6).
Example(The infinite well: energies, a beating superposition, and the current)Original
Take V=0 for 0<x<L and V=+∞ outside, which forces ψ(0)=ψ(L)=0.
The stationary states. Solving (3) as in the previous example and imposing the two boundary conditions gives ψ=sin(kx) with kL=nπ, so
ψn(x)=L2sinLnπx,En=2mL2n2π2ℏ2(n=1,2,…),
the normalization constant coming from ∫0Lsin2(nπx/L)dx=L/2. Here the quantization comes from the boundary condition rather than from decay at infinity, and the eigenfunctions are an orthonormal basis of L2(0,L) — the Fourier sine basis — so (5) applies with no continuous part.
Numbers. For an electron (m=9.1093837×10−31 kg) in a well of width L=1 nm, with ℏ=1.0545718×10−34 J s,
E1=2mL2π2ℏ2=6.024667×10−20J=0.376030eV,
and En=n2E1, so the first gap is E2−E1=3E1=1.128090 eV. A photon of that energy has wavelength hc/(E2−E1)=(1239.842eVnm)/(1.128090eV)=1099.06 nm, which is the near infrared, not the visible. The gap falls as L−2, so a narrower well is needed to reach visible light: the band from 1.77 eV (700 nm) to 3.10 eV (400 nm) corresponds to L between 0.798 nm and 0.603 nm. That L−2 dependence of the level spacing on size is the mechanism behind the colour of a quantum dot, which is set by its diameter — though a real dot is a well in a semiconductor, so one must add the bulk band gap to the confinement energy and use the carrier's effective mass in place of m.
A superposition, and its beat. Take the normalized initial state f=21(ψ1+ψ2), so c1=c2=1/2 in (5):
Ψ(x,t)=21[ψ1(x)e−iE1t/ℏ+ψ2(x)e−iE2t/ℏ].
Its density is ρ=21[ψ12+ψ22+2ψ1ψ2cosω21t] with
ω21=ℏE2−E1=2mL23π2ℏ=1.713871×1015s−1,
a period of 2π/ω21=3.666078 fs for the electron above. The mean position oscillates: since ∫0Lxψn2dx=L/2 for every n and
The particle sloshes from side to side with amplitude 0.18L — a classical-looking motion built from two stationary states, and the clearest demonstration that a superposition of different energies is not stationary.
The current, and an exact check of the continuity equation. The probability of being in the left half is
which is exactly dtd∫0L/2ρdx=−J(L/2,t), the one-dimensional case of Step 5 — the boundary x=0 contributing nothing because ψn(0)=0. The two sides agree identically in t, not just numerically, so this is a complete check of (6) and (7) on a nontrivial solution. Note also that a single stationary state here has J≡0, because ψn is real and Im(ψnψn′)=0: nothing flows in a standing wave.
Example(The harmonic oscillator: zero-point energy from the Hermite equation)AdaptedAdapted from Kreyszig (1978), §11.3, Example 11.3-1, Eqs. (7)–(12), pp. 585–588; Kreyszig (1978), §11.3, Prob. 5, p. 588; Kreyszig (1978), §11.3, Probs. 9–10, pp. 589–590
Take V=21mω02x2 in one dimension, Kreyszig's Example 11.3-1 in our notation. Equation (3) is
−2mℏ2ψ′′+21mω02x2ψ=Eψ.(8)
Non-dimensionalize. Set b=ℏ/mω0 (Kreyszig's b2=h/2πmω0, the same quantity), s=x/b and Λ=2E/ℏω0 (his Λ=4πE/ω0h). Multiplying (8) by 2/ℏω0 turns it into
ds2d2ψ~+(Λ−s2)ψ~=0.(9)
Strip the Gaussian. Substituting ψ~=v(s)e−s2/2 into (9): ψ~′′=(v′′−2sv′+(s2−1)v)e−s2/2, so (9) becomes
v′′−2sv′+(Λ−1)v=0,
which is Hermite's equation. It has a polynomial solution — and therefore ψ~ decays at infinity — exactly when
No other Λ is admissible. The displayed condition is what produces a polynomial; what makes these the only energies is that every other solution is thrown out by the normalization (2), and that needs one further fact. Substituting v=∑kaksk into Hermite's equation gives the two-term recursion
(k+2)(k+1)ak+2=(2k+1−Λ)ak,
which couples ak to ak+2 and so splits the solution space into an even series E(s) and an odd series O(s), with v=αE+βO. A series terminates only when Λ=2k+1 for some k of its own parity. Two things then have to be shown, and Kreyszig leaves the first of them as his Problem 5, p. 588.
A non-terminating series grows too fast. Suppose the even series, say, never terminates, so 2k+1−Λ=0 for every even k and no ak of that parity vanishes beyond the first. Fix c with 21<c<1. Comparing the recursion's ratio with the ratio of successive non-zero coefficients of ecs2=∑jcjs2j/j!, which is 2c/(k+2) at k=2j,
and since c<1 the right-hand inequality holds for all k beyond some k0. From k0 on the ratios are positive, so the coefficients keep one sign, and term by term they dominate a positive multiple of those of ecs2. Hence for large s>0 the series exceeds Cecs2 minus the finitely many terms below k0, and
∣ψ~(s)∣=∣v(s)∣e−s2/2≳e(c−1/2)s2⟶∞,
because c>21. Such a ψ~ is not in L2(R) and fails Kreyszig's criterion of staying finite and dying away at infinity.
The two series cannot rescue each other. One might hope that a growing even series and a growing odd series could cancel. They cannot, and parity is why: if ψ~=ve−s2/2 were in L2(R), then so would be its even and odd parts 21[ψ~(s)±ψ~(−s)], which are exactly αE(s)e−s2/2 and βO(s)e−s2/2. So each series would have to be admissible on its own, and by the previous paragraph a non-terminating one is not.
Hence an admissible Λ must terminate one of the two series, which is Λ=2n+1; and at such a Λ the series of the other parity still does not terminate, since 2k+1=Λ only at k=n, so its coefficient must be zero and the admissible solution is the polynomial one alone, up to a constant multiple.
Translating back, E=21ℏω0Λ gives
En=(n+21)ℏω0=(n+21)hν,ν=2πω0,
which is Kreyszig's formula (12). The levels are equally spaced, and the lowest is not zero: the zero-point energy is 21hν. Kreyszig remarks that this half-integral offset is characteristic of the oscillator and was not what Planck assumed in 1900, and he calls n the principal quantum number of the oscillator.
Check the ground state exactly. For n=0, H0=1, so ψ0=Ce−x2/2b2. Then
so E0=21ℏω0, confirming the formula at n=0 with no appeal to Hermite's equation. Normalizing, ∫e−x2/b2dx=bπ, so C=(πb2)−1/4 and
ψ0(x)=(πb2)1/41e−x2/2b2.
A number. For a diatomic vibration with ν=1014 Hz, the spacing is ℏω0=hν=6.62607×10−20 J =0.413567 eV, and the zero-point energy is half of that, 0.206783 eV — comparable with a chemical bond energy, which is why zero-point energy matters in reaction rates and in isotope effects.
Three dimensions. Kreyszig's Problems 9 and 10 take V=21m(ω12x12+ω22x22+ω32x32) and separate (3) in Cartesian coordinates, giving E=∑j(nj+21)ℏωj with a product of Hermite functions. For the isotropic case ω1=ω2=ω3=ω0 he notes that the ground level E0=23hν is non-degenerate while every higher level is degenerate: E=(N+23)ℏω0 is attained by every triple with n1+n2+n3=N, of which there are (2N+2). So the multiplicities of N=0,1,2,3,4 are 1,3,6,10,15.
ApplicationOriginal
Spectroscopy, which is where the energies are read off. Problem 2 of the previous remark is what makes (3) an experimental question: every line in an atomic or molecular spectrum is a difference En−En′, so the eigenvalues are read straight off a photographic plate or a detector. The equally spaced ladder of the oscillator gives the vibrational bands of a diatomic molecule and, through ℏω0=ℏk/μ, measures the bond stiffness k; the n2 ladder of the well gives the absorption edge of a quantum dot, which is why the colour of a cadmium-selenide nanocrystal is set by its diameter.
Hydrogen and the periodic table. Riley, Hobson and Bence's hydrogen equation is (3) with V=−e2/4πε0r. Because the potential is central, the general point of step 5 in separation of variables applies: the angular factors are the spherical harmonics Yℓm, so the quantum numbers ℓ and m and the 2ℓ+1-fold degeneracy are fixed by the Laplacian alone and are the same for every central potential. Only the radial equation knows about the Coulomb tail, and solving it gives the Laguerre functions and the energies En=−13.6eV/n2. Everything in the shell structure of the elements that does not depend on the particular potential — the s,p,d,f labels, the counts 2,6,10,14 — comes from the separation, not from the physics of the Coulomb force.
Tunnelling, and the sign of E−V. Read Kreyszig's derivation backwards: the local wavenumber is k(x)2=2m(E−V(x))/ℏ2, so (3) is a Helmholtz equation where E>V and a modified Helmholtz equation where E<V. In the first region the solution oscillates; in the second it grows or decays exponentially. A particle therefore has a nonzero amplitude inside a barrier it could not classically enter, and the transmission through a barrier of width d goes like e−2κd with ℏκ=2m(V−E). That exponential is the operating principle of the scanning tunnelling microscope — with a work function of about 4 eV, κ=1.02×1010m−1 and a change of one angstrom in the tip height changes the current by about a factor of ten, since ln10/2κ=1.12 Å — and of the flash memory cell.
The classical limit, and semiclassical computation. Hassani's Proposition 14.8.2 is the precise statement of the correspondence principle in this setting: substituting Ψ=AeiS/ℏ turns (1) into a pair of real equations, one of which is the continuity equation (7) for ρ=A2 and J=ρ∇S/m, and the other of which becomes the Hamilton–Jacobi equation for S when the term 2mℏ2∇2A/A is dropped. So as ℏ→0 a wave function behaves like a fluid of non-interacting classical particles following the characteristics of the Hamilton–Jacobi equation, with ∇S/m as the velocity field. Keeping one order further in ℏ gives the WKB approximation, and matching WKB solutions across a turning point produces the Bohr–Sommerfeld condition ∫ab2m(E−V)dx=(n+21)πℏ, which Hassani derives from the connection formulas and which reproduces the oscillator's (n+21) exactly.
Numerical quantum mechanics. Formula (5) is the design of a spectral time-stepper: expand the initial state in eigenfunctions, multiply each coefficient by the unit-modulus factor e−iEnt/ℏ, and sum. Because the factors have modulus one, the scheme is unconditionally stable and exactly norm-preserving, which an explicit finite-difference scheme for (1) is not — the standard alternative, the Crank–Nicolson or Cayley scheme, is chosen precisely because its propagator (1−iH^Δt/2ℏ)/(1+iH^Δt/2ℏ) is unitary and so inherits the conservation law (7). The split-operator method does the same with the kinetic part diagonalized by a fast Fourier transform (the discrete Fourier transform) and the potential part diagonal in position, alternating between the two.
Introduces the Klein-Gordon equation as the wave equation with a mass term, derives its dispersion relation and reads the relativistic energy-momentum relation off it, verifies its one-dimensional source function across the light cone and records the three-dimensional one, shows that the massless limit gives back the wave equation and the non-relativistic limit the Schrodinger equation, and computes the conserved density that fails to be positive.
Definition(The Klein-Gordon equation)AdaptedAdapted from Strauss (2008), §13.5, Eq. (1), pp. 373–374; Hassani (2013), §13.1, Eq. (13.6), p. 396; Gelfand & Fomin (1963), §38.1, Example 3, Eqs. (119)–(120), p. 181
The Klein-Gordon equation is the wave equation with a zeroth-order term,
utt−c2∇2u+μ2u=0,(1)
in which c is the propagation speed and μ is a positive constant. It is Strauss's equation (1) of §13.5, which he introduces as one of the simplest of the hyperbolic models of elementary particles and whose only difference from the wave equation, as he puts it, is the last term. Gelfand and Fomin identify the field it describes: the scalar field of uncharged spin-zero particles of mass M, the π0 meson being their example.
The constant in the mass term, and the three conventions in our sources. Strauss writes the last term as m2u with m called the mass of the particle. For the three terms of (1) to have the same dimensions, the coefficient must have the dimensions of an inverse time; in physical units it is
μ=ℏMc2,
M being the mass, and c/μ=ℏ/Mc is the reduced Compton wavelength of the particle (the Compton wavelength itself being 2π times that). Hassani avoids the issue by working in natural units ℏ=c=1, in which his equation (13.6) reads ∇2ϕ−m2ϕ=∂2ϕ/∂t2 — that is (1) with c=1 and μ=m. We write μ for the coefficient, state the relation μ=Mc2/ℏ once, and keep c explicit.
A sign to be careful with. Gelfand and Fomin display the equation as (□−M2)u=0 and give the Lagrangian density
L=21(ux02−ux12−ux22−ux32−M2u2),(2)
of which they say (1) is the Euler equation for a field (their Eq. (118), the multiple-integral Euler equation of several independent variables). Working out that Euler equation from (2) gives ux0x0−∇2u+M2u=0, which is (1) with c=1; so the sign of the mass term implied by their Lagrangian is the one Strauss and Hassani use, and their displayed (□−M2)u=0 matches it only for the metric signature in which □=∇2−∂x02. In our scanned copy the signs inside their d'Alembertian are not legible. We follow Strauss and Hassani, equation (1), whose mass term has the sign that makes the solutions oscillate rather than grow.
Application(Plane waves and the dispersion relation)AdaptedAdapted from Hassani (2013), §27.1.1, pp. 832–833; Strauss (2008), §13.5, Eq. (13), p. 376
The function u=Aei(k⋅x−ωt) with amplitude A=0 solves (1) if and only if
ω2=c2∣k∣2+μ2.(3)
Equation (3) is the dispersion relation of the Klein-Gordon equation. Writing k=∣k∣, it has the following consequences.
Phase and group velocity. For k=∣k∣>0,
vph=kω=c1+c2k2μ2>c,vgr=dkdω=ωc2k<c,vphvgr=c2.
The phase velocity exceeds c and the group velocity falls short of it; their product is exactly c2. At k=0 the phase velocity is undefined, since ω/k diverges, and the group velocity is 0, so the product identity needs k>0. For μ=0 both equal c and the medium is dispersionless.
The relativistic energy-momentum relation. Substituting E=ℏω, p=ℏk and μ=Mc2/ℏ turns (3) into
E2=p2c2+M2c4,
and then vgr=c2k/ω=pc2/E, which is the velocity of a relativistic particle of momentum p and energy E.
Both signs of ω occur. For each k there are two roots ω=±c2k2+μ2, and both give solutions of (1). There is no k for which ω lies between −μ and μ: the spectrum has a gap of width 2μ about zero.
Hassani reads the correspondence in the other direction, and this is how the equation was found: replacing E by iℏ∂/∂t and p by −iℏ∇ turns the non-relativistic relation E=p2/2M into the Schrodinger equation and the relativistic relation E2−p2c2=M2c4 into (1).
Idea of the proofOriginal
Substituting an exponential into a constant-coefficient equation turns each derivative into a multiplication, so the whole equation collapses to one algebraic relation among ω and k — the dispersion relation. Everything qualitative about the equation can be read off it. Here the mass term contributes a constant μ2 that is added to c2k2, which has three visible effects: long waves (ck≪μ) are strongly dispersive and hardly propagate, short waves (ck≫μ) behave as if the mass were absent, and there is a lowest frequency μ below which no travelling wave exists at all. The last of these is why a massive field cannot carry a static long-range force, and the first is why a Klein-Gordon pulse spreads where a wave-equation pulse does not.
ProofAdaptedAdapted from Hassani (2013), §27.1.1, p. 832
Step 1 (substitute). For u=Aei(k⋅x−ωt) with A=0,
utt=−ω2u,∇2u=−∣k∣2u,
so (1) becomes (−ω2+c2∣k∣2+μ2)u=0. Since u vanishes nowhere, this holds identically if and only if (3) does.
Step 2 (the velocities). From (3), ω=c2k2+μ2 for the positive root, so
kω=kc2k2+μ2=c1+c2k2μ2,
which exceeds c for every k>0 when μ>0. Differentiating (3) in k gives 2ωdω/dk=2c2k, hence vgr=c2k/ω, and multiplying the two,
vphvgr=kω⋅ωc2k=c2.
Since ω>ck by (3), vgr=c2k/ω<c.
Step 3 (the gap). By (3), ω2≥μ2 for every real k, with equality only at k=0. So ∣ω∣≥μ, and the two branches ω=±c2k2+μ2 are separated by the interval (−μ,μ), which no plane wave occupies. □
Application(The one-dimensional source function, and the light cone)AdaptedAdapted from Strauss (2008), §13.5, Eqs. (12)–(14) and (20), pp. 376–378
where J0 is the first-kind Bessel function of order zero. Strauss obtains it by Fourier transform: the transform satisfies S^tt=−(c2k2+μ2)S^ with S^(k,0)=0 and S^t(k,0)=1, hence
S^(k,t)=c2k2+μ2sin(tc2k2+μ2),(5)
and inverting (5) gives (4). Two features are worth naming.
The support is the same light cone as for the wave equation.S vanishes outside ∣x∣≤ct, and Strauss notes that it has the same jump discontinuity on the cone as the wave equation's source function: the mass term changes the solution inside the cone but moves no signal faster, and the domain of dependence is unchanged from d'Alembert's formula.
The massless limit is the wave equation. As μ→0, J0→J0(0)=1 and (4) becomes S=2c1 on ∣x∣<ct and 0 outside, which is the source function of the one-dimensional wave equation.
Carrying the same Fourier calculation out in three dimensions, Strauss arrives at
where r=∣x∣. The first term is an impulse carried on the cone itself; the second is spread through the cone's interior and is built from J1. At μ=0 the second term disappears and (6) is the wave equation's source function, concentrated on the cone alone. So Huygens' principle — a sharp pulse staying sharp — holds for the wave equation in three dimensions and fails for Klein-Gordon: a mass term makes the cone leak into its interior, and an impulse is followed by a decaying tail.
Idea of the proofOriginal
The Fourier transform diagonalizes (1), as it does any constant-coefficient equation: in k the problem is the oscillator S^tt+ω(k)2S^=0 with ω(k)2=c2k2+μ2, whose solution with the stated data is sin(ωt)/ω. All the information about the mass is in the single function ω(k), and the only question left is what the inverse transform of sin(tω(k))/ω(k) is.
That question has a pretty answer because of a coincidence of integral representations: J0 has one whose integrand is exactly sin(tc2k2+μ2)/c2k2+μ2. Rather than reproduce Strauss's table lookup, the proof below verifies (4) directly — it is a short computation, and it shows where the Bessel equation comes from. The key is that for a function of the invariant combination s=t2−x2/c2, the wave operator reduces to an ordinary differential operator in s which is precisely Bessel's operator of order zero.
Two things then remain that a pointwise substitution cannot give. The first is the cone itself: (4) is discontinuous there, and a jump across a surface normally leaves δ and δ′ terms in the wave operator, so solving on each side separately proves nothing about the union. They cancel here for two reasons worth separating — the cone is characteristic, which kills the δ′ term, and the Bessel profile is a function of t2−x2/c2, which vanishes on the cone to exactly the order that kills the δ term. The second is the initial data, and those follow from the shrinking support alone: S lives on [−ct,ct] and is bounded by 1/2c, so it tends to 0, while its integral comes out as t+O(t3), which is what makes its t-derivative at 0 a unit point mass.
ProofOriginal
We verify (4) by substitution. Set c=1 throughout; the general case follows by rescaling x↦x/c, which leaves (1) in the form utt−uxx+μ2u=0.
Step 1 (the invariant variable). Inside the forward cone ∣x∣<t put
s=t2−x2>0,
so that st=t/s and sx=−x/s. Differentiating once more,
stt=s1−s3t2,sxx=−s1−s3x2.
Step 2 (the wave operator on a function of s). For u=f(s), the chain rule gives utt=f′′st2+f′stt and uxx=f′′sx2+f′sxx, so
Step 3 (Bessel's equation). Take f(s)=J0(μs). Then f′=μJ0′(μs) and f′′=μ2J0′′(μs), so by (7) and Bessel's equation of order zero, J0′′(z)+J0′(z)/z+J0(z)=0 with z=μs,
Therefore utt−uxx+μ2u=0: the function 21J0(μt2−x2) solves the Klein-Gordon equation at every point strictly inside the cone. (With c restored, s=t2−x2/c2 and the constant is 1/2c, which is (4).)
Step 4 (across the cone). Outside the cone S≡0, which solves (1) trivially. The two pieces meet on ∣x∣=t, where the interior limit is 21J0(0)=21, so S jumps by 21 across the cone (by 1/2c once c is restored). Solving on each of the two open sets separately is not enough, because a jump across a surface in general leaves δ and δ′ terms in the wave operator, and what the statement asserts is that S solves the Cauchy problem on the whole line. We show those terms cancel.
Two preliminaries. First, J0 is even and entire, so
J0(μw)=j≥0∑(j!)2(−1)j(4μ2w)j
is an entire function of w with no branch point at w=0. Hence, putting G(w)=21J0(μw) and
g(x,t)=G(t2−x2),ϕ(x,t)=t−∣x∣,
the profile g is smooth on the whole plane — not merely inside the cone — and S=H(ϕ)g with H the Heaviside step. Second, Step 3 says exactly that gtt−gxx=−μ2g wherever ϕ>0.
Work on the open half-plane t>0. Near the axis x=0, the only place where ϕ fails to be differentiable, ϕ>0, so H(ϕ)≡1 and S=g on a neighbourhood; there the equation holds by Step 3. It therefore suffices to work where x=0, and there ϕ is smooth with
ϕt=1,ϕx=−sgnx,ϕtt=ϕxx=0.
Differentiating the product S=H(ϕ)g twice in the sense of distributions, ∂tS=δ(ϕ)ϕtg+H(ϕ)gt and
The δ′ term vanishes because ϕt2−ϕx2=1−1=0: the level sets of ϕ are characteristic for the wave operator. This is where the geometry of the light cone enters, and it is the precise content of the statement that the mass term moves no signal faster.
The second term vanishes because ϕ is affine on each of the two half-planes x>0 and x<0.
The third vanishes because g depends on x and t only through t2−x2. Indeed gt=2tG′(t2−x2) and gx=−2xG′(t2−x2), so
ϕtgt−ϕxgx=gt+(sgnx)gx=2(t−∣x∣)G′(t2−x2)=2ϕG′(t2−x2),
and δ(ϕ)ϕ=0. So the Bessel profile vanishes on the cone to exactly the order needed.
The last term is −μ2H(ϕ)g=−μ2S, by Step 3.
Hence Stt−Sxx+μ2S=0 in the sense of distributions on {t>0}, across the cone and not only on either side of it. The jump 1/2c is the same one the wave equation's source function has, which is Strauss's observation; the two computations above are why the mass term leaves it alone.
Step 5 (the initial data). It remains to show that S(⋅,t)→0 and St(⋅,0)=δ as distributions in x when t↓0. Everything needed is the shrinking support suppS(⋅,t)⊆[−t,t] and the bound ∣S∣≤21, which holds because ∣J0∣≤1.
First integrate (4) over its support. The series above converges uniformly on compact sets, so it may be integrated term by term:
Substituting x=tu and using the beta integral ∫−11(1−u2)jdu=2⋅4j(j!)2/(2j+1)!, the inner integral is t2j+1⋅2⋅4j(j!)2/(2j+1)!, and every factorial cancels:
So ⟨S(⋅,t),φ⟩=tφ(0)+O(t2) for every φ. Letting t↓0 gives S(⋅,0)=0; dividing by t first and using that value,
⟨St(⋅,0),φ⟩=t↓0limt1⟨S(⋅,t)−S(⋅,0),φ⟩=φ(0),
which says St(⋅,0)=δ. Neither step uses (5), which is the point: the Fourier route is not needed anywhere in this proof. (The two do agree: the integral just computed is S^(0,t), which is the value of (5) at k=0.)
Step 6 (the massless limit).J0 is continuous with J0(0)=1, and μt2−x2/c2→0 uniformly on compact subsets of the cone as μ→0. So S→2c1 pointwise inside the cone and 0 outside, which is the d'Alembert source function. □
Application(The conserved density, and the non-relativistic limit)Original
Let u be a complex-valued solution of (1). Define
ρ=2Mc2iℏ(uˉut−uuˉt),J=MℏIm(uˉ∇u).(8)
Then ∂tρ+∇⋅J=0. Moreover:
ρ is not positive. For the plane wave u=Aei(k⋅x−ωt),
ρ=Mc2ℏω∣A∣2,J=Mℏk∣A∣2,ρJ=ωc2k=vgrk^,
so ρ has the sign of ω and is negative on the branch ω<0 of the dispersion relation.
Non-relativistic limit. Write u=e−iω0tψ with ω0=Mc2/ℏ=μ. Then ψ satisfies
iℏ∂t∂ψ=−2Mℏ2∇2ψ+2Mc2ℏ2∂t2∂2ψ,(9)
exactly. Dropping the last term — legitimate when ∣ψtt∣≪ω0∣ψt∣, that is when the residual frequencies are small beside Mc2/ℏ — gives the free Schrodinger equation, and in the same limit ℏω→Mc2 makes ρ→∣ψ∣2, the Schrodinger probability density.
Hassani records the historical consequence of part 1: that the equation's second time derivative made the resulting probabilities capable of coming out negative, that it was set aside at first for that reason, and that it was nevertheless readmitted in the end as a proper equation for relativistic particles; and the wish for a first-order equation is what led Dirac to factor E2−p2−M2 and so to the Dirac equation. In Reed and Simon's framework the point is structural: their model demands a one-parameter unitary group generated by a self-adjoint H, with ∥Ψ∥=1 fixed; (1) is second order in t, so its natural phase space is the pair (u,ut), and the conserved quantity that does come out of it — the energy — is quadratic in the derivatives rather than in u alone.
Idea of the proofOriginal
The computation is the same two lines as for the Schrodinger current, run on a second-order equation. Multiply (1) by uˉ, multiply its conjugate by u, subtract: the terms in μ2∣u∣2 cancel because μ2 is real, and the second time derivatives assemble into ∂t(uˉut−uuˉt) while the Laplacians assemble into a divergence. So one gets a conservation law for free. The difference from the Schrodinger case is where the time derivatives sit: there the equation was first order, so the conserved density came out as ∣Ψ∣2, manifestly non-negative; here it involves ut, and nothing forces its sign.
For the non-relativistic limit, the trick is to peel off the rest-energy oscillation e−iMc2t/ℏ, which is the fastest thing in the solution. Doing so cancels the mass term against the −ω02 coming from differentiating the exponential twice — that cancellation is exactly the statement μ=ω0 — and leaves a first-order equation in t plus one leftover second derivative, suppressed by 1/c2.
ProofOriginal
Step 1 (the conservation law). Write (1) as utt=c2∇2u−μ2u and its conjugate as uˉtt=c2∇2uˉ−μ2uˉ, using that c and μ are real. Put
and J=MℏIm(ik∣A∣2)=Mℏk∣A∣2. Dividing, J/ρ=c2k/ω, whose magnitude is vgr by part 1 of the dispersion statement. Since both roots ω=±c2k2+μ2 are admissible, ρ takes both signs.
Step 3 (peel off the rest energy). Put u=e−iω0tψ with ω0=μ=Mc2/ℏ. Then
utt=e−iω0t(ψtt−2iω0ψt−ω02ψ),
and substituting into (1), after dividing by e−iω0t,
ψtt−2iω0ψt−ω02ψ−c2∇2ψ+μ2ψ=0.
Because μ=ω0, the two terms −ω02ψ and +μ2ψ cancel — this is the whole point of the substitution — leaving
ψtt−2iω0ψt−c2∇2ψ=0.
Step 4 (rearrange). Multiply by 2ω0ℏ=2Mc2ℏ2 and move the ψt term to the left:
iℏψt=2Mc2ℏ2ψtt−2Mℏ2∇2ψ,
which is (9). If ψ varies on a time scale T with ω0T≫1, the retained term iℏψt is of size ℏ∣ψ∣/T while the dropped term is of size ℏ2∣ψ∣/(2Mc2T2)=ℏ∣ψ∣/(2ω0T2), smaller by the factor 1/(2ω0T). Dropping it gives the free Schrodinger equation.
Step 5 (the density in the limit). In the same regime ut=e−iω0t(ψt−iω0ψ)≈−iω0e−iω0tψ, so
ρ≈2Mc2iℏ(−iω0−(−iω0))∣ψ∣2=Mc2ℏω0∣ψ∣2=∣ψ∣2,
using ℏω0=Mc2. So the Klein-Gordon density reduces to the Schrodinger probability density, and in particular becomes positive. □
RemarkOriginal
Where the hypotheses are used, and what changes without them.
μ real, in Step 1 of the last proof and in Step 3 of the dispersion proof. Both the cancellation of the mass term in the conservation law and the reality of ω need it. Strauss's Exercise 9 shows what an imaginary mass does: substituting u=e−t/2v into the telegraph equationutt−c2Δu+νut=0 with ν=1 gives vtt−c2Δv−41v=0, which is (1) with μ=i/2. (His exercise prints the conclusion in u where v is meant; it is the equation satisfied by v.) So a dissipative wave equation is a Klein-Gordon equation with imaginary mass — but the imaginary mass is not itself the damping. With μ=i/2 the dispersion relation (3) reads ω2=c2k2−41, whose roots are real for ck>21 and purely imaginary for ck<21, never genuinely complex: an imaginary mass closes the frequency gap of part 3 of the dispersion statement and replaces oscillation by the pair e±λt, λ=41−c2k2, in the long-wave band. One of those branches grows — the tachyonic instability — so (1) with imaginary μ is unstable, not damped. The telegraph equation's damping is the prefactor e−t/2 that the substitution strips off, and it is uniform in k; the two fit together because λ≤21 for every k, with equality only at k=0, so u=e−t/2v grows at no wavenumber, and decays strictly at every k=0. At k=0 exactly the two factors cancel and u is constant, as it must be: the k=0 telegraph equation is utt+ut=0, whose solutions are a+be−t, and the constant a never dies. Reading the correspondence the other way: a real mass term makes waves disperse but not decay, while an imaginary one destroys the oscillation of the long waves altogether.
The classification is unchanged by the mass term. Strauss's index lists Klein-Gordon under hyperbolic equations, and his point in §13.5 is that the term μ2u is of order zero, so the principal part of the operator — the part that decides the type, the characteristics and the domain of dependence (classification of second-order equations) — is that of the wave equation. That is why (4) and (6) are supported in the light cone and why Strauss's remark about the jump on the cone holds: the mass changes the solution inside the cone and nothing about the cone itself.
Huygens' principle fails. The second term of (6) is the price of the mass. In three dimensions the wave equation's source function is a pure delta on the cone, so a flash is seen as a flash; with μ>0 a flash is followed by a reverberating tail ∝J1, decaying but never ending. The same contrast appears in the wave equation in three dimensions, where the two-dimensional wave equation also has a tail.
What the "negative probability" means and what it does not. Part 1 of the last statement is a theorem: the conserved density (8) changes sign. Hassani reports the historical reading — the equation was initially abandoned because of it — together with the resolution, that the equation was in the end admitted as a proper description of relativistic particles after all. The modern reading is that ρ is a charge density and not a probability density, so its sign is a sign of charge; the equation is then a field equation rather than a one-particle wave equation, and its quantization describes particles and antiparticles. Gelfand and Fomin take the field view from the start: (1) is the Euler equation of the Lagrangian density (2), and §38.2 obtains its conservation laws from Noether's theorem, integrating the divergence identities over the boundary of a space-time cylinder.
The square-root equation. Strauss records the precise sense in which Dirac's equation is a square root of (1). Write D for the first-order operator built from his gamma matrices,
D=c1γ0∂t∂+γ1∂x∂+γ2∂y∂+γ3∂z∂,
whose square is D2=c−2∂t2−Δ, the gamma-matrix relations killing every cross term. Strauss's mass symbol here is a lowercase m of inverse-length dimension — it is μ/c in our notation, which is what lets m2 be added to −Δ, and it is not the M of μ=Mc2/ℏ — and with it
(D+im)(D−im)=D2+m2=c21∂t2∂2−Δ+m2,
which is c−2 times the operator in (1). Strauss displays this as the square(D+im)2; read as a square it would leave the cross term 2imD and reverse the sign of the mass, so it is the product of the two first-order factors, as displayed above, that is the Klein-Gordon operator. Hassani gives the algebraic version in exactly that product form: Dirac looked for β,αi making (βE+∑αipi−m)(βE+∑αipi+m) reduce to E2−∑pi2−m2, which is what forces the coefficients to be matrices satisfying a Clifford algebra.
The massless limit is not merely a limit of solutions. At μ=0 the dispersion relation (3) loses its gap, the phase and group velocities both become c, and the source functions (4) and (6) become the wave equation's. But the character of the field also changes qualitatively: only a massless field has a static solution falling off as 1/r and so a long-range force, which is why the photon is massless and the force carried by a massive field is screened at the length c/μ=ℏ/Mc — the same screening as in the modified Helmholtz equation of the Helmholtz equation, with e−μr/c/r in place of 1/r.
Example(Dispersion for a pion, and the exact relation vphvgr=c2)Original
Take the charged pion π±, with rest energy Mc2=139.57039 MeV. (A charged spin-zero particle is described by a complex-valued solution of (1), which is the case the conserved density (8) applies to, with ρ read as a charge density. Gelfand and Fomin's uncharged π0 is a real field and is lighter, Mc2=134.9768 MeV, giving c/μ=1.461933 fm in place of the 1.413817 fm below, and an ℏω column smaller by the factor 134.9768/139.57039; the velocity and product columns are the same for either particle, since they depend only on the ratio ck/μ.) Using ℏc=197.3269804 MeV fm, the mass constant and length scale of (1) are
the reduced Compton wavelength. Evaluate (3) at the three wavenumbers ck=21μ, ck=μ and ck=2μ, writing ℏω in MeV:
ck/μ
ω/μ=1+(ck/μ)2
ℏω (MeV)
vgr/c=ck/ω
vph/c=ω/(ck)
product
0.5
1.1180340
156.0444
0.4472136
2.2360680
1
1
1.4142136
197.3823
0.7071068
1.4142136
1
2
2.2360680
312.0889
0.8944272
1.1180340
1
The last column is exactly 1 in every row, which is the identity vphvgr=c2 of Step 2 — an algebraic identity, not a numerical coincidence. At ck=μ, that is at k=Mc/ℏ and so at the Compton wavelength 2πℏ/Mc, the group velocity is c/2=0.7071068c and the phase velocity 2c.
Check against relativistic kinematics. The row ck=μ corresponds to pc=Mc2, hence E=(pc)2+(Mc2)2=2Mc2=197.3823 MeV, which is the third column; the Lorentz factor is γ=E/Mc2=2 and the particle speed v=c1−γ−2=c1−21=c/2, which is the fourth column. So the group velocity of the Klein-Gordon wave packet is the velocity of the relativistic particle, as part 2 of the statement says, and the superluminal phase velocity carries no signal.
The gap. No solution of (1) oscillates more slowly than μ: the lowest frequency, at k=0, corresponds to ℏω=Mc2=139.57039 MeV. For the wave equation (μ=0) there is no such floor, and arbitrarily slow oscillations exist.
Example(The source function solves the equation, and its tail)Original
Verify (4) and watch the massless limit numerically. Put c=1 and μ=3, and take the candidate
S(x,t)=21J0(3t2−x2)(∣x∣<t).
Exact verification. By Step 3 of the proof, Stt−Sxx+9S=0 identically inside the cone; symbolic differentiation of J0(3t2−x2) and substitution confirms that the residual is exactly zero. Concretely at t=2, x=1, so s=3=1.7320508 and μs=5.1961524:
J0(5.1961524)=−0.1116107,S=−0.0558053.
The sign is worth noting: the source function is negative there, which cannot happen for the wave equation, where S≡1/2c inside the cone. The Bessel oscillation is the mass term at work.
The tail. Along the ray x=0 the source function is S(0,t)=21J0(3t), which oscillates and decays like t−1/2:
t
μt
S(0,t)=21J0(3t)
0.1
0.3
0.4888131
0.5
1.5
0.2559138
1
3
−0.1300260
2
6
0.0753226
5
15
−0.0071122
At small t the value is near 1/2c=0.5, the wave-equation answer, because J0(z)≈1−z2/4 for small z; at large t the solution rings. This is the one-dimensional face of the tail inside the light cone that formula (6) exhibits in three dimensions, and it is why Huygens' principle fails for a massive field.
The massless limit, quantitatively. At the fixed point t=1, x=0, so s=1:
μ
S=21J0(μ)
relative departure from 21
1
0.3825988
23.480%
0.3
0.4888131
2.2374%
0.1
0.4987508
0.24984%
0.01
0.4999875
0.0024999%
The departure is 41μ2 to leading order — 0.25, 0.0225, 0.0025, 0.000025 against the table's 0.23480, 0.022374, 0.0024984, 0.000024999 — so the solution approaches the d'Alembert source function quadratically in the mass, as Strauss's remark that (4) converges to the wave equation's source function at m→0 requires.
ApplicationOriginal
Why the range of a force is set by the mass of its carrier. Strip the time dependence from (1) with μ>0 and one gets the modified Helmholtz equation ∇2u−(μ/c)2u=0, whose decaying radial solution is e−μr/c/r rather than 1/r. So a static source of a massive field produces a potential screened at the length c/μ=ℏ/Mc. This is Yukawa's argument: from the observed range of the nuclear force, about 1.4 fm, one predicts a carrier of rest energy about 140 MeV, which is the pion — the number computed in the first example above. The massless limit removes the screening and restores the 1/r Coulomb and Newtonian potentials of the fundamental solution of Laplace's equation, which is why electromagnetism and gravity are long-ranged.
Waves in a plasma, and in a waveguide. The relation (3) is not special to particle physics. Transverse electromagnetic waves in a cold plasma obey exactly it, with μ the plasma frequency ωp: below ωp no wave propagates and the plasma reflects, which is how shortwave radio reaches over the horizon and why radio below a few megahertz cannot escape the ionosphere. A mode in a hollow waveguide obeys it too, with μ=cκ the cutoff set by the transverse eigenvalue of the Helmholtz equation; the superluminal phase velocity and subluminal group velocity of the table above are measured facts about waveguides. In each case the Klein-Gordon equation is the equation of the field in a medium with a gap, and the "mass" is whatever sets the gap.
Dispersion, and the spreading of pulses. Because vgr depends on k, a Klein-Gordon pulse spreads, unlike a wave-equation pulse in one dimension. Expanding (3) about a carrier wavenumber k0 gives ω≈ω0+vgr(k−k0)+21ω′′(k0)(k−k0)2, and the quadratic term ω′′=c2μ2/ω3 is the group-velocity dispersion that broadens the envelope — the same quantity that limits the bit rate of an optical fibre. The exact source function (4) exhibits the extreme case: an impulse does not stay an impulse but develops the ringing tail tabulated above.
The non-relativistic limit in practice. Part 2 of the conserved-density statement is the formal content of the statement that quantum mechanics is the low-energy limit of relativistic field theory, and equation (9) is used as written. Keeping the term 2Mc2ℏ2ψtt to first order and eliminating ψtt with the leading-order equation produces the relativistic kinetic-energy correction −p4/8M3c2, the first of the fine-structure corrections to atomic energy levels and one of the standard tests of the Dirac theory against the Schrodinger theory of the Schrodinger equation.
Numerics: why the two-sign spectrum matters. A time-stepping scheme for (1) must handle both branches of the dispersion relation, and because ω grows only linearly in k the stability constraint is the usual Courant condition cΔt≤Δx of a hyperbolic solver, with the mass term contributing a harmless O(μ2Δt2). The genuinely awkward regime is μΔt≫1, where the carrier oscillation e−iμt must be resolved even though nothing physical happens on that timescale; the standard remedy is the substitution of Step 3 — factor out e−iμt and integrate (9), whose time scale is the slow one. This is the same device that underlies the derivation of envelope equations such as the nonlinear Schrodinger equation from the nonlinear Klein-Gordon equation that Strauss writes down as his meson model.