Modelling Neutron-Star Ocean Dynamics

We re-visit the calculation of mode oscillations in the ocean of a rotating neutron star, which may be excited during thermonuclear X-ray bursts. Our present theoretical understanding of ocean modes relies heavily on the traditional approximation, commonly employed in geophysics. The approximation elegantly decouples the radial and angular sectors of the perturbation problem by neglecting the vertical contribution from the Coriolis force. However, as the implicit assumptions underlying it are not as well understood as they ought to be, we examine the traditional approximation and discuss the associated mode solutions. The results demonstrate that, while the approximation may be appropriate in certain contexts, it may not be accurate for rapidly rotating neutron stars. In addition, using the shallow-water approximation, we show analytically how the solutions that resemble r-modes change their nature in neutron-star oceans to behave like gravity waves. We also outline a simple prescription for lifting Newtonian results in a shallow ocean to general relativity, making the result more realistic.

In order to provide appropriate context, let us first consider the relevant observations.Thermonuclear (Type I) X-ray bursts are sudden eruptions due to the unstable burning of hydrogen and helium on the surface of an accreting neutron star in a low-mass X-ray binary [6][7][8].The fast timing capability of RXTE led to the discovery of narrow high-frequency features (mostly in the range 300 − 600 Hz) in the power spectrum of these X-ray bursts (see, e.g., Refs.[9,10]).These features became known as burst oscillations, either found during the rise/decay phase of the burst, or during both phases.The remarkable stability of the asymptotic oscillation frequency indicates that burst oscillations are linked to the neutron-star spin frequency (see, e.g., Ref. [11]).This notion was further supported by the detection of the oscillation frequency in the decay phase of thermonuclear bursts from the accreting millisecond X-ray pulsar SAX J1808.4−3658, which matched well with the known spin frequency (401 Hz) of this source [12].Currently, only a handful (about 20 − 21) of neutron-star low-mass X-ray binaries are known to exhibit burst oscillations [13].These systems are also known as nuclear-powered millisecond X-ray pulsars.The observed burst oscillations frequencies are, however, not always perfectly stable.Observations suggest a typical drift of 1 − 3 Hz.Conversely, it has been observed that the oscillation frequency does not significantly drift during the burst decay of persistent accreting millisecond X-ray pulsars and nuclear-powered millisecond X-ray pulsars, and typically remains close to the pulsar frequency in these cases.As an example, for XTE J1814−338, the burst oscillation frequency is extremely stable-within 10 −8 Hz of the pulsar spin frequency [14]-while for SAX J1808.4−3658,IGR J17511−3057 and IGR J17498−2921, the former frequency is within a few mHz, ∼ 0.05 Hz and ∼ 0.25 Hz of the latter, respectively [12,15,16].However, AstroSat observations of the neutron-star X-ray transient XTE J1739−385 during its 2019-2020 outburst revealed accretionpowered pulsations at 386 Hz during very short intervals (0.5 − 1 s) of X-ray flares, making it an intermittent accreting millisecond X-ray pulsar [17].One of the thermonuclear bursts observed during these observations showed the presence of coherent decay phase burst oscillations at 383 Hz, indicating a frequency offset of ∼ 3 Hz away from the spin frequency of the source.In addition, there have been observations of frequency drifts of ∼ 3.6 Hz for 4U 1916−053 [18] and ∼ 5 Hz for MXB 1659−298 [19].These frequency drifts-taking place in a few seconds-are significantly larger than the typical ≤ 1% feature.However, in both cases, the burst oscillation signal drops below the detection threshold in the middle of the burst.Since frequency cannot be tracked continuously throughout the burst, it is not clear whether the observed drift is real [20].One plausible explanation for the observed frequency drift during burst oscillations is the excitation of oscillation modes, the property of which evolves as the ocean settles down after the burst.Within this explanation, different families of modes (associated with different restoring forces) may lead to a range of observed frequencies.In the following, we aim to shed light on this notion.

