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

TL;DR

This paper develops a universal empirical formula for the atmospheric tidal torque of terrestrial planets, based on 3D climate simulations and linear theory, applicable across various orbital distances and atmospheric pressures.

Contribution

It introduces a generic, scalable model for atmospheric tidal torque that captures key resonances and frequency dependencies, improving planetary rotational evolution modeling.

Findings

Tidal torque shows a thermal peak near synchronization and a resonant Lamb mode peak.

The torque spectrum exhibits consistent behavior when scaled by power laws.

A single empirical formula accurately describes the torque across parameter space.

Abstract

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

Tables

Table 1. Table 1 : List of the main GCM parameters. The subscripts ♁ ♁ \Earth and ☉ ☉ \sun refer to the Earth and the Sun, respectively.

Parameter	Symbol	Value	Unit
Planet characteristicsa
Planet mass	$M_p$	0.815	$M_{♁}$
Planet radius	$R_p$	0.950	$R_{♁}$
Surface gravity	$g$	8.87	$m.s^{-2}$
Star characteristicsa
Mass of the star	$M_{\star }$	1.0	$M_{☉}$
Luminosity of the star	$L_{\star }$	1.0	$L_{☉}$
Atmospheric propertiesb
Mean molecular mass	$M_{\mathrm{atm}}$	28.0134	$g.{\mathrm{mol}}^{-1}$
Adiabatic exponent	$\kappa$	0.285	–
Surface parameters
Surface albedo	$A_s$	0.2	–
Surface thermal inertia	$I_{\mathrm{gr}}$	2000	$J.m^{-2}.s^{-1/2}.K^{-1}$

Table 2. Table 2 : Scaling laws of τ max subscript 𝜏 max \tau_{{\rm max}} and q max subscript 𝑞 max q_{{\rm max}} obtained using the LMDZ GCM for a dry terrestrial planet with a homogeneous N 2 subscript N 2 {\rm N_{2}} atmosphere. The scaling law of q max / τ max subscript 𝑞 max subscript 𝜏 max q_{{\rm max}}/\tau_{{\rm max}} is computed using the formers and should be compared to Eq. ( 35 ). Units: a 𝑎 a is given in AU, p s subscript 𝑝 s p_{{\rm s}} in bar, τ max subscript 𝜏 max \tau_{{\rm max}} in days, q max subscript 𝑞 max q_{{\rm max}} in Pa, the parameters of the linear fit α 𝛼 \alpha , β 𝛽 \beta and R 2 superscript 𝑅 2 R^{2} are dimensionless.

Case	Scaling laws of $q_{\max }$ and ${\tau }_{\max }$			$R^2$
	$\log \left(q_{\max }\right)=$	$-0.69\log \left(a\right)$	$+\mathrm{\hspace{0.17em}3.25}$	0.998
$p_s=10$ $\mathrm{bar}$	$\log \left({\tau }_{\max }\right)=$	$0.86\log \left(a\right)$	$+\mathrm{\hspace{0.17em}0.48}$	0.901
	$\log \left(\frac{q_{\max }}{{\tau }_{\max }}\right)=$	$-1.55\log \left(a\right)$	$+\mathrm{\hspace{0.17em}2.77}$	–
	$\log \left(q_{\max }\right)=$	$0.46\log \left(p_s\right)$	$+\mathrm{\hspace{0.17em}2.86}$	0.990
$a=a_{\mathrm{Venus}}$	$\log \left({\tau }_{\max }\right)=$	$0.30\log \left(p_s\right)$	$+\mathrm{\hspace{0.17em}0.13}$	0.959
	$\log \left(\frac{q_{\max }}{{\tau }_{\max }}\right)=$	$0.16\log \left(p_s\right)$	$+\mathrm{\hspace{0.17em}2.73}$	–

Table 3. Table 3 : Values of the imaginary part of the semidiurnal surface pressure anomaly δ ​ p s ; 2 2 , σ 𝛿 superscript subscript 𝑝 s 2 2 𝜎 \delta p_{{\rm s};2}^{2,\sigma} (Pa) obtained from GCM simulations in study 1, where p s = 10 subscript 𝑝 s 10 p_{{\rm s}}=10   bar bar {\rm bar} , and used to plot the spectra of Fig.  6 (top panels). The first column corresponds to the normalized tidal frequency ω = ( Ω − n ⋆ ) / n ⋆ 𝜔 Ω subscript 𝑛 ⋆ subscript 𝑛 ⋆ \omega=\left(\Omega-n_{\star}\right)/n_{\star} .

$\omega$	0.3 AU	0.4 AU	0.5 AU	0.6 AU	0.7 AU	$a_{\mathrm{Venus}}$	0.8 AU	0.9 AU
-30	-473.07	-545.09	-539.98	-533.95	-525.30	-517.29	-471.78	-434.51
-27	-511.43	-585.93	-607.25	-595.19	-512.35	-507.07	-465.63	-482.23
-24	-642.57	-644.54	-685.95	-580.62	-520.41	-513.39	-525.40	-651.49
-21	-704.17	-744.95	-703.56	-582.58	-592.16	-611.03	-772.29	-1035.44
-18	-804.42	-824.14	-666.45	-696.53	-1020.13	-1074.53	-1293.65	-1489.11
-15	-952.30	-778.00	-855.21	-1309.89	-1773.11	-1903.43	-1924.55	-1760.51
-12	-977.99	-1063.36	-1732.51	-2454.48	-2273.65	-2199.45	-2027.54	-1890.20
-9	-1472.86	-2614.20	-2777.14	-2534.91	-2237.35	-2243.82	-2102.46	-1926.23
-6	-3822.28	-3250.69	-2932.75	-2575.46	-2274.10	-2203.40	-2007.00	-1750.41
-3	-3588.82	-2847.24	-2059.37	-1753.74	-1530.60	-1486.34	-1389.68	-1213.63
0	334.42	234.58	151.64	89.75	55.15	49.50	32.02	15.13
3	3374.24	2785.68	2297.94	1964.36	1691.86	1638.27	1482.97	1241.96
6	3654.44	3058.04	2678.35	2433.68	2210.35	2145.61	1960.61	1728.26
9	1529.77	2386.06	2634.64	2356.54	2185.29	2135.41	2020.46	1839.10
12	973.74	1182.24	1720.70	2255.56	2092.64	2043.52	1925.20	1832.17
15	941.31	808.44	932.12	1302.36	1654.98	1732.38	1796.60	1618.30
18	797.69	812.74	694.15	762.96	969.00	1068.94	1276.22	1492.46
21	686.87	737.32	688.13	593.61	640.17	651.72	779.55	989.06
24	627.03	643.14	658.42	603.18	533.36	531.05	558.75	674.29
27	506.85	585.93	603.91	596.34	512.01	504.31	482.52	496.22
30	486.10	545.09	538.44	534.47	524.68	511.96	470.95	445.30

Table 4. Table 4 : Values of the imaginary part of the semidiurnal surface pressure anomaly δ ​ p s ; 2 2 , σ 𝛿 superscript subscript 𝑝 s 2 2 𝜎 \delta p_{{\rm s};2}^{2,\sigma} (Pa) obtained from GCM simulations in study 2, where a = a Venus 𝑎 subscript 𝑎 Venus a=a_{{\rm Venus}} , and used to plot the spectra of Fig.  6 (bottom panels). The first column corresponds to the normalized tidal frequency ω = ( Ω − n ⋆ ) / n ⋆ 𝜔 Ω subscript 𝑛 ⋆ subscript 𝑛 ⋆ \omega=\left(\Omega-n_{\star}\right)/n_{\star} .

$\omega$	1 bar	3 bar	10 bar	30 bar
-30	-354.77	-335.07	-517.29	-541.19
-27	-506.33	-486.33	-507.07	-547.23
-24	-611.39	-741.93	-513.39	-636.13
-21	-673.92	-1022.06	-611.03	-700.94
-18	-664.71	-1213.72	-1074.53	-737.02
-15	-678.58	-1242.91	-1903.43	-767.81
-12	-680.65	-1272.90	-2199.45	-1607.23
-9	-665.48	-1283.57	-2243.82	-2830.42
-6	-495.17	-1004.48	-2203.40	-3252.87
-3	-322.70	-657.66	-1486.34	-2989.54
0	25.37	37.06	49.50	51.20
3	358.98	712.00	1638.27	2989.06
6	504.55	1024.26	2145.61	3020.61
9	610.53	1249.16	2135.41	2598.07
12	679.89	1229.86	2043.52	1640.01
15	660.94	1198.29	1732.38	961.75
18	642.94	1179.86	1068.94	734.32
21	638.42	1003.31	651.72	703.06
24	586.54	726.95	531.05	632.88
27	489.53	488.95	504.31	551.61
30	352.31	371.83	511.96	539.69

Equations

$$\delta p_{\mathrm{s}}\left(t,\theta ,\phi \right)=\sum _{l=1}^{+\infty }\sum _{m=-l}^l\delta p_{\mathrm{s};l}^{m,\sigma }{Y}_l^m\left(\theta ,\phi \right){\mathrm{e}}^{i\sigma t}.$$

$$\mathcal{T}=\frac{1}{g}{\int }_{\mathcal{S}}{\partial }_{\phi }{U}_{\mathrm{T}}\delta p_{\mathrm{s}}\mathrm{d}\mathcal{S},$$

$$U_{\mathrm{T}}\left(t,\theta ,\phi \right)=\sum _{l=1}^{+\infty }\sum _{m=-l}^l{U}_{\mathrm{T};l}^{m,\sigma }{Y}_l^m\left(\theta ,\phi \right){\mathrm{e}}^{i\sigma t},$$

$$\mathcal{T}=\sqrt{\frac{24\pi }{5}}\frac{M_{\star }}{M_{\mathrm{p}}}\frac{R_{\mathrm{p}}^6}{a^3}\Im \left\{\delta p_{\mathrm{s};2}^{2,\sigma }\right\},$$

$${\tau }_{\mathrm{rad}}\propto \frac{p_{\mathrm{s}}}{g}\frac{C_p}{4{\sigma }_{\mathrm{SB}}{T}_{\mathrm{e}}^3},$$

$$f_{\mathrm{odd}}\left(\sigma \right)=\frac{1}{2}\left[f_{\mathrm{GCM}}\left(\sigma \right)-f_{\mathrm{GCM}}\left(-\sigma \right)\right],$$

$$d\left(\sigma \right)=\left|\left|f_{\mathrm{GCM}}\left(\sigma \right)-f\left(\sigma \right)\right|-\left|f_{\mathrm{GCM}}\left(-\sigma \right)-f\left(-\sigma \right)\right|\right|.$$

$$f_{\mathrm{even}}=\frac{1}{2}\left[f_{\mathrm{GCM}}\left(\sigma \right)+f_{\mathrm{GCM}}\left(-\sigma \right)\right],$$

$$H=\frac{{\mathcal{R}}_{\mathrm{s}}{T}_{\mathrm{s}}}{g},$$

$${\sigma }_{\mathrm{TP}}=\sqrt{\frac{4\kappa H{\Lambda }_{0}g}{R_{\mathrm{p}}^2}}.$$

