Precision cosmology with redshift-space bispectrum: a perturbation theory based model at one-loop order

TL;DR

This paper develops a perturbation theory-based model for the redshift-space matter bispectrum that accurately captures non-perturbative damping effects, improving cosmological parameter constraints from large-scale structure data.

Contribution

It extends previous models by incorporating a univariate damping function at one-loop order, validated against N-body simulations for better accuracy in bispectrum predictions.

Findings

Model matches N-body simulations for monopole and quadrupole moments

Univariate damping function effectively controls non-perturbative effects

Next-to-leading order calculations yield consistent damping parameters

Abstract

The large-scale matter distribution in the late-time Universe exhibits gravity-induced non-Gaussianity, and the bispectrum, three-point cumulant is expected to contain significant cosmological information. In particular, the measurement of the bispectrum helps to tighten the constraints on dark energy and modified gravity through the redshift-space distortions (RSD). In this paper, extending the work by Taruya, Nishimichi & Saito (2010, Phys.Rev.D 82, 063522), we present a perturbation theory (PT) based model of redshift-space matter bispectrum that can keep the non-perturbative damping effect under control. Characterizing this non-perturbative damping by a univariate function with single free parameter, the PT model of the redshift-space bispectrum is tested against a large set of cosmological $N$-body simulations, finding that the predicted monopole and quadrupole moments are in a good agreement with simulations at the scales of baryon acoustic oscillations (well beyond the range of agreement of standard PT). The validity of the univariate ansatz of the damping effect is also examined, and with the PT calculation at next-to-leading order, the fitted values of the free parameter is shown to consistently match those obtained from the PT model of power spectrum by Taruya, Nishimichi & Saito (2010).

Tables

Table 1. Table 1: Parameters of N 𝑁 N -body simulation data sets used in the paper.

	HR	LR
Set	B	C
Number of realizations	$96$	$512$
Box size [$h^{-1}$$\mathrm{Mpc}$]	$656.25$	$1312.5$
Number of particles	${1024}^3$	${512}^3$
Particle mass [$h^{-1}$$M_{\odot }$]	$1.88\times {10}^{10}$	$1.20\times {10}^{12}$
Initial redshift ($z_{\mathrm{ini}}$)	$190$	$106$
Output redshifts	$z=1,0$	$z=1,0.5,0$

Table 2. Table 2: Fitted values of σ v subscript 𝜎 v \sigma_{\rm v} given in Eq. ( 54 ) at various redshifts and triangular shapes. The results are obtained by fitting the monopole and quadrupole predictions of PT models to the measured results from the LR (HR) data of N 𝑁 N -body simulations, and are listed in units of h ℎ h\,  Mpc -1 . For reference, the rightmost column shows the linear theory predictions.

	$k_{\max }\left[h{\mathrm{Mpc}}^{-1}\right]$		$P_{\mathrm{model}}^{\left(s\right)}$: 1-loop	$B_{\mathrm{model}}^{\left(s\right)}$: 1-loop		$B_{\mathrm{model}}^{\left(s\right)}$: tree		${\sigma }_{v,\mathrm{lin}}$
	1-loop	tree		$k_1=k_2=k_3$	$k_1=2k_2=\left(2/\sqrt{3}\right)k_3$	$k_1=k_2=k_3$	$k_1=2k_2=\left(2/\sqrt{3}\right)k_3$
$z=1$	0.13	0.09	$4.3\pm 0.1$ ($4.2\pm 0.04$)	$4.9\pm 0.3$ ($5.9\pm 0.3$)	$4.0\pm 0.2$ ($4.3\pm 0.3$)	$6.1\pm 0.2$	$6.5\pm 0.1$	$3.8$
$z=0.5$	0.1	0.07	$5.9\pm 0.1$	$7.8\pm 0.4$	$6.3\pm 0.3$	$9.0\pm 0.3$	$9.5\pm 0.2$	$4.8$
$z=0$	0.08	0.06	$8.8\pm 0.3$ ($9.8\pm 0.3$)	$13\pm 0.8$ ($15\pm 1.0$)	$10\pm 0.6$ ($9.1\pm 1.1$)	$13\pm 0.6$	$15\pm 0.4$	$6.1$

Equations

$$\mathbf{s}=\mathbf{r}+\frac{(1+z)v_z(\mathbf{r})}{H(z)}\hat{\mathbf{z}},$$

$${\delta }^{(\mathrm{s})}(\mathbf{s})={\left|\frac{\partial \mathbf{s}}{\partial \mathbf{r}}\right|}^{-1}\{1+\delta (\mathbf{r})\}-1.$$

$${\delta }^{(\mathrm{s})}(\mathbf{k})=\int d{\mathbf{r}}^3\left\{\delta (\mathbf{r})-\frac{{\nabla }_z{v}_z(\mathbf{r})}{aH(z)}\right\}{e}^{i(k\mu v_{\mathrm{z}}/H+\mathbf{k}\cdot \mathbf{r})},$$

$$(2\pi {)}^3{\delta }_{\mathrm{D}}(\mathbf{k}+{\mathbf{k}}^{\prime })P^{(\mathrm{s})}(\mathbf{k})\equiv ⟨{\delta }^{(\mathrm{s})}(k){\delta }^{(\mathrm{s})}(k^{\prime })⟩,$$

