The traditional approximation of rotation for rapidly rotating stars and planets. II. Deformation and differential rotation

We examine the dynamics of low-frequency gravito-inertial waves (GIWs) in differentially rotating deformed radiation zones in stars and planets by generalising the traditional approximation of rotation (TAR). The TAR treatment was built on the assumptions that the star is spherical and uniformly rotating. However, it has been generalised in our previous work by including the effects of the centrifugal deformation using a non-perturbative approach in the uniformly rotating case. We aim to carry out a new generalisation of the TAR treatment to account for the differential rotation and the strong centrifugal deformation simultaneously. We generalise our previous work by taking into account the differential rotation in the derivation of our complete analytical formalism that allows the study of the dynamics of GIWs in differentially and rapidly rotating stars. We derived the complete set of equations that generalises the TAR, simultaneously taking the full centrifugal acceleration and the differential rotation into account. Within the validity domain of the TAR, we derived a generalised Laplace tidal equation for the horizontal eigenfunctions and asymptotic wave periods of the GIWs, which can be used to probe the structure and dynamics of differentially rotating deformed stars with asteroseismology. A new generalisation of the TAR, which simultaneously takes into account the differential rotation and the centrifugal acceleration in a non-perturbative way, was derived. This generalisation allowed us to study the detectability and the signature of the differential rotation on GIWs in rapidly rotating deformed stars and planets. We found that the effects of the differential rotation in early-type deformed stars on GIWs is theoretically largely detectable in modern space photometry using observations from $\textit{Kepler}$ and TESS.


