Baryon-CDM isocurvature galaxy bias with IllustrisTNG

TL;DR

This paper investigates how primordial baryon-CDM isocurvature perturbations influence galaxy bias and clustering, using IllustrisTNG simulations to quantify effects relevant for cosmological observations.

Contribution

The study introduces a novel simulation approach to measure baryon-CDM isocurvature galaxy bias and models its impact on galaxy clustering and cosmological parameter inference.

Findings

Galaxy bias parameter $b_{\delta}^{bc}$ is negative and decreases with mass and redshift.

Stellar-mass selected galaxies show weaker or opposite bias trends.

Impact on power spectrum and cosmological parameters is below 1% and 0.1%, respectively.

Abstract

We study the impact that baryon-CDM relative density perturbations $\delta_{bc}$ have on galaxy formation using cosmological simulations with the IllustrisTNG model. These isocurvature (non-adiabatic) perturbations can be induced primordially, if multiple fields are present during inflation, and are generated before baryon-photon decoupling when baryons did not comove with CDM. The presence of long-wavelength $\delta_{bc}$ perturbations in our simulations is mimicked by modifying the ratios of the cosmic densities of baryons $\Omega_b$ and CDM $\Omega_c$, at fixed total matter density $\Omega_m$. We measure the corresponding galaxy bias parameter $b_{\delta}^{bc}$ as the response of galaxy abundances to $\delta_{bc}$. When selecting by total host halo mass, $b_{\delta}^{bc}$ is negative and it decreases with mass and redshift. Stellar-mass selected simulated galaxies show a weaker or even the opposite trend because of the competing effects of $\delta_{bc}$ on the halo mass function and stellar-to-halo-mass relations. We show that simple modeling of the latter two effects describes $b_{\delta}^{bc}$ for stellar-mass-selected objects well. We find $b_{\delta}^{bc} =0.6$ for $M_* = 10^{11}\ M_{\odot}/h$ and $z=0.5$, which is representative of BOSS DR12 galaxies. For $\delta_{bc}$ modes generated by baryon-photon interactions, we estimate the impact on the DR12 power spectrum to be below $1\%$, and shifts on inferred distance and growth rate parameters should not exceed $0.1\%$.

Figures

Tables

Table 1. Table 1 : Parameters of the cosmologies simulated in this paper. The High and Low cosmologies describe the effect of Δ b = 0.05 subscript Δ 𝑏 0.05 \Delta_{b}=0.05 and Δ b = − 0.05 subscript Δ 𝑏 0.05 \Delta_{b}=-0.05 long-wavelength baryon-CDM density perturbations in the Fiducial cosmology (cf. Eqs. ( 2.2 ) and ( 2.3 )). Note that all of the cosmological parameters are the same, except Ω b subscript Ω 𝑏 \Omega_{b} and Ω c subscript Ω 𝑐 \Omega_{c} (this is the meaning of ”). We have simulated structure formation in these cosmologies at two particle resolutions: N p = 1250 3 subscript 𝑁 𝑝 superscript 1250 3 N_{p}=1250^{3} , L box = 75 ​ Mpc / h subscript 𝐿 box 75 Mpc ℎ L_{\rm box}=75\ {\rm Mpc}/h (called TNG100-1.5) and N p = 1250 3 subscript 𝑁 𝑝 superscript 1250 3 N_{p}=1250^{3} , L box = 205 ​ Mpc / h subscript 𝐿 box 205 Mpc ℎ L_{\rm box}=205\ {\rm Mpc}/h (called TNG300-2). Each simulation was also run without (dubbed Gravity) and with (dubbed Hydro) hydrodynamical physical processes (note that for the Hydro runs, the number of mass elements is twice the quoted values: N p subscript 𝑁 𝑝 N_{p} gas cells and N p subscript 𝑁 𝑝 N_{p} dark matter mass elements). The same random seed was used to generate the initial conditions of all the simulations. The value of the primordial power spectrum amplitude A s subscript 𝐴 𝑠 A_{s} is evaluated at a pivot scale k pivot = 0.05 ​ Mpc − 1 subscript 𝑘 pivot 0.05 superscript Mpc 1 k_{\rm pivot}=0.05\ {\rm Mpc}^{-1} ; this value yields σ 8 ​ ( z = 0 ) = 0.816 subscript 𝜎 8 𝑧 0 0.816 \sigma_{8}(z=0)=0.816 in the Fiducial cosmology.

Name	${\Omega }_m$	${\Omega }_b$	${\Omega }_c$	${\Omega }_{\Lambda }$	$h$	$n_s$	$A_s$
Fiducial	$0.3089$	$0.0486$	$0.2603$	$0.6911$	$0.6774$	$0.967$	$2.068\times {10}^{-9}$
High	”	$0.0510$	$0.2579$	”	”	”	”
Low	”	$0.0462$	$0.2627$	”	”	”	”

Equations

$${\delta }_g(\mathbf{x},z)=\sum _{\mathcal{O}}{b}_{\mathcal{O}}(z)\mathcal{O}(\mathbf{x},z),$$

$${\delta }_g(\mathbf{x},z)\supset b_{\delta }^{bc}(z){\delta }_{bc}(\mathbf{x})+b_{\theta }^{bc}(z){\theta }_{bc}(\mathbf{x},z),$$

$${\tilde{\Omega }}_b$$

$${\tilde{\Omega }}_c$$

$$1+{\delta }_{bc}$$

$$b_{\delta}^{bc}(z)=\frac{1}{n_{g}(z)}\frac{\partial n_{g}(z)}{\partial\delta_{bc}}\bigg{|}_{\delta_{bc}=0}.$$

$$b_{\delta }^{bc}(z,M)=\frac{b_{\delta }^{bc,\mathrm{High}}(z,M)+b_{\delta }^{bc,\mathrm{Low}}(z,M)}{2},$$

$$b_{\delta }^{bc,\mathrm{High}}(z,M)$$

$$b_{\delta }^{bc,\mathrm{Low}}(z,M)$$

$$n(M)$$

$$\frac{\mathrm{d}n(M)}{\mathrm{d}M}$$

$$f^{\mathrm{T}}(\sigma )$$

$${\sigma }^2=\frac{1}{2{\pi }^2}\int \mathrm{d}k{k}^2{P}_{{\delta }_m{\delta }_m}(k,z){\tilde{W}}^2(k,R(M)),$$

$$b_{\delta }^{bc,univ.}(z,M)=\frac{1}{{\delta }_{bc}}\left[\frac{n^{\mathrm{SepUni},univ.}(z,M)}{n^{\mathrm{Fiducial},univ.}(z,M)}-1\right],$$

$$b_{\delta }^{bc,univ.}(z,M_{\ast })=\frac{1}{{\delta }_{bc}}\left[\frac{n^{\mathrm{SepUni},univ.}\left(z,M_{\mathrm{h}}^{\mathrm{Fiducial}}\left(\left[1-ϵ{\Delta }_b\right]M_{\ast }\right)\right)}{n^{\mathrm{Fiducial},univ.}\left(z,M_{\mathrm{h}}^{\mathrm{Fiducial}}\left(M_{\ast }\right)\right)}-1\right].$$

$${\delta }_g(\mathbf{x},z)=b_1(z){\delta }_m(\mathbf{x},z)+b_{\delta }^{bc}(z){\delta }_{bc}(\mathbf{x}).$$

$$P_{gg}(k,z)$$

$$\frac{2b_1(z)b_{\delta }^{bc}(z)P_{{\delta }_m{\delta }_{bc}}(k,z)}{b_1(z{)}^2{P}_{{\delta }_m{\delta }_m}(k,z)}\approx \frac{2b_{\delta }^{bc}(z)}{b_1(z)}\frac{T_b(k,z)-T_c(k,z)}{T_m(k,z)},$$

$${\alpha }_{\perp }=\frac{D_A(z)r_d^{\ast }}{D_A^{\ast }(z)r_d};{\alpha }_{\parallel }=\frac{H^{\ast }(z)r_d^{\ast }}{H(z)r_d};f{\sigma }_8(z)=\frac{\mathrm{dln}D(z)}{\mathrm{dln}a}{\sigma }_8(z),$$

$${\dot{\delta }}_s(\mathbf{x},z)$$

$${\dot{\theta }}_s(\mathbf{x},z)+H(z){\theta }_s(\mathbf{x},z)$$

$${\ddot{\delta }}_r(\mathbf{x},z)+2H(z){\dot{\delta }}_r(\mathbf{x},z)$$

$${\dot{\theta }}_r(\mathbf{x},z)+H(z){\theta }_r(\mathbf{x},z)$$

$${\dot{\delta }}_r(\mathbf{x},z)$$

$${\delta }_r(\mathbf{x},z)={\delta }_{bc}(\mathbf{x})+R_{-}(\mathbf{x})D_r(z),$$

$$D_r(z)=H_0{\int }_{t(z)}^{\infty }\frac{d{t}^{\prime }}{a(t^{\prime }{)}^2}.$$

$${\delta }_r(\mathbf{x},z)$$

$${\theta }_{bc}(\mathbf{x},z)\equiv {\theta }_r(\mathbf{x},z)$$

$${\delta }_{bc}(\mathbf{k})={\delta }_r(\mathbf{k},z)-\frac{a{D}_r(z)}{H_0}{\theta }_{bc}(\mathbf{k},z).$$

$${\delta }_r(\mathbf{k},z)=\frac{T_b(k,z)-T_c(k,z)}{T_m(k,z)}{\delta }_m(\mathbf{k},z).$$

$${\theta }_{bc}(\mathbf{k},z)=i\mathbf{k}\cdot {\mathbf{v}}_{bc}=k\frac{T_{v_{bc}}(k,z)}{T_m(k,z)}{\delta }_m(\mathbf{k},z),$$

