Biased Tracers of Two Fluids in the Lagrangian Picture

TL;DR

This paper develops a Lagrangian perturbation theory framework for biased tracers in a two-fluid universe, quantifies corrections to bias and advection, and assesses their impact on BAO measurements, finding modest biases and manageable uncertainties.

Contribution

It introduces a novel two-fluid Lagrangian bias expansion, computes linear and nonlinear corrections, and evaluates their effects on BAO analysis.

Findings

Percent-level two-fluid corrections at low redshift.

Relative velocity effects can bias BAO scale by up to 0.5 sigma.

Fitting for two-fluid parameters increases BAO measurement errors by less than 30%.

Abstract

We explore Lagrangian perturbation theory (LPT) for biased tracers in the presence of two fluids, focusing on the case of cold dark matter (CDM) and baryons. The presence of two fluids induces corrections to the Lagrangian bias expansion and tracer advection, both of which we formulate as expansions in the three linear modes of the Lagrangian equations of motion. We compute the linear-order two-fluid corrections in the Zeldovich approximation, finding that modifications to the bias expansion and tracer advection both enter as percent-level corrections over a large range of wavenumbers at low redshift and draw parallels with the Eulerian formalism. We then discuss nonlinear corrections in the two-fluid picture, and calculate contributions from the relative velocity effect ($\propto v_r^2$) at one loop order. Finally, we conduct an exploratory Fisher analysis to assess the impact of two-fluid corrections on baryon acoustic oscillations (BAO) measurements, finding that while modest values of the relative bias parameters can introduce systematic biases in the measured BAO scale of up to $0.5\, \sigma$, fitting for these effects as additional parameters increases the error bar by less than $30\%$ across a wide range of bias values.

Equations

$${\text{x}}_{\sigma }(\text{q},\tau )=\text{q}+{\mathbf{\Psi }}_{\sigma }(\text{q},\tau ).$$

$$1+{\delta }_{\sigma }(\text{x},\tau )=\int d^{3}q{\delta }_D(\text{x}-\text{q}-{\mathbf{\Psi }}_{\sigma }(\text{q},\tau ))=\int d^{3}q\frac{d^{3}k}{(2\pi {)}^3}{e}^{i\text{k}\cdot (\text{x}-\text{q}-{\mathbf{\Psi }}_{\sigma }(\text{q},\tau ))},$$

$${\ddot{\mathbf{\Psi }}}_c+\mathcal{H}{\dot{\mathbf{\Psi }}}_c$$

$${\ddot{\mathbf{\Psi }}}_b+\mathcal{H}{\dot{\mathbf{\Psi }}}_b$$

$${\nabla }_x^2\Phi (\mathbf{x},\tau )=\frac{3}{2}{\Omega }_m(\tau ){\mathcal{H}}^2(\tau ){\delta }_m(\mathbf{x},\tau ),$$

$${\delta }_{\sigma }({\text{x}}_{\sigma })\leftrightarrow -\nabla \cdot {\mathbf{\Psi }}_{\sigma }(\text{q}),{\text{v}}_{\sigma }(x_{\sigma })\leftrightarrow {\dot{\mathbf{\Psi }}}_{\sigma }(\text{q}).$$

$${\ddot{\mathbf{\Psi }}}_m+\mathcal{H}{\dot{\mathbf{\Psi }}}_m$$

$${\ddot{\mathbf{\Psi }}}_r+\mathcal{H}{\dot{\mathbf{\Psi }}}_r$$

$${\mathbf{\Psi }}_m(\tau )=-{\text{m}}_{+}{D}_{+}(\tau )+{\text{m}}_{-}{D}_{-}(\tau )\approx -{\mathbf{m}}_{+}{D}_{+}(\tau ),$$

$$\dot{\boldsymbol{\Psi}}_{r}(\tau)=\textbf{v}_{r}(\tau_{i})\,\Big{(}\frac{a_{i}}{a}\Big{)}+\frac{1}{a}\int_{\tau_{i}}^{\tau}d\tau^{\prime}\ a(\tau^{\prime})\textbf{F}_{b}(\tau^{\prime}),$$

$$\dot{\boldsymbol{\Psi}}_{r}(\tau)=\dot{\boldsymbol{\Psi}}_{r}(\tau_{i})\frac{a(\tau_{i})}{a(\tau)}+\Bigg{[}\frac{1}{a}\int_{\ln(a(\tau_{i}))}^{\ln(a(\tau))}d\ln(a^{\prime})\Bigg{(}n_{e}(a^{\prime})\sigma_{T}\frac{\rho_{\gamma}(a^{\prime})}{\rho_{b}(a^{\prime})}\Bigg{)}\,\frac{f(a^{\prime})D_{+}(a^{\prime})}{a^{\prime 2}}\Bigg{]}\frac{\boldsymbol{\Psi}_{m}(\tau_{i})}{D_{+}(\tau_{i})}\;,$$

$${\mathbf{\Psi }}_r(\tau )=-{\mathbf{r}}_{+}+{\mathbf{r}}_{-}{D}_r(\tau ,{\tau }_i)+{\mathbf{m}}_{+}{D}_{\mathrm{CD}}(\tau ,{\tau }_i),D_r(\tau ,{\tau }_i)={\int }_{{\tau }_i}^{\tau }\frac{H_{0}d{\tau }^{\prime }}{a({\tau }^{\prime })}$$

$$T_{{\delta }_r}(k)\equiv T_{{\delta }_b}(k)-T_{{\delta }_c}(k)\mathrm{and}{T}_{{\theta }_r}(k)\equiv T_{{\theta }_b}(k)-T_{{\theta }_c}(k)$$

$$T_{\nabla \cdot {\mathbf{r}}_{+}}(k)$$

$$T_{\nabla \cdot {\mathbf{r}}_{-}}(k)$$

$$T_{\nabla \cdot {\mathbf{\Psi }}_r}(k,z)=T_{\nabla \cdot {\mathbf{r}}_{+}}(k)+D_r(z,z_i)T_{\nabla \cdot {\mathbf{r}}_{-}}(k)+D_{\mathrm{CD}}(z,z_i)T_{\nabla \cdot m_{+}}(k).$$

$$1+{\delta }_g(\text{x},\tau )=\int d^{3}q{F}_g[\mathbf{q}|\{{\mathbf{\Psi }}_{\sigma }(\text{q})\}]{\delta }_D\left[\text{x}-\text{q}-{\mathbf{\Psi }}_g(\text{q},\tau )\right].$$

$$\mathbf{q}\to \mathbf{q},{\mathbf{\Psi }}_{\sigma }\to {\mathbf{\Psi }}_{\sigma }+\text{n}(\tau ),\Phi \to \Phi \to \Phi -\text{x}\cdot (\ddot{\text{n}}+\mathcal{H}\dot{\text{n}}),$$

$$F_g=F_g\left[{\delta }_{\sigma },{\text{v}}_{\sigma },\Phi \right]\equiv F_g\left[{\delta }_m,{\delta }_r,\partial {\text{v}}_m,{\text{v}}_r,\partial \partial \Phi ,\cdots \right],$$

$$1+{\delta }_g(\mathbf{q})=1+b_1{\delta }_m+b_r{\delta }_r+b_{\theta }{\theta }_r+\cdots$$

$$F_g(\mathbf{q})=b_1{\delta }_m+b_{+}\nabla \cdot {\text{r}}_{+}+b_{-}\nabla \cdot {\text{r}}_{-}+...$$

$${\ddot{\mathbf{\Psi }}}_g+\mathcal{H}{\dot{\mathbf{\Psi }}}_g=-\nabla \Phi +{\mathbf{F}}_{b,g},$$

$${\mathbf{\Psi }}_g(\tau )={\mathbf{\Psi }}_m(\tau )+f_g[{\mathbf{\Psi }}_r(\tau ){]}_{\mathrm{CD}=0},$$

$$P_{gg}(k,z)=\int d^{3}q{e}^{i\mathbf{k}\cdot \mathbf{q}}{⟨F_g[{\mathbf{q}}_1]F_g[{\mathbf{q}}_2]e^{ik\cdot ({\mathbf{\Psi }}_g({\mathbf{q}}_1,z)-{\mathbf{\Psi }}_g({\mathbf{q}}_2,z))}⟩}_{q=|{\mathbf{q}}_2-{\mathbf{q}}_1|},$$

$$e^{i\mathcal{M}}\equiv ⟨\exp \left(i\mathbf{k}\cdot {\mathbf{\Delta }}_g(z)+{\lambda }_{{\delta }_{m,1}}{\delta }_{m,1}+{\lambda }_{+,1}\nabla {\mathbf{r}}_{+,1}+{\lambda }_{-,1}\nabla {\mathbf{r}}_{-,1}+(1\leftrightarrow 2)\right)⟩$$

