Open Access
Issue
A&A
Volume 624, April 2019
Article Number A17
Number of page(s) 22
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201834685
Published online 08 April 2019

© P. Auclair-Desrotour et al. 2019

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 N2-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 Rp and mass Mp, orbiting its host star, of mass M, circularly. The star–planet distance is denoted a, the mean motion of the system n, and the obliquity of the planet is set to zero. We assume that the planet rotates at the spin angular velocity Ω, which is positive if the spin rotation is along the same direction as the orbital motion, and negative otherwise.

The atmosphere of the planet undergoes both the tidal gravitational and thermal forcings of the host star. Below a certain orbital radius, the planet is sufficiently close to the star to make gravitational forces predominate. Thus, its rotation is driven towards spin–orbit synchronization (Ω = n), which is the unique possible final state of equilibrium for the planet rotation in the absence of obliquity and eccentricity. Conversely, the predominance of the thermal tide enables the existence of asynchronous final rotation states of equilibrium, as showed in the case of Venus (e.g. Gold & Soter 1969; Ingersoll & Dobrovolskis 1978; Dobrovolskis & Ingersoll 1980; Correia & Laskar 2001; Auclair-Desrotour et al. 2017a). As a consequence, we ignore here the action of gravitational forces on the atmosphere. We note however that the action of these forces on the atmospheric tidal bulge will be taken into account to compute the tidal torque, as seen in the following.

