Velocity bias and the nonlinear perturbation theory of peaks

TL;DR

This paper explores the origin of velocity bias in the peaks model of large-scale structure biasing using integrated perturbation theory, revealing how the flat constraint leads to velocity bias and deriving nonlinear velocity dispersion formulas.

Contribution

It provides a systematic derivation of velocity bias within the peaks model using nonlinear perturbation theory and develops a formalism to trace the nonlinear evolution of biased peaks.

Findings

Velocity bias arises from the flat constraint at peaks.

A formal perturbation theory for nonlinear evolution of biased peaks is developed.

A one-loop formula for nonlinear velocity dispersion of peaks is derived.

Abstract

The biasing in the large-scale structure of the universe is a crucial problem in cosmological applications. The peaks model of biasing predicts a linear velocity bias of halos, which is not present in a simple model of local bias. We investigate the origin of the velocity bias in the peaks model from the viewpoint of the integrated perturbation theory, which is a nonlinear perturbation theory in the presence of general Lagrangian bias. The presence of the velocity bias in the peaks model is a consequence of the "flat constraint," ${\nabla}\delta = 0$; i.e., all the first spatial derivatives should vanish at the locations of peaks. We show that the velocity bias in the peaks model is systematically derived in the framework of the integrated perturbation theory, and then develop a formal theory to perturbatively trace the nonlinear evolution of biased objects with the flat constraint. A formula for the nonlinear velocity dispersion of peaks with the one-loop approximation is also derived.

Equations

$${\mathbf{j}}_X(\mathbf{x})=[1+{\delta }_X(\mathbf{x})]\mathbf{v}(\mathbf{x}),$$

$${\mathbf{j}}_X^{\mathrm{L}}(\mathbf{q})=[1+{\delta }_X^{\mathrm{L}}(\mathbf{q})]{\mathbf{v}}^{\mathrm{L}}(\mathbf{q}),$$

$${\mathbf{j}}_X(\mathbf{x})=\int d^{3}q{\mathbf{j}}_X^{\mathrm{L}}(\mathbf{q}){\delta }_{\mathrm{D}}^3[\mathbf{x}-\mathbf{q}-\Psi (\mathbf{q})].$$

$${\mathbf{v}}^{\mathrm{L}}(\mathbf{q})=a\dot{\Psi }(\mathbf{q}),$$

$${\mathbf{j}}_X^{\mathrm{L}}(\mathbf{k})={\mathbf{v}}^{\mathrm{L}}(\mathbf{k})+{\int }_{{\mathbf{k}}_1+{\mathbf{k}}_2=\mathbf{k}}{\delta }_X^{\mathrm{L}}({\mathbf{k}}_1){\mathbf{v}}^{\mathrm{L}}({\mathbf{k}}_2),$$

$${\int }_{{\mathbf{k}}_1+\cdots +{\mathbf{k}}_2=\mathbf{k}}\cdots \equiv \int \frac{d^3{k}_1}{(2\pi {)}^3}\cdots \frac{d^3{k}_n}{(2\pi {)}^3}{\delta }_{\mathrm{D}}^3({\mathbf{k}}_1+\cdots +{\mathbf{k}}_n-\mathbf{k})\cdots .$$

$$\begin{gathered} ⟨{\delta }_{\mathrm{L}}(\mathbf{k}){\mathbf{j}}_X^{\mathrm{L}}({\mathbf{k}}^{\prime })⟩=⟨{\delta }_{\mathrm{L}}(\mathbf{k}){\mathbf{v}}^{\mathrm{L}}({\mathbf{k}}^{\prime })⟩ \\ +{\int }_{{\mathbf{k}}_1+{\mathbf{k}}_2={\mathbf{k}}^{\prime }}⟨{\delta }_{\mathrm{L}}(\mathbf{k}){\delta }_X^{\mathrm{L}}({\mathbf{k}}_1){\mathbf{v}}^{\mathrm{L}}({\mathbf{k}}_2)⟩. \end{gathered}$$

$$\Psi (\mathbf{k})=\sum _{n=1}^{\infty }\frac{i}{n!}{\int }_{{\mathbf{k}}_1+\cdots +{\mathbf{k}}_n=\mathbf{k}}{\bar{\mathbf{L}}}_n({\mathbf{k}}_1,\dots ,{\mathbf{k}}_n){\delta }_{\mathrm{L}}({\mathbf{k}}_1)\cdots {\delta }_{\mathrm{L}}({\mathbf{k}}_n),$$

$${\bar{\mathbf{L}}}_n({\mathbf{k}}_1,\dots ,{\mathbf{k}}_n)=W(k_{1}R)\cdots W(k_{n}R){\mathbf{L}}_n({\mathbf{k}}_1,\dots ,{\mathbf{k}}_n),$$

$${\bar{\mathbf{L}}}_n({\mathbf{k}}_1,\dots ,{\mathbf{k}}_n)=W(|{\mathbf{k}}_1+\cdots +{\mathbf{k}}_n|R){\mathbf{L}}_n({\mathbf{k}}_1,\dots ,{\mathbf{k}}_n).$$

$${\bar{\mathbf{L}}}_1(\mathbf{k})=W(kR){\mathbf{L}}_1(\mathbf{k}).$$

$${\mathbf{L}}^{(1)}(\mathbf{k})=\frac{\mathbf{k}}{k^2},$$

$${\mathbf{L}}^{(2)}({\mathbf{k}}_1,{\mathbf{k}}_2)=\frac{3}{7}\frac{{\mathbf{k}}_{12}}{{k_{12}}^2}\left[1-{\left(\frac{{\mathbf{k}}_1\cdot {\mathbf{k}}_2}{k_1{k}_2}\right)}^2\right],$$

