Waves in Thin Oceans on Oblate Neutron Stars
MMNRAS , 1–10 (2020) Preprint 12 June 2020 Compiled using MNRAS L A TEX style file v3.0
Waves in Thin Oceans on Oblate Neutron Stars
Bart F.A. van Baal, (cid:63) Frank R.N. Chambers and Anna L. Watts Anton Pannekoek Institute for Astronomy, University of Amsterdam, Postbus 94249, NL-1090 GE Amsterdam, the Netherlands
Accepted 2020 June 9. Received 2020 May 29; in original form 2020 February 21
ABSTRACT
Waves in thin fluid layers are important in various stellar and planetary problems.Due to rapid rotation such systems will become oblate, with a latitudinal variationin the gravitational acceleration across the surface of the object. In the case of ac-creting neutron stars, rapid rotation could lead to a polar radius smaller than theequatorial radius by a factor ∼ . . We investigate how the oblateness and a changinggravitational acceleration affect different hydrodynamic modes that exist in such fluidlayers through analytic approximations and numerical calculations. The wave vectorsof g -modes and Yanai modes increase for more oblate systems compared to sphericalcounterparts, although the impact of variations in the changing gravitational accel-eration is effectively negligible. We find that for increased oblateness, Kelvin modesshow less equatorial confinement and little change in their wave vector. For r -modes,we find that for more oblate systems the wave vector decreases. The exact manner ofthese changes for the r -modes depends on the model for the gravitational accelerationacross the surface. Key words: hydrodynamics - waves - stars: oscillations - stars: neutron - stars:rotation - X-rays: bursts
Waves in atmospheres and thin outer layers of stars arevery interesting in several astronomical applications. Previ-ous studies have mostly assumed spherical geometry; how-ever, rotating spheroids will become oblate. In particu-lar, accreting neutron stars can rotate up to several hun-dred times a second which could lead to eccentricities ofup to e = (cid:112) − c / a ≈ . , or alternatively a flattening f = − c / a ≈ . , where a is the semi-major axis and c thesemi-minor axis of the spheroid. In spite of this, the oblate-ness of neutron stars is not necessarily taken into consider-ation for all related phenomena. In this work, we considerthe effects of up to second order in the rotation frequency,meaning the neutron star can be modelled as an ellipsoid(Chandrasekhar 1969).The study of large-scale waves on a sphere comes witha very old and venerable history due to the applicationsto Earth’s oceans and atmosphere. Mathematically, thesewaves are described by Laplace’s tidal equation; a two-dimensional shallow-water model where the ocean is thinrelative to the radius of the sphere (Pedlosky 1987). Thesolutions to these equations consist of different families ofwaves which have different driving forces. One family are g -modes and arise due to the effects of buoyancy and theCoriolis force on displaced fluid packages (sometimes called (cid:63) E-mail: [email protected] gravito-inertial waves, or Poincar´e waves, Gill 1982; Ped-losky 1987, 2003). A different kind of low frequency grav-ity wave which instead arises from the conservation of spe-cific vorticity in combination with the stratification of theocean and are called Kelvin waves. Such waves are alwaysprograde and as such are sometimes referred to as low fre-quency prograde waves (Unno et al. 1989). The r -modes, orRossby waves, arise from the combination of the conserva-tion of specific vorticity combined with the curvature of thesurface and the variation of the Coriolis force that comeswith such curvatures. These waves can only be found on ro-tating spheres and are always retrograde. When the radialstructure of such modes are considered, they are often giventhe moniker ”buoyant” r -modes (Heyl 2004; Piro & Bildsten2005). Another wave was discovered due to pioneering workby Yanai & Maruyama (1966), called the Yanai mode af-ter this discovery. These modes are sometimes called mixedgravity-Rossby waves and mimic either the g -mode at lowfrequencies or r -mode at high frequencies (Gill 1982).Previous studies of large scale waves were done assum-ing a spherical geometry, which is a good approximation forthe Earth . In this work we are going to investigate the ef-fects of oblateness on these families of modes. We use the Our own Earth is also oblate due to rotation, although theflattening f = . × − is so low that it is often approximated asa sphere. © a r X i v : . [ a s t r o - ph . H E ] J un B.F.A. van Baal et al. work of Staniforth & White (2015) who derived a shallow-water model in a non-spherical geometry as was outlined byWhite & Wood (2012). We solve for linear perturbations intheir model and investigate the impact of rotationally in-duced oblateness on asymptotic approximations and numer-ical solutions for such perturbations, in a manner similar toTownsend (2003).An astrophysical example where the oblateness mightbe important is the phenomenon of thermonuclear burstoscillations, first discovered in X-ray observations byStrohmayer et al. (1996). These oscillations appear due toasymmetric brightness patches on the surface of accretingneutron stars during Type I thermonuclear X-ray bursts andare observed in about − of all bursts that are observedwith high time resolution instrumentation (Galloway et al.2008; Bilous & Watts 2019). Although it is currently un-clear what mechanism drives these asymmetries, global wavemodels were first suggested by Heyl (2004) as a possible ex-planation. The m = buoyant r -mode is a strong candidate.There have been concerns that this model overpredicts thechange in frequency throughout the course the burst (Piro &Bildsten 2005; Berkhout & Levin 2008), but recent calcula-tions that take into account more up to date nuclear physicsand relativistic effects (Chambers et al. 2019; Chambers &Watts 2020) show a reduced frequency change. A detailedoverview of burst oscillations, including their models, poten-tial origins and applications, is given by Watts (2012).These recent studies are the main motivation for thiswork, as the buoyant r -mode model for thermonuclear burstoscillations better matches observations when more realisticphysics is included. As previously mentioned, accreting neu-tron stars are expected to be oblate due to their rotationrates which range from Hz (IGR J17511-3057) to
Hz(4U 1608-522), but effects pertaining to this have not yetbeen included in the models .In Section 2 we give the mathematical derivation ofLaplace’s tidal equation for an oblate spheroid, and ourmodel for the gravitational acceleration across the surface;in Section 3 we investigate asymptotic approximations forthe solutions to these equations and compare to their spher-ical counterparts. In Section 4 we apply these new solutionsto the case of neutron star burst oscillations by investigatingthe r -modes in detail, before discussing the results in Section5. The strategy for our analysis is based on the work ofTownsend (2003), that investigated low-frequency pulsationmodes in rotating stars. By performing a separation of vari-ables on the perturbations used in the shallow-water equa-tions, it is possible to separate the partial differential equa-tions in spherical-polar coordinates into separate sets of or-dinary differential equations for the radial and latitudinaldependence of the perturbation. The important approxima-tion used is the so-called ‘traditional approximation’, which There is one source, IGR J17480-2446, which has a burst os-cillation frequency of Hz, but the magnetic field likely playsan important role for this source instead of rapid rotation (seeCavecchi et al. 2011). has its origins in geophysics (Eckart 1960) but has sincebeen applied to topics ranging from tidal forcing in mas-sive binaries (Papaloizou & Savonije 1997) to investigatinggravity modes on rotating neutron stars and their poten-tial link to quasi-periodic oscillations (Bildsten et al. 1996).This approximation amounts to neglecting the horizontalcomponent of the angular velocity vector when evaluatingthe Coriolis force (see e.g. Lee & Saio 1997, who discuss inwhich regimes this approximation is valid).The polar (colatitudinal, θ ) dependence is governed byLaplace’s tidal equations (Bildsten et al. 1996), which havea series of eigensolutions first discovered by Hough (1898).These Hough functions consist of a one-parameter familyof solutions in q = Ω / ω , where Ω is the angular velocityof the spheroid and ω the frequency of the wave, and havean eigenvalue λ . Only modes with positive λ will be consid-ered, although negative values of λ are possible as convectivemodes stabilized by the Coriolis force (Lee & Saio 1986).Generally, a numerical approach is necessary to find thesolutions to Laplace’s tidal equations. However, in the limitthat | q | becomes large, asymptotic approximations can bederived. This is well-established in the geophysics literature(see e.g. Longuet-Higgins 1968; Gill 1982), but the work byLee & Saio (1987) was the first time such approximationswere applied to stellar non-radial pulsations, albeit in onlythe radial direction. It was Bildsten et al. (1996) who firstpresented the asymptotic solutions to the angular equations,and Townsend (2003) who formally derived this approxima-tion. However, the effects of rapid rotation (large Ω ) werenot included even though rapid rotation is necessary in orderto reach the limit of large | q | .In this work we consider the effects of oblateness in-duced by rotation, which alters the problem in two mainways: firstly, the gravitational acceleration will vary withlatitude across the surface; the second, and more complex,effect is that by including the oblateness, the radial coor-dinate r is no longer constant for varying latitude, compli-cating the separation of variables used to solve the problem.Staniforth & White (2015) derived a set of shallow-waterequations in a zonally symmetric geometry, which allows fora latitudinal variation of gravity, known to be significantfor even slowly rotating neutron stars (AlGendy & Morsink2014).The separation of the radial and polar coordinates iskey in calculating the waves, and by using the parametricellipsoidal coordinates from Staniforth & White (2015) theequations can be written in such a way that the radial andpolar coordinates can be separated once again. This meansthat a new ODE can be recovered and we can follow manyof the steps taken by Townsend (2003), although there willbe new parameters introduced to track the eccentricity andvariation in gravitational acceleration which are present dueto the oblateness. Staniforth & White (2015) derive a set of equations appro-priate for an incompressible ocean that exists on the surfaceof a spheroid uniformly rotating about the z -axis with angu-lar velocity Ω . These shallow water equations describe thefluid flow in terms of azimuthal and latitudinal componentsof velocity ( u φ and u θ respectively) and the height of the MNRAS000
Hz(4U 1608-522), but effects pertaining to this have not yetbeen included in the models .In Section 2 we give the mathematical derivation ofLaplace’s tidal equation for an oblate spheroid, and ourmodel for the gravitational acceleration across the surface;in Section 3 we investigate asymptotic approximations forthe solutions to these equations and compare to their spher-ical counterparts. In Section 4 we apply these new solutionsto the case of neutron star burst oscillations by investigatingthe r -modes in detail, before discussing the results in Section5. The strategy for our analysis is based on the work ofTownsend (2003), that investigated low-frequency pulsationmodes in rotating stars. By performing a separation of vari-ables on the perturbations used in the shallow-water equa-tions, it is possible to separate the partial differential equa-tions in spherical-polar coordinates into separate sets of or-dinary differential equations for the radial and latitudinaldependence of the perturbation. The important approxima-tion used is the so-called ‘traditional approximation’, which There is one source, IGR J17480-2446, which has a burst os-cillation frequency of Hz, but the magnetic field likely playsan important role for this source instead of rapid rotation (seeCavecchi et al. 2011). has its origins in geophysics (Eckart 1960) but has sincebeen applied to topics ranging from tidal forcing in mas-sive binaries (Papaloizou & Savonije 1997) to investigatinggravity modes on rotating neutron stars and their poten-tial link to quasi-periodic oscillations (Bildsten et al. 1996).This approximation amounts to neglecting the horizontalcomponent of the angular velocity vector when evaluatingthe Coriolis force (see e.g. Lee & Saio 1997, who discuss inwhich regimes this approximation is valid).The polar (colatitudinal, θ ) dependence is governed byLaplace’s tidal equations (Bildsten et al. 1996), which havea series of eigensolutions first discovered by Hough (1898).These Hough functions consist of a one-parameter familyof solutions in q = Ω / ω , where Ω is the angular velocityof the spheroid and ω the frequency of the wave, and havean eigenvalue λ . Only modes with positive λ will be consid-ered, although negative values of λ are possible as convectivemodes stabilized by the Coriolis force (Lee & Saio 1986).Generally, a numerical approach is necessary to find thesolutions to Laplace’s tidal equations. However, in the limitthat | q | becomes large, asymptotic approximations can bederived. This is well-established in the geophysics literature(see e.g. Longuet-Higgins 1968; Gill 1982), but the work byLee & Saio (1987) was the first time such approximationswere applied to stellar non-radial pulsations, albeit in onlythe radial direction. It was Bildsten et al. (1996) who firstpresented the asymptotic solutions to the angular equations,and Townsend (2003) who formally derived this approxima-tion. However, the effects of rapid rotation (large Ω ) werenot included even though rapid rotation is necessary in orderto reach the limit of large | q | .In this work we consider the effects of oblateness in-duced by rotation, which alters the problem in two mainways: firstly, the gravitational acceleration will vary withlatitude across the surface; the second, and more complex,effect is that by including the oblateness, the radial coor-dinate r is no longer constant for varying latitude, compli-cating the separation of variables used to solve the problem.Staniforth & White (2015) derived a set of shallow-waterequations in a zonally symmetric geometry, which allows fora latitudinal variation of gravity, known to be significantfor even slowly rotating neutron stars (AlGendy & Morsink2014).The separation of the radial and polar coordinates iskey in calculating the waves, and by using the parametricellipsoidal coordinates from Staniforth & White (2015) theequations can be written in such a way that the radial andpolar coordinates can be separated once again. This meansthat a new ODE can be recovered and we can follow manyof the steps taken by Townsend (2003), although there willbe new parameters introduced to track the eccentricity andvariation in gravitational acceleration which are present dueto the oblateness. Staniforth & White (2015) derive a set of equations appro-priate for an incompressible ocean that exists on the surfaceof a spheroid uniformly rotating about the z -axis with angu-lar velocity Ω . These shallow water equations describe thefluid flow in terms of azimuthal and latitudinal componentsof velocity ( u φ and u θ respectively) and the height of the MNRAS000 , 1–10 (2020) blate Thin Ocean Waves free surface of the fluid layer ( H , measured from B ( θ ) whichis a rigid base). They are written:D hor u φ D t + (cid:18) u φ a sin θ + Ω (cid:19) u θ cos θσ ( θ ) + a sin θ ∂∂φ (cid:20) g ( θ ) H (cid:21) = , (1a)D hor u θ D t − (cid:18) u φ a sin θ + Ω (cid:19) u φ cos θσ ( θ ) + a σ ( θ ) ∂∂θ (cid:20) g ( θ ) H (cid:21) = , (1b)D hor [ H − B ( θ )] D t + H − B ( θ ) a sin θ (cid:20) ∂ u φ ∂φ + σ ( θ ) ∂∂θ (cid:0) u θ sin θ (cid:1)(cid:21) = , (1c)whereD hor D t ≡ ∂∂ t + u φ a sin θ ∂∂φ + u θ a σ ( θ ) ∂∂θ (2)is the horizontal material derivative . The eccentricity ofthe ellipse is e ≡ (cid:112) ( − c / a ) where a and c are the semi-major and semi-minor axes of the spheroid respectively. Thefunction σ ( θ ) depends on co-latitude and the eccentricity as √ − e sin θ . The surface B ( θ ) will be assumed to form ageopotential surface, thus ∂ θ [ g ( θ ) B ( θ )] = . See Staniforth &White (2015) for the full set of assumptions used in derivingthese equations.The gravitational acceleration on the surface of the el-lipsoid, g ( θ ) , is not constant, which is a departure from theshallow water equations on the surface of a sphere. The de-pendence of the gravitational acceleration on the co-latitudeis discussed further in Section 2.2. The product of the grav-itational acceleration and the height of the free surface, g H ,is constant along the surface of the ellipsoid; this result is re-lated to the fact that the pressure is assumed to be constantalong the free surface.We solve for perturbations of the form: g H = g B + g H + (cid:15) P p ( θ ) e i ( m φ + ω t ) , (3a) sin θ u θ = i (cid:15)ω P θ ( θ ) e i ( m φ + ω t ) , (3b) sin θ u φ = (cid:15)ω P φ ( θ ) e i ( m φ + ω t ) , (3c)where ω is the wave frequency, m denotes the azimuthal wavenumber (which is constrained to integer values), and (cid:15) is asmall factor used to keep track of perturbed terms. H isthe height of the fluid layer at rest. The Lagrangian fluiddisplacement is related to the fluid velocity perturbation as u ≡ D hor ξ / D t . The pressure along the free surface is relatedto the height as p = p + ρ g H where p is a reference surface.The wave propagates either prograde or retrograde depend-ing on the signs of m and ω . Prograde (retrograde) motioncorresponds to m ω < ( m ω > ).We apply the perturbations in Equations (3) to Equa-tions (1), retaining terms linear in (cid:15) . We define the param-eters q ≡ Ω / ω and λ ≡ a /( H g ρ ) and use the substitution Staniforth & White (2015), equations (48)-(51), are written interm of latitude, θ l , whereas we use to co-latitude, θ . Convertingbetween these coordinates introduces a sign change in the north-south velocity. Using their notation, the latitudinal velocity isrelated to the co-latitudinal velocity as u = − u θ . The eigenvalue λ is related to the wave vector k as k ∼ λ / / R in spherical geometry. µ ≡ cos ( θ ) and D ≡ ( − µ ) d / d µ to find the following set ofordinary differential equations: D P p = − σ P θ − q µ P φ , (4a) mP p = − P φ − q µ σ P θ , (4b) λσ ( − µ ) P p = − σ mP φ + D P θ − P θ g D g . (4c)Note that these equations reduce to their spherical counter-part when e = and the gravitational acceleration becomesconstant (equations (17)-(19) in Townsend 2003). The quan-tity P φ can be eliminated and the equations become: (cid:18) D − mq µ (cid:19) P p = σ (cid:18) q µ σ − (cid:19) P θ , (5) (cid:18) D + mq µ (cid:19) P θ − P θ D ln g = σ (cid:20) λ ( − µ ) − m (cid:21) P p . (6)It is possible to eliminate P θ and find one second order ordi-nary differential equation for P p which reduces to Laplace’stidal equation in the limit that e = and the gravitationalacceleration becomes constant. Further setting q = , theequation reduces to the associated Legendre equation. Including oblateness introduces an extra term in Equation(6) compared to the equation in spherical geometry. We re-quire a model for the variation in gravitational accelerationwith latitude across the surface of the spheroid. Both the ec-centricity e and any model for g ( θ ) are functions of the mass,radius, rotation rate and equation of state of the spheroid.Any model for g ( θ ) must be symmetrical about theequator. We choose to model this function as a polynomialin cos θ , g ( µ ) ≡ g E ( + χµ ) , (7)where g E is the gravitational acceleration at the equator,and χ is a new parameter which describes how gravity varieswith latitude. The gravitational acceleration at the pole is g E ( + χ ) . Through this equation, instead of dealing withparameters for radii, masses, rotation rates and equations ofstate, we can describe our variations through e and χ alone.This adds only two extra parameters to the shallow-waterequations. Equation (7) is just one way of describing g ( θ ) ,and more complicated choices are certainly possible.We now continue by making asymptotic approxima-tions, similar to Townsend (2003), while including these newparameters that account for oblateness. In this section we find asymptotic approximations for theeigenvalues and eigenfunctions of Equations (5) and (6) forlarge spin parameter and a range of values for eccentricity.We compare these approximations to their spherical coun-terparts and numerical solutions.
MNRAS , 1–10 (2020)
B.F.A. van Baal et al.
Yoshida (1960) recognised that for large values of thespin parameter, q , certain families of modes become equato-rially trapped due to the Coriolis force. This equatorial trap-ping was the principal property used by Townsend (2003) toderive asymptotic solutions, where it was noted that thespin parameter always appears together with the latitudinalcoordinate in Laplace’s tidal equation. Since eigenfunctionsmust remain finite over the whole range of latitude, they canonly be significantly different from zero around the equatorwhere terms q µ are or order unity. Similar arguments aretrue for these equations when the effect of oblateness is in-cluded. This leads to the simplification, D P i ≈ d P i / d µ . g -modes,Yanai modes and r -modes By assuming that λ is appreciably different from m , wemay also drop terms of order λµ P on the right hand sideof Equation (6). These two simplifications lead to the equa-tions: (cid:18) dd µ − mq µ (cid:19) P p = σ (cid:18) q µ σ − (cid:19) P θ , (8) (cid:18) dd µ + mq µ (cid:19) P θ − P θ dd µ ln g = σ (cid:18) λ − m (cid:19) P p . (9) P p can be eliminated to obtain a second-order differentialequation for P θ :d P θ d µ − (cid:18) d ln g d µ + e µσ (cid:19) d P θ d µ + (cid:18) mq − λ q µ + σ λ − σ m − e µ mq σ − d ln g d µ + mq µ d ln g d µ + e µσ d ln g d µ (cid:19) P θ = . (10)We now implement our model for the gravitational acceler-ation across the surface of the star, discussed in Section 2.2.Again, we take advantage of the fact that the eigenfunction isonly appreciably different from zero in the region around theequator. This leads to the approximations d µ ln g P θ ≈ χµ P θ and d µ ln g P θ ≈ χ P θ . For the same reason, we also make thesimplification that σ ( µ ) P i ≈ √ − e P i , which can be used tosimplify Equation (10):d P θ d µ − µ (cid:18) χ + e σ (cid:19) d P θ d µ + (cid:18) mq − λ q µ + σ λ − σ m − e µ mq σ + χ mq µ + χ e µ σ − χ (cid:19) P θ = . (11)In order to simplify the analysis of this equation, we intro-duce the following variables: L ≡ λ, (12) Υ ≡ e σ mq − χ mq − χ e σ , (13) α ≡ ( L q + Υ ) / µ, (14) A ≡ mq + σ ( L − m ) − χ ( L q + Υ ) / , (15) E ≡ − χ + e σ ( L q + Υ ) / . (16)With these new definitions comes an extra constraint; α must be real and thus L q (cid:62) − Υ . With the correct sub-stitutions and replacing µ as a variable with α , the ODEcan be written as:d P θ d α + α E d P θ d α + ( A − α ) P θ = . (17)In the case that E = this equation reduces to the time-independent Schr¨odinger equation for a quantum harmonicoscillator for a particle trapped in a potential well (Arfken& Weber 1999). In a similar manner, waves are trapped inthe region around the equator; E and Υ are responsible forthe effect of oblateness which alters this potential.A general solution to Equation (17) is given by a com-bination of a Hermite Polynomial ( H s ) and Kummer’s Con-fluent Hypergeometric Function ( F ). Although there is arelation between H s and F , fortunately we need not worryabout this because the boundary conditions of the problemsimplify this solution dramatically; either P p or d P p / d µ arezero at the equator depending on the parity of the mode(Bildsten et al. 1996). Since F will not be zero at thesepoints, this part of the general solution must vanish in or-der to satisfy the boundary conditions. Thus, the generalsolution is given by, P θ ( α ) = e − (cid:0) √ E + + E (cid:1) α H s (cid:18) α √ (cid:112) E + (cid:19) . (18)With this specific solution comes the integer index s whichis given by: s = − √ E + + E − A √ E + , (19)where the integer s (cid:62) . From Equation (19) the followingrelation can be recovered: A − AE − s ( s + ) E = ( s + ) . (20)Combining this relation with Equations (15) and (16), anexpression quadratic in λ can be obtained. From there it ispossible to find the two roots for λ ± : L ± = − mq σ − m σ − χσ σ − χσ + e σ + ( s + ) q σ (cid:40) ± (cid:20) − ( mq σ − m σ − χσ )( s + ) q − ( χσ + e )( s + ) q + σ Υ ( s + ) q + ( χσ + e ) ( s + ) q (cid:21) / (cid:41) , (21)and by using a Taylor expansion these roots can be approx-imated as λ + = ( s + ) q ( − e ) − (cid:20) mq ( − e ) − m − χ ( − e ) (cid:21) , (22)and λ − = (cid:2) mq − m ( − e ) (cid:3) ( s + ) q + mq (cid:20) e ( − e ) − χ (cid:21) (cid:2) + ( s + ) (cid:3) ( s + ) q . MNRAS000
Yoshida (1960) recognised that for large values of thespin parameter, q , certain families of modes become equato-rially trapped due to the Coriolis force. This equatorial trap-ping was the principal property used by Townsend (2003) toderive asymptotic solutions, where it was noted that thespin parameter always appears together with the latitudinalcoordinate in Laplace’s tidal equation. Since eigenfunctionsmust remain finite over the whole range of latitude, they canonly be significantly different from zero around the equatorwhere terms q µ are or order unity. Similar arguments aretrue for these equations when the effect of oblateness is in-cluded. This leads to the simplification, D P i ≈ d P i / d µ . g -modes,Yanai modes and r -modes By assuming that λ is appreciably different from m , wemay also drop terms of order λµ P on the right hand sideof Equation (6). These two simplifications lead to the equa-tions: (cid:18) dd µ − mq µ (cid:19) P p = σ (cid:18) q µ σ − (cid:19) P θ , (8) (cid:18) dd µ + mq µ (cid:19) P θ − P θ dd µ ln g = σ (cid:18) λ − m (cid:19) P p . (9) P p can be eliminated to obtain a second-order differentialequation for P θ :d P θ d µ − (cid:18) d ln g d µ + e µσ (cid:19) d P θ d µ + (cid:18) mq − λ q µ + σ λ − σ m − e µ mq σ − d ln g d µ + mq µ d ln g d µ + e µσ d ln g d µ (cid:19) P θ = . (10)We now implement our model for the gravitational acceler-ation across the surface of the star, discussed in Section 2.2.Again, we take advantage of the fact that the eigenfunction isonly appreciably different from zero in the region around theequator. This leads to the approximations d µ ln g P θ ≈ χµ P θ and d µ ln g P θ ≈ χ P θ . For the same reason, we also make thesimplification that σ ( µ ) P i ≈ √ − e P i , which can be used tosimplify Equation (10):d P θ d µ − µ (cid:18) χ + e σ (cid:19) d P θ d µ + (cid:18) mq − λ q µ + σ λ − σ m − e µ mq σ + χ mq µ + χ e µ σ − χ (cid:19) P θ = . (11)In order to simplify the analysis of this equation, we intro-duce the following variables: L ≡ λ, (12) Υ ≡ e σ mq − χ mq − χ e σ , (13) α ≡ ( L q + Υ ) / µ, (14) A ≡ mq + σ ( L − m ) − χ ( L q + Υ ) / , (15) E ≡ − χ + e σ ( L q + Υ ) / . (16)With these new definitions comes an extra constraint; α must be real and thus L q (cid:62) − Υ . With the correct sub-stitutions and replacing µ as a variable with α , the ODEcan be written as:d P θ d α + α E d P θ d α + ( A − α ) P θ = . (17)In the case that E = this equation reduces to the time-independent Schr¨odinger equation for a quantum harmonicoscillator for a particle trapped in a potential well (Arfken& Weber 1999). In a similar manner, waves are trapped inthe region around the equator; E and Υ are responsible forthe effect of oblateness which alters this potential.A general solution to Equation (17) is given by a com-bination of a Hermite Polynomial ( H s ) and Kummer’s Con-fluent Hypergeometric Function ( F ). Although there is arelation between H s and F , fortunately we need not worryabout this because the boundary conditions of the problemsimplify this solution dramatically; either P p or d P p / d µ arezero at the equator depending on the parity of the mode(Bildsten et al. 1996). Since F will not be zero at thesepoints, this part of the general solution must vanish in or-der to satisfy the boundary conditions. Thus, the generalsolution is given by, P θ ( α ) = e − (cid:0) √ E + + E (cid:1) α H s (cid:18) α √ (cid:112) E + (cid:19) . (18)With this specific solution comes the integer index s whichis given by: s = − √ E + + E − A √ E + , (19)where the integer s (cid:62) . From Equation (19) the followingrelation can be recovered: A − AE − s ( s + ) E = ( s + ) . (20)Combining this relation with Equations (15) and (16), anexpression quadratic in λ can be obtained. From there it ispossible to find the two roots for λ ± : L ± = − mq σ − m σ − χσ σ − χσ + e σ + ( s + ) q σ (cid:40) ± (cid:20) − ( mq σ − m σ − χσ )( s + ) q − ( χσ + e )( s + ) q + σ Υ ( s + ) q + ( χσ + e ) ( s + ) q (cid:21) / (cid:41) , (21)and by using a Taylor expansion these roots can be approx-imated as λ + = ( s + ) q ( − e ) − (cid:20) mq ( − e ) − m − χ ( − e ) (cid:21) , (22)and λ − = (cid:2) mq − m ( − e ) (cid:3) ( s + ) q + mq (cid:20) e ( − e ) − χ (cid:21) (cid:2) + ( s + ) (cid:3) ( s + ) q . MNRAS000 , 1–10 (2020) blate Thin Ocean Waves (23)Equations (22) and (23) reduce to the forms given byTownsend (2003) in spherical geometry.The two different branches of eigenvalues are associatedwith different types of waves. The λ + branch are known as g -modes, sometimes called Poincar´e waves (Gill 1982), whilethe λ − branch is associated with r -modes, since all valid solu-tions must be non-axisymmetric and retrograde (Saio 1982).Townsend (2003) furthermore finds solutions for a specialcase where s = , sometimes called the mixed gravity-Rossbywave since it can create both types of waves depending onwhether the mode is retrograde or prograde. However, thesekinds of waves were discovered by Yanai & Maruyama (1966)and as such are also known as Yanai waves in honour of thisdiscovery. In the limit of large | q | , the Yanai modes are bestapproximated as g -modes with s = in the case of eitherprograde or retrograde propagation. We now connect the value of the parameter dictating thegravitational acceleration across the surface of the star, χ ,with the eccentricity, e . We make specific assumptions forthe case of neutron star oceans using the models establishedby AlGendy & Morsink (2014). In order to investigate arelation between χ and e , we explore a range of values forcompactness, angular velocity, mass and radius .We are particularly interested in neutron stars withthermonuclear burst oscillations, which are known to havespin frequencies as high as Hz. For neutron stars rotat-ing at these high frequencies, e and χ follow the approximaterelation: χ = e . (24)Assuming a neutron star with equatorial radius between and . km, and mass between . and . M (cid:12) for a varietyof equations of state, we find e lies in the range e = . − . or e = . − . for a rotation rate of Hz and
Hz,respectively. In order to more thoroughly investigate the ef-fect of a larger variation of gravitational acceleration fromthe equator to the pole, we will also consider the relation χ = e . With new asymptotic approximations to eigenvalues includ-ing oblateness, we now investigate how well they compare tonumerical solutions, and how increasing oblateness changesthe values compared to their spherical counterparts as cal-culated by Townsend (2003). We use a shooting method tocalculate the eigenvalues and eigenfunctions, and use thesame normalisation condition as in Townsend (2003).In the non-rotating limit P p reduces to the associatedLegendre polynomials P ml (Lee & Saio 1997) (in the case of g -modes, as r -modes do not exist in non-rotating systems). See Equations (20) and (49) in AlGendy & Morsink (2014) forradius and gravitational acceleration, respectively.
Figure 1.
A comparison between the numerically calculated(black) and approximate (in colour) solutions to the m = − pro-grade and retrograde waves with k = for different eccentricities.Orange corresponds to a prograde g -mode and red to a retrograde g -mode. The solid lines correspond to e = , the dotted lines to e = . and the dashed lines to e = . , with the relation between e and χ as given by Equation (24). The analytic approximationsare given by Equation (22) for both modes. Figure 2.
A comparison between the numerically calculated(black) and approximate (in colour) solutions to the m = − pro-grade and retrograde waves with k = for different eccentricities.Purple corresponds to a retrograde g -mode and green to a pro-grade Yanai mode. The solid lines correspond to e = , the dottedlines to e = . and the dashed lines to e = . , with the rela-tion between e and χ as given by Equation (24). The analyticapproximations are given by Equation (22) for both modes. For these polynomials, m is the azimuthal order and the in-teger l (cid:62) is the harmonic degree. However, with this classi-fication it is impossible to describe retrograde Yanai modesand r -modes, as those do not exist in the non-rotating limit.Lee & Saio (1997) came up with an alternative classificationscheme which uses a unique integer index k for each solutionof the tidal equations. Positive and zero values of this indexindicates that the mode possesses a counterpart of harmonicdegree l = | m | + k , in the limit of no rotation, while nega-tive values of k denote r -modes and retrograde Yanai modes.Furthermore, odd values of k denote odd-parity modes whileeven values of k correspond to even-parity modes.It is straightforward to relate the k -index used by Lee &Saio (1997) with the s -index which appears in the approxi-mations of Equations (22) and (23), which was also used byTownsend (2003). For prograde modes, the relation is givenby s = k − , while for the retrograde modes there is a split, MNRAS , 1–10 (2020)
B.F.A. van Baal et al.
Figure 3.
A comparison between the numericallu calculated(black) and approximate (in colour) solutions to the m = − pro-grade and retrograde waves with k = for different eccentricities.Dark blue corresponds to a retrograde g -mode and light blue tothe prograde Kelvin mode. The solid lines correspond to e = , thedotted lines to e = . and the dashed lines to e = . , with therelation between e and χ as given by Equation (24). The analyticapproximations are given by Equations (22) for the retrograde(negative q ) and (29) for the prograde mode, respectively. Figure 4.
A comparison between the numerically calculated(black) and approximate (in colour) solutions to the m = − ret-rograde Yanai ( k = − ) and r -mode ( k = − ) waves for differenteccentricities of e = (solid lines), e = . (dotted) and e = . (dashed). The Yanai mode (green) is approximated by Equation(22) and matches well for | q | > . For increasing eccentricity,the analytic approximations of the r -mode (brown) do not matchwith the numerical solution. The approximation for the eigenval-ues are calculated using Equation (22) for the Yanai mode andEquation (23) for the r -mode. The relation between e and χ isgiven by Equation (24). as the relation is s = − k − for the r -modes and s = k + for the g -modes; both of these result in the correct relationfor the retrograde Yanai mode (which has k = − ⇔ s = ).The advantage of using the k -index over s is that the k -indexing scheme can uniquely identify all modes, while forthe s scheme additional information would be necessary inorder to differentiate between g -modes and r -modes.In Figures 1, 2, 3 and 4 we show the numerical solutionand analytic approximation to the eigenvalues for the m = − and k = , k = , k = and k = − , − modes. In each ofthe Figures, three different eccentricities e = , . , . areused. The other parameter we have to choose is χ , for whichthe relation in Equation (24) is used. From Equation (22) however, it can be seen that if | q | is large then the impactof χ will be marginal compared to the impact of e .For the k = modes (which are g -modes) in Figure 1it can be seen that beyond | q | = the analytic approxi-mation matches well with the numerical solution for all ec-centricities. For both the prograde and retrograde g -modes,the eigenvalues increase with q . Higher eccentricities (dottedand dashed lines) results in a steeper slope for λ against q ( λ increases more rapidly when e goes from . to . thanit does from to . ). For the retrograde g -mode in Figure2 the same effects can be observed. There is a good matchbetween the numerical and approximate values at | q | = .Eigenvalues scale with e as predicted by the asymptotic ap-proximation. On the prograde side of the Figure however, itcan be seen that the Yanai mode converges more slowly tothe approximation of Equation (22), but it does also followthe same trend of scaling with e .In Figure 3, the retrograde g -mode repeats the be-haviour of the other g -modes seen for k = , . At | q | = theanalytic approximations and the numerical solution agree,and for increasing e the eigenvalues increase. For the ret-rograde Yanai and the r -modes, we need to look at highervalues of | q | as there are no positive solutions for lower val-ues of q for these modes. In Figure 4 a comparison is shownbetween the numerical solutions and approximations to theeigenvalues. For the retrograde Yanai mode, the numericalsolutions and analytic approximations agree for values of | q | > , but converge quite slowly.Unlike the other wave families, the analytic approxima-tion for the r -mode does not match with the numerical solu-tion for higher eccentricity. It can be seen that for the zeroeccentricity (solid lines) the approximation and numericalsolutions do converge to the same value for high | q | , and theanalytic approximation changes little for increasing e . Thisis not the same behaviour shown by the numerical solutions.A higher value of | q | is required to find positive eigenvalues,and the asymptotic limit which these eigenvalues approachis much smaller for increasing e .We test the effect of χ on the eigenvalues by using adifferent relation between e and χ , such as χ = e . For the g -modes and Yanai modes, any differences are negligible,which is not unexpected given that terms with χ scale with q , while terms with e scale with q in the approximation for λ in Equation (22). We show the impact for the r -mode inSection 4. The final family of modes which we will investigate are theKelvin modes. In the spherical case, the eigenvalue of thesemodes may be approximated by λ ≈ m in the limit of q → ∞ (arbitrary but rapid rotation). We make the same approxi-mations associated with equatorial trapping for these modesas was performed for the g -modes and r -modes. However,we cannot drop the µ dependence on the right hand side ofEquation (26) since λ and m are not appreciably different.Furthermore, this fact implies that P p must be much largerthan P θ since we require the term on the right hand sideof the equation ( ∼ λµ P p ) to be the same order as mq µ P θ . MNRAS000
A comparison between the numerically calculated(black) and approximate (in colour) solutions to the m = − ret-rograde Yanai ( k = − ) and r -mode ( k = − ) waves for differenteccentricities of e = (solid lines), e = . (dotted) and e = . (dashed). The Yanai mode (green) is approximated by Equation(22) and matches well for | q | > . For increasing eccentricity,the analytic approximations of the r -mode (brown) do not matchwith the numerical solution. The approximation for the eigenval-ues are calculated using Equation (22) for the Yanai mode andEquation (23) for the r -mode. The relation between e and χ isgiven by Equation (24). as the relation is s = − k − for the r -modes and s = k + for the g -modes; both of these result in the correct relationfor the retrograde Yanai mode (which has k = − ⇔ s = ).The advantage of using the k -index over s is that the k -indexing scheme can uniquely identify all modes, while forthe s scheme additional information would be necessary inorder to differentiate between g -modes and r -modes.In Figures 1, 2, 3 and 4 we show the numerical solutionand analytic approximation to the eigenvalues for the m = − and k = , k = , k = and k = − , − modes. In each ofthe Figures, three different eccentricities e = , . , . areused. The other parameter we have to choose is χ , for whichthe relation in Equation (24) is used. From Equation (22) however, it can be seen that if | q | is large then the impactof χ will be marginal compared to the impact of e .For the k = modes (which are g -modes) in Figure 1it can be seen that beyond | q | = the analytic approxi-mation matches well with the numerical solution for all ec-centricities. For both the prograde and retrograde g -modes,the eigenvalues increase with q . Higher eccentricities (dottedand dashed lines) results in a steeper slope for λ against q ( λ increases more rapidly when e goes from . to . thanit does from to . ). For the retrograde g -mode in Figure2 the same effects can be observed. There is a good matchbetween the numerical and approximate values at | q | = .Eigenvalues scale with e as predicted by the asymptotic ap-proximation. On the prograde side of the Figure however, itcan be seen that the Yanai mode converges more slowly tothe approximation of Equation (22), but it does also followthe same trend of scaling with e .In Figure 3, the retrograde g -mode repeats the be-haviour of the other g -modes seen for k = , . At | q | = theanalytic approximations and the numerical solution agree,and for increasing e the eigenvalues increase. For the ret-rograde Yanai and the r -modes, we need to look at highervalues of | q | as there are no positive solutions for lower val-ues of q for these modes. In Figure 4 a comparison is shownbetween the numerical solutions and approximations to theeigenvalues. For the retrograde Yanai mode, the numericalsolutions and analytic approximations agree for values of | q | > , but converge quite slowly.Unlike the other wave families, the analytic approxima-tion for the r -mode does not match with the numerical solu-tion for higher eccentricity. It can be seen that for the zeroeccentricity (solid lines) the approximation and numericalsolutions do converge to the same value for high | q | , and theanalytic approximation changes little for increasing e . Thisis not the same behaviour shown by the numerical solutions.A higher value of | q | is required to find positive eigenvalues,and the asymptotic limit which these eigenvalues approachis much smaller for increasing e .We test the effect of χ on the eigenvalues by using adifferent relation between e and χ , such as χ = e . For the g -modes and Yanai modes, any differences are negligible,which is not unexpected given that terms with χ scale with q , while terms with e scale with q in the approximation for λ in Equation (22). We show the impact for the r -mode inSection 4. The final family of modes which we will investigate are theKelvin modes. In the spherical case, the eigenvalue of thesemodes may be approximated by λ ≈ m in the limit of q → ∞ (arbitrary but rapid rotation). We make the same approxi-mations associated with equatorial trapping for these modesas was performed for the g -modes and r -modes. However,we cannot drop the µ dependence on the right hand side ofEquation (26) since λ and m are not appreciably different.Furthermore, this fact implies that P p must be much largerthan P θ since we require the term on the right hand sideof the equation ( ∼ λµ P p ) to be the same order as mq µ P θ . MNRAS000 , 1–10 (2020) blate Thin Ocean Waves This leads to the set of equations: (cid:18) dd µ − mq µ (cid:19) P p = , (25) (cid:18) dd µ + mq µ (cid:19) P θ − P θ d ln g d µ = σ (cid:20) λ ( − µ ) − m (cid:21) P p , (26)Equation (25) is unchanged from spherical geometry, allow-ing the same solution: P p ( µ ) = e mq µ / . (27)We require mq < (prograde motion) since the mode mustdecay to zero at the pole. In the astrophysical literaturethey have, because of this, been referred to as low-frequencyprograde waves (Unno et al. 1989). Applying the solutionfor P p to Equation (26) we find the following first orderdifferential equation:dd µ (cid:18) e mq µ / P θ (cid:19) − e mq µ / P θ d ln g d µ = σ (cid:20) λ ( − µ ) − m (cid:21) e mq µ . (28)This equation can be solved using integrating factors. How-ever, we find that the solution for the eigenvalues in sphericalgeometry is valid for even large values of eccentricity. Thissolution is: λ = m mq mq + , (29)which satisfies the initial assumption that λ ≈ m .In Figure 3 the numerical solutions are shown for dif-ferent eccentricities e and compared to the analytic approx-imation of Equation (29). For all eccentricities the analyticapproximation is a good match, which shows that the Kelvinmodes on oblate spheroids can still be approximated usingthe spherical analysis. For increasing eccentricities a slightdecrease in the eigenvalue can be observed for this mode, al-though the difference from the spherical eigenvalue is small.We test the effect of χ on the eigenvalues by using adifferent relation between e and χ , namely χ = e . For theKelvin mode, the numerical solutions decrease slightly morefor this new relation but still match well to the sphericalapproximation.The eigenfunctions of the Kelvin wave depend on theoblateness in an interesting manner. Figure 5 shows that forincreasing eccentricity the peak of the P θ function moves tohigher latitudes, indicating that the wave is less confined tothe equatorial region for more oblate systems. This sets theKelvin modes apart from the g -modes and Yanai modes, asthose become more equatorially confined in more eccentricsystems. We find that even for the most oblate systems, thenew term in equation (26), related to gravitational variationacross the surface of the star, is small relative to the otherterms which appear in that equation. The lack of influenceof gravitational variation on the Kelvin wave is also reflectedin the minor differences in eigenvalues when including thiseffect. Rotationally induced oblateness is relevant for neutron stars,in particular for surface modes which may be responsible for
Figure 5.
A comparison between the numerically calculatedeigenfunctions P p , P θ and P φ for the m = − Kelvin mode at q = , for different eccentricities of e = (solid lines), e = . (dotted) and e = . (dashed). For more eccentric systems, theequatorial confinement of the Kelvin mode is weaker. The peakof P θ is found at higher latitudes. This behaviour is only observedfor the Kelvin mode, as the other modes all have stronger equato-rial confinement for more eccentric systems. The relation between e and χ is given by Equation (24). Figure 6.
Numerically calculated eigenvalues for several values ofeccentricity ( e = , . , . , . , . , . from top to bottomgoing light to dark) for the m = , k = − r -mode. Equation(24) is used to relate e and χ . The grey region reflects the valuesof q which are of interest for the study of thermonuclear burstoscillations as described in Section 4. the phenomenon of thermonuclear burst oscillations. Themost likely mode candidate to explain the observationalproperties of thermonuclear burst oscillations, suggested byHeyl (2004), is the m = r -mode with a small number of lat-itudinal nodes ( k = − , − ). The rotating frame frequency ofone of these modes must be small, ω / π ∼ − Hz (see Heyl2004, for further details), while the angular velocity of thesource is assumed to be large, between − Hz (Watts
MNRAS , 1–10 (2020)
B.F.A. van Baal et al.
Figure 7.
Same m = , k = − r -mode as Figure 6, using χ = e instead of Equation (24). Figure 8.
The eigenfunction of the m = , k = − retrograde r -mode with q = over the range µ = (the equator) to µ = . ( ≈ ◦ in latitude) for three different eccentricities. The relationbetween e and χ is given by Equation (24). For larger eccentrici-ties, the waves are more confined to the equator and the peak ofthe wave functions also moves towards lower latitudes. q falls in the range − .Sections 2.2 and 3.2 suggested a model for the varia-tion in gravitational acceleration with latitude, tailored toneutron stars based on calculations for a range of stellarmasses, radii, equations of state and rotation rates (AlGendy& Morsink 2014). We continue to use this model. Section 3showed that the asymptotic approximations deviate signifi-cantly from the numerical solutions for the r -modes, and assuch, for the rest of this analysis we rely purely on numericalsolutions.In Figure 6 the eigenvalues are shown for a range ofvalues of e , the region where we expect thermonuclear burstoscillation sources to be present is marked in grey. We wish to investigate the effect of a larger change in gravitationalacceleration across the surface of the neutron star. To dothis, we use the relation χ = e instead of χ = / e andshow results in Figure 7. This different relation means thatfor the same value of e , we now have a higher value of χ ,so the variation of the gravitation acceleration between theequator and the poles in greater than was shown in Figure6. A larger eccentricity increases the value of q for whichthe eigenvalues become positive. Using larger values of χ this effect becomes more extreme. Figure 6 shows that forthe relation χ = / e , this transition occurs at q = for e = . , and q = for e = . . In Figure 7, with therelation χ = e , it occurs at q = for e = . , and q = for e = . . In the spherical case, this transition occurs at q = (Lee & Saio 1997).The analytic approximation for the eigenvalue of r -modes, Equation (23), indicates that the r -modes in eccen-tric systems should reach the same asymptotic value as thosein spherical systems. The analytic approximations also pre-dict that the eigenvalue of an r -mode in an oblate system isalways smaller than that of a corresponding spherical sys-tem for the same value of q . This prediction, however, ismuch smaller than shown in the numerical results. The rateat which the eigenvalue converges to an asymptotic limit isdecreased significantly when including oblateness. For themost eccentric systems, the modes have not yet approachedtheir asymptotic limit by q = , while in spherical geome-try the asymptotic limit is reached at q ≈ , a much smallervalue. This discrepancy either indicates that the asymptoticlimit differs between spherical and oblate systems, or thatoblateness changes the rate at which the eigenvalue con-verges to an asymptotic limit more than is predicted by theanalytic approximations. Considering the very great discrep-ancy between the eigenvalues calculated using different re-lations between χ and e , it is possible that the approximateform for gravitational acceleration across the surface of thestar (the one used to make (11)) was not a reasonable choice,leading to a less accurate formula for the eigenvalues of r -modes.It can be seen that for the most eccentric system ( e = . ) the range of thermonuclear burst oscillation sources isclose to the point at which λ becomes positive. This ec-centricity, however, will likely only be reached by the mostrapidly rotating neutron stars at values of | q | (cid:29) . Since ω (and therefore | q | ) changes throughout the course of aburst, λ will also change during a burst, which was not thecase when using the spherical eigenvalue.Obtaining the correct eigenvalues is important for cal-culating the frequency of the wave when including the ra-dial structure of the mode. Equation (3) of Piro & Bildsten(2005) approximates the dependence of ω on the eigenvalueas λ / . By self-consistently solving ω = Ω / q = C λ ( q ) / ,where C is a constant, we can investigate the change in fre-quency of the wave when including oblateness. We considerthe eigenvalues at q = and for e = . , . as shownin Figure 6. We find that for e = . the wave frequencywould be reduced by − , while for e = . the wavefrequency would be reduced by − .Figure 8 shows eigenfunctions of the r -mode for severalvalues of eccentricity calculated using q = and Equation(24) to relate e and χ . For higher eccentricities, each com- MNRAS000
The eigenfunction of the m = , k = − retrograde r -mode with q = over the range µ = (the equator) to µ = . ( ≈ ◦ in latitude) for three different eccentricities. The relationbetween e and χ is given by Equation (24). For larger eccentrici-ties, the waves are more confined to the equator and the peak ofthe wave functions also moves towards lower latitudes. q falls in the range − .Sections 2.2 and 3.2 suggested a model for the varia-tion in gravitational acceleration with latitude, tailored toneutron stars based on calculations for a range of stellarmasses, radii, equations of state and rotation rates (AlGendy& Morsink 2014). We continue to use this model. Section 3showed that the asymptotic approximations deviate signifi-cantly from the numerical solutions for the r -modes, and assuch, for the rest of this analysis we rely purely on numericalsolutions.In Figure 6 the eigenvalues are shown for a range ofvalues of e , the region where we expect thermonuclear burstoscillation sources to be present is marked in grey. We wish to investigate the effect of a larger change in gravitationalacceleration across the surface of the neutron star. To dothis, we use the relation χ = e instead of χ = / e andshow results in Figure 7. This different relation means thatfor the same value of e , we now have a higher value of χ ,so the variation of the gravitation acceleration between theequator and the poles in greater than was shown in Figure6. A larger eccentricity increases the value of q for whichthe eigenvalues become positive. Using larger values of χ this effect becomes more extreme. Figure 6 shows that forthe relation χ = / e , this transition occurs at q = for e = . , and q = for e = . . In Figure 7, with therelation χ = e , it occurs at q = for e = . , and q = for e = . . In the spherical case, this transition occurs at q = (Lee & Saio 1997).The analytic approximation for the eigenvalue of r -modes, Equation (23), indicates that the r -modes in eccen-tric systems should reach the same asymptotic value as thosein spherical systems. The analytic approximations also pre-dict that the eigenvalue of an r -mode in an oblate system isalways smaller than that of a corresponding spherical sys-tem for the same value of q . This prediction, however, ismuch smaller than shown in the numerical results. The rateat which the eigenvalue converges to an asymptotic limit isdecreased significantly when including oblateness. For themost eccentric systems, the modes have not yet approachedtheir asymptotic limit by q = , while in spherical geome-try the asymptotic limit is reached at q ≈ , a much smallervalue. This discrepancy either indicates that the asymptoticlimit differs between spherical and oblate systems, or thatoblateness changes the rate at which the eigenvalue con-verges to an asymptotic limit more than is predicted by theanalytic approximations. Considering the very great discrep-ancy between the eigenvalues calculated using different re-lations between χ and e , it is possible that the approximateform for gravitational acceleration across the surface of thestar (the one used to make (11)) was not a reasonable choice,leading to a less accurate formula for the eigenvalues of r -modes.It can be seen that for the most eccentric system ( e = . ) the range of thermonuclear burst oscillation sources isclose to the point at which λ becomes positive. This ec-centricity, however, will likely only be reached by the mostrapidly rotating neutron stars at values of | q | (cid:29) . Since ω (and therefore | q | ) changes throughout the course of aburst, λ will also change during a burst, which was not thecase when using the spherical eigenvalue.Obtaining the correct eigenvalues is important for cal-culating the frequency of the wave when including the ra-dial structure of the mode. Equation (3) of Piro & Bildsten(2005) approximates the dependence of ω on the eigenvalueas λ / . By self-consistently solving ω = Ω / q = C λ ( q ) / ,where C is a constant, we can investigate the change in fre-quency of the wave when including oblateness. We considerthe eigenvalues at q = and for e = . , . as shownin Figure 6. We find that for e = . the wave frequencywould be reduced by − , while for e = . the wavefrequency would be reduced by − .Figure 8 shows eigenfunctions of the r -mode for severalvalues of eccentricity calculated using q = and Equation(24) to relate e and χ . For higher eccentricities, each com- MNRAS000 , 1–10 (2020) blate Thin Ocean Waves ponent of the eigenfunction is more confined to the equator.The area over which the P p component of the eigenfunc-tion is appreciably different from zero ( P p > . P p , max )changes by ≈ when the eccentricity is changed from to . . Using the relation χ = e , this area shrinks by .Greater equatorial confinement should result in a reducedpulsed amplitude and therefore the observed thermonuclearburst oscillation amplitudes since the contrasting pattern ofthe wave is contained in a smaller area Heyl (2005). In this work we investigated the impact of rotationally in-duced oblateness on waves that exist in a thin fluid layer onthe surface of a spheroid. Compared to their counterparts onthe surface of a slowly rotating sphere we found the eigen-functions and the eigenvalues change. Each family of modesis modified in a different manner. These results were ap-plied to the case of accreting neutron stars and specificallythe phenomenon of thermonuclear burst oscillations, whererotationally induced oblateness is expected to be relevantgiven the high spin rates. We investigate the m = r -modeand find that for more oblate systems, equatorial confine-ment increases and the wave frequency decreases.The numerical solutions and analytic approximationsfor m = − modes agree well with each other for the g -modes at | q | > and for the Yanai modes at | q | > . Nonew approximations are found for the Kelvin modes sincethe approximations from spherical geometry match well thesolutions calculated including oblateness. For the r -modeshowever, the analytic approximations no longer match thenumerical eigenvalues, especially for large oblateness andhigher variations in gravitational acceleration across the sur-face of the spheroid.We find that the eigenvalues of the g -modes, Yanaimodes and r -modes are altered by the presence of oblate-ness and do not depend particularly strongly on the modelfor gravitational acceleration across the surface of the star;however, they do depend on the degree of oblateness. This ef-fect is predicted by the asymptotic approximation Equation(22), where the leading term scales with e , the eccentricity,and not with χ , the parameter that characterises the changein gravitational acceleration across the surface of the star.The r -mode, on the other hand, depends strongly on χ ascan be seen from the differences in eigenvalues in Figures6 and 7. For the same values of e but higher χ , the eigen-values become smaller and only become positive towardshigher values of the spin parameter | q | . This behaviour isnot predicted from the asymptotic approximation, but Fig-ure 4 shows that the approximation no longer works whenmore eccentric systems are considered. For the Kelvin mode,the eigenvalues decrease slightly for larger oblateness but theeffect is small enough that the spherical approximations arestill valid. However, these modes do have one property thatdistinguishes them from the other families of modes; they be-come less equatorially confined for larger oblateness which isprecisely the opposite of what is seen for the g -mode, Yanaimode and r -mode (see Figures 5 and 8).For the m = r -mode we investigated how the eigenval-ues and eigenfunctions change due to oblateness. This hasimplications for the wave model for thermonuclear burst os- cillations. We find that for systems with an eccentricity be-tween . − . , the frequency of the wave drops by − ,while the area in which wave amplitude is significantly differ-ent from zero shrinks by ∼ . A less equatorially confinedpattern on the surface of the neutron star should lead to alarger burst oscillation amplitude; our findings suggest thatfor rapidly rotating (and thus more oblate) neutron stars,thermonuclear burst oscillations might (if caused by thismechanism) be more difficult to detect, as they would ex-hibit lower pulsed amplitudes. This might help explain whywe do not see burst oscillation sources which spin faster than Hz (if this is not due to the lack of more rapidly spinningneutron stars, Manchester et al. 2005). Note however thatthermonuclear burst oscillation amplitude also depends onother factors such as the accretion rate (Ootes et al. 2017),something that might confound efforts to isolate the effectsof rotation rate on amplitude in the current burst oscillationdata set (Bilous & Watts 2019).It is important to note that our results for the r -modes depend strongly on the parametrisation of the grav-itational acceleration across the neutron star surface. Withour simple parametrisation, Equation (7), significant differ-ences are present when increasing the difference in the grav-itational acceleration at the pole and the equator, changingthe asymptotic limit of the eigenvalues by up to a factor of2 and the area over which the eigenfunctions significantlydiffers from zero by as much as . Our parametrisationwas tailored for neutron stars specifically and based on thework of AlGendy & Morsink (2014), but other parametrisa-tions might fit better for applications other than for neutronstars.Inferring the rotation rate of neutron stars from thephenomenon of thermonuclear burst oscillations has impli-cations for continuous gravitational wave searches (Wattset al. 2008) and modelling the spacetime surrounding a neu-tron star (Riley et al. 2018). Any viable model for this phe-nomenon needs to provide a robust relationship between theburst oscillation frequency and the rotation rate of the star;in the mode model, this is the mode frequency in the ro-tating frame. Our study shows that the mode frequencycould decrease by as much as ∼ due to oblateness forthe most rapidly rotating systems. Maniopoulou & Ander-sson (2004) estimated that general relativistic effects candecrease the frequency of modes on the surface of a neutronstar by ∼ , and Chambers & Watts (2020) found similarreductions in frequency for r -modes when including the ra-dial structure. Other effects are also known to be important,such as nuclear burning throughout the course of the burst(Chambers et al. 2019), and perhaps magnetic fields (Heng& Spitkovsky 2009). A complete mode model for thermonu-clear burst oscillations should eventually include oblatenessalongside all of these other effects. ACKNOWLEDGEMENTS
The authors acknowledge support from ERC Starting GrantNo. 639217 CSINEUTRONSTAR (PI: Watts). This workbenefited from discussions at the ‘Bursting the Bubble’Lorentz Center workshop. The authors would like to ex-press their gratitude to the referee Richard Townsend forhis helpful and insightful comments, in particular with re-
MNRAS , 1–10 (2020) B.F.A. van Baal et al. gards to the impacts of the wave frequency, and in generalto improve the clarity of this work.
REFERENCES
AlGendy M., Morsink S. M., 2014, ApJ, 791, 78Arfken G. B., Weber H. J., 1999, Mathematical methods for physi-cistsBerkhout R. G., Levin Y., 2008, MNRAS, 385, 1029Bildsten L., Ushomirsky G., Cutler C., 1996, ApJ, 460, 827Bilous A. V., Watts A. L., 2019, ApJS, 245, 19Cavecchi Y., et al., 2011, ApJ, 740, L8Chambers F. R. N., Watts A. L., 2020, MNRAS, 491, 6032Chambers F. R. N., Watts A. L., Keek L., Cavecchi Y., GarciaF., 2019, ApJ, 871, 61Chandrasekhar S., 1969, Ellipsoidal figures of equilibrium. YaleUniv. PressEckart C., 1960, Hydrodynamics of Oceans and AtmospheresPergamonGalloway D. K., Muno M. P., Hartman J. M., Psaltis D.,Chakrabarty D., 2008, ApJS, 179, 360Gill A., 1982, Int. Geophys. Ser., 30, 662Heng K., Spitkovsky A., 2009, ApJ, 703, 1819Heyl J. S., 2004, ApJ, 600, 939Heyl J. S., 2005, MNRAS, 361, 504Hough S. S., 1898, Philosophical Transactions of the Royal Societyof London Series A, 191, 139Lee U., Saio H., 1986, MNRAS, 221, 365Lee U., Saio H., 1987, MNRAS, 224, 513Lee U., Saio H., 1997, ApJ, 491, 839Longuet-Higgins M. S., 1968, Philosophical Transactions of theRoyal Society of London Series A, 262, 511Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, AJ,129, 1993Maniopoulou A., Andersson N., 2004, MNRAS, 351, 1349Ootes L. S., Watts A. L., Galloway D. K., Wijnands R., 2017,ApJ, 834, 21Papaloizou J. C. B., Savonije G. J., 1997, MNRAS, 291, 651Pedlosky J., 1987, Geophysical fluid dynamics. Springer Science& Business MediaPedlosky J., 2003, Waves in the ocean and atmosphere: introduc-tion to wave dynamics. Springer Science & Business MediaPiro A. L., Bildsten L., 2005, ApJ, 629, 438Riley T. E., Raaijmakers G., Watts A. L., 2018, MNRAS, 478,1093Saio H., 1982, ApJ, 256, 717Staniforth A., White A., 2015, Quarterly Journal of the RoyalMeteorological Society, 141, 655Strohmayer T. E., Zhang W., Swank J. H., Smale A., TitarchukL., Day C., Lee U., 1996, ApJ, 469, L9Townsend R. H. D., 2003, MNRAS, 340, 1020Unno W., Osaki Y., Ando H., Saio H., Shibahashi H., 1989, Non-radial oscillations of stars, Tokyo: University of Tokyo Press,1989, 2nd ed.Watts A. L., 2012, ARA&A, 50, 609Watts A. L., Krishnan B., Bildsten L., Schutz B. F., 2008, MN-RAS, 389, 839White A. A., Wood N., 2012, Quarterly Journal of the RoyalMeteorological Society, 138, 980Yanai M., Maruyama T., 1966, Journal of the Meteorological So-ciety of Japan. Ser. II, 44, 291Yoshida K., 1960, Journal of the Oceanographical Society ofJapan, 15, 159This paper has been typeset from a TEX/L A TEX file prepared bythe author. MNRAS000