$$P_{{\delta }_m{\delta }_{bc}}(k,z)=\left[\frac{T_b(k,z)-T_c(k,z)}{T_m(k,z)}-\frac{k{D}_r(z)a}{H_0}\frac{T_{v_{bc}}(k,z)}{T_m(k,z)}\right]P_{{\delta }_m{\delta }_m}(k,z).$$

$$C^{\text{TNG300-2}}\left(M_{\mathrm{h}}\right)$$

Full text

1 Introduction

The study of the large-scale clustering of galaxies in the Universe is one of the most promising avenues to address long-standing open issues in fundamental physics such as the nature of dark matter and dark energy, what is the law of gravity on large scales and what is the value of neutrino masses. One of the main ingredients in the exercise of theoretically predicting the statistics of the galaxy distribution is a description of galaxy bias, i.e. a formalism that relates the observed distribution of galaxies to that of the underlying matter field (see Ref. [1] for a comprehensive review). Formally, the density contrast of the galaxies at position $\bm{x}$ and redshift $z$, $\delta_{g}(\bm{x},z)$, can be expanded as

[TABLE]

where the sum runs with all generality over all types of large-scale perturbations $\mathcal{O}(\bm{x},z)$ that can have an impact on galaxy formation. This expansion should also include stochastic terms, as well as projection and selection effects [2], but we refrain from writing these explicitly here. Physically, the bias parameters111The word “parameter” is a misnomer since they depend on time and also on galaxy properties (such as mass, luminosity, etc.). We nonetheless retain this nomenclature as it is usually adopted in the literature. encode all of the dependence of galaxy formation processes on the large-scale perturbations that each of them multiplies; more technically, the $b_{\mathcal{O}}(z)$ specify how the galaxy distribution “responds” to changes in the amplitude of the perturbations $\mathcal{O}(\bm{x},z)$.

A central question prior to any galaxy clustering study concerns therefore the number of terms that should be taken into account in the expansion of Eq. (1.1). For example, the so-called local-in-matter-density (LIMD) bias parameters correspond to including terms such as [3] $\delta_{g}(\bm{x},z)\supset\sum_{n}b_{n}(z)\delta^{n}_{m}(\bm{x},z)/n!$, where $\delta_{m}(\bm{x},z)\propto\nabla^{2}\Phi(\bm{x},z)$ is the total matter density contrast ($\Phi$ is the gravitational potential). This captures all of the dependency of galaxy formation on the amplitude of large-scale isotropic perturbations. The order $n$ up to which one should keep terms is determined by the order in perturbation theory [4] up to which one wishes to model the statistics of the galaxy distribution, which is in turn partly determined by how deep into the nonlinear regime of structure formation one wishes to analyse some given galaxy sample. For instance, to model the galaxy power spectrum at leading order (or tree-level) in perturbation theory, one would need to retain only the first-order term $\propto\delta_{m}$; next-to-leading order (or $1$-loop), one would need to include also the term $\propto\delta_{m}^{2}$, etc. In addition to the LIMD bias parameters, there are also important contributions from large-scale tidal fields [5, 6, 7, 8], e.g. $\delta_{g}(\bm{x},z)\supset b_{K^{2}}K^{2}_{ij}(\bm{x},z)$, where $K_{ij}(\bm{x},z)=(\partial_{i}\partial_{j}/\nabla^{2}-\delta_{ij}/3)\delta_{m}$; as well as from $\mathcal{O}(\bm{x},z)$ operators constructed from higher than second order derivatives of the potential such as $\nabla^{2}\delta_{m}(\bm{x},z)$ [9]. The bottom line is that, at a given order in perturbation theory, it is important to make sure that one enumerates all of the possible operators $\mathcal{O}(\bm{x},z)$ in Eq. (1.1) that can influence galaxy formation.

All the bias parameters mentioned so far have in common the fact that they are associated with perturbations that depend on the gravitational potential that is sourced by total matter density fluctuations (i.e., growing adiabatic perturbations); these are also the most widely studied bias parameters in the literature. In this paper, we turn our attention instead to relative perturbations between the baryon and cold dark matter (CDM) components, which have received less attention in the literature despite being potentially relevant. These isocurvature perturbations can be generated during inflation [10] in multifield scenarios [11, 12, 13, 14]. In addition, even for purely adiabatic primordial perturbations, they are sourced before baryon-photon decoupling, when baryons did not move along the same trajectories as dark matter because of their tight coupling to the photons. For both generation mechanisms, the consequences are long-wavelength modulations of the amount of baryons relative to CDM after decoupling [15, 16, 17, 18]. Given that galaxy formation and evolution is sensitive to the amount of baryons available to participate in processes such as baryonic accretion, star formation, black hole growth, and feedback, it is therefore important to study the impact of these baryon-CDM density perturbations. Likewise, after decoupling, there will also be regions in the Universe that exhibit relative baryon-CDM velocity perturbations, i.e. patches within which baryons and CDM move at different velocities [19, 20, 16]. A noteworthy aspect of baryon-CDM perturbations is that they exhibit strong baryon acoustic oscillation (BAO) features that are not completely in phase with those imprinted in the total matter fluctuations. This can lead to shifts in the BAO scale imprinted in the galaxy distribution, which should be accounted for to guard against potential biases on cosmological parameters.

More concretely, our main objective in this paper is to estimate the baryon-CDM density perturbation bias parameter $b_{\delta}^{bc}$. We shall do so by carrying out hydrodynamical cosmological simulations of galaxy formation with the AREPO code [21, 22] and the IllustrisTNG physics model [23, 24] in the presence of long-wavelength baryon-CDM density perturbations. These perturbations can be incorporated in the simulations by perturbing, relative to some fiducial cosmology, the cosmic fractions of baryons $\Omega_{b}$ and CDM $\Omega_{c}$, while keeping the total matter density $\Omega_{m}$ fixed. By way of the separate universe formalism, galaxy formation and evolution in this modified cosmology is equivalent to that taking place inside a long-wavelength baryon-CDM density perturbation in the fiducial cosmology. We will see that the effect of a baryon-CDM density perturbation impacts galaxy formation in two main ways: (i) the modified baryon-to-CDM ratio alters the shape of the initial matter power spectrum; and (ii) a modified baryon density alters the fuel supply for star formation, which results in modified stellar-to-halo-mass relations, as well as modified onset times for baryon feedback processes (such as active galactic nuclei (AGN) or supernovae feedback).

The amplitude of the baryon-CDM density bias parameter $b_{\delta}^{bc}$ has been previously estimated with simple analytical arguments in Refs. [15, 16]. The results that we present in this paper are, to the best of our knowledge, the first ever measurement of the baryon-CDM density bias from state-of-the-art galaxy formation simulations. Our numerical results will show that $b_{\delta}^{bc}$ can be sizeable for galaxy/halo mass scales and redshifts that are relevant to both current and future galaxy redshift surveys (while being well within the observational bounds reported by Ref. [25]; see also Refs. [26, 27]). This, together with the fact that $b_{\delta}^{bc}$ enters the bias expansion of Eq. (1.1) at leading order, makes it important to investigate the impact of including baryon-CDM density perturbations in models of galaxy clustering [17, 25, 18]. For adiabatic initial conditions after inflation, our results will show, however, that the impact of $\delta_{bc}$ on the galaxy power spectrum is not expected to exceed the $1\%$ level for mass scales and redshifts relevant for current and future galaxy surveys. These $\delta_{bc}$ perturbations are also not expected to bias distance, Hubble rate and growth rate measurements from galaxy samples like BOSS DR12 by more than $0.1\%$. Here, we do not study the impact of baryon-CDM velocity perturbations and we defer such an investigation to future work (see, however, Refs. [19, 28, 29, 30, 31, 32, 33, 20, 34, 16, 18] for such past studies, including Refs. [35, 36, 37, 38, 39, 40, 41, 42, 43, 44] for investigations with $N$-body simulations and Ref. [45] for a review).

This paper is organized as follows. In Sec. 2, we outline the main aspects of how baryon-CDM perturbations can impact large-scale structure formation and describe the numerical simulation setup we adopt to determine its associated density bias parameter $b_{\delta}^{bc}$. Section 3 contains our main numerical results: in Sec. 3.1, we present and discuss our results using simulated galaxies selected by their total host halo mass, whereas in Sec. 3.2, we do the same but selecting galaxies by their stellar mass. In Sec. 4, we estimate the impact that baryon-CDM density perturbations can have on the galaxy power spectrum. We summarize and conclude in Sec. 5. In App. A, we provide more details about the equations of baryon-CDM perturbations and in App. B, we describe a resolution correction strategy that we implement to self-consistently compare galaxy stellar mass values at different IllustrisTNG resolutions.

2 Baryon-CDM perturbations and Separate Universe simulations

In this section, we lay down the basics of the contribution of baryon-CDM perturbations to the general galaxy bias expansion and anticipate the main physical impact they will have on galaxy formation. We also describe the separate universe simulations that we perform with the IllustrisTNG model to obtain our numerical results.

2.1 Baryon-CDM perturbations in the galaxy bias expansion

Prior to the epoch of recombination (or more precisely, baryon-photon decoupling), the baryon and CDM components behave differently: the non-interacting CDM component is collisionless and experiences only gravity, but the baryons are tightly coupled to the photons and experience additional pressure forces that keep them from following the same trajectories as CDM. In other words, after baryons decouple from the photons, they are not comoving with the CDM component and cannot therefore be strictly treated as a single fluid [46, 47, 48, 49]. A practical consequence of this that is relevant for galaxy formation is that there will be regions in the Universe that exhibit relative velocity [19, 20, 16] and relative density perturbations [15, 16, 17] between the baryon and CDM components. These perturbations, which we call baryon-CDM perturbations here, are not normally taken into account in studies of large-scale structure formation, but it is important to move beyond (or at least assess the degree of validity of) this approximation since galaxy formation physics depends sensitively on the relative amounts of baryons and CDM. These perturbations are guaranteed to arise because of the baryon-photon interactions, but relative density perturbations can also be produced earlier if multiple fields are present during inflation; in the literature, these are often called compensated isocurvature perturbations (CIP) [50, 51, 10, 52, 53, 54, 27, 55].

