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

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_{\rm 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⁺, H₂, 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 Schauer et al. (2019), 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 $J_{21}$, the mean specific intensity just below the Lyman limit in units of $10^{-21}\>{\rm 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 $J_{21}$-dependent rates taken from Glover & Jappsen (2007) (for H₂ and HD), Glover & Savin (2009) (for H^-) and Glover (2015) (for H${}_{2}^{+}$).The most important of the new reactions is the photodissociation of H₂ 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 H₂ 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 H₂ self-shielding using the TreeCol algorithm (Clark, Glover & Klessen, 2012). This uses information on the H₂ mass of each Voronoi cell and tree node stored in the Barnes-Hut tree to compute an approximate map of the H₂ column densities surrounding each cell, pixelated using the healpix algorithm (Górski et al., 2005) with $N_{\rm pix}=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 H₂ in any given node or cell only contributes to the column density map if the centre-of-mass velocity of the node or cell, ${\bf v}_{\rm com}$, satisfies

where $v_{\rm curr}$ is the velocity of the current cell (i.e. the one for which we are computing the column density map), and where $v_{\rm th}$ is the thermal velocity of the H₂ in that cell. This modification to TreeCol accounts for the fact that H₂ 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 $\sim v_{\rm 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 H₂ column density map for each cell, we can then compute the effective self-shielding factor using a self-shielding function from Draine & Bertoldi (1996)

where $x_{i}=N_{\rm H_{2},i}/5\times 10^{14}\>{\rm cm^{-2}}$, with $N_{\rm H_{2},i}$ being the H₂ column density in pixel $i$ of the column density map, $y_{i}=(1+x_{i})^{0.5}$, and where $b_{5}=b/10^{5}\,{\rm 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 H₂ photodissociation rate in the cell is then simply given by

where $R_{\rm diss,thin}=1.38\times 10^{-12}\,J_{21}\>{\rm s^{-1}}$ is the value in the optically thin limit.

The Draine & Bertoldi (1996) self-shielding function assumes that H₂ 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 H₂ to provide effective cooling.

The basic setup of our simulations is the same as in Schauer et al. (2019). Briefly, we adopt a $\Lambda$CDM cosmological model, with parameters $h=0.6774$, $\Omega_{0}=0.3089$, $\Omega_{\rm b}=0.04864$, $\Omega_{\Lambda}=0.6911$, $n=0.96$ and $\sigma_{8}=0.8159$ (Planck Collaboration, 2016). The simulations are initialised at a redshift $z=200$ in a box of size $1h^{-1}\,{\rm 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 Schauer et al. 2019 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_{\rm dm}=99\>{\rm M_{\odot}}$. The gas is initially distributed across $1024^{3}$ Voronoi mesh cells, each of which has a mass $M_{\rm cell}\simeq 19\>{\rm M_{\odot}}$. 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 Schauer et al. (2019) we demonstrated that this resolution is high enough to yield converged results for the critical minihalo mass required for efficient H₂ cooling.

The initial chemical state of the gas in our simulations is specified in Table 2. The fractional abundances listed there are defined as $x_{i}=n_{i}/n$, where $n_{i}$ 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 $J_{21}$, and the size of the streaming velocity. For the strength of the Lyman-Werner background, we consider three values: no background (i.e. $J_{21}=0$), $J_{21}=0.01$ and $J_{21}=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\sim 15$ and above (see e.g. the models for the Lyman-Werner background presented in Haiman, Abel & Rees 2000, Trenti & Stiavelli 2009, Ahn et al. 2009, Wise et al. 2012 or Agarwal et al. 2012, all of which predict values for $J_{21}$ between these limits for redshifts $15<1+z<25$ for a broad range of input parameters). Although we expect $J_{21}$ to be higher at redshifts closer to the redshift of reionization, or in the immediate vicinity of massive star-forming galaxies (Ahn et al., 2009), these values are representative of the ones seen by most Pop III star-forming minihaloes during the range of redshifts in which Pop III star formation dominates. We also note that previous work (e.g. Haiman, Abel & Rees, 2000; Machacek, Bryan & Abel, 2001; Yoshida et al., 2003; O’Shea & Norman, 2008) suggests that these field strengths are large enough to yield significant suppression of H₂ cooling in high-redshift minihaloes.

For the size of the streaming velocity, we consider values $v_{\rm str}=0,6,12$ and $18\,{\rm km\,s^{-1}}$ at $z=200$, corresponding to 0, 1, 2 and 3 times $\sigma_{\rm 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 $J_{21}$, we represent the redshift dependence of the background with a simple step function,

where $J_{\rm LW}$ is the mean specific intensity of the LW background in units of ${\rm 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$.

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_{\rm cell}\simeq 19\>{\rm M_{\odot}}$. 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_{\rm th}=10^{5}\>{\rm 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 $r_{\rm acc}$ and above the threshold density $n_{\rm th}$ will be accreted by the sink provided that it is gravitationally bound to the sink particle. We adopt $r_{\rm acc}\approx 1$ physical pc and $n_{\rm 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_{\rm th}$. The mass and momentum of the accreted gas are added to the sink mass and momentum.

Our haloes are chosen the same way as in Schauer et al. (2019). 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_{\mathrm{halo}}\geq 5\times 10^{4}$ M_⊙ and investigate whether these haloes fulfil our star formation criteria. These are the same as in Schauer et al. (2019). We consider a halo to be capable of forming stars if at least one gas cell has:

• •

  a number density $n\geq 10^{2}$ cm^-2,

• •

  a temperature $T\leq 500$ K,

• •

  a molecular hydrogen abundance of $x_{\rm H_{2}}\geq 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 merger-trees. 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, $\Delta 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$.

In Figure 1, we show slices of the gas number density at redshift $z=1{}5$ 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$\sigma_{\mathrm{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; Schauer et al., 2019), as is the mean gas density within the virial radius. Even the dark matter power spectrum is slightly reduced (O’Leary & McQuinn, 2012; 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 H₂ 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.

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 $\sigma_{\mathrm{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 $J_{21}=0.1$, an H₂ abundance of more than $10^{-4}{}$ is reached. This immediately demonstrates that the H₂ in this halo is able to self-shield effectively. Nevertheless, the peak H₂ 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 $J_{21}>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 H₂ and harder to cool once it has formed H₂. 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 H₂. Instead, it suppresses cooling by destroying most of the H₂ that does form, leaving less available to cool the gas.

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 $J_{21}=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_{\rm bc}=3\sigma_{\mathrm{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\sigma_{\mathrm{rms}}$ runs is so small that this could just be a consequence of small-number statistics. For the stronger LW background with $J_{21}=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 $\sigma_{\mathrm{rms}}{}$ compared to the zero velocity case. This reduction reaches factors of 30 – 40 when going from zero to a streaming velocity of 2 $\sigma_{\mathrm{rms}}{}$, and can exceed values of 100 when increasing the streaming velocity to 3 $\sigma_{\mathrm{rms}}{}$.

To put these numbers into context, note that the most representative value for $v_{\rm bc}$ that we examine is $v_{\mathrm{bc}}\approx 1\,\sigma_{\mathrm{rms}}$, as the volume fraction of streaming velocities peaks at 0.8 $\sigma_{\mathrm{rms}}$ (compare Figure 12), while a typical value for the LW background at the redshifts studied here is $J_{21}\approx 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.

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 $\Delta 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\sim 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\geq 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 $\sigma_{\mathrm{rms}}$.

In the simulation with no LW background and no streaming velocity, we find the lowest mass threshold: $M_{\rm min}\simeq 3\times 10^{5}$ M_⊙. Increasing the strength of the LW background increases $M_{\rm min}$, 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.

In all of the simulations, typically only a few minihaloes with masses $M\sim M_{\rm min}$ 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 $M_{\rm min}$ (see e.g. Schauer et al. 2019 for a more detailed discussion of this point). Therefore, as well as looking at $M_{\rm min}$, 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, $M_{\rm ave}$. 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 $M_{\rm ave}$ than the LW background. 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}\>{\rm M_{\odot}}$ 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\,\sigma_{\rm rms}$ streaming, and the run with $2\,\sigma_{\rm rms}$ and $J_{21}=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, $M_{\mathrm{ave}}=\mu$, the standard deviation of the Gaussian fit to the mass distribution, $\sigma$, the standard error in the mean, $E=\sigma/\sqrt{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 $\sigma_{\mathrm{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 (Schauer et al., 2019), 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.

As we have already seen, neither Figure 5 nor Figure 6 display a strong redshift dependence. Similar to our previous work (Schauer et al., 2019), 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 Figure 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

where the intercept $\log_{10}M_{0}{}$ 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\,\sigma_{\mathrm{bc}}\leq v_{\mathrm{bc}}\leq 3\,\sigma_{\mathrm{bc}}$ and $0\leq J_{21}\leq 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.

with

and

Again, this relation is valid for the streaming velocity values and LW background values explored by these simulations, especially for $J_{21}\leq 0.1$.

We evaluate the fit to the data in Figures 10 and 11, 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 $\sigma_{\mathrm{rms}}$ streaming velocities, and moderately well for 2 $\sigma_{\mathrm{rms}}$ streaming. It does not fit well for the 3 $\sigma_{\mathrm{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_{\mathrm{bc}}\approx 2.5\sigma_{\mathrm{rms}}$, a result of the poor number statistics for $3\sigma_{\mathrm{rms}}$. An additional comparison of our fit functions to the data is given in the appendix, in Figures 13.

It is also important to note that the zero and 1 $\sigma_{\mathrm{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 $\sigma_{\mathrm{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 $J_{21}=0.1$), the medium grey (a LW background with $J_{21}=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%).

There have been a number of previous studies of the impact of a LW background on H₂ cooling in minihaloes (see e.g. Haiman, Abel & Rees, 2000; Machacek, Bryan & Abel, 2001; Yoshida et al., 2003; Wise & Abel, 2007; O’Shea & Norman, 2008). The majority of these studies found substantial suppression of H₂ cooling for LW background field strength $J_{21}=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 H₂ 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 $J_{21}=0.01$ is enough to increase $M_{\rm min}$ 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 $J_{21}=0$ and $J_{21}=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 self-consistent 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}\,\mathrm{M}_{\odot}$. However, it is difficult to directly compare the two simulations, as Hirano et al. (2015) do not report the values of $M_{\rm min}$ they obtain as a function of $J_{21}$. More recently, Skinner & Wise (2020) report results from a simulation of minihalo formation with a strongly time-varying LW background. They find values of $M_{\rm min}$ similar to those in our study, lying mostly in the range $3\times 10^{5}<M<10^{6}\>{\rm M_{\odot}}$. Their model includes the effects of H₂ self-shielding and they find an even weaker dependence of $M_{\rm min}$ 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 time-varying, 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.

Regarding the dependence of $M_{\rm min}$ and $M_{\rm ave}$, we find, unsurprisingly, that there is good agreement between our current results and those of our previous study (Schauer et al., 2019). 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 $M_{\rm min}$ 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.

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 H₂ and therefore could enhance Pop III star formation (Glover & Brand, 2003; Machacek, Bryan & Abel, 2003; Jeon et al., 2012, 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 approach the epoch of reionization (Ocvirk et al., 2016). However, both of these effects are beyond the scope of the current paper.