which is longitude-like, for the “toroidal” angle and θ, which is analogous to colatitude, for

the “poloidal” angle. All physical solutions must be periodic with a period of 2π in both

coordinates. As shown in the lower half of the ¬gure, it is trivial to “¬‚atten” the torus by

making a slit through a meridian, bending the cut torus into a cylinder, and then making a

second cut parallel to the “equator” to unwrap the cylinder into a rectangle.

CHAPTER 18. SPHERICAL & CYLINDRICAL GEOMETRY

392

Figure 18.3: Toroidal coordinates. (a) The torus. (b) The torus developed into a rectangle.

Because the geometry is periodic in both coordinates, the natural basis set on the torus

is a double Fourier series, i. e. any function on the surface of the torus can be expanded as

∞ ∞

acc acs cos(m») sin(nθ) (18.20)

f (», θ) = cos(m») cos(nθ) +

mn mn

m,n m,n

∞ ∞

asc sin(m») cos(nθ) + ass sin(m») sin(nθ)

+ mn mn

m,n m,n

If f (», θ) is free of singularities on the torus, then the series (18.20) will converge expo-

nentially. Schnack, Baxter, and Caramana (1984) is a good illustration of a code (in plasma

physics) that employs a two-dimensional Fourier series in toroidal coordinates.

The sphere is more complex because it is not a surface that can be unwrapped or “de-

veloped” into a cylinder. Its geometry is fundamentally different from that of the periodic

rectangle shown in Fig. 18.3b. To expand a scalar function on the sphere, we need only

half the terms in the general double Fourier series shown above. To explain which half,

however, requires us to look down on the north pole of the sphere and ¬‚atten the “polar

cap” region into a disk with a local polar coordinate system as shown in Fig. 18.4.

18.8. THE PARITY FACTOR FOR SCALARS: SPHERE VERSUS TORUS 393

Figure 18.4: Polar projections showing the positive (H) and negative (L; cross-hatched) re-

gions for a term gm (θ) sin(m») in the Fourier expansion of a function. The one-dimensional

graphs on the right show how gm (θ) must vary along the short, pole-crossing line segments

marked on the polar plots.

(a) Zonal wavenumber m is odd. The pole-crossing line segment has one end where sin(m»)

is negative and the other where sin(m») is positive. To avoid a discontinuity, gm (θ), must

be antisymmetric about the pole.

(b) Zonal wavenumber m is even, so sin(m») is positive everywhere along the pole-crossing

segment. To avoid a discontinuity, gm (θ) must be symmetric about the pole.

Note that f (», θ) can be written in the form

∞ ∞

(18.21)

f (», θ) = fm (θ) cos(m») + gm (θ) sin(m»)

m=0 m=1

The two diagrams in Fig. 18.4 schematically illustrate two typical terms in the sums: one

for m odd (top) and the other for m even (bottom). Let us follow gm (θ) sin(m») along a

meridian over the pole.

Fig. 18.4a shows that when m is odd, the component gm (θ) sin(m») must always change

sign as the pole is crossed. This sign change depends only on the fact that the zonal

wavenumber m is odd, and is quite independent of the details of gm (θ) and also of whether

m = 1, 3, 5, . . . so long as m is odd. This implies that gm (θ) sin(m») must have a jump dis-

continuity at the pole unless

[m odd] (18.22)

gm (0) = gm (π) = 0

The same diagram applies to fm (θ) cos(m») by rotating it through π/2 degrees, so the

same reasoning implies that fm (θ) must have zeros at both poles for all odd longitudinal

wavenumbers m.

CHAPTER 18. SPHERICAL & CYLINDRICAL GEOMETRY

394

This argument can be extended to higher latitudinal derivatives to show that ALL EVEN

derivatives of fm (θ) and gm (θ) must vanish at both poles. The Fourier series for these func-

tions could in principle include both sin(nθ) and cos(nθ) terms, but whereas the sine terms

individually have all even derivatives zero at the poles, the sum of the cosine terms would

be constrained by an in¬nite number of constraints. Reexpansion of a spherical harmonic

series into double Fourier series shows that the trivial solution to these constraints is in fact

the only solution: the colatitude cosine coef¬cients are always zero when the longitudinal

wavenumber m is odd. The series (18.20) (or one of its derivatives) will have a discontinuity

in f (», θ) unless fm (θ) and gm (θ) are the sums of Fourier sine series in θ when the zonal

wavenumber m is odd.

When m is even, Fig. 18.4b shows that the function must have the same value on a given

meridian on either side of the pole. The analytical proof is that if we compare two points

