From ray tracing to waves of topological origin in continuous media

Inhomogeneous media commonly support a discrete number of wave modes that are trapped along interfaces defined by spatially varying parameters. When they are robust against continuous deformations of parameters, such waves are said to be of topological origin. It has been realized over the last decades that such waves of topological origin can be predicted by computing a single topological invariant, the first Chern number, in a dual bulk wave problem that is much simpler to solve than the original wave equation involving spatially varying coefficients. The correspondence between the simple bulk problem and the more complicated interface problem is usually justified by invoking an abstract index theorem. Here, by applying ray tracing machinery to the paradigmatic example of equatorial shallow water waves, we propose a physical interpretation of this correspondence. We first compute ray trajectories in a phase space given by position and wavenumber of the wave packet, using Wigner-Weyl transforms. We then apply a quantization condition to describe the spectral properties of the original wave operator. We show that the Chern number emerges naturally from this quantization relation.

eigenvectors of the linear operator.When the coefficients of the partial differential equations (PDE) are homogeneous (constant in space), solutions are easily found in unbounded geometries using a Fourier transform.Such bulk wave problems are simple to compute, and are useful to understand many aspects of wave propagation.For instance the computation of this bulk problem allows to classify waves in different wave bands such as sound waves, gravity waves, among others.Yet, most practical situations involve spatially varying parameters, which are in general intractable analytically.
There are, however, three powerful complementary theoretical approaches that can be used to tackle wave problems in nonuniform media with smooth variations: the study of simple but emblematic normal forms, the use of ray tracing equations, and the topological analysis: • The first normal form approach is based on the idea that most salient features associated with spatial variations in the medium can be understood qualitatively by assuming linear variations of the parameters in one direction but constant in the other.A great success of this approach was the discovery of two unidirectional equatorially trapped waves by Matsuno in 1966 [1].The drawback of this approach is its lack of generality, and there is no guarantee that the problem will be solvable analytically, although numerical computation of the spectrum remains always a possibility.• The second ray tracing approach consists in computing the trajectories of wave packets, by assuming a scale separation between their wavelength and the spatial variations of the medium.This scale separation makes possible the use of standard asymptotic technique, such as the Wentzel-Kramers-Brillouin (WKB) approximation.This approach has long been used in fluids, and has even led to concrete applications in ocean dynamics [2,3], and atmospheric dynamics [4].• The third topological approach is more recent.The main idea is that global features of the dispersion relation in inhomogeneous media can be classified and predicted.In particular, when a spatially varying parameter defines an interface, the wave spectrum may exhibit modes that are trapped along the interface, and that fill the frequency gap between different wave bands.Those modes are often referred to as topological waves, which is a shorthand term to mean waves of topological origin, when their presence is robust to continuous changes in parameters.It turns out that the emergence of such modes in interface wave problems can be predicted by computing a topological invariant for a set of bulk wave problems that are much simpler to solve.The topological invariant is an integer, the Chern number, that predicts the number of trapped modes in interface wave problems.Born in condensed matter, those ideas have irrigated all fields of physics, including now astrophysical and geophysical flows.For instance, a Chern number 2 was computed for rotating shallow water waves, consistently with the presence of two unidirectional modes in the Matsuno spectrum [5].
The topological methods may be thought of as a way to justify the outcomes of normal form approaches to understand more complicated situations.For instance, the computation of a topological invariant related to the Matsuno wave problem guarantees that the presence of two unidirectional modes trapped at the equator is a robust feature of equatorial waves, independently from the details of the planet's curvature, or from continuous changes of other parameters of the problem.In most previous applications of topology to geophysical and astrophysical flows, and more generally in continuous media, the correspondence between the topological invariant of the bulk problem and the number of topological waves (more precisely defined as a spectral flow later on) in the interface wave problem was justified by invoking an abstract index theorem [6][7][8].The aim of these notes is to establish a connection between the ray tracing approach and the topological approach, in order to provide some physical intuition on this bulk-interface correspondence, building upon previous works on equatorial waves.
More precisely, we intend to relate the seminal work of Littlejohn and Flynn [9] on ray tracing equations for multicomponent wave systems, to more recent works on topological waves in continuous media, spectral flows and Berry monopole as presented for instance in the lecture notes of Faure [7] or Delplace [8].In both cases, a fundamental role is played by the Berry curvature, which can be computed from the knowledge of the bulk wave polarization relations.It was shown that ray trajectories in position/wavenumber phase space are modified by this term just as a charged particle is deviated by a magnetic field [9].This has found application from condensed matter [10][11][12] to geophysical fluid dynamics [13].
In the context of topological waves, the same Berry curvature (a geometrical quantity) is used to compute the Chern number (a topological quantity).This Chern number counts singularities in an abstract space of bulk wave polarization relations parameterized over a closed surface.The Chern number is related to the Berry curvature through a generalization of the Gauss-Bonnet formula that relates the genus of a surface to the integral of the Gaussian curvature over this surface [7,8].It is thus tempting to establish a connection between ray tracing, which describes the dynamics of a local wave packet influenced by the Berry curvature, to the global spectrum of a wave operator with spatially varying coefficient, which is ruled by a single Chern number.
The duality between ray tracing and spectral properties of an operator is reminiscent of the duality between classical and quantum mechanics.Classical mechanics deals with ordinary differential equations (ODE) that describe trajectories in phase space, while quantum mechanics deals with PDE that describe wave functions.Originally, quantum mechanics was derived by proposing quantization procedures that map functions of phase space variables such as the classical Hamiltonian to operators of quantum mechanics such as the Schrödinger operator [14].Semi-classical analysis is the reverse procedure that derives classical mechanics trajectories in phase space (and possible corrections) in the limit of vanishing Planck constant, as done for instance by using the celebrated WKB ansatz [15].Those techniques are not restricted to quantum mechanics, and a branch of mathematics called microlocal analysis was developed in the seventies to establish connections between ray tracing and spectral properties of operators, see e.g.[16].
The connection between ray tracing and spectral properties of an operator makes use of the Wigner-Weyl transform, that relates in a systematic way a wave operator to its symbol, which is a function playing for instance the role of the Hamiltonian in classical mechanics.Such tools have been used for instance recently by mathematicians to explain the emergence of singular wave patterns in density-stratified flows [17], or the dynamics of equatorial waves [18][19][20], including through the use of topological tools [7].This Weyl-Wigner formalism has also been used by theoretical physicists to describe wave transport phenomena in continuous media [9], including geophysical flows [21,22], see also the recent work by Onuki for a pedagogical introduction to the formalism in the context of geophysical flows [23].Our contribution will be to connect those works with previous studies on topological waves in continuous media such as geophysical flows [5,24], astrophysical flows [25,26], or plasmas [27,28].
We present in section 1 the equatorial shallow water model [1], and use this example to introduce the notion of spectral flow index [7,8], which is an integer that counts the number of topological modes that transit from one wave band to another in the spectrum of a wave operator, when a parameter is varied from −∞ to +∞.Those two limiting cases allow us to define the equivalent of a semi-classical limit.
We review in section 2 the standard procedure to diagonalize a multicomponent wave operator, taking advantage of a small parameter in the problem in the semi-classical limit, building upon [9,12,23].This derivation makes use of the Wigner-Weyl transform and notions of symbolic calculus that are presented in Appendix A. The diagonalization of the multicomponent wave operator allows us to define a scalar operator describing the dynamics of a wave packet in a given wave band.
We derive in section 3 the dynamics of those wave packets, which allows us to describe ray trajectories in the phase space defined by position and momentum of the wave packets.Those results were first derived by [9], who highlighted the central role of the Berry curvature, and a duality between two different possible definitions of position/wavenumber of the wave packets, one of them being independent from any gauge choice done during the diagonalization of the original wave operator.Our contribution is to propose a new physical interpretation of the transformation between two different phase space coordinates for the wave packets, and to explain that the Berry curvature appearing in the ray tracing equations is related to the presence of a topological invariant in the bulk wave problem: the first Chern number.
We recover in section 4 the spectral properties of the wave operator by applying a quantization condition (imposing that the phase picked up by a wave packet along a trajectory is a multiple of 2π), as proposed in [9].Our contribution is to show how this quantization can be used to recover the spectral flow result noticed in section 1.The derivation of the quantization condition will highlight terms than can be interpreted as integral of the Berry curvature over a closed surface, which is nothing but the Chern number, a topological invariant.Thus, we show how spectral properties are related to topological properties through ray tracing.We finally propose an alternative interpretation of the spectral flow result based on the modification of the phase space density of states by the Berry curvature, a standard result in non-canonical Hamiltonian systems.
Notations Any operator will involve a hat, as Ω.Operators or functions that are underlined once means that they involve multiple components.Operators underlined twice are matrices of operators.A character underlined twice but without a hat is a matrix of functions.Bold symbols such as z are vectors in phase space.A tilde is placed on a dimensional variable and it will be removed to define a scaled dimensionless one.