$$h\equiv \frac{R_{\mathrm{p}}^2{\sigma }^2}{{\Lambda }_{0}g},$$

$$\Im \left\{\delta p_{\mathrm{s}}\right\}={\sigma }^{-1}{p}_{\mathrm{s}}\frac{\kappa J_{\mathrm{s}}}{gH}\frac{\frac{H}{h}\left(b_{\mathrm{J}}+\frac{1}{2}+\kappa \right)-\frac{1}{2}\left(b_{\mathrm{J}}+1\right)}{\left[b_{\mathrm{J}}\left(b_{\mathrm{J}}+1\right)+\frac{\kappa H}{h}\right]\left(\frac{H}{h}-\frac{1}{{\Gamma }_1}\right)},$$

$$\Im \left\{\delta p_{\mathrm{s}}\right\}=\frac{{\sigma }^{-1}{p}_{\mathrm{s}}\frac{\kappa J_{\mathrm{s}}}{gh}\left[b_{\mathrm{J}}+\frac{1}{2}\left(1-\sqrt{1-\frac{4\kappa H}{h}}\right)\right]}{\left[b_{\mathrm{J}}\left(b_{\mathrm{J}}+1\right)+\frac{\kappa H}{h}\right]\left[\frac{H}{h}-\frac{1}{2}\left(1+\sqrt{1-\frac{4\kappa H}{h}}\right)\right]}.$$

$${\sigma }_{\mathrm{L}}=\sqrt{\frac{{\Lambda }_{0}g{h}_{\mathrm{L}}}{R_{\mathrm{p}}^2}}=\sqrt{\frac{{\Gamma }_1}{4\kappa }}{\sigma }_{\mathrm{TP}}.$$

$$\frac{\frac{H}{h}\left(b_{\mathrm{J}}+\frac{1}{2}+\kappa \right)-\frac{1}{2}\left(b_{\mathrm{J}}+1\right)}{\frac{H}{h}-\frac{1}{{\Gamma }_1}}\sim b_{\mathrm{J}}+\frac{1}{2}+\kappa .$$

$${\sigma }_{\mathrm{J}}=\frac{{\sigma }_{\mathrm{TP}}}{2\sqrt{b_{\mathrm{J}}\left(b_{\mathrm{J}}+1\right)}}\ll {\sigma }_{\mathrm{TP}}.$$

$$\Im \left\{\delta p_{\mathrm{s}}\right\}\approx \frac{2q_{\mathrm{J}}{\tau }_{\mathrm{J}}\sigma }{1+{\left({\tau }_{\mathrm{J}}\sigma \right)}^2},$$

$$\begin{array}[]{l l l}\displaystyle\tau_{{\rm J}}=\sigma_{{\rm J}}^{-1},&\mbox{and}&\displaystyle q_{{\rm J}}=\frac{\kappa J_{{\rm s}}\left(b_{{\rm J}}+\frac{1}{2}+\kappa\right)}{gH\sigma_{{\rm TP}}\sqrt{b_{{\rm J}}\left(b_{{\rm J}}+1\right)}}p_{{\rm s}}.\end{array}$$

$$\Im {\left\{\delta p_{\mathrm{s};2}^{2,\sigma }\right\}}_{\mathrm{M}}=\frac{2q_{\mathrm{M}}{\tau }_{\mathrm{M}}\sigma }{1+{\left({\tau }_{\mathrm{M}}\sigma \right)}^2},$$

$$\mathcal{T}=\frac{3\pi n_{\star }^2{R}_{\mathrm{p}}^4\epsilon F_{\star }}{8C_p{T}_{\mathrm{s}}}\frac{\sigma }{{\sigma }^2+{\tau }_{\mathrm{M}}^{-2}},$$

$$\frac{q_{\mathrm{M}}}{{\tau }_{\mathrm{M}}}=\frac{3}{128}\sqrt{\frac{5}{6\pi }}\frac{\mathcal{G}{M}_{\mathrm{p}}\epsilon L_{\star }}{C_p{T}_{\mathrm{s}}{R}_{\mathrm{p}}^2{a}^2},$$

$$T_{\mathrm{s}}={\left[\frac{\left(1-A_{\mathrm{s}}\right)L_{\star }}{16\pi {\sigma }_{\mathrm{SB}}{a}^2}\right]}^{\frac{1}{4}},$$

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

$${\mathcal{F}}_{\mathrm{par}}\left(\chi \right)\equiv$$

$$+b_{\mathrm{trans}}\left[1-{\mathcal{F}}_1\left(\chi \right)-{\mathcal{F}}_2\left(\chi \right)\right],$$

$$\begin{array}[]{l l l}\displaystyle\mathcal{F}_{1}\left(\chi\right)\equiv\frac{1}{1+{\rm e}^{\left(\chi-\chi_{1}\right)/d_{1}}},&\mbox{and}&\displaystyle\mathcal{F}_{2}\left(\chi\right)\equiv\frac{1}{1+{\rm e}^{-\left(\chi-\chi_{2}\right)/d_{2}}}.\end{array}$$

$${\mathcal{T}}_{\mathrm{par}}=10^{{\mathcal{F}}_{\mathrm{par}}\left(\log |\omega |\right)}\mathrm{sign}\left(\omega \right).$$

$$\begin{array}[]{llll}a_{1}=0.734,&b_{1}=2.85,&d_{1}=0.0100,&\chi_{1}=0.637,\\ a_{2}=-1,&b_{2}=4.23,&d_{2}=0.0232,&\chi_{2}=1.20,\\ b_{\rm trans}=3.33,&&&\end{array}$$

$$\delta T_{\mathrm{s}}=\frac{\delta F_{\mathrm{inc}}}{{\sigma }_{\mathrm{M}}C}\frac{1}{1+i\sigma /{\sigma }_{\mathrm{M}}},$$

$$\delta F_{\mathrm{inc}}-4{\sigma }_{\mathrm{SB}}{T}_{\mathrm{s}}^3\delta T_{\mathrm{s}}+\delta F_{\mathrm{atm}}-\delta Q_{\mathrm{gr}}-\delta Q_{\mathrm{atm}}=0,$$

$$\delta F_{\mathrm{atm}}={\int }_0^{+\infty }K\left(x,\theta ,\phi \right)\delta T\left(x,\theta ,\phi \right)\mathrm{d}x,$$

Full text

*Aims. *We aim at characterizing the dependence of the tidal torque resulting from the semidiurnal thermal tide on the tidal frequency, the planet orbital radius, and the atmospheric surface pressure.

*Methods. *The tidal torque is computed from full 3D simulations of the atmospheric climate and mean flows using a generic version of the LMDZ general circulation model (GCM) in the case of a nitrogen-dominated atmosphere. Numerical results are discussed with the help of an updated linear analytical framework. Power scaling laws governing the evolution of the torque with the planet orbital radius and surface pressure are derived.

*Results. *The tidal torque exhibits i) a thermal peak in the vicinity of synchronization, ii) a resonant peak associated with the excitation of the Lamb mode in the high frequency range, and iii) well defined frequency slopes outside these resonances. These features are well explained by our linear theory. Whatever the star-planet distance and surface pressure, the torque frequency spectrum – when rescaled with the relevant power laws – always presents the same behaviour. This allows us to provide a single and easily usable empirical formula describing the atmospheric tidal torque over the whole parameter space. With such a formula, the effect of the atmospheric tidal torque can be implemented in evolutionary models of the rotational dynamics of a planet in a computationally efficient, and yet relatively accurate way.

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_{\Earth}$ 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 to constrain 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 neighbors (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 toward 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 tide was computed using general circulation models (or 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 Maxwell model, which describes the forced response of a damped harmonic oscillator. It shows evidence of the important role played by dissipative processes such as radiative cooling in Venus-like configurations. In order 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 atmospheric tidal torque: an ab initio analytical model based upon the linear theory of atmospheric tides (e.g. Chapman & Lindzen, 1970), and a parametrized semi-analytical model derived from results obtained using GCM simulations. This later model describes in a realistic way the behaviour of the torque in the low-frequency range, where a thermal peak is observed. In addition, we investigate in this section the role played by the ground-atmosphere thermal coupling in the lag of the tidal bulge.

In Sect. 5, we examine the dependence of the tidal response on the planet orbital radius and surface pressure. We thus establish empirical scaling laws describing the evolution of the characteristic amplitude and timescale of the thermal peak with these two parameters. Combining together the obtained results, we finally derive a new generic formula to quantify the atmospheric tidal torque created by the thermal semidiurnal tide in the case of a ${\rm N_{2}}$-dominated atmosphere. We give our conclusions in Sect. 6.

2 Basic principle

We briefly recall in this section the main aspects of the mechanism of atmospheric tides in the case of terrestrial planets, and we introduce analytical expressions that will be used in the following to compute the resulting tidal torque. For the sake of simplification, we consider in this study the case of a spherical planet of radius $R_{\rm p}$ and mass $M_{\rm p}$, orbiting its host star, of mass $M_{\star}$, circularly. The star-planet distance is denoted $a$, the mean motion of the system $n_{\star}$, and the obliquity of the planet is set to zero. We assume that the planet rotates at the spin angular velocity $\Omega$, 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 ($\Omega=n_{\star}$), 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. Note however that the action of these forces on the atmospheric tidal bulge will be taken into account to compute the tidal torque, as seen in the following.

The thermal forcing results from the day-night periodic cycle. The atmosphere undergoes heating variations due to the time-varying component of the incoming stellar flux $F$, which scales as the equilibrium one, $F_{\star}=L_{\star}/\left(4\pi a^{2}\right)$, where $L_{\star}$ 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,

[TABLE]

where $l$ and $m$ designate the latitudinal and longitudinal degrees of a mode, $\theta$ and $\varphi$ the colatitude and longitude in the reference frame co-rotating with the planet, $t$ the time, $Y_{l}^{m}$ the normalized spherical harmonics (see Appendix A), $\delta p_{{\rm s};l}^{{m,\sigma}}$ the associated components, and $\sigma=m\left(\Omega-n_{\star}\right)$ the associated forcing frequencies (see e.g. Efroimsky, 2012; Ogilvie, 2014)111The expression of $\delta p_{{\rm s}}$ 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)). .

The tidal torque exerted on the atmosphere is obtained by integrating the gravitational force undergone by the tidal bulge over the sphere. Hence, denoting $U_{\rm T}$ the tidal gravitational potential at the planet surface, the atmospheric tidal torque is defined in the thin-layer approximation ($H\ll R_{\rm p}$) as (e.g. Zahn, 1966)

[TABLE]

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

Similarly as the surface pressure anomaly, $U_{\rm T}$ can be expanded in Fourier series of time and spherical harmonics,

[TABLE]