that are at equal distances from the poles, but on opposite sides, θ = δ for both points, but

» differs by π. When m is even, sin(m») = sin(m[» + π]), so

gm (δ) sin(m») = gm (δ) sin(m[» + π]) for any δ, » [m = 0, 2, 4, . . . ] (18.23)

as shown in the meridional slice graphed on the right side of Fig. 18.4b. However, if this

function is symmetric with respect to the pole on the line segment marked in the ¬gure, it

follows that its derivative must be antisymmetric. (Recall from Chapter 8 that differenti-

ation is a parity-reversing operation.) However, an antisymmetric function is zero at the

point of symmetry, in this case, the pole.

It follows that

d gm d gm

[m even] (18.24)

(0) = (π) = 0

dθ dθ

Extending this argument to higher derivatives shows that all the sine coef¬cients must be

zero when m is even.1

Thus, a scalar function f (», θ) that has no singularities on the sphere may be expanded

as an exponentially convergent Fourier series of the special form

∞

{ ac cos(m») + as sin(m») } cos(nθ)

f (», θ) = mn mn

m=0, 2, 4, ...

n=0

∞

{ bc cos(m») + bs sin(m») } sin(nθ) (18.25)

+ mn mn

m=1, 3, 5, ...

n=1

The identity

1

{sin([n ’ 1]θ) + sin([n + 1]θ)} (18.26)

sin(θ) cos(n θ) =

2

shows that equivalently, one can replace sin(nθ) by a basis whose elements are {sin(θ)

cos(nθ)}. [Exercise for the reader: Derive the relationship between the Fourier coef¬cients

for these two alternative forms of the sine series.] For historical reasons, this sin(θ) -times-

a-cosine series representation is often used to replace sin(nθ) in the terms in the second line

of (18.25). This extracted factor of sin(θ) for m odd is the “parity factor”.

argument that zero derivatives to all orders implies zero coef¬cients is not rigorous because a C ∞ func-

1 The

tion such as exp(’1/[θ(π ’ θ)]) has zero derivatives to all orders at the poles. However, such functions can

be represented on θ ∈ [0, π] by either a cosine series or sine series with a subgeometric but exponential rate of

convergence. Thus it is true, even when such exceptions are considered, that one needs only latitudinal Fourier

terms of a single parity to approximate fm (θ) and gm (θ).

18.9. PARITY II: HORIZONTAL VELOCITIES & OTHER VECTOR COMPONENTS 395

We added the restriction “for a scalar function” in stating (18.25) because the arguments

for a vector function such as the wind velocities are more complicated as explained in the

next section. The conclusions are similar, however: to represent the components of a vector,

one needs only half a general Fourier series, and which components in θ must be kept is

different for even and odd zonal wavenumber.

18.9 Parity II: Horizontal Velocities & Other Vector Compo-

nents

Computing the spherical harmonic expansions of vectors on the sphere is a little tricky.

As explained in Orszag (1974) the three Cartesian velocity components, ux , uy , and uz , all

transform like scalars, that is, may all be expanded directly in spherical harmonics2 . The

same is true for the radial velocity ur in spherical components. (A meteorologist would call

ur the “vertical” component.) However, the horizontal velocity components u» ≡ d»/dt

and uθ ≡ dθ/dt transform differently. As noted by Robert (1966), u» and uθ should be

expanded as 1/ sin(θ) times a spherical harmonic series.

For this reason, most spherical harmonics models replace the zonal and meridional

velocities by

∞

≡ u» sin(θ) = m

[modi¬ed east-west wind] (18.27)

U umn Yn (», θ)

m,n

∞

≡ uθ sin(θ) = m

[modi¬ed north-south wind] (18.28)

V vmn Yn (», θ)

m,n

Orszag (1974) gives an analytical proof of (18.27) and (18.28). We instead will offer a graph-

ical justi¬cation.

Fig. 18.5 illustrates a polar cap view of a typical ¬‚ow generated by a Cartesian velocity

component which is function of spherical components with even wavenumber ” in this

particular case, a ¬‚ow which is a constant, independent of all three spatial coordinates.

When we decompose this ¬‚ow into spherical velocity components as shown in the two

lower panels, we ¬nd that the zonal velocity and meridional velocity components are both

of odd zonal wavenumber. For this example of wavenumber one velocities, each has a

single nodal meridian, shown by the dotted line, where the ¬eld is zero. This interchange

of even wavenumber for odd wavenumber and vice versa when we decompose a Cartesian

vector component into spherical vector components explains why the parity factor for u»

and uθ for odd m is just what we expect for a scalar quantity for even m ” none.