Introduction
Understanding how angular momentum and chemicals are transported in the interiors of stars (and planets) along their evolution is one of the key challenges of modern stellar (and planetary) astrophysics. Indeed, rotation modifies their structure, their chemical stratification, their internal flows and magnetism, and their mass losses and winds (e.g . Maeder 2009;Mathis 2013;Aerts et al. 2019, and references therein). In this quest, asteroseismology has bought a fundamental breakthrough by demonstrating that all stars are the seat of a strong extraction of angular momentum during their evolution in comparison with the predictions by stellar models taking the rotation into account following the standard rotational transport and mixing theory (Eggenberger et al. 2012;Marques et al. 2013;Ceillier et al. 2013;Cantiello et al. 2014;Ouazzani et al. 2019). This was first obtained thanks to mixed pulsation modes splitted by rotation propagating in evolved low-and intermediate-mass stars (Beck et al. 2012;Mosser et al. 2012;Deheuvels et al. 2012;Beck et al. 2014;Deheuvels et al. 2014Deheuvels et al. , 2015Di Mauro et al. 2016;Triana et al. 2017;Gehan et al. 2018;Beck et al. 2018;Tayar et al. 2019;Deheuvels et al. 2020). Then, observations of oscillation modes in F-and A-type stars (Kurtz et al. 2014;Saio et al. 2015;Bedding et al. 2015;Keen et al. 2015;Van Reeth et al. 2015Schmid & Aerts 2016;Murphy et al. 2016;Sowicka et al. 2017;Guo et al. 2017;Van Reeth et al. 2018;Saio et al. 2018;Mombarg et al. 2019;Li et al. 2019Li et al. , 2020Ouazzani et al. 2020;Saio et al. 2021) and in B-type stars Triana et al. 2015;Moravveji et al. 2016;Pápics et al. 2017;Kallinger et al. 2017;Buysschaert et al. 2018;Szewczuk & Daszyńska-Daszkiewicz 2018;Pedersen et al. 2021;Szewczuk et al. 2021) provided us new Rosetta stones to constrain the transport of angular momentum in the whole Hertzsprung-Russell diagram. More particularly, this pushes gravity-and gravito-inertial mode pulsators such as γ-Doradus and SPB stars at the forefront of this research. For instance, recent theoretical developments have demonstrated how it is possible to probe stellar internal rotation in γ-Doradus stars from their surface to their convective core (Ouazzani et al. 2020;Saio et al. 2021). These stars are rapid rotators for a large proportion of them. Therefore, it is necessary to study gravito-inertial modes. These modes are gravity modes, which propagate only in stably stratified stellar radiation zones Article number, page 1 of 13 arXiv:2110.03619v1 [astro-ph.SR] 7 Oct 2021 A&A proofs: manuscript no. main under the action of the restoring buoyancy force in the absence of rotation, which are modified by rotation. If their frequency is super-inertial (i.e. above the inertial frequency 2Ω, Ω being the stellar angular velocity), they are propagating in stellar radiation zones and evanescent in convective regions. If their frequency is sub-inertial (below 2Ω) they propagate in an equatorial belt in radiation zones and they become propagative inertial waves in convective regions (e.g. Dintrans & Rieutord 2000;Mathis et al. 2014). The challenge of studying these waves is that the equation describing their dynamics are intrinsically bi-dimensional and non-separable (Dintrans et al. 1999;Prat et al. 2016;Mirouh et al. 2016;Prat et al. 2018). This makes the development of seismic diagnosis difficult analytically (Prat et al. 2017) or expansive in computation time when using 2D oscillation and stellar structure codes (e.g. Ouazzani et al. 2017;Reese et al. 2021) in the general case.
In this framework, the traditional approximation of rotation (TAR), which has been first introduced in geophysics (Eckart 1960) to treat the propagation of gravito-inertial waves (GIWs) in the case where the Coriolis acceleration can be neglected in front of the buoyancy force in the direction of the entropy and chemical stable stratifications, is very useful. Indeed, it allows one to consider that the velocity of gravito-inertial waves are mostly horizontal and to neglect the latitudinal component of the rotation vector in the momentum equation. That leads to decoupling the vertical and horizontal directions by keeping the non-perturbative action of the Coriolis acceleration in the latitudinal and azimuthal directions and neglecting it along the vertical one. Propagation equations become separable as in the non-rotating case (e.g. Bildsten et al. 1996;Lee & Saio 1997). In addition, the derivation of powerful and flexible seismic diagnostics using the period spacing between consecutive high radial order gravito-inertial modes becomes possible in uniformly rotating spherical stars (Bouabid et al. 2013). These period spacing provide constrains on properties of the chemical stratification and the rotation rate (through their slope) near the interface between the convective core and the radiative envelope in rapidly rotating intermediate-mass γ-Doradus (e.g. Van Reeth et al. 2015Ouazzani et al. 2017;Christophe et al. 2018;Li et al. 2019Li et al. , 2020Saio et al. 2021) and SPB stars (e.g. Pápics et al. 2015Pápics et al. , 2017Moravveji et al. 2016;Szewczuk & Daszyńska-Daszkiewicz 2018;Pedersen et al. 2021), thanks to the high precision of space-based photometric observations (e.g. Aerts et al. 2019;Aerts 2021, and references therein). This allows us to build intensive forward modelling of g-mode pulsating stars to constrain their internal structural and dynamical properties. This approximation, in its standard version, is only applicable for low-frequency GIWs propagating in strongly stably stratified zones of uniformly rotating spherical stars. In this case a set of assumptions should be verified: (i) the buoyancy force is stronger than the Coriolis acceleration (i.e. 2Ω N, where N is the Brunt-Väisälä frequency) in the direction of stable entropy or chemical stratification, (ii) the Brunt-Väisälä frequency is larger than the frequency of the waves in the rotating frame ω (ω N), (iii) the rotation is assumed to be uniform and (iv) the star is assumed to be spherical, in other words, the centrifugal deformation of the star is neglected (i.e. Ω Ω K ≡ GM/R 3 , where Ω 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).
Recently, the TAR has been examined in deformed stars. First,  have included the centrifugal acceleration for slightly deformed stars using a first-order perturbative approach. Then, this perturbative framework has been optimised for practical applications to one-dimensional stellar structure models by Henneco et al. (2021). Second, in Dhouib et al. (2021, hereafter Paper I), we have used a non-perturbative approach to include the centrifugal acceleration for strongly deformed stars and planets by studying the dynamics of low-frequency GIWs in a general spheroidal coordinate system defined by Bonazzola et al. (1998) which follows the shape of a deformed star. It is important to note that in these studies the rotation was assumed to be uniform. These semi-analytical studies demonstrated that centrifugal effects' signatures are potentially detectable. Moreover, their results are in good qualitative agreement with direct computations of gravito-inertial modes in centrifugally deformed stellar models using 2D oscillation codes that take into account the full Coriolis and centrifugal accelerations such as the Toulouse Oscillation Program (TOP) (Ballot et al. 2010(Ballot et al. , 2012 and the adiabatic code of oscillation including rotation (ACOR) code (Ouazzani et al. 2012(Ouazzani et al. , 2015(Ouazzani et al. , 2017. In the case where the centrifugal acceleration is treated as a perturbation, but where the full Coriolis acceleration was taken into account, we can refer to the Tohoku oscillation code (Lee & Saio 1987;Lee & Baraffe 1995). The chosen coordinate system is used in ACOR, TOP and in the Evolution STEllaire en Rotation (ESTER) code (Espinosa Lara & Rieutord 2013) that computes the structure and the stationary internal hydrodynamics of (differentially-)rotating early-type stars such as the observed gravito-inertial modes pulsators we are studying.
However, the solid-body rotation assumption has to be potentially abandoned along the evolution of real stars (and planets) where gradients of angular velocity can develop, both in the radial and in the latitudinal directions. They can be triggered by the redistribution of angular momentum induced by stellar winds (e.g. Zahn 1992), structural adjustments (Maeder & Meynet 2000;Decressin et al. 2009), and by transport processes (Mathis 2013). First, as shown in Charbonnel et al. (2013) and Gallet & Bouvier (2015), early main sequence low-mass stars are potentially subject to a strong differential rotation. Moreover, Aerts et al. (2019) and references therein (we refer the reader to those provided in the first paragraph for F, A, and B-type stars) pointed out that the radiative envelope of intermediate-mass stars are the seat of a weak differential rotation for those observed with a small difference of rotation between the stellar surface and the core. This is also predicted by numerical simulations and models of different transport mechanisms like internal gravity waves (Rogers et al. 2013;Rogers 2015) or meridional flows (Maeder & Meynet 2000;Decressin et al. 2009;Rieutord et al. 2016;Ouazzani et al. 2019). In this framework, Ogilvie & Lin (2004) and Mathis (2009) have generalised the TAR by including the effects of general differential rotation on low-frequency GIWs propagating in spherical stars (and planets). This new formalism has been successfully applied to gravito-inertial mode pulsators observed by the Kepler (Borucki et al. 2009) space telescope. Indeed, Van Reeth et al. (2018) have applied it to derive the variation of the asymptotic period spacing in the case of a weak radial differential rotation. They evaluated the sensitivity of GIWs period spacing to the effect of such a radial shear in spherical stars and they demonstrated that they can be potentially detected if observing several modes.
Since both effects of the centrifugal deformation and differential rotation are potentially observable, we generalise here these previous works, treating the case of a general differential rotation in deformed stars (and planets). We begin in Sect. 2 with the derivation of the system of linearised hydrodynamic equations for GIWs in spheroidal coordinates which follows the shape of a differentially rotating deformed body. We choose again the coordinate system that was introduced by Bonazzola et al. (1998) which is used in 2D numerical models of rotating stars and pulsation codes. In Sect. 3, we build the generalised TAR in this set of equations by adopting the adequate assumptions. In Sect. 4, we rewrite the oscillation equations in the form of a generalised Laplace tidal equation and deduce the asymptotic frequencies of low-frequency GIWs. In Sect. 5, we carry out a numerical exploration of the eigenvalues and Hough eigenfunctions of the generalised Laplace tidal equation within the domain of validity of the TAR. We model early-type stars with 2D ESTER models (Espinosa Lara & Rieutord 2013). In Sect. 6, we quantify the differential rotation and the centrifugal deformation combined effects on the asymptotic period spacing pattern and we discuss their detectability. Finally, we discuss and summarise our work and results in Sect. 8.

