ROTATIONAL SYNCHRONIZATION MAY ENHANCE HABITABILITY FOR CIRCUMBINARY PLANETS: KEPLER BINARY CASE STUDIES

Orbital stability is a prerequisite for habitability. To be stable, the semi-major axis of a circumbinary planet must be larger than a critical value a_c. For a large range of binary semi-major axes a_bin and eccentricities e, a_c is calculated from numerical fits, such as that provided by Equation (3) of Holman & Wiegert (1999).

Orbital stability within the RHZ is of greatest concern for similar mass binaries with large separations. For twin binaries in circular orbits, i.e., μ = 1/2 and e = 0, the stability criterion simplifies to a_c ∼ 2.4 a_bin. In this case, if a_bin > l_out/2.4, where l_out is the outer edge of the RHZ, planets throughout the habitable zone have unstable orbits, rendering the binary uninhabitable.

RHZ limits are found by calculating fluxes which allow the existence of liquid water. We follow the results of Kopparapu et al. (2013) that, for single stars, provide limiting fluxes S_eff (in units of present Earth solar flux) as a function of stellar effective temperature T_* = T_eff − 5780 K:

$$\begin{equation} S_{\rm eff}=S_{{\rm eff}\odot }+aT_{*}+bT^{2}_{*}+cT^{3}_{*}+dT^{4}_{*}. \end{equation} \tag{ 1 }$$

Here S_eff☉, a, b, c and d are constants which depend on the physical criterion defining a given limit and are tabulated in Table 3 of Kopparapu et al. (2013). The most generous limits are obtained by assuming that Venus had surface liquid water until a few Gyr ago (the recent Venus criterion) and that Mars was also habitable at the beginning of solar system evolution (the early Mars criterion).

Limits of the RHZ, either around single stars or binaries, are defined as the distance d where the averaged flux 〈S(d)〉 is equal to the critical flux:

$$\begin{equation} \langle S(d)\rangle = S_{\rm eff}. \end{equation} \tag{ 2 }$$

For single stars, circular orbits, and fast rotating planets, 〈S(d)〉 = L_*/d², where L_* is the stellar luminosity in solar units and d is in AU.

In order to apply Equation (1) to binaries, we first verify that an effective temperature can be associated with the combined stellar flux. For twins, the result is trivial, but for disparate binaries, a more complex situation arises. We have verified numerically that the peak flux in disparate systems remains close to that of the primary. So conservatively, we assume that the effective temperature is that of the most luminous star, since it has the dominant effect on the RHZ.

The flux at a distance d from the center of mass (CM), varies with the binary phase angle θ as:

$$\begin{equation} S_{\rm bin}(d,\theta)=\frac{L_1}{R_1^2(d,\theta)}+\frac{L_2}{R_2^2(d,\theta)} \end{equation} \tag{ 3 }$$

where R₁ and R₂ are planetary distances to the primary and secondary components. Assuming that the planetary orbit is in the same plane as the binary,

$$\begin{eqnarray*} R_1^2(d,\theta) & = & (d+r_1\sin \theta)^2+r_1^2\cos ^2\theta \\ R_2^2(d,\theta) & = & (d-r_2\sin \theta)^2+r_2^2\cos ^2\theta. \end{eqnarray*}$$

Here r₂ = a_bin/(1 + q) and r₁ = qr₂ are the average star–CM distances and q = M₂/M₁ is the binary mass ratio. To calculate RHZ limits, we compare the combined flux (Equation (3)) averaged over the binary orbit, with the effective flux for a given RHZ edge:

$$\begin{equation} \frac{1}{2\pi }\int _0^{2\pi }\left(\frac{L_1}{R_1^2(d,\theta)}+ \frac{L_2}{R_2^2(d,\theta)}\right)d\theta =S_{{\rm eff}}. \end{equation} \tag{ 4 }$$

The result of solving Equation (4) in the case of twins and disparate binaries is presented in Figure 1.

Zoom In Zoom Out Reset image size

Orbital stability and proper insolation do not guarantee the most important condition for habitability: the presence of a dense and wet atmosphere. Although important intrinsic factors are involved in the formation and preservation of planetary atmospheres, planet–star interactions play key roles in the fate of gaseous planetary envelopes, especially during the early, active, phases of stellar evolution. In the next section, we model several aspects of planet–star interactions and apply them to the survival of circumbinary planet atmospheres.

A rigorous treatment of binary stellar winds is challenging (Siscoe & Heinemann 1974). The problem has been extensively studied in the case of early type binaries (Stevens et al. 1992), but less attention has been paid to the case of low-mass stars in binaries. Here, we assume a simplified non-interacting model for the combined stellar wind.

Since for binary periods of P_bin ∼ 10–20 days, orbital velocities (v ∼ 80–100 km s⁻¹) are much lower than wind velocities as measured near the stellar surface (v_sw > 3000 km s⁻¹), we neglect orbital motion effects in calculating stellar wind properties.

