Multiple Climate States of Habitable Exoplanets: The Role of Obliquity and Irradiance

Stable, steady climate states on an Earth-size planet with no continents are determined as a function of the tilt of the planet’s rotation axis (obliquity) and stellar irradiance. Using a general circulation model of the atmosphere coupled to a slab ocean and a thermodynamic sea ice model, two states, the Aquaplanet and the Cryoplanet, are found for high and low stellar irradiance, respectively. In addition, four stable states with seasonally and perennially open water are discovered if comprehensively exploring a parameter space of obliquity from 0° to 90° and stellar irradiance from 70% to 135% of the present-day solar constant. Within 11% of today’s solar irradiance, we find a rich structure of stable states that extends the area of habitability considerably. For the same set of parameters, different stable states result if simulations are initialized from an aquaplanet or a cryoplanet state. This demonstrates the possibility of multiple equilibria, hysteresis, and potentially rapid climate change in response to small changes in the orbital parameters. The dynamics of the atmosphere of an aquaplanet or a cryoplanet state is investigated for similar values of obliquity and stellar irradiance. The atmospheric circulation substantially differs in the two states owing to the relative strength of the primary drivers of the meridional transport of heat and momentum. At 90° obliquity and present-day solar constant, the atmospheric dynamics of an Aquaplanet state and one with an equatorial ice cover is analyzed.


