<figure><img src="image_1.png" /></figure>

HAL
open science

# The traditional approximation of rotation for rapidly rotating stars and planets

H. Dhouib, V. Prat, T. van Reeth, S. Mathis

► To cite this version:

H. Dhouib, V. Prat, T. van Reeth, S. Mathis. The traditional approximation of rotation for rapidly rotating stars and planets: I. The impact of strong deformation. Astronomy and Astrophysics - A&A, 2021, 652, pp.A154. 10.1051/0004-6361/202140615 . cea-03330529

HAL Id: cea-03330529

https://cea.hal.science/cea-03330529

Submitted on 1 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# The traditional approximation of rotation for rapidly rotating stars and planets

## I. The impact of strong deformation

H. Dhouib<fnref n="1" />, V. Prat<fnref n="1" />, T. Van Reeth<fnref n="2" />, and S. Mathis<fnref n="1" />

Received 19 February 2021 / Accepted 17 May 2021

## ABSTRACT

**Context.** The traditional approximation of rotation (TAR) is a treatment of the hydrodynamic equations of rotating and stably stratified fluids in which the action of the Coriolis acceleration along the direction of the entropy and chemical stratifications is neglected because it is weak in comparison with the buoyancy Archimedean force. This leads to the neglect of the horizontal projection of the rotation vector in the equations for the dynamics of gravito-inertial waves (GIWs). The dependent variables in those equations then become separable into radial and horizontal parts as in the non-rotating case. The TAR is built on the assumptions that the star is spherical (i.e., its centrifugal deformation is neglected) and uniformly rotating. However, it has recently been generalised to include the effects of a moderate centrifugal deformation using a perturbative approach.

**Aims.** We study the feasibility of carrying out a new generalisation to account for the centrifugal acceleration in the case of strongly deformed uniformly and rapidly rotating stars (and planets), and to identify the validity domain of this approximation.

**Methods.** We built a complete formalism analytically that allows the study of the dynamics of GIWs in spheroidal coordinates which take the flattening of uniformly and rapidly rotating stars into account by assuming the hierarchies of frequencies adopted within the TAR in the spherical case.

**Results.** Using 2D stellar models, we determine the validity domain of the generalised TAR as a function of the rotation rate of the star normalised by its critical angular velocity and its pseudo-radius. Assuming the anelastic and the two-dimensional Jeffreys-Wentzel-Kramers-Brillouin approximations, we derive a generalised Laplace tidal equation for the horizontal eigenfunctions of the GIWs and their asymptotic wave periods, which can be used to probe the structure and dynamics of rotating deformed stars with asteroseismology. The generalised TAR where the centrifugal deformation of a star (or planet) is taken into account non-perturbatively allows us to identify, within the framework of 2D Evolution STEllaire en Rotation models, the validity domain of this approximation which is reduced by increasing the rate of rotation. We can affirm with a level of confidence of 90% that the TAR remains applicable in all the space domain of deformed stars rotating at a rotation rate lower than 20% of the critical rotation rate.

**Conclusions.** A new generalisation of the TAR, which takes the centrifugal acceleration into account in a non-perturbative way, is derived. This generalisation allows us to study the detectability and the signature of the centrifugal effects on GIWs in rapidly rotating deformed stars (and planets). We found that the effects of the centrifugal acceleration in rapidly rotating early-type stars on GIWs are theoretically detectable in modern space photometry using observations from *Kepler*. We found also, by comparing the period spacing pattern computed with the standard and the generalised TAR, that the centrifugal acceleration affects the period spacing by increasing its values for low radial orders and by decreasing them slightly for high radial orders.

**Key words.** hydrodynamics – waves – stars: rotation – stars: oscillations – methods: analytical – methods: numerical

## 1. Introduction

The traditional approximation of rotation (TAR) was introduced for the first time in geophysics by Eckart (1960) to study the dynamics of the shallow Earth atmosphere and oceans. Then, Berthomieu et al. (1978) used it for the first time in an astrophysical context. Afterwards, Bildsten et al. (1996) and Lee & Saio (1997) applied it to stellar pulsations to study low-frequency nonradial oscillations in rotating stars, more precisely the gravito-inertial waves (GIWs). These waves (often called g modes) propagate under the combined action of the buoyancy force and the Coriolis acceleration. A subclass of these waves are the Rossby waves (often called r modes) whose dynamics are driven by the conservation of the vorticity of the fluid.

The TAR is applicable for low-frequency waves propagating in strongly stratified zones. In this case, a hierarchy of frequencies is verified. Namely, the buoyancy force is stronger than the Coriolis force (i.e., $2\Omega \ll N$, where $\Omega$ is the angular velocity and $N$ is the Brunt-Vaïsälä frequency) in the direction of stable entropy or chemical stratification and the Brunt-Vaïsälä frequency is larger than the frequency of the waves in the co-rotating frame $\omega$ (i.e., $\omega \ll N$). Considering these two assumptions, we can neglect the horizontal projection of the rotation vector $\Omega_H = \Omega \sin\theta$ in the linearised hydrodynamic equations, where $\theta$ is the colatitude. Then, the vertical component of the Coriolis acceleration can be neglected and the wave velocity is almost horizontal within the description of the dynamics of GIWs. In this framework, the linearised hydrodynamic

<footnote n="1">¹ Département d'Astrophysique-AIM, CEA/DRF/IRFU, CNRS/INSU, Université Paris-Saclay, Université Paris-Diderot, Université de Paris, 91191 Gif-sur-Yvette, France
e-mail: hachem.dhouib@cea.fr</footnote>

<footnote n="2">² Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, 3001 Leuven, Belgium</footnote>
system, which is coupled in the general case (Dintrans et al. 1999; Dintrans & Rieutord 2000), becomes separable as in the non-rotating case and can be rewritten in the form of the Laplace tidal equation (Laplace 1799).

When the GIWs are observed with high-precision photometry, they allow us to study the internal structure and dynamics of stars with asteroseismology (Aerts et al. 2010; Aerts 2021). GIWs also trigger transport of angular momentum and mixing inside the stars, modifying their evolution (e.g., Mathis et al. 2008; Mathis 2009; Rogers et al. 2013; Rogers 2015; Aerts et al. 2019; Neiner et al. 2020). In this context, GIWs which propagate in stably stratified stellar radiation regions are of major interest. First, they are detected and analysed with asteroseismology in intermediate- and high-mass stars (e.g., Neiner et al. 2012; Van Reeth et al. 2015) allowing the characterisation of the internal stellar rotation (Van Reeth et al. 2016, 2018). In addition, GIWs are one of the candidates to explain the small rotation contrast observed between the stellar surface and core (Mathis 2009; Rogers 2015; Aerts et al. 2019). So, the TAR and its flexibility allow us to derive powerful seismic diagnostics, based on the period spacings between consecutive radial order *g* modes in uniformly and differentially rotating spherical stars (Bouabid et al. 2013; Ouazzani et al. 2017; Van Reeth et al. 2018). This enables us to determine the properties of the chemical stratification and the rotation rate near the convective core of early-type stars, which constrain our modelling of stellar structure, evolution and internal angular momentum transport (e.g., Pedersen et al. 2018; Aerts et al. 2018, 2019; Ouazzani et al. 2019).

In addition to the two conditions described above, three other assumptions are made when applying the TAR. First, the rotation is assumed to be uniform. Second, the star is assumed to be spherical, in other words, the centrifugal deformation of the star is neglected (i.e., $\Omega \ll \Omega_K \equiv \sqrt{GM/R^3}$, where $\Omega_K$ is the Keplerian critical (breakup) angular velocity, and *G*, *M*, and *R* are the universal constant of gravity, the mass of the star, and the stellar radius, respectively). Third, the Cowling approximation (Cowling 1941) is assumed where the fluctuations of the gravitational potential caused by the redistribution of mass by waves motion is neglected.

In addition, the anelastic approximation (Spiegel & Veronis 1960) is assumed where the high-frequency acoustic modes can be filtered out. Finally, for asymptotic low frequency modes which are rapidly oscillating in the vertical direction, we can assume the Jeffreys, Wentzel, Kramers, Brillouin (JWKB) approximation (e.g., Fröman & Fröman 1965).

The assumption of uniform rotation was abandoned by Ogilvie & Lin (2004) and Mathis (2009), who took into account the effects of general differential rotation (both in radius and latitude in a spherical star) on low-frequency GIWs within the framework of the TAR. The results of this new formalism have been successfully applied in Van Reeth et al. (2018) to derive the variation of the asymptotic period spacing in the case of a weak radial differential rotation and evaluate the sensitivity of GIWs to the effect of differential rotation.

The assumption of spherical symmetry was partially abandoned by Mathis & Prat (2019), who generalised the TAR for slightly deformed moderately rotating stars by considering the effects of the centrifugal acceleration using a perturbative approach. First, they derived the centrifugal perturbation of the spherical hydrostatic balance by determining the deformation of an isobar in the case of moderate uniform rotation. Afterwards, using this deformation they defined a mapping between the

spherical coordinates system and a spheroidal coordinates system which takes into account this weak centrifugal perturbation of the star. Then, by keeping only first-order terms in the deformation, they re-derived the linearised hydrodynamic equations (Saio 1981; Lee & Baraffe 1995) and deduced the perturbed Laplace tidal equation within this new coordinates system. Then, they carried out a numerical exploration of the eigenvalues and horizontal eigenfunctions (the so-called Hough functions; Hough 1898) of the perturbed Laplace tidal equation. Finally, they derived the asymptotic expression for the frequencies of GIWs, including the effect of the centrifugal acceleration. This perturbative framework was then optimised for practical intense seismic forward modelling (Aerts 2021) using one-dimensional stellar structure models by Henneco et al. (2021).

In this work, we consider the general case of rapidly and uniformly rotating stars and planets which will generalise for the first time the previous work of Lee & Saio (1997) and Mathis & Prat (2019). The shape of such objects becomes a strongly deformed spheroid because of the action of the centrifugal acceleration which cannot be modelled as a perturbation in this case. We study if the TAR could be generalised to include the effects of this deformation for any stellar (planetary) rotation including those close to their break-up velocity. The derived formalism could have several key applications. For instance, it could be used to build new seismic diagnostic and to study angular momentum transport (e.g., Lee & Saio 1993; Mathis et al. 2008; Mathis 2009; Lee et al. 2014) and tidal dissipation (e.g., Ogilvie & Lin 2004, 2007; Braviner & Ogilvie 2014) induced by low-frequency GIWs in deformed stars and planets.

