AN ANALYTIC METHOD TO DETERMINE HABITABLE ZONES FOR S-TYPE PLANETARY ORBITS IN BINARY STAR SYSTEMS

Fueled by the successes of current transit-observation incentives like Kepler (Welsh et al. 2012; Borucki & Koch 2011) and CoRoT (Baglin et al. 2009; Tingley et al. 2011) the quest for discovering the first Earth twin has lead to a considerable cross-disciplinary interest in the interplay between stellar and planetary properties necessary to produce habitable worlds. Even though opinions differ on what exactly to look for in a system harboring a terrestrial planet in order to declare it "habitable" (see, e.g., Buccino et al. 2006; Kaltenegger et al. 2007; Selsis et al. 2007; Lammer et al. 2009), the classical assumption investigated by Kasting et al. (1993), i.e., the capacity for water to stay liquid on the planet's surface, may still be considered a prerequisite for the development and sustainability of complex life as we know it (Kaltenegger & Sasselov 2011). This entails restrictions on planetological characteristics, such as mass, atmospheric, and bulk composition, and sets limits to the host star's activity as well as radiation properties (Lammer et al. 2009). Dynamical considerations are of equal importance, since changes in orbital stability or extreme variations in insolation due to large planetary eccentricities (e_p > 0.7) may also result in a hostile environment (Williams & Pollard 2002). It is therefore only natural that one would look for a copy of our solar system when searching for habitable worlds. Yet the study of exoplanetary systems so far has clearly shown that a broader perspective is required.

In fact, up to 70% of all stellar systems in our galaxy may not be single-stellar systems but multi-stellar systems (e.g., Kiseleva-Eggleton & Eggleton 2001 and references therein). Together with the approximately 60 planets that have already been discovered in systems harboring two stars (Schneider et al. 2011) this suggests that binary and multiple star systems should not be ignored in the search for habitable worlds.

Investigations of environments that permit planetary formation in binary star systems have progressed rapidly over the last decade (see, e.g., Thebault 2011 and references therein). Even though important questions regarding the early phases of planet formation in binary star systems—especially the transition from planetesimal to planetary embryos—still remain to be answered, late stage formation scenarios for terrestrial planets in such environments are available (Whitmire et al. 1998; Haghighipour & Raymond 2007; Haghighipour et al. 2010). Since previous studies did consider the extent of the classical habitable zone (KHZ) to be purely a function of the primary star's luminosity and spectral type as introduced by Kasting et al. (1993, hereafter KWR93), we aim to refine this definition to encompass the gravitational and radiative influence of a second star.

This article is structured as follows: after a short recapitulation of the main radiative aspects of habitability as defined in KWR93, Section 2 introduces three exemplary binary–planet configurations, which will serve as test-cases for habitability considerations. Section 3 briefly mentions the dynamical requirements that binary–planet configurations have to fulfill in order to ensure system stability. In Section 4, the maximum radiative influence of the second star on a terrestrial planet in the primary's HZ is estimated and compared to actual insolation simulations. The occurring differences are investigated in the following section. Finally, generalized, analytical estimates are developed and compared to numerical simulations in Sections 6 and 7, and the results concerning the behavior of HZs in binary star systems are presented in Section 8. A discussion of the results concludes this article.

Apart from planetological and dynamical aspects, the insolation a terrestrial planet receives from its host star is naturally the main driver determining the extent of the HZ. When considering planets within a binary star system it is therefore important to track the combined radiation of both stars arriving at the planet. KWR93 showed that not only the sheer amount of insolation but also the spectral distribution is essential to estimate limiting values for atmospheric collapse. In order to model the impact of diverse stellar spectral classes on an Earth-like planet's atmosphere, KWR93 introduced the so-called effective radiation values. These measure the relative impact a comparable amount of radiation (e.g., 1360 [W m⁻²]) with different spectral properties has on a planet's atmosphere. Taking the effects of different stellar spectra into account is especially important in binary star systems with different stellar components. Table 1 reproduces the effective radiation values for the inner (runaway greenhouse) and outer (maximum greenhouse) edges of a single star's HZ as given in KWR93. Notice how the onset of runaway greenhouse effects requires almost twice as much radiation for a spectral distribution akin to F0 class stars compared to M0 spectral types. Even though Kaltenegger & Sasselov (2011) assume similar effective radiation values for M and K spectral classes, actual calculations have only been done for F0, G2, and M0 zero-age main-sequence (ZAMS) stars. For this reason, we first investigate the following three stellar configurations:

• 1.  
  G2–M0 μ ≃ 0.3.
• 2.  
  G2–G2 μ ≃ 0.5.
• 3.  
  G2–F0 μ ≃ 0.6.

All stellar components are considered to be ZAMS stars, and μ = m₂/(m₁ + m₂) denotes the binary's mass ratio. A terrestrial planet is orbiting the Sun-like G2 host-star, hereafter referred to as primary. Such binary–planet configurations are considered to be of satellite type (S-Type), i.e., the planet revolves around one star (Dvorak 1984); see Figure 1. In fact, most of the planets in binary systems have been discovered to be of S-Type (Schneider et al. 2011), e.g., the system Gamma Cephei (Hatzes et al. 2003). In the following, we choose binary systems with semimajor axes between 10 and 50 AU for the comparison of numerical and analytical estimates on the extent of the HZ, as for closer binaries the HZs in G2–G2 and G2–F0 configurations are considerably reduced due to dynamical instability, whereas the qualitative differences for results beyond 50 AU are small. However, the methods presented in Section 6 are viable beyond those limits, as long as the assumptions given in Georgakarakos (2003) remain valid.

The terrestrial planet is assumed to be cloudless (Kaltenegger & Sasselov 2011) and orbits the G2 star only, resulting in configurations (1) and (2) to be classified as S-Type I (μ ⩽ 0.5) and (3) as S-Type II (μ > 0.5), respectively (see Figure 1).

Zoom In Zoom Out Reset image size

Binary–planet configurations that lead to highly chaotic orbits and eventual ejection of the planet cannot be considered habitable. Therefore, dynamical stability of the system is a basic requirement for our consideration. In order to assess the dynamical stability of a test planet's orbit within the given binary star systems, we applied results of numerical stability studies by Rabl & Dvorak (1988; hereafter RD88), Pilat-Lohinger & Dvorak (2002; hereafter PLD02), and Holman & Wiegert (1999; hereafter HW99), determining stable zones in a planar, normalized system (binary semimajor axis a_b = 1). In contrast to RD88 and HW99 who classified unstable orbits via ejections of test planets from the system, PLD02 applied the Fast Lyapunov Indicator (FLI) chaos detection method (Froeschlé et al. 1997) implemented in a Bulirsch–Stoer extrapolation code (Bulirsch & Stoer 1964). In Table 2 maximum planetary semimajor axes which still allow for a stable configuration in the normalized setup are shown for different binary orbits and mass ratios. Results gained via FLI by PLD02 are compared to those published in HW99. Even though the critical values are very similar, the onset of dynamical chaos (PLD02) in G2–G2 and G2–F0 configurations appears for smaller planetary semimajor axes than predicted by the ejection criterion in HW99. For the following investigations the more conservative FLI stability estimates by PLD02 are used.

The main question in dealing with habitability in binary star systems is doubtlessly: "How does the second star affect the HZ around the primary?" Let us focus on the dynamical aspects first, as stability is a prerequisite to habitability. The results from the previous section dictate that not all combinations of a binary's semimajor axis a_b and eccentricity e_b are viable, if orbital stability of planets in the primary's HZ is required. Figure 2 (top left) shows quadratic fits of the FLI stability data presented in Table 2. The secondary's allowed eccentricities are higher for planets orbiting the primary near the inner border of the KHZ than for cases where the planet's semimajor axis is close to the KHZ's outer rim. Such dynamical restrictions set limits on how close the secondary can approach the planet, which will in turn limit its insolation on the planet. Applying analytical approximations introduced later in this section Figure 2 (top right) demonstrates that the smallest possible separation between planet and secondary is always larger than the planet's aphelion distance to its primary. This lessens the second star's potential radiative contribution considerably for all but close S-Type II configurations. In the latter case, the enhanced luminosity of the secondary compensates larger distances to the planet. Using the minimum distances presented in Figure 2 (top right), one can estimate that in a G2–G2 S-Type I system with a_b = 20 AU, the secondary's radiative influence on a planet started at the inner edge of the KHZ is in the order of only 10% of the primary's contribution. As a consequence, the secondary was often considered not to have a significant radiative impact on the extent of the primary's HZ (Whitmire et al. 1998; Haghighipour & Raymond 2007; Haghighipour et al. 2010).

Zoom In Zoom Out Reset image size

