Please copy and paste this embed script to where you want to embed

A Student’s Guide to Geophysical Equations The advent of accessible student computing packages has meant that geophysics students can now easily manipulate datasets and gain ﬁrst-hand modeling experience – essential in developing an intuitive understanding of the physics of the Earth. Yet to gain a more in-depth understanding of the physical theory, and to be able to develop new models and solutions, it is necessary to be able to derive the relevant equations from ﬁrst principles. This compact, handy book ﬁlls a gap left by most modern geophysics textbooks, which generally do not have space to derive all of the important formulae, showing the intermediate steps. This guide presents full derivations for the classical equations of gravitation, gravity, tides, Earth rotation, heat, geomagnetism, and foundational seismology, illustrated with simple schematic diagrams. It supports students through the successive steps and explains the logical sequence of a derivation – facilitating self-study and helping students to tackle homework exercises and prepare for exams.

william lowrie was born in Hawick, Scotland, and attended the University of Edinburgh, where he graduated in 1960 with ﬁrst-class honors in physics. He achieved a masters degree in geophysics at the University of Toronto and, in 1967, a doctorate at the University of Pittsburgh. After two years in the research laboratory of Gulf Oil Company he became a researcher at the Lamont-Doherty Geological Observatory of Columbia University. In 1974 he was elected professor of geophysics at the ETH Zürich (Swiss Federal Institute of Technology in Zurich), Switzerland, where he taught and researched until retirement in 2004. His research in rock magnetism and paleomagnetism consisted of deducing the Earth’s magnetic ﬁeld in the geological past from the magnetizations of dated rocks. The results were applied to the solution of geologic-tectonic problems, and to analysis of the polarity history of the geomagnetic ﬁeld. Professor Lowrie has authored 135 scientiﬁc articles and a second edition of his acclaimed 1997 textbook Fundamentals of Geophysics was published in 2007. He has been President of the European Union of Geosciences (1987–9) and Section President and Council member of the American Geophysical Union (2000–2). He is a Fellow of the American Geophysical Union and a Member of the Academia Europaea.

A Student’s Guide to Geophysical Equations WILLIAM LOWRIE Institute of Geophysics Swiss Federal Institute of Technology Zurich, Switzerland

cambridge university press Cambridge, New York, Melbourne, Madrid, Cape Town, Singapore, São Paulo, Delhi, Tokyo, Mexico City Cambridge University Press The Edinburgh Building, Cambridge CB2 8RU, UK Published in the United States of America by Cambridge University Press, New York www.cambridge.org Information on this title: www.cambridge.org/9781107005846 © William Lowrie 2011 This publication is in copyright. Subject to statutory exception and to the provisions of relevant collective licensing agreements, no reproduction of any part may take place without the written permission of Cambridge University Press. First published 2011 Printed in the United Kingdom at the University Press, Cambridge A catalogue record for this publication is available from the British Library Library of Congress Cataloguing in Publication data Lowrie, William, 1939– A student’s guide to geophysical equations / William Lowrie. p. cm. Includes bibliographical references and index. ISBN 978-1-107-00584-6 (hardback) 1. Geophysics – Mathematics – Handbooks, manuals, etc. 2. Physics – Formulae – Handbooks, manuals, etc. 3. Earth – Handbooks, manuals, etc. I. Title. QC809.M37L69 2011 550.10 51525–dc22 2011007352 ISBN 978 1 107 00584 6 Hardback ISBN 978 0 521 18377 2 Paperback Cambridge University Press has no responsibility for the persistence or accuracy of URLs for external or third-party internet websites referred to in this publication, and does not guarantee that any content on such websites is, or will remain, accurate or appropriate.

This book is dedicated to Marcia

Contents

Preface Acknowledgments 1

2

page xi xiii

Mathematical background 1.1 Cartesian and spherical coordinates 1.2 Complex numbers 1.3 Vector relationships 1.4 Matrices and tensors 1.5 Conservative force, ﬁeld, and potential 1.6 The divergence theorem (Gauss’s theorem) 1.7 The curl theorem (Stokes’ theorem) 1.8 Poisson’s equation 1.9 Laplace’s equation 1.10 Power series 1.11 Leibniz’s rule 1.12 Legendre polynomials 1.13 The Legendre differential equation 1.14 Rodrigues’ formula 1.15 Associated Legendre polynomials 1.16 Spherical harmonic functions 1.17 Fourier series, Fourier integrals, and Fourier transforms Further reading Gravitation 2.1 Gravitational acceleration and potential 2.2 Kepler’s laws of planetary motion 2.3 Gravitational acceleration and the potential of a solid sphere 2.4 Laplace’s equation in spherical polar coordinates 2.5 MacCullagh’s formula for the gravitational potential Further reading

vii

1 1 1 4 8 17 18 20 23 26 28 32 32 34 41 43 49 52 58 59 59 60 66 69 74 85

viii

3

4

5

6

7

8

Contents

Gravity 3.1 The ellipticity of the Earth’s ﬁgure 3.2 The geopotential 3.3 The equipotential surface of gravity 3.4 Gravity on the reference spheroid 3.5 Geocentric and geographic latitude 3.6 The geoid Further reading The tides 4.1 Origin of the lunar tide-raising forces 4.2 Tidal potential of the Moon 4.3 Love’s numbers and the tidal deformation 4.4 Tidal friction and deceleration of terrestrial and lunar rotations Further reading Earth’s rotation 5.1 Motion in a rotating coordinate system 5.2 The Coriolis and Eötvös effects 5.3 Precession and forced nutation of Earth’s rotation axis 5.4 The free, Eulerian nutation of a rigid Earth 5.5 The Chandler wobble Further reading Earth’s heat 6.1 Energy and entropy 6.2 Thermodynamic potentials and Maxwell’s relations 6.3 The melting-temperature gradient in the core 6.4 The adiabatic temperature gradient in the core 6.5 The Grüneisen parameter 6.6 Heat ﬂow Further reading Geomagnetism 7.1 The dipole magnetic ﬁeld and potential 7.2 Potential of the geomagnetic ﬁeld 7.3 The Earth’s dipole magnetic ﬁeld 7.4 Secular variation 7.5 Power spectrum of the internal ﬁeld 7.6 The origin of the internal ﬁeld Further reading Foundations of seismology 8.1 Elastic deformation

86 86 88 91 96 102 106 115 116 116 119 124 130 136 137 138 140 142 155 157 169 170 171 172 176 178 179 182 197 198 198 200 205 213 214 217 225 227 227

Contents

ix

8.2 Stress 8.3 Strain 8.4 Perfectly elastic stress–strain relationships 8.5 The seismic wave equation 8.6 Solutions of the wave equation 8.7 Three-dimensional propagation of plane P- and S-waves Further reading

228 233 239 244 252 254 258

Appendix A Magnetic poles, the dipole ﬁeld, and current loops Appendix B Maxwell’s equations of electromagnetism References Index

259 265 276 278

Preface

This work was written as a supplementary text to help students understand the mathematical steps in deriving important equations in classical geophysics. It is not intended to be a primary textbook, nor is it intended to be an introduction to modern research in any of the topics it covers. It originated in a set of handouts, a kind of “do-it-yourself ” manual, that accompanied a course I taught on theoretical geophysics. The lecture aids were necessary for two reasons. First, my lectures were given in German and there were no comprehensive up-to-date texts in the language; the recommended texts were in English, so the students frequently needed clariﬁcation. Secondly, it was often necessary to explain classical theory in more detail than one ﬁnds in a multi-topic advanced textbook. To keep such a book as succinct as possible, the intermediate steps in the mathematical derivation of a formula must often be omitted. Sometimes the unassisted student cannot ﬁll in the missing steps without individual tutorial assistance, which is usually in short supply at most universities, especially at large institutions. To help my students in these situations, the “do-it-yourself ” text that accompanied my lectures explained missing details in the derivations. This is the background against which I prepared the present guide to geophysical equations, in the hope that it might be helpful to other students at this level of study. The classes that I taught to senior grades were largely related to potential theory and primarily covered topics other than seismology, since this was the domain of my colleagues and better taught by a true seismologist than by a paleomagnetist! Theoretical seismology is a large topic that merits its own treatment at an advanced level, and there are several textbooks of classical and modern vintage that deal with this. However, a short chapter on the relationship of stress, strain, and the propagation of seismic waves is included here as an introduction to the topic. Computer technology is an essential ingredient of progress in modern geophysics, but a well-trained aspiring geophysicist must be able to do more than xi

xii

Preface

apply advanced software packages. A fundamental mathematical understanding is needed in order to formulate a geophysical problem, and numerical computational skills are needed to solve it. The techniques that enabled scientists to understand much about the Earth in the pre-computer era also underlie much of modern methodology. For this reason, a university training in geophysics still requires the student to work through basic theory. This guide is intended as a companion in that process. Historically, most geophysicists came from the ﬁeld of physics, for which geophysics was an applied science. They generally had a sound training in mathematics. The modern geophysics student is more likely to have begun studies in an Earth science discipline, the mathematical background might be heavily oriented to the use of tailor-made packaged software, and some students may be less able to handle advanced mathematical topics without help or tutoring. To ﬁll these needs, the opening chapter of this book provides a summary of the mathematical background for topics handled in subsequent chapters.

Acknowledgments

In writing this book I have beneﬁted from the help and support of various people. At an early stage, anonymous proposal reviewers gave me useful suggestions, not all of which have been acted on, but all of which were appreciated. Each chapter was read and checked by an obliging colleague. I wish to thank Dave Chapman, Rob Coe, Ramon Egli, Chris Finlay, Valentin Gischig, Klaus Holliger, Edi Kissling, Emile Klingelé, Alexei Kuvshinov, Germán Rubino, Rolf Sidler, and Doug Smylie for their corrections and suggestions for improvement. The responsibility for any errors that escaped scrutiny is, of course, mine. I am very grateful to Derrick Hasterok and Dave Chapman for providing me with an unpublished ﬁgure from Derrick’s Ph.D. thesis. Dr. Susan Francis, Senior Commissioning Editor at Cambridge University Press, gave me constant support and friendly encouragement throughout the many months of writing, for which I am sincerely grateful. Above all, I thank my wife Marcia for her generous tolerance of the intrusion of this project into our retirement activities.

xiii

1 Mathematical background

1.1 Cartesian and spherical coordinates Two systems of orthogonal coordinates are used in this book, sometimes interchangeably. Cartesian coordinates (x, y, z) are used for a system with rectangular geometry, and spherical polar coordinates (r, θ, ) are used for spherical geometry. The relationship between these reference systems is shown in Fig. 1.1(a). The convention used here for spherical geometry is deﬁned as follows: the radial distance from the origin of the coordinates is denoted r, the polar angle θ (geographic equivalent: the co-latitude) lies between the radius and the z-axis (geographic equivalent: Earth’s rotation axis), and the azimuthal angle in the x–y plane is measured from the x-axis (geographic equivalent: longitude). Position on the surface of a sphere (constant r) is described by the two angles θ and . The Cartesian and spherical polar coordinates are linked as illustrated in Fig. 1.1(b) by the relationships x ¼ r sin θ cos y ¼ r sin θ sin z ¼ r cos θ

(1:1)

1.2 Complex numbers The numbers we most commonly use in daily life are real numbers. Some of them are also rational numbers. This means that they can be expressed as the quotient of two integers, with the condition that the denominator of the quotient must not equal zero. When the denominator is 1, the real number is an integer. Thus 4, 4/5, 123/456 are all rational numbers. A real number can also be irrational, which means it cannot be expressed as the quotient of two integers. 1

2

Mathematical background

z

(a)

(b)

z = r cosθ

θ r θ

φ

r

φ

y x

r sinθ

y = r sinθ sinφ

x = r sinθ cosφ

Fig. 1.1. (a) Cartesian and spherical polar reference systems. (b) Relationships between the Cartesian and spherical polar coordinates.

Imaginary axis z = x + iy

+y r sinθ

θ

r

r cosθ +x

Real axis

Fig. 1.2. Representation of a complex number on an Argand diagram.

Familiar examples are π, e (the base of natural logarithms), and some square roots, such as √2, √3, √5, etc. The irrational numbers are real numbers that do not terminate or repeat when expressed as decimals. In certain analyses, such as determining the roots of an equation, it is necessary to ﬁnd the square root of a negative real number, e.g. √(–y2), where y is real. The result is an imaginary number. The negative real number can be written as (–1)y2, and its square root is then √(–1)y. The quantity √(–1) is written i and is known as the imaginary unit, so that √(–y2) becomes ±iy. A complex number comprises a real part and an imaginary part. For example, z = x + iy, in which x and y are both real numbers, is a complex number with a real part x and an imaginary part y. The composition of a complex number can be illustrated graphically with the aid of the complex plane (Fig. 1.2). The real part is plotted on the horizontal axis, and the imaginary part on the vertical axis. The two independent parts are orthogonal on the plot and the complex number z

1.2 Complex numbers

3

is represented by their vector sum, deﬁning a point on the plane. The distance r of the point from the origin is given by r¼

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ x2 þ y2

(1:2)

The line joining the point to the origin makes an angle θ with the real (x-)axis, and so r has real and imaginary components r cos θ and r sin θ, respectively. The complex number z can be written in polar form as z ¼ rðcos θ þ i sin θÞ

(1:3)

It is often useful to write a complex number in the exponential form introduced by Leonhard Euler in the late eighteenth century. To illustrate this we make use of inﬁnite power series; this topic is described in Section 1.10. The exponential function, exp(x), of a variable x can be expressed as a power series as in (1.135). On substituting x = iθ, the power series becomes ðiθÞ2 ðiθÞ3 ðiθÞ4 ðiθÞ5 ðiθÞ6 þ þ þ þ þ 2! 3! 4! 5! 6! ðiθÞ2 ðiθÞ4 ðiθÞ6 ðiθÞ3 ðiθÞ5 ¼1þ þ þ þ ðiθÞ þ þ þ 2! 4! 6! 3! 5! θ2 θ4 θ6 θ3 θ 5 (1:4) ¼ 1 þ þ þ i θ þ þ 2! 4! 6! 3! 5!

expðiθÞ ¼ 1 þ ðiθÞ þ

Comparison with (1.135) shows that the ﬁrst bracketed expression on the right is the power series for cos θ; the second is the power series for sin θ. Therefore expðiθÞ ¼ cos θ þ i sin θ

(1:5)

On inserting (1.5) into (1.3), the complex number z can be written in exponential form as z ¼ r expðiθÞ

(1:6)

The quantity r is the modulus of the complex number and θ is its phase. Conversely, using (1.5) the cosine and sine functions can be deﬁned as the sum or difference of the complex exponentials exp(iθ) and exp(–iθ): expðiθÞ þ expðiθÞ 2 expðiθÞ expðiθÞ sin θ ¼ 2i

cos θ ¼

(1:7)

4

Mathematical background

1.3 Vector relationships A scalar quantity is characterized only by its magnitude; a vector has both magnitude and direction; a unit vector has unit magnitude and the direction of the quantity it represents. In this overview the unit vectors for Cartesian coordinates (x, y, z) are written (ex, ey, ez); unit vectors in spherical polar coordinates (r, θ, ) are denoted (er, eθ, e). The unit vector normal to a surface is simply denoted n.

1.3.1 Scalar and vector products The scalar product of two vectors a and b is deﬁned as the product of their magnitudes and the cosine of the angle α between the vectors: a · b ¼ ab cos α

(1:8)

If the vectors are orthogonal, the cosine of the angle α is zero and a·b ¼ 0

(1:9)

The vector product of two vectors is another vector, whose direction is perpendicular to both vectors, such that a right-handed rule is observed. The magnitude of the vector product is the product of the individual vector magnitudes and the sine of the angle α between the vectors: ja bj ¼ ab sin α

(1:10)

If a and b are parallel, the sine of the angle between them is zero and ab¼0

(1:11)

Applying these rules to the unit vectors (ex, ey, ez), which are normal to each other and have unit magnitude, it follows that their scalar products are ex · ey ¼ ey · ez ¼ ez · ex ¼ 0 e x · ex ¼ e y · ey ¼ e z · ez ¼ 1

(1:12)

The vector products of the unit vectors are ex ey ¼ ez ey ez ¼ ex ez ex ¼ ey ex ex ¼ ey ey ¼ ez ez ¼ 0

(1:13)

1.3 Vector relationships

5

A vector a with components (ax, ay, az) is expressed in terms of the unit vectors (ex, ey, ez) as a ¼ ax ex þ ay ey þ az ez

(1:14)

The scalar product of the vectors a and b is found by applying the relationships in (1.12): a · b ¼ ax ex þ ay ey þ az ez · bx ex þ by ey þ bz ez (1:15) ¼ ax bx þ ay by þ az bz The vector product of the vectors a and b is found by using (1.13): a b ¼ ax ex þ ay ey þ az ez bx ex þ by ey þ bz ez ¼ ay bz az by ex þ ðaz bx ax bz Þey þ ax by ay bx ez

(1:16)

This result leads to a convenient way of evaluating the vector product of two vectors, by writing their components as the elements of a determinant, as follows: ex ey ez a b ¼ ax ay az (1:17) bx by bz The following relationships may be established, in a similar manner to the above, for combinations of scalar and vector products of the vectors a, b, and c: a · ðb cÞ ¼ b · ðc aÞ ¼ c · ða bÞ

(1:18)

a ðb cÞ ¼ bðc · aÞ cða · bÞ

(1:19)

ða bÞ c ¼ bðc · aÞ aðb · cÞ

(1:20)

1.3.2 Vector differential operations The vector differential operator ∇ is deﬁned relative to Cartesian axes (x, y, z) as r ¼ ex

∂ ∂ ∂ þ ey þ ez ∂x ∂y ∂z

(1:21)

The vector operator ∇ determines the gradient of a scalar function, which may be understood as the rate of change of the function in the direction of each of the reference axes. For example, the gradient of the scalar function φ with respect to Cartesian axes is the vector

6

Mathematical background

rφ ¼ ex

∂φ ∂φ ∂φ þ ey þ ez ∂x ∂y ∂z

(1:22)

The vector operator ∇ can operate on either a scalar quantity or a vector. The scalar product of ∇ with a vector is called the divergence of the vector. Applied to the vector a it is equal to ∂ ∂ ∂ · ax ex þ ay ey þ az ez r · a ¼ ex þ ey þ ez ∂x ∂y ∂z ∂ax ∂ay ∂az (1:23) ¼ þ þ ∂x ∂y ∂z If the vector a is deﬁned as the gradient of a scalar potential φ, as in (1.22), we can substitute potential gradients for the vector components (ax, ay, az). This gives ∂ ∂φ ∂ ∂φ ∂ ∂φ r · rφ ¼ þ þ (1:24) ∂x ∂x ∂y ∂y ∂z ∂z By convention the scalar product (∇ · ∇) on the left is written ∇2. The resulting identity is very important in potential theory and is encountered frequently. In Cartesian coordinates it is r2 φ ¼

∂2 φ ∂2 φ ∂2 φ þ þ ∂x2 ∂y2 ∂z2

(1:25)

The vector product of ∇ with a vector is called the curl of the vector. The curl of the vector a may be obtained using a determinant similar to (1.17): ex ey ez r a ¼ ∂=∂x ∂=∂y ∂=∂z (1:26) ax ay az In expanded format, this becomes ∂az ∂ay ∂ax ∂az ∂ay ∂ax ex þ ey þ ez ra¼ ∂y ∂z ∂z ∂x ∂x ∂y

(1:27)

The curl is sometimes called the rotation of a vector, because of its physical interpretation (Box 1.1). Some commonly encountered divergence and curl operations on combinations of the scalar quantity φ and the vectors a and b are listed below: r · ðφaÞ ¼ ðrφÞ · a þ φðr · aÞ

(1:28)

1.3 Vector relationships

7

Box 1.1. The curl of a vector The curl of a vector at a given point is related to the circulation of the vector about that point. This interpretation is best illustrated by an example, in which a ﬂuid is rotating about a point with constant angular velocity ω. At distance r from the point the linear velocity of the ﬂuid v is equal to ω × r. Taking the curl of v, and applying the identity (1.31) with ω constant, r v ¼ r ðw rÞ ¼ wðr · rÞ ðw · rÞr

(1)

To evaluate the ﬁrst term on the right, we use rectangular coordinates (x, y, z): ∂ ∂ ∂ wðr · rÞ ¼ w ex þ ey þ ez · xex þ yey þ zez ∂x ∂y ∂z ¼ w ex · ex þ ey · ey þ ez · ez ¼ 3w (2) The second term is

ðw · rÞr ¼

ωx

∂ ∂ ∂ þ ωy þ ωz · xex þ yey þ zez ∂x ∂y ∂z

¼ ωx ex þ ωy ey þ ωz ez ¼ w

(3)

Combining the results gives r v ¼ 2w

(4)

1 w ¼ ð r vÞ 2

(5)

Because of this relationship between the angular velocity and the linear velocity of a ﬂuid, the curl operation is often interpreted as the rotation of the ﬂuid. When ∇ × v = 0 everywhere, there is no rotation. A vector that satisﬁes this condition is said to be irrotational.

r · ða bÞ ¼ b · ðr aÞ a · ðr bÞ

(1:29)

r ðφaÞ ¼ ðrφÞ a þ φðr aÞ

(1:30)

r ða bÞ ¼ aðr · bÞ bðr · aÞ ða · rÞb þ ðb · rÞa r ðrφÞ ¼ 0

(1:31) (1:32)

8

Mathematical background

m3 z0

n3 z

θ3

θ1

θ2 φ3

φ1

x0 m1

χ3

χ1 φ2

χ2

y

n2

y0 m 2

x n1

Fig. 1.3. Two sets of Cartesian coordinate axes, (x, y, z) and (x0, y0, z0), with corresponding unit vectors (n1, n2, n3) and (m1, m2, m3), rotated relative to each other.

r · ðr a Þ ¼ 0

(1:33)

r ðr aÞ ¼ rðr · aÞ r2 a

(1:34)

It is a worthwhile exercise to establish these identities from basic principles, especially (1.19) and (1.31)–(1.34), which will be used in later chapters.

1.4 Matrices and tensors 1.4.1 The rotation matrix Consider two sets of orthogonal Cartesian coordinate axes (x, y, z) and (x0, y0, z0) that are inclined to each other as in Fig. 1.3. The x0-axis makes angles (1, χ1, θ1) with each of the (x, y, z) axes in turn. Similar sets of angles (2, χ2, θ2) and (3, χ3, θ3) are deﬁned by the orientations of the y0- and z0-axes, respectively, to the (x, y, z) axes. Let the unit vectors along the (x, y, z) and (x0, y0, z0) axes be (n1, n2, n3) and (m1, m2, m3), respectively. The vector r can be expressed in either system, i.e., r = r(x, y, z) = r(x0, y0, z0), or, in terms of the unit vectors, r ¼ xn1 þ yn2 þ zn3 ¼ x0 m1 þ y0 m2 þ z0 m3

(1:35)

We can write the scalar product (r · m1) as r · m1 ¼ xn1 · m1 þ yn2 · m1 þ zn3 · m1 ¼ x0

(1:36)

The scalar product (n1 · m1) = cos 1 = α11 deﬁnes α11 as the direction cosine of the x0-axis with respect to the x-axis (Box 1.2). Similarly, (n2 · m1) = cos χ1 = α12

1.4 Matrices and tensors

9

and (n3 · m1) = cos θ1 = α13 deﬁne α12 and α13 as the direction cosines of the x0-axis with respect to the y- and z-axes, respectively. Thus, (1.36) is equivalent to x0 ¼ α11 x þ α12 y þ α13 z

(1:37)

On treating the y0- and z0-axes in the same way, we get their relationships to the (x, y, z) axes: y0 ¼ α21 x þ α22 y þ α23 z z0 ¼ α31 x þ α32 y þ α33 z

(1:38)

The three equations can be written as a single matrix equation 2 3 2 32 3 2 3 x0 α11 α12 α13 x x 4 y0 5 ¼ 4 α21 α22 α23 54 y 5 ¼ M4 y 5 z z z0 α31 α32 α33

(1:39)

The coefﬁcients αnm (n = 1, 2, 3; m = 1, 2, 3) are the cosines of the interaxial angles. By deﬁnition, α12 = α21, α23 = α32, and α31 = α13, so the square matrix M is symmetric. It transforms the components of the vector in the (x, y, z) coordinate system to corresponding values in the (x0, y0, z0) coordinate system. It is thus equivalent to a rotation of the reference axes. Because of the orthogonality of the reference axes, useful relationships exist between the direction cosines, as shown in Box 1.2. For example, ðα11 Þ2 þ ðα12 Þ2 þ ðα13 Þ2 ¼ cos2 1 þ cos2 χ 1 þ cos2 θ1 ¼

1 2 x þ y2 þ z2 ¼ 1 2 r (1:40)

and α11 α21 þ α12 α22 þ α13 α23 ¼ cos 1 cos 2 þ cos χ 1 cos χ 2 þ cos θ1 cos θ2 ¼ 0 (1:41) The last summation is zero because it is the cosine of the right angle between the x0-axis and the y0-axis. These two results can be summarized as 3 X 1; m¼n αmk αnk ¼ (1:42) 0; m 6¼ n k¼1

1.4.2 Eigenvalues and eigenvectors The transpose of a matrix X with elements αnm is a matrix with elements αmn (i.e., the elements in the rows are interchanged with corresponding elements in

10

Mathematical background

Box 1.2. Direction cosines The vector r is inclined at angles α, β, and γ, respectively, to orthogonal reference axes (x, y, z) with corresponding unit vectors (ex, ey, ez), as in Fig. B1.2. The vector r can be written r ¼ xex þ yey þ zez

(1)

where (x, y, z) are the components of r with respect to these axes. The scalar products of r with ex, ey, and ez are r

z

γ

α

x

β

ez ex

ey

y

Fig. B1.2. Angles α, β, and γ deﬁne the tilt of a vector r relative to orthogonal reference axes (x, y, z), respectively. The unit vectors (ex, ey, ez) deﬁne the coordinate system.

r · ex ¼ x ¼ r cos α r · ey ¼ y ¼ r cos β r · ez ¼ z ¼ r cos γ

(2)

Therefore, the vector r in (1) is equivalent to r ¼ ðr cos αÞex þ ðr cos βÞey þ ðr cos γÞez

(3)

The unit vector u in the direction of r has the same direction as r but its magnitude is unity: r u ¼ ¼ ðcos αÞex þ ðcos βÞey þ ðcos γÞez ¼ lex þ mey þ nez r where (l, m, n) are the cosines of the angles that the vector r makes with the reference axes, and are called the direction cosines of r. They are useful for describing the orientations of lines and vectors.

(4)

1.4 Matrices and tensors

11

The scalar product of two unit vectors is the cosine of the angle they form. Let u1 and u2 be unit vectors representing straight lines with direction cosines (l1, m1, n1) and (l2, m2, n2), respectively, and let θ be the angle between the vectors. The scalar product of the vectors is u1 · u2 ¼ cos θ ¼ l1 ex þ m1 ey þ n1 ez · l2 ex þ m2 ey þ n2 ez (5) Therefore, cos θ ¼ l1 l2 þ m1 m2 þ n1 n2

(6)

The square of a unit vector is the scalar product of the vector with itself and is equal to 1: u·u ¼

r·r ¼1 r2

(7)

On writing the unit vector u as in (4), and applying the orthogonality conditions from (2), we ﬁnd that the sum of the squares of the direction cosines of a line is unity: (8) lex þ mey þ nez · lex þ mey þ nez ¼ l2 þ m2 þ n2 ¼ 1

the columns). The transpose of a (3 × 1) column matrix is a (1 × 3) row matrix. For example, if X is a column matrix given by 2 3 x X ¼ 4y5 (1:43) z then its transpose is the row matrix X T, where XT ¼ ½ x T

The matrix equation X MX = K, where surface: 2 α11 XT MX ¼ ½ x y z 4 α21 α31

y

z

(1:44)

K is a constant, deﬁnes a quadric α12 α22 α32

32 3 x α13 α23 54 y 5 ¼ K z α33

The symmetry of the matrix leads to the equation of this surface:

(1:45)

12

Mathematical background

fðx; y; zÞ ¼ α11 x2 þ α22 y2 þ α33 z2 þ 2α12 xy þ 2α23 yz þ 2α31 zx ¼ K (1:46) When the coefﬁcients αnm are all positive real numbers, the geometric expression of the quadratic equation is an ellipsoid. The normal direction n to the surface of the ellipsoid at the point P(x, y, z) is the gradient of the surface. Using the relationships between (x, y, z) and (x0, y0, z0) in (1.39) and the symmetry of the rotation matrix, αnm = αmn for n ≠ m, the normal direction has components ∂f ¼ 2ðα11 x þ α12 y þ α13 zÞ ¼ 2x0 ∂x ∂f ¼ 2ðα21 x þ α22 y þ α23 zÞ ¼ 2y0 ∂y ∂f ¼ 2ðα31 x þ α32 y þ α33 zÞ ¼ 2z0 ∂z

(1:47)

and we write nðx; y; zÞ ¼ rf ¼ ex

∂f ∂f ∂f þ ey þ e z ∂x ∂x ∂x

(1:48)

nðx; y; zÞ ¼ 2ðx0 ex þ y0 ey þ z0 ez Þ ¼ 2rðx0 ; y0 ; z0 Þ (1:49) The normal n to the surface at P(x, y, z) in the original coordinates is parallel to the vector r at the point (x0, y0, z0) in the rotated coordinates (Fig. 1.4). The transformation matrix M has the effect of rotating the reference axes from one orientation to another. A particular matrix exists that will cause the directions of the (x0, y0, z0) axes to coincide with the (x, y, z) axes. In this case the normal to the surface of the ellipsoid is one of the three principal axes of the

tangent plane

z x

(x0,y0,z0) r

er n P(x,y,z)

y

Fig. 1.4. Location of a point (x, y, z) on an ellipsoid, where the normal n to the surface is parallel to the radius vector at the point (x0, y0, z0).

1.4 Matrices and tensors

13

ellipsoid. The component x0 is then proportional to x, y0 is proportional to y, and z0 is proportional to z. Let the proportionality constant be β. Then x0 = βx, y0 = βy, and z0 = βz, and we get the set of simultaneous equations ðα11 βÞx þ α12 y þ α13 z ¼ 0 α21 x þ ðα22 βÞy þ α23 z ¼ 0

(1:50)

α31 x þ α32 y þ ðα33 βÞz ¼ 0 which, in matrix form, is 2

α11 β 4 α21 α31

α12 α22 β α32

32 3 x α13 α23 54 y 5 ¼ 0 α33 β z

(1:51)

The simultaneous equations have a non-trivial solution only if the determinant of coefﬁcients is zero, i.e., α11 β α21 α31

α12 α22 β α32

α13 α23 ¼ 0 α33 β

(1:52)

This equation is a third-order polynomial in β. Its three roots (β1, β2, β3) are known as the eigenvalues of the matrix M. When each eigenvalue βn is inserted in turn into (1.50) it deﬁnes the components of a corresponding vector vn, which is called an eigenvector of M. Note that (1.51) is equivalent to the matrix equation 2

α11 4 α21 α31

α12 α22 α32

32 3 2 α13 x 1 α23 54 y 5 β4 0 α33 z 0

