Issue 
A&A
Volume 603, July 2017



Article Number  A107  
Number of page(s)  44  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201628252  
Published online  19 July 2017 
Atmospheric tides in Earthlike planets
^{1} IMCCE, Observatoire de Paris, CNRS UMR 8028, PSL, 77 Avenue DenfertRochereau, 75014 Paris, France
email: pierre.auclairdesrotour@ubordeaux.fr
^{2} Laboratoire AIM ParisSaclay, CEA/DSM – CNRS – Université Paris Diderot, IRFU/SAp Centre de Saclay, 91191 GifsurYvette, France
email: stephane.mathis@cea.fr
^{3} LESIA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC Univ. Paris 06, Univ. Paris Diderot, Sorbonne Paris Cité, 5 place Jules Janssen, 92195 Meudon, France
Received: 4 February 2016
Accepted: 8 November 2016
Context. Atmospheric tides can strongly affect the rotational dynamics of planets. In the family of Earthlike planets, which includes Venus, this physical mechanism coupled with solid tides makes the angular velocity evolve over long timescales and determines the equilibrium configurations of their spin.
Aims. Unlike the solid core, the atmosphere of a planet is subject to both tidal gravitational potential and insolation flux coming from the star. The complex response of the gas is intrinsically linked to its physical properties. This dependence has to be characterized and quantified for application to the wide variety of extrasolar planetary systems.
Methods. We develop a theoretical global model where radiative losses, which are predominant in slowly rotating atmospheres, are taken into account. We analytically compute the perturbation of pressure, density, temperature, and velocity field caused by a thermogravitational tidal perturbation. From these quantities, we deduce the expressions of atmospheric Love numbers and tidal torque exerted on the fluid shell by the star. The equations are written for the general case of a thick envelope and the simplified one of a thin isothermal atmosphere.
Results. The dynamics of atmospheric tides depends on the frequency regime of the tidal perturbation: the thermal regime near synchronization and the dynamical regime characterizing fastrotating planets. Gravitational and thermal perturbations imply different responses of the fluid, i.e. gravitational tides and thermal tides, which are clearly identified. The dependence of the torque on the tidal frequency is quantified using the analytic expressions of the model for Earthlike and Venuslike exoplanets and is in good agreement with the results given by global climate models (GCM) simulations.Introducing dissipative processes such as radiation regularizes the tidal response of the atmosphere, otherwise it is singular at synchronization.
Conclusions. We demonstrate the important role played by the physical and dynamical properties of a superEarth atmosphere (e.g. Coriolis, stratification, basic pressure, density, temperature, radiative emission) in its response to a tidal perturbation. We point out the key parameters defining tidal regimes (e.g. inertia, BruntVäisälä, radiative frequencies, tidal frequency) and characterize the behaviour of the fluid shell in the dissipative regime, which cannot be studied without considering the radiative losses.
Key words: hydrodynamics / waves / planets and satellites: atmospheres / planets and satellites: dynamical evolution and stability
© 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 Earthlike 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 nonnegligible fraction of its mass (about 20% in the case of small Neptunelike 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 secondorder Love number (k_{2}) introduced by Love (1911), which measures the quadrupolar hydrostatic elongation. The adiabatic Love number k_{2} 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 semidiurnal 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 Earthlike 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 spinorbit 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 biparameter 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 BeerLambert 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.
Fig. 1 Spherical coordinate system associated with the equatorial reference frame of the planet, and the geometrical parameters of the system. 
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_{0}, density ρ_{0}, and temperature T_{0}, 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 BruntVäisälä frequency (e.g. Lighthill 1978; Gerkema & Zimmerman 2008), is given by (1)where Γ_{1} = (∂lnp_{0}/∂lnρ_{0})_{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_{r}e_{r}; the corresponding displacement field ξ = ξ_{θ}e_{θ} + ξ_{ϕ}e_{ϕ} + ξ_{r}e_{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 sixequation system. We first introduce the NavierStokes equation. We assume the Cowling approximation in which the selfgravitational 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 NavierStokes 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 lefthand 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δρ/ρ_{0}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 BruntVä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 lefthand 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 BruntVä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 righthand 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 σ_{0}, 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_{0}, ρ_{0}, and T_{0}, 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_{0} = ρ_{0}/M, it is possible to express J_{rad} and σ_{0} 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 StefanBoltzmann 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^{1}. We also introduce the spin parameter (22)which defines the possible regimes of tidal gravitoinertial waves:

ν ≤ 1 corresponds to superinertial waves, for which the tidal frequency is greater than the inertia frequency;

ν > 1 corresponds to subinertial 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,σ}/ρ_{0}. Since V^{m,σ} = ∂_{t}ξ^{m,σ}, we also get the corresponding horizontal component of the displacement vector from Eqs. (25) and (26)
Fig. 2 Normalized Hough functions , , and corresponding to the Earth’s semidiurnal tide (m = 2 and ν = 1.0030). We note that ν = 2Ω /σ. Left: first gravity modes; 0 ≤ n ≤ 5. Right: first Rossby modes: − 6 ≤ n ≤ −1. In this case, ν is very close to 1 and the Rossby modes are trapped at the poles. The Earth’s semidiurnal tide is mainly described by gravity modes. For more details about the behaviour of Hough functions, see Fig. A.2. 
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)^{2}. 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.
Fig. 3 First gravity modes (0 ≤ n ≤ 5) of normalized Hough functions , , and corresponding to the semidiurnal tide of Venuslike exoplanets (m = 2 and ν = 0.4804). We note that ν = 2Ω /σ. For more details about the behaviour of Hough functions, see Figs. A.2. 
The system composed of Eqs. (31) to (33) can then be reduced to the pair of firstorder 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_{1}/ dr) /B_{1}, 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ödingerlike 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} = 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^{3}. 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 cutoff 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
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} = 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 mixedacoustic 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.
Fig. 4 Frequency spectrum of the atmospheric tidal response. The identified regimes are given as functions of the tidal frequency (σ). 
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 corotating 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 changeofbasis coefficients, (70)which allows us to write the atmospheric tidal potential (71)Let us treat the case of a simplified planetstar system, where the planet P circularly orbits around its host star, of mass M_{S}. The semimajor 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 secondorder 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^{)}^{4}(74)The tidal potential becomes (75)If the star can be considered a pointmass 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 secondorder 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 Earthlike planet, with g ~ 10 m s^{2} and R ~ 6.0 × 10^{3} 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_{0}, ρ_{0}, and T_{0} 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 BruntVä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_{1} (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 NavierStokes equation (Eq. (11)) and is negligible for usual values of (ε_{s;n} ≪ 1 and ς_{n} ∝ H^{2}/R^{2} ≪ 1). Therefore, by considering the case where σ_{0} ≪ σ, 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)
Fig. 5 Real and imaginary part of the vertical wavenumber as functions of the tidal (σ) and critical (σ_{c;n}) frequencies. Left: , with defined in Eq. (104), as a function of f^{(}σ/σ_{0}^{)} (horizontal axis) and (vertical axis), where f is the function defined by . The colour of the map, denoted c, is defined by . Hence, the purple region stands for slowly oscillating waves. The luminous (dark) regions correspond to strongly oscillating waves propagating downwards (upwards).Right: as a function of the same parameters. The luminous (dark) regions correspond to strongly (weakly) evanescent waves. The positions of the planets (Earth _{⊕}, Venus _{♀}) on the map are determined by the gravity mode of degree 0 in the case of the semidiurnal tide (Λ_{⊕;0} = 11.1596, Λ_{♀;0} = 7.1485). The tidal frequencies of the planets are given by σ_{⊕}/σ_{0} ≈ 194 and σ_{♀}/σ_{0} = 0.43, with σ_{0} = 7.5 × 10^{7} s^{1}. 
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 σ ≈ σ_{0}.
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 righthand 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)^{5}, (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 Earthlike 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 freesurface 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 waveguide 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 frequencyresonant. 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 freesurface 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)^{6}.
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 threedimensional 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 secondorder 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 secondorder 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_{2}/U_{2} and thus depend on the properties of the whole starplanet system, particularly the semimajor 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_{2} ≈ 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_{2} ∈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^{2}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 Venuslike 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 Jupiterlike 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 BruntVä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)
Fig. 6 Relative difference between the complex wavenumber obtained in this work (Eq. (103)) on the one hand, and the real one given by CL70 for Earth atmospheric tides ( ) on the other hand. The difference η, given by Eq. (158), is plotted in logarithmic scale as a function of the normalized tidal frequency σ/σ_{0} and of the ratio h_{c}/h_{n}. The colour of the map, denoted c, is defined by the expression c = log (η). The horizontal and vertical axis represent f^{(}σ/σ_{0}^{)} and , where f is the function defined by with X ∈R. The positions of the planets are determined by the gravity mode of degree 0 in the case of the semidiurnal tide. Dark (luminous) regions correspond to the dynamic (dissipative) regime. The vertical wavenumbers of Venus are strongly affected by dissipative processes, while those of the Earth are not. 
This difference, written (158)is plotted in Fig. 6 as a function of the reduced tidal frequency σ/σ_{0}, 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} = 0. This thermal regime corresponds to the middle whiteyellow region of the map. The two wavenumbers become similar for σ ≫ σ_{0}. 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 NavierStokes 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 Venuslike 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 BruntVä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 BruntVäisälä frequency can be expressed as(159)Therefore, the case N^{2} = 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 NavierStokes 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 NavierStokes 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_{0} the mean absorbed power per unit mass. At the upper boundary x_{atm} ≫ 1, we apply the freesurface 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. Secondorder 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 secondorder 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 Venuslike 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 nonsynchronized 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.
Fig. 7 Tidal torques exerted on the atmosphere in the case of the slowly rotating convective atmosphere (red line with the label “Conv”, Eq. (174)) and in the model by Correia & Laskar (2001) (green dashed line with the label “CL01”, Eq. (175)) as functions of the tidal frequency. Both functions are normalized by their respective maximum value, reached at σ = σ_{0}; the tidal frequency is normalized by the Newtonian cooling frequency (σ_{0}). The tidal forcing is quadrupolar, with U and J given by Sect. 4.3. The parameters a and b of Eq. (175) are adjusted so that the peaks of the two functions are superposed. 
5. Physical description of heat source terms
On the lefthand 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^{7}. 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 nonabsorbed 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 starplanet 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 starplanet 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_{0} by the atmosphere, we use the BeerLambert 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)
Fig. 8 Left: absorption of the flux emitted by the star along the starplanet direction. Right: absorption of the flux emitted by the ground. 
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_{0} = ρ_{0}/M. Therefore, the incident power flux and heat power per unit mass become (185)Near the planetstar 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.
Fig. 9 Vertical profiles of radiative fluxes (I), power per unit volume (J), and power per unit mass (J) for various values of the dimensionless optical depth τ_{λ} (Eq. (187)). As usual, the reduced altitude, x = z/H, is on the vertical axis. Left: incident flux, defined by Eq. (186), at ψ = 0. If τ_{λ} ≪ 1 the ground receives almost 100% of the total power. On the contrary, if τ_{λ} ≫ 1, all the power is absorbed by the atmosphere. Right: radiative flux emitted by the ground, defined by Eq. (206). For small τ_{λ}, the flux goes through the atmosphere without being absorbed. For high values of the parameter, it is absorbed by the lowest layers. 
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 starplanet 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)
Fig. 10 Left, from top to bottom: power per unit mass (J^{inc}, Eq. (184)) and power per unit volume (J^{inc}, Eq. 183) of a monochromatic incident flux normalized by their maxima and plotted as functions of the normalized spatial coordinates X (the direction of the star is along the horizontal axis) and Z (the spin axis is along the vertical axis). The planet is represented by a green disk of unit radius. The flux, of wavelength λ, propagates in an atmosphere assumed to be homogeneous in composition (C_{G} = C_{0}) and characterized by the ratio H/R = 0.1. The optical depth is given by τ_{λ} = 0.2. Right, from top to bottom: power per mass unit (J^{rad}, Eq. (206)) and power per volume unit (J^{rad}) of a monochromatic flux radiated by the ground, normalized by their maxima, and plotted as functions of X and Z, with τ_{λ} = 0.5. 
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 StefanBoltzmann law, (193)the parameter ϵ_{gr} being the emissivity of the ground and S= 5.670373 × 10^{8} W m^{2} K^{4} the StefanBoltzmann 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 counterradiation can be assumed to be proportional to I_{gr;BB}, and expressed by the semiempirical formula (194)which allows us to take it into account without studying the whole coupled system groundatmosphere. 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 frequencydependent 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)
Fig. 11 Nyquist diagram of the boundary layer transfer function in the frequency range − 10^{3} ≤ σ/σ_{bl} ≤ 10^{3} (σ being the tidal frequency and σ_{bl} the frequency characterizing the boundary layer introduced in Eq. (217)). The parametrized complex transfer function (Eq. (215)) reduced by its value at σ = 0 is plotted as a function of the tidal frequency (σ) in the complex plane, where the horizontal axis corresponds to the real part of the function, , and the vertical axis to the imaginary part . The arrow indicates the direction of increasing σ. The gain of the transfert function is given by the distance of a point P of the curve to the origin O of coordinates , and its phase by the angle of the vector OP to the horizontal axis. This plot shows that the response of the ground is always delayed with respect to the thermal forcing except at synchronization, where and the gain is maximum. It also highlights the two identified regimes of the ground response: the lag tends to 0 (response in phase with the perturbation) and the gain increases at σ → 0, while the lag tends to ± π/ 4 and the gain to 0 at σ → ± ∞, the critical transition frequency being σ_{bl}. 
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^{4} 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^{2} s^{1}, which had been prescribed before by Wilkes (1949) and which leads to σ_{bl} ~ 10^{6} s^{1}.
Fig. 12 Bichromatic case. Top: heat transfers between the atmosphere, space and the ground. Bottom: power balance of the interface groundatmosphere. 
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 wavelengthdependent 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_{0} 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 Venuslike 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 Venuslike 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^{8} is given by σ = 2(Ω−n_{orb}). One of the most crucial parameters of the model is the Newtonian cooling frequency (σ_{0}). 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 nondissipative (Earth) and dissipative (Venus) regimes. Therefore, following Lindzen & McKenzie (1967), we consider the interesting case where the radial profile σ_{0} is assumed to be constant. For both planets, we arbitrarily set σ_{0} = 7.5 × 10^{7} s^{1}, which is the effective value of σ_{0} computed by Leconte et al. (2015) with a GCM for Venus. This value is such that σ ~ σ_{0} for Venus and σ ≫ σ_{0} 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_{2} and J_{2} are fixed constant radial profiles. The tidal potential U_{2} is computed from the Kaula multipole expansion detailed in Appendix D, and explicitly given by the expression of Eq. (74). A zeroorder approximation of the thermal power J_{2} 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.
Fig. 13 Average daily variation in temperature in degrees (top) and pressure in mbar (bottom) over one year, plotted over solar time (in hours), the subsolar point corresponding to 12 h. The data was recorded every minute over the full year 2013 with a Vantage Pro 2 weather station at latitude 48.363°N. The time is converted to solar time and the data averaged over an exact count of 365 days. The average values of temperature (10.8973° C) and pressure (1015.83 mbar) have been removed to display only the variable part. In both figures, the raw data (in green) has been decomposed in a Fourier series over the 24 h period, limited to two harmonics (blue), or five harmonics (red). Temperature is dominated by the diurnal component, with a maximum value around 14h26mn. The largest harmonics of the pressure variations is the semidiurnal term, followed by the diurnal term. Two pressure maxima can be observed, at 9h38mn and 22h07mn. 
Fig. 14 Tidal torque applied to the atmosphere of Earth (left panel) and Venus (right panel) by the thermal semidiurnal tide and its components as functions of the reduced frequency χ = (Ω−n_{orb}) /n_{orb} (horizontal axis). The forcings consist of academic quadrupolar perturbations expressed as and with U_{2} = −0.985 m^{2} s^{2} and J_{2} = 1.0 × 10^{2} m^{2} s^{3} for the Earth and U_{2} = −2.349 m^{2} s^{2} and J_{2} = 1.0 × 10^{4} m^{2} s^{3} for Venus. Labels T, H, and V refer to the total response and its horizontal and vertical components, respectively (we note that ). The label CL01 refers to the model by Correia & Laskar (2001), given by Eq. (175). 
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.
Values of physical parameters used in simulations.
Fig. 15 Total tidal response of the Earth’s atmosphere to the solar semidiurnal perturbation caused by the quadrupolar academic forcings and with U_{2} = −0.985 m^{2} s^{2} and J_{2} = 1.0 × 10^{2} m^{2} s^{3}. Left, from top to bottom: amplitudes of δp, δρ, δT, V_{θ}, V_{ϕ}, and V_{r} in logarithmic scales as functions of the reduced latitude δ/π (horizontal axis) and altitude x = z/H in logarithmic scale (vertical axis). The colour function c is given by c = log (δf) for any quantity δf. Right, from top to bottom: Sine of the arguments of the same quantities as functions of the reduced latitude and altitude. In these plots, c = sin[arg(δf)]. 
Fig. 16 Total tidal response of the Earth’s atmosphere to the solar semidiurnal perturbation caused by the quadrupolar academic forcings and with U_{2} = −0.985 m^{2} s^{2} and J_{2} = 1.0 × 10^{2} m^{2} s^{3}. Left, from top to bottom: real parts of δp, δρ, δT, V_{θ}, and V_{ϕ} taken at the ground as functions of the reduced longitude ϕ/π (horizontal axis) and latitude δ/π (vertical axis). The colour function c is given by c = ℜ{δf} for any quantity δf. Right, from top to bottom: imaginary parts of the same quantities as functions of the reduced longitude and latitude. In these plots, c = ℑ{δf}. 
Fig. 17 Tidal surface oscillations of pressure, density, temperature, and their components at the equator of the Earth as functions of the reduced tidal frequency χ = (Ω−n_{orb}) /n_{orb} (horizontal axis). The perturbation is defined by the academic quadrupolar forcings and with U_{2} = −0.985 m^{2} s^{2} and J_{2} = 1.0 × 10^{2} m^{2} s^{3}. Labels T, H, and V refer to the total response and its horizontal and vertical components, respectively. 
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.
Fig. 18 Total tidal response of the atmosphere of Venus to the solar semidiurnal perturbation caused by the quadrupolar academic forcings and with U_{2} = −2.349 m^{2} s^{2} and J_{2} = 1.0 × 10^{4} m^{2} s^{3}. Left, from top to bottom: amplitudes of δp, δρ, δT, V_{θ}, V_{ϕ}, and V_{r} in logarithmic scales as functions of the reduced latitude δ/π (horizontal axis) and altitude x = z/H in logarithmic scale (vertical axis). The colour function c is given by c = log (δf) for any quantity δf. Right, from top to bottom: Sine of the arguments of the same quantities as functions of the reduced latitude and altitude. In these plots, c = sin[arg(δf)]. 
Fig. 19 Total tidal response of the atmosphere of Venus to the solar semidiurnal perturbation caused by the quadrupolar academic forcings and with U_{2} = −2.349 m^{2} s^{2} and J_{2} = 1.0 × 10^{4} m^{2} s^{3}. Left, from top to bottom: real parts of δp, δρ, δT, V_{θ}, and V_{ϕ} taken at the ground as functions of the reduced longitude ϕ/π (horizontal axis) and latitude δ/π (vertical axis). The colour function c is given by c = ℜ{δf} for any quantity δf. Right, from top to bottom: imaginary parts of the same quantities as functions of the reduced longitude and latitude. In these plots, c = ℑ{δf}. 
Fig. 20 Tidal surface oscillations of pressure, density, temperature, and their components at the equator of Venus as functions of the reduced tidal frequency χ = (Ω−n_{orb}) /n_{orb} (horizontal axis). The perturbation is defined by the academic quadrupolar forcing and with U_{2} = −2.349 m^{2} s^{2} and J_{2} = 1.0 × 10^{4} m^{2} s^{3}. Labels T, H, and V refer to the total response and its horizontal and vertical components, respectively. 
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^{2} ≈ 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^{37} kg m^{2} (NASA fact sheets^{10}), 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 Earthlike 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^{20} kg and M_{⊕} ~ 5.1 × 10^{18} 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_{2} ~ 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 Venuslike 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 σ ~ σ_{0}, 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 superinertial 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 crosssections 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^{37} kg m^{2} (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_{2} ~ 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, commonly used in the literature of geophysics, is very usefulbecause it allows us to separate the horizontal and vertical dependences in the NavierStokes 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 superinertial 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 gravitoinertial 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^{2} ≈ 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 σ_{0}, 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).
Fig. 21 Asymptotic domains of parameters where the analytical approach presented in this work can be applied (green) and where the numerical approach is necessary (red) as functions of the frequency ratios N^{2}/σ^{2} (horizontal axis) and (2Ω)^{2}/σ^{2} (vertical axis). Green regions: stably stratified atmospheres satisfying the hierarchy of frequencies (2Ω)^{2} ≪ N^{2} necessary to apply the traditional approximation a); weakly stratified or convective atmospheres of slowly rotating planets b). Red regions: stably stratified atmospheres where the hierarchy of frequencies (2Ω)^{2} ≪ N^{2} is not satisfied c); superadiabatic rotating convective atmospheres d). 
Fig. 22 Predictive aspects related to the dependence of the atmospheric tidal torque on stratification. The neutral convective case implies a strong torque, which can lead to nonsynchronized states of equilibrium. The stably stratified case implies a weak torque, which lets the planet evolve towards spinorbit synchronization. 
7. Conclusions
We have set the bases of a new theoretical framework to study the atmospheric tides of Earthlike 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,BruntVä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 nondissipative Earthlike dynamical regime characterizing fast rotating bodies, and a dissipative Venuslike 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 freesurface condition and a nomaterialescape 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 nonsynchronized 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.
The properties of Hough functions and the method used to compute them in this work are detailed in Appendix A.
To solve the vertical structure equation numerically in any case, we use the method described in Appendix B.
The parameters j, p, and q are the usual indexes of the Kaula’s expansion of the tidal gravitational potential; see Appendix C.
These coefficients are the same as those given by Eq. (60) to an H factor.
In the case of the semidiurnal tide treated by CL70 (no dissipative processes), this condition cannot be applied because . To address this particular situation, CL70 proposes applying a radiation condition at the upper boundary, i.e. considering that the energy is only propagating upwards and eliminating the terms of the solution corresponding to an energy propagating downwards. However, in the present work, when σ_{0}> 0, , so we do not have to consider this case.
For instance, see the Kaula expansion of the gravitational tidal potential (Kaula 1964) detailed in Appendix C and Mathis & Le PoncinLafitte (2009).
This frequency corresponds to the term defined by (l,m,j,p,q) = ^{(}2,2,2,0,0^{)} in the multipole expansion of the forcings detailed in Appendix C.
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 superEarth atmospheres from the point of view of numerical simulations. P. AuclairDesrotour 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 CEASaclay.
References
 Abramowitz, M., & Stegun, I. A. 1972, Handbook of Mathematical Functions (New York: Dover) [Google Scholar]
 Arras, P., & Socrates, A. 2010, ApJ, 714, 1 [Google Scholar]
 Auclair Desrotour, P., Mathis, S., & Le PoncinLafitte, C. 2015, A&A, 581, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Avduevsky, V. S., Marov, M. Y., Noykina, A. I., Polezhaev, V. I., & Zavelevich, F. S. 1970, J. Atmos. Sci., 27, 569 [NASA ADS] [CrossRef] [Google Scholar]
 Baker, R. D., Schubert, G., & Jones, P. W. 2000, J. Atmos. Sci., 57, 184 [NASA ADS] [CrossRef] [Google Scholar]
 Beer 1852, Annalen der Physik, 162, 78 [NASA ADS] [CrossRef] [Google Scholar]
 Bernard, E. A. 1962, Archives for Meteorology Geophysics and Bioclimatology Series A Meteorology and Atmopsheric Physics, 12, 502 [Google Scholar]
 Bouguer, P. 1729, Essai d’Optique sur la gradation de la Lumière (Paris: C. Jombert) [Google Scholar]
 Bruce, G. H., Peaceman, D., Rachford Jr, H., et al. 1953, Journal of Petroleum Technology, 5, 79 [CrossRef] [Google Scholar]
 Chapman, S., & Lindzen, R. 1970, Atmospheric tides. Thermal and gravitational (Dordrecht: Reidel) [Google Scholar]
 Correia, A. C. M., & Laskar, J. 2001, Nature, 411, 767 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Correia, A. C. M., & Laskar, J. 2003, J. Geophys. Res. (Planets), 108, 5123 [Google Scholar]
 Correia, A. C. M., Laskar, J., & de Surgy, O. N. 2003, Icarus, 163, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Correia, A. C. M., Levrard, B., & Laskar, J. 2008, A&A, 488, L63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Correia, A. C. M., Boué, G., Laskar, J., & Rodríguez, A. 2014, A&A, 571, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Covey, C., Dai, A., Marsh, D. R., & Lindzen, R. S. 2009, AGU Fall Meeting Abstracts [Google Scholar]
 Cowling, T. G. 1941, MNRAS, 101, 367 [NASA ADS] [Google Scholar]
 Dickinson, R. E., & Geller, M. A. 1968, J. Atmos. Sci., 25, 932 [NASA ADS] [CrossRef] [Google Scholar]
 Dobrovolskis, A. R., & Ingersoll, A. P. 1980, Icarus, 41, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Eckart, C. 1960, Physics of Fluids, 3, 421 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Efroimsky, M., & Lainey, V. 2007, J. Geophys. Res. (Planets), 112, 12003 [Google Scholar]
 Egbert, G. D., & Ray, R. D. 2000, Nature, 405, 775 [NASA ADS] [CrossRef] [Google Scholar]
 Fabrycky, D. C., Ford, E. B., Steffen, J. H., et al. 2012, ApJ, 750, 114 [Google Scholar]
 Gastineau, M., & Laskar, J. 2014, TRIP 1.3.8, TRIP Reference manual, IMCCE, Paris Observatory, http://www.imcce.fr/trip/ [Google Scholar]
 Gerkema, T., & Shrira, V. I. 2005, J. Fluid Mech., 529, 195 [NASA ADS] [CrossRef] [Google Scholar]
 Gerkema, T., & Zimmerman, J. 2008, Lecture Notes, Royal NIOZ, Texel [Google Scholar]
 Gold, T., & Soter, S. 1969, Icarus, 11, 356 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Nicholson, P. D. 1989, ApJ, 342, 1075 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Soter, S. 1966, Icarus, 5, 375 [Google Scholar]
 Green, J. S. A. 1965, Proc. of the Royal Society of London, Series A, Mathematical and Physical Sciences, 288, 564 [CrossRef] [Google Scholar]
 Greenberg, R. 2009, ApJ, 698, L42 [Google Scholar]
 Hagan, M., Forbes, J., & Richmond, A. 2003, Encyclopedia of Atmospheric Sciences, 1, 159 [CrossRef] [Google Scholar]
 Haurwitz, B. 1965, Archives for Meteorology Geophysics and Bioclimatology Series A Meteorology and Atmopsheric Physics, 14, 361 [NASA ADS] [Google Scholar]
 Henning, W. G., O’Connell, R. J., & Sasselov, D. D. 2009, ApJ, 707, 1000 [NASA ADS] [CrossRef] [Google Scholar]
 Hough, S. S. 1898, Phil. Trans. R. Soc. Lond. Ser. A, 191, 139 [NASA ADS] [CrossRef] [Google Scholar]
 Hut, P. 1980, A&A, 92, 167 [NASA ADS] [Google Scholar]
 Hut, P. 1981, A&A, 99, 126 [NASA ADS] [Google Scholar]
 Ingersoll, A. P., & Dobrovolskis, A. R. 1978, Nature, 275, 37 [NASA ADS] [CrossRef] [Google Scholar]
 Ioannou, P. J., & Lindzen, R. S. 1993a, ApJ, 406, 252 [NASA ADS] [CrossRef] [Google Scholar]
 Ioannou, P. J., & Lindzen, R. S. 1993b, ApJ, 406, 266 [NASA ADS] [CrossRef] [Google Scholar]
 Ioannou, P. J., & Lindzen, R. S. 1994, ApJ, 424, 1005 [NASA ADS] [CrossRef] [Google Scholar]
 Kato, S. 1966, J. Geophys. Res., 71, 3201 [NASA ADS] [CrossRef] [Google Scholar]
 Kaula, W. M. 1962, AJ, 67, 300 [NASA ADS] [CrossRef] [Google Scholar]
 Kaula, W. M. 1964, Reviews of Geophysics and Space Physics, 2, 661 [NASA ADS] [CrossRef] [Google Scholar]
 Kelvin, L. 1863, Phil. Trans. Roy. Soc. Lond., Treatise on Natural Philosophy, 2, 837 [Google Scholar]
 Kelvin, L. 1882, Proc. Roy. Soc. Edinb., 11, 396 [CrossRef] [Google Scholar]
 Klett, Witwe, E., Detleffsen, Peter, C., et al. 1760, IH Lambert... Photometria sive de mensura et gradibus luminis, colorum et umbrae (sumptibus viduae Eberhardi Klett) [Google Scholar]
 Konopliv, A. S., & Yoder, C. F. 1996, Geophys. Res. Lett., 23, 1857 [NASA ADS] [CrossRef] [Google Scholar]
 Lacis, A. A. 1975, J. Atmos. Sci., 32, 1107 [NASA ADS] [CrossRef] [Google Scholar]
 Lainey, V., Dehant, V., & Pätzold, M. 2007, A&A, 465, 1075 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lamb, H. 1911, Proc. Roy. Soc. Lond. Ser. A, 84, 551 [NASA ADS] [CrossRef] [Google Scholar]
 Lamb, H. 1932, Hydrodynamics (New York: Dover) [Google Scholar]
 Laplace, P. S. 1798, Traité de mécanique céleste (Duprat J. B. M.) [Google Scholar]
 Leconte, J., Wu, H., Menou, K., & Murray, N. 2015, Science, 347, 632 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, U., & Saio, H. 1997, ApJ, 491, 839 [NASA ADS] [CrossRef] [Google Scholar]
 Lighthill, J. 1978, Waves in fluids (Cambridge: Cambridge University Press) [Google Scholar]
 Lindzen, R. S. 1966, Monthly Weather Review, 94, 295 [NASA ADS] [CrossRef] [Google Scholar]
 Lindzen, R. S. 1967a, Nature, 215, 1260 [NASA ADS] [CrossRef] [Google Scholar]
 Lindzen, R. S. 1967b, Quart. J. Roy. Meteor. Soc., 93, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Lindzen, R. S. 1968, Proc. Roy. Soc. Lond. Ser. A, 303, 299 [NASA ADS] [CrossRef] [Google Scholar]
 Lindzen, R. S., & McKenzie, D. J. 1967, Pure and Applied Geophysics, 66, 90 [NASA ADS] [CrossRef] [Google Scholar]
 Linkin, V. M., Kerzhanovich, V. V., Lipatov, A. N., et al. 1986, Science, 231, 1417 [NASA ADS] [CrossRef] [Google Scholar]
 LonguetHiggins, M. S. 1968, Phil. Trans. R. Soc. Lond. Ser. A, 262, 511 [NASA ADS] [CrossRef] [Google Scholar]
 Love, A. E. H. 1911, Some Problems of Geodynamics (Cambridge University Press) [Google Scholar]
 Marov, M. Y., Avduevskij, V. S., Kerzhanovich, V. V., et al. 1973, J. Atmos. Sci., 30, 1210 [NASA ADS] [CrossRef] [Google Scholar]
 Mathis, S. 2009, A&A, 506, 811 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., & Le PoncinLafitte, C. 2009, A&A, 497, 889 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., Talon, S., Pantillon, F.P., & Zahn, J.P. 2008, Sol. Phys., 251, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Mayor, M., & Queloz, D. 1995, Nature, 378, 355 [Google Scholar]
 Mignard, F. 1979, Moon and Planets, 20, 301 [Google Scholar]
 Mignard, F. 1980, Moon and Planets, 23, 185 [NASA ADS] [CrossRef] [Google Scholar]
 Mohr, P. J., Taylor, B. N., & Newell, D. B. 2012, J. Phys. Chem. Ref. Data, 41, 043109 [NASA ADS] [CrossRef] [Google Scholar]
 Neron de Surgy, O., & Laskar, J. 1997, A&A, 318, 975 [NASA ADS] [Google Scholar]
 Ogilvie, G. I. 2014, ARA&A, 52, 171 [Google Scholar]
 Ogilvie, G. I., & Lin, D. N. C. 2004, ApJ, 610, 477 [NASA ADS] [CrossRef] [Google Scholar]
 Pekeris, C. L. 1937, Proc. Roy. Soc. Lond. Ser. A, 158, 650 [NASA ADS] [CrossRef] [Google Scholar]
 Perryman, M. 2011, The Exoplanet Handbook (Cambridge University Press) [Google Scholar]
 Planck, M. 1901, Annalen der Physik, 309, 553 [NASA ADS] [CrossRef] [Google Scholar]
 Pollack, J. B., & Young, R. 1975, J. Atmos. Sci., 32, 1025 [NASA ADS] [CrossRef] [Google Scholar]
 Prat, V., Lignières, F., & Ballot, J. 2016, A&A, 587, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Press, W. H. 1981, ApJ, 245, 286 [NASA ADS] [CrossRef] [Google Scholar]
 Press, W. H., Flannery, B. P., & Teukolsky, S. A. 1986, Numerical recipes. The art of scientific computing (Cambridge University Press) [Google Scholar]
 Ray, R. D., Eanes, R. J., & Lemoine, F. G. 2001, Geophys. J. Int., 144, 471 [NASA ADS] [CrossRef] [Google Scholar]
 Remus, F., Mathis, S., Zahn, J.P., & Lainey, V. 2012, A&A, 541, A165 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rossby, C.G., et al. 1939, Journal of Marine Research, 2, 38 [Google Scholar]
 Seiff, A., Kirk, D. B., Young, R. E., et al. 1980, J. Geophys. Res.: Space Phys., 85, 7903 [NASA ADS] [CrossRef] [Google Scholar]
 Shen, M., & Zhang, C. Z. 1990, Icarus, 85, 129 [NASA ADS] [CrossRef] [Google Scholar]
 Siebert, M. 1961, Advances in Geophysics, 7, 105 [NASA ADS] [CrossRef] [Google Scholar]
 Taylor, G. I. 1929, Proc. Roy. Soc. Lond. Ser. A, 126, 169 [NASA ADS] [CrossRef] [Google Scholar]
 Taylor, G. I. 1936, Proc. Roy. Soc. Lond. Ser. A, 156, 318 [NASA ADS] [CrossRef] [Google Scholar]
 Tort, M., & Dubos, T. 2014, Quart. J. Roy. Meteor. Soc., 140, 2388 [NASA ADS] [CrossRef] [Google Scholar]
 Townsend, R. H. D. 2003, MNRAS, 340, 1020 [NASA ADS] [CrossRef] [Google Scholar]
 Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1989, Nonradial oscillations of stars (Tokyo: University of Tokyo Press) [Google Scholar]
 Webb, D. J. 1980, Geophys. J., 61, 573 [NASA ADS] [CrossRef] [Google Scholar]
 White, A. A., Hoskins, B. J., Roulstone, I., & Staniforth, A. 2005, Quart. J. Roy. Meteor. Soc., 131, 2081 [NASA ADS] [CrossRef] [Google Scholar]
 Wilkes, M. V. 1949, Oscillations of the Earth’s Atmosphere (University Press Cambridge) [Google Scholar]
 Williams, J. G., Konopliv, A. S., Boggs, D. H., et al. 2014, J. Geophys. Res. (Planets), 119, 1546 [NASA ADS] [CrossRef] [Google Scholar]
 Zahn, J. P. 1966, Annales d’Astrophysique, 29, 313 [Google Scholar]
 Zahn, J.P., Talon, S., & Matias, J. 1997, A&A, 322, 320 [Google Scholar]
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 (LonguetHiggins 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), LonguetHiggins (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 subinertial 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):
Fig. A.1 Top: eigenvalues () of the gravity and Rossby modes of lowest degrees (− 6 ≤ n ≤ 5) as functions of the spin parameter (ν). The horizontal axis represents the function of ν given by , and the vertical axis the function of given by . Bottom: projection coefficients of the normalized Legendre polynomial on normalized Hough functions associated with these modes as functions of the spin parameter (see Eq. (70)). 
(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 SturmLiouville 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 subinertial 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 subinertial regime. To construct the set of Hough functions, we use as predictors the asymptotic results given by LonguetHiggins (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 nonrotating 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 subinertial 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 LonguetHiggins 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. LonguetHiggins (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 subinertial 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 subinertial 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.
Fig. A.2 Normalized Hough functions computed as functions of the colatitude (θ) for 0 ≤ n ≤ 5 (gravity modes) and various values of ν = (2Ω) /σ. Topleft: ν = 0; topright: ν = 0.8; middleleft: ν = 1.2; middleright: ν = 2; bottomleft: ν = 4; bottomright: ν = 12. The case ν = 0 corresponds to the normalized associated Legendre polynomials . The two plots at the top represent the regime of superinertial waves (ν ≤ 1) where Hough functions look like Legendre polynomials. In the regime of subinertial waves (ν > 1), functions become more and more evanescent when ν increases, and oscillations are trapped around the equator in the interval [θ_{ν},π−θ_{ν}], where . 
Fig. A.3 Normalized Hough functions computed as functions of the colatitude (θ) for − 6 ≤ n ≤ −1 (Rossby modes) and various values of ν = (2Ω) /σ. Topleft: ν = 1.2; topright: ν = 2; bottomleft: ν = 4; bottomright: ν = 12. In the regime of subinertial waves (ν > 1), Rossby modes are trapped around the poles when ν → 1. 
Appendix A.2: Spectral method
The early studies by LonguetHiggins (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 highdegree 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 nonnormalized Hough functions in nonnormalized 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 superinertial regime (ν ≤ 1), the Hough functions basis is composed of gravity modes alone. Therefore, the projection matrix is a bandmatrix. 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 subinertial 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 secondorder differential equations with nonconstant 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 secondorder 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 (α_{0} and β_{0}) are deduced from the boundary condition at x = x_{inf}, which can be written (B.7)where a_{0}, b_{0}, and c_{0} 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 PoncinLafitte (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.

Fig. C.1 Inertial reference, orbital and equatorial rotating frames, and the associated Euler’s angles of orientation. Figure taken from Mathis & Le PoncinLafitte (2009). 
Finally, we denote a the semimajor 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_{0}, ρ_{0}, T_{0} 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 PoncinLafitte 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 PoncinLafitte 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^{2}.
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_{0} 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 energybounded condition prescribed in the model, at the upper limit (x_{atm}) we set the commonly used stressfree 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 secondorder Love number reduced to the main contribution of the thermal semidiurnal tide, given by Eq. (76), can be expressed (E.9)
Fig. E.1 Tidal torques exerted on the atmosphere (in petaNewton metres (PN.m)) induced by the first gravity mode of the Earth’s semidiurnal tide (n = 0) as a function of the reduced tidal frequency (χ) obtained with two different conditions at the upper boundary: the energybounded condition chosen for the model (red line) and a stressfree condition (green dashed line). The equation for the structure in the vertical (resp. horizontal) direction has been integrated with 10 000 mesh points (resp. 400 associated Legendre polynomials). 
for a thin isothermal atmosphere. As a zeroorder approximation, we consider that the forcings are the constant profiles U_{2} = U_{2,2,2,0,0} and J_{2} = J_{2,2,2,0,0} used in Sect. 2 (with J_{2} ∈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 energybounded condition used in this work and the freesurface 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 freesurface 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 energybounded 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
Fig. F.1 Tidal response of the Earth’s atmosphere to a quadrupolar forcing ( and ) as a function of the reduced tidal frequency χ = (Ω−n_{orb})/n_{orb} (horizontal axis). For each quantity q_{n} associated with the Hough mode of degree n, except the coefficients (top left) and the scales of variations L_{V;n}, the function f(χ) = sign(q_{n})log (1 + q_{n}) is plotted (vertical axis) for n ∈ ˜−6,5¨. Left, from top to bottom: a) coefficients (Eq. (70)); b) eigenvalues of Laplace’s tidal equation ( defined by Eq. (35)); c) normalized variation scales of the vertical structure of the vertical component, L_{V;n} = 2π/ (Eq. (48)), plotted in logarithmic scale. Right, from top to bottom: d) real part of the vertical wavelengths ( , Eq. (103)); e) imaginary part of the vertical wavelengths; f) spin parameter ν = 2Ω /σ (Eq. (23)). 
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 lowfrequency 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 largescale contribution associated with the horizontal component (see Figs. 14).
Fig. G.1 Real part of the solution of Eq. (G.1) as a function of x in logarithmic scale computed numerically for various mesh resolutions. The number of points of the spatial discretization is set at N = 20 (green dashed line), N = 200 (blue dashed line), N = 2000 (pink dashed line), and N = 20 000 (cyan dashed line). These numerical solutions are compared to the analytical solution given by Eq. (129) (red line). 
Appendix H: Nomenclature
List of the symbols used throughout this work and their designations in order of introduction.
All Tables
List of the symbols used throughout this work and their designations in order of introduction.
All Figures
Fig. 1 Spherical coordinate system associated with the equatorial reference frame of the planet, and the geometrical parameters of the system. 

