Oct 14, 2013 - fined as the region surrounding a star which, if a terrestrial planet of Earth ..... fluxes. The planet was detected by the PlanetHunte...

0 downloads 2 Views 2MB Size

Printed 26 February 2018

(MN LATEX style file v2.2)

arXiv:1310.3611v1 [astro-ph.EP] 14 Oct 2013

Assessing Circumbinary Habitable Zones using Latitudinal Energy Balance Modelling Duncan Forgan 1⋆ 1

Scottish Universities Physics Alliance (SUPA), Institute for Astronomy, University of Edinburgh, Blackford Hill, Edinburgh, EH9 3HJ, Scotland, UK

Accepted

ABSTRACT

Previous attempts to describe circumbinary habitable zones have been concerned with the spatial extent of the zone, calculated analytically according to the combined radiation field of both stars. By contrast to these “spatial HZs”, we present a numerical analysis of the “orbital HZ”, a habitable zone defined as a function of planet orbital elements. This orbital HZ is better equipped to handle (for example) eccentric planet orbits, and is more directly connected to the data returned by exoplanet observations. Producing an orbital HZ requires a large number of climate simulations to be run to investigate the parameter space - we achieve this using Latitudinal Energy Balance Models (LEBMs), which handle the insolation of the planet by both stars (including mutual eclipses), as well as the planetary atmosphere’s ability to absorb, transfer and lose heat. We present orbital HZs for several known circumbinary planetary systems: Kepler-16, Kepler-34, Kepler-35, Kepler-47 and PH-1. Generally, the orbital HZs at zero eccentricity are consistent with spatial HZs derived by other authors, although we detect some signatures of variability that coincide with resonances between the binary and planet orbital periods. We confirm that Earthlike planets around Kepler-47 with Kepler-47c’s orbital parameters could possess liquid water, despite current uncertainties regarding its eccentricity. Kepler-16b is found to be outside the habitable zone, as well as the other circumbinary planets investigated. Key words: astrobiology, planets and satellites: general, methods:numerical

1 INTRODUCTION The Habitable Zone (HZ) is a useful conceptual tool in investigating the general habitability of planetary systems. It is usually defined as the region surrounding a star which, if a terrestrial planet of Earth mass and similar atmospheric composition were to reside within it, the water upon the planet’s surface would remain liquid (Huang 1959; Hart 1979). The boundaries of the HZ are subsequently defined by the properties of the host star: the outer boundary of the HZ is typically governed by the rate at which CO2 clouds maintain a sufficiently strong greenhouse effect, and the inner edge of the HZ is controlled by the rate of water loss via hydrogen escape and hydrolysis. The HZ is sensitive to the spectrum of the source of insolation - in particular, how strongly the source emits in the infrared (IR). As a result, the inner and outer boundaries of the HZ are a function of the effective temperature of the star. While the majority of the literature utilising the HZ concept has relied on the seminal atmospheric radiative transfer calculations of Kasting, Whitmire & Reynolds (1993), and subsequent parametrisations (e.g. Underwood, Jones & Sleep 2003; ⋆

E-mail: [email protected]

c RAS

Selsis et al. 2007; Kaltenegger & Sasselov 2011), it should be noted that Kopparapu et al. (2013) have since returned to these calculations, updating the atmospheric absorption models and extending the range of stellar effective temperatures calculated. This has the effect of moving the conservative HZ boundaries for the Solar System outwards slightly. For a single star, the HZ boundary conditions are spherically symmetric, and as a result, the single-star HZ is a circular annulus. Therefore, planets of Earth mass and atmospheric pressure/composition on circular orbits within the HZ are expected to possess liquid water, and hence be potentially habitable. If the planet’s orbit is elliptical, but impinges upon the HZ, then it can be habitable depending on the average flux received by the planet over the orbit, or equivalently how long it spends within the habitable zone (Williams & Pollard 2002; Kane & Gelino 2012a,b). Since the first detection of an exoplanet orbiting a main sequence star (Mayor & Queloz 1995), the science of exoplanet detection has quickly revealed a large number of exoplanets, with several residing in their local (single-star) HZs, e.g. Kepler-22b (Borucki et al. 2012), Kepler-62f (Borucki et al. 2013), or the three planets Gliese 667Cc, 667Ce and 667Cf, which occupy the same HZ (Anglada-Escud´e et al. 2013). However, these two Kepler plan-

2

D. Forgan