In terms of the bias expansion of Eq. (1.1), these baryon-CDM perturbations contribute to first order as

[TABLE]

where $\delta_{bc}(\bm{x})$ is a constant compensated baryon-CDM perturbation characterized by $\delta_{c}=-f_{b}\delta_{b}$, $\delta_{m}=0$ with $\delta_{b}$, $\delta_{c}$ the density contrasts of baryons and CDM, respectively, and $f_{b}=\Omega_{b}/\Omega_{c}$ (we neglect neutrino masses, so that neutrinos do not contribute to matter). The term $\theta_{bc}(\bm{x},z)$, on the other hand, corresponds to a relative velocity divergence $\theta_{bc}={\bm{\nabla}}\cdot\bm{v}_{bc}$, with $\bm{v}_{bc}=\bm{v}_{b}-\bm{v}_{c}$ the relative velocity between baryons and CDM; the amplitude of this mode decays with time (see App. A for more details). There is already a significant literature [15, 20, 16, 17] discussing these terms, as well as some observational constraints. For example, Ref. [25] using the galaxy power spectrum from the BOSS DR12 sample finds (assuming photon-baryon interactions as the production mechanism) $b_{\delta}^{bc}=-1.0\pm 6.2$ and $b_{\theta}^{bc}=-114\pm 175$ at the $95\%$ confidence level. In Refs. [26, 27] the authors looked for the impact of baryon-CDM density perturbations by comparing number- and luminosity-weighted galaxy statistics, which are expected to be sensitive to baryon-CDM perturbations; the analysis of Ref. [27] is also consistent with a null detection.

In this paper, we focus on the bias parameter $b_{\delta}^{bc}$, i.e. we are interested in studying galaxy formation within constant baryon-CDM density perturbations. According to the estimate of Ref. [16], this is also the term that is expected to be the most important in the bias expansion (see also Ref. [18]). We will further be interested in cases in which the wavelength of these perturbations is assumed to be sufficiently large compared to the typical scale of galaxy formation.222If $R_{*}$ denotes the size of this scale, then the requirement is for the wavenumber of the baryon-CDM perturbations to be sufficiently smaller than $2\pi/R_{*}$. For halos, this scale is of order the size of their Lagrangian radius, $R_{*}=R_{\rm Lag.}^{\rm halo}\sim 10\ {\rm Mpc}/h$, which implies $k\ll 2\pi/R_{*}\approx 0.6\ h/{\rm Mpc}$. The same scale for galaxies is however more uncertain and not necessarily the same. For instance, Ref. [56] noted that radiative-transfer effects during reionization introduce a new scale in the bias expansion for galaxies, of order the mean free path of ionizing photons. “Sufficiently” here means that the modulation of the baryon-CDM perturbations effectively acts as a modified background to the physical processes that govern galaxy formation. In this limit, one can then make use of the separate universe ansatz to gain intuition about the expected phenomenology, as well as to setup the numerical simulations needed to study the effect (see the next subsection). The separate universe ansatz states that local structure formation inside a long-wavelength perturbation in some fiducial cosmology is equivalent to global structure formation in an appropriately modified cosmology. For the case of baryon-CDM density perturbations, the modifications to the cosmology involve altering the background densities of baryons and CDM, $\Omega_{b}$, $\Omega_{c}$, while keeping the total matter density $\Omega_{m}=\Omega_{b}+\Omega_{c}$ the same. Physically, one can then anticipate two main ways through which baryon-CDM density perturbations impact large-scale structure formation:

1. (1)

Changes in the relative sizes of $\Omega_{b}$, $\Omega_{c}$ lead to a difference in shape of the linear matter power spectrum on scales $k\gtrsim k_{\rm eq}\approx 0.02\ h/{\rm Mpc}$ (larger $\Omega_{b}$ suppresses small-scale power). This effect impacts structure formation, even if the latter takes place in a purely gravitational, pressureless manner. 2. (2)

Different baryonic densities $\Omega_{b}$ will also modify the fraction of the total amount of non-relativistic matter that can experience non-gravitational forces, undergo radiative cooling, form stars and black holes, which in turn feedback onto the rest of the matter via supernovae explosions and gas ejected by AGN.

Effect (1) above can be split further into two main physical effects. First, at fixed total matter $\Omega_{m}$, an increase in the baryon density leads to a corresponding reduction in the dark matter density, and correspondingly less growth of structure between the end of radiation domination and baryon-photon decoupling. This suppresses the amplitude of the power spectrum on scales smaller than the horizon at the epoch of matter-radiation equality, $k>k_{\rm eq}$ (note that $k_{\rm eq}$ remains the same since $\Omega_{m}$ is fixed). A second physical effect is associated with the change in the sound speed of the photon-baryon plasma, which depends on the ratio of baryon to photon densities, and which impacts the total matter power spectrum via a modified BAO feature. This second effect has been studied recently in Ref. [55] in the context of baryon-CDM density perturbations generated during inflation, i.e. primordial CIPs.

For $\delta_{bc}$ perturbations generated during inflation and that are still outside the sound horizon at photon-baryon decoupling (cf. App. A), the size of effect (1) can be calculated directly using Einstein-Boltzmann codes with modified baryon and CDM cosmic fractions; this is the approach we take in this paper. For perturbations that are inside the horizon at the time of decoupling, including those generated by baryon-photon interactions, the calculation of the initial power spectrum is more involved, as $\delta_{bc}$ then evolves with time in a scale-dependent way before converging to a constant $\delta_{bc}$ value after recombination is complete. This means that the effect of a large-scale $\delta_{bc}$ mode on the evolution of small-scale modes before decoupling, which underlies the effect (1), depends on the wavelength of the mode and cannot be captured precisely by varying $\Omega_{b}/\Omega_{c}$; the latter can only describe the regime where $\delta_{bc}$ is constant in time. We argue in App. A that this fact actually leads to an overestimate of the impact of effect (1) above for $\delta_{bc}$ modes generated solely due to baryon-photon interactions. Effect (2), on the other hand, does not depend on the exact past evolution of $\delta_{bc}$, but only on its value at the starting time of the simulation. We will return to these points whenever relevant to the interpretation of our results below.

2.2 Separate universe simulations of baryon-CDM perturbations

As already mentioned above, the effects of baryon-CDM density perturbations on galaxy formation can be mimicked by changes in the $\Omega_{b}$ and $\Omega_{c}$ cosmological parameters.333The implementation of the separate universe technique for the case of baryon-CDM density perturbations is even more straightforward than for the case of total density perturbations [57, 58, 59, 60, 61, 62, 63, 64, 65, 66]. In the latter, the changes in cosmology alter the relation between comoving and physical distances, as well as the relation between redshift and physical time, which implies performing additional (although straightforward) adjustments to the simulation box size, desired output redshift values and structure-finding criteria. These are steps that we do not have to worry about here as baryon-CDM density perturbations keep the value of $\Omega_{m}$ unchanged. Here, we describe these changes with a parameter $\Delta_{b}$ as

[TABLE]

where a tilde indicates a quantity in the modified cosmology. The value of $\Delta_{b}$ is related to the $\delta_{bc}$ term that enters the bias expansion Eq. (2.1) as

[TABLE]

which follows from $\tilde{\Omega}_{b}/\tilde{\Omega}_{c}=[1+\delta_{bc}]\Omega_{b}/\Omega_{c}$. The equations of motion show that the relative density perturbation is a constant mode (cf. App. A). Physically, this is because the ratio of baryon to CDM densities is conserved under gravitational evolution. It is thus consistently absorbed in modified cosmological parameters $\Omega_{b}$ and $\Omega_{c}$. Then, given the number density $n_{g}$ of galaxies selected according to some property (e.g. total mass or stellar mass), the bias coefficient $b_{\delta}^{bc}$ can be computed as

[TABLE]

Table 1 summarizes the three cosmological scenarios we consider in this paper. The two separate universe cosmologies, which we call High and Low, are obtained from the Fiducial cosmology by considering $\Delta_{b}=\Delta_{b}^{\rm High}=0.05$ and $\Delta_{b}=\Delta_{b}^{\rm Low}=-0.05$, respectively. These numerical values of $\Delta_{b}$ are chosen from a compromise between having sizeable effects in the simulations (i.e., sufficiently high signal-to-noise in the evaluation of Eq. (2.4); see also Eqs. (2.6) and (2.7) below), while keeping higher-order corrections small. All cosmologies share the same numerical values of the matter and cosmological constant densities $\Omega_{m},\Omega_{\Lambda}$, dimensionless present-day Hubble rate $h$, spectral index $n_{s}$ and amplitude of the primordial power spectrum $A_{s}$ (evaluated at a pivot scale $k_{\rm pivot}=0.05\ {\rm Mpc}^{-1}$); the numerical value of $A_{s}$ is that which yields $\sigma_{8}(z=0)=0.816$ in the Fiducial cosmology.

We obtain our numerical results using the moving-mesh code AREPO [21, 22] together with the IllustrisTNG galaxy formation model [23, 24], which is an improved version of its precursor Illustris [67, 68]. We refer the interested reader to Refs. [69, 70, 71, 72, 73] for the first results from the IllustrisTNG simulations, as well as Ref. [66] for separate universe simulations of total matter density perturbations; see also Ref. [74] for an overview of the publicly available simulation data.

