CONSTRAINING THE RADIATION AND PLASMA ENVIRONMENT OF THE KEPLER CIRCUMBINARY HABITABLE-ZONE PLANETS

The outstanding work of the Kepler Observatory circumbinary team has resulted in the discovery of 10 confirmed planets around eight moderately separated binaries (all have binary periods in the range of 8–60 days). These are Kepler-16 (Doyle et al. 2011), Kepler-34 and Kepler-35 (Welsh et al. 2012), Kepler-38 (Orosz et al. 2012a), Kepler-47 (Orosz et al. 2012b), Kepler-64 (Kostov et al. 2013; Schwamb et al. 2013), Kepler-413 (Kostov et al. 2014), Kepler-47 (Hinse et al. 2015), and Kepler-453 (Welsh et al. 2015). Interestingly, three out of eight of these systems (38%), namely, Kepler-16, Kepler-47, and Kepler-453, host giant planets in the circumbinary habitable zone (BHZ)⁴ .

The study of conditions that planets face while orbiting a binary system has attracted the attention of the exoplanetary community. Several authors have studied the stability of their orbits (Dvorak 1986; Holman & Wiegert 1999) and the kind and level of insolation experienced by planets while illuminated by two stars (Harrington 1977; Haghighipour & Kaltenegger 2013). The conditions for the formation and migration of these planets have been studied by Pierens & Nelson (2008), Kley & Haghighipour (2014), and Kley & Haghighipour (2015). Efforts are being made to constrain the magnetic and plasma environment of circumbinary planets (Mason & Clark 2012; Mason et al. 2013; Sanz-Forcada et al. 2014), finding that binaries could potentially pose severe restrictions on the habitability of Earth-like planets and exomoons. We now know that stellar aggression, in the form of high-energy radiation and coronal-activity-related mass loss, would have an effect on the evolution of planetary atmospheres, especially their capacity to retain water (see, e.g., Lammer 2013).

In Mason et al. (2013) we proposed a mechanism for which some BHZ planets experience Earth-like or lower levels of stellar aggression and thereby have the potential to retain water. Tidal torquing of stellar rotation, especially for certain binary periods and initial eccentricities, aids in the reduction of stellar activity. The reduction in intensity of rotationally driven dynamo action results in a decrease in stellar coronal activity, potentially promoting life on planets. We refer to the sum of these and other beneficial effects for life on planets around binaries as the Binary Habitability Mechanism (BHM).

Recent observations lend support to the tidal astrophysics of BHM. For example, Wasp-18, a 0.63 Gyr old F6 star, is orbited by a nearby giant planet, with a 0.94-day orbital period. The star shows a low level of activity, normally expected for a much older star (Pillitteri et al. 2014). This effect has been described as "premature aging." In another example, the interaction of HD 189733 with a close-in planet has resulted in the opposite effect on a much larger timescale (Wolk et al. 2011). HD 189733, a relatively old star, displays anomalously high activity levels, lending support to our prediction that early tidal rotational torquing could result in a high level of activity. We call the latter phenomenon the "forever young" effect, in reference to the fact that tidal interaction of a star with a stellar or substellar companion could make the star very active (to appear young) for a longer time than expected from single-star evolution (Mason et al. 2013; Sanz-Forcada et al. 2014).

Although the specific physical mechanisms operating in the case of Wasp-18 and HD 189733 may be different from those predicted in stellar binary systems, these discoveries confirm the necessity of including tidal interaction in order to assess stellar activity levels of stars with companions and its effect on their planets.

In previous work, we used first-order estimation of tidal synchronization times (Mason et al. 2013). More recently (Mason et al. 2015), we explored habitable niches around hypothetical circumbinary environments, using a basic model for stellar rotation evolution including stellar tides but applied only to the early phases of main-sequence evolution.

Here our models are refined by (1) extending the time frame to the critical pre-main-sequence (PMS) phase where rotational and activity evolution is much richer; (2) applying a more physically motivated model for stellar rotational evolution, including but not restricted to angular momentum (AM) transport inside the star; and (3) estimating mass loss and X-ray emission using the latest solar-inspired models, relating basic stellar properties and rotation to chromospheric and coronal activity (Cranmer et al. 2013).

For these three critical refinements, we are applying models that have been developed and tested in the case of single stars. We are assuming here that the same models can be applied to describe low-mass stars in moderately separated binaries (orbital periods larger than 8 days and separations larger than several tens of stellar radii). We will justify and test this assumption later on in the paper.

We apply our improved and extended model to the interesting and well-characterized Kepler binaries with BHZ planets. The planets currently known in those systems are Neptune to Saturn sized, so no surface habitability is expected. However, it is reasonable to suspect that Earth-like planets may exist in the BHZ of similar binaries (Eggl et al. 2013). Additionally, since circumbinary giant planets appear to be common, some of them may harbor potentially habitable exomoons (Heller et al. 2014). These planets and putative moons are expected to experience harsh environments from aggressive host stars. Thus, understanding the radiation and plasma environment experienced by these bodies is a key effort to accessing their potential habitability (Zuluaga et al. 2012; Heller & Zuluaga 2013).

The present aim is to constrain the radiation and plasma environment of these HZs and to investigate the effect that these conditions have on the atmospheric evolution of hypothetical terrestrial planets and exomoons.

This model can readily be applied to other Kepler binaries or any other system with well-determined parameters. Future telescopes like PLATO (Rauer et al. 2013) will likely discover circumbinary planets exhibiting some potential for habitability. Our aim here is to inform such searches by showing how a more comprehensive theory of planetary habitability can be applied to circumbinary planets especially in the presence of interacting stars.

In order to provide the reader an opportunity to reproduce these results and to further explore the vast parameter space of circumbinary planets, a Web-based BHM Calculator, BHMcalc, is made available using the link http://bhmcalc.net. In Mason et al. (2015) we introduced the first version of the tool. Here we present a comprehensive version of the calculator, detailing all the new physical effects described in this paper, and providing an improved and user-friendly interface.

Other online tools have recently been made available (Cuntz & Bruntz 2014; Müller & Haghighipour 2014). Like these, BHMcalc provides renditions of the instantaneous and continuous BHZ. In addition, BHMcalc calculates time-dependent stellar and planetary properties and utilizes them to study HZ environments. Most critically, the rotational evolution of the stellar components is derived, allowing for estimates of the resulting stellar mass loss and magnetic activity. Stellar emission is used to evaluate effects experienced by planets. Specifically, we select the integrated X-ray and extreme-ultraviolet radiation (XUV) and stellar wind (SW) fluxes as derived for planets, located over a range of distances from the binary center of mass, as proxies for habitability. For this purpose we estimate the XUV- and SW-induced planetary atmospheric mass loss. In all cases, the evolution of each star is calculated independently and orbital averages are used to calculate the combined effect of the binary on the planet.

This paper is organized as follows: In Section 2, the sample of three Kepler binaries possessing BHZ planets is described. Section 3 presents an efficient model to calculate BHZs and their evolution in time. In Section 4, we revisit the argument for the need to access the radiation and plasma environment of planets as potential drivers of atmospheric evolution. High-energy radiation and winds are key to constraining habitability. The comprehensive model used to calculate the evolution of rotation and activity is described in detail in Section 5. The model is applied to constrain the radiation and plasma environments of the solar system as a reference in Section 6. The binary rotational evolution affects potential circumbinary worlds and so is presented in Section 7, and the radiation and plasma environments are investigated in Section 8. Finally, in Section 9, a summary and the conclusions of this investigation are presented.

