the Coriolis parameter. Time-evolution is determined by the material
conservation of surface potential vorticity (Bretherton, 1966),

$\theta=-\frac{1}{\sigma_{0}^{2}}\frac{\partial\psi}{\partial
z}\Big{|}_{z=0},$ (3.3)

at the upper boundary, $z=0$, where $\sigma_{0}=\sigma(0)$. Explicitly, the
time-evolution equation is

$\frac{\partial\theta}{\partial
t}+J(\psi,\theta)+\Lambda\,\partial_{x}\theta=F-D,$ (3.4)

at $z=0$, where
$J(\psi,\theta)=\partial_{x}\psi\,\partial_{y}\theta-\partial_{x}\theta\,\partial_{y}\,\psi$
represents the advection of $\theta$ by the geostrophic velocity,
$\bm{u}=\hat{\bm{z}}\times\bm{\nabla}\psi$. The frequency, $\Lambda$, is given
by

$\Lambda=\frac{1}{\sigma_{0}^{2}}\frac{\mathrm{d}U}{\mathrm{d}z}\Big{|}_{z=0},$
(3.5)

where $U(z)$ is a background zonal geostrophic flow. Without loss of
generality, we have assumed that $U(0)=0$ in the time-evolution equation (3.4)
to eliminate a constant advective term. The dissipation, $D$, consists of
linear damping and small-scale dissipation,

$D=r\,\theta+\mathrm{ssd},$ (3.6)

where $r$ is the damping rate. The forcing, $F$, and the small-scale
dissipation, $\mathrm{ssd}$, are described in section 3.4.

The surface buoyancy anomaly, $b|_{z=0}$, is related to the surface potential
vorticity, $\theta$, through

$b|_{z=0}=-f\,\sigma_{0}^{2}\,\theta.$ (3.7)

Therefore, the time-evolution equation (3.4) equivalently states that surface
buoyancy anomalies are materially conserved in the absence of forcing and
dissipation. In addition, the frequency, $\Lambda$, corresponds to a
meridional buoyancy gradient,

$\frac{\mathrm{d}B}{\mathrm{d}y}\Big{|}_{z=0}=-f\,\sigma_{0}^{2}\,\Lambda,$
(3.8)

where $B(y,z)$ is the buoyancy field that is in geostrophic balance with
background zonal velocity, $U(z)$.

If we further assume a doubly periodic domain in the horizontal, then we can
expand the streamfunction as

$\psi(\bm{r},z,t)=\sum_{\bm{k}}\hat{\psi}_{\bm{k}}(t)\,\Psi_{k}(z)\,\mathrm{e}^{\mathrm{i}\bm{k}\cdot\bm{x}},$
(3.9)

where $\bm{x}=(x,y)$ is the horizontal position vector, $z$ is the vertical
coordinate, $\bm{k}=(k_{x},k_{y})$ is the horizontal wavevector,
$k=\left\lvert\bm{k}\right\rvert$ is the horizontal wavenumber, and $t$ is the
time coordinate. The non-dimensional wavenumber-dependent vertical structure,
$\Psi_{k}(z)$, is determined by the boundary value problem (chapter 2)

$-\frac{\mathrm{d}}{\mathrm{d}z}\left(\frac{1}{\sigma^{2}}\frac{\mathrm{d}\Psi_{k}}{\mathrm{d}z}\right)+k^{2}\,\Psi_{k}(z)=0,$
(3.10)

with the upper boundary condition

$\Psi_{k}(0)=1,$ (3.11)

and lower boundary condition

$\Psi_{k}\rightarrow 0\quad\textrm{as }\quad z\rightarrow-\infty.$ (3.12)

The upper boundary condition (3.11) is a normalization for the vertical
structure, $\Psi_{k}(z)$, chosen so that

$\psi(\bm{r},z=0,t)=\sum_{\bm{k}}\hat{\psi}_{\bm{k}}(t)\,\mathrm{e}^{\mathrm{i}\bm{k}\cdot\bm{x}}.$
(3.13)

The corresponding Fourier expansion of the surface potential vorticity is
given by

$\theta(\bm{r},t)=\sum_{\bm{k}}\hat{\theta}_{\bm{k}}(t)\,\mathrm{e}^{\mathrm{i}\bm{k}\cdot\bm{x}},$
(3.14)

where

$\hat{\theta}_{\bm{k}}=-m(k)\,\hat{\psi}_{\bm{k}},$ (3.15)

and the function $m(k)$ is given by

$m(k)=\frac{1}{\sigma_{0}^{2}}\frac{\mathrm{d}\Psi_{k}(0)}{\mathrm{d}z}.$
(3.16)

The function $m(k)$ relates $\hat{\theta}_{\bm{k}}$ to $\hat{\psi}_{\bm{k}}$
in the Fourier space inversion relation (3.15) and so we call $m(k)$ the
_inversion function_.

To recover the well-known case of the uniformly stratified quasigeostrophic
model (Held et al., 1995), set $\sigma(z)=\sigma_{0}$. Then the vertical
structure equation (3.10) along with boundary conditions (3.11) and (3.12)
yield the exponentially decaying vertical structure
$\Psi_{k}(z)=\exp\left({\sigma_{0}\,k\,z}\right)$. On substituting
$\Psi_{k}(z)$ into equation (3.16), we obtain a linear inversion function

$m(k)=\frac{k}{\sigma_{0}}$ (3.17)

and hence [from the inversion relation (3.15)] a linear-in-wavenumber
inversion relation
$\hat{\theta}_{\bm{k}}=-(k/\sigma_{0})\,\hat{\psi}_{\bm{k}}$.

####  3.2.2 The inversion function and spatial locality

Figure 3.1: The inversion functions, $m(k)$ [in panel (a)] for two
stratification profiles [panel (b)] given by the piecewise stratification
profile (3.18). One stratification profile is increasing
[$\sigma^{\prime}(z)\geq 0$, blue], with $\sigma_{0}=1$,
$\sigma_{\mathrm{pyc}}=0.15$, $h_{\mathrm{mix}}=0.01$, and
$h_{\mathrm{lin}}=0.05$. The other stratification profile is decreasing
[$\sigma^{\prime}(z)\leq 0$, red] with $\sigma_{0}=1$,
$\sigma_{\mathrm{pyc}}=10$, $h_{\mathrm{mix}}=0.01$, and
$h_{\mathrm{lin}}=0.05$. The thin black line is given by $k/\sigma_{0}$ where
$\sigma_{0}=1$, whereas the blue and red lines are given by
$k/\sigma_{\mathrm{pyc}}$ with $\sigma_{\mathrm{pyc}}=0.15$ for the thin blue
line and $\sigma_{\mathrm{pyc}}=10$ for the thin red line.

The inversion function $m(k)$, which is determined by the stratification’s
vertical structure, controls the spatial locality of the resulting turbulence.
We illustrate this point with the following piecewise stratification profile,

$\sigma(z)=\begin{cases}\sigma_{0}\quad&\text{for}\quad-
h_{\mathrm{mix}}<z<0\\\
\sigma_{0}+\Delta\sigma\left(\frac{z+h_{\mathrm{mix}}}{h_{\mathrm{lin}}}\right)\quad&\text{for}\quad-(h_{\mathrm{mix}}+h_{\mathrm{lin}})<z<-h_{\mathrm{mix}}\\\
\sigma_{\mathrm{pyc}}\quad&\text{for}\quad-\infty<z<-(h_{\mathrm{mix}}+h_{\mathrm{lin}}),\end{cases}$
(3.18)

where $\Delta\sigma=\sigma_{0}-\sigma_{\mathrm{pyc}}$. At small horizontal
scales, where $k\gg k_{s}$, and

$k_{s}=1/\left(\sigma_{0}\,h_{\mathrm{mix}}\right),$ (3.19)

then $m(k)\approx k/\sigma_{0}$, as in the uniformly stratified model of Held
et al. (1995). Likewise, in the large-scale limit, where $k\ll
k_{\mathrm{pyc}}$, and

$k_{\mathrm{pyc}}=\begin{cases}1/\left(\sigma_{\mathrm{pyc}}\,h_{\mathrm{mix}}\right)\quad&\text{for}\quad\Delta\sigma\leq
0\\\
\sigma_{\mathrm{pyc}}/\left(\sigma_{0}^{2}\,h_{\mathrm{mix}}\right)\quad&\text{for}\quad\Delta\sigma>0,\end{cases}$
(3.20)

then $m(k)\approx k/\sigma_{\mathrm{pyc}}$. However, for wavenumbers between
$k_{\mathrm{pyc}}\lesssim k\lesssim k_{s}$, the inversion function takes an
approximate power law form

$m(k)\approx m_{0}\,k^{\alpha},$ (3.21)

where $m_{0}>0$ and $\alpha\geq 0$. The power $\alpha$ depends on the ratio
$\sigma_{\mathrm{pyc}}/\sigma_{0}$ between the deep and surface
stratification. If the stratification decreases towards the surface
[$\sigma^{\prime}(z)\leq 0$, or $\sigma_{\mathrm{pyc}}/\sigma_{0}>1$] then
$\alpha>1$, with $\sigma_{\mathrm{pyc}}/\sigma_{0}\rightarrow\infty$ sending
$\alpha\rightarrow 2$. In contrast, if the stratification increases towards
the surface [$\sigma^{\prime}(z)\geq 0$, or
$\sigma_{\mathrm{pyc}}/\sigma_{0}<1$] then $\alpha<1$, with
$\sigma_{\mathrm{pyc}}/\sigma_{0}\rightarrow 0$ sending $\alpha\rightarrow 0$.
Thus, for wavenumbers $k_{\mathrm{pyc}}\lesssim k\lesssim k_{s}$, the
inversion relation (3.15) has the approximate form

$\hat{\xi}_{\bm{k}}=-k^{\alpha}\,\hat{\psi}_{\bm{k}},$ (3.22)

where $\hat{\xi}_{\bm{k}}=\hat{\theta}_{\bm{k}}/m_{0}$, which is the inversion
relation for $\alpha$-turbulence (Pierrehumbert et al., 1994; Smith et al.,
2002; Sukhatme and Smith, 2009). Figure 3.1 provides two examples, one with
decreasing stratification (with $\alpha\approx 1.50$) and another with
increasing stratification (with $\alpha\approx 0.49$).

To see how the parameter $\alpha$ modifies the resulting dynamics, consider a
point vortex at the origin, given by
$\xi=\delta(\left\lvert\bm{x}\right\rvert)$, where
$\left\lvert\bm{x}\right\rvert$ is the horizontal distance from the vortex
centre, and $\delta(\left\lvert\bm{x}\right\rvert)$ is the Dirac delta. If
$\alpha=2$, then the streamfunction induced by the point vortex is
logarithmic,
$\psi(\left\lvert\bm{x}\right\rvert)=\log(\left\lvert\bm{x}\right\rvert)/(2\pi)$.
If $0<\alpha<2$, then
$\psi(\left\lvert\bm{x}\right\rvert)=-C_{\alpha}/\left\lvert\bm{x}\right\rvert^{2-\alpha}$
where $C_{\alpha}>0$ is a constant (Iwayama and Watanabe, 2010). Smaller
$\alpha$ leads to vortices with velocities decaying more quickly with the
horizontal distance $\left\lvert\bm{x}\right\rvert$, and hence a shorter
interaction range. Thus, the vertical stratification modifies the relationship
between a surface buoyancy anomaly and its induced velocity field: a surface
buoyancy anomaly over decreasing stratification [$\sigma^{\prime}(z)\leq 0$]
generates a longer range velocity field than an identical buoyancy anomaly
over increasing stratification [$\sigma^{\prime}(z)\geq 0$].

####  3.2.3 Wave dispersion in variable stratification

The background gradient term, $\Lambda$, in the time-evolution equation (3.4)
allows for the propagation of surface-trapped Rossby waves. Substituting a
wave solution of the form
$\psi(x,z,t)=\Psi_{k}(z)\,\exp\left[{\mathrm{i}\left(\bm{k}\cdot\bm{r}-\omega
t\right)}\right]$, where the vertical structure $\Psi_{k}(z)$ satisfies the
boundary value problem (3.10)–(3.12), into the time-evolution equation (3.4)
yields the angular frequency

$\omega(\bm{k})=-\frac{\Lambda\,k_{x}}{m(k)}.$ (3.23)

Given the relationship (3.8) between the meridional surface buoyancy gradient
$\mathrm{d}B/\mathrm{d}y|_{z=0}$ and the frequency $\Lambda$, a poleward
decreasing buoyancy gradient ($f\mathrm{d}{B}/\mathrm{d}{y}<0$) implies
westward propagating $(\omega<0)$ Rossby waves.

The dispersion relation (3.23) shows that Rossby wave dispersion is coupled to
the flow’s interaction range and hence the stratification’s vertical
structure. If we approximate the inversion function as a power law (3.21)
between $k_{\mathrm{pyc}}\lesssim k\lesssim k_{s}$, then the zonal phase
speed, $c=\omega/k_{x}$, becomes $c\sim 1/k^{\alpha}$. Therefore, at these
horizontal scales, Rossby waves are more dispersive over decreasing
stratification (with $\alpha>1$) than over increasing stratification (with
$\alpha<1$). In the limit that $\sigma_{0}\gg\sigma_{\mathrm{pyc}}$ in which
$\alpha\rightarrow 0$, then $c\approx$ constant, and so Rossby waves become
non-dispersive.

###  3.3 From edge waves to surface-trapped jets

The emergence of jets in barotropic $\beta$-plane turbulence is due to two
properties of the potential vorticity (Dritschel and McIntyre, 2008; Scott and
Dritschel, 2019). The first is the resilience of strong latitudinal potential
vorticity gradients to mixing (i.e., ”Rossby wave elasticity”, Dritschel and
McIntyre, 2008). Regions with weak latitudinal potential vorticity gradients
are preferentially mixed, weakening the gradient in these regions and
enhancing the gradient in regions where the latitudinal potential vorticity
gradient is already strong (Dritschel and Scott, 2011). The ultimate limit of
such latitudinally inhomogeneous mixing is a potential vorticity staircase
(Danilov and Gryanik, 2004; Dritschel and McIntyre, 2008; Scott and Dritschel,
2012), which consists of uniform regions of potential vorticity punctuated by
sharp potential vorticity gradients. The second property is that, through
potential vorticity inversion, strong (positive) latitudinal gradients in
potential vorticity correspond to eastward jets. Therefore, inverting a
potential vorticity staircase produces a flow with eastward zonal jets centred
at the sharp frontal zones, with weaker westward flows in between (Scott and
Dritschel, 2019).

However, the limit of a potential vorticity staircase is only achieved for
sufficiently large values of the non-dimensional number
$k_{\varepsilon}/k_{\mathrm{Rh}}$ (Scott and Dritschel, 2012), which is a
ratio of the forcing intensity wavenumber, $k_{\varepsilon}$, to the Rhines
wavenumber, $k_{\mathrm{Rh}}$. The forcing intensity wavenumber is given by
(Maltrud and Vallis, 1991)

$k_{\varepsilon}=(\beta^{3}/\varepsilon_{\mathcal{K}})^{1/5},$ (3.24)

where $\varepsilon_{\mathcal{K}}$ is the kinetic energy injection rate in the
barotropic model, and is obtained by setting the turbulent strain rate equal
to the Rossby wave frequency (Vallis and Maltrud, 1993). The Rhines wavenumber
is given by (Rhines, 1975)

$k_{\mathrm{Rh}}=\sqrt{\beta/U_{\mathrm{rms}}},$ (3.25)

where $U_{\mathrm{rms}}$ is the rms velocity. Scott and Dritschel (2012) found
that the ratio $k_{\varepsilon}/k_{\mathrm{Rh}}$ controls the structure of
zonal jets in barotropic $\beta$-plane turbulence; as
$k_{\varepsilon}/k_{\mathrm{Rh}}$ is increased, the zonal jet strength
increases and the potential vorticity gradient at the jet core becomes larger,
with the staircase limit approached as $k_{\varepsilon}/k_{\mathrm{Rh}}\sim
O(10)$.

Jet formation in surface quasigeostrophic turbulence proceeds similarly, with
the surface buoyancy (which is proportional to $\theta$) taking the role of
the potential vorticity and the frequency, $\Lambda$, taking the role of the
potential vorticity gradient, $\beta$. In this section, we first derive a non-
dimensional number analogous to $k_{\varepsilon}/k_{\mathrm{Rh}}$ for surface
quasigeostrophy. Then we consider how vertical stratification (and the non-
locality parameter $\alpha$) modifies jet structure in the buoyancy staircase
limit, as well as how it modifies the relationship between the Rhines
wavenumber and the jet spacing.

Before proceeding, we comment on two differences between two-dimensional
barotropic turbulence and its surface quasigeostrophic counterpart. First, in
the absence of forcing and dissipation, the kinetic energy,

$\mathcal{K}=-\frac{1}{2}\,\overline{\psi\nabla^{2}\psi}=\frac{1}{2}\,\overline{\left\lvert
u\right\rvert^{2}},$ (3.26)

is a conserved constant in two-dimensional barotropic turbulence (the overline
denotes an area average). With a constant kinetic energy injection rate,
$\varepsilon_{\mathcal{K}}$, and a linear damping rate, $r$, the equilibrium
kinetic energy is $\mathcal{K}=\varepsilon_{\mathcal{K}}/2r$. By definition,
the rms velocity is given by $U_{\mathrm{rms}}=\sqrt{2\mathcal{K}}$. Combining
this expression with the definition of the kinetic energy (3.26) and
substituting into the definition of the Rhines wavenumber (3.25) yields a
Rhines wavenumber expressed in terms of external parameters alone,

$k_{\mathrm{Rh}}=\beta^{1/2}(r/\varepsilon_{\mathcal{K}})^{1/4}.$ (3.27)

In contrast, in surface quasigeostrophy, the total energy,

$\mathcal{E}=-\frac{1}{2}\overline{\psi|_{z=0}\,\theta},$ (3.28)

is a conserved constant in the absence of forcing and dissipation and there is
no general relationship between the rms velocity, $U_{\mathrm{rms}}$, and the
equilibrium total energy, $\mathcal{E}=\varepsilon/2r$, where $\varepsilon$ is
the total energy injection rate in the surface quasigeostrophic model.
Therefore, we are not generally able to express the Rhines wavenumber in terms
of the external parameters $\varepsilon$, $\Lambda$, and $r$. Second, because
$\mathcal{E}$ and $\mathcal{K}$ have different dimensions, the kinetic energy
injection in the barotropic model, $\varepsilon_{\mathcal{K}}$, has different
dimensions than the total energy injection rate in the surface
quasigeostrophic model, $\varepsilon$. In particular, $\varepsilon$ has
dimensions of $L^{2}/T^{3}$.

####  3.3.1 The forcing intensity wavenumber

To obtain the forcing intensity wavenumber, $k_{\varepsilon}$, we compare the
Rossby wave frequency (3.23) to the turbulent strain rate, $\omega_{s}(k)$. If
the inversion function is not approximately constant (i.e., $\alpha\neq 0$)
then the strain rate is (chapter 2)

