© The Authors 2024

1 Introduction

For present-day Earth, the semi-diurnal atmospheric tide, driven by the thermal forcing of the Sun and generated via atmospheric pressure waves, describes the movement of atmospheric mass away from the substellar point. Consequently, mass culminates forming bulges on the nightside and the dayside, generating a torque that accelerates the Earth’s rotation. As such, this thermally generated torque counteracts the luni-solar gravitational torque associated with the Earth’s solid and oceanic tides. The latter components typically drive the closed system of the tidal players toward equilibrium states of orbital circularity, copla-narity, and synchronous rotation via dissipative mechanisms (e.g., Mignard 1980; Hut 1981). In contrast, the inclusion of the stellar flux as an external source of energy renders the system an open system where radiative energy is converted, by the atmosphere, into mechanical deformation and gravitational potential energy. Though this competition between the torques is established on Earth, the thermotidal torque remains, at least currently, orders of magnitude smaller.

Interestingly though, this dominance of the gravitational torque over the thermal counterpart admits exceptions. The question of the potential amplification of the atmospheric tidal response was initiated with Kelvin (1882), who invoked the theory of atmospheric tidal resonances, ushering a stream of theoretical studies investigating the normal modes’ spectrum of the Earth’s atmosphere (see Chapman & Lindzen 1970, for a neat and authoritative historical overview). Studies of the Earth’s tidal response spectrum advanced the theory of thermal tides for it to be applied to Venus (Goldreich & Soter 1966; Gold & Soter 1969; Ingersoll & Dobrovolskis 1978; Dobrovolskis & Ingersoll 1980; Correia & Laskar 2001, 2003b; Correia et al. 2003), hot Jupiters (e.g., Arras & Socrates 2010; Auclair-Desrotour & Leconte 2018; Gu et al. 2019; Lee 2020), and near-synchronous and Earth-like rocky exoplanets (Cunha et al. 2015; Leconte et al. 2015; Auclair-Desrotour et al. 2017a, Auclair-Desrotour et al. 2019). Namely, for planetary systems within the so-called habitable zone, the gravitational tidal torque diminishes in the regime near spin-orbit synchronization and becomes comparable in magnitude to the thermotidal torque. Consequently, the latter may actually prevent the planet from precisely reaching its destined synchronous state (Laskar & Correia 2004; Correia & Laskar 2010; Cunha et al. 2015; Leconte et al. 2015).

Going back to Earth, Holmberg (1952) suggested that the thermal tide at present is resonant, and the generated torque is equal in magnitude and opposite in sign to that generated by gravitational tides, thus placing the Earth into a rotational equilibrium with a stabilized spin rate. As this was proven to be untrue for the present Earth (Chapman & Lindzen 1970), Zahnle & Walker (1987) revived Holmberg’s hypothesis by applying the resonance scenario of thermal tides to the distant past. Their suggestion relied on two factors needed to close the gap between the competing torques. The first is the occurrence of a resonance in atmospheric Lamb waves (e.g., Lindzen & Blake 1972) – which we coin as a Lamb resonance – that characterizes the frequency overlap between the fundamental mode of atmospheric free oscillations and the semidiurnal forcing frequency. According to Zahnle & Walker (1987), this resonance occurred when the length of day (LOD) was around 21 h, exciting the thermotidal torque to large amplitudes. Secondly, the gravitational tidal torque must have been largely attenuated in the Precambrian. Recently, Bartlett & Stevenson (2016) revisited the equilibrium scenario and investigated the effect of temperature fluctuations on the stability of the resonance trapping and the Earth’s equilibrium. The authors concluded that the rotational stabilization could have lasted 1 billion years, only to be distorted by a drastic deglaciation event (on the scale that follows the termination of a snowball Earth), thus allowing the LOD to increase again from ~21 h to its present value. Evidently, the occurrence of such a scenario has very significant implications on paleoclimatic studies, with the growing evidence on links between the evolving LOD and the evolution of Precambrian benthic life (e.g., Klatt et al. 2021).

We are fresh out of a study on the tidal evolution of the Earth-Moon system (Farhat et al. 2022), where we focused on modeling tidal dissipation in the Earth’s paleo-oceans and solid interior. There we learned that the tidal response of the oceans, characterized by intermittent resonant excitations, is sufficient to explain the present rate of lunar recession and the estimated lunar age, and it is in good agreement with the geological proxies on the lunar distance and the LOD, leaving little to no place for an interval of a rotational equilibrium (Fig. 1).

On the other hand, major progress has been achieved in establishing the frequency spectrum of the thermotidal response of rocky planets with various approaches ranging from analytical models (Ingersoll & Dobrovolskis 1978; Dobrovolskis & Ingersoll 1980; Auclair-Desrotour et al. 2017a,b), to parameterized models that capture essential spectral features (e.g., Correia & Laskar 2001, 2003a), to fully numerical efforts that relied on the advancing sophistication of general circulation models (GCM; e.g., Leconte et al. 2015; Auclair-Desrotour et al. 2019). The latter work presents, to our knowledge, the first study to have numerically computed the planetary thermotidal torque in the high frequency regime, that is, around the Lamb resonance (Lindzen & Blake 1972). Of interest to us here are two perplexing results that Auclair-Desrotour et al. (2019) established: first, for planets near synchronization, the simplified Maxwellian models often used to characterize the thermotidal torque did not match the GCM simulated response; second, the torque at the Lamb resonance featured only a decelerating effect on the planet. Namely, it acts in the same direction of gravitational tides, and thus the effect required for the rotational stabilization disappeared.

More recently, two studies on the Precambrian LOD stabilization were published. Mitchell & Kirscher (2023) compiled various geological proxies on the Precambrian LOD and established the best piece-wise linear fit to this data compilation. The authors’ analysis depicts that a Precambrian LOD of 19 h was stabilized between 1 and 2 Ga. In parallel, Wu et al. (2023) also attempted to fit a fairly similar set of geological proxies, but using a simplified model of thermal tides. The authors conclude that the LOD was stabilized at ~19.5 h between 0.6 and 2.2 Ga, with a sustained very high mean surface temperature (40–55°C). Although using different approaches, the two studies have thus arrived at similar conclusions. A closer look at the subset of geological data that favored this outcome, however, indicates that both studies heavily rely on three stromatolitic records from the Paleoproterozoic that were originally studied by Pannella (1972a,b). These geological data have been, ever since, identified as unsuitable for precise quantitative interpretation (see e.g., Scrutton 1978; Lambeck 1980; Williams 2000). To this end, we provide a more detailed analysis of the geological proxies of the LOD, and of the model presented by Wu et al. (2023) in a parallel paper dedicated to the matter (Laskar et al. 2023).

With a view to greater physical realism, we aim here to study, analytically, the frequency spectrum of the thermotidal torque, from first principles, interpolating between the low and high frequency regimes. Our motivation is two-fold: first, to provide a novel physical model for the planetary thermotidal torque that better matches the GCM-computed response, and that can be used in planetary dynamical evolution studies; second, to apply this model to the Earth and attempt quantifying the amplitude of the torque at the Lamb resonance and explore the intriguing rotational equilibrium scenario.

2 Ab initio atmospheric dynamics

For an atmosphere enveloping a spherically symmetric planet, we define a reference frame corotating with the planet. In this frame, an atmospheric parcel is traced by its position vector r in spherical coordinates (r, θ, φ), such that θ is the colatitude, φ is the longitude, and the radial distance |r| = R_p + z, where R_p is the planet’s radius and z is the parcel’s atmospheric altitude. The atmosphere is characterized by the scalar fields of pressure p, temperature T, density ρ, and the three-dimensional vectorial velocity field 𝒱. Each of these fields varies in time and space, and is decomposed linearly into two terms: a background, equilibrium state field, subscripted with 0, and a tidally forced perturbation term of significantly smaller amplitude such that p = p₀ + δp, T = T₀ + δT, ρ = ρ₀ + δρ, and 𝒱 = V₀ + V.

Our fiducial atmosphere is subject to the perturbative gravitational tidal potential U and the thermal forcing per unit mass J. We shall define the latter component precisely in Sect. 2.2, but for now it suffices to say that J accounts for the net amount of heat, per unit mass, provided to the atmosphere, allowing for thermal losses driven by radiative dissipation. We take the latter effect into account by following the linear Newtonian cooling hypothesis¹ (Lindzen & McKenzie 1967), where radiative losses, J_rad, are parameterized by the characteristic frequency σ₀; namely J_rad = p₀σ₀/(κρ₀T₀)δT, where κ = (Γ₁ − 1)/Γ₁ = 0.285 and Γ₁ is the adiabatic exponent. Similar to Leconte et al. (2015), we associate with σ₀ a radiative cooling timescale τ_rad = 4τ/σ0.

2.1 The vertical structure of tidal dynamics

We are interested in providing a closed-form solution for the frequency² dependence of the thermotidal torque, which results from tidally driven atmospheric mass redistribution. By virtue of the hydrostatic approximation, this mass redistribution is encoded in the vertical profile of pressure. As such, it is required to solve for the vertical structure of tidal dynamics. With fellow non-theoreticians in mind, we delegate the detailed development of the governing system of equations describing the tidal response of the atmosphere to Appendix A. Therein, we employ the classical system of primitive equations describing momentum and mass conservation (e.g., Siebert 1961; Chapman & Lindzen 1970), atmospheric heat transfer augmented with linear radiative transfer à la Lindzen & McKenzie (1967), and the ideal gas law, all formulated in a dimensionless form.

Aided by the so-called traditional approximation (e.g., Unno et al. 1979, see also Appendix A), the analytical treatment of the said system is feasible as it decomposes into two parts describing, separately, the horizontal and vertical structures of tidal flows. The former part is completely described by the eigenvalue-eigenfunction problem defined as Laplace’s tidal equation (Laplace 1798; Lee & Saio 1997):

$${ℒ}^{m,v}{\text{Θ}}_n^{m,v}=-{\text{Λ}}_n^{m,v}{\text{Θ}}_n^{m,v},$$(1)

where the set of Hough functions $\left\{{\text{Θ}}_n^{m,v}\right\}$ serves as the solution (Hough 1898), $\left\{{\text{Λ}}_n^{m,v}\right\}$ is the associated set of eigenvalues, ℒ^m,v is a horizontal operator defined in Eq. (A.23), while v = 2Ω/σ, where Ω is the rotational velocity of the planet and σ is the tidal forcing frequency. In the tidal system under study, the variables and functions (δp, δρ, δT, V, J, Θ, Λ) are written in the Fourier domain using the longitudinal order m and frequency σ (Eq. (A.20)), and expanded in horizontal Hough modes with index n (Eq. (A.21)). We denote hereafter their coefficients $f_n^{m,v}$ by f_n to lighten the expressions. This horizontal structure of tidal dynamics is merely coupled to the vertical structure via the set of eigenvalues $\left\{{\text{Λ}}_n^{m,v}\right\}$. To construct these sets of eigenfunctions-eigenvalues we use the spectral method laid out by Wang et al. (2016).

The vertical structure on the other hand requires a more elaborate manipulation of the governing system, a procedure that we detail in Appendix B. The outcome is a wave-like equation that describes vertical thermotidal dynamics and reads as:

$$\frac{{\text{d}}^2{\text{Ψ}}_n}{\text{d}{x}^2}+{\widehat{k}}_n^2{\text{Ψ}}_n={\text{Φ}}^{-1}{C}_n.$$(2)

Here, as is the common practice (e.g., Siebert 1961; Chapman & Lindzen 1970), we use the reduced altitude $x={{{\displaystyle \int }}^{\text{}}}_0^z\text{\hspace{0.17em}d}z/H(z)$ as the vertical coordinate, where the pressure scale height H(z) = ℛ_sT₀(z)/g, with ℛ_s being the specific gas constant and g the gravitational acceleration. The quantity Ψ_n(x) is a calculation variable from which, once solved for, all the tidal scalar and vectorial quantities would flesh out (Appendix C). The vertical wave number ${\widehat{k}}_n(x)$ is defined via

$$\begin{array}{l}{\widehat{k}}_n^2(x)=-\frac{1}{4}\{{\left(1-\frac{i\kappa }{\alpha -i}(\gamma -1)\right)}^2+2\frac{\text{d}}{\text{d}x}\left(\frac{i\gamma }{\alpha -i}\right)\hfill \\ -\frac{4\alpha \kappa }{\alpha -i}\left[\beta \gamma {\text{Λ}}_n-\frac{i}{\alpha }(\gamma -1)\right]\}.\hfill \end{array}$$(3)

By virtue of the non-dimensionalization of the governing system of equations, dimensionless control parameters appear in the wavenumber definition. Namely, $\alpha =\sigma /{\sigma }_0,\beta ={\sigma }_{\text{w}}^2/{\sigma }^2$, and $\gamma =N_{\text{B}}^2/N_{\text{B};\text{ iso }}^2$, where we have introduced the characteristic frequency ${\sigma }_{\text{w}}=\sqrt{gH}/R_{\text{p}}$, a typical frequency of Lamb waves. Of significance to the computations that follow is the Brunt–Väisälä frequency, N_B, a measure of the vertical density stratification of the atmosphere against the strength of convection (e.g., Vallis 2017). Under hydrostatic equilibrium, N_B is defined as:

$$N_{\text{B}}^2=\frac{g}{H}\left(\kappa +\frac{\text{d}\mathrm{\hspace{0.17em}ln\hspace{0.17em}}\text{H}}{\text{dx}}\right).$$(4)

We define the Brunt–Väisälä frequency in the limit of an isothermal atmospheric profile, $N_{\text{B};\text{ iso }}=\sqrt{\kappa g/H}$. As such, the parameter γ measures the local atmospheric stability against convection with respect to the stability of an equivalent isothermal atmosphere. These key control parameters, among others, qualitatively determine the regime of the tidal response. We summarize these parameters in Table 1. The right hand side of the wave Eq. (2) combines the function Φ defined as

$$\text{Φ}(x)=\exp \left[\frac{x}{2}-\frac{1}{2}{{{\displaystyle \int }}^{\text{}}}_0^x\frac{i}{\alpha -i}\left(\frac{H{N}_{\text{B}}^2}{g}-\kappa \right)\text{d}{x}^{\prime }\right],$$(5)

and the forcing function C_n(x) which reads:

$$C_n(x)=\frac{\alpha \beta {\text{Λ}}_n{\tilde{J}}_n}{\alpha -i}+\left[\frac{\kappa (\gamma -1)}{\alpha -i}\left({\text{Λ}}_n+\frac{{\partial }_x}{\beta }\right)-i{\partial }_x\left(\frac{{\partial }_x}{\beta }\right)\right]{\tilde{U}}_n.$$(6)

Hereafter, we identify the dimensionless form of the tidal variables by the tilde diacritic. Namely in Eq. (6), the dimensionless thermal forcing function ${\tilde{J}}_n=\kappa /(g\sigma H)J_n$, while the gravitational analog ${\tilde{U}}_n=U_n/v_0^2$, where the reference velocity $v_0=\sigma R_{\text{p}}$ (see Eq. (A.11) for the rest of the dimensionless quantities). As we intend to solve the wave equation in the following sections, what is left for us to quantify the mass redistribution and compute the resulting tidal torque is to retrieve the vertical profile of pressure given the solution of the wave equation, Ψ_n(x). In Appendix C, we derive the vertical profiles of all the tidal variables, and specifically for the dimensionless pressure anomaly we obtain:

$$\tilde{\delta }{p}_n(x)=\frac{1}{i\beta {\text{Λ}}_n}\left(\frac{\text{d}{\tilde{G}}_n}{\text{d}x}-{\tilde{G}}_n\right)+\frac{1}{\beta }\left(1+\frac{1}{\beta {\text{Λ}}_n}\frac{\text{d}}{\text{dx}}\right){\tilde{U}}_n,$$(7)

where $\tilde{\delta }{p}_n(x)=\delta p_n(x)/p_0$, and the calculation variable $\tilde{G}n(x)={\text{Ψ}}_n(x)\text{Φ}(x)$.

2.2 The thermal forcing profile

To solve the nonhomogeneous wave Eq. (2), it is necessary to define a vertical profile for the tidal heating power per unit mass ${\tilde{J}}_n$ (or equivalently in dimensional form, J_n). We adopt a vertical tidal heating profile of the form

$$J_n(x)=J_{\text{s}}{\text{e}}^{-b_{\text{J}}x},$$(8)

where J_s is the heat absorbed at the surface and b_J is a decay rate that characterizes the exponential decay of heating along the vertical coordinate. As we are after a generic planetary model, this functional form of J_n allows the distribution of heat to vary between the Dirac distribution adopted by Dobrovolskis & Ingersoll (1980) where b_J → ∞, and a uniform distribution where the whole air column is uniformly heated (b_J = 0).

To determine J_s, we invoke its dependence on the total vertically propagating flux δF_tot by computing the energy budget over the air column. The net input of energy corresponds to the difference between the amount of flux absorbed by the column and associated with a local increase of thermal energy, and the amount that escapes into space or into the mean flows defining the background profiles. We quantify the fraction of energy transferred to the atmosphere and that is available for tidal dynamics by α_A, where 0 ≤ α_A ≤ 1; the rest of the flux amounting to 1 − α_A escapes the thermotidal interplay. We thus have

$${{{\displaystyle \int }}^{\text{}}}_0^{\infty }J(x){\rho }_0(x)H(x)\text{dx}={\alpha }_{\text{A}}\delta {\text{F}}_{\text{tot }}.$$(9)

