© ESO, 2017

1. Introduction

Since the pioneering work of Kelvin (1863) on the gravitationally forced elongation of a solid spherical shell, tides have been recognized as a phenomenon of major importance in celestial mechanics. They are actually one of the key mechanisms involved in the evolution of planetary systems over long timescales. Because they affect both stars and planets, tides introduce a complex dynamical coupling between all the bodies that compose a system. This coupling results from the mutual interactions of bodies and is tightly bound to their internal physics. Solid tides have long been the subject of studies (see the reference works by Love 1911; Goldreich & Soter 1966). The corresponding tidal dissipation, which is caused by viscous friction, can have a strong impact on the orbital elements of a planet (e.g. Mignard 1979, 1980; Hut 1980, 1981; Neron de Surgy & Laskar 1997; Henning et al. 2009; Remus et al. 2012). This evolution depends smoothly on the tidal frequency (e.g. Efroimsky & Lainey 2007; Greenberg 2009). Fluid tides cause internal dissipation as well as solid tides, and the energy dissipated is even more dependent on the tidal frequency than the solid tides are because of a highly resonant behaviour (Ogilvie & Lin 2004; Ogilvie 2014; Auclair Desrotour et al. 2015). As a consequence, in both cases the phase shifted elongation of the shell induces a tidal torque which modifies the rotational dynamics of the body. This torque drives the evolution of the spin of planets and determines its possible states of equilibrium (Gold & Soter 1969; Dobrovolskis & Ingersoll 1980; Correia & Laskar 2001; Correia et al. 2003, 2014; Correia & Laskar 2003; Arras & Socrates 2010).

The increased focus on exoplanets over the past two decades and the rapidly increasing number of discovered planetary systems(Mayor & Queloz 1995; Perryman 2011) have brought to light a huge diversity of existing orbital structures that can be very different from our own solar system (Fabrycky et al. 2012). Therefore, tidal effects now need to be studied in a general context to characterize these new systems, understand their orbital and rotational dynamics, and constrain the physics of exoplanets with the information provided by observational measurements.

In this work, we focus on Earth-like exoplanets, which are interesting for their complex internal structure and possible habitability. Indeed, a planet such as the Earth or Venus is composed of several solid and fluid layers with different physical properties, which implies that the tidal response of a layer cannot be simply correlated to its size. For instance, the Earth’s oceans, in spite of their very small depth compared to the Earth’s radius, are responsible for most of the tidal dissipation of the planet (Webb 1980; Egbert & Ray 2000; Ray et al. 2001). The atmospheric envelope of a terrestrial exoplanet sometimes represents a non-negligible fraction of its mass (about 20% in the case of small Neptune-like bodies). It is also a layer of great complexity because it is submitted to both tidal gravitational potentials due to other bodies and thermal forcing of the host star. Therefore, it is necessary to develop theoretical models of atmospheric tides with the smallest possible number of parameters. These models will be useful tools to explore the domain of parameters, and to quantify the tidal torque exerted on the spin axis of a planet and quantities used in the orbital evolution of planetary systems. One of these quantities is the second-order Love number (k₂) introduced by Love (1911), which measures the quadrupolar hydrostatic elongation. The adiabatic Love number k₂ has been estimated in the solar system (e.g. Konopliv & Yoder 1996; Lainey et al. 2007; Williams et al. 2014).

According to Wilkes (1949), the first theoretical global models of atmospheric tides were developed at the end of the eighteenth century by Laplace, who published them in Mécanique céleste (Laplace 1798). Laplace was interested in the case of a homogeneous atmosphere characterized by a constant pressure scale height and affected by tidal gravitational perturbations. He thus founded the classical theory of atmospheric tides. At the end of the nineteenth century, Lord Kelvin pointed out the unexpected importance of the Earth’s semidiurnal pressure oscillations in spite of the predominating diurnal perturbation (Kelvin 1882; Hagan et al. 2003; Covey et al. 2009, see also the measures of the daily variation of temperature and pressure, Fig. 13). This question motivated many further studies (e.g. Lamb 1911, 1932; Taylor 1929, 1936; Pekeris 1937) which contributed to enriching the classical theory of tides. In the 1960s, theoretical models were generalized in order to study other modes, such as diurnal oscillations (see Haurwitz 1965; Kato 1966; Lindzen 1966, 1967a,b, 1968), but still remained focused on the case of the Earth. They revealed the impact of thermal tides due to solar insolation by showing that these tides drive the Earth’s diurnal and semi-diurnal waves.Indeed, in the Earth atmosphere, the contribution of the tidal gravitational potential, which causes the elongation of solid layers, is negligible compared to the effect of the thermal forcing. A very detailed review of the general theory of thermal tides and of its history is given by Chapman & Lindzen (1970), hereafter CL70 (see also Wilkes 1949; Siebert 1961).

Leaning on this early work, here we generalize the existing theory to Earth-like planets and exoplanets. We introduce the mechanisms of radiative losses and thermal diffusion, which are usually not taken into account in geophysical models of atmospheric tides owing to their negligible effects on the Earth’s tidal response (radiation was, however, introduced with a Newtonian cooling by Lindzen & McKenzie 1967; Dickinson & Geller 1968, which inspired the present work). These mechanisms cannot be ignored because they play a crucial role in the case of planets like Venus, which are very close to spin-orbit synchronization. Models based on perfect fluid hydrodynamics, such as the one detailed by CL70, fail to describe the atmospheric tidal response of such planets because they are singular at synchronization. In these models, the amplitude of perturbed quantities like pressure, density, and temperature diverges at this point. The case of Venus was studied by Correia & Laskar (2001), who proposed an empirical bi-parameter model for the tidal torque describing two possible regimes: the regime of an atmosphere where dissipative processes can be ignored, far from synchronization, and its supposed regime at the vicinity of synchronization, where the perturbation is expected to annihilate itself. This behaviour was recently retrieved with global climate models (GCM) simulations by Leconte et al. (2015), who took into account dissipative mechanisms, thus showing that these mechanisms drive the response of planetary atmospheres in the neighbourhood of synchronization.

Our objective in this work is thus to propose a new theoretical global model controlled by a small number of physical parameters which could be used to characterize the response of any exoplanet atmosphere to a general tidal perturbation. We aim to answer the following questions:

• How does the amplitude of pressure, density, temperature, andvelocity oscillations depend on the parameters of the system (i.e.rotation, stratification, radiative losses)?

• What are the possible regimes for atmospheric tides?

• How do Love numbers and tidal torques vary with tidal frequency?

First, in Sect. 2, we derive the general theory for thick atmospheres, which are characterized by a depth comparable to the radius of the planet. In this framework, we expand the perturbed quantities in Fourier series in time and longitude, and in series of Hough functions (Hough 1898) in latitude. This allows us to compute the equations giving the vertical structure of the atmospheric response to the tidal perturbation. Then, from the hydrodynamics we derive the analytic expressions of the variations of pressure, density, temperature, velocity field, and displacement for any mode. At the end of the section, the density oscillations obtained are used to compute the atmospheric tidal gravitational potential, the Love numbers, and the tidal torque exerted on the spin axis of the planet. In Sect. 3, we apply the equations from Sect. 2 to the case of a thin isothermal stably stratified atmosphere. In this simplified case, the terms associated with the curvature of the planet disappear and most of the parameters depending on the altitude, such as the pressure scale height of the atmosphere, become constant. In Sect. 4, we treat the case of slowly rotating convective atmospheres, such as near the ground in Venus (Marov et al. 1973; Seiff et al. 1980). In the next section, we derive the terms of thermal forcing from the Beer-Lambert law (Bouguer 1729; Klett et al. 1760; Beer 1852) and the theory of temperature oscillations at the planetary surface. In Sect. 6, the model is applied to the Earth and Venus. The corresponding torque and atmospheric Love numbers are computed from analytical equations. We end the study with a discussion and by giving our conclusions and prospects. To facilitate the reading, several technical issues have been deferred to the Appendix where an index of notations is also provided.

2. Dynamics of a thick atmosphere forced thermally and gravitationally

The reference book CL70 has set the basis of analytical approaches for atmospheric tides. This pioneering work focuses on fast rotating telluric planets covered by a thin atmosphere, such as the Earth.

The goal of the present section is to generalise this work. We establish the equations that govern the dynamics of tides in a thick fluid shell in the whole range of tidal frequencies, from synchronization to fast rotation. We consider a spherical telluric planet of radius R covered by a stratified atmospheric layer of typical thickness H_atm. The atmosphere rotates uniformly with the rocky part at the angular velocity Ω (the spin vector of the planet being denoted Ω). Therefore, the dynamics are written in the natural equatorial reference frame rotating with the planet, R_E;T:{O,X_E,Y_E,Z_E}, where O is the centre of the planet, X_E and Y_E define the equatorial plane and Z_E = Ω/ |Ω|. We use the spherical basis and coordinates (r,θ,ϕ), where r is the radial coordinate, θ the colatitude, and ϕ the longitude. As usual, t stands for the time.

The atmosphere is assumed to be a perfect gas homogeneous in composition, of molar mass M, and stratified radially. Its pressure, density, and temperature are denoted p, ρ, and T respectively. For the sake of simplicity, all dissipative mechanisms (e.g. viscous friction and heat diffusion) are ignored except for the radiative losses of the gas, which play an important role in the vicinity of synchronization. Tides are considered to be a small perturbation around a global static equilibrium state. The basic pressure p₀, density ρ₀, and temperature T₀, are supposed to vary only with the radial coordinate. The tidal response is characterized by a combination of inertial waves, of typical frequency 2Ω, due to Coriolis acceleration, and gravity waves which are restored by the stable stratification. The typical frequency of these waves, the Brunt-Väisälä frequency (e.g. Lighthill 1978; Gerkema & Zimmerman 2008), is given by (1)where Γ₁ = (∂lnp₀/∂lnρ₀)_S is the adiabatic exponent (the index S being the specific macroscopic entropy) and g the gravity. We assume that the fluid follows the perfect gas law, (2)where R_s = R_GP/M is the specific gas constant of the atmosphere, R_GP = 8.3144621 J mol^-1 K^-1 (Mohr et al. 2012) being the molar gas constant. In the following, the notations ℜ, ℑ, and ^∗ are used to designate the real part, imaginary part, and conjugate of a complex number. The subscript symbols _⊕ and _♀ stand for the Earth and Venus, respectively, which are taken as examples to illustrate the theory. All functions and fields displayed in the figures are computed with the algebraic manipulator TRIP (Gastineau & Laskar 2014) and plotted with gnuplot.

2.1. Forced dynamics equations