ets possess radii 1.4 to 2.4 times larger than that of the Earth, and the Gliese 667C planets have masses greater than two Earth masses. Combined with the current ignorance as to their atmospheric composition, it is unclear if these objects are themselves habitable1 . Equally, these objects could possess exomoons which may themselves be habitable (Forgan & Kipping 2013), and the detection of Earth-mass exomoons is now possible with current observations (Kipping et al. 2013). The growing exoplanet population continues to challenge our preconceptions of what can constitute a stable, potentially habitable planetary system. Binary star systems are among the most recent of these exotic systems to be discovered. In S-type binary systems, such as Alpha Centauri, the binary typically has a sufficiently large semimajor axis (of order 10 - 50 AU) that stable planetary orbits exist around either of the two stars. It has been established by numerical simulation (Wiegert & Holman 1997; Quintana et al. 2002, 2007) that S-type binary systems can form planets in habitable regions around one or both stars. In this scenario, provided that the distance between the two stars remains sufficiently large, approximating the system’s habitable zone with two single star HZs placed around the binary components is usually acceptable. If this is not the case, e.g. if the binary eccentricity is large, then more detailed calculations are required (e.g. Forgan 2012; Eggl et al. 2012; Kaltenegger & Haghighipour 2013). In the P-type “circumbinary” systems, the stars orbit sufficiently closely that the planet orbits the system’s centre of mass, and the single-star approximation clearly fails. Kane & Hinkel (2013) produced analytical calculations which approximate the aggregrate stellar flux as a blackbody function, with a peak wavelength equal to that found by adding the flux from both stars. Applying Wien’s Law yields a combined effective temperature, which can then be used in conjunction with the bolometric flux to calculate HZ boundaries using the single star HZ prescriptions (Kasting, Whitmire & Reynolds 1993; Underwood, Jones & Sleep 2003). In a similar vein, Haghighipour & Kaltenegger (2013) also use the single-star HZ prescriptions, weighting the flux received from each star at a given location according to its effective temperature, and searching for the points where the weighted flux equals the flux received from a 1 M⊙ star at the inner and outer boundaries. Both methods produce similar calculations for the combined habitable zones, which can deviate strongly from the circular annuli depending on the binary mass ratio and orbital elements. While habitable zones can be defined spatially as described above, they can also be defined by the set of allowed planetary orbital elements that permit liquid water on their surface. Instead of analytically calculating what we might call “the spatial HZ”, and measuring the time that planets spend within the zone, we can attack the problem numerically, by evolving the climates of many planets on a multidimensional grid of orbital elements, mapping out an “ orbital HZ” in this parameter space. While the orbital HZ may not supply the same level of theoretical insight as a spatial HZ, it does possess two advantages: (i) The spatial HZs in multiple star systems are complex and time-dependent, and hence the time required to calculate a planet’s habitability using the spatial HZ increases quickly as the number of stars in the system increases. Conversely, numerical simulations 1 It should also be noted that liquid water is considered one of the primary necessary conditions for habitability, but it is unlikely to be a sufficient condition. Extrapolating astrobiological data from a single data point (the Earth) is demonstrably difficult (cf Spiegel & Turner 2012)

that produce an orbital HZ typically demonstrate a weaker scaling of compute time with star number. (ii) Simulations such as those used to generate the orbital HZ can incorporate the effect of stellar eclipses easily. For analytic calculations, some parametrisations are available (cf Heller 2012) but this has not yet been done for P-type binary systems. (iii) Exoplanet observations produce orbital parameters as output. As such, astrobiologists adopting an orbital HZ will have a more immediate and profound grasp on the habitability of an exoplanet than they might obtain by constructing a spatial HZ as an intermediate step. Identifying orbital HZs requires running a large number of individual simulations. The climate model used must therefore be fast, robust, and reliable. Latitudinal energy balance models (LEBMs) are well suited to this task (North, Cahalan & Coakley 1981; Williams & Kasting 1997; Williams & Pollard 2002; Spiegel, Menou & Scharf 2008, 2009; Dressing et al. 2010; Spiegel et al. 2010; Forgan 2012; Forgan & Kipping 2013; Vladilo et al. 2013). By splitting the planet into latitudinal strips, making some simplifying assumptions about atmospheric stratification, and the spectral energy distribution of the incoming radiation, LEBMs require very little CPU time to complete a climate simulation that faithfully reproduces climates on Earthlike planets (see e.g. Spiegel, Menou & Scharf 2008 or Vladilo et al. 2013 for examples of tests). In this work, we use LEBMs to assess the orbital HZs in several known circumbinary planetary systems: Kepler-16, Kepler-34, Kepler-35, Kepler-47 and PH1. We investigate the HZ as a function of planet semimajor axis ap and planet eccentricity ep , and compare the LEBM calculations to analytic calculations of circumbinary habitable zones. In section 2 we describe the construction of the LEBM and the initial conditions used; in section 3 we display the resulting orbital HZs produced using the LEBMs for the above circumbinary systems. In section 4, we investigate the dependence of the circumbinary HZs on the orbital parameters of the binary, and suggest routes for future improvement, and in section 5 we summarise the work.

2 METHOD 2.1 Latitudinal Energy Balance Models The LEBM is a one dimensional diffusion equation of surface temperature: ∂T (x, t) ∂ ∂T (x, t) C − D(1 − x2 ) = S(1−A(T ))−I(T ).(1) ∂t ∂x ∂x Rather than using the latitude, λ, directly, the variable x = sin λ is used instead for reasons of computational expediency (the (1 − x2 ) term being a geometric factor arising from the spherical geometry of the problem). This equation is evolved with the boundary condition dT = 0 at the poles (where λ = [−90, 90]◦ ). dx T (x, t) is the surface temperature, C is the effective heat capacity of the atmosphere, D is a diffusion coefficient that determines the efficiency of heat redistribution across latitudes, S is the insolation flux, I is the IR cooling and A is the albedo. In the above equation, C, S, I and A are functions of x (either explicitly, as S is, or implicitly through T ). The diffusion constant D is defined such that a planet at 1 au around a star of 1M⊙ , with rotation period of 1 day will reproduce the average temperature profile measured on Earth (see e.g. c RAS, MNRAS 000, 1–10

Circumbinary Habitable Zones Spiegel, Menou & Scharf 2008). Planets that rapidly rotate experience inhibited latitudinal heat transport, due to Coriolis effects, resulting in a D ∝ ωd−2 scaling, where ωd is the rotational angular velocity of the planet (see Farrell 1990). We therefore use: −2 ωd D = 5.394 × 102 , (2) ωd,⊕ where ωd,⊕ is the rotational angular velocity of the Earth. This expression is certainly too simple to describe the full effects of rotation, as more detailed global circulation modelling indicates (Del Genio 1993, 1996). A more rigorous expression would include the effects of atmospheric pressure and mean molecular weight (e.g. Williams & Kasting 1997, but see also Vladilo et al. 2013’s attempts to introduce a latitudinal dependence to D to mimic the Hadley convective cells on Earth). As in Forgan (2012) and Forgan & Kipping (2013), we solve the diffusion equation using an explicit forward time, centre space finite difference algorithm. A global timestep was adopted, with constraint δt <

(∆x)2 C . 2D(1 − x2 )

(3)

3

where q0 is the bolometric flux received from the star at a distance of 1 AU, and Z is the zenith angle: 4 M q0 = 1.36 × 106 erg s−1 cm−2 (12) M⊙ cos Z = µ = sin λ sin δ + cos λ cos δ cos h.

(13)

Here, we assume the luminosity can be determined from main sequence scaling (M⊙ represents one solar mass). In this form, the model cannot describe binary systems with post-main sequence components, but in principle it can be updated to do so, provided that the spectra of the stars are well described. The solar hour angle is h, and δ is the solar declination, which is calculated from the obliquity δ0 using: sin δ = − sin δ0 cos(φp − φperi − φa ),

(14)

where φ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. As we use diurnally averaged quantities, we must also diurnally average S: S = q0 µ ¯.

(15)

As the system is longitudinally averaged, a key assumption of the model (and its inputs) is that the planet rotates sufficiently quickly relative to its orbital period. We adopt the same input expressions for the atmospheric heat capacity, albedo, insolation and atmospheric cooling as was done in Forgan (2012), which we summarise here. The atmospheric heat capacity depends on what fraction of the planet’s surface is ocean, focean , what fraction is land fland = 1.0 − focean , and what fraction of the ocean is frozen fice :

We do this by integrating µ over the sunlit part of the day, i.e. h = [−H, +H], where H(x) is the radian half-day length at a given latitude. Multiplying by H/π (as H = π if a latitude is illuminated for a full rotation) gives the total diurnal insolation as H q0 S = q0 (H sin λ sin δ + cos λ cos δ sin H) . (16) µ ¯= π π

C = fland Cland + focean ((1 − fice )Cocean + fice Cice ) .

Both stars contribute to the total flux S. We calculate the orbital longitude, solar declination and radian half-day length for both stars, as well as the distance of the planet from both stars. If one star is eclipsed by the other, then we set its contribution to S to zero. We ensure that the simulation can accurately model a transit by adding an extra timestep criterion, ensuring that the transit’s duration will not be less than ten timesteps.

(4)

The heat capacities of land, ocean and ice covered areas are Cland = 5.25 × 109 erg cm−2 K−1

(5)

Cocean = 40.0Cland 9.2Cland Cice = 2Cland

(6) 263 K < T < 273 K T < 263 K.

σSB T 4 , 1 + 0.75τIR (T )

(9)

The albedo function is A(T ) = 0.525 − 0.245 tanh

T − 268 K . 5K

(10)

As the surface temperature drops and water freezes, the albedo increases rapidly and non-linearly. This sets up a positive feedback loop that can make the outer HZ extremely sensitive to small perturbations in dynamical or radiative properties. At any instant, for a single star, the insolation received at a given latitude at an orbital distance r is 2 1AU S = q0 cos Z , (11) r c RAS, MNRAS 000, 1–10

(17)

2.2 Determining the Habitable Zone - Classification of Model Outcomes (8)

where the optical depth of the atmosphere 3 T . τIR (T ) = 0.79 273 K

cos H = − tan λ tan δ.

(7)

The infrared cooling function is I(T ) =

The radian half day length is calculated as

When using LEBMs, it is common to calculate a “habitability function” ξ (see Spiegel, Menou & Scharf 2008): 1 273 K < T (λ, t) < 373 K (18) ξ(λ, t) = 0 otherwise. Strictly, this is a potential habitability function - it simply measures whether a given latitude lies in the temperature range where liquid water may exist. This paper relies heavily on this function, and discussions of habitability refer specifically to the fraction of surface where liquid water may exist. We average this function over latitude to calculate the fraction of potentially habitable surface at time t: Z π/2 ξ(t) = 1/2 ξ(λ, t) cos λ dλ. (19) −π/2

We will use this function to classify the planets we simulate in the following sections. Once each simulation has settled into a quasisteady state, we average ξ over the last ten years of the run, and use ¯ and its standard deviation σξ , to classify as follows: the mean, ξ,

4

D. Forgan 1.0

Table 1. Parameters used in this work to describe each binary system.

0.8

ep

0.6 0.4 0.2

M1 (M⊙ )

M2 (M⊙ )

abin (AU )

ebin

Kepler-16 Kepler-34 Kepler-35 Kepler-47 PH1

0.6897 1.0479 0.8877 1.043 1.384

0.2026 1.0208 0.8094 0.362 0.386

0.224 0.224 0.176 0.0836 0.144

0.15944 0.52087 0.1418 0.0234 0.0

2.3 Initial Conditions

0.0 −0.2

Name

0.6

0.8

ap

1.0 (AU)

1.2

1.4

Figure 1. The habitable zone for an Earth-like planet around a Sun-like star, as calculated from a LEBM using the classification system outlined above. We plot the results for each simulation according to the planet’s semimajor axis (x-axis) and the planet’s eccentricity (y-axis), and the colour of the point indicates its outcome. Red points are hot planets with no habitable surface; blue points are cold planets with no habitable surface; green points represent warm planets with at least ten percent of the surface habitable and low seasonal fluctuations; white circles represent warm planets with high seasonal fluctuations. The dashed lines indicate boundaries between classifications. The green dashed lines indicate the conventional habitable zone.

(i) Habitable Planets - these planets exhibit ξ¯ > 0.1, and σξ < ¯ i.e. the fluctuation in habitable surface is less than 10% of the 0.1ξ, mean. (ii) Hot Planets - these planets have temperatures above 373 K across all seasons, and are therefore completely uninhabitable (ξ¯ < 0.1). (iii) Snowball Planets - these planets are completely frozen and are therefore completely uninhabitable (ξ¯ < 0.1). (iv) Transient Planets - these planets possess a time-averaged ¯ i.e. the fluctuation in habitable surface is ξ¯ > 0.1, but σξ > 0.1ξ, greater than 10% of the mean. Figure 1 shows the single-star habitable zone for the Solar system as it would be classified by the above taxonomy. Note the extension of the habitable zone (as described by the green points) to low semimajor axis at low eccentricity. This is a symptom of only requiring ξ¯ > 0.1 for habitability. As the seasonal variations in climate around low eccentricity planets are relatively low, this allows a planet with a fairly inhospitable surface to maintain small habitable regions at the poles which do not vary greatly in extent. For comparison, the Earth’s parameters exhibits ξ¯ ∼ 0.85. This is much higher than the value required to classify a planet as habitable, and it might be suggested that requiring ξ¯ > 0.1 is not particularly demanding. The habitable zones we delineate here are quite generous, and planets at the edges of the zone will be largely inhospitable, but will still possess regions that remain habitable throughout the season, and as such sufficient to maintain a modest but limited biosphere.

Unless otherwise stated, the planets simulated are assumed to be Earthlike. The diurnal period is set equal to the Earth’s, the obliquity is set to 23.5◦ , and the surface ocean fraction focean is set to 0.7. We fix these parameters for expediency, but we should note that these parameters have their own effects on habitability. Increasing the rotation rate can suppress the latitudinal transport of heat (Farrell 1990). Planets with low surface ocean fractions will experience stronger seasonal temperature variations (Abe et al. 2011; Forgan 2012) which would have obvious consequences for the classification system used in this paper. Planets with larger obliquity appear to resist the “snowball” transition to a completely frozen state, even when rapid rotation would otherwise encourage it (Spiegel, Menou & Scharf 2009). The planets orbit in the binary plane, around the centre of mass of the binary system (with the exception of comparison simulations run without the secondary). The simulations begin at the northern winter solstice, which is assumed to occur at an orbital longitude of 0◦ . In the case of eccentric orbits, this is also the longitude of periastron2 . The planets’ initial temperature was set to 288 K at all latitudes. Each simulation is run for a sufficient length of time that the planet’s temperature profile reaches a periodic, steady state, such that the habitability classifications described earlier can be made. Table 1 lists the input parameters for all binary systems studied in this paper.

3 RESULTS 3.1 Kepler-16 Kepler-16 was the first circumbinary planetary system to be discovered during the Kepler mission (Doyle et al. 2011). Kepler-16b, with mass 0.3 MJup , orbits the binary with a period of 229 days, while the binary orbital period is 41 days. The left panel of Figure 2 shows the habitable zone in ap − ep space for the Kepler-16 binary system. For comparison, the right hand panel of the same figure shows the equivalent habitable zone in the absence of the secondary star. Given the large mass difference between the primary and secondary, it is not surprising that HZs produced with and without the secondary are so similar. What is more interesting is the switch in habitability classification for some parameters from habitable to transient as the secondary star is added. This might be expected for planets near the inner HZ edge on the right panel of Figure 2 - the extra insolation from the second star, coupled with eclipses of the primary by the secondary, can produce temperature variations of sufficient strength 2

Simulations were carried out where the longitude of periastron was varied. As the habitability calculations average over many orbits, the effect of changing the initial phase is minimal c RAS, MNRAS 000, 1–10

1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

ep

ep

Circumbinary Habitable Zones

0.2

0.2

0.0

0.0

−0.2

Kepler 16b 0.2

0.4

ap

0.6 (AU)

0.8

1.0

−0.2

5

Kepler 16b 0.2

0.4

ap

0.6 (AU)

0.8

1.0

Figure 2. Left: The habitable zone for an Earth-like planet in the Kepler-16 binary system; Right: the habitable zone for an Earth-like planet in the Kepler-16 system with the secondary removed. The colour of points represents the classification of each simulation run. As before, red points are hot planets with no habitable surface, blue points are cold planets with no habitable surface, green points are warm planets with at least ten percent of the surface habitable and low seasonal fluctuations, and white circles are warm moons with high seasonal fluctuations. Again, dashed lines denote boundaries between different classifications.

to periodically push large fractions of the planet’s surface above 373K. Indeed, this change is strongest at around 0.35 AU, which corresponds to a 2:1 resonance between the planet and binary orbital periods. However, we also see this reclassification at the outer HZ edge (approximately ap = 0.6AU ) at eccentricities above ep = 0.6, which is more surprising. The secondary insolation at this distance from the binary is less than a few percent of the primary insolation. So how is the surface temperature so strongly affected? At 0.6 AU, planets with eccentricities greater than 0.6 will have periastra located inside the binary’s orbit. These very close approaches to the secondary will produce climate variations that ensure the planet is classified as transient. The orbital HZ constructed here would suggest that Kepler16b, which has a close to circular orbit at 0.7 AU, is too cold to be in the habitable zone, as confirmed by Kane & Hinkel (2013) and Haghighipour & Kaltenegger (2013). It has been suggested that Kepler-16b could host a habitable Earth-mass captured planet in a satellite or Trojan orbit (Quarles, Musielak & Cuntz 2012), but this possibility is outside the scope of this work.

3.2 Kepler-34 Kepler-34 is a circumbinary planetary system possessing two G type stars, first reported in Welsh et al. (2012), with stellar masses very close to equal, in a 28 day orbit. The semi-major axis of the binary orbit is similar to that of Kepler-16, but the eccentricity is quite large. Figure 3 shows the habitable zones in the case where the secondary of Kepler-34 is either present (left panel) or absent (right panel). The effect of adding the secondary to the system is significant, pushing the outer HZ boundary from around 1.1 AU at zero eccentricity to around 1.5 AU. This shift is so extreme that there are very few simulations that are classified as habitable (green) in both cases (including the parameters corresponding to the planet Kepler-34b). The inner and outer edges of the habitable zone meet at a peak at e = 0.9 in the single star case - in the two star case, the height of this peak is reduced from e = 0.9 to 0.8. c RAS, MNRAS 000, 1–10

Kepler-34b has a semi-major axis of 1.0896 AU, with an eccentricity of 0.182 (Welsh et al. 2012). The simulation corresponding most closely to these parameters is classified as transient, although this is somewhat moot given that the planet mass is 0.22 MJup . This being the case, it could still be a promising host for a habitable exomoon, as the cooling effect of eclipses of the moon by Kepler-34b itself may help to make the surface more clement (Heller 2012; Forgan & Kipping 2013). Ironically, if the secondary is removed from the simulation (right panel), Kepler-34b is very much inside the habitable zone. 3.3 Kepler-35 Reported alongside Kepler-34 in Welsh et al. (2012), Kepler-35 also consists of two roughly equal mass G stars, but in a low eccentricity orbit with a period of 20 days. Figure 4 shows the habitable zones derived for this system (with and without the secondary star). Again, as the binary masses are close to equal, the HZ boundaries shift significantly, and the highest eccentricity orbits are no longer continuously habitable. Kepler-35b orbits with a period of 131 days at a semi-major axis of 0.6 AU, in a low eccentricity orbit (ep = 0.042). Our calculations indicate that Kepler-35b is too hot to be within the habitable zone, and it would seem unlikely that any moon it might possess would be habitable either. Again, on removal of the secondary, the exoplanet would be in the habitable zone (in this case near the inner edge). 3.4 Kepler-47 This binary system has the distinction of being the first P-type with multiple planets detected in orbit (Orosz et al. 2012). Consisting of a G and M star in a tight low eccentricity orbit, the system has two planets orbiting in the binary plane at 49 days (Kepler-47b) and 303 days (Kepler-47c). The outer planet is thought to be in the habitable zone, although with the eccentricity of the planet established only as an upper limit (ep < 0.411), it is unclear how long it will

D. Forgan 1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

ep

ep

6

0.2

0.2

0.0

0.0

Kepler 34b

−0.2 0.8

1.0

1.2

1.4

ap

(AU)

1.6

1.8

2.0

Kepler 34b

−0.2 0.8

2.2

1.0

1.2

1.4

ap

(AU)

1.6

1.8

2.0

1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

ep

ep

Figure 3. Left: The habitable zone for an Earth-like planet in the Kepler-34 binary system; Right: the habitable zone for an Earth-like planet in the Kepler-34 system with the secondary removed. The points are coloured according to the same classification system as previous figures.

0.2

0.2

0.0

0.0

−0.2

Kepler 35b 0.6

0.8

1.0 ap (AU)

1.2

1.4

−0.2

Kepler 35b 0.6

0.8

1.0 ap (AU)

1.2

1.4

Figure 4. Left: The habitable zone for an Earth-like planet in the Kepler-35 binary system; Right: the habitable zone for an Earth-like planet in the Kepler-35 system with the secondary removed. The points are coloured according to the same classification system as previous figures.

spend in the spatial HZ, as noted by both Kane & Hinkel (2013) and Haghighipour & Kaltenegger (2013). Figure 5 shows the orbital HZ constructed for the Kepler-47 binary system and for the Kepler-47 primary alone. The presence of the M star has little effect on the HZ - the increased flux allows low eccentricity planets to be more habitable at semi-major axes between 0.7 and 0.9 AU, but otherwise there is little else to report.

Kepler-47b has a near circular orbit at 0.2956 AU, and is clearly not in the habitable zone. Despite its currently uncertain eccentricity, Kepler-47c does indeed appear to be warm and habitable, with low climate variations. If Kepler-47c possessed an eccentricity larger than around 0.5, then it would fall into the region of parameter space occupied by transient classifications. Again, we should really only consider moons of Kepler-47c for Earthlike habitability, as the planet is Neptune-sized (Orosz et al. 2012).