Introduction
Habitability on exoplanets for known forms of life depends on the presence of liquid water. This puts severe constraints on the climatic conditions that result from the orbital configuration of the exoplanet, in particular stellar irradiation, the orientation of the planet, rotation rates, atmospheric composition, and surface albedo, all of which are drivers of planetary surface temperature (Pierrehumbert 2010). Some of the largest climate variations on a planetary surface are generated by the obliquity or tilt of the planet's rotation axis, that is, its orientation relative to the orbit around the star. For instance, on Earth, obliquity is responsible not only for the seasonal cycle but also for one of the strongest climate variations, manifested as glacial-interglacial cycles with an amplitude of 3°C-4°C in global mean surface temperature over the past 5 million years (Masson-Delmotte et al. 2013). Changes in obliquity of the Earth during this time have a periodicity of 41,000 years (Berger & Loutre 1991). Although they occur in a narrow range of 22°-24. 5 during the past 10 million years, they act as an effective pacemaker of substantial variations of the terrestrial ice mass and hence global sea level. This is evidenced in reconstructions of global sea level based on marine sediments (Huybers 2007;Lisiecki & Raymo 2007) and on the temperature reconstructed from an Antarctic ice core (Jouzel et al. 2007). While in these records the obliquity cycle persisted throughout the past 5 million years, in the past 800,000 years a longer cycle of 100,000 years, coinciding with the periodic variations of the eccentricity of the Earth's orbit cycle, has become the dominant one. The importance of planetary obliquity has also been recognized for the habitability of exoplanets. Armstrong et al. (2014) have shown that oscillations in obliquity, caused by the gravitational interaction with nearby planets, can expand the habitable zone of a planet considerably.
Planetary mean surface temperature and the atmospheric composition are fundamental quantities determining habitable conditions on a planetary surface. Surface temperature is determined by the net energy flux from the star reaching the planet's surface, which primarily depends on the stellar irradiance, the planetary albedo, and the composition of the atmosphere (H 2 O, CO 2 , CH 4 ). Through Planck's law, even small changes in stellar irradiance may cause fundamental climatic changes on a planet. Likewise, small changes in greenhouse gas content of the atmosphere will produce substantial climate changes through positive feedback effects such as the water vapor and ice-albedo feedbacks (Pierrehumbert 2010).
The purpose of this paper is to comprehensively explore the impact of changes in obliquity and in the stellar constant on the equilibrium climate of a planet in order to characterize possible climate states and determine their physical conditions, including habitability. Our study is therefore not limited to conditions that may have existed on Earth, but more generally it contributes to the current and challenging research on the habitability of exoplanets.
The key determinant of planetary climate is the stellar constant S; in the case of the Earth under current conditions, S 1361 0 = Wm −2 . The first energy balance models demonstrated that reductions in S will cause transitions in the climatic state from a habitable planet to a completely frozen planet (Budyko 1969;Sellers 1969). This is due to the temperaturedependent ice-albedo feedback, which is able to generate multiple equilibrium states of the planet. This phenomenon was first found in simple energy balance models (Budyko 1969;Sellers 1969;North et al. 1981) and is also simulated in stateof-the-art dynamical climate models (Rose & Marshall 2009;Ferreira et al. 2011;Linsenmeier et al. 2015). The occurrence of different climate states is therefore independent of model complexity. This insight has been offered as a physical explanation for geological indications of temporarily very cold conditions on Earth, in particular low-latitude glaciation some 700 million years ago, referred to as the Snowball Earth hypothesis (Hoffmann & Schrag 2002). An alternative, but still controversial, explanation is that the Earth may have resided in a high-obliquity state caused by a giant impact at the early stages of Earth with a subsequent rapid transition to current low obliquity values some 500 million years ago (Williams 2008).
Qualitatively different instabilities can also occur at moderately increased values of stellar irradiation, where a runaway greenhouse effect can be triggered. Using a onedimensional model, Kasting (1988) found the occurrence of runaway at an irradiation of S 1.4 0 · for cloud-free conditions and S 4.8 0 · for full cloud cover. A much lower threshold of S 1.4 0 · for full cloud cover was reported by Popp et al. (2015). Recent simulations using comprehensive atmospheric models show stable, moist greenhouse states (e.g., Popp et al. 2016) or conditions conducive to a runaway at about S 1.2 0 · (Wolf & Toon 2015). The threshold for a runaway greenhouse continues to be debated (Goldblatt et al. 2013;Leconte et al. 2013). Taken together, even using comprehensive atmospheric models, it does not appear settled that a moderate increase of irradiation, such as by 40%, necessarily results in a runaway greenhouse.
The second determinant of planetary climate is the orientation of the rotation axis of a planet relative to its orbit. It governs the amount of energy that a planet receives from the star at a given latitude in the course of a planet's orbit around the star. It is therefore responsible for both planetary meridional temperature gradients and seasonality. For small angles of obliquity, the maximum input of annual mean energy is at the equator, while at the poles, the irradiance is strongly reduced. Because of these spatial differences, the atmosphere and ocean must transport energy from the equator to the poles by means of the general circulation (Lorenz 1967). The latitudinal gradient of annual mean energy supply to the planet changes sign at an obliquity of about 54°: at higher values of obliquity more energy is delivered to the polar latitudes than to the tropics (Jenkins 2000). In this situation, the atmosphere and ocean will transport energy in the equatorward to achieve a steady state. With such fundamental changes of the spatial pattern of planetary energy supply, large climate changes on a planet are expected if the obliquity changes.
Additionally, seasonality depends on obliquity ε. At zero tilt of the rotation axis, there are almost uniform conditions throughout the year, while at ε = 90°the amplitude of the seasonal cycle is at a maximum, with each latitude receiving no energy from the star during some period of the year. Obliquity is highly variable among planets, as is evident in our solar system, where they range from ε = 0.01°for Mercury to 97.86°f or Uranus, and 177.4°for Venus. Unfortunately, the tilt of rotating exoplanets is not among the information that is currently available from astronomical observations (Winn & Fabrycky 2015), which leaves a key determinant of habitability unknown. With the recent discovery of entire planetary systems within the habitable zone (Gillon et al. 2017), such information would be particularly informative, and some research has shown that this is in principle possible (Carter & Winn 2010;Zhu et al. 2014).
Going beyond earlier studies that used highly simplified energy balance models to illustrate possible climate states of exoplanets (Spiegel et al. 2009), three-dimensional atmospheric models (Shields et al. 2013(Shields et al. , 2014, and even comprehensive climate models for a small number of parameter values, such as stellar irradiation or obliquities (Ferreira et al. 2014;Hu & Yang 2014), here we employ a computationally efficient climate model of intermediate complexity, the Planet Simulator (PlaSim) developed by Lunkeit et al. (2011). This allows us to carry out extensive parameter exploration and a robust search for stable states. The model consists of a coarsely resolved general circulation model for the atmosphere coupled to a slab ocean and a thermodynamic sea ice model. This model represents the relevant dynamical processes of the atmospheric meridional heat transport in a realistic manner and allows seasonal or perennial sea ice cover if sea surface temperatures fall below the freezing point. Here we use the model in aquaplanet configuration, that is, without any terrestrial surface. In a previous study with this model at lower resolution, the effect of gradual changes of the solar constant at ε = 0°, 60°, and 90°was investigated (Linsenmeier et al. 2015). Their results suggest the existence of several climate states, but it is difficult to establish whether these are true steady states because the model was integrated for only 50 years.
A more complex global model with a dynamical ocean component in aquaplanet configuration was used by Ferreira et al. (2014), who considered the cases of ε = 54°and 90°and present-day insolation. Only two equilibrium climates were found: an Aquaplanet and a Snowball Earth, and they concluded that a stable state with equatorial ice was unlikely. Also, Earth configurations with modern and past continental outlines were used to assess the effect of high obliquity (Williams & Pollard 2003). They found a significant increase in snow accumulation during the winter for extreme seasonal variations. However, this has not lead to a net annual accumulation and subsequent glaciation, or even a transition to a Snowball Earth.
Today little is known about habitable conditions on Earthlike exoplanets if the full range of obliquity is considered. It is therefore timely to explore in a comprehensive manner the entire parameter space of obliquity and stellar irradiance in order to characterize the different steady-state climates that specific astronomical conditions are able to support. Our model permits a more realistic investigation of the meridional transport of heat carried by the atmospheric circulation in a specific climate state. A further question that will be addressed here concerns the circulation changes in an exoplanet atmosphere that result from transitions from one state to another due to slowly changing astronomical parameters, that is, state changes in the presence of multiple equilibria. This implies hysteresis behavior and hence evolution-dependent climatic conditions on an exoplanet.
The paper is organized as follows. Section 2 describes the climate model and the experimental setup. In Section 3 we explore the parameter space and describe the various steadystate climates that are found. Section 4 presents transient experiments and an analysis of the general circulation of the atmosphere with a focus on the differences between a Cryoplanet and an Aquaplanet under nearly identical astronomical conditions. The special case of 90°obliquity is considered in Section 5, and discussions and conclusions are presented in Section 6.

Climate Model and Experiments
In this study, the Planet Simulator, PlaSim, is used (Lunkeit et al. 2011). It is a climate model of intermediate complexity consisting of a general circulation model for the atmosphere (AGCM) coupled to a slab ocean and a thermodynamic sea ice model. The spectral resolution of the AGCM is T42, corresponding to a grid resolution of approximately 2.8°, and 10 vertical levels are chosen. The radiation formulation resolves clear and cloudy sky for both long-and shortwave radiation (Stephens 1978;Stephens et al. 1984), and the presence of clouds and their effect on radiation are parameterized (Slingo & Slingo 1991). The model is run under Earth-like conditions (planetary mass and radius, atmospheric composition, rotation rate about the planetary axis, and form of the ecliptic). Surface pressure is 1.01 10 Pa 5 · , and the concentration of CO 2 is kept constant at 360 ppm. The dynamical atmosphere is coupled to a slab ocean in which heat is transported diffusively and which, through its heat capacity, serves as the slow component of the coupled climate system. At the ocean-atmosphere interface, sea ice can be formed if the sea surface temperature falls below the freezing point. This generates a positive feedback to perturbations through changes of the albedo. The albedo of open ocean surface is 0.069, whereas the albedo of the sea ice surface is temperature dependent with a maximum value of 0.7 (Lunkeit et al. 2011). The minimum snow albedo is 0.4 for 0°C surface temperature and increases linearly up to 0.8 for surface temperatures smaller than −10°C. The model represents a compromise between complexity, permitting the analysis of atmospheric circulation and planetary surface changes, and efficiency, allowing comprehensive exploration of the parameter space, here obliquity and stellar constant.
Our model has reduced complexity in several aspects. It uses classical, simplified cloud parameterizations, and the sea ice model is thermodynamic. With a slab ocean of 50 m depth, ocean heat transport is represented as a diffusive flux. Hence, potential effects of changing ocean circulation cannot be captured with our model. Although this could have important consequences for the meridional distribution of heat and the presence of multiple equilibria, recent simulations suggest that the overall effect of ocean dynamics is limited and that a slab ocean is a good approximation for the major feedback processes involving the storage and distribution of heat by the ocean (Ferreira et al. 2014).
To focus on the fundamental physical processes, we select an aquaplanet configuration for all simulations (Kilic et al. 2017). Three different types of sensitivity simulations are performed while changing the two parameters obliquity and stellar constant: (1) simulations starting from a warm initial climate state (Aquaplanet), (2) simulations starting from a cold initial climate state (Cryoplanet), and (3) transient simulations continuously changing one of the two parameters and keeping the other constant. The first two types of simulations are performed for different combinations of parameters and run into an equilibrium state. Obliquity varies from 0°, in which no seasonal variations are present, to 90°, where the planet's rotation axis lies in the orbital plane. For irradiance, we chose a range of S S S 0 =˜· with S 1361 W m 0 2 = -, and S˜varies from 0.7 to 1.35. Simulations were integrated for 100 years, by which time a steady state was diagnosed based on the global mean surface temperature and global mean sea ice thickness. For the "warm start" setup, 367 simulations were performed, for the "cold start" 200. The third setup consists of four transient simulations where the obliquity is fixed and the stellar constant is continuously increased from S 0.52 = to 1.27 for ε = 47°, and S 0.87 = to 1.35 for ε = 23.44°with a rate of change in stellar constant of 0.01% per year. Using the same obliquities, the stellar constant is decreased starting from S 1.27 = and S 1.3 = , respectively, for the other two simulations. We note that in the range of S˜considered here, all simulations are run to steady state. Therefore, the presence of a runaway greenhouse can be excluded in our model simulations. For example, the maximum values of annual mean surface temperature and surface humidity are about 65°C and 149 g kg −1 , respectively, for S 1.35 = and ε = 15°. This is characteristic of a moist greenhouse (Kasting 1988).

A Diversity of Stable Climate States Depending on Obliquity and Irradiation
Here we present the first systematic exploration of parameter space of the two determinants of planetary climate, stellar irradiance and obliquity. We find several qualitatively different steady states and identify areas in parameter space where, depending on the initial conditions, two stable states exist ( Figure 1). In total we have identified six stable climate states, as summarized in Table 1. These states are classified into two general types: Aquaplanets with a dominant water surface, and Cryoplanets with a predominantly frozen surface. These two types are well known from previous studies using a hierarchy of models. In fact, the first latitude-dependent energy balance models have already simulated these two stable states (Budyko 1969;Sellers 1969), and many later modeling studies have confirmed their existence.
Depending on the specific values of obliquity and irradiance, several variants of these two types emerge (Table 1). This is relevant because it significantly enlarges the area of habitability in this parameter space, here defined as the presence of water during some period of the year. It is evident that irrespective of obliquity a Cryoplanet occurs for low stellar irradiances, while an Aquaplanet is simulated for high values of irradiance. The simulations show that the exact location where either of these two types prevails depends on the initial conditions. However, we find a rich structure of stable climate states for intermediate values of these two parameters (Figure 1).
Starting from Earth-like conditions, with the relative stellar irradiance S S S 1 0 = = and ε = 23.44°, and using a cryoplanet state as the initial condition, the planet remains completely frozen for these values of irradiance and obliquity (Figure 1(a)). If, on the other hand, the initial condition is an aquaplanet state, we find the Capped Aquaplanet (Figure 1(b)). This state has a permanent ice cover in the polar areas with an ice edge position that varies seasonally (Figure 2). The model therefore exhibits two stable states for the same parameter values at this location in parameter space. This is indicated by the overlaid state boundaries of Figure 1 Figure 1(a).
If obliquity or irradiance are slightly increased, the permanent ice cover in the polar areas becomes seasonal because of the increased energy supply to the high latitudes. This state is termed Near Aquaplanet. It occurs in a very narrow region in parameter space (Figure 1(b)) and can only be reached from aquaplanet initial conditions. If obliquity is further increased, the polar ice cover disappears completely, resulting in a permanent ocean surface, or an Aquaplanet. The Aquaplanet is a state that occupies a wide range in parameter space. In particular for S 1.04 > , it occurs for all obliquities (Figure 1(b)) if simulations are initialized from an Aquaplanet, and for S 1.3 > if they are started from a Cryoplanet (Figure 1(a)). This is the consequence of multiple equilibria that we find to be a robust feature of our simulations. The Aquaplanet survives down to about S 0.83 = , provided there is sufficient seasonal energy input, which is the case of ε = 60°.
Starting again from Earth-like conditions in the Capped Aquaplanet state, but now reducing stellar irradiance, a Cryoplanet results, which is the steady state with a complete ice cover throughout the year (Figure 2). For S 0.74 < this state occurs independent of obliquity. Increasing seasonality will supply sufficient energy to the high latitudes so that the ice cover opens seasonally if S 0.73 > . This state is called Near Cryoplanet. It features a band of ice extending from the tropics to the high latitudes with an open polar ocean during the summer season ( Figure 2).
Finally, if starting a simulation from cryoplanet conditions with, for example, S 1 = and ε = 90°, we find the Uncapped Cryoplanet, which is the sixth stable state (Figure 2). It consists of a low-latitude ice cover and a perennial ocean at high latitudes. This state was also reported by Linsenmeier et al. (2015) and occurred in the range of S 0.95 1.05 < < in their model. Here we confirm the existence of this sixth steady state in our model configuration with increased grid resolution and longer integration time. The seasonal and latitudinal characteristics of the six states  the Uncapped Cryoplanet coexists with the Aquaplanet state. The area of habitability is on the right side of the state boundary of the Cryoplanet, and it differs in size depending on the initial conditions; that is, the evolution history of the exoplanet is crucial in determining the current climate. The white star indicates the location of current Earth conditions in parameter space. Notes. Stable climate states are found for ranges in obliquity from 0°to 90°and in stellar irradiance from 70% to 135% of today's solar constant, and simulations are started from aquaplanet or cryoplanet initial conditions. The states can be classified in two types: Aquaplanets and Cryoplanets, with variants depending on the presence of seasonal ice cover or ocean surface. The notion "perennial" does not suggest an absence of a seasonal cycle but the presence of ice or ocean throughout the year. Sea surface temperature or ice thickness will still exhibit a seasonal cycle in these cases (see Figure 2). are displayed in Figure 2 using Hovmöller diagrams. The symmetry of the two types of Cryoplanets and Aquaplanets and their respective stable variants should be noted: each state of one type has its counterpart in the other type.
Our finding of the Uncapped Cryoplanet state is evidence against the supposition of Ferreira et al. (2014) that a stable state with an equatorial ice cover is unlikely to exist. It is interesting to hypothesize why they did not find this state with their more comprehensive model that has a dynamical ocean. Orbital parameters were almost identical to those we use here. With ε = 90°and S 1 = and starting from an Aquaplanet, their model remains in this stable state. The results using a slab ocean component instead are very similar, leading them to conclude that the ocean is simply storing and releasing heat as in a one-dimensional column (Ferreira et al. 2014). Our extensive sampling of the parameter space does not indicate the presence of an Uncapped Cryoplanet if the simulations start from an aquaplanet. This is consistent with the finding of Ferreira et al. (2014). We find the Uncapped Cryoplanet state only if we start from a Cryoplanet, as is indicated in Figure 1, where the complete area of this state is shared by the aquaplanet state. We therefore speculate that the more complete climate model of Ferreira et al. (2014) would also sustain a stable equatorial ice cover for high obliquity if it had started from cryoplanet conditions. Note also that Figure 1 exhibits a certain degree of symmetry. The orientation of the boundary between Aquaplanet states and Cryoplanet states in parameter space indicates a general trade-off between irradiance and seasonality. For higher seasonality, a smaller amount of total energy from the star is necessary for a transition from a Cryoplanet to an Aquaplanet. The location of this transition depends on the initial climate state. For high values of obliquity and low irradiance, a seasonal opening of the ice cover at high latitudes is simulated. For low values of obliquity, a stable ice cover in polar areas throughout the year can be sustained on an Aquaplanet, generating the Capped Aquaplanet state.

Transition from an Aquaplanet to a Cryoplanet: Impact on Atmospheric Circulation
Multiple equilibria are a robust feature of our simulations. They are present in a large range of relative stellar irradiance for S 0.83 1.3 < < and for all obliquities, as evidenced by the white boundaries in Figure 1(a). This can also be seen in Figure 3, where the two panels of Figure 1 are laid on top of each other.
Depending on the exact values of the two parameters, different pairs of stable states are found. For example, under Earth-like conditions of S 1 = and ε = 23.44°, we find both the Cryoplanet and the Capped Aquaplanet, while if we increase obliquities to 90°, four more pairs of steady states are visited (Figure 3). Except for the parameter values where the Cryoplanet is found as a single steady state, the entire space produces conditions that permit at least the seasonal presence of water and hence are favorable for life in some regions on the planet.
Some of the stable states can be accessed only from specific initial states. The Capped and Near Aquaplanets cannot be reached from a cryoplanet initial state, since for these parameter values the Cryoplanet is the other stable equilibrium. The Uncapped Cryoplanet, on the other hand, can only be reached when simulations start from a Cryoplanet. A complete realization of possible steady states therefore requires simulations starting from different initial states.
We now present two examples of transient simulations across state boundaries by slowly varying stellar irradiance  ( Figure 4). At ε = 47°there is a single state boundary that is crossed between S 0.9 = and 0.89 (solid line). A smaller jump in planetary annual mean temperature already occurs around S 1.06 = . For ε = 23.44°, the simulation successively crosses three state boundaries (dashed line) if the simulations start from the aquaplanet, whereas only one boundary is crossed in the path back along the lower branch of the hysteresis.
Across a state boundary, only small changes in the parameters suffice to generate large changes in all climate variables, in particular the atmospheric circulation of an exoplanet. For illustration, we present the two steady states Aquaplanet and Cryoplanet that result at ε = 47°and 48°and at S 0.89 = , both initiated with aquaplanet conditions. Figure 5 shows that in both states strong seasonal dynamics characterizes the atmospheric circulation. The meridional overturning cell located at the equator changes sign during the year in both cases, while the zonal wind maxima change hemispheres in the course of a year.
In order to further assess the impact of strong seasonality on the climate at this obliquity, we investigate the general circulation of the atmosphere and its underlying physical processes using monthly means. We focus on the month of January, where the impact of the different physical processes is more evident for both climate states ( Figure 6).
Although the two states receive the same stellar irradiance with almost identical seasonal distribution, the temperature distribution at the surface and throughout the atmosphere is very different (Figures 6(a), (b)). On the Cryoplanet, a large meridional temperature gradient develops with minimum temperatures in the winter hemisphere and maximum temperatures at the south pole in January. The pole-to-pole surface air temperature difference is about 104 K. On the Aquaplanet, on the other hand, a meridional temperature gradient is still present, but it is much weaker: the pole-to-pole surface air temperature difference here is only about 11 K. The atmospheric temperature structure generates the general circulation of the planetary atmosphere. In both states, the meridional circulation essentially consists of a single overturning cell near the equator (Figures 6(c), (d)) that changes sign seasonally ( Figure 5). This cell is centered around the equator with the ascending branch in the summer hemisphere. The meridional circulation transports both mass and energy (Figures 6(e), (f)). While the mass transport in the Aquaplanet state has a much larger intensity than in the Cryoplanet, the opposite is true for the meridional energy transport because of the much larger temperature gradients in the Cryoplanet state.
The meridional temperature gradients in the atmosphere also generate zonal circulations through the thermal wind relationship (Vallis 2006). Isotherms exhibit large slopes at around 30°S for the Cryoplanet. This is where a strong jet is generated (Figure 6(c)). This zonal circulation is an easterly jet, in contrast to Earth conditions, because here under high-obliquity conditions the polar regions are warmer than the equatorial regions during January. The Aquaplanet state shows a much more Earth-like zonal circulation with two weak westerly jets in the subtropics (Figure 6(d)).
The sharp equatorward temperature decrease around 30°S in the summer hemisphere of the Cryoplanet results from the strong gradient in stellar irradiation. This generates a zone of high baroclinicity around 30-60°S that can be maintained because the surface is uniformly covered by sea ice. The importance of this baroclinic behavior is also evident in the Lorenz energy cycle (Lorenz 1955(Lorenz , 1967. This is shown in Figure 7, where the eddy conversion (CE) is significantly higher than the zonal mean conversion (CZ), although eddy available potential energy (AE) is dissipated by a greater loss of heat at the surface by melting a layer of the sea ice surface around 60-90°S. For a slightly larger obliquity, an Aquaplanet is simulated with a meridional temperature gradient that is much smaller, because in the absence of any sea ice surface the albedo is reduced and stellar energy is absorbed by the slab ocean, which generates much higher surface temperatures that persist in the winter hemisphere. This leads to a significantly less baroclinic atmosphere. Energy fluxes under such conditions are quantified in the Lorenz energy cycle (Figure 7), where all energy conversions are rather small compared to the Cryoplanet state.
To further examine reasons for the different characteristics of the meridional circulation and the temperature pattern, we analyze the meridional energy transport in the atmosphere by calculating the moist static energy where v is the meridional velocity and m is the MSE given by m c T g z L q p e = + + · · · with c 1004 p = J kg −1 K −1 the specific heat capacity of dry air, L 2.5 10 e 6 = · J kg −1 the latent heat of vaporization, T the temperature, z the height, and q the specific humidity. The overbars and primes denote the time mean and its deviation, so the two terms on the right-hand side of the equation quantify the transport by the mean flow and by the transient eddies.
The MSE flux around the equator is dominated by the transport of the mean flow for both states (Figures 6(e), (f)). This MSE flux is related to thermally direct cells (Figures 6(c), (d)). However, the driving mechanisms are different in the two states: under cryoplanet conditions, only radiative heating around the equator is able to drive the thermally direct cell. Latent heating is almost absent because there is no open water. For the aquaplanet state, both radiative and latent heating are the main drivers of this thermally direct cell (Figure 8). They are both stronger due to the presence of a warm ocean surface that acts as a source of heat and moisture.
The part of the MSE flux (Figures 6(e), (f)) related to eddies behaves differently for the two states. The Cryoplanet shows a peak contribution of eddy flux around 50°S, which overlaps with parts of a thermally direct cell (Figure 6(c)). This is in contrast to the Aquaplanet state, where we find a peak contribution around 30°N, which coincides with a thermally indirect cell (Figure 6(d)).
To assess this difference, we consider the zonal mean pattern of meridional eddy heat and momentum fluxes (Figure 9). For the cryoplanet state, the maximum eddy heat flux occurs in the lower troposphere at about 55-60°S (Figure 9(a)) with convergence of eddy heat flux north and divergence south of this maximum. Therefore, the eddy heat transport is northward and tends to reduce the pole-to-equator mean temperature gradient. Furthermore, the gradient of eddy momentum flux in the range 35-65°S is negative (Figure 9(c)) and thus illustrates a convergence of the eddy momentum flux. This convergence decreases with height for the lower troposphere up to 500 hPa, which implies a clockwise circulation (Figure 9(c)). Thus, a thermally indirect meridional cell in the southern hemisphere ( Figure 6(c)) results, which turns clockwise as the thermal direct cell between 30°S and 30°N. For the Aquaplanet, the maximum eddy heat flux is at about 30-40°N (Figure 9(b)), leading again to a northward heat transport. The gradient of the momentum flux (Figure 9(d)) is negative, showing that momentum flux converges, and this convergence increases with height in the lower troposphere up to 400 hPa around 30°-40°, resulting in a counterclockwise, thermally indirect mean meridional cell.
In summary, the main difference between the two states is that the dominant meridional circulation cell is a thermally direct cell, driven by radiative and latent heating under aquaplanet conditions (similar to Earth), whereas the dominant cell of the cryoplanet state is driven by a combination of radiative heating and eddy heat and momentum fluxes, and thus is a superposition of a thermally direct and a thermally indirect cell.

Two Stable States at 90°Obliquity
The stable state that features a permanent ice cover from about 30°S to 30°N and a permanent open ocean toward the poles was identified for S 1 = and ε = 90°. The Uncapped Cryoplanet occurs only when simulations are started from a cryoplanet state (Figures 1 and 3). This state does not occur as a single equilibrium but is always paired with the aquaplanet state. Figure 10 compares these two stable states at ε = 90°for the month of January. Although stellar radiation is identical, the strong ice-albedo feedback causes very cold temperatures above the equatorial ice cover. The presence of the ocean as an efficient heat storage maintains warm temperatures at high latitudes even when the stellar irradiation is zero during the winter season (Figure 10(a)). One overturning cell is located at the ice edge in the summer hemisphere and supports a subsidence at the equator that may stabilize the ice cover ( Figure 10(c)). The strong meridional temperature gradients throughout the atmospheric column at 30°S generate a strong easterly jet. Note that the meridional ocean heat flux under the ice cover vanishes (Figure 10(e)).
At the same parameter values, the Aquaplanet exhibits very warm surface temperatures at all latitudes (Figure 10(b)), with ocean heat storage maintaining ice-free conditions during the winter season. In the summer hemisphere, large meridional temperature gradients develop a strong easterly jet and a strong meridional overturning circulation, which is the primary mechanism for heat transport convergence (Figure 10(d)). The structures of the zonal wind and the circulation cell are remarkably similar in the two states, although the Aquaplanet's meridional circulation is much stronger. This is driven by the latent heat release over the open ocean surface, as was already seen in Figure 8(d). Related to the different strengths of the circulation and jets, an important deviation between the two states is illustrated in the meridional energy flux, as for latitudes south of 30°S the transient part dominates for the Uncapped Cryoplanet state due to increased baroclinicity (increased temperature gradient). This is in contrast to the aquaplanet state, where the mean circulation part of the energy flux drives a poleward flux between 30 and 50°S.

Discussion and Conclusions
Stellar irradiance and obliquity are key determinants of the surface climate conditions of a planet. We have demonstrated that a richer variety of stable states can be simulated by a coupled atmosphere-slab ocean-sea ice climate model than just the two classical types of an Aquaplanet and a Cryoplanet with perennial open water and sea ice surface, respectively. Each type comes in three variants with specific surface characteristics in the high latitudes. High obliquity causes a strong energy input in the high latitudes during the summer. By increasing the stellar irradiance of a Cryoplanet at high obliquity, the polar ice cap is melted during the summer, resulting in a Near Cryoplanet state. Further increasing the stellar irradiances then generates an Uncapped Cryoplanet with an open polar ocean throughout the year. Finally, upon even stronger irradiance, the ice is completely melted and an Aquaplanet results. Similarly, if starting from an Aquaplanet at strong irradiance and low obliquity, a reduction of irradiance produces first a Near Aquaplanet with seasonal polar ice cover and then a Capped Aquaplanet with a perennial polar sea ice, before transiting to a Cryoplanet. Hence a quite symmetric distribution of stable states results, which suggests that there are many more possibilities for habitability on an exoplanet than hitherto assumed on the basis of commonly used energy balance models.
We have also presented an examination of the atmospheric circulation and the meridional heat transport for the two extreme states of the Aquaplanet and the Cryoplanet. The presence of multiple equilibria is a robust feature of our simulations. It owes its existence to the ice-albedo feedback, which has been first described using zonally averaged energy balance models. As a consequence, hysteresis behavior exists, which implies that the current climate on an exoplanet depends on the evolution of a planet's orbital configuration, in particular the evolution of irradiance, or distance from the star, as well as of obliquity. The zone of habitability in the irradianceobliquity parameter space is more extended if the planet evolves from an aquaplanet than from a cryoplanet.
Our study has some limitations that are primarily associated with the type of climate model that we have employed. A dynamical ocean component, including a dynamic sea ice model, would be an obvious extension. Such a slowly responding component could generate a more complex solution structure by permitting different circulation modes in the ocean, which would modify the transport of heat meridionally. The chemical composition of the atmosphere, in particular the concentrations of greenhouse gases, would also strongly influence the results. By modifying the outgoing longwave radiation due to changed greenhouse gas concentrations, changes in annual mean energy fluxes in addition to the seasonal variations by obliquity would occur that would have an impact on the occurrence and characteristics of the steady states.
In summary, we have demonstrated that the tilt of the planet's rotation axis with respect to the orbit is a crucial parameter in determining habitability under multiple equilibria. While astronomical observations of exoplanetary systems, combined with model simulations, provide reliable estimates of stellar irradiance, the geometry of the planets' orbit, and their masses, some fundamental orbital information to determine habitability on an exoplanet is still elusive. Future efforts of astronomical observation should be directed to determine the obliquity of a planet in order to better assess the conditions of potential habitability on an exoplanet.