Atmospheric stability and collapse on tidally locked rocky planets

Over large timescales, a terrestrial planet may be driven towards spin-orbit synchronous rotation by tidal forces. In this particular configuration, the planet exhibits permanent dayside and nightside, which may induce strong day-night temperature gradients. The nightside temperature depends on the efficiency of the day-night heat redistribution and determines the stability of the atmosphere against collapse. To better constrain the atmospheric stability, climate, and surface conditions of rocky planets located in the habitable zone of their host star, it is thus crucial to understand the complex mechanism of heat redistribution. Building on early works and assuming dry thermodynamics, we developed a hierarchy of analytic models taking into account the coupling between radiative transfer, dayside convection, and large-scale atmospheric circulation in the case of slowly rotating planets. There are two types of these models: a zero-dimensional two-layer approach and a two-column radiative-convective-subsiding-upwelling (RCSU) model. They yield analytical solutions and scaling laws characterising the dependence of the collapse pressure on physical features, which are compared to the results obtained by early works using 3D global climate models (GCMs). The analytical theory captures (i) the dependence of temperatures on atmospheric opacities and scattering in the shortwave and in the longwave, (ii) the behaviour of the collapse pressure observed in GCM simulations at low stellar fluxes that are due to the non-linear dependence of the atmospheric opacity on the longwave optical depth at the planet's surface, (iii) the increase of stability generated by dayside sensible heating, and (iv) the decrease of stability induced by the increase of the planet size.


Introduction
Both observations and theoretical studies achieved over the last decade suggest that rocky planets can be found around stars of different masses and they represent a large fraction of the population of exoplanets (e.g.Cumming et al. 2008;Mordasini et al. 2009;Howard 2013;Ronco et al. 2017).Many of these planets are located in the habitable zone of their host stars (e.g.Kopparapu et al. 2013), which is basically the region where the incident stellar flux and greenhouse effect enable them to sustain surface liquid water.This is particularly the case for planets detected around brown dwarfs and very-low-mass stars (e.g.Payne & Lodato 2007;Raymond et al. 2007;Kopparapu et al. 2017).
Three out of seven of the Earth-sized rocky planets hosted by the TRAPPIST-1 ultra-cool dwarf star (planets e, f, and g) that were discovered in 2017 are within the habitable zone (Gillon et al. 2017;Grimm et al. 2018).Recently, two temperate Earth-mass planet candidates were detected around the Mdwarf Teegarden's Star by the CARMENES spectrograph using radial-velocity measurements (Zechmeister et al. 2019), while an Earth-sized planet orbiting in the habitable zone of the TOI-700 red dwarf was discovered by the TESS observatory (Gilbert et al. 2020;Rodriguez et al. 2020).This planet (TOI-700 d) is the first of its kind detected by TESS and it might harbor temperate surface conditions (Suissa et al. 2020).The number of discovered temperate rocky planets will keep growing in the future with the upcoming transit searches of the TESS (Barclay et al. 2018) and PLATO (Ragazzoni et al. 2016) observatories, which opens a new era of atmospheric characterisation and motivates a better understanding of the mechanisms that govern planetary climates and surface conditions.
Most planets located in the habitable zone of sub-Solarmass stars, and particularly those orbiting very-low-mass dwarfs such as TRAPPIST-1, are expected to be circularised and tidally locked in spin-orbit synchronous rotation (e.g.Kasting et al. 1993) or in a spin-orbit resonance (e.g.Correia et al. 2014) over short timescales.As rocky planets locked into synchronous rotation exhibit permanent dayside and nightside, their climate is affected by a strong gradient of thermal forcing along the axis connecting the sub-stellar and anti-stellar points.
Since it is continually heated by the stellar incident flux, the dayside is hotter than the nightside, which is cooled by infrared radiation towards space.If the nightside temperature goes below the condensation temperature of greenhouse gases present in the atmosphere, then the nightside becomes a cold trap for these gases, which condense and form an ice sheet at the planet surface.This triggers atmospheric collapse (e.g.Joshi et al. 1997;Heng & Kopparla 2012;Wordsworth 2015): the greenhouse effect diminishes and the atmosphere cools down, thereby accelerating the condensation process in the cold trap and leading the atmospheric composition, climate, and surface conditions to change radically.It is thus crucial to characterise preliminarily the atmospheric stability of rocky planets against collapse to bet-ter constrain the surface conditions of the Earth-like exoplanets found in the habitable zone of their host star.Kasting et al. (1993), who investigated the habitability conditions around main sequence stars, initially raised the question of atmospheric stability.This question was first addressed using 3D global climate models (GCMs), which are codes that solve the coupled angular momentum, mass conservation, radiative transfer, and energy equations in three dimensions on the surface of a sphere.In their pioneering works, Joshi et al. (1997) and Joshi (2003) examine the case of rocky planets orbiting M-Dwarfs and hosting CO 2 /H 2 O atmospheres, while Merlis & Schneider (2010) characterise the large-scale circulation patterns of Earthlike planets tidally locked in spin-orbit synchronisation in various regimes.
A series of studies continued on this path, such as Heng et al. (2011b) and Heng & Vogt (2011), who characterise the atmospheric dynamics of a hypothetical tidally locked Earth and treat the case of the Gliese 581g super-Earth as a scaled-up version of Earth.By considering the effects of moist thermodynamics, Leconte et al. (2013) study the circulation patterns of closein extrasolar terrestrial planets in the presence of clouds.They thus show evidence of the climate moist bistability, which is the fact that water can be either vaporised -which acts to increase the runaway greenhouse effect -or captured in permanent cold traps, leading the atmosphere to collapse.Carone et al. (2016) examine how surface friction affects the atmospheric stability, while Turbet et al. (2018) computed stability diagrams for TRAPPIST-1 planets assuming N 2 /CO 2 atmospheric mixtures.Similarly, Wordsworth (2015) and Koll & Abbot (2016) treat the case of CO, CO 2 , and N 2 -dominated dry atmospheres, and Kang & Wordsworth (2019) investigate the effect of shortwave absorption on collapse.
Meanwhile, with the rise of next-generation space telescopes, such as the James Webb Space Telescope (JWST; Deming et al. 2009), the mechanism of global heat redistribution itself was examined in order to decipher the upcoming thermal phase curves of tidally locked extrasolar rocky planets and hot Jupiters (e.g.Seager & Deming 2009;Cowan & Agol 2011;Selsis et al. 2011;Perez-Becker & Showman 2013;Koll & Abbot 2015;Koll & Komacek 2018;Koll 2019).Although their goals are slightly different, these works are complementary with those dealing with atmospheric collapse.
While they are convenient to treat the complex non-linear physics governing planetary climates, GCMs are far too costly to explore a broad range of the parameter space, with any single climate simulation requiring days or weeks of CPU time.Furthermore, the asymptotic regimes of the steady states reached after convergence may actually be characterised using a relatively small set of nondimensional control parameters (e.g.Koll & Abbot 2015).This motivated the development of simplified analytic theories reducing both the physics complexity and the dimensionality of the problem.
First, 1D analytic models were introduced to characterise the thermal structure of planetary atmospheres.Guillot (2010) derived temperature-pressure (T ´P) profiles in the limit of the 'double-gray' approximation and pure absorption.Robinson & Catling (2012) considered T ´P profiles in the presence of convection.Heng et al. (2012) generalised the work of Guillot (2010) to include isotropic scattering.Parmentier & Guillot (2014) used the 'picket fence' model of Chandrasekhar (1935) and Mihalas (1978) to generalise the work of Guillot (2010) by including the simplified non-grey radiative transfer due to spectral lines versus continua.Heng et al. (2014) generalised the work of Heng et al. (2012) to include non-isotropic/anisotropic scattering, and also demonstrated that the governing equations of these T ´P profiles and the two-stream solutions have a common origin.Mohandas et al. (2018) generalised the work of Parmentier & Guillot (2014) to include isotropic scattering.
Second, the analytic theory of global heat redistribution was developed incrementally through a diversity of approaches starting from analytical scalings of the atmospheric stability of super-Earths orbiting M-dwarf stars (e.g.Heng & Kopparla 2012).Laying the foundation of the 0-D theory, Wordsworth (2015) proposed a two-layer model that self-consistently includes both radiative transfer and sensible heat exchanges with planetary surface due to convection within the dayside planetary boundary layer (PBL).In this box model, the atmosphere is treated as a globally isothermal and optically thin layer, and the complex features of the three-dimensional general circulation are ignored.Despite these simplifications, the model approximately matches results obtained from GCM simulations for CO 2 -dominated planets.Particularly, it convincingly captures the behaviour of the collapse (or critical) pressure p crit -that is the minimum surface pressure for atmospheric stability -at high stellar fluxes.
Comparable progress was achieved by Koll & Abbot (2016) regarding the 1D theory.Building on the early work by Robinson & Catling (2012), these authors developed a two-column radiative-convective model that takes the effect of subsidence (i.e.downwelling flows) on the nightside temperature profile into account.This model is thus referred to as a radiative-convectivesubsiding (RCS) model in the following.By performing a series of simulations with a GCM, Koll & Abbot (2016) showed that the large-scale atmospheric circulation acts as a global heat engine (a similar work was done by Koll & Komacek 2018, in the case of hot Jupiters), and used this statement to derive scalings for the typical wind speeds parametrising heat fluxes.They thus recovered the stability diagrams of CO 2 -dominated atmospheres previously obtained by Wordsworth (2015) from GCM simulations.
In the present work, we aim to consolidate the analytic theory by building mainly on the developments made by Wordsworth (2015) and Koll & Abbot (2016).Particularly, our goal is to elucidate analytically how the collapse pressure is affected by the non-linear pressure dependence of the atmospheric optical thickness at low stellar fluxes, by dayside convection, and by largescale advection.Therefore, we adopt the 0-D approach proposed by Wordsworth (2015) and introduce incrementally in the model the effects of longwave and shortwave absorption and scattering, turbulent heat exchanges with the planet surface -which we call 'sensible heating' -and advective heat transport by a stellar and anti-stellar cell, which is the regime of slowly rotating planets.Following Koll & Abbot (2016), scalings of wind speeds are derived in the framework of heat engine theory.With every step, intermediate results are compared with those obtained in earlier studies.The full model -which encompasses radiative transfer, dayside convection, and large-scale circulation -is formally written as one single equation, given by Eq. ( 78), and controlled by a small set of dimensionless parameters.This model is used to compute stability diagrams, and finally compared to a 1D two-column radiative-convective-subsiding-upwelling (RCSU) model that we develop by introducing dayside upwelling flows in the RCS model of Koll & Abbot (2016).
In Sect.2, we detail the physical setup of the 0-D model and discuss the main assumptions.In Sect.3, we develop a twolayer radiative model of global heat redistribution including absorption and scattering in the two-stream, dual-band, and grey gas approximations.We then derive analytically the shortwave and longwave flux profiles, radiative transmission functions, and equilibrium temperatures.In Sect.4, the dayside sensible heating is included in the model.In Sect.5, we relax Wordsworth's globally isothermal atmosphere assumption and make the daynight temperature gradient depend on the heat transport due to stellar and anti-stellar circulation.In Sect.6, the model is used to compute stability diagrams for CO 2 -dominated atmospheres derived from the case treated by Wordsworth (2015).We thus show that the model captures the evolution of the collapse pressure at low stellar fluxes observed in GCM simulations, and the fact that the atmospheric stability increases with sensible heating and decreases as the planet's size increases.In Sect.7, we introduce our two-column RCSU model and discuss the role played by the atmospheric structure and, finally in Sect.8 we summarise the conclusions and future works.

Physical setup and main assumptions
We introduce in this section the main features of the 0-D model developed in the study as well as frequently used notations.This provides a global overview of the mechanisms involved in daynight heat redistribution of terrestrial planets that we take into account (Fig. 1).These mechanisms are detailed in the next sections.

Tidal locking in 1:1 spin-orbit resonance
We consider the simplified case of a tidally-locked rocky planet, of mass M p and radius R p , synchronously moving on a coplanar and circular orbit of period P ‹ .This regime is mainly relevant for planets orbiting small main-sequence stars, such as K and M dwarfs, since the habitable zone of the host star is within the tidal lock radius in this case (Kasting et al. 1993;Edson et al. 2011).This result may be derived from a simple scale analysis.
First, the tidal lock radius is usually estimated for dry planets from the formula (Peale 1977;Kasting et al. 1993;Dobrovolskis 2009;Edson et al. 2011), Here, M ‹ is the stellar mass, P rot;0 the original rotation period of the planet, t pf the time period from formation, and Q the tidal quality factor accounting for tidal dissipation in the planet's interior (the smaller Q and the larger tidally dissipated energy; Goldreich & Soter 1966).All quantities are in CGS units (centimetres, grams, seconds) in the formula.Second, the typical radius of the habitable zone r HZ can be roughly defined as the star-planet distance at which the black body equilibrium temperature of the planet has a given value T HZ .The equilibrium temperature is given as a function of the stellar flux received by the planet F ‹ , the parameter σ SB " 5.670367 ˆ10 ´8 W m ´2 K ´4 (Mohr et al. 2016) being the Stefan-Boltzmann constant, and the symbol " meaning 'defined by'.The above definition of r HZ thus leads to F ‹ " 4σ SB T 4 HZ .The stellar flux is also expressed as a function of the orbital radius a and the stellar bolometric luminosity L ‹ , By combining the two expressions of F ‹ , it follows that (4) The dependence of the stellar bolometric luminosity of mainsequence stars on the stellar mass has been quantified by empirical formulae, such as (Barnes et al. 2008) with µ " log 10 pM ‹ {M d q, the notations M d and L d referring to the Solar mass and luminosity, respectively.A linear regression of log 10 ´aL ‹ {L d ¯on the range ´1 ď µ ď 0 yields and thus r HZ 9 M 1.32 ‹ .By comparing this scaling law to r T 9 M 1{3 ‹ (Eq.( 1)), we observe that the tidal lock radius (r T ) increases slower with the stellar mass than the characteristic radius of the habitable zone (r HZ ).This means that the probability to be tidally locked in synchronous rotation for a planet orbiting in the habitable zone decays as the stellar mass increases, which is illustrated for instance by Kasting's plot (Kasting et al. 1993, Fig. 16).
In the circular coplanar regime, there is no seasonal variation of the local incident stellar flux because of the absence of obliquity and eccentricity.The stellar forcing at the top of the planet's atmosphere is thus invariant with time.The spin rotation of the planet is defined by the spin angular velocity Ω " n ‹ , where n ‹ " 2π{P ‹ is the orbital frequency.Hence, Ω is determined using the third Keplerian law, where G designates the universal gravitational constant.Besides, considering both the expression of the stellar flux received by the planet, given by Eq. ( 3), and the fact that the stellar bolometric luminosity L ‹ of main-sequence stars may be written as a function of the stellar mass (see Eq. ( 5)), we remark that the planet's angular velocity is actually fully constrained by the stellar mass and flux through the relationship and thus cannot be taken as a free parameter.Particularly, in simulations that we perform with general circulation models to benchmark the theory, the planet's rotation rate is specified as a function of M p , M ‹ , and F ‹ using the preceding expression with the scaling law of L ‹ given by Eq. ( 5).