It is interesting to note that while none of the currently known transiting circumbinary planets are Earth-like or even super-Earths, this may be due purely to three selection effects. The Kepler telescope is most sensitive to the detection of large planets, orbiting small stars, on short-period orbits. However, the fact that on average circumbinary planets have longer periods than those Kepler planets with single-star hosts is not a selection effect. Short-period circumbinary planets are limited by dynamical stability (Harrington 1977; Dvorak 1986; Holman & Wiegert 1999) and probably also because of formation (Pierens & Nelson 2008; Kley & Haghighipour 2014). Remarkably, the BHZ extends well beyond the orbital stability limit in many cases (see Mason et al. 2015). In addition, residence of planets in the BHZ appears to be common, as 3 of the 10 known transiting planets are located in the BHZ and three binaries out of eight harbor a BHZ planet.

In Table 1, properties of the three transiting BHZ planets and their host stars are listed. Hereafter, we will refer to this as the KBHZ sample. The binaries in the KBHZ sample are among the best-characterized binaries known, mainly owing to the extraordinary observational constraints provided by transiting circumbinary planets.

Table 1.  The KBHZ Sample, Kepler Binary Systems Harboring Planets in the Circumbinary Habitable Zone

Parameter	Kepler-16b	Kepler-47c	Kepler-453b
Stellar Properties
M1 (M⊙)	0.6897 ± 0.0034	1.043 ± 0.055	0.944 ± 0.010
R1 (R⊙)	0.6489 ± 0.0013	0.964 ± 0.017	0.833 ± 0.011
${T}_{{\mathrm{eff}}_{1}}$ (K)	4450 ± 150	5636 ± 100	5527 ± 100
Prot,1 (days)	35.1 ± 1.0	7.775 ± 0.022	20.31 ± 0.47
M2 (M⊙)	0.20255 ± 0.00066	0.362 ± 0.013	0.1951 ± 0.0020
R2 (R⊙)	0.22623 ± 0.00059	0.350 ± 0.006	0.2150 ± 0.0014
${T}_{{\mathrm{eff}}_{2}}$ (K)	⋯	3357 ± 100	3226 ± 100
Prot,2 (days)	⋯	⋯	⋯
[Fe/H] (dex)	−0.20	−0.25	−0.34
Closest Evolutionary Track
${M}_{{\rm{1}}}^{*}\ ({M}_{\odot })$	0.70	1.05	0.90
${M}_{{\rm{2}}}^{*}\ ({M}_{\odot })$	0.20	0.35	0.20
${Z}_{{\rm{1}}}^{*}\ ({M}_{\odot })$	0.010	0.010	0.008
${\rm{\Delta }}{\tau }_{{\rm{MS,1}}}^{*}$ (Gyr)	>13	5.53	10.1
Binary Properties
Tage (Gyr)	2.0–3.0	1.0–5.0	1.5–2.5
Pbin (days)	41.07922 ± 0.00007	7.44837695 ± 0.00000021	27.322037 ± 0.000017
abin (AU)	0.22431 ± 0.00035	0.0836 ± 0.0014	0.18539 ± 0.00066
ebin	0.15944 ± 0.00062	0.0234 ± 0.0010	0.0510 ± 0.0037
Planetary Orbit
Pp (days)	228.776 ± 0.029	303.158 ± 0.072	240.503 ± 0.053
ap (AU)	0.7048 ± 0,0011	0.989 ± 0.016	0.7903 ± 0.0028
ep	0.0069 ± 0.0012	<0.411	0.0359 ± 0.0088
References
	Doyle et al. (2011), Winn et al. (2011)	Orosz et al. (2012a)	Welsh et al. (2015)

Note. Quantities marked with an * are model parameters.

Download table as:  ASCIITypeset image

To study the binaries in the KBHZ sample, we focus on three different scenarios, selected to answer the following questions: (1) If the actual circumbinary planets were Earth-like planets instead of giant planets, then would those planets be habitable?⁵ From necessity, this must be based on current understanding of Earth-like habitability. (2) Assuming that the actual circumbinary planets, which have masses larger than or around that of Saturn, would be able to harbor a Mars-sized exomoon, then would these moons be able to retain a dense and moist atmosphere in order to sustain life? And maybe most importantly, binaries like these are very common, so (3) what are the constraints on habitability of currently unseen planets within the BHZ of similar binaries? Insight into these and other questions is gained by studying well-observed cases, where the environmental conditions experienced by planets under the reign of the binary may be investigated.

It is important to stress that the current investigation is not restricted by dynamical or formation constraints, which may or may not limit the likelihood of these scenarios (Chavez et al. 2013, 2015; Andrade-Ines & Michtchenko 2014; Forgan et al. 2015; Pilat-Lohinger et al. 2014). The basic dynamical constraint imposed by the critical limit surrounding the binary, beyond which planets have stable orbits, is included as it extends into the BHZ in some cases (see Section 3). We assume that exomoons and multiple planets could lie in stable orbits within the BHZ of these binaries, based on our own basic tests using numerical orbital integrations.

Estimates of BHZ boundaries quickly followed the Kepler discoveries (Mason & Clark 2012; Quarles et al. 2012; Haghighipour & Kaltenegger 2013; Kane & Hinkel 2013; Mason et al. 2013). Our approach, as that of others, is to use the models of Kasting et al. (1993), updated by Kopparapu et al. (2013a), derived for single stars to also estimate the more complex BHZ limits. An improved version of the first-order method presented in Mason et al. (2013) is developed and applied here.

HZ limits are determined by the critical fluxes at which planetary surface temperatures are compatible with the existence of liquid water. Using one-dimensional atmospheric models, Kopparapu et al. (2013a) calculated critical fluxes assuming that the planet is exposed to blackbody radiation with an effective temperature ${T}_{{\rm{eff}}}$:

$$\begin{eqnarray}&&{S}_{\mathrm{eff},i}\equiv {S}_{\mathrm{eff},i,\odot }+{a}_{i}{T}_{*}+{b}_{i}{T}_{*}^{2}+{c}_{i}{T}_{*}^{3}+{d}_{i}{T}_{*}^{4}.\end{eqnarray} \tag{ 1 }$$

Here T_* = T_eff − 5780 K, and ${S}_{\mathrm{eff},i}$ and ${a}_{i},{b}_{i},{c}_{i},{d}_{i}$ are the effective critical flux for a solar-mass star and the interpolation coefficients for the ith boundary, respectively (see Table 3 in Kopparapu et al. 2013b).

The limits of the HZ, around either single stars or a binary, are calculated by finding the radius of a circle centered in the barycenter of the system, where the average flux $\langle S\rangle$ received from the star or the binary is equal to the critical flux calculated with Equation (1).

In the case of a single star, $\langle S\rangle$ is equal to the nearly constant flux ${S}_{{\rm{sing}}}=L/{d}^{2}$ measured at distance d. Here L is the bolometric luminosity of the star measured in solar units and d is in astronomical units (AU). More precisely and as we will explain later, L is a corrected luminosity that takes into account the response of the atmosphere to the incident stellar radiation (see Equation (4)).

The radius of the ith HZ edge ${l}_{\mathrm{sing},i}$ will then be given by

$$\begin{eqnarray}&&{S}_{\mathrm{eff},i}=\displaystyle \frac{L}{{l}_{\mathrm{sing},i}^{2}}.\end{eqnarray} \tag{ 2 }$$

Around binaries, planetary atmospheres, even if planets are in a circular orbit, will be subject to a variable flux and spectrum of radiation. The instantaneous bolometric flux on a circumbinary planet is calculated as

$$\begin{eqnarray}&&{S}_{{\rm{bin}}}(d,\theta )=\displaystyle \frac{{L}_{1}}{{R}_{1}^{2}(d,\theta )}+\displaystyle \frac{{L}_{2}}{{R}_{2}^{2}(d,\theta )},\end{eqnarray} \tag{ 3 }$$

where R₁ and R₂ are the distances from the stars to a point in a circle of radius d. θ is the angle between the line joining the stars and a fixed point on the circle. This angle is a function of time, and it will depend on the relative angular velocity of the binary and a test particle on the circle.