$\omega_{s}(k)\sim\varepsilon^{1/3}\,k^{4/3}\,\left[m(k)\right]^{-1/3}.$
(3.29)

In particular, if $m(k)=m_{0}\,k^{\alpha}$, then $\omega_{s}(k)\sim
m_{0}^{1/3}\varepsilon^{1/3}\,k^{\left(4-\alpha\right)/3}$. Setting the
absolute value of the Rossby wave frequency for waves with $k=k_{x}$ equal to
the turbulent strain rate (3.29) yields the condition

$k_{\varepsilon}\left[m(k_{\varepsilon})\right]^{2}\sim\frac{\left\lvert\Lambda\right\rvert^{3}}{\varepsilon}.$
(3.30)

A solution to this equation always exists because $\mathrm{d}m/\mathrm{d}k\geq
0$. If the inversion function takes the power law form (3.21), then we obtain

$k_{\varepsilon}=\left(\frac{\left\lvert\Lambda\right\rvert^{3}}{m_{0}^{2}\,\varepsilon}\right)^{1/\left(2\alpha+1\right)},$
(3.31)

which is equivalent to a wavenumber derived in Smith et al. (2002).

####  3.3.2 The damping rate wavenumber and the Rhines wavenumber

Suppose the inversion function takes an approximate power law form,
$m(k)\approx m_{0}\,k^{\alpha}$, near the energy containing wavenumbers. Then
the generalization of the Rhines wavenumber at these wavenumbers is

$k_{\mathrm{Rh}}=\left(\frac{\Lambda}{m_{0}\,U_{\mathrm{rms}}}\right)^{1/\alpha}.$
(3.32)

However, unlike in two-dimensional barotropic turbulence where
$U_{\mathrm{rms}}=\sqrt{2\mathcal{K}}=\sqrt{\varepsilon_{\mathcal{K}}/r}$, we
do not have a general relationship between $U_{\mathrm{rms}}$ and the external
parameters $r$ and $\varepsilon$ in surface quasigeostrophic turbulence. To
obtain a second wavenumber that depends on the damping rate, $r$, we follow
Smith et al. (2002). From dimensional considerations, the energy spectrum at
small wavenumbers is

$E_{\Lambda}(k)\sim\Lambda^{2}\,k^{-(\alpha+3)}/m_{0}.$ (3.33)

Then, defining $k_{r}$ as the wavenumber at which the inverse cascade halts,
we obtain

$\frac{\varepsilon}{2r}\approx\int_{k_{r}}^{\infty}E(k)\,\mathrm{d}k\approx\left(\frac{\Lambda^{2}/m_{0}}{\alpha+2}\right)k_{r}^{-(\alpha+2)},$
(3.34)

where the second equality follows because the integral is dominated by its
peak at low wavenumbers. Solving for $k_{r}$ and neglecting any non-
dimensional coefficients, we obtain

$k_{r}=\left(\frac{\Lambda^{2}\,r}{m_{0}\,\varepsilon}\right)^{1/(\alpha+2)}.$
(3.35)

Note that the damping rate wavenumber, $k_{r}$, has the same dependence on
$\Lambda$, $\varepsilon$, and $r$ as the Rhines wavenumber, $k_{\mathrm{Rh}}$,
only if $\alpha=2$.

####  3.3.3 Surface potential vorticity inversion

A perfect surface potential vorticity staircase consists of mixed zones of
halfwidth $b$, where $\mathrm{d}\theta/\mathrm{d}y=-\Lambda$, separated by
jump discontinuities at which $\mathrm{d}\theta/\mathrm{d}y=\infty$. We find
it more conveniant to work with the relative surface potential vorticity,
$\theta$, rather than the total surface potential vorticity,
$\theta+\Lambda\,y$. In this case, if the total surface potential vorticity,
$\theta+\Lambda\,y$, is a perfect staircase with step width $2b$, then the
relative surface potential vorticity, $\theta$, is a $2b$-periodic sawtooth
wave.

Our first question is whether such a staircase is possible for general $m(k)$.
To answer this question, we consider the velocity field induced by a jump
discontinuity in $\theta$. For a jump discontinuity in an infinite domain,

$\theta=\begin{cases}\Delta\theta\quad&\text{for}\quad 0<y<\infty\\\
0\quad&\text{for}\quad-\infty<y<0,\end{cases}$ (3.36)

the zonal velocity is given by

$u=\frac{\Delta\theta}{2\pi}\int_{-\infty}^{\infty}\frac{\mathrm{e}^{\mathrm{i}\,k_{y}y}}{m\left(\left\lvert
k_{y}\right\rvert\right)}\mathrm{d}k_{y}.$ (3.37)

If $m(k)=m_{0}\,k^{\alpha}$, then this expression is proportional to
$\left\lvert y\right\rvert^{\alpha-1}$ if $\alpha\neq 1$ and logarithmic
otherwise, and so the zonal velocity diverges at $y=0$ if $\alpha\leq 1$.
Consequently, we expect that a perfect staircase should not be possible over
constant or increasing stratification due to the divergence of the zonal
velocity at a jump discontinuity.

Figure 3.2: Panel (a) shows a sloping sawtooth function (thick black line)
along with its derivative (thin black line). Panel (b) shows the normalized
zonal velocity induced by the sloping sawtooth function in panel (a) for
various values of the parameter $\alpha$. Panel (c) shows the normalized zonal
velocity induced by the sawtooth function in (a) in the increasing (blue line)
and decreasing stratifications (red line) shown in figure 3.1.

We therefore consider the more general case of a sloping staircase, where
there is a finite frontal zone of width $2a$ between the mixed zones. In this
case, $\theta$ is a $2(a+b)$-periodic sloping sawtooth wave (see figure 3.2),
and is given by the periodic extension of

$\theta=\Lambda\begin{cases}-\left[y-{(a+b)}\right]\quad&\text{for}\quad{a}<y<{a+b}\\\
\frac{\,b}{a}y\quad&\text{for}\quad\left\lvert y\right\rvert\leq{a}\\\
-\left[y+{(a+b)}\right]\quad&\text{for}\quad-{(a+b)}<y<-{a}.\end{cases}$
(3.38)

The meridional gradient $\mathrm{d}\theta/\mathrm{d}y$ is then a piecewise
constant $2(a+b)$-periodic function

$\frac{\mathrm{d}\theta}{\mathrm{d}y}=\Lambda\begin{cases}-1\quad&\text{for}\quad{a}<y<{a+b}\\\
\frac{\,b}{a}\quad&\text{for}\quad\left\lvert y\right\rvert\leq{a}\\\
-1\quad&\text{for}\quad-{(a+b)}<y<-{a}.\end{cases}$ (3.39)

Therefore the gradient in the frontal zones exceeds the gradient in the mixed
zones by a factor of $b/a$, which approaches infinity as
$b/a\rightarrow\infty$ in the sawtooth wave limit.

The zonal velocity, $u=-\partial_{y}\psi$, is obtained by using the inversion
relation (3.15) to solve for the streamfunction. Alternatively, taking the
meridional derivative of surface potential vorticity (3.3) gives

$\frac{\partial\theta}{\partial y}=\frac{1}{\sigma_{0}^{2}}\frac{\partial
u}{\partial z}\Big{|}_{z=0}.$ (3.40)

Then in Fourier space [$\partial_{y}\rightarrow\mathrm{i}k_{y}$ and
$\sigma_{0}^{-2}\partial_{z}|_{z=0}\rightarrow m(k)$] we obtain

$\hat{u}_{\bm{k}}=\frac{1}{m(k)}\left(\mathrm{i}\,k_{y}\,\hat{\theta}_{\bm{k}}\right),$
(3.41)

which shows that the induced zonal velocity is obtained by smoothing
$\mathrm{d}\theta/\mathrm{d}y$ by the function $m(k)$. An immediate
consequence is that the east-west asymmetry in the zonal velocity is
fundamentally due to the east-west asymmetry in the gradient
$\mathrm{d}\theta/\mathrm{d}y$.

Figure 3.2 shows an example of sloping sawtooth $\theta$ profile along with
the induced zonal velocities. For a power law inversion function,
$m(k)=m_{0}k^{\alpha}$, the parameter $\alpha$ modifies the zonal velocity in
two ways. First, in more local flows (with smaller $\alpha$), the zonal
velocity decays more rapidly away from the jet centre, as expected. Second,
the degree of smoothing increases with $\alpha$, and so more local regimes
(with smaller $\alpha$) are more east-west asymmetric, with the ratio
$\left\lvert u_{\mathrm{min}}\right\rvert/u_{\mathrm{max}}$ taking smaller
values for smaller $\alpha$. Figure 3.3(b) shows $\left\lvert
u_{\mathrm{min}}\right\rvert/u_{\mathrm{max}}$ as a function of $a/b$ for
$\alpha\in\\{1/2,\,1,\,3/2,\,2\\}$. For $\alpha=2$, we obtain $\left\lvert
u_{\mathrm{min}}\right\rvert/u_{\mathrm{max}}\rightarrow 1/2$ in the limit
$a/b\rightarrow 0$ so that eastward jets are only twice as strong as westward
flows in the perfect staircase limit (Danilov and Gurarie, 2004; Dritschel and
McIntyre, 2008). At $\alpha=3/2$, we find $\left\lvert
u_{\mathrm{min}}\right\rvert/u_{\mathrm{max}}\approx 0.29$ in the
$a/b\rightarrow 0$ limit so that eastward jets are now more than three time as
strong as westward flows. Once $\alpha\leq 1$, then the maximum jet velocity
diverges as $\alpha\rightarrow 0$ [figure 3.3(a)] and so $\left\lvert
u_{\mathrm{min}}\right\rvert/u_{\mathrm{max}}\rightarrow 0$ as $a/b\rightarrow
0$.

If $m(k)$ is not a power law, then the results are similar so long as $m(k)$
can be approximated by a power law at small wavenumbers. Figure 3.2 shows the
induced velocity for the inversion functions computed from idealized
stratifications profiles (shown in figure 3.1). Because these inversion
functions can be approximated by power laws $m(k)\approx k^{0.49}$ and
$m(k)\approx k^{1.50}$ at small wavenumbers, the induced velocity fields
nearly coincide with the velocity fields computed from power law inversion
functions with $\alpha=0.5$ and $\alpha=1.5$.

Figure 3.3: Properties of zonal velocity profiles induced by sloping sawtooth
profiles (3.38) of $\theta$ as a function of the non-dimensional frontal zone
width $a/b$ separating the mixed zones for four values of $\alpha$. Panel (a)
shows the maximum zonal velocity, panel (b) shows the ratio of westward speed
to eastward speed, panel (c) shows the rms zonal velocity, and panel (d) shows
the product $L_{j}k_{\mathrm{Rh}}$ where $L_{j}=a+b$ is the halfwidth
separation (the distance between $U_{\mathrm{mix}}$ and $U_{\mathrm{min}}$)
and $k_{\mathrm{Rh}}$ is the Rhines wavenumber (3.32).

Finally, we examine how the Rhines wavenumber, $k_{\mathrm{Rh}}$, relates to
jet spacing. Let

$L_{j}=a+b$ (3.42)

be the half-separation between the jets, i.e., the half distance between
consecutive zonal velocity maxima. For two-dimensional barotropic turbulence
(i.e., the $\alpha=2$ case), we have $L_{j}=45^{1/4}/k_{\mathrm{Rh}}\approx
2.59/k_{\mathrm{Rh}}$ in the staircase limit (i.e, for $a/b\rightarrow 0$,
Dritschel and McIntyre, 2008; Scott and Dritschel, 2012). This result is found
by solving for the zonal velocity induced by a staircase with halfwidth
$L_{j}=b$, taking the rms of the zonal velocity, and then substituting into
the definition of the generalized Rhines wavenumber (3.32). As figure 3.3(d)
shows, because the velocity field induced by a perfect staircase depends on
the inversion function, $m(k)$, the relationship between $L_{j}$ and
$k_{\mathrm{Rh}}$ also depends on the inversion function. For $m(k)=k^{3/2}$,
an analogous calculation gives $L_{j}\approx 2.35/k_{\mathrm{Rh}}$ in the
staircase limit. For $\alpha=1$, even though the maximum velocity diverges at
$a/b\rightarrow 0$, the rms velocity asymptotes to a constant value, and so we
obtain a half jet-separation of $L_{j}\approx 1.73/k_{\mathrm{Rh}}$ (figure
3.3). Finally in the $\alpha=1/2$ case, although the rms speed has not
converged by $a/b=10^{-6}$, the product $L_{j}\,k_{\mathrm{Rh}}$ is
approaching values close to zero.

###  3.4 Numerical Simulations

####  3.4.1 The numerical model

We use the pyqg pseudo-spectral model (Abernathey et al., 2019) which solves
the time-evolution equation (3.4) in a square domain with side length
$L=2\pi$. Time-stepping is through a third-order Adam-Bashforth scheme with
small-scale dissipation achieved through a scale-selective exponential filter
(Smith et al., 2002; Arbic and Flierl, 2003),

$\mathrm{ssd}=\begin{cases}1\quad\text{for}\quad k\leq k_{0}\\\
e^{-a(k-k_{0})^{4}}\quad\text{for}\quad k>k_{0},\end{cases}$ (3.43)

with $a=23.6$ and $k_{0}=0.65k_{\mathrm{Nyq}}$ where $k_{\mathrm{Nyq}}=\pi$ is
the Nyquist wavenumber. The forcing is isotropic, centred at wavenumber
$k_{f}=80$, and normalized so that the energy injection rate is
$\varepsilon=1$ (see appendix B in Smith et al., 2002). However, the effective
energy injection rate, $\varepsilon_{\mathrm{eff}}$, is smaller than
$\varepsilon$ due to dissipation. To determine $\varepsilon_{\mathrm{eff}}$
from numerical simulations, we use
$\varepsilon_{\mathrm{eff}}=2\,r\,\mathcal{E}$ where $\mathcal{E}$ is the
equilibrated total energy diagnosed from the model. In what follows, we report
values of $k_{\varepsilon}/k_{r}$ using $\varepsilon_{\mathrm{eff}}$ instead
of $\varepsilon$. The model is integrated forward in time until at least
$t=5/r$ to allow the fluid to reach equilibrium. All model runs use $1024^{2}$
horizontal grid points.

####  3.4.2 For what values of $k_{\varepsilon}/k_{r}$ do jets form?

Figure 3.4: Diagnostics from five series of simulations as a function of the
non-dimensional number $k_{\varepsilon}/k_{r}$. The first three series of
simulations have inversion function $m(k)=k^{\alpha}$ with
$\alpha\in\\{\nicefrac{{1}}{{2}},\,1,\,\nicefrac{{3}}{{2}}\\}$. For the other
two series, the inversion functions are shown in figure 3.8. Panel (a) shows
the ratio of energy in the zonal mode to total energy. Panel (b) shows the
ratio of domain averaged zonal speed to domain averaged meridional speed.
Panel (c) shows the ratio of westward zonal speed to eastward zonal speed.
Panel (d) shows the relationship between the halfwidth jet spacing, $L_{j}$,
and the Rhines wavenumber, $k_{\mathrm{Rh}}$.

For our first set of simulations, we vary $k_{\varepsilon}/k_{r}$ over the
values shown in figure 3.4. We do so by fixing $k_{r}=8$ and varying
$k_{\varepsilon}$. For a given value of $k_{\varepsilon}$, we choose $\Lambda$
and $r$ so as to maintain $k_{r}=8$ (the energy injection rate, $\varepsilon$,
is fixed at unity for all model runs). Given $k_{r}$ and $k_{\varepsilon}$, we
rearrange the definition of $k_{r}$ (3.35) to solve for
$\gamma=r\,\Lambda^{2}$,

$\gamma=m_{0}\,\varepsilon\,k_{r}^{\alpha+2},$ (3.44)

then solve for $r$ in the implicit equation (3.30) for $k_{\varepsilon}$ ,

$r=\frac{\gamma}{\left(\varepsilon\,k_{\varepsilon}\,[m(k_{\varepsilon})]^{2}\right)^{2/3}},$
(3.45)

and finally use the definition $\gamma=r\,\Lambda^{2}$ to solve for $\Lambda$.

##### Power law inversion functions

We first describe the results from three series of simulations with power law
inversion functions, $m(k)=k^{\alpha}$, with
$\alpha\in\\{\nicefrac{{1}}{{2}},\,1,\,\nicefrac{{3}}{{2}}\\}$. Summary
diagnostics from these simulations are shown in figure 3.4. In panel (a), we
observe that the ratio of energy in the zonal mode to total energy,
$\mathcal{E}_{\mathrm{zonal}}/\mathcal{E}$, increases with
$k_{\varepsilon}/k_{r}$, and that the majority of the total energy is in the
zonal mode for sufficiently large $k_{\varepsilon}/k_{r}$. For a fixed
$k_{\varepsilon}/k_{r}$, more of the total energy is zonal in more non-local
flows (with larger $\alpha$) than in more local flows (with smaller $\alpha$);
for $\alpha=\nicefrac{{3}}{{2}}$, we have
$\mathcal{E}_{\mathrm{zonal}}/\mathcal{E}\approx 1$ by
$k_{\varepsilon}/k_{r}\approx 6$ as compared to $k_{\varepsilon}/k_{r}\approx
12$ for $\alpha=1$. Moreover, for $\alpha=\nicefrac{{1}}{{2}}$, we find that
$\mathcal{E}_{\mathrm{zonal}}/\mathcal{E}$ asymptotes to approximately 0.9
once $k_{\varepsilon}/k_{r}\approx 18$ with little subsequent change for
larger values of $k_{\varepsilon}/k_{r}$. In panel (b), we observe a striking
contrast in the ratio $\overline{\left\lvert
u\right\rvert}/\overline{\left\lvert v\right\rvert}$ between different values
of $\alpha$ (the overline denotes a domain average). For
$\alpha=\nicefrac{{3}}{{2}}$, the domain averaged zonal speed,
$\overline{\left\lvert u\right\rvert}$, is approximately eight times larger
than the domain averaged meridional speed, $\overline{\left\lvert
v\right\rvert}$, for large $k_{\varepsilon}/k_{r}$. In contrast, for
$\alpha=\nicefrac{{1}}{{2}}$, $\overline{\left\lvert u\right\rvert}$ only
exceeds $\overline{\left\lvert v\right\rvert}$ by a multiple of two for large
$k_{\varepsilon}/k_{r}$.

Figure 3.5: Snapshots of the relative surface potential vorticity, $\theta$,
for simulations with power law inversion functions, $m(k)=k^{\alpha}$. In each
snapshot, the $\theta$ field is normalized by its maximum value in the
snapshot. Only one quarter of the domain is shown (i.e., $512^{2}$ grid
points).

Figure 3.6: The zonal mean total surface potential vorticity,
$\theta+\Lambda\,y$, in black and the zonal mean zonal velocity, $U$, in grey.