$${\mathbf{\Delta }}_g={\mathbf{\Psi }}_{g,1}-{\mathbf{\Psi }}_{g,2}={\mathbf{\Psi }}_{m,1}-{\mathbf{\Psi }}_{m,2}+f_g({\mathbf{\Psi }}_{r,1}-{\mathbf{\Psi }}_{r,2})\equiv {\mathbf{\Delta }}_m+f_g{\mathbf{\Delta }}_r$$

$$e^{i\mathcal{M}}=\exp \{$$

$$\displaystyle+ik\cdot\big{(}(\lambda_{\delta_{m,1}}+\lambda_{\delta_{m,2}})(U_{mm}+f_{g}U_{rm})$$

$$\displaystyle+(\lambda_{+,1}+\lambda_{+,2})(U_{m+}+f_{g}U_{r+})+(\lambda_{-,1}+\lambda_{-,2})(U_{m-}+f_{g}U_{r-})\big{)}$$

$$+({\lambda }_{{\delta }_{m,1}}{\lambda }_{+,2}+(1\leftrightarrow 2)){\xi }_{{\delta }_m\nabla {\mathbf{r}}_{+}}+({\delta }_m,\nabla {\mathbf{r}}_{-})+({\delta }_m,{\delta }_m)$$

$$+(\nabla {\mathbf{r}}_{+},\nabla {\mathbf{r}}_{+})+(\nabla {\mathbf{r}}_{+},\nabla {\mathbf{r}}_{-})+(\nabla {\mathbf{r}}_{-},\nabla {\mathbf{r}}_{-})\},$$

$$A_{ij}^{ab}=⟨{\Delta }_i^a(z){\Delta }_j^b(z)⟩,U_i^{a\pm }=⟨{\Delta }_i^a(z)\nabla \cdot {\mathbf{r}}_{\pm }({\mathbf{q}}_1)⟩,{\xi }_{ab}=⟨a({\mathbf{q}}_1)b({\mathbf{q}}_2)⟩,$$

$$({\delta }_m,\nabla {\mathbf{r}}_{-})\equiv \left({\lambda }_{{\delta }_{m,1}}{\lambda }_{-,2}+(1\leftrightarrow 2)\right){\xi }_{{\delta }_m\nabla {\mathbf{r}}_{-}}(q)$$

$$P_{gg}(k)=\int d^{3}q{e}^{i\mathbf{k}\cdot \mathbf{q}}{e}^{-\frac{1}{2}{k}_i{k}_j{A}_{ij}^{mm}}$$

Full text

1 Introduction

Observations of the large-scale structure (LSS) of the universe allow us to shed light on areas of physics ranging from galaxy formation and evolution to fundamental physics. A prime target of present and future LSS surveys is the measurement of baryon-acoustic oscillations (BAO) – the imprints of sound waves in the baryon-photon fluid observed in the cosmic microwave background (CMB) on the observed clustering of galaxies – which can be used as a standard ruler to constrain the expansion of the universe [1]. Upcoming surveys such as DESI [2], EUCLID [3] and WFIRST [4] will provide BAO measurements with higher-than-ever precision, and even more futuristic BAO surveys such as a Stage II 21-cm experiment [5] have been proposed. These next-generation observational campaigns will require us to model the LSS with unprecedented accuracy, at the sub-% level.

One area of recent interest in the field of LSS has been in accounting for the effects induced by the existence of multiple species (cold dark matter, baryons, neutrinos), with similar but distinct clustering properties, using analytic methods. Studies of the perturbative approach to structure formation have traditionally grouped all nonrelativistic species into a “total matter” fluid, whose gravitational collapse is the dominant source of structure on cosmological scales in the late-time universe, but many authors have recently extended these techniques to include neutrinos [6, 7, 8, 9, 10] and baryons [11, 12, 13, 14, 15, 16, 17] in the Eulerian framework of Standard Perturbation Theory (SPT). In parallel, the response of galaxy and halo formation to the existence of multiple fluid species has also been subject of extensive investigation [11, 18, 19, 20, 21, 22, 16, 17]. Of particular interest are the present-day imprint of relative perturbations between baryons and dark matter on large scales which, being seeded in the same epoch and at the same scales as the baryon acoustic oscillations, has the potential to confound future BAO measurements [11, 22, 23, 24]. While these relative perturbations do not grow significantly in time (and relative velocities in fact decay) and are thus small compared to the total-matter growing mode at late times, they amount to coherent supersonic flows post-recombination and could have significant effects on the formation of the first halos and galaxies [13, 11], which are the progenitors of the objects we observe today.

The goal of this work is to formulate perturbation theory and galaxy bias in the presence of multiple fluids within the Lagrangian framework, with a particular focus on the two-fluid baryon-dark matter scenario. Our work is a direct extension of the aforementioned SPT calculations. While Lagrangian Peturbation Theory (LPT) is order-by-order equivalent to SPT, it seamlessly allows a consistent treatment of large scales bulk flows, which are responsible for the final shape and position of the BAO features in the correlation functions or power spectrum [25, 26, 27, 28, 29, 30, 31]. The theory can also be extended to handle density field ‘reconstruction’ [32, 33, 34, 35, 36]. These features make LPT a natural language for investigating possible distortions to the BAO feature.

This paper is organized as follows. In Section 2, we introduce the linear Lagrangian equations of motion and discuss the role of non-gravitational forces such as Compton drag with the CMB. Modifications to Lagrangian galaxy bias and advection in the two-fluid limit are then introduced in Section 3. In Section 4, we employ the results of the preceding two sections and calculate the lowest-order two-fluid corrections to the galaxy power spectrum in the Zeldovich approximation. Cross spectra and subtleties in the IR resummation are briefly discussed in Section 4.2. In Section 5 we take up whether the calculated two-fluid corrections can significantly bias BAO measurements, arguing that any such biases can be mitigated by simultaneously fitting for these easily-characterizeable effects. Our conclusions are summarized in Section 6.

2 Linear Equations of Motion in Lagrangian Space

In the Lagrangian picture, fluid dynamics is encoded in the displacements $\Psi_{\sigma}(q)$ of fluid elements of each species, $\sigma$, originally situated at Lagrangian positions q, such that their Eulerian positions at conformal time $\tau$ ($d\tau=a^{-1}dt$) are given by [37, 28, 38]

[TABLE]

The subscript $\sigma=\{c,b\}$ denotes the species, either cold dark matter (CDM) or baryons, respectively, whose motion we are tracking. Assuming that initial displacements are infinitesimally small compared to those at the redshifts of interest, the overdensity, $\delta_{\sigma}$, of each species at Eulerian position $x$ can be solved for via mass conservation

[TABLE]

where $\delta_{D}$ is the Dirac delta function. Taylor expanding to first order in displacements yields the familiar result that $\delta_{\sigma}(x)=-\nabla\cdot\boldsymbol{\Psi}_{\sigma}(q)$, but, as seen in Equation 2.2, one feature of working in the Lagrangian picture is that the translation into Eulerian quantities, such as the density field, invariably involves nonlinear combinations of $\boldsymbol{\Psi}$ even when only the linear equations of motion are considered.

2.1 General Formalism

While CDM particles by assumption experience only the gravitational force, baryons are subject to non-gravitational effects, such as Compton drag and pressure gradients. These effects can be summarized in the equations of motion of the fluid elements

[TABLE]

where overdots signify derivatives with respect to $\tau$, $\mathcal{H}=d\ln a/d\tau$ is the conformal Hubble parameter, $\mathbf{F}_{b}$ is the non-gravitational force per unit mass felt by baryons, and $\Phi$ is the gravitational potential at Eulerian position $x$ satisfying Poisson’s equation

[TABLE]

where $\Omega_{m}$ is the total matter mass density and $\delta_{m}$ is the total matter overdensity (see below).

At the linear level, there is no difference between the Eulerian and Lagrangian positions in the above equations of motion, and we will neglect this distinction in the rest of this section unless otherwise stated. Indeed, taking the divergence of Equation 2.3 in the linear limit ($x_{\sigma}\approx q$) directly yields the Euler equation when we map overdensities to displacements and velocities to their derivatives:

[TABLE]

Note that the first mapping is correct only to linear order, while the second one is exact if the full $\mathbf{x}(\mathbf{q})$ is used. Assuming this translation, the solutions to the Lagrangian equations of motion as described below are essentially identical to those extracted from Boltzmann codes such as CAMB [39] or CLASS [40], provided one chooses post-recombination initial conditions for the Lagrangian displacements.

