Delayed Photons from Binary Evolution Help Reionize the Universe

Amy Secunda; Renyue Cen; Taysun Kimm; Ylva Götberg; Selma E. de Mink

doi:10.3847/1538-4357/abaefa

1. Introduction

High-redshift, star-forming, dwarf galaxies with virial masses that range from 10⁸ to 10^10.5 M_⊙ are the most plausible source of the hydrogen-ionizing radiation responsible for the reionization of the universe by z ∼6, provided that the escape fraction (f_esc) is sufficiently high (Haehnelt et al. 2001; Cowie et al. 2009; Fontanot et al. 2014; Madau & Haardt 2015; Madau & Fragos 2017). Models of cosmic reionization (e.g., Kuhlen & Faucher-Giguère 2012; Shull et al. 2012; Robertson et al. 2015) show that f_esc ≳ 20% is required for cosmic reionization. The Thompson optical depth Kuhlen & Faucher-Giguère (2012) calibrated their results against was a bit higher than current estimates by Planck, and the exact requirement for f_esc depends on uncertain parameters such as the clumping factor of the intergalactic medium (IGM) and the intrinsic ionizing luminosity density of the early universe. Nonetheless, f_esc ≥10%–20% is probably required.

However, observations of Lyman continuum (LyC) in the local universe, which is limited to starburst galaxies, generally suggest low escape fractions of ≲8% (Leitet et al. 2011, 2013; Borthakur et al. 2014; Izotov et al. 2016; Leitherer et al. 2016). Star-forming galaxies at z ∼ 1 with LyC detections also show low escape fractions of a few percent (Siana et al. 2007, 2010; Bridge et al. 2010; Rutkowski et al. 2016). Although observations at redshifts z ≳3 (e.g., Leethochawalit et al. 2016; Reddy et al. 2016; Steidel et al. 2018) seem to suggest that galaxies with f_esc ≳ 10% are more common than at lower redshifts, only a few of these galaxies with high values of f_esc have been shown to be robust detections uncontaminated by foreground sources (Mostardi et al. 2015; Shapley et al. 2016). Interestingly, the average relative escape fraction, the ratio of the escape fraction of 990–1500 Å photons, is found to be small (most often ≲5%) in nearly all observations, even at high redshift (Vanzella et al. 2010; Boutsia et al. 2011; Grazian et al. 2015; Mostardi et al. 2015; Siana et al. 2015; Marchi et al. 2017); however, see also Rivera-Thorsen et al. (2019).

Theory also struggles to provide a value of f_esc ≳ 10%–20%. Recent high-resolution (≲10pc) numerical simulations that include feedback (e.g., Kimm & Cen 2014; Ma et al. 2015; Xu et al. 2016) suggest average escape fractions of only 5%–11%. In these simulations, the escape fraction is lowest during periods of high star formation (i.e., the periods in which LyC is most often observed in the local universe) before feedback has the opportunity to clear the LyC-trapping gas from the birth clouds of the stars.

Wise et al. (2014) suggested that ionizing photons from low luminosity (M_UV > −13) mini-halos could play a large role in reionization. However, Kimm et al. (2017) found that, although the escape fractions of these mini-halos are generally large (∼20%–40%), the inefficiency of star formation in these mini-halos means they play only a minor role in reionization. Using a simple analytic argument, Conroy & Kratter (2012) showed that including runaway O/B stars may enhance f_esc by a factor of up to ∼4.5 for halos with virial masses between 10⁸ and 10⁹ M_⊙. However, Kimm & Cen (2014) and Ma et al. (2015) found that including a simple model for runaway O/B stars in their simulations only increased their mean values of f_esc to 14% and ∼6%, respectively.

One of the most important theoretical findings on dwarf galaxies at high redshift is that star formation is very episodic. Peaks in star formation and f_esc tend to be out of phase by 10–30 Myr, with the latter lagging the former (Kimm & Cen 2014). When star formation is most vigorous in dense regions, the produced ionizing photons do not easily escape. On the other hand, when supernova feedback has cleared out channels in the ISM, the O/B stars produced at the peak of star formation that dominate LyC production (e.g., Leitherer et al. 1999) are gone.

Thus, potentially the most promising additional sources of H-ionizing radiation are stars that interact in binaries. Spectroscopic observations of young massive stars in the Milky Way and Magellanic Clouds, show that a very large fraction (≳0.7) of stars are part of close binary systems (e.g., Kobulnicky & Fryer 2007; Mason et al. 2009; Chini et al. 2012; Sana et al. 2012; Almeida et al. 2017; Moe & Di Stefano 2017). In these binary systems interactions between the stars can lead to the exchange of mass and angular momentum through Roche-lobe overflow, common envelope evolution, and the merging of the two stars (Kippenhahn & Weigert 1967; Paczynski 1976; Wellstein & Langer 1999; Ivanova 2011; de Mink et al. 2013; Schneider et al. 2016; De Marco & Izzard 2017). These interactions can increase the number of high-mass stars at later times and can create envelope-stripped helium stars, which both emit ionizing photons tens to hundreds of Myr after a starburst (Van Bever et al. 1999; Götberg et al. 2019). These "delayed" ionizing photons could be particularly effective in increasing f_esc, because there would be more time for feedback from massive stars to remove most of the surrounding gas from the birth cloud that would normally trap LyC radiation (Götberg et al. 2020).

Ma et al. (2016, hereafter M16) studied the effect binary stellar evolution had on the escape fractions of three example halos from the Feedback in Realistic Environments project (fire; Hopkins et al. 2014) using models from the Binary Population and Spectral Synthesis code (bpass; Eldridge et al. 2008, 2017; Eldridge & Stanway 2009, 2011; Stanway et al. 2016), which accounts for the mass transfer and mergers of a binary star population. They found that, while their stellar population synthesis model that did not include binary stellar evolution produced an f_esc below 5% at most redshifts for their example halos, the inclusion of binaries in their stellar population synthesis model increased f_esc by factors of ∼3–5. Including binaries in their model also increased the amount of ionizing photons by a factor of 1.5, leading to a factor of ∼4–10 increase in the "effective" escape fraction. Rosdahl et al. (2018, hereafter R18), also used a model from bpass that includes binary stellar evolution in their sphinx suite of cosmological adaptive mesh refinement (AMR) simulations. R18 found a similar increase in their photon-weighted mean escape fraction by a factor of ∼3 from f_esc ∼ 2%–3% for a single-stellar population synthesis model, to f_esc ∼ 7%–10% for their run that included binary stellar evolution.

In this paper, we focus on, separately, two different products of interactions between stars in binary systems: namely, envelope-stripped helium stars (see Section 2.2.2), and stars that gain mass via mass exchange in binaries and become massive blue stragglers (see Section 2.2.3). A significant advantage of this approach is that we can investigate the effect on the escape fraction from different types of binary products, which allows us to better understand which sources are responsible for any change in the escape fraction. Equally important, it is easier to see how the uncertainties in the different aspects of our models produce uncertainties in the emission rate of ionizing photons.