R₁ and R₂ are given by

$$\begin{eqnarray*}{R}_{1}^{2}(d,\theta ) & = & {(d-{r}_{2}\mathrm{sin}\theta )}^{2}+{r}_{2}^{2}{\mathrm{cos}}^{2}\theta \\ {R}_{2}^{2}(d,\theta ) & = & {(d+{r}_{1}\mathrm{sin}\theta )}^{2}+{r}_{1}^{2}{\mathrm{cos}}^{2}\theta ,\end{eqnarray*}$$

where r₁ = qr/(1 + q) and r₂ = qr₁ are the instantaneous distances of the stars to the barycenter and q = M₂/M₁ is the binary mass ratio. The relative distance r between the stars is also a function of time for the general case of an eccentric binary orbit.

Since the two stars have different effective temperatures, ${S}_{{\rm{bin}}}$ and ${S}_{\mathrm{eff},i}$ cannot be directly compared as in the case of single stars. However, in the case of stellar twins, the two sources will have the same effective temperature and their fluxes can be summed directly. On the other hand, in disparate binaries the ratio of the flux of the secondary to that of the primary scales with the fourth power of the ratio of their effective temperatures. As a result, the combined spectrum, to first order, may be assumed to have an effective temperature equal to that of the primary and the total flux equals ${S}_{{\rm{bin}}}$ (Mason et al. 2013). A verification of this approximation is discussed below.

A more consistent method to take into account the differences in effective temperatures of the stellar components is to replace bolometric luminosities L₁ and L₂ in Equation (3) by weighted luminosities (Haghighipour & Kaltenegger 2013; Cuntz & Bruntz 2014):

$$\begin{eqnarray}&&L{^{\prime} }_{\mathrm{1,2}}={L}_{\mathrm{1,2}}(1+[{S}_{\mathrm{eff},i}({T}_{\mathrm{1,2}})-{S}_{\mathrm{eff},i,\odot }]{l}_{i,\odot }^{2}).\end{eqnarray} \tag{ 4 }$$

The weights in parentheses account for the different spectral distribution of each source and the response of planetary atmosphere to these spectra.

To obtain the limits of the BHZ, ${S}_{{\rm{bin}}}$ given by Equation (3), the equation should be averaged over a time long enough compared to the periods of the system. If the binary eccentricity is not too large, an analytical expression for $\langle S\rangle$ may be used:

$$\begin{eqnarray}&&\langle {S}_{{\rm{bin}}}\rangle (d)=\displaystyle \sum _{n\ne m}\displaystyle \frac{{L}_{n}^{\prime }}{\sqrt{{A}_{m}^{2}-{B}_{m}^{2}}}\left[1-\displaystyle \frac{1}{\pi }\mathrm{arctan}\left(\displaystyle \frac{{B}_{m}}{\sqrt{{A}_{m}^{2}-{B}_{m}^{2}}}\right)\right],\end{eqnarray} \tag{ 5 }$$

with n, m = 1, 2, ${A}_{m}\equiv {d}^{2}+{r}_{m}^{2}$, and ${B}_{m}\equiv 2\ d\ {r}_{m}$.

Finally, the BHZ limits are calculated by solving the equation

$$\begin{eqnarray}&&{S}_{\mathrm{eff},i}=\langle {S}_{{\rm{bin}}}\rangle ({l}_{\mathrm{bin},i}).\end{eqnarray} \tag{ 6 }$$

Although other authors have developed alternative methods to estimate the limits of the BHZ (Haghighipour & Kaltenegger 2013; Kane & Hinkel 2013; Mason et al. 2013; Cuntz & Bruntz 2014; Müller & Haghighipour 2014), all of them share the most important characteristics of the method presented here. In particular, our method is an improved version of the first-order method used in Mason et al. (2013) by incorporating some key elements of the methods presented by Haghighipour & Kaltenegger (2013) and Cuntz & Bruntz (2014). An interesting advantage of our method with respect to others is the fact that for low eccentricities, the calculation of the BHZ limits requires finding the roots of an analytical formula. No numerical integrations or calculations over a grid of points around the system are required. This feature greatly improves the speed of calculation, which is especially important for the exploration of the vast parameter space of circumbinary habitability.

We have compared the BHZ limits for all the Kepler binaries with known circumbinary planets calculated using our method and those used in three different investigations. The results are presented in Figure 1. We have found that the average of the relative difference between the BHZ limits using the new method and that of Haghighipour & Kaltenegger (2013) is less than 4%.⁶ The difference found by comparing our method with that of Mason et al. (2013) is less than 1%. The limits calculated by Cuntz & Bruntz (2014) are systematically more conservative than all of the other methods, with the inner limits ∼10% farther out and the outer limits ∼10% closer in than the other methods.

Zoom In Zoom Out Reset image size

3.1. Evolution of Insolation and BHZ

Generally, it is possible to assume that moderately separated stars in binaries are formed at the same time. Moreover, stellar evolution shows that stars with different masses evolve at different rates and in different ways. Thus, in order to investigate habitability around a binary, stellar evolution models for each star must be used to determine the time dependence of the BHZ.

It is important to ask whether the fundamental properties of low-mass stars, such as effective temperature and luminosity used to calculate the BHZ, evolve in binaries exactly as they do in single stars. For moderately separated binaries (periods larger than 8 days, stellar separations larger than tens of stellar radii), such as those studied in our sample, the Roche lobes of the stars are one order of magnitude larger than the stellar diameters. For instance, in the case of the tightest binary in our sample (Kepler-47), the filling factor of the primary (ratio of stellar radius to Roche lobe size) is ∼0.2 during PMS (when the star is largest) and ∼0.1 during main sequence (Eggleton 1983). At those filling factors the shape of the star is only slightly modified by the companion and its interior and atmospheric properties are practically the same as for a single star.

HZ evolution is also a natural outcome of the evolution of single stars. In the same way that the continuous HZ is defined for single stars (Kasting et al. 1993), it is generalized to binaries and a continuous circumbinary habitable zone (CBHZ) is defined.

In Figure 2, the evolution of the BHZ for Kepler-47 and its associated CBHZ is shown. We use stellar evolutionary models, at each age, to calculate the instantaneous properties of both stars. Then the limits of the HZ at that time are determined. Shown here are the the Recent Venus and Early Mars criteria. The outer CBHZ limit is defined by the instantaneous outer edge of the BHZ as calculated at the time when the binary flux is minimum. We have tested that this time is close to the time of zero-age main sequence of the primary star, which by definition is the more massive component of the system. Planets in the CBHZ remain inside the BHZ during the entire main-sequence stage of the primary as it is the star with the shortest lifetime. The definition of the CBHZ inner edge is trickier. We define it as the position of the inner edge of the instantaneous BHZ calculated when the primary star abandons the main sequence (close to the end of the nuclear hydrogen-burning phase).

Zoom In Zoom Out Reset image size

In order to assess the present insolation conditions in the KBHZ sample, we plot, in Figure 3, the instantaneous BHZ calculated in the middle of the range of estimated age (see Table 1 for ages and references). Along with BHZ limits, the insolation and photosynthetic photon flux density (PPFD) are shown as a function of planetary orbital phase. See Mason et al. (2015) for details on PPFD calculations.

Zoom In Zoom Out Reset image size

In the BHZ diagrams of Figure 3, the planetary orbit has been placed among different landmark distances: BHZ and CBHZ edges, binary orbit, equivalent Earth insolation distance, and, last but not least, the critical stability distance a_crit, defined as the minimum distance where test particles have stable circumbinary orbits. We compute this critical distance with the semi-empirical relationship of Holman & Wiegert (1999):

