Delayed Photons from Binary Evolution Help Reionize the Universe

High-resolution numerical simulations including feedback and aimed at calculating the escape fraction (fesc) of hydrogen-ionizing photons often assume stellar radiation based on single-stellar population synthesis models. However, strong evidence suggests the binary fraction of massive stars is 70%. Moreover, simulations so far yield values of fesc falling only on the lower end of the roughly 10-20% range, the amount presumed necessary to reionize the Universe. Analyzing a high-resolution (4 pc) cosmological radiation hydrodynamic simulation we study how fesc changes when we include two different products of binary stellar evolution - stars stripped of their hydrogen envelopes and massive blue stragglers. Both produce significant amounts of ionizing photons 10-200 Myr after each starburst. We find the relative importance of these photons are amplified with respect to escaped ionizing photons, because peaks in star formation rates (SFRs) and fesc are often out of phase by this 10-200 Myr. Additionally, low mass, bursty galaxies emit Lyman continuum radiation primarily from binary products when SFRs are low. Observations of these galaxies by the James Webb Space Telescope could provide crucial information on the evolution of binary stars as a function of redshift. Overall, including stripped stars and massive blue stragglers increases our photon-weighted mean escape fraction by around 13% and 10%, respectively, resulting in a mean fesc of 17%. Our results emphasize that using updated stellar population synthesis models with binary stellar evolution provides a more sound physical basis for stellar reionization.


INTRODUCTION
High-redshift, star-forming, dwarf galaxies with virial masses that range from 10 8 − 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 a high enough escape fraction (f esc ) (Haehnelt et al. 2001;Cowie et al. 2009;Fontanot et al. 2014;Madau & Fragos 2017;Madau & Haardt 2015). Models of cosmic reionization (e.g., Kuhlen & Faucher-Giguère 2012;Robertson et al. 2015) show that f esc 20% is required for cosmic reionization. The Thompson opti-cial 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(Leitet et al. , 2013Borthakur et al. 2014;Leitherer et al. 2016;Izotov 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(Siana et al. , 2010Bridge et al. 2010;Rutkowski et al. 2016). Although observations at redshifts z 3 (e.g., Reddy et al. 2016;Leethochawalit 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Å to 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;Mostardi et al. 2015;Siana et al. 2015;Grazian et al. 2015;Marchi et al. 2017), see however, 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 minihalos 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 8 and 10 9 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 source 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;Sana et al. 2012;Chini et al. 2012;Moe & Di Stefano 2017;Almeida et al. 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 since it would give time for feedback from massive stars to remove most of the surrounding gas from the birth cloud that would normally trap LyC radiation (Gotberg et al. 2019). 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 Environment project (fire, Hopkins et al. 2014) using models from the Binary Population and Spectral Synthesis code (bpass, Eldridge et al. 2008;Eldridge & Stanway 2009, 2011Stanway et al. 2016;Eldridge et al. 2017), 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 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 singlestellar 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 §2.2.2) and stars that gain mass via mass exchange in binaries and become massive blue stragglers (see §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 affect on the escape fraction. The outline of the paper is as follows: in §2.1 we briefly describe the cosmological simulations used. In §2.2 we describe the implementation of the various stellar populations included in our model, including massive stars ( §2.2.1), stripped stars ( §2.2.2), blue stragglers ( §2.2.3), and runaway stars ( §2.2.4). In §2.2.5 we discuss other sources of ionizing radiation that we have not included in our simulations, and in §2.2.6 we compare our stellar population synthesis models to the models M16 and R18 used from bpass. In §2.3 we describe our calculation of f esc . In §3 we show our results for a larger mass example halo ( §3.1) and all of the halos in our simulation ( §3.2). In §4 we compare our results with two previous studies, and in §5 we summarize our results.

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 .272, 0.728, 0.045, 0.702, 0.82, 0.96). They first run a dark matter only simulation with a sufficiently large box of size (25 Mpch −1 ) 3 with 256 3 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 zoom-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 5 M , effectively corresponding to 2048 3 over the entire box. The zoomed-in region is dynamically refined according to a pre-set 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 adaptive mesh refinement (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 multi-group radiative transfer (RT) module developed by Rosdahl et al. (2013) to follow radiative transfer of ionizing photons from stars through the ISM and the IGM. The module solves the moment equations for HI, HeI, and HeII 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).