To do so, we first introduce in Sect. 2 a new spheroidal coordinate system which follows the shape of a deformed star for any rotation. It is based on the work by Bonazzola et al. (1998) and was first used by Rieutord et al. (2005) to develop rotating polytropic models. Then Lignières et al. (2006) employed it to study acoustic oscillations in strongly deformed rotating stars. Afterwards, Reese et al. (2006, 2009) and Ouazzani et al. (2012) make use of it to develop the two-dimensional oscillation program (TOP) and the adiabatic code of oscillation including rotation (ACOR), respectively. They compute non-radial pulsations of rotating stars. Finally, Espinosa Lara & Rieutord (2013) utilised it to develop the Evolution STEllaire en Rotation (ESTER) code which computes the structure of an early-type star including its differential rotation and the associated meridional circulation. Here, we first consider the case of uniform rotation to distinguish the effects of the deformation from those of differential rotation associated to the centrifugal and Coriolis accelerations, respectively. From there, we derive the system of linearised hydrodynamic equations in spheroidal coordinates. In Sect. 3, we apply the TAR in this new coordinate system by adopting the adequate assumptions. In Sect. 4, we can then rewrite the oscillation equations in the form of a new generalised Laplace tidal equation (GLTE) and deduce the asymptotic frequencies of low-frequency GIWs. In Sect. 5, we identify the validity domain of the TAR when using ESTER models and carry out a numerical exploration of the eigenvalues and the eigenfunctions of the GLTE. In Sect. 6, we quantify the centrifugal deformation effect of the star on the asymptotic period spacing pattern and we discuss its detectability. In Sect. 7, we verify the validity of the hierarchy of the terms in the linearised hydrodynamics equations imposed by the TAR in deformed stars. Finally, we discuss and summarise our work and results in Sect. 8.
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Spheroidal coordinate system with $n = 13$ subdomains used to compute the equilibrium model of a star rotating at 60% of its break-up velocity corresponding to the case of the $3M_\odot$ ESTER model (with the central fraction in Hydrogen $X_\zeta = 0.7$) studied in Sect. 5. $\zeta$ and $\theta$ are the pseudo-radius and the colatitude, respectively. $\eta_i$ and $\Delta R_i(\theta)$ are defined in our mapping (Eq. (1)).</figcaption></figure>

## 2. Hydrodynamic equations in uniformly rotating deformed stars and planets

### 2.1. Spheroidal geometry

Because of the flattened shape of a rapidly rotating star, we use the spheroidal coordinate system ($\zeta, \theta, \varphi$) proposed by Bonazzola et al. (1998), where $\zeta$ is the radial coordinate (the pseudo-radius), $\theta$ the colatitude and $\varphi$ the azimuth. It was designed not only to fit the shape of the star at the surface but also to get closer to the spherical shape at the centre as illustrated in Fig. 1. Following Rieutord et al. (2013), this new coordinate system can be linked to the usual spherical one ($r, \theta, \varphi$) via the following mapping

$$r(\zeta, \theta) = a_i \xi \Delta\eta_i + R_i(\theta) + A_i(\xi) (\Delta R_i(\theta) - a_i \Delta\eta_i), \text{ for } \zeta \in [\eta_i, \eta_{i+1}], \quad (1)$$

where we split the spheroidal domain $\mathcal{D}$ in $n$ subdomains $\mathcal{D}_{i \in [0, n-1]} \in [R_i(\theta), R_{i+1}(\theta)]$ where $R_{i \in [0, n]}(\theta)$ are series of functions, such that $R_n(\theta) = R_s(\theta)$ is the outer boundary and $R_0(\theta) = 0$ is the centre. Additionally, $\eta_i = R_i(\theta = 0)$ are the polar radii of the interfaces between the subdomains, $\Delta\eta_i = \eta_{i+1} - \eta_i$, $\Delta R_i(\theta) = R_{i+1}(\theta) - R_i(\theta)$ and $\xi = (\zeta - \eta_i)/\Delta\eta_i$. The functions $A_{i \in [1, n-1]}(\xi) = -2\xi^3 + 3\xi^2$ and $A_0(\xi) = -1.5\xi^5 + 2.5\xi^3$ and the constants $a_i$ (in practice $a_i = 1$) are chosen to satisfy the boundary conditions between the different subdomains so that the mapping is continuous and derivable. The other coordinates ($\theta$ and $\varphi$) remain unchanged. Since we choose the polar radius as the unit length $R_s(\theta = 0) = 1$, $\zeta = 1$ coincides with the surface (i.e., $r(\zeta = 1, \theta) = R_s(\theta)$), while $r = 0$ is obtained at the centre when $\zeta = 0$ (i.e., $r(\zeta = 0, \theta) = 0$). Along the rotation axis ($\theta = 0$), the deformed radial coordinate $\zeta$ is the distance from the centre normalised by the polar radius and when the centrifugal deformation is zero, $\zeta$ is equivalent to the radius, hence the name pseudo-radius.

Once the spheroidal coordinate system is established, it is mandatory to specify the corresponding basis. To do so, we

define the spheroidal basis ($a_\zeta, a_\theta, a_\varphi$) from the natural covariant basis ($E_\zeta, E_\theta, E_\varphi$) defined as $E_i = \partial_i(\mathbf{r}e_r)$

$$\left\{
\begin{aligned}
a_{\zeta} &= \frac{\zeta^2}{r^2 r_{\zeta}} \mathbf{E}_{\zeta} = \frac{\zeta^2}{r^2} \mathbf{e}_r \\
a_{\theta} &= \frac{\zeta}{r^2 r_{\zeta}} \mathbf{E}_{\theta} = \frac{\zeta}{r^2 r_{\zeta}} (r_{\theta} \mathbf{e}_r + r \mathbf{e}_{\theta}), \\
a_{\varphi} &= \frac{\zeta}{r^2 r_{\zeta} \sin \theta} \mathbf{E}_{\varphi} = \frac{\zeta}{r r_{\zeta}} \mathbf{e}_{\varphi}
\end{aligned}
\right.
\quad (2)$$

where $r_\zeta \equiv \partial_\zeta r$, $r_\theta \equiv \partial_\theta r$ and ($e_r, e_\theta, e_\varphi$) is the usual spherical unit-vector basis. This basis has the advantage of being reduced to the usual spherical unit-vector basis when we approach the spherical geometry ($R_s(\theta)$ tends to the constant function 1).

### 2.2. Linearised hydrodynamic equations in spheroidal coordinates

Once the geometry is fixed, we can derive the complete adiabatic inviscid system of linearised equations to treat the wave dynamics in a uniformly rotating, strongly deformed star or planet. First, the linearised momentum equation for an inviscid fluid is written as (we refer the reader to Appendix A.1 for a detailed derivation)

$$\begin{aligned}
(\partial_t + \Omega \partial_\varphi) \left[ \frac{\zeta^2 r_\zeta}{r^2} v^\zeta + \frac{\zeta r_\theta}{r^2} v^\theta \right] \\
= 2\Omega \frac{\zeta \sin \theta}{r} v^\varphi - \frac{1}{\rho_0} \partial_\zeta P' + \frac{\rho'}{\rho_0^2} \partial_\zeta P_0 - \partial_\zeta \Phi',
\end{aligned}
\quad (3)$$

$$\begin{aligned}
(\partial_t + \Omega \partial_\varphi) \left[ \frac{\zeta^2 r_\theta}{r^2} v^\zeta + \frac{\zeta (r^2 + r_\theta^2)}{r^2 r_\zeta} v^\theta \right] \\
= 2\Omega \frac{\zeta (r_\theta \sin \theta + r \cos \theta)}{r r_\zeta} v^\varphi - \frac{1}{\rho_0} \partial_\theta P' + \frac{\rho'}{\rho_0^2} \partial_\theta P_0 - \partial_\theta \Phi',
\end{aligned}
\quad (4)$$

$$\begin{aligned}
\frac{\zeta}{r_\zeta} (\partial_t + \Omega \partial_\varphi) v^\varphi = -2\Omega \frac{\zeta^2 \sin \theta}{r} v^\zeta \\
-2\Omega \frac{\zeta (r_\theta \sin \theta + r \cos \theta)}{r r_\zeta} v^\theta - \frac{1}{\rho_0 \sin \theta} \partial_\varphi P' - \frac{1}{\sin \theta} \partial_\varphi \Phi',
\end{aligned}
\quad (5)$$

where $v^\zeta$, $v^\theta$ and $v^\varphi$ are the contravariant components of the velocity field $v = v^i a_i$ and $\Omega = \Omega(\cos \theta e_r - \sin \theta e_\theta)$ is the uniform angular velocity. $\rho$, $\Phi$ and $P$ are the fluid density, gravitational potential and pressure, respectively. Each of these scalar quantities has been expanded as

$$X(\zeta, \theta, \varphi, t) = X_0(\zeta, \theta) + X'(\zeta, \theta, \varphi, t),$$

where $X_0$ is the hydrostatic component of $X$ and $X'$ the wave's associated linear fluctuation. When assuming the Cowling (1941) approximation, the fluctuation of the gravitational potential is neglected ($\Phi' = 0$).

Subsequently, the linearised continuity equation is obtained (we refer the reader to Appendix A.2 for a detailed derivation)

$$\begin{aligned}
(\partial_t + \Omega \partial_\varphi) \rho' = - \frac{\zeta^2 \partial_\zeta \rho_0}{r^2 r_\zeta} v^\zeta - \frac{\zeta \partial_\theta \rho_0}{r^2 r_\zeta} v^\theta \\
- \frac{\zeta^2 \rho_0}{r^2 r_\zeta} \left[ \frac{\partial_\zeta (\zeta^2 v^\zeta)}{\zeta^2} + \frac{\partial_\theta (\sin \theta v^\theta)}{\zeta \sin \theta} + \frac{\partial_\varphi v^\varphi}{\zeta \sin \theta} \right].
\end{aligned}
\quad (6)$$
Then, the linearised energy equation in the adiabatic limit is derived (we focus in this work on adiabatic oscillations as in studies of the TAR in spherical stars)

$$ (\partial_t + \Omega \partial_\varphi) \left( \frac{1}{\Gamma_1} \frac{P'}{P_0} - \frac{\rho'}{\rho_0} \right) = \frac{N^2}{\|g_{\text{eff}}\|^2} \mathbf{v} \cdot \mathbf{g}_{\text{eff}}, \quad (7) $$

where $\Gamma_1 = (\partial \ln P_0 / \partial \ln \rho_0)_S$ ($S$ being the macroscopic entropy) is the adiabatic exponent, $g_{\text{eff}} = -\nabla\Phi_0 + \frac{1}{2}\Omega^2\nabla(r^2\sin^2\theta) = \nabla P_0 / \rho_0$ is the background effective gravity which includes the centrifugal acceleration and $N^2$ is the squared Brunt-Väisälä frequency given by

$$ N^2(\zeta, \theta) = g_{\text{eff}} \cdot \left( -\frac{1}{\Gamma_1} \frac{\nabla P_0}{P_0} + \frac{\nabla \rho_0}{\rho_0} \right). \quad (8) $$

Finally, since in this work we are interested in studying standing eigenmodes of oscillation we consider a discrete spectrum of eigenfrequencies rather than an integral over a continuous spectrum. So, we expand the wave's velocity and scalar quantities fluctuations ($X' = \{\rho', P'\}$) on discrete Fourier series both in time and in azimuth

$$ \mathbf{v}(\zeta, \theta, \varphi, t) \equiv \sum_{\omega^{\text{in},m}} \{\mathbf{u}(\zeta, \theta) \exp[i(\omega^{\text{in}}t - m\varphi)]\}, \quad (9) $$

$$ X'(\zeta, \theta, \varphi, t) \equiv \sum_{\omega^{\text{in},m}} \{\tilde{X}(\zeta, \theta) \exp[i(\omega^{\text{in}}t - m\varphi)]\}, \quad (10) $$