0 1 0

32 3 0 x 0 54 y 5 ¼ 0 1 z

(1:53)

which we can write in symbolic form ðM βIÞX ¼ 0

(1:54)

The matrix I, with diagonal elements equal to 1 and off-diagonal elements 0, is called a unit matrix: 2

1 0 I ¼ 40 1 0 0

3 0 05 1

(1:55)

14

Mathematical background

1.4.3 Tensor notation Equations describing vector relationships can become cumbersome when written in full or symbolic form. Tensor notation provides a succinct alternative way of writing the equations. Instead of the alphabetic indices used in the previous section, tensor notation uses numerical indices that allow summations to be expressed in a compact form. Let the Cartesian coordinates (x, y, z) be replaced by coordinates (x1, x2, x3) and let the corresponding unit vectors be (e1, e2, e3). The vector a in (1.14) becomes X a ¼ a1 e1 þ a2 e2 þ a3 e3 ¼ ai e i (1:56) i¼1;2;3

A convention introduced by Einstein drops the summation sign and tacitly assumes that repetition of an index implies summation over all values of the index, in this case from 1 to 3. The vector a is then written explicitly a ¼ ai e i

(1:57)

Alternatively, the unit vectors can be implied and the expression ai is understood to represent the vector a. Using the summation convention, (1.15) for the scalar product of two vectors a and b is a · b ¼ a1 b1 þ a2 b2 þ a3 b3 ¼ ai bi

(1:58)

Suppose that two vectors a and b are related, so that each component of a is a linear combination of the components of b. The relationship can be expressed in tensor notation as ai ¼ Tij bj

(1:59)

The indices i and j identify components of the vectors a and b; each index takes each of the values 1, 2, and 3 in turn. The quantity Tij is a second-order (or second-rank) tensor, representing the array of nine coefﬁcients (i.e., 32). A vector has three components (i.e., 31) and is a ﬁrst-order tensor; a scalar property has a single (i.e., 30) value, its magnitude, and is a zeroth-order tensor. To write the cross product of two vectors we need to deﬁne a new quantity, the Levi-Civita permutation tensor εijk. It has the value +1 when a permutation of the indices is even (i.e., ε123 = ε231 = ε312 = 1) and the value –1 when a permutation of the indices is odd (i.e., ε132 = ε213 = ε321 = –1). If any pair of indices is equal, εijk = 0. This enables us to write the cross product of two vectors in tensor notation. Let u be the cross product of vectors a and b: u ¼ a b ¼ ða2 b3 a3 b2 Þe1 þ ða3 b1 a1 b3 Þe2 þ ða1 b2 a2 b1 Þe3 (1:60)

1.4 Matrices and tensors

15

In tensor notation this is written ui ¼ εijk aj bk

(1:61)

This can be veriﬁed readily for each component of u. For example, u1 ¼ ε123 a2 b3 þ ε132 a3 b2 ¼ a2 b3 a3 b2

(1:62)

The tensor equivalent to the unit matrix deﬁned in (1.55) is known as Kronecker’s symbol, δij, or alternatively the Kronecker delta. It has the values 1; if i ¼ j δij ¼ (1:63) 0; if i 6¼ j Kronecker’s symbol is convenient for selecting a particular component of a tensor equation. For example, (1.54) can be written in tensor form using the Kronecker symbol: (1:64) αij βδij xj ¼ 0 This represents the set of simultaneous equations in (1.50). Likewise, the relationship between direction cosines in (1.42) simpliﬁes to αmk αnk ¼ δmn

(1:65)

in which a summation over the repeated index is implied.

1.4.4 Rotation of coordinate axes Let vk be a vector related to the coordinates xl by the tensor Tkl vk ¼ Tkl xl

(1:66)

A second set of coordinates x′n is rotated relative to the axes xl so that the direction cosines of the angles between corresponding axes are the elements of the tensor αnl: x0n ¼ αnl xl

(1:67)

Let the same vector be related to the rotated coordinate axes x′n by the tensor T ′kn: v0k ¼ T 0kn x0n

(1:68)

vk and v′k are the same vector, expressed relative to different sets of axes. Therefore, v0k ¼ αkn vn ¼ αkn Tnl xl

(1:69)

16

Mathematical background

Equating the expressions in (1.68) and (1.69) for v′k gives 0 0 T kn xn ¼ αkn Tnl xl

(1:70)

Using the relationships between the axes in (1.67), 0 0 0 xn ¼ T kn αnl xl T kn

(1:71)

0 αnl ¼ αkn Tnl T kn

(1:72)

Therefore,

On multiplying by αml and summing, 0 αml αnl T kn ¼ αml αkn Tnl

(1:73)

Note that in expanded form the products of direction cosines on the left are equal to αml αnl ¼ αm1 αn1 þ αm2 αn2 þ αm3 αn3 ¼ δmn

(1:74)

as a result of (1.42). Therefore the transformation matrix in the rotated coordinate system is related to the original matrix by the direction cosines between the two sets of axes: 0 ¼ αml αkn Tnl T km

(1:75)

The indices m and k can be interchanged without affecting the result. The sequence of terms in the summation changes, but its sum does not. Therefore, 0 ¼ αkl αmn Tnl T km

(1:76)

This relationship allows us to compute the elements of a matrix in a new coordinate system that is rotated relative to the original reference axes by angles that have the set of direction cosines αnl.

1.4.5 Vector differential operations in tensor notation In tensor notation the vector differential operator ∇ in Cartesian coordinates becomes r ¼ ei

∂ ∂xi

(1:77)

The gradient of a scalar function φ with respect to Cartesian unit vectors (e1, e2, e3) is therefore

1.5 Conservative force, ﬁeld, and potential

rφ ¼ e1

∂φ ∂φ ∂φ ∂φ þ e2 þ e3 ¼ ei ∂x1 ∂x2 ∂x3 ∂xi

17

(1:78)

Several shorthand forms of this equation are in common use; for example, ∂φ ¼ ðrφÞi ¼ φ;i ¼ ∂i φ ∂xi

(1:79)

The divergence of the vector a is written in tensor notation as r·a ¼

∂a1 ∂a2 ∂a3 ∂ai þ þ ¼ ¼ ∂ i ai ∂x1 ∂x2 ∂x3 ∂xi

The curl (or rotation) of the vector a becomes ∂a3 ∂a2 ∂a1 ∂a3 ∂a2 ∂a1 þ e2 þ e3 r a ¼ e1 ∂x2 ∂x3 ∂x3 ∂x1 ∂x1 ∂x2 ðr aÞi ¼ εijk

∂ak ¼ εijk ∂j ak ∂xj

(1:80)

(1:81)

(1:82)

1.5 Conservative force, ﬁeld, and potential If the work done in moving an object from one point to another against a force is independent of the path between the points, the force is said to be conservative. No work is done if the end-point of the motion is the same as the starting point; this condition is called the closed-path test of whether a force is conservative. In a real situation, energy may be lost, for example to heat or friction, but in an ideal case the total energy E is constant. The work dW done against the force F is converted into a gain dEP in the potential energy of the displaced object. The change in the total energy dE is zero: dE ¼ dEP þ dW ¼ 0

(1:83)

The change in potential energy when a force with components (Fx, Fy, Fz) parallel to the respective Cartesian coordinate axes (x, y, z) experiences elementary displacements (dx, dy, dz) is dEP ¼ dW ¼ Fx dx þ Fy dy þ Fz dz (1:84) The value of a physical force may vary in the space around its source. For example, gravitational and electrical forces decrease with distance from a source mass or electrical charge, respectively. The region in which a physical

18

Mathematical background

quantity exerts a force is called its ﬁeld. Its geometry is deﬁned by lines tangential to the force at any point in the region. The term ﬁeld is also used to express the value of the force exerted on a unit of the quantity. For example, the electric ﬁeld of a charge is the force experienced by a unit charge at a given point; the gravitational ﬁeld of a mass is the force acting on a unit of mass; it is therefore equivalent to the acceleration. In a gravitational ﬁeld the force F is proportional to the acceleration a. The Cartesian components of F are therefore (max, may, maz). The gravitational potential U is deﬁned as the potential energy of a unit mass in the gravitational ﬁeld, thus dEP = m dU. After substituting these expressions into (1.84) we get dU ¼ ax dx þ ay dy þ az dz (1:85) The total differential dU can be written in terms of partial differentials as dU ¼

∂U ∂U ∂U dx þ dy þ dz ∂x ∂y ∂z

(1:86)

On equating coefﬁcients of dx, dy, and dz in these equations: ax ¼

∂U ; ∂x

ay ¼

∂U ; ∂y

az ¼

∂U ∂z

(1:87)

These relationships show that the acceleration a is the negative gradient of a scalar potential U: a ¼ rU

(1:88)

Similarly, other conservative ﬁelds (e.g., electric, magnetostatic) can be derived as the gradient of the corresponding scalar potential. According to the vector identity (1.32) the curl of a gradient is always zero; it follows from (1.88) that a conservative force-ﬁeld F satisﬁes the condition rF¼0

(1:89)

1.6 The divergence theorem (Gauss’s theorem) Let n be the unit vector normal to a surface element of area dS. The ﬂux dΦ of a vector F across the surface element dS (Fig. 1.5) is deﬁned to be the scalar product dΦ ¼ F · n dS

(1:90)

If the angle between F and n is θ, the ﬂux across dS is dΦ ¼ F dS cos θ

(1:91)

1.6 The divergence theorem

19

F n

θ dS dS n

Fig. 1.5. The ﬂux of a vector F across a small surface dS, whose normal n is inclined to the vector, is equal to the ﬂux across a surface dSn normal to the vector.

z

dz Fx Fx + dFx

y x

dy x + dx x

Fig. 1.6. Figure for computing the change in the ﬂux of a vector in the x-direction for a small box with edges (dx, dy, dz).

where F is the magnitude of F. Thus the ﬂux of F across the oblique surface dS is equivalent to that across the projection dSn (=dS cos θ) of dS normal to F. Consider the net ﬂux of the vector F through a rectangular box with edges dx, dy, and dz parallel to the x-, y-, and z-axes, respectively (Fig. 1.6). The area dSx of a side normal to the x-axis equals dy dz. The x-component of the vector at x, where it enters the box, is Fx, and at x + dx, where it leaves the box, it is Fx + dFx. The net ﬂux in the x-direction is dΦx ¼ ððFx þ dFx Þ Fx ÞdSx ¼ dFx dy dz

(1:92)

If the distance dx is very small, the change in Fx may be written to ﬁrst order as dFx ¼

∂Fx dx ∂x

(1:93)

20

Mathematical background

The net ﬂux in the x-direction is therefore dΦx ¼

∂Fx ∂Fx dx dy dz ¼ dV ∂x ∂x

(1:94)

where dV is the volume of the small element. Similar results are obtained for the net ﬂux in each of the y- and z-directions. The total ﬂux of F through the rectangular box is the sum of these ﬂows: dΦ ¼ dΦx þ dΦy þ dΦz dΦ ¼

∂Fx ∂Fy ∂Fz þ þ dV ¼ ðr · FÞdV ∂x ∂y ∂z

(1:95) (1:96)

We can equate this expression with the ﬂux deﬁned in (1.90). The ﬂux through a ﬁnite volume V with a bounding surface of area S and outward normal unit vector n is ZZZ ZZ ðr · FÞdV ¼ F · n dS (1:97) V

S

This is known as the divergence theorem, or Gauss’s theorem, after the German mathematician Carl Friedrich Gauss (1777–1855). Note that the surface S in Gauss’s theorem is a closed surface, i.e., it encloses the volume V. If the ﬂux of F entering the bounding surface is the same as the ﬂux leaving it, the total ﬂux is zero, and so r·F ¼ 0

(1:98)

This is sometimes called the continuity condition because it implies that ﬂux is neither created nor destroyed (i.e., there are neither sources nor sinks of the vector) within the volume. The vector is said to be solenoidal.

1.7 The curl theorem (Stokes’ theorem) Stokes’ theorem relates the surface integral of the curl of a vector to the circulation of the vector around a closed path bounding the surface. Let the vector F pass through a surface S which is divided into a grid of small elements (Fig. 1.7). The area of a typical surface element is dS and the unit vector n normal to the element speciﬁes its orientation. First, we evaluate the work done by F around one of the small grid elements, ABCD (Fig. 1.8). Along each segment of the path we need to consider only the

1.7 The curl theorem

21

n S C

D dS

B

A

C Fig. 1.7. Conﬁguration for Stokes’ theorem: the surface S is divided into a grid of elementary areas dS and is bounded by a closed circuit C.

y

(x + dx, y + dy)

(x, y + dy) D

C

dy

Fy

F

A

(x, y)

B

dx

(x + dx, y)

x

Fx

Fig. 1.8. Geometry for calculation of the work done by a force F around a small rectangular grid.

vector component parallel to that segment. The value of F may vary with position, so, for example, the x-component along AB may differ from the x-component along CD. Provided that dx and dy are inﬁnitesimally small, we can use Taylor series approximations for the components of F (Section 1.10.2). To ﬁrst order we get ∂Fx dy ∂y ∂Fy dx Fy BC ¼ Fy DA þ ∂x ðFx ÞCD ¼ ðFx ÞAB þ

(1:99)

The work done in a circuit around the small element ABCD is the sum of the work done along each individual segment: xZþdx

I F·dl ¼ ABCD

yþdy Z

Fy BC dy þ

ðFx ÞAB dx þ x

y

Zx

Zy ðFx ÞCD dx þ

xþdx

Fy

DA

dy

yþdy

(1:100)

22

Mathematical background

I

xþdx Z

F·dl ¼

ðFx ÞAB ðFx ÞCD dx þ

x

ABCD

yþdy Z

Fy BC Fy DA dy

y

(1:101) Substituting from (1.99) gives xþdx Z

I

yþdy Z ∂Fx ∂Fy dy dx þ dx dy ∂y ∂x

F·dl ¼ ABCD

x

(1:102)

y

The mean-value theorem allows us to replace the integrands over the tiny distances dx and dy by their values at some point in the range of integration: I ∂Fy ∂Fx dx dy (1:103) F·dl ¼ ∂x ∂y ABCD

The bracketed expression is the z-component of the curl of F I F · d l ¼ ðr FÞz dx dy

(1:104)

ABCD

The normal direction n to the small area dS = dx dy is parallel to the z-axis (i.e., out of the plane of Fig. 1.8), and hence is in the direction of (∇ × F)z. Thus, I F · d l ¼ ðr FÞ·n dS (1:105) ABCD

The circuit ABCD is one of many similar grid elements of the surface S. When adjacent elements are compared, the line integrals along their common boundary are equal and opposite. If the integration is carried out for the entire surface S, the only surviving parts are the integrations along the bounding curve C (Fig. 1.7). Thus ZZ I ðr FÞ · n dS ¼ F · d l (1:106) S

C

This equation is known as Stokes’ theorem, after the English mathematician George Gabriel Stokes (1819–1903). It enables conversion of the surface integral of a vector to a line integral. The integration on the left is made over the surface S through which the vector F passes. The closed integration on the right is made around the bounding curve C to the surface S; d l is an inﬁnitesimal element of this

1.8 Poisson’s equation

23

boundary. The direction of dl around the curve is right-handed with respect to the surface S, i.e., positive when the path is kept to the right of the surface, as in Fig. 1.7. Note that the surface S in Stokes’ theorem is an open surface; it is like the surface of a bowl with the bounding curve C as its rim. The integration of F around the rim is called the circulation of F about the curve C. If the integral is zero, there is no circulation and the vector F is said to be irrotational. Comparison with the left-hand side shows that the condition for this is rF¼0

(1:107)

As shown in Section 1.5, this is also the condition for F to be a conservative ﬁeld.

1.8 Poisson’s equation The derivations in this and the following sections are applicable to any ﬁeld that varies as the inverse square of distance from its source. Gravitational acceleration is used as an example, but the electric ﬁeld of a charge may be treated in the same way. Let S be a surface enclosing an observer at P and a point mass m. Let dS be a small element of the surface at distance r in the direction er from the mass m, as in Fig. 1.9. The orientation of dS is speciﬁed by the direction n normal to the surface element. With G representing the gravitational constant (see Section 2.1), the gravitational acceleration aG at dS is given by aG ¼ G

m er r2

(1:108)

Let θ be the angle between the radius and the direction n normal to the surface element, and let the projection of dS normal to the radius be dSn. The solid angle dΩ with apex at the mass is deﬁned as the ratio of the normal surface element dSn to the square of its distance r from the mass (Box 1.3):

S dS m

dΩ

aG

θ

r

er n

dSn

Fig. 1.9. Representation of the ﬂux of the gravitational acceleration aG through a closed surface S surrounding the source of the ﬂux (the point mass m).

24

Mathematical background

Box 1.3. Deﬁnition of a solid angle A small element of the surface of a sphere subtends a cone with apex at the center of the sphere (Fig. B1.3(a)). The solid angle Ω is deﬁned as the ratio of the area A of the surface element to the square of the radius r of the sphere: (a)

r dθ

(b) dθ

dA r

Ω

r dA

θ

A

r sinθ dφ

dφ

r sinθ

Fig. B1.3. (a) Relationship of the solid angle Ω, the area A of an element subtended on the surface of a sphere, and the radius r of the sphere. (b) The surface of a sphere divided into rings, and each ring into small surface elements with sides r dθ and r sin θ d.

Ω¼

A r2

(1)

This deﬁnition can be used for an arbitrarily shaped surface. If the surface is inclined to the radial direction it must be projected onto a surface normal to the radius, as in Fig. 1.5. For example, if the normal to the surface A makes an angle α with the direction from the apex of the subtended cone, the projected area is A cos α and the solid angle Ω subtended by the area is Ω¼

A cos α r2

(2)

As an example, let the area on the surface of a sphere be enclosed by a small circle (Fig. B1.3(b)). Symmetry requires spherical polar coordinates to describe the area within the circle. Let the circle be divided into concentric rings, and let the half-angle subtended by a ring at the center of the sphere be θ. The radius of the ring is r sin θ and its width is r dθ. Let the angular position of a small surface element of the ring be ; the length of a side of the element is then r sin θ d. The area dA of a surface area element is equal to r2 sin θ dθ d. The solid angle subtended at the center of the sphere by the element of area dA is dΩ ¼

r2 sin θ dθ d ¼ sin θ dθ d r2

(3)

1.8 Poisson’s equation

25

This expression is also equivalent to the element of area on the surface of a unit sphere (one with unit radius). Integrating in the ranges 0 ≤ ≤ 2π and 0 ≤ θ ≤ θ0, we get the solid angle Ω0 subtended by a circular region of the surface of a sphere deﬁned by a half-apex angle θ0: Zθ0 Z2π Ω0 ¼

sin θ dθ d ¼ 2π ð1 cos θ0 Þ

(4)

θ¼0 ¼0

The unit of measurement of solid angle is the steradian, which is analogous to the radian in plane geometry. The maximum value of a solid angle is when the surface area is that of the complete sphere, namely 4πr2. The solid angle at its center then has the maximum possible value of 4π. This result is also obtained by letting the half-apex angle θ0 in (4) increase to its maximum value π.

dΩ ¼

dSn dS cos θ ðer · nÞdS ¼ ¼ r2 r2 r2

(1:109)

The ﬂux dN of the gravitational acceleration aG through the area element is m ðer · nÞdS r2

(1:110)

cos θ dS ¼ Gm dΩ r2

(1:111)

dN ¼ aG · n dS ¼ G

dN ¼ Gm

If we integrate this expression over the entire surface S we get the total gravitational ﬂux N, ZZ Z N¼ aG · n dS ¼ Gm dΩ ¼ 4πGm (1:112) Ω

S

Now we replace this surface integral by a volume integration, using the divergence theorem (Section 1.6) ZZZ ZZ ðr · aG ÞdV ¼ aG · n dS ¼ 4πGm (1:113) V

S

This is valid for any point mass m inside the surface S. If the surface encloses many point masses we may replace m with the sum of the point masses. If mass

26

Mathematical background

is distributed in the volume with mean density ρ, a volume integral can replace the enclosed mass: ZZZ ZZZ ðr · aG ÞdV ¼ 4πG ρ dV (1:114) V

V

ZZZ ðr · aG þ 4πGρÞdV ¼ 0

(1:115)

V

For this to be generally true the integrand must be zero. Consequently, r · aG ¼ 4πGρ

(1:116)

The gravitational acceleration is the gradient of the gravitational potential UG as in (1.88): r · ðrUG Þ ¼ 4πGρ r2 UG ¼ 4πGρ

(1:117) (1:118)

Equation (1.118) is known as Poisson’s equation, after Siméon-Denis Poisson (1781–1840), a French mathematician and physicist. It describes the gravitational potential of a mass distribution at a point that is within the mass distribution. For example, it may be used to compute the gravitational potential at a point inside the Earth.

1.9 Laplace’s equation Another interesting case is the potential at a point outside the mass distribution. Let S be a closed surface outside the point mass m. The radius vector r from the point mass m now intersects the surface S at two points A and B, where it forms angles θ1 and θ2 with the respective unit vectors n1 and n2 normal to the surface (Fig. 1.10). Let er be a unit vector in the radial direction. Note that the outward normal n1 forms an obtuse angle with the radius vector at A. The gravitational acceleration at A is a1 and its ﬂux through the surface area dS1 is Gm dN1 ¼ a1 · n1 dS1 ¼ 2 ðr · n1 ÞdS1 (1:119) r1 Gm cos θ1 dS1 (1:120) dN1 ¼ 2 cosðπ θ1 ÞdS1 ¼ Gm r1 r21

1.9 Laplace’s equation

27

S

B m

A

a1

dΩ

θ2

a2

dS1 n1

dS2

θ1

er n2

Fig. 1.10. Representation of the gravitational ﬂux through a closed surface S that does not enclose the source of the ﬂux (the point mass m).

dN1 ¼ Gm dΩ

(1:121)

The gravitational acceleration at B is a2 and its ﬂux through the surface area dS2 is dN2 ¼ a2 · n2 dS2 ¼ Gm

cos θ2 dS2 r22

dN2 ¼ Gm dΩ

(1:122) (1:123)

The total contribution of both surfaces to the gravitational ﬂux is dN ¼ dN1 þ dN2 ¼ 0

(1:124)

Thus, the total ﬂux of the gravitational acceleration aG through a surface S that does not include the point mass m is zero. By invoking the divergence theorem we have for this situation Z Z ðr · aG ÞdV ¼ aG · n dS ¼ 0 (1:125) V

S

For this result to be valid for any volume, the integrand must be zero: r · aG ¼ r · ðrUG Þ ¼ 0

(1:126)

r2 UG ¼ 0

(1:127)

Equation (1.127) is Laplace’s equation, named after Pierre Simon, Marquis de Laplace (1749–1827), a French mathematician and physicist. It describes the gravitational potential at a point outside a mass distribution. For example, it is applicable to the computation of the gravitational potential of the Earth at an external point or on its surface. In Cartesian coordinates, which are rectilinear, Laplace’s equation has the simple form

28

Mathematical background

∂2 UG ∂2 UG ∂2 UG þ þ ¼0 ∂x2 ∂y2 ∂z2

(1:128)

Spherical polar coordinates are curvilinear and the curvature of the angular coordinates results in a more complicated form: 1 ∂ 2 ∂UG 1 ∂ ∂UG 1 ∂2 UG þ þ ¼ 0 (1:129) r sin θ ∂r ∂θ r2 ∂r r2 sin θ ∂θ r2 sin2 θ ∂2

1.10 Power series A function ƒ(x) that is continuous and has continuous derivatives may be approximated by the sum of an inﬁnite series of powers of x. Many mathematical functions – e.g., sin x, cos x, exp(x), ln(1 + x) – fulﬁll these conditions of continuity and can be expressed as power series. This often facilitates the calculation of a value of the function. Three types of power series will be considered here: the MacLaurin, Taylor, and binomial series.

1.10.1 MacLaurin series Let the function ƒ(x) be written as an inﬁnite sum of powers of x: fðxÞ ¼ a0 þ a1 x þ a2 x2 þ a3 x3 þ a4 x4 þ þ an xn þ

(1:130)

The coefﬁcients an in this sum are constants. Differentiating (1.130) repeatedly with respect to x gives df ¼ a1 þ 2a2 x þ 3a3 x2 þ 4a4 x3 þ þ nan xn1 þ dx d 2f ¼ 2a2 þ ð3 2Þa3 x þ ð4 3Þa4 x2 þ þ nðn 1Þan xn2 þ dx2 d 3f ¼ ð3 2Þa3 þ ð4 3 2Þa4 x þ þ nðn 1Þðn 2Þan xn3 þ dx3 (1:131) After n differentiations, the expression becomes d nf ¼ ðnðn 1Þðn 2Þ . . . 3 2 1Þ an þ terms containing powers of x dxn (1:132) Now we evaluate each of the differentiations at x = 0. Terms containing powers of x are zero and

1.10 Power series df fð0Þ ¼ a0 ; ¼ a1 dx x¼0 2 3 d f d f ¼ 2a ; ¼ ð3 2Þa3 2 dx2 x¼0 dx3 x¼0 n d f ¼ ðnðn 1Þðn 2Þ . . . 3 2 1Þan ¼ n!an dxn x¼0

29

(1:133)

On inserting these values for the coefﬁcients into (1.130) we get the power series for ƒ(x): df x2 d 2 f x3 d 3 f þ þ þ fðxÞ ¼ fð0Þ þ x dx x¼0 2! dx2 x¼0 3! dx3 x¼0 (1:134) xn d n f þ þ n! dxn x¼0 This is the MacLaurin series for ƒ(x) about the origin, x = 0. It was derived in the eighteenth century by the Scottish mathematician Colin MacLaurin (1698– 1746) as a special case of a Taylor series. The MacLaurin series is a convenient way to derive series expressions for several important functions. In particular, x3 x5 x7 x2n1 þ þ ð1Þn1 3! 5! 7! ð2n 1Þ! x2 x 4 x6 x2n2 cos x ¼ 1 þ þ ð1Þn1 2! 4! 6! ð2n 2Þ! (1:135) x2 x3 xn1 x expðxÞ ¼ e ¼ 1 þ x þ þ þ þ þ 2! 3! ðn 1Þ! x2 x3 x4 xn lnð1 þ xÞ ¼ loge ð1 þ xÞ ¼ x þ þ ð1Þn1 2 3 4 n sin x ¼ x

1.10.2 Taylor series We can write the power series in (1.134) for ƒ(x) centered on any new origin, for example x = x0. To do this we substitute (x – x0) for x in the above derivation. The power series becomes df ðx x0 Þ2 d 2 f fðxÞ ¼fðx0 Þ þ ðx x0 Þ þ 2! dx x¼x0 dx2 x¼x0 ðx x0 Þ3 d 3 f ðx x0 Þn d n f þ þ þ þ dx3 x¼x0 dxn x¼x0 3! n! (1:136)

30

Mathematical background

This is called a Taylor series, after an English mathematician, Brooks Taylor (1685–1731), who described its properties in 1712. The MacLaurin and Taylor series are both approximations to the function ƒ(x). The remainder between the true function and its power series is a measure of how well the function is expressed by the series.

1.10.3 Binomial series Finite series An important series is the expansion of the function ƒ(x) = (a + x)n. If n is a positive integer, the expansion of ƒ(x) is a ﬁnite series, terminating after (n + 1) terms. Evaluating the series for some low values of n gives the following: n ¼ 0 : ða þ xÞ0 ¼ 1 n ¼ 1 : ða þ xÞ1 ¼ a þ x n ¼ 2 : ða þ xÞ2 ¼ a2 þ 2ax þ x2

(1:137)

n ¼ 3 : ða þ xÞ3 ¼ a3 þ 3a2 x þ 3ax2 þ x3 n ¼ 4 : ða þ xÞ4 ¼ a4 þ 4a3 x þ 6a2 x2 þ 4ax3 þ x4 The general expansion of ƒ(x) is therefore nðn 1Þ n2 2 a x þ 12 nðn 1Þ . . . ðn k þ 1Þ nk k a x þ xn þ k!

ða þ xÞn ¼ an þ nan1 x þ

(1:138)

The coefﬁcient of the general kth term is equivalent to nðn 1Þ . . . ðn k þ 1Þ n! ¼ k! k!ðn kÞ!

(1:139)

This is called the binomial coefﬁcient. When the constant a is equal to 1 and n is a positive integer, we have the useful series expansion ð1 þ xÞn ¼

n X

n! xk k! ð n k Þ! k¼0

(1:140)

Inﬁnite series If the exponent in (1.140) is not a positive integer, the series does not terminate, but is an inﬁnite series. The series for ƒ(x) = (1 + x)p, in which the exponent p is not a positive integer, may be derived as a MacLaurin series:

1.10 Power series

df dx

¼ pð1 þ xÞp1

x¼0

x¼0

31

¼p

d 2f ¼ pðp 1Þð1 þ xÞp2 ¼ pðp 1Þ 2 x¼0 dx x¼0 n d f ¼ ðpðp1Þ . . . ðpnþ1Þð1þxÞpn Þx¼0 ¼ pðp 1Þ . . . ðp n þ 1Þ dxn x¼0 (1:141) On inserting these terms into (1.134), and noting that ƒ(0) = 1, we get for the binomial series pðp 1Þ 2 pðp 1Þðp 2Þ 3 x þ x þ 12 123 pðp 1Þ . . . ðp n þ 1Þ n x þ þ n!

ð1 þ xÞp ¼ 1 þ px þ

(1:142)

If the exponent p is not an integer, or p is negative, the series is convergent in the range –1 < x < 1.

1.10.4 Linear approximations The variations in some physical properties over the surface of the Earth are small in relation to the main property. For example, the difference between the polar radius c and the equatorial radius a expressed as a fraction of the equatorial radius deﬁnes the ﬂattening ƒ, which is equal to 1/298. This results from deformation of the Earth by the centrifugal force of its own rotation, which, expressed as a fraction m of the gravitational force, is equal to 1/289. Both ƒ and m are less than three thousandths of the main property, so ƒ2, m2, and the product fm are of the order of nine parts in a million and, along with higher-order combinations, are negligible. Curtailing the expansion of small quantities at ﬁrst order helps keep equations manageable without signiﬁcant loss of geophysical information. In the following chapters much use will be made of such linear approximation. It simpliﬁes the form of mathematical functions and the usable part of the series described above. For example, for small values of x or (x – x0), the following ﬁrst-order approximations may be used: sin x x; expðxÞ 1 þ x; ð1 þ xÞp 1 þ px

cos x 1 lnð1 þ xÞ x

fðxÞ fðx0 Þ þ ðx x0 Þ

df dx

(1:143)

x¼x0

32

Mathematical background

1.11 Leibniz’s rule Assume that u(x) and v(x) are differentiable functions of x. The derivative of their product is d dvðxÞ duðxÞ ðuðxÞvðxÞÞ ¼ uðxÞ þ vðxÞ dx dx dx

(1:144)

If we deﬁne the operator D = d/dx, we obtain a shorthand form of this equation: DðuvÞ ¼ u Dv þ v Du

(1:145)

We can differentiate the product (uv) a second time by parts, D2 ðuvÞ ¼ DðDðuvÞÞ ¼ u D2 v þ ðDuÞðDvÞ þ ðDvÞðDuÞ þ v D2 u ¼ u D2 v þ 2ðDuÞðDvÞ þ v D2 u