Such a so-called single-star approach, however, stands in contrast to numerical experiments presented in the left panel of Figure 3. The effective insolation curves shown were generated solving the full Newtonian Three-Body Problem (3BP) numerically via Lie-Series (Hanslmeier & Dvorak 1984; Eggl & Dvorak 2010) and Gauss–Radau (Everhart 1974) integrators, where the actual amount of radiation the planet receives was calculated for each integration step. One can see that a planet started at the inner edge of the KHZ in a G2–G2 S-Type binary experiences an increase of more than 30% in insolation compared to a planet on a circular orbit with corresponding semimajor axis around a single G2 star. Consequently, an Earth-like planet would spend a considerable time outside the classical, circular HZ in a G2–G2 configuration. Is, therefore, the secondary's radiative impact more important than assumed? The fact that the variations in insolation perfectly correlate with the dynamical evolution of the planet's eccentricity (see Figure 3, right panel) permits an alternative explanation: changes in the planetary orbit induced via gravitational perturbation by the second star might be responsible for the increased insolation values. Even though planetary and binary orbits' semimajor axes are not expected to show any secular variation (e.g., Harrington 1968), it is known that even a distant companion would inject some eccentricity into the orbit of two bodies revolving around their common center of mass (e.g., Mazeh & Shaham 1979; Georgakarakos 2002). Elevated planetary eccentricities would entail smaller periastron distances—allowing for increased insolation by the primary.

In order to draw a clearer picture of whether the perturbation-induced rise in the planet's eccentricity or the secondary star's radiative influence are the dominant factors causing the increased insolation onto the planet encountered in Figure 3 (left panel), we will make three assumptions, which our analytical estimates will be based on:

• 1.  
  the binary–planet system is coplanar;
• 2.  
  stellar luminosities L are constant on timescales of the system's secular dynamics; and
• 3.  
  stellar occultation effects are negligible.

Since we assume coplanar orbits, we make extensive use of the analytic expressions in Georgakarakos (2003, 2005) to calculate the time-averaged-squared planetary eccentricity 〈e²_p〉_t. In Appendix B, we derive estimates on the planet's maximum injected eccentricity e^max_p. Unlike the recent ansatz by Giuppone et al. (2011) and earlier by Thébault et al. (2006), where the eccentricity evolution of a planet in a stellar binary was modeled by empirical formulae, Georgakarakos (2003, 2005) derived an entirely analytical formula for the eccentricity of the inner—in this case the planet's—orbit of a hierarchical triple system, which is valid for a wide range of mass ratios and initial conditions. Together with the estimates for e^max_p in Appendix B, the formulae presented in Georgakarakos (2003, 2005) are an analytical extension to the first order with respect to the perturbing mass secular perturbation theory, as given for example in the longstanding work of Heppenheimer (1978).

Assumptions 2 and 3 are reasonable for well-separated binary stars, where both components are on the main sequence. The combination of the dynamical stability results presented in Figure 2 (top left panel) and the analytic expressions for the injected planetary eccentricity allows us to estimate not only the minimum distances between the secondary and the planet as shown in Figure 2 (top right panel), but also the maximum contributions to the planetary insolation from the primary and the secondary star, respectively.

Figure 2 (bottom left panel) shows the primary's effective insolation onto a terrestrial planet during its closest perihelion passage d = a_p(1 − e^max_p). Here, the largest dynamically possible perturbation by the secondary is considered. In order to separate the secondary's radiative and gravitational effects, the second star's radiation has been excluded in this plot. Given a binary semimajor axis of 10 AU, Earth analogs started at the inner edge of the KHZ permit considerably higher binary eccentricities (Figure 2, top left panel) compared to planets near the outer edge. Therefore, higher insolations for planets near the inner edge of the KHZ are expected, since the injected planetary eccentricities are coupled strongly to the binary's eccentricity; see Appendix B. This effect abates around a_b = 20 AU, where the primary's insolation becomes almost equal for planets started near the inner and outer rims of the KHZ. For binary semimajor axes beyond a_b = 20 AU planetary orbits closer to the secondary are more severely perturbed than the ones near the inner edge of KHZ, as the binaries' maximum allowed eccentricities become similar for inner and outer borders of the KHZ and the perturbation is merely dependent on the binary to planet period ratios.

For planets on eccentric orbits in a G2–G2 configuration with a_b = 20 AU, the primary's insolation can increase up to almost 50% compared to circular orbits sharing the same semimajor axis. In contrast, even with the planet being as close to the second star as dynamically possible, the secondary's maximum radiative input only accounts for an additional 10% in the same setup. This can be seen in the bottom right panel of Figure 2. Hence, the primary is the main source of the additional insolation in S-Type I systems, which can only be explained via changes in the planet's eccentricity. Solely for S-Type II systems like G2–F0 binaries, the secondary's radiative contributions to planetary insolation can become comparable to the primary's (cf. Figure 2, bottom panels).