The atmosphere is affected by the tidal gravitational potential U and the heat power per mass unit J. In this approach, the equations of dynamics are linearized around the equilibrium state. The small variations in this state induced by the tidal perturbation are the velocity field V = V_θe_θ + V_ϕe_ϕ + V_re_r; the corresponding displacement field ξ = ξ_θe_θ + ξ_ϕe_ϕ + ξ_re_r defined as (3)and the pressure (δp), density (δρ), and temperature (δT) fluctuations where (4)This represents six unknown quantities to compute (V_r, V_θ, V_ϕ, δp, δρ, δT), and we therefore solve a six-equation system. We first introduce the Navier-Stokes equation. We assume the Cowling approximation in which the self-gravitational potential fluctuations induced by tides are neglected (Cowling 1941). It is well adapted to waves characterized by a high vertical wavenumber as are found in the vicinity of synchronization, as detailed farther. In the equatorial reference frame R_E;T, the linearized Navier-Stokes equation is thus written (5)which, projected on e_r, e_θ, and e_ϕ, gives (6)(7)(8)Equations (6)–(8) are coupled by the Coriolis terms (characterized by the factor 2Ω) on the left-hand side. To simplify the equations, we assume the traditional approximation, i.e. to neglect the terms 2ΩsinθV_r and 2ΩsinθV_ϕ in Eqs. (7) and (8), respectively. This hypothesis, commonly used in literature dealing with planetary atmospheres and star hydrodynamics (e.g. Eckart 1960; Mathis et al. 2008), can be applied to stratified fluids where the tidal flow satisfies two conditions: (i) the buoyancy (given by gδρ/ρ₀e_r in Eq. (5)) is strong compared to the Coriolis term (2Ω ∧ V) and (ii) the tidal frequency (introduced below) is greater than the inertia frequency (2Ω) and smaller than the Brunt-Väisälä frequency. These conditions are satisfied by the Earth’s atmosphere and fast rotating planets. The quantitative precision that can be expected with the traditional approximation may be lower in the vicinity of synchronization where the tidal frequency tends to zero; moreover, we note that it leads to issues with momentum and energy conservation in the case of deep atmospheres (for a discussion on the limitations of the traditional approximation, see Gerkema & Shrira 2005; White et al. 2005; Mathis 2009).

Hence, we obtain: (9)(10)(11)In strongly stratified fluids, the left-hand side of Eq. (11) is usually ignored because it is very small compared to the terms in δp and δρ in typical wave regimes. The radial acceleration will only play a role in regimes where the tidal frequency is close to or exceeds the Brunt-Väisälä frequency, which corresponds to fast rotators.

The second equation of our system is the conservation of mass, (12)which, in spherical coordinates, writes (13)The thermal forcing (J) appears on the right-hand side of the linearized heat transport equation (see CL70, Gerkema & Zimmerman 2008) given by (14)where and J_rad is the power per mass unit radiated by the atmosphere, supposed to behave as a grey body. We consider that J_rad∝δT. This hypothesis is known as “Newtonian cooling” and was used by Lindzen & McKenzie (1967) to introduce radiation analytically in the classical theory of atmospheric tides (see also Dickinson & Geller 1968). Physically, it corresponds to the case of an optically thin atmosphere in which the flux emitted by a layer propagates upwards or downwards without being absorbed by the other layers. In optically thick atmospheres, such as on Venus (Lacis 1975), this physical condition is not verified. Indeed, because of a stronger absorption, the power emitted by a layer is almost totally transmitted to the neighbourhood. Therefore, this significant thermal coupling within the atmospheric shell should ideally be taken into account in a rigorous way, which would lead to great mathematical difficulties (e.g. complex radiative transfers, Laplacian operators) in our analytical approach. However, recent numerical simulations of thermal tides in optically thick atmospheres by Leconte et al. (2015) show behaviour of the flow that is in good agreement with a model that uses radiative cooling. Therefore, in this work we assumeNewtonian cooling as a first model of the action of radiation on atmospheric tides. Newtonian cooling brings a new characteristic frequency, denoted σ₀, which we call radiative frequency and which depends on the thermal capacity of the atmosphere. The radiative power per unit mass is thus written (15)Like the basic fields p₀, ρ₀, and T₀, the radiative frequency varies with r and defines the transition between the dynamical regime, where the radiative losses can be ignored, and the radiative regime, where they predominate in the heat transport equation. Assuming that the radiative emission of the gas is proportional to the local molar concentration C₀ = ρ₀/M, it is possible to express J_rad and σ₀ as functions of the physical parameters of the fluid (Appendix D), (16)and (17)the parameter ϵ_a being an effective molar emissivity coefficient of the gas and S= 5.670373 × 10^-8 W m^-2 K^-4 the Stefan-Boltzmann constant (Mohr et al. 2012). The substitution of Eq. (15) in Eq. (14) yields (18)Finally, the system is closed by the perfect gas law (19)Substituting Eq. (19) in Eq. (18), we eliminate the unknown δT, and obtain (20)Because of the rotating motion of the perturber in the equatorial frame (R_E;T), a tidal perturbation is supposed to be periodic in time (t) and longitude (ϕ). So, any perturbed quantity f of our model can be expanded in Fourier series of t and ϕ(21)the parameter σ being the tidal frequency of a Fourier component and m its longitudinal degree¹. We also introduce the spin parameter (22)which defines the possible regimes of tidal gravito-inertial waves:

• |ν| ≤ 1 corresponds to super-inertial waves, for which the tidal frequency is greater than the inertia frequency;

• |ν| > 1 corresponds to sub-inertial waves, for which the tidal frequency is lower than the inertia frequency.

The parameters m and ν determine the horizontal operators

Here the operator is associated with and , while is associated with and (Lee & Saio 1997; Mathis et al. 2008). Hence, the traditional approximation allows us to write the horizontal component of the velocity field only as a function of the variations in the pressure and tidal gravitational potential

where y^m,σ = δp^m,σ/ρ₀. Since V^m,σ = ∂_tξ^m,σ, we also get the corresponding horizontal component of the displacement vector from Eqs. (25) and (26)

The substitution of Eqs. (25) and (26) in Eq. (13) yields (29)where ℒ^m,ν is the Laplace operator, parametrized by ν and m and is expressed by (30)with ν = 2Ω /σ (Eq. 22). When ν = 0 (i.e. Ω = 0), ℒ^m,ν is reduced to the horizontal Laplacian. Since the horizontal velocity field is not coupled with the other variables, only three unknowns remain. So, the new system is written

Under the traditional approximation, solutions of Eqs. (31)−(33) can be obtainedunder the form of series of functions of separated coordinates,

(34)where the functions are the eigenvectors of the Laplace operator, the corresponding eigenvalues being denoted (with ν = 2Ω /σ). Laplace’s tidal equation, expressed as (35)describes the horizontal structure of the perturbation. The solutions of Eq. (35), are called Hough functions (Hough 1898)². They form a complete set of continuous functions on the interval θ ∈ [0,π], parametrized by m and ν. They also verify the same boundary conditions as the associated Legendre polynomials (denoted ), which are solutions of Eq. (35) for ν = 0 with . In Figs. 2 and 3, the functions , , and are plotted for the Earth and Venus in the case of the solar semidiurnal tide defined by m = 2, σ = 2(Ω−n_orb), and ν = Ω/(Ω−n_orb), where n_orb stands for the mean motion of the planet. In the case of the Earth, the set is mainly composed of gravity modes. Since the spin parameter is slightly greater than 1, Rossby modes exist but are trapped at the poles. In the case of Venus, Hough functions are composed of gravity modes only because |ν| < 1. These functions are very similar to the associated Legendre polynomials , with l ≥ 2.

The system composed of Eqs. (31) to (33) can then be reduced to the pair of first-order differential equations (36)with the coefficients (37)and (38)where we introduce the sound velocity (39)The system of Eq. (36) can itself be reduced to a single equation in alone (40)with the set of coefficients (41)Equation (40) gives us the vertical structure of tidal waves generated in the fluid shell by both gravitational and thermal forcings. We have: (42)The last term, (dB₁/ dr) /B₁, is associated with curvature and can be ignored in the context of the “shallow atmosphere” approximation, as discussed in the next section. It is expressed as (see Eq. (50)) (43)We now transform Eq. (40) in a Schrödinger-like equation by introducing the radial function (44)and the variable change (45)where the volume is the wave function of tidal waves associated with the mode of degree n. We note that fixing σ₀ = 0 allows us to retrieve the usual variable change, (Press 1981; Zahn et al. 1997; Mathis et al. 2008). With the general variable change, Eq. (40) becomes (46)which gives us the vertical profile of the perturbation³. The form of this equation is very common in the literature because it describes the radial propagation of waves in a spherical shell. The typical wavenumber of these waves, denoted , is given by (47)We introduce the corresponding normalized length scale of variations of the mode of degree n, (48)This parameter is the typical height over which spatial variations of the mode can be observed.

2.2. Polarization relations

The perturbed quantities are readily deduced from Ψ. Before going further, let us introduce the Lamb frequency (the cut-off frequency of acoustic waves) associated with the mode (ν,m,n),(49)and the ratio (50)The parameter |ε_s,n| measures the relevance of the anelastic approximation, which consists in neglecting all the terms that correspond to an acoustic perturbation in Eq. (36). When |ε_s,n| ≪ 1 (i.e. |σ| ≪ |σ_s;n|), this hypothesis can be assumed. However, particular caution should be used when ignoring some terms. For instance, since the first gravity mode of the Earth’s semidiurnal tide is characterized by |σ| ≈ |σ_s;0| ≈ 10^-4 s^-1, the anelastic approximation cannot be done. In order to be able to study cases of the same kind, the polarization relations that follow will be given in their most general form. Hence, by substituting the series of Eq. (34) in Eqs. (25)–(28), and (31)–(33), we obtain the radial profiles , , , , , , , , , and as functions of and its first derivative. We obtain for the Lagrangian displacement

for the velocity

and for scalar quantities (57)

The coefficients in these expressions are given by

(60)where we have introduced the curvature term (61)which will be ignored in the thin atmosphere approximation. These relations and Eq. (46) define a linear operator, denoted Υ^σ, which determines the vertical structure of the atmospheric tidal response. Let be the output vector of the model and the input vector. Eqs. (51) to (59) can be reduced to (62)For any quantity f, the component of describing the vertical profile of f is denoted ( for the pressure, for the density, etc.).

Polarization relations can be naturally divided into two distinct contributions. The first contribution, which is called the “horizontal part” in this work, is a linear combination of the forcings and . It corresponds to the component of the response that does not depend on the tidal vertical displacement. The length scale of this horizontal component is the typical thickness of the atmosphere (H_atm). The other part, expressed as a function of and called “vertical part”, thus results from the tidal fluid motion along the vertical direction. In the following, we use the subscripts _H and _V to designate the horizontal and vertical components.

As can be noted, polarization relations point out the necessity of taking into account dissipative processes to study tidal regimes close to synchronization. Indeed, without dissipation (σ₀ = 0), synchronization (σ = 0) is characterized by a singularity. The terms associated with the horizontal component, which are directly proportional to and , tend to infinity at σ → 0. Terms in and behave similarly, as demonstrated in the next section. These terms are highly oscillating along the vertical direction over the whole depth of the atmosphere in this frequency range because when σ → 0. Hence, radiation regularizes the behaviour around σ = 0 and damps the vertical oscillations of waves associated with the vertical component. Diffusion acts in the same way by flattening the oscillations of smallest length scales (L_V;n) (Press 1981; Zahn et al. 1997; Mathis et al. 2008). We return to this point when we solve analytically the case of the isothermal atmosphere with constant profiles for the forcings (Sect. 3).

The whole spectrum of possible tidal regimes is represented in Fig. 4. The radiative regime (blue) characterizes a flow where thermal losses due to the linear Newtonian cooling (see the heat transport equation, Eq. 18) predominate. In the dynamic regime (orange), dissipation can be neglected because the system is governed by the Coriolis acceleration. Finally, the mixed-acoustic regime (red) corresponds to high tidal frequencies comparable to the Lamb frequency, given by Eq. (49). This regime marks the limit of the traditional approximation.

2.3. Tidal potential, Love numbers, and tidal torque