(1:146)

and, continuing in this way, D3 ðuvÞ ¼ u D3 v þ 3ðDuÞ D2 v þ 3 D2 u ðDvÞ þ v D3 u D4 ðuvÞ ¼ u D4 v þ 4ðDuÞ D3 v þ 6 D2 u D2 v þ 4 D2 u ðDvÞ þ v D4 u (1:147) The coefﬁcients in these equations are the binomial coefﬁcients, as deﬁned in (1.139). Thus after n differentiations we have Dn ðuvÞ ¼

n X

k nk n! D u D v k!ðn kÞ! k¼0

(1:148)

This relationship is known as Leibniz’s rule, after Gottfried Wilhelm Leibniz (1646–1716), who invented inﬁnitesimal calculus contemporaneously with Isaac Newton (1642–1727); each evidently did so independently of the other.

1.12 Legendre polynomials Let r and R be the sides of a triangle that enclose an angle θ and let u be the side opposite this angle (Fig. 1.11). The angle and sides are related by the cosine rule u2 ¼ r2 þ R2 2rR cos θ

(1:149)

Inverting this expression and taking the square root gives " ! 1 1 r ¼ 12 cos θ þ u R R

r R

!2 #1=2 (1:150)

1.12 Legendre polynomials

R

33

u

θ r

Fig. 1.11. Relationship of the sides r and R, which enclose an angle θ, and the side u opposite the angle, as used in the deﬁnition of Legendre polynomials.

Now let h ¼ r=R and x = cos θ, giving 1=2 1 1 1 ¼ 1 2xh þ h2 ¼ ð1 tÞ1=2 u R R

(1:151)

where t = 2xh – h2. The equation can be expanded as a binomial series 1 3 1 3 5 2 2 2 2 2 1 1=2 2 ð1tÞ ¼ 1þ ðtÞ þ ðtÞ þ ðtÞ3 þ 2 12 123 !2 !3 t 13 t 135 t ¼1þ þ þ þ 2 12 2 123 2 !n 1 3 5 . . . ð2n 1Þ t þ (1:152) þ 1 2 3...n 2 The inﬁnite series of terms on the right-hand side of the equation can be written ð1 tÞ1=2 ¼

1 X

an t n

(1:153)

n¼0

The coefﬁcient an is given by an ¼

1 3 5 . . . ð2n 1Þ 2n n!

(1:154)

Now, substitute the original expression for t,

1 2xh þ h2

1=2

¼

1 X n¼0

1 n X an 2xh h2 ¼ an hn ð2x hÞn

(1:155)

n¼0

This equation is an inﬁnite series in powers of h. The coefﬁcient of each term in the power series is a polynomial in x. Let the coefﬁcient of hn be Pn(x). The equation becomes 1 1=2 X Ψðx; hÞ ¼ 1 2xh þ h2 ¼ hn Pn ðxÞ n¼0

(1:156)

34

Mathematical background

Equation (1.156) is known as the generating function for the polynomials Pn(x). Using this result, and substituting h = r/R and x = cos θ, we ﬁnd that (1.151) becomes !n 1 1 1X r Pn ðcos θÞ (1:157) ¼ u R n¼0 R The polynomials Pn(x) or Pn(cos θ) are called Legendre polynomials, after the French mathematician Adrien-Marie Legendre (1752–1833). The deﬁning equation (1.157) is called the reciprocal-distance formula. An alternative formulation is given in Box 1.4.

1.13 The Legendre differential equation The Legendre polynomials satisfy an important second-order partial differential equation, which is called the Legendre differential equation. To derive this equation we will carry out a sequence of differentiations, starting with the generating function in the form 1=2 Ψ ¼ 1 2xh þ h2

(1:158)

Differentiating this function once with respect to h gives 3=2 ∂Ψ ¼ ðx hÞ 1 2xh þ h2 ¼ ðx hÞΨ3 ∂h

(1:159)

Differentiating Ψ twice with respect to x gives 3=2 ∂Ψ ¼ hΨ3 ¼ h 1 2xh þ h2 ∂x 1 ∂Ψ Ψ3 ¼ h ∂x ∂2 Ψ ∂Ψ ¼ 3h2 Ψ5 ¼ 3hΨ2 2 ∂x ∂x 1 ∂2 Ψ Ψ5 ¼ 2 2 3h ∂x

(1:160)

(1:161)

Next we perform successive differentiations of the product (hΨ) with respect to h. The ﬁrst gives ∂ ∂Ψ ðhΨÞ ¼ Ψ þ h ¼ Ψ þ hðx hÞΨ3 ∂h ∂h

(1:162)

1.13 The Legendre differential equation

35

Box 1.4. Alternative form of the reciprocal-distance formula The sides and enclosed angle of the triangle in Fig. 1.11 are related by the cosine rule u2 ¼ r2 þ R2 2rR cos θ

(1)

Instead of taking R outside the brackets as in (1.150), we can move r outside and write the expression for u as " 2 #1=2 1 1 R R (2) ¼ 12 cos θ þ u r r r Following the same treatment as in Section 1.12, but now with h ¼ R=r and x = cos θ, we get 1=2 1 1 1 ¼ ð1 tÞ1=2 ¼ 1 2xh þ h2 u r r

(3)

where t = 2xh – h2. The function (1 – t)−1/2 is expanded as a binomial series, which again gives an inﬁnite series in h, in which the coefﬁcient of hn is Pn(x). The deﬁning equation is as before: 1 1=2 X ¼ hn Pn ðxÞ Ψðx; hÞ ¼ 1 2xh þ h2

(4)

n¼0

On substituting h ¼ R=r and x = cos θ, we ﬁnd an alternative form for the generating equation for the Legendre polynomials: 1 n 1 1X R Pn ðcos θÞ (5) ¼ u r n¼0 r

On repeating the differentiation, and taking (1.159) into account, we get ∂2 ∂Ψ ∂Ψ ðhΨÞ ¼ þ hðx hÞ3Ψ2 þ ðx 2hÞΨ3 ∂h ∂h ∂h2 ¼ ðx hÞΨ3 þ 3hðx hÞ2 Ψ5 þ ðx 2hÞΨ3 ¼ ð2x 3hÞΨ3 þ 3hðx hÞ2 Ψ5 Now substitute for Ψ3, from (1.160), and Ψ5, from (1.161), giving

(1:163)

36

Mathematical background ∂2 1 ∂Ψ 1 ∂2 Ψ 2 ð hΨ Þ ¼ ð 2x 3h Þ þ 3h ð x h Þ ∂h2 h ∂x 3h2 ∂x2

(1:164)

Multiply throughout by h: 2 ∂2 ∂Ψ 2 ∂ Ψ h 2 ðhΨÞ ¼ ð2x 3hÞ þ ðx hÞ ∂h ∂x ∂x2 2 ∂Ψ ∂Ψ ∂ Ψ ¼ 2x 3h þ ð x hÞ 2 ∂x ∂x ∂x2

(1:165)

The second term on the right can be replaced as follows, again using (1.160) and (1.161): 3h

∂Ψ 1 ∂2 Ψ ¼ 3h2 Ψ3 ¼ 2 2 ∂x Ψ ∂x ∂2 Ψ ¼ 1 2xh þ h2 ∂x2

On substituting into (1.165) and gathering terms, we get h i ∂2 Ψ ∂2 ∂Ψ þ 2x h 2 ðhΨÞ ¼ ðx hÞ2 1 2xh þ h2 ∂x2 ∂x ∂h ∂2 Ψ ∂2 ∂Ψ h 2 ðhΨÞ ¼ x2 1 þ 2x ∂h ∂x2 ∂x

(1:166)

(1:167)

(1:168)

The Legendre polynomials Pn(x) are deﬁned in (1.156) as the coefﬁcients of hn in the expansion of Ψ as a power series. On multiplying both sides of (1.156) by h, we get hΨ ¼

1 X

hnþ1 Pn ðxÞ

(1:169)

n¼0

We differentiate this expression twice and multiply by h to get a result that can be inserted on the left-hand side of (1.168):

h

1 X ∂ ðn þ 1Þhn Pn ðxÞ ðhΨÞ ¼ ∂h n¼0

(1:170)

1 X ∂2 ðhΨÞ ¼ nðn þ 1Þhn Pn ðxÞ 2 ∂h n¼0

(1:171)

Using (1.156), we can now eliminate Ψ and convert (1.168) into a second-order differential equation involving the Legendre polynomials Pn(x),

1.13 The Legendre differential equation

37

Table 1.1. Some ordinary Legendre polynomials of low degree n

Pn(x)

Pn(cos θ)

0 1

1 x 1 2 3x 1 2 1 3 5x 3x 2 1 35x4 30x2 þ 3 8

1 cos θ 1 3 cos2 θ 1 2 1 5 cos3 θ 3 cos θ 2 1 35 cos4 θ 30 cos2 θ þ 3 8

2 3 4

1 X

hn

n¼0 1 X n¼0

X 1 2 d 2 Pn ðxÞ dPn ðxÞ x 1 þ 2x nðn þ 1Þhn Pn ðxÞ (1:172) ¼ dx2 dx n¼0

hn

2 d 2 Pn ðxÞ dPn ðxÞ x 1 þ 2x ðxÞ ¼0 nðn þ 1ÞPn dx2 dx

(1:173)

If this expression is true for every non-zero value of h, the quantity in curly brackets must be zero, thus d 2 Pn ðxÞ dPn ðxÞ 2x 1 x2 þ nðn þ 1ÞPn ðxÞ ¼ 0 2 dx dx

(1:174)

An alternative, simpler form for this equation is obtained by combining the ﬁrst two terms:

d 2 dPn ðxÞ (1:175) 1x þ nðn þ 1ÞPn ðxÞ ¼ 0 dx dx This is the Legendre differential equation. It has a family of solutions, each of which is a polynomial corresponding to a particular value of n. The Legendre polynomials provide solutions in potential analyses with spherical symmetry, and have an important role in geophysical theory. Some Legendre polynomials of low degree are listed in Table 1.1.

1.13.1 Orthogonality of the Legendre polynomials Two vectors a and b are orthogonal if their scalar product is zero: a · b ¼ ax bx þ ay by þ az bz ¼

3 X i¼1

ai bi ¼ 0

(1:176)

38

Mathematical background

By analogy, two functions of the same variable are said to be orthogonal if their product, integrated over a particular range, is zero. For example, the trigonometric functions sin θ and cos θ are orthogonal for the range 0 ≤ θ ≤ 2π, because Z2π

Z2π sin θ cos θ dθ ¼

θ¼0

θ¼0

2π 1 1 sinð2θÞdθ ¼ cosð2θÞ ¼ 0 2 4 θ¼0

(1:177)

The Legendre polynomials Pn(x) and Pl(x) are orthogonal over the range –1 ≤ x ≤ 1. This can be established as follows. First, we write the Legendre equation in short form, dropping the variable x for both Pn and Pl, and, for brevity, writing d Pn ðxÞ ¼ P0n dx

and

d2 Pn ðxÞ ¼ P00n dx2

(1:178)

Thus

1 x2 P00n 2xP0n þ nðn þ 1ÞPn ¼ 0 1 x2 P00l 2xP0l þ lðl þ 1ÞPl ¼ 0

(1:179) (1:180)

Multiplying (1.179) by Pl and (1.180) by Pn gives

1 x2 Pl P00n 2xPl P0n þ nðn þ 1ÞPl Pn ¼ 0 1 x2 Pn P00l 2xP0l Pn þ lðl þ 1ÞPl Pn ¼ 0

(1:181) (1:182)

Subtracting (1.182) from (1.181) gives

1 x2 Pl P00n Pn P00l 2x Pl P0n P0l Pn þ ½nðn þ 1Þ lðl þ 1ÞPl Pn ¼0 (1:183)

Note that d Pl P0n P0l Pn ¼ Pl P00n þ P0l P0n P0l P0n P00l Pn ¼ Pl P00n P00l Pn dx (1:184) and

d Pl P0n P0l Pn 2x Pl P0n P0l Pn dx d 1 x2 Pl P0n P0l Pn ¼ dx

1 x2

(1:185)

1.13 The Legendre differential equation

39

Thus d 1 x2 Pl P0n P0l Pn þ ½nðn þ 1Þ lðl þ 1ÞPl Pn ¼ 0 dx

(1:186)

Now integrate each term in this equation with respect to x over the range –1 ≤ x ≤ 1. We get

1x

2

Pl P0n

P0l Pn

þ1

x¼1

Zþ1 þ ½nðn þ 1Þ lðl þ 1Þ

Pl Pn dx ¼ 0 x¼1

(1:187) The ﬁrst term is zero on evaluation of (1 – x ) at x = ±1; thus the second term must also be zero. For n ≠ l the condition for orthogonality of the Legendre polynomials is Z þ1 Pn ðxÞPl ðxÞdx ¼ 0 (1:188) 2

x¼1

1.13.2 Normalization of the Legendre polynomials A function is said to be normalized if the integral of the square of the function over R þ1 its range is equal to 1. Thus we must evaluate the integral x¼1 ½Pn ðxÞ2 dx. We begin by recalling the generating function for the Legendre polynomials given in (1.156), which we rewrite for Pn(x) and Pl(x) individually: 1 X

1=2 hn Pn ðxÞ ¼ 1 2xh þ h2

(1:189)

1=2 hl Pl ðxÞ ¼ 1 2xh þ h2

(1:190)

n¼0 1 X l¼0

Multiplying these equations together gives 1 X 1 X

1 hnþl Pn ðxÞPl ðxÞ ¼ 1 2xh þ h2

(1:191)

l¼0 n¼0

Now let l = n and integrate both sides with respect to x, taking into account (1.188): 1 X n¼0

Zþ1 h

Zþ1 2

½Pn ðxÞ dx ¼

2n x¼1

x¼1

dx 1 þ h2 2xh

(1:192)

40

Mathematical background

The right-hand side of this equation is a standard integration that results in a natural logarithm: Z dx 1 ¼ lnða þ bxÞ (1:193) a þ bx b The right-hand side of (1.192) therefore leads to þ1 dx 1 2 ¼ ln 1 þ h 2xh 1 þ h2 2xh ð2hÞ x¼1

Zþ1 x¼1

¼

1 ln 1 þ h2 2h ln 1 þ h2 þ 2h 2h

(1:194)

and Zþ1 x¼1

dx 1 1 2 1 2 þ ¼ ln ð 1 h Þ ln ð 1 þ h Þ 1 þ h2 2xh h 2 2 1 ¼ ½lnð1 þ hÞ lnð1 hÞ h

(1:195)

Using the MacLaurin series for the natural logarithms as in (1.135), we get lnð1 þ hÞ ¼ h

h2 h3 h4 hn þ þ ð1Þn1 þ 2 3 4 n

lnð1 hÞ ¼ h

h2 h3 h4 ðhÞn ð1Þn1 þ 2 3 4 n

(1:196) (1:197)

Subtracting the second equation from the ﬁrst gives Zþ1 x¼1

1 dx 2 h3 h5 2X h2nþ1 ¼ þ þ ¼ h þ 2 3 5 1 þ h 2xh h h n¼0 2n þ 1

(1:198)

Inserting this result into (1.192) gives 1 X n¼0 1 X n¼0

0 h2n @

Zþ1 x¼1

Zþ1 h

½Pn ðxÞ2 dx ¼

2n x¼1

1 2X h2nþ1 h n¼0 2n þ 1

1 2 A¼0 ½Pn ðxÞ2 dx 2n þ 1

(1:199)

(1:200)

1.14 Rodrigues’ formula

41

This is true for every value of h in the summation, so we obtain the normalizing condition for the Legendre polynomials: Zþ1 ½Pn ðxÞ2 dx ¼ x¼1

2 2n þ 1

(1:201)

1=2 It follows that n þ 12 Pn ðxÞ is a normalized Legendre polynomial.

1.14 Rodrigues’ formula The Legendre polynomials can be easily computed with the aid of a formula derived by a French mathematician, Olinde Rodrigues (1795–1851). First, we deﬁne the function n fðxÞ ¼ x2 1 (1:202) Differentiating ƒ(x) once with respect to x gives n n1 df d 2 ¼ x 1 ¼ 2nx x2 1 dx dx

(1:203)

Multiplying the result by (x2 − 1) gives 2 d 2 n n x 1 x 1 ¼ 2nx x2 1 dx

(1:204)

2 df ¼ 2nxf x 1 dx

(1:205)

Now we use Leibniz’s rule (1.144) to differentiate both sides of this equation n + 1 times with respect to x. Writing D = d/dx as in Section 1.11, Dnþ1 ðuvÞ ¼

nþ1 X

ðn þ 1Þ! k nþ1k v D u D k!ðn þ 1 kÞ! k¼0

(1:206)

On the left-hand side of (1.205) let u(x) = (x2 − 1) and v(x) = dƒ/dx = Dƒ. Applying Leibniz’s rule, we note that after only three differentiations of (x2 − 1) the result is zero and the series is curtailed. On the right-hand side let u(x) = 2nx and v(x) = ƒ. Note that in this case the series is curtailed after two differentiations. Thus, using Leibniz’s rule to differentiate each side of (1.205) n + 1 times, we get

42

Mathematical background

ðn þ 1Þn n x2 1 Dnþ2 f þ 2xðn þ 1ÞDnþ1 f þ 2 D f ¼ 2nx Dnþ1 f 1·2 þ 2nðn þ 1ÞDn f (1:207)

On gathering terms and bringing all to the left-hand side, we have 2 x 1 Dnþ2 f þ 2x Dnþ1 f nðn þ 1ÞDn f ¼ 0

(1:208)

Now we deﬁne y(x) such that yðxÞ ¼ Dn f ¼

n dn 2 x 1 n dx

(1:209)

and we have

d 2y

x2 1

dx2

þ 2x

dy nðn þ 1Þy ¼ 0 dx

(1:210)

On comparing with (1.174), we see that this is the Legendre equation. The Legendre polynomials must therefore be proportional to y(x), so we can write Pn ðxÞ ¼ cn

n dn 2 x 1 n dx

(1:211)

The quantity cn is a calibration constant. To determine cn we ﬁrst write n d n dn 2 x 1 ¼ n ½ ð x 1Þ n ð x þ 1Þ n dxn dx

(1:212)

then we apply Leibniz’s rule to the product on the right-hand side of the equation: n n X dn 2 n! dm d nm x 1 ¼ ðx 1Þn nm ðx þ 1Þn n m dx dx m!ðn mÞ! dx m¼0

(1:213)

The successive differentiations of (x – 1)n give d ðx 1Þn ¼ nðx 1Þn1 dx d2 ðx 1Þn ¼ nðn 1Þðx 1Þn2 dx2 d n1 ðx 1Þn ¼ ðnðn 1Þðn 2Þ . . . 3·2·1Þðx 1Þ ¼ n!ðx 1Þ dxn1 dn ðx 1Þn ¼ n! dxn

(1:214)

1.15 Associated Legendre polynomials

43

Each differentiation in (1.214) is zero at x = 1, except the last one. Thus each term in the sum in (1.213) is also zero except for the last one, for which m = n. Substituting x = 1 gives n

n n d 2 nd n x 1 ¼ ð x þ 1 Þ ð x 1 Þ ¼ 2n n! (1:215) n dxn dx x¼1 x¼1 Putting this result and the condition Pn(1) = 1 into (1.211) gives n n d 2 x 1 ¼ cn 2n n! ¼ 1 Pn ð1Þ ¼ cn dxn x¼1

(1:216)

where cn ¼

1 2n n!

(1:217)

Rodrigues’ formula for the Legendre polynomials is therefore Pn ðxÞ ¼

1 dn 2n n! dxn

x2 1

n

(1:218)

1.15 Associated Legendre polynomials Many physical properties of the Earth, such as its magnetic ﬁeld, are not azimuthally symmetric about the rotation axis when examined in detail. However, these properties can be described using mathematical functions that are based upon the Legendre polynomials described in the preceding section. To derive these functions, we start from the Legendre equation, (1.174), which can be written in shorthand form as 1 x2 P 00n 2xP0n þ nðn þ 1ÞPn ¼ 0 (1:219) Now we differentiate this equation with respect to x: d 00 d d P 2xP00n 2x P0n 2P0n þ nðn þ 1Þ Pn ¼ 0 1 x2 dx n dx dx

(1:220)

On noting that we can equally write P 00n ¼ ðd=dxÞP0n and P0n ¼ ðd=dxÞPn , this can be written alternatively as

1 x2

d 00 d d P 4x P0n þ ½nðn þ 1Þ 2 Pn ¼ 0 dx n dx dx

which can be written, for later comparison,

(1:221)

44

Mathematical background d 00 d d Pn 2ð2Þx P0n þ ½nðn þ 1Þ 1ð2Þ Pn ¼ 0 1 x2 dx dx dx

(1:222)

Next, we differentiate this expression again, observing the same rules and gathering terms,

d 2 00 d 00 d 0 d2 0 d2 1x P 2x þ 4x P Pn ¼ 0 4 þ ½ n ð n þ 1 Þ2 P P n n n n dx2 dx2 dx2 dx dx 2

(1:223) d 2 00 d2 0 d2 0 d2 P 2x P 4x P þ ½ n ð n þ 1 Þ 2 4 Pn ¼ 0 1 x2 dx2 n dx2 n dx2 n dx2 (1:224)

1 x2

d 2 00 d2 0 d2 P 6x P þ ½ n ð n þ 1 Þ 6 Pn ¼ 0 n n dx2 dx2 dx2

(1:225)

which, as we did with (1.222), we can write in the form

1 x2

d 2 00 d2 d2 Pn 2ð3Þx 2 P0n þ ½nðn þ 1Þ 2ð3Þ 2 Pn ¼ 0 2 dx dx dx

(1:226)

On following step-by-step in the same manner, we get after the third differentiation

1 x2

d 3 00 d3 0 d3 P 2 ð 4 Þx P þ ½ n ð n þ 1 Þ 3 ð 4 Þ Pn ¼ 0 n n dx3 dx3 dx3

(1:227)

Equations (1.222), (1.226), and (1.227) all have the same form. The higherorder differentiation is accompanied by systematically different constants. By extension, differentiating (1.219) m times (where m ≤ n) yields the differential equation

1 x2

d m 00 dm dm Pn 2ðm þ 1Þx m P0n þ ½nðn þ 1Þ mðm þ 1Þ m Pn ¼ 0 m dx dx dx (1:228)

Now let the mth-order differentiation of Pn be written as dm QðxÞ Pn ðxÞ ¼ m=2 dxm ð1 x2 Þ

(1:229)

Substitution of this expression into (1.228) gives a new differential equation involving Q(x). We need to determine both ðd m =dxm ÞP0n and ðd m =dxm ÞP00n , so ﬁrst we differentiate (1.229) with respect to x:

1.15 Associated Legendre polynomials

dm 0 Q0 P ¼ n m=2 dxm ð1 x2 Þ

! m Q ð2xÞ m=2þ1 2 ð1 x 2 Þ

ðmþ2Þ=2 dm 0 P n ¼ 1 x2 1 x2 Q0 þ mxQ m dx

45

(1:230)

(1:231)

A further differentiation of (1.231) by parts gives ðmþ2Þ=2 d 1 x2 1 x2 Q0 þ mxQ dx ðmþ2Þ=2 d 1 x2 Q0 þ mxQ þ 1 x2 dx 2 ðmþ2Þ=21 ¼ ðm þ 2Þx 1 x 1 x2 Q0 þ mxQ ðmþ2Þ=2 1 x2 Q00 þ mxQ0 2xQ0 þ mQ þ 1 x2

d m 00 P ¼ dxm n

(1:232) n d m 00 2 ðmþ2Þ=2 P ¼ 1 x 1 x2 Q00 þ ðm 2ÞxQ0 þ mQ þ ðm þ 2ÞxQ0 n dxm mðm þ 2Þx2 Q (1:233) þ 1 x2 00 d m 00 mðm þ 2Þx2 Q 2 ðmþ2Þ=2 2 0 P ¼ 1 x 1 x þ 2mxQ þ mQ þ Q dxm n 1 x2 (1:234)

Now we substitute (1.231) and (1.234) into (1.228). Unless the multiplier (1 – x2)− (m + 2)/2 is always zero, Q must satisfy the following equation: 2 1 x2 Q00 þ 2mx 1 x2 Q0 þ m 1 x2 Q þ mðm þ 2Þx2 Q 2ðm þ 1Þx 1 x2 Q0 2mðm þ 1Þx2 Q þ ½nðn þ 1Þ mðm þ 1Þ 1 x2 Q ¼ 0

(1:235)

The remainder of the evaluation consists of gathering and reducing terms; we ﬁnally get

m2 Q¼0 1 x Q 2xQ þ nðn þ 1Þ 1 x2 2

00

0

(1:236)

46

Mathematical background

The functions Q(x) involve two parameters, the degree n and order m, and are written Pn,m(x). Thus

d2 d m2 Pn;m ðxÞ ¼ 0 P 1 x2 P ð x Þ 2x ð x Þ þ n ð n þ 1 Þ n;m n;m dx2 1 x2 dx (1:237) This is the associated Legendre equation. The solutions Pn,m(x) or Pn,m(cos θ), where x = cos θ, are called associated Legendre polynomials, and are obtained from the ordinary Legendre polynomials using the deﬁnition of Q in (1.229): m=2 d m Pn;m ðxÞ ¼ 1 x2 Pn ðxÞ dxm

(1:238)

Substituting Rodrigues’ formula (1.218) for Pn(x) into this equation gives

m=2 n 1 x2 d nþm 2 x 1 Pn;m ðxÞ ¼ n nþm 2 n! dx

(1:239)

The highest power of x in the function (x2 − 1)n is x2n. After 2n differentiations the result will be a constant, and a further differentiation will give zero. Therefore n + m ≤ 2n, and possible values of m are limited to the range 0 ≤ m ≤ n.

1.15.1 Orthogonality of associated Legendre polynomials For succinctness we again write Pn,m(x) as simply Pn,m. The deﬁning equations for the associated Legendre polynomials Pn,m and Pl,m are

0 m2 1 x Pn;m 2x Pn;m þ nðn þ 1Þ Pn;m ¼ 0 1 x2

00 0 m2 2 1 x Pl;m 2x Pl;m þ lðl þ 1Þ Pl;m ¼ 0 1 x2 2

00

(1:240)

(1:241)

As for the ordinary Legendre polynomials, we multiply (1.240) by Pl,m and (1.241) by Pn,m:

00 0 m2 Pn;m Pl;m ¼ 0 1 x2 Pn;m Pl;m 2x Pn;m Pl;m þ nðn þ 1Þ 1 x2 (1:242)

1.15 Associated Legendre polynomials

47

00 0 m2 1 x2 Pl;m Pn;m 2x Pl;m Pn;m þ lðl þ 1Þ Pn;m Pl;m ¼ 0 1 x2 (1:243)

On subtracting (1.243) from (1.242) we have i h i h 00 00 0 0 1 x2 Pn;m Pl;m Pl;m Pn;m 2x Pn;m Pl;m Pl;m Pn;m þ ½nðn þ 1Þ lðl þ 1ÞPn;m Pl;m ¼ 0

(1:244)

Following the method used to establish the orthogonality of the ordinary Legendre polynomials (Section 1.13.1), we can write this equation as

o 0 0 d n 1 x2 Pn;m Pl;m Pl;m Pn;m þ ½nðn þ 1Þ lðl þ 1ÞPn;m Pl;m dx ¼0 (1:245) On integrating each term with respect to x over the range –1 ≤ x ≤ 1, we get n

oþ1 0 0 1 x2 Pn;m Pl;m Pl;m Pn;m

x¼1

Zþ1 þ ½nðn þ 1Þ lðl þ 1Þ

Pn;m Pl;m dx ¼ 0

(1:246)

x¼1

The ﬁrst term is zero on evaluation of (1 – x2) at x = ±1; thus the second term must also be zero. Provided that n ≠ l, the condition of orthogonality of the associated Legendre polynomials is x¼þ1 Z

Pn;m ðxÞPl;m ðxÞdx ¼ 0

(1:247)

x¼1

1.15.2 Normalization of associated Legendre polynomials Squaring the associated Legendre polynomials and integrating over –1 ≤ x ≤ 1 gives x¼þ1 Z

x¼1

2 Pn;m ðxÞ dx ¼

2 ðn þ mÞ! 2n þ 1 ðn mÞ!

(1:248)

48

Mathematical background

The squared functions do not integrate to 1, so they are not normalized. If each polynomial is multiplied by a normalizing function, the integrated squared polynomial can be made to equal a chosen value. Different conditions for this apply in geodesy and geomagnetism. The Legendre polynomials used in geodesy are fully normalized. They are deﬁned as follows: Pm n ðxÞ ¼

2n þ 1 ðn mÞ! 1=2 Pn;m ðxÞ 2 ðn þ mÞ!

(1:249)

The Legendre polynomials used in geomagnetism are partially normalized (or quasi-normalized). Schmidt in 1889 deﬁned this method of normalization so that Pm n ðxÞ

¼

ðn mÞ! 2 ðn þ mÞ!

P0n ðxÞ ¼ Pn;0 ðxÞ;

1=2 Pn;m ðxÞ; m¼0

m 6¼ 0

(1:250) (1:251)

Integration of the squared Schmidt polynomials over the full range –1 ≤ x ≤ 1 gives the value 1 for m = 0 and 1/(2n + 1) for m > 0. Some fully normalized Legendre polynomials and partially normalized Schmidt polynomials are listed in Table 1.2. Table 1.2. Some fully normalized associated Legendre polynomials and partially normalized Schmidt polynomials of low degree and order Pm n ðcos θ Þ; Legendre, fully normalized

Pm n ðcos θ Þ; Schmidt, partially normalized

1

cos θ sin θ 1 ð3 cos2 θ 1Þ 2 3 sin θ cos θ

2

2

3 sin2θ

3

0

3

1

3

2

15 sin2θ cos θ

3

3

15 sin3θ

cos θ sin θ 1 ð3 cos2 θ 1Þ 2pﬃﬃﬃ pﬃﬃ3ﬃ sin θ cos θ 3 sin2 θ 2 1 cos θð5 cos2 θ 3Þ 2pﬃﬃﬃ 6 sin θð5 cos2 θ 1Þ 4 pﬃﬃﬃﬃﬃ 15 15 sin2 θ cos θ p2ﬃﬃﬃﬃﬃ 10 sin3 θ 4

n

m

1 1

0 1

2

0

2

1 cos θð5 cos2 θ 3Þ 2 3 sin θð5 cos2 θ 1Þ 2

1.16 Spherical harmonic functions

49

1.16 Spherical harmonic functions Several geophysical potential ﬁelds – for example, gravitation and geomagnetism – satisfy the Laplace equation. Spherical polar coordinates are best suited for describing a global geophysical potential. The potential can vary with distance r from the Earth’s center and with polar angular distance θ and azimuth (equivalent to co-latitude and longitude in geographic terms) on any concentric spherical surface. The solution of Laplace’s equation in spherical polar coordinates for a potential U may be written (see Section 2.4.5, (2.104)) U¼

1 X n X m Bn m An rn þ nþ1 am n cosðmÞ þ bn sinðmÞ Pn ðcos θÞ (1:252) r n¼0 m¼0

