Enlarging habitable zones around binary stars in hostile environments

TL;DR

This study shows that in dense star-forming regions, interactions can enlarge habitable zones around binary stars, suggesting such environments are not always hostile to planet habitability.

Contribution

The paper demonstrates through N-body simulations that binary star systems in dense regions can have enlarged habitable zones due to stellar interactions, challenging previous assumptions.

Findings

Binary systems can have their habitable zones enlarged through stellar interactions.

Enlarged habitable zones occur in 1-2 binaries per simulated region.

Dense star-forming regions are not necessarily hostile to planet habitability.

Abstract

Habitable zones are regions around stars where large bodies of liquid water can be sustained on a planet or satellite. As many stars form in binary systems with non-zero eccentricity, the habitable zones around the component stars of the binary can overlap and be enlarged when the two stars are at periastron (and less often when the stars are at apastron). We perform N-body simulations of the evolution of dense star-forming regions and show that binary systems where the component stars originally have distinct habitable zones can undergo interactions that push the stars closer together, causing the habitable zones to merge and become enlarged. Occasionally, overlapping habitable zones can occur if the component stars move further apart, but the binary becomes more eccentric. Enlargement of habitable zones happens to 1-2 binaries from an average initial total of 352 in each simulated star-forming region, and demonstrates that dense star-forming regions are not always hostile environments for planet formation and evolution.

Tables

Table 1. Table 1: Summary of the adopted binary fraction and semimajor axis distribution as a function of primary mass for stars in the Galactic field which we use to set up the binary populations in our star-forming regions.

Mass range	$f_{\mathrm{bin}}$	$\tilde{a}$	log${}_{10}{}^{}\sigma _a^{}$	refs.
	0.15	4 au	0.4	(1)
	0.34	16 au	0.8	(2),(3),(4)
	0.46	50 au	1.68	(5)
	0.48	230 au	0.79	(6)

Equations

$$p(m)\propto {\left(\frac{m}{\mu }\right)}^{-\alpha }{\left(1+{\left(\frac{m}{\mu }\right)}^{1-\alpha }\right)}^{-\beta },$$

$$f_{\mathrm{bin}}=\frac{B}{S+B}$$

$$l_{\mathrm{x}-\mathrm{Star}}=l_{\mathrm{x}-\odot }{\left[\frac{L_{\mathrm{Star}}/L_{\odot }}{1+{\alpha }_{\mathrm{x}}(T_i)l_{\mathrm{x}-\odot }^2}\right]}^{1/2},$$

$${\alpha }_{\mathrm{x}}(T_i)=a_{\mathrm{x}}{T}_i+b_{\mathrm{x}}{T}_i^2+c_{\mathrm{x}}{T}_i^3+d_{\mathrm{x}}{T}_i^4,$$

$$T_i(\mathrm{K})=T_{\mathrm{Star}}(\mathrm{K})-5870,$$

$$W_i(f,T_i)={\left[1+{\alpha }_{\mathrm{x}}(T_i)l_{\mathrm{x}-\odot }^2\right]}^{-1},$$

$$l_{\mathrm{x}-\mathrm{bin}}={\left[\frac{W_p{L}_p}{L_{\odot }/l_{\mathrm{x}-\odot }^2-W_s{L}_s/r_{\mathrm{pl}-\mathrm{s}}^2}\right]}^{1/2}.$$

Full text

††pagerange: Enlarging habitable zones around binary stars in hostile environments–Enlarging habitable zones around binary stars in hostile environments††pubyear: 2018

1 Introduction

One of the main drivers of exoplanet research is to find an Earth-like planet that orbits within the habitable zone of the planet’s host star. The exact definition and boundaries of habitable zones vary, but it is generally accepted that the habitable zone is the location where water can exist in liquid form (Kasting et al., 1993; Ramirez, 2018), as on Earth.

The calculation of the habitable zone for a single star like the Sun is relatively straightforward, with the flux (i.e. effective temperature) of the host star the main variable. To a lesser extent, the composition of the atmosphere, and the amount of cloud cover on a planet are also determining factors (Kasting et al., 1993). Habitable zones around single stars therefore only change if the luminosity of the host star changes (e.g. due to stellar evolution from the pre- to Main Sequence, and then once the star moves off the Main Sequence toward the end of its life), or the atmosphere of any planet changes drastically.