We perform post-processing on the cosmological radiation-hydrodynamic simulations from Kimm & Cen (2014) (hereafter, KC14) using the ionizing photon rates of stripped stars from Götberg et al. (2019) and a simple model for massive blue stragglers. Our aim is to better understand the role of stripped stars and massive blue stragglers during reionization and their effect on the escape fraction. The outline of the paper is as follows: in Section 2.1, we briefly describe the cosmological simulations used. In Section 2.2, we describe the implementation of the various stellar populations included in our model, including massive stars (Section 2.2.1), stripped stars (Section 2.2.2), blue stragglers (Section 2.2.3), and runaway stars (Section 2.2.4). In Section 2.2.5, we discuss other sources of ionizing radiation that we have not included in our simulations. In Section 2.2.6, we compare our stellar population synthesis models to the models M16 and R18 used from bpass. In Section 2.3, we describe our calculation of f_esc. In Section 3, we show our results for a larger mass example halo (Section 3.1) and all of the halos in our simulation (Section 3.2). In Section 4, we compare our results with two previous studies. Finally, in Section 5, we summarize our results.

2. Method

2.1. Cosmological Simulations

We perform our calculation of the escape fraction through post-processing of the "FRU" cosmological simulations of KC14, generated using ramses cosmological AMR code (Teyssier 2002). This enables us to make direct comparisons to the results presented in KC14 that do not include the two effects due to binary stellar evolution. Here, we briefly summarize the key components of these simulations.

KC14 use the music software (Hahn & Abel 2011) to generate the initial conditions with the WMAP7 cosmological parameters (Komatsu et al. 2011): (Ω_m, Ω_Λ, Ω_b, h, σ₈, n_s = 0.272, 0.728, 0.045, 0.702, 0.82, 0.96). They first run a simulation of dark matter only, using a sufficiently large box of size (25 Mpc h⁻¹)³ with 256³ dark matter particles in order to sample the large-scale gravitational field and include effects from the large-scale tidal field in the zoom-in simulation in the next step.

Next, a zoomed-in region of comoving size of 3.8 × 4.8 × 9.6 Mpc (comoving) is chosen, where a finer spacing is implemented in the initial condition to achieve a dark matter particle mass resolution of 1.6 × 10⁵ M_⊙, effectively corresponding to 2048³ over the entire box. The zoomed-in region is dynamically refined according to a preset criterion to better resolve the structure of the ISM, with up to 12 additional levels of refinement. This refinement results in a minimum physical cell size of 4.2 pc and a stellar mass resolution of approximately 49 M_⊙. The amiga halo finder (Gill et al. 2004; Knollmann & Knebe 2009) was used to identify dark matter (sub) halos and galaxies.

The AMR code, ramses (Teyssier 2002), is based on the fully threaded oct-tree structure (Khokhlov 1998) and uses the second-order Godunov scheme to solve the Euler equations. The hydrodynamic states reconstructed at the cell interface are limited using the MinMod method and then advanced using the Harten–Lax–van-Leer contact wave Riemann solver (Toro et al. 1994). KC14 use a typical Courant number of 0.8 and solve the Poisson equations using the adaptive particle-mesh method. Star formation and stellar feedback from supernova (SN) explosions are included in these simulations as outlined in KC14. Also included are radiative cooling (Sutherland & Dopita 1993; Rosen & Bregman 1995) and thermal stellar winds as in starburst99 (Leitherer et al. 1999).

KC14 also use the multigroup radiative transfer (RT) module developed by Rosdahl et al. (2013) to follow RT of ionizing photons from stars through the ISM and the IGM. The module solves the moment equations for H i, He i, and He ii ionizing photons using a first-order Godunov method with M1 closure for the Eddington tensor. KC14 adopt a Harten–Lax–van-Leer (Harten et al. 1983) intercell flux function, and use a photon production rate corresponding to a Kroupa IMF (Kroupa 2001) from the starburst99 library (Leitherer et al. 1999).

2.2. The Ionizing Emission from a Stellar Population

In addition to massive stars from a single-stellar population synthesis model, we account for two types of products of binary interaction: stars stripped of their hydrogen-rich envelopes (see Section 2.2.2) and stars that gain mass both from mass transfer and from coalescence (see Section 2.2.3). The stripped stars are the exposed hot and compact cores of their progenitors. In cases where their progenitors were high-mass stars, the stripped stars are Wolf-Rayet stars. The mass-gainers are rejuvenated and therefore appear bluer than the rest of the population. They are therefore often referred to as massive blue stragglers (e.g., Van Bever et al. 1999; Chen & Han 2009, 2010).

Note that we are replacing the stellar population synthesis models in post-processing, whereas the photoionization in the cosmological simulations is calculated using only the starburst99 models (see Section 2.1). If the photoionization in the cosmological simulation was calculated based on a stellar population that included binary stars, photons from these binary stars that failed to escape their host galaxy would ionize the gas that captures them. If these photons were able to ionize gas, it would increase the escape fraction in our simulations. We plan to include these effects on the fly in our future simulations.

The emission rates of hydrogen-ionizing photons from massive stars, stripped stars, and massive blue stragglers for a starburst of initially 10⁶ M_⊙ are shown as a function of time in Figure 1. We also show the ionizing photon emission rates from the model from bpass that includes binary stellar evolution that M16 and R18 used, for comparison (Stanway et al. 2016; Eldridge et al. 2017, version 2.0). Because most star particles in our simulation of high-redshift galaxies are very metal-poor, we show the ionizing photon emission rates for the lowest metallicity available for each stellar population. The lowest metallicity of each population is Z = 2 × 10⁻⁴ for the stripped stars, Z = 4 × 10⁻⁴ for the massive single stars and massive blue stragglers, and Z = 10⁻³ for the model from bpass. In our simulation, a star particle with Z ≲ 2 × 10⁻⁴ would be given the emission rates shown in Figure 1.

**Figure 1.** Photon production rate (s⁻¹) of LyC photons for a single 10⁶ M_⊙ starburst for the massive star component (in purple), stripped star component (in green), and massive blue straggler component (in orange). For comparison, we also include the LyC photon production rate of the model including binary stellar evolution that M16 and R18 used from BPASS (Stanway et al. 2016; Eldridge et al. 2017, version 2.0). We show the photon production rates for the lowest metallicity available for each population. In our simulation, a star particle with a metallicity ≲2 × 10⁻⁴ would be given the production rate shown here. After about 10 Myr, the photon production rate of the stripped stars becomes the dominant photon rate. The photon production rate of the massive blue stragglers peaks at about 5 Myr after the initial starburst and becomes larger than the photon production rate of the starburst99 component at around 15 Myr. Including photons from stripped stars increases the total amount of ionizing emission by 3.6%, and including photons from massive blue stragglers increases the total amount of ionizing emission by 3.3%.
Download figure:
Standard image High-resolution image

Each star particle in the ramses simulation has an age and metallicity. For each stellar population, we interpolate between the different available ages and metallicities to assign an ionizing photon production rate to each star particle. Below, we describe each model used in our simulations for the individual types of stars in more detail.

2.2.1. Massive Stars

We model the ionizing emission from massive stars using starburst99 (Leitherer et al. 1999, 2014). We use the stellar models from Padova (Bertelli et al. 1993, 1994; Marigo et al. 2008), cmfgen (Hillier & Miller 1998), and WM-basic (Pauldrach et al. 2001) atmospheric models, and assume a Kroupa (Kroupa 2001) initial mass function that stretches from 0.1 up to 100 M_⊙. We use the metallicities Z = 0.02 (solar), 0.008, 0.004, and 4 × 10⁻⁴ that are available in starburst99.