m Here An, Bn, am n , and bn are constants that apply to a particular situation. On the surface of the Earth, or an arbitrary sphere, the radial part of the potential of a point source at the center of the sphere has a constant value and the variation over the surface of the sphere is described by the functions in θ and . We are primarily interested in solutions outside the Earth, for which An is zero. Also we can set the constant Bn equal to Rn+1, where R is the Earth’s mean radius. The potential is then given by 1 X n nþ1 X m R m am U¼ n cosðmÞ þ bn sinðmÞ Pn ðcos θÞ r n¼0 m¼0

(1:253)

m Let the spherical harmonic functions Cm n ðθ; Þ and Sn ðθ; Þ be deﬁned as m Cm n ðθ; Þ ¼ cosðmÞ Pn ðcos θÞ m Sm n ðθ; Þ ¼ sinðmÞ Pn ðcos θÞ

(1:254)

The variation of the potential over the surface of a sphere may be described by these functions, or a more general spherical harmonic function Ym n ðθ; Þ that combines the sine and cosine variations: m Ym n ðθ; Þ ¼ Pn ðcos θ Þ

cosðmÞ sinðmÞ

(1:255)

Like their constituent parts – the sine, cosine, and associated Legendre functions – spherical harmonic functions are orthogonal and can be normalized.

50

Mathematical background

1.16.1 Normalization of spherical harmonic functions m Normalization of the functions Cm n ðθ; Þ and Sn ðθ; Þ requires integrating the squared value of each function over the surface of a unit sphere. The element of surface area on a unit sphere is dΩ = sin θ dθ d (Box 1.3) and the limits of integration are 0 ≤ θ ≤ π and 0 ≤ ≤ 2π. The integral is

ZZ

2 Cm n ðθ; Þ

Zπ

Z2π

dΩ ¼

Cm n ðθ; Þ

2

sin θ dθ d

θ¼0 ¼0

S

Zπ

Z2π

¼

cosðmÞPm n ðcos θ Þ

2

sin θ dθ d (1:256)

θ¼0 ¼0

Let x = cos θ in the associated Legendre polynomial, so that dx = –sin θ dθ and the limits of integration are –1 ≤ x ≤ 1. The integration becomes 8 Z1 > < Z2π x¼1

> :

¼0

9 > Zþ1 = 2 m 2 2 m cos ðmÞd Pn ðxÞ dx ¼ π Pn ðxÞ dx > ;

(1:257)

x¼1

Normalization of the associated Legendre polynomials gives the result in (1.248), thus ZZ

Cm n ðθ; Þ

2

dΩ ¼

S

2π ðn þ mÞ! 2n þ 1 ðn mÞ!

(1:258)

The normalization of the function Sm n ðθ; Þ by this method delivers the same result. Spherical harmonic functions make it possible to express the variation of a physical property (e.g., gravity anomalies, g(θ, )) on the surface of the Earth as an inﬁnite series, such as gðθ; Þ ¼

1 X n X

m m m am n Cn ðθ; Þ þ bn Sn ðθ; Þ

(1:259)

n¼0 m¼0 m The coefﬁcients am n and bn may be obtained by multiplying the function gðθ; Þ m m by Cn ðθ; Þ or Sn ðθ; Þ, respectively, and integrating the product over the surface of the unit sphere. The normalization properties give

1.16 Spherical harmonic functions

51

ZZ 2n þ 1 ðn mÞ! gðθ; Þ Cm n ðθ; ÞdΩ 2π ðn þ mÞ! S ZZ 2n þ 1 ð n m Þ! bm gðθ; Þ Sm n ¼ n ðθ; ÞdΩ 2π ðn þ mÞ!

am n ¼

(1:260)

S

1.16.2 Zonal, sectorial, and tesseral spherical harmonics The spherical harmonic functions Ym n ðθ; Þ have geometries that allow graphic representation of a potential on the surface of a sphere. Deviations of the potential from a constant value form alternating regions in which the potential is larger or smaller than a uniform value. Where the potential surface intersects the spherical surface a nodal line is formed. The appearance of any Ym n ðθ; Þ is determined by ð the distribution of its nodal lines. These occur where Ym n θ; Þ = 0. To simplify the discussion we will associate a constant value of the polar angle θ with a circle of latitude, and a constant value of the azimuthal angle with a circle of longitude. The deﬁnition of the associated Legendre polynomials in (1.239) shows that the equation Pm n ðxÞ = 0 has n – m roots, apart from the trivial solution x = ±1. The variation of the spherical harmonic Ym n ðθ; Þ with latitude θ thus has n – m nodal lines, each a circle of latitude, between the two poles. If additionally m = 0, the potential on the sphere varies only with latitude and there are n nodal lines separating zones in which the potential is greater or less than the uniform value. An example of a zonal spherical harmonic is Y02 ðθ; Þ, shown in Fig. 1.12(a). The solution of Laplace’s equation (1.253) shows that the variation in potential around any circle of latitude is described by the function m ΦðÞ ¼ am n cosðmÞ þ bn sinðmÞ

(1:261)

There are 2m nodal lines where Φ() = 0, corresponding to 2m meridians of longitude, or m great circles. In the special case in which n = m, there are no

(a) zonal, Y20

(b) sectorial, Y55

(c) tesseral, Y54

Fig. 1.12. Appearance of (a) zonal, (b) sectorial, and (c) tesseral spherical harmonics, projected on a meridian plane of the reference sphere.

52

Mathematical background

nodal lines of latitude and the longitudinal lines separate sectors in which the potential is greater or less than the uniform value. An example of a sectorial spherical harmonic is Y55 ðθ; Þ, shown in Fig. 1.12(b). In the general case (m ≠ 0, n ≠ m) the potential varies with both latitude and longitude. There are n – m nodal lines of latitude and m nodal great circles (2m meridians) of longitude. The appearance of the spherical harmonic resembles a patchwork of alternating regions in which the potential is greater or less than the uniform value. An example of a tesseral spherical harmonic is Y45 ðθ; Þ, which is shown in Fig. 1.12(c).

1.17 Fourier series, Fourier integrals, and Fourier transforms 1.17.1 Fourier series Analogously to the representation of a continuous function by a power series (Section 1.10), it is possible to represent a periodic function by an inﬁnite sum of terms consisting of the sines and cosines of harmonics of a fundamental frequency. Consider a periodic function ƒ(t) with period τ that is deﬁned in the interval 0 ≤ t ≤ τ, so that (a) ƒ(t) is ﬁnite within the interval; (b) ƒ(t) is periodic outside the interval, i.e., ƒ(t + τ) = ƒ(t); and (c) ƒ(t) is single-valued in the interval except at a ﬁnite number of points, and is continuous between these points. Conditions (a)–(c) are known as the Dirichlet conditions. If they are satisﬁed, ƒ(t) can be represented as fðtÞ ¼

1 a0 X þ ðan cosðnωtÞ þ bn sinðnωtÞÞ 2 n¼1

(1:262)

where ω = 2π/τ and the factor 12 in the ﬁrst term is included for reasons of symmetry. This representation of ƒ(t) is known as a Fourier series. The orthogonal properties of sine and cosine functions allow us to ﬁnd the coefﬁcients an and bn of the nth term in the series by multiplying (1.262) by sin(nωt) or cos(nωt) and integrating over a full period: 2 an ¼ τ

Zτ=2 fðtÞcosðnωtÞdt t¼τ=2

bn ¼

2 τ

(1:263)

Zτ=2 fðtÞsinðnωtÞdt t¼τ=2

1.17 Fourier series, integrals, and transforms

53

Instead of using trigonometric functions, we can replace the sine and cosine terms with complex exponentials using the deﬁnitions in (1.7), i.e., we write expðinωtÞ þ expðinωtÞ 2 expðinωtÞ expðinωtÞ sinðnωtÞ ¼ 2i

cosðnωtÞ ¼

(1:264)

Using these relationships in (1.262) yields bn ½expðinωtÞ expðinωtÞ 2 2i n¼0 1 1 X an ibn X an þ ibn ¼ expðinωtÞ þ expðinωtÞ (1:265) 2 2 n¼0 n¼0

fðtÞ ¼

1 X an

½expðinωtÞ þ expðinωtÞ þ

The summation indices are dummy variables, so in the second sum we can replace n by –n, and extend the limits of the sum to n = –∞; thus fðtÞ ¼ ¼

1 X an ibn n¼0 1 X

2

0 X an þ ibn expðinωtÞ þ expðinωtÞ 2 n¼1

ðan þ an Þ iðbn bn Þ expðinωtÞ 2 n¼1

(1:266)

If we deﬁne cn as the complex number cn ¼

ðan þ an Þ iðbn bn Þ 2

(1:267)

the Fourier series (1.262) can be written in complex exponential form as fðtÞ ¼

1 X

cn expðinωtÞ

(1:268)

n¼1

In this case the harmonic coefﬁcients cn are given by 1 cn ¼ τ

Zτ=2 fðtÞexpðinωtÞdt t¼τ=2

(1:269)

54

Mathematical background

1.17.2 Fourier integrals and Fourier transforms A Fourier series represents the periodic behavior of a physical property as an inﬁnite set of discrete frequencies. The theory can be extended to represent a function ƒ(t) that is not periodic and is made up of a continuous spectrum of frequencies, provided that the function satisﬁes the Dirichlet conditions speciﬁed above and that it has a ﬁnite energy: Z1

j fðtÞj2 dt51

(1:270)

t¼1

The inﬁnite sum in (1.268) is replaced by a Fourier integral and the complex coefﬁcients cn are replaced by an amplitude function g(ω): Z1 fðtÞ ¼

gðωÞexpðiωtÞdω

(1:271)

ω¼1

where g(ω) is a continuous function, obtained from the equation 1 gðωÞ ¼ 2π

Z1 fðtÞexpðiωtÞdt

(1:272)

t¼1

The transition from Fourier series to Fourier integral is explained in Box 1.5. The function g(ω) is called the forward Fourier transform of ƒ(t), and ƒ(t) is called the inverse Fourier transform of g(ω). Fourier transforms constitute a powerful mathematical tool for transforming a function ƒ(t) that is known in the time domain into a new function g(ω) in the frequency domain.

1.17.3 Fourier sine and cosine transforms A simple but important characteristic of a function is whether it is even or odd. An even function has the same value for both positive and negative values of its argument, i.e., ƒ(–t) = ƒ(t). The cosine of an angle is an example of an even function. The integral of an even function over a symmetric interval about the origin is equal to twice the integral of the function over the positive argument. The sign of an odd function changes with that of the argument, i.e., ƒ(–t) = –ƒ(t). For example, the sine of an angle is an odd function. The integral of an odd function over a symmetric interval about the origin is zero. The product of two odd functions or two even functions is an even function; the product of an odd function and an even function is an odd function.

1.17 Fourier series, integrals, and transforms

55

Box 1.5. Transition from Fourier series to Fourier integral The complex exponential Fourier series for a function ƒ(t) is 1 X

fðtÞ ¼

cn expðinωtÞ

(1)

n¼1

where the complex coefﬁcients cn are given by Zτ=2

1 cn ¼ τ

fðtÞexpðinωtÞdt

(2)

t¼τ=2

In these expressions ω = 2π/τ is the fundamental frequency and τ is the fundamental period. From one value of n to the next, the harmonic frequency changes by δω = 2π/τ, so the factor preceding the second equation can be replaced by 1/τ = δω/(2π). To avoid confusion when we insert (2) into (1), we change the dummy variable of the integration to u, giving Zτ=2

δω cn ¼ 2π

fðuÞexpðinωuÞdu

(3)

u¼τ=2

After insertion, (1) becomes 0 1 Zτ=2 1 X Bδω C fðtÞ ¼ fðuÞexpðinωuÞduA expðinωtÞ @ 2π n¼1 u¼τ=2

0

¼

1 X

Bδω @ 2π n¼1

Zτ=2

1 C fðuÞexpðinωðt uÞÞduA

(4)

u¼τ=2

We now deﬁne the function within the integral as Zτ=2 fðuÞexpðinωðt uÞÞdu

FðωÞ ¼ u¼τ=2

The initial Fourier series becomes

(5)

56

Mathematical background

fðtÞ ¼

1 1 X FðωÞδω 2π ω¼1

(6)

We now let the incremental frequency δω become very small, tending in the limit to zero; this is equivalent to letting the period τ become inﬁnite. The index n is dropped because ω is now a continuous variable; the discrete sum becomes an integral and the function f(t) is Z1

1 f ðt Þ ¼ 2π

FðωÞdω

(7)

ω¼1

while the function F(ω) from (5) becomes Z1 fðuÞexpðiωðt uÞÞdu

FðωÞ ¼

(8)

u¼1

On inserting F(ω) into (7) we get 2 1 3 Z Z1 1 4 fðtÞ ¼ fðuÞexpðiωðt uÞÞdu5dω 2π ω¼1 u¼1 2 1 3 Z Z1 1 4 fðuÞexpðiωuÞdu5 expðiωtÞdω ¼ 2π ω¼1

(9)

u¼1

The quantity in square brackets, on changing the variable from u back to t, is 1 gðωÞ ¼ 2π

Z1 fðtÞexpðiωtÞdt

(10)

t¼1

and the original expression can now be written Z1 fðtÞ ¼

gðωÞexpðiωtÞdω

(11)

ω¼1

The equivalence of these two equations is known as the Fourier integral theorem.

1.17 Fourier series, integrals, and transforms

57

Fourier series that represent odd or even functions consist of sums of sines or cosines, respectively. In the same way, there are sine and cosine Fourier integrals that represent odd and even functions, respectively. Suppose that the function ƒ(t) is even, and let us replace the complex exponential in (1.272) using (1.5): 1 gð ω Þ ¼ 2π

Z1 fðtÞ½cosðωtÞ i sinðωtÞdt

(1:273)

t¼1

The sine function is odd, so, if ƒ(t) is even, the product ƒ(t)sin(ωt) is odd, and the integral of the second term is zero. The product ƒ(t)cos(ωt) is even, and we can convert the limits of integration to the positive interval: Z1

1 gð ω Þ ¼ 2π ¼

1 π

fðtÞcosðωtÞdt

t¼1 Z1

fðtÞcosðωtÞdt

(1:274)

t¼0

Thus, if ƒ(t) is even, then g(ω) is also even. Similarly, one ﬁnds that, if ƒ(t) is odd, g(ω) is also odd. Now we expand the exponential in (1.271) and apply the same conditions of evenness and oddness to the products: Z1 fðtÞ ¼

gðωÞðcosðωtÞ þ i sinðωtÞÞdω ω¼1 Z1

¼2

gðωÞcosðωtÞdω

(1:275)

ω¼0

If we were to substitute (1.275) back into (1.274), the integration would be preceded by a constant 2/π, the product of the two constants in these equations. Equations (1.275) and (1.274) form a Fourier-transform pair, and it does not matter how the factor 2/π is divided between them. We will associate it here entirely with the second equation, so that we have the pair of equations Z1 gðωÞcosðωtÞdω

fðtÞ ¼ ω¼0

2 gðωÞ ¼ π

(1:276)

Z1 fðtÞcosðωtÞdt

t¼0

58

Mathematical background

The even functions ƒ(t) and g(ω) are Fourier cosine transforms of each other. A similar treatment for a function ƒ(t) that is odd leads to a similar pair of equations in which the Fourier transform g(ω) is also odd and Z1 fðtÞ ¼

gðωÞsinðωtÞdω ω¼0

2 gð ω Þ ¼ π

(1:277)

Z1 fðtÞsinðωtÞdt

t¼0

The odd functions ƒ(t) and g(ω) are Fourier sine transforms of each other. further reading Boas, M. L. (2006). Mathematical Methods in the Physical Sciences, 3rd edn. Hoboken, NJ: Wiley, 839 pp. James, J. F. (2004). A Student’s Guide to Fourier transforms, 2nd edn. Cambridge: Cambridge University Press, 135 pp.

2 Gravitation

2.1 Gravitational acceleration and potential The Universal Law of Gravitation deduced by Isaac Newton in 1687 describes the force of gravitational attraction between two point masses m and M separated by a distance r. Let a spherical coordinate system (r, θ, ) be centered on the point mass M. The force of attraction F exerted on the point mass m acts radially inwards towards M, and can be written F ¼ G

mM er r2

(2:1)

In this expression, G is the gravitational constant (6.674 21 × 10−11 m3 kg−1 s−2), er is the unit radial vector in the direction of increasing r, and the negative sign indicates that the force acts inwardly, towards the attracting mass. The gravitational acceleration aG at distance r is the force on a unit mass at that point: aG ¼ G

M er r2

(2:2)

The acceleration aG may also be written as the negative gradient of a gravitational potential UG aG ¼ rUG

(2:3)

The gravitational acceleration for a point mass is radial, thus the potential gradient is given by

∂UG M ¼ G 2 r ∂r

UG ¼ G 59

M r

(2:4) (2:5)

60

Gravitation

In Newton’s time the gravitational constant could not be veriﬁed in a laboratory experiment. The attraction between heavy masses of suitable dimensions is weak and the effects of friction and air resistance relatively large, so the ﬁrst successful measurement of the gravitational constant by Lord Cavendish was not made until more than a century later, in 1798. However, Newton was able to conﬁrm the validity of the inverse-square law of gravitation in 1687 by using existing astronomic observations of the motions of the planets in the solar system. These had been summarized in three important laws by Johannes Kepler in 1609 and 1619. The small sizes of the planets and the Sun, compared with the immense distances between them, enabled Newton to consider these as point masses and this allowed him to verify the inverse-square law of gravitation.

2.2 Kepler’s laws of planetary motion Johannes Kepler (1571–1630), a German mathematician and scientist, formulated his laws on the basis of detailed observations of planetary positions by Tycho Brahe (1546–1601), a Danish astronomer. The observations were made in the late sixteenth century, without the aid of a telescope. Kepler found that the observations were consistent with the following three laws ( Fig. 2.1). 1. The orbit of each planet is an ellipse with the Sun at one focus. 2. The radius from the Sun to a planet sweeps over equal areas in equal intervals of time.

Q

P* Q*

aphelion

r θ

a

S b

P(r, θ ) perihelion

p

Fig. 2.1. Illustration of Kepler’s laws of planetary motion. The orbit of each planet is an ellipse with the Sun at its focus (S); a, b, and p are the semi-major axis, semiminor axis, and semi-latus rectum, respectively. The area swept by the radius to a planet in a given time is constant (i.e., area SPQ equals area SP*Q*); the square of the period is proportional to the cube of the semi-major axis. After Lowrie ( 2007).

2.2 Kepler’s laws of planetary motion

61

3. The square of the period is proportional to the cube of the semi-major axis of the orbit. The fundamental assumption is that the planets move under the inﬂuence of a central, i.e., radially directed force. For a planet of mass m at distance r from the Sun the force F can be written F¼m

d 2r ¼ fðrÞer dt2

(2:6)

The angular momentum h of the planet about the Sun is h¼rm

dr dt

(2:7)

Differentiating with respect to time, the rate of change of angular momentum is dh d dr dr dr d 2r ¼m r ¼m þm r 2 dt dt dt dt dt dt

(2:8)

The ﬁrst term on the right-hand side is zero, because the vector product of a vector with itself (or with a vector parallel to itself) is zero. Thus dh d 2r ¼rm 2 dt dt

(2:9)

On substituting from ( 2.6) and applying the same condition, we have dh ¼ r fðrÞer ¼ fðrÞðr er Þ ¼ 0 dt

(2:10)

This equation means that h is a constant vector; the angular momentum of the system is conserved. On taking the scalar product of h and r, we obtain dr r·h ¼ r· r m (2:11) dt Rotating the sequence of the vectors in the triple product gives r·h ¼ m

dr · ðr rÞ ¼ 0 dt

(2:12)

This result establishes that the vector r describing the position of a planet is always perpendicular to its constant angular momentum vector h and therefore deﬁnes a plane. Every planetary orbit is therefore a plane that passes through the Sun. The orbit of the Earth deﬁnes the ecliptic plane.

62

Gravitation

2.2.1 Kepler’s Second Law Let the position of a planet in its orbit be described by polar coordinates (r, θ) with respect to the Sun. The coordinates are deﬁned so that the angle θ is zero at the closest approach of the planet to the Sun (perihelion). The angular momentum at an arbitrary point of the orbit has magnitude h ¼ mr2

dθ dt

(2:13)

In a short interval of time Δt the radius vector from the Sun to the planet moves through a small angle Δθ and deﬁnes a small triangle. The area ΔA of the triangle is 1 DA ¼ r2 Dθ 2 The rate of change of the area swept over by the radius vector is dA DA 1 2 Dθ ¼ lim ¼ lim r Dt!0 Dt Dt!0 2 dt Dt

(2:14)

(2:15)

dA 1 2 dθ ¼ r dt 2 dt

(2:16)

dA h ¼ dt 2m

(2:17)

On inserting from ( 2.13) we get

Thus the area swept over by the radius vector in a given time is constant. This is Kepler’s Second Law of planetary motion.

2.2.2 Kepler’s First Law If just the gravitational attraction of the Sun acts on the planet (i.e., we ignore the interactions between the planets), the total energy E of the planet is constant. The total energy E is composed of the planet’s orbital kinetic energy and its potential energy in the Sun’s gravitational ﬁeld: 2 2 1 dr 1 dθ S þ mr2 Gm ¼ E m 2 dt 2 dt r

(2:18)

The ﬁrst term here is the planet’s linear (radial) kinetic energy, the second term is its rotational kinetic energy (with mr2 being the planet’s moment of inertia

2.2 Kepler’s laws of planetary motion

63

about the Sun), and the third term is the gravitational potential energy. On writing dr dr dθ ¼ dt dθ dt

(2:19)

and rearranging terms we get 2 2 2 dr dθ S E 2 dθ þr 2G ¼ 2 dθ dt dt r m

(2:20)

Now, to simplify later steps, we make a change of variables, writing u¼

1 r

(2:21)

Then dr d 1 1 du 2 du ¼ ¼ 2 ¼ r dθ dθ u u dθ dθ Substituting from ( 2.22) into ( 2.20) gives 2 dθ 2 du 2 dθ S E r2 þ r2 2G ¼ 2 dt dθ dt r m

(2:22)

(2:23)

With the result of ( 2.13) we have dθ h ¼ r2 dt m dθ 1 h h r ¼ ¼u dt r m m On replacing these expressions, ( 2.23) becomes 2 2 2 h du E 2 h þu 2uGS ¼ 2 m dθ m m

du dθ

2 þ u2 2uGS

m2 Em ¼2 2 2 h h

(2:24)

(2:25)

(2:26)

The rest of the evaluation is straightforward, if painstaking. First we add a constant to each side,

du dθ

2 þ u2 2uGS

2 2 m2 m2 Em m2 þ GS ¼ 2 þ GS h2 h2 h2 h2

(2:27)

64

Gravitation

du dθ

2

m2 þ u GS 2 h

2

2 Em m2 ¼ 2 2 þ GS 2 h h

(2:28)

Next, we move the second term to the right-hand side of the equation, giving

du dθ

du dθ

2

2 2 Em m2 m2 ¼ 2 2 þ GS 2 u GS 2 h h h

2 ¼

2 2 m2 2Eh2 m2 GS 2 1 þ 2 2 3 u GS 2 h GS m h

(2:29)

(2:30)

Now, we deﬁne some combinations of these terms, as follows: u0 ¼ GS

m2 h2

(2:31)

2Eh2 G2 S2 m3

(2:32)

e2 ¼ 1 þ

Using these deﬁned terms, ( 2.30) simpliﬁes to a more manageable form: 2 du ¼ u20 e2 ðu u0 Þ2 (2:33) dθ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ du ¼ u20 e2 ðu u0 Þ2 (2:34) dθ The solution of this equation, which can be tested by substitution, is u ¼ u0 ð1 þ e cos θÞ

(2:35)

The angle θ is deﬁned to be zero at perihelion. The negative square root in ( 2.34) is chosen because, as θ increases, r increases and u must decrease. Let p¼

1 h2 ¼ u0 GSm2

(2:36)

r¼

p 1 þ e cos θ

(2:37)

This is the polar equation of an ellipse referred to its focus, and is the proof of Kepler’s First Law of planetary motion. The quantity e is the eccentricity of the ellipse, while p is the semi-latus rectum of the ellipse, which is half the length of a chord passing through the focus and parallel to the minor axis ( Fig. 2.1). These equations show that three types of trajectory around the Sun are possible, depending on the value of the total energy E in ( 2.18). If the kinetic

2.2 Kepler’s laws of planetary motion

65

energy is greater than the potential energy, the value of E in ( 2.32) is positive, and e is greater than 1; the path of the object is a hyperbola. If the kinetic energy and potential energy are equal, the total energy is zero and e is exactly 1; the path is a parabola. In each of these two cases the object can escape to inﬁnity, and the paths are called escape trajectories. If the kinetic energy is less than the potential energy, the total energy E is negative and the eccentricity is less than 1. In this case (corresponding to a planet or asteroid) the object follows an elliptical orbit around the Sun.

2.2.3 Kepler’s Third Law It is convenient to describe the elliptical orbit in Cartesian coordinates (x, y), centered on the mid-point of the ellipse, instead of on the Sun. Deﬁne the x-axis parallel to the semi-major axis a of the ellipse and the y-axis parallel to the semiminor axis b. The equation of the ellipse in Fig. 2.1 is x2 y2 þ ¼1 a2 b2

(2:38)

The semi-minor axis is related to the semi-major axis by the eccentricity e, so that b2 ¼ a2 1 e2

(2:39)

The distance of the focus of the ellipse from its center is by deﬁnition ae. The length p of the semi-latus rectum is the value of y for a chord through the focus. On setting y = p and x = ae in ( 2.38), we obtain p2 ðaeÞ2 ¼ 1 2 ¼ 1 e2 2 b a 2 p 2 ¼ a 2 1 e2

(2:40) (2:41)

Now consider the application of Kepler’s Second Law to an entire circuit of the elliptical orbit. The area of the ellipse is πab, and the period of the orbit is T, so dA πab ¼ dt T

(2:42)

h 2πab ¼ m T

(2:43)

Using (2.17),

66

Gravitation

From ( 2.36) and ( 2.43) we get the value of the semi-latus rectum, 1 h 2 1 2πab 2 ¼ p¼ GS m GS T

(2:44)

Substituting from ( 2.41) gives 4π 2 a2 b2 4π 2 a4 a 1 e2 ¼ ¼ 1 e2 GST 2 GST 2

(2:45)

After simplifying, we ﬁnally get T 2 4π 2 ¼ a3 GS

(2:46)

The quantities on the right-hand side are constant, so the square of the period is proportional to the cube of the semi-major axis, which is Kepler’s Third Law.

2.3 Gravitational acceleration and the potential of a solid sphere The gravitational potential and acceleration outside and inside a solid sphere may be calculated from the Poisson and Laplace equations, respectively.

2.3.1 Outside a solid sphere, using Laplace’s equation Outside a solid sphere the gravitational potential UG satisﬁes Laplace’s equation ( Section 1.9). If the density is uniform, the potential does not vary with the polar angle θ or azimuth . Under these conditions, Laplace’s equation in spherical polar coordinates ( 2.67) reduces to ∂ 2 ∂UG ¼0 (2:47) r ∂r ∂r This implies that the bracketed quantity that we are differentiating must be a constant, C, r2

∂UG ¼C ∂r

(2:48)

∂UG C ¼ 2 r ∂r

(2:49)

The gravitational acceleration outside the sphere is therefore

2.3 The potential of a solid sphere

67

∂UG C ¼ 2 er ∂r r

(2:50)

aG ðr4RÞ ¼

At its surface the gravitational acceleration has the value ∂UG C aG ðRÞ ¼ ¼ 2 er ∂r R

(2:51)

The boundary condition at the surface of the sphere is that the accelerations determined outside and inside the sphere must be equal there. We use this to derive the value of the constant C. On comparing ( 2.51) and ( 2.60) we have C ¼ GM

(2:52)

On inserting for C in ( 2.50), the gravitational acceleration outside the sphere is aG ðr4RÞ ¼ G

M er r2

(2:53)

The gravitational potential outside the solid sphere is obtained by integrating ( 2.53) with respect to the radius. This gives UG ðr4RÞ ¼ G

M r

(2:54)

2.3.2 Inside a solid sphere, using Poisson’s equation Inside a solid sphere with radius R and uniform density ρ the gravitational potential UG satisﬁes Poisson’s equation ( Section 1.8). Symmetry again requires the use of spherical polar coordinates, and, because the density is uniform, there is no variation of potential with the polar angle θ or azimuth . Poisson’s equation in spherical polar coordinates reduces to 1 ∂ 2 ∂UG ¼ 4πGρ r ∂r r2 ∂r

(2:55)

On multiplying by r2 and integrating with respect to r, we get ∂ 2 ∂UG ¼ 4πGρr2 r ∂r ∂r r2

∂UG 4 ¼ πGρr3 þ C1 ∂r 3

(2:56) (2:57)

This equation has to be valid at the center of the sphere where r = 0, so the constant C1 = 0 and

68

Gravitation ∂UG 4 ¼ πGρr ∂r 3 ∂UG ¼ aG ðr5RÞ ¼ ∂r

(2:58) 4 πGρr er 3

(2:59)

This shows that the gravitational acceleration inside a homogeneous solid sphere is proportional to the distance from its center. At the surface of the sphere, r = R, and the gravitational acceleration is 4 GM aG ðRÞ ¼ πGρR er ¼ 2 er 3 R

(2:60)

where the mass M of the sphere is 4 M ¼ πR3 ρ 3

(2:61)

To obtain the potential inside the solid sphere, we must integrate ( 2.58). This gives 2 UG ¼ πGρr2 þ C2 3

(2:62)

The constant of integration C2 is obtained by noting that the potential must be continuous at the surface of the sphere. Otherwise a discontinuity would exist and the potential gradient (and force) would be inﬁnite. Equating ( 2.54) and ( 2.62) at r = R gives 2 GM 4 πGρR2 þ C2 ¼ ¼ πGρR2 3 R 3

(2:63)

C2 ¼ 2πGρR2

(2:64)

The gravitational potential inside the uniform solid sphere is therefore given by 2 UG ¼ πGρr2 2πGρR2 3

(2:65)

2 UG ¼ πGρ r2 3R2 3

(2:66)

A schematic graph of the variation of the gravitational potential inside and outside a solid sphere is shown in Fig. 2.2.

2.4 Laplace’s equation

69

r /R 0

0

1 inside sphere

2

3

4

outside sphere

0.5 1

U G (R)

1.5 U G(r)

r=R

2

Fig. 2.2. Variation with radial distance r of the gravitational potential inside and outside a solid sphere of radius R. The potential of the surface of the sphere is UG(R).

2.4 Laplace’s equation in spherical polar coordinates In the above examples the sphere was assumed to have uniform density so that only the radial term in Laplace’s equation had to be solved. This is also the case when density varies only with radius. In the Earth, however, lateral variations of the density distribution occur, and the gravitational potential UG is then a solution of the full Laplace equation 1 ∂ 2 ∂UG 1 ∂ ∂UG 1 ∂2 UG ¼0 þ þ r sin θ ∂r ∂θ r2 ∂r r2 sin θ ∂θ r2 sin2 θ ∂2