where $m$ is the azimuthal order ($m$ is an integer because we are looking for a linear solution of a wave problem whose axisymmetric background depends only on $\zeta$ and $\theta$) and $\omega^{\text{in}}$ is the wave eigenfrequency in an inertial reference frame (we don't put an index for the eigenfrequencies in order to make the notations simpler). In a uniformly rotating region, the waves are Doppler shifted due to the rotation. We can define the wave eigenfrequency $\omega$ in the co-rotating reference frame as

$$ \omega = \omega^{\text{in}} - m\Omega. \quad (11) $$

## 3. Generalised TAR

### 3.1. Approximation on the stratification profile: $N^2(\zeta, \theta) \approx N^2(\zeta)$

To obtain a separable system of equations by applying the TAR, a first necessary condition is that the Brunt-Väisälä frequency $N^2$ depends mainly on the pseudo-radius $\zeta$, that is to say it has a small variation with respect to the colatitude $\theta$. This also implies that the hydrostatic pressure $P_0$ and the hydrostatic density $\rho_0$ depend mainly on $\zeta$ so the linearised energy equation simplifies into

$$ \frac{1}{\Gamma_1} \frac{\tilde{P}}{P_0} - \frac{\tilde{\rho}}{\rho_0} = \frac{1}{i\omega g_{\text{eff}}} \frac{N^2}{r^2 r_\zeta} u^\zeta. \quad (12) $$

In order to compute the relative error made by adopting this approximation and to define its validity domain, we compare in Sect. 5.2 an approximate value of $N^2$ where it depends only on $\zeta$ to the exact one depending on $\zeta$ and $\theta$ by fixing the value of $\theta$ and varying the angular velocity of the star using two-dimensional ESTER stellar models.

### 3.2. The TAR with centrifugal acceleration

Using this first approximation, and assuming the same hierarchies of frequencies ($2\Omega \ll N$ and $\omega \ll N$) and velocity scales ($|v^\zeta| \ll |v^\theta|, |v^\varphi|$)) that in the spherical and in the weakly deformed cases (this will be discussed in details in Sect. 7), we now built the generalised framework for the TAR in the case of a uniformly and rapidly rotating strongly deformed star (planet).

By adopting the approximation $N^2(\zeta, \theta) \approx N^2(\zeta)$ and the Cowling approximation, we can rewrite the radial momentum Eq. (3) as

$$ i\omega \left[ \frac{\zeta^2 r_\zeta}{r^2} u^\zeta + \frac{\zeta r_\theta}{r^2} u^\theta \right] = 2\Omega \frac{\zeta \sin \theta}{r} u^\varphi - \partial_\zeta \tilde{W} - \frac{\tilde{P}}{\rho_0^2} \partial_\zeta \rho_0 - \frac{N^2}{i\omega} \frac{\zeta^2}{r^2 r_\zeta} u^\zeta, \quad (13) $$

where $\tilde{W} = \tilde{P}/\rho_0$ is the normalised pressure. The propagation of low-frequency GWs can be studied within the anelastic approximation in which acoustic waves are filtered out (Dintrans & Rieutord 2000; Mathis 2009). Therefore, Eqs. (13) and (12) can be simplified accordingly by neglecting the terms $(\tilde{P}/\rho_0^2)\partial_\zeta\rho_0$ and $(1/\Gamma_1)\tilde{P}/P_0$. Within the TAR, we focus on low-frequency waves (i.e., $\omega \ll N$) propagating in strongly stratified regions (i.e., $2\Omega \ll N$). In this case, we can neglect the vertical component of the Coriolis acceleration because it is dominated by the buoyancy force and the vertical wave velocity since $|v^\zeta| \ll |v^\theta|, |v^\varphi|$ because of the strong stable stratification. Thus, the radial momentum equation can be simplified into

$$ i\omega \partial_\zeta \tilde{W} + N^2 \zeta^2 \mathcal{A} u^\zeta = 0, \quad (14) $$

where

$$ \mathcal{A}(\zeta, \theta) = \frac{1}{r^2 r_\zeta}. \quad (15) $$

Subsequently, we examine the latitudinal component of the momentum equation (Eq. (4)), which reduces, for the same reasons, to

$$ i\omega\zeta\mathcal{B}u^\theta - 2\Omega\zeta C u^\varphi + \partial_\theta \tilde{W} = 0, \quad (16) $$

where

$$ \mathcal{B}(\zeta, \theta) = \frac{r^2 + r_\theta^2}{r^2 r_\zeta}, \quad (17) $$

$$ C(\zeta, \theta) = \frac{r_\theta \sin \theta + r \cos \theta}{r r_\zeta}. \quad (18) $$

Finally, the azimuthal component of the momentum equation (Eq. (5)) simplifies onto

$$ i\omega\zeta\mathcal{D}u^\varphi + 2\Omega\zeta C u^\theta - \frac{im}{\sin\theta}\tilde{W} = 0, \quad (19) $$

where

$$ \mathcal{D}(\zeta, \theta) = \frac{1}{r_\zeta}. \quad (20) $$

We thus obtain a system of equations for strongly deformed rotating stars which has the same mathematical form that in the case of spherically symmetric stars (Lee & Saio 1997) and weakly deformed stars (Mathis & Prat 2019), and where we thus still manage to partially decouple the vertical and horizontal
components of the velocity. By solving the system formed by Eqs. (16) and (19) we can express $u^\theta$ and $u^\varphi$ as a function of $\tilde{W}$ as follows

$$u^{\theta}(\zeta, \theta) = i \frac{1}{\omega \zeta} \frac{1}{\mathcal{B}(\zeta, \theta)} \left[ \left( 1 + \nu^2 \frac{C^2(\zeta, \theta)}{\mathcal{E}(\zeta, \theta)\mathcal{B}(\zeta, \theta)} \right) \partial_{\theta} \tilde{W} - \frac{mv}{\sin \theta} \frac{C(\zeta, \theta)}{\mathcal{E}(\zeta, \theta)} \tilde{W} \right], \quad (21)$$

$$u^{\varphi}(\zeta, \theta) = -\frac{1}{\omega \zeta} \frac{1}{\mathcal{E}(\zeta, \theta)} \left[ \nu \frac{C(\zeta, \theta)}{\mathcal{B}(\zeta, \theta)} \partial_{\theta} \tilde{W} - \frac{m}{\sin \theta} \tilde{W} \right], \quad (22)$$

where we introduce

$$\mathcal{E}(\zeta, \theta) = \mathcal{D}(\zeta, \theta) - \nu^2 \frac{C^2(\zeta, \theta)}{\mathcal{B}(\zeta, \theta)}, \quad (23)$$

and the spin parameter

$$\nu = \frac{2\Omega}{\omega}. \quad (24)$$

The structure of the new equations that we obtain in the spheroidal case is similar to the one in the usual spherical case (e.g., Lee & Saio 1997; Townsend 2003; Pantillon et al. 2007). The most important difference is that the coefficients $\mathcal{A}, \mathcal{B}, \mathcal{C}, \mathcal{D}$ and $\mathcal{E}$ are function of $\zeta$ and $\theta$ through $r(\zeta, \theta)$ and their derivatives. As this is done in Lee & Saio (1997) and Mathis & Prat (2019), we introduce the reduced latitudinal coordinate $x = \cos \theta$, so that Eqs. (21) and (22) transform to

$$\begin{aligned} u^{\theta}(\zeta, x) &= \mathcal{L}_{\text{vm}}^{\theta} \left[ \tilde{W}(\zeta, x) \right] \\ &= -i \frac{1}{\omega \zeta} \frac{1}{\mathcal{B}(\zeta, x)} \frac{1}{\sqrt{1-x^2}} \left[ \left(1-x^2\right) \left(1+\nu^2 \frac{C^2(\zeta, x)}{\mathcal{E}(\zeta, x)\mathcal{B}(\zeta, x)}\right) \partial_x + m\nu \frac{C(\zeta, x)}{\mathcal{E}(\zeta, x)} \right] \tilde{W}, \quad (25) \end{aligned}$$

$$\begin{aligned} u^{\varphi}(\zeta, x) &= \mathcal{L}_{\text{vm}}^{\varphi} \left[ \tilde{W}(\zeta, x) \right] \\ &= \frac{1}{\omega \zeta} \frac{1}{\mathcal{E}(\zeta, x)} \frac{1}{\sqrt{1-x^2}} \left[ \left(1-x^2\right) \nu \frac{C(\zeta, x)}{\mathcal{B}(\zeta, x)} \partial_x + m \right] \tilde{W}. \quad (26) \end{aligned}$$

## 4. Dynamics of low-frequency gravito-inertial waves

Our goal in this section is to derive the generalised Laplace tidal equation (GLTE) for the normalised pressure $\tilde{W}$, which allows us to compute the frequencies and periods of low-frequency GIWs and to build the corresponding seismic diagnostics following the method by Bouabid et al. (2013), Van Reeth et al. (2018), and Mathis & Prat (2019). Applying the anelastic approximation and the approximation $N^2(\zeta, \theta) \approx N^2(\zeta)$ in the continuity equation (Eq. (6)) allows us to neglect the terms $\partial_t \rho'$ and $\partial_\theta \rho_0$ respectively, so it simplifies into

$$\zeta \partial_{\zeta} \rho_0 u^{\zeta} + \rho_0 \left[ \frac{\partial_{\zeta} (\zeta^2 u^{\zeta})}{\zeta} + \frac{1}{\sin \theta} \partial_{\theta} (\sin \theta u^{\theta}) - \frac{im}{\sin \theta} u^{\varphi} \right] = 0. \quad (27)$$

### 4.1. JWKB approximation

Assuming that $\omega \ll N$, each scalar field and each component of $u$ can be expanded using the two-dimensional Jeffreys-Wentzel-Kramers-Brillouin (JWKB) approximation which was introduced for strongly stratified differentially rotating fluids by Mathis (2009). In this case, the vertical wave number $k_V \equiv k_\zeta$ is assumed to be very large, and the associated wavelength is thus very small. Therefore, the spatial variation of the wave in the pseudo-radial direction is very fast compared to that of the hydrostatic background given by $\rho_0, P_0$ and $g_{\text{eff}}$. So the spatial structure of the waves can be described by the product of a rapidly oscillating plane-like wave function in the pseudo-radial direction multiplied by a bi-dimensional slowly varying envelope:

$$\tilde{W}(\zeta, \theta) = \sum_k \left\{ w_{vkm}(\zeta, \theta) \frac{A_{vkm}}{k_{V;vkm}^{1/2}} \exp \left[ i \int^{\zeta} k_{V;vkm} d\zeta \right] \right\}, \quad (28)$$

$$u^j(\zeta, \theta) = \sum_k \left\{ \hat{u}_{vkm}^j(\zeta, \theta) \frac{A_{vkm}}{k_{V;vkm}^{1/2}} \exp \left[ i \int^{\zeta} k_{V;vkm} d\zeta \right] \right\}, \quad (29)$$

with $j \equiv \{\zeta, \theta, \varphi\}$, $k$ is the index of a latitudinal eigenmode (cf. Sect. 4.3) and $A_{vkm}$ is the amplitude of the wave. Substituting the expansion given in Eqs. (28) and (29) into Eqs. (14), (25) and (26), the final pseudo-radial, latitudinal and azimuthal components of the velocity are obtained:

$$\hat{u}_{vkm}^{\zeta}(\zeta, x) = \frac{k_{V;vkm}(\zeta)}{N^2(\zeta)} \frac{\omega}{\zeta^2 \mathcal{A}(\zeta, x)} w_{vkm}(\zeta, x), \quad (30)$$

$$\hat{u}_{vkm}^{\theta}(\zeta, x) = \mathcal{L}_{\text{vm}}^{\theta} [w_{vkm}(\zeta, x)], \quad (31)$$

$$\hat{u}_{vkm}^{\varphi}(\zeta, x) = \mathcal{L}_{\text{vm}}^{\varphi} [w_{vkm}(\zeta, x)]. \quad (32)$$