To solve Equation 2.3 in the linear limit, it is convenient to rewrite the baryonic and CDM displacements in terms of a mass-weighted matter component ($\boldsymbol{\Psi}_{m}=w_{c}\boldsymbol{\Psi}_{c}+w_{b}\boldsymbol{\Psi}_{b}$), which sources the gravitational potential, and a relative component that characterizes the differential flows between baryons and CDM ($\boldsymbol{\Psi}_{r}=\boldsymbol{\Psi}_{b}-\boldsymbol{\Psi}_{c}$), where we have defined the mass fractions of each species, $w_{\sigma}=\rho_{\sigma}/\rho_{m}$. These are related to the Eulerian quantities $\delta_{m}=w_{b}\delta_{b}+w_{c}\delta_{c}$ and $v_{r}=v_{b}-v_{c}$ by $\delta_{a}=-\nabla\cdot\boldsymbol{\Psi}_{a}$ and $v_{a}=\dot{\boldsymbol{\Psi}}_{a}$, where $a=\{m,r\}$, again at the linear level. The equations of motion in terms of these components are

[TABLE]

If in addition non-gravitational forces are negligible, the matter and relative components decouple, such that Equation 2.6a can be solved as

[TABLE]

where $D_{+}$ is the usual linear-theory growth factor. In the last step we have neglected the decaying mode, $\textbf{m}_{-}$, since it is a tiny fraction of the total displacement at all redshifts of interest. For non-gravitational forces, like Compton drag or pressure gradients, direct integration of the linear equations of motion show that the non-gravitational terms make a negligible contribution to the matter component $\boldsymbol{\Psi}_{m}$, such that the transfer function at redshifts below $z=6$ agree with the linear solution in Equation 2.7 to within $0.2\%$, with even better agreement at the lower redshifts of interest in this paper. In the above we have included a minus sign for convenience such that $\delta_{m,0}=\nabla\cdot\textbf{m}_{+}$.

We end this subsection by discussing the full solution of the relative displacement when $\textbf{F}_{b}=\textbf{F}_{b}(\tau)$ is independent of $\boldsymbol{\Psi}_{r}$. In this case Equation 2.6b is linear and first order in $\dot{\boldsymbol{\Psi}}_{r}$ and can be solved as:

[TABLE]

where we have set the boundary conditions at initial time $\tau_{i}$ assuming the non-gravitational effects encoded in $F_{b}$ do not turn on until $\tau>\tau_{i}$. Equation 2.8 turns out to be an excellent approximation for the large-scale Compton drag electrons experience in the reionization era due to their relative motion with respect to the CMB rest frame, $\mathbf{F}_{b}=-n_{e}\sigma_{T}(\rho_{\gamma}/\rho_{b})a\mathbf{v}_{b}$, where $\sigma_{T}$ is the Thompson scattering cross section, $\rho_{\gamma}$ is the photon energy density and $n_{e}$ the free electron number density. Eq. (2.8) also applies baryonic pressure forces on small scales $\mathbf{F}_{b}\propto-\nabla\delta_{b}$— in both cases the total-matter component may be substituted for the baryonic component (i.e. $\delta_{b}\approx\delta_{m}$) at the sub-percent level [17]. In the case of the large-scale Compton drag, assuming $\mathbf{v}_{b}\simeq\mathbf{v}_{m}$ yields

[TABLE]

with $f=dD_{+}/d\ln(a)$ the linear theory growth factor. The Compton drag thus induces a mixing between the matter and relative components through a numerical prefactor dependent only on the linear growth factor $D_{+}$ and reionization history via $n_{e}$. Finally, we can integrate 2.9 to yield

[TABLE]

where we can identify $\boldsymbol{\Psi}_{r}(\tau_{i})=-\mathbf{r}_{+}$, $a(\tau_{i})\mathbf{v}_{r}=H_{0}\mathbf{r}_{-}$, and the Compton-drag kernel $D_{\rm CD}$ is defined as the conformal time integral of the square-bracketed function in 2.9. The linear solutions to both the total-matter and relative components are thus wholly specified by the three modes $\mathbf{m}_{+}$ and $\mathbf{r}_{\pm}$. Jeans instabilities and baryonic pressure forces affect much smaller scales and won’t be further discussed in the remainder of this work.

2.2 Initial conditions and transfer functions

The linear evolution of the density and velocity contrasts can be easily written in terms the CDM and baryon linear transfer functions (output from, e.g. CAMB) as

[TABLE]

where $\theta_{b,c}(k)\equiv-\dot{\delta}_{b,c}(k)$. It is worth noticing that while the velocity field is gauge dependent, velocity differences are not. The transfer function for $\nabla\cdot\mathbf{m}_{+}$ is simply the present-day matter transfer function $T_{m}$ and we can furthermore define

[TABLE]

These three functions specify the solution for the $\boldsymbol{\Psi}_{m}$, $\boldsymbol{\Psi}_{r}$ and $\dot{\boldsymbol{\Psi}}_{r}$ at any $z<z_{i}$. The choice of $z_{i}$ is somewhat arbitrary but choosing redshifts before the onset of reionization has the advantage of separating the effects of gravity from Compton drag. This choice also justifies the normalization in Eq. (2.12), since $\mathbf{r}_{-}$ is independent of redshift. In the remainder of the paper we assume $z_{i}=20$.

In addition to the above, we will show below that calculating the power spectrum at some redshift $z$ in the Lagrangian picture requires linear-theory spectra of the relative displacement at that redshift, which will typically include corrections from Compton drag. These can be calculated via Equations 2.7 and 2.8 to give

[TABLE]

Sample solutions of the equation of motion in Eq. (2.6) when $F_{b}$ is given by Compton drag with the CMB are shown in Figure 1. After reionization most of large scale power in the relative velocity transfer function, $T_{\nabla\cdot\dot{\boldsymbol{\Psi}}_{r}}$, is provided by the Compton drag, which in turn affects the evolution of the relative baryon-dark matter density at large scales (see top panels in Figure 1). Figure 1 also justifies the approximations we used to compute the drag forces, as one can see by the excellent agreement with the full CAMB output. Other non-gravitational effects like pressure terms (Jeans instability) and radiative transfer effects [41, 42, 43, 44], can be written in a similar form.

Ratios of the transfer functions to the total matter one are shown in Figure 2. We notice that the relative density perturbation is much larger than the relative velocity one, by a factor of a hundred at least, and the two relative components have the same behavior with wave-number $k$ at small and large scales. Nonetheless $\mathbf{r}_{+}$ and $\mathbf{r}_{-}$ have significant differences in shape around the BAO scales and therefore will have to be treated separately from the point of view of the galaxy bias expansion.

3 Lagrangian Bias in the Two-Fluid Dynamics

In the Lagrangian approach, galaxy bias is assumed to arise as the response of the overdensity of galaxies, or the precursors thereof, to the variation of the initial conditions encoded in the fields $\{\boldsymbol{\Psi}_{\sigma}(\textbf{q})\}$ of the various species, and then transported via advection to their present-day positions $\textbf{x}(q,t)=\textbf{q}+\boldsymbol{\Psi}_{g}(\textbf{q},t)$. Thus, when computing the density of a biased tracer the number-conservation Equation 2.2 is modified to

[TABLE]

The standard picture of (local) Lagrangian bias, outlined above, has been developed in the 1-fluid case by many authors, see for example [45, 46, 47, 48, 49, 50, 51, 52, 53] and [54] for a recent review on galaxy bias. In this section our focus will be on extending these arguments to the case of multiple fluids, and in particular to the two-fluid case. In the presence of two fluids, the form of Equation 3.1 raises two questions: (1) the form of the response $F_{g}$ and (2) whether biased tracers follow the dark matter, baryons, or a combination thereof. We address these in turn.

3.1 Bias Expansion

The initial tracer overdensity, $F_{g}[\mathbf{q}|\,\{\boldsymbol{\Psi}_{\sigma}(\textbf{q})\}]$, is defined to be a functional encoding the physics of gravitational collapse and galaxy formation at some Lagrangian position $\mathbf{q}$. Since the galaxy density field is a scalar quantity under rotations, $F_{g}$ will also be a scalar. We will assume this functional is local, in the sense gravitational collapse depends only on the value of the fields within a characteristic patch of size $R_{h}$, which then flows coherently on large scales with $\boldsymbol{\Psi}_{g}$ [54]. In the fluid limit, these conditions imply that the system is wholly specified – albeit in some complex, nonlinear way – by the species overdensities, $\delta_{\sigma}(\textbf{q})$, velocities, $\textbf{v}_{\sigma}(\textbf{q})$, and the gravitational potential111The gravitational potential $\Phi$, while not independent of $\delta_{m}$, depends on the total matter density in a very non-local way. To make our bias expansion local, and be able to truncate the derivative expansion at a reasonable order, we thus include it as a standalone quantity here., $\Phi(\textbf{q})$, at some initial time $\tau_{i}$. The condition that $F_{g}$ is local – or rather, nonlocal with width $R_{h}$ – can be equivalently (but more conveniently) expressed by requiring $F_{g}$ depend only on the initial fields and their spatial derivatives, with $n^{\rm th}$ derivatives suppressed by $n$ powers of $R_{h}$ [54].