We perform simulations at two mass resolutions: $N_{p}=1250^{3}$ particles on a cubic box with $L_{\rm box}=75\ {\rm Mpc}/h$ on a side, called TNG100-1.5 throughout, and $N_{p}=1250^{3}$, $L_{\rm box}=205\ {\rm Mpc/h}$, which we call TNG300-2. This resolution classification is similar to that adopted in the original IllustrisTNG runs; our TNG100-1.5 case falls in between the TNG100-1 ($N_{p}=1820^{3}$) and TNG100-2 ($N_{p}=910^{3}$) resolutions. For each resolution, we run gravity-only simulations (dubbed Gravity throughout), as well as hydrodynamical IllustrisTNG simulations (dubbed Hydro; for these, the number of ”particles” is doubled: $N_{p}$ gas cells and $N_{p}$ dark matter particles). The Gravity runs are sensitive to effect (1) mentioned in Sec. 2.1 in isolation, whereas the Hydro runs are sensitive to both effects (1) and (2). Comparing the two will thus indicate the relative importance of the two effects.

We generate the initial conditions at the same initial redshift $z_{i}=127$ for the Fiducial, High and Low cosmologies with the N-GenIC code [75] using the Zel’dovich approximation. We use the CAMB code [76] to compute the linear matter power spectra at $z=0$, which we scale back to $z_{i}$ assuming no cosmic radiation density and then give to N-GenIC as input. We also use the same random white-noise seed in N-GenIC to generate the initial conditions for the three cosmologies for each resolution. Note that we run CAMB with constant $\Omega_{b}$ and $\Omega_{c}$, which as discussed in the previous subsection, amounts to considering the impact on the initial power spectrum of baryon-CDM perturbations generated during inflation (CIPs) and that are still outside the sound horizon at decoupling.

There are two other points worth noting about our numerical methodology, which are important for the interpretation of our results below. One is that the gas and CDM mass elements at the starting redshift are initialized with the same density perturbations and velocities. That is, structure formation in our simulations begins with the baryons comoving with the CDM component (this is as in the original IllustrisTNG simulations; see, however, Refs. [77, 78, 79] for studies about the impact of initializing baryons and CDM with different transfer functions). The other point is that we do not modify any of the parameters of the IllustrisTNG physics model when we modify $\Omega_{b}$ and $\Omega_{c}$. Our separate universe simulations thus follow structure formation in the fiducial cosmology with the hydrodynamical processes as specified by the IllustrisTNG model, inside a region at cosmic mean total matter density where baryons and CDM are comoving, and whose abundances are not the cosmic mean ones. An interesting question that we leave to be addressed in future work concerns the dependence of the predicted galaxy bias values on the baryon physics implementation in hydrodynamical simulations of galaxy formation.

The galaxies in our simulations are identified as gravitationally bound groups of stars as determined by the SUBFIND algorithm [80] in data snapshots produced at redshifts $z=3,2,1,0.5,0$. Following Eq. (2.4), we measure the baryon-CDM density bias $b_{\delta}^{bc}$ as

[TABLE]

with

[TABLE]

where, recall, $\delta_{bc}=(1+f_{b})\Delta_{b}$ and $N(z,M)$ denotes the number of galaxies found in the corresponding cosmology at redshift $z$ in some mass bin centered at $M$. For each cosmology and resolution, we have only simulated a single realization of the initial conditions, which prevents us from quoting errors in a statistical-ensemble sense. Theoretically, the values of $b_{\delta}^{bc,\rm High}$ and $b_{\delta}^{bc,\rm Low}$ should be the same, and hence we shall use their difference as a rough guide for the error in our measurements. Note also that Eq. (2.5) corresponds to a central finite difference, and hence the numerical error is of order $\delta_{bc}^{2}$ ($\approx 0.35\%$ for $\Delta_{b}=0.05$).

3 Numerical results: baryon-CDM density galaxy bias

In this section, we show measurements of $b_{\delta}^{bc}$ for galaxies selected as a function of the mass of all particles that are gravitationally bound to the halo, $M_{\rm h}$, as well as a function of the total mass in stars enclosed within twice the stellar half-mass radius (the radius that encloses half of the mass in stars bound to the halo), $M_{*}$.

3.1 Dependence on total halo mass

Figure 1 shows the baryon-CDM density bias parameter $b_{\delta}^{bc}$ as a function of total halo mass $M_{\rm h}$ measured from the Gravity runs. The result is shown for the TNG100-1.5 and TNG300-2 resolutions, and for different redshifts, as labeled. Our results show that, overall, $b_{\delta}^{bc}$ is a decreasing function of halo mass (it becomes more negative with $M_{\rm h}$) and that the amplitude of the effect is stronger at higher redshift. More specifically, $b_{\delta}^{bc}$ is always negative at $z>2$ for the mass scales shown, being close to zero for $M_{\rm h}\approx 10^{10}\ M_{\odot}/h$ at both $z=2$ and $z=3$. For halo masses $M_{\rm h}\approx 10^{13}\ M_{\odot}/h$, we find $b_{\delta}^{bc}\approx-1.5\ (z=2)$ and $b_{\delta}^{bc}\approx-3\ (z=3)$. At lower redshifts $z<1$, the mass dependence becomes less pronounced and the amplitude gets overall closer to zero: $b_{\delta}^{bc}\approx-0.3$ and $b_{\delta}^{bc}\approx-0.5$ for $M_{\rm h}\approx 10^{13}M_{\odot}/h$ at redshifts $z=0.5$ and $z=1$, respectively.

The result shown in Fig. 1 comes from the gravity-only runs and it therefore captures only the effect due to the modified shape of the linear matter power spectrum (cf. effect (1) in Sec. 2.1). This makes it possible to predict it using semi-analytical universal halo abundance formulae such as the Tinker et al. halo mass function [81], in which the number density of halos in a given mass bin $\left[M_{\rm min},M_{\rm max}\right]$ is given by

[TABLE]

where $\bar{\rho}_{m0}=3\Omega_{m0}H_{0}^{2}/(8\pi G)$ is the total physical matter density today and $A=0.186$, $a=1.47$, $b=2.57$, $c=1.19$ are parameters fitted to $\Lambda{\rm CDM}$ gravity-only simulations at $z=0$ for spherical-overdensity halos with mass definition $M_{\rm 200}$ (see Tab. 2 of Ref. [81]); the superscript T refers to Tinker. In the equations above, $\sigma$ is the variance of the total density field defined as

[TABLE]

where $\tilde{W}(k,R(M))=3\left({\rm sin}(kR)-kR{\rm cos}(kR)\right)/\left(kR\right)^{3}$ and $R(M)=\left(3M/(4\pi\bar{\rho}_{m0})\right)^{1/3}$. The linear matter power spectrum $P_{\delta_{m}\delta_{m}}(k)$ is the only ingredient in the expressions above that depends on the relative abundance of $\Omega_{b}$ and $\Omega_{c}$. We can thus predict the baryon-CDM density bias parameter using the universal mass function as

[TABLE]

where $n^{\rm Fiducial,univ.}(z,M)$, $n^{\rm SepUni,univ.}(z,M)$ are the universal halo mass function predictions computed with the Tinker fitting function and the linear matter power spectrum of the fiducial and separate universe cosmologies, respectively. The prediction of Eq. (3.5) is depicted by the black solid line in Fig. 1, which shows a very good agreement with the simulation results at all of the redshifts and mass scales shown; the larger scatter of the simulation results at higher masses simply reflects the decreased statistical precision due to the smaller number of objects with those masses. The result is also in line with the physical expectation that a boost in $\Omega_{b}$ (at fixed $\Omega_{m}$), i.e., positive $\delta_{bc}$, lowers the amplitude of the linear matter power spectrum on scales $k\gtrsim 0.02\ h/{\rm Mpc}$, which suppresses the formation of the most massive objects. Recall that the linear power spectrum was obtained assuming constant $\delta_{bc}$ at all times, and thus this result corresponds to super-sound horizon $\delta_{bc}$ modes generated during inflation. For the case of $\delta_{bc}$ modes generated by baryon-photon interactions, $b_{\delta}^{bc}$ is expected to be smaller in absolute value (i.e., less negative) because the suppression in the amplitude of the linear power spectrum for $k\gtrsim 0.02\ h/{\rm Mpc}$ is not as pronounced (cf. App. A).

Figure 2 shows the baryon-CDM density bias $b_{\delta}^{bc}$ as a function of total mass $M_{\rm h}$, but now measured from the Hydro runs. The result shown is now due to both the changes in the shape of the linear matter power spectrum (effect (1) in Sec. 2.1) and the modified hydrodynamics and star formation processes that follow from the different amount of baryons (effect (2) in Sec. 2.1). In Fig. 2, the universal mass function result (solid black) is the same as in Fig. 1, where it is shown to agree very well with the Gravity results. Thus, comparing the Hydro results with the universal mass function prediction in Fig. 2 allows to visualize the impact of effect (2) on the baryon-CDM density bias parameter. Qualitatively, the values of $b_{\delta}^{bc}$ measured from the Hydro runs show a mass and redshift dependence that is very similar to that in the gravity-only results. Quantitatively, the impact of hydrodynamical processes in IllustrisTNG makes the baryon-CDM density bias slightly more negative. This is noticeable at $z=3$ ($z=2$) for $M_{\rm h}\lesssim 10^{11}\ M_{\odot}/h$ ($M_{\rm h}\lesssim 10^{12}\ M_{\odot}/h$), and effectively all mass scales shown at $z<1$. For instance, at $z=0.5$ for $M_{\rm h}=10^{13}\ M_{\odot}/h$, the value of $b_{\delta}^{bc}$ is reduced from $\approx-0.3$ in the Gravity run to $\approx-0.55$ in the Hydro run (the noise in the measurement also increases from the Gravity to the Hydro runs). The suppression in the size of $b_{\delta}^{bc}$ (i.e., more negative) caused by the hydrodynamical processes suggests that the increased amount of baryons ($\delta_{bc}>0$) effectively results in amplified feedback effects that suppress overall the number of objects that form at a given total mass $M_{\rm h}$. Finally, we note also that the TNG100-1.5 and TNG300-2 Hydro results are in good agreement, which indicates that our $b_{\delta}^{bc}$ measurements as a function of total halo mass are not strongly affected by limited numerical resolution.