$$\begin{eqnarray}{a}_{\mathrm{crit}} & \approx & (1.60+5.10{e}_{{\rm{bin}}}-2.22{e}_{{\rm{bin}}}^{2}\\ & & +4.12\mu -4.27{e}_{{\rm{bin}}}\mu -5.09{\mu }^{2}\\ & & +4.61{e}_{{\rm{bin}}}^{2}{\mu }^{2})\;{a}_{{\rm{bin}}}.\end{eqnarray} \tag{ 7 }$$

Here μ = M₂/(M₁ + M₂), and a_bin and e_bin are the binary semimajor axis and eccentricity, respectively. Although dynamical instability likely occurs outside the critical limit as well, as also discussed by Holman & Wiegert (1999), all the planets in our sample are safely beyond this critical distance. However, in Kepler-16, a_crit lies inside the BHZ, restricting the instantaneous extension of the BHZ as well as the CBHZ; see top left panel of Figure 3. When a_crit is larger than the inner edge of the continuous BHZ, we use its value as the inner edge distance.

In the instantaneous BHZ diagrams of Figure 3, distances at which orbital resonances with the binary occur are also shown. It is very interesting to notice, and likely not coincidental, that in both cases, Kepler-16b and Kepler-453b, the planet remains entirely in between two resonances throughout its orbit. This supports the idea that orbital binary resonances may play a role in the formation and final architecture of circumbinary systems. Notice also that Kepler-47c crosses many tightly spaced resonances, not plotted in its BHZ diagram in Figure 3 for clarity.

The planet Kepler-16b has a nearly circular orbit just inside the maximum greenhouse limit. Kepler-47c has a highly eccentric orbit and leaves the BHZ for part of its orbit. It does, however, currently spend most of its time within the HZ, especially given the slower velocity at apastron (see insolation diagram in the middle right panel of Figure 3). The planet Kepler-453b has a slightly eccentric orbit and always remains in the HZ, just inside the runaway greenhouse limit.

According to our calculations, of all the circumbinary Kepler planets, only Kepler-16b is in the CBHZ, i.e., it will lie inside the BHZ during the whole lifetime of the primary star, which is in turn larger than the age of the universe. This condition would seem at first very favorable for the emergence and evolution of complex life. However, other factors could limit the long-term habitability within environments surrounding Kepler-16 and similar binaries (see Section 8).

Kepler-47c and Kepler-453b do not remain inside the BHZ for the main-sequence lifetime of the primaries. However, both remain inside (although very close to the inner edge of) the BHZ for astrobiologically relevant times (3.5 and 4 Gyr, respectively). If we believe stellar evolution models, Kepler-47c is about to leave or has just left the most optimistic BHZ, being currently in a Venus-like state in terms of stellar insolation.

In order to calculate BHZ evolution and other key quantities for models, the stellar properties available in the PARSEC v1.2 evolutionary tracks are used (Bressan et al. 2012, 2013).⁷ In addition to providing good coverage in metallicity and mass, the PARSEC tracks sample the PMS stages of stellar evolution well (and of course also the main sequence), a key feature when attempting to study the early evolution of stellar rotation. Moreover, these models reflect improvements for low-mass stars (Chen et al. 2014). In the BHMcalc tool we provide access to other evolutionary model results.

So far we have calculated insolation as a way to assess habitability conditions in circumbinary planets. But insolation is not the only factor constraining the capability of a planet to sustain life (for a review of many other factors see Kasting 2009 and references therein). Among other endogenous and probably exogenous factors, the presence of a properly dense atmosphere with the right amount of water is mandatory.

The solar system is a good example of a planetary system where among three planets (and a sizeable moon), located close or inside the HZ, two of them, Venus and Mars (plus the moon), are uninhabitable mainly as a result of factors related to the composition and mass of their atmospheres.

The evolution of a planetary atmosphere is driven by the influx of high-energy radiation and plasma from the host star. This has been a matter of research in the past two decades (Kasting 1996; Lammer et al. 2007, p. 127; Lammer et al. 2009, p. 399; Lammer et al. 2010; Lammer 2013; Tian et al. 2013, p. 567). Although it is not yet clear which factors determine the survival of Earth-like atmospheres or drive their composition to favor habitable conditions, multiple lines of evidence show that the interaction of a planet with its young host star is a strong factor in determining the fate of potentially habitable planet atmospheres (Lammer et al. 2007, p. 127; Lammer et al. 2009, p. 399; Tian 2009; Zendejas et al. 2010; Zuluaga et al. 2013).

Two key factors are involved in this interaction. The first is the level of X-rays and extreme-ultraviolet (EUV) radiation, collectively called XUV radiation, incident on the planet. These photons heat up the upper atmosphere and under certain conditions could lead to massive loss of volatiles (Tian 2009; Tian et al. 2013, p. 567). The flux of charged particles, mainly protons, and their associated magnetic fields (SW) is the other factor. The direct interaction between this continuous flow of particles and fields and planetary atmospheres strongly impacts unmagnetized planets and moons (e.g., Venus, Mars, and the moon). Resulting nonthermal processes may remove large fractions of planetary atmospheres, or in extreme cases they can strip them almost completely (Zendejas et al. 2010). Mars is the best known example of the effect of this kind of aggression eliminating habitability. Recent data from the MAVEN Mission (Rahmati et al. 2014; Jakosky et al. 2015) are showing the increased level of importance that SW interaction with the atmosphere has in shaping the evolution of Mars's capability to have liquid water on its surface (habitability).

Modeling the interaction between the radiation and plasma fields of the star and, in this case, the binary system is key in assessing the habitability of all planets. However, since binaries are relatively novel environments, at least for planetary systems, we need to understand and learn to model aspects of stellar evolution that have been applied for years to single stars. In the following sections we develop a theoretical framework to model the evolution of stellar activity and its role in the changing conditions faced by the atmosphere of a circumbinary habitable planet.

5.1. Evolution of Stellar Rotation

The evolution of stellar rotation has been a matter of intense research over the past several decades, not only in single stars but also in binary stars (for a recent review see Bouvier 2013 and references therein; for the evolution of rotation in binaries see Claret et al. 2005).

The large diversity of stellar rotation periods observed in stars of all ages (especially PMS stars) has driven the development of a diverse array of semi-empirical and phenomenological models describing the evolution of stellar AM (Kawaler 1988; Chaboyer et al. 1995; Matt et al. 2012). These models have been extensively tested against samples of stars (both single and binary) in young stellar clusters and more recently in low-mass field stars (Irwin et al. 2011). The details of these models depend on weakly constrained parameters; however, a general consensus has been reached about the general forces and timescales driving stellar rotation evolution.

Here we apply and extend some of these models for stars in moderately separated binaries (i.e., binaries with similar masses and orbits to those of the KBHZ binaries). For details concerning the physical motivations, history of development, and success of these models at reproducing observations, please refer to Bouvier (2013) and references therein.

Although the applicability of these models to stars in binary systems has not been systematically tested, we assume that some of the underlying mechanisms will operate in the stellar components of these binaries as if they were isolated stars. A detailed test of this hypothesis should be pursued in future research efforts. Still, for the range of orbital separations and stellar masses considered here, the evolution in insolation hypothesis remains as a good one for the same arguments given before to support the applicability of single-star evolutionary models.

In general, the AM of the ith layer of a star evolves while obeying Newton's angular second law:

$$\begin{eqnarray}&&\displaystyle \frac{{{dJ}}_{i}}{{dt}}=\sum {{ \mathcal T }}_{i}.\end{eqnarray} \tag{ 8 }$$

Here J_i = I_iΩ_i is the instantaneous AM of the layer, I_i its moment of inertia, Ω_i its instantaneous angular velocity, and $\sum {{ \mathcal T }}_{i}$ the net torque over the layer.