In addition to the assumption of locality, the form of $F_{g}$ is strongly restricted by various symmetries. General relativity requires that all physical quantities be diffeomorphism invariant, which in our case reduces to rotational invariance and invariance under generalized Galilean transformations [55]:

[TABLE]

where $\mathbf{n}$ are time-dependent but spatially constant vector fields.

Rotational invariance simply requires that only contracted scalar quantities enter the bias; the restrictions placed on the bias expansion by general Galilean invariance are more subtle, and it is here that the two-fluid case diverges from the single-fluid case. Under this symmetry, densities remain unchanged— for instance that at first order $\delta_{\sigma}(\mathbf{q})=-\nabla\cdot\boldsymbol{\Psi}_{\sigma}(\mathbf{q})$— while velocities get boosted by a spatially constant amount (leaving $\partial\textbf{v}$ invariant) and the gravitational potential changes in a spatially linear way (leaving $\partial\partial\Phi$ invariant). In the single-fluid regime, where only one set of densities and velocities exist, this directly implies that velocities can only enter with at least one spatial derivative, and the gravitational potential can only enter as second (spatial) derivatives and beyond. The single-fluid overdensity, which is unchanged under the transformation, can enter at any order.

The presence of two or multiple fluids relaxes some of the above restrictions. In particular, since all particle velocities are boosted by the same amount ($\textbf{n}^{\prime}$) under a Galilean transformation, the relative velocity $\textbf{v}_{r}=\textbf{v}_{b}-\textbf{v}_{c}$ remains invariant and can thus enter the bias expansion at zeroth order in derivatives. The total matter velocity, $\textbf{v}_{m}$, on the other hand, is boosted and can thus still only enter at the derivative level. These two quantities form an equivalent basis to the individual species velocities and there is no loss of generality in defining the bias expansion in terms of them. We may similarly write terms involving species densities, which can enter separately, in the total matter and relative density basis. In general relativity the gravitational potential is unaffected by the number of species as a consequence of the equivalence principle, i.e. gravitational interactions are universal. The full set of physical fields that can enter $F_{g}$ in the two fluid case is then

[TABLE]

where the dots stand for higher derivative operators. To first order in the fields we can therefore write222A list of bias terms up to second order is given in Appendix B.2.

[TABLE]

which is similar to the Eulerian linear theory expression in [16]. This is not surprising, since at first order $\mathbf{q}\simeq\mathbf{x}$, however we will see below that differential advection can introduce further terms degenerate with the initial Lagrangian bias terms above, such that the Eulerian relative-component bias will in general be a combination of these terms.

Finally, since $F_{g}$ is defined as a functional on the initial conditions which can be chosen to be sufficiently early that they lie deep in the linear-theory regime, its form can be further simplified and expressed purely in terms of the initial modes $\textbf{m}_{+}$ and $\textbf{r}_{\pm}$. In the single fluid case, this restriction leads to the simplification that all bias terms can be written in terms of spatial derivatives of the total matter displacement $\textbf{m}_{+}\sim\boldsymbol{\Psi}_{m}$; this is a direct consequence that, up to time-dependent constant factors, $\delta_{m}\sim\partial\boldsymbol{\Psi}_{m}$, $\textbf{v}_{m}\sim\boldsymbol{\Psi}_{m}$ and $\partial\partial\Phi\sim\partial\boldsymbol{\Psi}$ in linear theory. In the two-fluid case these terms must be supplemented by those involving the relative modes. Specifically, including the $\textbf{v}_{r}$ dependence requires the inclusion of terms proportional to $\textbf{r}_{-}$ and including $\delta_{r}$ dependence similarly requires terms proportional to $\nabla\cdot\textbf{r}_{+}$. Equation 3.4 can thus be re-expressed as:

[TABLE]

We therefore have a direct correspondence in the bias expansion between the initial modes expressed in Eulerian and Lagrangian space. Notice that the bias expansion defined above is complete, in the sense that it contains all possible operators compatible with the symmetries of the problem. In particular, while $\mathbf{r}_{\pm}$ are defined at a particular initial redshift $z_{i}$, in the linear regime this dependence amounts to a simple linear transformation and can be absorbed into the definition of the bias parameters (Appendix A).

Finally, an additional complication arises when halo formation is affected by Compton drag. As pointed out by [17], by picking out the local CMB rest frame such that the drag force $\propto\mathbf{v}_{b}$, we lose the gauge redundancy of Galilean transformations. This will in general produce heretofore forbidden terms such as those proportional to the matter-component velocity $\mathbf{v}_{m}$. However, the terms thus generated are required by rotational invariance to enter at second order and beyond. For the remainder of this paper we will thus neglect these contributions, which are subdominant to the already sub-percent level contributions we study.

Whereas there exists quite a large literature on measuring and predicting, using approximate physical models, the value of the bias parameters in one-fluid scenarios, less attention has been devoted to the multi-fluid case. From an effective field theory perspective the dimensionless parameters should be of order unity, but in reality the actual value of the bias parameters is tracer-dependent and can be quite a bit larger or smaller. In this work we will assume, unless otherwise noted, that typical values are given by $b_{+}\simeq 1$ and $b_{-}\simeq 6.8$ derived in [16] using a spherical collapse model. These numbers are consistent with the non-detection of relative bias effects in BOSS DR12 by ref. [23], who find e.g. $b_{+}=-1.0\pm 2.5$ to within one sigma when fitting for $b_{+}$, $b_{-}$ and $c_{-}$ (Section 4.3) across all redshift bins, with large systematic biases measured in dark-matter only simulations that had to be subtracted.

3.2 Modifications to Tracer Advection

Once the initial, biased tracer overdensity is set, the overdensity at later times is set by the tracer “fluid” advecting from initial (Lagrangian) $q$ to final positions $q+\boldsymbol{\Psi}_{g}$ along trajectories described by the tracer equation of motion

[TABLE]

where we have included a non-gravitational term, $F_{b,\,g}$, to account for the possibility that tracers feel non-gravitational forces. Such non-gravitational contributions may arise, for example, from the Compton drag on the baryonic component of galaxies, or from various galaxy formation processes. Since such contributions are always local in space and time, we will assume the above equation satisfies the same symmetries of Eq. (3.2), i.e. the force acting on galaxies depends only on density fields and velocity gradients.

Equation 3.6 can be solved by subtracting the equation of motion of the total matter displacement (Equation 2.6a) and defining $\boldsymbol{\Psi}_{r,g}=\boldsymbol{\Psi}_{g}-\boldsymbol{\Psi}_{m}$. Neglecting the baryonic contributions such that the tracers’ dynamics are governed only by gravity, and assuming that the initial tracer displacements are a weighted average of the baryonic and CDM displacements, i.e. $\boldsymbol{\Psi}_{g,i}=\boldsymbol{\Psi}_{m,i}+f_{g}\boldsymbol{\Psi}_{r,i}$, this immediately yields the time evolution

[TABLE]

where the relative displacement is evaluated assuming zero Compton drag. Note that if we assume that the tracer field is made of objects composed of the same mass fractions of baryons and CDM as the total matter content of the universe, i.e. with $f_{g}=0$, Equation 3.7 reduces to the trajectory of the matter component. Similarly, objects composed purely of baryons or the CDM will (at the linear level) follow the baryon or CDM displacements, respectively.

We can alternatively think of Eq. (3.7) as a bias expansion of the galaxy displacements in terms of the underlying fields, since $\boldsymbol{\Psi}_{m}$ and $\boldsymbol{\Psi}_{r}$ are the only two linear operators allowed by symmetries at lowest order in spatial derivatives. If the tracer flow is purely gravitational, as assumed above, the equivalence principle further restricts the coefficient of the total matter displacement – which encapsulates the motion due to the gravitational potential – to be exactly 1 at all times. However, this restriction can be broken by baryonic contributions ($\propto\mathbf{F}_{b,\,g}$) such as the Compton drag. As seen in the second term on the right hand side of Equation 2.9, the acceleration due to Compton drag generates displacements proportional to $\boldsymbol{\Psi}_{m}$; this contribution, on top of the aforementioned gravitational displacements, can lead to an expansion $\boldsymbol{\Psi}_{g}=(1+\alpha_{\rm CD})\boldsymbol{\Psi}_{m}+f_{g}\boldsymbol{\Psi}_{r}+...$ for some nonzero coefficient $\alpha_{\rm CD}$ due to Compton drag, where the total-matter coefficient deviates from unity. Consequences of this modified expansion for the power spectrum are considered at the end of Section 4.1 and in Figure 6. Other baryonic forces, such as pressure forces at small scales, can similarly be included as further terms ($\boldsymbol{\Psi}_{g}\ni c_{s}^{2}\nabla\delta_{b}$) in this expansion.