In the text 
Fig. 2 Normalized Hough functions , , and corresponding to the Earth’s semidiurnal tide (m = 2 and ν = 1.0030). We note that ν = 2Ω /σ. Left: first gravity modes; 0 ≤ n ≤ 5. Right: first Rossby modes: − 6 ≤ n ≤ −1. In this case, ν is very close to 1 and the Rossby modes are trapped at the poles. The Earth’s semidiurnal tide is mainly described by gravity modes. For more details about the behaviour of Hough functions, see Fig. A.2. 

In the text 
Fig. 3 First gravity modes (0 ≤ n ≤ 5) of normalized Hough functions , , and corresponding to the semidiurnal tide of Venuslike exoplanets (m = 2 and ν = 0.4804). We note that ν = 2Ω /σ. For more details about the behaviour of Hough functions, see Figs. A.2. 

In the text 
Fig. 4 Frequency spectrum of the atmospheric tidal response. The identified regimes are given as functions of the tidal frequency (σ). 

In the text 
Fig. 5 Real and imaginary part of the vertical wavenumber as functions of the tidal (σ) and critical (σ_{c;n}) frequencies. Left: , with defined in Eq. (104), as a function of f^{(}σ/σ_{0}^{)} (horizontal axis) and (vertical axis), where f is the function defined by . The colour of the map, denoted c, is defined by . Hence, the purple region stands for slowly oscillating waves. The luminous (dark) regions correspond to strongly oscillating waves propagating downwards (upwards).Right: as a function of the same parameters. The luminous (dark) regions correspond to strongly (weakly) evanescent waves. The positions of the planets (Earth _{⊕}, Venus _{♀}) on the map are determined by the gravity mode of degree 0 in the case of the semidiurnal tide (Λ_{⊕;0} = 11.1596, Λ_{♀;0} = 7.1485). The tidal frequencies of the planets are given by σ_{⊕}/σ_{0} ≈ 194 and σ_{♀}/σ_{0} = 0.43, with σ_{0} = 7.5 × 10^{7} s^{1}. 