Spheroidal geometry
As in Paper I, here, we use the spheroidal coordinate system (ζ, θ, ϕ) proposed by Bonazzola et al. (1998), where ζ is the pseudo-radial coordinate, θ the colatitude and ϕ the azimuth. Following Espinosa Lara & Rieutord (2013), this new coordinate system, illustrated in Fig. 1, can be linked to the usual spherical one (r, θ, ϕ) via the following mapping in the ith subdomain D i∈ 0,n−1 ∈ [R i (θ), R i+1 (θ)] of the spheroidal domain D where R i∈ 0,n (θ) are series of functions, such that R n (θ) = R s (θ) is the outer boundary and R 0 (θ) = 0 is the centre. Additionally, η i = R i (θ = 0) are the polar radii of the interfaces between the subdomains, The functions A i∈ 1,n−1 (ξ) = −2ξ 3 + 3ξ 2 and A 0 (ξ) = −1.5ξ 5 + 2.5ξ 3 and the constants a i = 1 are chosen to satisfy the boundary conditions between the different subdomains. We recall also the spheroidal base (a ζ , a θ , a ϕ ) defined using the mapping (Eq. 1) (Rieutord et al. 2005;Lignières et al. 2006;Reese et al. 2006) where r ζ ≡ ∂ ζ r, r θ ≡ ∂ θ r, (e r , e θ , e ϕ ) is the usual spherical base, ε = 1−R pol /R eq is the flatness, R s (θ), R eq and R pol are the surface, the equatorial and the polar radii, respectively.

Linearised hydrodynamic equations in spheroidal coordinates
To treat the wave dynamics in differentially rotating, strongly deformed stars and planets, we derive the complete adiabatic inviscid system of equations in spheroidal coordinates. First, the linearised momentum equation for an inviscid fluid is written as where v ζ , v θ , and v ϕ are the contravariant components of the velocity field u = v i a i and Ω = Ω(ζ, θ)(cos θe r − sin θe θ ) is the differential angular velocity of the star. ρ, Φ and P are the fluid density, gravitational potential and pressure, respectively. Each of these scalar quantities has been expanded as where X 0 is the hydrostatic component of X and X the wave's associated linear fluctuation. Following Cowling (1941), we neglect the fluctuation of the gravitational potential.
Article number, page 3 of 13 A&A proofs: manuscript no. main Next, the linearised continuity equation is obtained Then, the linearised energy equation in the adiabatic limit is derived where Γ 1 = (∂ ln P 0 /∂ ln ρ 0 ) S (S being the macroscopic entropy) is the adiabatic exponent, g eff = −∇Φ 0 + 1 2 Ω 2 ∇(r 2 sin 2 θ) = ∇P 0 /ρ 0 is the background effective gravity which includes the centrifugal acceleration and N 2 the squared Brunt-Väisälä frequency (or buoyancy frequency) given by Finally, we expand the wave's velocity and fluctuations on Fourier series both in time and in azimuth where m is the azimutal order and ω in is the wave frequency in an inertial reference frame. In a differentially rotating region, the waves are Doppler-shifted due to the differential rotation so we can define the wave frequency ω in the rotating reference frame as

