Bias Loop Corrections to the Galaxy Bispectrum

TL;DR

This paper develops a comprehensive fourth-order galaxy bias model using Galilean invariants, enabling consistent inclusion of loop corrections in bispectrum analyses to improve cosmological parameter constraints.

Contribution

It introduces a fourth-order bias expansion with a Galilean invariant basis, facilitating loop correction calculations and bias evolution modeling for galaxy bispectrum analysis.

Findings

Bias parameters depend only on nonlocal operators in the split-basis.

Galilean invariance simplifies the computation of bias evolution.

Baseline model enables consistent loop correction inclusion in large-scale structure data.

Abstract

Combination of the power spectrum and bispectrum is a powerful way of breaking degeneracies between galaxy bias and cosmological parameters, enabling us to maximize the constraining power from galaxy surveys. Recent cosmological constraints treat the power spectrum and bispectrum on an uneven footing: they include one-loop bias corrections for the power spectrum but not the bispectrum. To bridge this gap, we develop the galaxy bias description up to fourth order in perturbation theory, conveniently expressed through a basis of Galilean invariants that clearly split contributions that are local and nonlocal in the second derivatives of the linear gravitational potential. In addition, we consider relevant contributions from short-range nonlocality (higher-derivative terms), stress-tensor corrections and stochasticity. To sidestep the usual renormalization of bias parameters that complicates predictions beyond leading order, we recast the bias expansion in terms of multipoint propagators, which take a simple form in our split-basis with loop corrections depending only on bias parameters corresponding to nonlocal operators. We show how to take advantage of Galilean invariance to compute the time evolution of bias and present results for the fourth-order parameters for the first time. We also discuss the possibilities of reducing the bias parameter space by using the evolution of bias and exploiting degeneracies between various bias contributions in the large-scale limit. Our baseline model allows to verify these simplifications for any application to large-scale structure data sets.

Figures

Tables

Table 1. Table 1: Overview of basis operators for galaxy bias, along with their associated bias parameters. Each single column presents a different order in SPT and we have categorized operators into groups, which are local and nonlocal in second derivatives of the linear gravitational potential (left and middle columns, respectively), and which contain higher than second derivatives of the potentials. Note that to simplify notation, we have relabeled some bias parameters compared to [ 9 ] — what was formerly γ 3 × superscript subscript 𝛾 3 \gamma_{3}^{\times} is now γ 2 × superscript subscript 𝛾 2 \gamma_{2}^{\times} and γ 3 − superscript subscript 𝛾 3 \gamma_{3}^{-} has become − 7 / 4 ​ γ 21 7 4 subscript 𝛾 21 -7/4\,\gamma_{21} . A superscript × \times denotes that a Galileon field has been multiplied by δ 𝛿 \delta . Operators related to noise terms are not listed here, see section  V.3 .

Local Evolution Operators				Nonlocal Evolution Operators		Higher-derivative
$1$st	$2$nd	$3$rd	$4$th	$3$rd	$4$th	$3$rd	$4$th
$b_1\delta$	$b_2{\delta }^2/2$	$b_3{\delta }^3/3!$	$b_4{\delta }^4/4!$	${\gamma }_{21}{𝒢}_2({\phi }_2,{\phi }_1)$	${\gamma }_{21}^{\times }\delta {𝒢}_2({\phi }_2,{\phi }_1)$	${\beta }_1{\nabla }^2\delta$	${\beta }_{2,1}{\nabla }^2{\delta }^2$
	${\gamma }_2{𝒢}_2\left({\Phi }_v\right)$	${\gamma }_2^{\times }\delta {𝒢}_2\left({\Phi }_v\right)$	${\gamma }_2^{\times \times }{\delta }^2{𝒢}_2\left({\Phi }_v\right)$		${\gamma }_{211}{𝒢}_3({\phi }_2,{\phi }_1,{\phi }_1)$		${\beta }_{2,2}{\left(\nabla \delta \right)}^2$
		${\gamma }_3{𝒢}_3\left({\Phi }_v\right)$	${\gamma }_3^{\times }\delta {𝒢}_3\left({\Phi }_v\right)$		${\gamma }_{22}{𝒢}_2({\phi }_2,{\phi }_2)$		${\beta }_{2,3}{\nabla }^2{𝒢}_2\left({\Phi }_v\right)$
			${\gamma }_2^{\mathrm{sq}}{𝒢}_2{\left({\Phi }_v\right)}^2$		${\gamma }_{31}{𝒢}_2({\phi }_3,{\phi }_1)$		${\beta }_{2,4}{𝒢}_2\left(\nabla {\Phi }_v\right)$

Table 2. Table 2: Number of bias parameters required for the galaxy power spectrum and bispectrum up to one-loop order in perturbation theory, organized according to the various contributions to the bias expansion. Of the 11 general bias parameters for P g + B g subscript 𝑃 𝑔 subscript 𝐵 𝑔 P_{g}+B_{g} description, 4 are enough to describe P g subscript 𝑃 𝑔 P_{g} alone ( b 1 , b 2 , γ 2 , γ 21 subscript 𝑏 1 subscript 𝑏 2 subscript 𝛾 2 subscript 𝛾 21 b_{1},b_{2},\gamma_{2},\gamma_{21} ), while adding B g subscript 𝐵 𝑔 B_{g} requires 7 more parameters ( γ 2 × , γ 21 × , b 3 , γ 3 , γ 22 , γ 31 , γ 211 superscript subscript 𝛾 2 superscript subscript 𝛾 21 subscript 𝑏 3 subscript 𝛾 3 subscript 𝛾 22 subscript 𝛾 31 subscript 𝛾 211 \gamma_{2}^{\times},\gamma_{21}^{\times},b_{3},\gamma_{3},\gamma_{22},\gamma_{31},\gamma_{211} ).

	General bias expansion	Stochastic		Stress-tensor & Higher-derivative
		$𝒪\left(k^0\right)$	$𝒪\left(k^2\right)$
$P_g$	11	1	1	1
$B_g$		2	3	5

Equations

$${\delta }_g=\sum _n\frac{1}{n!}{\left(\frac{{\partial }^n{\delta }_g}{\partial {\delta }^n}\right)}_0{\delta }^n\equiv$$

$$+\dots ,$$

$$⟨{\delta }_g({\mathbf{k}}_1){\delta }_g({\mathbf{k}}_2)⟩$$

$$⟨{\delta }_g({\mathbf{k}}_1){\delta }_g({\mathbf{k}}_2){\delta }_g({\mathbf{k}}_3)⟩$$

$${\delta }_g(\mathbf{x})={\int }_{\mathbf{k}}\exp (-i\mathbf{k}\cdot \mathbf{x}){\delta }_g(\mathbf{k}),$$

$$\begin{array}{rl}{P}_g(k)=& {\bar{b}}_1^2{P}_L(k)+\left[{\bar{b}}_1{\bar{b}}_3{\sigma }^2\right]P_L(k)\\ & +\frac{{\bar{b}}_2^2}{2}{\int }_{\mathbf{q}}{P}_L(|\mathbf{k}-\mathbf{q}|)P_L(q)+\dots ,\end{array}$$

$$P_g(k)=b_1^2{P}_L(k)+\frac{{\bar{b}}_2^2}{2}{\int }_{\mathbf{q}}{P}_L(|\mathbf{k}-\mathbf{q}|)P_L(q).$$

$$B_{g,123}=$$

$$+\frac{1}{2}{\bar{b}}_1{\bar{b}}_2{\bar{b}}_3{P}_1{\int }_{\mathbf{q}}{P}_L(|{\mathbf{k}}_2-\mathbf{q}|)P_L(q)$$

$$+{\bar{b}}_2^3{\int }_{\mathbf{q}}{P}_L(|{\mathbf{k}}_1-\mathbf{q}|)P_L(|{\mathbf{k}}_2+\mathbf{q}|)P_L(q)+cyc.$$

$$+\dots ,$$

$$B_{g,123}=$$

$$+{\bar{b}}_2^3{\int }_{\mathbf{q}}{P}_L(|{\mathbf{k}}_1-\mathbf{q}|)P_L(|{\mathbf{k}}_2+\mathbf{q}|)P_L(q)+cyc.$$

$$b_1=⟨\frac{\partial {\delta }_g}{\partial \delta }⟩={\bar{b}}_1+{\bar{b}}_3\frac{{\sigma }^2}{2}+\dots =\sum _{n=0}^{\infty }\frac{{\bar{b}}_{2n+1}}{n!}{\left(\frac{{\sigma }^2}{2}\right)}^n,$$

$$b_2=⟨\frac{{\partial }^2{\delta }_g}{\partial {\delta }^2}⟩={\bar{b}}_2+{\bar{b}}_4\frac{{\sigma }^2}{2}+\dots =\sum _{n=0}^{\infty }\frac{{\bar{b}}_{2n+2}}{n!}{\left(\frac{{\sigma }^2}{2}\right)}^n,$$

$${\left(\frac{{\partial }^n{\delta }_g}{\partial {\delta }^n}\right)}_0⟶⟨\frac{{\partial }^n{\delta }_g}{\partial {\delta }^n}⟩,$$

$${\delta }_g=$$

$$\displaystyle+{1\over 3!}\left<\frac{\partial^{3}\delta_{g}}{\partial\,\delta^{3}}\right>\ \Big{[}\delta^{3}-3\langle\delta^{2}\rangle\delta-\langle\delta^{3}\rangle\Big{]}+\ldots$$

$$=$$

$$\displaystyle+{b_{4}\over 4!}\Big{[}\delta^{4}-4\langle\delta^{3}\rangle\delta-6\langle\delta^{2}\rangle\delta^{2}+3\langle\delta^{2}\rangle^{2}-\langle\delta^{4}\rangle_{c}\Big{]}$$

$$+\dots$$

$$⟨{\delta }_g\delta ⟩$$

$$⟨{\delta }_g\delta \delta ⟩$$

$$P_g(k)=(2\pi {)}^3\sum _n\frac{b_n^2}{n!}\int {\text{displaylimits}}_{{\mathbf{q}}_1,\dots ,{\mathbf{q}}_n}[{\delta }_D{]}_n{P}_L(q_1)\dots P_L(q_n),$$

$${\xi }_g(r)=\sum _n\frac{b_n^2}{n!}{\left[\xi (r)\right]}^n,$$

$$\delta (\mathbf{x},\tau )=\sum _{n=1}^{\infty }{D}^n(\tau ){\delta }^{(n)}(\mathbf{x}),$$

$${\delta }^{(n)}(\mathbf{k})=(2\pi {)}^3\int {\text{displaylimits}}_{{\mathbf{k}}_1,\dots ,{\mathbf{k}}_n}$$

$$\times {\delta }^{(1)}({\mathbf{k}}_1)\cdots {\delta }^{(1)}({\mathbf{k}}_n).$$

$$⟨\frac{{\partial }^n{\delta }_g}{\partial {\delta }_L^n}⟩⟶⟨\frac{{\partial }^n{\delta }_g(\mathbf{k})}{\partial {\delta }_L({\mathbf{k}}_1)\cdots \partial {\delta }_L({\mathbf{k}}_n)}⟩,$$

$$\frac{1}{2}{⟨\frac{\partial {\delta }^2(\mathbf{k})}{\partial {\delta }_L({\mathbf{k}}_1)}⟩}^{\prime }$$

$$=2\underset{{\bar{F}}_2{\sigma }^2}{\underbrace{{\int }_{\mathbf{q}}{F}_2({\mathbf{k}}_1,\mathbf{q})P_L(q)}}+\dots =\frac{34}{21}{\sigma }^2+\dots ,$$

$$\frac{1}{3!}{⟨\frac{{\partial }^2{\delta }^3}{\partial {\delta }_1^L\partial {\delta }_2^L}⟩}^{\prime }$$

$$\displaystyle=\Big{[}4\bar{F_{2}}+F_{2}(\bm{k}_{1},\bm{k}_{2})\Big{]}\,\sigma^{2}+\ldots$$

$${\bar{b}}_1{⟨\frac{{\partial }^2\delta }{\partial {\delta }_1^L\partial {\delta }_2^L}⟩}^{\prime }={\bar{b}}_{1}2F_2({\mathbf{k}}_1,{\mathbf{k}}_2)+\dots ,$$

$$\frac{1}{2}{⟨\frac{{\partial }^2{\delta }^2}{\partial {\delta }_1^L\partial {\delta }_2^L}⟩}^{\prime }$$

$$=$$

$$+{⟨\frac{\partial {\delta }^{(3)}}{\partial {\delta }_1^L}+\frac{\partial {\delta }^{(3)}}{\partial {\delta }_2^L}⟩}^{\prime }+\dots$$

$$\begin{split}\left<{\partial\delta^{(2)}\over\partial\delta^{L}_{1}}\ast{\partial\delta^{(2)}\over\partial\delta^{L}_{2}}\right>^{\prime}+\left<\delta_{L}\ast{\partial^{2}\delta^{(3)}\over\partial\delta^{L}_{1}\partial\delta_{2}^{L}}\right>^{\prime}&=\int_{\bm{q}}\Big{[}4F_{2}(-\bm{q},\bm{k}_{1})F_{2}(\bm{q},\bm{k}_{2})+6F_{3}(\bm{q},\bm{k}_{1},\bm{k}_{2})\Big{]}P_{L}(q)\\ &={68\over 21}F_{2}(\bm{k}_{1},\bm{k}_{2})\sigma^{2}+{8126\over 2205}\sigma^{2}+{254\over 2205}K(\bm{k}_{1},\bm{k}_{2})\sigma^{2}+\int_{\bm{q}}{\cal F}(\bm{k}_{1},\bm{k}_{2},q)P_{L}(q)\end{split}$$

Full text

I Introduction

Any mismodeling of galaxy bias — the relation between galaxies (or any other luminous tracer) and the underlying matter distribution [1, 2, 3, 4, 5, 6] — risks an incorrect recovery of cosmological parameters from large-scale structure (LSS) surveys. Fortunately, although the formation of galaxies involves highly nonlinear, small-scale processes, recent developments [7, 8, 9, 10, 11, 12, 13, 14] have shown that a perturbative expansion provides a robust treatment of galaxy bias on sufficiently large scales. This comes at the price of a set of unknown bias parameters, which, once marginalized over, degrade the statistical power for constraining cosmological models. It is therefore imperative to combine traditional LSS two-point statistics with higher-order statistics, such as the bispectrum, which allows us to break the degeneracies that exist between cosmological parameters and galaxy bias [15, 16, 17, 18, 19, 20, 21]. For this program to succeed we require consistent models of the galaxy power spectrum and bispectrum at leading (“tree-level”) and next-to-leading (“one-loop”) order, which should improve the regime of validity and robustness of the results [22, 23, 24, 9]. This is particularly important since at present the analysis of the bispectrum in galaxy surveys [25, 26] is done in a manner that is inconsistent with the treatment of the power spectrum, in that bias is treated at one-loop order for the power spectrum but tree-level for the bispectrum. The main goal of this study is to present such a unified model that includes all relevant effects in real space, extending the results of [12].

