© P. Auclair-Desrotour et al. 2019

1 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, 2003; Correia et al. 2003). Over the past decade, the growing performances of computers have made full numerical approaches affordable, and the atmospheric torque created by the thermal tidewas 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 frequency-dependence of the tidal torque predicted by ab initio analytical models (Ingersoll & Dobrovolskis 1978; Auclair-Desrotour et al. 2017a, 2018). 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 Maxwellmodel, 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 nitrogen-dominated 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 atmospherictidal 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₂-dominated atmosphere.We give our conclusions in Sect. 6.

2 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_{\star }=L_{\star }/\left(4\pi a^2\right)$, 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, $$\delta p_{\text{s}}\left(t,\theta ,\phi \right)={\displaystyle \sum _{l=1}^{+\infty }}{\displaystyle \sum _{m=-l}^l}\delta p_{\text{s};l}^{m,\sigma }{Y}_l^m\left(\theta ,\phi \right){\text{e}}^{i\sigma t},$$(1)

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_l^m$ the normalized spherical harmonics (see Appendix A), $\delta p_{\text{s};l}^{m,\sigma }$ the associated components, and $\sigma =m\left(\Omega -n_{\star }\right)$ the associated forcing frequencies (see e.g. Efroimsky 2012; Ogilvie 2014)¹.

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) $$\mathcal{T}=\frac{1}{g}{\displaystyle {\int }_{\mathcal{S}}^{}{\partial }_{\phi }{U}_{\text{T}}\delta p_{\text{s}}}\text{d}\mathcal{S},$$(2)

where the notation g refers to the surface gravity of the planet, ∂_φ to the partial derivative in longitude, $\mathcal{S}$ to the sphere of radius R_p, and $\text{d}\mathcal{S}=R_{\text{p}}^2\sin \theta \text{d}\theta \text{d}\phi$ to the surface element.

Similarly as the surface pressure anomaly, U_T can be expanded in Fourier series of time and spherical harmonics, $$U_{\text{T}}\left(t,\theta ,\phi \right)={\displaystyle \sum _{l=1}^{+\infty }}{\displaystyle \sum _{m=-l}^l}{U}_{\text{T};l}^{m,\sigma }{Y}_l^m\left(\theta ,\phi \right){\text{e}}^{i\sigma t},$$(3)

where the $U_{\text{T};l}^{m,\sigma }$ 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=\left|m\right|=2$. Besides, since the $U_{\text{T};l}^{m,\sigma }$ scale as $U_{\text{T};l}^{m,\sigma }\propto {\left(R_{\text{p}}/a\right)}^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=\left|m\right|=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), $$\mathcal{T}=\sqrt{\frac{24\pi }{5}}\frac{M_{\star }}{M_{\text{p}}}\frac{R_{\text{p}}^6}{a^3}\Im \left\{\delta p_{\text{s};2}^{2,\sigma }\right\},$$(4)

