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

Mohammad Farhat<EMAIL_ADDRESS>Pierre Auclair-Desrotour  Gwenaël
Boué  Russell Deitrick  Jacques Laskar

###### Abstract

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 examine 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 the opposing torque of the oceans and solid
interior, driving the Earth into a rotational equilibrium. With this scenario
in mind, we construct 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 retrieve 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.

###### keywords:

Atmospheric dynamics , Thermal tides , Earth’s rotation , Precambrian Earth

††journal: Earth and Planetary Science Letters

[1]organization=IMCCE, CNRS, Observatoire de Paris, PSL University, Sorbonne
Université,addressline=77 Avenue Denfert-Rochereau, city=Paris,
postcode=75014, country=France [2]organization=School of Earth and Ocean
Sciences, University of Victoria, city= Victoria, British Columbia,
country=Canada

##  1 Introduction

For present day Earth, the semi-diurnal atmospheric tide, driven by the
thermal forcing of the Sun and generated via atmospheric pressure waves,
describes the movement of atmospheric mass away from the substellar point.
Consequently, mass culminates forming bulges on the nightside and the dayside,
generating a torque that accelerates the Earth’s rotation. As such, this
thermally generated torque counteracts the luni-solar gravitational torque
associated with the Earth’s solid and oceanic tides. The latter components
typically drive the closed system of the tidal players 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 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 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 and Laskar, 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).

Figure 1: 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 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${\sim}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.

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 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 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
${\sim}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
and Laskar, 2001, 2003a), to fully numerical efforts that relied on the
advancing sophistication of general circulation models (GCM; e.g., Leconte et
al., 2015; Auclair-Desrotour et al., 2019). The latter work presents, to our
knowledge, the first and, to-date111While this paper was under review, Wu et
al. (2023) presented another GCM-computed spectrum for the Earth in the high
frequency regime. We provide an elaborate discussion of their work’s results
in a separate dedicated paper (Laskar et al., 2023)., 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 piece-wise 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 $\sim 19.5$ hr between 0.6 and 2.2 Ga,
with a sustained very high mean surface temperature ($40-55^{\circ}$C).
Although using different approaches, the two studies have thus arrived at
similar conclusions. A closer look at the subset of geological data that
favored this outcome, however, indicates that both studies heavily rely on
three stromatolitic records from the Paleoproterozoic that were originally
studied by Pannella (1972a, b). These geological data have been, ever since,
identified as unsuitable for precise quantitative interpretation (see e.g.,
Scrutton, 1978; Lambeck, 1980; Williams, 2000). To this end, we provide a more
detailed analysis of the geological proxies of the LOD, and of the model
presented by Wu et al. (2023) in a parallel paper dedicated to the matter
(Laskar et al., 2023).

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

##  2 Ab initio atmospheric dynamics