Radiative transfer
The planet's atmosphere is considered in the standard shallow atmosphere framework (Vallis 2006) 1 , where the thickness of the fluid layer is assumed to be far smaller than horizontal scales, comparable to the planet's radius in order of magnitude.In this case, the horizontal propagation of radiation may be neglected and radiation only propagates upwards and downwards.Upwelling fluxes are denoted by F Ò and downwelling fluxes by F Ó .This sets the bases of the two-stream approximation (e.g.Heng 2017, Sect. 3.1), which is assumed to derive radiative net fluxes F ´" F Ò ´FÓ in the next section.
In addition to the two-stream approximation, we consider that the frequency spectra of the stellar and planetary radiative fluxes do not overlap, which is the so-called dual-band approximation (Heng 2017, Sect. 4.1).The radiative fluxes emitted by the star and the planet are thus split into two decoupled components, which are refereed to as the 'shortwave' and 'longwave' fluxes, and subscripted by S and L, respectively.While this simplification holds for rocky planets hosted by Sun-like stars, where the bodies' surface temperatures are separated by more than one decade in orders of magnitude, it tends to be less adapted to those hosted by ultra-cool dwarf stars, such as TRAPPIST-1, although the stellar temperature (« 3000 K) is still far higher than temperate planetary atmospheres (« 300 500 K).In these cases, the stellar and planetary radiative fluxes partly overlap.Cool stars are also those that are more likely to host tidally-locked exoplanets in their habitable zone than Sunlike stars, as stated in Sect.2.1.We should thus bear in mind the limitations of the dual-band approximation when applying the model to such systems.
Finally, grey gas opacities are assumed, meaning that atmospheric opacities in the longwave and shortwave are wavelengthindependent in the model (Heng 2017, Sect. 4.1).In each band, absorption is described by a unique parameter that we call 'effective opacity' to emphasise the fact that it does not really cor-1 The term 'shallow' is used here in the context of radiative transfer, where it just refers to the ratio between the vertical and horizontal scales of the fluid layer.Its meaning is more subtle in the context of fluid dynamics since it may designate different approximations in this case (see e.g.Vallis 2006, Chapter 3).respond to a mean opacity owing to the existing correlation between radiative spectra and opacity lines (see e.g.Wordsworth 2015, Fig. 10).Because of this correlation and of the strong wavelength dependence of opacity lines, this hypothesis appears as one of the strongest regarding radiative transfers.Particularly, as stated by Leconte et al. (2013) from the comparison of GCM simulations using different approaches, the cooling of the planet's surface is considerably underestimated when grey opacities are used, meaning that the nightside temperature is overestimated.
This is what motivated the strong effort made to include refined treatments of radiative transfer (e.g. the k-correlated distribution method; Lacis & Oinas 1991) in GCMs used to study the heat redistribution on tidally-locked exoplanets (e.g.Leconte et al. 2013;Wordsworth 2015).In the present work however, we have to admit the grey gas assumption as the price to pay for the simplification of the theoretical analysis.For exhaustive discussions of its limitations, as well as those of the shallow-water, two-stream, and dual-band approximations, one may consult the reference books by Seager (2010), Pierrehumbert (2010), andHeng (2017).We note that, in analytic developments, we follow the conventions and notations employed by Heng (2017) in Chapters 3 and 4.
As shown by Koll & Abbot (2016) through the development of a two-column 1D model, the self-consistent treatment of the coupling between radiative transfer and the atmospheric structure is a major source of complexity, mainly because of the degeneracies affecting boundary conditions.As our goal is to introduce simplified atmospheric dynamics in the analytical theory, we make the choice to favour the atmospheric circulation over the atmospheric structure, noting that the atmospheric dynamics already induces non-negligible mathematical complications.Hence, following the early work by Wordsworth (2015), we opt for a zero-dimensional two-layer model where the atmosphere is considered as isothermal across the vertical direction.However, we relax the well-mixed atmosphere approximation made by Wordsworth (2015), and assume that dayside and nightside atmospheric temperatures, denoted by T a;d and T a;n , are not the same in the general case.Similarly, the planet surface is characterised by its dayside and nightside temperatures, T s;d and T s;n (Fig. 1).
Article number, page 4 of 27 P. Auclair-Desrotour and K. Heng: Atmospheric stability and collapse on tidally locked rocky planets