3.5 PH-1 PH-1 (also designated Kepler-64) is a quadruple star system, with the planet PH1b orbiting a F and M binary system. The other two stars in the PH1 system form a separate binary which orbits at a distance of 1000 AU, which is sufficiently distant to neglect their fluxes. The planet was detected by the PlanetHunters citizen science program (Schwamb et al. 2013), and orbits with a period of 138 days. This system possesses the most extreme stellar mass ratio, and this is reflected in Figure 6, which shows the orbital HZs constructed both with and without the presence of the M star. The two figures are close to identical, with the exception of low eccentricity, low semimajor axis planets becoming slightly more habitable when the secondary is added, thanks to the cooling effect of eclipses. In the single star case (right panel), the inner and outer HZ edges meet at ep = 0.6 - with the addition of the second star, the inner and outer edges no longer meet, as the outer edge is pushed to 1.99 AU. The binary orbital period is 20 days (as the eccentricity is not c RAS, MNRAS 000, 1–10

2.2

1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

ep

ep

Circumbinary Habitable Zones

0.2

0.2

0.0

0.0

−0.2

Kepler47c 0.6

0.8

ap

1.0 (AU)

1.2

1.4

−0.2

7

Kepler47c 0.6

0.8

ap

1.0 (AU)

1.2

1.4

Figure 5. Left: The habitable zone for an Earth-like planet in the Kepler-47 binary system; Right: the habitable zone for an Earth-like planet in the Kepler-47 system with the secondary removed. The points are coloured according to the same classification system as previous figures. Kepler-47b is not marked as it is at lower semimajor axis than shown here. Kepler-47c is marked as a downward triangle to represent that the eccentricity measured is an upper limit.

constrained by observations, we assume the binary orbit is circular). Hence, as habitable planets will orbit with periods of 500 days or more, the effect of such frequent eclipses on the planetary climate is softened by the atmospheric thermal inertia of the planet. As a 0.5 MJup planet orbiting well within the inner HZ boundary, it is clear that PH1b is not habitable, and is unlikely to host habitable moons.

not considered here), then the dynamical landscape is quite rich (see e.g. Doolin & Blundell 2011). Kepler-16 does not have a stable orbital HZ according to this analysis, whereas the other systems have HZs well clear of the instability limit. We should note that the LEBM may underestimate the extent of the outer HZ boundary due to some missing physics (see section 4.3), but it remains clear that low mass circumbinary systems are not promising candidates for stable habitable zones.

4 DISCUSSION

4.2 Dependence on Binary Orbital Elements

4.1 Stability of Orbits in the Habitable Zone

The P-type systems investigated here give a useful indication of how the habitable zone changes as the binary mass ratio µ, binary semimajor axis abin and binary eccentricity ebin are altered. However, it is also instructive to select a single system and vary its parameters. We use Kepler-35 for this analysis, as the previous shows that systems with µ ∼ 0.5 tend to produce orbital HZs more sensitive to the binary orbital elements (in line with previous spatial HZ analyses), and we wish to consider HZs that are still orbitally stable as ebin and abin are changed. Figure 7 shows the Kepler-35 system orbital HZ as ebin is increased from 0.1418 to 0.3 (left panel) and 0.5 (right panel). The minimum stable semimajor axes are amin = 0.567 AU and 0.642 AU respectively. The habitability of planets at low eccentricity and low semi-major axis appears to increase as a result of increasing ebin , while remaining dynamically stable. The inner and outer edges of the habitable become extended in eccentricity, all the way to ep = 0.9. Objects with these orbital elements are likely to still possess strong temperature fluctuations on seasonal and dynamical timescales, but they appear to be just low enough to be less than 10% of the mean temperature, and are therefore conventionally habitable. If we now return ebin to the standard value for Kepler-35 and increase abin instead, we can see that the well-defined boundaries between habitable and non-habitable regions begin to blur. An increase of abin from 0.176 AU to 0.25 AU (left panel of Figure 8) prevents planets with high eccentricity being habitable, and deepens a region of variable habitability at ap = 0.82 AU, extending it downwards to ep = 0.05. This corresponds to a planetary orbital

In this work, we do not model the gravitational force of the binary upon the planet. Instead, we simply assume fixed Keplerian orbits for the planets around the barycentre of the binary system. By doing so, we ignore the complicating factors of forming habitable planets with the orbital elements within our parameter study, either in situ or through subsequent migration (cf Meschiari 2012 and Dunhill & Alexander 2013’s studies of Kepler-16b). More importantly, we have not yet considered if these orbits are expected to be stable. Holman & Wiegert (1999) simulate the motion of test particles in both P-type and S-type planetary systems, and produce fitting formulae for a critical semi-major axis for orbital stability. In the case of P-type systems, this critical semi-major axis is a minimum, and is sensitive to the binary’s orbital elements and the mass ratio of the binary µ = M2 /(M1 + M2 ): ap > amin = abin 1.6 + 5.1ebin + 4.12µ + 2.22e2bin −4.27µebin − 5.09µ2 + 4.61µ2 e2bin .

(20)

