The influence of streaming velocities and Lyman-Werner radiation on the formation of the first stars

The first stars in the Universe, the so-called Population III stars, form in small dark matter minihaloes with virial temperatures $T_{\rm vir}<10^{4}$~K. Cooling in these minihaloes is dominated by molecular hydrogen (H$_{2}$), and so Population III star formation is only possible in those minihaloes that form enough H$_{2}$ to cool on a short timescale. As H$_{2}$ cooling is more effective in more massive minihaloes, there is therefore a critical halo mass scale $M_{\rm min}$ above which Population III star formation first becomes possible. Two important processes can alter this minimum mass scale: streaming of baryons relative to the dark matter and the photodissociation of H$_{2}$ by a high redshift Lyman-Werner (LW) background. In this paper, we present results from a set of high resolution cosmological simulations that examine the impact of these processes on $M_{\rm min}$ and on $M_{\rm ave}$ (the average minihalo mass for star formation), both individually and in combination. We show that streaming has a bigger impact on $M_{\rm min}$ than the LW background, but also that both effects are additive. We also provide fitting functions quantifying the dependence of $M_{\rm ave}$ and $M_{\rm min}$ on the streaming velocity and the strength of the LW background.


INTRODUCTION
The transition from metal-free to metal-enriched star formation marks a fundamental period in the early Universe, taking place at redshifts z 10. The so-called Population III (Pop III) stars are the first stellar objects that formed out of a gas consisting only of hydrogen, helium, and negligible amounts of lithium, making their formation conditions different from present-day star formation (Yoshida, Hosokawa & Omukai 2012;Bromm 2013;Glover 2013;Klessen 2019).
These first stars are out of reach of current and upcoming telescopes. Although their emission peaks in the restframe ultraviolet, this is shifted into the near-infrared at the present day because of the high redshifts of these stars. This greatly hampers efforts to detect them with groundbased telescopes, owing to the high thermal background of the atmosphere, putting them out of reach not only of current 8 m class telescopes but also the next generation of 30 m class facilities. Space-based telescopes avoid the prob-E-mail: anna.schauer@utexas.edu † Hubble Fellow lems posed by the atmosphere but current and near-future examples (e.g. HST, JWST) are too small to observe the first stars directly (Zackrisson et al. 2011(Zackrisson et al. , 2015Jeon & Bromm 2019;Schauer, Drory & Bromm 2020), with only the supernova explosions that mark the end of the lives of massive Pop III stars potentially being detectable (Rydberg et al. 2020).
Another method for exploring the properties of the first stars comes from 21 cm radiation (Furlanetto, Oh & Briggs 2006). Before reionization, the hyperfine transition of atomic hydrogen can give insight into the global state of the gas, which depends on the properties of the first stars and galaxies. A number of different experiments around the world are currently attempting to measure this signal (Singh et al. 2017;Fialkov, Barkana & Cohen 2018;Price et al. 2018;Philip et al. 2019), with a recent success reported by EDGES (Bowman et al. 2018). Measurements of the high-redshift 21 cm signal allow us to infer a number of important quantities, such as the high z star formation rate density, the X-ray emissivity of the first star-forming systems, or the minimum halo mass required for Pop III formation (Fialkov &   . The parameter space, however, is large, and often allows for more than one solution. Cosmological hydrodynamical simulations are required to understand the nature of the first stars and star-forming regions better. These simulations have been performed for over twenty years. From these simulations, we have learned that the first stars are generally more massive than present day stars, with a clump mass of 100-1000M , (Bromm, Coppi & Larson 1999;Abel, Bryan & Norman 2002). Some more modern simulations see fragmentation and multiple lower-mass cores (Clark et al. 2011b;Greif et al. 2012;Stacy & Bromm 2014;Wollenberg et al. 2020), with a high binary fraction (Stacy, Greif & Bromm 2010;Stacy & Bromm 2013), while others still find the high-mass formation mode to be the most prominent one (Hirano et al. 2014;Susa, Hasegawa & Tominaga 2014).
An important open question that we can address with the help of cosmological simulations is the mass scale of the first dark matter haloes to host Pop III star formation. There is broad agreement that the first stars formed in H2-cooled dark matter "minihaloes", with masses of a few 10 5 M -10 7 M and formation redshifts z ∼ 15-30 (see e.g. Glover 2013. However, star formation in these systems is strongly influenced by the strength of the highredshift Lyman-Werner (LW) background and by the size of the relative streaming velocity between the baryons and the dark matter.
LW photons between 11.2 eV and 13.6 eV are produced in abundance by young massive stars and can dissociate H2 via the two-step Solomon process (Field, Somerville & Dressler 1966;Stecher & Williams 1967). The presence of a high-redshift LW background produced by the earliest Pop III stars results in a reduction in the amount of H2 present in minihaloes, suppressing cooling and star formation in them to an extent that depends on the mass of the minihalo and the strength of the LW background (see e.g. Machacek, Bryan & Abel 2001;O'Shea & Norman 2008;Hirano et al. 2015).
The streaming velocity is a large-scale difference between the motion of the gas and the dark matter in the Universe, originating from the coupling between baryons and radiation prior to the recombination epoch (Tseliakhovich & Hirata 2010a). It also acts to suppress star formation in the lowest mass minihaloes, shifting the onset of star formation to higher mass haloes (see e.g. Greif et al. 2011;Naoz, Yoshida & Gnedin 2012;Naoz, Yoshida & Gnedin 2013;Hirano et al. 2017;. With a one sigma velocity dispersion σrms ≈ 30 km s −1 , these velocities are supersonic in most of the Universe following the decoupling of baryons and radiation during the recombination epoch. As the Universe expands, the size of the streaming velocities decreases with decreasing redshift as σrms ∝ (1 + z). Therefore, at the onset of Pop III star formation (z ≈ 20), the one sigma velocity dispersion is approximately σrms ≈ 0.6 km s −1 . It is negligible in the present day Universe.
Although there have been many studies of both effects in isolation, until now there has been no direct numerical study of the combination of both effects, leaving their relative importance unclear. This is the gap that this paper attempts to fill. We carry out a series of high-resolution cosmological simulations in which we vary the strength of the LW background and the size of the initial streaming velocity. We consider three different values for the radiation field strength and four different values for the streaming, leading to 12 simulations in total. Each simulation is followed down to a redshift of z ≈ 14, by which time each has formed several thousand minhaloes, several hundred of which are able to form stars. This provides us with a large sample of haloes in which to examine the effects of the LW background and the baryon-dark matter streaming, both in isolation and together.
Our paper is structured as follows. Section 2 presents the numerical method used for our simulations, including details of our choice of chemical network, our treatment of LW radiation and the streaming velocity, and how we identify haloes. In this section, we also describe the initial conditions we adopt for our simulations. In Section 3, we show qualitatively how a LW background and streaming velocities impact the large-and small-scale physical properties of the early Universe. Section 4 is focused on the key quantity of our investigation: the typical mass for a star-forming minihalo and how this depends on the strength of the LW background and size of the streaming velocity. We also provide highly-accurate fitting functions that allows one to easily determine this mass scale as a function of both parameters. We compare our findings to previous results in the literature in Section 5 and conclude in Section 6.

Numerical method
We carry out our simulations using the arepo moving-mesh cosmological hydrodynamics code (Springel 2010). arepo is based on a quasi-Lagrangian approach in which the fluid is discretised with an unstructured mesh that is the Voronoi tessellation of a set of mesh-generating points that move with the flow of the gas. The flux of mass, momentum, energy etc. between cells in this mesh is determined by solving the Riemann problem at each interface between mesh cells. These can be refined according to various criteria simply by adding additional mesh generating points and appropriately partitioning the gas between the newly-created cells. The refinement criteria used in our study are discussed in Section 2.3 below. Dark matter is represented in the simulations by discrete collisionless particles for which the gravitational forces are computed using a Barnes & Hut (1986) oct tree. The gravitational force is softened as described in Springel (2010), with a fixed softening length l soft = 20 comoving pc for the dark matter and an adaptive approach that scales with the size of the mesh cells for the gas.
The version of arepo that we use in our study includes a model for primordial gas chemistry and cooling. Our chemical network tracks the non-equilibrium abundances of 9 species (H, H + , H2, D, D + , HD, He, He + and e − ) as well as the equilibrium abundances of the H − and H + 2 ions. It is based on the network used in , but has been supplemented with a treatment of the effects of a soft ultraviolet background. We assume that below 13.6 eV, this background has the spectral shape of a 10 5 K blackbody, as would be the case if it were produced purely by massive Pop III stars. Above 13.6 eV, we set the strength of the radiation field to zero to account for absorption by the intergalactic medium. The overall photon flux is quantified by J21, the mean specific intensity just below the Lyman limit in units of 10 −21 erg s −1 cm −2 sr −1 Hz −1 . The influence of the radiation background on the gas chemistry is accounted for by including the following reactions in our chemical network with J21-dependent rates taken from Glover & Jappsen (2007) (for H2 and HD), Glover & Savin (2009) (for H − ) and Glover (2015) (for H + 2 ).The most important of the new reactions is the photodissociation of H2 by LW photons in the 11.2 -13.6 eV energy range, and so for simplicity we will often refer to our adopted radiation background as a LW background, although in reality it extends over a wider range of energies than simply the LW bands.
The rate at which H2 is destroyed by LW photons is highly sensitive to the amount of self-shielding occurring in the gas (Draine & Bertoldi 1996). We model the effects of H2 self-shielding using the TreeCol algorithm . This uses information on the H2 mass of each Voronoi cell and tree node stored in the Barnes-Hut tree to compute an approximate map of the H2 column densities surrounding each cell, pixelated using the healpix algorithm (Górski et al. 2005) with Npix = 48 equal-area pixels. In this study, we use the modified version of the TreeCol algorithm introduced by Hartwig et al. (2015) that also accounts for the velocity of the gas. In this version of the algorithm, the H2 in any given node or cell only contributes to the column density map if the centre-of-mass velocity of the node or cell, vcom, satisfies where vcurr is the velocity of the current cell (i.e. the one for which we are computing the column density map), and where v th is the thermal velocity of the H2 in that cell. This modification to TreeCol accounts for the fact that H2 that is highly red or blue-shifted relative to the cell of interest will not contribute to the shielding of that cell, because its LW absorption lines will be significantly Doppler-shifted relative to those in the cell of interest. In low mass minihaloes with velocity dispersions ∼ v th , this effect is of limited importance, but it becomes increasingly relevant as we consider more massive minihaloes with higher characteristic velocity dispersions.
Once we have the TreeCol-derived H2 column density map for each cell, we can then compute the effective selfshielding factor using a self-shielding function from Draine & Bertoldi (1996) where xi = NH 2 ,i/5 × 10 14 cm −2 , with NH 2 ,i being the H2 column density in pixel i of the column density map, yi = (1+xi) 0.5 , and where b5 = b/10 5 cm s −1 is the scaled Doppler parameter of the molecular hydrogen in the cell, which we assume to be dominated by thermal broadening. The H2 photodissociation rate in the cell is then simply given by where R diss,thin = 1.38 × 10 −12 J21 s −1 is the value in the optically thin limit. The Draine & Bertoldi (1996) self-shielding function assumes that H2 is rotationally cold (i.e. that all of the molecules are in their ortho or para ground states). This is a good approximation at the low densities that we are primarily concerned with in this study, but may over-estimate the effectiveness of self-shielding in warm, dense gas (Wolcott-Green, Haiman & Bryan 2011;Wolcott-Green & Haiman 2019). However, we do not expect this to impact our results, since in order for the gas to collapse to the densities at which the Draine & Bertoldi (1996) function becomes inappropriate, it must already have formed enough H2 to provide effective cooling.

Initial conditions
The basic setup of our simulations is the same as in . Briefly, we adopt a ΛCDM cosmological model, with parameters h = 0.6774, Ω0 = 0.3089, Ω b = 0.04864, ΩΛ = 0.6911, n = 0.96 and σ8 = 0.8159 (Planck Collaboration 2016). The simulations are initialised at a redshift z = 200 in a box of size 1h −1 Mpc in comoving units. Initial conditions for the dark matter are generated using music (Hahn & Abel 2011) with the transfer function of Eisenstein & Hu (1998). The density distribution of the baryons is assumed to initially trace that of the dark matter. This approach is not entirely correct for runs with non-zero streaming velocities (see e.g. Naoz & Barkana 2005;Naoz, Barkana & Mesinger 2009;Naoz, Yoshida & Barkana 2011), but the impact of this simplification on our results should be small (see e.g. the discussion in  or the recent results of Park et al. 2020). In runs without streaming, the velocity field of the baryons is the same as that of the dark matter. In the runs with streaming, this is included as a constant velocity offset, arbitrarily chosen to be in the x-direction. Dark matter is represented by 1024 3 particles with a mass of M dm = 99M . The gas is initially distributed across 1024 3 Voronoi mesh cells, each of which has a mass M cell 19 M . In the subsequent evolution we allow for further mesh refinement in regions where gravitational contraction sets in, as discussed in more detail in Section 2.3. In  we demonstrated that this resolution is high enough to yield converged results for the critical minihalo mass required for efficient H2 cooling.
The initial chemical state of the gas in our simulations is specified in Table 2. The fractional abundances listed there are defined as xi = ni/n, where ni is the number density of the chemical species and n is the number density of hydrogen nuclei.
In the study presented in this paper, we vary two free parameters: the strength of the Lyman-Werner background, as quantified by J21, and the size of the streaming velocity. For the strength of the Lyman-Werner background, we consider three values: no background (i.e. J21 = 0), J21 = 0.01 and J21 = 0.1. Our choice of these particular non-zero values is motivated by the fact that we expect them to bracket the average value present in the Universe at redshifts 1 + z ∼ 15  For the size of the streaming velocity, we consider values vstr = 0, 6, 12 and 18 km s −1 at z = 200, corresponding to 0, 1, 2 and 3 times σrms, the root-mean-squared streaming velocity at that redshift. The combinations of parameters that we examine are summarised in Table 1. In the runs with non-zero J21, we represent the redshift dependence of the background with a simple step function, where JLW is the mean specific intensity of the LW background in units of erg s −1 cm −2 Hz −1 sr −1 . Although somewhat unrealistic, this choice allows for a much cleaner numerical experiment than would be the case for a more complex dependence on z.

Mesh refinement and sink particles
By default, arepo attempts to ensure that the gas mass in each Voronoi mesh cell remains within a factor of two of its initial value, which in our case is M cell 19M . If necessary, it will refine the mesh cells to ensure that this condition is met by adding additional mesh-generating points. We supplement this default refinement scheme with an additional Jeans criterion, and refine the grid as necessary to ensure that the local Jeans length is always resolved with at least 16 cells (for further details, see e.g. Federrath et al. 2010).
For efficiency reasons, we do not allow gravitationally collapsing gas to reach arbitrarily high densities. Instead, we replace regions of dense, collapsing gas with non-gaseous sink particles. The sink particle implementation we use here is the same as that introduced in Tress et al. (2020) and Wollenberg et al. (2020). Mesh cells denser than our sink particle creation threshold, n th = 10 5 cm −3 , are candidates for sink particle formation. However, in order to actually form a sink, the gas in the cell and its surroundings must pass an additional series of checks. It must be gravitationally bound and collapsing (as assessed based on the divergence of the local velocity and acceleration of the gas). It must also be located at a local minimum in the gravitational potential. Finally, sink formation is disabled in gas cells located within the accretion radius of an existing sink.
Once formed, sink particles can accrete gas from their surroundings. Gas within a given sink accretion radius racc and above the threshold density n th will be accreted by the sink provided that it is gravitationally bound to the sink particle. We adopt racc ≈ 1 physical pc and n th = 10 5 cm −3 , because these values roughly reflect the bulk properties of the star-forming central region of the halo (e.g. Bromm, Coppi & Larson 1999;Abel, Bryan & Norman 2002), and because we do not need to follow the full fragmentation process and the formation of individual stars in the current investigation. This would require resolving much smaller scales and higher densities (e.g. Haemmerlé et al. 2020). As in Tress et al. (2020) and Wollenberg et al. (2020), we do not remove all of the gas from cells satisfying this criterion, but instead simply remove enough gas to reduce the cell density to n th . The mass and momentum of the accreted gas are added to the sink mass and momentum.

Halo selection
Our haloes are chosen the same way as in . We use the standard friends-of-friends algorithm to identify dark matter haloes. Gas cells are included as secondary components. A subfind algorithm (Springel et al. 2001) is run in addition, which we use to determine the centre of the halo.
At every snapshot, we check all minihaloes with masses M halo 5 × 10 4 M and investigate whether these haloes fulfil our star formation criteria. These are the same as in . We consider a halo to be capable of forming stars if at least one gas cell has: • a number density n 10 2 cm −2 , • a temperature T 500 K, • a molecular hydrogen abundance of xH 2 10 −4 .
In contrast to our procedure in Schauer et al. (2019), we measure the mass of the star-forming halo only in the snapshot when these criteria have been met for the first time. As all DM particles have unique IDs, we are able to trace haloes through the entire simulation and can build mergertrees. With the help of these merger trees, we are able to identify when haloes fulfil the star formation criteria for the first time. If we did not use merger trees, our results would be biased to higher masses, as in the absence of feedback, almost all haloes retain their star-forming ability at lower redshifts but accumulate more mass through accretion and mergers.
Our sink particle formation criterion is more stringent than the star formation criterion described above, as sink formation occurs at a much higher density than the value we use to flag halos capable of star formation. Consequently, in any given halo that ultimately forms stars, the star formation criterion is fulfilled before sink particle formation occurs. This prompts the question of whether the two criteria actually identify the same set of halos: do all of the halos we flag as "capable of forming stars" go on to form sink particles? With the merger tree, we are able to trace back star-forming halos in time and associate the sink particles with their birth halo locations and times, allowing us to address this point. A full extended analysis that focuses on the spread of the star-forming halo masses, as well as the merger history, lies beyond the scope of this paper and will be addressed in future work. Nevertheless, a comparison of the two criteria for our most representative simulation (v1 lw-2) reveals that all all sink particles form in halos previously flagged as capable of star formation. The time delay between the fulfilment of the star formation criterion and the actual formation of a sink particle is on average one timestep, ∆z = 1, and all star-forming halos flagged between redshifts z = 20 and z = 17 have formed a sink particles before the end of the simulation at z = 14.

FROM LARGE TO SMALL SCALES: A QUALITATIVE ANALYSIS
As discussed in the Introduction, both the streaming velocity and the LW background suppress Pop III star formation through different mechanisms. We start by giving a qualitative analysis of our simulations, and proceed by discussing the effects that these processes have on large scales (comparable to the size of the simulation box) and small scales (comparable to the size of a halo).

The effects of a non-zero streaming velocity
In Figure 1, we show slices of the gas number density at redshift z = 15 in four simulations with different values for the streaming velocity, progressively zooming into the centre of a halo (top to bottom row). Our four simulations shown here, v0 lw0 (first column), v1 lw0 (second column), v2 lw0 (third column), and v3 lw0 (fourth column), all have a zero LW background. In the top row, a density slice through the entire computational volume can be seen for each of these simulations. While the density structure in the image on the left is well defined and crisp, it becomes more blurred going to the right. In the rightmost panel, which shows our highest streaming velocity value of 3σrms the gas distribution is quite fuzzy and washed out. The effect is especially strong in the x-direction, which corresponds to the direction of the streaming velocity in the initial conditions On these large scales, one can directly observe the effect of streaming velocities: they smooth out the gas distribution. The middle row of Figure 1 shows a close-up with a side length of 20 ckpc/h, centred around a minihalo (the size of the region is indicated in the top row). For this zoom-in view, we chose exactly the same x, y and z coordinates in all simulations to be able to better compare the effects of the different streaming velocities. As the streaming velocity increases, the highest density region shifts to the right, in the direction of the streaming velocity.
As seen on the large scales, streaming velocities lead to the washing out of the gas structures, whereas the dark matter structure remains largely unchanged. All structure in the Universe forms in a hierarchical manner: first, sheets and filaments emerge, before they collapse into haloes, with smaller haloes forming first. For a large time interval in the early Universe, these relative velocities between gas and dark matter were larger than the escape velocities of the dark matter potentials of the filaments and haloes that began emerging through gravitational interactions. As a result, the gas follows the dark matter structure later when the streaming velocity is higher, and we can see a lower density contrast in the slice plots (e.g. compare to Naoz & Narayan 2014, who find a phase shift between the linear collapse modes of gas and dark matter).
As a consequence of this behaviour, the gas fraction in the haloes is reduced (Tseliakhovich & Hirata 2010b;Popa et al. 2016;, as is the mean gas density within the virial radius. Even the dark matter power spectrum is slightly reduced Ahn & Smith 2018).
This directly impacts the ability of the gas in the haloes to cool and form stars. Haloes formed in regions of the Universe with high streaming velocities have less gas than those formed in regions with lower streaming velocities, and the gas that they do retain is less dense overall, and hence less able to form H2 and cool. This can be seen clearly in the bottom row of Figure 1, where we show a 2 ckpc/h close-up on the minihalo, centred around the highest density region (indicated by the white lines in the middle row). One immediately sees that the maximum density decreases with increasing streaming velocity. As star formation requires gas to reach high densities, it is not surprising that high streaming velocities hinder and delay the process. As we will show in the next subsection, only more massive haloes are able to retain enough high-density gas for Pop III star formation to proceed when the streaming velocity is large.

The effects of a Lyman-Werner background
LW radiation is another way to suppress star formation. By photodissociating molecular hydrogen, the main coolant at high redshifts is destroyed. We therefore investigate the abundance of molecular hydrogen in Figure 2. Similar to Figure 1, the top row shows a slice of the entire simulation box, progressively zoomed into one halo in the middle and bottom rows. As it is the most common of the streaming velocity values investigated here (see Section 4.3), we restrict our discussion to the 1 σrms case, and show simulations with no LW background (left column), a weak LW background (middle column), and a larger LW background (right column).
On large scales and in low density regions, molecular hydrogen is almost immediately destroyed. In the top row of Figure 2, the molecular hydrogen abundance in the intergalactic medium drops from a few 10 −6 in the run with no LW background to below 10 −9 in both runs with a non-zero LW background. Molecular hydrogen, however, can self-shield, so in high-density regions, the abundance stays higher. This is illustrated by the few pink and yellow regions in the middle and right top panels.
Zooming into one halo, as indicated by the white lines, one can identify the halo centres (compare the number density slice of these simulations in Figure 3) by their increased molecular hydrogen abundances. In case of a zero LW background, the abundance is highest and the most extended, but even for the strong LW background with J21 = 0.1, an H2 abundance of more than 10 −4 is reached. This immediately demonstrates that the H2 in this halo is able to selfshield effectively. Nevertheless, the peak H2 abundance in the runs with a non-zero LW background is clearly smaller than in the run with no LW background, reducing its effectiveness as a coolant. The impact of this on the density distribution inside the halo can been seen in Figure 3 (top row): the central density is slightly reduced in the runs with J21 > 0 compared to the case with no LW background (compare the second to the fifth and sixth panels).
Comparing our results here to those from the runs with high streaming but no LW, we see that the manner in which these two processes suppress cooling and star formation is quite different. Streaming reduces the gas density throughout the haloes, which has the knock on effect of making it harder for the gas to form H2 and harder to cool once it has formed H2. The LW background, on the other hand, has little effect on the gas density on halo scales and hence does not affect the ability of the halo to form H2. Instead, it suppresses cooling by destroying most of the H2 that does form, leaving less available to cool the gas.

PROPERTIES OF STAR-FORMING MINIHALOES
As our next step, we move to a more quantitative analysis. Our goal is to understand when minihaloes form Pop III stars and how the two effects, baryonic streaming and an elevated LW background, interact.

Suppression of Pop III star formation
We start by investigating the number of star-forming haloes as a function of redshift and show our results in Figure 4. As expected, the number of star-forming haloes increases with decreasing redshift, as haloes grow over time and the gas in the centres has more time to cool and collapse. For the first time, we show a direct comparison between the streaming velocity and LW radiation. Our key finding is that the LW background only plays a minor role in reducing the number of star-forming haloes compared to the streaming velocities. There is little difference between the results with no LW background and those for a weak background with J21 = 0.01. The presence of a weak LW background reduces the number of star-forming haloes, but only by around 25% in most simulations. The exception is the simulation with v bc = 3σrms, where we see a larger reduction in the star-forming halo number at z = 14. However, the number of star-forming haloes in all of the 3σrms runs is so small that this could just be a consequence of small-number statistics.
For the stronger LW background with J21 = 0.1, we see a larger effect: the number of star-forming haloes is reduced by a factor of 3-4 at most redshifts, leading to a significant decrease in the amount of star formation occurring in the simulation.
Compared to this, the presence or absence of streaming velocities has a much larger impact on the number of star-forming haloes. Keeping the LW background constant, this number decreases on average by a factor of 5 for 1 σrms compared to the zero velocity case. This reduction reaches factors of 30 -40 when going from zero to a streaming velocity of 2 σrms, and can exceed values of 100 when increasing the streaming velocity to 3 σrms.
To put these numbers into context, note that the most representative value for v bc that we examine is v bc ≈ 1 σrms, as the volume fraction of streaming velocities peaks at 0.8 σrms (compare Figure 12), while a typical value for the LW background at the redshifts studied here is J21 ≈ 0.01 (Wise et al. 2012). Neglecting the effects of streaming velocities therefore has a stronger influence on the number of star-forming halos than neglecting the presence of an LW background when considering a typical region of the Universe. . Slice plots of the gas number density (top row), dark matter density (second row), gas temperature (third row) and H 2 abundance (bottom row) of the 20 ckpc/h cut-out region shown in the middle rows of Figures 1 and 2 at redshift z = 15. From left to right, we show the simulations v0 lw0, v1 lw0, v2 lw0, v3 lw0, v1 lw-2, and v1 lw-1. While the number density and gas structure in general changes to lower maximum values for increasing streaming velocity values, it is only decreased a little bit for an increasing LW background. This is reflected in the temperature slices where the increased temperatures that indicate shock heating from the formation of the halo and filaments, largely follow the density structure. The dark matter density structure remains largely unchanged. For the two simulations with a non-zero LW background, the molecular hydrogen abundance drops to below 10 −7 , with only the centres of the haloes retaining a higher fraction. This is reflected in a slightly lower maximum gas number density, as seen in the top row.

Halo masses of star-forming minihaloes
We investigate the halo masses of star-forming minihaloes, and by how much these masses are shifted to larger values for increasing streaming velocities and a LW background. This is a more robust measurement than the number of haloes that are forming stars per redshift bin, as that number depends on the box size and on the details of setting up the initial conditions.
We further consider all halos only at the first snapshot (taken at equidistant redshift intervals with ∆z = 1) in which they are forming stars. Pop III stars are massive and short-lived (Schaerer 2002), and one or two supernova explosions are enough to enrich the halo gas to metallicities of Z ∼ 10 −3.5 (Rossi, Salvadori & Skúladóttir 2021), high enough for Pop III star formation to transition into Pop II star formation (Bromm et al. 2001). By only counting halos once, we hence avoid biasing against higher mass halos that grew through accretion and mergers, and that might take Pop II instead of Pop III star formation into account.
In Figure 5, we show the mass of the least massive halo that fulfils the SF criterion as a function of redshift for all 12 simulations. In this plot and those that follow, we exclude data at redshift z 24, as we initialise the LW background at that redshift. If no new halo is found to be star forming at a given snapshot, we don't show a data point at this redshift. After the onset of star formation, this only happens in the simulations with the highest streaming velocity, 3 σrms.
In the simulation with no LW background and no streaming velocity, we find the lowest mass threshold: Mmin 3 × 10 5 M . Increasing the strength of the LW background increases Mmin, as does increasing the streaming velocity. However, once again we see that changes in the streaming velocity have a much bigger impact than changes in the strength of the LW background, at least for the range of LW background strengths considered here. v0_lw0 v0_lw-2 v0_lw-1 v1_lw0 v1_lw-2 v1_lw-1 v2_lw0 v2_lw-2 v2_lw-1 v3_lw0 v3_lw-2 v3_lw-1 Figure 5. The minimum halo mass of star-forming haloes at the given redshift. Halos are only accounted for when they first fulfil our star formation criterion. Both streaming velocities and a LWBG lead to halos of higher masses forming stars.
In all of the simulations, typically only a few minihaloes with masses M ∼ Mmin actually form stars, while most remain starless. Star formation only becomes common in minihaloes with masses that are a factor of a few greater than Mmin (see e.g. Schauer et al. 2019 for a more detailed discussion of this point). Therefore, as well as looking at Mmin, it is also useful to look at the impact of streaming and LW radiation on the average mass that a minihalo has at the point at which it first forms stars, Mave. This is plotted in Figure 6 as a function of redshift for all twelve simulations.
Here again, we see that streaming velocities play a stronger role in increasing Mave than the LW background. v0_lw0 v0_lw-2 v0_lw-1 v1_lw0 v1_lw-2 v1_lw-1 v2_lw0 v2_lw-2 v2_lw-1 v3_lw0 v3_lw-2 v3_lw-1 Figure 6. Same as Figure 5, but for the average halo mass of star-forming minihaloes at that redshift. As in Figure 5, halos are only accounted for when they first fulfil the star formation criteria. Again, both streaming velocities and a LWBG increase the average halo mass for star formation. The behaviour is very stochastic and no clear trend with redshift can be observed.
To better understand this behavior, we show in Figure 7 the full distribution of the masses of star-forming haloes in each simulation, selected at the point at which they first start forming stars. Note that this is not the same as selecting the haloes at a fixed point in time. In particular, the fact that many of these distributions include few or no haloes with M > 10 7 M does not imply that haloes of this mass are unable to form stars. Rather, it implies that all haloes of this mass have at least one lower mass progenitor that was already able to form stars. We see from the Figure that in most cases, the halo mass distributions are well fit by Gaussians. The exceptions are the runs with 3 σrms streaming, and the run with 2 σrms and J21 = 0.1, all of which produce too few star-forming haloes (typically 1-4) to let us draw any conclusions about the shape of the halo mass distribution. In each panel in the Figure, we indicate the average halo mass, Mave = µ, the standard deviation of the Gaussian fit to the mass distribution, σ, the standard error in the mean, E = σ/ √ N , and the number of star-forming haloes selected from the simulation, N . One can see again that the average halo mass shifts to higher values for larger streaming velocities and a stronger LW background. We also note that standard deviation of the distribution is roughly 1/4 dex for the 0 and 1 σrms-streaming velocity simulations, but slightly smaller for higher values. This decline is unexpected and could be a result of our small sample size for large streaming velocities. A more detailed analysis requires significantly larger numerical simulations which are beyond the scope of our current investigation.
The number of halos in our simulation is limited due to our box size of 1cMpc/h. In previous work , we have seen that there is about a one order of magnitude spread between the lest massive halo forming stars and the most massive halo that does not form stars. Especially in simulations with high streaming velocity values or a LWBG of 0.1, the upper end of the halo mass function is influenced by our box size.

The average mass at which minihaloes first form stars
As we have already seen, neither Figure 5 nor Figure 6 display a strong redshift dependence. Similar to our previous work , we conclude that the mass thresholds for star formation are largely independent of redshift. Therefore we can stack the data for each simulation from all redshifts in order to obtain better statistics. Our goal is to identify a general formula that links the mass of the typical star-forming halo to its environmental parameters, here the size of the streaming velocity and the strength of the LW background it is exposed to. The resulting two parameter metric is illustrated in Fig-ure 8, where we present the average mass at which a minihalo first forms stars as a function of the streaming velocity, and in Figure 9, where we show the average halo mass as a function of the LW background. The average mass rises linearly with the streaming velocity, and the intercept scales as the square root of the strength of the LW background. We find that most of the data can be well fit by a relation log 10 Mave = log 10 M0 + 0.4159 × v bc σrms , where the intercept log 10 M0 is described by and all masses are given in solar masses. This relation is valid for the range of values explored in our simulations, i.e. 0 σ bc v bc 3 σ bc and 0 J21 0.1. The minimum halo mass can be fitted with a function of the same shape. This time, however, the slope s varies significantly with the LW background value, and we apply a fit to the slope as well.
We evaluate the fit to the data in Figures 10 and 11 Figure 11. Same as Figure 10, but for the minimum halo mass as a function of streaming velocity for our three values of the LW background.
where we display the average (minimum) halo mass as a function of streaming velocity in the top panel and the residual in the bottom panel. We find that the fit (solid lines) matches the data (dashed lines) very well (within 4%) for zero and 1 σrms streaming velocities, and moderately well for 2 σrms streaming. It does not fit well for the 3 σrms streaming runs, but as we have already seen, these runs are strongly affected by small number statistics. With our fit function, the minimum halo mass exceeds the average halo mass for v bc ≈ 2.5σrms, a result of the poor number statistics for 3σrms. An additional comparison of our fit functions to the data is given in the appendix, in Figures A1. It is also important to note that the zero and 1 σrms cases are by far the most representative of the Universe. This is quantified in Figure 12, where we show differential (cumulative) volume filling fractions of the Universe as a function of streaming velocity, illustrated by the blue (orange) line. One can see that there is a relatively sharp peak at 0.8 σrms. We also include the accuracy (the deviation of the residual fraction from the perfect one-to-one fit) in this Figure as described by the light grey (a LW background with J21 = 0.1), the medium grey (a LW background with J21 = 0.01) and the black (no LW background) lines. In the low streaming velocity regions, the fit is very accurate, and overall, the fit represents the data very well (at least 96%).  man 2008). The majority of these studies found substantial suppression of H2 cooling for LW background field strength J21 = 0.1, in some tension with our findings that a background with this strength has a relatively small impact. The main reason for this difference appears to be the fact that most of these previous studies neglected the effects of H2 self-shielding, which we find to have a significant effect. Particularly informative in this regard is the study by Yoshida et al. (2003). In the absence of self-shielding, they find that a field strength of only J21 = 0.01 is enough to increase Mmin by more than a factor of two. On the other hand, if they account approximately for the effects of self-shielding, they find little difference between their results with J21 = 0 and J21 = 0.01, in good agreement with our findings. Our masses appear to be slightly larger than the values reported by Hirano et al. (2015), who include a selfconsistent LW background and an approximate treatment of self-shielding, but no streaming velocities, and who find haloes with virial masses of a few 10 5 − 10 6 M . However, it is difficult to directly compare the two simulations, as Hirano et al. (2015) do not report the values of Mmin they obtain as a function of J21. More recently, Skinner & Wise (2020) report results from a simulation of minihalo formation with a strongly time-varying LW background. They find values of Mmin similar to those in our study, lying mostly in the range 3 × 10 5 < M < 10 6 M . Their model includes the effects of H2 self-shielding and they find an even weaker dependence of Mmin on the strength of the LW background than we do. However, they quantify the LW background field strength in terms of the instantaneous value at the point when a minihalo forms stars. Since their background is strongly timevarying, this can often be significantly larger than the mean value seen by the minihalo during its lifetime, making direct comparison with our time-independent values difficult.

Comparison with previous results
Regarding the dependence of Mmin and Mave, we find, unsurprisingly, that there is good agreement between our current results and those of our previous study . In that paper, we compared our results to a number of previous studies from the literature and showed that there was generally good agreement between our values of Mmin and those found in previous simulations. A similar result holds for our current study.
Finally, regarding the combined effect of the LW background and streaming, we remind the reader that our simulations represent the first comprehensive study that treats both effects simultaneously. Between the submission of the paper and answering the referee report, there has been one other publication, targeting the minimum and average halo mass for star formation in minihalos (Kulkarni, Visbal & Bryan 2020). The authors find a factor of 2-3 smaller halo masses and a mild redshift dependence. We assume the difference results from a different choice of halo mass (Kulkarni, Visbal & Bryan 2020 chose the virial mass instead of the halo mass of all particles and gas cells the minihalo is composed of), as well as an inherently redshift-dependent star formation criterion, as their temperature criterion for star formation in a minihalo depends on the virial temperature. The details of these difference will be explored in future work.

Other effects
In addition to the LW background and the streaming velocity, another large-scale effect at high redshift that can potentially influence the formation of Pop III stars is the presence of an X-ray background. X-rays will be emitted by sources such as high redshift quasars or high-mass X-ray binaries, but it is unclear whether their effect will be to enhance or delay star formation: on the one hand, they heat up their environment and therefore suppress small-scale star formation, while on the other hand, the additional ionization they provide can catalyze the formation of H2 and therefore could enhance Pop III star formation (Glover & Brand 2003;Machacek, Bryan & Abel 2003;Jeon et al. 2012Jeon et al. , 2014. Ionizing radiation can also delay star formation, if it is strong enough to break out from a high-redshift galaxy and affect neighbouring minihaloes. This can be the case for close pairs of low-mass haloes even at high redshift (Chen et al. 2017), but will become much more widespread as we ap-proach the epoch of reionization (Ocvirk et al. 2016). However, both of these effects are beyond the scope of the current paper.

CONCLUSIONS
In this paper, we have investigated how streaming velocities and a LW background influence Pop III star formation in minihaloes. We find that streaming velocities play a more important role in delaying star formation and in increasing the minimum and average mass at which minihaloes first form stars that the presence of a LW background of the strength that we expect to find in the Universe during the epoch dominated by Pop III star formation. The much higher LW radiation field strength found in the vicinity of massive star-forming galaxies (see e.g. Ahn et al. 2009) will likely have a much more significant effect relative to streaming, but this is relevant only for a small number of Pop III star-forming systems and is not the common case.
We find that, independent of redshift, the average mass at which minihaloes first form stars increases from ∼ 10 6 M for no streaming to ∼ 3 × 10 7 M for 3σrms streaming, i.e. beyond the atomic cooling threshold. We also provide two simple fitting function that describes how the minimum and the average mass vary as a function of LW background and streaming velocity, for use in future semi-analytical models or cosmological simulations. . Minimum (left panels) and average (right panels) halo masses for star formation. We show both the data points from our simulations, in the same way we presented them in Figures 5 and 6, as well as the values obtained by our fits in Equations (9)-(13).