Evaluating the effect of four unknown parameters included in a latitudinal energy balance model on the habitability of exoplanets

Abstract Among different models for determining the habitable zone (HZ) around a star, a Latitudinal Energy Balance Model (LEBM) is very beneficial due to its parametricity which keeps a good balance between complexity and simulation time. This flexibility makes the LEBM an excellent tool to assess the impact of some key physical parameters on the temperature and the habitability of a planet. Among different physical parameters, some of them, up until now, cannot be determined by any method such as the planet’s spin obliquity, diurnal period, ocean-land ratio, and pressure level. Here we apply this model to study the effect of these unknown parameters on the habitability of three exoplanets located in the inner, outer, and middle of their optimistic HZ. Among the examined parameters, the impact of pressure is more straightforward. It has a nearly direct relation with temperature and also with the habitability in the case of a cold planet. The effect of other parameters is discussed with details. To quantify the impact of all these unknown parameters we utilize a statistical interface which provides us with the conditional probability on habitability status of each planet.


Introduction
The first confirmation of an exoplanet orbiting a mainsequence star was made in 1995 when a giant planet was found in a four-day orbit around the nearby star 51 Pegasi. It was named 51 Pegasi b (Mayor and Queloz 1995). Following the discovery of 51 Pegasi b the field of exoplanet detection flourished and resulted in launching many projects, especially those based on transits. By now, over 2000 exoplanets have been detected by the Kepler telescope, and the numbers will expand further as a result of other projects like TESS (Transiting Exoplanet Survey Satellite) (Fischer et al. 2015).
Due to the crucial role of liquid water in the biochemistry of life on the Earth, a habitable zone is defined as the spatial extent around a main-sequence star in which the planet can maintain liquid water on its surface over an extended period of time (Kasting et al. 1993). There are different models of studying the habitability of exoplanets, each one with different assumptions. Based on the complexity and dimensions of a model, it contains its own unique set of parameters. Kasting et al. (1993) calculated HZ limits for stars with effective temperatures higher than 3700 K employing a 1-D radiative-convective, cloud-free climate model. Kopparapu et al. (2013) updated HZ limits for stars with effective temperatures between 2600 K and 7200 K to estimate the actual occurrence rate of the Earth-size planets in the HZs around M-dwarfs. Moreover, Kopparapu et al. (2014) have constructed both conservative and optimistic estimates of the HZ from 1D and 3D models. Heng et al. (2015) focused on atmospheric quantities of highly irradiated large exoplanets or Hot Jupiters. They concluded that the most important implication for the habitability of these planets is the equator-to-pole temperature difference response to the variation in orbital and atmospheric parameters. Kopparapu et al. (2016) estimated the inner edge of the HZ for synchronously rotating terrestrial planets around late-K and M-dwarf stars with selfconsistent relationships between stellar metallicity, stellar effective temperature, and the planetary orbital and rotational period. Furthermore, Kopparapu et al. (2017) studied ocean-covered planets at the inner edge of the HZ around late M to mid-K stars with a 3D climate system model. Yang et al. (2013) showed that the atmospheric circulation state affects horizontal heat transport and spatial distribution of clouds which both have an impact on habitability. In addition, planet rotation speed, and ,consequently, the Coriolis force have a significant effect on the planetary albedo and thus the climate of the exoplanet (Yang et al. 2014).
In addition to the Kasting definition of a spatial HZ, there is another way of assessing the concept by evolving a climate model for large numbers of planets on a multidimensional grid of orbital elements, mapping out an "orbital HZ" in this parameter space (Forgan 2013). An Energy Balance Model (EBM) is one of the above-mentioned climate models. After Budyko (1969) and Sellers (1969) a LEBM proved to be useful in climate science. After Williams and Kasting (1997); Franck et al. (2000), and Gaidos and Williams (2004) this model has been utilized in terrestrial exoplanets. Moreover, in Spiegel et al. (2008); Sellers (2011), and Vladilo et al. (2013) it was used for an Earth-like planet with different rotational speeds and various pressures, respectively.
A variety of other works has studied the dependence of climate on surface pressure in 1D or 3D models. Goldblatt et al. (2009) studied the effect of doubling the present atmospheric nitrogen level using a radiative-convective model. Kopparapu et al. (2014) showed that larger planets have wider HZs. Keles et al. (2018) demonstrated that increasing pressure would result in increasing the temperature due to the greenhouse effect. However, in around the pressure of 4 bar, Rayleigh scattering gradually becomes dominated and thus temperature decreases. Charnay et al. (2013) used a 3D global climate model (GCM) to investigate the effect of atmospheric pressure as different solutions to the faint young Sun problem over time. Wolf and Toon (2014) investigated the habitability of the Earth considering the increment of the Sun's brightness. Chemke et al. (2016) showed that higher pressures tend to increase the nearsurface temperature. Rushby et al. (2019) determined the effect of varying fractional and latitudinal distribution of land as a function of host star spectral energy distribution on the overall planetary albedo, climate, and ice-albedo feedback response. Colose et al. (2019) showed that increasing obliquity results in increasing the temperature of the planet because of ice-albedo feedbacks in cold climates and water vapor warm climates. Some of the physical orbit parameters like semi-major axis and eccentricity, which affect habitability, can be achieved by observation. However, many parameters are still unknown. Among them are atmospheric pressure, diurnal period, spin obliquity, and ocean coverage.
For an orbital HZ, we are free to manipulate some parameters which makes this a suitable tool to speculate on cases having different ocean coverages or different spin obliquities. It may create an issue on how changing unknown parameters can change the habitable fraction of a planet and whether this can change the status of a planet from habitable to unhabitable or vice versa. This paper is organized as follows: Section 2 introduces the LEBM, section 3 presents our exoplanet case studies, section 4 discusses the effect of each parameter on the temperature and habitability, and section 5 concludes our results.