### 4.2. Approximation: $\mathcal{A}(\zeta, \theta) \approx \mathcal{A}(\zeta)$

To be able to derive a generalised Laplace tidal equation as in the case of slightly deformed spheres (Mathis & Prat 2019) and to introduce its horizontal eigenvalues $\Lambda_{vkm}$, which should depend only on $\zeta$, we have to do a partial separation between the pseudo-radial and latitudinal variables in the pseudo-radial component of the velocity $\hat{u}^\zeta$. In other words, we have to write Eq. (30) as the product of a function which depends only on $\zeta$ and the normalised pressure $w_{vkm}$. To do so, we have to assume that the coefficient $\mathcal{A}(\zeta, \theta)$ depends mainly on $\zeta$. To evaluate the relative error made by adopting this approximation and to define its validity domain, we compare in Sect. 5.2 an approximate value of $\mathcal{A}$ where it depends only on $\zeta$ to the exact one depending on $\zeta$ and $\theta$ by fixing the value of $\theta$ and varying the angular velocity of the star using two-dimensional ESTER stellar models.

### 4.3. Generalised Laplace tidal equation (GLTE)

Substituting the radial component of the velocity (Eq. (30)) into the continuity equation (Eq. (27)) we get:

$$\frac{\omega_{km}^2 k_{V;vkm}^2}{\mathcal{A}N^2} w_{vkm} + \omega_{km}\zeta \left( i \partial_x \left( \sqrt{1-x^2}\hat{u}^\theta \right) - \frac{m}{\sqrt{1-x^2}} \hat{u}^\varphi \right) = 0, \quad (33)$$

where the JWKB approximation allows us to neglect $\zeta \partial_\zeta \rho_0 / \rho_0 u^\zeta$ and $2u^\zeta$ in front of the dominant term $\zeta \partial_\zeta u^\zeta$. Using Eqs. (31) and
<table><caption>**Table 1.** Terms involved in the derivation of the GLTE in the general case of spheroidal geometry and in the particular case of spherical geometry (we refer the reader to Appendix B for a visual presentation of these terms).</caption><thead><tr><th>Terms</th><th>Spheroidal geometry</th><th>Spherical geometry</th></tr></thead><tbody><tr><td>A</td><td><sup>1</sup>&frasl;<sub>r<sup>2</sup>r<sub>ζ</sub></sub></td><td><sup>1</sup>&frasl;<sub>ζ<sup>2</sup></sub></td></tr><tr><td>B</td><td><sup>r<sup>2</sup> + (1 - x<sup>2</sup>)r<sub>x</sub><sup>2</sup></sup>&frasl;<sub>r<sup>2</sup>r<sub>ζ</sub></sub></td><td>1</td></tr><tr><td>C</td><td><sup>-(1 - x<sup>2</sup>)r<sub>x</sub> + rx</sup>&frasl;<sub>r<sub>ζ</sub></sub></td><td>x</td></tr><tr><td>D</td><td><sup>rr<sub>ζ</sub></sup>&frasl;<sub>r<sub>ζ</sub></sub></td><td>1</td></tr><tr><td>E</td><td><sup>D - v<sup>2</sup>C<sup>2</sup></sup>&frasl;<sub>B</sub> = <sup>1</sup>&frasl;<sub>r<sub>ζ</sub></sub> - v<sup>2</sup><sup>-(1 - x<sup>2</sup>)r<sub>x</sub> + rx</sup>&frasl;<sub>r<sup>2</sup> + (1 - x<sup>2</sup>)r<sub>x</sub><sup>2</sup></sub></td><td>1 - v<sup>2</sup>x<sup>2</sup></td></tr><tr><td>∂<sub>x</sub>B</td><td><sup>2rr<sub>x</sub> - 2xr<sub>x</sub><sup>2</sup> + 2(1 - x<sup>2</sup>)r<sub>x</sub>r<sub>xx</sub> - B(2rr<sub>x</sub>r<sub>ζ</sub> + r<sup>2</sup>r<sub>ζx</sub>)</sup>&frasl;<sub>r<sup>2</sup>r<sub>ζ</sub></sub></td><td>0</td></tr><tr><td>∂<sub>x</sub>C</td><td><sup>2xr<sub>x</sub> - (1 - x<sup>2</sup>)r<sub>xx</sub> + r + xr<sub>x</sub> - C(r<sub>x</sub>r<sub>ζ</sub> + rr<sub>ζx</sub>)</sup>&frasl;<sub>rr<sub>ζ</sub></sub></td><td>1</td></tr><tr><td>∂<sub>x</sub>D</td><td><sup>-D<sup>2</sup>r<sub>ζx</sub></sup>&frasl;<sub>D<sup>2</sup>r<sub>ζ</sub></sub></td><td>0</td></tr><tr><td>∂<sub>x</sub>(C/EB)</td><td><sup>∂<sub>x</sub>C(DB - v<sup>2</sup>r<sub>ζ</sub>C<sup>2</sup>) - C(B∂<sub>x</sub>D + D∂<sub>x</sub>B - v<sup>2</sup>r<sub>ζx</sub>C<sup>2</sup> - 2v<sup>2</sup>r<sub>ζ</sub>C∂<sub>x</sub>C)</sup>&frasl;<sub>(DB - v<sup>2</sup>r<sub>ζ</sub>C<sup>2</sup>)<sup>2</sup></sub></td><td><sup>1 + v<sup>2</sup>x<sup>2</sup></sup>&frasl;<sub>(1 - v<sup>2</sup>x<sup>2</sup>)<sup>2</sup></sub></td></tr></tbody></table>

Notes. $r_{\zeta x}$ and $r_{xx}$ denote $\partial_{\zeta x}^2 r$ and $\partial_{xx}^2 r$, respectively and $f$ is a function of $x$ given by the stellar model.

(32), we obtain the GLTE for the normalised pressure $w_{vkm}$

$$ \mathcal{L}_{vm}[w_{vkm}] = \partial_x \left[ \frac{1}{\mathcal{B}} \left( 1 + \frac{\nu^2 C^2}{\mathcal{E}\mathcal{B}} \right) (1 - x^2) \partial_x w_{vkm} \right] \\ + \left[ m\nu \partial_x \left( \frac{C}{\mathcal{E}\mathcal{B}} \right) - \frac{m^2}{\mathcal{E}(1 - x^2)} \right] w_{vkm} \\ = -\Lambda_{vkm}(\zeta) w_{vkm}. \qquad (34) $$

When the centrifugal acceleration is taken into account, the eigenvalues $\Lambda_{vkm}$ and the generalised Hough functions (eigenfunctions) $w_{vkm}$ of the GLTE vary with the pseudo-radius $\zeta$. We choose to define our latitudinal ordering number $k$ to enumerate, for each $(\nu, m)$, the infinite set of solutions as in Lee & Saio (1997) and Mathis & Prat (2019) by considering the eigenvalues and eigenfunctions at the centre where they are not affected by the centrifugal acceleration since our mapping (Eq. (1)) is such that $r(\zeta \to 0, \theta) \to 0$.

We can check the validity of our model analytically by calculating the coefficients of the GLTE (Eq. (34)) in the spherical case (cf. Table 1) and by substituting these values into the GLTE. We then recover the standard Laplace tidal equation (SLTE) of Lee & Saio (1997):

$$ \mathcal{L}_{vm}^{\text{stand.}}[w_{vkm}^{\text{stand.}}] = \partial_x \left[ \frac{1 - x^2}{1 - \nu^2 x^2} \partial_x w_{vkm}^{\text{stand.}} \right] \\ + \left[ m\nu \frac{1 + \nu^2 x^2}{(1 - \nu^2 x^2)^2} - \frac{m^2}{(1 - x^2)(1 - \nu^2 x^2)} \right] w_{vkm}^{\text{stand.}} \\ = -\Lambda_{vkm}^{\text{stand.}} w_{vkm}^{\text{stand.}}. \qquad (35) $$

In this case, the eigenvalues $\Lambda_{vkm}^{\text{stand.}}$ are independent of the pseudo-radius $\zeta$. We note that at the limit $\nu = 0$, which corresponds to the non-rotating case, the SLTE reduces to the differential equation for the associated Legendre polynomials $P_{\ell,m}$ of

$$ \begin{gathered} \text{degree } \ell \text{ and order } m \text{ with } \Lambda_{0km}^{\text{stand.}} = \ell(\ell + 1) \text{ and } w_{0km}(x) = a_{\ell,m}P_{\ell,m}(x) \text{ where } \ell = |m| + k \text{ and} \\ a_{\ell,m} = (-1)^{\frac{m+|m|}{2}} \sqrt{\frac{2l+1}{4\pi} \frac{(l-|m|)!}{(l+|m|)!}}, \end{gathered} $$

is the normalisation coefficient.

In our work, we solve for each spheroid only one linear ordinary differential equation (the generalised Laplace tidal equation) which depends on the colatitude and in a parametric way on the pseudo-radius. Indeed, since we focus on low-frequency gravito-inertial waves, which are rapidly oscillating along the pseudo-radial direction, we adopt the JWKB approximation. It allows us to describe the spatial structure of the waves by the product of a rapidly oscillating plane-like wave function in the pseudo-radial direction multiplied by a bi-dimensional slowly varying envelope (Eqs. (28) and (29)) which is obtained by solving the GLTE (Eq. (34)) for each spheroidal shell for the normalised pressure and polarisation relationships for the different components of the velocity (Eqs. (30)–(32)). Therefore, we avoid having an ordinary differential equation to solve in the pseudo-radial direction. This method has been introduced and used by Mathis (2009) and Van Reeth et al. (2018) in the case of the generalisation of the TAR to the case of spherical differentially rotating stars and by Mathis & Prat (2019) and Henneco et al. (2021) in the case of the generalisation of the TAR for slightly deformed rotating stars. Finally, solving the GLTE allows us to compute the horizontal eigenvalues and the corresponding eigenfrequencies.

### 4.4. Asymptotic frequency and period spacing of low-frequency GWs

From Eqs. (33) and (34), we can identify the dispersion relation for low-frequency GWs within the TAR for strongly deformed
stars

$$k_{V;vkm}^2(\zeta) = \frac{N^2(\zeta)\mathcal{A}(\zeta)}{\omega_{km}^2}\Lambda_{vkm}(\zeta). \qquad (36)$$

To be able to derive the eigenfrequencies of low-frequency GIWs, we applied the quantisation relation in the vertical (pseudo-radial) direction used by Mathis (2009) and Mathis & Prat (2019) and introduced by Unno et al. (1989), Gough et al. (1993), and Christensen-Dalsgaard (1997)

$$\int_{\zeta_1}^{\zeta_2} k_{V;vkm} d\zeta = (n + 1/2)\pi, \qquad (37)$$

where $\zeta_1$ and $\zeta_2$ are the turning points of the Brunt–Väisälä frequency $N$ and $n$ is the radial order.

So, by substituting the dispersion relation (Eq. (36)) into the vertical quantisation one (Eq. (37)), we get the asymptotic expression for the frequencies of low-frequency GIWs

$$\omega_{nkm} = \frac{\int_{\zeta_1}^{\zeta_2} N(\zeta) \sqrt{\mathcal{A}(\zeta) \Lambda_{v_n km}(\zeta)} d\zeta}{(n + 1/2)\pi}, \qquad (38)$$

and the corresponding period

$$P_{nkm} = \frac{2(n + 1/2)\pi^2}{\int_{\zeta_1}^{\zeta_2} N(\zeta) \sqrt{\mathcal{A}(\zeta) \Lambda_{v_n km}(\zeta)} d\zeta}, \qquad (39)$$