What are the essential elements? First, we need to consider all contributions from the general bias expansion up to fourth order in perturbation theory [6, 13, 12]. These are generated by the gravitational evolution of the dark matter field and include the common linear, and nonlinear bias parameters [1, 27, 28], as well as the second order nonlocal (or tidal field) bias [29, 8, 9, 10]. Second, it has been argued in [30, 31] that so-called higher-derivative terms, which have been known for a long time [3] and can be understood as deriving from a dependence of galaxy formation on the spatial distribution of matter [8], can be as important or even dominate over the general bias terms in certain circumstances. Third, on scales of the weakly nonlinear regime stress-tensor corrections to the evolution of dark matter [32, 33, 34, 35, 36, 37], might potentially become relevant. And finally, the impact of very small scales, in the absence of primordial non-Gaussianity largely uncorrelated with the previously mentioned large-scale effects, leads to an additional stochastic bias [38, 39].

However, a challenge that afflicts all theoretical predictions of galaxy clustering beyond leading order is a mismatch between the bias parameters from the perturbative expansion, and those an observer would define through the measurement of correlation functions. This was first pointed out in [7], which showed that appropriate renormalizations of the perturbative bias parameters restore agreement with the measurements. These renormalizations must be done on a statistic-by-statistic and parameter-by-parameter basis, and while tractable for the power spectrum [7, 8], this becomes increasingly complicated for higher-order statistics [12]. As a way of circumventing this procedure altogether, we advocate the use of multipoint propagators.

The multipoint propagator formalism was originally introduced in the context of renormalized perturbation theory (RPT) and its generalizations [40, 41, 42, 43, 44, 45] to obtain an accurate description of the dark matter field in the quasi-linear regime. Multipoint propagators have also been extended to include redshift space distortions and galaxy bias in [46, 47, 48], but without the clear connection to the renormalization issue that we aim to establish in this paper. In fact, we will demonstrate that the multipoint propagators correspond to the bias parameters that are commonly identified through cross-correlations of galaxy and matter fields [49, 22]. Thus, they are observable and have a well-defined physical meaning, based on which we argue that they provide the most natural approach towards galaxy bias. As in RPT they further act as building blocks for the general $N$-point correlation function, and determination of the first three propagators already fixes the power spectrum and bispectrum at the one-loop level, simplifying the computations significantly.

Modeling higher-order statistics to the same precision as their lower-order counterparts brings about a quickly inflating number of terms from the various bias contributions. In order to reduce the growing parameter space, a further goal of this article is to examine 1) degeneracies between the bias terms, and 2) relations between bias parameters arising from evolution, which are studied here for the first time at fourth order.

The paper is organized as follows. In Sec. II we present the central idea of the multipoint propagator expansion based on the simple and widely known model of local galaxy bias. A complete basis for the general bias and higher-derivative terms is provided in Sec. III. Following [9] this basis makes intuitive use of Galilean invariants of Lagrangian perturbation theory (LPT) potentials, and is particularly well suited for the computation of the multipoint propagators that we perform in Sec. IV. We also determine their time evolution and present the relations between initial and evolved bias parameters. Finally, Sec. V uses these results to compute the power spectrum and bispectrum, and demonstrates how we incorporate stress-tensor corrections and stochasticity. Sec. VI summarizes the final model components and gives our conclusions. The appendices A and B give further details and include relations between our basis for galaxy bias and others in the literature. Appendix C shows how to take advantage of Galilean invariance for determining the evolution of galaxy bias, while Appendix D demonstrates that the corresponding results are unaffected by renormalization.

II Basics of the renormalized bias expansion

We are interested in the statistical properties of the observed galaxy distribution. These are commonly quantified by a hierarchy of correlation functions of the density perturbations, which we write as $\delta_{g}(\bm{x})=n_{g}(\bm{x})/\bar{n}_{g}-1$, with $n_{g}(\bm{x})$ denoting the galaxy number density and $\bar{n}_{g}$ its average.111Note that while the subscript ‘g’ here stands for galaxies, it could equally well denote any other tracer of the mass field, such as quasars, galaxy clusters or the Ly-$\alpha$ forest. Model predictions of the correlation functions require a relation between the galaxy perturbations and matter fluctuations $\delta(\bm{x})$, which is usually written as some functional $\delta_{g}[\delta(\bm{x})]$ that is then Taylor expanded.

In order to illustrate the main idea pursued in this paper, we start with the simplest and most well-known bias relation, which considers $\delta_{g}$ to be a local function of $\delta$ [27, 28] and can be written as a Taylor series around $\delta=0$. Dropping all position arguments we have

[TABLE]

where the bias parameters are identified as $\bar{b}_{n}=(\partial^{n}\delta_{g}/\partial\,\delta^{n})_{0}$ ($n>0$) with $()_{0}$ denoting evaluation at $\delta=0$, and $\bar{b}_{0}$ enforces $\langle\delta_{g}\rangle=0$. However, the bias parameters so defined are not observables, as we measure correlators (expectation values or ensemble averages of fields), not quantities evaluated at $\delta=0$. As we will detail below, this is the reason for difficulties in the computation of galaxy correlation functions that ultimately require us to redefine (or “renormalize” [7]) the bias parameters above. We will further show that a more natural bias expansion, in the sense that its coefficients are actually measurable, can be defined in terms of cross-correlations between the galaxy and matter fields. In the language of renormalized perturbation theory (RPT) and its generalizations [40, 41] these cross-correlations correspond to the so-called multipoint propagators.

II.1 Galaxy clustering statistics in the

standard approach

We define the two- and three-point correlation functions of the galaxy perturbations in Fourier space — the power spectrum and bispectrum — as follows

[TABLE]

where $\bm{k}_{1\ldots n}\equiv\bm{k}_{1}+\ldots+\bm{k}_{n}$ and the appearance of the delta distribution is a manifestation of statistical homogeneity. Statistical isotropy further demands that the power spectrum only depends on the magnitude of the two wave vectors participating in the correlator, while the bispectrum is a function of three wave numbers $k_{1}$, $k_{2}$ and $k_{3}$. Finally, to pass from configuration to Fourier space we have adopted the convention

[TABLE]

using a short-hand notation for $k$-space integrals, i.e. $\int_{\bm{k}_{1},\ldots,\bm{k}_{n}}\equiv\int\text{d}^{3}k_{1}/(2\pi)^{3}\cdots\text{d}^{3}k_{n}/(2\pi)^{3}$.

Analogous definitions hold for the dark matter field, and to begin with let us assume that the bias expansion is done in Lagrangian space, so that the matter fluctuations can be considered linear and Gaussian (we neglect any possible primordial non-Gaussianities for the remainder of this paper). In that case all statistical information is contained in the linear power spectrum, which we denote by $\langle\delta_{L}(\bm{k}_{1})\,\delta_{L}(\bm{k}_{2})\rangle=(2\pi)^{3}\,P_{L}(k_{1})\,\delta_{D}(\bm{k}_{12})$, while the bispectrum and all other higher-order statistics vanish. Under this assumption, plugging Eq. (II) into Eq. (2) and using Wick’s theorem leads to

[TABLE]

where the dots denote contributions of two-loop and higher order, i.e. terms that involve five or more powers of $\delta$. The variance $\sigma^{2}\equiv\langle\delta^{2}(\bm{x})\rangle$, which appears in the first term of the loop contribution in square brackets, is formally infinite or at least highly sensitive to the nonlinear regime, depending on the shape of the linear power spectrum222Imposing a high-$k$ cutoff for the power spectrum would only mask the problem, as this would lead to all observables being dependent on an arbitrary choice of scale.. This implies that the large-scale galaxy power spectrum would be heavily influenced by scales where our perturbative approach is not expected to hold. However, we notice that this term is also proportional to $P_{L}(k)$, such that if we redefine the linear bias parameter to be $b_{1}\equiv\bar{b}_{1}+\bar{b}_{3}\,\sigma^{2}/2$, we retain the form

[TABLE]

Clearly, as more bias loops are included, the expression for $b_{1}$ keeps changing, but the principle remains the same — the observed linear bias is defined as (the square root of) the coefficient in front of $P_{L}$.

A similar situation holds for the quadratic bias parameter $b_{2}$. To see this let us now consider the bispectrum up to one-loop order (up to six powers of $\delta$), under the same assumption as above. We then get from Eq. (3):

[TABLE]

where $P_{i}\equiv P_{L}(k_{i})$ and cyc. denotes cyclic permutations of each term over the three wave vectors. Again, the first term in square brackets corresponds to the renormalization of the linear bias seen in the power spectrum, but in addition now there is also a second term that can be absorbed by a renormalization of the quadratic bias $b_{2}\equiv\bar{b}_{2}+\bar{b}_{4}\,\sigma^{2}/2$, so that we can write

[TABLE]

A number of questions arise from this procedure. 1) Is the $b_{2}$ renormalization that follows from the one-loop bispectrum consistent with the one that follows from the two-loop power spectrum? If so, is that true for all other bias parameters? 2) Is there a way to do all these renormalizations “automatically”, instead of calculating statistic by statistic, and order by order? 3) Is there a simple connection between renormalizations of different $N$-point correlators? As we will see next, the answer to all these questions is ‘yes’ [40, 50, 41, 43].

II.2 A bias expansion based on observables

First, we notice that the terms that contribute to the renormalization of bias parameters must derive from factorizable loop corrections, i.e. they can be written as products of lower-order terms and new contributions (in the example of Sec. II.1 these new contributions are functions of bias parameters and $\sigma^{2}$). In general, any loop correction can be portrayed as a diagram with a fixed number of external lines and a variable number of internal lines depending on the loop order. Those diagrams that lead to factorizable contributions are so-called reducible, meaning they can be decomposed into two or more connected diagrams by cutting one or more internal lines, while all other ones that do not share this property are classified as irreducible.

In order to automatically include renormalizations to arbitrary loop order, it is therefore desirable to construct the bias expansion in terms of the sum over all reducible diagrams with a given number of external lines. Such objects are known as multipoint propagators (see Fig. 1), and correspond to the expectation value of functional derivatives of the galaxy field with respect to the matter field: the derivatives produce the external lines, while the expectation value generates the loops [40].

An arbitrarily complicated loop diagram for any statistic can be decomposed easily and uniquely into multipoint propagators: the reducible subdiagrams are absorbed into the multipoint propagators while the irreducible subdiagrams are generated by connecting multipoint propagators among themselves. Since a given multipoint propagator appears in various statistics, this gives us a connection between renormalizations of different $N$-point correlation functions; e.g. the renormalization of $b_{2}$ that follows from the one-loop bispectrum is consistent with the one that follows from the two-loop power spectrum, because both are incorporated into the two-point propagator (i.e. they correspond to the same subdiagram in the power spectrum and bispectrum). As a result of this, $\bar{b}_{2}$ in Eq. (6) can be replaced by $b_{2}$ and similarly $\bar{b}_{1}\bar{b}_{2}\bar{b}_{3}\to b_{1}b_{2}b_{3}$ and $\bar{b}_{2}^{3}\to b_{2}^{3}$ in Eq. (8).

The observed bias parameters of order $n$ that appear in correlators thus correspond to the sum over all reducible diagrams with $n$ external legs. For example when $n=1$, i.e. linear bias, we simply need to consider the sum of all reducible diagrams with a single external leg (corresponding to $P_{L}$ at $k$):

[TABLE]

where we use the symbol $\partial$ to denote a functional derivative. Similarly, for quadratic bias ($n=2$) we have from Eq. (II)

[TABLE]

that is, the sum over all reducible diagrams with two external legs (corresponding to $P_{L}$ at $k_{1}$ and $k_{2}$ as seen in the leading order bispectrum, Eq. 8). Clearly, the calculations in Eqs. (9-10) are significantly easier than performing the renormalization procedure order by order and statistic by statistic leading to Eqs. (6, 8).

To conclude, we can remove the disconnect between the parameters appearing in the standard bias expansion (Eq. II) and those in correlators (Eqs. 6, 8), and therefore stop thinking about renormalization altogether, if we construct the bias expansion in terms of reducible diagrams with a given number of external legs. In effect, we trade kernels for multipoint propagators,

[TABLE]

which leads to a new expansion of the form (denoted as Gamma expansion in [41, 43])

[TABLE]

The structure of the square brackets is given by $\delta^{n}$ minus all possible actions of $\langle\rangle$ on $\delta^{n}$ with a constant term that respects that the expectation value is zero for non-Gaussian $\delta$, and we note that Eq. (13) automatically satisfies $\langle\delta_{g}\rangle=0$. The second equality assumes local galaxy bias as we have done so far, but we stress that in this expansion the bias parameters $b_{n}$ are the renormalized ones, that is, replacing the propagators by numbers gives precisely the renormalized local bias expansion.

Furthermore, if the expansion in Eq. (13) is done with respect to a Gaussian $\delta$, for instance when using the linear fluctuations, terms such as $\langle\delta^{3}\rangle$ and $\langle\delta^{4}\rangle_{c}$ will vanish and the expansion of $\delta_{g}$ will be given in terms of Hermite polynomials — a suggestion that was already put forward in [51]. Another crucial property of the multipoint propagators in this case is that they can be shown to be the cross-correlation bias between galaxies and matter fluctuations [50], e.g. for linear and quadratic bias we have

[TABLE]

That is, the observables corresponding to the multipoint propagators are no other than the standard cross-correlation Lagrangian bias coefficients routinely measured in N-body simulations, e.g.  [52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63].

Using the expansion in Eq. (13) allows us to rederive the renormalization procedure by writing it in the form of Eq. (II) and matching coefficients of $\delta^{n}$, e.g. $\bar{b}_{1}=b_{1}-b_{3}\sigma^{2}/2$ and $\bar{b}_{2}=b_{2}-b_{4}\sigma^{2}/2$, again, without doing any calculations of correlators. This can become useful if such relations are desired for the more general case of nonlocal bias, as we shall discuss below.

Finally, the multipoint propagator expansion has the additional advantage that it simplifies the computation of correlators considerably. For instance, for the power spectrum there is only one irreducible diagram present at each loop order, as opposed to an increasing number of diagrams at each loop order for the standard approach; indeed, assuming local bias for now we simply have:

[TABLE]

where $[\delta_{D}]_{n}\equiv\delta_{D}(\bm{k}-\sum_{i=1}^{n}\bm{q}_{i})$. For the correlation function this leads to,

[TABLE]

a result first obtained in [51]. If bias is not local, the functions of wave vectors that replace $b_{n}$ are just placed inside the integral in Eq. (16).

II.3 Renormalized bias expansion with

nonlinear evolution