In the previous section, the rise in the planet's eccentricity was identified as a main driver for the strong variance in simulated insolation curves in S-Type I systems. Eccentric planetary orbits entail strong variations in insolation over an orbital period. Nevertheless, Williams & Pollard (2002) concluded that this might not necessarily be prohibitive for habitability. As long as the average insolation values lie within habitable limits and e_p < 0.7, the atmosphere should be able to act as a buffer, preventing immediate carbon freeze-out or instantaneous water evaporation. In order to distinguish between cases where the planet remains within the radiative boundaries of the KHZ for all times, configurations where the planet is "on average" within the HZ, and non-habitability, we introduce the following three categories.

Permanently Habitable Zone (PHZ). The PHZ is the region where a planet always stays within the insolation limits (A, B) of the corresponding KHZ, i.e.,

$$\begin{equation} A \ge S_{{\rm eff}} \ge B, \end{equation} \tag{ 1 }$$

where A denotes the inner and B the outer effective radiation limit for a given spectral type; see Table 1.

Extended Habitable Zone (EHZ). In contrast to the PHZ—where the planet stays within the KHZ for all times—parts of the planetary orbit lie outside the HZ due to, e.g., the planet's eccentricity. Yet the binary–planet configuration is still considered to be habitable when most of its orbit remains inside the HZ boundaries:

$$\begin{equation} \left\langle S_{{\rm eff}}\right\rangle _t +\sigma \le A \wedge \left\langle S_{{\rm eff}}\right\rangle _t -\sigma \ge B , \end{equation} \tag{ 2 }$$

where 〈S_eff〉_t denotes the time-averaged effective insolation from both stars and σ² is the effective insolation variance.

Averaged Habitable Zone (AHZ). Following the argument of Williams & Pollard (2002), this category encompasses all configurations which allow for the planet's time-averaged effective insolation to be within the limits of the KHZ:

$$\begin{equation} A \ge \left\langle S_{{\rm eff}} \right\rangle _t \ge B. \end{equation} \tag{ 3 }$$

We now propose analytical estimates to achieve a classification of planetary habitability as was suggested in the previous section. The aim is to circumvent time-consuming numerical integrations when global parameter scans are required to check systems for possible habitability. Even though the following analytical estimates are presented utilizing the S_eff values developed in KWR93, more advanced atmospheric models for exoplanets can easily be introduced by exchanging the effective insolation values A and B.

6.1. Estimates for the PHZ

Let the second star move on a fixed Keplerian orbit with semimajor axes a_b and eccentricity e_b. Accordingly, the planet's orbit has the semimajor axis a_p and acquires a maximum eccentricity of e^max_p due to the secondary's gravitational perturbations (cf. Figure 3, right panel). This permits us to estimate the maximum and minimum insolation conditions for the planet to permanently remain within the KHZ.

Zoom In Zoom Out Reset image size

Insolation minimum condition. Both planet and secondary are assumed to be in apocenter position and opposition. The additional normalization of the stellar luminosities per solid angle (L = L_bol/(4π)) with regard to the respective outer insolation limits for each star (B₁, B₂) ensures that different spectral properties are taken into account:

$$\begin{equation*} 1 \le \frac{L_1}{B_1} \big(a_p \big(1+e^{{\rm max}}_p\big)\big)^{-2} + \frac{L_2}{B_2} \big(a_b (1+e_b)+a_p\big(1+e^{{\rm max}}_p\big)\big)^{-2}. \end{equation*}$$

Insolation maximum condition. Again the luminosities are normalized, but this time with regard to the inner insolation limits (A₁, A₂). Since we consider S-Type I and II systems, it is possible that the secondary at pericenter may produce a higher insolation on the planet than the primary star. If this is the case, the maximum insolation configuration will have the planet at apocenter with regard to the primary and the secondary at pericenter:

$$\begin{eqnarray*} &&1 \ge {\rm max} \left\lbrace \begin{array}{@{}l@{}} \dsty \frac{L_1}{ A_1} \big(a_p \big(1-e^{{\rm max}}_p\big)\big)^{-2} \\ \dsty \quad +\,\frac{L_2}{ A_2} \big(a_b (1-e_b)-a_p \big(1-e^{{\rm max}}_p\big)\big)^{-2}, \vspace{8.5359pt} \\ \dsty \frac{L_1}{ A_1} \big(a_p\big(1+ e^{{\rm max}}_p\big)\big)^{-2}\\ \dsty \quad +\,\frac{L_2}{ A_2} \big(a_b (1-e_b)-a_p \big(1+e^{{\rm max}}_p\big)\big)^{-2}. \end{array} \right. \end{eqnarray*}$$

