Cradles of the first stars: self-shielding, halo masses, and multiplicity
MMNRAS , 1–13 (2019) Preprint 15 January 2020 Compiled using MNRAS L A TEX style file v3.0
Cradles of the first stars: self-shielding, halo masses, andmultiplicity
Danielle Skinner (cid:63) and John H. Wise Center for Relativistic Astrophysics, School of Physics, Georgia Institute of Technology, Atlanta, GA 30332, USA
15 January 2020
ABSTRACT
The formation of Population III (Pop III) stars is a critical step in the evolutionof the early universe. To understand how these stars affected their metal-enricheddescendants, the details of how, why and where Pop III formation takes place needsto be determined. One of the processes that is assumed to greatly affect the formationof Pop III stars is the presence of a Lyman-Werner (LW) radiation background, thatdestroys H , a necessary coolant in the creation of Pop III stars. Self-shielding canalleviate the effect the LW background has on the H within haloes. In this work,we perform a cosmological simulation to study the birthplaces of Pop III stars, usingthe adaptive mesh refinement code Enzo . We investigate the distribution of host halomasses and its relationship to the LW background intensity. Compared to previouswork, haloes form Pop III stars at much lower masses, up to a factor of a few, dueto the inclusion of H self-shielding. We see no relationship between the LW intensityand host halo mass. Most haloes form multiple Pop III stars, with a median numberof four, up to a maximum of 16, at the instance of Pop III formation. Our resultssuggest that Pop III star formation may be less affected by LW radiation feedbackthan previously thought and that Pop III multiple systems are common. Key words: star formation – cosmology: first stars – stars: Population III – methods:numerical
The first generation of stars transformed a cold and darkuniverse by illuminating and enriching up their neighbor-hoods with radiation and heavy elements. The formation ofthese first stars is a crucial step in the cosmological evolu-tion of the universe because the metals and feedback theydeliver to their local environments are necessary for furtherstar production and chemical enrichment. Without these ini-tial stars, heavier metals would not have been produced,and a different universe would be observed than what isobserved today. These stars are, by definition, metal-free(Population III; Pop III) and are traditionally thought tobe generally massive (Abel et al. 2002; Bromm et al. 2002;Turk et al. 2009; Hosokawa et al. 2011, 2016; Hirano et al.2015), although more recent work has shown that lower mass(
M < M (cid:12) ) Pop III stars can be formed (Greif et al. 2011b;Clark et al. 2011; Stacy et al. 2016). A substantial fractionof these stars will generate prodigious amounts of ionizingphotons and will end their lives in some form of a supernova(e.g. Schaerer 2002; Heger & Woosley 2002a). The supernova (cid:63) E-mail: [email protected] will spread the enriched guts of the Pop III star out acrossits local environment, providing the area with elements theuniverse has not yet seen.Once a halo becomes chemically enriched by the deathof Pop III stars, then by definition, it can no longer formmore Pop III stars. This marks the end of metal-free star for-mation in that pregalactic object. Understanding the mix-ing of metals into the environment of these haloes is nec-essary to constrain the reach of these metals and the ef-fects they have for future star formation. Numerical studieshave shown that turbulence within haloes can mix metalswell down to their resolution limit (Wise & Abel 2008; Greifet al. 2010; Smith et al. 2015). Stars will continue to form inhaloes, but with an increased metal abundance. These starsstill are metal-poor, having metallicities of 10 − to 10 − ofthe solar abundance (Chiaki et al. 2016, 2018; Ritter et al.2016), have a direct chemical connection to Pop III stars,and can survive until the present day (Gnedin & Kravtsov2006; Tumlinson 2010; Griffen et al. 2018; Magg et al. 2018;Ezzeddine et al. 2019). Understanding the formation andchemical abundance of second generation stars can providemore evidence and insight for the earlier generation of PopIII stars. c (cid:13) a r X i v : . [ a s t r o - ph . GA ] J a n D. Skinner and J. H. Wise
Without metals and dust to facilitate efficient radia-tive cooling, Pop III stars rely on H formation in the gas-phase for cooling. These molecules are however fragile todissociation from Lyman-Werner (LW) radiation in the en-ergy range 11.2–13.6 eV, through the Solomon process (Fieldet al. 1966; Stecher & Williams 1967). This is a two-step pro-cess through which H is excited to a higher state, H ∗ , viaabsorption of a LW photon: H + γ → H ∗ (1)This excited state then has a probability of dissociating intotwo hydrogen atoms: H ∗ → H + H (2)Diffuse gas is optically thin to LW radiation, thus abackground builds over time and can suppress H forma-tion, delaying Pop III star formation. Furthermore, nearbysources of LW radiation can boost the intensity above thebackground value which facilitates further H dissociationin minihaloes. This process can be counteracted with a suf-ficient amount of H already present within a halo, via H self-shielding. Halos with an H column density of N H2 ≥ cm − can suppress the photodissociation of H by LWphotons (Draine & Bertoldi 1996), giving the halo a chanceto form Pop III stars. Other sources of the suppression of PopIII star formation are streaming baryonic velocities (Tseli-akhovich et al. 2011; Greif et al. 2011a; Naoz et al. 2012;O’Leary & McQuinn 2012; Schauer et al. 2019; Hirano et al.2017), arising from recombination, and dynamical heating,occurring during the virialization of haloes (Yoshida et al.2003; Fernandez et al. 2014). X-ray heating from high-massX-ray binaries also has the potential to suppress Pop III starformation (Jeon et al. 2014a).Given the paradigm of hierarchical structure formation,haloes grow through smooth accretion and successive merg-ers of smaller haloes. But in which haloes do these firstgenerations of stars form? Their formation rates and loca-tions are important to constrain because they influence thevery beginnings of galaxy formation and cosmic reionization.Various semi-analytic investigations have been conducted tolearn more about the halo collapse criterion, and thus, whichhaloes can host Pop III stars. Tegmark et al. (1997) dis-covered that haloes can have a strongly redshift dependent,minimum baryonic mass of 10 M (cid:12) at z ≈
15. In particu-lar, they derive an analytic expression for the fraction of H needed in a halo for efficient cooling, and determine whichhaloes can cool in a Hubble time. Visbal et al. (2018) deviseda semi-analytic model for the formation of Pop III stars andthe transition to metal-enriched stars. They find that vary-ing the Pop III star formation efficiency, the time from aPop III supernova to metal-enriched star formation, the ex-ternal enrichment of the halo, and the ionizing escape frac-tion, leads to large differences in the star formation historyof Pop III and metal-enriched stars. This method is usefulfor exploring the wide parameter space of star formation inthe early universe, and could lead to further constraints inthe future.Previous work have mainly investigated and establishedthe lower limit for a halo to host Pop III stars, or for a halo tocollapse. However, subsequent simulations have shown thatthey do not necessarily form at this minimum due to theaforementioned physical processes playing a role. Mebane et al. (2018) used a semi-analytic model of early star forma-tion and found that the LW background coming from therapidly increasing supply of Pop III stars becomes respon-sible for suppressing Pop III star formation. Griffen et al.(2018) also found that the LW background can significantlysuppress the amount of potential sites for Pop III star for-mation.Numerical simulations have also been employed to in-vestigate collapse thresholds and Pop III star formation.Machacek et al. (2001, hereafter M01) found similar resultsto those mentioned previously, in that the LW feedback cansuppress the collapse of small mass haloes. They fit a simpleanalytic expression for the mass threshold of haloes given aparticular LW intensity J LW : M TH ( M (cid:12) ) = 1 . × + 8 . × (cid:18) πJ LW − (cid:19) . . (3)We will compare our results with this relation in future sec-tions. Yoshida et al. (2003) also found that cooling is inef-ficient in haloes with a LW background > J , where J is a specific intensity of 10 − erg s − cm − Hz − sr − ,although with sufficient H shielding, haloes are able to coolin the given LW background radiation. In fact, Wise & Abel(2007) found that central shocks drive H formation, allow-ing the halo to cool in a LW background intensity of 1 J .They find that H cooling is always a dominant process, evenin large LW background fluxes. O’Shea & Norman (2008)also found that Pop III star formation can occur in rela-tively high LW backgrounds, implying that the LW back-ground may not be a complete indicator of whether or notPop III stars will form in a given halo.In this paper, we focus on the distribution of host halomasses, not just the minimum, and its dependence on red-shift and the LW background at the instance of star forma-tion, augmenting results from prior work. We aim to providefurther insight into the host haloes of Pop III stars. In § self-shielding are described. In § § § We run and analyze a cosmological simulation with theadaptive mesh refinement (AMR) code
Enzo (Bryan et al.2014) and the toolkit yt (Turk et al. 2011). Enzo uses an N-body adaptive particle-mesh solver (Efstathiou et al. 1985;Couchman 1991) to follow the dark matter (DM) dynam-ics. We use the nine-species (H i , H ii , He i , He ii , He iii ,e − , H , H +2 , H − ) non-equilibrium chemistry model (Abelet al. 1997; Anninos et al. 1997). This simulation is similarto the RP simulation in Wise et al. (2012, hereafter W12),but with updated cosmological and Pop III parameters, andthe inclusion of H shielding.We simulate a 1 Mpc comoving box with a 256 basegrid resolution and a dark matter particle mass of 2001 M (cid:12) .It has a maximum refinement level of 12 which provides amaximal comoving resolution of ∼ MNRAS , 1–13 (2019) radles of the first stars when one of the following criteria are met: (1) the baryonor dark matter overdensity exceeds 3¯ ρ b , DM × (cid:96) (3+ (cid:15) ) , or (2)the local Jeans length is less than eight cell widths. Here ρ b , DM is the baryon or dark matter cosmic mean density, (cid:96) is the AMR level, and (cid:15) = − . Mu-sic (Hahn & Abel 2011) at z = 130 and uses the cosmologi-cal parameters from the Planck collaboration best fit PlanckCollaboration et al. (2014): Ω M = 0.3175, Ω Λ = 0.6825, Ω DM = 0.2685, and h = 0 . z >
20. The sim-ulation is run until z = 9 .
32, when it becomes too compu-tationally expensive to continue. At this point, about 20%of the volume is over 10% ionized. We output 918 datasets,roughly 0.5 Myr apart. In this paper, we focus only on out-puts where Pop III stars have just formed, from z = 27.23to 9.39. Throughout the rest of the paper, all numbers arein physical units, unless otherwise specified.A time-dependent LW optically thin radiation back-ground is applied in the simulation. This was fit in W12(see their Eq. 16) and is consistent with the values in Trenti& Stiavelli (2009). The background evolution of the specificintensity takes the following form:log ( J LW /J ) = A + Bz − Cz + Dz − Ez (4)where ( A, B, C, D, E ) = ( − . , . , − . , . × − , − . × − ). This form only varies by a factor offive in z = 9 −
25 and has a maximum of J LW /J = 0 .
97 ata redshift of 13.5 and slowly declines to 0.60 J at z = 9 . dissociating radiationfield with an optically thin, inverse square profile. TheseLW sources are added on top of the background intensitydescribed above. In this section, we briefly discuss our implementation of PopIII and metal-enriched star formation. We do not considerthe formation and feedback from stellar remnants or asymp-totic giant branch (AGB) stars. For further details, we referthe reader to W12.
The utilized Pop III star formation model is the same asin Wise & Abel (2008) but with updated parameters. Eachstar particle represents a single massive star, and is formedin a cell when the following criteria are met:(i) a metallicity Z ≤ × − Z (cid:12) (ii) a gas number density n > cm − (iii) converging gas flow, ∇ · v gas < H2 > − The critical metallicity marks where dust cooling becomesefficient enough to cause fragmentation at high densities (e.g. Schneider et al. 2006). In practice, we find about 10per cent of Pop III star particles having a non-zero metal-licity below this value. This scenario occurs when a halo isexternally enriched by a nearby supernova. Because of un-certainties with turbulent mixing and collapse timescales,we choose not to predict the final metallicity at zero-agemain sequence and assign the star particle the metallicity ofthe densest cell. The chosen threshold density are consistentwith values found in previous simulations of Pop III starformation (e.g. Hirano et al. 2015) at our resolution limit.Additionally, it is similar to a density that would trigger an-other mesh refinement past level 12. At these scales, Jeanslength refinement is the most common because of the coldtemperatures in the cloud core. The molecular hydrogen cri-terion is also consistent with previous work in the “loitering”phase of molecule formation between 10 and 10 cm − be-fore three-body reactions become important (e.g. Omukaiet al. 2010).If within 10 pc, multiple cells meet this criteria, then asingle Pop III star forms at the center of mass of these cells.In a parameter study, we have found that the number of PopIII star particles is insensitive to this algorithmic mergingdistance. We ran another 1 Mpc simulation with the samesetup as the main simulation with the exception that it was azoom-in simulation focusing on a single minihalo of mass 5 × M (cid:12) collapsing at z (cid:39)
23. Starting from 1 Myr before thecollapse, we followed Pop III star formation in this one halowith a merging distance of 10, 5, 2.5, 1.25, and 0.625 pc. Wefound that even when this number was decreased to 0.625 pc,the number of star particles that formed was unchanged. Wenote that our density threshold for star particle formationis equivalent to a molecular core. Fragmentation of metal-free gas into massive stellar precursors can occur at scalesdown to 1000 AU within molecular cores (e.g. Turk et al.2009; Stacy et al. 2016), therefore our results on multiplicityshould be considered as lower limits.The mass of the Pop III star is randomly sampled froman initial mass function (IMF) of the form: f (log M ) dM = M − . exp (cid:34) − (cid:18) M char M (cid:19) . (cid:35) dM (5)where M char = 20 M (cid:12) (Hirano & Bromm 2017). This char-acteristic mass is different from the W12 choice of 100 M (cid:12) .This results in an exponential cutoff below M char and apower-law IMF above. The Pop III IMF mass ranges from1 ≤ M/M (cid:12) ≤ ph = 29.6 eV, which is appro-priate for the nearly mass-independent surface temperatureof Pop III stars, at 10 K. A Type II supernova results ifthe Pop III star dies with a mass between 11 ≤ M (cid:63) /M (cid:12) ≤ MNRAS000
23. Starting from 1 Myr before thecollapse, we followed Pop III star formation in this one halowith a merging distance of 10, 5, 2.5, 1.25, and 0.625 pc. Wefound that even when this number was decreased to 0.625 pc,the number of star particles that formed was unchanged. Wenote that our density threshold for star particle formationis equivalent to a molecular core. Fragmentation of metal-free gas into massive stellar precursors can occur at scalesdown to 1000 AU within molecular cores (e.g. Turk et al.2009; Stacy et al. 2016), therefore our results on multiplicityshould be considered as lower limits.The mass of the Pop III star is randomly sampled froman initial mass function (IMF) of the form: f (log M ) dM = M − . exp (cid:34) − (cid:18) M char M (cid:19) . (cid:35) dM (5)where M char = 20 M (cid:12) (Hirano & Bromm 2017). This char-acteristic mass is different from the W12 choice of 100 M (cid:12) .This results in an exponential cutoff below M char and apower-law IMF above. The Pop III IMF mass ranges from1 ≤ M/M (cid:12) ≤ ph = 29.6 eV, which is appro-priate for the nearly mass-independent surface temperatureof Pop III stars, at 10 K. A Type II supernova results ifthe Pop III star dies with a mass between 11 ≤ M (cid:63) /M (cid:12) ≤ MNRAS000 , 1–13 (2019)
D. Skinner and J. H. Wise
40. A pair instability supernova (PISNe) results if the stardies with a mass between 140 ≤ M (cid:63) /M (cid:12) ≤ ≤ M (cid:63) /M (cid:12) <
20, the star dies as a normalType II supernova with an energy of 10 erg. From 20 ≤ M (cid:63) /M (cid:12) ≤
40, the star dies as a hypernova with energiestaken from Nomoto et al. (2006). The PISNe have explosionenergies from Heger & Woosley (2002b), of which the ana-lytic function fit to their models can be seen in Eq. 12 ofW12. With our chosen IMF (Eq. 5), 59, 36, and 5 per centof Pop III stars produce Type II SN, black holes with noSN, and PISNe, respectively.
The Pop II star formation model is the same as Wise &Cen (2009), which is equivalent to the above criteria forPop III star formation but without the molecular hydrogenfraction requirement, and with a metallicity greater thanthe above value. Star formation is only allowed for gas withtemperatures
T < min =1000 M (cid:12) for these star particles.The lifetime of these particles is 20 Myr, during whichthey radiate 6000 hydrogen ionising photons per stellarbaryon, or 1.13 × photons s − M − (cid:12) , which is appro-priate for a Salpeter IMF (Schaerer 2003). Their spectra areapproximated with a monochromatic spectrum with an en-ergy of 21.6 eV at a constant luminosity. After living for 4Myr, the stars begin to return supernova energies of 6.8 × erg s − M − (cid:12) back to their surroundings. Self-Shielding
We use the Sobolev-like approximation from Wolcott-Greenet al. (2011) to model H self-shielding. To determine theH shielding factor, the column density of H needs to becalculated: N H2 = n H2 L char , (6)where n H2 is the number density of H and L char is a char-acteristic length over which n H2 is assumed to be constant.The method employed in this simulation is to define thecharacteristic length as: L char = ρ |∇ ρ | . (7)This Sobolev-like method determines the length over whichgas with density ρ is diminished. They found that the mostaccurate, non-ray tracing, method is to use a single Sobolev-like length determined from the mean Sobolev-like lengthover all Cartesian directions.Upon numerically calculating the shielding factor of H for simulated protogalaxies, Wolcott-Green et al. determinedthat a slight adjustment to the shielding factor from Draine& Bertoldi (1996) is sufficient to account for inaccuraciesat higher temperatures. The implemented shielding factor is then (see Eq. 10 in Wolcott-Green et al. (2011)): f sh ( N H2 , T ) = 0 . x/b ) . + 0 . x ) . × exp[ − . × − (1 + x ) . ] (8)Here, x ≡ N H2 / × cm − , b ≡ b/ cm s − , and b is the Doppler broadening parameter. This shielding fac-tor acts as a multiplier to the H dissociation rate, whereif the shielding factor equals one, there is no shielding, andzero if there is maximum shielding. Hartwig et al. (2015)use the same shielding factor equation from Wolcott-Greenet al. (2011) but implement an improved calculation of thecolumn density to incorporate the relative gas velocities andgeometries of the halo. This method takes into account theDoppler-shifting of spectral lines, which can affect H shield-ing. With the improved calculation of the column density,they find that the shielding factor is larger for gas at highertemperatures near 10 K than when the shielding factor iscalculated for an average column density for the halo, mean-ing that self-shielding is generally less efficient for highertemperatures. In this work, we do not take into account theDoppler-shifting of spectral lines, but as will be seen below,our haloes are less massive and contain gas at lower temper-atures. Our calculation of the shielding factor is still validfor the temperature range of our haloes, and therefore, theDoppler-shifting of spectral lines is not expected to affectour results.
To demonstrate how the shielding factor changes as a func-tion of halo mass and redshift, we can examine N H2 betweenthe core and halo virial radius as a function of halo mass foran isothermal halo, with a realistic H halo profile, and useEquation 8 to calculate the shielding factor. The H moleculenumber density within the halo is given by n H2 ( r ) = f H2 ρ R µm H r (9)where µ = 1 .
22 is the mean molecular weight, ρ = 200 ρ c / ρ c = 3H / π G(1 +z) is the critical density. Here, H is the Hubble constant,G is the gravitational constant, and z is the redshift. Thecolumn density is then equation 9 integrated from the radiusof the core to the virial radius, where the core has a radiusof f c R vir : N H2 = f H2 R vir ρ c µm H (cid:18) f c − (cid:19) (10)For H = 70 km s − Mpc − , and assuming the Dopplerbroadening parameter, b , is equal to the circular velocity ofthe halo, Equation 10 and Equation 8 will give the shieldingfactor for a given core radius, f c , redshift, and for a givenH fraction halo profile. Figure 1 shows the shielding factoras a function of halo mass, for redshifts z = 9, 15, and 25 inred, blue and black, respectively. The solid lines indicate acore radius of f c = 10 − and the dashed lines a core radiusof f c = 10 − . To have an accurate depiction of the shield-ing factor for a given halo radius, the H fraction is givenby the H fraction halo profiles at z = 24 . J = 0, which was MNRAS , 1–13 (2019) radles of the first stars Halo Mass [M ]10 S h i e l d i n g F a c t o r z = 9z = 15z = 2510 H2 Column Density [cm ] Figure 1.
Shielding factor as a function of halo mass. Solid linesindicate a core radius of f c = 10 − and the dashed lines indi-cate a core radius of f c = 10 − . For these parameters, assumingthe Doppler broadening parameter b equals the circular velocity,and assuming H fraction halo profiles from O’Shea & Norman(2007), the shielding factor is below 1 percent. The green lineshows the shielding factor as a function of H column density for b = 10 km s − . Note that a shielding factor of 1 (0) correspondsto no (complete) shielding. extracted directly from their plot. These profiles were ap-proximated with a piecewise power-law for simplicity givenby log ( f H2 ) = − . × log ( r ) − .
039 for r < . ( f H2 ) = − . × log ( r ) − .
737 for r > . is being added to the calculation of the column density.The more H present within a halo, the smaller the shield-ing factor will become. Across the top axis and plotted asa green line, the shielding factor as a function of H col-umn density is shown for a Doppler broadening parameter b = 10 km s − . At each redshift, all halo masses, and foreach core radius, we see the shielding factor is significantlysmaller than one, implying that H shielding for these haloesis almost at a maximum. Self-shielding will then significantlydecrease the dissociation rate of H in the halo core, allowingfor the presence of a significant H fraction even in a strongLW background. In this section, we will present the distribution of host halomasses of Pop III stars, the relation between the mass dis-tribution and the LW background radiation, and the distri-bution of various Pop III properties.To ensure that our box is representative of a typicalpiece of the universe, Figure 2 shows the simulated halomass function along with the analytic Sheth-Tormen massfunction at z = 9 .
3. The analytic function closely matchesour simulation until the mass resolution of our simulation isunable to resolve below 10 M (cid:12) haloes that contain ≈ Halo Mass [M ]10 n ( > M H a l o ) [ M p c ] Sheth-TormenSimulated Halos
Figure 2.
Halo mass function of the last output of the simulationat z = 9.3. The thick, black line is the analytic Sheth-Tormenmass function. The distribution of halo masses matches well withanalytical expectations above 10 M (cid:12) , where haloes are capturedby 50 particles. The vertical blue line indicates our halo massresolution limit of 10 M (cid:12) . particles. After a sufficient amount of time after the begin-ning of the simulation, enough matter has collapsed to be-gin forming Pop III stars throughout the box. When the firstPop III star forms in our simulation at z = 27, there are only6 haloes that are above 10 M (cid:12) , which are resolved by 50DM particles, consistent with analytical expectations froman ellipsoidal variant of Press-Schechter formalism (Press &Schechter 1974; Sheth et al. 2001). These Pop III stars eitherbecome black holes or supernovae at the end of their lives,with the latter resulting in an expulsion of metals into thestar’s surroundings. Metal-enriched star formation begins totake place as soon as enough metals are present within ahalo. In the meantime, haloes form and merge with one an-other, resulting in a wide variety of galaxies, and reionizationstarts to take place. By the end of the simulation, 20% ofthe box has an ionized fraction greater than 10 per cent, ascan be seen in Figure 3, which also shows the mass weightedLW intensity coming from stars within the simulation as afunction of time, the LW background given by Equation 4in units of J , and temperature projections of the box atcertain points in time.In our simulation, 688 Pop III stars form between 27.23 ≥ z ≥ . × − M (cid:12) yr − Mpc − at z = 21 . z ≤ − M (cid:12) yr − Mpc − at z = 15and z = 20 respectively. Jaacks et al. (2019) find a SFRD MNRAS000
Halo mass function of the last output of the simulationat z = 9.3. The thick, black line is the analytic Sheth-Tormenmass function. The distribution of halo masses matches well withanalytical expectations above 10 M (cid:12) , where haloes are capturedby 50 particles. The vertical blue line indicates our halo massresolution limit of 10 M (cid:12) . particles. After a sufficient amount of time after the begin-ning of the simulation, enough matter has collapsed to be-gin forming Pop III stars throughout the box. When the firstPop III star forms in our simulation at z = 27, there are only6 haloes that are above 10 M (cid:12) , which are resolved by 50DM particles, consistent with analytical expectations froman ellipsoidal variant of Press-Schechter formalism (Press &Schechter 1974; Sheth et al. 2001). These Pop III stars eitherbecome black holes or supernovae at the end of their lives,with the latter resulting in an expulsion of metals into thestar’s surroundings. Metal-enriched star formation begins totake place as soon as enough metals are present within ahalo. In the meantime, haloes form and merge with one an-other, resulting in a wide variety of galaxies, and reionizationstarts to take place. By the end of the simulation, 20% ofthe box has an ionized fraction greater than 10 per cent, ascan be seen in Figure 3, which also shows the mass weightedLW intensity coming from stars within the simulation as afunction of time, the LW background given by Equation 4in units of J , and temperature projections of the box atcertain points in time.In our simulation, 688 Pop III stars form between 27.23 ≥ z ≥ . × − M (cid:12) yr − Mpc − at z = 21 . z ≤ − M (cid:12) yr − Mpc − at z = 15and z = 20 respectively. Jaacks et al. (2019) find a SFRD MNRAS000 , 1–13 (2019)
D. Skinner and J. H. Wise
Figure 3.
For the entire box, the mass-weighted average LW in-tensity from stars within the simulation is shown in solid red,and the fraction of the volume ionized above 10 per cent ver-sus time is shown in solid black. The red dashed line representsthe LW background implemented in the simulation (see Equation4). The solid red line in the legend refers to both the solid anddashed lines in the figure. Temperature projections are shown atthe points indicated by the arrows. of 10 . M (cid:12) yr − Mpc − at z = 20, which flattens out to10 − M (cid:12) yr − Mpc − by z = 7. Our SFRD rates are con-sistent with these other studies, but it is likely to changedepending on the particular modes captured by the initialconditions. To determine the host haloes of new Pop III stars, we use thehalo finding code,
Rockstar (Behroozi et al. 2013), an algo-rithm that adaptively refines haloes found with a friends-of-friends method in position and velocity phase-space and alsotemporally. The mass determined by
Rockstar is the darkmatter mass, so to determine the total mass of the haloes, wemultiply the masses by Ω M / Ω DM . The total mass is subse-quently used in our analysis. For stars that form in subhaloesidentified by Rockstar , we assign the parent halo mass asthe host halo mass. There is one interesting case that we ex-clude from our analysis, where the parent halo is above theatomic cooling limit, because the radiative cooling physicsare different in this regime and is not sensitive to a LW back-ground. Here seven Pop III stars form in a pristine 10 M (cid:12) subhalo inside of a 3 × M (cid:12) atomic cooling parent halo.Because this object only contains 1 per cent of the numberof Pop III stars, it only affects our average halo masses bya similar amount. We also take the center and virial radiusfrom the Rockstar halo catalog. About 1% of Pop III starsdo not form in a halo identified by
Rockstar . Assumingthat Pop III stars must form in collapsed haloes, we calcu-late the virial mass and radius of a halo centered on the Pop
100 200 300 400 500Time [Myr]10 S F R D [ M y r M p c ]
20 15 12 10Redshift
Figure 4.
The comoving star formation rate density of PopIII stars. The SFRD peaks at a redshift of 20, and slowly fallsoff. Note that metal-enriched star formation continually increasesthroughout time, although that is not plotted here. Due to thesmall box size of our simulation, this SFRD is particular for ourbox.
III star directly from the simulation data by determiningthe sphere that contains an average overdensity of 200 ρ c .We calculate the LW intensity impinging on each halo bysumming up the contribution coming from each radiatingstar particle outside the virial radius of the host halo. Thisis then added to the constant LW background implementedin the simulation (Eq. 4).In redshift bins of ∆ z = 1, the mean halo mass host-ing Pop III stars is determined, and plotted as black dotsagainst the redshift bins in Figure 5. The median for eachbin is represented by a cross, and the 15.9 and 84.1 per-centiles are plotted below and above the mean, respectively.The M01 relation (Equation 3) for our LW background in-tensity in Eq. 4 is shown as the black solid line. Notably,the mean halo mass falls well below the relation, at < M (cid:12) until z = 15 (270 Myr), at which point, the mean halomass steadily rises above 3 × M (cid:12) . The large discrepancybetween our data and the mass threshold from M01 is indica-tive of the H shielding which is included in our simulation.H shielding allows for haloes to cool at much lower masses,and therefore, Pop III stars are forming in these haloes atearlier times. In the M01 analysis, H shielding is neglectedin their calculations for a variety of reasons, including theDoppler broadening of LW bands and large column densitiesof H only beginning to form at late times. Interestingly, afull radiation hydrodynamics simulation, such as ours, doesnot produce identical results as M01 predicts.Throughout the literature, there appears to be a con-sensus that H self-shielding will help smaller mass haloescollapse in the presence of a LW background radiation (e.g.Yoshida et al. 2003; Ricotti et al. 2001; Glover & Brand2001; Hartwig et al. 2015). As can be seen by the red solidline in Figure 5, nearly 100% of the haloes hosting Pop IIIstars fall below the M01 relation, until around a redshift of15, at which point, the mean halo mass begins to rise abovethe relation. The increase in the mean halo mass above the MNRAS , 1–13 (2019) radles of the first stars P e r c e n t a g e o f H a l o s B e l o w M a c h a c e k +
20 15 12 10Redshift200 300 400 500Time [Myr]5.505.756.006.256.506.757.007.25 L o g ( M H a l o / M ) Figure 5.
For haloes hosting new Pop III stars, the mean halomass for redshift bins of ∆ z = 1 is plotted along the left hand sideversus time. The black scatter points indicate the mean halo massfor that redshift bin and the x indicate the median halo mass. Theerror bars indicate the 15.9 and 84.1 percentiles. The right handside shows the percentage of haloes that fall below the M01 massthreshold as a function of redshift. The M01 mass threshold (Eq.3), is plotted for the constant LW background in Eq. 4. Almostall haloes fall below the M01 relation, until ≈
380 Myr, when themean halo mass rises above the relation.
M01 relation is indicative of the increased LW intensity per-meating the box due to a galaxy that becomes very activein star formation at late times. This galaxy has a halo massof 3 . × M (cid:12) , a total stellar mass of 3 . × M (cid:12) , and hasa peak Pop III SFR of 2 . × − M (cid:12) yr − at z = 15 (250Myr) and a peak metal enriched SFR of 5 . × − M (cid:12) yr − at z = 9 . J . The LW background intensity at the instance of Pop IIIstar formation versus the host halo mass for each Pop IIIhalo in our dataset is shown in Figure 6. Because the LWbackground only varies by a factor of five as Pop III starsform, most points lie within a small intensity range, justbelow 1 J . Each point is colored by the redshift where thePop III star forms and the mass threshold from M01 (Eq.3) for the given LW background (Eq. 4) is the solid blackline. We see again that almost all haloes form Pop III starsbelow the M01 threshold. There are also a few haloes thatform Pop III stars in very high J LW , below the relation.This situation arises when there are multiple Pop III starsforming at about the same time, within ≈
100 pc of eachother. Star formation can occur in a neighboring halo of aPop III star whose LW radiation does not have ample timeto photodissociate enough H to completely suppress star Halo Mass [M ]10 J L W / J Machacek et al. 2001Local Sources + BG (W12) 101214161820222426 R e d s h i f t Figure 6.
The average LW background for host haloes of new PopIII stars is plotted versus the host halo mass, colored by redshift.The Machacek et al. relation is plotted given the background LWin Eq. 4. Almost all haloes fall below the relation, across a rangeof redshifts. There are extreme cases where Pop III stars form insmall haloes located in a very intense LW radiation field or largehaloes in relatively weak radiation fields. formation. The high J LW is generally an indicator that thehaloes will likely have their star formation suppressed, butin cases where Pop III star formation occurs at about thesame time, star formation will not be suppressed, so longas they are a suitable distance away from each other. Asthere are only a few data points in this region, this situationis rare. It should be noted that there are duplicate haloeswithin this plot, since haloes are allowed to form multiplePop III stars if the conditions are sufficient and each pointrepresents an instance of Pop III formation. The groupedpoints in the higher end of J LW are representative of suchhaloes.We find that a total of 84% of the haloes forming Pop IIIstars lie below the M01 mass threshold over the entire sim-ulation redshift range. Given that the background intensityonly varies by a factor of five, previous work suggested thatmost haloes should have formed stars around or above theM01 relationship. However, we see a large spread in host halomasses, over an order of magnitude, and no clear relation-ship between the instantaneous LW intensity and the hosthalo mass. This behavior suggests that the LW intensity isnot the primary indicator of the mass at which a metal-freehalo collapses. H self-shielding is the main cause behindthis conclusion because the internal dynamics and chemo-thermal evolution are now isolated from its large-scale envi-ronment and mostly depends on the local conditions withinthe halo. Alternatively, it is possible that there could be acorrelation between the time-integrated LW background andthe host halo mass. These calculations are outside the scopeof this paper and require further study. We now inspect the number of Pop III stars and the totalmass of Pop III stars per halo. We are only tracking star-
MNRAS000
MNRAS000 , 1–13 (2019)
D. Skinner and J. H. Wise forming regions where the star particles form at densitiesequivalent to molecular cores that may form multiple mas-sive stars. Furthermore, we are not modeling the formationof low-mass stars ( < M (cid:12) ) that may form out of fragmentedgas at higher densities and smaller scales. Thus our resultsin this section should be considered as lower limits to thenumber of Pop III stars.The left panel of Figure 7 shows the number of Pop IIIstars per halo for a given halo mass, just after star formation.The histogram of the number of Pop III stars for all massesis projected on the right hand side. We find that a mediannumber of four Pop III stars form per halo, with a maximumof 16 Pop III stars forming per halo. Since we did not restrictthe number of Pop III stars that can form in a halo, we findthat the conditions are often sufficient for multiple Pop IIIstars to form in a single halo. Out of the haloes formingPop III stars, only 16% of them form single Pop III stars,whereas 54% form between two and five Pop III stars, similarto the radiation hydrodynamic simulations of Susa (2013);Susa et al. (2014). It should be noted that these haloes arenot necessarily forming all of the Pop III stars at once. Forexample, a halo can form stars again if it did not host asupernova. The number of Pop III stars that form within ahalo is indicative of the star formation history of that halo.Figure 8 shows the total mass of Pop III stars per halo fora given halo mass, where the histogram of the total PopIII mass for all haloes is projected on the right hand side.Lines of constant star formation efficiency are overplotted.We find that most of our Pop III stars are forming between10 − < f (cid:63) < − . The mean total mass of Pop III stars forall haloes is 195 M (cid:12) . Note that the mass of the Pop III starsdepends on the chosen IMF.At the end of the simulation, the right panel of Figure7 shows the distribution of the number of Pop III stars perhalo for a given halo mass, for all haloes in the simulationthat host either a living or a remnant Pop III star. The redline shows the median number of Pop III stars for each halomass bin. At low masses, the number of Pop III stars per halosits around 2 Pop III stars, and quickly rises once the halomass reaches 10 M (cid:12) . Xu et al. (2013) found that at z = 15,haloes with masses at about 10 M (cid:12) will contain between 10and 100 Pop III stellar remnants, similar to our results. Wecan also look at the spread in creation times of these PopIII stars per halo, shown in Figure 9. The haloes that have azero spread in their Pop III creation time, meaning the PopIII stars formed at the same time, lie at 10 − Myr. The redline indicates the median spread in creation times for eachhalo mass bin. There are two groups within this plot. Thefirst group contains haloes that form their Pop III stars ina spread of <
10 Myr, with a median spread of about 0.1Myr. The second group contains haloes that form their PopIII stars in a spread of >
10 Myr, with a median spreadof about 100 Myr. The second group represents typicallylarger mass haloes that have acquired older Pop III stars bythe end of the simulation through their long merger history.These haloes have acquired these Pop III stars from outsidehaloes, accounting for the large spread in Pop III creationtime as well as the increased number of Pop III stars, as seenin the right panel of Figure 7.
Previous studies have focused on the minimum halo mass ofPop III hosts as a function of LW intensity. However, theyhave also found that Pop III stars often form in haloes up toan order of magnitude greater than the minimum. We alsohave found that Pop III star formation occurs in a similarrange of host halo masses throughout our simulation.There are three main causes of this significant variationin mass. First given our assumed IMF (Equation 5), 36 percent of metal-free stars directly form a black hole withoutmetal ejecta (Heger et al. 2003), leaving the halo chemicallypristine. Second-generation stars forming later in that halocould still be metal-free if no external enrichment occurs.Thus they would form in haloes substantially larger thanthe minimum, as it takes tens of millions of years for thehalo to recover its gas after radiative feedback (Muratovet al. 2013; Jeon et al. 2014b).Second, dynamical heating from mergers and accretionprovide a heating source inside the growing halo, prevent-ing the gas from efficiently cooling via H . This turbulentstirring can delay star formation that could cool throughH in an ideal situation but otherwise forms in haloes moremassive than the minimum (Yoshida et al. 2003; Wise et al.2019).Finally, temporal fluctuations in the local LW radia-tion field can greatly affect the amount of H within a haloand thus, how efficiently the halo can cool. While this mayonly be a small effect in some haloes where the local LWradiation fluctuations are relatively small and star forma-tion is distant, Pop III star formation may be significantlydelayed or completely prevented in haloes that are close toactive star formation sites. In summary, the timing of PopIII star formation depends on both local – halo histories ofstar formation and growth – and environmental properties.The combination of these three processes determines whichhaloes will be able to form Pop III stars, and thus results ina wide range of host halo masses. There has been a large amount of work done throughoutthe literature investigating the relationship between the LWbackground intensity and halo masses. Tegmark et al. (1997)numerically integrated a chemical network to calculate theminimum mass a halo must have in order to cool, by con-sidering H formation and cooling rates. They conclude thatcooling by H is efficient and leads to the first formation ofstructures. They found a minimum halo mass needed to col-lapse which depends on the virial temperature of the haloand the virial redshift. At z = 15, this minimum halo massis 10 M (cid:12) . Trenti et al. (2009) uses a similar argument asTegmark et al. (1997), and found that metal-free haloes canexist until z ≈ from Trenti & Stiavelli(2009) as one part of their halo mass model (see the blue MNRAS , 1–13 (2019) radles of the first stars Halo Mass [M ]246810121416 N P o p III p e r H a l o At Pop III Formation 0 100 Halo Mass [M ]0510152025 N P o p III p e r H a l o Final Redshift 024681012
Figure 7.
The number of Pop III stars per halo versus halo mass at the instance of Pop III star formation (left) and the final redshift(right). At the instance of formation, Pop III stars will form in a halo which contains a median number of four Pop III stars, with somecontaining as many as 16 Pop III stars. On the right, the red line indicates the median number of Pop III stars in each halo mass bin.By the end of the simulation, high mass haloes contain a large number of Pop III stars due to their long merger history with smallerhaloes hosting Pop III stars. Halo Mass [M ]10 M P o p III p e r H a l o [ M ] f = f = f = Figure 8.
Total mass of Pop III stars in haloes hosting newPop III stars versus halo mass. Lines of constant star formationefficiencies are overplotted. Most haloes form Pop III stars atefficiencies between 10 − and 10 − . solid line in Figure 10, given our LW background). NeitherTrenti et al. (2009) nor Tegmark et al. (1997) included theeffects of H self-shielding, which can explain the discrep-ancy between our data and their mass threshold relationship(Figure 10). For a LW background of J LW /J = 0 .
2, we seehaloes hosting Pop III stars down to 4 × M (cid:12) , while therelation from Trenti & Stiavelli (2009) would expect to seeminimum halo masses capable of cooling, and thus capableof hosting Pop III stars, at 10 M (cid:12) , almost two orders ofmagnitude larger. Since we include H self-shielding, PopIII stars are allowed to form in smaller mass haloes for agiven LW background. Halo Mass [M ]10 Sp r e a d i n P o p III C r e a t i o n T i m e [ M y r ] Figure 9.
The spread in Pop III creation times per halo versushalo mass at the end of the simulation. The red line indicates themedian spread in creation times for each halo mass bin. Thereare two groups that appear in this plot, haloes that have a spreadof less than 10 Myr, representing haloes that form Pop III starsquickly, and haloes that have a spread of greater than 10 Myr,representing larger mass haloes that acquire old Pop III starsthrough mergers.
Visbal et al. (2014) used one-zone models to study thegaseous centers of dark matter haloes and how they reactto a LW background. They use a particular parametrizationof the central gas density, which then follows their analyticequation for critical mass (see their Equation 4) to withina factor of two. This critical mass is plotted in our Figure10 at redshift 20. Here almost all of our haloes fall belowtheir critical mass relationship, with only the highest masshaloes rising above. Visbal et al. also compare their results to
MNRAS000
MNRAS000 , 1–13 (2019) D. Skinner and J. H. Wise three-dimensional hydrodynamical simulations, where theyfind similar threshold halo masses as Machacek et al. (2001).Mebane et al. (2018) used a semi-analytic model of star for-mation, including feedback properties such as a LW back-ground, photoionization due to Pop III stars, supernovae ofPop III stars and metal-enriched stars, and chemical enrich-ment, to determine for how long Pop III stars will survivefor and found that Pop III stars can continue to form untilz ≈
6. They found a minimum halo mass for hosting PopIII star formation, but with the caveat that they did not in-clude H self-shielding. They found that when Pop III starscontribute to the LW background, the minimum halo massfor Pop III star formation is 4 × M (cid:12) at z = 20, simi-lar to results from Tegmark et al. (1997). The results fromTrenti et al. (2009) and Mebane et al. (2018) have higherhalo masses in comparison with M01, although this moreclosely matches results from simulations (see Wise & Abel2007; O’Shea & Norman 2008).It should be noted that because our box is fairly small,we cannot capture the cosmic variance of rare haloes andgalaxies. For example, at late times, a single galaxy dom-inates the LW radiation field, and drives up the host halomasses. However, we would expect this to happen at dif-ferent times in other cosmological volumes. Therefore, wecannot directly compare the time dependence, however, acomparison as a function of LW intensity is still valid.Yoshida et al. (2003) used cosmological simulations tostudy the formation of primordial star-forming clouds. Theyfollowed the growth of structure to find where gas cools andcondenses which would form the first stars, and how theeffects of LW radiation may affect these gas clouds. Impor-tantly, they included H self-shielding within haloes. In aseries of simulations, they found a minimum halo mass forthose haloes hosting gas clouds which may result in activestar formation (see their Figure 12). They found that in thepresence of a LW background of 0.01 J , the minimum halomass is 7 × M (cid:12) . When H self-shielding is taken intoaccount along with a LW background, their minimum halomass decreases to 4 . × M (cid:12) . This value lies close to thecase where there is no LW background applied, where theminimum halo mass is 3 . × M (cid:12) . They found that H self-shielding does appear to be an efficient mechanism en-abling primordial gas cooling. In comparison with our work,we find similar minimum halo masses for each redshift, al-though we do see a wider range of halo masses, ranging from2 . × M (cid:12) to 1 . × M (cid:12) . While our haloes never expe-rience a LW intensity as low as 0.01 J , we find that theirminimum mass for this LW intensity lies in the same massrange as our haloes (see blue open circle in Figure 10).Wise & Abel (2007) used cosmological simulations toinvestigate H cooling in a LW background. They foundthat H cooling is dominant even when there is a large LWbackground present. They did not include self-shielding andsubsequently found halo masses at their collapse that liewell above the M01 relation (see red open triangles in Fig-ure 10). Their J LW = 0 control is plotted in Figure 10 at J LW /J = 10 − . O’Shea & Norman (2008) used cosmologi-cal simulations to investigate Pop III star formation in var-ious LW backgrounds. They found that due to an increasedLW background, there is a delay in star formation, and thusthere is an increase in the halo masses at collapse. Theyalso ignored H self-shielding in their calculations which can J LW /J H a l o M a ss [ M ] T r e n t i + V i s b a l + M a c h a c e k + This workYoshida+ 2003O'Shea+ 2008Wise+ 2007 v Naoz =1 vbc v Naoz =1.7 vbc v Naoz =3.4 vbc v Sch. =0 vbc v Sch. =1.1 vbc v Sch. =2.2 vbc v Sch. =3.3 vbc
Figure 10.
The halo mass versus the average LW intensity nor-malized by J21 versus for this work and a variety of other works.The black solid line is from M01 given our LW background, bluesolid line from Trenti & Stiavelli (2009) at z = 20, the purplesolid line from Visbal et al. (2014) at z = 20, blue open circlefrom Yoshida et al. (2003), green open squares from O’Shea &Norman (2008), and red open triangles from Wise & Abel (2007). J LW = 0 for the various sources is plotted here at J LW / J =10 − . The right panel represents halo masses corresponding todifferent streaming velocities from Naoz et al. (2013) and Schaueret al. (2019). The solid bands of halo masses correspond to thecharacteristic masses of haloes to contain a baryon fraction thatis half of the cosmic mean for various redshifts and for the shownstreaming velocities from Naoz et al. (2013). For each band, themaximum halo mass occurs at z = 15. For streaming velocitiesof 3 . σ vbc , 1 . σ vbc , and 1 σ vbc (5.8 km s − at z = 199), the min-imum halo mass occurs at z = 17 , , and 27, respectively. Theminimum halo mass hosting cold dense gas for different streamingvelocities, given at their starting redshift z = 200, from Schaueret al. (2019) are shown as horizontal black lines. again account for the increased halo masses they found com-pared to this work, as can be seen in Figure 10 (open greensquares). Their halo masses are similar to the results of Wise& Abel (2007). Their control of J = 0 control is plottedat J LW /J = 10 − . Hirano et al. (2015) used cosmologicalsimulations to study 1540 collapsing metal-free gas cloudsin the early universe to derive the corresponding primor-dial stellar mass distribution. They find that stars form inhaloes above virial masses M vir = 2 . × M (cid:12) at z = 30and M vir = 9 . × M (cid:12) at z = 10. These masses are consis-tent with our results with our mean halo mass lying betweenthese values. There are a few minor shortcomings to this work that shouldbe noted. As discussed in § MNRAS , 1–13 (2019) radles of the first stars results are only strengthened by the possibility of havingmore Pop III stars forming at these locations.Schauer et al. (2017) investigated the escape fractionof LW photons in the near and far-field and found that theLW escape fraction of atomic cooling haloes can vary signif-icantly depending on the ionisation front within the halo. Inthe near-field with an ionisation front breaking out of a halo,they find that the LW escape fraction is greater than 95%.But when the ionisation front is contained within the halo,they found that the LW escape fraction can range from 3%to 88%. In general, they find that the LW escape fractionsin the far-field are higher than 75%. Their work was basedon previous work, done by Schauer et al. (2015), where theystudy LW radiation coming from single Pop III stars, ratherthan stellar populations. They find similar results as Schaueret al. (2017), with the exception that H shielding by neu-tral hydrogen decreases LW escape fractions as compared towhen only H shielding is considered. In their more recentwork, they find that this difference in LW escape fractionsbetween when H shielding by neutral hydrogen is includedor not is negligible.Kitayama et al. (2004) investigated the transition be-tween D-type and R-type ionization fronts within haloesto determine the escape fractions of ionizing and photodis-sociating photons using a one-dimensional Lagrangian hy-drodynamics code. They found that in haloes with masses < M (cid:12) , escape fractions are greater than 80 %, and inlarger mass haloes, with masses > M (cid:12) , the escape frac-tions are essentially zero. Since we are not accounting forthe reduced escape fraction within our haloes, we are over-estimating the LW radiation coming from point sources inour simulation. We could account for this by decreasing theintensity of the stars, since the overall LW radiation wouldbe reduced, although this should not significantly affect ourresults as we find little dependence on LW intensity.We also do not include streaming velocities betweenbaryons and dark matter within our simulation (Tseli-akhovich et al. 2011; Greif et al. 2011a; Naoz et al. 2012;O’Leary & McQuinn 2012; Schauer et al. 2019; Hirano et al.2017). Streaming velocities can suppress star formation atvery high ( z > ∼
20) redshifts only, because the streaming ve-locities decrease inversely with the scale factor. Pop III starformation is delayed by an average of δz = 4 by increas-ing the halo mass needed to overcome the bulk velocity byabout a factor of three (Greif et al. 2011a). The streamingvelocities will generally result in lower gas densities withinhaloes as well as an offset of the peak density from the centerof the halo (O’Leary & McQuinn 2012). Naoz et al. (2013)estimated the minimum mass needed to retain the bulk ofthe baryons within the dark matter halo using a series ofcosmological simulations, the results of which can be seen inFigure 10 (see shaded bands). These haloes only reach up toabout 2 − × M (cid:12) when the streaming velocity is about1 σ vbc = 3 km s − at z = 99, which is smaller than the starforming haloes in our simulation. Naoz et al. (2013) mentionthat at the largest streaming velocity, their simulations donot match expected values due to a small sample size of highmass haloes.A more recent study by Schauer et al. (2019), who usedcosmological simulations to study cold gas content withinhaloes for varying streaming velocities, showed that highstreaming velocities could restrict first star formation to atomic cooling haloes in the most extreme cases. Their min-imum halo masses which host cold dense gas for variousstreaming velocities are shown as black horizontal lines inFigure 10. Hirano et al. (2017) investigated the effect thatstreaming velocities have on massive black hole seed forma-tion and find that in regions with high streaming velocitiesand for z >
30, first generation star formation can be sup-pressed in minihaloes up to the atomic cooling limit.The results of Naoz et al. (2013) and (Schauer et al.2019) in Figure 10 show that typical streaming velocities( ≤ σ vbc = 5 . − at z = 199) would affect haloesbelow 10 M (cid:12) , especially at redshifts z >
20 when stream-ing velocities are the strongest. Our results sample mini-haloes that experience such streaming velocities primarilyat z <
20, whereas Hirano et al. (2017) studied much ear-lier structure formation. In the evolution of a halo affectedby streaming velocities, it would have started to accumulategas once above these critical masses while also being irra-diated by a LW background. Afterwards, the collapse willbe controlled by the processes discussed in Section 4.1, andthus streaming velocities would likely not affect our mainconclusion that H plays an important role in Pop III starformation.Finally with most studies of metal-free stars, there isthe caveat of its uncertain IMF that affects the individualstellar luminosities and surface temperatures of the Pop IIIstars. Changes in the IMF could alter the multiplicity andtotal stellar mass of the system because the radiative feed-back from the earlier stars can affect the thermodynamics ofother molecular clouds within the same halo. Ultimately theIMF determines the fraction of stars going supernova andthe metal production in pre-galactic objects that is the rootcause of suppressing and ending Pop III star formation inthe early universe. In this work, we presented the analysis of the host halo massdistribution of Pop III stars and how it relates to the LWbackground radiation. From our results, we conclude the fol-lowing: • We find that Pop III stars are forming in haloes with amean mass below 10 M (cid:12) , which falls well below the M01relation, until z = 15, at which point the mean halo massrises above the M01 relation to a mean value of 3 × M (cid:12) due to the domination of LW radiation produced by metal-enriched stars in the most massive galaxy in the simulationbox. • H self-shielding allows haloes below the M01 mass rela-tion, as small as half as previously thought, to cool and formPop III stars. Self-shielding creates a disconnect between theinstantaneous LW background and the collapse, resulting ina broad mass distribution of metal-free star forming haloesand no strict correlation between the two quantities. • Halos are likely to form multiple Pop III stars. At theinstance of Pop III star formation, a halo will have formeda median number of four Pop III stars, up to a maximum of16. • By the end of the simulation, haloes with masses > M (cid:12) acquire multiple Pop III stars due to mergers with MNRAS000
20, whereas Hirano et al. (2017) studied much ear-lier structure formation. In the evolution of a halo affectedby streaming velocities, it would have started to accumulategas once above these critical masses while also being irra-diated by a LW background. Afterwards, the collapse willbe controlled by the processes discussed in Section 4.1, andthus streaming velocities would likely not affect our mainconclusion that H plays an important role in Pop III starformation.Finally with most studies of metal-free stars, there isthe caveat of its uncertain IMF that affects the individualstellar luminosities and surface temperatures of the Pop IIIstars. Changes in the IMF could alter the multiplicity andtotal stellar mass of the system because the radiative feed-back from the earlier stars can affect the thermodynamics ofother molecular clouds within the same halo. Ultimately theIMF determines the fraction of stars going supernova andthe metal production in pre-galactic objects that is the rootcause of suppressing and ending Pop III star formation inthe early universe. In this work, we presented the analysis of the host halo massdistribution of Pop III stars and how it relates to the LWbackground radiation. From our results, we conclude the fol-lowing: • We find that Pop III stars are forming in haloes with amean mass below 10 M (cid:12) , which falls well below the M01relation, until z = 15, at which point the mean halo massrises above the M01 relation to a mean value of 3 × M (cid:12) due to the domination of LW radiation produced by metal-enriched stars in the most massive galaxy in the simulationbox. • H self-shielding allows haloes below the M01 mass rela-tion, as small as half as previously thought, to cool and formPop III stars. Self-shielding creates a disconnect between theinstantaneous LW background and the collapse, resulting ina broad mass distribution of metal-free star forming haloesand no strict correlation between the two quantities. • Halos are likely to form multiple Pop III stars. At theinstance of Pop III star formation, a halo will have formeda median number of four Pop III stars, up to a maximum of16. • By the end of the simulation, haloes with masses > M (cid:12) acquire multiple Pop III stars due to mergers with MNRAS000 , 1–13 (2019) D. Skinner and J. H. Wise smaller mass haloes, resulting in galaxies with young andold metal-free stars.Our results provide another piece of the Pop III puzzle.With the inclusion of H self-shielding, we find Pop III starsforming in smaller mass haloes than previously predicted.These smaller mass haloes provide earlier and additionalsites for Pop III star formation, and therefore, more metal-enrichment within haloes. We also find that Pop III starsdo not necessarily have to form in isolation, in fact, theyare more likely to form in a halo that has already formedmultiple other Pop III stars. Since these stars are gener-ally massive and are likely to produce stellar mass blackholes, their remnants may be gravitational wave candidatesdetected through LIGO (Hartwig et al. 2016). The exactdetection signature of these gravitational waves is a topicwhich needs to be studied further. While we do not see anyrelationship between the LW intensity and Pop III host halomass, in order to determine exactly how the LW intensityinfluences the host halo mass, a more systematic approachneeds to be taken. The time-integrated LW background lead-ing up to the formation of a Pop III star would be an in-teresting calculation to follow up with, and may result in acorrelation with the host halo mass.The assumptions often made about the formation ofPop III stars, like how they will only form in isolation andhow the LW background mainly determines the collapsemass, appear to break down when self-shielding is taken intoaccount and restrictions about their formation are lifted.Further work needs to be done to determine why some haloesform Pop III stars and others do not, to provide insightabout the formation process of these stars and thus how theywill affect their surrounding pre-galactic environment. Theeffect of the LW intensity on host halo mass should be sys-tematically investigated with H self-shielding to determinewhether or not this background radiation really has a pro-found effect on Pop III formation or not. Further constraintson these primordial objects will assist multi-messenger as-tronomers in identifying these ancient systems. ACKNOWLEDGMENTS
JHW is supported by National Science Foundationgrants AST-1614333 and OAC-1835213, and NASA grantNNX17AG23G. The simulation was performed on Blue Wa-ters operated by the National Center for SupercomputingApplications (NCSA) with PRAC allocation support by theNSF (awards ACI-1514580 and OAC-1810584). The subse-quent analysis was performed with NSF’s XSEDE allocationAST-120046 on the Comet resource and also on the Geor-gia Tech PACE compute system. This research is part of theBlue Waters sustained-petascale computing project, which issupported by the NSF (awards OCI-0725070, ACI-1238993)and the state of Illinois. Blue Waters is a joint effort of theUniversity of Illinois at Urbana-Champaign and its NCSA.The freely available plotting library matplotlib (Hunter2007) was used to construct numerous plots within this pa-per. Computations and analysis described in this work wereperformed using the publicly-available enzo and yt codes,which is the product of a collaborative effort of many in-dependent scientists from numerous institutions around theworld. REFERENCES
Abel T., Wandelt B. D., 2002, MNRAS, 330, L53Abel T., Anninos P., Zhang Y., Norman M. L., 1997, New As-tronomy, 2, 181Abel T., Bryan G. L., Norman M. L., 2002, Science, 295, 93Anninos P., Zhang Y., Abel T., Norman M. L., 1997, New As-tronomy, 2, 209Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013, ApJ, 762, 109Bromm V., Coppi P. S., Larson R. B., 2002, ApJ, 564, 23Bryan G. L., et al., 2014, ApJS, 211, 19Chiaki G., Yoshida N., Hirano S., 2016, MNRAS, 463, 2781Chiaki G., Susa H., Hirano S., 2018, MNRAS, 475, 4378Clark P. C., Glover S. C. O., Smith R. J., Greif T. H., KlessenR. S., Bromm V., 2011, Science, 331, 1040Couchman H. M. P., 1991, ApJ, 368, L23Draine B. T., Bertoldi F., 1996, ApJ, 468, 269Efstathiou G., Davis M., White S. D. M., Frenk C. S., 1985, ApJS,57, 241Ezzeddine R., et al., 2019, ApJ, 876, 97Fernandez R., Bryan G. L., Haiman Z., Li M., 2014, MNRAS,439, 3798Field G. B., Somerville W. B., Dressler K., 1966, Annual Reviewof Astronomy and Astrophysics, 4, 207Glover S. C. O., Brand P. W. J. L., 2001, MNRAS, 321, 385Gnedin N. Y., Kravtsov A. V., 2006, ApJ, 645, 1054Greif T. H., Glover S. C. O., Bromm V., Klessen R. S., 2010, ApJ,716, 510Greif T. H., White S. D. M., Klessen R. S., Springel V., 2011a,ApJ, 736, 147Greif T. H., Springel V., White S. D. M., Glover S. C. O., ClarkP. C., Smith R. J., Klessen R. S., Bromm V., 2011b, ApJ,737, 75Griffen B. F., Dooley G. A., Ji A. P., O’Shea B. W., G´omez F. A.,Frebel A., 2018, MNRAS, 474, 443Hahn O., Abel T., 2011, MNRAS, 415, 2101Hartwig T., Glover S. C. O., Klessen R. S., Latif M. A., VolonteriM., 2015, MNRAS, 452, 1233Hartwig T., Volonteri M., Bromm V., Klessen R. S., Barausse E.,Magg M., Stacy A., 2016, MNRAS, 460, L74Heger A., Woosley S. E., 2002a, ApJ, 567, 532Heger A., Woosley S. E., 2002b, ApJ, 567, 532Heger A., Fryer C. L., Woosley S. E., Langer N., Hartmann D. H.,2003, ApJ, 591, 288Hirano S., Bromm V., 2017, MNRAS, 470, 898Hirano S., Hosokawa T., Yoshida N., Omukai K., Yorke H. W.,2015, MNRAS, 448, 568Hirano S., Hosokawa T., Yoshida N., Kuiper R., 2017, Science,357, 1375Hosokawa T., Omukai K., Yoshida N., Yorke H. W., 2011, Science,334, 1250Hosokawa T., Hirano S., Kuiper R., Yorke H. W., Omukai K.,Yoshida N., 2016, ApJ, 824, 119Hunter J. D., 2007, Computing In Science & Engineering, 9, 90Jaacks J., Finkelstein S. L., Bromm V., 2019, MNRAS, 488, 2202Jeon M., Pawlik A. H., Bromm V., Milosavljevi´c M., 2014a, MN-RAS, 440, 3778Jeon M., Pawlik A. H., Bromm V., Milosavljevi´c M., 2014b, MN-RAS, 444, 3288Kitayama T., Yoshida N., Susa H., Umemura M., 2004, ApJ, 613,631Machacek M. E., Bryan G. L., Abel T., 2001, ApJ, 548, 509Magg M., Hartwig T., Glover S. C. O., Klessen R. S., WhalenD. J., 2016, MNRAS, 462, 3591Magg M., Hartwig T., Agarwal B., Frebel A., Glover S. C. O.,Griffen B. F., Klessen R. S., 2018, MNRAS, 473, 5308Mebane R. H., Mirocha J., Furlanetto S. R., 2018, MNRAS, 479,4544 MNRAS , 1–13 (2019) radles of the first stars Muratov A. L., Gnedin O. Y., Gnedin N. Y., Zemp M., 2013,ApJ, 772, 106Naoz S., Yoshida N., Gnedin N. Y., 2012, ApJ, 747, 128Naoz S., Yoshida N., Gnedin N. Y., 2013, ApJ, 763, 27Nomoto K., Tominaga N., Umeda H., Kobayashi C., Maeda K.,2006, Nuclear Physics A, 777, 424O’Leary R. M., McQuinn M., 2012, ApJ, 760, 4O’Shea B. W., Norman M. L., 2007, ApJ, 654, 66O’Shea B. W., Norman M. L., 2008, ApJ, 673, 14Omukai K., Hosokawa T., Yoshida N., 2010, ApJ, 722, 1793Planck Collaboration et al., 2014, A&A, 571, A16Press W. H., Schechter P., 1974, ApJ, 187, 425Ricotti M., Gnedin N. Y., Shull J. M., 2001, ApJ, 560, 580Ritter J. S., Safranek-Shrader C., Milosavljevi´c M., Bromm V.,2016, MNRAS, 463, 3354Schaerer D., 2002, A&A, 382, 28Schaerer D., 2003, A&A, 397, 527Schauer A. T. P., Whalen D. J., Glover S. C. O., Klessen R. S.,2015, MNRAS, 454, 2441Schauer A. T. P., et al., 2017, MNRAS, 467, 2288Schauer A. T. P., Glover S. C. O., Klessen R. S., Ceverino D.,2019, MNRAS, 484, 3510Schneider R., Omukai K., Inoue A. K., Ferrara A., 2006, MNRAS,369, 1437Sheth R. K., Mo H. J., Tormen G., 2001, MNRAS, 323, 1Smith B. D., Wise J. H., O’Shea B. W., Norman M. L., KhochfarS., 2015, MNRAS, 452, 2822Stacy A., Bromm V., Lee A. T., 2016, MNRAS, 462, 1307Stecher T. P., Williams D. A., 1967, ApJ, 149, L29Susa H., 2013, ApJ, 773, 185Susa H., Hasegawa K., Tominaga N., 2014, ApJ, 792, 32Tegmark M., Silk J., Rees M. J., Blanchard A., Abel T., Palla F.,1997, ApJ, 474, 1Trenti M., Stiavelli M., 2009, ApJ, 694, 879Trenti M., Stiavelli M., Shull J. M., 2009, ApJ, 700, 1672Tseliakhovich D., Barkana R., Hirata C. M., 2011, MNRAS,p. 1501Tumlinson J., 2010, ApJ, 708, 1398Turk M. J., Abel T., O’Shea B., 2009, Science, 325, 601Turk M. J., Smith B. D., Oishi J. S., Skory S., Skillman S. W.,Abel T., Norman M. L., 2011, ApJS, 192, 9Visbal E., Haiman Z., Terrazas B., Bryan G. L., Barkana R., 2014,MNRAS, 445, 107Visbal E., Haiman Z., Bryan G. L., 2018, MNRAS, 475, 5246Wise J. H., Abel T., 2007, ApJ, 671, 1559Wise J. H., Abel T., 2008, ApJ, 685, 40Wise J. H., Abel T., 2011, MNRAS, 414, 3458Wise J. H., Cen R., 2009, ApJ, 693, 984Wise J. H., Abel T., Turk M. J., Norman M. L., Smith B. D.,2012, MNRAS, 427, 311Wise J. H., Regan J. A., O’Shea B. W., Norman M. L., DownesT. P., Xu H., 2019, Nature, 566, 85Wolcott-Green J., Haiman Z., Bryan G. L., 2011, MNRAS,p. 1673Xu H., Wise J. H., Norman M. L., 2013, ApJ, 773, 83Yoshida N., Abel T., Hernquist L., Sugiyama N., 2003, ApJ, 592,645This paper has been typeset from a TEX/L A TEX file prepared bythe author.MNRAS000