Motivation from equatorial shallow water waves 1.Wave equation and definition of the spectral flow
The shallow water model describes the dynamics of a thin layer of fluid with homogeneous density under gravity (constant g), with a solid boundary at the bottom and a free surface at the top [29].For the sake of simplicity, we assume that the dynamics takes place in a plane tangent to Earth near the equator, over a flat bottom, and the fluid layer thickness is denoted by H (see Fig. 1).
In this geometry, the dynamics is conveniently expressed in Cartesian coordinates (x, y), where x points eastward (zonal direction) and y points northward (meridional direction).Metric terms are neglected, but planet's sphericity is taken into account by considering linear variations of the Coriolis parameter in the meridional direction y.This Coriolis parameter denoted by βy is twice the projection of the planet's rotation vector on the local vertical axis, as in the Foucault pendulum experiment.The linearized shallow water dynamics around a state of rest reads where (u, v) is a depth-independent two-dimensional velocity field and η is the variation of surface elevation with respect to the rest state.1A one-dimensional wave problem can be derived from this model.We take advantage of the invariance with x to write (ũ, ṽ, η) = (cu, c v, Hη) e ik x , c = g H . ( From now on we chose β = 1 and c = 1, which can always be done by adimensionalizing time and length units with 1/ β c and c/β, respectively.The wave problem is now The operator Ĥk depends on a single parameter k.This problem was fully solved by Matsuno, and the corresponding dispersion relation is displayed in Fig. (2).One can look at how the wave spectrum changes when this parameter is varied: When varying the parameter k from −∞ to +∞, two modes are gained by the positivefrequency wave band, called the Poincaré wave band.
Similarly, two modes are lost by the negative frequency Poincaré wave band.By contrast, the central wave band called the Rossby wave band has a net gain of modes equal to 0: there are as many branches that leave the Rossby wave band as branches that enter the Rossby wave band when k is varied from −∞ to +∞.
The bottom line is that some of the branches in the dispersion relation transit from one wave band to another when k is varied, which defines a spectral flow, 2 and k is called the spectral flow parameter.The spectral flow is quantified by an integer, the spectral flow index, that counts how many modes (branches in the dispersion relation) are gained or lost as k is varied from −∞ to +∞.More precise and rigorous definitions of this index are given for instance in [7], and its relation with another integer named analytical index is explained in [8].
Matsuno's computation leads to a spectral flow index +2 associated with positive-frequency Poincaré wave band, 0 for the Rossby wave band, and −2 for the negative frequency Poincaré wave band.
Note that in many cases, the spectral flow is directly related to the difference between the number of right-moving (positive group velocity) and left-moving modes in a given range of frequencies, most often taken within the gap between different bulk wavebands when it exists.In Matsuno's case, at any nonzero frequency, there are two more modes with eastward group velocity than modes with westward group velocity, and this is related to the spectral flow index ±2 of the upper and low wave bands [5].Thus, since group velocity corresponds to energy transport, spectral flow is related to unidirectional wave transport in that case.

The semi-classical limit for equatorial shallow water waves
Our aim is to understand this spectral flow from a semi-classical perspective, which is related to the familiar WKB approximation commonly used in geophysical fluid dynamics.For that purpose we take advantage of the existence of a small parameter in the limit k → ±∞, by rescaling the wave equation as follows: which leads to The parameter ε appears only in front of the y-derivative, which is suggestive of the reduced Planck constant ħ h in quantum mechanics.This is why the limit ε → 0 will from now on be referred to as a semi-classical limit, and ε will be called the semi-classical parameter. 3Given that the length unit used to adimensionalize equations is c/β, the semi-classical parameter compares this intrinsic length scale, called the equatorial radius of deformation, to the zonal wavelength 1/k.Since c/β corresponds to a trapping length scale in the meridional direction, the semi-classical limit is a limit of large meridional to zonal aspect ratio for the waves.The upper-script index ± is used to distinguish the limit k → +∞ from the limit k → −∞, which both correspond to the semi-classical scaling ε → 0. An important remark follows: Understanding the spectral flow can be tackled by comparing the spectral properties of the operators Ĥ+ and Ĥ− in the semi-classical limit ε → 0.More precisely, we will pair together modes of the two operators that share common properties.By continuity of the eigenvalues, the modes that can not be paired will be those belonging to a branch that transits from one wave band to another as k varies.

Strategy to interpret the spectral flow
We will show in section 2 how to project the original multi-component wave operator Ĥ± into scalar operators for each wave band, denoted by Ω(j)± , where j is the band index.Those scalar operators are different for the three wave bands, just as their dispersion relation are different.In fact, we will highlight a formal duality between such operators and the corresponding dispersion relations.
For that purpose, we will follow [9] and [23] and rely on the Wigner-Weyl transform that is introduced in Appendix A.1.After this first step, the spectral flow index for each wave band will be understood loosely as the difference between two infinite numbers N ( j) = # of modes for Ω(j)+ − # of modes for Ω(j)− , (6) where # means "number".We will derive in section 3 the ray tracing equation in wavenumber/physical space (phase space), for wave packets described by operators Ω(j)± .We will highlight the central role played by a quantity called the Berry curvature in ray tracing equations, as noted first in [9].We will explain how to compute this Berry curvature from the knowledge of the wave polarization relations.Using an integration over a surface in parameter space, we will relate this curvature, which is a local geometrical quantity, to a singularity described by a topological invariant named the first Chern number, which describes global properties of continuous families of eigenvectors.
We will show in section 4 how to recover the spectrum of Ω(j)± from the use of a quantization procedure of ray trajectories in phase space, following [9].The quantization procedure amounts to assuming that the phase picked-up by a wave packet along a trajectory is an integer multiple of 2π.This is analogous to the Bohr-Sommerfeld quantization in quantum mechanics.This procedure will make possible an explicit computation of the mode imbalance (6).This will also highlight the topological origin of this mode imbalance, that is related to an integral of the Berry curvature over a closed surface, and thus to the Chern number.
The index ( j) for the wave bands will be dropped in the following when discussing general properties of the operator Ω(j)± .In many steps of the derivations, the formula will be sufficiently general to apply for a variety of flow models.However, for pedagogical reasons we chose to focus on the particular case of shallow water waves.
2 Scalar wave operator and symbols

From multicomponent to scalar wave operator
The initial multicomponent wave problem is written formally as where Ψ is a multicomponent wave field and Ĥ a linear operator involving spatially varying coefficients and spatial derivatives, together with an external parameter k to be varied.In the shallow water case, this spectral flow parameter k is the wavenumber in zonal (West-East) direction.Here we will focus only on 1D problems involving the spatial (meridional) coordinate y, with the corresponding derivative ∂ y .The parameter to be varied will be the wavenumber k in the x (zonal) direction.An example is given in the previous section, where Ĥ is to be replaced by Ĥ± .
We recall here general results that are presented in a pedagogical way by Onuki [23] (for fluid dynamicists) and Reijnders et al [12] (for condensed matter physicists).We would like to express the multicomponent wave dynamics on the form of a scalar equation for a single wavefield denoted by ψ, assumed to be written on the following form: We also assume that reconstruction of the multicomponent wavefield from the scalar field is formally expressed with a multicomponent operator χ (a vector whose each component is an operator depending on the y coordinate): where 1 is the identity operator.The second equality guarantees that both Ψ and ψ are normalized: In the case of shallow water waves, this norm represents the total energy of the flow, 4 which is the sum of kinetic plus available potential energy.At this stage the operators Ω and χ are not known.A useful equation relating those operators is obtained by combining ( 7), (8), and (9): We will see in the next section that the expression for Ω and χ can be obtained in the semiclassical limit, and are different for each wave bands of the original multi-component wave problem.

Expression of the scalar operator in the semi-classical limit
Our aim here is to make an educated use of Wigner transform, symbolic calculus, and Weyl transform to find the expression of the scalar operator Ω, by exploiting (11).Those tools are introduced in Appendix A. If an operator such as Ω is known, then the Wigner transform defined by (A.2) yields its symbol Ω( y, p).This symbol is a function of y and p, the conjugate momentum in the direction y, which is a scaled wavenumber.This procedure generalizes the Fourier transform to problems with spatially varying coefficients.As such, the symbol Ω( y, p) can be interpreted as a local dispersion relation for a plane wave oscillating much faster than the spatial variations of the parameter in the medium.The inverse transform that builds an operator Ω from the knowledge of a function Ω( y, p) is called the Weyl transform, and is defined in Appendix A.1, equation (A.3).The operator Ω is called a pseudo-differential operator as it cannot always be expressed as an explicit polynomial in ∂ y [30].
Formally, all the operators in ( 11) can be expanded as In the following, removing the hat (ˆ) from an operator will refer to its symbol.Symbols can also be expanded in the semi-classical limit as ) Products of operators âb are in general different from the operator a b obtained by a Weyl transform of the standard product between symbols a( y, p) and b( y, p).As explained in Appendix A.3, equation (A.15), we can however define a particular symbol product denoted by a ⋆ b such that the Weyl transform of this star product corresponds to the standard product of operators (i.e.their composition).Let us apply this procedure to develop the operator products in (11).The r.h.s. is expanded as The star product is computed by using its definition (A.15) together with the symbol expansions ( 16) and ( 17) in the semi-classical limit, which yields where we have introduced the Poisson bracket for two symbols χ 0 ( y, p) and Ω 0 ( y, p): The corresponding operator is finally obtained by applying the Weyl transform (A.3) to the star product.

