Issue 
A&A
Volume 684, April 2024



Article Number  A49  
Number of page(s)  23  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202348625  
Published online  04 April 2024 
Thermal tides in neutrally stratified atmospheres: Revisiting the Earth’s Precambrian rotational equilibrium
^{1}
IMCCE, CNRS, Observatoire de Paris, PSL University, Sorbonne Université,
77 Avenue DenfertRochereau,
75014
Paris, France
email: mohammad.farhat@obspm.fr
^{2}
School of Earth and Ocean Sciences, University of Victoria,
Victoria, British Columbia, Canada
Received:
15
November
2023
Accepted:
15
January
2024
Rotational dynamics of the Earth, over geological timescales, have profoundly affected local and global climatic evolution, probably contributing to the evolution of life. To better retrieve the Earth’s rotational history, and motivated by the published hypothesis of a stabilized length of day during the Precambrian, we examined the effect of thermal tides on the evolution of planetary rotational motion. The hypothesized scenario is contingent upon encountering a resonance in atmospheric Lamb waves, whereby an amplified thermotidal torque cancels out the opposing torque generated by the oceans and solid interior, driving the Earth into rotational equilibrium. With this scenario in mind, we constructed an ab initio model of thermal tides on rocky planets describing a neutrally stratified atmosphere. The model takes into account dissipative processes with Newtonian cooling and diffusive processes in the planetary boundary layer. We retrieved, from this model, a closedform solution for the frequencydependent tidal torque, which captures the main spectral features previously computed using 3D general circulation models. In particular, under longwave heating, diffusive processes near the surface and the delayed thermal response of the ground prove to be responsible for attenuating, and possibly annihilating, the accelerating effect of the thermotidal torque at the resonance. When applied to the Earth, our model prediction suggests the occurrence of the Lamb resonance in the Phanerozoic, but with an amplitude that is insufficient for the rotational equilibrium. Interestingly, though our study was motivated by the Earth’s history, the generic tidal solution can be straightforwardly and efficiently applied in exoplanetary settings.
Key words: Earth / planets and satellites: atmospheres / planets and satellites: dynamical evolution and stability / planets and satellites: terrestrial planets / planet–star interactions
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1 Introduction
For presentday Earth, the semidiurnal 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 lunisolar 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, coplanarity, 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; AuclairDesrotour & Leconte 2018; Gu et al. 2019; Lee 2020), and nearsynchronous and Earthlike rocky exoplanets (Cunha et al. 2015; Leconte et al. 2015; AuclairDesrotour et al. 2017a, AuclairDesrotour et al. 2019). Namely, for planetary systems within the socalled habitable zone, the gravitational tidal torque diminishes in the regime near spinorbit 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 EarthMoon system (Farhat et al. 2022), where we focused on modeling tidal dissipation in the Earth’s paleooceans 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; AuclairDesrotour 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; AuclairDesrotour 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 AuclairDesrotour 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 piecewise 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 twofold: first, to provide a novel physical model for the planetary thermotidal torque that better matches the GCMcomputed 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.
Fig. 1 Modeled histories of the rotational motion of the Earth. The Earth’s LOD evolution is plotted in time over geological timescales for three models: (i) the model of Farhat et al. (2022), where the evolution is driven solely by oceanic and solid tidal dissipation; and (ii) the model of Zahnle & Walker (1987), where the Lamb resonance is encountered for LOD~21 h, forcing a rotational equilibrium on the Earth; (iii) the model of Bartlett & Stevenson (2016), which also adopts the equilibrium scenario, but further studies the effect of thermal noise, and the required temperature variation to escape the equilibrium. Three curves of the latter model correspond to different parameterizations of the gravitational tide. Plotted on top of the modeled histories are geological proxies of the LOD evolution which can be retrieved from http://astrogeo.eu. 
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 threedimensional 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_{0} + δp, T = T_{0} + δT, ρ = ρ_{0} + δρ, and 𝒱 = V_{0} + 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^{1} (Lindzen & McKenzie 1967), where radiative losses, J_{rad}, are parameterized by the characteristic frequency σ_{0}; namely J_{rad} = p_{0}σ_{0}/(κρ_{0}T_{0})δT, where κ = (Γ_{1} − 1)/Γ_{1} = 0.285 and Γ_{1} is the adiabatic exponent. Similar to Leconte et al. (2015), we associate with σ_{0} a radiative cooling timescale τ_{rad} = 4τ/σ0.
2.1 The vertical structure of tidal dynamics
We are interested in providing a closedform solution for the frequency^{2} 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 nontheoreticians 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 socalled 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 eigenvalueeigenfunction problem defined as Laplace’s tidal equation (Laplace 1798; Lee & Saio 1997):
where the set of Hough functions serves as the solution (Hough 1898), 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 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 . To construct these sets of eigenfunctionseigenvalues 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 wavelike equation that describes vertical thermotidal dynamics and reads as:
Here, as is the common practice (e.g., Siebert 1961; Chapman & Lindzen 1970), we use the reduced altitude as the vertical coordinate, where the pressure scale height H(z) = ℛ_{s}T_{0}(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 is defined via
By virtue of the nondimensionalization of the governing system of equations, dimensionless control parameters appear in the wavenumber definition. Namely, , and , where we have introduced the characteristic frequency , 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:
We define the Brunt–Väisälä frequency in the limit of an isothermal atmospheric profile, . 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
and the forcing function C_{n}(x) which reads:
Hereafter, we identify the dimensionless form of the tidal variables by the tilde diacritic. Namely in Eq. (6), the dimensionless thermal forcing function , while the gravitational analog , where the reference velocity (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:
where , and the calculation variable .
Dimensionless control parameters determining the regime of the atmospheric tidal response and defined throughout the text.
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 (or equivalently in dimensional form, J_{n}). We adopt a vertical tidal heating profile of the form
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
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 , 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 σ_{0}, and the atmospheric opacity parameter α_{A}. In analogy to a = σ/σ_{0} we define the dimensionless parameter for the boundary layer .
By virtue of the power budget balance established in Appendix D, we define the total propagating flux δF_{tot} as
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 thermotidal 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 Venuslike planets, as is the fiducial setting studied using GCMs in Leconte et al. (2015) and AuclairDesrotour et al. (2019).
Finally, as we are interested in the semidiurnal tidal response, we decompose the thermal forcing in Appendix E to obtain the amplitude of the quadrupolar component as , where , L_{⋆} being the stellar luminosity, and a_{p} the starplanet 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; AuclairDesrotour et al. 2019), in which case the scale height H becomes independent of the vertical coordinate, and by virtue of Eq. (4), .
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 . In fact, can be expressed in terms of the potential temperature Θ_{0} (e.g., Sect. 2.10 of Vallis 2017):
whereby the stability of the atmosphere is controlled by the slope of Θ_{0}. 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 buoyancyrestored 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^{3}
It also follows that the background profiles of the scalar variables read as (AuclairDesrotour et al. 2017a):
We thus obtain for the heating profile (using Eqs. (8), (9), and (10))
As such, the wave Eq. (2) is rewritten as
where the complex functions 𝒜_{n} and ℬ are defined as:
The wave Eq. (15) admits the general solution
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 . Accordingly, the nondivergence of the flow condition is satisfied if one sets c_{2} = 0 and takes the proper sign of the wavenumber (Eq. (12)), namely:
– The second condition is the natural wall condition imposed by the ground interface, which enforces . 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_{1} in the form:
Under these boundary conditions, we are now fully geared to analytically compute the solution of the wave equation, Ψ_{n}(x) (or equivalently ), but we are specifically interested in retrieving a closedform solution of the quadrupolar tidal torque. The latter takes the general form^{4} (Appendix G):
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 . 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_{2} = 0 and c_{1} 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:
where we have defined the complex functions 𝒳 and 𝒴 as
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 humiditydependent analogs.
Before embarking on any results, we pause here for a few remarks on the provided closedform 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 thermotidal 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 Sunlike 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 Λ_{2}. Noting that we have dropped the superscripts, we remind the reader that for the semidiurnal (m = 2) response, , thus A is frequencydependent 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 Λ_{2} 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 closedform 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 βΛ_{2} = 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 Venuslike dry planet (M_{p} = 0.815M_{⊕}, R_{p} = 0.95 R_{⊕}, a_{p} = 0.73 au, g = 8.87 m s^{−2}), with a 10 bar atmosphere and a scale height at the surface H_{0} = 10 km, thermally forced by a solarlike 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 spinorbit 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 GCMgenerated spectra when varying the planetary setup (e.g., AuclairDesrotour 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.
Fig. 2 Spectrum of semidiurnal atmospheric thermal tides. Plotted is the imaginary part of the normalized pressure anomaly (; Eq. (22)) as a function of the normalized forcing frequency ω = (Ω − n_{⋆})/n_{⋆} = σ/2n_{⋆}, where n_{⋆} is the mean motion of the stellar perturber. The planetarystellar parameters are those of the fiducial planetary system defined in Sect. 3.1. 
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 (AuclairDesrotour 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^{−2} s^{−1}, 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^{−5} s^{−1} 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 AuclairDesrotour 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 massbased 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 smallamplitude nonmigrating tidal components they produce (e.g., Schindelegger & Ray 2014).
Fig. 3 (A)symmetry of the Lamb resonance. Similar to Fig. 2, the imaginary part of the normalized pressure anomaly (Eq. (22)), which is associated with the semidiurnal tide, is plotted as a function of the normalized forcing frequency ω = (Ω − n_{⋆})/n_{⋆} = σ/2n_{⋆}, for the same planetarystellar parameters. We focus here on the high frequency regime around the Lamb resonance. Different panels correspond to different values of σ_{bl}, or different thermal inertias in the ground and the atmosphere. Allowing for thermal inertia results in a delayed ground response, of which the signature is clear in inducing an asymmetry in the spectral behavior around the resonance. 
Fig. 4 Thermally induced tidal bulge revealed. Shown are polar snapshots of the radial and longitudinal variations of the tidal pressure anomaly δp(x) in the equatorial plane. The snapshots are shown from a top view, and the troposphere is puffed in size by virtue of the used massbased vertical coordinate ϛ = p/p_{s}. The longitudinal axes are shown in increments of 30° with 0° at the substellar point, while the radial axes are in increments of 0.25. The profile of the pressure perturbation is also normalized by the exponentially decaying pressure background profile. Snapshots are taken at different spectral positions that cover the passage through the Lamb resonance, which specifically occurs at ω = 262.6. In the top row, the response describes the limit of a planet with a synchronous atmospheric absorption, mimicking the Earth’s direct absorption by ozone and water vapor, and it shows the continuous movement of the bulge, function of ω, from lagging to leading the substellar point. In contrast, in the bottom row, and for the prescribed value of σ_{bl}, the delayed response of the ground forces the bulge to always lag the substellar point, thus acting to decelerate the planetary rotation. 
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 closedform 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 modeldependent. 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 Precambrian 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 σ_{0} 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 σ_{0} ≳ 10^{−4}s^{−1}, we severely suppress the amplified pressure response around the resonance. Conversely, for values of σ_{0} ≲ 10^{−6.5} s^{−1}, 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^{16} Nm) is roughly a quarter of the present gravitational torque (~4.51 × 10^{16} Nm) (e.g., Zahnle & Walker 1987; Farhat et al. 2022), thus requiring, via Eq. (21), Im{δp_{s}} on the order of 880 Pa^{5}. 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 σ_{0} and α_{A} can place the Earth’s atmosphere in the identified regions. Stringent constraints on σ_{0} are hard to obtain for the Earth since σ_{0} is an effective parameter that in reality is dependent on altitude. Furthermore, in the linear theory of tides, we are forced to ignore the layertolayer 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 kdistributions (e.g., Lacis & Oinas 1991) as performed in Leconte et al. (2015) and in AuclairDesrotour 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 σ_{0} 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^{6}. 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_{2}) 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 semidiurnal 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 semidiurnal 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_{2} constraint we obtain. Finally, it is also noteworthy how the plotted S_{2} constraint is insensitive to variations in σ_{0} over a wide interval, which prohibits the determination of the present value of σ_{0} 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^{7} 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 (; AuclairDesrotour 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 σ_{0}. 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 quasiadiabatic regime, the model ignores the associated nonlinearities of dissipative mechanisms. The remaining question is therefore: when did the Lamb resonance actually occur?
Fig. 5 Parametric study of the tidal response. Plotted is a contoured surface of the amplitude of the imaginary part of the positive semidiurnal pressure anomaly at the Lamb resonance, over a grid of values of our free parameters σ_{0} and α_{A}. The solid black isoline marks the level curve of Im{δp_{s}} = 880 Pa, and defines from below a region in (α_{A}, σ_{0})space where the thermotidal response is sufficient to cancel the gravitational counterpart in the Precambrian. Analogously, the dashed isoline defines the threshold needed in the late Paleozoic/early Mesozoic, 350−250 Ma. The horizontal shaded area corresponds to typical values of the radiative cooling rate as described in the main text. The other shaded area defines the region of parameter space that yields the presently observed semidiurnal tidal bulge. The gray area on the left covers the parametric region where the resonance features a lower pressure amplitude than the present. 
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 σ_{0}, but the effect of the latter is unnoticeable (if one varies σ_{0} 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:
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; KrissansenTotton 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 KrissansenTotton 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 timedependence to LODdependence, and we further identify the corresponding geological eras of the LOD evolution with the color shadings. Given the presentday equilibrium surface temperature, the resonance occurs at LOD = 22.8 h. This value is in agreement with the 11.38 ± 0.16 h semidiurnal 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 KrissansenTotton 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 σ_{0} by almost 19% (as we discuss in the previous section, ; 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.
Fig. 6 Dependence of the resonant rotational period on the mean surface temperature. By virtue of Eq. (24), the LOD at which the Lamb resonance occurs scales as the inverse square root of the mean surface temperature (solid curve). In Appendix I, we show that the same dependence holds if one considers a moist adiabatic profile. The shaded area of temperature variations highlights 95% confidence intervals for the past temperature evolution according to the carbon cycle model of KrissansenTotton et al. (2018). The identified geological eras correspond to the LOD evolution model of Farhat et al. (2022). The overlap between the modeled temperature evolution and the black curve places the resonance occurrence in the late Paleozoic/early Mesozoic. 
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 closedform 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 lunisolar 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 ANR19CE31000201) and by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Advanced Grant AstroGeo885250). This work was granted access to the HPC resources of MesoPSL financed by the Region ÎledeFrance and the project Equip@Meso (reference ANR10EQPX2901) 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:
(A.1) (A.2) (A.3) (A.4) (A.5) (A.6) (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 of any variable y is defined as
and we denote by χ the divergence of the velocity vector field V, and by the notation ∇_{h} · V its horizontal divergence such that
We also introduce the calculation variable
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 planetaryscale 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 υ_{0}, and we nondimensionalize all the variables; namely:
Under these definitions, the material derivative now reads
and the calculation variable becomes
We next set υ_{0} = σR_{p} and introduce the ratios of length scale η = R_{p} /H, and frequency α = σ/σ_{0}. Allowing for these changes of variables in the governing system of equations, one directly obtains a dimensionless system in the form:
(A.14) (A.15) (A.16) (A.17) (A.18) (A.19)
Due to the periodic nature of the tidal forcing, the tidal response is Fourierdecomposed in time and longitude, allowing us to write all the physically varying quantities in the form
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:
Unless stated otherwise, we denote throughout the paper the coefficients by f_{n} for simplicity. As such, the Laplace tidal equation given in the main text by Eq. (1) is rewritten as:
where the operator ℒ is defined as
Appendix B Retrieving the wave equation of vertical dynamics
The governing system of dimensionless equations (A.14A.19) decouples into equations describing the horizontal structure (A.14A.15), and those describing the vertical structure. In the Fourier domain, the latter equations read:
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 . First, we combine Eq. (B.1) with the heat transfer equation (B.4) to eliminate the vertical component of velocity, . Then we use the ideal gas law (B.5) to replace by and in the resulting equation. It follows that
Next, we use the hydrostatic equilibrium condition (B.2) to express in terms of and Ũ_{n} in the above equation to obtain
Now our aim is to obtain another first order equation in the calculation variable . We start with the hydrostatic equilibrium condition (B.2), and we substitute by its expression (B.1). Then we use the continuity equation (B.3) to express as a function of and , and we finally eliminate using equation (B.1) to obtain
We thus have two first order partial differential equations (B.7  B.8) of the form
Taking the derivative of Eq.(B.9), it is straightforward to obtain a second order equation of the form:
where
Next, we implement a change of variable of the form 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
which is satisfied by
This gives us the form of Φ(x) defined in Eq. (5), and consequently the wave equation (2), with a vertical wavenumber defined as
With the definitions of the dimensionless control parameters given in Table 1, the wave number (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 .
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
Given , the horizontal component of the velocity field is straightforwardly obtained via (e.g., Eqs.(2526) of AuclairDesrotour et al. 2017b):
Fig. D.1 Schematic of the bichromatic radiative and diffusive flux exchange mechanisms at the surface interface between the lower troposphere and the ground. The power balance between the fluxes in given by Eq. (D.1). 
which gives
The vertical component of the velocity is established from equation (B.1), replacing by its solution (C.1):
Finally, the solution for the density perturbation is obtained by replacing the solution of (C.1) in equation (B.6) to get:
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 such that . 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^{8}:
To define each of these terms explicitly, we start with the surface radiative emission term δF_{rad}, which can be obtained by differentiating the StefanBoltzmann law as a function of the surface temperature assuming that the planetary surface radiates as a black body. Namely we have
with σ_{SB} = 5.670373 × 10^{−8} W m^{−2} K^{−4} being the StephanBoltzmann 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):
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)
with K_{gr} = k_{gr}/(ρ_{0}(0^{−})C_{gr}) and K_{atm} = k_{atm}/(ρ_{0}(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}(0^{−}) and ρ_{0}(0^{+}) are the densities in the vicinity of the surface interface. These equations have solutions of the form:
Here, we remind the reader that s = sgn(σ), and we denote by and the skin thicknesses of heat transport by thermal diffusion in the ground and the atmosphere defined as
To obtain the explicit form of , we first rewrite the power budget balance of Eq.(D.1) as
and we substitute the solutions of δT(x) of Eq. (D.7) to obtain
Here we have defined the boundary layer timescale τ_{bl} as
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:
Finally, with δQ_{gr} written as
we are fully geared to define the total propagating flux feeding the atmosphere, δF_{tot}, as
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}(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^{6} J m^{−3} K^{−1} (Van Wijk & De Vries 1963). The value for the density ρ_{0} is dependent on the planetary setting understudy and is computed via ρ_{0} = p_{0}/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^{−1} K^{−1} and K_{gr} ~ 10^{−6} – 10^{−8} m^{2} s^{−1} (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^{2} s^{−1} (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 . 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 selfconsistently varying the values of the thermal properties. For instance, the value of σ_{bl} = 10^{−5} s^{−1}, for which the accelerating torque at the Lamb resonance completely vanishes, corresponds to k_{gr} = 1.25 W m^{−1} K^{−1}, K_{gr} = 5 × 10^{−7} m^{2} s^{−1}, and .
Appendix E The quadrupolar component of the incident stellar flux
As we are interested in the semidiurnal tidal response, we decompose here the incoming stellar flux harmonically to isolate the quadrupolar component. We start by expressing the daynight periodically varying incident flux δF_{inc} (denoted hereafter by F for convenience) as:
where Θ is the zenith angle. To obtain the quadrupolar component of the thermal forcing we first expand F in Legendre Polynomials:
where the expansion coefficients are given by
Next, using the addition theorem (e.g., Abramowitz et al. 1988), we write the Legendre Polynomials as series of spherical harmonics:
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:
where we have defined C_{l} = 4π/(2l + 1). Using the definition of spherical harmonics
with
and the associated Legendre functions
we rewrite Eq.(E.5) as:
where , and σ_{11} = Ω – n_{⋆}. Next we note that
by virtue of the equality
Finally, F is written as the real part of the thermal forcing function in the Fourier domain:
from which we retrieve the quadrupolar component (l = m = 2) of the forcing in the form
or in spherical harmonics as
where σ_{22} = 2(Ω – n_{⋆}) is the semidiurnal tidal forcing frequency.
Appendix F Energy carried by tidal waves
Essential to the closedform 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 nondivergence 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
By substituting by its expression in Eq.(B.1), we notice that
Therefore, using Eq.(7), we define the tidal energy flow as a function of the wave equation solution as
Ignoring the gravitational tidal potential Ũ, and using the variables Ψ and Φ of the wave equation in the main text, we rewrite the energy flow as
which scales as Ψ^{2}Φ^{2} 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 AuclairDesrotour 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
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
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
As such, the torque is rewritten as
The Fourier coefficients of the tidal potential and the pressure anomaly are further expanded in spherical harmonics (Eq.E.6):
Substituting these expansions in Eq.(G.5) and using the orthogonality property of spherical harmonics, we obtain:
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 , and that the quadrupolar component of the potential is given by (e.g., Ogilvie 2014):
we straightforwardly obtain the expression of the torque in Eq.(21) (see also Leconte et al. 2015; AuclairDesrotour et al. 2019):
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
where 𝒳 and 𝒴 are defined by Eq. (23), and 𝒳_{R} and 𝒴_{R} are defined by
Next, we reformulate this complex solution as
where we introduced , 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
and the dimensionless constant factor K defined as
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.
Fig. H.1 Sample of the tidal response, plotted in terms of the modulus, real part, and imaginary part of the surface pressure anomaly, all given by Eqs. (H.7). The function f(x) on the yaxis indicates the functional form of the pressure variation, that is, the pressure function divided by the constant factor K. The xaxis is centered at the position of the Lamb resonance. Notable is the difference in the resonant amplification of the response between the real part and the imaginary part. 
Next we define the normalized frequency X = σ/σ_{L} and the dissipation parameter ɛ = σ_{0}/σ_{L} ≪ 1, and we rewrite Eq. (H.3) as
The modulus, real part, and the imaginary part of this complex surface pressure oscillation are then straightforwardly deduced as:
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
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 BruntVäisälä frequency is expressed as
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
equivalent to the κ parameter characterizing the dry adiabat. We also introduce the ratio q_{M}, such that 0 < q_{M} ≤ 1, defined as
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 BruntVäisälä frequency is thus rewritten as
With these definitions, assuming neutral stratification (N_{B} = 0) leads to
The other background fields are deduced from H by making use of the hydrostatic approximation and the ideal gas law,
with H_{s} = R_{s}T_{s}/g, and ρ_{s} = p_{s}/ (gH_{s}). Written as a function of the altitude, the normalized altitude coordinate x is expressed as
and the atmospheric vertical profiles, given in Eq. (13) for the dry adiabat, in the case of a moist adiabat become
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 thermodynamic equation thus becomes
By introducing the calculation variable G, we rewrite the above equation as
Finally, we normalize the equations and we switch to the frequency domain. This yields,
and the normalized equation for ,
with η = R_{p}/H and α = σ/σ_{0}. The structure factors in the tidal equations are expressed as
and
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
Finally, we use Eq. (B.2) to express as a function of and , and we end up with
The equation for 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
with the operators and coefficients
This system is reduced to a single second order equation of in the form given by Eq. (B.11) but with the modified operators A, B, and C given by
The wave equation is finally written in its classical form by substituting ℒ = −Λ_{n} and introducing the change of variable ,
with ϕ being expressed as
and the squared vertical wavenumber as
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
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
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
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_{2}O) and ozone (O_{3}).
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. (143144))
with the coefficients J^{m,σ} given by
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 . Following the authors (their Eq. 163), we set H = 7.6 km. The semidiurnal distributions of for water vapor and ozone can be approximated by the functions
with and . 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 of Eq.(J.3) is set to
where . Similarly, the vertical profile of ozone absorption can be approximated by a truncated sine function,
with , x_{1} = z_{1}/H, and x_{2} = z_{2}/H, the altitudes z_{1} and z_{2} being set to z_{1} = 18 km and z_{2} = 78 km. At a given location on the sphere (θ, φ), the total flux absorbed by the atmosphere for the spherical harmonic Y_{22} is defined as
The efficiency parameter α_{A} of our model is basically the ratio between δℱ_{22} and the corresponding component of the incident Solar flux, which we retrieve in Appendix E (Eq. E.13) as
In order to simplify calculations, we assume that the latitudinal profile of ozone can be approximated by the function sin^{2} θ, similar to the latitudinal profile of water vapor. Under this assumption, we rewrite Eqs.(J.3) and (J.4) as
and the corresponding heating rates as
where the imaginary number appearing in Eq.(J.2) has been ignored since it does not affectthe amplitude of the thermal forcing^{9}. For a given component, the total input power is obtained by integrating the heating rate over the air column
noting that ρ_{0}H = (p_{0}/g). This requires computing the integral of the vertical distribution of the heating rates,
We end up with and . Therefore, the total fluxes absorbed by water vapor and ozone are given by
with the components
Using σ = 1.45 × 10^{−4} s^{−1}, C_{p} = 1004.7 J kg^{−1} K^{−1}, p_{s} = 1013.25 hPa, g = 9.81 m s^{−2}, and F_{⋆} = 1366 W m^{−2}, we obtain δF_{22} = 828.8 W m^{−2}, and
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 semidiurnal atmospheric tide on Earth. Without elucidating again the copious details of the procedure, we obtain for the profiles in Vichare & Rajaram (2013) the following:
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 presentday semidiurnal surface pressure variations.
Appendix K Benchmarking our model against GCM simulations
The work of AuclairDesrotour 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 Venuslike planet.
Fig. K.1 Similar to Figure 2, the tidal spectra, in terms of the imaginary part of the surface pressure anomaly, are plotted as a function of the normalized tidal frequency ω. The blue curve is that computed numerically in AuclairDesrotour et al. (2019) using 3D GCM simulations, while the orange curve is recovered from our model by tuning the two free frequencies σ_{bl} and σ_{0} as discussed in Appendix K. 
Imposing on our analytical solution (Eq.22) the same planetary parameters defined in Table 1 of AuclairDesrotour et al. (2019), for a 10 bar atmosphere, and setting σ_{bl} = 10^{−5} s^{−1} and σ_{0} = 4 × 10^{−6} s^{−1}, 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 𝒯 ∝ σ^{−1}; while for our model, 𝒯 ∝ σ^{0.53} in the former regime, and 𝒯 ∝ σ^{−1} 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 σ_{0} ∝ T^{3}. Henceforth, the temperature increase required to place the resonance occurrence in the Precambrian would also attenuated the resonant amplitude of the torque by increasing σ_{0}. 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
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 Λ_{2} are indeed constant, the rest of the parameters would modify the expression of the amplification factor ŋ (Eq. L.1) into:
The resonant frequency in the adiabatic limit scales as . Moreover, the radiative cooling frequency is also a function of the surface temperature and it scales as (e.g., Eq. 10 of Showman & Guillot 2002). Including these scaling relationships in the latter expression of the amplification factor one finds that
The equilibrium surface temperature varies as T_{s} ∝ (α_{h}F_{⋆})^{1/4}, where α_{h} is the fraction of flux that is absorbed by the atmosphere, mainly in the infrared absorption. It follows that
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_{2}), 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 adiabatic 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; AuclairDesrotour 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 AuclairDesrotour et al. (2019) that yield a frequencydependence for the pressure surface anomaly in the form (their Eq. 12):
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} = Γ_{1} H (see also Lindzen & Blake 1972). Hence, introducing the definition of the equivalent depth
we obtain, for the fundamental mode, the resonant frequency and consequently the resonant rotational period in the form:
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 Γ_{1} factor multiplied by the temperature. To compare the two equations, we assign for the isothermal model a constant value of temperature corresponding to a massaverage of the temperature profile in the neutrally stratified case assuming a fixed total enthalpy of the air column. Namely, using Eq. 13, we obtain , thus we rewrite Eq. M.3 as:
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
 Abramowitz, M., Stegun, I. A., & Romer, R. H. 1988, Am. J. Phys., 56, 958 [Google Scholar]
 Andújar Márquez, J. M., Martínez Bohórquez, M. Á., & Gómez Melgar, S. 2016, Sensors, 16, 306 [Google Scholar]
 Arkhangelskaya, T., & Lukyashchenko, K. 2018, Biosyst. Eng., 168, 83 [Google Scholar]
 Arras, P., & Socrates, A. 2010, ApJ, 714, 1 [Google Scholar]
 AuclairDesrotour, P., Laskar, J., & Mathis, S. 2017a, A&A, 603, A107 [Google Scholar]
 AuclairDesrotour, P., Laskar, J., Mathis, S., & Correia, A. 2017b, A&A, 603, A108 [Google Scholar]
 AuclairDesrotour, P., & Leconte, J. 2018, A&A, 613, A45 [Google Scholar]
 AuclairDesrotour, P., Leconte, J., & Mergny, C. 2019, A&A, 624, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bartlett, B. C., & Stevenson, D. J. 2016, Geophys. Res. Lett., 43, 5716 [Google Scholar]
 Berger, B. W., & Grisogono, B. 1998, Bound. Layer Meteorol., 87, 363 [Google Scholar]
 Bernard, E. A. 1962, Arch. Meteorol. Geophys. Bioklimatol. A, 12, 502 [Google Scholar]
 Blake, R. E., Chang, S. J., & Lepland, A. 2010, Nature, 464, 1029 [Google Scholar]
 Catling, D. C., & Zahnle, K. J. 2020, Sci. Adv., 6, eaax1420 [Google Scholar]
 Chapman, S., & Lindzen, R. S. 1970, Atmospheric Tides: Thermal and Gravitational, 15 (Springer Science & Business Media) [Google Scholar]
 Charnay, B., Le Hir, G., Fluteau, F., Forget, F., & Catling, D. C. 2017, Earth Planet. Sci. Lett., 474, 97 [Google Scholar]
 Charnay, B., Wolf, E. T., Marty, B., & Forget, F. 2020, Space Sci. Rev., 216, 1 [Google Scholar]
 Correia, A., & Laskar, J. 2001, Nature, 411, 767 [Google Scholar]
 Correia, A. C., & Laskar, J. 2003a, J. Geophys. Res. Planets, 108 [Google Scholar]
 Correia, A. C., & Laskar, J. 2003b, Icarus, 163, 24 [Google Scholar]
 Correia, A. C. M., & Laskar, J. 2010, in Exoplanets (Tucson, AZ: University of Arizona Press), 239 [Google Scholar]
 Correia, A. C., Laskar, J., & de Surgy, O. N. 2003, Icarus, 163, 1 [Google Scholar]
 Covey, C., Dai, A., Lindzen, R. S., & Marsh, D. R. 2014, J. Atmos. Sci., 71, 1905 [Google Scholar]
 Cunha, D., Correia, A. C., & Laskar, J. 2015, Int. J. Astrobiol., 14, 233 [Google Scholar]
 Dai, A., & Wang, J. 1999, J. Atmos. Sci., 56, 3874 [Google Scholar]
 de Wit, M. J., & Furnes, H. 2016, Sci. Adv., 2, e1500368 [Google Scholar]
 Dobrovolskis, A. R., & Ingersoll, A. P. 1980, Icarus, 41, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Farhat, M., AuclairDesrotour, P., Boué, G., & Laskar, J. 2022, A&A, 665, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Feiccabrino, J., Graff, W., Lundberg, A., Sandström, N., & Gustafsson, D. 2015, Hydrology, 2, 266 [Google Scholar]
 Garratt, J. R. 1994, EarthSci. Rev., 37, 89 [Google Scholar]
 Gerkema, T., & Zimmerman, J. 2008, Lecture Notes, Royal NIOZ, Texel, 207 [Google Scholar]
 Gold, T., & Soter, S. 1969, Icarus, 11, 356 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Soter, S. 1966, Icarus, 5, 375 [Google Scholar]
 Gough, D. 1981, in Physics of Solar Variations: Proceedings of the 14th ESLAB Symposium held in Scheveningen, The Netherlands, 16–19 September, 1980 (Springer), 21 [Google Scholar]
 Gu, P.G., Peng, D.K., & Yen, C.C. 2019, ApJ, 887, 228 [Google Scholar]
 Hagan, M., & Forbes, J. 2002, J. Geophys. Res. Atmos., 107, ACL [Google Scholar]
 Haurwitz, B., & Cowley, A. D. 1973, Pure Appl. Geophys., 102, 193 [Google Scholar]
 Hilsenrath, J. 1955, Tables of thermal properties of gases: comprising tables of thermodynamic and transport properties of air, argon, carbon dioxide, carbon monoxide, hydrogen, nitrogen, oxygen, and steam, 564 (US Department of Commerce, National Bureau of Standards) [Google Scholar]
 Holmberg, E. 1952, Geophys. J. Int., 6, 325 [Google Scholar]
 Holtslag, A., & Boville, B. 1993, J. Climate, 6, 1825 [Google Scholar]
 Hough, S. S. 1898, Philos. Trans. Roy. Soc. Lond. Ser. A, 191, 139 [Google Scholar]
 Hourdin, F., Musat, I., Bony, S., et al. 2006, Climate Dyn., 27, 787 [Google Scholar]
 Hut, P. 1981, A&A, 99, 126 [NASA ADS] [Google Scholar]
 Ingersoll, A. P., & Dobrovolskis, A. R. 1978, Nature, 275, 37 [CrossRef] [Google Scholar]
 Kelvin, L. O. K. A. W. T. 1882, Proc. Roy. Soc. Edinb., 11, 396 [Google Scholar]
 Klatt, J. M., Chennu, A., Arbic, B. K., Biddanda, B., & Dick, G. J. 2021, Nat. Geosci., 14, 564 [Google Scholar]
 Knauth, L. P. 2005, in Geobiology: Objectives, Concepts, Perspectives (Elsevier), 53 [Google Scholar]
 KrissansenTotton, J., Arney, G. N., & Catling, D. C. 2018, PNAS, 115, 4105 [Google Scholar]
 Lacis, A. A., & Oinas, V. 1991, J. Geophys. Res. Atmos., 96, 9027 [Google Scholar]
 Lambeck, K. 1980, The Earth’s Variable Rotation: Geophysical Causes and Consequences (Cambridge University Press) [Google Scholar]
 Lantink, M. L., Davies, J. H., Ovtcharova, M., & Hilgen, F. J. 2022, PNAS, 119, e2117146119 [Google Scholar]
 Laplace, P. 1798, Traite de Mecanique Celeste, five volumes (Paris: Chez JBM Duprat) [Google Scholar]
 Laskar, J., & Correia, A. C. 2004, Extrasolar Planets: Today and Tomorrow, 321, 401 [Google Scholar]
 Laskar, J., Farhat, M., Lantink, M. L., et al. 2023, arXiv eprints [arXiv:2309.11479] [Google Scholar]
 Leconte, J., Wu, H., Menou, K., & Murray, N. 2015, Science, 347, 632 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, U. 2020, MNRAS, 494, 3141 [Google Scholar]
 Lee, U., & Saio, H. 1997, ApJ, 491, 839 [Google Scholar]
 Lindzen, R. S. 1978, Monthly Weather Rev., 106, 526 [Google Scholar]
 Lindzen, R. S., & Blake, D. 1972, J. Geophys. Res., 77, 2166 [Google Scholar]
 Lindzen, R. S., & McKenzie, D. J. 1967, Pure Appl. Geophys., 66, 90 [Google Scholar]
 Madsen, O. S. 1977, J. Phys. Oceanogr., 7, 248 [Google Scholar]
 Mignard, F. 1980, Moon Planets, 23, 185 [Google Scholar]
 Mitchell, R. N., & Kirscher, U. 2023, Nat. Geosci., 1, 567 [Google Scholar]
 Nieuwstadt, F. 1983, Bound. Layer Meteorol., 26, 377 [Google Scholar]
 Ogilvie, G. I. 2014, ARA&A, 52, 171 [Google Scholar]
 Pannella, G. 1972a, Astrophys. Space Sci., 16, 212 [Google Scholar]
 Pannella, G. 1972b, Int. Geol. Congr. Rep. Sess., 24th, 50 [Google Scholar]
 Pierrehumbert, R. T. 2010, Principles of Planetary Climate (Cambridge University Press) [Google Scholar]
 Robert, F., & Chaussidon, M. 2006, Nature, 443, 969 [Google Scholar]
 Sakazaki, T., & Hamilton, K. 2020, J. Atmos. Sci., 77, 2519 [Google Scholar]
 Sakazaki, T., Hamilton, K., Zhang, C., & Wang, Y. 2017, Geophys. Res. Lett., 44, 1998 [Google Scholar]
 Schindelegger, M., & Ray, R. D. 2014, Monthly Weather Rev., 142, 4872 [Google Scholar]
 Scrutton, C. 1978, in Tidal Friction and the Earth’s Rotation (Springer), 154 [Google Scholar]
 Showman, A. P., & Guillot, T. 2002, A&A, 385, 166 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Siebert, M. 1961, in Advances in Geophysics, 7 (Elsevier), 105 [Google Scholar]
 Sleep, N. H., & Zahnle, K. 2001, J. Geophys. Res.: Planets, 106, 1373 [Google Scholar]
 Tiesinga, E., Mohr, P. J., Newell, D. B., & Taylor, B. N. 2021, J. Phys. Chem. Ref. Data, 50, 033105 [Google Scholar]
 Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1979, Nonradial oscillations of stars (Tokyo: University of Tokyo Press) [Google Scholar]
 Vallis, G. K. 2017, Atmospheric and Oceanic Fluid Dynamics (Cambridge University Press) [Google Scholar]
 Van Wijk, W. R., & De Vries, D. 1963, Phys. Plant Environ., 1, 103 [Google Scholar]
 Vichare, G., & Rajaram, R. 2013, J. Earth Syst. Sci., 122, 1207 [Google Scholar]
 Walker, J. C., & Zahnle, K. J. 1986, Nature, 320, 600 [Google Scholar]
 Wang, H., Boyd, J. P., & Akmaev, R. A. 2016, Geosci. Model Dev., 9, 1477 [Google Scholar]
 Wilkes, M. V. 1949, Oscillations of the Earth’s Atmosphere (Cambridge University Press) [Google Scholar]
 Williams, G. E. 1990, J. Phys. Earth, 38, 475 [Google Scholar]
 Williams, G. E. 2000, Rev. Geophys., 38, 37 [Google Scholar]
 Wu, H., Murray, N., Menou, K., Lee, C., & Leconte, J. 2023, Sci. Adv., 9, eadd2499 [Google Scholar]
 Zahnle, K., & Walker, J. C. 1987, Precamb. Res., 37, 95 [Google Scholar]
Noting that surface friction is another dissipative mechanism as discussed by Lindzen & Blake (1972).
We note that this form corresponds to the quadrupolar component of the torque about the spin axis, and it is only valid assuming a thin atmospheric layer under the hydrostatic approximation. In the case of a thick atmosphere, one should integrate the mass redistribution over the radial direction.
This value does not correspond to the maximum amplitude of the surface pressure oscillation at the equator, but to the coefficient of its quadrupolar spherical harmonic component, which via Appendix E and Appendix G is normalized by a factor of .
We note that in Leconte et al. (2015), the tidal frequency under study is that of the diurnal component, thus we multiply their ω_{0} value by 2; i.e., σ_{0} = 2ω_{0}.
All Tables
Dimensionless control parameters determining the regime of the atmospheric tidal response and defined throughout the text.
All Figures
Fig. 1 Modeled histories of the rotational motion of the Earth. The Earth’s LOD evolution is plotted in time over geological timescales for three models: (i) the model of Farhat et al. (2022), where the evolution is driven solely by oceanic and solid tidal dissipation; and (ii) the model of Zahnle & Walker (1987), where the Lamb resonance is encountered for LOD~21 h, forcing a rotational equilibrium on the Earth; (iii) the model of Bartlett & Stevenson (2016), which also adopts the equilibrium scenario, but further studies the effect of thermal noise, and the required temperature variation to escape the equilibrium. Three curves of the latter model correspond to different parameterizations of the gravitational tide. Plotted on top of the modeled histories are geological proxies of the LOD evolution which can be retrieved from http://astrogeo.eu. 

In the text 
Fig. 2 Spectrum of semidiurnal atmospheric thermal tides. Plotted is the imaginary part of the normalized pressure anomaly (; Eq. (22)) as a function of the normalized forcing frequency ω = (Ω − n_{⋆})/n_{⋆} = σ/2n_{⋆}, where n_{⋆} is the mean motion of the stellar perturber. The planetarystellar parameters are those of the fiducial planetary system defined in Sect. 3.1. 

In the text 
Fig. 3 (A)symmetry of the Lamb resonance. Similar to Fig. 2, the imaginary part of the normalized pressure anomaly (Eq. (22)), which is associated with the semidiurnal tide, is plotted as a function of the normalized forcing frequency ω = (Ω − n_{⋆})/n_{⋆} = σ/2n_{⋆}, for the same planetarystellar parameters. We focus here on the high frequency regime around the Lamb resonance. Different panels correspond to different values of σ_{bl}, or different thermal inertias in the ground and the atmosphere. Allowing for thermal inertia results in a delayed ground response, of which the signature is clear in inducing an asymmetry in the spectral behavior around the resonance. 

In the text 
Fig. 4 Thermally induced tidal bulge revealed. Shown are polar snapshots of the radial and longitudinal variations of the tidal pressure anomaly δp(x) in the equatorial plane. The snapshots are shown from a top view, and the troposphere is puffed in size by virtue of the used massbased vertical coordinate ϛ = p/p_{s}. The longitudinal axes are shown in increments of 30° with 0° at the substellar point, while the radial axes are in increments of 0.25. The profile of the pressure perturbation is also normalized by the exponentially decaying pressure background profile. Snapshots are taken at different spectral positions that cover the passage through the Lamb resonance, which specifically occurs at ω = 262.6. In the top row, the response describes the limit of a planet with a synchronous atmospheric absorption, mimicking the Earth’s direct absorption by ozone and water vapor, and it shows the continuous movement of the bulge, function of ω, from lagging to leading the substellar point. In contrast, in the bottom row, and for the prescribed value of σ_{bl}, the delayed response of the ground forces the bulge to always lag the substellar point, thus acting to decelerate the planetary rotation. 

In the text 
Fig. 5 Parametric study of the tidal response. Plotted is a contoured surface of the amplitude of the imaginary part of the positive semidiurnal pressure anomaly at the Lamb resonance, over a grid of values of our free parameters σ_{0} and α_{A}. The solid black isoline marks the level curve of Im{δp_{s}} = 880 Pa, and defines from below a region in (α_{A}, σ_{0})space where the thermotidal response is sufficient to cancel the gravitational counterpart in the Precambrian. Analogously, the dashed isoline defines the threshold needed in the late Paleozoic/early Mesozoic, 350−250 Ma. The horizontal shaded area corresponds to typical values of the radiative cooling rate as described in the main text. The other shaded area defines the region of parameter space that yields the presently observed semidiurnal tidal bulge. The gray area on the left covers the parametric region where the resonance features a lower pressure amplitude than the present. 

In the text 
Fig. 6 Dependence of the resonant rotational period on the mean surface temperature. By virtue of Eq. (24), the LOD at which the Lamb resonance occurs scales as the inverse square root of the mean surface temperature (solid curve). In Appendix I, we show that the same dependence holds if one considers a moist adiabatic profile. The shaded area of temperature variations highlights 95% confidence intervals for the past temperature evolution according to the carbon cycle model of KrissansenTotton et al. (2018). The identified geological eras correspond to the LOD evolution model of Farhat et al. (2022). The overlap between the modeled temperature evolution and the black curve places the resonance occurrence in the late Paleozoic/early Mesozoic. 

In the text 
Fig. D.1 Schematic of the bichromatic radiative and diffusive flux exchange mechanisms at the surface interface between the lower troposphere and the ground. The power balance between the fluxes in given by Eq. (D.1). 

In the text 
Fig. H.1 Sample of the tidal response, plotted in terms of the modulus, real part, and imaginary part of the surface pressure anomaly, all given by Eqs. (H.7). The function f(x) on the yaxis indicates the functional form of the pressure variation, that is, the pressure function divided by the constant factor K. The xaxis is centered at the position of the Lamb resonance. Notable is the difference in the resonant amplification of the response between the real part and the imaginary part. 

In the text 
Fig. K.1 Similar to Figure 2, the tidal spectra, in terms of the imaginary part of the surface pressure anomaly, are plotted as a function of the normalized tidal frequency ω. The blue curve is that computed numerically in AuclairDesrotour et al. (2019) using 3D GCM simulations, while the orange curve is recovered from our model by tuning the two free frequencies σ_{bl} and σ_{0} as discussed in Appendix K. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.