Let us now consider the renormalized bias expansion in Eulerian space, which requires us to take into account the nonlinear evolution of the matter field. This computation serves two main objectives. Firstly, it highlights the complexity in determining contributions to the renormalized bias parameters from evolved fields, and thus motivates our approach taken in Sec. IV, where we instead evolve the multipoint propagators from their initial conditions. Secondly, it demonstrates that the bias expansion in Eq. (II) cannot be completely renormalized by the local parameters $b_{n}$ alone [see also 12], calling for additional terms in the expansion of $\delta_{g}$, which we address in Sec. III.

On large scales, where we can take the dark matter velocity field333Note that we work with the scaled velocity field $\bm{v}$, which is related to the peculiar velocity via $\bm{v}\equiv-\bm{u}/f{\cal H}$, with $f$ the logarithmic growth rate and ${\cal H}$ the comoving Hubble rate. $\bm{v}(\bm{x})$ to be irrotational, the time evolution is governed by the coupled equations of motion for the density perturbations and velocity divergence $\theta(\bm{x})\equiv\bm{\nabla}\cdot\bm{v}(\bm{x})$ under influence of the gravitational potential (see [64] for a detailed review). In standard Eulerian perturbation theory (SPT) these equations are solved as series expansions about the linear density field $\delta^{(1)}$, which we still assume to be Gaussian. At time $\tau$ we have (neglecting transients)

[TABLE]

where, to very good accuracy, all cosmology dependence is encoded in the linear growth factor $D(\tau)$ [65, 66, 19]. Note that what we denote as $\delta_{L}$ is actually the linear density field extrapolated to some final time (whose argument we will usually drop), i.e. $\delta_{L}(\bm{x})\equiv D(\tau)\,\delta^{(1)}(\bm{x})$. In Fourier space the $n$-th order solution is constructed out of $n$ powers of the linear density field, which are coupled via the SPT kernels $F_{n}$:

[TABLE]

The velocity divergence can be expanded in a similar manner and the $n$-th order solutions are obtained by replacing $F_{n}$ with $G_{n}$ in the equation above. Explicit expressions for these kernels can be found in [67, 68, 64]. We also note that at linear order, we have $\delta=\theta$.

Due to the mode-coupling in Eq. (II.3) the nonlinear (or time-evolved) density field becomes non-Gaussian. In order to retain the usual linear (as opposed to nonlinear) spectra in the loop integrals when computing clustering statistics, we however still define the multipoint propagators with respect to the linear Gaussian matter fluctuations, i.e. derivatives are taken with respect to $\delta_{L}$. However, this implies that the kernels $(\partial^{n}\delta_{g}/\partial\,\delta_{L}^{n})_{0}$ are not just numbers but acquire a scale dependence, in which case it becomes more convenient to express the multipoint propagators in Fourier space,

[TABLE]

as was already indicated in Fig. 1. While we present the resulting Fourier space equivalent of Eq. (12) in Sec. IV.1, in what follows it suffices to notice that the Fourier space multipoint propagators are functions of the $n$ momenta $\bm{k}_{1},\ldots,\bm{k}_{n}$ and include the factor $(2\pi)^{3}\,\delta_{D}(\bm{k}-\bm{k}_{1\ldots n})$. Primed angle brackets indicate that we have dropped this common factor.

Returning now to the calculation of the renormalized linear bias parameter, we see that because of second-order SPT (2SPT) we have an effective cubic term inside the quadratic bias contribution, $\bar{b}_{2}\,\delta^{2}/2$, from Eq. (II). Accordingly, the one-point propagator has the following extra term,

[TABLE]

where $\ast$ denotes a convolution and $\bar{F_{2}}$ is the angular average of the 2SPT kernel (i.e. the contribution from spherical collapse to $F_{2}$). This leads to $b_{1}=\bar{b}_{1}+(34\,\bar{b}_{2}/21+\bar{b}_{3}/2)\,\sigma^{2}+\ldots$ [7]. More broadly, since each term in the bias expansion of Eq. (II) appears at all orders in perturbation theory when $\delta$ corresponds to the Eulerian matter field, the renormalization of each bias parameter depends on all $\bar{b}_{n}$ ($n\geq 2$). Thus, the simplicity of Eqs. (9-10) is lost.

Let us consider what that implies for the renormalization of the quadratic bias, where we now have an effective quartic term (again due to 2SPT) inside the $\bar{b}_{3}\,\delta^{3}/3!$ contribution in Eq. (II), extending the two-point propagator by

[TABLE]

with $\delta^{L}_{i}\equiv\delta_{L}(\bm{k}_{i})$. The second term depending on the full $F_{2}$ kernel simply gives the desired contribution to the renormalization of the linear bias parameter by $\bar{b}_{3}$ as found for the one-point propagator before. That is because

[TABLE]

such that when calculating $\langle\partial^{2}\delta_{g}/\partial\delta^{L}_{1}\partial^{L}_{2}\rangle$ the second term in Eq. (22) corresponds to changing $\bar{b}_{1}\to\bar{b}_{1}+\bar{b}_{3}\,\sigma^{2}/2$ in Eq. (23). The first term in Eq. (22) on the other hand is the contribution to the renormalization of the quadratic bias parameter from nonlinear evolution, which now reads $b_{2}=\bar{b}_{2}+(68\,\bar{b}_{3}/21+\bar{b}_{4}/2)\,\sigma^{2}+\ldots$. In addition, quadratic bias contributes to the renormalization of $b_{2}$ itself, since

[TABLE]

The last term in this expression dresses the external lines to include the one-loop propagator due to nonlinear evolution (i.e. the standard $P_{13}$ contribution from the one-loop matter power spectrum [64]), while the remaining two terms give rise to the renormalization of the quadratic bias we are interested in, and they read

[TABLE]

with ($\mu\equiv\hat{\bm{k}}_{1}\cdot\hat{\bm{k}}_{2}$, $K(\bm{k}_{1},\bm{k}_{2})\equiv\mu^{2}-1$, and ${\cal L}_{\ell}$ denoting Legendre polynomials)

[TABLE]

While the integral over the function ${\cal F}$ in Eq. (26) represents the finite part (going as $1/q^{2}$ as $q\to\infty$) of the one-loop two-point propagator due to quadratic bias, the three terms proportional to $\sigma^{2}$ in the second line of Eq. (25) must be absorbed by renormalizations. We recognize that the first of these corresponds to a linear bias renormalization by $\bar{b}_{2}$, leading to $\bar{b}_{1}\to b_{1}=\bar{b}_{1}+(34\bar{b}_{2}/21+\bar{b}_{3}/2)\,\sigma^{2}$ in Eq. (23), which is consistent with the result found for the one-point propagator in Eq. (21). The second term renormalizes the quadratic bias, giving $b_{2}=\bar{b}_{2}+(8126\,\bar{b}_{2}/2205+68\,\bar{b}_{3}/21+\bar{b}_{4}/2)\,\sigma^{2}$ [12]. Finally, the third term has a different scale dependence compared to the previous ones, encoded by $K(\bm{k}_{1},\bm{k}_{2})$, such that it cannot be absorbed by a redefinition of any local bias parameter $b_{n}$. This illustrates that nonlinear evolution gives rise to additional effects that enter the relation between the galaxy and matter density and thus renders the local-in-matter expansion incomplete. In particular, we will see in Appendix D that the term proportional to $K(\bm{k}_{1},\bm{k}_{2})\,\sigma^{2}$ corresponds to a renormalization of the nonlocal quadratic bias $\gamma_{2}=\bar{\gamma}_{2}+127\,\bar{b}_{2}\sigma^{2}/2205$ [12], which is associated to the gravitational tidal field.

As we alluded to at the beginning of this section, compared to the calculation when the bias relation is written in Lagrangian space, the steps taken here are a lot more complicated (and we just scratched the surface, since we only assumed local galaxy bias so far). On the other hand, if we would like to maintain the simplicity of the Lagrangian expansion we would instead have to evolve the Lagrangian propagators by nonlinear evolution, which is where the complications may surface again. However, this is not the case for the main reason that the evolution from Lagrangian to Eulerian bias by conserving the number of objects cannot generate unphysical terms going like $\sigma^{2}$. Therefore, the time-evolved propagators do not require extra renormalizations from nonlinear dynamics, and we obtain the ‘finite parts’ such as Eq. (26) automatically.

It is further worth noting that the evolved propagators connect the Lagrangian bias parameters with the (late time) Eulerian ones, i.e. they give us the time evolution of the observable bias parameters, as opposed to the evolution of the bare bias parameters. We discuss this in some detail in Appendix D, as to whether renormalization affects the time evolution of bias parameters.

III A complete Galilean invariant basis for galaxy bias

We have already seen how the renormalization of local quadratic bias, when applied in Eulerian space, requires the existence of at least one additional term in the bias expansion, which was not initially included in Eq. (II). That is because gravitational evolution automatically generates a variety of terms nonlocal in $\delta$, such as the tidal field, which in principle can affect the formation of galaxies and should consequently be taken into account in the perturbative description of galaxy bias [69, 29, 8, 9, 12].

We will refer to any such terms (both local and nonlocal in $\delta$) as operators, and it is our aim in this section to provide a complete basis of operators required by the bispectrum at one-loop order. By basis we mean a set of linearly independent operators at each order of perturbation theory. We will largely follow up on the earlier work of [9, 12, 13, 6], but distinguish between operators local in second derivatives of the linear potential and those which are not. The former will be generated even if nonlinear evolution is entirely local (Zel’dovich approximation), whereas the latter derive from corrections at second-order Lagrangian perturbation theory (LPT) and beyond. Accordingly, we denote these as “local evolution” (LE) and “nonlocal evolution” (NLE) operators.

Our choice of basis will make this distinction explicit and therefore differs in the type of operators from that given in [8, 70, 13, 6]. The rationale for this choice is:

1. to give the multipoint propagators a particularly simple form, and 2) to provide a hierarchy of approximations based on which one can reduce the total number of bias parameters. While the usual local Lagrangian bias approximation (where only LE operators which are also local in $\delta$ are present at the initial conditions) may well be too restrictive for cosmological parameter estimation, it might prove useful in practice to relax this at least to all LE operators. This approximation is more accurate but still reduces the number of free bias parameters compared to the case when we also allow for NLE operators at the initial conditions. Other bases in the literature mix our two sets of operators, so one might not be able to see these subtleties when measuring bias parameters from simulations or data. We discuss their relation to ours in Appendix A.2.

III.1 Galileons as general basis operators

Let us begin with two physical scales important for the process of galaxy formation: 1) the spatial extent $R_{*}$ on which this process depends on the precise distribution of matter, and 2) the typical time $T_{*}$ it takes for this matter distribution to collapse into a bound object. While the latter is a significant fraction of the Hubble time $H^{-1}$, the scale $R_{*}$ usually corresponds to the Lagrangian radius of the galaxy’s host halo, which is of the order $\sim 1~{}\text{Mpc}$. If we are interested in the clustering of galaxies on scales $r\gg R_{*}$, then we can consider galaxy formation as essentially local in space. For now we will take this to be the case, before relaxing this assumption in Sec. III.4.

As stated at the beginning of this section, other properties than density of the matter field, such as the tidal field, must enter the bias relation [69, 29, 8, 9]. More generally, we can argue that these should depend on the (scaled) gravitational and velocity potentials, defined by

[TABLE]

as these drive the time evolution in the regime where the dark matter flow is irrotational. According to the equivalence principle all leading local gravitational effects must stem from second derivatives, which we write as $\nabla_{ij}\Phi(\bm{x},\tau)\equiv\partial_{i}\partial_{j}\Phi(\bm{x},\tau)$. Similarly, if we assume that dark matter and galaxies are comoving (i.e. no velocity bias), Galilean invariance of the equations of motion [71] implies that only second derivatives of the velocity potential are allowed to appear.

Furthermore, as $\delta_{g}$ is a scalar and therefore invariant under spatial coordinate transformations, we can limit ourselves to all scalar invariants of the tensors $\nabla_{ij}\Phi$ and $\nabla_{ij}\Phi_{v}$. In three dimensions the Cayley-Hamilton theorem [72] guarantees that there can only be three such invariants, which can be expressed by the so-called Galileons [9] (repeated indices are summed over):

[TABLE]

and similarly for $\Phi_{v}$. We note that their leading SPT expressions are of first, second, and third order, respectively. By inverting the Poisson equation (Eq. 27) we can derive their Fourier space analogs and for the latter two we obtain:

[TABLE]

where we have defined the following two kernel functions,

[TABLE]

with $\mu_{ij}\equiv\bm{k}_{i}\cdot\bm{k}_{j}/k_{i}\,k_{j}$. Similar expressions hold for $\Phi_{v}$ by replacing $\delta$ with $\theta$ in Eqs. (32) and (33). Note that these kernels vanish as $k^{2}$ when $k=|\sum_{i}\bm{k}_{i}|\to 0$, and therefore $\langle{\cal G}_{2}\rangle=\langle{\cal G}_{3}\rangle=0$ [see e.g. 9, 12]. In addition $\langle\partial{\cal G}_{3}/\partial\delta\rangle=0$ since $L(\bm{k},\bm{q},-\bm{q})=0$, which means that ${\cal G}_{3}$ cannot contribute to the one-loop galaxy propagator, as we shall discuss in section IV.2.1 below. All these properties make using Galileons as basis functions better suited for calculations compared to other choices that do not obey these relations, e.g. [8].

III.2 Local evolution (LE) operators

Let us consider the bias relation on some initial time slice in the far past. In that case we are only dealing with linear fluctuations and the single degree of freedom is $\Phi_{L}$, as at linear order we have $\Phi=\Phi_{v}$. If we assume that objects can be identified by a local procedure on $\nabla_{ij}\Phi_{L}$, the only terms that can appear in the bias relation for objects at the initial conditions are ${\cal G}_{m}(\Phi_{L})$ ($m\leq 3$), such that there will be $n$ basis operators at $n$-th order in the expansion ($n\leq 5$), i.e.

[TABLE]

and we assign a free bias parameter to each of these operators. We stress that the corresponding tracer density at initial time will thus be a local function of $\nabla_{ij}\Phi_{L}$, which is similar in spirit with more phenomenological approaches, such as the peak and excursion set bias models  [3, 4, 73, 74]. In reality though, tracers are not identified at the initial time, but in the late-time, nonlinear field, which means that the basis in Eq. (36) will be insufficient if gravitational instability produces terms nonlocal in $\nabla_{ij}\Phi_{L}$. This is the case starting with second-order corrections to the Zel’dovich approximation [75], as we will see in Sec. III.3. Therefore, even when traced back to the initial conditions such terms can in principle not be ignored. If we do so, the nonlocal terms are still produced during the evolution process, but with their bias coefficients fixed in terms of the parameters associated to the operators in Eq. (36) [9]. For that reason, this provides us with a very useful approximation (with fewer free parameters) that can be tested against numerical simulations [76].

III.3 Nonlocal evolution (NLE) operators