From these equations, it becomes clear that the frequency of planetary systems with potentially habitable planets will depend in the first instance on the Galactic binary fraction and their orbital statistics (see e.g. Parker & Quanz 2013). Note that this stability prescription ignores the potential for instability islands at ap > amin due to mean motion resonances, and if orbits out of the binary plane are permitted (which we have c RAS, MNRAS 000, 1–10

D. Forgan 1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

ep

ep

8

0.2

0.2

0.0

0.0

−0.2

1.2

1.4

ap

1.6 (AU)

1.8

−0.2

2.0

1.2

1.4

ap

1.6 (AU)

1.8

2.0

1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

ep

ep

Figure 6. Left: The habitable zone for an Earth-like planet in the PH1 binary system; Right: the habitable zone for an Earth-like planet in the PH1 system with the secondary removed. The points are coloured according to the same classification system as previous figures. PH1b is not plotted as its semimajor axis is less than 1.2 AU.

0.2

0.2

0.0

0.0

−0.2

Kepler35b 0.6

0.8

1.0 ap (AU)

1.2

1.4

−0.2

Kepler35b 0.6

0.8

1.0 ap (AU)

1.2

1.4

Figure 7. Left: The habitable zone for an Earth-like planet in the Kepler-35 binary system if the binary eccentricity is increased to ebin = 0.3. Right: the same, with the binary eccentricity increased to ebin = 0.5.

period of 208 days, which is close to a 6:1 resonance with the binary (which now has an orbital period of 35 days). The minimum stable orbit for this configuration is 0.705 AU, implying a small fraction of the inner HZ at low ap , low ep cannot be considered habitable.

Increasing abin to 0.3 AU (right panel of Figure 8) serves only to exacerbate these issues. The variable habitability at 0.82 AU now extends to circular planetary orbits, despite no longer corresponding to any orbital resonance. The inner and outer boundaries of the HZ at high ep are no longer smooth, with the ability to determine whether a planet is continuously or variably habitable (green or clear) depending sensitively on the time period used to carry out the averaging. The minimum stable orbit moves to 0.845 AU, rendering most of the inner HZ dynamically unstable.

4.3 Limitations of the Model Using a LEBM by definition requires some concessions to simplicity, especially if the goal is to run a large number of simulations. However, we acknowledge that there are some potential improvements that should be considered in future work. As previously mentioned, the orbits of the planets are fixed and Keplerian. A more accurate representation would involve specifying a Keplerian orbit as initial conditions, and allowing the orbit to evolve under the gravitational influence of both stars in the system, either via full N Body calculations (e.g. Meschiari 2012) or using analytic expressions which assume the planet mass to be negligible (e.g. Leung & Lee 2013). The non-Keplerian orbits that result from these calculations will add important variations in climate, which may make planets previously classified as “habitable” into “transient” planets (especially the low semi-major axis, low eccentricity planets that are commonly classified as habitable), or set up long term Milankovitch cycles (Spiegel et al. 2010). c RAS, MNRAS 000, 1–10

1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

ep

ep

Circumbinary Habitable Zones

0.2

0.2

0.0

0.0

−0.2

Kepler35b 0.6

0.8

ap

1.0 (AU)

1.2

1.4

−0.2

9

Kepler35b 0.6

0.8

ap

1.0 (AU)

1.2

1.4

Figure 8. Left: The habitable zone for an Earth-like planet in the Kepler-35 binary system if the binary semimajor axis is increased to abin = 0.25 AU. Right: the same, with the binary semimajor axis increased to abin = 0.3 AU (right).

Also, the precession of periastron in circumbinary systems will strongly affect the range of seasonal variations eccentric planets will experience, as the longitude of solstices shifts further from the longitude of periastron (Doolin & Blundell 2011; Armstrong et al. 2013). Kepler-34b is expected to undergo a complete cycle of periapse precession in around 20,000 years, whereas Kepler-35b is expected to do so in less than 10,000 years (Welsh et al. 2012). These timescales are similar to the aforementioned Milankovitch cycles measured on Earth. When comparing this paper to analytical calculations of the spatial HZ, we find that our calculation of the outer HZ boundary (at zero eccentricity) is typically lower than that of the other authors. This is most likely due to the rapid snowball albedo effect present due to the freezing of ice, which does not have a counteropposing mechanism to suppress it (excluding adding more radiation to melt the ice). In reality, more accurate modelling of the carbon-silicate cycle (Williams & Kasting 1997) would allow cooler planets to modify their atmospheric CO2 levels. As such, we do not fully model the “maximum greenhouse” conditions that are a standard of spatial HZ calculations, and this is an important feature that must be added to future models.

the stars are of approximately equal mass. If the primary is much more massive than the secondary, then the single star HZ and circumbinary HZs are very similar, although we note that in the case of Kepler-16, which contains two low mass stars with a relatively large mass difference, eclipses can become important. Of the circumbinary planets orbiting the binaries we investigated, Kepler-47c was the only planet found to reside within the habitable zone. We are able to make this determination despite the uncertainty of the planet’s eccentricity, an advantage of orbital HZ modelling. Kepler-47c is therefore an interesting target for future exomoon detections, as while Kepler-47c is not Earthlike, it may possess terrestrial moons. Kepler-34b would be marginally habitable if it were of Earth mass - if it possesses an Earthlike moon, it may be able to sustain a biosphere, but it would need to be robust against strong oscillations in the moon’s climate. While we have not explicitly simulated the orbital stability of these planets, previous analytical calculations of the minimum semimajor axis for a stable circumbinary orbit indicate that with the exception of Kepler-16b, the HZs produced in this work should be amenable to terrestrial planets on stable orbits. However, the dynamical complexity of circumbinary systems warrants further investigation with more appropriate gravitational physics included.

5 CONCLUSIONS We have used one dimensional latitudinal energy balance modelling (LEBM) to investigate the habitable zones (HZs) of planets orbiting P-type star systems. By running many models for each star system, an “orbital HZ” can be produced, which maps out the HZ in terms of the planet’s orbital elements. This numerical analysis is complementary to the common practice of mapping out the HZ in terms of its spatial extent using analytical calculations. With the use of LEBMs, the orbital HZ allows the effects of stellar eclipses and planet eccentricity to be more simply incorporated. We apply this technique to the circumbinary planetary systems Kepler-16, Kepler-34, Kepler-35, Kepler-47 and PH1. In general, our orbital HZs are consistent with the spatial HZs derived by other authors (e.g. Kane & Hinkel 2013 and Haghighipour & Kaltenegger 2013). As has been found previously, the habitable zone strongly deviates from the single star HZ when c RAS, MNRAS 000, 1–10

ACKNOWLEDGMENTS DF gratefully acknowledges support from STFC grant ST/J001422/1. The author would like to thank the referee for their invaluable comments which greatly improved this paper.

REFERENCES Abe Y., Abe-Ouchi A., Sleep N. H., Zahnle K. J., 2011, Astrobiology, 11, 443 Anglada-Escud´e G. et al., 2013, A&A, 556, A126 Armstrong D. et al., 2013, MNRAS, in press Borucki W. J. et al., 2013, Science (New York, N.Y.), 340, 587 Borucki W. J. et al., 2012, ApJ, 745, 120 Del Genio A., 1993, Icarus, 101, 1

10

D. Forgan

Del Genio A., 1996, Icarus, 120, 332 Doolin S., Blundell K. M., 2011, MNRAS, 418, 2656 Doyle L. R. et al., 2011, Science (New York, N.Y.), 333, 1602 Dressing C. D., Spiegel D. S., Scharf C. A., Menou K., Raymond S. N., 2010, ApJ, 721, 1295 Dunhill A., Alexander R., 2013, MNRAS, in press Eggl S., Pilat-Lohinger E., Georgakarakos N., Gyergyovits M., Funk B., 2012, ApJ, 752, 74 Farrell B. F., 1990, Journal of Atmospheric Sciences, 47, 2986 Forgan D., 2012, MNRAS, 422, 1241 Forgan D., Kipping D., 2013, MNRAS, 432, 2994 Haghighipour N., Kaltenegger L., 2013, ApJ, in press Hart M. H., 1979, Icarus, 37, 351 Heller R., 2012, A&A, 545, L8 Holman M. J., Wiegert P. A., 1999, The Astronomical Journal, 117, 621 Huang S.-S., 1959, PASP, 71, 421 Kaltenegger L., Haghighipour N., 2013, ApJ, in press Kaltenegger L., Sasselov D., 2011, ApJ, 736, L25 Kane S. R., Gelino D. M., 2012a, Astrobiology, 12, 940 Kane S. R., Gelino D. M., 2012b, PASP, 124, 323 Kane S. R., Hinkel N. R., 2013, ApJ, 762, 7 Kasting J., Whitmire D., Reynolds R., 1993, Icarus, 101, 108 Kipping D. M., Forgan D., Hartman J., Nesvorny D., Bakos G. A., Schmitt A. R., Buchhave L. A., 2013, ApJ, in press Kopparapu R. K. et al., 2013, ApJ, 765, 131 Leung G. C. K., Lee M. H., 2013, ApJ, 763, 107 Mayor M., Queloz D., 1995, Nature, 378, 355 Meschiari S., 2012, ApJ, 752, 71 North G., Cahalan R., Coakley J., 1981, Rev. Geophys. Space Phys., 19, 91 Orosz J. A. et al., 2012, Science (New York, N.Y.), 337, 1511 Parker R. J., Quanz S. P., 2013, MNRAS, in press Quarles B., Musielak Z. E., Cuntz M., 2012, ApJ, 750, 14 Quintana E. V., Adams F. C., Lissauer J. J., Chambers J. E., 2007, ApJ, 660, 807 Quintana E. V., Lissauer J. J., Chambers J. E., Duncan M. J., 2002, ApJ, 576, 982 Schwamb M. E. et al., 2013, ApJ, 768, 127 Selsis F., Kasting J. F., Levrard B., Paillet J., Ribas I., Delfosse X., 2007, A&A, 476, 1373 Spiegel D. S., Menou K., Scharf C. A., 2008, ApJ, 681, 1609 Spiegel D. S., Menou K., Scharf C. A., 2009, ApJ, 691, 596 Spiegel D. S., Raymond S. N., Dressing C. D., Scharf C. A., Mitchell J. L., 2010, ApJ, 721, 1308 Spiegel D. S., Turner E. L., 2012, Proceedings of the National Academy of Sciences of the United States of America, 109, 395 Underwood D., Jones B., Sleep P., 2003, International Journal of Astrobiology, 2, 289 Vladilo G., Murante G., Silva L., Provenzale A., Ferri G., Ragazzini G., 2013, ApJ, 767, 65 Welsh W. F. et al., 2012, Nature, 481, 475 Wiegert P. A., Holman M. J., 1997, The Astronomical Journal, 113, 1445 Williams D., Kasting J. F., 1997, Icarus, 129, 254 Williams D. M., Pollard D., 2002, International Journal of Astrobiology, 1, 61

c RAS, MNRAS 000, 1–10