Circumbinary habitable zones in the presence of a giant planet

Determining habitable zones in binary star systems can be a challenging task due to the combination of perturbed planetary orbits and varying stellar irradiation conditions. The concept of"dynamically informed habitable zones"allows us, nevertheless, to make predictions on where to look for habitable worlds in such complex environments. Dynamically informed habitable zones have been used in the past to investigate the habitability of circumstellar planets in binary systems and Earth-like analogs in systems with giant planets. Here, we extend the concept to potentially habitable worlds on circumbinary orbits. We show that habitable zone borders can be found analytically even when another giant planet is present in the system. By applying this methodology to Kepler-16, Kepler-34, Kepler-35, Kepler-38, Kepler-64, Kepler-413, Kepler-453, Kepler-1647 and Kepler-1661 we demonstrate that the presence of the known giant planets in the majority of those systems does not preclude the existence of potentially habitable worlds. Among the investigated systems Kepler-35, Kepler-38 and Kepler-64 currently seem to offer the most benign environment. In contrast, Kepler-16 and Kepler-1647 are unlikely to host habitable worlds.


INTRODUCTION
Over the past three decades exoplanet researchers have discovered more than four thousands planets outside our Solar System 1 . Improvements in detection techniques have now reached a point where finding planets of similar size to our Earth has become a reality. Four of the planets in the Trappist-1 system [17] or TOI-700d [16,36] are prime examples. Several planets in those systems orbit their host star inside the habitable zone. The habitable zone is the region where a terrestrial planet on a circular orbit about its host star can support liquid water on its surface [24]. A number of planets have been found to reside in binary stars systems, some of which even orbit both stars [e.g. 4,43,27]. We refer to the latter as circumbinary planets. Unresolved questions regarding the formation and dynamical evolution of such systems [30] have motivated a number of studies in recent years, particularly on whether or not such systems could host potentially habitable worlds [23,18,21,2,8,3,39,1,46]. Some of the challenges in assessing habitability in binary star systems arises from the fact that one has to account for two sources of radiation, possibly of different spectral type. Moreover, the distance of the planet to each star keeps changing over time in a non-trivial manner due to gravitational interactions between the planet and the two stars.
The introduction of "dynamically informed habitable zones" allowed [7] to study the prospects for habitability of planets orbiting a single star in a binary star system. Dynamically informed habitable zones for systems with a potentially habitable world on a circumbinary orbit were developed in [5] and [6]. In this work, we improve on previous analytic estimates for circumbinary dynamically informed habitable zones and extend the concept to systems that are known to host an additional giant planet. Here we refer to giant planets as bodies with masses ranging from Neptune mass up to a few Jupiter masses.
The structure of this paper is as follows: in the next section we explain the general principles behind dynamically informed habitable zones and construct the required tools to extend the concept to include a giant planet in the system. In section 3 we investigate the potential of binary star systems with known circumbinary planets observed during the Kepler mission [e.g. 4,43,27] to host additional potentially habitable worlds. Finally, in section 4 we summarize and discuss the work presented here.

DYNAMICALLY INFORMED HABITABLE ZONES
The nearly circular orbit of the Earth around the Sun ensures that the planet receives an almost constant amount of radiation on a permanent basis. That assumption falters for a circumbinary planet, however. The second star provides an additional source of radiation, and more importantly, it is also a source of gravitational perturbations for the planetary orbit. Even if a planet is on an initially circular orbit around the binary, the orbit will become elliptic over time [see e.g. 12,14]. As a consequence, the planet experiences time dependent irradiation. Depending on how effectively the climate on the planet can buffer changes in incoming radiation its response to radiative forcing can differ widely [see e.g. 34,42,19,45]. In order to capture the various responses defined by a planet's climate inertia, we make use of so-called 'dynamically informed habitable zones' (DIHZs) [15]. DIHZs not only take the orbital evolution of the planet around the binary into account, but they can even trace habitable zone limits for different climate inertia values a planet may have. Climate inertia is defined as the time it takes climate parameters, such as the mean surface temperature, to react to radiative forcing. The faster the mean surface temperature changes, the lower the climate inertia of a planet is.
In order to account for the effect of climate inertia, we follow the general methodology as outlined in [7] and [15] and introduce three different DIHZs: the permanently habitable zone (PHZ), the averaged habitable zone (AHZ) and the extended habitable zone (EHZ). The PHZ is the most conservative region. For a planet to reside in the PHZ means that it stays continuously within habitable insolation limits in spite of orbit-induced variability. In other words, the PHZ is the region where a planet with essentially zero climate inertia could remain habitable on stellar evolution timescales. On the other end of the spectrum is the averaged habitable zone, where a planet is assumed to buffer all variations in irradiation and remains habitable as long as the insolation average stays within habitable limits. This scenario corresponds to a planet with a very high climate inertia.
The extended habitable zone lies between the above extremes by assuming that the planetary climate has limited buffering capabilities. The EHZ is defined as the region where the planet stays on average plus minus one standard deviation within habitable insolation limits.
In order to calculate DIHZ borders for circumbinary systems, we need to understand a) whether or not a configuration is dynamical stable, b) how the orbital evolution of the system affects the amount of radiation the planet receives and c) how the combined quantity and spectral distribution of the star light influences the climate of a potentially habitable world.