The dominating source of ionizing photons in starburst99 are the massive O/B-type main-sequence stars. The most massive stars die after a few Myr. Their deaths lead to the decline in the emission rate of ionizing photons seen in Figure 1, as the remaining stars are cooler and less luminous, making them less efficient ionizing sources.

Although starburst99 is created for stellar populations containing only single stars, it does a reasonable job representing the emission from a population that also contains binary stars on the main sequence. It does a reasonable job because binary interaction primarily occurs when the most massive star in the system has depleted hydrogen in its center and is expanding to become a red (super)giant (de Mink et al. 2008). During this phase, the star's radius expands significantly more than during the star's previous main-sequence evolution, making interactions with companions more likely to occur. However, mass gain through accretion or coalescence can affect main-sequence stars as well. For example, we do not subtract the contribution to the binary product from the remaining main-sequence star, meaning that the contribution from starburst99 is somewhat overestimated. However, the binary product from these processes is a significantly brighter and hotter main-sequence star (see Section 2.2.3), and we therefore consider the overestimate of ionizing emission from starburst99 to be negligible.

Our default calculation uses only the ionizing emission predicted by starburst99. We refer to this simulation as the SB99 run.

2.2.2. Stars Stripped in Binaries

We account for the ionizing emission from stripped stars in stellar populations of different metallicities by adopting the emission rates of ionizing photons presented by Götberg et al. (2019).⁵ Götberg et al. (2019) calculated the number and type of stripped stars present in stellar populations as a function of time. Convolving with the emission rate of ionizing photons from individual stripped stars of given masses (Götberg et al. 2018), they then computed the total emission rate of ionizing photons from stripped stars in stellar populations. The models are created using an initial mass function from Kroupa (2001) with lower and upper mass limits of 0.1 and 100 M_⊙ and for the metallicities Z = 0.014, 0.006, 0.002, and 2 × 10⁻⁴.

In Figure 1, the stripped stars (green, dashed curve) start contributing LyC emission about 10 Myr after a starburst, and become the dominant source of LyC photons soon after. This delay corresponds to the main-sequence lifetime of a 20 M_⊙ star, which is the most massive and thus shortest-lived star that Götberg et al. (2019) consider as progenitors for stripped stars. Over 1 Gyr, stripped stars at this metallicity (Z = 2 × 10⁻⁴) contribute an additional 3.6% of the total LyC emission from massive stars. Stripped stars contribute a higher percentage of the total number of LyC photons at higher metallicities. For example, at solar metallicity over 1 Gyr, stripped stars contribute an additional ∼6% of the total LyC emission from massive stars. The relative size of the stripped star contribution changes because, while massive stars emit fewer ionizing photons at higher metallicities, in our models the emission rates of the stripped star stellar population remains reasonably constant as a function of metallicity.

In one of our simulations, we account for the ionizing emission from massive stars and stars stripped in binaries. We refer to this run as the SSB run.

2.2.3. Massive Blue Stragglers

We take a simple approach when estimating the ionizing contribution from stars that accrete mass or merge during interaction with a binary companion. During mass transfer, the accreting star evolves on its main sequence. When it gains mass from the donor star, it gets rejuvenated as the convective core grows into the hydrogen-rich layers when the star increases in mass. This means that the emission rate of ionizing photons from a star that has accreted mass can be modeled as a younger and more massive main-sequence star.

The result of a merging binary star is most probably a more massive and rejuvenated main-sequence star in the case in which the merger occurred during the main-sequence evolution of both stars (Tout et al. 1997; Schneider et al. 2019). If the merger occurred when one of the stars is in a later evolutionary phase, the outcome of the merger is uncertain, but it is possible that the star expands and becomes a red supergiant. Since red supergiants are inefficient emitters of ionizing photons, we only account for mergers that occur during the main-sequence evolution of the two stars in the binary. We take the same approach when modeling the ionizing emission from blue stragglers resulting from mass accretion and coalescence.

We make similar assumptions to previous work on blue stragglers (e.g. Leigh et al. 2007, 2013; Leigh & Sills 2011). First, we assume that the interaction occurs when the main-sequence lifetime of the most massive star in the system has passed. Second, we assume that the blue straggler is rejuvenated such that it is halfway through its main-sequence evolution after the rejuvenation. Third, we assume that the blue straggler becomes 50% more massive than what the most massive star in the system was initially. Last, we assume that 10% of the stars with initial primary masses between 2 and 50 M_⊙ in a population go through this evolutionary phase and become blue stragglers.

The ionizing emission in starburst99 comes primarily from main-sequence stars, but also with some contribution from Wolf-Rayet stars. It is a reasonable assumption that mass-gaining stars evolve in similar ways as single stars, at least for the duration of the main sequence. We therefore use the predictions from starburst99 to estimate the emission rates of ionizing photons from mass-gaining stars via the following method. Using lifetimes of massive stars from evolutionary models (e.g., Figure 1 of Zapartas et al. 2017), we infer what the emission rate of ionizing photons from each mass bin of stars is via interpolation and subtracting the contribution from lower-mass stars. We calculate what the contribution from each mass bin should be following the assumptions described above. We then shift it with the time delay corresponding to the time of interaction.

The resulting emission rate of LyC photons can be seen in Figure 1, marked with a dashed–dotted, yellow line. The mass-gaining stars begin contributing ionizing photons after several Myr have passed. After around 20 Myr, the ionizing contribution of these mass-gainers becomes more significant than that of the stars in the population that did not interact. Since the mass-gainers are more massive than the single stars that are present in the stellar population, they also mildly harden the ionizing emission from the stellar population. Over 1 Gyr at a metallicity of Z = 4 × 10⁻⁴, our blue straggler population contributes a total of 3.3% of the ionizing radiation contributed by the single-stellar population. The contribution from blue stragglers decreases slightly at higher metallicities. For example, at solar metallicity, blue stragglers contribute a total of 2.1% of the ionizing radiation contributed by the single-stellar population.

We refer to the run where we include ionizing radiation from massive stars and massive blue stragglers as MBS.

2.2.4. Runaway Stars

Tetzlaff et al. (2011) found that roughly 30% of O/B stars are runaway stars with peculiar velocities greater than 28 km s⁻¹. They found that these peculiar velocities fall along a Maxwellian distribution with a dispersion of 24.4 km s⁻¹. Motivated by Tetzlaff et al. (2011), in their "FRU" simulations KC14 divided 30% of their star particles into runaway O/B stars. They drew the peculiar velocities of these runaway stars from a Maxwellian distribution with a dispersion of 24.4 km s⁻¹ and a minimum space velocity of 28 km s⁻¹, and added them to the initial velocity of the star particle. The directions of motion for the runaway stars were chosen at random.

