Topological Origin of Equatorial Waves

TL;DR

This paper reveals that equatorial waves like Kelvin and Yanai modes have a topological origin, with their existence guaranteed by the Earth's rotation and the associated Chern number, linking climate dynamics to topological insulators.

Contribution

It demonstrates a topological explanation for equatorially trapped waves, connecting geophysical fluid dynamics with topological phases of matter.

Findings

Kelvin and Yanai waves are topologically protected due to Earth's rotation.

The bulk Poincaré wave modes have a non-trivial Chern number of 2.

Topology underpins the robustness of these oceanic and atmospheric waves.

Abstract

Topology sheds new light on the emergence of unidirectional edge waves in a variety of physical systems, from condensed matter to artificial lattices. Waves observed in geophysical flows are also robust to perturbations, which suggests a role for topology. We show a topological origin for two celebrated equatorially trapped waves known as Kelvin and Yanai modes, due to the Earth's rotation that breaks time-reversal symmetry. The non-trivial structure of the bulk Poincar\'e wave modes encoded through the first Chern number of value $2$ guarantees existence for these waves. This invariant demonstrates that ocean and atmospheric waves share fundamental properties with topological insulators, and that topology plays an unexpected role in the Earth climate system.

Equations

$${\partial }_{t}h+\nabla \cdot \left(h\mathbf{u}\right)$$

$${\partial }_t\mathbf{u}+(\mathbf{u}\cdot \nabla )\mathbf{u}$$

$$L_d=\frac{c}{2\Omega },c=\sqrt{gH}.$$

$$\tilde{\eta }=\frac{h-H}{H},\tilde{\mathbf{u}}=\frac{\mathbf{u}}{c},\tilde{t}=2\Omega t,\tilde{f}=\frac{f}{2\Omega },\tilde{\mathbf{x}}=\frac{\mathbf{x}}{L_d},$$

$${\partial }_{\tilde{t}}\tilde{\eta }$$

$${\partial }_{\tilde{t}}\tilde{u}$$

$${\partial }_{\tilde{t}}\tilde{v}$$

$$\omega \left(\begin{array}{c}\hat{\eta }\\ \hat{u}\\ \hat{v}\end{array}\right)=\left(\begin{array}{ccc}0& k_x& k_y\\ k_x& 0& -if\\ k_y& if& 0\end{array}\right)\left(\begin{array}{c}\hat{\eta }\\ \hat{u}\\ \hat{v}\end{array}\right).$$

$${\Psi }_{+}(k_x,k_y,f)=\frac{1}{\sqrt{2}}\left(\begin{array}{c}\frac{k}{\sqrt{k^2+f^2}}\\ \frac{k_x}{k}-i\frac{f{k}_y}{k\sqrt{k^2+f^2}}\\ \frac{k_y}{k}+i\frac{f{k}_x}{k\sqrt{k^2+f^2}}\end{array}\right)$$

$${\Psi }_0(k_x,k_y,f)=\frac{1}{\sqrt{k^2+f^2}}\left(\begin{array}{c}f\\ i{k}_y\\ -i{k}_x\end{array}\right)$$

$${\Psi }_{+}=\frac{1}{\sqrt{2}}\left(\begin{array}{c}\sin \theta \\ \cos \phi -i\cos \theta \sin \phi \\ \sin \phi +i\cos \theta \cos \phi \end{array}\right).$$

$${\Psi }_{+}^N={\Psi }_{+}{e}^{i\phi },{\Psi }_{+}^S={\Psi }_{+}{e}^{-i\phi }$$

$${\mathbf{A}}_{+}^N=-i⟨{\Psi }_{+}^N,{\nabla }_s{\Psi }_{+}^N⟩,{\mathbf{A}}_{+}^S=-i⟨{\Psi }_{+}^S,{\nabla }_s{\Psi }_{+}^S⟩$$

$${\mathbf{A}}_{+}^S={\mathbf{A}}_{+}^N-2{\nabla }_s\phi .$$

$${\mathbf{B}}_{+}={\nabla }_s\times {\mathbf{A}}_{+}^N={\nabla }_s\times {\mathbf{A}}_{+}^S.$$

$$\Delta {\mathcal{C}}_{+}=\frac{1}{2\pi }{\int }_0^{\pi }\mathrm{d}\theta {\int }_0^{2\pi }\mathrm{d}\phi {\mathbf{B}}_{+}\cdot {\mathbf{e}}_k$$

$$\Delta {\mathcal{C}}_{+}=\frac{1}{2\pi }{\int }_0^{2\pi }\mathrm{d}\phi \left({\mathbf{A}}_{+}^N-{\mathbf{A}}_{+}^S\right)\cdot {\mathbf{e}}_{\phi }.$$

$$\omega u_{m,n}$$

$$\omega v_{m,n}$$

$$\omega h_{m,n}$$

$$f_n=\tanh \left((n-1)\frac{\delta y}{R_f}\right)-\tanh \left((n-N_y)\frac{\delta y}{R_f}\right)+\tanh \left((n-2N_y)\frac{\delta y}{R_f}\right)$$