Dynamical stability
Dynamical stability is a necessary condition for habitability of a circumbinary planet. If a planet is ejected from a system, water that may be present on its surface will ultimately freeze. There are a number of ways to predict whether or not a potentially habitable world is on a stable orbit [e.g. see 11]. In this work, we make use of the empirical condition developed in [20] which provides critical semi-major axis values below which planetary orbits around binary stars become unstable. The critical semi-major axis in a circumbinary system depends on the eccentricity, semi-major axis and mass ratio of the binary as can be seen from the following equation: where a b is the semi-major axis of the binary, e b is the eccentricity of the binary and µ = m 2 /(m 1 + m 2 ); m 1 and m 2 are the masses of the two stars (m 1 > m 2 without loss of generality). A value of planetary semi-major axis below a c indicates that the planet will escape from the system or collide with one of the stars.
Once we confirm that a potentially habitable planet is on a stable orbit, we can proceed to investigate how much radiation it receives from the two stars. The latter depends on the orbit of the planet which evolves over time. By modeling the evolution of the orbit of the planet we can estimate the actual amount and spectral composition of the radiation the planet receives. To this end we make use of an analytic orbit propagation technique for circumbinary planets developed in [14].

The Classical Habitable Zone
For systems consisting of a star and a terrestrial planet on a fixed circular orbit, the limits of the classical habitable zone (CHZ) read where r is the distance of the planet to its host star in astronomical units, L is the host star's luminosity in solar luminosities, and S I and S O are effective insolation values or "spectral weights" [24]. The latter correspond to the number of solar constants that trigger a runaway greenhouse process (subscript I) evaporating surface oceans, or a snowball state (subscript O) freezing oceans on a global scale. Spectral weights are functions of the effective temperature of the host star and, therefore, they take the specific wavelength distribution of a star's light into account. They can be calculated by the expressions below which can be found in [26]: S O = 0.356 + 6.171 · 10 −5 T c + 1.698 · 10 −9 T 2 c − 3.198 · 10 −12 T 3 c − 5.575 · 10 −16 T 4 c .
T c = T ef f /1 K − 5780, with T ef f being the effective temperature of the star, while the value of 5780 used in the above fit formulae corresponds to the effective temperature T of the Sun. The coefficients in equations (3) and (4) refer to a planet of one Earth mass but similar functions are available for different terrestrial planet masses. In case stellar luminosities have not been observed directly, one can use stellar radii R * and effective temperatures T ef f instead: where R is the radius of the Sun.