In the text 
Fig. 6 Relative difference between the complex wavenumber obtained in this work (Eq. (103)) on the one hand, and the real one given by CL70 for Earth atmospheric tides ( ) on the other hand. The difference η, given by Eq. (158), is plotted in logarithmic scale as a function of the normalized tidal frequency σ/σ_{0} and of the ratio h_{c}/h_{n}. The colour of the map, denoted c, is defined by the expression c = log (η). The horizontal and vertical axis represent f^{(}σ/σ_{0}^{)} and , where f is the function defined by with X ∈R. The positions of the planets are determined by the gravity mode of degree 0 in the case of the semidiurnal tide. Dark (luminous) regions correspond to the dynamic (dissipative) regime. The vertical wavenumbers of Venus are strongly affected by dissipative processes, while those of the Earth are not. 

In the text 
Fig. 7 Tidal torques exerted on the atmosphere in the case of the slowly rotating convective atmosphere (red line with the label “Conv”, Eq. (174)) and in the model by Correia & Laskar (2001) (green dashed line with the label “CL01”, Eq. (175)) as functions of the tidal frequency. Both functions are normalized by their respective maximum value, reached at σ = σ_{0}; the tidal frequency is normalized by the Newtonian cooling frequency (σ_{0}). The tidal forcing is quadrupolar, with U and J given by Sect. 4.3. The parameters a and b of Eq. (175) are adjusted so that the peaks of the two functions are superposed. 