with

$$v_n = v_{nkm} = \frac{2\Omega}{\omega_{nkm}}. \qquad (40)$$

We can thus compute the period spacing $\Delta P_{km} = P_{n+1km} - P_{nkm}$. We note that by applying a first-order Taylor development we can generalise the result obtained by Bouabid et al. (2013) in the spherical case

$$\Delta P_{km} \approx \frac{2\pi^2}{\int_{\zeta_1}^{\zeta_2} N \sqrt{\mathcal{A} \Lambda_{v_{n+1}km}} d\zeta \left(1 + \frac{1}{2} \frac{\int_{\zeta_1}^{\zeta_2} N \sqrt{\mathcal{A} \Lambda_{v_{nkm}} \frac{d\ln \Lambda_{v_{nkm}}}{d\ln \nu}} d\zeta}{\int_{\zeta_1}^{\zeta_2} N \sqrt{\mathcal{A} \Lambda_{v_{nkm}}}}\right)}, \qquad (41)$$

which will potentially allow us to probe the internal rotation of a large number of rapidly rotating stars taking into account their possibly strong flattening by the centrifugal acceleration and perform their detailed seismic modelling.

This result generalises and is very similar to the one obtained in the spherical case (Lee & Saio 1997; Bouabid et al. 2013) and in the weakly deformed case using a perturbative approach (Mathis & Prat 2019). The major interest of the eigenfrequencies, periods and periods spacing (Eqs. (38)–(41)) is that they can be easily evaluated for a large number of rapidly rotating stars and compared to observations which is a great asset to perform detailed seismic modelling. This flexibility is what makes the TAR a very interesting treatment.

## 5. Application to rapidly rotating early-type stars

In order to implement our equations and solve the GLTE, we must first compute an equilibrium model which gives the structure of the rotating star. The most advanced model of stellar evolution that we have and which is adapted to the formalism that we developed here is the ESTER model (Espinosa Lara & Rieutord 2013). It is unique of its kind because it calculates the effect of the centrifugal acceleration in a non-perturbative and a self-consistent way at any rotation rate.

### 5.1. ESTER models

The ESTER code is a 2D stellar structure code which takes into account the effects of rotation in a non-perturbative and coherent manner. The ESTER code makes it possible to calculate the equilibrium state of a star with a mass greater than 2 $M_\odot$ by solving the stationary equations of hydrodynamics for pressure, density, temperature and angular velocity for the entire volume. These realistic two-dimensional models of rapidly rotating stars are computed with the actual baroclinic meridional flows and differential rotation for intermediate-mass and massive stars. In this work, as a first step we work only with a uniform rotation (the mean value of the differential rotation over the pseudo-radius and the colatitude). The inclusion of differential rotation in the TAR formalism will be done in Paper II. The code uses spectral methods in the radial and latitudinal directions to transform the original system of non-linear partial differential equations in a system of non-linear algebraic equations (Chebyshev polynomials in the pseudo-radial direction and with spherical harmonics in the horizontal direction) (Rieutord et al. 2013).

In this study, we use $3M_\odot$ stellar models with a hydrogen mass fraction in the core $X_c = 0.7$ at [0%, 90%] of the Keplerian break-up rotation rate $\Omega_K$. Such a model is representative of observed stars near the zero-age main-sequence with asteroseismic probes of their angular velocity (Van Reeth et al. 2015, 2016, 2018; Pápics et al. 2015, 2017; Li et al. 2019a,b).

### 5.2. Domain of validity of the TAR

In this section, we study the validity domain of the two approximations (3.1 and 4.2) that are necessary to build the TAR in the case of rapidly rotating deformed bodies:

$$N^2(\zeta, \theta) \approx N^2(\zeta), \qquad (42)$$

$$\mathcal{A}(\zeta, \theta) \approx \mathcal{A}(\zeta), \qquad (43)$$

To visualise the problem, we first illustrate in Figs. 2 and 3 the functions $N^2(\zeta, \theta)$ and $\mathcal{A}(\zeta, \theta)$ computed with an ESTER model rotating at $\Omega/\Omega_K = 60\%$ (we use $r$ instead of $\zeta$ in these figures so we can show the spheroidal shape of the star caused by the centrifugal deformation). We can see from Fig. 2 that the convective region of the star is located between $\zeta = 0$ and $\zeta = 0.153$ and that the pseudo-radial variation of the Brunt–Väisälä frequency is the highest at the surface. From Fig. 3 we can see that the coefficient $\mathcal{A}$ varies the most, according to the pseudo-radius, at the centre of the star where we approach its singularity. In order to be able to apply our approximations (Eqs. (42) and (43)), we are interested in quantifying the latitudinal variation of this two quantities.

So to evaluate the committed error by adopting Eqs. (42) and (43) we represent in Figs. 4 and 5 the Brunt–Väisälä frequency and the coefficient $\mathcal{A}$ pseudo-radial profiles for different values of the colatitude $\theta_f$, at a fixed rotation rate $\Omega = 0.6\Omega_K$. In the bottom panels, we show the corresponding relative error

$$\delta X_{\theta_f}(\zeta) = \frac{X_{\text{approx}}(\zeta) - X(\zeta, \theta_f)}{X(\zeta, \theta_f)}, \qquad (44)$$

between the exact value $X(\zeta, \theta_f)$ given by the model and the approximated model $X_{\text{approx}}(\zeta)$ given by the weighted average of the exact value over the colatitude $\theta$

$$X_{\text{approx}}(\zeta) = \bar{X}(\zeta) = \int_{0}^{\pi/2} X(\zeta, \theta) \sin \theta d\theta, \qquad (45)$$
<figure><img src="image_3.png" /><figcaption>Fig. 2. Squared Brunt–Väisälä frequency $N^2(\zeta, \theta)$ profile for $\Omega/\Omega_K = 60\%$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Coefficient $A(\zeta, \theta)$ profile for $\Omega/\Omega_K = 60\%$.</figcaption></figure>

with

$$X \equiv \{N^2, A\},$$

where $\theta_f$ is a fixed value of the colatitude $\theta$.

Then, in order to determine the validity domain of these approximations as a function of the rotation rate and the

<figure><img src="image_5.png" /><figcaption>Fig. 4. Profile of the squared Brunt–Väisälä frequency $N^2$ and the relative error $\delta N_{\theta_f}^2$ of the adopted approximation (Eq. (42)) as a function of $\zeta$ at different colatitudes $\theta_f$ using an ESTER model rotating at 60% of the Keplerian breakup rotation rate (The light orange area indicates the margin of error which we allow and the red hatched area represents the convective region of the star).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Profile of the coefficient $A$ and the relative error $\delta A_{\theta_f}$ of the adopted approximation (Eq. (43)) as a function of $\zeta$ at different colatitudes $\theta_f$ using an ESTER model rotating at 60% of the Keplerian breakup rotation rate (The light orange area indicates the margin of error which we allow).</figcaption></figure>

pseudo-radius, we vary the rotation rate of the star $\Omega/\Omega_K$ (we use here 10 sample points for the rotation rate) and we calculate for each case the associate maximum relative error $\max_i \delta X_{\theta_i}(\zeta)$. Afterwards, we fix a maximum error rate equal to 10% and we deduce the pseudo-radius limit $\zeta_{\text{limit}}$ where the maximum relative error exceeds this threshold and we decide that the approximations become invalid. In Fig. 6, we display the pseudo-radius limit $\zeta_{\text{limit}}$ as a function of the rotation rate $\Omega/\Omega_K$. This curve defines the validity domains of the two approximations (Eqs. (42) and (43)) which consequently define the validity domain of the TAR.

We notice that the validity domains of the two approximations are more and more restricted (the pseudo-radius limit $\zeta_{\text{limit}}$ decrease) by increasing the rate of rotation and approaching the critical angular velocity. The tendency of the reduction of the pseudo-radius limit with the increasing of the rotation rate is expected since the action of the centrifugal acceleration is stronger further away from the axis of rotation. Thus, to be
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Validity domain of the TAR within the framework of ESTER models ($3 M_\odot$, $X_c = 0.7$) as a function of the rotation rate $\Omega/\Omega_K$ and the pseudo-radius $\zeta$ (with 90% degree of confidence).</figcaption></figure>

able to apply the TAR in the whole space domain, we shouldn't exceed a rotation rate $\Omega = 0.2\Omega_K$. It is therefore possible to apply the generalised TAR to early-type stars ($3 M_\odot$, $X_c = 0.7$) with an external radiative zone rotating up to 20% of the Keplerian critical angular velocity. This limit is lower than the one proposed by Mathis & Prat (2019) which was $\Omega = 0.4\Omega_K$. The limit proposed here should be more realistic since it is derived from a two-dimensional model which takes into account the centrifugal acceleration in a non-perturbative way whereas Mathis & Prat (2019) treated it in a perturbative manner. In the case where $\zeta_{\text{limit}} < 1$, generalised TAR should also be applicable to the radiative core of rapidly rotating (young) solar-type stars for which $\zeta_{\text{rad}} < \zeta_{\text{limit}}$, where $\zeta_{\text{rad}}$ is the pseudo radius of the radiation-convection interface. So, this formalism can be applied to all rapidly rotating bodies (stars of all types and planets) as long as we have a 2D model with access to the quantities that we used to derive our equations.

### 5.3. Eigenvalues and Hough functions

We solve the GLTE for different pseudo-radii, spin parameters, and rotation rates within the defined validity domain using the implementation based on Chebyshev polynomials presented and discussed by Wang et al. (2016) and then used by Mathis & Prat (2019) and Henneco et al. (2021).

Figure 7 shows the eigenvalues $\Lambda_{vkm}$ of the GLTE at $\zeta = \zeta_{\text{limit}}$ as a function of the spin parameter $v$ for $m = 1$ at two different rotation rates: $\Omega = 0.2\Omega_K$ the maximum rotation rate for which the generalised TAR is valid in all the radiative region ($\zeta_{\text{limit}} = 1$) and $\Omega = 0.9\Omega_K$ the highest rotation rate that we consider in this study for which the generalised TAR is valid in a restricted space domain ($\zeta_{\text{limit}} = 0.57$). Since $\Omega$ points in the direction of $\theta = 0$ and the oscillations are proportional to $e^{i(\omega t - m\phi)}$, a prograde (retrograde) oscillation has a positive (negative) value of the product $mv$. Therefore, since $m = 1$ the prograde and retrograde modes correspond respectively to the positive and negative values of $v$. We see here that the centrifugal deformation of the star causes a modest gradual shift in the eigenvalues as already identified in previous studies where the centrifugal acceleration has been treated in a perturbative way (Mathis & Prat 2019; Henneco et al. 2021).

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Spectrum of the GLTE as a function of the spin parameter $v$ at $\zeta = \zeta_{\text{limit}}$ and $m = 1$ for $\Omega = 0.2\Omega_K$ (above) and $\Omega = 0.9\Omega_K$ (below). Blue (respectively, orange) dots correspond to even (respectively, odd) eigenfunctions.</figcaption></figure>