(2:67)

This equation is solved using the method of separation of variables. This is a valuable mathematical technique, which allows the variables in a partial differential equation to be separated so that only terms in one variable are on one side of the equation and terms in other variables are on the opposite side. A trial solution for UG is UG ðr; θ; Þ ¼

View more...
william lowrie was born in Hawick, Scotland, and attended the University of Edinburgh, where he graduated in 1960 with ﬁrst-class honors in physics. He achieved a masters degree in geophysics at the University of Toronto and, in 1967, a doctorate at the University of Pittsburgh. After two years in the research laboratory of Gulf Oil Company he became a researcher at the Lamont-Doherty Geological Observatory of Columbia University. In 1974 he was elected professor of geophysics at the ETH Zürich (Swiss Federal Institute of Technology in Zurich), Switzerland, where he taught and researched until retirement in 2004. His research in rock magnetism and paleomagnetism consisted of deducing the Earth’s magnetic ﬁeld in the geological past from the magnetizations of dated rocks. The results were applied to the solution of geologic-tectonic problems, and to analysis of the polarity history of the geomagnetic ﬁeld. Professor Lowrie has authored 135 scientiﬁc articles and a second edition of his acclaimed 1997 textbook Fundamentals of Geophysics was published in 2007. He has been President of the European Union of Geosciences (1987–9) and Section President and Council member of the American Geophysical Union (2000–2). He is a Fellow of the American Geophysical Union and a Member of the Academia Europaea.

A Student’s Guide to Geophysical Equations WILLIAM LOWRIE Institute of Geophysics Swiss Federal Institute of Technology Zurich, Switzerland

cambridge university press Cambridge, New York, Melbourne, Madrid, Cape Town, Singapore, São Paulo, Delhi, Tokyo, Mexico City Cambridge University Press The Edinburgh Building, Cambridge CB2 8RU, UK Published in the United States of America by Cambridge University Press, New York www.cambridge.org Information on this title: www.cambridge.org/9781107005846 © William Lowrie 2011 This publication is in copyright. Subject to statutory exception and to the provisions of relevant collective licensing agreements, no reproduction of any part may take place without the written permission of Cambridge University Press. First published 2011 Printed in the United Kingdom at the University Press, Cambridge A catalogue record for this publication is available from the British Library Library of Congress Cataloguing in Publication data Lowrie, William, 1939– A student’s guide to geophysical equations / William Lowrie. p. cm. Includes bibliographical references and index. ISBN 978-1-107-00584-6 (hardback) 1. Geophysics – Mathematics – Handbooks, manuals, etc. 2. Physics – Formulae – Handbooks, manuals, etc. 3. Earth – Handbooks, manuals, etc. I. Title. QC809.M37L69 2011 550.10 51525–dc22 2011007352 ISBN 978 1 107 00584 6 Hardback ISBN 978 0 521 18377 2 Paperback Cambridge University Press has no responsibility for the persistence or accuracy of URLs for external or third-party internet websites referred to in this publication, and does not guarantee that any content on such websites is, or will remain, accurate or appropriate.

This book is dedicated to Marcia

Contents

Preface Acknowledgments 1

2

page xi xiii

Mathematical background 1.1 Cartesian and spherical coordinates 1.2 Complex numbers 1.3 Vector relationships 1.4 Matrices and tensors 1.5 Conservative force, ﬁeld, and potential 1.6 The divergence theorem (Gauss’s theorem) 1.7 The curl theorem (Stokes’ theorem) 1.8 Poisson’s equation 1.9 Laplace’s equation 1.10 Power series 1.11 Leibniz’s rule 1.12 Legendre polynomials 1.13 The Legendre differential equation 1.14 Rodrigues’ formula 1.15 Associated Legendre polynomials 1.16 Spherical harmonic functions 1.17 Fourier series, Fourier integrals, and Fourier transforms Further reading Gravitation 2.1 Gravitational acceleration and potential 2.2 Kepler’s laws of planetary motion 2.3 Gravitational acceleration and the potential of a solid sphere 2.4 Laplace’s equation in spherical polar coordinates 2.5 MacCullagh’s formula for the gravitational potential Further reading

vii

1 1 1 4 8 17 18 20 23 26 28 32 32 34 41 43 49 52 58 59 59 60 66 69 74 85

viii

3

4

5

6

7

8

Contents

Gravity 3.1 The ellipticity of the Earth’s ﬁgure 3.2 The geopotential 3.3 The equipotential surface of gravity 3.4 Gravity on the reference spheroid 3.5 Geocentric and geographic latitude 3.6 The geoid Further reading The tides 4.1 Origin of the lunar tide-raising forces 4.2 Tidal potential of the Moon 4.3 Love’s numbers and the tidal deformation 4.4 Tidal friction and deceleration of terrestrial and lunar rotations Further reading Earth’s rotation 5.1 Motion in a rotating coordinate system 5.2 The Coriolis and Eötvös effects 5.3 Precession and forced nutation of Earth’s rotation axis 5.4 The free, Eulerian nutation of a rigid Earth 5.5 The Chandler wobble Further reading Earth’s heat 6.1 Energy and entropy 6.2 Thermodynamic potentials and Maxwell’s relations 6.3 The melting-temperature gradient in the core 6.4 The adiabatic temperature gradient in the core 6.5 The Grüneisen parameter 6.6 Heat ﬂow Further reading Geomagnetism 7.1 The dipole magnetic ﬁeld and potential 7.2 Potential of the geomagnetic ﬁeld 7.3 The Earth’s dipole magnetic ﬁeld 7.4 Secular variation 7.5 Power spectrum of the internal ﬁeld 7.6 The origin of the internal ﬁeld Further reading Foundations of seismology 8.1 Elastic deformation

86 86 88 91 96 102 106 115 116 116 119 124 130 136 137 138 140 142 155 157 169 170 171 172 176 178 179 182 197 198 198 200 205 213 214 217 225 227 227

Contents

ix

8.2 Stress 8.3 Strain 8.4 Perfectly elastic stress–strain relationships 8.5 The seismic wave equation 8.6 Solutions of the wave equation 8.7 Three-dimensional propagation of plane P- and S-waves Further reading

228 233 239 244 252 254 258

Appendix A Magnetic poles, the dipole ﬁeld, and current loops Appendix B Maxwell’s equations of electromagnetism References Index

259 265 276 278

Preface

This work was written as a supplementary text to help students understand the mathematical steps in deriving important equations in classical geophysics. It is not intended to be a primary textbook, nor is it intended to be an introduction to modern research in any of the topics it covers. It originated in a set of handouts, a kind of “do-it-yourself ” manual, that accompanied a course I taught on theoretical geophysics. The lecture aids were necessary for two reasons. First, my lectures were given in German and there were no comprehensive up-to-date texts in the language; the recommended texts were in English, so the students frequently needed clariﬁcation. Secondly, it was often necessary to explain classical theory in more detail than one ﬁnds in a multi-topic advanced textbook. To keep such a book as succinct as possible, the intermediate steps in the mathematical derivation of a formula must often be omitted. Sometimes the unassisted student cannot ﬁll in the missing steps without individual tutorial assistance, which is usually in short supply at most universities, especially at large institutions. To help my students in these situations, the “do-it-yourself ” text that accompanied my lectures explained missing details in the derivations. This is the background against which I prepared the present guide to geophysical equations, in the hope that it might be helpful to other students at this level of study. The classes that I taught to senior grades were largely related to potential theory and primarily covered topics other than seismology, since this was the domain of my colleagues and better taught by a true seismologist than by a paleomagnetist! Theoretical seismology is a large topic that merits its own treatment at an advanced level, and there are several textbooks of classical and modern vintage that deal with this. However, a short chapter on the relationship of stress, strain, and the propagation of seismic waves is included here as an introduction to the topic. Computer technology is an essential ingredient of progress in modern geophysics, but a well-trained aspiring geophysicist must be able to do more than xi

xii

Preface

apply advanced software packages. A fundamental mathematical understanding is needed in order to formulate a geophysical problem, and numerical computational skills are needed to solve it. The techniques that enabled scientists to understand much about the Earth in the pre-computer era also underlie much of modern methodology. For this reason, a university training in geophysics still requires the student to work through basic theory. This guide is intended as a companion in that process. Historically, most geophysicists came from the ﬁeld of physics, for which geophysics was an applied science. They generally had a sound training in mathematics. The modern geophysics student is more likely to have begun studies in an Earth science discipline, the mathematical background might be heavily oriented to the use of tailor-made packaged software, and some students may be less able to handle advanced mathematical topics without help or tutoring. To ﬁll these needs, the opening chapter of this book provides a summary of the mathematical background for topics handled in subsequent chapters.

Acknowledgments

In writing this book I have beneﬁted from the help and support of various people. At an early stage, anonymous proposal reviewers gave me useful suggestions, not all of which have been acted on, but all of which were appreciated. Each chapter was read and checked by an obliging colleague. I wish to thank Dave Chapman, Rob Coe, Ramon Egli, Chris Finlay, Valentin Gischig, Klaus Holliger, Edi Kissling, Emile Klingelé, Alexei Kuvshinov, Germán Rubino, Rolf Sidler, and Doug Smylie for their corrections and suggestions for improvement. The responsibility for any errors that escaped scrutiny is, of course, mine. I am very grateful to Derrick Hasterok and Dave Chapman for providing me with an unpublished ﬁgure from Derrick’s Ph.D. thesis. Dr. Susan Francis, Senior Commissioning Editor at Cambridge University Press, gave me constant support and friendly encouragement throughout the many months of writing, for which I am sincerely grateful. Above all, I thank my wife Marcia for her generous tolerance of the intrusion of this project into our retirement activities.

xiii

1 Mathematical background

1.1 Cartesian and spherical coordinates Two systems of orthogonal coordinates are used in this book, sometimes interchangeably. Cartesian coordinates (x, y, z) are used for a system with rectangular geometry, and spherical polar coordinates (r, θ, ) are used for spherical geometry. The relationship between these reference systems is shown in Fig. 1.1(a). The convention used here for spherical geometry is deﬁned as follows: the radial distance from the origin of the coordinates is denoted r, the polar angle θ (geographic equivalent: the co-latitude) lies between the radius and the z-axis (geographic equivalent: Earth’s rotation axis), and the azimuthal angle in the x–y plane is measured from the x-axis (geographic equivalent: longitude). Position on the surface of a sphere (constant r) is described by the two angles θ and . The Cartesian and spherical polar coordinates are linked as illustrated in Fig. 1.1(b) by the relationships x ¼ r sin θ cos y ¼ r sin θ sin z ¼ r cos θ

(1:1)

1.2 Complex numbers The numbers we most commonly use in daily life are real numbers. Some of them are also rational numbers. This means that they can be expressed as the quotient of two integers, with the condition that the denominator of the quotient must not equal zero. When the denominator is 1, the real number is an integer. Thus 4, 4/5, 123/456 are all rational numbers. A real number can also be irrational, which means it cannot be expressed as the quotient of two integers. 1

2

Mathematical background

z

(a)

(b)

z = r cosθ

θ r θ

φ

r

φ

y x

r sinθ

y = r sinθ sinφ

x = r sinθ cosφ

Fig. 1.1. (a) Cartesian and spherical polar reference systems. (b) Relationships between the Cartesian and spherical polar coordinates.

Imaginary axis z = x + iy

+y r sinθ

θ

r

r cosθ +x

Real axis

Fig. 1.2. Representation of a complex number on an Argand diagram.

Familiar examples are π, e (the base of natural logarithms), and some square roots, such as √2, √3, √5, etc. The irrational numbers are real numbers that do not terminate or repeat when expressed as decimals. In certain analyses, such as determining the roots of an equation, it is necessary to ﬁnd the square root of a negative real number, e.g. √(–y2), where y is real. The result is an imaginary number. The negative real number can be written as (–1)y2, and its square root is then √(–1)y. The quantity √(–1) is written i and is known as the imaginary unit, so that √(–y2) becomes ±iy. A complex number comprises a real part and an imaginary part. For example, z = x + iy, in which x and y are both real numbers, is a complex number with a real part x and an imaginary part y. The composition of a complex number can be illustrated graphically with the aid of the complex plane (Fig. 1.2). The real part is plotted on the horizontal axis, and the imaginary part on the vertical axis. The two independent parts are orthogonal on the plot and the complex number z

1.2 Complex numbers

3

is represented by their vector sum, deﬁning a point on the plane. The distance r of the point from the origin is given by r¼

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ x2 þ y2

(1:2)

The line joining the point to the origin makes an angle θ with the real (x-)axis, and so r has real and imaginary components r cos θ and r sin θ, respectively. The complex number z can be written in polar form as z ¼ rðcos θ þ i sin θÞ

(1:3)

It is often useful to write a complex number in the exponential form introduced by Leonhard Euler in the late eighteenth century. To illustrate this we make use of inﬁnite power series; this topic is described in Section 1.10. The exponential function, exp(x), of a variable x can be expressed as a power series as in (1.135). On substituting x = iθ, the power series becomes ðiθÞ2 ðiθÞ3 ðiθÞ4 ðiθÞ5 ðiθÞ6 þ þ þ þ þ 2! 3! 4! 5! 6! ðiθÞ2 ðiθÞ4 ðiθÞ6 ðiθÞ3 ðiθÞ5 ¼1þ þ þ þ ðiθÞ þ þ þ 2! 4! 6! 3! 5! θ2 θ4 θ6 θ3 θ 5 (1:4) ¼ 1 þ þ þ i θ þ þ 2! 4! 6! 3! 5!

expðiθÞ ¼ 1 þ ðiθÞ þ

Comparison with (1.135) shows that the ﬁrst bracketed expression on the right is the power series for cos θ; the second is the power series for sin θ. Therefore expðiθÞ ¼ cos θ þ i sin θ

(1:5)

On inserting (1.5) into (1.3), the complex number z can be written in exponential form as z ¼ r expðiθÞ

(1:6)

The quantity r is the modulus of the complex number and θ is its phase. Conversely, using (1.5) the cosine and sine functions can be deﬁned as the sum or difference of the complex exponentials exp(iθ) and exp(–iθ): expðiθÞ þ expðiθÞ 2 expðiθÞ expðiθÞ sin θ ¼ 2i

cos θ ¼

(1:7)

4

Mathematical background

1.3 Vector relationships A scalar quantity is characterized only by its magnitude; a vector has both magnitude and direction; a unit vector has unit magnitude and the direction of the quantity it represents. In this overview the unit vectors for Cartesian coordinates (x, y, z) are written (ex, ey, ez); unit vectors in spherical polar coordinates (r, θ, ) are denoted (er, eθ, e). The unit vector normal to a surface is simply denoted n.

1.3.1 Scalar and vector products The scalar product of two vectors a and b is deﬁned as the product of their magnitudes and the cosine of the angle α between the vectors: a · b ¼ ab cos α

(1:8)

If the vectors are orthogonal, the cosine of the angle α is zero and a·b ¼ 0

(1:9)

The vector product of two vectors is another vector, whose direction is perpendicular to both vectors, such that a right-handed rule is observed. The magnitude of the vector product is the product of the individual vector magnitudes and the sine of the angle α between the vectors: ja bj ¼ ab sin α

(1:10)

If a and b are parallel, the sine of the angle between them is zero and ab¼0

(1:11)

Applying these rules to the unit vectors (ex, ey, ez), which are normal to each other and have unit magnitude, it follows that their scalar products are ex · ey ¼ ey · ez ¼ ez · ex ¼ 0 e x · ex ¼ e y · ey ¼ e z · ez ¼ 1

(1:12)

The vector products of the unit vectors are ex ey ¼ ez ey ez ¼ ex ez ex ¼ ey ex ex ¼ ey ey ¼ ez ez ¼ 0

(1:13)

1.3 Vector relationships

5

A vector a with components (ax, ay, az) is expressed in terms of the unit vectors (ex, ey, ez) as a ¼ ax ex þ ay ey þ az ez

(1:14)

The scalar product of the vectors a and b is found by applying the relationships in (1.12): a · b ¼ ax ex þ ay ey þ az ez · bx ex þ by ey þ bz ez (1:15) ¼ ax bx þ ay by þ az bz The vector product of the vectors a and b is found by using (1.13): a b ¼ ax ex þ ay ey þ az ez bx ex þ by ey þ bz ez ¼ ay bz az by ex þ ðaz bx ax bz Þey þ ax by ay bx ez

(1:16)

This result leads to a convenient way of evaluating the vector product of two vectors, by writing their components as the elements of a determinant, as follows: ex ey ez a b ¼ ax ay az (1:17) bx by bz The following relationships may be established, in a similar manner to the above, for combinations of scalar and vector products of the vectors a, b, and c: a · ðb cÞ ¼ b · ðc aÞ ¼ c · ða bÞ

(1:18)

a ðb cÞ ¼ bðc · aÞ cða · bÞ

(1:19)

ða bÞ c ¼ bðc · aÞ aðb · cÞ

(1:20)

1.3.2 Vector differential operations The vector differential operator ∇ is deﬁned relative to Cartesian axes (x, y, z) as r ¼ ex

∂ ∂ ∂ þ ey þ ez ∂x ∂y ∂z

(1:21)

The vector operator ∇ determines the gradient of a scalar function, which may be understood as the rate of change of the function in the direction of each of the reference axes. For example, the gradient of the scalar function φ with respect to Cartesian axes is the vector

6

Mathematical background

rφ ¼ ex

∂φ ∂φ ∂φ þ ey þ ez ∂x ∂y ∂z

(1:22)

The vector operator ∇ can operate on either a scalar quantity or a vector. The scalar product of ∇ with a vector is called the divergence of the vector. Applied to the vector a it is equal to ∂ ∂ ∂ · ax ex þ ay ey þ az ez r · a ¼ ex þ ey þ ez ∂x ∂y ∂z ∂ax ∂ay ∂az (1:23) ¼ þ þ ∂x ∂y ∂z If the vector a is deﬁned as the gradient of a scalar potential φ, as in (1.22), we can substitute potential gradients for the vector components (ax, ay, az). This gives ∂ ∂φ ∂ ∂φ ∂ ∂φ r · rφ ¼ þ þ (1:24) ∂x ∂x ∂y ∂y ∂z ∂z By convention the scalar product (∇ · ∇) on the left is written ∇2. The resulting identity is very important in potential theory and is encountered frequently. In Cartesian coordinates it is r2 φ ¼

∂2 φ ∂2 φ ∂2 φ þ þ ∂x2 ∂y2 ∂z2

(1:25)

The vector product of ∇ with a vector is called the curl of the vector. The curl of the vector a may be obtained using a determinant similar to (1.17): ex ey ez r a ¼ ∂=∂x ∂=∂y ∂=∂z (1:26) ax ay az In expanded format, this becomes ∂az ∂ay ∂ax ∂az ∂ay ∂ax ex þ ey þ ez ra¼ ∂y ∂z ∂z ∂x ∂x ∂y

(1:27)

The curl is sometimes called the rotation of a vector, because of its physical interpretation (Box 1.1). Some commonly encountered divergence and curl operations on combinations of the scalar quantity φ and the vectors a and b are listed below: r · ðφaÞ ¼ ðrφÞ · a þ φðr · aÞ

(1:28)

1.3 Vector relationships

7

Box 1.1. The curl of a vector The curl of a vector at a given point is related to the circulation of the vector about that point. This interpretation is best illustrated by an example, in which a ﬂuid is rotating about a point with constant angular velocity ω. At distance r from the point the linear velocity of the ﬂuid v is equal to ω × r. Taking the curl of v, and applying the identity (1.31) with ω constant, r v ¼ r ðw rÞ ¼ wðr · rÞ ðw · rÞr

(1)

To evaluate the ﬁrst term on the right, we use rectangular coordinates (x, y, z): ∂ ∂ ∂ wðr · rÞ ¼ w ex þ ey þ ez · xex þ yey þ zez ∂x ∂y ∂z ¼ w ex · ex þ ey · ey þ ez · ez ¼ 3w (2) The second term is

ðw · rÞr ¼

ωx

∂ ∂ ∂ þ ωy þ ωz · xex þ yey þ zez ∂x ∂y ∂z

¼ ωx ex þ ωy ey þ ωz ez ¼ w

(3)

Combining the results gives r v ¼ 2w

(4)

1 w ¼ ð r vÞ 2

(5)

Because of this relationship between the angular velocity and the linear velocity of a ﬂuid, the curl operation is often interpreted as the rotation of the ﬂuid. When ∇ × v = 0 everywhere, there is no rotation. A vector that satisﬁes this condition is said to be irrotational.

r · ða bÞ ¼ b · ðr aÞ a · ðr bÞ

(1:29)

r ðφaÞ ¼ ðrφÞ a þ φðr aÞ

(1:30)

r ða bÞ ¼ aðr · bÞ bðr · aÞ ða · rÞb þ ðb · rÞa r ðrφÞ ¼ 0

(1:31) (1:32)

8

Mathematical background

m3 z0

n3 z

θ3

θ1

θ2 φ3

φ1

x0 m1

χ3

χ1 φ2

χ2

y

n2

y0 m 2

x n1

Fig. 1.3. Two sets of Cartesian coordinate axes, (x, y, z) and (x0, y0, z0), with corresponding unit vectors (n1, n2, n3) and (m1, m2, m3), rotated relative to each other.

r · ðr a Þ ¼ 0

(1:33)

r ðr aÞ ¼ rðr · aÞ r2 a

(1:34)

It is a worthwhile exercise to establish these identities from basic principles, especially (1.19) and (1.31)–(1.34), which will be used in later chapters.

1.4 Matrices and tensors 1.4.1 The rotation matrix Consider two sets of orthogonal Cartesian coordinate axes (x, y, z) and (x0, y0, z0) that are inclined to each other as in Fig. 1.3. The x0-axis makes angles (1, χ1, θ1) with each of the (x, y, z) axes in turn. Similar sets of angles (2, χ2, θ2) and (3, χ3, θ3) are deﬁned by the orientations of the y0- and z0-axes, respectively, to the (x, y, z) axes. Let the unit vectors along the (x, y, z) and (x0, y0, z0) axes be (n1, n2, n3) and (m1, m2, m3), respectively. The vector r can be expressed in either system, i.e., r = r(x, y, z) = r(x0, y0, z0), or, in terms of the unit vectors, r ¼ xn1 þ yn2 þ zn3 ¼ x0 m1 þ y0 m2 þ z0 m3

(1:35)

We can write the scalar product (r · m1) as r · m1 ¼ xn1 · m1 þ yn2 · m1 þ zn3 · m1 ¼ x0

(1:36)

The scalar product (n1 · m1) = cos 1 = α11 deﬁnes α11 as the direction cosine of the x0-axis with respect to the x-axis (Box 1.2). Similarly, (n2 · m1) = cos χ1 = α12

1.4 Matrices and tensors

9

and (n3 · m1) = cos θ1 = α13 deﬁne α12 and α13 as the direction cosines of the x0-axis with respect to the y- and z-axes, respectively. Thus, (1.36) is equivalent to x0 ¼ α11 x þ α12 y þ α13 z

(1:37)

On treating the y0- and z0-axes in the same way, we get their relationships to the (x, y, z) axes: y0 ¼ α21 x þ α22 y þ α23 z z0 ¼ α31 x þ α32 y þ α33 z

(1:38)

The three equations can be written as a single matrix equation 2 3 2 32 3 2 3 x0 α11 α12 α13 x x 4 y0 5 ¼ 4 α21 α22 α23 54 y 5 ¼ M4 y 5 z z z0 α31 α32 α33

(1:39)

The coefﬁcients αnm (n = 1, 2, 3; m = 1, 2, 3) are the cosines of the interaxial angles. By deﬁnition, α12 = α21, α23 = α32, and α31 = α13, so the square matrix M is symmetric. It transforms the components of the vector in the (x, y, z) coordinate system to corresponding values in the (x0, y0, z0) coordinate system. It is thus equivalent to a rotation of the reference axes. Because of the orthogonality of the reference axes, useful relationships exist between the direction cosines, as shown in Box 1.2. For example, ðα11 Þ2 þ ðα12 Þ2 þ ðα13 Þ2 ¼ cos2 1 þ cos2 χ 1 þ cos2 θ1 ¼

1 2 x þ y2 þ z2 ¼ 1 2 r (1:40)

and α11 α21 þ α12 α22 þ α13 α23 ¼ cos 1 cos 2 þ cos χ 1 cos χ 2 þ cos θ1 cos θ2 ¼ 0 (1:41) The last summation is zero because it is the cosine of the right angle between the x0-axis and the y0-axis. These two results can be summarized as 3 X 1; m¼n αmk αnk ¼ (1:42) 0; m 6¼ n k¼1

1.4.2 Eigenvalues and eigenvectors The transpose of a matrix X with elements αnm is a matrix with elements αmn (i.e., the elements in the rows are interchanged with corresponding elements in

10

Mathematical background

Box 1.2. Direction cosines The vector r is inclined at angles α, β, and γ, respectively, to orthogonal reference axes (x, y, z) with corresponding unit vectors (ex, ey, ez), as in Fig. B1.2. The vector r can be written r ¼ xex þ yey þ zez

(1)

where (x, y, z) are the components of r with respect to these axes. The scalar products of r with ex, ey, and ez are r

z

γ

α

x

β

ez ex

ey

y

Fig. B1.2. Angles α, β, and γ deﬁne the tilt of a vector r relative to orthogonal reference axes (x, y, z), respectively. The unit vectors (ex, ey, ez) deﬁne the coordinate system.

r · ex ¼ x ¼ r cos α r · ey ¼ y ¼ r cos β r · ez ¼ z ¼ r cos γ

(2)

Therefore, the vector r in (1) is equivalent to r ¼ ðr cos αÞex þ ðr cos βÞey þ ðr cos γÞez

(3)

The unit vector u in the direction of r has the same direction as r but its magnitude is unity: r u ¼ ¼ ðcos αÞex þ ðcos βÞey þ ðcos γÞez ¼ lex þ mey þ nez r where (l, m, n) are the cosines of the angles that the vector r makes with the reference axes, and are called the direction cosines of r. They are useful for describing the orientations of lines and vectors.

(4)

1.4 Matrices and tensors

11

The scalar product of two unit vectors is the cosine of the angle they form. Let u1 and u2 be unit vectors representing straight lines with direction cosines (l1, m1, n1) and (l2, m2, n2), respectively, and let θ be the angle between the vectors. The scalar product of the vectors is u1 · u2 ¼ cos θ ¼ l1 ex þ m1 ey þ n1 ez · l2 ex þ m2 ey þ n2 ez (5) Therefore, cos θ ¼ l1 l2 þ m1 m2 þ n1 n2

(6)

The square of a unit vector is the scalar product of the vector with itself and is equal to 1: u·u ¼

r·r ¼1 r2

(7)

On writing the unit vector u as in (4), and applying the orthogonality conditions from (2), we ﬁnd that the sum of the squares of the direction cosines of a line is unity: (8) lex þ mey þ nez · lex þ mey þ nez ¼ l2 þ m2 þ n2 ¼ 1

the columns). The transpose of a (3 × 1) column matrix is a (1 × 3) row matrix. For example, if X is a column matrix given by 2 3 x X ¼ 4y5 (1:43) z then its transpose is the row matrix X T, where XT ¼ ½ x T

The matrix equation X MX = K, where surface: 2 α11 XT MX ¼ ½ x y z 4 α21 α31

y

z

(1:44)

K is a constant, deﬁnes a quadric α12 α22 α32

32 3 x α13 α23 54 y 5 ¼ K z α33

The symmetry of the matrix leads to the equation of this surface:

(1:45)

12

Mathematical background

fðx; y; zÞ ¼ α11 x2 þ α22 y2 þ α33 z2 þ 2α12 xy þ 2α23 yz þ 2α31 zx ¼ K (1:46) When the coefﬁcients αnm are all positive real numbers, the geometric expression of the quadratic equation is an ellipsoid. The normal direction n to the surface of the ellipsoid at the point P(x, y, z) is the gradient of the surface. Using the relationships between (x, y, z) and (x0, y0, z0) in (1.39) and the symmetry of the rotation matrix, αnm = αmn for n ≠ m, the normal direction has components ∂f ¼ 2ðα11 x þ α12 y þ α13 zÞ ¼ 2x0 ∂x ∂f ¼ 2ðα21 x þ α22 y þ α23 zÞ ¼ 2y0 ∂y ∂f ¼ 2ðα31 x þ α32 y þ α33 zÞ ¼ 2z0 ∂z

(1:47)

and we write nðx; y; zÞ ¼ rf ¼ ex

∂f ∂f ∂f þ ey þ e z ∂x ∂x ∂x

(1:48)

nðx; y; zÞ ¼ 2ðx0 ex þ y0 ey þ z0 ez Þ ¼ 2rðx0 ; y0 ; z0 Þ (1:49) The normal n to the surface at P(x, y, z) in the original coordinates is parallel to the vector r at the point (x0, y0, z0) in the rotated coordinates (Fig. 1.4). The transformation matrix M has the effect of rotating the reference axes from one orientation to another. A particular matrix exists that will cause the directions of the (x0, y0, z0) axes to coincide with the (x, y, z) axes. In this case the normal to the surface of the ellipsoid is one of the three principal axes of the

tangent plane

z x

(x0,y0,z0) r

er n P(x,y,z)

y

Fig. 1.4. Location of a point (x, y, z) on an ellipsoid, where the normal n to the surface is parallel to the radius vector at the point (x0, y0, z0).

1.4 Matrices and tensors

13

ellipsoid. The component x0 is then proportional to x, y0 is proportional to y, and z0 is proportional to z. Let the proportionality constant be β. Then x0 = βx, y0 = βy, and z0 = βz, and we get the set of simultaneous equations ðα11 βÞx þ α12 y þ α13 z ¼ 0 α21 x þ ðα22 βÞy þ α23 z ¼ 0

(1:50)

α31 x þ α32 y þ ðα33 βÞz ¼ 0 which, in matrix form, is 2

α11 β 4 α21 α31

α12 α22 β α32

32 3 x α13 α23 54 y 5 ¼ 0 α33 β z

(1:51)

The simultaneous equations have a non-trivial solution only if the determinant of coefﬁcients is zero, i.e., α11 β α21 α31

α12 α22 β α32

α13 α23 ¼ 0 α33 β

(1:52)

This equation is a third-order polynomial in β. Its three roots (β1, β2, β3) are known as the eigenvalues of the matrix M. When each eigenvalue βn is inserted in turn into (1.50) it deﬁnes the components of a corresponding vector vn, which is called an eigenvector of M. Note that (1.51) is equivalent to the matrix equation 2

α11 4 α21 α31

α12 α22 α32

32 3 2 α13 x 1 α23 54 y 5 β4 0 α33 z 0

0 1 0

32 3 0 x 0 54 y 5 ¼ 0 1 z

(1:53)

which we can write in symbolic form ðM βIÞX ¼ 0

(1:54)

The matrix I, with diagonal elements equal to 1 and off-diagonal elements 0, is called a unit matrix: 2

1 0 I ¼ 40 1 0 0

3 0 05 1

(1:55)

14

Mathematical background

