Relative baryon-dark matter velocities in cosmological zoom simulations

Supersonic relative motion between baryons and dark matter due to the decoupling of baryons from the primordial plasma after recombination affects the growth of the first small-scale structures. Large box sizes (greater than a few hundred Mpc) are required to sample the full range of scales pertinent to the relative velocity, while the effect of the relative velocity is strongest on small scales (less than a few hundred kpc). This separation of scales naturally lends itself to the use of `zoom' simulations, and here we present our methodology to self-consistently incorporate the relative velocity in zoom simulations, including its cumulative effect from recombination through to the start time of the simulation. We apply our methodology to a large-scale cosmological zoom simulation, finding that the inclusion of relative velocities suppresses the halo baryon fraction by $46$--$23$ per cent between $z=13.6$ and $11.2$, in qualitative agreement with previous works. In addition, we find that including the relative velocity delays the formation of star particles by $\sim 20 {~\rm Myr}$ Myr on average (of the order of the lifetime of a $\sim 9~{\rm M}_\odot$ Population III star) and suppresses the final stellar mass by as much as $79$ per cent at $z=11.2$.


INTRODUCTION
The cosmic microwave background (CMB) radiation carries an image of the Universe at the moment of recombination, when the first neutral atoms formed at rec ≈ 1089. Prior to recombination, photons and baryons were tightly coupled and oscillated as a single plasma until they decoupled at dec ≈ 1020. These oscillations, referred to as the baryon acoustic oscillations (BAO), are observed today as fluctuations of the CMB temperature (e.g. Planck Collaboration et al. 2020). Acoustic oscillations in the baryons' velocity at the moment of their decoupling lead to clumping in the distribution of baryons at later times, resulting in over and underdense regions. The initially tiny perturbations grew under the effect of gravity and are detected today in the distribution of galaxies on the largest cosmological scales (e.g. Alam et al. 2017).
Not only did the plasma oscillations lead to the BAO features in the post-recombination distribution of baryons, they also impacted the baryons' velocities (Sunyaev & Zeldovich 1970). Tseliakhovich & Hirata (2010) were the first to point out that the BAO pattern is imprinted in the magnitude of the relative velocity between baryons and dark matter, because the latter was not coupled to the primordial plasma at the time of recombination. At decoupling, the relative velocity had a root-mean-square (RMS) of ⟨ 2 bc ⟩ 1 2 ≈ 30 km s −1 , or ∼ 10 −4 with the speed of light. There is a vast separation of scales relevant to the relative velocity. Over scales smaller than a few ★ E-mail: luke.conaboy@nottingham.ac.uk Mpc (the coherence scale), the relative velocity is roughly constant; however, box sizes greater than a few hundred Mpc are required to properly sample the relative velocity (see Fig. 1 of this work and also fig. 1 in Tseliakhovich & Hirata 2010). At recombination the sound speed of the baryonic fluid dropped from being relativistic, ∼ / √ 3, to the thermal velocities of hydrogen atoms, ∼ 2 × 10 −5 , which is much smaller than the RMS value of bc . Therefore, on average, at decoupling baryons were travelling with supersonic velocities relative to the underlying potential wells generated by dark matter haloes (Tseliakhovich & Hirata 2010). The relative velocity remained supersonic all the way down to ∼ 15, sourcing shocks and entropy generation (Tseliakhovich & Hirata 2010;O'Leary & Mc-Quinn 2012). The amplitude of the velocity field decayed with time as (1+ ), and thus the effect weakened as the Universe expanded. For instance, the signature of bc in the low-power spectrum of BOSS galaxies (Yoo & Seljak 2013;Beutler et al. 2017) and the threepoint correlation function of BOSS CMASS galaxies was found to be negligible (Slepian & Eisenstein 2015;Slepian et al. 2018).
In the post-recombination Universe, growth of structure on large cosmological scales is generally described by linear perturbation theory, which follows the evolution of density and velocity fields to the leading order in perturbations. Despite being formally second-order contributions, terms that involve the supersonic relative velocity can actually be as large as the first order terms. Moreover, on scales below the coherence scale, bc is position-independent and the secondorder terms become effectively linear (Tseliakhovich & Hirata 2010). Using such a 'quasi-linear' approach, analytical methods were employed to explore implications of bc in the cosmological context. Supersonic relative velocities modulate the abundance of minihaloes and their gas content on the BAO scale (e.g. Tseliakhovich & Hirata 2010;Tseliakhovich et al. 2011;Fialkov et al. 2012;Ahn 2016;Ahn & Smith 2018), affecting fluctuations of the 21-cm signal of neutral hydrogen (e.g. Dalal et al. 2010;McQuinn & O'Leary 2012;Visbal et al. 2012;Fialkov et al. 2013;Cohen et al. 2016;Fialkov et al. 2018;Muñoz 2019;Cain et al. 2020;Muñoz et al. 2022;Long et al. 2022).
Finally, a hybrid approach was used to incorporate the non-linear effects into the large-scale cosmological picture by tiling regions of fixed bc together (e.g. Visbal et al. 2012;Fialkov et al. 2013). In such studies, the distribution of bc on scales larger than the 'pixel' size was generated from the corresponding density field using the continuity equation, while star formation in each 'pixel' was calibrated to numerical simulations (Maio et al. 2011;Stacy et al. 2011;Greif et al. 2011;Naoz et al. 2012Naoz et al. , 2013. This method was applied to the 21-cm signal of neutral hydrogen revealing enhanced BAO patterns (Visbal et al. 2012;Fialkov et al. 2013).
To fully capture the non-linear effect of bc in the cosmological context, it is necessary to properly include the velocity effect in the initial conditions of -body and hydrodynamical simulations. Such a task would require an accurate non-linear treatment of dark matter, baryons, and radiation, starting at rec and following the growth of structure all the way down to the lowest simulated redshift. This treatment is not possible due to the large dynamical range: the amplitude of density fluctuations at recombination is smaller than the precision of many commonly used integration schemes. Today, stateof-the-art numerical simulations are typically initialised at ini ∼ 200 with fields that do not incorporate the effect of the relative velocity.
Recently, Ahn & Smith (2018) introduced a new cosmological initial condition generator bccomics, based on a code that solves the quasi-linear equations between rec and ini for fixed values of largescale density, , and bc at decoupling (Ahn 2016). Next, the solver is applied to a larger cosmic volume divided in regions of fixed and bc . The code simulates the growth of small-scale structure inside density peaks and voids, by treating each patch of fixed and bc as a separate universe. They do this because the evolution equations from Tseliakhovich & Hirata (2010) are only valid at mean density, and Ahn (2016) found that structure formation was boosted (suppressed) in overdense (underdense) patches. They estimated the effect of bc and on the abundance of small haloes and, in some cases, reported qualitative disagreement with prior works.
Here we take an independent approach and develop a new initial conditions generator. Our approach differs from that of bccomics as, instead of treating small patches as separate universes, we follow small patches embedded in a larger cosmological volume through the use of zoom simulations. We compensate for the lacking effect of bc between rec and ini = 200 through the use of a bias factor, described in Section 3.3. After this compensation, the effect of bc from ini = 200 is naturally included by the simulation, through initialising the simulations with separate transfer functions for the baryon and dark matter velocities. Unlike Ahn (2016), we do not account for the largescale density environment, since our simulation focuses on a patch with very near to mean density and so we can safely ignore the impact on structure formation due to non-zero overdensity. Our methodology employs the widely used code music (Hahn & Abel 2011) to generate high-resolution 'zoom' initial conditions (Bertschinger 2001) in large cosmic volumes and by design is well-matched to AMR simulations. We demonstrate the performance by generating initial conditions in a 400 ℎ −1 Mpc box before extracting a 100 ℎ −1 Mpc subbox, which is used to run a simulation from ini = 200 to the final redshift of 11.2 with the AMR code ramses (Teyssier 2002). We explore the performance by computing the effects of bc on the number of haloes formed, baryon fraction, and star formation.
The paper is organised as follows: in Section 2, we recap the theoretical background and discuss why large simulation box sizes are needed; in Section 3, we discuss the simulation setup and our methodology for incorporating the effect of the bc through a scaledependent bias parameter ( , bc ); in Section 4, we present the results of a first demonstration of our methodology and discuss the findings in Sections 5 and 6. We present a comparison of our methodology to some previous works in Appendix A. Throughout this paper, we assume a flat ΛCDM cosmology consistent with the Planck 2018 results (Planck Collaboration et al. 2020) with parameters: Ω m = 0.314, Ω Λ = 0.686, Ω b = 0.049, = 0.965, 8 = 0.812 and ℎ = 0.673 1 .

THEORY
In this section, we briefly review the relevant theoretical background, directing the reader to Fialkov (2014) and Barkana (2016) for more comprehensive reviews.
In the non-linear regime, the evolution of the density and velocity of baryons and dark matter is governed by: where b and c are dimensionless perturbations in baryonic and dark matter densities respectively, b and c are the velocities of baryonic and dark matter respectively, is the scale factor, ≡ / is the Hubble parameter, Φ is the gravitational potential, is the baryonic pressure, and¯b and¯m are the average densities of baryons and total matter, respectively.
Following Tseliakhovich & Hirata (2010), we split the velocities into a coherent bulk motion, bc (of magnitude bc ), and random velocity perturbations, b and c , so that in the cold dark matter frame, the velocities can be written as b = bc + b and c = c .
Though they are second-order terms, components involving bc are large for typical values of bc at high redshifts. In addition, these terms become effectively first order on scales where bc is coherent. In this quasi-linear regime, perturbations in density ( b and c ) and velocities ( b and c ) evolve according to the following set of equations: where = −1 ∇ · is the velocity divergence in comoving coordinates (though contrary to Tseliakhovich & Hirata 2010, we work in the rest frame of the dark matter; see also the appendix of O'Leary & McQuinn 2012), is the fluctuation in the baryons' temperature, B is the Boltzmann constant, is the mean molecular weight, and p is the mass of the proton. Both the baryon and dark matter density parameters Ω b and Ω c are functions of (we drop the explicit dependence for clarity).
The importance of baryonic pressure for the growth of density modes was stressed by Naoz & Barkana (2005, 2007, which requires solving an extra equation to track fluctuations in the temperature . We follow Bovy & Dvorkin (2013) and Ahn (2016) in neglecting tracking fluctuations in photon density and temperature within the evolution equations, since they are subdominant at most of our scales and redshifts of interest. We then add the equation for the temperature fluctuations and¯= 2.726 K/ is the mean photon temperature, e ( ) is the electron fraction out of the total number density, of gas particles at time ,¯is the mean gas temperature,¯, 0 is the mean photon energy density at = 0, T is the Thomson scattering cross-section for an electron and e is the mass of an electron. Both e ( ) and are calculated using recfast++ (Seager et al. 1999;Chluba et al. 2010;Chluba & Thomas 2011). The initial conditions for are set as in Naoz & Barkana (2005) by requiring that ( − )/ = 0 at the initial redshift = 1000, where / is calculated from camb (Lewis et al. 2000). The above set of linearised equations can be solved using a publicly available code, e.g. cicsass (O'Leary & McQuinn 2012). We reimplement a python version of cicsass, which complements our methodology. In the current implementation, for simplicity, we ignore the directionality of bc when solving this set of equations, instead taking bc · = bc cos = bc (i.e. assuming bc is parallel to ). Future implementations will take the directionality of bc into account. Tseliakhovich & Hirata (2010) showed that most of the contributions to the variance of bc come from scales between 0.005 ℎ Mpc . RMS bc at = 200 as a function of box size, calculated by integrating the bc power spectrum Δ 2 bc (computed for the cosmological parameters listed in Section 1) from 2 / to 3365 ℎ Mpc −1 . The RMS bc converges for ≳ 400 ℎ −1 Mpc, since Δ 2 bc drops to zero at small-, and so box sizes of this or larger are needed to capture all the relevant scales. power spectrum of bc fluctuations from the fundamental mode of the box min = 2 / to infinity. The mean square bc in a box of size is given by where Δ 2 bc is the dimensionless power spectrum of the bc , taken from camb and, in theory, the upper limit max of the integral in equation (5) should be the maximum wavenumber of the box, dictated by the number of simulation elements. In practice, however, any upper limit max ≫ 0.5 ℎ Mpc −1 is sufficient, since the bc power spectrum drops off rapidly above this value. Fig. 1 shows the RMS bc , calculated as the square root of equation (5), where the oscillatory nature of Δ 2 bc at low-(cf. fig. 1 in Tseliakhovich & Hirata 2010) is clearly visible. From Fig. 1, we can see that even in a box size of 100 ℎ −1 Mpc, we do not capture all of the scales relevant to bc . The curve only begins to plateau around ∼ 400 ℎ −1 Mpc, so using a box size smaller than this means that we may miss out on some of the effect, for example by not sampling extreme values of bc . Simultaneously simulating this large-scale box and the very high-resolution zoom region needed to observe the effect would be computationally infeasible. In Section 3.2, we discuss our solution to this problem.

Simulations
We follow the evolution of dark matter, gas, and stars in the cosmological context using ramses 2 , which employs a second-order Godunov method to solve the equations of hydrodynamics. Gas states are computed at cell interfaces using the Harten-Lax-van Leer-contact Riemann solver, with a MinMod slope limiter. Dark matter and stars are modelled as a collisionless -body system, described by the Vlasov-Poisson equations. Grid refinement is performed whenever a cell contains more than eight high-resolution dark matter particles, or has the equivalent amount of baryonic mass scaled by Ω b /Ω m . We allow the AMR grid to refine from the coarsest level ℓ min = 8 to the finest level ℓ max = 23, but in practice grid hold-back within ramses means that the finest level reached is ℓ = 21, corresponding to a maximum comoving resolution of 47.7 ℎ −1 pc 3 .
Star formation is allowed whenever the gas density of a cell is greater than ★ = 1 cm −3 in units of the number density of hydrogen atoms and when the local overdensity is greater than 200 cr , where the latter condition prevents spurious star formation at extremely high redshift. We impose a polytropic temperature function with index ★ = 2 and 0 = 1050 K, which ensures that the Jeans length is always resolved by at least eight cells. We do not rigorously calibrate the star formation parameters to reproduce any stellar mass-halo mass relation, since we are interested only in the differences between simulations. Star particles, which represent a population of stars, form with a mass of 108.0 ℎ −1 M ⊙ . Supernova feedback is included using the kinetic feedback model of Dubois & Teyssier (2008), with a mass fraction SN = 0.1 and a metal yield of 0.1. We allow gas cooling and follow the advection of metals. We do not include molecular hydrogen in this simulation so, to attempt to compensate for this missing cooling channel, we initialise the zoom region with a metallicity of = 10 −3 Z ⊙ , where Z ⊙ = 0.02 in ramses.

Initial conditions
As described in Section 2, large box sizes of > ∼ 400 ℎ −1 Mpc are required in order to capture all of the scales pertaining to bc . By performing calibration runs, we found that very high resolution (a cell size of Δ < ∼ 2 ℎ −1 kpc) is needed in the ICs in order to properly resolve the effect. To this end, we employ 'zoom' initial conditions (ICs), generating density and velocity fields at ini = 200 first in a 400 ℎ −1 Mpc box using music (Hahn & Abel 2011). The ICs are refined from the base level ℓ min = 10 (1,024 3 ) up to ℓ = 18 (262,144 3 effective) in a cube of side length 543 ℎ −1 kpc at the finest level. Such extremely high resolution is required because bc suppresses structure formation on very small scales. We found that the resolution we used in this work is the minimum necessary to observe the effect of bc , and using lower resolution largely misses the effect. Since the zoom region is very small compared to the box size, we use extra padding between zoom levels, increasing the number of padding cells on each side for each dimension from the typical value of 4 to 32. We use transfer functions from camb (Lewis et al. 2000), which gives distinct density and velocity fields for the baryons and dark matter.
In order to make the simulation tractable, we extract a 100 ℎ −1 Mpc 3 Releasing higher levels of AMR grids at specified steps in scale factor ('grid hold-back') is a technique employed in ramses to ensure that the physical (as opposed to comoving) resolution remains roughly constant over an entire simulation, which is desirable for e.g. ensuring that cells at high redshift do not over-refine and end up containing too little mass to form stars. Snaith et al. (2018) performed a detailed study of the effects of grid holdback on simulation properties, finding, among other results, that the sudden release of high-resolution grids can lead to spikes in the star formation rate. Our simulation reaches its highest refinement level ℓ = 21 before any star formation occurs, and so is not impacted by these grid hold-back effects on star formation (ℓ min = 8, 256 3 ) base grid from the 400 ℎ −1 Mpc (ℓ min = 10, 1,024 3 ) box and use this as our coarsest level, meaning that the maximum refinement level in the zoom region also drops two levels from ℓ = 18 to ℓ = 16 (65,536 3 effective). In principle, this methodological choice could introduce some error around the edges of the box due to using periodic boundary conditions with non-periodic ICs, though in practice we expect the impact of this to be negligible since we are concerned with a sub-ℎ −1 Mpc region in the centre of the box. Additionally, these errors will be common between runs, so their effect will wash out when comparing between runs. Table 1 details the sets of ICs used in this work. We selected a region for zoom-in with bc,ini = 20.09 km s −1 at = 200, corresponding to bc,rec = 100.07 km s −1 , or ∼ 3.3 bc , at recombination. The no bc case is often used in cosmological simulations, for example when using transfer functions that do not have separate amplitudes for the baryon and dark matter velocity fields (in fact, it is the default behaviour for older versions of ramses, where the dark matter velocity field is used to initialise both the dark matter and baryon velocities). The bc -ini case is where the simulation is initialised using separate transfer functions for the baryon and dark matter velocity fields, such as by generating ICs using music with transfer functions from camb. In this case, bc is included from the start time of the simulation ini , but the effect of bc on density and velocity perturbations between recombination and ini is missed. In the final, and most realistic, case, bc -rec, we include the contributions from bc across all by computing a bias factor which is applied to the ICs. The methodology for computing the bias factor is detailed in Section 3.3.

Bias factor
Using transfer functions that have distinct amplitudes for the baryon and dark matter velocity fluctuations naturally yields the bc field at the start time of the simulation ini . First, we interpolate the dark matter particle velocities onto the same grid as the baryons, then take the difference of these two fields to calculate the magnitude as bc = | b − c |. A 0.39 ℎ −1 Mpc thick slice through the resultant bc field is shown in Fig. 2.
With the bc field in hand, we split our ICs into cubic patches, aiming for a patch extent of 0.5 ℎ −1 Mpc, though the actual extent depends upon how many patches can be fit in each level of the grafic files. The size of these patches is chosen to be smaller than the scale over which bc is coherent (Tseliakhovich & Hirata 2010). Within each patch, the average value of bc is calculated and used as bc in equations (2). The initial values for equations (2) are set using the transfer functions from camb at = 1000, and the equations are integrated from = 1000 to = 200 using the lsoda ordinary differential equation solver. Equations (2) are solved for the average patch value of bc and also for bc = 0 km s −1 , which yields power spectra for the baryon perturbations both with and without bc . We use these power spectra to calculate a 'bias' factor at ini = 200 that depends both upon scale and the magnitude of the relative velocity bc where the square root arises from ∝ 2 . In Fig. 3, we show the bias factor for the baryon and dark matter densities ( b and c ) and velocities ( b and c ), computed for the average bc in our zoom region. The strongest suppression is seen in the baryons and in particular the baryon density, while the dark matter is hardly affected. We do not expect the oscillatory features in ( , bc ) at the very small scales to have much, if any, impact since the power spectrum of fluctuations in the baryon density contrast begins to fall rapidly for ≳ 300 ℎ Mpc −1 , while for the velocity most of the power is at much larger scales. Ali-Haïmoud et al. (2014) also found oscillatory features in the small-scale baryon perturbations, and we have checked that we find similar oscillations for typical values of bc and also find that increasing the magnitude of bc increases the frequency of oscillations for the larger-scale (≳ 100 ℎ Mpc −1 ) modes too. For a detailed study into the origin of these small-scale oscillations, we defer the reader to Ali-Haïmoud et al. (2014). This factor is then convolved with the Fourier transform of the corresponding patch of baryon overdensitŷ to give individual patches of biased overdensityˆb, which are then stitched together to generate the bc -rec set of ICs. In this way, the bias factor compensates for the suppression of baryon perturbations  Figure 3. The bias factors ( bc ) at ini = 200 for the average bc in the zoom region, scaled back to its value at recombination bc,rec = 100.07 km s −1 , which is a ∼ 3.3 bc value. We show ( , bc ) for baryon and dark matter overdensities and peculiar velocities. These are the bias factors that are applied to the perturbations in the zoom region. Note how the perturbations in the baryon overdensity (dark blue, short-dashed) are strongly suppressed for > 40 ℎ Mpc −1 , while the perturbations in the dark matter overdensity (dark red, long-dashed) are largely unaffected, to the few per cent level. Note that we show ( , bc ) for the peculiar velocities for the full range of , but the velocity field is dominated by large-scale (small-) modes. Therefore, ( , bc ) will have little if any impact on the small scales (large-) of the velocity field. between = 1000 and ini that is missing if bc is included only from ini . We only modify the baryons, since as discussed earlier, they are much more strongly affected than the dark matter, as can be seen from Fig. 3.
We deal with the peculiar velocity field for the baryons in a similar way, by first converting the velocity divergence to peculiar velocities as b ( ) = −i b ( )/ 2 . Note again that we do not include directionality when solving the evolution equations, and therefore, the bias factor is applied to each direction of b equally. In reality, there would be preferential directions for the bias factor, depending on the direction of bc , but we defer that implementation to future work. Fig. 4 shows a 1.53 ℎ −1 kpc thick slice through the highest resolution level of the zoom ICs directly from music ('unmodified', left column) and after the bias factor ( , bc ) has been applied ('modified', right column). For b (top row), the unmodified ICs contain a lot of small-scale structure, which is almost totally washed out after applying ( , bc ). Most of what remains is in the form of lower amplitude, larger scale fluctuations. For b, (bottom rows) 4 , there is less small-scale structure to begin with, since the peculiar velocity fields are dominated by large scales. The effect of ( , bc ) on b, is therefore much less striking than on the b , with the main effect being smoothing and a slight reduction in amplitude.

Haloes
After the ICs have been correctly initialised with bc , we can characterise the effect of bc on structure formation, principally by exploring  (7), right column) baryon overdensity (top row), peculiar velocity in the -direction (second row), -direction (third row) and -direction (bottom row) in the high-resolution zoom region, of side length 543 ℎ −1 kpc. Each pixel corresponds to a cell width of 1.53 ℎ −1 kpc and the slice has a thickness of one cell width. The effect of applying ( , bc ), defined in equation (6), can be clearly seen in the baryon overdensity, in that it washes out the small-scale fluctuations. The effect is less pronounced in the peculiar velocities, which are dominated by large-scale modes. how haloes are affected. Haloes are identified using ahf (Gill et al. 2004;Knollmann & Knebe 2009), which supports multi-resolution datasets and calculates baryonic properties of haloes. In order to use ahf with a ramses dataset, we use the supplied ramses2gadget tool to convert leaf cells of the AMR hierarchy into gas pseudoparticles placed at the centre of each cell. This conversion is a standard recipe in the analysis of AMR datasets (e.g. Wu et al. 2015) and allows us to conveniently assign gas to haloes. We ignore the internal energy of the gas and do not allow ahf to perform any unbinding, as this has been shown to remove most of the gas from subhaloes in a manner that is dependent upon the choice of halo finder (Knebe et al. 2013). We define the halo overdensity with respect to the critical density cr , such that the average density inside the halo is 200 cr . Guided by the resolution study of Naoz & Barkana (2007), we perform our analysis on haloes that have ≥ 500 particles, as they found that this is the level at which the baryon fraction is resolved with a scatter of ∼ 20 per cent when compared to higher resolution simulations.
We only use haloes comprised entirely of high-resolution dark matter particles, since contaminant low-resolution particles can disrupt the dynamics of haloes (Oñorbe et al. 2014). Due to the small size of the zoom region, it is often the case that haloes which were initially solely composed of high-resolution particles can become contaminated by low-resolution particles as the simulation progresses. If contaminated haloes are removed at each timestep, then haloes can 'disappear' if they become contaminated between one timestep and the next. To counter this effect, we generate merger trees using consistent-trees (Behroozi et al. 2013) and only keep haloes that have never been contaminated at any point in the simulation. Then, when exploring the effect on halo properties (such as when looking at their gas content), we go one step further and match haloes between the simulations. We do this using ahf's MergerTree tool to correlate the particle IDs between each run. This allows us to isolate the effect on identical haloes, as opposed to also capturing the impact on the global population of haloes. The difficulties in exploring the impact on the global population of haloes (e.g. through the cumulative number of haloes) are explored in Section 4.1.

Halo abundances
Throughout this section, we calculate the cumulative number of haloes (> ) as the number of haloes with a total mass greater than and estimate the Poisson uncertainty on (> ) for each case as √︁ (> ). To begin with, we ignore the conditions described in Section 3.4, instead considering the cumulative number of all haloes formed in the simulation all (> ), shown in Fig. 5. Aside from a slight suppression in all (> ) for the bc -rec case below 10 6 ℎ −1 M ⊙ at = 14.2, the bc -ini and bc -rec cases are consistent both with each other and with the no bc case at the 1 level. Analysing the haloes in this fashion (i.e. by performing no cleaning on the catalogue) gives a picture of the global impact of bc on halo abundance, however this picture is inaccurate because a significant fraction of the haloes are contaminated by, or formed entirely of, lower-resolution particles, affecting the accuracy of their properties (see Section 3.4 for a discussion of contamination).
If we now apply the conditions described in Section 3.4, namely that we do not include in our analysis any haloes that are contaminated at any point in the simulation, we are left with a reduced catalogue. In Fig. 6, we show the cumulative number of haloes in this cleaned catalogue clean (> ), which shows a much more striking suppression of clean (> ) for the bc -rec case and an apparent increase in the abundance of low-mass ( < 10 6 ℎ −1 M ⊙ ) haloes for the bc -ini case at = 11.2 (though still consistent with no difference to the no bc case at the 1 level).
To highlight the effect of the cleaning procedure we also show the raw cumulative number of haloes at = 14.2 and = 11.2 in Fig. 7. The impact of the cleaning process described in Section 3.4 can be clearly seen in the bottom panel of Fig. 7, which shows the ratio of the cleaned catalogue to the full catalogue clean (> )/ all (> ). Comparing the different runs, it can be clearly seen that a different fraction of haloes is removed between each run at each of the times shown. that have been selected by the process detailed in Section 3.4. Selecting never-contaminated haloes in this fashion leads to the removal of an unequal and inconsistent number of haloes between each run-this leads to the odd behaviour of the bc -ini curve, which jumps from a suppression to a boost in the number of haloes with < 2 × 10 6 ℎ −1 M ⊙ between = 12.6 and 11.2.

Baryon fraction
We allow star formation in these runs, so the total baryon fraction of a given halo is defined as where g is the gas, ★ the stellar, and d the dark matter mass in each halo. We upweight the best resolved (i.e. most massive) haloes by calculating the mass-weighted average baryon fraction as and the associated mass-weighted standard deviation as where the sum is over all haloes that satisfy the conditions in Section 3.4. Fig. 8 shows ⟨ b ⟩ and associated 1 errorbars as function of . We show each where ≥ 30 haloes have formed that satisfy We show (> ) at = 14.2 (thin) and 11.2 (thick). From the top panel, it can be seen that all (> ) is very similar between all runs at each shown, except for a slight suppression in the abundance of haloes in the bc -rec case at = 14.2. The huge reduction in the number of haloes after cleaning the catalogue is abundantly clear, as is the difference in fraction of haloes removed between each run. the criteria in Section 3.4, starting from = 13.6 where we are able to match 38 haloes between the three cases. The gas fraction is suppressed at all for the bc -ini and bc -rec cases compared to the no bc case. At earlier , the suppression is stronger, though even by the final snapshot at = 11.2, ⟨ b ⟩ for both the bc -rec and bc -ini cases are not within 1 of the no bc case. Notably, at all , ⟨ b ⟩ in bc -ini and bc -rec cases are almost indistinguishable from, and certainly consistent with, one another.

Star formation
In Fig. 9, we show the cumulative ★ formed in the simulation, not accounting for mass loss due to supernovae, and the corresponding number of stellar particles ★ , which each have a mass of 108.0 ℎ −1 M ⊙ . In each case, all of the star particles in the simulation formed inside a single halo. In total, 29 star particles formed by = 11.2 in the no bc case, 10 in the bc -ini case, and 7 in 11.5 12.0 12.5 13.0 13.5 z  (10). We find a suppression in both the bc -ini and bc -rec cases compared to the no bc run, though between the two runs with bc there is little difference. . Cumulative stellar mass ★ formed as function of redshift . Also shown is the corresponding number of stellar particles ★ . Despite the small number of stars formed in this simulation, there is a clear hierarchy in how many stars each forms, with the no bc run forming the most, followed by bc -ini and then bc -rec. There is also a clear delay in the onset of star formation which follows a similar trend, with the bc -rec run starting to form stars the latest. the bc -rec case. This hierarchy persists across all , with more star particles having formed in the no bc case than in the bc -ini case and fewer still in the bc -rec case.
In the following, all times quoted in Myr are measured relative to the Big Bang. To quantify how bc impacts star formation in our simulations, we can calculate the delay in formation time of the th star for each run with bc compared to the no bc run , no bc as Δ = − , no bc . For the bc -ini case, the smallest (largest) delay of 8.6 Myr (35.5 Myr) occurs for the formation of the fourth (second) star particle, while for the bc -rec case we find a delay of 23.4 Myr (49.4 Myr) for the sixth (second) star particle to form.
Star formation in ramses (like in other simulation codes) is a stochastic process, hence the large fluctuations in delay times. We can reduce the impact this stochasticity has on our results by averaging over delay times, as opposed to considering the delay times for individual star particles. We find a mean delay time of ⟨Δ ⟩ = 19.4 Myr and 34.9 Myr for the bc -ini and bc -rec cases, respectively.

DISCUSSION
To understand the perceived uptick in bc -ini haloes in Fig. 6, we need to explore the raw cumulative number of haloes, shown in Fig. 7 at = 14.2 and = 11.2. The key point is that while only keeping haloes that are never contaminated ensures that haloes do not disappear between timesteps of the same run, it also means that different numbers of haloes will be removed between different runs at any given timestep. This discrepancy in the number of haloes removed is responsible for the apparent low-mass boost in the bc -ini case, because over the mass range 3×10 5 < /ℎ −1 M ⊙ < 10 6 fewer haloes are removed in the bc -ini case than in the no bc case. Since the cumulative numbers of haloes are already very similar, removing a smaller number of haloes in the bc -ini case manifests itself as a boost relative to the no bc case, when in fact it is simply an artefact of the cleaning procedure.
Given the difficulties in extracting a sample of haloes that is both free from contamination and comparable between runs, we will not draw any quantitative conclusions on the impact of bc on global properties like the number of haloes formed. We retain this discussion of the well-known impact of contamination and, crucially, the importance of verifying any mitigation techniques as it may prove helpful for other works employing zoom simulations.
Including bc also significantly affects the baryon fraction b , where we see that the mass-weighted baryon fraction ⟨ b ⟩ is suppressed at all redshifts in both cases, with the suppression stronger at higher redshift. Even by = 11.2, the bc -ini and bc -rec cases are still not in agreement with the no bc case, though the difference between the two populations has decreased. Again, this is likely due to the decay in the magnitude of bc , which allows the haloes to accrete more gas. Interestingly, ⟨ b ⟩ is almost indistinguishable between the bc -ini and bc -rec cases, suggesting that including bc from ini is sufficient to observe its impact on halo baryon fraction, though a larger sample is required to confirm this. The suppression in ⟨ b ⟩ when bc is included is in qualitative agreement with previous studies.
This decrement in baryon fraction for the bc -ini and bc -rec cases is reflected in the cumulative stellar mass formed, as fewer star particles formed in both cases than in the no bc case. Not only do they form fewer star particles, they also start forming star particles later since the effect of bc is to wash out the peaks (and troughs) in the baryon density contrast, meaning that it takes longer for gas to reach the densities required for star formation. The mean delay for forming star particles is 19.4 Myr for the bc -ini case and 34.9 Myr for the bc -rec case. From Schaerer (2002), we find that these delays are all of the order of the lifetime of a 9 M ⊙ first-generation Population (Pop) III star, which has a lifetime of 20.02 Myr (table 3 in Schaerer 2002). More massive Pop III stars have even shorter lifetimes, for example a 120 M ⊙ Pop III star lives for only 2.52 Myr. Pop III stars form from initially pristine gas, and their death pollutes their immediate surroundings with metals, introducing new cooling channels into the high-redshift Universe. Any delay in this introduction of metals will delay the transition between Pop III to Pop II (i.e. from metalenriched gas), which can, for example, affect the 21 cm signal (Magg et al. 2022). In our case, though we do not form Pop III stars, chemical enrichment is still vitally important for star formation to get properly underway, particularly as all of the star particles form in the same halo.
Despite there being almost no difference in ⟨ b ⟩ between the bc -ini and bc -rec cases at most redshifts, there is a clear hierarchy in the amount of stars formed-no bc forms the most, bc -ini forms fewer, and bc -rec forms the least-albeit on the order of a few star particles. This effect is expected, since the bias factor washes out baryonic density peaks, and there are slightly more haloes (i.e. star formation locations) present in the bc -ini case than in the bc -rec case.

CONCLUSIONS
We have performed the first cosmological zoom simulations that self-consistently sample the relative baryon-dark matter velocity bc from a large 400 ℎ −1 Mpc box. This relative velocity naturally arises when simulations are initialised using transfer functions that have separate amplitudes for the baryon and dark matter velocities, and we have shown that a box roughly as large as this is required to properly sample all of the scales associated with the relative velocity. However, solely initialising simulations in this manner misses out on the effect of the relative velocities from = 1000 to the start time of the simulation, ini . We developed a methodology that compensates for the effect of bc on baryon density and velocity perturbations by computing a 'bias' factor ( , bc ), which is convolved with the ICs. We verified that our methodology performs as expected by comparing to previous works (see Appendix A).
As a first demonstration of our methodology, we applied it to an extremely high-resolution zoom region in a 100 ℎ −1 Mpc subbox, extracted from the main 400 ℎ −1 Mpc box. The zoom region is centred on the region with the largest relative velocity in the 400 ℎ −1 Mpc box, which has an RMS value of bc = 100.07 km s −1 at = 1000, corresponding to ∼ 3.3 bc . We find qualitative agreement with previous works, namely a reduction in the halo baryon fraction, a delay in the onset of star formation at high redshift, and a suppression of the final stellar mass. The strength of the effect decreases with redshift, but the two simulations still exhibit some differences by = 11.2. We find that the delay in the onset of star formation is of the order of the lifetime of a ∼ 9 M ⊙ Pop III star. We also test the effect of incorporating the bias factor by running a simulation that includes the relative velocity from the start time of the simulation only. In this case, we find that more stars are formed when compared to the simulation that includes the bias factor, but there is almost no change in the average baryon fraction, except at the earliest redshift. Due to the small size of our zoom region, the vast majority of haloes in our simulation are contaminated with low-resolution particles, and we are thus unable to draw any robust conclusions regarding the halo mass function.
Our code for producing these compensated ICs is publicly available 5 , and we hope will be of use for studying this effect in the full cosmological context.

ACKNOWLEDGEMENTS
We thank the anonymous referee for constructive feedback which improved this paper. LC thanks Antony Lewis and Joakim Rosdahl for helpful discussions. LC was supported by an STFC studentship and acknowledges the indirect support of a Royal  (Virtanen et al. 2020); yt (Turk et al. 2011); ytree (Smith & Lang 2019); hmf (Murray et al. 2013) and cmasher (van der Velden 2020).

DATA AVAILABILITY
The code for computing and applying the bias factor ( , bc ) is available at https://github.com/lconaboy/drft. The data underlying this article will be shared on reasonable request to the corresponding author.
Beutler F., Seljak U., Vlah Z., 2017, Monthly Notices of the Royal Astronom-a cosmology consistent with Komatsu et al. (2009) with parameters: Ω m = 0.28, Ω Λ = 0.72, Ω b = 0.046, ℎ = 0.7 7 and a boosted 8 = 1.4. We adopt the boosted 8 = 1.4 as used in the original studies. We initialise the simulations at ini = 200 8 both with and without a relative velocity of 1.7 bc = 10 km s −1 at ini . We also compute and apply the bias factor ( , bc ) to the baryonic component of the ICs, while the Naoz et al. (2012) simulations are initialised by computing transfer functions which explicitly include the effect of the bc . Following Naoz et al. (2012), we set both of the velocity fields to that of the dark matter and apply the bc to the -component of the baryon velocity as Including the bc in this manner is justified by the box size being much smaller than the coherence scale of the bc . We do not allow star formation in these runs. Haloes are identified as described in Section 3.4, noting that although Naoz et al. (2012) define their halo overdensity with respect to the background matter density¯m, we are working in the period of matter domination so the difference will be small. While the final halo masses in Naoz et al. (2012) are defined using a spherical overdensity method, the initial halo finding is done using a friend-of-friends method, which could be susceptible to 'overlinking' (e.g. Davis et al. 1985), whereby disparate groups of particles are spuriously connected by a diffuse particle bridge. Their redefinition of the halo mass using the spherical overdensity method will go some way to alleviating the impact (if any) of overlinking.
To calculate the effect of the bc , we compare to simulations without bc , where the velocity field of the baryons is equal to that of the dark matter. In order to quantify this effect, we calculate the fractional difference of a quantity as First, we look at the effect on the cumulative halo mass function (> ), as in Naoz et al. (2012). Fig. A1 shows the decrement in (> ) for the case with bc compared to the case without, both for our simulations and for the Naoz et al. (2012) run. We see qualitatively similar behaviour, observing a decrement between ∼ 0 per cent and 50 per cent at all redshifts shown and for almost all masses.
However, the overall shape of our Δ is slightly different to Naoz et al. (2012); we match well below ∼ 3 × 10 5 ℎ −1 M ⊙ but show more relative suppression above this mass. This discrepancy is due, at least in part, to the different simulation codes used and the different white noise fields in the initial conditions. One further significant source of difference is the cosmologies used. Fig. A2 shows the difference expected at = 15 by comparing the analytic Watson et al. (2013) (> ) mass functions (magenta solid). From this, we would expect the Naoz et al. (2012) simulation to have ∼ 9 per cent more haloes with > 3 × 10 5 ℎ −1 M ⊙ . This increase in the number of haloes increases with mass and for > 1 × 10 7 ℎ −1 M ⊙ , we would expect ∼ 11 per cent more haloes in the Naoz et al. for the bc,ini = 10 km s −1 case, calculated with equation (A2). We show Δ at = 25 (top), 19 (centre), and = 15 (bottom) for our simulations (red solid) and for the Naoz et al. (2012) work (grey dashed). The shaded regions indicate the 1 Poisson uncertainty on the ratio. At most masses we are consistent with Naoz et al. (2012), though there is some disagreement for < 6 × 10 4 ℎ −1 M ⊙ , which could arise because of the friends-of-friends algorithm used to initially find haloes in Naoz et al. (2012). at all masses. At higher masses, Δ diverges as the absolute number of haloes becomes small. As discussed earlier, Naoz et al. (2012) use the friends-of-friends groups as a starting point, which could be a source of the high-mass discrepancy in Fig. A2.
Next, we turn our attention to the gas fraction of haloes, as studied in Naoz et al. (2013). Since we do not include star formation in these runs, the baryon fraction is simply the halo gas mass divided by the total halo mass b = g g + d .
(A3) Fig. A3 shows the binned gas fractions (top panel) for our and for the Naoz et al. (2013) simulations, each normalised to the cosmic mean Ω b /Ω m for the appropriate cosmology, and the decrement (bottom panel) as defined in equation (A2). We take the midpoint of the mass bin to be the mean of all the mass values in that bin. The binned gas fractions for Naoz et al. (2013) are slightly higher than in this work, though they exhibit roughly the same mass dependence. from the simulations.For ≲ 4 × 10 6 ℎ −1 M ⊙ the difference in (> ) between the simulations is mostly consistent with the expected difference due to cosmology and halo mass definitions (i.e. the difference between the Watson et al. 2013 curves). Potential sources of the high-mass ( ≳ 4 × 10 6 ℎ −1 M ⊙ ) discrepancy are discussed in the text.
The agreement between the two simulations for the decrement is striking -they have an extremely similar mass dependence. There is some difference in the binned baryon fractions, in particular we find slightly more suppression at lower masses. This is likely due to differences in code used since, as mentioned previously, Naoz et al. (2013) used gadget2 (Springel 2005), where we use ramses. There are well documented differences between Lagrangian (e.g. SPH) and Eulerian (e.g. AMR) codes (e.g. Agertz et al. 2007), and indeed it has been shown that numerical diffusion due to baryon-grid relative velocities can artificially smooth densities in Eulerian codes (Pontzen et al. 2020). In any case, we are not interested in comparing the merits of different codes, so by calculating the difference between the runs with and without bc , we can remove artefacts due to the choice of code.  Figure A3. Binned baryon fraction b (top panels) and relative fractional difference Δ b between each run with bc to the run without, calculated with equation (A2) (bottom panels). We show data from our simulations (red triangles, points and plusses) and from the Naoz et al. (2013) work (grey squares, stars and crosses). The panels show = 25 (left), 19 (centre) and 15 (right). The errorbars indicate the 1 standard deviation in each mass bin (top panels) and the combined 1 uncertainty (bottom panels). We find excellent agreement in the relative difference between our work and Naoz et al. (2013), and broad agreement in the absolute value of b . Some slight disagreement in the absolute value could be due to the different choice of simulation methodology, since we employ an Eulerian code where Naoz et al. (2013) use a Lagrangian code.