Due to nonlinear evolution, the gravitational and velocity potentials will no longer be identical. Following our previous arguments, it thus seems obvious to simply double the number of operators at each order by including a set of Galileons for both, the evolved $\Phi$ and $\Phi_{v}$, and also allow for their combinations, in order to complete the basis in Eq. (36). Unfortunately, this produces a lot of redundancy as many of these operators are degenerate, so our task will be to identify those, which are linearly independent. We follow the strategy first developed in [9].

At first order in SPT we have already established that ${\cal G}_{1}(\Phi)={\cal G}_{1}(\Phi_{v})$, and we choose the former, i.e. the matter fluctuation itself, as our first basis operator. Likewise, the two second-order Galileons are degenerate at second order in SPT and furthermore, we have

[TABLE]

proving that the basis in Eq. (36) is complete up to that order. The need for an additional operator occurs for the first time at third order. Using the notation $\Delta_{n}{\cal G}_{m}\equiv{\cal G}_{m}^{(n)}(\Phi)-{\cal G}_{m}^{(n)}(\Phi_{v})$ for the difference between the $m$-th Galileons evaluated at $n$-th order in SPT, we see that in addition to the ones already written in Eq. (36) there are the following four combinations:

[TABLE]

From Eq. (37) follows that the latter two are degenerate with $\delta\,{\cal G}_{2}(\Phi_{v})$, while $\Delta_{3}{\cal G}_{1}$ contains a contribution that cannot be written in terms of second derivatives of $\Phi_{v}$ and is thus not Galilean invariant. This only leaves the second combination, which gives

[TABLE]

demonstrating, as claimed above, that the additional basis operators induced by gravity can no longer be expressed as local functions of the linear gravitational potential, i.e. $\nabla_{ij}\Phi_{L}$. As these effects are precisely captured by second-order LPT and beyond, instead of explicitly calculating the differences between Galileons of $\Phi$ and $\Phi_{v}$, we can follow an alternative strategy (which builds on [9]) that will prove particularly useful for extending the basis beyond third order. Let us consider LPT, which summarizes all of the dynamics in its Lagrangian displacement field

[TABLE]

that moves particles from their initial positions $\bm{q}$ to their final destinations $\bm{x}=\bm{q}+\bm{\Psi}(\bm{q},\tau)$ at conformal time $\tau$. The functions $\bm{\Psi}^{(n)}(\bm{q})$ are the $n$-th order contributions and $D_{n}(\tau)$ are the corresponding growth factors ($D_{1}\equiv D$ is the linear growth factor). At any order both the gravitational and velocity potentials can be expressed in terms of these $\bm{\Psi}^{(n)}$, which are in turn given by the LPT potentials $\varphi_{n}$ [77, 66], e.g.

[TABLE]

Any set of linearly independent operators induced by gravity must therefore be connected to combinations of the LPT potentials. In order to guarantee that these are still Galilean invariant, we can generalize the definition of the Galileons to [9]

[TABLE]

and similarly for ${\cal G}_{3}(A,B,C)$. From Eq. (41) we have $\varphi_{1}=-\Phi_{L}$, so that the first new combination appears at third order of perturbation theory: ${\cal G}_{2}(\varphi_{2},\varphi_{1})$. Evaluating this Galileon using Eq. (42) shows that it is precisely related to the only gravity induced operator that we previously identified at third-order, i.e. [9]

[TABLE]

Following this line of argument we can now easily determine the additional operators at fourth order: apart from the combination $\delta\,{\cal G}_{2}(\varphi_{2},\varphi_{1})$, we can construct the following three invariants out of the LPT potentials

[TABLE]

However, beyond second order the LPT solutions are no longer purely potential anymore and at third order in particular it consists of two scalar and a transverse vector potential, all with different time dependencies:

[TABLE]

where in the Einstein-de Sitter (EdS) approximation the growth factors are given by $D_{3}^{(a)}(\tau)=1/18\,D(\tau)^{3}$, $D_{3}^{(b)}(\tau)=5/42\,D(\tau)^{3}$ and $D_{3}^{(c)}(\tau)=1/14\,D(\tau)^{3}$. The functions $\Psi^{(3,a/b/c)}$ and their corresponding potentials satisfy the following Poisson equations [78]

[TABLE]

where $\epsilon_{ijk}$ denotes the fully anti-symmetric Levi-Civita symbol and $\hat{e}_{i}$ the unit vector in direction $i$. The combination of the third and first order LPT potentials is thus made up of three pieces and factoring out $D(\tau)^{3}$ from the EdS solutions, we define

[TABLE]

Due to the different time dependencies we should in principle allow these three pieces to enter the bias basis individually, but in practice the departures from the EdS approximation to growth factors are below the percent-level, so for all purposes of this paper we are safe to ignore this complication.

To conclude, our choice of a complete Galilean invariant basis for the evolved galaxy perturbations is given by a set of 15 operators up to fourth order, which are summarized in the first two columns of Table 1. We have separated what we denote as local evolution operators (first column), which are local in $\nabla_{ij}\,\Phi_{L}$, from those that correspond to nonlinear corrections to the gravitational and velocity potentials that are nonlocal in $\nabla_{ij}\,\Phi_{L}$ (middle column). With respect to the number of basis operators we are thus in agreement with [6], the relation between our and their basis up to fourth order is given in Appendix A.2.

III.4 Higher-derivative galaxy bias

Although some of the basis operators we derived in Sec. III.2 and III.3 are nonlocal in the matter fluctuations and gravitational potential, we made the central assumption that the formation of galaxies is spatially-local, meaning it is determined by the value of these operators at a single point in space. On small scales this approximation must break down because galaxies form due to matter collapsing from a finite region of size $R_{*}$, which is of the order of the Lagrangian radius of the host dark matter halo. As shown in [8, 70, 6] this gives rise to additional operators in the bias expansion that contain higher than second derivatives of the gravitational and velocity potentials444These have been known for a very long time, e.g. [3] since they naturally appear in peak models of biased tracers..

As these operators must still be scalars, the simplest one involves four derivatives of $\Phi$ and is given by $\nabla^{2}\,\delta$. Its effect can be interpreted as an emerging scale dependence of the linear bias parameter, because upon Fourier transformation and grouping all terms linear in $\delta$, we get

[TABLE]

where $\beta_{1}$ is the bias parameter associated with $\nabla^{2}\,\delta$ and ${\cal O}(k^{4})$ stands for contributions involving even higher derivatives. Note that $\beta_{1}$ has units of length squared, and as we can expect it to scale with $R_{*}$, we see that on scales much larger than the Lagrangian radius, i.e. $k\,R_{*}\ll 1$, the higher-derivative contributions are suppressed by powers of $(k\,R_{*})^{2}$ compared to the spatially-local linear term. This is similar to the nonlinear bias terms, where an increase in nonlinear order is suppressed by powers of $\delta\simeq(k/k_{\rm nl})^{(n_{\rm eff}+3)/2}$ at large scales, with $k_{nl}^{3}P(k_{nl})/(2\pi^{2})\equiv 1$ and $n_{\rm eff}$ the effective spectral index at the nonlinear scale. However, comparing the nonlinear and higher-derivative contributions against each other a priori is difficult, since this depends on the relative size of $R_{*}$ (which in turn depends on the biased tracer in question) and $k_{nl}$ and $n_{\rm eff}$ (which depend on the linear spectrum and redshift), in addition to the size of the bias parameters. We therefore take the following approach: we consider a priori each higher-derivative factor in correlators as a nonlinear bias loop, counting each additional derivative acting on the gravitational or velocity potentials as an increase of the SPT order by one. The operator $\nabla^{2}\,\delta$ would thus be considered as third order. After performing a measurement of bias parameters from clustering data, one can reassess a posteriori whether bias loops of derivatives are more important for a particular biased tracer. In [31], the authors conclude that higher-derivative biases are more important than loops for tracers of a very wide range of masses. We revisit this issue in [76], but Figs. 4 and 5 below already suggest that loop corrections are as important for the bispectrum as for the power spectrum and that they matter on scales commonly used in the analysis of galaxy surveys.

Based on our counting, the fourth-order higher-derivative terms will simply be given by acting with two derivatives on the spatially-local quadratic bias operators. This allows only for the following four independent combinations

[TABLE]

where $\delta\,\nabla^{2}\,\delta$ as well as ${\cal G}_{2}(\Phi_{v},\nabla^{2}\Phi_{v})$ do not enter individually, as they can be expressed through combinations of the operators above. In total we thus obtain a set of five higher-derivative operators, which are summarized in Table 1 according to their classification in terms of SPT order as discussed above. Note that Ref. [6] gives a total of ten operators up to the same order. We find that the set presented here is equivalent with the first five operators in their Eq. (2.74), while the last five can be expressed in terms of the other ones, making them superfluous.

IV Multipoint propagators in Lagrangian and Eulerian space

We now return to the main idea presented in Sec. II.2: in order to guarantee that the bias parameters corresponding to all of the Galilean basis operators are observable quantities, we should construct the galaxy density field out of multipoint propagators,

[TABLE]

as opposed to the usual kernel functions that we obtain from $\left[\partial^{n}\delta_{g}(\bm{k})/\partial\,\delta_{L}(\bm{k}_{1})\cdots\partial\,\delta_{L}(\bm{k}_{n})\right]_{0}$. Our goal is to compute the complete multipoint propagators at initial and final time, while paying particular attention to the evolution of the bias parameters. Before that we illustrate their relation to the Wiener-Hermite functionals, which formalizes the multipoint propagator expansion in Fourier space and will prove useful for determining their time evolution.

IV.1 Relating multipoint propagators and the Wiener-Hermite expansion

For a linear Gaussian dark matter density field, the probability density function (PDF) for a mode in Fourier space is given by

[TABLE]

with normalization factor $N$. The $n$-th generalized Wiener-Hermite functional ${\cal H}_{n}$ is then defined by taking $n$ functional derivatives of the PDF [79]:

[TABLE]

Using this definition, we obtain the following first three functionals (suppressing the momentum arguments):

[TABLE]

where a superscript $*$ denotes complex conjugation. Like the standard Hermite polynomials they satisfy an orthogonality relation, which can be shown to be (see [79] and App. B.1):

[TABLE]

where $\delta_{nm}^{\rm K}$ is the Kronecker delta, and “sym.” stands for the remaining $(n!-1)$ combinations of the arguments. Expanding the Fourier space galaxy overdensity in terms of Hermite functionals gives a convolution over ${\cal H}_{n}$ at each order, such that

[TABLE]

where the coefficients $\Gamma_{g}^{(n)}(\bm{k}_{1},\ldots,\bm{k}_{n})$ are scale-dependent functions (kernels), that can be interpreted as the corresponding bias parameters. Multiplying both sides of Eq. (IV.1) with ${\cal H}_{m}$ and using the orthogonality relation (Eq. IV.1), we see that

[TABLE]

which generalizes Eqs. (14-15). Cross-correlating the galaxy density field with Hermite polynomials to measure local bias parameters has been proposed in a similar way in [80]. Moreover, this is also how propagators for the mapping from linear to nonlinear fluctuations have been measured in [50, 41, 43]. Note that in general for nonlocal bias this cross-correlation has to be further decomposed in terms of the structures inside $\Gamma_{g}^{(n)}$ to end up with parameter estimates (i.e. one has to separate $b_{2}$ from $\gamma_{2}$ contributions etc.), e.g. [81, 23, 82, 60, 61], but that procedure is basis-dependent.

We can derive a different relation by plugging in Eq. (55) into Eq. (IV.1) and replacing the ensemble average by its definition — the functional integral of all modes $\delta_{L}$ over their joint PDF:

[TABLE]

where we have integrated by parts $n$ times in going from the second to the third line and ignored any surface terms as ${\cal P}[\delta_{L}]\to 0$ (same as its derivatives) for $\delta_{L}\to\pm\infty$. Thus, it follows that the kernels of the Wiener-Hermite expansion are indeed the multipoint propagators:

[TABLE]

This equation is equivalent to the definition of multipoint propagators in the previous work of [41, 43], which considered the nonlinear matter fluctuations in place of the galaxy fluctuations. As they showed that the multipoint propagators function as basic building blocks for constructing arbitrary $N$-point spectra, the same holds for the $\Gamma_{g}^{(n)}$ and we will use this fact to compute the galaxy power spectrum and bispectrum in Sec. V.

IV.2 Initial conditions

We now determine the first three multipoint propagators at an initial time where nonlinearities in the dark matter field can be ignored, which implies we can set all SPT kernels $F_{n}$ and $G_{n}$ for $n\geq 2$ to zero. We include all of the basis operators summarized in Table 1 and work up to fourth order in SPT, as required by the one-loop bispectrum. For reasons explained below, the general structure of the multipoint propagators in our basis is very simple: they are given by a tree-level contribution consisting of all basis operators that correspond to the order of the propagator itself, in addition to loop corrections involving only NLE operators.

IV.2.1 The one-point propagator

To compute the one-point propagator, we need to take a single derivative of $\delta_{g}$ and take the expectation value, which implies due to Gaussianity of $\delta_{L}$ that only odd orders in the bias expansion enter. More generally, as each derivative cancels exactly one factor of $\delta_{L}$, we see that odd (even) numbered propagators can only contain terms stemming from odd (even) orders of the bias expansion.

The first term in the bias expansion gives just linear bias, thus to compute $\Gamma_{g}^{(1)}$ to one loop we need the derivative of a generic third order term ${\cal O}_{B}^{(3)}$, which can be written as the convolution

[TABLE]

with $B\in\{\delta^{3},\,\delta\,{\cal G}_{2},\,{\cal G}_{3},\,{\cal G}_{2}(\varphi_{2},\varphi_{1})\}$ and the kernels ${\cal K}_{B}^{(3)}$ are given in Eqs. (152) to (155). Making use of the symmetry of ${\cal K}_{B}^{(3)}$, we obtain:

[TABLE]

using the same notation for ensemble averages as in Sec. II.3. We note that Eq. (62) is analogous to the result obtained for the one-loop contribution to the matter propagator, which gives rise to $P_{13}/2P_{L}$ [40].

As pointed out already, most of the possibilities for third-order kernel contributions are trivial or zero in our choice of basis functions for bias. In fact, no basis function local in $\nabla_{ij}\Phi$ (denoted as LE operators in Table 1) can give a non-trivial contribution to the loop corrections of $\Gamma_{g}^{(1)}$ (with trivial contributions meaning renormalizations of the linear bias parameter, as discussed in section II.2). The reason is that in propagator loop integrals no wavevector angles can appear inside power spectra (due to being reducible diagrams), thus one can always perform angular integrations over momenta of the resulting kernels, which for operators local in $\nabla_{ij}\Phi$ simply give rise to numerical prefactors (many of them zero in our choice of basis). Therefore such loops are only functions of $\sigma^{2}$, e.g. for $B=\delta\,{\cal G}_{2}$,

[TABLE]

similarly to the case $B=\delta^{3}$, while $\langle\partial{\cal G}_{3}/\partial\delta_{L}\rangle=0$ since $L(\bm{k},\bm{q},-\bm{q})=0$. On the other hand, for the only NLE operator (see Table 1) at third order $B={\cal G}_{2}(\varphi_{2},\varphi_{1})$, we get