$${\mathbf{L}}^{(3)}({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3)=\frac{1}{3}\left[{\mathbf{L}}^{(3\mathrm{a})}({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3)+perm.\right];$$

$${\mathbf{L}}^{(3\mathrm{a})}({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3)$$

$$=\frac{{\mathbf{k}}_{123}}{{k_{123}}^2}\{\frac{5}{7}\left[1-{\left(\frac{{\mathbf{k}}_1\cdot {\mathbf{k}}_2}{k_1{k}_2}\right)}^2\right]\left[1-{\left(\frac{{\mathbf{k}}_{12}\cdot {\mathbf{k}}_3}{k_{12}{k}_3}\right)}^2\right]$$

$$-\frac{1}{3}\left[1-3{\left(\frac{{\mathbf{k}}_1\cdot {\mathbf{k}}_2}{k_1{k}_2}\right)}^2+2\frac{({\mathbf{k}}_1\cdot {\mathbf{k}}_2)({\mathbf{k}}_2\cdot {\mathbf{k}}_3)({\mathbf{k}}_3\cdot {\mathbf{k}}_1)}{{k_1}^2{k_2}^2{k_3}^2}\right]\}$$

$$+\frac{3}{7}\frac{{\mathbf{k}}_{123}\times ({\mathbf{k}}_1\times {\mathbf{k}}_{23})}{{k_{123}}^2{k_1}^2{k_{23}}^2}({\mathbf{k}}_1\cdot {\mathbf{k}}_{23})\left[1-\left(\frac{{\mathbf{k}}_2\cdot {\mathbf{k}}_3}{k_2{k}_3}\right)\right].$$

$${\delta }_X^{\mathrm{L}}(\mathbf{k})$$

$${\mathbf{v}}^{\mathrm{L}}(\mathbf{k})$$

$$\times {\int }_{{\mathbf{k}}_1+\cdots +{\mathbf{k}}_n=\mathbf{k}}{\bar{\mathbf{L}}}_n({\mathbf{k}}_1,\dots ,{\mathbf{k}}_n){\delta }_{\mathrm{L}}({\mathbf{k}}_1)\cdots {\delta }_{\mathrm{L}}({\mathbf{k}}_n),$$

$$⟨{\delta }_{\mathrm{L}}({\mathbf{k}}_1){\mathbf{j}}_X^{\mathrm{L}}(\mathbf{k})⟩=⟨{\delta }_{\mathrm{L}}({\mathbf{k}}_1){\mathbf{v}}^{\mathrm{L}(1)}(\mathbf{k})⟩+⟨{\delta }_{\mathrm{L}}({\mathbf{k}}_1){\mathbf{v}}^{\mathrm{L}(3)}(\mathbf{k})⟩+\cdots$$

$$+{\int }_{{\mathbf{k}}^{\prime }+{\mathbf{k}}^{\prime \prime }=\mathbf{k}}[⟨{\delta }_{\mathrm{L}}({\mathbf{k}}_1){\delta }_X^{\mathrm{L}(2)}({\mathbf{k}}^{\prime }){\mathbf{v}}^{\mathrm{L}(1)}({\mathbf{k}}^{\prime \prime })⟩$$

$$+⟨{\delta }_{\mathrm{L}}({\mathbf{k}}_1){\delta }_X^{\mathrm{L}(1)}({\mathbf{k}}^{\prime }){\mathbf{v}}^{\mathrm{L}(2)}({\mathbf{k}}^{\prime \prime })⟩+\cdots ]$$

$$=iaHf(2\pi {)}^3{\delta }_{\mathrm{D}}^3({\mathbf{k}}_1+\mathbf{k})P_{\mathrm{L}}(k_1)$$

$$\displaystyle\hskip 36.0pt\times\Biggl{[}\bar{\bm{L}}_{1}(\bm{k})+\frac{3}{2}\int\frac{d^{3}p}{(2\pi)^{3}}\bar{\bm{L}}_{3}(\bm{k},\bm{p},-\bm{p})P_{\mathrm{L}}(p)+\cdots$$

$$+\int \frac{d^{3}p}{(2\pi {)}^3}{c}_X^{(2)}(\mathbf{k},\mathbf{p}){\bar{\mathbf{L}}}_1(-\mathbf{p})P_{\mathrm{L}}(p)$$

$$\displaystyle\hskip 48.0pt+2\int\frac{d^{3}p}{(2\pi)^{3}}c_{X}^{(1)}(p)\bar{\bm{L}}_{2}(\bm{k},\bm{p})P_{\mathrm{L}}(p)+\cdots\Biggr{]},$$

$$c_X^{(n)}({\mathbf{k}}_1,\dots ,{\mathbf{k}}_n)$$

$$\equiv \sum _{m=0}^{\infty }\frac{1}{m!}\int \frac{d^3{p}_1}{(2\pi {)}^3}\cdots \frac{d^3{p}_m}{(2\pi {)}^3}{b}_X^{\mathrm{L}(n+m)}({\mathbf{k}}_1,\dots ,{\mathbf{k}}_n,{\mathbf{p}}_1,\dots ,{\mathbf{p}}_m)$$

$$\times ⟨{\delta }_{\mathrm{L}}({\mathbf{p}}_1)\cdots {\delta }_{\mathrm{L}}({\mathbf{p}}_m)⟩$$

$$=\sum _{m=0}^{\infty }\frac{1}{2^{m}m!}\int \frac{d^3{p}_1}{(2\pi {)}^3}\cdots \frac{d^3{p}_m}{(2\pi {)}^3}$$

$$\times b_X^{\mathrm{L}(n+2m)}({\mathbf{k}}_1,\dots ,{\mathbf{k}}_n,{\mathbf{p}}_1,-{\mathbf{p}}_1,\dots ,{\mathbf{p}}_m,-{\mathbf{p}}_m)$$

$$\times P_{\mathrm{L}}(p_1)\cdots P_{\mathrm{L}}(p_m).$$

$$c_X^{(1)}(k)=b_X^{\mathrm{L}(1)}(k)+\frac{1}{2}\int \frac{d^{3}p}{(2\pi {)}^3}{b}_X^{\mathrm{L}(3)}(\mathbf{k},\mathbf{p},-\mathbf{p})P_{\mathrm{L}}(p)$$

$$+\frac{1}{8}\int \frac{d^{3}p}{(2\pi {)}^3}\frac{d^3{p}^{\prime }}{(2\pi {)}^3}{b}_X^{\mathrm{L}(5)}(\mathbf{k},\mathbf{p},-\mathbf{p},{\mathbf{p}}^{\prime },-{\mathbf{p}}^{\prime })$$

$$\times P_{\mathrm{L}}(p)P_{\mathrm{L}}(p^{\prime })+\cdots ,$$

$$c_X^{(2)}({\mathbf{k}}_1,{\mathbf{k}}_2)=b_X^{\mathrm{L}(2)}({\mathbf{k}}_1,{\mathbf{k}}_2)$$

$$+\frac{1}{2}\int \frac{d^{3}p}{(2\pi {)}^3}{b}_X^{\mathrm{L}(4)}({\mathbf{k}}_1,{\mathbf{k}}_2,\mathbf{p},-\mathbf{p})P_{\mathrm{L}}(p)+\cdots ,$$

$$c_X^{(1)}(k)=\left(b_{10}+b_{01}{k}^2\right)W(kR),$$

$$c_X^{(2)}({\mathbf{k}}_1,{\mathbf{k}}_2)=\{b_{20}+b_{11}\left({k_1}^2+{k_2}^2\right)+b_{02}{k_1}^2{k_2}^2$$

$$-2{\chi }_1({\mathbf{k}}_1\cdot {\mathbf{k}}_2)+{\omega }_{10}\left[3({\mathbf{k}}_1\cdot {\mathbf{k}}_2{)}^2-{k_1}^2{k_2}^2\right]\}$$

Full text

I

Introduction

The large-scale structure (LSS) of the universe is one of the main sources of cosmological information. The statistical properties of LSS constrain the physics of the early universe as well as the nature of dark matter, dark energy, etc. The large-scale structure is mainly probed by observable objects such as galaxies, quasars, and other astronomical objects. The spatial distribution of objects is not the same as that of dark matter; the observable objects are biased tracers of matter distributions. The precise relation between the objects and the dark matter is complicated because of, e.g., nonlinear, nongravitational physics, etc. On large scales where linear theory is applicable, it is a common practice to simply apply a linear bias. However, the linear bias is too simplistic to extract maximum information from the LSS which experiences nonlinear effects of gravitational evolution.

Modeling the nonlinear biasing is a theoretically nontrivial problem in general. One of the natural formulations of the bias modeling is provided by the peak approach BBKS . In this approach, the density peaks of the initial density field are identified as formation sites of astronomical objects. The peaks model predicts a linear velocity bias of halos Des08 ; DS10 ; DCSS10 ; Bal15 . Unlike the prediction of the coupled-fluids approximation for the coevolution of dark matter and halos Fry96 ; TP98 ; CSS12 , the velocity bias in the peaks remains constant in time. This property is explained by the modification of the halo momentum conservation equation in the presence of peak constraints BDKR14 .

The time evolution of the statistical properties of peaks with the velocity bias has been investigated mainly by applying the Zel’dovich approximation Zel70 . The Zel’dovich approximation corresponds to the first-order approximation in terms of the Lagrangian perturbation theory Buc89 ; Mou91 ; Buc92 . A Lagrangian perturbation theory in the presence of general bias is given by the integrated perturbation theory (iPT) Mat11 ; Mat14 . In this formalism, any form of the bias model in Lagrangian space can be taken into account MD16 , including the halo bias MW96 ; MJW97 , peaks model BBKS ; LMD15 , and excursion set peaks AP90 ; PS12 . Therefore, the dynamical evolution of peaks can be described by iPT including higher-order effects of perturbations beyond the Zel’dovich approximation.

In this paper, we apply the iPT to investigate the properties of peak evolution, focusing on the origin of velocity bias in the peaks model. We see how the velocity bias appears in the formalism of iPT with peak bias and derived expressions of two- and three-point propagators which are ingredients to predict the one-loop approximation of the power spectrum. In the limit of the Zel’dovich approximation, previously known results of the velocity bias of peaks are reproduced. We give an analytic expression of the velocity dispersion of peaks in the one-loop approximation. Then we consider the velocity bias from a more formal point of view and show that the velocity bias originates from the “flat constraint” in the peaks model. The flat constraint is the condition that the first spatial derivatives of the field should vanish at the location of peaks. That is, any bias model with the flat constraint, such as the peaks model and excursion set peaks, should predict the velocity bias. We consider a resummation theory of the flat constraint and derive resummed propagators of density and velocity.

In this way, the origin and properties of velocity bias in the peaks model are theoretically investigated in this work. Our paper is organized as follows. In Sec. II, the presence of linear velocity bias is derived in the framework of iPT. The two- and three-point propagators of density and velocity, which are ingredients of predicting the power spectrum, are derived. In Sec. III, the one-loop expression of velocity dispersion is derived. The resulting expression is given by combinations of one-dimensional fast Fourier transforms (FFT). In Sec. IV, a formal aspect of the higher-order perturbation theory in the presence of velocity bias is investigated. Conclusions are summarized in Sec. V. Detailed derivations of necessary functions are given in Appendix A.

II

Velocity bias in the peaks

In this section, we show how the velocity bias is derived in the iPT of lower-order perturbations.

II.1

Velocity of peaks and iPT

In the bias models such as the peaks model, the biased objects are discretely distributed in space. In this case, the velocity field of biased objects is not a continuous field. There is not any value of velocity at every location where discrete objects do not reside; the velocities of the discrete objects are defined only at the locations of objects. Therefore, a natural quantity to describe the velocity of discrete objects are the number-weighted velocity field, or the momentum,

[TABLE]

where $\delta_{X}(\bm{x})$ is the number density contrast of objects $X$ at a location $\bm{x}$, and $\bm{v}(\bm{x})$ is the peculiar velocity at the same location. Throughout this paper, time dependencies are omitted in the arguments of various variables, although they actually depend on time.

For the Lagrangian-space counterpart, we have

[TABLE]

where the Lagrangian velocity field $\bm{v}^{\mathrm{L}}$ is defined by $\bm{v}^{\mathrm{L}}(\bm{q})\equiv\bm{v}(\bm{x})$, where $\bm{x}=\bm{q}+\bm{\varPsi}(\bm{q})$ and $\bm{\varPsi}(\bm{q})$ is the displacement field. The Eulerian and Lagrangian momenta are related by

[TABLE]

The above relation is derived by noting a relation $[1+\delta_{x}(\bm{x})]d^{3}x=[1+\delta^{\mathrm{L}}_{X}(\bm{q})]d^{3}q$. The Lagrangian velocity $\bm{v}^{\mathrm{L}}$ in Eq. (2) is given by

[TABLE]

where $a$ is the scale factor, and the dot denotes a derivative with respect to time, $\dot{\bm{\varPsi}}=\partial\bm{\varPsi}/\partial t$. The momentum $\bm{j}_{X}(\bm{x})$ is nonzero only at the locations of discrete objects, so Eq. (4) is evaluated only at locations of discrete objects as well (at locations of peaks, in the case of the peaks model). Thus the velocity of discrete objects is statistically biased with respect to the velocity of mass, even when one assumes that velocities of matter and objects are the same at locations of discrete objects.

Throughout this paper, it is assumed that the peculiar velocity of the peaks is the same as for the dark matter field, at the peak location. While this is a critical assumption in this paper, it is not necessarily held for actual objects such as halos with relatively small masses LP11 , in which case the velocities of tracers can depend on the environment, as well as on the underlying smoothing procedure in general. The inclusion of such dependences is beyond the scope of this paper, and we assume the velocities of biased objects are well approximated by those of peaks.

Using Eqs. (2)–(4), statistics of momentum field (and cross correlations to the density field as well) can be calculated in the way of iPT Mat11 . The full formulation of extending the iPT including the momentum field is developed in Sec. IV. In the rest of this section, we focus on its implication to the Lagrangian velocity bias.

II.2 Velocity bias in Lagrangian space

The Fourier transform of Eq. (2) is given by

[TABLE]

where we use a notation

[TABLE]

Taking a cross correlation with linear density field $\delta_{\mathrm{L}}(\bm{k})$, we have

[TABLE]

Although the second term on the right-hand side (RHS) corresponds to a higher-order correction in the usual sense of perturbation theory, what we see below is that this term gives a contribution which linearly biases the velocities of halos in the peaks model.

Equation (7) is to be evaluated by means of iPT. For the purpose that the results can be compared with those in the literature, we consider a model where the displacement field is smoothed by a filtering scale $R$, in accordance with most of the literature in which smoothed underlying density and velocity fields are considered. In this case, the perturbative expansion of the displacement field is given by

[TABLE]

where $\bar{\bm{L}}_{n}$ is the smoothed, symmetrized kernel of $n$th-order Lagrangian perturbation theory (LPT).

There are at least two schemes of smoothing the LPT kernels. The first one is to smooth the linear density contrast $\delta_{\mathrm{L}}$ in Lagrangian space, and the second one is to smooth the nonlinear displacement field $\bm{\varPsi}$. Accordingly, the relation between the LPT kernel $\bm{L}_{n}$ and the smoothed kernel $\bar{\bm{L}}_{n}$ is given by either

[TABLE]

or

[TABLE]

The first-order kernel does not have any difference between the two smoothing schemes above, and is given by

[TABLE]

In the following, we assume this relation for the first-order kernel, and we leave the expression $\bar{\bm{L}}_{n}$ for higher-order kernels. The LPT kernels up to third order are given by Cat95 ; CT96 ; Mat15

[TABLE]

where

[TABLE]

Perturbative expansions of other variables are given by

[TABLE]

where $Hf=\dot{D}/D$, and $D$ is the linear growth factor. We have used the approximate time dependence of the kernel function, $\bm{L}_{n}\propto D^{n}$ Mat15 , and an identity $d(D^{n})/dt=nHfD^{n}$. Substituting the expansion of Eqs. (16) and (17) into Eq. (7), and assuming Gaussian initial conditions, we have

[TABLE]

where the higher-order self-loops of the bias functions are renormalized according to the iPT. The functions $c_{X}^{(n)}$ are the renormalized bias functions given by Mat11

[TABLE]

The first line is derived from the general definition [Eq. (64) below], and the second line is derived with Gaussian initial conditions. For example,

[TABLE]

and so forth. We have used the fact that the self-contractions vanish $\langle\bm{v}^{\mathrm{L}(n)}\rangle=\langle\delta_{X}^{\mathrm{L}(n)}\rangle=0$ for every order. The renormalized bias functions of the peaks model up to the second order are given by Des13 ; LMD15 ; MD16

[TABLE]

where coefficients $b_{ij}$, $\chi_{1}$ and $\omega_{10}$ are scale-independent constants which depend on the threshold of the peaks (see Refs. LMD15 ; MD16 for their definitions).

The second and fourth terms in the square bracket of Eq. (18) are the usual mode-coupling terms. In the third term, however, we have $\bar{\bm{L}}_{1}(\bm{p})=W(pR)\bm{p}/p^{2}$, and the integration over the angle of $\bm{p}$ leaves the term with only $\chi_{1}$ of Eq. (II.2), resulting in

[TABLE]

where we have used $(4\pi)^{-1}\int d\varOmega_{p}p_{i}p_{j}/p^{2}=\delta_{ij}/3$ and $\chi_{1}=-3/(2{\sigma_{1}}^{2})$, ${R_{\mathrm{v}}}^{2}={\sigma_{0}}^{2}/{\sigma_{1}}^{2}$, and spectral parameters $\sigma_{n}(R)$ are defined by

[TABLE]

The above Eq. (24) has the same functional form as the linear term. Thus, Eq. (18) reduces to

[TABLE]

where “$+\,\mbox{m.c.}$” represents mode-coupling terms. The correlator for the matter momentum field $\langle\delta_{\mathrm{L}}\bm{j}_{\mathrm{m}}^{\mathrm{L}}\rangle$ is just given by putting ${R_{\mathrm{v}}}^{2}=0$ and $W(kR)=1$ in the above expression. Thus we have

[TABLE]

In this way, the $\chi_{1}$ term in $c_{X}^{(2)}$ of the peaks model introduces a linear velocity bias. The scale-dependent velocity bias factor, $(1-{R_{\mathrm{v}}}^{2}k^{2})W(kR)$, is consistent with the finding of Ref. Des08 . Equation (27) corresponds to Eqs. (2) and (9) of Ref. Bal15 . In the following, we use the notation,

[TABLE]

which represents the lowest-order effect of velocity bias in Lagrangian space.

II.3 Relation to the one-loop density power spectra of iPT

The same type of the seemingly linear term is also contained in a one-loop correction of iPT power spectrum. In the iPT, the normalized two-point propagator is given by Mat14

[TABLE]

where the displacement field is filtered by a smoothing scale $R$ in accordance with the model of the previous subsection, and we drop a term which vanishes for an isotropic function $c^{(1)}_{X}(p)$. The second and third terms in the integrand are mode-coupling terms. The first term has essentially the same form as the third term of Eq. (18), which is the very source of the velocity bias. After the integration, this term is equal to $-{R_{\mathrm{v}}}^{2}k^{2}\bm{k}\cdot\bar{\bm{L}}_{1}(\bm{k})$. Hence, the first term in the curly bracket of Eq. (29) gives a term without mode coupling, and we have

[TABLE]

In real space without the redshift-space distortions, we have $\bm{L}_{1}(\bm{k})=\bm{k}/k^{2}$. Consequently, the Eulerian linear bias factor is effectively given by

[TABLE]

Therefore, the $\chi_{1}$ term again modifies the linear bias, not only for the velocity field, but also for the density field.

In Eq. (31), the magnitude of the scale-dependent part $\propto k^{2}$ of the linear bias factor is reduced compared with the prediction of the linear theory. This effect can be interpreted from the presence of velocity bias, because the velocities of halos are slower than those of dark matter on average, and halos are less clustered than in the case without velocity bias. More quantitatively, one can see this as follows: consider the simplest case when the halos are not linearly biased, $b_{10}=b_{01}=0$, and formations and destructions of halos can be neglected at the lowest order. In this ideal case, the number density of halos are approximately conserved at the lowest order, $\dot{\delta}_{X}+a^{-1}i\bm{k}\cdot\bm{j}_{X}\sim 0$. Because of the velocity bias, we have $\bm{k}\cdot\bm{j}_{X}\sim(1-{R_{\mathrm{v}}}^{2}k^{2})W(kR)\bm{k}\cdot\bm{j}_{\mathrm{m}}$. Combining these relations with the mass conservation equation $\dot{\delta}_{\mathrm{m}}+a^{-1}i\bm{k}\cdot\bm{j}_{\mathrm{m}}=0$, we have a relation, $\dot{\delta}_{X}\sim(1-{R_{\mathrm{v}}}^{2}k^{2})W(kR)\dot{\delta}_{\mathrm{m}}$, which indicates $b(k)\sim(1-{R_{\mathrm{v}}}^{2}k^{2})W(kR)$, and this is consistent to Eq. (31) within our assumption $b_{10}=b_{01}=0$. In this way, the amplitude of the scale-dependent part of the linear bias factor is reduced due to the velocity bias. All such effects are already included in the iPT with the peak constraint.

II.4 Propagators in Lagrangian space

The correlator of Eq. (18) is related to the two-point propagators BCS08 . We define the two-point propagator of the momentum field in Lagrangian space, $\bm{\varGamma}^{\mathrm{vL}(1)}_{X}(\bm{k})$, by

[TABLE]

where $\delta/\delta\delta_{\mathrm{L}}(\bm{k})$ is the functional derivative with respect to the linear density field $\delta_{\mathrm{L}}(\bm{k})$ in Fourier space, and the appearance of the Dirac’s delta function is a consequence of translation invariance of space. The imaginary unit $i$ is put in the definition, so as to make the propagator real vector. In Gaussian initial conditions, we have a relation,

[TABLE]

where $\langle\cdots\rangle_{\mathrm{c}}$ represents the cumulants of the second order, and we have $\bm{\varGamma}^{\mathrm{vL}(1)}_{X}(-\bm{k})=-\bm{\varGamma}^{\mathrm{vL}(1)}_{X}(\bm{k})$ from the parity symmetry.

Including the mode-coupling terms of Eqs. (18) and (26), the two-point correlators of the momentum field in Lagrangian space is given by

[TABLE]

Similarly, we consider the three-point propagator in Lagrangian space, $\bm{\varGamma}^{\mathrm{vL}(2)}_{X}(\bm{k}_{1},\bm{k}_{2})$, which is defined by

[TABLE]

In Gaussian initial conditions, we have a relation,

[TABLE]

Thus, it is sufficient to evaluate the left-hand side (LHS) of the above equation to obtain an expression of the three-point propagator.

From Eq. (5), we have

[TABLE]

The cumulants in the above equation are all third order, where the square bracket in the cumulant of the last term is considered as a single variable. Substituting the expansion of Eqs. (16) and (17) into Eq. (37), we have

[TABLE]

where the higher-order self-loops of the bias functions are renormalized according to the iPT.

The fourth term in the square bracket is similar to the integral of Eq. (24). The third-order renormalized bias functions of the peaks model is given in Ref. LMD15 . The terms of odd parity with respect to the third argument in their expression are given by111In the $c_{00001}$ term of Eq. (91) of Ref. LMD15 , the term $(\bm{k}_{1}\cdot\bm{k}_{2})k_{3}^{2}$ should be replaced by $(\bm{k}_{1}\cdot\bm{k}_{2})^{2}k_{3}^{2}$, i.e., this term has even parity.

[TABLE]

The coefficients can be calculated by their definition, and the results are given by $c_{10100}=-3b_{10}/2{\sigma_{1}}^{2}$, $c_{01100}=-3b_{01}/2{\sigma_{1}}^{2}$. Other terms with even parity do not contribute to the integral. Performing the similar calculation of Eq. (24), we have

[TABLE]

Consequently, the three-point propagator of Eq. (35) is given by

[TABLE]

up to the lowest-order approximation.

II.5 The displacement field

For the purpose of comparison, we consider the statistics of the displacement field. In the calculations above, the momentum field can be replaced by a displacement field, $a\dot{\bm{\varPsi}}\rightarrow\bm{\varPsi}$, and we can follow almost the same steps. The number-weighted displacement field $\bm{\psi}_{X}^{\mathrm{L}}$ in Lagrangian space is defined by

[TABLE]

and the perturbative expansions of the displacement field in Lagrangian space is given by Eq. (8). Comparing these equations with Eqs. (2), (4) and (17), we see that replacing the Lagrangian kernels $\bar{\bm{L}}_{n}\rightarrow\bar{\bm{L}}_{n}/naHf$ in the expression of the above results for $\bm{j}_{X}^{\mathrm{L}}$ can give the results for $\bm{\psi}_{X}^{\mathrm{L}}$. Following the calculations to derive Eqs. (34) and (41), the two- and three-point propagators of displacement field are derived as

[TABLE]

III

The velocity dispersion of peaks

III.1 Nonlinear corrections to the velocity dispersion of peaks

As an example of applications of the velocity bias derived from the framework of iPT in the previous subsection, we consider the velocity dispersion of peaks in this section. The velocity dispersion of biased objects is defined by

[TABLE]

where $\bm{v}(\bm{x}_{a})$ is the Eulerian velocity at Eulerian coordinates $\bm{x}_{a}$ of an object labeled by $a$, and $N_{\mathrm{obj}}$ is the total number of biased objects in a sample. In the last expression, $\bm{v}^{\mathrm{L}}(\bm{q}_{a})$ is the Lagrangian velocity at Lagrangian coordinates $\bm{q}_{a}$. The last equality holds because the Lagrangian velocity field is defined by $\bm{v}^{\mathrm{L}}(\bm{q})\equiv\bm{v}(\bm{x})$, where $\bm{x}=\bm{q}+\bm{\varPsi}(\bm{q})$. Therefore, there are not any differences between velocity dispersions of objects in Eulerian space and in Lagrangian space.

The velocity dispersion of objects is not identical to the dispersion of velocity field $\langle|\bm{v}(\bm{x})|^{2}\rangle$. Instead, it is given by a number-weighted dispersion of velocity field, $\langle n_{X}(\bm{x})|\bm{v}(\bm{x})|^{2}\rangle/\bar{n}_{X}$. The velocity dispersion of Eq. (45) is given by an expression,

[TABLE]

This expression can be evaluated by applying the perturbative expansions of Eqs. (16) and (17), and by following the similar calculations to derive the correlators of Eqs. (18) and (41). Just as in that case, only terms of odd parity in the higher-order bias functions $c^{(2)}_{X}$ and $c^{(3)}_{X}$ survive and they are represented by ${R_{\mathrm{v}}}^{2}$ and $c^{(1)}_{X}$ because of Eqs. (24) and (40). Alternatively, Eq. (46) is equivalent to an expression $\sigma_{\mathrm{nv}}^{2}=\langle\bm{v}^{\mathrm{L}}\cdot\bm{j}^{\mathrm{L}}_{X}\rangle$, and one can evaluate this expression by applying the perturbative expansion of Eq. (17) and using the results of the correlators, Eqs. (34) and (41). In either way, one can derive the same result,

[TABLE]

where

[TABLE]

This parameter $\sigma_{\mathrm{dpk}}^{2}$ is the “peak displacement dispersion” defined in Ref. Bal15 .

III.2 Reducing the dimensions of integral

The six-dimensional integrals in Eq. (47) can be evaluated by a combination of one-dimensional integrals as shown below. Recently it is pointed out that the multidimensional integrations appeared in the perturbation theory of the nonlinear power spectrum can be evaluated by combining only one-dimensional (1D) integrations of Hankel transforms with FFT SVM16 ; SV16 ; MFHB16 ; FBMH16 . The essence of the method is the realization that angular parts of the multidimensional integrations can be analytically performed, and all the remaining integrations can be represented by a set of 1D Hankel transforms. Essentially the same method can be applied here.

The particular forms of LPT kernels in Eq. (47) are given by

[TABLE]

where we assume the smoothing scheme of linear density field, Eq. (9). When the smoothing scheme of displacement, Eq. (10), is adopted, the products of smoothing kernels in Eqs. (49)–(51) are, respectively, replaced by $W^{2}(kR)W^{2}(pR)\rightarrow W^{2}(|\bm{k}+\bm{p}|R)$, $W(kR)W(pR)\rightarrow W(|\bm{k}+\bm{p}|R)$, and $W(kR)W^{2}(pR)\rightarrow W(kR)$.

The angular integration of the variable $\bm{p}$ can be evaluated by a following formula for an arbitrary function $F(k)$ and Legendre polynomials $\mathsf{P}_{l}(x)$:

[TABLE]

which can be shown, e.g., by applying a 3D Fourier transform of $F(k)$, the plane wave expansion of $e^{i\bm{k}\cdot\bm{r}}$, the addition formula of Legendre polynomials, $\mathsf{P}_{l}(x)$, and the orthogonality relation of spherical harmonics, $Y_{l}^{m}(\varOmega)$. The last integral is a Hankel transform of $F(k)$. For example, for $F(k)=1/k^{2}$, we have

[TABLE]

Equations (52) and (53) are sufficient to evaluate the angular integrations of Eqs. (49)–(51), by changing the integration variable as $\bm{p}\rightarrow-\bm{p}$ when necessary. The angular variable $\mu=\bm{k}\cdot\bm{p}/kp$ in those equations are represented by Legendre polynomials as $(1-\mu^{2})^{2}=(8/15)\mathsf{P}_{0}(\mu)-(16/21)\mathsf{P}_{2}(\mu)+(8/35)\mathsf{P}_{4}(\mu)$, $1-\mu^{2}=(2/3)\mathsf{P}_{0}(\mu)-(2/3)\mathsf{P}_{2}(\mu)$ and $\mu-\mu^{3}=(2/5)\mathsf{P}_{1}(\mu)-(2/5)\mathsf{P}_{3}(\mu)$.

Using the above formula, angular integrations of Eqs. (49)–(51) are evaluated as

[TABLE]

where the first smoothing scheme, Eq. (9) is adopted. Substituting the above results into Eq. (47), we have

[TABLE]

where

[TABLE]

The number-weighted dispersion of displacement field, $\sigma_{\mathrm{nd}}^{2}\equiv\langle(1+\delta^{\mathrm{L}}_{X})|\bm{\varPsi}|^{2}\rangle$, is just obtained by replacements $\bar{\bm{L}}_{n}\rightarrow\bar{\bm{L}}_{n}/naHf$ in the expression of the number-weighted dispersion of the velocity field. Corresponding to Eq. (48), we have

[TABLE]

Correspondingly, Eq. (57) can easily be modified to give the formula for $\sigma_{\mathrm{nd}}^{2}$ by replacing the coefficients in front of three integrals as $18/49\rightarrow 9/98$, $12/7\rightarrow 6/7$ and $30/21\rightarrow 10/21$.

IV General formulation of velocity bias with

the flat constraint

IV.1

A class of bias models with the flat constraint

The appearance of the linear velocity bias from a second term in Eq. (5) is due to the $\chi_{1}$ term in the second-order renormalized bias function of Eq. (II.2), which is an odd function of the wave vector $\bm{p}$. The same is also the case for the linear bias factor of density field in Eq. (31). Thus it is crucial that the $\chi_{1}$ term is present in the renormalized bias function to have the effects of velocity bias.

The $\chi_{1}$ factor arises from the flat constraint of peaks, $\bm{\eta}=\bm{0}$, where $\bm{\eta}=\bm{\nabla}\delta_{R}/\sigma_{1}$ Des13 ; LMD15 ; MD16 . Therefore, the flat constraint is an essential ingredient for the appearance of velocity bias. This is simply interpreted as follows: at a point with the flat constraint, density gradients are zero, and the magnitude of velocity at the same point is expected to be smaller than the average value. Thereby the velocity is statistically biased in the presence of the flat constraint.

The peak constraint contains the flat constraint, and the velocity bias appears from the $\chi_{1}$ term in the lowest-order approximation in the last section. Any other model with the flat constraint, such as the excursion set peaks (ESP) AP90 ; PS12 , is naturally expected to have the velocity bias of the same kind. In the following, we consider a general situation in the presence of the flat constraint, in order to investigate the origin and properties of the velocity bias.

For the sake of generality, we assume a general constraint in the Lagrangian number density of objects which contains the flat constraint,

[TABLE]

In the above notation, $\nu$ is a scalar, $\bm{\eta}$ is a three-dimensional vector, and $\bm{\zeta}$ is a $3\times 3$ tensor, which components are defined by

[TABLE]

The function $F$ is an arbitrary function of the above variables at a location $\bm{x}$ where the number density of biased objects $n_{X}(\bm{x})$ is defined. From the rotational symmetry, the variable $\bm{\eta}$ is uncorrelated with variables $\nu$ and $\bm{\zeta}$: $\langle\nu\eta_{i}\rangle=\langle\eta_{i}\zeta_{jk}\rangle=0$. We assume the function $F$ can depend on other variables which are assumed to be uncorrelated with $\bm{\eta}$. For example, the variable $\bm{\eta}$ is correlated with the third-order derivatives $\partial_{i}\partial_{j}\partial_{k}\delta_{R}$, which is assumed to be absent in the function $F$. The ESP model has an additional variable $\partial\delta_{R}/\partial R$, and this variable is really uncorrelated with $\bm{\eta}$ because of the rotational symmetry. Thus the ESP model is an example with the general constraint of Eq. (62).

In the general formalism below, different smoothing kernels can be adopted for every variable, $\nu$, $\bm{\eta}$, $\bm{\zeta}$, $\ldots$. For example, both the top-hat kernel and the Gaussian kernel are adopted in a certain version of the ESP model. The essential assumption below is that $\bm{\eta}$ is the only variable which has odd parity. Other variables are inevitably uncorrelated to $\bm{\eta}$ at a single point by parity symmetry.

IV.2

Implications of the flat constraint and the renormalized bias functions

The observable quantities, such as the power spectrum, depends on the series of renormalized bias functions. The renormalized bias function $c_{X}^{(n)}$ is defined by Mat12

[TABLE]

where $\delta_{X}^{\mathrm{L}}(\bm{k})$ is the Fourier transform of the density contrast of the biased objects in Lagrangian space, and $\bm{k}_{1\cdots n}\equiv\bm{k}_{1}+\cdots+\bm{k}_{n}$. The density contrast in Lagrangian space is given by $n_{X}(\bm{x})/\bar{n}_{X}-1$, where $\bar{n}_{X}=\langle n_{X}\rangle$ is the mean number density of biased objects.

In a class of models in which the number density $n_{X}(\bm{x})$ is given by a function of finite number of variables, $y_{\alpha}(\bm{x})$, which are linearly related to the linear density field $\delta_{\mathrm{L}}$, the renormalized bias function in Eq. (64) is given by Mat11

[TABLE]

where $U_{\alpha}(\bm{k})$ are the Fourier coefficients of the variables $y_{\alpha}(\bm{x})$,

[TABLE]

Defining an operator

[TABLE]

Eq. (65) can be represented by

[TABLE]

In Eqs. (65) and (68), the average $\langle\cdots\rangle$ is taken over the random variables $\bm{y}=(y_{\alpha})$. We assume Gaussian initial conditions throughout this paper, and the distribution function is given by Dor70 ; BBKS ; PGP09 ; GPP12 ; Des13 ; LMD15

[TABLE]

In the above equation, we introduce notations

[TABLE]

where repeated indices are summed over, and

[TABLE]

is the traceless part of $\zeta_{ij}$. The factor $N_{0}$ is a normalization constant to ensure that the total probability is equal to one, but its actual value is not used for our applications in this paper.

Although the RHS of Eq. (69) is represented by rotationally invariant variables $\nu$, $J_{1}$ and $J_{2}$, the probability distribution function $P(\bm{y})$ is still for the linear variables $\bm{y}$, and not for the invariant variables.

In our class of models, Eq. (62), the variables $y_{\alpha}$ are given by

[TABLE]

and

[TABLE]

where $i\leq j$. When different smoothing functions are applied to each variable, the window function $W(kR)$ in the above equation is replaced by corresponding functions. Since we assume the variable $\bm{\eta}$ is uncorrelated to other variables, the only components which have odd parity, $U_{\alpha}(-\bm{k})=-U_{\alpha}(\bm{k})$, are the components $ik_{i}W(kR)/\sigma_{1}$ in Eq. (73). All the other components are assumed to have even parity, $U_{\alpha}(-\bm{k})=+U_{\alpha}(\bm{k})$.

It is an essential assumption in the bias models analyzed in this paper that the variable $\bm{\eta}$ is the only quantity which has odd parity in the set of variables $(y_{\alpha})$. This is the case for the peaks model and is the reason why only the term with $\chi_{1}$ in $c_{X}^{(2)}$ survives and causes the velocity bias as observed in Sec. II. It is essential that the term with $\chi_{1}$ in $c_{X}^{(2)}(\bm{k},\bm{p})$ of Eq. (II.2) is the only term which is an odd function of $\bm{p}$. All the other terms are even functions of $\bm{p}$, so that they vanish in the angular integral of $c_{X}^{(2)}(\bm{k},\bm{p})\bm{L}_{1}(-\bm{p})$. In order to give a general formalism in this section, we pursue similar mechanisms in higher-order terms.

For that purpose, we decompose the operator $\hat{D}(\bm{k})$ of Eq. (67) as

[TABLE]

where

[TABLE]

corresponds to an operator of even parity, $\hat{D}_{0}(-\bm{k})=+\hat{D}_{0}(\bm{k})$, and

[TABLE]

corresponds to an operator of odd parity, $\hat{D}_{1}(-\bm{k})=-\hat{D}_{1}(\bm{k})$.

We consider the following sum of integrals:

[TABLE]

This form of integral appears in the applications of iPT, and a corresponding section of the iPT diagram is given in the left diagram of Fig. 1. For the diagrammatic rules of iPT, see Refs. Mat11 ; Mat14 . The corresponding iPT diagram is denoted by a double square as shown in Fig. 1. This iPT diagram has another resummation factor $\varPi(\bm{k})k_{i_{1}}\cdots k_{i_{m}}$ in addition to the function $T_{X}^{(n)}$, where $\varPi(\bm{k})$ is the resummation factor of the displacement field222The resummation of the iPT explicitly violates the “extended Galilean invariance” PP16 in the small-scale limit. Since the iPT is the theory for the weakly nonlinear regime and should not be applied on the scales below the shell-crossing scale, the violation is not considered to be a practical issue of iPT. Mat11 . The exact diagrammatic rule for this vertex of the double square is also shown in Fig. 1.

Because $\bm{L}_{1}(\bm{p})=\bm{p}/p^{2}$ is an odd function of $\bm{p}$, components of even parity in $c^{(n+m)}_{X}$ with respect to $\bm{p}_{1},\ldots,\bm{p}_{m}$ vanish. Substituting the expression of Eq. (68) into Eq. (77), $\hat{D}(\bm{p}_{i})$’s are replaced by $\hat{D}_{1}(\bm{p}_{i})$’s. Using an identity,

[TABLE]

we have

[TABLE]

This expression reduces to a more useful form. In Appendix A, it is shown that the function $T_{X}^{(n)}(\bm{k};\ldots)$ generally has a factor $\exp(k^{2}{\sigma_{0}}^{4}/6{\sigma_{1}}^{2})$, and explicit forms for $n=0,1,2$ are derived. We define the normalized function $\hat{T}_{X}^{(n)}$ by

[TABLE]

and the results for $n=0,1,2$ are given by Eqs. (129), (132), and (135). Substituting $W(k_{1}R)\bm{k}_{1}={k_{1}}^{2}\bar{\bm{L}}_{1}(\bm{k}_{1})$, $W(k_{2}R)\bm{k}_{2}={k_{2}}^{2}\bar{\bm{L}}_{1}(\bm{k}_{2})$ in these equations, they can be represented by

[TABLE]

The above results are derived in real space. In redshift space, the LPT kernels are replaced by Mat08

[TABLE]

where $\hat{\bm{z}}$ is the unit vector along the line of sight, and $f=d\ln D/d\ln a$ is the linear growth rate. In particular, we have

[TABLE]

where $k_{z}=\bm{k}\cdot\hat{\bm{z}}$ is the line-of-sight component of the vector $\bm{k}$. Since the $\bm{k}$-dependence and LPT kernels in Eq. (77) are included only in the form of Eq. (85), the function $T_{X}^{(n)}$ in redshift space is given by Eq. (77), with the replacement $\bm{k}\rightarrow\bm{k}+fk_{z}\hat{\bm{z}}$, i.e., $T_{X}^{(n)}(\bm{k}+fk_{z}\hat{\bm{z}};\bm{k}_{1},\ldots,\bm{k}_{n})$. On the one hand, this result is equivalent to just replacing the first-order kernel $\bm{L}_{1}$ by $\bm{L}^{\mathrm{s}}_{1}$ in Eqs. (81)–(83). On the other hand, $k^{2}$ is replaced by $k^{2}+f(f+2){k_{z}}^{2}$ in the exponential factor of Eq. (80). Therefore, the function $T_{X}^{(n)}$ in redshift space is given by

[TABLE]

where $\mu_{\bm{k}}\equiv k_{z}/k$, and the function $\hat{T}_{X}^{\mathrm{s}(n)}(\bm{k};\bm{k}_{1},\ldots,\bm{k}_{n})$ is given by Eqs. (81)–(83) with a replacement $\bar{\bm{L}}_{1}\rightarrow\bar{\bm{L}}^{\mathrm{s}}_{1}$.

IV.3

Propagators of density field with the flat constraint

In this section, we consider propagators $\varGamma_{X}^{(n)}(\bm{k}_{1},\ldots,\bm{k}_{n})$ of the density field with the flat constraint in the formalism of iPT. Using the propagators, the density power spectrum $P_{X}(\bm{k})$ of biased object $X$ is given by

[TABLE]

in the case of the Gaussian initial condition. In more general cases with primordial non-Gaussianity, other terms with the linear bispectrum, trispectrum, etc., are added to the above expansion Mat11 ; Mat12 ; YM13 . The effects of primordial non-Gaussianity also modify the the calculations of the renormalized bias functions in Sec. IV.2. In this paper, we only consider the case of Gaussian initial conditions throughout.

In the course of evaluating the propagators, special types of loop corrections, depicted by Fig. 1 with the flat constraint, are taken into account. First we consider the two-point propagator. Using a diagrammatic notation of Fig. 1 for the iPT vertex, the two-point propagator, $\varGamma^{(1)}_{X}(\bm{k})$, is given by Fig. 2 up to the one-loop approximation.

It is straightforward to apply the diagrammatic rules of iPT Mat11 and a new rule of Fig. 1 with the flat constraint for a resummed vertex. The two-point propagator is represented by $\varGamma_{X}^{(1)}(\bm{k})=\varPi(\bm{k})\hat{\varGamma}_{X}^{(1)}(\bm{k})$, where $\varPi(\bm{k})$ is the resummation factor of the displacement field, and $\hat{\varGamma}_{X}^{(1)}(\bm{k})$ is the normalized two-point propagator. The resummation factor in real space, $\varPi(k)$, is given by Mat08 ; Mat11

[TABLE]

up to the one-loop order.

Applying diagrammatic rules to Fig. 2, the normalized two-point propagator up to the one-loop approximation is given by

[TABLE]

Because the function $T_{X}^{(n)}$ also has a common factor as in Eq. (80), the propagators in general have a common factor

[TABLE]

where $\sigma_{\mathrm{dpk}}^{2}$ is the peak displacement dispersion defined by Eq. (48). Substituting Eqs. (81)–(83) into Eq. (89), we obtain

[TABLE]

Comparing the above equation with Eq. (29), the resummation of the flat constraint in the bias replaces the first-order kernel $\bar{\bm{L}}_{1}(\bm{k})$ by $(1-{R_{\mathrm{v}}}^{2}k^{2})\bar{\bm{L}}_{1}(\bm{k})$, removes the second-order bias function $c^{(2)}_{X}$ in the integrand of the one-loop term, and converts $\sigma_{-1}^{2}$ into $\sigma_{\mathrm{dpk}}^{2}$ in the exponential prefactor.

Neglecting mode-coupling terms in Eq. (91), the effective bias factor in Eulerian space is given by

[TABLE]

The two-point propagator $\varGamma_{X}^{(1)}(k)$ is effectively a bias factor $b_{X}^{\mathrm{eff}}(k)$ in real space. This result agrees with a previous result, Eq. (10) of Ref. Bal15 in the peaks model.

When the factor $\exp(k^{2}{\sigma_{0}}^{4}/6{\sigma_{1}}^{2})$ is dropped, this result is consistent with that of the previous section, Eq. (29). The reason for not showing this factor in the previous section is that this contribution comes from two- and higher-loop corrections in terms of the original diagram of iPT. As one can explicitly show, the diagram of Fig. 1 with even numbers $m=2l$ gives a factor of $(k^{2}{\sigma_{0}}^{4}/6{\sigma_{1}}^{2})^{l}/l!$. The diagram with odd numbers $m=2l+1$ only contributes when $n\geq 1$, and gives a factor of ${R_{\mathrm{v}}}^{2}(k^{2}{\sigma_{0}}^{4}/6{\sigma_{1}}^{2})^{l}/l!$. After the summation over $l$, the factor $\exp(k^{2}{\sigma_{0}}^{4}/6{\sigma_{1}}^{2})$ or ${R_{\mathrm{v}}}^{2}\exp(k^{2}{\sigma_{0}}^{4}/6{\sigma_{1}}^{2})$ emerges. If we take only the cases $m=0$ and $1$ into account, the exponential factor does not appear. Since $m$ is the number of loops in the diagram of Fig. 1, the exponential factor comes from two-, or higher-loop corrections. In this way, loop corrections with the flat constraint contains contributions without mode couplings, which look like a linear term.

The above results are explicitly derived for real space. As noted in the previous subsection, extending the results to redshift space is straightforward. In redshift space, the resummation factor $\varPi(\bm{k})$ is given by Mat08 ; Mat11

[TABLE]

while the exponential factor in $T_{X}^{(n)}$ is given by Eq. (86). Accordingly, the factor $G_{\mathrm{d}}(\bm{k})$ of Eq. (90) in real space is replaced by

[TABLE]

in redshift space. The two-point propagator in redshift space $\varGamma^{\mathrm{s}(n)}_{X}$ is given by the form of Eq. (91) with replacements $\bar{\bm{L}}_{n}\rightarrow\bar{\bm{L}}^{\mathrm{s}}_{n}$ and $G_{\mathrm{d}}(k)\rightarrow G_{\mathrm{d}}^{\mathrm{s}}(\bm{k})$.

In redshift space, the corresponding factor of Eq. (92) is given by

[TABLE]

In the large-scale limit, $k\rightarrow 0$, this expression reduces to the Kaiser’s factor, $b+f{\mu_{\bm{k}}}^{2}=b(1+f{\mu_{\bm{k}}}^{2}/b)$, for the linear redshift-space distortion Kai87 . It is a common practice to define the redshift-space distortion parameter $\beta=f/b$ in the Kaiser’s limit. In Eq. (95), this parameter is effectively scale dependent. The corresponding parameter is given by

[TABLE]

and Eq. (95) is represented by

[TABLE]

The exponential prefactor corresponds to the damping factor along the line of sight, which represents the fingers-of-God effect of biased objects within the Zel’dovich approximation.

Next we consider the three-point propagator, $\varGamma_{X}^{(2)}(\bm{k}_{1},\bm{k}_{2})=\varPi(\bm{k}_{12})\hat{\varGamma}_{X}^{(2)}(\bm{k}_{1},\bm{k}_{2})$, where $\bm{k}_{12}=\bm{k}_{1}+\bm{k}_{2}$. The diagram for this propagator is given in Fig. 3.

In calculating the power spectrum, the three-point propagator is always accompanied by loop diagrams. Thus, we do not need to consider loop diagrams for the three-point propagator for the purpose of evaluating the power spectrum in the one-loop approximation. Applying diagrammatic rules to Fig. 3, the normalized three-point propagator is given by

[TABLE]

where $\bm{k}=\bm{k}_{1}+\bm{k}_{2}$. Substituting Eqs. (81)–(83), we have

[TABLE]

Comparing the above equation with Eq. (16) of Ref. Mat14 , the resummation of the flat constraint in the bias replaces the first-order kernel $\bar{\bm{L}}_{1}(\bm{k})$ by $(1-{R_{\mathrm{v}}}^{2}k^{2})\bar{\bm{L}}_{1}(\bm{k})$ and converts $\sigma_{-1}^{2}$ into $\sigma_{\mathrm{dpk}}^{2}$ in the exponential prefactor. The redshift-space counterpart is again obtained by replacements, $\bar{\bm{L}}_{n}\rightarrow\bar{\bm{L}}^{\mathrm{s}}_{n}$ and $G_{\mathrm{d}}(k)\rightarrow G_{\mathrm{d}}^{\mathrm{s}}(\bm{k})$, in the above expression.

Substituting Eqs. (91) and (99) into Eq. (87), the expression for the power spectrum $P_{X}(\bm{k})$ is obtained. In the case of Zel’dovich approximation, $\bar{\bm{L}}_{2}=\bar{\bm{L}}_{3}=\cdots=0$, the resulting expression is consistent with a previous result of Ref. BD17 in the peaks model.

IV.4

Propagators of momentum field with the flat constraint

Extending the method of the previous subsection, we consider the propagators of the velocity field in this section. For this purpose, we need to have the diagrammatic rules of iPT for the velocity momentum, $\bm{j}_{X}$. One can straightforwardly follow the derivation of Ref. Mat11 for the diagrammatic rules of iPT. The Fourier transform of the momentum field $\bm{j}_{X}(\bm{x})$, Eq. (3), is given by

[TABLE]

where the expansion of the number density field, Eq. (16) is applied. The expansion of the displacement field $\tilde{\bm{\varPsi}}(\bm{k})$ is given by Eq. (8). From these expansions, the diagrammatic rules for the momentum field are given by Fig. 4.

The rules in Fig. 4 are extensions of those for the number density field $\delta_{X}(\bm{x})$ given in Fig. 7 of Ref. Mat11 . An important difference of the momentum field from the density field is the existence of the dotted lines in the top and bottom graphs. The external, double dotted line corresponds to the biased momentum, $\bm{j}_{X}$. The internal, single dotted line corresponds to the velocity field $\bm{v}=a\dot{\bm{\varPsi}}$. The internal, wavy line corresponds to the displacement field $\bm{\varPsi}$. The solid line corresponds to the linear density field $\delta_{\mathrm{L}}$. The only one dotted line should be attached to the external vertex. Numbers of solid and wavy lines attached to the external vertex are arbitrary, including zero. The index $j$ in the top graph corresponds to the spatial index of the momentum field $\bm{j}_{X}$, and this index appears only in the internal vertex of the bottom graph. Thus the corresponding factor in the rule for the external vertex does not depend on the index $j$ and the external vertex of the top graph only transmits the index to the internal vertex of the bottom graph. All these properties are the consequences of the expression of Eq. (100). For the velocity moment, we do not consider the redshift-space distortions. In the bottom rule of Fig. 4, we apply an approximation that the LPT kernel is independent of the time; i.e., perturbations of $n$th-order $\bm{\varPsi}^{(n)}$ are approximately proportional to $D^{n}$. If this approximation is not assumed, $nHf\bar{\bm{L}}_{n}$ should be replaced by $nHf\bar{\bm{L}}_{n}+\partial\bar{\bm{L}}_{n}/\partial t$.

The vertex resummation for the velocity momentum can be derived in the same way as that for the density field in Ref. Mat11 . The diagrammatic rule for the resummed vertex is given in Fig. 5.

This rule is an extension of that for the density field given in Fig. 15 of Ref. Mat11 . An important difference of the momentum field from density field is the existence of the dotted line. Only one dotted line should be attached to the resummed vertex. The grey ellipse with wavy lines represents the summation of every kind of diagrams which is attached to the external vertex by any number of wavy lines. This piece of diagram can be either connected or disconnected. The same is true for the grey ellipse with solid lines. Any diagram which is attached to the external vertex by both solid and wavy lines is not included in the resummation vertex (see Ref. Mat11 for details).

A further resummation of the vertex in the presence of a flat constraint is given in Fig. 6.

This resummation is an extension of Fig. 1 for the number density field. The rules are almost the same as those for the density field, except that they carry the spatial index $j$ for the momentum $\bm{j}_{X}$ which should be transmitted to the internal vertex of the bottom graph in Fig. 4. We also define the double square vertex without a dotted line by Fig. 7.

This diagram can be calculated by the rules of Figs. 6 and 4. As a result, the factor of $U^{(n)}_{X,j}$ in Fig. 7 is given by

[TABLE]

The first term is proportional to the derivative $\partial T^{(n)}_{X}/\partial k_{j}$, because we have

[TABLE]

which is straightforwardly shown by substituting the expression of Eq. (77). The second term can be calculated by noting

[TABLE]

As a result, Eq. (101) reduces to

[TABLE]

The functions $U^{(n)}_{X,j}$ commonly have an exponential factor, and we define normalized functions $\hat{U}^{(n)}_{X,j}$ by

[TABLE]

Substituting Eqs. (80) and (105) into Eq. (104), we have

[TABLE]

For $n=0,1,2$, Eqs. (81)–(83) and (106) give

[TABLE]

Using a somewhat lengthy preparation above, propagators of momentum field with the flat constraint are straightforwardly obtained by almost the same method as that in the previous subsection. The $(n+1)$-point propagator of the momentum field, $\bm{\varGamma}_{X}^{\mathrm{v}(n)}(\bm{k}_{1},\ldots,\bm{k}_{n})$, is defined by

[TABLE]

The diagrams for the two-point propagator for the momentum field are given in Fig. 8.

Applying diagrammatic rules of Figs. 6 and 7, and substituting Eqs. (81), (82), (107), and (108), we obtain a result,

[TABLE]

where “$+$ [higher-order terms]” indicates the terms which are proportional to $\sigma_{\mathrm{dpk}}^{2}P_{\mathrm{L}}$. These terms correspond to two-loop corrections in the usual sense, and we drop them in the following. Neglecting mode-coupling terms of Eq. (111), the two-point propagator reduces to

[TABLE]

This result agrees with a previous result, Eq. (11) of Ref. Bal15 in the peaks model ($\sigma^{2}_{\mathrm{d,pk}}$ in their notation is equal to $\sigma^{2}_{\mathrm{dpk}}/3$ in this paper).

The diagrams for the three-point propagator for the momentum field are given in Fig. 9.

Applying the diagrammatic rules as before, the result is given by

[TABLE]

where $\bm{k}=\bm{k}_{1}+\bm{k}_{2}$ and “$+$ [higher-order terms]” indicates the terms proportional to $\sigma_{\mathrm{dpk}}^{2}$, which corresponds to one-loop corrections in the usual sense, and we drop them in the following.

We define the power spectrum of number-weighted velocity field $P^{\mathrm{v}}_{X,ij}(\bm{k})$ of biased object $X$ by

[TABLE]

where $j_{X,i}$ is the $i$ component of $\bm{j}_{X}$. Using the propagators derived above, the number-weighted velocity power spectrum of biased object $X$ is given by

[TABLE]

where Gaussian initial conditions are assumed.

IV.5

Propagators of displacement field with the flat constraint

It is also interesting to find statistics of displacement field $\bm{\varPsi}$ for the biased objects. The number-weighted displacement field of the biased objects $X$ in Eulerian space can be defined by

[TABLE]

where $\bm{\varPsi}^{\mathrm{E}}(\bm{x})$ is the displacement field in Eulerian space. The latter is defined by the corresponding value of the displacement field in Lagrangian space, i.e., $\bm{\varPsi}^{\mathrm{E}}[\bm{x}=\bm{q}+\bm{\varPsi}(\bm{q})]=\bm{\varPsi}(\bm{q})$. The number-weighted displacement field in Lagrangian space is defined by

[TABLE]

The above displacement fields are related by

[TABLE]

The formalism in the previous subsection to calculate the propagators of momentum field can be similarly applied to the case of the number-weighted displacement field. The only difference is the diagrammatic rule in the bottom graph of Fig. 4: the corresponding factor should be replaced as $inaHfL_{n,j}\rightarrow iL_{n,j}$. This replacement is explained by the fact that the Fourier transform of the number-weighted displacement field in Eulerian space, $\tilde{\bm{\psi}}^{\mathrm{E}}_{X}$, is just given by replacing $a\dot{\tilde{\bm{\varPsi}}}\rightarrow\tilde{\bm{\varPsi}}$ in Eq. (100). Accordingly, the rules of Figs. 5 and 6 do not change, and the factor $aHf$ in the rule of Fig. 7 is just removed.

The $(n+1)$-point propagator of the number-weighted displacement field, $\bm{\varGamma}_{X}^{\mathrm{d}(n)}(\bm{k}_{1},\ldots,\bm{k}_{n})$, is defined by

[TABLE]

Following the same procedures in the previous subsection, and applying the above changes, the two-point propagator is given by

[TABLE]

The three-point propagator is given by

[TABLE]

where $\bm{k}=\bm{k}_{1}+\bm{k}_{2}$. The number-weighted displacement power spectrum is given by a similar expression of Eq. (115) with replacements, $P^{\mathrm{v}}_{X,ij}(k)\rightarrow P^{\mathrm{d}}_{X,ij}(k)$ and $\varGamma^{\mathrm{v}(n)}_{X,i}(k)\rightarrow\varGamma^{\mathrm{dv}(n)}_{X,i}(k)$.

V

Conclusion

In this paper, we investigate the effect of velocity bias in the formalism of iPT. The iPT is a formalism which can evaluate dynamical evolution of biased objects, based on the Lagrangian perturbation theory. The peaks model predicts the existence of the velocity bias. Even though the velocities of peaks and matter are assumed to be the same at peak locations, the velocities of peaks are statistically biased with respect to those of matter. With the formalism of iPT, the dynamical evolution of a biased object and its velocity momentum can be evaluated by higher-order Lagrangian perturbation theory beyond the Zel’dovich approximation.

In the first half of this paper, we see how the linear velocity bias emerges in the framework of iPT. The effects of velocity bias are already present in the formalism of iPT for the peaks model. Two- and three-point propagators for the velocity of peaks are calculated. As an example of their applications, we derive a formula for the one-loop approximation of velocity dispersion of peaks, which is also biased with respect to that of matter.

In the second half of this paper, a formal development of iPT in the presence of velocity bias are presented. We show that the emergence of velocity bias is a consequence of the flat constraint in general. Assuming that the flat constraint is the only element with odd parity, we generally derive the two- and three-point propagators of density in the one-loop approximation, with a resummation technique regarding the flat constraint. Diagrammtic rules for the propagators of momentum are generally derived, and a resummation technique regarding the flat constraint is also applied.

The formal development presented in this paper gives the basis of future applications of the iPT formalism. The one-loop power spectra from the two- and three-point propagators are given by Eqs. (87) and (115), and numerical evaluations of them can be performed by a similar method given in Ref. Mat14 . The purpose of this paper is to provide a theoretical understanding of the velocity bias with the formalism of iPT. Numerical evaluations of the power spectra, and comparisons with numerical simulations, etc., will be addressed in future work.

Acknowledgements.

The author gratefully thanks Vincent Desjacques and Tobias Baldauf for valuable discussions. This work is supported by JSPS KAKENHI Grants No. JP16H03977 and No. JP19K03835.

Appendix A

Evaluation of the function $T^{(n)}_{X}$

In this appendix, explicit expressions of $T_{X}^{(n)}(\bm{k};\bm{k}_{1},\ldots,\bm{k}_{n})$, defined by Eq. (77), are evaluated. With our assumption on the flat constraint, this function reduces to the form of Eq. (79), i.e.,

[TABLE]

To evaluate the above equation, we introduce a vector, $\bm{J}={\sigma_{0}}^{2}\bm{k}/3\sigma_{1}$. Using Eqs. (74) and (76), the operator which appears in Eq. (122) reduces to

[TABLE]

Substituting Eqs. (62) and (123) into (122), we have

[TABLE]

With our assumption, the variables $\bm{\eta}$ are independent of the other variables $\nu,\bm{\zeta},\ldots$, and thus the average over $\bm{\eta}$ is independently evaluated in the above equation. The probability distribution function of $\bm{\eta}$ is given by $P(\bm{\eta})=(3/2\pi)^{3/2}e^{-3\eta^{2}/2}$. Expressing the delta function by a Fourier integral, and applying the multidimensional formula of Gaussian integration, we derive

[TABLE]

Temporarily putting $\bm{J}=0$ in Eq. (125), Eq. (62) indicates

[TABLE]

Thus we have

[TABLE]

From the above expression, the functions $T_{X}^{(n)}$ have a common factor of $\exp(3J^{2}/2)=\exp(k^{2}{\sigma_{0}}^{4}/6{\sigma_{1}}^{2})$ for every $n$. We thus define normalized functions $\hat{T}_{X}^{(n)}$ by

[TABLE]

One can evaluate Eq. (127) as follows. For $n=0$, we have $T_{X}^{(0)}(\bm{k})=e^{3J^{2}/2}$, and therefore

[TABLE]

For $n=1$, we have

[TABLE]

When $\bm{k}=\bm{0}$ and $\bm{J}=\bm{0}$, we have $T_{X}^{(1)}(\bm{0};\bm{k}_{1})=c_{X}^{(1)}(\bm{k}_{1})$ from the definition, Eq. (77), and the above equation in this case reduces to

[TABLE]

Substituting this equation into Eq. (130), the normalized function is given by

[TABLE]

For $n=2$, we have

[TABLE]

When $\bm{k}=\bm{0}$ and $\bm{J}=\bm{0}$, we have $T_{X}^{(2)}(\bm{0};\bm{k}_{1},\bm{k}_{2})=c_{X}^{(2)}(\bm{k}_{1},\bm{k}_{2})$ from the definition, Eq. (77), and the above equation in this case reduces to

[TABLE]

Substituting the above equation and Eq. (131) into Eq. (133), the normalized function is given by

[TABLE]

The expressions of the function $T_{X}^{(n)}$ for $n\geq 3$ can be derived by following the similar procedures above.