4 Galaxy Power Spectra in the Zeldovich Approximation

4.1 Analytic Form

From Equation 3.1, the power spectrum at redshifts $z$ for a biased tracer can be computed as

[TABLE]

where the subscripts denote quantities evaluated at two points separated by $q$ in Lagrangian space. It is important to note that the bias functions $F_{g}$ are evaluated in terms of the linear modes $m_{+},r_{\pm}$ defined at the initial redshift $z_{i}$. In the Zeldovich approximation displacements are solved to linear order but the full mapping between initial and final times is kept. This amounts to keeping the displacement correlators exponentiated in what follows [46]. We will adopt the bias expansion in Equation 3.5. We evaluate integrals involving $F_{g}$ by functional differentiation in the usual manner [46, 47, 56]: we include a term (e.g. $\lambda X$) in the exponential for each argument, $X$, of $F_{g}$ and evaluate terms like $X^{n}$ via $\partial^{n}/\partial\lambda^{n}$ of $\exp[\lambda X]$.

Under the above assumptions our task reduces to evaluating

[TABLE]

with numerical subscripts referring to Lagrangian coordinates, $q_{1}$ and $q_{2}$, and

[TABLE]

The function $e^{i\mathcal{M}}$ can be evaluated using the cumulant theorem as the exponential of the connected components. The Zeldovich approximation assumes linear dynamics, such that only quadratic terms survive

[TABLE]

where we have defined

[TABLE]

noting that the $\Delta$’s carry an implicit redshift dependence while the other fields do not. For the total-matter component this redshift dependence is a direct growth factor scaling and we will for convenience take the linear field’s value as evaluated at the observed redshift $\delta_{m}=-D_{m}(z)\nabla\cdot\mathbf{m}_{+}$. The paired parentheses denote terms similar to the preceding except with the indicated pair of variables. For example, in the third line

[TABLE]

and when the elements of a pair are repeated the term should be divided by a symmetry factor of two.

Figure 3 shows the different correlation functions entering the above calculation. Since the correlation function of the different displacements fields, $A^{ab}_{ij}(q)$, is a tensor, we can decompose it as $A^{ab}_{ij}(q)=X^{ab}(q)\delta^{K}_{ij}+Y^{ab}(q)\hat{q}_{i}\hat{q}_{j}$, and the functions $X(q)$’s and $Y(q)$’s are shown in the left panel of Figure 3. Clearly the galaxy displacements are dominated by the total matter component, with the relative terms contributing much less than a % to the bulk flows. This fact will enable us to treat the terms proportional to $f_{g}$ perturbatively, as they will be much smaller than one for wavenumbers below the nonlinear scale defined by $k^{2}\Sigma^{2}\lesssim 1$, where the Zeldovich r.m.s. displacement is $\Sigma\propto X_{mm}(q\rightarrow\infty)$. The same conclusions apply for the bias weighted displacements $U(q)$’s, shown on the right hand panels in Figure 3, where $U_{m}(q)\gg U_{\pm}(q)$.

Working to linear order in the power spectrum we then have that the galaxy-galaxy power spectrum is given by

[TABLE]

Figure 4 shows the different contributions to the galaxy power spectrum in the Zeldovich approximation at $z=1.2$. The leading corrections to the total-matter power spectrum come at the roughly percent level from terms in Equation 4.7 linear in $\mathbf{r}_{+}$, i.e. in $b_{+}$ and $f_{g}$. These contributions are essentially degenerate, with differences due to the dynamical evolution of $\boldsymbol{\Psi}_{r}$ in the $f_{g}$ term, as we will discuss in the next paragraph. Corrections quadratic in $\mathbf{r}_{+}$ or linear in $\mathbf{r}_{-}$ enter at roughly the same size four orders of magnitude below the total-matter contributions.

An interesting consequence of the advection of biased tracers with $|f_{g}|>0$ is the appearance of relative bias terms even if none were present in the initial Lagrangian bias expansion. To see this, we can take the low-$k$ limit of Eq. (4.7), neglecting for the moment non-gravitational contributions to $\boldsymbol{\Psi}_{r}(q)$, and obtain up to $\mathcal{O}(P(k))$

[TABLE]

We immediately recognize the familiar expression for the Eulerian linear bias, $b^{E}_{1}=1+b_{1}$, and that the relative density and velocity bias terms get renormalized by terms proportional to $f_{g}$. To make further contact with the existing literature employing the Eulerian formulation of the equations of motion [16, 17], we can identify the relative baryon dark-matter density perturbation $\delta_{r}$ with the divergence of $\mathbf{r}_{+}$, $\delta_{r}\equiv\nabla\cdot\mathbf{r}_{+}$, and the relative baryon dark-matter velocity divergence $\theta_{r}$ with the divergence of $\mathbf{r}_{-}$, $\theta_{r}\equiv(1+z)H_{0}\nabla\cdot\mathbf{r}_{-}$. This implies that the bias parameters in [16, 17] associated to the Eulerian fields are $b^{E}_{\delta_{r}}=b_{+}+f_{g}$ and $b^{E}_{\theta_{r}}=(1+z)^{-1}H_{0}^{-1}(b_{-}+f_{g}D_{r}(z))$. Note that the referenced overdensities and velocities are those defined at the initial redshift $z_{i}$ so should not be directly substituted for their Eulerian counterparts; for more details about the mapping of bias parameters from some initial time $z_{i}$ to Eulerian coordinates see Appendix A.

A final caveat occurs when the non-gravitational forces on the tracer, $\mathbf{F}_{b,g}$ are nonzero. The integrated effect of such forces on $\boldsymbol{\Psi}_{r,g}$ must then be accounted for. For example, when dealing with baryons and dark matter, the effects of Compton drag on large scales are non-negligible. In this case, since the Compton drag force is proportional to the total-matter displacement, the two-point functions in Eq. 4.5 involving $\mathbf{\Delta}^{r}$ will gain a contribution proportional to $\mathbf{\Delta}^{m}$ (Fig. 5). Such contributions can be non-negligible at large scales and can dominate in the contributions to the power spectrum proportional to $f_{g}$ at low wavenumber (Fig. 6). Importantly, terms proportional to $b_{\pm}$ are unaffected since they are related only to the primordial modes $\mathbf{r}_{\pm}$, breaking the degeneracy between $f_{g}$ and $b_{+}$. Since the difference between these terms is proportional to the total-matter component, this difference can alternatively be absorbed into the total-matter bias $b_{m}$ [17]. Comparisons of these terms with and without Compton drag are shown in Figure 6. Comparing the $f_{g}$ contribution with and without Compton drag we see, as expected, that renormalizing the linear total-matter bias $b_{1}$ to include a contribution proportional to $f_{g}D_{\rm CD}(z)$ (purple dotted curve) is sufficient to account for the non-gravitational Compton drag contributions.

4.2 Cross-Spectra of different tracers and IR Resummation

So far we have dealt only with tracer auto-spectra. The situation for cross-spectra is complicated by the non-cancellation of the IR-exponent at small separations, $q$. For two generic fluids, $X$ and $Y$, such that $\boldsymbol{\Psi}_{X,Y}=\boldsymbol{\Psi}_{m}+f_{X,Y}\boldsymbol{\Psi}_{r}$, the cross spectrum will take the form as in Equation 2.2:

[TABLE]

where the exponentiated two-point function $A^{XY}$ is given by

[TABLE]

and expectation values of point operators are displayed without arguments. Both terms in parentheses on the RHS of Equation 4.10 are well-defined and invariant under generalized Galilean transformations; however the second term vanishes as $q\rightarrow 0$ while the first does not333A similar non-cancellation occurs in the modeling of BAO reconstruction, where the cross-term between the ‘displaced’ and ‘shifted’ fields exhibits the same behavior [33, 35].. As first noted in Ref. [15], this is in contrast to the single-fluid case where $A_{ij}$ had to vanish at small scales due to Galilean invariance.