We can distinguish two families of eigenvalues. Namely, the first one is gravity-like solutions ($\Lambda_{vkm} \ge 0$) which exist for any value of $v$ and we attach positive $k$'s to them. They correspond to gravity waves ($g$ modes) modified by rotation. The second one is Rossby-like solutions which exist only when $|v| > 1$ and we attach negative $k$'s to them. They only propagate in rotating stars. They correspond to Rossby modes ($r$ modes) if they are retrograde and have positive eigenvalues ($\Lambda_{vkm} > 0$) (Saio et al. 2018) and to overstable convective modes if they are prograde and have negative eigenvalues ($\Lambda_{vkm} < 0$). These modes are able to propagate in convective regions stabilised by the rotation under the action of the Coriolis acceleration (Lee & Saio 1997, 2020; Lee 2019). Even though the TAR is not a valid approximation in convective zones (Ogilvie & Lin 2004), we do find such propagative modes as solutions from our approach. Contrary to the spherical case we obtain here eigenvalues which depend slightly on the pseudo-radius $\zeta$ as shown in Fig. 8 where we represent the spectrum of the GLTE as a function of the pseudo-radius at two different rotation rates $\Omega/\Omega_K = \{20\%, 90\%\}$ for $m = 1$ and $v = 10$.

We focus now on the generalised Hough functions $w_{vkm}$ (normalised to unity) which vary with the pseudo-radius $\zeta$ and the horizontal coordinate $x$. This dependence is illustrated at $v = 10$ in Fig. 9 at two different rotation rates $\Omega/\Omega_K = \{20\%, 90\%\}$.
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Spectrum of the GLTE as a function of the pseudo-radius for *m* = 1 and *v* = 10 at Ω/Ω<sub>K</sub> = 0.2 (above) and Ω/Ω<sub>K</sub> = 0.9 (below). Blue (respectively, orange) dots correspond to even (respectively, odd) Hough functions.</figcaption></figure>

for gravity like-solutions (prograde dipole {k = 0, m = 1} modes) and Rossby like-solutions (retrograde Rossby modes with {k = -2, m = -1}). We can see that increasing the rotation rate from 20% to 90% of the critical angular velocity does hardly affect the Hough functions. This is caused by the fact that increasing the rotation rate implies decreasing the pseudo-radius limit (ζ<sub>limit</sub>) so we can no longer trust the computed Hough functions near the surface where the centrifugal effect is maximal. We can see clearly in this figure that the centrifugal acceleration slightly influences and modifies the eigenfunctions of the GLTE by introducing a new dependency on the pseudo-radial coordinate ζ, whereas in the spherically symmetric case they were only dependant on the latitudinal coordinate θ. In fact, the gravity-like solutions and the Rossby-like solutions shift as the distance from the border of the radiative zone of the star (ζ = 0.153) to its surface (ζ = 1) increases. More specifically, we observe a non-monotonous behaviour in Rossby-like solutions: these solutions migrate at first onwards, away from the equator (x = 0), causing a broadening of its shape. Next, by increasing the reduced latitudinal coordinate x this behaviour becomes less and less visible until it changes to the opposite one where these solutions migrate inwards, towards the equator, causing a narrowing of its shape. However, in gravity-like solutions we found a monotonous behaviour for which these

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Generalised Hough functions (normalised to unity) as a function of the horizontal coordinate *x* at different pseudo-radii from *ζ* = 0.153 (blue) to *ζ* = *ζ*<sub>limit</sub> (yellow) for *v* = 10, Ω/Ω<sub>K</sub> = 0.2 (above), and Ω/Ω<sub>K</sub> = 0.9 (below). The solid lines correspond to gravity-like solutions with {*k* = 1, *m* = 1}, while the dotted lines correspond to Rossby-like solutions with {*k* = −2, *m* = −1}.</figcaption></figure>

solutions migrate inwards, towards the equator, causing a narrowing of its shape. Contrary to the non-monotonous behaviour found in Rossby-like solutions, Henneco et al. (2021) found using the perturbative TAR a monotonous one with a migration towards the equator in the two types of solutions.

# 6. Asymptotic seismic diagnosis

## 6.1. Asymptotic period spacing pattern