Generalised TAR
3.1. Approximation on the stratification profile: As in Paper I, in order to obtain a separable system of equations when applying the TAR, we assume here that the Brunt-Väisälä frequency N 2 depends mainly on the pseudo-radius ζ. This implies also that the hydrostatic pressure P 0 and the hydrostatic density ρ 0 depend mainly on ζ. We present the validity domain of this approximation in Sect. 5.1 using two-dimensional ESTER stellar models.

The TAR with centrifugal acceleration in differentially rotation stars
By assuming the following frequencies hierarchy imposed by the TAR: ω N and 2Ω N (we verify this hierarchy in Sect. 7) which leads to a mostly horizontal velocity equation because of the strong buoyancy restoring force action in the vertical direction (|u ζ | {|u θ |, |u ϕ |}), the approximation N 2 (ζ, θ) ≈ N 2 (ζ), the Cowling approximation (Φ = 0), and the anelastic approximation (acoustic waves are filtered out), we can rewrite accordingly the energy equation (Eq. 7) and the three components of the momentum equation (Eqs. 3,4 & 5). First, we rewrite the energy equation as where we neglect the term (1/Γ 1 )P /P 0 in Eq. (7) using the anelastic approximation and we simplify the expression of the squared Brunt-Väisälä frequency (Eq. 8) using the approximation 3.1. Then, we simplify the radial momentum equation where we neglect the term (P /ρ 2 0 )∂ ζ ρ 0 thanks to the anelastic approximation and the vertical component of the acceleration and of the Coriolis acceleration, which are dominated by the buoyancy term since we assume that ω N and 2Ω N within the TAR. Subsequently, the latitudinal component of the momentum equation (Eq. 4) reduces to where we neglect the terms involving the vertical wave velocity since |u ζ | {|u θ |, |u ϕ |} because of the strong stable stratification. Finally, the azimuthal component of the momentum equation (Eq. 5) simplifies, for the same reasons, into The coefficients A, B, C, and D are defined in Table 1 and W = P/ρ 0 is the normalised pressure. We thus obtain a system of equations which has the same mathematical form as in the case of: (i) spherically symmetric uniformly rotating (Bildsten et al. 1996;Lee & Saio 1997) and differentially rotating (Mathis 2009) stars, (ii) weakly deformed uniformly rotating ) stars and (iii) strongly deformed uniformly rotating (Paper I) stars. We thus still manage to partially decouple the vertical and horizontal components of the velocity. By solving the system formed by Eqs. (14) and (15) we can express u θ and u ϕ as a function of W as follows where E and F are defined in Table 1, x = cos θ is the reduced latitudinal coordinate and : Terms involved in the derivation of the generalised Laplace tidal equation in the general case of spheroidal geometry and in the particular case of spherical geometry with differential rotation.

Terms Spheroidal Spherical
is the spin parameter. The structure of the new equations that we obtain in the spheroidal differentially rotating case is similar to the one with uniform rotation (Paper I). But we can point out two major differences. First, the wave frequency ω and the spin parameter ν are no longer uniform but they now depend on the pseudo-radius ζ and the reduced latitudinal coordinate x. That is why it is more convenient to index the operators and the variables by ω in (independent of ζ and θ) instead of using the spin parameter ν as in Paper I. Second, a new coefficient F is involved in the derivation of the dynamics of GIWs that models the action of differential rotation.

Dynamics of low-frequency gravito-inertial waves
We now derive the generalised Laplace tidal equation for the normalised pressure W, which allows us to compute the frequencies and periods of low-frequency GIWs and to build the corresponding seismic diagnostics. Applying the anelastic approximation and the approximation N 2 (ζ, θ) ≈ N 2 (ζ) in the continuity equation (Eq. 6) simplifies it into

JWKB approximation
Under the assumption that ω N, each scalar field and each component of u can be expanded using the two-dimensional Jeffreys-Wentzel-Kramers-Brillouin JWKB approximation (Fröman & Fröman 1965). In this case, 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 slowly varying envelope: with j ≡ {ζ, θ, ϕ}, k is the index of a latitudinal eigenmode (cf. Sect. 4.4) and A ω in km is the amplitude of the wave. Substituting the expansion given in Eqs. (20) and (21) into Eqs. (13), (16) and (17), the final pseudo-radial, latitudinal and azimuthal components of the velocity are obtained: 4.2. Approximation: A(ζ, θ) ≈ A(ζ) As in Paper I, to be able to introduce the eigenvalues Λ ω in km which depends only on ζ and derive the generalised Laplace tidal equation, we have to do a partial separation between the pseudo-radial and latitudinal variables in the pseudo-radial component of the velocity (Eq. 22). So, we have to assume that the coefficient A depends mainly on the pseudo-radius ζ. The validity domain of this approximation is presented in Sect. 5.1 using two-dimensional ESTER stellar models. Unlike the uniformly rotating case (Paper I), this assumption alone is insufficient to partially separate the variables in Eq. (22) in the differentially rotating case), so it is mandatory to assume a third approximation.