The ionizing emission from a stellar population
In addition to massive stars from a single-stellar population syntheses model, we account for two types of products of binary interaction: stars stripped of their hydrogen-rich envelopes (see §2.2.2) and stars that gain mass both from mass transfer and from coalescence (see §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. 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 §2.1). If the . For comparison, we also include the LyC photon production rate of the model including binary stellar evolution that M16 and R18 used from BPASS (Eldridge et al. 2017;Stanway et al. 2016, 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 −4 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 star-burst99 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%. 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 6 M are shown as a function of time in Fig. 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 −4 for the stripped stars, Z = 4 × 10 −4 for the massive single stars and massive blue stragglers, and Z = 10 −3 for the model from bpass. In our simulation a star particle with Z 2 × 10 −4 would be given the emission rates shown in Fig. 1.
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.
The dominating source of ionizing photons in star-burst99 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 Fig. 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 §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.

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) 1 . 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 −4 .
In Fig. 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 −4 ) 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.

Massive blue stragglers
1 see also the online interface for starburst99 http://www. stsci.edu/science/starburst99/docs/default.html 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.
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 main sequence duration. We therefore use the predictions from starburst99 to estimate the emission rates of ionizing photons from mass gaining stars. It is done as follows. Using lifetimes of massive stars from evolutionary models (see e.g., Fig.  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 Fig. 1, marked with a dash-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 −4 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.

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 −1 . They found that these peculiar velocities fall along a Maxwellian distribution with a dispersion of 24.4 km s −1 . 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 −1 and a minimum space velocity of 28 km s −1 and added them to the initial velocity of the star particle. The directions of motion of runaway stars were chosen at random.
The star particles which 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 a 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 towards 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. 2000Hoogerwerf et al. , 2001Guseinov et al. 2005), we consider this simple approach adequate.

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 2013Chen et al. 2015), and X-ray binaries (Fragos et al. 2013;Schaerer et al. 2019;Senchyna et al. 2019). 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 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 (Eldridge & Stanway 2011;de Mink et al. 2013;Dufton et al. 2011). 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;Ekström et al. 2012;Cantiello et al. 2007). Although many effects of rotation are interesting since they predict an increased emission rate of ionizing photons (e.g., Topping & Shull 2015;Levesque et al. 2012;Kubátová et al. 2019), there is only circumstantial observational evidence for the existence of these processes (see however, Abdul-Masih et al. 2019; and circumstantial evidence from e.g., Martins et al. 2008;Hainich et al. 2015; see also Shenar et al. 2019;Schootemeijer & Langer 2018).

Comparison to bpass
Both M16 and R18 use models from bpass version 2.0, which have very similar initial conditions to 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 possibly 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 −3 . 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 −4 or Z=2 × 10 −4 .

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, ν, iṡ where Ω is the solid angle; N HI (Ω), N HeI (Ω), and N HeII (Ω), are the neutral hydrogen, and 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.Ṅ i (ν) 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Ṅ i (ν) is extrapolated from the age-and metallicity-dependent starburst99 spectrum, with the age-and metallicitydependent photon rates from the stripped star or mas-sive blue straggler spectrum added on for the SSB and MBS runs, respectively.
The escape fraction for each star particle, f esc,i , in each halo snapshot is then computed as the ratio ( 2) For each ramses snapshot of each halo in each run, where f 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 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.