The star particles that represent runaway stars are given the same ionizing photon rates as the other star particles, including the photon rates for stripped stars and massive blue stragglers in the SSB and MBS runs. There are currently two proposed channels for the formation of runaway stars. The first way a star could become a runaway is through a three-body interaction with other stars in a young cluster (Leonard & Duncan 1988), and the second way is from an SN explosion of a companion in a binary system (Blaauw 1961). The second channel would likely result primarily in runaway massive blue stragglers or main-sequence stars, depending on the level of interaction between the stars in the binary system before ejection. It could also potentially result in a runaway stripped star (Pols 1994; Renzo et al. 2019). Renzo et al. (2019) even found that stripped stars are just as likely to become runaways as massive blue stragglers. On the other hand, if fewer stripped stars become runaways, then fewer stripped stars will end up toward the outskirts of a galaxy where their emission would more easily be able to escape into the IGM. Therefore, fewer runaway stripped stars would result in a lower escape fraction. However, because the first channel for the creation of runaway stars would be impartial to stellar type, and because the relative importance of each channel is still not well-constrained (e.g., Hoogerwerf et al. 2000, 2001; Guseinov et al. 2005), we consider this simple approach to be adequate.

2.2.5. Other Sources of Ionizing Radiation

There are other stellar sources of ionizing radiation that we do not yet account for, such as post-AGB stars (Stasińska et al. 2015; Byler et al. 2019), accreting white dwarfs (van den Heuvel et al. 1992; Woods & Gilfanov 2013, 2014; Chen et al. 2015), and X-ray binaries (Fragos et al. 2013; Schaerer et al. 2019; Senchyna et al. 2020). The ionizing emission from post-AGB stars appears after about 100 Myr, at a rate that is about five orders of magnitudes lower than that of the most massive stars (Byler et al. 2019). The accreting compact objects provide ionizing emission at a much lower rate than living stars, but their emission is also significantly harder (e.g., Lepo & van Kerkwijk 2013). Since massive stars, stripped stars, and blue stragglers are expected to emit ionizing photons at the highest rates, compared to these of other stellar ionizing sources, we choose to start with including only these stars.

We also do not account for the stellar rotation of massive stars in our stellar population synthesis models (Huang et al. 2010; Ramírez-Agudelo et al. 2013) or stars spun up from binary interaction (Dufton et al. 2011; Eldridge & Stanway 2011; de Mink et al. 2013). Stellar rotation has been predicted to lead to, for example, more luminous stars, more WR stars, and stars evolving chemically homogeneously (e.g., Maeder 1987; Cantiello et al. 2007; Ekström et al. 2012). Although many effects of rotation are interesting, since they predict an increased emission rate of ionizing photons (e.g., Levesque et al. 2012; Topping & Shull 2015; Kubátová et al. 2019), there is only circumstantial observational evidence for the existence of these processes. However, interested readers should refer to the work of Abdul-Masih et al. (2019) and circumstantial evidence from, e.g., Martins et al. (2008) and Hainich et al. (2015); see also Schootemeijer & Langer (2018) and Shenar et al. (2019).

2.2.6. Comparison to bpass

Both M16 and R18 use models from bpass version 2.0, which have initial conditions very similar to those of our models. However, Figure 1 shows that, even with the additional emission from stripped stars and massive blue stragglers combined, the model from bpass produces more LyC photons over the course of a Gyr than the models used in this paper. It does not appear as dramatically because of the logarithmic scale of Figure 1, but the most significant difference between bpass and our model occurs between 3 and 20 Myr. This boost in ionizing emission occurs because bpass assumes that high-mass accretor stars evolve chemically homogeneously at Z ≤ 0.004 (Eldridge & Stanway 2011).

bpass also predicts a higher, almost flat emission rate at very late times (>200 Myr) (see Stanway & Eldridge 2018). This difference could be related to the inclusion of post-AGB stars in bpass, and also possibly to our model's inclusion of gravitational settling in the atmospheres of the low-mass stripped stars that form in older stellar populations (see Götberg et al. 2018). Interestingly, in the most recent model from bpass (Eldridge et al. 2017, version 2.2.1), fewer delayed LyC photons are produced through binary interactions.

In addition, the models from bpass version 2.0 only go down to a metallicity of Z = 1 × 10⁻³. This metallicity is still high for very metal-poor high-redshift dwarf galaxies, and the metallicity of a model can have a significant affect on the amount of ionizing photons. For example, in our models the total fraction of ionizing photons that come from binary products over 1 Gyr increases from our lowest-metallicity models to our next lowest-metallicity models from 6% to 9%. The models we use here go down to either Z = 4 × 10⁻⁴ or Z = 2 × 10⁻⁴.

2.3. Calculating the Escape Fraction

The rate of photons that escape the virial sphere of each halo for each star particle in that halo, when accounting for absorption by hydrogen, helium, and dust, as a function of frequency, ν, is

$\begin{eqnarray}\begin{array}{rcl}{\dot{N}}_{\mathrm{esc},{\rm{i}}}(\nu ) & = & {\int }_{0}^{4\pi }d{\rm{\Omega }}\dot{{N}_{{\rm{i}}}}(\nu )\exp [-{\sigma }_{{\rm{H}}{\rm{\small{I}}}}(\nu ){N}_{{\rm{H}}{\rm{\small{I}}}}({\rm{\Omega }})\\ & & -{\sigma }_{\mathrm{He}{\rm{\small{I}}}}(\nu ){N}_{\mathrm{He}{\rm{\small{I}}}}({\rm{\Omega }})-{\sigma }_{\mathrm{He}{\rm{\small{II}}}}(\nu ){N}_{\mathrm{He}{\rm{\small{II}}}}({\rm{\Omega }})\\ & & -{k}_{\mathrm{ext}}(\nu ){\rm{\Sigma }}({\rm{\Omega }})],\end{array}\end{eqnarray} \tag{ 1 }$

where Ω is the solid angle; N_{H i}(Ω), N_{He i}(Ω), and N_{He ii}(Ω), are the neutral hydrogen and the neutral and singly ionized helium column densities output from the ramses snapshots; k_ext is dust extinction opacity; and Σ(Ω) is the surface density of dust along the line of sight also output from the ramses simulations. Here, $\dot{{N}_{{\rm{i}}}}(\nu )$ is the age- and metallicity-dependent number of photons per second emitted by an individual star particle, i, in a ramses snapshot at frequency, ν. The value of $\dot{{N}_{{\rm{i}}}}(\nu )$ is extrapolated from the age- and metallicity-dependent starburst99 spectrum, with the age- and metallicity-dependent photon rates from the stripped star or massive blue straggler spectrum added on for the SSB and MBS runs, respectively.

The neutral hydrogen, σ_{H i}(ν), neutral helium, σ_{He i}(ν), and singly ionized helium, σ_{He ii}(ν), cross sections are calculated as in Osterbrock (1989). The frequency-dependent dust extinction opacity, k_ext(ν), is extrapolated from the dust model for the metal-poor Small Magellanic Cloud in Weingartner & Draine (2001) and Li & Draine (2001).

The escape fraction for each star particle, ${f}_{\mathrm{esc},{\rm{i}}}$ , in each halo snapshot is then computed as the ratio

$\begin{eqnarray}&&{f}_{\mathrm{esc},{\rm{i}}}=\displaystyle \frac{{\int }_{13.6\mathrm{eV}/h}^{\infty }{\dot{N}}_{\mathrm{esc},{\rm{i}}}(\nu )d\nu }{{\int }_{13.6\mathrm{eV}/h}^{\infty }{\dot{N}}_{{\rm{i}}}(\nu )d\nu }.\end{eqnarray} \tag{ 2 }$

For each ramses snapshot of each halo in each run,