B. Theory
The effort to build theory models that match the observed burst oscillation phenomenology draws heavily on the broader effort to understand neutron-star oscillations, with a significant body of work motivated by issues from, in particular, gravitational-wave astronomy [21].The equation of state of the low-density neutron-star ocean is naturally different from that of the star's high-density core, but the mathematical machinery required for the problem of ocean waves is pretty much the same as that used for global oscillations.Having said that, the shallow-ocean problems lends itself to further assumptions.In particular, there is a natural connection with geophysics and work aimed at understanding the Earth's atmosphere and oceans [22].This, in turn, has led to models for neutron-star oceans commonly drawing on the so-called traditional approximation.The first study in this direction was conducted by Bildsten et al. [23], who used the approximation to compute low-multipole g-modes in the oceans of fairly slowly rotating accreting neutron stars.
Following this, the connection between burst oscillations and ocean modes was first made by Heyl [24].He arguedbased on the results from Longuet-Higgins [25], which further assume the shallow-water approximation-that the observations were best explained by an |m| = 1 "buoyant" r -mode, where m is an integer that characterises the mode's azimuthal dependence.This explanation seemed natural because the r -mode solutions occupy a wide band near the equator, resulting in a large variability (required for the oscillations to impact on observations), while these modes have low inertial-frame frequencies (as required by the proximity of the burst oscillation frequency to the star's spin frequency).
Further work on the problem showed that the evolution of the ocean (effectively, cooling after the thermonuclear surface explosion) led to the r -mode frequency drifting towards the spin frequency, thus adding further support for the model [26][27][28].However, Newtonian gravity models suggested that the frequency drift ought to be larger than observed (by perhaps a factor of 2).This issue was addressed by the work of Chambers and Watts [5], who extended the traditional approximation calculation to general relativity (drawing on the results from Ref. [29]).The results from Chambers and Watts [5] indicate mode frequencies of ∼ 2 Hz and frequency drifts of ∼ 2 Hz.This is more in line with the observations, indicating that relativistic gravity is an important part of any quantitative model (as one would have expected).
At face value, the r -mode explanation for the observed quasi-periodic oscillations seems reasonable and the available results are promising.Yet, there is a slight disconnect between the argument and the observations.In particular, it seems difficult to explain systems where the burst oscillation frequency is very close to the spin frequency (as in the case of XTE J1814−33, for example).The problem is that, any mode that has a finite frequency in the rotating frame must lead to an offset from the spin frequency in the (observed) inertial frame.This is easy to see from the relationship between the frequencies: where ω and σ are the mode frequencies as measured by an observer in the frame of the star (rotating with angular spin Ω) and in an inertial frame outside the star, respectively.Clearly, a mode with m = ±1 and ω Ω can lead to |σ| being close to Ω, as in the results of Chambers and Watts [5], but (and this is a significant but) we cannot have σ → Ω.The only way this can happen is if the ocean dynamics relaxes to geostrophic motion [30,31], and this is not represented by the r -mode solutions.
As things stand, it seems evident that we do not yet have final answers to the questions raised by the observations.The mode explanation is appealing and may well turn out to be correct, but the available results do not allow us to explain all observed systems (some: yes; all: no).The natural next step then is to try to understand to what extent the existing calculations are precise enough to warrant matching against the observations.Given that these models rely (universally) on the traditional approximation it is natural to point the bright light of the inquisitor in this direction.This is our aim with this paper.We want to understand to what extent the traditional approximation is robust for the relatively fast-spinning neutron stars we are interested in.This turns out to be an interesting question in its own right and the answer may (in the extension) help us make progress on the questions that motivated the discussion in the first place.We will, however, not be providing the final answers here.

