Influence of stellar structure, evolution and rotation on the tidal damping of exoplanetary spin-orbit angles
Cilia Damiani Stéphane Mathis
 
   
   
   Université Paris-Sud, CNRS, Institut d'Astrophysique Spatiale, UMR8617, 91405 Orsay Cedex, France
            Max Planck Institute for Solar System Research, Justus-von-Liebig-Weg 3, 37077 Göttingen, Germany ;   
          damiani@mps.mpg.de
          IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France
          Université Paris Diderot, AIM, Sorbonne Paris Cité, CEA, CNRS, F-91191 Gif-sur-Yvette, France
             
             
Influence of stellar evolution on the tidal damping of obliquities

   
  
  It is debated whether close-in giant planets can form in-situ and if not, which mechanisms are responsible for their migration. One of the observable tests for migration theories is the current value of the obliquity, i.e. the angle between the stellar equatorial plane and the orbital plane. But after the main migration mechanism has ended, the obliquity and the semi-major axis keep on evolving due to the combined effects of tides and magnetic braking. The observed correlation between effective temperature and measured projected obliquity in well characterised systems has been taken as evidence of such mechanisms being at play.
  
 Our aim is to produce an improved model for the tidal evolution of the obliquity, including all the components of the dynamical tide for circular misaligned systems. This model is developed to take into account the strong variations in structure and rotation of stars during their evolution, and their consequences for the efficiency of tidal dissipation.
   
Our model uses an analytical formulation for the frequency-averaged dissipation in convective layers for each mode, depending only on global stellar parameters and rotation. It also includes the effect of magnetic braking in the framework of a double zone stellar model.
  
For typical hot-Jupiters orbital configurations, the obliquity is generally damped on a much shorter timescale than the semi-major axis. The final outcome of tidal evolution is also very sensitive to the initial conditions, with Jupiter-mass planets being either quickly destroyed or brought on more distant orbit, depending on the initial ratio of planetary orbital momentum to stellar spin momentum. However we find that everything else being the same, the evolution of the obliquity around low-mass stars with a thin convective zone is not slower than around those with a thicker convective zone. On the contrary, we find that more massive stars, remaining faster rotator throughout their main-sequence, produce more efficient dissipation.
  
 

 

   Influence of stellar structure, evolution and rotation 
 on the tidal damping of exoplanetary spin-orbit angles
    C. Damiani
          1,2
          S. Mathis 3,4
          
    Received xxx, 2017; accepted xxx
===============================================================================================================






§ INTRODUCTION


	Over the last twenty years, tremendous observational efforts have resulted in the complete characterisation and determination of orbital properties of a few hundreds of exoplanetary systems <cit.>. For about a hundred of them, this includes measurements of the Rossiter-MacLaughlin effect <cit.>, a signature in radial velocities produced during a transit that contains information about the obliquity Θ of the system, i.e. the angle between the stellar spin axis and the orbital normal. Most of the time, what is observed is the projection on the sky of this angle, usually called λ. The constraint obtained on Θ through the observation of λ has a complex nature, but when λ is observed to be large, then it is more likely that Θ is also large. Moreover, when the orbit is seen edge-on, the value of λ can distinguish wether the orbit is prograde or retrograde <cit.>. Those measurements are especially valuable, because the obliquity of exoplanetary systems is one of the main observables allowing to put constraints on their formation and evolution.

	In the case of close-in planets however, which constitute the vast majority of systems with known obliquities, the orbital elements keep evolving after the end of the main migration mechanism under the influence of tides. According to <cit.> and <cit.>, observations seem to show that hot stars with hot Jupiters have high obliquities. They interpret this as evidence of tidal evolution, acting upon a population of initially randomly distributed obliquities. Indeed, tidal dissipation is expected to be more efficient in stars having important convective envelopes, because of the turbulent friction acting on the tidal flows <cit.>. Thus cooler star's spin would realign with the orbit more rapidly. Because the mass of the outer convection zone on the main sequence decreases rapidly with increasing stellar mass after 1.2 M⊙, which corresponds to an effective temperature of T_ eff≈ 6250 K, this would naturally explain the sharp increase in λ that occurs around this effective temperature. This interpretation would favour a single migration mechanism, capable of producing randomly distributed obliquities.

	There is a major flaw in this reasoning: if tides are efficient at damping the obliquity of hot-Jupiters, they would also lead to significant orbital decay, and the planet would be destroyed before re-alignment. This result is obtained with the classical theory of the hydrostatic equilibrium tides <cit.>, an approximation that seems to be incompatible with other orbital properties of hot Jupiters <cit.>. Corrections to the equilibrium tide are provided by the dynamical tide, i.e. the response of a rotating fluid to the harmonic forcing of a time-varying perturbing potential. This includes excitation and damping of inertial waves due to the Coriolis acceleration in convective regions and excitation and damping of gravity waves due to buoyancy modified by rotation (then called gravito-inertial waves) in radiative regions <cit.>. 

	In the framework of the dynamical tide, <cit.> suggested that tidal dissipation efficiency may not be the same for different processes, e.g. spin-orbit alignment and orbital decay. <cit.>  pointed out that in the case of a misaligned system, the tidal potential has a finite number of components, corresponding to all combinations of the order m and degree l (-l≤ m ≤ l) of the spherical harmonics decomposition, and that in general each of the distinct components will generate a tidal perturbation with its own dissipation efficiency. The efficiency of tidal dissipation is usually parameterised by a dimensionless quality factor Q proportional to the ratio of the total kinetic energy of the tidal distortion to the energy dissipated in one tidal period <cit.>. In the theory, for the quadrupolar term, Q always appears in the combination Q' ≡ (3/2)(Q/k_2), where k_2 is the static Love number of the star that measures its density stratification. Therefore, Q' is the most widely used parameter to estimate tidal dissipation efficiency. The lower the value of Q' the stronger the tidal dissipation. In general, Q' depends on l, m, the tidal frequency (a linear combination of stellar rotation frequency and orbital frequency), and the structure of the body under consideration. Thus a rigorous treatment of the tidal dissipation should consider the sum of the effects associated with the different tidal components each having its specific Q' <cit.>.  

	Using the tidal prescription proposed by <cit.>, <cit.> have computed the coupled evolution of the orbital elements and stellar spin of five representative systems, taking into account the combined effects of tides, stellar wind, mass loss, magnetic braking and stellar evolution. Their results show that, accounting for all the relevant physical mechanisms, the current properties of the systems they consider can indeed be naturally explained. However, an important limitation of this study is the absence of a reliable estimation for the tidal quality factor Q'_ 10 (in the notations introduced by <cit.>) of the obliquity tide (corresponding to the component associated with a torque but with no dissipated power). In their study, Q'_ 10 is a free parameter which is allowed to vary within a broad range of values, instead of being imposed by stellar properties, to reproduce the observations. While this approach validates 's basic idea, it is not enough to accurately constrain the initial distribution of the obliquity of the population of known hot Jupiters.
	
		
	Another fundamental limitation of their approach, inherited from 's formulation, is that they only consider one of the components of the dynamical tide to have a possibly different tidal quality factor than the others. Indeed, this parametrisation may be reasonable for the current observed orbital state of known hot-Jupiters, for which the orbital frequency is usually much greater than the stellar rotation frequency. Then in this case, there is only one mode that is capable of exciting inertial waves in the star. However, when the orbital frequency is shorter than the rotation frequency, which can in general be the case in early stages, other components of the tide may fall in the frequency range of inertial waves, and the resulting dissipation can also contribute to the evolution of the semi-major axis. As already shown by other studies, this clearly stresses the need to take into account the rotational history of the star to compute the tidal evolution of the system <cit.>.
	
	In this paper, we propose a tidal evolution model that follows Lai's theory and take into account the dynamical tide. We expand on this and <cit.> works by including all the seven components of the quadrupolar term of the tidal potential into the orbital evolution. Moreover, in a similar approach as <cit.>, we also include an analytical formulation of the tidal dissipation efficiency associated with each mode, as a function of stellar structure and rotation, by using the frequency-averaged formulation for dissipation in convective zones derived by <cit.>. In a series of paper <cit.>, this approach has already shown the importance of the evolution of stellar structure and rotation, for the frequency-averaged dissipation efficiency. This compels the use of a double-zone model, one zone representing the radiative core and the other the convective envelope, accounting for the evolution of stellar structure in time. We also use a double-zone semi-empirical model for magnetic braking <cit.>, to account for the evolution of stellar rotation and the internal redistribution of angular momentum within the star.
	

The paper is organised as follows. In Section <ref>, we formulate the double-zone tidal model for the evolution of the obliquity in exoplanetary systems. We also describe and motivate our approximations. In Section <ref>, we compute the frequency averaged-dissipation for different modes involved in the evolution of the semi-major axis and obliquity, and we discuss its dependence  with mass and rotation. In Section <ref>, we compute the coupled evolution of semi-major axis, stellar rotation and obliquity for a sample of initially misaligned systems and discuss the different characteristic time-scales. In Section <ref>, we discuss the main results and perspectives of this work.
	




§ THE DOUBLE-ZONE WEAK FRICTION MODEL FOR TIDAL INTERACTION

We consider a star of mass M_⋆ orbited by a planet of mass M_ p.  The star is modelled as a deformable core of homogenous density ρ_ c and mean radius R_ c in uniform rotation, surrounded by an envelope of homogeneous lesser density ρ_ e and mean radius R_⋆ also rotating uniformly at the angular frequency Ω_ e, but not necessarily the same as that of the core Ω_ c. We only consider the tides raised by the planet on the star. As shown by other authors <cit.>, for typical hot-Jupiters, tides raised on the planet result in the evolution of the planetary spin magnitude and direction towards their equilibrium value within a time scale of about 10^5 yr. This time is sufficiently short compared to the typical stellar evolution timescales to safely neglect the effects of tidal dissipation in the planet on the overall evolution in the case of a circular orbit. We also neglect the mutual interaction of the tidal bulges and consider that the tide-generating potential of the planet is that of a mass-point. This hypothesis should be relaxed only if the semi-major axis a < 5 R_⋆ and if the gravitational quadrupolar moment J_2 > 10^-2 for the planet <cit.>. The planet generates a time-varying tidal potential Ψ that changes the shape and as a result, the exterior potential Φ of the star. In polar coordinates, the time-varying tidal potential produced by the point-like planet can be expanded in terms of spherical harmonics Y_l^m of degree l and order m. For an eccentric orbit, Ψ must be expanded for an infinite number of terms. But for a circular orbit, not necessarily in the equatorial plane of the reference frame, there is a finite number of tidal components for any given value of l. Here we are interested in the evolution of the obliquity of the orbit resulting from the dissipation of the tides in the star, so we limit our study to the case of circular orbits, and produce a tidal evolution model that is valid for arbitrary values of the obliquity.



 §.§ Dynamical tidal potential and response in the circular case

 We use spherical polar coordinates (r, θ, ϕ) , where r is the radial distance from the center of the star, θ is the colatitude, and ϕ is the azimuthal angle. We adopt an inertial coordinate system with polar axis aligned with the star's angular momentum L_⋆. We note Θ, the angle between the orbital angular momentum L_ o and L_⋆, a, the semi-major axis of the relative orbit of the planet and Ω_0, its mean orbital angular velocity. We chose the y-axis in the direction L_⋆×L_ o. Neglecting precession, the tide-raising potential of the planet can be written as

    Ψ(r, θ,ϕ, t) = Re[∑_l = 2^∞∑_m=-l^l∑_p=0^lG M_ p R_⋆^l/a^l+1. 
    . W_l^l-2p d^(l)_m, l-2p(Θ) e^-i(l-2p) Ω_0 t(r/R_⋆)^l Y_l^m(θ, ϕ) ],

where G is the gravitational constant, d^(l)_m, l-2p(Θ) are the Wigner d-matrix elements that only depend on Θ, and 

    W_l^l-2p = 4π/2l+1 Y_l^l-2p(π/2, 0).

Each tidal component rotates with the angular velocity ω given by

    ω = (l-2p)Ω_0,

where p is an integer with 0≤ p≤ l that arises from the inclination of the orbital plane. This defines the tidal frequency in the inertial frame. In a frame that rotates with the spin angular velocity Ω_ e, the tidal frequency must be Doppler shifted and we thus define ω̂, the tidal frequency in the fluid frame

    ω̂ = (l-2p)Ω_0 - m Ω_ e.

While the tidal frequency can be positive or negative, the physical forcing frequency is positive and is thus simply defined as χ = |ω̂|. Since we have ω̂_l, m, p=-ω̂_l, -m, l-p, the (m, p) component is physically identical to the  (-m, l-p) component.

When the amplitude of the tidal disturbance is small, as this is generally the case for close-in exoplanetary systems, and for axisymmetric bodies, the tidal response of the perturbed body can be determined by treating each component of the tidal potential independently, and considering that the total response is simply the sum of each component. Taking the Fourier transform of the equations of motion, mass conservation and gravitational potential leads to a purely spatial problem in which the frequency ω appears as a parameter. 

Considering a tidal potential of the form Ψ=∑Ψ_l^m(t) (r/R_⋆)^l Y_l^m (θ, ϕ), we write the perturbed stellar gravitational potential in the form Φ=∑Φ_l^m(t) (r/R_⋆)^-(l+1)Y_l^m(θ, ϕ). We denote the temporal Fourier transforms with a tilde and we associate to each component of the perturbed gravitational potential an array of potential Love numbers k_l, l'^m, m'(ω), that quantify the frequency-dependent response of the star to tidal forcing : 

    Φ̃_l^m(ω)= ∑_l', m'k_l, l'^m, m'(ω) Ψ̃^l'_m'(ω).

For an axisymmetric body, k_l, l'^m, m'(ω)=0 unless m'=m. The part of the response that has the same form as the applied potential is thus quantified by the potential Love numbers k_l^m(ω) :

    k_l^m(ω) = k_l, l^m, m(ω).




 §.§ Tidal dissipation
 
 Stars are fluid bodies that cannot be considered as purely frictionless. Dissipation of mechanical energy occurs even in adiabatically stratified fluid layers through turbulent friction or partial wave reflection at the interface with stably stratified zones for example. The fluid dynamics of tidally forced bodies is a very active field of research <cit.>, but whatever the mechanism actually responsible for dissipation, it suffices to say here that the response of a fluid body to tidal excitation is not frictionless and some energy is transferred to the perturbed body. In the presence of dissipation, the perturbed external gravitational potential of the star involves complex Love numbers and the imaginary part of the Love numbers Im[k_l^m(ω)] quantifies the component of the response that is out of phase with the tidal forcing, and is associated with transfers of energy and angular momentum. 

The rates of transfer of energy and of the axial component of angular momentum from the orbit to the body, measured in an inertial frame, define the tidal power P, and the tidal torque T, respectively. So for each component of the tidal potential of the form of Eq. (<ref>), averaged over azimuth in the case m=0, it can be shown <cit.> that P=ω𝒯 and T=m𝒯, where 𝒯 depends on the amplitude of the tidal component and on Im[k_l^m(ω)]. To compute the evolution of orbital elements associated with these transfers of energy and angular momentum, we use an approximation similar to the one made in the equilibrium theory of tides. We consider that the tidal response resembles the hydrostatic one in the absence of dissipation, and that dissipation introduces a small frequency-dependent phase lag in the response of the body. However, contrary to the usual equilibrium theory of tides, we do not impose that all components have the same phase-lag. So for each component of the tidal potential, we assume that dissipation introduces separate frequency-dependent phase-lags Δ_l^m(ω) which can be related to the imaginary part of the Love numbers <cit.> with 

    Im[k_l^m(ω)] =  |k_l^m(ω)| sinΔ_l^m(ω)

where Im[k_l^m(ω)] has the same sign as ω̂.
As detailed in <cit.>, typically Re[k_l^m(ω)] is a quantity of order unity, only weakly dependent on m and ω, and can be well approximated by its hydrostatic value. The hydrostatic Love number k_l is a real quantity and does not depend on m. Its evaluation for a fluid body is a classical problem involving the Clairaut's equation (see Appendix. <ref>). In the weak friction approximation <cit.>, we thus consider that |Δ_l^m(ω)| << 1 so that 

    |k_l^m(ω)| ≈ k_l

and 

    Im[k_l^m(ω)] ≈  k_l Δ_l^m(ω).


In general, the tidal response involves resonances with stellar oscillations when the tidal frequency matches that of an appropriate mode. In convection zones, for solar-like slow rotating stars, the frequencies of tidal oscillations typically lie in the range of inertial modes, for which the Coriolis acceleration provides the restoring force. The eigenfrequencies of inertial modes in the frame that rotates with angular velocity Ω_ e are dense in the interval [-2Ω_ e, 2Ω_ e], but do not exist outside of this range. The spatial structure of the inertial modes imply that any spherical harmonic potential resonates with only a finite number of inertial modes. As shown for example by <cit.>, enhanced tidal dissipation rates are expected when the tidal potential can resonantly excite those inertial modes. Dissipation has then a complex dependence on the tidal frequency, and thus on both the rotational and orbital period. As dissipation depends on the internal structure of the body and on its rotation, the computation of Im[k_l^m(ω)] becomes very heavy when one is interested in the secular effects of the tidal torque in the evolution of the orbital element as the star itself is evolving <cit.>. 

Here, we use a simplified model for tidal dissipation, which does not involve the complicated details of the frequency-dependence of the response functions, and neglect the enhancement of dissipation due to resonances. Indeed, as shown by <cit.>, the typical level of dissipation can be approximated using a simple analytical formulation of the frequency-averaged dissipation ∫^∞_-∞ Im[k_l^m(ω)]  dω/ω, obtained by means of an impulse calculation. Their solution depends on the internal structure of the body, either homogeneous or double-layered, but makes no assumption on the details of dissipation mechanisms, so it is smooth and free of surface effects (like wave reflection at the interface between convective and radiative zones). Their derivations are done for arbitrary degree and order of the tidal components. Taking-up on their work, we use here ∫^∞_-∞ Im[k_l^m(ω)]  dω/ω as a measure of the tidal response in the low-frequency part of the spectrum where inertial waves are found, when the details of the frequency-dependent response curve are filtered out. Thus, we define an average phase-shift Δ̅_l^m for the tidal component of degree l and order m obtained using the frequency-averaged formulation of <cit.>[We recall the reasoning behind the formulation of ∫^∞_-∞ Im[k_l^m(ω)]   dω/ω and derive its expression in terms of global stellar parameters for arbitrary values of l and m in Appendix <ref>.]

    ∫^∞_-∞ Im[k_l^m(ω)]   dω/ω≈ k_l ∫^∞_-∞Δ_l^m(ω) dω/ω≡ k_l Δ̅_l^m.


This averaged phase-shift Δ̅_l^m can be related to an average positive-definite time-lag τ̅_l^m. They are  obtained by assuming that, in the fluid frame, components of different p but same l and m have the same time-lags, which is equivalent to assume that the relevant tidal frequency is ω̂ = m Ω_ e. Since the tidal dissipation rate is a positive-definite quantity, by construction  Δ̅_l^m is also positive-definite. Because ω̂_l, m, p=-ω̂_l, -m ,l-p,  we must have τ̅_l^-m = τ̅_l^m, and thus

    k_l  τ̅_l^m = k_l Δ̅_l^m/|m| Ω_ e.

Subsequently, we compute the tidal power and torque associated to every component of the tidal potential, but instead of using a frequency dependent time-lag associated to Im[k_l^m(ω)] = k_l τ_l, m, pω̂, we approximate it to Im[k_l^m(ω)] ≈ k_l τ̅_l^mω̂. For the axisymmetric m=0 component, the axial component of the tidal torque vanishes, but not the tidal power averaged over azimuths, because tidal deformation modulates the moment of inertia of the star. Hence, we compute the average time-lag for m=0 considering that the relevant frequency is the orbital frequency

    k_l τ̅_l^0  = k_l Δ̅_l^0/Ω_0.

The frequency averaged formulation derived from 's work corresponds to the dissipation of tidally excited inertial waves. Yet those waves exist only when |ω̂| < 2 Ω_ e. Outside of this range, the dissipation efficiency corresponds to the one of the equilibrium tide. Hence, we introduce a frequency dependence for the value of τ̅_l^m by considering

    k_l τ̅_l^m, p=
    
    k_l τ̅_l^m,  if  |ω̂| <2 Ω_ e
    3/2 Ω_0 Q_ eq',  otherwise

where the value of Q_ eq' remains uncertain and we perform calculations considering a range of value to study the sensitivity of the model (see Sec. <ref>).

Following <cit.>, we thus write the average energy transfer ⟨Ė⟩ associated to the component of degree l as

    ⟨Ė⟩ = Ω_0 ∑_m, p (l-2p)(2l+1)/8 πG  M_ p^2 R_⋆^2l+1/a^2l+2
    (W_l^l-2p d^(l)_m, l-2p(Θ))^2 k_l  τ̅_l^m, pω̂_l, m, p,

and the corresponding axial component of the tidal torque T_z

    T_z = ∑_m, p m (2l+1)/8 πG  M_ p^2 R_⋆^2l+1/a^2l+2
    (W_l^l-2p d^(l)_m, l-2p(Θ))^2 k_l  τ̅_l^m, pω̂_l, m, p ,

where we calculate the value of τ̅_l^m, p using Eqs. (<ref>), (<ref>), (<ref>) and (<ref>).

Up to this point, all our derivations are valid for arbitrary value of l and m, and in principle, they allow the computation of the tidal response for any spherical harmonic of the multipole expansion of the tidal potential. As can be seen in Eq. (<ref>), the contribution of the spherical harmonic of degree l is proportional to r^l, so usually only the quadrupolar terms are retained in the formulation of the temporal evolution of orbital elements. Explicitly then, for the quadrupolar components,  Eq. (<ref>) yields

    ⟨Ė⟩ = 3/64G M_ p^2 R_⋆^5/a^6Ω_0  [ 6 sin^4 Θ (2Ω_0)   k_2  τ̅_2^0,2.
    
    +4 sin^2 Θ(1+cosΘ)^2 (2Ω_0-Ω_ e) k_2 τ̅_2^1,0
    
    - 4 sin^2 Θ(1-cosΘ)^2  (-2Ω_0-Ω_ e)  k_2 τ̅_2^1,2
    
    + (1+cosΘ)^4 (2Ω_0-2Ω_ e)  k_2 τ̅_2^2,0
    . - (1 - cosΘ)^4 (-2Ω_0-2Ω_ e)  k_2 τ̅_2^2,2],

and for Eq. (<ref>) we have

    T_z = 3/64G M_ p^2 R_⋆^5/a^6[ 2 sin^2 Θ (1+cosΘ)^2  (2Ω_0 - Ω_ e)  k_2   τ̅_2^1,0. 
    
    + 8 sin^2 Θcos^2 Θ (-Ω_ e)  k_2   τ̅_2^1,1 + 2 sin^2 Θ  (1-cosΘ)^2(-2Ω_0 - Ω_ e)  k_2   τ̅_2^1,2
     
    +  (1+cosΘ)^4  (2Ω_0 - 2Ω_ e) k_2   τ̅_2^2,0 +  4  sin^4 Θ k_2   τ̅_2^2,1  (- 2Ω_ e) 
    .+  (1-cosΘ)^4  k_2   τ̅_2^2,2  (-2Ω_0 - 2Ω_ e) ]


The y-component only contributes to precession that we are neglecting here. The perpendicular component of the torque is computed following the formulation of <cit.>[The computations of the tidal power in <cit.>  (Eq. 25) and <cit.>  (Eq. 41) are equivalent to a factor (2l+1)/8π with the notations we have introduced.]

    T_x =   3/64G M_ p^2 R_⋆^5/a^6sinΘ[ 6 sin^2 Θ (2Ω_0) k_2 τ̅_2^0,2
    
    +  2 (1+cosΘ)^2 (2-cosΘ) (2Ω_0 - Ω_ e) k_2  τ̅_2^1,0 - 8  cos^3 Θ (-Ω_ e) k_2  τ̅_2^1,1
    
    - 2 (1-cosΘ)^2 (2+cosΘ) (- 2Ω_0 - Ω_ e) k_2  τ̅_2^1,2
    
    + (1+cosΘ)^3 (2Ω_0 - 2Ω_ e) k_2  τ̅_2^2,0 -  4 sin^2 ΘcosΘ (- 2Ω_ e) k_2  τ̅_2^2,1
    . - ( 1-cosΘ)^3 (-2Ω_0 - 2 Ω_ e)  k_2 τ̅_2^2,2]




 §.§ Stellar angular momentum modeling

The angular momentum of the star L_⋆ = L_⋆L̂_⋆ is modelled with core-envelope decoupling under the assumptions of the double-zone model introduced by <cit.>;

    L_⋆ = L_ c + L_ e = I_ cΩ_ c + I_ eΩ_ e

where I_ c and I_ e are the moment of inertia for the core and the envelope respectively. The core and the envelope are assumed to rotate as solid bodies with different angular velocities. As noted in <cit.>, in a differentially rotating body there is no simple global relation between the tidal dissipation rate and the rates of energy and angular momentum transfer as computed from Eqs. (<ref>), (<ref>), (<ref>). However, those quantities depend on stellar parameters only through the mean stellar radius at the surface, the static Love number and the average time-lag. The static (real) Love numbers link the interior structure of the star to its external potential through the surface value of the logarithmic derivative of the deformation <cit.>. Here we neglect the effect of the centrifugal force on the shape of the star, so the internal discontinuity in rotation speed is irrelevant. The impulsing forcing approach used to estimate the average time-lag also neglects the effect the Coriolis acceleration on the non-wavelike part of the displacement, so the boundary condition specifying the solution remains the same even if we consider that the core rotates at a different angular velocity than the envelope. Moreover, as is clearly apparent in Eq. (96), (97) and (98) of <cit.>, only the convective envelope of the star is involved in the tidal energy transfer, since we are only considering the dissipation of tidally-excited inertial wave. Here we suppose that the envelope rotates rigidly and as such our assumptions are fully compatible with the derivations of <cit.>, setting the angular velocity to the one of the envelope. Furthermore, we consider that the tidal torque is applied on the envelope of the star, but we allow angular momentum exchanges between the two zones. As a matter of fact, in the piece-wise homogenous rigidly rotating model considered in <cit.>, the non-barotropic nature of the density-jump results in a singularity of vorticity at the interface. In the double-zone model considered here, considering our simplifying assumptions regarding rotation, and other overlooked mechanisms that may be involved, we simply parametrise the transfer of angular momentum between the zones by a quantity Δ L defined as

    Δ L = I_ c I_ e/I_ c+ I_ e(Ω_ c - Ω_ e)

at a rate determined by a coupling time-scale τ_c. For the dependence of τ_c on stellar rotation, we follow <cit.> and consider a two-valued function with

    τ_c =
    
    10  Myr ,   if Ω_ e(t_0) ≥Ω_ crit, 
    
    128  Myr  ,  otherwise.


We follow the evolution of the star from its formation until the end of its main sequence, i.e. before the core starts contracting during the subgiant phase. Initially the star is fully convective and is assumed to rotate rigidly. We consider a phase of disc-locking, where the net effect of interactions with the disc is that of keeping the surface angular velocity of the star constant for the disc lifetime, i.e. 

    Ω_ e/t = 0,  while   t≤τ_ disc.

The disc lifetime τ_ disc may vary, but it is in general shorter than the PMS, as observations show that most primordial discs have disappeared by the first 10 Myr <cit.>. As the radiative core develops, a quantity of material contained in a thin shell at the base of the convective zone, with a velocity of Ω_ e, becomes radiative, producing an angular momentum transfer towards the core

    .L/t|_ growth = (2/3 R^2_ cM_ c/t) Ω_ e.

From the moment the disc disappears, the star experiences angular momentum loss through its magnetised wind. We assume that because of wind braking, the envelope loses angular momentum at a rate given by the following parametric formula <cit.>

    .L_ e/t|_ wind = K_ w(R_⋆/ R_⊙)^1/2(M_⋆/ M_⊙)^-1/2 min(Ω^3_ e, Ω^2_ satΩ_ e),

where K_ w determine the braking intensity and Ω_ sat is the angular frequency threshold defining a saturated (when Ω_ e≥Ω_ sat) regime of angular momentum loss.

In fine, the evolution of the stellar angular momentum is modelled with the following differential equations

    L_ c/t   = -Δ L /τ_c + .L/t|_ growth
    L_ e/t   = + Δ L /τ_c - .L/t|_ growth -.L_ e/t|_ wind + T_z.




 §.§ Tidal evolution of orbital elements

 

The transfer of energy from the orbit to the rotation of the star due to tides results in the variation of the semi-major axis through the relationship

    ⟨Ė⟩ = /tG M_⋆ M_ p/2a = - G M_⋆ M_ p/2a^2a/t.

For hot-Jupiters, the quadrupolar components of the tidal potential is strongly dominant due to the small ratio between stellar radius and semi-major axis, even during the PMS. The evolution of the semi-major axis is straightforwardly obtained using Eq. (<ref>). As can be seen in Eq. (<ref>), tidal components with p=1 do not contribute to the transfer of energy. Moreover, for m=0, the components with p=0 and p=2 are physically identical. So there are only 5 components that contribute to the evolution of the semi-major axis. The importance of each of them depends on their respective tidal frequency, the average time-lag, and a weighting factor that is a function of the instantaneous value of Θ. This quantity is shown in Fig. <ref> (top panel). Note that for Θ = 0, only the l=2, m=2, p=0 component contributes to the evolution of the semi-major axis.

The magnitude of the orbital angular momentum L_ o= || L_ o || is 

    L_ o=M_⋆ M_ p/M_⋆ + M_ pΩ_0 a^2.

Because we consider here that only the envelope is involved in tidal energy transfer, we will consider the evolution of the orientation of the direction of the spin of the envelope so 

    cosΘ = L_ e . L_ o/L_ e L_ o.

By definition, 0 ≤Θ≤π, the orbit is polar when Θ=π/2 and retrograde when Θ> π/2. Supposing that in the absence of magnetic braking, the total angular momentum is conserved we get 

    Θ/t= - T_x/L_ e -T_x/L_ ocosΘ + T_z/L_ osinΘ,

The differential equation governing the evolution of the obliquity can thus be obtained using Eq. (<ref>) and (<ref>). Here, all the 7 physical components of the tidal forcing will contribute to the evolution of the obliquity. Their relative importance is a function of the instantaneous value of Θ and of L_ o/L_ e. It is shown in Fig. <ref> (bottom panel) for L_ o/L_ e = 0.5.





§ EVOLUTION OF TIDAL DISSIPATION AS A FUNCTION OF STELLAR MASS AND EVOLUTIONARY STAGE

        
As can be seen from Eq.  (<ref>), (<ref>), (<ref>), (<ref>), (<ref>) and (<ref>), the coupled evolution of the orbital elements depends on the dissipation efficiency of the tides raised by the three components of the quadrupole moment l=2 of the tidal potential : the m=0 zonal harmonic, the m=1 tesseral harmonic and the m=2 sectoral harmonic. The so-called obliquity-tide (l=2,m=1,p=1), which has the frequency ω̂=-Ω_ e regardless of the orbital frequency, contributes to the tidal torque but has no associated variation of mechanical energy. However, it is important to note that other terms associated with the zonal or tesseral harmonics participate to both angular momentum and energy transfer, but their contribution to the evolution of the semi-major axis is non-null only when sinΘ≠ 0, i.e. when the system is not aligned or anti-aligned. 

The effect of stellar structure and evolution on the frequency-averaged dissipation of the sectoral (l =m = 2) harmonic have already been computed in other studies <cit.>. In the impulsive forcing problem, the value of m affects the energy dissipated by inertial waves through the coupling of the spheroidal and toroidal wavelike velocity components due to the Coriolis force.  The only difference between the energy transfer from the sectoral harmonic and the zonal and tesseral ones stems from the contribution of the toroidal part, related to the coupling coefficient q̃_̃l̃,

    q̃_̃l̃ = 1/l(l^2-m^2/4l^2 -1)^1/2.



  
To compare the evolution of ∫^∞_-∞ Im[k_l^m(ω)]  dω/ω for the sectoral and tesseral harmonics, we computed it for models of a 1 M_⊙ and a 1.2 M_⊙ and two values of initial stellar rotation frequency 3 Ω_⊙ and 25 Ω_⊙. Those values correspond to rotation periods of ∼ 8.5 and 1 days respectively, and are representative of the range of surface rotation periods measured in clusters younger than ∼ 2 Myr for stars in this mass range <cit.>. The stellar structure quantities and their time derivatives are taken from evolutionary models by <cit.>. To obtain the variables of our piecewise-homogeneous fluid model, we compute the mean density in each layer simply using the mass and thickness of the corresponding zone. We plot in Fig. <ref> the temporal evolution of α the core to envelope radius ratio (top panel), and f the envelope to core density ratio ( bottom panel; see also Appendix <ref>) that we get from the stellar models for the two masses considered here.

 Our evolution begins at t_0 ∼ 1.5 Myr, when stars are almost fully convective, justifying our assumption of solid-body rotation as our initial condition. After the disk-locking phase, the angular momentum of the star evolves as a consequence of both structural evolution and magnetic braking (see Sec. <ref>). Here the parameters are taken from <cit.> and summarised in Table <ref>. There is a lack of reliable calibration of the magnetic braking laws in the framework of the double-zone model for stars with M_⋆ > 1.1 M_⊙. Here, we simply assume the same values apply also to our  1.2 M_⊙ model for all parameters, except for K_w, that we divide by a factor 10 as in <cit.>. 


   

To allow easy comparison with the values found in the literature for the equilibrium tide, we also give the conversion of our estimation of the dissipation rate as the traditional equivalent modified tidal quality factor Q', using the relationship

    <Q_l^'m> = 3/2 k_l Δ̅_l^m.

Here we get a value of Q' that is not exactly frequency-averaged, but rather corresponds to a frequency average of the phase lag <cit.>. The results are shown in Fig. <ref>. 

 
 