To define δF_tot, we establish the flux budget for a small thermal perturbation at the planetary surface. We start with δF_inc, a variation of the effective incident stellar flux, after the reflected component has been removed. δF_inc generates a variation δT_s in the surface temperature T_s. The proportionality between δF_inc and δT_s is parameterized by τ_bl, a characteristic diffusion timescale of the ground and atmospheric surface thermal responses. We detail on this proportionality in Appendix D, but for now it suffices to state that τ_bl is a function of the thermal inertia budgets in the ground, I_gr, and the atmosphere I_atm.

We associate with τ_bl the frequency ${\sigma }_{\text{bl}}={\tau }_{\text{bl}}^{-1}$, a characteristic frequency that reflects the thermal properties of the diffusive boundary layer. It will serve as another free parameter of our tidal model, besides the Newtonian cooling frequency σ₀, and the atmospheric opacity parameter α_A. In analogy to a = σ/σ₀ we define the dimensionless parameter for the boundary layer $\zeta =\sqrt{\left|\sigma \right|{\tau }_{\text{b}1}}=\sqrt{\left|\sigma \right|/{\sigma }_{\text{bl}}}$.

By virtue of the power budget balance established in Appendix D, we define the total propagating flux δF_tot as

$$\delta F_{\text{tot}}=\delta F_{\text{inc}}\left[1-{\mu }_{\text{gr}}\zeta \frac{1+si}{1+(1+si)\zeta }\right].$$(10)

Here, s = sign(σ), and µ_gr is a dimensionless characteristic function weighing the relative contribution of ground thermal inertia to the total inertia budget; namely µ_gr = I_gr/(I_gr + I_atm).

The generic form of the flux in Eq. (10) clearly depicts two asymptotic regimes of thermotidal forcing:

• (i)

  Ignoring the surface layer effects associated with the term on the right, that is, setting ζ = µ_gr = 0, leaves us with ther-motidal heating that is purely attributed to the direct atmospheric absorption of the incident flux. This limit can be used to describe the present understanding of thermotidal forcing on Earth where, to first order, direct insolation absorption in the shortwave by ozone and water vapor appears sufficient to explain the observed tidal amplitudes in barometric measurements (e.g., Chapman & Lindzen 1970; Schindelegger & Ray 2014). Nevertheless, it is noteworthy that the observed tidal phases of pressure maxima could not be explained by this direct absorption, a discrepancy later attributed to an additional semidiurnal forcing, namely that of latent heat release associated with cloud and raindrop formation (e.g., Lindzen 1978; Hagan & Forbes 2002; Sakazaki et al. 2017).

• (ii)

  Allowing for the surface layer term on the other hand (ζ ≠ 0, µ_gr ≠ 0) places us in the limit where the ground radiation in the infrared and heat exchange processes occurring in the vicinity of the surface would dominate the thermotidal heating. The total tidal forcing in this case is nonsynchronous with the incident flux due to the delayed thermal response of the ground, which here is a function of τ_bl. This limit better describes dry Venus-like planets, as is the fiducial setting studied using GCMs in Leconte et al. (2015) and Auclair-Desrotour et al. (2019).

Finally, as we are interested in the semi-diurnal tidal response, we decompose the thermal forcing in Appendix E to obtain the amplitude of the quadrupolar component as $\delta F_{\text{inc }}=\delta F_{22}=(\sqrt{30\pi }/16)F_{\star }$, where $F_{\star }=L_{\star }/4\pi a_{\text{p}}^2,L_{\star }$, L_⋆ being the stellar luminosity, and a_p the star-planet distance.

3 The tidal response

3.1 The tidal torque in the neutral stratification limit

Under the defined forcing in the previous section, to solve the wave equation analytically, a choice has to be made on the Brunt–Väisälä frequency, N_B (Eq. (4)), which describes the strength of atmospheric buoyancy forces and consequently the resulting vertical temperature profile. Earlier analytical solutions have been obtained in the limit of an isothermal atmosphere (Lindzen & McKenzie 1967; Auclair-Desrotour et al. 2019), in which case the scale height H becomes independent of the vertical coordinate, and by virtue of Eq. (4), $N_{\text{B}}^2=\kappa g/H=\text{const}$.

Motivated by the Earth’s atmosphere, where the massive troposphere (~80% of atmospheric mass) controls the tidal mass redistribution, we derive next an analytical solution in a different, and perhaps more realistic limit. Namely, the limit corresponding to the case of a neutrally stratified atmosphere, where $N_{\text{B}}^2=0$. In fact, $N_{\text{B}}^2$ can be expressed in terms of the potential temperature Θ₀ (e.g., Sect. 2.10 of Vallis 2017):

$$N_{\text{B}}^2=\frac{g}{H}\frac{\text{d}\mathrm{\hspace{0.17em}ln\hspace{0.17em}}{\text{Θ}}_0}{\text{dx}},$$(11)

whereby the stability of the atmosphere is controlled by the slope of Θ₀. To this end, atmospheric temperature measurements (e.g., Figs. 2.1–2.3 of Pierrehumbert 2010) clearly depict that the troposphere is characterized by a negative temperature gradient, and a weak potential temperature gradient, which is closer to an idealized adiabatic profile than it is to an idealized isothermal profile. Moreover, the heating in the troposphere generates strong convection and efficient turbulent stirring, thus enhancing energy transfer and driving the layer toward an adiabatic temperature profile, away from the isothermal profile. Consequently, the adiabatic temperature profile prohibits the propagation of buoyancy-restored gravity waves, which compose the baroclinic component of the atmospheric tidal response (e.g., Gerkema & Zimmerman 2008). This leaves the atmosphere with the barotropic component of the tidal flow, a feature that is consistent with tidal dynamics under the shallow water approximation (Appendix A).

Hereafter, we focus on the thermotidal heating as the only tidal perturber, and we ignore the much weaker gravitational potential Ũ. It follows, in the neutral stratification limit, that γ = 0 (Table 1), and the vertical wavenumber (Eq. (3)) reduces to³

$${\widehat{k}}^2={\left[\frac{1+\kappa +i\alpha }{2(\alpha -i)}\right]}^2.$$(12)

It also follows that the background profiles of the scalar variables read as (Auclair-Desrotour et al. 2017a):

$$\begin{array}{l}{p}_0(x)=p_0(0)e^{-x},{\rho }_0(x)=\frac{p_0(0)}{gH(0)}{e}^{(\kappa -1)x},\hfill \\ T_0(x)=\frac{gH(0)}{{ℛ}_{\text{s}}}{\text{e}}^{-\kappa x}.\hfill \end{array}$$(13)

We thus obtain for the heating profile (using Eqs. (8), (9), and (10))

$$J_{\text{s}}=\delta F_{22}\frac{{\alpha }_{\text{A}}g\left(b_{\text{J}}+1\right)}{p_0(0)}\left[1-{\mu }_{\text{gr}}\zeta \frac{1+si}{1+(1+si)\zeta }\right].$$(14)

As such, the wave Eq. (2) is rewritten as

$$\frac{{\text{d}}^2{\text{Ψ}}_{\text{n}}}{{\text{dx}}^2}+{\widehat{k}}^2{\text{Ψ}}_n={\mathcal{𝒜}}_n{\text{e}}^{-ℬx},$$(15)

where the complex functions 𝒜_n and ℬ are defined as:

$${\mathcal{𝒜}}_n=\frac{\kappa {\text{Λ}}_n}{R_{\text{p}}^2{\sigma }^3}\frac{{\alpha }^2+i\alpha }{{\alpha }^2+1}{J}_{\text{s}}\text{, }$$(16) $$ℬ=\frac{1}{2\left({\alpha }^2+1\right)}\left[\left(2b_{\text{J}}+1\right)\left({\alpha }^2+1\right)-\kappa +i\kappa \alpha \right].$$(17)

The wave Eq. (15) admits the general solution

$${\text{Ψ}}_n\left(x\right)=c_1{\text{e}}^{i\widehat{k}x}+c_2{\text{e}}^{-i\widehat{k}x}+\frac{{\mathcal{𝒜}}_n}{{ℬ}^2+{\widehat{k}}^2}{\text{e}}^{-ℬx}.$$(18)

We consider the following two boundary conditions:

– First, the energy of tidal flows, 𝒲, should be bounded as x → ∞. In Appendix F, we derive the expression of the tidal energy following Wilkes (1949), and it scales as $\mathcal{W}\propto \left|\text{Ψ}{|}^2\right|\text{Φ}{|}^2$. Accordingly, the non-divergence of the flow condition is satisfied if one sets c₂ = 0 and takes the proper sign of the wavenumber (Eq. (12)), namely:

$$\widehat{k}=\frac{1}{2\left({\alpha }^2+1\right)}\left[\kappa \alpha +i\left(1+\kappa +{\alpha }^2\right)\right].$$(19)

– The second condition is the natural wall condition imposed by the ground interface, which enforces ${\tilde{V}}_{r;n}(x=0)=0$. We derive the expression of the profile of the vertical velocity in Appendix C, and by virtue of Eq. (C.6), this condition allows us to write c₁ in the form:

$$c_1\left(x\right)=\frac{{\mathcal{𝒜}}_n}{ℬ+{\widehat{k}}^2}\times \frac{ℬ-\frac{1}{2}\left(1+\frac{i\kappa }{\alpha -i}\right)-\beta {\text{Λ}}_n+1}{i\widehat{k}+\frac{1}{2}\left(1+\frac{i\kappa }{\alpha -i}\right)+\beta {\text{Λ}}_n-1}.$$(20)

Under these boundary conditions, we are now fully geared to analytically compute the solution of the wave equation, Ψ_n(x) (or equivalently ${\tilde{G}}_n(x)$), but we are specifically interested in retrieving a closed-form solution of the quadrupolar tidal torque. The latter takes the general form⁴ (Appendix G):

$$\mathcal{𝒯}=\sqrt{\frac{6\pi }{5}}\frac{M_{\star }}{M_{\text{p}}}\frac{R_{\text{p}}^6}{a_{\text{p}}^3}\text{Im}\left\{\delta p_{\text{s}}\right\}.$$(21)

Here M_⋆ and M_p designate the stellar and planetary masses respectively, and Im refers to the imaginary part of a complex number, the latter in this case being the quadrupolar pressure anomaly at the surface $\delta p_{\text{s}}=\delta p_2^{2,v}(x=0)$. We further note that while this torque is computed for the atmosphere, it does act on the whole planet since the atmosphere is a thin layer that features no differential rotation with respect to the rest of the planet.

Taking the solution Ψ(0) of Eq. (18) (with c₂ = 0 and c₁ defined in Eq. (20)), we retrieve δp_s from Eq. (7). After straightforward, but rather tedious manipulations, we extract the imaginary part of the pressure anomaly and write it in the simplified form:

$$\begin{array}{l}\text{Im}\left\{\delta p_{\text{s}}\right\}={\alpha }_{\text{A}}\delta F_{22}\frac{\kappa g{\text{Λ}}_2}{R_{\text{p}}^2{\sigma }^3}\frac{\left(𝒳\alpha +𝒴\right)\alpha }{\left(1+2\zeta +2{\zeta }^2\right)}\hfill \\ \times \underset{\text{position of the Lamb resonance}}{\underbrace{{\left[{\left(\kappa -\beta {\text{Λ}}_2+1\right)}^2+{\alpha }^2{\left(\beta {\text{Λ}}_2-1\right)}^2\right]}^{-1}}},\hfill \end{array}$$(22)

where we have defined the complex functions 𝒳 and 𝒴 as

$$\begin{array}{l}𝒳=\left(\beta {\text{Λ}}_2-1\right)\left[2{\zeta }^2\left(1-{\mu }_{\text{gr}}\right)+\zeta \left(2-{\mu }_{\text{gr}}\right)+1\right],\hfill \\ 𝒴=-s{\mu }_{\text{gr}}\zeta \left(\kappa -\beta {\text{Λ}}_2+1\right).\hfill \end{array}$$(23)

We note that we provide the full complex transfer function of the surface pressure anomaly, along with further analysis on its functional form in Appendix H. In Appendix I, we show that the same form of the tidal solution given by Eq. (22) holds for moist and dry adiabatic profiles, subject to interchanging few effective parameters by their humidity-dependent analogs.

Before embarking on any results, we pause here for a few remarks on the provided closed-form solution of the torque.

• The parameter α_A, defined earlier (Eq. (9)) as the fraction of radiation actually absorbed by the atmosphere, can evidently be correlated with the typical transmission function of the atmosphere and therefore its optical depth. Presuming that ther-motidal heating on Earth is driven by ozone and water vapor, α_A can then characterize atmospheric opacity parameter in the visible. Explicitly showing this dependence now takes us too far afield, though we compute and infer estimates of α_A in Sect. 4.1 and Appendix J.

• The quadrupolar component of the equilibrium stellar flux, entering through a fraction of F_⋆ (Appendix E), is directly proportional to the stellar luminosity L_⋆. Standard models suggest that the Sun’s luminosity was around 80% of its present value ~3 Ga (Gough 1981). Such luminosity evolution of Sun-like stars can be directly accommodated in the model if one were to study the evolution of the tidal torque with time.

• As we mentioned earlier, upon separating the horizontal and vertical structure of tidal dynamics, the only remaining coupling factor between the two structures is the eigenvalue of horizontal flows, A_n, which reduces in our case to the dominant fundamental mode Λ₂. Noting that we have dropped the superscripts, we remind the reader that for the semidiurnal (m = 2) response, ${\text{Λ}}_2={\text{Λ}}_2^{2,v}$, thus A is frequency-dependent in the general case. The Earth however, over its lifetime, lives in the asymptotic regime of v ≈ 1 since 2Ω ≫ n_⋆, thus it is safe to assume that Λ₂ is invariant over the geological history with a value of 11.159 that we compute using the spectral method of Wang et al. (2016).

• Of significance to us in the Precambrian rotational equilibrium hypothesis is the tidal frequency, and consequently the LOD, at which the Lamb resonance occurs. It is evident from the closed-form solution (Eq. (22)) that the position of the resonance is controlled by the highlighted term. Had it not been for the introduced radiative losses, entering here through a, this term would have encountered a singularity at the spectral position of the resonance, namely for βΛ₂ = 1. Here, however, the amplitude of the tidal peak is finite, and its position is a function of the planetary radius, gravitational acceleration, average surface temperature, eigenvalue of the fundamental Hough mode of horizontal flows, and the Newtonian cooling frequency. We detail further on this dependence in Sect. 4.2.

In Fig. 2, we plot the spectrum of the tidal response for a fiducial system in terms of the normalized surface pressure anomaly over a wide range of tidal frequencies covering the low and high frequency regimes. The system describes a Venus-like dry planet (M_p = 0.815M_⊕, R_p = 0.95 R_⊕, a_p = 0.73 au, g = 8.87 m s⁻²), with a 10 bar atmosphere and a scale height at the surface H₀ = 10 km, thermally forced by a solar-like star (M_⋆ = 1 M_⊙, L_⋆ = 1 L_⊙). We further ignore here the thermal inertia in the ground and the atmosphere by taking σ_bl → ∞, or ζ → 0, thus assuming an synchronous response of the ground with the thermal excitation.

Key tidal response features are recovered in this spectrum: first, we obtain a tidal peak near synchronization (ω = 0) that generates a positive torque for σ > 0 and a negative torque for σ < 0, driving the planet in both cases away from its destined spin-orbit synchronization due to the effect of solid tides (e.g., Gold & Soter 1969; Correia & Laskar 2001; Leconte et al. 2015). The peak has often been modeled by a Maxwellian functional form, though this form does not always capture GCM-generated spectra when varying the planetary setup (e.g., Auclair-Desrotour et al. 2019). Second, we recover the Lamb resonance in the high frequency regime. The resonance is characterized here by two symmetric peaks of opposite signs. Thus upon passage through the resonance, the thermotidal torque shifts from being a rotational brake to being a rotational pump. In this work, we are more interested in the high frequency regime, thus we delegate further discussion and analysis on the low frequency tidal response to a forthcoming work, and we focus next on the Lamb resonance.

3.2 The longwave heating limit: breaking the symmetry of the Lamb resonance

We now allow for variations of the characteristic time scale associated with the boundary layer diffusive processes, τ_bl (Eq. (D.12)), or equivalently σ_bl. Variations in σ_bl are physically driven by variations in the thermal conductive capacities of the ground and the atmosphere, and are significant when infrared ground emission and boundary layer turbulent processes contribute significantly to the thermotidal heating.

In such a case, the value of σ_bl plays a significant role in the tidal response of the planet. Namely, the ratio σ/ σ_bl determines the angular delay of the ground temperature variations. For our study of the global tidal response, this frequency ratio determines whether the ground response is synchronous with the thermal excitation (when σ ≪ σ_bl), meaning thermal inertias vanish, the ground and the surface layer do not store energy, and the ground response is instantaneous; or if due to the combination of thermal inertias, the energy reservoir of the ground is huge, and the ground response lags the excitation, imposing another angular shift on the generated tidal bulge (when σ ≳ σ_bl). We now reap from the analytical model the signature of σ_bl in order to explain the Lamb resonance asymmetry – as opposed to its symmetry in Fig. 2 – observed in GCM simulations of an atmosphere forced by a longwave flux (Auclair-Desrotour et al. 2019).

In Fig. 3, we plot the tidal spectrum around the Lamb resonance, in terms of the normalized pressure anomaly at the surface, for different values of σ_bl. For σ_bl = 5 × 10⁻² s⁻¹, the almost instantaneous response of the ground leaves us with two pressure peaks that are symmetric around the resonant frequency. Decreasing σ_bl and allowing for a delayed ground response, the two pressure peaks of the resonance are attenuated in amplitude, but not with the same magnitude; namely, the amplitude damping is stronger against the positive pressure peak. Decreasing σ_bl to 10⁻⁵ s⁻¹ in the panel in the middle, the positive pressure peak completely diminishes, leaving only the negative counterpart. Decreasing σ_bl further, both peaks are amplified, thus the positive peak emerges again. However, the spectral position of the peaks is now opposite to what it was in the limit of an instantaneous ground response.