However, most stars are not born as singles, but rather in binary or higher-order multiple systems where the semimajor axis can be anywhere between $\sim 0.01-10^{4}$ au (Duchêne & Kraus, 2013). The binary fraction of solar-mass primary stars in the Sun’s local Galactic environment is just less than 50 per cent (Raghavan et al., 2010), though lower-mass primaries (e.g. $m_{p}<0.5\,$M⊙) likely have a lower binary fraction (30 – 40 per cent, Bergfors et al., 2010; Janson et al., 2012; Ward-Duong et al., 2015) whereas higher mass primary stars (>1.5 M⊙) have a higher binary fraction (increasing towards 100 per cent for the most massive stars, De Rosa et al., 2014; Sana et al., 2013).

The binary fraction is more difficult to measure in Giant Molecular Clouds where stars (and planets) are forming, and the values in the Galactic disc are thought to be a lower-limit to the initial fraction of stars that are found in binaries. In star-forming regions the stellar density can exceed the value in the Galactic disc by several orders of magnitude (King et al., 2012), and at these densities ($>$100 M⊙ pc*-3*) dynamical interactions can change the initial orbits of binary systems, or destroy the system altogether to produce two single stars. Conversely, it is difficult to form all but the most extreme wide binary ($a_{\rm bin}>1000$ au) systems through dynamical interactions (Kouwenhoven et al., 2010; Moeckel & Bate, 2010), and systems with semimajor axes lower than these values are usually primordial.

Observational and theoretical work has shown that binary stars do not preclude planet formation. Radial velocity surveys (e.g. Raghavan et al., 2006; Bonavita & Desidera, 2007) have found planets orbiting the primary component of binary systems (so-called ‘S-type’, or satellite orbits, Dvorak, 1986) and more recently the Kepler mission has shown that planets can orbit both components of a binary system (Doyle et al., 2011) in a so-called ‘P-type’, or primary orbit (Dvorak, 1986).

Furthermore, dynamical stability calculations have shown that planets can exist on long-lived orbits in S-type systems, usually if they orbit their host star at a distance of $<0.3a_{\rm bin}$ and the eccentricity of the binary is not too high ($e_{\rm bin}<0.9$, Holman & Wiegert, 1999). In this situation, each component of the binary could feasibly host planets, and each component of the binary will have its own habitable zone, unless the component stars are so close that their habitable zones merge. This varies depending on the masses of (and hence fluxes from) the component stars of the binary, but for a mildly eccentric ($e=0.5$) binary with semimajor axis $\sim 5$ au it is possible that the binary habitable zone will merge and be larger than if the two stars were single or widely separated (Kaltenegger & Haghighipour, 2013; Cuntz, 2014; Jaime et al., 2014).

This then raises the interesting possibility that the orbit of a binary star system can be altered due to dynamical interactions in its birth star-forming region to the extent that the previously separate habitable zones merge, and the enhanced incident flux from both stars enlarges the width of the habitable zone around both stars. In star-forming regions, binaries with a binding energy higher than the local average energy of passing stars will not break apart, and are also likely to have their semimajor axis decreased or “hardened". (A binary with binding energy less than the surrounding stars is “soft" and will likely break apart or increase its semimajor axis after an interaction – Heggie, 1975; Hills, 1975.) The boundary between the hard and soft regimes depends on the stellar density, interaction rate and local velocity dispersion in a star-forming region, but the latter two variables are more or less constant.

In very dense star-forming regions ($>$1000M⊙ pc*-3*) the hard-soft boundary is several au, which is the regime in which the habitable zones of main sequence stars in binaries can overlap. In this Letter, we determine how often binary stars are hardened in dense star-forming regions such as the Orion Nebula Cluster (ONC), and whether this causes the habitable zones around the component stars to merge and become enlarged. We describe our simulations of the evolution of star-forming regions and outline our method to calculate the habitable zones in Section 2. We present our results in Section 3 and we conclude in Section 4.

2 Method

2.1 $N$-body simulations

We set up $N$-body simulations of star-forming regions with spatial and kinematic substructure in an attempt to mimic observations that indicate that young star-forming regions appear spatially clumpy and filamentary. The local velocity dispersions are small, meaning that stellar velocities are correlated on local scales. We create substructured star-forming regions containing 1500 stars arranged in a fractal distribution according to the prescription in Goodwin & Whitworth (2004), with fractal dimension $D=1.6$, resulting in a high degree of spatial and kinematic substructure.

The star-forming regions have initial radii of 1 pc, and primary masses drawn from a Maschberger (2013) Initial Mass Function (IMF), which has a probability distribution of the form

[TABLE]