To compute the period spacing patterns, we first calculated the asymptotic frequencies by following the method developed by Henneco et al. (2021) (the detailed steps of this computation can be found in their Appendix B). First, we compute for each radial order *n* of a given mode (*k*, *m*), the corresponding *v*<sub>*nkm*</sub> = 2Ω/*ω*<sub>*nkm*</sub> using the implicit equation of the asymptotic frequencies of low-frequency GIWs (Eq. (38)). Then, by taking their inverse and multiplying them with 2Ω, we retrieve the asymptotic frequencies *ω*<sub>*nkm*</sub>. The corresponding asymptotic frequencies in the inertial (observer's) frame *ω*<sub>*nkm*</sub><sup>in</sup> are then found using
the angular Doppler shift (Eq. (11))

$$ \omega_{nkm}^{\text{in}} = m\Omega + \underbrace{\omega_{nkm}}_{=2\Omega/\nu_{nkm}}. \qquad (46) $$

Now, we can easily calculate the asymptotic periods in the inertial frame ($P_{nkm}^{\text{in}} = 2\pi/\omega_{nkm}^{\text{in}}$) and the corresponding period spacing which is shown in Fig. 10 for the prograde dipole sectoral $\{k = 0, m = 1\}$ and the retrograde Rossby $\{k = -2, m = -1\}$ modes using the standard TAR (in spherical geometry) and the generalised TAR. The asymptotic periods within the standard TAR are computed using the following asymptotic frequencies deduced from the SLTE (Eq. (35))

$$ \omega_{nkm}^{\text{stand}} = \frac{\Lambda_{vkm}^{1/2}}{(n + 1/2)\pi} \int_{\zeta_1}^{\zeta_2} \frac{N(\zeta)}{\zeta} d\zeta. \qquad (47) $$

By comparing the period spacing pattern computed using the standard TAR and the generalised TAR in each of bottom panels of Fig. 10, we find that the spacing values of the $\{k = 0, m = 1\}$ and $\{k = -2, m = -1\}$ modes increase with the decrease in radial orders (low (high) radial orders equivalent to short (long) pulsation periods for g modes and long (short) pulsation periods for r modes) under the influence of the centrifugal deformation. The effect of the centrifugal acceleration on the period spacing that we found here is different from the one discussed in Henneco et al. (2021) who found an increase in the period spacing of the two modes. Another notable difference between the generalised TAR and the standard (Lee & Saio 1997) and the perturbative (Henneco et al. 2021) TAR is in the domain of validity and applicability, where the general treatment adopted in this work should be more reliable than the standard and the perturbative ones. Figure 10 reveals the difference in the period spacing values computed with the standard and the generalised formalisms and also offers a comparison of this difference with typical observational uncertainties of measured period spacing values from modern space photometric observations. In fact, Van Reeth et al. (2015) shows that for a sample of 40 γ Dor stars the error margins of the period spacings have values up to 25 seconds for the prograde dipole modes (larger than our period spacing differences) and between 20 and 50 seconds for the retrograde Rossby modes (smaller than our period spacing differences).

## 6.2. Detectability of the centrifugal deformation effect

To evaluate the detectability of the centrifugal effect with space-based photometric observations, we first compute the frequency differences $\Delta f_{\text{centrifugal}}$ between asymptotic frequencies calculated in the standard TAR (TARs) and those computed in the generalised TAR (TARg) through Eq. (38):

$$ \Delta f_{\text{centrifugal}}(n) = |f_{\text{TARg}}(n) - f_{\text{TARs}}(n)|. \qquad (48) $$

Then, by comparing the obtained frequency differences with the frequency resolutions ($f_{\text{res}} = 1/T_{\text{obs}}$) of nominal 4-year Kepler (Borucki et al. 2009) and TESS CVZ (Transiting Exoplanet Survey Satellite Continuous Viewing Zone; Ricker et al. 2014) light curves covering quasi-continuously observation times of $T_{\text{obs}} = 4$ yrs and $T_{\text{obs}} = 351$ d, respectively, we can deduce the radial orders $n$ for which the frequency differences are, at least in principle, expected to be detectable. In Fig. 11, we display these results for the $\{k = 0, m = 1\}$ and $\{k = -2, m = -1\}$ modes for a stellar model rotating at $0.2\Omega_K$. We can see that the centrifugal effect for the $\{k = 0, m = 1\}$ mode is detectable for all radial

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 10. Period spacing pattern in the inertial frame for $\{k = 0, m = 1\}$ (above) and $\{k = -2, m = -1\}$ (below) modes at $\Omega/\Omega_K = 20\%$ within the standard (red line) and the generalised (blue line) TAR frameworks. The bottom panel for each mode shows the differences with respect to the standard TAR period spacing pattern (the sudden fluctuations in the period spacing pattern are caused by the numerical noise which is introduced by the numerical derivatives of the mapping with respect to $\zeta$ and $\theta$ used in the resolution of the GLTE).</figcaption></figure>

orders using Kepler observations, but detectable only for $n < 41$ using TESS observations (yellow dashed line in Fig. 11). Indeed even if for $n \in [20, 40]$ using nominal TESS CVZ light curves and $n \in [75, 100]$ using nominal Kepler light curves we have $\Delta f_{\text{centrifugal}}/f_{\text{res}} \gtrsim 1$, so we would expect that the detectability of the centrifugal effect for these radial orders is in principle guaranteed. However, since the frequency differences is nearly equal to the resolution one ($\Delta f_{\text{centrifugal}} \approx f_{\text{res}}$), this effect become difficult to detect in this case. Regarding the $\{k = -2, m = -1\}$ mode, we can see that the centrifugal effect is less detectable in comparison with the $\{k = 0, m = 1\}$ mode but remains in principle observable (yet difficult to observe because $\Delta f_{\text{centrifugal}} \approx f_{\text{res}}$)
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. Detectability of {k = 0, m = 1} (above) and {k = -2, m = -1} (below) modes as a function of the radial order n at a rotation rate Ω/Ω<sub>K</sub> = 20% based on the frequency resolution of TESS (brown line) and Kepler (indigo line).</figcaption></figure>

using Kepler observations for small radial orders (n ∈ [5, 8]) and for n ∈ [16, 45].

Our results in Figs. 10 and 11 are formal comparisons in the sense that uncalibrated astrophysical effects occurring in stellar models (e.g., the presence of rotational mixing and atomic diffusion effects) usually imply larger uncertainties in the prediction of the period spacing values than the differences shown in these figures, as outlined in Aerts et al. (2018). For this reason, our theoretical findings in the current work offer a key result in terms of reliability of future asteroseismic analyses of observed pulsators. The centrifugal effect is not significant when assessed within the validity domain of the generalised TAR which is consistent with the results of Ouazzani et al. (2017) and Henneco et al. (2021).

## 7. Evaluation of the terms hierarchy imposed by the TAR within uniformly rotating deformed stars

The hierarchy of terms imposed by the TAR can be summarised by the following frequency hierarchy:

$$2\Omega \ll N, \qquad (49)$$

$$\omega \ll N, \qquad (50)$$

which ensures the other hierarchies. In fact, the vertical wave vector is larger than the horizontal one ($|k_H| \ll |k_V|$) which comes from the dispersion relation of GIWs (Lee & Saio 1997)

$$\omega^2 \approx \frac{N^2 k_H^2 + (2\Omega \cdot k)^2}{k^2}. \qquad (51)$$

In addition, the latitudinal component of the rotation vector in the momentum equation Ω<sub>H</sub> = Ω sin θ is neglected for all colattitudes because (2Ω · k) ≈ 2Ω<sub>V</sub>k<sub>V</sub>. And finally, the wave velocity is almost horizontal ($|v^\perp| \ll \{|v^\theta|, |v^\psi|\}$). Indeed, by adopting the anelastic approximation we obtain

$$\nabla \cdot (\rho v) \approx 0 \Rightarrow k \cdot v = k_V v_V + k_H \cdot v_H \approx 0 \\ \Rightarrow \frac{v_V}{v_H} \approx -\frac{k_H}{k_V} \ll 1. \qquad (52)$$

The anelastic approximation filters out acoustic waves which have higher frequencies, and is justified by the fact that we study only low frequency waves (Eq. (50)) within the framework of the TAR treatment.

So in order to verify whether the TAR is still valid in uniformly rotating deformed stars, we discuss the frequency hierarchy (Eqs. (49) and (50)) using the Brunt–Väisälä frequency and rotation profiles from ESTER models and using the asymptotic frequencies calculated in Sect. 6.

### 7.1. The strong stratification assumption: $2\Omega \ll N$

We evaluate the term N/2Ω using the Brunt–Väisälä frequency and the angular velocity computed with ESTER models (3 M<sub>☉</sub>, X<sub>c</sub> = 0.7) for rotation rates Ω/Ω<sub>K</sub> ∈ [0.1, 0.9]. As shown in Fig. 12, the strong stratification assumption is valid only in the radiative zone far from the border between convection and radiation (ζ ≥ 0.2 at Ω/Ω<sub>K</sub> ≤ 0.2). Beyond this critical value, the Brunt–Väisälä frequency and the frequency of rotation have a very close order of magnitude so we can no longer adopt the strong stratification approximation for stars close to the zero-age main-sequence. As stars evolve, however, they build up a strong gradient of molecular weight, increasing the value of N(ζ) in the region near their convective core. Therefore, the strong stratification assumption becomes better fulfilled as stars evolve along the main sequence because the spin-up of their core remains modest (Li et al. 2020; Aerts 2021), while N(ζ) increases.

### 7.2. The low frequency assumption: $\omega \ll N$

We evaluate the term N/ω using the Brunt–Väisälä frequency computed by ESTER models and the asymptotic frequencies for {k = 0, m = 1} and {k = -2, m = -1} modes calculated in Sect. 6. As shown in Fig. 13, the low frequency assumption is valid for {k = -2, m = -1} mode for all radial orders and pseudoradii while for {k = 0, m = 1} mode this assumption is valid for all radial orders only if we are far from the transition layer between the convective core and the radiative envelope. If this is not the case, we have to work with high radial orders (n > 20). Below this critical value, the Brunt–Väisälä frequency and the wave frequency have a very close order of magnitude so we can no longer adopt the low frequency approximation.

## 8. Discussion and conclusions

The aim of this work was to study the possibility of carrying out a new generalisation of the TAR that takes into account the
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** N/2Ω term as a function of the pseudo-radius $\zeta$ at different rotation rates $\Omega/\Omega_K$ (the hatched area represents the convective region of the star).</figcaption></figure>

centrifugal acceleration in the case of strongly deformed, rapidly rotating stars and planets. We approached this issue by deriving the generalised Laplace Tidal equation in a general spheroidal coordinate system. We relied on two necessary assumptions (described by Eqs. (42) and (43)) that are applicable in a specific domain, which will then define the validity domain of the generalised TAR. The equation that we derive has a similar form to the one obtained when the TAR is applied to weakly rotating spherical stars (Lee & Saio 1997), differentially rotating spherical stars (Mathis 2009; Van Reeth et al. 2018) and moderately rotating, weakly deformed stars (Mathis & Prat 2019). This new formalism can be applied in the radiative region of all types of stars and planets.

Using this theoretical framework with 2D ESTER models we can define the validity domain of the TAR with a chosen confidence level as a function of the rotation rate and the pseudo-radius. We apply this general formalism to rapidly rotating early-type stars ($M = 3M_\odot$, $X_c = 0.7$), where we found that the generalised TAR is only applicable for all rotation rates for values of the pseudo-radius smaller than 0.57 ($\zeta_{\text{limit}} < 0.57$). To calculate the asymptotic pulsation frequencies we then have to integrate the dispersion relation on the pulsation mode cavity (the radiative zone) (Eq. (37)) where $\zeta \in [\zeta_1, \zeta_2]$ (in our considered case $\zeta_1 = 0.153$ and $\zeta_2 = 1$). This condition will limit the rotation rates where the TAR is valid. We can therefore deduce a rotation rate limit ($\Omega_{\text{limit}} = 0.2\Omega_K$ in our case) where the TAR ceases to be applicable with a fixed level of confidence (90% in our study).

We found a difference between the $\text{Islimit}$ of the TAR applicability defined in Mathis & Prat (2019) which was $\Omega_{\text{limit}} = 0.4\Omega_K$ and the one defined here. The current new limit should be more realistic and accurate since it is first derived from a 2D model which takes the impact of the centrifugal acceleration into account in a non-perturbative way. Next, the validity domain of the adopted approximations is studied carefully in more details and with a level of confidence of 90%.

The centrifugal acceleration affects two aspects within our new generalised theoretical formalism compared to the standard application of the TAR (Lee & Saio 1997). On the one hand, the stellar structure as a whole becomes deformed which causes the usual radial coordinate $r$, the hydrostatic background and thus the Brunt-Väisälä frequency to have a dependency not only on the pseudo-radius $\zeta$ but also on the colatitude $\theta$. This introduces

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 13.** N/ω term as a function of the pseudo-radius $\zeta$ and the radial order $n$ for $\{k = 0, m = 1\}$ (top panel) and $\{k = -2, m = -1\}$ (bottom panel) modes at $\Omega/\Omega_K = 20\%$.</figcaption></figure>

a new challenge in the decoupling of the vertical and the horizontal dynamics and the applicability of the quantisation relation unless we adopt adequate approximations which will restrict the validity domain of our formalism. On the other hand, the dynamics of the GIWs represented by the GLTE, is altered and get a radial dependence compared to the one in the spherical case represented by the SLTE. Then, we apply our new formalism to a rapidly rotating early-type star ($M = 3M_\odot$, $X_c = 0.7$, $\Omega/\Omega_K = 20\%$). In this case, we observe a non-monotonous geometrical behaviour in the Rossby-like solutions and a monotonous one (only an inward migration) in the gravity-like solutions. In fact, near the equator the Rossby-like solutions migrate onwards, causing a broadening of its shape. But far from the equator this behaviour becomes less and less visible until it changes to the opposite one where these solutions migrate inwards, causing a narrowing of its shape. We found also that the centrifugal acceleration increases the period spacing values for low radial orders (equivalent to short pulsation periods for $g$ modes and long pulsation periods for $r$ modes) and decrease slightly for high radial orders (equivalent to long pulsation periods for $g$ modes and short pulsation periods for $r$ modes). But yet the effects of the centrifugal acceleration (within the validity domain of the generalised TAR) are of the order of or smaller than the effects of other missing stellar physical processes in our models (e.g., accurate rotational mixing, atomic diffusion,
etc.). This can mask the effect of the centrifugal acceleration in
forward asteroseismic modelling analyses, even when it is theo-
retically detectable.

This new generalisation of the TAR can also be used to study
the tidal dissipation induced by low-frequency GIWs in rapidly
rotating deformed stars and planets (Braviner & Ogilvie 2014)
and the angular momentum transport with a formalism directly
implantable in ESTER models. In addition, we can apply our
formalism to different 2D ESTER model grids and then compare
the obtained results with direct computations using 2D oscilla-
tion codes TOP (Reese et al. 2021) and ACOR (Ouazzani et al.
2017).

The next steps will be first to take into account the general
differential rotation in the deformed case. So far, Ogilvie & Lin
(2004) and Mathis (2009) included the effects of general differ-
ential rotation on low-frequency GIWs within the framework of
the TAR for spherical stars. Then, their formalism has been suc-
cessfully applied in Van Reeth et al. (2018) to derive the varia-
tion of the asymptotic period spacing in the case of a weak radial
differential rotation and evaluate the sensitivity of GIWs to the
effect of differential rotation in spherical stars. Next, we will
take into account the magnetic field in a non-perturbative way
(Mathis & de Brye 2011, 2012). So far, Prat et al. (2019, 2020)
and Van Beeck et al. (2020) have focused on the case where
magnetic fields are weak enough to be treated within a pertur-
bative treatment to study the effects of a magnetic field on the
seismic parameters of *g* modes which become magneto-gravito
inertial modes.

So in the near future, magneto-gravito inertial waves have
to be studied in the general case of rapidly and differentially
rotating stars and planets. As a first step, in a forthcoming paper
(Dhouib et al. 2021, Paper II) we will examine if the TAR can be
generalised to the case of strongly deformed differentially rotat-
ing stars and planets.

*Acknowledgements.* We thank the referee, Pr. H.Shibahashi, for his construc-
tive report that allowed us to improve the quality of our work and of the article.
We thank the ESTER code developers for their efforts and for making their codes
publicly available. We are also grateful to Pr. C.Aerts who kindly commented on
the manuscript and suggested improvements. H.D., V.P. and S.M. acknowledge
support from the European Research Council through ERC grant SPIRE 647383.
H.D. and S.M. acknowledge support from the CNES PLATO grant at CEA/DAp.
T.V.R. gratefully acknowledges support from the Research Foundation Flanders
(FWO) under grant agreement No. 12ZB620N.

References

Aerts, C. 2021, *Rev. Mod. Phys.*, **93**, 015001

Aerts, C., Christensen-Dalsgaard, J., & Kurtz, D. W. 2010, *Asteroseismology* (Berlin: Springer)

Aerts, C., Molenberghs, G., Michielsen, M., et al. 2018, ApJS, **237**, 15

Aerts, C., Mathis, S., & Rogers, T. M. 2019, ARA&A, **57**, 35

Berthomieu, G., Gonczi, G., Graff, P., Provost, J., & Rocca, A. 1978, A&A, **70**, 597

Bildsten, L., Ushomirsky, G., & Cutler, C. 1996, ApJ, **460**, 827

Bonazzola, S., Gourgoulhon, E., & Marck, J.-A. 1998, Phys. Rev. D, **58**, 104020

Borucki, W., Koch, D., Batalha, N., et al. 2009, in *Transiting Planets*, eds. F. Pont, D. Sasselov, M. J. Holman, et al., IAU Symp., **253**, 289

Bouabid, M. P., Dupret, M. A., Salmon, S., et al. 2013, MNRAS, **429**, 2500

Braviner, H. J., & Ogilvie, G. I. 2014, MNRAS, **441**, 2321

Christensen-Dalsgaard, J. 1997, *Stellar Oscillations*

Cowling, T. G. 1941, MNRAS, **101**, 367

Dhouib, H., Prat, V., Van Reeth, T., & Mathis, S. 2021, A&A, submitted
(Paper II)

Dintrans, B., & Rieutord, M. 2000, A&A, **354**, 86

Dintrans, B., Rieutord, M., & Valdettaro, L. 1999, J. Fluid Mech., **398**, 271

Eckart, C. 1960, *Hydrodynamics of Oceans and Atmospheres* (Oxford: Pergamon Press)

Espinosa Lara, F., & Rieutord, M. 2013, A&A, **552**, A35

Fröman, N., & Fröman, P. O. 1965, *JWKB approximation: contributions to the theory* (Amsterdam: North-Holland Publishing Company)

Gough, D. 1993, in *Les Houches Session XLVII*, eds. J. P. Zahn, & J. Zinn-Justin (Amsterdam: Elsevier), 399

Henneco, J., Van Reeth, T., Prat, V., et al. 2021, A&A, **648**, A97

Hough, S. S. 1898, *Philos. Trans. R. Soc. London, Ser. A*, **191**, 139

Laplace, P. S. 1799, *Traité de Mécanique Céleste* (Paris: Imprimerie de Crapelet)

Lee, U. 2019, MNRAS, **484**, 5845

Lee, U., & Baraffe, I. 1995, A&A, **301**, 419

Lee, U., & Saio, H. 1993, MNRAS, **261**, 415

Lee, U., & Saio, H. 1997, ApJ, **491**, 839

Lee, U., & Saio, H. 2020, MNRAS, **497**, 4117

Lee, U., Neiner, C., & Mathis, S. 2014, MNRAS, **443**, 1515

Li, G., Bedding, T. R., Murphy, S. J., et al. 2019a, MNRAS, **482**, 1757

Li, G., Van Reeth, T., Bedding, T. R., Murphy, S. J., & Antoci, V. 2019b, MNRAS, **487**, 782

Li, G., Van Reeth, T., Bedding, T. R., et al. 2020, MNRAS, **491**, 3586

Lignières, F., Rieutord, M., & Reese, D. 2006, A&A, **455**, 607

Mathis, S. 2009, A&A, **506**, 811

Mathis, S., & de Brye, N. 2011, A&A, **526**, A65

Mathis, S., & de Brye, N. 2012, A&A, **540**, A37

Mathis, S., & Prat, V. 2019, A&A, **631**, A26

Mathis, S., Talon, S., Pantillon, F. P., & Zahn, J. P. 2008, Sol. Phys., **251**, 101

Neiner, C., Floquet, M., Samadi, R., et al. 2012, A&A, **546**, A47

Neiner, C., Lee, U., Mathis, S., et al. 2020, A&A, **644**, A9

Ogilvie, G. I., & Lin, D. N. C. 2004, ApJ, **610**, 477

Ogilvie, G. I., & Lin, D. N. C. 2007, ApJ, **661**, 1180

Ouazzani, R. M., Dupret, M. A., & Reese, D. R. 2012, A&A, **547**, A75

Ouazzani, R.-M., Salmon, S. J. A. J., Antoci, V., et al. 2017, MNRAS, **465**, 2294

Ouazzani, R. M., Marques, J. P., Goupil, M. J., et al. 2019, A&A, **626**, A121

Pantillon, F. P., Talon, S., & Charbonnel, C. 2007, A&A, **474**, 155

Pápics, P. I., Tkachenko, A., Aerts, C., et al. 2015, ApJ, **803**, L25

Pápics, P. I., Tkachenko, A., Van Reeth, T., et al. 2017, A&A, **598**, A74

Pedersen, M. G., Aerts, C., Pápics, P. I., & Rogers, T. M. 2018, A&A, **614**, A128

Prat, V., Mathis, S., Buysschaert, B., et al. 2019, A&A, **627**, A64

Prat, V., Mathis, S., Neiner, C., et al. 2020, A&A, **636**, A100

Reese, D., Lignières, F., & Rieutord, M. 2006, A&A, **455**, 621

Reese, D. R., MacGregor, K. B., Jackson, S., Skumanich, A., & Metcalfe, T. S. 2009, A&A, **506**, 189

Reese, D. R., Mirouh, G. M., Espinosa Lara, F., Rieutord, M., & Putigny, B. 2021, A&A, **645**, A46

Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, in *Space Telescopes and Instrumentation 2014: Optical, Infrared, and Millimeter Wave*, eds. J. Oschmann, M. Jacobus, M. Clampin, G. G. Fazio, H. A. MacEwen, et al., SPIE Conf. Ser., **9143**, 914320

Rieutord, M., & Espinosa Lara, F. 2013, in *Ab Initio Modelling of Steady Rotating Stars*, eds. M. Goupil, K. Belkacem, C. Neiner, F. Lignières, & J. J. Green, 865, 49

Rieutord, M., Dintrans, B., Lignières, F., Corbard, T., & Pichon, B. 2005, in *Semaine de l'Astrophysique Francaise*, eds. F. Casoli, T. Contini, J. M. Hameury, & L. Pagani SF2A-2005, 759

Rogers, T. M. 2015, AJ, **815**, L30

Rogers, T. M., Lin, D. N. C., McElwaine, J. N., & Lau, H. H. B. 2013, ApJ, **772**, 21

Saio, H. 1981, ApJ, **244**, 299

Saio, H., Kurtz, D. W., Murphy, S. J., Antoci, V. L., & Lee, U. 2018, MNRAS, **474**, 2774

Spiegel, E. A., & Veronis, G. 1960, ApJ, **131**, 442

Townsend, R. H. D. 2003, MNRAS, **340**, 1020

Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1989, *Nonradial oscillations of stars* (Tokyo: University of Tokyo Press)

Van Beeck, J., Prat, V., Van Reeth, T., et al. 2020, A&A, **638**, A149

Van Reeth, T., Tkachenko, A., Aerts, C., et al. 2015, A&A, **574**, A17

Van Reeth, T., Tkachenko, A., & Aerts, C. 2016, A&A, **593**, A120

Van Reeth, T., Mombarg, J. S. G., Mathis, S., et al. 2018, A&A, **618**, A24

Wang, H., Boyd, J. P., & Akmaev, R. A. 2016, Geosci. Model Dev., **9**, 1477
# Appendix A: Derivation of the linearised momentum and continuity equations in the spheroidal basis

## A.1. The momentum equation

The linearised momentum equation for an inviscid uniformly rotating fluid is written as

$$ (\partial_t + \Omega \partial_\varphi) \boldsymbol{v} + 2\Omega \wedge \boldsymbol{v} = -\frac{1}{\rho_0} \nabla P' + \frac{\rho'}{\rho_0^2} \nabla P_0 - \nabla \Phi'. \quad (A.1) $$

Using the spheroidal basis ($\mathbf{a}_\zeta, \mathbf{a}_\theta, \mathbf{a}_\varphi$) defined in Eq. (2) we can express the velocity field as follows

$$ \boldsymbol{v} = v^i \mathbf{a}_i = v^i \alpha_i \mathbf{E}_i = v^i \alpha_i g_{ij} \mathbf{E}^j, \quad (A.2) $$

with $\alpha_\zeta = \frac{\zeta^2}{r^2 r_\zeta}$, $\alpha_\theta = \frac{\zeta}{r^2 r_\zeta}$, $\alpha_\varphi = \frac{\zeta}{r^2 r_\zeta \sin\theta}$, ($\mathbf{E}^\zeta, \mathbf{E}^\theta, \mathbf{E}^\varphi$) is the dual basis (the contravariant basis) and $g_{ij}$ is the covariant metric tensor defined as follows

$$ g_{ij} = \mathbf{E}_i \cdot \mathbf{E}_j = \begin{pmatrix} r_\zeta^2 & r_\zeta r_\theta & 0 \\ r_\zeta r_\theta & r^2 + r_\theta^2 & 0 \\ 0 & 0 & r^2 \sin^2 \theta \end{pmatrix}. \quad (A.3) $$

We write also the uniform angular velocity in the natural covariant basis as

$$ \Omega = \Omega(\cos\theta e_r - \sin\theta e_\theta) \quad (A.4) $$

$$ = \Omega \underbrace{\left( \frac{\cos \theta}{r_\zeta} + \frac{\sin \theta r_\theta}{r r_\zeta} \right)}_{=\Omega^\zeta} \mathbf{E}_\zeta + \underbrace{\left( -\frac{\Omega \sin \theta}{r} \right)}_{=\Omega^\theta} \mathbf{E}_\theta. \quad (A.5) $$

We then express the Coriolis acceleration in the spheroidal basis as follows

$$ 2\Omega \wedge \boldsymbol{v} = 2\varepsilon_{ilk} \Omega^l \alpha_k v^k \mathbf{E}^i, \quad (A.6) $$

where $\varepsilon^{ilk}$ is the Levi-Civita tensor

$$ \varepsilon_{ilk} = \sqrt{|g|}[i, l, k], \quad (A.7) $$

with $|g| = \det(g_{ij}) = r^4 r_\zeta^2 \sin^2 \theta$ and

$$ [i, l, k] = \begin{cases} +1 & \text{if } (i, l, k) \text{ is an even permutation of } (\zeta, \theta, \varphi) \\ -1 & \text{if } (i, l, k) \text{ is an odd permutation of } (\zeta, \theta, \varphi) \\ 0 & \text{otherwise} \end{cases}. \quad (A.8) $$

Finally, we express the gradient of the scalar function $f = \{P', P_0, \Phi'\}$ in the spheroidal basis as

$$ \nabla f = \partial_i f \mathbf{E}^i. \quad (A.9) $$

With these definitions, it is now possible to give an explicit expression of the momentum equation (Eq. (A.1)) in the spheroidal basis

$$ \begin{aligned} & (\partial_t + \Omega \partial_\varphi) \alpha_j v^j g_{ji} + 2\varepsilon_{ilk} \Omega^l \alpha_k v^k \\ & = -\frac{1}{\rho_0} \partial_i P' + \frac{\rho'}{\rho_0^2} \partial_i P_0 - \partial_i \Phi', \end{aligned} \quad (A.10) $$

for $i = \zeta, i = \theta$ and $i = \varphi$ we get the pseudo-radial (Eq. (3)), the latitudinal (Eq. (4)), and the azimuthal (Eq. (5)) components of the momentum equation, respectively.

## A.2. The continuity equation

The linearised continuity equation is given by

$$ (\partial_t + \Omega \partial_\varphi) \rho' + \nabla \cdot (\rho_0 \boldsymbol{v}) = 0, \quad (A.11) $$

then by expressing the divergence in the spheroidal basis as

$$ \nabla \cdot (\rho_0 \boldsymbol{v}) = \frac{1}{\sqrt{|g|}} \partial_j (\sqrt{|g|} \rho_0 \alpha_j v^j), \quad (A.12) $$

we get the Eq. (6).

# Appendix B: Visual presentation of the terms involved in the derivation of the GLTE

We represent in Fig. B.1 the profiles of the 8 terms ($\mathcal{A}$ is already represented in the top panel of Fig. 5) involved in the derivation of the GLTE (Eq. (34)) in the general case of spheroidal geometry presented in the Table 1 as a function of the pseudo-radius $\zeta$ at different reduced colatitudes $x$. In the main text (Sect. 5.2), we mainly focused on the coefficient $\mathcal{A}$ and the Brunt-Vaisälä frequency $N$ because they are the important quantities to be able to build the generalised TAR formalism. The oscillating behaviour of some terms is caused by the fact that the mapping (Eq. (1)) is different in each subdomain of the model which is continuous and derivable, but not twice differentiable with respect to $\zeta$ between subdomains.
<figure><img src="image_18.png" /><figcaption>**Fig. B.1.** Profiles of the coefficients $\mathcal{B}, \mathcal{C}, \mathcal{D}, \mathcal{E}, \partial_\alpha \mathcal{B}, \partial_\alpha \mathcal{C}, \partial_\alpha \mathcal{D}$ and $\partial_\alpha \left(\frac{C}{\mathcal{E}\mathcal{B}}\right)$ involved in the derivation of the GLTE (Eq. (34)) as a function of $\zeta$ at different reduced colatitudes $\alpha$ using an ESTER model ($M = 3M_\odot$, $X_c = 0.7$) rotating at 20% of the Keplerian break-up rotation rate for a spin parameter $\nu = 10$.</figcaption></figure>