Relative baryon-dark matter velocities in cosmological zoom simulations

ABSTRACT

1 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 z_rec ≈ 1089. Prior to recombination, photons and baryons were tightly coupled and oscillated as a single plasma until they decoupled at z_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 VI 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 |$\langle v_{\rm bc}^2\rangle ^{\frac{1}{2}}\approx 30~\mathrm{km} \,\mathrm{s}^{-1}$|⁠, or |$\sim 10^{-4}\, c$| with c the speed of light. There is a vast separation of scales relevant to the relative velocity. Over scales smaller than a few 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, |$\sim c/\sqrt{3}$|⁠, to the thermal velocities of hydrogen atoms, ∼2 × 10⁻⁵ c, which is much smaller than the RMS value of v_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 z ∼ 15, sourcing shocks and entropy generation (Tseliakhovich & Hirata 2010; O’Leary & McQuinn 2012). The amplitude of the velocity field decayed with time as (1 + z), and thus the effect weakened as the Universe expanded. For instance, the signature of v_bc in the low-z power spectrum of BOSS galaxies (Yoo & Seljak 2013; Beutler, Seljak & Vlah 2017) and the three-point correlation function of BOSS CMASS galaxies was found to be negligible (Slepian & Eisenstein 2015; Slepian et al. 2018).

Figure 1.

RMS v_bc at z = 200 as a function of box size, calculated by integrating the v_bc power spectrum |$\Delta _{{v_{\rm bc}}}^2$| (computed for the cosmological parameters listed in Section 1) from 2π/L to 3365 h Mpc⁻¹. The RMS v_bc converges for L ≳ 400 h⁻¹ Mpc, since |$\Delta ^2_{v_{\rm bc}}$| drops to zero at small-k, and so box sizes of this or larger are needed to capture all the relevant scales.

Open in new tabDownload slide

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, v_bc is position-independent and the second-order terms become effectively linear (Tseliakhovich & Hirata 2010). Using such a ‘quasi-linear’ approach, analytical methods were employed to explore implications of v_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, Barkana & Hirata 2011; Fialkov et al. 2012; Ahn 2016; Ahn & Smith 2018), affecting fluctuations of the 21-cm signal of neutral hydrogen (e.g. Dalal, Pen & Seljak 2010; McQuinn & O’Leary 2012; Visbal et al. 2012; Fialkov et al. 2013; Cohen, Fialkov & Barkana 2016; Fialkov, Barkana & Cohen 2018; Muñoz 2019; Cain et al. 2020; Long, Givans & Hirata 2022; Muñoz et al. 2022).

Numerical simulations were used to explore non-linear effects of v_bc on scales well-below its coherence scale. Such simulations typically employ boxes of several comoving Mpc or less and assume a position-independent velocity field. These simulations demonstrated that v_bc suppresses formation of small dark matter haloes (Naoz, Yoshida & Gnedin 2012, 2013; O’Leary & McQuinn 2012), induces shocks (O’Leary & McQuinn 2012), affects the formation of first stars (e.g. Greif et al. 2011; Maio, Koopmans & Ciardi 2011; Stacy, Bromm & Loeb 2011; Schauer et al. 2019, 2021) and black holes (e.g. Hirano et al. 2017; Schauer et al. 2017), and may even influence shaping globular clusters (Naoz & Narayan 2014; Chiou et al. 2018, 2019, 2021; Druschke et al. 2020; Lake et al. 2021).