where the $U_{{\rm 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 center. Thus, the main components of the expansion are those associated with the quadrupolar semidiurnal tide, that is with degrees $l=\left|m\right|=2$. Besides, since the $U_{{\rm T};l}^{m,\sigma}$ scale as $U_{{\rm T};l}^{m,\sigma}\,\propto\,\left(R_{\rm p}/a\right)^{l}$, terms of greater order in $l$ can be neglected with respect to the quadrupolar components if the radius of the planet $R_{\rm p}$ is assumed to be small compared to the star-planet distance, which is the case in the present study.

Thus, by substituting $U_{\rm T}$ and $\delta p_{{\rm s}}$ by their expansions in spherical harmonics in Eq. (2), we note that the quadrupolar terms $l=\left|m\right|=2$ only remain, as for the tidal potential $U_{\rm T}$, and we end up with the well-known expression of the semidiurnal quadrupolar torque in the thin layer approximation (e.g. Leconte et al., 2015),

[TABLE]

with $\sigma=2\left(\Omega-n_{\star}\right)$, and where the notation $\Im$ refers to the imaginary part of a complex number ($\Re$ referring to the real part). In this expression, $\delta p_{{\rm s};2}^{2,\sigma}$ designates the component of degrees $l=2$ and $m=2$ in the expansion on spherical harmonics given by Eq. (1). This complex quantity is the most important one since it encompasses the whole physics of the atmospheric tidal response. In the following, it will be calculated using a GCM.

The action of the torque on the planet is fully determined by the sign of the product $\eta={\rm sign}\left(\sigma\right)\Im\left\{\delta p_{{\rm s};2}^{2,\sigma}\right\}$. When $\eta<0$ ($\eta>0$), the atmospheric tidal torque pushes the planet toward (away from) spin-orbit synchronization, $\left|\Omega-n_{\star}\right|$ decays (increases). Positions for which $\eta=0$ correspond to the stable ($\left.{\rm d}\eta/{\rm d}\sigma\right|_{\rm eq}<0$) or unstable ($\left.{\rm d}\eta/{\rm d}\sigma\right|_{\rm eq}>0$) equilibrium rotation rates that the planet would reach if it were subject to atmospheric tides only, that is if solid tides were ignored in the case of a dry terrestrial planet.

3 Method

As mentioned in the previous section, the guideline of the method is to compute the quadrupolar component of the surface pressure anomaly from 3D GCM simulations. We detail here the basic physical setup of these simulations in a first time, and the way $\delta p_{{\rm s};2}^{2,\sigma}$ is extracted from pressure snapshots in a second time. In the whole study, we focus on a Venus-sized planet orbiting a Sun-like star.

A “reference case” of fixed surface pressure and star-planet distance is defined. Specifically, the surface pressure is set in this case to $p_{{\rm s}}=10$ ${\rm bar}$ and we assume that the planet is located at the Venus-Sun distance, that is $a_{{\rm Venus}}=0.723$ ${\rm AU}$. This configuration, characterized in Sect. 4, corresponds to the case illustrated by Fig. 1 of Leconte et al. (2015), and seems thereby a convenient choice for comparisons with this early work.

In Sect. 5, two families of configurations will be studied, both including the reference case. In the first family, the surface pressure is set to $p_{{\rm s}}=10$ ${\rm bar}$ and the semi-major axis varies. Conversely, in the second family, planets have the same orbital radius, $a=a_{{\rm Venus}}$, and various surface pressures.

3.1 Physical setup of the 3D simulations

Apart from the surface pressure and star-planet distance, all simulations are based on a common physical setup. For the stellar incoming flux, the emission spectrum of the Sun is used. The planet is assumed to be dry, with no surface liquid water or water vapour, which allows us to filter out effects associated with the formation of clouds in the study of its atmospheric tidal response. The atmosphere is arbitrarily assumed to be nitrogen-dominated. However, a pure ${\rm N_{2}}$-atmosphere would be an extreme case for radiative transfer owing to the absence of radiator. Hence, we have to set a non-zero volume mixing ratio for carbon dioxyde to avoid numerical issues in the treatment of the radiative transfer with the LMDZ, which was originally designed to study the Earth atmosphere. Although any value could be used, we choose to set the value of the ${\rm CO_{2}}$ volume mixing ratio to that of the Earth atmosphere at the beginning of the twenty-first century, that is $\sim 370$ ${\rm ppmv}$ (e.g. Etheridge et al., 1996). The mass ratio corresponding to this volume mixing ratio being negligible, we use the value of ${\rm N_{2}}$ for the mean molecular mass of the atmosphere, $M_{{\rm atm}}=28.0134$ ${\rm g.mol^{-1}}$ (Meija et al., 2016).

For a perfect diatomic gas, the ratio of heat capacities (also called first adiabatic exponent) is $\Gamma_{1}=1.4$, and it follows that $\kappa=\left(\Gamma_{1}-1\right)/\Gamma_{1}=0.285$ (the parameter $\kappa$ can also be written $\kappa=\mathcal{R}_{\rm GP}/\left(M_{{\rm atm}}C_{p}\right)$, where $\mathcal{R}_{\rm GP}$ and $C_{p}$ stand for the perfect gaz constant and the thermal capacity per unit mass of the atmosphere, respectively). The effects of topography are ignored and the surface of the planet is thus considered as an isotropic sphere of albedo $A_{{\rm s}}=0.2$ and thermal inertia $I_{{\rm gr}}=2000$ ${\rm J.m^{-2}.s^{-1/2}.K^{-1}}$, which is a typical value for Venus-like soils (see e.g. Lebonnois et al., 2010)222This is the value prescribed for nonporous basalts by Zimbelman (1986).. 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 general circulation model (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 of 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 hypotheses are assumed:

• •

Hydrostatic approximation. The pressure and gravitational forces compensate each other exactly along the vertical direction.

• •

Traditional approximation. Components of the Coriolis approximation associated with a vertical motion of fluid particules or generating a force along the vertical direction are ignored.

• •

Thin layer approximation. The thickness of the atmosphere is considered as small with respect to the radius of the planet.

A spatial resolution of $32\times 32\times 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 interpolated every radiative timestep during simulations to determine local radiative transfers. The method is commonly used and has been thoroughly discussed in past studies (e.g. Leconte et al., 2013). We thus refer the readers to these works for a detailed description.

3.2 Extraction of the quadrupolar surface pressure anomaly

For a given planet, of fixed rotation, semi-major axis and surface pressure, the calculation of the quadrupolar surface pressure anomaly follows several steps.

First, the GCM is run for a period $P_{{\rm conv}}$ corresponding to the convergence timescale necessary to reach a steady cycle. Note that this period has to be specified for each doublet $\left(a,p_{{\rm s}}\right)$. As a first approximation, it depends on the radiative timescale of the deepest layers of the atmosphere $\tau_{{\rm rad}}$, which scales as

[TABLE]

where $T_{\rm e}$ stands for the mean effective, or black body, temperature of the atmosphere (see e.g. Showman & Guillot, 2002, Eq. 10). In the reference case, we observe that the atmospheric state has converged toward a steady cycle after $\sim 5800$ Earth Solar days, and we thus use this value for calculations in this section.

After this first step, a simulation is run for 300 Solar days of the planet, defined by $P_{{\rm sol}}=2\pi/\left|\Omega-n_{\star}\right|$, except in the case of the spin-orbit synchronization ($\Omega=n_{\star}$), 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.

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 centered on the substellar point. Since meteorological fluctuations can be considered as a perturbation varying randomly over short timescales, we filter them by folding the surface pressure anomaly over one Solar day.

We finally apply a spherical harmonics transform to the resulting averaged surface pressure snapshot in order to get the complex coefficient $\delta p_{{\rm s};2}^{2,\sigma}$ associated with the semidiurnal tidal mode of degrees $l=2$ and $m=2$ (see Eq. (1)). The method is illustrated by Fig. 2 in the reference case ($a=a_{{\rm Venus}}$ and $p_{{\rm s}}=10$ ${\rm bar}$).

This procedure provides the value of the tidal torque for a given forcing frequency. In practice, the torque is computed over an interval of the normalized frequency $\omega=\left(\Omega-n_{\star}\right)/n_{\star}$ centered on synchronization ($\omega=0$) with $n_{\star}$ fixed, the planet rotation rate being deduced from $\omega$ (the normalized frequency $\omega$ is employed here instead of $\sigma$ to follow along the line by Leconte et al., 2015). Typically, we use $-30\leq\omega\leq 30$ to study the low-frequency regime of the atmospheric tidal response and $-300\leq\omega\leq 300$ to study the high-frequency regime.

The frequency range is thus divided into $N$ intervals, meaning that the whole above procedure has to be repeated $N+1$ times to construct a frequency-spectrum of the tidal torque. The size of an interval is defined as $\Delta\omega\equiv\left(\omega_{\rm sup}-\omega_{\rm inf}\right)/N$. For instance, for the exploration of the parameters space detailed in Sect. 5, $N=20$, $\omega_{\rm inf}=-30$, $\omega_{\rm sup}=30$, and thus $\Delta\omega=3$.

4 Frequency behaviour of the atmospheric tidal torque

The apparent complexity of the physics involved in thermal atmospheric tides requires that we opt for a graduated approach of the problem. Hence, before investigating the dependence of the tidal torque on the planet orbital radius and atmospheric surface pressure as mentioned above, we have to preliminary characterize how it varies with the tidal frequency. To address this question, we consider the reference case ($p_{{\rm s}}=10$ ${\rm bar}$ and $a=a_{{\rm Venus}}$).

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 tide are computed in low-frequency and high-frequency ranges. For convenience, we introduce the function $f_{{\rm GCM}}\left(\sigma\right)$, which is the interpolating function of GCM results with cubic splines. Noting that the tidal torque should be an odd function of the tidal frequency in the absence of rotation (or if the effect of rotation on the tidal response were negligible), we also introduce the function $f_{{\rm odd}}$, defined by

[TABLE]

which is the odd function $f$ minimizing for any $\sigma$ the distance defined by

[TABLE]

The complementary function $f_{{\rm even}}$, such that $f_{{\rm GCM}}=f_{{\rm odd}}+f_{{\rm even}}$, is defined by

[TABLE]

and provides a measure of the impact of Coriolis effects on the tidal torque. The data, the interpolating function $f_{\rm GCM}$ and its components $f_{{\rm odd}}$ and $f_{{\rm even}}$ are plotted in Fig. 3 as functions of the normalized tidal frequency $\omega=\sigma/\left(2n_{\star}\right)$ in linear and logarithmic scales. Additional functions of the frequency are plotted in dashed lines. They correspond to the ab initio analytical (“Ana.”), Maxwell, and parametrized (“Param.”) models that will be introduced and discussed further.

We first consider the low-frequency range ($-30\leq\omega\leq 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 $\omega\sim 5$. We introduce here the maximal value of the peak $q_{{\rm max}}\equiv\max\left\{f_{{\rm odd}}\left(\sigma\right)\right\}$ and the associated frequency $\sigma_{{\rm max}}$, such that $f_{{\rm odd}}\left(\sigma_{{\rm max}}\right)=q_{{\rm max}}$, timescale $\tau_{{\rm max}}\equiv\sigma_{{\rm max}}^{-1}$ and normalized frequency $\omega_{{\rm max}}=\sigma_{{\rm max}}/\left(2n_{\star}\right)$.

The tidal torque is negative for $\sigma<0$ and positive otherwise, which corresponds to the typical behaviour of the thermally induced atmospheric tidal response in the the vicinity of synchronization, as discussed in Sect. 2. As shown by early studies (Gold & Soter, 1969; Ingersoll & Dobrovolskis, 1978; Dobrovolskis & Ingersoll, 1980; Correia & Laskar, 2001), thermal atmospheric tides thus tends to drive the planet away from synchronous rotation and determines its non-synchronized rotation states of equilibrium.

In the zero-frequency limit, the torque scales as $\mathcal{T}\,\propto\,\sigma^{\alpha}$ , with $\alpha\approx 0.73$. In the high-frequency range ($20\lesssim\left|\omega\right|\leq 300$), it scales as $\mathcal{T}\,\propto\,\sigma^{-1}$ with a remarkable regularity (see Fig. 3, bottom left panel) and exhibits a resonance at $\omega\approx 260$. We will see in the next section that these features can be explained using the linear theory of atmospheric tides (Wilkes, 1949; Siebert, 1961; Lindzen & Chapman, 1969).

We note that the spectrum of $f_{{\rm GCM}}$ exhibits a slight systematic asymmetry with respect to the synchronization. This feature is obvious in the low-frequency range, where $\left|f_{{\rm GCM}}\left(-\sigma\right)\right|>\left|f_{{\rm GCM}}\left(\sigma\right)\right|$, and tends to vanish while $\left|\sigma\right|$ increases. Particularly, a small departure between $f_{{\rm GCM}}$ and $f_{{\rm odd}}$ can be observed around the extrema of the tidal torque, and we note that the atmosphere undergoes a non-negligible tidal torque at synchronization ($\sigma=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.

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 atmospheric tides (Wilkes, 1949; Siebert, 1961; Lindzen & Chapman, 1969). In Appendix B, by using an ab initio approach, we compute analytically the atmospheric tidal torque created by the semidiurnal thermal tide in the idealized case of an isothermal atmosphere undergoing the tidal heating of the planet surface. The atmospheric structure is here characterized by the constant pressure height

[TABLE]

where $\mathcal{R}_{\rm s}$ and $T_{{\rm s}}$ designate the specific gas constant and the surface temperature, respectively. It allows to renormalizes the altitude $z$ with the introduction of the pressure height scale $x=z/H$. In the analytic model, we choose for the heat per unit mass inducing the tidal response the vertical profile $J=J_{{\rm s}}{\rm e}^{-b_{{\rm J}}x}$, where $J_{{\rm s}}$ is the heat per unit mass at the planet surface, and $b_{{\rm J}}$ a dimensionless optical depth corresponding to the inverse of the characteristic thickness of the heated layer. Note that the limit $b_{{\rm J}}\rightarrow+\infty$ 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, note that dissipative processes are ignored since they are associated to 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 $\sigma$ compares to the frequency characterizing the turning point, where the vertical wavenumber annihilates (see Appendix B),

[TABLE]

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

Hence, introducing the equivalent depth of the mode,

[TABLE]

we obtain, for $\left|\sigma\right|\leq\sigma_{{\rm TP}}$,

[TABLE]

and, for $\left|\sigma\right|>\sigma_{{\rm TP}}$,

[TABLE]

We recall that $\kappa=\mathcal{R}_{\rm s}/C_{p}$, where $C_{p}$ designates the heat capacity per unit mass, and $\Gamma_{1}=1/\left(1-\kappa\right)$ the adiabatic exponent at constant entropy (Gerkema & Zimmerman, 2008). The solution given by Eqs. (12) and (13) provides a useful diagnosis about the frequency-behaviour of the torque in the high-frequency range.

The most striking feature of this behaviour is the peak that can be observed in Fig. 3 (top and bottom left panels) at the normalized frequency $\omega\approx 260$. This peak correspond to the fundamental resonance of the atmospheric vertical structure associated with the propagation of the Lamb mode (e.g. Lindzen et al., 1968; Bretherton, 1969; Lindzen & Blake, 1972; Platzman, 1988; Unno et al., 1989), which is an acoustic type wave of long horizontal wavelengh. In an inviscid, isothermal atmosphere, the Lamb mode is characterized by the equivalent depth $h_{\rm L}=\Gamma_{1}H$ (Lindzen & Blake, 1972). In the asymptotic regime, where $n_{\star}\ll\left|\Omega\right|$, the characteristic Lamb frequency follows from Eq. (11),

[TABLE]

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

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

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

The other interesting feature highlighted by Fig. 3 is the scaling law of the torque $\mathcal{T}\,\propto\,\sigma^{-1}$ in the range of intermediate frequencies, that is between the thermal and Lamb resonances, typically. This behaviour is described by the analytical model. As discussed before (see Eq. (14)), $\sigma_{{\rm TP}}$ and $\sigma_{{\rm L}}$ are close to each other. The intermediate-frequency range thus corresponds to the case $\left|\sigma\right|<\sigma_{{\rm TP}}$, which leads us to consider the solution given by Eq. (12). We place ourselves in the configuration where characteristic timescales are clearly separated, that is $\left|\sigma\right|\ll\sigma_{{\rm TP}}$ and $n_{\star}\ll\left|\Omega\right|$ in the meantime. As $H/h\,\propto\,\sigma^{-2}$, the preceding condition implies that $H/h\gg 1$. It follows that

[TABLE]

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

[TABLE]

Basically, $\sigma_{{\rm J}}$ is the frequency for which the vertical wavelength of the mode and the characteristic depth of the heated layer are of the same order of magnitude.

From the moment that $\left|\sigma\right|\ll\sigma_{{\rm L}}$ (or $H/h\gg 1$), Eq. (88) can be approximated by the function

[TABLE]

where the associated characteristic timescale $\tau_{{\rm J}}$ and maximal amplitude of the pressure anomaly $q_{{\rm J}}$ are

[TABLE]

We recognize in the form of the function given by Eq. (17) the well-known Maxwell model, which is commonly used to describe the dependence of the tidally dissipated energy on the forcing frequency in the case of solid bodies (e.g. Efroimsky, 2012; Correia et al., 2014). Its use in the case of thermal atmospheric tides is discussed in the next section.

4.3 Discussion on the Maxwell model

Analytic ab initio approaches based on a linear analysis of the atmospheric tidal response – including this work (cf. previous section) – predict that the imaginary part of surface pressure variations can be expressed as a function of the forcing frequency $\sigma=2\left(\Omega-n_{\star}\right)$ as (e.g. Ingersoll & Dobrovolskis, 1978; Auclair-Desrotour et al., 2017a)

[TABLE]

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

Note that other works based upon different approaches converged toward the functional form of the Maxwell model. For instance, Correia & Laskar (2001) used the parametrized function $f\left(\sigma\right)=\sigma^{-1}\left(1-{\rm e}^{-\gamma\sigma^{2}}\right)$ ($\gamma$ 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 analyzing 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 $\sigma\approx 0$). To our knowledge, most of early works using the classical tidal theory to study the spin rotation of Venus and ignoring dissipative processes obtained a torque scaling as $\mathcal{T}\,\propto\,\sigma^{-1}$, and thus singular at the synchronization (e.g. Dobrovolskis & Ingersoll, 1980; Correia & Laskar, 2001, 2003). This is precisely the reason that led Correia & Laskar (2001) to introduce the regular ad hoc parametrized function mentioned above. Conversely, Ingersoll & Dobrovolskis (1978) and, later, Auclair-Desrotour et al. (2017a), derived a Maxwell-like tidal torque analytically by introducing a characteristic thermal time associated with boundary layer processes and radiative cooling. These early results may let think that dissipative processes are a necessary ingredient for a regular tidal torque to exist at the synchronization.

Although dissipative processes definitely regularize the atmospheric tidal torque at the synchronization (e.g. Auclair-Desrotour et al., 2017a), we showed in Sect. 4.2 that regularity also naturally emerges from approaches ignoring them when the vertical structure equation is solved in a self-consistent way. For a sufficiently small frequency, namely $\left|\sigma\right|\ll\sigma_{{\rm J}}$, the torque derived from our analytic solution in the absence of dissipative mechanisms scales as $\mathcal{T}\,\propto\,\sigma$. Therefore, it seems that the singularity at $\sigma=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, note that our analytical model asymptotically converges towards the function obtained by Dobrovolskis & Ingersoll (1980) when the vertical profile of tidal heating tends towards the Dirac distribution used by these authors (i.e. when $b_{{\rm J}}\rightarrow+\infty$).

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 behaviors 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 $\tau_{{\rm 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),

[TABLE]

where $\varepsilon$ stands for the effective fraction of the incoming flux absorbed by the atmosphere. Substituting $\Im\left\{\delta p_{{\rm 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

[TABLE]

the notation $\mathcal{G}$ referring to the gravitational constant.

Assuming that the atmosphere is optically thin in the visible frequency range and that the surface temperature corresponds to a black body equilibrium, we write the mean surface temperature as

[TABLE]

where we have introduced the Stefan-Boltzmann constant $\sigma_{\rm SB}$ and the surface albedo $A_{{\rm s}}$. By substituting $T_{{\rm s}}$ by Eq. (22) in Eq. (21), we obtain that the ratio $q_{{\rm M}}/\tau_{{\rm M}}$ does not depend on the surface pressure and scales as

[TABLE]

with $\varepsilon=1-A_{{\rm s}}$ if the atmosphere is optically thick in the infrared. This relationship between $\tau_{{\rm M}}$ and $q_{{\rm M}}$ means that the two parameters of the Maxwell model (Eq. (19)) can theoretically be reduced to the effective thermal timescale only, which is determined by complex boundary layer and dissipative processes in the general case. The scaling law given by Eq. (23) will be tested using GCM simulations in Sect. 5.

We now compare the Maxwell model to numerical results by assimilating the Maxwell amplitude and timescales to the maximum value of $f_{{\rm odd}}$ and its associated timescale, respectively. The ab initio analytical solution given by Eqs. (12) and (13) (“Ana.”) and its Maxwell-like form, derived for $\left|\sigma\right|\ll\sigma_{{\rm L}}$ and given by Eq. (17) (“Maxwell”), are both plotted in Fig. 3 as functions of the normalized forcing frequency ($\omega$). The numerical values of $\sigma_{{\rm L}}$ and $\sigma_{{\rm TP}}$ used for the plot are determined by the eigenfrequency of the resonance associated with the Lamb mode in GCM simulations, that is $\omega_{{\rm L}}\approx 260$. We arbitrarily choose to set $\sigma_{{\rm J}}=\sigma_{{\rm max}}$ (correspondence between the numerically-derived and the Maxwell maxima), which determines the value of $b_{{\rm J}}$ (i.e. $b_{{\rm J}}\approx 14$). Finally, the maximum $q_{{\rm J}}$ is obtained by fitting the slope in the intermediate frequency-range to numerical results ($q_{{\rm J}}\approx 1042$ ${\rm Pa}$), and provides the value of the parameter $J_{{\rm s}}$ (i.e. $J_{{\rm s}}\approx 0.05$ ${\rm W.kg^{-1}}$).

Figure 3 highlights the fact that the Maxwell model does not allow 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. Note that a smaller departure between the Maxwell and numerical maxima would certainly be obtained by fitting the Maxwell function to the whole spectrum of numerical results, and not only to the peak. However, this would also lead to overestimate the Maxwell timescale, and the fit would not be satisfactory either. A a consequence, a novel parametrized model has to be introduced to better describe the behaviour of the tidal torque in the low-frequency range. This is the purpose of the next section.

4.4 Introduction of a new parametrized model

It has been shown that the ab initio analytic model described in Sect. 4.2 and Appendix B reproduces the main features of the tidal torque in the high-frequency range, namely the resonance associated with the Lamb mode and the asymptotic scaling law $\mathcal{T}\,\propto\,\sigma^{-1}$. However, in the low-frequency range, the behaviour of the torque appears to be a little bit more complex than that predicted by the model, which reduces to a simple Maxwell function. This is not surprising since the atmospheric tidal response at low tidal frequencies involves complex non-linear mechanisms, interactions with mean flows, and dissipative processes, which are clearly outside of the scope of the classical tidal theory used to establish the solution given by Eqs. (12) and (13).

Yet, the frequency dependence of the tidal torque has to be characterized in the vicinity of synchronization as this is where its action of the planetary rotation is the strongest. Our effort has thus to be concentrated on the low-frequency regime and the transition with the high-frequency regime. As they treat the full non-linear 3D dynamics of the atmosphere in a self-consistent way, GCM simulations are particularly useful in this prospect.

To make oneself an intuition of the behaviour of the torque, it is instructive to look at the logarithmic plot of Fig. 3 (bottom left panel), which enables us to identify the different regimes at first glance. We basically observe two tendencies, highlighted in the plot by slopes taking the form of a straight line, in the zero-frequency limit ($\log\left(\omega\right)\lesssim 0.5$) and the high-frequency asymptotic regime ($\log\left(\omega\right)\gtrsim 1.5$). In the interval $0.5\lesssim\log\left(\omega\right)\lesssim 1.5$, the tidal torque reaches a maximum and undergoes an abrupt decay.

Considering these observations, it seems relevant to approximate the logarithm of the torque by linear functions corresponding to the low and high-frequency regimes, and multiplied by sigmoid activation functions. By introducing the notation $\chi\equiv\log\omega$, we thus define the parametrized function as

[TABLE]

where $b_{\rm trans}\approx\log\left(q_{{\rm max}}\right)$ is the level of the transition plateau, $a_{1}$, $b_{1}$, $a_{2}$ and $b_{2}$ the dimensionless coefficients of linear functions describing asymptotic regimes, and $\mathcal{F}_{1}$ and $\mathcal{F}_{2}$ two sigmoid activation functions expressed as

[TABLE]

In these expressions, the dimensionless parameters $\chi_{1}$ and $\chi_{2}$ designate the cutoff frequencies of $\mathcal{F}_{1}$ and $\mathcal{F}_{2}$ in logarithmic scale, and $d_{1}$ and $d_{2}$ the widths of transition intervals. The corresponding tidal torque is given by

[TABLE]

As the scaling law $\mathcal{T}\,\propto\,\sigma^{-1}$ was derived from the ab initio model of Sect. 4.2 in the high-frequency range, we enforce it by setting $a_{2}=-1$. The eight left parameters are then obtained by fitting the function given by Eq. (24) to numerical results (as done previously, the odd function $f_{{\rm odd}}$ is used). We thus end up with

[TABLE]

and plot the model function $\mathcal{F}_{{\rm 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 composition has to be discussed. In Appendix C, we treat the case of a ${\rm CO_{2}}$-dominated atmosphere with a mixture of water and sulfuric acid (${\rm H_{2}SO_{4}}$) comparable to that hosted by the Venus planet. The obtained spectrum and the associated functions introduced above are plotted in Fig. 4, and shall be compared to those computed for the ${\rm N_{2}}$-dominated atmosphere, plotted in Fig. 3 (right panel). Several interesting features may be noted.

First, the tidal torque exerted on the ${\rm CO_{2}}$-dominated atmosphere is more than twice weaker than that exerted on the ${\rm N_{2}}$-dominated atmosphere. Particularly, peaks are strongly attenuated. This results from the vertical distribution of tidal heating. Because of the optical thickness of carbon dioxide in the visible frequency range, an important part of the incoming stellar flux is absorbed above clouds. This is not the case of the ${\rm N_{2}}$-dominated atmosphere, where most of the flux reaches the planet surface and and is re-emited in the infrared frequency range, leading to the thermal forcing of dense atmospheric layers located at high pressure levels.

Second, we observe a greater asymmetry between the negative and positive frequency ranges, the function $f_{{\rm even}}$ being not negligible with respect to $f_{{\rm odd}}$. This is also an effect of the vertical distribution of tidal heating. In the case of the ${\rm N_{2}}$-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 ${\rm N_{2}}$ to ${\rm CO_{2}}$ decreases the contribution of these layers, while it increases the contribution of layers located at pressure levels where the strong zonal jets mentioned above are generated.

Despite the clear interest there is to study the tidal response of ${\rm CO_{2}}$-dominated atmospheres for the similarity of configurations they offer with the Venus planet, we choose to focus in this work on ${\rm N_{2}}$-dominated atmospheres owing to their simpler frequency-behaviour.

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 ${\rm N_{2}}$-dominated atmosphere belongs to the second category. Thus, the role played by the thermal response of the ground should be considered in the present study to explain the observed difference between the obtained tidal torque and the Maxwell model.

Concerning this point, we note that Leconte et al. (2015) included the heat capacity of the surface $C_{{\rm s}}$ in the simplified model they used to establish the Maxwell-like behaviour of the tidal torque (see Sect. 4 in the Material and Methods of their article). Hence, by introducing the heat capacity of the atmosphere/surface system $C=C_{p}p_{{\rm s}}/g+C_{{\rm s}}$ and the emission temperature ($T_{\rm e}$), they expressed the relationship between surface temperature variations $\delta T$ and the variations of the incoming stellar flux $\delta\!\!\>F_{{\rm inc}}$ as

[TABLE]

with $\sigma_{{\rm M}}=4\sigma_{\rm SB}T_{\rm e}^{3}/C$ (the subscript ${\rm M}$ refers to the Maxwell-like form of the function given by Eq. (28)). As we generally observe that $T_{\rm e}\approx T_{{\rm s}}$ in our GCM simulations of a 10 bar atmosphere (the mean surface temperature of the planet is well approximated by the black body equilibrium temperature, given by Eq. (22), in this case), this model implies that $\sigma_{{\rm M}}$ should be always less than $\sigma_{{\rm M}}^{{\rm sup}}=4\sigma_{\rm SB}T_{{\rm s}}^{3}g/\left(C_{p}p_{{\rm s}}\right)$. However, in light of typical values of $\tau_{{\rm M}}$ obtained with the GCM (see Table 2), it appears that the above formula for $\sigma_{{\rm M}}^{{\rm sup}}$ leads to underestimate $\sigma_{{\rm 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,

[TABLE]

where we have introduced the small variations of the incoming stellar flux $\delta F$, surface temperature $\delta T_{{\rm s}}$, radiative heating by the atmosphere $\delta\!\!\>F_{{\rm atm}}$ and diffusive losses in the ground $\delta Q_{{\rm gr}}$ and in the atmosphere $\delta Q_{{\rm atm}}$. Owing to the absence of water, latent heats associated with changes of states are ignored.

In the general case, $\delta T_{{\rm s}}$ and $\delta Q_{{\rm atm}}$ are coupled with the atmospheric tidal response. Particularly, in the Newtonian cooling approximation (i.e. variations of the emitted flux are proportional to temperature variations), $\delta\!\!\>F_{{\rm atm}}$ can be expressed as

[TABLE]

where $K$ designates an effective coefficient of Newtonian cooling. In order to avoid mathematical complications, we ignore this coupling by assuming either that $\left|\delta\!\!\>F_{{\rm atm}}\right|\ll 4\sigma_{\rm SB}T_{{\rm s}}^{3}\left|\delta T_{{\rm s}}\right|$, or, following Bernard (1962), that the variation of the atmospheric flux scales as $\delta\!\!\>F_{{\rm atm}}\,\propto\,\delta T_{{\rm s}}$ in a similar way as the variation of the flux emitted by the ground. This allows us to simplify radiative terms by writing

[TABLE]

where $\epsilon_{{\rm s}}\approx 1$ stands for the effective emissivity of the surface.

With the above approximations, surface temperature variations can be written for a given mode as $\delta T_{{\rm s}}^{\sigma}=\mathfrak{B}_{{\rm s}}^{\sigma}\delta\!\!\>F_{{\rm inc}}^{\sigma}$. We thus end up with (see detailed calculations in Appendix D)

[TABLE]

where $\mathfrak{B}_{{\rm s}}^{0}=\left(4\sigma_{\rm SB}T_{{\rm s}}^{3}\epsilon_{{\rm s}}\right)^{-1}$, and $\tau_{{\rm s}}$ designates the characteristic timescale of the surface thermal response, which depends on the thermal inertia of the ground $I_{{\rm gr}}$ and of the atmosphere $I_{{\rm atm}}$ at the interface, and is expressed as

[TABLE]

We compare this model to numerical results by extracting the $Y_{2}^{2}$ component of the surface temperature distribution $\delta T_{{\rm s};2}^{2,\sigma}$ provided by GCM simulations, as previously done for the surface pressure distribution. The obtained values are plotted in the complex plane in Fig. 5. In this plot, the horizontal and vertical axes correspond to the real and imaginary parts of the normalized transfer function $\mathfrak{B}_{{\rm s}}^{\sigma}/\mathfrak{B}_{{\rm s}}^{0}$ (such that $\delta T_{{\rm s};2}^{2,\sigma}=\mathfrak{B}_{{\rm s}}^{\sigma}\delta\!\!\>F_{{\rm 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\leq\omega<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 $\tau_{{\rm s}}\sim 0.3$ ${\rm 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), $\tau_{{\rm s}}\approx 4.6$ ${\rm days}$ (we use values given by Table 1, set $\epsilon_{{\rm s}}=1$ and neglect $I_{{\rm atm}}$), which shows the limitations of the approach detailed above.

While the forcing frequency increases, the behaviour of the function interpolated using numerical results starts to change radically. In the vicinity of the resonance ($\sigma\sim\sigma_{{\rm max}}$), the imaginary part of $\mathfrak{B}_{{\rm 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 $\delta\!\!\>F_{{\rm atm}}$ (see Eq. (30)) can no longer be neglected, as done in the model. The abrupt variation of the surface thermal lag around the resonance partially explain the behavior of the tidal torque in this range. Nevertheless, to better understand it, one should study the whole dynamics of the atmospheric tidal response, which is beyond the scope of this work.

In the high-frequency range, that is for $\sigma\gg\tau_{{\rm s}}^{-1}$, the model predicts that the amplitude of temperature variations should tend to zero. Yet, we observe that $\left|\delta T_{{\rm 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).

5 Exploration of the parameter space

We now examine the evolution of the tidal torque with the planet semi-major axis ($a$) and surface pressure ($p_{{\rm s}}$).

5.1 Frequency spectra of the tidal torque

Considering the planet defined in Sect. 3.1, we carry out two studies. In study 1, we set $p_{{\rm s}}=10$ ${\rm bar}$ and we compute frequency spectra of the imaginary part of the $Y_{2}^{2}$-surface pressure component in the low-frequency range for $a$ varying from 0.3 to 0.9 ${\rm AU}$. In study 2, we set $a=a_{{\rm Venus}}$, that is $a=0.723$ ${\rm AU}$, and frequency spectra are computed for $p_{{\rm s}}$ varying from 1 to 30 ${\rm bar}$. The reference case characterized in the previous section, and parametrized by $a=a_{{\rm Venus}}$ and $p_{{\rm s}}=10$ ${\rm bar}$, is located at the intersection of the two studies.

Limitations concerning the lower bound of the orbital radius range and the upper bound of the surface pressure range come from the spectra of optical properties used in simulations to compute radiative transfers (see Sect. 3.1), which were produced for temperatures below 710 ${\rm K}$. Indeed, for $a<0.3$ ${\rm AU}$ or $p_{{\rm s}}>30$ ${\rm bar}$, the planet surface temperature exceeds this maximum. As this might lead to erroneous estimations of radiative transfers, we choose not to treat extremal cases, although there is no formal limitation for the GCM to run normally in these conditions.

Radiative transfers also determine the convergence timescale necessary for the atmosphere to reach a steady state, $P_{{\rm conv}}$. For study 1, we use the timescale obtained in the reference case, that is 5800 Earth Solar days, considering that the steady state is reached more rapidly in mosts cases, where the planet is closer to the star (see Eq. 5). Similarly, to take into account the dependence of $P_{{\rm conv}}$ on the planet surface pressure in study 2, we set $P_{{\rm conv}}$ to 1100, 2300, 5800 and 14000 Earth Solar days for $p_{{\rm s}}=1,3,10,30$ ${\rm bar}$, respectively.

The obtained frequency spectra are plotted in Fig. 6 in linear (left) and logarithmic scales (right) for study 1 (top) and study 2 (bottom). In all plots, points designate the results of GCM simulations with the method described in Sect. 3, while solid lines correspond to the associated cubic splines interpolations. The reference case ($a=a_{{\rm Venus}}$ and $p_{{\rm s}}=10$ ${\rm bar}$) is designated by the solid grey line. Numerical values used to produce these plots are given by Table 3 for study 1 and Table 4 for study 2.

We retrieve here the features identified in Sect. 4. The tidal torque exhibits maxima located at the transition between the low-frequency and high-frequency asymptotic regimes. The corresponding peaks are slightly higher in the negative-frequency range than in the positive-frequency range owing to Coriolis effects and the impact of zonal jets on the angular lag of the tidal bulge. As expected, the amplitude of peaks increases with both the incoming stellar flux and the planet surface pressure. Interestingly, the evolution of $q_{{\rm max}}$ and $\sigma_{{\rm max}}$ with $a$ and $p_{{\rm s}}$ looks very regular. This suggests that the dependences of the peak maximum and characteristic timescale on the planet surface pressure and distance to star are well approximated by simple power scaling laws, and it is the case indeed, as shown in Sect. 5.2.

As previously noticed in the study of the reference case, the asymptotic behaviour of the tidal torque in the zero-frequency limit differs from that described by the Maxwell model. Particularly, the logarithmic plot of study 2 (bottom right panel) shows that the torque follows the scaling law $f_{{\rm GCM}}\left(\sigma\right)\,\propto\,\sigma^{1/2}$ in cases characterized by low surface pressures, that is 1 and 3 bar. These cases correspond to the thin-atmosphere asymptotic limit, where thermal tides are driven by diffusion in the ground in the vicinity of the surface. We note that the simplified linear model of the surface thermal response detailed in Sect. 4.6 and Appendix D leads to a surface-generated radiative heating scaling as $\Im\left\{\delta T_{{\rm s}}\right\}\,\propto\,\sigma^{1/2}$ in the zero-frequency limit, which is precisely the dependence observed in Fig. 6.

5.2 Evolution of the thermal peak with the planet semi-major axis and atmospheric surface pressure

Let us now quantify the regular dependence of the peak of tidal torque on the planet orbital radius and surface pressure observed in the preceding section. We have thereby to determine how the two parameters defining the peak – namely its maximum value $q_{{\rm max}}$ and associated timescale $\tau_{{\rm max}}$ – vary with $a$ and $p_{{\rm s}}$.

Thus, for each study, we fit numerical values of $q_{{\rm max}}$ and $\tau_{{\rm max}}$ using a linear regression, formulated as

[TABLE]

where $Y$ designates the logarithm of $q_{{\rm max}}$ (Pa) or $\tau_{{\rm max}}$ (days), $X$ the logarithm of $a$ (AU) or $p_{{\rm s}}$ (bar), and $\alpha$ and $\beta$ the dimensionless parameters of the fit. The values of these parameters are given by Table 2, as well as those of the corresponding coefficients of determination $R^{2}$. We also compute $\log\left(q_{{\rm max}}/\tau_{{\rm max}}\right)$ for comparison with the theoretical scaling law given by Eq. (23) in the case where the dependence of the tidal torque on the forcing frequency is approximated by a Maxwell function.

Linear regressions are plotted in Fig. 7 (blue solid line). In order to provide an estimation of the variability of numerical results, error bars are given for the reference case. These error bars do not literally correspond to a margin of error, but indicate the resolution of the sampling for the frequency and maximum of the tidal torque. For $q_{{\rm max}}$, the amplitude of the error bar is the departure between the maxima of the interpolating function and data. For $\tau_{{\rm max}}$, the two bounds of the error bar are the values associated with the nearest points of the sampling, designated by the subscripts ${\rm inf}$ and ${\rm sup}$, such that $\tau_{{\rm inf}}\leq\tau_{{\rm max}}\leq\tau_{{\rm 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$  ${\rm UA}$, which makes the fit less reliable in this case.

Comparing coefficients of determination in Table 2, we observe that a better fit is systematically obtained for $q_{{\rm max}}$ than for $\tau_{{\rm max}}$. This difference may be explained by the aspect of spectra displayed in Figs. 6 and 3. Since the peak of the tidal torque computed with the GCM is both flatter and larger than that of the Maxwell function, the position of the maximum is more sensitive to small fluctuations than the maximum itself. As a consequence, the variability of $q_{{\rm max}}$ is less than the variability of $\tau_{{\rm max}}$.

Hence, the linear regression fits particularly well the dependence of $q_{{\rm max}}$ on $a$, while the plot of $\tau_{{\rm max}}$ exhibits a relatively important variability with respect to the linear tendency. 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 $\tau_{{\rm max}}$ on $a$, one may also observe that the slope, given by $\alpha=0.86$, is almost twice smaller than that predicted by the scaling law of the radiative timescale given by Eq. (5), that is $\tau_{{\rm max}}\,\propto\,n_{\star}^{-1}\,\propto\,a^{3/2}$.

As regards the ratio $q_{{\rm max}}/\tau_{{\rm 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

[TABLE]

if we assume that $\varepsilon=1-A_{{\rm s}}$ (i.e. the flux reemitted by the ground is entirely absorbed by the atmosphere).

As may be seen, the dependence of $q_{{\rm max}}/\tau_{{\rm max}}$ on the surface pressure is small ($\alpha=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 $\alpha$ (i.e. 1.55 and $3/2$, respectively) is around 3%. However, the value of $\beta$ computed from GCM simulations (2.77) is higher than that predicted by the model (2.49), despite the fact that this latter is an upper estimation. This difference illustrates the limitations of the Maxwell model, which fails to describe the sharp variations of the tidal torque with the tidal frequency when $\left|\sigma\right|\sim\sigma_{{\rm max}}$.

5.3 Scaling laws and generic formula for the tidal torque

By proceeding to a quantitative study of the evolution of the tidal torque maximum with the planet orbital radius and surface pressure, we demonstrated in the preceding section the regularity observed in Fig. 3. The scaling laws given by Table 2 and plotted in Fig. 7 show that frequency-spectra have the same aspect from the moment that the horizontal and vertical axes are rescaled following the obtained dependences on $a$ and $p_{{\rm s}}$. In this section, our purpose is to compute this rescaling in a robust way, by taking into account the whole set of data at our disposal rather than the maximal value of the torque and the associated timescale only. Combining this rescaling with the parametrized model given by Eq. (24), we will obtain a novel generic formula for the frequency-behaviour of the thermally generated atmospheric tidal torque.

The parameter with respect to which axes are rescaled, $a$ or $p_{{\rm s}}$, is denoted by $\mathfrak{p}$, and the considered case is subscripted $j$. A given family is thus composed of $N_{\mathfrak{p}}$ couples of numerical vectors $\left(\sigma_{j},\mathcal{T}_{j}\right)$, with $1\leq j\leq N_{\mathfrak{p}}$ (see Table 3 and Table 4), associated with the value $\mathfrak{p}_{j}$ of the varying parameter $\mathfrak{p}$. For a given couple of vector, one may introduce the associated interpolating function

[TABLE]

We also introduce the renormalized vectors $\hat{\sigma}_{j}$ and $\hat{\mathcal{T}}_{j}$, defined by

[TABLE]

where $\alpha_{1}$ and $\alpha_{2}$ are the exponents characterizing the renormalization.

The size of frequency domains covered by the $\hat{\sigma}_{j}$ vectors varies with $\mathfrak{p}$ in the general case. As a consequence, rescaling axes requires to define the bounds of the largest common interval,

[TABLE]

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

The values of $\alpha_{1}$ and $\alpha_{2}$ are obtained by minimizing the squared difference function

[TABLE]

Note that the parameters derived from these calculations, $q_{0}$ and $\tau_{0}$, slightly differ from $q_{{\rm max}}$ and $\tau_{{\rm max}}$. They stand for the characteristic amplitude and timescale of the peak and not for its maximum value and corresponding forcing period, as $q_{{\rm max}}$ and $\tau_{{\rm max}}$. These parameters are defined by functions of $\mathfrak{p}$ (that is $a$ or $p_{{\rm s}}$), by

[TABLE]

where $\mathfrak{p}_{{\rm ref}}$ designates the value of $\mathfrak{p}$ in the reference case (typically $a_{{\rm Venus}}$ or 10 bar), $q_{{\rm ref}}=2.24$ ${\rm kPa}$ the corresponding maximum of the peak and $\tau_{{\rm ref}}=2.18$ ${\rm days}$ the associated timescale. Hence, we end up for the first family ($p_{{\rm s}}=10$ ${\rm bar}$ and variable $a$) with the scaling laws

[TABLE]

and for the second family ($a=a_{{\rm Venus}}$ and variable $p_{{\rm s}}$), with

[TABLE]

Let us remind here the used units: $a$ is given in AU, $p_{{\rm s}}$ in bar, $q_{0}$ in Pa, and $\tau_{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

[TABLE]

We remind here that $\mathcal{F}_{{\rm par}}$ is the function

[TABLE]

where $\mathcal{F}_{1}$ and $\mathcal{F}_{2}$ are the activation functions defined by

[TABLE]

The parameters characterizing the generic formula given by Eq. (46) take the values

[TABLE]

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 $\tau_{0}\sigma$ (dashed black line). Fig. 8 clearly shows the relevance of the rescaling as regards the first family, where the dependence of the torque on the star-planet distance is investigated. After rescaling, spectra look similar and the model matches them fairly well. As regards the second family, we observe a greater variability of $q_{0}$ and $\tau_{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 atmospheric torque created by the semidiurnal thermal tide as a function of the tidal frequency by extracting the $Y_{2}^{2}$ component of the surface pressure anomaly in simulations.

As a first step, we characterized the variation of the torque with the forcing frequency for a reference case ($p_{{\rm s}}=10$ ${\rm bar}$ and $a=a_{{\rm Venus}}$), and explained its various features with an independent analytical model. As a second step, we explored the parameters space by focusing on the dependence of the tidal torque on the planet orbital radius and atmospheric surface pressure. The obtained results were then used to derive scaling laws characterizing the torque, renormalize the pressure anomaly and forcing period, and finally propose a novel generic parametrized function to model the frequency-behaviour of the torque in a realistic way in the case of a nitrogen-dominated atmosphere.

The first investigation confirmed and extended the results obtained by Leconte et al. (2015). We showed that the torque follows two different asymptotic regimes. In the high-frequency range, the torque decays inversely proportionally to the tidal frequency until it exhibits a resonance peak. These two features are both explained by the analytical solution derived using the ab initio linear theory of atmospheric tides. Particularly, the peak corresponds to a resonance associated with the Lamb mode, an acoustic type wave of wavelength comparable with the planet radius.

In the low-frequency range the torque, which is zero at synchronization, increases following a power law of index ranging from $0.5$ to $0.7$ until it reaches a maximum. While the increase and presence of a maximum are predicted by the analytical solution, the exponent and the value of amplitude of the peak differ significantly. These discrepancies result from the complex interactions between mechanisms laying beyond the scope of standard analytical treatments but resolved in 3D GCM simulations, such as the non-linear effects inherent to the atmospheric dynamics in the vicinity of synchronization, and the strong radiative coupling between the atmosphere and the planet surface. Typically, the low-frequency asymptotic regime of the tidal response is characterized by diurnal oscillations of large amplitude. The resulting differences in the day- and night-side temperature profiles significantly affect the stratification of the atmosphere. This clearly violates the small perturbation approximation upon which the analytic approach is based, and induces a non-linear coupling between the diurnal and semidiurnal oscillations that is important enough to modify the dependence of the tidal torque on the forcing frequency.

The parametrized function that we propose in the present work (given by Eq. (47)) appears as a good compromise as it matches numerical results in a more satisfactory way than the Maxwell model while being defined by a reasonably small number of parameters. It is thus perfectly suited to be implemented in evolutionary models of the rotational dynamics of a planet. Nevertheless, the Maxwell-like analytic solution derived by early studies (e.g. Ingersoll & Dobrovolskis, 1978; Auclair-Desrotour et al., 2017a) provides a first order of magnitude approximation of the torque. It also predicts a relationship between the maximum of the thermal peak and the associated characteristic timescale. By establishing scaling laws governing the evolution of these features with the planet orbital radius and surface pressure, we retrieved numerically this relationship, which is $q_{{\rm max}}/\tau_{{\rm max}}\,\propto\,a^{-3/2}$.

The fact that scaling laws match well numerical results reveals that the torque and the tidal frequency can be normalized by the characteristic amplitude and frequency associated with the low-frequency regime. This was confirmed by the rescaling of spectra, which shows that numerical results obtained in all of the treated cases actually describe the same frequency dependence, whatever the star-planet distance and surface pressure. The combination of the parametrized function and scaling laws derived in this work thus leads to a generic empirical formula for the atmospheric tidal torque in the vicinity of synchronous rotation.

In spite of its limitations in the low-frequency regime, the analytic approach remains complementary with GCM calculations owing to the high computational cost of this method (several days of parallel computation on 80 processors are necessary to produce a spectrum with a sampling of 21 points in frequency). Results obtained from simulations can be used to improve the linear analysis, which provides in return a diagnosis of the physical and dynamical mechanisms involved in the tidal response.

As the study showed evidence of the interest of the numerical method using GCMs in characterizing the atmospheric tidal torque of terrestrial planets, several prospects can be considered for future works. First, the effects of clouds and optical thickness should be investigated owing to their strong impact on the tidal response. The case of an exo-Earth hosting a cloudy atmosphere may be treated in a similar way as the idealized planet of the present study. Second, it would be interesting to better characterize the dependence of the tidal torque on the atmospheric structure using ab initio analytic models. Third, numerical results and the derived generic parametrized function may be coupled to evolutionary models in order to quantify in a realistic way the contribution of the atmosphere to the evolution of the planet rotation over long timescales.

Acknowledgements.

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

Appendix A Normalized spherical harmonics

With the notations introduced in Sect. 2, the normalized spherical harmonics associated with the degrees $l$ and $m$ ($l\in\mathbb{N}$ and $-l\leq m\leq l$) are defined by (e.g. Arfken & Weber 2005)

[TABLE]

where the $P_{l}^{m}$ designate the associated Legendre functions, expressed, for $x\in\left[-1,1\right]$, as

[TABLE]

the $P_{l}$ being the Legendre polynomials, defined by

[TABLE]

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 $\boldsymbol{V}$, pressure $p$, density $\rho$ and temperature $T$ are written as the sum of a spherically symmetrical constant component, subscripted [math], and a time-dependent small perturbation, identified by the symbol $\delta$. 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 $\theta$, and longitude $\varphi$. The constant component of $\boldsymbol{V}$ is ignored, which enforces a solid rotation condition and allows us to simply denote by $\boldsymbol{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)

[TABLE]

the equation of mass conservation,

[TABLE]

the conservation of energy,

[TABLE]

and the perfect gas law,

[TABLE]

the symbol $\partial_{t}$ designating the partial time derivative, $\nabla$ the gradient operator, $\nabla\!\cdot\!$ the divergence and $\dfrac{{\rm D}}{{\rm D}t}\equiv\partial_{t}+\boldsymbol{V}\!\cdot\!\nabla$ the material derivative.

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

Other approximations are convenient to simplify analytic calculations. First, considering that $H\ll R_{\rm p}$, we assume the hydrostatic approximation, as discussed in Sect. 2. Second, we neglect Coriolis components associated with a vertical displacement. This approximation, known as the traditional approximation (e.g. Unno et al. 1989), is appropriate provided that the fluid layer is stably-stratified or thin with respect to the planet radius, which is the case here. Thus the preceding equations reduce to

[TABLE]

where the horizontal part of the velocity divergence $\nabla_{\rm h}\!\cdot\!\boldsymbol{V}$ is defined by

[TABLE]

Beside, we introduce the atmospheric pressure heigh scale,

[TABLE]

and the pressure altitude,

[TABLE]

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 $\sigma$ and a longitudinal degree $m$. Fourier coefficients are functions of the colatitude ($\theta$) and altitude ($z$). Any perturbed quantity $q$ is thus expressed as

[TABLE]

where the $q^{m,\sigma}$ are the Fourier coefficients of the expansion. Under the assumed approximations, Eqs. (57) and (58) are decoupled from the radial momentum equation, which allows us to write $V_{\theta}^{{m,\sigma}}$ and $V_{\varphi}^{{m,\sigma}}$ as functions of $\delta p_{n}^{{m,\sigma}}/\rho_{0}$. By substituting horizontal winds by the obtained expressions in $\nabla_{\rm h}\!\cdot\!\boldsymbol{V}$ (Eq. (63)), and introducing the variable

[TABLE]

the whole system of governing equations given by Eqs. (57 - 62) can be put after some manipulations into the form (see Lindzen & Chapman 1969)

[TABLE]

Here, $\mathcal{F}^{{m,\sigma}}$ is an operator depending on the $x$ coordinate only and $\mathcal{L}^{{m,\nu}}$ the Laplace’s tidal operator, which depends on the $\theta$ coordinate only and is formulated as (e.g. Lee & Saio 1997)

[TABLE]

the quantity $\nu\equiv 2\Omega/\sigma$ designating the so-called spin parameter.

The above separation of coordinates allows us to expand the Fourier coefficients of $G$ as

[TABLE]

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

[TABLE]

and determine the equivalent depth of the mode associated with the triplet $\left(n,m,\sigma\right)$ (e.g. Taylor 1936),

[TABLE]

In the absence of resonances, the semidiurnal tidal response is generally dominated by the fundamental gravity mode, indicated by $n=0\,$333We follow here the indexing notation by Lee & Saio (1997), which associates g-modes with positive $n$ and r-modes to strictly negative $n$., which corresponds to the associated Legendre function $P_{2}^{2}$ in the static case (e.g. Auclair-Desrotour & Leconte 2018). In the high-frequency regime, $\Lambda_{0}^{2,\nu}\approx\Lambda_{0}^{2,1}\approx 11.1$, from the moment that $n_{\star}\ll\left|\Omega\right|$.

Note that this value, denoted by $\Lambda_{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 $\sigma_{\rm R}$, one may show that the eigenvalue of the modes tends to the value of the static case, that is $\Lambda_{0}^{2,\nu}\approx 6$, if $\sigma/\sigma_{\rm R}\rightarrow 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 $\sigma$ to lighten notations. The function $G_{n}^{{m,\sigma}}$ is now simply denoted by $G$, and so on for the tidal heat source, pressure, density, temperature, wind velocity components, eigenvalues and equivalent depths. The usual change of variable $G={\rm e}^{x/2}y$ leads to the vertical structure equation in its canonical form,

[TABLE]

where we have introduced the dimensionless vertical wavenumber $\hat{k}_{x}$, defined by

[TABLE]

The vertical structure equation describe the behaviour of a forced harmonic oscillator, and $\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, $\hat{k}_{x}^{2}\in\mathbb{R}$, and its sign directly determines the nature of waves across the vertical axis. The condition $\hat{k}_{x}^{2}>0$ indicates a propagating mode. Conversely, $\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. Note that we acknowledge the limitations of this academic atmospheric structure regarding real ones, where convective instability leads to a strong temperature gradient near the planet surface. However, this approach appears to be sufficient for the purpose of this appendix.

In the isothermal approximation, the temperature profile is supposed to be invariant with the radial coordinate. In light of Eq. (64), it immediately follows that $H$ is a constant, ${\rm d}H/{\rm d}x=0$, and

[TABLE]

The above expression shows the existence of a turning point for $h=4\kappa H$, where the sign of $\hat{k}_{x}^{2}$ changes. This turning point occurs at the frequency

[TABLE]

In the reference case of the study, GCM simulations provide $T_{{\rm s}}\approx 316$ ${\rm K}$, which, combined with $\mathcal{R}_{\rm s}\approx 297$ ${\rm J.kg^{-1}.K^{-1}}$, gives $H\approx 10.6$ ${\rm km}$ in the isothermal approximation. An estimation of the normalized frequency $\omega_{{\rm TP}}=\sigma_{{\rm TP}}/\left(2n_{\star}\right)$ using Eq. (64) thus gives $\omega_{{\rm TP}}\approx 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 $\hat{k}_{x}^{2}>0$ ($\left|\sigma\right|<\sigma_{{\rm TP}}$) corresponds to an oscillatory regime, while $\hat{k}_{x}^{2}<0$ ($\left|\sigma\right|>\sigma_{{\rm 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

[TABLE]

where $J_{{\rm s}}$ stands for the heat absorbed at the planet surface and $b_{{\rm J}}$ is a dimensionless optical depth characterizing the decay of heating across the vertical coordinate. This profile is derived from the Beer’s law (e.g. Heng 2017) applied to an isothermal atmosphere of constant extinction coefficient, where the variation of flux is supposed to be proportional to the local pressure,

[TABLE]

Integrating this expression from the surface to $x$ and introducing the surface flux $F_{{\rm s}}$, we obtain

[TABLE]

which, linearized in the vicinity of the surface ($x\ll 1$) and assuming $b_{{\rm J}}\gg 1$, can be approximated by

[TABLE]

Since the absorbed heat per unit mass is given by

[TABLE]

and $\rho_{0}\left(x\right)\,\propto\,{\rm e}^{-x}$, we retrieve Eq. (77), with

[TABLE]

With the preceding choices and approximations, the general solution of Eq. (73) can be written as

[TABLE]

the notations $A$ and $B$ designating integrating constants, and $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 $\hat{k}_{x}^{2}>0$. If the convention $\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=+\infty$, 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 $\hat{k}_{x}^{2}<0$, the radiation condition is not appropriate any more since $\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 $\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 $V_{r}=0$ at $x=0$. As the vertical velocity is given by (Lindzen & Chapman 1969)

[TABLE]

the rigid wall condition is expressed as

[TABLE]

at $x=0$. Thus the solution is

[TABLE]

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

The surface pressure variation is readily deduced from $y$ using Eq. (67),

[TABLE]

Taking the imaginary part of the preceding expression, we finally obtain, for $\hat{k}_{x}^{2}>0$ (i.e. $\left|\sigma\right|<\sigma_{{\rm TP}}$),

[TABLE]

and, for $\hat{k}_{x}^{2}<0$ (i.e. $\left|\sigma\right|>\sigma_{{\rm TP}}$),

[TABLE]

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. (88). By assuming that the condition $1\ll H/h\ll\kappa^{-1}b_{{\rm J}}^{2}$ is satisfied, that is for $\left|\sigma\right|\ll\sigma_{{\rm TP}}$ and a small thickness of the heated layer, the torque scales as $\mathcal{T}\,\propto\,\sigma^{-1}$. This corresponds to what may be observed in Fig. 3 in the interval $30<\left|\omega\right|<200$. In the zero-frequency limit, $h\rightarrow 0$ by scaling as $h\,\propto\,\sigma^{2}$ if the dependence of the eigenvalue $\Lambda_{0}$ on the forcing frequency is ignored (see Eq. (72)). It follows that $\mathcal{T}\,\propto\,\sigma$ in this asymptotic regime. The transition between the two regimes occurs for $\left|\sigma\right|\approx\sigma_{{\rm J}}$ with

[TABLE]

Since $H/h\gg 1$ when $\left|\sigma\right|\ll\sigma_{{\rm L}}$ (let us remind that $\sigma_{{\rm L}}$ designates the Lamb frequency given by Eq. (14)), Eq. (88) can be put into the form of the Maxwell function,

[TABLE]

where the associated characteristic timescale $\tau_{{\rm J}}$ and maximal amplitude $q_{{\rm J}}$ are

[TABLE]

Let us now consider the case where $\hat{k}_{x}^{2}<0$, described by Eq. (89). We notice that the surface pressure anomaly diverges for $h=\Gamma_{1}H$. This feature was discussed by early 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 $\omega_{{\rm TP}}$, we obtain $\omega_{{\rm L}}\approx 308$. The resonance can be observed in Fig. 3 for a smaller frequency ($\omega_{{\rm L}}\approx 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 ${\rm N_{2}}$-dominated atmosphere with a small amount of ${\rm CO_{2}}$. Hence, we ignore the effects of clouds and compute the thermal tide of an optically thin atmosphere in the visible frequency range, where the major part of the stellar flux reaches the planet surface without being absorbed.

In this appendix, we consider the case of a planet hosting a Venus-like ${\rm CO_{2}}$-dominated atmosphere with a mixture of water and sulfuric acid (${\rm H_{2}SO_{4}}$) in the same reference configuration ($a=a_{{\rm Venus}}$ and $p_{{\rm s}}=10$ ${\rm 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 ($C_{p}$) to 1000 ${\rm J.kg^{-1}.K^{-1}}$, which is the typical value of $C_{p}$ in the case of Venus (e.g. Seiff et al. 1985), where the parameter decreases from 1181 ${\rm J.kg^{-1}.K^{-1}}$ near the surface to 904 ${\rm J.kg^{-1}.K^{-1}}$ at an altitude of 50 km.

Similarly, the mean molecular mass is set to $43.45$ ${\rm g.mol^{-1}}$, and the volume mixing ratio of water vapor to 20 ${\rm ppm}$ (Moroz et al. 1979). We set the diameter of water particles to 3 ${\rm\mu m}$, which is a typical value in the lower cloud (e.g. Knollenberg & Hunten 1980). To take into account the impact of sulfuric acid on the saturation pressure of water vapor $p_{{\rm H_{2}O}}$, we use the prescription given by Gmitro & Vermeulen (1964) for aqueous sulfuric acid (see Eq. (24) of their article). This prescription is written as a function of the local temperature $T$, as

[TABLE]

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. 9 with the spectrum of the ${\rm N_{2}}$ reference case for comparison.

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

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

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

As mentioned in Sect. 4.6, we follow along the line by Bernard (1962) to study the thermal response of the planet surface. In this approach, we consider the surface-atmosphere interface, located at the altitude $z=0$, and write the power flux budget for a small perturbation in the framework of a frequency linear analysis. Hence, any quantity $q$ can be expressed as $q=q^{\sigma}{\rm e}^{i\sigma t}$, where $\sigma$ is the forcing frequency introduced in Sect. 2. In the following, we omit the superscript $\sigma$ and use $q$ in place of $q^{\sigma}$, 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 $\delta\!\!\>F_{{\rm inc}}$, is absorbed by the planet surface. A fraction $\delta Q_{{\rm gr}}$ of this power is transmitted to the ground by thermal conduction, and an other fraction, $\delta Q_{{\rm atm}}$, is transmitted to the atmosphere through turbulent thermal diffusion. Finally, the increase of surface temperature $\delta T_{{\rm s}}$ generated by $\delta\!\!\>F_{{\rm inc}}$ induces a radiative emission, $\delta\!\!\>F_{{\rm rad}}$, which is expressed as $\delta\!\!\>F_{{\rm rad}}=4\sigma_{\rm SB}T_{{\rm s}}^{3}\delta T_{{\rm s}}$ in the black body approximation (we recall that $\sigma_{\rm SB}$ and $T_{{\rm s}}$ are the Stefan-Boltzmann constant and mean surface temperature introduced in Sect. 3.2, respectively). Since the atmosphere is heated by both the incoming stellar flux and the surface thermal forcing, it undergoes a radiative cooling, similarly as the surface. The flux emitted downward to the surface is denoted $\delta\!\!\>F_{{\rm atm}}$. Thus, the power budget of the thermal perturbation at the interface is expressed as

[TABLE]

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 $\delta\!\!\>F_{{\rm atm}}$. This amounts to assuming either that the emission of the atmosphere toward the planet surface is negligible, or that it is proportional to $\delta T_{{\rm s}}$ (see Bernard 1962). Thus, introducing the surface effective emissivity $\epsilon_{{\rm s}}$, radiative terms can be reduced to

[TABLE]

The next step consists in defining the thermal exchanges resulting from diffusive processes, $\delta Q_{{\rm gr}}$ and $\delta Q_{{\rm atm}}$. These flux are directly proportional to the gradient of the temperature profile anomaly in the vicinity of the interface, and are expressed as

[TABLE]

where $\partial_{z}$ designates the partial derivative in altitude, $\delta T$ the profile of temperature variations, and $k_{{\rm gr}}$ and $k_{{\rm atm}}$ the thermal conductivities of the ground and of the atmosphere at $z=0$, respectively. By introducing the mean density profile of the planet $\rho_{0}$ and the thermal capacity per unit mass of the ground $C_{{\rm gr}}$ (the analogous parameter for the atmosphere being $C_{p}$; see Sect. 3.1), the corresponding diffusivities can be defined by

[TABLE]

Temperature variations in the vicinity of the interface are described by the heat transport equation. We assume that diffusive processes predominates in the $z\rightarrow 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 $\partial_{z}$ in both the solid and atmospheric regions. It follows that

[TABLE]

Solving these two equations with constant $K_{{\rm gr}}$ and $K_{{\rm atm}}$, and ignoring the diverging term in solutions, we end up with

[TABLE]

where we have introduced the frequency-dependent skin thicknesses of heat transport by thermal diffusion in the ground ($h_{{\rm gr}}^{\sigma}$) and atmosphere ($h_{{\rm atm}}^{\sigma}$), expressed as

[TABLE]

The transfer function $\mathfrak{B}_{{\rm s}}^{\sigma}$ defined by Eq. (32) comes straightforwardly when the profiles of $\delta T$ are substituted by Eq. (101) in Eqs. (96) and (94) successively. Introducing $\mathfrak{B}_{{\rm s}}^{0}=\left(4\sigma_{\rm SB}T_{{\rm s}}^{3}\epsilon_{{\rm s}}\right)^{-1}$, we obtain $\delta T_{{\rm s}}=\mathfrak{B}_{{\rm s}}^{\sigma}\delta\!\!\>F$, with

[TABLE]

where $\tau_{{\rm s}}$ can be interpreted as the characteristic timescale of the surface thermal response. The parameter $\tau_{{\rm s}}$ is a function of the thermal inertia of the ground $I_{{\rm gr}}\equiv\rho_{0}\left(0^{-}\right)C_{{\rm gr}}\sqrt{K_{{\rm gr}}}$ and of the atmosphere $I_{{\rm atm}}\equiv\rho_{0}\left(0^{+}\right)C_{p}\sqrt{K_{{\rm atm}}}$,

[TABLE]

The above expression shows that $\tau_{{\rm 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 $\tau_{{\rm s}}\,\propto\,T_{{\rm s}}^{-6}$.

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

As discussed in Sect. 4.6, the transfer function obtained using GCM simulations is well approximated by Eq. (103) in the low-frequency regime (see Fig. 5). But it diverges from the model when the forcing frequency increases, typically for $\left|\sigma\right|\gtrsim\tau_{{\rm 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 $\delta\!\!\>F_{{\rm atm}}$, although it may be very strong. Particularly, this is the case for resonances of the atmospheric tidal response, where $\delta\!\!\>F_{{\rm atm}}$ is increased similarly as the amplitude of pressure and temperature oscillations.

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

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