II. OCEAN OSCILLATIONS
In the first instance, let us model the neutron star as a Newtonian, uniformly rotating, perfect fluid.The rotation is assumed to be slow such that the angular frequency Ω Ω 0 ≡ GM/R 3 is small, where M is the star's mass and R is its corresponding (non-rotating) radius. 1 In this limit, the star is well approximated as a sphere.Moreover, our focus will be on the dynamics of the ocean, representing a thin outer layer, which will be small in mass and contribute little to the star's gravitational potential.Therefore, it is appropriate to adopt the Cowling approximation, neglecting variations in the star's gravitational potential.
We are interested in mode solutions of the form e i(ωt+mφ) (i.e., working in the frequency domain and using a Fouriermode decomposition for the φ-dependence).Hence, in a spherical coordinate basis, the linearised Euler equation in the rotating frame has components −ω 2 rξ θ − 2iωΩr cos θ sin −ω 2 r sin θξ φ + 2iωΩ(r cos where ξ j represents the components of the Lagrangian displacement vector, ρ is the equilibrium mass density, g is the local gravitational acceleration and δρ and δp are the Eulerian perturbations of the mass density and pressure, respectively.The mode frequency ω is as measured by a rotating observer.The perturbations must also satisfy the continuity equation To close the system of equations, we require an equation of state.For an adiabatic ocean with frozen composition, it is appropriate to introduce where p is the equilibrium pressure, Γ 1 is the adiabatic index and N is the local Brunt-Väisälä frequency given by Although we have made a number of assumptions up to this point, we find ourselves in a situation that plagues much of hydrodynamics.We need to solve a set of coupled partial differential equations ( 2)-( 4).The problem is no less complicated by the presence of rotation.In order to proceed and determine the mode solutions, we need to separate the radial and angular dependencies.
A. The traditional approximation The historical approach-with origins in geophysics-is to use the traditional approximation.The aim of the traditional approximation is simple: We want to effect the separation of the radial dependence from the angular behaviour.This is achieved by disregarding the vertical component of the Coriolis force Ω sin θ in the Euler equation (2).This involves two assumptions.The first is to discard the term with the radial displacement in Eq. (2c).This is valid when and should apply for modes that are predominantly horizontal.This is appropriate for some low-frequency oscillations (see Ref. [21] for an overview of neutron-star oscillations), like g-and r -modes, but in principle excludes p-modes and general inertial modes. 2 Then, the φ-component of the Euler equation becomes The second aspect of the approximation removes the angular displacement from the r-component of the Euler equation.The motivation for this step is less clear, but one may argue for it based on the expected ocean stratification [23,33].The equation of state (4) combines with Eq. (2a) to remove δρ/ρ, Based on the first assumption ( 6), we have the scalings |ξ r | ∼ h and |rξ θ | ∼ |r sin θξ φ | ∼ R, where h R is the depth of the ocean.As indicative values for a neutron-star ocean, we will take h ∼ 10 − 100 m [34].Now, to neglect the Coriolis term in favour of the buoyancy force in Eq. ( 8), we clearly require Due to the relative smallness of h and the assumed low frequencies of the modes we are interested in (recall that the rotating-frame frequency must be small in order to explain the observed phenomenon), this necessitates a combination of a highly stratified layer with an additional restriction on the rotation rate.For example, for accreting neutron stars, Ref. [23] argue that this limits the spin to Ω/(2π) < ∼ 200 Hz.However, to retain the first term in Eq. ( 8), we must also have Since we have been focusing on modes with low frequencies, this constraint requires the star to be very slowly rotating indeed, certainly beyond what is relevant for the sources we are interested in.In essence, the assumptions associated with the traditional approximation are unlikely to be satisfied for the objects we want to model.Nevertheless, assuming that the two inequalities are satisfied, we arrive at Equations ( 7) and ( 11) are the φ-and r-component of the perturbed Euler equation, respectively, in the traditional approximation.The θ-component (2b) is unaltered.
Before we proceed, it is worth commenting on the character of the wave-solutions expected within the traditional approximation.Focusing on short wavelengths, a plane-wave analysis sheds light on this issue (discussed in detail in Appendix A).The key take-home message is that-contrary to the impression one might get from much of the literature-the traditional approximation is not (strictly) a low-frequency approximation.The equations also support high-frequency sound waves.However, if we insist on consistency and impose Eq. ( 9) [and hence omit the frequency dependence from the radial Euler equation (11)], then the approximation filters out the sound waves and permits only low-frequency gravity modes (supported by buoyancy).It is worth keeping these observations in mind when making use of the traditional approximation.
In the classic geophysics literature, further approximations are made.The outer layer is modelled as a shallow, uniform-density, incompressible fluid (water, obviously), thereby circumventing (as we discuss later in Sec.III A) radial derivatives and reducing the problem to one dimension [25,[35][36][37][38][39] .It turns out that, in this context, the traditional approximation becomes necessary in order to ensure that the perturbations conserve angular momentum [40][41][42].However, this seems to be a pathology of neglecting the radial behaviour, as the radial component of the Euler equation is ignored entirely.
The conceptual issues that arise with the traditional approximation are important, since much of our theoretical understanding of oscillations in neutron-star oceans relies on it [5,23,24,27,28].Based on our arguments, it would seem-as also suggested by Ref. [22]-that the traditional approximation is more of an assumption than an approximation.This suggests that care should be taken when applying it to realistic scenarios.
With this caveat in mind, we note that traditional approximation neatly decouples the radial and angular behaviour in the following way.The perturbation equations now show that ξ r , δρ and δp share the same polar dependence such that3 where Θ is known as a Hough function and µ = cos θ is the latitudinal coordinate.That this separation of variables is possible is a direct consequence of the approximation and is mathematically quite appealing.Given Θ, one is in a position to describe the entire angular behaviour of the perturbations.To see this, the angular components of the Euler equation ( 2b) and ( 7) can be combined to give where To determine Θ, we consider the continuity equation ( 3), supplemented by Eq. ( 4), to obtain where we have introduced the eigenvalue equation with separation constant λ, operator L q defined by and spin parameter q = 2Ω/ω.Equation ( 16) is known as Laplace's tidal equation and fully characterises the angular sector of the perturbation problem.It possesses some noteworthy features.In the non-rotating limit q = 0, Laplace's tidal equation reduces to the associated Legendre equation, where the eigenfunction becomes an associated Legendre polynomial Θ = P m l and the eigenvalue is λ = l(l + 1).This reduction is consistent with the familiar fact that a mode of a spherical star may be uniquely described by a single associated Legendre polynomial with a given degree l and order m.
All the rotational features are captured by Laplace's tidal equation.However, it is interesting to note that the equation knows nothing about the stratification of the star.This information sits in the radial sector, given by Eqs. ( 11) and (15), which is identical to that of a non-rotating star when λ = l(l + 1) and describes an eigenvalue problem for the mode frequency ω, accompanied by the relevant boundary conditions-zero radial displacement at the ocean base and a vanishing Lagrangian variation of the pressure at the surface.This similarity does allow for a correspondence to the case of spherical symmetry, where useful results can be drawn (see, e.g., Refs.[23,33]).

III. SOLUTIONS OF LAPLACE'S TIDAL EQUATION
Due to the mathematical usefulness of the traditional approximation, quite some effort has been dedicated to solving Laplace's tidal equation.Notable mentions include Longuet-Higgins's series of seminal papers [25,37,38], which discusses the topic extensively.Fortunately, with modern computational techniques and machinery, it is rather straightforward to obtain solutions of Laplace's tidal equation through direct numerical integration [23,33].We will discuss the solutions of Laplace's tidal equation here, paying particular attention to features that do (and, conversely, do not) capture the expected physics.
Provided a choice of (m, q), one obtains from Eq. ( 16) an infinite set of eigenvalues λ.As a representative example, we show a selection of eigenvalues for m = −2 in Fig. 1. (Due to the symmetry of Laplace's tidal equation with respect to mq, the sign of q can be reversed to give solutions with m = 2.) These agree with the results of Refs.[25,33].The eigenfunctions Θ are even or odd, depending on their symmetry about µ = 0 (for more information, see Ref. [33]).  1.The eigenvalues λ of Laplace's tidal equation for m = −2 as functions of the spin parameter q = 2Ω/ω.Solutions with q > 0 are prograde, whereas q < 0 correspond to retrograde perturbations.The q = 0 axis corresponds to a non-rotating star.The integer k labels sequences of solutions, where solid and dotted lines indicate even and odd solutions, respectively.The modes from these solutions arise as a consequence of the traditional approximation.The solutions at q = 0 correspond to the oscillation modes of a spherical star.With finite q, these modes become rotationally modified.The other class of solutions only exists for |q| > 1 and can have negative λ.For λ ≥ 0, this class of solutions capture features of the r -modes.
As can be seen in Fig. 1, we can categorise the solutions of Laplace's tidal equation into two classes: (i) perturbations that can exist on the non-rotating star (with label k ≥ 0) and (ii) those that require rotation (k < 0).Since Eq. ( 16) limits to the associated Legendre equation when q = 0, we note that class (i) limit to λ = l(l + 1) and can be identified as the familiar polar modes of a spherical star. 4In general, this includes f -, p-and (in the presence of stratification) g-modes.Although in principle, we would anticipate the high-frequency solutions to be excluded via Eq.( 6), the perturbation equations in the traditional approximation do permit such oscillations (as we have shown in Appendix A).Class (i) are simply the rotationally modified polar modes.
We follow the convention introduced in Ref. [33] and label the sequences with the integer k.The integers k ≥ 0 of class (i) have been chosen so that, in the limiting case q = 0, the solutions correspond to associated Legendre polynomials of degree l = |m| + k.Therefore, for m = −2, the lowest eigenvalue corresponds to l = 2.
Solutions with different signs of q travel in different directions around the star.Modes with mq < 0 are progradethey travel in the same direction as the rotation of the star.On the other hand, modes with mq > 0 are retrograde.As we would expect, there is a splitting of the modes when rotation is introduced, with respect to the q = 0 axis of Fig. 1.This splitting is related to the deviation from polar symmetry that arises through the Coriolis force.For a given |q|, each mode has a partner that moves in the opposite sense around the star.As is shown in Fig. 1 for k ≥ 0, the retrograde solutions have larger λ than their prograde partners.
We note that the k = 0 solution seems to be special among the class (i) modes, in that it does not reach a minimum λ when q = 0.For q > 0, it asymptotes to a fixed value relatively rapidly.Longuet-Higgins [25] described this solution (when q > 0) as a prograde Kelvin wave.This seems to be a feature of the assumptions associated with the traditional approximation.
Class (ii) solutions only arise when |q| > 1.Since these solutions require rotation to exist, they belong to the inertial-mode family.One unique aspect of these solutions is that they may acquire negative values of λ.There has been some discussion of the meaning of such negative eigenvalues.Longuet-Higgins [25] proposed that these solutions may be relevant when the system experiences an external force at some fixed frequency.It has also been suggested that they correspond to convective modes that become stabilised by rotation [33].Whatever the explanation, the λ < 0 solutions are of no particular relevance for our discussion.They are simply included in the figure for completeness.
The λ ≥ 0 solutions of class (ii) are all retrograde.This is a well-known feature of the r -modes [43][44][45].[Recall that we expect general inertial modes to be excluded through the assumption (6).]Furthermore, each of these solutions has the property that when λ = 0, mq = (|m| + |k + 1|)(|m| + |k + 1| + 1).This is precisely the r -mode frequency in the slowly rotating limit, where l = |m| + |k + 1|.We show this analytically below.Therefore, the region near to λ = 0 must be closely related to the slow-rotation approximation.
Finally, there exists a special solution k = −1 that has r -mode character for small |q|, but joins the class (i) family of solutions when |q| becomes large.This feature was observed by Longuet-Higgins [25].In the slow-rotation limit, this solution corresponds to the l = |m| = 2 r -mode (with overtones, if the star is stratified).That this r -mode joins the rotationally modified polar modes at large |q| is surprising and could also be a feature of the particular assumptions in the traditional approximation.
It is also worth noting that we have, from the outset, assumed the star to be slowly rotating.In theory, this will impact the range of q that Laplace's tidal equation is valid for.By the constraint (10), we arrive at For any shallow layer, where h R, this severely limits the values that q can take.This implies that solutions of class (ii) are unlikely to accurately describe the inertial modes they approximate.This further suggests that care should be taken with interpreting results in the traditional approximation.

A. The shallow-water problem
As we have alluded to above, the first attempts of calculating oscillations in the oceans of rotating bodies neglected the radial dependence entirely, approximating the fluid layer as shallow.In this context, the modes are completely described by Laplace's tidal equation (16).Let us show how this comes out.
In the absence of radial structure, r = R and the r-component of the Euler equation ( 11) is ignored. 5Further, the shallow ocean is treated as a uniform-density, incompressible fluid with ρ = const and δρ = 0. Considering the Lagrangian pressure perturbation at the surface, we then find which is assumed to hold throughout the ocean with g = const.This expression may be used to remove δp/ρ from the system of equations.It is worth noting that this violates Eq. ( 11), which has been discarded, as well as the boundary condition on ξ r at the base of the ocean.The shallowness of the layer implies the approximation With these assumptions, Eqs.(2b), ( 3) and ( 7) may be expressed in the forms that are found in the classic literature (see, e.g., Refs.[25,36]).Equation ( 15) readily simplifies to show that where is the Lamb parameter and λ is the separation constant we introduced earlier with Laplace's tidal equation (16).Effectively, Eq. ( 21) shows how the two eigenvalues λ and ω become linked in the shallow-water approximation and the problem reduces to a single eigenvalue.However, it is important to keep in mind that the relation (21) only holds in the shallow-water approximation.
In the historical shallow-water literature, is the favoured eigenvalue.This is accomplished simply by absorbing a factor of q 2 into the definition of Eq. ( 16) and solving for instead of λ. Figure 2 shows the results (for λ > 0) of Fig. 1 expressed in terms of the Lamb parameter.These are in complete agreement with Longuet-Higgins [25] (see his Fig. 3).The left panel shows the prograde perturbations, which only include solutions of class (i).The right panel shows the retrograde solutions, which include members of class (i) and (ii).The solutions that asymptote to constant values at large 1/ 1/2 approach the expected r -mode frequencies.
FIG. 2. The inverse spin parameter 1/q = ω/(2Ω) against 1/ 1/2 = (gh) 1/2 /(2ΩR) for m = −2.The Lamb parameter is related to the eigenvalue λ by = q 2 λ when the fluid is assumed to be uniform density, incompressible and shallow.We only show > 0 solutions.The left panel shows prograde solutions and the right panel shows retrograde solutions.The integer k labels the sequences of solutions, where solid and dotted lines indicate even and odd solutions, respectively.
In the shallow-water approximation, where all the mode features are determined from Laplace's tidal equation, progress can be made analytically, which gives insight into the solutions.(Indeed, this is almost certainly one of the reasons the two approximations, shallow-water and traditional, were adopted in the first place!)With the relation ( 21), Laplace's tidal equation ( 16) can be brought into the form [25] where L 0 is the operator L q , defined in Eq. ( 17), when q = 0 (which is simply the horizontal Laplacian) and y = 1 − µ 2 Θ.Laplace's tidal equation in the form of Eq. ( 23) is convenient to consider limiting cases.
For small values of /|q|, Eq. ( 23) reduces to We further assume small and an immediate solution is found with associated Legendre polynomials for mq = l(l +1), so This is the r -mode frequency in the slow-rotation limit.In this regime, the angular dependence of the displacement vector becomes axial [see Eqs.(14)].These asymptotic solutions can be seen in the right panel of Fig. 2.However, this regime will not be relevant for fast-rotating neutron stars.This is evident from the Lamb parameter, Clearly, the sources we are interested in lie firmly in the large-regime.Suppose we now examine large /|q|, then Eq. ( 23) becomes [38] Assuming that we only need to retain the highest derivative to balance the large terms, close to the equator we have where A = mq + /q 2 and we have changed variables to x = 1/4 µ.This approximation connects directly to Longuet-Higgins's discussion of the asymptotics [25,37,38].In short, the solutions of Eq. ( 28) are where and H ν is a Hermite polynomial of integer order ν.Therefore, by Eq. ( 30), the mode frequencies are now determined by the cubic equation For large , one root is clearly 1/q ≈ m/[(2ν + 1) 1/2 ], leading to Notably, this frequency does not depend on the spin of the star.As discussed by Longuet-Higgins [25] and shown in Figs. 1 and 2, all but one (k = −1) of the solutions of class (ii) tend to this frequency.
In terms of the nature of the modes, it follows immediately from Eq. ( 28) that these waves are trapped in the region close to the equator.A standard WKB argument shows that the solution is exponentially small beyond the turning points and it is easy to see that the modes will be equatorially trapped through a Taylor expansion, which leads to Heyl [24] used these asymptotic results from Longuet-Higgins [25] to argue, based on the structure of the bursts, that the |m| = 1 manifestly retrograde solutions of class (ii) (excluding the k = −1 solution that changes character) best explained the observed features.These he termed "buoyant" r -modes.The terminology seems somewhat confused, however, because the frequency (32) has the scaling one would expect for a surface gravity wave (and bears no resemblance at all to the r -modes, which are linear in the rotation rate).We get a clear hint of what is going on from the right panel in Fig. 2. The modes change nature-essentially undergoing an avoided crossing-as the parameters change from the slow-rotation regime (where the solution behaves like the classic r -mode) to the fast-rotation regime [where Eq. ( 32) applies].
Additional hints that the r -modes must change nature follows by repeating the constant-density calculation from Ref. [47] with a zero boundary condition imposed at depth h.A straightforward calculation leads to the frequency correction (see the original paper for definitions) This clearly diverges as h R, thus violating the slow-rotation ordering assumed in the r -mode calculation (an issue paid close attention to in Ref. [48]).In particular, the l = |m| r -mode assumptions require so we need to have At the surface of the star this leads to so the Lamb parameter cannot be very large.In essence, the l = |m| r -modes cannot exist in a shallow neutronstar ocean.This divergence in the frequency correction ω2 was also observed in numerical calculations we performed following Ref.[48], which incorporate rotational shape corrections.This is in accordance with the asymptotic behaviour outlined by Longuet-Higgins [25,37,38].
Returning to the mode frequency (32), let us ask what kind of values we expect.We have so, in the inertial frame we get [using Eq. ( 1)] In order to be close to the spin frequency we need |m| = 1 and the result then predicts an offset between the spin frequency and the burst oscillation.For the fastest spinning stars, we have Ω ≈ 0.3Ω 0 , so it follows that (for m = 1) If we expect the offset to be 1 %-in order to connect with typical observations of the bursts (see, e.g., Ref. [12])-then we need or so we need ν = 4 − 5 for h ≈ 10 m.This does not seem unreasonable.
It is important to note that we have not considered (since it is not the focus of this work) the impact of the ocean cooling after the burst.As demonstrated by Chambers et al. [28], the effect of nuclear reactions in the ocean, accompanied by varying composition and heat fluxes from the crust, is required to explain the drift of the oscillation modes.In particular, Chambers and Watts [5] brought this discussion into the arena of general relativity, showing a further reduction in the frequency drifts.Clearly, a realistic calculation of oscillations in neutron-star oceans will need to incorporate all of these features.

IV. RELATIVISTIC OCEAN MODES
Even though we know that accurate theoretical predictions for neutron stars require general relativity to be taken into account, the discussion we have provided so far has been entirely in a Newtonian setting.Given this, it is natural to consider how relativity would enter and modify the results we have discussed.A first answer to this question was provided by Chambers and Watts [5], who built on the relativistic version of the traditional approximation from Ref. [29], and calculated oscillations for a cooling ocean following an X-ray burst.The calculation focused on the "most likely suspects"; retrograde solutions of class (ii) with λ > 0 (excluding k = −1).Chambers and Watts [5] chose m = 1 and assumed large (positive) q.In the regime of large |q|, the eigenvalues of these mode-solutions asymptote to fixed values and become independent of the spin parameter (see Fig. 1 and Ref. [27]).Using these eigenvalues-effectively turning the problem into that of a single eigenvalue-Chambers and Watts [5] solved the relativistic analogues of Eqs.(11) and (15) to obtain the (rotating-frame) frequencies ω.Thus they obtained inertial-frame mode frequencies of σ ∼ 2 Hz and, by incorporating a model for the ocean cooling, observed drifts of ∼ 2 Hz, in reasonable agreement with (at least some of) the observations.Intuitively, we expect three relativistic effects to enter the mode problem for rotating stars.First, frequencies measured at infinity will involve the gravitational redshift.Second, the rotation will introduce frame-dragging which, in turn, also affects locally measured frequencies.Third, any variation in mass density associated with a wave will generate gravitational waves.Now, for dynamics in the star's low-density ocean, there will not be a significant gravitational-wave contribution.Moreover, as long as the ocean is shallow, we may account for the redshift and the frame-dragging through constant multiplicative factors.This considerably simplifies the problem.
Notably, for modes confined in a thin shell of fluid, the general-relativistic corrections to the frequencies can be quickly obtained following the strategy outlined in Ref. [49].Let us briefly discuss this strategy and apply it to the shallow ocean modes under examination here. 6or the spacetime of a slowly rotating star, it is generally convenient to work in the coordinate system given by the line element [50] ds 2 = −e ν(r) dt 2 + e λ(r) dr 2 − 2 (r)r 2 sin 2 θ dt dφ + r 2 (dθ 2 + sin 2 θ dφ 2 ). ( As we are interested in ocean modes, we specialise to a thin shell on the star's surface r = R.Because the framedragging is small ∼ O(Ω/Ω 0 ), the line element on a thin shell at the surface of the star can also be written as where the metric potentials ν R and R are approximated by their values at the surface-and it is worth noting that this is consistent with Ref. [5].It is then easy to see that, via the simple coordinate transformation the line element (45) takes the Minkowski form (on a thin shell and in spherical coordinates) Hence, we have transformed the coordinates to a (flat) Lorentz frame with r = R = const.This is a convenience of the shallowness of the layer we are considering.Moreover, the slow-rotation assumption implies that the fluid elements in the star move at non-relativistic velocities [50].The upshot to this is that we can take Newtonian results (provided they are confined to a thin shell, the case of interest here) and use the transformation (46) to derive the corresponding result in relativity.
Let us show how we can take advantage of this in practice.The first thing to note is that we first need to convert from rotating frame to inertial-frame frequencies as the line element (47) is of Minkowskian form.Taking the shallow-water result as a starting point, this means that the inertial-frame frequency is [by Eq. ( 32)] where we have explicitly introduced the "N" subscript to emphasise that the quantities are calculated assuming Newtonian gravity.Next, using the fact that the phase ωt + mφ is an invariant quantity, together with the transformation (46), we readily obtain the useful correspondence where, as before, ν R and R are the redshift factor and the frame-dragging evaluated at the surface of the star, respectively, and we have the "GR" label to denote relativistic quantities.The last step is to note that also the angular velocity as measured at infinity is going to change from Newtonian theory to relativity.Noting that this is defined as dφ/dt and using the transformation (46), we obtain Finally, using Eqs.( 48)-( 50), we arrive at the final result We see that the frame-dragging has no effect whatsoever, while the redshift factor enters explicitly.This result is consistent with what we might have expected from the beginning, as noted by, for example, Ref. [27].It is important to stress, however, that this is due to the fact that the rotating frame frequency in the large-limit does not depend on the spin frequency.It is not true in general.As an explicit demonstration, we may follow the same steps, this time instead starting from Eq. ( 25)-although we should (obviously) keep in mind that this result follows in the smalllimit, which does not seem to be appropriate for the systems we are interested in here.Anyway, in this case we arrive at Notably, the frame dragging now enters the relativistic expression explicitly, while the redshift factor does not.Finally, let us estimate the magnitude of the relativistic effects in the large-limit.For a fiducial neutron star of M = 1.4M and R = 10 km, we have so that, in order for the frequency offset to be approximately 1% we need Not surprisingly, accounting for general-relativistic corrections changes the Newtonian prediction by 30%, in good agreement with the numerical results from Chambers and Watts [5].

V. CONCLUDING REMARKS
The traditional approximation finds its roots in the classic geophysics literature and has also been used to calculate the dynamics of neutron-star oceans.However, the conditions under which the approximation holds are not particularly well understood.In an attempt to demystify it, we examined the implicit assumptions that underlie this approximation.For the low-frequency modes that we are interested in (in order to describe the small offsets in some X-ray detections), we found that the traditional approximation necessitates the neutron star to be rotating far slower than the observed systems.Thus, we suggest that the approximation may be inaccurate for the low-frequency oscillations of neutron-star oceans.
The central appeal of the traditional approximation is that it decouples the radial and angular dependencies of the perturbations.The angular sector of the perturbation problem becomes solely described by Laplace's tidal equation, associated with the eigenvalue λ, while the radial sector holds information about the stratification of the fluid and corresponds to the mode frequency ω.We discussed the eigenvalues of Laplace's tidal equation, which belong to two classes.Class (i) are the rotationally modified polar modes, which exist on the non-rotating star, and class (ii) include particular inertial modes, which limit to the r -modes in the slowly rotating regime.We noted the interesting behaviour of two solutions.(1) The k = 0 perturbation, previously known as a prograde Kelvin wave for q > 0, asymptotes rapidly to a minimum, which is contrast with other class (i) solutions.(2) At small λ, the k = −1 solution corresponds to the l = |m| r -mode.However, as |q| becomes large, this perturbation joins the class (i) solutions.
To relate to the classic literature, we discussed the shallow-water approximation, where the layer is assumed thin, uniform density and incompressible.Under these assumptions, the mode problem reduces to one eigenvalue, the Lamb parameter , and is solely angular.We examined the large-limit, appropriate for neutron-star oceans.We showed how the character of the class (ii) solutions (excluding the special k = −1 perturbation) change nature in this regime and no longer behave like inertial modes.This is supported by the standard mode calculation, which shows that the classic r -modes cannot exist in a shallow neutron-star ocean.
Finally, taking a step towards a more realistic setting for neutron stars, we demonstrated a simple prescription for lifting Newtonian results in a shallow layer to general relativity.While being a useful procedure for straightforwardly generalising Newtonian calculations, it also illustrates the importance of relativity.We found that relativistic corrections adjust calculations in Newtonian gravity by 30%.
Moving forward, an obvious improvement would be to ensure consistency within the traditional approximation, thereby making it an appropriate low-frequency filter.Fortunately, the resolution is basically trivial: One simply neglects the mode frequency in the radial component of the Euler equation (11).This will spoil the precise correspondence between the traditional approximation and the perturbation equations in spherical symmetry, but this seems like a small price to pay.To what extent this impacts on the numerical mode-solutions remains to be seen.
Additionally, it is worth recalling that one of the assumptions within the traditional approximation is to assume the fluid is rotating sufficiently slowly such that it may safely be approximated as spherical in shape.For bodies like our Earth, this is a perfectly reasonable constraint.However, we know that neutron stars can rotate very rapidly indeed.Along this vein, there has been work developing the perturbation equations in the traditional approximation to include the centrifugal deformation of the star [51,52].Other aspects of neutron-star physics will inevitably be important to incorporate into any detailed model of these systems.This includes (but is not limited to) the crust and magnetic field.However, in both of these contexts, we will need to depart from the traditional approximation.It would, of course, be interesting to progress these ideas further.
The issues we have discussed were important to raise, since our quantitative understanding of the natural oscillation modes of neutron-star oceans relies almost entirely on the traditional approximation.The interest in this direction has largely stemmed from using these mode solutions to explain the timing phenomenology of X-ray bursts on accreting neutron stars.Ultimately, we have shown that one should be careful in applying the traditional approximation to this context.We must also re-emphasise that a number of the observations exhibit bursts at the same frequency as the star's spin and a mode with a finite rotating-frame frequency will necessarily result in an offset.Given this tension, modes can describe some observations, but certainly not all.
Assuming we can Taylor expand the square-root this leads to the two roots ω 2 = λ(q)c 2