The new mass distribution resulting from tidal waves generates a tidal gravitational potential , which is a linear perturbation of the spherical gravitational potential of the planet. This atmospheric potential presents disymmetries with respect to the direction of the perturber. The tidal torque thus induced will affect the rotationnal dynamics of the planet over secular timescales and determine the possible equilibrium states of the spin, as demonstrated by Correia & Laskar (2001). This torque is deduced from the Poisson’s equation expressed in the co-rotating frame (63)the notation G designating the gravitational constant. Like δρ, the potential is expanded in series of functions of separated variables (64)where the (with l ∈N such as l ≥ |m|) are the normalized associated Legendre polynomials. Decomposing δρ^m,σ on this basis, we get the equation describing the vertical structure of , (65)For the upper boundary condition, the tidal potential is required to remain bounded at r → ∞. At r = 0, we impose the same condition. Therefore, the solution of Eq. (65) is (66)where and are the functions (67)Love numbers are defined as the ratio between the gravitational potential due to the tidal response of the atmosphere and the forcing potential at the upper boundary of the layer. Let us denote R_atm = R + H_atm this upper boundary. Then, at the upper boundary, the atmospheric tidal potential given by Eq. (66) is simply expressed: (68)Considering Eq. (72) and using the linearity property of , the density component can be expanded (69)with = and the change-of-basis coefficients, (70)which allows us to write the atmospheric tidal potential (71)Let us treat the case of a simplified planet-star system, where the planet P circularly orbits around its host star, of mass M_S. The semi-major axis, obliquity, and mean motion of the planet are denoted a, ε, and n_orb. In this case, the modes of the gravitational and thermal forcings applied on the atmosphere can be written (Appendix C) (72)Hence the complex Love numbers can be written generically (73)In celestial dynamics, the second-order Love number (l = 2) is commonly used to quantify the tidal response of a body. Therefore, we illustrate the expression given by Eq. (73) by computing this coefficient for m = 2. For the sake of simplicity, we set the obliquity ε = 0 (the spin of the planet is supposed to be aligned with its orbital angular momentum). The tidal frequency is σ_SD = 2(Ω−n_orb) and the series of Eq. (C.14) are reduced to the terms characterized by (l,m,j,p,q) = ⁽2,2,2,0,0⁾⁴(74)The tidal potential becomes (75)If the star can be considered a point-mass perturber (R ≪ a), then rapidly decays with k. So, the terms of orders k> 2 are only generated by the thermal forcing. We also assume that the horizontal pattern of J is well represented by the function , which allows us to ignore other terms. Finally, we consider the case |ν| ~ 1 (|n_orb| ≪ |Ω|, see Eq. (22)), where and for n> 0 (). In this simplified framework, the second-order Love number can be approximated in magnitude order by (76)To conclude this section, we compute the tidal torque exerted on the atmosphere with respect to the spin axis. The monochromatic torque associated with the frequency σ writes (Zahn 1966) (77)where the notation V designates the volume of the atmospheric shell. Denoting the phase lag between the forcing and the response, we obtain (78)and the total tidal torque, (79)As previously done for the Love numbers, is expanded using Hough modes (80)with l ≥ |m|, n ∈Z, and k ≥ |m|. Indeed, we note here that is the component of density variations represented by caused by the component of the excitation represented by and projected on Θ_n. Hence, the parameter is the phase difference between the component of degree of the response and the component of degree l of the excitation. When the tidal gravitational potential is quadrupolar (R ≪ a), the terms of orders higher than l = 2 can be neglected. Denoting and the quantities and for l = 2 (which should not be confused with the projections and of the forcings on the set of Hough functions), it follows (81)If the forcing is a positive real function, then the previous expression becomes (82)Substituting the polarization relation of density in Eq. (82), we obtain the torque as a sum of two contributions (83)where (84)(85)with and (86)Assuming |ε_s;n| ≪ 1 and noticing that , we establish that does not depend on the eigenvalues of Laplace’s tidal equation and can be written (87)The torque induced by this component is expanded in series of Hough functions and characterized by the solutions of the vertical structure equation. The typical scales of the vertical component are given by the vertical wavenumbers (see Eq. (47)) and can be very small compared to the scale of the horizontal ones. For example, in the vicinity of synchronization, if the radiative damping is ignored (CL70), which means that the characteristic scales of variations rapidly decay when σ → 0. Finally, we note that the series of modes in Eq. (81) can often be reduced to one dominating term. For instance, in the case of the Earth’s semidiurnal tide, and the expression above (Eq. (81)) can be simplified (88)

3. Tides in a thin stably stratified isothermal atmosphere

We establish in this section the analytical equations that describe the tidal perturbation in the case where H_atm ≪ R. This theoretical approach has been exhaustively formalized for the Earth in “atmospheric tides” (CL70) and corresponds to the case of thin atmospheres. The formalism developed in the previous section allows us to go beyond this pioneer work, and more particularly to study regimes near synchronization thanks to the inclusion of radiative losses.

3.1. Equilibrium state

For a thin atmosphere (r ≈ R), the equilibrium structure can be determined analytically. We introduce the altitude z, such that r = R + z, which is here a more appropriate coordinate than r.

The fluid is supposed to be at hydrostatic equilibrium, which leads to (89)The local acceleration in R_E;T, denoted a_RE;T, is not equal to g because of the rotating motion. The total acceleration can be written (90)where a_c is the centrifugal acceleration: (91)Here, Coriolis acceleration does not intervene because the global circulation of the atmosphere is ignored. The centrifugal acceleration (Eq. (90)) can be neglected with respect to the gravity if the rotation rate of the planet satisfies the condition Ω ≪ Ω_c, where is the critical Keplerian angular velocity. In the case of an Earth-like planet, with g ~ 10 m s^-2 and R ~ 6.0 × 10³ km, Ω_c ~ 20 rotations per day. From Eqs. (2) and (89), a typical depth can be identified, (92)which is the local characteristic pressure height scale. The parameter H represents the vertical scale of variations of basic fields. It is of the same order of magnitude as H_atm and allows us to introduce the reduced altitude (93)which is used in the following. Assuming that p₀, ρ₀, and T₀ only depends on x, we deduce their expressions from Eqs. (2) and (89) (94)In the fluid shell, the gravitational acceleration does not vary very much (95)which means that g can be considered as a constant. Assuming that H also does not vary with the altitude, we obtain a constant profile of the temperature and exponentially decaying profiles of the pressure and density: (96)In this case, which corresponds to an isothermal atmosphere, 95% of the mass of the gas is contained within the interval x ∈ [0,3]. Therefore, we can write H_atm ≈ 3H. For the Earth, H ~ 8 km (CL70) and H_atm ~ 24 km.

3.2. Wave equation and polarization relations

The constant H hypothesis is very useful in order to simplify the expressions of Sect. 2. Indeed, with this approximation, the typical frequencies of the system no longer vary with the radial coordinate. The Brunt-Väisälä frequency simply writes (97)and the acoustic frequency associated with the mode , (98)with . Moreover, the horizontal structure of tidal waves does not change because Laplace’s tidal equation is not modified. So, the vertical structure equation and the radial polarization relations only are affected by the constant H hypothesis. The coefficients defined by Eq. (41) become (99)where (100)The parameter h_n, which writes (101)is usually called the “equivalent depth” in the literature (e.g. Taylor 1936; Chapman & Lindzen 1970) and represents a typical length scale associated with the mode of degree n. We note that the curvature term in Eq. (42) vanishes in the shallow atmosphere approximation because the coefficient B₁ (Eq. (37)) is a constant in this case. This allows us to simplify the vertical structure equation (Eq. (46)). Unlike the vertical wavenumber in the thick shell (Eq. (47)), the vertical wavenumber does not vary with x in this case. It is expressed (102)with the scale ratio . It can also be written (103)The term ς_n(ε_s;n−1) comes from the radial acceleration in the Navier-Stokes equation (Eq. (11)) and is negligible for usual values of (|ε_s;n| ≪ 1 and ς_n ∝ H²/R² ≪ 1). Therefore, by considering the case where σ₀ ≪ |σ|, we recover the wavenumber obtained by CL70 for the Earth’s semidiurnal tide. In the general case, the expression of Eq. (103) can be approximated by (104)which gives for the real and imaginary parts, introducing the characteristic depth h_c = 4κH and fixing , (105)with (106)

The behaviour of the vertical wavenumber defines the possible regimes of the perturbation. These regimes are represented by the map in Fig. 5, which shows the real and imaginary part of the vertical wavenumber as colour functions of the tidal frequency and critical frequency, σ_c;n, given by (the frequency σ_c;n is proportional to the Lamb frequency defined by Eq. (98)). Propagative modes are characterized by (right panel, dark regions), evanescent modes by (yellow regions). The transition between the radiative regime and the dynamic regime corresponds to |σ| ≈ σ₀.

In the case of the thin isothermal atmosphere, the equation giving the vertical profiles (Eq. (46)) writes (107)where C is given by Eq. (99). Moreover, if the condition (108)is satisfied, the forcing in the right-hand member is dominated by thermal tides and the contribution of the gravitational tidal potential can be ignored in Eq. (107), which is rewritten (109)The polarization relations of the thin atmospheric shell are deduced from Eqs. (51) to (59)

with (119)and where , ℬ_n, and are the constant coefficients of Eq. (60)⁵, (120)

3.3. Boundary conditions

Solving Eq. (107) requires that we choose two boundary conditions. Following CL70, we fix ξ_r = 0 at the ground. This corresponds to a smooth rigid wall and is equivalent to at x = 0. In the case of Earth-like exoplanets where the telluric surface is well defined, this condition is relevant. Thus, given the form of Eq. (107), can be written (121)where A and B are the complex functions (122)The parameter K ∈C is an integration constant which is fixed by the upper boundary condition.

The upper border should be treated very carefully, as has been discussed by Green (1965). The fluid envelopes of stars and giant gaseous planets, which are considered as bounded fluid shells, are usually treated with a free-surface condition applied at the upper limit x = x_atm (e.g. Unno et al. 1989), i.e. (123)Hence, (124)with the complex coefficients (125)The atmosphere then behaves as a wave-guide of typical thickness x_atm (Appendix E). In this case, atmospheric tides are analogous to ocean tides (see Webb 1980), the amplitude of the perturbation being highly frequency-resonant. However, this condition appears to be inappropriate for the thin isothermal atmosphere because the fluid is not homogeneous and there is no discontinuous interface in this case. Setting the free-surface condition would give birth to unrealistic resonances depending on x_atm (see Fig. E.1 in Appendix F). Therefore, we have to choose a condition consistent with the exponentially decaying density and pressure. Following Shen & Zhang (1990), we consider that there is no material escape when x → + ∞, i.e. V_r|_{x → ∞} = 0. This means that the amplitude of oscillations decrease with the altitude and, consequently, that we have to eliminate the diverging term in Eq. (121)⁶.

Let us assume . As shown by CL70 (Eq. (88)), the above condition requires that (126)at the upper boundary (x = x_atm). To make this condition relevant, x_atm must be chosen so that x_atm/x_crit ≫ 1, where x_crit is the typical damping depth computed from the case where C is a constant, and given by (127)we usually fix x_atm = 10 in computations. The integration constant of Eq. (122) is then obtained as (128)

3.4. Analytical expressions of the semidiurnal tidal torque, lag angle, and amplitude of the pressure bulge

If and are assumed to be constant (with ), C does not vary with x and the expression of Eq. (121) can be simplified significantly. Indeed, in this case it can be written (129)We note that we use the convention , which implies at infinity. Substituting this simplified solution in polarization relations (Eqs. (110)to (118)), we compute three-dimensional variations of the velocity field, displacement, pressure, density, and temperature caused by the semidiurnal tide as explicit functions of the physical parameters of the atmosphere and tidal frequency. Let us introduce the function of the altitude, parametrized by (which can be , ℬ_n, or ), and defined by (130)Then, polarization relations write