For twins, we assume a wind source with a coronal temperature and wind velocity profile equal to that of a single star. Mass-loss rates are assumed to be double of that calculated for single stars of the same type. For disparate binaries, we simply sum the stellar wind pressure from each component.

Stellar wind properties are calculated using Parker's model (Parker 1958). It has been shown (Zuluaga et al. 2013) that planetary magnetospheric properties computed with this model are only ∼10% different than those obtained with more realistic models.

The stellar wind's average particle velocity, v_sw(d), at a distance d from the host star is obtained by solving Parker's equation (Equation (35) in Zuluaga et al. 2013). The density profile, n_sw(d) is obtained by equating particle velocity and mass-loss rate $\dot{M_\star }$:

$$\begin{equation} n_{\rm sw}=\frac{\dot{M_\star }}{4\pi d^2 v_{\rm sw} m}. \end{equation} \tag{ 5 }$$

A procedure to estimate Parker model parameters as a function of stellar age for single stars was devised by Grießmeier et al. (2007). It relies on empirical relationships between age, X-ray flux, and rotation of single stars (Wood et al. 2005). According to these relations, the product of mass-loss rate and wind velocity $\dot{M}v_{\rm sw}$ scale with stellar rotation period P_rot following:

$$\begin{equation} \dot{M}v_{\rm sw}\propto P_{\rm rot}^{-3.3}. \end{equation} \tag{ 6 }$$

We adapt these results to binaries by introducing the so-called rotational age, τ_rot defined as the equivalent stellar age at which rotational period of a single star P_{rot, single} will be equal to rotational period of a star in the binary P_{rot, bin}:

$$\begin{equation} P_{\rm rot,single}(\tau _{\rm rot})=P_{\rm rot,bin}. \end{equation} \tag{ 7 }$$

For non-negligible binary eccentricity, stars eventually reach a pseudo-synchronous state with P_{rot, bin} = P_bin/n_sync (see Section 3.2) where n_sync is a real number depending on eccentricity (Hut 1981).

Applying these isolated star relationships to binaries is reasonable since physical mechanisms connecting rotation, age, and activity would not be different for stars in binaries with separations of tens to hundreds of stellar radii, i.e., for P_bin > 5 days.

Rotational ages for stars in binaries are depicted in Figure 2. We see that F, G, and late K stars (M_⋆ > 0.8) in binaries with orbital periods P_{bin, rot} < 20 days, could experience premature aging if they are quickly tidally locked. They would appear as old as single stars with τ > 3 Gyr, in terms of magnetic and stellar activity. Premature aging might be an advantage for circumbinary terrestrial planet atmospheres. On the other hand, lower mass stars in binaries with similar periods exhibit a forever-young effect, i.e., components freeze at rotational ages less than approximately 2 Gyr, thereby reducing habitability.

Zoom In Zoom Out Reset image size

Stellar XUV luminosity depends on chromospheric and coronal activity, which in turn depends on rotation. Since rotation of single main sequence stars slows down with age, XUV luminosity should also decrease with time (Garces et al. 2011). The rotational aging mechanism will reduce XUV luminosities and their potential harmful effects on planetary atmospheres.

Despite large uncertainties in measured XUV stellar emission (Pizzolato et al. 2003), several authors have developed simple empirical laws expressing XUV luminosities, or its proxy, X-ray luminosity, as a function of age. To be conservative, we use the law by Garces et al. 2011 providing X-ray luminosity of GKM-types as a power-law of stellar age:

$$\begin{equation} L_X=\left\lbrace \begin{array}{@{}ll}6.3\times 10^{-4} L_\star & {\rm if}\;\tau <\tau _{\rm i} \\ 1.89\times 10^{28}\;\tau ^{-1.55} & {\rm otherwise} \end{array} \right. \end{equation} \tag{ 8 }$$

where L_⋆ is the bolometric luminosity and τ_i is the so-called saturation time scaling with L_⋆ according to:

$$\begin{equation} \tau _i=0.06\,{\rm Gyr}\;\left(\frac{L_\star }{L_\odot }\right)^{-0.65}. \end{equation} \tag{ 9 }$$

For high X-ray luminosities, we approximate L_XUV ≈ L_X (Guinan & Eagle 2009). Using this model, we verify that XUV Present Earth Level (PEL) is 0.88 erg cm⁻² s⁻¹, in agreement with the observed value (Judge et al. 2003). We also predict that at τ_☉ ≈ 1 Gyr the Earth XUV flux was F_XUV = 8 PEL in agreement with previous estimates (see, e.g., Kulikov et al. 2006).

Benefits are maximized if tidal locking occurs, at least for the primary component, before the rise of the secondary planetary atmosphere. For solar system terrestrial planets, this time is estimated as τ_atm ∼ 0.3–0.8 Gyr (Hart 1978; Hunten 1993). Therefore, in order to evaluate if circumbinary planets benefit from rotational aging we estimate synchronization times.