where $\mu=0.2$ M⊙ is the average stellar mass, $\alpha=2.3$ is the Salpeter (1955) power-law exponent for higher mass stars, and $\beta=1.4$ describes the slope of the slope of the IMF for low-mass objects. We sample this distribution in the mass range 0.01 – 50 M⊙. The average local stellar density in these simulations is $\simeq 10^{4}$M⊙ pc*-3*, consistent with inferred initial densities of some star-forming regions such as the ONC (Marks & Kroupa, 2012; Parker, 2014).

We determine whether a star is the primary component of a binary system depending on its mass. Observations of binary stars where the primary component is roughly Solar mass in the local Galactic neighbourhood suggest a binary fraction of $f_{\rm bin}=0.46$, where

[TABLE]

and $S$ and $B$ are the number of single and binary systems, respectively. This fraction decreases for lower-mass primary systems and increases for higher-mass systems. The distribution of orbital semimajor axes of Solar-type primaries can be fit with a log-normal distribution with mean 50 au and a variance log${}_{\rm 10}\sigma_{a}=1.68$. The distribution of mass ratios $q=m_{s}/m_{p}$ is observed to be flat in the field (Reggiani & Meyer, 2011) and the distribution of eccentricities is also flat for binaries with a semimajor axis greater than 0.1 au (Raghavan et al., 2010), with systems on shorter orbits tending to zero eccentricity. We draw mass ratios and eccentricities from these distributions and then draw semimajor axes according to the primary mass (we summarise the adopted binary fraction, peak of the semimajor axis distribution and variance of this distribution for a given primary mass range in Table 1).

We then randomly distribute our primordial binary systems and remaining single stars within the fractal distributions. We use the kira integrator (Portegies Zwart et al., 2001), a 4th-order Hermite scheme $N$-body integrator within the Starlab environment to evolve the simulations for 10 Myr which is long enough for most dynamical processing to act on the binary systems in dense stellar environments. We do not include internal stellar or binary evolution as we do not expect theses processes to be important for stars with $m<3$M⊙ in the 10 Myr timeframe of the simulations.

2.2 Calculating the habitable zones in binaries

In our analysis we take the primary and secondary masses of the binary system, as well as the semimajor axis and eccentricity, to calculate the extent of the habitable zones before dynamical evolution ($t$ = 0 Myr) and at the end-point of the simulation ($t$ = 10 Myr). We use the mass–luminosity–temperature relations for Main Sequence stars in Cox (2000), where $L/L_{\odot}=\left(M/M_{\odot}\right)^{3.5}$ to derive an effective temperature for each star in a given binary system. We do not use pre-main sequence temperatures and luminosities because we are interested in the habitability of the star once it reaches the main sequence, and how different this habitability would have been had the binary system in question not had its orbit altered by a dynamical interaction in the birth star-forming region.

We follow the method adopted by Kaltenegger & Haghighipour (2013) to calculate the narrow and empirical habitable zones around each component star in the binary. A more generalised version of these calculations was derived by Cuntz (2014, 2015) – see also Jaime et al. (2014). For individual stars, the position of the outer or inner boundary of the habitable zone is given by

[TABLE]

where $l_{\rm x-\odot}$ is the inner or outer boundary of the Sun’s habitable zone, $L_{\rm Star}$ is the luminosity of the star, $L_{\odot}$ the luminosity of the Sun. $\alpha_{\rm x}(T_{i})$ is given by

[TABLE]

where

[TABLE]

and $a_{\rm x}$, $b_{\rm x}$, $c_{\rm x}$ and $d_{\rm x}$ are coefficients which depend on the type of habitable zone under consideration and are found in Kopparapu et al. (2013) and also summarised in table 1 in Kaltenegger & Haghighipour (2013). The subscript i refers to either primary component p or secondary component s.

We consider a narrow habitable zone, which varies depending on the Greenhouse Effect from CO2 with boundary values for the Sun of $l_{\rm x-\odot}=0.97$ au for a ‘Runaway Greenhouse’ and $l_{\rm x-\odot}=1.67$ au for a ‘Maximum Greenhouse’ (Kopparapu et al., 2013; Kaltenegger & Haghighipour, 2013). We also consider an empirical habitable zone, where $l_{\rm x-\odot}=0.75$ au for a ‘Recent Venus’ and $l_{\rm x-\odot}=1.77$ au for an ‘Early Mars’; both planets show no evidence for liquid water on their surfaces at these times, indicating they would not be within the habitable zone (Kopparapu et al., 2013).

We then calculate the spectral weight factors, $W_{i}(f,T_{i})$, which determine the relative contribution of the flux from each component of the binary star in determining the habitable zone:

[TABLE]

which depend on the assumed cloud cover on the planet in the habitable zone $f$ and the effective temperature of the star $T_{i}$.