The Permanently Habitable Zone
We will now use the above concepts to determine the permanently habitable zone. In order to find the borders of the region wherein the planet stays always within habitable insolation limits, i.e. the PHZ, we need to find the effective insolation extrema a circumbinary planet is likely to encounter. In hierarchical systems of two stars and a circumbinary planet, the planetary semi-major axis remains practically constant over time. In addition, if a system is coplanar, the time evolution of the eccentricity vectors determines the geometric configuration at any given moment. Assuming furthermore that the gravitational effect of the planet on the stellar binary is negligible, maximum and minimum insolation configurations are determined through the maximum planetary eccentricity e max p . We use the expressions derived in [14] that allow us to calculate to e max p as a function of initial conditions and system parameters. The apocenter and pericenter distance between the planet and the barycenter of the binary star can in turn be expressed through the maximum eccentricity as well as the semi-major axis of the planetary orbit a p , i.e. Q p = a p (1 + e max p ) and q p = a p (1 − e max p ). Figure 1 is a schematic representation of our system, meant to help the reader visualize configurations that lead to insolation extrema. The planet receives maximum insolation when it is at pericenter with respect to the binary barycenter and closer to the brighter star, i.e. when the angle φ between the line that connects the stars and the line that connects the planet to the binary barycenter is φ = 0 • . Thus, the three bodies are aligned and the mathematical expression that determines the inner edge of the PHZ is P HZ I : where the subscript p refers to the planet and the subscript b refers to the binary star orbit. Note that we have normalized the stellar luminosities with their respective spectral weights in order to account for the effect of each star's individual spectral energy distribution on the planetary climate. The minimum insolation configuration is not as straight forward to determine. If one star is substantially more massive and luminous than the other a minimum can be reached when the brighter star is farthest from the planet and the planet is at apocenter with respect to the barycenter of the binary star. This is the case when φ = 180 • which corresponds to the following condition for the outer edge of the PHZ P HZ O : If we have two stars of equal mass, then the minimum radiation configuration is reached when φ = 90 • . Hence, the minimum insolation condition in this case is When the mass ratio of the binary is close to -but not exactly -equal to one (with m 1 > m 2 ) and the planet is located at a suitable distance, the minimum insolation geometry occurs between φ = 90 • and φ = 180 • . Figure 2 illustrates this behavior for a variety of binary mass ratios. We note, however that the difference in insolation received by the planet between said configurations is minute. Hence, we limit our approach to comparing the perpendicular and straight line configurations and use the one that provides the smallest value. Consequently the minimum insolation condition for the outer border of the PHZ reads One can find the numeric values for the borders of the PHZ by solving equations (5) through (8) for a p .
We would like to point out here that for all the above insolation extrema configurations we have assumed that the stars are point masses. In reality, however, the stars have finite sizes. Depending on the distance between the two stars and on the distance between the planet and the binary, it is possible that when the three bodies are aligned the planet may receive reduced insolation due to the eclipsed star. In such a scenario, the minimum insolation configuration (panel (c) in Figure 1) will still be valid. On the other hand, the maximum insolation configuration (panel (d) of Figure 1) would provide a smaller insolation value. [45] have shown, however, that eclipses have little effect on the overall stability of the climate. But even if the insolation change is considerable, as soon as the three bodies get out of alignment, the planet will then receive near maximum insolation. Hence, equation (5) constitutes a reasonable approximation for our purposes.

The Averaged Habitable Zone
The averaged habitable zone is the region around a binary star where a planet remains habitable inspite of variations in irradiation. That is, as long as the insolation average is compatible with habitable limits a planet with a high climate inertia can remain potentially habitable. The averaged over time radiation that a planet receives when orbiting a single star is where L is the stellar luminosity, P the period of the planetary motion, n p is the mean motion of the planet, f p is its true anomaly and r p the distance between the source of radiation and the planet; we made use of the well known relation r 2 p df p /dt = n p a 2 p 1 − e 2 p . We can now extend this simple relation to circumbinary orbits. Assuming that the distance between the two stars is small compared to the distance between the planet and the binary star barycenter, i.e. the orbital period of the binary pair is much smaller than that of the planet, we can approximate the average over time insolation around the binary by placing the two stars at their barycenter and make use of equation (9). In other words we approximate the three body system as a two-body problem with a central "hybrid star" that has the combined mass, luminosity and spectral energy distribution of the stellar binary. This leads to where e 2 p is the averaged square planetary eccentricity over time and initial angles as given in [14]. Under those assumptions the borders of the AHZ are defined through the following inequalities and Note that e 2 p generally depends on a p in a non-trivial way. The full expression for e 2 p along with the expression for e max p is provided in the Appendix. Once more, solving equations (11) and (12) for a p yields the numeric values for the borders of the AHZ.

The Extended Habitable Zone
The definition of the extended habitable zone in section 2 translates into the following equations: The standard deviation σ can be found via the insolation variance We already have an expression for S from equation (10), but we are yet to find S 2 . For a planet around a single star, S 2 is Using the same approach that lead to equation (10), namely combining the stellar binary into a "hybrid star" we can construct a closed analytic expression for S 2 of a circumbinary planet, namely Combining equations (10) and (16)  System Kepler - The inner border of the EHZ is then defined through while the outer border is defined via Numerical values for the inner and outer EHZ borders can be found via solving the above equations for a p .