We see that the m=0 and m=1 components are almost identical due to similar contributions of the toroidal part to the impulse energy, and they are more than one order of magnitude greater than the m=2 component during the main-sequence. This clearly shows that assuming that different components have the same time-lag, as is done in the equilibrium tide theory, severely underestimates the participation of the tesseral harmonics of the tide, when the tidal frequency is within the range allowing the existence of inertial waves.  Overall, we find that the value of  Q^'0_2 and Q^'1_2 assumes value between 1 and 10^6 throughout stellar evolution across slow and fast rotators. For 1.2 M_⊙ mass stars, fast rotators can have Q^'0_2 and Q^'1_2 smaller than 100 throughout their evolution. This is extremely efficient for fluid dissipation, and is comparable to the value of Q' estimated for the rocky bodies of the solar system, estimated in the framework of the equilibrium tide <cit.>. 

We must stress that we are using here a simplified model which considers that the densities in the radiative core and the convective envelope are both homogeneous. For the sectoral harmonic l=m=2, <cit.> has shown, using polytropes of different indexes and in the presence of a solid core, that using various density profiles can change the frequency-averaged dissipation-rate by orders of magnitude depending on the core to envelope radius ratio α <cit.>. The general trend, for α≳ 0.3, is that the dissipation becomes less efficient as the star gets more centrally condensed. In addition, <cit.> have shown numerically that the dissipation rate of the obliquity tide, for a fluid body containing a solid core, is enhanced by 1 to 2 orders of magnitude in the homogenous case compared to the more realistic case of a polytropic density of index 1(we refer the reader to Figs. 4 & 7 in their article). Here, we use a piecewise homogeneous fluid model. Therefore, we cannot exclude that modelling the convective envelope as homogeneous may artificially enhance the dissipation and this could potentially significantly modify our results regarding the timescales of tidal evolution. In addition, we know from previous works <cit.> that the dissipation of tidal inertial waves can vary over several orders of magnitude as a function of the tidal frequency. At a given frequency, it may thus differ from its frequency-averaged value. However, we expect this latter to provide us a reasonable order of magnitude of the dissipation value that we can use in planetary systems secular evolution models.

Dissipation is strongly affected by stellar rotation. In the magnetic braking model considered here, initially fast and slow rotators converge to a common rotation rate after the first Gyr of evolution, but they have very different behaviours during the firsts hundreds of Myr. During this phase, the angular momentum evolution is driven by the evolution of the structure of the star, dominated by its global contraction. Thus, the value of Q' during this phase is very similar for the two stellar masses considered here, with Q^'0_2 and Q^'1_2 ∼ 10^1 - 10^3 and  Q^'2_2 ∼ 10^3 - 10^5 depending on the value of the initial rotation. Those values are maintained during the first 20-30 Myr of evolution. This implies that any hot-Jupiter that would reach a 3-day period orbit in the first tens of Myr would spiral into the star and be destructed by tides within a few 10^5 years, if the tidal frequency falls within the range of inertial waves and if the rotation period of the star is greater than 3 days, regardless of the mass of the host. Yet, one should keep in mind that the values obtained here may be considered as upper bounds on the dissipation because of the hypothesis of homogeneous density in the convective envelope.

 On the main sequence, the stellar structure stops evolving and the evolution of the dissipation is mainly controlled by the rotation frequency of the envelope. For the 1 M_⊙ mass star, we see a steady decrease in dissipation efficiency, with Q^'0_2 and Q^'1_2 ∼ 10^4 - 10^5 and  Q^'2_2 ∼ 10^7 - 10^9. For the 1.2 M_⊙ mass star, dissipation efficiency appears to plateau around Q^'0_2 and Q^'1_2 ∼ 10^1 - 10^3 and  Q^'2_2 ∼ 10^5 - 10^7, because angular momentum loss is reduced and the star spins down less efficiently. After about 4 Gyrs, the mass of the convective envelope becomes negligible and our computation of Q' suffers from numerical instabilities inherent to the stellar evolutionary track we are using. 
These order of magnitudes are in agreement with those obtained by <cit.>, who also computed the impact of the simultaneous structural and rotational evolution of low-mass stars on the dissipation of tidal inertial waves in their convective envelope (we refer the reader to their Figs. 4 and 5 - top left panel) for the l=m=2 components.
 
Interestingly here, unlike what is found in <cit.>, we see that despite the less important convective zone of the 1.2 M_⊙ mass star, the energy dissipated by the inertial waves for this stellar mass is about the same order of magnitude, if not greater, than for the 1 M_⊙ mass star, for a given set of l, m and initial rotation frequency. 
 
This may seem to be also in contradiction with previous numerical simulations which showed that this dissipation rate should scale with the extent of the convection zone <cit.>. But one must remember that those simulations were computed at a fixed rotation rate while here we compute the evolution of the rotation of the star as its internal structure evolves. To better compare our results with these numerical simulations, we show in Fig. <ref> the frequency-averaged dissipation obtained with our model along the evolution of the stellar structure, but fixing the rotation frequency of the star at 8.5Ω_⊙, which corresponds to the rotation period of 3 days, a value chosen arbitrarily by <cit.> and <cit.> to illustrate their results.

Firstly, let us consider the case of the l=m=2 component. For this component, the frequency dependent dissipation rate of inertial waves for a solar mass star rotating with the 3-day period has been computed in  <cit.>. Besides the complex dependance upon tidal frequency, their equivalent Q' ranges between 10^6 and 10^8. Using the same numerical set up, <cit.> obtained, for a 1.2 M_⊙ mass star rotating at the same angular velocity, a significantly lower Q', ranging between 10^8 and 10^10 for the same tidal component. Here, on the main sequence, we obtain Q'∼ 10^6 and Q'∼ 10^7 for the respective masses. Considering again that our model assumes piecewise homogenous density, which leads to upper bounds on the dissipation, our results are thus in reasonable agreement with the trends of these more realistic numerical simulations, albeit the difference in dissipation as a function of mass is less contrasted. Again, this is in broad agreement with <cit.> (compare their Fig. 4 top-left and bottom-left panels) where the difference in dissipation between a 1 M_⊙ and a 1.2 M_⊙ mass star is found to be greater when considering only the effect of structure, compared to the case when the evolution of both the structure and rotation is taken into account. Comparing Fig. <ref> and <ref> (and also top-left and bottom-left panels of Fig. 4 in <cit.>, although in their case the 1.2 M_⊙ mass star always dissipates less), it is clear that the main factor affecting the difference in dissipation on the main sequence between the two stellar models considered here is thus the evolution of stellar rotation : the 1.2 M_⊙ mass star spins down less rapidly and thus maintains a dissipation rate throughout its main sequence similar to those of the 1 M_⊙ mass star independently from the thickness of their respective convective envelope.

 


Nevertheless, for a given initial rotation rate, as previously noted, the dissipation of the m=0 and m=1 components is always orders of magnitude greater than for the m=2 component. This seems to disagree with previous results from numerical simulations computed in <cit.>, who showed that the dissipation evaluated for different values of m are similar in magnitude. Such a discrepancy is also observed between the frequency-averaged estimation and numerical results obtained by <cit.> for the dissipation rate of the obliquity tide in the case of an homogeneous fluid containing a solid core. Indeed, they show in their Fig. 4 the dimensionless dissipation rate D̃ of the obliquity tide as defined in their Eq. (69) for a homogeneous fluid. It ranges between 10^-5 and 1 depending on the value of the fractional core radius. However, Fig. 6 of <cit.> (second red line from the top, also partly reproduced in Fig. <ref>) shows that the frequency-averaged, dimensionless, rotation normalised, equivalent dissipation for the same homogenous fluid is greater than 10 and diverges as α→ 1. The dissipation rate computed in <cit.> is not exactly equivalent to the quantity obtained in the impulsive response formulation in <cit.>. However, the forcing frequency is independent of the orbital frequency of the planet for the obliquity tide, and both quantities, if we relate them through an equivalent quality factor Q', should agree to a factor of the order of unity. This does not seem to be the case and may be related to the effect of precessional forces.

In addition, compared to the fluid case, the boundary conditions imposed by a solid core have large repercussions on the functional dependence of the frequency-averaged tidal dissipation with the core size. This can be appreciated on Fig. <ref>, where we plot the result of the frequency-averaged dissipation, divided by ϵ^2=Ω^2 (G M_⋆/R_⋆), in the cases of a homogenous fluid with a solid core and of the piecewise fluid model considered here, for all components of the quadrupolar tidal forcing. 

First, we can clearly see that the piecewise homogenous fluid has overall less efficient dissipation for all values of m. This is expected since the core in this model is not perfectly rigid but instead undergoes some tidal deformation. This produces slower horizontal flow than in the solid core case, and through the Coriolis force, a smaller wavelike velocity. 

Next, in the range of parameters corresponding to main sequence values (see Fig. <ref>), the dissipation of the zonal and tesseral components clearly dominates that of the sectoral component in the case of the piecewise fluid model. This quantitative difference exists also in the solid-core case but is mainly important for thick convection zones. 

Here in the piecewise homogeneous fluid model, the integrated response to a tesseral harmonic is always much stronger than the one of the sectoral harmonic. This seems to be a robust result for the dissipation computed for an impulsive forcing. However, the question of the effect of the density profile on the relative dissipations of the different tidal components has not yet been investigated beyond the two simple models discussed here. This requires further theoretical investigations that are beyond the scope of this article. Let us note here that, in the present model, the difference in stellar structure between a 1 M_⊙ and a 1.2 M_⊙ mass star on the main sequence (mainly stemming from the density of the envelope) has little effect on the tidal dissipation for the m=0 and m=1 components, but has a greater impact on the m=2 component.

Again, the predominant factor here will be the evolution of stellar rotation. Since the 1.2 M_⊙ mass star remains a fast rotator throughout its main-sequence, it has overall a more efficient dissipation, if we consider that our impulsive response formulation provides a robust measure of the dissipative properties of the star.


This would tend to produce an anti-correlation of obliquity with stellar effective temperature, contrary to what is observed. Nevertheless, this result is very sensitive to the model chosen for angular momentum evolution. For example, <cit.> use a different prescription for the angular momentum evolution of the star. They use stellar evolution models that take into account rotation, with angular momentum transport operated through meridional circulation and diffusion by shear turbulence. This results in a much less efficient coupling between the radiative zone and convective envelope than the model used here <cit.>. Moreover, in their model, the magnetic braking law involves a parametrisation based on the Rossby number. Consequently, for the l=m=2 mode, they find that dissipation scales with the extent of the convective zone for stellar mass between 1 and 1.4 M_⊙ on the main sequence, with <Q_2^'2> increasing by about an order of magnitude over this mass range. Here our results differ because we allow the coupling of the radiative and convective zones, which affects the rotation rate of the envelope and thus the value of the dissipation. We estimate tidal dissipation efficiency using a simplified model, valid in the convective zone only and which assumes rigid body rotation, but our computations show that even in this framework, the internal angular momentum transport throughout the star may have a strong impact on tidal dissipation efficiency.

To conclude, we find different dissipation efficiencies associated to different modes, and this supports the possibility of having different timescales for the tidal evolution of the semi-major axis and the obliquity. But the absolute value of the dissipation is also very sensitive to the evolution of the rotation of the star. With the magnetic braking law used here, this would tend to produce an anti-correlation of obliquity with stellar effective temperature, contrary to what is observed. In any cases, the value we find here would imply typical timescales for tidal evolution that could be much shorter than the main-sequence life-time of the star, depending on the orbital configuration. Thus the final state of Hot-Jupiter systems is not easily inferred from the relative values of Q'^l_m alone, but must also consider the initial rotation period of the star and the coupled evolution of orbital elements.