Approximation
: The wave frequency in the rotating frame (Eq. 11) depends on ζ and θ. This dependency comes from the differential rotation Ω(ζ, θ). Therefore, in order to preform the partial separation of variables in Eq. (22), the angular velocity Ω should have a dependency only on ζ. We examine in detail the validity of this approximation for two-dimensional ESTER stellar models in Sect. 5.1.1.

Generalised Laplace tidal equation
Substituting the expansion given in Eqs. (20) and (21) into Eqs. (13), (16) and (17), then replacing them into the continuity equation (Eq. 19) we obtain the generalised Laplace tidal operator (GLTO) thus the generalised Laplace tidal equation (GLTE) for the normalised pressure w ω in km is written as where Λ ω in km are the eigenvalues deduced from the following dispersion relation: and w ω in km the generalised Hough functions (eigenfunctions) of the GLTE. We choose to define our latitudinal ordering number k to enumerate, for each (ω in , m), the infinite set of solutions as in Lee & Saio (1997), Mathis (2009) and .
The GLTO (Eq. 25) reduces to the Laplace tidal operator in the case of a uniformly rotating deformed star (Paper I). In this case, ω and ν are constants so we can take them out of the derivatives and F = 0; thus we obtain . (28) Furthermore, the GLTO (Eq. 25) reduces to the Laplace tidal operator derived by Mathis (2009) in the case of a spherical differentially rotating star. In this case, we can replace coefficients A, B, C, D, and E by their analytical expressions in the particular case of a spherical geometry presented in Table 1 of Paper I, so the new coefficient F defined here can be written in this particular case as thus we obtain with

Asymptotic frequency and period spacing of low-frequency GIWs
By substituting the dispersion relation (Eq. 27) into the following quantisation condition in the vertical direction defined by Unno et al. (1989), Gough (1993), and Christensen-Dalsgaard (1997) and used in Paper I ζ 2 ζ 1 k V;ω in nkm dζ = (n + 1/2)π, where ζ 1 and ζ 2 are the turning points of the Brunt-Väisälä frequency N and n is the radial order, we compute numerically the asymptotic frequencies of low-frequency GIWs in deformed differentially rotating stars in the inertial frame. We describe in detail the used numerical method in Sec. 6. In the case of a weak differential rotation, we can extract, analytically, the expression of the asymptotic frequencies in the rotating frame with ω in n = ω in nkm = ω nkm + mΩ av , where Ω av is here the averaged rotation. Then by applying a firstorder Taylor development following Bouabid et al. (2013), we can generalise his expression of the corresponding period spacing obtained in the spherical case

Application to rapidly and differentially rotating early-type stars
As in Paper I, we use here also ESTER models (Espinosa Lara & Rieutord 2013) to implement our equations. We use 3M stellar models with a hydrogen mass fraction in the core X c = 0.7 rotating at the equator at [0%, 90%] of the Keplerian break-up rotation rate.

Domain of validity of the generalised TAR
We study, within the framework of ESTER models, the validity domain of the three approximations (3.1, 4.2 and 4.3) that are necessary to build the TAR in the case of rapidly differentially rotating deformed bodies: The validity of the approximations (37) and (38) is well discussed in Paper I (the N 2 and A profiles are the same in the uniformly and differentially rotating cases). Here, we focus on determining the validity domain of the approximation (39) by applying the same method used to evaluate the other two assumptions that we recall here.

Validity of Ω(ζ, θ) ≈ Ω(ζ)
To visualise the problem, first we illustrate in Fig. 2 the function Ω(ζ, θ) computed with an ESTER model (3 M , X c = 0.7) rotating at the equator at Ω/Ω K = 20%. We can see that the angular velocity is equal to 14.66 µHz at the core, and at the surface we get 11.41 µHz at the pole and 12.22 µHz at the equator. Then, to evaluate the committed error by adopting the approximation we represent in Fig. 3 the angular velocity profile as a function of the pseudo-radius ζ for different values of the colatitude θ f , at a fixed rotation rate at the equator Ω = 0.2Ω K . Visually, we can notice that the latitudinal gradient of the angular velocity is less important than the radial gradient. In order to confirm this observation, we show, in the bottom panels, the corresponding relative error ∆Ω θ f (ζ) between the exact value Ω(ζ, θ f ) given by the model and the approximated model Ω approx (ζ) (the weighted average of the exact value over the colatitude θ,Ω(ζ)). We define these quantities as follows: with where θ f is a fixed value of the colatitude θ and we used the symmetry propriety of Ω with respect to the equator. We can conclude that this approximation is valid with a level of confidence of 90% for all pseudo-radii ζ and for all rotation rates Ω/Ω K ∈ ]0%, 90%] since the error rate is always lower than the maximal error rate fixed to 10%. We also recall here that within the physics treated in the ESTER code the differential rotation results from the baroclinic meridional circulation (Zahn 1992;Maeder & Zahn 1998;Rieutord 2006) and an isotropic viscous transport. A weaker dependence of Ω on θ can be enforced by horizontal turbulent transport (Zahn 1992;Maeder 2003;Mathis et al. , 2018Park et al. 2021) or magnetic fields (Mestel & Weiss 1987).