For an atmosphere enveloping a spherically symmetric planet, we define a
reference frame co-rotating with the planet. In this frame, an atmospheric
parcel is traced by its position vector $\boldsymbol{r}$ in spherical
coordinates ($r,\theta,\varphi)$, such that $\theta$ is the colatitude,
$\varphi$ is the longitude, and the radial distance $|\boldsymbol{r}|=R_{\rm
p}+z$, where $R_{\rm 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 $\rho$, and the three-dimensional
vectorial velocity field $\boldsymbol{\mathcal{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}+\delta p,\,T=T_{0}+\delta T,\,\rho=\rho_{0}+\delta\rho,\,$and
$\boldsymbol{\mathcal{V}}=\boldsymbol{V}_{0}+\boldsymbol{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
hypothesis222noting that surface friction is another dissipative mechanism as
discussed by Lindzen and Blake (1972). (Lindzen and McKenzie, 1967), where
radiative losses, $J_{\rm rad}$, are parameterized by the characteristic
frequency $\sigma_{0}$; namely $J_{\rm
rad}=p_{0}\sigma_{0}/(\kappa\rho_{0}T_{0})\delta T$, where
$\kappa=(\Gamma_{1}-1)/\Gamma_{1}=0.285$ and $\Gamma_{1}$ is the adiabatic
exponent. Similar to Leconte et al. (2015), we associate with $\sigma_{0}$ a
radiative cooling timescale $\tau_{\rm rad}=4\pi/\sigma_{0}.$

###  2.1 The vertical structure of tidal dynamics

We are interested in providing a closed form solution for the frequency333The
frequency in this case being the tidal forcing frequency $\sigma$, typically a
linear function of the planet’s spin rate $\Omega$ and the stellar perturber’s
mean motion $n_{\star}$. The semi-diurnal tidal frequency, for instance,
$\sigma_{22}=2(\Omega-n_{\star}).$ dependence of the thermotidal torque, which
results from tidally driven atmospheric mass redistribution. By virtue of the
hydrostatic approximation, this mass redistribution is encoded in the vertical
profile of pressure. As such, it is required to solve for the vertical
structure of tidal dynamics. With fellow non-theoreticians in mind, we
delegate the detailed development of the governing system of equations
describing the tidal response of the atmosphere to 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 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):

$\mathcal{L}^{m,\nu}\Theta^{m,\nu}_{n}=-\Lambda^{m,\nu}_{n}\Theta^{m,\nu}_{n},$
(1)

where the set of Hough functions $\\{\Theta^{m,\nu}_{n}\\}$ serves as the
solution (Hough, 1898), $\\{\Lambda^{m,\nu}_{n}\\}$ is the associated set of
eigenvalues, $\mathcal{L}^{m,\nu}$ is a horizontal operator defined in Eq.(A)
of A, while $\nu=2\Omega/\sigma$, where $\Omega$ is the rotational velocity of
the planet and $\sigma$ is the tidal forcing frequency. In the tidal system
under study, the variables and functions $(\delta p,\delta\rho,\delta
T,\boldsymbol{V},J,\Theta,\Lambda)$ are written in the Fourier domain using
the longitudinal order $m$ and frequency $\sigma$ (Eq.44), and expanded in
horizontal Hough modes with index $n$ (Eq.45). We denote hereafter their
coefficients $f_{n}^{m,\nu}$ 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 $\\{\Lambda^{m,\nu}_{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 B. The
outcome is a wave-like equation that describes vertical thermotidal dynamics
and reads as:

$\frac{d^{2}\Psi_{n}}{dx^{2}}+\hat{k}_{n}^{2}\Psi_{n}=\Phi^{-1}C_{n}.$ (2)

Here, as is the common practice (e.g., Siebert, 1961; Chapman and Lindzen,
1970), we use the reduced altitude $x=\int_{0}^{z}dz/H(z)$ as the vertical
coordinate, where the pressure scale height $H(z)=\mathcal{R}_{\rm
s}T_{0}(z)/g$; $\mathcal{R}_{\rm s}$ being the specific gas constant and $g$
the gravitational acceleration. The quantity $\Psi_{n}(x)$ is a calculation
variable from which, once solved for, all the tidal scalar and vectorial
quantities would flesh out (C). The vertical wave number $\hat{k}_{n}(x)$ is
defined via

$\displaystyle\hat{k}_{n}^{2}(x)$
$\displaystyle=-\frac{1}{4}\Bigg{\\{}\\!\\!\\!\left(\\!1\\!-\\!\frac{i\kappa}{\alpha-i}(\gamma-1)\right)^{2}\\!\\!\\!\\!+2\frac{d}{dx}\\!\\!\left(\frac{i\gamma}{\alpha-i}\right)$
$\displaystyle-\frac{4\alpha\kappa}{\alpha-i}\\!\left[\beta\gamma\Lambda_{n}-\\!\frac{i}{\alpha}(\gamma-1)\right]\Bigg{\\}}.$
(3)

Table 1:  Dimensionless control parameters determining the regime of the atmospheric tidal response and defined throughout the text.  Parameter | Scale significance of
---|---
$\eta={R_{\rm p}}/{H}$ | horizontal structure of dynamics
$\nu={2\Omega}/{\sigma}$ | Coriolis effects
$\alpha={\sigma}/{\sigma_{0}}$ | radiative cooling
$\beta={\sigma_{\rm w}^{2}}/{\sigma^{2}}$ | Lamb waves
$\gamma={N_{\rm B}^{2}}/{N_{\rm B;iso}^{2}}$ | buoyancy forces
$\zeta=\sqrt{{|\sigma|}/\sigma_{\rm bl}}$ | ground and atmospheric inertia
$\mu_{\rm gr}=\frac{I_{\rm gr}}{I_{\rm gr}+I_{\rm atm}}$ | ground thermal inertia

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

$N_{\rm B}^{2}=\frac{g}{H}\left(\kappa+\frac{d\ln H}{dx}\right).$ (4)

We define the Brunt–Väisälä frequency in the limit of an isothermal
atmospheric profile, $N_{\rm B;iso}=\sqrt{\kappa g/H}$. As such, the parameter
$\gamma$ 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 $\Phi$ defined as

$\Phi(x)=\exp\left[\frac{x}{2}-\frac{1}{2}\int_{0}^{x}\frac{i}{\alpha-i}\left(\frac{HN_{\rm
B}^{2}}{g}-\kappa\right)dx^{\prime}\right],$ (5)

and the forcing function $C_{n}(x)$ which reads:

$\displaystyle C_{n}(x)$
$\displaystyle=\frac{\alpha\beta\Lambda_{n}\tilde{J}_{n}}{\alpha-i}+\Bigg{[}\frac{\kappa(\gamma-1)}{\alpha-i}\left(\Lambda_{n}+\frac{\partial_{x}}{\beta}\right)$
$\displaystyle-i\partial_{x}\left(\frac{\partial_{x}}{\beta}\right)\\!\Bigg{]}\tilde{U}_{n}.$
(6)

Hereafter, we identify the dimensionless form of the tidal variables by the
tilde diacritic. Namely in Eq.(6), the dimensionless thermal forcing function
$\tilde{J}_{n}=\kappa/(g\sigma H)J_{n}$, while the gravitational analogue
$\tilde{U}_{n}=U_{n}/v_{0}^{2}$, where the reference velocity $v_{0}=\sigma
R_{\rm p}$ (see Eq.(35) 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, $\Psi_{n}(x)$. In C, we derive the vertical profiles of all the
tidal variables, and specifically for the dimensionless pressure anomaly we
obtain:

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

where $\tilde{\delta p}_{n}(x)=\delta p_{n}(x)/p_{0}$, and the calculation
variable $\tilde{G}_{n}(x)=\Psi_{n}(x)\Phi(x)$.

###  2.2 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 $\tilde{J}_{n}$ (or
equivalently in dimensional form, $J_{n}$). We adopt a vertical tidal heating
profile of the form

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

where $J_{\rm s}$ is the heat absorbed at the surface and $b_{\rm 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_{\rm
J}\rightarrow\infty$, and a uniform distribution where the whole air column is
uniformly heated ($b_{\rm J}=0)$.

To determine $J_{\rm s}$, we invoke its dependence on the total vertically
propagating flux $\delta F_{\rm 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
$\alpha_{\rm A}$, where $0\leq\alpha_{\rm A}\leq 1$; the rest of the flux
amounting to $1-\alpha_{\rm A}$ escapes the thermotidal interplay. We thus
have

$\int_{0}^{\infty}J(x)\rho_{0}(x)H(x)dx=\alpha_{\rm A}\delta F_{\rm tot}.$ (9)

To define $\delta F_{\rm tot}$, we establish the flux budget for a small
thermal perturbation at the planetary surface. We start with $\delta F_{\rm
inc}$, a variation of the effective incident stellar flux, after the reflected
component has been removed. $\delta F_{\rm inc}$ generates a variation $\delta
T_{\rm s}$ in the surface temperature $T_{\rm s}$. The proportionality between
$\delta F_{\rm inc}$ and $\delta T_{\rm s}$ is parameterized by $\tau_{\rm
bl}$, a characteristic diffusion timescale of the ground and atmospheric
surface thermal responses. We detail on this proportionality in D, but for now
it suffices to state that $\tau_{\rm bl}$ is a function of the thermal inertia
budgets in the ground, $I_{\rm gr}$, and the atmosphere $I_{\rm atm}$. We
associate with $\tau_{\rm bl}$ the frequency $\sigma_{\rm bl}=\tau_{\rm
bl}^{-1}$, a characteristic frequency that reflects the thermal properties of
the diffusive boundary layer. It will serve as another free parameter of our
tidal model, besides the Newtonian cooling frequency $\sigma_{0}$, and the
atmospheric opacity parameter $\alpha_{\rm A}$. In analogy to
$\alpha=\sigma/\sigma_{0}$, we define the dimensionless parameter for the
boundary layer $\zeta=\sqrt{|\sigma|\tau_{\rm bl}}=\sqrt{|\sigma|/\sigma_{\rm
bl}}$.

By virtue of the power budget balance established in D, we define the total
propagating flux $\delta F_{\rm tot}$ as

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

Here, $s={\rm sign}(\sigma)$, and $\mu_{\rm gr}$ is a dimensionless
characteristic function weighing the relative contribution of ground thermal
inertia to the total inertia budget; namely $\mu_{\rm gr}={I_{\rm
gr}}/({I_{\rm gr}+I_{\rm 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 $\zeta=\mu_{\rm 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 ($\zeta\neq 0$,
$\mu_{\rm gr}\neq 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 $\tau_{\rm 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 E to obtain the amplitude of the quadrupolar component
as $\delta F_{\rm inc}=\delta F_{22}=({\sqrt{30\pi}}/{16})F_{\star}$ , where
$F_{\star}=L_{\star}/4\pi a_{\rm p}^{2}$, $L_{\star}$ being the stellar
luminosity, and $a_{\rm p}$ the star-planet distance.

##  3 The tidal response

###  3.1 The tidal torque in the neutral stratification limit

Under the defined forcing in the previous section, to solve the wave equation
analytically, a choice has to be made on the Brunt–Väisälä frequency, $N_{\rm
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_{\rm B}^{2}=\kappa g/H={\rm const}.$

Motivated by the Earth’s atmosphere, where the massive troposphere
(${\sim}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_{\rm B}^{2}=0.$ In fact, $N_{\rm B}^{2}$ can be expressed
in terms of the potential temperature $\Theta_{0}$ (e.g., Section 2.10 of
Vallis, 2017):

$N_{\rm B}^{2}=\frac{g}{H}\frac{d\ln\Theta_{0}}{dx}.$ (11)

whereby the stability of the atmosphere is controlled by the slope of
$\Theta_{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 (A).

Hereafter, we focus on the thermotidal heating as the only tidal perturber,
and we ignore the much weaker gravitational potential $\tilde{U}$. It follows,
in the neutral stratification limit, that $\gamma=0$ (Table 1), and the
vertical wavenumber (Eq. 3) reduces to444It is noteworthy that the wavenumber
in the neutral stratification limit is not longer dependent on the horizontal
structure.

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

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

$\displaystyle p_{0}(x)$ $\displaystyle=p_{0}(0)e^{-x},\hskip
14.22636pt\rho_{0}(x)=\frac{p_{0}(0)}{gH(0)}e^{(\kappa-1)x},\hskip 8.5359pt$
$\displaystyle T_{0}(x)$ $\displaystyle=\frac{gH(0)}{\mathcal{R}_{\rm
s}}e^{-\kappa x}.$ (13)

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

$J_{\rm s}=\delta F_{22}\frac{\alpha_{\rm A}g(b_{\rm
J}+1)}{p_{0}(0)}\\!\\!\left[1-\mu_{\rm
gr}\zeta\frac{1+si}{1+(1+si)\zeta}\right].$ (14)

As such, the wave equation (2) is rewritten as

$\frac{d^{2}\Psi_{n}}{dx^{2}}+\hat{k}^{2}\Psi_{n}=\mathcal{A}_{n}e^{-\mathcal{B}x},$
(15)

where the complex functions $\mathcal{A}_{n}$ and $\mathcal{B}$ are defined
as:

$\displaystyle\mathcal{A}_{n}=\frac{\kappa\Lambda_{n}}{R_{\rm
p}^{2}\sigma^{3}}\frac{\alpha^{2}+i\alpha}{\alpha^{2}+1}J_{\rm s},$ (16)
$\displaystyle\mathcal{B}=\frac{1}{2(\alpha^{2}+1)}\left[(2b_{\rm
J}+1)(\alpha^{2}+1)-\kappa+i\kappa\alpha\right].$ (17)

The wave equation (15) admits the general solution

$\Psi_{n}(x)=c_{1}e^{i\hat{k}x}+c_{2}e^{-i\hat{k}x}+\frac{\mathcal{A}_{n}}{\mathcal{B}^{2}+\hat{k}^{2}}e^{-\mathcal{B}x}.$
(18)

We consider the following two boundary conditions:

  * 1.

First, the energy of tidal flows, $\mathcal{W}$, should be bounded as
$x\\!\rightarrow\\!\infty$. In F, we derive the expression of the tidal energy
following Wilkes (1949), and it scales as
$\mathcal{W}\propto|\Psi|^{2}|\Phi|^{2}$. Accordingly, the non-divergence of
the flow condition is satisfied if one sets $c_{2}=0$ and takes the proper
sign of the wavenumber (Eq. 12), namely:

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

  * 2.

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

$c_{1}(x)=\frac{\mathcal{A}_{n}}{\mathcal{B}+{\hat{k}^{2}}}\times\frac{\mathcal{B}-\frac{1}{2}\left(1+\frac{i\kappa}{\alpha-i}\right)-\beta\Lambda_{n}+1}{i\hat{k}+\frac{1}{2}\left(1+\frac{i\kappa}{\alpha-i}\right)+\beta\Lambda_{n}-1}.$
(20)

Under these boundary conditions, we are now fully geared to analytically
compute the solution of the wave equation, $\Psi_{n}(x)$ (or equivalently
$\tilde{G}_{n}(x))$, but we are specifically interested in retrieving a closed
form solution of the quadrupolar tidal torque. The latter takes the general
form555We 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.  (G):

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

Here $M_{\star}$ and $M_{\rm p}$ designate the stellar and planetary masses
respectively, and $\Im$ refers to the imaginary part of a complex number, the
latter in this case being the quadrupolar pressure anomaly at the surface
$\delta p_{\rm s}=\delta p_{2}^{2,\nu}(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 $\Psi(0)$ of Eq. (18) (with $c_{2}=0$ and $c_{1}$ defined
in Eq. 20), we retrieve $\delta p_{\rm 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:

$\displaystyle\Im\\{\delta p_{\rm s}\\}=$ $\displaystyle\alpha_{\rm A}\delta
F_{22}\frac{\kappa g\Lambda_{2}}{R_{\rm
p}^{2}\sigma^{3}}\frac{(\mathcal{X}\alpha+\mathcal{Y})\alpha}{\left(1+2\zeta+2\zeta^{2}\right)}$
$\displaystyle\times{\underbrace{\left[(\kappa-\beta\Lambda_{2}+1)^{2}+\alpha^{2}(\beta\Lambda_{2}-1)^{2}\right]}_{\text{position
of the Lamb resonance}}}^{-1},$ (22)

where we have defined the complex functions $\mathcal{X}$ and $\mathcal{Y}$ as

$\displaystyle\mathcal{X}$
$\displaystyle=(\beta\Lambda_{2}-1)\left[2\zeta^{2}(1-\mu_{\rm
gr})+\zeta(2-\mu_{\rm gr})+1\right],$ $\displaystyle\mathcal{Y}$
$\displaystyle=-s\mu_{\rm gr}\zeta(\kappa-\beta\Lambda_{2}+1).$ (23)

We note that we provide the full complex transfer function of the surface
pressure anomaly, along with further analysis on its functional form in H.
Before embarking on any results, we pause here for a few remarks on the
provided closed form solution of the torque.

Figure 2: The spectrum of semi-diurnal atmospheric thermal tides. Plotted is
the imaginary part of the normalized pressure anomaly ($\tilde{\delta
p}=\delta p/p_{\rm s};$ Eq.3.1) as a function of the normalized forcing
frequency $\omega=(\Omega-n_{\star})/n_{\star}=\sigma/2n_{\star}$, where
$n_{\star}$ is the mean motion of the stellar perturber. The planetary-stellar
parameters are those of the fiducial planetary system defined in Section 3.1.

  * 1.

The parameter $\alpha_{\rm 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, $\alpha_{\rm 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 $\alpha_{\rm A}$ in
Section 4.1 and I.

  * 2.

The quadrupolar component of the equilibrium stellar flux, entering through a
fraction of $F_{\star}$ (E), is directly proportional to the stellar
luminosity $L_{\star}$. Standard models suggest that the Sun’s luminosity was
around 80% of its present value ${\sim}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.

  * 3.

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, $\Lambda_{n}$, in our case
reducing to the dominant fundamental mode $\Lambda_{2}$. Noting that we have
dropped the superscripts, we remind the reader that for the semidiurnal
($m=2)$ response, $\Lambda_{2}=\Lambda^{2,\nu}_{2}$, thus $\Lambda$ is
frequency-dependent in the general case. The Earth however, over its lifetime,
lives in the asymptotic regime of $\nu\approx 1$ since $2\Omega\gg n_{\star}$,
thus it is safe to assume that $\Lambda_{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).

  * 4.

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 (3.1) that the position of
the resonance is controlled by the highlighted term. Had it not been for the
introduced radiative losses, entering here through $\alpha$, this term would
have encountered a singularity at the spectral position of the resonance, i.e.
for $\beta\Lambda_{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_{\rm p}=0.815M_{\oplus},\,R_{\rm
p}=0.95R_{\oplus},\,a_{\rm p}=0.73\,{\rm au},\,g=8.87\,{\rm 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_{\star}=1M_{\odot},\,L_{\star}=1L_{\odot})$.
We further ignore here the thermal inertia in the ground and the atmosphere by
taking $\sigma_{\rm bl}\rightarrow\infty$, or $\zeta\rightarrow 0$, thus
assuming an synchronous response of the ground with the thermal excitation.

Figure 3: The (a-)symmetry of the Lamb resonance. Similar to Figure 2, plotted
is the imaginary part of the normalized pressure anomaly (Eq.3.1), associated
with the semidiurnal tide, as a function of the normalized forcing frequency
$\omega=(\Omega-n_{\star})/n_{\star}=\sigma/2n_{\star}$, 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
$\sigma_{\rm 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.

Key tidal response features are recovered in this spectrum: First, we obtain a
tidal peak near synchronization $(\omega=0)$ that generates a positive torque
for $\sigma>0$ and a negative torque for $\sigma<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.

###  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, $\tau_{\rm bl}$ (Eq.84), or
equivalently $\sigma_{\rm bl}$. Variations in $\sigma_{\rm 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.

Figure 4: The thermally induced tidal bulge revealed. Shown are polar
snapshots of the radial and longitudinal variations of the tidal pressure
anomaly $\delta 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 $\varsigma=p/p_{\rm s}$. The longitudinal axes are
shown in increments of 30∘ with $0^{\circ}$ 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
$\omega=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 $\omega$, from lagging to leading the substellar point.
In contrast, in the bottom row, and for the prescribed value of $\sigma_{\rm
bl}$, the delayed response of the ground forces the bulge to always lag the
substellar point, thus acting to decelerate the planetary rotation.

In such a case, the value of $\sigma_{\rm bl}$ plays a significant role in the
tidal response of the planet. Namely, the ratio $\sigma/\sigma_{\rm 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
$\sigma\ll\sigma_{\rm 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
$\sigma\gtrsim\sigma_{\rm bl})$. We now reap from the analytical model the
signature of $\sigma_{\rm 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
$\sigma_{\rm bl}$. For $\sigma_{\rm bl}=5\times 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 $\sigma_{\rm 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 $\sigma_{\rm bl}$ to $10^{-5}\,s^{-1}$ in the panel in the middle,
the positive pressure peak completely diminishes, leaving only the negative
counterpart. Decreasing $\sigma_{\rm 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_{\Omega}$, due to
solid or oceanic gravitational tides, $\omega$ 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_{\Omega}$-pumping part of the atmospheric torque is attenuated, and for
some values of $\sigma_{\rm 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 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 vertical profile of the pressure anomaly
from Eq. (7), and augment it with the latitudinal and longitudinal
dependencies from Eqs. (44-45). As the massive troposphere dominates the tidal
mass redistribution, we use the mass-based vertical coordinate
$\mathcal{\varsigma}=p/p_{\rm s}$ (i.e. $x=-\ln\varsigma$, and $\varsigma$
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 $\sigma_{\rm 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 $\omega=253$, i.e. before encountering the
resonance, the bulge leads the substellar point and acts to accelerate the
planet’s rotation. Increasing $\omega$ 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 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
$\sigma_{\rm 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 small-amplitude non-migrating tidal components they
produce (e.g., Schindelegger and Ray, 2014).

##  4 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.

###  4.1 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, 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 ${\sim}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 $\zeta\rightarrow 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 $\sigma_{0}$ and
$\alpha_{\rm 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 $\sigma_{0}\gtrsim 10^{-4}\,s^{-1}$, we severely
suppress the amplified pressure response around the resonance. Conversely, for
values of $\sigma_{0}\lesssim 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.

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 $\sigma_{0}$ and $\alpha_{\rm A}$. The solid black isoline marks
the level curve of $\Im\\{\delta p_{\rm s}\\}=880$ Pa, and defines from below
a region in ($\alpha_{\rm A},\sigma_{0})$-space where the thermotidal response
is sufficient to cancel the gravitational counterpart in the Precambrian.
Analogously, the dashed isoline defines the threshold ($\Im\\{\delta p_{\rm
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.

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 (${\sim}1.13\times
10^{16}$ N m) is roughly a quarter of the present gravitational torque
(${\sim}4.51\times 10^{16}$ N m) (e.g., Zahnle and Walker, 1987; Farhat et
al., 2022), thus requiring, via Eq.(21), $\Im\\{\delta p_{\rm s}\\}$ on the
order of 880 Pa666This 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 E and G is normalized by a
factor of $\sqrt{16\pi/15}$. . 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 ($\Im\\{\delta p_{\rm 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 $\Im\\{\delta p_{\rm 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 $\sigma_{0}$ and $\alpha_{\rm A}$ can place the Earth’s
atmosphere in the identified regions. Stringent constraints on $\sigma_{0}$
are hard to obtain for the Earth since $\sigma_{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 $\sigma_{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 $\tau_{\rm 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 values777We note that in Leconte
et al. (2015), the tidal frequency under study is that of the diurnal
component, thus we multiply their $\omega_{0}$ value by 2; i.e.
$\sigma_{0}=2\omega_{0}$..

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 $\alpha_{\rm A}{\sim}14\%$. In I, we compute estimates
of the present value of $\alpha_{\rm 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 $\alpha_{\rm A}{\sim}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 $\sigma_{0}$ over a
wide interval, which prohibits the determination of the present value of
$\sigma_{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\times 10^{15}$ N m) needs to
be amplified by a factor of 3.9 to reach the absolute minimum of the opposing
gravitational torque888subject to the uncertainty of the present measurement
of the semi-diurnal surface pressure oscillation discussed earlier 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 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 constraints 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 ($\sigma_{0}\propto T_{0}^{3}$; Auclair-Desrotour et al.,
2017a, Eq. 17). As such, a warmer planet in the past would yield a shorter
cooling timescale, and consequently, more efficient damping of the resonant
amplitude (see 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 $\alpha_{\rm A}=24\%$
can indeed place the response beyond the Precambrian threshold for some values
of $\sigma_{0}$. An increase to four times the present value of $\alpha_{\rm
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 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?

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

###  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.(3.1). The latter is a function of $\sigma$ 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_{\rm s}$, and is independent of
$\sigma_{\rm 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_{\rm s}$. The apparent single curve is actually a
bundle of curves with different values of $\sigma_{0}$, but the effect of the
latter is unnoticeable (if one varies $\sigma_{0}$ by two orders of magnitude,
the resonant rotational period varies by few minutes). As such, the resonant
frequency is predominantly controlled by $T_{\rm s}$, which allows us to take
the adiabatic limit of Eq.(3.1), and straightforwardly derive the tidal
frequency that minimizes the denominator. In terms of the rotational period,
the position of the resonance then reads:

${\rm LOD_{res}}=\frac{4\pi R_{\rm p}}{\sqrt{\mathcal{R_{\rm
s}}\Lambda_{n}T_{\rm s}}+2R_{\rm p}n_{\star}}.$ (24)

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^{\circ}$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\,{\rm hr}$. This value
is in agreement with the $11.38\pm 0.16\,{\rm hr}$ 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 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\,{\rm 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 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 ${\sim}21\,{\rm hr}$, but
this requires an increase in the equilibrium temperature of at least
$55^{\circ}$C, which is inconsistent with the studies mentioned above. Such an
increase in temperature would also increase $\sigma_{0}$ by almost $19\%$ (as
we discuss in the previous section, $\sigma_{0}\propto T_{0}^{3}$; see also
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.

##  5 Summary and Outlook

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

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

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

  * 1.

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.

  * 2.

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 ($\geq 55^{\circ}$C), a long enough
radiative cooling timescale ($\geq 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.

## Acknowledgments

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

## References

  * Abramowitz et al. (1988) Abramowitz, M., Stegun, I.A., Romer, R.H., 1988.  Handbook of mathematical functions with formulas, graphs, and mathematical tables.
  * Andújar Márquez et al. (2016) Andújar Márquez, J.M., Martínez Bohórquez, M.Á., Gómez Melgar, S., 2016.  Ground thermal diffusivity calculation by direct soil temperature measurement. application to very low enthalpy geothermal energy systems.  Sensors 16, 306\.
  * Arkhangelskaya and Lukyashchenko (2018) Arkhangelskaya, T., Lukyashchenko, K., 2018\.  Estimating soil thermal diffusivity at different water contents from easily available data on soil texture, bulk density, and organic carbon content.  Biosystems Engineering 168, 83–95.
  * Arras and Socrates (2010) Arras, P., Socrates, A., 2010\.  Thermal tides in fluid extrasolar planets.  The Astrophysical Journal 714, 1\.
  * Auclair-Desrotour et al. (2017a) Auclair-Desrotour, P., Laskar, J., Mathis, S., 2017a.  Atmospheric tides in earth-like planets.  Astronomy & Astrophysics 603, A107.
  * Auclair-Desrotour et al. (2017b) Auclair-Desrotour, P., Laskar, J., Mathis, S., Correia, A., 2017b.  The rotation of planets hosting atmospheric tides: from venus to habitable super-earths.  Astronomy & Astrophysics 603, A108.
  * Auclair-Desrotour and Leconte (2018) Auclair-Desrotour, P., Leconte, J., 2018\.  Semidiurnal thermal tides in asynchronously rotating hot jupiters.  Astronomy & Astrophysics 613, A45.
  * Auclair-Desrotour et al. (2019) Auclair-Desrotour, P., Leconte, J., Mergny, C., 2019.  Generic frequency dependence for the atmospheric tidal torque of terrestrial planets.  Astronomy & Astrophysics 624, A17.
  * Bartlett and Stevenson (2016) Bartlett, B.C., Stevenson, D.J., 2016\.  Analysis of a precambrian resonance-stabilized day length.  Geophysical Research Letters 43, 5716–5724.
  * Berger and Grisogono (1998) Berger, B.W., Grisogono, B., 1998\.  The baroclinic, variable eddy viscosity ekman layer.  Boundary-layer meteorology 87, 363–380.
  * Bernard (1962) Bernard, E.A., 1962.  Théorie des oscillations annuelles et diurnes de la température à la surface des continents et des océans.  Archiv für Meteorologie, Geophysik und Bioklimatologie, Serie A 12, 502–543.
  * Blake et al. (2010) Blake, R.E., Chang, S.J., Lepland, A., 2010.  Phosphate oxygen isotopic evidence for a temperate and biologically active archaean ocean.  Nature 464, 1029–1032.
  * Catling and Zahnle (2020) Catling, D.C., Zahnle, K.J., 2020\.  The archean atmosphere.  Science advances 6, eaax1420.
  * Chapman and Lindzen (1970) Chapman, S., Lindzen, R.S., 1970\.  Atmospheric tides: thermal and gravitational. volume 15.  Springer Science & Business Media.
  * Charnay et al. (2017) Charnay, B., Le Hir, G., Fluteau, F., Forget, F., Catling, D.C., 2017.  A warm or a cold early earth? new insights from a 3-d climate-carbon model.  Earth and Planetary Science Letters 474, 97–109.
  * Charnay et al. (2020) Charnay, B., Wolf, E.T., Marty, B., Forget, F., 2020\.  Is the faint young sun problem for earth solved?  Space Science Reviews 216, 1–29.
  * Correia and Laskar (2001) Correia, A., Laskar, J., 2001\.  The four final rotation states of venus.  Nature 411, 767–770.
  * Correia and Laskar (2003a) Correia, A.C., Laskar, J., 2003a.  Different tidal torques on a planet with a dense atmosphere and consequences to the spin dynamics.  Journal of Geophysical Research: Planets 108\.
  * Correia and Laskar (2003b) Correia, A.C., Laskar, J., 2003b.  Long-term evolution of the spin of venus: Ii. numerical simulations.  Icarus 163, 24–45.
  * Correia et al. (2003) Correia, A.C., Laskar, J., de Surgy, O.N., 2003.  Long-term evolution of the spin of venus: I. theory.  Icarus 163, 1–23.
  * Correia and Laskar (2010) Correia, A.C.M., Laskar, J., 2010\.  Tidal Evolution of Exoplanets, in: Exoplanets. Tucson, AZ: University of Arizona Press, pp. 239–266.
  * Covey et al. (2014) Covey, C., Dai, A., Lindzen, R.S., Marsh, D.R., 2014\.  Atmospheric tides in the latest generation of climate models.  Journal of the Atmospheric Sciences 71, 1905–1913.
  * Cunha et al. (2015) Cunha, D., Correia, A.C., Laskar, J., 2015.  Spin evolution of earth-sized exoplanets, including atmospheric tides and core–mantle friction.  International Journal of Astrobiology 14, 233–254.
  * Dai and Wang (1999) Dai, A., Wang, J., 1999.  Diurnal and semidiurnal tides in global surface pressure fields.  Journal of the atmospheric sciences 56, 3874–3891.
  * Dobrovolskis and Ingersoll (1980) Dobrovolskis, A.R., Ingersoll, A.P., 1980\.  Atmospheric tides and the rotation of venus i. tidal theory and the balance of torques.  Icarus 41, 1–17.
  * Farhat et al. (2022) Farhat, M., Auclair-Desrotour, P., Boué, G., Laskar, J., 2022\.  The resonant tidal evolution of the earth-moon distance.  Astronomy & Astrophysics 665, L1.
  * Garratt (1994) Garratt, J.R., 1994.  The atmospheric boundary layer.  Earth-Science Reviews 37, 89–134.
  * Gerkema and Zimmerman (2008) Gerkema, T., Zimmerman, J., 2008\.  An introduction to internal waves.  Lecture Notes, Royal NIOZ, Texel 207\.
  * Gold and Soter (1969) Gold, T., Soter, S., 1969.  Atmospheric tides and the resonant rotation of venus.  Icarus 11, 356–366.
  * Goldreich and Soter (1966) Goldreich, P., Soter, S., 1966\.  Q in the solar system.  icarus 5, 375–389.
  * Gough (1981) Gough, D., 1981.  Solar interior structure and luminosity variations, in: Physics of Solar Variations: Proceedings of the 14th ESLAB Symposium held in Scheveningen, The Netherlands, 16–19 September, 1980, Springer. pp. 21–34.
  * Gu et al. (2019) Gu, P.G., Peng, D.K., Yen, C.C., 2019.  Modeling the thermal bulge of a hot jupiter with the two-stream approximation.  The Astrophysical Journal 887, 228\.
  * Hagan and Forbes (2002) Hagan, M., Forbes, J., 2002\.  Migrating and nonmigrating diurnal tides in the middle and upper atmosphere excited by tropospheric latent heat release.  Journal of Geophysical Research: Atmospheres 107, ACL–6.
  * Haurwitz and Cowley (1973) Haurwitz, B., Cowley, A.D., 1973\.  The diurnal and semidiurnal barometric oscillations global distribution and annual variation.  pure and applied geophysics 102, 193–222.
  * Hilsenrath (1955) 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. volume 564\.  US Department of Commerce, National Bureau of Standards.
  * Holmberg (1952) Holmberg, E., 1952.  A suggested explanation of the present value of the velocity of rotation of the earth.  Geophysical Journal International 6, 325–330.
  * Holtslag and Boville (1993) Holtslag, A., Boville, B., 1993\.  Local versus nonlocal boundary-layer diffusion in a global climate model.  Journal of climate 6, 1825–1842.
  * Hough (1898) Hough, S.S., 1898.  On the Application of Harmonic Analysis to the Dynamical Theory of the Tides. Part II: On the General Integration of Laplace’s Dynamical Equations.  Philosophical Transactions of the Royal Society of London Series A 191, 139–185.
  * Hourdin et al. (2006) Hourdin, F., Musat, I., Bony, S., Braconnot, P., Codron, F., Dufresne, J.L., Fairhead, L., Filiberti, M.A., Friedlingstein, P., Grandpeix, J.Y., et al., 2006.  The lmdz4 general circulation model: climate performance and sensitivity to parametrized physics with emphasis on tropical convection.  Climate Dynamics 27, 787–813.
  * Hut (1981) Hut, P., 1981.  Tidal evolution in close binary systems.  Astronomy and Astrophysics 99, 126–140.
  * Ingersoll and Dobrovolskis (1978) Ingersoll, A.P., Dobrovolskis, A.R., 1978\.  Venus’ rotation and atmospheric tides.  Nature 275, 37–38.
  * Kelvin (1882) Kelvin, L.o.k.a.W.T., 1882.  On the thermodynamic acceleration of the earth’s rotation.  Proceedings of the Royal Society of Edinburgh 11, 396–405.
  * Klatt et al. (2021) Klatt, J.M., Chennu, A., Arbic, B.K., Biddanda, B., Dick, G.J., 2021.  Possible link between earth’s rotation rate and oxygenation.  Nature Geoscience 14, 564–570.
  * Knauth (2005) Knauth, L.P., 2005.  Temperature and salinity history of the precambrian ocean: implications for the course of microbial evolution, in: Geobiology: Objectives, concepts, perspectives. Elsevier, pp. 53–69.
  * Krissansen-Totton et al. (2018) Krissansen-Totton, J., Arney, G.N., Catling, D.C., 2018.  Constraining the climate and ocean ph of the early earth with a geological carbon cycle model.  Proceedings of the National Academy of Sciences 115, 4105–4110.
  * Lacis and Oinas (1991) Lacis, A.A., Oinas, V., 1991\.  A description of the correlated k distribution method for modeling nongray gaseous absorption, thermal emission, and multiple scattering in vertically inhomogeneous atmospheres.  Journal of Geophysical Research: Atmospheres 96, 9027–9063.
  * Lambeck (1980) Lambeck, K., 1980.  The Earth’s variable rotation: geophysical causes and consequences.  Cambridge University Press.
  * Lantink et al. (2022) Lantink, M.L., Davies, J.H., Ovtcharova, M., Hilgen, F.J., 2022\.  Milankovitch cycles in banded iron formations constrain the earth–moon system 2.46 billion years ago.  Proceedings of the National Academy of Sciences 119, e2117146119.
  * Laplace (1798) Laplace, P., 1798.  Traite de mecanique celeste, five volumes.  Paris: Chez JBM Duprat .
  * Laskar and Correia (2004) Laskar, J., Correia, A.C., 2004\.  The rotation of extra-solar planets, in: Extrasolar Planets: Today and Tomorrow, p. 401\.
  * Laskar et al. (2023) Laskar, J., Farhat, M., Lantink, M., Auclair-Desrotour, P., Boué, G., Sinnesael, M., 2023\.  Did atmospheric thermal tides cause a daylength locking in the Precambrian? A review on recent results. arXiv:2309.11479, doi:10.48550/arXiv.2309.11479.
  * Leconte et al. (2015) Leconte, J., Wu, H., Menou, K., Murray, N., 2015\.  Asynchronous rotation of earth-mass planets in the habitable zone of lower-mass stars.  Science 347, 632–635.
  * Lee (2020) Lee, U., 2020.  Tidal oscillations of rotating hot jupiters.  Monthly Notices of the Royal Astronomical Society 494, 3141–3155.
  * Lee and Saio (1997) Lee, U., Saio, H., 1997.  Low-frequency nonradial oscillations in rotating stars. i. angular dependence.  The Astrophysical Journal 491, 839\.
  * Lindzen (1978) Lindzen, R.S., 1978.  Effect of daily variations of cumulonimbus activity on the atmospheric semidiurnal tide.  Monthly Weather Review 106, 526–533.
  * Lindzen and Blake (1972) Lindzen, R.S., Blake, D., 1972\.  Lamb waves in the presence of realistic distributions of temperature and dissipation.  Journal of Geophysical Research 77, 2166–2176.
  * Lindzen and Chapman (1968) Lindzen, R.S., Chapman, S., 1968\.  The application of classical atmospheric tidal theory.  Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 303, 299–316.
  * Lindzen and McKenzie (1967) Lindzen, R.S., McKenzie, D.J., 1967\.  Tidal theory with newtonian cooling.  pure and applied geophysics 66, 90–96.
  * Madsen (1977) Madsen, O.S., 1977.  A realistic model of the wind-induced ekman boundary layer.  Journal of Physical Oceanography 7, 248–255.
  * Mignard (1980) Mignard, F., 1980.  The evolution of the lunar orbit revisited, ii.  The Moon and the planets 23, 185–201.
  * Mitchell and Kirscher (2023) Mitchell, R., Kirscher, U., 2023\.  Mid-Proterozoic day length stalled by tidal resonance.  Nature Geoscience 9, 1–3.
  * Nieuwstadt (1983) Nieuwstadt, F., 1983.  On the solution of the stationary, baroclinic ekman-layer equations with a finite boundary-layer height.  Boundary-layer meteorology 26, 377–390.
  * Ogilvie (2014) Ogilvie, G.I., 2014.  Tidal dissipation in stars and giant planets.  Annual Review of Astronomy and Astrophysics 52, 171–210.
  * Pannella (1972a) Pannella, G., 1972.  Paleontological evidence on the Earth’s rotational history since early Precambrian.  Astrophysics and Space Science 16, 212–237.
  * Pannella (1972b) Pannella, G., 1972.  Precambrian stromatolites as paleontological clocks.  Int. Geol. Congr. Rep. Sess., 24th 50–57.
  * Pierrehumbert (2010) Pierrehumbert, R.T., 2010.  Principles of planetary climate.  Cambridge University Press.
  * Robert and Chaussidon (2006) Robert, F., Chaussidon, M., 2006\.  A palaeotemperature curve for the precambrian oceans based on silicon isotopes in cherts.  Nature 443, 969–972.
  * Sakazaki et al. (2017) Sakazaki, T., Hamilton, K., Zhang, C., Wang, Y., 2017\.  Is there a stratospheric pacemaker controlling the daily cycle of tropical rainfall?  Geophysical Research Letters 44, 1998–2006.
  * Sakazaki and Hamilton (2020) Sakazaki, T., Hamilton, K., 2020\.  An array of ringing global free modes discovered in tropical surface pressure data  Journal of the Atmospheric Sciences 77, 2519–2539.
  * Schindelegger and Ray (2014) Schindelegger, M., Ray, R.D., 2014\.  Surface pressure tide climatologies deduced from a quality-controlled network of barometric observations.  Monthly Weather Review 142, 4872–4889.
  * Scrutton (1978) Scrutton, C., 1978.  Periodic growth features in fossil organisms and the length of the day and month, in: Tidal friction and the earth’s rotation. Springer, pp. 154–196.
  * Showman and Guillot (2002) Showman, A.P., Guillot, T., 2002\.  Atmospheric circulation and tides of “51 pegasus b-like” planets.  Astronomy & Astrophysics 385, 166–180.
  * Siebert (1961) Siebert, M., 1961.  Atmospheric tides, in: Advances in geophysics. Elsevier. volume 7, pp. 105–187.
  * Sleep and Zahnle (2001) Sleep, N.H., Zahnle, K., 2001\.  Carbon dioxide cycling and implications for climate on ancient earth.  Journal of Geophysical Research: Planets 106, 1373–1399.
  * Tiesinga et al. (2021) Tiesinga, E., Mohr, P.J., Newell, D.B., Taylor, B.N., 2021\.  Codata recommended values of the fundamental physical constants: 2018.  Journal of Physical and Chemical Reference Data 50, 033105.
  * Unno et al. (1989) Unno, W., Osaki, Y., Ando, H., Saio, H., Shibahashi, H., 1989\.  Nonradial oscillations of stars, ed.  Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H 755.
  * Vallis (2017) Vallis, G.K., 2017.  Atmospheric and oceanic fluid dynamics.  Cambridge University Press.
  * Van Wijk and De Vries (1963) Van Wijk, W.R., De Vries, D., 1963\.  Periodic temperature variations in a homogeneous soil.  Physics of plant environment 1, 103–143.
  * Vichare and Rajaram (2013) Vichare, G., Rajaram, R., 2013\.  Diurnal and semi-diurnal tidal structures due to o2, o3 and h2o heating.  Journal of Earth System Science 122, 1207–1217.
  * Walker and Zahnle (1986) Walker, J.C., Zahnle, K.J., 1986\.  Lunar nodal tide and distance to the moon during the precambrian.  Nature 320, 600–602.
  * Wang et al. (2016) Wang, H., Boyd, J.P., Akmaev, R.A., 2016.  On computation of hough functions.  Geoscientific Model Development 9, 1477–1488.
  * Wilkes (1949) Wilkes, M.V., 1949.  Oscillations of the earth’s atmosphere.
  * Williams (1990) Williams, G.E., 1990.  Tidal rhythmites: key to the history of the earth’s rotation and the lunar orbit.  Journal of Physics of the Earth 38, 475–491.
  * Williams (2000) Williams, G.E., 2000.  Geological constraints on the precambrian history of earth’s rotation and the moon’s orbit.  Reviews of Geophysics 38, 37–59.
  * de Wit and Furnes (2016) de Wit, M.J., Furnes, H., 2016\.  3.5-ga hydrothermal fields and diamictites in the barberton greenstone belt—paleoarchean crust in cold environments.  Science advances 2, e1500368.
  * Wu et al. (2023) Wu, H., Murray, N., Menou, K., Lee, C., Leconte, J., 2023\.  Why the day is 24 hours long: The history of earth’s atmospheric thermal tide, composition, and mean temperature.  Science Advances 9, eadd2499.
  * Zahnle and Walker (1987) Zahnle, K., Walker, J.C., 1987\.  A constant daylength during the precambrian era?  Precambrian Research 37, 95–105.

##  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 and Lindzen,
1970), and we augment the heat transfer equation by the additional radiative
term $J_{\rm rad}$ described in the main text. These equations read:

$\displaystyle\partial_{{t}}{V}_{\theta}-2\Omega\cos\theta{V}_{\varphi}=-\frac{1}{R_{\rm
p}}\partial_{\theta}\left(\frac{\delta p}{\rho_{0}}-{U}\right)\,,$ (25)
$\displaystyle\partial_{{t}}{V}_{\varphi}\\!\\!+2\Omega\cos\theta{V}_{\theta}=\\!-\frac{1}{R_{\rm
p}\sin\theta}\partial_{\varphi}\\!\left(\\!\frac{\delta
p}{\rho_{0}}\\!-\\!{U}\\!\right)\\!,$ (26) $\displaystyle
0=-\partial_{z}\delta p-g\delta\rho+\rho_{0}\partial_{z}U\,,$ (27)
$\displaystyle\frac{D\rho}{Dt}=-\rho_{0}\chi\,,$ (28)
$\displaystyle\frac{1}{\Gamma_{1}}\frac{Dp}{Dt}-gH\frac{D\rho}{Dt}=\kappa\rho_{0}J-\rho_{0}\mathcal{R}_{\rm
s}\sigma_{0}\delta T\,,$ (29) $\displaystyle\chi-\partial_{z}V_{r}=\nabla_{\rm
h}\cdot\boldsymbol{V}\,,$ (30) $\displaystyle\frac{\delta
p}{p_{0}}=\frac{\delta\rho}{\rho_{0}}+\frac{\delta T}{T_{0}}.$ (31)

This system describes atmospheric momentum conservation (Eqs. 25 \- 27), mass
conservation (Eq. 28), heat transfer with Newtonian cooling (Eq. 29), and the
ideal gas law (Eq. 31). Furthermore, we have adopted the traditional
approximation (e.g., Unno et al., 1989), 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 (27), the traditional approximation amounts to imposing the
hydrostatic equilibrium for the tidal perturbation in addition to the
background profiles, which are in turn averaged over the sphere. The
hydrostatic approximation would allow us later to define the tidal torque as a
function of the pressure anomaly at the surface.

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

$\frac{Dy}{Dt}=\partial_{t}y+V_{r}\frac{dy_{0}}{dz},$ (32)

and we denote by $\chi$ the divergence of the velocity vector field
$\boldsymbol{V}$, and by the notation $\nabla_{\rm h}\cdot\boldsymbol{V}$ its
horizontal divergence such that

$\nabla_{\rm h}\cdot\boldsymbol{V}=\frac{1}{R_{\rm
p}\sin\theta}\left[\partial_{\theta}(\sin\theta
V_{\theta})+\partial_{\varphi}V_{\varphi}\right].$ (33)

We also introduce the calculation variable

$G=-\frac{1}{\Gamma_{1}p_{0}}\frac{Dp}{Dt}.$ (34)

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\delta\rho$
term in Eq.(27) is supposed to vanish since it is often referred to as the
buoyancy term. However, this is not exactly the case since internal gravity
waves are not the only source of density fluctuations. In fact, density
fluctuations can also result from the planetary-scale compressibility waves,
read Lamb modes. As such, the atmosphere can be at the same time neutrally
stratified and still feature strong density variations across the horizontal
direction. With this subtlety clarified, the tidal perturbations are then
expanded in Fourier series of time and longitude, such that the tidal
excitation frequency is denoted by $\sigma$. We introduce the reference
velocity $v_{0}$, and we non-dimensionalize all the variables; namely:

$\displaystyle dz=Hdx,\hskip 11.38092ptt=\sigma^{-1}\tilde{t},\hskip
11.38092pt\boldsymbol{V}=v_{0}\tilde{\boldsymbol{V}},\hskip 11.38092pt\delta
p=p_{0}\tilde{\delta p},$
$\displaystyle\delta\rho=\rho_{0}\tilde{\delta\rho},\hskip 8.5359pt\delta
T=T_{0}\tilde{\delta T},\hskip 8.5359ptU=v_{0}^{2}\tilde{U},\hskip
8.5359ptG=\frac{\sigma}{\Gamma_{1}}\tilde{G},$ $\displaystyle\nabla_{\rm
h}\cdot=R_{\rm p}^{-1}\tilde{\nabla}\cdot,\hskip 11.38092pt\nabla_{\rm
h}=R_{\rm p}^{-1}\tilde{\nabla},\hskip 11.38092ptJ=\frac{\sigma
gH}{\kappa}\tilde{J}\,,$ (35)

Under these definitions, the material derivative now reads

$\frac{Dy}{Dt}=\sigma y_{0}\left[\partial_{\tilde{t}}\tilde{\delta
y}+\frac{v_{0}}{\sigma H}\frac{d\ln y_{0}}{dx}\tilde{V}_{r}\right]\,,$ (36)

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

$\partial_{\tilde{t}}\tilde{\delta p}+\frac{v_{0}}{\sigma H}\frac{d\ln
p_{0}}{dx}\tilde{V}_{r}=-\tilde{G}.$ (37)

We next set $v_{0}=\sigma R_{\rm p}$ and introduce the ratios of length scale
$\eta=R_{\rm p}/H$, and frequency $\alpha=\sigma/\sigma_{0}$. Allowing for
these changes of variables in the governing system of equations, one directly
obtains a dimensionless system in the form:

$\displaystyle\partial_{\tilde{t}}\tilde{V}_{\theta}-\nu\cos\theta\tilde{V}_{\varphi}=-\partial_{\theta}\left(\frac{gH}{v_{0}^{2}}\tilde{\delta
p}-\tilde{U}\right)\,,$ (38)
$\displaystyle\partial_{\tilde{t}}\tilde{V}_{\varphi}+\nu\cos\theta\tilde{V}_{\theta}=-\frac{1}{\sin\theta}\partial_{\varphi}\left(\frac{gH}{v_{0}^{2}}\tilde{\delta
p}-\tilde{U}\right)\,,$ (39) $\displaystyle(\partial_{x}-1)\tilde{\delta
p}=-\tilde{\delta\rho}+\frac{v_{0}^{2}}{gH}\partial_{x}\tilde{U}\,,$ (40)
$\displaystyle\partial_{\tilde{t}}\tilde{\delta\rho}+\tilde{\nabla}\cdot\tilde{\boldsymbol{V}}=-\eta\left(\partial_{x}+\frac{d\ln{\rho_{0}}}{dx}\right)\tilde{V}_{r}\,,$
(41)
$\displaystyle\frac{\tilde{G}}{\Gamma_{1}}+\partial_{\tilde{t}}\tilde{\delta\rho}+\eta\frac{d\ln\rho_{0}}{dx}\tilde{V}_{r}=-\tilde{J}+\frac{\delta\tilde{T}}{\alpha}\,,$
(42) $\displaystyle\tilde{\delta p}=\tilde{\delta\rho}+\tilde{\delta T}.$ (43)

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

$f(x,\theta,\varphi,\tilde{t})=\Re\left\\{\sum_{m}f^{m,\nu}(x,\theta)e^{i(\tilde{t}+m\varphi)}\right\\},$
(44)

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. (44)
into series of Hough functions (Hough, 1898) describing the horizontal tidal
flow; namely:

$f^{m,\nu}(x,\theta)=\sum_{n}f_{n}^{m,\nu}(x)\Theta_{n}^{m,\nu}(\theta).$ (45)

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

$\mathcal{L}\Theta_{n}=-\Lambda_{n}\Theta_{n},$ (46)

where the operator $\mathcal{L}$ is defined as

$\displaystyle\mathcal{L}$
$\displaystyle=\frac{1}{\sin\theta}\frac{d}{d\theta}\left(\frac{\sin\theta}{1-{\nu}^{2}\cos^{2}\theta}\frac{d}{d\theta}\right)$
$\displaystyle-\frac{1}{1-{\nu}^{2}\cos^{2}\theta}\left(m{\nu}\frac{1+{\nu}^{2}\cos^{2}\theta}{1-{\nu}^{2}\cos^{2}\theta}+\frac{m^{2}}{\sin^{2}\theta}\right).$
(47)

##  Appendix B Retrieving the wave equation of vertical dynamics

The governing system of dimensionless equations (38-43) decouples into
equations describing the horizontal structure (38-39), and those describing
the vertical structure. In the Fourier domain, the latter equations read:

$\displaystyle i\tilde{\delta p}_{n}-\eta\tilde{V}_{r;n}=-\tilde{G}_{n},$ (48)
$\displaystyle\\!\\!\left(\partial_{x}-1\right)\tilde{\delta
p}_{n}=-\tilde{\delta\rho}_{n}+\frac{v_{0}^{2}}{gH}\partial_{x}\tilde{U}_{n},$
(49)
$\displaystyle\tilde{\delta\rho}_{n}\\!-\\!i\eta\left(\\!\partial_{x}\\!+\\!\frac{d\ln\rho_{0}}{dx}\\!\right)\\!\tilde{V}_{r;n}\\!=\\!-\mathcal{L}\\!\left(\frac{gH}{v_{0}^{2}}\tilde{\delta
p}_{n}-\tilde{U}_{n}\\!\right)\\!,$ (50)
$\displaystyle\frac{\tilde{G}_{n}}{\Gamma_{1}}+i\tilde{\delta\rho}_{n}+\eta\frac{d\ln\rho_{0}}{dx}\tilde{V}_{r;n}=-\tilde{J}_{n}+\frac{\tilde{\delta
T}_{n}}{\alpha},$ (51) $\displaystyle\tilde{\delta
p}_{n}=\tilde{\delta\rho}_{n}+\tilde{\delta T}_{n}.$ (52)

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

$\tilde{\delta p}_{n}\\!\left(\\!1\\!+\\!\frac{d\ln
H}{dx\\!}-\\!\frac{i}{\alpha}\right)-\tilde{\delta\rho}_{n}\left(\\!1\\!-\\!\frac{i}{\alpha}\right)\\!=\frac{iN^{2}_{\rm
B}H}{g}\tilde{G}_{n}-i\tilde{J}_{n}.$ (53)

Next, we use the hydrostatic equilibrium condition (49) to express
$\tilde{\delta\rho}_{n}$ in terms of $\tilde{\delta p}_{n}$ and
$\tilde{U}_{n}$ in the above equation to obtain

$\displaystyle\partial_{x}\tilde{\delta p}_{n}$ $\displaystyle=i\frac{\alpha
HN_{\rm B}^{2}}{g(\alpha-i)}\tilde{G}_{n}-\frac{\alpha}{\alpha-i}\frac{d\ln
H}{dx}\tilde{\delta p}_{n}-\frac{i\alpha}{\alpha-i}\tilde{J}_{n}$
$\displaystyle+\frac{v_{0}^{2}}{gH}\partial_{x}\tilde{U}_{n}.$ (54)

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

$\partial_{x}\tilde{G}_{n}=\tilde{G}_{n}-i\frac{gH}{v_{0}^{2}}\mathcal{L}\tilde{\delta
p}_{n}+i\left(\mathcal{L}-\frac{v_{0}^{2}}{gH}\partial_{x}\right)\tilde{U}_{n}.$
(55)

We thus have two first order partial differential equations (54 \- 55) of the
form

$\displaystyle\partial_{x}\tilde{G}_{n}=A_{1}\tilde{G}_{n}+B_{1}\tilde{\delta
p}_{n}+C_{1},$ (56) $\displaystyle\partial_{x}\tilde{\delta
p}_{n}=A_{2}\tilde{G}_{n}+B_{2}\tilde{\delta p}_{n}+C_{2}.$ (57)

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

$\left[\partial_{x}^{2}+A\partial_{x}+B\right]\tilde{G}_{n}=C_{n},$ (58)

where

$\displaystyle A(x)=-1+\frac{i}{\alpha-i}\frac{d\ln H}{dx},$ (59)
$\displaystyle B(x)=-\frac{\alpha}{\alpha-i}\left[\frac{i}{\alpha}\frac{d\ln
H}{dx}+\frac{H^{2}N_{\rm B}^{2}}{v_{0}^{2}}\mathcal{L}\right],$ (60)
$\displaystyle
C_{n}(x)=-\frac{gH\alpha}{v_{0}^{2}(\alpha-i)}\mathcal{L}\tilde{J}_{n}-\Bigg{[}i\frac{d}{dx}\left(\frac{v_{0}^{2}}{gH}\frac{d}{dx}\right)$
(61) $\displaystyle\hskip 22.76228pt+\frac{1}{\alpha-i}\frac{d\ln
H}{dx}\left(\mathcal{L}-\frac{v_{0}^{2}}{gH}\frac{d}{dx}\right)\Bigg{]}\tilde{U}_{n}.$
(62)

Next, we implement a change of variable of the form
$\tilde{G}_{n}(x)=\Phi(x)\Psi_{n}(x)$ to write Eq. (58) in the standard form
of wave equations. We find that the first order term in Eq. (58) cancels when

$2\frac{d\Phi(x)}{dx}+A(x)\Phi(x)=0,$ (63)

which is satisfied by

$\displaystyle\Phi(x)$
$\displaystyle=\exp\left(-\frac{1}{2}\int_{0}^{x}A(x^{\prime})dx^{\prime}\right).$
(64)

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

$\displaystyle\hat{k}_{n}^{2}$
$\displaystyle=-\frac{1}{4}\left(A^{2}+2\frac{dA}{dx}-4B\right)$
$\displaystyle=-\frac{1}{4}\Bigg{\\{}\\!\\!\\!\left(\\!1\\!-\\!\frac{i}{\alpha-i}\frac{d\ln
H}{dx}\right)^{2}\\!\\!\\!\\!+2\frac{d}{dx}\\!\\!\left(\frac{i}{\alpha-i}\frac{d\ln
H}{dx}\right)$ $\displaystyle-\frac{4\alpha}{\alpha-i}\\!\left[\frac{gHN_{\rm
B}^{2}\Lambda_{n}}{v_{0}^{2}}-\\!\frac{i}{\alpha}\frac{d\ln
H}{dx}\right]\Bigg{\\}}.$ (65)

With the definitions of the dimensionless control parameters given in Table 1,
the wave number $\hat{k}_{n}$ (Eq.65) and the forcing term (Eq.62) 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 $\Psi(x)_{n}$,
or equivalently the variable $\tilde{G}_{n}$.

First, the expression of pressure variations is readily obtained from Eq.(55),
which for a given Hough mode and using the definition of $\beta$ (Table 1)
reads

$\tilde{\delta
p}_{n}=\frac{1}{i\beta\Lambda_{n}}\left(\frac{d\tilde{G}_{n}}{dx}-\tilde{G}_{n}\right)+\frac{1}{\beta}\left(1+\frac{1}{\beta\Lambda_{n}}\frac{d}{dx}\right)\tilde{U}_{n}.$
(66)

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

$\displaystyle\tilde{V}_{\theta;n}$
$\displaystyle=\frac{\partial_{\theta}+m\nu\cot\theta}{1-\nu^{2}\cos^{2}\theta}\left(\frac{gH}{v_{0}^{2}}\tilde{\delta
p}_{n}-\tilde{U}_{n}\right),$ (67) $\displaystyle\tilde{V}_{\varphi;n}$
$\displaystyle=-\frac{\nu\cos\theta\partial_{\theta}+{m}{\sin^{-1}\theta}}{1-\nu^{2}\cos^{2}\theta}\left(\frac{gH}{v_{0}^{2}}\tilde{\delta
p}_{n}-\tilde{U}_{n}\right),$ (68)

which gives

$\displaystyle\tilde{V}_{\theta;n}$
$\displaystyle=\frac{1}{\Lambda_{n}}\left(\frac{d\tilde{G}_{n}}{dx}-\tilde{G}_{n}+i\beta^{-1}\frac{d\tilde{U}_{n}}{dx}\right),$
(69) $\displaystyle\tilde{V}_{\varphi;n}$
$\displaystyle=\frac{i}{\Lambda_{n}}\left(\frac{d\tilde{G}_{n}}{dx}-\tilde{G}_{n}+i\beta^{-1}\frac{d\tilde{U}_{n}}{dx}\right).$
(70)

The vertical component of the velocity is established from equation (48),
replacing $\tilde{\delta p}_{n}$ by its solution (66):

$\displaystyle\tilde{V}_{r;n}$
$\displaystyle=\frac{1}{\eta\beta\Lambda_{n}}\left[\frac{d\tilde{G}_{n}}{dx}+\left(\beta\Lambda_{n}-1\right)\tilde{G}_{n}\right]$
$\displaystyle+\frac{i}{\eta\beta}\left(1+\frac{1}{\beta\Lambda_{n}}\frac{d}{dx}\right)\tilde{U}_{n}.$
(71)

Finally, the solution for the density perturbation is obtained by replacing
the solution of $\tilde{\delta p}_{n}$ (66) in equation (53) to get:

$\displaystyle\tilde{\delta\rho}_{n}$
$\displaystyle=\frac{1}{i\beta\Lambda_{n}}\left(1+\frac{\kappa\alpha(\gamma-1)}{\alpha-i}\right)\Bigg{[}\frac{d\tilde{G}_{n}}{dx}+$
$\displaystyle\left(\frac{\alpha\beta\gamma\kappa{\Lambda}}{\alpha-i+\kappa\alpha(\gamma-1)}-1\right)\tilde{G}_{n}\Bigg{]}+\frac{i\alpha\tilde{J}_{n}}{\alpha-i}$
$\displaystyle+\frac{1}{\beta}\left(1+\frac{\kappa\alpha(\gamma-1)}{\alpha-i}\right)\left(1+\frac{1}{\beta\Lambda_{n}}\frac{d}{dx}\right)\tilde{U}_{n},$
(72)

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

##  Appendix D Thermal exchange mechanisms and the total propagating flux

Figure 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. (73).

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, $\delta F_{\rm inc}$
generates a variation $\delta T_{\rm s}$ in the surface temperature $T_{\rm
s}$. We denote the proportionality by the complex transfer function
$\mathcal{B}_{\rm gr}^{\sigma}$ such that $\delta T_{\rm s}=\mathcal{B}_{\rm
gr}^{\sigma}\delta F_{\rm inc}$. The variation $\delta T_{\rm s}$ induces
radiative emission with magnitude $\delta F_{\rm rad}$. A fraction of the
incident power, $\delta Q_{\rm gr}$, is then transmitted to the ground by
thermal conduction, and another fraction, $\delta Q_{\rm atm}$, is transmitted
to the atmosphere through turbulent thermal diffusion (see Figure 7).

Having isolated these thermal mechanisms, the total power balance of the
thermal perturbation at the surface interface is expressed as999we ignore
latent heat exchanges associated with changes of state.:

$\delta F_{\rm inc}-\delta F_{\rm rad}=\delta Q_{\rm gr}+\delta Q_{\rm atm}.$
(73)

To define each of these terms explicitly, we start with the surface radiative
emission term $\delta F_{\rm rad}$, which can be obtained by differentiating
the Stefan-Boltzmann law as a function of the surface temperature assuming
that the planetary surface radiates as a black body. Namely we have

$\delta F_{\rm rad}=4\sigma_{\rm SB}T_{\rm s}^{3}\delta T_{\rm s},$ (74)

with $\sigma_{\rm SB}=5.670373\times 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,
$\delta Q_{\rm gr}$ and $\delta Q_{\rm atm}$, are derived using the gradient
flux theory (Garratt, 1994):

$\delta Q_{\rm gr}(t)=\frac{k_{\rm gr}}{H}\frac{\partial\delta T}{\partial
x}\Bigg{|}_{x=0^{-}},$ (75) $\delta Q_{\rm atm}(t)=-\frac{k_{\rm
atm}}{H}\frac{\partial\delta T}{\partial x}\Bigg{|}_{x=0^{+}}.$ (76)

Here, $k_{\rm gr}$ is the thermal conductivity of the ground, and $k_{\rm
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)

$i\sigma\delta T=\frac{K_{\rm gr}}{H^{2}}\frac{\partial^{2}\delta T}{\partial
x^{2}}\,\,\,\,\,{\rm for}\,\,x<0\,,$ (77) $i\sigma\delta T=\frac{K_{\rm
atm}}{H^{2}}\frac{\partial^{2}\delta T}{\partial x^{2}}\,\,\,\,\,{\rm
for}\,\,x>0.$ (78)

with $K_{\rm gr}=k_{\rm gr}/(\rho_{0}(0^{-})C_{\rm gr})$ and $K_{\rm
atm}=k_{\rm atm}/(\rho_{0}(0^{+})C_{\rm p})$ being the thermal diffusivities
near the surface in the ground and the atmosphere respectively, $C_{\rm gr}$
and $C_{\rm p}$ being the respective thermal capacities per unit mass, while
$\rho_{0}(0^{-})$ and $\rho_{0}(0^{+})$ are the densities in the vicinity of
the surface interface. These equations have solutions of the form:

$\delta T(x)=\delta T_{\rm s}\exp\left\\{(1+is)x/h_{\rm
gr}^{\sigma}\right\\},\,\,\,\,{\rm for}\,\,x\leq 0,$ (79) $\delta T(x)=\delta
T_{\rm s}\exp\left\\{-(1+is)x/h_{\rm atm}^{\sigma}\right\\},\,\,\,\,{\rm
for}\,\,x>0.$ (80)

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

$h_{\rm gr}^{\sigma}=\frac{1}{H}\frac{2K_{\rm gr}}{|\sigma|}\hskip
17.07182pt{\rm and}\hskip 17.07182pth_{\rm
atm}^{\sigma}=\frac{1}{H}\frac{2K_{\rm atm}}{|\sigma|}.$ (81)

To obtain the explicit form of $\mathcal{B}_{\rm gr}^{\sigma}$, we first
rewrite the power budget balance of Eq.(73) as

$\delta F_{\rm inc}-\frac{k_{\rm gr}}{H}\frac{\partial\delta T}{\partial
x}\Bigg{|}_{x=0^{-}}\\!\\!\\!\\!\\!\\!\\!\\!\\!=-\frac{k_{\rm
atm}}{H}\frac{\partial\delta T}{\partial
x}\Bigg{|}_{x=0^{+}}\\!\\!\\!\\!\\!\\!\\!\\!\\!+4\sigma_{\rm SB}T_{\rm
s}^{3}\delta T_{\rm s}\,,$ (82)

and we substitute the solutions of $\delta T(x)$ of Eq. (79) to obtain

$\mathcal{B}_{\rm gr}^{\sigma}=\left\\{\left(4\sigma_{\rm SB}T_{\rm
s}^{3}\right)\left[1+(1+si)\sqrt{\tau_{\rm bl}|\sigma|}\right]\right\\}^{-1}.$
(83)

Here we have defined the boundary layer timescale $\tau_{\rm bl}$ as

$\tau_{\rm bl}=\frac{1}{2}\left(\frac{I_{\rm gr}+I_{\rm atm}}{4\sigma_{\rm
SB}T_{\rm s}^{3}}\right)^{2}\,,$ (84)

where the functions $I_{\rm gr}$ and $I_{\rm atm}$ denote, respectively, the
thermal inertias of the ground and the atmospheric surface layers, and are
specifically given by:

$\displaystyle I_{\rm gr}=$ $\displaystyle\rho_{0}(0^{-})C_{\rm
gr}\sqrt{K_{\rm gr}}\,,$ $\displaystyle I_{\rm atm}=$
$\displaystyle\rho_{0}(0^{+})C_{\rm p}\sqrt{K_{\rm atm}}\,.\,$ (85)

Finally, with $\delta Q_{\rm gr}$ written as

$\delta Q_{\rm gr}=\frac{k_{\rm gr}}{\sqrt{2K_{\rm
gr}}}\sqrt{|\sigma|}\mathcal{B}_{\rm gr}^{\sigma}(1+si)\delta F_{\rm inc},$
(86)

we are fully geared to define the total propagating flux feeding the
atmosphere, $\delta F_{\rm tot}$, as

$\delta F_{\rm tot}=\delta F_{\rm inc}-\delta Q_{\rm gr},$ (87)

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_{\rm p}=1005$ J/kg$\cdot$K
(Hilsenrath, 1955). The volumetric heat capacity in the ground
$\rho_{0}(0^{-})C_{\rm 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\times 10^{6}$ J m3 K-1 (Van Wijk and De Vries, 1963). The value
for the density $\rho_{0}$ is dependent on the planetary setting understudy
and is computed via $\rho_{0}=p_{0}/gH$. The values for ground thermal
conductivities and diffusivities vary among rock types and water content
taking typical values of $k_{\rm gr}\sim 0.5$ – $4$ W m-1 K-1 and $K_{\rm gr}$
$\sim$ $10^{-6}$ – $10^{-8}$ m2 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_{\rm 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$ m2 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
$\overline{K}_{\rm atm}=3.6$ m2 s-1. This variability of the ground and
atmospheric thermal parameters propagates to our model free parameter
$\tau_{\rm bl}$, or equivalently the frequency $\sigma_{\rm bl}$, via the
thermal ground and atmospheric thermal inertias defined in Eq. (85). Thus in
sweeping the limits of $\sigma_{\rm 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 $\sigma_{\rm
bl}=10^{-5}$ s-1, for which the accelerating torque at the Lamb resonance
completely vanishes, corresponds to $k_{\rm gr}=1.25$ W m-1 K-1, $K_{\rm
gr}=5\times 10^{-7}$ m2 s-1, and $K_{\rm atm}=\overline{K}_{\rm atm}=3.6$ m2
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 $\delta
F_{\rm inc}$ (denoted hereafter by $F$ for convenience) as:

$F=\begin{cases}F_{\star}\cos\Theta,&\text{for}\ 0\leq\Theta\leq\pi/2\\\
0,&\text{otherwise}\,,\end{cases}$ (88)

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

$F(\Theta)=\sum_{l=0}^{\infty}F_{l}P_{l}(\cos\Theta),$ (89)

where the expansion coefficients are given by

$\displaystyle F_{l}$
$\displaystyle=\frac{2l+1}{2}\int_{0}^{\pi}F(\Theta)P_{l}(\cos\Theta)\sin\Theta
d\Theta\,,$
$\displaystyle=\frac{2l+1}{2}F_{\star}\int_{0}^{\pi/2}P_{l}(\cos\Theta)\sin\Theta\cos\Theta
d\Theta.$ (90)

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

$P_{l}(\cos\Theta)=\frac{4\pi}{2l+1}\sum_{m=-l}^{m=l}Y_{lm}(\theta,\varphi)Y_{lm}^{*}(\theta_{\rm
S},\varphi_{\rm S})\,,$ (91)

where the asterisk corresponds to complex conjugation, and the coordinates
$(\theta_{\rm S},\varphi_{\rm S})=\left[\pi/2,(\Omega-n_{\star})t\right]$
define the position of the star when ignoring planetary obliquity. We thus
have:

$F=\sum_{l=0}^{\infty}\sum_{m=-l}^{m=l}C_{l}F_{l}Y_{lm}(\theta,\varphi)Y_{lm}^{*}(\theta_{\rm
S},\varphi_{\rm S})\,,$ (92)

where we have defined $C_{l}={4\pi}/(2l+1)$. Using the definition of spherical
harmonics

$Y_{lm}(\theta,\varphi)=N_{lm}P_{lm}(\cos\theta)e^{im\varphi},$ (93)

with

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

and the associated Legendre functions

$P_{lm}(\cos\theta)=(-1)^{m}(1-\cos^{2}\theta)^{\frac{m}{2}}\frac{d^{m}}{dx^{m}}P_{l}(\cos\theta),$
(95)

we re-write Eq.(92) as:

$\displaystyle
F=\sum_{l=0}^{\infty}\sum_{m=-l}^{m=l}\gamma_{lm}(\theta)\cos\left[m(\sigma_{11}t+\varphi)\right]$
$\displaystyle+i\sum_{l=1}^{\infty}\sum_{m=1}^{m=l}\sin\left[m(\sigma_{11}t+\varphi)\right]\left[\gamma_{lm}(\theta)-\gamma_{l-m}(\theta)\right]\,,$
(96)

where $\gamma_{lm}(\theta)=C_{l}F_{l}N_{lm}^{2}P_{lm}(\cos\theta)P_{lm}(0)$,
and $\sigma_{11}=\Omega-n_{\star}$. Next we note that

$\displaystyle\gamma_{lm}(\theta)-\gamma_{l-m}(\theta)=C_{l}F_{l}\Big{[}N_{lm}^{2}P_{lm}(\cos\theta)P_{lm}(0)$
$\displaystyle-N_{l-m}^{2}P_{l-m}(\cos\theta)P_{l-m}(0)\Big{]}=0,$

by virtue of the equality

$P_{l-m}=\frac{(l-m)!}{(l+m)!}P_{lm}.$ (97)

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

$F=\Re\left\\{\sum_{l=0}^{\infty}\sum_{m=0}^{l}(2-\delta_{0m})\gamma_{lm}(\theta)e^{im(\sigma_{11}t+\varphi)}\right\\},$
(98)

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

${F}_{22}(\theta,\varphi,t)=\frac{5}{64}F_{\star}P_{22}(\cos\theta)e^{i(\sigma_{22}t+2\varphi)},$
(99)

or in spherical harmonics as

${F}_{22}(\theta,\varphi,t)=\frac{\sqrt{30\pi}}{16}F_{\star}Y_{22}(\theta,\varphi)e^{i\sigma_{22}t},$
(100)

where $\sigma_{22}=2(\Omega-n_{\star})$ is the semi-diurnal tidal forcing
frequency.

##  Appendix F Energy carried by tidal waves

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

$\mathcal{W}=\frac{1}{2}p_{0}v_{0}\Re\left\\{\tilde{\delta
p}\tilde{V}_{r}^{*}\right\\}.$ (101)

By substituting $\tilde{V}_{r}$ by its expression in Eq.(48), we notice that

$\Re\left\\{\tilde{\delta
p}\tilde{V}_{r}^{*}\right\\}=\eta^{-1}\Re\left\\{\tilde{\delta
p}\tilde{G}^{*}\right\\}.$ (102)

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

$\displaystyle\mathcal{W}$
$\displaystyle=\frac{p_{0}v_{0}}{2\eta\beta\Lambda}\Im\left\\{\frac{d\tilde{G}}{dx}\tilde{G}^{*}\right\\}$
$\displaystyle+\frac{p_{0}v_{0}}{2\eta\beta}\Re\left\\{\left(1+\frac{1}{\beta\Lambda}\frac{d}{dx}\right)\tilde{U}\tilde{G}^{*}\right\\}.$
(103)

Ignoring the gravitational tidal potential $\tilde{U}$, and using the
variables $\Psi$ and $\Phi$ of the wave equation in the main text, we re-write
the energy flow as

$\mathcal{W}=\frac{p_{0}v_{0}}{2\eta\beta\Lambda}\Im\left\\{\frac{d\Psi}{dx}\Psi^{*}-\frac{A}{2}|\Psi|^{2}\right\\}|\Phi|^{2},$

which scales as $|\Psi|^{2}|\Phi|^{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 $\mathcal{V}$ can be expressed as

$\mathcal{T}=\int_{\mathcal{V}}\partial_{\varphi}U_{\rm T}d\mu,$ (104)

where $U_{\rm T}(\vec{r},t)$ represents the tidal potential, $\vec{r}$ being
the separation between the tidal players, and $d\mu$ 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

$\mathcal{T}=\lim_{P\rightarrow+\infty}\frac{1}{P}\int_{0}^{P}\\!\\!\int_{\mathcal{S}}\partial_{\varphi}U_{\rm
T}\frac{R_{\rm p}^{2}\delta p_{\rm s}}{g}d{S}dt,$ (105)

where $\mathcal{S}$ is a sphere of unit radius, and $d{S}=\sin\theta d\theta
d\varphi$. Due to the periodic nature of the forcing and the response, the
tidal potential and the pressure anomaly are expanded in Fourier series with
$\sigma$ being the tidal frequency such that

$\displaystyle U_{\rm T}(\vec{r},t)=\Re\left\\{\sum_{\sigma}U_{\rm
T}^{\sigma}(\vec{r})e^{i\sigma t}\right\\},$ (106) $\displaystyle\delta p_{\rm
s}(\theta,\varphi,t)=\Re\left\\{\sum_{\sigma}\delta p_{\rm
s}^{\sigma}(\theta,\varphi)e^{i\sigma t}\right\\}.$ (107)

As such, the torque is rewritten as

$\mathcal{T}=\frac{R_{\rm
p}^{2}}{2g}\sum_{\sigma}\Re\left\\{\int_{\mathcal{S}}{\partial_{\varphi}^{*}U_{\rm
T}^{\sigma*}}\delta p_{\rm s}^{\sigma}dS\right\\}.$ (108)

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

$\displaystyle U_{\rm
T}^{\sigma}(\vec{r})=\sum_{l=2}^{+\infty}\sum_{m=-l}^{l}U_{{\rm
T};l}^{m,\sigma}(r)Y_{l,m}(\theta,\varphi),$ (109) $\displaystyle\delta p_{\rm
s}^{\sigma}(\theta,\varphi)=\sum_{l=2}^{+\infty}\sum_{m=-l}^{l}\delta p_{{\rm
s};l}^{m,\sigma}Y_{l,m}(\theta,\varphi).$ (110)

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

$\mathcal{T}=\frac{R_{\rm
p}^{2}}{2g}\sum_{\sigma}\sum_{l=2}^{+\infty}\sum_{m=-l}^{l}m\Im\left\\{{U_{{\rm
T};l}^{m,\sigma}}\delta p_{{\rm s};l}^{m,\sigma}\right\\}.$

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

$U_{{\rm T};2}^{2,\sigma}=\sqrt{\frac{6\pi}{5}}\frac{GM_{\star}}{a_{\rm
p}}\left(\frac{R_{\rm p}}{a_{\rm p}}\right)^{2},$ (111)

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

$\mathcal{T}=\sqrt{\frac{6\pi}{5}}\frac{M_{\star}}{M_{\rm p}}\frac{R_{\rm
p}^{6}}{a_{\rm p}^{3}}\Im\\{\delta p_{{\rm s};2}^{2,\sigma}\\}.$ (112)

##  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. (3.1) 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. (3.1), is given by

$\displaystyle\tilde{\delta p}_{\rm s}=$ $\displaystyle\frac{\alpha_{\rm
A}\delta F_{22}}{p_{s}}\frac{\kappa g\Lambda_{2}}{R_{\rm
p}^{2}\sigma^{3}}\frac{\left[\mathcal{X}_{\rm R}+\mathcal{Y}_{\rm
R}\alpha+i(\mathcal{X}\alpha+\mathcal{Y})\right]\alpha}{\left(1+2\zeta+2\zeta^{2}\right)}$
$\displaystyle\times\left[(\kappa-\beta\Lambda_{2}+1)^{2}+\alpha^{2}(\beta\Lambda_{2}-1)^{2}\right]^{-1},$
(113)

where $\mathcal{X}$ and $\mathcal{Y}$ are defined by Eq. (23), and
$\mathcal{X}_{\rm R}$ and $\mathcal{Y}_{\rm R}$ are defined by

$\displaystyle\mathcal{X}_{\rm R}$
$\displaystyle=(\kappa\\!-\\!\beta\Lambda_{2}\\!+\\!1)\\!\left[2\zeta^{2}(1-\mu_{\rm
gr})\\!+\\!\zeta(2\\!-\\!\mu_{\rm gr})+1\right],$
$\displaystyle\mathcal{Y}_{\rm R}$ $\displaystyle=s\mu_{\rm
gr}\zeta(\beta\Lambda_{2}-1).$ (114)

Next, we reformulate this complex solution as

$\tilde{\delta p}_{\rm s}=K\left(\frac{\sigma_{\rm
L}}{\sigma}\right)^{3}\frac{\mathcal{G}(\sigma)}{\alpha^{-1}\left(\kappa+1-\beta\Lambda_{2}\right)+i\left(\beta\Lambda_{2}-1\right)}\,,$
(115)

where we introduced $\sigma_{\rm L}=\sigma_{\rm w}\sqrt{\Lambda_{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
$\mathcal{G}(\sigma)$ defined as

$\mathcal{G}(\sigma)=\frac{2\zeta^{2}(1-\mu_{\rm gr})+\zeta(2-\mu_{\rm
gr})+1-is\mu_{\rm gr}\zeta}{1+2\zeta+2\zeta^{2}}\,,$ (116)

and the dimensionless constant factor $K$ defined as

$K=\frac{\alpha_{\rm A}\delta F_{22}}{p_{s}}\frac{\kappa g\Lambda_{2}}{R_{\rm
p}^{2}\sigma_{\rm L}^{3}}\,.$ (117)

The transfer function $\mathcal{G}$ thus describes the delayed forcing of the
atmosphere due to thermal inertia near the planetary surface. Therefore, this
function reduces to $\mathcal{G}\left(\sigma\right)=1$ in the case of direct
absorption of the incident flux and the synchronous response of the surface
(namely, in the limit of $\zeta=\mu_{\rm gr}=0$, which is the limit with which
we describe the Earth in the main text). With $\mathcal{G}(\sigma)$ 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.

Figure 8: A sample of the tidal response, plotted in terms of modulus, real
part, and imaginary part of the surface pressure anomaly, all given by Eqs.
(119). The function $f(x)$ on the y-axis indicates the functional form of the
pressure variation, that is, the pressure function divided by the constant
factor K. The x-axis 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=~{}\sigma/\sigma_{\rm L}$ and the
dissipation parameter $\varepsilon=~{}\sigma_{0}/\sigma_{\rm L}\ll 1$, and we
rewrite Eq. (115) as

$\tilde{\delta p}_{\rm
s}=K\frac{\mathcal{G}(X)}{\varepsilon\left[(\kappa+1)X^{2}-1\right]+iX\left(X^{2}-1\right)}\,.$
(118)

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

$|\tilde{\delta p}_{\rm
s}|=K\frac{|\mathcal{G}|}{\sqrt{\varepsilon^{2}\left[(\kappa+1)X^{2}-1\right]^{2}+X^{2}\left(X^{2}-1\right)^{2}}}\,,$
$\Re\\{\tilde{\delta p}_{\rm
s}\\}=K\frac{\varepsilon\left[(\kappa+1)X^{2}-1\right]\mathcal{G}_{\rm
R}+X(X^{2}-1)\mathcal{G}_{\rm
I}}{{\varepsilon^{2}\left[(\kappa+1)X^{2}-1\right]^{2}+X^{2}\left(X^{2}-1\right)^{2}}}\,,$
$\Im\\{\tilde{\delta p}_{\rm
s}\\}=K\frac{\varepsilon\left[(\kappa+1)X^{2}-1\right]\mathcal{G}_{\rm
I}-X(X^{2}-1)\mathcal{G}_{\rm
R}}{{\varepsilon^{2}\left[(\kappa+1)X^{2}-1\right]^{2}+X^{2}\left(X^{2}-1\right)^{2}}}\,,$
(119)

where $\mathcal{G}_{\rm R}=\Re\\{\mathcal{G}\\}$ and $\mathcal{G}_{\rm
I}=\Im\\{\mathcal{G}\\}$. In Figure 8 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. (119) in the low frequency limit
($X\rightarrow 0)$ and setting $\mathcal{G}=1$ we obtain

$\Im\\{\tilde{\delta p}_{\rm s}\\}=K\frac{X}{\varepsilon^{2}+X^{2}}\,,$ (120)

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
$\alpha_{\rm A}$. As defined in the main text (Section 2.2), $\alpha_{\rm 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 (${\rm H_{2}O}$) and
ozone ($\rm 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 $\tau_{\rm G}$, which
approximately corresponds to the temperature oscillation produced by
$J^{m,\sigma}$. The input power is thus defined as (Chapman and Lindzen, 1970,
Eqs. (143-144))

$J=\Re\left\\{\sum_{\sigma}\sum_{m}J^{m,\sigma}\left(x,\theta,\varphi\right)e^{i\sigma
t}\right\\},$ (121)

with the coefficients $J^{m,\sigma}$ given by

$J^{m,\sigma}\left(x,\theta,\varphi\right)=i\sigma C_{\rm p}\tau_{\rm
G}^{m}{\sigma}\left(x,\theta,\varphi\right),$ (122)

where the heat capacity of the gas per unit mass $C_{\rm p}=\mathcal{R}_{\rm
s}/\kappa$. Figure 3.2 of Chapman and Lindzen (1970) shows the normalised
vertical distribution of $\tau_{\rm G}^{2,\sigma}$. Following the authors
(their Eq. 163), we set $H=7.6$ km. The semidiurnal distributions of
$\tau_{\rm G}^{2,\sigma}$ for water vapour and ozone can be approximated by
the functions

$\displaystyle\tau_{\rm G,H_{2}O}\left(x,\theta,\varphi\right)$
$\displaystyle=T_{{\rm H_{2}O}}V_{{\rm H_{2}O}}\left(x\right)\sin\theta
e^{i2\varphi},$ (123) $\displaystyle\tau_{\rm
G,O_{3}}\left(x,\theta,\varphi\right)$ $\displaystyle=T_{\rm O_{3}}V_{\rm
O_{3}}\left(x\right)\sin^{2}\theta e^{i2\varphi},$ (124)

with $T_{{\rm H_{2}O}}=0.035$ K and $T_{\rm 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_{{\rm
H_{2}O}}$ of Eq.(123) is set to

$V_{\rm H_{2}O}\left(x\right)=e^{-k_{\rm H_{2}O}x},$ (125)

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

$V_{\rm O_{3}}\left(x\right)=\left\\{\begin{array}[]{ll}K_{\rm
O_{3}}\sin\left(\frac{\pi\left(x-x_{1}\right)}{x_{2}-x_{1}}\right),&x_{1}\leq
x\leq x_{2},\\\ 0&\mbox{elsewhere},\end{array}\right.$ (126)

with $K_{\rm 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 ($\theta,\varphi$), the total flux absorbed by the
atmosphere for the spherical harmonic $Y_{22}$ is defined as

$\delta\mathcal{F}\left(\theta,\varphi\right)=\delta\mathcal{F}_{22}Y_{22}\left(\theta,\varphi\right).$
(127)

The efficiency parameter $\alpha_{\rm A}$ of our model is basically the ratio
between $\delta\mathcal{F}_{22}$ and the corresponding component of the
incident Solar flux, which we retrieve in E (Eq. 100) as

$\delta F_{22}=\frac{\sqrt{30\pi}}{16}F_{\star}.$ (128)