Given the direct proportionality between the tidal torque and the surface pressure anomaly (Eq. (21)), the effect of thermal inertia thus contributes to the rotational dynamics when encountering the Lamb resonance. If a planet is decelerating and is losing rotational angular momentum, L_Ω, due to solid or oceanic gravitational tides, ω decreases, and the planet encounters the resonance from the right. In the first panel of Fig. 3, the thermotidal torque in this regime is also negative, thus it complements the effect of gravitational tides. When the resonance is encountered, the thermotidal torque shifts its sign to counteract the effect of gravitational tides, with an amplified effect in the vicinity of the resonance. However, with the introduction of thermal inertia into the linear theory of tides, the L_Ω-pumping part of the atmospheric torque is attenuated, and for some values of σ_bl, it completely disappears. This modification of the analytical theory allows us to explain the asymmetry of the Lamb resonance depicted in the 3D GCM simulations of Auclair-Desrotour et al. (2019). In Appendix K, we show that we are able to recover from our model the essential features of the tidal spectrum computed in the mentioned simulations.

To understand the signature of the surface response further, in Fig. 4, we generate snapshots of the tidal pressure variation in the equatorial plane, seen from a top view. The snapshots thus show the thermally induced atmospheric mass redistribution and the resulting tidal bulge, if any. To generate these plots, we compute the vertical profile of the pressure anomaly from Eq. (7), and augment it with the latitudinal and longitudinal dependencies from Eqs. (A.20)–(A.21). As the massive troposphere dominates the tidal mass redistribution, we use the mass-based vertical coordinate ϛ = p/p_s (i.e., x = – ln ϛ, and ϛ ranges between 1 at the surface and 0 in the uppermost layer).

In Fig. 4, we show the tidal bulge as the planet passes through the Lamb resonance, for two values of σ_bl that correspond to the limits of synchronous atmospheric absorption (top row), and a delayed thermal response in the ground (bottom row). First, the accumulation of mass and its culmination on a tidal bulge is indicated by the color red, with varying intensity depicting varying pressure amplitudes. In the case of synchronous atmospheric absorption, for ω = 253, that is, before encountering the resonance, the bulge leads the substellar point and acts to accelerate the planet’s rotation. Increasing ω and encountering the resonance, the bulge reorients smoothly toward lagging the substellar point thus decelerating the planet’s rotation. This behavior is consistent with the established response spectrum in the first panel of Fig. 3, and is relevant to the Earth’s case, assuming that thermotidal heating is predominantly driven by direct synchronous absorption. In the bottom row, the delayed response of the ground imposes another shift on the bulge: for the prescribed value of σ_bl, the passage through the resonance only amplifies the response, but the bulge barely leads the tidal vector, leaving us with a tidal torque that mainly complements the gravitational counterpart, as seen in the fourth panel of Fig. 3.

From what preceded, the reader can find it quite natural that the effect of thermal inertias in the ground and the boundary layer should be accounted for when studying planetary rotational dynamics using the linear theory, especially under longwave forcing. The results also make it tempting to revisit these effects in the case of the dominant shortwave forcing on Earth, as they have been often ignored from the theory (e.g., Chapman & Lindzen 1970) on the basis of the small-amplitude non-migrating tidal components they produce (e.g., Schindelegger & Ray 2014).

4 The stabilized Precambrian LOD hypothesis

So where does all this leave us with the Precambrian rotational equilibrium hypothesis? The occurrence of this scenario straddles several factors, the most significant of which is that the Lamb resonance amplifies the thermotidal response when the opposing gravitational tide is attenuated. Consequently, to investigate the scenario, the two essential quantities that need to be well constrained are the amplitude of the thermotidal torque when the resonance is encountered, and the geological epoch of its occurrence. Having provided a closed-form solution for the tidal torque, it is straightforward for us to investigate these elements.

4.1 The amplitude of the resonant thermal tide: a parametric study

Constraints on the amplitude of the gravitational tide during the Precambrian are model-dependent. The study in Zahnle & Walker (1987), and later in Bartlett & Stevenson (2016), relied on rotational deceleration estimates fitted to match the distribution of geological proxies available at the time (e.g., Lambeck 1980). Specifically, the estimate of the Precambrian gravitational torque relied on the tidal rhythmite record preserved in the Weeli–Wolli Banded Iron formation (Walker & Zahnle 1986). The record is fraught with multiple interpretations featuring different inferred values for the LOD (Williams 1990, 2000), altogether different from a recent cyclostratigraphic inference that roughly has the same age (Lantink et al. 2022, see Fig. 1 for the geological data points ~2.45 Ga). Nevertheless, the claim of an attenuated Precambrian torque still holds, as the larger interval of the Pre-cambrian is associated with a “dormant” gravitational torque phase, lacking any significant amplification in the oceanic tidal response, in contrast with the present state where the oceanic response lives in the vicinity of a spectral resonance (e.g., Farhat et al. 2022).

That said, we explore the atmospheric parameter space of our analytical model to check the potential outcomes of the torques’ competition. Given that the dominant thermotidal forcing on Earth is the direct absorption of the incident flux, we consider the synchronous limit of ζ → 0, whereby the Lamb resonance is symmetric (first panel of Fig. 3; top row of Fig. 4). In Fig. 5, on a grid of values of our free parameters σ₀ and α_A, we contour the surface of the maximum value of the imaginary part of the positive pressure anomaly that is attained when the Lamb resonance is encountered. The two parameters have a similar signature on the tidal response. Moving vertically upward and increasing the rate of Newtonian cooling typically attenuates the amplitude of the peak. For very high cooling rates corresponding to σ₀ ≳ 10⁻⁴s⁻¹, we severely suppress the amplified pressure response around the resonance. Conversely, for values of σ₀ ≲ 10^−6.5 s⁻¹, we approach the adiabatic limit of the tidal model where the Lamb resonance becomes a singularity. A similar signature is associated with increasing the opacity parameter of the atmosphere.

On the contour surface, we highlight with the solid blackisoline the pressure anomaly value required to generate a thermotidal torque of equal magnitude to the Precambrian gravitational tidal torque. The latter (~1.13 × 10¹⁶ Nm) is roughly a quarter of the present gravitational torque (~4.51 × 10¹⁶ Nm) (e.g., Zahnle & Walker 1987; Farhat et al. 2022), thus requiring, via Eq. (21), Im{δp_s} on the order of 880 Pa⁵. This isoline bounds from below a cornered region of the parameter space where the thermal tide is sufficiently amplified upon the resonance encounter. It is noteworthy that this Precambrian value of the torque is the minimum throughout the Earth’s history. We mark by the dashed isoline, for comparison, the threshold needed if the Lamb resonance is encountered in the late Paleozoic/early Mesozoic. The solid gray region on the left side of the parameter space is bounded by the isoline corresponding to the present value Im{δp_s} = 224 Pa (Schindelegger & Ray 2014). Thus it defines to the left an area where the present thermal tide is stronger than it would be around the resonance.

We take this parametric study one step further to study whether typical values of the parameters σ₀ and α_A can place the Earth’s atmosphere in the identified regions. Stringent constraints on σ₀ are hard to obtain for the Earth since σ₀ is an effective parameter that in reality is dependent on altitude. Furthermore, in the linear theory of tides, we are forced to ignore the layer-to-layer radiative transfer and assume a gray body atmospheric radiation directly into space. However, radiative transfer can be consistently accommodated in numerical GCMs using the method of correlated k-distributions (e.g., Lacis & Oinas 1991) as performed in Leconte et al. (2015) and in Auclair-Desrotour et al. (2019), both studies using the LMD GCM (Hourdin et al. 2006). In fact, Leconte et al. (2015) fitted their numerically obtained atmospheric torques to effective values of σ₀ for various atmospheric parameters (see Table 1 of Leconte et al. 2015). The closest of these settings to the Earth yields a radiative cooling timescale τ_rad = 32 days. We consider a range of τ_rad for our model between half and twice this value, and we highlight this range in Fig. 5 with the horizontal shaded area⁶. It is noteworthy here that these dissipative timescales are only relevant for radiative cooling. Other dissipative mechanisms, such as friction (e.g., Lindzen & Blake 1972), could be as efficient as Newtonian cooling and could act over shorter timescales. This said, the parametric analysis that follows, restricted to radiative cooling as the only dissipative mechanism, can be considered to lie on the optimistic limit of resonant tidal amplification.

Another constraint on the free parameters emerges from present in situ barometric observations of the semidiurnal (S₂) tidal response. We use the analysis of compilations of measurements performed in Haurwitz & Cowley (1973); Dai & Wang (1999); Covey et al. (2014) and Schindelegger & Ray (2014), which constrain the amplitude of the semi-diurnal surface pressure oscillation to within 107–150 Pa, occurring around 0945 LT. The narrow shaded area defines the region of parameter space with which we can explain these observables using the present semi-diurnal frequency, placing the opacity parameter in the region α_A ~ 14%. In Appendix J, we compute estimates of the present value of α_A by studying distributions of heating rates that are obtained either by direct measurements of the Earth’s atmosphere (Chapman & Lindzen 1970), or using GCM simulations (Vichare & Rajaram 2013). Our analysis of the data suggests that the efficiency parameter is around α_A ~ 17–18%, which is consistent with the S₂ constraint we obtain. Finally, it is also noteworthy how the plotted S₂ constraint is insensitive to variations in σ₀ over a wide interval, which prohibits the determination of the present value of σ₀ using this constraint.

The overlap of the parametric constraints allows for a narrow region in parameter space where the thermotidal response is sufficient to satisfy the Precambrian rotational equilibrium condition. The present thermotidal torque needs to be amplified by a factor of 3.9 to reach the absolute minimum of the opposing gravitational torque⁷ in the Precambrian, and by a factor of 12.3 to reach the Paleozoic/Mesozoic threshold. Our parametric exploration partially allows for the former but suggests that the latter level of amplification is unlikely. It is important to also note that larger amplification factors would be required if one were focused on the modulus of the pressure oscillation, or its real value, rather than its imaginary part. This derives from Fig. H.1 where we show that the amplification in the imaginary part is almost half that of the modulus of the surface pressure oscillation.

One can argue, however, that the used constraints derive from present measurements, and the likelihood of the scenario still hinges on possible atmospheric variations as we go backward in time. Nonetheless, the radiative cooling timescale exhibits a strong dependence on the equilibrium temperature of the atmosphere (${\sigma }_0\propto T_0^3$; Auclair-Desrotour et al. 2017a, Eq. (17)). As such, a warmer planet in the past would yield a shorter cooling timescale, and consequently, more efficient damping of the resonant amplitude (see Appendix L). On the other hand, atmospheric compositional variations can change the opacity parameter of the atmosphere in the visible and the infrared. An increase of the opacity in the visible to double its present value, that is, to α_A = 30%, can indeed place the response beyond the Precambrian threshold for all the considered typical values of σ₀. An increase to three times the present value of α_A is required to cross the threshold in the late Paleozoic/early Mesozoic when the longest dissipative timescales are considered. In contrast, an increase in atmospheric opacity in the infrared, which accompanies the abundance of Precambrian greenhouse gases (e.g., Catling & Zahnle 2020), delivers the opposite effect by attenuating the resonant tidal response, as we elaborate in Appendix L. Furthermore, the latter effect would also trigger the contribution of asynchronous tidal heating, which further attenuates the amplitude of the positive peak as we show in Sect. 3.2. Thus, with these analyses, it is unlikely that compositional variations in the past would have rendered the resonant thermotidal response certainly sufficient to cross the required threshold. This conclusion can be further regarded as conservative, since the employed linear model tends to overestimate the resonant amplification of the tidal response. This derives from the fact that, in the quasi-adiabatic regime, the model ignores the associated non-linearities of dissipative mechanisms. The remaining question is therefore: when did the Lamb resonance actually occur?

4.2 The spectral position of the Lamb resonance

The spectral position of the Lamb resonance, or equivalently, the geological time of its occurrence, is identified in the analytical model via the denominator highlighted in Eq. (22). The latter is a function of σ and is dependent on the planetary radius, gravitational acceleration, eigenvalue of the fundamental Hough mode, the radiative frequency, and the equilibrium surface temperature at the surface T_s, and is independent of σ_bl. Thus for the Earth, the resonance position is merely dependent on the equilibrium temperature at the surface and the radiative cooling frequency.

In Fig. 6, we plot the dependence of the spectral position of the resonance, in terms of LOD, on T_s. The apparent single solid curve of the dry adiabatic profile is actually a bundle of curves with different values of σ₀, but the effect of the latter is unnoticeable (if one varies σ₀ by two orders of magnitude, the resonant rotational period varies by few minutes). As such, the resonant frequency is predominantly controlled by T_s, which allows us to take the adiabatic limit of Eq. (22), and straightforwardly derive the tidal frequency that minimizes the denominator. In terms of the rotational period, the position of the resonance then reads:

$${\text{LOD}}_{\text{res}}=\frac{4\pi R_{\text{p}}}{\sqrt{{ℛ}_{\text{s}}{\text{Λ}}_n{T}_{\text{s}}}+2R_{\text{p}}{n}_{\star }}.$$(24)

In Appendix I, we compute the tidal solution in the limit of a moist adiabatic profile by modifying the thermodynamic equation to allow for the effect of humidity. There we show that the same expression of the resonant period holds in the moist adiabatic limit. In Appendix M, we derive the analogous expression in the limit of an isothermal atmospheric profile, and we show that the resonant period is reduced by almost one hour compared to that in the adiabatic limit when a fixed enthalpy is considered for an air column in the two limits. For all profiles, nonetheless, the resonant rotational period scales as the inverse square root of the surface equilibrium temperature. However, the evolution of the latter for the early Earth is widely debated. For instance, marine oxygen isotopes have been interpreted to indicate Archean ocean temperatures around 60–80° C (e.g., Knauth 2005; Robert & Chaussidon 2006). This interpretation is in contrast with geochemical analysis using phosphates (e.g., Blake et al. 2010), geological evidence of Archean glacial deposits (e.g., de Wit & Furnes 2016), geological carbon cycle models (e.g., Sleep & Zahnle 2001; Krissansen-Totton et al. 2018), numerical results of 3D GCMs (e.g., Charnay et al. 2017), and the fact that solar luminosity was 10–25% lower during the Precambrian (e.g., Charnay et al. 2020), altogether predicting a temperate climate and moderate temperatures throughout the Earth’s history.

We highlight with the gray shading on top of the curves modeled mean surface temperature variations adopted from Krissansen-Totton et al. (2018). As the latter temperature evolution is established in the time domain, we use the LOD evolution in Farhat et al. (2022) to map from time-dependence to LOD-dependence, and we further identify the corresponding geological eras of the LOD evolution with the color shadings. Given the present-day equilibrium surface temperature, the resonance occurs at LOD = 22.8 h. This value is in agreement with the 11.38 ± 0.16 h semi-diurnal period obtained by analyzing the spectrum of normal modes using pressure data on global scales (see Table 1 of Sakazaki & Hamilton 2020, first symmetric gravity mode of wavenumber k = −2). The resonant rotational period assuming an isothermal profile of the atmosphere is roughly one hour less than that in the neutrally stratified limit, placing it closer to 21.3 h estimate of Zahnle & Walker (1987) and Bartlett & Stevenson (2016). We emphasize here, however, that the resonant period does not exactly mark the period at which the thermotidal torque is maximum. The latter occurs at the peaks surrounding the resonance (see Figs. 2 and H.1), the difference between the two being dependent on the radiative cooling frequency.

Taking the LOD evolution model of Farhat et al. (2022) at face value, the temperature variations predicted in Krissansen-Totton et al. (2018) locate the resonance encounter in the early Mesozoic (or late Paleozoic if one considers wider temperature variations), and not in the Precambrian. In fact, for the resonance to be encountered in the Precambrian, even in the latest eras of it, the resonant period should move to less than ~21 h, but this requires an increase in the equilibrium temperature of at least 55°C, which is inconsistent with the studies mentioned above. Such an increase in temperature would also increase σ₀ by almost 19% (as we discuss in the previous section, ${\sigma }_0\propto T_0^3$; see also Appendix L), reducing the radiative cooling timescale and prompting more efficient damping of the tidal amplitude at the resonance. Moreover, such an increase in temperature would accompany increased greenhouse effects in the past, which in turn would increase the atmospheric absorption and thermotidal heating in the infrared. The latter would then place the Earth’s atmosphere in the regime of asynchronous thermotidal heating studied in Sect. 3.2, whereby the accelerative peak of the torque is further attenuated.

5 Summary and outlook

We were drawn to the problem of atmospheric thermal tides by the hypothesized scenario of a constant length of day on Earth during the Precambrian. Our motivation in investigating the scenario lies in its significant implications on paleoclimatic evolution, and the evident mismatch between LOD geological proxies and the predicted LOD evolution if this rotational equilibrium is surmised. The scenario hinges on the occurrence of a Lamb resonance in the atmosphere whereby an amplified thermotidal torque would cancel the opposing torque generated by solid and oceanic gravitational tides. Naturally, the atmospheric tidal torque is that of two flavors: it can either pump or deplete the rotational angular momentum budget of the planet, depending on the orientation of the generated tidal bulge.