We combine these equations to determine the inner and outer boundaries of the narrow habitable zone, and the empirical habitable zone in the binary system, $l_{\rm x-bin}$ using the following equation

[TABLE]

Here, $W_{p}$ and $W_{s}$ are the spectral weight factors for the primary and secondary, $L_{p}$ and $L_{s}$ are the luminosities of the component stars and $r_{\rm pl-s}$ is the position of the planet with respect to the secondary component star (recall that we are interested in how the habitable zone of the lower-mass secondary star may be enlarged, so we use the distance from the hypothetical planet to the secondary star). If the binary is eccentric ($e>0$), $l_{\rm x-bin}$ will change as the two component stars orbit each other.

3 Results

We first determine the number of hard, or close, binary systems that have their orbital semimajor axis significantly hardened during dynamical interactions in our simulated star-forming regions and examine the distribution of these altered systems. In Fig. 1 we show the change in semimajor axis of binaries against their initial semimajor axis. Systems that are hardened are shown in red, whereas softened systems – those whose semimajor axis increases due to a dynamical interaction – are shown in blue. In the following we will focus on binaries with semimajor axes that are typically in the range 1 – 20 au as these are systems that can (a) host stable planets in their habitable zones and (b) could be altered such that the habitable zones of one or both stars could be enlarged. However, we have performed our calculations for all surviving binaries in our star-forming regions.

Clearly, significantly more binary systems in this range are hardened to smaller semimajor axes than are softened to larger semimajor axes. Indeed, for every one binary that has its semimajor axis softened by more than 10 per cent of the initial value, two binaries are hardened by more than 10 per cent. When binary stars are subject to dynamical processing it is often the orbital eccentricity that changes the most. This is apparent in Fig. 1, where we show the change in eccentricity and semimajor axis in one simulation. The initial $t$ = 0 Myr values are shown by the asterisks and the final $t$ = 10 Myr values are shown by the circles. We also show intermediate values at 0.1, 1, 2 and 5 Myr by the filled symbols; however because these systems are dynamically hard, most binaries experience only one perturbative interaction, and this interaction usually occurs within the first 0.1 Myr when the star-forming region is at its densest. The coloured lines denote whether the binary is hardened or softened. Most systems that are dynamically hardened have their eccentricities increased, which results in the two components of the binary becoming even closer at periastron.

In Fig. 1 we show the change in the size of the habitable zone for binaries that are hardened or softened in all simulations, as a function of the change in the semimajor axis of the binary. From a total of 7032 binaries in twenty realisations of the same star-forming region (352 per region), 4745 survive (237 per region) and of these, 354 are hardened or softened so that the habitable zone of the binary increases (18 per region).

We then consider systems where the habitable zones are enlarged such that the habitable zones around the binary components overlap. Across our ensemble of simulations, we find that 31 systems (i.e. 1–2 per region, or 0.44 per cent of the total number of initial binaries [0.65 per cent of the surviving binaries]) are dynamically altered such that their habitable zones expand and overlap. Most of these systems only overlap at periastron (denoted by the symbols within the black circles in Fig. 1), but we found two systems whose habitable zones even overlap at apastron (the symbols within the black squares). In the majority of systems, the habitable zone increases in size. Across all of our simulations, we find two systems where the habitable zone decreases in size; these systems are shown by the downward arrow symbols in Fig. 1.

We show one example of an overlapping habitable zone in Fig. 2, and typically one to two binaries in each star forming region undergo significant habitable zone enlargement. The binary system shown in Fig. 2 consists of a primary star of mass $m_{p}=0.99$M⊙ and a secondary star of mass $m_{s}=0.63$M⊙, initially on an orbit with semimajor axis $a=6.4$ au and an orbital eccentricity of $e=0.32$. The respective effective temperatures of these two stars are $T_{{\rm eff},p}=5780$K and $T_{{\rm eff},s}=4410$K. The binary then experiences an interaction that hardens the system to a final semimajor axis of $a=5.4$ au, but which also increases the eccentricity to $e=0.59$.

Before this interaction, the binary component stars each have distinct, separate habitable zones, with the lower-mass component having a noticeable smaller habitable zone than the higher-mass (primary) component (Fig. 2(a)). Following the interaction, when the binary is at periastron, the habitable zones are enlarged and also merge together (Fig. 2(a)).