1.4.3 Tensor notation Equations describing vector relationships can become cumbersome when written in full or symbolic form. Tensor notation provides a succinct alternative way of writing the equations. Instead of the alphabetic indices used in the previous section, tensor notation uses numerical indices that allow summations to be expressed in a compact form. Let the Cartesian coordinates (x, y, z) be replaced by coordinates (x1, x2, x3) and let the corresponding unit vectors be (e1, e2, e3). The vector a in (1.14) becomes X a ¼ a1 e1 þ a2 e2 þ a3 e3 ¼ ai e i (1:56) i¼1;2;3

A convention introduced by Einstein drops the summation sign and tacitly assumes that repetition of an index implies summation over all values of the index, in this case from 1 to 3. The vector a is then written explicitly a ¼ ai e i

(1:57)

Alternatively, the unit vectors can be implied and the expression ai is understood to represent the vector a. Using the summation convention, (1.15) for the scalar product of two vectors a and b is a · b ¼ a1 b1 þ a2 b2 þ a3 b3 ¼ ai bi

(1:58)

Suppose that two vectors a and b are related, so that each component of a is a linear combination of the components of b. The relationship can be expressed in tensor notation as ai ¼ Tij bj

(1:59)

The indices i and j identify components of the vectors a and b; each index takes each of the values 1, 2, and 3 in turn. The quantity Tij is a second-order (or second-rank) tensor, representing the array of nine coefﬁcients (i.e., 32). A vector has three components (i.e., 31) and is a ﬁrst-order tensor; a scalar property has a single (i.e., 30) value, its magnitude, and is a zeroth-order tensor. To write the cross product of two vectors we need to deﬁne a new quantity, the Levi-Civita permutation tensor εijk. It has the value +1 when a permutation of the indices is even (i.e., ε123 = ε231 = ε312 = 1) and the value –1 when a permutation of the indices is odd (i.e., ε132 = ε213 = ε321 = –1). If any pair of indices is equal, εijk = 0. This enables us to write the cross product of two vectors in tensor notation. Let u be the cross product of vectors a and b: u ¼ a b ¼ ða2 b3 a3 b2 Þe1 þ ða3 b1 a1 b3 Þe2 þ ða1 b2 a2 b1 Þe3 (1:60)

1.4 Matrices and tensors

15

In tensor notation this is written ui ¼ εijk aj bk

(1:61)

This can be veriﬁed readily for each component of u. For example, u1 ¼ ε123 a2 b3 þ ε132 a3 b2 ¼ a2 b3 a3 b2

(1:62)

The tensor equivalent to the unit matrix deﬁned in (1.55) is known as Kronecker’s symbol, δij, or alternatively the Kronecker delta. It has the values 1; if i ¼ j δij ¼ (1:63) 0; if i 6¼ j Kronecker’s symbol is convenient for selecting a particular component of a tensor equation. For example, (1.54) can be written in tensor form using the Kronecker symbol: (1:64) αij βδij xj ¼ 0 This represents the set of simultaneous equations in (1.50). Likewise, the relationship between direction cosines in (1.42) simpliﬁes to αmk αnk ¼ δmn

(1:65)

in which a summation over the repeated index is implied.

1.4.4 Rotation of coordinate axes Let vk be a vector related to the coordinates xl by the tensor Tkl vk ¼ Tkl xl

(1:66)

A second set of coordinates x′n is rotated relative to the axes xl so that the direction cosines of the angles between corresponding axes are the elements of the tensor αnl: x0n ¼ αnl xl

(1:67)

Let the same vector be related to the rotated coordinate axes x′n by the tensor T ′kn: v0k ¼ T 0kn x0n

(1:68)

vk and v′k are the same vector, expressed relative to different sets of axes. Therefore, v0k ¼ αkn vn ¼ αkn Tnl xl

(1:69)

16

Mathematical background

Equating the expressions in (1.68) and (1.69) for v′k gives 0 0 T kn xn ¼ αkn Tnl xl

(1:70)

Using the relationships between the axes in (1.67), 0 0 0 xn ¼ T kn αnl xl T kn

(1:71)

0 αnl ¼ αkn Tnl T kn

(1:72)

Therefore,

On multiplying by αml and summing, 0 αml αnl T kn ¼ αml αkn Tnl

(1:73)

Note that in expanded form the products of direction cosines on the left are equal to αml αnl ¼ αm1 αn1 þ αm2 αn2 þ αm3 αn3 ¼ δmn

(1:74)

as a result of (1.42). Therefore the transformation matrix in the rotated coordinate system is related to the original matrix by the direction cosines between the two sets of axes: 0 ¼ αml αkn Tnl T km

(1:75)

The indices m and k can be interchanged without affecting the result. The sequence of terms in the summation changes, but its sum does not. Therefore, 0 ¼ αkl αmn Tnl T km

(1:76)

This relationship allows us to compute the elements of a matrix in a new coordinate system that is rotated relative to the original reference axes by angles that have the set of direction cosines αnl.

1.4.5 Vector differential operations in tensor notation In tensor notation the vector differential operator ∇ in Cartesian coordinates becomes r ¼ ei

∂ ∂xi

(1:77)

The gradient of a scalar function φ with respect to Cartesian unit vectors (e1, e2, e3) is therefore

1.5 Conservative force, ﬁeld, and potential

rφ ¼ e1

∂φ ∂φ ∂φ ∂φ þ e2 þ e3 ¼ ei ∂x1 ∂x2 ∂x3 ∂xi

17

(1:78)

Several shorthand forms of this equation are in common use; for example, ∂φ ¼ ðrφÞi ¼ φ;i ¼ ∂i φ ∂xi

(1:79)

The divergence of the vector a is written in tensor notation as r·a ¼

∂a1 ∂a2 ∂a3 ∂ai þ þ ¼ ¼ ∂ i ai ∂x1 ∂x2 ∂x3 ∂xi

The curl (or rotation) of the vector a becomes ∂a3 ∂a2 ∂a1 ∂a3 ∂a2 ∂a1 þ e2 þ e3 r a ¼ e1 ∂x2 ∂x3 ∂x3 ∂x1 ∂x1 ∂x2 ðr aÞi ¼ εijk

∂ak ¼ εijk ∂j ak ∂xj

(1:80)

(1:81)

(1:82)

1.5 Conservative force, ﬁeld, and potential If the work done in moving an object from one point to another against a force is independent of the path between the points, the force is said to be conservative. No work is done if the end-point of the motion is the same as the starting point; this condition is called the closed-path test of whether a force is conservative. In a real situation, energy may be lost, for example to heat or friction, but in an ideal case the total energy E is constant. The work dW done against the force F is converted into a gain dEP in the potential energy of the displaced object. The change in the total energy dE is zero: dE ¼ dEP þ dW ¼ 0

(1:83)

The change in potential energy when a force with components (Fx, Fy, Fz) parallel to the respective Cartesian coordinate axes (x, y, z) experiences elementary displacements (dx, dy, dz) is dEP ¼ dW ¼ Fx dx þ Fy dy þ Fz dz (1:84) The value of a physical force may vary in the space around its source. For example, gravitational and electrical forces decrease with distance from a source mass or electrical charge, respectively. The region in which a physical

18

Mathematical background

quantity exerts a force is called its ﬁeld. Its geometry is deﬁned by lines tangential to the force at any point in the region. The term ﬁeld is also used to express the value of the force exerted on a unit of the quantity. For example, the electric ﬁeld of a charge is the force experienced by a unit charge at a given point; the gravitational ﬁeld of a mass is the force acting on a unit of mass; it is therefore equivalent to the acceleration. In a gravitational ﬁeld the force F is proportional to the acceleration a. The Cartesian components of F are therefore (max, may, maz). The gravitational potential U is deﬁned as the potential energy of a unit mass in the gravitational ﬁeld, thus dEP = m dU. After substituting these expressions into (1.84) we get dU ¼ ax dx þ ay dy þ az dz (1:85) The total differential dU can be written in terms of partial differentials as dU ¼

∂U ∂U ∂U dx þ dy þ dz ∂x ∂y ∂z

(1:86)

On equating coefﬁcients of dx, dy, and dz in these equations: ax ¼

∂U ; ∂x

ay ¼

∂U ; ∂y

az ¼

∂U ∂z

(1:87)

These relationships show that the acceleration a is the negative gradient of a scalar potential U: a ¼ rU

(1:88)

Similarly, other conservative ﬁelds (e.g., electric, magnetostatic) can be derived as the gradient of the corresponding scalar potential. According to the vector identity (1.32) the curl of a gradient is always zero; it follows from (1.88) that a conservative force-ﬁeld F satisﬁes the condition rF¼0

(1:89)

1.6 The divergence theorem (Gauss’s theorem) Let n be the unit vector normal to a surface element of area dS. The ﬂux dΦ of a vector F across the surface element dS (Fig. 1.5) is deﬁned to be the scalar product dΦ ¼ F · n dS

(1:90)

If the angle between F and n is θ, the ﬂux across dS is dΦ ¼ F dS cos θ

(1:91)

1.6 The divergence theorem

19

F n

θ dS dS n

Fig. 1.5. The ﬂux of a vector F across a small surface dS, whose normal n is inclined to the vector, is equal to the ﬂux across a surface dSn normal to the vector.

z

dz Fx Fx + dFx

y x

dy x + dx x

Fig. 1.6. Figure for computing the change in the ﬂux of a vector in the x-direction for a small box with edges (dx, dy, dz).

where F is the magnitude of F. Thus the ﬂux of F across the oblique surface dS is equivalent to that across the projection dSn (=dS cos θ) of dS normal to F. Consider the net ﬂux of the vector F through a rectangular box with edges dx, dy, and dz parallel to the x-, y-, and z-axes, respectively (Fig. 1.6). The area dSx of a side normal to the x-axis equals dy dz. The x-component of the vector at x, where it enters the box, is Fx, and at x + dx, where it leaves the box, it is Fx + dFx. The net ﬂux in the x-direction is dΦx ¼ ððFx þ dFx Þ Fx ÞdSx ¼ dFx dy dz

(1:92)

If the distance dx is very small, the change in Fx may be written to ﬁrst order as dFx ¼

∂Fx dx ∂x

(1:93)

20

Mathematical background

The net ﬂux in the x-direction is therefore dΦx ¼

∂Fx ∂Fx dx dy dz ¼ dV ∂x ∂x

(1:94)

where dV is the volume of the small element. Similar results are obtained for the net ﬂux in each of the y- and z-directions. The total ﬂux of F through the rectangular box is the sum of these ﬂows: dΦ ¼ dΦx þ dΦy þ dΦz dΦ ¼

∂Fx ∂Fy ∂Fz þ þ dV ¼ ðr · FÞdV ∂x ∂y ∂z

(1:95) (1:96)

We can equate this expression with the ﬂux deﬁned in (1.90). The ﬂux through a ﬁnite volume V with a bounding surface of area S and outward normal unit vector n is ZZZ ZZ ðr · FÞdV ¼ F · n dS (1:97) V

S

This is known as the divergence theorem, or Gauss’s theorem, after the German mathematician Carl Friedrich Gauss (1777–1855). Note that the surface S in Gauss’s theorem is a closed surface, i.e., it encloses the volume V. If the ﬂux of F entering the bounding surface is the same as the ﬂux leaving it, the total ﬂux is zero, and so r·F ¼ 0

(1:98)

This is sometimes called the continuity condition because it implies that ﬂux is neither created nor destroyed (i.e., there are neither sources nor sinks of the vector) within the volume. The vector is said to be solenoidal.

1.7 The curl theorem (Stokes’ theorem) Stokes’ theorem relates the surface integral of the curl of a vector to the circulation of the vector around a closed path bounding the surface. Let the vector F pass through a surface S which is divided into a grid of small elements (Fig. 1.7). The area of a typical surface element is dS and the unit vector n normal to the element speciﬁes its orientation. First, we evaluate the work done by F around one of the small grid elements, ABCD (Fig. 1.8). Along each segment of the path we need to consider only the

1.7 The curl theorem

21

n S C

D dS

B

A

C Fig. 1.7. Conﬁguration for Stokes’ theorem: the surface S is divided into a grid of elementary areas dS and is bounded by a closed circuit C.

y

(x + dx, y + dy)

(x, y + dy) D

C

dy

Fy

F

A

(x, y)

B

dx

(x + dx, y)

x

Fx

Fig. 1.8. Geometry for calculation of the work done by a force F around a small rectangular grid.

vector component parallel to that segment. The value of F may vary with position, so, for example, the x-component along AB may differ from the x-component along CD. Provided that dx and dy are inﬁnitesimally small, we can use Taylor series approximations for the components of F (Section 1.10.2). To ﬁrst order we get ∂Fx dy ∂y ∂Fy dx Fy BC ¼ Fy DA þ ∂x ðFx ÞCD ¼ ðFx ÞAB þ

(1:99)

The work done in a circuit around the small element ABCD is the sum of the work done along each individual segment: xZþdx

I F·dl ¼ ABCD

yþdy Z

Fy BC dy þ

ðFx ÞAB dx þ x

y

Zx

Zy ðFx ÞCD dx þ

xþdx

Fy

DA

dy

yþdy

(1:100)

22

Mathematical background

I

xþdx Z

F·dl ¼

ðFx ÞAB ðFx ÞCD dx þ

x

ABCD

yþdy Z

Fy BC Fy DA dy

y

(1:101) Substituting from (1.99) gives xþdx Z

I

yþdy Z ∂Fx ∂Fy dy dx þ dx dy ∂y ∂x

F·dl ¼ ABCD

x

(1:102)

y

The mean-value theorem allows us to replace the integrands over the tiny distances dx and dy by their values at some point in the range of integration: I ∂Fy ∂Fx dx dy (1:103) F·dl ¼ ∂x ∂y ABCD

The bracketed expression is the z-component of the curl of F I F · d l ¼ ðr FÞz dx dy

(1:104)

ABCD

The normal direction n to the small area dS = dx dy is parallel to the z-axis (i.e., out of the plane of Fig. 1.8), and hence is in the direction of (∇ × F)z. Thus, I F · d l ¼ ðr FÞ·n dS (1:105) ABCD

The circuit ABCD is one of many similar grid elements of the surface S. When adjacent elements are compared, the line integrals along their common boundary are equal and opposite. If the integration is carried out for the entire surface S, the only surviving parts are the integrations along the bounding curve C (Fig. 1.7). Thus ZZ I ðr FÞ · n dS ¼ F · d l (1:106) S

C

This equation is known as Stokes’ theorem, after the English mathematician George Gabriel Stokes (1819–1903). It enables conversion of the surface integral of a vector to a line integral. The integration on the left is made over the surface S through which the vector F passes. The closed integration on the right is made around the bounding curve C to the surface S; d l is an inﬁnitesimal element of this

1.8 Poisson’s equation

23

boundary. The direction of dl around the curve is right-handed with respect to the surface S, i.e., positive when the path is kept to the right of the surface, as in Fig. 1.7. Note that the surface S in Stokes’ theorem is an open surface; it is like the surface of a bowl with the bounding curve C as its rim. The integration of F around the rim is called the circulation of F about the curve C. If the integral is zero, there is no circulation and the vector F is said to be irrotational. Comparison with the left-hand side shows that the condition for this is rF¼0

(1:107)

As shown in Section 1.5, this is also the condition for F to be a conservative ﬁeld.

1.8 Poisson’s equation The derivations in this and the following sections are applicable to any ﬁeld that varies as the inverse square of distance from its source. Gravitational acceleration is used as an example, but the electric ﬁeld of a charge may be treated in the same way. Let S be a surface enclosing an observer at P and a point mass m. Let dS be a small element of the surface at distance r in the direction er from the mass m, as in Fig. 1.9. The orientation of dS is speciﬁed by the direction n normal to the surface element. With G representing the gravitational constant (see Section 2.1), the gravitational acceleration aG at dS is given by aG ¼ G

m er r2

(1:108)

Let θ be the angle between the radius and the direction n normal to the surface element, and let the projection of dS normal to the radius be dSn. The solid angle dΩ with apex at the mass is deﬁned as the ratio of the normal surface element dSn to the square of its distance r from the mass (Box 1.3):

S dS m

dΩ

aG

θ

r

er n

dSn

Fig. 1.9. Representation of the ﬂux of the gravitational acceleration aG through a closed surface S surrounding the source of the ﬂux (the point mass m).

24

Mathematical background

Box 1.3. Deﬁnition of a solid angle A small element of the surface of a sphere subtends a cone with apex at the center of the sphere (Fig. B1.3(a)). The solid angle Ω is deﬁned as the ratio of the area A of the surface element to the square of the radius r of the sphere: (a)

r dθ

(b) dθ

dA r

Ω

r dA

θ

A

r sinθ dφ

dφ

r sinθ

Fig. B1.3. (a) Relationship of the solid angle Ω, the area A of an element subtended on the surface of a sphere, and the radius r of the sphere. (b) The surface of a sphere divided into rings, and each ring into small surface elements with sides r dθ and r sin θ d.

Ω¼

A r2

(1)

This deﬁnition can be used for an arbitrarily shaped surface. If the surface is inclined to the radial direction it must be projected onto a surface normal to the radius, as in Fig. 1.5. For example, if the normal to the surface A makes an angle α with the direction from the apex of the subtended cone, the projected area is A cos α and the solid angle Ω subtended by the area is Ω¼

A cos α r2

(2)

As an example, let the area on the surface of a sphere be enclosed by a small circle (Fig. B1.3(b)). Symmetry requires spherical polar coordinates to describe the area within the circle. Let the circle be divided into concentric rings, and let the half-angle subtended by a ring at the center of the sphere be θ. The radius of the ring is r sin θ and its width is r dθ. Let the angular position of a small surface element of the ring be ; the length of a side of the element is then r sin θ d. The area dA of a surface area element is equal to r2 sin θ dθ d. The solid angle subtended at the center of the sphere by the element of area dA is dΩ ¼

r2 sin θ dθ d ¼ sin θ dθ d r2

(3)

1.8 Poisson’s equation

25

This expression is also equivalent to the element of area on the surface of a unit sphere (one with unit radius). Integrating in the ranges 0 ≤ ≤ 2π and 0 ≤ θ ≤ θ0, we get the solid angle Ω0 subtended by a circular region of the surface of a sphere deﬁned by a half-apex angle θ0: Zθ0 Z2π Ω0 ¼

sin θ dθ d ¼ 2π ð1 cos θ0 Þ

(4)

θ¼0 ¼0

The unit of measurement of solid angle is the steradian, which is analogous to the radian in plane geometry. The maximum value of a solid angle is when the surface area is that of the complete sphere, namely 4πr2. The solid angle at its center then has the maximum possible value of 4π. This result is also obtained by letting the half-apex angle θ0 in (4) increase to its maximum value π.

dΩ ¼

dSn dS cos θ ðer · nÞdS ¼ ¼ r2 r2 r2

(1:109)

The ﬂux dN of the gravitational acceleration aG through the area element is m ðer · nÞdS r2

(1:110)

cos θ dS ¼ Gm dΩ r2

(1:111)

dN ¼ aG · n dS ¼ G

dN ¼ Gm

If we integrate this expression over the entire surface S we get the total gravitational ﬂux N, ZZ Z N¼ aG · n dS ¼ Gm dΩ ¼ 4πGm (1:112) Ω

S

Now we replace this surface integral by a volume integration, using the divergence theorem (Section 1.6) ZZZ ZZ ðr · aG ÞdV ¼ aG · n dS ¼ 4πGm (1:113) V

S

This is valid for any point mass m inside the surface S. If the surface encloses many point masses we may replace m with the sum of the point masses. If mass

26

Mathematical background

is distributed in the volume with mean density ρ, a volume integral can replace the enclosed mass: ZZZ ZZZ ðr · aG ÞdV ¼ 4πG ρ dV (1:114) V

V

ZZZ ðr · aG þ 4πGρÞdV ¼ 0

(1:115)

V

For this to be generally true the integrand must be zero. Consequently, r · aG ¼ 4πGρ

(1:116)

The gravitational acceleration is the gradient of the gravitational potential UG as in (1.88): r · ðrUG Þ ¼ 4πGρ r2 UG ¼ 4πGρ

(1:117) (1:118)

Equation (1.118) is known as Poisson’s equation, after Siméon-Denis Poisson (1781–1840), a French mathematician and physicist. It describes the gravitational potential of a mass distribution at a point that is within the mass distribution. For example, it may be used to compute the gravitational potential at a point inside the Earth.

1.9 Laplace’s equation Another interesting case is the potential at a point outside the mass distribution. Let S be a closed surface outside the point mass m. The radius vector r from the point mass m now intersects the surface S at two points A and B, where it forms angles θ1 and θ2 with the respective unit vectors n1 and n2 normal to the surface (Fig. 1.10). Let er be a unit vector in the radial direction. Note that the outward normal n1 forms an obtuse angle with the radius vector at A. The gravitational acceleration at A is a1 and its ﬂux through the surface area dS1 is Gm dN1 ¼ a1 · n1 dS1 ¼ 2 ðr · n1 ÞdS1 (1:119) r1 Gm cos θ1 dS1 (1:120) dN1 ¼ 2 cosðπ θ1 ÞdS1 ¼ Gm r1 r21

1.9 Laplace’s equation

27

S

B m

A

a1

dΩ

θ2

a2

dS1 n1

dS2

θ1

er n2

Fig. 1.10. Representation of the gravitational ﬂux through a closed surface S that does not enclose the source of the ﬂux (the point mass m).

dN1 ¼ Gm dΩ

(1:121)

The gravitational acceleration at B is a2 and its ﬂux through the surface area dS2 is dN2 ¼ a2 · n2 dS2 ¼ Gm

cos θ2 dS2 r22

dN2 ¼ Gm dΩ

(1:122) (1:123)

The total contribution of both surfaces to the gravitational ﬂux is dN ¼ dN1 þ dN2 ¼ 0

(1:124)

Thus, the total ﬂux of the gravitational acceleration aG through a surface S that does not include the point mass m is zero. By invoking the divergence theorem we have for this situation Z Z ðr · aG ÞdV ¼ aG · n dS ¼ 0 (1:125) V

S

For this result to be valid for any volume, the integrand must be zero: r · aG ¼ r · ðrUG Þ ¼ 0

(1:126)

r2 UG ¼ 0

(1:127)

Equation (1.127) is Laplace’s equation, named after Pierre Simon, Marquis de Laplace (1749–1827), a French mathematician and physicist. It describes the gravitational potential at a point outside a mass distribution. For example, it is applicable to the computation of the gravitational potential of the Earth at an external point or on its surface. In Cartesian coordinates, which are rectilinear, Laplace’s equation has the simple form

28

Mathematical background

∂2 UG ∂2 UG ∂2 UG þ þ ¼0 ∂x2 ∂y2 ∂z2

(1:128)

Spherical polar coordinates are curvilinear and the curvature of the angular coordinates results in a more complicated form: 1 ∂ 2 ∂UG 1 ∂ ∂UG 1 ∂2 UG þ þ ¼ 0 (1:129) r sin θ ∂r ∂θ r2 ∂r r2 sin θ ∂θ r2 sin2 θ ∂2

1.10 Power series A function ƒ(x) that is continuous and has continuous derivatives may be approximated by the sum of an inﬁnite series of powers of x. Many mathematical functions – e.g., sin x, cos x, exp(x), ln(1 + x) – fulﬁll these conditions of continuity and can be expressed as power series. This often facilitates the calculation of a value of the function. Three types of power series will be considered here: the MacLaurin, Taylor, and binomial series.

1.10.1 MacLaurin series Let the function ƒ(x) be written as an inﬁnite sum of powers of x: fðxÞ ¼ a0 þ a1 x þ a2 x2 þ a3 x3 þ a4 x4 þ þ an xn þ

(1:130)

The coefﬁcients an in this sum are constants. Differentiating (1.130) repeatedly with respect to x gives df ¼ a1 þ 2a2 x þ 3a3 x2 þ 4a4 x3 þ þ nan xn1 þ dx d 2f ¼ 2a2 þ ð3 2Þa3 x þ ð4 3Þa4 x2 þ þ nðn 1Þan xn2 þ dx2 d 3f ¼ ð3 2Þa3 þ ð4 3 2Þa4 x þ þ nðn 1Þðn 2Þan xn3 þ dx3 (1:131) After n differentiations, the expression becomes d nf ¼ ðnðn 1Þðn 2Þ . . . 3 2 1Þ an þ terms containing powers of x dxn (1:132) Now we evaluate each of the differentiations at x = 0. Terms containing powers of x are zero and

1.10 Power series df fð0Þ ¼ a0 ; ¼ a1 dx x¼0 2 3 d f d f ¼ 2a ; ¼ ð3 2Þa3 2 dx2 x¼0 dx3 x¼0 n d f ¼ ðnðn 1Þðn 2Þ . . . 3 2 1Þan ¼ n!an dxn x¼0

29

(1:133)

On inserting these values for the coefﬁcients into (1.130) we get the power series for ƒ(x): df x2 d 2 f x3 d 3 f þ þ þ fðxÞ ¼ fð0Þ þ x dx x¼0 2! dx2 x¼0 3! dx3 x¼0 (1:134) xn d n f þ þ n! dxn x¼0 This is the MacLaurin series for ƒ(x) about the origin, x = 0. It was derived in the eighteenth century by the Scottish mathematician Colin MacLaurin (1698– 1746) as a special case of a Taylor series. The MacLaurin series is a convenient way to derive series expressions for several important functions. In particular, x3 x5 x7 x2n1 þ þ ð1Þn1 3! 5! 7! ð2n 1Þ! x2 x 4 x6 x2n2 cos x ¼ 1 þ þ ð1Þn1 2! 4! 6! ð2n 2Þ! (1:135) x2 x3 xn1 x expðxÞ ¼ e ¼ 1 þ x þ þ þ þ þ 2! 3! ðn 1Þ! x2 x3 x4 xn lnð1 þ xÞ ¼ loge ð1 þ xÞ ¼ x þ þ ð1Þn1 2 3 4 n sin x ¼ x

1.10.2 Taylor series We can write the power series in (1.134) for ƒ(x) centered on any new origin, for example x = x0. To do this we substitute (x – x0) for x in the above derivation. The power series becomes df ðx x0 Þ2 d 2 f fðxÞ ¼fðx0 Þ þ ðx x0 Þ þ 2! dx x¼x0 dx2 x¼x0 ðx x0 Þ3 d 3 f ðx x0 Þn d n f þ þ þ þ dx3 x¼x0 dxn x¼x0 3! n! (1:136)

30

Mathematical background

This is called a Taylor series, after an English mathematician, Brooks Taylor (1685–1731), who described its properties in 1712. The MacLaurin and Taylor series are both approximations to the function ƒ(x). The remainder between the true function and its power series is a measure of how well the function is expressed by the series.

1.10.3 Binomial series Finite series An important series is the expansion of the function ƒ(x) = (a + x)n. If n is a positive integer, the expansion of ƒ(x) is a ﬁnite series, terminating after (n + 1) terms. Evaluating the series for some low values of n gives the following: n ¼ 0 : ða þ xÞ0 ¼ 1 n ¼ 1 : ða þ xÞ1 ¼ a þ x n ¼ 2 : ða þ xÞ2 ¼ a2 þ 2ax þ x2

(1:137)

n ¼ 3 : ða þ xÞ3 ¼ a3 þ 3a2 x þ 3ax2 þ x3 n ¼ 4 : ða þ xÞ4 ¼ a4 þ 4a3 x þ 6a2 x2 þ 4ax3 þ x4 The general expansion of ƒ(x) is therefore nðn 1Þ n2 2 a x þ 12 nðn 1Þ . . . ðn k þ 1Þ nk k a x þ xn þ k!

ða þ xÞn ¼ an þ nan1 x þ

(1:138)

The coefﬁcient of the general kth term is equivalent to nðn 1Þ . . . ðn k þ 1Þ n! ¼ k! k!ðn kÞ!

(1:139)

This is called the binomial coefﬁcient. When the constant a is equal to 1 and n is a positive integer, we have the useful series expansion ð1 þ xÞn ¼

n X

n! xk k! ð n k Þ! k¼0

(1:140)

Inﬁnite series If the exponent in (1.140) is not a positive integer, the series does not terminate, but is an inﬁnite series. The series for ƒ(x) = (1 + x)p, in which the exponent p is not a positive integer, may be derived as a MacLaurin series:

1.10 Power series

df dx

¼ pð1 þ xÞp1

x¼0

x¼0

31

¼p

d 2f ¼ pðp 1Þð1 þ xÞp2 ¼ pðp 1Þ 2 x¼0 dx x¼0 n d f ¼ ðpðp1Þ . . . ðpnþ1Þð1þxÞpn Þx¼0 ¼ pðp 1Þ . . . ðp n þ 1Þ dxn x¼0 (1:141) On inserting these terms into (1.134), and noting that ƒ(0) = 1, we get for the binomial series pðp 1Þ 2 pðp 1Þðp 2Þ 3 x þ x þ 12 123 pðp 1Þ . . . ðp n þ 1Þ n x þ þ n!

ð1 þ xÞp ¼ 1 þ px þ

(1:142)

If the exponent p is not an integer, or p is negative, the series is convergent in the range –1 < x < 1.

1.10.4 Linear approximations The variations in some physical properties over the surface of the Earth are small in relation to the main property. For example, the difference between the polar radius c and the equatorial radius a expressed as a fraction of the equatorial radius deﬁnes the ﬂattening ƒ, which is equal to 1/298. This results from deformation of the Earth by the centrifugal force of its own rotation, which, expressed as a fraction m of the gravitational force, is equal to 1/289. Both ƒ and m are less than three thousandths of the main property, so ƒ2, m2, and the product fm are of the order of nine parts in a million and, along with higher-order combinations, are negligible. Curtailing the expansion of small quantities at ﬁrst order helps keep equations manageable without signiﬁcant loss of geophysical information. In the following chapters much use will be made of such linear approximation. It simpliﬁes the form of mathematical functions and the usable part of the series described above. For example, for small values of x or (x – x0), the following ﬁrst-order approximations may be used: sin x x; expðxÞ 1 þ x; ð1 þ xÞp 1 þ px

cos x 1 lnð1 þ xÞ x

fðxÞ fðx0 Þ þ ðx x0 Þ

df dx

(1:143)

x¼x0

32

Mathematical background

1.11 Leibniz’s rule Assume that u(x) and v(x) are differentiable functions of x. The derivative of their product is d dvðxÞ duðxÞ ðuðxÞvðxÞÞ ¼ uðxÞ þ vðxÞ dx dx dx

(1:144)

If we deﬁne the operator D = d/dx, we obtain a shorthand form of this equation: DðuvÞ ¼ u Dv þ v Du

(1:145)

We can differentiate the product (uv) a second time by parts, D2 ðuvÞ ¼ DðDðuvÞÞ ¼ u D2 v þ ðDuÞðDvÞ þ ðDvÞðDuÞ þ v D2 u ¼ u D2 v þ 2ðDuÞðDvÞ þ v D2 u

(1:146)

and, continuing in this way, D3 ðuvÞ ¼ u D3 v þ 3ðDuÞ D2 v þ 3 D2 u ðDvÞ þ v D3 u D4 ðuvÞ ¼ u D4 v þ 4ðDuÞ D3 v þ 6 D2 u D2 v þ 4 D2 u ðDvÞ þ v D4 u (1:147) The coefﬁcients in these equations are the binomial coefﬁcients, as deﬁned in (1.139). Thus after n differentiations we have Dn ðuvÞ ¼

n X

k nk n! D u D v k!ðn kÞ! k¼0

(1:148)

This relationship is known as Leibniz’s rule, after Gottfried Wilhelm Leibniz (1646–1716), who invented inﬁnitesimal calculus contemporaneously with Isaac Newton (1642–1727); each evidently did so independently of the other.

1.12 Legendre polynomials Let r and R be the sides of a triangle that enclose an angle θ and let u be the side opposite this angle (Fig. 1.11). The angle and sides are related by the cosine rule u2 ¼ r2 þ R2 2rR cos θ

(1:149)

Inverting this expression and taking the square root gives " ! 1 1 r ¼ 12 cos θ þ u R R

r R

!2 #1=2 (1:150)

1.12 Legendre polynomials

R

33

u

θ r

Fig. 1.11. Relationship of the sides r and R, which enclose an angle θ, and the side u opposite the angle, as used in the deﬁnition of Legendre polynomials.

Now let h ¼ r=R and x = cos θ, giving 1=2 1 1 1 ¼ 1 2xh þ h2 ¼ ð1 tÞ1=2 u R R

(1:151)

where t = 2xh – h2. The equation can be expanded as a binomial series 1 3 1 3 5 2 2 2 2 2 1 1=2 2 ð1tÞ ¼ 1þ ðtÞ þ ðtÞ þ ðtÞ3 þ 2 12 123 !2 !3 t 13 t 135 t ¼1þ þ þ þ 2 12 2 123 2 !n 1 3 5 . . . ð2n 1Þ t þ (1:152) þ 1 2 3...n 2 The inﬁnite series of terms on the right-hand side of the equation can be written ð1 tÞ1=2 ¼

1 X

an t n

(1:153)

n¼0

The coefﬁcient an is given by an ¼

1 3 5 . . . ð2n 1Þ 2n n!

(1:154)

Now, substitute the original expression for t,

1 2xh þ h2

1=2

¼

1 X n¼0

1 n X an 2xh h2 ¼ an hn ð2x hÞn

(1:155)

n¼0

This equation is an inﬁnite series in powers of h. The coefﬁcient of each term in the power series is a polynomial in x. Let the coefﬁcient of hn be Pn(x). The equation becomes 1 1=2 X Ψðx; hÞ ¼ 1 2xh þ h2 ¼ hn Pn ðxÞ n¼0

(1:156)

34

Mathematical background

Equation (1.156) is known as the generating function for the polynomials Pn(x). Using this result, and substituting h = r/R and x = cos θ, we ﬁnd that (1.151) becomes !n 1 1 1X r Pn ðcos θÞ (1:157) ¼ u R n¼0 R The polynomials Pn(x) or Pn(cos θ) are called Legendre polynomials, after the French mathematician Adrien-Marie Legendre (1752–1833). The deﬁning equation (1.157) is called the reciprocal-distance formula. An alternative formulation is given in Box 1.4.

1.13 The Legendre differential equation The Legendre polynomials satisfy an important second-order partial differential equation, which is called the Legendre differential equation. To derive this equation we will carry out a sequence of differentiations, starting with the generating function in the form 1=2 Ψ ¼ 1 2xh þ h2

(1:158)

Differentiating this function once with respect to h gives 3=2 ∂Ψ ¼ ðx hÞ 1 2xh þ h2 ¼ ðx hÞΨ3 ∂h

(1:159)

Differentiating Ψ twice with respect to x gives 3=2 ∂Ψ ¼ hΨ3 ¼ h 1 2xh þ h2 ∂x 1 ∂Ψ Ψ3 ¼ h ∂x ∂2 Ψ ∂Ψ ¼ 3h2 Ψ5 ¼ 3hΨ2 2 ∂x ∂x 1 ∂2 Ψ Ψ5 ¼ 2 2 3h ∂x

(1:160)

(1:161)

Next we perform successive differentiations of the product (hΨ) with respect to h. The ﬁrst gives ∂ ∂Ψ ðhΨÞ ¼ Ψ þ h ¼ Ψ þ hðx hÞΨ3 ∂h ∂h

(1:162)

1.13 The Legendre differential equation

35

Box 1.4. Alternative form of the reciprocal-distance formula The sides and enclosed angle of the triangle in Fig. 1.11 are related by the cosine rule u2 ¼ r2 þ R2 2rR cos θ

(1)

Instead of taking R outside the brackets as in (1.150), we can move r outside and write the expression for u as " 2 #1=2 1 1 R R (2) ¼ 12 cos θ þ u r r r Following the same treatment as in Section 1.12, but now with h ¼ R=r and x = cos θ, we get 1=2 1 1 1 ¼ ð1 tÞ1=2 ¼ 1 2xh þ h2 u r r

(3)

where t = 2xh – h2. The function (1 – t)−1/2 is expanded as a binomial series, which again gives an inﬁnite series in h, in which the coefﬁcient of hn is Pn(x). The deﬁning equation is as before: 1 1=2 X ¼ hn Pn ðxÞ Ψðx; hÞ ¼ 1 2xh þ h2

(4)

n¼0

On substituting h ¼ R=r and x = cos θ, we ﬁnd an alternative form for the generating equation for the Legendre polynomials: 1 n 1 1X R Pn ðcos θÞ (5) ¼ u r n¼0 r

On repeating the differentiation, and taking (1.159) into account, we get ∂2 ∂Ψ ∂Ψ ðhΨÞ ¼ þ hðx hÞ3Ψ2 þ ðx 2hÞΨ3 ∂h ∂h ∂h2 ¼ ðx hÞΨ3 þ 3hðx hÞ2 Ψ5 þ ðx 2hÞΨ3 ¼ ð2x 3hÞΨ3 þ 3hðx hÞ2 Ψ5 Now substitute for Ψ3, from (1.160), and Ψ5, from (1.161), giving

(1:163)

36

Mathematical background ∂2 1 ∂Ψ 1 ∂2 Ψ 2 ð hΨ Þ ¼ ð 2x 3h Þ þ 3h ð x h Þ ∂h2 h ∂x 3h2 ∂x2

(1:164)

Multiply throughout by h: 2 ∂2 ∂Ψ 2 ∂ Ψ h 2 ðhΨÞ ¼ ð2x 3hÞ þ ðx hÞ ∂h ∂x ∂x2 2 ∂Ψ ∂Ψ ∂ Ψ ¼ 2x 3h þ ð x hÞ 2 ∂x ∂x ∂x2

(1:165)

The second term on the right can be replaced as follows, again using (1.160) and (1.161): 3h

∂Ψ 1 ∂2 Ψ ¼ 3h2 Ψ3 ¼ 2 2 ∂x Ψ ∂x ∂2 Ψ ¼ 1 2xh þ h2 ∂x2

On substituting into (1.165) and gathering terms, we get h i ∂2 Ψ ∂2 ∂Ψ þ 2x h 2 ðhΨÞ ¼ ðx hÞ2 1 2xh þ h2 ∂x2 ∂x ∂h ∂2 Ψ ∂2 ∂Ψ h 2 ðhΨÞ ¼ x2 1 þ 2x ∂h ∂x2 ∂x

(1:166)

(1:167)

(1:168)

The Legendre polynomials Pn(x) are deﬁned in (1.156) as the coefﬁcients of hn in the expansion of Ψ as a power series. On multiplying both sides of (1.156) by h, we get hΨ ¼

1 X

hnþ1 Pn ðxÞ

(1:169)

n¼0

We differentiate this expression twice and multiply by h to get a result that can be inserted on the left-hand side of (1.168):

h

1 X ∂ ðn þ 1Þhn Pn ðxÞ ðhΨÞ ¼ ∂h n¼0

(1:170)

1 X ∂2 ðhΨÞ ¼ nðn þ 1Þhn Pn ðxÞ 2 ∂h n¼0

(1:171)

Using (1.156), we can now eliminate Ψ and convert (1.168) into a second-order differential equation involving the Legendre polynomials Pn(x),

1.13 The Legendre differential equation

37

Table 1.1. Some ordinary Legendre polynomials of low degree n

Pn(x)

Pn(cos θ)

0 1

1 x 1 2 3x 1 2 1 3 5x 3x 2 1 35x4 30x2 þ 3 8

1 cos θ 1 3 cos2 θ 1 2 1 5 cos3 θ 3 cos θ 2 1 35 cos4 θ 30 cos2 θ þ 3 8

2 3 4

1 X

hn

n¼0 1 X n¼0

X 1 2 d 2 Pn ðxÞ dPn ðxÞ x 1 þ 2x nðn þ 1Þhn Pn ðxÞ (1:172) ¼ dx2 dx n¼0

hn

2 d 2 Pn ðxÞ dPn ðxÞ x 1 þ 2x ðxÞ ¼0 nðn þ 1ÞPn dx2 dx

(1:173)

If this expression is true for every non-zero value of h, the quantity in curly brackets must be zero, thus d 2 Pn ðxÞ dPn ðxÞ 2x 1 x2 þ nðn þ 1ÞPn ðxÞ ¼ 0 2 dx dx

(1:174)

An alternative, simpler form for this equation is obtained by combining the ﬁrst two terms:

d 2 dPn ðxÞ (1:175) 1x þ nðn þ 1ÞPn ðxÞ ¼ 0 dx dx This is the Legendre differential equation. It has a family of solutions, each of which is a polynomial corresponding to a particular value of n. The Legendre polynomials provide solutions in potential analyses with spherical symmetry, and have an important role in geophysical theory. Some Legendre polynomials of low degree are listed in Table 1.1.

1.13.1 Orthogonality of the Legendre polynomials Two vectors a and b are orthogonal if their scalar product is zero: a · b ¼ ax bx þ ay by þ az bz ¼

3 X i¼1

ai bi ¼ 0

(1:176)

38

Mathematical background

By analogy, two functions of the same variable are said to be orthogonal if their product, integrated over a particular range, is zero. For example, the trigonometric functions sin θ and cos θ are orthogonal for the range 0 ≤ θ ≤ 2π, because Z2π

Z2π sin θ cos θ dθ ¼

θ¼0

θ¼0

2π 1 1 sinð2θÞdθ ¼ cosð2θÞ ¼ 0 2 4 θ¼0

(1:177)

The Legendre polynomials Pn(x) and Pl(x) are orthogonal over the range –1 ≤ x ≤ 1. This can be established as follows. First, we write the Legendre equation in short form, dropping the variable x for both Pn and Pl, and, for brevity, writing d Pn ðxÞ ¼ P0n dx

and

d2 Pn ðxÞ ¼ P00n dx2

(1:178)

Thus

1 x2 P00n 2xP0n þ nðn þ 1ÞPn ¼ 0 1 x2 P00l 2xP0l þ lðl þ 1ÞPl ¼ 0

(1:179) (1:180)

Multiplying (1.179) by Pl and (1.180) by Pn gives

1 x2 Pl P00n 2xPl P0n þ nðn þ 1ÞPl Pn ¼ 0 1 x2 Pn P00l 2xP0l Pn þ lðl þ 1ÞPl Pn ¼ 0

(1:181) (1:182)

Subtracting (1.182) from (1.181) gives

1 x2 Pl P00n Pn P00l 2x Pl P0n P0l Pn þ ½nðn þ 1Þ lðl þ 1ÞPl Pn ¼0 (1:183)

Note that d Pl P0n P0l Pn ¼ Pl P00n þ P0l P0n P0l P0n P00l Pn ¼ Pl P00n P00l Pn dx (1:184) and

d Pl P0n P0l Pn 2x Pl P0n P0l Pn dx d 1 x2 Pl P0n P0l Pn ¼ dx

1 x2

(1:185)

1.13 The Legendre differential equation

39

Thus d 1 x2 Pl P0n P0l Pn þ ½nðn þ 1Þ lðl þ 1ÞPl Pn ¼ 0 dx

(1:186)

Now integrate each term in this equation with respect to x over the range –1 ≤ x ≤ 1. We get

1x

2

Pl P0n

P0l Pn

þ1

x¼1

Zþ1 þ ½nðn þ 1Þ lðl þ 1Þ

Pl Pn dx ¼ 0 x¼1

(1:187) The ﬁrst term is zero on evaluation of (1 – x ) at x = ±1; thus the second term must also be zero. For n ≠ l the condition for orthogonality of the Legendre polynomials is Z þ1 Pn ðxÞPl ðxÞdx ¼ 0 (1:188) 2

x¼1

1.13.2 Normalization of the Legendre polynomials A function is said to be normalized if the integral of the square of the function over R þ1 its range is equal to 1. Thus we must evaluate the integral x¼1 ½Pn ðxÞ2 dx. We begin by recalling the generating function for the Legendre polynomials given in (1.156), which we rewrite for Pn(x) and Pl(x) individually: 1 X

1=2 hn Pn ðxÞ ¼ 1 2xh þ h2

(1:189)

1=2 hl Pl ðxÞ ¼ 1 2xh þ h2

(1:190)

n¼0 1 X l¼0

Multiplying these equations together gives 1 X 1 X

1 hnþl Pn ðxÞPl ðxÞ ¼ 1 2xh þ h2

(1:191)

l¼0 n¼0

Now let l = n and integrate both sides with respect to x, taking into account (1.188): 1 X n¼0

Zþ1 h

Zþ1 2

½Pn ðxÞ dx ¼

2n x¼1

x¼1

dx 1 þ h2 2xh

(1:192)

40

Mathematical background

The right-hand side of this equation is a standard integration that results in a natural logarithm: Z dx 1 ¼ lnða þ bxÞ (1:193) a þ bx b The right-hand side of (1.192) therefore leads to þ1 dx 1 2 ¼ ln 1 þ h 2xh 1 þ h2 2xh ð2hÞ x¼1

Zþ1 x¼1

¼

1 ln 1 þ h2 2h ln 1 þ h2 þ 2h 2h

(1:194)

and Zþ1 x¼1

dx 1 1 2 1 2 þ ¼ ln ð 1 h Þ ln ð 1 þ h Þ 1 þ h2 2xh h 2 2 1 ¼ ½lnð1 þ hÞ lnð1 hÞ h

(1:195)

Using the MacLaurin series for the natural logarithms as in (1.135), we get lnð1 þ hÞ ¼ h

h2 h3 h4 hn þ þ ð1Þn1 þ 2 3 4 n

lnð1 hÞ ¼ h

h2 h3 h4 ðhÞn ð1Þn1 þ 2 3 4 n

(1:196) (1:197)

Subtracting the second equation from the ﬁrst gives Zþ1 x¼1

1 dx 2 h3 h5 2X h2nþ1 ¼ þ þ ¼ h þ 2 3 5 1 þ h 2xh h h n¼0 2n þ 1

(1:198)

Inserting this result into (1.192) gives 1 X n¼0 1 X n¼0

0 h2n @

Zþ1 x¼1

Zþ1 h

½Pn ðxÞ2 dx ¼

2n x¼1

1 2X h2nþ1 h n¼0 2n þ 1

1 2 A¼0 ½Pn ðxÞ2 dx 2n þ 1

(1:199)

(1:200)

1.14 Rodrigues’ formula

41

This is true for every value of h in the summation, so we obtain the normalizing condition for the Legendre polynomials: Zþ1 ½Pn ðxÞ2 dx ¼ x¼1

2 2n þ 1

(1:201)

1=2 It follows that n þ 12 Pn ðxÞ is a normalized Legendre polynomial.

1.14 Rodrigues’ formula The Legendre polynomials can be easily computed with the aid of a formula derived by a French mathematician, Olinde Rodrigues (1795–1851). First, we deﬁne the function n fðxÞ ¼ x2 1 (1:202) Differentiating ƒ(x) once with respect to x gives n n1 df d 2 ¼ x 1 ¼ 2nx x2 1 dx dx

(1:203)

Multiplying the result by (x2 − 1) gives 2 d 2 n n x 1 x 1 ¼ 2nx x2 1 dx

(1:204)

2 df ¼ 2nxf x 1 dx

(1:205)

Now we use Leibniz’s rule (1.144) to differentiate both sides of this equation n + 1 times with respect to x. Writing D = d/dx as in Section 1.11, Dnþ1 ðuvÞ ¼

nþ1 X

ðn þ 1Þ! k nþ1k v D u D k!ðn þ 1 kÞ! k¼0

(1:206)

On the left-hand side of (1.205) let u(x) = (x2 − 1) and v(x) = dƒ/dx = Dƒ. Applying Leibniz’s rule, we note that after only three differentiations of (x2 − 1) the result is zero and the series is curtailed. On the right-hand side let u(x) = 2nx and v(x) = ƒ. Note that in this case the series is curtailed after two differentiations. Thus, using Leibniz’s rule to differentiate each side of (1.205) n + 1 times, we get

42

Mathematical background

ðn þ 1Þn n x2 1 Dnþ2 f þ 2xðn þ 1ÞDnþ1 f þ 2 D f ¼ 2nx Dnþ1 f 1·2 þ 2nðn þ 1ÞDn f (1:207)

On gathering terms and bringing all to the left-hand side, we have 2 x 1 Dnþ2 f þ 2x Dnþ1 f nðn þ 1ÞDn f ¼ 0

(1:208)

Now we deﬁne y(x) such that yðxÞ ¼ Dn f ¼

n dn 2 x 1 n dx

(1:209)

and we have

d 2y

x2 1

dx2

þ 2x

dy nðn þ 1Þy ¼ 0 dx

(1:210)

On comparing with (1.174), we see that this is the Legendre equation. The Legendre polynomials must therefore be proportional to y(x), so we can write Pn ðxÞ ¼ cn

n dn 2 x 1 n dx

(1:211)

The quantity cn is a calibration constant. To determine cn we ﬁrst write n d n dn 2 x 1 ¼ n ½ ð x 1Þ n ð x þ 1Þ n dxn dx

(1:212)

then we apply Leibniz’s rule to the product on the right-hand side of the equation: n n X dn 2 n! dm d nm x 1 ¼ ðx 1Þn nm ðx þ 1Þn n m dx dx m!ðn mÞ! dx m¼0

(1:213)

The successive differentiations of (x – 1)n give d ðx 1Þn ¼ nðx 1Þn1 dx d2 ðx 1Þn ¼ nðn 1Þðx 1Þn2 dx2 d n1 ðx 1Þn ¼ ðnðn 1Þðn 2Þ . . . 3·2·1Þðx 1Þ ¼ n!ðx 1Þ dxn1 dn ðx 1Þn ¼ n! dxn

(1:214)

1.15 Associated Legendre polynomials

43

Each differentiation in (1.214) is zero at x = 1, except the last one. Thus each term in the sum in (1.213) is also zero except for the last one, for which m = n. Substituting x = 1 gives n

n n d 2 nd n x 1 ¼ ð x þ 1 Þ ð x 1 Þ ¼ 2n n! (1:215) n dxn dx x¼1 x¼1 Putting this result and the condition Pn(1) = 1 into (1.211) gives n n d 2 x 1 ¼ cn 2n n! ¼ 1 Pn ð1Þ ¼ cn dxn x¼1

(1:216)

where cn ¼

1 2n n!

(1:217)

Rodrigues’ formula for the Legendre polynomials is therefore Pn ðxÞ ¼

1 dn 2n n! dxn

x2 1

n

(1:218)

1.15 Associated Legendre polynomials Many physical properties of the Earth, such as its magnetic ﬁeld, are not azimuthally symmetric about the rotation axis when examined in detail. However, these properties can be described using mathematical functions that are based upon the Legendre polynomials described in the preceding section. To derive these functions, we start from the Legendre equation, (1.174), which can be written in shorthand form as 1 x2 P 00n 2xP0n þ nðn þ 1ÞPn ¼ 0 (1:219) Now we differentiate this equation with respect to x: d 00 d d P 2xP00n 2x P0n 2P0n þ nðn þ 1Þ Pn ¼ 0 1 x2 dx n dx dx

(1:220)

On noting that we can equally write P 00n ¼ ðd=dxÞP0n and P0n ¼ ðd=dxÞPn , this can be written alternatively as

1 x2

d 00 d d P 4x P0n þ ½nðn þ 1Þ 2 Pn ¼ 0 dx n dx dx

which can be written, for later comparison,

(1:221)

44

Mathematical background d 00 d d Pn 2ð2Þx P0n þ ½nðn þ 1Þ 1ð2Þ Pn ¼ 0 1 x2 dx dx dx

(1:222)

Next, we differentiate this expression again, observing the same rules and gathering terms,

d 2 00 d 00 d 0 d2 0 d2 1x P 2x þ 4x P Pn ¼ 0 4 þ ½ n ð n þ 1 Þ2 P P n n n n dx2 dx2 dx2 dx dx 2

(1:223) d 2 00 d2 0 d2 0 d2 P 2x P 4x P þ ½ n ð n þ 1 Þ 2 4 Pn ¼ 0 1 x2 dx2 n dx2 n dx2 n dx2 (1:224)

1 x2

d 2 00 d2 0 d2 P 6x P þ ½ n ð n þ 1 Þ 6 Pn ¼ 0 n n dx2 dx2 dx2

(1:225)

which, as we did with (1.222), we can write in the form

1 x2

d 2 00 d2 d2 Pn 2ð3Þx 2 P0n þ ½nðn þ 1Þ 2ð3Þ 2 Pn ¼ 0 2 dx dx dx

(1:226)

On following step-by-step in the same manner, we get after the third differentiation

1 x2

d 3 00 d3 0 d3 P 2 ð 4 Þx P þ ½ n ð n þ 1 Þ 3 ð 4 Þ Pn ¼ 0 n n dx3 dx3 dx3

(1:227)

Equations (1.222), (1.226), and (1.227) all have the same form. The higherorder differentiation is accompanied by systematically different constants. By extension, differentiating (1.219) m times (where m ≤ n) yields the differential equation

1 x2

d m 00 dm dm Pn 2ðm þ 1Þx m P0n þ ½nðn þ 1Þ mðm þ 1Þ m Pn ¼ 0 m dx dx dx (1:228)

Now let the mth-order differentiation of Pn be written as dm QðxÞ Pn ðxÞ ¼ m=2 dxm ð1 x2 Þ

(1:229)

Substitution of this expression into (1.228) gives a new differential equation involving Q(x). We need to determine both ðd m =dxm ÞP0n and ðd m =dxm ÞP00n , so ﬁrst we differentiate (1.229) with respect to x:

1.15 Associated Legendre polynomials

dm 0 Q0 P ¼ n m=2 dxm ð1 x2 Þ

! m Q ð2xÞ m=2þ1 2 ð1 x 2 Þ

ðmþ2Þ=2 dm 0 P n ¼ 1 x2 1 x2 Q0 þ mxQ m dx

45

(1:230)

(1:231)

A further differentiation of (1.231) by parts gives ðmþ2Þ=2 d 1 x2 1 x2 Q0 þ mxQ dx ðmþ2Þ=2 d 1 x2 Q0 þ mxQ þ 1 x2 dx 2 ðmþ2Þ=21 ¼ ðm þ 2Þx 1 x 1 x2 Q0 þ mxQ ðmþ2Þ=2 1 x2 Q00 þ mxQ0 2xQ0 þ mQ þ 1 x2

d m 00 P ¼ dxm n

(1:232) n d m 00 2 ðmþ2Þ=2 P ¼ 1 x 1 x2 Q00 þ ðm 2ÞxQ0 þ mQ þ ðm þ 2ÞxQ0 n dxm mðm þ 2Þx2 Q (1:233) þ 1 x2 00 d m 00 mðm þ 2Þx2 Q 2 ðmþ2Þ=2 2 0 P ¼ 1 x 1 x þ 2mxQ þ mQ þ Q dxm n 1 x2 (1:234)

Now we substitute (1.231) and (1.234) into (1.228). Unless the multiplier (1 – x2)− (m + 2)/2 is always zero, Q must satisfy the following equation: 2 1 x2 Q00 þ 2mx 1 x2 Q0 þ m 1 x2 Q þ mðm þ 2Þx2 Q 2ðm þ 1Þx 1 x2 Q0 2mðm þ 1Þx2 Q þ ½nðn þ 1Þ mðm þ 1Þ 1 x2 Q ¼ 0

(1:235)

The remainder of the evaluation consists of gathering and reducing terms; we ﬁnally get

m2 Q¼0 1 x Q 2xQ þ nðn þ 1Þ 1 x2 2

00

0

(1:236)

46

Mathematical background

The functions Q(x) involve two parameters, the degree n and order m, and are written Pn,m(x). Thus

d2 d m2 Pn;m ðxÞ ¼ 0 P 1 x2 P ð x Þ 2x ð x Þ þ n ð n þ 1 Þ n;m n;m dx2 1 x2 dx (1:237) This is the associated Legendre equation. The solutions Pn,m(x) or Pn,m(cos θ), where x = cos θ, are called associated Legendre polynomials, and are obtained from the ordinary Legendre polynomials using the deﬁnition of Q in (1.229): m=2 d m Pn;m ðxÞ ¼ 1 x2 Pn ðxÞ dxm

(1:238)

Substituting Rodrigues’ formula (1.218) for Pn(x) into this equation gives

m=2 n 1 x2 d nþm 2 x 1 Pn;m ðxÞ ¼ n nþm 2 n! dx

(1:239)

The highest power of x in the function (x2 − 1)n is x2n. After 2n differentiations the result will be a constant, and a further differentiation will give zero. Therefore n + m ≤ 2n, and possible values of m are limited to the range 0 ≤ m ≤ n.

1.15.1 Orthogonality of associated Legendre polynomials For succinctness we again write Pn,m(x) as simply Pn,m. The deﬁning equations for the associated Legendre polynomials Pn,m and Pl,m are

0 m2 1 x Pn;m 2x Pn;m þ nðn þ 1Þ Pn;m ¼ 0 1 x2

00 0 m2 2 1 x Pl;m 2x Pl;m þ lðl þ 1Þ Pl;m ¼ 0 1 x2 2

00

(1:240)

(1:241)

As for the ordinary Legendre polynomials, we multiply (1.240) by Pl,m and (1.241) by Pn,m:

00 0 m2 Pn;m Pl;m ¼ 0 1 x2 Pn;m Pl;m 2x Pn;m Pl;m þ nðn þ 1Þ 1 x2 (1:242)

1.15 Associated Legendre polynomials

47

00 0 m2 1 x2 Pl;m Pn;m 2x Pl;m Pn;m þ lðl þ 1Þ Pn;m Pl;m ¼ 0 1 x2 (1:243)

On subtracting (1.243) from (1.242) we have i h i h 00 00 0 0 1 x2 Pn;m Pl;m Pl;m Pn;m 2x Pn;m Pl;m Pl;m Pn;m þ ½nðn þ 1Þ lðl þ 1ÞPn;m Pl;m ¼ 0

(1:244)

Following the method used to establish the orthogonality of the ordinary Legendre polynomials (Section 1.13.1), we can write this equation as

o 0 0 d n 1 x2 Pn;m Pl;m Pl;m Pn;m þ ½nðn þ 1Þ lðl þ 1ÞPn;m Pl;m dx ¼0 (1:245) On integrating each term with respect to x over the range –1 ≤ x ≤ 1, we get n

oþ1 0 0 1 x2 Pn;m Pl;m Pl;m Pn;m

x¼1

Zþ1 þ ½nðn þ 1Þ lðl þ 1Þ

Pn;m Pl;m dx ¼ 0

(1:246)

x¼1

The ﬁrst term is zero on evaluation of (1 – x2) at x = ±1; thus the second term must also be zero. Provided that n ≠ l, the condition of orthogonality of the associated Legendre polynomials is x¼þ1 Z

Pn;m ðxÞPl;m ðxÞdx ¼ 0

(1:247)

x¼1

1.15.2 Normalization of associated Legendre polynomials Squaring the associated Legendre polynomials and integrating over –1 ≤ x ≤ 1 gives x¼þ1 Z

x¼1

2 Pn;m ðxÞ dx ¼

2 ðn þ mÞ! 2n þ 1 ðn mÞ!

(1:248)

48

Mathematical background

The squared functions do not integrate to 1, so they are not normalized. If each polynomial is multiplied by a normalizing function, the integrated squared polynomial can be made to equal a chosen value. Different conditions for this apply in geodesy and geomagnetism. The Legendre polynomials used in geodesy are fully normalized. They are deﬁned as follows: Pm n ðxÞ ¼

2n þ 1 ðn mÞ! 1=2 Pn;m ðxÞ 2 ðn þ mÞ!

(1:249)

The Legendre polynomials used in geomagnetism are partially normalized (or quasi-normalized). Schmidt in 1889 deﬁned this method of normalization so that Pm n ðxÞ

¼

ðn mÞ! 2 ðn þ mÞ!

P0n ðxÞ ¼ Pn;0 ðxÞ;

1=2 Pn;m ðxÞ; m¼0

m 6¼ 0

(1:250) (1:251)

Integration of the squared Schmidt polynomials over the full range –1 ≤ x ≤ 1 gives the value 1 for m = 0 and 1/(2n + 1) for m > 0. Some fully normalized Legendre polynomials and partially normalized Schmidt polynomials are listed in Table 1.2. Table 1.2. Some fully normalized associated Legendre polynomials and partially normalized Schmidt polynomials of low degree and order Pm n ðcos θ Þ; Legendre, fully normalized

Pm n ðcos θ Þ; Schmidt, partially normalized

1

cos θ sin θ 1 ð3 cos2 θ 1Þ 2 3 sin θ cos θ

2

2

3 sin2θ

3

0

3

1

3

2

15 sin2θ cos θ

3

3

15 sin3θ

cos θ sin θ 1 ð3 cos2 θ 1Þ 2pﬃﬃﬃ pﬃﬃ3ﬃ sin θ cos θ 3 sin2 θ 2 1 cos θð5 cos2 θ 3Þ 2pﬃﬃﬃ 6 sin θð5 cos2 θ 1Þ 4 pﬃﬃﬃﬃﬃ 15 15 sin2 θ cos θ p2ﬃﬃﬃﬃﬃ 10 sin3 θ 4

n

m

1 1

0 1

2

0

2

1 cos θð5 cos2 θ 3Þ 2 3 sin θð5 cos2 θ 1Þ 2

1.16 Spherical harmonic functions

49

1.16 Spherical harmonic functions Several geophysical potential ﬁelds – for example, gravitation and geomagnetism – satisfy the Laplace equation. Spherical polar coordinates are best suited for describing a global geophysical potential. The potential can vary with distance r from the Earth’s center and with polar angular distance θ and azimuth (equivalent to co-latitude and longitude in geographic terms) on any concentric spherical surface. The solution of Laplace’s equation in spherical polar coordinates for a potential U may be written (see Section 2.4.5, (2.104)) U¼

1 X n X m Bn m An rn þ nþ1 am n cosðmÞ þ bn sinðmÞ Pn ðcos θÞ (1:252) r n¼0 m¼0

m Here An, Bn, am n , and bn are constants that apply to a particular situation. On the surface of the Earth, or an arbitrary sphere, the radial part of the potential of a point source at the center of the sphere has a constant value and the variation over the surface of the sphere is described by the functions in θ and . We are primarily interested in solutions outside the Earth, for which An is zero. Also we can set the constant Bn equal to Rn+1, where R is the Earth’s mean radius. The potential is then given by 1 X n nþ1 X m R m am U¼ n cosðmÞ þ bn sinðmÞ Pn ðcos θÞ r n¼0 m¼0