With this rotational equilibrium scenario in mind, we have developed a novel analytical model that describes the tidal response of thermally forced atmospheres on rocky planets. The model derivation is based on the secure ground of the first principles of linear atmospheric dynamics, studied under classical approximations that are commonly drawn in earlier analytical works and in more recent numerical frameworks. The distinct feature that we imposed in this model is that of neutral atmospheric stratification, which presents a more realistic description of the Earth’s troposphere than the isothermal profile imposed in earlier analytical studies. In this limit, we derive from the model a closed-form solution of the tidal torque that can be efficiently used to study the evolution of planetary rotational dynamics. We accommodate into the model dissipative thermal radiation via linear Newtonian cooling, and turbulent and diffusive processes related to thermal inertia budgets in the boundary layer and the ground. As such, the model can be used to study a planetary thermotidal response when heated either by direct synchronous absorption of the incident stellar flux, or by a delayed infrared radiation from the ground.

We probed the spectral behavior of the tidal torque using this developed model in the two aforementioned limits. In the limit of longwave heating flux, the inherently delayed thermal response in the planetary boundary layer maneuvers the tidal bulge in such a way that, for typical values of thermal inertia in the ground and atmosphere, the accelerating effect of the tidal torque at the Lamb resonance is attenuated, and possibly annihilated. In the case of the Earth, where we apply the opposite limit of shortwave thermotidal heating and ignore the attenuating effect of asynchronous forcing, while the encounter of the resonance in the atmosphere is guaranteed, the epoch of its occurrence and the tidal amplitude it generates are uncertain. As such, we attempted a cautious incursion on constraining them and learned that:

• Assuming that temperate climatic conditions have prevailed over the Earth’s history, the resonance is likely to have occurred in the late Paleozoic/early Mesozoic, and not in the Precambrian. Unlike the Precambrian, this era is characterized by an amplified decelerating luni-solar gravitational torque;

• For judiciously constrained estimates of our atmospheric model parameters, the resonance does not amplify the accelerating thermotidal torque to a level comparable in magnitude to the gravitational counterpart in the said era.

It is important to note that these model predictions presume that thermotidal heating in the Earth has always been dominated by the shortwave. Compositional variations however, namely those associated with increased greenhouse contributions in the past would amplify the asynchronous thermotidal forcing in the longwave. The latter in turn, as we show in this work, further attenuates the accelerating flavor of the resonant torque. Exploring this end is certainly worthy of future efforts, but with the present indications at hand, we conclude that the occurrence of the rotational equilibrium is contingent upon a drastic increase in the Earth’s surface temperature (≥ 55°C), a long enough radiative cooling timescale (≥40 days), an increase in the shortwave flux opacity of the atmosphere, and that infrared thermotidal heating remained negligible in the past. We cannot completely preclude these requirements when considered separately, especially given the uncertainty in reconstructing the Earth’s temperature evolution in the Proterozoic. However, a warmer paleoclimate goes hand in hand with a shorter radiative cooling timescale, along with increased greenhouse gases that amplify the asynchronous thermotidal forcing. Both effects damp the accelerating flavor of the thermotidal torque. Put together, these indications suggest that the occurrence of the rotational equilibrium for the Earth is unlikely. To this end, future GCM simulations that properly model the Precambrian Earth to provide stringent constraints on our analytical predictions of the resonant amplification are certainly welcome.

Ultimately though, even if the locking into the resonance did not occur, the effect of the thermotidal torque at the resonance remains a robust and significant feature, and it should be accommodated in future modeling attempts of the Earth’s rotational evolution. Our model sets the table for efficiently studying such a complex interplay between several tidal players, both for the Earth and duly for its analogs. Interestingly, the question of the climatic response to the Lamb resonance, or similarly to oceanic tidal resonances, where abrupt and significant astronomical variations occur, largely remains an unexplored territory, perhaps requiring an armada of rigorous GCM simulations. This only leaves us with anticipated pleasure in weaving yet another thread in the rich tidal history of the Earth. Furthermore, we anticipate that the growing abundance of geological proxies, especially robust inferences associated with cyclostratigraphy, may help detect the whereabouts of these resonances and provide further constraints to our modeling efforts.

Acknowledgements

The authors are thankful to the referee, Dorian Abbot, for insightful comments which helped improve the manuscript considerably. M.F. expresses his gratitude to Kevin Heng for his hospitality at the LMU Munich Observatory where part of this work was completed. This work has been supported by the French Agence Nationale de la Recherche (AstroMeso ANR-19-CE31-0002-01) and by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Advanced Grant AstroGeo-885250). This work was granted access to the HPC resources of MesoPSL financed by the Region Île-de-France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche.

Appendix A The dimensionless governing equations of atmospheric tidal dynamics

Using the same definition of atmospheric variables as in the main text, we develop here the dimensionless governing system of equations describing tidal dynamics, which shall be used to recover the wave equation (2). We start with the classical primitive equations describing the atmospheric tidal response under the thin shell approximation (e.g., Siebert 1961; Chapman & Lindzen 1970), and we augment the heat transfer equation by the additional radiative term J_rad described in the main text. These equations read:

$${\partial }_t{V}_{\theta }-2\text{Ω\hspace{0.17em}}\cos \theta V_{\phi }=-\frac{1}{R_{\text{p}}}{\partial }_{\theta }\left(\frac{\delta p}{{\rho }_0}-U\right),$$(A.1) $${\partial }_t{V}_{\phi }+2\text{Ω\hspace{0.17em}}\cos \theta V_{\theta }=-\frac{1}{R_{\text{p}}\sin \theta }{\partial }_{\phi }\left(\frac{\delta p}{{\rho }_0}-U\right),$$(A.2) $$0=-{\partial }_z\delta p-g\delta \rho +{\rho }_0{\partial }_{z}U,$$(A.3) $$\frac{D\rho }{Dt}=-{\rho }_0\chi ,$$(A.4) $$\frac{1}{{\text{Γ}}_1}\frac{Dp}{Dt}-gH\frac{D\rho }{Dt}=\kappa {\rho }_{0}J-{\rho }_0{ℛ}_{\text{s}}{\sigma }_0\delta T,$$(A.5) $$\chi -{\partial }_z{V}_r={\nabla }_{\text{h}}\cdot V,$$(A.6) $$\frac{\delta p}{p_0}=\frac{\delta \rho }{{\rho }_0}+\frac{\delta T}{T_0}.$$(A.7)

This system describes atmospheric momentum conservation (Eqs. A.1 - A.3), mass conservation (Eq. A.4), heat transfer with Newtonian cooling (Eq. A.5), and the ideal gas law (Eq. A.7). Furthermore, we have adopted the traditional approximation (e.g., Unno et al. 1979), where we ignore the Coriolis coupling term between the vertical and horizontal parts of the momentum equations, thus allowing for the analytical treatment of the system. For the vertical momentum equation (A.3), the traditional approximation amounts to imposing the hydrostatic equilibrium for the tidal perturbation in addition to the background profiles, which are in turn averaged over the sphere. The hydrostatic approximation would allow us later to define the tidal torque as a function of the pressure anomaly at the surface.

In this system, the material (read Lagrangian) time derivative $\frac{D}{Dt}$ of any variable y is defined as

$$\frac{Dy}{Dt}={\partial }_{t}y+V_r\frac{d{y}_0}{dz},$$(A.8)

and we denote by χ the divergence of the velocity vector field V, and by the notation ∇_h · V its horizontal divergence such that

$${\nabla }_{\text{h}}\cdot V=\frac{1}{R_{\text{p}}\sin \theta }\left[{\partial }_{\theta }\left(\sin \theta V_{\theta }\right)+{\partial }_{\phi }{V}_{\phi }\right].$$(A.9)

We also introduce the calculation variable

$$G=-\frac{1}{{\text{Γ}}_1{p}_0}\frac{Dp}{Dt}.$$(A.10)

It is noteworthy here, as we impose neutral stratification on the atmosphere later in the model (Section 3.1), that one can argue that the gôp term in Eq.(A.3) is supposed to vanish since it is often referred to as the buoyancy term. However, this is not exactly the case since internal gravity waves are not the only source of density fluctuations. In fact, density fluctuations can also result from the planetary-scale compressibility waves, read Lamb modes. As such, the atmosphere can be at the same time neutrally stratified and still feature strong density variations across the horizontal direction. With this subtlety clarified, the tidal perturbations are then expanded in Fourier series of time and longitude, such that the tidal excitation frequency is denoted by σ. We introduce the reference velocity υ₀, and we non-dimensionalize all the variables; namely:

$$\begin{array}{l}dz=Hdx,t={\sigma }^{-1}\tilde{t},V=v_0\tilde{V},\delta p=p_0\tilde{\delta }p,\hfill \\ \delta \rho ={\rho }_0\tilde{\delta }\rho ,\delta T=T_0\tilde{\delta }T,U=v_0^2\tilde{U},G=\frac{\sigma }{{\text{Γ}}_1}\tilde{G},\hfill \\ {\nabla }_{\text{h}}=R_{\text{p}}^{-1}\tilde{\nabla }\cdot ,{\nabla }_{\text{h}}=R_{\text{p}}^{-1}\tilde{\nabla },J=\frac{\sigma gH}{\kappa }\tilde{J},\hfill \end{array}$$(A.11)

Under these definitions, the material derivative now reads

$$\frac{Dy}{Dt}=\sigma y_0\left[{\partial }_{\tilde{t}}\tilde{\delta }y+\frac{v_0}{\sigma H}\frac{d\ln y_0}{dx}{\tilde{V}}_r\right]\text{,}$$(A.12)

and the calculation variable $\tilde{G}$ becomes

$${\partial }_{\tilde{t}}\tilde{\delta }p+\frac{{\upsilon }_0}{\sigma H}\frac{d\ln p_0}{dx}{\tilde{V}}_r=-\tilde{G}.$$(A.13)

We next set υ₀ = σR_p and introduce the ratios of length scale η = R_p /H, and frequency α = σ/σ₀. Allowing for these changes of variables in the governing system of equations, one directly obtains a dimensionless system in the form:

$${\partial }_{\tilde{t}}{\tilde{V}}_{\theta }-v\cos \theta {\tilde{V}}_{\phi }=-{\partial }_{\theta }\left(\frac{gH}{{\upsilon }_0^2}\tilde{\delta }p-\tilde{U}\right),$$(A.14) $${\partial }_{\tilde{t}}{\tilde{V}}_{\phi }+v\cos \theta {\tilde{V}}_{\theta }=-\frac{1}{\sin \theta }{\partial }_{\phi }\left(\frac{gH}{{\upsilon }_0^2}\tilde{\delta }p-\tilde{U}\right),$$(A.15) $$\left({\partial }_x-1\right)\tilde{\delta }p=-\tilde{\delta }\rho +\frac{{\upsilon }_0^2}{gH}{\partial }_x\tilde{U},$$(A.16) $${\partial }_{\tilde{t}}\tilde{\delta }+\tilde{\nabla }\cdot \tilde{V}=-\eta \left({\partial }_x+\frac{d\ln {\rho }_0}{dx}\right){\tilde{V}}_r,$$(A.17) $$\frac{\tilde{G}}{{\text{Γ}}_1}+{\partial }_{\tilde{t}}\tilde{\delta }\rho +\eta \frac{d\ln {\rho }_0}{dx}{\tilde{V}}_r=-\tilde{J}+\frac{\delta \tilde{T}}{\alpha },$$(A.18) $$\tilde{\delta }p=\tilde{\delta }\rho +\tilde{\delta }T.$$(A.19)

Due to the periodic nature of the tidal forcing, the tidal response is Fourier-decomposed in time and longitude, allowing us to write all the physically varying quantities in the form

$$f\left(x,\theta ,\phi ,\tilde{t}\right)=\mathrm{Re}\left\{{\displaystyle \underset{m}{{\sum }^{\text{}}}}{f}^{m,v}\left(x,\theta \right)e^{i(\tilde{t}+m\phi )}\right\},$$(A.20)

where m is the longitudinal degree. Furthermore, due to the traditional approximation, the decoupling of the vertical and horizontal structures of tidal dynamics allows the expansion of the Fourier coefficients of Eq. (A.20) into series of Hough functions (Hough 1898) describing the horizontal tidal flow; namely:

$$f^{m,v}(x,\theta )={\displaystyle \underset{n}{{\sum }^{\text{}}}}{f}_n^{m,v}(x){\text{Θ}}_n^{m,v}(\theta ).$$(A.21)

Unless stated otherwise, we denote throughout the paper the coefficients $f_n^{m,v}$ by f_n for simplicity. As such, the Laplace tidal equation given in the main text by Eq. (1) is rewritten as:

$$ℒ{\text{Θ}}_n=-{\text{Λ}}_n{\text{Θ}}_n,$$(A.22)

where the operator ℒ is defined as

$$\begin{array}{l}ℒ=\frac{1}{\sin \theta }\frac{d}{d\theta }\left(\frac{\sin \theta }{1-v^2{\cos }^2\theta }\frac{d}{d\theta }\right)\hfill \\ -\frac{1}{1-v^2{\cos }^2\theta }\left(mv\frac{1+v^2{\cos }^2\theta }{1-v^2{\cos }^2\theta }+\frac{m^2}{{\sin }^2\theta }\right).\hfill \end{array}$$(A.23)

Appendix B Retrieving the wave equation of vertical dynamics

The governing system of dimensionless equations (A.14-A.19) decouples into equations describing the horizontal structure (A.14-A.15), and those describing the vertical structure. In the Fourier domain, the latter equations read:

$$i\tilde{\delta }{p}_n-\eta {\tilde{V}}_{r;n}=-{\tilde{G}}_n,$$(B.1) $$\left({\partial }_x-1\right)\tilde{\delta }{p}_n=-\tilde{\delta }{\rho }_n+\frac{v_0^2}{gH}{\partial }_x{\tilde{U}}_n,$$(B.2) $$\tilde{\delta }{\rho }_n-i\eta \left({\partial }_x+\frac{d\ln {\rho }_0}{dx}\right){\tilde{V}}_{r;n}=-ℒ\left(\frac{gH}{v_0^2}\tilde{\delta }{p}_n-{\tilde{U}}_n\right)\text{, }$$(B.3) $$\frac{{\tilde{G}}_n}{{\text{Γ}}_1}+i\tilde{\delta }{\rho }_n+\eta \frac{d\ln {\rho }_0}{dx}{\tilde{V}}_{r;n}=-{\tilde{J}}_n+\frac{\delta T_n}{\alpha },$$(B.4) $$\tilde{\delta }{p}_n=\tilde{\delta }{\rho }_n+\tilde{\delta }{T}_n.$$(B.5)

To recover the wave equation (2) from this system, we aim to reduce this system to a single second order partial differential equation in the calculation variable ${\tilde{G}}_n$. First, we combine Eq. (B.1) with the heat transfer equation (B.4) to eliminate the vertical component of velocity, ${\tilde{V}}_{r;n}$. Then we use the ideal gas law (B.5) to replace $\tilde{\delta }{T}_n$ by $\tilde{\delta }{p}_n$ and $\tilde{\delta }{\rho }_n$ in the resulting equation. It follows that

$$\tilde{\delta }{p}_n\left(1+\frac{d\ln H}{dx}-\frac{i}{\alpha }\right)-\tilde{\delta }{\rho }_n\left(1-\frac{i}{\alpha }\right)=\frac{i{N}_{\text{B}}^{2}H}{g}{\tilde{G}}_n-i{\tilde{J}}_n\text{.}$$(B.6)

Next, we use the hydrostatic equilibrium condition (B.2) to express $\tilde{\delta }{\rho }_n$ in terms of $\tilde{\delta }{p}_n$ and Ũ_n in the above equation to obtain

$${\partial }_x\tilde{\delta }{p}_n=i\frac{\alpha H{N}_{\text{B}}^2}{g(\alpha -i)}{\tilde{G}}_n-\frac{\alpha }{\alpha -i}\frac{d\ln H}{dx}\tilde{\delta }{p}_n-\frac{i\alpha }{\alpha -i}{\tilde{J}}_n+\frac{{\upsilon }_0^2}{gH}{\partial }_x{\tilde{U}}_n.$$(B.7)

Now our aim is to obtain another first order equation in the calculation variable ${\tilde{G}}_n$. We start with the hydrostatic equilibrium condition (B.2), and we substitute $\tilde{\delta }{p}_n$ by its expression (B.1). Then we use the continuity equation (B.3) to express $\tilde{\delta }{\rho }_n$ as a function of ${\tilde{V}}_{r;n}$ and $\tilde{\delta }{p}_n$, and we finally eliminate ${\tilde{V}}_{r;n}$ using equation (B.1) to obtain

$${\partial }_x{\tilde{G}}_n={\tilde{G}}_n-i\frac{gH}{v_0^2}ℒ{\tilde{\delta }}_n+i\left(ℒ-\frac{{\upsilon }_0^2}{gH}{\partial }_x\right){\tilde{U}}_n.$$(B.8)

We thus have two first order partial differential equations (B.7 - B.8) of the form

$${\partial }_x{\tilde{G}}_n=A_1{\tilde{G}}_n+B_1\tilde{\delta }{p}_n+C_1,$$(B.9) $${\partial }_x\tilde{\delta }p{p}_n=A_2{\tilde{G}}_n+B_2\tilde{\delta }{p}_n+C_2.$$(B.10)

Taking the derivative of Eq.(B.9), it is straightforward to obtain a second order equation of the form:

$$\left[{\partial }_x^2+A{\partial }_x+B\right]{\tilde{G}}_n=C_n,$$(B.11)

where

$$A(x)=-1+\frac{i}{\alpha -i}\frac{d\ln H}{dx},$$(B.12) $$B(x)=-\frac{\alpha }{\alpha -i}\left[\frac{i}{\alpha }\frac{d\ln H}{dx}+\frac{H^2{N}_{\text{B}}^2}{v_0^2}ℒ\right]\text{, }$$(B.13) $$C_n(x)=-\frac{gH\alpha }{v_0^2(\alpha -i)}ℒ{\tilde{J}}_n-[i\frac{d}{dx}\left(\frac{v_0^2}{gH}\frac{d}{dx}\right)$$(B.14) $$+\frac{1}{\alpha -i}\frac{d\ln H}{dx}\left(ℒ-\frac{{\upsilon }_0^2}{gH}\frac{d}{dx}\right)]{\tilde{U}}_n.$$(B.15)

