. 81
( 136 .)


Fig. 18.3 illustrates toroidal geometry. To stress the analogue with the sphere, we use »,
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.

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.

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
∞ ∞
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.

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, ...

{ bc cos(m») + bs sin(m») } sin(nθ) (18.25)
+ mn mn
m=1, 3, 5, ...

The identity
{sin([n ’ 1]θ) + sin([n + 1]θ)} (18.26)
sin(θ) cos(n θ) =
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 (θ).

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-
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 (», θ)

≡ uθ sin(θ) = m
[modi¬ed north-south wind] (18.28)
V vmn Yn (», θ)

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.

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)
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
 »=
all θ [y-axis] (18.30)
u» =
 π
 ux »=’
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

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.

18.10 The Pole Problem: Spherical Coordinates
The Courant-Friedrichs-Levy criterion for the stability of explicit time-stepping algorithms
is that
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


. 81
( 136 .)