APPLICATION TO KEPLER SYSTEMS WITH KNOWN CIRCUMBINARY PLANETS
In order to demonstrate the merit of dynamically informed habitable zones, we apply our method to the Kepler circumbinary planets. For that purpose we select Kepler-16 [4], Kepler-34 and Kepler-35 [43], Kepler-38 [31], Kepler-64 [38], Kepler-413 [27], Kepler-453 [44], Kepler-1647 [28] and Kepler-1661 [40]. Kepler-47 [32] has three planets and is, therefore, beyond the limits of our model. The masses of Kepler-38b, Kepler-64b, Kepler-453b are not well defined. For those cases, we have decided to use the upper limit provided in the relevant publications. The parameters necessary for our calculations of the selected systems can be found in Tables 1 (binary star) and 2 (planet).
First, we calculate all dynamically informed habitable zones assuming no giant planets are present in the systems. That provides us with a first idea of the location of the habitable zones and how the presence of the  [18]. Finally, the column CHZ-K13 provides values for our classical zone using [25] as that specific version of the climate model is also used by [18].  In both stages, we allowed the eccentricities of the binary and the existing planet to vary. That way we get a better picture of the effect of orbital eccentricity, an important quantity that regulates distances between bodies, on the extent of habitable zones in the system. In order to simplify the complex dynamics in the presence of the giant planet, we acknowledge the double hierarchical structure of the problem that allows us to consider the binary as one massive body located at the barycenter. The two stars at their barycenter are considered as one body of mass m b = m 1 + m 2 , therefore reducing the four-body problem to a three-body one. To describe the dynamical evolution of such a system we made use of the relevant equations in [14] when a gp < a p and [7,9,10] when a gp > a p , where a gp is the semi-major axis of the orbit of the giant planet. The stability of that particular kind of triple system was assessed using the criterion developed in [33]. The corresponding empirical formula is based on numerical simulations of a star and two planets. The planet-star mass ratios investigated in [14] ranged from 10 −4 to 10 −2 , the ratio of the planetary semi-major axes was in the interval [3,10], while the planetary eccentricities took values in the interval [0, 0.9].

Frontiers
A terrestrial planet in a star-planet-planet system is stable against either ejections or collisions with the central object when Here, m gp is the mass of the giant planet and e gp its orbital eccentricity. The above expression is valid when the giant planet is closer to the binary than the terrestrial planet. If the giant planet orbits externally to the terrestrial planet, then the criterion becomes Strictly speaking, the planetary systems investigated here do not fall in the planet-star mass ratios investigated in [14]. The planet to planet mass ratio, however, remains within those limits which is ultimately more important for the validity of the stability criterion at hand. This hypothesis can be supported by the reasonable agreement between the stability estimates given by equations (20) (21) and those in [13], where hierarchical triple systems with a wide range of masses and on initially circular orbits where integrated numerically and their stability limit was determined.
We now proceed to calculating habitable zones for the aforementioned Kepler systems. In a first step, we ignore the presence of residing giant planets in the respective systems. Considering a potential terrestrial planet orbiting the stellar binary, we can determine the borders of the classical and dynamically informed habitable zones for all the Kepler systems under investigation. As we can see from the left column plots of Figures 3, 4 and 5, all systems but one have well defined habitable zones and would be capable of hosting a broad range of potentially habitable worlds. The only exception to that is Kepler-16, where more than half of the habitable zone is truncated due to dynamical instability arising from the gravitational perturbations of the stellar binary. A terrestrial planet could only survive near the outer border of the habitable zone with the additional requirement that it does not have a low climate inertia since the PHZ has been completely eliminated.
When we add the giant planet in our model (right column plots of Figures 3, 4 and 5), we notice that in Kepler-16 and Kepler-1647 the added gravitational effect of the giant planet renders the entire habitable zone dynamically unstable. Regarding Kepler-16, this is in agreement with other studies such as [35,29]. Note that these results do not exclude the presence of potentially habitable moons around Kepler-16b. In the Kepler-1647 system, the giant planet is located near the middle of the habitable zone. Although Kepler-453b and Kepler-1661b are located inside the classical habitable zones like the other two giant planets mentioned above, they allow a terrestrial planet to exist near the outer border of the habitable zone and they may even allow that planet to have a partial PHZ. This is because of the relatively small masses (∼Neptune mass) and eccentricities of those planets. In contrast, Kepler-1647b is 1.5M J and resides near the center of the habitable zone.
Regarding the remaining systems, the entire classical habitable zone is essentially dynamically stable. The difference between dynamically informed habitable zones with and without the giant planet perturbers shows, however, that the influence of giant planets goes beyond dynamical instability. In Kepler-34, potentially habitable worlds with low climate inertia could only remain habitable in a tiny region centered around 2.17 au. This can be seen from a comparison between the extent of the PHZ and the black vertical lines in the middle row right panel of Figure 3. Relatively dry planets with a low concentration of greenhouse gases would fall into this category. On the other hand, Kepler-35, Kepler-38 and Kepler-64 could host potentially habitable worlds with low climate inertia over a significant fraction (more than 75%) of their classical habitable zone (Figures 4 and 5). Finally, the extent of the PHZ of Kepler-413 is around 40% of its classical habitable zone. While not prohibitive, the presence of giant planets in Kepler-34, Kepler-413, Kepler-453 and Kepler-1661 requires additional terrestrial worlds to buffer significant variations in irradiation in order to remain habitable.
For comparison, some results from [18,2,3,41] are also provided. Our results seem to be in good agreement with those of [2,3,41], with our classical habitable zone being a bit wider. In order to compare with [18], we use [25] to calculate our classical habitable zone as that specific climate model was used by [18] for their habitable zone calculations. In this case, however, it is more difficult to make a direct comparison as [18] used time variable habitable zones centered at the primary star. Habitable zones borders for the Kepler systems under investigation are presented in Table 3.