Atmospheric heat transport
The four temperatures of the system are coupled together by the mechanisms of heat transport taken into account in the model, which are of three types: (i) radiative transfers along the vertical direction (shortwave and longwave absorption, radiation, and scattering); (ii) surface-atmosphere turbulent exchanges in the dayside convective boundary layer; and (iii) day-night horizontal heat transport.Each of these components may defined separately.Owing to the above simplifications, radiative transfers are described by simple analytic solutions, which are derived in Sect.3.For the sensible heating due to turbulent heat transport, among the different prescriptions existing in literature, we choose to follow that proposed by Koll & Abbot (2016), which treats the atmosphere as a heat engine where the convective flow acts against friction in the surface boundary layer.This part of the theory is detailed in Sect. 4.
The large-scale day-night heat transport turns out to be the most complex effect to include since there is no general theory of atmospheric heat flux to our knowledge.This is due to the fact that the mechanisms contributing to the heat transport cannot be decoupled and are spanning over the three dimensions at planetary scales.Nevertheless, as demonstrated by Koll & Abbot (2015), the involved processes may be disentangled by making use of the Buckingham-Pi theorem (Buckingham 1914), which unravels the nondimensional parameters governing the system.
Among these parameters, we find the nondimensional Rossby deformation length, the length L Ro being the Rossby radius of deformation (see e.g.Vallis 2006, Sect. 3.8.2),and c wave the characteristic speed of gravity waves, defined hereafter.The nondimensional Rossby deformation length governs the circulation regime of the atmosphere.In the case of slow rotators, LRo Á 1, meaning that waves can propagate planetwide.Because the effect of rotation is weak, the circulation that develops in this regime is approximately symmetric with respect to the axis connecting the substellar and anti-stellar points, the day-night flow being driven by the balance between advection and pressure-gradient acceleration (e.g.Leconte et al. 2013).In this stellar and anti-stellar circulation, high-altitude winds blow from the dayside to the nightside and low-altitude wind from the nightside to the dayside (e.g.Merlis & Schneider 2010;Leconte et al. 2013).
We illustrate the slow rotators regime by performing simulations with the Open-source General Circulation Model THOR (Mendonça et al. 2016;Deitrick et al. 2019), which solves the general non-hydrostatic equations within a spherical shell using the icosahedral grid.In these calculations, simple grey gas radiative transfers are used and two cases treated by early works are reproduced: a 0.1-bar CO 2 -dominated atmosphere at T eq " 279 K (Wordsworth 2015, Fig. 2), and a 0.5-bar N 2 -dominated atmosphere at T eq " 400 K (Koll & Abbot 2016, Fig. 4). Figure 2 shows instantaneous snapshots of the zonally-averaged temperature and vertical wind speeds in both cases.These quantities are plotted as functions of latitude and pressure in the reference frame where the north pole is located at the substellar point and the south pole at the anti-stellar point.We emphasise that this latitude coordinate is different from the usual latitude, which is defined from the planet axis of rotation.As we observe that a steady cycle has been reached in both simulations after 1000 days, snapshots are taken at this date.The plots particularly highlight the day-night temperature gradient (left panels) and the dayside convective cell (right panels), which is characterised by strong upward flows in the substellar region.
When the Rossby deformation radius becomes less than the planet radius ( LRo ă 1), the circulation regime changes and the atmospheric dynamics starts developing super-rotation, that is strong eastward equatorial jets (Showman & Guillot 2002;Showman & Polvani 2011;Thrastarson & Cho 2010;Heng & Vogt 2011;Heng et al. 2011a;Leconte et al. 2013;Mendonça 2019).As shown by Showman & Polvani (2011), the emergence of super-rotation is due to the formation of standing, planetaryscale equatorial Rossby and Kelvin waves (i.e.waves restored by the Coriolis acceleration; see e.g. Lee & Saio 1997), which feed zonal flows by continually pumping angular momentum from the mid-latitudes towards the equator.The latitude width of the equatorial band where super-rotating flows are accelerated through this mechanism scales as the equatorial Rossby radius (L Ro ), and thus decays as the planet angular velocity increases.As the Rossby and Kelvin waves responsible for the super-rotation are caused by the latitudinal variations in radiative heating, the strength of equatorial jets increases with the day-night temperature difference (Showman & Polvani 2011).
The speed of gravity waves that intervenes in the nondimensional Rossby deformation length depends on the atmospheric structure.Denoting by g the surface gravity, R GP the ideal gas constant, and M a the mean molecular weight of the atmosphere, we first introduce the specific gas constant R s " R GP {M a and the pressure scale height, given as a function of temperature T by Gravity waves are restored by the Archimedean force associated with the fluid buoyancy (Gerkema & Zimmerman 2008).The strength of this force is quantified by the Brunt-Väisälä frequency N, which, in a dry stably stratified atmosphere, is given by where C p designates the heat capacity per unit mass of the gas at constant pressure, and z the altitude.In terms of H and N, the typical speed of gravity waves is simply written as c wave " HN (e.g.Koll & Abbot 2015;Leconte et al. 2013).Thus, substituting these parameters by the above expressions (Eqs.( 10) and ( 11)) in Eq. ( 9) and assuming an isothermal atmosphere yields LRo " In the present work, we only consider the heat advected by a day-night stellar and anti-stellar circulation, which corresponds to the asymptotic regime of slowly rotating planets ( LRo " 1).We do not include other mechanisms of day-night heat transport, such as that due to the propagation of gravity waves, which was already studied in early works (Koll & Abbot 2015, 2016).Thus, in Sect.5, a rough approximation of the day-night heat flow is derived from a scale analysis in the heat engine framework (Koll & Abbot 2016) using for wind speeds the prescription given by Article number, page 5 of 27 Fig. 2. Instantaneous snapshots of zonally-averaged temperature (left panels) and vertical wind speed distributions (right panels) as functions of latitude and pressure with north pole at substellar point.Top: 0.1-bar CO 2 -dominated atmosphere at T eq " 279 K (Wordsworth 2015, Fig. 2).Bottom: 0.5 bar N 2 -dominated atmosphere at T eq " 400 K (Koll & Abbot 2016, Fig. 4).Latitude is given in degrees (horizontal axis).This coordinate does not correspond to the usual latitude but to a latitude coordinate in the system of spherical coordinates where the north pole is the substellar point and the south pole the anti-stellar point.The associated colatitude is the stellar zenithal angle in this system of coordinates.Pressure is given in bars (vertical axis), temperature in K, and vertical velocity in m s ´1.The white area at the substellar point indicates the region where surface pressure is lower than the spherically averaged surface pressure due to the updraft.The dotted red line indicates the transition between rising and subsiding flows, where the vertical velocity is zero.Koll & Komacek (2018) in the case where circulation is balancing Rayleigh drag.As predicted by dimensional analyses (e.g.Leconte et al. 2013, Sect. 3.2.1),this heat flow is scaled by the ratio t rad {t adv between the radiative and advective timescales (e.g.Showman & Guillot 2002) t rad " p s C p gσ SB T 3 eq , and t adv " where p s designates the atmospheric surface pressure and V adv the typical speed of stellar and anti-stellar advective flows.This accounts for the fact that the heat transported by weak flows (t rad {t adv ! 1) is radiated towards space before being advected to the nightside.Conversely, strong flows (t rad {t adv " 1) efficiently transport heat since the advected fluid parcels reach the nightside before being radiatively cooled.In the present work, the scaling law chosen for the advected heat flux is different from that used to describe thermal exchanges with the surface on the dayside in order to emphasise the distinction between the two involved mechanisms in the general case.These mechanisms could nevertheless be linked to each other by assuming the same scalings for the associated heat fluxes.
Finally, we adopt dry thermodynamics as a first convenient step to set the basis of the model, ignoring thereby the role played by an ocean or moisture on the day-night heat redistribution.We note however that this role may be non-negligible.Edson et al. (2011) showed for instance that the presence of a slab ocean considerably affects the stability of the planet's atmosphere with respect to collapse, by inducing an additional heat transport by latent heat flux in the atmosphere and by heat diffusion in the ocean itself.Water vapour also leads to the formation of clouds, which affects the albedo of the atmosphere as well as its thermal structure and general circulation (Leconte et al. 2013).While not taken into account in the present work, these aspects may be included in the model in future studies.

A two-layer grey radiative model
Radiative transfer is the first building block of the model.We establish analytically in this section the vertical profiles of net radiative fluxes in the short-and longwave for the isothermal atmosphere.These profiles provide the integrated transfer functions that are used to derive the atmospheric and surface power budget equations.
Article number, page 6 of 27 P. Auclair-Desrotour and K. Heng: Atmospheric stability and collapse on tidally locked rocky planets

Two-stream analytic solution
In the two-stream, dual-band, and gray-gas approximations (see Sect. 2.2), the longwave and shortwave fluxes decouple and only propagate upwards or downwards.As it is governed by the same set of equations, the shortwave case is readily deduced from the longwave case that we treat hereafter.The longwave total and net fluxes are defined as F `L " F ÒL `FÓL and F ´L " F ÒL ´FÓL , respectively.They are governed by the two-stream Schwarzschild equations, which may be written as the system of first order partial differential equations (Heng 2017) where B X designates the partial derivative with respect to X, τ L the optical depth of the atmosphere in the longwave, and β L0 the scattering parameter.The optical depth characterises the optical thickness of the atmosphere and varies between 0 (top of the atmosphere) and τ s;L (planet's surface).The scattering parameter indicates the fraction of absorbed flux with respect to the scattered component and takes its values between 0 and 1 (β L0 " 1 corresponds to pure absorption and β L0 " 0 to pure scattering; see e.g.Heng 2017, Sect. 4.3).Eliminating the total flux in the preceding system yields the second order ordinary differential equation the notation d dX designating the total derivative of some quantity X.
As the isothermal approximation is assumed, the atmospheric temperature profile is approximated by the constant T a .In this particular case, the right-hand member of Eq. ( 16) vanishes, which leads to the analytic expressions of the total and net fluxes, and of the upwelling and downwelling fluxes, where A L and B L are the two integration constants of the solution, and ζ Ĺ and ζ L the coupling coefficients defined by The parameters ζ L characterise the coupling of radiative fluxes due to scattering.In the case of pure absorption (β L0 " 1), there is no coupling (ζ Ĺ " 0 and ζ L " 1).Conversely, in the case of pure scattering (β L0 " 0), the coupling is strong (ζ Ĺ " ζ L " 1{2).Similar expressions may be derived for the shortwave band by assuming T a " 0 in Eqs.(17)(18)(19)(20).The shortwave fluxes (F `S, F ´S, F ÒS , and F ÓS ) are parametrised by the as functions of optical depths (vertical axis) in case of isothermal atmosphere.Fluxes profiles are computed from the analytic formulae given by Eqs. ( 22) and ( 23) and normalised by the black body equilibrium flux F eq " σ SB T 4 eq (the equilibrium temperature T eq being given by Eq. ( 2)).Optical depths are normalised by their value at the planet's surface, namely τ s;S for the shortwave and τ s;L for the longwave.Shortwave (blue lines) and longwave (red lines) net fluxes are plotted for various optical thicknesses, the optical depths at the planet's surface taking the values τ s;S " 0.1, 1, 10 (from light to dark blue lines) and τ s;L " 0.1, 1, 10 (from light to dark red lines).Pure absorption is assumed (β L0 " β S0 " 1).The incident, atmospheric, and surface fluxes are arbitrarily set to F i " 4F eq , F a " 0.7F eq , and F s " 2F eq , respectively.
shortwave optical depth τ S (such that 0 ď τ S ď τ s;S , τ s;S being the shortwave optical depth at the planet's surface), scattering parameter β S0 , coupling coefficients ζ S " p1 ˘βS0 q {2, and integration constants A S and B S .
The four integration constants of the solution are determined by assuming two boundary conditions for each wavelength domain.At the top of the atmosphere, the shortwave downwelling flux is the stellar incident flux F i and the longwave downwelling flux is zero because of the absence of thermal source at the infinity, which is written as F ÓS " F i and F ÓL " 0. At the planet surface, a fraction of the shortwave downwelling flux is reflected while the remaining part is absorbed by the ground and re-emitted in the longwave.Denoting by A s the albedo of the planet's surface in the shortwave, and F s the flux emitted by the planet's surface in the longwave, these conditions are expressed as F ÒS " A s F ÓS and F ÒL " F s .We thus end up with the net fluxes in the short-and longwave, and Article number, page 7 of 27 where we have introduced the atmospheric and surface fluxes, F a and F s , and the vertically-integrated transmission functions T L " e ´τs;L , and T S " e ´τs;S . (24) The vertical profiles of the short-and longwave net fluxes given by Eqs. ( 22) and ( 23) are plotted in Fig. 3 for various atmospheric optical depths.One may observe here that increasing the opacity tends to confine the absorption region in the upper atmosphere in the shortwave band, while it tends to neutralise the net flux within the atmosphere in the longwave.

Surface and atmospheric energy budgets
The derived analytic solution is now used to establish the local energy budgets of the planet's surface and atmosphere at thermodynamical equilibrium.By considering the upper (τ L " τ S " 0) and lower (τ L " τ s;L and τ S " τ s;S ) boundaries, we obtain the fluxes radiated towards space at the top of the atmosphere and towards the atmosphere at the planet's surface, respectively.These fluxes are expressed as and which brings out the coefficient parameterising the surface and atmospheric energy balances, with C L " K L `AL ´1 (similarly, in the short wavelength domain, we may define C S " K S `AS ´1).These coefficients are expressed as and Thus, introducing the heat power per unit area due to surface turbulent exchanges on the dayside (F sen ) and day-night advection (F adv ) -specified in Sect. 4 and Sect.5, respectively -the local energy balances of the planet's surface and atmosphere are formulated as p1 ´AS q F i `CL F a ´KL F s ´Fsen " 0 (surf.;day), (40) C S F i ´2C L F a `CL F s ´Fadv `Fsen " 0 (atm.; day), (41) C L F a ´KL F s " 0 (surf.;night), (42) ´2C L F a `CL F s `Fadv " 0 (atm.; night).( 43) We note that Eqs. ( 40) and ( 41) are the counterparts of Eqs. ( 14) and (15) in Wordsworth (2015), respectively.Following Pierrehumbert (2011) and Wordsworth (2015), the notation F sen refers to 'sensible heat flux', which designates an energy flux corresponding to a change of temperature of the fluid2 .
As illustrated by Fig. 1, the dayside and nightside hemispheres of the planet are characterised by different surface and atmospheric temperatures (T s;d , T a;d , T s;n , and T a;n ).This leads to average the preceding energy budgets over each of the two hemispheres.Hence, we introduce the dayside and nightside hemisphere-averaged surface fluxes, denoted by B s;d and B s;n , the blackbody fluxes radiated by the atmosphere, B a;d and B a;n , and the averaged heat flows, these quantities being defined by (see Wordsworth 2015, Eqs. ( 16) and (17)) Besides, the incident stellar flux is given as a function of the stellar zenithal angle φ by Thus, averaged over the dayside and nightside hemispheres, the system of Eqs.(40-43) becomes This set of equations is the counterpart of Eqs.(21-23) in Wordsworth (2015) extended outside of the optically thin limit and with an additional energy flux between the dayside and nightside hemispheres (F adv ).

Purely radiative case
Before investigating the role played by sensible heating and advection, it is worth examining the purely radiative case.For simplification, we place ourselves in the regime treated by Wordsworth (2015), where day-night heat transfers are strong enough to homogenise the atmospheric temperature (T a;d " T a;n " T a and B a;d " B a;n " B a ).In this framework, the sensible heating is neglected (F sen " 0), and the advected heat flow does not intervene in the equilibrium, meaning that the dayside and nightside atmospheric budgets may be combined into a globally averaged atmospheric budget, Solving this equation together with Eqs. ( 46) and ( 48) yields the surface and atmospheric equilibrium temperatures which capture the non-linear dependence of radiative transfer on the atmospheric optical thickness in the general case, and highlight asymptotic regimes.Assuming pure absorption in the longwave (β L0 " 1) and transparency in the shortwave (τ s;S " 0), we recover the behaviour derived by Wordsworth (2015) in the optically thin limit (T L « 1 ´τs;L ).Conversely, as the optical depth increases, both the atmosphere and surface temperatures converge towards the black body temperature, as showed by Table 1.We note that the surface temperature of rocky planets in the optically thick limit are expected to be higher than that predicted by the model, because the atmospheric structure cannot be approximated by an isothermal temperature profile in this limit.For instance, Venus' troposphere is characterised by a strong temperature gradient (e.g.Seiff et al. 1985), which is far closer to the adiabatic profile (N " 0) than to the isothermal profile.Consequently, Venus's mean surface temperature is T s « 740 K whereas its blackbody equilibrium temperature in the absence of albedo is T eq " 327.3 K (we have taken the value F ‹ " 2601.3W m ´2 for the Solar bolometric flux).
Table 1.Optically thin and thick asymptotic limits in the pure radiative (dayside sensible heating and day-night heat transfers ignored) and pure absorption approximations (no scattering).We adopt the notations τ ‚ ! 1 and τ ‚ " 1 to designate the optically thin (τ s;S !1; τ s;L ! 1) and thick (τ s;S " 1; τ s;L " 1) limits, respectively.The case treated by Wordsworth (2015) is obtained by setting τ s;S " 0 in the optically thin limit.

Parameter
Thin (τ T s;d T eq r2 p1 ´As qs The effect that absorption and scattering have on the steady state is explored in Fig. 4, where the atmospheric and surface temperatures are plotted against the longwave opacity for various scattering parameters and shortwave opacities.We first consider the quasi-transparent limit in the shortwave (τ s;S " 10 ´4), which corresponds to τ s;S !τ s;L .As the longwave optical thickness increases, the equilibrium state switches from the optically thin regime, where T s;n 9 τ 1{4 s;L , to the optically thick one, where the three temperatures reach a plateau.The interplay between short-and longwave optical depths defines another dimension in the parameter space.When τ s;L À τ s;S , the atmosphere has to increase its temperature to evacuate the power absorbed in the shortwave.As a consequence, the atmospheric temperature increases as the optical thickness in the longwave decays, which is the behaviour observed in Fig. 4 (right panels).
Scattering in the shortwave and in the longwave have opposite effects on temperatures.When scattering in the shortwave is present (β S0 " 0.5, middle panels), it induces an anti-greenhouse cooling.A fraction of the incident stellar flux is scattered back to space, meaning that the planet absorbs less energy than in the case of pure absorption.As a consequence, the atmosphere and surface temperatures are smaller.Conversely, the presence of scattering in the longwave (β L0 " 0.5, bottom panels) generates an effect known as scattering greenhouse effect (e.g.Heng 2017, Sect. 4.7, p. 71).Thermal emissions of the surface and atmosphere towards space are scattered back to the planet, leading to the observed increase of the atmospheric, dayside and nightside surface temperatures with respect to the case of pure absorption.

Inclusion of sensible heating
As a second step, the dayside sensible heating is included in the model.While we still assume a uniform atmospheric temperature (T a;d " T a;n " T a ), the surface turbulent heat flux F sen is now expressed as a function of the dayside surface and atmospheric temperatures.We introduce first the expression of turbulent heat exchanges between the surface and the atmosphere within the planetary boundary layer (PBL), and second the scaling law of the horizontal velocity parametrising this flux, following Koll & Abbot (2016), where the atmospheric circulation is treated as a heat engine (Fig. 5).

Scaling of turbulent heat exchanges
The turbulent heat exchanges due to convection in the dayside PBL may be quantified in the framework of the mixing length theory (Holton 1973, 5 th edition, Sect.8.3.5, p. 268), as detailed in Appendix A. The hemisphere-averaged sensible heat flux is written as (e.g.Pierrehumbert 2011, Eq. (6.11), p. 396) where ρ a;d " p s { pR s T a;d q is the surface density, C D the bulk drag coefficient accounting for the strength of friction, and V sen the typical horizontal wind speed quantifying the strength of the convective cell that develops in the vicinity of the substellar point.If the ground is warmer than air, heat is carried away from the ground (F sen ą 0) at a rate proportional to the temperature difference.Conversely, if the ground is cooler than the air, the atmosphere warms the ground (F sen ă 0).We note that the convective cell vanishes in case of stable stratification, which implies V sen " 0 and F sen " 0. Thus, we ignore the sensible flux on the nightside.On the dayside however, although radiative transfer is described in the isothermal atmosphere approximation for simplification, we assume that the planetary boundary layer is convective, so that F sen ‰ 0 in the general case.
The bulk drag coefficient weakly depends on the surface and boundary layer properties.In the case of a neutrally buoyant boundary layer, following the empirical scaling law proposed by von Kármán for the mixing length (Appendix A), it is expressed as a function of the surface roughness length z 0 , the thickness of the turbulent boundary layer z h , and the von Kármán constant K « 0.4, as (e.g.Esau 2004, Eq. ( 7)) The bulk drag coefficient is of order C D " 10 ´3 (a typical value of C D over oceans is C D « 1.5 ˆ10 ´3; e.g.Holton 1973, Sect. 8.3.1),although it may be larger over rough ground.In the present work, this parameter is set to C D " 3.4 ˆ10 ´3, which is the value used by Wordsworth (2015) and corresponding to z 0 " 10 ´2 m and z h " 10 m.The characteristic horizontal wind speed varies with the surface and atmospheric temperatures.To take this dependence into account, Wordsworth (2015) derived a relationship between V sen , T a , and T s;d from the thermodynamic equation by assuming the Weak Temperature Gradient (WTG) approximation, that is considering that winds are driven by small horizontal temperature gradients.This equation (Eq.( 43) in the article) was combined with Eq. ( 54) to close the system of equations derived from hemisphere-averaged energy budgets.Later, by performing a series of numerical simulations using a GCM, Koll & Abbot (2016) showed that the atmospheric circulations of rocky planets resemble heat engines.Particularly, they established that the heat engine theory better fit horizontal wind speeds than the relationship derived by Wordsworth (2015) in the explored region of the parameter space.We thus choose to follow here the prescription proposed by Koll & Abbot (2016).
In this framework, the dayside convection is idealised as a single overturning cell between the substellar point and cooler regions.This cell is driven by the surface-atmosphere temperature difference, and the heat engine is fed by the absorption of radiative fluxes.The atmosphere absorbs heat near the dayside surface at the temperature T s;d (hot reservoir) and emits it to space at the temperature T a;d (cold reservoir) in the upper regions of the atmosphere.During the cycle, a fluid parcel works against friction in the boundary layer, which is the place where the energy is dissipated.This work is formulated as W " C D ρ a;d V 3 sen (Bister & Emanuel 1998;Koll & Abbot 2016).In the case of an isentropic cycle (i.e.composed of adiabatic reversible processes solely), W is related to the amount of power per unit area available to drive atmospheric motion Q in by Carnot's theorem (e.g.Schroeder 2000) the factor η sen " pT s;d ´Ta;d q {T s;d being the atmosphere's thermodynamic efficiency.The typical horizontal wind speed of the convective cell is thus expressed as where we have introduced the efficiency coefficient ε sen accounting for the non-isentropic nature of the thermodynamic cycle.
For an isentropic cycle, ε sen " 1, which corresponds to the idealised Carnot's heat engine.In reality, the thermodynamic cycle is not isentropic because additional dissipative processes, such as diffusion, induce an irreversible production of entropy, which decreases the efficiency of the heat engine.As a consequence, the effective amount of energy available to drive atmospheric motion is smaller than its theoretical estimate, and ε sen ă 1. Koll & Abbot (2016) hence noted that the typical wind speeds given by their GCM simulations could be twice smaller than those predicted by the theory, that is ε sen " 1{2.
The power per unit area available to drive atmospheric motion (Q in ) has to be specified from an ad hoc prescription.As it was benchmarked against GCM simulations, we follow that proposed by Koll & Abbot (2016), which is formulated as with F eq " σ SB T 4 eq , the black body equilibrium temperature T eq being defined in Eq. ( 2).
The first factor of the expression given by Eq. ( 58) is just the hemisphere-averaged power received by the planet on the dayside.The second factor is the vertically-integrated atmospheric transmission function in the shortwave.The third factor is a function of optical depth scaling the fraction of flux emitted by the surface in the longwave that is absorbed by the atmosphere.Assuming τ s;L ! 1 leads to 1 ´e´τ s;L « τ s;L and we recover the usual transfer function in the optically thin limit (Pierrehumbert 2011).

Role played by sensible heating
The formulation of the sensible heat flux given by Eqs. ( 54) and (57) closes the system formed by energy budget equations and allows us to determine the state of equilibrium of the system.It is convenient here to adopt the formalism of Wordsworth (2015) since we follow the same approach as this early work.
Article number, page 11 of 27 τ s;L " 10 ´2 τ s;L " 10 ´1 τ s;L " 10 0 Fig. 6.Atmospheric, dayside and nightside surface temperatures as functions of logarithm of control parameter L sen .Left: τ s;L " 0.01 (quasitransparent).Middle: τ s;L " 0.1 (small optical thickness).Right: τ s;L " 1 (Earth-like optical thickness).In all cases, the atmosphere is quasitransparent in the shortwave (τ s;S " 10 ´4), pure absorption is assumed (β L0 " β S0 " 1), A s " 0.2, and F ‹ " 1366 Wm ´2.The atmospheric temperature is designated by the solid black line, the dayside surface temperature by the dashed red line, the nightside surface temperature by the dotted blue line, and the black body equilibrium temperature (Eq.( 2)) by the dotted grey line.The value of L sen for a 1-bar CO 2 -dominated atmosphere is indicated by the green triangle.
Thus, introducing the normalised temperature, radiative flux, and sensible heat flux, , and Fsen " we transform Eqs. ( 46), (50), and (54) into Fsen " L sen F11{12 where the dimensionless parameter L sen is defined as The parameter L sen controls the intensity of sensible heating with respect to radiative heating on the dayside.The case L sen " 0 corresponds to pure radiative equilibrium (no turbulent heat exchanges), while L sen " 1 indicates a very strong turbulent mixing within the boundary layer, which leads to T a " T s;d .The system of Eqs.(60-62) is finally reduced to one single equation.By combining Eqs. ( 60) and ( 61) we first express the normalised stellar flux as a function of the normalised temperature, Substituted in Eqs. ( 61) and ( 62), this expression of F‹ yields ´Lsen which is the counterpart of Eq. ( 45) in Wordsworth (2015) with the heat engine approach proposed by Koll & Abbot (2016).This equation can be solved by using a combination of the dichotomy and secant methods (e.g.Press et al. 2007, p. 449).The obtained Td is then substituted in Eq. ( 64) and Eqs.(60-62) to compute successively F‹ , Fsen , the atmospheric temperature, and the dayside and nightside surface temperatures.Numerical calculations are performed using TRIP (Gastineau & Laskar 2011).
Figure 6 shows the evolution of the atmospheric temperature (T a ), and the dayside (T s;d ) and nightside (T s;n ) surface temperatures with the atmospheric opacity in the longwave -quantified by τ s;L -and the control parameter of sensible heating L sen .In this example, pure absorption is assumed (β L0 " β S0 " 1), the atmosphere is quasi-transparent in the shortwave (τ s;S " 10 ´4), and parameters are set to A s " 0.2 and F ‹ " 1366 W m ´2, which is the stellar flux received by the Earth.The value of L sen in the typical case of a 1-bar CO 2 -dominated atmosphere (green triangle) is estimated by taking p s " 1 bar, R s " 189 J kg ´1 K ´1, C p " 650 J kg ´1 K ´1, C D " 3.4 ˆ10 ´3, and ε sen " 1{2.
The figure highlights the role played by sensible heating.As L sen increases, the thermodynamic state of equilibrium switches from the purely radiative regime (L sen ! 1) to a strongly convective regime (L sen " 1), where atmospheric and dayside surface temperatures are homogenised by sensible exchanges.Both in the optically thin (τ s;L ! 1) and thick (τ s;L " 1) limits, we retrieve for Td a similar evolution with L sen as that shown by Fig. 7 of Wordsworth (2015), which indicates that the two different prescriptions used to estimate the typical horizontal wind speed (heat engine and scaling analysis using the thermodynamic equation) give comparable results.As expected, increasing the atmospheric opacity tends to accentuate the greenhouse effect, and thus to increase the nightside surface temperature.

Inclusion of large-scale advection
The last physical ingredient that has to be introduced in the model is a mechanism responsible for heat transport from the dayside to the nightside.As discussed in Sect.2.3, heat can be advected by large-scale atmospheric flows or transported by gravity waves (e.g.Koll & Abbot 2015, 2016), which propagate in stably stratified fluid layers (Gerkema & Zimmerman 2008).
Article number, page 12 of 27 P. Auclair-Desrotour and K. Heng: Atmospheric stability and collapse on tidally locked rocky planets In the present work, we assume that heat transport is due to atmospheric circulation solely and we focus on the regime of slow rotators ( LRo " 1), where mean flows are symmetric with respect to the axis connecting the stellar and anti-stellar points.From now on, we relax the uniform temperature approximation and consider that T a;d ‰ T a;n in the general case.First, we derive a scaling of the heat flux associated with advection F adv as a function of the dayside and nightside atmospheric temperatures by using the heat engine theory (Fig. 7).Second, we use this expression to derive and solve the equilibrium equation in the general case.

Radiative cooling
Large-scale circulation

Heating
T a;n T a;d Fig. 7.A diagram of large-scale stellar and anti-stellar atmospheric circulation as planetary heat engine.The heat engine is driven by dayside heating and nightside cooling to space.In this study, the atmospheric circulation is balancing Rayleigh drag, which causes frictional energy dissipation.

Scaling of the advected heat flux
To quantify the atmospheric heat transport, we proceed to a standard scale analysis (e.g.Pierrehumbert 2011, Sect.9.2.3, p. 610).
The heat per unit mass contained in a fluid parcel starting from the dayside is C p T a;d .When the parcel moves to the nightside, it decays to C p T a;n .The energy released by the parcel thus corresponds to C p pT a;d ´Ta;n q.Moreover, the mass flow that goes across the annulus separating the day-and nightside is proportional to the air column mass p s {g, to the size of the annulus 2πR p , and to a typical wind speed V adv .According to these considerations, the hemisphere-averaged heat flux due to atmospheric circulation may be written as where ε circ is an efficiency coefficient related to the threedimensional geometry of the flow.The wind speed V adv has now to be specified as a function of the atmospheric temperatures and of the system parameters.In order to clearly separate the mechanisms responsible for day-night heat transport on the one hand, and for sensible heat exchanges with the surface on the other hand, we choose to adopt a scaling law for the advection velocity different from that used to parametrise the convective cell (Eq.( 57)).We note however that the two mechanisms could be linked to each other by adopting the same scaling for both velocities.Koll & Komacek (2018) showed that the large-scale atmospheric circulation of hot Jupiters can be modelled as planetary heat engines and derived the typical wind speeds of mean flows in this framework.The reasoning applied to hot Jupiters in their work holds for rocky planets, and is very similar to that followed in Sect.4.1 to derive the expression of the sensible heating flux used in the model (Eq.( 62)).Here, the hot and cold reservoirs of the heat engine are the dayside and nightside parts of the atmosphere, and the temperature gradient driving the flow is thus horizontal instead of being vertical.In steady state, the work generated by differential heating and cooling is balanced by dissipation of kinetic energy, which determines the typical wind speed of the flow.In the present study, we assume that the atmospheric circulation is balancing Rayleigh drag in the terminator region.As a consequence, wind speeds scale as (e.g.Koll & Komacek 2018, Eq. ( 11)) where η adv " pT a;d ´Ta;n q {T a;d designates Carnot's thermodynamic efficiency for the large-scale circulation, ε adv the efficiency coefficient associated with the non-isentropic nature of the cycle, Q in the theoretical amount of power available to drive atmospheric motion, given by Eq. ( 58), and t drag the drag timescale, that is the timescale over which winds are linearly damped by Rayleigh drag.Substituting V adv by Eq. ( 67) in Eq. ( 66), introducing the total efficiency factor ε eff " ε circ ε adv , and extracting the dependence on the dayside and nightside atmospheric temperatures yields where we retrieve the scaling F adv {F eq 9 t rad {t adv given by early studies (e.g.Leconte et al. 2013).We remind ourselves here that t rad and t adv are the radiative and advective timescales defined by Eq. ( 13).The ratio t rad {t adv quantifies the ability of the stellar and anti-stellar circulation to transport heat from the dayside to the nightside.The case t rad {t adv " `8 corresponds to the regime treated by Wordsworth (2015) and Sect. 4 of the present work, where advection is so strong that the atmosphere is horizontally well mixed and its temperature homogenised.Conversely, if t rad {t adv ! 1, advection is not efficient in warming the nightside, which leads to strong day-night temperature gradients.We note that the drag timescale (t drag ) and the efficiency factor (ε eff ) appearing in Eq. ( 68) are unknown a priori.Both these parameters and the scaling law used for F adv should eventually be constrained with the help of more sophisticated models solving the coupled momentum and thermodynamic equations, such as GCMs typically.

Steady state equation in the general case
We can now write the system of equations describing the steady state in the general case, where radiative transfer, sensible heating and large-scale circulation are taken into account (Fig. 1).In (69) The equations Eqs. ( 46), (50), and (49) thus become where the normalised sensible heating and advection fluxes are expressed as Fadv In the preceding equation, we have introduced the nondimensional parameter which controls the intensity of heating due to atmospheric circulation with respect to radiative cooling.This parameter is complementary with L sen , defined by Eq. ( 63), which controls the intensity of sensible heating with respect to radiative cooling on dayside.which controls advective heat transport by the stellar and anti-stellar circulation, varies in the range log 10 pL adv q " ´2, . . ., 2 (from dark blue to green).The atmosphere is quasi-transparent in the shortwave (τ s;S " 10 ´4), pure absorption is assumed (β L0 " β S0 " 1), A s " 0.2, and F ‹ " 1366 Wm ´2, similarly as in Fig. 6.Here, the atmosphere is optically thin in the longwave: τ s;L " 0.1, which corresponds to the middle panels of Fig. 6.
The system of Eqs.(70-74) may be reduced to one single equation by proceeding similarly as in Sect.4.2.First, F‹ is expressed as a function of Ta and Td by substituting Fadv by Eq. ( 74) in Eq. ( 72), Second, Eqs. ( 70) and ( 71) are combined to eliminate Fsen , and the resulting equation is combined with Eq. ( 74) to eliminate Fadv .This yields the expression of Td as a function of Ta , Third, Td is substituted by Eq. ( 77) in Eq. ( 76), and both Td and F‹ are substituted by their expressions as functions of Ta in the normalised fluxes (Eqs.( 73) and ( 74)), which allows us to reduce the system of Eqs.(70-74) to a single equation of Ta .By using Eq. ( 70), we get This equation can be solved using a combination of the dichotomy and secant methods (Press et al. 2007), as previously done for Eq. ( 65).However, we note that the denominator of Eq. ( 77) is negative if Ta is less than a minimum value Ta;inf , meaning that Td is not defined in this range.This minimum value simply results from the fact that the advected heat flux cannot be greater than the heat absorbed on the dayside hemisphere.As the heat transport by large-scale advection increases as Ta decays (see Eq. ( 74)), a maximum of Fadv corresponds to a minimum of Ta , which is Ta;inf .We note that Ta;inf ! 1 if the efficiency of advection is low (L adv !1), while Ta;inf « 1 if the efficiency is high (L adv " 1).As a consequence, one has to preliminarily calculate Ta;inf by finding the zero of the denominator of Td in Eq. ( 77).As a second step, the solution of Eq. ( 78) is sought within the interval Ta;inf ă Ta ď 1 considering the fact that the nightside atmospheric temperature cannot be greater than the dayside atmospheric temperature in this model.Although there is no evidence for the existence and unicity of solutions from an analytical point of view, we verify them numerically in the studied range of the parameter space (see Appendix B).
However, as L adv Ñ `8, the dayside and nightside atmospheric temperatures tend to join together, meaning that very small temperature variations lead to huge variations of the flux associated with atmospheric circulation ( Fadv ).In other words, the day-night atmospheric temperature difference do not contain information any more in the strong circulation limit, and one should drop Eqs. ( 72) and ( 74), which corresponds to the case treated in Sect.4.2 and by Wordsworth (2015).This degeneracy Article number, page 14 of 27 P. Auclair-Desrotour and K. Heng: Atmospheric stability and collapse on tidally locked rocky planets L adv " 10 ´2 L adv " 10 0 L adv " 10 2 τ s;L " 0.1 T s;n T eq Fig. 9. Atmospheric, dayside and nightside surface temperatures as functions of control parameters L sen and L adv for optically thin (τ s;L " 0.1, top panels) and thick (τ s;L " 1.0, bottom panels) atmospheres in longwave.Left: L adv " 0.01 (weak circulation).Middle: L adv " 1.0 (mediumstrength circulation).Right: L adv " 100 (strong circulation).In all cases, the atmosphere is quasi-transparent in the shortwave (τ s;S " 10 ´4), pure absorption is assumed (β L0 " β S0 " 1), A s " 0.2, and F ‹ " 1366 Wm ´2, similarly as in Fig. 6.The dayside atmospheric (solid pink line) and surface (dashed red line) temperatures, as well as the nightside atmospheric (dotted cyan line) and surface (dashed blue line) temperatures are plotted as a function of the logarithm of L sen .The notation log 10 designates the decimal logarithm.has repercussions on the numerical solution of Eq. ( 78), which cannot be performed beyond a certain upper bound of L adv " 1 (typically L adv " 10 2 ´10 3 ).We also note that the logarithm of the distance to Ta;inf is a better appropriate coordinate than Ta itself in the root-finding procedure because of the strong variation of the function given by Eq. ( 78) in the vicinity of Ta;inf (see Fig.

B.1).
Figure 8 shows the evolution of the normalised temperature Td " T a;d {T s;d (Eq.( 59)) as a function of the parameters controlling sensible heating (L sen ) and large-scale heat transport (L adv ) in the optically thin case of Fig. 6 (τ s;L " 0.1, middle panels).We observe here the transition between the asymptotic limits of weak (L adv ! 1) and strong (L adv " 1) atmospheric circulations.In the regime of strong advection, the solution resembles that obtained by Wordsworth (2015) (see Fig. 7, top panel).Conversely, while the circulation weakens, a larger amount of heat remains on the planet's dayside.As a consequence, the temperature difference between the surface and the atmosphere decays.
In Fig. 9, the four temperatures of the system in steady state are plotted as functions of L sen for L adv " 0.01, 1, 100 in the optically thin (τ s;L " 0.1) and thick (τ s;L " 1.0) regimes.This figure shows that the nightside surface temperature strongly depends on L adv , while it is not very sensitive to L sen .As discussed above, increasing the efficiency of heat transport by the atmospheric circulation leads the nightside atmospheric and surface temperature to increase.In the regime of strong circulation (L adv " 1), the atmospheric nightside and dayside temperatures join together and we recover the results obtained in Fig. 6, where L adv " `8.The dependence of the nightside temperature on the atmospheric cir-culation is increased by greenhouse effect.In the optically thick case the temperature evolution between the cases L adv " 0.01 and L adv " 100 is larger than in the optically thick one.

Atmospheric stability
The analytical developments made in the preceding sections allow us to characterise the stability of the atmosphere, that is its ability to conserve its greenhouse gases and, more generally, to remain unchanged in composition.In the following, we use the model (Eq.( 78)) in parallel with that proposed by Wordsworth (2015) (Eq.( 45) of the article) to quantify the effects of the different processes participating to global heat redistribution on atmospheric stability.
The collapse is triggered by the condensation of one of the gases present in the atmospheric mixture, which occurs from the moment that the temperature of the air goes below the condensation temperature of the gas.The gas then condensates and forms an ice sheet at the planet's surface in the cold trap.This modifies the radiative properties of the atmosphere, as well as the atmospheric circulation, which are both coupled to the atmospheric gas mixture.In the standard case of greenhouse gases, collapse acts as a positive feedback by leading the atmosphere to cool down, and favouring thereby the condensation of other gases.
In the simplified approach of the present work, the complex dynamics of the collapse is not treated and the cold trap is assumed to be the whole nightside surface.We thus consider that, because of atmospheric circulation, a fluid parcel will end up cooling down at the nightside surface temperature, and we do Article number, page 15 of 27 not investigate the timescale of the process consequently.Under these hypotheses, the atmospheric stability is simply determined by the ratio between the nightside surface temperature and the lowest condensation temperature among the greenhouse gases present in the mixture, following early analytical studies (Wordsworth 2015;Koll & Abbot 2016).
In order to better understand the role played by the different involved processes, we use as a reference case the case studied by Wordsworth (2015), which is an Earth-sized planet hosting a CO 2 -dominated atmosphere.In this case, the atmosphere is assumed to be stable with respect to atmospheric collapse if T s;n ą T cond,CO 2 pp s q, the condensation temperature of CO 2 in K being given, below the triple point (p ă 5.18 ˆ10 5 Pa), by (Fanale et al. 1982;Wordsworth et al. 2010 To specify the relation between the short-and longwave optical thicknesses and surface pressure, we use a linear law of the form (Wordsworth 2015, Eq. ( 12)) where κ S and κ L are the supposed constant shortwave and longwave opacities, respectively.The longwave opacity is set to κ L " 5 ˆ10 ´5 m 2 kg ´1, which is the value used by Wordsworth (2015) in GCM simulations.This value corresponds to τ s;L « 1 for a 1-bar atmosphere.The shortwave opacity is set to κ S " 10 ´9 m 2 kg ´1, that is a value corresponding to the quasitransparent regime (κ S !κ L ).We note that pure absorption is assumed in the reference case (β S0 " β L0 " 1).Finally, the parameters characterising heat engines are specified from the observations made by Koll & Abbot (2016) and Koll & Komacek (2018) using GCM simulations.As mentioned above, the efficiency of the heat engine responsible for sensible heating was noticed to be twice smaller than the theoretical value (Koll & Abbot 2016).Thus, it is set to ε sen " 0.5.As the timescale of Rayleigh drag and the efficiency parameter of advection both scale the inter-hemispheric heat transport due to atmospheric circulation (see Eq. ( 75)), only one of these two parameters is needed to quantify the flux.Thus we arbitrarily fix t drag " 10 days, and choose the efficiency parameter ε adv so that L adv " 1, which corresponds to the asymptotic regime of the thermally homogenised atmosphere studied by Wordsworth (2015).The parameters values used in the reference case are gathered in Table 2.

Optical opacity favours collapse at low stellar fluxes
Our investigations start with the role played by longwave opacity on atmospheric collapse.By performing simulations with a GCM with correlated-k radiative transfer, Wordsworth (2015) computed a diagram of the atmospheric stability as a function of stellar flux and surface pressure in the reference case of the CO 2dominated atmosphere for an Earth-sized planet and a 10 M Csuper-Earth (Wordsworth 2015, Fig. 12).This diagram shows that the thresholds between stable and unstable regions of the parameter space predicted by GCM simulations and the analytic theory in the optically thin limit diverge from each other at low stellar fluxes, while they match well at high stellar fluxes.The divergence is patent around F ‹ « 500 W m ´2, where the critical pressure increases drastically as the stellar flux decays.Following Wordsworth (2015), we interpret this feature as a limitation of the optically thin approximation assumed in this early work.In order to better understand it, we compute the nightside temperature both with Wordsworth's model (designated by the acronym 'W2015'), and our model ('This work'), which captures in a simplified way the non-linear dependence of the atmospheric thickness on the surface pressure.Figure 10 shows the results of these calculations in the purely radiative case (L sen " 0, L adv " `8) and in the case where sensible heating is introduced (L sen ‰ 0), refereed to as 'General case'.As mentioned above, this later case corresponds to the asymptotic regime treated by Wordsworth (2015), where the day-night heat transport is efficient (L adv " 1), leading to T a;n « T a;d .In this regime, the values chosen for the timescale of Rayleigh drag (t drag ) and the efficiency parameter of advection (ε adv ) do not affect the planet's thermal state of equilibrium and the associated atmospheric stability from the moment that L adv " 1. Besides, the computed nightside temperature and atmospheric stability correspond to upper bounds since they can only decay if the efficiency of the day-night heat transport is modified3 .
The values of parameters used in the two cases are given by Table 2. Nightside temperatures as well as the corresponding stability diagrams are plotted as function of the stellar flux and surface pressure in logarithmic scales.For comparison, the incident stellar fluxes of TRAPPIST-1 f (Gillon et al. 2017 78) and Eq. ( 45) of Wordsworth (2015) with parameters values given by Table 2.
We first focus on stability diagrams (Fig. 10, bottom panels).Both in the purely radiative and in the general cases, the two models diverge at low stellar fluxes.While the evolution of the critical pressure with the stellar flux is well described by a power law in the optically thin limit (W2015), it undergoes a radical change of behaviour around F ‹ « F C when the nonlinear dependence of the atmospheric optical thickness on surface pressure is taken into account (this work).This change may be explained by comparing the scalings of the nightside temperature with surface pressure given by the two models (Fig. 10, top and middle panels).Below p s « 1 bar, the atmospheric opacity is linear with κ L , as predicted by Eq. ( 53) in the purely radiative case (see Table 1).Combining the fact that T s;n 9 T eq τ 1{4 s;L in the optically thin limit with the expression of the condensation temperature of CO 2 given by Eq. ( 79) yields the scaling where T 1mbar cond,CO 2 " T cond,CO 2 p100 Paq " 136.37 K is the condensation temperature of carbon dioxide at 1 mbar.We remark that this scaling law depends on the relation between the longwave optical thickness and pressure.In the case where this relation is a power law of the form τ s;L 9 p n s with n ą 0, the preceding scaling becomes p crit 9 F ´1{n ‹ . For instance, if the opacity increases due to pressure broadening, then n " 2 (Robinson & Catling 2012;Pierrehumbert 2011), meaning that the dependence of the critical pressure on the stellar flux is weaker in this case than in the studied case.
Beyond p s « 1 bar, the atmospheric opacity stops growing linearly with κ L .As a consequence, T s;n { pT s;n q rad.Ñ 0 as surface pressure increases, which leads to the observed dependence inversion.We note that the change of behaviour predicted by the model exceeds that derived from GCM simulations since the critical pressure tends to decay as surface pressure increases for p s Á 3 bar.This may be related to the limitations of the isothermal approximation made to derive radiative transfer functions in Sect.3, which is not appropriate to model the structure of thick atmospheric layers as discussed in the case of Venus.In spite of these limitations, the isothermal approximation turns out to be sufficient to capture the effect of large longwave opacities on the atmospheric stability at low stellar fluxes, and recover the associated behaviour of the critical pressure highlighted by GCM simulations.
Besides, we observe that the prediction of the model for atmospheric stability in Fig. 10 corresponds to an upper limit since stability diagrams were computed in the asymptotic regime previously treated by Wordsworth (2015), which maximises daynight heat transport (L adv " 1).Decreasing the efficiency parameter of advection would decrease the atmospheric stability and widen the region of the parameter space where collapse may occur as detailed in Sect.6.3.In the light of these considerations, the obtained stability diagrams (Fig. 10, bottom panels) predict atmospheric collapse on TRAPPIST-1 f, which agrees with the results obtained by Wolf (2017) 65) (L adv " `8) and Eq. ( 45) of Wordsworth (2015) with parameters values given by Table 2. Violet dots (top right panel) indicate the region of atmospheric stability computed by Wordsworth (2015) from GCM simulations (Fig. 12), while small blue dots indicate the region of collapse.

Dayside sensible heating increases stability
As a second step, we investigate how the dayside sensible heating affects the atmospheric stability.To do so, we place ourselves in the asymptotic regime studied by Wordsworth (2015), where the atmosphere has the same temperature on the dayside and nightside hemispheres.This allows us to ignore the effects of advection for the moment (L adv " `8) and to use the simplified equation given by Eq. ( 65) in the root-finding procedure instead of the general one (Eq.( 78)).We solve this equation for various values of the efficiency parameter controlling the intensity of sensible heating (ε sen ) and proceed similarly with Wordsworth's model (Wordsworth 2015, Eq. 45).
Figure 11 shows the stability diagrams derived from these calculations.The general case (i.e., with sensible heating) is compared with the purely radiative case treated in Sect.3.3 (Fig. 11, left column), and with the general case of W2015 (Fig. 11, right column).In addition, outcomes computed from 3D GCM simulations for the 1 M C -planet in this early study (Wordsworth 2015, Fig. 12, top panel) are included (Fig. 11, top right panel).Violet dots indicate simulations where the atmo-sphere remained stable, and blue dots simulations where collapse occurred.
As discussed in Sect.4, increasing the intensity of sensible heating acts to warm the atmosphere up.Consequently, the critical pressure decays as the efficiency of the dayside convective heat engine increases.However, this evolution is bounded.From the moment that L sen " 1, T a;d « T s;d and p crit reaches a minimum that is not sensitive to L sen , as observed for the nightside temperature in Figs. 6 and 9.One may show analytically that the maximum amplitude of critical pressure variations is of 0.6 decades in the studied case.
First, assuming that T a;d " T s;d in the strongly convective asymptotic regime (L sen " 1), we obtain from Eq. ( 64) that which yields the nightside temperature in this regime, Article number, page 18 of 27 P. Auclair-Desrotour and K. Heng: Atmospheric stability and collapse on tidally locked rocky planets If the atmosphere is transparent in the shortwave and optically thin in the longwave, then this expression simplifies to T s;n « T eq r2 p1 ´As q τ s;L s 1{4 .
The critical optical depth below which collapse occurs is determined by the equality T s;n pτ s;L q " T cond,CO 2 .This critical optical depth is denoted by τ rad s;L in the purely radiative regime and τ sen s;L in the regime dominated by sensible heating, and the corresponding critical pressures are denoted by p rad crit and p sen crit , respectively.Thus, assuming that T cond,CO 2 « T 1mbar cond,CO 2 , and using both Eqs. ( 53) and ( 85), we get the ratio between these two boundaries that is roughly 0.6 decades.This corresponds to what is observed for ε sen " 1 in Fig. 11 (top left panel).In W2015, no efficiency parameter such as ε sen intervene to modulate the sensible heat flux except a geometrical coefficient (χ), which is implicitly set to 1 here.This is why the critical pressure predicted by W2015 is the same in all panels of Fig. 11.Interestingly, setting ε sen " 1.0 reproduces the asymptotic behaviour of the critical pressure derived from W2015 at high stellar fluxes, although the two models use different approaches to scale the horizontal wind speed at planet's surface.Moreover, we note that the behaviour of the critical pressure predicted by Eq. ( 78) at low fluxes matches GCM simulations fairly well in this case.

Increasing the planet's size decreases stability
In the text above, the atmosphere is assumed to be homogeneous in temperature (T a;n " T a;d ).We now assume that the day-night atmospheric temperature distribution is controlled by the mechanism of heat transport by atmospheric circulation, which depends on dayside and nightside temperatures in return.This allows us to study the link between the intensity of advective heat transport -scaled by the efficiency parameter ε eff -and the atmospheric stability.Particularly, we aim to show how the 0-D theory, in spite of its limitations, may give us some insight to better understand the decrease of atmospheric stability that is observed as the planet size increases (e.g.Wordsworth 2015, Fig. 12).
In order to benchmark the predictions of the model against outcomes computed from GCM simulations, we reproduce the case treated by Wordsworth (2015) and thus consider two planets: the Earth-sized planet characterised by Table 2 and a super-Earth of 10 M C and 1.88 R C .As GCM simulations were performed using correlated-k distributions for radiative transfers in this early work, there is no evident value for the equivalent gray-gas longwave opacity.We thus set the absorptivity to κ L " 8.5 ˆ10 ´5 m 2 kg ´1 for a pedagogical purpose, this value leading the analytic model to match GCM simulations for a simple value of ε eff , as shown in the following.The equation of the steady state in the general case (Eq.( 78)) is then solved for various values of the efficiency parameters controlling sensible heating (ε sen ) and heat transport by atmospheric circulation (ε eff ).
Diagrams comparing the atmospheric stability of the two planets are plotted in Fig. 12 in the absence of sensible heating (ε sen " 0) and in the regime of strong dayside convection (ε sen " 1.0), for log 10 pε eff q " ´4, ´2.5, `8 (this later case is treated by solving Eq. ( 65) instead of Eq. ( 78)).Besides, we include the data computed from 3D GCM simulations (Wordsworth 2015, Fig. 12) in the figure (middle right panel).Violet dots designate simulations where the atmosphere remained stable in both cases, while red dots indicate those where the atmosphere of the 1 M Cplanet only remained stable.The figure shows that, in all cases, the atmospheric stability decays as the planet size increases.
We first consider the purely radiative regime (Fig. 12, top left panel).In this case, the observed difference between the Earthsized planet and the super-Earth comes from the dependence of the longwave atmospheric optical depth on the planet mass and radius that is hidden in the surface gravity (Eq.( 81)).As the planet size increases, the atmospheric optical depth decays, and so does the greenhouse effect.Consequently, the nightside surface temperature decays.The difference between the two planets does not change when sensible heating is included (top right panel) because this flux does not depend on the planet size at first order in the scaling given by Eq. ( 63) (there is a dependence hidden in Q in but it is weakened by the 1{3 exponent in the scaling of the horizontal wind speed).However, we observe that the gap of collapse pressure between the two planets predicted by GCM simulations is larger than that associated with the aforementioned dependence.Among the multiple reasons that may be invoked to explain this difference, we examine the hypothesis where it is related to a decrease of the efficiency of heat redistribution mechanisms through the example of large-scale advection.
Except in the asymptotic limit of strong horizontal mixing (L adv " 1), the intensity of the heat flux due to large-scale advection is modulated by the efficiency parameter ε eff , leading the atmospheric stability to decay with ε eff (Fig. 12, from top to bottom).Furthermore, in the derived scaling of the heat transport by stellar and anti-stellar atmospheric circulation, the heat flux is quantified by L adv 9 g ´1{2 R ´1 p (Eq. ( 75)) and thus decays as the planet size increases.As a consequence, the intensity of the heat flux from the dayside to the nightside decays as one switches from the Earth-sized planet to the super-Earth, which tends to widen the gap between the two planets (Fig. 12, bottom panels).This effect is stronger in the transition regime (L adv " 1) than in asymptotic ones, as shown by the evolution of the nightside and dayside temperatures ratio plotted in Fig. 8.
The case log 10 pε eff q " ´2.5 (i.e.ε eff « 3 ˆ10 ´3, middle panels) highlights the divergence between the analytic model and GCM simulations (Wordsworth 2015, Fig. 12).First, we note that the model tends to underestimate atmospheric stability at low stellar fluxes and high pressures, which is partly due to the simplified temperature profile used to derive radiative fluxes, as discussed in Sect.6.1.Second, the critical pressure predicted by the analytic model scales similarly with the stellar flux for both planets, while it is not the case in GCM results.This is obviously due to the numerous simplifications made in the present work, where the processes responsible for heat redistribution are described using rough scalings.Despite these limitations, the model approximately matches GCM simulations and captures the decrease of atmospheric stability caused by the increase of the planet size.It should be possible to improve the analytical theory with better scalings of the fluxes in future studies.

Role played by atmospheric structure
In our zero-dimensional approach, the implications of the atmospheric structure on the nightside temperature are ignored, given that atmospheric temperature profiles are reduced to dayside and nightside bulk atmospheric temperatures.This approximation Article number, page 19 of 27 12. Comparative stability diagrams for 1 M C and 10 M C planets as functions of logarithms of stellar flux (horizontal axis) and surface pressure (vertical axis).The efficiency parameter controlling advective heat transport ε eff take the values log 10 pε eff q " ´4, ´2.5, `8 (from bottom to top).Left column: Case without sensible heating (ε sen " 0).Right column: Case with efficient sensible heating (ε sen " 1.0).The stellar flux is normalised by the modern Earth's value F C " 1366 W m ´2, and the surface pressure is given in bar.The acronyms 'E' and 'SE' are employed for 'Earth' and 'Super-Earth', and designate the 1 M C and 10 M C planets, respectively.Colours indicate the stability of the steady state obtained in the two cases, similarly as in Fig. 10.These plots are obtained by solving Eq. ( 65) (L adv " `8, top panels) and Eq. ( 78) (other panels) with κ L " 8.5 ˆ10 ´5 m 2 kg ´1 and parameters values given by Table 2 for the Earth-sized planet.In the case of the 10 M C planet, the planet's radius and mass are set to M p " 10M C and R p " 1.88R C , respectively.Data computed by Wordsworth (2015) from GCM simulations (Fig. 12) are included in the middle right panel for comparison and indicate the region of atmospheric stability for the 1 M C planet (red dots) and both planets (violet dots), and the region of atmospheric collapse (small blue dots).
enabled us to derive radiative fluxes analytically, and to capture in a simplified way the non-linear dependence of the atmospheric opacity on the shortwave and longwave optical depths.In reality, the atmospheric structure is not the same on the dayside and on the nightside, and differs from the idealised isothermal temperature profile, as highlighted by analytical models (e.g.Robinson & Catling 2012;Koll & Abbot 2016) and GCM simulations (e.g.Leconte et al. 2013;Wordsworth 2015).Particularly, the evolution of the temperature profile with the stellar zenithal angle appears clearly in simulations performed using THOR (Fig. 2, left panels).The substellar region is characterised by a quasiadiabatic temperature gradient, which is due to the strong convection generated by the absorption of stellar heating near the ground.On the nightside, the structure inversion that may be observed indicates that the atmosphere is stable with respect to convection.
As mentioned above, the role played by the atmospheric structure was already investigated in the framework of the analytical theory.For instance Robinson & Catling (2012) proposed a one-dimensional radiative convective model to predict the atmospheric structure of a wide range of planets.Koll & Abbot (2016) refined this approach in the case of tidally locked rocky planets by including in a simplified way the effect of subsidence on the nightside temperature profile and surface temperature, which led them to develop a two-column radiative-convectivesubsiding (RCS) model.In this early work, the dayside temperature profile was supposed to be shaped by convection, and was thus specified in calculations.This assumption holds in the case of convective tropospheres, which corresponds to the regime of strongly irradiated atmospheres transparent in the shortwave.However, it may break down if the absorption near the ground is not large enough to generate strong convection.In view of these limitations, the prospect of extending the theory to optically thicker atmospheres in the shortwave motivates us to go past this assumption.
Therefore, although the core of the present study is the zerodimensional model, we make a first attempt to relax the specification of the dayside atmospheric structure in one-dimensional models for the purpose of this section.We aim to give here some insight about the transition from zero-dimensional to one-dimensional models where both the dayside and nightside temperature profiles are computed from atmospheric opacities.This leads us to develop a simplified two-column radiativeconvective-subsiding-upwelling (RCSU) model by introducing dayside upwelling winds in the RCS model of Koll & Abbot (2016).The RCSU model is detailed as a first step, and used as a second step to compute dayside and nightside atmospheric temperature profiles and study atmospheric stability against collapse.

A two-column RCSU model
To treat the coupling between the atmospheric structure and mean flows, one should integrate self-consistently the momentum, mass conservation, thermodynamic and radiative transfer equations, which is achieved by three-dimensional GCM Article number, page 20 of 27 but cannot be envisaged in the simplified framework of onedimensional models.We thus introduce vertical winds by using the thermodynamic equation solely, following the method by Koll & Abbot (2016).Written in pressure coordinates, this equation reads (Vallis 2006) the notation ∇ p designating the horizontal gradient over an isobar, V the horizontal velocity vector, ω " Dp Dt the pressure velocity (i.e. the vertical velocity in pressure coordinates; D Dt stands for the material derivative), and D the vertical diffusive power flux.To relate pressure to optical depth consistently with the early work by Koll & Abbot (2016), we use the standard power law where the exponent n specifies how optical depths increases with pressure.Introducing the dimensionless parameter β " R s { pnC p q and substituting the pressure coordinate by the optical depth in the longwave, Eq. ( 87) thus becomes In steady state, B t T " 0. Besides, we ignore the diffusion term B τ L D assuming that radiative transfers and advective heat transport predominates.In the general case, the horizontal advection term (V ¨∇p T ) is comparable to the vertical advection term (ωB τ L T ), and should be taken into account.However, we ignore it for simplicity following Koll & Abbot (2016), since including horizontal advection would require to specify ad hoc horizontal wind speeds and temperature distributions, which goes beyond the scope of a simplified two-column approach.Nevertheless, one may think about clever ways to introduce this component in the modelling in future studies.The above approximations yield where the complex atmospheric circulation is now reduced to one single quantity, ω, scaling the strength of vertical flows.We notice that ω ą 0 corresponds to subsidence, and ω ă 0 to updraft, since pressure and optical depth decay with altitude.
In this one-dimensional model, we consider a regime similar to that treated in Sect. 4 in the framework of the zerodimensional RC model.The heat transport from the dayside to the nightside is supposed to be efficient, and the stratosphere is consequently globally isothermal.However, the two cells introduced in the zero-dimensional model to characterise sensible heating and large-scale heat transport by atmospheric circulation separately are now encompassed into one single stellar and anti-stellar cell (see Fig. 5), so that updraft takes place within the dayside substellar region and subsidence on the nightside.Moreover, we assume uniform vertical profiles for pressure velocity both in the updraft and subsidence regions.The corresponding amplitudes are denoted by ω up and ω sub , respectively.The conservation of mass thus yields the relationship the notations A up and A sub designating the updraft and subsidence areas, respectively.As shown by snapshots of vertical wind speeds obtained from GCM simulations in typical cases (Fig. 2, right panels), there is an asymmetry between rising and sinking air motions.Particularly, air rises rapidly near the substellar point where heat absorption is the strongest, while it sinks slowly over a large area in cooler atmospheric regions.This asymmetry is quantified by the ratio A up {A sub , which is a fixed parameter of the model.As it characterises the strength of convection, the pressure velocity of rising air in the substellar region is related to the typical horizontal wind speed V sen introduced previously in the expression of the sensible heating flux (Eq.( 54)) and quantified using the heat engine theory (Eq.( 57)).
A simple scale analysis based on the conservation of mass yields With a similar analysis, we derive the typical time for a fluid parcel to rise up and to subside, Because of the asymmetry between rising and sinking motions, t sub " t up .For simplification, we assume in the following that A sub " 2πR 2 p and A sub {A up " t sub {t up " 10 and ignore the subsidence happening on the dayside (see Fig. 2, right panels).The steady state thermodynamic equation thus reads, for the dayside updraft region, and, for the dayside subsidence region, where t rad is the timescale of radiative cooling defined by Eq. ( 13).We note that the dayside net flux in the right-hand member of Eq. ( 94) is the sum of the longwave and shortwave components.For the later, we use the analytical solution given by Eq. ( 22), which is not limited to the case of isothermal profiles since it does not depend on temperature.
Combined with the day-and nightside net flux equations (Eq.( 16)), Eqs. ( 94) and (95) form the system of equations governing the global heat redistribution and dayside and nightside temperature profiles.For both the dayside and nightside, two boundary conditions are required to integrate the radiative transfer equation, and one to integrate the thermodynamic equation, meaning that six boundary conditions have to be specified.As may be noticed, the simplified dynamics adopted to derive Eqs. ( 94) and ( 95) generates a singularity at τ L " 0 in their Article number, page 21 of 27 A&A proofs: manuscript no.auclair-desrotour left-hand members.For this reason, upper boundary conditions cannot be applied at τ L " 0, and we have to set arbitrarily the optical depth of the top of the atmosphere τ top;L ą 0. We verify a posteriori that results do not vary with τ top;L from the moment that τ top;L !τ s;L .
We remark that Koll & Abbot (2016) could preliminarily calculate the level of the tropopause to which they applied nightside boundary conditions by assuming a convective dayside temperature profile.As the dayside temperature profile is not fixed in the present work, the pressure level of the nightside upper boundary cannot be determined in a similar way.Thus, upper boundary conditions have to be applied to the fixed level τ top;L .
On the dayside substellar area of rising air, we assume that there is no radiative source at infinity in the longwave, and subtract the amount of power available to drive atmospheric motion from the outgoing longwave radiation, assuming that dissipation occurs outside of the area.In the scaling of the typical wind speed V sen given by Eq. ( 57), the power available to drive atmospheric motion is quantified by ε 3 sen Q in , where Q in is the theoretical power defined by Eq. ( 58), and ε sen " 1{2 the efficiency coefficient evaluated by Koll & Abbot (2016) for the wind speed using GCM simulations 4 .We remind here that the day-night heat transport by atmospheric circulation is not taken into account in this section.The stratosphere is considered as globally isothermal, meaning that dayside and nightside atmospheric temperatures are equal at the upper boundary.The dayside surface temperature is defined by the blackbody radiation of the surface, meaning that T a;d " T s;d at surface.On the nightside, we assume the absence of downwelling longwave flux at infinity, as on the dayside.Besides, following Koll & Abbot (2016), we consider that the net flux in the longwave vanishes at the planet surface on the nightside -which is a consequence of stable stratification (no sensible heat flux) -and that the stratosphere is horizontally isothermal.Mathematically, these conditions read, for the dayside, and, for the nightside, The boundary value problem is solved over a fixed range of optical depths by means of a relaxation method.Starting from initial temperature profiles, the flux and temperature equations are integrated iteratively until convergence, an adaptive relaxation coefficient depending on the history of residuals being applied every iteration to stabilise the convergence process.At a given step, the dayside and nightside radiative transfer equations are integrated first with a finite difference scheme, by means of the shooting method (Press et al. 2007, Sect. 18.1) and Thomas's algorithm (Press et al. 2007, Sect. 2.4), respectively.Then, quantities parametrising the dynamics are calculated using auxiliary 4 We remark that including or not the term ε 3 sen Q in in the power budget has little impact on the obtained results given that the efficiency factor ε 3 sen is small.
equations.Finally, the dayside and nightside thermodynamic equations are integrated, which yields new temperature profiles.
The old temperature profiles are then incremented by the temperature difference between the old and new solutions weighted by the relaxation coefficient.This defines the temperature profiles for the next iteration, and so on.

Temperatures and fluxes profiles
We benchmark the two-column RCSU model against simulations performed by Wordsworth (2015), who used the threedimensional LMD GCM (Hourdin et al. 2006) with grey gas opacities.These simulations correspond to the reference case described by parameters values of Table 2 with n " 1 and p s " 0.01, 0.1, 1 bar (see Wordsworth 2015, Fig. 2, right panel).
We remind ourselves that the associated optical depths in the longwave are τ s;L " 0.01, 0.1, 1, respectively, owing to the chosen effective longwave opacity (κ L " 5 ˆ10 ´5 m 2 kg ´1).As may be noted, the substellar region is far smaller than the whole dayside hemisphere.Thus, the dayside atmospheric structure in this region is calculated by taking the value of the incident stellar flux at the substellar point, that is F ‹ .This differs both from Wordsworth (2015) and Koll & Abbot (2016), who calculated dayside averaged profiles (with an averaged incident flux of F ‹ {2), and leads consequently to higher dayside temperatures and stronger fluxes.
Temperatures and fluxes profiles computed using our RCSU model are plotted in Fig. 13, which is the counterpart of Fig. 2 (right panel) in the study by Wordsworth (2015).For comparison, we indicate in plots the radiative equilibrium temperature T rad , which corresponds to F ´L " 0, and is expressed in the absence of scattering as (e.g.Pierrehumbert 2011) T rad " T skin p1 `τL q 1{4 . (102) In this equation, T skin designates the skin temperature, that is temperature the outer regions of the atmosphere would have in the absence of in situ heating by stellar absorption (Pierrehumbert 2011, Sect.3.6, p.169).In the case where the atmosphere is transparent to incident stellar radiation, T skin " 2 ´1{4 p1 ´As q 1{4 T eq , which is typically the temperature of the stratosphere.
We recover in Fig. 13 the behaviour described by the zerodimensional RC model for surface temperatures, which are both increasing with the atmospheric optical thickness in the longwave because of greenhouse effect.Quantitatively, the predictions of the RCSU model for the nightside surface temperature match rather well the results obtained from GCM simulations, although the RCSU model tends to overestimate the nightside temperature with respect to the GCM in the case p s " 1 bar.The RCSU model also captures the variation of atmospheric structure between the dayside and the nightside.While convection on the dayside leads to a quasi-adiabatic temperature profile, the atmosphere is cooler and stably stratified on the nightside, which induces a smaller temperature gradient.Particularly, we retrieve here the nightside structure inversion observed in GCM simulations (e.g.Leconte et al. 2013;Koll & Abbot 2016) and due to subsidence, as discussed by Koll & Abbot (2016).
Because of the assumed vertically-uniform pressure velocity, both dayside and nightside atmospheric temperatures decay with pressure in order to compensate the increase of the singular term of the simplified thermodynamic equations (Eqs.( 94) and ( 95)).In reality, vertical winds vanish above the troposphere, and there .Temperature structure and fluxes profiles calculated with radiative-convective-subsiding-upwelling (RCSU) model in reference case for various surface pressures as functions of pressure levels in logarithmic scale.Left: p s " 0.01 bar and τ s;L " 0.01.Middle: p s " 0.1 bar and τ s;L " 0.1.Right: p s " 1 bar and τ s;L " 1. Dayside atmospheric (solid pink line), nightside atmospheric (dashed sky blue line), dayside surface (red square), nightside surface (blue square), and radiative equilibrium (dotted grey line) temperatures (K) are plotted in top panels.Dayside (shades of red) and nightside (shades of blue) radiative net, total, upwelling, and downwelling fluxes in the longwave, and net flux in the shortwave (dashed cyan line) are plotted in bottom panels.Fluxes are normalised by the black body equilibrium flux F eq " σ SB T 4 eq , where T eq is the black body equilibrium temperature defined by Eq. (2).Calculations are performed using parameters values given by Table 2 with n " 1, which corresponds to Fig. 2 (right panel) of Wordsworth (2015).Complementary parameters of the RCSU model are set to A sub " 2πR 2 p and A sub {A up " 10.
is no singularity.The tropopause approximately corresponds to the pressure level where T a;d and T a;n join together.Above this pressure level, the temperature profile should approach the radiative equilibrium temperature profile T rad (Eq.( 102)) instead of decaying, as shown by GCM simulations (e.g.Wordsworth 2015, Fig. 2, right panel).In spite of the unrealistic tendency observed in the stratosphere, the nightside temperature profile derived from the RCSU model in the troposphere is qualitatively similar to that obtained by Koll & Abbot (2016) with their radiative-convective-subsiding model (Fig. 4 of the article).Thus we may now compare the predictions of the two models regarding atmospheric stability and test thereby the RCSU model.

Atmospheric stability
Considering the reference case of the study, we run grid calculations using the RCSU model to compute the nightside temperature of the planet, and determine how the stability of the steady state against collapse evolves with the incident stellar flux and surface pressure.Consistent with previous sections, the atmospheric stability is simply obtained by comparing the nightside surface temperature to the condensation temperature of CO 2 , given by Eqs. ( 79) and ( 80).We note that the boundary value problem becomes difficult to solve accurately in the optically thick limit because tiny changes of the upper boundary condition (τ L " τ top;L ) in Schwarzschild equation (Eq.( 16)) substantially affect radiative fluxes at large τ L .As we start encountering convergence issues for τ s;L Á 10, we set the upper bound of the surface pressure range to log 10 pp s q " 0.5 (bar).
Figure 14 shows the obtained stability diagram, on which data computed by Wordsworth (2015) from GCM simulations are also plotted.This diagram may be compared to those calculated using the zero-dimensional RC model (Fig. 11) and to the results obtained by Koll & Abbot (2016) with their RCS model (Koll & Abbot 2016, Fig. 12, left panel).Although it includes the day-night difference of atmospheric structure, the RCSU model does not match GCM simulations, which were well approached by the RCS model of Koll & Abbot (2016).Particularly, the model does not capture the behaviour caused by the non-linear dependence of temperature on the atmospheric optical thickness, whereas this behaviour was captured by the 0-D approach, as shown by Fig. 11 (top right panel).The scaling of the collapse pressure given by the RCSU model actually corresponds to the purely radiative case of the 0-D model, that is p crit 9 F ´1 ‹ (Eq.( 82)).These limitations of the theory are apparently due to the oversimplified dynamics assumed to derive the differential equations of dayside and nightside temperatures (Eqs.( 94) and ( 95)).
Article number, page 23 of 27 Particularly, the pressure velocity (i.e. the velocity in pressure coordinates) is supposed uniform over the whole air column while it should tend to zero above the tropopause.As a consequence, the model does not describe properly the stratosphere.The inclusion of shortwave opacities in the two-column RCS approach thus still remain an open question for future studies.Moreover, temperature and winds profiles strongly vary with the stellar zenithal angle on the dayside, which cannot be properly modelled in the framework of a 1D theory.This encourages the development of 2D intermediate models treating the coupling between the large-scale atmospheric circulation and radiative transfers in a self-consistent way in the regime of slowly rotating planets.

Discussion and conclusions
In order to better understand the mechanism of atmospheric collapse, which is of fundamental importance to characterise surface conditions of tidally locked rocky planets, we developed a hierarchy of models building on pioneering works of the analytical theory in the spirit of Held (2005).We aimed thereby to enrich and to consolidate the theory by means of cross-comparisons between multiple approaches ranging from 0-D scalings to 3D GCM simulations.Models of increasing complexity were thus developed following a step-by-step method to treat the case of dry atmospheres hosted by slowly rotating rocky planets by including incrementally the couplings between radiative transfer, convection and large-scale circulation.In these models, radiative transfer is described using the twostream, dual-band and grey gas approximations, while the heat transport by atmospheric circulation is treated in the framework of heat engine theory following Koll & Abbot (2016).
We first revisited the early study by Wordsworth (2015), which is based on a 0-D approach.Assuming an isothermal atmosphere, we derived analytically the short-and longwave net flux profiles in presence of absorption and scattering, as well as the atmospheric, dayside and nightside surface temperatures characterising the steady state in the purely radiative case.This enabled us to quantify at first order how scattering and shortwave absorption modify the equilibrium temperatures and atmospheric stability against collapse by acting on greenhouse effect.
Typically, longwave and shortwave absorption both increase greenhouse effect, and thus the atmospheric stability, which is not the case of scattering.While scattering in the shortwave induces anti-greenhouse cooling and favours collapse, scattering in the longwave tends to stabilise the atmosphere by generating scattering greenhouse effect.The model captures these dependences in a simplified way.
It also captures the particular behaviour of the collapse pressure (or critical pressure) at low stellar fluxes due to the nonlinear dependence of the atmospheric thickness on the longwave opacity.Although highlighted by Wordsworth (2015) using GCM simulations, this behaviour was absent of the theory proposed by the author since the regime of optically thin atmospheres solely was considered in this early study.It results from the fact that the absorbed power ceases to grow linearly with the atmospheric optical depth at the planet's surface as one leaves the optically thin regime.The present study shows that, in spite of being a rough approximation of the atmospheric structure, the isothermal atmosphere assumption is sufficient to capture the induced scaling of collapse pressure with stellar flux that may be observed in GCM simulations.
As a second step, we included in the model the turbulent exchanges due to convection within the dayside boundary layer in the asymptotic regime where the atmosphere is globally isothermal (efficient large-scale heat transport).The sensible heat flux between the atmosphere and surface was scaled considering the atmosphere as a heat engine controlled by the atmospheric and dayside surface temperatures.In this framework, we derived the steady-state equation that governs the global heat redistribution, and solved it numerically to recover the early results obtained by Wordsworth (2015) and Koll & Abbot (2016) in the case of CO 2 -dominated atmospheres.Besides, we derived analytically a lower bound for the collapse pressure.This lower bound corresponds to the regime of strong convection, which maximises greenhouse effect.
As a third step, we ceased to treat the atmosphere as a globally isothermal layer and examined the case where the dominating heat transport mechanism is coupled with dayside and nightside temperatures.By way of an example, we supposed the heat to be advected from the dayside to the nightside by a stellar and anti-stellar atmospheric circulation, which is the typical dynamical regime of slow rotators.The associated heat flux was thus quantified from a scale analysis by treating the stellar and antistellar cell as a heat engine, as done for dayside convection.This simplified modelling accounts for the fact that the efficiency of heat transport by large-scale advection decays as the planet size increases.It allowed us to derive a generalised version of the steady state equation including the coupled effects of radiative transfer, sensible heating, and large-scale atmospheric circulation.
In GCM simulations, the atmospheric stability decays as the planet-size increases.Our generalised 0-D model approximately captures this behaviour and relate it to two conjugated effects: (i) for given gas opacity and surface pressure, the optical depth at the planet's surface decays as the surface gravity increases, which weakens greenhouse effect; and (ii) the ratio of the advective and radiative timescales decay, which decreases the hemisphere-averaged amount of heat transported from the dayside to the nightside.
The 0-D approach may certainly be improved in future studies by introducing refined scalings for heat fluxes and mean flows, which would better capture the dependence of the collapse pressure on the physical properties of the star-planet system.However, the model already allows us to characterise asymptotic Article number, page 24 of 27 P. Auclair-Desrotour and K. Heng: Atmospheric stability and collapse on tidally locked rocky planets regimes using a small number of dimensionless control parameters, such as -for instance -L sen and L adv , which compare dayside sensible heating and large-scale heat transport with radiative cooling, respectively.Furthermore, the fact that the root-finding procedure solving the steady state equation can be massively parallelised makes it possible to widely explore the parameter space.This clearly encourages the development of analytic 0-D models as a complement of studies based on GCM simulations.
As a last step, we investigated the transition between 0-D and 1D models by studying how the atmospheric circulation affects the dayside and nightside atmospheric structure.Starting from the two-column radiative-convective-subsiding (RCS) model proposed by Koll & Abbot (2016), we developed a radiative-convective-subsiding-upwelling (RCSU) model.In this approach, the dayside and nightside temperature profiles are modified by vertical winds.While upwelling winds due to convection in the substellar region drive the temperature gradient towards the adiabatic profile, subsidence on the nightside causes a structure inversion leading to stable stratification.In order to extend the theory to the case of non-transparent atmospheres in the shortwave, we let the dayside temperature profile be integrated simultaneously with the nightside temperature profile, instead of specifying its shape, which was done by Koll & Abbot (2016).
From a quantitative point of view, the RCSU model does not match well the behaviour of the collapse pressure obtained from GCM simulations, and rather predict a scaling similar to that given by the 0-D model in the purely radiative case.This apparently results from the simplifications made in the dynamics to include the effect of vertical winds.Particularly, in both subsidence and updraft areas, the vertical wind speed was supposed to be uniform over the whole air column, although the nature of the dynamics is not the same in tropospheric and stratospheric layers.Moreover, horizontal advection was ignored notwithstanding the fact that it is comparable to vertical advection in the general case.Future studies will have to remedy to both of these rough approximations to derive self-consistently realistic temperature profiles and collapse pressure at the same time with the 1D two-column approach, which still remains an open question.Nevertheless, in spite of its limitations, the RCSU model is a promising first attempt that allowed us to recover the features identified by early works, and particularly the nightside structure inversion previously captured by the RCS model (Koll & Abbot 2016).
For future works, it is crucial to improve the treatment of the dynamics in the analytic theory since it is a major aspect of the problem.This suggests to develop an intermediate class of 2D atmospheric models describing the complex coupling between mean flows and the thermodynamics in the asymptotic regime of slowly rotating tidally locked planets, the circulation being symmetric with respect to the star-planet axis in this regime.These simplified 2D GCMs would also be useful to quantify timescales associated with the atmospheric collapses, which cannot be achieved with steady states models, and requires a substantial computational effort when 3D GCMs are used.
In addition to the dynamics, one should evolve towards improved modellings of radiative transfers including wavelengthdependent opacities and accounting for the relationship between atmospheric absorption and stellar spectra, which is a blind spot of the dual-band approximation.This aspect is of fundamental importance when the stellar and planetary radiative fluxes overlap, that is typically in the case of temperate exoplanets orbiting cool dwarf stars.TRAPPIST-1 planets are a representative example of such a configuration.This question could be addressed for instance by introducing simplified dynamics in 1D radiative transfer codes.
Similarly as their predecessors, the simplified 0-D and 1D analytical models detailed in the present work are limited to dry thermodynamics and slow rotation, which corresponds to the simplest physics.The presence of water is however a feature of uppermost interest in the study of Earth-sized exoplanets located in the habitable zone of their host star because of its implications on climate and prebiotic chemistry.Thus, this feature should be introduced in the analytic theory in order to consolidate the conclusions of works essentially based on GCM simulations (e.g.Merlis & Schneider 2010;Leconte et al. 2013).
Particularly, the presence of a slab ocean induces an additional heat transport both by latent heat flux in the atmosphere and by heat diffusion in the ocean itself (Edson et al. 2011), and tends thereby to stabilise the atmosphere against collapse.Additionally, water vapour leads to the formation of clouds, which affects the atmospheric albedo, thermal structure, and general circulation in a complex way.To complete the theory, similar inclusions may be done for the heat transport by propagating gravity waves, and super-rotation in the regime of rapid rotators.Most of these aspects were discussed in early works (see e.g.Leconte et al. 2013).
Finally, we shall emphasise that steady states do not necessarily exist, and if so, are not necessarily stable.For example, Edson et al. (2011) showed evidence of the bistability characterising the dependence of the wind speeds on the planet spin by highlighting an abrupt transition between two regimes.This transition leads to a significant change of velocity.Regarding this type of questions, the analytic theory has a strong asset with respect to GCMs since it offers the opportunity to characterise formally the existence and stability of steady states.

Fig. 1 .
Fig. 1.Diagram of the two-layer radiative-convective-advective model detailed in the present study.The diagram shows the processes involved in dayside-nightside heat redistribution that are taken into account in the model.Radiative net fluxes (F ´" F Ò ´FÓ ) are upwards in the convention used, meaning that the shortwave net flux F ´S is negative.

Fig. 3 .
Fig.3.Net radiative fluxes in the short-and longwave (horizontal axis) as functions of optical depths (vertical axis) in case of isothermal atmosphere.Fluxes profiles are computed from the analytic formulae given by Eqs.(22) and (23) and normalised by the black body equilibrium flux F eq " σ SB T 4 eq (the equilibrium temperature T eq being given by Eq. (2)).Optical depths are normalised by their value at the planet's surface, namely τ s;S for the shortwave and τ s;L for the longwave.Shortwave (blue lines) and longwave (red lines) net fluxes are plotted for various optical thicknesses, the optical depths at the planet's surface taking the values τ s;S " 0.1, 1, 10 (from light to dark blue lines) and τ s;L " 0.1, 1, 10 (from light to dark red lines).Pure absorption is assumed (β L0 " β S0 " 1).The incident, atmospheric, and surface fluxes are arbitrarily set to F i " 4F eq , F a " 0.7F eq , and F s " 2F eq , respectively.
e 3 p A F b P 8 c 5 C M 7 X i 9 Z G c e t k s 1 D a z 0 Y 9 j i U s Y 5 X m u Y 0 S D n G M M n n f 4 B F P e D Z i 4 9 a 4 M + 7 7 q U Y u 0 y z i 2 z I e P g D p F p 1 z < / l a t e x i t > / l a t e x i t > R b x b R n z L n e k = < / l a t e x i t > Fig. 5.A diagram of the dayside atmospheric heat engine.The heat engine is driven by dayside heating and cooling to space.A fluid parcel works against friction in the boundary layer, which limits the speed of the flow.

Fig. 8 .
Fig.8.Normalised temperature Td " T a;d {T s;d as function of logarithm of parameter L sen for various values of L adv .The parameter L adv , which controls advective heat transport by the stellar and anti-stellar circulation, varies in the range log 10 pL adv q " ´2, . . ., 2 (from dark blue to green).The atmosphere is quasi-transparent in the shortwave (τ s;S " 10 ´4), pure absorption is assumed (β L0 " β S0 " 1), A s " 0.2, and F ‹ " 1366 Wm ´2, similarly as in Fig.6.Here, the atmosphere is optically thin in the longwave: τ s;L " 0.1, which corresponds to the middle panels of Fig.6.
Fig.10.Nightside temperature (K) predicted by analytical theory in this work (top panels) and inWordsworth (2015) (middle panel), and associated comparative stability diagrams (bottom panels) as functions of logarithms of stellar flux (horizontal axis) and surface pressure (vertical axis).Left: Purely radiative case as defined in Sect.3.3 (L sen " 0, L adv " `8).Right: General case (i.e., with sensible heating and heat transport by stellar and anti-stellar atmospheric circulation).The stellar flux is normalised by the modern Earth's value F C " 1366 W m ´2, and the surface pressure is given in bar.In stability diagrams, the acronyms 'TW' and 'W2015' are used for 'This work' and 'Wordsworth (2015)', respectively.Colours indicate the stability of the steady state predicted by the two models: both models predict stability (orange-red areas); the present model predicts stability while W2015 predicts collapse (orange areas); the present model predicts collapse while W2015 predicts stability (light blue areas); both models predict collapse (blue areas).The incident stellar fluxes of TRAPPIST-1 f and a hypothetical tidally locked Venus are designated by green triangles.These plots are obtained by solving Eq. (78) and Eq.(45) ofWordsworth (2015) with parameters values given by Table2.

Fig. 11 .
Fig. 11.Comparative stability diagrams as functions of logarithms of stellar flux (horizontal axis) and surface pressure (vertical axis) for various intensities of sensible heating, ε sen " 10 ´4, 10 ´2, 10 0 (from bottom to top).Left: Comparison between the general case (L sen ą 0) and the purely radiative case characterised in Sect.3.3 (L sen " 0).Right: Comparison between this work and the model developed by Wordsworth (2015) in the general case.The stellar flux is normalised by the modern Earth's value F C " 1366 W m ´2, and the surface pressure is given in bar.The acronyms 'TW', 'RC', and 'W2015' are used for 'This work' (with sensible heating), 'Radiative case' (without sensible heating), and 'Wordsworth (2015)', respectively.Colours indicate the stability of the steady state predicted by the two models, similarly as in Fig. 10.These plots are obtained by solving Eq. (65) (L adv " `8) and Eq.(45) ofWordsworth (2015) with parameters values given by Table2.Violet dots (top right panel) indicate the region of atmospheric stability computed byWordsworth (2015) from GCM simulations (Fig.12), while small blue dots indicate the region of collapse.
Fig.13.Temperature structure and fluxes profiles calculated with radiative-convective-subsiding-upwelling (RCSU) model in reference case for various surface pressures as functions of pressure levels in logarithmic scale.Left: p s " 0.01 bar and τ s;L " 0.01.Middle: p s " 0.1 bar and τ s;L " 0.1.Right: p s " 1 bar and τ s;L " 1. Dayside atmospheric (solid pink line), nightside atmospheric (dashed sky blue line), dayside surface (red square), nightside surface (blue square), and radiative equilibrium (dotted grey line) temperatures (K) are plotted in top panels.Dayside (shades of red) and nightside (shades of blue) radiative net, total, upwelling, and downwelling fluxes in the longwave, and net flux in the shortwave (dashed cyan line) are plotted in bottom panels.Fluxes are normalised by the black body equilibrium flux F eq " σ SB T 4 eq , where T eq is the black body equilibrium temperature defined by Eq. (2).Calculations are performed using parameters values given by Table2with n " 1, which corresponds to Fig.2(right panel) ofWordsworth (2015).Complementary parameters of the RCSU model are set to A sub " 2πR 2 p and A sub {A up " 10.
Fig. 14.Stability diagram computed with one-dimensional radiativeconvective-subsiding-upwelling (RCSU) model in reference case.The orange area indicate a stable atmosphere, and the blue area collapse in the RCSU model.Data computed byWordsworth (2015) from GCM simulations (Fig.12) are included for comparison.Violet dots indicate simulations where the atmosphere remained stable, while small blue dots indicate atmospheric collapse.Calculations are performed using parameters values given by Table2with n " 1, which corresponds to Fig.2(right panel) ofWordsworth (2015).Complementary parameters of the RCSU model are set to A sub " 2πR 2 p and A sub {A up " 10.

Table 2 .
Wordsworth 2015)ters used in the reference case of the present work (case treated byWordsworth 2015).The acronyms SW and LW are used in place of 'shortwave' and 'longwave', respectively.