with (140)From the analytical solution and the polarization relation of density given by Eq. (138), we then compute the second-order Love number and the tidal torque. Integrating δρ over the whole thickness of the atmosphere, we obtain (141)where the expressions of ℬ_n and are given by Eqs. (120) and (100), respectively. The semidiurnal tide is supposed to be represented by spatial forcings of the form and with the tidal frequency σ = 2(Ω−n_orb). The second-order Love number, deduced from Eq. (75), is thus expressed (142)Assuming |ε_s;n| ≪ 1, we identify the different contributions. It follows (143)where the subscripts _therm and _grav indicate the origin of a contribution (thermally or gravitationally forced) and _H and _V its horizontal and vertical components. The terms of this expansion are expressed as functions of the parameters of the atmosphere and tidal frequency

the parameter being a dimensionless constant given by (148)Contrary to gravitational Love numbers ( and ) which are intrinsic to the planet, thermal Love numbers ( and ) are proportional to the forcings ratio J₂/U₂ and thus depend on the properties of the whole star-planet system, particularly the semi-major axis and the stellar luminosity (e.g. Correia et al. 2008).

In the same way as we obtained Love numbers, the torque exerted on the atmosphere may be computed by substituting Eq. (141) in Eq. (80). This torque writes (149)Finally, assuming |ε_s;n| ≪ 1 and introducing the factor (150)allows us to write under the form (151)with

Let us consider the case of the pure thermal tide (U₂ ≈ 0) at the vicinity of synchronization. If σ → 0, then , |ℬ_n| → + ∞, and . As a consequence, and . The tidal torque exerted on the stably stratified isothermal atmosphere is thus very small compared to the horizontal thermal component,(156)where we have assumed that the thermal forcing is in phase with the perturber (J₂ ∈R). This behaviour had been identified before, for instance by Arras & Socrates (2010) who studied thermal tides in hot Jupiters. According to Arras & Socrates (2010), at the vicinity of synchronization, the vertical displacement of fluid due to the restoring effect of the Archimedean force (N²V_r term in the heat transport equation, Eq. (14)) exactly compensates the local density variations generated by the thermal forcing, which annihilates the quadrupolar tidal torque. We discuss this effect in Sect. 6 when considering the isothermal and stably stratified atmospheres of the Earth and a Venus-like planet (see Fig. 14). For a stably stratified structure, the tidal torque exerted on the atmosphere is weak and cannot balance the solid torque, which leads the planet’s spin to the synchronization configuration. We note here that the importance of stable stratification has been pointed out for the case of Jupiter-like planets by Ioannou & Lindzen (1993b) (see also Ioannou & Lindzen 1993a, 1994), who showed that the atmospheric tidal response was very sensitive to variations in the Brunt-Väisälä frequency (N).

3.5. Comparison with CL70

The vertical wavenumber is a key parameter of the tidal perturbation. Hence, to highlight the interest of taking into account dissipative mechanisms such as radiation, we consider the relative difference between the given by Eq. (103) and the value established via the classical theory of tides without dissipation (e.g. Wilkes 1949, CL70), denoted in reference to CL70 and given by (157)

This difference, written (158)is plotted in Fig. 6 as a function of the reduced tidal frequency σ/σ₀, and of the height ratio h_c/h_n. The supplementary term of Newtonian cooling in the heat transport equation induces an additional regime with respect to CL70. At low frequencies, around the radiation frequency, the vertical profiles are damped, while they are highly oscillating in the case where σ₀ = 0. This thermal regime corresponds to the middle white-yellow region of the map. The two wavenumbers become similar for |σ| ≫ σ₀. This is the regime of fast rotating planets studied by CL70. Beyond the threshold corresponding to |σ| ~ N, the traditional approximation assumed in Sect. 2 is no longer valid and the coupling between the three components of the Navier-Stokes equation is considered. Finally, the area located at h_c/h_n ≈ 1 represents the discontinuous transition between propagative and evanescent waves in CL70, this transition being regular in our model.

4. Tides in a slowly rotating convective atmosphere

In this section, we treat the asymptotic case of a slowly rotating convective atmosphere for which the equations of tidal hydrodynamics can be simplified drastically. This situation typically corresponds to Venus-like planets. Indeed, the rotation period of Venus, 243 days, is of the same order of magnitude as its orbital period. It follows that |ν| < 1 for the semidiurnal tide of Venus. As a consequence, the Coriolis acceleration can be neglected as the first step. The atmospheric layer located below 60 km is characterized by a strongly negative temperature gradient (Seiff et al. 1980). In this layer, the temperature decreases from 750 K to 250 K, which make it subject to convective instability (Baker et al. 2000). Therefore, the Brunt-Väisälä frequency of Venus is far lower than the value given by the isothermal approximation. As this frequency represents the “stiffness” of the Archimedean force which restores gravity waves, these waves cannot propagate in the unstable region above the surface.

4.1. Equilibrium distributions of pressure, density, and temperature

We assume Ω = 0 and N = 0. According to Eq. (1), the Brunt-Väisälä frequency can be expressed as(159)Therefore, the case N² = 0 corresponds to the adiabatic temperature gradient, given by(160)For this structure, the reduced altitude is expressed by (161)and we thus obtain the distributions of the basic pressure height, temperature, pressure, and density, respectively (162)with .

4.2. Tidal response

For the sake of simplicity, we neglect the contribution of the gravitational forcing in the Navier-Stokes equation and examine the case of pure thermal tides. When σ → 0, the tidal response tends to its hydrostatic component, the equilibrium tide (Zahn 1966). We thus assume the hydrostatic approximation. This allows us to reduce the radial projection of the Navier-Stokes equation (Eq. (8)) to(163)Now, by expanding the perturbed quantities in Fourier series and substituting Eq. (163) in the heat transport equation (Eq. (20)), we get the simplified vertical structure equation of pressure (164)represents the complex damping rate of the perturbation. Following Dobrovolskis & Ingersoll (1980) and Shen & Zhang (1990), we consider that thermal tides are generated by an increase in heat at the ground and apply the profile of thermal forcing (165)where 1 /τ_J ≪ 1 is the thickness of the heated region (see Eq. (213)for a physical expression of this parameter) and J₀ the mean absorbed power per unit mass. At the upper boundary x_atm ≫ 1, we apply the free-surface condition δp = 0 and get the pressure profile(166)from which we deduce the vertical profiles of the Lagrangian horizontal displacements (167)(168)of the horizontal velocities(169)(170)and of the density and temperature (171)(172)

4.3. Second-order Love number and tidal torque

The previous results enable us to compute the tidal Love numbers and torque associated with the atmospheric tidal response. Owing to the hydrostatic approximation, these quantities are directly proportional to δp⁽x = 0⁾. We consider that the semidiurnal thermal tide can be reduced to its quadrupolar component, given by and . Hence, using the expressions Eq. (76) for the second-order Love number and Eq. (82) for the tidal torque, we obtain (173)and(174)

4.4. Comparison with previous models

The tidal torque plotted in Fig. 7 is identical to the horizontal component of the torque exerted on the stably stratified isothermal atmosphere given by Eq. (156) because convection eliminates the component resulting from the fluid vertical displacement. In agreement with the early result of Ingersoll & Dobrovolskis (1978) (Eq. (4)), it is of the same form as the result given by the Maxwell model (Correia et al. 2014), where we identify the Maxwell time . It also corresponds to the torque computed by Leconte et al. (2015) for Venus-like planets with numerical simulations using GCM, which suggests that the slow rotation and convective instability approximations are appropriate for planets of this kind. This torque is stronger than that computed in the case of the stably stratified atmosphere. Consequently, it can lead to non-synchronized rotation states of equilibrium by counterbalancing the solid tidal torque. Finally, it must be compared to the model introduced by Correia & Laskar (2001) in early theoretical studies of the tidal torque exerted on Venus atmosphere. In this model, the tidal torque is given by (175)where a and b are two empirical positive real parameters. As can be observed in Fig 7 where (174) and (175) are plotted as functions of the tidal frequency, looks like from Eq. (174) with amplitude maxima located at (176)This strengthens the statement of Correia & Laskar (2001), who expected that “further studies of atmospheres of extra solar synchronous planets [ ...] may provide an accurate solution for the case σ ≈ 0”. In Fig. 14, the expression of Eq. (175) is plotted for the Earth and Venus in addition to the analytical results of the model presented here with parameters a and b adjusted numerically.

5. Physical description of heat source terms

On the left-hand side of Eq. (46) and in polarization relations, the perturbation is driven by and . The tidal gravitational potential has been expanded in spherical harmonics for many years and there is no need to come back to it⁷. Nevertheless, it is necessary to establish physical expressions of the heat power per unit mass to compute a solution. The different contributions must be clearly identified. We detail them in the case of the optically thin atmosphere. The most obvious heat source, but not necessarily the one with the largest amplitude, is the absorption of the light flux coming from the star. The non-absorbed part of this flux reaches the ground and causes temperature oscillations at x = 0. This implies two contributions from the ground: a radiative emission in infrared and diffusive heating due to turbulence within the surface boundary layer.

5.1. Insolation

The heating fluxes coming from the ground are all oriented radially and proportional to the local insolation flux at x = 0. Therefore, they can be easily written as a product of separable functions: . This is not the case of the insolation flux, which goes through the spherical shell along the star-planet direction (see Fig. 10). Particularly, the atmosphere is partly enlightened by stellar rays on the dark side. We note I_S, the bolometric flux transmitted to the atmosphere. If we assume that the star radiates like a black body, this flux is given by the integral of the spectral radiance I_S;λ over a wide range of wavelengths (177)where I_S;λ is given by Plancks’s law (Planck 1901): (178)the parameter c being the speed of light in vacuum, h the Planck constant, k_B the Boltzmann constant, d the star-planet distance, R_S the radius of the star, and T_S its surface temperature. The flux is partly reflected by the atmosphere with the albedo A_atm;λ for the wavelength λ, which gives the effective heating flux (179)To describe the absorption of I₀ by the atmosphere, we use the Beer-Lambert law (Bouguer 1729; Klett et al. 1760; Beer 1852). Light rays propagate along a straight line defined by the linear spatial coordinate l. So, the absorption of the incident spectral density by a gas G of molar concentration C_G, illustrated by Fig. 8, can be written (180)the parameter ε_G;λ being the molar extinction density coefficient of the gas, which is supposed to depend on the light wavelength (λ) alone. This dependence should be taken into account because ε_G;λ can vary over several orders of magnitude, typically, ε_G;λ ∝ λ^-4 in the regime of Rayleigh diffusion (molecules are small compared to the wavelength). From Eq. (180), we get (181)and the incident flux (182)

Then, performing the derivation of Eq. (181), the heat power per unit volume is obtained, (183)and the heat power per unit mass comes straightforwardly (184)We retrieve the formulation of the thermal forcing proposed by CL70. Considering that the gas is homogeneous in composition leads us to C_G = C₀ = ρ₀/M. Therefore, the incident power flux and heat power per unit mass become (185)Near the planet-star axis, they are simply expressed as (186)where τ_λ is the dimensionless optical depth for the radiance of wavelength λ, (187)The parameter 1 /τ_λ corresponds to an absorption depth of the flux in a homogeneous fluid normalized by the pressure height scale (H). If τ_λ ≪ 1, the heating power of a wavelength λ is almost totally transmitted to the ground. On the contrary, if τ_λ ≫ 1, the flux does not reach the surface because it is entirely absorbed by the atmosphere (see Fig. 9). In this case, the atmosphere is not optically thin, and it cannot be assumed that its radiative losses are proportional to δT as was done in Sect. 2, and thermal diffusion has to be considered. Therefore, the case treated by the present work corresponds to τ_λ ≪ 1.