In principle, the non-cancellation discussed above will introduce a large scale damping in the power spectrum at scales proportional to the difference $|\boldsymbol{\Psi}^{X}-\boldsymbol{\Psi}^{Y}|^{2}$. However, since $\boldsymbol{\Psi}^{X,Y}$ are both expected to have the same coefficient in the total-matter component (i.e. unity) this difference squared will generically be proportional to $(f_{X}-f_{Y})^{2}\mathcal{O}(\boldsymbol{\Psi}_{r}^{2})$, and thus is suppressed by about four orders of magnitude relative to the Zeldovich displacement, $\Sigma^{2}$, at the redshifts with which we are concerned ($z<10$). On the other hand, while differential streaming is expected to damp cross spectra negligibly even if $f_{X}$ is of order unity, as discussed in the previous section it will still generate an observable effect degenerate with the relative bias $b_{+}$.

4.3 Higher Order Bias

Thus far we have not discussed the fact that any perturbative model should be considered an effective field theory, working up to some scale $\Lambda$ [57, 58, 59]. This forces us to introduce a set of counterterms that remove the small scale sensitivities of the perturbative calculations. For instance all the $A_{ij}(q)$ terms contain a zero-lag piece computed at zero separation, i.e. $q=0$, where perturbation theory breaks down. In the single fluid case, this UV-sensitivity is renormalized to lowest order in the power spectrum by a counterterm $c_{s}k^{2}P_{ZA}(k)$ [29, 60], where the free parameter $c_{s}$ has to be matched to simulations or data. The same structure of the counterterms appears in the two fluid scenario: for instance, the $A^{ab}_{ij}(q)$ required to calculate auto and cross spectra feature the same UV-sensitive contributions as $q\rightarrow 0$, requiring one value of $c^{a}_{s}$ for each species. In principle, terms in the equations of motion due to the relative component will add additional UV sensitivities to our predictions; in practice, however, such contributions are subdominant in the dynamics of the relative component and negligible for the total-matter component (Appendix B). To the extent that these contributions can be ignored, then, the two-fluid equations of motion can be renormalized identically to the single fluid case with one set of counterterms for each species or tracer. As counterterms have minor impact on BAO scales, and are anyway fitted to the data in both the single and multiple fluid cases, we do not include them in the Fisher calculation in the next section.

We have equally refrained from discussing bias beyond linear order. As in the equations of motion, contributions beyond first order in the linear power spectrum proportional only to the total-matter component can be added consistently as in the single-fluid case, and we will ignore small nonlinear contributions proportional to one or more powers of the relative component444A proper accounting of such terms would in addition require solving the relative-component equations of motion to beyond linear order, which is beyond our present scope.. However, one exception must be made: operators involving the relative-velocity between the baryon and dark matter squared, which, despite being at second order in the relative component, can be non-negligible due to their distinct dimensional scaling [11, 18, 22, 16]. Such contributions were the focus of the first studies of bias [11, 18, 22] in the two-fluid picture, and we will show how their calculation fits naturally into the Lagrangian framework. For a discussion of other second order bias operators see Appendix B.2.

At second order in the bias expansion we can write

[TABLE]

where $\sigma_{v_{r}}^{2}$ is the 1-point variance of the relative velocities and $\sigma_{r_{-}}^{2}=(1+z)^{-2}H_{0}^{-2}\sigma_{v_{r}}^{2}$. As several authors [11, 16] have pointed out, baryon-dark matter relative velocities can be quite large at the time when the first halos and galaxies form, which could result in a large value of $b_{v}$ for their late time descendants. The value of $b_{v}\sigma_{v_{r}}^{2}$ can be as large as 0.01, which will make this contribution at second order in the power spectrum larger than the $b_{-}$ terms, even on linear scales. It is however worth remembering that a value of $b_{v}\sigma_{v_{r}}^{2}\simeq 10^{-5}$ is also plausible, which would substantially reduce the importance of this contribution.

To consistently compute the power spectrum contributions due to $c_{-}\sim b_{v^{2}}$ we must go beyond the Zeldovich approximation. Up to 1-loop in Lagrangian perturbation theory we have to compute 4 new terms to properly include the new bias parameter $c_{-}$. Beyond these, terms proportional to $c_{-}^{2}$ can be safely neglected as they are $\mathcal{O}(P_{\nabla\mathbf{r}_{-}^{2}}^{2})$. For the same reason we drop all the terms proportional to $b_{\pm}c_{-}$, as well as contributions of the relative component to the equations of motion. This leaves us with contributions proportional to $c_{-}$, $b_{1}c_{-}$, $b_{2}c_{-}$, and $b_{s^{2}}c_{-}$.

The first of these, proportional to $c_{-}$, contains a 1-loop contribution and is given by

[TABLE]

where we have defined $A^{m-}_{ij}=\left\langle\,\Delta_{m,i}\ (\mathbf{r}_{-,2}-\mathbf{r}_{-,1})_{j}\,\right\rangle$ and the 1-loop contribution from the second-order Lagrangian displacement $\boldsymbol{\Psi}^{(2)}$ enters as

[TABLE]

The kernel $Q_{v^{2}}$ is derived in Appendix C.

The remaining terms do not contain loop contributions and follow straightforwardly from evaluating the second and third cumulants in Eq. (4.1) within the Zeldovich approximation. These are those proportional to the first order bias:

[TABLE]

second order bias:

[TABLE]

and shear

[TABLE]

where we have defined the 2-point functions $U^{-m}\equiv\left\langle\,\mathbf{r}_{-}(q)\delta_{m}(0)\,\right\rangle=U^{m-}$ and $W^{s-}_{ijk}(q)=\left\langle\,s_{ij}(q)\mathbf{r}_{-,k}(0)\,\right\rangle$. Details of the above calculation can be found in Appendix C.

The contributions proportional to $c_{-}$ and their comparison with the 1-piece in Eq. (4.7) and with the $b_{\pm}$ ones computed in the previous section is shown in Fig 7 for $z=1.2$, assuming $b_{v}\sigma_{v_{r}^{2}}^{2}=0.01$. The $c_{-}$ terms are indeed larger than the $b_{-}$ terms on most scales, but still subdominant compared to the $b_{+}$ terms. Notably, the $c_{-}$ terms feature significantly larger oscillatory features than contributions from $b_{\pm}$, with minima that differ from maxima by more than an order of magnitude.

5 Degeneracies and bias to BAO

Baryon acoustic oscillations (BAO) in the photon-baryon fluid before combination imprint a characteristic clustering scale in the distribution of galaxies that can be used as a standard ruler to constrain the cosmic expansion history [61]. In general this method is regarded as highly robust as it probes very large scales which are largely unaffected by astrophysical processes. However, relative component contributions to the two-point function also occur on very large scales and their oscillatory features, although arising from the same physical process of the standard BAO features in the matter density power spectrum, could bias our estimates of the distance scale if not properly taken into account [22, 23, 24]. Indeed, as shown in the left panel of Fig. 8, all the relative component contributions we have considered show distinct features around the BAO peak.

The extent to which contributions from the relative component can contaminate measurements of the BAO scale can be estimated using the Fisher matrix formalism [62]. The galaxy overdensity has a covariance that is diagonal in Fourier space and given by the power spectrum plus shot noise, $\hat{P}_{gg}=P_{gg}(k)+\bar{n}^{-1}$; for the parameters $\{\theta_{i}\}$, the Fisher matrix is given by

[TABLE]

where $V_{obs}$ is the observed volume. For simplicity we neglect redshift space distortions and focus only on the isotropic BAO signal, though we will comment on how our Lagrangian analysis can be naturally extended to redshift space in the final paragraph. We model the power spectrum using the two-fluid Zeldovich terms derived above and include matter contributions up to one loop (see e.g. [52]), including contributions from the quadratic Lagrangian bias $b_{2}$. We consider only scales between $k_{\rm min}=10^{-2}\,\,h\,{\rm Mpc}^{-1}$ and $k_{\rm max}=0.25\,\,h\,{\rm Mpc}^{-1}$, and fiducial value of $b_{1}=0.53$ and $b_{2}=0.2$. The number density of galaxies is $\bar{n}=4.2\times 10^{-4}\ \,h^{3}\,{\rm Mpc}^{-3}$ and we assume $V=5\,h^{-3}\,{\rm Gpc}^{3}$. These numbers are chosen to be similar to what galaxy surveys like DESI [2] or Euclid [3] are expected to measure, and in particular are based off the expected DESI ELG population at $z=1.25$ in a bin of width $\Delta z=0.1$ and 14,000 square degrees of observation.

To quantify the potential impact of the relative component on standard BAO analyses, we will compare two models of the power spectrum within the Fisher formalism: the “correct” model $\textbf{M}_{1}$, which is a function of all total-matter and relative component biases, and the nested “standard” model $\textbf{M}_{0}$, wherein the relative component biases are set to zero (i.e. $b_{\pm},\ c_{-}=0$). The observed power spectrum is in addition a function of the BAO scaling parameter $\alpha$ such that