When we inspect a polar view of an odd wavenumber portion of a scalar quantity as

in Fig. 18.4a, we observe that it changes sign as we cross the poles while moving on a

meridian. We argued that the colatitude dependence of the scalar had to be antisymmetric

about θ = 0 ” that is, like sin(θ) ” so the scalar would not have a jump discontinuity at

the poles. Why does this argument fail for u» and uθ for odd m? The answer is that these

spherical vector components actually are discontinuous at the poles!

If we look again at the upper left of Fig. 18.5, we note that what is interpreted as a

westward ¬‚ow in the top of the disk is an eastward ¬‚ow in the bottom half (opposite

hemisphere). The wind velocity, however, does not go to zero at the pole ” the magnitude

and direction of the ¬‚ow are independent of coordinate. If we embed a ¬xed Cartesian

2 Warning: in this section, subscripts will be used to denote a particular component of a vector and not differ-

entiation with respect to that coordinate.

CHAPTER 18. SPHERICAL & CYLINDRICAL GEOMETRY

396

Figure 18.5: Polar cap view of a ¬‚ow parallel to the Cartesion x-axis, chosen to be represen-

tative of even zonal wavenumber. In this case, the zonal wavenumber is zero ” the ¬‚ow

is independent of position.

When decomposed into vector components, however, both the zonal and meridional ve-

locity have odd zonal wavenumber, in this case, m = 1. To avoid a discontinuity at the pole,

both spherical velocity components must be symmetric about the pole even though a scalar

quantity of odd wavenumber would have to be antisymmetric.

coordinate system in the globe so that the origin is at the north pole, then along the y-axis

the north-south ¬‚ow is zero and

π

|u» | = |ux | for » = ± , all θ [y-axis] (18.29)

2

so that the magnitude of the zonal velocity is independent of θ along the whole meridian.

The only way that this constancy of magnitude can be reconciled with the change of sign

is if u» has a jump discontinuity as the pole is crossed:

± π

’ux

»=

2

all θ [y-axis] (18.30)

u» =

π

ux »=’

2

This in turn requires that we omit the sin(θ) parity that is needed for well-behaved quanti-

ties of odd zonal wavenumber (such as scalars).

It is important to note, however, that this discontinuity is caused by a change in math-

ematical interpretation rather than a physical variation in the ¬‚ow itself. In Fig. 18.5, the

wind as measured by an anemometer never changes either magnitude or direction. What

does change is that we interpret the same wind as an eastward wind on one side of the

pole but as a westward current on the other side. Thus, the prescription against a physical

18.9. PARITY II: HORIZONTAL VELOCITIES & OTHER VECTOR COMPONENTS 397

discontinuity remains in force, and nothing in this section should be (mis)interpreted as

implying otherwise.

The meridional velocity component has a similar jump discontinuity on the x-axis (and

every other meridian except the y-axis where uθ ≡ 0). The modi¬ed velocity components

U and V in contrast, are continuous across the poles ” another reason for preferring them

to u» and uθ . However, the radial velocity is not discontinuous as we cross the pole: up

away from the center of the earth is still up.

Figure 18.6: Same as previous ¬gure except that the Cartesian velocity is now a function

of odd zonal wavenumber, in this case, m = 1. The spherical velocity components are

the sums of terms of even zonal wavenumber. To avoid polar discontinuities, the velocity

components must vanish at the poles. A scalar quantity of similar zonal wavenumber would

be symmetric with respect to the poles.

Fig. 18.6 illustrates the situation when the spherical velocity components are of even

zonal wavenumber. For the case shown, the zonal velocity is eastward over the whole

globe. The meridional velocity alternates in sign, but if we follow a meridian over the

pole, we ¬nd that a poleward ¬‚ow is a poleward ¬‚ow on both sides of the pole. However,

this in turn requires that uθ decrease smoothly to zero as the pole is approached along any

meridian. This in turn requires the parity factor, sin(θ). This may be seen in the analytical

expressions for uθ and uθ at the bottom of Fig. 18.6.

Swarztrauber (1981, 1993) presents some good illustrations of discontinuities in hori-

zontal wind ¬elds. His proposed remedy, vector spherical harmonics, will be discussed in

Sec. 18.22.

CHAPTER 18. SPHERICAL & CYLINDRICAL GEOMETRY

398

18.10 The Pole Problem: Spherical Coordinates

The Courant-Friedrichs-Levy criterion for the stability of explicit time-stepping algorithms

is that

x

(18.31)

t<

c

where c is the speed of the fastest waves allowed by the differential equations being solved

and x is the smallest spatial grid interval. One reason why nearly-uniform ¬nite difference