with $\sigma =2\left(\Omega -n_{\star }\right)$, and where the notation ℑ refers to the imaginary part of a complex number (ℜ referring to the real part). In this expression, $\delta p_{\text{s};2}^{2,\sigma }$ 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 $\eta =\text{sign}\left(\sigma \right)\Im \left\{\delta p_{\text{s};2}^{2,\sigma }\right\}$. When η < 0 (η > 0), the atmospheric tidal torque pushes the planet towards (away from) spin-orbit synchronization, $\left|\Omega -n_{\star }\right|$ decays (increases). Positions for which η = 0 correspond to the stable (${\text{d}\eta /\text{d}\sigma |}_{\text{eq}}<0$) or unstable (${\text{d}\eta /\text{d}\sigma |}_{\text{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.

3 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 $\delta p_{\text{s};2}^{2,\sigma }$ 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.

3.1 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 atmospherictidal response. The atmosphere is arbitrarily assumed to be nitrogen-dominated. However, a pure N₂-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₂ 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₂ for the mean molecular mass of the atmosphere, M_atm = 28.0134 g mol⁻¹ (Meija et al. 2016).

For a perfect diatomic gas, the ratio of heat capacities (also called first adiabatic exponent) is Γ₁ = 1.4, and it follows that $\kappa =\left({\Gamma }_1-1\right)/{\Gamma }_1=0.285$ (the parameter κ can also be written $\kappa ={\mathcal{R}}_{\text{GP}}/\left(M_{\text{atm}}{C}_p\right)$, where ${\mathcal{R}}_{\text{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⁻² s^−1∕2 K⁻¹, which is a typical value for Venus-like soils (see e.g. Lebonnois et al. 2010)². 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. 2010, 2011, 2013; Forget 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 wide range of temperatures and pressures using the HITRAN 2008 database (Rothman et al. 2009). These spectra are interpolatedevery 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.

3.2 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 $\left(a,p_{\text{s}}\right)$. As a first approximation, it depends on the radiative timescale of the deepest layers of the atmosphere τ_rad, which scalesas $${\tau }_{\text{rad}}\propto \frac{p_{\text{s}}}{g}\frac{C_p}{4{\sigma }_{\text{SB}}{T}_{\text{e}}^3},$$(5)

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_{\text{sol}}=2\pi /\left|\Omega -n_{\star }\right|$, except in the case of the spin-orbit synchronization (Ω = n_⋆), where there is no day-night 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 (seee.g. Fig. 1).

The third step consists in post-processing these data. We first remove the constant component, that is the mean surface 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 $\delta p_{\text{s};2}^{2,\sigma }$ 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 $\omega =\left(\Omega -n_{\star }\right)/n_{\star }$ 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 $\Delta \omega \equiv \left({\omega }_{\text{sup}}-{\omega }_{\text{inf}}\right)/N$. For instance, for the exploration of the parameters space detailed in Sect. 5, N = 20, ω_inf = − 30, ω_sup = 30, and thus Δω = 3.

Fig. 1 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).

4 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 this question, we consider the reference case (p_s = 10 bar and a = a_Venus).

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 $\omega =\left(\Omega -n_{\star }\right)/n_{\star }$ is increased from 0 (spin-orbit synchronization) to 24 (this corresponds to the length of the Solar day Psol = 9.36 days) for the reference case of the study (a = aVenus and ps = 10 bar).

4.1 Characterization of the reference case

In order to characterize the reference case, frequency-spectra of the atmospheric torque created by the semidiurnal thermal tideare computed in low-frequency and high-frequency ranges. For convenience, we introduce the function $f_{\text{GCM}}\left(\sigma \right)$, which is theinterpolating 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 $$f_{\text{odd}}\left(\sigma \right)=\frac{1}{2}\left[f_{\text{GCM}}\left(\sigma \right)-f_{\text{GCM}}\left(-\sigma \right)\right],$$(6)

which is the odd function f minimizing for any σ the distance defined by $$d\left(\sigma \right)=\left|\left|f_{\text{GCM}}\left(\sigma \right)-f\left(\sigma \right)\right|-\left|f_{\text{GCM}}\left(-\sigma \right)-f\left(-\sigma \right)\right|\right|.$$(7)

The complementary function f_even, such that f_GCM = f_odd + f_even, is defined by $$f_{\text{even}}=\frac{1}{2}\left[f_{\text{GCM}}\left(\sigma \right)+f_{\text{GCM}}\left(-\sigma \right)\right],$$(8)

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 $\omega =\sigma /\left(2n_{\star }\right)$ 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_{\text{max}}\equiv \max \left\{f_{\text{odd}}\left(\sigma \right)\right\}$ and the associated frequency σ_max, such that $f_{\text{odd}}\left({\sigma }_{\text{max}}\right)=q_{\text{max}}$, timescale ${\tau }_{\text{max}}\equiv {\sigma }_{\text{max}}^{-1}$ and normalized frequency ${\omega }_{\text{max}}={\sigma }_{\text{max}}/\left(2n_{\star }\right)$.

The tidal torque is negative for σ < 0 and positiveotherwise, 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 $\mathcal{T}\propto {\sigma }^{\alpha }$, with α ≈ 0.73. In the high-frequency range ($20\lesssim \left|\omega \right|\le 300$), it scales as $\mathcal{T}\propto {\sigma }^{-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 $\left|f_{\text{GCM}}\left(-\sigma \right)\right|>\left|f_{\text{GCM}}\left(\sigma \right)\right|$, and tends to vanish while $\left|\sigma \right|$ 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 $\left|\Omega \left(-\sigma \right)\right|\ne \left|\Omega \left(\sigma \right)\right|$ (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.

Fig. 3 Imaginary part of the $Y_2^2$ component of surface pressure variations as a function of the normalized forcing frequency $\omega =\left(\Omega -n_{\star }\right)/n_{\star }$ in the reference case (a = aVenus and ps = 10 bar). Top left panel: spectrum over the high-frequency range (− 300 < ω < 300) in linear scales. Top right panel: spectrum over the low-frequency range (− 30 < ω < 30) in linear scales. Bottom left panel: spectrum in logarithmic scales (for ω > 0). Numerical results obtained with GCM simulations are plotted (grey points), as well as the interpolating function fGCM (grey solid line), the odd and even functions of σ defined by Eqs. (6) and (8) (black and purple solid line, respectively), the function derived from the abinitio analytical model and given by Eqs. (12) and (13) (orange dashed line), its Maxwell-like approximation given by Eq. (17) (red dotted line), and the parametrized function given by Eq. (26) (blue dashed line).

4.2 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 atmospherictides (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 $$H=\frac{{\mathcal{R}}_{\text{s}}{T}_{\text{s}}}{g},$$(9)

where ${\mathcal{R}}_{\text{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_se^−b_Jx, 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), $${\sigma }_{\text{TP}}=\sqrt{\frac{4\kappa H{\Lambda }_{0}g}{R_{\text{p}}^2}}.$$(10)

The notation Λ₀ 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 Λ₀ can be approximated as a constant provided that $n_{\star }\ll |\Omega |$.

Hence, introducing the equivalent depth of the mode, $$h\equiv \frac{R_{\text{p}}^2{\sigma }^2}{{\Lambda }_{0}g},$$(11)

we obtain, for $\left|\sigma \right|\le {\sigma }_{\text{TP}}$, $$\Im \left\{\delta p_{\text{s}}\right\}={\sigma }^{-1}{p}_{\text{s}}\frac{\kappa J_{\text{s}}}{gH}\frac{\frac{H}{h}\left(b_{\text{J}}+\frac{1}{2}+\kappa \right)-\frac{1}{2}\left(b_{\text{J}}+1\right)}{\left[b_{\text{J}}\left(b_{\text{J}}+1\right)+\frac{\kappa H}{h}\right]\left(\frac{H}{h}-\frac{1}{{\Gamma }_1}\right)},$$(12)

and, for $\left|\sigma \right|>{\sigma }_{\text{TP}}$, $$\Im \left\{\delta p_{\text{s}}\right\}=\frac{{\sigma }^{-1}{p}_{\text{s}}\frac{\kappa J_{\text{s}}}{gh}\left[b_{\text{J}}+\frac{1}{2}\left(1-\sqrt{1-\frac{4\kappa H}{h}}\right)\right]}{\left[b_{\text{J}}\left(b_{\text{J}}+1\right)+\frac{\kappa H}{h}\right]\left[\frac{H}{h}-\frac{1}{2}\left(1+\sqrt{1-\frac{4\kappa H}{h}}\right)\right]}.$$(13)

We recall that $\kappa ={\mathcal{R}}_{\text{s}}/C_p$, where C_p designates the heat capacity per unit mass, and ${\Gamma }_1=1/\left(1-\kappa \right)$ 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 = Γ₁H (Lindzen & Blake 1972). In the asymptotic regime, where $n_{\star }\ll |\Omega |$, the characteristic Lamb frequency follows from Eq. (11), $${\sigma }_{\text{L}}=\sqrt{\frac{{\Lambda }_{0}g{h}_{\text{L}}}{R_{\text{p}}^2}}=\sqrt{\frac{{\Gamma }_1}{4\kappa }}{\sigma }_{\text{TP}}.$$(14)

By noticing that σ_L > σ_TP in the case of a diatomic gas (Γ₁ = 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 $\left|\sigma \right|={\sigma }_{\text{L}}$. The resonance hence occurs when the phase velocity of the forced mode equalizes the characteristic Lamb velocity $V_{\text{L}}=\sqrt{g{h}_{\text{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, Λ₀ ≈ 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 $\mathcal{T}\propto {\sigma }^{-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 $\left|\sigma \right|<{\sigma }_{\text{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 $\left|\sigma \right|\ll {\sigma }_{\text{TP}}$ and $n_{\star }\ll |\Omega |$ in the meantime. As H∕h ∝ σ⁻², the preceding condition implies that H∕h ≫ 1. It follows that $$\frac{\frac{H}{h}\left(b_{\text{J}}+\frac{1}{2}+\kappa \right)-\frac{1}{2}\left(b_{\text{J}}+1\right)}{\frac{H}{h}-\frac{1}{{\Gamma }_1}}~b_{\text{J}}+\frac{1}{2}+\kappa .$$(15)

By invoking the strong optical thickness of the atmosphere in the infrared (b_J ≫ 1), we remark that we recover analytically the scaling law $\mathcal{T}\propto {\sigma }^{-1}$ observed in Fig. 3 from the moment that the condition $1\ll H/h\ll {\kappa }^{-1}{b}_{\text{J}}^2$ is satisfied. This provides a definition for the intermediate frequency-range, which is now the range corresponding to ${\sigma }_{\text{J}}\ll \left|\sigma \right|\ll {\sigma }_{\text{L}}$, where we have introduced the thermal frequency $${\sigma }_{\text{J}}=\frac{{\sigma }_{\text{TP}}}{2\sqrt{b_{\text{J}}\left(b_{\text{J}}+1\right)}}\ll {\sigma }_{\text{TP}}.$$(16)

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 $\left|\sigma \right|\ll {\sigma }_{\text{L}}$ (or H∕h ≫ 1), Eq. (12) can be approximated by the function $$\Im \left\{\delta p_{\text{s}}\right\}\approx \frac{2q_{\text{J}}{\tau }_{\text{J}}\sigma }{1+{\left({\tau }_{\text{J}}\sigma \right)}^2},$$(17)

where the associated characteristic timescale τ_J and maximal amplitude of the pressure anomaly q_J are $$\begin{array}{lll}{\tau }_{\text{J}}={\sigma }_{\text{J}}^{-1},\hfill & \text{and}\hfill & q_{\text{J}}=\frac{\kappa J_{\text{s}}\left(b_{\text{J}}+\frac{1}{2}+\kappa \right)}{gH{\sigma }_{\text{TP}}\sqrt{b_{\text{J}}\left(b_{\text{J}}+1\right)}}{p}_{\text{s}}.\hfill \end{array}$$(18)

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.

4.3 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 $\sigma =2\left(\Omega -n_{\star }\right)$ as (e.g. Ingersoll & Dobrovolskis 1978; Auclair-Desrotour et al. 2017a) $$\Im {\left\{\delta p_{\text{s};2}^{2,\sigma }\right\}}_{\text{M}}=\frac{2q_{\text{M}}{\tau }_{\text{M}}\sigma }{1+{\left({\tau }_{\text{M}}\sigma \right)}^2},$$(19)

the notations τ_M and q_M referring to an effective thermal time constant and the amplitude of the maximum (located at $\sigma ={\tau }_{\text{M}}^{-1}$), 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\left(\sigma \right)={\sigma }^{-1}\left(1-{\text{e}}^{-\gamma {\sigma }^2}\right)$ (γ 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 $\mathcal{T}\propto {\sigma }^{-1}$, and thus singular at the synchronization (e.g. Dobrovolskis & Ingersoll 1980; Correia & Laskar 2001, 2003). 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 $\left|\sigma \right|\ll {\sigma }_{\text{J}}$, the torque derived from our analytic solution in the absence of dissipative mechanisms scales as $\mathcal{T}\propto \sigma$. 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)), $$\mathcal{T}=\frac{3\pi n_{\star }^2{R}_{\text{p}}^4\epsilon F_{\star }}{8C_p{T}_{\text{s}}}\frac{\sigma }{{\sigma }^2+{\tau }_{\text{M}}^{-2}},$$(20)

where ε stands for the effective fraction of the incoming flux absorbed by the atmosphere. Substituting $\Im \left\{\delta p_{\text{s};2}^{2,\sigma }\right\}$ 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 $$\frac{q_{\text{M}}}{{\tau }_{\text{M}}}=\frac{3}{128}\sqrt{\frac{5}{6\pi }}\frac{\mathcal{G}{M}_{\text{p}}\epsilon L_{\star }}{C_p{T}_{\text{s}}{R}_{\text{p}}^2{a}^2},$$(21)

the notation $\mathcal{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 $$T_{\text{s}}={\left[\frac{\left(1-A_{\text{s}}\right)L_{\star }}{16\pi {\sigma }_{\text{SB}}{a}^2}\right]}^{\frac{1}{4}},$$(22)

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 $$\frac{q_{\text{M}}}{{\tau }_{\text{M}}}\approx \frac{3}{128}\sqrt{\frac{5}{6\pi }}\frac{\mathcal{G}{M}_{\text{p}}\epsilon L_{\star }}{C_p{R}_{\text{p}}^2}{\left[\frac{\left(1-A_{\text{s}}\right)L_{\star }}{16\pi {\sigma }_{\text{SB}}}\right]}^{-\frac{1}{4}}{a}^{-3/2}.$$(23)

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 $\left|\sigma \right|\ll {\sigma }_{\text{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⁻¹).

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.

4.4 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 $\mathcal{T}\propto {\sigma }^{-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 self-consistent 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 zero-frequency limit ($\log \left(\omega \right)\lesssim 0.5$) and the high-frequency asymptotic regime ($\log \left(\omega \right)\gtrsim 1.5$). In the interval $0.5\lesssim \log \left(\omega \right)\lesssim 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 $$\begin{array}{ll}{\mathcal{F}}_{\text{par}}\left(\chi \right)\equiv \hfill & \left(a_1\chi +b_1\right){\mathcal{F}}_1\left(\chi \right)+\left(a_2\chi +b_2\right){\mathcal{F}}_2\left(\chi \right)\hfill \\ \hfill & +b_{\text{trans}}\left[1-{\mathcal{F}}_1\left(\chi \right)-{\mathcal{F}}_2\left(\chi \right)\right],\hfill \end{array}$$(24)

where $b_{\text{trans}}\approx \log \left(q_{\text{max}}\right)$ is the level of the transition plateau, a₁, b₁, a₂ and b₂ the dimensionless coefficients of linear functions describing asymptotic regimes, and ${\mathcal{F}}_1$ and ${\mathcal{F}}_2$ two sigmoid activation functions expressed as $$\begin{array}{lll}{\mathcal{F}}_1\left(\chi \right)\equiv \frac{1}{1+{\text{e}}^{\left(\chi -{\chi }_1\right)/d_1}},\hfill & \text{and}\hfill & {\mathcal{F}}_2\left(\chi \right)\equiv \frac{1}{1+{\text{e}}^{-\left(\chi -{\chi }_2\right)/d_2}}.\hfill \end{array}$$(25)

In these expressions, the dimensionless parameters χ₁ and χ₂ designate the cutoff frequencies of ${\mathcal{F}}_1$ and ${\mathcal{F}}_2$ in logarithmic scale, and d₁ and d₂ the widths of transition intervals. The corresponding tidal torque is given by $${\mathcal{T}}_{\text{par}}={10}^{{\mathcal{F}}_{\text{par}}\left(\log \left|\omega \right|\right)}\text{sign}\left(\omega \right).$$(26)

As the scaling law $\mathcal{T}\propto {\sigma }^{-1}$ was derived from the ab initio model of Sect. 4.2 in the high-frequency range, we enforce it by setting a₂ = − 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 $$\begin{array}{llll}{a}_1=0.734,\hfill & b_1=2.85,\hfill & d_1=0.0100,\hfill & {\chi }_1=0.637,\hfill \\ a_2=-1,\hfill & b_2=4.23,\hfill & d_2=0.0232,\hfill & {\chi }_2=1.20,\hfill \\ b_{\text{trans}}=3.33,\hfill & \hfill & \hfill & \hfill \end{array}$$(27)

and plot the model function ${\mathcal{F}}_{\text{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.

4.5 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 compositionhas to be discussed. In Appendix C, we treat the case of a CO₂-dominated atmosphere with a mixture of water and sulphuric acid (H₂SO₄) 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₂-dominated atmosphere, plotted in Fig. 3 (right panel). Several interesting features may be noted.

First, the tidal torque exerted on the CO₂-dominated atmosphere is more than twice weaker than that exerted on the N₂-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₂-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 respectto f_odd. This is also an effect of the vertical distribution of tidal heating. In the case of the N₂-dominated 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₂ to CO₂ decreases the contribution of these layers, while it increases the contribution of layers located at pressure levels where the strong zonaljets mentioned above are generated. Despite the clear interest there is to study the tidal response of CO₂-dominated atmospheresfor the similarity of configurations they offer with the Venus planet, we choose to focus in this work on N₂-dominated atmospheres owing to their simpler frequency-behaviour.

Fig. 4 Imaginary part of the $Y_2^2$ component of surface pressure variations as a function of the normalized forcing frequency $\omega =\left(\Omega -n_{\star }\right)/n_{\star }$ in the reference case (a = aVenus and ps = 10 bar) with a CO2-dominated atmosphere.Numerical results obtained with GCM simulations are plotted (grey points), as well as the interpolating function fGCM (dashed greyline), the odd and even functions of σ defined by Eqs. (6) and (8) (solid black and purple lines, respectively), the Maxwell function given by Eq. (17) (red dotted line), and the parametrized function given by Eq. (24) (blue dashed line) with the parameters: a1 = 0.549, b1 = 2.52, a2 = − 1, b2 = 3.94, χ1 = 0.798, d1 = 0.072, χ2 = 1.35, d2 = 0.0098, and btrans = 2.59.

4.6 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 Venus-like 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₂-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.

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_pp_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 $$\delta T_{\text{s}}=\frac{\delta \text{}\text{}{F}_{\text{inc}}}{{\sigma }_{\text{M}}C}\frac{1}{1+i\sigma /{\sigma }_{\text{M}}},$$(28)

with ${\sigma }_{\text{M}}=4{\sigma }_{\text{SB}}{T}_{\text{e}}^3/C$ (the subscript M refers to theMaxwell-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 $\sigma M^{\text{sup}}=4{\sigma }_{\text{SB}}{T}_{\text{s}}^{3}g/\left(C_p{p}_{\text{s}}\right)$. However, in light of typical values of τ_M obtained with the GCM (see Table 2), it appears that the above formula for σM^sup 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, $$\delta \text{}\text{}{F}_{\text{inc}}-4{\sigma }_{\text{SB}}{T}_{\text{s}}^3\delta T_{\text{s}}+\delta \text{}\text{}{F}_{\text{atm}}-\delta Q_{\text{gr}}-\delta Q_{\text{atm}}=0,$$(29)

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 temperaturevariations), δF_atm can be expressed as $$\delta \text{}\text{}{F}_{\text{atm}}={\displaystyle {\int }_0^{+\infty }K\left(x,\theta ,\phi \right)\delta T\left(x,\theta ,\phi \right)}\text{d}x,$$(30)

where K designates an effective coefficient of Newtonian cooling. In order to avoid mathematical complications, we ignore this coupling by assuming either that $\left|\delta \text{}\text{}{F}_{\text{atm}}\right|\ll 4{\sigma }_{\text{SB}}{T}_{\text{s}}^3\left|\delta T_{\text{s}}\right|$, 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 $$4{\sigma }_{\text{SB}}{T}_{\text{s}}^3\delta T_{\text{s}}-\delta \text{}\text{}{F}_{\text{atm}}=4{\sigma }_{\text{SB}}{T}_{\text{s}}^3{\varepsilon }_{\text{s}}\delta T_{\text{s}},$$(31)

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 $\delta T_{\text{s}}^{\sigma }={\mathfrak{B}}_{\text{s}}^{\sigma }\delta \text{}\text{}{F}_{\text{inc}}^{\sigma }$. We thus end up with (see detailed calculations in Appendix D) $${\mathfrak{B}}_{\text{s}}^{\sigma }=\frac{{\mathfrak{B}}_{\text{s}}^0}{1+\left[1+\text{sign}\left(\sigma \right)i\right]\sqrt{{\tau }_{\text{s}}\left|\sigma \right|}},$$(32)

where ${\mathfrak{B}}_{\text{s}}^0={\left(4{\sigma }_{\text{SB}}{T}_{\text{s}}^3{\varepsilon }_{\text{s}}\right)}^{-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 $${\tau }_{\text{s}}=\frac{1}{2}{\left(\frac{I_{\text{gr}}+I_{\text{atm}}}{4{\sigma }_{\text{SB}}{T}_{\text{s}}^3{\varepsilon }_{\text{s}}}\right)}^2.$$(33)

We compare this model to numerical results by extracting the $Y_2^2$ component of the surface temperature distribution $\delta T_{\text{s};2}^{2,\sigma }$ 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 ${\mathfrak{B}}_{\text{s}}^{\sigma }/{\mathfrak{B}}_{\text{s}}^0$ (such that $\delta T_{\text{s};2}^{2,\sigma }={\mathfrak{B}}_{\text{s}}^{\sigma }\delta \text{}\text{}{F}_{\text{inc};2}^{2,\sigma }$), 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 showsthe 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 ${\mathfrak{B}}_{\text{s}}^{\sigma }$ 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 $\sigma \gg {\tau }_{\text{s}}^{-1}$, the model predicts that the amplitude of temperature variations should tend to zero. Yet, we observe that $\left|\delta T_{\text{s};2}^{2,\sigma }\right|$ 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).

Fig. 5 Nyquist plot of the transfer function ${\mathfrak{B}}_{\text{s}}^{\sigma }$ associated with the $Y_2^2$ component of the semidiurnal tide, that is such that $\delta T_{\text{s};2}^{2,\sigma }={\mathfrak{B}}_{\text{s}}^{\sigma }\delta \text{}\text{}{F}_{\text{inc};2}^{2,\sigma }$. The imaginary part of the normalized function ${\mathfrak{B}}_{\text{s}}^{\sigma }/\left|{\mathfrak{B}}_{\text{s}}{}^0\right|$ is plotted in absolute value as a function of the real part in the reference case of the study (a = aVenus and ps = 10 bar). Values obtained using GCM simulations are indicated by points. They are interpolated by a green line in the range 0 ≤ ω < 3 (step Δω = 0.3), a black line in the range 3 ≤ ω < 30 (step Δω = 3), and a blue line in the range 30 ≤ ω ≤ 300 (step Δω = 30), $\omega =\left(\Omega -n_{\star }\right)/n_{\star }$ being the normalized tidal frequency. Similarly, values obtained in these ranges with the model given by Eq. (32) for the surface thermal time τs = 0.3 days are designated by crosses. The red line corresponds to the function itself.

5 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).

5.1 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 atthe 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 Table E.1 for study 1 and Table E.2 for study 2.

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_{\text{GCM}}\left(\sigma \right)\propto {\sigma }^{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 $\Im \left\{\delta T_{\text{s}}\right\}\propto {\sigma }^{1/2}$ in the zero-frequency limit, which is precisely the dependence observed in Fig. 6.

Fig. 6 Imaginary part of the $Y_2^2$ harmonic of surface pressure oscillations (Pa) as a function of the normalized tidal frequency $\omega =\left(\Omega -n_{\star }\right)/n_{\star }$ in linear (leftpanels) and logarithmic scales (for ω > 0). Top panels: spectra obtained with GCM numerical simulations for a fixed surface pressure, ps = 10 bar, and various values of the star–planet distance: a = 0.3–0.9 au with a step Δa = 0.1 au. Bottom panels: spectra obtained for a fixed star planet–planet distance, a = aVenus and various values of the surface pressure : ps = 1, 3, 10, 30 bar. Numerical results are designated by points and interpolated with cubic splines. The reference case of the study (a = aVenus and ps = 10 bar) corresponds to the grey line in all plots.

5.2 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 $$Y=\alpha X+\beta ,$$(34)

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². We also compute $\log \left(q_{\text{max}}/{\tau }_{\text{max}}\right)$ for comparison with the theoretical scaling law given by Eq. (23) in the case where the dependence of the tidal torque on the forcing frequency is approximated by a Maxwell function.

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 ${\tau }_{\text{max}}\propto n_{\star }^{-1}\propto 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 $$\log \left(\frac{q_{\text{M}}}{{\tau }_{\text{M}}}\right)=-\frac{3}{2}\log \left(a\right)+2.49,$$(35)

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 $\left|\sigma \right|~{\sigma }_{\text{max}}$.

Fig. 7 Parameters of the Maxwell model given by Eq. (19) computed from GCM simulations as functions of the star-planet distance a (left panels, in au) and ps (right panels, in bar) in logarithmic scales. Top panels: amplitude qmax of the maximum (Pa). Bottom panels: associated characteristic thermal time τmax (days). Numerical results obtained with GCM simulations are designated by black points, the corresponding linear regressions (see Table 2) by solid blue lines. For the reference case, error bars are given to indicate the resolution ofthe sampling (for details about how they are constructed, see Sect. 5.2).

5.3 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 𝔭, and the considered case is subscripted j. A given family is thus composed of N_𝔭 couples of numerical vectors $\left({\sigma }_j,{\mathcal{T}}_j\right)$, with 1 ≤ j ≤ N_𝔭 (see Tables E.1 and E.2), associated with the value 𝔭_j of the varying parameter 𝔭. For a given couple of vector, one may introduce the associated interpolating function $$f_j\left(\sigma \right)\equiv f\left({\sigma }_j,{\mathcal{T}}_j,\sigma \right).$$(36)

We also introduce the renormalized vectors ${\widehat{\sigma }}_j$ and ${\widehat{\mathcal{T}}}_j$, defined by $$\begin{array}{lll}{\widehat{\sigma }}_j\equiv {\mathfrak{p}}_j^{-{\alpha }_1}{\sigma }_j,\hfill & \text{and}\hfill & {\widehat{\mathcal{T}}}_j\equiv {\mathfrak{p}}_j^{-{\alpha }_2}{\mathcal{T}}_j.\hfill \end{array}$$(37)

where α₁ and α₂ are the exponents characterizing the renormalization.

The size of frequency domains covered by the ${\widehat{\sigma }}_j$ vectors varies with 𝔭 in the general case. As a consequence, rescaling axes requires to define the bounds of the largest common interval, $$\begin{array}{l}{\widehat{\sigma }}_{\text{inf}}\left({\alpha }_1\right)=\max {\left\{{\mathfrak{p}}_j^{-{\alpha }_1}{\sigma }_{j,1}\right\}}_j,\\ {\widehat{\sigma }}_{\text{sup}}\left({\alpha }_1\right)=\min {\left\{{\mathfrak{p}}_j^{-{\alpha }_1}{\sigma }_{j,N_{\sigma }}\right\}}_j,\end{array}$$

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 ${\sigma }_{j,N_{\sigma }}$ to the lower and upper bounds of the interval sampled by σ_j, respectively.

The values of α₁ and α₂ are obtained by minimizing the squared difference function $$F\left({\alpha }_1,{\alpha }_2\right)=\frac{{\displaystyle \sum _{j<k}}{\displaystyle {\int }_{{\widehat{\sigma }}_{\text{inf}}\left({\alpha }_1\right)}^{{\widehat{\sigma }}_{\text{sup}}\left({\alpha }_1\right)}\text{}\text{}{\left[{\mathfrak{p}}_j^{-{\alpha }_2}{f}_j\left({\mathfrak{p}}_j^{{\alpha }_1}\widehat{\sigma }\right)-{\mathfrak{p}}_k^{-{\alpha }_2}{f}_k\left({\mathfrak{p}}_k^{{\alpha }_1}\widehat{\sigma }\right)\right]}^2\text{}\text{}}\text{d}\widehat{\sigma }}{{\widehat{\sigma }}_{\text{sup}}\left({\alpha }_1\right)-{\widehat{\sigma }}_{\text{inf}}\left({\alpha }_1\right)}.$$(40)

We note that the parameters derived from these calculations, q₀ and τ₀, 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 𝔭 (that is a or p_s), by $$\begin{array}{lll}{\tau }_0\equiv {\tau }_{\text{ref}}{\left(\frac{\mathfrak{p}}{{\mathfrak{p}}_{\text{ref}}}\right)}^{-{\alpha }_1},\hfill & \text{and}\hfill & q_0\equiv q_{\text{ref}}{\left(\frac{\mathfrak{p}}{{\mathfrak{p}}_{\text{ref}}}\right)}^{{\alpha }_2},\hfill \end{array}$$(41)

where 𝔭_ref designates the value of 𝔭 in the reference case (typically a_Venus or 10 bar), q_ref = 2.24 kPa the corresponding maximum of the peak and τ_ref = 2.18 days the associated timescale. Hence, we end up for the first family (p_s = 10 bar and variable a) with the scaling laws $$\begin{array}{l}\log \left(q_0\right)=-0.61\log \left(a\right)+3.26,\\ \log \left({\tau }_0\right)=0.68\log \left(a\right)+0.43,\end{array}$$

and for the second family (a = a_Venus and variablep_s), with $$\begin{array}{l}\log \left(q_0\right)=0.48\log \left(p_{\text{s}}\right)+2.87,\\ \log \left({\tau }_0\right)=0.30\log \left(p_{\text{s}}\right)+\mathrm{0.038.}\end{array}$$

Let us remind here the used units: a is given in au, p_s in bar, q₀ in Pa, and τ₀ 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 $$\Im {\left\{\delta p_{\text{s};2}^{2,\sigma }\right\}}_{\text{par}}\equiv q_0{10}^{{\mathcal{F}}_{\text{par}}\left(\log \left|{\tau }_0\sigma \right|\right)}\text{sign}\left(\sigma \right).$$(46)

We remind here that ${\mathcal{F}}_{\text{par}}$ is the function $$\begin{array}{ll}{\mathcal{F}}_{\text{par}}\left(\chi \right)\equiv \hfill & \left(a_1\chi +b_1\right){\mathcal{F}}_1\left(\chi \right)+\left(a_2\chi +b_2\right){\mathcal{F}}_2\left(\chi \right)\hfill \\ \hfill & +b_{\text{trans}}\left[1-{\mathcal{F}}_1\left(\chi \right)-{\mathcal{F}}_2\left(\chi \right)\right],\hfill \end{array}$$(47)

where ${\mathcal{F}}_1$ and ${\mathcal{F}}_2$ are the activation functions defined by $$\begin{array}{lll}{\mathcal{F}}_1\left(\chi \right)\equiv \frac{1}{1+{\text{e}}^{\left(\chi -{\chi }_1\right)/d_1}},\hfill & \text{and}\hfill & {\mathcal{F}}_2\left(\chi \right)\equiv \frac{1}{1+{\text{e}}^{-\left(\chi -{\chi }_2\right)/d_2}}.\hfill \end{array}$$(48)

The parameters characterizing the generic formula given by Eq. (46) take the values $$\begin{array}{llll}{a}_1=0.734,\hfill & b_1=0.171,\hfill & d_1=0.010,\hfill & {\chi }_1=-0.277,\hfill \\ a_2=-1,\hfill & b_2=-0.031,\hfill & d_2=0.023,\hfill & {\chi }_2=0.290,\hfill \\ b_{\text{trans}}=-\mathrm{0.020.}\hfill & \hfill & \hfill & \hfill \end{array}$$(49)

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 τ₀ σ (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₀ and τ₀ 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.

6 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 Venus-sized 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 atmospherictorque 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 night-side 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.

Fig. 8 Normalized imaginary part of the $Y_2^2$ harmonic of surface pressure oscillations as a function of the normalized tidal frequency in linear (left panels) and logarithmic scales (for ω > 0). Top panels: spectra obtained with GCM numerical simulations for a fixed surface pressure, ps = 10 bar, and various values of the star–planet distance : a = 0.3–0.9 au with a step Δa = 0.1 au. Bottom panels: spectra obtained for a fixed star planet–planet distance, a = aVenus and various values of the surface pressure : ps = 1, 3, 10, 30 bar. Results are designated by points and interpolated with cubic splines. The reference case of the study (a = aVenus and ps = 10 bar) corresponds to the grey line in all plots. The generic formula of the tidal torque derived in this study and given by Eq. (46) is designated by the black dashed line. The characteristic torque q0 and timescale τ0 used for the rescaling are given by Eqs. (42)–(45) as functions of a and ps.

Acknowledgements

The authors acknowledge funding by the European Research Council through the ERC grant WHIPLASH 679030. They also thank the anonymous reviewer for helpful comments that improved the manuscript.

Appendix A Normalized spherical harmonics

With the notations introduced in Sect. 2, the normalized spherical harmonics associated with the degrees l and m ($l\in \mathbb{N}$ and − l ≤ m ≤ l) are defined by (e.g. Arfken & Weber 2005) $$Y_l^m\left(\theta ,\phi \right)\equiv {\left(-1\right)}^m\sqrt{\frac{2l+1}{4\pi }\frac{\left(l-m\right)!}{\left(l+m\right)!}}{P}_l^m\left(\cos \theta \right){\text{e}}^{im\phi },$$(A.1)

where the $P_l^m$ designate the associated Legendre functions, expressed, for $x\in \left[-1,1\right]$, as $$P_l^m\left(x\right)\equiv {\left(-1\right)}^m{\left(1-x^2\right)}^{m/2}\frac{{\text{d}}^m}{\text{d}{x}^m}{P}_l\left(x\right),$$(A.2)

the P_l being the Legendre polynomials, defined by $$P_l\left(x\right)\equiv \frac{1}{2^{l}l!}\frac{{\text{d}}^1}{\text{d}{x}^1}\left[{\left(x^2-1\right)}^l\right].$$(A.3)

Appendix B Linear analytical model for the high-frequency regime

In Sect. 4.2, we give an expression of the surface pressure anomaly as a function of the tidal frequency (see Eqs. (12) and (13)). This expression is a solution of the thermally generated tidal response derived in the framework of the linear theory of atmospheric tides, as described for instance in Chapman & Lindzen (1970). We detail here the calculations that allowed us to obtain it.

In the linear analysis, the wind velocity V, pressure p, density ρ and temperatureT are written as the sum of a spherically symmetrical constant component, subscripted 0, and a time-dependent small perturbation, identified bythe symbol δ. Hence, background equilibrium quantities depend on the radial coordinate only, this later being here the altitude with respect to the planet surface z. Perturbed quantities are functions of the time t, altitude z, colatitude θ, and longitude φ. The constant component of V is ignored, which enforces a solid rotation condition and allows us to simply denote by V the velocity vector of tidal winds.

The tidal response of the atmosphere in the accelerated frame rotating with the planet is governed by the perturbed momentum equation (e.g. Siebert 1961, p. 128) $${\partial }_{t}V+2\Omega \wedge V=-\frac{1}{\rho }\nabla \delta p+\frac{\delta \rho }{\rho }g,$$(B.1)

the equation of mass conservation, $$\frac{\text{D}\rho }{\text{D}t}+{\rho }_0÷V=0,$$(B.2)

the conservation of energy, $$\frac{{\mathcal{R}}_{\text{s}}}{{\Gamma }_1-1}\frac{\text{D}T}{\text{D}t}=\frac{gH}{{\rho }_0}\frac{\text{D}\rho }{\text{D}t}+J,$$(B.3)

and the perfect gas law, $$p=\rho {\mathcal{R}}_{\text{s}}T,$$(B.4)

the symbol ∂_t designating the partial time derivative, ∇ the gradient operator, ∇⋅ the divergence and $\frac{\text{D}}{\text{D}t}\equiv {\partial }_t+V\text{}\cdot \text{}\nabla$ the materialderivative.

We note that we have neglected the force resulting from the tidal gravitational potential in the momentum equation as we focus on the thermally generated component of the tidal response. The only source term is thus the net absorbed heat per unit mass J in the right-hand side of the conservation of energy. Moreover, dissipative processes are not taken into account since they play a role in the low-frequency regime mainly, where their characteristic associated timescales are comparable or greater than the tidal period. The considered regime is thus adiabatic, which is an appropriate approximation for studying the atmospheric tidal response in the high-frequency range, as demonstrated by early studies of the Earth’s case (see Lindzen & Chapman 1969, for a review).

Other approximations are convenient to simplify analytic calculations. First, considering that H ≪ R_p, we assume the hydrostatic approximation, as discussed in Sect. 2. Second, we neglect Coriolis components associated with a vertical displacement. This approximation, known as the traditional approximation (e.g. Unno et al. 1989), is appropriate provided that the fluid layer is stably-stratified or thin with respect to the planet radius, which is the case here. Thus the preceding equations reduce to $$\begin{array}{l}{\partial }_t{V}_{\theta }-2\Omega \cos \theta V_{\phi }=-\frac{1}{R_{\text{p}}}{\partial }_{\theta }\left(\frac{\delta p}{{\rho }_0}\right),\\ {\partial }_t{V}_{\phi }+2\Omega \cos \theta V_{\theta }=-\frac{1}{R_{\text{p}}\sin \theta }{\partial }_{\phi }\left(\frac{\delta p}{{\rho }_0}\right),\\ {\partial }_z\delta p=-g\delta \rho ,\\ {\partial }_t\delta \rho +\frac{\text{d}{\rho }_0}{\text{d}z}{V}_r+{\rho }_0{\partial }_z{V}_r=-{\rho }_0{\nabla }_{\text{h}}\text{}\cdot \text{}V\\ \frac{{\mathcal{R}}_{\text{s}}}{{\Gamma }_1-1}\left({\partial }_t\delta T+\frac{\text{d}{T}_0}{\text{d}z}{V}_r\right)=\frac{gH}{{\rho }_0}\left({\partial }_t\delta \rho +\frac{\text{d}{\rho }_0}{\text{d}z}{V}_r\right)+J,\\ \frac{\delta p}{p_0}=\frac{\delta \rho }{{\rho }_0}+\frac{\delta T}{T_0},\end{array}$$

where the horizontal part of the velocity divergence ∇_h ⋅ V is defined by $${\nabla }_{\text{h}}\text{}\cdot \text{}V\equiv \frac{1}{R_{\text{p}}\sin \theta }\left[{\partial }_{\theta }\left(V_{\theta }\sin \theta \right)+{\partial }_{\phi }{V}_{\phi }\right].$$(B.11)

Beside, we introduce the atmospheric pressure height scale, $$H\equiv \frac{{\mathcal{R}}_{\text{s}}{T}_0}{g},$$(B.12)

and the pressure altitude, $$x\equiv {\displaystyle {\int }_0^z\frac{\text{d}z}{H}},$$(B.13)

which will be used instead of z for convenience in the following.

Since it is periodical in time and longitude, the perturbation can be written as a combination of Fourier modes, each one being associated with a forcing frequency σ and a longitudinal degree m. Fourier coefficients are functions of the colatitude (θ) and altitude (z). Any perturbed quantity q is thus expressed as $$q\left(x,\theta ,\phi \right)={\displaystyle \sum _{\sigma ,m}}{q}^{m,\sigma }\left(z,\theta \right){\text{e}}^{i\left(\sigma t+m\phi \right)},$$(B.14)

where the q^m,σ are the Fourier coefficients of the expansion. Under the assumed approximations, Eqs. (B.5) and (B.6) are decoupled from the radial momentum equation, which allows us to write $V_{\theta }{}^{m,\sigma }$ and $V_{\phi }{}^{m,\sigma }$ as functions of $\delta p_n^{m,\sigma }/{\rho }_0$. By substituting horizontal winds by the obtained expressions in ∇_h ⋅ V (Eq. (B.11)), and introducing the variable $$G=-\frac{1}{{\Gamma }_1{p}_0}\frac{\text{D}p}{\text{D}t},$$(B.15)

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) $${\mathcal{F}}^{m,\sigma }{G}^{m,\sigma }={\mathcal{L}}^{m,\nu }\left[\left(\frac{\text{d\hspace{0.17em}}\ln H}{\text{d}x}+\kappa \right)G^{m,\sigma }-\frac{\kappa J^{m,\sigma }}{{\Gamma }_{1}gH}\right].$$(B.16)

Here, ${\mathcal{F}}^{m,\sigma }$ is an operator depending on the x coordinate only and ${\mathcal{L}}^{m,\nu }$ the Laplace’s tidal operator, which depends on the θ coordinate only and is formulated as (e.g. Lee & Saio 1997) $$\begin{array}{ll}{\mathcal{L}}^{m,\nu }\equiv \hfill & \frac{1}{\sin \theta }{\partial }_{\theta }\left(\frac{\sin \theta }{1-{\nu }^2{\cos }^2\theta }{\partial }_{\theta }\right)\hfill \\ \hfill & -\frac{1}{1-{\nu }^2{\cos }^2\theta }\left(m\nu \frac{1+{\nu }^2{\cos }^2\theta }{1-{\nu }^2{\cos }^2\theta }+\frac{m^2}{{\sin }^2\theta }\right),\hfill \end{array}$$(B.17)

the quantity ν ≡ 2Ω∕σ designating the so-called spin parameter.

The above separation of coordinates allows us to expand the Fourier coefficients of G as $$G^{m,\sigma }\left(x,\theta \right)={\displaystyle \sum _n}{G}_n^{m,\sigma }\left(x\right){\Theta }_n^{m,\nu }\left(\theta \right).$$(B.18)

The Fourier coefficients of δp, δρ and δT are written likewise. In Eq. (B.18), the integer n corresponds to the latitudinal wavenumber of a spherical mode modified by the planet spin rotation (in the static case, where ν = 0, $n=l-\left|m\right|$). Similarly, the Hough functions ${\Theta }_n^{m,\nu }$ correspond to the associated Legendre functions (see Appendix A) modified by rotation. Hough functions are the eigenvectors of the Laplace operator (e.g. Lee & Saio 1997). They are associated with the eigenvalues ${\Lambda }_n^{m,\nu }$ through the relationship $${\mathcal{L}}^{m,\nu }{\Theta }_n^{m,\nu }={\Lambda }_n^{m,\nu }{\Theta }_n^{m,\nu },$$(B.19)

and determine the equivalent depth of the mode associated with the triplet $\left(n,m,\sigma \right)$ (e.g. Taylor 1936), $$h_n^{m,\sigma }\equiv \frac{R_{\text{p}}^2{\sigma }^2}{g{\Lambda }_n^{m,\nu }}.$$(B.20)

In the absence of resonances, the semidiurnal tidal response is generally dominated by the fundamental gravity mode, indicated by n = 0 ³, 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, ${\Lambda }_0^{2,\nu }\approx {\Lambda }_0^{2,1}\approx 11.1$, from the moment that $n_{\star }\ll |\Omega |$.

We note that this value, denoted by Λ₀ 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 ${\Lambda }_0^{2,\nu }\approx 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_n^{m,\sigma }$ is now simply denoted by G, and so on for the tidal heat source, pressure, density, temperature, wind velocity components, eigenvalues and equivalent depths. The usual change of variable G = e^x∕2y leads to the vertical structure equation in its canonical form, $$\frac{{\text{d}}^{2}y}{\text{d}{x}^2}+{\widehat{k}}_x^{2}y=\frac{\kappa J}{{\Gamma }_{1}gh}{\text{e}}^{-x/2},$$(B.21)

where we have introduced the dimensionless vertical wavenumber ${\widehat{k}}_x$, defined by $${\widehat{k}}_x^2\equiv \frac{1}{4}\left[\frac{4}{h}\left(\kappa H+\frac{\text{d}H}{\text{d}x}\right)-1\right].$$(B.22)

The vertical structure equation describe the behaviour of a forced harmonic oscillator, and ${\widehat{k}}_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, ${\widehat{k}}_x^2\in ℝ$, and its signdirectly determines the nature of waves across the vertical axis. The condition ${\widehat{k}}_x^2>0$ indicates a propagating mode. Conversely, ${\widehat{k}}_x^2<0$ corresponds to an evanescent mode.