Validity domain of the TAR
We determine the validity domain of the three approximations (37), (38), and (39) as a function of the rotation rate and the pseudo-radius by varying systematically the normalised rotation rate of the star Ω/Ω K and by calculating for each case the associated maximum relative error. Afterwards, we fix a maximum error rate equal to 10% and we deduce the pseudo-radius limit ζ limit where the maximum relative error exceeds this threshold and we decide that the approximations become invalid. Physically, it means that the variation of the structural properties and rotation profile are weak in the latitudinal direction and vary mostly with the pseudo-radius that allows us to build the TAR. The best method in the future to improve the control of the adopted assumptions would be to compute asymptotic frequencies using the latitudinally averaged quantities within the TAR and to compare them to frequencies computed using 2D oscillation codes Fig. 3: Profile of the angular velocity Ω and the relative error δΩ θ f of the approximation Ω(ζ, θ) ≈ Ω(ζ) as a function of ζ at different colatitudes θ f using an ESTER model rotating at the equator at 20% (above) and 90% (below) of the Keplerian breakup rotation rate (the light orange area indicate the margin of error which we allow).
applied to deformed models (e.g. Ouazzani et al. 2015Ouazzani et al. , 2017Reese et al. 2021). In Fig. 4, we display the pseudo-radius limit ζ limit as a function of the rotation rate Ω/Ω K for 3 solar masses ESTER models with X c = 0.7 and X c = 0.2. This curve defines the validity domains of the approximations (37), (38), and (39) which consequently define the validity domain of the TAR. We note that the influence of the hydrogen mass fraction in the core X c on the validity domain of the TAR is very weak. It is therefore possible to apply the generalised TAR to differentially rotating early-type stars rotating up to 20% of the Keplerian critical angular velocity. This limit is the same as the one found in Paper I, so the differential rotation have no influence on the validity domain of the TAR. Chemical gradients created near the core along the evolution are not treated in 2D ESTER models yet. Their treatment can be foreseen by adding them 'by hand' to explore their effects as this has been done for instance for acoustic glitches by Reese et al. (2014Reese et al. ( , 2021.

Eigenvalues and Hough functions
We solve the GLTE for different pseudo-radii, wave frequency in the inertial frame, and rotation rates within the scope of the defined validity domain using an implementation based on Chebyshev polynomials (Wang et al. 2016) as in Paper I. Here, we are no longer able to calculate the spectrum of the GLTE as a function of the spin parameter ν because in the differentially rotating case, ν is no longer a constant value but it becomes a function of ζ and θ. So to make the comparison between the new results and our previous results clear, we define the weighted average of the spin parameter Fig. 5 shows the eigenvalues Λ ω in km as a function of ω in for m = 1 and ζ = ζ limit = 1 at Ω/Ω K = 0.2. Since Ω points in the direction of θ = 0 and the oscillations are proportional to e i(ωt−mϕ) , a prograde (retrograde) oscillation correspond to a positive (negative) value of the product mν which depends on the spatial coordinates. This figure reveals that the differential rotation of the star causes at ζ = 1 a modest gradual shift in the eigenvalues. We Fig. 5: Spectrum of the GLTE as a function of the wave frequency in the inertial frame ω in at ζ = 1 and m = 1 for Ω = 0.2Ω K in the case of a uniformly rotating star (above) and of a differentially rotating star (below). Blue (respectively, orange) dots correspond to even (respectively, odd) eigenfunctions. The vertical lines correspond to the mean values of the spin parameterν.
recall here the two main families of eigenvalues. The first one being gravity-like solutions (Λ ω in km 0). They exist for any value of ω in 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. They exist only for specific values of wave frequency in the inertial reference frame (ω in ∈ [ω in − , ω in + ]) and we attach negative k's to them. They appear only in rotating stars and they correspond to Rossby modes (r modes) if they are retrograde and have positive eigenvalues (Λ ω in km > 0) and to overstable convective modes if they are prograde and have negative eigenvalues (Λ ω in km < 0).
Contrary to the uniformly rotating case we obtain here eigenvalues which depend considerably on the pseudo-radius ζ as shown in Fig. 6 where we represent the spectrum of the GLTE as a function of the pseudo-radius at Ω/Ω K = 0.2, m = 1 and ω in = 1.12 × 10 −4 rad/s. More precisely we can see that the differential rotation influence especially the region close to the convective core (0.153 ≤ ζ ≤ 0.35) where the rotation is maximal and the radial gradient of the angular velocity is the highest (cf. Fig. 2). We can see that this dependency becomes very low in the external region 0.5 ≤ ζ ≤ 1 where we obtain eigenvalues very close to the ones in the uniformly rotating case.
We focus now on the influence of the differential rotation on the generalised Hough functions w km which varies with the pseudo-radius ζ and the horizontal coordinate x. This dependence is illustrated in Fig. 7 at Ω/Ω K = 20% for prograde dipole {k = 0, m = 1} modes with ω in = 10 −4 rad/s and retrograde Rossby {k = −2, m = −1} modes with ω in = −5.5×10 −5 rad/s for four different cases: Fig. 7a represents our starting point where Fig. 6: Spectrum of the GLTE as a function of the pseudo-radius for m = 1 and ω in = 1.12 × 10 −4 rad/s at Ω/Ω K = 0.2 in the case of a uniformly rotating star (above) and of a differentially rotating star (below). Blue (respectively, orange) dots correspond to even (respectively, odd) Hough functions. the rotation is uniform (Paper I); in Fig. 7b, we take into account the pseudo-radial differential rotation; in Fig. 7c, we take into account the latitudinal differential rotation and Fig. 7d represents our final point where we take into account the full differential rotation profile. As the pseudo-radius decreases from the surface (ζ = 1) to the edge of the radiative zone (ζ = 0.153), Rossby-like solutions are slightly modified, whereas gravity-like solutions considerably change. We can also see that the major modification comes from the differential rotation gradient along the pseudoradius ζ (Fig. 7b) while the differential rotation according to the colatitude θ modifies marginally the gravity-like solutions and has almost no impact on the Rossby-like solution (Fig. 7c). On top of that, we can see that the Hough functions at the surface are not too affected in all the cases and that the strongest variation occur in the interior of the star where the pseudo-radial gradient of the rotation is high.
Overall, we can see clearly in Fig. 7 that the differential rotation strongly influences and modifies the eigenfunctions of the GLTE by introducing a new dependency on the pseudo-radial coordinate ζ, whereas in the spherically symmetric uniformly rotating case they were only dependant on the latitudinal coordinate θ while they were also slightly dependent on ζ in the uniformly rotating deformed case. In fact, the gravity-like 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 that the gravity-like solutions migrate onwards, away from the equator (x = 0), causing a broadening of its shape. This corresponds to the modification of the region of propagation of gravito-inertial waves by the differential rotation (Mathis 2009;Mirouh et al. 2016;Prat et al. 2018). Indeed, since the angular velocity Ω(ζ) decreases with the distance to the centre, the equatorial trapping of sub-inertial (ν > 1) gravito-inertial modes become less and less important.

Asymptotic period spacing pattern
In order to compute the period spacing patterns, we adapt the method developed in Henneco et al. (2021) and used in Paper I. First, we calculate for each radial order n of a given mode (k, m), the corresponding asymptotic frequencies in the inertial frame ω in nkm using Eqs. (11) and (33). Then, we can calculate the asymptotic periods in the inertial frame (P in nkm = 2π/ω in nkm ) and the corresponding period spacing which is shown in Fig. 8 for {k = 0, m = 1}, {k = 1, m = 0}, and {k = 0, m = −1} modes. The periods represented here are calculated for radial orders between n = 5 and n = 50 but the Cowling approximation that we adopted here to develop our formalism is valid only for high radial orders (Bouabid et al. 2013;Ouazzani et al. 2017). That is why we hatched the region with low radial order (n < 20), so we can identify where the Cowling approximation is potentially invalid. We find a net decrease in the period spacing ∆P in nkm caused by the differential rotation. The centrifugal and the differential rotation effects are moderately significant when assessed within the validity domain of the generalised TAR nonetheless the global characteristics of the period spacing pattern are conserved. This is in good agreement with the results obtained by Ouazzani et al. (2017). They computed gravito-inertial modes and their period spacing (represented in Fig. 6 of their article) using the ACOR 2D oscillation code with deformed stellar models (see Ouazzani et al. (2015) for details on their method), which take the centrifugal acceleration into account following the method of deformation of acoustic models by Roxburgh (2006). Ballot et al. (2012) also studied the effect of the centrifugal acceleration on the period spacing. They found, for a sectoral mode (k = 0), a discrepancy between the period spacing (represented in Fig. 2 of their article) computed using the spherical TAR and the complete computations using TOP with spheroidal models. They demonstrate that this discrepancy originates in the centrifugal distortion of the 2D models. But for other modes, they do not see a considerable effect of the centrifugal deformation on the period spacing pattern in agreement with our results.
Even though the centrifugal deformation effects predicted by the TAR are weak in intermediate-mass stars, they are theoretically detectable for some radial orders of several modes using observations from TESS (Transiting Exoplanet Survey Satellite; Ricker et al. 2014) and Kepler (Paper I). This is why we do not neglect them and we study the differential rotation and the deformation effects simultaneously. To quantify the induced variation and evaluate its detectability, we compute the frequency differences.

Detectability of the differential rotation effect
To evaluate the detectability of the differential rotation effect with space-based photometric observations, we compute first the frequency differences ∆ f centrifugal between asymptotic frequencies calculated in the standard TAR (TAR s ) and those calculated in the generalised TAR with a uniform rotation (TAR g u ) and a differential rotation (TAR g d ) through Eq. (33): with x and y two parameters that allow us to choose the effect to be evaluated: (g u , s) for the effect of the centrifugal acceleration; (g u , g d ) for the effect of the differential rotation; (g d , s) for the two effects simultaneously.
Then, by comparing the obtained frequency differences with the frequency resolutions ( f res = 1/T obs ) of Kepler and TESS light curves covering quasi-continuously observation times of T obs = 4 years and T obs = 351 days, respectively, we can deduce the radial orders n's for which the frequency differences are expected to be theoretically detectable. In Fig. 9, we display these results for {k = 0, m = 1}, {k = 1, m = 0}, and {k = 0, m = −1} modes rotating at 0.2 Ω K . In this figure, we show the detectability of the centrifugal effect, the differential rotation effect and of the two effects taken simultaneously into account. Unlike the centrifugal effect, we can see clearly that the effect of the differential rotation for these modes is, theoretically, largely detectable for all radial order higher than 15 using TESS and Kepler observations (not represented here because its observation time is larger than the one of TESS so the detectability using Kepler is guaranteed). The detectability of the centrifugal effect is discussed in details in Paper I.
As we pointed out in Paper I, this is a formal evaluation due to the presence of large correlations between the different parameters of our stellar models. This will mask the effect of the centrifugal acceleration and the differential rotation in forward asteroseismic modelling analyses, even when it is theoretically detectable:

Evaluation of the terms hierarchy imposed by the TAR within differentially rotating deformed stars
The hierarchy of terms imposed by the TAR can be summarised by the following frequency hierarchy: which ensures the other hierarchies. Using the Brunt-Väisälä frequency and rotation profiles from ESTER models and using the asymptotic frequencies calculated in Sect. 6, we compare these frequencies and we discuss whether the TAR is still valid in differentially rotating deformed stars as in Paper I.
7.1. The strong stratification assumption: 2Ω N We evaluate the term N/2Ω using 3 M , X c = 0.7 ESTER models for rotation rates Ω/Ω K ∈ [0.1, 0.9]. As shown in Fig. 10, the strong stratification assumption is valid only in the radiative zone away from the border between the convective core and the radiative envelope (ζ ≥ 0.2) at Ω/Ω K ≤ 0.2. Beyond this critical value, the Brunt-Väisälä frequency and the frequency of rotation have a close order of magnitude so the strong stratification approximation is no longer valid.
7.2. The low frequency assumption: ω N We evaluate the term N/ω using the asymptotic frequencies for {k = 0, m = 1}, {k = 1, m = 0}, and {k = 0, m = −1} modes calculated in Sect. 6. As shown in Fig. 11, the low frequency assumption is valid in all the space domain for the {k = 0, m = 1} mode only for high radial orders (n > 20). Below this critical value, the Brunt-Väisälä frequency and the wave frequency can have a very close order of magnitude near the transition layer between the convective core and the radiative envelope. But far from this interface, the low frequency approximation is valid for all the modes.

Discussion and conclusions
This work is a continuation of our previous study where we generalised the TAR in the case of strongly deformed, rapidly and uniformly rotating stars and planets. In this perspective, we study the possibility of carrying out a new generalisation of the TAR that abandons the assumption of uniform rotation and takes into account the radial and latitudinal differential rotation. We approached this exploration by deriving the generalised Laplace Tidal equation in a spheroidal coordinate system. We relied mainly on two assumptions that will define the validity domain of the generalised TAR (Paper I). The equation that we derive has a similar form to the one that is 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), moderately rotating weakly deformed stars Henneco et al. 2021) and uniformly rotating strongly deformed stars (Paper I). So with this new formalism, we can study GIWs in the radiative region of all types of stars and planets. We apply this general formalism to rapidly and differentially rotating early-type stars using 2D ESTER models (M = Fig. 9: Detectable radial orders n for {k = 0, m = 1} (top) {k = 1, m = 0} (middle), {k = 0, m = −1} (bottom) modes at a rotation rate Ω/Ω K = 20% based on the frequency resolution of TESS. Blue, yellow and red dots represent respectively the centrifugal effect, the differential rotation effect and the sum of the two effects (the hatched pink area indicates where the Cowling approximation is potentially invalid) .
3M , X c = 0.7) where we found that the signature of the differential rotation effect in the period spacing patterns is stronger than the signature of the centrifugal effect for the {k = 0, m = 1} mode at Ω/Ω K = 20%, typically up to a factor twenty. This is caused essentially by the presence of a strong pseudo-radial Ωgradients in stellar interior of early-type stars. This new generalisation of the TAR can be used to study the dissipation of stellar and planetary tides induced by low-frequency GIWs in rapidly and differentially rotating deformed stars and planets and the angular momentum transport with a formalism that can be directly implemented in ESTER models. Fig. 10: N/2Ω term as a function of the pseudo-radius ζ at different rotation rates Ω/Ω K (the hatched area represents the convective region of the star).
The next step will be the inclusion of the magnetic field in a non-perturbative way within the generalised TAR formalism. So far, Prat et al. (2019Prat et al. ( , 2020 and Van Beeck et al. (2020) have focused on the case where magnetic fields are weak enough to be treated within a perturbative treatment to study the effects of a magnetic field on the seismic parameters of g modes which become magneto-gravito inertial modes (Mathis & de Brye 2011. So a possible follow-up of this work (Paper III) is to take into account the stellar magnetic fields and generalise the TAR framework to the case of differentially rotating strongly deformed magnetic stars.