3.2 Dependence on galaxy stellar mass

We now turn our attention to the stellar mass dependence of the baryon-CDM density bias parameter. This measurement has to be made with some care since, at fixed total halo mass, different resolutions return slightly different stellar mass values. This is explored in detail in Ref. [70] (see their Appendix A), in which the authors devise a resolution correction strategy to ensure a more trustworthy comparison of quantities that use stellar masses (like the stellar mass function) at different IllustrisTNG resolutions. Specifically, Ref. [70] proposes using the ratio of stellar masses found at TNG100-1 and TNG100-2 resolutions as a multiplicative correction factor of the stellar masses found at TNG300-1 resolution, all at fixed total halo mass (see Eq. (A1) in Ref. [70]).444Important for this to work is the fact that TNG300-1 ($L_{\rm box}=205\ {\rm Mpc}/h$, $N_{p}=2500^{3}$) and TNG100-2 ($L_{\rm box}=75\ {\rm Mpc}/h$, $N_{p}=910^{3}$) have approximately the same mass resolution. The higher resolution TNG100-1 simulation has $L_{\rm box}=75\ {\rm Mpc}/h$, $N_{p}=1820^{3}$.

We have applied a similar resolution correction scheme using the original TNG100-2 ($L_{\rm box}=75\ {\rm Mpc}/h$, $N_{p}=910^{3}$) and TNG100-3 ($L_{\rm box}=75\ {\rm Mpc}/h$, $N_{p}=455^{3}$) simulations of the Fiducial cosmology to scale the stellar masses of our TNG300-2 simulations to values representative of TNG100-2 resolutions. Our resolution correction scheme in described in App. B, which shows also that the bias measured from the corrected TNG300-2 catalogues is nearly the same as that measured from the uncorrected ones (cf. right panel of Fig. 8). The reason is that our results are sensitive to relative differences between cosmologies with different $\Omega_{b}$, $\Omega_{c}$ values, which are less affected by numerical resolution issues compared to absolute values. In this section, we thus opt to show the results measured from our TNG300-2 simulations without any correction.

Figure 3 shows the stellar mass dependence of the baryon-CDM density bias parameter, $b_{\delta}^{bc}(M_{*})$. The figure shows that the two resolutions are in good agreement, although to a slightly lesser extent compared to the same measurement made in terms of total halo mass (cf. Fig. 2). The difference between $b_{\delta}^{bc,{\rm High}}$ and $b_{\delta}^{bc,{\rm Low}}$ (shaded area) is also larger, but one can nonetheless discern the main differences relative to the measurement made in terms of $M_{\rm h}$. Namely, the shape of $b_{\delta}^{bc}(M_{*})$ is appreciably flatter, with the simulation results suggesting even a slight increase with $M_{*}$ at $z<1$, as opposed to the monotonic decrease with $M_{\rm h}$ displayed in Fig. 2. Further, at redshifts $z<1$, our results are consistent with $b_{\delta}^{bc}(M_{*})$ being positive for most of the mass scales shown, whereas $b_{\delta}^{bc}(M_{\rm h})$ is always negative in Fig. 2.

As we describe next, the shape of $b_{\delta}^{bc}(M_{*})$ shown in Fig. 3 can be explained with the aid of a simple model based on universal mass function formulae and the halo-to-stellar-mass relation of the Fiducial cosmology. Figure 4 shows the impact that changes in $\Delta_{b}$ have on galaxy stellar masses as a function of halo mass. Specifically, the figure shows the ratio of the median stellar mass in total halo mass bins in the High and Low cosmologies to the Fiducial one. As expected, at fixed halo mass, an increase (decrease) in the amount of baryons in the High (${\rm Low}$) cosmology translates into an increase (decrease) in stellar mass. Further, the magnitude of these changes remains approximately constant across the total halo mass values shown; the two resolutions TNG100-1.5 and TNG300-2 also agree on this result. More quantitatively, the horizontal bands bracket $M_{\rm h}$-independent changes of $50\%$ to $100\%$ of the values of $\Delta_{b}$ ($\pm\left[0.025,0.05\right]$), which corresponds roughly to the changes observed in the simulations. This reflects the physical expectation that star formation efficiency is not entirely independent of the cosmic baryon density, and consequently not all of the extra baryons are turned into stars.

The result of Fig. 4 indicates that if $M_{\rm h}^{\rm Fiducial}\left(M_{*}\right)$ is the median halo-to-stellar-mass relation in the Fiducial cosmology, then the same relation in the separate universe cosmologies can be approximated as $M_{\rm h}^{\rm Sep.Uni.}(M_{*})\approx M_{\rm h}^{\rm Fiducial}\left(\left[1-\epsilon\Delta_{b}\right]M_{*}\right)$ where $\epsilon$ is a number between $0.5$ and $1$. One can then plug these halo-to-stellar-mass relations into the universal mass function formulae and generalize Eq. (3.5) as

[TABLE]

This equation’s prediction is shown by the shaded grey band in Fig. 3, which brackets the result from using $\epsilon=0.5$ and $\epsilon=1$. We evaluate $M_{\rm h}^{\rm Fiducial}\left(M_{*}\right)$ using a polynomial fit to the median relations found in the TNG100-1.5 simulation of the Fiducial cosmology. Figure 3 shows that the simple model of Eq. (3.6) succeeds at explaining the overall amplitude of the $b_{\delta}^{bc}(z,M_{*})$ measured in the simulations, as well as its dependence on $z$ and $M_{*}$. The level of agreement is not perfect, but this is as expected for at least three reasons: (i) the universal mass function with the Tinker fitting formula is already not a perfect description of the halo mass function in hydrodynamical simulations (cf. Fig. 1 vs. Fig. 2); (ii) the relation $M_{\rm h}^{\rm Fiducial}\left(M_{*}\right)$ we use is fitted to the median total halo mass found in a given stellar mass bin, which fails to capture the shape of the distribution within the bin; (iii) we have assumed that $\epsilon$ is constant with mass, which is an approximation valid only to the degree shown in Fig. 4.

To gain more insight on how much of the shape of $b_{\delta}^{bc}(M_{*})$ is due to the changes in the halo mass function and the halo-to-stellar-mass relation, we compare in Fig. 5 the result of Eq. (3.6) (grey band) with that obtained by considering only the changes to the $M_{\rm h}(M_{*})$ relation, i.e., by replacing $n^{{\rm SepUni},{\rm univ.}}$ with $n^{{\rm Fiducial},{\rm univ.}}$ in the numerator of Eq. (3.6) (orange band). Noting that the halo abundances are predominantly controlled by the shape of the power spectrum (effect (1) in Sec. 2.1; cf. Figs. 1 vs. 2), the orange band thus captures effect (2) in Sec. 2.1 in isolation. The figure shows that the lower the redshift, the greater the importance of effect (2). This result could have been anticipated from Figs. 1 and 4, which show that the impact on halo abundances decreases with redshift, but the changes to the stellar-to-halo-mass relations remain effectively constant. Specifically, at $z=0,0.5$ and $1$, effect (1) has barely any impact for $M_{*}\lesssim 10^{12},10^{11}$ and $10^{10}\ M_{\odot}/h$, respectively. On the other hand, for $z>2$, both effects (1) and (2) contribute sizeably, and $b_{\delta}^{bc}$ is given by the result of their competition.

Figure 5 allows us to also comment on the expected result for $\delta_{bc}$ modes generated by photon-baryon interactions. In the previous section, we noted that the impact on halo abundances registered in our simulations of constant $\delta_{bc}$ modes can be regarded as upper bounds to the impact of $\delta_{bc}$ modes generated by baryon-photon interactions (because the impact on the initial power spectrum is not expected to be as large). Further, the stellar-to-halo-mass relation is sensitive mostly to the relative amount of baryons and CDM at the onset of nonlinear structure formation; in other words, if $\delta_{bc}$ is already a constant at the starting time of the simulation, then the result of Fig. 4 is independent of the exact past evolution of $\delta_{bc}$. Putting these arguments together, we expect the $b_{\delta}^{bc}$ values for $\delta_{bc}$ generated by photon-baryon interactions to be bracketed by the orange and grey bands. Further, noting that at low redshift the two practically overlap, we conclude that for these redshifts the stellar mass dependence of $b_{\delta}^{bc}$ is fairly independent of the exact origin of the $\delta_{bc}$ modes (inflation vs. baryon-photon interactions).

4 Impact on the galaxy power spectrum

We can use the values of the baryon-CDM density galaxy bias parameter $b_{\delta}^{bc}$ measured from the simulations to estimate the corresponding impact on the galaxy power spectrum. We work to linear order, do not include redshift-space distortions and consider the following galaxy bias expansion (we keep the mass dependence of the bias parameters implicit to ease the notation):

[TABLE]

We focus only on the baryon-CDM density contribution, i.e. we disregard the contribution from the relative baryon-CDM velocity divergence term $\theta_{bc}(\bm{x},z)$ (whose bias parameter we cannot measure with our separate universe simulations). After Fourier-transforming, the galaxy power spectrum $P_{gg}(k,z)$ can be written as

[TABLE]