Next, we implement a change of variable of the form ${\tilde{G}}_n(x)=\text{Φ}(x){\text{Ψ}}_n(x)$ to write Eq. (B.11) in the standard form of wave equations. We find that the first order term in Eq. (B.11) cancels when

$$2\frac{d\text{Φ}(x)}{dx}+A(x)\text{Φ}(x)=0,$$(B.16)

which is satisfied by

$$\text{Φ}(x)=\exp \left(-\frac{1}{2}{{{\displaystyle \int }}^{\text{}}}_0^{x}A\left(x^{\prime }\right)d{x}^{\prime }\right).$$(B.17)

This gives us the form of Φ(x) defined in Eq. (5), and consequently the wave equation (2), with a vertical wavenumber $\widehat{k}$ defined as

$$\begin{array}{l}{\widehat{k}}_n^2=-\frac{1}{4}\left(A^2+2\frac{dA}{dx}-4B\right)\hfill \\ =-\frac{1}{4}\{{\left(1-\frac{i}{\alpha -i}\frac{d\ln H}{dx}\right)}^2+2\frac{d}{dx}\left(\frac{i}{\alpha -i}\frac{d\ln H}{dx}\right)\hfill \\ -\frac{4\alpha }{\alpha -i}\left[\frac{gH{N}_{\text{B}}^2{\text{Λ}}_n}{v_0^2}-\frac{i}{\alpha }\frac{d\ln H}{dx}\right]\}.\hfill \end{array}$$(B.18)

With the definitions of the dimensionless control parameters given in Table 1, the wave number ${\widehat{k}}_n$ (Eq.B.18) and the forcing term (Eq.B.15) can be rewritten in the form given in the main text by Eqs. (3) and (6) respectively.

Appendix C Solutions of the vertical profiles of tidal variables

In this section we provide the solutions to the vertical profiles of the tidally varying scalar fields of pressure, density, and temperature, and the velocity vector field, given the solution of the wave equation Ψ(x)_n, or equivalently the variable ${\tilde{G}}_n$.

First, the expression of pressure variations is readily obtained from Eq.(B.8), which for a given Hough mode and using the definition of β (Table 1) reads

$$\tilde{\delta }{p}_n=\frac{1}{i\beta {\text{Λ}}_n}\left(\frac{d{\tilde{G}}_n}{dx}-{\tilde{G}}_n\right)+\frac{1}{\beta }\left(1+\frac{1}{\beta {\text{Λ}}_n}\frac{d}{dx}\right){\tilde{U}}_n.$$(C.1)

Given $\tilde{\delta }{p}_n(x)$, the horizontal component of the velocity field is straightforwardly obtained via (e.g., Eqs.(25-26) of Auclair-Desrotour et al. 2017b):

$${\tilde{V}}_{\theta ;n}=\frac{{\partial }_{\theta }+mv\cot \theta }{1-v^2{\cos }^2\theta }\left(\frac{gH}{{\upsilon }_0^2}\tilde{\delta }{p}_n-{\tilde{U}}_n\right),$$(C.2)

$${\tilde{V}}_{\phi ;n}=-\frac{v\cos \theta {\partial }_{\theta }+m{\sin }^{-1}\theta }{1-v^2{\cos }^2\theta }\left(\frac{gH}{v_0^2}\tilde{\delta }{p}_n-{\tilde{U}}_n\right),$$(C.3)

which gives

$${\tilde{V}}_{\theta ;n}=\frac{1}{{\text{Λ}}_n}\left(\frac{d{\tilde{G}}_n}{dx}-{\tilde{G}}_n+i{\beta }^{-1}\frac{d{\tilde{U}}_n}{dx}\right),$$(C.4)

$${\tilde{V}}_{\phi ;n}=\frac{i}{{\text{Λ}}_n}\left(\frac{d{\tilde{G}}_n}{dx}-{\tilde{G}}_n+i{\beta }^{-1}\frac{d{\tilde{U}}_n}{dx}\right).$$(C.5)

The vertical component of the velocity is established from equation (B.1), replacing $\tilde{\delta }{p}_n$ by its solution (C.1):

$${\tilde{V}}_{r;n}=\frac{1}{\eta \beta {\text{Λ}}_n}\left[\frac{d{\tilde{G}}_n}{dx}+\left(\beta {\text{Λ}}_n-1\right){\tilde{G}}_n\right]+\frac{i}{\eta \beta }\left(1+\frac{1}{\beta {\text{Λ}}_n}\frac{d}{dx}\right){\tilde{U}}_n.$$(C.6)

Finally, the solution for the density perturbation is obtained by replacing the solution of $\tilde{\delta }{p}_n$ (C.1) in equation (B.6) to get:

$$\begin{array}{l}\tilde{\delta }{\rho }_n=\frac{1}{i\beta {\text{Λ}}_n}\left(1+\frac{\kappa \alpha (\gamma -1)}{\alpha -i}\right)\left[\frac{d{\tilde{G}}_n}{dx}+\left(\frac{\alpha \beta \gamma \kappa \text{Λ}}{\alpha -i+\kappa \alpha (\gamma -1)}-1\right){\tilde{G}}_n\right]\hfill \\ +\frac{i\alpha {\tilde{J}}_n}{\alpha -i}+\frac{1}{\beta }\left(1+\frac{\kappa \alpha (\gamma -1)}{\alpha -i}\right)\left(1+\frac{1}{\beta {\text{Λ}}_n}\frac{d}{dx}\right){\tilde{U}}_n,\hfill \end{array}$$(C.7)

and the vertical profile of temperature is then readily obtained by the ideal gas law (B.5).

Appendix D Thermal exchange mechanisms and the total propagating flux

We define here explicitly the thermal exchange mechanisms used to establish the thermal energy budget at the surface and to construct the thermal forcing profile given by Eq.(10) of Section 2.2.

As mentioned in the main text, the incident stellar flux, δF_inc generates a variation δT_s in the surface temperature T_s. We denote the proportionality by the complex transfer function ${ℬ}_{\text{gr}}^{\sigma }$ such that $\delta T_{\text{s}}={ℬ}_{\text{gr}}^{\sigma }\delta F_{\text{inc}}$. The variation δT_s induces radiative emission with magnitude δF_rad. A fraction of the incident power, δQ_gr, is then transmitted to the ground by thermal conduction, and another fraction, δQ_atm, is transmitted to the atmosphere through turbulent thermal diffusion (see Figure D.1).

Having isolated these thermal mechanisms, the total power balance of the thermal perturbation at the surface interface is expressed as⁸:

$$\delta F_{\text{inc}}-\delta F_{\text{rad}}=\delta Q_{\text{gr}}+\delta Q_{\text{atm}}.$$(D.1)

To define each of these terms explicitly, we start with the surface radiative emission term δF_rad, which can be obtained by differentiating the Stefan-Boltzmann law as a function of the surface temperature assuming that the planetary surface radiates as a black body. Namely we have

$$\delta F_{\text{rad}}=4{\sigma }_{\text{SB}}{T}_{\text{s}}^3\delta T_{\text{s}},$$(D.2)

with σ_SB = 5.670373 × 10⁻⁸ W m⁻² K⁻⁴ being the Stephan-Boltzmann constant (Tiesinga et al. 2021).

Next, the terms associated with heat diffusion at the boundary interface, δQ_gr and δQ_atm, are derived using the gradient flux theory (Garratt 1994):

$$\delta Q_{\text{gr}}(t)={\frac{k_{\text{gr}}}{H}\frac{\partial \delta T}{\partial x}|}_{x=0^{-}},$$(D.3)

$$\delta Q_{\text{atm}}(t)=-{\frac{k_{\text{atm}}}{H}\frac{\partial \delta T}{\partial x}|}_{x=0^{+}}.$$(D.4)

Here, k_gr is the thermal conductivity of the ground, and k_atm is its atmospheric analog. We also denote by x = 0⁻ and x = 0⁺ the vicinity below and above the surface interface, respectively. Temperature variations near the surface can be traced by the heat transport equation, which in the frequency domain reads as (e.g., Chapman & Lindzen 1970)

$$i\sigma \delta T=\frac{K_{\text{gr}}}{H^2}\frac{{\partial }^2\delta T}{\partial x^2}\text{\hspace{0.17em}for\hspace{0.17em}}x<0,$$(D.5)

$$i\sigma \delta T=\frac{K_{\text{atm}}}{H^2}\frac{{\partial }^2\delta T}{\partial x^2}\text{\hspace{0.17em}for\hspace{0.17em}}x>0,$$(D.6)

with K_gr = k_gr/(ρ₀(0⁻)C_gr) and K_atm = k_atm/(ρ₀(0⁺)C_p) being the thermal diffusivities near the surface in the ground and the atmosphere respectively, C_gr and C_p being the respective thermal capacities per unit mass, while ρ₀(0⁻) and ρ₀(0⁺) are the densities in the vicinity of the surface interface. These equations have solutions of the form:

$$\delta T(x)=\delta T_{\text{s}}\exp \left\{(1+is)x/h_{\text{gr}}^{\sigma }\right\},\text{\hspace{0.17em}for\hspace{0.17em}}x\le 0,$$(D.7)

$$\delta T(x)=\delta T_{\text{s}}\exp \left\{-(1+is)x/h_{\text{atm}}^{\sigma }\right\},\text{\hspace{0.17em}for\hspace{0.17em}}x>0.$$(D.8)

Here, we remind the reader that s = sgn(σ), and we denote by $h_{\text{gr}}^{\sigma }$ and $h_{\text{atm}}^{\sigma }$ the skin thicknesses of heat transport by thermal diffusion in the ground and the atmosphere defined as

$$h_{\text{gr}}^{\sigma }=\frac{1}{H}\sqrt{\frac{2K_{\text{gr}}}{\left|\sigma \right|}}\text{\hspace{0.17em}and\hspace{0.17em}\hspace{0.17em}}{h}_{\text{atm}}^{\sigma }=\frac{1}{H}\sqrt{\frac{2K_{\text{atm}}}{\left|\sigma \right|}}.$$(D.9)

To obtain the explicit form of ${ℬ}_{\text{gr}}^{\sigma }$, we first rewrite the power budget balance of Eq.(D.1) as

$$\delta F_{\text{inc}}-{\frac{k_{\text{gr}}}{H}\frac{\partial \delta T}{\partial x}|}_{x=0^{-}}=-{\frac{k_{\text{atm}}}{H}\frac{\partial \delta T}{\partial x}|}_{x=0^{+}}+4{\sigma }_{\text{SB}}{T}_{\text{s}}^3\delta T_{\text{s}},$$(D.10)

and we substitute the solutions of δT(x) of Eq. (D.7) to obtain

$${ℬ}_{\text{gr}}^{\sigma }={\left\{\left(4{\sigma }_{\text{SB}}{T}_{\text{s}}^3\right)\left[1+(1+si)\sqrt{{\tau }_{\text{bl}}\left|\sigma \right|}\right]\right\}}^{-1}.$$(D.11)

Here we have defined the boundary layer timescale τ_bl as

$${\tau }_{\text{bl}}=\frac{1}{2}{\left(\frac{I_{\text{gr}}+I_{\text{atm}}}{4{\sigma }_{\text{SB}}{T}_{\text{s}}^3}\right)}^2,$$(D.12)

where the functions I_gr and I_atm denote, respectively, the thermal inertias of the ground and the atmospheric surface layers, and are specifically given by:

$$\begin{array}{c}{I}_{\text{gr}}={\rho }_0\left(0^{-}\right)C_{\text{gr}}\sqrt{K_{\text{gr}}},\\ I_{\text{atm}}={\rho }_0\left(0^{+}\right)C_{\text{p}}\sqrt{K_{\text{atm}}}.\end{array}$$(D.13)

Finally, with δQ_gr written as

$$\delta Q_{\text{gr}}=\frac{k_{\text{gr}}}{\sqrt{2K_{\text{gr}}}}\sqrt{\left|\sigma \right|}{ℬ}_{\text{gr}}^{\sigma }(1+si)\delta F_{\text{inc}},$$(D.14)

we are fully geared to define the total propagating flux feeding the atmosphere, δF_tot, as

$$\delta F_{\text{tot}}=\delta F_{\text{inc}}-\delta Q_{\text{gr}},$$(D.15)

from which we retrieve Eq.(10) of the main text. As for the numerical values of the material properties considered in the formalism above, we first set the gas heat capacity to its nominal value C_p = 1005 J/kg·K (Hilsenrath 1955). The volumetric heat capacity in the ground ρ₀(0⁻)C_gr varies between the land and the oceans, but since we are considering an effective value over a rocky planet, we adopt the typical value of 2.5 × 10⁶ J m⁻³ K⁻¹ (Van Wijk & De Vries 1963). The value for the density ρ₀ is dependent on the planetary setting understudy and is computed via ρ₀ = p₀/gH. The values for ground thermal conductivities and diffusivities vary among rock types and water content taking typical values of k_gr ~ 0.5 – 4 W m⁻¹ K⁻¹ and K_gr ~ 10⁻⁶ – 10⁻⁸ m² s⁻¹ (e.g., Andújar Márquez et al. 2016; Arkhangelskaya & Lukyashchenko 2018). In contrast, for the atmosphere, the turbulent eddy diffusivity K_atm is known to have a strong altitude dependence, but the profile of this dependence varies significantly in the literature between linear, exponentially decaying, and higher order polynomial profiles with several inflection points, with near surface estimates ranging between 0.1 and 100 m² s⁻¹ (e.g., Bernard 1962; Madsen 1977; Nieuwstadt 1983; Holtslag & Boville 1993; Berger & Grisogono 1998). Specifically, Bernard (1962) computes an average effective value over continental areas ${\overline{K}}_{\text{atm}}=3.6{\text{m}}^2{\text{s}}^{-1}$. This variability of the ground and atmospheric thermal parameters propagates to our model free parameter τ_bl, or equivalently the frequency σ_bl, via the thermal ground and atmospheric thermal inertias defined in Eq. (D.13). Thus in sweeping the limits of σ_bl in Figure 3, we are moving between strong and weak thermal inertia limits by self-consistently varying the values of the thermal properties. For instance, the value of σ_bl = 10⁻⁵ s⁻¹, for which the accelerating torque at the Lamb resonance completely vanishes, corresponds to k_gr = 1.25 W m⁻¹ K⁻¹, K_gr = 5 × 10⁻⁷ m² s⁻¹, and $K_{\text{atm}}={\overline{K}}_{\text{atm}}=3.6{\text{m}}^2{\text{s}}^{-1}$.

Appendix E The quadrupolar component of the incident stellar flux

As we are interested in the semi-diurnal tidal response, we decompose here the incoming stellar flux harmonically to isolate the quadrupolar component. We start by expressing the day-night periodically varying incident flux δF_inc (denoted hereafter by F for convenience) as:

$$F=\{\begin{array}{ll}{F}_{\star }\cos \text{Θ},\hfill & \text{for }0\le \text{Θ}\le \pi /2\hfill \\ 0,\hfill & \text{otherwise},\hfill \end{array}$$(E.1)

where Θ is the zenith angle. To obtain the quadrupolar component of the thermal forcing we first expand F in Legendre Polynomials:

$$F(\text{Θ})={\displaystyle \underset{l=0}{\overset{\infty }{{\sum }^{\text{}}}}}{F}_l{P}_l(\cos \text{\hspace{0.17em}Θ}),$$(E.2)

where the expansion coefficients are given by

$$\begin{array}{l}{F}_l=\frac{2l+1}{2}{{{\displaystyle \int }}^{\text{}}}_0^{\pi }F(\text{Θ})P_l(\cos \text{\hspace{0.17em}Θ})\sin \text{\hspace{0.17em}Θ}d\text{\hspace{0.17em}Θ},\hfill \\ =\frac{2l+1}{2}{F}_{\star }{{{\displaystyle \int }}^{\text{}}}_0^{\pi /2}{P}_l(\cos \text{\hspace{0.17em}Θ})\sin \text{\hspace{0.17em}Θ\hspace{0.17em}}\cos \text{\hspace{0.17em}Θ}d\text{Θ}.\hfill \end{array}$$(E.3)

Next, using the addition theorem (e.g., Abramowitz et al. 1988), we write the Legendre Polynomials as series of spherical harmonics:

$$P_l(\cos \text{\hspace{0.17em}Θ})=\frac{4\pi }{2l+1}{\displaystyle \underset{m=-l}{\overset{m=l}{{\sum }^{\text{}}}}}{Y}_{lm}(\theta ,\phi )Y_{lm}^{\ast }\left({\theta }_{\text{S}},{\phi }_{\text{S}}\right),$$(E.4)

where the asterisk corresponds to complex conjugation, and the coordinates (θ_S, φ_S) = [π/2, (Ω – n_⋆)t] define the position of the star when ignoring planetary obliquity. We thus have:

$$F={\displaystyle \underset{l=0}{\overset{\infty }{{\sum }^{\text{}}}}}{\displaystyle \underset{m=-l}{\overset{m=l}{{\sum }^{\text{}}}}}{C}_l{F}_l{Y}_{lm}(\theta ,\phi )Y_{lm}^{\ast }\left({\theta }_{\text{S}},{\phi }_{\text{S}}\right),$$(E.5)

where we have defined C_l = 4π/(2l + 1). Using the definition of spherical harmonics

$$Y_{lm}(\theta ,\phi )=N_{lm}{P}_{lm}(\cos \theta )e^{im\phi },$$(E.6)

with

$$N_{lm}=\sqrt{\frac{(2l+1)(l-m)!}{4\pi (l+m)!}},$$(E.7)

and the associated Legendre functions

$$P_{lm}(\cos \theta )={(-1)}^m{\left(1-{\cos }^2\theta \right)}^{\frac{m}{2}}\frac{d^m}{d{x}^m}{P}_l(\cos \theta ),$$(E.8)

we rewrite Eq.(E.5) as:

$$\begin{array}{l}F={\displaystyle \underset{l=0}{\overset{\infty }{{\sum }^{\text{}}}}}{\displaystyle \underset{m=-l}{\overset{m=l}{{\sum }^{\text{}}}}}{\gamma }_{lm}(\theta )\cos \left[m\left({\sigma }_{11}t+\phi \right)\right]\hfill \\ +i{\displaystyle \underset{l=1}{\overset{\infty }{{\sum }^{\text{}}}}}{\displaystyle \underset{m=1}{\overset{m=l}{{\sum }^{\text{}}}}}\sin \left[m\left({\sigma }_{11}t+\phi \right)\right]\left[{\gamma }_{lm}(\theta )-{\gamma }_{l-m}(\theta )\right],\hfill \end{array}$$(E.9)

where ${\gamma }_{lm}(\theta )=C_l{F}_l{N}_{lm}^2{P}_{lm}(\cos \theta )P_{lm}(0)$, and σ₁₁ = Ω – n_⋆. Next we note that

$$\begin{array}{l}{\gamma }_{lm}(\theta )-{\gamma }_{l-m}(\theta )=C_l{F}_l[N_{lm}^2{P}_{lm}(\cos \theta )P_{lm}(0)\hfill \\ -N_{l-m}^2{P}_{l-m}(\cos \theta )P_{l-m}(0)]=0,\hfill \end{array}$$

by virtue of the equality

$$P_{l-m}=\frac{(l-m)!}{(l+m)!}{P}_{lm}.$$(E.10)

Finally, F is written as the real part of the thermal forcing function in the Fourier domain:

$$F=\mathrm{Re}\left\{{\displaystyle \underset{l=0}{\overset{\infty }{{\sum }^{\text{}}}}}{\displaystyle \underset{m=0}{\overset{l}{{\sum }^{\text{}}}}}\left(2-{\delta }_{0m}\right){\gamma }_{lm}(\theta )e^{im({\sigma }_{{11}^t+\phi )}}\right\},$$(E.11)

from which we retrieve the quadrupolar component (l = m = 2) of the forcing in the form

$$F_{22}(\theta ,\phi ,t)=\frac{5}{64}{F}_{\star }{P}_{22}(\cos \theta )e^{i\left({\sigma }_{22}t+2\phi \right)},$$(E.12)

or in spherical harmonics as

$$F_{22}(\theta ,\phi ,t)=\frac{\sqrt{30\pi }}{16}{F}_{\star }{Y}_{22}(\theta ,\phi )e^{i{\sigma }_{22}t},$$(E.13)

where σ₂₂ = 2(Ω – n_⋆) is the semi-diurnal tidal forcing frequency.

Appendix F Energy carried by tidal waves

Essential to the closed-form solution we provide for the tidal response is the boundary condition imposed against tidal waves at the uppermost layer of the atmosphere. There we enforce a non-divergence condition on the tidal flow, namely, the energy of tidal flows 𝒲 should be bounded as x → ∞. It follows from Wilkes (1949) that the energy flow associated with the tidal (n, m, v)-mode across the vertical direction is given by

$$\mathcal{𝒲}=\frac{1}{2}{p}_0{v}_0\mathrm{Re}\left\{\tilde{\delta }p{\tilde{V}}_r^{\ast }\right\}.$$(F.1)

By substituting ${\tilde{V}}_r$ by its expression in Eq.(B.1), we notice that

$$\mathrm{Re}\left\{\tilde{\delta }p{\tilde{V}}_r^{\ast }\right\}={\eta }^{-1}\mathrm{Re}\left\{\tilde{\delta }p{\tilde{G}}^{\ast }\right\}.$$(F.2)

Therefore, using Eq.(7), we define the tidal energy flow as a function of the wave equation solution as

$$\mathcal{𝒲}=\frac{p_0{v}_0}{2\eta \beta \text{Λ}}Im\left\{\frac{d\tilde{G}}{dx}{\tilde{G}}^{\ast }\right\}+\frac{p_0{v}_0}{2\eta \beta }\mathrm{Re}\left\{\left(1+\frac{1}{\beta \text{Λ}}\frac{d}{dx}\right)\tilde{U}{\tilde{G}}^{\ast }\right\}.$$(F.3)

Ignoring the gravitational tidal potential Ũ, and using the variables Ψ and Φ of the wave equation in the main text, we rewrite the energy flow as

$$\mathcal{𝒲}=\frac{p_0{v}_0}{2\eta \beta \text{Λ}}Im\left\{\frac{d\text{Ψ}}{dx}{\text{Ψ}}^{\ast }-\frac{A}{2}|\text{Ψ}{|}^2\right\}|\text{Φ}{|}^2,$$

which scales as |Ψ|²|Φ|² as we discuss in the main text.

Appendix G The tidal torque exerted about the spin axis

At first glance, the definitions of the atmospheric tidal torque, as a function of the pressure anomaly variation, present in various works of the literature may seem inconsistent [see eg., Eq.(27) of Dobrovolskis & Ingersoll (1980); Eq.(3) of Zahnle & Walker (1987); Eq.(4) of Auclair-Desrotour et al. (2019)]. Upon a more careful inspection, one can relate the evident mismatch to the adopted convention in defining the pressure anomaly, which when accounted for reveals the overlap between the different expressions. Here we recover the expression of the tidal torque that we use in Eq.(21).

The instantaneous torque exerted by tidal force about the spin axis of a planet of volume 𝒱 can be expressed as

$$\mathcal{𝒯}={\displaystyle {\int }_{\mathcal{𝒱}}{\partial }_{\phi }{U}_{\text{T}}d\mu },$$(G.1)

where U_T(r, t) represents the tidal potential, r being the separation between the tidal players, and dµ being an infinitesimal mass parcel of the tidally induced mass redistribution, that is, the tidal bulge. Under the hydrostatic approximation, where the pressure and gravity forces are balanced, the distribution of mass can be quantified by the pressure variation instead of the density variation (e.g., Leconte et al. 2015). Furthermore, since the atmospheric layer is relatively thin, the mass distribution can be further quantified by the surface pressure anomaly. In the case of a thick atmospheric layer, however, one should integrate the mass distribution over the full radial extent of the layer. That said, the torque, averaged over an infinite time period P, is written as

$$\mathcal{𝒯}=\underset{P\to +\infty }{\lim }\frac{1}{P}{{{\displaystyle \int }}^{\text{}}}_0^P{\displaystyle {\int }_{\mathcal{𝒮}}{\partial }_{\phi }{U}_{\text{T}}\frac{R_{\text{p}}^2\delta p_{\text{s}}}{g}dSdt}.$$(G.2)

where 𝒮 is a sphere of unit radius, and dS = sin θdθdφ. Due to the periodic nature of the forcing and the response, the tidal potential and the pressure anomaly are expanded in Fourier series with σ being the tidal frequency such that

$$U_{\text{T}}(r,t)=\mathrm{Re}\left\{{\displaystyle \underset{\sigma }{{\sum }^{\text{}}}}{U}_{\text{T}}^{\sigma }(r)e^{i\sigma t}\right\},$$(G.3) $$\delta p_{\text{s}}(\theta ,\phi ,t)=\mathrm{Re}\left\{{\displaystyle \underset{\sigma }{{\sum }^{\text{}}}}\delta p_{\text{s}}^{\sigma }(\theta ,\phi )e^{i\sigma t}\right\}.$$(G.4)

As such, the torque is rewritten as

$$\mathcal{𝒯}=\frac{R_{\text{p}}^2}{2g}{\displaystyle \underset{\sigma }{{\sum }^{\text{}}}}\mathrm{Re}\left\{{\displaystyle {\int }_{\mathcal{𝒮}}{\partial }_{\phi }^{\ast }{U}_{\text{T}}^{\sigma \ast }\delta p_{\text{s}}^{\sigma }dS}\right\}.$$(G.5)

The Fourier coefficients of the tidal potential and the pressure anomaly are further expanded in spherical harmonics (Eq.E.6):

$$U_{\text{T}}^{\sigma }(r)={\displaystyle \underset{l=2}{\overset{+\infty }{{\sum }^{\text{}}}}}{\displaystyle \underset{m=-l}{\overset{l}{{\sum }^{\text{}}}}}{U}_{\text{T};l}^{m,\sigma }(r)Y_{l,m}(\theta ,\phi ),$$(G.6) $$\delta p_{\text{s}}^{\sigma }(\theta ,\phi )={\displaystyle \underset{l=2}{\overset{+\infty }{{\sum }^{\text{}}}}}{\displaystyle \underset{m=-l}{\overset{l}{{\sum }^{\text{}}}}}\delta p_{\text{s};l}^{m,\sigma }{Y}_{l,m}(\theta ,\phi ).$$(G.7)

Substituting these expansions in Eq.(G.5) and using the orthogonality property of spherical harmonics, we obtain:

$$\mathcal{𝒯}=\frac{R_{\text{p}}^2}{2g}{\displaystyle \underset{\sigma }{{\sum }^{\text{}}}}{\displaystyle \underset{l=2}{\overset{+\infty }{{\sum }^{\text{}}}}}{\displaystyle \underset{m=-l}{\overset{l}{{\sum }^{\text{}}}}}m\text{\hspace{0.17em}Im}\left\{U_{\text{T};l}^{m,\sigma }\delta p_{\text{s};l}^{m,\sigma }\right\}.$$

Considering that the harmonic terms of the tidal potential are components of a series in the ratio R_p/a_p, it suffices to truncate the series at the quadrupolar order (l = m = 2) to describe the semidiurnal response since R_p ≪ a_p. As such, noting that $g=G{M}_{\text{p}}/R_{\text{p}}^2$, and that the quadrupolar component of the potential is given by (e.g., Ogilvie 2014):

$$U_{\text{T};2}^{2,\sigma }=\sqrt{\frac{6\pi }{5}}\frac{G{M}_{\star }}{a_{\text{p}}}{\left(\frac{R_{\text{p}}}{a_{\text{p}}}\right)}^2,$$(G.8)

we straightforwardly obtain the expression of the torque in Eq.(21) (see also Leconte et al. 2015; Auclair-Desrotour et al. 2019):

$$\mathcal{𝒯}=\sqrt{\frac{6\pi }{5}}\frac{M_{\star }}{M_{\text{p}}}\frac{R_{\text{p}}^6}{a_{\text{p}}^3}\text{Im}\left\{\delta p_{\text{s};2}^{2,\sigma }\right\}.$$(G.9)

Appendix H The functional form of the tidal response

Our goal in this sections is to provide a more elaborate description of the functional form of the tidal solution obtained from our model. Namely, while Eq. (22) quantifies the imaginary part of the surface pressure anomaly – as required for computing the tidal torque – expressed in physical parameters, here we provide the full complex solution in a different form that facilitates future model modifications. First, the obtained expression of the complex surface pressure anomaly, from which we extract the imaginary part given by Eq. (22), is given by

$$\begin{array}{l}\tilde{\delta }{p}_{\text{s}}=\frac{{\alpha }_{\text{A}}\delta F_{22}}{p_s}\frac{\kappa g{\text{Λ}}_2}{R_{\text{p}}^2{\sigma }^3}\frac{\left[{\mathcal{𝒳}}_{\text{R}}+{\mathcal{𝒴}}_{\text{R}}\alpha +i(\mathcal{𝒳}\alpha +\mathcal{𝒴})\right]\alpha }{\left(1+2\zeta +2{\zeta }^2\right)}\hfill \\ \times {\left[{\left(\kappa -\beta {\text{Λ}}_2+1\right)}^2+{\alpha }^2{\left(\beta {\text{Λ}}_2-1\right)}^2\right]}^{-1},\hfill \end{array}$$(H.1)

where 𝒳 and 𝒴 are defined by Eq. (23), and 𝒳_R and 𝒴_R are defined by

$$\begin{array}{l}{𝒳}_{\text{R}}=\left(\kappa -\beta {\text{Λ}}_2+1\right)\left[2{\zeta }^2\left(1-{\mu }_{\text{gr}}\right)+\zeta \left(2-{\mu }_{\text{gr}}\right)+1\right],\hfill \\ {\mathcal{𝒴}}_{\text{R}}=s{\mu }_{\text{gr}}\zeta \left(\beta {\text{Λ}}_2-1\right).\hfill \end{array}$$(H.2)

Next, we reformulate this complex solution as

$$\tilde{\delta }{p}_{\text{s}}=K{\left(\frac{{\sigma }_{\text{L}}}{\sigma }\right)}^3\frac{\mathcal{𝒢}(\sigma )}{{\alpha }^{-1}\left(\kappa +1-\beta {\text{Λ}}_2\right)+i\left(\beta {\text{Λ}}_2-1\right)},$$(H.3)

where we introduced ${\sigma }_{\text{L}}={\sigma }_{\text{w}}\sqrt{{\text{Λ}}_2}$, the frequency which marks the peak of the modulus of the resonant pressure variation in the adiabatic limit, the transfer function of the thermal forcing 𝒢(σ) defined as

$$\mathcal{𝒢}(\sigma )=\frac{2{\zeta }^2\left(1-{\mu }_{\text{gr}}\right)+\zeta \left(2-{\mu }_{\text{gr}}\right)+1-is{\mu }_{\text{gr}}\zeta }{1+2\zeta +2{\zeta }^2},$$(H.4)

and the dimensionless constant factor K defined as

$$K=\frac{{\alpha }_{\text{A}}\delta F_{22}}{p_s}\frac{\kappa g{\text{Λ}}_2}{R_{\text{p}}^2{\sigma }_{\text{L}}^3}.$$(H.5)

The transfer function 𝒢 thus describes the delayed forcing of the atmosphere due to thermal inertia near the planetary surface. Therefore, this function reduces to 𝒢 (σ) = 1 in the case of direct absorption of the incident flux and the synchronous response of the surface (namely, in the limit of ζ = µ_gr = 0, which is the limit with which we describe the Earth in the main text). With 𝒢(σ) defined as such, the complex solution is structured in a convenient form where the effects of atmospheric dynamics and the thermal forcing are clearly separated. This facilitates future modifications of the model if one intends, for instance, to replace the thermal forcing function by another.

Next we define the normalized frequency X = σ/σ_L and the dissipation parameter ɛ = σ₀/σ_L ≪ 1, and we rewrite Eq. (H.3) as

$$\tilde{\delta }{p}_{\text{s}}=K\frac{\mathcal{𝒢}(X)}{\epsilon \left[(\kappa +1)X^2-1\right]+iX\left(X^2-1\right)}.$$(H.6)

The modulus, real part, and the imaginary part of this complex surface pressure oscillation are then straightforwardly deduced as:

$$\begin{array}{l}\left|\tilde{\delta }{p}_{\text{s}}\right|=K\frac{\left|\mathcal{𝒢}\right|}{\sqrt{{\epsilon }^2{\left[(\kappa +1)X^2-1\right]}^2+X^2{\left(X^2-1\right)}^2}},\hfill \\ \mathrm{Re}\left\{{\widetilde{\delta p}}_{\text{s}}\right\}=K\frac{\epsilon \left[(\kappa +1)X^2-1\right]{\mathcal{𝒢}}_{\text{R}}+X\left(X^2-1\right){\mathcal{𝒢}}_{\text{I}}}{{\epsilon }^2{\left[(\kappa +1)X^2-1\right]}^2+X^2{\left(X^2-1\right)}^2},\hfill \\ \text{Im}\left\{{\widetilde{\delta p}}_{\text{s}}\right\}=K\frac{\epsilon \left[(\kappa +1)X^2-1\right]{\mathcal{𝒢}}_{\text{I}}-X\left(X^2-1\right){\mathcal{𝒢}}_{\text{R}}}{{\epsilon }^2{\left[(\kappa +1)X^2-1\right]}^2+X^2{\left(X^2-1\right)}^2},\hfill \end{array}$$(H.7)

where 𝒢_R = Re{𝒢} and 𝒢_I = Im{𝒢}. In Figure H.1 we plot a sample of the tidal response in terms of the three components, the modulus, the imaginary, and the real part, centered at the resonance. Notable in the figure is the difference between the resonant amplification obtained for the absolute value of the pressure anomaly, or its real part, compared to its imaginary part. Thus one has to be careful in analyzing the tidal amplification induced by the resonance in terms of the surface pressure anomaly compared to its imaginary part, which goes into the computation of the tidal torque.

It is noteworthy that taking Eq. (H.7) in the low frequency limit (X → 0) and setting 𝒢 = 1 we obtain

$$\text{Im}\left\{\tilde{\delta }{p}_{\text{s}}\right\}=K\frac{X}{{\epsilon }^2+X^2},$$(H.8)

which is similar to the functional form of the imaginary part of the pressure oscillation given in Ingersoll & Dobrovolskis (1978), and obtained empirically using GCMs in Leconte et al. (2015). This provides the equivalence between the radiative cooling frequency used in our model and the dissipative frequency used in Leconte et al. (2015), as we did in Section 4.1.

Appendix I Moist adiabatic tidal model

In this section, we establish the tidal solution in the case of a moist adiabat rather than the dry adiabat assumed thus far. We show that the moist and dry adiabatic profiles exhibit the same functional form of the tidal response.

For a moist adiabatic temperature profile, the Brunt-Väisälä frequency is expressed as

$$N_{\text{B}}^2=\frac{g}{T_0}\left(\frac{\text{d}{T}_0}{\text{d}z}+\frac{g}{C_p+\text{Δ}{C}_p}\right),$$(I.1)

where ΔC_p designates the increase of the heat capacity per unit mass due to humidity. As a first approximation, this parameter can be considered to be uniform over the air column. This allows us to introduce the constant parameter κ_M, defined as

$${\kappa }_{\text{M}}\equiv \frac{R_{\text{s}}}{C_p+\text{Δ}{C}_p},$$(I.2)

equivalent to the κ parameter characterizing the dry adiabat. We also introduce the ratio q_M, such that 0 < q_M ≤ 1, defined as

$$q_{\text{M}}\equiv \frac{{\kappa }_{\text{M}}}{\kappa },$$(I.3)

which quantifies the effect of humidity on the atmospheric structure and tidal waves. Taking the ratio of the slopes of moist (5.9ºC/km) and dry (9.8ºC/km) adiabatic temperature profiles, we get q_M = 0.602 (e.g., Feiccabrino et al. 2015). The Brunt-Väisälä frequency is thus rewritten as

$$N_{\text{B}}^2=\frac{g}{H}\left(\frac{\text{d\hspace{0.17em}}\ln H}{\text{d}x}+{\kappa }_{\text{M}}\right).$$(I.4)

With these definitions, assuming neutral stratification (N_B = 0) leads to

$$H=H_{\text{s}}{\text{e}}^{-{\kappa }_{\text{M}}x}.$$(I.5)

The other background fields are deduced from H by making use of the hydrostatic approximation and the ideal gas law,

$$p_0(x)=p_{\text{s}}{\text{e}}^{-x},T_0(x)=T_{\text{s}}{\text{e}}^{-{\kappa }_{\text{M}}x},{\rho }_0(x)={\rho }_{\text{s}}{\text{e}}^{\left({\kappa }_{\text{M}}-1\right)x},$$(I.6)

with H_s = R_sT_s/g, and ρ_s = p_s/ (gH_s). Written as a function of the altitude, the normalized altitude coordinate x is expressed as

$$x=\frac{1}{{\kappa }_{\text{M}}}\ln \left(\frac{H_{\text{s}}}{H_{\text{s}}-{\kappa }_{\text{M}}z}\right)\text{,}$$(I.7)

and the atmospheric vertical profiles, given in Eq. (13) for the dry adiabat, in the case of a moist adiabat become

$$\begin{array}{ll}H(z)=H_{\text{s}}-{\kappa }_{\text{M}}z,\hfill & p_0(z)=p_{\text{s}}{\left(1-\frac{{\kappa }_{\text{M}}z}{H_{\text{s}}}\right)}^{\frac{1}{{\kappa }_{\text{M}}}},\hfill \\ T_0(z)=T_{\text{s}}-\frac{{\kappa }_{\text{M}}g}{R_{\text{s}}}z,\hfill & {\rho }_0(z)=\frac{p_{\text{s}}}{g{H}_{\text{s}}}{\left(1-\frac{{\kappa }_{\text{M}}z}{H_{\text{s}}}\right)}^{\frac{1-{\kappa }_{\text{M}}}{{\kappa }_{\text{M}}}}.\hfill \end{array}$$(I.8)

I.1 The wave equation and the tidal solution in the moist adiabatic limit

The humidity of the atmosphere only affects the thermodynamic equation through the modified heat capacity. The thermody-namic equation thus becomes

$$\left(q_{\text{M}}^{-1}-\kappa \right)\frac{\text{d}p}{\text{d}t}-q_{\text{M}}^{-1}gH\frac{\text{d}\rho }{\text{d}t}=\kappa {\rho }_{0}J-{\rho }_0{R}_{\text{s}}{\sigma }_0\delta T.$$(I.9)

By introducing the calculation variable G, we rewrite the above equation as

$$\left(q_{\text{M}}^{-1}-\kappa \right)G+q_{\text{M}}^{-1}\frac{gH}{p_0}\frac{\text{d}\rho }{\text{d}t}=-\kappa \frac{{\rho }_0}{p_0}J+\frac{{\rho }_0{R}_{\text{s}}}{p_0}{\sigma }_0\delta T.$$(I.10)

Finally, we normalize the equations and we switch to the frequency domain. This yields,

$$\frac{1}{\sigma {\rho }_0}\frac{\text{d}\rho }{\text{d}t}={\partial }_{\tilde{t}}\tilde{\delta }\rho +\frac{V_0}{\sigma H}\frac{\text{d\hspace{0.17em}}\ln {\rho }_0}{\text{d}x}{\tilde{V}}_r,$$(I.11)

and the normalized equation for $\tilde{G}$,

$$\left(1-\kappa q_{\text{M}}\right)\frac{\tilde{G}}{{\text{Γ}}_1}+i\tilde{\delta }\rho +\eta \frac{\text{d\hspace{0.17em}}\ln {\rho }_0}{\text{d}x}{\tilde{V}}_r=q_{\text{M}}\left(-\tilde{J}+\frac{\delta \tilde{T}}{\alpha }\right),$$(I.12)

with η = R_p/H and α = σ/σ₀. The structure factors in the tidal equations are expressed as

$$\frac{\text{d\hspace{0.17em}}\ln H}{\text{d}x}=\frac{H{N}_{\text{B}}^2}{g}-{\kappa }_{\text{M}},\frac{\text{d\hspace{0.17em}}\ln p_0}{\text{d}x}=-1,$$(I.13)

and

$$\frac{\text{d\hspace{0.17em}}\ln {\rho }_0}{\text{d}x}=-\left(1+\frac{\text{d}H}{\text{d}x}\right)=-\left[1-{\kappa }_{\text{M}}+\frac{H{N}_{\text{B}}^2}{g}\right].$$(I.14)

The horizontal structure of the tidal perturbation is obtained similarly as in the other tidal solutions. For the vertical structure, we follow the same steps as in Appendix B. First, the vertical velocity and the temperature variation are eliminated from the thermodynamic equation by making use of Eq. (B.1) and Eq. (B.5), respectively. It follows

$$\left[1+\frac{\text{d\hspace{0.17em}}\ln H}{\text{d}x}-i\frac{q_{\text{M}}}{\alpha }\right]\tilde{\delta }p-\left(1-i\frac{q_{\text{M}}}{\alpha }\right)\tilde{\delta }\rho =i\frac{H{N}_{\text{B}}^2}{g}\tilde{G}-i{q}_{\text{M}}\tilde{J}.$$(I.15)

Finally, we use Eq. (B.2) to express $\tilde{\delta }\rho$ as a function of $\tilde{\delta }p$ and $\tilde{U}$, and we end up with

$${\partial }_x\tilde{\delta }p=\frac{\alpha }{{\left(\alpha -i{q}_{\text{M}}\right)}^{-1}}\left[i\frac{H{N}_{\text{B}}^2}{g}\tilde{G}-\frac{\text{d\hspace{0.17em}}\ln H}{\text{d}x}\tilde{\delta }p-i{q}_{\text{M}}\tilde{J}\right]+\frac{V_0^2}{gH}{\partial }_x\tilde{U}.$$(I.16)

The equation for ${\partial }_x\tilde{G}$ is the same as in the dry adiabatic case and is given by Eq. (B.8). We thus reformulate the tidal equations as the first order system given by

$${\partial }_x\tilde{G}=A_1\tilde{G}+B_1\tilde{\delta }p+C_1,$$(I.17)

$${\partial }_x\tilde{\delta }p=A_2\tilde{G}+B_2\tilde{\delta }p+C_2,$$(I.18)

with the operators and coefficients

$$\begin{array}{ll}{A}_1=1,\hfill & A_2=i\frac{\alpha H{N}_{\text{B}}^2}{g\left(\alpha -i{q}_{\text{M}}\right)},\hfill \end{array}$$(I.19)

$$\begin{array}{ll}{B}_1=-i\frac{gH}{V_0^2}ℒ,\hfill & B_2=-\frac{\alpha }{\left(\alpha -i{q}_{\text{M}}\right)}\frac{\text{d\hspace{0.17em}}\ln H}{\text{d}x},\hfill \end{array}$$(I.20)

$$C_1=i\left(ℒ-\frac{V_0^2}{gH}{\partial }_x\right)\tilde{U},C_2=\frac{V_0^2}{gH}{\partial }_x\tilde{U}-\frac{i\alpha q_{\text{M}}\tilde{J}}{\alpha -i{q}_{\text{M}}}.$$(I.21)

This system is reduced to a single second order equation of $\tilde{G}$ in the form given by Eq. (B.11) but with the modified operators A, B, and C given by

$$A=-1+\frac{i{q}_{\text{M}}}{\alpha -i{q}_{\text{M}}}\frac{\text{d\hspace{0.17em}}\ln H}{\text{d}x},$$(I.22)

$$B=-\frac{\alpha }{\left(\alpha -i{q}_{\text{M}}\right)}\left[\frac{i}{\alpha }{q}_{\text{M}}\frac{\text{d\hspace{0.17em}}\ln H}{\text{d}x}+\frac{H^2{N}_{\text{B}}^2}{V_0^2}ℒ\right],$$(I.23)

$$\begin{array}{l}C=-\frac{gH\alpha q_{\text{M}}}{V_0^2\left(\alpha -i{q}_{\text{M}}\right)}ℒ\tilde{J}\hfill \\ -\left[\frac{q_{\text{M}}}{\alpha -i{q}_{\text{M}}}\frac{\text{d\hspace{0.17em}}\ln H}{\text{d}x}\left(ℒ-\frac{V_0^2}{gH}\frac{\text{d}}{\text{d}x}\right)+i\frac{\text{d}}{\text{d}x}\left(\frac{V_0^2}{gH}\frac{\text{d}}{\text{d}x}\right)\right]\tilde{U}.\hfill \end{array}$$(I.24)

The wave equation is finally written in its classical form by substituting ℒ = −Λ_n and introducing the change of variable $ℒ=-{\text{Λ}}_n$,

$$\frac{d^2\text{Ψ}}{d{x}^2}+k_n^2\text{Ψ}={\varphi }^{-1}C,$$(I.25)

with ϕ being expressed as

$$\varphi (x)=\exp \left[\frac{x}{2}-\frac{1}{2}{\displaystyle {\int }_0^x\frac{i{q}_{\text{M}}}{\alpha -i{q}_{\text{M}}}\left(\frac{H{N}_{\text{B}}^2}{g}-\kappa q_{\text{M}}\right)dx}\right],$$(I.26)

and the squared vertical wavenumber as

$$\begin{array}{l}{k}_n^2=-\frac{1}{4}\{{\left(1-\frac{i{q}_{\text{M}}}{\alpha -i{q}_{\text{M}}}\frac{\text{d\hspace{0.17em}}\ln H}{\text{d}x}\right)}^2+2\frac{\text{d}}{\text{d}x}\left(\frac{i{q}_{\text{M}}}{\alpha -i{q}_{\text{M}}}\frac{\text{d\hspace{0.17em}}\ln H}{\text{d}x}\right)\hfill \\ -\frac{4\alpha }{\alpha -i{q}_{\text{M}}}\left[\frac{gH{\text{Λ}}_n}{V_0^2}\left(\frac{\text{d\hspace{0.17em}}\ln H}{\text{d}x}+\kappa q_{\text{M}}\right)-i\alpha q_{\text{M}}\frac{\text{d\hspace{0.17em}}\ln H}{\text{d}x}\right]\},\hfill \end{array}$$(I.27)

Thus far, we have not implemented the neutral stratification assumption (N_B = 0). In doing so, and considering the thermal tide solely (Ũ = 0), the components of the wave equation simplify to

$$A=-\frac{\alpha +i{q}_{\text{M}}\left(q_{\text{M}}\kappa -1\right)}{\alpha -i{q}_{\text{M}}},$$(I.28)

$$B=\frac{i{q}_{\text{M}}^2\kappa }{\alpha -i{q}_{\text{M}}},$$(I.29)

$$C=\frac{g\alpha H{q}_{\text{M}}{\text{Λ}}_n}{V_0^2\left(\alpha -i{q}_{\text{M}}\right)}\tilde{J},$$(I.30)

$$\varphi (x)=\exp \left[-\frac{A}{2}x\right],$$(I.31)

$$k_n^2={\left[\frac{i\alpha +q_{\text{M}}+q_{\text{M}}\kappa }{2\left(\alpha -i{q}_{\text{M}}\right)}\right]}^2.$$(I.32)

It is now straightforward to obtain the tidal solution in the moist adiabatic limit following the same recipe we implemented in the dry adiabatic limit. For the vertical profile of the pressure perturbation we obtain

$$\begin{array}{l}\tilde{\delta }p(x)=\frac{iK{q}_{\text{M}}{\text{e}}^{q_{\text{M}}\kappa x}}{b_{\text{J}}\left(X-i\epsilon q_{\text{M}}\right)+i\epsilon q_{\text{M}}^2\kappa }\{\left(b_{\text{J}}+1\right){\text{e}}^{-b_{\text{J}}x}\hfill \\ +\frac{\left[1-\left(1+b_{\text{J}}\right)X^2\right]\left[X-i\epsilon q_{\text{M}}\left(1+q_{\text{M}}\kappa \right)\right]}{\left(X-i\epsilon q_{\text{M}}\right)\left(X^2-1\right)-i\epsilon q_{\text{M}}^2\kappa X^2}{\text{e}}^{\left(i{k}_n-\frac{A}{2}\right)x}\}.\hfill \end{array}$$(I.33)

with K being the dimensionless constant defined by Eq. (H.5) for a dry adiabat, and e and X are defined as in Appendix H. Thus the appearance of q_M in the common factor of the latter equation renders the dependence of the tidal response on humidity explicit. Finally, the surface pressure perturbation (x = 0) is readily deduced from the above expression in the form

$$\tilde{\delta }{p}_{\text{s}}=\frac{K{q}_{\text{M}}}{q_{\text{M}}\epsilon \left[\left(1+\kappa q_{\text{M}}\right)X^2-1\right]+iX\left(X^2-1\right)}.$$(I.34)

We observe that switching from the dry adiabat to the moist adiabat amounts to replacing everywhere κ by κ_M, and ɛ by q_Mɛ. The functional form of the solution is the same in both cases. Consequently, the resonant frequency in the moist adiabatic limit is exactly equal to that in the dry adiabatic limit, and so is the resonant LOD. Moreover, the normalizing factor at the numerator of Eq. (I.34) is reduced by the factor q_M as is the dissipation parameter ϵ. Therefore, the amplification factor which scales as the ratio of the two is barely affected by the humidity of the atmosphere.

Appendix J Estimates of the absorbed tidal energy

In this section, we compute estimates of the atmospheric opacity parameter α_A. As defined in the main text (Section 2.2), α_A quantifies the fraction of thermal energy that is absorbed by the atmosphere and is actually available for semidiurnal thermotidal dynamics. To our knowledge, this input tidal power is rarely evaluated in the literature. Here we show that it can be estimated from the distributions of heating rates that are measured in the Earth’s atmosphere or computed using GCM simulations. We refer in what follows to Chapman & Lindzen (1970) and Vichare & Rajaram (2013) who provide such distributions for water vapor (H₂O) and ozone (O₃).

First, we consider the vertical and latitudinal profiles shown by Chapman & Lindzen (1970). These profiles do not directly quantify the heating rate J, but rather the equivalent temperature variation τ_G, which approximately corresponds to the temperature oscillation produced by J^m,σ. The input power is thus defined as (Chapman & Lindzen 1970, Eqs. (143-144))

$$J=\mathrm{Re}\left\{{\displaystyle \sum _{\sigma }}{\displaystyle \sum _m}{J}^{m,\sigma }(x,\theta ,\phi )e^{i\sigma t}\right\},$$(J.1)

with the coefficients J^m,σ given by

$$J^{m,\sigma }(x,\theta ,\phi )=i\sigma C_{\text{p}}{\tau }_{\text{G}}^m\sigma (x,\theta ,\phi ),$$(J.2)

where the heat capacity of the gas per unit mass C_p = ℛ_s/κ. Figure 3.2 of Chapman & Lindzen (1970) shows the normalized vertical distribution of ${\tau }_{\text{G}}^{2,\sigma }$. Following the authors (their Eq. 163), we set H = 7.6 km. The semidiurnal distributions of ${\tau }_{\text{G}}^{2,\sigma }$ for water vapor and ozone can be approximated by the functions

$${\tau }_{\text{G},{\text{H}}_2\text{O}}(x,\theta ,\phi )=T_{{\text{H}}_2\text{O}}{V}_{{\text{H}}_2\text{O}}(x)\sin \theta e^{i2\phi },$$(J.3)

$${\tau }_{\text{G},{\text{O}}_3}(x,\theta ,\phi )=T_{{\text{O}}_3}{V}_{{\text{O}}_3}(x){\sin }^2\theta e^{i2\phi },$$(J.4)

with $T_{{\text{H}}_2\text{O}}=0.035\text{\hspace{0.17em}\hspace{0.17em}K}$ and $T_{{\text{O}}_3}=0.25\text{\hspace{0.17em}\hspace{0.17em}K}$. As the vertical profile of water vapor shown by Figure 3.2 of Chapman & Lindzen (1970) is an exponentially decaying function of the altitude, the function $V_{{\text{H}}_2\text{O}}$ of Eq.(J.3) is set to

$$V_{{\text{H}}_2\text{O}}(x)=e^{-k_{{\text{H}}_2\text{O}}x},$$(J.5)

where $k_{{\text{H}}_2\text{O}}\approx 1/3$. Similarly, the vertical profile of ozone absorption can be approximated by a truncated sine function,

$$V_{{\text{O}}_3}\left(x\right)=\{\begin{array}{ll}{K}_{{\text{O}}_3}\sin \left(\frac{\pi \left(x-x_1\right)}{x_2-x_1}\right),\hfill & x_1\le x\le x_2,\hfill \\ 0\hfill & \text{elsewhere},\hfill \end{array}$$(J.6)

with $K_{{\text{O}}_3}=1.4$, x₁ = z₁/H, and x₂ = z₂/H, the altitudes z₁ and z₂ being set to z₁ = 18 km and z₂ = 78 km. At a given location on the sphere (θ, φ), the total flux absorbed by the atmosphere for the spherical harmonic Y₂₂ is defined as

$$\delta ℱ(\theta ,\phi )=\delta {ℱ}_{22}{Y}_{22}(\theta ,\phi ).$$(J.7)

The efficiency parameter α_A of our model is basically the ratio between δℱ₂₂ and the corresponding component of the incident Solar flux, which we retrieve in Appendix E (Eq. E.13) as

$$\delta F_{22}=\frac{\sqrt{30\pi }}{16}{F}_{\star }.$$(J.8)

In order to simplify calculations, we assume that the latitudinal profile of ozone can be approximated by the function sin² θ, similar to the latitudinal profile of water vapor. Under this assumption, we rewrite Eqs.(J.3) and (J.4) as

$${\tau }_{\text{G},{\text{H}}_2\text{O}}(x,\theta ,\phi )=\frac{4}{3}\sqrt{\frac{6\pi }{5}}{T}_{{\text{H}}_2\text{O}}{V}_{{\text{H}}_2\text{O}}(x)Y_{22}(\theta ,\phi ),$$(J.9)

$${\tau }_{\text{G},{\text{O}}_3}(x,\theta ,\phi )=\frac{4}{3}\sqrt{\frac{6\pi }{5}}{T}_{{\text{O}}_3}{V}_{{\text{O}}_3}(x)Y_{22}(\theta ,\phi ),$$(J.10)

and the corresponding heating rates as

$$J_{{\text{H}}_2\text{O}}(x,\theta ,\phi )=\frac{4}{3}\sqrt{\frac{6\pi }{5}}\sigma C_{\text{p}}{T}_{{\text{H}}_2\text{O}}{V}_{{\text{H}}_2\text{O}}(x)Y_{22}(\theta ,\phi ),$$(J.11)

$$J_{{\text{O}}_3}(x,\theta ,\phi )=\frac{4}{3}\sqrt{\frac{6\pi }{5}}\sigma C_{\text{p}}{T}_{{\text{O}}_3}{V}_{{\text{O}}_3}(x)Y_{22}(\theta ,\phi ),$$(J.12)

where the imaginary number appearing in Eq.(J.2) has been ignored since it does not affectthe amplitude of the thermal forc-ing⁹. For a given component, the total input power is obtained by integrating the heating rate over the air column

$$\delta ℱ={\displaystyle {\int }_0^{+\infty }J{\rho }_{0}Hdx=\frac{p_{\text{s}}}{g}}{\displaystyle {\int }_0^{+\infty }J{e}^{-x}dx},$$(J.13)

noting that ρ₀H = (p₀/g). This requires computing the integral of the vertical distribution of the heating rates,

$$A_{{\text{H}}_2\text{O}}={\displaystyle {\int }_0^{+\infty }{V}_{{\text{H}}_2\text{O}}{e}^{-x}dx=\frac{1}{1+k_{{\text{H}}_2\text{O}}}},$$(J.14)

$$\begin{array}{l}{A}_{{\text{O}}_3}={\displaystyle {\int }_0^{+\infty }{V}_{{\text{O}}_3}{e}^{-x}dx}\hfill \\ =K_{{\text{O}}_3}\pi \frac{\left(x_2-x_1\right)\left(e^{-x_1}+e^{-x_2}\right)}{{\pi }^2+{\left(x^2-x_1\right)}^2}.\hfill \end{array}$$(J.15)

We end up with $A_{{\text{H}}_2\text{O}}=0.75$ and $A_{{\text{O}}_3}=4.5\times {10}^{-2}$. Therefore, the total fluxes absorbed by water vapor and ozone are given by

$$\delta F_{{\text{H}}_2\text{O}}(\theta ,\phi )=\delta {ℱ}_{{\text{H}}_2\text{O};22}{Y}_{22}(\theta ,\phi ),$$(J.16)

$$\delta F_{{\text{O}}_3}(\theta ,\phi )=\delta {ℱ}_{{\text{O}}_3;22}{Y}_{22}(\theta ,\phi ),$$(J.17)

with the components

$$\delta {ℱ}_{{\text{H}}_2\text{O};22}=\frac{4}{3}\sqrt{\frac{6\pi }{5}}\sigma C_{\text{p}}\frac{p_{\text{s}}}{g}{T}_{{\text{H}}_2\text{O}}{A}_{{\text{H}}_2\text{O}},$$(J.18)

$$\delta {ℱ}_{{\text{O}}_3;22}=\frac{4}{3}\sqrt{\frac{6\pi }{5}}\sigma C_{\text{p}}\frac{p_{\text{s}}}{g}{T}_{{\text{O}}_3}{A}_{{\text{O}}_3}.$$(J.19)

Using σ = 1.45 × 10⁻⁴ s⁻¹, C_p = 1004.7 J kg⁻¹ K⁻¹, p_s = 1013.25 hPa, g = 9.81 m s⁻², and F_⋆ = 1366 W m⁻², we obtain δF₂₂ = 828.8 W m⁻², and

$$\begin{array}{ll}\delta {ℱ}_{{\text{H}}_2\text{O};22}\approx 102\text{\hspace{0.17em}W\hspace{0.17em}}{\text{m}}^{-2},\hfill & \frac{\delta {ℱ}_{{\text{H}}_2\text{O};22}}{\delta F_{22}}\approx 12.3\text{\%},\hfill \\ \delta {ℱ}_{{\text{O}}_3;22}\approx 44\text{\hspace{0.17em}W\hspace{0.17em}}{\text{m}}^{-2},\hfill & \frac{\delta {ℱ}_{{\text{O}}_3;22}}{\delta F_{22}}\approx 5.3\text{\%},\hfill \\ \delta {ℱ}_{\text{tot};22}\approx 146\text{\hspace{0.17em}W\hspace{0.17em}}{\text{m}}^{-2},\hfill & \frac{\delta {ℱ}_{\text{tot};22}}{\delta F_{22}}\approx 17.6\%.\hfill \end{array}$$(J.20)

Similar estimates can be computed from the profiles shown by Vichare & Rajaram (2013) (see their Figure 4), where the heating rate is plotted as a function of altitude for the main Hough functions structuring the semi-diurnal atmospheric tide on Earth. Without elucidating again the copious details of the procedure, we obtain for the profiles in Vichare & Rajaram (2013) the following:

$$\begin{array}{ll}\delta {ℱ}_{{\text{H}}_2\text{O};22}\approx 116\text{\hspace{0.17em}W\hspace{0.17em}}{\text{m}}^{-2},\hfill & \frac{\delta {ℱ}_{{\text{H}}_2\text{O};22}}{\delta F_{22}}\approx 14.0\text{\%},\hfill \\ \delta {ℱ}_{{\text{O}}_3;22}\approx 34\text{\hspace{0.17em}W\hspace{0.17em}}{\text{m}}^{-2},\hfill & \frac{\delta {ℱ}_{{\text{O}}_3;22}}{\delta F_{22}}\approx 4.2\text{\%},\hfill \\ \delta {ℱ}_{\text{tot};22}\approx 150\text{\hspace{0.17em}W\hspace{0.17em}}{\text{m}}^{-2},\hfill & \frac{\delta {ℱ}_{\text{tot};22}}{\delta F_{22}}\approx 18.1\%.\hfill \end{array}$$(J.21)

Put together, these computed estimates are consistent and suggest that the efficiency parameter is around α_A~17–18%, which is close to the constraint on α_A (14%, Section 4.1) obtained by fitting our model parameters to reproduce the present-day semidiurnal surface pressure variations.

Appendix K Benchmarking our model against GCM simulations

The work of Auclair-Desrotour et al. (2019) presents, to our knowledge, the only attempt on recovering the tidal torque exerted about the spin axis of a thermally forced atmosphere in the high frequency regime, namely, the regime where the Lamb resonance is encountered, using GCM simulations. Specifically, the mentioned work computed the tidal torque from 3D simulations of atmospheric dynamics using the generic version of the LMDZ GCM (Hourdin et al. 2006) in the case of a nitrogen dominated atmosphere around a dry rocky Venus-like planet.

Imposing on our analytical solution (Eq.22) the same planetary parameters defined in Table 1 of Auclair-Desrotour et al. (2019), for a 10 bar atmosphere, and setting σ_bl = 10⁻⁵ s⁻¹ and σ₀ = 4 × 10⁻⁶ s⁻¹, we recover in Figure K.1 a tidal spectrum that reasonably fits the GCM computed spectrum. In particular, we recover the two amplified responses of the torque in the vicinity of synchronization and at the Lamb resonance, and we recover the asymmetry in the Lamb resonance where the positive peak of the torque is annihilated. Furthermore, we recover a similar scaling of the torque with the forcing frequency σ. Namely, for the GCM spectrum, the tidal torque 𝒯 ∝ σ^0.75 near synchronization, and in the high frequency regime 𝒯 ∝ σ⁻¹; while for our model, 𝒯 ∝ σ^0.53 in the former regime, and 𝒯 ∝ σ⁻¹ in the latter.

Appendix L Atmospheric compositional variations and the possible interplay between shortwave and longwave heating

When studying the case of the Earth in the main text, we have assumed that the atmosphere is always heated by the shortwave. We have thus ignored the asynchronous heating in the infrared, which, as we showed in Section 3.2, breaks the symmetry of the Lamb resonance by attenuating the accelerative peak of the tidal torque. The results obtained for the amplitude of the accelerative tidal torque at the resonance can thus be considered as the optimistic limit of the resonant amplitude, as the latter can be attenuated once the longwave heating starts to play a significant role. However, an additional attenuation effect is delivered by the more efficient radiative cooling on a warmer Earth. Namely, the radiative cooling frequency features a strong temperature dependence that scales as σ₀ ∝ T³. Henceforth, the temperature increase required to place the resonance occurrence in the Precambrian would also attenuated the resonant amplitude of the torque by increasing σ₀. Here we intend to elaborate further on this signature, which is driven by atmospheric compositional variations, by studying the functional form of the tidal response obtained in Appendix H.

The form of the pressure anomaly in Eqs. (H.7) allows for an analytical investigation of the Lamb resonance. Namely, we can characterize the positive and negative peaks of the resonance by studying the derivative of the imaginary part of the pressure anomaly. We focus on the limit studied for the Earth in the main text and set 𝒢 = 1, that is, we assume thermotidal heating is dominated by direct shortwave absorption. Computing the roots of the derivative of Eq. (H.7) and substituting them in the expression, we obtain the amplitude of the peaks around the resonance. Namely, we denote by η the ratio between the value of the tidal torque reached at the peaks around resonance to that at present. Considering the predominant terms of the obtained roots, this amplification ratio reads as

$$\eta \approx \frac{K_{\oslash }}{K_{\oplus }}\frac{{\epsilon }_{\oplus }}{{\epsilon }_{\oslash }}{\eta }_{\oplus },$$(L.1)

where we have denoted by the subscript ⊘ the parameter values at the resonance, in contrast with the present values denoted by ⊕.

Now we allow for variations in K that are induced by the parameters entering into the definition of Eq. (H.5). Namely, presuming that κ, R_p, g, and Λ₂ are indeed constant, the rest of the parameters would modify the expression of the amplification factor ŋ (Eq. L.1) into:

$$\eta ={\eta }_{\oplus }\left(\frac{{\alpha }_{\text{A};\oslash }}{{\alpha }_{\text{A};\oplus }}\right)\left(\frac{F_{\star ;\oslash }}{F_{\star :\oplus }}\right)\left(\frac{p_{\text{s};\oplus }}{p_{\text{s};\oslash }}\right){\left(\frac{{\sigma }_{\text{L};\oplus }}{{\sigma }_{\text{L};\oslash }}\right)}^2\left(\frac{{\sigma }_{0;\oplus }}{{\sigma }_{0;\oslash }}\right).$$(L.2)

The resonant frequency in the adiabatic limit scales as ${\sigma }_{\text{L}}\propto T_{\text{s}}^{1/2}$. Moreover, the radiative cooling frequency is also a function of the surface temperature and it scales as ${\sigma }_0\propto p_{\text{s}}^{-1}{T}_{\text{s}}^3$ (e.g., Eq. 10 of Showman & Guillot 2002). Including these scaling relationships in the latter expression of the amplification factor one finds that

$$\eta ={\eta }_{\oplus }\left(\frac{{\alpha }_{\text{A};\oslash }}{{\alpha }_{\text{A};\oplus }}\right)\left(\frac{F_{\star ;\oslash }}{F_{\star :\oplus }}\right){\left(\frac{T_{\text{s};\oplus }}{T_{\text{s};\odot }}\right)}^4.$$(L.3)

The equilibrium surface temperature varies as T_s ∝ (α_hF_⋆)^1/4, where α_h is the fraction of flux that is absorbed by the atmosphere, mainly in the infrared absorption. It follows that

$$\eta ={\eta }_{\oplus }\left(\frac{{\alpha }_{\text{A};\oslash }}{{\alpha }_{\text{A};\oplus }}\right)\left(\frac{{\alpha }_{\text{h};\oplus }}{{\alpha }_{\text{h};\oslash }}\right).$$(L.4)

This interesting and simplified version of the amplification factor shows that η depends on the opacity of the atmosphere in the visible (via α_A) and in the infrared (via α_h), which are functions of the composition. Namely, the two absorption coefficients have opposite effects on η: a past atmosphere that is more opaque in the visible features a greater amplification factor, and thus prompts a stronger tidal response upon the resonance occurrence. In contrast, a larger past opacity in the infrared tends to attenuate the resonance amplified response. Thus the bottom line is, increased greenhouse gases in the past (mainly CO₂), tend to decrease the amplitude of the accelerative resonant response of the atmosphere via, first, Eq. (L.4), and via increasing the contribution of the asynchronous thermotidal heating that induces a phase lag on the tidal bulge and consequently attenuates the accelerative torque around the resonance.

Appendix M The resonant period in the isothermal limit

Here we discuss the limit where the temperature profile of the atmosphere is modeled by an isotherm, in contrast with the adi-abatic profile we assume in our neutrally stratified model. The isothermal limit translates to a positive and constant value of the Brunt–Väisälä frequency N_B, and it has been studied earlier using the linear theory of atmospheric tides (e.g., Lindzen 1978; Lindzen & Blake 1972; Auclair-Desrotour et al. 2019). We are specifically interested in estimating the difference between the two models in terms of the Lamb resonance period. Namely, we aim to define the spectral position of the resonance in the isothermal limit and compare with what we obtained for our model in Eq. 24.

We adopt the analytical computations in Auclair-Desrotour et al. (2019) that yield a frequency-dependence for the pressure surface anomaly in the form (their Eq. 12):

$$\text{Im}\left\{\delta p_{\text{s}}\right\}=p_{\text{s}}\frac{\kappa J_{\text{s}}}{gH\sigma }\frac{\frac{H}{h}\left(b_{\text{J}}+\frac{1}{2}+\kappa \right)-\frac{1}{2}\left(b_{\text{J}}+1\right)}{\left[b_{\text{J}}\left(b_{\text{J}}+1\right)+\frac{\kappa H}{h}\right]\left(\frac{H}{h}-\frac{1}{{\text{Γ}}_1}\right)}.$$(M.1)

Here the parameters and variables are equivalent in definition to those we use in our work, and h is the equivalent depth. This equation defines the position of the Lamb resonance by the equivalent depth h_L = Γ₁ H (see also Lindzen & Blake 1972). Hence, introducing the definition of the equivalent depth

$$h_n=\frac{R_{\text{p}}^2{\sigma }^2}{{\text{Λ}}_{n}g},$$(M.2)

we obtain, for the fundamental mode, the resonant frequency and consequently the resonant rotational period in the form:

$${\text{LOD}}_{\text{res; iso}}=\frac{4\pi R_{\text{p}}}{\sqrt{{ℛ}_{\text{s}}{\text{Λ}}_2{\text{Γ}}_1\overline{T}}+2R_{\text{p}}{n}_{\star }}.$$(M.3)

This equation differs from that in the neutrally stratified limit (Eq. 24) in the fact that temperature here cannot be assumed to be the surface temperature T_s, and the Γ₁ factor multiplied by the temperature. To compare the two equations, we assign for the isothermal model a constant value of temperature corresponding to a mass-average of the temperature profile in the neutrally stratified case assuming a fixed total enthalpy of the air column. Namely, using Eq. 13, we obtain $\overline{T}=T_{\text{s}}/(\kappa +1)$, thus we rewrite Eq. M.3 as:

$${\text{LOD}}_{\text{res;\hspace{0.17em}so}}=\frac{4\pi R_{\text{p}}}{\sqrt{{ℛ}_{\text{s}}{\text{Λ}}_n{T}_{\text{s}}/\left(1-{\kappa }^2\right)}+2R_{\text{p}}{n}_{\star }}.$$(M.4)

This equation yields a resonant rotational period that is roughly one hour smaller than that obtained in the neutrally stratified case; therefore, the curve shown in Figure 6 shifts vertically downward by ~ 1 hr. This places the resonance for present conditions around 21.6 hr, which is closer to the 21.3 hr value obtained by Zahnle & Walker (1987).

References

---

All Tables

All Figures