6.2. Estimates for the AHZ

The combined stellar insolation S_tot on the planet is, of course, a function of time. In order to calculate time-averaged insolation values, we will use that

$$\begin{equation} \left\langle S_{{\rm tot}} \right\rangle _t = \left\langle S_{1} \right\rangle _t + \left\langle S_{2} \right\rangle _t, \end{equation} \tag{ 4 }$$

where 〈S₁〉_t is the time average of the planetary insolation due to its host star and 〈S₂〉_t is the time-averaged contribution of the second star. Let us focus on the two-body problem planet–host star first. Here, 〈S₁〉_t = 〈L₁/δ²(t)〉_t where δ(t) denotes the scalar distance between planet and primary. The average insolation a planet on an unperturbed Keplerian orbit experiences can be calculated using the planet's angular momentum $h=\delta ^2\dot{f}$, f being the true anomaly:

$$\begin{eqnarray} \langle S_1 \rangle _t &=&\frac{L_1}{P}\int _0^P \frac{1}{\delta ^2(t)}dt \\ \nonumber &=&\frac{L_1 n}{2\pi } \int _0^{2\pi } \frac{1}{h}df \\ \nonumber &=& \frac{L_1 n}{ h}. \end{eqnarray} \tag{ 5 }$$

This expression states that the time-averaged insolation a planet receives over one orbit depends only on the star's luminosity, the planet's mean motion n, and its orbital angular momentum (see, e.g., Seager 2010, p. 18). Now let us construct a circular orbit so that h_circle = h_ellipse. A planet moving on such an orbit with "equivalence radius" r = a_p(1 − 〈e²_p〉_t) will experience the same amount of insolation per unit time as a planet on an elliptic orbit sharing the same angular momentum. The advantage of considering "equivalent circular orbits" is that the insolation remains constant for any orbital position of the planet. Of course, the reduction of elliptic orbits to circular orbits with common angular momentum will decrease the planet's orbital period, as r ⩽ a and r = a only if 〈e²_p〉_t = 0. However, as we have chosen the "equivalent circular orbit" to have the same angular momentum, it will share the same constant rate of change of insolation with the true orbit. Therefore, even an average over the longer period of the elliptic orbit will still yield valid results.

We will now apply the same line of argumentation to construct the equivalence radius R for the secondary. This allows us to extend the pseudostatic radiation environment to include the average radiative influence of the secondary. Since the distinction between averaged and initial eccentricity is small for the secondary—its secular variance is negligible for the cases investigated—we can safely use the secondary's initial eccentricity to calculate R. The suggested configuration is shown in Figure 4. The secondary circles the primary in a fictitious orbit with equivalence radius R, and the planet orbits the primary in a circle with radius r.

Zoom In Zoom Out Reset image size

The insolation on the planet depends on the relative distances of the planet to both stars. In order to estimate the insolation on the planet caused by the secondary (S₂), we simply apply law of cosine, and average over all possible geometric configurations (Figure 4):

$$\begin{eqnarray} \rho ^2 &= &R^2+r^2-2 r R \cos (\phi) \nonumber \\ \left\langle S_2 \right\rangle _t &= & \left\langle \frac{L_2}{ \rho ^2(\phi)} \right\rangle _\phi \\ &= & \frac{L_2}{2\pi }\int _0^{2\pi }\frac{1}{R^2+r^2-2 r R \,\cos (\phi)}\,d\phi \nonumber \\ &= &\frac{L_2}{(R^2-r^2)}\quad {\rm if} \; R>r, \nonumber \end{eqnarray} \tag{ 6 }$$

where ρ is the planet's distance to the secondary and ϕ denotes the angle between the distance vectors of planet and secondary to the host star. Since we do not allow for orbit crossings of the planet and the secondary, R > r will always hold. Consequently, the total, time-averaged insolation onto the planet is given by

$$\begin{equation} \left\langle S_{{\rm tot}} \right\rangle _t = \frac{L_1}{r^2}+\frac{L_2}{(R^2-r^2)}. \end{equation} \tag{ 7 }$$

Equation (7) does not yet take the different spectral properties of the binary's components into account. Therefore, we note the following conditions for the planet's averaged insolation being within habitable limits.

Average insolation minimum condition:

$$\begin{eqnarray*} &&1 \le \langle S_{{\rm eff}, B} \rangle _t = \frac{L_1}{ B_1}\frac{1}{ r^2}+\frac{L_2}{ B_2}\frac{1}{R^2-r^2}. \end{eqnarray*}$$

Average insolation maximum condition:

$$\begin{eqnarray*} &&1 \ge \langle S_{{\rm eff}, A} \rangle _t = \frac{L_1}{ A_1}\frac{1}{r^2}+\frac{L_2}{ A_2}\frac{1}{R^2-r^2}. \end{eqnarray*}$$

Here, the indices 1, 2 indicate the respective star's KHZ boundary values A and B, R = a_b(1 − e²_b), and r = a_p(1 − 〈e²_p〉_t), where the averaged-squared planetary eccentricity was calculated following Georgakarakos (2005).

6.3. Estimates for the EHZ

Having derived the insolation time averages in the previous section, the expected insolation variance (σ²) still remains to be determined:

$$\begin{equation} \sigma ^2=\big\langle S_{{\rm tot}}^2\big\rangle _t- \big\langle S_{{\rm tot}}\big\rangle _t^2. \end{equation} \tag{ 8 }$$

Using Equation (7), and analytic estimates for 〈S²_tot〉_t, the effective insolation variance can be calculated as follows (see Appendix A):

$$\begin{eqnarray} \sigma ^2_{X} &=&\frac{L_1^2}{X_1^2r^4}\big(-1+3\big\langle e_p^2\big\rangle -3\big\langle e_p^2\big\rangle ^2+\big\langle e_p^2\big\rangle ^3\big) \nonumber \\ &&+\frac{L_1^2}{X_1^2r^4}\sqrt{1-\big\langle e_p^2\big\rangle }\left(1-\frac{\big\langle e_p^2\big\rangle }{2} - \frac{\big\langle e_p^2\big\rangle ^2}{2} \right) \nonumber \\ &&-\frac{2L_1L_2}{X_1X_2(r^4-r^2R^2)}\big(1-\sqrt{1-\big\langle e_p^2\big\rangle }\big(1+\big\langle e_p^2\big\rangle \big)\big) \nonumber \\ &&-\frac{2L_2^2r^2}{X_2^2(r^2-R^2)^3}, \end{eqnarray} \tag{ 9 }$$

where X_i ∈ {A_i, B_i} and the index i denotes the respective star. The minimum and maximum conditions for a planet to be within the EHZ are given as follows.

Extended insolation minimum condition:

$$\begin{eqnarray*} &&1 \le \langle S_{{\rm eff}, B} \rangle _t -\sigma _B. \end{eqnarray*}$$

Extended insolation maximum condition:

$$\begin{eqnarray*} &&1 \ge \langle S_{{\rm eff}, A} \rangle _t + \sigma _A. \end{eqnarray*}$$

In order to test the reliability of the analytical estimates presented in Section 6, we used high-precision numerical integration methods based on Gauss–Radau quadrature (Everhart 1974) and Lie-Series (Hanslmeier & Dvorak 1984; Eggl & Dvorak 2010) to determine the actual positions of both stars with respect to the Earth-like planet. Assuming that the stellar luminosities will not change significantly on the timescale of the planet's secular period in eccentricity, good approximations of insolation patterns can be obtained. Figure 6 shows a comparison of analytic habitability classifications versus results gained via numerical orbit integration and direct insolation calculation. The setup consists of a terrestrial planet (1 M_⊕) in S-Type orbit around a G2 host star, with three different spectral types as secondary: F0 (top), G2 (middle), and M0 (bottom). The terrestrial planet was started on circular orbits with semimajor axes between 0.6 AU ⩽ a_p ⩽ 2 AU with the secondaries' semimajor axes being a_b = 50 AU. The time span of the numerical integrations encompassed at least two secular periods in the planet's eccentricity. It is evident that for small binary eccentricities, all three types of HZ coincide well with the borders defined in KWR93 indicated by the vertical lines at 0.84 and 1.67 AU. For e_b > 0.1, however, a splitting into the HZ categories defined in Sections 5 and 6 becomes eminent. The PHZ (black; blue in the online journal) shrinks considerably with growing eccentricity of the binary's orbit. This is due to the perturbation-induced elevation of the planet's eccentricity. In contrast, the region defined in KWR93 is best approximated by the AHZ (light gray; yellow in the online journal), which remains virtually unaffected by the secondary's eccentricity. In this setup all analytically calculated HZs are in excellent agreement with the numerical ones. Only close to the stability limit (shaded region; purple in the online journal) the correspondence between simulation and analytical estimates decreases. This can be seen more clearly when the secondary's influence becomes stronger, e.g., in the cases of a_b = 10 AU; see Figure 7. In general the analytical approach is producing more conservative results compared to the numerical data (cf. Figures 5 and 7). In order to determine whether these deviations in the determination of the PHZ are due to

• 1.  
  inaccurate analytical estimates for e^max_p,
• 2.  
  insufficient time resolution in numerical simulations with regard to determining e^max_p values, or
• 3.  
  insufficient total integration time to reach minimum and maximum insolation conditions,

Zoom In Zoom Out Reset image size

we constructed semi-analytical PHZs using numerically determined e^max_p values in the analytic equations to determine the PHZ presented in Section 6.1. The borders of the semi-analytically derived PHZs are depicted as white dash-dotted lines in Figures 5 and 7. As they are nearly identical with the fully analytic estimates, we can exclude (1) and (2), which would have led to significantly different results for semi-analytic and analytic approaches. Therefore, (3) seems most likely to cause the differences between the numerical and analytical PHZs, since encountering an exact alignment of planetary aphelion and secondary perihelion at the moment where e_p = e^max_p may take far longer than two secular periods. As the computational efforts required to ensure that a simulation's time resolution as well as total integration time is sufficient to identify the correct PHZ boundaries are enormous, the necessity to have analytical methods at hand becomes evident.

As far as EHZ and AHZ regions are concerned, clear differences can be seen in high perturbation environments close to the transition to instability (Figures 5 and 7, shaded regions). In these cases, the authors favor the numerical results, as single configurations are not critical for the more statistically oriented measures.

A quantitative overview of the correspondence between numerical and analytical results for all system configurations investigated is given in Table 3. Here, similar maps as presented in Figures 6 and 7 were generated with a resolution of Δa_p = 0.01 AU, and Δe_b = 0.01 and evaluated statistically. In spite of their shortcomings in determining the PHZ, the numerical classification results have been used as reference values and are compared to the analytical estimates given in Section 6. The total correspondence percentages are calculated as the number of all orbits below the stability limit that were classified identically via numerics and analytics, divided by the total number of orbits simulated. The number of orbits classified as PHZs analytically divided by the number of orbits classified as PHZs numerically yields the percentage of PHZs, etc. Table 3 shows that the global correspondence between both methods is quite convincing, which can be considered a strong indicator that the behavior of the respective HZs is modeled correctly.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Also, Figures 5–7 indicate that most of the significant deviations between numerical and analytical results occur near the border of orbital instability, especially for high mass and small period ratios. This might be expected for AHZs and EHZs given the approximations involved in determining the analytic estimates. In the case of PHZs, however, semi-analytical results suggest caution in using simulation outcomes as reference values.

Diverse trends in the behavior of the different types of HZs can be seen in Figures 6 and 7. While the AHZ is almost independent of the binary's eccentricity and coincides well with the KHZ by KWR93 for distant stellar companions, the PHZ and EHZ shrink with higher binary eccentricities. The fact that the PHZ and EHZ contract around the center of the KHZ emphasizes the importance of the changes in the planet's eccentricity, as the secondary's radiation alone could only account for one-sided features toward the outer edge of the KHZ. The PHZs seem most affected by strong perturbations and shrink to almost half their size before the systems investigated become unstable. Interestingly, in close S-Type II binaries with low eccentricity the extent of the EHZs and AHZs can reach beyond the predicted values by KWR93. This can be seen in Figure 5, where a zoom on the outer border of the KHZ in a G2–F0 configuration is shown. Here, an extension toward the second star of about 0.1 AU seems possible for the system's AHZ if the binary's eccentricity remains below the system's stability limit of e_b ≃ 0.2.

In this work, the impact of the second star on habitable zones in S-Type binary star systems with different stellar constituents has been investigated analytically as well as numerically. The radiative contribution of the secondary on a terrestrial planet is negligible in all but S-Type II systems, if orbital stability of the planet is required. The gravitational influence of the second star, on the other hand, perturbs in the planet's eccentricity, which in turn can lead to substantial changes in planetary insolation. Therefore, the secondary indeed has to be taken into account when calculating the extent of habitable zones.