where $P_{\delta_{m}\delta_{bc}}$ is the cross power spectrum of $\delta_{m}$ and $\delta_{bc}$ and $P_{\delta_{bc}\delta_{bc}}$ is the auto power spectrum of $\delta_{bc}$. In this section, we assume that the $\delta_{bc}$ modes are generated by photon-baryon interactions before decoupling; the calculation of the corresponding spectra is described in App. A. In order to evaluate $P_{gg}(k,z)$ using Eq. (4.2), we need to specify the values of $b_{1}$ and $b_{\delta}^{bc}$. We consider galaxy samples like that of BOSS DR12 with typical stellar masses of order $M_{*}\approx 10^{11}\ M_{\odot}/h$ [82] at median effective redshifts $z_{\rm eff}\approx 0.5$. The value of $b_{\delta}^{bc}=0.6$ can be read off from the $z=0.5$ panel of Fig. 2 (black star). For $b_{1}$ we use the halo bias semi-analytical formulae of Ref. [83] assuming a typical host halo mass for BOSS galaxies with $M_{h}\approx 10^{13}\ M_{\odot}/h$ (in accordance with halo abundance matching analyses [82]). This yields $b_{1}=1.5$, which is also in line with the constraints on $b_{1}\sigma_{8}\approx 1.3$ obtained by Ref. [84] for the same sample.

The left panel of Fig. 6 shows the contribution from the $P_{\delta_{m}\delta_{m}}(k,z)$ and $P_{\delta_{m}\delta_{bc}}(k,z)$ terms in Eq. (4.2); the right panel shows the relative contribution of the $P_{\delta_{m}\delta_{bc}}(k,z)$ term to the total $P_{gg}(k,z)$. Overall, the figure makes apparent that the $b_{1}^{2}P_{\delta_{m}\delta_{m}}$ term is the dominant contribution for BOSS DR12-like galaxies, with the $P_{\delta_{m}\delta_{bc}}(k,z)$ accounting for a subpercent fraction that oscillates between $0.3\%$ and $0.7\%$ for the scales shown. The contribution from the term $\propto P_{\delta_{bc}\delta_{bc}}(k,z)$ in Eq. (4.2) is not shown, but we have explicitly confirmed that it is smaller in absolute value than the $P_{\delta_{m}\delta_{bc}}(k,z)$ term by over two orders of magnitude.

The relatively small impact of the baryon-CDM density perturbations on the galaxy power spectrum can be understood by inspecting the ratio of the $P_{\delta_{m}\delta_{bc}}$ to the $P_{\delta_{m}\delta_{m}}$ contributions in Eq. (4.2), which, in the case of adiabatic primordial perturbations considered in this section, is given by

[TABLE]

where $T_{b}$, $T_{c}$ and $T_{m}$ are the transfer functions of baryons, CDM and total matter, respectively (defined as $\delta_{s}(k,z)=(2k^{2}/(5\Omega_{m}H_{0}^{2}))\mathcal{R}(k)T_{s}(k,z)$, $s\in\{b,c,m\}$, where $\mathcal{R}(k)$ is the adiabatic primordial scalar perturbation generated by inflation); the approximation in Eq. (4.3) is valid at late times, as explained in App. A. The main suppression is due to the small absolute value of $(T_{b}-T_{c})/T_{m}$, which Ref. [16] had used to predict already $\approx 1\%$-level contributions using $b_{\delta}^{bc}=1$ (see also Ref. [15] for a similar earlier prediction and Ref. [18] for the impact of larger bias values). An additional small suppression comes from the ratio $2b_{\delta}^{bc}/b_{1}\approx 0.8$. Also, note that the bias parameters are positive (for stellar mass selection at $z=0.5$), but $T_{b}-T_{c}$ is negative, and hence, the baryon-CDM density perturbations contribute negatively to the galaxy power spectrum. Further, baryon-CDM perturbations also contribute to the cross spectrum of galaxies with total matter (which is relevant for galaxy-galaxy lensing measurements), but their relative contribution is suppressed by a similar amount following the same arguments.

The oscillatory behavior of the relative contribution of $P_{\delta_{m}\delta_{bc}}(k,z)$ is indicative of an offset of the phases of the BAO oscillations. In configuration space, i.e. at the level of the galaxy two-point correlation function, this can potentially lead to shifts in the position of the BAO peak, which is one of the main geometrical probes of the expansion history of the Universe. References [15, 16] estimated that the BAO peak position is likely shifted by less than $1\%$ by the baryon-CDM density term. Our results sharpen this conclusion, since we now have a likely range for the bias parameter $b_{\delta}^{bc}$. Figure 5 of Ref. [25] shows the expected biases on the distance and Hubble-rate parameters $\alpha_{\perp}$ and $\alpha_{\parallel}$, and growth rate of structure $f\sigma_{8}$ that would arise from ignoring $\delta_{bc}$ in observational analyses of the BOSS DR12 galaxy sample.555These parameters are defined as

$\displaystyle\alpha_{\perp}=\frac{D_{A}(z)r^{*}_{d}}{D^{*}_{A}(z)r_{d}}\ \ \ ;\ \ \ \alpha_{\parallel}=\frac{H^{*}(z)r^{*}_{d}}{H(z)r_{d}}\ \ \ ;\ \ \ f\sigma_{8}(z)=\frac{{\rm dln}D(z)}{{\rm dln}a}\sigma_{8}(z),$

(4.4) where $D_{A}^{*}(z)$ is the angular diamater distance to redshift $z$, $r^{*}_{d}$ is the sound horizon at photon-baryon decoupling and $H^{*}(z)$ is the Hubble rate, all evaluated in a reference cosmology assumed in the analysis; the same quantities without the superscript ∗ denote the true/unknown values that are to be inferred from the data. Further, $D(z)$ is the linear growth factor of total matter perturbations. At $95\%$ confidence level, the analysis of Ref. [25] constrains $b_{\delta}^{bc}=-1.0\pm 6.2$, which keeps the biases on $\alpha_{\perp}$, $\alpha_{\parallel}$ and $f\sigma_{8}$ below the $0.5\%$, $0.5\%$ and $2\%$ levels, respectively (cf. left middle panel of Fig. 5 of Ref. [25]). Our simulation results show that $b_{\delta}^{bc}\approx 0.6$ for galaxy stellar masses and redshifts representative of BOSS DR12. Taking the estimates of Ref. [25] as a guide, the expected biases caused by our $\delta_{bc}$ values are reduced to $\approx 0.1\%$ for all $\alpha_{\perp}$, $\alpha_{\parallel}$ and $f\sigma_{8}$. Furthermore, according to the forecasts of Ref. [18], the bias values measured from our simulations ($-1\lesssim b_{\delta}^{bc}\lesssim 1$; cf. Fig. 5) should remain below $1\sigma$ detection thresholds in surveys like DESI [85].

Before summarizing, we note that there is an interesting application of our results in the context of constraints on primordial compensated isocurvature perturbations (CIP) generated during inflation (see e.g. Refs. [50, 51, 10, 52, 53, 54, 27, 55, 86] for a number of existing observational constraints and forecasts). These contribute to the galaxy power spectrum as $\propto b_{1}b_{\delta}^{bc}A_{\rm CIP}P_{\delta_{m}\mathcal{R}}(k,z)$, where $A_{\rm CIP}=\Delta_{b}/\mathcal{R}$ is a parameter that determines the size of isocurvature perturbations $\Delta_{b}$ relative to the adiabatic curvature perturbations $\mathcal{R}$. This contribution has the same scale dependence as the one coming from primordial non-Gaussianity of the local type, whose amplitude is parametrized by the $f_{\rm NL}$ parameter [87]. Noting that our measurements give $b_{\delta}^{bc}$ values that are $\mathcal{O}(0.1-1)$ (the sign and exact value can depend on redshift and mass selection), it follows that galaxy surveys like SPHEREx [88] that aim to place $\mathcal{O}(1)$ constraints on $f_{\rm NL}$ should be able to place similarly tight constraints on $A_{\rm CIP}$, since both $f_{\rm NL}$ and $A_{\rm CIP}$ impact the galaxy power spectrum in the same way. More recently, Ref [86] presented also forecast constraints on CIP from contributions $\propto b_{\delta}^{bc}A_{\rm CIP}$ in the tomographic kinetic Sunyaev-Zel’dovich signal.

5 Summary and conclusions

In this paper, we have used hydrodynamical cosmological simulations to study the impact of baryon-CDM perturbations on galaxy formation. These are compensated isocurvature perturbations between baryons and CDM that can be generated during inflation if multiple fields are present. However, even for adiabatic initial conditions after inflation, these perturbations will also be produced during the epoch prior to baryon-photon decoupling, when baryons and CDM did not comove because of the tight coupling of the baryons to the photons (and lack thereof for CDM). Once the baryons decoupled from the photons, they were able to collapse gravitationally, but there were regions within which the amount of baryons relative to CDM differed from that at cosmic mean and baryons were moving with a different velocity relative to the CDM component. These baryon-CDM density and velocity perturbations can naturally have an impact on galaxy formation, but are not customarily included in hydrodynamical simulations, as well as in theoretical models of galaxy clustering statistics. Importantly, strong BAO features are imprinted on the statistics of the baryon-CDM perturbations, so assessing their contribution to galaxy clustering is important to ensure robust cosmology inference from the BAO feature.

We have focused on long-wavelength compensated baryon-CDM density perturbations $\delta_{bc}$ characterized by $\delta_{c}=-f_{b}\delta_{b}$ and $\delta_{m}=0$. Under the separate universe ansatz, galaxy formation taking place sufficiently inside these baryon-CDM density perturbations is effectively equivalent to structure formation taking place in a cosmology with modified baryon ($\Omega_{b}$) and CDM ($\Omega_{c}$) density parameters and fixed total matter ($\Omega_{m}$; cf. Tab. 1). There are two ways in which this change in cosmology impacts galaxy formation. The first is that an increased (decreased) value of $\Omega_{b}/\Omega_{c}$ lowers (enhances) the linear matter power spectrum on scales $k\gtrsim 0.02\ h/{\rm Mpc}$. In this paper, we have assumed that the value of $\Omega_{b}/\Omega_{c}$ is constant at all times after inflation, which maximizes this suppression on the matter power spectrum. For baryon-CDM perturbations generated by photon-baryon interactions, the long-wavelength perturbations in $\Omega_{b}/\Omega_{c}$ are time-dependent before decoupling and are expected to lead to a weaker modulation of the small-scale modes. The second effect is that modified relative amounts of baryons and CDM also result in different fractions of the total matter that can experience hydrodynamical forces and form stars.