In the text 
Fig. 8 Left: absorption of the flux emitted by the star along the starplanet direction. Right: absorption of the flux emitted by the ground. 

In the text 
Fig. 9 Vertical profiles of radiative fluxes (I), power per unit volume (J), and power per unit mass (J) for various values of the dimensionless optical depth τ_{λ} (Eq. (187)). As usual, the reduced altitude, x = z/H, is on the vertical axis. Left: incident flux, defined by Eq. (186), at ψ = 0. If τ_{λ} ≪ 1 the ground receives almost 100% of the total power. On the contrary, if τ_{λ} ≫ 1, all the power is absorbed by the atmosphere. Right: radiative flux emitted by the ground, defined by Eq. (206). For small τ_{λ}, the flux goes through the atmosphere without being absorbed. For high values of the parameter, it is absorbed by the lowest layers. 

In the text 
Fig. 10 Left, from top to bottom: power per unit mass (J^{inc}, Eq. (184)) and power per unit volume (J^{inc}, Eq. 183) of a monochromatic incident flux normalized by their maxima and plotted as functions of the normalized spatial coordinates X (the direction of the star is along the horizontal axis) and Z (the spin axis is along the vertical axis). The planet is represented by a green disk of unit radius. The flux, of wavelength λ, propagates in an atmosphere assumed to be homogeneous in composition (C_{G} = C_{0}) and characterized by the ratio H/R = 0.1. The optical depth is given by τ_{λ} = 0.2. Right, from top to bottom: power per mass unit (J^{rad}, Eq. (206)) and power per volume unit (J^{rad}) of a monochromatic flux radiated by the ground, normalized by their maxima, and plotted as functions of X and Z, with τ_{λ} = 0.5. 

In the text 
Fig. 11 Nyquist diagram of the boundary layer transfer function in the frequency range − 10^{3} ≤ σ/σ_{bl} ≤ 10^{3} (σ being the tidal frequency and σ_{bl} the frequency characterizing the boundary layer introduced in Eq. (217)). The parametrized complex transfer function (Eq. (215)) reduced by its value at σ = 0 is plotted as a function of the tidal frequency (σ) in the complex plane, where the horizontal axis corresponds to the real part of the function, , and the vertical axis to the imaginary part . The arrow indicates the direction of increasing σ. The gain of the transfert function is given by the distance of a point P of the curve to the origin O of coordinates , and its phase by the angle of the vector OP to the horizontal axis. This plot shows that the response of the ground is always delayed with respect to the thermal forcing except at synchronization, where and the gain is maximum. It also highlights the two identified regimes of the ground response: the lag tends to 0 (response in phase with the perturbation) and the gain increases at σ → 0, while the lag tends to ± π/ 4 and the gain to 0 at σ → ± ∞, the critical transition frequency being σ_{bl}. 

In the text 
Fig. 12 Bichromatic case. Top: heat transfers between the atmosphere, space and the ground. Bottom: power balance of the interface groundatmosphere. 

In the text 
Fig. 13 Average daily variation in temperature in degrees (top) and pressure in mbar (bottom) over one year, plotted over solar time (in hours), the subsolar point corresponding to 12 h. The data was recorded every minute over the full year 2013 with a Vantage Pro 2 weather station at latitude 48.363°N. The time is converted to solar time and the data averaged over an exact count of 365 days. The average values of temperature (10.8973° C) and pressure (1015.83 mbar) have been removed to display only the variable part. In both figures, the raw data (in green) has been decomposed in a Fourier series over the 24 h period, limited to two harmonics (blue), or five harmonics (red). Temperature is dominated by the diurnal component, with a maximum value around 14h26mn. The largest harmonics of the pressure variations is the semidiurnal term, followed by the diurnal term. Two pressure maxima can be observed, at 9h38mn and 22h07mn. 

In the text 
Fig. 14 Tidal torque applied to the atmosphere of Earth (left panel) and Venus (right panel) by the thermal semidiurnal tide and its components as functions of the reduced frequency χ = (Ω−n_{orb}) /n_{orb} (horizontal axis). The forcings consist of academic quadrupolar perturbations expressed as and with U_{2} = −0.985 m^{2} s^{2} and J_{2} = 1.0 × 10^{2} m^{2} s^{3} for the Earth and U_{2} = −2.349 m^{2} s^{2} and J_{2} = 1.0 × 10^{4} m^{2} s^{3} for Venus. Labels T, H, and V refer to the total response and its horizontal and vertical components, respectively (we note that ). The label CL01 refers to the model by Correia & Laskar (2001), given by Eq. (175). 

In the text 
Fig. 15 Total tidal response of the Earth’s atmosphere to the solar semidiurnal perturbation caused by the quadrupolar academic forcings and with U_{2} = −0.985 m^{2} s^{2} and J_{2} = 1.0 × 10^{2} m^{2} s^{3}. Left, from top to bottom: amplitudes of δp, δρ, δT, V_{θ}, V_{ϕ}, and V_{r} in logarithmic scales as functions of the reduced latitude δ/π (horizontal axis) and altitude x = z/H in logarithmic scale (vertical axis). The colour function c is given by c = log (δf) for any quantity δf. Right, from top to bottom: Sine of the arguments of the same quantities as functions of the reduced latitude and altitude. In these plots, c = sin[arg(δf)]. 

In the text 
Fig. 16 Total tidal response of the Earth’s atmosphere to the solar semidiurnal perturbation caused by the quadrupolar academic forcings and with U_{2} = −0.985 m^{2} s^{2} and J_{2} = 1.0 × 10^{2} m^{2} s^{3}. Left, from top to bottom: real parts of δp, δρ, δT, V_{θ}, and V_{ϕ} taken at the ground as functions of the reduced longitude ϕ/π (horizontal axis) and latitude δ/π (vertical axis). The colour function c is given by c = ℜ{δf} for any quantity δf. Right, from top to bottom: imaginary parts of the same quantities as functions of the reduced longitude and latitude. In these plots, c = ℑ{δf}. 

In the text 
Fig. 17 Tidal surface oscillations of pressure, density, temperature, and their components at the equator of the Earth as functions of the reduced tidal frequency χ = (Ω−n_{orb}) /n_{orb} (horizontal axis). The perturbation is defined by the academic quadrupolar forcings and with U_{2} = −0.985 m^{2} s^{2} and J_{2} = 1.0 × 10^{2} m^{2} s^{3}. Labels T, H, and V refer to the total response and its horizontal and vertical components, respectively. 

In the text 
Fig. 18 Total tidal response of the atmosphere of Venus to the solar semidiurnal perturbation caused by the quadrupolar academic forcings and with U_{2} = −2.349 m^{2} s^{2} and J_{2} = 1.0 × 10^{4} m^{2} s^{3}. Left, from top to bottom: amplitudes of δp, δρ, δT, V_{θ}, V_{ϕ}, and V_{r} in logarithmic scales as functions of the reduced latitude δ/π (horizontal axis) and altitude x = z/H in logarithmic scale (vertical axis). The colour function c is given by c = log (δf) for any quantity δf. Right, from top to bottom: Sine of the arguments of the same quantities as functions of the reduced latitude and altitude. In these plots, c = sin[arg(δf)]. 

In the text 
Fig. 19 Total tidal response of the atmosphere of Venus to the solar semidiurnal perturbation caused by the quadrupolar academic forcings and with U_{2} = −2.349 m^{2} s^{2} and J_{2} = 1.0 × 10^{4} m^{2} s^{3}. Left, from top to bottom: real parts of δp, δρ, δT, V_{θ}, and V_{ϕ} taken at the ground as functions of the reduced longitude ϕ/π (horizontal axis) and latitude δ/π (vertical axis). The colour function c is given by c = ℜ{δf} for any quantity δf. Right, from top to bottom: imaginary parts of the same quantities as functions of the reduced longitude and latitude. In these plots, c = ℑ{δf}. 

In the text 
Fig. 20 Tidal surface oscillations of pressure, density, temperature, and their components at the equator of Venus as functions of the reduced tidal frequency χ = (Ω−n_{orb}) /n_{orb} (horizontal axis). The perturbation is defined by the academic quadrupolar forcing and with U_{2} = −2.349 m^{2} s^{2} and J_{2} = 1.0 × 10^{4} m^{2} s^{3}. Labels T, H, and V refer to the total response and its horizontal and vertical components, respectively. 

In the text 
Fig. 21 Asymptotic domains of parameters where the analytical approach presented in this work can be applied (green) and where the numerical approach is necessary (red) as functions of the frequency ratios N^{2}/σ^{2} (horizontal axis) and (2Ω)^{2}/σ^{2} (vertical axis). Green regions: stably stratified atmospheres satisfying the hierarchy of frequencies (2Ω)^{2} ≪ N^{2} necessary to apply the traditional approximation a); weakly stratified or convective atmospheres of slowly rotating planets b). Red regions: stably stratified atmospheres where the hierarchy of frequencies (2Ω)^{2} ≪ N^{2} is not satisfied c); superadiabatic rotating convective atmospheres d). 

In the text 
Fig. 22 Predictive aspects related to the dependence of the atmospheric tidal torque on stratification. The neutral convective case implies a strong torque, which can lead to nonsynchronized states of equilibrium. The stably stratified case implies a weak torque, which lets the planet evolve towards spinorbit synchronization. 

In the text 
Fig. A.1 Top: eigenvalues () of the gravity and Rossby modes of lowest degrees (− 6 ≤ n ≤ 5) as functions of the spin parameter (ν). The horizontal axis represents the function of ν given by , and the vertical axis the function of given by . Bottom: projection coefficients of the normalized Legendre polynomial on normalized Hough functions associated with these modes as functions of the spin parameter (see Eq. (70)). 

In the text 
Fig. A.2 Normalized Hough functions computed as functions of the colatitude (θ) for 0 ≤ n ≤ 5 (gravity modes) and various values of ν = (2Ω) /σ. Topleft: ν = 0; topright: ν = 0.8; middleleft: ν = 1.2; middleright: ν = 2; bottomleft: ν = 4; bottomright: ν = 12. The case ν = 0 corresponds to the normalized associated Legendre polynomials . The two plots at the top represent the regime of superinertial waves (ν ≤ 1) where Hough functions look like Legendre polynomials. In the regime of subinertial waves (ν > 1), functions become more and more evanescent when ν increases, and oscillations are trapped around the equator in the interval [θ_{ν},π−θ_{ν}], where . 

In the text 
Fig. A.3 Normalized Hough functions computed as functions of the colatitude (θ) for − 6 ≤ n ≤ −1 (Rossby modes) and various values of ν = (2Ω) /σ. Topleft: ν = 1.2; topright: ν = 2; bottomleft: ν = 4; bottomright: ν = 12. In the regime of subinertial waves (ν > 1), Rossby modes are trapped around the poles when ν → 1. 

In the text 
Fig. C.1 Inertial reference, orbital and equatorial rotating frames, and the associated Euler’s angles of orientation. Figure taken from Mathis & Le PoncinLafitte (2009). 

In the text 
Fig. E.1 Tidal torques exerted on the atmosphere (in petaNewton metres (PN.m)) induced by the first gravity mode of the Earth’s semidiurnal tide (n = 0) as a function of the reduced tidal frequency (χ) obtained with two different conditions at the upper boundary: the energybounded condition chosen for the model (red line) and a stressfree condition (green dashed line). The equation for the structure in the vertical (resp. horizontal) direction has been integrated with 10 000 mesh points (resp. 400 associated Legendre polynomials). 

In the text 
Fig. F.1 Tidal response of the Earth’s atmosphere to a quadrupolar forcing ( and ) as a function of the reduced tidal frequency χ = (Ω−n_{orb})/n_{orb} (horizontal axis). For each quantity q_{n} associated with the Hough mode of degree n, except the coefficients (top left) and the scales of variations L_{V;n}, the function f(χ) = sign(q_{n})log (1 + q_{n}) is plotted (vertical axis) for n ∈ ˜−6,5¨. Left, from top to bottom: a) coefficients (Eq. (70)); b) eigenvalues of Laplace’s tidal equation ( defined by Eq. (35)); c) normalized variation scales of the vertical structure of the vertical component, L_{V;n} = 2π/ (Eq. (48)), plotted in logarithmic scale. Right, from top to bottom: d) real part of the vertical wavelengths ( , Eq. (103)); e) imaginary part of the vertical wavelengths; f) spin parameter ν = 2Ω /σ (Eq. (23)). 

In the text 
Fig. G.1 Real part of the solution of Eq. (G.1) as a function of x in logarithmic scale computed numerically for various mesh resolutions. The number of points of the spatial discretization is set at N = 20 (green dashed line), N = 200 (blue dashed line), N = 2000 (pink dashed line), and N = 20 000 (cyan dashed line). These numerical solutions are compared to the analytical solution given by Eq. (129) (red line). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.