[TABLE]

Collecting all these results, the first galaxy propagator at initial time is given by

[TABLE]

where

[TABLE]

and we have added the higher-derivative term. As noted in Sec. II.2, Eq. (66) can also be derived from the renormalized expansion in Eq. (13) without explicitly computing the loops for the terms local in $\nabla_{ij}\Phi$.

In the limit $k\to 0$, the $\gamma_{21}$ integral scales as

[TABLE]

which ensures that $\Gamma_{g}^{(1)}$ indeed corresponds to the linear bias parameter on large scales. Moreover, we note that if terms of order $k^{4}$ and beyond are negligible, the $\gamma_{21}$ contribution is entirely degenerate with the higher-derivative term [see also 8, 83, 24]. As long as we restrict ourselves to sufficiently large scales, this implies it is not necessary to include the higher-derivative parameter in excess of $\gamma_{21}$. However, a priori it is difficult to tell what “sufficient” means for a given data set, and therefore requires careful testing using mock data.

IV.2.2 The two-point propagator

In a similar manner we can now derive all remaining multipoint propagators. The two-point propagator gets contributions from second and fourth order bias operators, which, analogous to Eq. (IV.2.1), we write as the symmetric kernels ${\cal K}_{B}^{(2)}$ and ${\cal K}_{B}^{(4)}$. The corresponding expressions are given in Eqs. (149) and (150), and Eqs. (159) to (166), respectively. Differentiating a generic second or fourth order contribution twice results in

[TABLE]

and

[TABLE]

As explained above, the non-trivial loop corrections to the two-point galaxy propagator can only consist of the NLE operators and so we get:

[TABLE]

where $k_{12}^{2}\equiv|\bm{k}_{1}+\bm{k}_{2}|^{2}$ and the bracket inside the loop integral is evaluated at $(\bm{k}_{1},\bm{k}_{2},\bm{q},-\bm{q})$. Furthermore, the kernel ${\cal K}_{\delta{\cal G}_{2}(\varphi_{2},\varphi_{1})}^{(4,\text{F})}$ denotes the finite part of ${\cal K}_{\delta{\cal G}_{2}(\varphi_{2},\varphi_{1})}^{(4)}$ (i.e. removing contributions proportional to $\sigma^{2}$) and we have added the higher-derivative terms in the last two lines. For the same reasons mentioned above, all of the initial operators again only contribute to renormalizations of

[TABLE]

These are the observable quadratic and tidal tensor bias parameters, which match those in Eqs. (3.16-3.17) in [12] up to the signs of the $\gamma_{21}^{\times}$ and $\gamma_{3}^{\times}$ terms in the expression for $\gamma_{2}$ (note that in their notation $b_{\Gamma_{3}\delta}=-4/7\,\gamma_{21}^{\times}$ and $b_{{\cal G}_{3}\delta}=\gamma_{3}^{\times}$). A more direct comparison with [12] including renormalization due to nonlinear time evolution is discussed in Appendix D. Again, the relations in Eqs. (71-72) can be obtained directly from Eq. (13) without computing the loops explicitly.

The large-scale behaviour of the nonlocal evolution terms that we found for the one-point propagator in the previous section also applies to the two-point propagator. More precisely, it can be shown that the four loop integrals in Eq. (IV.2.2) are fully expressable in terms of combinations of the four higher-derivative operators (and vice versa), and we present the corresponding relations in App. A.3. As noted above, this means that the higher-derivative operators become degenerate with 4th-order contributions from the general bias expansion, and their impact is automatically accounted for by the latter, as long as terms of order $k_{1}^{4}$, $k_{2}^{4}$, etc. are negligible.

IV.2.3 The three-point propagator

Finally, for the three-point propagator, we compute three derivatives of ${\cal O}_{B}^{(3)}$, giving

[TABLE]

At the order of SPT we are working in, we only need the tree-level expression for $\Gamma_{g}^{(3)}$ and we simply get:

[TABLE]

IV.3 Time evolution

As we already discussed in Sec. II.3 for the case of local galaxy bias, until the time when galaxies are observed, all of the basis operators will have evolved and thus developed nonlinear corrections that are of higher SPT orders. Consequently, if we were to compute the galaxy multipoint propagators at late times, we would have to account for these corrections, which lead to numerous new renormalizations (see discussion related to Eqs. 25-26). Especially when pushed to fourth order in SPT this approach becomes very cumbersome.

For that reason this section presents an alternative — we start from the multipoint propagator expansion at initial time using the propagators derived in Sec. IV.2 and then evolve this expansion instead of resorting to the usual SPT solutions. Because the mapping from the initial time to the final time conserves the number of tracers, no new divergences going like $\sigma^{2}$ can arise, and this bypasses the need to deal with these complexities. Therefore, the only bias renormalizations done in our approach are trivial ones, done at initial time, and in fact dealt simply by using the multipoint propagators: in this sense, in our approach one does not have to think about renormalization at all.

Since time evolution from the initial conditions also gives rise to an evolution of the bias parameters, we first consider them separately from the propagators. This allows us to illustrate precisely how various assumptions made about the initial bias relation affect the late-time galaxy overdensity, and thus generalizes previous results in the literature up to fourth order. In order to distinguish initial from evolved quantities, we will make the following changes to notation:

[TABLE]

for the initial bias parameters and propagators, respectively.

IV.3.1 Bias parameters

One of the key principles underpinning the bias relation is that it must be Galilean invariant if there is no velocity bias. We can exploit this symmetry to determine the evolution of the bias parameters in a significantly simpler way by restricting all quantities that appear in intermediate steps of the calculation to a Galilean invariant basis of operators, i.e. the basis we presented in Sec. III. To illustrate this novel technique we start from the evolution equations for conserved tracers, which in the absence of velocity bias can be directly integrated to give (see Eq. 48 in [9])

[TABLE]

where Lagrangian fields (with argument $\bm{q}$) are evaluated at the initial time when $\delta=\delta_{L}$. Since we are not interested in decaying modes, we can neglect $\delta^{\rm IC}(\bm{q})$ in the denominator because it is suppressed by one growth factor compared to the Eulerian fields, leading to the well-known expression for the evolution of bias [84]. Similarly, we can ignore decaying modes and write the initial bias relation in terms of the extrapolated (to final time) linear density fluctuations and Lagrangian bias parameters (identical to the parameters in Sec. IV.2),

[TABLE]

where all RHS fields are linear Gaussian evaluated at the final time. To proceed we would then have to use $\bm{q}=\bm{x}-\bm{\Psi}$ to relate initial to final positions and expand all quantities around $\bm{x}$. However, this can be done trivially by realizing that in order for the LHS of Eq. (76) to be Galilean invariant, so must be the RHS, which implies that the dipole terms that arise from the $\bm{q}$ to $\bm{x}$ mapping must cancel against those in the nonlinear matter density. All we need to do is thus to consider the Galilean invariant contributions to the nonlinear density and velocity expressed in terms of our bias basis, which is done in Appendix C.1. Indicating this by the superscript “GI”, it follows that

[TABLE]

and by using our bias basis at each order we can then find the coefficients (bias parameters) for the Eulerian expansion in terms of those in the initial expansion. For instance, at linear order we have $\delta_{g}^{\text{IC}}(\bm{x})=b_{1,{\cal L}}\,\delta_{L}(\bm{x})$ and plugging this into Eq. (78) gives

[TABLE]

which agrees with the result from full evolution [69] when neglecting transients. For higher-order bias evolution we proceed as in [9], i.e. to find bias parameters at a given order we subtract the contributions expected from the lower order bias parameters. The full details of this computation are given in Appendix C.2 and at second order we get the following well known relations

[TABLE]

while we obtain

[TABLE]

for cubic bias (note that there was a typo in Eq. (116) of [9], which should have had a minus sign for its second term), and finally

[TABLE]

for quartic bias. Appendix C.3 shows how to recover the familiar spherically symmetric results from these expressions, which serves as a robust consistency check. Although the relations above have been derived for the bare bias parameters, they are equally valid for the renormalized ones. This is demonstrated explicitly in Appendix D, which considers what happens if one applies Eq. (78) to the renormalized bias expansion. However, as we mentioned above, in this paper we advocate a different, simpler route which evolves the initial multipoint propagators, as we discuss next.

IV.3.2 Propagators

As a first step we write a combined evolution equation for matter, its velocity divergence and galaxies, and for that purpose we define the three-component vector

[TABLE]

In terms of $\bm{\Psi}$ and by using the logarithm of the linear growth rate as our new time variable, i.e. $\eta\equiv\ln{D(\tau)}$, the evolution equations can be recast as (see [9])

[TABLE]

where we follow the convention that repeated indices are summed over, and we assume galaxies move with the dark matter, i.e. no velocity bias. The matrix

[TABLE]

describes the coupling between densities and velocities, while $\gamma_{abc}$ encodes the nonlinear interactions between different Fourier modes. Its only non-zero components are given by

[TABLE]

and $\gamma_{112}(\bm{k}_{1},\bm{k}_{2})=\gamma_{121}(\bm{k}_{2},\bm{k}_{1})$. Given some arbitrary initial conditions $\Psi_{a}(\bm{k},\eta=0)\equiv\phi_{a}(\bm{k})$, Eq. (IV.3.2) has the integral solution [85]

[TABLE]

This expression depends on the linear propagator $g_{ab}(\eta)$, which solves the linearized equations of motion (i.e. setting the righthand side of Eq. (IV.3.2) to zero) and presents a mixture of growing and decaying, as well as time independent modes [85, 9]:

[TABLE]

As we did for the galaxy overdensity in Sec. IV.1, we now expand $\Psi(\bm{k},\eta)$ in terms of generalized Wiener-Hermite functionals,

[TABLE]

where $\Gamma_{3}^{(n)}\equiv\Gamma_{g}^{(n)}$. An equivalent expansion holds at initial time $\eta=0$ and we denote the corresponding propagators by the symbol $\Gamma_{{\cal L}}$. Since the dark matter and velocity fields are linear at that time, we have $\Gamma_{1,{\cal L}}^{(1)}=1=\Gamma_{2,{\cal L}}^{(1)}$, while all higher-order propagators must vanish, i.e.

[TABLE]

The initial galaxy propagators, on the other hand, are given by the expressions from Sec. IV.2. As mentioned in the introduction to Sec. III, it is here where one might take advantage of some simplifying assumptions, for instance, removing fourth-order NLE operators from the initial conditions. This sets

[TABLE]

in Eqs. (83-84), which does not imply that we ignore such NLE operators completely because, as demonstrated by these equations, time evolution will generate amplitudes $\gamma_{21}^{\times},\gamma_{31},\gamma_{22},\gamma_{211}$ fixed through lower order bias parameters (note this approximation is more general than the local Lagrangian approximation, which sets all nonlocal in $\delta$ Lagrangian bias parameters to zero).

By multiplying both sides of Eq. (IV.3.2) with ${\cal H}_{n}$ and taking the ensemble average, we can derive a recursion relation for the time evolved multipoint propagators. Following the steps detailed in App. B, we arrive at the expression:

[TABLE]

where the quantity $\Gamma^{(n,m,r)}_{a}$ represents $r$ loop integrals over propagators of orders $m$ and $n-m+2r$:

[TABLE]

with $\text{sym.}(\bm{k}_{i})$ standing for the symmetrization over all possibilities of building a subset of $m-r$ $\bm{k}$-modes from a total group of $n$, i.e. $\binom{n}{m-r}$ terms.

Furthermore, Eq. (IV.3.2) illustrates the point we made at the beginning of this section — the time evolved multipoint propagators are free of potentially divergent contributions proportional to $\sigma^{2}$, which means we do not have to perform any additional renormalization steps. This is a consequence of the scale dependence of the mode coupling kernels $\gamma_{abc}$, which contribute to $\Gamma_{g}^{(n)}$ only via the symmetrized combination $\alpha^{s}(\bm{k}_{1},\bm{k}_{2})=1/2[\alpha(\bm{k}_{1},\bm{k}_{2})+\alpha(\bm{k}_{2},\bm{k}_{1})]$ coming from the continuity equation (conservation of tracers).

Plugging in the wavevectors from Eq. (IV.3.2) and expanding in inverse powers of $q_{1\cdots r}$, we find that to leading order

[TABLE]

where $\bm{K}=\bm{k}_{1\ldots m-r}$. In this paper we are interested in corrections only up to the one-loop level, so we are led to consider expressions of the form

[TABLE]

where both propagators are evaluated at tree-level. That means the propagators remain finite when the loop momentum $q$ becomes large, which in turn implies that the overall integrand scales, at most, as $1/q^{3}$ in this limit (using that $P_{L}(q)\sim 1/q^{3}$ for $q\rightarrow\infty$). This guarantees that any integral of the above type is quickly convergent.

For higher than one-loop corrections this argument no longer holds, as each loop adds an additional power spectrum to the integrand, while the scaling with the loop momenta remains the same. However, the scale dependence on the external momenta, i.e. $\sim k_{1\cdots n}^{2}$, suggests that these terms can be absorbed by redefinitions of the higher-derivative bias parameters, which display the exact same scale dependence, as shown in Sec. III.4. Using higher-derivative terms to absorb sensitivities to the highly nonlinear regime was already presented in a slightly different manner in [86, 87]. We note that this behavior of the loop integrals is well known in the context of matter perturbations [43, 88, 89, 90, 91], where the small scale sensitivity can be understood in terms of non-zero stress tensor corrections. Indeed, we will show in Sec. V.2 that these are completely degenerate with the contributions from higher-derivative bias.

Using Eq. (114) it is straightforward to compute the evolved multipoint propagators by starting from lowest order at tree-level and constructing all higher-order solutions recursively. One part of the of the solution is due to the evolution of the bias parameters and can be simply obtained by replacing the Lagrangian parameters in the initial multipoint propagators by their Eulerian analogs (according to the relations given in Sec. IV.3.1). Here we focus on the remaining corrections from nonlinear evolution, which we write as $\Delta\Gamma_{g}^{(n)}$, and neglecting all but the fastest growing mode and expressing the result back in terms of the SPT kernels, we find that these corrections for the one-point propagator at one-loop ($r=1$) are given by

[TABLE]

The two- and three-point propagators are already affected at tree-level by nonlinear evolution and the resulting corrections are

[TABLE]

For the one-loop bispectrum we also require nonlinear corrections of $\Gamma_{g}^{(2)}$ at one-loop order, which are, however, not easily expressed in terms of SPT kernels, so we give the full result in terms of the initial multipoint propagators in Appendix B.3.

V Power Spectrum and Bispectrum

This section details the final step of this paper: the computation of the galaxy power spectrum and bispectrum in terms of the multipoint propagators. In addition, we focus on residual sensitivities of the loop corrections on the highly nonlinear regime, where our perturbative approach breaks down. These affect the statistics not only on small scales, but also on asymptotically large scales, and we discuss how they can be regularized by the addition of physically motivated terms.