The proportion of flux density transmitted to the ground α_λ is deduced from Eq. (186) taken at x = 0, (188)Like the atmosphere, the ground has an albedo denoted A_gr;λ, which implies that the upcoming flux is partly reflected. Therefore, we can write the absorption equation (189)and compute the reflected surface power (190)the coordinate ψ being the angle between the star-planet direction and the position vector of the current point. The heat power per unit mass and wavelength provided by the reflected flux is derived from the previous equation (191)So, the contribution of the reflected flux is finally obtained (192)

5.2. Radiative heating from the ground

The part of the incident flux which is not reflected causes surface temperature oscillations. This effect was studied in meteorological works throughout the twentieth century. We will follow here the theoretical approach proposed by Bernard (1962) for its simplicity and the physical landmarks that it brings. The ground is considered to be a black body of temperature T_gr, emitting the power flux I_gr;BB given by the Stefan-Boltzmann law, (193)the parameter ϵ_gr being the emissivity of the ground and S= 5.670373 × 10^-8 W m^-2 K^-4 the Stefan-Boltzmann constant (Mohr et al. 2012). Moreover, the atmosphere emits a counter radiation to the surface, which implies that the effective flux is less than I_gr;BB. According to Bernard (1962), this counter-radiation can be assumed to be proportional to I_gr;BB, and expressed by the semi-empirical formula (194)which allows us to take it into account without studying the whole coupled system ground-atmosphere. The factor ϵ_atm corresponds to an effective emissivity. So, introducing the spectral radiance of the atmosphere I_atm;λ, the albedo of the ground in the infrared A_gr;IR can be defined by the relationship (195)and the effective flux emitted by the telluric surface can be written (196)Given that , Eq. (196) becomes (197)where ϵ = ϵ_gr⁽1−ϵ_atm⁾ ≤ 1 is the effective emissivity of the ground. The terrestrial atmosphere being rather thin and transparent, ϵ ≈ 1 for the Earth. This expression will be used further to establish the heating by turbulent diffusion. We also introduce here the corresponding spectral radiance of the ground, (198)The radiative flux emitted by the ground, denoted I^rad, is the solution of the absorption equation (similar to the one used to compute I^inc), (199)with the boundary condition at x = 0. Therefore, (200)and finally, (201)The temperature of the ground T_gr depends on the power reaching the surface, (202)Both quantities are linearized near the equilibrium, and T_gr = T_gr;0 + δT_gr, and the perturbation is expanded in series (203)The spatial functions and are related by a transfer function of the tidal frequency, denoted B_gr, that is given explicitly in the next subsection (204)So, the perturbation of the spectral radiance emitted by the ground can be written (205)and the flux and heat power per mass unit (206)The partial derivative of is explicitly given by (207)Here, we note that the response of the ground to the tidal forcing induces a dependence on the tidal frequency through the oscillations of T_gr and the transfer function . Thus, the spatial functions describing the contribution of the ground are parametrized by σ contrary to those associated with the incident flux.

5.3. Boundary layer turbulent heat diffusion

Near the ground, heat transfers are dominated by turbulent mechanisms and diffusion. This is known as the planetary boundary layer and it is illustrated in Fig. 12. Thus, for x ≪ 1, thermal diffusion cannot be ignored as it is in Sect. 2 and drives the oscillations of temperature. However, we note that it is only valid at the vicinity of the surface, where friction with the ground is important. After CL70 and Siebert (1961), we neglect as a first step the heat transport through the troposphere due to advection to provide a first analytical treatment of the problem. At the interface, the diffusive fluxes in the ground and in the atmosphere write (208)the parameters k_gr and k_atm representing the thermal conductivities of the ground and of the lowest layer of the atmosphere. In addition, denoting the associated thermal capacities c_gr and c_atm, we introduce the thermal diffusivities (209)the parameter K_atm being the vertical turbulent thermal diffusivity. Therefore, temperature variations near the surface are described by the equations of heat transport (e.g. CL70), (210)where only the diffusive terms are taken into account. The power balance of the surface is deduced from Eqs. (208) and (197), (211)and provides the temperature T_gr;0 of the ground at equilibrium. Linearized with variables of the form given by Eq. (21), Eq. (210) can be written (212)where and are the tidal frequency-dependent skin thicknesses of diffusive heat transport in the solid and fluid parts, respectively. Their expressions are given by (213)It should be noted that and both decay when the tidal frequency increases. Since corresponds to the typical depth of the boundary layer, it has to satisfy , otherwise the diffusive effects cannot be ignored in the dynamical core of the model (Sect. 2). Let us denote the radiative impedance of the ground considered as a perfect black body and submitted to a perturbation in temperature. The transfer function introduced in Eq. (204) can be deduced from the linearized power balance, (214)in which the expressions of Eq. (212) are substituted. We finally get (see Fig. 11) (215)

an expression that can be decomposed in its modulus and argument (216)The frequency σ_bl, which parametrizes , is a characteristic frequency reflecting the thermal properties of the diffusive boundary layer, and is expressed (217)where the parameters β_gr and β_atm are the thermal conductive capacities of the ground and of the atmosphere near the telluric surface (218)According to Bernard (1962), turbulent diffusion plays a more important role than thermal conduction in the ground, which means β_atm ≫ β_gr. However, the turbulent diffusivity K_atm may vary over a wide range of orders of magnitude, taking extremal values above oceans (typically β_atm ~ 10⁴ J m^-2 K^-1 s^{− 1/2}). For these reasons, Bernard (1962) prescribes for the Earth the effective mean value K_atm ~ 10 m² s^-1, which had been prescribed before by Wilkes (1949) and which leads to σ_bl ~ 10^-6 s^-1.

As shown by Eq. (216), the frequency ratio |σ| /σ_bl determines the angular delay of the ground temperature variations. If |σ| ≪ σ_bl, the response is in phase with the excitation. It corresponds to low conductive capacities. On the contrary, if |σ| ≫ σ_bl, then the delay tends to its asymptotical limit, π/ 4. By using the frequency σ_bl ~ 10^-6 s^-1 derived from the prescriptions of Wilkes (1949) and Bernard (1962), we note that the Earth’s solar diurnal tide belongs to the second category. This explains the position of the diurnal peak observed in surface temperature oscillations, three hours late with respect to the subsolar (midday) point. The gain of is damped by the conductive capacities of the material at the interface. Therefore, a very diffusive interface tends to attenuate the temperature oscillations of the ground, as expected.

We now consider the profile of temperature variations in the boundary layer obtained in Eq. (212) (219)where is the vertical damping ratio of the diffusive perturbation. Chapman & Lindzen (1970) shows that a heat power due to diffusion J^diff can be deduced from the previous expression. Indeed, by writing the heat transport equation as (220)and linearizing it, we compute the expression of the tidal diffusive forcing term as (221)

5.4. Simplified bichromatic case

The expressions of the thermal forcings obtained in the previous section may appear of some complexity because of the dependence of albedos and optical depths on the wavelength. They can be simplified notably if each source is assumed to be monochromatic. In this case, the star emits in the frequency range of visible light, and the planet in the range of infrared. Consequently, the wavelength-dependent A_gr;λ, A_atm;λ, ε_λ, ϵ_λ, α_λ, and τ_λ are reduced to the constant coefficients A_gr;V, A_atm;V, ε_V, ϵ_V, α_V, and τ_V in the range of visible light and A_gr;IR, A_atm;IR, ε_IR, ϵ_IR, α_IR, and τ_IR in the infrared. Hence, the bolometric flux is given by (222)where T_S is the temperature of the star. This flux is partly reflected back to space by the atmosphere so that the effective flux I₀ getting through the layer can be written (223)The power of the incident flux I^inc decays along the path of a light beam, absorbed by the gas. The power absorbed per unit mass is denoted J^inc and the expressions of Eq. (185) become (224)From I^inc, we deduce the flux reflected by the ground I^ref and the power absorbed per mass unit J^ref, given by (225)Finally, the expressions of Eq. (206) are simply replaced by (226)The expressions of the powers per unit mass of the incident and radiated flux (J^inc and J^rad), and the corresponding powers per volume unit (J^inc and J^rad), are plotted in Fig. 10. The heating due to the turbulent boundary layer remains unchanged (see Eq. (221)).

6. Application to the Earth and Venus

To illustrate the regimes identified and the dependence of the atmospheric response on the forcing frequency, we apply the model to the cases of the semidiurnal thermal tides on Earth and Venus assuming first for both planets isothermal stably stratified atmospheres, as described by the equations in Sect. 3. We note here that the atmosphere of Venus is not stably stratified but convective in regions close to the surface (Seiff et al. 1980; Baker et al. 2000) where the tidal mass redistribution is the most important (Dobrovolskis & Ingersoll 1980; Shen & Zhang 1990). Nevertheless, examining the case of a stably stratified isothermal atmosphere for a Venus-like planet will allow us to better understand in an academic framework the role of stratification and radiative losses in the tidal response. To make a clear difference between Venus and the studied stably stratified Venus-like planet, we call the latter “VenusX”. For the sake of simplicity, we assume that planets orbit circularly in their equatorial plane. Hence, the tidal frequency⁸ is given by σ = 2(Ω−n_orb). One of the most crucial parameters of the model is the Newtonian cooling frequency (σ₀). This parameter varies over several orders of magnitude with the altitude (Pollack & Young 1975). However, our purpose here is not to compute a quantitative tidal perturbation, but to illustrate qualitatively the non-dissipative (Earth) and dissipative (Venus) regimes. Therefore, following Lindzen & McKenzie (1967), we consider the interesting case where the radial profile σ₀ is assumed to be constant. For both planets, we arbitrarily set σ₀ = 7.5 × 10^-7 s^-1, which is the effective value of σ₀ computed by Leconte et al. (2015) with a GCM for Venus. This value is such that σ ~ σ₀ for Venus and σ ≫ σ₀ for the Earth. In the case of an optically thin atmosphere, the Newtonian cooling frequency can be estimated using Eq. (17). Moreover, we impose academic quadrupolar perturbations of the form (227)where U₂ and J₂ are fixed constant radial profiles. The tidal potential U₂ is computed from the Kaula multipole expansion detailed in Appendix D, and explicitly given by the expression of Eq. (74). A zero-order approximation of the thermal power J₂ can be quantified by writing the total surface power absorbed by the atmosphere as a function of θ and ϕ, which is then expanded in Fourier series of longitude and associated Legendre polynomials. Given that atmospheric tides are considered in this model as linear perturbations around the equilibrium state, the amplitudes of the perturbed quantities are proportional to the forcing.

The horizontal structure equation is solved using the spectral method described in Appendix B. The vertical structure equation is integrated numerically on the domain x ∈ [0,x_atm] using a regular mesh with element of size δx. As pointed out by CL70, the number of points in the mesh, denoted N, must be sufficiently large to obtain vertical profiles of the perturbation with a good accuracy (see Appendix G). To fix N correctly in the case of the Earth, CL70 suggests using a criterion that we adapt for other configurations, i.e. (228)considering that the scale of the vertical variations is defined by the module of the complex wavenumber given in Eq. (107).

For each planet, we give here the spatial distribution of perturbed quantities and the evolution of the tidal torque with the tidal frequency. The numerical values used in these simulations are summarized in Table 1.

6.1. The Earth

The Earth’s semidiurnal thermal tide corresponds to one of the cases detailed in CL70. Here, the radiation frequency is clearly negligible compared to the tidal frequency σ ≈ 2Ω. Therefore, the tidal response of the Earth’s atmosphere is driven by dynamical effects only and the radiative losses added in the heat transport equation can be ignored. This case is typical of a pure dynamic behaviour (see Figs. 5 and 6). Given that ν is slightly larger than 1, the horizontal structure of the perturbations is essentially described by gravity modes (Fig. 2). Rossby modes are trapped at the poles. Figures 15 and 16 show maps of the perturbed quantities in latitude (denoted δ), longitude (ϕ), and altitude (x), the subsolar point being indicated by the coordinates (δ,ϕ) = ⁽0,0⁾. In these figures, we can observe that the pressure and density variations present very similar spatial distributions. Particularly, they are in phase with each other and the phase lag is approximately the same in the whole layer, which causes the tidal bulge to materialize. The lag angle γ = π/ 4 of the bulge with respect to the subsolar point is given by Eq. (137) taken at x = 0 for the gravity mode of lowest degree (n = 0) and is not very sensitive to temperature structure (Lindzen 1968). For the Earth, these semidiurnal pressure peaks are measured around 9h38mn and 22h07mn (see Fig. 13).

The patterns of the velocities and are also explained well by the first horizontal Hough functions, and (see Fig. 2, middle and bottom left). In contrast, the behaviour of the temperature involves other modes.

We then compute the evolution of the tidal torque with the tidal frequency by changing the value of the spin frequency (Ω) in simulations. Hence, the reduced tidal frequency χ = (Ω−n_orb) /n_orb varies within the interval [− 15,15]. The torque is computed using the analytical formulae of Eqs. (141) and (149). We also compute the horizontal and vertical components isolated in Sect. 3. The results are plotted in Fig. 14.

At first, we note that the horizontal and vertical components of the torque are of the same order of magnitude and in opposition to one another, as shown by Eq. (141) taken at the limit σ → 0. This behaviour was identified previously in Sect. 3.4. It results from the fact that vertical displacement of fluid induced by the stable stratification compensates the local density variations in this frequency regime and strongly flattens the total tidal torque. As demonstrated in Sect. 4, this is not the case in general and it depends on the strength of the stratification. In the case of a slowly rotating convective atmosphere (N² ≈ 0) forced by a heating located at the ground, only the horizontal component remains. This leads to a much stronger tidal torque, similar to those obtained by Leconte et al. (2015) with a GCM (given by Eq. (174)) and Correia & Laskar (2001) (see Eq. (175)).

In Fig. 17, we plot the surface variations of pressure, density, temperature, and their components as functions of the tidal frequency. For each quantity, the horizontal component varies smoothly while the vertical one is discontinuous at the synchronization. This discontinuity is a consequence of the dissymmetry between prograde and retrograde Hough modes due to the Coriolis acceleration (Fig. A.1), particularly as regards the gravity mode of the lowest degree, which is the most important one.

Assuming that the tidal torque is entirely transmitted to the telluric core of the planet, it is possible to estimate the timescale of the spin evolution induced by the atmospheric semidiurnal tide. Introducing the moment of inertia of the planet, ℐ, we write (229)Thus, the timescale of the spin evolution is given by (230)For the Earth’s semidiurnal tide, χ ≈ 365 and the horizontal component of the corresponding torque is . With ℐ = 8.02 × 10³⁷ kg m² (NASA fact sheets¹⁰), T_⊕ ~ 180 Gyr. Therefore, the spin frequency of the Earth is not affected by atmospheric tides.

6.2. Venus

Geometrically speaking, Venus can be seen as an Earth-like planet, with scales of the same magnitude order for the telluric core and the thickness of the atmosphere. However, there are some important differences compared to the case of the Earth. The Venusian atmosphere does not have the same properties as the Earth’s atmosphere. First, the fluid layer is a hundred times more massive (denoting M the mass of the atmosphere, M_♀ ~4.8 × 10²⁰ kg and M_⊕ ~ 5.1 × 10¹⁸ kg). Thus, the surface pressure is around 93 bar. Second, it is optically thicker. The major part of the heating flux is thus deposited at high altitudes. The absorbed flux near the planet surface at the subsolar point F_abs is estimated at F_abs ~ 100−200 W m^-2 (Avduevsky et al. 1970; Lacis 1975; Dobrovolskis & Ingersoll 1980), which implies J₂ ~ 10^-4 W kg^-1. Finally, layers below 60 km are characterized by a strongly negative temperature gradient (Seiff et al. 1980) and are, therefore, weakly stratified or convective. As discussed in Sects. 3 and 4, these properties have a strong impact on the nature of the tidal response. Particularly, the tidal torque may be much stronger in the convective case than in the stably stratified one. In the present section, we study a Venus-like planet with a stably stratified isothermal atmosphere (VenusX) in order to understand the effects of stratification and radiative losses in a general context.

The other important difference between the Earth and Venus is in their rotational dynamics. While the spin frequency of the Earth is far higherthan its mean motion and the radiation frequency of its atmosphere, these three frequencies are of the same order of magnitude for Venus. Therefore, the Venusian semidiurnal thermal tide is characterized by a tidal frequency |σ| ~ σ₀, which is typical of the thermal regime. The term describing radiative losses in the heat transport equation (Eq. (18)) plays an important role in this regime by damping the profiles in altitude of the tidal response; see Fig. 5. This figure shows that the imaginary part of , which is responsible for the damping, is always comparable to the real part in this frequency range.

Contrary to the case of the Earth, the semidiurnal thermal tide of VenusX belongs to the family of super-inertial waves, with ν ≈ 0.48, which means that the horizontal component of the tidal response is composed of gravity modes only. We note that the surface variations of perturbed quantities plotted in Fig. 19 are very well represented by the first mode. The patterns of the Hough functions , , and (see Fig. 3, red curves) clearly appear on the maps. Moreover, the lags of δp and δρ are different from that observed in the Earth’s semidiurnal tide. In particular, it is interesting to note that the pressure and density peaks are not superposed in this case. The effect of the damping caused by radiative losses can be observed in vertical cross-sections in Fig. 18, where the variations of pressure and density are located close to the ground. Given that the real and imaginary part of the vertical wavenumber are both very high, vertical component is both strongly oscillating and strongly damped.

As was previously done for the Earth, we draw the variations of the tidal torque (Fig. 14) and perturbed quantities at the ground (Fig. 20) around the synchronization for χ ∈ [−15,15]. The behaviour of the torque and its components is the same as previously observed for the Earth. The maxima of horizontal and vertical components are located at χ = ± χ_♀ with χ_♀<χ_⊕ owing to the difference of mean motions between the two planets. Since the tidal torque is proportional to the surface density of the atmosphere, it is a hundred times stronger in the case of Venus than in the case of the Earth for a given thermal forcing. This is also verified by the surface variations of the perturbed quantities, which – except for their amplitudes – are comparable to those of the Earth.

To conclude this section, we compute the timescale of the spin evolution. The Venusian semidiurnal tide is characterized by the frequency χ ≈ −1.92. With the moment of inertia ℐ = 5.88 × 10³⁷ kg m² (NASA fact sheets) and the horizontal component of the torque given by Fig. 14, i.e. , we obtain T_♀ ~ 10 Myr. The thermal forcing corresponds approximately to the mean level of the heating profile used by Dobrovolskis & Ingersoll (1980) (J₂ ~ 10^-4 W kg^-1) so the obtained torque is of the same order of magnitude as the value given by this early work, . We note here that the amplitude of the response is directly proportional to the amplitude of the excitation if one of its two components, J or U, can be neglected (Sect. 2).

6.3. Discussion of physical assumptions

We have developed a new analytic model for atmospheric tides by generalizing the classical tidal theory detailed in CL70 to dissipative atmospheres.

This work leans on several assumptions which contribute to make computations simpler. The most important of them are listed here.

• The traditional approximation. This hypothesis, com-monly used in the literature of geophysics, is very usefulbecause it allows us to separate the horizontal and verti-cal dependences in the Navier-Stokes equation. Hence,the horizontal and vertical structure can be computedseparately. This approximation is relevant in stronglystably stratified fluids where 2Ω < |σ| ≪ N, but is less appropriate if|σ| ~ 2Ω, as discussed in Mathis (2009). Its validity domain, not well defined, remains essentially limited to the regime of super-inertial waves (Prat et al. 2016). It is represented in Fig. 21, which shows the cases that can be treated analytically (regions a and b). Concerning the importance of the neglected Coriolis terms, we refer the reader to Gerkema & Shrira (2005), who explored the physics of gravito-inertial waves without the traditional approximation (see also Tort & Dubos 2014). The traditional approximation appears to be a relatively strong hypothesis because it does not take into consideration the case of rotating convective atmospheres (N² ≈ 0, e.g. Ogilvie & Lin 2004). Convection is an important feature of the Earth’s nearly adiabatic layer (altitudes below 2 km) as well as the Venus’ atmosphere (Linkin et al. 1986). It has strong implications on the propagations of waves: gravity waves cannot exist in convective layers, which may significantly modify the tidal response, particularly if the convective turnover time is close to the period of the tidal forcing. In the general case, the effect of convection on the structure of waves can be studied either numerically with GCM (Fig. 21, d) or analytically with simplified models (e.g. based on a Cartesian geometry). The only exception for a global analytical study of a convective atmosphere is the case treated in Sect. 4 where stable stratification and rotation are too weak to be neglected (Fig. 21, b).

• The linearization around the equilibrium. In Sects. 2 and 3, the equations have been linearized around the equilibrium. This method is only relevant in cases where the tidal response of the atmosphere can be considered as a small perturbation with respect to the equilibrium state. This condition is verified for the Earth and Venus.

• The solid rotation approximation. We assumed that the atmosphere rotates with the solid core uniformly. Thus, the atmospheric circulation (both meridional and zonal) is ignored. This approximation is justified by the fact that the most important pressure and density variations are concentrated in the densest layers, near the ground. Given that these layers are bound to the solid core by viscous friction, they rotate at the same spin frequency. However, tidal waves might interact with the mean flow in cases where σ ≈ 0 (Goldreich & Nicholson 1989) and it will be interesting to take into account the global atmospheric circulation in further works, as is done in GCM simulations.

• The use of simplified atmospheric vertical structures. In order to lighten calculations, we considered the cases of uniform background temperature (isothermal stably stratified atmosphere) and temperature gradient (convective atmosphere), and supposed that some parameters, such as σ₀, do not vary with altitude. These simplified structures are a first step towards more sophisticated ones. Owing to the importance of the stratification in the tidal response, it will be necessary to treat the cases of observed structures in future works, for instance those of the atmospheres of Earth, Venus, and Mars.

• The optically thin layer hypothesis. Radiative processes are introduced in the dynamics of atmospheric tides through Newtonian cooling. The linear dependence of the corresponding radiative term on the temperature perturbation can be interpreted as an optically thin layer approximation. Indeed, this term corresponds to the thermal behaviour of an atmosphere that would emit grey body radiation towards space and the planet surface. Actually, the flux emitted is partly absorbed by the atmosphere itself, which induces a complex thermal coupling between the different layers of the atmosphere. This is true particularly for optically thick atmospheres. In these cases, the Newtonian cooling should be seen as a convenient approximation of thermal radiation for an analytical approach where the dissipative processes cannot be treated without approximations. This approach seems to be well validated by numerical simulations using GCM (Leconte et al. 2015).

7. Conclusions

We have set the bases of a new theoretical framework to study the atmospheric tides of Earth-like exoplanets analytically. This work is complementary to numerical simulations. The main objective is to explain the qualitative behaviour of an atmosphere submitted to both thermal and gravitational tidal perturbations and to explore the domain of key physical parameters. The parameters highlighted by the model are

• the characteristic frequencies of the system (e.g. inertia,Brunt-Väisälä, radiative, boundary layer, tidal frequencies);

• the involved characteristic length scales (e.g. planetary radius, pressure height scale of the atmosphere, length scales of the vertical variations of tidal waves);

• the physical parameters of the planet (e.g. the molar mass of the atmosphere, the adiabatic exponent, the surface gravity, the surface pressure).

Using the equations of hydrodynamics, we have adapted the theoretical approach detailed in CL70 to thick dissipative atmospheric shells. We computed the analytic expressions of the horizontal and vertical structure equations, of the perturbed quantities (pressure, density, temperature, velocity, displacement), and of Love numbers and tidal torques exerted on the atmosphere. These expressions have been used to identify the possible regimes of atmospheric tides: a non-dissipative Earth-like dynamical regime characterizing fast rotating bodies, and a dissipative Venus-like thermal regime for planets near synchronization. The behaviours involved in the two cases are very different. The dynamical regime corresponds to the case studied by CL70. In this regime, tidal waves are either oscillating or evanescent depending of the equivalent depth (h_n) associated with the mode. In the radiative regime, they are oscillating and evanescent, and the perturbation can be strongly damped.

The tidal response of the atmosphere will also be affected by the nature of the tidal excitation. A predominating thermal forcing generates thermal tides. On the contrary, if the thermal energy brought to the fluid layer is negligible compared to the effect of the tidal gravitational potential, the perturbation belongs to the family of gravitational tides. This case has not been studied in this work, but the vertical structure equation and polarization relations of the model clearly attest to behaviour that is qualitatively different from that observed for thermal tides. Finally, the tidal response of the atmospheric shell is obviously affected by lower and upper boundary conditions. We showed for instance that a free-surface condition and a no-material-escape condition would not give the same vertical profiles of the perturbation. Similarly, the ground may affect temperature oscillations by introducing a delayed thermal forcing locally.

Finally, the results obtained for the tidal torque are in agreement with the GCM simulations by Leconte et al. (2015). This torque appears to be strongly dependent both on the tidal frequency and on the atmospheric structure (Fig. 22). In the case of a stably stratified atmosphere, the tidal bulge vanishes at the vicinity of synchronization. As a consequence, the tidal torque exerted on the atmosphere is very weak and cannot balance the solid tidal torque, which should lead the planet to synchronization. On the contrary, if the stable stratification is sufficiently weak (weakly stratified or convective atmosphere), the atmospheric tidal torque can be strong. In this case we recover the Maxwellian behaviour obtained numerically by Leconte et al. (2015) and described by the early model of Correia & Laskar (2001). This would lead planets to non-synchronized states. These results highlight the predictive interest of the modelling of atmospheric tides by showing that it is possible to constrain the atmospheric structure of terrestrial planets from their rotation states of equilibrium.

---

Acknowledgments

We thank the anonymous referee for comments that helped to improve the paper. We thank Richard Lindzen for the very interesting discussions we had about the analytic modelling of the Earth’s atmospheric tides. We are also very grateful to Umin Lee for his help in understanding how to compute Hough functions efficiently. Finally, we thank Jeremy Leconte, who shared with us his helpful approach to the physical modelling of super-Earth atmospheres from the point of view of numerical simulations. P. Auclair-Desrotour and S. Mathis acknowledge funding by the European Research Council through ERC grant SPIRE 647383. This work was also supported by the Programme National de Planétologie (CNRS/INSU) and CoRoT/Kepler and PLATO CNES grant at CEA-Saclay.

References

Appendix A: Hough functions

Hough functions, named for Hough’s initial work (Hough 1898), are the solutions of Eq. (35) integrated as an eigenvalue problem of with appropriate boundary conditions (Longuet-Higgins 1968; Chapman & Lindzen 1970; Lee & Saio 1997). It is possible to compute these functions in several ways. Here, we present two possible approaches: the finite differences method proposed in the Numerical recipes (Press et al. 1986) and detailed in Lee & Saio (1997), and a spectral method developed from the works of Hough (1898), Longuet-Higgins (1968), and Chapman & Lindzen (1970).

Appendix A.1: Finite differences method

The operator ℒ^m,ν of Eq. (30) writes (A.1)where μ = cosθ. The new Laplacian tidal equation is not modified when μ is replaced by − μ. Therefore, the domain over which this equation is integrated can be reduced to 0 ≤ μ ≤ 1. Like the associated Legendre polynomials, the solutions are either even or odd with respect to the equatorial plane μ = 0. The even solutions satisfy , the odd solutions . For the sake of simplicity, the parity of the subscript used in the notations corresponds to the parity of the associated function. The solutions belonging to one of these two families are all defined by the boundary condition applied at μ = 0: for even solutions and for odd solutions. The operator ℒ^m,ν has two regular singular points, a permanent one at μ = 1 and another one at μ = 1/|ν| in the sub-inertial regime, characterized by |ν| > 1 and, therefore, depending on the spin parameter. Both are sources of numerical difficulties. To tackle the problem of the singular point at μ = 1, we introduce the variable y(μ) defined by Press et al. (1986) and Lee & Saio (1997):

(A.2)Laplace’s tidal equation becomes (A.3)where (A.4)The boundary condition applied at μ = 0 depends on the parity of the sought solution, i.e. dy/ dμ = 0 at the equator for even modes and y = 0 for odd modes. We also apply at μ = 0 the normalization condition y = 1 for even modes and dy/ dμ = 1 for odd modes. Since and are not defined at μ = 1, the boundary condition at the poles is applied at μ = 1−ϵ, where 0 <ϵ ≪ 1. This condition is derived from the Sturm-Liouville theory, which requires that should be regular at μ = 1, and is given explicitly by (see Lee & Saio 1997).

In order to solve Eq. (A.3) numerically, the domain is discretized. Hence, the equation is reduced to a linear system of dimension N_S = 2N, where N is the number of elements of the mesh. In the regime of sub-inertial waves, the singular point at μ = 1/|ν| corresponds to a turning point, in the vicinity of which the resolution of the mesh must be sufficiently high. Lee & Saio (1997) showed that a regular solution could be locally expanded in series of the form , where x = μ−1/|ν|, and satisfying the condition (A.5)This allows us to compute the solution properly in the sub-inertial regime. To construct the set of Hough functions, we use as predictors the asymptotic results given by Longuet-Higgins (1968) and Townsend (2003). Two types of modes can be identified.

The modes of the first type, called gravity modes, are defined for ν ∈R and are characterized by positive eigenvalues (Fig. A.1, top panel). In the non-rotating case, where ν = 0, they correspond to the associated Legendre polynomials (Fig. A.1, bottom panel). The gravity modes of Hough functions are ordered with positive n, by analogy with this particular case, for which the degrees l of Legendre polynomials are given by l = |m| + n. In the asymptotic cases where |ν| → + ∞, analytic solutions using Hermite polynomials can be computed (Townsend 2003). In the sub-inertial regime, the oscillations of gravity modes are trapped in an equatorial belt in the interval given by 0 ≤ μ ≤ 1/|ν| (see Fig. A.2). This explains why these modes are not well represented by associated Legendre polynomials for high values of |ν|.

The modes of the second type have various denominations in literature (see Longuet-Higgins 1968). We call them Rossby waves in reference to the early studies by Rossby et al. (1939). These modes are defined for |ν| > 1 only. Most of them are characterized by negative and are identified in this work by subscripts n< 0, e.g. Λ_-1 ≥ Λ_-2 ≥ Λ_-3 ≥ ..., the parity of n corresponding to the parity of the associated function. Longuet-Higgins (1968) established the analytic expressions of the functions in the vicinity of the boundaries ν = ± 1 using Laguerre polynomials. In contrast with gravity modes, the oscillations of Rossby modes are concentrated around the poles in a region defined by 0 ≤ μ ≤ 1/|ν|. As pointed out by Lindzen (1966), Rossby modes have to be taken into account in the regime of sub-inertial waves, otherwise the set of eingenfunctions is not complete.

In Fig. (A.1), we can observe that the eigenvalues () strongly depend on the spin parameter. This dependence will affect the propagation of tidal waves along the vertical direction through their vertical wavenumbers. In addition, the Coriolis acceleration induces a dissymmetry between prograde (ν< 0) and retrograde (ν> 0) modes. Particularly, the gravity mode of lowest degree present two different asymptotic behaviours in the sub-inertial regime. The projection of the associated Legendre polynomial , which is representative of the semidiurnal tidal forcing, is also directly related to ν. We retrieve for the dissymmetry observed for the eigenvalues and note that the variations of the decomposition with ν are not simple to describe. In the present study, this decomposition is computed numerically with the spectral method presented below.

Appendix A.2: Spectral method

The early studies by Longuet-Higgins (1968) and CL70 provide a detailed theoretical approach of Hough functions and present the very powerful general method introduced by Hough (1898) to compute the eigenvalues and eigenvectors of ℒ^m,ν. In this method, which is not detailed here, Hough functions are expanded in series of normalized associated Legendre polynomials (Abramowitz & Stegun 1972). In the previous subsection, Eq. (A.1) highlights the symmetry between the two parameters defining these functions, m and ν. Hence, for m> 0 and ν ∈R, (A.6)Therefore, Hough functions and eigenvalues for m< 0 are readily deduced from those obtained for m> 0. This allows us to write (A.7)where s = |m|. A first possible way to compute the coefficients and eigenvalues , described by CL70, consists in finding the roots of a high-degree polynomial (typically, the degree N of this polynomial corresponds to the number of coefficients of the series given by Eq. (A.7)). With the current computing resources, N = 500 can be reached easily. However, this method rapidly becomes tricky to apply beyond N ~ 500 because of very large polynomial coefficients. Since the size of the range depends on N, this can limit the number of gravity modes or Rossby modes found. For instance, in the case of the Earth’s semidiurnal tide, where ν is very close to 1, the first Rossby modes are characterized by strongly negative eigenvalues. Therefore, only a few functions of this type can be computed. To address this issue, the previous problem may be written as an explicit eigenvalue problem. Hence, following Hough (1898) and CL70, we expand the non-normalized Hough functions in non-normalized associated Legendre Polynomials, denoted : (A.8)To lighten the expressions, the are simply denoted in the following analytic development. Let us introduce X = 1/ and the coefficients (A.9)Hough (1898) demonstrated the following equation, which corresponds to Eq. (60) in CL70 (A.10)with (A.11)The system of equations then reduces to the two eigenvalue problems (A.12)and (A.13)with the matrices(A.14)(A.15)They are solved numerically by linear algebra algorithms with great efficiency. For example, the case N = 1000 is computed on a laptop instantaneously. To make the comparison with Lee & Saio (1997) easier (see the previous subsection), we set at the equatorial point (θ = π/ 2) for even solutions and for odd solutions. Hence, the coefficients of Eq. (A.7) are explicitly given by (A.16)where (A.17)The solutions and projection matrices given by the spectral method are plotted in Figs. 2, 3, A.2, A.3. In the super-inertial regime (|ν| ≤ 1), the Hough functions basis is composed of gravity modes alone. Therefore, the projection matrix is a band-matrix. In this regime Hough functions are very similar to the associated Legendre polynomials. In the case ν = 0 (top left), the bases are exactly the same and the projection matrix is diagonal.

The sub-inertial regime is characterized by modes of both types, the transition being identified by Hough functions of the lowest degrees: Rossby modes on the left of the map (lowest eigenvalues) and gravity modes on the right (highest eigenvalues). The maps reveal the side effect inherent to the spectral method. Hough functions of the highest degrees are obviously not well represented by the family of associated Legendre Polynomials; more polynomials are necessary. This is explained by the asymptotic behaviour of Hough functions. At |ν| → 1⁺, Rossby modes are trapped around the poles and are thus well represented by highly oscillating associated Legendre polynomials (Fig. A.3). At |ν| → + ∞, gravity modes behave similarly and are trapped around the equatorial belt (Figs. 2, A.2).

Since the spectral method has been proved to be very efficient at computing a large number of eigenvalues and eigenfunctions at once, it is used for all calculations in this work. We note that some particular cases cannot be treated using Eq. (A.10). These cases have been discussed by CL70 and will not be developed in this work.

Appendix B: Numerical scheme used to integrate spatial differential equations

We detail here the numerical scheme used to integrate Laplace’s tidal equation (with the finite differences method) and the vertical structure equation. This method, introduced by Bruce et al. (1953) and described by Chapman & Lindzen (1970) for a particular case, is very convenient when solving second-order differential equations with non-constant coefficients. Those equations are written (B.1)where P, Q, and R are given functions of x. The domain on which Eq. (B.1) is integrated is divided into N intervals, i.e. N + 1 points, identified by the index n. Hence, for a regular mesh with elements of size δx, the first- and second-order derivative of a solution f are expressed (B.2)with n ∈ ˜1,N−1¨. Substituting Eq. (B.2) in Eq. (B.1), we obtain the numerical relations (B.3)with the coefficients (B.4)the P_n, Q_n, and R_n coefficients being the components of the numerical vectors corresponding to P, Q, and R. Then, we introduce the coefficients α_n and β_n, such that (B.5)These coefficients are defined straightforwardly by the recursive relations (B.6)The first terms (α₀ and β₀) are deduced from the boundary condition at x = x_inf, which can be written (B.7)where a₀, b₀, and c₀ are real or complex coefficients. Thus, we obtain (B.8)and, step by step, all the following terms. At x = x_sup, we apply the condition (B.9)with , which gives (B.10)The solution is finally obtained by computing the f_n backwards using Eq. (B.5).

Appendix C: Decomposition of the thermal forcing in Keplerian elements

We establish here the general multipole expansion in spherical harmonics of the thermal power given in Eq. (72) using the notations of Mathis & Le Poncin-Lafitte (2009). We consider a planet (A) orbiting circularly around its star (B). Three reference frames, all centred on the centre of mass of the planet and represented in Fig. C.1, must be defined:

• R_R:{O_A,X_R,Y_R,Z_R}, time independent inertial frame, with Z_R the direction of the total angular momentum of the system.

• R_O:{O_A,X_O,Y_O,Z_O}, the orbital frame linked to R_R by the Euler angles:

  • I_B, the inclination of the orbital frame with respect to ;

  • ω_B, the argument of the pericentre;

  • , the longitude of the ascending node.

• R_E;T:{O_A,X_E,Y_E,Z_E}, the spin equatorial frame rotating with the angular velocity Ω_A, and linked to R_R by three Euler angles:

  • ε_A, the obliquity, i.e. the inclination of the equatorial plane with respect to ;

  • Θ_A, the mean sidereal angle such that Ω_A = dΘ_A/ dt, which is the angle between the minimal axis of inertia and the straight line due to the intersection of the planes and ;

  • φ_A, the general precession angle.

Finally, we denote a the semi-major axis and the mean anomaly with , n_B being the mean motion. Neglecting the action of the centrifugal acceleration, the hydrostatic structure of the atmosphere presents a spherical symmetry, p₀, ρ₀, T₀ being functions only of the radial coordinate. Hence, any heat forcing caused by the star can be written (C.1)where h is a function of r and ψ, the angle between the points M(r,θ,ϕ) and B(a,θ_AB,ϕ_AB). We decompose J in normalized Legendre polynomials (C.2)with . Using the addition theorem (C.3)with (C.4)we obtain (C.5)Kaula’s transform is finally introduced to express as a function of Keplerian elements (Kaula 1962; Mathis & Le Poncin-Lafitte 2009). Since the motion of the planet is circular (e = 0), we get (C.6)with (j,p,q) ∈ ˜−l,l¨ × ˜0,l¨ ×Z, where the κ_l,j coefficients are given by (C.7)The functions denoted , F_l,j,p, and G_l,p,q are called the obliquity, inclination, and eccentricity functions (for more details, see Mathis & Le Poncin-Lafitte 2009). The angle Ψ_l,m,j,p,q can be written (C.8)the parameter σ_l,m,p,q = mΩ_A−(l−2p + q)n_B being the tidal frequency of the mode (l,m,p,q) and the phase lag. Substituting Eq. (C.6) in Eq. (C.5), we finally get(C.9)where (C.10)with ϕ_l,m,j,p,q = ψ_l,m,j,p,q/m.

With the same notations, the gravitational potential is written (C.11)The tidal frequency, written σ = mΩ−sn_orb, is defined by the doublet ⁽m,s⁾ ∈ Z².

The perturbing tidal potential and thermal power finally write (C.12)where the spatial functions U^m,σ and J^m,σ are expressed (C.13)with (C.14)the notation δ standing for the Kronecker symbol.

Appendix D: Radiated power as a function of physical parameters

The power per unit mass emitted by an optically thin atmosphere (J_rad) can be expressed as an analytic function of physical parameters in our simplified framework. Indeed, assuming that the flux radiated by a layer of thickness dx, denoted δI, propagates along the vertical direction without being absorbed by the other layers, we write (D.1)where φ_rad is the total power emitted per unit mass. Then, the atmosphere is treated as a grey body of molar emissivity per unit mass ϵ_a, which – using the molar concentration C₀ introduced in Sect. 2 – gives (D.2)Combining Eqs. (D.1) and (D.2) allows us to get (D.3)Finally, we derive the perturbed radiative loss J_rad associated with the tidal perturbation as a function of δT by linearizing the previous expression around the equilibrium state, (D.4)which gives the radiation frequency introduced in Eq. (17): (D.5)

Appendix E: Analytic solution for constant profiles and upper boundary condition

We show here that the solution of the equation describing the vertical structure strongly depends on the choice of the upper boundary condition. As discussed in Sect. 3, (E.1)at the ground level x = 0.

Instead of the energy-bounded condition prescribed in the model, at the upper limit (x_atm) we set the commonly used stress-free condition, (E.2)We obtain, at x = x_atm, (E.3)with the complex coefficients (E.4)The parameter x_atm must be fixed so that , which corresponds to x_atm ≫ 1. A simple analytic solution can be computed assuming that and are constant profiles. The second member of the vertical equation, C, given by Eq. (99), becomes (E.5)the solution is written (E.6)with and , and the complex integration constant (E.7)Now, the condition of Eq. (108) characterizing thermal tides, i.e. (E.8)is assumed. The second-order Love number reduced to the main contribution of the thermal semidiurnal tide, given by Eq. (76), can be expressed (E.9)

for a thin isothermal atmosphere. As a zero-order approximation, we consider that the forcings are the constant profiles U₂ = U_2,2,2,0,0 and J₂ = J_2,2,2,0,0 used in Sect. 2 (with J₂ ∈R). Therefore, introducing the complex impedance Ξ_atm, defined by (E.10)and substituting the analytic solution (Eq. E.6) in Eq. (117), we get (E.11)where Ξ_atm is explicitly given by (E.12)with (E.13)Similarly, the torque of Eq. (88) associated with the -contribution of the semidiurnal tide can be written (E.14)To illustrate the difference induced by the choice of the upper boundary condition, we compute the variation of the tidal torque generated by the Earth’s semidiurnal tide with the reduced tidal frequency (χ = (Ω−n_orb) /n_orb) for two conditions: the energy-bounded condition used in this work and the free-surface condition studied here. To simplify calculations, we focus on the gravity mode of lowest degree (n = 0). Results are plotted in Fig. E.1. In the first case, we obtain the smooth evolution observed in Fig. 14 (with slight differences in the vicinity of the synchronization because of the absence of Rossby modes), but the free-surface condition obviously leads to a highly resonant behaviour. This can be explained by the form of the solution given by Eq. (E.6), where the integration constant depends on and x_atm. By fixing x_atm, we implicitly define values of for which reaches a maximum. On the contrary, the integration constant obtained by applying the energy-bounded condition, given by Eq. (128), does not depend on the upper limit far beyond the critical altitude given by Eq. (127).

Appendix F: Supplementary materials regarding the dependence of the tidal response on the tidal frequency

As a complement to the tidal torques (Figs. 14) and perturbed pressure, density, and temperature (Figs. 17) resulting from the Earth’s solar semidiurnal tide, we plot several quantities as a function of the tidal frequency that can be helpful for the diagnostic of results (Fig. F.1). These quantities are the coefficients introduced in Eq. (81) which weight the contribution of Hough modes to the tidal torque; the eigenvalues of Laplace’s tidal equation (); the length scale of vertical variations associated to Hough modes (L_V;n); the spin parameter (ν); and the real and imaginary parts of the complex vertical wavenumber ( denoted k_n in plots), which determines the nature (propagative or evanescent) of tidal waves.

These plots highlight the difference between prograde (ν< 0) and retrograde (ν> 0) modes. Indeed, it can be noted that the projection of the associated Legendre polynomial on Hough functions is not the same before and after the synchronization: the weight of the gravity mode of degree 0, which is the most representative one, varies and Rossby modes are predominating in a very short interval of positive tidal frequencies. This dissymmetry is retrieved in the variations of the horizontal eigenvalues, where Rossby modes only appear for σ> 0.

As predicted by its analytic expression (Eq. (102)), the real and imaginary part of the vertical wavenumber diverge at the synchronization. Therefore, Hough modes become both highly oscillating and strongly damped in this region. Their spatial variation scale along the vertical direction decays, which corresponds to what is observed in the variations of L_V;n. Actually, below a given typical length scale, these modes are damped by dissipative processes such as thermal diffusion. Finally, thegraph of ν shows that this parameter strongly varies around the synchronization. The transition at σ = 0 is discontinuous, ν switching from − ∞ to + ∞.

Appendix G: Consequences of an insufficient spatial resolution on the numerical computation of the vertical structure

As pointed out by CL70 and detailed in Sect. 5, solving numerically the vertical structure equation (Eq. (107)) with the numerical scheme of Appendix B requires using a sufficiently high spatial resolution. This condition is mathematically expressed by the criterion given by Eq. (228) for a regular mesh and depends on the vertical wavenumbers that characterize Hough modes. To compute the response of a mode, it is necessary to use a spatial step δx lower than the corresponding variations length scale (L_V;n, defined by Eq. (48)), otherwise the signal does not appear. We illustrate this point in here by considering the simplified forced wave equation (G.1)where k_V ∈C is such that ℑ{k_V} > 0, with the previously chosen boundary conditions, (G.2)We arbitrarily set k_V = 100⁽1 + i⁾ in order to satisfy the condition L_V ≪ 1 (keeping in mind that x_sup = 10). This corresponds to a strongly damped tidal mode similar to those predicted by the model in a frequency range close to synchronization. The critical number of points necessary to compute the response with our scheme is then N_crit ~ 400. We solve Eq. (G.1) for various resolutions. The corresponding results are plotted in Fig. G.1 where they are also compared to the analytical solution given by Eq. (129). The flattening observed for N<N_crit gives an idea of the kind of numerical difficulties that need to be addressed to compute tidal low-frequency waves with a GCM. Particularly, this shows that if we compute the tidal torque numerically by using the scheme of Appendix B with an insufficient resolution, we will only be able to get the large-scale contribution associated with the horizontal component (see Figs. 14).

Appendix H: Nomenclature

All Tables

All Figures