Next, we examine the jet structure for different $\alpha$ as a function of
$k_{\varepsilon}/k_{r}$. Figure 3.5 shows $\theta$-snapshots from model runs
with $m(k)=k^{\alpha}$. For each value of $\alpha$, two model runs are shown:
one where jets have just become visible in the $\theta$-snapshot and another
with the largest value of $k_{\varepsilon}/k_{r}$, which we expect to be
closest to the staircase limit. The jets are visible in these snapshots as the
regions with strong gradients. Because these are $\theta$-snapshots rather
than $(\theta+\Lambda\,y)$-snapshots, the $(\theta+\Lambda\,y)$-staircase is
instead a $\theta$-sawtooth, and the mixed zones between the jets are
approximately linear in $\theta$. We confirm this to be the case in figure
3.6, where the zonal averages of the total surface potential vorticity,
$\theta+\Lambda\,y$, and the zonal velocity are shown. For the
$\alpha=\nicefrac{{3}}{{2}}$ and $\alpha=1$ cases, we observe an approximate
staircase structure with nearly uniform mixed zones separated by frontal
zones, and with jets centred at sharp $\theta$ gradients. As expected from the
idealized staircases of section 3.3, close to the staircase limit, the
$\alpha=1$ jets are narrower than the $\alpha=\nicefrac{{3}}{{2}}$ jets, and
the ratio of maximum westward speed to maximum eastward speed,
$|U_{\mathrm{min}}|/|U_{\mathrm{max}}|$, is smaller at $\alpha=1$ than at
$\alpha=\nicefrac{{3}}{{2}}$.

In contrast to the $\alpha=\nicefrac{{3}}{{2}}$ and the $\alpha=1$ series, the
$\alpha=\nicefrac{{1}}{{2}}$ series approaches the staircase limit slowly with
$k_{\varepsilon}/k_{r}$. The $\alpha=\nicefrac{{1}}{{2}}$ staircase remains
smooth even at $k_{\varepsilon}/k_{r}=42$ [figure 3.6(c)]. The ratio of
frontal zone width to mixed zone width, $a/b$, is between $0.5$ and $0.65$ for
$\alpha=\nicefrac{{1}}{{2}}$ jets. In contrast, this ratio is between $0.15$
and $0.2$ for the $\alpha=\nicefrac{{3}}{{2}}$ and $\alpha=1$ jets. In part,
the broadness of the $\alpha=\nicefrac{{1}}{{2}}$ frontal zones is a
consequence of zonal averaging in the presence of large amplitude undulations.
However, it is evident from the $\theta$-snapshots of figure 3.5 that the
$\alpha=\nicefrac{{1}}{{2}}$ frontal zones are indeed broader than the
$\alpha=\nicefrac{{3}}{{2}}$ and $\alpha=1$ frontal zones [e.g., compare
panels (a) and (d) with (f) in figure 3.5], even without zonal averaging.

We now examine how the generalized Rhines wavenumber, $k_{\mathrm{Rh}}$,
relates to the jet spacing. From figure 3.3(d), a ratio of $a/b\approx 0.2$
leads to a $L_{j}\,k_{\mathrm{Rh}}\approx 2.2$ for
$\alpha=\nicefrac{{3}}{{2}}$ and $L_{j}\,k_{\mathrm{Rh}}\approx 2.0$ for
$\alpha=1$. But as figure 3.4(d) shows, we find values closer to
$L_{j}\,k_{\mathrm{Rh}}\approx 3$ for both of these cases. In contrast, for
the $\alpha=\nicefrac{{1}}{{2}}$ jets, figure 3.3(d) predicts $1.98\lesssim
L_{j}\,k_{\mathrm{Rh}}\lesssim 2.5$ for the observed range of $0.5\lesssim
a/b\lesssim 0.65$, but we find $L_{j}\,k_{\mathrm{Rh}}\approx 1.5$ for
$k_{\varepsilon}/k_{r}\geq 18$, which is smaller than predicted.

Returning to figure 3.5, we observe that there are undulations along the jets,
with smaller values of $\alpha$ corresponding to larger amplitude undulations.
These undulations propagate as waves and are less dispersive for smaller
$\alpha$, propagating eastward for $\alpha=\frac{3}{2}$, westward for
$\alpha=\nicefrac{{1}}{{2}}$, and are nearly stationary for $\alpha=1$.
Moreover, the waves in the $\alpha=\nicefrac{{1}}{{2}}$ case maintain their
shape as they propagate for a significant fraction of the domain, although
they eventually disperse or merge with other along jet waves. That we obtain
larger amplitude along jet undulations for smaller $\alpha$ is a consequence
of the more local inversion operator (3.15) at smaller $\alpha$. A jet in a
highly local flow (with small $\alpha$) is “a coherent structure that hangs
together strongly while being easy to push sideways” (McIntyre, 2008, in the
context of equivalent barotropic jets). However, although both an equivalent
barotropic jet and an $\alpha=\nicefrac{{1}}{{2}}$ jet exhibit large
meridional undulations, the undulations in the equivalent barotropic case are
frozen in place (because of a vanishing group velocity at large scales,
McIntyre, 2008) and so the equivalent barotropic jet behaves like a meandering
river with a fixed shape. In contrast, the $\alpha=\nicefrac{{1}}{{2}}$ jet
behaves like a flexible string whose shape evolves in time with the
propagation of weakly dispersive waves. Another difference between the two
cases is that an equivalent barotropic jet has a width given by the
deformation radius. In contrast, there is no analogous characteristic scale
for $\alpha=\nicefrac{{1}}{{2}}$ jets and, in principle, the jets should
become infinitely thin as $k_{\varepsilon}/k_{r}\rightarrow\infty$.

Figure 3.7: The total energy spectrum, $E(k)$, as a function of the
wavenumber, $k=k_{x}^{2}+k_{x}^{2}$, for three simulations with power law
inversion functions, $m(k)=k^{\alpha}$. The values of $k_{\varepsilon}/k_{r}$
are 18.0 for panel (a), 21.0 for panel (b), and 42.3 for panel (c).

Energy spectra for the three power law simulations are shown in figure 3.7.
The energy spectrum obtained from dimensional analysis (3.33) gives a
$k^{-\alpha-3}$ wavenumber dependence, which leads to the familiar $k^{-5}$
spectrum for beta-plane barotropic turbulence ($\alpha=2$). Although early
investigations (Chekhlov et al., 1996; Huang et al., 2000; Danilov and
Gryanik, 2004) found a $k^{-5}$ spectrum in barotropic $\beta$-plane
turbulence, Scott and Dritschel (2012) instead found a shallower $k^{-4}$
spectrum in the staircase limit (suggested earlier by Danilov and Gryanik,
2004; Danilov and Gurarie, 2004), which they explained as a consequence of the
sharp discontinuities of the staircase. Generalizing their argument to the
present case, a one dimensional $\theta(y)$ series with discontinuities
implies a Fourier series with coefficients decaying as $k^{-1}$, leading to a
$\theta^{2}$ spectrum of $k^{-2}$, and hence an energy spectrum

$E(k)\sim k^{-2}\left[m(k)\right]^{-1}.$ (3.46)

If $m(k)\sim k^{\alpha}$, then we obtain a spectrum $E(k)\sim k^{-\alpha-2}$,
which yields the $k^{-4}$ spectrum observed in Scott and Dritschel (2012),
where $\alpha=2$. For $\alpha=\nicefrac{{3}}{{2}}$, $\alpha=1$, and
$\alpha=\nicefrac{{1}}{{2}}$, the predicted spectrum is proportional to
$k^{-3.5}$, $k^{-3}$, and $k^{-2.5}$, respectively. The diagnosed spectra
shown in figure 3.7 are consistent with these shallow spectra, instead of
energy spectrum (3.33) obtained from dimensional considerations.

##### Inversion functions from $\sigma(z)$

Figure 3.8: Inversion functions [panel (a) and (c)] along with their
corresponding stratification profiles [panels (b) and (d), respectively]. The
stratification profiles are given by the piecewise function (3.18). For panel
(a), we have $\sigma_{0}=1$, $\sigma_{\mathrm{pyc}}=0.1$,
$h_{\mathrm{mix}}=0.01$, and $h_{\mathrm{lin}}=0.05$. For panel (c), we have
$\sigma_{0}=0.133$, $\sigma_{\mathrm{pyc}}=1$, $h_{\mathrm{mix}}=0.125$, and
$h_{\mathrm{lin}}=0.2$. The thin grey lines in panels (a) and (c) are given by
$k/\sigma_{0}$ and $k/\sigma_{\mathrm{pyc}}$.

We also ran two series of simulations where we specified a piecewise
stratification profile (3.18), and then obtained $m(k)$ by solving the
boundary value problem (3.10)–(3.12) at each wavenumber. The stratification
profiles and the resulting inversion functions are shown in figure 3.8. One
case consists of an increasing stratification profile [$\sigma^{\prime}(z)\geq
0$] with $\sigma_{0}=1$, $\sigma_{\mathrm{pyc}}=0.1$, $h_{\mathrm{mix}}=0.01$
and $h_{\mathrm{lin}}=0.05$. The resulting $m(k)$ is approximately linear for
$k\gtrsim 70$ and transitions to an approximate sub-linear wavenumber
dependence $m(k)\sim k^{0.40}$ for wavenumbers $5\lesssim k\lesssim 50$. The
second case consists of a decreasing stratification profile
[$\sigma^{\prime}(z)\leq 0$] with $\sigma_{0}=0.13$,
$\sigma_{\mathrm{pyc}}=1$, $h_{\mathrm{mix}}=0.125$ and
$h_{\mathrm{lin}}=0.2$. The resulting $m(k)$ is approximately linear at
wavenumbers $k\gtrsim 60$ and transitions to an approximate super linear
wavenumber dependence $m(k)\sim k^{1.50}$ between $3\lesssim k\lesssim 60$.

As seen in figure 3.4, the $\sigma^{\prime}(z)\leq 0$ case is similar to the
$\alpha=\nicefrac{{3}}{{2}}$ case, with the various diagnostics close to the
$\alpha=\nicefrac{{3}}{{2}}$ counterpart. In contrast, there are significant
differences between the $\sigma^{\prime}(z)\geq 0$ simulations and the
$\alpha=\nicefrac{{1}}{{2}}$ simulations. In the $\sigma^{\prime}(z)\geq 0$
series, the ratio of energy in the zonal mode to total energy continues to
increase as $k_{\varepsilon}/k_{r}$ is increased, whereas it asymptotes to a
constant in the $\alpha=\nicefrac{{1}}{{2}}$ series. Moreover, the ratio of
domain average zonal speed to domain averaged meridional speed,
$\overline{|u|}/\overline{|v|}$, is generally larger in the
$\sigma^{\prime}\geq 0$ series than in the $\alpha=\nicefrac{{1}}{{2}}$
series. Finally, for the largest values of $k_{\varepsilon}/k_{r}$, the
product $L_{j}\,k_{\mathrm{Rh}}$ reaches smaller values in the
$\sigma^{\prime}\geq 0$ simulations than in the $\alpha=\nicefrac{{1}}{{2}}$
simulations.

Figure 3.9: Snapshots of the relative surface potential vorticity, $\theta$,
normalized by its maximum value in the snapshot, for simulations with
inversion functions shown in figure 3.8. Only one quarter of the domain is
shown (i.e., $512^{2}$ grid points).

These differences can be explained by the snapshots of figure 3.9 as well as
the zonal averages of figure 3.6. As expected from the model diagnostics, both
the snapshots and the zonal average from the $\sigma^{\prime}\leq 0$
simulation are qualitatively similar to the $\alpha=\nicefrac{{3}}{{2}}$
simulation. In contrast, the $\sigma^{\prime}\geq 0$ snapshot is evidently
closer to the staircase limit than the $\alpha=\nicefrac{{1}}{{2}}$ snapshot:
the mixed zones are more homogeneous and the frontal zones are sharper. The
zonal average of the $\sigma^{\prime}\geq 0$ simulation in figure 3.6 also
shows how the $\sigma^{\prime}\geq 0$ simulation is closer to the staircase
limit than the $\alpha=\nicefrac{{1}}{{2}}$ simulation, although, again, zonal
averaging in the presence of large amplitude undulations is artificially
smoothing the jets. Therefore, the differences in the diagnostics between the
$\sigma^{\prime}\geq 0$ series and the $\alpha=\nicefrac{{1}}{{2}}$ series
stem from the more rapid approach (i.e., at smaller $k_{\varepsilon}/k_{r}$)
of the $\sigma^{\prime}\geq 0$ series to the staircase limit.

####  3.4.3 Simulations with fixed parameters

The dependence of the non-dimensional number $k_{\varepsilon}/k_{r}$ on the
external parameters $\varepsilon,\Lambda$, and $r$ depends on the functional
form of $m(k)$. For example, if $m(k)\sim k^{\alpha}$, then

$k_{\varepsilon}/{k_{r}}=|\Lambda|^{\frac{4-\alpha}{(2\alpha+1)(\alpha+2)}}\,\varepsilon^{\frac{\alpha-1}{(1+2\alpha)(\alpha+2)}}\,r^{\frac{-1}{\alpha+2}}.$
(3.47)

Because the forcing intensity wavenumber, $k_{\varepsilon}$, is obtained by
solving the implicit equation for $k_{\varepsilon}$ (3.30), an analogous
expression for $k_{\varepsilon}/k_{r}$ is not possible for general $m(k)$.
However, at sufficiently large $k_{\varepsilon}$, the inversion function
asymptotes to $m(k_{\varepsilon})\approx k_{\varepsilon}/\sigma_{0}$ and so,
using $\alpha$-turbulence expression for $k_{\varepsilon}$ (3.31) with
$\alpha=1$, we obtain

$k_{\varepsilon}/k_{r}\approx|\Lambda|^{\frac{\alpha}{\alpha+2}}\,\varepsilon^{\frac{1-\alpha}{3\alpha+6}}\,r^{\frac{-1}{\alpha+2}}\,m_{0}^{\frac{1}{\alpha+2}}\,\sigma_{0}^{2/3}$
(3.48)

for large $k_{\varepsilon}$, where $\alpha$ is the approximate power law
dependence of $m(k)$ near $k_{r}$.

Therefore, simulations with identical $k_{\varepsilon}/k_{r}$ but distinct
inversion functions cannot be directly compared because they have different
values of $\Lambda$ and $r$. Here, we investigate how the stratification
modifies jet structure as all other parameters are held fixed. We therefore
run two additional series of simulations with the stratification profiles and
inversion functions shown in figure 3.1. The stratification profiles were
chosen so that they both have identical stratification at the upper boundary.
One case corresponds to an increasing stratification profile,
$\sigma^{\prime}\geq 0$, with an approximate power law dependence of $m(k)\sim
k^{0.49}$ at small wavenumbers. The second case consists of a decreasing
stratification profile, $\sigma^{\prime}\leq 0$, with a $m(k)\sim k^{1.50}$ at
small wavenumbers. Aside from the different stratification profiles, these two
series of simulations are run under the same conditions as the constant
stratification ($\alpha=1$) simulations of section 3.4.2, with identical
values of $\Lambda$, $\varepsilon$, and $r$.

Figure 3.10: As in figure 3.4, but the $\sigma^{\prime}\geq 0$ and
$\sigma^{\prime}\leq 0$ series now only differ from the $\sigma^{\prime}=0$
(i.e., $\alpha=1$) series only in the vertical stratification (and hence the
inversion function).

Summary diagnostics are shown in figure 3.10. We see that, at a fixed value of
$\Lambda$ and $r$, more of the total energy is in the zonal mode in the
$\sigma^{\prime}(z)\leq 0$ simulation than in the constant stratification
simulation, which in turn is larger than the $\sigma^{\prime}(z)\geq 0$
simulation (and similarly for the ratio of area averaged zonal to meridional
speeds, $\overline{\left\lvert u\right\rvert}/\overline{\left\lvert
v\right\rvert}$). Therefore, increased non-locality (larger $\alpha$) promotes
anisotropy in the velocity field and leads to larger zonal velocities relative
to meridional velocities. Indeed, figure 3.11 shows $\theta$ snapshots from
these simulations; the more local, $\sigma^{\prime}\geq 0$, simulations have
larger meridional undulations along the jets. Moreover, compared to the
$k_{\varepsilon}/k_{r}=15$ constant stratification simulation in figure
3.5(c), the $\sigma^{\prime}(z)\leq 0$ simulation in figure 3.11(a) is closer
to the staircase limit whereas the frontal zones in the
$\sigma^{\prime}(z)\geq 0$ simulation [figure 3.11(c)] remain broad. Finally,
we show values of the product $L_{j}\,k_{\mathrm{Rh}}$, relating the Rhines
wavenumber to the half spacing between the jets, in figure 3.10(c). These
values are similar to those in shown in figure 3.4(c).

Figure 3.11: Snapshots of relative surface potential vorticity, $\theta$,
where $\theta$ is normalized by its maximum value in the snapshot. Panels (a)
and (c) are from simulations with identical $\Lambda$, $r$, and $\varepsilon$
as the $\alpha=1$ simulation shown in figure 3.5(c), whereas panels (b) and
(d) are from simulations with identical $\Lambda$, $r$, and $\varepsilon$ as
the $\alpha=1$ simulation shown in figure 3.5(d). Only one quarter of the
domain is shown (i.e., $512^{2}$ grid points).

###  3.5 Conclusion

We have examined the emergence of staircase-like buoyancy structures in
surface quasigeostrophic turbulence with a mean background buoyancy gradient.
We found that the stratification’s vertical structure controls the locality of
the inversion operator and the dispersion of surface-trapped Rossby waves. As
we go from decreasing stratification profiles [$\sigma^{\prime}(z)\leq 0$] to
increasing stratification profiles [$\sigma^{\prime}(z)\geq 0$], the inversion
operator becomes more local and Rossby wave less dispersive. In all cases, we
find that the non-dimensional ratio, $k_{\varepsilon}/k_{r}$, controls the
extent of inhomogeneous buoyancy mixing. Larger $k_{\varepsilon}/k_{r}$
correspond to sharper buoyancy gradients at jet centres with larger peak jet
velocities that are separated by more homogeneous mixed-zones. Moreover, we
found that the staircase limit is reached at smaller $k_{\varepsilon}/k_{r}$
in more non-local flows; the staircase limit is reached by
$k_{\varepsilon}/k_{r}\approx 15$ for our $\sigma\leq 0$ simulations compared
to $k_{\varepsilon}/k_{r}\approx 25$ for our $\sigma\geq 0$ simulations.

In addition, once the staircase limit is reached, the dynamics of the jets
depends on the locality of the inversion operator and, hence, on the
stratification’s vertical structure. In flows with a more non-local inversion
operator [or decreasing stratification, $\sigma^{\prime}(z)\leq 0$], we obtain
straight jets that are perturbed by dispersive, eastward propagating, along
jet waves. In contrast, for more local flows [or over increasing
stratification, $\sigma^{\prime}(z)\geq 0$], we obtain jets with latitudinal
meanders on the order of the jet spacing. The shape of these jets evolves in
time as these meanders propagate westwards as weakly dispersive waves.