V.1 Reconstructing correlators from

multipoint propagators

The multipoint propagators serve as the basic building blocks for computing $N$-point spectra. This follows easily from the orthogonality relations of the generalized Wiener-Hermite functionals (see App. B.1), and was already shown for the dark matter density in [41], whose results we can directly apply to the present case of galaxy clustering.

In particular, by evaluating $\langle\Psi_{a}(\bm{k})\Psi_{b}(\bm{k}^{\prime})\rangle$ and using Eq. (IV.1) one finds that the power spectrum is given by a series of two contracted multipoint propagators of the same order. Diagramatically this can be represented by glueing together two of the objects shown in Fig. 1, where each combination of the incoming lines gives rise to a (linear) power spectrum. This is demonstrated in Fig. 2 and as the shaded circles include all vertex loop corrections (i.e. vertex renormalizations) we only need to consider one distinct diagram for the power spectrum at one-loop level, compared to the usual two in the standard treatment. The galaxy power spectrum is thus given by [41]

[TABLE]

where $\Gamma_{g}^{(1)}$ is evaluated up to one-loop order, while tree-level terms are sufficient for $\Gamma_{g}^{(2)}$.

We proceed in a similar manner for the bispectrum, which is obtained from $\langle\Psi_{a}(\bm{k}_{1})\Psi_{b}(\bm{k}_{2})\Psi_{c}(\bm{k}_{3})\rangle$ and application of Eq. (B.1). The complete solution can be found in [41], but its diagrammatic depiction in Fig. 3 is straightforward — a combination of three multipoint propagators with a varying number of connecting lines between each pair (maximally one such pair is allowed to be disconnected to avoid an overall disconnected graph). Taking care of the appropriate symmetry factors that arise in these various combinations, it follows that [41]:

[TABLE]

Therefore, for the one-loop galaxy bispectrum we require both, $\Gamma_{g}^{(1)}$ and $\Gamma_{g}^{(2)}$, up to one-loop order in the first term of Eq. (V.1), but tree-level expressions for them and $\Gamma_{g}^{(3)}$ are enough in the loop integrals, i.e. in the second and third line. Note that there are only two one-loop diagrams instead of the four in SPT; this is because there are two diagrams in SPT that are reducible, and thus incorporated into the one-loop $\Gamma_{g}^{(1)}$ (giving the 321-II diagram in the notation of [92]) and the one-loop $\Gamma_{g}^{(2)}$ diagram (giving the 411 diagram).

From comparing Eqs. (V.1) and (V.1) we note that the two-point propagator contributes to the leading order bispectrum, while showing up as a loop correction for the power spectrum. This structure extends to consecutively higher orders, for instance, the three-point propagator which appears as a one-loop expression in the bispectrum, will enter at tree-level for the trispectrum. That suggests that constraints on bias parameters required to fit the small-scale behaviour of a given correlator (and thus cosmological parameters, too) will already highly benefit from the large-scale information of the next order correlator.

When splitting each expression in Eq. (V.1) into its individual contributions, one finds that in total there are 40 additional terms in the one-loop galaxy bispectrum compared to tree-level. Each of these terms is multiplied by a combination of various bias parameters, making a comparison of the individual terms not particularly meaningful. However, given a set of bias parameters, it is interesting to consider whether the bispectrum loop contributions become relevant on similar scales as for the power spectrum. To that end, in Fig. 4 we show a comparison of the one-loop galaxy power spectrum evaluated from Eq. (V.1) and its tree-level prediction $P_{g}^{\mathrm{tree}}=b_{1}^{2}\,P_{L}$, where different colors correspond to a different choice of the linear bias parameter. In order to adopt representative values for the higher-order bias parameters, we fix them in terms of $b_{1}$ by making use of the local Lagrangian approximation (obtained by setting all parameters with subscript ’${\cal L}$’ in Eqs. (80-84) to zero), as well as the peak-background split relations for $b_{2}$ and $b_{3}$.555Note that the power spectrum and bispectrum models in Fig. 4 and 5 include the stochastic corrections in the low-$k$ limit to be discussed in Sec. V.3, but with all noise parameters set to zero. These were calibrated using separate universe simulations in [93], yielding

[TABLE]

As is demonstrated by the lower panel of Fig. 4, the one-loop corrections to the galaxy power spectrum become relevant on scales $k\gtrsim 0.1\,h/\text{Mpc}$. A similar trend can be observed for the galaxy bispectrum, whose relative difference between the tree-level and one-loop models for $b_{1}=1.8$ is plotted in Fig. 5. To show its configuration dependence we evaluate the bispectrum as a function of $x_{2}=k_{2}/k_{1}$ and $x_{3}=k_{3}/k_{1}$, averaged over a thin shell centered at a given value of $k_{1}$, i.e.

[TABLE]

where we have chosen $\Delta k=0.02\,h/\text{Mpc}$ and three different values for $k_{1}$ (see Fig. 5). We notice that for all triangle configurations the relative difference surpasses $10\,\%$ for $k_{1}\gtrsim 0.1\,h/\text{Mpc}$, with the strongest deviations to be found for equilateral triangles. Of course, the precise numerical values obtained in this demonstration depend sensitively on the set of bias parameters, but we do not expect any big impact on the overall conclusion: bias loop corrections affect the power spectrum and bispectrum starting from comparable scales.

V.2 Stress-tensor corrections and their degeneracy with higher-derivative bias

While the introduction of multipoint propagators has automatically removed potentially divergent contributions proportional to $\sigma^{2}$ (i.e. $(k/q)^{0}$ with $k$ the external momentum and $q$ the loop momentum), the loop integrals remain sensitive to the nonlinear regime through terms scaling as powers of $(k/q)^{2}$ (and $(k/q)^{4}$, etc). This is true for loop corrections of individual multipoint propagators, as discussed in Sec. IV.2 and at the end of Sec. IV.3, but also for the loop integrals over tree-level propagators appearing in Eqs. (V.1) and (V.1). Sensitivity of loop integrals to the nonlinear regime is obviously a problem since perturbation theory does not hold at small scales.

At such small scales the dark matter field can no longer be treated as a pressureless perfect fluid, as is assumed in SPT. Initially, or on large scales, this is a good approximation as dark matter particles tend to move within single coherent flows, which implies a vanishing stress tensor $\sigma_{ij}$. At later times, multi-streaming induces non-zero stresses, which can have an impact on quasi-linear scales, with the same $(k/q)^{2}$ scaling as one-loop SPT and one-loop bias terms [32].

These stress tensor corrections have been computed in the framework of the effective field theory (EFT) [34, 35, 36, 37] and directly from the Vlasov equation [94] following [32, 95]. This leads to additional terms in the power spectrum and bispectrum that at lowest order scale with powers of $k^{2}$, identically to one-loop corrections sensitive to the nonlinear regime. In the particularly economic parametrization of [94] (which is otherwise equivalent to the EFT calculations mentioned above) they read,

[TABLE]

where $\beta_{P}$ and $\beta_{B,i}$ are numbers that result from integrating over stress tensor components weighted by functions of time (growth factors). Since the time-dependence of the stress tensor is not calculable in a generic $\Lambda$CDM universe, these numbers are free parameters. Crucially though, they can absorb and regularize all contributions from SPT one-loop integrals that are sensitive to a range of scales where our perturbative approach breaks down and have the same $k^{2}$ scaling.

As we have seen in Sec. III.4, the higher-derivative bias terms display the same momentum scaling and thus give rise to very similar terms. Ignoring all precoefficients, evaluation of Eqs. (V.1) and (V.1) shows that the higher-derivative contribution to the power spectrum is exactly degenerate with Eq. (126), while for the bispectrum we get five different terms:

[TABLE]

The contributions iii) and iv) are clearly degenerate with terms in Eq. (127) and using that $\bm{k}_{1}\cdot\bm{k}_{2}=1/2\left(k_{3}^{2}-k_{1}^{2}-k_{2}^{2}\right)$ we see that v) can be written as a combination where $\beta_{B,d}=1/2=-\beta_{B,c}$. Furthermore, one can show that

[TABLE]

such that a combination of all four terms in Eq. (127) with $\beta_{B,a}=1$, $\beta_{B,b}=-1$, $\beta_{B,c}=5/7$ and $\beta_{B,d}=2/7$ can also accommodate for i). Only ii) cannot be expressed through the stress tensor terms and must consequently enter the galaxy bispectrum as an independent contribution. In total, this demonstrates that the stress tensor corrections to the bispectrum are also completely degenerate with those from higher-derivative bias. For that reason we consider them collectively, using the following basis:

[TABLE]

which reduces the number of free parameters to five, in addition to another one for the power spectrum. Note that if we had considered only higher-derivative bias, then we would have only five parameters in total for the power spectrum and bispectrum, since in that case $2\beta_{P}=\beta_{B,a}+\beta_{B,b}$. However that is not the case for stress tensor contributions since $\beta_{P}$ and $\beta_{B}$’s result from integrating the stress tensor weighted with different powers of the growth factor [94]. In addition, it’s worth noting that some of the contributions of velocity bias are also degenerate with higher-derivative bias (see [96]). Finally, note that the contributions in Eq. (135) are degenerate with one-loop galaxy bias contributions in the low-$k$ limit, see Appendix A.3 for an explicit discussion of this.

V.3 Stochastic contributions

The effects discussed in the last section become relevant towards smaller scales, but the bias loop corrections can also have an impact on the large-scale power spectrum and bispectrum. As already pointed out in [97, 7], that is because a subset of the bias loop integrals do not vanish in the large-scale limit, and thus come to dominate for values of $k^{-1}$ above a certain scale. Taking the limit $k\rightarrow 0$ of Eq. (V.1) shows that the galaxy power spectrum approaches a constant:

[TABLE]

where we have assumed that the linear power spectrum falls off to zero on large scales and used that $K(-\bm{q},\bm{q})=0$. Crucially, corrections from successively higher orders of SPT add in comparable measures to this large-scale limit, meaning that this low-$k$ limit is not controlled by any particular order in perturbation theory and thus one must introduce a new parameter that describes its size.

The same extends to the loop corrections of the galaxy bispectrum, and in the limit of all three triangle sides approaching zero, Eq. (V.1) equally becomes constant:

[TABLE]

As for the power spectrum the value of this constant receives corrections from higher-order SPT terms, which are shown schematically up to the three-loop level by the diagrams in Fig. 6. Although diagrams which contain a one-point propagator at one of their external legs do not contribute to this limit, they are affected by large-scale loop corrections in a similar manner. The corresponding terms can be identified as follows:

[TABLE]

where B_{g}(k_{1},k_{2},k_{3})\big{|}_{\propto P_{L}(k_{1})} denotes all terms in the galaxy bispectrum that are proportional to $P_{L}(k_{1})$. After cyclic permutations their total contribution to the one-loop bispectrum is thus given by

[TABLE]

These effects can be interpreted as the impact of small-scale perturbations on the formation of galaxies that cannot be captured by any perturbative bias model, and were originally introduced as a “stochastic” bias in the literature [38, 39]. Its defining property is that it must be mostly uncorrelated on large scales (assuming Gaussian initial conditions as we do), meaning it will manifest in the same way as shot noise. Inspecting Eq. (136) and Eqs. (137, 139), we notice that this is indeed the case for the large-scale limit of the one-loop power spectrum and bispectrum. We could have chosen to incorporate these effects by including stochastic fields in the expansion for $\delta_{g}$ (see for instance [98, 99, 6]), instead we now introduce a posteriori three effective shot noise parameters [7, 8] for the power spectrum and bispectrum, such that

[TABLE]

Similar to the stress tensor parameters, these parameters are able to absorb any residual large-scale contributions stemming from the bias loops. We stress that the values of $N_{0}$, $\epsilon_{0}$ and $\eta_{0}$ are typically not given by their Poissonian shot noise predictions, i.e. $N_{0}=1/\overline{n}$, $\epsilon_{0}=1/\overline{n}^{2}$ and $\eta_{0}=1/\overline{n}$ [100] for an average number density of galaxies $\overline{n}$, but must be determined from the data itself.

Finally, in order to take into account a slight correlation of the stochastic bias on large scales, we can think of Taylor expanding its contributions in powers of $k^{2}$, where in general the expansion coefficients must also be considered as free parameters. This is motivated by requiring that not only the leading order terms in the large-scale limits above can be absorbed, but also their next-to-leading order (NLO) contributions. These terms exclusively scale as powers of $k^{2}$ (or $k_{1}^{2}$, $k_{2}^{2}$ and $k_{3}^{2}$ in case of the bispectrum), and it is straightforward to show that the allowed terms can be summarized as follows:

[TABLE]

which leads to an additional four parameters.

VI Discussion and conclusions

In this paper, we presented the complete set of one-loop (next-to-leading order) galaxy bias corrections for the real-space galaxy bispectrum. These corrections serve to increase the range of scales where the bispectrum can be robustly used to extract information from the clustering of galaxies (or any other tracer), and brings the bispectrum to the same state-of-the-art as the galaxy power spectrum, meaning joint analyses can be performed consistently. We carry this out in detail and compare against numerical simulations of biased tracers in a follow-up paper [76]

Our perturbative bias model systematically combines a variety of effects that become relevant on scales of the weakly non-linear regime. These effects include contributions from: 1) the general bias expansion, generated by gravitational evolution of the matter field; 2) higher-derivative terms, due to the spatial nonlocality of galaxy formation, and a non-zero stress tensor; and 3) stochasticity, resulting from the impact of non-perturbative physics on the large-scale galaxy distribution.

In all generality, we found that a joint one-loop power spectrum and bispectrum analysis in real-space requires 24 free parameters. In agreement with [12], we showed that 11 parameters are due to the general bias expansion, which are shared among the power spectrum and bispectrum. They include the usual linear bias parameter at first order ($b_{1}$), the local and nonlocal quadratic parameters at second ($b_{2},\gamma_{2}$) etc. To account for higher-derivative effects, which were not considered by [12], we demonstrated that we need to introduce an additional five parameters for the bispectrum and one for the power spectrum, five less than given in [6]. We also found that the stress-tensor corrections are entirely degenerate with those from higher derivatives and therefore do not consider them separately. Finally, at leading order stochasticity contributes with two terms (one term) to the bispectrum (power spectrum), whereas at next-to-leading order we identified a further three (one). An overview of these numbers is given in Table 2.

Obviously, it would be bad news if a joint analysis of power spectrum and bispectrum demands that each of the 24 parameters enters the model with a freely adjustable amplitude. In that case the gain from including smaller scales in the analysis could be easily cancelled out by the loss in constraining power from having to marginalize over so many nuisance parameters. From that point of view the work presented in this paper is best regarded as the theory on which we can analyze various ways of reducing the parameter space, as will be discussed in detail in a follow-up paper [76]. For example, it is known in the case of the one-loop galaxy power spectrum that higher-derivatives are degenerate with bias loops [8, 24, 101, 102], and as now shown here (see Appendix A.3) the structure of one-loop bias at low-$k$ is also degenerate with higher-derivative bias for the bispectrum. However, it remains to be seen how the higher-derivative terms compare to those from the general bias expansion. While [30, 31] suggest that the former dominate for tracers residing in very massive halos, this is not necessarily correct for certain types of tracers that are targeted by upcoming galaxy surveys such as DESI. This includes, for instance, emission line galaxies, which are usually found in less massive halos.