We check the dynamical stability of the binary system before and after the interaction that causes the habitable zones to merge, using the criterion derived by Holman & Wiegert (1999) for ‘S-type’ orbits. For the system shown in Fig. 2, a planet will be stable if it is within 1.24 au of the star before the interaction, which decreases to 0.55 au following the interaction. The narrow habitable zone for the secondary star before the interaction ranges between 0.32 – 0.79 au, and is then 0.46 – 1.26 au afterwards, which means a planet on an orbit between 0.46 and 0.55 au would both be stable and orbit in the enlarged habitable zone. These stability criteria are shown in Fig. 2 by the dotted grey lines (planets will be stable if they orbit interior to these lines).

We determine the Holman & Wiegert (1999) stability criterion for each binary where the habitable zones are enlarged, and find that for 20 per cent of these systems a planet would only be stable if it lay interior to the habitable zone around either star. However, we have assumed that the planets are on circular orbits, whereas planets may be on eccentric orbits if they had already formed before the binary underwent its perturbative interaction. Because most of the hardening interactions occur in the first 0.1 Myr, a planet that is forced onto an eccentric orbit by a perturbative interaction may subsequently be circularised due to interactions with the remaning protoplanetary disc.

4 Discussion and conclusions

We have shown that dynamical interactions in dense star-forming regions can dynamically alter the orbits of relatively close ($\sim$5 – 10 au) binary stars, pushing the two component stars closer together and often increasing the eccentricity of the system. If the binary orbit is sufficiently ‘hardened’, i.e.  to around 5 au, the habitable zone around a low-mass secondary star ($\sim$0.5 M⊙) can be enlarged due to the increase in flux from the (typically) solar-mass primary star.

Depending on the orbit and the individual masses of the component star, the habitable zones overlap and are enlarged when the binary is at periastron (i.e. the stars are closest). Almost all of our systems do not have overlapping habitable zones at apastron, but the habitable zone of the lower mass star is still enlarged by the flux from the higher mass star. It is worth emphasising that any planet in the enlarged habitable zone would experience periods of lower habitability when the component stars were widely separated, unless the atmosphere of the planet can act as a buffer for the flux incident at periastron (Kaltenegger & Haghighipour, 2013). When the habitable zones do overlap at apastron, any planets would spend more time in the habitable zone.

We assume that the planets are on circular orbits, but they would likely be pushed onto eccentric orbits if they formed before the interaction that hardened the binary’s orbit. We are unable to quantify the number of systems that this could occur in, but we note that in such a scenario eccentric planets could pass in and out of the habitable zone.

Dense star-forming regions of the type we model here are usually considered to be extremely hostile to planet formation. Extreme and far-ultraviolet (EUV and FUV) radiation fields generated by OB-type stars have been shown to photoevaporate the gas content of protoplanetary discs in dense star-forming regions (Adams et al., 2004; Scally & Clarke, 2001). We would therefore expect the discs from which any planets form to be primarily composed of dust and unlikely to form gas or ice giants. However, given the closeness of the component stars in the binary, we would anyway expect that the planetary system would consist solely of terrestrial planets (Haworth et al., 2018).

We also note that we have assumed that the protoplanetary discs are circumprimary and circumsecondary, whereas in close binaries they may be circumbinary.

Interactions that would alter the orbit of binary systems would also likely disrupt the orbits of fledgling planetary systems (e.g. Parker & Quanz, 2012; Cai et al., 2017), and planets orbiting binary stars would likely be most susceptible due to a higher cross section for collisions (Adams et al., 2006). However, in our simulations the majority of the perturbative interactions occur in the first 1–2 Myr of the star-forming regions’ evolution, which is shorter than the timescale for terrestrial planet formation.

This is helped by our simulated star-forming regions rapidly dissolving in the first few Myrs, placing the binary systems into the Galactic field. As such, our hard binaries typically only experience one hardening encounter. In less realistic models for the initial conditions of star forming regions (e.g. Plummer (1911) or King (1966) profiles), the star-forming reigons would be longer lived and so the binaries would experience further hardening or softening encounters.

Soft binary systems are highly likely to be disrupted in star-forming regions; however, habitable zones around the component stars in these systems would be distinct, and for our purposes would be treated as the habitable zones around single stars. In the dense star-forming regions we simulate hard binaries (5 – 10 au) are much more likely to be hardened, increasing the size of the habitable zone around both stars.

In summary, we find that dense stellar environments can enlarge the habitable zones around the components of close binary star systems, without detrimentally affecting the formation process for terrestrial planets that may reside in the habitable zone.

Acknowledgements

We are grateful to the referee, Maxwell Xu Cai, for helpful comments and suggestions which improved the original manuscript. We thank Stuart Littlefair, Chris Watson and Belinda Nicholson for helpful discussions. RJP acknowledges support from the Royal Society in the form of a Dorothy Hodgkin Fellowship.