$\begin{eqnarray}&&{f}_{\mathrm{esc}}=\displaystyle \frac{{\displaystyle \sum }_{i}{f}_{\mathrm{esc},{\rm{i}}}{\dot{N}}_{{\rm{i}}}}{{\displaystyle \sum }_{i}{\dot{N}}_{{\rm{i}}}},\end{eqnarray} \tag{ 3 }$

where ${f}_{\mathrm{esc}}$ is the photon production rate–weighted mean escape fraction for each halo snapshot.

Note that our calculation of the number of escaping ionizing photons is done in post-processing using a simulation that did not include the effects of binary evolution. As a result, we predict that our value of f_esc is likely lower than it would be if it were calculated directly in the cosmological simulation. We expect that f_esc is lower because when we include photons from binary stars only in post-processing; photons from this population that fail to escape also do not ionize the gas that absorbs them. If these photons were included in the cosmological simulation, they would ionize the gas that absorbs them, helping to clear a path for future photons. It is also useful to note that, because our calculation is done in post-processing, the total number of escaping LyC photons from both additional sources (stripped stars and massive blue stragglers) combined is simply the sum of the number of escaping LyC photons from each source calculated separately.

3. Results

3.1. An Example Halo

Before presenting statistical results, we use one halo to illustrate (in Figures 2 and 3) some basic effects on the escape fraction of ionizing photon due to binary evolution. A relatively massive halo, for which we are able to construct a merger tree to a very high redshift, is used here as an example. The virial and stellar masses of this halo at z = 7 are 3.16 × 10¹⁰ M_⊙ and 2 × 10⁸ M_⊙, respectively.

**Figure 2.** Evolution of the escape fraction (f_esc) for an example halo as a function of the age of the universe in Gyr. SB99 run is shown in purple (solid line), SSB run in green (dashed line), and MBS run in orange (dashed–dotted line). Purple shaded region shows star formation rate (SFR) in ${M}_{\odot }\,{\mathrm{yr}}^{-1}$ . Redshifts corresponding to the ages of the universe shown on the lower x-axis are given on the upper x-axis along with the log of the stellar mass (in ${M}_{\odot }$ ) of the halo at that redshift. As noted in KC14, the peaks in f_esc are out of phase with the SFR.
Download figure:
Standard image High-resolution image

**Figure 2.** Evolution of the escape fraction (f_esc) for an example halo as a function of the age of the universe in Gyr. SB99 run is shown in purple (solid line), SSB run in green (dashed line), and MBS run in orange (dashed–dotted line). Purple shaded region shows star formation rate (SFR) in ${M}_{\odot }\,{\mathrm{yr}}^{-1}$ . Redshifts corresponding to the ages of the universe shown on the lower x-axis are given on the upper x-axis along with the log of the stellar mass (in ${M}_{\odot }$ ) of the halo at that redshift. As noted in KC14, the peaks in f_esc are out of phase with the SFR.
Download figure:
Standard image High-resolution image

KC14 found that more massive halos, on average, have lower escape fractions. The photon production rate–weighted escape fraction of this halo for the SB99 run is 11%. When photons from either stripped stars or massive blue stragglers are included, we see that the escape fraction increases significantly.

Figure 2 shows the escape fraction as a function of time, in Gyr on the bottom axis and redshift on the top axis. The SB99 run is shown in purple (solid line), the SSB run in green (dashed line), and the MBS run in orange (dashed–dotted line). This color scheme (and line type) is consistent for all figures in this section. The shaded purple region shows the star formation rate (SFR) in M_⊙ yr⁻¹ as a function of time, which is calculated by summing over star particles formed in the last 1 Myr. The logarithm of the stellar mass of the halo at each redshift is also shown on the upper axis above the redshift.

As shown first in KC14, the escape fraction decreases significantly when the SFR is highest due to the presence of large amounts of gas during star formation, which occurs deep in the cores of giant molecular clouds. There is therefore a delay between peak star formation and the increase in the escape fraction that occurs once much of the birth cloud has been cleared through supernova-driven blastwaves. This delay is commonly seen in dwarf galaxy simulations at high redshift (e.g., Wise & Cen 2009; Wise et al. 2014; Ma et al. 2015). For this halo, the escape fraction is always higher for both runs that include products of binary interactions than for the run that does not.

The largest increases in the escape fraction in the SSB and MBS runs compared to the SB99 run occur about 10 Myr after a major starburst has ended, for example at t_H ∼ 0.705 Gyr in Figure 2. At first, the escape fractions of all three runs increase, because supernova feedback has significantly cleared out the birth clouds of the starburst. However, shortly after the escape fraction for the SB99 run increases, it will rapidly decrease, because the most massive single stars are gone after 10 Myr. The escape fractions of the SSB and MBS runs decrease less, because LyC photon production for stripped stars and massive blue stragglers peaks at around 10 Myr.

Figure 3 compares the number of escaped photons as a function of redshift for the three runs. The top panel shows the cumulative number of escaped photons for SB99 (solid purple), SSB (dashed green), and MBS (dotted–dashed orange). The middle panel shows the ratio of the cumulative number of photons of SSB+SB99 to SB99 (dashed green), and that of MBS+SB99 to SB99 (dotted–dashed orange), respectively. The bottom panel is similar to the middle panel but for the instantaneous number of photons. The shaded purple region once again shows the SFR in M_⊙ yr⁻¹. We see that the instantaneous ratio of escaped photons changes significantly over time, often reverting back to being close to order unity directly after a starburst. For this example halo, we see that when stripped stars are included, 27% more LyC photons escape into the IGM by z = 7, compared to when only massive stars are accounted for. When massive blue stragglers are included, 17% more LyC photons escape. Combined LyC emission increases by 42% for this halo when both stripped stars and massive blue stragglers are included. This would increase the escape fraction of this halo to ∼14% when stripped stars are included, and ∼16% when both stripped stars and massive blue stragglers are included. The dip present at z ∼ 7.7 is from the boost in the starburst99 ionizing emission rate right after a new starburst.

3.2. Statistical Results

The cosmological simulations of KC14 include hundreds of halos ranging in virial mass from 10⁸ M_⊙ to 10¹¹ M_⊙ at z ∼ 7. The left panels of Figure 4 show the escape fraction (f_esc), in percent, as a function of virial mass. The right panels show the production rate of ionizing photons that escape the virial sphere (in number per second) as a function of virial mass. In the top three panels, each point represents an individual halo. The purple, green, and yellow lines show the median values for all halos within a mass range as a function of the median of those mass ranges for the SB99, SSB, and MBS runs, respectively. The error bars represent the interquartile range. The gray lines show the boxcar-smoothed, photon-weighted mean values for the three runs. The bottom of the gray line represents the mean for the SB99 run, and the top of the line represents the mean of the SSB run. The bottom panel is a zoom-in to show the boxcar-smoothed, photon-weighted means for the three runs at z = 7, with the same color scheme as in Figure 2. To increase the statistical significance, results from six consecutive snapshots are combined to determine the medians, interquartile ranges, and means for each redshift specified. The color scheme and line types are the same as above. The top panels are at z = 11, the panels directly below them are at z = 9, and the bottom two panels are at z = 7.