Apart from testing the significance of the higher-derivative contributions, another clear possibility of reducing the parameter space is to check which bias operators determine the galaxy perturbations at the initial time. Our bias basis is designed so that we separate operators which are local in second derivatives of the linear potential (what we call local evolution operators, or LE), and those which are nonlocal (induced by nonlocal evolution, hence NLE operators) that involve the nonlinear Lagrangian potentials. Other bases in the literature mix these properties, see Appendix A.2.

The most drastic reduction in parameter space is achieved by putting initially all NLE operators to zero and all LE operators nonlocal in $\delta$ to zero, this is the local Lagrangian approximation: it reduces the 11 parameters in Table 2 to just 3 free parameters ($b_{1},b_{2},b_{3}$) and the other eight parameters get determined from them by the time evolution arguments in Section IV.3.1 (putting all of the eight Lagrangian values to zero). A hierarchy of more accurate approximations (but with more free parameters) can be made where one takes e.g. only the NLE operators to be zero in the far past (as in peak or excursion set models of bias), or just the 4th-order ones that only enter through the loop correction of a single diagram in the bispectrum. The validity of such assumptions must be tested numerically and compared to the full model, which we will present in [76]. For practical applications to galaxy surveys, given finite error bars, one should be able to reduce the parameter space sequentially going from most general to least general checking the final results are not impacted by the assumptions.

On the more technical side, the central concept explored in this work is the use of multipoint propagators, which correspond to the sum over all reducible diagrams, as a way to writing the perturbative expansion in terms of explicit observables (that corresponds for Gaussian fluctuations to the commonly measured cross-correlation bias between galaxies and matter). As a result of this reformulation we were able to circumvent the renormalization procedure — a tedious redefinition of the standard bias parameters required to absorb diverging contributions in the computation of correlation functions at one-loop and beyond. To achieve this we showed that it is best to proceed in two steps: first, compute the multipoint propagators at initial time where nonlinearities in the matter field can be ignored, then evolve these initial propagators to the time of observation. Evolving the propagators conserving the number of tracers (the so-called “coevolution” first given by [69]) cannot generate any diverging terms proportional to $\sigma^{2}=\langle\delta(\bm{x})^{2}\rangle$ and therefore the late-time propagators are already renormalized (and any $k^{2}$ or higher-order renormalizations are simple to handle by adding the corresponding terms). Once the evolved multipoint propagators have been computed, they serve as simple building blocks for the general $N$-point correlation function.

We also discussed how our choice of basis, written in terms of Galileons and with a clear split between local and nonlocal functions of second derivatives of the linear potential, simplifies the calculation of the multipoint propagators: they are given by the tree-level expressions plus loop integrals over the NLE operators contributions only. Other bases that mix our operators (see Appendix A.2) complicate the calculation, as one has to separate terms that get absorbed by the renormalization procedure.

A significant part of the evolution of the propagators is absorbed into the evolution of the bias parameters in them. Our calculation of this time evolution was carried out using a new approach that takes full advantage of Galilean invariance (see Appendix C) and we present results for the fourth-order parameters for the first time (see Section IV.3.1), generalizing the third-order results originally given in [9]. Up to quadratic bias we also show explicitly that renormalization by fourth-order operators does not change the time evolution of bias parameters obtained at tree level (see Appendix D), as expected.

Finally, to summarize the steps necessary for applying the results of this paper in practice, one would proceed as follows: first, one computes the Lagrangian multipoint propagators given in Eqs. (65), (IV.2.2) and (IV.2.3). These serve as initial conditions for the evolution equations yielding the Eulerian multipoint propagators, i.e. Eqs. (209-B.3), which in turn determine the power spectrum and bispectrum at the time of observation via Eqs. (V.1) and (V.1). The final expressions are obtained after adding the combination of stress-tensor and higher-derivative corrections from Eqs. (126) and (135), as well as the stochastic contributions up to next-to-leading order (Eqs. 140-143).

This will be demonstrated explicitly in our follow-up paper [76]. In particular, we carry out a detailed likelihood analysis of the power spectrum and bispectrum in numerical simulations of biased tracers to ascertain the performance of the results derived here, the importance of the different terms in the bias expansion, the accuracy of theoretical assumptions that reduce the number of free parameters, and the improvement of our results compared to current models in the literature.

Acknowledgements.

We thank M. Crocce, A. Sanchez, R. Sheth for useful discussions. AE acknowledges support from the UK Science and Technology Facilities Council (STFC) via Research Training Grant (grant number ST/M503836/1), as well as from the European Research Council (grant number ERC-StG-716532-PUNCA). AE also thanks the Center for Cosmology and Particle Physics of New York University for hospitality while part of this work was undertaken, and the STFC for supporting this research stay. RES acknowledges support from the STFC (grant number ST/P000525/1).

Appendix A Further notes on Galilean basis for galaxy bias

A.1 Basis operators in Fourier space

Here we briefly summarize Fourier space expressions for our basis operators given in Table1, which are being used in the computation of the multipoint propagators. In genereal, we write any $n$-th order operator ${\cal O}^{(n)}$ as an integral over $n$ linear matter perturbations:

[TABLE]

where $B$ stands for any of the basis operators at that order and we use the notation $\left[\delta_{D}\right]_{n}\equiv\delta_{D}(\bm{k}-\bm{k}_{1\ldots n})$. For $n=2$ we have from Eq. (32):

[TABLE]

and the kernel $K(\bm{k}_{1},\bm{k}_{2})$ was already defined in Eq. (34). At third order the only nontrivial operator is ${\cal G}_{2}(\varphi_{2},\varphi_{1})$, which becomes upon Fourier transformation:

[TABLE]

where we have made use of Eq. (32) and made the redefinition $\bm{k}_{2}\,\rightarrow\,\bm{k}_{23}$ in the second step. After symmetrization we then obtain:

[TABLE]

Next, let us consider the most complicated combination that appears at fourth order, ${\cal G}_{2}(\varphi_{3},\varphi_{1})$, all other operators will follow in a very similar manner. Starting from the definition in Eq. (III.3) and using the relations (47) to (49) for the LPT potentials, we have

[TABLE]

Plugging in the Fourier expressions for the remaining potentials and Galileons (using Eq. A.1), and replacing $\bm{k}_{2}\,\rightarrow\,\bm{k}_{234}$ we get

[TABLE]

where we have introduced the new kernel

[TABLE]

which is symmetric under cyclic permutations of its four momenta. The fully symmetric kernels for the basis operators at fourth order are thus given by

[TABLE]

where $\text{sym.}(n)$ denotes the total number of terms the expressions have to be symmetrized over with respect to the four wave vectors.

A.2 Relation to other bias bases in the literature

In [6], the authors extend the bias basis from [13]. Let us establish the connection between their basis of operators and ours up to fourth order. They use operators that correspond to tracing objects denoted as $\Pi^{[n]}$ defined from convective derivatives. The calculation of the connection between both bases is conceptually simple but in practice long and tedious beyond second order, so we only provide the final results. To linear order, we simply have ${\rm Tr}[\Pi^{[1]}]=\delta$. The only non-trivial (i.e. beyond $\delta^{2}$) at quadratic order is ${\rm Tr}[(\Pi^{[1]})^{2}]=\delta^{2}+{\cal G}_{2}$. At cubic order we have for the non-trivial ones,

[TABLE]

whereas at fourth order we have,

[TABLE]

Note that these relations show that this basis mixes operators that are local and nonlocal in second derivatives of the linear potential.

Another set of papers use a yet different basis [70, 103, 30, 31]. The comparison to these cases is more complicated since they also include velocity bias operators mixed together with standard ones. The simplest to compare with is [31], which when constrained to the no-velocity bias case shows one more operator (five in total) at third order than in our basis. In fact, this set is precisely the same set used in [8]; however, as pointed out in [9] not all of these five operators are independent (only four of them are). In terms of our basis, apart from $\delta^{3}$ they have at cubic order the following operators

[TABLE]

which means that we can express $\delta s^{2},s^{3},st,\psi$ each in terms of the four other basis elements in their basis, e.g.

[TABLE]

Again, it is worth noting that Eq. (172) shows that this basis mixes operators that are local and nonlocal in second derivatives of the linear potential.

The papers [103, 30, 31] do not consider fourth-order operators, but they do include higher-derivative operators. For example, in the bispectrum [31] include terms corresponding to our $\beta_{1},\beta_{2,1},\beta_{2,2}$ (but they don’t include $\beta_{2,3},\beta_{2,4}$) and for the noise terms they include $\eta_{2,1},\eta_{2,2}$ and seem to be missing $\epsilon_{2}$ (though they do include such terms in the power spectrum, which for reasons that are unclear also includes a $k^{4}$ noise term); in addition they have the usual $\epsilon_{0},\eta_{0}$ shot noise terms. See our Table 1 and Eqs. (141) and (143) for reference of the operators associated with these coefficients.

Finally, the approach described in [47, 48] uses a number of basis of operators obtained from different phenomenological models of biasing (e.g. halos, peaks). Note that the ’renormalized bias functions’ they define correspond to our initial multipoint propagators, whereas what they refer to as the multipoint propagators correspond to our evolved multipoint propagators. The LPT formalism of [46, 47] was later extended by [104] and in the context of the two-point correlation function the authors of [105] introduced second order bias operators equivalent to those given in [8], as well as a higher-derivative term. However, they neglect third order bias parameters and thus do not cover the complete bias model at one-loop order.

A.3 Large-scale limit of 4th-order NLE operators

As shown in Sec. IV.2, loop corrections of the multipoint propagators are exclusively given in terms of the NLE operators. In the large-scale limit these loop integrals can be Taylor expanded in powers of the participating wave vectors. For the 4th-order contributions to the two-point propagator we thus obtain:

[TABLE]

where ${\cal O}(k_{i}^{4})$ denotes terms which are of order $k_{1}^{4}$, $k_{2}^{4}$ or $k_{12}^{4}$ and higher. These expressions serve to illustrate that the nonlocal evolution terms in the large-scale limit are given by combinations of the higher-derivative operators. The converse is equally true as can be easily verified by inverting the relations above.

Appendix B Time evolution of multipoint propagators

In this appendix we give a detailed derivation for the time evolution of the multipoint propagators. We proceed in two steps: first, we evaluate expectation values of products of two or three Wiener-Hermite functionals, and second, by using these results we directly show how to obtain the recursion relations reported in Eq. (IV.3.2).

B.1 Orthogonality relations for generalized Wiener-Hermite functionals

Let us consider the PDF of $\delta_{L}$, shifted by a generic source term $\alpha(\bm{k})$, which we take to be an arbitrary function of wavenumber $\bm{k}$. A Taylor expansion around $\alpha=0$ yields:

[TABLE]

where $\alpha_{i}\equiv\alpha(\bm{k}_{i})$. Swapping the derivatives from $\alpha$ to $\delta_{L}$ and using the definition of the Wiener-Hermite functionals from Eq. (55) we get

[TABLE]

and thus:

[TABLE]

To derive the orthogonality relation between two Wiener-Hermite functionals of orders $m$ and $n$, we first compute the following integral

[TABLE]

where we have plugged in Eq. (54). However, using Eq. (180) to replace the PDF’s, we must also have:

[TABLE]

and by Taylor expanding the righthand side of the expression above, we see that we need to match up all $\bm{k}$ and $\bm{q}$ modes, which is only possible if $m=n$. From that observation it immediately follows that

[TABLE]

where we have used the short-hand notation,

[TABLE]

for writing all possible ways ($m!$ in total) of matching up the two sets of modes.

Let us now compute the expectation value of three Wiener-Hermite functionals, i.e. $\langle{\cal H}_{m}\,{\cal H}_{n}\,{\cal H}_{l}\rangle$. Similar to the above procedure, we first evaluate an integral over PDF’s, now with the three different sources $\alpha$, $\beta$ and $\gamma$:

[TABLE]

Expanding both sides of Eq. (185), we get

[TABLE]

and

[TABLE]

In order for Eqs. (B.1) and (187) to be equal, they need to contain the same number of source terms, which requires that $l+m+n=2(a+b+c)$ and thus, $l+m+n\in 2\mathbb{N}$. Moreover, the indices must satisfy the conditions

[TABLE]

from which follows that $m+n\geq l$, and cyclic permutations thereof. According to Eq. (194), we can divide all $\bm{k}$-, $\bm{q}$-, and $\bm{p}$-modes into two subsets each, either of size $a$, $b$ or $c$. For instance, the $\bm{k}$-modes will be split into a group containing $a$ modes, and another containing $b$, such that $a+b=m$. Each mode in the former group can then be assigned a mode from an equally sized group of $\bm{q}$’s, i.e. $\bm{q}_{1}=-\bm{k}_{1},\,\ldots,\,\bm{q}_{a}=-\bm{k}_{a}\,,$ while the latter are matched in a similar manner with a set of $\bm{p}$-modes, which guarantees that we obtain Eq. (187). Altogether, this means we need to require that the expectation value of the three Wiener-Hermite functionals is given by

[TABLE]

and it must vanish if the conditions above are not satisfied. The number of terms that have to be added in order to symmetrize Eq. (B.1) is given by the number of possibilities of selecting subsets of $\bm{k}$-, $\bm{q}$-, and $\bm{p}$-modes that are of size $b$, $a$ and $c$, respectively. The total number of terms in the square brackets is thus:

[TABLE]

B.2 The $\Gamma$-recursion relation

We are interested in the time evolution of the $n$-th order multipoint propagator. Exploiting the orthogonality of Wiener-Hermite functionals, we can single out its contribution to the series expansion in Eq. (IV.3.2) by multiplying both sides with ${\cal H}_{n}$ and taking the ensemble average:

[TABLE]

Applying the same procedure to the integral solution for $\Psi_{a}$ from Eq. (IV.3.2) we get

[TABLE]

where $\Gamma_{a,{\cal L}}$ denotes multipoint propagators at initial time $\eta=0$. Upon inserting Eq. (B.1) we perform the first $a=(n+m-l)/2$ integrations over $\bm{q}$-modes as well as over all of the $\bm{p}$-modes, such that the second term above reduces to

[TABLE]

The symmetrization is carried out by summing over all $\binom{n}{a}$ subsets of $\bm{k}$-modes, which implies that the combinatorial factor $\kappa$ is given by:

[TABLE]

Next, let us change the summation index from $l$ to $r=(l+m-n)/2=m-a$ and relabel the remaining integrations over $\bm{q}$’s into $\bm{q}_{1},\ldots,\bm{q}_{r}$. The conditions $n+m\geq l$ (and cyclic permutations) then transform into $r_{*}\leq r\leq m$, where

