Issue 
A&A
Volume 638, June 2020



Article Number  A77  
Number of page(s)  27  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202037513  
Published online  16 June 2020 
Atmospheric stability and collapse on tidally locked rocky planets
^{1}
ASD/IMCCE, CNRSUMR 8028, Observatoire de Paris, PSL, UPMC,
77 avenue DenfertRochereau,
75014
Paris, France
email: pierre.auclairdesrotour@obspm.fr
^{2}
University of Bern, Center for Space and Habitability,
Gesellschaftsstrasse 6,
3012
Bern,
Switzerland
^{3}
University of Warwick, Department of Physics, Astronomy & Astrophysics Group,
Coventry
CV4 7AL, UK
Received:
16
January
2020
Accepted:
13
April
2020
Context. Over large timescales, a terrestrial planet may be driven towards spinorbit synchronous rotation by tidal forces. In this particular configuration, the planet exhibits permanent dayside and nightside, which may induce strong daynight temperature gradients. The nightside temperature depends on the efficiency of the daynight heat redistribution and determines the stability of the atmosphere against collapse.
Aims. 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.
Methods. 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 largescale atmospheric circulation in the case of slowly rotating planets. There are two types of these models: a zerodimensional twolayer approach and a twocolumn radiativeconvectivesubsidingupwelling 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).
Results. 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 nonlinear 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.
Key words: astrobiology / planets and satellites: atmospheres / planets and satellites: terrestrial planets / methods: analytical / methods: numerical
© P. AuclairDesrotour and K. Heng 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 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 verylowmass stars (e.g. Payne & Lodato 2007; Raymond et al. 2007; Kopparapu et al. 2017).
Three out of seven of the Earthsized rocky planets hosted by the TRAPPIST1 ultracool 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 Earthmass planet candidates were detected around the Mdwarf Teegarden’s Star by the CARMENES spectrograph using radialvelocity measurements (Zechmeister et al. 2019), while an Earthsized planet orbiting in the habitable zone of the TOI700 red dwarf was discovered by the TESS observatory (Gilbert et al. 2020; Rodriguez et al. 2020). This planet (TOI700 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 subSolarmass stars, and particularly those orbiting verylowmass dwarfs such as TRAPPIST1, are expected to be circularised and tidally locked in spinorbit synchronous rotation (e.g. Kasting et al. 1993) or in a spinorbit 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 substellar and antistellar 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 better constrain the surface conditions of the Earthlike 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 atmosphericstability. 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 Mdwarfs and hosting CO_{2} /H_{2}O atmospheres, while Merlis & Schneider (2010) characterise the largescale circulation patterns of Earthlike planets tidally locked in spinorbitsynchronisation in various regimes.
A series of studies continued on this path, such as Heng et al. (2011a) and Heng & Vogt (2011), who characterise the atmosphericdynamics of a hypothetical tidally locked Earth and treat the case of the Gliese 581g superEarth as a scaledup 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 atmosphericstability, while Turbet et al. (2018) computed stability diagrams for TRAPPIST1 planets assuming N_{2} /CO_{2} atmosphericmixtures. 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 nextgeneration 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; PerezBecker & 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 nonlinear 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 temperaturepressure (T−P) profiles in the limit of the “doublegray” 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 nongrey radiative transfer due to spectral lines versus continua. Heng et al. (2014) generalised the work of Heng et al. (2012) to include nonisotropic/anisotropic scattering, and also demonstrated that the governing equations of these T−P profiles and the twostream 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 superEarths orbiting Mdwarf stars (e.g. Heng & Kopparla 2012). Laying the foundation of the 0D theory, Wordsworth (2015) proposed a twolayer model that selfconsistently 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 threedimensional 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 twocolumn radiativeconvective 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 radiativeconvectivesubsiding (RCS) model in the following. By performing a series of simulations with a GCM, Koll & Abbot (2016) showed that the largescale 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 nonlinear pressure dependence of the atmospheric optical thickness at low stellar fluxes, by dayside convection, and by largescale advection. Therefore, we adopt the 0D 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 antistellar 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 largescale 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 twocolumn radiativeconvectivesubsidingupwelling (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 0D 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 twostream, dualband, 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 antistellar circulation. In Sect. 6, the model is used to compute stability diagrams for CO_{2} dominated atmospheresderived 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 twocolumn RCSU model and discuss the role played by the atmospheric structure and, finally in Sect. 8 we summarise the conclusions and future works.
Fig. 1 Diagram of the twolayer radiativeconvectiveadvective model detailed in the present study. The diagram shows the processes involved in daysidenightside 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. 
2 Physical setup and main assumptions
We introduce in this section the main features of the 0D 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.
2.1 Tidal locking in 1:1 spinorbit resonance
We consider the simplified case of a tidallylocked 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 mainsequence 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), (1)
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 starplanet 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_{⋆}, (2)
the parameter σ_{SB} = 5.670367 × 10^{−8} W m^{−2} K^{−4} (Mohr et al. 2016) being the StefanBoltzmann constant, and the symbol ≡ meaning “defined by”. The above definition of r_{HZ} thus leads to . The stellar flux is also expressed as a function of the orbital radius a and the stellar bolometric luminosity L_{⋆}, (3)
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) (5)
with , the notations M_{⊙} and L_{⊙} referring to the Solar mass and luminosity, respectively. A linear regression of on the range − 1 ≤ μ ≤ 0 yields (6)
and thus . By comparing this scaling law to (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, (7)
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 mainsequence 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 (8)
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).
2.2 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 twostream 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 twostream approximation, we consider that the frequency spectra of the stellar and planetary radiative fluxes do not overlap, which is the socalled dualband 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 Sunlike 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 ultracool dwarf stars, such as TRAPPIST1, 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 tidallylocked exoplanets in their habitable zone than Sunlike stars, as stated in Sect. 2.1. We should thus bear in mind the limitations of the dualband 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 correspond 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 approximation 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 kcorrelated distribution method; Lacis & Oinas 1991) in GCMs used to study the heat redistribution on tidallylocked 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 shallowwater, twostream, and dualband approximations, one may consult the reference books by Seager (2010), Pierrehumbert (2010), and Heng (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 twocolumn 1D model, the selfconsistent 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 nonnegligible mathematical complications. Hence, following the early work by Wordsworth (2015), we opt for a zerodimensional twolayer model where the atmosphere is considered as isothermal across the vertical direction. However, we relax the wellmixed 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).
2.3 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 ofthree types: (i) radiative transfers along the vertical direction (shortwave and longwave absorption, radiation, and scattering); (ii) surfaceatmosphere turbulent exchanges in the dayside convective boundary layer; and (iii) daynight 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 largescale daynight heat transport turns out to be the most complex effect to include since there is no general theory of atmosphericheat 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 BuckinghamPi theorem (Buckingham 1914), which unravels the nondimensional parameters governing the system.
Among these parameters, we find the nondimensional Rossby deformation length, (9)
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, , 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 antistellar points, the daynight flow being driven by the balance between advection and pressuregradient acceleration (e.g. Leconte et al. 2013). In this stellar and antistellar circulation, highaltitude winds blow from the dayside to the nightside and lowaltitude 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 Opensource General Circulation Model THOR (Mendonça et al. 2016; Deitrick et al. 2019), which solves the general nonhydrostatic 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.1bar CO_{2} dominated atmosphere at T_{eq} = 279 K (Wordsworth 2015, Fig. 2), and a 0.5bar N_{2}dominated atmosphere at T_{eq} = 400 K (Koll & Abbot 2016, Fig. 4). Figure 2 shows instantaneous snapshots of the zonallyaveraged temperatureand 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 antistellar 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 daynight 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 (), the circulation regime changes and the atmospheric dynamics starts developing superrotation, that is strong eastward equatorial jets (Showman & Guillot 2002; Showman & Polvani 2011; Thrastarson & Cho 2010; Heng & Vogt 2011; Heng et al. 2011b; Leconte et al. 2013; Mendonça 2020). As shown by Showman & Polvani (2011), the emergence of superrotation 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 midlatitudes towards the equator. The latitude width of the equatorial band where superrotating 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 superrotation are caused by the latitudinal variations in radiative heating, the strength of equatorial jets increases with the daynight 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 the mean molecular weight of the atmosphere, we first introduce the specific gas constant and the pressure scale height, given as a function of temperature T by (10)
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 BruntVäisälä frequency N, which, in a dry stably stratified atmosphere, is given by (11)
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 (12)
In the present work, we only consider the heat advected by a daynight stellar and antistellar circulation, which corresponds to the asymptotic regime of slowly rotating planets (). We do not include other mechanisms of daynight 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 daynight heat flow is derived from a scale analysis in the heat engine framework (Koll & Abbot 2016) using for wind speeds the prescription given by 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) (13)
where p_{s} designates the atmospheric surface pressure and V_{adv} the typical speed of stellar and antistellar 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 daynight heat redistribution. We note however that this role may be nonnegligible. 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.
Fig. 2 Instantaneous snapshots of zonallyaveraged temperature (left panels) and vertical wind speed distributions (right panels) as functions of latitude and pressure with north pole at substellar point. Top: 0.1bar 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 antistellar 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. 
3 A twolayer 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.
3.1 Twostream analytic solution
In the twostream, dualband, and graygas 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 twostream Schwarzschild equations, which may be written as the system of first order partial differential equations (Heng 2017)
where ∂_{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 (β_{L 0} = 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 (16)
the notation 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 righthand 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 and are the two integration constants of the solution, and and the couplingcoefficients defined by (21)
The parameters characterise the coupling of radiative fluxes due to scattering. In the case of pure absorption (β_{L 0} = 1), there is no coupling ( and ). Conversely, in the case of pure scattering (β_{L0} = 0), the couplingis strong (). Similar expressions may be derived for the shortwave band by assuming T_{a} = 0 in Eqs. (17–20). The shortwave fluxes (F_{+S}, F_{−S}, F_{↑S}, and F_{↓S}) are parametrised by the shortwave optical depth τ_{S} (such that 0 ≤ τ_{S} ≤ τ_{s;S}, τ_{s;S} being the shortwave optical depth at the planet’s surface), scattering parameter β_{S 0}, coupling coefficients , and integration constants and .
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 reemitted 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, (22)
where we have introduced the atmospheric and surface fluxes, F_{a} and F_{s}, and the verticallyintegrated transmission functions (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.
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 (the equilibrium temperature T_{eq} being givenby 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. 
3.2 Surface and atmospheric energy budgets
The derivedanalytic 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 (25) (26)
which brings out the coefficient parameterising the surface and atmospheric energy balances,
with (similarly, in the short wavelength domain, we may define ). These coefficients are expressed as (34) (35) (36) (37) (38)
Thus, introducing the heat power per unit area due to surface turbulent exchanges on the dayside (F_{sen}) and daynight advection (F_{adv}) – specified in Sects. 4 and 5, respectively – the local energy balances of the planet’s surface and atmosphere are formulated as
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 fluid^{2}.
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 hemisphereaveraged 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)) (44)
Besides, the incident stellar flux is given as a function of the stellar zenithal angle ϕ by (45)
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 ().
3.3 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 daynight 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 (), 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, (50)
Solving this equation together with Eqs. (46) and (48) yields the surface and atmospheric equilibrium temperatures (51) (52) (53)
which capture the nonlinear 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 (). 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.3 W m^{−2} for the Solar bolometric flux).
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 quasitransparent 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 , 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 temperatureincreases 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 (β_{S 0} = 0.5, middle panels), it induces an antigreenhouse 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 (β_{L 0} = 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.
Optically thin and thick asymptotic limits in the pure radiative (dayside sensible heating and daynight heat transfers ignored) and pure absorption approximations (no scattering).
4 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 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).
4.1 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, 5th edition, Sect. 8.3.5, p. 268), as detailed in Appendix A. The hemisphereaveraged sensible heat flux is written as (e.g. Pierrehumbert 2011, Eq. (6.11), p. 396) (54)
where 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 () at a rate proportional to the temperature difference. Conversely, if the ground is cooler than the air, the atmosphere warms the ground (). We note that the convective cell vanishes in case of stable stratification, which implies V_{sen} = 0 and . 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 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 (1930) 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 (1930) constant , as (e.g. Esau 2004, Eq. (7)) (55)
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 hemisphereaveraged 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 surfaceatmosphere 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 (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) (56)
the factor being the atmosphere’s thermodynamic efficiency. The typical horizontal wind speed of the convective cell is thus expressed as (57)
where we have introduced the efficiency coefficient ε_{sen} accounting for the nonisentropic 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 (58)
with , 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 hemisphereaveraged power received by the planet on the dayside. The second factor is the verticallyintegrated 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).
Fig. 4 Atmospheric, dayside and nightside surface temperatures in pure radiative regime. Temperatures are plotted as functions of the logarithm of τ_{s;L} in the case of pure absorption (top panels), scattering in the shortwave (middle panels) and scattering in the longwave (bottom panels). Left: quasitransparence in the shortwave (τ_{s;S} =10^{−4}). Middle: small shortwave opacity (τ_{s;S} = 10^{−2}). Large shortwave opacity (τ_{s;S} = 10^{0}). The atmospheric temperature (solid black line), dayside surface temperature (dashed red line), nightside surface temperature (dashed blue line), and blackbody equilibrium temperature (dotted grey line) are plotted using Eqs. (51)–(53) and (2), respectively. Parameters values: A_{s} = 0.2, and F_{⋆} = 1366 W m^{−2}. 
Fig. 5 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. 
4.2 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. Thus, introducing the normalised temperature, radiative flux, and sensible heat flux, (59)
we transform Eqs. (46), (50), and (54) into
where the dimensionless parameter L_{sen} is defined as (63)
The parameter L_{sen} controls the intensity of sensible heating with respect to radiative heating on the dayside. The case L_{sen} = 0 correspondsto 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, (64)
Substituted in Eqs. (61) and (62), this expression of yields (65)
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 is then substituted in Eq. (64) and Eqs. (60)–(62) to compute successively , , 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 quasitransparent 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 1bar 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 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.
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 (Earthlike 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 W m^{−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 1bar CO_{2}dominated atmosphere is indicated by the green triangle. 
5 Inclusion of largescale 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 largescale atmospheric flows or transported by gravity waves (e.g. Koll & Abbot 2015, 2016), which propagate in stably stratified fluid layers (Gerkema & Zimmerman 2008). In the present work, we assume that heat transport is due to atmospheric circulation solely and we focus on the regime of slow rotators (), where mean flows are symmetric with respect to the axis connecting the stellar and antistellar 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.
Fig. 7 Diagram of largescale stellar and antistellar 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. 
5.1 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 . 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 windspeed V_{adv}. According to these considerations, the hemisphereaveraged heat flux due to atmospheric circulation may be written as (66)
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 daynight 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 largescale 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)) (67)
where designates Carnot’s thermodynamic efficiency for the largescale circulation, ε_{adv} the efficiency coefficient associated with the nonisentropic 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 (68)
where we retrieve the scaling given by earlier 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 antistellar circulation to transport heat from the dayside to the nightside. The case t_{rad} ∕t_{adv} = +∞ 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 daynight 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 should eventually be constrained with the help of more sophisticated models solving the coupled momentum and thermodynamic equations, such as GCMs typically.
5.2 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 largescale circulation are taken into account (Fig. 1). In accordance with the formalism employed in Sect. 4.2, we introduce the normalised temperature and heat flux, (69)
Equations (46), (50), and (49) thus become
where the normalised sensible heating and advection fluxes are expressed as
In the preceding equation, we have introduced the nondimensional parameter (75)
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.
The system of Eqs. (70)–(74) may be reduced to one single equation by proceeding similarly as in Sect. 4.2. First, is expressed as a function of and by substituting by Eq. (74) in Eq. (72), (76)
Second,Eqs. (70) and (71) are combined to eliminate , and the resulting equation is combined with Eq. (74) to eliminate . This yields the expression of as a function of , (77)
Third, is substituted by Eq. (77) in Eq. (76), and both and are substituted by their expressions as functions of in the normalised fluxes (Eqs. (73) and (74)), which allows us to reduce the system of Eqs. (70)–(74) to a single equation of . By using Eq. (70), we get (78)
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 is less than a minimum value , meaning that 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 largescale advection increases as decays (see Eq. (74)), a maximum of corresponds to a minimum of , which is . We note that if the efficiency of advection is low (L_{adv} ≪ 1), while if the efficiency is high (L_{adv} ≫ 1). As a consequence, one has to preliminarily calculate by finding the zero of the denominator of in Eq. (77). As a second step, the solution of Eq. (78) is sought within the interval 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} → +∞, 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 (). In other words, the daynight 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 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 is a better appropriate coordinate than itself in the rootfinding procedure because of the strong variation of the function given by Eq. (78) in the vicinity of (see Fig. B.1).
Figure 8 shows the evolution of the normalised temperature (Eq. (59)) as a function of the parameters controlling sensible heating (L_{sen}) and largescale 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) atmosphericcirculations. 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 atmosphericnightside and dayside temperatures join together and we recover the results obtained in Fig. 6, where L_{adv} = +∞. The dependence of the nightside temperature on the atmospheric circulation 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.
Fig. 8 Normalised temperature 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 antistellar circulation, varies in the range (from dark blue to green). 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 W m^{−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. 
6 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 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 Earthsized planet hosting a CO_{2}dominated atmosphere. In this case, the atmosphere is assumed to be stable with respect to atmospheric collapse if , 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; Wordsworth 2015) (79)
and, beyond the triple point (p ≥ 5.18 × 10^{5} Pa), by (80)
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)) (81)
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 1bar 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 interhemispheric 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.
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 quasitransparent in the shortwave (τ_{s;S} =10^{−4}), pure absorption is assumed (β_{L0} = β_{S0} = 1), A_{s} = 0.2, and F_{⋆} = 1366 W m^{−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. 
6.1 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 correlatedk 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_{2}dominated atmosphere for an Earthsized planet and a 10 M_{⊕}superEarth (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 nonlinear dependence of the atmospheric thickness onthe surface pressure. Figure 10 shows the results of these calculations in the purely radiative case (L_{sen} = 0, L_{adv} = +∞) 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 daynight 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 daynight heat transport is modified^{3}.
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 TRAPPIST1 f (Gillon et al. 2017) and a hypothetical tidally locked Venus are indicated in the bottom right panel of the figure (green triangles).
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_{⊕} 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 in the optically thin limit with the expression of the condensation temperature of CO_{2} given by Eq. (79) yields the scaling (82)
where K is the condensation temperatureof 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 with n >0, the preceding scaling becomes . 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, 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 TRAPPIST1 f, which agrees with the results obtained by Wolf (2017) using 3D GCM simulations.
Fig. 10 Nightside temperature (K) predicted by analytical theory in this work (top panels) and in Wordsworth (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} = +∞). Right: general case (i.e., with sensible heating and heat transport by stellar and antistellar atmospheric circulation). The stellar flux is normalised by the modern Earth’s value F_{⊕} = 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 (orangered 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 TRAPPIST1 f and a hypothetical tidally locked Venus are designated by green triangles. These plots are obtained by solving Eqs. (78) and (45) of Wordsworth (2015) with parameters values given by Table 2. 
Values of parameters used in the reference case of the present work (case treated by Wordsworth 2015).
6.2 Dayside sensible heating increases stability
As a secondstep, 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 temperatureon the dayside and nightside hemispheres. This allows us to ignore the effects of advection for the moment (L_{adv} = +∞) and to use the simplified equation given by Eq. (65) in the rootfinding 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 proceedsimilarly 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_{⊕}planet in this early study (Wordsworth 2015, Fig. 12, top panel) are included (Fig. 11, top right panel). Violet dots indicate simulations where the atmosphere 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 (83)
which yields the nightside temperature in this regime, (84)
If the atmosphere is transparent in the shortwave and optically thin in the longwave, then this expression simplifies to (85)
The critical optical depth below which collapse occurs is determined by the equality . This critical optical depth is denoted by in the purely radiative regime and in the regime dominated by sensible heating, and the corresponding critical pressures are denoted by and , respectively. Thus, assuming that , and using both Eqs. (53) and (85), we get the ratio between these two boundaries (86)
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.
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 thegeneral case. The stellar flux is normalised by the modern Earth’s value F_{⊕} = 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} = +∞) 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. 
6.3 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 daynight 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 0D 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 Earthsized planet characterised by Table 2 and a superEarth of 10 M_{⊕} and 1.88 R_{⊕}. As GCM simulations were performed using correlatedk distributions for radiative transfers in this early work, there is no evident value for the equivalent graygas 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 inthe 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 (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_{⊕}planet 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 superEarth 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 largescale advection.
Except in the asymptotic limit of strong horizontal mixing (L_{adv} ≫ 1), the intensity of the heat flux due to largescale 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 antistellar atmospheric circulation, the heat flux is quantified by (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 Earthsized planet to the superEarth, 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 (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.
Fig. 12 Comparative stability diagrams for 1 M_{⊕} and 10 M_{⊕} 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 (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_{⊕} = 1366 W m^{−2}, and the surface pressure is given in bar. The acronyms “E” and “SE” are employed for “Earth” and “SuperEarth”, and designate the 1 M_{⊕} and 10 M_{⊕} 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} = +∞, 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 Earthsized planet. In the case of the 10 M_{⊕} planet, the planet’s radius and mass are set to M_{p} = 10 M_{⊕} and R_{p}= 1.88 R_{⊕}, 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_{⊕} planet (red dots) and both planets (violet dots), and the region of atmospheric collapse (small blue dots). 
7 Role played by atmospheric structure
In our zerodimensional 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 enabled us to derive radiative fluxes analytically, and to capture in a simplified way the nonlinear 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 onedimensional 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 twocolumn radiativeconvectivesubsiding (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 onedimensional models for the purpose of this section. We aim to give here some insight about the transition from zerodimensional to onedimensional models where both the dayside and nightside temperature profiles are computed from atmospheric opacities. This leads us to develop a simplified twocolumn radiativeconvectivesubsidingupwelling (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.
7.1 A twocolumn RCSU model
To treat the coupling between the atmospheric structure and mean flows, one should integrate selfconsistently the momentum, mass conservation, thermodynamic and radiative transfer equations, which is achieved by threedimensional GCM 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) (87)
the notation∇_{p} designating the horizontal gradient over an isobar, V the horizontal velocity vector, the pressure velocity (i.e. the vertical velocity in pressure coordinates; stands for the material derivative), and 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 (88)
where the exponent n specifies how optical depths increases with pressure. Introducing the dimensionless parameter β ≡ and substituting the pressure coordinate by the optical depth in the longwave, Eq. (87) thus becomes (89)
In steady state, ∂_{t}T = 0. Besides, we ignore the diffusion term 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 (), 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 twocolumn approach. Nevertheless, one may think about clever ways to introduce this component in the modelling in future studies. The above approximations yield (90)
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 onedimensional 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 zerodimensional model to characterise sensible heating and largescale heat transport by atmospheric circulation separately are now encompassed into one single stellar and antistellar 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 (91)
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 (92)
With a similar analysis, we derive the typical time for a fluid parcel to rise up and to subside, (93)
Because of the asymmetry between rising and sinking motions, t_{sub} ≫ t_{up}. For simplification, we assume in the following that 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, (94)
and, for the dayside subsidence region, (95)
where t_{rad} is the timescale of radiative cooling defined by Eq. (13). We note that the dayside net flux in the righthand 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 lefthand 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 , 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 daynight heat transport by atmospheric circulation is not taken into account in this section. The stratosphereis 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 consequenceof stable stratification (no sensible heat flux) – and that the stratosphere is horizontally isothermal. Mathematically, these conditions read, for the dayside,
The boundary value problem is solved over a fixed range of optical depths by means of a relaxation method. Starting from initial temperatureprofiles, 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 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.
Fig. 13 Temperature structure and fluxes profiles calculated with radiativeconvectivesubsidingupwelling (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 , 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 and A_{sub}∕A_{up} = 10. 
7.2 Temperatures and fluxes profiles
We benchmark the twocolumn 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, which 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) (102)
In this equation, T_{skin} designates the skin temperature, which 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, , 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 quasiadiabatic 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 verticallyuniform 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 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 radiativeconvectivesubsiding 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.
7.3 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 withthe 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 (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 zerodimensional 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 daynight 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 nonlinear dependence of temperature on the atmospheric optical thickness, whereas this behaviour was captured by the 0D 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 0D model, that is (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)). 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 twocolumn 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 largescale atmospheric circulation and radiative transfers in a selfconsistent way in the regime of slowly rotating planets.
Fig. 14 Stability diagram computed with onedimensional radiativeconvectivesubsidingupwelling (RCSU) model in reference case. The orange area indicate a stable atmosphere, and the blue area collapse in the RCSU model. Data computed by Wordsworth (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 Table 2 with n = 1, which corresponds to Fig. 2 (right panel) of Wordsworth (2015). Complementary parameters of the RCSU model are set to and A_{sub}∕A_{up} = 10. 
8 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 analytical theory in the spirit of Held (2005). We aimed thereby to enrich and to consolidate the theory by means of crosscomparisons between multiple approaches ranging from 0D scalings to 3D GCM simulations. Models of increasing complexity were thus developed following a stepbystep method to treat the case of dry atmospheres hosted by slowly rotating rocky planets by including incrementally the couplings between radiative transfer, convection and largescale circulation. In these models, radiative transfer is described using the twostream, dualband 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 0D 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 atmosphericstability against collapse by acting on greenhouse effect. Typically, longwave and shortwave absorption both increase the greenhouse effect, and thus the atmospheric stability, which is not the case of scattering. While scattering in the shortwave induces antigreenhouse 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 from the theory proposed by the author since the regime of optically thin atmospheres was solely 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 largescale 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 steadystate 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 antistellar 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 largescale 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 largescale atmospheric circulation.
In GCM simulations, the atmospheric stability decays as the planetsize increases. Our generalised 0D 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 hemisphereaveraged amount of heat transported from the dayside to the nightside.
The 0D 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 starplanet system. However, the model already allows us to characterise asymptotic regimes using a small number of dimensionless control parameters, such as – for instance – L_{sen} and L_{adv}, which compare dayside sensible heating and largescale heat transport with radiative cooling, respectively. Furthermore, the fact that the rootfinding 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 0D models as a complement of studies based on GCM simulations.
As a last step, we investigated the transition between 0D and 1D models by studying how the atmospheric circulation affects the dayside and nightside atmospheric structure. Starting from the twocolumn radiativeconvectivesubsiding (RCS) model proposed byKoll & Abbot (2016), we developed a radiativeconvectivesubsidingupwelling (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 nontransparent 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 0D 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 selfconsistently realistic temperature profiles and collapse pressure at the same time with the 1D twocolumn 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 complexcoupling between mean flows and the thermodynamics in the asymptotic regime of slowly rotating tidally locked planets, the circulation being symmetric with respect to the starplanet 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 modelling of radiative transfers including wave lengthdependent opacities and accounting for the relationship between atmospheric absorption and stellar spectra, which is a blind spot of the dualband approximation. This aspect is of fundamental importance when the stellar and planetary radiative fluxes overlap, which is typically the case for temperate exoplanets orbiting cool dwarf stars. TRAPPIST1 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 0D 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 utmost interest in the study of Earthsized 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 superrotation 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.
Acknowledgements
The authors thank the reviewer, D. D. B. Koll, for helpful comments that contributed to improve the manuscript. They acknowledge financial support from the European Research Council via the Consolidator grant EXOKLEIN (grant number 771620). K.H. also acknowledges partial financial support from the National Swiss Foundation, the PlanetS National Center of Competence in Research, the Center for Space & Habitability and the MERAC Foundation. This research has made use of NASA’s Astrophysics Data System.
Appendix A Turbulent heat exchanges
At the planet surface, the heat flux resulting from turbulent friction is written as (A.1)
where is the density of the atmosphere at the planet’s surface, T′ and are fluctuations of temperature and radial velocity, respectively.
Fluctuations of the radial flow are induced by turbulence due to the interaction between the dominating horizontal flow and the planet surface. Thus, we can write the energy transport in Eq. (54) as , where V_{*} and T_{*} correspond to frictional velocity and temperature, respectively. Using the flux gradient theory with the mixing length hypothesis, these quantities are defined as (A.2)
where z is the altitude with respect to the surface, ℓ the turbulence mixing length, and and the average horizontal velocity and temperature, respectively.
The turbulence mixing length is calculated from the empirical scaling law of Blackadar (1962), (A.3)
the parameter being the von Kármán constant and ℓ_{0} the maximum attainable mixing length in the boundary layer. By substituting Eq. (A.3) into Eq. (A.2), and integrating the equations, we obtain
where z_{0} is the roughness height. Thus, considering that z is the thickness of the boundary layer and introducing the bulk drag coefficient (A.6)
we can write the vertical turbulent heat flux (from the surface to the atmosphere) as (A.7)
By assuming ℓ_{0} = +∞ (no fixed maximum for the turbulence mixing length), we recover the scaling law of the mixing length initially proposed by von Kármán on the basis of laboratory experiments, (von Kármán 1930), and the classical expression for the bulk drag coefficient given by Eq. (55). The parameters z and z_{0} are free parameters in the general case.They are set to z = 10 m and z_{0} = 1 cm in Wordsworth (2015).
Appendix B Solving the steady state equation in the general case
Equation (78) defines a rootfinding problem. This problem is formulated as (B.1)
where f is the function defined by the lefthand member of Eq. (78), that is
To solve Eq. (B.1), we first have to determine the range of where is defined, which corresponds to the range where the denominator of the expression given by Eq. (77) is strictly positive. This is a preliminary zerofinding problem, which is treated using the standard secant method (e.g. Press et al. 2007). We thus obtain the lower bound of the range, namely . Equation (B.1) is then solved by using a combination of the dichotomy and secant method within the interval .
We remark that f is subject to very sharp variations inthe vicinity of while it smoothly varies with in the rest of the interval. As a consequence, the coordinate is not well appropriate to the rootfinding procedure. We rather use the logarithm of the normalised distance, (B.2)
the notation log_{10} designating the decimal logarithm. The problem thus becomes , where . Finally, to deal with the possible issues resulting from the variations of f over several orders of magnitude, we apply the root finding method to the problem (B.3)
where ○ is the function composition operator and h the smoothing function defined as a function of the variable X by (B.4)
the notation sign referring to the sign function ( if X < 0 or X > 0, respectively, and ).
For an illustrative purpose, f is plotted in Fig. B.1 as a function of for differentvalues of the nondimensional parameters controlling atmospheric circulation, L_{sen} and L_{adv}. In these calculations, the atmosphere is optically thick in the longwave (τ_{s;L} = 10) and thin in the shortwave (τ_{s;S} = 0.1). The figure shows that there is one unique solution to the steady state equation in each of the treated cases. Besides, roots are found to be very close to with values ranging within the interval . This interval depends on the atmospheric opacity, and increases when the optical depth in the longwave decays.
Fig. B.1 Lefthand member of Eq. (78) as function of normalised temperature difference in logarithm scale for L_{sen}, L_{adv} = 0.01, 1, 100. The functionis denoted by f. The parameters L_{sen} and L_{adv} are thenondimensional control parameters of sensible heating and advected heat transport defined by Eqs. (63) and (75), respectively. Parameters values: β_{L 0} = β_{S0} = 1, τ_{s;L} = 10; τ_{s;S} = 0.1; F_{⋆} = 1366 W m^{−2}; A_{s} = 0.2. 
References
 Barclay, T., Pepper, J., & Quintana, E. V. 2018, ApJS, 239, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Barnes, R., Raymond, S. N., Jackson, B., & Greenberg, R. 2008, Astrobiology, 8, 557 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Bister, M., & Emanuel, K. A. 1998, Meteorol. Atmos. Phys., 65, 233 [Google Scholar]
 Blackadar, A. K. 1962, J. Geophys. Res., 67, 3095 [Google Scholar]
 Buckingham, E. 1914, Phys. Rev., 4, 345 [Google Scholar]
 Carone, L., Keppens, R., & Decin, L. 2016, MNRAS, 461, 1981 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1935, MNRAS, 96, 21 [NASA ADS] [Google Scholar]
 Correia, A. C. M., Boué, G., Laskar, J., & Rodríguez, A. 2014, A&A, 571, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cowan, N. B., & Agol, E. 2011, ApJ, 726, 82 [Google Scholar]
 Cumming, A., Butler, R. P., Marcy, G. W., et al. 2008, PASP, 120, 531 [CrossRef] [Google Scholar]
 Deitrick, R., Mendonça, J. M., Schroffenegger, U., et al. 2019, ApJS, submitted [arXiv:1911.13158] [Google Scholar]
 Deming, D., Seager, S., Winn, J., et al. 2009, PASP, 121, 952 [NASA ADS] [CrossRef] [Google Scholar]
 Dobrovolskis, A. R. 2009, Icarus, 204, 1 [Google Scholar]
 Edson, A., Lee, S., Bannon, P., Kasting, J. F., & Pollard, D. 2011, Icarus, 212, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Esau, I. 2004, Ann. Geophys., 22, 3353 [Google Scholar]
 Fanale, F. P., Salvail, J. R., Banerdt, W. B., & Saunders, R. S. 1982, Icarus, 50, 381 [Google Scholar]
 Gastineau, M., & Laskar, J. 2011, ACM Commun. Comput. Algebra, 44, 194 [Google Scholar]
 Gerkema, T., & Zimmerman, J. 2008, Lecture Notes, Royal NIOZ, Texel [Google Scholar]
 Gilbert, E. A., Barclay, T., Schlieder, J. E., et al. 2020, AAS J., submitted [arXiv:2001.00952] [Google Scholar]
 Gillon, M., Triaud, A. H. M. J., Demory, B.O., et al. 2017, Nature, 542, 456 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Goldreich, P., & Soter, S. 1966, Icarus, 5, 375 [Google Scholar]
 Grimm, S. L., Demory, B.O., Gillon, M., et al. 2018, A&A, 613, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Guillot, T. 2010, A&A, 520, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Held, I. M. 2005, BAMS, 86, 1609 [Google Scholar]
 Heng, K. 2017, Exoplanetary Atmospheres: Theoretical Concepts and Foundations (Princeton: Princeton University Press) [Google Scholar]
 Heng, K., & Kopparla, P. 2012, ApJ, 754, 60 [Google Scholar]
 Heng, K., & Vogt, S. S. 2011, MNRAS, 415, 2145 [NASA ADS] [CrossRef] [Google Scholar]
 Heng, K., Menou, K., & Phillipps, P. J. 2011a, MNRAS, 413, 2380 [NASA ADS] [CrossRef] [Google Scholar]
 Heng, K., Frierson, D. M. W., & Phillipps, P. J. 2011b, MNRAS, 418, 2669 [Google Scholar]
 Heng, K., Hayek, W., Pont, F., & Sing, D. K. 2012, MNRAS, 420, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Heng, K., Mendonça, J. M., & Lee, J.M. 2014, ApJS, 215, 4 [Google Scholar]
 Holton, J. R. 1973, Am. J. Phys., 41, 752 [Google Scholar]
 Hourdin, F., Musat, I., Bony, S., et al. 2006, Clim. Dyn., 27, 787 [Google Scholar]
 Howard, A. W. 2013, Science, 340, 572 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Joshi, M. 2003, Astrobiology, 3, 415 [NASA ADS] [CrossRef] [Google Scholar]
 Joshi, M. M., Haberle, R. M., & Reynolds, R. T. 1997, Icarus, 129, 450 [NASA ADS] [CrossRef] [Google Scholar]
 Kang, W., & Wordsworth, R. 2019, ApJ, 885, L18 [Google Scholar]
 Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icarus, 101, 108 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Koll, D. D. B. 2019, ApJ, submitted [arXiv:1907.13145] [Google Scholar]
 Koll, D. D. B., & Abbot, D. S. 2015, ApJ, 802, 21 [Google Scholar]
 Koll, D. D. B., & Abbot, D. S. 2016, ApJ, 825, 99 [Google Scholar]
 Koll, D. D. B., & Komacek, T. D. 2018, ApJ, 853, 133 [NASA ADS] [CrossRef] [Google Scholar]
 Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013, ApJ, 765, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Kopparapu, R. k., Wolf, E. T., Arney, G., et al. 2017, ApJ, 845, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Lacis, A. A., & Oinas, V. 1991, J. Geophys. Res., 96, 9027 [NASA ADS] [CrossRef] [Google Scholar]
 Leconte, J., Forget, F., Charnay, B., et al. 2013, A&A, 554, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lee, U., & Saio, H. 1997, ApJ, 491, 839 [NASA ADS] [CrossRef] [Google Scholar]
 Mendonça, J. M. 2020, MNRAS, 491, 1456 [CrossRef] [Google Scholar]
 Mendonça, J. M., Grimm, S. L., Grosheintz, L., & Heng, K. 2016, ApJ, 829, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Merlis, T. M., & Schneider, T. 2010, J. Adv. Model. Earth Syst., 2, 13 [Google Scholar]
 Mihalas, D. 1978, Stellar Atmospheres (New York: John Wiley & Sons) [Google Scholar]
 Mohandas, G., Pessah, M. E., & Heng, K. 2018, ApJ, 858, 1 [Google Scholar]
 Mohr, P. J., Newell, D. B., & Taylor, B. N. 2016, Rev. Mod. Phys., 88, 035009 [NASA ADS] [CrossRef] [Google Scholar]
 Mordasini, C., Alibert, Y., & Benz, W. 2009, A&A, 501, 1139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Parmentier, V., & Guillot, T. 2014, A&A, 562, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Payne, M. J., & Lodato, G. 2007, MNRAS, 381, 1597 [NASA ADS] [CrossRef] [Google Scholar]
 Peale, S. J. 1977, IAU Colloq. 28: Planetary Satellites, 87 [Google Scholar]
 PerezBecker, D., & Showman, A. P. 2013, ApJ, 776, 134 [NASA ADS] [CrossRef] [Google Scholar]
 Pierrehumbert, R. T. 2010, Principles of Planetary Climate (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
 Pierrehumbert, R. T. 2011, ApJ, 726, L8 [NASA ADS] [CrossRef] [Google Scholar]
 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2007, Numerical Recipes, 3rd edition: The Art of Scientific Computing (Cambridge: Cambridge University Press) [Google Scholar]
 Ragazzoni, R., Magrin, D., Rauer, H., et al. 2016, SPIE Conf. Ser., 9904, 990428 [NASA ADS] [Google Scholar]
 Raymond, S. N., Scalo, J., & Meadows, V. S. 2007, ApJ, 669, 606 [Google Scholar]
 Robinson, T. D., & Catling, D. C. 2012, ApJ, 757, 104 [NASA ADS] [CrossRef] [Google Scholar]
 Rodriguez, J. E., Vanderburg, A., Zieba, S., et al. 2020, AAS J., submitted [arXiv:2001.00954] [Google Scholar]
 Ronco, M. P., Guilera, O. M., & de Elía, G. C. 2017, MNRAS, 471, 2753 [NASA ADS] [CrossRef] [Google Scholar]
 Schroeder, D. V. 2000, An Introduction to Thermal Physics (Boston: Addison Wesley Longman) [Google Scholar]
 Seager, S. 2010, Exoplanet Atmospheres: Physical Processes (Princeton: Princeton University Press) [Google Scholar]
 Seager, S., & Deming, D. 2009, ApJ, 703, 1884 [Google Scholar]
 Seiff, A., Schofield, J. T., Kliore, A. J., et al. 1985, Adv. Space Res., 5, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Selsis, F., Wordsworth, R. D., & Forget, F. 2011, A&A, 532, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Showman, A. P., & Guillot, T. 2002, A&A, 385, 166 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Showman, A. P., & Polvani, L. M. 2011, ApJ, 738, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Suissa,G., Wolf, E. T., Kopparapu, R. k., et al. 2020, AJ, submitted [arXiv:2001.00955] [Google Scholar]
 Thrastarson, H. T., & Cho, J. Y.K. 2010, ApJ, 716, 144 [NASA ADS] [CrossRef] [Google Scholar]
 Turbet, M., Bolmont, E., Leconte, J., et al. 2018, A&A, 612, A86 [CrossRef] [EDP Sciences] [Google Scholar]
 Vallis, G. K. 2006, Atmospheric and Oceanic Fluid Dynamics (Cambridge: Cambridge Universsity Press), 770 [Google Scholar]
 von Kármán, T. 1930, Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, MathematischPhysikalische Klasse, 1930, 58 [Google Scholar]
 Wolf, E. T. 2017, ApJ, 839, L1 [Google Scholar]
 Wordsworth, R. 2015, ApJ, 806, 180 [CrossRef] [Google Scholar]
 Wordsworth, R. D., Forget, F., Selsis, F., et al. 2010, A&A, 522, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zechmeister, M., Dreizler, S., Ribas, I., et al. 2019, A&A, 627, A49 [NASA ADS] [EDP Sciences] [Google Scholar]
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).
Sensible heat is thus distinct from latent heat, which designates the energy associated with the change of phase of a condensable substance (Pierrehumbert 2011, Sect. 6.3).
If the value chosen for ε_{adv} were small enough, one would leave the asymptotic regime of efficient daynight heat transport and get in the transition regime showed by Figs. 8 and 9, which corresponds to a weaker heat transport, and thus leads to lower nightside temperatures and decreases the atmospheric stability.
All Tables
Optically thin and thick asymptotic limits in the pure radiative (dayside sensible heating and daynight heat transfers ignored) and pure absorption approximations (no scattering).
Values of parameters used in the reference case of the present work (case treated by Wordsworth 2015).
All Figures
Fig. 1 Diagram of the twolayer radiativeconvectiveadvective model detailed in the present study. The diagram shows the processes involved in daysidenightside 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. 

In the text 
Fig. 2 Instantaneous snapshots of zonallyaveraged temperature (left panels) and vertical wind speed distributions (right panels) as functions of latitude and pressure with north pole at substellar point. Top: 0.1bar 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 antistellar 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. 

In the text 
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 (the equilibrium temperature T_{eq} being givenby 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. 

In the text 
Fig. 4 Atmospheric, dayside and nightside surface temperatures in pure radiative regime. Temperatures are plotted as functions of the logarithm of τ_{s;L} in the case of pure absorption (top panels), scattering in the shortwave (middle panels) and scattering in the longwave (bottom panels). Left: quasitransparence in the shortwave (τ_{s;S} =10^{−4}). Middle: small shortwave opacity (τ_{s;S} = 10^{−2}). Large shortwave opacity (τ_{s;S} = 10^{0}). The atmospheric temperature (solid black line), dayside surface temperature (dashed red line), nightside surface temperature (dashed blue line), and blackbody equilibrium temperature (dotted grey line) are plotted using Eqs. (51)–(53) and (2), respectively. Parameters values: A_{s} = 0.2, and F_{⋆} = 1366 W m^{−2}. 

In the text 
Fig. 5 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. 

In the text 
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 (Earthlike 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 W m^{−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 1bar CO_{2}dominated atmosphere is indicated by the green triangle. 

In the text 
Fig. 7 Diagram of largescale stellar and antistellar 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. 

In the text 
Fig. 8 Normalised temperature 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 antistellar circulation, varies in the range (from dark blue to green). 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 W m^{−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. 

In the text 
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 quasitransparent in the shortwave (τ_{s;S} =10^{−4}), pure absorption is assumed (β_{L0} = β_{S0} = 1), A_{s} = 0.2, and F_{⋆} = 1366 W m^{−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. 

In the text 
Fig. 10 Nightside temperature (K) predicted by analytical theory in this work (top panels) and in Wordsworth (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} = +∞). Right: general case (i.e., with sensible heating and heat transport by stellar and antistellar atmospheric circulation). The stellar flux is normalised by the modern Earth’s value F_{⊕} = 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 (orangered 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 TRAPPIST1 f and a hypothetical tidally locked Venus are designated by green triangles. These plots are obtained by solving Eqs. (78) and (45) of Wordsworth (2015) with parameters values given by Table 2. 

In the text 
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 thegeneral case. The stellar flux is normalised by the modern Earth’s value F_{⊕} = 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} = +∞) 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. 

In the text 
Fig. 12 Comparative stability diagrams for 1 M_{⊕} and 10 M_{⊕} 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 (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_{⊕} = 1366 W m^{−2}, and the surface pressure is given in bar. The acronyms “E” and “SE” are employed for “Earth” and “SuperEarth”, and designate the 1 M_{⊕} and 10 M_{⊕} 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} = +∞, 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 Earthsized planet. In the case of the 10 M_{⊕} planet, the planet’s radius and mass are set to M_{p} = 10 M_{⊕} and R_{p}= 1.88 R_{⊕}, 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_{⊕} planet (red dots) and both planets (violet dots), and the region of atmospheric collapse (small blue dots). 

In the text 
Fig. 13 Temperature structure and fluxes profiles calculated with radiativeconvectivesubsidingupwelling (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 , 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 and A_{sub}∕A_{up} = 10. 

In the text 
Fig. 14 Stability diagram computed with onedimensional radiativeconvectivesubsidingupwelling (RCSU) model in reference case. The orange area indicate a stable atmosphere, and the blue area collapse in the RCSU model. Data computed by Wordsworth (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 Table 2 with n = 1, which corresponds to Fig. 2 (right panel) of Wordsworth (2015). Complementary parameters of the RCSU model are set to and A_{sub}∕A_{up} = 10. 

In the text 
Fig. B.1 Lefthand member of Eq. (78) as function of normalised temperature difference in logarithm scale for L_{sen}, L_{adv} = 0.01, 1, 100. The functionis denoted by f. The parameters L_{sen} and L_{adv} are thenondimensional control parameters of sensible heating and advected heat transport defined by Eqs. (63) and (75), respectively. Parameters values: β_{L 0} = β_{S0} = 1, τ_{s;L} = 10; τ_{s;S} = 0.1; F_{⋆} = 1366 W m^{−2}; A_{s} = 0.2. 

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