Thermal tides in neutrally stratified atmospheres: Revisiting the Earth’s Precambrian rotational equilibrium

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 closed-form solution for the frequency-dependent 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.


Introduction
For present day Earth, the semi-diurnal atmospheric tide, driven by the thermal forcing of the Sun and generated via atmospheric pressure waves, describes the movement of atmospheric mass away from the substellar point.Consequently, mass culminates forming bulges on the nightside and the dayside, generating a torque that accelerates the Earth's rotation.As such, this thermally generated torque counteracts the luni-solar gravitational torque associated with the Earth's solid and oceanic tides.The latter components typically drive the closed system of the tidal players towards 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 me-chanical 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 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 and 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 and Soter, 1966;Gold and Soter, 1969;Ingersoll and Dobrovolskis, 1978;Dobrovolskis and Ingersoll, 1980;Correia andLaskar, 2001, 2003b;Correia et al., 2003), hot Jupiters (e.g., Arras and Socrates, 2010;Auclair-Desrotour and Leconte, 2018;Gu et al., 2019;Lee, 2020), and near-synchronous and Earth-like rocky exoplanets (Cunha et al., 2015;Leconte et al., 2015;Auclair-Desrotour et al., 2017a, 2019).Namely, for planetary systems within the so-called habitable zone, the gravitational tidal torque diminishes in the regime near spin-orbit synchronization and becomes comparable in magnitude to the thermotidal torque.Consequently, the latter may actually prevent the planet from precisely reaching its destined synchronous state (Laskar and Correia, 2004;Correia and 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 present Earth (Chapman and Lindzen, 1970), Zahnle and 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 Plotted is the Earth's LOD evolution 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; ii) the model of Zahnle and Walker (1987), where the Lamb resonance is encountered for LOD∼21 hr, forcing a rotational equilibrium on the Earth; iii) the model of Bartlett and 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 modelled histories are geological proxies of the LOD evolution that can be retrieved from http://astrogeo.eu.
competing torques.The first is the occurrence of a resonance in atmospheric Lamb waves (e.g., Lindzen and 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 and Walker (1987), this resonance occurred when the length of day (LOD) was around 21 hrs, exciting the thermotidal torque to large amplitudes.Secondly, the gravitational tidal torque must have been largely attenuated in the Precambrian.Recently, Bartlett and Stevenson (2016) revisited the equilibrium scenario and investigated the ef-fect 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 hr to its present value.Evidently, the occurrence of such a scenario has very significant implications on paleoclimatic studies, with the growing evidence on links between the evolving LOD and the evolution of Precambrian benthic life (e.g., Klatt et al., 2021).
We are fresh out of a study on the tidal evolution of the Earth-Moon system (Farhat et al., 2022), where we focused on modelling tidal dissipation in the Earth's paleo-oceans and solid interior.There we learned that the tidal response of the oceans, characterized by intermittent resonant excitations, is sufficient to explain the present rate of lunar recession and the estimated lunar age, and 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 (Figure 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 and Dobrovolskis, 1978;Dobrovolskis and Ingersoll, 1980;Auclair-Desrotour et al., 2017a,b), to parameterized models that capture essential spectral features (e.g., Correia andLaskar, 2001, 2003a), to fully numerical efforts that relied on the advancing sophistication of general circulation models (GCM; e.g., Leconte et al., 2015;Auclair-Desrotour et al., 2019).The latter work presents, to our knowledge, the first and, to-date1 , the only study to have numerically computed the planetary thermotidal torque in the high frequency regime, i.e. around the Lamb resonance (Lindzen and Blake, 1972).Of interest to us here are two perplexing results that Auclair-Desrotour et al. (2019) established: first, for planets near synchronization, the simplified Maxwellian models often used to characterize the thermotidal torque did not match the GCM simulated response; second, the torque at the Lamb resonance featured only a decelerating effect on the planet.Namely, it acts in the same direction of gravitational tides, and thus the effect required for the rotational stabilization disappeared.
More recently, while this work was under review, two studies on the Precambrian LOD stabilization were published.Mitchell and 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 hr 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 hr between 0.6 and 2.2 Ga, with a sustained very high mean surface temperature (40 − 55 • C).Although using different approaches, the two studies have thus arrived at similar conclusions.A closer look at the subset of geological data that favored this outcome, however, indicates that both studies heavily rely on three stromatolitic records from the Paleoproterozoic that were originally studied by Pannella (1972a,b).These geological data have been, ever since, identified as unsuitable for precise quantitative interpretation (see e.g., Scrutton, 1978;Lambeck, 1980;Williams, 2000).To this end, we provide a more detailed analysis of the geological proxies of the LOD, and of the model presented by Wu et al. (2023) in a parallel paper dedicated to the matter (Laskar et al., 2023).
With a view to greater physical realism, we aim here to study, analytically, the frequency spectrum of the thermotidal torque, from first principles, interpolating between the low and high frequency regimes.Our motivation is two-fold: first, to provide a novel physical model for the planetary thermotidal torque that better matches the GCM-computed response, and that can be used in planetary dynamical evolution studies; second, to apply this model to the Earth and attempt quantifying the amplitude of the torque at the Lamb resonance and explore the intriguing rotational equilibrium scenario.

Ab initio atmospheric dynamics
For an atmosphere enveloping a spherically symmetric planet, we define a reference frame corotating with the planet.In this frame, an atmospheric parcel is traced by its position vector r in spherical coordinates (r, θ, φ), such that θ is the colatitude, φ is the longitude, and the radial distance |r| = R p +z, where R p is the planet's radius and z is the parcel's atmospheric altitude.The atmosphere is characterized by the scalar fields of pressure p, temperature T , density ρ, and the three-dimensional vectorial velocity field V.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 = 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 Section 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 2 (Lindzen and 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 noting that surface friction is another dissipative mechanism as discussed by Lindzen and Blake (1972).

The vertical structure of tidal dynamics
We are interested in providing a closed form solution for the frequency3 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 and Lindzen, 1970), atmospheric heat transfer augmented with linear radiative transfer à la Lindzen and McKenzie (1967), and the ideal gas law, all formulated in a dimensionless form.
Aided by the so-called traditional approximation (e.g., Unno et al., 1989, see also Appendix A), the analytical treatment of the said system is feasible as it decomposes into two parts describing, separately, the horizontal and vertical structures of tidal flows.The former part is completely described by the eigenvalue-eigenfunction problem defined as Laplace's tidal equation (Laplace, 1798;Lee and Saio, 1997): where the set of Hough functions {Θ m,ν n } serves as the solution (Hough, 1898), {Λ m,ν n } is the associated set of eigenvalues, L m,ν is a horizontal operator defined in Eq.(A.23) of Appendix A, while ν = 2Ω/σ, where Ω is the rotational velocity of the planet and σ is the tidal forcing frequency.In the tidal system under study, the variables and functions (δp, δρ, δT, V , J, Θ, Λ) are written in the Fourier domain using the longitudinal order m and frequency σ (Eq.A.20), and expanded in horizontal Hough modes with index n (Eq.A.21).We denote hereafter their coefficients f m,ν n 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 {Λ m,ν n }.To construct these sets of eigenfunctions-eigenvalues we use the spectral method laid out by Wang et al. (2016).
The vertical structure on the other hand requires a more elaborate manipulation of the governing system, a procedure that we detail in Appendix B. The outcome is a wave-like equation that describes vertical thermotidal dynamics and reads as: Here, as is the common practice (e.g., Siebert, 1961;Chapman and Lindzen, 1970), we use the reduced altitude x = z 0 dz/H(z) as the vertical coordinate, where the pressure scale height H(z) = R s T 0 (z)/g; R 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 kn (x) is defined via By virtue of the non-dimensionalization of the governing system of equations, dimensionless control parameters appear in the wavenumber definition.Namely, α = σ/σ 0 , β = σ 2 w /σ 2 , and γ = N 2 B /N 2 B;iso , where we have introduced the characteristic frequency σ w = √ gH/R p , a typical frequency of Lamb waves.Of significance to the computations that follow is the Brunt-Väisälä frequency, N B , a measure of the vertical density stratification of the atmosphere against the strength of convection (e.g, Vallis, 2017).Under hydrostatic equilibrium, N B is defined as: We define the Brunt-Väisälä frequency in the limit of an isothermal atmospheric profile, N B;iso = κg/H.As such, the parameter γ measures the local atmospheric stability against convection with respect to the stability of an equivalent isothermal atmosphere.These key control parameters, among others, qualitatively determine the regime of the tidal response.We summarize these parameters in Table 1.The right hand side of the wave equation (2) combines the function Φ defined as (5) 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 Jn = κ/(gσH)J n , while the gravitational analogue Ũn = U n /v 2 0 , where the reference velocity v 0 = σR p (see Eq.(A.11) for the rest of the dimensionless quantities).As we intend to solve the wave equation in the following sections, what is left for us to quantify the mass redistribution and compute the resulting tidal torque is to retrieve the vertical profile of pressure given the solution of the wave equation, Ψ n (x).In Appendix C, we derive the vertical profiles of all the tidal variables, and specifically for the dimensionless pressure anomaly we obtain: where δp n (x) = δp n (x)/p 0 , and the calculation variable Gn (x) = Ψ n (x)Φ(x).

The thermal forcing profile
To solve the non-homogeneous wave equation (2), it is necessary to define a vertical profile for the tidal heating power per unit mass Jn (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 and 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 σ bl = τ −1 bl , 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 α = σ/σ 0 , we define the dimensionless parameter for the boundary layer ζ = |σ|τ bl = |σ|/σ bl .
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, i.e. 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 and Lindzen, 1970;Schindelegger and 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 and 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 non-synchronous with the incident flux due to the delayed thermal response of the ground, which here is a function of τ bl .This limit better describes dry Venus-like planets, as is the fiducial setting studied using GCMs in Leconte et al. (2015) and Auclair-Desrotour et al. (2019).
Finally, as we are interested in the semi-diurnal tidal response, we decompose the thermal forcing in Appendix E to obtain the amplitude of the quadrupolar component as δF inc = δF 22 = ( √ 30π/16)F ⋆ , where F ⋆ = L ⋆ /4πa 2 p , L ⋆ being the stellar luminosity, and a p the star-planet distance.

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 and McKenzie, 1967;Auclair-Desrotour et al., 2019), in which case the scale height H becomes independent of the vertical coordinate, and by virtue of Eq. (4), N 2 B = κg/H = const.Motivated by the Earth's atmosphere, where the massive troposphere (∼80% of atmospheric mass) controls the tidal mass redistribution, we derive next an analytical solution in a different, and perhaps more realistic limit.Namely, the limit corresponding to the case of a neutrally stratified atmosphere, where N 2 B = 0.In fact, N 2 B can be expressed in terms of the potential temperature Θ 0 (e.g., Section 2.10 of Vallis, 2017): whereby the stability of the atmosphere is controlled by the slope of Θ 0 .That said, atmospheric temperature measurements (e.g., Figures 2.1-2.3 of Pierrehumbert, 2010) clearly depict that the troposphere is characterized by a negative temperature gradient, and a very weak potential temperature gradient, which is closer to an idealised adiabatic profile than it is to an idealised isothermal profile.Moreover, the heating in the troposphere generates strong convection and efficient turbulent stirring, thus enhancing energy transfer and driving the layer towards an adiabatic temperature profile.As such, the temperature profile being adiabatic would prohibit the propagation of buoyancy-restored gravity waves, which compose the baroclinic component of the atmospheric tidal response (e.g., Gerkema and Zimmerman, 2008).This leaves the atmosphere with the barotropic component of the tidal flow, a feature 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 It also follows that the background profiles of the scalar variables read as (Auclair-Desrotour et al., 2017a): We thus obtain for the heating profile (using Eqs. 8, 9, and 10) (14) As such, the wave equation ( 2) is rewritten as where the complex functions A n and B are defined as: The wave equation ( 15) admits the general solution We consider the following two boundary conditions: • First, the energy of tidal flows, W, should be bounded as x → ∞.In Appendix F, we derive the expression of the tidal energy following Wilkes (1949), and it scales as W ∝ |Ψ| 2 |Φ| 2 .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 Ṽr;n (x = 0) = 0. We derive the expression of the profile of the vertical velocity in Appendix C, and by virtue of Eq.(C.6), this condition allows us to write c 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 Gn (x)), but we are specifically interested in retrieving a closed form solution of the quadrupolar tidal torque.The latter takes the general form5 (Appendix G): Here M ⋆ and M p designate the stellar and planetary masses respectively, and ℑ refers to the imaginary part of a complex number, the latter in this case being the quadrupolar pressure anomaly at the surface δp s = δp 2,ν 2 (x = 0).We further note that while this torque is computed for the atmosphere, it does act on the whole planet since the atmosphere is a thin layer that features no differential rotation with respect to the rest of the planet.
Taking the solution Ψ(0) of Eq. ( 18) (with c 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 X and Y 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. Before embarking on any results, we pause here for a few remarks on the provided closed form solution of the torque.
• The parameter α A , defined earlier (Eq.9) as the fraction of radiation actually absorbed by the atmosphere, can evidently be correlated with the typical transmission function of the atmosphere and therefore its optical depth.Presuming that 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 Section 4.1 and Appendix I.
• The quadrupolar component of the equilibrium stellar flux, entering through a fraction of F ⋆ (Appendix E), is directly proportional to the stellar luminosity L ⋆ .Standard models suggest that the Sun's luminosity was around 80% of its present value ∼3 Ga (Gough, 1981).Such luminosity evolution of Sun-like stars can be directly accommodated in the model if one were to study the evolution of the tidal torque with time.
• As we mentioned earlier, upon separating the horizontal and vertical structure of tidal dynamics, the only remaining coupling factor between the two structures is the eigenvalue of horizontal flows, Λ n , in our case reducing to the dominant fundamental mode Λ 2 .Noting that we have dropped the superscripts, we remind the reader that for the semidiurnal (m = 2) response, Λ 2 = Λ 2,ν 2 , thus Λ is frequency-dependent in the general case.The Earth however, over its lifetime, lives in the asymptotic regime of ν ≈ 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 closed form solution ( 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 α, this term would have encountered a singularity at the spectral position of the resonance, i.e. 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 Section 4.2.
In Fig. 2, we plot the spectrum of the tidal response for a fiducial system in terms of the normalized surface pressure anomaly over a wide range of tidal frequencies covering the low and high frequency regimes.The system describes a Venus-like dry planet (M p = 0.815M ⊕ , R p = 0.95R ⊕ , 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 solar-like star (M ⋆ = 1M ⊙ , L ⋆ = 1L ⊙ ).We further ignore here the thermal inertia in the ground and the atmosphere by taking σ bl → ∞, or ζ → 0, thus assuming an synchronous response of the ground with the thermal excitation.Key tidal response features are recovered in this spectrum: First, we obtain a tidal peak near synchronization (ω = 0) that generates a positive torque for σ > 0 and a negative torque for σ < 0, driving the planet in both cases away from its destined spin-orbit synchronization due to the effect of solid tides (e.g., Gold and Soter, 1969;Correia and Laskar, 2001;Leconte et al., 2015).The peak has often been modelled by a Maxwellian functional form, though this form does not always capture GCM-generated spectra when varying the planetary setup (e.g., Auclair-Desrotour et al., 2019).Second, we recover the Lamb resonance in the high frequency regime.The resonance is characterized here by two symmetric peaks of opposite signs.Thus upon passage through the resonance, the thermotidal torque shifts from being a rotational brake to being a rotational pump.In this work, we are more interested in the high frequency regime, thus we delegate further discussion and analysis on the low frequency tidal response to a forthcoming work, and we focus next on the Lamb resonance.

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 Figure 2 -observed in GCM simulations of an atmosphere forced by a longwave flux (Auclair-Desrotour et al., 2019).
In Figure 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 Figure 3, the thermotidal torque in this regime is also negative, thus it complements the effect of gravitational tides.When the resonance is encountered, the thermotidal torque shifts its sign to counteract the effect of gravitational tides, with an amplified effect in the vicinity of the resonance.However, with the introduction of thermal inertia into the linear theory of tides, the L Ω -pumping part of the atmospheric torque is attenuated, and for some values of σ bl , it completely disappears.This modification of the analytical theory allows us to explain the asymmetry of the Lamb resonance depicted in the 3D GCM simulations of Auclair-Desrotour et al. (2019).In Appendix J, 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 Figure 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 verti-cal profile of the pressure anomaly from Eq. ( 7), and augment it with the latitudinal and longitudinal dependencies from Eqs. (A.20-A.21).As the massive troposphere dominates the tidal mass redistribution, we use the mass-based vertical coordinate ς = p/p s (i.e.x = − ln ς, and ς ranges between 1 at the surface and 0 in the uppermost layer).
In Figure 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, i.e. 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 towards lagging the substellar point thus decelerating the planet's rotation.This behavior is consistent with the established response spectrum in the first panel of Figure 3, and is relevant to the Earth's case, assuming that thermotidal heat-  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.
ing 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 Figure 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 and Lindzen, 1970) on the basis of the smallamplitude non-migrating tidal components they produce (e.g., Schindelegger and Ray, 2014).

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

Was the resonance resonant enough? A parametric study
Constraints on the amplitude of the gravitational tide during the Precambrian are model-dependent.The study in Zahnle and Walker (1987), and later in Bartlett and 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 and Zahnle, 1986).The record is fraught with multiple interpretations featuring different inferred values for the LOD (Williams, 1990(Williams, , 2000)), altogether different from a recent cyclostratigraphic inference that roughly has the same age (Lantink et al., 2022, see Figure 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 Figure 3; top row of Figure 4).In Figure 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 upwards 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 black isoline 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 N m) is roughly a quarter of the present gravitational torque (∼4.51 × 10 16 N m) (e.g., Zahnle and Walker, 1987;Farhat et al., 2022), thus requiring, via Eq.( 21), ℑ{δp s } on the order of 880 Pa 6 .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 Mesozoic (ℑ{δp s } = 2275 Pa).The solid gray region on the left side of the parameter space is bounded by the isoline corresponding to the present value ℑ{δp s } = 224 Pa (Schindelegger and 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 layer-to-layer radiative transfer and assume a gray body atmospheric radiation directly into space.However, radiative transfer can be consistently accommodated in numerical GCMs using the method of correlated k-distributions (e.g., Lacis and Oinas, 1991) as performed in Leconte et al. (2015) and in Auclair-Desrotour et al. (2019), both studies using the LMD GCM (Hourdin et al., 2006).In fact, Leconte et al. (2015) fitted their numerically obtained atmospheric torques to effective values of σ 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.In contrast, Lindzen and Chapman (1968) and later Lindzen and Blake (1972) estimated the timescale to be on the order of 1 day.We presume that these estimates should encompass the possible effective values for the Earth's atmosphere, and we highlight with the horizontal shaded area the range of these values 7 .
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 and Cowley (1973); Dai and Wang (1999); Covey et al. (2014) and Schindelegger and Ray (2014), which constrain the amplitude of the semi-diurnal surface pressure oscillation to within 107 − 150 Pa, occurring around 0945 LT.The narrow shaded area defines the region of parameter space that can explain these observables using the present semi-diurnal frequency, placing the opacity parameter in the region α A ∼14%.In Appendix I, 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 and Lindzen, 1970), or using GCM simulations (Vichare and 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.
Evidently, the overlap of the parametric constraints lives outside the region where the thermotidal response is sufficient for the rotational equilibrium condition.The present thermotidal torque (2.89 × 10 15 N m) needs to be amplified by a factor of 3.9 to reach the absolute minimum of the opposing gravitational torque8 in the Precambrian, and by a factor of 12.3 to reach the Mesozoic value.Our parametric exploration precludes these levels of amplification.It is important to also note that larger amplification factors would be required if one were focused on the modulus of the pressure oscillation, rather than its imaginary part.This derives from Figure H.8 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 con-straints derive from present measurements, and the likelihood of the scenario still hinges on possible atmospheric variations as we go backwards in time.Nonetheless, the radiative cooling timescale exhibits a strong dependence on the equilibrium temperature of the atmosphere (σ 0 ∝ T 3 0 ; Auclair-Desrotour et al., 2017a, Eq. 17).As such, a warmer planet in the past would yield a shorter cooling timescale, and consequently, more efficient damping of the resonant amplitude (see Appendix K).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 α A = 24% can indeed place the response beyond the Precambrian threshold for some values of σ 0 .An increase to four times the present value of α A is required to cross the threshold in the Mesozoic.These increases, however, can be precluded, based on the fact that the Archean lacked a stratospheric ozone layer (e.g., Catling and Zahnle, 2020).In contrast, an increase in atmospheric opacity in the infrared, which accompanies the abundance of Precambrian greenhouse gases, delivers the opposite effect by attenuating the resonant tidal response, as we elaborate in Appendix K. Furthermore, the latter increase would also trigger the contribution of asynchronous tidal heating, which further attenuates the amplitude of the positive peak as we show in Section 3.2.Thus, with these analyses, it is unlikely that the resonance could have amplified the thermotidal response beyond the required threshold.This conclusion can be further regarded as conservative, since the employed linear model tends to overestimate the resonant amplification of the tidal response.This derives from the fact that, in the quasi-adiabatic regime, the model ignores the associated non-linearities of dissipative mechanisms.The remaining question is therefore: when did the Lamb resonance actually occur?

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  24), the LOD at which the Lamb resonance occurs scales as the inverse square root of the mean surface temperature.The gray shaded area highlights 95% confidence intervals for the past temperature evolution according to the carbon cycle model of Krissansen-Totton et al. (2018).The identified geological eras correspond to the LOD evolution model of Farhat et al. (2022).The overlap between the modelled temperature evolution and the black curve places the resonance occurrence in the early Mesozoic.
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 Figure 6, we plot the dependence of the spectral position of the resonance, in terms of LOD, on T s .The apparent single curve 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: The resonant rotational period thus 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 and 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 and Furnes, 2016), geological carbon cycle models (e.g., Sleep and Zahnle, 2001;Krissansen-Totton et al., 2018), numerical results of 3D GCMs (e.g., Charnay et al., 2017), and the fact that solar luminosity was 10-25% lower during the Precambrian (e.g., Charnay et al., 2020), altogether predicting a temperate climate and moderate temperatures throughout the Earth's history.
We highlight with the gray shading on top of the curve modelled mean surface temperature variations adopted from Krissansen-Totton et al. (2018).As the latter temperature evolution is established in the time domain, we use the LOD evolution in Farhat et al. (2022) to map from time-dependence to LOD-dependence, and we further identify the corresponding geological eras of the LOD evolution with the color shadings.Given the present day equilibrium surface temperature, the resonance occurs at LOD = 22.8 hr.This value is in agreement with the 11.38±0.16hr semi-diurnal period obtained by analyzing the spectrum of normal modes using pressure data on global scales (see Table 1 of Sakazaki and Hamilton, 2020, first symmetric gravity mode of wavenumber k = −2).In Appendix L, we compute the resonant rotational period assuming an isothermal profile of the atmosphere, and we show that it is roughly one hour less than that in the neutrally stratified limit, placing it closer to 21.3 hr estimate of Zahnle and Walker (1987) and Bartlett and 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 Figures 2 and H.8), the difference between the two being dependent on the radiative cooling frequency.
Taking the LOD evolution model of Farhat et al. (2022) at face value, the temperature variations predicted in Krissansen-Totton et al. (2018) locate the resonance encounter in the Triassic, 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 hr, 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, σ 0 ∝ T 3 0 ; see also Appendix K), reducing the radiative cooling timescale and prompting more efficient damping of the tidal amplitude at the resonance.Moreover, such an increase in temperature would most probably 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 Section 3.2, whereby the accelerative peak of the torque is further attenuated.