[TABLE]

The derivative of the baseline galaxy power spectrum with respect to the parameters is shown in Figure 8. These templates all show oscillatory features of roughly the same frequency as the BAO scale but exhibit distinguishable scale dependence. For reference, applying Eq. (5.1) returns sub-% error on the BAO scale, with $\sigma_{\alpha}=0.9\%$, for the standard analysis using $\textbf{M}_{0}$.

We can now compute the systematic shifts in $\alpha$ that would be incurred by neglecting the relative component, i.e. by fitting to $\textbf{M}_{0}$. For convenience, we will split the parameters in $\textbf{M}_{1}$ into $\theta=(\phi_{a},\psi_{\sigma}),$ where $\phi_{a}$ with Latin indices are the BAO scale and total-matter parameters and $\psi_{\sigma}$ with Greek indices are the relative component biases, such that $\textbf{M}_{0}$ is given by $\theta=(\phi_{a},\psi_{\sigma}=0).$ In this language the shift in $\alpha$ and $b_{1}$ due to using the standard model can be calculated to first order as [63]

[TABLE]

Here $F_{0}$ and $G$ are respectively diagonal and off-diagonal blocks of the full Fisher matrix $F=F(\theta_{0})$ calculated at the best fit parameters $\theta_{0}$ for the full model $\textbf{M}_{1}$, such that $F_{0,ab}=F_{ab}$ and $G_{b\sigma}=F_{b\sigma}$, and $\delta\psi$ is the deviation of $\psi$ in the standard cold dark matter only model $\textbf{M}_{0}$ from $\mathbf{M}_{1}$, i.e. $\delta\psi=-\psi_{0}.$

As a simple first example, we consider a toy-model Universe in which the only relative contribution is $b_{+}$. Figure 9 compares the “true” power spectrum, $P_{1}(k)$, assuming $b_{+}=5$, with best fits to the power spectrum in a dark matter only universe $P_{0}(k)$, described by the model parameters $\textbf{M}_{0}$, where the values of $\alpha$, $b_{1}$, $b_{2}$ are shifted from their true values according to Equation 5.3. Different values of the maximum wave number $k_{\rm max}$ included in the Fisher calculation are shown with different lines. For $k_{\rm max}=0.1\,h\,{\rm Mpc}^{-1}$, we find a significant departure in phase between the two models, compared to higher limiting wavenumbers, as evident from the phase of the residual in the bottom panel. Beyond $k_{\rm max}=0.15\,h\,{\rm Mpc}^{-1}$ there are sufficient BAO wiggles that the phase of the residuals are essentially locked. We caution that the same exercise repeated with both matter and relative terms in the Zeldovich approximation can lead to wide swings in the BAO scale $\delta\alpha$ as a function of $k_{\rm max}$. This can be understood as follows: at $k\gtrsim 0.1\,h\,{\rm Mpc}^{-1}$, $b_{+}$ contributes both oscillatory behavior and a broadband shape identical to the total matter component. The latter is essentially an amplitude change and can be roughly cancelled by a shift $\delta b_{m},$ which it is thus fixed independently of $k_{\rm max}$. This then requires $\delta\alpha$ to shift with $k_{\rm max}$ as more oscillations are included until the oscillations in $\mathbf{r}_{+}$ relative to $\mathbf{m}_{+}$ are damped at large $k$ (Figure 2). This broadband effect is ameliorated by including nonlinear terms for BAO measurements, but the partial degeneracy of $b_{+}$ with the power spectrum amplitude likely implies that ignoring two-fluid effects may affect measurement of the amplitude of the power spectrum (though this effect will also be partially mitigated by redshift-space distortions).

The same formalism can be applied to more realistic bias models. In the upper left panel of Figure 10 we consider the case when the observed power spectrum contains nonzero values $b_{\pm}$ and $c_{-}=0,$ and forecast the shifts in $\alpha$ due to the wrong assumption of $b_{\pm}=0$. Due to the small size of the $b_{-}$ contributions (see Figure 4), we expect shifts in BAO inferred distances to be dominated by $b_{+}$, and this is indeed what we find, contours of constant $\delta\alpha$ are almost independent of $b_{-}$ even when $|b_{-}|=10.$ On the other hand, we see that values of $b_{+}\sim 5$ shift the measured $\alpha$ by up to $0.4\%$, close to half of the error on $\alpha$ expected when using $\textbf{M}_{0}$.

However, the physics behind the relative components is quite well understood and can be easily included in Fisher forecasts or power spectrum analyses. Indeed, as seen in Figure 8, the templates for the various relative biases and $\alpha$ have distinct shape and could be distinguishable depending on the noise level of the measurements. The upper right plot in Figure 10 shows the increase in $\sigma_{\alpha}$ induced by marginalizing over $b_{\pm}$ in universes where $b_{\pm}$ and $b_{v}$ are not necessarily nonzero 555The nonzero $b_{\pm,v}$ produce shifts in the measured $\alpha,b_{m}$ when using $\textbf{M}_{0}$, which must be taken account when computing $\sigma_{\alpha}.$ To first order, the shifted Fisher matrix is given by $F_{0}$.. The total loss of constraining power is modest, with less than 10% worse error bars even after marginalizing over two extra parameters. In both the computations of the shifts in $\alpha$ and the increase of $\sigma(\alpha)$, the volume of the survey does not enter, and the final results depend only on the shot noise levels.

In the lower set of plots in Figure 8, we repeat the same exercise described above including $c_{-}$ as an extra free parameter. Since $b_{-}$ is irrelevant for the final results we set it to zero (but still marginalized over it). We find that $b_{+}$ and $c_{-}$ are anti-correlated, with larger shifts compared to the $b_{\pm}$ case, but $\delta\alpha/\alpha\leq 0.5\%$ in all cases. Marginalizing over the extra parameter $c_{-}$ results in a 20-30% increase in $\sigma(\alpha)$, which is still benign for BAO constraints. Our results therefore advocate for the implementation of relative component biases, at least of $b_{+}$ and $c_{-}$, in standard BAO data analysis of the galaxy power spectrum or correlation function.

Finally, in Figure 11 we investigate the detectability of the two-fluid effects in the same setup. On their own, both $b_{+}$ and $c_{-}$ become 1$\sigma$ detectable at the upper end of our explored parameter ranges, shown as the red lines in Figure 11. However, once all three relative bias parameters are marginalized over, the black set of curves in Figure 11, neither will be detectable within our fiducial volumes, with $c_{-}$ in particular at $0.1\sigma$, well out of reach even if all the DESI redshift bins are combined.

6 Conclusions

The large scale structure of the universe, whose formation is dominated by the dynamics of gravitational collapse, is one of the premier probes into fundamental physics. At subleading order, the presence of multiple particle species, broadly categorized into cold dark matter, baryons and neutrinos, with distinct properties beyond their shared gravitational attraction, can present additional features in this structure, which will become increasingly important as future surveys push to higher precision. In particular, relative perturbations between baryons and cold dark matter are prominent at the same scale as baryon acoustic oscillations and have the potential to cause systematic biases in future BAO measurements.

In this paper, we develop the Lagrangian formalism to calculate the clustering of biased tracers in the presence of multiple fluids, focusing specifically on the two-fluid scenario with dark matter and baryons. The Eulerian description of two-fluid dynamics has been studied extensively in the past and we make contact with previous work as appropriate throughout the text. LPT includes an automatic resummation over long-wavelength bulk flows and is thus able to accurately capture the shape of BAO features for biased tracers. In addition, LPT naturally maps bias terms from their initial Lagrangian positions to advected Eulerian positions, in contrast to Eulerian theory in which advective terms must be put in by hand, thereby simplifying the treatment of bias as responses to linear initial perturbations.

The presence of two fluids introduces terms beyond those encountered in traditional single fluid cosmological perturbation theory, with modifications in both the bias expansion and tracer advection. In the former, the generalized Galilean invariance that restricted the bias to contain only second derivatives of the gravitational potential in the single fluid case, allows terms including relative overdensities and velocities between different species. In the latter, initial relative displacements between various species are preserved under free fall and present an additional source of bias. Large scale non-gravitational forces such as Compton drag induced by the CMB can introduce additional corrections. We formulate modifications to tracer bias and advection in terms of three initial modes, constants of motion in the linear equations of motion, which roughly correspond to the initial total-matter displacement field and the relative displacement and velocity fields between dark matter and baryons.