An Example Halo
Before presenting statistical results, we use one halo to illustrate some basic effects on the escape fraction of ionizing photon due to binary evolution in Figures 2  and 3. 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 mass of this halo at z = 7 is 3.16 × 10 10 M and 2 × 10 8 M , respectively.
KC14 found that more massive halos, on average, have lower escape fractions. The photon production rateweighted 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 −1 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 supernovadriven 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 after 10 Myr the most massive single-stars are gone. 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 (dotdashed 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 (dotdashed 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 −1 . We see that the instantaneous ratio of escaped photons changes significantly over time, often reverting back to 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  Figure 2. The middle panel shows the ratio of cumulative, escaped LyC photons for the SSB run to the SB99 run, in green, and that of escaped LyC photons for the MBS run to the SB99, in orange. The bottom panel is analogous to the middle panel but for the instantaneous escape fraction. We have also plotted the star formation rate (SFR) as in Figure 2 to show the offset between when the increase in the number of escaping photons is greatest, and when the SFR is highest. Looking forward, we note that these ratios are larger than the mean values shown in Figure 6. 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.

Statistical Results
The cosmological simulations of KC14 include hundreds of halos ranging in virial mass from 10 8 M to 10 11 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 grey lines show the boxcar smoothed photon-weighted mean values for the three runs. The bottom of the grey 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 box-car 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 photon production rate-weighted mean escape fraction ( fesc ), in percent, as a function of redshift. The top panel shows fesc for all halos in the simulation. The middle panel shows fesc for halos with virial masses less than 10 9 M and the lower panel shows fesc for halos with virial masses greater than 10 9 M . The color-scheme and line-styles are the same as in Figure 2. fesc is greater for both runs that include binary effects than for the SB99 run over all redshifts in every panel. 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 vir < 10 9 M ) 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 as 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 large for the production rate of escaping photons too, but the median values range from roughly 10 45 s −1 at lower virial masses to 10 52 s −1 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 9 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 9 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, and 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 6 M starburst with an age of 100 Myr will be only ∼ 10 47 s −1 . Meanwhile, at that age stripped stars will have a photon production rate roughly 200 times greater, which accounts for the two orders of magnitude difference 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 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 which were not emitting any hydrogen-ionizing photons previously, now emit some radiation that can still be trapped, which lowers the escape fraction.
The box-car smoothed photon-weighted mean escape fractions, shown as the grey line in the top three lefthand 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 49 s −1 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 Ṅ ph f esc is largest for the most massive halos.
Ṅ ph f esc 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 com-pletely 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 9 M at that redshift (middle panel), and halos with a virial mass greater than 10 9 M at that redshift (lower panel). The color-scheme and linetypes are the same as in Figure 2. Overall, the mean escape fractions stay mainly within 10 and 20 percent, 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 is greater at all redshifts for both the SSB and the MBS runs than the SB99 run. Figure 6 compares the total number of photons which 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 9 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 68 photons and 8 × 10 67 photons, respectively, that is by a factor of ∼ 1.125 and ∼ 1.10, respectively. Including both binary products would increase the total number of escaped photons by ∼ 2 × 10 68 , or 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 −4 the additional photons from stripped stars and massive blue stragglers would lead to an increase in the total number of photons by a factor 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 §4 using Equation 5 we 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 9 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 staying 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 Ṅ ph f esc 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 which 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 high enough photon rates to influence the total number of escaped photons, as well as significant and even dominant contributions from stripped stars and massive blue stragglers.
Additionally, KC14 found 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 which will have more time to undergo their own stellar evolution that will generate feedback.
The massive example galaxy shown in Figures 2 and 3 (see §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.

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 10 M ) halo, similar in mass to our example halo in §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 singlestellar population synthesis runs in M16 and R18. One reason for this is that our single-stellar population model includes runaway O/B stars (see §2.2.4), which KC14 showed increased 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 −3 in theirs vs 100 cm −3 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 has 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 in part 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 §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 esc,bin , are more heavily weighted 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 esc,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, [f esc ] R18 = 0.6f esc,ss + 0.4f esc,bin . (4) On the other hand, in this paper only 7% of the ionizing radiation comes from interacting binaries. The overall escape fraction here is, f esc = 0.93f esc,ss + 0.07f esc,bin .
Clearly f esc,bin is more heavily weighted in R18. It is therefore useful to compare f esc,bin , among the three papers. Physically, f esc,bin represents the escape fraction of the additional LyC radiation from interacting binaries. We can calculate f esc,bin using a rearranged generalized version of equations 4 and 5, where x is the fraction of LyC radiation from the singlestellar component in simulations that include binary interaction. From each paper we know x; the escape fraction of the run which includes only the single-stellar population, which we take as f esc,ss ; and the overall escape fraction of the run that includes binary evolution, f esc . This information will give accurate values of f esc,bin for this paper and M16 because both apply a ray tracer in post-processing. However, using Equation 6 to calculate f esc,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 esc,bin for R18 in mind, f esc,bin = 32%, 17%, 57% for M16, R18, and this paper, respectively. Both values of f esc,bin for M16 and R18 are greater than their values of f esc,ss by a factor of ∼ 6, which one would expect given that the radiation from single stars is so well trapped in M16 and R18. Our value increases from f esc,ss to f esc,bin by a slightly smaller factor of ∼4. f esc,bin in this paper is much greater than the values 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. . Photon production rate-weighted mean escape fractions calculated from the ramses cosmological simulation run by KC14 when using a starburst99 spectrum, and when adding in runaway O/B stars, then stripped stars, then massive blue stragglers. The final two steps are what we examine in this paper and they increase the effective escape fraction from 14% to 17%.

CONCLUSIONS
Using ultra-high resolution cosmological radiationhydrodynamic 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 time scales 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 6 M starburst at a metallicity of 2 × 10 −4 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 we see in our SSB and MBS runs, respectively. The physical reason is that the escape fraction and the instantaneous SFR are often strongly anti-correlated with a phase shift of 10 − 200 Myr, reflecting the fact 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 typical massive stars, potentially masquerading as Pop III stars. Figure 7 summarizes the increase of the photonweighted 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, 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;Robertson et al. 2015). While 17% is still a few points lower than the quoted 20%, the most important implication is that our modeling that 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 their lower redshift counterparts.
AS would like to thank Charles Emmett Maher for useful discussions. AS is supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. #DGE1656466. TK 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). YG acknowledges the funding from the Alvin E. Nashman fellowship for Theoretical Astrophysics.