$$(2\pi {)}^3{\delta }_{\mathrm{D}}({\mathbf{k}}_1+{\mathbf{k}}_2+{\mathbf{k}}_3)B^{(\mathrm{s})}({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3)$$

$$\equiv ⟨{\delta }^{(\mathrm{s})}(k_1){\delta }^{(\mathrm{s})}(k_2){\delta }^{(\mathrm{s})}(k_3)⟩,$$

$$P^{(\mathrm{s})}(\mathbf{k})=\int d\mathbf{r}{e}^{i\mathbf{k}\cdot {\mathbf{r}}_{12}}$$

$$\displaystyle\times\biggl{\langle}e^{-ifk_{z}u_{12z}}\left\{\delta(\bm{r}_{1})+f\nabla_{z}u_{z}(\bm{r}_{1})\right\}\left\{\delta(\bm{r}_{2})+f\nabla_{z}u_{z}(\bm{r}_{2})\right\}\biggr{\rangle},$$

$$B^{(\mathrm{s})}({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3)=\int d{\mathbf{r}}_{13}{\mathbf{r}}_{23}{e}^{i({\mathbf{k}}_1\cdot {\mathbf{r}}_{13}+{\mathbf{k}}_2\cdot {\mathbf{r}}_{23})}$$

$$\displaystyle\times\biggl{\langle}e^{-if(k_{1z}u_{13z}+k_{2z}u_{23z})}\left\{\delta(\bm{r}_{1})+f\nabla_{z}u_{z}(\bm{r}_{1})\right\}$$

$$\displaystyle\times\left\{\delta(\bm{r}_{2})+f\nabla_{z}u_{z}(\bm{r}_{2})\right\}\left\{\delta(\bm{r}_{3})+f\nabla_{z}u_{z}(\bm{r}_{3})\right\}\biggr{\rangle},$$

$${\delta }^{(\mathrm{s})}(\mathbf{k})=\sum _{n=0}\int \frac{d^3{\mathbf{q}}_1}{(2\pi {)}^3}\cdots \int \frac{d^3{\mathbf{q}}_n}{(2\pi {)}^3}{\delta }_{\mathrm{D}}(\mathbf{k}-{\mathbf{q}}_{1\cdots n})$$

$$\displaystyle\times\bigl{\{}\delta(\bm{k})+f\mu^{2}\theta(\bm{k})\bigr{\}}\,\frac{(f\mu k)^{n}}{n!}\frac{\mu_{1}}{q_{1}}\theta(\bm{q}_{1})\cdots\frac{\mu_{n}}{q_{n}}\theta(\bm{q}_{n})$$

$$\delta (\mathbf{k})=\sum _{n=1}\int \frac{d^3{\mathbf{p}}_1\cdots d^3{\mathbf{p}}_n}{(2\pi {)}^{3n}}{\delta }_{\mathrm{D}}(\mathbf{k}-{\mathbf{p}}_{1\cdots n})F_n({\mathbf{p}}_1,\cdots ,{\mathbf{p}}_n)$$

$$\times {\delta }_{\mathrm{L}}({\mathbf{p}}_1)\cdots {\delta }_{\mathrm{L}}({\mathbf{p}}_n),$$

$$\theta (\mathbf{k})=\sum _{n=1}\int \frac{d^3{\mathbf{p}}_1\cdots d^3{\mathbf{p}}_n}{(2\pi {)}^{3n}}{\delta }_{\mathrm{D}}(\mathbf{k}-{\mathbf{p}}_{1\cdots n})G_n({\mathbf{p}}_1,\cdots ,{\mathbf{p}}_n)$$

$$\times {\delta }_{\mathrm{L}}({\mathbf{p}}_1)\cdots {\delta }_{\mathrm{L}}({\mathbf{p}}_n),$$

$${\delta }^{(\mathrm{s})}(\mathbf{k})=\sum _{n=1}\int \frac{d^3{\mathbf{q}}_1}{(2\pi {)}^3}\cdots \int \frac{d^3{\mathbf{q}}_n}{(2\pi {)}^3}{\delta }_D(\mathbf{k}-{\mathbf{q}}_{1\cdots n})$$

$$\times Z_n({\mathbf{q}}_1,\cdots ,{\mathbf{q}}_n){\delta }_{\mathrm{L}}({\mathbf{q}}_1)\cdots {\delta }_{\mathrm{L}}({\mathbf{q}}_n).$$

$$P^{(\mathrm{s})}(\mathbf{k})$$

$$B^{(\mathrm{s})}({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3)$$

$$+B_{411}^{1\text{mbox}-\mathrm{loop}}+\cdots .$$

$$P^{\mathrm{lin}}(\mathbf{k})$$

$$P^{1\text{mbox}-\mathrm{loop}}(\mathbf{k})$$

$$B^{\mathrm{tree}}({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3)$$

$$B_{222}^{1\text{mbox}-\mathrm{loop}}\left({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3\right)$$

$$+2\mathrm{perms}({\mathbf{k}}_1\leftrightarrow {\mathbf{k}}_2\leftrightarrow {\mathbf{k}}_3),$$

$$B_{321\text{mbox}-\mathrm{I}}^{1\text{mbox}-\mathrm{loop}}\left({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3\right)$$

$$+5\mathrm{perms}({\mathbf{k}}_1\leftrightarrow {\mathbf{k}}_2\leftrightarrow {\mathbf{k}}_3),$$

$$B_{321\text{mbox}-\mathrm{II}}^{1\text{mbox}-\mathrm{loop}}\left({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3\right)$$

$$B_{411}^{1\text{mbox}-\mathrm{loop}}\left({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3\right)$$

$$B^{(\mathrm{s})}({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3)$$

$$A_1$$

$$A_2$$

$$A_3$$

$$A_4$$

$$A_5$$

$$j_4$$

$$j_5$$

$$⟨e^{\mathbf{j}\cdot \mathbf{A}}⟩=\exp \left\{⟨e^{\mathbf{j}\cdot \mathbf{A}}{⟩}_{\mathrm{c}}\right\},$$

$$⟨A_1{A}_2{A}_3{e}^{j_4{A}_4+j_5{A}_5}⟩$$

$$+⟨A_2{e}^{j_4{A}_4+j_5{A}_5}{⟩}_{\mathrm{c}}⟨A_1{A}_3{e}^{j_4{A}_4+j_5{A}_5}{⟩}_{\mathrm{c}}+⟨A_1{e}^{j_4{A}_4+j_5{A}_5}{⟩}_{\mathrm{c}}⟨A_2{A}_3{e}^{j_4{A}_4+j_5{A}_5}{⟩}_{\mathrm{c}}$$

$$\displaystyle+\langle A_{1}e^{j_{4}A_{4}+j_{5}A_{5}}\rangle_{\rm c}\langle A_{2}e^{j_{4}A_{4}+j_{5}A_{5}}\rangle_{\rm c}\langle A_{3}e^{j_{4}A_{4}+j_{5}A_{5}}\rangle_{\rm c}\biggr{]}.$$

$$B^{(\mathrm{s})}({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3)⟶D_{\mathrm{FoG}}(k_1{\mu }_1,k_2{\mu }_2,k_3{\mu }_3)\int d{\mathbf{r}}_{13}d{\mathbf{r}}_{23}$$

Full text

cosmology, large-scale structure

††preprint: YITP-17-46

I Introduction

The galaxy clustering observed via the spectroscopic measurements appears distorted due to the peculiar velocity of galaxies. This systematic effect is known as the redshift-space distortions (RSD) (e.g., Refs. Davis and Peebles (1983); Kaiser (1987); Hamilton (1992); Scoccimarro (2004)), which manifestly break the statistical isotropy. While the RSD complicates the interpretation of the galaxy clustering data, it now attracts much attention as a probe of gravity on cosmological scales. This is because, on large scales, the effect is accounted by the linear theory, and the size of the anisotropies is related to the growth of structure (e.g., Refs. Kaiser (1987); Hamilton (1992)). To be precise, the measurement of RSD offers an opportunity to determine the parameter $f\,\sigma_{8}$ (e.g., Ref. Percival and White (2009); Song and Percival (2009); Yamamoto et al. (2008)). Here, $f$ is the linear growth rate defined by $f\equiv\,d\ln D_{+}/d\ln a$ with $D_{+}$ and $a$ being respectively the linear growth factor and scale factor of the Universe, and $\sigma_{8}$ is the normalization amplitude of the linear power spectrum smoothed over $8\,h^{-1}$ Mpc. Combining the RSD with the measurement of baryon acoustic oscillations (BAO) imprinted on galaxy clustering (e.g., Eisenstein and Hu (1998); Meiksin et al. (1998); Eisenstein et al. (2005)), we can simultaneously estimate both the geometric distances and growth of structure (e.g., Taruya et al. (2010); Reid et al. (2012); Oka et al. (2014); Beutler et al. (2013)), from which we are able to scrutinize the theory of gravity that accounts for cosmic expansion and structure formation (e.g., Refs. Linder (2008); Seo and Eisenstein (2003); Blake and Glazebrook (2003); Glazebrook and Blake (2005); Song et al. (2015)).

The measurements of RSD taken from various observations now cover a wide range of redshifts out to $z\sim 1.4$ (e.g., Refs. Guzzo et al. (2008); Okumura et al. (2016); de la Torre et al. (2013); Beutler et al. (2016); Blake et al. (2011), see also Sec. 5.5 of Ref. Planck Collaboration et al. (2016)), and the results are broadly consistent with general relativity. But a closer look at the constrained values of $f\,\sigma_{8}$ suggests a mild tension with the Planck $\Lambda$CDM model Planck Collaboration et al. (2015), indicating some systematics or potentially new physics. A further improvement on the RSD measurement is thus important, mitigating both the statistical errors and systematics. This is one of the main reasons why there are various projects aiming at precisely measuring RSD which will uncover a large cosmic volume.

With a gigantic survey volume, the next-generation galaxy redshift surveys (e.g. Euclid, LSST) will therefore offer an opportunity to precisely measure the power spectrum or correlation function at an unprecedented precision, and provided a high-precision template, a robust and tight constraint on $f\,\sigma_{8}$ will be expected. Furthermore, one promising point is that taking advantage of a large-volume data, a high-statistical measurement of the higher-order statistics such as bispectrum or three-point correlation function is made available, and combining it with power spectrum measurement further helps to tighten the constraint on $f\,\sigma_{8}$.

It is widely known that the nonlinearity of gravity generically incorporates the non-Gaussian statistical nature into the large-scale structures, and a part of the cosmological information initially encoded in the two-point statistics is leaked to the higher-order statistics. In this respect, the combination of the bispectrum data with power spectrum in the cosmological data analysis is a natural first step to efficiently extract the cosmological information from the large-scale structures. There have been various studies based on the Fisher matrix analysis to quantify the statistical impact of the bispectrum on the cosmological parameter estimation, showing that adding the bispectrum indeed plays an important role to break parameter degeneracy (e.g., Refs. Sefusatti et al. (2006); Kayo et al. (2013); Sato and Nishimichi (2013); Greig et al. (2013)). In particular, the impact of combining bispectrum measurement is demonstrated for the RSD measurement by Refs. Song et al. (2015); Gagrani and Samushia (2017), showing that even if we restrict the data to the large-scale modes, the constraint on $f\,\sigma_{8}$ will be improved by more than a factor of two.

However, most of these analysis has been demonstrated in a very simplified setup, assuming that the observed bispectrum is described by the perturbation theory (PT) at leading order Scoccimarro et al. (1999), on top of a simple prescription for galaxy bias. While such a simplified forecast study helps to understand the ability and/or potential of the planned galaxy surveys, there are a number of remarks and cautions in a practical data analysis. One important point to be noted is the theoretical template of redshift-space bispectrum. While the perturbative description is supposed to be adequate in real space at the weakly nonlinear stage of gravitational clustering, the situation becomes more subtle in redshift space, because the redshift space is related to the real space through nonlinear mapping. As a result, in terms of the real-space quantities, even the redshift-space power spectrum cannot be simply expressed as the large-scale two-point statistics of the underlying fields, and is significantly affected by the small-scale physics (e.g., Ref. Scoccimarro (2004)). Hence, a sophisticated treatment is needed for a reliable theoretical template even in the weakly nonlinear regime Taruya et al. (2010); Reid and White (2011); Vlah et al. (2012); Matsubara (2008a, b); Carlson et al. (2013); Matsubara (2014); Taruya et al. (2013).

The aim of this paper is to address this issue in the case of bispectrum, and based on the next-to-leading order calculations, we present a perturbation-theory model of redshift-space bispectrum, which consistently incorporates the effect coming from the small-scale virial motion. While modeling the redshift-space bispectrum has been already investigated based on fitting formula Gil-Marín et al. (2012); Scoccimarro and Couchman (2001) and halo model (e.g., Refs. Smith et al. (2008); Yamamoto et al. (2016)), analytical treatment would be certainly powerful in characterizing the anisotropies of the redshift-space bispectrum (see e.g., Rampf and Wong (2012) for discussion on resummed PT treatment). In particular, with perturbation theory, one can give an accurate description for the large scales of our interest. Albeit its limitation, the PT-based modeling therefore deserves further consideration, and we present, for the first time, monopole and quadrupole moments of the redshift-space bispectrum at one-loop order, which agree well with $N$-body simulations at BAO scales.

This paper is organized as follows. In Sec. II, we begin by briefly reviewing the relation between real and redshift spaces, and derive the exact formulas for power spectrum and bispectrum in redshift space. Based on the standard PT, the expressions for redshift-space bispectrum at one-loop order are also presented. Then, in Sec. III, extending the work by Ref. Taruya et al. (2010), we give the PT-based model of redshift-space bispectrum, which consistently includes both the next-to-leading order PT corrections and non-perturbative damping effect coming from the small-scale virial motion. Rewriting the exact formulas with the cumulants, we derive the expressions for our PT model valid at one-loop order. In Sec. IV, the proposed PT model is quantitatively tested against a large suite of cosmological $N$-body simulations, and the validity of the assumption is checked in detail. Finally, Sec. V is devoted to discussion and conclusion on our important findings. Throughout the paper, we assume $\Lambda$CDM model and adopt the cosmological parameters determined from seven-year WMAP results Komatsu et al. (2011): $h=0.72$, $\Omega_{\rm m}=0.2573$, $\Omega_{b}=0.04356$, $\sigma_{8}=0.801$ and $n_{s}=0.963$.

II Power spectrum and Bispectrum in redshift space

In this section, we begin by writing down the exact expressions of power spectrum and bispectrum for matter density field in redshift space. We then consider the standard PT calculations, and based on the explicit expressions for bispectrum up to the next-to-leading order called one-loop, we discuss its limitation and applicability.

II.1 Exact formula of power spectrum and bispectrum

Throughout the paper, we consider the distant-observer case, in which observer’s line-of-sight direction is described by a specific direction given by the unit vector $\hat{\bm{z}}$. The mapping between real-space position $\bm{r}$ and redshift-space position $\bm{s}$ is then given by

[TABLE]

where $z$ is the redshift, $v_{z}$ is the line-of-sight component of the velocity field, $\bm{v}$ (i.e., $v_{z}=\bm{v}\cdot\bm{z}$). The quantity $H$ is the Hubble constant. Denoting respectively the density fields in real and redshift spaces by $\delta(\bm{r})$ and $\delta^{({\rm s})}(\bm{s})$, the mass conservation implies $\{1+\delta^{({\rm s})}(\bm{s})\}d^{3}\bm{s}=\{1+\delta(\bm{r})\}d^{3}\bm{r}$, which leads to

[TABLE]

Since we are interested in the statistical quantities in Fourier-space, we rewrite the above relation in terms of the Fourier counterpart:

[TABLE]

where the variable $\mu$ is defined as the directional cosine of the angle between the line-of-sight $\hat{\bm{z}}$ and wave vector $\bm{k}$.

Using Eq. (3), the exact formulas for the power spectrum and bispectrum in redshift space, expressed in terms of the real-space quantities, are derived. Recalling that the power spectrum and bispectrum are respectively defined as

[TABLE]

substituting Eq. (3) into the above leads to

[TABLE]

for the redshift-space power spectrum, and

[TABLE]

for the redshift-space bispectrum. Here, we define the normalized peculiar velocity $\bm{u}(\bm{r})=-\bm{v}(\bm{r})/\{faH\,(z)\}$ and denote the pairwise normalized velocity of the separation $\bm{r}_{ij}=\bm{r}_{i}-\bm{r}_{j}$ by $\bm{u}_{ij}\equiv\bm{u}_{i}-\bm{u}_{j}$. The function $f$ is the linear growth rate defined by $f\equiv d\ln{D(z)}/d\ln{a}$, with $D(z)$ being the linear growth factor.

The above expressions show that albeit the simple relation (1), the power spectrum and bispectrum in redshift space are rather intricate statistical relation. Qualitatively, the amplitude of the power spectrum is known to be boosted by the additional term of the velocity field at large scales (Kaiser effect), and to be suppressed at small scales by an exponential damping factor (Fingers-of-God effect). Since the structure of the expressions is of the same form in both the power spectrum and bispectrum, we expect that the redshift-space bispectrum possesses similar qualitative features. The additional complexity in bispectrum is, however, that it is no longer characterized simply by the shape of the triangle, i.e., length of three wave vectors $k_{1}$, $k_{2}$ and $k_{3}$, or length of vectors $k_{1}$ and $k_{2}$ and their angle $\theta_{12}\equiv\cos^{-1}(\bm{\hat{k}}_{1}\cdot\bm{\hat{k}}_{2})$. On top of these three variables, one needs two more variables to describe the orientation of the triangular shape with respect to the line-of-sight direction. In this respect, the identification and separation of the Kaiser and Finger-of-God effects are rather non-trivial, and a careful treatment is required for an accurate modeling of redshift-space bispectrum.

II.2 Standard perturbation theory in redshift space

In this paper, we employ the perturbation theory technique to construct an analytic model of redshift-space bispectrum relevant at large scales. To do this, we here present the standard PT results at next-to-leading order called one-loop, and discuss its limitation and the issues to be improved. We will then present, in next section, a PT model of redshift-space bispectrum which properly incorporates the effect of FoG damping in a non-perturbative manner.

To derive the standard PT results of redshift-space bispectrum, let us first expand the redshift-space density field in Fourier space. Regarding the real-space density and velocity fields as perturbative quantities, a systematic expansion of Eq. (3) leads to

[TABLE]

with $\bm{q}_{1\cdots n}$ being $\bm{q}_{1}+\cdots+\bm{q}_{n}$. In the above, we assumed the irrotational flow of the velocity field, and introduced the velocity-divergence field defined by $\theta(\bm{x})=-\nabla\cdot\bm{v}/(f\,aH)$. Then, we further expand the real-space quantities in terms of the standard PT kernels:

[TABLE]

where $\delta_{\rm L}$ is the linear density field, and $F_{n}$ and $G_{n}$ are the standard PT kernels (see for example Bernardeau et al. (2002)). Plugging the above expressions into Eq. (8), reorganizing the perturbative expansion in powers of $\delta_{\rm L}$ leads to

[TABLE]

Here, the $Z_{n}$ are the so-called redshift-space PT kernel, and these are expressed in terms of the real-space PT kernels, $F_{n}$ and $G_{n}$ in a rather intricate way. The explicit expression for $Z_{n}$ is presented in Appendix A. With the PT expression of the density field in Eq. (11), the redshift-space power spectrum and bispectrum are expanded as

[TABLE]

Each term at the right-hand-side is given by

[TABLE]

for the power spectrum, and

[TABLE]

for the bispectrum. In deriving the equations above, we assume the Gaussianity of linear density field $\delta_{\rm L}$, and denote its power spectrum by $P_{\rm L}$. These expressions can be reduced to the real-space power spectrum and bispectrum of density field if we replace the kernels $Z_{n}$ with $F_{n}$. A notable point is that the redshift-space kernel implicitly depends on the line-of-sight direction. As a result, the statistical isotropy is broken in each term of power spectrum and bispectrum, and we need one and two more additional variables to characterize the redshift-space power spectrum and bispectrum, respectively. Another important point is that the kernels $Z_{n}$ at $n\geq 2$ includes the mode-coupling contributions from the velocity fields, which basically come from the perturbative expansion of the exponential factor in Eq. (3). This implies that the Finger-of-God damping effect cannot be reproduced in a naive standard PT treatment, and we need to resum the infinite series of PT expansions. As it has been shown in Ref. Taruya et al. (2010), the standard PT prediction in redshift space largely overestimates the power spectrum amplitude at one-loop order, and cannot accurately describe the BAO in redshift space. The applicable range of one-loop prediction thus becomes narrower than that of the real-space prediction. Since the expression of bispectrum also uses the redshift PT kernels and it even includes the higher-order, the situation must be similar or rather worse than that in the power spectrum case. We will see in Sec. IV that without accounting the Finger-of-God damping, standard PT prediction of bispectrum starts to deviate from simulation even at large scales. A proper account of the damping effect is essential, and we will consider how to incorporates the damping effect into the PT calculation in next section.

III An improved model prescription for bispectrum at one-loop order

In this section, we present the PT model of redshift-space bispectrum which keeps the non-perturbative damping effect. Our strategy is to decompose the contributions into non-perturbative part and the terms which can be evaluated with PT calculation, starting with the exact expression, Eq. (7). For this purpose, we follow the treatment by Ref. Taruya et al. (2010), and rewrite the exact expression in terms of cumulants. We then identify the non-perturbative part responsible for the FoG damping. Based on the simple proposition similarly made by Ref. Taruya et al. (2010), the non-perturbative damping term is separated out from the rest of the contributions, for which we apply the PT calculation. We derive the expression valid at one-loop order.

Let us begin by rewriting Eq. (7) in the form

[TABLE]

where the quantities $A_{i}$, $j_{i}$ are respectively defined by

[TABLE]

with $\mu_{i}=\bm{k}_{i}\cdot\hat{\bm{z}}/k_{i}$. To express the moment given above in terms of the cumulants, we use the relation between moment and cumulant generating functions (e.g., Taruya et al. (2010); Scoccimarro (2004); Matsubara (2008a)). For the stochastic vector field $\bm{A}$, we have

[TABLE]

with $j$ being an arbitrary constant vector, $\bm{j}$. To be specific, we assume that the vector fields given above are five components, i.e., $\bm{A}=\{A_{1},A_{2},A_{3},A_{4},A_{5}\}$ and $\bm{j}=\{j_{1},j_{2},j_{3},j_{4},j_{5}\}$. Then, taking the derivative three times with respect to $j_{1}$, $j_{2}$ and $j_{3}$, we set $j_{1}=j_{2}=j_{3}=0$. We obtain

[TABLE]

This equation is indeed what we want to derive, and the left hand side is exactly the same one as in the integrand of Eq. (21).

The above expression shows that the pairwise velocity fields $A_{4}$ and $A_{5}$ in the exponent produce non-trivial correlations with density and velocity gradient fields. Further, these fields appear in the overall prefactor. This is indeed the same structure as seen in the expression of power spectrum Taruya et al. (2010). In power spectrum, the exponential prefactor is known to give a suppression of the amplitude due to the large-scale coherent and small-scale virialized motions, and at the large-scale of our interest, it mainly affects the broadband shape of the power spectrum. We expect that the overall prefactor in Eq. (30) similarly behaves like the one in the power spectrum, and it can alter the broadband shape of the bispectrum. Because of its functional form, the prefactor is likely to be affected by the small-scale nonlinearity even at large scales, and we may thus take it as non-perturbative part. On the other hand, the terms in the square bracket of Eq. (30) include the density fields and are responsible for reproducing the real-space bispectrum in the absence of redshift-space distortions [this is simply obtained by setting all the velocity fields in Eq. (30) to zero]. Thus, these terms basically carry the cosmological information, and imprints the acoustic feature of BAO. Although each term in the square bracket contains the exponential factor $e^{j_{4}A_{4}+j_{5}A_{5}}$, the contribution can be small as long as we consider the BAO scales, and the perturbative expansion may work well.

Based on the considerations, we adopt the proposition by Ref. Taruya et al. (2010) to derive the PT model of redshift-space bispectrum valid at weakly nonlinear scales. That is,

(i)

The overall prefactor, $\exp\left\{\langle e^{j_{4}A_{4}+j_{5}A_{5}}\rangle_{\rm c}\right\}$, is kept as a non-perturbative contribution, and is replaced with general functional form $D_{\rm FoG}$, which is assumed to be given as a function of $k_{1}\mu_{1}$, $k_{2}\mu_{2}$ and $k_{3}\mu_{3}$, ignoring the spatial correlation of $A_{4}$ and $A_{5}$. This means that the zero-lag correlation is the only dominant contribution. The relevant functional form of $D_{\rm FoG}$ will be discussed in Sec. IV.

(ii)

The terms in the square bracket are treated perturbatively, and regarding the variables $j_{4}$ and $j_{5}$ as expansion parameters, we collect the contributions valid at one-loop order in standard PT.

From the proposition (i), the overall exponential factor is factorized outside the integral over $\bm{r}_{13}$ and $\bm{r}_{23}$. We have

[TABLE]

We then expand the terms in the square bracket. Up to the third order in $j_{n}$, we obtain

[TABLE]

where the functions $C_{n}$ are defined by

[TABLE]

with the integrands $S_{n}$ given below:

[TABLE]

In the above, the relevant terms for one-loop PT calculations, which are of the order of $\mathcal{O}(P_{\rm L}^{3})$, appear at $n\leq 6$, and the rest of the terms turns out to be higher-order. Hence, keeping the terms valid at the one-loop level, we model the redshift-space bispectrum as

[TABLE]

The above model is compared with the standard PT results in Sec. II.2, which can be also derived from Eqs. (21) and (30) by a naive expansion of all the exponential factors, $e^{i(j_{4}A_{4}+j_{5}A_{5})}$. Collecting the relevant contributions at one-loop order, we have

[TABLE]

with the functions $D_{n}$ given by

[TABLE]

where the function $S_{n}$ $(n=1,2)$ is defined by Eqs. (34) and (35). These terms come from the expansion of the overall prefactor $\exp\left\{\langle e^{j_{4}A_{4}+j_{5}A_{5}}\rangle_{\rm c}\right\}$. Hence, at one-loop order, Eq. (45) is recast as

[TABLE]

In what follows, we use Eq. (48) to compute the PT model of redshift-space bispectrum, and compare the predictions with $N$-body simulations. To be precise, we first compute $B^{\rm(s)}_{\rm SPT,1\mbox{-}loop}$ based on the standard PT calculations summarized in Sec. II.2 [Eqs. (13) with (16), (17), (18), (19), and (20)]. Then, we subtract $D_{1}$ and $D_{2}$ terms from the standard PT bispectrum and, we take into account Fingers-of-God effect. The explicit expressions for $D_{1}$ and $D_{2}$ relevant for the one-loop calculations are presented in Appendix B.

Before closing this section, we look at the significance of the difference between the standard PT bispectrum and the model given in Eq. (48) or (45). In Fig. 2, ignoring the $D_{\rm FoG}$ contribution, the monopole and quadrupole moments of the bispectrum are computed at $z=1$ for equilateral (top) and isosceles (bottom) configurations, and the results are plotted as function of $k$ and angle $\theta_{12}\equiv\cos^{-1}(\hat{\bm{k}}_{1}\cdot\hat{\bm{k}}_{2})$, respectively. Here, the multipole moments of the bispectrum, $B^{\rm(s)}_{\ell}$, are defined by:

[TABLE]

where the function $\mathcal{P}_{\ell}(\mu)$ is the Legendre polynomials with $\mu$ being the directional cosine given by $\mu=\cos\omega$. The angles $\omega$ and $\phi$ characterize the orientation of the triangles (i.e., $\bm{k}_{1}$, $\bm{k}_{2}$, and $\bm{k}_{3}$) with respect to the line-of-sight direction (see Fig. 1). Precise definition of the angles $\omega$ and $\phi$, together with the properties of multipole expansion, is described in Appendix C. Note that our definition differs from the one frequently used in the literature (e.g., Refs. Scoccimarro et al. (1999); Scoccimarro (2015); Yamamoto et al. (2016); Gagrani and Samushia (2017)), but a nice property is that the bispectrum multipoles become fully symmetric under the permutation of the order of $k_{1}$, $k_{2}$ and $k_{3}$.

Fig. 2 shows that both the one-loop bispectra, $B^{\rm(s)}_{\rm SPT,1\mbox{-}loop}$ (red) and $B_{\rm model}^{\rm(s)}$ (magenta), have a larger amplitude than the tree-level prediction. While these predictions lead to similar scale and shape dependencies, differences appear manifest at smaller scales due to the the $D_{1}$ and $D_{2}$ terms. Though these two terms are basically small and weakly depend on scales, we will see below that subtracting these from $B^{\rm(s)}_{\rm SPT,1\mbox{-}loop}$ indeed plays an important role to get a consistent damping behavior of $D_{\rm FoG}$ from both power spectrum and bispectrum (see Sec. IV.4).

IV Comparison with N-body simulations

In this section, we present a detailed comparison of the redshift-space bispectrum between PT predictions and $N$-body simulations. After briefly describing the $N$-body data set used in the analysis in Sec. IV.1, we first compare the real-space results in Sec IV.2 to see the applicable range of PT as well as the quality of $N$-body data. We then move to the redshift space, and compare the monopole and quadrupole bispectra from $N$-body simulations with those obtained by PT calculations in Sec. IV.3. Sec. IV.4 discusses the validity and consistency of the ansatz for the damping function $D_{\rm FoG}$ in our PT model of bispectrum.

IV.1 N-body simulations and measurement of the bispectrum

We use the simulation ensembles Set B and Set C from Dark Energy Universe Simulation - Parallel Universe Runs (DEUS-PUR) introduced by Blot et al. (2015) 111Set C was not directly presented in Blot et al. (2015) but it was performed at the same time for the same project. Each simulation started from a given realization of the initial matter density field. The initial conditions were generated with an improved version of MPGRAFIC Prunet et al. (2008) while the particles evolution were computed with an optimized version of the RAMSES N-body code Teyssier (2002). The main characteristics of the two ensembles of simulations are summarized in Table 1. Set B consists in 96 simulations with $1024^{3}$ particles in a cosmological volume of ($656.25$ $h^{-1}$Mpc)3. The total effective volume is 27 ($h^{-1}$Gpc)3 and the mass resolution is $1.88\times 10^{10}$ $h^{-1}$M*⊙. We call this set HR (“High Resolution”). Set C consists in 512 simulations with $512^{3}$ particles in a cosmological volume of ($1312.5$ $h^{-1}$Mpc)3. The total effective volume of 1158 ($h^{-1}$Gpc)3 is larger but the mass resolution of $1.20\times 10^{12}$ $h^{-1}$M⊙* is more limited. We call this set LR (“Low Resolution”). The two sets are complementary because they are affected at different level by numerical effects such as sample variance, finite volume and mass resolution effects (see Refs Blot et al. (2015); Rasera et al. (2014) for a study of these effects).

The matter bispectrum is estimated using the BISP_MES code Scoccimarro et al. (1998) kindly provided by S.Colombi. The code has been updated to take into account redshift space distortions (RSD) and projections onto multipoles. We provide here a short summary of the numerical methods, for more details see Scoccimarro et al. (1998). Particles position and velocity from a given snapshot are provided as an input of the code. Particles position are then displaced along the $\hat{\bm{z}}$-direction of the box using Eq. (1). When a particle falls outside of the simulation box, periodic boundary conditions are assumed to ensure conservation of the total number of particles. The density in Fourier-space is then computed using Cloud-In-Cell (CIC) mass assignment followed by a Fast Fourier Transform (FFT). The density field is further deconvolved with the CIC window function. We call the resulting density field from a given snapshot of the j-th simulation of a set, $\delta_{\rm sim,j}$. For triangles with sides $k_{1}$, $k_{2}$ and $k_{3}$, the multipole moments of the bispectrum are computed by a projection onto Legendre polynomials followed by an averaging over an homogeneous sample of modes within a bin of size $\Delta k_{1}$, $\Delta k_{2}$ and $\Delta k_{3}$ centered on $k_{1}$, $k_{2}$ and $k_{3}$. The averaging procedure is performed by randomly picking $N_{\rm mode}$ possible orientations and sizes of triangles within this interval (Monte-Carlo method). The estimated multipole projection of the bispectrum for a given snapshot of the j-th simulation is then given by

[TABLE]

where the subscript $i$ indicates the rank of the orientation/size of the triangle.

The number of orientations/sizes $N_{\rm mode}$ used for the averaging procedure is set to $10^{7}$. With such a large value, the results are insensitive to the exact choice of $N_{\rm mode}$. The size of the bin is chosen to be equal to the fundamental frequency of the box. We also use an FFT grid with $512^{3}$ elements. The corresponding Nyquist frequency is 2.45 $h$ Mpc*-1* (1.23 $h$ Mpc*-1*) for the HR (LR) simulation. As a consequence, we do not apply any shot noise corrections nor any aliasing corrections since such effects are negligible at the scale of interest in this paper ($k<0.3$ $h$ Mpc*-1*). However, at the very large scales of interest for comparison to perturbation theory ($k=0.01-0.1$ $h$ Mpc*-1*), HR simulations might suffer from non-negligible finite-mode sampling Takahashi et al. (2008). In the linear regime, the bispectrum of the density field should be zero for perfect ensemble average. Because the total number of independent modes in the simulation is finite, the resulting bispectrum is non-zero. To mitigate this effect at large scales, we compute the initial bispectrum of each snapshot $B^{(s)}_{\ell,{\rm ini,j}}$. The linearly evolved bispectrum is then subtracted from the snapshot bispectrum to obtain a corrected bispectrum $B^{(s)}_{\ell,{\rm corr,j}}$. The linearly evolved bispectrum is computed using Eq. (50) but instead of using the snapshot density $\delta_{\rm sim,j}^{(s)}(\bm{k}^{i})$, we use the linearly evolved density field

[TABLE]

where $f(z)$ is the linear growth rate and $D(z)$ is the linear growth factor. This correction plays a role at the percent level at small k ($k<0.1$ $h$ Mpc*-1*) for the HR run.

For each snapshot we consider isosceles triangles ($k_{1}=k_{2}$) and scalene triangles ($k_{1}=2k_{2}$). For each type of triangle, we explore the scale dependence by fixing $\theta_{12}$ and varying $k_{1}$ and the shape dependence by fixing $k_{1}$ and varying $\theta_{12}$. Once the bispectrum is computed for all simulations of a given set, we perform an ensemble average of $B^{(s)}_{\ell,{\rm corr,j}}(k_{1},k_{2},k_{3})$ to get the mean bispectrum of the set

[TABLE]

where $N_{\rm sim}$ is the number of simulations of the set. This bispectrum estimate is the one used in the rest of the paper. We also compute the standard deviation of the bispectrum and we estimate statistical error bars assuming

[TABLE]

The same analysis is repeated in comoving space (i.e. no redshift space distortions) by setting the velocity field and linear growth rate to zero.

IV.2 Results in real space

Let us first look at the results in real space and check the applicable range of PT as well as the quality of $N$-body data.

Fig. 4 shows the power spectra at $z=1$ (top), $0.5$ (middle) and [math] (bottom). Left panel plots the power spectra multiplied by $k^{3/2}$, while the right panel summarizes their fractional difference, for which we take the predictions by the Emulator code Heitmann et al. (2010, 2009); Lawrence et al. (2010); Heitmann et al. (2014) as the base value, and evaluate $[P(k)-P_{\rm Emu}(k)]/P_{\rm Emu}(k)$ with $P_{\rm Emu}$ being the Emulator power spectrum. Note here that the error bars in the $N$-body results indicate the standard error of the averaged power spectrum over the number of realizations. The claimed error bars of the Emulator are $1\%$.

The PT predictions at one-loop order (red) reasonably agree with those of the emulator code, and the agreement is at the 3 percent level for the scales of $k\lesssim 0.18$, $0.12$ and $0.1\,h$ Mpc*-1* at $z=1$, $0.5$ and [math], respectively. This is consistent with what has been found in the literature. The high-resolution $N$-body data (HR) also shows a reasonable agreement with the emulator (1 percent level over all the studied range of wavenumber) and one-loop PT predictions, but the low-resolution data (LR) systematically deviates from others at small scales. The deviation gradually increases from about 0.5 percent near $k=0.05\,h$ Mpc*-1*) to 5 percent near $k=0.2\,h$ Mpc*-1*. This is a well known mass-resolution effect which tends to decrease the power at small scale Heitmann et al. (2010); Rasera et al. (2014). Indeed, in the case of high-resolution setup, we found that the simulations starting at a lower initial redshift gives a better agreement with PT predictions.

On the other hand, turning to look at the real-space bispectrum, we do not clearly see the systematic difference between HR and LR simulations. Fig. 4 shows the results for the equilateral triangular shapes plotted as function of $k$ (left two panels), and those for the isosceles triangles plotted as function of $\theta_{12}\equiv\cos^{-1}(\hat{\bm{k}}_{1}\cdot\hat{\bm{k}}_{2})$ (right two panels). Note that in panels showing the amplitude of bispectrum, the results are all multiplied by $k^{3}$. Also, in plotting the fractional difference (second left and rightmost panels), we take the one-loop PT predictions as the base model, and evaluate the ratio, $(B-B_{\rm 1\mbox{-}loop})/B_{\rm 1\mbox{-}loop}$.

Compared to the power spectrum case, the statistical errors in $N$-body simulations is larger, and the difference between the two data set can be seen only at large scales $k\lesssim 0.05\,h$ Mpc*-1*. Rather, there seems to be a systematic difference between simulations and PT prediction, and because of this, the agreement between simulations and PT prediction look somewhat worse, and is at $5\%$ level even at the scales where the reasonable agreement at a few percent level can be seen in the power spectrum. Perhaps, these results might be partly ascribed to the setup of initial conditions (early starting redshift) or to the overestimate of non-linear effects on the bispectrum by standard PT, but without any reference, it is difficult to further clarify the origin of systematics. We will leave it to future investigation. Nevertheless, one important point is that the one-loop PT predictions can capture the major trend in the $N$-body simulations at small scales; scale-dependent enhancement of the bispectrum amplitude at $k\gtrsim 0.05\,h$ Mpc*-1*. Because of this, the agreement still remains at the $5-10$% level even at small scales at all redshifts.

Keeping the systematics and a level of agreement in the real space in mind, we will move to the redshift space, and continue the comparison in next subsection.

IV.3 Results in redshift space

To see the impact of FoG damping on the bispectrum, we first compare the standard PT predictions with results of $N$-body simulations.

Figs. 6 and 6 respectively show the scale- and shape-dependencies of the bispectrum for equilateral and isosceles configurations. Based on the definition of multipole expansion in Eq. (49) and Appendix C, we compute and measure the monopole and quadrupole moments of the bispectrum, which are respectively plotted in left and right panels. Overall, the $N$-body results both from the LR and HR data agree well each other. A closer look at $k\gtrsim 0.1\,h$ Mpc*-1* reveals that the amplitude of the LR data is systematically larger than that of the HR data, but within the errorbars there is no significant discrepancy in both monopole and quadrupole moments. Rather, a discrepancy between the $N$-body results and PT predictions is manifest, and compared to the results in real space, the range of agreement is fairly narrower for standard PT one-loop, and is restricted to a low-$k$ region. This is rather manifest at higher redshifts, and the prediction generally overestimates the simulations. On the other hand, shifting the overall amplitude, the shape dependence predicted by standard PT one-loop seems to reasonably match the measured results of $N$-body simulations (see Fig. 6). The results clearly manifest that a naive standard PT fails to describe the damping behaviors seen in the $N$-body simulations, and an appropriate prescription for the damping effect is important for prediction even at large scales.

Let us then see how the PT model presented in Sec. III works well. In doing so, an appropriate functional form of $D_{\rm FoG}$ needs to be first specified. The function $D_{\rm FoG}$ is generally expressed as function of $k_{1}\mu_{1}$, $k_{2}\mu_{2}$, and $k_{3}\mu_{3}$, and it contains the non-perturbative damping behavior arising from the exponential factor, $\exp\{\langle e^{j_{4}A_{4}+j_{5}A_{5}}\rangle_{c}\}$ [see Eq. (30) and proposition (i) in Sec. III]. Here, we adopt the following Gaussian form as a simple and relevant ansatz, similar to what has been frequently used in the power spectrum cases Taruya et al. (2010):

[TABLE]

The function $f$ is the linear growth rate, and $\sigma_{v}$ is the constant parameter corresponding to the one-dimensional velocity dispersion, which controls the strength of the FoG damping. Note that similar functional form is obtained by expanding the exponential pre-factor and truncating it at the second-order in cumulants, just ignoring the spatial correlation. A non-trivial point may be whether $D_{\rm FoG}$ is still expressed as univariate function of $(k_{1}\mu_{1})^{2}+(k_{2}\mu_{2})^{2}+(k_{3}\mu_{3})^{2}$ or not beyond the scales relevant for tree-level predictions. We will discuss and check it in Sec. IV.4.

Adopting Eq. (54), Figs. 8-10 compare the prediction of the PT model with $N$-body simulations. While Figs. 8 and 8 show the scale-dependence of the bispectrum amplitudes for equilateral ($k_{1}=k_{2}=k_{3}$) and scalene triangular configuration with $k_{1}=2k_{2}=(2/\sqrt{3})k_{3}$, Figs. 10 and 10 respectively plot the shape-dependence of the bispectrum as function of $\theta_{12}=\cos^{-1}(\hat{\bm{k}}_{1}\cdot\hat{\bm{k}}_{2})$ for the triangles of $k_{1}=k_{2}$ and $k_{1}=2k_{2}$. In each figure, the monopole and quadrupole moments of the bispectrum are computed/measured according to the definition in Appendix C, and the results are presented in left and right panels, respectively. Here, the measured results of the bispectra are shown only for LR data, since no notable difference has been found in both LR and HR data.

Overall, the one-loop PT model depicted as magenta solid lines better agree with simulations over a wider range of $k$ as well as for a wide range of shapes. Note that the free parameter $\sigma_{\rm v}$ in Eq. (54) is determined at each redshift by fitting the predicted monopole and quadrupole moments with measured results of $N$-body simulations at the range $[k_{\rm min},\,k_{\rm max}]$. While $k_{\rm min}$ is set to $0.05\,h$ Mpc*-1*, we adopt $k_{1\%}$ defined by Ref. Nishimichi et al. (2009) as the maximum wavenumber $k_{\rm max}$, indicated by the vertical arrow in each panel of Figs. 8 and 8 (blue for tree-level PT and red for one-loop PT). The $k_{1\%}$ indicates the maximum wavenumber below which the predicted power spectrum is shown to well reproduce the N-body result within $1\%$ accuracy in the real space, and from Fig. 4, we see that the $k_{1\%}$ also gives a good indicator for the applicable range of one-loop bispectrum in real space. Because the deviations between PT and simulations do not behave as a monotonous function it is however hard to extract an exact scale. This is why we rely on more accurate power spectrum measurements to define $k_{1\%}$. The fitted result of the parameter $\sigma_{\rm v}$ is given in left panel of each figure, which are close to the linear theory prediction of one-dimensional velocity dispersion (see also Table 2, and Figs. 12 and 12).

A notable point may be that the one-loop PT model reproduces the $N$-body results even beyond the fitting range of $\sigma_{\rm v}$. For comparison, in Figs. 8-10, we have also plotted the tree-level standard PT predictions multiplied by the damping function of Eq. (54)222In terms of the descriptions given in Sec. III, the tree-level standard PT multiplied by $D_{\rm FoG}$ corresponds to the leading-order PT calculations of Eq. (45) with the functions $C_{n}$ summing up to $n=2$. , depicted as blue short-dashed lines, but the agreement with $N$-body simulations is restricted to the fitting range indicated by the blue vertical arrow. In this respect, the one-loop corrections play an important role, together with damping function, to better describe the redshift-space bispectrum at weakly non-linear regime. A closer look at the equilateral case in Fig. 8 reveals that the agreement of the one-loop model is bit degraded compared to the scalene case (Fig. 8), especially for the quadrupole moment. Similar trend is also found for the shape dependence in Figs. 10 and 10. This is partly because for a fixed wavenumber $k$, one of the side length for the triangle becomes smaller or larger than $k$, and the results can be less or more sensitive to the nonlinearity of the gravitational clustering and RSD. Indeed, for the cases shown in Figs. 8 and 10, the tree-level PT predictions get closer to the one-loop PT results at large scales.

IV.4 Testing the ansatz of damping function

Adopting the Gaussian form of the non-perturbative damping function, we have seen that the PT model successfully describes the measured results of the bispectrum. However, the validity and consistency of the treatment have to be checked and/or verified, at least in the following two aspects. One is the consistency of the fitted values of $\sigma_{\rm v}$. Since the present PT model is constructed in similar manner to the model of power spectrum in Ref. Taruya et al. (2010), the fitted values of the parameter $\sigma_{\rm v}$ derived from bispectrum have to be consistent with those from the power spectrum. The other aspect is the functional form of the damping function. We have assumed in Eq. (54) that the damping function is expressed as the univariate function of $(k_{1}\mu_{1})^{2}+(k_{2}\mu_{2})^{2}+(k_{3}\mu_{3})^{2}$, but this can be verified only at the leading-order, and have to be checked at the scales relevant for the one-loop PT.

Let us first discuss the consistency of the fitted value, $\sigma_{\rm v}$. Table 2 summarizes the results derived both from the bispectrum and power spectrum. The estimated results are based on the $N$-body simulations of the LR (HR) data, and we used the one-loop model for power spectrum, while for the bispectrum, the fitted results from two different configurations are presented in both tree-level and one-loop PT cases. The best-fitted values both from the power spectrum and bispectrum one-loop models reasonably agree with each other at $z=1$. Although a deviation is manifest at lower redshifts, this is small compared to the cases with tree-level predictions of bispectrum.

To see the robustness of the fitted values, we extend the analysis in Table 2, and using the LR data, we examine the fitting in various range of $k$. The results are shown in Figs. 12 and 12. In each figure, top and bottom panels respectively plot the fitted values $\sigma_{\rm v}$ and the reduced $\chi^{2}$ as function of $k_{\rm max}$. Note that in estimating $\chi^{2}_{\rm red}$ and the errors of the fitted $\sigma_{\rm v}$, we took account of the systematics in the $N$-body simulations. That is, at each data point, we added the systematic errors $\Delta P_{\ell,{\rm sys}}^{\rm(s)}$ and $\Delta B_{\ell,{\rm sys}}^{\rm(s)}$ to the statistical errors of the power spectrum and bispectrum multipoles, $\Delta P_{\ell,{\rm stat}}^{\rm(s)}$ and $\Delta B_{\ell,{\rm stat}}^{\rm(s)}$, as shown in Figs. 8-10 for the bispectrum case (solid error bars: $\Delta B^{\rm(s)}_{\ell,{\rm stat}}$, dotted error bars: $\Delta B^{\rm(s)}_{\ell,{\rm stat}}+\Delta B^{\rm(s)}_{\ell,{\rm sys}}$). The size of the systematics is estimated from the real-space results in Figs. 4 and 4. We assume that the relative systematics in redshift space for the power spectrum and bispectrum are proportional to the relative systematics in real space for the matter power spectrum, $\Delta P_{\rm sys}/P_{\rm sim}=|P_{\rm sim}-P_{\rm ref}|/P_{\rm sim}$, where $P_{\rm sim}$ represents the measured real-space power spectrum in $N$-body simulations, and $P_{\rm ref}$ is our reference for the real-space power spectrum. For the reference, we choose the emulator power spectrum (although the estimated systematics is of order 1 percent, we neglect this contribution). The systematics for the power spectrum and bispectrum in redshift space are estimated as,

[TABLE]

where $\alpha$ and $\beta$ are two fudge factors that we have fixed to $\alpha\sim 1$ and $\beta\sim 1$. We have checked that at large-k (where statistical error bars are small compared to systematics), Eq. (56) reproduces the order of magnitude of the relative difference between the bispectrum of the HR set and the LR set (which is affected by mass resolution effect).

The rough systematic errors adopted here may result in a rather crude estimate of the goodness of fit, and thus the derived $\chi_{\rm red}^{2}$ can only be used for a comparison between tree-level and one-loop results in Fig. 12 and 12. Nevertheless, we see that the fitting results in one-loop PT cases are basically stable against the variation of $k_{\rm max}$, and the estimated values of $\chi_{\rm red}^{2}$ are smaller than those in the tree-level cases. Further, we checked that the (best-)fitted values of $\sigma_{\rm v}$ are robust against the systematic errors, and the resultant values in one-loop PT reasonably agree well with each other, especially at $z=1$. Although the deviation becomes manifest at lower redshift, this would be probably due to the break down of the one-loop predictions. In fact, the $\chi_{\rm red}^{2}$ systematically increases with $k_{\rm max}$, indicating that the fitting starts to fail. Thus, at least at the redshift $z=1$, the one-loop PT models work fine, and the FoG damping is described with a single parameter.

Next consider the validity of the ansatz for $D_{\rm FoG}$ at Eq. (54). To clarify whether the non-perturbative part is described by the univariate function or not, we directly measure the bispectrum, not applying the multipole expansion. In this case, the bispectrum in redshift space is described by the five variables. For each shape and orientation of the bispectrum, we compute the corresponding PT prediction based on Eq. (45) or (48), but setting $D_{\rm FoG}$ to $1$. Taking the ratio gives

[TABLE]

At the scales where the one-loop PT is applicable, this ratio directly quantifies the functional form of $D_{\rm FoG}$, and thus we can check whether it is expressed as univariate function of $(k_{1}\mu_{1})^{2}+(k_{2}\mu_{2})^{2}+(k_{3}\mu_{3})^{2}$ or not.

Fig. 13 shows the measured results of the ratio, Eq. (57), for various orientations at $z=1$, plotted against $(k_{1}\mu_{1})^{2}+(k_{2}\mu_{2})^{2}+(k_{3}\mu_{3})^{2}$. To be precise, what is shown here is the dispersion of the measured ratio depicted as shaded color region, and the four different color strengths imply the length of $k_{1}$ used to estimate the ratio: $k_{1}\leq 0.05$, $0.15$, $0.5$, and $1\,h$ Mpc*-1* from dark to light. The results are compared with the univariate damping function of the Gaussian (dashed) and Lorentzian (dotted) form. Clearly, the scatter of the ratio for the one-loop PT model, given in top panel, is small, and its mean values fairly trace the univariate damping function. This is in marked contrast to the results for tree-level PT shown in bottom panel, where we see a large scatter. Further, the results seem robust irrespective of the shape of the bispectrum triangle, as seen in both left and right panels, where we respectively show the results for the equilateral case ($k_{1}=k_{2}=k_{3}$) and the scalene triangle with $k_{1}=2k_{2}=(2/\sqrt{3})k_{3}$. A closer look at results suggests that Lorentzian form describes the measured ratio reasonably well at the high-$k$ tail, although it is mostly the boundary where we can apply one-loop PT prediction. Hence, we conclude that the univariate ansatz for $D_{\rm FoG}$ is validated at least in the applicable range of one-loop predictions.

V Discussion and conclusion

In this paper, we have studied the matter bispectrum in redshift space, and presented a perturbation theory (PT) model that can keep the non-perturbative damping effect of the redshift-space distortions (RSD) under control. Starting with the exact formula for redshift-space bispectrum, we rewrite the expression in terms of the cumulants to identify the non-perturbative term. Separating the non-perturbative term responsible for the so-called Fingers-of-God (FoG) damping effect, we derive the perturbative expressions for bispectrum valid at one-loop order. The resultant model has been constructed similarly to the power spectrum model in Ref. Taruya et al. (2010), and it incorporates the non-perturbative damping term on top of the terms that can be computed with standard PT.

Adopting the Gaussian form of damping function, we have performed a detailed comparison between the predictions of PT model with the measured results of the bispectrum from a suite of cosmological $N$-body simulations. Incorporating a single free parameter into the damping function, the one-loop PT model reproduces the simulation results fairly well at weakly nonlinear scales at $z=0-1$. The fitted results of the parameter $\sigma_{\rm v}$ are found to agree well with those obtained from the power spectrum, and the agreement generally holds irrespective of the shape of the triangles. On the other hand, even if we incorporate the damping function into the model, the tree-level PT predictions start to deviate from $N$-body results at rather low-$k$. Also, the fitted value $\sigma_{\rm v}$ does not match the one obtained from the power spectrum, and varies with triangular shapes.

We have further examined the validity of the ansatz imposed in the functional form of the damping function [Eq. (54)]. Combining the simulation data with standard PT results, we confirmed that the univariate ansatz for damping function $D_{\rm FoG}$ indeed hold for one-loop PT model, and its functional form is shown to be very close to the Gaussian, although the Lorentzian form looks slightly better. Note, however, that in the case of the tree-level PT model, univariate ansatz does not give a good description, and this can be a part of the reason why the failure of tree-level PT prediction appears at rather low-$k$. Hence, even at large scales, a careful modeling with one-loop correction is essential, and together with the model of power spectrum proposed by Ref. Taruya et al. (2010), the present one-loop model of bispectrum gives a coherent description for RSD.

Finally, we note that toward the practical application to real applications, there still remain several issues to be addressed. One is the improvement of the PT prediction. Including the higher-order (two-loop) corrections or applying the resummation technique, the applicable range of PT is expected to become wider, and a more tighter test of gravity will be made. Effective field theory approach may also help to improve the prediction (e.g., Refs. Baumann et al. (2012); Carrasco et al. (2012); Hertzberg (2014); Baldauf et al. (2015)). Another issue is the estimation of statistical error covariance of bispectrum, which is crucial and necessary for unbiased and robust cosmological data analysis. The $N$-body measurement of the covariance is, however, known to be computationally extensive (e.g., Refs. Kayo et al. (2013); Sato and Nishimichi (2013) for weak lensing case), and a clever approach, involving the analytic treatment, may have to be developed, for instance. One final big issue is the galaxy biasing. Throughout the paper, we have focused on the matter bispectrum, but the real observable is the biased object. Incorporating the prescription of galaxy biasing is thus very crucial. Toward the practically useful model, the study with halo/subhalo catalogs may help a lot, and we will tackle this issue in the near future.

Acknowledgements.

This work is in part supported by MEXT/JSPS KAKENHI Grant Number JP15H05899 and JP16H03977. It is also supported by JSPS Grant L16519. This work was granted access to HPC resources of IDRIS through allocations made by GENCI (Grand Equipement National de Calcul Intensif) under the allocations 2015-042287 and 2016-042287. We acknowledge support from the DIM ACAV of the Region Ile-de-France. We also thanks Stéphane Colombi for providing his bispectrum code.

Appendix A Redshift-space perturbation theory kernels

In this Appendix, we summarize the explicit expressions for redshift-space kernels $Z_{n}$ defined in Eq. (11). Substituting the higher-order PT solutions given by Eqs. (9) and (10) into the Taylor-expanded form of the redshift-space density field in Eq. (8), the kernel $Z_{n}$ can be read off from the reorganized perturbative expansion in powers of $\delta_{\rm L}$. The kernels up to fourth order become (e.g., Scoccimarro et al. (1999)):

[TABLE]

Here, the fourth-order kernels $Z_{\rm 4a}$ and $Z_{\rm 4b}$ are respectively given by

[TABLE]

In the above, the vector $\bm{k}$ in the $n$-th order kernel implies $\bm{k}=\bm{k}_{1}+\cdots+\bm{k}_{n}$. The quantities $\mu$, $\mu_{i}$, $\mu_{ij}$, and $\mu_{ijk}$ are defined by:

[TABLE]

In applying these results to the statistical calculations described in Sec. II.2, the kernels $Z_{n}$ have to be symmetrized under the exchange of each argument. One important remark is that even with fully symmetrized kernels $F_{n}$ and $G_{n}$, the resultant redshift-space kernels $Z_{n}$ at $n\geq 3$ only preserve partial symmetry. For instance, while the kernel $Z_{3}$ given at Eq. (60) is symmetric under $\bm{k}_{1}\leftrightarrow\bm{k}_{3}$, the expression for $Z_{4a}$ preserve the symmetry of $\bm{k}_{2}\leftrightarrow\bm{k}_{3}\leftrightarrow\bm{k}_{4}$. These kernels become fully symmetrized if we take the cyclic permutations.

Appendix B Perturbative calculations for $D_{1}$ and $D_{2}$ terms

In this Appendix, we present the perturbative expressions for the $D_{1}$ and $D_{2}$ terms defined in Eq. (47), relevant at one-loop order calculations.

First note that Eq. (47) is expressed in terms of the real-space quantities, including the auto- and cross-power spectra and bispectra of density and velocity fields. For convenience, we introduce the two-component multiplet:

[TABLE]

with $\theta$ being the dimensionless velocity divergence defined in real space by $\theta(\bm{x})=-\nabla\cdot\bm{v}/(f\,aH)$. Then the auto- and cross-power spectra and bispectra are given by

[TABLE]

where subscripts $a,b,c$ run from $1$ to $2$.

Since the integrand of the $D_{1}$ and $D_{2}$ terms respectively involves $B_{abc}P_{de}$ and $P_{ab}P_{cd}P_{ef}$, the leading-order non-vanishing contributions become $\mathcal{O}(P_{\rm L}^{3})$. This is the same order as in the one-loop redshift-space bispectrum. Hence, the tree-level calculations of $P_{ab}$ and $B_{abc}$ are sufficient for a consistent one-loop treatment of redshift-space bispectrum. That is,

[TABLE]

where the kernel ${\mathcal{F}}_{a}$ is the standard PT kernels at second-order, given by ${\mathcal{F}}_{a}=(F_{2},\,G_{2})$. Below, we will separately present the explicit expressions for $D_{1}$ and $D_{2}$ terms. In doing so, we use the following expressions for the cumulants:

[TABLE]

with the quantities $\mu_{i}$ and $\mu_{p}$ respectively defined by $\mu_{i}=(\bm{k}_{i}\cdot\hat{\bm{z}})/k_{i}$ and $\mu_{p}=(\bm{p}\cdot\hat{\bm{z}})/p$.

B.1 $D_{1}$ term

According to Eq. (47), the $D_{1}$ term is explicitly written as

[TABLE]

Substituting Eq. (70) and (71) into the above, we obtain

[TABLE]

with the quantity $\sigma_{\rm v,lin}$ being the linear theory estimate of the one-dimensional velocity dispersion, given by

[TABLE]

Here, the function $C_{1}$ is the same one as defined in Eq. (33 ) with (34), and is explicitly given by

[TABLE]

B.2 $D_{2}$ term

From Eq. (47), the $D_{2}$ term is expressed in terms of the cumulants below:

[TABLE]

Substituting Eqs. (71)-(73) into the above expression, after lengthy calculation, we obtain

[TABLE]

Here, the function $C_{2}$ is the same one as defined in Eq. (33) with (35). The explicit expressions for the functions $C_{2}$, $J_{2}$ and $K_{2}$ are given below:

[TABLE]

with the functions $X$ and $Y$ defined by

[TABLE]

Appendix C Multipole expansion of redshift-space bispectrum

In Sec. III and IV, we have applied the multipole expansion to the redshift-space bispectrum, and evaluated its monopole and quadrupole moments. In this Appendix, we present the definition of our multipole expansion which differs from the one frequently used in the literature (e.g., Scoccimarro et al. (1999); Scoccimarro (2015); Yamamoto et al. (2016); Gagrani and Samushia (2017)). The newly defined bispectrum multipoles have several nice properties, which we will discuss below.

The redshift-space bispectrum is characterized as a function of five variables. Three of them characterize the shape of triangle, i.e., the length of two wave vectors $\bm{k}_{1}$ and $\bm{k}_{2}$, and the angle between them, $\theta_{12}\equiv\cos^{-1}(\hat{\bm{k}}_{1}\cdot\hat{\bm{k}}_{2})$. The two remaining variables describe the orientation of the triangle with respect to the line-of-sight. We denote them by $\omega$ and $\phi$. In Ref. Scoccimarro et al. (1999), the dependance of the orientation of the triangle are conveniently described by decomposing into spherical harmonics:

[TABLE]

Similar to Ref. Scoccimarro et al. (1999), we focus on the $m=0$ multipoles, which correspond to averaging over $\phi$. Then, the above equation leads to

[TABLE]

were $P_{\ell}(\mu)$ is the Legendre polynomials. The variable $\mu$ represents the directional cosine of the orientation, $\mu=\cos{\omega}$. This is rewritten with

[TABLE]

In the above, the bispectrum multipole, $B_{\rm s}^{(\ell)}$, is the quantity of our interest, but at this moment, Eq. (87) is ambiguous because we do not yet specify what is $\omega$ and $\phi$. To describe the orientation of the triangle, a simple way is to choose $\bm{k}_{1}$ specifically, and define the orientation angle as $\mu=\cos{\omega}=\hat{\bm{k}}_{1}\cdot\hat{\bm{z}}$ Scoccimarro et al. (1999). Then, we set $\phi$ to the azimuthal angle around $\bm{k}_{1}$. In this case, however, the resultant bispectrum multipoles $B_{\rm s}^{(\ell)}(k_{1},k_{2},\theta_{12})$ are not fully symmetric under the permutation of order, $k_{1}$, $k_{2}$ and $k_{3}$.

Here, we give alternative definitions of $\omega$ and $\phi$ to preserve the symmetry of bispectrum multipoles. This is shown in Fig. 1. Given two vectors $\bm{k}_{1}$ and $\bm{k}_{2}$, we define the orientation of the triangle by the angle between line-of-sight direction and the vector normal to the triangle, i.e., $\hat{\bm{k}}_{1}\times\hat{\bm{k}}_{2}$. That is,

[TABLE]

The remaining angle $\phi$ may be defined as the azimuthal angle around the vector $\hat{\bm{k}}_{1}\times\hat{\bm{k}}_{2}$. Since this is perpendicular to the plane of the triangle, we have

[TABLE]

One may suspect that the above definition still breaks the symmetry. Indeed, the vectors $\bm{k}_{1}$ and $\bm{k}_{2}$ can be arbitrarily chosen among three vectors, and we can even exchange $\hat{\bm{k}}_{1}\,\longleftrightarrow\,\hat{\bm{k}}_{2}$. Then the orientation angle is changed to $\omega\to\pi-\omega$. Nevertheless, the bispectrum multipole $B_{\rm s}^{(\ell)}$ is invariant. This is because only the $\ell=$ even modes of the bispectrum become non-vanishing, and we have $P_{\ell}(\cos\omega)=P_{\ell}(\cos(\pi-\omega))$ for $\ell=$ even.

Once accepting the new definition, the measurement of the bispectrum multipoles is straightforward:

• •

In harmonic space, we first pick up the three density fields, $\delta(\bm{k}_{1})$, $\delta(\bm{k}_{2})$, and $\delta(\bm{k}_{3})$, with the vectors satisfying $\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3}=\bm{0}$. Here, labels of $\bm{k}_{1}$, $\bm{k}_{2}$ and $\bm{k}_{3}$ are arbitrary.

• •

Choosing the two vectors among the three, we calculate the directional cosine $\mu$ according to Eq. (88).

• •

Multiplying the bispectrum estimator by the Legendre polynomial $P_{\ell}(\mu)$.

• •

Repeating the above three steps for the same triangle configuration but with different $\mu$, we average the weighted bispectrum estimator over $\mu$. Further multiplying the averaged bispectrum by the factor $(2\ell+1)/2$, we finally obtain the bispectrum multipole, $B_{\rm s}^{(\ell)}$, which is characterized by the length of two vectors, and the angle between them [see Eq. (87)].