Summary and Outlook
We were drawn to the problem of atmospheric thermal tides by the hypothesized scenario of a constant length of day on Earth during the Precambrian.Our motivation in investigating the scenario lies in its significant implications on paleoclimatic evolution, and the evident mismatch between LOD geological proxies and the predicted LOD evolution if this rotational equilibrium is surmised.The scenario hinges on the occurrence of a Lamb resonance in the atmosphere whereby an amplified thermotidal torque would cancel the opposing torque generated by solid and oceanic gravitational tides.Naturally, the atmospheric tidal torque is that of two flavors: it can either pump or deplete the rotational angular momentum budget of the planet, depending on the orientation of the generated tidal bulge.
With this rotational equilibrium scenario in mind, we have developed a novel analytical model that describes the tidal response of thermally forced atmospheres on rocky planets.The model derivation is based on the secure ground of the first principles of linear atmospheric dynamics, studied under classical approximations that are commonly drawn in earlier analytical works and in more recent numerical frameworks.The distinct feature that we imposed in this model is that of neutral atmospheric stratification, which presents a more realistic description of the Earth's troposphere than the isothermal profile imposed in earlier analytical studies.In this limit, we derive from the model a closed form solution of the tidal torque that can be efficiently used to study the evolution of planetary rotational dynamics.We accommodate into the model dissipative thermal radiation via linear Newtonian cooling, and turbulent and diffusive processes related to thermal inertia budgets in the boundary layer and the ground.As such, the model can be used to study a planetary thermotidal response when heated either by direct synchronous absorption of the incident stellar flux, or by a delayed infrared radiation from the ground.
We probed the spectral behavior of the tidal torque using this developed model in the two aforementioned limits.In the limit of longwave heating flux, the inherently delayed thermal response in the planetary boundary layer maneuvers the tidal bulge in such a way that, for typical values of thermal inertia in the ground and atmosphere, the accelerating effect of the tidal torque at the Lamb resonance is attenuated, and possibly annihilated.In the case of the Earth, -where we apply the opposite limit of shortwave thermotidal heating and ignore the attenuating effect of asynchronous forcing -while the encounter of the resonance in the atmosphere is guaranteed, the epoch of its occurrence and the tidal amplitude it generates are uncertain.As such, we attempted a cautious incursion on constraining them and learned that: • Assuming that temperate climatic conditions have prevailed over the Earth's history, the resonance is likely to have occurred in the early Mesozoic, and not in the Precambrian.The early Mesozoic, unlike the Precambrian, is characterized by an amplified decelerating luni-solar gravitational torque.
• For judiciously constrained estimates of our atmospheric model parameters, the resonance does not amplify the accelerating thermotidal torque to a level comparable in magnitude to the gravitational counterpart.
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 that 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 modelling 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 analogues.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.
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 Gn .First, we combine Eq. (B.1) with the heat transfer equation (B.4) to eliminate the vertical component of velocity, Ṽr;n .Then we use the ideal gas law (B.5) to replace δT n by δp n and δρ n in the resulting equation.It follows that (B.6) Next, we use the hydrostatic equilibrium condition (B.2) to express δρ n in terms of δp n and Ũn in the above equation to obtain Now our aim is to obtain another first order equation in the calculation variable Gn .We start with the hydrostatic equilibrium condition (B.2), and we substitute δp n by its expression (B.1).Then we use the continuity equation (B.3) to express δρ n as a function of Ṽr;n and δp n , and we finally eliminate Ṽr;n using equation (B.1) to obtain (B.8)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 Gn (x) = Φ(x)Ψ n (x) to write Eq. (B.11) in the standard form of wave equations.We find that the first order term in Eq. (B.11) cancels when 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 k defined as With the definitions of the dimensionless control parameters given in Table 1, the wave number kn (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 Gn .
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 (C.1) Given δp n (x), the horizontal component of the velocity field is straightforwardly obtained via (e.g., Eqs.(25-26) of Auclair-Desrotour et al., 2017b): The vertical component of the velocity is established from equation (B.1), replacing δp n by its solution (C.1): Finally, the solution for the density perturbation is obtained by replacing the solution of δp n (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 B σ gr such that δT s = B σ gr δF inc .The variation δT s induces radiative emission with magnitude δF rad .A fraction of the incident power, δQ gr , is then transmitted to the ground by thermal conduction, and another fraction, δQ atm , is transmitted to the atmosphere through turbulent thermal diffusion (see Figure D.7).
Having isolated these thermal mechanisms, the total power balance of the thermal perturbation at the surface interface is expressed as9 : δF inc − δF rad = δQ gr + δQ atm . (D.1) To define each of these terms explicitly, we start with the surface radiative emission term δF rad , which can be obtained by differentiating the Stefan-Boltzmann law as a function of the surface temperature assuming that the planetary surface radiates as a black body.Namely we have with σ SB = 5.670373 × 10 −8 W m −2 K −4 being the Stephan-Boltzmann constant (Tiesinga et al., 2021).
Next, the terms associated with heat diffusion at the boundary interface, δQ gr and δQ atm , are derived using the gradient flux theory (Garratt, 1994): Here, k gr is the thermal conductivity of the ground, and k atm is its atmospheric analogue.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 and 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: (D.8) Here, we remind the reader that s = sgn(σ), and we denote by h σ gr and h σ atm the skin thicknesses of heat transport by thermal diffusion in the ground and the atmosphere defined 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 δF tot = δF inc − δQ gr , (D.15) from which we retrieve Eq.( 10) of the main text.
As for the numerical values of the material properties considered in the formalism above, we first set the gas heat capacity to its nominal value C p = 1005 J/kg•K (Hilsenrath, 1955).The volumetric heat capacity in the ground ρ 0 (0 − )C gr varies between the land and the oceans, but since we're 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 and 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 and 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 and Boville, 1993;Berger and Grisogono, 1998).Specifically, Bernard (1962) computes an average effective value over continental areas K atm = 3.6 m 2 s −1 .This variability of the ground and atmospheric thermal parameters propagates to our model free parameter τ bl , or equivalently the frequency σ bl , via the thermal ground and atmospheric thermal inertias defined in Eq. (D.13).
Thus in sweeping the limits of σ bl in Figure 3, we are moving between strong and weak thermal inertia limits by self-consistently varying the values of the thermal properties.For instance, the value of σ bl = 10 −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 K atm = K atm = 3.6 m 2 s −1 .

Appendix E. The quadrupolar component of the incident stellar flux
As we are interested in the semi-diurnal tidal response, we decompose here the incoming stellar flux harmonically to isolate the quadrupolar component.We start by expressing the day-night periodically varying incident flux δF inc (denoted hereafter by F for convenience) as: 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: (E.4) where the asterisk corresponds to complex conjugation, and the coordinates (θ S , φ S ) = [π/2, (Ω − n ⋆ )t] define the position of the star when ignoring planetary obliquity.We thus have: Finally, F is written as the real part of the thermal forcing function in the Fourier domain: (2 − δ 0m )γ lm (θ)e im(σ 11 t+φ) , (E.11) from which we retrieve the quadrupolar component (l = m = 2) of the forcing in the form F 22 (θ, φ, t) = 5 64 F ⋆ P 22 (cos θ)e i(σ 22 t+2φ) , (E.12) or in spherical harmonics as where σ 22 = 2(Ω − n ⋆ ) is the semi-diurnal tidal forcing frequency.

Appendix F. Energy carried by tidal waves
Essential to the closed-form solution we provide for the tidal response is the boundary condition imposed against tidal waves at the uppermost layer of the atmosphere.There we enforce a nondivergence condition on the tidal flow, namely, the energy of tidal flows W should be bounded as x → ∞.It follows from Wilkes (1949) that the energy flow associated with the tidal (n, m, ν)-mode across the vertical direction is given by By substituting Ṽr 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 re-write 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 and Ingersoll (1980); Eq.( 3) of Zahnle and Walker (1987); Eq.( 4) of Auclair-Desrotour et al. (2019)].Upon a more careful inspection, one can relate the evident mismatch to the adopted convention in defining the pressure anomaly, which when accounted for reveals the overlap between the different expressions.Here we recover the expression of the tidal torque that we use in Eq.( 21).
The instantaneous torque exerted by tidal force about the spin axis of a planet of volume V 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.i.e. 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 S 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 g = GM p /R 2 p , 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;Auclair-Desrotour 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 X and Y are defined by Eq. ( 23), and X R and Y R are defined by Next, we reformulate this complex solution as 3) where we introduced σ L = σ w √ Λ 2 , the frequency which marks the peak of the modulus of the resonant pressure variation in the adiabatic limit, the transfer function of the thermal forcing G(σ) defined as and the dimensionless constant factor K defined as The transfer function G thus describes the delayed forcing of the atmosphere due to thermal inertia near the planetary surface.Therefore, this function reduces to G (σ) = 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 G(σ) defined as such, the complex solution is structured in a convenient form where the effects of atmospheric dynamics and the thermal forcing are clearly separated.This facilitates future modifications of the model if one intends, for instance, to replace the thermal forcing function by another.
Next we define the normalized frequency X = σ/σ L and the dissipation parameter ε = σ 0 /σ L ≪ 1, and we rewrite Eq. (H.3) as . (H.6) The modulus, real part, and the imaginary part of this complex surface pressure oscillation are then straightforwardly deduced as: 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 G = 1 we obtain which is similar to the functional form of the imaginary part of the pressure oscillation given in Ingersoll and 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. 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 and Lindzen (1970) and Vichare and Rajaram (2013) who provide such distributions for water vapour (H 2 O) and ozone (O 3 ).
First, we consider the vertical and latitudinal profiles shown by Chapman and 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 and Lindzen, 1970, Eqs. (143-144)) J m,σ (x, θ, φ) e iσt , (I.1) with the coefficients J m,σ given by J m,σ (x, θ, φ) = iσC p τ m G σ (x, θ, φ) , (I.2) where the heat capacity of the gas per unit mass C p = R s /κ. Figure 3.2 of Chapman and Lindzen (1970) shows the normalised vertical distribution of τ 2,σ G .Following the authors (their Eq. 163), we set H = 7.6 km.The semidiurnal distributions of τ 2,σ G for water vapour and ozone can be approximated by the functions where k H 2 O ≈ 1/3.Similarly, the vertical profile of ozone absorption can be approximated by a truncated sine function, x 2 −x 1 , x 1 ≤ x ≤ x 2 , 0 elsewhere, (I.6) with K O 3 = 1.4,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 δF (θ, φ) = δF 22 Y 22 (θ, φ) . (I.7) The efficiency parameter α A of our model is basically the ratio between δF 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 cooling frequency is also a function of the surface temperature and it scales as σ 0 ∝ p −1 s T 3 s (e.g., Eq. 10 of Showman and 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. (K.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 L. The resonant period in the isothermal limit Here we discuss the limit where the temperature profile of the atmosphere is modelled 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 and Blake, 1972;Auclair-Desrotour et al., 2019).We are specifically interested in estimating the difference between the two models in terms of the Lamb resonance period.Namely, we aim to define the spectral position of the resonance in the isothermal limit and compare with what we obtained for our model in Eq. 24.
We adopt the analytical computations in Auclair-Desrotour et al. ( 2019) that yield a frequency-dependence for the pressure surface anomaly in the form (their Eq. 12): 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 and 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 T = T s /(κ + 1), thus we rewrite Eq.L.3 as: LOD res; iso = 4πR p R s Λ n T s /(1 − κ 2 ) + 2R p n ⋆ . (L.4) This equation yields a resonant rotational period that is roughly one hour smaller than that obtained in the neutrally stratified case; i.e., the curve shown in Figure 6 shifts vertically downwards 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 and Walker (1987).

Figure 1 :
Figure1: Modeled histories of the rotational motion of the Earth.Plotted is the Earth's LOD evolution in time over geological timescales for three models: i) the model ofFarhat et al. (2022), where the evolution is driven solely by oceanic and solid tidal dissipation; ii) the model ofZahnle and Walker (1987), where the Lamb resonance is encountered for LOD∼21 hr, forcing a rotational equilibrium on the Earth; iii) the model ofBartlett and 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 modelled histories are geological proxies of the LOD evolution that can be retrieved from http://astrogeo.eu.

Figure 2 :
Figure 2: The spectrum of semi-diurnal atmospheric thermal tides.Plotted is the imaginary part of the normalized pressure anomaly ( δp = δp/p s ; Eq.22) as a function of the normalized forcing frequency ω = (Ω − n ⋆ )/n ⋆ = σ/2n ⋆ , where n ⋆ is the mean motion of the stellar perturber.The planetary-stellar parameters are those of the fiducial planetary system defined in Section 3.1.

Figure 3 :
Figure 3: The (a-)symmetry of the Lamb resonance.Similar to Figure 2, plotted is the imaginary part of the normalized pressure anomaly (Eq.22), associated with the semidiurnal tide, as a function of the normalized forcing frequency ω = (Ω − n ⋆ )/n ⋆ = σ/2n ⋆ , for the same planetary-stellar 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.

Figure 4 :
Figure 4: The 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 mass-based 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.

Figure 5 :
Figure 5: A 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 ℑ{δ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 (ℑ{δp s } = 2275 Pa) needed in the early Mesozoic, 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 semi-diurnal tidal bulge.The gray area on the left covers the parametric region where the resonance features a lower pressure amplitude than the present.

Figure 6 :
Figure 6: The 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.The gray shaded area highlights 95% confidence intervals for the past temperature evolution according to the carbon cycle model ofKrissansen-Totton et al. (2018).The identified geological eras correspond to the LOD evolution model ofFarhat et al. (2022).The overlap between the modelled temperature evolution and the black curve places the resonance occurrence in the early Mesozoic.

Figure
Figure D.7: A 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).
sin θe i2φ , (I.3)τ G,O 3 (x, θ, φ) = T O 3 V O 3 (x)sin 2 θe i2φ , (I.4) with T H 2 O = 0.035 K and T O 3 = 0.25 K.As the vertical profile of water vapour shown by Figure 3.2 of Chapman and Lindzen (1970) is an exponentially decaying function of the altitude, the function V H 2 O of Eq.(I.3) is set to V H 2 O (x) = e −k H 2 O x , (I.5) Figure J.9: Similar to Figure 2, plotted are the tidal spectra, in terms of the imaginary part of the surface pressure anomaly, as a function of the normalized tidal frequency ω.The blue curve is that computed numerically in Auclair-Desrotour 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 J.

Table 1 :
Dimensionless control parameters determining the regime of the atmospheric tidal response and defined throughout the text.