In the case of low-mass stars (M_⋆ < 1.3 M_⊙), it is customary to assume that the star has two distinct interior layers: a radiative inner core (i = core) and a convective envelope (i = conv). It is also assumed that both layers rotate independently as rigid bodies (i.e., no differential rotation).

The convective envelope is subject to two different exogenic effects: (1) tidal breaking, ${{ \mathcal T }}_{{\rm{tid}}}$, and (2) loss of AM via a magnetized SW ${{ \mathcal T }}_{{\rm{SW}}}$.

To calculate ${{ \mathcal T }}_{{\rm{tid}}}$, we use the formalism developed by Hut (1981). Accordingly, a target star of mass ${M}_{{\rm{targ}}}$ and radius ${R}_{{\rm{targ}}}$, rotating at an instantaneous rate Ω and affected by a secondary star of mass ${M}_{{\rm{field}}}$ in an orbit with semimajor axis a, eccentricity e, and mean angular velocity n, feels an effective tidal torque given by

$$\begin{eqnarray}{{ \mathcal T }}_{{\rm{tid}}} & = & -\displaystyle \frac{3{k}_{2}I}{{\kappa }^{2}{t}_{{\rm{diss}}}}{\left(\displaystyle \frac{{M}_{{\rm{field}}}}{{M}_{{\rm{targ}}}}\right)}^{2}{\left(\displaystyle \frac{{R}_{{\rm{targ}}}}{a}\right)}^{6}\\ & & \times \displaystyle \frac{n}{{(1-{e}^{2})}^{6}}\left[{f}_{2}({e}^{2})-{(1-{e}^{2})}^{3/2}{f}_{5}({e}^{2})\displaystyle \frac{{\rm{\Omega }}}{n}\right].\end{eqnarray} \tag{ 9 }$$

Here k₂ is the apsidal motion constant that quantifies the response of the star to tidal distortion, and it depends in general on the internal mass distribution of the star; ${t}_{{\rm{diss}}}={f}_{{\rm{diss}}}{({M}_{{\rm{targ}}}{R}_{{\rm{targ}}}^{2}/{L}_{{\rm{targ}}})}^{1/3}$ is the timescale for dissipation of tidal energy inside the target star, and $\kappa \equiv \sqrt{I/({M}_{{\rm{targ}}}{R}_{{\rm{targ}}}^{2})}$ is the gyration radius. Eccentricity functions f₂(e²) and f₅(e²) are defined by Equations (2) and (3) in Hut (1981).

We assume that while tides are raised over the whole star, only tides within the convective layer contribute to AM dissipation (Zahn 2008). As a result, the torque in Equation (9) only affects the envelope of the star. Alternatively, we can assume that tides raised in the radiative core are dissipated on timescales much larger than that in the convective zone, rendering the resulting torque comparatively small.

The parameter f_diss quantifies the efficiency of tidal energy dissipation. Although a fundamental theory of this process capable of correctly explaining observations of close binaries is still lacking (see Claret 2011 for a discussion), we assume for simplicity that f_diss ≈ 1, as is used, for example, in Claret (2012). Alternative theories predict values of f_diss within one order of magnitude. Thus, for instance, the turbulent dissipation theory of Zahn (2008) predicts a value f_diss ≈ 3.5.

For nonnegligible eccentricities, tidal braking could drive stars to a pseudosynchronous final state where rotational angular velocity is related to, but not exactly, the average orbital angular velocity (Hut 1981),

$$\begin{eqnarray}&&{n}_{{\rm{sync}}}\equiv \displaystyle \frac{{{\rm{\Omega }}}_{{\rm{sync}}}}{\omega }=\displaystyle \frac{1+\frac{15}{2}{e}^{2}+\frac{45}{8}{e}^{4}+\frac{5}{16}{e}^{6}}{{(1-{e}^{2})}^{\frac{3}{2}}(1+3{e}^{2}+\frac{3}{8}{e}^{4})}.\end{eqnarray} \tag{ 10 }$$

For binary eccentricities e in the range 0–0.5, $1\lt {n}_{{\rm{sync}}}\lt 2.8$, implying that provided enough time, stellar components of a binary will be locked in a rotational state where their periods are close to that of the binary period. This result will not depend on the mass-loss history of each star. Once locked, the magnetic activity of the stars becomes frozen at a given value, which in some cases could mimic the behavior of a star that is much younger or older (Mason et al. 2013).

The second exogenic effect removing AM from the convective layer is mass loss via magnetized SWs. Not only does the mass leaving the surface of the star carry AM, but in magnetically active stars the flow of charged particles is also coupled with the corotating magnetic field, up to few to several tens of stellar radii. This leads to a loss of AM that is several orders of magnitude larger than that in stars without a convective envelope (Schatzman 1962).

Although a successful, comprehensive, and self-consistent model of magnetized SW AM loss has not yet been developed, a semi-empirical formula, developed by Kawaler (1988) and modified by Chaboyer et al. (1995), succeeds at producing results consistent with observations for a large range of stellar masses and ages. Accordingly, the torque exerted by the magnetized SW scales with angular velocity Ω, stellar mass M_⋆, and radius R_⋆, following

$$\begin{eqnarray}&&{{ \mathcal T }}_{{\rm{SW}}}={K}_{{\rm{C}}}\;{\rm{\Omega }}\ {\rm{min}}{({\rm{\Omega }},{{\rm{\Omega }}}_{{\rm{sat}}})}^{2}\;{\left(\displaystyle \frac{{R}_{\star }}{{R}_{\odot }}\right)}^{1/2}\;{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{-1/2}.\end{eqnarray} \tag{ 11 }$$

K_C and Ω_sat are two free parameters that can be adjusted to reproduce the scale of observed rotational rates and their distribution, respectively. In order to reduce the number of free parameters in our comprehensive model, Ω_sat is calculated following the scaling relationship

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{sat}}}={{\rm{\Omega }}}_{\odot }\left(\displaystyle \frac{{\tau }_{{\rm{conv}},\odot }}{{\tau }_{{\rm{conv}}}}\right),\end{eqnarray} \tag{ 12 }$$

where τ_conv is the convective overturn time (see Section 5.2.2). This scaling relationship has proven to be useful in reproducing the distribution of rotation periods for low-mass field stars (Irwin et al. 2011).

More recently, success in modeling the rotational rate distribution of stars in young clusters was achieved (Matt et al. 2012; Gallet & Bouvier 2013) using a semianalytical fit of MHD simulations to the original model of Kawaler (1988) along with a recent model for stellar mass loss (Cranmer & Saar 2011). However, its application is still restricted to stars with mass nearly identical to the Sun. Our interest here is in stars ranging from 0.2 to 1.05 M_⊙; see Table 1. A comparison between the results obtained with the semi-empirical formula used here (Equation (11)) and those obtained with the formula by Matt et al. (2012) is left for future work, but it is not expected to alter the major conclusions of the current work.

As the convective envelope loses AM, a difference in angular velocity between the envelope and the radiative core is established. A variety of physical processes, such as magnetic fields, hydrodynamical instabilities, and gravity waves, may redistribute AM in the stellar interior, reducing this difference with time. To model AM redistribution, we follow, as is customary, the prescription by MacGregor & Brenner (1991), by estimating the mutual torque due to the difference in rotation rates as given by

$$\begin{eqnarray}&&| {{ \mathcal T }}_{{\rm{\Delta }}J}| =\displaystyle \frac{{\rm{\Delta }}J}{{\tau }_{{\rm{ce}}}}.\end{eqnarray} \tag{ 13 }$$

The AM exchanged is ΔJ = (I_convJ_core − I_coreJ_conv)/I_tot, and τ_ce is a free parameter called the core–envelope coupling timescale.

In summary, the AM of the convective envelope and the radiative core of each star in the binary evolve following the equations