(1:253)

m Let the spherical harmonic functions Cm n ðθ; Þ and Sn ðθ; Þ be deﬁned as m Cm n ðθ; Þ ¼ cosðmÞ Pn ðcos θÞ m Sm n ðθ; Þ ¼ sinðmÞ Pn ðcos θÞ

(1:254)

The variation of the potential over the surface of a sphere may be described by these functions, or a more general spherical harmonic function Ym n ðθ; Þ that combines the sine and cosine variations: m Ym n ðθ; Þ ¼ Pn ðcos θ Þ

cosðmÞ sinðmÞ

(1:255)

Like their constituent parts – the sine, cosine, and associated Legendre functions – spherical harmonic functions are orthogonal and can be normalized.

50

Mathematical background

1.16.1 Normalization of spherical harmonic functions m Normalization of the functions Cm n ðθ; Þ and Sn ðθ; Þ requires integrating the squared value of each function over the surface of a unit sphere. The element of surface area on a unit sphere is dΩ = sin θ dθ d (Box 1.3) and the limits of integration are 0 ≤ θ ≤ π and 0 ≤ ≤ 2π. The integral is

ZZ

2 Cm n ðθ; Þ

Zπ

Z2π

dΩ ¼

Cm n ðθ; Þ

2

sin θ dθ d

θ¼0 ¼0

S

Zπ

Z2π

¼

cosðmÞPm n ðcos θ Þ

2

sin θ dθ d (1:256)

θ¼0 ¼0

Let x = cos θ in the associated Legendre polynomial, so that dx = –sin θ dθ and the limits of integration are –1 ≤ x ≤ 1. The integration becomes 8 Z1 > < Z2π x¼1

> :

¼0

9 > Zþ1 = 2 m 2 2 m cos ðmÞd Pn ðxÞ dx ¼ π Pn ðxÞ dx > ;

(1:257)

x¼1

Normalization of the associated Legendre polynomials gives the result in (1.248), thus ZZ

Cm n ðθ; Þ

2

dΩ ¼

S

2π ðn þ mÞ! 2n þ 1 ðn mÞ!

(1:258)

The normalization of the function Sm n ðθ; Þ by this method delivers the same result. Spherical harmonic functions make it possible to express the variation of a physical property (e.g., gravity anomalies, g(θ, )) on the surface of the Earth as an inﬁnite series, such as gðθ; Þ ¼

1 X n X

m m m am n Cn ðθ; Þ þ bn Sn ðθ; Þ

(1:259)

n¼0 m¼0 m The coefﬁcients am n and bn may be obtained by multiplying the function gðθ; Þ m m by Cn ðθ; Þ or Sn ðθ; Þ, respectively, and integrating the product over the surface of the unit sphere. The normalization properties give

1.16 Spherical harmonic functions

51

ZZ 2n þ 1 ðn mÞ! gðθ; Þ Cm n ðθ; ÞdΩ 2π ðn þ mÞ! S ZZ 2n þ 1 ð n m Þ! bm gðθ; Þ Sm n ¼ n ðθ; ÞdΩ 2π ðn þ mÞ!

am n ¼

(1:260)

S

1.16.2 Zonal, sectorial, and tesseral spherical harmonics The spherical harmonic functions Ym n ðθ; Þ have geometries that allow graphic representation of a potential on the surface of a sphere. Deviations of the potential from a constant value form alternating regions in which the potential is larger or smaller than a uniform value. Where the potential surface intersects the spherical surface a nodal line is formed. The appearance of any Ym n ðθ; Þ is determined by ð the distribution of its nodal lines. These occur where Ym n θ; Þ = 0. To simplify the discussion we will associate a constant value of the polar angle θ with a circle of latitude, and a constant value of the azimuthal angle with a circle of longitude. The deﬁnition of the associated Legendre polynomials in (1.239) shows that the equation Pm n ðxÞ = 0 has n – m roots, apart from the trivial solution x = ±1. The variation of the spherical harmonic Ym n ðθ; Þ with latitude θ thus has n – m nodal lines, each a circle of latitude, between the two poles. If additionally m = 0, the potential on the sphere varies only with latitude and there are n nodal lines separating zones in which the potential is greater or less than the uniform value. An example of a zonal spherical harmonic is Y02 ðθ; Þ, shown in Fig. 1.12(a). The solution of Laplace’s equation (1.253) shows that the variation in potential around any circle of latitude is described by the function m ΦðÞ ¼ am n cosðmÞ þ bn sinðmÞ

(1:261)

There are 2m nodal lines where Φ() = 0, corresponding to 2m meridians of longitude, or m great circles. In the special case in which n = m, there are no

(a) zonal, Y20

(b) sectorial, Y55

(c) tesseral, Y54

Fig. 1.12. Appearance of (a) zonal, (b) sectorial, and (c) tesseral spherical harmonics, projected on a meridian plane of the reference sphere.

52

Mathematical background

nodal lines of latitude and the longitudinal lines separate sectors in which the potential is greater or less than the uniform value. An example of a sectorial spherical harmonic is Y55 ðθ; Þ, shown in Fig. 1.12(b). In the general case (m ≠ 0, n ≠ m) the potential varies with both latitude and longitude. There are n – m nodal lines of latitude and m nodal great circles (2m meridians) of longitude. The appearance of the spherical harmonic resembles a patchwork of alternating regions in which the potential is greater or less than the uniform value. An example of a tesseral spherical harmonic is Y45 ðθ; Þ, which is shown in Fig. 1.12(c).

1.17 Fourier series, Fourier integrals, and Fourier transforms 1.17.1 Fourier series Analogously to the representation of a continuous function by a power series (Section 1.10), it is possible to represent a periodic function by an inﬁnite sum of terms consisting of the sines and cosines of harmonics of a fundamental frequency. Consider a periodic function ƒ(t) with period τ that is deﬁned in the interval 0 ≤ t ≤ τ, so that (a) ƒ(t) is ﬁnite within the interval; (b) ƒ(t) is periodic outside the interval, i.e., ƒ(t + τ) = ƒ(t); and (c) ƒ(t) is single-valued in the interval except at a ﬁnite number of points, and is continuous between these points. Conditions (a)–(c) are known as the Dirichlet conditions. If they are satisﬁed, ƒ(t) can be represented as fðtÞ ¼

1 a0 X þ ðan cosðnωtÞ þ bn sinðnωtÞÞ 2 n¼1

(1:262)

where ω = 2π/τ and the factor 12 in the ﬁrst term is included for reasons of symmetry. This representation of ƒ(t) is known as a Fourier series. The orthogonal properties of sine and cosine functions allow us to ﬁnd the coefﬁcients an and bn of the nth term in the series by multiplying (1.262) by sin(nωt) or cos(nωt) and integrating over a full period: 2 an ¼ τ

Zτ=2 fðtÞcosðnωtÞdt t¼τ=2

bn ¼

2 τ

(1:263)

Zτ=2 fðtÞsinðnωtÞdt t¼τ=2

1.17 Fourier series, integrals, and transforms

53

Instead of using trigonometric functions, we can replace the sine and cosine terms with complex exponentials using the deﬁnitions in (1.7), i.e., we write expðinωtÞ þ expðinωtÞ 2 expðinωtÞ expðinωtÞ sinðnωtÞ ¼ 2i

cosðnωtÞ ¼

(1:264)

Using these relationships in (1.262) yields bn ½expðinωtÞ expðinωtÞ 2 2i n¼0 1 1 X an ibn X an þ ibn ¼ expðinωtÞ þ expðinωtÞ (1:265) 2 2 n¼0 n¼0

fðtÞ ¼

1 X an

½expðinωtÞ þ expðinωtÞ þ

The summation indices are dummy variables, so in the second sum we can replace n by –n, and extend the limits of the sum to n = –∞; thus fðtÞ ¼ ¼

1 X an ibn n¼0 1 X

2

0 X an þ ibn expðinωtÞ þ expðinωtÞ 2 n¼1

ðan þ an Þ iðbn bn Þ expðinωtÞ 2 n¼1

(1:266)

If we deﬁne cn as the complex number cn ¼

ðan þ an Þ iðbn bn Þ 2

(1:267)

the Fourier series (1.262) can be written in complex exponential form as fðtÞ ¼

1 X

cn expðinωtÞ

(1:268)

n¼1

In this case the harmonic coefﬁcients cn are given by 1 cn ¼ τ

Zτ=2 fðtÞexpðinωtÞdt t¼τ=2

(1:269)

54

Mathematical background

1.17.2 Fourier integrals and Fourier transforms A Fourier series represents the periodic behavior of a physical property as an inﬁnite set of discrete frequencies. The theory can be extended to represent a function ƒ(t) that is not periodic and is made up of a continuous spectrum of frequencies, provided that the function satisﬁes the Dirichlet conditions speciﬁed above and that it has a ﬁnite energy: Z1

j fðtÞj2 dt51

(1:270)

t¼1

The inﬁnite sum in (1.268) is replaced by a Fourier integral and the complex coefﬁcients cn are replaced by an amplitude function g(ω): Z1 fðtÞ ¼

gðωÞexpðiωtÞdω

(1:271)

ω¼1

where g(ω) is a continuous function, obtained from the equation 1 gðωÞ ¼ 2π

Z1 fðtÞexpðiωtÞdt

(1:272)

t¼1

The transition from Fourier series to Fourier integral is explained in Box 1.5. The function g(ω) is called the forward Fourier transform of ƒ(t), and ƒ(t) is called the inverse Fourier transform of g(ω). Fourier transforms constitute a powerful mathematical tool for transforming a function ƒ(t) that is known in the time domain into a new function g(ω) in the frequency domain.

1.17.3 Fourier sine and cosine transforms A simple but important characteristic of a function is whether it is even or odd. An even function has the same value for both positive and negative values of its argument, i.e., ƒ(–t) = ƒ(t). The cosine of an angle is an example of an even function. The integral of an even function over a symmetric interval about the origin is equal to twice the integral of the function over the positive argument. The sign of an odd function changes with that of the argument, i.e., ƒ(–t) = –ƒ(t). For example, the sine of an angle is an odd function. The integral of an odd function over a symmetric interval about the origin is zero. The product of two odd functions or two even functions is an even function; the product of an odd function and an even function is an odd function.

1.17 Fourier series, integrals, and transforms

55

Box 1.5. Transition from Fourier series to Fourier integral The complex exponential Fourier series for a function ƒ(t) is 1 X

fðtÞ ¼

cn expðinωtÞ

(1)

n¼1

where the complex coefﬁcients cn are given by Zτ=2

1 cn ¼ τ

fðtÞexpðinωtÞdt

(2)

t¼τ=2

In these expressions ω = 2π/τ is the fundamental frequency and τ is the fundamental period. From one value of n to the next, the harmonic frequency changes by δω = 2π/τ, so the factor preceding the second equation can be replaced by 1/τ = δω/(2π). To avoid confusion when we insert (2) into (1), we change the dummy variable of the integration to u, giving Zτ=2

δω cn ¼ 2π

fðuÞexpðinωuÞdu

(3)

u¼τ=2

After insertion, (1) becomes 0 1 Zτ=2 1 X Bδω C fðtÞ ¼ fðuÞexpðinωuÞduA expðinωtÞ @ 2π n¼1 u¼τ=2

0

¼

1 X

Bδω @ 2π n¼1

Zτ=2

1 C fðuÞexpðinωðt uÞÞduA

(4)

u¼τ=2

We now deﬁne the function within the integral as Zτ=2 fðuÞexpðinωðt uÞÞdu

FðωÞ ¼ u¼τ=2

The initial Fourier series becomes

(5)

56

Mathematical background

fðtÞ ¼

1 1 X FðωÞδω 2π ω¼1

(6)

We now let the incremental frequency δω become very small, tending in the limit to zero; this is equivalent to letting the period τ become inﬁnite. The index n is dropped because ω is now a continuous variable; the discrete sum becomes an integral and the function f(t) is Z1

1 f ðt Þ ¼ 2π

FðωÞdω

(7)

ω¼1

while the function F(ω) from (5) becomes Z1 fðuÞexpðiωðt uÞÞdu

FðωÞ ¼

(8)

u¼1

On inserting F(ω) into (7) we get 2 1 3 Z Z1 1 4 fðtÞ ¼ fðuÞexpðiωðt uÞÞdu5dω 2π ω¼1 u¼1 2 1 3 Z Z1 1 4 fðuÞexpðiωuÞdu5 expðiωtÞdω ¼ 2π ω¼1

(9)

u¼1

The quantity in square brackets, on changing the variable from u back to t, is 1 gðωÞ ¼ 2π

Z1 fðtÞexpðiωtÞdt

(10)

t¼1

and the original expression can now be written Z1 fðtÞ ¼

gðωÞexpðiωtÞdω

(11)

ω¼1

The equivalence of these two equations is known as the Fourier integral theorem.

1.17 Fourier series, integrals, and transforms

57

Fourier series that represent odd or even functions consist of sums of sines or cosines, respectively. In the same way, there are sine and cosine Fourier integrals that represent odd and even functions, respectively. Suppose that the function ƒ(t) is even, and let us replace the complex exponential in (1.272) using (1.5): 1 gð ω Þ ¼ 2π

Z1 fðtÞ½cosðωtÞ i sinðωtÞdt

(1:273)

t¼1

The sine function is odd, so, if ƒ(t) is even, the product ƒ(t)sin(ωt) is odd, and the integral of the second term is zero. The product ƒ(t)cos(ωt) is even, and we can convert the limits of integration to the positive interval: Z1

1 gð ω Þ ¼ 2π ¼

1 π

fðtÞcosðωtÞdt

t¼1 Z1

fðtÞcosðωtÞdt

(1:274)

t¼0

Thus, if ƒ(t) is even, then g(ω) is also even. Similarly, one ﬁnds that, if ƒ(t) is odd, g(ω) is also odd. Now we expand the exponential in (1.271) and apply the same conditions of evenness and oddness to the products: Z1 fðtÞ ¼

gðωÞðcosðωtÞ þ i sinðωtÞÞdω ω¼1 Z1

¼2

gðωÞcosðωtÞdω

(1:275)

ω¼0

If we were to substitute (1.275) back into (1.274), the integration would be preceded by a constant 2/π, the product of the two constants in these equations. Equations (1.275) and (1.274) form a Fourier-transform pair, and it does not matter how the factor 2/π is divided between them. We will associate it here entirely with the second equation, so that we have the pair of equations Z1 gðωÞcosðωtÞdω

fðtÞ ¼ ω¼0

2 gðωÞ ¼ π

(1:276)

Z1 fðtÞcosðωtÞdt

t¼0

58

Mathematical background

The even functions ƒ(t) and g(ω) are Fourier cosine transforms of each other. A similar treatment for a function ƒ(t) that is odd leads to a similar pair of equations in which the Fourier transform g(ω) is also odd and Z1 fðtÞ ¼

gðωÞsinðωtÞdω ω¼0

2 gð ω Þ ¼ π

(1:277)

Z1 fðtÞsinðωtÞdt

t¼0

The odd functions ƒ(t) and g(ω) are Fourier sine transforms of each other. further reading Boas, M. L. (2006). Mathematical Methods in the Physical Sciences, 3rd edn. Hoboken, NJ: Wiley, 839 pp. James, J. F. (2004). A Student’s Guide to Fourier transforms, 2nd edn. Cambridge: Cambridge University Press, 135 pp.

2 Gravitation

2.1 Gravitational acceleration and potential The Universal Law of Gravitation deduced by Isaac Newton in 1687 describes the force of gravitational attraction between two point masses m and M separated by a distance r. Let a spherical coordinate system (r, θ, ) be centered on the point mass M. The force of attraction F exerted on the point mass m acts radially inwards towards M, and can be written F ¼ G

mM er r2

(2:1)

In this expression, G is the gravitational constant (6.674 21 × 10−11 m3 kg−1 s−2), er is the unit radial vector in the direction of increasing r, and the negative sign indicates that the force acts inwardly, towards the attracting mass. The gravitational acceleration aG at distance r is the force on a unit mass at that point: aG ¼ G

M er r2

(2:2)

The acceleration aG may also be written as the negative gradient of a gravitational potential UG aG ¼ rUG

(2:3)

The gravitational acceleration for a point mass is radial, thus the potential gradient is given by

∂UG M ¼ G 2 r ∂r

UG ¼ G 59

M r

(2:4) (2:5)

60

Gravitation

In Newton’s time the gravitational constant could not be veriﬁed in a laboratory experiment. The attraction between heavy masses of suitable dimensions is weak and the effects of friction and air resistance relatively large, so the ﬁrst successful measurement of the gravitational constant by Lord Cavendish was not made until more than a century later, in 1798. However, Newton was able to conﬁrm the validity of the inverse-square law of gravitation in 1687 by using existing astronomic observations of the motions of the planets in the solar system. These had been summarized in three important laws by Johannes Kepler in 1609 and 1619. The small sizes of the planets and the Sun, compared with the immense distances between them, enabled Newton to consider these as point masses and this allowed him to verify the inverse-square law of gravitation.

2.2 Kepler’s laws of planetary motion Johannes Kepler (1571–1630), a German mathematician and scientist, formulated his laws on the basis of detailed observations of planetary positions by Tycho Brahe (1546–1601), a Danish astronomer. The observations were made in the late sixteenth century, without the aid of a telescope. Kepler found that the observations were consistent with the following three laws ( Fig. 2.1). 1. The orbit of each planet is an ellipse with the Sun at one focus. 2. The radius from the Sun to a planet sweeps over equal areas in equal intervals of time.

Q

P* Q*

aphelion

r θ

a

S b

P(r, θ ) perihelion

p

Fig. 2.1. Illustration of Kepler’s laws of planetary motion. The orbit of each planet is an ellipse with the Sun at its focus (S); a, b, and p are the semi-major axis, semiminor axis, and semi-latus rectum, respectively. The area swept by the radius to a planet in a given time is constant (i.e., area SPQ equals area SP*Q*); the square of the period is proportional to the cube of the semi-major axis. After Lowrie ( 2007).

2.2 Kepler’s laws of planetary motion

61

3. The square of the period is proportional to the cube of the semi-major axis of the orbit. The fundamental assumption is that the planets move under the inﬂuence of a central, i.e., radially directed force. For a planet of mass m at distance r from the Sun the force F can be written F¼m

d 2r ¼ fðrÞer dt2

(2:6)

The angular momentum h of the planet about the Sun is h¼rm

dr dt

(2:7)

Differentiating with respect to time, the rate of change of angular momentum is dh d dr dr dr d 2r ¼m r ¼m þm r 2 dt dt dt dt dt dt

(2:8)

The ﬁrst term on the right-hand side is zero, because the vector product of a vector with itself (or with a vector parallel to itself) is zero. Thus dh d 2r ¼rm 2 dt dt

(2:9)

On substituting from ( 2.6) and applying the same condition, we have dh ¼ r fðrÞer ¼ fðrÞðr er Þ ¼ 0 dt

(2:10)

This equation means that h is a constant vector; the angular momentum of the system is conserved. On taking the scalar product of h and r, we obtain dr r·h ¼ r· r m (2:11) dt Rotating the sequence of the vectors in the triple product gives r·h ¼ m

dr · ðr rÞ ¼ 0 dt

(2:12)

This result establishes that the vector r describing the position of a planet is always perpendicular to its constant angular momentum vector h and therefore deﬁnes a plane. Every planetary orbit is therefore a plane that passes through the Sun. The orbit of the Earth deﬁnes the ecliptic plane.

62

Gravitation

2.2.1 Kepler’s Second Law Let the position of a planet in its orbit be described by polar coordinates (r, θ) with respect to the Sun. The coordinates are deﬁned so that the angle θ is zero at the closest approach of the planet to the Sun (perihelion). The angular momentum at an arbitrary point of the orbit has magnitude h ¼ mr2

dθ dt

(2:13)

In a short interval of time Δt the radius vector from the Sun to the planet moves through a small angle Δθ and deﬁnes a small triangle. The area ΔA of the triangle is 1 DA ¼ r2 Dθ 2 The rate of change of the area swept over by the radius vector is dA DA 1 2 Dθ ¼ lim ¼ lim r Dt!0 Dt Dt!0 2 dt Dt

(2:14)

(2:15)

dA 1 2 dθ ¼ r dt 2 dt

(2:16)

dA h ¼ dt 2m

(2:17)

On inserting from ( 2.13) we get

Thus the area swept over by the radius vector in a given time is constant. This is Kepler’s Second Law of planetary motion.

2.2.2 Kepler’s First Law If just the gravitational attraction of the Sun acts on the planet (i.e., we ignore the interactions between the planets), the total energy E of the planet is constant. The total energy E is composed of the planet’s orbital kinetic energy and its potential energy in the Sun’s gravitational ﬁeld: 2 2 1 dr 1 dθ S þ mr2 Gm ¼ E m 2 dt 2 dt r

(2:18)

The ﬁrst term here is the planet’s linear (radial) kinetic energy, the second term is its rotational kinetic energy (with mr2 being the planet’s moment of inertia

2.2 Kepler’s laws of planetary motion

63

about the Sun), and the third term is the gravitational potential energy. On writing dr dr dθ ¼ dt dθ dt

(2:19)

and rearranging terms we get 2 2 2 dr dθ S E 2 dθ þr 2G ¼ 2 dθ dt dt r m

(2:20)

Now, to simplify later steps, we make a change of variables, writing u¼

1 r

(2:21)

Then dr d 1 1 du 2 du ¼ ¼ 2 ¼ r dθ dθ u u dθ dθ Substituting from ( 2.22) into ( 2.20) gives 2 dθ 2 du 2 dθ S E r2 þ r2 2G ¼ 2 dt dθ dt r m

(2:22)

(2:23)

With the result of ( 2.13) we have dθ h ¼ r2 dt m dθ 1 h h r ¼ ¼u dt r m m On replacing these expressions, ( 2.23) becomes 2 2 2 h du E 2 h þu 2uGS ¼ 2 m dθ m m

du dθ

2 þ u2 2uGS

m2 Em ¼2 2 2 h h

(2:24)

(2:25)

(2:26)

The rest of the evaluation is straightforward, if painstaking. First we add a constant to each side,

du dθ

2 þ u2 2uGS

2 2 m2 m2 Em m2 þ GS ¼ 2 þ GS h2 h2 h2 h2

(2:27)

64

Gravitation

du dθ

2

m2 þ u GS 2 h

2

2 Em m2 ¼ 2 2 þ GS 2 h h

(2:28)

Next, we move the second term to the right-hand side of the equation, giving

du dθ

du dθ

2

2 2 Em m2 m2 ¼ 2 2 þ GS 2 u GS 2 h h h

2 ¼

2 2 m2 2Eh2 m2 GS 2 1 þ 2 2 3 u GS 2 h GS m h

(2:29)

(2:30)

Now, we deﬁne some combinations of these terms, as follows: u0 ¼ GS

m2 h2

(2:31)

2Eh2 G2 S2 m3

(2:32)

e2 ¼ 1 þ

Using these deﬁned terms, ( 2.30) simpliﬁes to a more manageable form: 2 du ¼ u20 e2 ðu u0 Þ2 (2:33) dθ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ du ¼ u20 e2 ðu u0 Þ2 (2:34) dθ The solution of this equation, which can be tested by substitution, is u ¼ u0 ð1 þ e cos θÞ

(2:35)

The angle θ is deﬁned to be zero at perihelion. The negative square root in ( 2.34) is chosen because, as θ increases, r increases and u must decrease. Let p¼

1 h2 ¼ u0 GSm2

(2:36)

r¼

p 1 þ e cos θ

(2:37)

This is the polar equation of an ellipse referred to its focus, and is the proof of Kepler’s First Law of planetary motion. The quantity e is the eccentricity of the ellipse, while p is the semi-latus rectum of the ellipse, which is half the length of a chord passing through the focus and parallel to the minor axis ( Fig. 2.1). These equations show that three types of trajectory around the Sun are possible, depending on the value of the total energy E in ( 2.18). If the kinetic

2.2 Kepler’s laws of planetary motion

65

energy is greater than the potential energy, the value of E in ( 2.32) is positive, and e is greater than 1; the path of the object is a hyperbola. If the kinetic energy and potential energy are equal, the total energy is zero and e is exactly 1; the path is a parabola. In each of these two cases the object can escape to inﬁnity, and the paths are called escape trajectories. If the kinetic energy is less than the potential energy, the total energy E is negative and the eccentricity is less than 1. In this case (corresponding to a planet or asteroid) the object follows an elliptical orbit around the Sun.

2.2.3 Kepler’s Third Law It is convenient to describe the elliptical orbit in Cartesian coordinates (x, y), centered on the mid-point of the ellipse, instead of on the Sun. Deﬁne the x-axis parallel to the semi-major axis a of the ellipse and the y-axis parallel to the semiminor axis b. The equation of the ellipse in Fig. 2.1 is x2 y2 þ ¼1 a2 b2

(2:38)

The semi-minor axis is related to the semi-major axis by the eccentricity e, so that b2 ¼ a2 1 e2

(2:39)

The distance of the focus of the ellipse from its center is by deﬁnition ae. The length p of the semi-latus rectum is the value of y for a chord through the focus. On setting y = p and x = ae in ( 2.38), we obtain p2 ðaeÞ2 ¼ 1 2 ¼ 1 e2 2 b a 2 p 2 ¼ a 2 1 e2

(2:40) (2:41)

Now consider the application of Kepler’s Second Law to an entire circuit of the elliptical orbit. The area of the ellipse is πab, and the period of the orbit is T, so dA πab ¼ dt T

(2:42)

h 2πab ¼ m T

(2:43)

Using (2.17),

66

Gravitation

From ( 2.36) and ( 2.43) we get the value of the semi-latus rectum, 1 h 2 1 2πab 2 ¼ p¼ GS m GS T

(2:44)

Substituting from ( 2.41) gives 4π 2 a2 b2 4π 2 a4 a 1 e2 ¼ ¼ 1 e2 GST 2 GST 2

(2:45)

After simplifying, we ﬁnally get T 2 4π 2 ¼ a3 GS

(2:46)

The quantities on the right-hand side are constant, so the square of the period is proportional to the cube of the semi-major axis, which is Kepler’s Third Law.

2.3 Gravitational acceleration and the potential of a solid sphere The gravitational potential and acceleration outside and inside a solid sphere may be calculated from the Poisson and Laplace equations, respectively.

2.3.1 Outside a solid sphere, using Laplace’s equation Outside a solid sphere the gravitational potential UG satisﬁes Laplace’s equation ( Section 1.9). If the density is uniform, the potential does not vary with the polar angle θ or azimuth . Under these conditions, Laplace’s equation in spherical polar coordinates ( 2.67) reduces to ∂ 2 ∂UG ¼0 (2:47) r ∂r ∂r This implies that the bracketed quantity that we are differentiating must be a constant, C, r2

∂UG ¼C ∂r

(2:48)

∂UG C ¼ 2 r ∂r

(2:49)

The gravitational acceleration outside the sphere is therefore

2.3 The potential of a solid sphere

67

∂UG C ¼ 2 er ∂r r

(2:50)

aG ðr4RÞ ¼

At its surface the gravitational acceleration has the value ∂UG C aG ðRÞ ¼ ¼ 2 er ∂r R

(2:51)

The boundary condition at the surface of the sphere is that the accelerations determined outside and inside the sphere must be equal there. We use this to derive the value of the constant C. On comparing ( 2.51) and ( 2.60) we have C ¼ GM

(2:52)

On inserting for C in ( 2.50), the gravitational acceleration outside the sphere is aG ðr4RÞ ¼ G

M er r2

(2:53)

The gravitational potential outside the solid sphere is obtained by integrating ( 2.53) with respect to the radius. This gives UG ðr4RÞ ¼ G

M r

(2:54)

2.3.2 Inside a solid sphere, using Poisson’s equation Inside a solid sphere with radius R and uniform density ρ the gravitational potential UG satisﬁes Poisson’s equation ( Section 1.8). Symmetry again requires the use of spherical polar coordinates, and, because the density is uniform, there is no variation of potential with the polar angle θ or azimuth . Poisson’s equation in spherical polar coordinates reduces to 1 ∂ 2 ∂UG ¼ 4πGρ r ∂r r2 ∂r

(2:55)

On multiplying by r2 and integrating with respect to r, we get ∂ 2 ∂UG ¼ 4πGρr2 r ∂r ∂r r2

∂UG 4 ¼ πGρr3 þ C1 ∂r 3

(2:56) (2:57)

This equation has to be valid at the center of the sphere where r = 0, so the constant C1 = 0 and

68

Gravitation ∂UG 4 ¼ πGρr ∂r 3 ∂UG ¼ aG ðr5RÞ ¼ ∂r

(2:58) 4 πGρr er 3

(2:59)

This shows that the gravitational acceleration inside a homogeneous solid sphere is proportional to the distance from its center. At the surface of the sphere, r = R, and the gravitational acceleration is 4 GM aG ðRÞ ¼ πGρR er ¼ 2 er 3 R

(2:60)

where the mass M of the sphere is 4 M ¼ πR3 ρ 3

(2:61)

To obtain the potential inside the solid sphere, we must integrate ( 2.58). This gives 2 UG ¼ πGρr2 þ C2 3

(2:62)

The constant of integration C2 is obtained by noting that the potential must be continuous at the surface of the sphere. Otherwise a discontinuity would exist and the potential gradient (and force) would be inﬁnite. Equating ( 2.54) and ( 2.62) at r = R gives 2 GM 4 πGρR2 þ C2 ¼ ¼ πGρR2 3 R 3

(2:63)

C2 ¼ 2πGρR2

(2:64)

The gravitational potential inside the uniform solid sphere is therefore given by 2 UG ¼ πGρr2 2πGρR2 3

(2:65)

2 UG ¼ πGρ r2 3R2 3

(2:66)

A schematic graph of the variation of the gravitational potential inside and outside a solid sphere is shown in Fig. 2.2.

2.4 Laplace’s equation

69

r /R 0

0

1 inside sphere

2

3

4

outside sphere

0.5 1

U G (R)

1.5 U G(r)

r=R

2

Fig. 2.2. Variation with radial distance r of the gravitational potential inside and outside a solid sphere of radius R. The potential of the surface of the sphere is UG(R).

2.4 Laplace’s equation in spherical polar coordinates In the above examples the sphere was assumed to have uniform density so that only the radial term in Laplace’s equation had to be solved. This is also the case when density varies only with radius. In the Earth, however, lateral variations of the density distribution occur, and the gravitational potential UG is then a solution of the full Laplace equation 1 ∂ 2 ∂UG 1 ∂ ∂UG 1 ∂2 UG ¼0 þ þ r sin θ ∂r ∂θ r2 ∂r r2 sin θ ∂θ r2 sin2 θ ∂2

(2:67)

This equation is solved using the method of separation of variables. This is a valuable mathematical technique, which allows the variables in a partial differential equation to be separated so that only terms in one variable are on one side of the equation and terms in other variables are on the opposite side. A trial solution for UG is UG ðr; θ; Þ ¼

Thank you for interesting in our services. We are a non-profit group that run this website to share documents. We need your help to maintenance this website.

To keep our site running, we need your help to cover our server cost (about $400/m), a small donation will help us a lot.