$$\omega u_{k_x,n}$$

$$\omega v_{k_x,n}$$

$$\omega h_{k_x,n}$$

$$\frac{\partial }{\partial t}\delta =\dots +A{g}_{10^{\circ }}(\theta -\pi /2)g_{\pi /m}(\varphi -\omega t/m-\pi /m)\sin (\omega t-m\varphi )$$

Full text

Topological Origin of Equatorial Waves

Pierre Delplace1∗ J. B. Marston2∗ Antoine Venaille1∗

1Univ Lyon, Ens de Lyon, Univ Claude Bernard, CNRS, Laboratoire de Physique

AF-69342 Lyon, France

2 Department of Physics, Box 1843, Brown University, Providence, RI 02912-1843 USA.

∗ E-mails: pierre.delplace@ens-lyon.fr, marston@brown.edu, antoine.venaille@ens-lyon.fr

Topology sheds new light on the emergence of unidirectional edge waves in a variety of physical systems, from condensed matter to artificial lattices. Waves observed in geophysical flows are also robust to perturbations, which suggests a role for topology. We show a topological origin for two celebrated equatorially trapped waves known as Kelvin and Yanai modes, due to the Earth’s rotation that breaks time-reversal symmetry. The non-trivial structure of the bulk Poincaré wave modes encoded through the first Chern number of value $2$ guarantees existence for these waves. This invariant demonstrates that ocean and atmospheric waves share fundamental properties with topological insulators, and that topology plays an unexpected role in the Earth climate system.

One sentence summary: Equatorial Kelvin and Yanai waves, important components of Earth’s climate system, have an origin in topology.

Symmetries and topology are central to an understanding of physics. In condensed matter, topology explains the precise quantization of the Hall effect  (?), where a magnetic field breaks the discrete symmetry of time reversal. Interest in topological properties was reinvigorated following the discovery of the quantum spin Hall effect, and the subsequent classification of different states of matter according to discrete symmetries  (?). Recently topologically protected edge excitations have been found in artificial lattices of various types  (?, ?, ?). A correspondence between topological properties of waves in the bulk and the existence of unidirectional edge modes along boundaries exists in all these systems  (?, ?). The edge modes fill frequency or energy gaps found in the bulk and are immune to various types of disorder. We show that topologically protected edge waves also manifest in atmospheres and oceans.

Equatorial Kelvin and mixed Rossby-gravity (Yanai) waves are edge modes that propagate energy along the Earth’s equator with eastward group velocity  (?). Remarkably, the dispersion relations for these waves (Figure 1a) were derived within the framework of the rotating shallow water model  (?) just prior to their first observation in the 1960s. Since then, observations of the atmosphere have revealed a robust signature of these trapped modes in wavenumber-frequency spectra  (?) (Figure 1b). Equatorial Kelvin and Yanai waves have been shown to play a crucial role in several aspects of climate dynamics. For instance, Kelvin waves are a key component of the El Niño Southern Oscillation, traveling across the waters of the Pacific ocean  (?). The waves are also part of the quasi-biennal oscillation in the stratosphere, and are thought to be an important component of the Madden Julian Oscillation in the troposphere  (?).

The fact that Yanai and Kelvin waves are equatorially trapped unidirectional modes filling a frequency gap between low frequency planetary (Rossby) wave band and the high frequency inertia-gravity (Poincaré) wave band  (?), as shown in Figure 1a, suggests they can be interpreted as topological boundary states, similar to those emerging in various topological insulating media. More precisely, bulk (Poincaré and/or Rossby) waves possess a topological property, that should be directly related to the existence of these two unidirectional boundary waves, by virtue of the bulk-boundary correspondence  (?, ?). According to this correspondence, the number of states inherited by a band when the zonal (directed along the equator) wavenumber $k_{x}$ varies from $-\infty$ to $+\infty$ is given by an integer-valued topological number called the first Chern number. The first Chern number is an integer that quantifies the number of phase singularities in a bundle of eigenmodes parameterized on a closed manifold. These singularities are somewhat analogous to amphidromic points ($\pm 2\pi$ phase vortices of tidal modes), but they occur in parameter space rather than in physical space. We demonstrate the existence of a non-trivial global structure in the bulk Poincaré modes as being encoded through the first Chern number of value $\pm 2$, thus ensuring the existence of $2$ unidirectional edge modes at the equator that fill the two frequency gaps, in agreement with the existence of Kelvin and Yanai waves. The existence of the frequency gap originates from a broken time-reversal symmetry of the flow model due to Earth’s rotation. The structure of tidal modes   (?) and bifurcations in large scale geophysical flow  (?) have previously invoked the effect of breaking time-reversal symmetry. Our study shows that another far reaching consequence of this broken symmetry is to confer non-trivial topological properties to bundles of fluid waves, giving rise to robust edge states.

The rotating shallow water equations  (?) that describe the dynamics of a thin layer of fluid on a two-dimensional surface of height $h(\mathbf{x},t)$ and horizontal velocity $\mathbf{u}(\mathbf{x},t)$ furnish a minimal model for equatorial waves:

[TABLE]

The Coriolis parameter $f=2\boldsymbol{\Omega}\cdot\hat{\mathbf{n}}$ is twice the projection of the planetary angular rotation vector on the local vertical axis $\hat{\mathbf{n}}$ and $g$ is the constant of gravitational acceleration. When linearized about a state of rest ($\mathbf{u}=0$) and mean height ($h=H$), this dynamical system may be rewritten as $i\partial_{t}\Psi=\mathcal{H}\Psi$, where $\Psi=(\mathbf{u},\eta)$ is a triplet of fields describing the two components of the perturbed velocity field and the perturbed height field $\eta$, and where $\mathcal{H}$ is a Hermitian operator  (?). Because the fields $(\mathbf{u},\eta)$ are real, the operator $\mathcal{H}$ is equal to the negative of its complex conjugate: $\Xi\mathcal{H}\Xi^{-1}=-\mathcal{H}$ where $\Xi$ is the operator that effects complex conjugation, with $\Xi^{2}=1$. In the quantum context, the operation is referred to as a particle-hole transformation because it inverts the spectrum. Time reversal symmetry $t\rightarrow-t$, $\mathbf{x}^{\prime}\rightarrow\mathbf{x}$, $\eta\rightarrow\eta$, $\mathbf{u}\rightarrow-\mathbf{u}$ is broken by non-zero Coriolis parameter $f\neq 0$ in Eq. (2). The broken symmetry generates gaps in the shallow water spectrum  (?).

The $f$-plane approximation commonly used in geophysics  (?) amounts to the neglect of Earth sphericity by assuming that the dynamics take place on a tangent plane with constant $f$ (Figure 2a). Translational symmetry ensures that eigenmodes of the linearized dynamics in this geometry are of the form $\hat{\Psi}\exp\left(i\omega t-ik_{x}x-ik_{y}y\right)$ where $\hat{\Psi}$ has three components. Viewing $f/c$ as an external parameter, where $c=\sqrt{gH}$ is the speed of non-rotating shallow water gravity waves, the eigenmodes may be easily found at each point in the space $(k_{x},k_{y},f/c)$ as depicted in Figure 2b. There are 3 bands with frequencies $\omega_{\pm}=\pm\left(f^{2}+c^{2}k^{2}\right)^{1/2}$ and $\omega_{0}=0$ where $k^{2}\equiv k_{x}^{2}+k_{y}^{2}$, with corresponding wavefunctions $\{\Psi_{\pm},~{}\Psi_{0}\}$. For $f\neq 0$, the bands separate by gaps of frequency $f$ (Figure 3). The zero-frequency modes are in geostrophic balance; the other two modes are Poincaré waves with dispersions $\omega_{\pm}$ that are symmetric with respect to the origin in ($k_{x},k_{y},\omega$) space.

Eigenmodes depend on the triplet of parameters $(k_{x},k_{y},f/c)$ that correspond to the set of waves in all possible $f$-plane configurations. The eigenmodes do not vary with the distance from the origin in $(k_{x},k_{y},f/c)$-space and can therefore be parameterized on the surface of a sphere $\mathcal{S}$ that encloses the singular band-crossing point at the origin $(k_{x},k_{y},f/c)=(0,0,0)$ [see Figure 2b and  (?)]. Each of the eigenstates $\left\{\Psi_{-},\Psi_{0},\Psi_{+}\right\}$ defines a fiber bundle over $\mathcal{S}$ that may possess topological defects. The singularities reflect the impossibility of continuously defining the eigenmodes everywhere on the sphere, and in particular over both of Earth’s two hemispheres simultaneously. They are quantified by the first Chern number $\Delta\mathcal{C}$ that can be calculated for each bulk mode $n$ as the flux of the Berry curvature $B_{n}=-{\rm i}\nabla_{p}\times\left(\Psi_{n}^{\dagger}\nabla_{p}\Psi_{n}\right)$ through the sphere $\mathcal{S}$, with $\Psi^{\dagger}_{n}$ the conjugate transpose of $\Psi_{n}$, and $\nabla_{p}=\left(\partial_{k_{x}},\partial_{k_{y}},\partial_{f/c}\right)$. In other words, there exists a quantized Berry flux generated by a (Berry) monopole located at the center of $\mathcal{S}$, where the three bands cross  (?, ?). The singularities are analogous to the one exhibited by an electron wavefunction that cannot be defined continuously around a Dirac magnetic monopole  (?). We find $\{\Delta\mathcal{C}_{-},\Delta\mathcal{C}_{0},\Delta\mathcal{C}_{+}\}=\{-2,0,2\}$  (?), namely only the Poincaré modes $\Psi_{\pm}$ are topologically non-trivial as the geostrophic modes $\Psi_{0}$ have zero Chern index, in agreement with the bulk-boundary correspondence  (?, ?).

To understand qualitatively the correspondence between these bulk properties and the emergence of unidirectional edge states in the presence of an equator, is it worth considering the case of a planar flow in an unbounded domain with $f$ varying in the $y$-direction from $-2\Omega$ to $2\Omega$ (see Figure 3). Far from the interface, the eigenmodes are given by delocalized solutions, i.e. by those computed in the case of constant $f$. If one could continuously deform the whole set of positive frequency eigenmodes from one hemisphere to other, for instance by varying $f$ slowly with y, then the eigenmodes would be given by solutions close to those calculated for constant $f$. Our previous calculation shows that this continuous deformation is prohibited by the occurance of $\Delta C_{+}=2$ phase singularities (positive vortices) when the plane $f=0$ is crossed. In order to remove these two singularities, the positive frequency band and the negative frequency bands must be connected to each other as the sum of their Chern numbers is zero. This connection happens through the emergence of two edge states that fill the frequency gaps. For any frequency that lies within the bulk gaps, the number of topological edge states is fixed by the set of Chern numbers  (?). Because $\Delta\mathcal{C}_{\pm}=\pm 2$ there are two extra unidirectional edge modes in the frequency gaps  (?).

It is instructive to examine the Berry curvature for the Poincaré modes. As shown in Figure 3, the curvature is mainly concentrated around $k=0$ where it reaches extremal values, and importantly, changes sign with $f$. It follows that its flux for each Poincaré mode $\mathcal{C}_{\pm}=\frac{1}{2\pi}\int_{-\infty}^{\infty}\rm{d}k_{x}\rm{d}k_{y}B_{\pm}=\pm\text{sgn}(f)$ is an integer that only depends on the hemisphere. It is thus tempting to say that the Poincaré eigenmodes on the two hemispheres are topologically distinct by interpreting $\mathcal{C}_{\pm}$ also as a Chern number, as the difference $\mathcal{C}_{\pm}(f>0)-\mathcal{C}_{\pm}(f<0)=\pm 2$ coincides with the first Chern number $\Delta\mathcal{C}_{\pm}$. This would be rigorously true if the two-dimensional manifold through which this Berry flux is computed at fixed $f$ were closed, for instance when the wave-numbers $(k_{x},k_{y})$ live on a Brillouin zone that reflects an underlying lattice. For continuous fluids, only $\Delta\mathcal{C}_{\pm}$ is a well defined topological number, but this suffices to characterize the topological property of the bulk modes, and thus the existence of the two equatorial unidirectional modes.

We stress one important point concerning the role of the spherical geometry of the planet in our approach. We removed this sphericity with the $f$-plane approximation, equivalent to holding the Coriolis parameter constant in space. However, through the construction of the sphere $\mathcal{S}$ in parameter space $(k_{x},k_{y},f/c)$, we recover the effect of a varying Coriolis parameter $f$ on the shallow water eigenmodes. In this way, sphericity works its way back into the problem. The detailed geometry of the Earth is no longer needed as topology itself requires the existence of Yanai and Kelvin waves. Even a misshapen sphere would support the waves.

Topology guarantees the existence of equatorial Yanai and Kelvin waves, obviating the need to carry out the classic but more complex calculation on the equatorial beta plane  (?). On the equatorial beta plane, Rossby and Poincaré waves can also be equatorially trapped. However, this trapping depends on the precise longitudinal variation of $f(y)$, as may be demonstrated numerically. In contrast, the topological origin of Kelvin and Yanai modes makes them insensitive to the details of the interface, such as the detailed shape of $f(y)$  (?). We also performed numerical scattering experiment showing that there is no possibility for Kelvin or Yanai wave excited within the bulk frequency gap, away from the other bands, to exchange energy with other modes that propagate energy westward  (?). Consequently, there is no energy backscattering in the presence of topography. This robustness against disorder can now be understood as a consequence of topology.

Other ideas from topology have been applied to hydrodynamics  (?, ?, ?). However, the appearance of singularities in the set of eigenmodes that arises from the breaking time-reversal symmetry has so far been overlooked in this context, as well as the striking physical consequence of unidirectional edge modes filling the frequency gaps. The general principle of bulk-boundary correspondence may now be applied to other fluid systems of interest.

The shallow water system exhibits particle-hole symmetry stemming real-valued velocity and displacement fields. More generally, any linearized fluid flow model that can be written in terms of an hermitian operator that breaks time reversal symmetry belong to the symmetry class with Cartan label D, which means that non trivial topological properties may arise  (?, ?). Other physical systems that may belong to class D are chiral p-wave superconductors  (?, ?) and superfluid 3He-A  (?). The linear operator of flow dynamics can be non-Hermitian in the presence of mean-flows and dissipation, in which case other topological properties may appear  (?). We expect that topology may ultimately shape the global structure of a number of other astrophysical and geophysical wave spectra, where similar gaps opened in the presence of symmetry breaking fields are known to exist. For instance, Lamb waves are edge states that fill the gap between acoustic and gravity waves because gravity breaks another discrete symmetry, that of inversion. Hall magnetohydrodynamics is another possible setting for topological edge waves  (?). It will also be also interesting to study in more detail the resilience of topological waves against dissipation, and non-linear wave-wave scattering processes.

Acknowledgments

We thank David Carpentier, Baylor Fox-Kemper, Thibaud Louvet, Jim Sauls, and Steve Tobias for discussions. P.D. was supported by the French Agence Nationale de la Recherche (ANR) under Grant TopoDyn (ANR-14-ACHN-0031).

Supplementary materials

1- First Chern Number of Shallow Water Eigenmodes

2- Numerical Calculation of Eigenmodes With $f(y)$: From the Shallow Water Model to a Haldane-like Model on the Lieb lattice

3 -Scattering of Equatorial Kelvin Waves by a Static Perturbation

Figs. S1 to S4

References (28-33)

Supplementary Material for “Topological Origin of Equatorial Waves”

1 First Chern Number of Shallow Water Eigenmodes

In the classical dispersion of equatorial waves on a beta plane ($f=\beta y$, Figure 1a in the Main Text), one recognizes the three bands obtained within the $f$-plane approximation, except that part of the degeneracy of geostrophic modes is lifted, with the emergence of Planetary (Rossby) waves due to the spatial variation of the Coriolis parameter. In addition, when varying the zonal (directed along the equator) wavenumber $k_{x}$ from $-\infty$ to $+\infty$, the Poincaré band of positive (negative) frequency acquires (is depleted of) two states whereas the total number of states acquired by the Rossby band vanishes. Thus, the occurence of one-way equatorial waves in the dispersion relation can be seen as a spectral flow from one band to another that inherits $\Delta\nu$ states, with $\Delta\nu_{\pm}=\pm 2$ and $\Delta\nu_{0}=0$. According to the bulk-boundary correspondence  (?, ?), such a spectral flow suggests a topological property of the bulk (Poincaré and Rossby) waves, obtained by assuming that the Coriolis , that encodes the impossibility to smoothly define each of these waves for any $k_{x}$, $k_{y}$ when the $f$ can change sign. This correspondence lies on the index theorem that relates the index of the bulk Hamiltonian both to the spectral flow and to a topological index  (?). More precisely, we expect that the spectral flow to the band $n$ is given by $\Delta\nu_{n}=\Delta\mathcal{C}_{n}$ where $\mathcal{C}_{n}$ is a topological invariant called the first Chern number assigned to the band $n$  (?, ?). In the following, we compute the Chern numbers of the Poincaré and Rossby modes and show that their value indeed satisfy the bulk-boundary correspondence, thus unveiling the topological origin of the Kelvin and Yanai equatorial waves.

We consider flows on a $f$-plane with $x$ in the zonal direction (positive towards the East) and $y$ the meridional direction (positive towards to North) [see Figure 1(a) in the Main Text]. The velocity field is denoted by $\mathbf{u}(\mathbf{x},t)=(u(\mathbf{x},t),v(\mathbf{x},t))$, with $\mathbf{x}=(x,y)$. The value of the Coriolis parameter $f$ is spatially constant, but will be considered as an external parameter that can be varied between $\pm 2\Omega$, where $\Omega$ is the planet rotation rate. The extremal values correspond to the value of $f$ at the Northern and Southern pole. The system possess one intrinsic time scale $1/2\Omega$ and one intrinsic length scale given by the global Rossby radius of deformation $L_{d}$ along with a speed $c$:

[TABLE]

The dynamics of Equations (1)-(2) can then be expressed in terms of dimensionless quantities

[TABLE]

and the linearized equations of motion around a state at rest ($\mathbf{u}=0$, $h=H$) read

[TABLE]

For convenience, we drop the $\tilde{}$ tilde in the following.

Eigenmodes (equivalently, normal modes) are propagating planar waves. Let $\hat{\eta}$, $\hat{u}$, $\hat{v}$ be the amplitudes of the fields ${\eta}$, ${u}$, ${v}$ for mode $e^{i{\omega}{t}-ik_{x}x-ik_{y}y}$. The amplitudes obey the eigenvalue equation:

[TABLE]

Diagonalization of this matrix leads to three eigenvalues corresponding to three wave bands in wavenumber-frequency space. The normalized eigenvectors associated with the positive frequency band ${\omega}_{+}=\left(k^{2}+{f}^{2}\right)^{1/2}$ are

[TABLE]

The ones associated with the zero frequency band $\omega_{0}=0$ are

[TABLE]

Eigenvectors for the negative frequency band ${\omega}_{-}=-\left(k^{2}+{f}^{2}\right)^{1/2}$ are easily obtained by taking advantage of the system symmetry: $\Psi_{-}(k_{x},k_{y},{f})=\Psi_{+}(-k_{x},-k_{y},-{f})$.

We calculate the Chern number of the fiber bundle defined by eigenvectors of the positive frequency wave band parameterized on the unit sphere in the space $\left(k_{x},k_{y},{f}\right)$ [see Figure 1(b) in the Main Text]. Importantly, normalized eigenvectors $\Psi_{0}$ and $\Psi_{\pm}$ are defined only up to a phase. Here we show that this phase cannot be defined continuously for $\Psi_{\pm}$ over all possible values of $\left(k_{x},k_{y},{f}\right)$. To see this, first notice that the eigenstates can be parametrized with the use of spherical coordinates. Introducing the colatitude $\theta$ and the longitude $\varphi$ (that should not be confused with colatitude and longitude on the planetary sphere in physical space), the normalized eigenvectors of the positive frequency wave band may be rewritten as:

[TABLE]

The eigenvectors are multivalued at $\theta=0$ and $\theta=\pi$, respectively the north and south poles of the parameter sphere $\mathcal{S}$ shown in Fig. 1(b) of the Main Text. The multivaluedness is a topological property of the vector bundle that is defined as the set of $\Psi_{+}$ over the base space $\left(\theta,\varphi\right)$. The multivaluedness can be quantified by non-vanishing first Chern number that we evaluate below.

Recalling that the eigenmodes are defined only up to a phase, a gauge transformation can be used to regularize the eigenmodes at one of the poles. In particular, we define

[TABLE]

so that $\Psi_{+}^{N}$ is a single-valued function of $(\theta,\varphi)$ everywhere but at the south pole. Likewise $\Psi_{+}^{S}$ is a single-valued function of $(\theta,\varphi)$ everywhere but at the north pole. The Berry connections associated with these vector bundles are

[TABLE]

with the inner product defined for the 3-component amplitudes as $\left<A,B\right>\equiv A_{1}^{*}B_{1}+A^{*}_{2}B_{2}+A_{3}^{*}B_{3}$ and $\nabla_{s}$ is the gradient operator on the unit sphere in $(k_{x},k_{y},f)$ parameter space. These Berry connections are related by the gauge transformation through

[TABLE]

At any point where $\mathbf{A}_{+}^{N}$ and $\mathbf{A}_{+}^{S}$ are not singular, the Berry curvature is given by

[TABLE]

The Chern number for the fiber bundle of eigenvectors with positive frequencies on the unit sphere in $(k_{x},k_{y},f)$-space is defined as the flux of the Berry curvature through $\mathcal{S}$

[TABLE]

where $\mathbf{e}_{k}(\theta,\varphi)$ is the vector normal to the sphere pointing in the outward direction. Stokes theorem can be applied separately for each hemisphere of $\mathcal{S}$ where the Berry connection is well defined, e.g. along the equator circle ($\theta=\pi/2$), with $\mathbf{B}_{+}=\nabla_{s}\times\mathbf{A}_{+}^{N}$ in the northern hemisphere and $\mathbf{B}_{+}=\nabla_{s}\times\mathbf{A}_{+}^{S}$ in the southern hemisphere. The calculation yields

[TABLE]

Using Eq. (S12), we find $\Delta\mathcal{C}_{+}=2$. This Chern number can be interpreted as the topological invariant of the transition from one hemisphere to another for the case of the positive frequency eigenmodes of the shallow water model in the $f$-plane approximation, with $f$ taken as an external parameter that tunes this transition. Note that we have derived this result by considering the fiber bundle of eigenmodes parameterized on the unit sphere, but it is straightforward to show that it holds on any sphere of radius $\omega_{+}>0$ in $(k_{x},k_{y},f)$ space since the normalized eigenvectors $\Psi_{+}$ do not depend on $\omega_{+}$. Similar calculations show that the Chern number for the negative frequency eigenmodes is $\Delta\mathcal{C}_{-}=-2$, and $\Delta\mathcal{C}_{0}=0$ for the zero-frequency eigenmodes.

According to bulk-edge correspondence  (?, ?), the number of topological edge states at the interface between two domains of opposite $f$ is fixed by the set of Chern numbers $\left\{\Delta\mathcal{C}_{-},\ \Delta\mathcal{C}_{0},\ \Delta\mathcal{C}_{+}\right\}$. Let $N^{\text{above}}_{\text{right}}-N^{\text{above}}_{\text{left}}$ be the difference in the number of right-moving and left-moving modes in the gap above the central geostrophic band, and $N^{\text{below}}_{\text{right}}-N^{\text{below}}_{\text{left}}$ be the difference in the number of right-moving and left-moving modes in the gap below. Bulk-edge correspondence then states that $|\Delta\mathcal{C}|=\left|\left(N^{\text{above}}_{\text{right}}-N^{\text{above}}_{\text{left}}\right)-\left(N^{\text{below}}_{\text{right}}-N^{\text{below}}_{\text{left}}\right)\right|$. The numbers can vary with frequency within the gap, but the difference $|\Delta\mathcal{C}|$ is the number of topologically protected chiral edge modes. In the present case, Chern number $\Delta\mathcal{C}_{0}=0$ means that there are equal numbers of topological edge states with the same chirality in the gaps above and below the central band of geostrophic modes. Since there is no gap above the positive frequency band, and no gap below the negative frequency band, the number of chiral edge states in each gap surrounding the central band is given by $\left|\Delta\mathcal{C}_{\pm}\right|=2$.

2 Numerical Calculation of Eigenmodes With $f(y)$: From the Shallow Water Model to a Haldane-like Model on the Lieb lattice

The linearized rotating shallow water model on a doubly periodic domain with $f$ varying in the $y$ direction is diagonalized numerically using a standard C-grid discretization procedure  (?). This approach establishes a direct link to lattice models used in condensed matter physics.

The domain is discretized into a uniform grid with $N_{x}\times 2N_{y}$ unit cells located at $(x,y)=(m\,\delta x,n\,\delta y)$ where $\delta x$ and $\delta y$ are the size of the grid cells, and $m$ and $n$ are integers. Assuming time-harmonic solutions, eigenmodes of the discretized shallow water system are described at each grid point by the real part of the triplet $e^{{\rm i}\omega t}(h_{m,n},\mathbf{u}_{m,n})$. Following the C-grid discretization procedure, the system of equations (S3)-(S4)-(S5) then turns into an eigenvalue problem for the fields $(h(m,n,\omega),\mathbf{u}(m,n,\omega))$:

[TABLE]

where the change of variable ${\rm i}h\rightarrow h$ was performed for convenience, and where we have defined $t_{x}=\frac{1}{\delta x}$ and $t_{y}=\frac{1}{\delta y}$. We consider a meridional dependence of the Coriolis parameter through the map $f_{y}$ that we specify below. This allows us to treat arbitrary variations of $f$, and in particular any equator-like interface from $f<0$ to $f>0$. For concreteness, we consider periodic boundary conditions in the $y$ direction with $n$ running from $1$ to $2N_{y}$, and we choose

[TABLE]

where $R_{f}$ gives the scale of variation for the gradient of the Coriolis parameter in the $y$ direction. There are two regions with opposite $f$ and two interfaces between the regions (see Figure S1). The doubly periodic geometry has the advantage of simplifying numerical computations. The price to pay is the existence of tow equators.

The frequency $\omega$ is the eigenvalue of a tight-binding model on the Lieb lattice shown in Figure S2. Each site is coupled to its nearest neighbors by hopping parameters $\pm t_{x}$ and $\pm t_{y}$. The minus sign that appears in these matrix elements can be interpreted as a phase of $\pi$ picked up when hopping from site to site. While this phase does not break time-reversal symmetry, the coupling $\pm{\rm i}f/4$ that accounts for a phase of $\pm\pi/2$ accumulated between second nearest neighbors, does break it. The total phase accumulated around a unit cell is a multiple of $2\pi$, so that the net flux through the system vanishes. This behavior is analogous to the seminal Haldane model on the honeycomb lattice, the archetypal example of a topological Chern insulator in the absence of an applied magnetic field  (?).

Because of translational invariance, the Fourier transformed fields $(h(k_{x},n,\omega),\mathbf{u}(k_{x},n,\omega))$ satisfy

[TABLE]

There are four dimensionless parameters given by $R_{f}/L_{y}$ $L_{d}/L_{y}$, $L_{x}/L_{y}$, and $\delta y/L_{y}$. The first parameter controls the interface thickness between regions of constant (but opposite) $f$. It is of order $1$ for a spherical geometry, but the difference between non-topological and topological edge modes is enhanced by considering a sharper interface, with $R_{f}/L_{y}\approx 0.1$. The second non-dimensional parameter, $L_{d}/L_{y}$, controls the intrinsic length scale of the domain in comparison to the domain size. It is of order one in the atmosphere, and of order $0.1$ in the ocean. The third parameter $L_{x}/L_{y}$ is a geometric parameter, of order one for the atmosphere and oceans on Earth. For convenience, we consider $L_{x}/L_{y}\gg 1$, which yields continuous lines $\omega(k_{x})$ in the dispersion relation. The last parameter is related to the discretization of the continuous problem and is chosen such that $\delta y/L_{y}\ll R_{f}/L_{y}$ and $\delta y/L_{y}\ll L_{d}/L_{y}$. The precise form of the dispersion relation for $k_{x}\sim\pm 2\pi/\delta_{y}$ depends on the choice of the discrete geometry. However, for a given $k_{x}$, the dispersion relation of the continuous shallow water system is recovered in the limit $\delta x\rightarrow 0$.

The system (S18) is diagonalized numerically. The frequency spectrum shown in Fig. S3 is obtained for $N_{y}=50$, $R_{f}/\delta y=10$ and Rossby radius of deformation $L_{d}=R_{f}/2$. We have checked that regardless of the value of these dimensionless parameters, the gaps obtained for $f\neq 0$ are filled by two equatorially trapped modes, and that for each equator, the difference in the number of right-moving and left-moving edge modes within the frequency gaps is equal to $\pm 2$: in the earth-like equator ($\partial_{y}f>0$), these unidirectional modes have eastward group velocity. By contrast, the amplitude of the Rossby waves, and the emergence of equatorially trapped Poincaré waves, depend on the shape of $f(y)$.

3 Scattering of Equatorial Kelvin Waves by a Static Perturbation

Equatorial Kelvin waves were predicted in Refs.  (?) and subsequently observed in Earth’s oceans  (?) and atmosphere  (?). The waves are reproduced by models more general than the single-layer shallow water equations considered in the Main Text  (?). Conversion of low-frequency eastward-moving Kelvin waves to westward Rossby waves can occur at the eastern margins of ocean basins  (?) or by scattering from topography  (?, ?). Here we investigate the frequency dependence of the scattering of equatorial Kelvin waves by a region of enhanced stratification in a 2-layer primitive equation modelto demonstrate that Kelvin waves that oscillates at an intermediate frequency in the gap between the Rossby and Poincaré waves are not backscattered by the perturbation.

The model is defined on a spherical geodesic grid as described in Ref.  (?). This minimal model of the general circulation describes fluid motion in the horizontal within each layer, flow between the layers, and the variation of temperature within each layer. The model has rigid upper and lower lids and supports internal waves at the interface between the two layers. The temperature field is advected by the fluid and acts back upon it via pressure as determined by the equation of state for an ideal gas and hydrostatic balance. There are in total 5 independent dynamical fields once the incompressibility of the fluid has been taken into account. The idealized model can be used to describe either a planetary atmosphere or an ocean that covers the whole planet. Code that implements the model is publicly available111The application “GCM” is available for macOS 10.9 and higher on the Apple Mac App Store at URL http://appstore.com/mac/gcm.

The vertical mean temperature is initially uniform over the sphere, and can be set to zero without loss of generality. Stable stratification is imposed by setting the temperature of the upper layer to be $1$ K greater than that of the lower layer, except in a region centered on the equator at longitude $90^{\circ}$ W (indicated by the dashed circles in Figure S4) where the temperature difference increases up to a maximum of $10$ K further with a Gaussian profile of width $20^{\circ}$. The lower layer is thus denser than the upper layer and supports internal Kelvin and Rossby waves that move at speeds comparable to those observed in the equatorial Pacific ocean. Waves can scatter off the region of enhanced stability because the wave speed is higher in that region. Kelvin waves in the ocean can be observed from satellite via surface altimetry as water that warms also expands. The imposed model stratification means that the simulated vertical mean temperature likewise serves as a proxy for the interface height.

Wavepackets are generated early in the simulations by adding a forcing term to the tendency of the difference $\delta$ in the divergence of the velocity field between the two layers

[TABLE]

where $\ldots$ denotes the other terms in the equation of motion (Equation 5 of Ref.  (?)) and $g_{\Delta}(\alpha)\equiv\exp[-(\alpha/\Delta)^{2}]$ is a Gaussian of unit height centered at $\alpha=0$ and of width $\Delta$. The added driving acts for only half of an oscillation period, from $t=0$ to $t=\pi/\omega$. We choose $A=10^{-4}$ (day)-1, an amplitude sufficiently small for the waves to behave linearly. Apart from the wave driving term at the early stage of simulation, all other forcing terms and friction are are turned off for the numerical experiments performed here.

An intermediate frequency wave train that consists largely of an eastward moving Kelvin wave is produced by setting the angular frequency $\omega=1$ (day)-1 and wavenumber $m=10$ [Figure S4(a)]. A low frequency, long wavelength Kelvin wave is created by choosing instead $\omega=0.2$ (day)-1 and $m=2$ [Figure S4(d)]. (Yanai waves can be generated by modifying the forcing term in Eq. S19 to be antisymmetric with respect to the equator.)222Videos of the two simulations may be viewed at URLs https://vimeo.com/204597701 and https://vimeo.com/204598091.

The localized region of enhanced stability breaks translational symmetry; therefore wavenumber is not conserved. Because the perturbation is static in time, however, the frequencies of the waves passing through the region remain unchanged. The intermediate frequency Kelvin wave packet propagates at a frequency that lies in the gap between the Rossby and Poincaré waves as indicated in Figure 1a in the Main Text. There are thus no available scattering channels that match its frequency, and it emerges from the interaction intact [Figure S4(b, c)]. (Reflection symmetry of the Kelvin wave about the equator blocks scattering into the Yanai wave.) Because the waves in the frequency gap have a topological origin, the absence of scattering is a manifestation of topological protection. By contrast there is a band of Rossby waves at low frequencies, and the low-frequency Kelvin wave is strongly disturbed by its interaction with the enhanced stratification [Figure S4(e)], producing Rossby waves. At late times, the low-frequency Kelvin wave has largely been converted to Rossby waves moving in the opposite (westward) direction [Figure S4(f)].