$$\begin{eqnarray}&&\displaystyle \frac{{{dJ}}_{{\rm{conv}}}}{{dt}}={{ \mathcal T }}_{{\rm{tid}}}+{{ \mathcal T }}_{{\rm{SW}}}+{{ \mathcal T }}_{{\rm{\Delta }}J}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{dJ}}_{{\rm{rad}}}}{{dt}}=-{{ \mathcal T }}_{{\rm{\Delta }}J}.\end{eqnarray} \tag{ 15 }$$

The moment of inertia of both the convective envelope and the radiative core evolves in time and affects the rotational history of the whole star. PMS evolution is especially important in this respect. During PMS, stars contract the most. A pragmatic approach to account for this early contraction is to integrate the angular velocity instead of the AM using

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Omega }}}_{i}}{{dt}}=\displaystyle \frac{1}{{I}_{i}}\displaystyle \frac{{{dJ}}_{i}}{{dt}}-{{\rm{\Omega }}}_{i}\displaystyle \frac{{{dI}}_{i}}{{dt}},\end{eqnarray} \tag{ 16 }$$

where the term $-{{\rm{\Omega }}}_{i}({{dI}}_{i}/{dt})$ accounts for the contraction (expansion) of the respective layer. It is important to stress that although the larger changes in I_i happen during the PMS, we are also integrating Equation (16) during the main sequence.

From a purely dynamical point of view this contraction term is associated with two discrepancies. The rapid contraction of the star during the first stages of stellar evolution is compounded with the gain of specific AM by accreting gas. Angular acceleration of the star may continue even to near breakup velocities (P_rot ∼ 0.1 days). Observations, however, show that the rotation rates of PMS stars are orders of magnitude slower (with periods of 2–10 days). Although still debated, the solution to this discrepancy may rest on the assumption that stars transfer this AM excess to their accretion disk. This disk locking persists for a timescale τ_disk ∼ 2–5 Myr, ending with disk dissipation.

The second discrepancy arises when considering the value and sign of the contraction term for the radiative core. In stellar evolution simulations, the core suddenly forms at around 10–100 Myr and grows to its final size over a time an order of magnitude smaller. As a result, the moment of inertia of the core increases rapidly, potentially producing a sudden rotational breaking of this layer. Even after considering the exchange of AM during the growth of the core, as is done in Allain (1998), this sudden braking still occurs. Although unobservable, a sudden deceleration of the radiative core would produce a signature in the evolution of the envelope and would create a discrepancy with observations.

To be consistent with previous results, we assume that the contraction term for the core in Equation (16) is the same as that of the convective envelope. An investigation of how this apparent deceleration in the radiative core seems to be prevented is a critical matter in need of further research.

For each system studied here, numerical integrations of Equations (16) were performed for both stars in the binary including their mutual tidal interaction. For that purpose we take as inputs the values of key stellar structure properties as a function of time (e.g., apsidal motion constants, moments of inertia, and convective zone thickness), from stellar evolutionary models (Bressan et al. 2012, 2013) as well as from interior structure models (I. Baraffe 2014, personal communication).

5.2. Magnetic Activity

The ultimate purpose of calculating instantaneous stellar rotational periods is to calculate the magnetic activity as a function of time. Magnetic phenomena in the atmospheres of cool stars are driven by the action of a dynamo in their convective zones that in turn is responsible for producing SWs and high-energy radiation (XUV emission).

Progress in understanding and characterizing these processes in the Sun allowed S. Cranmer and S. Saar to develop a self-consistent model of SW acceleration applicable to other cool stars (Cranmer & Saar 2011). Their model accurately predicts stellar activity observables such as the average photospheric magnetic field strength and mass-loss rates as a function of several basic stellar properties such as mass, radius, luminosity, and rotational period.

Their publicly available routine BOREAS is applied here in order to calculate the instantaneous mass-loss rate of each star in the binary system. For details on the input physics, assumptions, and free parameters used by the models of Cranmer and Saar, refer to Cranmer & Saar (2011).

5.2.1. Stellar Wind

Armed with the mass-loss rate $\dot{M}$, the properties of the SW are calculated in a straightforward manner. We assume that the wind reaches its terminal velocity v_SW at a distance r_t from the stellar center, being equal to the local escape velocity:

$$\begin{eqnarray}&&{v}_{{\rm{SW}}}=\sqrt{\displaystyle \frac{2{{GM}}_{\star }}{{r}_{t}}}.\end{eqnarray} \tag{ 17 }$$

A distance r_t = 2.5 R_⋆ is assumed in order to reproduce the current average solar wind condition as observed at Earth distance.

By mass conservation, the particle mean density as measured at a distance r from the star is given by

$$\begin{eqnarray}&&{n}_{{\rm{SW}}}=\displaystyle \frac{\dot{M}}{\mu {m}_{p}4\pi {r}^{2}{v}_{{\rm{SW}}}},\end{eqnarray} \tag{ 18 }$$

where μ is the mean molecular weight of the wind particles and m_p is the proton mass. For a fully magnetized SW μ = 0.6 (Cranmer 2008; Cranmer & Saar 2011).

This stream of charged particles continuously impacts the atmosphere of planets, at a rate given by

$$\begin{eqnarray}&&{F}_{{\rm{SW}}}={n}_{{\rm{SW}}}{v}_{{\rm{SW}}}.\end{eqnarray} \tag{ 19 }$$

At the distance of the HZ, the wind is supersonic and has a negligible thermal pressure. As a result, it exerts a dynamic pressure over planetary magnetospheres given by

$$\begin{eqnarray}&&{P}_{{\rm{SW}}}=\mu {m}_{p}{n}_{{\rm{SW}}}({v}_{{\rm{SW}}}^{2}+{v}_{{\rm{orb}}}^{2}),\end{eqnarray} \tag{ 20 }$$

where ${v}_{{\rm{orb}}}$ is the orbital velocity of the planet.

The Cranmer & Saar model has been developed in order to describe activity in single stars. The previous assumptions concerning how to convert mass-loss rates into SW properties may work well in describing SWs for moderately spaced binaries, as well as isolated stars, but the applicability of these models to binaries has not been tested. Additionally, other effects could arise, modifying the underlying assumptions on which this model is built. For example, a star could illuminate the atmosphere of the companion, thereby modifying the vertical atmospheric structure. Coronal X-ray heating may play a critical role in some cases. The acceleration due to gravity is also modified by the presence of a close companion. Stellar wind shocks could also arise, modifying the simple model described herein. We assume that the acceleration mechanism working in single stars provides a valid constraint for stars in the moderately separated binaries under consideration.

At distances large enough from the binary or in very disparate binaries (q ≪ 1), it is reasonable to assume that the SW flux and dynamic pressure are simply the sum of the quantities produced by each star, i.e.,

$$\begin{eqnarray}&&{F}_{{\rm{SW,bin}}}\approx {F}_{{\rm{SW,1}}}+{F}_{{\rm{SW,2}}}\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{P}_{{\rm{SW,bin}}}\approx {P}_{{\rm{SW,1}}}+{P}_{{\rm{SW,2}}}.\end{eqnarray} \tag{ 22 }$$

For close orbiting planets, these simple assumptions are expected to break down under realistic conditions in some binaries. Recently, the conditions under which the interaction between the winds of binaries could strongly interact, rendering this approximation imprecise, have been studied using MHD simulations. The results of these simulations suggest that even in an extreme case of solar-mass twins, the SW flow is only strongly modified at close circumbinary distances. Within the BHZ the simple approximation given in Equation (21) is applicable for the purposes of studying the plasma environment of potentially habitable planets.

5.2.2. XUV Emission