We explicitly calculate the galaxy auto-power spectrum in the Zeldovich approximation within this formalism. Cross correlations between the relative modes introduce eight terms linear in the power spectrum—however, those quadratic in the relative component are suppressed by four orders of magnitude relative to the single fluid terms at low redshifts relevant for the next generation of galaxy surveys. Comparing to the Eulerian result explicitly to first order in the power spectrum, we find that the Eulerian relative component bias corresponds to linear mixtures of the Lagrangian bias, with modifications to the tracer advection entering both the Eulerian relative overdensity bias and the Eulerian relative velocity divergence bias. We then take up the calculation of cross spectra, finding a large scale damping due to an IR noncancellation in the relative component that is nonetheless negligibly small on perturbative scales. We also briefly discuss higher order corrections to the equations of motion in the presence of two fluids from an effective theory point of view, and perform an example one loop calculation for the relative velocity effect ($\propto\mathbf{v}_{r}^{2}$).

We conduct an exploratory analysis into whether two-fluid effects can cause systematic biases in measurements of the BAO scale. Taking the example of DESI ELGs at $z=1.25,$ we show that while ignoring two-fluid effects can lead to systematic shifts in the measured BAO scale as large as half a sigma, properly marginalizing over these effects induces less than ten percent loss in precision for a wide range of bias values. Since the scale dependence of the underlying physics is well understood, these results advocate for including two-fluid terms at linear order in future analyses. The dominant relative bias term ($\propto b_{+}$) does not fall quadratically with the growth factor like the total-matter contributions, and we therefore expect the relative bias signal as a fraction of total power to scale with redshift as $D^{-1}_{+}(z)$ and become proportionally more significant for surveys (such as the proposed Stage II 21-cm survey [5]) at higher redshifts. Studies of more highly biased tracers such as DESI quasars [2], for which the total-matter contributions are correspondingly larger, will on the other hand be less influenced by the relative bias for similar reasons.

While the Lagrangian picture is a natural playground for their study, in this paper we have opted not to study redshift space distortions (RSD). We note, however, that of the two relative components, $\mathbf{r}_{+}$ is dominant but stationary while $\mathbf{r}_{-}$ is so small as to be essentially negligible— two-fluid impacts should thus have a relatively small impact on RSD. However, as noted in the previous section, since the dominant relative component contribution $b_{+}$ is somewhat degenerate with the overall power spectrum amplitude, it is possible that two-fluid effects could hinder the accuracy of $f\sigma_{8}$ measurements beyond the percent level. We will return to this issue in future work.

Acknowledgments

We thank Jonathan Blazek and Jahmour Givans for helpful comments. SC is supported by the National Science Foundation Graduate Research Fellowship (Grant No. DGE 1106400) and by the UC Berkeley Theoretical Astrophysics Center Astronomy and Astrophysics Graduate Fellowship. M.W. is supported by the U.S. Department of Energy and by NSF grant number 1713791. This work made extensive use of the NASA Astrophysics Data System and of the astro-ph preprint archive at arXiv.org.

Appendix A Redshift dependence and size of the of bias parameters

The bias expansion in Equation 3.5 has an implicit dependence on the initial redshift $z_{i}$ that must be taken into account to reach consistent conclusions. Since the initial conditions mix at most linearly, no information can be lost by choosing one initial time $\tau_{i}$ over another; for example, the sensitivity of halos to the relative velocity divergence after reionization, which contains a contribution from the total matter overdensity (Eq. 2.9), can be directly accounted for by calibrating the bias parameter for $\delta_{m}$ at an earlier redshift.

As a simple example we consider the redshift dependence of the relative components in the sourceless ($F_{b}=0$) case. If we set our initial time at $\tau_{i}^{\prime}$ instead of $\tau_{i}$ we will get

[TABLE]

Re-expanding $F_{g}$ at $\tau_{i}^{\prime}$ thus yields

[TABLE]

Since $b^{\prime}$ and $b$ apply to the same field configurations at different times, they must yield the same initial overdensity $F_{g}$ — this requirement can be satisfied by enforcing the differential equations

[TABLE]

Intriguingly, the presence of a relative overdensity bias can “generate” a relative velocity bias at later times. This can be understood as follows: the relative overdensity at late times is a linear combination of the relative overdensity and velocities at earlier times. Similar, though more complicated, versions of this relation hold when $F_{b}\propto\mathbf{m}_{+}$, in which case mixing of all three initial fields must be taken into account.

Appendix B Beyond Linear Order

B.1 Equations of motions

In this appendix we derive the equations of motion beyond linear order in the two-fluid scenario, and show that the nonlinear contribution of $\mathbf{r}_{\pm}$ to the total-matter component are quadratically suppressed, and that the nonlinear relative component is always sourced by at least one component of $\mathbf{r}_{\pm}$.

The Lagrangian equations of motion at higher order can be found by taking the real-space divergence of both equations in 2.3. To do so we make use of the identities

[TABLE]

where the negative powers in the second line denote matrix inverses, to account for the coordinate transformations between Lagrangian coordinates $q$ and the fluid trajectories for each species $x^{a}=q+\boldsymbol{\Psi}^{a},$ as well as the standard matrix identities $(I+A)^{-1}=I-A+A^{2}-A^{3}+\mathcal{O}(A^{4})$ and $\det(I+A)=1+\text{Tr}[A]+\frac{1}{2}(\text{Tr}[A]^{2}-\text{Tr}[A^{2}])+\mathcal{O}(A^{3})$. We will neglect the effects of Compton drag, which affects the relative displacement at a few percent level even at late times and on linear scales and thus enter into our final power spectra at the same order of magnitude as the relative component squared, and assume potential flow.

The above equations imply that the Lagrangian equations of motion for the fluid displacements of a species $a$ at nth order takes the generic form

[TABLE]

where the superscript $(a,n)$ denotes the species and order of each term, the derivative operator $D$ is defined such that $DX=X^{\prime\prime}+\mathcal{H}X^{\prime}$, and the $F^{(a,n)}$’s are kernels composed of displacements of the species $a$ at order $n$ and below. Switching to the total matter and relative components, we have

[TABLE]

We can derive some elementary properties of these equations without solving for their particular forms using symmetry arguments. Noting that the RHS of B.4 is symmetric under species index exchange ($b\leftrightarrow c$) while B.5 is antisymmetric, we can conclude that (1) the first relative contribution to the total matter EOM at each order must be of order $\mathcal{O}(\mathbf{r}_{\pm}^{2})$, and all subsequent contributions suppressed by further even powers of the linear relative component and (2) the relative component is always sourced by at least one power of $\mathbf{r}_{\pm}$, since the total matter component is even under this swap while the relative component is even. Note that (1) implies that the dynamics of the total-matter displacement are affected by the relative component only at the percent-of-a-percent level, and (2) implies that the nonlinear relative displacement is never less suppressed in $\mathbf{r}_{\pm}$ than the linear solution.. A similar result occurs in Eulerian theory, as described in [15].

For completeness, the explicit second and third order equations of motion for the relative component are, up to first order in the linear relative perturbation

[TABLE]

which are in-line with the symmetry arguments outlined above. The equations of motion for the total matter component at second and third order can be similarly verified to be simply the equations of motion in the one-fluid case to this order in the relative component.

B.2 Biasing at second order

Below we list all contributions to the bias expansion up to second order in the initial fields omitting derivative corrections:

[TABLE]

In the main body of this paper we consider relative bias terms up to first order in the power spectrum, since even these represent only percent level effects, with the exception of the relative velocity effect $\propto\mathbf{v}_{r}^{2}$, which has a distinct scaling. As noted in the text, we note that the presence of Compton drag can introduce additional terms due to loss of gauge redundancy; we refer readers to the extensive discussion in [17].

Appendix C Relative Velocity Bias Terms

In this appendix we provide details for the contributions of the relative velocity bias $b_{v^{2}}$ at $\mathcal{O}(P^{2})$ to the galaxy power spectrum. These contributions require the calculation of two new 2-point functions, the one-loop correlation between matter displacements and the squared relative velocity, and the correlation function between the shear field $s_{ij}$ and the relative velocity. We describe these in turn.

The second order solution to the total-matter displacement (correct up to first order in the relative component) is given by

[TABLE]

and more simply the “normalized” relative velocity at first order is given by

[TABLE]

From this we can calculate the two-point function

[TABLE]

which can be simplified to give

[TABLE]

where the kernel is defined as

[TABLE]

Next, the shear-velocity correlation function $W^{s-}_{ijk}$ is given in Fourier space by

[TABLE]

where in the last equality we have split $W^{s-}_{ijk}$ into a totally-symmetric piece and a familiar piece proportional to $U^{m-}$. The former can be decomposed into scalar components

[TABLE]

with the scalar components defined as spherical Bessel transformations:

[TABLE]

[TABLE]