Asymptotic expansion up to order one
Using symbolic calculus and the asymptotic expansion introduced in section 2.2, one can now exploit the equality (11) to find expressions of χ 0 , Ω 0 and Ω 1 .At the leading order, we obtain a matrix eigenvalue equation for H 0 : We see that the symbols Ω 0 and χ 0 are just the eigenvalue and the eigenvector of H 0 , the leading order component of the symbol associated with the multicomponent wave operator.This zeroth order result may be understood as the outcome of a "local" plane wave solution: in this framework, χ 0 is the local polarization vector associated with the local dispersion relation Ω 0 ( y, p).
Collecting the first order terms in (11), multiplying on the left by χ † 0 , and using owing to the hermiticity of the wave operator5 ) leads to ) More detailed on this computation can be found for instance in [9,12].
Let us now explain why we have separated the first order expression into two compontents Ω 1A and Ω 1B .The important point noticed by Littlejohn and Flynn is that the eigenvectors χ 0 of the zeroth order equation are defined up to a phase factor [9].In physical jargon, choosing this phase amounts to a gauge choice.It appears that the first order expression of the scalar symbol (and the corresponding operator Ω1 ) depends on this gauge choice.It is however possible to split the symbol into a part Ω 1A that is gauge independent and a part Ω 1B that is not gauge independent.This can be checked by applying the transformation χ 0 → χ 0 e i g( y,p) , where g( y, p) an arbitrary real-valued function.The term Ω 1A is left unchanged, while the term Ω 1B is shifted as Ω 1B → Ω 1B + {g, Ω 0 }.As we shall see later, the term Ω 1B is related to the Berry curvature describing local variations of the eigenvectors χ 0 in parameter space (i.e.local variations of the polarization relation in a given wave band).As such, the term Ω 1B is sometimes called the Berry term.

Application to the shallow water wave problem
Let us come back to the shallow water wave problem introduced in (5).The symbol of the wave operator Ĥ± is of order 0 in ε: The three eigenvalues of the symbol are, at leading order: Notice that the eigenvalues of H + 0 and H − 0 are the same.The difference between the spectra of both operators manifests at next order (equations ( 23)-( 24)), for which the expression of the eigenvectors at zeroth order is needed: Using those expressions in ( 23)-( 24) yields to the first order corrections: The correction Ω 1A for the middle wave band is known as the dispersion relation of Rossby wave.The Berry correction Ω 1B is zero in that case, which is consistent with ( 24), together with Ω 0 = 0 for this wave band.As explained previously, the correction Ω 1B for the other bands depends on the phase choice made in (27) for χ 0 .
Note on the literature.The term Ω 1A will play a crucial role in ray tracing equations, in the next section.In that context, it was computed by N. Perez and called gradient correction, see (B15) in [13].Using the same gauge choice, Y. Onuki obtained the expression of Ω 1 as a particular case of a more general computing including horizontal variations of the layer thickness (equation (3.9) in [23]), without splitting the expression into a Berry part and a gauge-independent part.Gallagher and coworkers also obtained similar expressions in a series of paper on the mathematics of equatorial waves [18], including more general configurations with a prescribed horizontal mean flow field [19,20].

Wave packet dynamics in the semi-classical regime
This section builds upon previous works on the manifestation of the Berry curvature in ray tracing as reviewed by Niu's group for electronic waves [10], Perez et al for geophysical waves [13], see also the pioneering work of Littlejohn and Flynn [9].The only novelty here is a physical interpretation of the formal change of variable appearing in the paper by Littlejohn and Flynn, which allows us to make connections with Onuki's recent work on wave transport in geophysics [23].

Wave packet center of mass and wave packet momentum
We define the wave packet's location y v and wavenumber p v as the corresponding operators averaged with the energy density Ψ † • Ψ as weight, since the multicomponent wave field Ψ is the physical field to be observed: where the subscript v stands for vectorial, and where Ψ is normalized according to (10).In the quantum mechanical context, y v and p v are just the expectation values of the position and momentum operators.Recall that the normalization constraint ( 10) is equivalent to the energy conservation for shallow water waves.This is why y v and p v are respectively interpreted in this context as an averaged energy-weighted position and momentum (wavenumber) for a given wavepacket.The weight |Ψ| 2 corresponds indeed to the sum of local kinetic and available potential energy, which, in dimensional units, is 0.5(Hu 2 + H v 2 + gη 2 ).
It will be useful to introduce similar quantities defined formally from the scalar wavefield ψ, although their physical interpretation is less straightforward: where the subscript s stands for scalar.The quantities ( y s , p s ) in ( 31) are averaged position and momentum, just as ( y v , p v ) in ( 30), albeit with a different weight.In the case of ( y v , p v ), the weight was the local energy density, a physically meaningful quantity.This is not the case for ( y s , p s ).In fact, we will see in subsection 3.3 that ( y s , p s ) are not physical observables for the wave-packet since they depend on the gauge choice for the wave vector reconstruction operator χ that relates the scalar field ψ to the vector field Ψ through (9).The two different definitions ( 30) and ( 31) for an averaged position and momentum of the wave packet can be related together.As we shall see, they are not equivalent, due to the non-commutation between the multicomponent projection operator χ and the operators y and using the commutation relations (A.12-A.13)established in Appendix A.3, together with the Poisson bracket definition in (20), we find From ( 30) to (33), no approximation is involved, as no assumption on the form of Ψ or ψ is needed.