Baryon-CDM density perturbations enter the galaxy bias expansion at leading order (cf. Eq. (2.1)) and our main objective was to measure the corresponding bias parameter $b_{\delta}^{bc}$. We have done so by measuring the response of galaxy number counts to changes in the amplitude of $\delta_{bc}$ (cf. Eqs. (2.4), (2.5), (2.6) and (2.7)). We have performed our simulations using the AREPO code with the IllustrisTNG galaxy formation model at TNG100-1.5 ($L_{\rm box}=75\ {\rm Mpc}/h$, $N_{p}=1250^{3}$) and TNG300-2 ($L_{\rm box}=205\ {\rm Mpc}/h$, $N_{p}=1250^{3}$) resolutions. For each resolution, we have run simulations without (referred to as “Gravity”) and with (“Hydro”) hydrodynamical processes taken into account. We have studied the redshift and mass dependence of $b_{\delta}^{bc}$ for two mass definitions: (i) the total mass due to all mass elements that belong to the host halo of the galaxy, $M_{\rm h}$; (ii) the mass in stars within twice the stellar half-mass radius, $M_{*}$. A summary of our main results is as follows:

• •

The bias parameter $b_{\delta}^{bc}$ becomes more negative with increasing host halo mass $M_{\rm h}$, and it does so more pronouncedly at higher redshift (cf. Fig. 2). Quantitatively, we find for $M_{\rm h}=10^{13}\ M_{\odot}/h$ that $b_{\delta}^{bc}\approx-0.5,-0.7,-1.5$ at $z=0.5,1,2$, respectively. At $z=1$, we have $b_{\delta}^{bc}\approx-0.1,-0.3,-0.7$ for $M_{\rm h}=10^{11},10^{12},10^{13}\ M_{\odot}/h$, respectively.

• •

The mass dependence of $b_{\delta}^{bc}$ becomes markedly different when galaxies are selected according to stellar mass $M_{*}$, instead of $M_{\rm h}$ (cf. Fig. 3). In particular, our results suggest a slight increase of $b_{\delta}^{bc}$ with $M_{*}$ and the redshift dependence is also less pronounced. At redshifts $z<1$, $b_{\delta}^{bc}$ is positive for $M_{*}\gtrsim 2\times 10^{10}\ M_{\odot}/h$; for $M_{*}\approx 10^{11}\ M_{\odot}/h$ at $z=0.5$, which is representative of BOSS galaxies, we find $b_{\delta}^{bc}=0.6$.

• •

The values of $b_{\delta}^{bc}$ for both halo-mass- and stellar-mass-selected tracers can be explained with a simple empirical model that takes into account (1) the effect of the modified linear power spectrum on the halo mass function via a universal mass function (Figs. 1 and 2); and (2) the modification of the stellar mass-halo mass relation $M_{*}(M_{\rm h})$ due to the changed baryon fraction (Fig. 4). For halo-mass-selected tracers, only the first effect is relevant in this simple model, while the bias $b_{\delta}^{bc}$ of stellar-mass-selected tracers is determined by the combination of both effects, with effect (2) being dominant at low redshifts (cf. Fig. 5).

• •

The contribution $2b_{1}b_{\delta}^{bc}P_{\delta_{m}\delta_{bc}}(k)$ from baryon-CDM density perturbations generated by baryon-photon interactions to the galaxy power spectrum is expected to lie below the $1\%$ level for galaxy samples like BOSS DR12 (cf. Fig. 6). Even taking into account the stronger BAO feature in $\delta_{bc}$, the impact of $b_{\delta}^{bc}$ is expected to be of order $0.1\%$ on the inferred distance and Hubble-rate parameters $\alpha_{\perp}$, $\alpha_{\parallel}$, and growth rate of structure $f\sigma_{8}$ (cf. Sec. 4), at least for an analysis pipeline similar to that of Ref. [25].

The impact of baryon-CDM perturbations is expected to be more pronounced at earlier times because the size of the constant $\delta_{bc}$ mode relative to the growing total matter density $\delta_{m}$ increases toward higher redshifts. This provides motivation to extend our baryon-CDM density perturbation study to higher redshifts (we limited ourselves to $z<3$ here), which would allow to study the impact of baryon-CDM density perturbations on the formation of the first halos, stars and galaxies. Further, at higher redshift, it would be interesting to go beyond the study of the impact of baryon-CDM perturbations on the total number of galaxies and measure the response of observables associated with the gas distribution. An example is the study of the response of the distribution of neutral Hydrogen, which could provide insights on the epoch of reionization ($z\approx 10-7$) and the signals of $21$-cm emission and Lyman-$\alpha$ absorption spectra.

It would also be interesting to extend our implementation of baryon-CDM perturbations in IllustrisTNG to include not only $\delta_{bc}$, but also the relative velocity divergence $\theta_{bc}$ and uniform relative velocity $\bm{v}_{bc}$ terms. Past simulation work focusing on $z\gtrsim 10$ [35, 36, 37, 38, 39, 40, 41, 42, 45, 43, 44] has shown that the presence of a uniform relative velocity suppresses the formation of low-mass haloes, as well as the rate at which gas can cool and form stars at early times. These simulations require however appreciably higher resolution than the ones currently performed with the IllustrisTNG model. The imprints that these early-time effects can leave on later-time structure formation (e.g., via modified reionization histories) are therefore currently uncertain, which motivates further investigations; we defer these to future work.

Finally, we note that our baryon-CDM bias measurements are specific to the IllustrisTNG galaxy formation model, and that it is plausible to expect the result to vary for varying implementations of baryon physics in cosmological simulations of galaxy formation. The bias parameters describe formally the environmental dependence of galaxy formation, and hence, a comparison of galaxy bias predictions from different baryon physics models can be used to take interesting steps in studies of galaxy formation.

Acknowledgments

We thank Raul Angulo, Simeon Bird, Jonathan Blazek, Emanuele Castorina, Bruno Henriques, Titouan Lazeyras, Kaloian Lozanov, Rüdiger Pakmor, Annalisa Pillepich, Volker Springel and Simon White for useful comments and discussions. The simulations used in this work were run on the Cobra supercomputer at the Max Planck Computing and Data Facility (MPCDF) in Garching near Munich. AB, GC and FS acknowledge support from the Starting Grant (ERC-2015-STG 678652) “GrInflaGal” of the European Research Council.

Appendix A Evolution of baryon-CDM perturbations

In this appendix, we present the equations that describe the evolution of baryon-CDM perturbations. We evaluate first the general solution of the baryon-CDM density and velocity divergence perturbations, and then we describe the calculation of the power spectra for those generated by baryon-photon interactions.

A.1 The constant and decaying modes

After photon-baryon decoupling, the linear density contrast $\delta$ and velocity divergence $\theta={\bm{\nabla}}\cdot\bm{v}$ of baryons and CDM obey the following equations

[TABLE]

where the subscripts b, c and m denote baryons, CDM and total matter respectively, $H(z)$ is the Hubble rate, $\bar{\rho}_{m}(z)$ is the background density of total matter and an overdot denotes a derivative w.r.t. physical time $t$; we use $t$, the scale factor $a$ and the redshift $z$ interchangeably as time variables. With these equations, we can write the following three equations (not independent) for the relative density, $\delta_{r}=\delta_{b}-\delta_{c}$, and relative velocity divergence $\theta_{r}=\theta_{b}-\theta_{c}$ between baryons and CDM:

[TABLE]

From Eq. (A.3), we know that the solution for $\delta_{r}$ is of the form (we follow closely the notation of Ref. [16])

[TABLE]

where $\delta_{bc}(\bm{x})$ and $R_{-}(\bm{x})$ are time-independent and $D_{r}(z)$ is a decaying function of time that satisfies $\ddot{D}_{r}+2H(z)\dot{D}_{r}=0$ and admits the following integral solution (the $H_{0}$ factor merely makes $D_{r}(z)$ dimensionless)

[TABLE]

Strictly speaking, $\delta_{bc}(\bm{x})$ and $R_{-}(\bm{x})$ should be evaluated at the Lagrangian position associated with $\bm{x}$, $\bm{q}[{\bm{x}},z]$; here, we always work with linear theory for which these complications can be ignored (actually, it should be $\bm{q}[{\bm{x}},z_{\rm dec}]$, where $z_{\rm dec}$ is the redshift at decoupling, but the difference is small and of the same order as other nonlinear terms at recombination that would not be considered anyway). Further, from Eq. (A.4), one knows that the solution is $\theta_{r}(\bm{x},z)={\theta_{bc,0}(\bm{x})}/{a}$, where $\theta_{bc,0}$ is the present-day value ($a=1$, $z=0$) of $\theta_{r}(\bm{x},z)$. Finally, plugging Eq. (A.6) and this $\theta_{r}(\bm{x},z)$ solution into Eq. (A.5) yields $R_{-}=\theta_{bc,0}/H_{0}$ (note that $\dot{D}_{r}=-H_{0}/a^{2}$ and $R_{-}\propto D_{r}/\dot{D}_{r}$, so the $H_{0}$ normalization in Eq. (A.7) cancels).

Putting it all together, the relative density and velocity divergence between baryons and CDM are given by

