pdz

Chapter 18 PartialDifferentialEquationsin Three–DimensionalSpace At last we have ascended the dimensional ladder to its ultima te rung (at least for those of us living in a three-dimensional universe): partia l differential equations in physical space. As in the one and two-dimensional settings developed in the preceding chapters, the three key examples are the three-dimensional Laplace eq uation, modeling equilibrium configurationsof solidbodies, thethree-dimensional wave equation, governing vibrationsof solids, liquids, gasses, andelectromagneticwaves, andth ethree-dimensional heat equation, modeling basic spatial diffusion processes. Fortunately, almost everything of importance has already a ppeared in the one- and two-dimensional situations, and appending a third dimensi on is, for the most part, sim- ply a matter of appropriately adapting the constructions. W e have already seen the basic underlying solutiontechniques: separationofvariables a ndGreen’s functions or fundamen- tal solutions. (Unfortunately, the most powerful of our pla nar tools, conformal mapping, doesnotcarry over to higher dimensions.) In three-dimensional pro blems, separation of variables is applicable in rectangular, cylindrical and sp herical coordinates. The first two do not produce anything fundamentally new, and are therefor e relegated to the exercises. Separation in spherical coordinates leads to spherical har monics and spherical Bessel func- tions, whose properties are investigated in some detail. Th ese new special functions play important roles in a number of physical systems, including t he quantum theory of atomic structure that underlies the spectral and chemical propert ies of atoms. The Green’s function for the three-dimensional Poisson equ ation in space can be identified as the classic Newtonian (and Coulomb) 1 /rpotential. The fundamental so- lution for the three-dimensional heat equation can be easil y guessed from its one- and two-dimensional versions. The three-dimensional wave equ ation, surprisingly, has an ex- plicit, although more intricate, solution formula of d’Ale mbert form, due to Poisson. Para- doxically, the best way to treat the two-dimensional versio n is by “descending” from the simpler three-dimensional formula. This result highlight s a remarkable difference between waves in planar and spacial media. In three-dimensions, Huy gens’ principle states that waves emanating from a localized initial disturbance remai n localized as they propagate through space. In contrast, in two dimensions, initially co ncentrated pulses leave a slowly decaying remnant that never entirely disappears. 18.1. The Laplace and Poisson Equations. We begin our investigations, as usual, with systems in equil ibrium, deferring dynamics 12/11/12 972 c/ci∇cleco√y∇t2012 Peter J. Olver until later. The prototypical equilibrium system is the thr ee-dimensional Laplace equation ∆u=∂2u ∂x2+∂2u ∂y2+∂2u ∂z2= 0, (18.1) in which x= (x,y,z)Trepresents rectangular coordinates on R3. The solutions u(x,y,z) continue to be known as harmonic functions . The Laplace equation models unforced equilibria; Poisson’s equation is the inhomogeneous version −∆u=f(x,y,z), (18.2) where the inhomogeneity frepresents some form of external forcing. The basic boundary value problem for the Laplace or the Poiss on equation seeks a solution inside a bounded domain Ω ⊂R3subject to either Dirichlet boundary conditions , prescribing the function values u=hon∂Ω, (18.3) orNeumann boundary conditions prescribing its normal derivative or flux ∂u ∂n=kon∂Ω, (18.4) ormixed boundary conditions in which one imposes Dirichlet conditions on part of the boundary and Neumann conditions on the remainder. Keep in mi nd that the boundary of the solid domain Ω consists of one or more piecewise smooth cl osed surfaces, which will be oriented by the outwards unit normal n; Theboundaryvalueproblemsforthethree-dimensionalLapl aceandPoissonequations govern a wide variety of physical systems, including: (a)Heat conduction : In this application, urepresents the equilibrium temperature in a solid body. Dirichlet conditions correspond to fixing the t emperature on the bounding surface(s), whereas homogeneous Neumann conditi ons correspond to an insulated boundary, i.e., one which does not allow any heat fl ux. The inhomogene- ityfrepresents some form of internal heat source. (b)Ideal fluid flow : Hereurepresents the velocity potential for an incompressible, i rro- tational steady state fluid flow inside a container, Ω, with ve locity vector field v=∇u. Homogeneous Neumann boundary conditions correspond to a s olid boundary which the fluid cannot penetrate. (c)Elasticity : In certain restricted situations, urepresents an equilibrium deformation of a solid body, e.g., the radial deformation of a solid ball. Fu lly three-dimensional elasticityisgovernedbyamorecomplicatedsystemofparti aldifferentialequations, which can be found in Example 21.9. (d)Electrostatics : In applications to electromagnetism, urepresents the electric potential in a conducting medium; its gradient ∇uprescribes the electromotive force on a charged particle. The inhomogeneity represents an electos tatic force field. 12/11/12 973 c/ci∇cleco√y∇t2012 Peter J. Olver (e)Gravitation : The Newtonian gravitational potential in flat empty space i s also pre- scribed by the Laplace equation. (In contrast, general rela tivity requires a vastly more complicated nonlinear system of partial differential e quations, [ 132].) Self–Adjoint Formulation and Minimum Principle The Laplace and Poisson equations naturally fit into the gene ral self-adjoint equi- librium framework summarized in Section 14.7. The construc tion is a straightforward adaptation of the planar version of Section 15.4. We introdu ce the L2inner products /an}b∇acketle{tu;/tildewideu/an}b∇acket∇i}ht=/integraldisplay/integraldisplay/integraldisplay Ωu(x,y,z)/tildewideu(x,y,z)dxdydz, /an}b∇acketle{tv;/tildewidev/an}b∇acket∇i}ht=/integraldisplay/integraldisplay/integraldisplay Ωv(x,y,z)·/tildewidev(x,y,z)dxdydz,(18.5) between scalar fields u,/tildewideu, and between vector fields v,/tildewidevdefined on the domain Ω ⊂R3. We assume that the functions in question are sufficiently nice that these inner products are well-defined; ifΩisunbounded, thisrequiresthat, atlarge distances, theydecay reasonably rapidly to zero. When subject to suitable homogeneous boundary conditions, the three-dimensional Laplace equation can be placed in our standard self-adjoint form −∆u=−∇·∇u=∇∗◦∇u. (18.6) This relies on the fact that the adjoint of the gradient opera tor with respect to the L2 inner products (18.5) is minus the divergence operator: ∇∗v=−∇·v. (18.7) As usual, the determination of the adjoint rests on an integr ation by parts formula, which, in three-dimensional space, follows from the Divergence Th eorem B.35. The first step is to establishthethree-dimensional analogofGreen’s formula (15.88). Weapplythedivergence identity (B.85) to the product uvof a scalar field uand a vector field v, leading to /integraldisplay/integraldisplay/integraldisplay Ω/parenleftbig u∇·v+∇u·v/parenrightbig dxdydz =/integraldisplay/integraldisplay/integraldisplay Ω∇·(uv)dxdydz =/integraldisplay/integraldisplay ∂Ωu(v·n)dS.(18.8) Rearrangingthetermsproducesthedesiredintegrationbyp artsformulafortripleintegrals: /integraldisplay/integraldisplay/integraldisplay Ω(∇u·v)dxdydz =/integraldisplay/integraldisplay ∂Ωu(v·n)dS−/integraldisplay/integraldisplay/integraldisplay Ωu(∇·v)dxdydz. (18.9) The boundary integral will vanish provided either u= 0 orv·n=0at each point on ∂Ω. Whenu= 0 on all of ∂Ω, we have homogeneous Dirichlet conditions. Setting v·n=0 everywhere on ∂Ω results in the homogeneous Neumann boundary value problem ; see Section 15.4 for a detailed explanation. Finally, when u= 0 on part of ∂Ω andv·n=0on the rest leads to the mixed boundary value problem. Thus, sub ject to one of these choices, the integration by parts formula (18.9) reduces to /an}b∇acketle{t∇u;v/an}b∇acket∇i}ht=/an}b∇acketle{tu;−∇·v/an}b∇acket∇i}ht, (18.10) which suffices to prove the adjoint formula (18.7). 12/11/12 974 c/ci∇cleco√y∇t2012 Peter J. Olver Remark: Adopting more general weighted inner products results in a more general elliptic boundary value problem. See Exercise for details. According to the abstract Theorem 7.60, the self-adjoint fo rmulation (18.6) implies positivesemi-definitenessoftheboundaryvalueproblem, a ndpositivedefinitenessprovided ker∇={0}. Since, on a connected domain, only constant functions are a nnihilated by the gradient operator — see Theorem B.28 — both the Dirichlet and mixed boundary value problems are positive definite, while the Neumann boun dary value problem is only semi-definite. Finally, in the positive definite cases, the solution can be c haracterized by the three- dimensional version of the Dirichlet minimization princip le (15.103). Theorem 18.1. The solution u(x,y,z)to the Poisson equation (18.2)subject to homogeneous Dirichlet or mixed boundary conditions (18.3)is characterized as the unique function that minimizes the Dirichlet integral 1 2/ba∇dbl∇u/ba∇dbl2−/an}b∇acketle{tu;f/an}b∇acket∇i}ht=/integraldisplay/integraldisplay/integraldisplay Ω/bracketleftbig1 2(u2 x+u2 y+u2 z)−fu/bracketrightbig dxdydz (18.11) among all C1functions that satisfy the prescribed boundary conditions . As in the two-dimensional version discussed in Section 15.4 , he minimization principle continuestoholdwithoutmodificationinthecaseoftheinho mogeneousDirichletboundary value problem. Modifications for the inhomogeneous mixed bo undary value problem are discussed in Exercise . The three-dimensional finite element method for construct ing numerical solutions to such boundary value problems rests o n the associated minimization principle; see [ 153,109] for details. 18.2. Separation of Variables. With conformal mapping no longer a viable option in three dim ensional space, sep- aration of variables reasserts its primacy for generating e xplicit solutions to the Laplace equation. As always, its applicability is unfortunately re stricted to rather special, but important, geometrical configurations. In three-dimensio nal space, the simplest separable cases are problems formulated on rectangular, cylindrical or spherical domains. Since the first two are straightforward extensions of their two-dimen sional counterparts, we will only discuss spherically separable solutions in any detail. The simplest domain to which the separation of variables method applies is a rectangular box : B=/braceleftbig 0<x<a, 0<y<b, 0<z <c/bracerightbig . For functions of three variables, one begins the separation process by splitting off one of them, by setting u(x,y,z) =v(x)w(y,z), say. The function v(x) satisfies a simple second order ordinary differential equation, while w(y,z) solves the two-dimensional Helmholtz equation, which is then separated by writing w(y,z) =p(y)q(z). The resulting fully separated solutions u(x,y,z) =v(x)p(y)q(z) are (mostly) products of trigonometric and hyperbolic functions. Complete details of the technique an d the resulting series solution are relegated to Exercise . 12/11/12 975 c/ci∇cleco√y∇t2012 Peter J. Olver In the case when the domain is a cylinder, one passes to cylind rical coordinates r,θ,z to effect the separation. The resulting separable solutions u(r,θ,z) =v(r,θ),w(z) = p(r)q(θ),w(z) are products of Bessel functions of the cylindrical radius r, trigonometric functions of the polar angle θ, and hyperbolic functions of z. Details are outlined in Exercise . The most interesting case is that of spherical coordinates, which we proceed to analyze in detail in the following subsection. Remark: Beyond these three well-known cases, there are, in fact, a t otal of eleven different coordinate systems in which the three-dimensiona l Laplace equation separates. See [131,134,136] for details on the more exotictypes of separation, includi ng ellipsoidal, toroidal, and parabolic spheroidal coordinates. The resul ting separable solutions lead to new classes of special functions. Laplace’s Equation in a Ball Suppose a solid ball (e.g., the earth), is subject a specified steady temperature distri- bution on its spherical boundary. Our task is to determine th e equilibrium temperature within the ball. To simplify matters, we assume that the body is composed of an isotropic, homogeneous medium, and shall choose units in which its radi us equals 1. Then, to find the equilibrium temperature within the ball, we must solve t he Dirichlet boundary value problem ∂2u ∂x2+∂2u ∂y2+∂2u ∂z2= 0, x2+y2+z2<1, u(x,y,z) =h(x,y,z), x2+y2+z2= 1.(18.12) Problems in spherical geometries tend to simplify when re-e xpressed in terms of spherical coordinates r,ϕ,θ, as defined by the usual formulae x=rsinϕcosθ, y =rsinϕsinθ, z =rcosϕ. (18.13) Here 0≤θ <2πmeasures the azimuthal angle orlongitude , while 0 ≤ϕ≤πmeasures thezenith angle orlatitude. Warning : We use the mathematician’s convention for spherical coord inates. Physi- cists often interchange the notation for the azimuthal and z enith angles; see Example B.8 for a detailed discussion. In spherical coordinates, the Laplace equation for†u(r,ϕ,θ) takes the form ∆u=∂2u ∂r2+2 r∂u ∂r+1 r2∂2u ∂ϕ2+cosϕ r2sinϕ∂u ∂ϕ+1 r2sin2ϕ∂2u ∂θ2= 0. (18.14) This important formula is the final result of a fairly nasty ch ain rule computation, whose details are left to the motivated reader. (Set aside lots of p aper and keep an eraser handy!) †Warning : See Section 15.2 for our convention on rewriting functions in new coord inates. 12/11/12 976 c/ci∇cleco√y∇t2012 Peter J. Olver Toconstructseparablesolutionstothesphericalcoordina teform(18.14)oftheLaplace equation, we begin by separating off the radial part of the sol ution, setting u(r,ϕ,θ) =v(r)w(ϕ,θ). (18.15) Substituting this ansatz into (18.14), multiplying the res ulting equation through byr2 vw, and then placing all the terms involving ron one side yields 1 v/parenleftbigg r2d2v dr2+2rdv dr/parenrightbigg =−1 w∆Sw=µ, (18.16) whereµis the separation constant, and ∆Sw=∂2w ∂ϕ2+cosϕ sinϕ∂w ∂ϕ+1 sin2ϕ∂2w ∂θ2. (18.17) The second order differential operator ∆S, which contains only the angular components of the full Laplacian operator ∆, is of particular significan ce. It is known as the spherical Laplacian , and governs the equilibrium and dynamics of thin spherical shells, cf. Exam- ple 18.13. Returning to equation (18.16), our usual separation argume nt applies. The left hand side depends only on r, while the right hand side depends only on the angles ϕ,θ. This can only occur when both sides are equal to a common separatio n constant, denoted by µ. As a consequence, the radial component v(r) satisfies the ordinary differential equation r2v′′+2rv′−µv= 0, (18.18) which is of Euler type (7.51), and hence can be readily solved . We will put this equation aside for the time being, and concentrate our efforts on the mo re complicated part. The angular components in (18.16) assume the form ∆S[w]+µw=∂2w ∂ϕ2+cosϕ sinϕ∂w ∂ϕ+1 sin2ϕ∂2w ∂θ2+µw= 0. (18.19) This second order partial differential equation can be regar ded as the eigenvalue equation for the spherical Laplacian operator ∆S, and is known as the spherical Helmholtz equation . To solve it, we adopt a further separation of angular variabl es, w(ϕ,θ) =p(ϕ)q(θ), (18.20) which we substitute into (18.19). Dividing the result by the productw=pq, multiplying by sin2ϕ, and then rearranging terms, we are led to the separated syst em sin2ϕ pd2p dϕ2+cosϕsinϕ pdp dϕ+µsin2ϕ=−1 qd2q dθ2=ν. The left hand side depends only on the zenith coordinate ϕwhile the right hand side depends only on the azimuthal coordinate θ. Since these angles are independent, the only way this could hold is when the two sides equal a common separa tion constant, denoted 12/11/12 977 c/ci∇cleco√y∇t2012 Peter J. Olver byν. The spherical Helmholtz equation thereby splits into a pai r of ordinary differential equations sin2ϕd2p dϕ2+cosϕsinϕdp dϕ+(µsin2ϕ−ν)p= 0,d2q dθ2+νq= 0. The equation for q(θ) is easy to solve. As one circumnavigates the sphere from wes t to east, the azimuthal angle θincreases from 0 to 2 π, soq(θ) must be a 2 πperiodic function. Thus, q(θ) satisfies the well-studied periodic boundary value proble m treated, for instance, in (15.33). Up to constant multiple, non-zero periodic solutions occur only when the separation constant assumes one of the values ν=m2, wherem= 0,1,2,...is an integer, with q(θ) = cosmθor sinmθ, m = 0,1,2,.... (18.21) Each positive ν=m2>0 admits two linearly independent 2 πperiodic solutions, while whenν= 0, only the constant solutions are periodic. With this information, we endeavor to solve the zenith equat ion sin2ϕd2p dϕ2+cosϕsinϕdp dϕ+(µsin2ϕ−m2)p= 0. (18.22) This is not easy, and constructing analytic formulas for its solutions requires some effort. The motivation behind the following steps will not be immedi ately apparent to the reader, since they are the result of a long, detailed study of this imp ortant differential equation by mathematicians over the last 200 years. As an initial simplification, we will eliminate the trigonom etric functions. To this end, we invoke the change of variables t= cosϕ, with p(ϕ) =P(cosϕ) =P(t). (18.23) Since 0≤ϕ≤π, we have 0 ≤/radicalbig 1−t2= sinϕ≤1. According to the chain rule, dp dϕ=−sinϕdP dt=−/radicalbig 1−t2dP dt, d2p dϕ2= sin2ϕd2P dt2−cosϕdP dt= (1−t2)d2P dt2−tdP dt. Substituting these expressions into (18.22), we conclude t hatP(t) must satisfy (1−t2)2d2P dt2−2t(1−t2)dP dt+/bracketleftbig µ(1−t2)−m2/bracketrightbig P= 0. (18.24) Unfortunately, the resulting differential equation is stil l not so easy to solve, but at least its coefficients are polynomials. Equation (18.24) is known a s theLegendre differential equation of orderm, and its solutions are known as Legendre functions , having first been employed by Legendre to study the gravitational attraction of ellipsoidal bodies. 12/11/12 978 c/ci∇cleco√y∇t2012 Peter J. Olver Power series solutions to the Legendre equation can be const ructed by the standard techniques presented in Appendix C. The most general soluti on is a new type of special function, known as a Legendre function , [3,145]. However, the solutions we are actually interested in can all be written in terms of elementary algeb raic functions. First of all, sincet= cosϕ, the solution only needs to be defined on the interval −1≤t≤1. The endpoints of this interval, t=±1, correspond to the sphere’s north pole, ϕ= 0, and south pole, ϕ=π. Both endpoints are singular points for the Legendre equati on since the coefficient (1 −t2)2of the leading order derivative vanishes when t=±1. In fact, both are regular singular points, as shown in Exercise C.3.4 4. Since ultimately we need the separable solution (18.15) to be a well-defined function ofx,y,z(even at points where the spherical coordinates degenerate, i.e., on the zaxis), we need p(ϕ) to be well-defined atϕ= 0 andπ, and this requires P(t) to be bounded at the singular points: |P(−1)|<∞, |P(+1)|<∞. (18.25) The combined boundary value problem (18.24–25) takes the fo rm of an eigensystem, in which the separation constant µis the eigenvalue and the non-zero solutions P(t)/ne}ationslash≡0 are the associated eigenfunctions. Mathematical justifications of the following statements ca n be found in Appendix C. Consider first the case m= 0, which assumes the simpler form (1−t2)d2P dt2−2tdP dt+µP= 0. (18.26) In this case, it turns out that the eigenfunctions, i.e., sol utions to the Legendre boundary value problem (18.26,25), are the Legendre polynomials Pn(t) =1 2nn!dn dtn(t2−1)n(18.27) thatfirst aroseinChapter 5 asoursimplest exampleoforthog onalpolynomials. Indeed, we can now finally comprehend the reason for the orthogonality o f the Legendre polynomials: they are the common eigenfunctions of a self-adjoint bounda ry value problem! Explicit formulas for the first few Legendre polynomials appear in (5. 46). Whenm>0, the eigenfunctions of the Legendre boundary value proble m (18.24–25) are not always polynomials. They are known as the associated Legendre functions , and have the explicit formula† Pm n(t) = (1−t2)m/2dm dtmPn(t) = (−1)n(1−t2)m/2 2nn!dn+m dtn+m(1−t2)n,n=m,m+1,..., (18.28) †Warning : Some authors include ( −1)min the formula, resulting in the opposite sign when m is odd. Another source of confusion is that many tables define the associ ated Legendre functions using the alternative initial factor ( t2−1)m/2. But we are solely interested in values of tlying between −1≤t≤1, which would result in complex values for odd m. 12/11/12 979 c/ci∇cleco√y∇t2012 Peter J. Olver which generalizes the Rodrigues formula (5.48) for the clas sical Legendre polynomials. Its proof is similar, and done in Exercise . Here is a list of the first few Legendre polynomials and associated Legendre functions: P0 0(t) = 1, P0 1(t) =t, P1 1(t) =/radicalbig 1−t2, P0 2(t) =−1 2+3 2t2, P1 2(t) = 3t/radicalbig 1−t2, P2 2(t) = 3(1−t2), P0 3(t) =−3 2t+5 2t3, P1 3(t) =/parenleftbig −3 2+15 2t2/parenrightbig/radicalbig 1−t2, P2 3(t) = 15t(1−t2), P3 3(t) = 15(1 −t2)3/2, (18.29) P0 4(t) =3 8−15 4t2+35 8t4, P1 4(t) =/parenleftbig −15 2+35 2t2/parenrightbig/radicalbig 1−t2, P2 4(t) =/parenleftbig −15 2+105 2t2/parenrightbig (1−t2), P3 4(t) = 105t(1−t2)3/2, P4 4(t) = 105(1 −t2)2. Whenm= 2k≤nis an even integer, Pm n(t) is a polynomial function, while when m= 2k+1≤nis odd, there is an extra factor of√ 1−t2. Keep in mind that the square root is real and positive since we are restricting our attention t o the interval −1≤t≤1. If m>n, formula (18.28) reduces to the zero function, and is not nee ded in the final tally. Warning : Even though half of the associated Legendre functions are p olynomials, only those with m= 0, i.e.,Pn(t) =P0 n(t), are called Legendre polynomials . Graphs of the first few Legendre polynomials can be found in Fi gure 5.4. In addition, Figure 18.1 displays the graphs of the associated Legendre f unctionsPm n(t) for 1≤m≤ n≤4. Pay particular attention to the fact that, owing to the cho ice of normalization factor, their graphs have very different vertical scales. The following result states that the Legendre polynomials a nd associated Legendre functions are a complete list of solutions to the Legendre bo undary value problem (18.24– 25). A proof can be found in [ 26]. Theorem 18.2. Letm≥0be a non-negative integer. Then the mthorder Legendre boundary value problem prescribed by (18.24–25)has eigenvalues µn=n(n+1)forn= 0,1,2,..., and associated eigenfunctions Pm n(t)wherem= 0,...,n. Returning totheoriginalvariable ϕvia(18.23),Theorem 18.2impliesthatouroriginal boundary value problem sin2ϕd2p dϕ2+cosϕsinϕdp dϕ+(µsin2ϕ−m2)p= 0,|p(0)|,|p(π)|<∞,(18.30) has its eigenvalues and eigenfunctions expressed in terms o f the Legendre functions: µn=n(n+1), pm n(ϕ) =Pm n(cosϕ),for 0 ≤m≤n. (18.31) The eigenfunction pm n(ϕ) is, in fact, a trigonometric polynomial of degree n; here are the 12/11/12 980 c/ci∇cleco√y∇t2012 Peter J. Olver -1 -0.5 0.5 10.20.40.60.81 P1 1(t)-1 -0.5 0.5 1 -1.5-1-0.50.511.5 P1 2(t)-1 -0.5 0.5 10.511.522.53 P2 2(t) -1 -0.5 0.5 1 -1.5-1-0.50.511.52 P1 3(t)-1 -0.5 0.5 1 -6-4-2246 P2 3(t)-1 -0.5 0.5 12468101214 P3 3(t) -1 -0.5 0.5 1 -2-112 P1 4(t)-1 -0.5 0.5 1 -7.5-5-2.52.557.510 P2 4(t)-1 -0.5 0.5 1 -30-20-10102030 P3 4(t)-1 -0.5 0.5 120406080100 P4 4(t) Figure 18.1. Associated Legendre Functions. first few, written in Fourier form: p0 0(ϕ) = 1, p0 1(ϕ) = cosϕ, p1 1(ϕ) = sinϕ, p0 2(ϕ) =1 4+3 4cos2ϕ, p1 2(ϕ) =3 2sin2ϕ, p2 2(ϕ) =3 2−3 2cos2ϕ, p0 3(ϕ) =3 8cosϕ+5 8cos3ϕ, p1 3(ϕ) =3 8sinϕ+15 8sin3ϕ, p2 3(ϕ) =15 4cosϕ−15 4cos3ϕ, p3 3(ϕ) =45 4sinϕ−15 4sin3ϕ, p0 4(ϕ) =9 64+5 16cos2ϕ+35 64cos4ϕ, p1 4(ϕ) =5 8sin2ϕ+35 16sin4ϕ, p2 4(ϕ) =45 16+15 4cos2ϕ−105 16cos4ϕ, p3 4(ϕ) =105 4sin2ϕ−105 8sin4ϕ, p4 4(ϕ) =315 8−105 2cos2ϕ+105 8cos4ϕ.(18.32) It is also instructive to plot the eigenfunctions in terms of the zenith angle ϕ; see Fig- ure 18.2. As in Figure 18.1, the vertical scales are not the sa me. At this stage, we have determined both angular components of our separable solutions 12/11/12 981 c/ci∇cleco√y∇t2012 Peter J. Olver 0.5 11.5 22.5 30.511.52 p0 0(ϕ)0.5 11.5 22.5 3 -1-0.50.51 p0 1(ϕ)0.5 11.5 22.5 30.20.40.60.81 p1 1(ϕ) 0.5 11.5 22.5 3 -0.4-0.20.20.40.60.81 p0 2(ϕ)0.5 11.5 22.5 3 -1.5-1-0.50.511.5 p1 2(ϕ)0.5 11.5 22.5 30.511.522.53 p2 2(ϕ)0.5 11.5 22.5 3 -1-0.50.51 p0 3(ϕ) 0.5 11.5 22.5 3 -1.5-1-0.50.511.52 p1 3(ϕ)0.5 11.5 22.5 3 -6-4-2246 p2 3(ϕ)0.5 11.5 22.5 32468101214 p3 3(ϕ)0.5 11.5 22.5 3 -0.4-0.20.20.40.60.81 p0 4(ϕ) 0.5 11.5 22.5 3 -2-112 p1 4(ϕ)0.5 11.5 22.5 3 -7.5-5-2.52.557.510 p2 4(ϕ)0.5 11.5 22.5 3 -30-20-10102030 p3 4(ϕ)0.5 11.5 22.5 320406080100 p4 4(ϕ) Figure 18.2. Trigonometric Legendre Functions. (18.20). Multiplying the two parts together results in the s pherical angle functions Ym n(ϕ,θ) =Pm n(cosϕ)cosmθ, /tildewideYm n(ϕ,θ) =Pm n(cosϕ)sinmθ,n= 0,1,2,..., m= 0,1,...,n,(18.33) known as spherical harmonics . They satisfy the spherical Helmholtz equation ∆SYm n+n(n+1)Ym n= 0 = ∆S/tildewideYm n+n(n+1)/tildewideYm n, (18.34) and so are eigenfunctions for the spherical Laplacian opera tor, (18.17), with associated 12/11/12 982 c/ci∇cleco√y∇t2012 Peter J. Olver eigenvalues µn=n(n+1) forn= 0,1,2,.... Thentheigenvalueµnadmits a (2 n+1)– dimensional eigenspace, spanned by the spherical harmonic s Y0 n(ϕ,θ), Y1 n(ϕ,θ), ... Yn n(ϕ,θ),/tildewideY1 n(ϕ,θ), .../tildewideYn n(ϕ,θ). (The omitted function /tildewideY0 n(ϕ,θ)≡0 is trivial, and so does not contribute.) In Figure 18.3 we plot the first few spherical harmonic surfaces r=Ym n(ϕ,θ). In these graphs, in view of the spherical coordinate formula (18.13), points with a neg ative value of rappear on the opposite side of the origin from the point on the unit sphere w ith anglesϕ,θ. Incidentally, the graphs of their counterparts r=/tildewideYm n(ϕ,θ), whenm/ne}ationslash= 0, are obtained by rotation around the zaxis by 90◦. On the other hand, the graphs of Y0 nare cylindrically symmetric (why?), and hence unaffected by such a rotation. Self-adjointness of the spherical Laplacian, cf. Exercise , implies that the spherical harmonics are orthogonal with respect to the L2inner product /an}b∇acketle{tf;g/an}b∇acket∇i}ht=/integraldisplay/integraldisplay S1fgdS=/integraldisplayπ 0/integraldisplay2π 0f(ϕ,θ)g(ϕ,θ)sinϕ dθdϕ (18.35) given by integrating the product of the functions with respe ct to surface area over the unit sphere S1={/ba∇dblx/ba∇dbl= 1}, cf. (B.42). More correctly, self-adjointness only guaran - tees orthogonality for the harmonics corresponding to dist inct eigenvalues. However, the orthogonality relations /an}b∇acketle{tYm n;Yk l/an}b∇acket∇i}ht=/integraldisplay/integraldisplay S1Ym nYk ldS= 0, (m,n)/ne}ationslash= (k,l), /an}b∇acketle{tYm n;/tildewideYk l/an}b∇acket∇i}ht=/integraldisplay/integraldisplay S1Ym n/tildewideYk ldS= 0,for all ( m,n),(k,l), /an}b∇acketle{t/tildewideYm n;/tildewideYk l/an}b∇acket∇i}ht=/integraldisplay/integraldisplay S1/tildewideYm n/tildewideYk ldS= 0 ( m,n)/ne}ationslash= (k,l),(18.36) do, in fact, hold in full generality; Exercise contains the details. Their norms can be explicitly computed: /ba∇dblY0 n/ba∇dbl2=4π 2n+1, /ba∇dblYm n/ba∇dbl2=/ba∇dbl/tildewideYm n/ba∇dbl2=2π(n+m)! (2n+1)(n−m)!.(18.37) A proof of this formula appears in Exercise . Withsomefurther work, itcanbeshown thatthespherical har monicsform a complete orthogonalsystemoffunctionsontheunitsphere. Thismean sthatanyreasonablefunction h:S1→R, e.g., piecewise C1, can be expanded into a convergent spherical Fourier series h(ϕ,θ) =c0,0 2+∞/summationdisplay n=1/parenleftBigg c0,n 2Y0 n(ϕ)+n/summationdisplay m=1/bracketleftBig cm,nYm n(ϕ,θ)+/tildewidecm,n/tildewideYm n(ϕ,θ)/bracketrightBig/parenrightBigg (18.38) in the spherical harmonics. Applying the orthogonality rel ations (18.36), we find that the spherical Fourier coefficients are given by the inner product s c0,n=2/an}b∇acketle{th;Y0 n/an}b∇acket∇i}ht /ba∇dblY0 n/ba∇dbl2, cm,n=/an}b∇acketle{th;Ym n/an}b∇acket∇i}ht /ba∇dblYm n/ba∇dbl2,/tildewidecm,n=/an}b∇acketle{th;/tildewideYm n/an}b∇acket∇i}ht /ba∇dbl/tildewideYm n/ba∇dbl2,0≤n, 1≤m≤n, 12/11/12 983 c/ci∇cleco√y∇t2012 Peter J. Olver Y0 0(ϕ,θ) Y0 1(ϕ,θ) Y1 1(ϕ,θ) Y0 2(ϕ,θ) Y1 2(ϕ,θ) Y2 2(ϕ,θ) Y0 3(ϕ,θ) Y1 3(ϕ,θ) Y2 3(ϕ,θ) Y3 3(ϕ,θ) Y0 4(ϕ,θ) Y1 4(ϕ,θ) Y2 4(ϕ,θ) Y3 4(ϕ,θ) Y4 4(ϕ,θ) Figure 18.3. Spherical Harmonics. 12/11/12 984 c/ci∇cleco√y∇t2012 Peter J. Olver or, explicitly, using (18.35) and the formulae (18.37) for t he norms, cm,n=(2n+1)(n−m)! 2π(n+m)!/integraldisplay2π 0/integraldisplayπ 0h(ϕ,θ)Pm n(cosϕ) cosmθsinϕ dϕdθ, /tildewidecm,n=(2n+1)(n−m)! 2π(n+m)!/integraldisplay2π 0/integraldisplayπ 0h(ϕ,θ)Pm n(cosϕ) sinmθsinϕ dϕdθ.(18.39) As with an ordinary Fourier series, the extra1 2was appended to the c0,nterms in the series (18.38) so that the formulae (18.39) are valid for all m,n. In particular, the constant term in the spherical harmonic series is the mean of the function hover the unit sphere: c0,0 2=1 4π/integraldisplay/integraldisplay S1hdS=1 4π/integraldisplay2π 0/integraldisplayπ 0h(ϕ,θ) sinϕdϕdθ (18.40) Remark: An alternative approach is to replace the real trigonometr ic functions by complex exponentials, and work with the complex spherical harmonics Ym n(θ,ϕ) =Ym n(θ,ϕ)+ i/tildewideYm n(θ,ϕ) =Pm n(cosϕ)eimθ,n= 0,1,2,... , m=−n,−n+1,...,n.(18.41) The complex orthogonality and expansion formulas are releg ated to the exercises. To complete our solution to the Laplace equation on the solid ball, we still need to analyze the ordinary differential equation (18.18) for the r adial component v(r). In view of our analysis of the spherical Helmholtz equation, the ori ginal separation constant is µ=n(n+ 1) for some non-negative integer n≥0, and so the radial equation takes the form r2v′′+2rv′−n(n+1)v= 0. (18.42) As noted earlier, to solve such a second order linear equatio n of Euler type (3.84), we substitutethepower ansatz v(r) =rα. Theexponent αmust satisfythequadraticequation α2+α−n(n+1) = 0,and hence α=norα=−(n+1). Therefore, the two linearly independent solutions are v1(r) =rnandv2(r) =r−n−1. (18.43) Since here we are only interested in solutions that remain bo unded atr= 0 — the center of the ball — we should just retain the first solution v(r) =rnin our subsequent analysis. At this stage, we have solved all three ordinary differential equations for the sepa- rable solutions. We combine the results (18.21,33,43) toge ther to produce the following spherically separable solutions to the Laplace equation: Hm n=rnYm n(ϕ,θ) =rnPm n(cosϕ)cosmθ /tildewideHm n=rn/tildewideYm n(ϕ,θ) =rnPm n(cosϕ)sinmθn= 0,1,2,..., m= 0,1,...,n.(18.44) 12/11/12 985 c/ci∇cleco√y∇t2012 Peter J. Olver Although apparently complicated, these solutions are, sur prisingly, elementary polynomial functions of the rectangular coordinates x,y,z, and hence harmonic polynomials . The first few are H0 0= 1, H0 1=z, H0 2=z2−1 2x2−1 2y2H0 3=z3−3 2x2z−3 2y2z H1 1=x, H1 2= 3xz, H1 3= 6xz2−3 2x3−3 2xy2 /tildewideH1 1=y,/tildewideH1 2= 3yz, /tildewideH1 3= 6yz2−3 2x2y−3 2y3 H2 2= 3x2−3y2, H2 3= 15x2z−15y2z /tildewideH2 2= 6xy, /tildewideH2 3= 30xyz H3 3= 15x3−45xy2 /tildewideH3 3= 45x2y−15y3.(18.45) The polynomials H0 n, H1 n, ... , Hn n,/tildewideH1 n, .../tildewideHn n form a basis for the vector space H(n)of all homogeneous harmonic polynomials of degree n, which therefore has dimension 2 n+1. The harmonic polynomials form a complete orthogonal system , and therefore the general solution to the Laplace equation inside the unit bal l can be written as a harmonic polynomial series: u(x,y,z) =c0,0 2+∞/summationdisplay n=1/parenleftBigg c0,n 2H0 n(x,y,z)+n/summationdisplay m=1/bracketleftBig cm,nHm n(x,y,z)+/tildewidecm,n/tildewideHm n(x,y,z)/bracketrightBig/parenrightBigg , (18.46) or, equivalently in spherical coordinates, u(r,ϕ,θ) =c0,0 2+∞/summationdisplay n=1/parenleftBigg c0,n 2rnY0 n(ϕ)+n/summationdisplay m=1/bracketleftBig cm,nrnYm n(ϕ,θ)+/tildewidecm,nrn/tildewideYm n(ϕ,θ)/bracketrightBig/parenrightBigg . (18.47) The coefficients cm,n,/tildewidecm,nare uniquely prescribed by the boundary conditions. Indeed , substituting (18.47) into the Dirichlet boundary conditio ns on the unit sphere r= 1 yields u(1,ϕ,θ) =c0,0 2+∞/summationdisplay n=1/parenleftBigg c0,n 2Y0 n(ϕ)+n/summationdisplay m=1/bracketleftBig cm,nYm n(ϕ,θ)+/tildewidecm,n/tildewideYm n(ϕ,θ)/bracketrightBig/parenrightBigg =h(ϕ,θ). (18.48) Thus, the coefficients cm,n,/tildewidecm,nare given by the orthogonality formulae (18.39). If the terms in the resulting series are uniformly bounded — which o ccurs for all integrable functionsh, as well as also certain generalized functions including th e delta function — then the harmonic polynomial series (18.47) converges ever ywhere, and, in fact, uniformly on any smaller ball /ba∇dblx/ba∇dbl=r≤r0<1. In rectangular coordinates, the nthsummand of the series (18.46) is a homogeneous polynomial of degree n. Therefore, repeating the argument used on the two-dimensi onal polar coordinate solution (15.39), we conclude that the har monic polynomial series is, 12/11/12 986 c/ci∇cleco√y∇t2012 Peter J. Olver in fact, a power series, and hence hence provides the Taylor expansion for the harmonic functionu(x,y,z)at the origin ! In particular, this implies that the harmonic function u(x,y,z) is analytic at 0. The constant term in such a Taylor series can be identified wit h the value of the function at the origin: u(0,0,0) =1 2c0,0. On the other hand, since u=honS1=∂Ω, the coefficient formula (18.40) tells us that u(0,0,0) =c0,0 2=1 4π/integraldisplay/integraldisplay S1udS. (18.49) Therefore, we have established the three-dimensional coun terpart of Theorem 15.8: the value of the harmonic function at the center of the sphere is e qual to the average of its values on the sphere’s surface. In addition, the higher orde r coefficients cm,n,/tildewidecm,nserve to prescribe the partial derivatives∂i+j+ku ∂xi∂yj∂zk(0,0,0). In this way, the orthogonality formulae (18.39) can be re-interpreted as three-dimension al counterparts of the Cauchy formulae (16.139) for the derivatives of the real and imagin ary parts of a complex analytic function; see Exercise for details. So far, we have restricted our attention to the sphere of unit radius. A simple scaling argument serves to establish the general result. Theorem 18.3. Ifu(x)is a harmonic function defined on a domain Ω⊂R3, thenu is analytic inside Ω. Moreover, its value at any point x0∈Ωis obtained by averaging its values on any sphere centered at x0, so u(x0) =1 4πa2/integraldisplay/integraldisplay /bardblx−x0/bardbl=au dS, (18.50) provided the enclosed ball {/ba∇dblx−x0/ba∇dbl ≤a} ⊂Ωlies entirely within the domain of defini- tion. Proof: It is easily checked that, under the hypothesis of the theor em, the rescaled and translated function U(y) =u(ay+x0) =u(x),where y=x−x0 a, (18.51) isharmoniconthe unit ball /ba∇dbly/ba∇dbl ≤1, and hence solvesthe boundary valueproblem (18.12) with boundary values h(y) =U(y) =u(ay+x0) on/ba∇dbly/ba∇dbl= 1. By the preceding remarks, U(y) is analytic at y=0, and sou(x) =U/parenleftBigx−x0a/parenrightBig is analytic at x=x0. Sincex0can be any point inside Ω, this proves analyticity of ueverywhere in Ω. Moreover, according to the integral formula (18.49), u(x0) =U(0) =1 4π/integraldisplay/integraldisplay /bardbly/bardbl=1UdS=1 4πa2/integraldisplay/integraldisplay /bardblx−x0/bardbl=au dS, since the change of variables (18.51) has the effect of rescal ing the spherical surface inte- gral. Q.E.D. 12/11/12 987 c/ci∇cleco√y∇t2012 Peter J. Olver Arguing as in the planar case of Theorem 15.9, we readily esta blish the corresponding Maximum Principle for harmonic functions of three variable s. Theorem 18.4. A non-constant harmonic function cannot have a local maximu m or minimum at any interior point of its domain of definition. Mor eover, its global maximum or minimum (if any)can only occur on the boundary of the domain. For instance, the Maximum Principle implies that the maximu m and minimum tem- peratures in a solid body in thermal equilibrium are to be fou nd only on its boundary. In physical terms, since heat energy must flow away from any inte rnal maximum and towards any internal minimum, any local temperature extremum insid e the body would preclude it from being in thermal equilibrium. Example 18.5. In this example, we shall determine the electrostatic poten tial inside a hollow sphere when the upper and lower hemispheres are held at different constant potentials. This device is called a spherical capacitor and is realized experimentally by separating the two charged conducting hemispherical shell s by a thin insulating ring at the equator. A straightforward scaling argument allows us t o choose our units so that the sphere has radius 1, while the potential is set equal to 1 on th e upper hemisphere and equal to 0, i.e., grounded, on the lower hemisphere. The resu lting electrostatic potential ssatisfies the Laplace equation ∆u= 0 inside a solid ball /ba∇dblx/ba∇dbl<1, and is subject to Dirichlet boundary conditions u(x,y,z) =/braceleftbigg1, z> 0, 0, z< 0,on the unit sphere /ba∇dblx/ba∇dbl= 1. (18.52) The solution will be prescribed by a harmonic polynomial ser ies (18.46) whose coeffi- cientsarefixedby theboundaryvalues(18.52). Beforetakin gontherequiredcomputation, let us first note that since the boundary data does not depend u pon the azimuthal angle θ, the solution u=u(r,ϕ) will also be independent of θ. Therefore, we need only consider theθ-independent spherical harmonic polynomials (18.33), whi ch are those with m= 0, and hence u(r,ϕ) =1 2∞/summationdisplay n=0cnH0 n(x,y,z) =1 2∞/summationdisplay n=0cnrnPn(cosϕ), where we abbreviate cn=c0,n. The boundary conditions (18.52) require u(1,ϕ) =1 2∞/summationdisplay n=0cnPn(cosϕ) =h(ϕ) =/braceleftBigg 1,0≤ϕ<1 2π, 0,1 2π<ϕ≤π. The coefficients are given by (18.39), which, in the case m= 0, reduce to cn=2n+1 2π/integraldisplay/integraldisplay S1f Y0 ndS= (2n+1)/integraldisplayπ/2 0Pn(cosϕ)sinϕ dϕ= (2n+1)/integraldisplay1 0Pn(t)dt, (18.53) 12/11/12 988 c/ci∇cleco√y∇t2012 Peter J. Olver sincef= 0 when1 2π<ϕ≤π. The first few are c0= 1, c1=3 2, c2= 0, c3=−7 8, c4= 0, ... . Therefore, the solution has the explicit Taylor expansion u(x,y,z) =1 2+3 4rcosϕ−21 128r3cosϕ−35 128r3cos3ϕ+··· =1 2+3 4z+21 32(x2+y2)z−7 16z3+···.(18.54) Note in particular that the value u(0,0,0) =1 2at the center of the sphere is the average of its boundary values, in accordance with Theorem 18.3. Obs erve that the solution only depends upon the cylindrical coordinates r,z. This follows from the invariance of the Laplace equation under general rotations, coupled with the invariance of the boundary data under rotations around the zaxis. Remark: The same solution u(x,y,z) describes the thermal equilibrium in a solid sphere whose upper hemisphere is held at temperature 1◦and lower hemisphere at 0◦. Example 18.6. A closely related problem is to determine the electrostatic potential outsidea spherical capacitor. As in the preceding example, we take o ur capacitor of radius 1, with electrostatic charge of 1 on the upper hemisphere and 0 on the lower hemisphere. Here, we need to solve the Laplace equation ∆ u= 0 in the unbounded domain Ω = {/ba∇dblx/ba∇dbl>1}— the exterior of the unit sphere — subject to same Dirichlet b oundary conditions (18.52). We anticipate that the potential will b e vanishingly small at large distances away from the capacitor: r=/ba∇dblx/ba∇dbl ≫1. Therefore, the harmonic polynomial solutions (18.44) will not help us solve this problem, since (except for the constant case) they become unboundedly large far away from the origin. However, reconsideration of our original separation of var iables argument will produce a different class ofsolutionshaving thedesired decay prope rties. When we solvedtheradial equation (18.42), we discarded the solution v2(r) =r−n−1because it had a singularity at the origin. In the present situation, the behavior of the fun ction atr= 0 is irrelevant; our requirement is that the solution decays as r→ ∞, andv2(r) has this property. Therefore, we will utilize the complementary harmonic functions Km n(x,y,z) =r−2n−1Hm n(x,y,z) =r−n−1Ym n(ϕ,θ) =r−n−1Pm n(cosϕ) cosmθ, /tildewideKm n(x,y,z) =r−2n−1/tildewideHm n(x,y,z) =r−n−1/tildewideYm n(ϕ,θ) =r−n−1Pm n(cosϕ) sinmθ,(18.55) for solving such exterior problems. For the capacitor probl em, we need only those that are independent of θ, which have m= 0. We write the resulting solution as a series u(r,ϕ) =1 2∞/summationdisplay n=0cnK0 n(x,y,z) =1 2∞/summationdisplay n=0cnr−n−1Pn(cosϕ). The boundary conditions u(1,ϕ) =1 2∞/summationdisplay n=0cnPn(cosϕ) =f(ϕ) =/braceleftBigg 1,0≤ϕ<1 2π, 0,1 2π<ϕ≤π, 12/11/12 989 c/ci∇cleco√y∇t2012 Peter J. Olver are identical with those in the previous example. Therefore , the coefficients are given by (18.53), leading to the series expansion u=1 2r+3cosϕ 4r2−21cosϕ+35cos3ϕ 128r4+···=1 2r+3z 4r3+21(x2+y2)z−14z3 32r7+···, (18.56) wherer=/radicalbig x2+y2+z2. Interestingly, at large distances, the higher order terms become negligible, and the potential looks like that associated wi th a point charge of magnitude1 2 — the average of the potential over the sphere — that is concen trated at the origin. This is indicative of a general fact, to be explored in Exercise . 18.3. The Green’s Function. We now turn to the inhomogeneous form of the three-dimension al Laplace equation: thePoisson equation −∆u=ffor all x∈Ω (18 .57) on a solid domain Ω ⊂R3. In order to uniquely specify the solution, we must impose appropriate boundary conditions: Dirichlet or mixed. We on ly need to discuss the case of homogeneous boundary conditions, since, by linearity, an i nhomogeneous boundary value problemcanbesplitupintoahomogeneousboundaryvaluepro blemfortheinhomogeneous Poisson equation and an inhomogeneous boundary value probl em for the homogeneous Laplace equation. As in Chapters 11 and 15, we begin by analyzing the case of a del ta function inhomo- geneity that is concentrated at a single point in the domain. Thus, for each ξ= (ξ,η,ζ)∈ Ω, theGreen’s function G(x;ξ) =G(x,y,z;ξ,η,ζ) is the unique solution to the Poisson equation −∆u=δ(x−ξ) =δ(x−ξ)δ(y−η)δ(z−ζ) for all x∈Ω, (18.58) subject to the chosen homogeneous boundary conditions. The solution to the general Poisson equation (18.57) is then obtained by superposition : We write the forcing function f(x,y,z) =/integraldisplay/integraldisplay/integraldisplay Ωf(ξ,η,ζ)δ(x−ξ)δ(y−η)δ(z−ζ)dξdηdζ as a linear superposition of delta functions. By linearity, the solution u(x,y,z) =/integraldisplay/integraldisplay/integraldisplay Ωf(ξ,η,ζ)G(x,y,z;ξ,η,ζ)dξdηdζ (18.59) to the homogeneous boundary value problem for the Poisson eq uation (18.57) is then given as the corresponding superposition of the Green’s function solutions. The Green’s Function in Space Only in a few specific instances is the explicit formula for th e Green’s function known. Nevertheless, certain general guiding features can be read ily established. The starting point is to investigate the Poisson equation (18.58) when th e domain Ω = R3is all of 12/11/12 990 c/ci∇cleco√y∇t2012 Peter J. Olver three-dimensional space. We impose boundary constraints b y seeking a solution that goes to zero,u(x)→0, at large distances /ba∇dblx/ba∇dbl → ∞. Since the Laplacian is invariant under translations we can, without loss of generality, place our d elta impulse at the origin, and concentrate on solving the particular case −∆u=δ(x), x∈R3. Sinceδ(x) = 0 for all x/ne}ationslash=0, the desired solution will, in fact, be a solution to the homogeneous Laplace equation ∆u= 0, x/ne}ationslash=0, save, possibly, for a singularity at the origin. The Laplace equation models the equilibria of a homogeneous , isotropic medium, and so is also invariant under three-dimensional rotations; de tails can be found in Exercise . This suggests that, in any radially symmetric configuratio n, the solution should only depend upon the distance r=/ba∇dblx/ba∇dblfrom the origin. Referring to the spherical coordinate form (18.14) of the Laplacian operator, if uonly depends upon r, its derivatives with respect to the angular coordinates ϕ,θare zero, and so u(r) solves the ordinary differential equation d2u dr2+2 rdu dr= 0. (18.60) This equation is, in effect, a first order linear ordinary diffe rential equation for v=du/dr and hence is particularly easy to solve: du dr=v(r) =−b r2,and hence u(r) =a+b r, wherea,bare arbitrary constants. The constant solution u(r) =adoes not die away at large distances, nor does it have a singularity at the origin . Therefore, if our intuition is valid, the desired solution should be of the form u=b r=b /ba∇dblx/ba∇dbl=b/radicalbig x2+y2+z2. (18.61) Indeed, this function is harmonic — solves Laplace’s equati on — everywhere away from the origin, and has a singularity at x=0. The solution (18.61) is, up to constant multiple, the three- dimensional Newtonian gravitational potential due to a point mass at the origin. It s gradient g(x) =∇/parenleftbiggb /ba∇dblx/ba∇dbl/parenrightbigg =−bx /ba∇dblx/ba∇dbl3. (18.62) defines the gravitational force vector at the point x. Whenb >0, the force g(x) points towards the mass at the origin. Its magnitude /ba∇dblg/ba∇dbl=b /ba∇dblx/ba∇dbl2=b r2 12/11/12 991 c/ci∇cleco√y∇t2012 Peter J. Olver is proportional to one over the squared distance, which is th e well-known inverse square law of three-dimensional Newtonian gravity. Thus, (18.61) can also be interpreted as the electrostatic Coulomb potential on a charged mass at positi onxdue to a concentrated electric charge at the origin, with (18.62) the correspondi ng electrostatic force. The con- stantbis positive when the charges are of opposite signs, leading t o an attractive force, and negative in the repulsive case of like charges. Returning to our problem, the remaining task is to fix the mult iplebsuch that the Laplacian of our candidate solution (18.61) has a delta func tion singularity at the origin; equivalently, we must determine c= 1/bsuch that −∆r−1=cδ(x). (18.63) This equation is certainly valid away from the origin, since δ(x) = 0 when x/ne}ationslash=0. To investigate near the singularity, we integrate both sides o f (18.63) over a small solid ball Bε={/ba∇dblx/ba∇dbl ≤ε}of radiusε: −/integraldisplay/integraldisplay/integraldisplay Bε∆r−1dxdydz =/integraldisplay/integraldisplay/integraldisplay Bεcδ(x)dxdydz =c, (18.64) where we used the definition of the delta function to evaluate the right hand side. On the other hand, since ∆ r−1=∇·∇r−1, we can use the divergence theorem (B.85) to evaluate the left hand integral, whence /integraldisplay/integraldisplay/integraldisplay Bε∆r−1dxdydz =/integraldisplay/integraldisplay/integraldisplay Bε∇·∇r−1dxdydz =/integraldisplay/integraldisplay Sε∂ ∂n/parenleftbigg1 r/parenrightbigg dS, where the surface integral is over the bounding sphere Sε=∂Bε={/ba∇dblx/ba∇dbl=ε}. The sphere’s unit normal npoints in the radial direction, and hence the normal derivat ive coincides with differentiation with respect to r; in particular, ∂ ∂n/parenleftbigg1 r/parenrightbigg =∂ ∂r/parenleftbigg1 r/parenrightbigg =−1 r2. The surface integral can now be explicitly evaluated: /integraldisplay/integraldisplay Sε∂ ∂n/parenleftbigg1 r/parenrightbigg dS=−/integraldisplay/integraldisplay Sε1 r2dS=−/integraldisplay/integraldisplay Sε1 ε2dS=−4π, sinceSεhas surface area 4 πε2. Substituting this result back into (18.64), we conclude that c= 4π,and hence −∆r−1= 4πδ(x). (18.65) This is our desired formula! We conclude that the solution to Poisson’s equation for a delta function impulse at the origin is G(x,y,z) =1 4πr=1 4π/ba∇dblx/ba∇dbl=1 4π/radicalbig x2+y2+z2, (18.66) which is the three-dimensional Newtonian potential due to a unit point mass situated at the origin. 12/11/12 992 c/ci∇cleco√y∇t2012 Peter J. Olver If the singularity is concentrated at some other point ξ= (ξ,η,ζ), then we merely translate the preceding solution. This leads immediately t o the Green’s function G(x;ξ) =G(x−ξ) =1 4π/ba∇dblx−ξ/ba∇dbl=1 4π/radicalbig (x−ξ)2+(y−η)2+(z−ζ)2.(18.67) The superposition principle (18.59) implies the following integral formula for the solutions to the Poisson equation on all of three-dimensional space. Theorem 18.7. A particular solution to the Poisson equation −∆u=ffor x∈R3(18.68) is given by u⋆(x) =1 4π/integraldisplay/integraldisplay/integraldisplay R3f(ξ) /ba∇dblx−ξ/ba∇dbldξ=1 4π/integraldisplay/integraldisplay/integraldisplay R3f(ξ,η,ζ)dξdηdζ/radicalbig (x−ξ)2+(y−η)2+(z−ζ)2.(18.69) The general solution is u(x,y,z) =u⋆(x,y,z)+w(x,y,z), wherew(x,y,z)is an arbitrary harmonic function. Example 18.8. In this example, we compute the gravitational (or electrost atic) potential in three-dimensional space due to a uniform solid ball, e.g., a spherical planet such as the earth. By rescaling, it suffices to consider the cas e when the forcing function f(x) =/braceleftbigg1,/ba∇dblx/ba∇dbl<1, 0,/ba∇dblx/ba∇dbl>1, is equal to 1 inside a solid ball of radius 1 and zero outside. T he particular solution to the resulting Poisson equation (18.68) is given by the integral u(x) =1 4π/integraldisplay/integraldisplay/integraldisplay /bardblξ/bardbl<11 /ba∇dblx−ξ/ba∇dbldξdηdζ. (18.70) Clearly, since the forcing function is radially symmetric, the solution u=u(r) is also radially symmetric. To evaluate the integral, then, we can t akex= (0,0,z) to lie on the zaxis, so that r=/ba∇dblx/ba∇dbl=|z|. We use cylindrical coordinates ξ= (ρcosθ,ρsinθ,ζ), so that /ba∇dblx−ξ/ba∇dbl=/radicalbig ρ2+(z−ζ)2. The integral in (18.70) can then be explicitly computed: 1 4π/integraldisplay1 −1/integraldisplay√ 1−ζ2 0/integraldisplay2π 0ρdθdρdζ/radicalbig ρ2+(z−ζ)2= =1 2/integraldisplay1 −1/parenleftBig/radicalbig 1+z2−2zζ− |z−ζ|/parenrightBig dζ=  1 3|z|,|z| ≥1, 1 2−z2 6,|z| ≤1. 12/11/12 993 c/ci∇cleco√y∇t2012 Peter J. Olver 1 2 3 40.10.20.30.40.5 Figure 18.4. Solution to Poisson’s Equation in a Solid Ball. Therefore, by radial symmetry, the solution is u(x) =  1 3r, r =/ba∇dblx/ba∇dbl ≥1, 1 2−r2 6, r=/ba∇dblx/ba∇dbl ≤1,(18.71) plotted, as a function of r=/ba∇dblx/ba∇dblin Figure 18.4. Note that, outside the solid ball, the solutionisa Newtonianpotentialcorresponding to a concen trated point mass of magnitude 4 3π— the total mass of the planet. We have thus demonstrated a wel l-known result in gravitation and electrostatics: the exterior potential du e to a spherically symmetric mass (or electric charge) is the same as if all the mass (charge) we re concentrated at its center. In outer space if you can’t see a spherical planet, you can onl y determine its mass, not its size, by measuring its external gravitational force. Bounded Domains and the Method of Images Suppose we now wish to solve the inhomogeneous Poisson equat ion (18.57) on a bounded domainΩ ⊂R3. Toconstruct thedesired Green’s function, weproceed asfo llows. The Newtonian potential (18.67) is a particular solution to the underlying inhomogeneous equation −∆u=δ(x−ξ), x∈Ω, (18.72) but it almost surely does not have the proper boundary values on∂Ω. By linearity, the general solution to such an inhomogeneous linear equation i s of the form u(x) =1 4π/ba∇dblx−ξ/ba∇dbl−v(x), (18.73) where the first summand is a particular solution, which we now know, while†v(x) is an arbitrary solution to the homogeneous equation ∆ v= 0, i.e., an arbitrary harmonic function. The solution (18.73) satisfies the homogeneous bo undary conditions provided †The minus sign is for later convenience. 12/11/12 994 c/ci∇cleco√y∇t2012 Peter J. Olver ξηx Figure 18.5. Method of Images for the Unit Ball. the boundary values of v(x) match those of the Green’s function. Let us explicitly stat e the result in the Dirichlet case. Theorem 18.9. The Green’s functionfor thehomogeneous Dirichlet boundar y value problem −∆u=f, x∈Ω, u = 0,x∈∂Ω, for the Poisson equation in a domain Ω⊂R3has the form G(x;ξ) =1 4π/ba∇dblx−ξ/ba∇dbl−v(x;ξ),x,ξ∈Ω, (18.74) wherev(x;ξ)is the harmonic function of xthat satisfies v(x;ξ) =1 4π/ba∇dblx−ξ/ba∇dblfor all x∈∂Ω. In this manner, we have reduced the determination of the Gree n’s function to the solutiontoaparticularfamilyofLaplaceboundaryvaluepr oblems, whichareparametrized by the point ξ∈Ω. Incertaindomainswithsimplegeometry,theMethodofImage scanbeusedtoproduce an explicit formula for the Green’s function. As in Section 1 5.3, the idea is to match the boundary values of the free space Green’s function due to a de lta impulse at a point inside the domain with one or more additional Green’s functions cor responding to impulses at points outside the domain — the “image points”. The case of a solid ball of radius 1 with Dirichlet boundary co nditions is the easiest to handle. Indeed, the samegeometrical construction that we used for a planar disk, red rawn in Figure 18.5, applies here. Although the same as Figure 15. 8, we are re-interpreting it as a three-dimensional diagram, with the circle representing the unit sphere, while the lines remain lines. The required image point is given by inversion : η=ξ /ba∇dblξ/ba∇dbl2,whereby /ba∇dblξ/ba∇dbl=1 /ba∇dblη/ba∇dbl. 12/11/12 995 c/ci∇cleco√y∇t2012 Peter J. Olver By the similar triangles argument used before, we find /ba∇dblξ/ba∇dbl /ba∇dblx/ba∇dbl=/ba∇dblx/ba∇dbl /ba∇dblη/ba∇dbl=/ba∇dblx−ξ/ba∇dbl /ba∇dblx−η/ba∇dbl, and therefore /ba∇dblx/ba∇dbl= 1. As a result, the function v(x,ξ) =1 4π/ba∇dblη/ba∇dbl /ba∇dblx−η/ba∇dbl=1 4π/ba∇dblξ/ba∇dbl /ba∇dblξ−/ba∇dblξ/ba∇dbl2x/ba∇dbl has the same boundary values on the unit sphere as the Newtoni an potential: 1 4π/ba∇dblη/ba∇dbl /ba∇dblx−η/ba∇dbl=1 4π/ba∇dblx−ξ/ba∇dblwhenever /ba∇dblx/ba∇dbl= 1. We conclude that their difference G(x;ξ) =1 4π/parenleftbigg1 /ba∇dblx−ξ/ba∇dbl−/ba∇dblξ/ba∇dbl /ba∇dblξ−/ba∇dblξ/ba∇dbl2x/ba∇dbl/parenrightbigg (18.75) has the required properties of the Green’s function: it sati sfies the Laplace equation inside the unit ball except at the delta function singularity x=ξ, and, moreover, G(x;ξ) = 0 has homogeneous Dirichlet conditions on the spherical boun dary/ba∇dblx/ba∇dbl= 1. With the Green’s function in hand, we can apply the general su perposition for- mula (18.59)to arriveat a solutionto theDirichlet boundar y value problem for the Poisson equation in the unit ball. Theorem 18.10. The solution to the homogeneous Dirichlet boundary value pr ob- lem −∆u=f,for/ba∇dblx/ba∇dbl<1, u = 0,for/ba∇dblx/ba∇dbl= 1 is u(x) =1 4π/integraldisplay/integraldisplay/integraldisplay /bardblξ/bardbl≤1/parenleftbigg1 /ba∇dblx−ξ/ba∇dbl−/ba∇dblξ/ba∇dbl /ba∇dblξ−/ba∇dblξ/ba∇dbl2x/ba∇dbl/parenrightbigg f(ξ)dξdηdζ. (18.76) The Green’s function can also be used to solve the inhomogene ous boundary value problem −∆u= 0,x∈Ω, u =h,x∈∂Ω. (18.77) The same argument used in the two-dimensional situation pro duces the solution u(x) =−/integraldisplay/integraldisplay ∂Ω∂G(x;ξ) ∂nh(ξ)dS. (18.78) In the case when Ω is a solid ball, this integral formula effect ively sums the spherical harmonic series (18.46). 12/11/12 996 c/ci∇cleco√y∇t2012 Peter J. Olver 18.4. The Heat Equation in Three-Dimensional Media. Thermal diffusion in a homogeneous, isotropic solid body Ω ⊂R3is governed by the three-dimensional heat equation ∂u ∂t=γ∆u=γ/parenleftbigg∂2u ∂x2+∂2u ∂y2+∂2u ∂z2/parenrightbigg , (x,y,z)∈Ω.(18.79) The positivity of the body’s thermal diffusivity γ >0 is reaquired on both physical and mathematical grounds. The physical derivation is exactly t he same as the two-dimensional version (17.1), and does not need to be repeated in detail. Br iefly, the heat flux vector is proportional to the temperature gradient, w=−κ∇u, while its divergence is proportional to the rate of change of temperature: ∇·w=−σut. Combining these two physical laws and assuming homogeneity, whereby κandσare constant, produces (18.79) with γ=κ/σ. As always, we must impose suitable boundary conditions: Dir ichlet conditions u=h that specify the boundary temperature; (homogeneous) Neum ann conditions ∂u/∂n= 0 corresponding to an insulated boundary; or a mixture of the two. Given the initial temperature of the body u(t0,x,y,z) =f(x,y,z) (18 .80) at the initial time t0, it can be proved, [ 48], that the resulting initial-boundary value problem is well-posed, and so there is a unique solution u(t,x,y,z) that is defined for all subsequent times t≥t0and depends continuously on the initial data. As in the one- and two-dimensional versions, we do not lose ge nerality by restricting our attention to homogeneous boundary conditions. Separat ion of variables method works as usual, and we quickly review the basic ideas. One begins by imposing an exponential ansatzu(t,x) =e−λtv(x). Substituting into the differential equation and cancelin g the exponentials, it follows that vsatisfies the Helmholtz eigenvalue problem γ∆v+λv= 0, subject to the relevant boundary conditions. For Dirichlet and mixedboundary conditions, the Laplacian is a positive definite operator, and hence the e igenvalues are all strictly positive, 0<λ1≤λ2≤ ···,withλn−→ ∞,asn→ ∞. Moreover, on a bounded domain, the Helmholtz eigenfunction s are complete, and so linear superposition implies that the solution can be written as an eigenfunction series u(t,x) =∞/summationdisplay n=1cne−λntvn(x). (18.81) The coefficients cnare uniquely prescribed by the initial condition (18.80): u(t0,x) =∞/summationdisplay n=1cne−λnt0vn(x) =f(x). (18.82) 12/11/12 997 c/ci∇cleco√y∇t2012 Peter J. Olver Self-adjointness of the boundary value problem implies ort hogonalityof the eigenfunctions, and hence the coefficients are given by the usual orthogonalit y formulae cn=eλnt0/an}b∇acketle{tf;vn/an}b∇acket∇i}ht /ba∇dblvn/ba∇dbl2=e−λnt0/integraldisplay/integraldisplay/integraldisplay Ωf(x)vn(x)dxdydz /integraldisplay/integraldisplay/integraldisplay Ωvn(x)2dxdydz. (18.83) The resulting solution u(t,x)→0 decays exponentially fast to thermal equilibrium, at a rate equal to the smallest positive eigenvalue λ1>0. Since the higher modes — the terms with n≫0 — go to zero extremely rapidly with increasing t, the solution can be well approximated by the first few terms in its Fourier expans ion. As a consequence, the heat equation rapidly smoothes out discontinuities and eli minates high frequency noise in the initial data, and so can be used to process three-dimen sional images and video — although better nonlinear techniques are now available, [ 162]. Unfortunately, the explicit formulae for the eigenfunctio ns and eigenvalues few and far between, [ 134]. Most explicit eigensolutions of the Helmholtz boundary v alue problem require a further separation of variables. In a rectangular box, one separates into a product of functions depending upon the individual Cartesian coord inates, and the eigenfunctions arewrittenasproductsoftrigonometricandhyperbolicfun ctions; seeExercise fordetails. In a cylindrical domain, the separation is effected in cylind rical coordinates, and leads to eigensolutions involving trigonometric and Bessel functi ons, as outlined in Exercise . The most interesting and enlightening case is a spherical domai n, and we treat this particular problem in complete detail. Heating of a Ball Our goal is to study heat propagation in a solid spherical bod y, e.g., the earth†. For simplicity, we take the diffusivity γ= 1, and consider the heat equation on a solid spherical ballB1={/ba∇dblx/ba∇dbl<1}of radius 1, subject to homogeneous Dirichlet boundary cond itions. Once we know how to solve this particular case, an easy scalin g argument, as outlined in Exercise , will allow us to find the solution for a ball of arbitrary radi us and with a general diffusion coefficient. As usual, when dealing with a spherical geometry, we adopt sp herical coordinates r,ϕ,θas in (18.13), in terms of which the heat equation takes the fo rm ∂u ∂t= ∆u=∂2u ∂r2+2 r∂u ∂r+1 r2∂2u ∂ϕ2+cosϕ r2sinϕ∂u ∂ϕ+1 r2sin2ϕ∂2u ∂θ2, (18.84) where we have used our handy formula (18.14) for the Laplacia n in spherical coordinates. The standard diffusive separation of variables ansatz u(t,r,ϕ,θ) =e−λtv(r,ϕ,θ) †In this simplified model, we are assuming that the earth is composed of a completely homo- geneous and isotropic solid material. 12/11/12 998 c/ci∇cleco√y∇t2012 Peter J. Olver requires us to analyze the spherical coordinate form of the H elmholtz equation ∆v+λv=∂2v ∂r2+2 r∂v ∂r+1 r2∂2v ∂ϕ2+cosϕ r2sinϕ∂v ∂ϕ+1 r2sin2ϕ∂2v ∂θ2+λv= 0 (18.85) on the unit ball Ω = {r<1}with homogeneous Dirichlet boundary conditions. To make further progress, we invoke a second variable separation, s plitting off the radial coordinate by setting v(r,ϕ,θ) =p(r)w(ϕ,θ). The function wmust be 2πperiodic in θand well-defined at the poles ϕ= 0,π. Substi- tuting our ansatz into (18.85), and separating all the r-dependent terms from those terms depending upon the angular variables ϕ,θleads to a pair of differential equations: The first is an ordinary differential equation r2d2p dr2+2rdp dr+(λr2−µ)p= 0, (18.86) for the radial component p(r), while the second is a familiar partial differential equati on ∆Sw+µw=∂2w ∂ϕ2+cosϕ sinϕ∂w ∂ϕ+1 sin2ϕ∂2w ∂θ2+µw= 0, (18.87) for its angular counterpart w(ϕ,θ). The operator ∆Sis thespherical Laplacian from (18.17). In Section 18.2, we showed that its eigenvalues are µm=m(m+1) for m= 0,1,2,3,.... Themtheigenvalue admits 2 m+ 1 linearly independent eigenfunctions — the spherical harmonicsY0 m,...,Ym m,/tildewideY1 m,...,/tildewideYm mdefined in (18.33). The radial ordinary differential equation (18.86) can be sol ved by setting p(r) =√rq(r). We use the product rule to relate their derivatives p=1√rq,dp dr=1√rdq dr−1 2r3/2q,d2p dr2=1√rd2q dr2−1 r3/2dq dr+3 4r5/2q. Substituting these expressions back into (18.86) with µ=µm=m(m+1), and multiplying the resulting equation by√r, we discover that q(r) must solve the differential equation r2d2q dr2+rdq dr+/bracketleftBig λr2−/parenleftbig m+1 2/parenrightbig2/bracketrightBig q= 0, (18.88) which turns out to be the rescaled Bessel equation (17.52) of half integer order m+1 2. As a result, the solution to (18.88) that remains bounded at r= 0 is (up to scalar multiple) the rescaled Bessel function q(r) =Jm+1/2/parenleftbig√ λr./parenrightbig The corresponding solution p(r) =r−1/2Jm+1/2/parenleftbig√ λr/parenrightbig (18.89) to (18.86) is important enough to warrant a special name. 12/11/12 999 c/ci∇cleco√y∇t2012 Peter J. Olver Definition 18.11. Thespherical Bessel function of orderm≥0 is defined by the formula Sm(x) =/radicalbiggπ 2xJm+1/2(x). (18.90) Remark: The multiplicative factor/radicalbig π/2 is included in the definition so as to avoid annoying factors of√πand√ 2 in subsequent formulae. Surprisingly, unlikethe Bessel functions of integer order , the spherical Bessel functions are elementary functions! According to formula (C.67), the spherical Bessel function of order 0 is S0(x) =sinx x. (18.91) Thehigher orderspherical Besselfunctionscanbeobtained byuseofthegeneralrecurrence relation Sm+1(x) =−dSm dx+m xSm(x), (18.92) which is a consequence of Proposition C.13. The next few are, therefore, S1(x) =−dS0 dx=−cosx x+sinx x2, S2(x) =−dS1 dx+S1 x=−sinx x−3cosx x2+3sinx x3, S3(x) =−dS2 dx+2S1 x=cosx x−6sinx x2−15cosx x3+15sinx x4,(18.93) and so on. Graphs can be found in Figure 18.6. Our radial solut ion (18.89) is, apart from an inessential constant multiple, a rescaled spherical Bes sel function of order m: vm(r) =Sm/parenleftbig√ λr/parenrightbig . So far, we have not taken into account the (homogeneous) Diri chlet boundary condi- tion atr= 1. This requires p(1) = 0,and hence Sm/parenleftbig√ λ/parenrightbig = 0. Therefore,√ λmust be a root of the mthorder spherical Bessel function. We introduce the notation 0<σm,1<σm,2<σm,e<··· to denote the successive (positive) spherical Bessel roots , satisfying Sm(σm,n) = 0 for n= 1,2, ... . (18.94) In particular the roots of the zerothorder spherical Bessel function S0(x) =x−1sinxare just the integer multiples of π: σ0,n=nπ forn= 1,2, ... . 12/11/12 1000 c/ci∇cleco√y∇t2012 Peter J. Olver 2 4 6 8 10 12 14 -0.20.20.40.60.81 S0(x)2 4 6 8 10 12 14 -0.20.20.40.60.81 S1(x) 2 4 6 8 10 12 14 -0.20.20.40.60.81 S2(x)2 4 6 8 10 12 14 -0.20.20.40.60.81 S3(x) Figure 18.6. Spherical Bessel Function. Spherical Bessel Roots σm,n n/backslashbigg m0 1 2 3 4 5 6 7 1 3.1416 4.4934 5.7635 8.1826 9.3558 10.5128 11.6570 12.7908... 2 6.2832 7.7253 9.0950 11.7049 12.9665......... 3 9.4248 10.9041 12.3229...... 4 12.5664...... ...... A table of all spherical Bessel roots that are <13 appears above. The columns of the table are indexed by m, the order, while the rows are indexed by n, the root number. Re-assembling the individual constituents, we have now dem onstrated that the sep- arable eigenfunctions of the Helmholtz equation on a solid b all of radius 1, when subject to homogeneous Dirichlet boundary conditions, are product s of spherical Bessel functions and spherical harmonics, vk,m,n(r,ϕ,θ) =Sm(σm,nr)Yk m(ϕ,θ),/tildewidevk,m,n(r,ϕ,θ) =Sm(σm,nr)/tildewideYk m(ϕ,θ).(18.95) 12/11/12 1001 c/ci∇cleco√y∇t2012 Peter J. Olver The corresponding eigenvalues λm,n=σ2 m,n, m = 0,1,2,..., n = 1,2,3,... , (18.96) are given by the squared spherical Bessel roots. Since there are 2m+ 1 independent spherical harmonics of order m, the eigenvalue λm,nadmits 2m+1 linearly independent eigenfunctions, namely v0,m,n,...,vm,m,n,/tildewidev1,m,n,...,/tildewidevm,m,n. In particular, the radially symmetric solutions are the eigenfunctions with k=m= 0, namely vn(r) =v0,0,n(r) =S0(σ0,nr) =sinnπr nπr, n = 1,2,... . (18.97) Further analysis demonstrates that the separable solution s (18.95) form a complete system of eigenfunctions for the Helmholtz equation on the unit bal l subject ot homogeneous Dirichlet boundary conditions, cf. [ 47]. Wehavethuscompletely determinedthebasicseparablesolu tionstotheheatequation on a solid unit ball subject to homogeneous Dirichlet bounda ry conditions. They are products of exponential functions of time, spherical Besse l functions of the radius and spherical harmonics: uk,m,n(t,r,ϕ,θ) =e−σ2 m,ntSm(σm,nr)Yk m(ϕ,θ), /tildewideuk,m,n(t,r,ϕ,θ) =e−σ2 m,ntSm(σm,nr)/tildewideYk m(ϕ,θ).(18.98) The general solution can be written as an infinite “Fourier–B essel–spherical harmonic” series in these fundamental modes u(t,r,θ,ϕ) =∞/summationdisplay m=0∞/summationdisplay n=1e−σ2 m,ntSm(σm,nr)/parenleftBigg c0,m,n 2Y0 m(ϕ,θ) + +m/summationdisplay k=1/bracketleftBig ck,m,nYk m(ϕ,θ)+/tildewideck,m,n/tildewideYk m(ϕ,θ)/bracketrightBig/parenrightBigg .(18.99) The series’ coefficients ck,m,n,/tildewideck,m,nare uniquely prescribed by the initial data; explicit formulae follow from the usual orthogonality relations amo ng the eigenfunctions. Detailed formulae are relegated to the exercises. In particular, the slowest decaying mode is the spherically symmetric function u0,0,1(t,r) =e−π2tsinπr r(18.100) corresponding to the smallest eigenvalue λ0,1=σ2 0,1=π2. Therefore, typically, the decay to thermal equilibrium of a unit sphere is at an exponential r ate ofπ2≈9.8696, or, to a very rough approximation, 10. The Fundamental Solution to the Heat Equation For the heat equation (as well as more general diffusion equat ions), the fundamental solution measures the response of the body to a concentrated unit heat source. Thus, given a pointξ= (ξ,η,ζ)∈Ω within the body, the fundamental solution u(t,x) =F(t,x;ξ) =F(t,x,y,z;ξ,η,ζ) 12/11/12 1002 c/ci∇cleco√y∇t2012 Peter J. Olver solves the initial-boundary value problem ut= ∆u, u (0,x) =δ(x−ξ),forx∈Ω, t>0, (18.101) subject to the selected homogeneous boundary conditions — D irichlet, Neumann or mixed. In general, the fundamental solution has no explicit formul a, although in certain domains it is possible to construct it as an eigenfunction se ries. The one case amenable to a complete analysis is when the heat is distributed over al l of three-dimensional space, so Ω =R3. We recall that Lemma 17.1 showed how to construct solutions of the two- dimensional heat equation as products of one-dimensional s olutions. In a similar manner, ifv(t,x),w(t,x) andq(t,x) are any three solutions to the one-dimensional heat equati on ut=γuxx, then their product u(t,x,y,z) =p(t,x)q(t,y)r(t,z) (18 .102) is a solution to the three-dimensional heat equation ut=γ(uxx+uyy+uzz). In particular, choosing p(t,x) =e−(x−ξ)2/4γt 2√πγt, q (t,y) =e−(y−η)2/4γt 2√πγt, r (t,z) =e−(z−ζ)2/4γt 2√πγt, toallbeone-dimensionalfundamentalsolutions, weareimm ediatelyledtothefundamental solution in the form of a three-dimensional Gaussian kernel . Theorem 18.12. The fundamental solution F(t,x;ξ) =F(t,x−ξ) =e−/bardblx−ξ/bardbl2/4γt 8(πγt)3/2(18.103) solves the three-dimensional heat equation ut=γ∆uonR3with an initial temperature equal to a delta function concentrated at the point x=ξ. Thus, the initially concentrated heat energy immediately b egins to spread out in a radially symmetric manner, with a minuscule, but nonzero eff ect felt at arbitrarily large distances away from the initial concentration. At each indi vidual point x∈R3, after an initial warm-up, the temperature decays back to zero at a rat e proportional to t−3/2— even more rapidly than in two dimensions because, intuitive ly, there are more directions for the heat energy to disperse. Tosolveamoregeneralinitialvalueproblemwiththeinitia ltemperature u(0,x,y,z) = f(x,y,z) distributed over all of space, we first write f(x,y,z) =/integraldisplay/integraldisplay/integraldisplay f(ξ)δ(x−ξ)dξdηdζ as a linear superposition of delta functions. By linearity, the solution to the initial value problem is given by the corresponding superposition u(t,x) =1 8(πγt)3/2/integraldisplay/integraldisplay/integraldisplay f(ξ)e−/bardblx−ξ/bardbl2/4γtdξdηdζ. (18.104) 12/11/12 1003 c/ci∇cleco√y∇t2012 Peter J. Olver of the fundamental solutions. Since the fundamental soluti on has exponential decay as /ba∇dblx/ba∇dbl → ∞, the superposition formula is valid even for initial temper ature distributions which are moderately increasing at large distances. We rema rk that the integral (18.104) has the form of a three-dimensional convolution u(t,x) =F(t,x)∗f(x) =/integraldisplay/integraldisplay/integraldisplay f(ξ)F(t,x−ξ)dξdηdζ (18.105) of the initial data with a one-parameter family of increasin gly spread out Gaussian filters. Consequently, convolution with the Gaussian kernel has a sm oothing effect on the initial temperature distribution. 18.5. The Wave Equation in Three-Dimensional Media. The three-dimensional wave equation utt=c2∆u=c2(uxx+uyy+uzz), (18.106) in whichcdenotes the velocity of light, governs the propagation of el ectromagnetic waves (light, radio, X-rays, etc.) in a homogeneous medium, inclu ding (in the absence of gravita- tional effects) empty space. While the electric and magnetic vector fields E,Bare intrin- sically coupled by the more complicated system of Maxwell’s equations, each individual component satisfies the wave equation; see Exercise for details. The wave equation also models certain classes†of vibrations of a uniform solid body. The solution u(t,x) =u(t,x,y,z) represents a scalar-valued displacement of the body at timetand position x= (x,y,z)∈Ω⊂R3. For example, u(t,x) might represent the radial displacement of the body. One imposes suitable bound ary conditions, e.g., Dirichlet, Neumann or mixed, on ∂Ω, along with a pair of initial conditions u(0,x) =f(x),∂u ∂t(0,x) =g(x), x∈Ω, (18.107) that specify the body’s initial displacement and initial ve locity. As long as the initial and boundary data are reasonably nice, there exists a unique sol ution to the initial-boundary value problem for all −∞< t<∞, cf. [47]. Thus, in contrast to the heat equation, one can follow solutions to the wave equation both forwards and b ackwards in time; see also Exercise . Let us fix our attnetion on the homogeneous boundary value pro blem. The funda- mental vibrational modes are found by imposing our usual tri gonometric ansatz u(t,x,y,z) = cos(ωt)v(x,y,z). †Since the solution u(t,x) to the wave equation is scalar-valued, it cannot measure the full range of possible three-dimensional motions of a solid body. The more comp licated dynamical systems governing the elastic motions of solids are discussed in Exe rcise. 12/11/12 1004 c/ci∇cleco√y∇t2012 Peter J. Olver Substituting into the wave equation (18.106), we discover ( yet again) that v(x,y,z) must be an eigenfunction solving the associated Helmholtz eigen value problem ∆v+λv= 0,where λ=ω2 c2, (18.108) coupled to the relevant boundary conditions. In the positiv e definite cases, i.e., Dirichlet and mixed boundary conditions, the eigenvalues λk=ω2 k/c2>0 are all positive; each eigenfunction vk(x,y,z) yields two normal vibrational modes uk(t,x,y,z) = cos(ωkt)vk(x,y,z),/tildewideuk(t,x,y,z) = sin(ωkt)vk(x,y,z), offrequency ωk=c√ λkequaltothesquarerootofthecorresponding eigenvaluemul tiplied by the wave speed. The general solution is a quasi-periodic l inear combination u(t,x,y,z) =∞/summationdisplay k=1/parenleftbig akcosωkt+bksinωkt/parenrightbig vk(x,y,z) (18.109) of these fundamental vibrational modes. The coefficients ak,bkare uniquely prescribed by the initial conditions (18.107). Thus, u(0,x,y,z) =∞/summationdisplay k=1akvk(x,y,z) =f(x,y,z), ∂u ∂t(0,x,y,z) =∞/summationdisplay k=1ωkbkvk(x,y,z) =g(x,y,z). The explicit formulas follow immediately from the orthogon ality of the eigenfunctions: ak=/an}b∇acketle{tf;vk/an}b∇acket∇i}ht /ba∇dblvk/ba∇dbl2=/integraldisplay/integraldisplay/integraldisplay Ωfvkdxdydz /integraldisplay/integraldisplay/integraldisplay Ωv2 kdxdydz, bk=1 ωk/an}b∇acketle{tg;vk/an}b∇acket∇i}ht /ba∇dblvk/ba∇dbl2=/integraldisplay/integraldisplay/integraldisplay Ωgvkdxdydz ωk/integraldisplay/integraldisplay/integraldisplay Ωv2 kdxdydz. (18.110) In the positive semi-definite Neumann boundary value proble m, there is an additional zero eigenvalue λ0= 0 corresponding to the constant null eigenfunction v0(x,y,z)≡1. This results in two additional terms in the eigenfunction ex pansion — a constant term a0=1 volΩ/integraldisplay/integraldisplay/integraldisplay Ωf(x,y,z)dxdydz that equals the average initial displacement, and an unstab le modeb0tthat grows linearly in time, whose speed b0=1 volΩ/integraldisplay/integraldisplay/integraldisplay Ωg(x,y,z)dxdydz isthe average ofthe initialvelocityover theentire body. T he unstable mode willbe excited if and only if there is a non-zero net initial velocity: b0/ne}ationslash= 0. 12/11/12 1005 c/ci∇cleco√y∇t2012 Peter J. Olver Most of the basic solution techniques we learned in the two-d imensional case apply here, and we will not dwell on the details. The case of a rectan gular box is a particularly straightforward application of the method of separation of variables, and is outlined in the exercises. A similar analysis, now in cylindrical coordina tes, can be applied to the case of a vibrating cylinder. The most interesting case is that of a s olid spherical ball, which is the subject of the next subsection. Vibrations of a Ball Letusfocusontheradialvibrationsofasolidball,asmodel edbythethree-dimensional wave equation (18.106). The solution u(t,x,y,z) represents the radial displacement of the “atom” that is situated at position ( x,y,z) when the ball is at rest. For simplicity, we look at the Dirichlet boundary value prob lem on the unit ball B1={/ba∇dblx/ba∇dbl<1}. The normal modes of vibration are governed by the Helmholtz equation (18.108) subject to homogeneous Dirichlet boundary condit ions. According to (18.95), the eigenfunctions are vk,m,n(r,ϕ,θ) =Sn(σn,mr)Yk m(ϕ,θ), /tildewidevk,m,n(r,ϕ,θ) =Sm(σm,nr)/tildewideYk m(ϕ,θ),forn= 1,2,3,... , m= 0,1,2,... , k= 0,1,...,m.(18.111) HereSmdenotes the mthorder spherical Bessel function (18.90), σm,nis itsnthroot, while Ym n,/tildewideYm nare the spherical harmonics (18.33). Each eigenvalue λm,n=σ2 m,n, m = 0,1,2,..., n = 1,2,3,..., corresponds to 2 m+1 independent eigenfunctions, namely vk,m,0(r,ϕ,θ), vk,m,1(r,ϕ,θ), ... vk,m,m(r,ϕ,θ),/tildewidevk,m,1(r,ϕ,θ), .../tildewidevk,m,m(r,ϕ,θ). Consequently, the fundamental vibrational frequencies of a solid ball ωm,n=c/radicalbig λm,n=cσm,n, m = 0,1,2,..., n = 1,2,3,..., (18.112) are equal to the spherical Bessel roots σm,nmultiplied by the wave speed. There are a total of 2(2 m+1) independent vibrational modes associated with each dis tinct frequency (18.112), namely uk,m,n(t,r,ϕ,θ) = cos(cσm,nt)Sm(σm,nr)Yk m(ϕ,θ), /hatwideuk,m,n(t,r,ϕ,θ) = sin(cσm,nt)Sm(σm,nr)Yk m(ϕ,θ), /tildewideuk,m,n(t,r,ϕ,θ) = cos(cσm,nt)Sm(σm,nr)/tildewideYk m(ϕ,θ), /hatwide/tildewideuk,m,n(t,r,ϕ,θ) = sin(cσm,nt)Sm(σm,nr)/tildewideYk m(ϕ,θ).n= 1,2,3,... , m= 0,1,2,... , k= 0,1,...,m.(18.113) 12/11/12 1006 c/ci∇cleco√y∇t2012 Peter J. Olver Relative Spherical Bessel Roots σk,m/σ0,1 n/backslashbigg m0 1 2 3 4 6 7 8 ... 1 1.0000 1.4303 1.8346 2.2243 2.6046 2.9780 3.3463 3.7105... 2 2.0000 2.4590 2.8950 3.3159 3.7258......... 3 3.0000 3.4709 3.9225...... 4 4.0000...... ...... In particular, the radially symmetric modes of vibration ha ve, according to (18.91), the elementary form u0,0,n(r,ϕ,θ) = cos(cnπt)S0(nπr) =coscnπtsinnπr r, /hatwideu0,0,n(r,ϕ,θ) = sin(cnπt)S0(nπr) =sincnπtsinnπr r,k= 1,2,3,... .(18.114) Their vibrational frequencies, ω0,n=cnπ, are integral multiples of the lowest frequency ω0,1=π. Therefore, interestingly, if you only excite the radially symmetric modes, the resulting motion of the ball is periodic. More generally, adopting the same scaling argument as in (17 .111), we conclude that the fundamental frequencies for a solid ball of radius Rand wave speed care given by ωm,n=cσm,n/R. The relative vibrational frequencies ωm,n ω0,1=σm,n σ0,1=σm,n π(18.115) are independent of the size of the ball Ror the wave speed c. In the accompanying table, we display all relative vibrational frequencies that are le ss than 4 in magnitude. The purely radial modes of vibration (18.114) have individu al frequencies ω0,n=nπc R,soω0,n ω0,1=n, and appear in the first column of the table. The lowest frequen cy isω0,1=πc/R, cor- responding to a vibration with period 2 π/ω0,1= 2R/c. In particular, for the earth, the radiusR≈6,000 km and the wave speed in rock is, on average, c≈5 km/sec, so that the fundamental mode of vibration has period 2 R/c≈2400 seconds, or 40 minutes. Vibra- tions of the earth are also known as seismic waves and, of course, earthquakes are their most severe manifestation. Understanding the modes of vibr ation is an issue of critical im- portance in geophysics and civil engineering, including th e design of structures, buildings and bridges and the avoidance of resonant frequencies. 12/11/12 1007 c/ci∇cleco√y∇t2012 Peter J. Olver Of course, we have suppressed almost all interesting terres trial geology in this very crude approximation, which has been based on the assumption that the earth is a uniform body, vibrating only in its radial direction. A more realist ic modeling of the vibrations of the earth requires an understanding of the basic partial d ifferential equations of linear and nonlinear elasticity, [ 87]. Nonuniformities in the earth lead to scattering of the vi- brational waves, which are then used to locate subterranean geological structures, e.g., oil and gas deposits. We refer the interested reader to [ 6] for a comprehensive introduction to mathematical seismology. Example 18.13. The radial vibrations of a hollow spherical shell (e.g., an e lastic balloon) are governed by the differential equation utt=c2∆S[u] =c2/parenleftbigg∂2u ∂ϕ2+cosϕ sinϕ∂u ∂ϕ+1 sin2ϕ∂2u ∂θ2/parenrightbigg , (18.116) where ∆Sdenotes the spherical Laplacian (18.17). The radial displa cementu(t,ϕ,θ) of a point on the sphere only depends on time tand the angular coordinates ϕ,θ. The solution u(t,ϕ,θ) is required to be 2 πperiodic in the azimuthal angle θand bounded at the poles ϕ= 0,π. According to (18.33), the ntheigenvalueλn=n(n+ 1) of the spherical Laplacian possesses 2 n+1 linearly independent eigenfunctions, namely, the spher ical harmonics Y0 n(ϕ,θ), Y1 n(ϕ,θ), ..., Yn n(ϕ,θ),/tildewideY1 n(ϕ,θ), ...,/tildewideYn n(ϕ,θ). As a consequence, the fundamental frequencies of vibration for a spherical shell are ωn=c/radicalbig λn=c/radicalbig n(n+1), n = 1,2,... . (18.117) The vibrational solutions are quasi-periodic combination s of the fundamental spherical harmonic modes cos/parenleftbig/radicalbig n(n+1)t/parenrightbig Ym n(ϕ,θ), sin/parenleftbig/radicalbig n(n+1)t/parenrightbig Ym n(ϕ,θ), cos/parenleftbig/radicalbig n(n+1)t/parenrightbig/tildewideYm n(ϕ,θ), sin/parenleftbig/radicalbig n(n+1)t/parenrightbig/tildewideYm n(ϕ,θ).(18.118) Representative graphs can be seen in Figure 18.3. The smalle st positive eigenvalue is λ1= 2, yielding a lowest tone of frequency ω1=c√ 2. The higher order frequencies are irrational multiples of the fundamental frequency, implyi ng that a vibrating spherical bell sounds percussive to our ears. One further remark isin order. The spherical Laplacian oper ator is only positivesemi- definite, since the lowest mode has eigenvalue λ0= 0, which corresponds to the constant null eigenfunction v0(ϕ,θ) =Y0 0(ϕ,θ)≡1. Therefore, the wave equation (18.116) admits an unstable mode b0,0t, corresponding to a uniform radial inflation; its coefficient b0,0=3 4π/integraldisplay/integraldisplay S1∂u ∂t(0,ϕ,θ)dS represents the sphere’s average initial velocity. The exis tence of such an unstable mode is an artifact of the simplified linear model we are using, that f ails to account for nonlinearly elastic effects that serve to constrain the inflation of a sphe rical balloon. 12/11/12 1008 c/ci∇cleco√y∇t2012 Peter J. Olver 18.6. Spherical Waves and Huygens’ Principle. For any dynamical (time-varying) partial differential equa tion, the fundamental solu- tion measures the effect of applying an instantaneous concen trated unit impulse at a single point. Two representative physical effects to keep in mind ar e the light waves emanat- ing from a sudden concentrated blast, e.g., a lightning bolt or a stellar supernova, and the sound waves due to an explosionor thunderclap, propagating in air at a much slower speed. Linear superposition utilizes the fundamental solution to build up more general solutions to initial value problems. For the wave and other second orde r vibrational equations, the impulse can be applied either to the initial displacement or to the initial velocity, resulting in two types of fundamental solution. In a uniform isotropic medium, an initial concentrated blas t results in a spherically expanding wave, moving away at the speed of light (or sound) i n all directions. Invoking translation invariance, we will assume that the source of th e disturbance is at the origin, andsothesolution u(t,x)shouldonlydependonthedistance r=/ba∇dblx/ba∇dblfromthesource. We adopt spherical coordinates and look for a solution u=u(t,r) to the three-dimensional wave equation with no angular dependence. Substituting the formula (18.14) for the spherical Laplacian and setting both angular derivatives t o 0, we are led to the following partial differential equation ∂2u ∂t2=c2/parenleftbigg∂2u ∂r2+2 r∂u ∂r/parenrightbigg , (18.119) that governs the propagation of spherically symmetric wave s in three-dimensional space. Surprisingly, we can explicitly solve this partial differen tial equation. The secret is to multiply both sides of the equation by r: ∂2(ru) ∂t2=r∂2u ∂t2=c2/parenleftbigg r∂2u ∂r2+2∂u ∂r/parenrightbigg =c2∂2 ∂r2(ru), Thus, the function w(t,r) =ru(t,r) solves the one-dimensional wave equation ∂2w ∂t2=c2∂2w ∂r2. (18.120) AccordingtoTheorem14.9,thegeneralsolutionto(18.120) canbewrittenind’Alembert form w(t,r) =p(r−ct)+q(r+ct), wherep(ξ) andq(η) are arbitrary functions of a single characteristic variab le. Therefore, spherically symmetric solutions to the three-dimensional wave equation assume the form u(t,r) =p(r−ct) r+q(r+ct) r. (18.121) The first term u(t,r) =p(r−ct) r(18.122) 12/11/12 1009 c/ci∇cleco√y∇t2012 Peter J. Olver represents a wave moving at speed cin the direction of increasing r, and so describes the effect of a variable light source that is concentrated at the o rigin, e.g., a pulsating quasar in interstellar space. To highlight this interpretation, l et us concentrate on the case when p(ξ) =δ(ξ−a) is a delta function, keeping in mind that more general solut ions can then be assembled by linear superposition. The induced solution u(t,r) =δ(r−ct−a) r=δ/parenleftbig r−c(t−t0)/parenrightbig r,where t0=−a c.(18.123) represents a spherical wave propagating through space. At t he instantt=t0, the light is entirely concentrated at the origin r= 0. The signal then moves away from the origin in all directions at speed c. At each later time t>t0, the wave is concentrated on the surface of a sphere of radius r=c(t−t0). Its intensity at each point on the sphere, however, has decreased by a factor 1 /r, and so, the farther from the source, the dimmer the light. A stationary observer sitting at a fixed point in space will onl y see an instantaneous flash of light of intensity 1 /ras the spherical wave passes by at time t=t0+r/c, whereris the observer’s distance from the light source. A similar statem ent holds for sound waves — the sound of the explosion will only last momentarily. Thund er and lightning are the most familiar examples of this everyday phenomenon. On the other hand, for t < t0, the impulse is concentrated at a negative radius r=c(t−t0)<0. To interpret this, note that, for a given value of the spher ical angles ϕ,θ, the point x=rsinϕcosθ, y =rsinϕsinθ, z =rcosϕ, forr <0 lies on the antipodal point of the sphere of radius |r|, so that replacing rby −rhas the same effect as changing xto−x. Thus, the solution (18.123) represents a concentrated spherically symmetric light wave arriving fr om the edges of the universe at speedc, that strengthens inintensity as it collapsesinto the orig inatt=t0. After collapse, it immediately reappears and expands back out into the unive rse. The second solution in the d’Alembert formula (18.121) has, in fact, exactly the same physical form. Indeed, if we set /tildewider=−r,/tildewidep(ξ) =−q(−ξ),thenq(r+ct) r=/tildewidep(/tildewider−ct) /hatwider. Thus, the second d’Alembert solution is redundant, and we on ly need to consider solutions of the form (18.122) from now on. To effectively utilize such spherical wave solutions, we nee d to understand the nature of their originating singularity. For simplicity, we set a= 0 in (18.123) and concentrate on the particular solution u(t,r) =δ(r−ct) r, (18.124) which has a singularity at the origin r= 0 whent= 0. We need to pin down precisely which sort of distribution this solutionrepresents. Invok ing the limiting definition is tricky, 12/11/12 1010 c/ci∇cleco√y∇t2012 Peter J. Olver and it will be easier to work with the dual characterization o f a distribution as a linear functional. Thus, at a fixed time t≥0, we must evaluate the inner product /an}b∇acketle{tu;f/an}b∇acket∇i}ht=/integraldisplay/integraldisplay/integraldisplay u(t,x,y,z)f(x,y,z)dxdydz. of the solution with a smooth test function f(x) =f(x,y,z). We rewrite the triple integral in spherical coordinates using the change of variables form ula (B.68), whereby /an}b∇acketle{tu;f/an}b∇acket∇i}ht=/integraldisplay2π 0/integraldisplayπ 0/integraldisplay∞ 0δ(r−ct) rf(r,ϕ,θ)r2sinϕ drdϕdθ Whent/ne}ationslash= 0, therintegration can be immediately computed, and so /an}b∇acketle{tu;f/an}b∇acket∇i}ht=ct/integraldisplay2π 0/integraldisplayπ 0f(ct,ϕ,θ) sinϕ dϕdθ= 4πctM0 ct[f], (18.125) where, according to (B.43), M0 ct[f] =1 4π/integraldisplay2π 0/integraldisplayπ 0f(ct,ϕ,θ) sinϕ dϕdθ=1 4πc2t2/integraldisplay/integraldisplay Sctf dS (18.126) is the mean or average value of the function fon the sphere Sct=/braceleftbig /ba∇dblx/ba∇dbl=ct/bracerightbig of radius r=ct. In particular, the mean over the limiting sphere of radius r= 0 reduces to the value of the function at the origin: M0 0[f] =f(0). (18.127) Thus, in the limit as t→0, (18.125) implies that /an}b∇acketle{tu;f/an}b∇acket∇i}ht= 0 for allfunctions f, and henceu(0,r)≡0 represents a zero initial displacement. In the absense of any intial diplacement, how, then, can the s olution (18.124) be non- zero? Clearly, this must be the result of a nonzero initial ve locity. To find ut(0,r), we differentiate (18.125) with respect to t, whereby /angbracketleftbigg∂u ∂t;f/angbracketrightbigg =∂ ∂t/an}b∇acketle{tu;f/an}b∇acket∇i}ht=∂ ∂t/parenleftbigg ct/integraldisplay2π 0/integraldisplayπ 0f(ct,ϕ,θ) sinϕ dϕdθ/parenrightbigg =c/integraldisplay2π 0/integraldisplayπ 0f(ct,ϕ,θ) sinϕ dϕdθ+c2t/integraldisplay2π 0/integraldisplayπ 0∂f ∂r(ct,ϕ,θ) sinϕ dϕdθ = 4πcM0 ct[f]+4πc2tM0 ct/bracketleftbigg∂f ∂r/bracketrightbigg . (18.128) The result is a linear combination of the mean of fand of its radial derivative frover the sphere of radius ct. In particular, lim t→0/an}b∇acketle{tut;f/an}b∇acket∇i}ht= 4πcM0 0[f] = 4πcf(0). 12/11/12 1011 c/ci∇cleco√y∇t2012 Peter J. Olver Since this holds for all test functions, we conclude that the initial velocity ut(0,r) = 4πcδ(x) is a multiple of a delta function at the origin! Dividing thro ugh by 4πc, we conclude that the spherical expanding wave u(t,r) =δ(r−ct) 4πcr(18.129) solves the initial value problem u(0,x)≡0,∂u ∂t(0,x) =δ(x), corresponding to an initial unit velocity impulse concentr ated at the origin. This solution can be viewed as the three-dimensional version of the hammer -blow solution (14.128) to the one-dimensional wave equation. A significant difference is that, in three dimensions, there is no residual effect after the wave passes by. More generally, if the unit impulse is concentrated at the po intξ, we invoke transla- tional symmetry to conclude that the function G(t,x;ξ) =δ/parenleftbig /ba∇dblx−ξ/ba∇dbl−ct/parenrightbig 4πc/ba∇dblx−ξ/ba∇dbl, t ≥0, (18.130) isthefundamental solution tothewaveequationresultingfromaconcentratedunitvelo city at the initial time t= 0: G(0,x;ξ) = 0,∂G ∂t(0,x;ξ) =δ(x−ξ). (18.131) We can then apply linear superposition to solve the initial v alue problem u(0,x,y,z) = 0,∂u ∂t(0,x,y,z) =g(x,y,z), (18.132) with zero initial displacement. Namely, we write the initia l velocity g(x) =/integraldisplay/integraldisplay/integraldisplay g(ξ)δ(x−ξ)dξdηdζ as a superposition of impulses, and immediately conclude th at the relevant solution is the self-same superposition of spherical waves: u(t,x) =1 4πc/integraldisplay/integraldisplay/integraldisplay g(ξ)δ/parenleftbig /ba∇dblx−ξ/ba∇dbl−ct/parenrightbig /ba∇dblx−ξ/ba∇dbldξdηdζ=1 4πc2t/integraldisplay/integraldisplay /bardblξ−x/bardbl=ctg(ξ)dS. (18.133) Its value u(t,x) =tMx ct[g], (18.134) at a point xand timet≥0, isttimes the average of the initial velocity function gon a sphere of radius r=ctcentered at the point x. 12/11/12 1012 c/ci∇cleco√y∇t2012 Peter J. Olver 1t rαxSx t B1 Figure 18.7. A Sphere Intersecting a Ball. Example 18.14. Let us set the wave speed c= 1 for simplicity. Suppose that the initial velocity g(x) =/braceleftbigg1,/ba∇dblx/ba∇dbl<1, 0,/ba∇dblx/ba∇dbl>1 is 1 inside the unit ball B1centered at the origin, and 0 outside. To solve the initial va lue problem, we must compute the average value of gover a sphere Sx tof radiust>0 centered at a point x∈R3. Sinceg= 0 outside the unit ball, its average will be equal to the surf ace area of that part of the sphere that is contained inside the un it ball,Sx t∩B1, divided by the total surface area of Sx t, namely 4πt2. To compute this quantity, let r=/ba∇dblx/ba∇dbl. Ift>r+1 or 0<t<r−1, then the sphere of radiustlies entirely outside the unit ball, and so the average is 0; i f 0<t<1−r, then the sphere lies entirely within the unit ball and so the avera ge is 1. Otherwise, referring to Figure 18.7, and referring to Exercise B.5.3, we see that t he area of the spherical cap Sx t∩B1is, by the Law of Cosines, 2πt2(1−cosα) = 2πt2/parenleftbigg 1−r2+t2−1 2rt/parenrightbigg =πt r[1−(t−r)2], whereαdenotes the angle between the line joining the centers of the two spheres and the circle formed by their intersection. Assembling the differe nt subcases, we conclude that Mx ct[g] =  1, 0≤t≤1−r, 1−(t−r)2 4rt,|r−1| ≤t≤r+1, 0, 0≤t≤r−1 ort≥r+1.(18.135) The solution (18.134) is obtained by multiplying by t, and hence for t≥0, u(t,x) =  t, 0≤t≤1−/ba∇dblx/ba∇dbl, 1−/parenleftbig t−/ba∇dblx/ba∇dbl/parenrightbig2 4/ba∇dblx/ba∇dbl,/vextendsingle/vextendsingle/ba∇dblx/ba∇dbl−1/vextendsingle/vextendsingle≤t≤ /ba∇dblx/ba∇dbl+1, 0, 0≤t≤ /ba∇dblx/ba∇dbl−1 ort≥ /ba∇dblx/ba∇dbl+1.(18.136) 12/11/12 1013 c/ci∇cleco√y∇t2012 Peter J. Olver 0.5 11.5 22.5 30.20.40.60.811.2 r= 00.5 11.5 22.5 30.20.40.60.8 r=.3 0.5 11.5 22.5 30.20.40.60.8 r=.70.5 11.5 22.5 30.20.40.60.8 r= 1.3 Figure 18.8. Solution to the Wave Equation Solution due to an Initial Concentrated Velocity. Figure 18.8 plots the solution as a function of time for sever al fixed values of r=/ba∇dblx/ba∇dbl. An observer sitting at the origin will see a linearly increas ing light intensity followed by a sudden decrease to 0. At other points inside the sphere, the d ecrease follows a parabolic arc; if the observer is closer to the edge than the center, the parabolic portion will continue toincreaseforawhilebeforeeventuallytaperingoff. Onthe otherhand, anobserver sitting outside the sphere will experience, after an initially dark period, a symmetrical, parabolic increase to a maximal intensity and then decrease back to dar k after a total time lapse of 2. We also show a plot of uas a function of rfor various fixed times in Figure 18.9. Note that, up until time t= 1, the light spreads out while increasing in intensity near the origin, after which the solution is of gradually decreasing magnitude, supported within the domain lying between two concentric spheres of respective r adiit−1 andt+1. The solution described by formula (18.133)only handles ini tialvelocities. What about solutions reuslting from a nonzero initial displacement? S urprisingly, the answer is differ- entiation! The key observation is that if u(t,x) is any (sufficiently smooth) solution to the wave equation, so is its time derivative v(t,x) =∂u ∂t(t,x). This follows at once from differentiating both sides of the wa ve equation with respect to t andusingtheequalityofmixedpartialderivatives. Physic ally,thisimpliesthatthevelocity of a wave obeys the same evolutionary principle as the wave it self, which is a manifestation of the linearity and time-independence (autonomy) of the eq uation. Suppose uhas initial conditions u(0,x) =f(x), ut(0,x) =g(x). What are the initial conditions for its derivative v=ut? Clearly, its initial displacement v(0,x) =ut(0,x) =g(x) 12/11/12 1014 c/ci∇cleco√y∇t2012 Peter J. Olver 0.5 11.5 22.5 30.20.40.60.811.2 t= 00.5 11.5 22.5 30.20.40.60.811.2 t=.50.5 11.5 22.5 30.20.40.60.811.2 t=.9 0.5 11.5 22.5 30.20.40.60.811.2 t= 1.00.5 11.5 22.5 30.20.40.60.811.2 t= 1.10.5 11.5 22.5 30.20.40.60.811.2 t= 1.5 Figure 18.9. Solution to the Wave Equation Solution due to an Initial Concentrated Velocity. equals the initial velocity of u. As for its initial velocity, we have ∂v ∂t=∂2u ∂t2=c2∆u because we are assuming that usolves the wave equation. Thus, at the initial time ∂v ∂t(0,x) =c2∆u(0,x) =c2∆f(x) equalsc2times the Laplacian of the initial displacement†. In particular, if usatisfies the initial conditions u(0,x) = 0, ut(0,x) =g(x), (18.137) thenv=utsatisfies the initial conditions v(0,x) =g(x), vt(0,x) = 0. (18.138) Thus, paradoxically, to solve the initial displacement pro blem we differentiate the initial velocity solution (18.133) with respect to t, and hence v(t,x) =∂u ∂t(t,x) =∂ ∂t/parenleftbig tMx ct[g]/parenrightbig = Mx ct[g]+ctMx ct/bracketleftbigg∂g ∂n/bracketrightbigg , (18.139) using our computation in (18.128). Therefore, v(t,x) is a linear combination of the mean of the function gand the mean of its normal or radial derivative ∂g/∂n=∂g/∂r, taken †In Section 14.6, a similar device was used to initiate the numerical sol utions to the wave equation. 12/11/12 1015 c/ci∇cleco√y∇t2012 Peter J. Olver 0.5 11.5 22.5 30.20.40.60.811.2 r= 00.5 11.5 22.5 3 -2-1.5-1-0.50.51 r=.3 0.5 11.5 22.5 3 -1-0.50.51 r=.70.5 11.5 22.5 3 -1-0.50.51 r= 1.3 Figure 18.10. Solution to the Wave Equation due to an Initial Concentrated Displacement. over a sphere of radius ctcentered at the point x. In particular, to obtain the solution corresponding to a concentrated initial displacement, F(0,x;ξ) =δ(x−ξ),∂F ∂t(0,x;ξ) = 0, (18.140) we differentiate the solution (18.130), resulting in F(t,x;ξ) =∂G ∂t(t,x;ξ) =−δ′/parenleftbig /ba∇dblx−ξ/ba∇dbl−ct/parenrightbig 4π/ba∇dblξ−x/ba∇dbl, (18.141) which represents a spherically expanding doublet, cf. Figu re 11.10. Thus, interestingly, a concentrated initial displacement spawns an expanding sph erical doublet wave, whereas a concentrated initial velocity spawns a spherical singlet o r delta wave. Example 18.15. Letc= 1 for simplicity. Consider the initial displacement u(0,x) =f(x) =/braceleftbigg1,/ba∇dblx/ba∇dbl<1, 0,/ba∇dblx/ba∇dbl>1 along with zero initial velocity, modeling the effect of an in stantaneously illuminated solid ball. To obtain the solution, we differentiate (18.136) with respect tot, leading to u(t,x) =  1, 0≤t<1−/ba∇dblx/ba∇dbl, /ba∇dblx/ba∇dbl−t 2/ba∇dblx/ba∇dbl,/vextendsingle/vextendsingle/ba∇dblx/ba∇dbl−1/vextendsingle/vextendsingle≤t≤ /ba∇dblx/ba∇dbl+1, 0, 0≤t</ba∇dblx/ba∇dbl−1 ort>1+/ba∇dblx/ba∇dbl.(18.142) As illustrated in Figure 18.10, an observer sitting at the ce nter of the ball will see a constant light intensity until t= 1, at which time the solution suddenly goes dark. At 12/11/12 1016 c/ci∇cleco√y∇t2012 Peter J. Olver 0.5 11.5 22.5 3 -3-2-11 t= 00.5 11.5 22.5 3 -3-2-11 t=.50.5 11.5 22.5 3 -3-2-11 t=.9 0.5 11.5 22.5 3 -3-2-11 t= 1.00.5 11.5 22.5 3 -3-2-11 t= 1.30.5 11.5 22.5 3 -3-2-11 t= 1.8 Figure 18.11. Solution to the Wave Equation due to an Initial Concentrated Displacement. other points inside the ball, 0 < r <1, the downwards jump in intensity arrives sooner, and even goes below 0, followed by a further linear decrease, and finally a jump back to quiescent. An observer placed outside the ball will exper ience, after an initially dark period, a sudden increase in the light intensity, followed b y a linear decrease to negative, followed by a jump back up to darkness. The farther away from t he source, the fainter the light. In Figure 18.11 we plot the same solution as a funct ion ofrfor different values oft. Note the sudden appearnace of a 1 /rsingularity at the origin at time t= 1, due to the focussing of the initial discontinuities in uover the entire unit sphere. Afterwards, the residual disturbance moves off to ∞while gradually decreasing in intensity. Linearly combining the two solutions (18.134,139) establi shesKirchhoff’s formula — although it was first discovered by Poisson — which is the thre e-dimensional counterpart to the d’Alembert’s solution formula for the wave equation. Theorem 18.16. The solution to the initial value problem utt=c2∆u, u (0,x) =f(x),∂u ∂t(0,x) =g(x),x∈R3,(18.143) for the wave equation in three-dimensional space is given by u(t,x) =∂ ∂t/parenleftbig tMx ct[f]/parenrightbig +tMx ct[g] = Mx ct[f]+ctMx ct/bracketleftbigg∂f ∂n/bracketrightbigg +tMx ct[g].(18.144) Here,Mx ct[f]denotes the average of the function fover a sphere of radius ctcentered at position x. A crucially important consequence of the Kirchhoff solution formula is the cele- bratedHuygens’ Principle , which was first highlighted the pioneering seventeenth cen - tury Dutch scientist Christiaan Huygens. Roughly, Huygens ’ Principle states that, in 12/11/12 1017 c/ci∇cleco√y∇t2012 Peter J. Olver three-dimensional space, localized solutions to. the wave equation remain localized. More concretely, (18.144) implis that the value of the solution a t a point xand timetonly depends upon the values of the initial displacements and vel ocities at a distance ctaway. Thus, all signals propagate along the light cone c2t2=x2+y2+z2 in four-dimensional Minkowski space-time. For electromag netic waves, this fact lies at the foundationofspecialrelativity. Physically,Huygens’ Pr inciplemeansthatthelightthatwe see at a given time tarrived from points at a distance exactly d=ctaway at time t= 0. In particular, a sharp, localized initial signal — whether i nitial displacement or initial velocity — that is concentrated near a point produces a sharp , localized response that remains concentrated on an ever expanding sphere surroundi ng the point. In our three- dimensional universe, we only witness the light from an expl osion for a brief moment, after which if there is no subsequent light source, the view return s to darkness. Similarly, a sharp sound remains sharply concentrated, with diminishin g magnitude, as it propagates through space. Remarkably, as we will show next, Huygens’ Pr inciple does not hold in a two dimensional universe! In the plane, concentrated impul ses will be spread out as time progresses. Descent to Two Dimensions So far, we have explicitly determined the solution to the wav e equation in one- and three-dimensional space. The two-dimensional case utt=c2∆u=c2(uxx+uyy). (18.145) is, counter-intuitively, more complicated! For instance, seeking a radially symmetric solu- tionu(t,r) requires solving the partial differential equation ∂2u ∂t2=c2/parenleftbigg∂2u ∂r2+1 r∂u ∂r/parenrightbigg (18.146) which, unlike its three-dimensional cousin (18.119), is no t so easily integrated. However, our solution to the three-dimensional problem can be easily adapted to construct a solution using the so-called Method of Descent . Any solution u(t,x,y) to the two-dimensional wave equation (18.145) can be viewed as a solution to the three- dimensional wave equation (18.106) that does not depend upo n the vertical zcoordinate, whence∂u/∂z= 0. Clearly, if the three-dimensional initial data does not depend on z, then the resulting solution u(t,x,y) will also be independent of z. Consider first the zero initial displacement initial condit ions u(0,x,y) = 0,∂u ∂t(0,x,y) =g(x,y). (18.147) We rewrite the solution formula (18.133) in the form of a surf ace integral over the sphere Sct=/braceleftbig /ba∇dblξ/ba∇dbl=ct/bracerightbig centered at the origin: u(t,x) =1 4πc2t/integraldisplay/integraldisplay Sctg(ξ)dS=1 4πc2t/integraldisplay/integraldisplay /bardblξ/bardbl=ctg(x+ξ)dS. (18.148) 12/11/12 1018 c/ci∇cleco√y∇t2012 Peter J. Olver Imposing the condition that g(x,y) does not depend upon the zcoordinate, we see that the integrals over the upper and lower hemispheres S+ ct=/braceleftbig /ba∇dblξ/ba∇dbl=ct, ζ≥0/bracerightbig , S− ct=/braceleftbig /ba∇dblξ/ba∇dbl=ct, ζ≤0/bracerightbig , are identical. As in (B.47), to evaluate the upper hemispher ical integral, we parametrize the upper hemisphere as the graph of ζ=/radicalbig c2t2−ξ2−η2over the disk Dct=/braceleftbig ξ2+η2≤c2t2/bracerightbig , We conclude that u(t,x,y) =1 2πc2t/integraldisplay/integraldisplay S+ ctg(x+ξ)dS=1 2πc/integraldisplay/integraldisplay Dctg(x+ξ,y+η)/radicalbig c2t2−ξ2−η2dξdη(18.149) solves the initial value problem (18.147). In particular, i f we take the initial velocity g(x,y)=δ(x−ξ)δ(y−η) to be a concentrated impulse, then the resulting solution is G(t,x,y;ξ,η)=  1 2πc/radicalbig c2t2−(x−ξ)2−(y−η)2,(x−ξ)2+(y−η)2<ct, 0, (x−ξ)2+(y−η)2>c2t2. (18.150) An observer placed at position xwill first experience a concentrated displacement singu- larity at time t=/ba∇dblx−ξ/ba∇dbl/c. However, in contrast to the three-dimensional solution, e ven after the impulse passes by, the observer will continue to ex perience a decreasing, but non- zero signal of magnitude roughly proportional to 1 /t. In Figure 18.12, we plot the solution corresponding to a concentrated impulse at the origin, with unit wave speed c= 1. The first line shows the displacement at three different times as a function of r=/ba∇dblx/ba∇dbl; note the initial singularity, indicated by a spike in the graph, is fo llowed by a progressively smaller residual displacement. The second line plots intensity as a function of tat three different radii; the further away from the initial impulse, the faster the residual displacement decays back to 0 — although it never entirely disappears. Similarly, the solution to the initial displacement condit ions u(0,x,y) =f(x,y),∂u ∂t(0,x,y) = 0, (18.151) can be obtained by differentiation with respect to t, and so u(t,x,y) =∂ ∂t/parenleftBigg 1 2πc/integraldisplay/integraldisplay Dctf(x+ξ,y+η)/radicalbig c2t2−ξ2−η2dξdη/parenrightBigg . (18.152) Again, for a concentrated impulse in the initial displaceme nt, an observer will witness, after a certain time lapse, an abrupt impulse passing by that is followed by a progressively decaying residual effect. The general solution to the two-di mensional wave equation on all ofR2is a linear combination of these two types of solutions (18.1 49,152). 12/11/12 1019 c/ci∇cleco√y∇t2012 Peter J. Olver 0.5 11.5 22.5 30.20.40.60.811.21.4 t=.50.5 11.5 22.5 30.20.40.60.811.21.4 t= 10.5 11.5 22.5 30.20.40.60.811.21.4 t= 2 0.5 11.5 22.5 30.20.40.60.811.21.4 r=.50.5 11.5 22.5 30.20.40.60.811.21.4 r= 10.5 11.5 22.5 30.20.40.60.811.21.4 r= 1.5 Figure 18.12. Solution to the Two-Dimensional Wave Equation for a Concentrated Impulse. Thus, Huygens’ Principle is notvalid in a two-dimensional universe. The solution to the two-dimensional wave equation at a point xat timetdepends upon the initial displacement and velocity on the entire disk of radius rtcentered at the point, and not just on the points a distance ctaway. So a two-dimensional creature would experience not only a initial effect of any sound or light wave but also an “ afterglow” with slowly diminishing magnitude. It would be like living in a permanen t echo chamber, and so understanding and acting upon sensory phenomena would more challenging. In general, Huygens’ principle isonlyvalidinodd-dimensional spaces ; see also[ 16] for recent advances in the classification of partial differential equations that admit a Huygens’ principle. Remark: Since the solutions to the two-dimensional wave equation c an be interpreted as three-dimensional solutions with no zdependence, a concentrated delta impulse in the two-dimensional waveequationwouldcorrespondtoa concen trated impulsealonganentire verical line in three dimensions. If light starts propagati ng from the line at t= 0, after the initial signal reaches us, we will continue to receive light from points that lie progressively farther away along the line, and this accounts for the two-di mensional afterglow. 18.7. The Schr¨ odinger Equation and the Hydrogen Atom. Ahydrogen atom consists of a single electron, of mass mand charge e, circling an atomicnucleus containing a single protoninthree-dimensi onal space. Asa result ofquanti- zation, the Schr¨ odinger equation governing the dynamical behavior of the electron around 12/11/12 1020 c/ci∇cleco√y∇t2012 Peter J. Olver the nucleus takes the explicit form i/planckover2pi1∂ψ ∂t=−/planckover2pi12 2m∆ψ−α rψ=−/planckover2pi12 2m/parenleftbigg∂2ψ ∂x2+∂2ψ ∂y2+∂2ψ ∂y2/parenrightbigg −α/radicalbig x2+y2+z2ψ.(18.153) Hereψ(t,x,y,z) denotes the electron’s time-dependent wave function givi ng its quantum probability density as it moves around the nucleus, which, o wing to it relatively small size, is assumed to be concentrated at the origin. The coefficient of the Laplacian depends on Planck’sconstant /planckover2pi1andtheelectronmass m. Thefinaltermrepresents theelectromagnetic (or Newtonian) potential function V=α/rattracting the electron to the nucleus, where α=e2is the square of the electron’s (and proton’s) charge. More g enerally, if the nucleus containsZprotons, then one needs to multiply the potential according ly:α=Ze2. Incidentally, the Schr¨ odinger equation for multi-electr on atoms or even molecules is not hard to write down, but its solution, even for, say, the heliu m atom, is muchmore difficult, and is still a major challenge for numerical approximations on today’s supercomputers. Thus, we will only consider a single electron atom in this sec tion. AccordingtotheanalysisinSection14.7,thenormalmodeso lutionstotheSchr¨ odinger equation are of the form ψ(t,x,y,z) =eiλt//planckover2pi1u(x,y,z), whereuis an eigenfunction of the Hamiltonian operator with eigenv alueλ, and hence /planckover2pi12 2m∆u+/parenleftBig λ+α r/parenrightBig u= 0. (18.154) Thebound states of the atom, in which the electron remains trapped by the nucl eus, are represented the non-zero solutions to the eigenvalue probl em with bounded L2norm: /ba∇dblu/ba∇dbl2=/integraldisplay/integraldisplay/integraldisplay |v(x,y,z)|dxdydz< ∞. The eigenvalue λspecifies the state’s energy, which is necessarily negative :λ<0. Unlike the Laplace equation, the bound states do notform a complete system of eigenfunctions, andsonotevery wavefunction ϕ∈L2(R3)canbeapproximatedbyaneigenfunctionseries. The missing data are the so-called scattering states arising from the continuous spectrum of the Schr¨ odinger operator; these represent electrons th at scatter off of the nucleus, and so do not remain bounded or trapped. We will leave the discuss ion of the scattering states and continuous spectrum to a more advanced treatment, [ 130,154]. Tounderstandtheboundstates,webeginbyrewritingtheeig envalueproblem(18.154) in spherical coordinates: /planckover2pi12 2m/parenleftbigg∂2u ∂r2+2 r∂u ∂r+1 r2∂2u ∂ϕ2+cosϕ r2sinϕ∂u ∂ϕ+1 r2sin2ϕ∂2u ∂θ2/parenrightbigg +/parenleftBig λ+α r/parenrightBig u= 0. (18.155) We then separate off the radial coordinate, setting u(r,ϕ,θ) =v(r)w(ϕ,θ). 12/11/12 1021 c/ci∇cleco√y∇t2012 Peter J. Olver The angular component satisfies the spherical Helmholtz equ ation ∆Sw+µw=∂2w ∂ϕ2+cosϕ sinϕ∂w ∂ϕ+1 sin2ϕ∂2w ∂θ2+µw= 0, thatwehavealreadysolved. Theeigensolutionsarespheric alharmonicswhich, because the quantum mechanical solutions are intrinsically complex-v alued, we take in their complex form (18.41). The associated eigenvalue µ=l(l+1),where the integer l= 0,1,2,... , (18.156) known as the angular quantum number , admits a total of 2 l+ 1 linearly independent eigenfunctions Ym l(θ,ϕ) =Pm l(cosϕ)eimθ, m =−l,−l+1,...,l−1,l. (18.157) The radial equation associagted with the separation consta nt (18.156) is /planckover2pi12 2m/parenleftbiggd2v dr2+2 rdv dr/parenrightbigg +/parenleftbigg λ+α r−l(l+1) r2/parenrightbigg v= 0. (18.158) To eliminate the physical parameters, let’s rescale the rad ial coordinate by setting s=σr, where σ=/radicalbigg −8mλ /planckover2pi12, (18.159) where we use the fact that λ <0. The resulting ordinary differential equation for the rescaled function P(s) =v/parenleftBigs σ/parenrightBig is d2P ds2+2 sdP ds−/parenleftbigg1 4−n s+l(l+1) s2/parenrightbigg P= 0, (18.160) where n=2mα σ/planckover2pi12=2mα σ/planckover2pi12/radicalbigg −m 2λ. (18.161) Since we are searching for bound states, the relevant soluti on(s) of this second order ordi- nary differential equation should be defined on 0 ≤s<∞, remain bounded at s= 0, and go to zero as s→ ∞. The proof of the key result is outlined in the exercises in Ch apter C. Theorem 18.17. The bound state solutions of (18.160) subject to the boundary conditionsP(0+)<∞,lim s→∞P(s) = 0, only occur when n≥l+1is an integer, and are given by P(s) =sle−s/2L2l+1 n−l−1(s), (18.162) where Lk j(s) =k/summationdisplay i=0(−1)i i!/parenleftbiggj+k j−i/parenrightbigg xi, j,k = 0,1,2,... , (18.163) is a certain polynomial function, known as an associated Laguerre polynomial . 12/11/12 1022 c/ci∇cleco√y∇t2012 Peter J. Olver The resulting integer n, whose physical value was noted in (18.161), is known as the principle quantum number . We further note that the scaling factor in (18.159) can be written as σ=2mα n/planckover2pi12=2 nawhere a=/planckover2pi12 mα≈.529×10−10meter is called the Bohr radius , in honor of the pioneering Danish quantum physicist Niels Bohr, and approximates the radius of the first atomic energy l evel. Reverting to physical coordinates, the bound states (18.162) become, up to an ines sential constant multiple, the radial wave functions vn l(r) =/parenleftbigg2r na/parenrightbiggl e−r/(na)L2l+1 n−l−1/parenleftbigg2r na/parenrightbigg . (18.164) Comining them with the spherical harmonics (18.157) yields the atomic eigenfunctions Ulmn(r,ϕ,θ) =vn l(r)Ym l(θ,ϕ). (18.165) These eigenstates depend upon three integers: •l= 0,1,2,3,...: the angular quantum number; •n=l+1,l+2,l+3,...: the principle quantum number; •m=−l,−l+1,...,l−1,l: the magnetic quantum number. The energy of the eigenstate is the associated eigenvalue λn=−α2m 2/planckover2pi121 n2, n = 1,2,3,... . The fact that the ratios λn/λ1= 1/n2between the higher and lowest energy levels are inverse squares of integers was first established by Bohr. Th enthenergy level has a total of n−1/summationdisplay l=0(2l+1) =n2 bound states with that energy, indicating the number of orbi tal shells in the atom associ- ated with that energy level. In chemistry, the electron levels are indexed by the angular quantum number, i.e., the orderlof the spherical harmonic, and traditionally labeled by a le tter in the sequence p,s,d,f,... . Thus, the two order l= 0 spherical harmonics correspond to the pshells; the six harmonics of order l= 1 are the sshells, and so on. Since electrons are allowed to have one of two possible spins, the Pauli exclusion principle tel ls us that each energy shell can be occupied by at most two electrons. Thus, the number of elec trons that can reside in the atomic shell with angular quantum number lis, in fact, 2(2 l+1). The configuration of energy shells and electrons in atoms are responsible for the periodic table. Thus, hydrogen has a single electron in the pshell. Helium has two electrons in the pshell. Lithium has 3 electrons, with two of them filling the first pshell and the third in the second pshell. Neon has 10 electrons filling the two pand first three sshells. And so on. The chemical properties of the elements are, to a very large extent, deter mined by the placement of the electrons within the different shells. See [ Chem] for further details. 12/11/12 1023 c/ci∇cleco√y∇t2012 Peter J. Olver