We find a median unweighted escape fraction of roughly 10%–50% for all halo masses, although there is a large interquartile range. Less massive halos ( ${M}_{\mathrm{vir}}\lt {10}^{9}\,{M}_{\odot }$ ) overall have a higher median escape fraction than higher-mass halos, as noted earlier in KC14. There is not a large variance in the escape fraction over different redshifts, except an increased smoothness in the median values because there are more halos at later redshifts.

The median production rate of escaping photons increases with virial mass, as more massive galaxies tend to have higher SFRs. The median production rates also decrease slightly with decreasing redshift, because for a given halo, the SFR decreases and the metallicity increases with redshift. The interquartile range is also large for the production rate of escaping photons, but the median values range from roughly 10⁴⁵ s⁻¹ at lower virial masses to 10⁵² s⁻¹ at higher virial masses.

In Figure 4, the median unweighted escape fraction is slightly lower for the SSB run than the SB99 run for halos with virial masses less than about 10⁹ M_⊙, but higher for the most massive halos. The median unweighted escape fraction for our MBS run is very similar to the median unweighted escape fraction for the SB99 run for halos with virial masses less than about 10⁹ M_⊙, and then slightly higher for the more massive halos. The rate of escaping photons is greater at all redshifts for all mass halos for the SSB and MBS runs, but by a larger factor for less massive halos than more massive halos. At low masses, this factor is around 40–200 for the SSB run and around 5–10 for the MBS run; at high mass, this factor is around order unity to 20 for the SSB run and 1–5 for the MBS run. These factors are slightly greater at later redshifts than earlier ones.

The differences in the median rate of escaping photons between the SSB and MBS runs and the SB99 run are greater at lower masses because the starbursts of low-mass halos are more episodic, meaning that many halos at these low masses have undergone very little recent star formation. During periods of low star formation, stripped stars—and to a lesser extent, massive blue stragglers—become the dominant source of LyC radiation, because of the relatively long time range over which these sources have a significantly higher LyC photon rate. For example, in Figure 1, the starburst99 photon rate of a 10⁶ M_⊙ starburst with an age of 100 Myr will be only ∼10⁴⁷ s⁻¹. Meanwhile, at that age, stripped stars will have a photon production rate roughly 200 times greater, which accounts for the difference of two orders of magnitude between the escaping photon rates of the SSB and SB99 runs in low-mass halos.

Stripped stars in particular have harder ionizing emission than massive stars, and so these low-mass galaxies whose LyC emission is dominated by stripped stars would also emit harder spectra. These low-mass galaxies emitting harder spectra may be observable in the future by the James Webb Space Telescope (JWST), or they could even produce a unique signature in future observations of reionization bubbles by the Nancy Grace Roman Space Telescope. Interestingly, if JWST observes these low-mass halos, then the observed ionizing emission would almost exclusively come from binary products, dominated by stripped stars (see right panel of Figure 4), which means that these low-mass halos could be very useful for studying binary interactions at high redshift. However, the dramatic increase in the photon rate for the SSB run for the bursty low-mass halos is also responsible for the slight decrease in the median escape fractions for the SSB run. Star particles that previously were not emitting any hydrogen-ionizing photons now emit some radiation that can still be trapped, which lowers the escape fraction.

The boxcar-smoothed, photon-weighted mean escape fractions, shown as the gray line in the top three left-hand panels of Figure 4, are mostly lower than the median escape fractions. The most massive halos occasionally have similar mean and median escape fractions because there are not many of them. The difference between the unweighted median and weighted mean escape fractions is not surprising since, as KC14 found and as can be seen in Figure 2, star formation peaks are out of phase with peaks in the escape fraction by 10–200 Myr. The effects of the star formation peak and escape fraction peak phase difference are somewhat mitigated by the delayed photons in the SSB and MBS runs leading to the higher weighted mean escape fractions for these runs. These higher escape fractions can be seen in the bottom panel of Figure 4. Photons from the stripped star component are particularly delayed, with the photon rate for this component remaining above 10⁴⁹ s⁻¹ for over ∼100 Myr after the initial starburst.

The bottom right panel of Figure 4 shows that, at z = 7, the photon-weighted mean production rate of escaping photons is greater for both runs that include binary effects, although the dramatic increase in the production rate for lower-mass halos is gone. In fact, the difference in $\langle {\dot{N}}_{\mathrm{ph}}{f}_{\mathrm{esc}}\rangle$ is largest for the most massive halos. The value of $\langle {\dot{N}}_{\mathrm{ph}}{f}_{\mathrm{esc}}\rangle$ for the lower-mass halos does not increase as dramatically between runs because these halos' photon rates are so much higher during the 20 Myr after a starburst that photons from young stellar populations completely dominate the mean rate of escaping LyC photons. Nonetheless, the fact that most low-mass galaxies will have their LyC emission rate so enhanced during a majority of their evolution while their SFRs are low provides an exciting opportunity to learn about binary products, if these galaxies could be observed.

Figure 5 shows the photon production rate–weighted mean escape fractions as a function of redshift, for all halos in our simulation (upper panel), halos with a virial mass less than 10⁹ M_⊙ at that redshift (middle panel), and halos with a virial mass greater than 10⁹ M_⊙ at that redshift (lower panel). The color scheme and line types are the same as in Figure 2. Overall, the mean escape fractions stay mainly within 10%–20%, with higher-mass halos having a somewhat larger variance. The large temporal variations are dominated by a handful of large starburst events. The weighted mean escape fraction at all redshifts for both the SSB and the MBS runs is greater than that for the SB99 run.

$\langle {f}_{\mathrm{esc}}\rangle $ — **Figure 5.** Photon production rate–weighted mean escape fraction ( $\langle {f}_{\mathrm{esc}}\rangle$ ), in percent, as a function of redshift. Top panel shows $\langle {f}_{\mathrm{esc}}\rangle$ for all halos in the simulation. Middle panel shows $\langle {f}_{\mathrm{esc}}\rangle$ for halos with virial masses less than 10⁹ M_⊙. Lower panel shows $\langle {f}_{\mathrm{esc}}\rangle$ for halos with virial masses greater than 10⁹ M_⊙. Color scheme and line styles are the same as in Figure 2. Both runs that include binary effects have an $\langle {f}_{\mathrm{esc}}\rangle$ greater than than that for the SB99 run over all redshifts in every panel.
Download figure:
Standard image High-resolution image

Figure 6 compares the total number of photons that escape the virial sphere for all halos as a function of redshift for the various runs. The top panel shows the cumulative number of escaping LyC photons for all halos in each run. The panel directly below it shows the ratio of these numbers of escaping photons for the SSB and MBS runs to the SB99 run in green and orange, respectively. The bottom two panels show the instantaneous ratios of the total number of escaping photons for all halos (third panel) and halos more massive than 10⁹ M_⊙ (bottom panel). The color scheme and line types are the same as in Figure 2. The dotted line in the bottom panel shows the average instantaneous ratio for all the halos, for comparison.