The inversion operator’s locality is also reflected in two more aspects of the
dynamics. First, the domain-averaged zonal speed exceeds the domain-averaged
meridional speed by approximately a factor of eight in our most non-local
simulations, whereas this ratio is merely two in our most local simulations.
This observation is consistent with the fact that jets are narrower and
exhibit larger latitudinal meanders in more local flows. Second, for a given
Rhines wavenumber, jets in more local flows are closer together. Indeed, we
found $L_{j}\,k_{\mathrm{Rh}}\approx 3-4$ in our most non-local simulations,
where $L_{j}$ is the jet half spacing, as compared to
$L_{j}\,k_{\mathrm{Rh}}\approx 0.5-1.5$ in our most local simulations.

Several open questions remain. First, we have not examined the dynamics of the
along jet waves. As we observed, these waves propagate eastwards in our most
non-local simulations [with $\sigma^{\prime}(z)\leq 0$] but westwards for our
most local simulations [with $\sigma^{\prime}(z)\geq 0$]. These waves are not
described by the dispersion relation (3.23); rather, the relevant model is
that of freely propagating edge waves along a buoyancy discontinuity
(McIntyre, 2008). However, the difficulty here is that a jump discontinuity in
the buoyancy field results in infinite velocities over constant or increasing
stratification. In addition, the relationship of the along jet waves in the
staircase limit to the non-linear zonons found by Sukoriansky et al. (2008)
remains unclear.

The divergence of the velocity at a buoyancy discontinuities raises a second
question. Is there a limit to how close the staircase limit can be approached?
In barotropic dynamics, the velocity remains finite at a jump continuity in
the vorticity, and, in this case, Scott and Dritschel (2012) report that a
vorticity staircase case can be approached arbitrarily. Whether this result
continues to hold for arbitrarily sharp buoyancy gradients and arbitrarily
large zonal velocities is not clear. Because the rms velocity seems to
converge for arbitrarily sharp staircases, even for the most local inversion
relations we considered, there may not be any energetic reason precluding
arbitrarily sharp buoyancy gradients.

Finally, there remains the question of how relevant these results are for the
upper ocean, which, in addition to surface buoyancy gradients, has interior
potential vorticity gradients as well. In particular, our neglect of the
$\beta$-effect limits the direct relevance of this model to the upper ocean.
Whether surface buoyancy staircases can emerge under more realistic oceanic
conditions requires further investigation.

##  Chapter 4 Normal Modes With Boundary Dynamics in Geophysical Fluids

Abstract

  * Three-dimensional geophysical fluids support both internal and boundary-trapped waves. To obtain the normal modes in such fluids we must solve a differential eigenvalue problem for the vertical structure (for simplicity, we only consider horizontally periodic domains). If the boundaries are dynamically inert (e.g., rigid boundaries in the Boussinesq internal wave problem, flat boundaries in the quasigeostrophic Rossby wave problem) the resulting eigenvalue problem typically has a Sturm-Liouville form and the properties of such problems are well-known. However, when restoring forces are also present at the boundaries, then the equations of motion contain a time-derivative in the boundary conditions and this leads to an eigenvalue problem where the eigenvalue correspondingly appears in the boundary conditions. In certain cases, the eigenvalue problem can be formulated as an eigenvalue problem in the Hilbert space $L^{2}\oplus\mathbb{C}$ and this theory is well-developed. Less explored is the case when the eigenvalue problem takes place in a Pontryagin space, as in the Rossby wave problem over sloping topography. This article develops the theory of such problems and explores the properties of wave problems with dynamically-active boundaries. The theory allows us to solve the initial value problem for quasigeostrophic Rossby waves in a region with sloping bottom (we also apply the theory to two Boussinesq problems with a free-surface). For a step-function perturbation at a dynamically-active boundary, we find that the resulting time-evolution consists of waves present in proportion to their projection onto the dynamically-active boundary.

###  4.1 Introduction

An important tool in the study of wave motion near a stable equilibrium is the
separation of variables. When applicable, this elementary technique transforms
a linear partial differential equation into an ordinary differential
eigenvalue problem for each coordinate (e.g., Hillen et al., 2012). Upon
solving the differential eigenvalue problems, one obtains the normal modes of
the physical system. The normal modes are the fundamental wave motions for the
given restoring forces, each mode represents an independent degree of freedom
in which the physical system can oscillate, and any solution of the wave
problem may be written as a linear combination of these normal modes.

To derive the normal modes, we must first linearize the dynamical equations of
motion about some equilibrium state. We then encounter linearized restoring
forces of two kinds:

  * _1._

volume-permeating forces experienced by fluid particles in the interior, and

  * _2._

boundary-confined forces only experienced by fluid particles at the boundary.

Examples of volume-permeating forces include the restoring forces resulting
from continuous density stratification and continuous volume potential
vorticity gradients. These restoring forces respectively result in internal
gravity waves (Sutherland, 2010) and Rossby waves (Vallis, 2017). Examples of
boundary-confined restoring forces include the gravitational force at a free-
surface (i.e., at a jump discontinuity in the background density), forces
arising from gradients in surface potential vorticity (Schneider et al.,
2003), and the molecular forces giving rise to surface tension. These
restoring forces respectively result in surface gravity waves (Sutherland,
2010), topographic/thermal waves (Hoskins et al., 1985), and capillary waves
(Lamb, 1975).

In the absence of boundary-confined restoring forces, we can often apply
Sturm-Liouville theory (e.g., Hillen et al., 2012; Zettl, 2010) to the
resulting eigenvalue problem. We thus obtain a countable infinity of waves
whose vertical structures form a basis of $L^{2}$, the space of square-
integrable functions (see section 4.2), and, given some initial vertical
structure, we know how to solve for the subsequent time-evolution as a linear
combination for linearly independent waves. Moreover, a classic result of
Sturm-Liouville theory is that the $n$th mode has $n$ internal zeros.

In the presence of boundary-confined restoring forces, the governing equations
have a time-derivative in the boundary conditions. The resulting eigenvalue
problem correspondingly contains the eigenvalue parameter in the boundary
conditions. Sturm-Liouville theory is inapplicable to such problems.

In this chapter, we present a general method for solving these problems by
delineating a generalization of Sturm-Liouville theory. Some consequences of
this theory are the following. There is a countable infinity of waves whose
vertical structures form a basis of $L^{2}\oplus\mathbb{C}^{s}$, where $s$ is
the number of dynamically-active boundaries; thus, each boundary-trapped wave,
in mathematically rigorous sense, provides an additional degree of freedom to
the problem. The modes satisfy an orthogonality relation involving boundary
terms, the modes may have a negative norm, and the modes may have finite jump
discontinuities at dynamically-active boundaries (although the _solutions_ are
always continuous, see section 4.3.3). When negative norms are possible (as in
quasigeostrophic theory), there is a new expression for the Fourier
coefficients that one must use to solve initial value problems [see equation
(4.23)]. We can also expand boundary step-functions (representing some
boundary localized perturbation) as a sum of modes. Moreover, the $n$th mode
may not have $n$ internal zeros; indeed, depending on physical parameters in
the problem, two or three linearly independent modes with an identical number
of internal zeros may be present.

We also show that the eigenfunction expansion of a function is term-by-term
differentiable, with the derivative series converging uniformly on the whole
interval, regardless of the boundary condition the function satisfies at the
dynamically-active boundaries. This property is in contrast with a traditional
Sturm-Liouville eigenfunction expansion where the term-by-term derivative
converges uniformly only if the function satisfies the same boundary condition
as the eigenfunctions.

We apply the theory to three geophysical wave problems. The first is that of a
Boussinesq fluid with a free-surface; we find that the $n$th mode has $n$
internal zeros. The second example is that of a rotating Boussinesq fluid with
a free-surface where we assume that the stratification suppresses rotational
effects in the interior but not at the upper boundary. We find that there are
two linearly independent modes with $M$ internal zeros, where the integer $M$
depends on the ratio of the Coriolis parameter to the horizontal wavenumber,
and that the eigenfunctions have a finite jump discontinuity at the upper
boundary. The third application is to a quasigeostrophic fluid with a sloping
lower boundary. We find that modes with an eastward phase speed have a
negative norm whereas modes with a westward phase speed have a positive norm
(the sign of the norm has implications for the relative phase of a wave and
for series expansions). Moreover, depending on the propagation direction,
there can be two linearly independent modes with no internal zeros. For all
three examples, we outline the properties of the resulting series expansions
and provide the general solution. We also consider the time-evolution
resulting from a vertically localized perturbation at a dynamically-active
boundary; we idealize such a perturbation as a boundary step-function. The
step-function perturbation induces a time-evolution in which the amplitude of
each constituent wave is proportional to the projection of that wave onto the
boundary.

To our knowledge, most of the above results cannot be found in the literature
[however, the gravity wave orthogonality relation has been noted before, e.g.,
Gill (1982) and Kelly (2016) for the hydrostatic case and Olbers (1986) and
Early et al. (2020) for the non-hydrostatic case]. For instance, we provide
the only solution to the initial value problem for Rossby waves over
topography in the literature [equation (4.78)]. Moreover, many of the
properties we discuss arise in practical problems in physical oceanography.
The number of internal zeros of Rossby waves is also a useful quantity in
observational physical oceanography [e.g., Clément et al. (2014) and de La
Lama et al. (2016)]. In addition, the question of whether the quasigeostrophic
baroclinic modes are complete is a controversial one. Lapeyre (2009) has
suggested that the baroclinic modes are incomplete because they assume a
vanishing surface buoyancy anomaly. Consequently, Smith and Vanneste (2012)
address this issue by deriving an $L^{2}\oplus\mathbb{C}^{2}$ basis for
quasigeostrophic theory. Yet many authors, citing completeness theorems from
Sturm-Liouville theory, insist that the baroclinic modes are indeed complete
and can represent all quasigeostrophic states (Ferrari and Wunsch, 2010;
LaCasce, 2012; Rocha et al., 2015). This chapter shows that, by including
boundary-confined restoring forces, we obtain a set of modes with additional
degrees-of-freedom. These degrees-of-freedom manifest in the behaviour of
eigenfunction expansions at the boundaries. In addition, the distinction
between $L^{2}$ and $L^{2}\oplus\mathbb{C}^{s}$ bases that we present here is
useful for equilibrium statistical mechanical calculations where one must
decompose fluid motion onto a complete set of modes (Bouchet and Venaille,
2012; Venaille et al., 2012).

The plan of the chapter is the following. We formulate the mathematical theory
in section 4.2. We then apply the theory to the two Boussinesq wave problems,
in section 4.3, and to the quasigeostrophic wave problem, in section 4.4. We
consider the time-evolution of a localized perturbation at a dynamically-
active boundary in section 4.5. We then conclude in section 4.6.

###  4.2 The eigenvalue problem

In this section, we outline the theory of the differential eigenvalue problem,

$\displaystyle-(p\,\phi^{\prime})^{\prime}+q\,\phi$
$\displaystyle=\lambda\,r\,\phi\quad\text{for }\quad
z\in\left(z_{1},z_{2}\right)$ (4.1)
$\displaystyle-\left[a_{1}\,\phi(z_{1})-b_{1}\,(p\,\phi^{\prime})(z_{1})\right]$
$\displaystyle=\lambda\left[c_{1}\,\phi(z_{1})-d_{1}\,(p\,\phi^{\prime})(z_{1})\right]$
(4.2)
$\displaystyle-\left[a_{2}\,\phi(z_{2})-b_{2}\,(p\,\phi^{\prime})(z_{2})\right]$
$\displaystyle=\lambda\left[c_{2}\,\phi(z_{2})-d_{2}\,(p\,\phi^{\prime})(z_{2})\right],$
(4.3)

where $p^{-1},q$, and $r$ are real-valued integrable functions;
$a_{i},b_{i},c_{i}$, and $d_{i}$ are real numbers with $i\in\\{1,2\\}$; and
where $\lambda\in\mathbb{C}$ is the eigenvalue parameter. We further assume
that $p>0$ and $r>0$, that $p$ and $r$ are twice continuously differentiable,
that $q$ is continuous, and that $(a_{i},b_{i})\neq(0,0)$ for $i\in\\{1,2\\}$.
The system of equations (4.1)–(4.3) is an eigenvalue problem for the
eigenvalue $\lambda\in\mathbb{C}$ and differs from a regular Sturm-Liouville
problem in that $\lambda$ appears in the boundary conditions (4.2) and (4.3).
That is, setting $c_{i}=d_{i}=0$ recovers the traditional Sturm-Liouville
problem. The presence of $\lambda$ as part of the boundary condition leads to
some fundamentally new mathematical features that are the subject of this
section and fundamental to the physics of this chapter.

It is useful to define the two boundary parameters

$D_{i}=(-1)^{i+1}\left(a_{i}\,d_{i}-b_{i}\,c_{i}\right)\quad i=1,2.$ (4.4)

Just as the function $r$ acts as a weight for the interval $(z_{1},z_{2})$ in
traditional Sturm-Liouville problems, the constants $D_{i}^{-1}$ will play
analogous roles for the boundaries $z=z_{i}$ when $D_{i}\neq 0$.

##### Outline of the mathematics

The right-definite case, when the $D_{i}\geq 0$ for $i\in\\{1,2\\}$, is well-
known in the mathematics literature; most of the right-definite results in
this section are due to Evans (1970), Walter (1973), and Fulton (1977). In
contrast, the left-definite case, defined below, is much less studied. In this
section, we generalize the right-definite results of Fulton (1977) to the
left-definite problem as well as provide an intuitive formulation (in terms of
functions rather than vectors, for a vector formulation see Fulton, 1977) of
the eigenvalue problem.

In section 4.2.1 we state the conditions under which we obtain real
eigenvalues and a basis of eigenfunctions. We proceed, in section 4.2.2, to
explore the properties of eigenfunctions and eigenfunction expansions.
Finally, in section 4.2.3, we discuss oscillation properties of the
eigenfunctions. Additional properties of the eigenvalue problem are found in
appendix 4.A and a literature review, along with various technical proofs, is
found in appendix 4.B.

####  4.2.1 Formulation of the problem

##### The function space of the problem

We denote by $L^{2}$ the Hilbert space of square-integrable “functions” $\phi$
on the interval $(z_{1},z_{2})$ satisfying

$\int_{z_{1}}^{z_{2}}\left\lvert\phi\right\rvert^{2}\,r\,\mathrm{d}z<\infty.$
(4.5)

To be more precise, the elements of $L^{2}$ are not functions but rather
equivalence classes of functions (e.g., Reed and Simon, 1980, section I.3).
Two functions, $\phi$ and $\psi$, are equivalent in $L^{2}$ (i.e., $\phi=\psi$
in $L^{2}$) if they agree in a mean-square sense on $[z_{1},z_{2}]$,

$\int_{z_{1}}^{z_{2}}\left\lvert\phi(z)-\psi(z)\right\rvert^{2}\,r\,\mathrm{d}z=0.$
(4.6)

Significantly, we can have $\phi=\psi$ in $L^{2}$ but $\phi\neq\psi$
pointwise.

Furthermore, as a Hilbert space, $L^{2}$ is endowed with a positive-definite
inner product

$\left<\phi,\psi\right>_{\sigma}=\int_{z_{1}}^{z_{2}}{\phi}^{\ast}\,\psi\,\mathrm{d}\sigma=\int_{z_{1}}^{z_{2}}{\phi}^{\ast}\,\psi\
r\,\mathrm{d}z,$ (4.7)

where the symbol ${\ }^{\ast}$ denotes complex conjugation and the measure
$\sigma$ associated $L^{2}$ induces a differential element
$\mathrm{d}\sigma=r\,\mathrm{d}z$ (see appendix 4.A). The positive-
definiteness is ensured by our assumption that $r>0$ (i.e.,
$\left<\phi,\phi\right>_{\sigma}>0$ for $\phi\neq 0$ when $r>0$).

It is well-known that traditional Sturm-Liouville problems [i.e., equations
(4.1)–(4.3) with $c_{i}=d_{i}=0$ for $i=1,2$] are eigenvalue problems in some
subspace of $L^{2}$ (Debnath and Mikusinski, 2005). For the more general case
of interest here, the eigenvalue problem occurs over a “larger” function space
denoted by $L^{2}_{\mu}$ which we construct in appendix 4.A.

Let the integer $s\in\\{0,1,2\\}$ denote the number of $\lambda$-dependent
boundary conditions and let $S$ denote the set

$S=\\{j\ |\ j\in\\{1,2\\}\text{ and }(c_{j},d_{j})\neq(0,0)\\}.$ (4.8)

$S$ is one of $\emptyset,\\{1\\},\\{2\\},\\{1,2\\}$ and $s$ is the number of
elements in the set $S$. In appendix 4.A, we show that $L^{2}_{\mu}$ is
isomorphic to the space $L^{2}\oplus\mathbb{C}^{s}$ and is thus “larger” than
$L^{2}$ by $s$ dimensions.

We denote elements of $L^{2}_{\mu}$ by upper case letters $\Psi$; we define
$\Psi(z)$ for $z\in[z_{1},z_{2}]$ by

$\Psi(z)=\begin{cases}\Psi(z_{i})\quad&\textrm{at }z=z_{i},\textrm{ for
}\,i\in S,\\\ \psi(z)\quad&\textrm{otherwise},\end{cases}$ (4.9)

where $\Psi(z_{i})\in\mathbb{C}$ are constants, for $i\in S$, and the
corresponding lower case letter $\psi$ denotes an element of $L^{2}$. Two
elements $\Phi$ and $\Psi$ of $L^{2}_{\mu}$ are equivalent in $L^{2}_{\mu}$ if
and only if

  * _1._

$\Phi(z_{i})=\Psi(z_{i})$ for $i\in S$, and

  * _2._

$\phi(z)$ and $\psi(z)$ are equivalent in $L^{2}$ [i.e., as in equation
(4.6)].

Here, $\Phi$, as an element of $L^{2}_{\mu}$, is defined as in equation (4.9).
The primary difference between $L^{2}$ and $L^{2}_{\mu}$ is that $L^{2}_{\mu}$
discriminates between functions that disagree at $\lambda$-dependent
boundaries.

The measure $\mu$ associated with $L^{2}_{\mu}$ (see appendix 4.A) induces a
differential element

$\mathrm{d}\mu(z)=\left[r(z)+\sum_{i\in
S}D_{i}^{-1}\,\delta(z-z_{i})\right]\mathrm{d}z,$ (4.10)

where $\delta(z)$ is the Dirac delta. The induced inner product on
$L^{2}_{\mu}$ is

$\left<\Phi,\Psi\right>=\int_{z_{1}}^{z_{2}}{\Phi}^{\ast}\,{\Psi}\,\mathrm{d}\mu=\int_{z_{1}}^{z_{2}}{\Phi}^{\ast}\,{\Psi}\,r\,\mathrm{d}z+\sum_{i\in
S}D_{i}^{-1}\,{\Phi(z_{i})}^{\ast}\,\Psi(z_{i}).$ (4.11)

If $D_{i}>0$ for $i\in S$ then this inner product is positive-definite and
$L^{2}_{\mu}$ is a Hilbert space. However, this is not the case in general.

Let $\kappa$ denote the number of negative $D_{i}$ for $i\in S$ (the possible
values are $\kappa=0,1,2$). Then $L^{2}_{\mu}$ has a $\kappa$-dimensional
subspace of elements $\Psi$ satisfying

$\left<\Psi,\Psi\right><0.$ (4.12)

This makes $L^{2}_{\mu}$ a Pontryagin space of index $\kappa$ (Bognár, 1974).
If $\kappa=0$ then $L^{2}_{\mu}$ is again a Hilbert space. In the present
case, $L^{2}_{\mu}$ also has an infinite-dimensional subspace of elements
$\psi$ satisfying

$\left<\Psi,\Psi\right>>0.$ (4.13)

##### Reality and completeness

In appendix 4.A.2, we reformulate the eigenvalue problem (4.1)–(4.3) as an
eigenvalue problem of the form,

$\mathcal{L}\,\Phi=\lambda\,\Phi,$ (4.14)

in a subspace of $L^{2}_{\mu}$, where $\mathcal{L}$ is a linear operator and
$\Phi$ an element of $L^{2}_{\mu}$. We also define the notions of right- and
left-definiteness that are required for the reality and completeness theorem
below. The following two propositions can be considered to define right- and
left-definiteness for applications of the theory. Both propositions are
obtained through straightforward manipulations (see appendix 4.A).

######  Proposition 4.2.1 (Criterion for right-definiteness).

The eigenvalue problem (4.1)–(4.3) is right-definite if $r>0$ and $D_{i}>0$
for $i\in S$.

######  Proposition 4.2.2 (Criterion for left-definiteness).

The eigenvalue problem (4.1)–(4.3) is left-definite if the following
conditions hold:

  * (i)

the functions $p,q$ satisfy $p>0,q\geq 0$,

  * (ii)

for the $\lambda$-dependent boundary conditions, we have

$\frac{a_{i}\,c_{i}}{D_{i}}\leq 0,\quad\frac{b_{i}\,d_{i}}{D_{i}}\leq
0,\quad(-1)^{i}\frac{a_{i}\,d_{i}}{D_{i}}\geq 0\quad\text{for }i\in S.$ (4.15)

  * (iii)

for the $\lambda$-independent boundary conditions, we have

$b_{i}=0\quad\text{or}\quad(-1)^{i+1}\frac{a_{i}}{b_{i}}\geq 0\quad\text{ if
}b_{i}\neq 0\quad\text{ for }i\in\\{1,2\\}\setminus S.$ (4.16)

The notions of right and left-definiteness are not mutually exclusive. Namely,
a problem can be neither right- or left-definite; both right- and left-
definite; only right-definite; or only left-definite. In this chapter, we
always assume that $p>0$ and $r>0$.

The reality of the eigenvalues and the completeness of the eigenfunctions in
the space $L^{2}_{\mu}$ is given by the following theorem.

######  Theorem 4.2.3 (Reality and completeness).

Suppose the eigenvalue problem (4.1)–(4.3) is either right-definite or left-
definite. Moreover, if the problem is not right-definite, we assume that
$\lambda=0$ is not an eigenvalue. Then the eigenvalue problem (4.1)–(4.3) has
a countable infinity of real simple eigenvalues $\lambda_{n}$ satisfying

$\lambda_{0}<\lambda_{1}<\dots<\lambda_{n}<\dots\rightarrow\infty,$ (4.17)

with corresponding eigenfunctions $\Phi_{n}$. Furthermore, the set of
eigenfunctions $\\{\Phi_{n}\\}_{n=0}^{\infty}$ is a complete orthonormal basis
for $L^{2}_{\mu}$ satisfying

$\left<\Phi_{m},\Phi_{n}\right>=\pm\delta_{mn}.$ (4.18)

###### Proof.

See appendix 4.B.3. ∎

Recall that $\kappa$ denotes the number of negative $D_{i}$ for $i\in S$. We
then have the following corollary of the proof of theorem 4.2.3.

######  Proposition 4.2.4.

Suppose the eigenvalue problem (4.1)–(4.3) is left-definite and that
$\lambda=0$ is not an eigenvalue. Then there are $\kappa$ negative eigenvalues
and their eigenfunctions satisfy

$\left<\Phi,\Phi\right><0.$ (4.19)

The remaining eigenvalues are positive and their eigenfunctions satisfy

$\left<\Phi,\Phi\right>>0.$ (4.20)

In other words, proposition 4.2.4 states that we have the relationship

$\lambda_{n}\left<\Phi_{n},\Phi_{n}\right>>0$ (4.21)

for left-definite problems.

####  4.2.2 Properties of the eigenfunctions

For the remainder of section 4.2, we assume that the eigenvalue problem
(4.1)–(4.3) satisfies the requirements of theorem 4.2.3.

##### Eigenfunction expansions

The eigenvalue problem (4.1)–(4.3) has _eigenfunctions_
$\\{\Phi_{n}\\}_{n=0}^{\infty}$ as well as corresponding _solutions_
$\\{\phi_{n}\\}_{n=0}^{\infty}$. In other words, while the $\phi_{n}$ are the
solutions to the differential equation defined by equations (4.1)–(4.3) with
$\lambda=\lambda_{n}$, the eigenfunctions required by the operator formulation
of the problem [equation (4.14)] are $\Phi_{n}$. The functions $\Phi_{n}$ and
$\phi_{n}$ are related by equation (4.9), with the boundary values
$\Phi_{n}(z_{i})$ of $\Phi_{n}$ determined by

$\Phi_{n}(z_{i})=\left[c_{i}\,\phi(z)-d_{i}\,(p\,\phi^{\prime})(z)\right]\quad\text{for
}i\in S.$ (4.22)

Thus, while the solutions $\phi_{n}$ are continuously differentiable over the
closed interval $[z_{1},z_{2}]$, the eigenfunctions $\Phi_{n}$ are
continuously differentiable over the open interval $(z_{1},z_{2})$ but
generally have finite jump discontinuities at the $\lambda$-dependent
boundaries. The eigenfunctions $\Phi_{n}$ are continuous in the closed
interval $[z_{1},z_{2}]$ only if $c_{i}=1$ and $d_{i}=0$ for $i\in S$. In this
case, the eigenfunctions $\Phi_{n}$ coincide with the solutions $\phi_{n}$ on
the closed interval $[z_{1},z_{2}]$.

The boundary conditions of the eigenvalue problem (4.1)–(4.3) are not unique.
One can multiply each boundary condition by an arbitrary constant to obtain an
equivalent problem. To uniquely specify the eigenfunctions in physical
applications, the boundary coefficients $\\{a_{i},b_{i},c_{i},d_{i}\\}$ of
equations (4.1)–(4.3) must be chosen so that $r\,\mathrm{d}z$ has the same
dimensions as $D_{i}^{-1}\,\delta(z-z_{i})\,\mathrm{d}z$ [recall that
$\delta(z)$ has the dimension of inverse length]. In the quasigeostrophic
problem, we must also invoke continuity and set $c_{i}=1$.

Since $\\{\Phi_{n}\\}_{n=0}^{\infty}$ is a basis for $L^{2}_{\mu}$, then any
$\Psi\in L^{2}_{\mu}$ may be expanded in terms of the eigenfunctions (Bognár,
1974, thereom IV.3.4),

$\Psi=\sum_{n=0}^{\infty}\frac{\left<\Psi,\Phi_{n}\right>}{\left<\Phi_{n},\Phi_{n}\right>}\,\Phi_{n}.$
(4.23)

We emphasize that the above equality is an equality in $L^{2}_{\mu}$ and not a
pointwise equality [see the discussion following equation (4.9)]. Some
properties of $L^{2}_{\mu}$ expansions are given in appendix 4.A.3.

An important property that distinguishes the basis
$\\{\Phi_{n}\\}_{n=0}^{\infty}$ of $L^{2}_{\mu}$ from an $L^{2}$ basis is its
“sensitivity” to function values at boundary points $z=z_{i}$ for $i\in S$.
See section 4.5 for a physical application.

A natural question is whether the basis $\\{\Phi_{n}\\}_{n=0}^{\infty}$ of
$L^{2}_{\mu}$ is also a basis of $L^{2}$. Recall that the set
$\\{\Phi_{n}\\}_{n=0}^{\infty}$ is a basis of $L^{2}$ if every element
$\psi\in L^{2}$ can be written _uniquely_ in terms of the functions
$\\{\Phi_{n}\\}_{n=0}^{\infty}$. However, in general, this is not true. If
$s>0$, the $L^{2}_{\mu}$ basis $\\{\Phi_{n}\\}_{n=0}^{\infty}$ is overcomplete
in $L^{2}$ (Walter, 1973; Russakovskii, 1997).

##### Uniform convergence and term-by-term differentiability

Along with the eigenfunction expansion (4.23) in terms of the eigenfunctions
$\\{\Phi_{n}\\}_{n=0}^{\infty}$, we also have the expansion

$\sum_{n=0}^{\infty}\frac{\left<\Psi,\Phi_{n}\right>}{\left<\Phi_{n},\Phi_{n}\right>}\,\phi_{n}$
(4.24)

in terms of the solutions $\phi_{n}$. The two expansions differ in their
behaviour at $\lambda$-dependent boundaries, $z=z_{i}$ for $i\in S$, but are
otherwise equal. In particular, the $\Phi_{n}$ eigenfunction expansion (4.23)
must converge to $\Psi(z_{i})$ at $z=z_{i}$ for $i\in S$ as this equality is
required for $\Psi$ to be equal to the series expansion (4.23) in
$L^{2}_{\mu}$ [see the discussion following equation (4.9)]. Some properties
of both expansions are given in appendix 4.A.4. In particular, theorem 4.A.4
shows that the $\phi_{n}$ solution series (4.24) does not generally converge
to $\Psi(z_{i})$ at $z=z_{i}$.

Figure 4.1: Convergence to a function $F(z)=1+2z+(3/2)\sin(2\pi
z)\cos(\pi^{2}z^{2}+3)$ for $z\in[-1,0]$, shown in panel (a), by various
eigenfunction expansions of $-\phi^{\prime\prime}=\lambda\,\phi$ with fifteen
terms, as discussed in section 4.2.2. Panel (b) shows the Fourier sine
expansion of $F$. Since the sine eigenfunctions vanish at the boundaries
$z=-1,0$, the series expansion will not converge to $F$ at the boundaries.
Panel (c) shows the cosine expansion of $F$ which converges uniformly to $F$
on the closed interval $[-1,0]$. Panel (d) shows an expansion with boundary
coefficients in equations (4.2)–(4.3) given by
$(a_{1},b_{1},c_{1},d_{1})=(-0.5,-5,1,0)$ and
$(a_{2},b_{2},c_{2},d_{2})=(0.5,-5,1,0)$. Since the $c_{i}=1$ and $d_{i}=0$,
then $\Phi_{n}=\phi_{n}$ and the series expansions (4.23) and (4.24) coincide.
As with the cosine series, the expansion converges uniformly to $F$ on
$[-1,0]$. The derivative of $F$ is shown in panel (e). Panel (f) show the
derivative of the sine series expansion. In panel (g), we show the
differentiated cosine series which does not converge to the derivative
$F^{\prime}$ at the boundaries $z=z_{1},z_{2}$. In contrast, in panel (h), the
differentiated series obtained from a problem with $\lambda$-dependent
boundary conditions converges uniformly to the derivative $F^{\prime}$.

The following theorem is of central concern for physical applications.

######  Theorem 4.2.5 (Uniform convergence).

Let $\psi$ be a twice continuously differentiable function on $[z_{1},z_{2}]$
satisfying all $\lambda$-independent boundary conditions of the eigenvalue
problem (4.1)–(4.3). Define the function $\Psi$ on $z\in[z_{1},z_{2}]$ by

$\Psi(z)=\begin{cases}c_{i}\,\psi(z)-d_{i}\,(p\,\psi^{\prime})(z)\quad&\textrm{at
}z=z_{i},\textrm{ for }\,i\in S,\\\
\psi(z)\quad&\textrm{otherwise}.\end{cases}$ (4.25)

Then

$\psi(z)=\sum_{n=0}^{\infty}\frac{\left<\Psi,\Phi_{n}\right>}{\left<\Phi_{n},\Phi_{n}\right>}\,\phi_{n}(z)\quad\textrm{and}\quad\psi^{\prime}(z)=\sum_{n=0}^{\infty}\frac{\left<\Psi,\Phi_{n}\right>}{\left<\Phi_{n},\Phi_{n}\right>}\,\phi_{n}^{\prime}(z)$
(4.26)

with both series converging uniformly and absolutely on $[z_{1},z_{2}]$.

###### Proof.

See appendix 4.B.4. ∎

If $c_{i}=1$ and $d_{i}=0$ for $i\in S$ then we can replace $\Phi_{n}$ by
$\phi_{n}$ and $\Psi$ by $\psi$ in equation (4.26).

In addition, if both boundary conditions of the eigenvalue problem (4.1)–(4.3)
are $\lambda$-dependent, then both expansions in equation (4.26) converge
uniformly on $[z_{1},z_{2}]$ regardless of the boundary conditions $\psi$
satisfies. As discussed in appendix 4.A.4, for traditional Sturm-Liouville
expansions, an analogous result holds only if $\psi$ satisfies the same
boundary conditions as the eigenfunctions. Figure 4.1 contrasts the
convergence behaviour of such a problem (with continuous eigenfunctions, so
$c_{i}=1$ and $d_{i}=0$ for $i\in S$) with the convergence behaviour of sine
and cosine series. All numerical solutions in this chapter are obtained using
a pseudo-spectral code in Dedalus (Burns et al., 2020).

Figure 4.2: Convergence to a function $\tilde{F}$ with finite jump
discontinuities at the boundaries by two eigenfunction expansions (with
$\lambda$-dependent boundary conditions) of
$-\phi^{\prime\prime}=\lambda\,\phi$ with fifteen terms, as discussed in
section 4.2.2. The function $\tilde{F}(z)$ is defined by $\tilde{F}(z)=F(z)$
for $z\in(z_{1},z_{2})$ where $F(z)$ is the function defined in figure 4.1,
$F(-1)=0.5$ at the lower boundary, and $F(0)=-0.7$ at the upper boundary. The
function $\tilde{F}$ is shown in panel (a). In panel (b), the boundary
coefficients in equations (4.2)–(4.3) are given by
$(a_{1},b_{1},c_{1},d_{1})=(-0.5,-5,1,0)$ and
$(a_{2},b_{2},c_{2},d_{2})=(0.5,-5,1,0)$ as in figure 4.1. In panel (c), the
boundary coefficients are $(a_{1},b_{1},c_{1},d_{1})=(-0.5,-5,1,0.1)$ and
$(a_{2},b_{2},c_{2},d_{2})=(0.5,-5,1,-0.1)$. The $\Phi_{n}$ expansion (4.23)
and the $\phi_{n}$ expansion (4.24) are not generally equal at the boundaries
$z=-1,0$; this figure shows the $\Phi_{n}$ expansion. The $\Phi_{n}$ series
(4.23) converges pointwise to $\tilde{F}$ on $[-1,0]$, however, the
convergence will not be uniform if $d_{i}=0$ for $i\in S$, as in panel (b).
The boundary values of the $\Phi_{n}$ series (4.23) are shown with a black
dot. In panel (b), the eigenfunctions $\Phi_{n}$ are continuous and a large
number of terms are required for the series to converge to the discontinuous
function $\tilde{F}$. Panel (c) shows that the discontinuous eigenfunction
$\Phi_{n}$ have almost converged to the $\tilde{F}$—including at the jump
discontinuities; the black dot in panel (c) overlap with the grey dots, which
represent the boundary values of $\tilde{F}$. Although the $\phi_{n}$ series
(4.24) converges to $\tilde{F}$ in the interior $(-1,0)$, the $\phi_{n}$
series does not generally converge to $\tilde{F}$ at the boundaries but
instead converges to the values given in theorem 4.A.4.

Another novel property of the eigenfunction expansions is that we obtain
pointwise convergence to functions that are smooth in the interior of the
interval, $(z_{1},z_{2})$, but have finite jump discontinuities at
$\lambda$-dependent boundaries (see appendix 4.A.4). If $d_{i}\neq 0$ for
$i\in S$, the convergence is even uniform (Fulton, 1977, corollary 2.1).
Figure 4.2 illustrates the convergence behaviour for eigenfunction expansions
with $\lambda$-dependent boundary conditions in the two cases $d_{i}=0$ and
$d_{i}\neq 0$. Note the presence of Gibbs-like oscillations in the case
$d_{i}=0$ shown in panel (b). Although the $\Phi_{n}$ eigenfunction series
(4.23) converges pointwise to the discontinuous function, the $\phi_{n}$
solution series (4.24) converges to the values given in theorem 4.A.4 at the
$\lambda$-dependent boundaries. The ability of these series expansions to
converge to functions with boundary jump discontinuities is related to their
ability to expand distributions in the Bretherton (1966) “$\delta$-function
formulation” of a problem.

####  4.2.3 Oscillation theory

Recall that for regular Sturm-Liouville problems [i.e., equations (4.1)–(4.3)
with $c_{i}=d_{i}=0$] we obtain a countable infinity of real simple
eigenvalues, $\lambda_{n}$, that may be ordered as

$\lambda_{0}<\lambda_{1}<\lambda_{2}<\dots\rightarrow\infty,$ (4.27)

with associated eigenfunctions $\phi_{n}$. The $n$th eigenfunction $\phi_{n}$
has $n$ internal zeros in the interval $(z_{1},z_{2})$ so that no two
eigenfunctions have the same number of internal zeros.

However, once the eigenvalue $\lambda$ appears in the boundary conditions,
there may be up to $s+1$ linearly independent eigenfunctions with the same
number of internal zeros. The crucial parameters deciding the number of zeros
is $-b_{i}/d_{i}$ for $i\in S$, where $b_{i}$ and $d_{i}$ are the boundary
coefficients appearing in the boundary conditions (4.2)–(4.3). The following
lemma outlines the possibilities when only one boundary condition is
$\lambda$-dependent.

######  Lemma 4.2.6 (Location of double oscillation count).

Suppose that $s=1$, $i\in S$, and let $\kappa$ be the number of negative
$D_{i}$ for the eigenvalue problem (4.1)–(4.3). We have the following
possibilities.

  * (i)

Right-definite, $d_{i}\neq 0$: The eigenfunction $\Phi_{n}$ corresponding to
the eigenvalue $\lambda_{n}$ has $n$ internal zeros if
$\lambda_{n}<-b_{i}/d_{i}$ and $n-1$ internal zero if
$-b_{i}/d_{i}\leq\lambda_{n}$.

  * (ii)

Right-definite, $d_{i}=0$: The $n$th eigenfunction has $n$ internal zeros.

  * (iii)

Left-definite: If $\kappa=0$ then all eigenvalues are positive, the problem is
right-definite, and either (i) or (ii) applies. Otherwise, if $\kappa=1$, then
the eigenvalues may be ordered as

$\lambda_{0}<0<\lambda_{1}<\lambda_{2}<\dots\rightarrow\infty.$ (4.28)

Both eigenfunctions $\Phi_{0}$ and $\Phi_{1}$ have no internal zeros. The
remaining eigenfunctions $\Phi_{n}$, for $n>1$, have $n-1$ internal zeros.

###### Proof.

Parts (i), (ii) and (ii) are due to Linden (1991), Binding et al. (1994), and
Binding and Browne (1999), respectively. ∎

When both boundary conditions are $\lambda$-dependent, the situation is
similar. See Binding et al. (1994) and Binding and Browne (1999) for further
discussion.

###  4.3 Boussinesq gravity-capillary waves

Consider a rotating Boussinesq fluid on an $f$-plane with a reference
Boussinesq density of $\rho_{0}$. The fluid is subject to a constant
gravitational acceleration $g$ in the downwards, $-\hat{\bm{z}}$, direction,
and to a surface tension $T$ (with dimensions of force per unit length, see
Lamb, 1975) at its upper boundary. The upper boundary of the fluid, given by
$z=\eta$, is a free-surface defined by the function $\eta(\bm{x},t)$, where
$\bm{x}=\hat{\bm{x}}\,x+\hat{\bm{y}}\,y$ is the horizontal position vector.
The lower boundary of the fluid is a flat rigid surface given by $z=-H$. The
fluid region is periodic in both horizontal directions $\hat{\bm{x}}$ and
$\hat{\bm{y}}$.

####  4.3.1 Linear equations of motion

The governing equations for infinitesimal perturbations about a background
state of no motion, characterized by a prescribed background density of
$\rho_{B}=\rho_{B}(z)$, are

$\displaystyle\partial_{t}^{2}\nabla^{2}w+f_{0}^{2}\,\partial_{z}^{2}w+N^{2}\,\nabla^{2}_{z}w=0\quad$
$\displaystyle\text{for }z\in\left(-H,0\right)$ (4.29) $\displaystyle
w=0\quad$ $\displaystyle\text{for }z=-H$ (4.30)
$\displaystyle-\partial_{t}^{2}\partial_{z}w-f_{0}^{2}\,\partial_{z}w+g_{b}\,\nabla^{2}_{z}w-\tau\,\nabla^{4}_{z}w=0\quad$
$\displaystyle\text{for }z=0,$ (4.31)

where $w$ is the vertical velocity, $f_{0}$ is the constant value of the
Coriolis frequency, the prescribed buoyancy frequency $N^{2}$ is given by

$N^{2}(z)=-\frac{g}{\rho_{0}}\frac{\mathrm{d}\rho_{B}(z)}{\mathrm{d}z},$
(4.32)

the acceleration $g_{b}$ is the effective gravitational acceleration at the
upper boundary

$g_{b}=-\frac{g}{\rho_{0}}\left[\rho_{a}-\rho_{B}(0-)\right]$ (4.33)

where $\rho_{a}$ is the density of the overlying fluid, and the parameter
$\tau$ is given by

$\tau=\frac{T}{\rho_{0}}$ (4.34)

where $T$ is the surface tension. The three-dimensional Laplacian is denoted
$\nabla^{2}=\partial_{x}^{2}+\partial_{y}^{2}+\partial_{z}^{2}$, the
horizontal Laplacian is denoted by
$\nabla^{2}_{z}=\partial_{x}^{2}+\partial_{y}^{2}$, and the horizontal
biharmonic operator is given by
$\nabla^{4}_{z}=\nabla^{2}_{z}\,\nabla^{2}_{z}$. See equation (1.37) in
Dingemans (1997) for the surface tension term in (4.31). The remaining terms
in equation (4.29)–(4.31) are standard (Gill, 1982). Consistent with our
assumption that $\eta(\bm{x},t)$ is small, we evaluate the upper boundary
condition at $z=0$ in equation (4.31).

####  4.3.2 Non-rotating Boussinesq fluid

Figure 4.3: The vertical velocity eigenfunctions $\hat{W}_{n}=\hat{w}_{n}$ of
the non-rotating Boussinesq eigenvalue problem (4.36)–(4.38) for two distinct
wavenumbers with constant stratification, as discussed in section 4.3.2. For
both wavenumbers, the $n$th eigenfunction has $n$ internal zeros as in regular
Sturm-Liouville theory. The zeroth mode ($n=0$) corresponds to a surface
gravity wave and is trapped to the upper boundary for large horizontal
wavenumbers. In contrast to the internal wave problem with a rigid lid, the
modes $\hat{w}_{n}$ now depend on the horizontal wavenumber $k$ through the
boundary condition (4.38), however, this dependence is weak for $n\gg 1$, as
can be observed in this figure; for $n>2$, the modes for $k=0.01$ (in black)
and for $k=10$ (in grey) nearly coincide. The horizontal wavenumbers $k$ are
non-dimensionalized by $H$.

We assume wave solutions of the form

$w(\bm{x},z,t)=\hat{w}(z)\,\mathrm{e}^{\mathrm{i}\left(\bm{k}\cdot\bm{x}-\omega
t\right)}$ (4.35)

where $\bm{k}=\hat{\bm{x}}\,k_{x}+\hat{\bm{y}}\,k_{y}$ is the horizontal
wavevector and $\omega$ is the angular frequency. Substituting the wave
solution (4.35) into equations (4.29)–(4.31) and setting $f_{0}=0$ yields

$\displaystyle-\hat{w}^{\prime\prime}+k^{2}\,\hat{w}=\sigma^{-2}\,N^{2}\,\hat{w}\quad$
$\displaystyle\text{for }z\in(-H,0)$ (4.36) $\displaystyle\hat{w}=0\quad$
$\displaystyle\text{for }z=-H$ (4.37)
$\displaystyle(g_{b}+\tau\,k^{2})^{-1}\hat{w}^{\prime}=\sigma^{-2}\,\hat{w}\quad$
$\displaystyle\text{for }z=0,$ (4.38)

where $\sigma=\omega/k$ is the phase speed and
$k=\left\lvert\bm{k}\right\rvert$ is the horizontal wavenumber. Equations
(4.36)–(4.38) are an eigenvalue problem for the eigenvalue
$\lambda=\sigma^{-2}$.

##### Definiteness & the underlying function space

Equations (4.36)–(4.38) form an eigenvalue problem with one
$\lambda$-dependent boundary condition, namely, the upper boundary condition
(4.38). The underlying function space is then

$L^{2}_{\mu}\cong L^{2}\oplus\mathbb{C}.$ (4.39)

We write $\hat{W}_{n}$ for the eigenfunctions and $\hat{w}_{n}$ for the
solutions of the eigenvalue problem (4.36)–(4.38) [see the paragraph
containing equation (4.22)]. The eigenfunctions $\hat{W}_{n}$ are related to
the solutions $\hat{w}_{n}$ by equation (4.9) with boundary values
$\hat{W}_{n}(0)$ given by equation (4.22). However, since $c_{2}=1$ and
$d_{2}=0$ in equation (4.38) [compare with equations (4.1)–(4.3)] then
$\hat{W}_{n}=\hat{w}_{n}$ on the closed interval $[-H,0]$; thus, the solutions
$w_{n}$ are also the eigenfunctions.

By theorem 4.2.3, the eigenfunctions $\\{\hat{w}_{n}\\}_{n=0}^{\infty}$ form
an orthonormal basis of $L^{2}_{\mu}$. For functions $\varphi$ and $\phi$, the
inner product is

$\displaystyle\left<\varphi,\phi\right>$
$\displaystyle=\frac{1}{N_{0}^{2}\,H}\left[\int_{-H}^{0}\varphi\,\phi\,N^{2}\,\textrm{d}z+(g_{b}+\tau\,k^{2})\varphi(0)\,\phi(0)\right]$
(4.40)

obtained from equations (4.11) and equation (4.4); we have introduced the
factor $1/(N^{2}_{0}\,H)$ in the above expression for dimensional consistency
in eigenfunction expansions ($N^{2}_{0}$ is a typical value of $N^{2}$).
Orthonormality is then given by

$\displaystyle\delta_{mn}$
$\displaystyle=\left<\hat{w}_{m},\hat{w}_{n}\right>$ (4.41)

and we have chosen the solutions $\hat{w}_{n}$ to be non-dimensional (so the
Kronecker delta is non-dimensional as well).

One verifies that the eigenvalue problem (4.36)–(4.38) is right-definite using
proposition 4.2.1 and left-definite using proposition 4.2.2. Right-
definiteness implies that $L^{2}_{\mu}$, with the inner product (4.40), is a
Hilbert space. That is, all eigenfunctions $\hat{w}_{n}$ satisfy

$\left<\hat{w}_{n},\hat{w}_{n}\right>>0.$ (4.42)

Left-definiteness, along with proposition 4.2.4, ensures that all eigenvalues
$\lambda_{n}=\sigma_{n}^{-2}$ are positive. Indeed, the phase speeds
$\sigma_{n}$ satisfy

$\sigma_{0}^{2}>\sigma_{1}^{2}>\dots>\sigma^{2}_{n}>\dots\rightarrow 0.$
(4.43)

##### Properties of the eigenfunctions

By lemma 4.2.6, the $n$th eigenfunction $\hat{w}_{n}$ has $n$ internal zeros
in the interval $(-H,0)$. See figure 4.3 for an illustration of the first six
eigenfunctions.

The eigenfunctions $\\{\hat{w}_{n}\\}_{n=0}^{\infty}$ are complete in $L^{2}$
but do not form a basis in $L^{2}$; in fact, the basis is overcomplete in
$L^{2}$. The presence of a free-surface provides an additional degree of
freedom over the usual rigid-lid $L^{2}$ basis of internal wave
eigenfunctions. Indeed, the $n=0$ wave in figure 4.3 corresponds to a surface
gravity wave, while the remaining modes are internal gravity waves (with some
surface motion).

##### Expansion properties

Given a twice continuously differentiable function $\chi(z)$ satisfying
$\chi(-H)=0$, then, from theorem 4.2.5, we have

$\chi(z)=\sum_{n=0}^{\infty}\left<\chi,\hat{w}_{n}\right>\,\hat{w}_{n}(z)\quad\textrm{and}\quad\chi^{\prime}(z)=\sum_{n=0}^{\infty}\left<\chi,\hat{w}_{n}\right>\,w_{n}^{\prime}(z),$
(4.44)

with both series converging uniformly on $[-H,0]$ (note that $\chi$ is not
required to satisfy any particular boundary condition at $z=0$). If $\chi$ is
the vertical structure at time $t=0$ (and at some wavevector $\bm{k}$) and we
assume $\partial_{t}w(\bm{x},z,t=0)=0$, then the subsequent time-evolution is
given by

$w(\bm{x},z,t)=\sum_{n=0}^{\infty}\left<\chi,\hat{w}_{n}\right>\,w_{n}(z)\,\cos\left(\sigma_{n}kt\right)\mathrm{e}^{\mathrm{i}\bm{k}\cdot\bm{x}}.$
(4.45)

##### The $f$-plane hydrostatic problem

Suppose we have hydrostatic gravity waves on an $f$-plane with free surface at
the upper boundary, as in Kelly (2016). The appropriate inner product is
obtained by setting $\tau=0$ in the inner product (4.40). All the above
results on the eigenfunctions of gravity-capillary waves carry over to the
hydrostatic $f$-plane problem provided we set

$\sigma^{2}=\frac{\omega^{2}-f_{0}^{2}}{k^{2}}.$ (4.46)

####  4.3.3 A Boussinesq fluid with a rotating upper boundary

Figure 4.4: The vertical velocity eigenfunctions $\hat{W}_{n}$ of a Boussinesq
fluid with a rotating upper boundary—eigenvalue problem (4.48)–(4.50). This
figure is discussed in section 4.3.3. The wavenumbers $k$ in the figure are
non-dimensionalized by the depth $H$. The dots represent the values of the
eigenfunctions at the boundaries. Note that the eigenfunctions have a finite
jump discontinuity at $z=0$. For $k\,H=0.01$ (given by the black line) there
are two modes with no internal zeros. As $k$ increases, we obtain two modes
with one internal zero (at $k\,H=0.05$, the thick grey line) and then two
modes with three internal zeros (at $k\,H=0.11$, the thin grey line).

Although this next example is not geophysically relevant, it has the curious
property that the resulting eigenfunctions are discontinuous.

Let $N^{2}_{0}$ be a typical value of $N^{2}(z)$. Consider the situation where
$f_{0}^{2}/N_{0}^{2}\ll 1$ but

$\frac{g_{b}+\tau\,k^{2}}{f_{0}^{2}\,H}\sim O(1).$ (4.47)

Accordingly, we may neglect the Coriolis parameter in the interior equation
(4.29) but not at the upper boundary condition (4.31). Substituting the wave
solution (4.35) into equations (4.29)–(4.31) yields

$\displaystyle-\hat{w}^{\prime\prime}+k^{2}\,\hat{w}=\sigma^{-2}\,N^{2}\,\hat{w}\quad$
$\displaystyle\text{for }z\in(-H,0)$ (4.48) $\displaystyle\hat{w}=0\quad$
$\displaystyle\text{for }z=-H$ (4.49)
$\displaystyle(g_{b}+\tau\,k^{2})^{-1}\hat{w}^{\prime}=\sigma^{-2}\left[\hat{w}+\frac{f_{0}^{2}}{k^{2}}(g_{b}+\tau\,k^{2})^{-1}\hat{w}^{\prime}\right]\quad$
$\displaystyle\text{for }z=0,$ (4.50)

where $\sigma=\omega/k$ is the phase speed. Equations (4.48)–(4.50) form an
eigenvalue problem for the eigenvalue $\lambda=\sigma^{-2}$.

##### Definiteness & the underlying function space

As in the previous case, the eigenvalue problem is both right-definite and
left-definite, the underlying function space $L^{2}_{\mu}$ is given by
equation (4.39), and the appropriate inner product is equation (4.40). By
right-definiteness, the space $L^{2}_{\mu}$, equipped with the inner product
(4.40), is a Hilbert space; thus, all eigenfunctions $\hat{W}_{n}$ satisfy

$\left<\hat{W}_{m},\hat{W}_{n}\right>>0.$ (4.51)

By theorem 4.2.3, all eigenvalues $\lambda_{n}=\sigma_{n}^{-2}$ are real and
the corresponding eigenfunctions $\\{\hat{W}_{n}\\}_{n=0}^{\infty}$ form an
orthonormal basis of the Hilbert space $L^{2}_{\mu}$. By proposition 4.2.4,
all eigenvalues $\lambda_{n}=\sigma_{n}^{-2}$ are positive and satisfy
equation (4.43).

##### Boundary jump discontinuity of the eigenfunctions

The main difference between the previous non-rotating problem (4.36)–(4.38)
and the above problem (4.48)–(4.50) is that, in the present problem, if
$f_{0}\neq 0$ then $d_{2}\neq 0$ [see equation (4.3)]. Thus, by equation
(4.22), the eigenfunctions $\hat{W}_{n}$ generally have a jump discontinuity
at the upper boundary $z=0$ (see figure 4.4) and so are not equal to the
solutions $\hat{w}_{n}$. The eigenfunctions $\hat{W}_{n}$ are defined by
$\hat{W}_{n}(z)=\hat{w}_{n}(z)$ for $z\in[-H,0)$ and

$\hat{W}_{n}(0)=\hat{w}_{n}(0)+\frac{f_{0}^{2}}{k^{2}}(g_{b}+\tau\,k^{2})^{-1}\,\hat{w}_{n}^{\prime}(0)$
(4.52)

[see equation (4.22)]. It is not difficult to show that

$\hat{W}_{n}(0)\approx 0\quad\text{ for }n\text{ sufficiently large,}$ (4.53)

as can be seen in figure 4.4.

Physical motion is given by the solutions $\hat{w}_{n}$ which are continuous
over the closed interval $[-H,0]$. The jump discontinuity in the
eigenfunctions $\hat{W}_{n}$ does not correspond to any physical motion;
instead, the eigenfunctions $\hat{W}_{n}$ are convenient mathematical aids
used to obtain eigenfunction expansions in the function space $L^{2}_{\mu}$.

##### Number of internal zeros of the eigenfunctions

Another consequence of $d_{2}\neq 0$ is that by, lemma 4.2.6, there are two
distinct solutions $\hat{w}_{M}$ and $\hat{w}_{M+1}$ with the same number of
internal zeros (i.e., $M$) in the interval $(-H,0)$. Noting that

$-\frac{b_{2}}{d_{2}}=\frac{k^{2}}{f_{0}^{2}}$ (4.54)

the integer $M$ is determined by

$\sigma_{0}^{2}>\sigma_{1}^{2}>\dots>\sigma_{M}^{2}>\frac{f_{0}^{2}}{k^{2}}\geq\sigma^{2}_{M+1}>\dots>0.$
(4.55)

A smaller $f_{0}$ or a larger $k$ implies a larger $M$ and hence that
$\hat{w}_{M}$ and $\hat{w}_{M+1}$ have a larger number of internal zeros, as
shown in figure 4.4.

##### Expansion properties

As in the previous problem, the eigenfunctions are complete in $L^{2}_{\mu}$
but overcomplete in $L^{2}$ due to the additional surface gravity-capillary
wave.

Given a twice continuously differentiable function $\chi(z)$ satisfying
$\chi(-H)=0$, we define the discontinuous function $X(z)$ by

$X(z)=\begin{cases}\chi(z)&\quad\textrm{for }z\in[-H,0)\\\
\chi(0)+\frac{f_{0}^{2}}{k^{2}}\left(g_{b}+\tau\,k^{2}\right)^{-1}\chi^{\prime}(0)&\quad\textrm{for
}z=0\end{cases}$ (4.56)

as in theorem 4.2.5. Then, by theorem 4.2.5, we have the expansions

$\chi(z)=\sum_{n=0}^{\infty}\left<X,\hat{W}_{n}\right>\,\hat{w}_{n}(z)\quad\textrm{and}\quad\chi^{\prime}(z)=\sum_{n=0}^{\infty}\left<X,\hat{W}_{n}\right>\,w_{n}^{\prime}(z).$
(4.57)

Moreover, if $\chi(z)$ is the vertical structure at $t=0$ (and at some
wavevector $\bm{k}$) and we assume $\partial_{t}w(\bm{x},z,t=0)=0$, then the
subsequent time-evolution is given by

$w(\bm{x},z,t)=\sum_{n=0}^{\infty}\left<X,\hat{W}_{n}\right>\,\hat{w}_{n}(z)\,\cos\left(\sigma_{n}kt\right)\,\mathrm{e}^{\mathrm{i}\bm{k}\cdot\bm{x}}.$
(4.58)

###  4.4 Quasigeostrophic waves

####  4.4.1 Linear equations

Linearizing the quasigeostrophic equations about a quiescent background state
with an infinitesimally sloping lower boundary, at $z=-H$, and a rigid flat
upper boundary, at $z=0$, renders

$\displaystyle\partial_{t}\left[\nabla^{2}_{z}\psi+\partial_{z}\left(S^{-1}\,\partial_{z}\psi\right)\right]+\hat{\bm{z}}\cdot\left(\bm{\nabla}_{z}\psi\times\bm{\nabla}_{z}f\right)$
$\displaystyle=0\quad\text{for }z\in(-H,0)$ (4.59)
$\displaystyle\partial_{t}\left(S^{-1}\,\partial_{z}\psi\right)+\hat{\bm{z}}\cdot\left(\bm{\nabla}_{z}\psi\times
f_{0}\,\bm{\nabla}_{z}h\right)$ $\displaystyle=0\quad\text{for }z=-H$ (4.60)
$\displaystyle\partial_{t}\left(S^{-1}\,\partial_{z}\psi\right)$
$\displaystyle=0\quad\text{for }z=0.$ (4.61)

See Rhines (1970), Charney and Flierl (1981), Straub (1994) for details. The
streamfunction $\psi$ is defined through
$\bm{u}=\hat{\bm{z}}\times\bm{\nabla}_{z}\psi$ where $\bm{u}$ is the
horizontal velocity and
$\bm{\nabla}_{z}=\hat{\bm{x}}\,\partial_{x}+\hat{\bm{y}}\,\partial_{y}$ is the
horizontal Laplacian. The stratification parameter $S$ is given by

$S(z)=\frac{N^{2}(z)}{f_{0}^{2}},$ (4.62)

where $N^{2}$ is the buoyancy frequency and $f_{0}$ is the reference Coriolis
parameter. The latitude dependent Coriolis parameter $f$ is defined by

$f(y)=f_{0}+\beta\,y.$ (4.63)

Finally, $h(\bm{x})$ is the height of the topography at the lower boundary and
is a linear function of the horizontal position vector $\bm{x}$. Consistent
with quasigeostrophic theory, we assume that topography $h$ is small and so we
evaluate the lower boundary condition at $z=-H$ in equation (4.60).

####  4.4.2 The streamfunction eigenvalue problem

We assume wave solutions of the form

$\psi(\bm{x},z,t)=\hat{\psi}(z)\,\mathrm{e}^{\mathrm{i}(\bm{k}\cdot\bm{x}-\omega
t)}$ (4.64)

where $\bm{k}=\hat{\bm{x}}\,k_{x}+\hat{\bm{y}}\,k_{y}$ is the horizontal
wavevector and $\omega$ is the angular frequency.

We denote by $\Delta\theta_{f}$ the angle between the horizontal wavevector
$\bm{k}$ and the gradient of Coriolis parameter $\bm{\nabla}_{z}f$,

$\sin{(\Delta\theta_{f})}=\frac{1}{k\,\beta}\,\hat{\bm{z}}\cdot\left(\bm{k}\times{\bm{\nabla}_{z}f}\right),$
(4.65)

where $k=\left\lvert\bm{k}\right\rvert$ is the horizontal wavenumber. Positive
angles are measured counter-clockwise relative to $\bm{k}$. Thus,
$\Delta\theta_{f}>0$ indicates that $\bm{k}$ points to the right of
$\bm{\nabla}_{z}f$ while $\Delta\theta_{f}<0$ indicates that $\bm{k}$ points
to the left of $\bm{\nabla}_{z}f$.

We define the topographic parameter $\alpha$ by

$\alpha=\left\lvert f_{0}\,\bm{\nabla}_{z}h\right\rvert.$ (4.66)

In analogy with $\Delta\theta_{f}$, we define the angle $\Delta\theta_{h}$ by

$\sin{(\Delta\theta_{h})}=\frac{1}{k\,\alpha}\,\hat{\bm{z}}\cdot\left(\bm{k}\times{f_{0}\bm{\nabla}_{z}h}\right)$
(4.67)

with a similar interpretation assigned to $\Delta\theta_{h}>0$ and
$\Delta\theta_{h}<0$.

Substituting the wave solution (4.64) into the linear quasigeostrophic
equations (4.59)–(4.61) and assuming that $\alpha\,\sin(\Delta\theta_{h})\neq
0$, $\omega\neq 0$, and $k\neq 0$, we obtain

$\displaystyle-(S^{-1}\,\hat{\psi}^{\prime})^{\prime}+k^{2}\,\hat{\psi}=\lambda\,\hat{\psi}\quad$
$\displaystyle\text{for }z\in(-H,0)$ (4.68)
$\displaystyle-\frac{\beta}{\alpha}\,\frac{\sin{(\Delta\theta_{f})}}{\sin{(\Delta\theta_{h})}}\,S^{-1}\,\hat{\psi}^{\prime}=\lambda\,\psi\quad$
$\displaystyle\text{for }z=-H$ (4.69) $\displaystyle
S^{-1}\,\psi^{\prime}=0\quad$ $\displaystyle\text{for }z=0,$ (4.70)

where we have defined the eigenvalue $\lambda$ by

$\lambda=-\frac{k\,\beta\,\sin{(\Delta\theta_{f})}}{\omega}.$ (4.71)

Since $k\neq 0$ then $\lambda=0$ is not an eigenvalue. The above problem
(4.68)–(4.70) was recently considered in LaCasce (2017).

Figure 4.5: The streamfunction eigenfunctions $\hat{\psi}_{n}$ of the
quasigeostrophic eigenvalue problem with a sloping bottom from section 4.4.2.
Two cases are shown. The first is with $\Delta\theta_{f}=-90^{\circ}$ and
$\Delta\theta_{1}=-30^{\circ}$ and is both right-definite and left-definite.
The second is with $\Delta\theta_{f}=-45^{\circ}$ and
$\Delta\theta_{1}=15^{\circ}$ and is only left-definite. In the right-definite
case the $n$th eigenfunction has $n$ internal zero whereas in the left-
definite only case there are two eigenfunctions ($n=0,1$) with no internal
zeros.

##### Definiteness & the underlying function space

The eigenvalue problem has one $\lambda$-dependent boundary condition and so
the underlying function space is

$L^{2}_{\mu}\cong L^{2}\oplus\mathbb{C}.$ (4.72)

The appropriate inner product is obtained from equations (4.11) and (4.4)

$\displaystyle\left<\varphi,\phi\right>=\frac{1}{H}\left[\int_{-H}^{0}\varphi\,\phi\,\mathrm{d}z+\frac{\alpha}{\beta}\,\frac{\sin\left(\Delta\theta_{h}\right)}{\sin\left(\Delta\theta_{f}\right)}\,\varphi(-H)\,\phi(-H)\right]$
(4.73)

where we have introduced the factor $1/H$ for dimensional consistency in
eigenfunction expansions. By proposition 4.2.1, the problem is right-definite
for horizontal wavevectors $\bm{k}$ satisfying

$\frac{\sin{(\Delta\theta_{h})}}{\sin{(\Delta\theta_{f})}}>0$ (4.74)

and, in such cases, $L^{2}_{\mu}$ equipped with the inner product (4.11) is a
Hilbert space. However, $L^{2}_{\mu}$ is not a Hilbert space for all
wavevectors $\bm{k}$. By proposition 4.2.2, the problem is left-definite for
all wavevectors $\bm{k}$ and so $L^{2}_{\mu}$, equipped with the inner product
(4.11), is generally a Pontryagin space.

We write $\hat{\Psi}_{n}$ for the eigenfunctions and $\hat{\psi}_{n}$ for the
solutions of equations (4.68)–(4.70). The eigenfunctions $\hat{\Psi}_{n}$ are
related to the solutions $\hat{\psi}_{n}$ by (4.9) with boundary values
$\hat{\Psi}_{n}(0)$ given by equation (4.22). However, since $c_{1}=1$ and
$d_{1}=0$ in equation (4.69) [compare with equations (4.1)–(4.3)] then
$\hat{\Psi}_{n}=\hat{\psi}_{n}$ on the closed interval $[-H,0]$. Thus, the
solutions $\psi_{n}$ are also the eigenfunctions.

With theorem 4.2.3, we deduce that all eigenvalues $\lambda_{n}$ are real and
the corresponding eigenfunctions $\\{\hat{\psi}_{n}\\}_{n=0}^{\infty}$ form an
orthonormal basis for $L^{2}_{\mu}$. Orthonormality is defined with respect to
the inner product given by equation (4.73) and takes the form

$\pm\delta_{mn}=\left<\hat{\psi}_{m},\hat{\psi}_{n}\right>$ (4.75)

where we have taken the eigenfunctions $\hat{\psi}_{m}$ and $\hat{\psi}_{n}$
to be non-dimensional.

##### Properties of the eigenfunctions

By lemma 4.2.6, the number of internal zeros of the eigenfunctions
$\\{\hat{\psi}_{n}\\}_{n=0}^{\infty}$ depends on the propagation direction and
hence [by equation (4.74)] on the definiteness of the problem (see figure
4.5):

  * _1._

if the problem is right-definite then the $n$th eigenfunction has $n$ internal
zeros,

  * _2._

if the problem is not right-definite then both $\psi_{0}$ and $\psi_{1}$ have
no internal zeros; the remaining eigenfunctions $\psi_{n}$, for $n>1$, have
$n-1$ internal zeros.

As the problem is left-definite for all wavevectors $\bm{k}$, we can use
proposition 4.2.4 to determine the sign of the eigenvalues. Proposition 4.2.4
informs us that

$\lambda_{n}\left<\hat{\psi}_{n},\hat{\psi}_{n}\right>>0.$ (4.76)

In the first case, when the problem is right-definite, all eigenvalues are
positive and all eigenfunctions $\hat{\psi}_{n}$ satisfy
$\left<\hat{\psi}_{n},\hat{\psi}_{n}\right>>0$. In the second case, when the
problem is only left-definite, then there is one negative eigenvalue
$\lambda_{0}$ and the corresponding eigenfunction $\hat{\psi}_{0}$ satisfies
$\left<\hat{\psi}_{0},\hat{\psi}_{0}\right><0$. The remaining eigenvalues are
positive and their corresponding eigenfunctions satisfy
$\left<\hat{\psi}_{n},\hat{\psi}_{n}\right>>0$. In fact, from equation (4.71),
we see that waves with $\left<\hat{\psi}_{n},\hat{\psi}_{n}\right>>0$ have
westward phase speeds $\omega_{n}/k<0$ while waves with
$\left<\hat{\psi}_{n},\hat{\psi}_{n}\right><0$ have eastward phase speeds
$\omega_{n}/k>0$.

##### Expansion properties

The eigenfunctions $\\{\hat{\psi}_{n}\\}_{n=0}^{\infty}$ are complete in
$L^{2}_{\mu}$ but overcomplete in $L^{2}$. Physically, there is now an
additional eigenfunction corresponding to a topographic Rossby wave ($n=0$ in
figure 4.5).

Given a twice continuously differentiable function $\phi(z)$ satisfying
$\phi^{\prime}(0)=0$, then from theorem 4.2.5, we have

$\phi(z)=\sum_{n=0}^{\infty}\frac{\left<\phi,\hat{\psi}_{n}\right>}{\left<\hat{\psi}_{n},\hat{\psi}_{n}\right>}\,\hat{\psi}_{n}(z)\quad\textrm{and}\quad\phi^{\prime}(z)=\sum_{n=0}^{\infty}\frac{\left<\phi,\hat{\psi}_{n}\right>}{\left<\hat{\psi}_{n},\hat{\psi}_{n}\right>}\,\hat{\psi}^{\prime}_{n}(z),$
(4.77)

with both series converging uniformly on $[-H,0]$ (note that $\phi$ is not
required to satisfy any particular boundary condition at $z=-H$). If the
vertical structure at time $t=0$ (and at some wavevector $\bm{k}$) is given by
$\phi$, then the subsequent time-evolution is given by

$\psi(\bm{x},z,t)=\sum_{n=0}^{\infty}\frac{\left<\phi,\hat{\psi}_{n}\right>}{\left<\hat{\psi}_{n},\hat{\psi}_{n}\right>}\,\hat{\psi}_{n}(z)\,\cos\left(\omega_{n}t\right)\,\mathrm{e}^{\mathrm{i}\bm{k}\cdot\bm{x}},$
(4.78)

where the angular frequency $\omega_{n}$ is given by equation (4.71).

###  4.5 A localized perturbation at the boundary

We now consider a localized perturbation at a dynamically-active boundary; we
idealize such a perturbation by a boundary step-function $\Theta_{i}$ (for
$i\in S$) given by

$\Theta_{i}(z)=\begin{cases}1\quad\text{if }z=z_{i}\\\
0\quad\text{otherwise.}\end{cases}$ (4.79)

Using equation (4.23), the series expansion of $\Theta_{i}$ is found to be

$\Theta_{i}=\frac{1}{D_{i}}\sum_{n=0}^{\infty}\frac{\,\Phi_{n}(z_{i})}{\left<\Phi_{n},\Phi_{n}\right>}\,\Phi_{n}(z).$
(4.80)

For the non-rotating Boussinesq problem of section 4.3.2, a step-function
perturbation with amplitude $w_{0}$ (at some wavevector $\bm{k}$) yields the
time-evolution

$w(\bm{x},z,t)=w_{0}\left(\frac{g_{b}+\tau\,k^{2}}{N_{0}^{2}\,H}\right)\sum_{n=0}^{\infty}\hat{w}_{n}(0)\,\hat{w}_{n}(z)\,\cos\left(\sigma_{n}kt\right)\,\mathrm{e}^{\mathrm{i}\bm{k}\cdot\bm{x}}.$
(4.81)

Analogously, for the quasigeostrophic problem of section 4.4.2, a step-
function perturbation with amplitude $\psi_{0}$ (at some wavevector $\bm{k}$)
yields the time-evolution

$\psi(\bm{x},z,t)=\psi_{0}\left[\frac{\alpha\,\sin(\Delta\theta_{h})}{H\,\beta\,\sin(\Delta\theta_{f})}\right]\sum_{n=0}^{\infty}\frac{\hat{\psi}_{n}(-H)}{\left<\hat{\psi}_{n},\hat{\psi}_{n}\right>}\,\hat{\psi}_{n}(z)\,\cos\left(\omega_{n}t\right)\,\mathrm{e}^{\mathrm{i}\bm{k}\cdot\bm{x}}.$
(4.82)

That both the above series converge to a step-function at $t=0$ (and
$\bm{x}=\bm{0}$) is confirmed by theorem 4.A.3 along with theorem 2 in Fulton
(1977).

We thus see that a step-function perturbation induces wave motion with an
amplitude that is proportional to the boundary-confined restoring force (at
wavevector $\bm{k}$). Moreover, the amplitude of each constituent wave in the
resulting motion is proportional to the projection of that wave onto the
dynamically-active boundary.

###  4.6 Summary and conclusions

We have developed a mathematical framework for the analysis of three-
dimensional wave problems with dynamically-active boundaries (i.e., boundaries
where time derivatives appear in the boundary conditions). The resulting waves
have vertical structures that depend on the wavevector $\bm{k}$: For
Boussinesq gravity waves, the dependence is only through the wavenumber $k$
whereas the dependence for quasigeostrophic Rossby waves is on both the
wavenumber $k$ and the propagation direction $\bm{k}/k$. Moreover, the
vertical structures of the waves are complete in a space larger than $L^{2}$,
namely, they are complete in $L^{2}_{\mu}\cong L^{2}\oplus\mathbb{C}^{s}$
where $s$ is the number of dynamically active boundaries (and the number of
boundary-trapped waves). Each dynamically active boundary contributes an
additional boundary-trapped wave and hence an additional degree of freedom to
the problem. Mathematically, the presence of boundary-trapped waves allows us
to expand a larger collection of functions (with a uniformly convergent
series) in terms of the modes. The resulting series are term-by-term
differentiable and the differentiated series converges uniformly. In fact, the
normal modes have the intriguing property converging pointwise to functions
with finite jump discontinuities at the boundaries, a property related to
their ability to expand distributions in the Bretherton (1966)
“$\delta$-function formulation” of a physical problem. By considering a step-
function perturbation at a dynamically-active boundary, we find that the
subsequent time-evolution consists of waves whose amplitude is proportional to
their projection at the dynamically-active boundary. Within the mathematical
formulation is a qualitative oscillation theory relating the number of
internal zeros of the eigenfunctions to physical quantities; indeed, for the
quasigeostrophic problem, the number of zeros of the topographic Rossby wave
depends on the propagation direction while, for the rotating Boussinesq
problem, the ratio of the Coriolis parameter to the horizontal wavenumber
determines at which integer $M$ we obtain two modes with $M$ zeros.

Our results also clarify the difference between the traditional
quasigeostrophic baroclinic modes and the the $L^{2}\oplus\mathbb{C}^{2}$
eigenfunctions of Smith and Vanneste (2012). Namely, the series expansion of a
function in terms of the Smith and Vanneste (2012) eigenfunctions has a term-
by-term derivative that converges uniformly over the whole interval regardless
of the boundary conditions satisfied by the function. In contrast, an
eigenfunction expansion in terms of the baroclinic modes only has this
property if the function satisfies the same boundary conditions as the
baroclinic modes. One consequence is the following. Suppose we expand an
arbitrary quasigeostrophic state, with boundary buoyancy anomalies, in terms
of the baroclinic modes. The presence of these boundary buoyancy anomalies
implies that this state does not satisfy the same boundary conditions as the
baroclinic modes. The resulting series expansion in term of the baroclinic
modes is then not differentiable at the boundaries. We are thus unable to
recover the value of the boundary buoyancy anomalies from the series expansion
and so we have lost information in the expansion process. This loss of
information does not occur with $L^{2}\oplus\mathbb{C}^{2}$ expansions.

Normal mode decompositions of quasigeostrophic motion play an important role
in physical oceanography (e.g., Wunsch, 1997; Lapeyre, 2009; LaCasce, 2017).
Other applications include the extension of equilibrium statistical mechanical
calculations (e.g., Bouchet and Venaille, 2012; Venaille et al., 2012) to
three-dimensional systems with dynamically-active boundaries. Moreover, the
mathematical framework developed here is useful for the development of weakly
non-linear wave turbulence theories (e.g., Fu and Flierl, 1980; Smith and
Vallis, 2001; Scott, 2014) in systems with both internal and boundary-trapped
waves.

##  Appendix 4.A Additional properties of the eigenvalue problem

####  4.A.1 Construction of $L^{2}_{\mu}$

First, define the weighted Lebesgue measure $\sigma$ by

$\sigma([a,b])=\int_{a}^{b}r\,\mathrm{d}z\quad\text{where
}a,b\in[z_{1},z_{2}].$ (4.83)

The measure $\sigma$ induces the differential element

$\mathrm{d}\sigma(z)=r(z)\,\mathrm{d}z$ (4.84)

and is the measure associated with $L^{2}$ [see equations (4.6) and (4.7)].

Now, for $i\in S$ [see equation (4.8)], define the pure point measure
$\nu_{i}$ by (e.g., Reed and Simon, 1980, section I.4, example 2)

$\nu_{i}([a,b])=\begin{cases}D_{i}^{-1}\ \ \ &\text{if }z_{i}\in[a,b]\\\ 0\ \
\ &\text{otherwise,}\end{cases}$ (4.85)

where $D_{i}$ is the combination of boundary condition coefficients given by
equation (4.4). The pure point measure $\nu_{i}$ induces the differential
element

$\mathrm{d}\nu_{i}(z)=D_{i}^{-1}\,\delta(z-z_{i})\,\mathrm{d}z,$ (4.86)

where $\delta(z)$ is the Dirac distribution.

Consider now the space $L^{2}_{\nu_{i}}$ of “functions” $\phi$ satisfying

$\left\lvert\int_{z_{1}}^{z_{2}}\left\lvert\phi\right\rvert^{2}\,\mathrm{d}\nu_{i}\right\rvert=\left\lvert
D_{i}^{-1}\right\rvert\,\int_{z_{1}}^{z_{2}}\left\lvert\phi\right\rvert^{2}\,\delta(z-z_{i})\,\mathrm{d}z=\left\lvert
D_{i}^{-1}\right\rvert\,\left\lvert\phi(z_{i})\right\rvert^{2}<\infty.$ (4.87)

Elements of $L^{2}_{\nu_{i}}$ are not functions, but rather equivalence
classes of functions. Two functions, $\phi$ and $\psi$, on the interval
$[z_{1},z_{2}]$ are equivalent in $L^{2}_{\nu_{i}}$ if
$\phi(z_{i})=\psi(z_{i})$. In particular, $L^{2}_{\nu_{i}}$ is a one-
dimensional vector space and is hence isomorphic to the field of complex
numbers $\mathbb{C}$

$L^{2}_{\nu_{i}}\cong\mathbb{C}.$ (4.88)

Now define the measure $\mu$ by

$\mu=\sigma+\sum_{i\in S}\nu_{i}$ (4.89)

with an induced differential element of

$\mathrm{d}\mu(z)=\left[r(z)+\sum_{i\in
S}D_{i}^{-1}\,\delta(z-z_{i})\right]\mathrm{d}z.$ (4.90)

Then $L^{2}_{\mu}$ is the space of equivalence classes of functions that are
square-integrable with respect to the measure $\mu$.

Since the measures $\sigma$ and $\nu_{i}$, for $i\in S$, are mutually
singular, we have (Reed and Simon, 1980, section II.1, example 5)

$L^{2}_{\mu}\cong L^{2}\oplus\sum_{i\in S}L^{2}_{\nu_{i}}\cong
L^{2}\oplus\mathbb{C}^{s}$ (4.91)

from which we see that $L^{2}_{\mu}$ is “larger” by $s$ dimensions.

####  4.A.2 The eigenvalue problem in $L^{2}_{\mu}$

We construct here an operator formulation of (4.1)–(4.3) as an eigenvalue
problem in the Pontryagin space $L^{2}_{\mu}$.

Define the differential operator $\ell$ acting on a function $\phi$ by

$\ell\,\phi=\frac{1}{r}\,\left[(p\,\phi^{\prime})^{\prime}-q\,\phi\right].$
(4.92)

We also define the following boundary operators for $i\in S$,

$\displaystyle\mathcal{B}_{i}\phi$
$\displaystyle=\left[a_{i}\,\phi(z_{i})-b_{i}\,(p\,\phi^{\prime})(z_{i})\right]$
(4.93) $\displaystyle\mathcal{C}_{i}\phi$
$\displaystyle=\left[c_{i}\,\phi(z_{i})-d_{i}\,(p\,\phi^{\prime})(z_{i})\right].$
(4.94)

Let $\Phi$ be an element of $L^{2}_{\mu}$, as in equation (4.9), with boundary
values $\Phi(z_{i})=\mathcal{C}_{i}\phi$ for $i\in S$ and equal to $\phi$
elsewhere. We then define the operator $\mathcal{L}$, acting on functions
$\Phi$, by

$\mathcal{L}\,\Phi=\begin{cases}-\ell\,\phi\quad&\text{for
}z\in(z_{1},z_{2})\\\ -\mathcal{B}_{i}\,\phi\quad&\text{for }z=z_{i}\text{
where }i\in S\end{cases}$ (4.95)

with a domain $D(\mathcal{L})\subset L^{2}_{\mu}$ defined by

$\displaystyle D(\mathcal{L})=\\{\Phi\in L^{2}_{\mu}\ |$ $\displaystyle\
\phi\text{ is continuously differentiable, }\ell\,\phi\in L^{2},\
\Phi(z_{i})=\mathcal{C}_{i}\,\phi$ (4.96) $\displaystyle\text{ for }i\in S$
$\displaystyle\text{ and }\mathcal{B}_{i}\phi=0\text{ for
}i\in\\{1,2\\}\setminus S\\}.$

Recall that $S$ contains indices of the $\lambda$-dependent boundary
conditions, and therefore, $\\{1,2\\}\setminus S$ contains the indices of the
$\lambda$-independent boundary conditions.

Then, on the subspace $D(\mathcal{L})$ of $L^{2}_{\mu}$, the eigenvalue
problem (4.1)–(4.3) may be written as

$\mathcal{L}\,\Phi=\lambda\,\Phi.$ (4.97)

As shown in Russakovskii (1975, 1997), $\mathcal{L}$ is a self-adjoint
operator in the space $L^{2}_{\mu}$.

There is a natural quadratic form $Q$, induced by the eigenvalue problem
(4.1)–(4.3), given by

$Q(\Phi,\Psi)=\left<\Phi,\mathcal{L}\,\Psi\right>.$ (4.98)

For elements $\Phi,\Psi\in D(\mathcal{L})$, we obtain

$\displaystyle Q(\Phi,\Psi)$
$\displaystyle=\int_{z_{1}}^{z_{2}}\left[p\,{\phi^{\prime}}^{\ast}\,\psi^{\prime}+q\,{\phi}^{\ast}\,\psi\right]\,\mathrm{d}z+\sum_{i\in\\{1,2\\}\setminus
S}(-1)^{i+1}\,\frac{a_{i}}{b_{i}}\,{\phi(z_{i})}^{\ast}\,\psi(z_{i})\,$ (4.99)
$\displaystyle\quad-\sum_{i\in
S}\frac{1}{D_{i}}{\left(\begin{matrix}\psi(z_{i})\\\
-(p\,\psi^{\prime})(z_{i})\end{matrix}\right)}^{\ast}\cdot\left(\begin{matrix}a_{i}\,c_{i}&a_{i}\,d_{i}\\\
a_{i}\,d_{i}&b_{i}\,d_{i}\end{matrix}\right)\left(\begin{matrix}\phi(z_{i})\\\
-(p\,\phi^{\prime})(z_{i})\end{matrix}\right)$

for $b_{i}\neq 0$ for $i\in\\{1,2\\}\setminus S$. If $b_{i}=0$ for
$i\in\\{1,2\\}\setminus S$ then we replace the term $a_{i}/b_{i}$ with zero.

To develop the reality and completeness theorem 4.2.3, we provide the
following definitions.

######  Definition 4.A.1 (Right-definite).

The eigenvalue problem (4.1)–(4.3) is said to be right-definite if
$L^{2}_{\mu}$ is a Hilbert space or, equivalently, if

$\left<\Phi,\Phi\right>>0$ (4.100)

for all non-zero $\Phi\in L^{2}_{\mu}$.

######  Definition 4.A.2 (Left-definite).

The eigenvalue problem (4.1)–(4.3) is said to be left-definite if

$Q(\Phi,\Phi)\geq 0$ (4.101)

for all $\Phi\in D(\mathcal{L})$.

One can then prove propositions 4.2.1 and 4.2.2 through straightforward
manipulations.

####  4.A.3 Properties of eigenfunction expansions

The following theorem features some of the novel properties of the basis
$\\{\Phi_{n}\\}_{n=0}^{\infty}$ of $L^{2}_{\mu}$. Theorem 4.A.3 below is a
generalization of a theorem first formulated, in the right-definite case, by
Walter (1973) and Fulton (1977).

######  Theorem 4.A.3 (Eigenfunction expansions).

Let $\\{\Phi_{n}\\}_{n=0}^{\infty}$ be the set of eigenfunctions of the
eigenvalue problem (4.1)–(4.3). Then the following properties hold.

  * (i)

Null series: For $i\in S$, we have

$0=D_{i}^{-1}\,\sum_{n=0}^{\infty}\frac{1}{\left<\Phi_{n},\Phi_{n}\right>}\,\Phi_{n}(z_{i})\,\phi_{n}(z)$
(4.102)

with equality in the sense of $L^{2}$.

  * (ii)

Unit series: For $i\in S$, we have

$1=D_{i}^{-1}\,\sum_{n=0}^{\infty}\frac{1}{\left<\Phi_{n},\Phi_{n}\right>}\,\left\lvert\Phi_{n}(z_{i})\right\rvert^{2}.$
(4.103)

  * (iii)

$L^{2}$-expansion: Let $\psi\in L^{2}$, then

$\psi=\sum_{n=0}^{\infty}\frac{1}{\left<\Phi_{n},\Phi_{n}\right>}\,\left(\int_{z_{1}}^{z_{2}}{\psi}^{\ast}\,\phi_{n}\,r\,\mathrm{d}z\right)\,\phi_{n}.$
(4.104)

with equality in the sense of $L^{2}$.

  * (iv)

Interior-boundary orthogonality: Let $\psi\in L^{2}$, then for $i\in S$, we
have

$0=\sum_{n=0}^{\infty}\frac{1}{\left<\Phi_{n},\Phi_{n}\right>}\left(\int_{z_{1}}^{z_{2}}{\psi}^{\ast}\,\phi_{n}\,r\,\mathrm{d}z\right)\Phi_{n}(z_{i}).$
(4.105)

###### Proof.

The proof is similar to the proof of corollary 1.1 in Fulton (1977). ∎

####  4.A.4 Pointwise convergence and Sturm-Liouville series

Theorem 3 in Fulton (1977) states that the $\Phi_{n}$ series expansion (4.23)
behaves like a Fourier series in the interior of the interval $(z_{1},z_{2})$
(see appendix 4.B for why this theorem applies in the left-definite case).
Since the expansions (4.23) and (4.24) in terms of $\Phi_{n}$ and $\phi_{n}$
are equal in the interior, then the above theorem applies to the $\phi_{n}$
series (4.24) as well. It is at the boundaries points, $z=z_{1},z_{2}$, where
the novel behaviour of the series expansions (4.23) and (4.24) appears.

For traditional Sturm-Liouville expansions [with eigenfunctions of problem
(4.1)-(4.3) with $c_{i},d_{i}=0$ for $i=1,2$], eigenfunction expansions behave
like the analogous Fourier series on $[z_{1},z_{2}]$ [page 16 in Titchmarsh
(1962) or chapter 1, section 9, in Levitan and Sargsjan (1975)]. In
particular, for a twice continuously differentiable function $\psi$, the
eigenfunction expansion of $\psi$ converges uniformly to $\psi$ on
$[z_{1},z_{2}]$ so long as the eigenfunctions $\phi_{n}$ do not vanish at the
boundaries. If the eigenfunctions vanish at one of the boundaries, then we
only obtain uniform convergence if $\psi$ vanishes at the corresponding
boundary as well (Brown and Churchill, 1993, section 22). Under these
conditions, the resulting expansion will be differentiable in the interior of
the interval, $(z_{1},z_{2})$, but not at the boundaries $z=z_{1},z_{2}$ [see
chapter 8, section 3, in Levitan and Sargsjan (1975) for the equiconvergence
of differentiated Sturm-Liouville series with Fourier series and see section
23 in Brown and Churchill (1993) for the convergence behaviour of
differentiated Fourier series].

Returning to the case of eigenfunction expansions for the eigenvalue problem
(4.1)–(4.3) with $\lambda$-dependent boundaries, the following theorem
provides pointwise (as well as uniform, in the case $d_{i}\neq 0$) convergence
conditions for the $\phi_{n}$ series (4.24).

######  Theorem 4.A.4 (Pointwise convergence).

Let $\psi$ be a twice continuously differentiable function on the interval
$[z_{1},z_{2}]$ satisfying any $\lambda$-independent boundary conditions in
the eigenvalue problem (4.1)–(4.3). Define the function $\Psi$ on
$[z_{1},z_{2}]$ by

$\Psi(z)=\begin{cases}\Psi(z_{i})\quad&\textrm{at }z=z_{i},\textrm{ for }i\in
S,\\\ \psi(z)\quad&\textrm{otherwise}.\end{cases}$ (4.106)

where $\Psi(z_{i})$ are constants for $i\in S$ (the $\lambda$-dependent
boundaries). Then we have the following.

  * (i)

If $d_{i}\neq 0$ for $i\in S$, then the $\phi_{n}$ series expansion (4.24)
converges uniformly to $\psi(z)$ on the closed interval $[z_{1},z_{2}]$,

$\sum_{n=0}^{\infty}\frac{\left<\Psi,\Phi_{n}\right>}{\left<\Phi_{n},\Phi_{n}\right>}\,\phi_{n}(z)=\psi(z).$
(4.107)

Furthermore, for the differentiated series, we have

$\sum_{n=0}^{\infty}\frac{\left<\Psi,\Phi_{n}\right>}{\left<\Phi_{n},\Phi_{n}\right>}\,\phi_{n}^{\prime}(z)=\begin{cases}\left(c_{i}\,\psi(z_{i})-\Psi(z_{i})\right)/d_{i}\quad&\textrm{at
}z=z_{i},\textrm{ for }i\in S\\\
\psi^{\prime}(z)\quad&\textrm{otherwise}.\end{cases}$ (4.108)

  * (ii)

If $d_{i}=0$, then we have

$\sum_{n=0}^{\infty}\frac{\left<\Psi,\Phi_{n}\right>}{\left<\Phi_{n},\Phi_{n}\right>}\,\phi_{n}=\begin{cases}\Psi(z_{i})/c_{i}\quad&\textrm{at
}z=z_{i},\textrm{ for }i\in S\\\ \psi(z)\quad&\textrm{otherwise}.\end{cases}$
(4.109)

###### Proof.

This theorem is a generalization of corollary 2.1 in Fulton (1977). We provide
the extension of the corollary to the left-definite problem in appendix 4.B.4.
∎

The $\Phi_{n}$ series (4.23) converges to $\Psi(z_{i})$ at $z=z_{i}$ for $i\in
S$ (i.e., at $\lambda$-dependent boundaries) but otherwise behaves as in
theorem 4.A.4.

##  Appendix 4.B Literature survey and mathematical proofs

####  4.B.1 Literature survey

There is an extensive literature associated with the eigenvalue problem
(4.1)–(4.3) with $\lambda$-dependent boundary conditions (see Schäfke and
Schneider, 1966; Fulton, 1977, and citations within). One can use the
$S$-hermitian theory of Schäfke and Schneider (1965, 1966, 1968) to show that
one obtains real eigenvalues when the problem is either right-definite or
left-definite (see section 4.2) but completeness results in $L^{2}_{\mu}$ are
unavailable in this theory.

The right-definite theory is well-known (Evans, 1970; Walter, 1973; Fulton,
1977). In particular, Fulton (1977) applies the residue calculus techniques of
Titchmarsh (1962) to the right-definite problem and, in the process, extends
some well-known properties of Fourier series to eigenfunction expansions
associated with (4.1)–(4.3). A recent Hilbert space approach to the right-
definite problem, in the context of obtaining a projection basis for
quasigeostrophic dynamics, is given by Smith and Vanneste (2012).

The left-definite problem is less examined. As we show in this chapter, the
eigenvalue problem is naturally formulated in a Pontryagin space, and, in such
a setting, one can prove, in the left-definite case, that the eigenvalues are
real and that the eigenfunctions form a basis for the underlying function
space. We prove this result, stated in theorem 4.2.3, in appendix 4.B.3.

With these completeness results, we may apply the residue calculus techniques
of Titchmarsh (1962) to extend the results of Fulton (1977) to the left-
definite problem. Indeed, Fulton (1977) uses a combination of Hilbert space
methods as well as residue calculus techniques to prove various convergence
results for the right-definite problem. However, only theorem 1 of Fulton
(1977) makes use of Hilbert space methods. If we extend Fulton’s theorem 1 to
the left-definite problem, then all the results of Fulton (1977) will apply
equally to the left-definite problem. A left-definite analogue of theorem 1 of
Fulton (1977), along with its proof, is given in appendix 4.B.4.

####  4.B.2 A Pontryagin space theorem

A Pontryagin space $\Pi_{\kappa}$, for a finite non-negative integer $\kappa$,
is a Hilbert space with a $\kappa$-dimensional subspace of elements satisfying

$\displaystyle\left<\phi,\phi\right><0.$ (4.110)

An introduction to the theory of Pontryagin spaces can be found in Iohvidov
and Krein (1960) as well as in the monograph of Bognár (1974). Another
resource is the monograph of Azizov and Iokhvidov (1989) on linear operators
in indefinite inner product spaces.

Pontryagin spaces admit a decomposition

$\Pi_{\kappa}=\Pi^{+}\oplus\Pi^{-}$ (4.111)

into orthogonal subspaces $(\Pi^{+},+\left<\cdot,\cdot\right>)$ and
$(\Pi^{-},-\left<\cdot,\cdot\right>)$. Moreover, one can associate with a
Pontryagin space $(\Pi_{\kappa},\left<\cdot,\cdot\right>)$ a corresponding
Hilbert space $(\Pi,\left<\cdot,\cdot{}\right>_{+})$ where the positive-
definite inner product $\left<\cdot,\cdot\right>_{+}$ is defined by

$\left<\phi,\psi\right>_{+}=\left<\phi_{+},\psi_{+}\right>-\left<\phi_{-},\psi_{-}\right>,\quad\phi,\psi\in\Pi,$
(4.112)

where $\phi=\phi_{+}+\phi_{-}$ and $\psi=\psi_{+}+\psi_{-}$, with
$\phi_{\pm},\psi_{\pm}\in\Pi^{\pm}$ (Azizov and Iokhvidov, 1981).

As a prerequisite to proving theorem 4.2.3, we require the following.

######  Theorem 4.B.1 (Positive compact Pontryagin space operators).

Let $\mathcal{A}$ be a positive compact operator in a Pontryagin space
$\Pi_{\kappa}$ and suppose that $\lambda=0$ is not an eigenvalue. Then all
eigenvalues are real and the corresponding eigenvectors form an orthonormal
basis for $\Pi_{\kappa}$. There are precisely $\kappa$ negative eigenvalues
and the remaining eigenvalues are positive. Moreover, positive eigenvalues
have positive eigenvectors and negative eigenvalues have negative
eigenvectors.

###### Proof.

By theorem VII.1.3 in Bognár (1974) the eigenvalues are all real. Moreover,
since $\lambda=0$ is not an eigenvalue, then all eigenspaces are definite