Latitudinal Energy Balance Model
The simplest method of climate models, the zerodimensional model, considers the planet as a particle satisfying energy conservation which could be used to solve the planet's effective temperature. This temperature could be different from the surface temperature based on the greenhouse effect of the atmospheric gases that absorb thermal radiation. To have a more realistic view of the climate, the temperature should be latitudinally resolved as a one-dimensional climate model. A LEBM or zonally averaged EBM treats each zonal strip separately and, thus, the rate of transporting of energy between different layers is also considered.
Following previous works, Williams and Kasting (1997) and Spiegel et al. (2008) the main equation of a LEBM is a one-dimensional heat diffusion equation which relates incoming, reflecting and outgoing radiation to the temperature of each layer.
C is the effective heat capacity, T is temperature, D is the diffusion coefficient, which determines the efficacy of zonal exchange heat, x is related to the latitude such that λ = sin −1 x, I is Outgoing Long wave Radiation (OLR), S is instellation, and A represents the albedo of the layer. The (1 − x 2 ) term arises from solving diffusion equation in spherical geometry.

Model Parameters
To Include the pressure, we follow the definitions of Vladilo et al. (2013) for heat capacity, diffusion coefficient, albedo, and OLR, which are as follows.

Heat Capacity
The atmosphere thermal inertia depends on the fraction of the planet's surface which is covered with ocean, land, and ice. Considering different levels of pressure, the heat capacity of the atmosphere becomes: where cp and P are specific heat capacity and pressure of the atmosphere. Index o in this relation corresponds to the values for the Earth. Heat capacity of the atmosphere, C atm,o is set to be 10.1 × 10 6 J m −2 K −1 , and heat capacities of land, C l , ocean, Co, and ice, C i are defined as follows: Co = 210 × 10 6 + C atm,o 1.0 × 10 6 + C atm T < 263 43 × 10 6 + C atm 263 < T < 273 The total heat capacity is equal to The relation between f l , fraction of land, and fo, fraction of ocean, is The fraction of the surface covered by ice is determined by

Diffusion Coeflcient
The diffusion coefficient D which approximate the heat transport by atmospheric circulation is defined such that a planet at 1 A.U. around a star of mass 1M ⊙ with a rotational period of 1 day will reproduce the average temperature profile measured on the Earth. For a planet with a different rotational period, D is proportional to the rotational velocity ωo −2 and for different pressures, D is linearly proportional to the pressure. The diffusion coefficient is thus: ω is the rotational velocity of the planet, and m is the mean molecular weight. Index o corresponds to the values of the Earth (Williams and Kasting 1997).

Albedo
Albedo is defined based on the fraction of land, f l , ocean, fo, ice on them, f i , the fraction of clouds on land, f cl , clouds on water, fcw, and clouds on ice, f ci . ao, a l , and ac are the albedos of ocean, land, and clouds respectively. a io and a il are the albedos of the ice on the ocean and ice on the land. Z * is the zenith distance of the star and α and β are estimated such that the relation α + βZ * yield a good fit to Cess data (Cess 1976). For a more detailed description of each parameter see Vladilo et al. (2013).
Implementing this relation for albedo instead of the temperature dependent relation in Spiegel et al. (2009) prevents the Earth from freezing when the simulation starts at the northern winter solstice or when the solar distance is changed from 1 AU to 1.025 AU. In each configuration, the Earth freezes as a consequence of choosing an approximation for albedo which fails to predict some situations correctly (Appendix A).

Outgoing Long wave Radiation
The OLR is defined as: Where σ is the Boltzmann constant and τ IR is the optical depth of the atmosphere which for atmospheric pressure is defined as: Changing the pressure represents itself in the optical depth which governs the greenhouse effect. It is expected that τ IR has a strong dependence on P since τ IR = κ IR p/g where κ IR is the total absorption coefficient and g the gravitational acceleration. Furthermore, in the range of planetary surface pressures, κ IR is dominated by collisional broadening, which introduces a linear dependence on P of the widths of the absorption lines (Pierrehumbert 2010;Salby 2012).
In absence of line saturation, κ IR should therefore increase linearly with P. Therefore, one might expect that τ IR ∝ P 2 ; however, it is milder in reality and deriving an analytical function is not possible. It is required to perform a series of radiative calculations and tabulate the results as a function of T and P. The detailed description of the procedure can be found in Vladilo et al. (2013).

Incoming Solar Radiation
Assuming a simple main sequence scaling for the luminosity, the total averaged diurnal instellation becomes: µ is the mean diurnal value of µ = cos Z * when the star is above the horizon. The bolometric flux received at a distance of 1 A.U. becomes: The radian half-day length (H) is: δ is the solar declination which is calculated as follows: δo is the obliquity, ϕp is the current orbital longitude of the planet, ϕ peri is the longitude of periastron, and ϕa is the longitude of winter solstice relative to the longitude of periastron.

Simulation Procedure and the Habitability
In order to solve equation 1 we used a staggered grid in which variables are calculated at the center of the grid cells and their derivatives at the cells borders. The spatial resolution is 1.25 ∘ (145 grid points from the north to the south pole), and a global time step was adopted with a stable time step constraint for equation 1: The boundary condition is dT dλ = 0 at the poles. The habitability function of each latitude is defined as: The fraction of the potentially habitable surface at time t is: To define habitability status, we use equation 17 which calculates the habitable fraction of the surface. When a steady-state is reached, we calculate mean H and the standard deviation σ H of the habitability fraction. Finally, the habitability status of the planet is determined according to Table 1 in which we classify the status of the planet as habitable, hot, snowball, or transient (Forgan 2013

Exoplanet case study
Our goal is to study the impact of different unknown parameters of an exoplanet on its temperature and habitability. These parameters include pressure, the fraction of ocean and land, obliquity, and diurnal period. Usually, these parameters are set to the values of the Earth parameters. However, it is just one case among numerous cases which have different values. We want to study as many cases as possible in order to see how planet habitability behaves with a change in those parameters.
To show this effect, we choose three planetsτ Ceti e, Kepler-22b, and Kepler-62f -among the discovered exoplanets. These three planets according to Kopparapu et al. (2013Kopparapu et al. ( , 2014 are respectively on the inner, middle, and outer edge of the optimistic habitable zone. Therefore, each one serves as a representation of other exoplanets located in the same approximate location in their star HZ. The planet τ Ceti e is on the inner edge of the optimistic HZ of a star with a mass of 0.783±0.012 M ⊙ , and a temperature T=5344 ± 50 K. The planet has a semi-major axis of 0.552 AU, eccentricity of 0.05, and orbital period of 162.87 days (Teixeira et al. 2009). The minimum mass of this planet is 3.94 M ⊕ or 0.0124 M j (⊕ and j refer to the Earth and Jupiter respectively) and its radius is not defined (Feng et al. 2017). The planet Kepler-22b, designated as KOI-087.01, (Borucki et al. 2011) orbits within the optimistic HZ of the Sun-like star Kepler-22 with mass 0.970 M ⊙ and effective temperature T e = 5581 ± 44K. The planet is located in an orbit with a semi-major axis of 0.849 AU, zero eccentricity, and orbital period of 286.89 days. This planet has a maximum mass of 35.87 M ⊕ (0.113 M j ) and radius of 2.38 R ⊕ (0.212 R j ) (Borucki et al. 2012;Torres et al. 2015).
The planet Kepler-62f, designated as KOI-701, ) is located on the outer edge of the optimistic HZ of a star with mass 0.69 M ⊙ and effective temperature T e = 4925 ± 70K. The planet has a semi-major axis of 0.718 AU, zero eccentricity and orbital period of 267 days. Kepler-62f has a maximum mass of 34.92 M ⊕ (0.110 M j ) and radius of 1.41 R ⊕ (0.126 R j ) (Borucki et al. 2013;Akeson et al. 2013). Kepler-62f is probably a rocky world because of its Earth-like radius. The radius of Kepler-22b is in the range of a gas dwarf or mini-Neptune, however, being a rocky world without a very dense atmosphere is also possible (Buchhave et al. 2014;Kane et al. 2016).

Results & Discussion
In the simulation, the spin obliquity (S.O.) is varied to (0.0, 22.5, 45.0, 67.5, 90.0) degrees, diurnal period (D.P.) to (0.5, 0.75, 1.0, 1.25, 1.5) days, ocean fraction (O.F.) to (0.1, 0.325, 0.55, 0.775, 1.0), pressure (P) to (0.1, 0.4, 0.7, 1.0, 1.3) atm for τ Ceti e and Kepler-22b, and (0.5, 1.5, 2.5, 3.5, 4.5) atm for Kepler-62f (The effect of changing these variables has been demonstrated in Figures 1 to 4). Since the model is longitudinally averaged, it can only be utilized in situations where the diurnal period is small compared to the orbital period. Considering this limitation, the range of diurnal period is relatively small to not to approach critical values (Forgan 2013). The reason for choosing different values of pressure for Kepler-62f was that in the range of other planets' pressures there was not a notable change in the temperature or the habitability of Kepler-62f. Then, raising the pressure to the aforementioned values showed its effect on temperature and habitability To present the effect of changing the parameters on the temperature and habitability, we monitored the changes at the end of the determined interval by fixing each parameter to its minimum and maximum and letting our parameter of interest to vary across the entire range of values.
The effect of each parameter of interest on minimum, mean, and maximum temperatures and habitability is generally non-linear. Thus, the interpretation is not straight forward. However, an overall description tolerating a few deviations is possible. We explain this in the following sections.

Changing Pressure
In LEBM, the main effect of increasing pressure is increasing the diffusion. Therefore, heat transfers between zonal strips more effectively. The result is that the process of glaciation would be reduced and it would lower the albedo which in turn increases the temperature. This effect is in accordance with the plots of Figure 1. A planet like τ Ceti e is prone to be hot due to its location on the inner edge of the HZ. Therefore, increasing the pressure could decrease the habitability. In Figure 1a, the habitability fraction of this planet starts with a value of 0.67 at a pressure of 0.1 atm, increases to nearly one at a pressure of 0.7 atm, and decreases again to the value of 0.78 at a pressure of 1.3 atm. In the other two cases , what is observed is increasing the habitability, especially in the case of Kepler-62f, where the planet becomes unhabitable at low pressures. In Figures 1a  1b 1c, Kepler-62f has a zero habitability at pressures below 3.5 atm. However, when its pressure increases to 4.5 atm, its habitability has a significant increase and reaches the value of one.
In more accurate models, various mechanism are included. In fact, The effect of surface temperature could not easily be estimated. Increasing surface pressure would result in the increment of albedo via Rayleigh scattering and also increases the greenhouse effect which implementing them is beyond the scope of this paper (Zsom et al. 2013). Furthermore, since increasing the pressure has a high impact on the temperature, other effects such as the moist and runaway greenhouse, which are not implemented in the standard LEBM, should be studied with more cautious. Those two effects are expected to occur on a planet with

Changing Ocean Fraction
According to equation 4 , as the water has a higher amount of heat capacity in comparison to the land, increasing the ocean fraction will increase the total heat capacity. Regarding the required time for a planet to reach a semi steadystate temperature profile, the main effect of heat capacity is to slow down the changes in the temperature profile of the planet. However, if we allow enough time for the planet to evolve, it will finally reach its semi steady-state temperature profile regardless of the value of its heat capacity. This fact has been demonstrated in an Earth-Sun like system in Figure 8. In low spin obliquities, increasing the ocean fraction generally results in slightly increasing the temperature due to the capability of maintaining more heat. This temperature increment makes planets more habitable in two ways. First, the habitability fraction increases for Kepler-22b, located in the middle region of the HZ (See Figure 2 in which this increment ranges from 0.25 in state c to 0.3 in state b). Second, in addition to some increase in the habitability fraction of τ Ceti e, its habitability status changes from transient to habitable in some cases. Look at Figure 5 for these special cases. When the obliquity is very high, increasing the ocean fraction decreases the difference between maximum and minimum temperatures resulting in a flatter temperature profile. This narrowing of temperature differences is most evident for Kepler-22b in the case of Figure 2a where the difference between these two temperatures changes from 176 to 68 degrees. In these cases, because of the particular configuration of the planet in its orbit, the ice belts created on the equator are weaker than ice caps in the case of zero obliquity. The reason is that while the star altitude at the poles of the planet in a zero obliquity case is always zero, it changes from 0-to-90 degrees at the equator of the planet with a 12-hour day cycle in the case of 90-degree obliquity based on the phase of the planet in its orbit. When the planet is semi-habitable (having a habitability fraction be-tween zero and one), this ocean fraction increment leads to an increase in the habitability of the planet as in the case of Kepler-22b, where its habitability increases from 0.11 to 0.22 while it decreases for τ Ceti e as it is prone to get warmer.

Changing Spin Obliquity
When spin obliquity is zero, high latitudes receive a negligible amount of solar energy. Therefore, ice poles are created and spread into other latitudes until the received instellation compensates for the negative effect of albedo and OLR on the stored energy. By increasing the spin obliquity, the poles receive more instellation, which in turn decreases the area of ice caps. In the median obliquities around 45 degrees, since all latitudes receive a moderate amount of instellation, the temperature profile experiences less diversity. When this increment reaches more radical values such as 90 degrees, an ice belt is created on the equator since this is the region with the least instellation. However, the equator of a zero degress spin obliquity planet receives more yearly instellation in comparison to the poles of a 90 degrees spin obliquity planet. Therefore, the ice belt on the equator is narrower and warmer in general. If we compare the maximum, minimum, and mean temperatures of the cases in Figure 3, the general trend of higher temperatures in 90 degrees spin obliquity compared with zero degrees obliquity is evident. Regarding the numerical difference between maximum and minimum temperatures for the planets at 3a, this difference for τ Ceti e changes from 135 K at zero obliquity to 40 K at 45 degrees obliquity, and 82 K at 90 degrees obliquity. These values are 109, 56, and 128 Kelvin for Kepler-22b and 27, 13, and 22 Kelvin for Kepler-62f, respectively.
The habitability fraction follows a more complicated process. When the planet has a relatively moderate or high temperature, adjusting the temperature profile to a more uniform one by increasing the obliquity causes the habitability to increase since there are more regions with a moderate temperature. Reaching more radical values of obliquity makes temperature profiles again more diverse with a higher density of frozen and boiling regions; therefore, habitability decreases. As an example, the habitability of τ Ceti e in the case of Figure 3a starts with a value of 0.52, When the planet has a relatively low temperature, a reverse process occurs. In low obliquities, parts of the planet which receive more instellation are habitable. When spin obliquity increases, it again makes the temperature profile more uniform. This uniformity results in lowering the temperature in above freezing regions and raising the temperature in frozen regions. However, as the instellation is not high enough the overall effect is to turn more regions into frozen rather than turning others to becoming more habitable. Therefore, habitability decreases. In high obliquities there is a similar situation to low obliquities; there are frozen parts with high instellation which begin to melt. These areas are responsible for increasing the habitability. The planet Kepler-22b in the case of Figure 3d is a good example, where its habitability starts from 0.34, decreases to zero, and then increases back to 0.22.
There are some cases for which changing obliquity has no significant effect on habitability. These cases might be completely frozen in any spin condition or the dominant parameter is a parameter other than obliquity (e.g. pressure or diurnal period). Moreover, changing the status of some cases from transient to habitable or vice versa is considerable in Figures 5 and 6. By increasing obliquity, planet τ Ceti e experiences transformation from habitable to transient in low pressures and vice versa in high pressures. The planet Kepler-22b experiences transformation from habitable to transient in moderate pressures.

Changing the Diurnal Period
As mentioned earlier, due to the limitation of the model, the diurnal period should be small enough with respect to the orbital period. According to equation 7, increasing the rotational velocity results in decreasing the diffusion coefficient. This means that heat cannot transfer effectively between different latitudes. This is a result of the combination of the facts that as rotational rate increases, the decreases in eddy length scale leads to less eddy transport poleward. On the contrary, as the rotation slows down, large planetary scale Hadley cells increase the heat transport by the mean meridional circulation (Kaspi and Showman 2015). In fact, we expect that the main impact of increasing the diurnal period should be in changing the latitudinal temperature profile where the difference between the maximum and minimum temperatures decreases. However, the effect of changing this parameter on mean temperature is negligible.  . Simulation for a Sun-Earth system with different diurnal period of the Earth shows that increasing this parameter results in a more homogeneous temperature distribution.
Since increasing the diurnal period converges the maximum and minimum temperatures, if there is an effect on habitability, that effect would be increasing the habitability as a result of making hotter regions colder and colder regions warmer (It is displayed for an Earth-Sun system like simulation in Figure 9). For instance, in Figure 4b, the difference between maximum and minimum temperatures of the planet τ Ceti e starts with a value of 134 K and finally reaches 31 Kelvin. The two other planets follow nearly the same process. The habitability for τ Ceti e increases from 0.53 to 1.0.

Statistical interface to the unknown parameters
Up until now, we investigated the signature of the mentioned unknown parameters on the habitability of an exoplanet. As we see, they have a variety of effects on the temperature and habitability and should be taken into account carefully. The problem with these parameters is that we just know or even suspect a possible range for them but not their exact values. Therefore, to quantify their effect we utilized a statistical approach.
With the values for the parameters introduced in section 2, the simulation's phase space contains 625 points.
To have a better perception of how the temperature and habitability of a planet vary when a selected parameter or a combination of parameters changes, we encapsulated the result of all of the diagrams in one figure called "Map of parameters". There are three of these figures (Figures 5 to  7), one for each of the three planets. In these figures, by moving to the right, the pressure increases in each little square and the obliquity increases by jumping between the large squares. Similarly, by moving to the top, the diurnal period increases in each little square and the ocean fraction increases by jumping between the large squares. The habitability status and mean temperature of each case could be determined using the guidance next to each figure.
From these maps a planet could have different states based on the value of its four unknown parameters. Therefore, setting these values to those of the Earth does not give us an accurate understanding of the habitability of an exoplanet. For instance, Figure 6 shows that 333 out of 625 cases (slightly more than 50 percent) have a habitability more than 0.4. However, it is still an inaccurate inference to say that the planet is likely to be habitable rather than unhabitable.
To assess this issue more accurately, we need to meticulously determine our phase space. For the spin obliquity, the possible range is covered in the simulation. Ocean fraction takes the values from 0 to 1 with the exception of 0.0 since the planet with no water is unhabitable in our assumption of habitability. The pressure could not be 0.0 atm which indicates no atmosphere on the planet and should be excluded from the possible range. The upper limit of 1,100 atm for the pressure is a good choice based on the endurance of extremophiles in an extreme environment (Stan-Lotter 2007). To make a statistical interpretation and comparison, it is essential to specify a single phase space for all three planets. In this regard, the phase space on the pressure coordinate will be confined to the interval (0.0, 1]. Choosing this interval helps us avoid boundaries problem during the extrapolations. Moreover, it is a representation of different atmospheric pressures on the Earth. The value of diurnal period depends on the formation history of the planet and has no apparent bounds. However, using the statistical approach based on a simulation of planet formation with different initial planetary system parameters it is possible to determine the distribution of the diurnal period (Miguel and Brunini 2010). According to those results, the most probable values for diurnal period are larger than 10 hours but the number of planets with diurnal periods between 0.1 and 10 hours is also considerable. Values less than 0.1 hours are really rare, because at those periods the spin angular velocities exceed the critical instability angular velocity value. Therefore, we set the lower limit of 0.1 hours for the diurnal period. The upper limit could not be bounded. Using the probability distribution for a primordial diurnal period introduced in the aforementioned study we can discuss the habitability probability with the help of conditional probability by computing probability of habitability given that the probability of diurnal period lies in a specific range; this is the only way to reconcile the probability and LEBM inherent limitations on diurnal period.
Considering the above discussion, we set the range (0.0 -1.0] atm for pressure, (0.0 -1.0] for ocean fraction, [0.0 -90.0] degrees for spin obliquity, and (0.0042 -1.5] days for diurnal period. Using Miguel and Brunini (2010) the probability of diurnal period being within this mentioned range can be calculated and is equal to 0.1226.
Since the parameters are continuous, to calculate the probabilities we need the area of the region which gives us the desired outcome. Running the simulation for each point in phase space is not possible. To overcome the situation we used the Monte Carlo method. We generated many random coordinates in the domain of possible values. Then, we used a weighted average according to the simulated point to estimate the habitable fraction and habitability status.
To estimate those values in the 4D space of the parameters by using a weighted average, we first determine the nearest neighbor simulated coordinate and the 4D cube in which the random value is located. Then, we calculated the average according to this equation: where vxyzw is the simulated value for each corner of the 4D cube and dxyzw is the normalized-to-step-size Euclidean distance from each nearest neighbors: The coefficient a is set to be the Napierian logarithm of 2.220446 × 10 −16 (the smallest positive number in floatingpoint representation with double precision). With this value for a distance equal to the edge of the cube, the exponential value becomes technically zero in the calculation; if the random value is located on the vertices, the average yields the simulated value. Habitability fraction is a continuous variable which we can estimate using equation 19. Habitability status, on the other hand, is a categorical variable; to average it, we assign +1, 0, and −1 respectively for habitable, transient and unhabitable status. The averaged value will be rounded to the nearest category number.
There are 30000 generated random points which is high enough such that repeating the Monte Carlo process changes the value of outcomes by less than 0.1%. The results for our three surveyed planets are summarized in table 2. To calculate the probability of habitability status (habitable -transient -unhabitable) we need to multiply the habitability probability given that the diurnal period lies in the range of (0.0042 -1.5] by the probability of the diurnal period being in the same range.  (North and Coakley 1979) global tem-perature data shows that the model, despite its simplicity, is precise enough to make it suitable for our survey. The main improvement of our simulation compared with Forgan (2013) and Spiegel et al. (2009) is defining some parameters based on Vladilo et al. (2013). When the albedo is defined by relating it only to the temperature, a global freezing problem occurs which is dependent on the starting point of the planet in its orbit. For instance, by increasing the semimajor axis of the Earth to 1.025 AU, based on the defined starting point, the whole surface of the Earth could turn into ice. This problem is solved when equation 8, which uses the albedo of land, ocean, ice, and cloud portions on each of those surfaces, is applied (see Appendix A). To study the impact of each parameter on the temperature and habitability, three maps of parameters for each exoplanet consisting of 625 different sets of parameters were formed. The effect of a parameter is shown in four different sets of parameters extracted from the maps in which one parameter is free, changing from a minimum to a maximum value. These effects are shown in Figures 1 to 4 . The pressure has a strong, nearly linear, impact on temperature. It is the dominant parameter in the examined range of our simulations. When it is increased, the temperature of the planet is also increased. Habitability, on the other hand, could undergo different scenarios. Depending on the mean temperature of the planet at the starting point, it could increase when more frozen regions start to melt or decrease when more regions go beyond the boiling point of water.
Since the ocean has a much higher thermal inertia, the effect of increasing its ratio is to increase the total heat capacity of the surface. In a short period, it prevents abrupt changes in the temperature; however, over the long term, this feature does not play a role as the planet reaches a stable condition. The other impact of increasing the ocean fraction is reducing the difference between the maximum and minimum temperatures. These effects make the planet more habitable by increasing the habitability fraction or changing the status of the planet from transient to habitable.
Spin obliquity is changed from 0 to 90 degrees. In low obliquities, ice caps are created and therefore the difference between the maximum and minimum temperatures is high. As the obliquity increases, ice caps start to melt and these extreme temperatures approach each other. In very high obliquities, where equators receive the least instellation, ice belts are formed on the equator and thus the difference between the maximum and minimum temperatures starts to increase.
The main impact of the increasing diurnal period is a better heat transfer between latitudes which in turn de-creases the difference between the maximum and minimum temperatures of the planet. The value of the diurnal period does not have a noticeable effect on the mean temperature, and so the habitability does not change by varying this parameter.
Since the LEBM does not depend on the mass and radius of a planet, it is expected that nearly the same results for other exoplanets located in the inner, middle, and outer regions of the HZ would be obtained. However, it should be mentioned that the radius and gravity of a planet could have a significant effect on its climate for an even broader range. These effects are shown by Kopparapu et al. (2014); Kaspi and Showman (2015); Yang et al. (2019).
A better understanding of the changing of unknown parameters simultaneously could be achieved by probing the "Map of Parameters" for each planet in figures 5 to 7. To make any statistical inference from these maps, we need to have the area of the phase space corresponding to our desired outcome. Utilizing a Monte Carlo process, we created a large number of points in the phase space and estimated the value of each point with a weighted averaging method. This approach is inevitable since the random value is continuous and it is not possible to run a simulation for each point.
As mentioned in section 4, there is no upper bound for the diurnal period. Thus, we utilized conditional probability using Miguel and Brunini (2010). We set the probability of the diurnal period to be in the range of (0.0042 -1.5] days, and calculated the conditional probability of habitability of the planet, having the ranges of (0.0 -1.0] atm for pressure, (0.0 -1.0] for ocean fraction, [0.0 -90.0] degrees for spin obliquity, given that the diurnal period is in the aforementioned range. To calculate the probability of habitability status we needed to multiply the conditional probability with the probability of diurnal period being in that specific range.
For each planet, we simulated 30000 points of phase space for the diurnal period to be in the range of (0.0042 -1.5] days, having the ranges of (0.0 -1.0] atm for pressure, (0.0 -1.0] for ocean fraction, and [0.0 -90.0] degrees for spin obliquity. Then, we computed the conditional probability habitability status given that the diurnal period is the in range of (0.0042 -1.5] days. As it is presented in table II, τ Ceti e is mostly habitable, the probability of Kepler-22b to be habitable and transient is almost equal, and Kepler-62 is unhabitable in those aforementioned ranges.

A The effect of changing the simulation's starting point and the planet-star distance in an Earth-Sun like system
The LEBM implemented by Forgan (2013) and Spiegel et al. (2009) works well for the Earth-Sun system. However, it is highly dependent on the starting point of the simulation and the exact 1 A.U. distance of the Sun-Earth system. The problem arises from the fact that their choice for the Albedo is based on only temperature. Equation A1 is implemented by both Forgan (2013) and Spiegel et al. (2009).
]︂ (A1) Figure A1 shows that by implementing A1, if the simulation starts at northern solstice, the planet will freeze totally after about 40 years. Also, Figure A2 shows that if the Earth-Sun distance increases by only 0.025 AU, the planet will freeze after about 5 years. However, if we use the relation implemented by Vladilo et al. (2013), which is described in the Albedo subsection 2.1.3, both problems in turning the planet in an Earth-Sun like system to a frozen one would be solved. The bottom rows of Figures A1 and A2 show this remedy.

B Model validation
We should assess the validity of the model by applying it to an Earth/Sun-like system, comparing its result with real Earth data. In Figure B1a, the green line shows the average taken from the NCEP/NCAR global temperature data and the dashed line is our model simulation for an Earth-like planet orbiting a Sun-like star. It is a good match except in the south polar region (North and Coakley 1979). Figure B1b shows the temperature profile in different latitudes for a one hundred year simulation run. The decreasing temperature from the equator to the poles and the effect of seasonal changes are visible in these figures.