For a target star with initial rotational and orbital angular velocities Ω = 2π/P_rot and ω = 2π/P_bin, subject to secondary star tides placed in an orbit with eccentricity e, synchronization time t_sync is calculated following Zahn (2008), see Figure 3,

$$\begin{eqnarray} \frac{1}{t_{{\rm sync}}} &=& \frac{1}{t_{{\rm diss}}} \frac{f_2(e^2)}{(1-e^2)^6}\left[1-\frac{(1-e^2)^{3/2}f_5(e^2)}{f_2(e^2)}\frac{\omega }{\Omega }\right] \nonumber\\ && \times \left(\frac{M_{\rm field}}{M_{\rm targ}}\right)^{2}\frac{M_{\rm targ}R_{\rm targ}^{2}}{I}\left(\frac{R}{a}^{6}\right) \end{eqnarray} \tag{ 10 }$$

where:

$$\begin{eqnarray*} f_2(e^2) & = & 1+\frac{15}{2}e^2+\frac{45}{8}e^4+\frac{5}{16}e^6\\ f_5(e^2) & = & 1+3e^2+\frac{3}{8}e^4. \end{eqnarray*}$$

Zoom In Zoom Out Reset image size

Here I is the moment of inertia of the target star, and M_targ and R_targ are its mass and radius. M_field is the mass of the star producing the tidal field. Moments of inertia MoI ≡ I/MR² have been calculated from stellar evolution models (Claret & Gimenez 1990). Zero-age main sequence values of MoI = 0.08 for solar-mass stars, MoI = 0.1 for a 0.8 M_☉ star and MoI = 0.14 for 0.6 M_☉ stars have been used to interpolate the value of this parameter for other masses. For less massive stars, we use values close to 0.23, which is the limit for less centrally concentrated substellar objects (Leconte et al. 2011).

Since we are interested in low-mass binaries, i.e., M_⋆ < 1.5M_☉, for which convection happens in the whole star or in the outer envelope (Chabrier & Baraffe 1997), we assume that turbulent convection is the dominant tidal dissipation mechanism. We compute the viscous dissipation time t_diss directly from convection overturn times following Equation (2.8) in Zahn (2008).

For non-negligible eccentricities, tidal breaking drives stars to a pseudo-synchronous final state where rotational angular velocity is not exactly the average orbital angular velocity. Final rotational velocity in eccentric binaries is obtained when the term in the square brackets in Equation (10) becomes zero:

$$\begin{equation} n_{\rm sync}\equiv \frac{\Omega _{\rm sync}}{\omega }= \frac{1+\frac{15}{2}e^{2}+\frac{45}{8}e^{4}+\frac{5}{16}e^{6}}{(1-e^{2})^{\frac{3}{2}} \left(1+3e^{2}+\frac{3}{8}e^{4} \right)}. \end{equation} \tag{ 11 }$$

For eccentricities in the range 0–0.5, 1 < n_sync < 2.8. Note that this equation is also Equation (42) of Hut (1981), where we use average angular velocity rather than instantaneous periastron velocity.

For magnetic protection, we use the models of Zuluaga et al. (2013). Magnetosphere sizes, quantified by stand-off distance R_S, are scaled with stellar wind dynamical pressure $P_{\rm sw}=m n_{\rm sw} v_{\rm sw}^2$ and the planetary dipole moment ${\cal M}$ according to:

$$\begin{equation} \frac{R_S}{R_\oplus } = 9.75 \left(\frac{\cal M}{{\cal M}_\oplus }\right)^{1/3} \left(\frac{P_{\rm sw}}{P_{\rm sw\odot }}\right)^{-1/6} \end{equation} \tag{ 12 }$$

where ${\cal M}_\oplus = 7.768\times 10^{22}$ A m² and P_sw☉ = 2.24 × 10⁻⁹ Pa are the present Earth dipole moment and the average dynamical pressure of the solar wind at Earth. We note that during the first two Gyr, our thermal evolution models predict a dipole moment for Earth analogs of around 0.6 ${\cal M}_\oplus$ (see Figure 4 in Zuluaga et al. 2013).

Strong magnetic fields were probably required to create magnetosphere cavities large enough to protect early bloated atmospheres. For a magnetosphere comparable in size, or smaller than, that predicted for early Venus (R_{S, Venus} = 3.8 R_p), subject to similar levels of XUV radiation, we assume magnetic protection is insufficient to prevent water loss similar to Venus. On the other hand, if magnetospheres are larger than that of early Earth R_{S, Earth} = 4.5 R_p, the planet is magnetically protected.

Earth and Venus limits are drawn in the contour plots for R_S in Figures 4 and 5. We stress that while R_S depends on the dipole moment, an intrinsic planetary property, it also depends on the stellar wind pressure and velocity. This is the reason why we can display R_S in the e–P_bin plane for Earth-like magnetic properties.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size