Finally, a hybrid approach was used to incorporate the non-linear effects into the large-scale cosmological picture by tiling regions of fixed v_bc together (e.g. Visbal et al. 2012; Fialkov et al. 2013). In such studies, the distribution of v_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 (Greif et al. 2011; Maio et al. 2011; Stacy et al. 2011; Naoz et al. 2012, 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 v_bc in the cosmological context, it is necessary to properly include the velocity effect in the initial conditions (ICs) of N-body and hydrodynamical simulations. Such a task would require an accurate non-linear treatment of dark matter, baryons, and radiation, starting at z_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, state-of-the-art numerical simulations are typically initialized at z_ini ∼ 200 with fields that do not incorporate the effect of the relative velocity.

Recently, Ahn & Smith (2018) introduced a new cosmological IC generator bccomics, based on a code that solves the quasi-linear equations between z_rec and z_ini for fixed values of large-scale density, δ, and v_bc at decoupling (Ahn 2016). Next, the solver is applied to a larger cosmic volume divided in regions of fixed δ and v_bc. The code simulates the growth of small-scale structure inside density peaks and voids, by treating each patch of fixed δ and v_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 v_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 ICs 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 v_bc between z_rec and z_ini = 200 through the use of a bias factor, described in Section 3.3. After this compensation, the effect of v_bc from z_ini = 200 is naturally included by the simulation, through initializing the simulations with separate transfer functions for the baryon and dark matter velocities. Unlike Ahn (2016), we do not account for the large-scale 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’ ICs (Bertschinger 2001) in large cosmic volumes and by design is well-matched to AMR simulations. We demonstrate the performance by generating ICs in a 400 h⁻¹ Mpc box before extracting a 100 h⁻¹ Mpc subbox, which is used to run a simulation from z_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 v_bc on the number of haloes formed, baryon fraction, and star formation.

The paper is organized 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 v_bc through a scale-dependent bias parameter b(k, v_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 VI 2020) with parameters: Ω_m = 0.314, Ω_Λ = 0.686, Ω_b = 0.049, n_s = 0.965, σ₈ = 0.812, and h = 0.673.¹

2 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, |$\boldsymbol { v}_{\rm b}$| and |$\boldsymbol { v}_{\rm c}$| are the velocities of baryonic and dark matter, respectively, a is the scale factor, |$H \equiv \dot{a}/a$| is the Hubble parameter, Φ is the gravitational potential, p is the baryonic pressure, and |$\bar{\rho }_{\rm b}$| and |$\bar{\rho }_{\rm m}$| are the average densities of baryons and total matter, respectively.

Following Tseliakhovich & Hirata (2010), we split the velocities into a coherent bulk motion, |$\boldsymbol {v}_{\rm bc}$| (of magnitude v_bc), and random velocity perturbations, |$\boldsymbol {u}_{\rm b}$| and |$\boldsymbol {u}_{\rm c}$|⁠, so that in the cold dark matter frame, the velocities can be written as |$\boldsymbol {v}_{\rm b} = \boldsymbol {v}_{\rm bc}+ \boldsymbol {u}_{\rm b}$| and |$\boldsymbol {v}_{\rm c} = \boldsymbol {u}_{\rm c}$|⁠.

Though they are second-order terms, components involving |$\boldsymbol {v}_{\rm bc}$| are large for typical values of v_bc at high redshifts. In addition, these terms become effectively first order on scales, where v_bc is coherent. In this quasi-linear regime, perturbations in density (δ_b and δ_c) and velocities (⁠|$\boldsymbol {u}_{\rm b}$| and |$\boldsymbol {u}_{\rm c}$|⁠) evolve according to the following set of equations:

where |$\theta _i = a^{-1}\nabla \cdot {\boldsymbol {u}_i}$| 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), δ_T is the fluctuation in the baryons’ temperature, k_B is the Boltzmann constant, μ is the mean molecular weight, and m_p is the mass of the proton. Both the baryon and dark matter density parameters Ω_b and Ω_c are functions of t (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 δ_T. We follow Bovy & Dvorkin (2013) and Ahn (2016) in neglecting tracking fluctuations in photon density and temperature within the evolution equations, since they are sub-dominant at most of our scales and redshifts of interest. We then add the equation for the temperature fluctuations

to equations (2), where

and |$\bar{T}_\gamma =2.726\ {\mathrm{K}}/a$| is the mean photon temperature, x_e(t) is the electron fraction out of the total number density, of gas particles at time t, |$\bar{T}$| is the mean gas temperature, |$\bar{\rho }_{\gamma ,0}$| is the mean photon energy density at z = 0, σ_T is the Thomson scattering cross-section for an electron and m_e is the mass of an electron. Both x_e(t) and |$\bar{T}$| are calculated using recfast+ + (Seager, Sasselov & Scott 1999; Chluba, Vasil & Dursi 2010; Chluba & Thomas 2011). The ICs for δ_T are set as in Naoz & Barkana (2005) by requiring that |$\partial (\delta _T-\delta _{T_\gamma })/\partial t = 0$| at the initial redshift z = 1000, where |$\partial \delta _{T_\gamma }/\partial t$| is calculated from camb (Lewis, Challinor & Lasenby 2000). The above set of linearized 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 |$\boldsymbol {v}_{\rm bc}$| when solving this set of equations, instead taking |$\boldsymbol {v}_{\rm bc}\cdot \boldsymbol { k}=v_{\rm bc}k~\mathrm{cos}\theta =v_{\rm bc}k$| (i.e. assuming |$\boldsymbol {v}_{\rm bc}$| is parallel to |$\boldsymbol {k}$|⁠). Future implementations will take the directionality of |$\boldsymbol {v}_{\rm bc}$| into account.

Tseliakhovich & Hirata (2010) showed that most of the contributions to the variance of v_bc come from scales between 0.005 and 0.5 h Mpc⁻¹. In a similar fashion to Pontzen et al. (2020), we can compute the RMS v_bc inside a box of size L by integrating the power spectrum of v_bc fluctuations from the fundamental mode of the box k_min = 2π/L to infinity. The mean square v_bc in a box of size L is given by

where |$\Delta _{{v_{\rm bc}}}^2$| is the dimensionless power spectrum of the v_bc, taken from camb and, in theory, the upper limit k_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 k_max ≫ 0.5 h Mpc⁻¹ is sufficient, since the v_bc power spectrum drops off rapidly above this value. Fig. 1 shows the RMS v_bc, calculated as the square root of equation (5), where the oscillatory nature of |$\Delta _{v_{\rm bc}}^2$| at low-k (cf. fig. 1 in Tseliakhovich & Hirata 2010) is clearly visible. From Fig. 1, we can see that even in a box size of 100 h⁻¹ Mpc, we do not capture all of the scales relevant to v_bc. The curve only begins to plateau around ∼400 h⁻¹ 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 v_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.

3 METHODS

3.1 Simulations

We follow the evolution of dark matter, gas, and stars in the cosmological context using ramses,² 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 N-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 h⁻¹ pc.³

Star formation is allowed whenever the gas density of a cell is greater than n_⋆ = 1 cm⁻³ 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 g_⋆ = 2 and |$T_0=1050\, \mathrm{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 h⁻¹ 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 initialize the zoom region with a metallicity of Z = 10⁻³ Z_⊙, where Z_⊙ = 0.02 in ramses.

3.2 ICs

As described in Section 2, large box sizes of |$\gtrsim 400~h^{-1} \mathrm{Mpc}$| are required in order to capture all of the scales pertaining to v_bc. By performing calibration runs, we found that very high resolution (a cell size of |$\Delta x \lesssim 2~h^{-1} \mathrm{kpc}$|⁠) is needed in the ICs in order to properly resolve the effect. To this end, we employ ‘zoom’ ICs, generating density and velocity fields at z_ini = 200 first in a 400 h⁻¹ Mpc box using music (Hahn & Abel 2011). The ICs are refined from the base level ℓ_min = 10 (1024³) up to ℓ = 18 (262 144³ effective) in a cube of side length 543 h⁻¹ kpc at the finest level. Such extremely high resolution is required because v_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 v_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 h⁻¹ Mpc (ℓ_min = 8256³) base grid from the 400 h⁻¹ Mpc (ℓ_min = 10, 1024³) 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³ 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-h⁻¹ 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 v_{bc, ini} = 20.09 km s⁻¹ at z = 200, corresponding to v_{bc, rec} = 100.07 km s⁻¹, or |$\sim 3.3\sigma _{v_{\rm bc}}$|⁠, at recombination. The no v_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 initialize both the dark matter and baryon velocities). The v_bc–ini case is where the simulation is initialized 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, v_bc is included from the start time of the simulation z_ini, but the effect of v_bc on density and velocity perturbations between recombination and z_ini is missed. In the final, and most realistic, case, v_bc–rec, we include the contributions from v_bc across all z by computing a bias factor, which is applied to the ICs. The methodology for computing the bias factor is detailed in Section 3.3.

3.3 Bias factor

Using transfer functions that have distinct amplitudes for the baryon and dark matter velocity fluctuations naturally yields the v_bc field at the start time of the simulation z_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 |$v_{\rm bc}= \left| \boldsymbol {v}_{\rm b}-\boldsymbol {v}_{\rm c} \right|$|⁠. A 0.39 h⁻¹ Mpc thick slice through the resultant v_bc field is shown in Fig. 2.

Figure 2.

A slice through the large-scale |$v_{\rm bc}=|\boldsymbol {v}_{\rm b}-\boldsymbol {v}_{\rm c}|$| field in the full 400 h⁻¹ Mpc ICs at z_ini = 200. Each pixel corresponds to a cell width of 0.39 h⁻¹ Mpc and the slice has a thickness of one cell width. Also shown (white square) is the position of the extracted 100 h⁻¹ Mpc subbox, centred on the peak v_bc in the ICs. The zoom region is located in the centre of the subbox. The colour shows the magnitude of the v_bc at z_ini = 200, where light pink is high v_bc and dark blue is low v_bc.

Open in new tabDownload slide

With the v_bc field in hand, we split our ICs into cubic patches, aiming for a patch extent of 0.5 h⁻¹ 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 v_bc is coherent (Tseliakhovich & Hirata 2010). Within each patch, the average value of v_bc is calculated and used as v_bc in equations (2). The initial values for equations (2) are set using the transfer functions from camb at z = 1000, and the equations are integrated from z = 1000 to z = 200 using the lsoda ordinary differential equation solver. Equations (2) are solved for the average patch value of v_bc and also for v_bc = 0 km s⁻¹, which yields power spectra for the baryon perturbations both with and without v_bc. We use these power spectra to calculate a ‘bias’ factor at z_ini = 200 that depends both upon scale k and the magnitude of the relative velocity v_bc

where the square root arises from P ∝ |δ²|. In Fig. 3, we show the bias factor for the baryon and dark matter densities (δ_b and δ_c) and velocities (v_b and v_c), computed for the average v_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 b(k, v_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 k ≳ 300 h Mpc⁻¹, while for the velocity most of the power is at much larger scales. Ali-Haïmoud, Meerburg & Yuan (2014) also found oscillatory features in the small-scale baryon perturbations, and we have checked that we find similar oscillations for typical values of v_bc and also find that increasing the magnitude of v_bc increases the frequency of oscillations for the larger-scale (≳ 100 h Mpc⁻¹) 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 overdensity

to give individual patches of biased overdensity |$\hat{\delta }_{\rm b}$|⁠, which are then stitched together to generate the v_bc–rec set of ICs. In this way, the bias factor compensates for the suppression of baryon perturbations between z = 1000 and z_ini that is missing if v_bc is included only from z_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.

Figure 3.

The bias factors |$b(k\, v_{\rm bc})$| at z_ini = 200 for the average v_bc in the zoom region, scaled back to its value at recombination v_{bc, rec} = 100.07 km s⁻¹, which is a |$\sim 3.3\sigma _{v_{\rm bc}}$| value. We show b(k, v_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 k > 40 h Mpc⁻¹, 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 b(k, v_bc) for the peculiar velocities for the full range of k, but the velocity field is dominated by large-scale (small-k) modes. Therefore, b(k, v_bc) will have little if any impact on the small scales (large-k) of the velocity field.

Open in new tabDownload slide

We deal with the peculiar velocity field for the baryons in a similar way, by first converting the velocity divergence to peculiar velocities as |$\boldsymbol {v}_{\rm b}(k) = -{\rm i}a\boldsymbol {k}\theta _{\rm b}(k)/k^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 |$\boldsymbol {v}_{\rm b}$| equally. In reality, there would be preferential directions for the bias factor, depending on the direction of |$\boldsymbol {v}_{\rm bc}$|⁠, but we defer that implementation to future work.

Fig. 4 shows a 1.53 h⁻¹ kpc thick slice through the highest resolution level of the zoom ICs directly from music (‘unmodified’, left-hand column) and after the bias factor b(k, v_bc) has been applied (‘modified’, right-hand column). For δ_b (top row), the unmodified ICs contain a lot of small-scale structure, which is almost totally washed out after applying b(k, v_bc). Most of what remains is in the form of lower amplitude, larger scale fluctuations. For v_{b, i} (bottom rows),⁴ there is less small-scale structure to begin with, since the peculiar velocity fields are dominated by large scales. The effect of b(k, v_bc) on v_{b, i} is therefore much less striking than on the δ_b, with the main effect being smoothing and a slight reduction in amplitude.

Figure 4.

Slices of the unmodified (left-hand column) and modified (through convolution with the bias factor, as in equation (7), right-hand column) baryon overdensity (top row), peculiar velocity in the x-direction (second row), y-direction (third row), and z-direction (bottom row) in the high-resolution zoom region, of side length 543 h⁻¹ kpc. Each pixel corresponds to a cell width of 1.53 h⁻¹ kpc and the slice has a thickness of one cell width. The effect of applying b(k, v_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.

Open in new tabDownload slide

3.4 Haloes

After the ICs have been correctly initialized with v_bc, we can characterize the effect of v_bc on structure formation, principally by exploring how haloes are affected. Haloes are identified using ahf (Gill, Knebe & Gibson 2004; Knollmann & Knebe 2009), which supports multiresolution data sets and calculates baryonic properties of haloes. In order to use ahf with a ramses data set, we use the supplied ramses2gadget tool to convert leaf cells of the AMR hierarchy into gas pseudo-particles placed at the centre of each cell. This conversion is a standard recipe in the analysis of AMR data sets (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.

4 RESULTS

4.1 Halo abundances

Throughout this section, we calculate the cumulative number of haloes N(> M) as the number of haloes with a total mass greater than M and estimate the Poisson uncertainty on N(> M) for each case as |$\sqrt{N(\gt M)}$|⁠. To begin with, we ignore the conditions described in Section 3.4, instead considering the cumulative number of all haloes formed in the simulation N_all(> M), shown in Fig. 5. Aside from a slight suppression in N_all(> M) for the v_bc–rec case below 10⁶ h⁻¹ M_⊙ at z = 14.2, the v_bc–ini and v_bc–rec cases are consistent both with each other and with the no v_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 v_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).

Figure 5.

Ratio of the cumulative number of haloes N_all(> M) in the v_bc–ini (pink short-dashed) and the v_bc–rec (blue dotted) cases to the no v_bc case at z = 14.2 (top), 12.6 (middle), and 11.2 (bottom). All the haloes found by ahf are included in N_all(> M). The shaded regions indicate the 1σ Poisson uncertainty on the ratio. N_all(> M) is broadly the same between all three sets of simulations, except for a suppression in the number of haloes with M < 10⁶ h⁻¹ M_⊙ at z = 14.2 in the v_bc–rec case compared to the no v_bc case.

Open in new tabDownload slide

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 N_clean(> M), which shows a much more striking suppression of N_clean(> M) for the v_bc–rec case and an apparent increase in the abundance of low-mass (M < 10⁶ h⁻¹ M_⊙) haloes for the v_bc–ini case at z = 11.2 (though still consistent with no difference to the no v_bc case at the 1σ level).

Figure 6.

As in Fig. 5, but this time N_clean(> M) includes only haloes 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 v_bc–ini curve, which jumps from a suppression to a boost in the number of haloes with M < 2 × 10⁶ h⁻¹ M_⊙ between z = 12.6 and 11.2.

Open in new tabDownload slide

To highlight the effect of the cleaning procedure we also show the raw cumulative number of haloes at z = 14.2 and z = 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 N_clean(> M)/N_all(> M). 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.

Figure 7.

Cumulative number of haloes (top panel) for the raw catalogue N_all(> M) (solid) and the catalogue cleaned according to the prescription in Section 3.4 N_clean(> M) (long-dashed) and the ratio of N_clean(> M) to N_all(> M) (bottom panel) at each z. We show N(> M) at z = 14.2 (thin) and 11.2 (thick). From the top panel, it can be seen that N_all(> M) is very similar between all runs at each z shown, except for a slight suppression in the abundance of haloes in the v_bc–rec case at z = 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.

Open in new tabDownload slide

4.2 Baryon fraction

We allow star formation in these runs, so the total baryon fraction of a given halo is defined as follows:

where M_g is the gas, M_⋆ the stellar, and M_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 〈f_b〉_M and associated 1σ errorbars as function of z. We show each z where ≥30 haloes have formed that satisfy the criteria in Section 3.4, starting from z = 13.6 where we are able to match 38 haloes between the three cases. The gas fraction is suppressed at all z for the v_bc–ini and v_bc–rec cases compared to the no v_bc case. At earlier z, the suppression is stronger, though even by the final snapshot at z = 11.2, 〈f_b〉_M for both the v_bc–rec and v_bc–ini cases are not within 1σ of the no v_bc case. Notably, at all z, 〈f_b〉_M in v_bc–ini and v_bc–rec cases are almost indistinguishable from, and certainly consistent with, one another.

Figure 8.

Mass-weighted average baryon fraction f_b as a function of redshift z, normalized to the cosmic mean Ω_b/Ω_m. The average is shown for all z where more than 30 haloes have formed. The errorbars show the 1σ standard deviation as calculated using equation (10). We find a suppression in both the |$\rm{v_{\rm bc}--ini}$| and |$\rm{v_{\rm bc}--rec}$| cases compared to the no v_bc run, though between the two runs with v_bc there is little difference.

Open in new tabDownload slide

4.3 Star formation

In Fig. 9, we show the cumulative M_⋆ formed in the simulation, not accounting for mass loss due to supernovae, and the corresponding number of stellar particles N_⋆, which each have a mass of 108.0 h⁻¹ M_⊙. In each case, all of the star particles in the simulation formed inside a single halo. In total, 29 star particles formed by z = 11.2 in the no v_bc case, 10 in the v_bc–ini case, and seven in the v_bc–rec case. This hierarchy persists across all z, with more star particles having formed in the no v_bc case than in the v_bc–ini case and fewer still in the v_bc–rec case.

Figure 9.

Cumulative stellar mass M_⋆ formed as function of redshift z. Also, shown is the corresponding number of stellar particles N_⋆. Despite the small number of stars formed in this simulation, there is a clear hierarchy in how many stars each forms, with the no v_bc run forming the most, followed by |$\rm{v_{\rm bc}--ini}$|⁠, and then |$\rm{v_{\rm bc}--rec}$|⁠. There is also a clear delay in the onset of star formation, which follows a similar trend, with the |$\rm{v_{\rm bc}--rec}$| run starting to form stars the latest.

Open in new tabDownload slide

In the following, all times quoted in Myr are measured relative to the Big Bang. To quantify how v_bc impacts star formation in our simulations, we can calculate the delay in formation time of the nth star t_n for each run with v_bc compared to the no v_bc run |$t_{n,\, \text{no}~v_{\rm bc}}$| as |$\Delta t_n = t_n-t_{n,\, \text{no}~v_{\rm bc}}$|⁠. For the v_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 v_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 〈Δt〉 = 19.4 Myr and 34.9 Myr for the v_bc–ini and v_bc–rec cases, respectively.

5 DISCUSSION

To understand the perceived uptick in v_bc–ini haloes in Fig. 6, we need to explore the raw cumulative number of haloes, shown in Fig. 7 at z = 14.2 and z = 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 v_bc–ini case, because over the mass range 3 × 10⁵ < M h⁻¹ M_⊙ < 10⁶ fewer haloes are removed in the v_bc–ini case than in the no v_bc case. Since the cumulative numbers of haloes are already very similar, removing a smaller number of haloes in the v_bc–ini case manifests itself as a boost relative to the no v_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 v_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 v_bc also significantly affects the baryon fraction f_b, where we see that the mass-weighted baryon fraction 〈f_b〉_M is suppressed at all redshifts in both cases, with the suppression stronger at higher redshift. Even by z = 11.2, the v_bc–ini and v_bc–rec cases are still not in agreement with the no v_bc case, though the difference between the two populations has decreased. Again, this is likely due to the decay in the magnitude of v_bc, which allows the haloes to accrete more gas. Interestingly, 〈f_b〉_M is almost indistinguishable between the v_bc–ini and v_bc–rec cases, suggesting that including v_bc from z_ini is sufficient to observe its impact on halo baryon fraction, though a larger sample is required to confirm this. The suppression in 〈f_b〉_M when v_bc is included is in qualitative agreement with previous studies.

This decrement in baryon fraction for the v_bc–ini and v_bc–rec cases is reflected in the cumulative stellar mass formed, as fewer star particles formed in both cases than in the no v_bc case. Not only do they form fewer star particles, they also start forming star particles later since the effect of v_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 v_bc–ini case and 34.9 Myr for the v_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 and Pop II (i.e. from metal-enriched 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 〈f_b〉_M between the v_bc–ini and v_bc–rec cases at most redshifts, there is a clear hierarchy in the amount of stars formed – no v_bc forms the most, v_bc–ini forms fewer, and v_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 v_bc–ini case than in the v_bc–rec case.

6 CONCLUSIONS

We have performed the first cosmological zoom simulations that self-consistently sample the relative baryon-dark matter velocity v_bc from a large 400 h⁻¹ Mpc box. This relative velocity naturally arises when simulations are initialized 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 initializing simulations in this manner misses out on the effect of the relative velocities from z = 1000 to the start time of the simulation, z_ini. We developed a methodology that compensates for the effect of v_bc on baryon density and velocity perturbations by computing a ‘bias’ factor b(k, v_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 h⁻¹ Mpc subbox, extracted from the main 400 h⁻¹ Mpc box. The zoom region is centred on the region with the largest relative velocity in the 400 h⁻¹ Mpc box, which has an RMS value of v_bc = 100.07 km s⁻¹ at z = 1000, corresponding to |$\sim 3.3\sigma _{v_{\rm 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 z = 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,⁵ 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 a Science and Technology Facilities Council studentship and acknowledges the indirect support of a Royal Society Dorothy Hodgkin Fellowship and a Royal Society Enhancement Award. ITI was supported by the Science and Technology Facilities Council (grant numbers ST/I000976/1 and ST/T000473/1) and the Southeast Physics Network. AF is supported by the Royal Society University Research Fellowship URF/R1/180523. This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1, ST/R002371/1, and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time through the John von Neumann Institute for Computing (NIC) on the GCS Supercomputer JUWELS at Jülich Supercomputing Centre (JSC). We acknowledge that the results of this research have been achieved using the PRACE Research Infrastructure resource Galileo based in Italy at CINECA. We acknowledge PRACE for awarding us access to Beskow/Dardel hosted by the PDC Center for High Performance Computing, KTH Royal Institute of Technology, Sweden. This work made use of the following software packages: numpy (Harris et al. 2020), matplotlib (Hunter 2007), scipy (Virtanen et al. 2020), yt (Turk et al. 2011); ytree (Smith & Lang 2019), hmf (Murray, Power & Robotham 2013) and cmasher (van der Velden 2020).

DATA AVAILABILITY

The code for computing and applying the bias factor b(k, v_bc) is available at https://github.com/lconaboy/drft. The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

References

APPENDIX A: COMPARISON TO PREVIOUS WORKS

We run a series of test simulations set up as in Naoz et al. (2012, 2013).⁶ The specific case shown here has a base resolution ℓ_min = 9 (512³ dark matter particles and, initially, cells), in a 471.1 h⁻¹ kpc periodic box. The simulations in Naoz et al. (2012, 2013) were performed using the SPH code gadget2 (Springel 2005), while we use the AMR code ramses. We allow the AMR grid to refine freely up to ℓ_max = 14, corresponding to a maximum comoving resolution of 28.9 h⁻¹ pc, comparable to the gravitational softening length of 45.8 h⁻¹ pc comoving used in Naoz et al. (2012). We use our fiducial cosmology (Section 3.1), whereas Naoz et al. (2012) used a cosmology consistent with Komatsu et al. (2009) with parameters: Ω_m = 0.28, Ω_Λ = 0.72, Ω_b = 0.046, h = 0.7,⁷ and a boosted σ₈ = 1.4. We adopt the boosted σ₈ = 1.4 as used in the original studies. We initialize the simulations at z_ini = 200⁸ both with and without a relative velocity of |$1.7\sigma _{v_{\rm bc}}=10~\text{km~s}^{-1}$| at z_ini. We also compute and apply the bias factor b(k, v_bc) to the baryonic component of the ICs, while the Naoz et al. (2012) simulations are initialized by computing transfer functions, which explicitly include the effect of the v_bc. Following Naoz et al. (2012), we set both of the velocity fields to that of the dark matter and apply the v_bc to the x-component of the baryon velocity as follows:

Including the v_bc in this manner is justified by the box size being much smaller than the coherence scale of the v_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 |$\bar{\rho }_{\rm 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 v_bc, we compare to simulations without v_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 A as follows:

First, we look at the effect on the cumulative halo mass function N(> M), as in Naoz et al. (2012). Fig. A1 shows the decrement in N(> M) for the case with v_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.

Figure A1.

The fractional difference Δ_N in halo mass function N(> M) for the v_{bc, ini} = 10 km s⁻¹ case, calculated with equation (A2). We show Δ_N at z = 25 (top), 19 (centre), and z = 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 M < 6 × 10⁴ h⁻¹ M_⊙, which could arise because of the friends-of-friends algorithm used to initially find haloes in Naoz et al. (2012).

Open in new tabDownload slide

However, the overall shape of our Δ_N is slightly different to Naoz et al. (2012); we match well below ∼3 × 10⁵ h⁻¹ 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 ICs. One further significant source of difference is the cosmologies used. Fig. A2 shows the difference expected at z = 15 by comparing the analytic Watson et al. (2013) N(> M) mass functions (magenta solid). From this, we would expect the Naoz et al. (2012) simulation to have ∼9 per cent more haloes with M > 3 × 10⁵ h⁻¹ M_⊙. This increase in the number of haloes increases with mass and for M > 1 × 10⁷ h⁻¹ M_⊙, we would expect ∼11 per cent more haloes in the Naoz et al. (2012) simulation. Indeed this is borne out by the simulations (cyan dashed), which show that the Naoz et al. (2012) simulations do produce more haloes at all masses. At higher masses, Δ_N 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.

Figure A2.

Comparison of the fractional difference Δ_N in halo mass function N(> M) between our simulations and the simulations from Naoz et al. (2012) (cyan short-dashed) and between the Watson et al. (2013) curve, which is fit to N-body simulations, for the cosmology used in our work and the one used in Naoz et al. (2012) (magenta solid). The fractional difference is calculated as Δ_N = (N₁ − N₀)/N₀, where N₀ are the data corresponding to our work and N₁ to Naoz et al. (2012), so Δ_N > 0 means there are more haloes in Naoz et al. (2012). Neither run includes v_bc. The shaded cyan region indicates the combined 1σ Poisson uncertainty on N(> M) from the simulations. For M ≲ 4 × 10⁶ h⁻¹ M_⊙ the difference in N(> M) 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 (M ≳ 4 × 10⁶ h⁻¹ M_⊙) discrepancy are discussed in the text.

Open in new tabDownload slide

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 as follows:

Fig. A3 shows the binned gas fractions (top panel) for our and for the Naoz et al. (2013) simulations, each normalized to the cosmic mean Ω_b/Ω_m for the appropriate cosmology, and the decrement (bottom panel) as defined in equation (A2). We take the mid-point 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. 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 v_bc, we can remove artefacts due to the choice of code.

Figure A3.

Binned baryon fraction f_b (top panels) and relative fractional difference |$\Delta _{f_{\rm b}}$| between each run with v_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 z = 25 (left-hand panel), 19 (centre panel), and 15 (right-hand panel). 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 f_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.

Open in new tabDownload slide