One of the most detrimental effects that planets face is early stellar activity. During the early history of planets, the high-energy radiation flux far exceeds that which they receive later on. In the case of Earth, for example, it has been estimated (Guinan & Engle 2009) that 4.5 Gyr ago the X-ray flux was 10–100 times larger than the present Earth level (hereafter PEL) owing to solar coronal activity.

The effect that such flux has over time on planetary atmospheres ranges from simply a modification of the upper atmospheric chemistry (Tian et al. 2008; Segura et al. 2010; Hu et al. 2012, 2013; Hu & Seager 2014) to rapid photoevaporation via induced heating of the exosphere (Lammer et al. 2003, 2009, p. 399; Lammer et al. 2010; Zendejas et al. 2010). Both phenomena could have very detrimental effects on planetary habitability. Estimating the flux of high-energy radiation is thus an essential element in the assessment of habitability of any planetary environment.

The emission of X-ray and EUV radiation (1–91.2 nm) and their dependence on stellar age and rotational period have been extensively studied in the case of single stars. A strong anticorrelation between age and XUV luminosity has been observed (Ribas et al. 2005; Penz & Micela 2008; Penz et al. 2008).

Binary stars follow a different set of rules. First, rotation and age are not necessarily correlated, especially when the binary period is short (of the order of days) or if the binary eccentricity is larger than about 0.1, owing to tidal torques. Second, and as explained before, strong magnetic interaction between the stellar components in the first tens of Myr could drive rotational evolution and hence activity toward values that are uncommon in single stars.

In spite of the most common recipes used to estimate XUV levels as a function of stellar age (see, e.g., Zuluaga et al. 2012), which will not be applicable in the case of binaries for the reasons given above, we use here an approach consistent with the methods used to calculate mass-loss rates, i.e., using the activity proxies provided by the models of Cranmer & Saar (2011).

Using a catalog of ∼1600 stars whose rotation and X-ray brightness have been measured, Vardavas (2005) and more recently Wright et al. (2013) showed that X-ray luminosity is strongly correlated with the Rossby number (the ratio of inertial to the Coriolis force in the convection zone). The Rossby number is calculated as

$$\begin{eqnarray}&&\mathrm{Ro}={P}_{{\rm{rot}}}/{\tau }_{{\rm{conv}}}.\end{eqnarray} \tag{ 23 }$$

Here P_rot is the rotational period (which determines Coriolis forces) and τ_conv (the convective overturn time) is the typical timescale for convective motions (inertial forces).

Wright et al. (2013) found the following empirical relationship:

$$\begin{eqnarray}{R}_{{\rm{X}}}=\displaystyle \frac{{L}_{{\rm{X}}}}{{L}_{{\rm{bol}}}}={R}_{{\rm{X,sat}}}\times \left\{\begin{array}{ll}1.0 & {\rm{if}}\ {\rm{Ro}}\lt {{\rm{Ro}}}_{{\rm{sat}}}\\ {({\rm{Ro}}/{{\rm{Ro}}}_{{\rm{sat}}})}^{-2.7} & {\rm{otherwise}}\end{array}\right.\end{eqnarray} \tag{ 24 }$$

Here Ro_sat and R_X,sat are fitting constants quantifying the Rossby number at which the emission is saturated and the constant ratio of X-ray to bolometric luminosity observed at saturation levels. Using their catalog, Wright et al. (2013) found that (Ro_sat, R_X,sat) are between a high-level range of (0.17, 0.001075) and a low-level range of (0.09, 0.00051).

It is interesting to notice that at a given value of Ro the X-ray luminosity could be predicted using the empirical law in Equation (24) with an uncertainty of a factor of 2 even among a large range of ages, stellar masses, and rotational periods.

In order to incorporate this result into our model, we need to predict the Rossby number of a star at a given stage of its rotational evolution. For that purpose, the routines of Cranmer & Saar (2011) that provide the convective overturn times as a function of stellar effective temperature are used.

In order to calibrate binary star models, we have calculated the evolution of rotation and activity for an isolated solar-mass star. The results are presented in Figure 4. We call this the Solar Reference Model (SRM).

Zoom In Zoom Out Reset image size

The free parameters of the model that were described in previous sections have been chosen to reproduce the observed properties of the Sun, namely, the period of rotation, mass-loss rate, and XUV luminosity at present times.

Since activity and rotation of the Sun during the early phases of solar system evolution are poorly constrained, three scenarios are considered: (1) a nominal scenario (initial rotation period P_rot,ini = 7 days, τ_ce = 7 Myr), (2) an initially fast-rotating Sun (P_rot,ini = 1.5 days, τ_ce = 15 Myr), and (3) an initially slowly rotating Sun (P_rot,ini = 12 days, τ_ce = 1 Myr). In each case we use a disk locking time τ_disk = 2 Myr, in agreement with the estimations of the solar nebula dissipation time (Krot et al. 2005). For each scenario the value of the braking constant K_C has been set so as to reproduce the properties of the Sun at τ = 4.56 Gyr (see legend in upper panel of Figure 4).

Using rotation periods, the mass-loss rate, calculated from the BOREAS routine, as a function of time for each scenario is shown in the lower panel of Figure 4. The XUV luminosity has been calculated in each case using Equation (24) and by assuming a critical Rossby number Ro_sat = 0.16 and a saturation X-ray luminosity of R_X,sat = 7.4 × 10⁻⁴. These values are chosen in order to reproduce the observed X-ray luminosity of the Sun.

With the calculated XUV luminosities and mass-loss rates in hand, the SW and XUV fluxes measured at the distance of Venus, Earth, and Mars are determined. From these fluxes, we calculate the integrated effect on a planetary atmosphere:

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{\rm{XUV,SW}}}(\tau ;{\tau }_{{\rm{ini}}})={\displaystyle \int }_{{\tau }_{{\rm{ini}}}}^{\tau }{F}_{{\rm{XUV,SW}}}(t){dt}.\end{eqnarray} \tag{ 25 }$$

Integrated values of XUV and SW fluxes are shown in Figure 5. Integrated fluxes rapidly reach a maximum around 200–400 Myr in the case of SW and ∼1 Gyr in the case of XUV radiation. At about those times, any hypothetical primitive planetary atmosphere (protoatmosphere) has received most of the incoming flux that it will ever receive.

Zoom In Zoom Out Reset image size

If we assume that a "secondary" atmosphere is degassed or left over from the formation process (see, e.g., Lammer 2013 and Lammer et al. 2013 and references therein), we can also calculate the cumulative effect that SW and XUV radiation has on these secondary envelopes. In the inset panels of Figure 5, we plot the integrated flux starting at 700 Myr (around 3.8 Gyr before present) and up to 8 Gyr. Again, most of the integrated effects reach an asymptotic value after 200–300 Myr (around 1 Gyr after formation).

Values of the integrated particle and photon fluxes can be converted into total atmospheric mass loss from unmagnetized planets. However, converting one quantity into the other is nontrivial given the complexity of the nonthermal and thermal mass-loss processes involved. For present purposes, first-order estimations can be obtained with the phenomenological formalisms of Zendejas et al. (2010) (for atmospheric stripping) and Sanz-Forcada et al. (2011) (for XUV-driven erosion).

Using the asymptotic values of the integrated fluxes in Figure 5, a total solar-driven mass loss of the order of 0.1–0.2 bars for a hypothetical "primary" atmosphere of Mars is obtained. Venus could have lost 1.5–3 bars from any protoatmosphere as a result of the same mechanism. In both cases, and to be conservative, we assumed that both planets started with CO₂-rich atmospheres.⁸

For secondary atmospheres (i.e., any atmosphere remaining or degassed after several hundreds of Myr), we estimate that Mars and Venus were subject to mass losses of the order of 0.06–0.1 bars and 0.8–1.5 bars, respectively. These numbers are in agreement with previous estimations (see, e.g., Kulikov et al. 2006).