DISCUSSION AND SUMMARY
We present an analytical approach to determine dynamically informed habitable zones in binary star systems with a circumbinary giant planet. The method takes into consideration the orbital evolution of the giant and terrestrial planet as well as different responses of planetary climates to variations in the quantity and spectral energy distribution of incoming radiation. It does not apply, however, during the planet formation stage, when we can have planets migrating due to interactions with the protoplanetary disk or planet-planet scattering events during late stage formation. In addition, we do not consider systems where there is significant tidal interaction between the two stars which may lead to changes in the orbit and rotation rates of the stars, as well as to changes in the emission of XUV radiation that can affect the atmosphere of a planet within the habitable zone [e.g. 37,22].
As the method mainly relies on analytical equations, it can provide a quick assessment of the capability of terrestrial circumbinary planets in complex dynamical environments to retain liquid water on their surface. The construction and comparison of dynamically informed habitable zones, i.e. the PHZ, the EHZ and the AHZ allows us to better understand where potentially habitable worlds with different climate characteristics can exist in binary star systems. The method presented here is very versatile as it has been constructed in such a way that it does not depend on the dynamical model and the insolation limits.
In this work, we investigated the effects of stellar binarity and circumbinary giant planets on the habitable zones of nine systems observed by the Kepler mission. We confirm earlier studies that suggest Kepler-16 is not suitable for hosting a terrestrial planet within its classical habitable zone. The situation is similar for Kepler-1647. In contrast, Kepler-34, Kepler-35, Kepler-38, Kepler-64 and Kepler-413 seemed more promising with Kepler-38 being the best candidate in this respect. Kepler-453 and Kepler-1661 stand between the previous two categories of systems. We find that nearly equal binary mass ratios and small eccentricities of the perturbing bodies provide favorable, from the orbital evolution point of view, conditions for an Earth-like planet to exist in the habitable zone. We show, furthermore, that the presence of a giant planet can have a significant effect on the potential habitability of terrestrial worlds in the same system. We, thus, recommend gravitational perturbations of known giant planets to be taken into account in future studies regarding habitability in binary star systems. where The above equations can be used for nearly coplanar systems and systems that are not close to mean motion resonances. Also, the equations become less reliable as the maximum planetary eccentricity gets larger than 0.2-0.25. Figure 1. A schematic representation of possible binary star -planet geometries representing irradiation minima: a), b), c) and maxima d) for various mass ratios of the binary star (m 2 /m 1 ). A standard massluminosity relation for main sequence stars is assumed. The more massive star is m 1 , the less massive star is m 2 , and the planet is represented by m p . Figure 2. Angle φ min for the minimum planetary insolation configuration against binary mass ratio. We assume a planet on a circular orbit at a distance of 2.37 au from the binary barycenter, while the binary stars are evolved on a circular orbit with a semi-major axis of 1 au. Plots in the left column show the different types of habitable zones without the presence of the known giant planets. The right column includes the influence of the known giant planets. Red colored regions correspond to uninhabitable areas, blue, green, yellow and purple colors denote the PHZ, the EHZ, the AHZ, and unstable areas according to [20] stability criterion, respectively. Violet colored areas mark regions of dynamical instability caused by the giant planet in the system ( [33] dynamical stability criterion). The vertical black lines denote the classical habitable zone limits, while the horizontal white line in the left column plots marks the current eccentricity of the binary star orbit. In the right column graphs, the white line marks the current eccentricity of the giant planet orbit. Finally, the black dot in the right column plots shows the position of the giant planet in the presented parameter space.