Our analytical estimates for planetary eccentricities in binary star systems introduced in Appendix B are an extension to secular perturbation theory as used in, e.g., Heppenheimer (1978). Together with methods presented in Section 6 suggesting an analytic determination of habitable zones, they allow us to paint a global picture of habitability in S-Type binary star systems without having to rely on time consuming numerical orbit integrations. Our approach is quite flexible in the sense that different planetary atmospheric models and average stellar luminosities can be integrated via adaption of the S_eff values. Thus, the formulae presented in this article grant access to calculating habitable zones for a large set of possible binary–planet configurations.

For the three stellar configurations investigated it could be shown that the permanently habitable zone, i.e., the zone where the planet never exceeds the classical insolation limits for habitability shrinks considerably with the binary's eccentricity. If one considers average insolation values only, the extent of the average habitable zone coincides well with predictions by Kasting et al. (1993) for wide binaries, whereas a significant extension toward the secondary is possible for close, eccentric binary systems. The overall correspondence between numerical and analytical results presented is excellent, as 93%–99% of all investigated orbits were classified identically. The computational efforts required to calculate the true extent of permanently habitable zones numerically, however, can be enormous and might in fact be prohibitive in some cases. In contrast, the analytical method presented offers immediate, reliable estimates.

A more careful approach than the one proposed in this work is advisable when multi-planetary systems, systems close to the stability limit, or resonant configurations are being investigated. In a next step we plan to extend our classification methods to mutually inclined systems.

The authors acknowledge the support of FWF projects AS11608-N16 (E.P.-L., S.E.), P20216-N16 (S.E., M.G., and E.P.-L.), and P22603-N16 (E.P.-L. and B.F.). S.E. acknowledges the support of University of Vienna's Forschungsstipendium 2012.

The maximum possible eccentricity a terrestrial planet's orbit can acquire in an S-Type setup is composed of

$$\begin{equation} e^{{\rm max}}_p=e^{{\rm sp}}_p+e^{{\rm sec}}_p, \end{equation} \tag{ B1 }$$

where e^sp_p denotes amplitude of short-period terms and e^sec_p the planetary eccentricity's secular amplitude. Using expressions derived in Georgakarakos (2003) via the Laplace–Runge–Lenz vector, we can estimate the maximum short-period contributions for the planet's eccentricity by taking only terms including the dominant frequencies into account. The amplitude of the secular part of the planet's eccentricity is given by 2C/(B − A) (Georgakarakos 2003) resulting in the following expressions:

$$\begin{eqnarray} e^{{\rm sp}}_p&=&\alpha \Bigg(\frac{15}{64}\frac{\beta }{X^{5/3}}\frac{\big(4+11 e_b^2\big)}{\big(1-e_b^2\big)^{5/2}}+\frac{11}{4}\frac{1}{X^{2}}\frac{(1+e_b)^3}{\big(1-e_b^2\big)^3}\nonumber\\ &&+\frac{3}{4}\frac{1}{X^{3}}\frac{(1+e_b)^4(6+11 e_b)}{\big(1-e_b^2\big)^{9/2}}\Bigg) \end{eqnarray} \tag{ B2 }$$

$$\begin{eqnarray} e^{{\rm sec}}_p&=& e_b\beta \Bigg(\frac{5}{4}\frac{\alpha }{X^{1/3}}\frac{3+2e_b^2}{\big(1-e_b^2\big)^{1/2}} -\frac{2}{5}\gamma X^{1/3}\big(1-e_b^2\big)^{1/2}\nonumber\\ && +\frac{2}{5} X^{2/3}\big(1-e_b^2\big)\Bigg)^{-1}. \end{eqnarray} \tag{ B3 }$$

The mass parameters α, β, and γ are defined as

$$\begin{eqnarray*} &&\alpha =\frac{m_2}{M}\quad \beta =\frac{m_1-m_p}{(m_1+m_p)^{2/3} M^{1/3}} \quad \gamma =\frac{m_1 m_p M^{1/3}}{m_2(m_1+m_p)^{4/3}}, \end{eqnarray*}$$

where m_p is the planetary mass, and m₁ and m₂ are the stellar masses of primary and secondary, respectively.

M = m₁ + m₂ + m_p is the total mass of the system. Finally, X denotes the secondary-to-planet period ratio:

$$\begin{eqnarray*} &&X=\frac{P_b}{P_p}=\left(\frac{m_1+m_p}{M}\right)^{1/2}\left(\frac{a_b}{a_p}\right)^{3/2}. \end{eqnarray*}$$