WKB ansatz for the scalar wavefield
We now consider the traditional WKB ansatz for the scalar wave field ψ: where a 0 ( y, t), φ 0 ( y, t) and φ 1 ( y, t) are real fields of order ∼ 1.The normalization condition (10) for ψ leads to the constraint The scalar wave field ψ is related to the multicomponent wave field Ψ through (9).Using the order zero expansion of operators acting on the WKB ansatz (as detailed in Appendix A.5) leads to Using this WKB ansatz, the operators in (33) can be replaced by the expression of their symbols at leading order, in order to obtain an approximation at order ε: where the symbol χ 0 and the Poisson bracket are evaluated at point ( y, p), with p( y, t) given by (37).

Assumption of a localized wave packet
We now assume that the wave packet is localized at y v and extends over a scale ∆ y ≪ 1, keeping ∆ y ≫ ε.Using (35), equation ( 38) is further simplified into where the terms in the r.h.s. are evaluated at z s .The expression remains actually correct when evaluated at z v , but corrections are then of order O(ε∆ y).From now on, we remove those correction terms.It is worth insisting on the importance of (39) by writing independently the two components, and by expanding the Poisson bracket: The correction terms involve the components of a vector called the Berry connection: Here, the Berry connection is a measure of how the polarization vector χ 0 ( y s , p s ) varies in phase space ( y s , p s ).This term is gauge dependent: it is not invariant for a change of phase choice in χ 0 .Since y v and p v are quantities built from the initial multicomponent wave problem that does not depend on any phase choice, they are gauge invariant.Thus, the breaking of gauge invariance by the Berry connection implies that the coordinates ( y s , p s ) based on the scalar wave equation are not gauge-invariant.Those observations will be confirmed by the ray tracing equations obtained in both choices of coordinates.
Equations ( 40)-( 41) correspond to the change of variable proposed by Littlejohn and Flynn to obtain gauge-independent phase space coordinates.[9].We showed here that this change of variable has a wave packet interpretation. 6

Gauge dependent, canonical Hamiltonian form
A time differentiation of (31), together with the Hermiticity of the operator Ω and the commutation rule (A.12), leads to Using the WKB ansatz in the limit ε → 0 followed by the localized wave packet assumption leads to the canonical ray tracing equations, where the frequency Ω 0 + εΩ 1 plays the role of the Hamiltonian, where Ω s is the symbol of the scalar operator.An unpleasant situation occurs: the ray dynamics in phase space ( y s , p s ) depends on the term Ω 1B which is gauge-dependent, as Ω 1B depends on the phase choice for χ 0 .

Gauge independent, non-canonical Hamiltonian form
Using the change of variables ( 40)-( 41) in ray tracing equations ( 44)-( 45)-( 46) leads to a new set of equations: where ) is called the Berry curvature: The Berry curvature is the curl of the Berry connection A( y, p) = (A y , A p ) introduced in (42): By contrast with the Berry connection, the Berry curvature is gauge-independent, as it is left unchanged by a change of phase choice for χ 0 .In fact, all terms in the ray tracing equations ( 47)-( 48)-( 49) are gauge-independent.In particular, the term Ω 1B present in (46) has been cancelled out in (49) by a contribution induced by the change of variable from ( y s , p s ) to ( y v , p v ).This was expected as we already noticed that y v and p v are physical observable interpreted as averaged position and momentum with a local energy density weight, and the temporal evolution of such physical observables cannot be gauge-dependent.The price to pay for being gauge-independent in this new set of coordinates is that the presence of Berry curvature renders the Hamiltonian dynamics non-canonical [9].Indeed, the system of equations (47)-( 48) is Hamiltonian as it can be written at order ε in the form: and where J is the antisymmetric matrix It is not a canonical Hamiltonian system as y v and p v are not conjugate variables.The conserved "Hamiltonian energy" in phase space y v , p v is the frequency Ω v .This conserved quantity is indeed left invariant by the change of phase space coordinates: Ω v ( y v , p v ) = Ω s ( y s , p s ), up to order ε, as noted by Littlejohn and Flynn [9].
The same set of gauge-invariant ray tracing equations can actually be obtained with a different derivation.Starting from a variational formulation of the linearized dynamics with a similar scaling, Perez et al recovered ( 47)- (49), and discussed geophysical applications.Our contribution here is to clarify the reason why the term Ω 1B that contributes to the symbol of Ω does not enter into the ray tracing equations with gauge-independent phase space coordinates.Note that Ω v is denoted by Ω in Perez et al [13].
We conclude that the ray tracing equations describing the trajectories of the center of mass of the multicomponent wave involve only gauge-independent terms.The presence of Berry curvature corrections to those trajectories makes the Hamiltonian dynamics non-canonical.

Application to shallow water waves
For shallow water waves described by the operators Ĥ± , the Berry curvature is computed by using the expression ( 27) for χ ± 0 in (50).This yields an explicit expression of the Berry curvature for each of the three wave bands As noticed in [13], there is no Berry curvature correction for the flat geostrophic wave band, but inertia-gravity wave packet trajectories are influenced by those corrections.We show in the next section that this Berry curvature is related to peculiar topological properties of the symbol of the rotating shallow water wave operator.

Topological properties of Matsuno symbol
The aim of this section is to interpret the Berry curvature (54) as a limiting case of a more general computation performed in [5], that highlights the topological origin of this quantity.
For that purpose, we need to step back to the original Matsuno shallow wave operator Ĥk defined in (3), and to compute its symbol.This is done by noting that each term of the matrix operator in ( 3) is an elementary operator ĝ( ỹ, ∂ ỹ ) whose symbol has already been computed in (A.4): Note that we have introduced the wavenumber l = p/ε in zonal direction, that will be more convenient to manipulate than the rescaled wavenumber p. Recall also that ỹ = |k| y.One can then check that the symbols of the original shallow water wave operator and of the rescaled shallow water operators used previously in the semi-classical computations are related by . Note that there is no higher order corrections to the symbol, due to its simple form that involves no products of derivatives by functions of y.Such higher order corrections would emerge for instance when considering the effect of a varying bottom topography in this model, as detailed in [23,24].
Our strategy now is (i) to compute the eigenvectors of this symbol matrix, (ii) to introduce the Berry curvature describing how the polarization relations of those eigenvectors change when parameters (k, l, ỹ) are varied, (iii) to show that such Berry curvature is generated by a topological singularity in parameter space (iv) to recover the Berry curvature terms (54) as a limiting case.Those steps will be essential to bridge the gap between spectral flow properties of the operator Ĥk , and ray tracing equations deduced from operators Ĥ± in the semi-classical limit k → ±∞.

Matsuno symbol and Kelvin plane waves.
The symbol matrix (55) can be interpreted as the Fourier representation of the wave operator considered in the original work of Kelvin for shallow water waves on the unbounded tangent plane to the sphere, assuming ỹ = f as a constant [31], see Fig.  of those shallow water plane waves considered originally by Kelvin [5].The Kelvin wave problem is simpler to solve than the Matsuno wave problem as it only involves the diagonalization of the 3 × 3 matrix (55).
Existence of a degeneracy point.We now assume that ỹ is a constant parameter.Except in the particular case (k, l, ỹ) = (0, 0, 0) The symbol matrix (55) admits three eigenvalues, namely 0 and ± k 2 + l 2 + ỹ2 .These three eigenvalues define three distinct wave bands when (k, l, ỹ) are varied, except at the particular point (k, l, ỹ) = (0, 0, 0).This point is peculiar as it corresponds to a degeneracy point for the eigenvalues, where the three wave bands touch each other.This band-touching point will play an important role later on.

Berry curvature for eigenvectors of Matsuno symbol
Berry curvature in (k, l, ỹ)-parameter space.An explicit expression for the normalized eigenvectors χ(k, l, ỹ) of the symbol matrix (55) is given in [13].From the knowledge of these eigenvectors, one can compute the Berry curvature F, which is conveniently expressed as a vector in parameter space (k, l, ỹ): where each component is obtained by applying the formula The Berry curvature of the three shallow water wave bands was computed in [5]: with r = (k, l, ỹ) and r its norm.Those vector fields are displayed in Figs.4a and 4b for negative and positive Poincaré wavebands, respectively.
Let us now explain how this Berry curvature vector is related to the Berry curvature introduced previously in (50).
From (k, l, ỹ)-parameter space to phase space ( y, p).The eigenvectors χ ± 0 ( y, p) of H ± 0 are given by the eigenvectors χ(k, l, ỹ) of the matrix H k defined in (55), using the change of variable The Berry curvature introduced in (50) for ray tracing equation is then recovered from the first component F l ỹ = −F ỹ l of the Berry curvature introduced through (57), by expressing the semi-classical limit ε → 0 as a limiting case for k: This expression will play a crucial role in the next section, to relate the spectral properties found by quantization of ray trajectories to the topological properties in parameter space.