In the case of Mars, which currently has a 6 mbar atmosphere, an estimated mass loss of ∼100 mbar is fairly compatible with an early loss of "habitable" conditions (the presence of surface water). Venus, having a current ∼100 bar atmosphere, seems to have barely lost 1% of its mass via direct stripping from solar wind.

If, on the other hand, we assume that terrestrial planets started with massive hydrogen/helium envelopes (see Lammer 2013 and references therein), subject to photoevaporation from the absorption of XUV radiation, we estimate (from the integrated XUV fluxes) that the protoatmosphere of Venus could have lost 0.3%–1% of its total mass in the first 300 Myr. This is of the same order as the mass a planet with its size could have accumulated from the planetary nebula or from other early outgassing processes (Elkins-Tanton & Seager 2008).

In summary, and recognizing that other complex endogenous and exogenous factors are involved in the evolution of the atmosphere of terrestrial planets and moons (vulcanism, degassing, impacts, etc.), integrated SW and XUV fluxes calculated for the reference solar model, as well as for the KBHZ sample, are fairly compatible with the state of affairs in the solar system. In consequence, our model can be used to place first-order constraints on the critical level of stellar aggression that planets around single and binary stars could potentially withstand, before losing their water inventory and becoming uninhabitable.

The model developed and tested in previous sections is applied in order to calculate the evolution of rotation and activity of the KBHZ sample. Figure 6 shows rotation rate as a function of time for each star in the three binary systems.

Zoom In Zoom Out Reset image size

In each case, the equation of motion (Equation (16)) including the effect of tidal interaction is integrated. For primary stars, we use the same parameters of the rotational evolution model (${{\rm{\Omega }}}_{{\rm{sat}}}$, ${K}_{{\rm{C}}}$, ${P}_{{\rm{ini}}}$) as those of the solar reference nominal model. Using the same parameters for the secondaries (all of them are low-mass stars, M < 0.36 M_⊙) produces anomalously low rotational rates at late times. In order to be consistent with observed periods or rotation of field M dwarfs, we use a much lower braking parameter, K_C = 10⁴⁰, for companions. With this adjustment the rotation periods of Proxima Centauri and Barnard's Star are both accurately predicted.

Although there are discrepancies between the rotation periods of the primary star measured photometrically and spectroscopically (blue and cyan error bars in Figure 6, respectively), the model (blue line) predicts the observed values reasonably well considering the large uncertainty box.

The observed period of Kepler-16A, P_rot = 35.1 (35.8) days, is very close to that expected if the star is synchronized with the binary orbit. For the present eccentricity of the system the pseudosynchronous period is P_sync = 35.6 days (see Equation (10)). However, the tidal torque, even for the shortest expected value of tidal dissipation, is too low to achieve synchronization during the first couple of Gyr (notice that the rotational evolution calculated without tidal interaction is shown, using dashed lines, in Figure 6), for reference. Interestingly, the rotation evolution model also predicts a rotation period very close to the pseudosynchronous period for the estimated median age. This may be considered a major accomplishment for the proposed model.

The case of Kepler-47A is also very interesting. The period of rotation measured photometrically, P_rot,photo = 7.7 days, is close to although larger than the pseudosyncronization period of P_sync = 7.4 days. Independently, a spectroscopically determined velocity of rotation predicts a slightly larger rotation period P_rot,vel = 11.9 days (provided that the inclination angle is close to 90°). Stars in this system are close enough to have their rotation periods substantially affected by tides. Our rotation evolution model predicts that within the estimated range of ages for the system (1–5 Gyr), the rotation period will reach a maximum of ∼10 days. This value is in between the independently measured periods. The primary star in the model is not synchronized as early as first-order models suggest. Instead, once all the effects modifying the rotational rate during the PMS are taken into account, tides become relevant only after several hundreds of Myr.

The case of Kepler-453 is more conventional, with tides having a negligible effect on rotation during the relevant stellar evolution stages. The nominal model predicts a period of rotation slightly larger than the photometrically measured one, but nearly identical to the spectroscopically measured period. These three cases provide evidence that rotational evolution is complex yet predictable by the current model.

The focus of this paper is on the question of how the evolution of rotation and activity translates into a variable plasma and radiation environment within the BHZ of three well-known circumbinary planets. Both effects impact the evolution of potential planets and moons and their habitability.

Figures 7–9 show the XUV radiation and SW fluxes as calculated within the BHZ of the binaries (gray shaded region). Instantaneous and integrated fluxes are shown. In order to evaluate the level of stellar aggression experienced around each Kepler binary, we compare these results with those calculated for the SRM (solid thin lines). We have compared fluxes around binaries with those for Venus, Earth, and Mars in the solar system. For each planet we calculate the range of flux experienced by those planets for two extreme solar activity evolution models (fast and slow initial rotation).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

8.1. Kepler-16

8.2. Kepler-47

8.3. Kepler-453

We have developed a comprehensive and largely consistent model to calculate the evolution of stellar activity in binaries, aimed at constraining the radiation and plasma environment of BHZ planets. We apply the model to the Kepler binaries with known planets in the BHZ, the KBHZ sample, namely, Kepler-16b, Kepler-47c, and Kepler-453b.

We find that Kepler-16 probably has a hospitable plasma environment for the retention of atmospheres of Mars-sized planets and exomoons in the BHZ. Our model predicts that the integrated SW flux anywhere in the BHZ of the system is lower than that measured at the distance of Mars in the solar system. The levels of high-energy radiation are, however, similar to those observed in the solar system.

Kepler-47 is the only system among the three in the KBHZ sample where tidal interaction has shaped the recent history of stellar activity. High-energy (XUV) radiation levels are affected the most by the evolution of rotation. The historical levels of XUV radiation predicted by our model for this system are lower than those measured in the solar system. This result hints at an intriguing possibility that Earth-like atmospheres could survive desiccation at distances analogous to that of Venus in the solar system. The predicted plasma flux in this system is, however, too high for Mars-sized planets to retain their primary atmospheres. If, however, other mechanisms are responsible for the degassing of a secondary atmosphere, the late integrated fluxes seem to be lower than those predicted in the solar system.

Finally, the case of Kepler-453 seems to be rather conventional. The circumbinary plasma and radiation environment is quite similar to that found in the solar system. However, at late times, the integrated high-energy radiation as measured at the inner edge of the HZ (where the actual planet Kepler-453b lies) could be better than that experienced by Venus.

These results are not definitive, but they provide a solid starting point to consistently assess the key problem of the radiation and plasma environment of planets in binaries. Rotational evolution models of stars in binaries need to be further improved and, more importantly, compared with a larger sample of observed rotations in binaries. Models of the relationship between stellar activity and fundamental properties and rotation period of stars, across a wide range of stellar masses, need to be further developed and tested. We cannot rely only on age–activity relationships found for single stars in order to assess the evolution of activity of stars in binaries. Stellar rotation is modified by mutual tidal interaction. Any future improvement of these aspects of the models can be tested using the computational tool we have developed and made avaiable at http://bhmcalc.net.

We thank W. Welsh for discussions concerning the Kepler planets and N. Haghighapor for pointing out the importance of resonances. We also thank S. Cranmer, G. Torres, and A. Claret for discussion about their models applied in this work. We appreciate the help of I. Baraffe, who provided detailed stellar evolutionary models, including moment-of-inertia calculations, required in the rotational evolution calculations. Special thanks to the anonymous referee, who provided useful input for improving the BHZ models and suggestions for improving the clarity of the manuscript. J.I.Z. and P.A.C. are supported by Estrategia de Sostenibilidad 2014-2015 de la Universidad de Antioquia and by COID/UdeA. P.A.M. is supported by NMSU-DACC. J.I.Z. is also supported by the Fulbright Commission, Colombia. He also thanks the Harvard-Smithsonian Center for Astrophysics for their hospitality.