The inclusion of stripped stars and massive blue stragglers increases the total number of escaped photons at a redshift of ∼7 by around 10⁶⁸ photons and 8 × 10⁶⁷ photons, respectively (i.e., by respective factors of ∼1.125 and ∼1.10). Including both binary products would increase the total number of escaped photons by ∼2 × 10⁶⁸ (i.e., by 22.5%). After the initial high-redshift starbursts, the overall trend is for the ratio of escaped photons between the SSB and MBS runs and the SB99 run to increase as redshift decreases. This trend is also present in the ratios of the median production rate of escaping photons, as mentioned above in the discussion of Figure 4. This ratio increases because the SFR decreases with redshift. A lower SFR means there will be more older stellar populations at lower redshifts, and stripped stars and massive blue stragglers are the dominant sources of LyC radiation in these populations. If the escape fractions for the various runs were identical, then for star particles at a metallicity of 2 × 10⁻⁴, the additional photons from stripped stars and massive blue stragglers would lead to increases in the total numbers of photons by factors of 1.036 and 1.033, respectively. Therefore, LyC photons from stripped stars and massive blue stragglers are significantly more effective at escaping their host galaxies. In Section 4, we use Equation (5) to calculate that the mean escape fraction for photons from these binary products is 57%, versus 14% for the massive stars.

The instantaneous ratios between the SSB/MBS and SB99 runs of the total number of escaping ionizing photons at a given redshift vary between just above unity to 1.3 for all halos. When only considering halos more massive than 10⁹ M_⊙, these same ratios vary between just below 1.1 and 1.8, from z = 11 to z = 7. The instantaneous ratios of the more massive halos often stay above the mean value of the instantaneous ratio for all the halos (shown in the bottom panel as the dotted line). In the bottom right panel of Figure 4, $\langle {\dot{N}}_{\mathrm{ph}}{f}_{\mathrm{esc}}\rangle$ also increases the most for runs including binary products compared to the SB99 run for higher-mass halos. One reason these more massive halos will have a larger increase in the number of escaping photons is because larger halos that undergo more starbursts will have more stellar populations at ages of around 10–30 Myr since the starburst. Stellar populations at these ages will still have photon rates high enough to influence the total number of escaped photons, as well as significant and even dominant contributions from stripped stars and massive blue stragglers.

Furthermore, KC14 found that when a single-stellar population synthesis model is used, the most massive halos tend to have lower escape fractions. They attribute this to the fact that, in larger halos, young massive stars can be buried in many star-forming clouds that are more resilient to SN feedback arising in neighboring star clusters. Our results suggest that delayed photons have a greater chance of escaping these halos because they allow more time for additional SN feedback to occur—not just in the cloud where the radiation is coming from, but also in adjacent clouds that will have more time to undergo their own stellar evolution, which will generate feedback.

The massive example galaxy shown in Figures 2 and 3 (see Section 3.1) clearly has a greater-than-average increase in escaped LyC photons when stripped stars and massive blue stragglers are included. Because it is one of the more massive galaxies in the simulation, it fits the overall trend that more massive galaxies have a larger additional amount of escaping photons in our SSB and MBS runs.

4. Comparison to Previous Studies

M16 used the stellar population synthesis models from bpass that include binary evolution to calculate the escape fraction of LyC for three example mock halos. For a high-mass (M_vir = 10¹⁰ M_⊙) halo, similar in mass to our example halo in Section 3.1, M16 found that including binary evolution increased their escape fraction by a factor of ∼3 from ∼5% to ∼14%. For our SB99 run, the more massive halos, like our example halo, had smaller-than-average escape fractions of ∼11%. Combining the increase in the number of escaped photons from both the SSB and MBS runs for our example halo increases the total effective escape fraction by a factor of 1.44 from 11% to 16%.

R18 also used models from bpass. They found that adding binary stellar evolution to their simulations increased the luminosity-weighted mean escape fraction for their large sample of halos by a similar factor to M16 of ∼3 from f_esc = 2.7% to f_esc = 8.5%. The photon rate–weighted mean escape fraction over all of our simulated halos increased by a factor 1.225 from 14% to 17%.

Our single-stellar population synthesis run clearly produces significantly larger values of f_esc than the single-stellar population synthesis runs in M16 and R18. One reason for this is that our single-stellar population model includes runaway O/B stars (see Section 2.2.4), which KC14 showed to increase their overall value of f_esc from 11% to 14%. However, 11% is still a factor of 2–5 larger than the percent of photons escaping in the M16 and R18 simulations. This discrepancy is likely due to the different star formation and feedback schemes among the different simulations.

R18, in particular, use a varying star formation efficiency (star formation efficiency per dynamic time) that leads to preferential star formation in higher-density regions. Because our simulations use a fixed star formation efficiency of 2% per dynamical time, even though the density threshold for star formation in their simulation is lower than ours (10 cm⁻³ in theirs versus 100 cm⁻³ in ours), stars in our simulations are less abundant in the densest regions. These denser regions of gas where the majority of stars are forming in the R18 sphinx simulations are more difficult to clear with stellar feedback, which is probably the reason R18 have much lower values for their escape fractions with and without binary stellar populations, exacerbated by the omission of runaway O/B stars in their simulations.

Importantly, the galaxy stellar mass and luminosity functions in R18's simulations are consistent with observations but lie near the upper end of observational constraints and other recent models. Their results would still fall within the observational constraints even if the amplitudes of the stellar masses were decreased by as much as a factor of two through increased SN feedback. This stronger SN feedback would lead to an increase in f_esc and a decrease in star formation—and therefore ionizing luminosities. These changes would make the results of R18 more similar to ours here.

The different star formation treatment in R18 is also partially responsible for their much larger increase in f_esc when including binary stellar evolution in their simulation. Because the birth clouds in R18 are denser, it takes longer for feedback and LyC photons to clear them. Therefore, most massive stars will be gone by the time these clouds are cleared, and so very few LyC photons from this stellar population are able to escape. Delayed LyC photons from binary stellar evolution, on the other hand, will be able to escape. As a result, the increase in the escape fraction due to binary stellar evolution will be greater. The denser birth clouds in R18 also allow for a concerted launch of SN-driven blastwaves from a more concentrated, larger stellar cluster with nearly coeval star formation, leading to a "cleaner" sweep of the remaining gas in star birth clouds. This can be seen indirectly in the rather large changes in escape fraction for singles and binaries in Figure 13 of R18.

Most significantly, the escape fractions in M16 and R18 both increase with the inclusion of binary evolution by a greater factor than in this paper because the number of additional LyC photons from binary stellar evolution is greater in both of those papers than in this paper. The number of LyC photons increases by roughly 50% and 70% for M16 and R18, respectively, while here we only have a 6% increase. Figure 1 shows that the LyC photon ionization rate of the lowest-metallicity model from bpass that M16 and R18 use for their runs that include binary products is significantly greater than ours. As discussed in Section 2.2.6, this difference is likely mainly due to the chemically homogeneously evolving accretor stars in bpass. The high LyC photon rate in R18, may be the reason why they find that their simulations are reionized slightly earlier than observations would predict, despite a low escape fraction of 7%.

Because the M16 and R18 simulations have a greater increase in the number of LyC photons due to binary interaction, the escape fractions of the additional ionizing photons from these interacting binaries, ${f}_{\mathrm{esc},\mathrm{bin}}$ , are weighted more heavily in their overall escape fractions than in ours. For example, 40% of the ionizing radiation in the R18 simulation that uses the spectrum from bpass comes from interacting binaries. Therefore, if ${f}_{\mathrm{esc},\mathrm{ss}}$ is the escape fraction for the single star component, the overall escape fraction for the R18 simulation that uses the spectrum from bpass is:

$\begin{eqnarray}&&{[{f}_{\mathrm{esc}}]}_{{\rm{R}}18}=0.6{f}_{\mathrm{esc},\mathrm{ss}}+0.4{f}_{\mathrm{esc},\mathrm{bin}}.\end{eqnarray} \tag{ 4 }$

On the other hand, in this paper only 7% of the ionizing radiation comes from interacting binaries. The overall escape fraction here is:

$\begin{eqnarray}&&{f}_{\mathrm{esc}}=0.93{f}_{\mathrm{esc},\mathrm{ss}}+0.07{f}_{\mathrm{esc},\mathrm{bin}}.\end{eqnarray} \tag{ 5 }$

Clearly, ${f}_{\mathrm{esc},\mathrm{bin}}$ is more heavily weighted in R18.

It is therefore useful to compare ${f}_{\mathrm{esc},\mathrm{bin}}$ among the three papers. Physically, ${f}_{\mathrm{esc},\mathrm{bin}}$ represents the escape fraction of the additional LyC radiation from interacting binaries. We can calculate ${f}_{\mathrm{esc},\mathrm{bin}}$ using a rearranged generalized version of Equations (4) and (5),

$\begin{eqnarray}&&{f}_{\mathrm{esc},\mathrm{bin}}=\displaystyle \frac{{f}_{\mathrm{esc}}-{{xf}}_{\mathrm{esc},\mathrm{ss}}}{1-x},\end{eqnarray} \tag{ 6 }$

where x is the fraction of LyC radiation from the single-stellar component in simulations that include binary interaction. From each paper, we know x (the escape fraction of the run that includes only the single-stellar population), which we take as ${f}_{\mathrm{esc},\mathrm{ss}}$ ; we also know the overall escape fraction of the run that includes binary evolution, f_esc. This information will give accurate values of ${f}_{\mathrm{esc},\mathrm{bin}}$ for this paper and M16 because both apply a ray tracer in post-processing. However, using Equation (6) to calculate ${f}_{\mathrm{esc},\mathrm{bin}}$ for R18 is not entirely accurate, because they use a model from bpass in the cosmological simulation itself. Implementing a model from bpass in the cosmological simulation itself means that the extra photons from binary products can also help LyC photons from single stars to escape.

With this caveat for the value of ${f}_{\mathrm{esc},\mathrm{bin}}$ for R18 in mind, ${f}_{\mathrm{esc},\mathrm{bin}}=32 \% ,17 \% ,57 \%$ for M16, R18, and this paper, respectively. Both values of ${f}_{\mathrm{esc},\mathrm{bin}}$ for M16 and R18 are greater than their values of ${f}_{\mathrm{esc},\mathrm{ss}}$ by a factor of ∼6, which one would expect given that the radiation from single stars is trapped so well in M16 and R18. Our value increases from ${f}_{\mathrm{esc},\mathrm{ss}}$ to ${f}_{\mathrm{esc},\mathrm{bin}}$ by a slightly smaller factor of ∼4. In this paper, the value of ${f}_{\mathrm{esc},\mathrm{bin}}$ is much greater than those from both M16 and R18, once again because of the different feedback and star formation schemes of the various simulations. Interestingly, Götberg et al. (2020) found that if the escape fraction of stripped stars is high (∼50%–80%), they can provide most of the LyC emission required for reionization, despite massive stars having typical escape fractions of ∼10%–20%. Our simulations support that the high escape fraction Götberg et al. (2020) set in their simulations is achievable for binary products.

5. Conclusions

Using ultra-high-resolution cosmological radiation-hydrodynamic simulations, we study the effects of interacting binaries on the fraction of LyC photons that escape from their host galaxy. Binary evolution modeling inevitably remains uncertain due to lack of direct information of stellar binarity and the stellar initial mass function in galaxies at the epoch of reionization. However, utilizing results from our simple model of blue stragglers and state-of-the-art model of stripped stars based on local stellar observations, we demonstrate that the extra LyC photons from these two sources increase f_esc significantly.

LyC photons from these two sources are delayed relative to the initial starbursts, on timescales of ∼10 Myr. As a result, the photons released from stripped stars and massive blue stragglers increase the number of escaped photons by a more significant factor than just the number of additional photons they produce. For example, for a single 10⁶ M_⊙ starburst at a metallicity of ≲2 × 10⁻⁴, our SSB model produces a factor of 1.036 more photons than the SB99 model, and the MBS model produces a factor of 1.033 more photons than the SB99 model. These factors are significantly less than the factor of 1.125 and 1.10 additional escaping photons by z = 7 that we see in our SSB and MBS runs, respectively. The physical reason is that the escape fraction and the instantaneous SFR are often strongly anticorrelated with a phase shift of 10–200 Myr, reflecting the fact that star formation tends to occur in regions of high obscuration and that it takes roughly 10–200 Myr to clear the birth giant molecular cloud.

Our results also suggest that LyC photons from massive blue stragglers—and in particular, stripped stars—significantly dominate the ionizing emission of low-mass galaxies for a majority of these galaxies' histories. If these galaxies can be detected by JWST, they would make excellent laboratories for studying interacting binaries as a function of redshift. Future work is needed to better predict the spectral shape of these galaxies—as, for example, LyC radiation from stripped stars tends to be harder than that from typical massive stars, potentially masquerading as Population III stars.

Figure 7 summarizes the increase of the photon-weighted mean escape fraction due to the addition of photons from binary products, from the 14% found by KC14, who used a starburst99 spectrum and included runaway O/B stars in their simulations, up to 16% and 17% when stripped stars and both stripped stars and massive blue stragglers are included.

The requirement for cosmological reionization remains somewhat unclear due to a combination of uncertain factors, including the initial mass function and clumping factor of the IGM as a function of redshift. R18 found that an escape fraction as low as 7% would be enough to reionize the universe by z ∼ 6. However, 20% is often considered to be the minimum required escape fraction to balance the recombination rate at z ∼ 6, if the initial mass function at z ∼ 6 is not too different from its local counterpart and a clumping factor of ∼6 is adopted (e.g., Kuhlen & Faucher-Giguère 2012; Shull et al. 2012; Robertson et al. 2015). While 17% is still a few points lower than the quoted 20%, the most important implication is that our modeling, which includes new stellar physics, high resolution, and advanced treatment of supernova feedback, puts the stellar reionization picture on a more solid footing. Future work is still needed to better understand the connection between the escape fractions of high-redshift dwarf galaxies with those of their lower-redshift counterparts.

We would like to thank the anonymous referee for thoughtful feedback that strengthened and improved the clarity of the paper. A. S. would like to thank Charles Emmett Maher for useful discussions. A. S. is supported by the National Science Foundation Graduate Research Fellowship Program under grant No. #DGE1656466. T. K. was supported in part by the National Research Foundation of Korea (NRF-2017R1A5A1070354 and NRF-2020R1C1C100707911) and in part by the Yonsei University Future-leading Research Initiative (RMS2-2019-22-0216). Y. G. acknowledges funding from the Alvin E. Nashman fellowship for Theoretical Astrophysics. This research is supported in part by NASA grant 80NSSC18K1101. Computing resources were provided in part by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center, and in part by Horizon-UK program through DiRAC-2 facilities.

Delayed Photons from Binary Evolution Help Reionize the Universe

Article metrics

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction