A generic frequency dependence for the atmospheric tidal torque of terrestrial planets

Thermal atmospheric tides have a strong impact on the rotation of terrestrial planets. They can lock these planets into an asynchronous rotation state of equilibrium. We aim at characterizing the dependence of the tidal torque resulting from the semidiurnal thermal tide on the tidal frequency, the planet orbital radius, and the atmospheric surface pressure. The tidal torque is computed from full 3D simulations of the atmospheric climate and mean flows using a generic version of the LMDZ general circulation model (GCM) in the case of a nitrogen-dominated atmosphere. Numerical results are discussed with the help of an updated linear analytical framework. Power scaling laws governing the evolution of the torque with the planet orbital radius and surface pressure are derived. The tidal torque exhibits i) a thermal peak in the vicinity of synchronization, ii) a resonant peak associated with the excitation of the Lamb mode in the high frequency range, and iii) well defined frequency slopes outside these resonances. These features are well explained by our linear theory. Whatever the star-planet distance and surface pressure, the torque frequency spectrum -- when rescaled with the relevant power laws -- always presents the same behaviour. This allows us to provide a single and easily usable empirical formula describing the atmospheric tidal torque over the whole parameter space. With such a formula, the effect of the atmospheric tidal torque can be implemented in evolutionary models of the rotational dynamics of a planet in a computationally efficient, and yet relatively accurate way.


Introduction
Understanding the evolution of planetary systems has become a crucial question with the rapidly growing number of exoplanets discovered up to now. Terrestrial planets particularly retain our attention as they offer a fascinating diversity of orbital configurations, and possible climates and surface conditions. This diversity is well illustrated by Proxima-b, an exo-Earth with a minimum mass of 1.3 M ⊕ orbiting Proxima Centauri (Anglada-Escudé et al. 2016;Ribas et al. 2016), and the TRAPPIST-1 system, which is a tightly-packed system of seven Earth-sized planets orbiting an ultracool dwarf star (Gillon et al. 2017;Grimm et al. 2018).
Characterizing the atmospheric dynamics and climate of these planets is a topic that motivated numerous theoretical works, both analytical and numerical (e.g. Pierrehumbert 2011;Heng & Kopparla 2012;Leconte et al. 2013;Heng & Workman 2014;Wolf et al. 2017;Wolf 2017;Turbet et al. 2018). This tendency will be reinforced in the future by the rise of forthcoming space observatories such as the James Webb Space Telescope (JWST), which will unravel features of the planetary atmospheric structure by performing high resolution spectroscopy over the infrared frequency range (Lagage 2015).
Constraining the climate and surface conditions of the observed terrestrial planets requires constraint of their rotation rate first because of the key role played by this parameter in the equilibrium atmospheric dynamics (Vallis 2006;Pierrehumbert 2010). Particularly, it is important to know whether a planet is locked into the configuration of spin-orbit synchronization with its host star and the extent to which asynchronous rotation states of equilibrium might exist. Over long timescales, the planet rotation is driven by tidal effects, that is the distortion of the planet by its neighbours (star, planets and satellites) resulting from mutual distance interactions. Tides are a source of internal dissipation inducing a variation of mass distribution delayed with respect to the direction of the perturber. As a consequence, the planet undergoes a tidal torque, which modifies its rotation by establishing a transfer of angular momentum between the orbital and spin motions.
Tides can be generated by forcings of different natures. First, the whole planet is distorted by the gravitational tidal potential generated by the perturber, and is driven by the resulting tidal torque towards spin-orbit synchronous rotation and a circular orbital configuration. Second, if the perturber is the host star, the atmosphere of the planet undergoes a heating generated by the day-night cycle of the incoming stellar flux. The variations of the atmospheric mass distribution generated by this forcing are the so-called thermal atmospheric tides (Chapman & Lindzen 1970).
As demonstrated by the pioneering study by Gold & Soter (1969) in the case of Venus, thermal tides are able to drive a terrestrial planet away from spin-orbit synchronization since they induce a tidal torque in opposition with that resulting from solid tides in the low frequency range. Hence, the competition between the two effects locks the planet into an asynchronous rotation state of equilibrium, which explains the departure of the rotation rate of Venus to spin-orbit synchronization.
The understanding of this mechanism has been progressively consolidated by analytical works based upon the classical tidal theory (e.g. Ingersoll & Dobrovolskis 1978;Dobrovolskis & Ingersoll 1980;Auclair-Desrotour et al. 2017a,b) or using parametrized models (Correia & Laskar 2001. Over the past decade, the growing performances of computers have made full numerical approaches affordable, and the atmospheric torque created by the thermal tide was computed using general circulation models (GCM; Leconte et al. 2015). This approach remains complementary with analytical models owing to its high computational cost. However, it is particularly interesting since it allows to characterize the atmospheric tidal response of a planet by taking into account the atmospheric structure, mean flows and other internal processes by solving the primitive equations of fluid dynamics in a self-consistent way.
By using a generic version of the LMDZ GCM (Hourdin et al. 2006), Leconte et al. (2015) retrieved the frequencydependence of the tidal torque predicted by ab initio analytical models (Ingersoll & Dobrovolskis 1978;Auclair-Desrotour et al. 2017a. The torque increases linearly with the tidal frequency in the vicinity of synchronization. It reaches a maximum associated with a thermal time of the atmosphere and then decays in the high-frequency range. This behaviour is approximated at the first order by the Maxwell model, which describes the forced response of a damped harmonic oscillator. It shows evidence of the important role played by dissipative processes such as radiative cooling in Venus-like configurations. To better understand the action of the thermal tide on the planet rotation, this frequency-dependent behaviour has to be characterized. Thus, our purpose in this study is to investigate the dependences of the tidal torque created by the semidiurnal tide on the tidal frequency and on key control parameters. We follow along the line by Leconte et al. (2015) for the method, and treat the case of an idealized dry terrestrial planet hosting a nitrogendominated atmosphere and orbiting a Sun-like star. Hence, we recall in Sect. 2 the mechanism of the thermal atmospheric tide. In Sect. 3, we detail the method and the physical setup of the treated case.
In Sect. 4, we compute the tidal torque exerted on the atmosphere from simulations using the LMDZ GCM and examine its dependence on the tidal frequency. We introduce in this section two new models for the thermally generated atmospheric tidal torque: an ab initio analytical model based upon the linear theory of atmospheric tides (e.g. Chapman & Lindzen 1970), and a parametrized semi-analytical model derived from results obtained using GCM simulations. This later model describes in a realistic way the behaviour of the torque in the low-frequency range, where a thermal peak is observed. In addition, we investigate in this section the role played by the ground-atmosphere thermal coupling in the lag of the tidal bulge.
In Sect. 5, we examine the dependence of the tidal response on the planet orbital radius and surface pressure. We thus establish empirical scaling laws describing the evolution of the characteristic amplitude and timescale of the thermal peak with these two parameters. Combining together the obtained results, we finally derive a new generic formula to quantify the atmospheric tidal torque created by the thermal semidiurnal tide in the case of a N 2 -dominated atmosphere. We give our conclusions in Sect. 6.

Basic principle
We briefly recall in this section the main aspects of the mechanism of atmospheric tides in the case of terrestrial planets, and we introduce analytical expressions that will be used in the following to compute the resulting tidal torque. For the sake of simplification, we consider in this study the case of a spherical planet of radius R p and mass M p , orbiting its host star, of mass M , circularly. The star-planet distance is denoted a, the mean motion of the system n , and the obliquity of the planet is set to zero. We assume that the planet rotates at the spin angular velocity Ω, which is positive if the spin rotation is along the same direction as the orbital motion, and negative otherwise.
The atmosphere of the planet undergoes both the tidal gravitational and thermal forcings of the host star. Below a certain orbital radius, the planet is sufficiently close to the star to make gravitational forces predominate. Thus, its rotation is driven towards spin-orbit synchronization (Ω = n ), which is the unique possible final state of equilibrium for the planet rotation in the absence of obliquity and eccentricity. Conversely, the predominance of the thermal tide enables the existence of asynchronous final rotation states of equilibrium, as showed in the case of Venus (e.g. Gold & Soter 1969;Ingersoll & Dobrovolskis 1978;Dobrovolskis & Ingersoll 1980;Correia & Laskar 2001;Auclair-Desrotour et al. 2017a). As a consequence, we ignore here the action of gravitational forces on the atmosphere. We note however that the action of these forces on the atmospheric tidal bulge will be taken into account to compute the tidal torque, as seen in the following.
The thermal forcing results from the day-night periodic cycle. The atmosphere undergoes heating variations due to the time-varying component of the incoming stellar flux F, which scales as the equilibrium one, F = L / 4πa 2 , where L is the luminosity of the star. Hence, the absorbed energy induces a delayed variation of the atmospheric mass distribution. Let us assume the hydrostatic approximation (i.e. that pressure and gravitational forces compensate each other exactly in the vertical direction) and consider that the surface of the planet is rigid enough to support the atmospheric pressure variations with negligible distortions. It follows that the variation of mass distribution is directly proportional to the surface pressure anomaly, where l and m designate the latitudinal and longitudinal degrees of a mode, θ and ϕ the colatitude and longitude in the reference frame co-rotating with the planet, t the time, Y m l the normalized spherical harmonics (see Appendix A), δp m,σ s;l the associated components, and σ = m (Ω − n ) the associated forcing frequencies (see e.g. Efroimsky 2012; Ogilvie 2014) 1 .
The tidal torque exerted on the atmosphere is obtained by integrating the gravitational force undergone by the tidal bulge over the sphere. Hence, denoting U T the tidal gravitational potential at the planet surface, the atmospheric tidal torque is defined in the thin-layer approximation (H R p ) as (e.g. Zahn 1966) where the notation g refers to the surface gravity of the planet, ∂ ϕ to the partial derivative in longitude, S to the sphere of radius R p , and dS = R 2 p sin θdθdϕ to the surface element.
Similarly as the surface pressure anomaly, U T can be expanded in Fourier series of time and spherical harmonics, where the U m,σ T;l are the amplitudes of the different modes. Terms associated with l = 1 do not contribute to the tidal torque since they just correspond to a displacement of the planet gravity centre. Thus, the main components of the expansion are those associated with the quadrupolar semidiurnal tide, that is with degrees l = |m| = 2. Besides, since the U m,σ T;l scale as U m,σ T;l ∝ R p /a l , terms of greater order in l can be neglected with respect to the quadrupolar components if the radius of the planet R p is assumed to be small compared to the star-planet distance, which is the case in the present study. Thus, by substituting U T and δp s by their expansions in spherical harmonics in Eq. (2), we note that the quadrupolar terms l = |m| = 2 only remain, as for the tidal potential U T , and we end up with the well-known expression of the semidiurnal quadrupolar torque in the thin layer approximation (e.g. Leconte et al. 2015), with σ = 2 (Ω − n ), and where the notation refers to the imaginary part of a complex number ( referring to the real part). In this expression, δp 2,σ s;2 designates the component of degrees l = 2 and m = 2 in the expansion on spherical harmonics given by Eq. (1). This complex quantity is the most important one since it encompasses the whole physics of the atmospheric tidal response. In the following, it will be calculated using a GCM.
The action of the torque on the planet is fully determined by the sign of the product η = sign (σ) δp 2,σ s;2 . When η < 0 (η > 0), the atmospheric tidal torque pushes the planet towards (away from) spin-orbit synchronization, |Ω − n | decays (increases). Positions for which η = 0 correspond to the stable ( dη/dσ| eq < 0) or unstable ( dη/dσ| eq > 0) equilibrium rotation rates that the planet would reach if it were subject to atmospheric tides only, that is if solid tides were ignored in the case of a dry terrestrial planet.

Method
As mentioned in the previous section, the guideline of the method is to compute the quadrupolar component of the surface pressure anomaly from 3D GCM simulations. We detail here the basic physical setup of these simulations in a first time, and the way δp 2,σ s;2 is extracted from pressure snapshots in a second time. In the whole study, we focus on a Venus-sized planet orbiting a Sun-like star.
A 'reference case' of fixed surface pressure and star-planet distance is defined. Specifically, the surface pressure is set in this case to p s = 10 bar and we assume that the planet is located at the Venus-Sun distance, that is a Venus = 0.723 au. This configuration, characterized in Sect. 4, corresponds to the case illustrated by Fig. 1 of Leconte et al. (2015), and seems thereby a convenient choice for comparisons with this early work.
In Sect. 5, two families of configurations will be studied, both including the reference case. In the first family, the surface pressure is set to p s = 10 bar and the semi-major axis varies.
Conversely, in the second family, planets have the same orbital radius, a = a Venus , and various surface pressures.

Physical setup of the 3D simulations
Apart from the surface pressure and star-planet distance, all simulations are based on a common physical setup. For the stellar incoming flux, the emission spectrum of the Sun is used. The planet is assumed to be dry, with no surface liquid water or water vapour, which allows us to filter out effects associated with the formation of clouds in the study of its atmospheric tidal response. The atmosphere is arbitrarily assumed to be nitrogendominated. However, a pure N 2 -atmosphere would be an extreme case for radiative transfer owing to the absence of radiator. Hence, we have to set a non-zero volume mixing ratio for carbon dioxyde to avoid numerical issues in the treatment of the radiative transfer with the LMDZ, which was originally designed to study the Earth atmosphere. Although any value could be used, we choose to set the value of the CO 2 volume mixing ratio to that of the Earth atmosphere at the beginning of the twenty-first century, that is ∼370 ppmv (e.g. Etheridge et al. 1996). The mass ratio corresponding to this volume mixing ratio being negligible, we use the value of N 2 for the mean molecular mass of the atmosphere, M atm = 28.0134 g mol −1 (Meija et al. 2016).
For a perfect diatomic gas, the ratio of heat capacities (also called first adiabatic exponent) is Γ 1 = 1.4, and it follows that κ = (Γ 1 − 1) /Γ 1 = 0.285 (the parameter κ can also be written κ = R GP / M atm C p , where R GP and C p stand for the perfect gaz constant and the thermal capacity per unit mass of the atmosphere, respectively). The effects of topography are ignored and the surface of the planet is thus considered as an isotropic sphere of albedo A s = 0.2 and thermal inertia I gr = 2000 J m −2 s −1/2 K −1 , which is a typical value for Venus-like soils (see e.g. Lebonnois et al. 2010) 2 . All of these parameters remain unchanged for the whole study and are summarized in Table 1. Our simulations are performed with an upgraded version of the LMD GCM specifically developed for the study of extrasolar planets and paleoclimates (see e.g. Wordsworth et al. 2010Wordsworth et al. , 2011Wordsworth et al. , 2013Forget et al. 2013;Leconte et al. 2013), and used previously by Leconte et al. (2015) for the study of atmospheric tides. The model is based on the dynamical core of the LMDZ 4 GCM (Hourdin et al. 2006), which uses a finite-difference formulation of the primitive equations of geophysical fluid dynamics. Particularly, the following approximations are assumed.
The main one is the hydrostatic approximation (e.g. Vallis 2006), meaning that the pressure and gravitational forces compensate each other exactly along the vertical direction. The second approximation is the traditional approximation (e.g. Unno et al. 1989), which consists in ignoring the components of the Coriolis acceleration associated with a vertical motion of fluid particles or generating a force along the vertical direction. The third important assumption in the code is the thin layer approximation, meaning that the thickness of the atmosphere is considered as small with respect to the radius of the planet (e.g. Vallis 2006).
A spatial resolution of 32 × 32 × 26 in longitude, latitude, and altitude is used for the simulations.
The radiative transfer is computed in the model using a method similar to Wordsworth et al. (2011) and Leconte et al. (2013). High-resolution spectra characterizing optical properties were preliminary produced for the chosen gas mixture over a A&A 624, A17 (2019) wide range of temperatures and pressures using the HITRAN 2008 database (Rothman et al. 2009). These spectra are interpolated every radiative timestep during simulations to determine local radiative transfers. The method is commonly used and has been thoroughly discussed in past studies (e.g. Leconte et al. 2013). We thus refer the readers to these works for a detailed description.

Extraction of the quadrupolar surface pressure anomaly
For a given planet, of fixed rotation, semi-major axis and surface pressure, the calculation of the quadrupolar surface pressure anomaly follows several steps.
First, the GCM is run for a period P conv corresponding to the convergence timescale necessary to reach a steady cycle. We note that this period has to be specified for each doublet (a, p s ). As a first approximation, it depends on the radiative timescale of the deepest layers of the atmosphere τ rad , which scales as where T e stands for the mean effective, or black body, temperature of the atmosphere (see e.g. Showman & Guillot 2002, Eq. (10)). In the reference case, we observe that the atmospheric state has converged towards a steady cycle after ∼5800 Earth Solar days, and we thus use this value for calculations in this section. After this first step, a simulation is run for 300 Solar days of the planet, defined by P sol = 2π/ |Ω − n |, except in the case of the spin-orbit synchronization (Ω = n ), where there is no daynight cycle (in this case, the simulation is simply run for 3000 Earth Solar days). At the end of the simulation, we have at our disposal a time series of snapshots of the surface pressure given as a function of the longitude and latitude (see e.g. Fig. 1).
The third step consists in post-processing these data. We first remove the constant component, that is the mean surface Surface pressure and horizontal winds computed with the LMDZ GCM for a Venus-sized terrestrial planet hosting a 10 bar atmosphere (reference case). In this study, the surface pressure anomaly is folded over one Solar day and expanded in spherical harmonics to calculate the atmospheric tidal torque using the formula given by Eq. (4).
pressure. Then, we proceed to a change of variable: the time coordinate is replaced by the Solar zenith angle, so that snapshots are all centred on the substellar point. Since meteorological fluctuations can be considered as a perturbation varying randomly over short timescales, we filter them out by folding the surface pressure anomaly over one Solar day.
We finally apply a spherical harmonics transform to the resulting averaged surface pressure snapshot in order to get the complex coefficient δp 2,σ s;2 associated with the semidiurnal tidal mode of degrees l = 2 and m = 2 (see Eq. (1)). The method is illustrated by Fig. 2 in the reference case (a = a Venus and p s = 10 bar).
This procedure provides the value of the tidal torque for a given forcing frequency. In practice, the torque is computed over an interval of the normalized frequency ω = (Ω − n ) /n centred on synchronization (ω = 0) with n fixed, the planet rotation rate being deduced from ω (the normalized frequency ω is employed here instead of σ to follow along the line by Leconte et al. 2015). Typically, we use −30 ≤ ω ≤ 30 to study the low-frequency regime of the atmospheric tidal response and −300 ≤ ω ≤ 300 to study the high-frequency regime.
The frequency range is thus divided into N intervals, meaning that the whole above procedure has to be repeated N + 1 times to construct a frequency-spectrum of the tidal torque. The size of an interval is defined as ∆ω ≡ ω sup − ω inf /N. For instance, for the exploration of the parameters space detailed in Sect. 5, N = 20, ω inf = − 30, ω sup = 30, and thus ∆ω = 3.

Frequency behaviour of the atmospheric tidal torque
The apparent complexity of the physics involved in thermal atmospheric tides requires that we opt for a graduated approach of the problem. Hence, before investigating the dependence of the tidal torque on the planet orbital radius and atmospheric surface pressure as mentioned above, we have to preliminary characterize how it varies with the tidal frequency. To address Mean pressure anomaly [Pa] Semidiurnal component [Pa] Fig. 2. Surface pressure anomaly created by the thermal tide. Left panels: daily averaged spatial distribution of the departure of the surface pressure from its mean value created by the thermal tide. Right panels: spatial distribution of the semidiurnal component only. The surface pressure anomaly is computed for 300 Solar days and folded over one Solar day centred on the substellar point, whose location and direction of motion are shown with a white arrow. From top to bottom panels: normalized forcing frequency ω = (Ω − n ) /n is increased from 0 (spin-orbit synchronization) to 24 (this corresponds to the length of the Solar day P sol = 9.36 days) for the reference case of the study (a = a Venus and p s = 10 bar).
this question, we consider the reference case (p s = 10 bar and a = a Venus ).

Characterization of the reference case
In order to characterize the reference case, frequency-spectra of the atmospheric torque created by the semidiurnal thermal tide are computed in low-frequency and high-frequency ranges. For convenience, we introduce the function f GCM (σ), which is the interpolating function of GCM results with cubic splines.
Noting that the tidal torque should be an odd function of the tidal frequency in the absence of rotation (or if the effect of rotation on the tidal response were negligible), we also introduce the function f odd , defined by which is the odd function f minimizing for any σ the distance defined by The complementary function f even , such that f GCM = f odd + f even , is defined by " 10 bar). Top left: and provides a measure of the impact of Coriolis effects on the tidal torque. The data, the interpolating function f GCM and its components f odd and f even are plotted in Fig. 3 as functions of the normalized tidal frequency ω = σ/ (2n ) in linear and logarithmic scales. Additional functions of the frequency are plotted in dashed lines. They correspond to the ab initio analytical ('Ana.'), Maxwell, and parametrized ('Param.') models that will be introduced and discussed further.
We first consider the low-frequency range (−30 ≤ ω ≤ 30). The reference case of our study exactly reproduces the results plotted in Fig. 1 of Leconte et al. (2015), with a maximum slightly greater that 2000 Pa located around ω ∼ 5. We introduce here the maximal value of the peak q max ≡ max { f odd (σ)} and the associated frequency σ max , such that f odd (σ max ) = q max , timescale τ max ≡ σ −1 max and normalized fre- The tidal torque is negative for σ < 0 and positive otherwise, which corresponds to the typical behaviour of the thermally induced atmospheric tidal response in the vicinity of synchronization, as discussed in Sect. 2. As shown by early studies (Gold & Soter 1969;Ingersoll & Dobrovolskis 1978;Dobrovolskis & Ingersoll 1980;Correia & Laskar 2001), thermal atmospheric tides thus tends to drive the planet away from synchronous rotation and determines its non-synchronized rotation states of equilibrium.
In the zero-frequency limit, the torque scales as T ∝ σ α , with α ≈ 0.73. In the high-frequency range (20 |ω| ≤ 300), it scales as T ∝ σ −1 with a remarkable regularity (see Fig. 3, bottom left panel) and exhibits a resonance at ω ≈ 260. We will see in the next section that these features can be explained using the linear theory of atmospheric tides (Wilkes 1949;Siebert 1961;Lindzen & Chapman 1969).
We note that the spectrum of f GCM exhibits a slight systematic asymmetry with respect to the synchronization. This feature is obvious in the low-frequency range, where | f GCM (−σ)| > | f GCM (σ)|, and tends to vanish while |σ| increases. Particularly, a small departure between f GCM and f odd can be observed around the extrema of the tidal torque, and we note that the atmosphere undergoes a non-negligible tidal torque at synchronization (σ = 0), although the perturber does not move in the reference frame co-rotating with the planet.
This asymmetry is an effect of the Coriolis acceleration, which comes from the fact that |Ω (−σ)| |Ω (σ)| (in the low-frequency range, the spin rotation rate is not proportional to the tidal frequency). The Coriolis acceleration affects the atmospheric general circulation by generating strong zonal jets through the mechanism of non-linear Rossby waves pumping angular momentum equatorward (e.g. Showman & Polvani 2011). These jets induce a Doppler-like angular lag of the tidal bulge with respect to the direction of the perturber.

Ab initio analytical model
The behaviour of the torque in the high-frequency range can be explained with the help of the linear theory of thermal atmospheric tides (Wilkes 1949;Siebert 1961;Lindzen & Chapman 1969). In Appendix B, by using an ab initio approach, we compute analytically the atmospheric tidal torque created by the semidiurnal thermal tide in the idealized case of an isothermal atmosphere undergoing the tidal heating of the planet surface. The atmospheric structure is here characterized by the constant pressure height where R s and T s designate the specific gas constant and the surface temperature, respectively. It allows to renormalizes the altitude z with the introduction of the pressure height scale x = z/H. In the analytic model, we choose for the heat per unit mass inducing the tidal response the vertical profile J = J s e −b J x , where J s is the heat per unit mass at the planet surface, and b J a dimensionless optical depth corresponding to the inverse of the characteristic thickness of the heated layer. We note that the limit b J → +∞ corresponds to the case studied by Dobrovolskis & Ingersoll (1980), where the vertical profile of heat is approximated by a Dirac distribution. The surface pressure anomaly is obtained by solving the vertical structure equation of the dominating mode with the above profile of the forcing. We refer the reader to the appendix for the detail of approximations and calculations made to get this result. Particularly, we note that dissipative processes are ignored since they are associated with timescales that are supposed to far exceed typical tidal periods in the high-frequency range.
The solution takes two different forms depending on the way σ compares to the frequency characterizing the turning point, where the vertical wavenumber annihilates (see Appendix B), The notation Λ 0 designates here the eigenvalue of the predominating mode in the expansion of perturbed quantities on the basis of Hough functions (see Eq. (B.17)). This mode is the gravity mode of latitudinal wavenumber n = 0 in the indexing notation used by Lee & Saio (1997). Its eigenvalue Λ 0 can be approximated as a constant provided that n |Ω|. Hence, introducing the equivalent depth of the mode, we obtain, for |σ| ≤ σ TP , and, for |σ| > σ TP , We recall that κ = R s /C p , where C p designates the heat capacity per unit mass, and Γ 1 = 1/ (1 − κ) the adiabatic exponent at constant entropy (Gerkema & Zimmerman 2008). The solution given by Eqs. (12) and (13) provides a useful diagnosis about the frequency-behaviour of the torque in the high-frequency range.
The most striking feature of this behaviour is the peak that can be observed in Fig. 3 (top and bottom left panels) at the normalized frequency ω ≈ 260. This peak correspond to the fundamental resonance of the atmospheric vertical structure associated with the propagation of the Lamb mode (e.g. Lindzen et al. 1968;Bretherton 1969;Lindzen & Blake 1972;Platzman 1988;Unno et al. 1989), which is an acoustic type wave of long horizontal wavelengh. In an inviscid, isothermal atmosphere, the Lamb mode is characterized by the equivalent depth h L = Γ 1 H (Lindzen & Blake 1972). In the asymptotic regime, where n |Ω|, the characteristic Lamb frequency follows from Eq. (11), By noticing that σ L > σ TP in the case of a diatomic gas (Γ 1 = 1.4) and substituting h by h L into the corresponding expression of the solution -that is Eq. (13) -we can easily observe that the tidal torque is singular at |σ| = σ L . The resonance hence occurs when the phase velocity of the forced mode equalizes the characteristic Lamb velocity V L = gh L .
With the numerical values given by Table 1 and the mean surface temperature computed from GCM simulations (T s ≈ 316 K), the isothermal approximation leads to H ≈ 10.6 km and h L ≈ 15 km for the reference case. Besides, Λ 0 ≈ 11.1 in the adiabatic asymptotic regime of high rotation rates. It thus follows that ω L ≈ 308, and we recover the order of magnitude of the frequency identified in Fig. 3 (top left panel) using GCM simulations (i.e. ω L ≈ 260).
The observed departure between values of ω L obtained in analytical and numerical approaches can be explained by the dependence of the resonance on the atmospheric vertical structure (see e.g. Bretherton 1969;Lindzen & Blake 1972). The analytical value corresponds to the case of an isothermal atmosphere of temperature T s . In reality, the mean temperature vertical profile is characterized by a strong gradient in the troposphere, the temperature decaying linearly from ∼316 K at z = 0 to ∼160 K at z ≈ 25 km in GCM simulations. As a consequence, the mean pressure height scale of the tidally heated layer is less than the surface pressure heigh scale, which leads to smaller equivalent depth and resonance frequency for the Lamb mode.
The other interesting feature highlighted by Fig. 3 is the scaling law of the torque T ∝ σ −1 in the range of intermediate frequencies, that is between the thermal and Lamb resonances, typically. This behaviour is described by the analytical model. As discussed before (see Eq. (14)), σ TP and σ L are close to each other. The intermediate-frequency range thus corresponds to the case |σ| < σ TP , which leads us to consider the solution given by Eq. (12). We place ourselves in the configuration where characteristic timescales are clearly separated, that is |σ| σ TP and n |Ω| in the meantime. As H/h ∝ σ −2 , the preceding condition implies that H/h 1. It follows that By invoking the strong optical thickness of the atmosphere in the infrared (b J 1), we remark that we recover analytically the scaling law T ∝ σ −1 observed in Fig. 3 from the moment that the condition 1 H/h κ −1 b 2 J is satisfied. This provides a definition for the intermediate frequency-range, which is now the range corresponding to σ J |σ| σ L , where we have introduced the thermal frequency Basically, σ J is the frequency for which the vertical wavelength of the mode and the characteristic depth of the heated layer are of the same order of magnitude.
From the moment that |σ| σ L (or H/h 1), Eq. (12) can be approximated by the function where the associated characteristic timescale τ J and maximal amplitude of the pressure anomaly q J are We recognize in the form of the function given by Eq. (17) the well-known Maxwell model, which is commonly used to describe the dependence of the tidally dissipated energy on the forcing frequency in the case of solid bodies (e.g. Efroimsky 2012; Correia et al. 2014). Its use in the case of thermal atmospheric tides is discussed in the next section.

Discussion on the Maxwell model
Analytic ab initio approaches based on a linear analysis of the atmospheric tidal response -including this work (cf. previous section) -predict that the imaginary part of surface pressure variations can be expressed as a function of the forcing frequency σ = 2 (Ω − n ) as (e.g. Ingersoll & Dobrovolskis 1978;Auclair-Desrotour et al. 2017a) the notations τ M and q M referring to an effective thermal time constant and the amplitude of the maximum (located at σ = τ −1 M ), respectively (the factor 2 sets the maximal amplitude to q M ). This functional form corresponds to the so-called Maxwell model mentioned above. It describes the behaviour of an idealized forced oscillator composed of a string and a damper arranged in series (Greenberg 2009;Efroimsky 2012;Correia et al. 2014).
We note that other works based upon different approaches converged towards the functional form of the Maxwell model. For instance, Correia & Laskar (2001) used the parametrized function f (σ) = σ −1 1 − e −γσ 2 (γ being a real parameter, see Eq. (26) of the article) to mimic the behaviour of the atmospheric tidal torque, while Leconte et al. (2015) retrieved Eq. (19) empirically by analysing results obtained from simulations run with the LMDZ GCM.
An important remark should be made here concerning the behaviour of the tidal torque in the vicinity of the synchronization (i.e. for σ ≈ 0). To our knowledge, most of early works using the classical tidal theory to study the spin rotation of Venus and ignoring dissipative processes obtained a torque scaling as T ∝ σ −1 , and thus singular at the synchronization (e.g. Dobrovolskis & Ingersoll 1980;Correia & Laskar 2001. This is precisely the reason that led Correia & Laskar (2001) to introduce the regular ad hoc parametrized function mentioned above. Conversely, Ingersoll & Dobrovolskis (1978) and, later, Auclair-Desrotour et al. (2017a), derived a Maxwell-like tidal torque analytically by introducing a characteristic thermal time associated with boundary layer processes and radiative cooling. These early results may let think that dissipative processes are a necessary ingredient for a regular tidal torque to exist at the synchronization.
Although dissipative processes definitely regularize the atmospheric tidal torque at the synchronization (e.g. Auclair-Desrotour et al. 2017a), we showed in Sect. 4.2 that regularity also naturally emerges from approaches ignoring them when the vertical structure equation is solved in a self-consistent way. For a sufficiently small frequency, namely |σ| σ J , the torque derived from our analytic solution in the absence of dissipative mechanisms scales as T ∝ σ. Therefore, it seems that the singularity at σ = 0 obtained by early works could result from oversimplifying hypotheses, such as neglecting the three-dimensional aspect of the tidal response or tidal winds. For instance, we note that our analytical model asymptotically converges towards the function obtained by Dobrovolskis & Ingersoll (1980) when the vertical profile of tidal heating tends towards the Dirac distribution used by these authors (i.e. when b J → +∞).
The above statement means that the analytical solutions given by Eqs. (12) and (13) can be used in practice over the whole range of tidal frequencies without leading to unrealistic behaviours at the vicinity of synchronization, notwithstanding the fact that they were derived assuming that characteristic timescales associated with dissipative processes far exceed the tidal period.
In studies taking into account dissipative processes (e.g. Ingersoll & Dobrovolskis 1978;Auclair-Desrotour et al. 2017a), the parameter τ M of Eq. (19) can be interpreted as an effective timescale associated with the radiative cooling of the atmosphere in the Newtonian cooling approximation, where radiative losses are assumed to be proportional to temperature variations (Lindzen & McKenzie 1967;Auclair-Desrotour et al. 2017a;Auclair-Desrotour & Leconte 2018). These early analytical works established the following expression of the tidal torque (see e.g. Ingersoll & Dobrovolskis 1978, Eq. (2)), where ε stands for the effective fraction of the incoming flux absorbed by the atmosphere. Substituting δp 2,σ s;2 by Eq. (19) in Eq. (4) and comparing the obtained result with the preceding expression leads to a relationship between the Maxwell thermal time and maximum, which is the notation G referring to the gravitational constant.
Assuming that the atmosphere is optically thin in the visible frequency range and that the surface temperature corresponds to a black body equilibrium, we write the mean surface temperature as where we have introduced the Stefan-Boltzmann constant σ SB and the surface albedo A s . By substituting T s by Eq. (22) in Eq. (21), we obtain that the ratio q M /τ M does not depend on the surface pressure and scales as with ε = 1 − A s if the atmosphere is optically thick in the infrared. This relationship between τ M and q M means that the two parameters of the Maxwell model (Eq. (19)) can theoretically be reduced to the effective thermal timescale only, which is determined by complex boundary layer and dissipative processes in the general case. The scaling law given by Eq.
(23) will be tested using GCM simulations in Sect. 5. We now compare the Maxwell model to numerical results by assimilating the Maxwell amplitude and timescales to the maximum value of f odd and its associated timescale, respectively. The ab initio analytical solution given by Eqs. (12) and (13) ('Ana.') and its Maxwell-like form, derived for |σ| σ L and given by Eq. (17) ('Maxwell'), are both plotted in Fig. 3 as functions of the normalized forcing frequency (ω). The numerical values of σ L and σ TP used for the plot are determined by the eigenfrequency of the resonance associated with the Lamb mode in GCM simulations, that is ω L ≈ 260. We arbitrarily choose to set σ J = σ max (correspondence between the numerically-derived and the Maxwell maxima), which determines the value of b J (i.e. b J ≈ 14). Finally, the maximum q J is obtained by fitting the slope in the intermediate frequency-range to numerical results (q J ≈ 1042 Pa), and provides the value of the parameter J s (i.e. J s ≈ 0.05 W.kg −1 ). Figure 3 highlights the fact that the Maxwell model does not allow us to recover the behaviour of the torque in the low-frequency regime. The functional form given by numerical results and the Maxwell function clearly differ in this regime. Particularly, the maximal amplitude obtained from GCM simulations is about twice larger than that given by the model. We note that a smaller departure between the Maxwell and numerical maxima would certainly be obtained by fitting the Maxwell function to the whole spectrum of numerical results, and not only to the peak. However, this would also lead to overestimate the Maxwell timescale, and the fit would not be satisfactory either. A a consequence, a novel parametrized model has to be introduced to better describe the behaviour of the tidal torque in the low-frequency range. This is the purpose of the next section.

Introduction of a new parametrized model
It has been shown that the ab initio analytic model described in Sect. 4.2 and Appendix B reproduces the main features of the tidal torque in the high-frequency range, namely the resonance associated with the Lamb mode and the asymptotic scaling law T ∝ σ −1 . However, in the low-frequency range, the behaviour of the torque appears to be a little bit more complex than that predicted by the model, which reduces to a simple Maxwell function. This is not surprising since the atmospheric tidal response at low tidal frequencies involves complex non-linear mechanisms, interactions with mean flows, and dissipative processes, which are clearly outside of the scope of the classical tidal theory used to establish the solution given by Eqs. (12) and (13).
Yet, the frequency dependence of the tidal torque has to be characterized in the vicinity of synchronization as this is where its action of the planetary rotation is the strongest. Our effort has thus to be concentrated on the low-frequency regime and the transition with the high-frequency regime. As they treat the full non-linear 3D dynamics of the atmosphere in a selfconsistent way, GCM simulations are particularly useful in this prospect.
To make oneself an intuition of the behaviour of the torque, it is instructive to look at the logarithmic plot of Fig. 3 (bottom  left panel), which enables us to identify the different regimes at first glance. We basically observe two tendencies, highlighted in the plot by slopes taking the form of a straight line, in the zerofrequency limit (log (ω) 0.5) and the high-frequency asymptotic regime (log (ω) 1.5). In the interval 0.5 log (ω) 1.5, the tidal torque reaches a maximum and undergoes an abrupt decay.
Considering these observations, it seems relevant to approximate the logarithm of the torque by linear functions corresponding to the low and high-frequency regimes, and multiplied by sigmoid activation functions. By introducing the notation χ ≡ log ω, we thus define the parametrized function as where b trans ≈ log (q max ) is the level of the transition plateau, a 1 , b 1 , a 2 and b 2 the dimensionless coefficients of linear functions describing asymptotic regimes, and F 1 and F 2 two sigmoid activation functions expressed as In these expressions, the dimensionless parameters χ 1 and χ 2 designate the cutoff frequencies of F 1 and F 2 in logarithmic scale, and d 1 and d 2 the widths of transition intervals. The corresponding tidal torque is given by As the scaling law T ∝ σ −1 was derived from the ab initio model of Sect. 4.2 in the high-frequency range, we enforce it by setting a 2 = − 1. The eight left parameters are then obtained by fitting the function given by Eq. (24) to numerical results (as done previously, the odd function f odd is used). We thus end up with and plot the model function F par in Fig. 3 using these numerical values ('Param.').
As shown by Fig. 3, the parametrized function defined by Eq. (24) describes important features that escaped the Maxwell function, such as the fact that the tidal torque does not scale linearly with the forcing frequency in the zero-frequency limit, and the rapid decay characterizing the transition between low and high-frequency regimes.

Dependence of the tidal torque on the atmospheric composition
As it clearly has a strong impact, the dependence of the tidal torque on the atmospheric composition has to be discussed. In Appendix C, we treat the case of a CO 2 -dominated atmosphere with a mixture of water and sulphuric acid (H 2 SO 4 ) comparable to that hosted by the Venus planet. The obtained spectrum and the associated functions introduced above are plotted in Fig. 4, and shall be compared to those computed for the N 2 -dominated atmosphere, plotted in Fig. 3 (right panel). Several interesting features may be noted. First, the tidal torque exerted on the CO 2 -dominated atmosphere is more than twice weaker than that exerted on the N 2 -dominated atmosphere. Particularly, peaks are strongly attenuated. This results from the vertical distribution of tidal heating. Because of the optical thickness of carbon dioxide in the visible frequency range, an important part of the incoming stellar flux is absorbed above clouds. This is not the case of the N 2 -dominated atmosphere, where most of the flux reaches the planet surface and is re-emited in the infrared frequency range, leading to the thermal forcing of dense atmospheric layers located at high pressure levels.
Second, we observe a greater asymmetry between the negative and positive frequency ranges, the function f even being not negligible with respect to f odd . This is also an effect of the vertical distribution of tidal heating. In the case of the N 2 -dominated A17, page 9 of 22  atmosphere, most of the tidal torque is generated by density variations occurring at low altitudes, where the fluid is well coupled to the solid part of the planet by frictional forces. Switching from N 2 to CO 2 decreases the contribution of these layers, while it increases the contribution of layers located at pressure levels where the strong zonal jets mentioned above are generated. Despite the clear interest there is to study the tidal response of CO 2 -dominated atmospheres for the similarity of configurations they offer with the Venus planet, we choose to focus in this work on N 2 -dominated atmospheres owing to their simpler frequency-behaviour.

The surface-atmosphere coupling
The specific role played by the surface thermal response is not taken into account in linear models used to establish the Maxwell-like behaviour of the tidal torque described by Eq. (19) (e.g. Dobrovolskis & Ingersoll 1980;Auclair-Desrotour et al. 2017a). In these early works, the thermal forcing is assumed to be in phase with the stellar incoming flux, which amounts to considering that thermal tides are caused by the direct absorption of the flux. This approximation seems realistic in the case of Venuslike planets given that their atmospheres are optically thick in the visible range, and sufficiently dense to neglect their interactions with the surface.
However, it appears as a rough approximation in the case of optically thin atmospheres, where most of the stellar flux reaches the surface. In this case, thermal tides are mainly caused by the absorption of the flux emitted by the surface in the infrared range, which is delayed with respect to the stellar incoming flux owing to the surface inertia and dissipative processes such as thermal diffusion. Our N 2 -dominated atmosphere belongs to the second category. Thus, the role played by the thermal response of the ground should be considered in the present study to explain the observed difference between the obtained tidal torque and the Maxwell model. Table 2. Scaling laws of τ max and q max obtained using the LMDZ GCM for a dry terrestrial planet with a homogeneous N 2 atmosphere.

Case
Scaling laws of q max and τ max R 2 log (q max ) = −0.69 log ( Notes. The scaling law of q max /τ max is computed using the formers and should be compared to Eq. (35). Units: a is given in au, p s in bar, τ max in days, q max in Pa, the parameters of the linear fit α, β and R 2 are dimensionless.
Concerning this point, we note that Leconte et al. (2015) included the heat capacity of the surface C s in the simplified model they used to establish the Maxwell-like behaviour of the tidal torque (see Sect. 4 in the Material and Methods of their article). Hence, by introducing the heat capacity of the atmosphere/surface system C = C p p s /g + C s and the emission temperature (T e ), they expressed the relationship between surface temperature variations δT and the variations of the incoming stellar flux δF inc as with σ M = 4σ SB T 3 e /C (the subscript M refers to the Maxwell-like form of the function given by Eq. (28)). As we generally observe that T e ≈ T s in our GCM simulations of a 10 bar atmosphere (the mean surface temperature of the planet is well approximated by the black body equilibrium temperature, given by Eq. (22), in this case), this model implies that σ M should be always less than σ sup M = 4σ SB T 3 s g/ C p p s . However, in light of typical values of τ M obtained with the GCM (see Table 2), it appears that the above formula for σ sup M leads to underestimate σ M by a factor 10 to 100 for the case treated in the present study.
To understand the role played by the ground in the atmospheric tidal response, we adopt an ab initio approach describing thermal exchanges at the surface-atmosphere interface. Following along the line by Bernard (1962; see also Auclair-Desrotour et al. 2017a), we write the local budget of perturbative power inputs and losses, where we have introduced the small variations of the incoming stellar flux δF, surface temperature δT s , radiative heating by the atmosphere δF atm and diffusive losses in the ground δQ gr and in the atmosphere δQ atm . Owing to the absence of water, latent heats associated with changes of states are ignored.
In the general case, δT s and δQ atm are coupled with the atmospheric tidal response. Particularly, in the Newtonian cooling approximation (i.e. variations of the emitted flux are proportional to temperature variations), δF atm can be expressed as where K designates an effective coefficient of Newtonian cooling. In order to avoid mathematical complications, we ignore this A17, page 10 of 22  coupling by assuming either that |δF atm | 4σ SB T 3 s |δT s |, or, following Bernard (1962), that the variation of the atmospheric flux scales as δF atm ∝ δT s in a similar way as the variation of the flux emitted by the ground. This allows us to simplify radiative terms by writing where s ≈ 1 stands for the effective emissivity of the surface. With the above approximations, surface temperature variations can be written for a given mode as δT σ s = B σ s δF σ inc . We thus end up with (see detailed calculations in Appendix D) where B 0 s = 4σ SB T 3 s s −1 , and τ s designates the characteristic timescale of the surface thermal response, which depends on the thermal inertia of the ground I gr and of the atmosphere I atm at the interface, and is expressed as We compare this model to numerical results by extracting the Y 2 2 component of the surface temperature distribution δT 2,σ s;2 provided by GCM simulations, as previously done for the surface pressure distribution. The obtained values are plotted in the complex plane in Fig. 5. In this plot, the horizontal and vertical axes correspond to the real and imaginary parts of the normalized transfer function B σ s /B 0 s (such that δT 2,σ s;2 = B σ s δF 2,σ inc;2 ), respectively. Normalization is obtained by fitting numerical results with the function given by Eq. (32) in the low-frequency range (0 ≤ ω < 3). Figure 5 shows a good agreement between the functional form of the model and numerical results in the zero-frequency limit. However, we observe that the value of the thermal time τ s ∼ 0.3 days obtained by fitting Eq. (32) to numerical results in the low-frequency range is a decade smaller than the theoretical value given by Eq. (33), τ s ≈ 4.6 days (we use values given by Table 1, set s = 1 and neglect I atm ), which shows the limitations of the approach detailed above.
While the forcing frequency increases, the behaviour of the function interpolated using numerical results starts to change radically. In the vicinity of the resonance (σ ∼ σ max ), the imaginary part of B σ s decays abruptly whereas its theoretical analogous keeps growing. This divergence suggests a strong radiative coupling between the surface and the atmosphere, which comes from the fact that the emission of the atmosphere to the surface δF atm (see Eq. (30)) can no longer be neglected, as done in the model. The abrupt variation of the surface thermal lag around the resonance partially explain the behaviour of the tidal torque in this range. Nevertheless, to better understand it, one should study the whole dynamics of the atmospheric tidal response, which is beyond the scope of this work.
In the high-frequency range, that is for σ τ −1 s , the model predicts that the amplitude of temperature variations should tend to zero. Yet, we observe that δT 2,σ s;2 increases until reaching a maximum before decaying. This maximum corresponds here to a resonance whose frequency coincides with that of the main Lamb mode identified previously, in Sect. 4.2 (see Lamb 1917;Vallis 2006).

Exploration of the parameter space
We now examine the evolution of the tidal torque with the planet semi-major axis (a) and surface pressure (p s ).

Frequency spectra of the tidal torque
Considering the planet defined in Sect. 3.1, we carry out two studies. In study 1, we set p s = 10 bar and we compute frequency spectra of the imaginary part of the Y 2 2 -surface pressure component in the low-frequency range for a varying from 0.3 to 0.9 au. In study 2, we set a = a Venus , that is a = 0.723 au, and frequency spectra are computed for p s varying from 1 to 30 bar. The reference case characterized in the previous section, and parametrized by a = a Venus and p s = 10 bar, is located at the intersection of the two studies.
Limitations concerning the lower bound of the orbital radius range and the upper bound of the surface pressure range come from the spectra of optical properties used in simulations to compute radiative transfers (see Sect. 3.1), which were produced for temperatures below 710 K. Indeed, for a < 0.3 au or p s > 30 bar, the planet surface temperature exceeds this maximum. As this might lead to erroneous estimations of radiative transfers, we choose not to treat extremal cases, although there is no formal limitation for the GCM to run normally in these conditions. Radiative transfers also determine the convergence timescale necessary for the atmosphere to reach a steady state, P conv . For study 1, we use the timescale obtained in the reference case, that is 5800 Earth Solar days, considering that the steady state is reached more rapidly in mosts cases, where the planet is closer to the star (see Eq. (5)). Similarly, to take the dependence of P conv on the planet surface pressure into account in study 2, we set P conv to 1100, 2300, 5800 and 14000 Earth Solar days for p s = 1, 3, 10, 30 bar, respectively.
The obtained frequency spectra are plotted in Fig. 6 in linear (left) and logarithmic scales (right) for study 1 (top) and study 2 (bottom). In all plots, points designate the results of GCM simulations with the method described in Sect. 3, while solid  lines correspond to the associated cubic splines interpolations. The reference case (a = a Venus and p s = 10 bar) is designated by the solid grey line. Numerical values used to produce these plots are given by We retrieve here the features identified in Sect. 4. The tidal torque exhibits maxima located at the transition between the low-frequency and high-frequency asymptotic regimes. The corresponding peaks are slightly higher in the negative-frequency range than in the positive-frequency range owing to Coriolis effects and the impact of zonal jets on the angular lag of the tidal bulge. As expected, the amplitude of peaks increases with both the incoming stellar flux and the planet surface pressure. Interestingly, the evolution of q max and σ max with a and p s looks very regular. This suggests that the dependences of the peak maximum and characteristic timescale on the planet surface pressure and distance to star are well approximated by simple power scaling laws, and it is the case indeed, as shown in Sect. 5.2.
As previously noticed in the study of the reference case, the asymptotic behaviour of the tidal torque in the zero-frequency limit differs from that described by the Maxwell model. Particularly, the logarithmic plot of study 2 (bottom right panel) shows that the torque follows the scaling law f GCM (σ) ∝ σ 1/2 in cases characterized by low surface pressures, that is 1 and 3 bar. These cases correspond to the thin-atmosphere asymptotic limit, where thermal tides are driven by diffusion in the ground in the vicinity of the surface. We note that the simplified linear model of the surface thermal response detailed in Sect. 4.6 and Appendix D leads to a surface-generated radiative heating scaling as {δT s } ∝ σ 1/2 in the zero-frequency limit, which is precisely the dependence observed in Fig. 6.

Evolution of the thermal peak with the planet semi-major axis and atmospheric surface pressure
Let us now quantify the regular dependence of the peak of tidal torque on the planet orbital radius and surface pressure observed in the preceding section. We have thereby to determine how the two parameters defining the peak -namely its maximum value q max and associated timescale τ max -vary with a and p s . Thus, for each study, we fit numerical values of q max and τ max using a linear regression, formulated as where Y designates the logarithm of q max (Pa) or τ max (days), X the logarithm of a (au) or p s (bar), and α and β the dimensionless parameters of the fit. The values of these parameters are given by Table 2, as well as those of the corresponding coefficients of determination R 2 . We also compute log (q max /τ max ) for comparison with the theoretical scaling law given by Eq. (23)   Linear regressions are plotted in Fig. 7 (blue solid line). In order to provide an estimation of the variability of numerical results, error bars are given for the reference case. These error bars do not literally correspond to a margin of error, but indicate the resolution of the sampling for the frequency and maximum of the tidal torque. For q max , the amplitude of the error bar is the departure between the maxima of the interpolating function and data. For τ max , the two bounds of the error bar are the values associated with the nearest points of the sampling, designated by the subscripts inf and sup, such that τ inf ≤ τ max ≤ τ sup . These error bars depend on the ratio between the size of a frequency interval and the width of the thermal peak. For example, the thermal peak is undersampled for a = 0.3 au, which makes the fit less reliable in this case.
Comparing coefficients of determination in Table 2, we observe that a better fit is systematically obtained for q max than for τ max . This difference may be explained by the aspect of spectra displayed in Figs. 6 and 3. Since the peak of the tidal torque computed with the GCM is both flatter and larger than that of the Maxwell function, the position of the maximum is more sensitive to small fluctuations than the maximum itself. As a consequence, the variability of q max is less than the variability of τ max .
Hence, the linear regression fits particularly well the dependence of q max on a, while the plot of τ max exhibits a relatively important variability with respect to the linear tendency. We note however that differences with the fit are not significative since they remain small compared to the width of the peak. Concerning the dependence of τ max on a, one may also observe that the slope, given by α = 0.86, is almost twice smaller than that predicted by the scaling law of the radiative timescale given by Eq. (5), that is τ max ∝ n −1 ∝ a 3/2 .
As regards the ratio q max /τ max however, we recover numerically the scaling law predicted by the theoretical model (Eq. (23)) with a good approximation. This scaling law is numerically expressed in the units of Table 2 as if we assume that ε = 1 − A s (i.e. the flux reemitted by the ground is entirely absorbed by the atmosphere).
As may be seen, the dependence of q max /τ max on the surface pressure is small (α = 0.13) for want of being zero, as predicted by the model. Regarding the dependence on a, the relative difference between the numerical and theoretical values of α (i.e. 1.55 and 3/2, respectively) is around 3%. However, the value of β computed from GCM simulations (2.77) is higher than that predicted by the model (2.49), despite the fact that this latter is an upper estimation. This difference illustrates the limitations of the Maxwell model, which fails to describe the sharp variations of the tidal torque with the tidal frequency when |σ| ∼ σ max .

Scaling laws and generic formula for the tidal torque
By proceeding to a quantitative study of the evolution of the tidal torque maximum with the planet orbital radius and surface pressure, we demonstrated in the preceding section the regularity observed in Fig. 3. The scaling laws given by Table 2 and plotted in Fig. 7 show that frequency-spectra have the same aspect from the moment that the horizontal and vertical axes are rescaled following the obtained dependences on a and p s . In this section, our purpose is to compute this rescaling in a robust way, by taking into account the whole set of data at our disposal rather than the maximal value of the torque and the associated timescale only. Combining this rescaling with the parametrized model given by Eq. (24), we will obtain a novel generic formula for the frequency-behaviour of the thermally generated atmospheric tidal torque.
The parameter with respect to which axes are rescaled, a or p s , is denoted by p, and the considered case is subscripted j. A given family is thus composed of N p couples of numerical vectors σ j , T j , with 1 ≤ j ≤ N p (see Tables E.1 and E.2), associated with the value p j of the varying parameter p. For a given couple of vector, one may introduce the associated interpolating function We also introduce the renormalized vectorsσ j andT j , defined bŷ A17, page 13 of 22 A&A 624, A17 (2019) where α 1 and α 2 are the exponents characterizing the renormalization. The size of frequency domains covered by theσ j vectors varies with p in the general case. As a consequence, rescaling axes requires to define the bounds of the largest common interval, the notation N σ referring to the size of the frequency sampling (typically N σ = 21, see Tables E.1 and E.2), and σ j,1 and σ j,N σ to the lower and upper bounds of the interval sampled by σ j , respectively. The values of α 1 and α 2 are obtained by minimizing the squared difference function We note that the parameters derived from these calculations, q 0 and τ 0 , slightly differ from q max and τ max . They stand for the characteristic amplitude and timescale of the peak and not for its maximum value and corresponding forcing period, as q max and τ max . These parameters are defined by functions of p (that is a or p s ), by and for the second family (a = a Venus and variable p s ), with log (q 0 ) = 0.48 log (p s ) + 2.87, (44) log (τ 0 ) = 0.30 log (p s ) + 0.038.
Let us remind here the used units: a is given in au, p s in bar, q 0 in Pa, and τ 0 in days. The last step consists in combining these scaling laws with the parametrized model derived in the reference case (Eq. (24)). Proceeding to the change of variables associated with the renormalization, we obtain the generic parametrized model function δp 2,σ s;2 par ≡ q 0 10 F par (log|τ0σ|) sign (σ) .
We remind here that F par is the function where F 1 and F 2 are the activation functions defined by The parameters characterizing the generic formula given by Eq. (46) take the values The spectra of Fig. 6 are replotted in Fig. 8 using the normalized variables derived from the axes rescaling. In addition to numerical results and their interpolating functions, the tidal torque described by the generic parametrized model (Eq. (46)) is plotted as a function of the normalized tidal frequency τ 0 σ (dashed black line). Figure 8 clearly shows the relevance of the rescaling as regards the first family, where the dependence of the torque on the star-planet distance is investigated. After rescaling, spectra look similar and the model matches them fairly well. As regards the second family, we observe a greater variability of q 0 and τ 0 with a net separation between the reference and 30 bar cases. However, the frequency behaviour of the torque does not change much from one case to another and the parametrized function given by Eq. (46) remains a reasonable approximation of its main features.

Conclusions
In order to better understand the behaviour of the atmospheric torque created by the thermal tide, we computed the tidal response of the atmosphere hosted by a terrestrial planet using the LMDZ general circulation model. This work builds on both the early study by Leconte et al. (2015), which was a first attempt to characterize the atmospheric tidal response with this approach, and the early analytical works based upon the linear theory of atmospheric tides (e.g. Auclair-Desrotour et al. 2017a,b). It is motivated by the need to merge these two different approaches together in a self-consistent picture. Our aim was to proceed to a methodic comparison of their predictions while exploring the parameter space.
Hence, we considered the simplified case of a dry Venussized terrestrial planet orbiting a Sun-like star circularly and hosting a nitrogen-dominated atmosphere. Following the method by Leconte et al. (2015), we computed the atmospheric torque created by the semidiurnal thermal tide as a function of the tidal frequency by extracting the Y 2 2 component of the surface pressure anomaly in simulations.
As a first step, we characterized the variation of the torque with the forcing frequency for a reference case (p s = 10 bar and a = a Venus ), and explained its various features with an independent analytical model. As a second step, we explored the parameters space by focusing on the dependence of the tidal torque on the planet orbital radius and atmospheric surface pressure. The obtained results were then used to derive scaling laws characterizing the torque, renormalize the pressure anomaly and forcing period, and finally propose a novel generic parametrized function to model the frequency-behaviour of the torque in a realistic way in the case of a nitrogen-dominated atmosphere.
The first investigation confirmed and extended the results obtained by Leconte et al. (2015). We showed that the torque follows two different asymptotic regimes. In the high-frequency range, the torque decays inversely proportionally to the tidal frequency until it exhibits a resonance peak. These two features are both explained by the analytical solution derived using the ab initio linear theory of atmospheric tides. Particularly, the peak corresponds to a resonance associated with the Lamb mode, an acoustic type wave of wavelength comparable with the planet radius.  In the low-frequency range the torque, which is zero at synchronization, increases following a power law of index ranging from 0.5 to 0.7 until it reaches a maximum. While the increase and presence of a maximum are predicted by the analytical solution, the exponent and the value of amplitude of the peak differ significantly. These discrepancies result from the complex interactions between mechanisms laying beyond the scope of standard analytical treatments but resolved in 3D GCM simulations, such as the non-linear effects inherent to the atmospheric dynamics in the vicinity of synchronization, and the strong radiative coupling between the atmosphere and the planet surface. Typically, the low-frequency asymptotic regime of the tidal response is characterized by diurnal oscillations of large amplitude. The resulting differences in the day-and nightside temperature profiles significantly affect the stratification of the atmosphere. This clearly violates the small perturbation approximation upon which the analytic approach is based, and induces a non-linear coupling between the diurnal and semidiurnal oscillations that is important enough to modify the dependence of the tidal torque on the forcing frequency.
The parametrized function that we propose in the present work (given by Eq. (47)) appears as a good compromise as it matches numerical results in a more satisfactory way than the Maxwell model while being defined by a reasonably small number of parameters. It is thus perfectly suited to be implemented in evolutionary models of the rotational dynamics of a planet. Nevertheless, the Maxwell-like analytic solution derived by early studies (e.g. Ingersoll & Dobrovolskis 1978;Auclair-Desrotour et al. 2017a) provides a first order of magnitude approximation of the torque. It also predicts a relationship between the maximum of the thermal peak and the associated characteristic timescale. By establishing scaling laws governing the evolution of these features with the planet orbital radius and surface pressure, we retrieved numerically this relationship, which is q max /τ max ∝ a −3/2 .
The fact that scaling laws match well numerical results reveals that the torque and the tidal frequency can be normalized by the characteristic amplitude and frequency associated with the low-frequency regime. This was confirmed by the rescaling of spectra, which shows that numerical results obtained in all of the treated cases actually describe the same frequency dependence, whatever the star-planet distance and surface pressure. The combination of the parametrized function and scaling laws derived in this work thus leads to a generic empirical formula for the atmospheric tidal torque in the vicinity of synchronous rotation.
In spite of its limitations in the low-frequency regime, the analytic approach remains complementary with GCM calculations owing to the high computational cost of this method (several days of parallel computation on 80 processors are necessary to produce a spectrum with a sampling of 21 points in frequency). Results obtained from simulations can be used to improve the linear analysis, which provides in return a diagnosis of the physical and dynamical mechanisms involved in the tidal response.
As the study showed evidence of the interest of the numerical method using GCMs in characterizing the atmospheric tidal torque of terrestrial planets, several prospects can be considered for future works. First, the effects of clouds and optical thickness should be investigated owing to their strong impact on the tidal response. The case of an exo-Earth hosting a cloudy atmosphere may be treated in a similar way as the idealized planet of the present study. Second, it would be interesting to better characterize the dependence of the tidal torque on the atmospheric structure using ab initio analytic models. Third, numerical results and the derived generic parametrized function may be coupled to evolutionary models in order to quantify in a realistic way the contribution of the atmosphere to the evolution of the planet rotation over long timescales.
us to write V m,σ θ and V m,σ ϕ as functions of δp m,σ n /ρ 0 . By substituting horizontal winds by the obtained expressions in ∇ h ·V (Eq. (B.11)), and introducing the variable the whole system of governing equations given by Eqs. (B.5)-(B.10) can be put after some manipulations into the form (see Lindzen & Chapman 1969) Here, F m,σ is an operator depending on the x coordinate only and L m,ν the Laplace's tidal operator, which depends on the θ coordinate only and is formulated as (e.g. Lee & Saio 1997) L m,ν ≡ 1 sin θ ∂ θ sin θ 1 − ν 2 cos 2 θ ∂ θ (B.17) − 1 1 − ν 2 cos 2 θ mν 1 + ν 2 cos 2 θ 1 − ν 2 cos 2 θ + m 2 sin 2 θ , the quantity ν ≡ 2Ω/σ designating the so-called spin parameter. The above separation of coordinates allows us to expand the Fourier coefficients of G as and determine the equivalent depth of the mode associated with the triplet (n, m, σ) (e.g. Taylor 1936), (B.20) In the absence of resonances, the semidiurnal tidal response is generally dominated by the fundamental gravity mode, indicated by n = 0 3 , which corresponds to the associated Legendre function P 2 2 in the static case (e.g. Auclair-Desrotour & Leconte 2018). In the high-frequency regime, Λ 2,ν 0 ≈ Λ 2,1 0 ≈ 11.1, from the moment that n |Ω|. We note that this value, denoted by Λ 0 in the following, can be modified by dissipative processes. For instance, by including friction with the planet surface using a Raleigh drag of constant characteristic frequency σ R , one may show that the eigenvalue of the modes tends to the value of the static case, that is Λ 2,ν 0 ≈ 6, if σ/σ R → 0 (see e.g. Volland 1974;Auclair-Desrotour et al. 2017b).
As we focus on the n = 0 mode, we can drop the subscripts and superscripts n, m and σ to lighten notations. The function G m,σ n is now simply denoted by G, and so on for the tidal 3 We follow here the indexing notation by Lee & Saio (1997), which associates g-modes with positive n and r-modes to strictly negative n. heat source, pressure, density, temperature, wind velocity components, eigenvalues and equivalent depths. The usual change of variable G = e x/2 y leads to the vertical structure equation in its canonical form, where we have introduced the dimensionless vertical wavenumberk x , defined bŷ The vertical structure equation describe the behaviour of a forced harmonic oscillator, andk x thus corresponds to the inverse of a length scale of the variation of perturbed quantities across the vertical coordinate. Since the tidal response is adiabatic,k 2 x ∈ R, and its sign directly determines the nature of waves across the vertical axis. The conditionk 2 x > 0 indicates a propagating mode. Conversely,k 2 x < 0 corresponds to an evanescent mode. Computing analytic solutions turns out to be a very challenging problem except for a few simplified configurations. Therefore, we treat here the idealized case of the isothermal atmosphere, which is one of these configurations. We acknowledge the limitations of this academic atmospheric structure regarding real ones, where convective instability leads to a strong temperature gradient near the planet surface. However, this approach appears to be sufficient for the purpose of this appendix.
In the isothermal approximation, the temperature profile is supposed to be invariant with the radial coordinate. In light of Eq. (B.12), it immediately follows that H is a constant, dH/dx = 0, and The above expression shows the existence of a turning point for h = 4κH, where the sign ofk 2 x changes. This turning point occurs at the frequency In the reference case of the study, GCM simulations provide T s ≈ 316 K, which, combined with R s ≈ 297 J kg −1 K −1 , gives H ≈ 10.6 km in the isothermal approximation. An estimation of the normalized frequency ω TP = σ TP / (2n ) using Eq. (B.12) thus gives ω TP ≈ 270, showing that the turning point occurs in the high-frequency range and must therefore be taken into account in the calculation of an analytical solution. The conditionk 2 x > 0 (|σ| < σ TP ) corresponds to an oscillatory regime, whilek 2 x < 0 (|σ| > σ TP ) corresponds to an evanescent one. To solve the vertical structure equation, we have to choose a vertical profile for the tidal heat power per unit mass J. Following Lindzen et al. (1968), we opt for a profile of the form where J s stands for the heat absorbed at the planet surface and b J is a dimensionless optical depth characterizing the decay of heating across the vertical coordinate. This profile is derived from the Beer's law (e.g. Heng 2017) applied to an isothermal atmosphere The atmospheric torque generated by the thermal tide depends on the atmospheric composition, which has a strong impact on the vertical distribution of the tidal heating through clouds formation and the optical thickness of the gas mixture. In the study, we treat the case of a terrestrial planet hosting a cloudless N 2dominated atmosphere with a small amount of CO 2 . Hence, we ignore the effects of clouds and compute the thermal tide of an optically thin atmosphere in the visible frequency range, where the major part of the stellar flux reaches the planet surface without being absorbed.
In this appendix, we consider the case of a planet hosting a Venus-like CO 2 -dominated atmosphere with a mixture of water and sulphuric acid (H 2 SO 4 ) in the same reference configuration (a = a Venus and p s = 10 bar). We do not attempt to reproduce exactly the composition and dynamics of the Venus atmosphere, which is a complex problem beyond the scope of this study (see e.g. Lebonnois et al. 2010Lebonnois et al. , 2016, but to simply retrieve its main features (optical opacity, clouds absorption, etc.). As a consequence, we opt for a generic approach excluding a fine tuning of the atmospheric properties. We set the thermal capacity per unit mass of the gas (C p ) to 1000 J kg −1 K −1 , which is the typical value of C p in the case of Venus (e.g. Seiff et al. 1985), where the parameter decreases from 1181 J kg −1 K −1 near the surface to 904 J kg −1 K −1 at an altitude of 50 km.
Similarly, the mean molecular mass is set to 43.45 g mol −1 , and the volume mixing ratio of water vapour to 20 ppm (Moroz et al. 1979). We set the diameter of water particles to 3 µm, which is a typical value in the lower cloud (e.g. Knollenberg & Hunten 1980). To take into account the impact of sulphuric acid on the saturation pressure of water vapour p H 2 O , we use the prescription given by Gmitro & Vermeulen (1964) for aqueous sulphuric acid (see Eq. (24) of their article). This prescription is written as a function of the local temperature T , as where A, B, C, D, and E are empirical constants. The optical properties of the atmosphere used to compute radiative transfers are pre-computed using the HITRAN 2008 database (Rothman et al. 2009) for the Venus atmospheric mixture (instead of the Earth mixture used in the study). The spectrum of the atmospheric tidal torque due to the semidiurnal tide is plotted in Fig. C.1 with the spectrum of the N 2 reference case for comparison.
Because of the opacity of the atmosphere in the visible range, the fraction of the incoming stellar flux reaching the planet surface is less in the case of the CO 2 atmosphere than in the case of the N 2 atmosphere. Particularly, the resonance peak is strongly attenuated. We also observe a greater impact of Coriolis effects, the asymmetry of the tidal torque between negative and positive frequency ranges being more significant.
This difference can be explained by the vertical distribution of tidal heating. As mentioned above, the major part of the stellar flux reaches the planet surface if the atmosphere is composed of N 2 , which means that the tidal torque is mainly due to density variations occurring in the vicinity of the ground, where friction predominate over Coriolis forces. In the case of the CO 2 atmosphere, an important fraction of the incoming energy flux is absorbed at the cloud level. The contribution of this fraction is thus strongly affected by Coriolis effects through zonal mean flows characterizing the equilibrium dynamical state.

Appendix D: Simplified ab initio analytical model for the ground thermal response
As mentioned in Sect. 4.6, we follow along the line by Bernard (1962) to study the thermal response of the planet surface. In this approach, we consider the surface-atmosphere interface, located at the altitude z = 0, and write the power flux budget for a small perturbation in the framework of a frequency linear analysis. Hence, any quantity q can be expressed as q = q σ e iσt , where σ is the forcing frequency introduced in Sect. 2. In the following, we omit the superscript σ and use q in place of q σ , given that we work in the frequency domain. A variation of the effective incoming stellar flux (i.e. where the reflected component has been removed), denoted δF inc , is absorbed by the planet surface. A fraction δQ gr of this power is transmitted to the ground by thermal conduction, and an other fraction, δQ atm , is transmitted to the atmosphere through turbulent thermal diffusion. Finally, the increase of surface temperature δT s generated by δF inc induces a radiative emission, δF rad , which is expressed as δF rad = 4σ SB T 3 s δT s in the black body approximation (we recall that σ SB and T s are the Stefan-Boltzmann constant and mean surface temperature introduced in Sect. 3.2, respectively). Since the atmosphere is heated by both the incoming stellar flux and the surface thermal forcing, it undergoes a radiative cooling, similarly as the surface. The flux emitted downward to the surface is denoted δF atm . Thus, the power budget of the thermal perturbation at the interface is A17, page 20 of 22 P. Auclair-Desrotour et al.: Generic frequency dependence for the atmospheric tidal torque of terrestrial planets expressed as δF inc − 4σ SB T 3 s δT s + δF atm − δQ gr − δQ atm = 0. (D.1) To study the surface thermal response without having to consider the full atmospheric tidal response in its whole complexity, it is necessary to ignore the coupling induced by δF atm . This amounts to assuming either that the emission of the atmosphere towards the planet surface is negligible, or that it is proportional to δT s (see Bernard 1962). Thus, introducing the surface effective emissivity s , radiative terms can be reduced to 4σ SB T 3 s δT s − δF atm = 4σ SB T 3 s s δT s . (D. 2) The next step consists in defining the thermal exchanges resulting from diffusive processes, δQ gr and δQ atm . These flux are directly proportional to the gradient of the temperature profile anomaly in the vicinity of the interface, and are expressed as δQ gr = k gr (∂ z δT ) z=0 − , δQ atm = −k atm (∂ z δT ) z=0 + , (D.3) where ∂ z designates the partial derivative in altitude, δT the profile of temperature variations, and k gr and k atm the thermal conductivities of the ground and of the atmosphere at z = 0, respectively. By introducing the mean density profile of the planet ρ 0 and the thermal capacity per unit mass of the ground C gr (the analogous parameter for the atmosphere being C p ; see Sect. 3.1), the corresponding diffusivities can be defined by K gr ≡ k gr ρ 0 (0 − ) C gr and K atm ≡ k atm Temperature variations in the vicinity of the interface are described by the heat transport equation. We assume that diffusive processes predominates in the z → 0 limit. Moreover, since the typical horizontal length scale is far greater than the vertical one in the thin layer approximation, the horizontal component of the Laplacian describing diffusive processes can be neglected with respect to the vertical component in both the solid and atmospheric regions. It follows that iσδT = K gr ∂ zz δT, for z ≤ 0, (D.5) iσδT = K atm ∂ zz δT, for z > 0. (D.6) Solving these two equations with constant K gr and K atm , and ignoring the diverging term in solutions, we end up with δT (z) = δT s e [1+sign(σ)i]z/h σ gr , for z ≤ 0, (D.7) δT (z) = δT s e −[1+sign(σ)i]/h σ atm , for z > 0, (D.8) where we have introduced the frequency-dependent skin thicknesses of heat transport by thermal diffusion in the ground (h σ gr ) and atmosphere (h σ atm ), expressed as where τ s can be interpreted as the characteristic timescale of the surface thermal response. The parameter τ s is a function of the thermal inertia of the ground I gr ≡ ρ 0 (0 − ) C gr K gr and of the atmosphere I atm ≡ ρ 0 (0 + ) C p √ K atm , τ s ≡ 1 2 I gr + I atm 4σ SB T 3 s s 2 . (D.11) The above expression shows that τ s compares the efficiency of diffusive processes to that of the radiative cooling of the surface. The thermal time increases with the interface thermal inertia and decays when the surface temperature increases, scaling as τ s ∝ T −6 s . The expression of B σ s given by Eq. (D.10) highlights two asymptotic regimes. In the low-frequency regime, where |σ| τ −1 s , the surface responds instantaneously to the forcing δF, leading to a surface temperature oscillation in phase with the incoming stellar flux. At σ = 0, the incoming flux is equal to the radiative flux (δF inc = 4σ SB T 3 s s ), and B σ s = B 0 s . In the hightfrequency regime, where |σ| τ −1 s , the amplitude of the surface temperature variations decays and tends to zero in the limit |στ s | → +∞. At the transition, that is |σ| = τ −1 s , B σ s = (2/5) B 0 s and {B s } = − (1/5) B 0 s . As discussed in Sect. 4.6, the transfer function obtained using GCM simulations is well approximated by Eq. (D.10) in the low-frequency regime (see Fig. 5). But it diverges from the model when the forcing frequency increases, typically for |σ| τ −1 s . This divergences seems to result mainly from the fact that we ignored the radiative coupling between the atmosphere and the surface associated with δF atm , although it may be very strong. Particularly, this is the case for resonances of the atmospheric tidal response, where δF atm is increased similarly as the amplitude of pressure and temperature oscillations.

Appendix E: Tables of values obtained with GCM simulations for the exploration of the parameter space
The values used to plot the frequency-spectra of Fig. 6 are given by Table E.1 for study 1 (dependence on the star-planet distance) and Table E.2 for study 2 (dependence on the planet surface pressure). In both cases, the first column corresponds to the normalized tidal frequency ω = (Ω − n ) /n .  Notes. The first column corresponds to the normalized tidal frequency ω = (Ω − n ) /n . Notes. The first column corresponds to the normalized tidal frequency ω = (Ω − n ) /n .