From Berry curvature to the first Chern number and Berry monopoles
The first Chern number.The Berry curvature introduced previously is a geometrical quantity, as it describes local properties of eigenvectors in parameter space.In loose terms, it describes how twisted is the eigenvector field locally, independently from any phase choice of the eigenvectors.Let us now introduce a topological invariant, the first Chern number, an integer that counts singularities in families of eigenvectors parameterized over a closed surface.Although more detailed and formal definitions of this invariant exist, we give below a definition of this number as an integral quantity involving the flux of Berry curvature across a closed surface.
Chern-Gauss-Bonnet formula.When considering eigenvectors parameterized over a 2D closed surface embedded in the 3D parameter space, one can compute a Berry flux induced by this curvature across any surface elements da as F • nda, with n a unit vector normal to the surface, as displayed in Fig. 4. When integrated over the whole surface and normalized by 2π, we get an integer, which is the first Chern number: (61) is a generalization of the more familiar Gauss-Bonnet formula that relates the integral of the (geometrical) Gaussian curvature over a closed surface to the (topological) gender of this surface (the number of holes).Here, we consider a simply connected surface that can be continuously deformed to a sphere, and the Chern number appearing in (61) counts the number of phase singularities associated with the bundle of eigenvectors that are parameterized over the surface.This first Chern number is a topological index as it can not be changed under continuous deformation of the surface S, provided that no degeneracy point is crossed during the deformation, just as the number of holes of a given closed 2D surface is not changed by continuously deforming this surface.Since the first Chern number is a topological invariant, singularities can be moved on a given surface (for instance by changing the phase choice for the eigenvectors), but can not be removed.

Application to shallow water waves.
Coming back to the particular case of the shallow water wave problem, it turns out that only two configurations are possible depending on the surface S, when computing (61) using the Berry curvature definition in Eqs (56)-(57): If the surface S encloses a degeneracy point, the Chern number can be nonzero.Otherwise, the Chern number is always zero.For this reason, it is sometimes said that the band degeneracy point carries a topological charge for a given wave band, whose value is given by the Chern number that can be computed either from (61), as done originally in [5], or by other methods [7].Considering the parameter space (k, l, ỹ) the result of the computation7 is a triplet of Chern numbers associated with the three wave bands (by increasing order of frequency): Thus, the positive-frequency Poincaré wave band is described by the first Chern number C = 2.
Our aim in the next section will be to explain how this topological invariant is related to the observed spectral flow of 2 in Matsuno spectrum through ray tracing dynamics in the semiclassical limit.
Physical interpretation: analogy with a magnetic monopole.As noted previously, the Berry curvature does not depend on the phase choice for the normalized eigenvector χ.It describes how fast the polarization relation changes when parameters are varied in the vicinity of a given point in parameter space.A direct consequence of (61) together with the existence of a degeneracy point carrying a non-zero Chern number is that the Berry curvature diverges close to the band-degeneracy points.In fact, one can interpret this curvature F as being generated by a Berry monopole, in the same way as a divergent or convergent magnetic field would be generated by a positive or negative magnetic monopole, whose charge must be quantized, according to a celebrated work by Dirac [6,8].In this analogy the first Chern number plays the role of the magnetic charge of the monopole.The Berry monopole can be interpreted as a hedgehogs topological point defect for the Berry curvature vector field F. This whole section on Berry curvature and Berry monopoles can now be summarized as follows: To conclude, the presence of a non-zero Berry curvature term in ray tracing equations is related to the presence of a Berry monopole at the origin of parameter space (k, l, ỹ), where the three bands touch each others.This Berry monopole generates the Berry curvature away from the degeneracy point, and, in particular, generates the component F y p involved in ( 47)- (48) for ray trajectories in phase space ( y, p).

Bohr-Sommerfeld Quantization
Ray tracing can be used to find the set of eigenfunctions of the wave equation in the limit ε → 0. The standard semi-classical procedure is to look for closed orbits indexed by the frequency ω in phase space, and to select those orbits such that the phase gained by the wave after a period (in phase space) is a multiple of 2π, taking into account additional phase jumps of ±π/2 picked up at turning points [15,32].Turning points occur where the wavenumber vanishes (p = 0).WKB expansion fails at such point, but the solution can be patched with another ansatz, and matching those two solutions in the asymptotic limit ε → 0 yields to a phase jump π/2 (see [32,33]).Those phase jumps are related to a Maslov index [34].To describe eigenmodes within the ray tracing framework, the phase φ originally defined in (34) is now assumed to be separated like φ( y, t) = φ ′ ( y) − ωt.In the following discussion we describe φ ′ ( y) while dropping the prime to simplify notations.
In phase space with canonical coordinates ( y s , p s ), ray trajectories with frequency ω are found by solving This is the equivalent of the Hamilton-Jacobi equation in classical mechanics, with the phase φ playing the role of an action.The phase picked up along a segment d y is We explained in the previous section that computing trajectories in phase space ( y s , p s ) is not easy and awkward, as the computation involves gauge-dependent terms.Following Littlejohn and Flynn [9], it is more convenient to consider phase space variable involving gauge-invariant coordinates where Ω v is defined in (49).Using the relations ( 40)-( 41), the integral transforms into where d g is the differential of the function g = −p v iχ † 0 • ∂ p v χ 0 in phase space.A proof of this can be found in [9] (p.5249).
Let us now consider the case of a closed ray trajectory of frequency ω, with two turning points.This will be the case of all ray trajectories for equatorial beta plane shallow water waves to be considered later.Let us compute the total phase ∆φ picked up along a closed trajectory winding the origin once clockwise.First, we note that the term involving dg vanishes.Second, we add the contributions from the WKB solutions to the phase jump picked up at turning points: where the contour integral are taken along the trajectory given by (65).For more complex trajectories with µ ∈ turning points, the additional term επ should be replaced by εµπ/2.The second term of the r.h.s.involves the Berry connection defined in (42).Using the Stokes theorem, it can be expressed as a flux of Berry curvature across the surface delimited by the closed trajectory at frequency ω: where the integral is performed in the region delimited by the closed trajectory solution of (65).Finally, the quantization condition is obtained by stating that the phase φ/ε picked up along the trajectory must be an integer multiple of 2π to ensure that the WKB ansatz is singlevalued.The quantization condition is finally expressed as Note that some values of m may not be associated with solution of this equation.In fact, we show in the next section that m admits a lower bound for shallow water waves.Note on the literature.This quantization relation has long been used in condensed matter context to address the role of Berry curvature generated by a topological charge in shaping electronic waves, see e.g.[11].Yet there seems to be so far little discussion on how to use those relation to describe the spectral flow.This is the aim of the next subsection, in the case of shallow water wave.