[TABLE]

In Eq. (A.9), we have introduced the notation $\theta_{r}\equiv\theta_{bc}$ that is used in the main text (and across most literature). These two equations make apparent that there are two baryon-CDM perturbation modes that can impact galaxy formation: the constant compensated baryon-CDM density perturbation that we studied in the main text, $\delta_{bc}(\bm{x})$, and the baryon-CDM velocity divergence $\theta_{bc}(\bm{x},z)$. The latter decays as $a^{-1}$ and its contribution to the $\delta_{r}$ perturbation becomes also smaller with time because the function $D_{r}(z)$ is decaying. Note that the results derived here are valid in general long after baryon-photon decoupling, regardless of what physical effect (primordial or pre-decoupling) set the initial conditions for $\delta_{bc}$ and $\theta_{bc}$.

A.2 The power spectra of baryon-CDM perturbations

We now focus on baryon-CDM perturbations generated by baryon-photon interactions before decoupling. Using Eqs. (A.8) and (A.9), we can write down the Fourier transform of $\delta_{bc}$ as (we distinguish Fourier- and real-space quantities by their arguments)

[TABLE]

Using that $\delta_{s}(k,z)=(2k^{2}/(5\Omega_{m}H_{0}^{2}))\mathcal{R}(k)T_{s}(k,z),s\in\{b,c,m\}$, we have

[TABLE]

Further, the velocity divergence term can be worked out as

[TABLE]

where the second equality defines the transfer function $T_{v_{bc}}(k,z)$ of the velocity difference between baryons and CDM $\bm{v}_{bc}=\bm{v}_{b}-\bm{v}_{c}$. We evaluate all of the transfer functions with the CAMB code [76]. The cross power spectrum between $\delta_{m}$ and $\delta_{bc}$, $\langle\delta_{m}(\bm{k},z)\delta_{bc}({\bm{k}^{\prime}})\rangle=(2\pi)^{3}P_{\delta_{m}\delta_{bc}}(k,z)\delta_{D}(\bm{k}+\bm{k}^{\prime})$ is then given by

[TABLE]

At the low redshifts we have considered in this paper, the term $\propto T_{v_{bc}}(k,z)$ is only a small contribution ($\approx 4\%$ at $z=3$ and $\approx 2\%$ at $z=1$, for $k=0.1\ h/{\rm Mpc}$). Hence, it is a good approximation to discuss the importance of $P_{\delta_{m}\delta_{bc}}(\bm{k},z)$ using only the first term (as we did in the discussion below Eq. (4.3)); numerically, however, we evaluate $P_{\delta_{m}\delta_{bc}}$ using all terms in Eq. (A.13). As a side remark, we note that according to Eqs. (21) and (26) of Ref. [25], their $b_{\delta}^{bc}$ parameter multiplies only the $(T_{b}-T_{c})/T_{m}$ contribution in Eq. (A.13). Strictly speaking, the constraints quoted in Ref. [25] do not correspond to the exact same bias parameter definition, but this does not have any practical consequence given the unimportance of the $T_{v_{bc}}$ term at late times.

The auto power spectrum of $\delta_{bc}$, $P_{\delta_{bc}\delta_{bc}}(k,z)$, can be calculated analogously by evaluating $\langle\delta_{bc}(\bm{k})\delta_{bc}({\bm{k}^{\prime}})\rangle$. Noting that $\delta_{m}$ and $\delta_{bc}$ are fully correlated, $P_{\delta_{bc}\delta_{bc}}$ can also be obtained from $P_{\delta_{m}\delta_{bc}}=\sqrt{P_{\delta_{m}\delta_{m}}P_{\delta_{bc}\delta_{bc}}}$, which shows that $P_{\delta_{bc}\delta_{bc}}$ is an even smaller contribution to the galaxy power spectrum than $P_{\delta_{m}\delta_{bc}}$.

A.3 The generation of $\delta_{bc}$ modes

Equations (A.1) and (A.2) admit a constant mode solution $\delta_{bc}(\bm{x})$ because they ignore the pressure forces that the baryons feel due to their coupling to the photons, i.e., they are valid only sufficiently after the epoch of photon-baryon decoupling. Prior to this epoch, these forces are sizeable and are what is in fact responsible for generating non-zero $\delta_{bc}$ for adiabatic perturbations after inflation. Figure 7 shows the time evolution of $\delta_{bc}$ as defined in Eq. (A.10), but with the coupling between baryons and photons appropriately taken into account using the CAMB code. The figure shows that, as expected, the modes that enter the horizon earlier (higher $k$), begin evolving earlier than larger-scale (lower $k$) modes. For example, the $k=0.1\ h/{\rm Mpc}$ mode has a wavenumber that is larger than the inverse sound horizon at decoupling, so it undergoes more than one full oscillation. Further, sufficiently large-scale modes, specifically modes that are still super-sound horizon by the time of decoupling, display negligible evolution in comparison (cf. blue line, $k=0.001\ h/{\rm Mpc}$). Importantly, however, after decoupling, the pressure that drives the generation of $\delta_{bc}$ becomes negligible and all $\delta_{bc}$ modes approach constant values.

For the case of baryon-CDM density perturbations generated during inflation, the picture remains qualitatively the same, except that the initial conditions after inflation are not necessarily zero. Specifically, a $\delta_{bc}(k)$ mode with $k\lesssim 0.001\ h/{\rm Mpc}$ generated during inflation will retain approximately the same amplitude at all times, much in the same way as the $k=0.001\ h/{\rm Mpc}$ mode in Fig. 7 remains small throughout. It is this case that our separate universe setup strictly applies to, since we generated the initial conditions for the simulations assuming constant modified $\Omega_{b}$ and $\Omega_{c}$ at all times up to the starting redshift of the simulation (corresponding to $a_{\rm ics}$ in the figure). The reason why this overestimates the impact on the amplitude of the initial power spectrum of the simulations for $\delta_{bc}$ modes generated by baryon-photon interactions is because these are still growing in between $a_{\rm eq.}$ and $a_{\rm dec.}$, which is when the modified relative abundances of baryons and CDM modify the amplitude of the power spectrum.

Appendix B Stellar mass resolution correction factors

In this appendix, we describe a stellar mass resolution correction scheme that we apply to the galaxies in our TNG300-2 simulations to test the robustness of our $b_{\delta}^{bc}(M_{*})$ measurements against numerical resolution. The numerical results we show here are for $z=1$, but they are representative of the other redshift values analysed in the paper.

The left panel of Figure 8 shows the ratio $M_{*}/M_{\rm h}$ as a function of $M_{\rm h}$ for the Fiducial cosmology at TNG100-1.5 (green) and TNG300-2 (blue) resolutions at $z=1$. The solid and dot-dashed lines indicate the corresponding median relation in bins of total halo mass. The observed difference between the TNG100-1.5 and TNG300-2 curves reflects the varying levels of convergence of $M_{*}$ at different resolutions. In Ref. [70], the authors demonstrate that a correction factor constructed using the stellar-to-halo-mass relations of two same-volume IllustrisTNG resolutions works well in bringing the results of simulations done at different volumes together in other stellar-mass-related quantities such as the stellar mass function or stellar mass radial profiles (cf. Appendix A of Ref. [70] for more details). Here, we apply the same resolution correction strategy to the TNG300-2 catalogues. Specifically, we define a stellar mass correction factor as

[TABLE]

where the numerator and denominator on the right-hand side are the median relations in the original TNG100-2 and TNG100-3 simulations, respectively; the superscripts ${}^{\text{oTNG100-2}}$ and ${}^{\text{oTNG100-3}}$ stress that we use the original simulations, which have matching phases of the initial conditions. It is important to note also that the TNG300-2 ($N_{p}=1250^{3}$, $L_{\rm box}=205\ {\rm Mpc/h}$) and TNG100-3 ($N_{p}=455^{3}$, $L_{\rm box}=75\ {\rm Mpc/h}$) simulations have approximately the same mass resolution. Our corrected stellar mass catalogues are then subsequently obtained by multiplying all of the stellar mass values of the TNG300-2 galaxies by $C^{\text{TNG300-2}}\left(M_{\rm h}\right)$; specifically, we interpolate over the values of $M_{\rm h}$ in Eq. (B.1), which is defined only at a finite number of total mass bins. We label the results of the corrected TNG300-2 catalogues as rTNG300-2, whose stellar mass values should be representative of a TNG100-2 resolution. While this still falls short of our higher resolution TNG100-1.5 results, it nonetheless allows us to check whether numerical convergence at the level of stellar masses plays a critical role in our measurements.

The dashed black line in the left panel of Fig. 8 shows the resulting median relation of the rTNG300-2 galaxy catalogues. Further, the middle panel of Fig. 8 shows the TNG300-2 stellar mass function measured with (dashed blue) and without (solid blue) the stellar mass correction. As expected, for both the stellar-to-halo-mass relation and the stellar mass function, the rTNG300-2 results are in closer agreement with TNG100-1.5, compared to TNG300-2. More importantly, the right panel of Fig. 8 shows the same as the $z=1$ panel of Fig. 3, but with the result of the corrected rTNG300-2 galaxies shown as well (dashed line). The $b_{\delta}^{bc}(M_{*})$ measured from the TNG300-2 and rTNG300-2 catalogues display only small differences with one another; in particular, both agree well with the higher-resolution TNG100-1.5 result. The discussion in the main body of the paper (in Sec. 3.2) about the TNG300-2 catalogues thus holds equally to the case of the corrected rTNG300-2 ones. This is not surprising since we have corrected the ${\rm Fiducial}$, High and Low cosmologies using the same correction factors, which preserve the relative difference between cosmologies that is effectively measured by $b_{\delta}^{bc}$; the observed small differences are caused by galaxies moving to different stellar mass bins.