The thermal forcing results from the day-night periodic cycle. The atmosphere undergoes heating variations due to the time-varying component of the incoming stellar flux F, which scales as the equilibrium one, F=L/(4πa2)$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, δps(t,θ,φ)=l=1+m=llδps;lm,σYlm(θ,φ)eiσt,\begin{equation*} \delta {p_{\textrm{s}}} \left(t , \theta , \varphi \right) = \sum_{l = 1}^{+ \infty} \sum_{m =- l}^l {{\delta p}}_{\textrm{s};l}^{{m,\sigma}} {Y_{l}^{m}} \left( \theta , \varphi \right) \textrm{e}^{i \sigma t }{,}\end{equation*}(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, Ylm${Y_{l}^{m}} $ the normalized spherical harmonics (see Appendix A), δps;lm,σ$ {{\delta p}}_{\textrm{s};l}^{{m,\sigma}} $ the associated components, and σ=m(Ωn)$\sigma\,{=}\,m \left(\Omega - n_{\star} \right)$ the associated forcing frequencies (see e.g. Efroimsky 2012; Ogilvie 2014)1.

The tidal torque exerted on the atmosphere is obtained by integrating the gravitational force undergone by the tidal bulge over the sphere. Hence, denoting UT the tidal gravitational potential at the planet surface, the atmospheric tidal torque is defined in the thin-layer approximation (HRp) as (e.g. Zahn 1966) T=1gSφUTδps dS,\begin{equation*} {\mathcal{T}} = \frac{1}{g} {\int_{\mathcal{S}}^{} {{\partial_{\varphi} {U_{\textrm{T}}}} \, \delta {p_{\textrm{s}}}} \textrm{d} {\mathcal{S}}},\end{equation*}(2)

where the notation g refers to the surface gravity of the planet, φ to the partial derivative in longitude, S$ \mathcal{S} $ to the sphere of radius Rp, and dS=Rp2sinθdθdφ$\textrm{d} \mathcal{S}\,{=}\,R_{\textrm{p}}^2 \sin \theta {\textrm{d}} \theta {\textrm{d}} \varphi$ to the surface element.

Similarly as the surface pressure anomaly, UT can be expanded in Fourier series of time and spherical harmonics, UT(t,θ,φ)=l=1+m=llUT;lm,σYlm(θ,φ)eiσt,\begin{equation*} U_{\textrm{T}} \left( t , \theta , \varphi \right) = \sum_{l = 1}^{+ \infty} \sum_{m=-l}^l {U}_{\textrm{T};l}^{m,\sigma} {Y_{l}^{m}} \left( \theta , \varphi \right) \textrm{e}^{i \sigma t}, \end{equation*}(3)

where the UT;lm,σ$ {U}_{\textrm{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=|m|=2$ l\,{=}\,{\left| m \right|}\,{=}\,2 $. Besides, since the UT;lm,σ$ {U}_{\textrm{T};l}^{m,\sigma} $ scale as UT;lm,σ(Rp/a)l$ {U}_{\textrm{T};l}^{m,\sigma} \,{\propto}\, \left( R_{\textrm{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 Rp is assumed to be small compared to the star-planet distance, which is the case in the present study.

Thus, by substituting UT and δps by their expansions in spherical harmonics in Eq. (2), we note that the quadrupolar terms l=|m|=2$ l\,{=}\, {\left| m \right|}\,{=}\,2$ only remain, as for the tidal potential UT, 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), T=24π5MMpRp6a3{δps;22,σ},\begin{equation*} {\mathcal{T}} = \sqrt{\frac{24 \pi}{5}} \frac{M_{\star}}{M_{\textrm{p}}} \frac{R_{\textrm{p}}^6}{a^3} {\Im \left\{{{\delta p}}_{\textrm{s};2}^{2,\sigma} \right\}},\end{equation*}(4)

with σ=2(Ωn)$\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, δps;22,σ${{\delta p}}_{\textrm{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 η=sign(σ){δps;22,σ}$\eta\,{=}\,\textrm{sign} \left( \sigma \right) {\Im \left\{{{\delta p}}_{\textrm{s};2}^{2,\sigma} \right\}}$. When η < 0 (η > 0), the atmospheric tidal torque pushes the planet towards (away from) spin-orbit synchronization, |Ωn|${\left| {\Omega - n_{\star} } \right|} $ decays (increases). Positions for which η = 0 correspond to the stable (dη/dσ|eq<0$ \left. \textrm{d} \eta / \textrm{d} \sigma \right|_{\textrm{eq}} < 0 $) or unstable (dη/dσ|eq>0$ \left. \textrm{d} \eta / \textrm{d} \sigma \right|_{\textrm{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 δps;22,σ${{\delta p}}_{\textrm{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 ps = 10 bar and we assume that the planet is located at the Venus-Sun distance, that is aVenus = 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 ps = 10 bar and the semi-major axis varies. Conversely, in the second family, planets have the same orbital radius, a = aVenus, 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 N2-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 CO2 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 N2 for the mean molecular mass of the atmosphere, Matm = 28.0134 g mol−1 (Meija et al. 2016).

For a perfect diatomic gas, the ratio of heat capacities (also called first adiabatic exponent) is Γ1 = 1.4, and it follows that κ=(Γ11)/Γ1=0.285$\kappa\,{=}\,\left( \Gamma_1 - 1 \right) / \Gamma_1\,{=}\,0.285$ (the parameter κ can also be written κ=RGP/(MatmCp)$\kappa\,{=}\,\mathcal{R}_{\textrm{GP}} / \left( M_{\textrm{atm}} C_{p} \right)$, where RGP$\mathcal{R}_{\textrm{GP}}$ and Cp 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 As = 0.2 and thermal inertia Igr = 2000 J m−2 s−1∕2 K−1, which is a typical value for Venus-like soils (see e.g. Lebonnois et al. 2010)2. All of these parameters remain unchanged for the whole study and are summarized in Table 1.

Our simulations are performed with an upgraded version of the LMD GCM specifically developed for the study of extrasolar planets and paleoclimates (see e.g. Wordsworth et al. 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.

Table 1

Main GCM parameters.

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 Pconv corresponding to the convergence timescale necessary to reach a steady cycle. We note that this period has to be specified for each doublet (a,ps)$\left( a , {p_{\textrm{s}}} \right)$. As a first approximation, it depends on the radiative timescale of the deepest layers of the atmosphere τrad, which scalesas τradpsgCp4σSBTe3,\begin{equation*} \tau_{\textrm{rad}} \,{\propto}\, \frac{{p_{\textrm{s}}}}{g} \frac{C_{p}}{4 \sigma_{\textrm{SB}} T_{\textrm{e}}^3},\end{equation*}(5)

where Te 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 Psol=2π/|Ωn|$P_{\textrm{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 δps;22,σ${{\delta p}}_{\textrm{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 = aVenus and ps = 10 bar).

This procedure provides the value of the tidal torque for a given forcing frequency. In practice, the torque is computed over an interval of the normalized frequency ω=(Ωn)/n$\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 Δω(ωsupωinf)/N$\Delta \omega \equiv \left( \omega_{\textrm{sup}} - \omega_{\textrm{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.

thumbnail 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 (ps = 10 bar and a = aVenus).

thumbnail Fig. 2

Surface pressure anomaly created by the thermal tide. Left panels: daily averaged spatial distribution of the departure of the surface pressure from its mean value created by the thermal tide. Right panels: spatial distribution of the semidiurnal component only. The surface pressure anomaly is computed for 300 Solar days and folded over one Solar day centred on the substellar point, whose location and direction of motion are shown with a white arrow. From top to bottom panels: normalized forcing frequency ω=(Ωn)/n$\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 fGCM(σ)$f_{\textrm{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 fodd, defined by fodd(σ)=12[fGCM(σ)fGCM(σ)],\begin{equation*} f_{\textrm{odd}} \left( \sigma \right)\,{=}\,\frac{1}{2} \left[ f_{\textrm{GCM}}\left( \sigma \right) - f_{\textrm{GCM}} \left( -\sigma \right) \right],\end{equation*}(6)

which is the odd function f minimizing for any σ the distance defined by d(σ)=||fGCM(σ)f(σ)||fGCM(σ)f(σ)||.\begin{equation*} d \left( \sigma \right)\,{=}\,\left| \left| f_{\textrm{GCM}} \left( \sigma \right) - f \left( \sigma \right) \right| - \left| f_{\textrm{GCM}} \left( - \sigma \right) - f \left(- \sigma \right) \right| \right|. \end{equation*}(7)

The complementary function feven, such that fGCM = fodd + feven, is defined by feven=12[fGCM(σ)+fGCM(σ)],\begin{equation*} f_{\textrm{even}} = \frac{1}{2} \left[ f_{\textrm{GCM}} \left( \sigma \right) &#x002B; f_{\textrm{GCM}} \left( - \sigma \right) \right],\end{equation*}(8)

and provides a measure of the impact of Coriolis effects on the tidal torque. The data, the interpolating function fGCM and its components fodd and feven are plotted in Fig. 3 as functions of the normalized tidal frequency ω=σ/(2n)$\omega\,{=}\,\sigma / \left( 2 n_{\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 qmaxmax{fodd(σ)}$q_{\textrm{max}} \equiv \max \left\{ f_{\textrm{odd}} \left( \sigma \right) \right\} $ and the associated frequency σmax, such that fodd(σmax)=qmax$f_{\textrm{odd}} \left( \sigma_{\textrm{max}} \right)\,{=}\,q_{\textrm{max}}$, timescale τmaxσmax1$\tau_{\textrm{max}} \equiv \sigma_{\textrm{max}}^{-1}$ and normalized frequency ωmax=σmax/(2n)$\omega_{\textrm{max}}\,{=}\,\sigma_{\textrm{max}} / \left( 2 n_{\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 Tσα${\mathcal{T}} \,{\propto}\, \sigma^{\alpha}$, with α ≈ 0.73. In the high-frequency range (20|ω|300$20 \lesssim {\left| \omega\right|} \leq 300$), it scales as Tσ1$ {\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 fGCM exhibits a slight systematic asymmetry with respect to the synchronization. This feature is obvious in the low-frequency range, where |fGCM(σ)|>|fGCM(σ)|$\left|{ f_{\textrm{GCM}}\left( -\sigma \right)}\right|> \left|{ f_{\textrm{GCM}} \left( \sigma \right) }\right|$, and tends to vanish while |σ|${\left| {\sigma}\right|} $ increases. Particularly, a small departure between fGCM and fodd 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| \neq \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.

thumbnail Fig. 3

Imaginary part of the Y22${Y_{2}^{2}}$ component of surface pressure variations as a function of the normalized forcing frequency ω=(Ωn)/n$\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=RsTsg,\begin{equation*} H = \frac{\mathcal{R}_{\textrm{s}} {T_{\textrm{s}}}}{g}, \end{equation*}(9)

where Rs$\mathcal{R}_{\textrm{s}}$ and Ts 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 = zH. In the analytic model, we choose for the heat per unit mass inducing the tidal response the vertical profile J = JsebJx, where Js is the heat per unit mass at the planet surface, and bJ a dimensionless optical depth corresponding to the inverse of the characteristic thickness of the heated layer. We note that the limit bJ → + 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), σTP=4κHΛ0gRp2.\begin{equation*} \sigma_{\textrm{TP}} = \sqrt{\frac{4 \kappa H \Lambda_0 g}{R_{\textrm{p}}^2}}.\end{equation*}(10)

The notation Λ0 designates here the eigenvalue of the predominating mode in the expansion of perturbed quantities on the basis of Hough functions (see Eq. (B.17)). This mode is the gravity mode of latitudinal wavenumber n = 0 in the indexing notation used by Lee & Saio (1997). Its eigenvalue Λ0 can be approximated as a constant provided that n|Ω|$n_{\star} \ll {\left|{\Omega}\right|}$.

Hence, introducing the equivalent depth of the mode, hRp2σ2Λ0g,\begin{equation*} h \equiv \frac{R_{\textrm{p}}^2 \sigma^2}{\Lambda_0 g},\end{equation*}(11)

we obtain, for |σ|σTP${\left| {\sigma}\right|} \leq \sigma_{\textrm{TP}}$, {δps}=σ1psκJsgHHh(bJ+12+κ)12(bJ+1)[bJ(bJ+1)+κHh](Hh1Γ1),\begin{equation*} {\Im \left\{{\delta {p_{\textrm{s}}}}\right\}}\,{=}\,\sigma^{-1} {p_{\textrm{s}}} \frac{\kappa J_{\textrm{s}}}{g H} \frac{\frac{H}{h} \left( b_{\textrm{J}} &#x002B; \frac{1}{2} &#x002B; \kappa \right) - \frac{1}{2} \left( b_{\textrm{J}} &#x002B; 1 \right) }{\left[ b_{\textrm{J}} \left( b_{\textrm{J}} &#x002B; 1 \right) &#x002B; \frac{\kappa H}{h} \right] \left( \frac{H}{h} - \frac{1}{\Gamma_1} \right) } ,\end{equation*}(12)

and, for |σ|>σTP${\left| {\sigma}\right|} > \sigma_{\textrm{TP}}$, {δps}=σ1psκJsgh[bJ+12(114κHh)][bJ(bJ+1)+κHh][Hh12(1+14κHh)].\begin{equation*} {\Im \left\{{\delta {p_{\textrm{s}}}}\right\}} = \frac{ \sigma^{-1} {p_{\textrm{s}}} \frac{\kappa J_{\textrm{s}}}{g h} \left[ b_{\textrm{J}} &#x002B; \frac{1}{2} \left( 1 - \sqrt{1 - \frac{4 \kappa H}{h}} \right) \right]}{ \left[ b_{\textrm{J}} \left( b_{\textrm{J}} &#x002B; 1 \right) &#x002B; \frac{\kappa H}{h} \right] \left[ \frac{H}{h} - \frac{1}{2} \left( 1 &#x002B; \sqrt{1 - \frac{4 \kappa H}{h}} \right) \right] }.\end{equation*}(13)

We recall that κ=Rs/Cp$\kappa\,{=}\,\mathcal{R}_{\textrm{s}} / C_{p}$, where Cp designates the heat capacity per unit mass, and Γ1=1/(1κ)$\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 hL = Γ1H (Lindzen & Blake 1972). In the asymptotic regime, where n|Ω|$n_{\star} \ll {\left|{\Omega}\right|}$, the characteristic Lamb frequency follows from Eq. (11), σL=Λ0ghLRp2=Γ14κσTP.\begin{equation*} \sigma_{\textrm{L}}\,{=}\,\sqrt{\frac{\Lambda_0 g h_{\textrm{L}}}{R_{\textrm{p}}^2}}\,{=}\,\sqrt{\frac{\Gamma_1}{4 \kappa}} \sigma_{\textrm{TP}}.\end{equation*}(14)

By noticing that σL > σTP in the case of a diatomic gas (Γ1 = 1.4) and substituting h by hL into the corresponding expression of the solution – that is Eq. (13) – we can easily observe that the tidal torque is singular at |σ|=σL${\left| {\sigma}\right|}\,{=}\,\sigma_{\textrm{L}}$. The resonance hence occurs when the phase velocity of the forced mode equalizes the characteristic Lamb velocity VL=ghL$V_{\textrm{L}}\,{=}\,\sqrt{g h_{\textrm{L}}}$.

With the numerical values given by Table 1 and the mean surface temperature computed from GCM simulations (Ts ≈ 316 K), the isothermal approximation leads to H ≈ 10.6 km and hL ≈ 15 km for the reference case. Besides, Λ0 ≈ 11.1 in the adiabatic asymptotic regime of high rotation rates. It thus follows that ωL ≈ 308, and we recover the order of magnitude of the frequency identified in Fig. 3 (top left panel) using GCM simulations (i.e. ωL ≈ 260).

The observed departure between values of ωL obtained in analytical and numerical approaches can be explained by the dependence of the resonance on the atmospheric vertical structure (see e.g. Bretherton 1969; Lindzen & Blake 1972). The analytical value corresponds to the case of an isothermal atmosphere of temperature Ts. In reality, the mean temperature vertical profile is characterized by a strong gradient in the troposphere, the temperature decaying linearly from ~316 K at z = 0 to ~160 K at z ≈ 25 km in GCM simulations. As a consequence, the mean pressure height scale of the tidally heated layer is less than the surface pressure heigh scale, which leads to smaller equivalent depth and resonance frequency for the Lamb mode.

The other interesting feature highlighted by Fig. 3 is the scaling law of the torque Tσ1${\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 |σ|<σTP${\left| {\sigma}\right|} < \sigma_{\textrm{TP}}$, which leads us to consider the solution given by Eq. (12). We place ourselves in the configuration where characteristic timescales are clearly separated, that is |σ|σTP${\left| {\sigma}\right|} \ll \sigma_{\textrm{TP}}$ and n|Ω|$n_{\star} \ll {\left|{\Omega}\right|}$ in the meantime. As Hhσ−2, the preceding condition implies that Hh ≫ 1. It follows that Hh(bJ+12+κ)12(bJ+1)Hh1Γ1~bJ+12+κ.\begin{equation*} \frac{\frac{H}{h} \left( b_{\textrm{J}} &#x002B; \frac{1}{2} &#x002B; \kappa \right) - \frac{1}{2} \left( b_{\textrm{J}} &#x002B; 1 \right) }{ \frac{H}{h} - \frac{1}{\Gamma_1} } \sim b_{\textrm{J}} &#x002B; \frac{1}{2} &#x002B; \kappa. \end{equation*}(15)

By invoking the strong optical thickness of the atmosphere in the infrared (bJ ≫ 1), we remark that we recover analytically the scaling law Tσ1${\mathcal{T}} \,{\propto}\, \sigma^{-1}$ observed in Fig. 3 from the moment that the condition 1H/hκ1bJ2$1 \ll H / h \ll \kappa^{-1} b_{\textrm{J}}^2$ is satisfied. This provides a definition for the intermediate frequency-range, which is now the range corresponding to σJ|σ|σL$\sigma_{\textrm{J}} \ll {\left| {\sigma}\right|} \ll \sigma_{\textrm{L}} $, where we have introduced the thermal frequency σJ=σTP2bJ(bJ+1)σTP.\begin{equation*} \sigma_{\textrm{J}} = \frac{\sigma_{\textrm{TP}}}{2 \sqrt{ b_{\textrm{J}} \left( b_{\textrm{J}} &#x002B; 1 \right)}} \ll \sigma_{\textrm{TP}}. \end{equation*}(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 |σ|σL${\left| {\sigma}\right|} \ll \sigma_{\textrm{L}}$ (or Hh ≫ 1), Eq. (12) can be approximated by the function {δps}2qJτJσ1+(τJσ)2,\begin{equation*} {\Im \left\{{\delta {p_{\textrm{s}}}}\right\}} \approx \frac{2 q_{\textrm{J}} \tau_{\textrm{J}} \sigma }{1 &#x002B; \left( \tau_{\textrm{J}} \sigma \right)^2},\end{equation*}(17)

where the associated characteristic timescale τJ and maximal amplitude of the pressure anomaly qJ are τJ=σJ1,andqJ=κJs(bJ+12+κ)gHσTPbJ(bJ+1)ps. \begin{equation*} \begin{array}{lll} \hspace*{-4.5pt}\displaystyle {\tau_{\textrm{J}} = \sigma_{\textrm{J}}^{-1}} , & \mbox{and} & \displaystyle {q_{\textrm{J}} = \frac{\kappa J_{\textrm{s}} \left( b_{\textrm{J}} &#x002B; \frac{1}{2} &#x002B; \kappa \right)}{ g H \sigma_{\textrm{TP}} \sqrt{b_{\textrm{J}} \left( b_{\textrm{J}} &#x002B; 1 \right)} }{p_{\textrm{s}}} .} \end{array} \end{equation*}(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 σ=2(Ωn)$\sigma\,{=}\,2 \left( \Omega - n_{\star} \right) $ as (e.g. Ingersoll & Dobrovolskis 1978; Auclair-Desrotour et al. 2017a) {δps;22,σ}M=2qMτMσ1+(τMσ)2,\begin{equation*} \Im \left\{ {{\delta p}}_{\textrm{s};2}^{2,\sigma} \right\}_{\textrm{M}}= \frac{2 q_{\textrm{M}} \tau_{\textrm{M}} \sigma}{1 &#x002B; \left( \tau_{\textrm{M}} \sigma \right)^2},\end{equation*}(19)

the notations τM and qM referring to an effective thermal time constant and the amplitude of the maximum (located at σ=τM1$ \sigma\,{=}\,\tau_{\textrm{M}}^{-1}$), respectively (the factor 2 sets the maximal amplitude to qM). This functional form corresponds to the so-called Maxwell model mentioned above. It describes the behaviour of an idealized forced oscillator composed of a string and a damper arranged in series (Greenberg 2009; Efroimsky 2012; Correia et al. 2014).

We note that other works based upon different approaches converged towards the functional form of the Maxwell model. For instance, Correia & Laskar (2001) used the parametrized function f(σ)=σ1(1eγσ2)$ f \left( \sigma \right)\,{=}\,\sigma^{-1} \left(1 - \textrm{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 Tσ1${\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 |σ|σJ${\left| {\sigma}\right|} \ll \sigma_{\textrm{J}}$, the torque derived from our analytic solution in the absence of dissipative mechanisms scales as Tσ${\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 bJ → +).

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)), T=3πn2Rp4εF8CpTsσσ2+τM2,\begin{equation*} {\mathcal{T}} = \frac{3 \pi n_{\star}^2 R_{\textrm{p}}^4 \varepsilon F_{\star}}{8 C_{p} {T_{\textrm{s}}}} \frac{\sigma}{\sigma^2 &#x002B; \tau_{\textrm{M}}^{-2}}, \end{equation*}(20)

where ε stands for the effective fraction of the incoming flux absorbed by the atmosphere. Substituting {δps;22,σ}$\Im \left\{ {{\delta p}}_{\textrm{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 qMτM=312856πGMpεLCpTsRp2a2,\begin{equation*} \frac{q_{\textrm{M}}}{\tau_{\textrm{M}}} = \frac{3}{128} \sqrt{\frac{5}{6 \pi}} \frac{\mathcal{G} M_{\textrm{p}} \varepsilon L_{\star}}{C_{p} {T_{\textrm{s}}} R_{\textrm{p}}^2 a^2},\end{equation*}(21)

the notation G$\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 Ts=[(1As)L16πσSBa2]14,\begin{equation*} {T_{\textrm{s}}} = \left[ \frac{\left( 1 - A_{\textrm{s}} \right) L_{\star} }{16 \pi \sigma_{\textrm{SB}} a^2} \right]^{\frac{1}{4}},\end{equation*}(22)

where we have introduced the Stefan-Boltzmann constant σSB and the surface albedo As. By substituting Ts by Eq. (22) in Eq. (21), we obtain that the ratio qMτM does not depend on the surface pressure and scales as qMτM312856πGMpεLCpRp2[(1As)L16πσSB]14a3/2.\begin{equation*} \frac{q_{\textrm{M}}}{\tau_{\textrm{M}}} \approx \frac{3}{128} \sqrt{\frac{5}{6 \pi}} \frac{\mathcal{G} M_{\textrm{p}} \varepsilon L_{\star}}{C_{p} R_{\textrm{p}}^2 } \left[ \frac{\left( 1 - A_{\textrm{s}} \right) L_{\star} }{16 \pi \sigma_{\textrm{SB}} } \right]^{-\frac{1}{4}} a^{-3/2}.\end{equation*}(23)

with ε = 1 − As if the atmosphere is optically thick in the infrared. This relationship between τM and qM 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 fodd and its associated timescale, respectively. The ab initio analytical solution given by Eqs. (12) and (13) (‘Ana.’) and its Maxwell-like form, derived for |σ|σL${\left| {\sigma}\right|} \ll \sigma_{\textrm{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 bJ (i.e. bJ ≈ 14). Finally, the maximum qJ is obtained by fitting the slope in the intermediate frequency-range to numerical results (qJ ≈ 1042 Pa), and provides the value of the parameter Js (i.e. Js ≈ 0.05 W.kg−1).

Figure 3 highlights the fact that the Maxwell model does not allow us to recover the behaviour of the torque in the low-frequency regime. The functional form given by numerical results and the Maxwell function clearly differ in this regime. Particularly, the maximal amplitude obtained from GCM simulations is about twice larger than that given by the model. We note that a smaller departure between the Maxwell and numerical maxima would certainly be obtained by fitting the Maxwell function to the whole spectrum of numerical results, and not only to the peak. However, this would also lead to overestimate the Maxwell timescale, and the fit would not be satisfactory either. A a consequence, a novel parametrized model has to be introduced to better describe the behaviour of the tidal torque in the low-frequency range. This is the purpose of the next section.

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 Tσ1${\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(ω)0.5$\log \left( \omega \right) \lesssim 0.5$) and the high-frequency asymptotic regime (log(ω)1.5$\log \left( \omega \right) \gtrsim 1.5$). In the interval 0.5log(ω)1.5$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 Fpar(χ)(a1χ+b1)F1(χ)+(a2χ+b2)F2(χ)+btrans[1F1(χ)F2(χ)], \begin{align*}{\mathcal{F}}_{\textrm{par}} \left( \chi \right) \equiv & \left( a_1 \chi &#x002B; b_1 \right) {\mathcal{F}}_{1} \left( \chi \right) &#x002B; \left( a_2 \chi &#x002B; b_2 \right) {\mathcal{F}}_{2} \left( \chi \right) \\ & &#x002B; b_{\textrm{trans}} \left[ 1 - {\mathcal{F}}_{1} \left( \chi \right) - {\mathcal{F}}_{2} \left( \chi \right) \right], \nonumber \end{align*}(24)

where btranslog(qmax)$b_{\textrm{trans}} \approx \log \left( q_{\textrm{max}} \right) $ is the level of the transition plateau, a1, b1, a2 and b2 the dimensionless coefficients of linear functions describing asymptotic regimes, and F1${\mathcal{F}}_{1}$ and F2${\mathcal{F}}_{2}$ two sigmoid activation functions expressed as F1(χ)11+e(χχ1)/d1,andF2(χ)11+e(χχ2)/d2. \begin{equation*} \begin{array}{l l l} \hspace*{-4.5pt}\displaystyle {{\mathcal{F}}_{1} \left( \chi \right) \equiv \frac{1}{1 &#x002B; \textrm{e}^{\left( \chi - \chi_1 \right)/d_1}} } , & \mbox{and} & \displaystyle {{\mathcal{F}}_{2} \left( \chi \right) \equiv \frac{1}{1 &#x002B; \textrm{e}^{- \left( \chi - \chi_2 \right) / d_2}}.} \end{array} \end{equation*}(25)

In these expressions, the dimensionless parameters χ1 and χ2 designate the cutoff frequencies of F1${\mathcal{F}}_{1}$ and F2${\mathcal{F}}_{2}$ in logarithmic scale, and d1 and d2 the widths of transition intervals. The corresponding tidal torque is given by Tpar=10Fpar(log|ω|)sign(ω).\begin{equation*} {\mathcal{T}}_{\textrm{par}} = 10^{{\mathcal{F}}_{\textrm{par}} \left( \log {\left| \omega\right|} \right) } \textrm{sign} \left( \omega \right).\end{equation*}(26)

As the scaling law Tσ1${\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 a2 = − 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 fodd is used). We thus end up with a1=0.734,b1=2.85,d1=0.0100,χ1=0.637,a2=1,b2=4.23,d2=0.0232,χ2=1.20,btrans=3.33, \begin{equation*} \begin{array}{llll} \hspace*{-4.5pt} a_1 = 0.734, & b_1 = 2.85, & d_1 = 0.0100 , & \chi_1 = 0.637 , \\ \hspace*{-4.5pt} a_2 = -1 , & b_2 = 4.23 , & d_2 = 0.0232, & \chi_2 = 1.20, \\ \hspace*{-4.5pt} b_{\textrm{trans}} = 3.33, & & & \end{array} \end{equation*}(27)

and plot the model function Fpar${\mathcal{F}}_{\textrm{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 CO2-dominated atmosphere with a mixture of water and sulphuric acid (H2SO4) 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 N2-dominated atmosphere, plotted in Fig. 3 (right panel). Several interesting features may be noted.

First, the tidal torque exerted on the CO2-dominated atmosphere is more than twice weaker than that exerted on the N2-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 N2-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 feven being not negligible with respectto fodd. This is also an effect of the vertical distribution of tidal heating. In the case of the N2-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 N2 to CO2 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 CO2-dominated atmospheresfor the similarity of configurations they offer with the Venus planet, we choose to focus in this work on N2-dominated atmospheres owing to their simpler frequency-behaviour.

thumbnail Fig. 4

Imaginary part of the Y22${Y_{2}^{2}}$ component of surface pressure variations as a function of the normalized forcing frequency ω=(Ωn)/n$\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 N2-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 Cs 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 = Cppsg + Cs and the emission temperature (Te), they expressed the relationship between surface temperature variations δT and the variations of the incoming stellar flux δFinc as δTs=δFincσMC11+iσ/σM,\begin{equation*} \delta {T_{\textrm{s}}} = \frac{{\delta \! \! \: F}_{\textrm{inc}}}{\sigma_{\textrm{M}} C} \frac{1}{1 &#x002B; i \sigma/ \sigma_{\textrm{M}}},\end{equation*}(28)

with σM=4σSBTe3/C$\sigma_{\textrm{M}}\,{=}\,4 \sigma_{\textrm{SB}} T_{\textrm{e}}^3 / C$ (the subscript M refers to theMaxwell-like form of the function given by Eq. (28)). As we generally observe that TeTs 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 σMsup=4σSBTs3g/(Cpps)$\sigma M^{\textrm{sup}}\,{=}\,4 \sigma_{\textrm{SB}} {T_{\textrm{s}}}^3 g / \left( C_{p} {p_{\textrm{s}}} \right)$. However, in light of typical values of τM obtained with the GCM (see Table 2), it appears that the above formula for σMsup 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, δFinc4σSBTs3δTs+δFatmδQgrδQatm=0,\begin{equation*} {\delta \! \! \: F}_{\textrm{inc}} - 4 \sigma_{\textrm{SB}} {T_{\textrm{s}}}^3 \delta {T_{\textrm{s}}} &#x002B; {\delta \! \! \: F}_{\textrm{atm}} - \delta Q_{\textrm{gr}} - \delta Q_{\textrm{atm}} = 0, \end{equation*}(29)

where we have introduced the small variations of the incoming stellar flux δF, surface temperature δTs, radiative heating by the atmosphere δFatm and diffusive losses in the ground δQgr and in the atmosphere δQatm. Owing to the absence of water, latent heats associated with changes of states are ignored.

In the general case, δTs and δQatm are coupled with the atmospheric tidal response. Particularly, in the Newtonian cooling approximation (i.e. variations of the emitted flux are proportional to temperaturevariations), δFatm can be expressed as δFatm=0+K(x,θ,φ)δT(x,θ,φ) dx,\begin{equation*} {\delta \! \! \: F}_{\textrm{atm}} = {\int_{0}^{&#x002B; \infty} {K \left( x,\theta,\varphi \right) \delta T \left( x,\theta,\varphi \right)} \textrm{d}{x}},\end{equation*}(30)

where K designates an effective coefficient of Newtonian cooling. In order to avoid mathematical complications, we ignore this coupling by assuming either that |δFatm|4σSBTs3|δTs|$\left| {\delta \! \! \: F}_{\textrm{atm}} \right| \ll 4 \sigma_{\textrm{SB}} {T_{\textrm{s}}}^3 \left|{\delta {T_{\textrm{s}}}}\right|$, or, following Bernard (1962), that the variation of the atmospheric flux scales as δFatmδTs in a similar way as the variation of the flux emitted by the ground. This allows us to simplify radiative terms by writing 4σSBTs3δTsδFatm=4σSBTs3ɛsδTs,\begin{equation*} 4 \sigma_{\textrm{SB}} {T_{\textrm{s}}}^3 \delta {T_{\textrm{s}}} - {\delta \! \! \: F}_{\textrm{atm}} = 4 \sigma_{\textrm{SB}} {T_{\textrm{s}}}^3 \epsilon_{\textrm{s}} \delta {T_{\textrm{s}}}, \end{equation*}(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 δTsσ=BsσδFincσ$\delta {T_{\textrm{s}}}^{\sigma}\,{=}\,{\mathfrak{B}_{\textrm{s}}}^{\sigma} {\delta \! \! \: F}_{\textrm{inc}}^{\sigma}$. We thus end up with (see detailed calculations in Appendix D) Bsσ=Bs01+[1+sign(σ)i]τs|σ|,\begin{equation*} {\mathfrak{B}_{\textrm{s}}}^{\sigma} = \frac{{\mathfrak{B}_{\textrm{s}}}^{0}}{1 &#x002B; \left[ 1 &#x002B; \textrm{sign} \left( \sigma \right) i \right] \sqrt{\tau_{\textrm{s}} {\left| {\sigma}\right|}}},\end{equation*}(32)

where Bs0=(4σSBTs3ɛs)1${\mathfrak{B}_{\textrm{s}}}^{0}\,{=}\,\left( 4 \sigma_{\textrm{SB}} {T_{\textrm{s}}}^3 \epsilon_{\textrm{s}} \right)^{-1}$, and τs designates the characteristic timescale of the surface thermal response, which depends on the thermal inertia of the ground Igr and of the atmosphere Iatm at the interface, and is expressed as τs=12(Igr+Iatm4σSBTs3ɛs)2.\begin{equation*} \tau_{\textrm{s}} = \frac{1}{2} \left( \frac{I_{\textrm{gr}} &#x002B; I_{\textrm{atm}}}{4 \sigma_{\textrm{SB}} {T_{\textrm{s}}}^3 \epsilon_{\textrm{s}}} \right)^2.\end{equation*}(33)

We compare this model to numerical results by extracting the Y22${Y_{2}^{2}}$ component of the surface temperature distribution δTs;22,σ${\delta T}_{\textrm{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 Bsσ/Bs0${\mathfrak{B}_{\textrm{s}}}^{\sigma} / {\mathfrak{B}_{\textrm{s}}}^{0}$ (such that δTs;22,σ=BsσδFinc;22,σ${\delta T}_{\textrm{s};2}^{2,\sigma}\,{=}\,{\mathfrak{B}_{\textrm{s}}}^{\sigma} {{\delta \! \! \: F}}_{\textrm{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 Iatm), 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 Bsσ${\mathfrak{B}_{\textrm{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 δFatm (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 στs1$ \sigma \gg \tau_{\textrm{s}}^{-1}$, the model predicts that the amplitude of temperature variations should tend to zero. Yet, we observe that |δTs;22,σ|$\left|{{\delta T}_{\textrm{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).

Table 2

Scaling laws of τmax and qmax obtained using the LMDZ GCM for a dry terrestrial planet with a homogeneous N2 atmosphere.

thumbnail Fig. 5

Nyquist plot of the transfer function Bsσ${\mathfrak{B}_{\textrm{s}}}^{\sigma}$ associated with the Y22${Y_{2}^{2}}$ component of the semidiurnal tide, that is such that δTs;22,σ=BsσδFinc;22,σ${\delta T}_{\textrm{s};2}^{2,\sigma}\,{=}\, {\mathfrak{B}_{\textrm{s}}}^{\sigma} {{\delta \! \! \: F}}_{\textrm{inc};2}^{2,\sigma} $. The imaginary part of the normalized function Bsσ/|Bs0|${\mathfrak{B}_{\textrm{s}}}^{\sigma}/ \left|{{\mathfrak{B}_{\textrm{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), ω=(Ωn)/n$\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 (ps).

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 ps = 10 bar and we compute frequency spectra of the imaginary part of the Y22${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 = aVenus, that is a = 0.723 au, and frequency spectra are computed for ps varying from 1 to 30 bar. The reference case characterized in the previous section, and parametrized by a = aVenus and ps = 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 ps > 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, Pconv. 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 Pconv on the planet surface pressure into account in study 2, we set Pconv to 1100, 2300, 5800 and 14000 Earth Solar days for ps = 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 = aVenus and ps = 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 qmax and σmax with a and ps 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 fGCM(σ)σ1/2$f_{\textrm{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 {δTs}σ1/2${\Im \left\{ \delta {T_{\textrm{s}}} \right\}}\,{\propto}\, \sigma^{1/2} $ in the zero-frequency limit, which is precisely the dependence observed in Fig. 6.

thumbnail Fig. 6

Imaginary part of the Y22${Y_{2}^{2}}$ harmonic of surface pressure oscillations (Pa) as a function of the normalized tidal frequency ω=(Ωn)/n$ \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.30.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 qmax and associated timescale τmax – vary with a and ps.

Thus, for each study, we fit numerical values of qmax and τmax using a linear regression, formulated as Y=αX+β,\begin{equation*} Y = \alpha X &#x002B; \beta, \end{equation*}(34)

where Y designates the logarithm of qmax (Pa) or τmax (days), X the logarithm of a (au) or ps (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 R2. We also compute log(qmax/τmax)$\log \left( q_{\textrm{max}} / \tau_{\textrm{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 qmax, 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 qmax 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 qmax is less than the variability of τmax.

Hence, the linear regression fits particularly well the dependence of qmax 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 τmaxn1a3/2$\tau_{\textrm{max}} \,{\propto}\, n_{\star}^{-1} \,{\propto}\, a^{3/2}$.

As regards the ratio qmaxτ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(qMτM)=32log(a)+2.49,\begin{equation*} \log \left( \frac{q_{\textrm{M}}}{\tau_{\textrm{M}}} \right) = -\frac{3}{2} \log \left( a \right) &#x002B; 2.49,\end{equation*}(35)

if we assume that ε = 1 − As (i.e. the flux reemitted by the ground is entirely absorbed by the atmosphere).

As may be seen, the dependence of qmaxτmax on the surface pressure is small (α = 0.13) for want of being zero, as predicted by the model. Regarding the dependence on a, the relative difference between the numerical and theoretical values of α (i.e. 1.55 and 3∕2, respectively) is around 3%. However, the value of β computed from GCM simulations (2.77) is higher than that predicted by the model (2.49), despite the fact that this latter is an upper estimation. This difference illustrates the limitations of the Maxwell model, which fails to describe the sharp variations of the tidal torque with the tidal frequency when |σ|~σmax$ {\left| {\sigma}\right|} \sim \sigma_{\textrm{max}}$.

thumbnail 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 ps. 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 ps, is denoted by 𝔭, and the considered case is subscripted j. A given family is thus composed of N𝔭 couples of numerical vectors (σj,Tj)$\left( \sigma_{j} , {\mathcal{T}}_{j} \right)$, with 1 ≤ jN𝔭 (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 fj(σ)f(σj,Tj,σ).\begin{equation*} f_{j} \left( \sigma \right) \equiv f \left( \sigma_{j} , {\mathcal{T}}_{j} , \sigma \right) . \end{equation*}(36)

We also introduce the renormalized vectors σ^j$ {\hat{\sigma}}_{j} $ and T^j$ {\hat{{\mathcal{T}}}}_{j} $, defined by σ^jpjα1σj,andT^jpjα2Tj. \begin{equation*} \begin{array}{lll} \displaystyle {{\hat{\sigma}}_{j} \equiv {\mathfrak{p}}_{j}^{-\alpha_1} \sigma_{j} } , & \mbox{and} & \displaystyle {{\hat{{\mathcal{T}}}}_{j} \equiv {\mathfrak{p}}_{j}^{-\alpha_2} {\mathcal{T}}_{j} .} \end{array} \end{equation*}(37)

where α1 and α2 are the exponents characterizing the renormalization.

The size of frequency domains covered by the σ^j${\hat{\sigma}}_{j}$ vectors varies with 𝔭 in the general case. As a consequence, rescaling axes requires to define the bounds of the largest common interval, σ^inf(α1)=max{pjα1σj,1}j,σ^sup(α1)=min{pjα1σj,Nσ}j,\begin{align}& {\hat{\sigma}}_{\textrm{inf}} \left( \alpha_1 \right) = \max \left\{ {\mathfrak{p}}_{j}^{-\alpha_1} \sigma_{j,1} \right\}_{j}, \\ & {\hat{\sigma}}_{\textrm{sup}} \left( \alpha_1 \right) = \min \left\{ {\mathfrak{p}}_{j}^{-\alpha_1} \sigma_{j,{N_{\sigma}}} \right\}_{j}, \end{align}

the notation Nσ referring to the size of the frequency sampling (typically Nσ = 21, see Tables E.1 and E.2), and σj,1 and σj,Nσ$\sigma_{j,{N_{\sigma}}}$ to the lower and upper bounds of the interval sampled by σj, respectively.

The values of α1 and α2 are obtained by minimizing the squared difference function F(α1,α2)=j<kσ^inf(α1)σ^sup(α1)[pjα2fj(pjα1σ^)pkα2fk(pkα1σ^)]2 dσ^σ^sup(α1)σ^inf(α1).\begin{equation*} F \left( \alpha_1 , \alpha_2 \right) = \frac{ \displaystyle \sum_{j < k} {\int_{{\hat{\sigma}}_{\textrm{inf}} \left( \alpha_1 \right)}^{{\hat{\sigma}}_{\textrm{sup}} \left( \alpha_1 \right)} { \! \! \left[ {\mathfrak{p}}_{j}^{-\alpha_2} f_{j} \left( {\mathfrak{p}}_{j}^{\alpha_1} {\hat{\sigma}} \right) - {\mathfrak{p}}_{k}^{-\alpha_2} f_{k} \left( {\mathfrak{p}}_{k}^{\alpha_1} {\hat{\sigma}} \right) \right]^2 \! \!} \textrm{d} {{\hat{\sigma}}} } }{{\hat{\sigma}}_{\textrm{sup}} \left( \alpha_1 \right) - {\hat{\sigma}}_{\textrm{inf}} \left( \alpha_1 \right)} . \end{equation*}(40)

We note that the parameters derived from these calculations, q0 and τ0, slightly differ from qmax and τmax. They stand for the characteristic amplitude and timescale of the peak and not for its maximum value and corresponding forcing period, as qmax and τmax. These parameters are defined by functions of 𝔭 (that is a or ps), by τ0τref(ppref)α1,andq0qref(ppref)α2, \begin{equation*} \begin{array}{l l l} \hspace*{-5pt} \displaystyle {\tau_{0} \equiv \tau_{\textrm{ref}} \left( \frac{{\mathfrak{p}}}{{\mathfrak{p}}_{\textrm{ref}}} \right)^{-\alpha_1}} , & \mbox{and} & \displaystyle {q_{0} \equiv q_{\textrm{ref}} \left( \frac{{\mathfrak{p}}}{{\mathfrak{p}}_{\textrm{ref}}} \right)^{\alpha_2},} \end{array} \end{equation*}(41)

where 𝔭ref designates the value of 𝔭 in the reference case (typically aVenus or 10 bar), qref = 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 (ps = 10 bar and variable a) with the scaling laws log(q0)=0.61log(a)+3.26,log(τ0)=0.68log(a)+0.43,\begin{align}\hspace*{-2pt} & \log \left( q_{0} \right) = - 0.61 \log \left( a \right) &#x002B; 3.26, \\ \hspace*{-2pt} & \log \left( \tau_{0} \right) = 0.68 \log \left( a \right) &#x002B; 0.43, \end{align}

and for the second family (a = aVenus and variableps), with log(q0)=0.48log(ps)+2.87,log(τ0)=0.30log(ps)+0.038.\begin{align}\hspace*{-2pt} & \log \left( q_{0} \right) = 0.48 \log \left( {p_{\textrm{s}}} \right) &#x002B; 2.87, \\\hspace*{-2pt} & \log \left( \tau_{0} \right) = 0.30 \log \left( {p_{\textrm{s}}} \right) &#x002B; 0.038. \end{align}

Let us remind here the used units: a is given in au, ps in bar, q0 in Pa, and τ0 in days.

The last step consists in combining these scaling laws with the parametrized model derived in the reference case (Eq. (24)). Proceeding to the change of variables associated with the renormalization, we obtain the generic parametrized model function {δps;22,σ}parq010Fpar(log|τ0σ|)sign(σ).\begin{equation*} \Im \left\{ {{\delta p}}_{\textrm{s};2}^{2,\sigma} \right\}_{\textrm{par}} \equiv q_{0} 10^{{\mathcal{F}}_{\textrm{par}} \left( \log \left|{\tau_{0} \sigma}\right| \right) } \textrm{sign} \left( \sigma \right).\end{equation*}(46)

We remind here that Fpar${\mathcal{F}}_{\textrm{par}}$ is the function Fpar(χ)(a1χ+b1)F1(χ)+(a2χ+b2)F2(χ)+btrans[1F1(χ)F2(χ)], \begin{align*}{\mathcal{F}}_{\textrm{par}} \left( \chi \right) \equiv & \left( a_1 \chi &#x002B; b_1 \right) {\mathcal{F}}_{1} \left( \chi \right) &#x002B; \left( a_2 \chi &#x002B; b_2 \right) {\mathcal{F}}_{2} \left( \chi \right) \\ & &#x002B; b_{\textrm{trans}} \left[ 1 - {\mathcal{F}}_{1} \left( \chi \right) - {\mathcal{F}}_{2} \left( \chi \right) \right], \nonumber \end{align*}(47)

where F1${\mathcal{F}}_{1}$ and F2${\mathcal{F}}_{2}$ are the activation functions defined by F1(χ)11+e(χχ1)/d1,andF2(χ)11+e(χχ2)/d2. \begin{equation*} \begin{array}{l l l} \hspace*{-4.5pt}\displaystyle {{\mathcal{F}}_{1} \left( \chi \right) \equiv \frac{1}{1 &#x002B; \textrm{e}^{\left( \chi - \chi_1 \right)/d_1}} } , & \mbox{and} & \displaystyle {{\mathcal{F}}_{2} \left( \chi \right) \equiv \frac{1}{1 &#x002B; \textrm{e}^{- \left( \chi - \chi_2 \right) / d_2}}.} \end{array} \end{equation*}(48)

The parameters characterizing the generic formula given by Eq. (46) take the values a1=0.734,b1=0.171,d1=0.010,χ1=0.277,a2=1,b2=0.031,d2=0.023,χ2=0.290,btrans=0.020. \begin{equation*} \begin{array}{llll} \hspace*{-4.5pt}a_1 = 0.734, & b_1 = 0.171, & d_1 = 0.010 , & \chi_1 = -0.277 , \\ \hspace*{-4.5pt}a_2 = -1 , & b_2 = -0.031, & d_2 = 0.023, & \chi_2 = 0.290, \\ \hspace*{-4.5pt} b_{\textrm{trans}} = -0.020. & & &\vspace*{-12pt} \end{array} \end{equation*}(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 τ0 σ (dashed black line). Figure 8 clearly shows the relevance of the rescaling as regards the first family, where the dependence of the torque on the star–planet distance is investigated. After rescaling, spectra look similar and the model matches them fairly well. As regards the second family, we observe a greater variability of q0 and τ0 with a net separation between the reference and 30 bar cases. However, the frequency behaviour of the torque does not change much from one case to another and the parametrized function given by Eq. (46) remains a reasonable approximation of its main features.

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 Y22$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 (ps = 10 bar and a = aVenus), 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 qmaxτmaxa−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.

thumbnail Fig. 8

Normalized imaginary part of the Y22${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.30.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$l \in \mathbb{N}$ and − lml) are defined by (e.g. Arfken & Weber 2005) Ylm(θ,φ)(1)m2l+14π(lm)!(l+m)!Plm(cosθ)eimφ,\begin{equation*} {Y_{l}^{m}} \left( \theta, \varphi \right) \equiv \left( - 1 \right)^m \sqrt{\frac{2 l &#x002B;1}{4 \pi} \frac{\left( l - m \right) !}{\left( l &#x002B; m \right) !}} {P_{l}^{m}} \left( \cos \theta \right) \textrm{e}^{i m \varphi}, \end{equation*}(A.1)

where the Plm${P_{l}^{m}}$ designate the associated Legendre functions, expressed, for x[1,1]$x \in \left[ -1,1 \right]$, as Plm(x)(1)m(1x2)m/2dmdxmPl(x),\begin{equation*} {P_{l}^{m}} \left( x \right) \equiv \left( - 1 \right)^m \left( 1 - x^2 \right)^{m/2} {\dfrac{\textrm{d}^{m}}{\textrm{d} x^{m}}} P_{l} \left( x \right), \end{equation*}(A.2)

the Pl being the Legendre polynomials, defined by Pl(x)12ll!d1dx1[(x21)l].\begin{equation*} P_{l} \left( x \right) \equiv \frac{1}{2^l l !} {\dfrac{\textrm{d}^{1}}{\textrm{d} x^{1}}} \left[ \left( x^2 - 1 \right)^l \right]. \end{equation*}(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) tV+2ΩV=1ρδp+δρρg,\begin{equation*} {\partial_{t} \vec{V}} &#x002B; 2 \vec{\Omega} \wedge \vec{V} = - \frac{1}{\rho} \nabla {\delta p} &#x002B; \frac{{\delta \rho}}{\rho} \vec{g}, \end{equation*}(B.1)

the equation of mass conservation, DρDt+ρ0÷V=0,\begin{equation*} {\dfrac{\textrm{D} \rho}{\textrm{D} t}} &#x002B; \rho_{0} \div \vec{V} = 0, \end{equation*}(B.2)

the conservation of energy, RsΓ11DTDt=gHρ0DρDt+J,\begin{equation*} \frac{\mathcal{R}_{\textrm{s}}}{\Gamma_1 - 1} {\dfrac{\textrm{D} T}{\textrm{D} t}} = \frac{g H}{\rho_{0}} {\dfrac{\textrm{D} \rho}{\textrm{D} t}} &#x002B; J, \end{equation*}(B.3)

and the perfect gas law, p=ρRsT,\begin{equation*} p = \rho \mathcal{R}_{\textrm{s}} T , \end{equation*}(B.4)

the symbol t designating the partial time derivative, ∇ the gradient operator, ∇⋅ the divergence and DDtt+V${\dfrac{\textrm{D}}{\textrm{D} t}} \equiv {\partial_{t}} &#x002B; \vec{V} \! \cdot \! \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 HRp, 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 tVθ2ΩcosθVφ=1Rpθ(δpρ0),tVφ+2ΩcosθVθ=1Rpsinθφ(δpρ0),zδp=gδρ,tδρ+dρ0dzVr+ρ0zVr=ρ0hVRsΓ11(tδT+dT0dzVr)=gHρ0(tδρ+dρ0dzVr)+J,δpp0=δρρ0+δTT0,\begin{align}&{\partial_{t}{V_{\theta}}} - 2 \Omega \cos \theta {V_{\varphi}} = - \frac{1}{R_{\textrm{p}}} {\partial_{\theta}} \left( \frac{{\delta p}}{\rho_{0}} \right), \\&{\partial_{t} {V_{\varphi}}} &#x002B; 2 \Omega \cos \theta {V_{\theta}} = - \frac{1}{R_{\textrm{p}} \sin \theta} {\partial_{\varphi}} \left( \frac{{\delta p}}{\rho_{0}} \right), \\&{\partial_{z} {\delta p}} = - g {\delta \rho}, \\&{\partial_{t} {\delta \rho}} &#x002B; {\dfrac{\textrm{d} \rho_{0}}{\textrm{d} z}} V_r &#x002B; \rho_{0} {\partial_{z} V_r} = - \rho_{0} \nabla_{\textrm{h}} \! \cdot \! \vec{V} \\&\frac{\mathcal{R}_{\textrm{s}}}{\Gamma_1 - 1} \left( {\partial_{t}\delta T} &#x002B; {\dfrac{\textrm{d} T_{0}}{\textrm{d} z}} V_r \right) = \frac{g H}{\rho_{0}} \left({\partial_{t} {\delta \rho}} &#x002B; {\dfrac{\textrm{d} \rho_{0}}{\textrm{d} z}} V_r \right) &#x002B; J, \\&\frac{{\delta p}}{p_{0}} = \frac{{\delta \rho}}{\rho_{0}} &#x002B; \frac{\delta T}{T_{0}}, \end{align}

where the horizontal part of the velocity divergence ∇hV is defined by hV1Rpsinθ[θ(Vθsinθ)+φVφ].\begin{equation*} \nabla_{\textrm{h}} \! \cdot \! \vec{V} \equiv \frac{1}{R_{\textrm{p}} \sin \theta} \left[ {\partial_{\theta}} \left( {V_{\theta}} \sin \theta \right) &#x002B; {\partial_{\varphi} {V_{\varphi}}} \right].\end{equation*}(B.11)

Beside, we introduce the atmospheric pressure height scale, HRsT0g,\begin{equation*} H \equiv \frac{\mathcal{R}_{\textrm{s}} T_{0}}{g } ,\end{equation*}(B.12)

and the pressure altitude, x0zdzH ,\begin{equation*} x \equiv \int_0^z \frac{\textrm{d} z}{H},\end{equation*}(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(x,θ,φ)=σ,mqm,σ(z,θ)ei(σt+mφ),\begin{equation*} q \left( x , \theta , \varphi \right) = \sum_{\sigma,m} q^{m , \sigma} \left( z , \theta \right) \textrm{e}^{i \left( \sigma t &#x002B; m \varphi \right)}, \end{equation*}(B.14)

where the qm,σ 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θm,σ${V_{\theta}}^{{m,\sigma}}$ and Vφm,σ${V_{\varphi}}^{{m,\sigma}}$ as functions of δpnm,σ/ρ0${\delta p}_{n}^{{m,\sigma}}/\rho_{0}$. By substituting horizontal winds by the obtained expressions in ∇hV (Eq. (B.11)), and introducing the variable G=1Γ1p0DpDt,\begin{equation*} G = - \frac{1}{\Gamma_1 p_{0}} {\dfrac{\textrm{D} p}{\textrm{D}t}},\end{equation*}(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) Fm,σGm,σ=Lm,ν[(d lnHdx+κ)Gm,σκJm,σΓ1gH].\begin{equation*} \mathcal{F}^{{m,\sigma}} G^{{m,\sigma}} = \mathcal{L}^{{m,\nu}} \left[ \left( {\dfrac{\textrm{d} \ln H}{\textrm{d} x}} &#x002B; \kappa \right) G^{{m,\sigma}} - \frac{\kappa J^{{m,\sigma}}}{\Gamma_1 g H} \right]. \end{equation*}(B.16)

Here, Fm,σ$\mathcal{F}^{{m,\sigma}}$ is an operator depending on the x coordinate only and Lm,ν$\mathcal{L}^{{m,\nu}}$ the Laplace’s tidal operator, which depends on the θ coordinate only and is formulated as (e.g. Lee & Saio 1997) Lm,ν1sinθθ(sinθ1ν2cos2θθ)11ν2cos2θ(mν1+ν2cos2θ1ν2cos2θ+m2sin2θ), \begin{align*}\mathcal{L}^{{m,\nu}} \equiv & \frac{1}{\sin \theta} {\partial_{\theta}} \left( \frac{\sin \theta}{1 - \nu^2 \cos^2 \theta} {\partial_{\theta}} \right) \\ & - \frac{1}{1 - \nu^2 \cos^2 \theta} \left( m \nu \frac{1 &#x002B; \nu^2 \cos^2 \theta}{1 - \nu^2 \cos^2 \theta} &#x002B; \frac{m^2}{\sin^2 \theta} \right), \nonumber \end{align*}(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 Gm,σ(x,θ)=nGnm,σ(x)Θnm,ν(θ).\begin{equation*} G^{{m,\sigma}} \left( x , \theta \right) = \sum_n G_{n}^{{m,\sigma}} \left( x \right) \Theta_{n}^{m,\nu} \left( \theta \right).\end{equation*}(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|m|$n\,{=}\,l - \left| m \right|$). Similarly, the Hough functions Θnm,ν$\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 Λnm,ν$\Lambda_{n}^{m,\nu}$ through the relationship Lm,νΘnm,ν=Λnm,νΘnm,ν,\begin{equation*} \mathcal{L}^{{m,\nu}} \Theta_{n}^{m,\nu} = \Lambda_{n}^{m,\nu} \Theta_{n}^{m,\nu}, \end{equation*}(B.19)

and determine the equivalent depth of the mode associated with the triplet (n,m,σ)$\left( n , m , \sigma \right)$ (e.g. Taylor 1936), hnm,σRp2σ2gΛnm,ν.\begin{equation*} h_{n}^{{m,\sigma}} \equiv \frac{R_{\textrm{p}}^2 \sigma^2}{g \Lambda_{n}^{m,\nu}}.\end{equation*}(B.20)

In the absence of resonances, the semidiurnal tidal response is generally dominated by the fundamental gravity mode, indicated by n = 0 3, which corresponds to the associated Legendre function P22${P_{2}^{2}}$ in the static case (e.g. Auclair-Desrotour & Leconte 2018). In the high-frequency regime, Λ02,νΛ02,111.1${\Lambda_{0}^{2,\nu}} \approx {\Lambda_{0}^{2,1}} \approx 11.1$, from the moment that n|Ω|$ n_{\star} \ll {\left|{\Omega}\right|} $.

We note that this value, denoted by Λ0 in the following, can be modified by dissipative processes. For instance, by including friction with the planet surface using a Raleigh drag of constant characteristic frequency σR, one may show that the eigenvalue of the modes tends to the value of the static case, that is Λ02,ν6${\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 Gnm,σ$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 = ex∕2y leads to the vertical structure equation in its canonical form, d2ydx2+k^x2y=κJΓ1ghex/2,\begin{equation*} {\dfrac{\textrm{d}^{2} y}{\textrm{d} x^{2}}} &#x002B; \hat{k}_{x}^2 \, y = \frac{\kappa J}{\Gamma_1 g h} \textrm{e}^{-x / 2},\end{equation*}(B.21)

where we have introduced the dimensionless vertical wavenumber k^x$\hat{k}_{x}$, defined by k^x214[4h(κH+dHdx)1].\begin{equation*} \hat{k}_{x}^2 \equiv \frac{1}{4} \left[ \frac{4}{h} \left( \kappa H &#x002B; {\dfrac{\textrm{d} H}{\textrm{d} x}} \right) -1 \right]. \end{equation*}(B.22)

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

Computing analytic solutions turns out to be a very challenging problem except for a few simplified configurations. Therefore, we treat here the idealized case of the isothermal atmosphere, which is one of these configurations. We acknowledge the limitations of this academic atmospheric structure regarding real ones, where convective instability leads to a strong temperaturegradient near the planet surface. However, this approach appears to be sufficient for the purpose of this appendix.

In the isothermal approximation, the temperature profile is supposed to be invariant with the radial coordinate. In light of Eq. (B.12), it immediately follows that H is a constant, dH∕dx = 0, and k^x2=14[4κHh1].\begin{equation*} \hat{k}_{x}^2 = \frac{1}{4} \left[ \frac{4 \kappa H}{h} - 1 \right]. \end{equation*}(B.23)

The above expression shows the existence of a turning point for h = 4κH, where the sign of k^x2$\hat{k}_{x}^2$ changes. This turning point occurs at the frequency σTP=4κHΛ0gRp2.\begin{equation*} \sigma_{\textrm{TP}} = \sqrt{\frac{4 \kappa H \Lambda_0 g}{R_{\textrm{p}}^2}}. \end{equation*}(B.24)

In the reference case of the study, GCM simulations provide Ts ≈ 316 K, which, combined with Rs297 Jkg1K1$\mathcal{R}_{\textrm{s}} \approx 297~{\textrm{J\,kg}^{-1}\,\textrm{K}^{-1}}$, gives H ≈ 10.6 km in the isothermal approximation. An estimation of the normalized frequency ωTP=σTP/(2n)$\omega_{\textrm{TP}}\,{=}\,\sigma_{\textrm{TP}} / \left( 2 n_{\star} \right)$ using Eq. (B.12) thus gives ωTP ≈ 270, showing that the turning point occurs in the high-frequency range and must therefore be taken into account in the calculation of an analytical solution. The condition k^x2>0$\hat{k}_{x}^2 > 0$ (|σ|<σTP${\left| {\sigma}\right|} < \sigma_{\textrm{TP}}$) corresponds to an oscillatory regime, while k^x2<0$\hat{k}_{x}^2 < 0$ (|σ|>σTP${\left| {\sigma}\right|} > \sigma_{\textrm{TP}}$) corresponds to an evanescent one.

To solve the vertical structure equation, we have to choose a vertical profile for the tidal heat power per unit mass J. Following Lindzen et al. (1968), we opt for a profile of the form J=JsebJx,\begin{equation*} J = J_{\textrm{s}} {\textrm{e}}^{- b_{\textrm{J}} x},\end{equation*}(B.25)

where Js stands for the heat absorbed at the planet surface and bJ is a dimensionless optical depth characterizing the decay of heating across the vertical coordinate. This profile is derived from the Beer’s law (e.g. Heng 2017) applied to an isothermal atmosphere of constant extinction coefficient, where the variation of flux is supposed to be proportional to the local pressure, dFF=bJexdx.\begin{equation*} \frac{\textrm{d} F}{F} = - b_{\textrm{J}} {\textrm{e}}^{-x} \textrm{d} x. \vspace*{-1pt}\end{equation*}(B.26)

Integrating this expression from the surface to x and introducing the surface flux Fs, we obtain F=FsebJ(ex1),\begin{equation*} F = F_{\textrm{s}} {\textrm{e}}^{b_{\textrm{J}} \left( \textrm{e}^{- x} - 1 \right)}, \vspace*{-1pt}\end{equation*}(B.27)

which, linearized in the vicinity of the surface (x ≪ 1) and assuming bJ ≫ 1, can be approximated by FFsebJx.\begin{equation*} F \approx F_{\textrm{s}} {\textrm{e}}^{- b_{\textrm{J}} x}. \vspace*{-1pt}\end{equation*}(B.28)

Since the absorbed heat per unit mass is given by J=1Hρ0dFdx,\begin{equation*} J = - \frac{1}{H \rho_{0}} {\dfrac{\textrm{d} F}{\textrm{d} x}}, \vspace*{-1pt}\end{equation*}(B.29)

and ρ0(x)ex$\rho_{0} \left( x \right) \,{\propto}\, \textrm{e}^{- x}$, we retrieve Eq. (B.25), with Js=gbJpsFs.\begin{equation*} J_{\textrm{s}} = \frac{g b_{\textrm{J}}}{{p_{\textrm{s}}}} F_{\textrm{s}}. \vspace*{-1pt} \end{equation*}(B.30)

With the preceding choices and approximations, the general solution of Eq. (B.21) can be written as y=Aeik^xx+Beik^xx+y(0)e(bJ+12)x,\begin{equation*} y = A {\textrm{e}}^{i \hat{k}_{x} x} &#x002B; B {\textrm{e}}^{ - i \hat{k}_{x} x} &#x002B; y^{\left( 0 \right)} {\textrm{e}}^{- \left( b_{\textrm{J}} &#x002B; \frac{1}{2} \right) x}, \end{equation*}(B.31)

the notations A and B designating integrating constants, and y(0)$y^{\left( 0 \right)}$ the constant factor of the particular solution. The integrating constant B can be easily eliminated by setting the appropriate upper boundary condition. First, consider the oscillatory case, where k^x2>0$\hat{k}_{x}^2 > 0$. If the convention k^x>0$\hat{k}_{x} > 0$ is adopted, the first term of the solution is associated with upward propagation of energy and the second term with downward propagation (Wilkes 1949; Lindzen & Chapman 1969).

Since there is no energy source at x = + , we expect that the energy propagation is only upward at the upper boundary. It immediately follows that B = 0. This boundary condition is known as the radiation condition (Lindzen & Chapman 1969). When k^x2<0$\hat{k}_{x}^2 < 0 $, the radiation condition is not appropriate any more since k^x$\hat{k}_{x}$ is now a pure imaginary number. In this case the wave is evanescent, and the first and second term of the solution correspond to decaying and diverging components, respectively, from the moment that the convention {k^x}>0${\Im \left\{ {\hat{k}_{x}} \right\}} > 0$ is adopted. In this case, a non-divergence condition is applied at the upper boundary on the energy flux (e.g. Wilkes 1949; Lindzen et al. 1968), which leads to B = 0 again.

At the lower boundary, the fact that fluid particles cannot penetrate through the planet surface is enforced by a rigid wall condition. The vertical velocity is set to Vr = 0 at x = 0. As the vertical velocity is given by (Lindzen & Chapman 1969) Vr=Γ1hex/2[dydx+(Hh12)y],\begin{equation*} V_r\,{=}\,\Gamma_1 h \textrm{e}^{x / 2} \left[ {\dfrac{\textrm{d} y}{\textrm{d} x}} &#x002B; \left( \frac{H}{h} - \frac{1}{2} \right) y \right], \end{equation*}(B.32)

the rigid wall condition is expressed as dydx+(Hh12)y=0,\begin{equation*} {\dfrac{\textrm{d} y}{\textrm{d} x}} &#x002B; \left( \frac{H}{h} - \frac{1}{2} \right) y = 0 , \end{equation*}(B.33)

at x = 0. Thus the solution is y=κJsΓ1gh(bJ(bJ+1)+κHh)[1+bJHhik^x+Hh12eik^xx+e(bJ+12)x],\begin{equation*} y = \frac{\kappa J_{\textrm{s}}}{\Gamma_1 g h \left( b_{\textrm{J}} \left( b_{\textrm{J}} &#x002B; 1 \right) &#x002B; \frac{\kappa H}{h} \right) } \left[ \frac{1 &#x002B; b_{\textrm{J}} - \frac{H}{h}}{i \hat{k}_{x} &#x002B; \frac{H}{h} - \frac{1}{2} } \textrm{e}^{i \hat{k}_{x} x} &#x002B; \textrm{e}^{- \left( b_{\textrm{J}} &#x002B; \frac{1}{2} \right) x } \right] , \end{equation*}(B.34)

with the conventions k^x>0$\hat{k}_{x} > 0$ if k^x2>0$\hat{k}_{x}^2 > 0$ and {k^x}>0${\Im \left\{{\hat{k}_{x}}\right\}} > 0$ else. We recover here the solution previously obtained by Lindzen et al. (1968) for different values of bJ.

The surface pressure variation is readily deduced from y using Eq. (B.15), δpsps=iΓ1σy(0).\begin{equation*} \frac{\delta {p_{\textrm{s}}}}{{p_{\textrm{s}}}} = i \frac{\Gamma_1}{\sigma} y \left( 0 \right). \end{equation*}(B.35)

Taking the imaginary part of the preceding expression, we finally obtain, for k^x2>0$\hat{k}_{x}^2 >0 $ (i.e. |σ|<σTP${\left| {\sigma}\right|} < \sigma_{\textrm{TP}}$), {δps}=σ1psκJsgHHh(bJ+12+κ)12(bJ+1)[bJ(bJ+1)+κHh](Hh1Γ1),\begin{equation*} {\Im \left\{{\delta {p_{\textrm{s}}}}\right\}} = \sigma^{-1} {p_{\textrm{s}}} \frac{\kappa J_{\textrm{s}}}{g H} \frac{\frac{H}{h} \left( b_{\textrm{J}} &#x002B; \frac{1}{2} &#x002B; \kappa \right) - \frac{1}{2} \left( b_{\textrm{J}} &#x002B; 1 \right) }{\left[ b_{\textrm{J}} \left( b_{\textrm{J}} &#x002B; 1 \right) &#x002B; \frac{\kappa H}{h} \right] \left( \frac{H}{h} - \frac{1}{\Gamma_1} \right) } ,\end{equation*}(B.36)

and, for k^x2<0$\hat{k}_{x}^2 <0 $ (i.e. |σ|>σTP${\left| {\sigma}\right|} > \sigma_{\textrm{TP}}$), {δps}=σ1psκJsgh[bJ+12(114κHh)][bJ(bJ+1)+κHh][Hh12(1+14κHh)].\begin{equation*} {\Im \left\{{\delta {p_{\textrm{s}}}}\right\}} = \frac{ \sigma^{-1} {p_{\textrm{s}}} \frac{\kappa J_{\textrm{s}}}{g h} \left[ b_{\textrm{J}} &#x002B; \frac{1}{2} \left( 1 - \sqrt{1 - \frac{4 \kappa H}{h}} \right) \right]}{ \left[ b_{\textrm{J}} \left( b_{\textrm{J}} &#x002B; 1 \right) &#x002B; \frac{\kappa H}{h} \right] \left[ \frac{H}{h} - \frac{1}{2} \left( 1 &#x002B; \sqrt{1 - \frac{4 \kappa H}{h}} \right) \right] }.\end{equation*}(B.37)

Although they have been obtained using an idealized atmospheric structure, these two expressions of the semidiurnal surface pressure anomaly inform us about the frequency-behaviour of the thermally generated atmospheric tidal torque in the high-frequency range. First consider Eq. (B.36). By assuming that the condition 1H/hκ1bJ2$ 1 \ll H / h \ll \kappa^{-1} b_{\textrm{J}}^2 $ is satisfied, that is for |σ|σTP${\left| {\sigma}\right|} \ll \sigma_{\textrm{TP}}$ and a small thickness of the heated layer, the torque scales as Tσ1${\mathcal{T}} \,{\propto}\, \sigma^{-1}$. This corresponds to what may be observed in Fig. 3 in the interval 30<|ω|<200$30 < {\left| \omega\right|} < 200$. In the zero-frequency limit, h → 0 by scaling as hσ2 if the dependence of the eigenvalue Λ0 on the forcing frequency is ignored (see Eq. (B.20)). It follows that Tσ${\mathcal{T}} \,{\propto}\, \sigma$ in this asymptotic regime. The transition between the two regimes occurs for |σ|σJ${\left| {\sigma}\right|} \approx \sigma_{\textrm{J}} $ with σJ=σTP2bJ(bJ+1)σTP.\begin{equation*} \sigma_{\textrm{J}} = \frac{\sigma_{\textrm{TP}}}{2 \sqrt{b_{\textrm{J}} \left( b_{\textrm{J}} &#x002B; 1 \right)}} \ll \sigma_{\textrm{TP}}. \end{equation*}(B.38)

Since Hh ≫ 1 when |σ|σL${\left| {\sigma}\right|} \ll \sigma_{\textrm{L}} $ (let us remind that σL designates the Lamb frequency given by Eq. (14)), Eq. (B.36) can be put into the form of the Maxwell function, {δps}2qJτJσ1+(τJσ)2,\begin{equation*} {\Im \left\{{\delta {p_{\textrm{s}}}}\right\}} \approx \frac{2 q_{\textrm{J}} \tau_{\textrm{J}} \sigma }{1 &#x002B; \left( \tau_{\textrm{J}} \sigma \right)^2}, \end{equation*}(B.39)

where the associated characteristic timescale τJ and maximal amplitude qJ are τJ=σJ1,andqJ=κJs(bJ+12+κ)gHσTPbJ(bJ+1)ps. \begin{equation*} \begin{array}{l l l} \hspace*{-3pt}\displaystyle {\tau_{\textrm{J}} = \sigma_{\textrm{J}}^{-1}} , & \mbox{and} & \displaystyle {q_{\textrm{J}} = \frac{\kappa J_{\textrm{s}} \left( b_{\textrm{J}} &#x002B; \frac{1}{2} &#x002B; \kappa \right)}{ g H \sigma_{\textrm{TP}} \sqrt{b_{\textrm{J}} \left( b_{\textrm{J}} &#x002B; 1 \right)} } {p_{\textrm{s}}} .} \end{array} \end{equation*}(B.40)

Let us now consider the case where k^x2<0$\hat{k}_{x}^2 < 0$, described by Eq. (B.37). We notice that the surface pressure anomaly diverges for h = Γ1H. This feature was discussed byearly studies (e.g. Lindzen & Blake 1972). It corresponds to the resonance associated with an horizontally propagating acoustic mode of large wavelength known as the Lamb mode (e.g. Lindzen et al. 1968; Bretherton 1969; Lindzen & Blake 1972; Platzman 1988; Unno et al. 1989). In the reference case, with the values previously used to estimate ωTP, we obtain ωL ≈ 308. The resonance can be observed in Fig. 3 for a smaller frequency (ωL≈ 260) because of the departure between the isothermal atmospheric structure and the realistic one computed in GCM simulations, as discussed in Sect. 4.2.

Appendix C Dependence of the tidal torque on the atmospheric composition

The atmospheric torque generated by the thermal tide depends on the atmospheric composition, which has a strong impact on the vertical distribution of the tidal heating through clouds formation and the optical thickness of the gas mixture. In the study, we treat the case of a terrestrial planet hosting a cloudless N2-dominated atmosphere with a small amount of CO2. Hence, we ignore the effects of clouds and compute the thermal tide of an optically thin atmosphere in the visible frequency range, where the major part of the stellar flux reaches the planet surface without being absorbed.

In this appendix, we consider the case of a planet hosting a Venus-like CO2-dominated atmosphere with a mixture of water and sulphuric acid (H2SO4) in the same reference configuration (a = aVenus and ps = 10 bar). We do not attempt to reproduce exactly the composition and dynamics of the Venus atmosphere, which is a complex problem beyond the scope of this study (see e.g. Lebonnois et al. 2010, 2016), but to simply retrieve its main features (optical opacity, clouds absorption, etc.). As a consequence, we opt for a generic approach excluding a fine tuning of the atmospheric properties. We set the thermal capacity per unit mass of the gas (Cp) to 1000J kg−1 K−1, which is the typical value of Cp in the case of Venus (e.g. Seiff et al. 1985), where the parameter decreases from 1181J kg−1 K−1 near the surface to 904J kg−1 K−1 at an altitude of 50 km.

Similarly, the mean molecular mass is set to 43.45 g mol−1, and the volume mixing ratio of water vapour to 20 ppm (Moroz et al. 1979). We set the diameter of water particles to 3 μ m, which is a typical value in the lower cloud (e.g. Knollenberg & Hunten 1980). To take into account the impact of sulphuric acid on the saturation pressure of water vapour pH2O$p_{{\textrm{H}_2 \textrm{O}}}$, we use the prescription given by Gmitro & Vermeulen (1964) for aqueous sulphuric acid (see Eq. (24) of their article). This prescription is written as a function of the local temperature T, as ln(pH2O)=Aln(298T)+BT+C+DT+ET2,\begin{equation*} \ln \left( p_{{\textrm{H}_2 \textrm{O}}} \right) = A \ln \left( \frac{298}{T} \right) &#x002B; \frac{B}{T} &#x002B; C &#x002B; D T &#x002B; E T^2, \end{equation*}(C.1)

where A, B, C, D, and E are empirical constants. The optical properties of the atmosphere used to compute radiative transfers are pre-computed using the HITRAN 2008 database (Rothman et al. 2009) for the Venus atmospheric mixture (instead of the Earth mixture used in the study). The spectrum of the atmospheric tidal torque due to the semidiurnal tide is plotted in Fig. C.1 with the spectrum of the N2 reference case for comparison.

Because of the opacity of the atmosphere in the visible range, the fraction of the incoming stellar flux reaching the planet surface is less in the case of the CO2 atmosphere than in the case of the N2 atmosphere.Particularly, the resonance peak is strongly attenuated. We also observe a greater impact of Coriolis effects, the asymmetry of the tidal torque between negative and positive frequency ranges being more significant.

This difference can be explained by the vertical distribution of tidal heating. As mentioned above, the major part of the stellar flux reaches the planet surface if the atmosphere is composed of N2, which means that the tidal torque is mainly due to density variations occurring in the vicinity of the ground, where friction predominate over Coriolis forces. In the case of the CO2 atmosphere, an important fraction of the incoming energy flux is absorbed at the cloud level. The contribution of this fraction is thus strongly affected by Coriolis effects through zonal mean flows characterizing the equilibrium dynamical state.

thumbnail Fig. C.1

Imaginary part of the surface pressure variation created by the semidiurnal thermal tide as a function of the normalized tidal frequency ω=(Ωn)/n$\omega\,{=}\, \left( \Omega - n_{\star} \right) / n_{\star}$ for a Venus-sized planet of surface pressure ps = 10 bar and orbital radius a = aVenus. Results obtained using GCM simulations are designated by points. Spectra are plotted for two cases: the reference case of the study, where the atmosphere is mainly composed of N2 (black line), and the case treated in Appendix C, where the atmosphere is composed of CO2 with a mixture of water vapour and sulphuric acid.

Appendix D Simplified ab initio analytical model for the ground thermal response

As mentioned in Sect. 4.6, we follow along the line by Bernard (1962) to study the thermal response of the planet surface. In this approach, we consider the surface-atmosphere interface, located at the altitude z = 0, and write the power flux budget for a small perturbation in the framework of a frequency linear analysis. Hence, any quantity q can be expressed as q = qσeiσt, where σ is the forcing frequency introduced in Sect. 2. In the following, we omit the superscript σ and use q in place of qσ, given that we work in the frequency domain.

A variation of the effective incoming stellar flux (i.e. where the reflected component has been removed), denoted δFinc, is absorbed by the planet surface. A fraction δQgr of this power is transmitted to the ground by thermal conduction, and an other fraction, δQatm, is transmitted to the atmosphere through turbulent thermal diffusion. Finally, the increase of surface temperature δTs generated by δFinc induces a radiative emission, δFrad, which is expressed as δFrad=4σSBTs3δTs${\delta \! \! \: F}_{\textrm{rad}}\,{=}\,4 \sigma_{\textrm{SB}} {T_{\textrm{s}}}^3 \delta {T_{\textrm{s}}}$ in the black body approximation (we recall that σSB and Ts are the Stefan-Boltzmann constant and mean surface temperature introduced in Sect. 3.2, respectively). Since the atmosphere is heated by both the incoming stellar flux and the surface thermal forcing, it undergoes a radiative cooling, similarly as the surface. The flux emitted downward to the surface is denoted δFatm. Thus, the power budget of the thermal perturbation at the interface is expressed as δFinc4σSBTs3δTs+δFatmδQgrδQatm=0.\begin{equation*} {\delta \! \! \: F}_{\textrm{inc}} - 4 \sigma_{\textrm{SB}} {T_{\textrm{s}}}^3 \delta {T_{\textrm{s}}} &#x002B; {\delta \! \! \: F}_{\textrm{atm}} -\delta Q_{\textrm{gr}} - \delta Q_{\textrm{atm}} = 0.\end{equation*}(D.1)

To study the surface thermal response without having to consider the full atmospheric tidal response in its whole complexity, it is necessary to ignore the coupling induced by δFatm. This amounts to assuming either that the emission of the atmosphere towards the planet surface is negligible, or that it is proportional to δTs (see Bernard 1962). Thus, introducing the surface effective emissivity ɛs, radiative terms can be reduced to 4σSBTs3δTsδFatm=4σSBTs3ɛsδTs.\begin{equation*} 4 \sigma_{\textrm{SB}} {T_{\textrm{s}}}^3 \delta {T_{\textrm{s}}} - {\delta \! \! \: F}_{\textrm{atm}} = 4 \sigma_{\textrm{SB}} {T_{\textrm{s}}}^3 \epsilon_{\textrm{s}} \delta {T_{\textrm{s}}}. \end{equation*}(D.2)

The next step consists in defining the thermal exchanges resulting from diffusive processes, δQgr and δQatm. These flux are directly proportional to the gradient of the temperature profile anomaly in the vicinity of the interface, and are expressed as δQgr=kgr(zδT)z=0,δQatm=katm(zδT)z=0+,\begin{equation*} \delta Q_{\textrm{gr}} = k_{\textrm{gr}} \left( {\partial_{z} \delta T} \right)_{z = 0^-},\; \delta Q_{\textrm{atm}} = - k_{\textrm{atm}} \left( {\partial_{z}\delta T} \right)_{z = 0^&#x002B;},\end{equation*}(D.3)

where z designates the partial derivative in altitude, δT the profile of temperature variations, and kgr and katm the thermal conductivities of the ground and of the atmosphere at z = 0, respectively. By introducing the mean density profile of the planet ρ0 and the thermal capacity per unit mass of the ground Cgr (the analogous parameter for the atmosphere being Cp; see Sect. 3.1), the corresponding diffusivities can be defined by Kgrkgrρ0(0)CgrandKatmkatmρ0(0+)Cp. \begin{equation*} \begin{array}{rll} \hspace*{-5pt} K_{\textrm{gr}} \equiv \dfrac{k_{\textrm{gr}}}{\rho_{0} \left( 0^- \right) C_{\textrm{gr}}} & \mbox{and} & K_{\textrm{atm}} \equiv \dfrac{k_{\textrm{atm}}}{\rho_{0} \left( 0^&#x002B; \right) C_{p}}. \end{array} \end{equation*}(D.4)

Temperature variations in the vicinity of the interface are described by the heat transport equation. We assume that diffusive processes predominates in the z → 0 limit. Moreover, since the typical horizontal length scale is far greater than the vertical one in the thin layer approximation, the horizontal component of the Laplacian describing diffusive processes can be neglected with respect to the vertical component in both the solid and atmospheric regions. It follows that iσδT=KgrzzδT,forz0,iσδT=KatmzzδT,forz>0. \begin{align}& i \sigma \delta T = K_{\textrm{gr}} {\partial_{z z} \delta T}, & \textrm{for} \ z \leq 0, \\ & i \sigma \delta T = K_{\textrm{atm}} {\partial_{z z} \delta T}, & \textrm{for} \ z >0. \end{align}

Solving these two equations with constant Kgr and Katm, and ignoring the diverging term in solutions, we end up with δT(z)=δTse[1+sign(σ)i]z/hgrσ,forz0,δT(z)=δTse[1+sign(σ)i]/hatmσ,forz>0, \begin{align}& \delta T \left( z \right) = \delta {T_{\textrm{s}}} {\textrm{e}}^{\left[ 1 &#x002B; \textrm{sign} \left( \sigma \right) i \right] z / h_{\textrm{gr}}^{\sigma}}, & \textrm{for} \ z \leq 0, \\ & \delta T \left( z \right) = \delta {T_{\textrm{s}}} {\textrm{e}}^{- \left[ 1 &#x002B; \textrm{sign} \left( \sigma \right) i \right] / h_{\textrm{atm}}^{\sigma} }, & \textrm{for} \ z > 0,\end{align}

where we have introduced the frequency-dependent skin thicknesses of heat transport by thermal diffusion in the ground (hgrσ$h_{\textrm{gr}}^{\sigma}$) and atmosphere (hatmσ$h_{\textrm{atm}}^{\sigma}$), expressed as hgrσ=2Kgr|σ|andhatmσ=2Katm|σ|. \begin{equation*} \begin{array}{rcl} \hspace*{-5pt} h_{\textrm{gr}}^{\sigma} = \sqrt{ \dfrac{2 K_{\textrm{gr}}}{{\left| {\sigma}\right|}}} & \mbox{and} & h_{\textrm{atm}}^{\sigma} = \sqrt{ \dfrac{2 K_{\textrm{atm}}}{{\left| {\sigma}\right|}}}. \end{array} \end{equation*}(D.9)

The transfer function Bsσ${\mathfrak{B}_{\textrm{s}}}^{\sigma}$ defined by Eq. (32) comes straight forwardly when the profiles of δT are substituted by Eq. (D.8) in Eqs. (D.3) and (D.1) successively. Introducing Bs0=(4σSBTs3ɛs)1${\mathfrak{B}_{\textrm{s}}}^{0}\,{=}\,\left( 4 \sigma_{\textrm{SB}} {T_{\textrm{s}}}^3 \epsilon_{\textrm{s}} \right)^{-1}$, we obtain δTs=BsσδF$\delta {T_{\textrm{s}}}\,{=}\,{\mathfrak{B}_{\textrm{s}}}^{\sigma} {\delta \! \! \: F}$, with Bsσ=Bs01+[1+sign(σ)i]τs|σ|,\begin{equation*} {\mathfrak{B}_{\textrm{s}}}^{\sigma} = \frac{{\mathfrak{B}_{\textrm{s}}}^{0}}{1 &#x002B; \left[ 1 &#x002B; \textrm{sign} \left( \sigma \right) i \right] \sqrt{\tau_{\textrm{s}} {\left| {\sigma}\right|}}},\vspace*{3pt}\end{equation*}(D.10)

where τs can be interpreted as the characteristic timescale of the surface thermal response. The parameter τs is a function of the thermal inertia of the ground Igrρ0(0)CgrKgr$I_{\textrm{gr}} \equiv \rho_{0} \left( 0^- \right) C_{\textrm{gr}} \sqrt{K_{\textrm{gr}}}$ and of the atmosphere Iatmρ0(0+)CpKatm$I_{\textrm{atm}} \equiv \rho_{0} \left( 0^&#x002B; \right) C_{p} \sqrt{K_{\textrm{atm}}} $, τs12(Igr+Iatm4σSBTs3ɛs)2.\begin{equation*} \tau_{\textrm{s}} \equiv \frac{1}{2} \left( \frac{I_{\textrm{gr}} &#x002B; I_{\textrm{atm}}}{4 \sigma_{\textrm{SB}} {T_{\textrm{s}}}^3 \epsilon_{\textrm{s}}} \right)^2. \vspace*{3pt}\end{equation*}(D.11)

The above expression shows that τs compares the efficiency of diffusive processes to that of the radiative cooling of the surface. The thermal time increases with the interface thermal inertia and decays when the surface temperature increases, scaling as τsTs6$\tau_{\textrm{s}} \,{\propto}\, {T_{\textrm{s}}}^{-6}$.

The expression of Bsσ${\mathfrak{B}_{\textrm{s}}}^{\sigma} $ given by Eq. (D.10) highlights two asymptotic regimes. In the low-frequency regime, where |σ|τs1${\left| {\sigma}\right|} \ll \tau_{\textrm{s}}^{-1}$, the surface responds instantaneously to the forcing δF, leading to a surface temperature oscillation in phase with the incoming stellar flux. At σ = 0, the incoming flux is equal to the radiative flux (δFinc=4σSBTs3ɛs${\delta \! \! \: F}_{\textrm{inc}}\,{=}\,4 \sigma_{\textrm{SB}} {T_{\textrm{s}}}^3 \epsilon_{\textrm{s}}$), and Bsσ=Bs0${\mathfrak{B}_{\textrm{s}}}^{\sigma}\,{=}\,{\mathfrak{B}_{\textrm{s}}}^{0}$. In the hight-frequency regime, where |σ|τs1${\left| {\sigma}\right|} \gg \tau_{\textrm{s}}^{-1}$, the amplitude of the surface temperature variations decays and tends to zero in the limit |στs|+$\left|{\sigma \tau_{\textrm{s}}}\right| \rightarrow &#x002B; \infty $. At the transition, that is |σ|=τs1${\left| {\sigma}\right|}= \tau_{\textrm{s}}^{-1} $, Bsσ=(2/5)Bs0$\real{ {\mathfrak{B}_{\textrm{s}}}^{\sigma}}= \left( 2/5 \right) {\mathfrak{B}_{\textrm{s}}}^{0} $ and {Bs}=(1/5)Bs0${\Im \left\{ {{\mathfrak{B}_{\textrm{s}}}}\right\}}\,{=}\,- \left( 1/5 \right) {\mathfrak{B}_{\textrm{s}}}^{0} $.

As discussed in Sect. 4.6, the transfer function obtained using GCM simulations is well approximated by Eq. (D.10) in the low-frequency regime (see Fig. 5). But it diverges from the model when the forcing frequency increases, typically for |σ|τs1${\left| {\sigma}\right|} \gtrsim \tau_{\textrm{s}}^{-1}$. This divergences seems to result mainly from the fact that we ignored the radiative coupling between the atmosphere and the surface associated with δFatm, although it may be very strong. Particularly, this is the case for resonances of the atmospheric tidal response, where δFatm is increased similarly as the amplitude of pressure and temperature oscillations.

Appendix E Tables of values obtained with GCM simulations for the exploration of the parameter space

The values used to plot the frequency-spectra of Fig. 6 are given by Table E.1 for study 1 (dependence on the star–planet distance) and Table E.2 for study 2 (dependence on the planet surface pressure). In both cases, the first column corresponds to the normalized tidal frequency ω=(Ωn)/n$\omega\,{=}\,\left( \Omega - n_{\star} \right) / n_{\star}$.

Table E.1

Values of the imaginary part of the semidiurnal surface pressure anomaly δps;22,σ${{\delta p}}_{\textrm{s};2}^{2,\sigma}$ (Pa) obtained from GCM simulations in study 1, where ps = 10 bar, and used to plot the spectra of Fig. 6 (top panels).

Table E.2

Values of the imaginary part of the semidiurnal surface pressure anomaly δps;22,σ${{\delta p}}_{\textrm{s};2}^{2,\sigma}$ (Pa) obtained from GCM simulations in study 2, where a = aVenus, and used to plot the spectra of Fig. 6 (bottom panels).

References

  1. Anglada-Escudé, G., Amado, P. J., Barnes, J., et al. 2016, Nature, 536, 437 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  2. Arfken,G. B., & Weber, H. J. 2005, Mathematical Methods for Physicists, 6th ed. (Amstradam: Elsevier) [Google Scholar]
  3. Auclair-Desrotour, P., & Leconte, J. 2018, A&A, 613, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Auclair-Desrotour, P., Laskar, J., & Mathis, S. 2017a, A&A, 603, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Auclair-Desrotour, P., Laskar, J., Mathis, S., & Correia, A. C. M. 2017b, A&A, 603, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Auclair-Desrotour, P., Mathis, S., & Laskar, J. 2018, A&A, 609, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Bernard, E. A. 1962, Arch. Meteorol. Geophys. Bioclimatol. A Meteorol. Atmops. Phys., 12, 502 [NASA ADS] [CrossRef] [Google Scholar]
  8. Bretherton, F. P. 1969, QJRMS, 95, 754 [NASA ADS] [CrossRef] [Google Scholar]
  9. Bureau Des Longitudes & Institut de Mécanique Céleste Et de Calcul Des Ephémérides Observatoire de Paris. 2011, Ephémérides Astronomiques 2011: Connaissance des Temps (Paris: EDP Sciences) [Google Scholar]
  10. Chapman, S., & Lindzen, R. 1970, Atmospheric Tides: Thermal and Gravitational (Dordrecht: Reidel) [Google Scholar]
  11. Correia, A. C. M., & Laskar, J. 2001, Nature, 411, 767 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  12. Correia, A. C. M., & Laskar, J. 2003, J. Geophys. Res. Planets, 108, 5123 [NASA ADS] [CrossRef] [Google Scholar]
  13. Correia, A. C. M., Laskar, J., & de Surgy O. N. 2003, Icarus, 163, 1 [NASA ADS] [CrossRef] [Google Scholar]
  14. Correia, A. C. M., Boué, G., Laskar, J., & Rodríguez, A. 2014, A&A, 571, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Dobrovolskis, A. R., & Ingersoll, A. P. 1980, Icarus, 41, 1 [NASA ADS] [CrossRef] [Google Scholar]
  16. Efroimsky, M. 2012, ApJ, 746, 150 [NASA ADS] [CrossRef] [Google Scholar]
  17. Etheridge, D. M., Steele, L. P., Langenfelds, R. L., et al. 1996, J. Geophys. Res., 101, 4115 [NASA ADS] [CrossRef] [Google Scholar]
  18. Forget, F., Wordsworth, R., Millour, E., et al. 2013, Icarus, 222, 81 [NASA ADS] [CrossRef] [Google Scholar]
  19. Gerkema, T., & Zimmerman, J. 2008, An Introduction to Internal Waves, Lecture Notes (Texel: Royal NIOZ), 207 [Google Scholar]
  20. Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Nature, 542, 456 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  21. Gmitro, J. I., & Vermeulen, T. 1964, AIChE J., 10, 740 [CrossRef] [Google Scholar]
  22. Gold, T., & Soter, S. 1969, Icarus, 11, 356 [NASA ADS] [CrossRef] [Google Scholar]
  23. Greenberg, R. 2009, ApJ, 698, L42 [NASA ADS] [CrossRef] [Google Scholar]
  24. Grimm, S. L., Demory, B.-O., Gillon, M., et al. 2018, A&A, 613, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Heng, K. 2017, Exoplanetary Atmospheres: Theoretical Concepts and Foundations (Princeton: Princeton University Press) [Google Scholar]
  26. Heng, K., & Kopparla, P. 2012, ApJ, 754, 60 [Google Scholar]
  27. Heng, K., & Workman, J. 2014, ApJS, 213, 27 [NASA ADS] [CrossRef] [Google Scholar]
  28. Hourdin, F., Musat, I., Bony, S., et al. 2006, Clim. Dyn., 27, 787 [Google Scholar]
  29. Ingersoll, A. P., & Dobrovolskis, A. R. 1978, Nature, 275, 37 [NASA ADS] [CrossRef] [Google Scholar]
  30. Knollenberg, R. G., & Hunten, D. M. 1980, J. Geophys. Res., 85, 8039 [NASA ADS] [CrossRef] [Google Scholar]
  31. Lagage, P.-O. 2015, European Planetary Science Congress, 10, EPSC2015 [Google Scholar]
  32. Lamb, H. 1917, Proc. R. Soc. London Ser. A, 93, 114 [Google Scholar]
  33. Lebonnois, S., Hourdin, F., Eymet, V., et al. 2010, J. Geophys. Res. Planets, 115, E06006 [NASA ADS] [CrossRef] [Google Scholar]
  34. Lebonnois, S., Sugimoto, N., & Gilli, G. 2016, Icarus, 278, 38 [NASA ADS] [CrossRef] [Google Scholar]
  35. Leconte, J., Forget, F., Charnay, B., et al. 2013, A&A, 554, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Leconte, J., Wu, H., Menou, K., & Murray, N. 2015, Science, 347, 632 [NASA ADS] [CrossRef] [Google Scholar]
  37. Lee, U., & Saio, H. 1997, ApJ, 491, 839 [NASA ADS] [CrossRef] [Google Scholar]
  38. Lindzen, R. S., & Blake, D. 1972, J. Geophys. Res., 77, 2166 [NASA ADS] [CrossRef] [Google Scholar]
  39. Lindzen, R. S., & Chapman, S. 1969, Space Sci. Rev., 10, 3 [NASA ADS] [CrossRef] [Google Scholar]
  40. Lindzen, R. S., & McKenzie, D. J. 1967, Pure Appl. Geophys., 66, 90 [Google Scholar]
  41. Lindzen,R. S., Batten, E. S., & Kim, J.-W. 1968, Mon. Weather Rev., 96, 133 [NASA ADS] [CrossRef] [Google Scholar]
  42. Meija, J., Coplen, T. B., Berglund, M., et al. 2016, Pure Appl. Chem., 88, 265 [Google Scholar]
  43. Moroz, V. I., Moshkin, B. E., Ekonomov, A. P., et al. 1979, Sov. Astron. Lett., 5, 118 [NASA ADS] [Google Scholar]
  44. Ogilvie, G. I. 2014, ARA&A, 52, 171 [Google Scholar]
  45. Pierrehumbert, R. T. 2010, Principles of Planetary Climate (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
  46. Pierrehumbert, R. T. 2011, ApJ, 726, L8 [NASA ADS] [CrossRef] [Google Scholar]
  47. Platzman, G. W. 1988, Meteorol. Atmos. Phys., 38, 70 [NASA ADS] [CrossRef] [Google Scholar]
  48. Ribas, I., Bolmont, E., Selsis, F., et al. 2016, A&A, 596, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Rothman, L. S., Gordon, I. E., Barbe, A., et al. 2009, J. Quant. Spectr. Rad. Transf., 110, 533 [NASA ADS] [CrossRef] [Google Scholar]
  50. Seiff, A., Schofield, J. T., Kliore, A. J., et al. 1985, Adv. Space Res., 5, 3 [NASA ADS] [CrossRef] [Google Scholar]
  51. Showman, A. P., & Guillot, T. 2002, A&A, 385, 166 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Showman, A. P., & Polvani, L. M. 2011, ApJ, 738, 71 [NASA ADS] [CrossRef] [Google Scholar]
  53. Siebert, M. 1961, Adv. Geophys., 7, 105 [Google Scholar]
  54. Taylor, G. I. 1936, R. Soc. London Proc. Ser. A, 156, 318 [Google Scholar]
  55. Turbet, M., Bolmont, E., Leconte, J., et al. 2018, A&A, 612, A86 [CrossRef] [EDP Sciences] [Google Scholar]
  56. Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1989, Nonradial Oscillations of Stars (Tokyo: University of Tokyo Press) [Google Scholar]
  57. Vallis, G. K. 2006, Atmospheric and Oceanic Fluid Dynamics (Cambridge: Cambridge University Press), 770 [Google Scholar]
  58. Volland, H. 1974, J. Atmos. Terr. Phys., 36, 445 [NASA ADS] [CrossRef] [Google Scholar]
  59. Wilkes, M. V. 1949, Oscillations of the Earth’s Atmosphere (University Press Cambridge) [Google Scholar]
  60. Wolf, E. T. 2017, ApJ, 839, L1 [Google Scholar]
  61. Wolf, E. T., Shields, A. L., Kopparapu, R. K., Haqq-Misra, J., & Toon, O. B. 2017, ApJ, 837, 107 [NASA ADS] [CrossRef] [Google Scholar]
  62. Wordsworth, R. D., Forget, F., Selsis, F., et al. 2010, A&A, 522, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Wordsworth, R. D., Forget, F., Selsis, F., et al. 2011, ApJ, 733, L48 [NASA ADS] [CrossRef] [Google Scholar]
  64. Wordsworth, R., Forget, F., Millour, E., et al. 2013, Icarus, 222, 1 [NASA ADS] [CrossRef] [Google Scholar]
  65. Zahn, J. P. 1966, Ann. Astrophys., 29, 313 [NASA ADS] [Google Scholar]
  66. Zimbelman, J. R. 1986, Icarus, 68, 366 [NASA ADS] [CrossRef] [Google Scholar]

1

The expression of δps given by Eq. (1) is similar to that of the tidal gravitational potential, which is derived by applying the addition theorem to the components of the potential expanded in Legendre polynomials (see e.g. Efroimsky 2012, Eq. (47)).

2

This is the value prescribed for nonporous basalts by Zimbelman (1986).

3

We follow here the indexing notation by Lee & Saio (1997), which associates g-modes with positive n and r-modes to strictly negative n.

All Tables

Table 1

Main GCM parameters.

Table 2

Scaling laws of τmax and qmax obtained using the LMDZ GCM for a dry terrestrial planet with a homogeneous N2 atmosphere.

Table E.1

Values of the imaginary part of the semidiurnal surface pressure anomaly δps;22,σ${{\delta p}}_{\textrm{s};2}^{2,\sigma}$ (Pa) obtained from GCM simulations in study 1, where ps = 10 bar, and used to plot the spectra of Fig. 6 (top panels).

Table E.2

Values of the imaginary part of the semidiurnal surface pressure anomaly δps;22,σ${{\delta p}}_{\textrm{s};2}^{2,\sigma}$ (Pa) obtained from GCM simulations in study 2, where a = aVenus, and used to plot the spectra of Fig. 6 (bottom panels).

All Figures

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

In the text
thumbnail Fig. 2

Surface pressure anomaly created by the thermal tide. Left panels: daily averaged spatial distribution of the departure of the surface pressure from its mean value created by the thermal tide. Right panels: spatial distribution of the semidiurnal component only. The surface pressure anomaly is computed for 300 Solar days and folded over one Solar day centred on the substellar point, whose location and direction of motion are shown with a white arrow. From top to bottom panels: normalized forcing frequency ω=(Ωn)/n$\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).

In the text
thumbnail Fig. 3

Imaginary part of the Y22${Y_{2}^{2}}$ component of surface pressure variations as a function of the normalized forcing frequency ω=(Ωn)/n$\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).

In the text
thumbnail Fig. 4

Imaginary part of the Y22${Y_{2}^{2}}$ component of surface pressure variations as a function of the normalized forcing frequency ω=(Ωn)/n$\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.

In the text
thumbnail Fig. 5

Nyquist plot of the transfer function Bsσ${\mathfrak{B}_{\textrm{s}}}^{\sigma}$ associated with the Y22${Y_{2}^{2}}$ component of the semidiurnal tide, that is such that δTs;22,σ=BsσδFinc;22,σ${\delta T}_{\textrm{s};2}^{2,\sigma}\,{=}\, {\mathfrak{B}_{\textrm{s}}}^{\sigma} {{\delta \! \! \: F}}_{\textrm{inc};2}^{2,\sigma} $. The imaginary part of the normalized function Bsσ/|Bs0|${\mathfrak{B}_{\textrm{s}}}^{\sigma}/ \left|{{\mathfrak{B}_{\textrm{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), ω=(Ωn)/n$\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.

In the text
thumbnail Fig. 6

Imaginary part of the Y22${Y_{2}^{2}}$ harmonic of surface pressure oscillations (Pa) as a function of the normalized tidal frequency ω=(Ωn)/n$ \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.30.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.

In the text
thumbnail 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).

In the text
thumbnail Fig. 8

Normalized imaginary part of the Y22${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.30.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.

In the text
thumbnail Fig. C.1

Imaginary part of the surface pressure variation created by the semidiurnal thermal tide as a function of the normalized tidal frequency ω=(Ωn)/n$\omega\,{=}\, \left( \Omega - n_{\star} \right) / n_{\star}$ for a Venus-sized planet of surface pressure ps = 10 bar and orbital radius a = aVenus. Results obtained using GCM simulations are designated by points. Spectra are plotted for two cases: the reference case of the study, where the atmosphere is mainly composed of N2 (black line), and the case treated in Appendix C, where the atmosphere is composed of CO2 with a mixture of water vapour and sulphuric acid.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.