[TABLE]

Defining the quantity

[TABLE]

we finally obtain the desired recursion relation for multipoint propagators:

[TABLE]

where we have arranged the two summations in the second contribution such that the outer sum runs over the number of loop integrals, $r$, and used Eq. (204) to determine the allowed range of the index $m$.

B.3 Computation of time-evolved propagators up to one-loop order

We now use the recursion relation from Sec. B.2 to derive explicit expressions for the building blocks of the one-loop bispectrum, which requires us to evaluate the first two multipoint propagators at one-loop order and the three-point propagator at tree-level. In order to make the connection to the initial conditions evident, where convenient, we write them completely in terms of the Lagrangian multipoint propagators, which we computed in Sec. IV.2.

To simplify the resulting expressions it is useful to define the following two quantities that arise through contractions of the linear propagator $g_{ab}$ with the vertices $\gamma_{abc}$:

[TABLE]

where we let the superscript indicate the number of external lines that are attached to these objects, i.e. two for $g^{(2)}_{abc}$ etc. The tree-level expression ($r=0$) for $\Gamma_{a}^{(1)}$ is trivial since the second piece in Eq. (206) makes no contribution. For that reason we start with the two-point propagator, and by using the kernel from Eq. (207) we obtain

[TABLE]

with all initial multipoint propagators being evaluated at tree-level. The tree-level three-point propagator receives two contributions from the second term in Eq. (206), such that

[TABLE]

where we have plugged in the solution for $\Gamma^{(2)}_{a}$ from above. Similarly, Eq. (209) can be used to determine the nonlinear correction to the one-point propagator at one-loop order (or $r=1$), giving

[TABLE]

Equations (209-211) are the full results including all transients and the evolution of the Lagrangian bias parameters. For comparison to observational data only the fastest growing mode is relevant and we can further absorb all terms due to bias evolution into the final Eulerian bias parameters (see Sec. IV.3.1), allowing us to greatly simplify these expressions. Moreover, if we separate the contributions that get newly generated by nonlinear evolution (i.e. they are not due to bias evolution of terms already present in the initial multipoint propagators), and expressing them in terms of the SPT kernels, we arrive at Eqs. (IV.3.2-120).

While Eqs. (IV.3.2-120) can be more easily obtained from the nonlinear evolution of the matter and velocity fields directly, starting from $\Gamma^{(2)}$ at one-loop order the evolution of the multipoint propagators becomes useful, since we do not have to deal with cumbersome renormalizations as shown in Sec. II.3. Because the result is still complicated we have chosen to present it in terms of the previously derived expressions for the evolved propagators,

[TABLE]

where we have indicated at which order the propagators are evaluated.

Appendix C Using Galilean Invariance to compute time evolution of bias

C.1 SPT in Galilean basis for galaxy bias

Let us consider the Galilean invariant contributions to SPT at each order in perturbation theory (PT). This will allow us to calculate in a more straightforward way the time evolution of bias parameters up to any order in PT (we will compute up to fourth-order here). Since our bias basis contains the LPT potentials, we start from the expression for density fluctuations in terms of the Lagrangian displacement field:

[TABLE]

where $\Psi_{i,j}\equiv\partial\Psi_{i}/\partial q_{j}$ ($\delta_{ij}^{\rm K}$ denotes the Kronecker delta) and the $G_{i}$’s are the Galileons built out of the deformation tensor, working up to fourth-order in PT they read:

[TABLE]

By using these results in Eq. (213), and using that $\bm{x}=\bm{q}+\bm{\Psi}$, one can derive the standard Eulerian PT expansion. However, we are interested in finding out the bias relation in the absence of velocity bias, therefore any dipole terms coming from displacement from $\bm{q}$ to $\bm{x}$ will not appear in any of the final expressions, and our basis being Galilean invariant does not include any such operators. Therefore, we can simply work with the Galilean invariant SPT expression at each order, since dipole terms will cancel at the end of the calculation. Thus we have,

[TABLE]

Using Eqs. (214-216) in Eq. (217) we obtain the Galilean invariant piece of the SPT expansion at each order in our basis of bias operators,

[TABLE]

[TABLE]

A simple check on these expressions is provided by the spherical collapse dynamics, which is also Galilean invariant. To do this all we need is the spherical average of our basis operators, which is straightforward

[TABLE]

Using these results in Eqs. (218-219) we recover the well-known vertices $\nu_{n}$ from the spherical collapse dynamics [106], i.e.

[TABLE]

Now we proceed to derive the analogous results for the velocity field, which is what is needed for the nonlinear evolution of $\Phi_{v}$ or the normalized divergence $\theta=\nabla^{2}\Phi_{v}$. By definition, the Eulerian velocity field agrees with the velocity of particles calculated from the time derivative of the displacement field at the same location. Taking the normalized divergence $-\nabla/{\cal H}f$ we then have that

[TABLE]

where ${}^{\prime}\equiv\partial/\partial{\ln D}$, with $D$ the linear growth factor. Again, we are only interested in the Galilean invariant piece of this exact relation. In order to obtain that, all we need is $\bm{\Psi}^{\prime}(\bm{q})$ as a function of $\bm{x}$, with at most one displacement free from derivatives acting on it, because other contributions cannot be rendered Galilean invariant by taking a single derivative when calculating $\theta_{\text{GI}}(\bm{x})$. A straightforward Taylor expansion gives then,

[TABLE]

where all quantities on the RHS of this equation are now evaluated at $\bm{x}$. The result for $\theta_{\text{GI}}$ follows by acting with the divergence in Eq. (222) on the derivative-free displacement field in each term of Eq. (223) giving us the simple result,

[TABLE]

The traces here can be easily rewritten in terms of the deformation tensor Galileons in Eq. (214), and the derivative with respect to $\ln D$ just multiplies each contribution by its PT order. Therefore we obtain the Galilean invariant piece of the normalized divergence at each order,

[TABLE]

[TABLE]

Again, we can double check these results comparing them to the spherical collapse dynamics, by performing the spherical average using Eq. (C.1) in Eqs. (225-226) to compute the vertices $\mu_{n}$ for the normalized velocity divergence,

[TABLE]

which agree with the standard spherical collapse dynamics values [106].

C.2 Time Evolution of Bias

As already explained in Sec. IV.3.1, in order to determine the time evolution of bias, we can use a simple trick that exploits the Galilean invariance of the galaxy overdensity in the absence of velocity bias. That means all dipole terms in the matter density must cancel, such that

[TABLE]

where $\delta^{\text{GI}}$ is the Galilean invariant part of the (nonlinear) matter density as determined in the last section, and $\delta_{g}^{\text{IC}}$ is the galaxy density at initial time (where all fluctuations are linear and Gaussian). At linear order it is easy to check that Eq. (228) simply gives $b_{1}=1+b_{1,{\cal L}}$, which agrees with the full evolution [69] when neglecting transients. To find the evolution of higher-order bias parameters, we must first subtract all contributions due to lower order operators (see [9]). For instance, by computing (suppressing $\bm{x}$ arguments)

[TABLE]

and writing the RHS of this equation in our bias basis using Eqs. (77) and (218) at second order, we are able to identify the coefficients of the $\delta^{2}$ and ${\cal G}_{2}$ operators, which leads to

[TABLE]

At cubic order we proceed in the same way and compute,

[TABLE]

and we use Eqs. (77) and (218) plus that

[TABLE]

which follows from Eq. (93) in [9] when constrained to Galilean invariant operators (for which we need Eqs. 218 and 225). We thus identify the coefficients (i.e. bias parameters) in front of the operators in our bias basis at third order,

[TABLE]

Equations (230)-(233) agree with those derived in [9] up to decaying modes (which we neglect here) and a typo in their Eq. (116) which should have a minus sign for its second term.

Finally, at quartic order we need to subtract all fourth-order contributions from up to cubic bias,

[TABLE]

and then, using that (again, constrained to Galilean invariant operators)

[TABLE]

[TABLE]

we can identify the evolved bias parameters from the coefficients of each of the fourth-order operators in our bias basis. We thus obtain,

[TABLE]

C.3 Spherical Average Limit

Let us now briefly consider the spherical-average limit of these results, which is somewhat subtle due to the fact that operators beyond local bias do not vanish under spherical average. As a result of this, the spherically averaged bias parameters $b_{n}^{\text{sph}}$ are linear combinations of all the bias parameters (local and nonlocal) in our basis.

Consider quadratic bias, where the impact of spherical averaging first appears (linear bias is invariant under spherical average). From Eq. (C.1) one can take the spherical average of the quadratic bias relation $(b_{2}/2)\,\delta^{2}+\gamma_{2}\,{\cal G}_{2}$ and write the spherically averaged quadratic bias parameter in terms of our bias parameters as,

[TABLE]

where in the second equality we have used our time-evolution result from Eq. (230) and then rewritten it in terms of the Lagrangian spherical averaged quadratic bias parameter, giving the well-known spherical collapse model bias evolution at quadratic order [73].

At cubic order, we see for the first time something subtler. The reason is that the spherical average of Eq. (231) does not correspond to

[TABLE]

due to the fact that at third-order the spherically averaged Galileon operator obeys from Eqs. (C.1) and (232)

[TABLE]

Therefore, at cubic order in our bias basis part of what would be considered $b_{3}^{\text{sph}}$ is encoded by the deviation of the spherical average of $[{\cal G}_{2}]^{(3)}$ from that of $(-2/3)[\delta^{2}]^{(3)}$, the precise relation is then

[TABLE]

Note this apparent dependence of $b_{3}^{\text{sph}}$ on a quadratic bias parameter such as $\gamma_{2}$ only happens at the final conditions, since at the initial conditions matter fluctuations are linear and thus $[{\cal G}_{2}]^{(3)}=[\delta^{2}]^{(3)}=0$. Therefore $\gamma_{2,{\cal L}}$ cannot appear in $b_{3,{\cal L}}^{\text{sph}}$, and in fact all this means that the only role of the $\gamma_{2}$ term in Eq. (243) is to cancel the appearance of such terms in the evolved values of $\gamma_{3}$ and $\gamma_{21}$ in Eq. (233) and in $b_{3}$ and $\gamma_{2}^{\times}$ through $b_{2}=b_{2}^{\text{sph}}+4\gamma_{2}/3$. Indeed, evaluating Eq. (243) yields,

[TABLE]

where $b_{3,{\cal L}}^{\text{sph}}=b_{3,{\cal L}}-4\gamma_{2,{\cal L}}^{\times}+(4/3)\gamma_{3,{\cal L}}+(8/3)\gamma_{21,{\cal L}}$. At quartic order, a similar subtlety arises with all nonlocal bias parameters of lower order than fourth. However, for the same reasons as in the cubic bias calculation, all such terms must cancel at the end. Therefore, a fast way of doing the calculation is to simply ignore all nonlocal quadratic and cubic bias parameters contributions to Eqs. (237-239) and in the relation between spherical bias parameters and local bias parameters. A quick calculation then gives,

[TABLE]

where $b_{4,{\cal L}}^{\text{sph}}=b_{4,{\cal L}}-16\gamma_{2,{\cal L}}^{\times\times}+(16/3)\gamma_{3,{\cal L}}^{\times}+(32/3)\gamma_{2,{\cal L}}^{\text{sq}}+(32/3)\gamma_{21,{\cal L}}^{\times}-(32/9)\gamma_{211,{\cal L}}-(64/9)\gamma_{22,{\cal L}}-(368/567)\gamma_{31,{\cal L}}$. The results in Eqs. (244) and (245) agree with those in [46] who corrected the cubic and quartic spherical collapse bias results in [73].

Appendix D Evolution of Renormalized Bias

In Appendix C we discussed time evolution of bias parameters in terms of the standard expansion, as opposed to the bias parameters that appear in the multipoint-propagator expansion, therefore we ignored renormalization. One may ask whether those results are somehow changed due to renormalization. In fact, [12] contrasted the relation between renormalized and bare (final time) bias parameters to the time evolution of bare bias parameters found in [9], and appeared to suggest that the results of [9] may be affected by renormalization.

What matters of course is the time evolution of the observable bias parameters, i.e. how the renormalized late-time bias parameters are related to the renormalized initial-time bias parameters (this is what is often measured in numerical simulations). We therefore start from the initial conditions written using the multipoint propagator expansion up to fourth order,

[TABLE]

where, unlike Eq. (77) the bias parameters here are renormalized, and we have neglected ${\cal O}(\sigma^{4})$ contributions. Also, we have used that the loops in the multipoint propagators coming from nonlocal evolution operators cancel against the corresponding terms in the Gamma expansion to make the expression more compact.

The next step is to use Eq. (228) to evolve this to find the late-time evolved $\delta_{g}$, and find the evolved renormalized bias parameters. To compare with and address the issue suggested by [12] we proceed as they did and write the non-linear renormalization formulae for linear and quadratic bias parameters to ${\cal O}(\sigma^{2})$, following the calculations we discussed in Section II.3, that is

[TABLE]

for linear bias [8] and

[TABLE]

for the quadratic bias parameters. Note that all bias parameters in Eqs. (247-249) are final-time quantities. The first parenthesis contributions in each of these expressions correspond to the renormalizations by operators that depend on the linear Gaussian fluctuations (that is, those that are straightforward to compute) and the second parenthesis contributions describe renormalizations induced by non-linear evolution (discussed in Section II.3) that were bypassed in the method we advocate elsewhere in the main text. Apart from some sign typos and the renormalizations induced by $\bar{\gamma}_{2}^{\times}$ (what they call $b_{{\cal G}_{2}\delta}^{(0)}$), Eqs. (248-249) agree with Eqs. (3.17-18) in [12].

Using Eq. (246) in Eq. (228), we can identify the evolved bare bias parameters as the coefficients of the operators in our bias basis up to quadratic order, and then renormalize according to Eqs. (247)-(249), finding the evolved renormalized bias parameters in terms of the initial renormalized bias parameters. In doing so we have taken into account all renormalizations induced by operators up to fourth order (both in the initial bias parameters and in the final ones). As expected, all $\sigma^{2}$ dependencies go away (a nontrivial check on the nonlinear renormalization calculations in Eqs. 248-249), and we recover

[TABLE]

i.e. the same results obtained from the evolution of the bare bias parameters, showing explicitly that the bare evolution results in [9] are inmune to renormalization. Therefore final-time renormalization formulae such as Eqs. (247-249) are not in conflict with time evolution of bare bias parameters, but rather guarantee that the time evolution results also hold for the renormalized bias parameters. In our approach advocated in the main text, the initial-time renormalization (equivalent to the first parenthesis in Eqs. 247-249) is automatically carried out by the multipoint propagator expansion, while the ‘second parenthesis’ renormalizations are never needed because time evolution of conserved tracers cannot generate such divergences (i.e. the evolved propagators already avoid spurious contributions proportional to $\sigma^{2}$).