Application to shallow water waves
Let us now consider the two scalar operators Ω± describing the positive-frequency Poincaré wave band of Ĥ± in the semi-classical limit ε → 0. To find the eigenvalues ω ± of the operators Ω± , we apply the quantization procedure introduced in subsection 5.1.Ray trajectories in phase space ( y v , p v ) are found by solving (65), using where Ω 0 is given by (26) and Ω ± 1A is given by ( 28).This computation leads to circular trajectories of radius ϱ(ω), and eigenvalues ω ± are solutions of Admissible values of ω ± are then obtained by applying the quantization relation ( 69), which leads to where the integral in the r.h.s. of ( 71) is just the area inside the (clockwise) trajectory, and where the Berry flux is obtained by injecting in (68) the expression of the Berry curvature F p y = −F y p given in (54): The functions ϱ(ω) and Γ ± (ω) for the positive frequency Poincaré wave band are displayed in Fig. 5a and 5b, respectively.Our aim is now to compare the spectra of Ω+ and Ω− , taking advantage of the quantization condition in (71).
We need to choose a convention to index the eigenfunctions of both operators, in such a way that it may be possible to pair them together based on some physical criterion.For that purpose, it is useful to consider the limit of large frequency ω → +∞.In this limit, the Berry curvature terms Γ ± tend to ±2π, and the phase space trajectories ω = Ω ± v ( y v , p v ) become identical at order ε.Indeed, the symbols Ω ± are expressed as given in ( 28)- (29).A direct inspection of those first-order correction terms for the positive-frequency Poincaré wave band shows that they vanish in the large frequency limit, for which ϱ(ω) → +∞.This implies Because of the convergence of their symbol to a common expression, spectral properties of the operators Ω± are also expected to converge in the large frequency limit.The pairing procedure between eigenmodes of Ω+ and Ω− can thus formally be performed by assigning a common index n to those modes in the large frequency limit.This is done by choosing The term ±1 cancels the Berry curvature terms in (71) in the large frequency limit.While the WKB solution is valid in the limit ε → 0 for a given value of ϱ, the quantization condition (71) is not guaranteed in the limit ϱ → 0, for a given ε, as one can not distinguish anymore phase jumps at the turning points from the phase gained by the WKB solution.More precisely, (71) is valid with the scaling ϱ ≫ ε, that makes possible a scale separation between the WKB part of the solution and the phase picked up at the turning points.Fortunately, in the limit ϱ → 0 the operators Ω± are described by shifted quantum harmonic oscillator operators, which allows us to find an explicit computation for their eigenmodes at finite n ± , as shown in Appendix B. This part of the spectrum can then be matched to the semi-classical computation presented in this section, as the quantum harmonic oscillator soution remains valid for ϱ ∼ ε α with 0 < α < 0.5.An important outcome of the quantum harmonic oscillator computation is the condition: In fact, these results could have directly been deduced from Bohr-Sommerfeld rule (71) by noting that the area term must be positive and that Γ ± = 0 in the limit ϱ → 0, as displayed Fig. 5.The constraints (75)-( 76) mean that two modes labelled by n = −1 and n = 0 are unpaired!This was precisely expected from the bulk-boundary correspondence [5,7].In fact, one can check that • n = −1 is the Kelvin mode, • n = 0 is the Yanai (or mixed Rossby-gravity) mode, and that n ≥ 0 counts the number of zeros for the field v, which is an invariant of a given branch in the dispersion relation of equatorial shallow water waves [1].
To conclude, ray tracing and quantization condition in the limit ε → 0 allow us to recover a semi-classical version of the spectral flow result: just as two modes are gained by the positivefrequency Poincaré wave band as k is varied from −∞ to +∞, the operator Ω+ admits two more eigenmodes than the operator Ω− in the semi-classical limit ε → 0.

Mode imbalance interpreted by Chern-Gauss-Bonnet formula
The analysis of subsection 5.2 shows that the imbalance of 2 modes between operators Ω+ and Ω− in the semi-classical limit is due to the Berry curvature flux corrections Γ ± involved in the Bohr-Sommerfeld quantization conditions, with as illustrated graphically in Fig. 5b.We explained in subsection 5.2 that the Berry curvature term F ± y p involved in the expression of Γ ± in (68) is related to the presence of a Berry monopole in parameter space (k, l, ỹ) for the eigenvectors of the symbol matrix (55).
In fact, (77) can be interpreted as the direct outcome of the Chern-Gauss-Bonnet formula (61).To show this, let us consider the parameter space (k, l, ỹ), and the closed cylindrical surface S depicted Fig. 6.The length of the cylinder is 2k = 2/ ε, while its circular ends are delimited by a closed circular ray trajectory with a radius ϱ(ω) in phase space (±1, y, p), which is related to parameters (k, l, ỹ) through the relation (59), where the index v for phase space variables ( y, p) has been dropped for convenience.The surface S encloses the degeneracy point (0, 0, 0).Thus, according to Chern-Gauss-Bonnet formula (61), the Berry flux across this surface normalized by 2π is equal to the Chern number carried by this degeneracy point for each wave bands, which are given by (62) for the shallow water case.The normalized Berry flux can be decomposed into three contributions: where S c y l represents the open cylindrical surface, and where we have used F ỹ l d f dl = F y p d ydp, consistently with the definition of Berry curvature in (57).
To estimate those three contributions to the Berry curvature flux, we consider a double limit whose order matters: first, the limit of infinite radius ϱ → +∞ for a given value of k second, the semi-classical limit k → ±∞.The motivation for this double limit comes from the fact that we want to interpret how global properties of the full spectrum are changed when k is varied.The first limit allows to get properties of the full spectrum for a given value of k, since this limit allows to scan all possible trajectories in phase space ( y, p).The second limit allows us to get asymptotic properties of this full spectrum in the semi-classical limit.
The first limit ϱ → +∞ implies that the contribution of Berry flux across the open cylinder surface S c y l vanishes, as the Berry curvature vanishes at infinite distance from the origin in parameter space [5].More precisely, the Berry curvature decreases as ϱ −3/2 , which is faster than the increase in the area S c y l ∼ ϱk.Once this limit has been taken, the only non-zero contribution to the Berry flux across S comes from the two circular surfaces at the ends of the cylinder.The Berry flux across those surfaces only involve the component F y p which tends to F ± y p computed in (54) when taking the limit k → ±∞.We note that the limit ρ → +∞ is equivalent to the limit ω → +∞ for positive-frequency Poincaré waves, and that k → ±∞ is equivalent to the semi-classical limit ε → 0. Finally, we get which is equal to a first Chern number C = 2, according to the Chern-Gauss-Bonnet formula (61) and expression (62) for the Chern number in shallow water case.This shows the topological origin of the integer number 2 in the r.h.s. of (77).A last remark is in order: Our point here is to stress that ray tracing followed by quantization in an appropriate semiclassical limit offers a physically appealing intuitive explanation on the relation between the topological index and the spectral properties of the operator.In that respect, it may complement previous lecture notes on this topic [7,8].While the derivation has been focused on the shallow water wave problem, the method illustrated with this example applies to a much broader class of wave problems.
Figure 6: Parameter space and phase space for ray tracing.A Berry monopole is located at the origin of (k, l, ỹ)-parameter space, which corresponds to a degeneracy point for the eigenvalues of matrix H k .Phase space for ray tracing associated with the scalar wave operators Ω± are recovered in the semi-classical limit ε → 0, with the change of variable from (l, ỹ) to (p v , y v ) given in (59).The surface S in parameter space is a cylinder where both circular ends are delimited by closed phase space trajectories of radius ϱ at positions k = ±1/ ε.

Spectral flow and phase space density of states
Liouville theorem.We noticed that ray tracing equations ( 47)- (48) in phase space ( y v , p v ) are an example of non-canonical Hamiltonian systems, which are commonly encountered in fluid dynamics [35].For such systems, the equivalent of standard Liouville theorem holds [35].To see this, let us start with the canonical Hamiltonian system (44)- (45) for trajectories in phase space ( y s , p s ).In that case phase space volume element δV = d y s dp s is trivially conserved owing to the canonical structure of the ray tracing equations.Now, the non-canonical ray tracing equations ( 47)- (48) were derived from this canonical set of equations by performing the change of variable ( 40)- (41).Applying this change of variable to the volume element then yields to the conservation of This means that a Liouville theorem holds in phase space ( y v , p v ), albeit with a modified density of states that depends on the Berry curvature.
Counting the number of states.Owing to the uncertainty principle, a given state of the system occupies at least the elementary phase space volume 2πε.According to standard arguments, this elementary phase space volume just reflects the impossibility to have a perfectly localized wave packet with a given wavenumber [36].One can then count the total number of states that can be hosted in phase space as the integral of δV over phase space divided by the elementary phase space volume 2πε.Here, this number is infinite.However, just as one can compute mode imbalance between operators Ω+ and Ω− , we can compare the number of modes that are hosted in phase space ( y v , p v ) for ray trajectories ruled by Ω+ and Ω− .The corresponding volume elements are denoted by δV ± , and a direct application of (80) leads to: where the integrals are performed over the whole phase space ( y v , p v ).The last equality follows from previous computations involving (54).
We see that the mode imbalance between operators Ω+ and Ω− can be interpreted as a consequence of the continuous deformation of the density of state (1 − εF ± y p )/(2πε) induced by the Berry curvature in phase space ( y v , p v ).In condensed matter context, this phenomenon has been realized and popularized by Niu and coworkers [10].In fact, the modification of phase space density has long been known by physicists working with non canonical hamiltonian systems, e.g.[35].Yet the relation with a spectral flow was not explicit in those previous studies.This relation can now be stated as follows for the shallow wave problem: The Berry curvature induced by a Chern-Berry monopole of charge 2 in ( y, p, k)-space modifies phase space density in ( y c , p c )-space such that 2 more states (among an infinite number) are hosted in the limit k → +∞ than in the limit k → −∞.

Main physical ideas underlying the bulk-interface correspondence
By considering the particular case of equatorial waves, we have proposed a novel physical interpretation of the interface-bulk correspondence that is commonly invoked to interpret the wave spectrum of continuous media with spatially varying coefficients.Let us summarize the main points of the argument to emphasize the physical concepts that have been introduced along the way and to show how they are related to previous work.
In the case of equatorial shallow water waves, Mastuno discovered in 1966 a spectral flow of index 2, corresponding to two modes gained by the upper frequency waveband when the zonal wavenumber k was varied from −∞ to +∞ [1].He found that those modes were more localized than the others along the equator, and that they only propagate energy eastward.In 2017, [5] noticed that this spectral flow index is related to a Chern number describing singularities in families of plane wave solutions for the same shallow water flow model albeit in a simpler f -plane configuration, where the Coriolis parameter is held constant in space, and thus considered as an external parameter.Concretely, [5] computed this Chern number carried by a three-fold degeneracy point for the f -plane wave bands in (k, f , p) parameter space, with p the wavenumber in meridional direction, and interpreted it as the charge of a Berry monopole.The authors concluded that global properties of the dispersion relation in the the beta-plane (complicated) problem could be predicted just by computing the polarization relation of the (simpler) f -plane wave problem, as expected from an interface-bulk correspondence and index theorems [7].
Mastuno's computation was performed for arbitrary wavenumber k with brute force analytical computation.Here, we have recovered Matsuno's result in the particular limit |k| → +∞.The interest of considering this particular limit is to take advantage of the small zonal wavelength parameter to perform a multiple scale development akin to semi-classical or WKB analysis.In this particular situation, the semi-classical computation is less elegant that the direct solution of the problem valid for all k, but is extremely useful to get a physical interpretation of the result, and to be applied to a much wider class of wave problems.
In the semi-classical limit, the computation of the wave spectrum boils down to computing trajectories of wave-packets in phase space ( y, p), and to select trajectories such that the phase picked up by a wave-packet along a closed contour is an integer multiple of 2π.
At lowest order, local properties of the wave-packet trajectory at point ( y, p) in phase space are obtained by assuming a local plane-wave solution, with f ( y) considered constant at the scale of the wave-packet.This local plane wave solution comes along with a polarization relation, and a dispersion relation, which are obtained by solving a matrix problem, referred to as the bulk problem -or the symbol in mathematics.As explained by [7] for a more general class of wave problems exhibiting spectral flows, this procedure makes a direct connection between the parameter space (k, f , p) considered in [5], and the phase space ( y, p) for wavepackets with zonal wavenumber k sufficiently large.Wigner-Weyl transforms and symbolic calculus is just a standard technical way to formalize and perform this translation from the operator to the phase space with local plane wave solutions.
Our contribution here is to explicitly compute the ray trajectories in phase space for a wave-packet in a given wave band, including first order corrections within the semi-classical expansion framework.The key point of our analysis is thus to compute the wave trajectory taking into account a small correction in the polarization that results in asymmetry for phases of left-ward and right-ward propagating waves.The formal development leading to those first order corrections can be found in previous studies following the seminal work of Littlejohn and Flynn [9,10,12,13].
Once trajectories are computed, the second key step is to use the Bohr-Sommerfeld selection rule to find the discrete set of wave-packet trajectories that correspond to eigenmodes of the equatorial beta-plane shallow water wave problem, and to compare the sets of eigenmodes in the limit k → +∞ and k → −∞, for the upper frequency Poincaré wave band.Any mismatch is interpreted as a spectral flow for varying k.Indeed, wave branches in the dispersion relation can neither be created nor annihilated as k is varied [37].The reason is that the wave operator is self-adjoint for any given k and admits only discrete spectrum.Thus it is possible to follow the different dispersion branches by varying k.The set of eigenmodes for any value of k defines a complete orthonormal basis for triplet of 1D fields in y direction.Thus, if a wave-band has more modes for k → +∞ than for k → −∞, it means that the additional modes have transited from other wave-bands as k was varied.
The phase picked up by the wave-packet along a closed trajectory in phase space ( y, p) involves two contributions, which can be interpreted as the dynamical and geometrical (Berry) phase.
The dynamical phase is just the phase that one would expect from inspection of the dispersion relation.This dispersion relation is at lowest order the dispersion of the f -plane problem with the value of f at the averaged location of the wave packet.This dispersion relation is the same for k > 0 and k < 0, and thus involve no asymmetry in the selection rule.
The geometrical phase corresponds to a mismatch of the (complex) polarization phase of the eigenvector after circulating along the closed orbit.This is very much like the angle gained by the pendulum in Foucault experiment after a day, after the system has completed a closed trajectory along a latitude circle.In Foucault pendulum case, the angle increases from 0 to 2π when varying the latitude from the 0 to the North pole.Similarly, for shallow water wave problem in phase space with k > 0, the geometrical phase picked up by the wave packets also increases from 0 to 2π when the radius of ray trajectory increases from 0 to infinity.The opposite phase is picked up for k < 0. The geometrical corrections in phase space ( y, p) are directly related to the Berry monopole described and computed in 2017 by [5] in parameter space (k, f ( y), p).The key message here is that the phase difference of 4π corresponds both to this Berry-Chern monopole and to the mismatch of two modes in the Poincaré wave band for k → +∞ and k → −∞, and ray tracing bridge the gap between those two point of views.
It should be stressed here that gradients of the Coriolis parameter f involve additional corrections to the dynamical phase that we computed, building upon [9,10,13], and that are of the same order as the corrections due to the geometrical phase.It turns out that such corrections do not lead to a net gain or loss of modes when counting them at a given value of k, even if they play a prominent role in setting the value of the frequency levels.Such corrections may break the k ↔ −k symmetry, but we showed that they actually vanish for closed orbits with a sufficiently large radius.
Another subtlety of the approach is that the Bohr-Sommerfeld rule applies to trajectory with sufficiently large radius.We explained that this regime could be matched asymptotically with another regime for which the wave operator admits explicit solutions, corresponding to the eigenmodes of a shifted scalar quantum harmonic oscillator.Again, the precise value of the frequency levels depends on this procedure, but not the number of modes.
The ray tracing argument thus explains how a topological defect in parameter space leads to a mismatch of wavemode numbers for the spectrum in the semi-classical limit.It also predicts that the modes involved in the spectral flow have a lower index than the others, meaning that the modes are more localized than the others at a given k, since the trajectories in phase space associated with those modes are closer to the origin.It would be interesting to show in this ray tracing framework that a deformation of the beta plane into a smooth step profile for f ( y) would indeed lead to the delocalization of all the modes but those two additional modes of topological origin, as described for instance in [38] by computing explicitly the spectrum in a solvable case.
We also showed along the way that the noncanonical structure of the Hamiltonian ray equations allows for a second, physically appealing interpretation of the spectral flow: in the limit of large wavenumber |k|, the density of state in phase space ( y, p) is modified in such a way that two more states can be accommodated for k > 0 than for k < 0. Again, this density of state is governed by the presence of a Berry monopole in parameter space.
In that respect, our study provides a physical interpretation of the index theorem, complementary to more rigorous approaches [7].In fact, our study relates a spectral flow index to a topological index, while Atiyah-Singer theorem involves an analytical index, which, in the shallow water case and other relevant physical systems, has been suitably defined and interpreted in a recent study by Delplace [8].In the present paper, a semi-classical approach was used to reconstruct formally spectral properties of a wave problem with a parameter varying in a given direction.On a more rigorous side, semi-classical technique has been recently used in mathematical context to describe how a given wave packet evolving between two topologically distinct materials splits into a "bulk" part that spreads and eventually collapses and "edge" part that remains coherent and propagates at the interface between the two topological phases [39].

Perspectives
We stress that we focused on a semi-classical interpretation of a bulk-interface correspondence, in problems without boundaries.In such problems, waves of topological origin emerge along an interface of a parameter that opens a frequency gap for the bulk waves, as the Coriolis parameter for equatorial waves.While the method was presented in the case of equatorial shallow water waves, it is sufficiently general to be applied to a much wider class of problems including more general shallow water waves [24,52], plasma [27,28], continuously stratified and compressible flows [25,26], among others.In all those problems with Hermitian wave operator, topological features are found by looking for Berry monopole in parameter space, and we have shown how such monopole may affect wave-packet dynamics in phase space.The bulk-interface correspondence discussed in this paper is different than the bulk-boundary correspondence that is commonly encountered in condensed matter.In the case of a bulkboundary correspondence, a bulk Chern number can be defined on each side of the interface (e.g. each side of the equator for shallow water waves).Assigning a topological invariant on each side of the interface for continuous media is not always possible.For instance, for rotating shallow water waves, it requires the introduction of a regularization parameter such as odd-viscosity [38,[41][42][43].When this parameter is added into the model, one can assign a topological index to the f -plane problem, i.e. in each hemisphere, and predict the number of unidirectional modes at the equator [38].The case of sharp (discontinuous) interfaces or hard boundaries such as solid walls is much more complicated for continuous media: unidirectional modes exist [43], yet apparent violation of standard bulk-boundary correspondence have been noted and explained [44,45].Our present work applies to the unbounded case and as such provide an explanation for the bulk-interface correspondence only.However, hardboundary cases may sometimes be interpreted as limiting cases of an interfaces.For instance, unidirectional trapped modes along coasts found by Kelvin in 1880 [31] can be recovered from the study of a shallow water with varying bottom topography defining an interface at the coast [24].In that context, the ray tracing approach may help to bridge the gap between the modern topological point of view and the more traditional textbook skipping orbit picture explaining the emergence of the unidirectional mode along the wall.In that context of coastal waves, it will be interesting to relate the ray tracing approach to the interpretation of rotating shallow water dynamics as a Chern-Simons topological field theory [46].
It will also be interesting to see in future work how non-Hermitian wave problems can be tackled with the ray tracing approach.Such problems naturally occur in the context of flow instabilities [47], or in the presence of dissipative terms in the dynamics.On the one hand, some of the spectral flow properties found in the Hermitian case seem to be robust to the presence of non-Hermitian effects, as reported in the context of shallow water dynamics in the presence of a mean velocity shear [52].On the other hand, intriguing new unidirectional trapped modes have been reported in rotating convection [48,53,54].The ray tracing point of view may help to bring new hindsight on the topological origin of those phenomena.where all functions inside the integral are evaluated at y.The term in the exponential is expanded up to the second order because of the 1/ε prefactor.Powers of y ′′ in the integrand can be replaced by derivatives with respect to p in front of the exponential term.Keeping only terms up to order ε yields to:

B Quantum harmonic oscillator limit
We show in this subsection that eigenmodes with sufficiently small n index can be approximated by solutions of a differential equation analogous to the quantum harmonic oscillator problem [36], in the semi-classical limit ε → 0. This amounts to consider the limit of small radius ϱ for trajectories in phase space, or equivalently, the limit of frequencies asymptotically close to one: |ω − 1| ≪ 1.Before presenting the computation, let us stress that operators Ω± are not differential operators in general.By contrast, the original multicomponent wave operator Ĥ± is a differential operator.Differential operators involve only terms such as ∂ n /∂ y n with n a non-negative integer.Consequently, their symbols involve only polynomial terms in the wavenumber p.However, the scalar operators Ω± are obtained after a diagonalization procedure of the symbols H ± , and this diagonalization leads to the presence of non-polynomial terms in the ex- pression of Ω ± ( y, p), see for instance the r.h.s. of (70).Owing to the presence of those nonpolynomial terms, the corresponding operators Ω± belongs to the family of pseudo-differential operators.
In the limit of trajectories that are close to the origin in phase space ( y, p), a Taylor expansion allows us to turn the pseudo-differential operators into differential operators analogous to the quantum harmonic oscillator operator, up to a constant.In the shallow water case, taking the limit ϱ ≪ 1 in (70) leads to The term εϱ 2 can also be dropped in the small ε limit.Then, using Ω ± = ω ± + O(ε 2 ), ϱ 2 = y 2 + p 2 , and applying Weyl quantization to this symbol leads to The eigenvalues of this operator are ω ± = (m ± + 1/2)ε + 1 ∓ ε/2 with m ± ∈ [36].Now, we want to match this dispersion relation with the semi-classical result obtained in the main text under the condition ϱ ≫ ε, while keeping ϱ ≪ 1 to stay in the harmonic oscillator limit.This is done by choosing m ± = n ± 1: The expression n(ω) is drawn as a dashed line in Fig. 5d.

Figure 1 :
Figure 1: a) The Coriolis parameter f as a function of latitude on a rotating planet, and the shallow water model with H the fluid layer thickness, much smaller than the typical horizontal scales of motion.b) Beta-plane approximation: the flow takes place in a plane (x, y) tangent to the equator, with f = βy.

Figure 2 :
Figure 2: Eigenvalues of the beta plane shallow water wave operator defined in (3), as a function of the zonal wavenumber k, for β = 1, c = 1, often referred to as the Matsuno spectrum [1].Two modes (in red) are gained by the positive-frequency Poincaré wave band as k increases.The aim of this paper is to interpret this spectral flow using semi-classical analysis.

3 .
Delplace et al showed that the spectral flow associated with Matsuno wave operator for shallow water waves on the beta plane is encoded in the topological properties of eigenvectors

Figure 3 :
Figure 3: Dispersion relation of shallow water waves on a plane tangent to the planet, assuming constant Coriolis parameter f = ỹ (Kelvin wave problem).Geostrophic modes correspond to the Rossby modes of the Matsuno wave problem.The pink dot corresponds to the band degeneracy point that plays a central role in the topology of shallow water waves [5].

Figure 4 :
Figure 4: a) Berry curvature F vector field of the Matsuno symbol (55) for the negative Poincaré band.b) The same for the positive Poincaré band.Both diverge in amplitude at (k, l, ỹ) = (0, 0, 0).c) An integral of the flux of this vector field on any closed surface S enclosing this point evaluates to −4π and 4π, for the negative and positive Poincaré bands respectively.The pink dot represents a Berry-Chern monopole of charge ±2 that can be interpreted as the source of observed Berry curvature.

Figure 5 :
Figure 5: (a) Variations of the circular trajectory radius ϱ as a function of frequency ω in phase space ( y v , p v ).(b) Variation of the Berry flux across the area delimited by a circular ray trajectory having a radius ϱ.(c) A circular trajectory in phase space.(d) Quantization condition for shallow water waves in the semiclassical limit, as expressed in (71)-(74).Red and blue curves are associated respectively with operators Ω+ and Ω− .The sign ± is due to the definition n ± = 1 Dashed lines correspond to explicit computations performed in the limit ε → 0 at finite n, see Appendix B.