§ COUPLED EVOLUTION OF SEMI-MAJOR AXIS AND OBLIQUITY

Using Eqs. (<ref>), (<ref>), (<ref>),  (<ref>), (<ref>) and (<ref>), we compute the coupled temporal evolution of the orbital elements for different values of initial obliquity, semi-major axis and stellar rotation. Our evolution begins again at t_0 ∼ 1.5 Myr, when stars are almost fully convective, and thus rotating as a solid-body. It stops at the end of the main-sequence phase (i.e. after 10 Gyr for the 1 M_⊙ star and 5 Gyr for the 1.2 M_⊙ star) or when the planet reaches the Roche limit (a_ R=2.422 R_ p√(M_⋆/M_ p), with R_ p = 1.2 R_ Jup for a Jupiter-mass planet), or alternatively when the star is spinning faster than the break-up velocity (Ω_ B = √(G M_⋆/ R_⋆^3)). During the disk-locking phase, we consider that the planet is maintained on its orbit by its interaction with the disk. Since the rotational evolution of the star is clearly important for tidal evolution, and the distinction between initially slow and fast rotator may be paramount, we opted to relate the disk-lifetime τ_ disk and the coupling time-scale τ_ c to the initial rotation rate of the star following the prescription of <cit.>. However, keeping the prescription of <cit.> for the other parameters as in the previous section resulted in problematic evolutions in the case of initially fast rotating 1.2 M_⊙ mass stars. For most of the initial parameter space explored, tidal evolution would lead either to the star spinning up to the break-up velocity or the planet reaching the Roche limit before the first few Myrs. This suggests that the chosen angular momentum loss model may not apply to this kind of host. Therefore, and in the lack of calibration for angular momentum loss rate in the framework of the double-zone model for 1.2 M_⊙ mass stars, we opted for a parametrisation akin to the one proposed by <cit.> for the solid rotation case. Namely, we take the same value for the scaling factor K_ w for all masses, and we include the mass dependence  through the value of the saturation threshold  Ω_ sat. In this prescription,  Ω_ sat  is inversely proportional to the convective turnover timescale in the star τ_ cz

    Ω_ sat =  Ω_ sat_⊙τ_ cz_⊙/τ_ cz,

where  Ω_ sat_⊙ and τ_ cz_⊙ are normalisation values and where τ_ cz is computed from the effective temperature on the ZAMS according to the prescription of <cit.>.  Table <ref> and <ref> lists the values of all adopted parameters for stellar angular momentum loss in the present section. This slightly different magnetic braking model produces a qualitatively similar evolution of dissipation as the one presented in Fig. <ref>, but now, most of the 1.2 M_⊙ hosts do not enter the saturated regime of the magnetic braking law. The peak of their rotation speed is not as high during the initial spin up phase, so even though they do not enter the saturated regime, their instantaneous rate of angular speed loss is weaker than for the 1 M_⊙ stars, during the first Gyr of evolution.





   
          
Finally, we performed orbital evolution calculations considering a range of the value of Q_ eq' defined in Eq. (<ref>) (corresponding to the frequency range where inertial waves cannot exist). It turns out that in most case, the results do not change significantly for Q_ eq'≃ 10^6 - 10^10 <cit.>. Thus we present here the results obtained with Q_ eq'=10^7. We integrate the set of equations presented in Section <ref> using a fourth-order Runge-Kutta integrator with an adaptive time-step routine <cit.>. Examples of evolutions are given in Fig. <ref> and for a 1 M_⊙ mass star and Fig. <ref>.

Let us first consider the case of a solar-mass star orbited by a Jupiter-mass planet at 0.1 AU at the start of the evolution (Fig. <ref>). On the top panel, we see that the evolution of the obliquity depends mainly on the initial rotation frequency of the star, whatever the obliquity's initial value. This is expected in our model where the efficiency of tidal dissipation is directly related to the rotation frequency of the star, knowing that in this orbital configuration, inertial waves can be excited throughout the first Gyr of evolution, as can be seen on Fig. <ref>. 

For initially slow rotators (solid lines), tidal dissipation is not efficient and this result in very little evolution of both the semi-major axis and the obliquity over the lifetime of the system. For initially moderately rotating stars (dotted-lines), tidal dissipation becomes significant and the value of sinΘ has an important contribution to typical timescale of evolution of the semi-major axis. For both quasi-aligned prograde and retrograde system, sinΘ≈ 1 and both the obliquity and the semi-major axis do not change significantly in time. But for 80^∘≲Θ≲ 120^∘,  sinΘ << 1, the semi-major axis decreases very rapidly and the planet reaches the Roche Lobe during the first 10 Myrs of evolution. However in this case, tidal evolution does not change significantly the value of the obliquity before the final engulfment. 

For initially fast rotators (dashed lines) the prograde and retrograde orbits do not behave in a mirrored way. For systems starting almost aligned (Θ < 60^∘), the typical timescale for evolution of the semi-major axis and the obliquity appear to be similar, and faster than for the initially moderately rotating star. Since at the beginning of the evolution, Ω_0 << cosΘΩ_ e, the tides tend to push the planet on a wider orbit, and the obliquity is damped as the start contracts and Ω_ e increases. When  Ω_ e reaches its maximum, the star enters its phase of spin down in the saturated regime, and tidal dissipation becomes less efficient. Essentially, both the semi-major axis and the obliquity stop evolving and they keep the value they reached after the first ∼ 40 Myrs of evolution. The almost polar orbit starting with Θ=80^∘ has a somewhat different behaviour starting with a faster evolution of the semi-major axis which, in this case, decreases in time, because cosΘ≈ 0 and thus Ω_0 >> cosΘΩ_ e. When the star starts spinning down, the tidal regime changes and the evolution of the obliquity becomes faster than that of the semi-major axis. This changes again when the frequency ratio between the star and the orbit does not allow the excitation of the dynamical tide any more and tidal evolution resumes, damping the obliquity but with little effect on the semi-major axis. 

On the other hand, systems with fast rotators starting retrograde with Θ > 100^∘ are much more unstable and the planet is immediately engulfed as negative angular momentum is transferred (from the orbit to the star because |Ω_0| << |cosΘΩ_ e|) into the axial component of the orbital angular momentum through an increasing value of Θ because the system is retrograde.

The evolution starting with Θ=160^∘ is a remarkable example. As the obliquity reaches perfect retrograde alignement (Θ=180^∘) the m=0 and m=1 components of the tides cease to dissipate energy, and the rate of evolution of the semi-major axis reduces dramatically. The system is now capable of maintaining a retrograde very close-in orbit for about two Gyrs, the star continuously being spun down to the point of reaching sub-solar rotation by the end of the evolution. Within the space of initial conditions considered here, this is the only track for which tidal evolution affect significantly the spin of the star on secular timescales.






Finally Fig. <ref> shows similar evolutionary track but this time for a 1.2 M_⊙ star. Qualitatively, we find again the same respective differences between the outcome of tidal evolution of initially slow, moderate and fast rotators. However, the star is now rotating faster throughout its main sequence, and consequently maintains a strong dissipation for most of its lifetime. This results in a faster evolution of the obliquity even around slowly rotating stars. In this set of simulations, all systems reach values of obliquity inferior to 10^∘ after the first Gyr of evolution, except some initially retrograde systems that settle on well-aligned retrograde orbits.










The same general trends are observed if the planet is starting its evolution on a shorter orbit. Fig.<ref> show the same kind of simulations but with now a Jupiter-mass planet starting on a 0.05 AU orbit, close to where hot-Jupiters are observed. Now in these orbital configurations, inertial waves can be only be excited episodically during the evolution, especially for initially slow rotators, as can be seen on Fig. <ref>. However, there is little qualitative difference compared to the previous case. Typical timescales of evolution are faster, but the obliquity tend to evolve faster than the semi-major axis, in all cases. However, the 1.2 M_⊙ star remains more efficient at damping any initial non-null obliquity. This suggests that tidal dissipation in the convective zone alone is not enough to interpret the observed correlation between effective temperature and obliquity. Clearly others physical mechanisms are at play.



§ CONCLUSIONS

We have implemented for the first time a model of tidal evolution for the obliquity in circular exoplanetary systems that uses an ab-initio estimation of tidal dissipation based on the physical properties of the host star, and tied to a known mechanism, the excitation and damping of inertial waves in convection zones. Even though we don't apply a full frequency-dependent treatment of the dynamical tide, our formulation retains some dependence on tidal frequency, through the use of different tidal dissipation efficiency for different components of the tidal potential. Moreover, we also condition the value of dissipation on the range of tidal frequencies where inertial waves are excited <cit.>. We use here a simplified piecewise homogeneous fluid body model. This may enhance the obtained value of the dissipation compared to a case where the realistic density profile is taken into account <cit.>. Moreover, the actual dissipation at a given frequency could differ from the frequency-averaged dissipation by several orders of magnitude <cit.>, so that we cannot exclude that real systems behave differently from this simplified model <cit.>.

We find that the frequency-averaged tidal dissipation is very sensitive on the angular velocity of the convective envelope, as it scales as Ω_ e^2 for all modes. We also find that for the quadrupole moment of the tidal potential, the dissipation due to the sectoral harmonic (l = 2, m = 2) is about one to two orders of magnitude less efficient than the dissipation due to both the zonal (l = 2, m=0) and the tesseral (l = 2, m = 1) harmonic, for both models of 1 M_⊙ and 1.2 M_⊙ mass stars. We predict that the dissipation in the 1.2 M_⊙ mass star is not significantly reduced compared to the one of the 1 M_⊙ mass star for all modes considered here, despite a smaller convective zone. This is due to the fact that the evolution of tidal dissipation in convective envelopes of low-mass stars is mainly driven by their rotational evolution on the main-sequence <cit.>. More massive stars remain fast rotators until the end of the main sequence, as a consequence of a less efficient magnetic braking. This explains why their dissipation is comparable to the one of less massive stars.

In this framework, typical timescales of evolution for the orbital elements cannot be simply formulated in terms of global parameters. However our simulations show that, for typical orbital configurations of hot-Jupiters, the damping timescale of the obliquity is generally shorter than that of the semi-major axis, except for the case of rapidly rotating stars with planets on very inclined orbits. This is a formal proof of concept of the conjecture put forward by <cit.>, where the dynamical tide does indeed lead to different speed of evolution for the semi-major axis and the obliquity.

Nevertheless, our simulations fail to produce slower tidal evolution for more massive stars, in the range of mass corresponding to low-mass stars developing a radiative core and a convective envelope. The obliquity and/or the semi-major axis are damped before the star reaches the main-sequence for both stellar masses considered here, when the initial rotation period is shorter than ∼ 3 days (i.e. Ω_ e≳ 8.5 Ω_⊙) and the initial semi-major axis is smaller than ∼ 0.05 AU (i.e. P∼ 4 days). For planets starting with Θ≳ 150^∘, tidal dissipation can bring the planet on a meta-stable retrograde aligned orbit (Θ =180^∘), as found in other studies <cit.>. For initially slow rotators, an initial misalignment can be maintained throughout the life of the system, but in this case, the obliquity of planets orbiting more massive host is consistently damped on shorter time-scales.

This would thus fail to explain the observed correlation between effective temperature and obliquity. However, many factors are at play here, and our simulations show that the final outcome of tidal evolution is very sensitive to the initial conditions. In addition, the few first Myr of evolution are determinant for the final outcome. This is also where our model is more affected by the parameters of the magnetic braking law, such as the duration of the disc-locking phase, the coupling timescale for exchange of angular momentum between the radiative core and the convective envelope, or the saturation threshold. 


One of the most severe limitations of magnetic braking laws using the double-zone model is the paucity of calibration points. To this day, there is only a handful of main sequence Sun-like stars for which constraints on internal differential rotation and age can be obtained using asteroseismology <cit.>. Precise measurements of internal differential rotation are only available for the Sun <cit.>, where it is measured to be weak through most of the interior, except maybe in the innermost region of the core. There is also a general lack of stars of known age, rotation period and mass for masses smaller than 1 M_⊙ and ages greater than 600 Myr <cit.>.

Clearly those parameters play a role as important as the actual mechanism responsible for tidal dissipation in the global outcome of tidal evolution. A better understanding of the fundamental underlying mechanisms is required to be able to retrace the tidal evolution of hot-Jupiter, and eventually understand their migration and formation.

But even with a perfect understanding of magnetic braking, and internal angular momentum transport <cit.>, this study would just be a first step in the theoretical effort to produce tidal models accurate enough to constrain migration scenarios. Indeed, in all scenarios of migration that do not involve interactions with the protoplanetary disk, the eccentricity of the proto-hot Jupiter must be excited to high values to allow its subsequent circularisation on a close-in orbit <cit.>. Our model here is not suitable for the treatment of eccentricity. In the classical equilibrium tide, it is easy to show that the circularisation of the orbit in a star-planet system is mainly produced by the dissipation of tides raised in the planet. Here we neglect those tides, so that our model would need to include them to also account for the evolution of the eccentricity. Since we are neglecting the interactions of tidal bulges between them, the orbital evolution due to the tides raised and dissipated in the planet can simply be obtained by adding the torque resulting from the response of the planet, to the point-mass-like tidal potential exerted by the star. The efficiency of tidal dissipation on the planet could also use the frequency-averaged formulation of <cit.>, assuming that the planet could be considered as a homogeneous body with a solid core <cit.>. Because this formulation is very sensitive to the core size, this raises the question of the mass of the solid core in gas-giants, for which there is still no definite observational evidence, even in the case of our own Jupiter. Finally, in the case of an eccentric orbit, the dynamical tide assumes an infinite number of terms, and non-dominant harmonics can be important through resonant contributions. In this case, the frequency-averaged dissipation would be needed for arbitrary value of l and m, and the analytical formulation we have made explicit in Appendix can be used. The introduction of the eccentricity in the model is thus the next step to improve models of tidal evolution and constrain migration scenarios. 

To get a complete picture, it would also be necessary to take into account the dissipation of tidal gravity waves <cit.> <cit.> in the radiative cores of low-mass stars through linear thermal diffusion <cit.> or non-linear breaking mechanisms <cit.>. This dissipation can be of the same order of magnitude that the one of tidal inertial waves in the external convective envelope <cit.>. Therefore, tidal gravity (gravito-inertial) waves can also play a key role to shape the orbital dynamics of short-period systems and the rotational evolution of their components <cit.>. As in the case of tidal inertial waves, the dissipation of tidal gravito-inertial waves in the radiative core of low-mass stars depends on their mass, their age, and their rotation. This dependence, combined with observational constraints on the orbital state of short-period systems, can also allow us to constrain the mass and the evolutionary state of the host star <cit.>. Besides, the effects of differential rotation, both in convective and in radiative zones, should also be taken into account, as it is expected that the range of frequencies at which tidal (gravito-)inertial waves may propagate is broader than in the case of solid-body rotation <cit.> while corotation resonances may enhance their dissipation <cit.>. Finally, one should also consider the effects of the elliptic instability <cit.> and of magnetic fields on tidal inertial and gravito-inertial waves <cit.>.

Efforts to confront observations of hot-Jupiters and migration theories require also better observational constraints. In a very near future, these studies will benefit from the observations of the ground-based near-IR SPIRou instrument <cit.> and of the space missions TESS <cit.> and CHEOPS <cit.>. Currently, one of the most stringent limitation to test migration theories is the poor constrain we have on the age of most known hot-Jupiter systems. Besides, measuring the internal rotation profile of stars, possibly down to the core <cit.>, would help to better understand the mechanisms involved in the evolution of their angular momentum <cit.>. In this context, it is clear that a mission like PLATO <cit.>, which will provide the complete characterisation of host stars using asteroseismology, will enable a breakthrough for dynamical studies, thanks to the precise and accurate determination of stellar radii and masses, differential rotation and, above all, ages.


The authors are grateful to the referee, Adrian Barker, for several valuable comments on the first version of the manuscript that helped improved the quality of the work. CD wishes to thank J.P. Marques, M. Janvier, T. Appourchaux, R. Cameron and E. Bolmont for helpful discussions. CD acknowledges funding from the ANR (Agence Nationale de la Recherche, France) program IDEE (ANR-12-BS05-0008) "Interaction Des Étoiles et des Exoplanètes", and from the German Space Agency (Deutsches Zentrum für Luftind Raumfahrt) under PLATO grant 50OO1501. SM acknowledges funding from European Research Council (ERC) under SPIRE grant (647383), and from CNES PLATO grant at CEA Saclay.


 aa 
  






§ ANALYTICAL INTEGRATED RESPONSE FOR A PIECEWISE-HOMOGENEOUS FLUID BODY

In <cit.>, the computation of the frequency-averaged dissipation ∫^∞_-∞ Im[k_l^m(ω)]  dω/ω for inertial waves in a piece-wise homogeneous body is given in Appendix B for the l=m=2 mode. It is the only mode required to compute the tidal evolution of coplanar circular systems. Here we need to compute this value for other values of m, namely m=0 and m=1. In this appendix, we recall the main steps leading to the formulation of ∫^∞_-∞ Im[k_l^m(ω)]  dω/ω as derived by <cit.>. We derive the formulation of ∫^∞_-∞ Im[k_l^m(ω)]  dω/ω for arbitrary values of l and m, which is not provided explicitly in  <cit.> but is needed to study the obliquity tides.



 §.§ Low-frequency tidal forcing of a rotating barotropic fluid

The star is considered as a steady axisymmetric body weakly perturbed by a tidal gravitational potential Ψ (see Eq. <ref>). We focus on its convective zone that rotates with uniform angular velocity Ω. It is also considered to be a barotropic ideal fluid of pressure p uniquely related to the density ρ. Finally, we assume isentropy so we can write d p /ρ = h, h being the specific enthalpy. In the reference frame rotating with the envelope, the linearised equations of motion, mass conservation and gravitational potential, are 

    ξ̈ + 2Ω×ξ̇   =    - ∇ W, 
    
    h'+Φ'+Ψ   =   W, 
    ρ'    =    - ∇· (ρξ),
    ∇^2 Φ'    =    4 π G ρ' ,

where the prime denotes an Eulerian perturbation, the dot denotes ∂/∂ t, ξ is the displacement. The external tidal potential is Ψ and Φ' represent the internal self-gravitational potential perturbation. Thus W is a reduced pressure perturbation. In the low-frequency limit, this quantity can be neglected to leading order and thus the right-hand side of Eq. (<ref>) can be set to 0. Even with this simplification, solving the system (<ref>)-(<ref>) is a task that requires further approximations if progress is to be made analytically.

Thus solutions are sought by decomposition : all quantities are supposed to be the sum of a non-wavelike part and a wavelike part. The non-wavelike part is an instantaneous hydrostatic response to the tidal potential. It is calculated by neglecting the Coriolis force, but involves the correct Eulerian perturbations of internal potential, pression and density induced by the tidal forcing. The wavelike part is the component of the response that is driven by the Coriolis force, and satisfies the anelastic approximation where acoustic waves are filtered out, and rigid boundary conditions on surfaces with both vacuum and the core. The non-wavelike part therefore satisfies the following equations:

    ξ̈_ nw   =    - ∇ W_ nw, 
    
    h'_ nw+Φ'_ nw+Ψ   =   0, 
    ρ'_ nw   =    - ∇· (ρξ_ nw),
    ∇^2 Φ' _ nw   =    4 π G ρ'_ nw,

together with the wavelike part that verifies

    ξ̈_ w + 2Ω×ξ̇_̇ ̇ẇ   =    - ∇ W_ w + f,
    ∇· (ρξ_ w)    =   0, 
    
    h'_ w   =   0,
    Φ'_ w   =   0, 
    ρ'_ w   =   0,

where f = - 2 Ω×ξ̇_ nw is the Coriolis force per unit mass acting on the non-wavelike part which drives the wavelike-like part of the solution. The boundary conditions are also decomposed as

    ξ_ nw,r   =   0   at   r=α R_⋆, 
    ξ_ nw,r   =   -Φ' + Ψ/g  at   r=R_⋆,
    ξ_ w,r   =   0   at   r=α R_⋆,
    ξ_ w,r   =   0   at   r=R_⋆,

where R_⋆ is the mean radius of the free surface of the star and α =R_ c/R_⋆ with R_ c the mean radius of the radiative zone.

The system of ordinary differential equations for the wavelike and non-wavelike parts, together with their respective boundary conditions at interfaces, form our initial-value problem. This can be solved by means of Fourier (or Laplace) transforms methods. Let us consider the case where the tidal forces derives from a potential Ψ =  Re[A_l(r/R_⋆)^l Y_l^m(θ, ϕ) e^-iω̂ t], in the fluid frame with coordinates (r, θ, ϕ) defined in Sec.<ref>, and A_l an amplitude of appropriate dimensions. Taking the Fourier transform of Eqs. (<ref>) and (<ref>) is formally equivalent to assume that the forcing and the response depend harmonically in time through a common factor e^-iω̂ t.



 §.§ Impulsive forcing

It is important to stress that even if the equations (<ref>)-(<ref>) contain time derivatives, they do in fact imply that the response is instantaneously related to the tidal potential, and the corresponding boundary conditions are instantaneous in time. Taking the Fourier transforms of those equations for an impulsive forcing thus brings back the problem to an hydrostatic setting with only spatial dependence . If we consider a tidal potential with the temporal dependence 

    Ψ = Ψ̂(r, θ, ϕ) H(t),

where H(t) is the Heaviside step function, we place ourself in the special case where the wave-like part of the response experiences an impulsive effective force that derives from this potential. In this case, the spatial structure of the tidal potential

    Ψ̂ =  Re[A_l(r/R_⋆)^l Y_l^m(θ, ϕ) ]=  Re[Ψ̂_̂l̂(r) Y_l^m(θ, ϕ) ],

and the associated perturbed internal gravitational potential 

    Φ̂'_ nw =  Re[Φ̂'_ nw,l(r) Y_l^m(θ, ϕ) ],

are the only components intervening in the Fourier-transformed problem. The Fourier transforms of Eq. (<ref>) and (<ref>) combine to produce the inhomogeneous Helmholtz equation

    ∇^2 Φ̂_ nw' + 4 π G ρ/v_s^2 (Φ̂_ nw' + Ψ̂) = 0,

where v_s is the sound speed, which projects on the considered spherical harmonics

    1/r^2 d/ dr(r^2  dΦ̂'_ nw,l/ dr) - l(l+1)/r^2Φ̂'̂_ nw,l+ 4 π G ρ/v_s^2 (Φ̂'̂_ nw,l+ Ψ̂_l) = 0

 in R_ c < r< R_⋆. The known tidal potential Ψ̂_l is of the form Ψ̂_l = A_l(r/R_⋆)^l Y_l^m(θ, ϕ) for the impulse problem and Φ̂'̂_ nw,l is the internal perturbed gravitational potential to be solved for. For the two-layered model, the solutions to this equation can be written under the form 

    Φ̂'̂_ nw,l =
    
     (B_1 + B_2)(r/R_⋆)^l,     0 < r < R_ c, 
    
     B_1 (r/R_⋆)^l + B_2 α^2l+1 (r/R_⋆)^-(l+1),     R_ c < r< R_⋆,
    
    (B_1 + B_2 α^2l+1) (r/R_⋆)^-(l+1),    r>R_⋆.

The coefficients must satisfy the matching conditions to ensure that the boundaries are equipotentials

    (2l+1) B_2    =    3(1-f) (A_l +B_1+B_2)
    
    (2l+1) B_1    =   3f/f+(1-f)α^3 (A _l+B_1+B_2α^2l+1),

with f=ρ_ e/ρ_ c and α =R_ c/R_⋆, where ρ_ c and ρ_ e denote the density of the core and of the envelope respectively.
We give here the explicit algebraic solutions for B_1 and B_2 as a function of the stellar parameters, which are not made explicit in <cit.>

    B_1 = A_l  1+q_l/p_l-q_l


    B_2 = A_l  (1+p_l/p_l-q_l)( 1- f  /2l-2/3 +f )

where we have introduced

    p_l    =2l-2/3(1+(2l+1/2l-2)(1/f-1)α^3)
    
    q_l    = 1- f /2l-2/3 + f α^2l+1.

Next, we use Eq. (<ref>) to find the non-wavelike tidal displacement by using the change of variable W_ nw = Ẍ so that its Fourier transform can be re-written as 

    ξ̂_ nw = - ∇X̂,

where X̂ =  Re [X̂_l (r) Y_l^m(θ, ϕ)]. It projects on the spherical harmonics Y_l^m as 

    1/r^2 d/ dr(r^2 ρ dX̂_l/ dr) - l(l+1)/r^2ρX̂_l = - ρ/v_s^2 (Φ̂'̂_ nw,l+ Ψ̂_l)

for R_ c < r< R_⋆. The solution satisfies Poisson's equation and can be written in the form

    X̂_l =
    
    C_1(r/R_⋆)^l,     0 < r < R_ c, 
    
    C_2 (r/R_⋆)^l + C_3 α^2l+1 (r/R_⋆)^-(l+1),     R_ c < r< R_⋆.

The regularity of X̂_l at each interface requires dX̂_l/ dr = (Φ̂'̂_ nw,l+ Ψ̂_l)/g. This implies the following relationships between the coefficients C_1, C_2, C_3, B_1 and B_2.

    l C_1 = l C_2 - (l+1) C_3    =   3 (A_l + B_1 + B_2)/4 π G ρ_ c,
    
    l C_2 - (l+1) C_3 α^2l+1   =   3 (A _l+ B_1 + B_2 α^2l+1)/4 π G ρ_ c (f + (1-f)α^3).

Again this can be written in terms of densities and radii for the two-layered model as

    C_1 = A_l/4 π G l ρ_ c(2l+1/2l-2/3+f)( 1+p_l/p_l-q_l),


    C_2 = A_l/4 π G l ρ_ c(1/f+(1-f)α^3) 
    (3 - (f/2l-2/3+f)(α^2l+1 + (2l+1)(1-α^3)/1-α^2l+1)) ( 1+p_l/p_l-q_l),


    C_3 = A_l/4 π G (l+1) ρ_ c(3α^2l+1 - (2l+1)α^3 + 2l-2/1-α^2l+1) 
    (1-f/(2l-2/3+f) (f+(1-f)α^3))( 1+p_l/p_l-q_l).

This gives the explicit form for X̂ that determines the non-wavelike displacement ξ̂_ nw. From ξ̂_ nw we deduce the impulsive effective force f̂ = - 2 Ω×ξ̂_ nw that allows us to solve for the wavelike part of the response.


Assuming that the fluid is at rest before the impulse, and considering that the Coriolis force does not act on the wavelike velocity during the impulse process (besides providing the effective restoring force from its action on the non-wavelike velocity), the wavelike velocity û_ w is given by

    û_ w = f̂ - ∇Ŵ_ w,

where the pressure perturbation Ŵ_ w satisfies the anelastic approximation

    ∇· (ρû_ w) = 0

and the boundary conditions û_ w,r=0. The wavelike velocity can be projected on to vector spherical harmonics, using standard methods used to describe stellar oscillations. We note â_l(r) and b̂_l(r) the coefficients representing the spheroidal part and ĉ_l± 1(r) the toroidal part. The wavelike velocity is written as

    û_ w =  Re [e_r â_l^m(r) Y_l^m + r^2 b̂_l^m(r) ∇ Y_l^m 
    
    -r^2 ĉ_l-1^m(r) e_r  ×∇ Y_l-1^m - r^2 ĉ_l+1^m(r)e_r  ×∇ Y_l+1^m ].


These spheroidal and toroidal components satisfy ordinary differential equations according to Eq. (<ref>), which involves X̂_l found previously from the non-wavelike part, and Ŵ_ w of the wavelike part. If we write Ŵ_ w =  Re [Ŵ_l^m(r) Y_l^m(θ, ϕ)], it can be shown that Ŵ_l^m satisfies

    1/r^2 d/ dr(r^2 ρ dŴ_l^m/ dr) - l(l+1)/r^2ρŴ_l^m = - 2 i m Ω/r dρ/ drX̂_l

with boundary conditions dŴ_l^m/ dr = -2imΩX̂_l /r at r=R_ c and r=R_⋆. Consequently, Ŵ_l^m satisfies Poisson's equation and the solution is 

    Ŵ_l^m =
    
    -2 i m ΩC_1/l(r/R_⋆)^l,     0 < r < R_ c, 
    
    -2 i m Ω[C_2/l (r/R_⋆)^l - C_3/l+1α^2l+1 (r/R_⋆)^-(l+1)],     R_ c < r< R_⋆.


Using Eqs. (<ref>) and (<ref>), the spheroidal and toroidal parts of the initial velocity just after the impulse can be made explicit. This yields null values for the spheroidal component for all values of l, â_l^m=b̂_l^m=0, and a toroidal part given by


    ĉ_l-1^m =
    
    -2 Ωq̃_l^m r^-2 (2l +1) C_1 (r/R_⋆)^l,     0 < r < R_ c, 
    
    -2 Ωq̃_l^m r^-2 (2l +1) C_2 (r/R_⋆)^l,     R_ c < r < R_⋆,

and 

    ĉ_l+1^m =
    
    0,     0 < r < R_ c, 
    
    -2 Ωq̃_l+1^m r^-2 (2l +1) C_3 α^2l+1 (r/R_⋆)^-(l+1),     R_ c < r < R_⋆,

where q̃_l^m is a coefficient arising from the coupling of spheroidal and toroidal velocity components by the Coriolis force, given in Eq. (<ref>). The energy transfer associated to the impulsive response is given by

    Ê=1/2∫ρ |û_ w|^2  dV,

which, thanks to the projection given in Eq. (<ref>), is simply found to be

    Ê=Ê_l^m+Ê_l-1^m+ Ê_l+1^m

with 

    Ê_l^m   =   0,
    Ê_l-1^m   =   1/4∫_R_ c^R_⋆ρ r^2 [l(l-1)r^2 |ĉ_l-1^m|^2] dr, 
    Ê_l+1^m   =   1/4∫_R_ c^R_⋆ρ r^2 [(l+1)(l+2)r^2 |ĉ_l+1^m|^2] dr.

Thus inserting Eqs. (<ref>) and (<ref>) into Eq. (<ref>) and (<ref>), we get 

    Ê_l^m   =   0,
    Ê_l-1^m   =   l (l-1) (2l+1) (q̃_l^m)^2  Ω^2 R_⋆
          [ρ_ cα^2l+1 |C_1|^2 + ρ_ e (1- α^2l+1) |C_2|^2],
    Ê_l+1^m   =    (l+1)(l+2) (2l+1) (q̃_l+1^m)^2 Ω^2 R_⋆
          ρ_ eα^2l+1 (1-α^2l+1) |C_3|^2 ,

which are readily expressed in terms of stellar quantities using Eqs. (<ref>) – (<ref>):

    Ê_l-1^m=l (l-1) (2l+1) (q̃_l^m)^2 Ω^2  R_⋆(A_l/4 π G l ρ_ c)^2 ( 1+p_l/p_l-q_l)^2 
    [ρ_ cα^2l+1(2l+1/2l-2/3+f)^2 + ρ_ e(1- α^2l+1) (1/f+(1-f)α^3)^2 .
    .×(3 - (f/2l-2/3+f)(α^2l+1 + (2l+1)(1-α^3)/1-α^2l+1))^2 ],



    Ê_l+1^m= (l+1)(l+2) (2l+1 (q̃_l+1^m)^2 Ω^2  R_⋆
    ρ_ eα^2l+1(A_l/4 π G (l+1) ρ_ c)^2 ( 1+p_l/p_l-q_l)^2 
    (1-f )^2/1-α^2l+1(3α^2l+1 - (2l+1)α^3 + 2l-2/(2l-2/3+f) (f+(1-f)α^3))^2 .


It can be shown that for a tidal potential of the form Ψ =  Re[A_l(r/R_⋆)^l Y_l^m (θ, ϕ) e^-i ω t], the rate of energy transfer given by P=∫ρξ̇·( - ∇Ψ)  dV
is related to the potential Love numbers through

    P = (2l+1) R_⋆/ 4 π G1/2|A_l|^2ω Im[k_l^m(ω)].


By analogy, in the Fourier-transformed problem, again assuming that the force derives from a tidal potential of the form Ψ =  Re[A_l(r/R_⋆)^l Y_l^m (θ, ϕ) e^-i ω t], we define the complex dimensionless response functions K_l^m(ω) as

    K_l^m(ω) ≡4 π G/|A_l|^2 R_⋆ (2l+1)∫ρξ̃_̃w̃·f̅̃̅  dV.

So that the total energy transferred to the wavelike tide is, using the same analogy,

    Ê= 2 (2l+1)R_⋆/4 π G∫_-∞^∞ω Im[K_l^m(ω)] | Ψ̃_l^m (ω)|^2   dω/ω.

In the case of an impulsive tidal potential, the energy transfer is simply

    Ê= 2 (2l+1)R_⋆/4 π GA_l^2/8π∫_-∞^∞ Im[K_l^m(ω)]  dω/ω,

and using Eq. (<ref>), we finally get that

    ∫_-∞^∞ Im[K_l^m(ω)]  dω/ω = Ê_l-1^m+Ê_l+1^m /(2l+1) R_⋆ G(4π G/A_l)^2,

so that

    ∫_-∞^∞ Im[K_l^m(ω)]  dω/ω = 4/3(l-1/l)  ϵ^2   ( 1+p_l/p_l-q_l)^2  ( f/f+(1-f)α^3)   
    [(q̃_l^m)^2 (2l+1/2l-2/3+f)^2 (f+(1-f)α^3)^2 (α^2l+1/f) . 
    
    + (q̃_l^m)^2  (1- α^2l+1) (3 - (f/2l-2/3+f)(α^2l+1 + (2l+1)(1-α^3)/1-α^2l+1))^2 
    
    + (q̃_l+1^m)^2 (l/l-1)(l+2/l+1) (1-f/2l-2/3+f)^2 (α^2l+1/1-α^2l+1) 
    . (3α^2l+1 - (2l+1)α^3 + 2l-2 )^2 ],

where ϵ is a dimensionless parameter measuring stellar rotation normalised to the breakup velocity

    ϵ = Ω_⋆(G M_⋆/R_⋆)^-1/2.


Lastly, to the level of accuracy of our approximations, we have Im[K_l^m(ω)] =   Im[k_l^m(ω)]. So in this way we get an expression of ∫^∞_-∞ Im[k_l^m(ω)]  dω/ω that depends only on the angular velocity and the masses and radii of the radiative core and the convective envelope, valid for arbitrary values of l and m. 



§ STATIC LOVE NUMBERS

We stress here that the Helmholtz-like equation (<ref>) for the non-wavelike part in the Fourier-transform problem, neglecting the effect of the Coriolis force, is exactly equivalent to the one that would be obtained by combining Newton's second law of action and Poisson's equation, for a non-rotating fluid in hydrostatic equilibrium experiencing the static perturbing gravitational potential from a point-like companion. The general form of the solutions given in Eq. (<ref>) also describe the gravitational potential arising from the mass of the tidally distorted star in hydrostatic equilibrium at any point r in space.

Indeed, as it is well-known, in a frame of reference centred on the star, the gravitational potential produced by a point-like companion of mass M_p placed at r_p=(r_p, θ', ϕ') with r_p>R_⋆ can be expanded over the Legendre polynomials. The tidal potential Ψ is simply the potential from which derives the variable part of the perturbing force, so 


    Ψ(r, θ,ϕ) =  G M_p/r_p∑_l = 2^∞(r/r_p)^l P_l(cosγ),

where  γ denotes the angular separation between the vectors r and  r_p <cit.>. We have, by the addition theorem for spherical harmonics,

    P_l(cosγ) = 4π/2l+1∑_m=-l^l Y_l^m(θ, ϕ)Y̅_l^m(θ', ϕ').


In hydrostatic equilibrium, any surface of a fluid body characterised by equal density and equal pressure corresponds to a surface of equal total potential of forces acting upon the body.  Because of the tidal force, the star in hydrostatic equilibrium cannot assume a perfectly spherical shape. Following <cit.>, without any loss of generality, we can assume that on such a surface the radial distance r can be written as 

    r = r_*( 1 + ∑_l=1^∞∑_m=-l^l  S_l^m(r_*)Y_l^m (θ, ϕ) ),

where the Y^m_l(θ, ϕ) are the normalised spherical harmonics following the Condon-Shortley phase convention, and S_l^m(r_*) are complex functions of the mean radius r_* that must verify

    S_l^-m =(-1)^m S̅_l^m,

for r to be a real quantity. By virtue of the uniqueness of the potential function, and because the density must remain constant over an equipotential, the density can be regarded as a function of the single variable r_*, denoting the mean radius of the corresponding equipotential. Thus, for an arbitrary point (r,θ,ϕ) in the interior of the star,  we can write the interior self-potential U of this distorted configuration  as 

    U(r,θ,ϕ) = U_0 + 
    ∑_l=1^∞4π G r^l/2l+1∑_m=-l^l Y_l^m(θ, ϕ)  ∫_r_*^R_⋆ρ∂ (r_*^2-l S_l^m(r_*))/∂ r_*dr_*,

where ρ is the density, R_⋆ denotes the mean radius at the surface of the body, and we have introduced

    U_0 =  4 π G ∫_r_*^R_⋆ρ r_* dr_*.

In the same way, the exterior potential V(r, θ, ϕ) can be expressed as

    V(r,θ,ϕ) = V_0/r + 
    ∑_l=2^∞4π G/(2l+1)r^l+1∑_m=-l^l Y_l^m(θ, ϕ)  ∫_0^r_*ρ∂ (r_*^l+3 S_l^m(r_*))/∂ r_*dr_*

with 

    V_0 =  4 π G ∫_0^r_*ρ r_*^2 dr_*.


Now let Φ  be the total potential of forces acting on any arbitrary point. In hydrostatic equilibrium, the function Φ is to remain constant over any surface of equal pressure of density, over which r verifies Eq. (<ref>). By binomial expansion of r, neglecting squares and higher order of Y_l^m, the shape of an equipotential is described by the spherical harmonic expansion with coefficients S_l^m that are solution of the following equation

    (2l+1) S_l^m(r_*) ∫_0^r_*ρ r_*^2 dr_* - r_*^l+1∫_r_*^R_⋆ρ∂ (r_*^2-l S_l^m(r_*))/∂ r_*dr_*  
    
    - 1/r_*^l∫_0^r_*ρ∂ (r_*^l+3 S_l^m(r_*))/∂ r_*dr_* =  M_p  (r_*/r_p)^l+1Y̅_l^m(θ', ϕ'),

for any l≥1 and -l≤ m≤ l. This is the well-known Clairaut's equation. It implicitly specifies the spherical harmonics describing the form of equipotential level surfaces as distorted by an external force derived from the potential (<ref>). Solving Clairaut's equation allows us to derive the explicit form of the radial distance r_S describing the free surface of the perturbed body of mass M_⋆ in terms of the mass M_p of the perturber 

    r_S = R_⋆( 1 + ∑_l=1^∞M_p/M_⋆(R_⋆/r_p)^l+12l +1/ l+ η_l (R_⋆) P_l(cosγ) ),

where η_l(R_⋆) is the surface value of the logarithmic derivative of the shape coefficients

    η_l (r_*) = r_*/S_l^m(r_*)∂ S_l^m(r_*)/∂ r_*.


The potential produced by this distorted configuration at any point (r, θ, ϕ), with r >> R_⋆ using the center of mass of the extended body as the origin, is generally written as

    Φ(r,θ,ϕ) = G M_⋆/r[1 + ∑_l=2^∞ k_l M_p/M_⋆(R_⋆/r_p)^l+1(R_⋆/r)^l  P_l(cosγ)   ]

where k_l are called the static Love numbers <cit.>. Using equations Eqs. (<ref>) and (<ref>), the static Love numbers are found to be

    k_l =  l+1 - η_l (R_⋆)/ l+ η_l (R_⋆).
 
Thus the Love numbers link the external potential to the internal structure through the surface value of the logarithmic derivative of the shape coefficients. 

For certain density distributions, Clairaut's equation admit simple solutions in a closed form and the Love number can take explicit forms. For example, we have the well-known result

    k_l = 3/2l-2

if the density remains uniform throughout the interior of the body, and 

    k_l = l+2/l-1

if the whole mass of the configuration were confined to an infinitesimally thin surface shell. If the whole mass of the configuration is condensed at its centre, obviously, k_l = 0 for all l , and the exterior potential (<ref>) is that of a point mass.


Now in this paper, we consider a two-layer model that simplifies the actual density distribution within the star. To our knowledge, the expression of the Love numbers in this case is not easily available in recent literature so we provide it here for convenience, since it is a direct by-product of our derivations. Since we consider that ρ_ c and ρ_ e are uniform within their respective layers, Clairaut's equation, taken at the core boundary,  simplifies to

    ( 2l-2/3 + ρ_ e/ρ_ c) ρ_ c  R_ c^3  S_l^m(R_ c) - (R_ c/R_⋆)^l+1ρ_ e R_⋆^3 S_l^m(R_⋆)
    
     =  M_ p(R_ c/r_p)^l+1Y̅_l^m(θ', ϕ'),

while it gives at the free surface

    ( 2l+1/3(ρ_ c/ρ_ e -1 ) (R_ c/ R_⋆)^3  + 2l-2/3)  ρ_ e R_⋆^3S_l^m(R_⋆)  
    
    +(ρ_ e/ρ_ c -  1 ) (R_ c/R_⋆)^lρ_ c R_ c^3 S_l^m(R_ c) 
     =  M_ p(R_⋆/r_p)^l+1Y̅_l^m(θ', ϕ').

Using those two equations to solve for S_l^m(R_ c) and S_l^m(R_⋆) and letting f=ρ_ e/ρ_ c, we find

    S_l^m(R_⋆)    =   M_ p/ρ_ e R_⋆^3(R_⋆/r_p)^l+11+ q_l/p_l - q_l Y̅_l^m(θ', ϕ'),
    
    S_l^m(R_ c)    =   M_ p/ρ_ c R_ c^3(R_ c/r_p)^l+11/2l-2/3 + f1+ p_l/p_l - q_l Y̅_l^m(θ', ϕ'),

where p_l and q_l are defined in Eq. (<ref>) and  (<ref>), respectively. Those expressions are equivalent to the ones derived by <cit.> for a two-layer planet. They considered the case of a purely elastic solid core, which is equivalent to the fluid case when setting the shear modulus to zero. For the two-layer model considered here, the potential at any point (r, θ, ϕ) with r >> R_⋆ produced by the distorted configuration can be then written in the form

    Φ(r,θ,ϕ) =   G M_⋆/r +∑_l=2^∞∑_m=-l^l B_1,  eq(R_⋆/r)^l+1  Y_l^m(θ, ϕ) 
    
    +∑_l=2^∞∑_m=-l^l B_2,  eq(R_ c/R_⋆)^2l+1(R_⋆/r)^l+1    Y_l^m(θ, ϕ),

where we have abbreviated

    B_1, eq = 4π/(2l+1)G M_ p/R_⋆(R_⋆/r_p)^l+11+q_l/p_l-q_lY̅_l^m(θ', ϕ')


    B_2,  eq = 4π/(2l+1)G M_ p/R_⋆(R_⋆/r_p)^l+11+p_l/p_l-q_l( 1- f ) /(2l-2/3 +f )Y̅_l^m(θ', ϕ').


Finally the potential produced by the distorted star Eq.(<ref>) can also be rearranged to take the form of Eq. (<ref>) where the static Love number now has the explicit expression

    k_l = (2+p_l) q_l +1/p_l-q_l.

This yields the expected k_l = 3/2l-2 in the limit case of a homogenous body with ρ_ e = ρ_ c.

Note that the correspondence with Eqs. (<ref>) and (<ref>) is easily made by considering that, in the hydrostatic case, the tidal potential expressed as Eq. (<ref>) yields

    A_l ≡4π/(2l+1)G M_ p/R_⋆(R_⋆/r_p)^l+1Y̅_l^m(θ', ϕ'),

with A_l a constant since the perturber is considered here to be static.

In the dynamical case, if one neglects the Coriolis force, the reasoning presented here still holds, but the coordinates of the perturber now must assume a temporal dependence . Therefore, the non-wavelike part of the response (which is equivalent to the instantaneous adjustment of the star when neglecting the Coriolis force) can easily be solved by considering an adequate setting where the temporal dependence  of the coordinates of the perturber is removed and A_l is indeed a constant. This is done in <cit.> by using a frame of reference that rotates with the star, considering a tidal forcing of impulsive nature, and taking the problem into the Fourier domain <cit.>. For the circular misaligned problem, the same method can be used. In this case, the radial coordinate of the companion is constant over an orbit with r_p = a. Its latitudinal coordinate θ' intervenes only through the spherical harmonics Y̅_l^m(θ', ϕ'), which can be expressed as a function of the obliquity Θ by a simple rotation of the axis of the reference frame (hence the Wigner-d matrix coefficient in Eq. <ref>). Thus, the temporal dependence  only persists in the exponential term of Y̅_l^m(θ', ϕ'), which contributes to the harmonic forcing, and defines the tidal frequency in the rotating frame.

