<figure><img src="image_1.png" /></figure>

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

Ichihiko Hashimoto, Yann Rasera, Atsushi Taruya

► To cite this version:

Ichihiko Hashimoto, Yann Rasera, Atsushi Taruya. Precision cosmology with redshift-space bispectrum: a perturbation theory based model at one-loop order. Physical Review D, 2017, 96 (4), pp.043526. 10.1103/PhysRevD.96.043526. hal-01645374

HAL Id: hal-01645374

https://hal.science/hal-01645374v1

Submitted on 5 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Precision cosmology with redshift-space bispectrum: A perturbation theory based model at one-loop order

Ichihiko Hashimoto,¹ Yann Rasera,²,³ and Atsushi Taruya³,⁴

¹Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan

²LUTH, Observatoire de Paris, PSL Research University, CNRS, Université Paris Diderot, Sorbonne Paris Cité, 5 place Jules Janssen, 92195 Meudon, France

³Center for Gravitational Physics, Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan

⁴Kavli Institute for the Physics and Mathematics of the Universe (WPI), Todai Institute for Advanced Study, University of Tokyo, Kashiwa, Chiba 277-8568, Japan
(Received 17 May 2017; published 31 August 2017)

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. In this paper, extending the work by Taruya, Nishimichi, and Saito [Phys. Rev. D **82**, 063522 (2010)], we present a perturbation theory (PT) based model of redshift-space matter bispectrum that can keep the nonperturbative damping effect under control. Characterizing this nonperturbative damping by a univariate function with a 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 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, and Saito.

DOI: 10.1103/PhysRevD.96.043526

## 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. [1–5]), 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 for by the linear theory, and the size of the anisotropies is related to the growth of structure (e.g., Refs. [2,3]). To be precise, the measurement of RSD offers an opportunity to determine the parameter $f\sigma_8$ (e.g., Refs. [6–8]). 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 $8h^{-1}$ Mpc. Combining the RSD with the measurement of baryon acoustic oscillations (BAO) imprinted on galaxy clustering (e.g., [9–11]), we can simultaneously estimate both the geometric distances and growth of structure (e.g., [12–15]), from which we are able to scrutinize the theory of gravity that accounts for cosmic expansion and structure formation (e.g., Refs. [16–20]).

The measurements of RSD taken from various observations now cover a wide range of redshifts out to $z \sim 1.4$ (e.g., Refs. [21–25]; see also Sec. 5.5 of Ref. [26]), 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 ΛCDM model [26], 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. [27–31]). In particular, the impact of combining bispectrum measurement is demonstrated for the RSD measurement by Refs. [32,33], 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 2.

However, most of these analyses have been demonstrated in a very simplified setup, assuming that the observed bispectrum is described by the perturbation theory (PT) at leading order [34], 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 the 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 it is significantly affected by the small-scale physics (e.g., Ref. [4]). Hence, a sophisticated treatment is needed for a reliable theoretical template even in the weakly nonlinear regime [12,35–41].

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 the redshift-space bispectrum, which consistently incorporates the effect coming from the small-scale virial motion. While modeling the redshift-space bispectrum has already been investigated based on the fitting formula [42,43] and the halo model (e.g., Refs. [44,45]), analytical treatment would certainly be powerful in characterizing the anisotropies of the redshift-space bispectrum (see, e.g., [46] 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 we derive the exact formulas for power spectrum and bispectrum in redshift space. Based on the standard PT, the expressions for the redshift-space bispectrum at one-loop order are also presented. Then, in Sec. III, extending the work by Ref. [12], we give the PT-based model of the redshift-space bispectrum, which consistently includes both the next-to-leading order PT corrections and the nonperturbative 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 the $\Lambda$CDM model and adopt the cosmological parameters determined from the seven-year WMAP results [47]: $h = 0.72$, $\Omega_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 the matter density field in redshift space. We then consider the standard PT calculations, and based on the explicit expressions for the bispectrum up to the next-to-leading order called one-loop, we discuss its limitation and applicability.

### A. Exact formula of power spectrum and bispectrum

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

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

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

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

Since we are interested in the statistical quantities in Fourier space, we rewrite the above relation in terms of the Fourier counterpart,
$$ \delta^{(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_z / aH + k \cdot \mathbf{r})}, \quad (3) $$

where the variable $\mu$ is defined as the directional cosine of the angle between the line-of-sight $\hat{z}$ and wave vector $\mathbf{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

$$ (2\pi)^3 \delta_D(\mathbf{k} + \mathbf{k}') P^{(s)}(\mathbf{k}) \equiv \langle \delta^{(s)}(\mathbf{k}) \delta^{(s)}(\mathbf{k}') \rangle, \quad (4) $$

$$ \begin{aligned} & (2\pi)^3 \delta_D(\mathbf{k}_1 + \mathbf{k}_2 + \mathbf{k}_3) B^{(s)}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) \\ & \equiv \langle \delta^{(s)}(\mathbf{k}_1) \delta^{(s)}(\mathbf{k}_2) \delta^{(s)}(\mathbf{k}_3) \rangle, \end{aligned} \quad (5) $$

substituting Eq. (3) into the above leads to

$$ P^{(s)}(\mathbf{k}) = \int dr e^{i\mathbf{k}\cdot\mathbf{r}_{12}} \times \langle e^{-if\mathbf{k}_z u_{12z}} \{\delta(\mathbf{r}_1) + f\nabla_z u_z(\mathbf{r}_1)\} \times \{\delta(\mathbf{r}_2) + f\nabla_z u_z(\mathbf{r}_2)\} \rangle \quad (6) $$

for the redshift-space power spectrum, and

$$ \begin{aligned} B^{(s)}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = & \int dr \, \mathbf{r}_{13} \mathbf{r}_{23} e^{i(\mathbf{k}_1 \cdot \mathbf{r}_{13} + \mathbf{k}_2 \cdot \mathbf{r}_{23})} \\ & \times \langle e^{-if(\mathbf{k}_1 u_{13z} + \mathbf{k}_2 u_{23z})} \{\delta(\mathbf{r}_1) + f \nabla_z u_z(\mathbf{r}_1)\} \\ & \times \{\delta(\mathbf{r}_2) + f \nabla_z u_z(\mathbf{r}_2)\} \\ & \times \{\delta(\mathbf{r}_3) + f \nabla_z u_z(\mathbf{r}_3)\} \rangle \end{aligned} \quad (7) $$

for the redshift-space bispectrum. Here, we define the normalized peculiar velocity $\mathbf{u}(\mathbf{r}) = -\mathbf{v}(\mathbf{r})/\{faH(z)\}$ and denote the pairwise normalized velocity of the separation $\mathbf{r}_{ij} = \mathbf{r}_i - \mathbf{r}_j$ by $\mathbf{u}_{ij} \equiv \mathbf{u}_i - \mathbf{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 relations. 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., the length of three wave vectors $k_1, k_2$, and $k_3$, or the length of vectors $k_1$ and $k_2$ and their angle $\theta_{12} \equiv \cos^{-1}(\hat{\mathbf{k}}_1 \cdot \hat{\mathbf{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 (FoG) effects are rather nontrivial, and a careful treatment is required for an accurate modeling of the redshift-space bispectrum.

## B. Standard perturbation theory in redshift space

In this paper, we employ the perturbation theory technique to construct an analytic model of the 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 the next section, a PT model of the redshift-space bispectrum which properly incorporates the effect of FoG damping in a nonperturbative manner.

To derive the standard PT results of the 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

$$ \begin{aligned} \delta^{(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_D(\mathbf{k} - \mathbf{q}_{1\cdots n}) \\ & \times \{\delta(\mathbf{k}) + f\mu^2 \theta(\mathbf{k})\} \frac{(f\mu k)^n}{n!} \frac{\mu_1}{q_1} \theta(\mathbf{q}_1) \cdots \frac{\mu_n}{q_n} \theta(\mathbf{q}_n) \end{aligned} \quad (8) $$

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

$$ \delta(\mathbf{k}) = \sum_{n=1} \int \frac{d^3 \mathbf{p}_1 \cdots d^3 \mathbf{p}_n}{(2\pi)^{3n}} \delta_D(\mathbf{k} - \mathbf{p}_{1\cdots n}) F_n(\mathbf{p}_1, \cdots, \mathbf{p}_n) \\ \times \delta_L(\mathbf{p}_1) \cdots \delta_L(\mathbf{p}_n), \quad (9) $$

$$ \theta(\mathbf{k}) = \sum_{n=1} \int \frac{d^3 \mathbf{p}_1 \cdots d^3 \mathbf{p}_n}{(2\pi)^{3n}} \delta_D(\mathbf{k} - \mathbf{p}_{1\cdots n}) G_n(\mathbf{p}_1, \cdots, \mathbf{p}_n) \\ \times \delta_L(\mathbf{p}_1) \cdots \delta_L(\mathbf{p}_n), \quad (10) $$

where $\delta_L$ is the linear density field, and $F_n$ and $G_n$ are the standard PT kernels (see, for example, [48]). Plugging the above expressions into Eq. (8), reorganizing the perturbative expansion in powers of $\delta_L$ leads to

$$ \begin{aligned} \delta^{(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_L(\mathbf{q}_1) \cdots \delta_L(\mathbf{q}_n). \end{aligned} \quad (11) $$
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

$$P^{(s)}(\mathbf{k}) = P^{\text{lin}} + P^{\text{1-loop}} + \dots, \quad (12)$$

$$B^{(s)}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = B^{\text{tree}} + B_{222}^{\text{1-loop}} + B_{321-\text{I}}^{\text{1-loop}} + B_{321-\text{II}}^{\text{1-loop}} \\ + B_{411}^{\text{1-loop}} + \dots \quad (13)$$

Each term at the right-hand side is given by

$$P^{\text{lin}}(\mathbf{k}) = \{Z_1(\mathbf{k})\}^2 P_L(\mathbf{k}), \quad (14)$$

$$P^{\text{1-loop}}(\mathbf{k}) = \int \frac{d^3 p}{(2\pi)^3} \{Z_2(\mathbf{p}, \mathbf{k}-\mathbf{p})\}^2 P_L(p) P_L(|\mathbf{k}-\mathbf{p}|) + 2Z_1(\mathbf{k}) P_L(k) \int \frac{d^3 p}{(2\pi)^3} \{Z_3(\mathbf{p}, -\mathbf{p}, \mathbf{k})\} P_L(p) \quad (15)$$

for the power spectrum, and

$$B^{\text{tree}}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = 2Z_2(\mathbf{k}_1, \mathbf{k}_2)Z_1(\mathbf{k}_1)Z_1(\mathbf{k}_2)P_L(\mathbf{k}_1)P_L(\mathbf{k}_2) + 2 \text{perms}(\mathbf{k}_1 \leftrightarrow \mathbf{k}_2 \leftrightarrow \mathbf{k}_3), \quad (16)$$

$$B_{222}^{\text{1-loop}}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = \int \frac{d^3 p}{(2\pi)^3} Z_2(\mathbf{p}, \mathbf{k}_1 - \mathbf{p}) Z_2(-\mathbf{p}, \mathbf{k}_2 + \mathbf{p}) Z_2(-\mathbf{k}_1 + \mathbf{p}, -\mathbf{k}_2 - \mathbf{p}) P_L(p) P_L(|\mathbf{k}_1 - \mathbf{p}|) P_L(|\mathbf{k}_2 + \mathbf{p}|) \\ + 2 \text{perms}(\mathbf{k}_1 \leftrightarrow \mathbf{k}_2 \leftrightarrow \mathbf{k}_3), \quad (17)$$

$$B_{321-\text{I}}^{\text{1-loop}}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = Z_1(\mathbf{k}_1) P_L(\mathbf{k}_1) \int \frac{d^3 p}{(2\pi)^3} Z_2(\mathbf{p}, \mathbf{k}_2 - \mathbf{p}) Z_3(-\mathbf{k}_1, -\mathbf{p}, -\mathbf{k}_2 + \mathbf{p}) P_L(p) P_L(|\mathbf{k}_2 - \mathbf{p}|) \\ + 5 \text{perms}(\mathbf{k}_1 \leftrightarrow \mathbf{k}_2 \leftrightarrow \mathbf{k}_3), \quad (18)$$

$$B_{321-\text{II}}^{\text{1-loop}}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = Z_1(\mathbf{k}_1)Z_2(\mathbf{k}_1, \mathbf{k}_2)P_L(\mathbf{k}_1)P_L(\mathbf{k}_2) \int \frac{d^3 p}{(2\pi)^3} Z_3(\mathbf{k}_2, \mathbf{p}, -\mathbf{p})P_L(p) + 5 \text{perms}(\mathbf{k}_1 \leftrightarrow \mathbf{k}_2 \leftrightarrow \mathbf{k}_3), \quad (19)$$

$$B_{411}^{\text{1-loop}}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = Z_1(\mathbf{k}_1)Z_1(\mathbf{k}_2)P_L(\mathbf{k}_1)P_L(\mathbf{k}_2) \int \frac{d^3 p}{(2\pi)^3} Z_4(-\mathbf{k}_1, -\mathbf{k}_2, \mathbf{p}, -\mathbf{p})P_L(p) + 2 \text{perms}(\mathbf{k}_1 \leftrightarrow \mathbf{k}_2 \leftrightarrow \mathbf{k}_3) \quad (20)$$

for the bispectrum. In deriving the equations above, we assume the Gaussianity of linear density field $\delta_L$ and denote its power spectrum by $P_L$. These expressions can be reduced to the real-space power spectrum and bispectrum of the 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 the 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 \ge 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. [12], 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 the one-loop prediction thus becomes narrower than that of the real-space prediction. Since the expression of the bispectrum also uses the redshift PT kernels and it even includes the higher order, the situation must be similar to or rather worse than that in the power spectrum case. We will see in Sec. IV that without accounting for the finger-of-God damping, the standard PT prediction of the 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 incorporate the damping effect into the PT calculation in the next section.

### III. AN IMPROVED MODEL PRESCRIPTION FOR BISPECTRUM AT ONE-LOOP ORDER

In this section, we present the PT model of the redshift-space bispectrum which keeps the nonperturbative damping effect. Our strategy is to decompose the contributions into a nonperturbative 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. [12], and rewrite the exact expression in terms of cumulants. We then identify the nonperturbative part responsible for the FoG damping. Based on the simple proposition similarly made by Ref. [12], the nonperturbative 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

$$B^{(s)}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = \int d\mathbf{r}_{13} d\mathbf{r}_{23} e^{i(\mathbf{k}_1 \cdot \mathbf{r}_{13} + \mathbf{k}_2 \cdot \mathbf{r}_{23})} \langle A_1 A_2 A_3 e^{j_4 A_4 + j_5 A_5} \rangle, \quad (21)$$

where the quantities $A_i$, $j_i$ are, respectively, defined by

$$A_1 = \delta(\mathbf{r}_1) + f \nabla_z u_z(\mathbf{r}_1), \quad (22)$$

$$A_2 = \delta(\mathbf{r}_2) + f \nabla_z u_z(\mathbf{r}_2), \quad (23)$$

$$A_3 = \delta(\mathbf{r}_3) + f \nabla_z u_z(\mathbf{r}_3), \quad (24)$$

$$A_4 = u_z(\mathbf{r}_1) - u_z(\mathbf{r}_3), \quad (25)$$

$$A_5 = u_z(\mathbf{r}_2) - u_z(\mathbf{r}_3), \quad (26)$$

$$j_4 = -ik_1 \mu_1 f, \quad (27)$$

$$j_5 = -ik_2 \mu_2 f, \quad (28)$$

with $\mu_i = \mathbf{k}_i \cdot \hat{\mathbf{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., [4,12,37]). For the stochastic vector field $\mathbf{A}$, we have

$$\langle e^{j\mathbf{A}} \rangle = \exp \{\langle e^{j\mathbf{A}} \rangle_c\}, \quad (29)$$

with $j$ being an arbitrary constant vector, $\mathbf{j}$. To be specific, we assume that the vector fields given above are five components, i.e., $\mathbf{A} = \{A_1, A_2, A_3, A_4, A_5\}$ and $\mathbf{j} = \{j_1, j_2, j_3, j_4, j_5\}$. Then, taking the derivative 3 times with respect to $j_1, j_2$, and $j_3$, we set $j_1 = j_2 = j_3 = 0$. We obtain

$$\langle A_1 A_2 A_3 e^{j_4 A_4 + j_5 A_5} \rangle = \exp \{\langle e^{j_4 A_4 + j_5 A_5} \rangle_c \} [\langle A_1 A_2 A_3 e^{j_4 A_4 + j_5 A_5} \rangle_c + \langle A_1 A_2 e^{j_4 A_4 + j_5 A_5} \rangle_c \langle A_3 e^{j_4 A_4 + j_5 A_5} \rangle_c \\ + \langle A_2 e^{j_4 A_4 + j_5 A_5} \rangle_c \langle A_1 A_3 e^{j_4 A_4 + j_5 A_5} \rangle_c + \langle A_1 e^{j_4 A_4 + j_5 A_5} \rangle_c \langle A_2 A_3 e^{j_4 A_4 + j_5 A_5} \rangle_c \\ + \langle A_1 e^{j_4 A_4 + j_5 A_5} \rangle_c \langle A_2 e^{j_4 A_4 + j_5 A_5} \rangle_c \langle A_3 e^{j_4 A_4 + j_5 A_5} \rangle_c]. \quad (30)$$

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 nontrivial 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 the power spectrum [12]. 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 as 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 a nonperturbative part. On the other hand, the terms in the square brackets 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 imprint the acoustic feature of BAO. Although each term in the square brackets 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. [12] to derive the PT model of the redshift-space bispectrum valid at weakly nonlinear scales. That is,

(i) The overall prefactor, $\exp\{\langle e^{j_4 A_4 + j_5 A_5} \rangle_c\}$, is kept as a nonperturbative contribution and is replaced with general functional form $D_{\text{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_{\text{FoG}}$ will be discussed in Sec. IV.

(ii) The terms in the square brackets 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 proposition (i), the overall exponential factor is factorized outside the integral over $\mathbf{r}_{13}$ and $\mathbf{r}_{23}$. We have
$$B^{(s)}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) \rightarrow D_{\text{FoG}}(k_1\mu_1, k_2\mu_2, k_3\mu_3) \int dr_{13} dr_{23} \\ \times e^{i(\mathbf{k}_1 r_{13} + \mathbf{k}_2 r_{23})} [\langle A_1 A_2 A_3 e^{j_4 A_4 + j_5 A_5} \rangle_c \\ + \langle A_1 A_2 e^{j_4 A_4 + j_5 A_5} \rangle_c \langle A_3 e^{j_4 A_4 + j_5 A_5} \rangle_c \\ + \langle A_2 e^{j_4 A_4 + j_5 A_5} \rangle_c \langle A_1 A_3 e^{j_4 A_4 + j_5 A_5} \rangle_c \\ + \langle A_1 e^{j_4 A_4 + j_5 A_5} \rangle_c \langle A_2 A_3 e^{j_4 A_4 + j_5 A_5} \rangle_c \\ + \langle A_1 e^{j_4 A_4 + j_5 A_5} \rangle_c \langle A_2 e^{j_4 A_4 + j_5 A_5} \rangle_c \\ \times \langle A_3 e^{j_4 A_4 + j_5 A_5} \rangle_c]. \quad (31)$$

We then expand the terms in the square brackets. Up to the third order in $j_n$, we obtain

$$B^{(s)}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) \rightarrow D_{\text{FoG}}(k_1\mu_1, k_2\mu_2, k_3\mu_3) \sum_{n=1}^{11} C_n(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3), \quad (32)$$

where the functions $C_n$ are defined by

$$C_n(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) \equiv \int dr_{13} dr_{23} e^{i(\mathbf{k}_1 r_{13} + \mathbf{k}_2 r_{23})} S_n \quad (33)$$

with the integrands $S_n$ given below,

$$S_1 = \langle A_1 A_2 A_3 \rangle_c, \quad (34)$$

$$S_2 = \langle A_1 A_2 \rangle_c \langle (j_4 A_4 + j_5 A_5) A_3 \rangle_c + \text{cyc}, \quad (35)$$

$$S_3 = \langle (j_4 A_4 + j_5 A_5) A_1 A_2 A_3 \rangle_c, \quad (36)$$

$$S_4 = \langle (j_4 A_4 + j_5 A_5) A_1 A_2 \rangle_c \langle (j_4 A_4 + j_5 A_5) A_3 \rangle_c + \text{cyc}, \quad (37)$$

$$S_5 = \langle A_1 A_2 \rangle_c \langle (j_4 A_4 + j_5 A_5)^2 A_3 \rangle_c + \text{cyc}, \quad (38)$$

$$S_6 = \langle (j_4 A_4 + j_5 A_5) A_1 \rangle \langle (j_4 A_4 + j_5 A_5) A_2 \rangle \\ \times \langle (j_4 A_4 + j_5 A_5) A_3 \rangle, \quad (39)$$

$$S_7 = \langle (j_4 A_4 + j_5 A_5)^2 A_1 A_2 A_3 \rangle_c, \quad (40)$$

$$S_8 = \langle (j_4 A_4 + j_5 A_5)^3 A_1 A_2 A_3 \rangle_c, \quad (41)$$

$$S_9 = \langle (j_4 A_4 + j_5 A_5)^2 A_1 A_2 \rangle_c \langle (j_4 A_4 + j_5 A_5) A_3 \rangle_c + \text{cyc}, \quad (42)$$

$$S_{10} = \langle (j_4 A_4 + j_5 A_5) A_1 A_2 \rangle_c \langle (j_4 A_4 + j_5 A_5)^2 A_3 \rangle_c + \text{cyc}, \quad (43)$$

$$S_{11} = \langle A_1 A_2 \rangle_c \langle (j_4 A_4 + j_5 A_5)^3 A_3 \rangle_c + \text{cyc}. \quad (44)$$

In the above, the relevant terms for one-loop PT calculations, which are of the order of $O(P_L^3)$, appear at $n \le 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

$$B_{\text{model}}^{(s)}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = D_{\text{FoG}}(k_1\mu_1, k_2\mu_2, k_3\mu_3) \\ \times \sum_{n=1}^{6} C_n(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3). \quad (45)$$

The above model is compared with the standard PT results in Sec. II B, which can also be 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

$$B_{\text{SPT,1-loop}}^{(s)}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = \sum_{n=1}^{6} C_n(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) \\ + \sum_{n=1}^{2} D_n(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) \quad (46)$$

with the functions $D_n$ given by

$$D_n(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = \frac{1}{2} \int dr_{13} dr_{23} e^{i(\mathbf{k}_1 r_{13} + \mathbf{k}_2 r_{23})} \\ \times S_n \langle (j_4 A_4 + j_5 A_5)^2 \rangle_c, \quad (47)$$

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\{\langle e^{j_4 A_4 + j_5 A_5} \rangle_c\}$. Hence, at one-loop order, Eq. (45) is recast as

$$B_{\text{model}}^{(s)}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = D_{\text{FoG}}(k_1\mu_1, k_2\mu_2, k_3\mu_3) \\ \times \left\{ B_{\text{SPT,1-loop}}^{(s)}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) \right. \\ \left. - \sum_{n=1}^{2} D_n(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) \right\}. \quad (48)$$

In what follows, we use Eq. (48) to compute the PT model of the redshift-space bispectrum and compare the predictions with N-body simulations. To be precise, we first compute $B_{\text{SPT,1-loop}}^{(s)}$ based on the standard PT calculations summarized in Sec. II B [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 the 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. 1, ignoring the $D_{\text{FoG}}$
<figure><img src="image_3.png" /><figcaption>FIG. 1. Monopole (left) and quadrupole (right) moments of the redshift-space bispectrum from PT. The results at $z=1$ are shown. While top panels show the scale dependence of the equilateral configuration (i.e., plotted as a function of $k = k_1 = k_2 = k_3$ fixing $\theta_{12}$ to $2\pi/3$), bottom panels plot the shape dependence for isosceles configuration fixing the wave numbers to $k1 = k2 = k = 0.19h \text{ Mpc}^{-1}$. In each panel, green dotted and red dashed lines are, respectively, the standard PT predictions at tree-level ($B_{\text{SPT,tree}}^{(s)}$) and one-loop order ($B_{\text{SPT,1-loop}}^{(s)}$). Magenta solid lines are the predictions based on Eq. (48), $B_{\text{model}}^{(s)}$, for which the term $D_{\text{FoG}}$ is set to 1. This is essentially $B_{\text{SPT,1-loop}}^{(s)}$, but the terms, $D_1$ (orange lines) and $D_2$ (cyan lines), are subtracted.</figcaption></figure>

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 a function of $k$ and angle $\theta_{12} \equiv \cos^{-1}(\hat{k}_1 \cdot \hat{k}_2)$, respectively. Here, the multipole moments of the bispectrum, $B_{\ell}^{(s)}$, are defined by

$$B_{\ell}^{(s)}(k_1, k_2, \theta_{12}) = \int_{0}^{2\pi} \frac{d\phi}{2\pi} \int_{0}^{1} d\mu B^{(s)}(k_1, k_2, k_3) P_{\ell}(\mu), \quad (49)$$

<figure><img src="image_4.png" /><figcaption>FIG. 2. Definition of the angles, $\omega$ and $\phi$, given in Eqs. (C4) and (C5), which characterize the orientation of the triangle against the line-of-sight direction, $\hat{z}$. Here, the unit vectors, $\hat{r}_1$ and $\hat{r}_2$, are expressed in terms of the quantities indicated in the figure by $\hat{r}_1 = \hat{z} \times (k_1 \times k_2)/(k_1 k_2 \sin \theta_{12} \sin \omega)$ and $\hat{r}_2 = k_1 \times k_2/(k_1 k_2 \sin \theta_{12})$.</figcaption></figure>

where the function $P_{\ell}(\mu)$ is the Legendre polynomial with $\mu$ being the directional cosine given by $\mu = \cos \omega$. The angles $\omega$ and $\phi$ characterize the orientation of the triangles (i.e., $k_1, k_2$, and $k_3$) with respect to the line-of-sight direction (see Fig. 2). 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. [33,34,45,49]), 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$.

Figure 1 shows that both the one-loop bispectra, $B_{\text{SPT,1-loop}}^{(s)}$ (red lines) and $B_{\text{model}}^{(s)}$ (magenta lines), 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 $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_{\text{SPT,1-loop}}^{(s)}$ indeed plays an important role to get a consistent damping behavior of $D_{\text{FoG}}$ from both power spectrum and bispectrum (see Sec. IV D).

## 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. IVA, we first compare the real-space results in Sec. IV B 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 C. Section IV D discusses the validity and consistency of the ansatz for the damping function $D_{\text{FoG}}$ in our PT model of bispectrum.

A. 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 [50].<fnref n="1" /> Each simulation started from a given realization of the initial matter density field. The initial conditions were generated with an improved version of MPGRAFIC [51] while the particles evolution were computed with an optimized version of the RAMSES N-body code [52]. The main characteristics of the two ensembles of simulations are summarized in Table I. Set B consists in 96 simulations with $1024^3$ particles in a cosmological volume of $(656.25 \ h^{-1} \text{Mpc})^3$. The total effective volume is $27 \ (h^{-1} \text{Gpc})^3$, and the mass resolution is $1.88 \times 10^{10} \ h^{-1} \text{M}_\odot$. We call this set high resolution (HR). Set C consists in 512 simulations with $512^3$ particles in a cosmological volume of $(1312.5 \ h^{-1} \text{Mpc})^3$. The total effective volume of $1158 \ (h^{-1} \text{Gpc})^3$ is larger but the mass resolution of $1.20 \times 10^{12} \ h^{-1} \text{M}_\odot$ is more limited. We call this set low resolution (LR). The two sets are complementary because they are affected at different levels by numerical effects such as sample variance, finite volume, and mass resolution effects (see Refs. [50,53] for a study of these effects).

The matter bispectrum is estimated using the BISP_MES code [54] kindly provided by 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 [54]. Particle positions and velocity from a given snapshot are provided as an input of the code. Particle positions are then displaced along the $\hat{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 jth simulation of a set, $\delta_{\text{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$.

<table><caption>TABLE I. Parameters of N-body simulation data sets used in the paper.</caption><thead><tr><th></th><th>HR</th><th>LR</th></tr></thead><tbody><tr><td>Set</td><td>B</td><td>C</td></tr><tr><td>Number of realizations</td><td>96</td><td>512</td></tr><tr><td>Box size [$h^{-1}$ Mpc]</td><td>656.25</td><td>1312.5</td></tr><tr><td>Number of particles</td><td>1024<sup>3</sup></td><td>512<sup>3</sup></td></tr><tr><td>Particle mass [$h^{-1}$M<sub>☉</sub>]</td><td>1.88 × 10<sup>10</sup> h<sup>-1</sup>M<sub>☉</sub></td><td>1.20 × 10<sup>12</sup> h<sup>-1</sup>M<sub>☉</sub></td></tr><tr><td>Initial redshift ($z_{\text{ini}}$)</td><td>190</td><td>106</td></tr><tr><td>Output redshifts</td><td>$z = 1, 0$</td><td>$z = 1, 0.5, 0$</td></tr></tbody></table>

The averaging procedure is performed by randomly picking $N_{\text{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 jth simulation is then given by

$$B_{\ell, \text{sim,j}}^{(s)}(k_1, k_2, k_3) = \frac{1}{N_{\text{mode}}} \times \sum_{i=1}^{N_{\text{mode}}} \delta_{\text{sim,j}}^{(s)}(k_1^i) \delta_{\text{sim,j}}^{(s)}(k_2^i) \delta_{\text{sim,j}}^{(s)}(k_3^i) P_{\ell}(\mu^i), \quad (50)$$

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

The number of orientations/sizes $N_{\text{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_{\text{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 \ \text{Mpc}^{-1}$ ($1.23 \ h \ \text{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 \ \text{Mpc}^{-1}$). However, at the very large scales of interest for comparison to perturbation theory ($k = 0.01-0.1 \ h \ \text{Mpc}^{-1}$), HR simulations might suffer from non-negligible finite-mode sampling [55]. In the linear regime, the bispectrum of the density field should be zero for a perfect ensemble average. Because the total number of independent modes in the simulation is finite, the resulting bispectrum is nonzero. To mitigate this effect at large scales, we compute the initial bispectrum of each snapshot $B_{\ell, \text{ini,j}}^{(s)}$. The linearly evolved bispectrum is then subtracted from the snapshot bispectrum to obtain a corrected bispectrum $B_{\ell, \text{corr,j}}^{(s)}$. The linearly evolved bispectrum is computed using Eq. (50), but instead of using the snapshot density $\delta_{\text{sim,j}}^{(s)}(k^i)$, we use the linearly evolved density field

$$\delta_{\text{linevol,j}}^{(s)}(\mathbf{k}^i, z) = \frac{1 + f(z)\mu^2}{1 + f(z_{\text{ini}})\mu^2} \frac{D(z)}{D(z_{\text{ini}})} \times \delta_{\text{sim,j}}^{(s)}(\mathbf{k}^i, z_{\text{ini}}), \quad (51)$$

<footnote n="1">¹Set C was not directly presented in [50] but it was performed at the same time for the same project.</footnote>
<figure><img src="image_5.png" /><figcaption>FIG. 3. Real-space power spectrum as a function of wave number $k$ at $z = 1$ (top), 0.5 (middle), and 1 (bottom). The left panel plots the power spectrum multiplied by $k^{3/2}$, while the right panel shows the fractional difference between the emulator and other predictions, $(P - P_{\text{Emu}})/P_{\text{Emu}}$. In each panel, black and gray points with error bars indicate the results of LR and HR simulations, respectively. The green dashed and red solid lines are the linear theory and the standard PT one-loop predictions. The black dashed lines in the left panels are the predictions based on the cosmic emulator [56–59].</figcaption></figure>

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 \text{ 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_{\ell, \text{corr}, j}^{(s)}(k_1, k_2, k_3)$ to get the mean bispectrum of the set

$$B_{\ell, \text{sim}}^{(s)}(k_1, k_2, k_3) = \frac{1}{N_{\text{sim}}} \sum_{j=1}^{N_{\text{sim}}} B_{\ell, \text{corr}, j}^{(s)}(k_1, k_2, k_3), \quad (52)$$

where $N_{\text{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 error of the bispectrum, and we estimate statistical error bars assuming<fnref n="2" />

$$\Delta B_{\ell, \text{stat}}^{(s)} = \frac{\sqrt{\frac{1}{N_{\text{sim}}} \sum_{j=1}^{N_{\text{sim}}} B_{\ell, \text{corr}, j}^{(s)}(k_1, k_2, k_3)^2 - B_{\ell, \text{sim}}^{(s)}(k_1, k_2, k_3)^2}}{\sqrt{N_{\text{sim}}}} \quad (53)$$

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.

## B. 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.

Figure 3 shows the power spectra at $z = 1$ (top), 0.5 (middle), and 0 (bottom). The 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 [56–59] as the base value and evaluate $[P(k) - P_{\text{Emu}}(k)]/P_{\text{Emu}}(k)$ with $P_{\text{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 lines) reasonably agree with those of the emulator code, and the agreement is at the 3% level for the scales of $k \lesssim 0.18$, 0.12, and $0.1 \text{ h Mpc}^{-1}$ at $z = 1$, 0.5, and 0, 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% level over all the studied ranges of wave number) and one-loop PT predictions, but the low-resolution data (LR) systematically deviate from others at small scales. The deviation gradually increases from about 0.5% near $k = 0.05 \text{ h Mpc}^{-1}$ to 5%

<footnote n="2">²Multiplying the expression of Eq. (53) by $\sqrt{N_{\text{sim}}/(N_{\text{sim}} - 1)}$ decreases the bias in the estimate of the standard error. However, we think that our realizations are enough to neglect $-1$.</footnote>
<figure><img src="image_6.png" /></figure>

FIG. 4. Real-space bispectrum in real space at $z = 1$ (top), 0.5 (middle), and 0 (bottom). The first two panels from the left show the results for equilateral triangular shape (i.e., $k_1 = k_2 = k_3 = k$), plotted as a function of wave number $k$. The leftmost panel is the bispectrum amplitude multiplied by $k^3$, while the second panel plots the fractional difference between standard PT one-loop predictions and others, i.e., $(B - B_{1\text{-loop}})/B_{1\text{-loop}}$. The right two panels also present the results similar to the left two panels, but the cases with isosceles triangular shape ($k_1 = k_2$) are plotted as a function of the angle $\theta_{12} \equiv \cos^{-1}(\hat{k}_1 \cdot \hat{k}_2)$, fixing the scale $k_1 = k_2$ to 0.096 $h$ Mpc$^{-1}$. In each panel, meanings of the symbols and line types are the same as in Fig. 3, except for the green dashed lines, which indicate the standard PT prediction at tree-level order.

near $k = 0.2 h$ Mpc$^{-1}$. This is a well known mass-resolution effect which tends to decrease the power at a small scale [53, 56, 60]. We refer to Ref. [53] for a dedicated section about statistical and systematics errors on the matter power spectrum with a similar simulation setup as the one used in this paper.

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. Figure 4 shows the results for the equilateral triangular shapes plotted as a function of $k$ (left two panels), and those for the isosceles triangles plotted as a function of $\theta_{12} \equiv \cos^{-1}(\hat{k}_1 \cdot \hat{k}_2)$ (right two panels). Note that in panels showing the amplitude of the 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_{1\text{-loop}})/B_{1\text{-loop}}$.

Compared to the power spectrum case, the statistical errors in N-body simulations are larger, and the difference between the two data sets 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 looks somewhat worse, and is at the 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 nonlinear 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 the next subsection.

### C. 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.

Figures 5 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 from both the LR and the 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 error bars 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
<figure><img src="image_7.png" /><figcaption>FIG. 5. Monopole (left) and quadrupole (right) moments of the redshift-space bispectrum for the equilateral triangles, plotted as a function of wave number $k$ at $z = 1$ (top), 0.5 (middle), and 0 (bottom). The standard PT predictions at tree-level (green dashed line) and one-loop (red solid line) order are compared with the measured bispectrum from LR (black lines) and HR (gray lines) data of N-body simulations. The results for equilateral shape triangles are particularly shown.</figcaption></figure>

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

<figure><img src="image_8.png" /><figcaption>FIG. 6. Same as Fig. 5, but the shape dependence of the bispectrum is plotted as a function of $\theta_{12} = \cos^{-1}(\hat{k}_1 \cdot \hat{k}_2)$, fixing both $k_1$ and $k_2$ to 0.096 h Mpc$^{-1}$. The meanings of the line types are the same as in Fig. 5.</figcaption></figure>
$D_{FoG}$ needs to be first specified. The function $D_{FoG}$ is generally expressed as a function of $k_1\mu_1$, $k_2\mu_2$, and $k_3\mu_3$, and it contains the nonperturbative damping behavior arising from the exponential factor, $\exp\{\langle e^{j4A_4+j5A_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 [12]:

$$D_{FoG} = \exp \left[ -\frac{1}{2} (f \sigma_v)^2 \{(k_1 \mu_1)^2 + (k_2 \mu_2)^2 + (k_3 \mu_3)^2\} \right]. \quad (54)$$

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 are obtained by expanding the exponential prefactor and truncating it at the second-order in cumulants, just ignoring the spatial correlation. A nontrivial point may be whether $D_{FoG}$ is still expressed as a 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 D.

Adopting Eq. (54), Figs. 7–10 compare the prediction of the PT model with N-body simulations. While Figs. 7 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. 9 and 10, respectively, plot the shape dependence of the bispectrum as a function of $\theta_{12} = \cos^{-1}(\hat{k}_1 \cdot \hat{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_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_{min}, k_{max}$]. While $k_{min}$ is set to $0.05 \text{ h Mpc}^{-1}$, we adopt $k_{1\%}$ defined by Ref. [61] as the maximum wave number $k_{max}$, indicated by the vertical arrow in each panel of Figs. 7 and 8 (blue for tree-level PT and red for one-loop PT). The $k_{1\%}$ indicates the maximum wave number 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

<figure><img src="image_9.png" /><figcaption>FIG. 7. Monopole (left) and quadrupole (right) moments of the redshift-space bispectrum for the equilateral triangles, plotted as a function of wave number $k$ at $z = 1$ (top), 0.5 (middle), and 0 (bottom). Here, the PT models involving the damping function $D_{FoG}$ are compared with the measured results obtained from the LR data of N-body simulations. The predictions of tree- and one-loop PT models are depicted as blue dashed and magenta solid lines, respectively. In plotting the predictions, the free parameter $\sigma_v$ in the damping function is determined by fitting the monopole and quadrupole moments to the N-body data at $0.05 \text{ h Mpc}^{-1} \le k \le k_{max}$, with $k_{max}$ indicated by the vertical arrows in each panel. Note that the error bars depicted as solid and dotted lines, respectively, represent the statistical error averaged over the number of realizations and the one including both the statistical and systematic errors (see Sec. IV D).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 8. Same as Fig. 7, but the results for the scalene triangle of $k_1 = 2k_2 = (2/\sqrt{3})k_3 = k$ are shown.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 9. Same as Fig. 7, but here the shape dependence of the bispectrum for isosceles triangle ($k_1 = k_2$) is plotted as a function of $\theta_{12} = \cos^{-1}(\hat{k}_1 \cdot \hat{k}_2)$, fixing the scale of $k_1$ and $k_2$ to 0.096 h Mpc$^{-1}$. Here, in plotting the predictions of PT models, we adopt the fitted results of $\sigma_v$ obtained from Fig. 7.</figcaption></figure>

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_v$ is given in the left panel of each figure, which is close to the linear theory prediction of one-dimensional velocity dispersion (see also Table II and Figs. 11 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_v$. For comparison, in Figs. 7–10, we have also plotted the tree-level standard PT predictions multiplied by the damping function of Eq. (54),<fnref n="3" /> depicted as blue short-dashed lines, but the agreement with N-body simulations is

<footnote n="3">³In terms of the descriptions given in Sec. III, the tree-level standard PT multiplied by $D_{\text{FOG}}$ corresponds to the leading-order PT calculations of Eq. (45) with the functions $C_n$ summing up to $n=2$.</footnote>
<figure><img src="image_12.png" /><figcaption>FIG. 10. Same as Fig. 9, but the results for another triangular shape with $k_1 = 2k_2$ are shown, again fixing $k_1$ and $2k_2$ to 0.096 h Mpc$^{-1}$. Here, in plotting the predictions of PT models, we adopt the fitted results of $\sigma_v$ obtained from Fig. 8.</figcaption></figure>

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 the weakly non-linear regime. A closer look at the equilateral case in Fig. 7 reveals that the agreement of the one-loop model is a bit degraded compared to the scalene case (Fig. 8), especially for the quadrupole moment. A similar trend is also found for the shape dependence in Figs. 9 and 10. This is partly because for a fixed wave number $k$, one of the side lengths 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.

### D. Testing the ansatz of damping function

Adopting the Gaussian form of the nonperturbative 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_v$. Since the present PT model is constructed in a manner similar to the model of power spectrum in Ref. [12], the fitted values of the parameter $\sigma_v$ derived from the 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 has to be checked at the scales relevant for the one-loop PT.

Let us first discuss the consistency of the fitted value, $\sigma_v$. Table II summarizes the results derived from both the bispectrum and the 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 the power spectrum, while for the bispectrum, the fitted results from two different configurations are presented in both tree-level

<table><caption>TABLE II. Fitted values of $\sigma_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-body simulations and are listed in units of h Mpc$^{-1}$. For reference, the rightmost column shows the linear theory predictions.</caption><thead><tr><th rowspan="2"></th><th colspan="2">k<sub>max</sub> [h Mpc<sup>-1</sup>]</th><th rowspan="2">P<sup>(s)</sup><sub>model</sub>: 1-loop</th><th colspan="2">B<sup>(s)</sup><sub>model</sub>: 1-loop</th><th colspan="2">B<sup>(s)</sup><sub>model</sub>: tree</th><th rowspan="2">σ<sub>v,lin</sub></th></tr><tr><th>1-loop</th><th>Tree</th><th>k<sub>1</sub> = k<sub>2</sub> = k<sub>3</sub></th><th>k<sub>1</sub> = 2k<sub>2</sub> = (2/√3)k<sub>3</sub></th><th>k<sub>1</sub> = k<sub>2</sub> = k<sub>3</sub></th><th>k<sub>1</sub> = 2k<sub>2</sub> = (2/√3)k<sub>3</sub></th></tr></thead><tbody><tr><td>z = 1</td><td>0.13</td><td>0.09</td><td>4.3 ± 0.1 (4.2 ± 0.04)</td><td>4.9 ± 0.3 (5.9 ± 0.3)</td><td>4.0 ± 0.2 (4.3 ± 0.3)</td><td>6.1 ± 0.2</td><td>6.5 ± 0.1</td><td>3.8</td></tr><tr><td>z = 0.5</td><td>0.1</td><td>0.07</td><td>5.9 ± 0.1</td><td>7.8 ± 0.4</td><td>6.3 ± 0.3</td><td>9.0 ± 0.3</td><td>9.5 ± 0.2</td><td>4.8</td></tr><tr><td>z = 0</td><td>0.08</td><td>0.06</td><td>8.8 ± 0.3 (9.8 ± 0.3)</td><td>13 ± 0.8 (15 ± 1.0)</td><td>10 ± 0.6 (9.1 ± 1.1)</td><td>13 ± 0.6</td><td>15 ± 0.4</td><td>6.1</td></tr></tbody></table>
<figure><img src="image_13.png" /><figcaption>FIG. 11. Dependence of the fitted results of parameter $\sigma_v$ (top) and the reduced chi-square (bottom) on the maximum wave number $k_{\text{max}}$ at $z=1$ (left), 0.5 (middle), and 0 (right). Red and green symbols, respectively, represent the results from the bispectrum for equilateral and scalene shapes with $k_1 = 2k_2 = (2/\sqrt{3})k_3$. On the other hand, blue symbols are obtained from the power spectrum, based on the PT model by Ref. [12]. Note that all the results are obtained by fitting the one-loop PT model predictions to the measured results from LR data of N-body simulations.</figcaption></figure>

and one-loop PT cases. The best-fitted values from both the power spectrum and the 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 the bispectrum.

To see the robustness of the fitted values, we extend the analysis in Table II, and using the LR data, we examine the fitting in various ranges of $k$. The results are shown in Figs. 11 and 12. In each figure, top and bottom panels, respectively, plot the fitted values $\sigma_v$ and the reduced $\chi^2$ as a function of $k_{\text{max}}$. Note that in estimating $\chi^2_{\text{red}}$ and the errors

of the fitted $\sigma_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,\text{sys}}^{(s)}$ and $\Delta B_{\ell,\text{sys}}^{(s)}$ to the statistical errors of the power spectrum and bispectrum multipoles, $\Delta P_{\ell,\text{stat}}^{(s)}$ and $\Delta B_{\ell,\text{stat}}^{(s)}$, as shown in Figs. 7–10 for the bispectrum case (solid error bars: $\Delta B_{\ell,\text{stat}}^{(s)}$; dotted error bars: $\Delta B_{\ell,\text{stat}}^{(s)} + \Delta B_{\ell,\text{sys}}^{(s)}$). The size of the systematics is estimated from the real-space results in Figs. 3 and 4. We assume that the relative systematics in redshift space for the power spectrum and bispectrum are proportional to the

<figure><img src="image_14.png" /><figcaption>FIG. 12. Same as in Fig. 11, but the results obtained by fitting the tree-level PT model predictions to the N-body data are shown. Meanings of the lines and symbols are the same as in Fig. 11.</figcaption></figure>
relative systematics in real space for the matter power spectrum, $\Delta P_{\text{sys}}/P_{\text{sim}} = |P_{\text{sim}} - P_{\text{ref}}|/P_{\text{sim}}$, where $P_{\text{sim}}$ represents the measured real-space power spectrum in N-body simulations, and $P_{\text{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%, we neglect this contribution). The systematics for the power spectrum and bispectrum in redshift space are estimated as

$$ \frac{\Delta P_{\ell, \text{sys}}^{(s)}(k_a)}{P_{\ell, \text{sim}}^{(s)}(k_a)} = \alpha \frac{\Delta P_{\text{sys}}(k_a)}{P_{\text{sim}}(k_a)}, \qquad (55) $$

$$ \frac{\Delta B_{\ell, \text{sys}}^{(s)}(k_1, k_2, k_3)}{B_{\ell, \text{sim}}^{(s)}(k_1, k_2, k_3)} = \beta \sqrt{\frac{1}{3} \sum_{a=1}^{3} \left( \frac{\Delta P_{\text{sys}}(k_a)}{P_{\text{sim}}(k_a)} \right)^2}, \qquad (56) $$

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 the 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^2_{\text{red}}$ can only be used for a comparison between tree-level and one-loop results in Figs. 11 and 12. Nevertheless, we see that the fitting results in one-loop PT cases are basically stable against the variation of $k_{\text{max}}$, and the

estimated values of $\chi^2_{\text{red}}$ are smaller than those in the tree-level cases. Further, we checked that the (best-)fitted values of $\sigma_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 probably be due to the breakdown of the one-loop predictions. In fact, the $\chi^2_{\text{red}}$ systematically increases with $k_{\text{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_{\text{FoG}}$ at Eq. (54). To clarify whether the nonperturbative part is described by the univariate function, 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_{\text{FoG}}$ to 1. Taking the ratio gives

$$ \frac{B_{\text{sim}}^{(s)}(k_1, k_2, \theta_{12}, \omega, \phi)}{B_{\text{model}}^{(s)}(k_1, k_2, \theta_{12}, \omega, \phi)|_{D_{\text{FoG}}=1}}. \qquad (57) $$

At the scales where the one-loop PT is applicable, this ratio directly quantifies the functional form of $D_{\text{FoG}}$, and thus we can check whether it is expressed as a univariate function of $(k_1\mu_1)^2 + (k_2\mu_2)^2 + (k_3\mu_3)^2$.

Figure 13 shows the measured results of the ratio, Eq. (57), for various orientations at $z = 1$, plotted against

<figure><img src="image_15.png" /><figcaption>FIG. 13. Ratio of the redshift-space bispectrum measured from the LR data and to the one-loop (top) and tree-level (bottom) PT models, setting $D_{\text{FoG}}$ to unity (or equivalently, taking the limit $\sigma_v \to 0$). Here, without applying multipole expansion, the ratio is evaluated for various orientations of the triangle, and the scatters of the estimated results are shown as color shaded regions, fixing the equilateral shape (left) and the scalene shape with $k_1 = 2k_2 = (2/\sqrt{3})k_3$ (right). The resultant ratio is then plotted against $(k_1\mu_1)^2 + (k_2\mu_2)^2 + (k_3\mu_3)^2$. Note that the four different color strengths imply the length of $k_1$ used to estimate the ratio, indicated at the upper horizontal axis. For reference, Gaussian (dashed lines) and Lorentzian (dotted lines) forms of the damping function are also plotted.</figcaption></figure>
$(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 the shaded color region, and the four different color strengths imply the length of $k_1$ used to estimate the ratio: $k_1 \le 0.05, 0.15, 0.5$, and $1 \text{ h Mpc}^{-1}$ from dark to light. The results are compared with the univariate damping function of the Gaussian (dashed lines) and Lorentzian (dotted lines) form. Clearly, the scatter of the ratio for the one-loop PT model, given in the 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 the 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 the Lorentzian form describes the measured ratio reasonably well at the high-$k$ tail, although it is mostly the boundary where we can apply the one-loop PT prediction. Hence, we conclude that the univariate ansatz for $D_{\text{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 PT model that can keep the nonperturbative damping effect of the RSD under control. Starting with the exact formula for the redshift-space bispectrum, we rewrite the expression in terms of the cumulants to identify the nonperturbative term. Separating the nonperturbative term responsible for the so-called FoG damping effect, we derive the perturbative expressions for bispectrum valid at one-loop order. The resultant model has been constructed similar to the power spectrum model in Ref. [12], and it incorporates the nonperturbative damping term on top of the terms that can be computed with standard PT.

Adopting the Gaussian form of the damping function, we have performed a detailed comparison between the predictions of the 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_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_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_{\text{FoG}}$ indeed holds for the 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 the one-loop correction is essential, and together with the model of power spectrum proposed by Ref. [12], 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 tighter test of gravity will be made. The effective field theory approach may also help to improve the prediction (e.g., Refs. [62–65]). Another issue is the estimation of statistical error covariance of the 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. [28,66] for the 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.

## ACKNOWLEDGMENTS

This work is in part supported by MEXT/JSPS KAKENHI Grants No. JP15H05899 and No. JP16H03977. It is also supported by JSPS Grant No. 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, as well as Takahiro Nishimichi and Shun Saito for discussion and useful comments.
## 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_L$. The kernels up to fourth order become (e.g., [34])

$$Z_1(\mathbf{k}) = 1 + f\mu^2, \qquad (A1)$$

$$Z_2(\mathbf{k}_1, \mathbf{k}_2) = F_2(\mathbf{k}_1, \mathbf{k}_2) + f\mu^2 G_2(\mathbf{k}_1, \mathbf{k}_2) + \frac{f\mu k}{2} \left[ \frac{\mu_1}{k_1} (1 + f\mu_2^2) + \frac{\mu_2}{k_2} (1 + f\mu_1^2) \right], \qquad (A2)$$

$$Z_3(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = F_3(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) + f\mu^2 G_3(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) + f\mu k [F_2(\mathbf{k}_2, \mathbf{k}_3) + f\mu_{23}^2 G_2(\mathbf{k}_2, \mathbf{k}_3)] \frac{\mu_1}{k_1} \\ + f\mu k (1 + f\mu_1^2) \frac{\mu_{23}}{k_{23}} G_2(\mathbf{k}_2, \mathbf{k}_3) + \frac{(f\mu k)^2}{2} (1 + f\mu_1^2) \frac{\mu_2 \mu_3}{k_2 k_3}, \qquad (A3)$$

$$Z_4(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3, \mathbf{k}_4) = Z_{4a}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3, \mathbf{k}_4) + Z_{4b}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3, \mathbf{k}_4). \qquad (A4)$$

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

$$Z_{4a}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3, \mathbf{k}_4) = F_4(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3, \mathbf{k}_4) + f\mu^2 G_4(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3, \mathbf{k}_4) + f\mu k [F_3(\mathbf{k}_2, \mathbf{k}_3, \mathbf{k}_4) + f\mu_{23}^2 G_3(\mathbf{k}_2, \mathbf{k}_3, \mathbf{k}_4)] \frac{\mu_1}{k_1} \\ + f\mu k (1 + f\mu_1^2) \frac{\mu_{234}}{k_{234}} G_3(\mathbf{k}_2, \mathbf{k}_3, \mathbf{k}_4) + \frac{(f\mu k)^3}{6} (1 + f\mu_1^2) \frac{\mu_2 \mu_3 \mu_4}{k_2 k_3 k_4}, \qquad (A5)$$

$$Z_{4b}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3, \mathbf{k}_4) = f\mu k [F_2(\mathbf{k}_1, \mathbf{k}_2) + f\mu_{12}^2 G_2(\mathbf{k}_1, \mathbf{k}_2)] \frac{\mu_{34}}{k_{34}} G_2(\mathbf{k}_3, \mathbf{k}_4) \\ + \frac{(f\mu k)^2}{2} [F_2(\mathbf{k}_1, \mathbf{k}_2) + f\mu_{12}^2 G_2(\mathbf{k}_1, \mathbf{k}_2)] \frac{\mu_3 \mu_4}{k_3 k_4} + \frac{(f\mu k)^2}{2} \frac{\mu_{12}}{k_{12}} G_2(\mathbf{k}_1, \mathbf{k}_2) \left[ (1 + f\mu_3^2) \frac{\mu_4}{k_4} + (1 + f\mu_4^2) \frac{\mu_3}{k_3} \right]. \qquad (A6)$$

In the above, the vector $\mathbf{k}$ in the nth order kernel implies $\mathbf{k} = \mathbf{k}_1 + \cdots + \mathbf{k}_n$. The quantities $\mu$, $\mu_i$, $\mu_{ij}$, and $\mu_{ijk}$ are defined by

$$\mu \equiv \frac{\mathbf{k} \cdot \hat{\mathbf{z}}}{\mathbf{k}}, \qquad \mu_i \equiv \frac{\mathbf{k}_i \cdot \hat{\mathbf{z}}}{\mathbf{k}_i}, \qquad \mu_{ij} \equiv \frac{(\mathbf{k}_i + \mathbf{k}_j) \cdot \hat{\mathbf{z}}}{|\mathbf{k}_i + \mathbf{k}_j|}, \qquad \mu_{ijk} \equiv \frac{(\mathbf{k}_i + \mathbf{k}_j + \mathbf{k}_k) \cdot \hat{\mathbf{z}}}{|\mathbf{k}_i + \mathbf{k}_j + \mathbf{k}_k|}. \qquad (A7)$$

In applying these results to the statistical calculations described in Sec. II B, 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 \ge 3$ only preserve partial symmetry. For instance, while the kernel $Z_3$ given at Eq. (A3) is symmetric under $\mathbf{k}_1 \leftrightarrow \mathbf{k}_3$, the expression for $Z_{4a}$ preserves the symmetry of $\mathbf{k}_2 \leftrightarrow \mathbf{k}_3 \leftrightarrow \mathbf{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,

$$\Psi_a(\mathbf{k}) = (\delta(\mathbf{k}), \theta(\mathbf{k})), \qquad (B1)$$

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

$$(2\pi)^3 \delta_D(\mathbf{k}_1 + \mathbf{k}_2) P_{ab}(k_1) \equiv \langle \Psi_a(\mathbf{k}_1) \Psi_b(\mathbf{k}_2) \rangle, \qquad (B2)$$

$$(2\pi)^3 \delta_D(\mathbf{k}_1 + \mathbf{k}_2 + \mathbf{k}_3) B_{abc}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) \\ \equiv \langle \Psi_a(\mathbf{k}_1) \Psi_b(\mathbf{k}_2) \Psi_c(\mathbf{k}_3) \rangle, \qquad (B3)$$

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 $O(P_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,

$$P_{ab}(k) \approx P_L(k), \qquad (B4)$$

$$B_{abc}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) \approx 2 \left\{ \begin{aligned} & \mathcal{F}_a(\mathbf{k}_2, \mathbf{k}_3) P_L(\mathbf{k}_2) P_L(\mathbf{k}_3) \\ & + \mathcal{F}_b(\mathbf{k}_1, \mathbf{k}_3) P_L(\mathbf{k}_1) P_L(\mathbf{k}_3) \\ & + \mathcal{F}_c(\mathbf{k}_1, \mathbf{k}_2) P_L(\mathbf{k}_1) P_L(\mathbf{k}_2) \end{aligned} \right\}, \qquad (B5)$$

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:

$$\langle A_1 A_2 A_3 \rangle_c = \int \frac{d^3 \mathbf{p}_1 d^3 \mathbf{p}_2}{(2\pi)^6} e^{i(\mathbf{p}_1 \cdot \mathbf{r}_{13} + \mathbf{p}_2 \cdot \mathbf{r}_{23})} \sum_{a,b,c=1}^{2} f^{a+b+c-3} \times \mu_{p_1}^{2(a-1)} \mu_{p_2}^{2(b-1)} \mu_{p_3}^{2(c-1)} B_{abc}(\mathbf{p}_1, \mathbf{p}_2, \mathbf{p}_3), \qquad (B6)$$

$$\langle (j_4 A_4 + j_5 A_5)^2 \rangle_c = -f^2 \int \frac{d^3 \mathbf{p}}{(2\pi)^3} \frac{\mu_p^2}{p^2} P_{22}(p) \times \{ (k_1 \mu_1)^2 + (k_2 \mu_2)^2 + (k_3 \mu_3)^2 + 2(k_1 k_2 \mu_1 \mu_2 e^{i \mathbf{p} \cdot \mathbf{r}_{12}} + k_1 k_3 \mu_1 \mu_3 e^{i \mathbf{p} \cdot \mathbf{r}_{13}} + k_2 k_3 \mu_2 \mu_3 e^{i \mathbf{p} \cdot \mathbf{r}_{23}}) \}, \qquad (B7)$$

$$\langle A_1 A_2 \rangle_c = \int \frac{d^3 \mathbf{p}}{(2\pi)^3} e^{i \mathbf{p} \cdot \mathbf{r}_{12}} \times \{ P_{11}(p) + 2f \mu_p^2 P_{12}(p) + f^2 \mu_p^4 P_{22}(p) \}, \qquad (B8)$$

$$\langle (j_4 A_4 + j_5 A_5) A_3 \rangle_c = -f \int \frac{d^3 \mathbf{p}}{(2\pi)^3} \frac{\mu_p}{p} \times \{ P_{12}(p) + f \mu_p^2 P_{22}(p) \} \times \{ k_1 \mu_1 e^{i \mathbf{p} \cdot \mathbf{r}_{13}} + k_2 \mu_2 e^{i \mathbf{p} \cdot \mathbf{r}_{23}} + k_3 \mu_3 \} \qquad (B9)$$

with the quantities $\mu_i$ and $\mu_p$, respectively, defined by $\mu_i = (\mathbf{k}_i \cdot \hat{\mathbf{z}})/k_i$ and $\mu_p = (\mathbf{p} \cdot \hat{\mathbf{z}})/p$.

## 1. $D_1$ term

According to Eq. (47), the $D_1$ term is explicitly written as

$$D_1(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = \frac{1}{2} \int d^3 \mathbf{r}_{13} d^3 \mathbf{r}_{23} e^{i(\mathbf{k}_1 \cdot \mathbf{r}_{13} + \mathbf{k}_2 \cdot \mathbf{r}_{23})} \times \langle A_1 A_2 A_3 \rangle_c \langle (j_4 A_4 + j_5 A_5) \rangle_c^2. \qquad (B10)$$

Substituting Eq. (B6) and (B7) into the above, we obtain

$$D_1(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = -\frac{f^2}{2} \{(k_1\mu_1)^2 + (k_2\mu_2)^2 + (k_3\mu_3)^2\} \times \sigma_{v,\text{lin}}^2 C_1(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) - f^2 \int \frac{d^3\mathbf{p}}{(2\pi)^3} \frac{\mu_p^2}{p^2} P_{22}(p) \times \{(k_1 k_2 \mu_1 \mu_2) C_1(\mathbf{k}_1 + \mathbf{p}, \mathbf{k}_2 - \mathbf{p}, \mathbf{k}_3) + (k_1 k_3 \mu_1 \mu_3) C_1(\mathbf{k}_1 + \mathbf{p}, \mathbf{k}_2, \mathbf{k}_3 - \mathbf{p}) + (k_2 k_3 \mu_2 \mu_3) C_1(\mathbf{k}_1, \mathbf{k}_2 + \mathbf{p}, \mathbf{k}_3 - \mathbf{p})\} \qquad (B11)$$

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

$$\sigma_{v,\text{lin}}^2 = \int \frac{d^3 \mathbf{p}}{(2\pi)^3} \frac{\mu_p^2}{p^2} P_{22}(p) \approx \int \frac{dp}{6\pi^2} P_L(p). \qquad (B12)$$

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

$$C_1(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = \sum_{a,b,c=1}^{2} \mu_1^{2(a-1)} \mu_2^{2(b-1)} \mu_3^{2(c-1)} B_{abc}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3). \qquad (B13)$$

## 2. $D_2$ term

From Eq. (47), the $D_2$ term is expressed in terms of the cumulants below,

$$D_2(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = \frac{1}{2} \int d^3 \mathbf{r}_{13} d^3 \mathbf{r}_{23} e^{i(\mathbf{k}_1 \cdot \mathbf{r}_{13} + \mathbf{k}_2 \cdot \mathbf{r}_{23})} \times \{ \langle A_1 A_2 \rangle_c \langle (j_4 A_4 + j_5 A_5) A_3 \rangle_c + \langle A_1 A_3 \rangle_c \langle (j_4 A_4 + j_5 A_5) A_2 \rangle_c + \langle A_2 A_3 \rangle_c \langle (j_4 A_4 + j_5 A_5) A_1 \rangle_c \} \times \langle (j_4 A_4 + j_5 A_5)^2 \rangle_c. \qquad (B14)$$

Substituting Eqs. (B7)–(B9) into the above expression, after lengthy calculation, we obtain
$$D_2(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = -\frac{f^2}{2} \{(k_1\mu_1)^2 + (k_2\mu_2)^2 + (k_3\mu_3)^2\} \\ \times \sigma_{\text{v,lin}}^2 C_2(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) \\ + f^3(k_1\mu_1)(k_2\mu_2)(k_3\mu_3)J_2(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) \\ + f^3 K_2(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3). \quad (B15)$$

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,

$$C_2(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = -f(k_1\mu_1)\{X(\mathbf{k}_2)Y(\mathbf{k}_3) + X(\mathbf{k}_3)Y(\mathbf{k}_2)\} \\ + \text{cyc}, \quad (B16)$$

$$J_2(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = \int \frac{d^3\mathbf{p}}{(2\pi)^3} \frac{\mu_p^2}{p^2} P_{22}(p) \{X(\mathbf{k}_2 - \mathbf{p})Y(\mathbf{k}_1 + \mathbf{p}) \\ + X(\mathbf{k}_1 + \mathbf{p})Y(\mathbf{k}_2 - \mathbf{p})\} + \text{cyc}, \quad (B17)$$

$$K_2(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = \int \frac{d^3\mathbf{p}}{(2\pi)^3} \frac{\mu_p^2}{p^2} P_{22}(p) \\ \times \left[ \{(k_1\mu_1)^2(k_2\mu_2)X(\mathbf{k}_2 - \mathbf{p}) \\ + (k_1\mu_1)(k_2\mu_2)^2X(\mathbf{k}_1 - \mathbf{p})\}Y(\mathbf{k}_3) \\ + (X \leftrightarrow Y) \right] + \text{cyc}, \quad (B18)$$

with the functions X and Y defined by

$$X(\mathbf{p}) = P_{11}(p) + 2f\mu_p^2 P_{12}(p) + f^2\mu_p^4 P_{22}(p), \quad (B19)$$

$$Y(\mathbf{p}) = \frac{\mu_p}{p} \{ P_{12}(p) + f\mu_p^2 P_{22}(p) \}. \quad (B20)$$

## 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., [33,34,45,49]). 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 a triangle, i.e., the length of two wave vectors $\mathbf{k}_1$ and $\mathbf{k}_2$, and the angle between them, $\theta_{12} \equiv \cos^{-1}(\hat{\mathbf{k}}_1 \cdot \hat{\mathbf{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. [34], the dependence of the orientation of the triangle is conveniently described by decomposing into spherical harmonics,

$$B_s(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = \sum_{\ell=0}^{\infty} \sum_{m=-\ell}^{\ell} B_s^{(\ell,m)}(\mathbf{k}_1, \mathbf{k}_2, \theta_{12}) Y_{\ell m}(\omega, \phi). \quad (C1)$$

Similar to Ref. [34], we focus on the $m = 0$ multipoles, which correspond to averaging over $\phi$. Then, the above equation leads to

$$\int_0^{2\pi} \frac{d\phi}{2\pi} B_s(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3) = \sum_{\ell=0}^{\infty} B_s^{(\ell)}(k_1, k_2, \theta_{12}) P_\ell(\mu), \quad (C2)$$

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

$$B_s^{(\ell)}(k_1, k_2, \theta_{12}) = \frac{2\ell + 1}{2} \int_{-1}^{1} d\mu P_\ell(\mu) \int_0^{2\pi} \frac{d\phi}{2\pi} B_s(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3). \quad (C3)$$

In the above, the bispectrum multipole, $B_s^{(\ell)}$, is the quantity of our interest, but at this moment, Eq. (C3) is ambiguous because we do not yet specify what are $\omega$ and $\phi$. To describe the orientation of the triangle, a simple way is to choose $\mathbf{k}_1$ specifically, and to define the orientation angle as $\mu = \cos \omega = \hat{\mathbf{k}}_1 \cdot \hat{\mathbf{z}}$ [34]. Then, we set $\phi$ to the azimuthal angle around $\mathbf{k}_1$. In this case, however, the resultant bispectrum multipoles $B_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. 2. Given two vectors $\mathbf{k}_1$ and $\mathbf{k}_2$, we define the orientation of the triangle by the angle between the line-of-sight direction and the vector normal to the triangle, i.e., $\hat{\mathbf{k}}_1 \times \hat{\mathbf{k}}_2$. That is,

$$\mu = \cos \omega = \frac{(\hat{\mathbf{k}}_1 \times \hat{\mathbf{k}}_2) \cdot \hat{\mathbf{z}}}{\sin \theta_{12}}. \quad (C4)$$

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

$$\cos \phi = \frac{\{\hat{\mathbf{z}} \times (\hat{\mathbf{k}}_1 \times \hat{\mathbf{k}}_2)\} \cdot \hat{\mathbf{k}}_1}{\sin \omega}. \quad (C5)$$

One may suspect that the above definition still breaks the symmetry. Indeed, the vectors $\mathbf{k}_1$ and $\mathbf{k}_2$ can be arbitrarily chosen among three vectors, and we can even exchange $\hat{\mathbf{k}}_1 \leftrightarrow \hat{\mathbf{k}}_2$. Then the orientation angle is changed to $\omega \to \pi - \omega$. Nevertheless, the bispectrum multipole $B_s^{(\ell)}$ is invariant. This is because only the $\ell = 0$ modes of the bispectrum become nonvanishing, and we have $P_\ell(\cos \omega) = P_\ell(\cos(\pi - \omega))$ for $\ell = 0$.

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

(i) In harmonic space, we first pick up the three density fields, $\delta(\mathbf{k}_1)$, $\delta(\mathbf{k}_2)$, and $\delta(\mathbf{k}_3)$, with the vectors satisfying $\mathbf{k}_1 + \mathbf{k}_2 + \mathbf{k}_3 = \mathbf{0}$. Here, labels of $\mathbf{k}_1, \mathbf{k}_2$, and $\mathbf{k}_3$ are arbitrary.
(ii) Choosing the two vectors among the three, we calculate the directional cosine $\mu$ according to Eq. (C4).

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

(iv) 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_s^{(\ell)}$, which is characterized by the length of two vectors, and the angle between them [see Eq. (C3)].

[1] M. Davis and P. J. E. Peebles, Astrophys. J. **267**, 465 (1983).

[2] N. Kaiser, Mon. Not. R. Astron. Soc. **227**, 1 (1987).

[3] A. J. S. Hamilton, Astrophys. J. **385**, L5 (1992).

[4] R. Scoccimarro, Phys. Rev. D **70**, 083007 (2004).

[5] J. C. Jackson, Mon. Not. R. Astron. Soc. **156**, 1P (1972).

[6] W. J. Percival and M. White, Mon. Not. R. Astron. Soc. **393**, 297 (2009).

[7] Y.-S. Song and W. J. Percival, J. Cosmol. Astropart. Phys. **10** (2009) 004.

[8] K. Yamamoto, T. Sato, and G. Huetsi, Prog. Theor. Phys. **120**, 609 (2008).

[9] D. J. Eisenstein and W. Hu, Astrophys. J. **496**, 605 (1998).

[10] A. Meiksin, M. J. White, and J. A. Peacock, Mon. Not. R. Astron. Soc. **304**, 851 (1999).

[11] D. J. Eisenstein et al. (SDSS Collaboration), Astrophys. J. **633**, 560 (2005).

[12] A. Taruya, T. Nishimichi, and S. Saito, Phys. Rev. D **82**, 063522 (2010).

[13] B. A. Reid et al., Mon. Not. R. Astron. Soc. **426**, 2719 (2012).

[14] A. Oka, S. Saito, T. Nishimichi, A. Taruya, and K. Yamamoto, Mon. Not. R. Astron. Soc. **439**, 2515 (2014).

[15] F. Beutler et al. (BOSS Collaboration), Mon. Not. R. Astron. Soc. **443**, 1065 (2014).

[16] E. V. Linder, Astropart. Phys. **29**, 336 (2008).

[17] H.-J. Seo and D. J. Eisenstein, Astrophys. J. **598**, 720 (2003).

[18] C. Blake and K. Glazebrook, Astrophys. J. **594**, 665 (2003).

[19] K. Glazebrook and C. Blake, Astrophys. J. **631**, 1 (2005).

[20] Y.-S. Song, A. Taruya, E. Linder, K. Koyama, C. G. Sabiu, G.-B. Zhao, F. Bernardeau, T. Nishimichi, and T. Okumura, Phys. Rev. D **92**, 043522 (2015).

[21] L. Guzzo et al., Nature (London) **451**, 541 (2008).

[22] T. Okumura et al., Publ. Astron. Soc. Jpn. **68**, 38 (2016).

[23] S. de la Torre et al., Astron. Astrophys. **557**, A54 (2013).

[24] F. Beutler et al. (BOSS Collaboration), Mon. Not. R. Astron. Soc. **466**, 2242 (2017).

[25] C. Blake et al., Mon. Not. R. Astron. Soc. **415**, 2876 (2011).

[26] P. A. R. Ade et al. (Planck Collaboration), Astron. Astrophys. **594**, A13 (2016).

[27] E. Sefusatti, M. Crocce, S. Pueblas, and R. Scoccimarro, Phys. Rev. D **74**, 023522 (2006).

[28] I. Kayo, M. Takada, and B. Jain, Mon. Not. R. Astron. Soc. **429**, 344 (2013).

[29] M. Sato and T. Nishimichi, Phys. Rev. D **87**, 123538 (2013).

[30] B. Greig, E. Komatsu, and J. S. B. Wyithe, Mon. Not. R. Astron. Soc. **431**, 1777 (2013).

[31] J. Byun, A. Eggemeier, D. Regan, D. Seery, and R. E. Smith, arXiv:1705.04392.

[32] Y.-S. Song, A. Taruya, and A. Oka, J. Cosmol. Astropart. Phys. **08** (2015) 007.

[33] P. Gagrani and L. Samushia, Mon. Not. R. Astron. Soc. **467**, 928 (2017).

[34] R. Scoccimarro, H. M. P. Couchman, and J. A. Frieman, Astrophys. J. **517**, 531 (1999).

[35] B. A. Reid and M. White, Mon. Not. R. Astron. Soc. **417**, 1913 (2011).

[36] Z. Vlah, U. Seljak, P. McDonald, T. Okumura, and T. Baldauf, J. Cosmol. Astropart. Phys. **11** (2012) 009.

[37] T. Matsubara, Phys. Rev. D **77**, 063530 (2008).

[38] T. Matsubara, Phys. Rev. D **78**, 083519 (2008); **78**, 109901 (E) (2008).

[39] J. Carlson, B. Reid, and M. White, Mon. Not. R. Astron. Soc. **429**, 1674 (2013).

[40] T. Matsubara, Phys. Rev. D **90**, 043537 (2014).

[41] A. Taruya, T. Nishimichi, and F. Bernardeau, Phys. Rev. D **87**, 083509 (2013).

[42] H. Gil-Marín, C. Wagner, F. Fragkoudi, R. Jimenez, and L. Verde, J. Cosmol. Astropart. Phys. **02** (2012) 047.

[43] R. Scoccimarro and H. M. P. Couchman, Mon. Not. R. Astron. Soc. **325**, 1312 (2001).

[44] R. E. Smith, R. K. Sheth, and R. Scoccimarro, Phys. Rev. D **78**, 023523 (2008).

[45] K. Yamamoto, Y. Nan, and C. Hikage, Phys. Rev. D **95**, 043528 (2017).

[46] C. Rampf and Y. Y. Y. Wong, J. Cosmol. Astropart. Phys. **06** (2012) 018.

[47] E. Komatsu et al., Astrophys. J. Suppl. Ser. **192**, 18 (2011).

[48] F. Bernardeau, S. Colombi, E. Gaztañaga, and R. Scoccimarro, Phys. Rep. **367**, 1 (2002).

[49] R. Scoccimarro, Phys. Rev. D **92**, 083532 (2015).

[50] L. Blot, P. S. Corasaniti, J.-M. Alimi, V. Reverdy, and Y. Rasera, Mon. Not. R. Astron. Soc. **446**, 1756 (2015).

[51] S. Prunet, C. Pichon, D. Aubert, D. Pogosyan, R. Teyssier, and S. Gottloeber, Astrophys. J. Suppl. Ser. **178**, 179 (2008).

[52] R. Teyssier, Astron. Astrophys. **385**, 337 (2002).

[53] Y. Rasera, P.-S. Corasaniti, J.-M. Alimi, V. Bouillot, V. Reverdy, and I. Balmès, Mon. Not. R. Astron. Soc. **440**, 1420 (2014).
[54] R. Scoccimarro, S. Colombi, J. N. Fry, J. A. Frieman, E. Hivon, and A. Melott, Astrophys. J. **496**, 586 (1998).

[55] R. Takahashi et al., Mon. Not. R. Astron. Soc. **389**, 1675 (2008).

[56] K. Heitmann, M. White, C. Wagner, S. Habib, and D. Higdon, Astrophys. J. **715**, 104 (2010).

[57] K. Heitmann, D. Higdon, M. White, S. Habib, B. J. Williams, and C. Wagner, Astrophys. J. **705**, 156 (2009).

[58] E. Lawrence, K. Heitmann, M. White, D. Higdon, C. Wagner, S. Habib, and B. Williams, Astrophys. J. **713**, 1322 (2010).

[59] K. Heitmann, E. Lawrence, J. Kwan, S. Habib, and D. Higdon, Astrophys. J. **780**, 111 (2014).

[60] A. Schneider, R. Teyssier, D. Potter, J. Stadel, J. Onions, D. S. Reed, R. E. Smith, V. Springel, F. R. Pearce, and R. Scoccimarro, J. Cosmol. Astropart. Phys. 04 (2016) 047.

[61] T. Nishimichi et al., Publ. Astron. Soc. Jpn. **61**, 321 (2009).

[62] D. Baumann, A. Nicolis, L. Senatore, and M. Zaldarriaga, J. Cosmol. Astropart. Phys. 07 (2012) 051.

[63] J. J. M. Carrasco, M. P. Hertzberg, and L. Senatore, J. High Energy Phys. 09 (2012) 82.

[64] M. P. Hertzberg, Phys. Rev. D **89**, 043521 (2014).

[65] T. Baldauf, L. Mercolli, and M. Zaldarriaga, Phys. Rev. D **92**, 123007 (2015).

[66] M. Sato and T. Nishimichi, Phys. Rev. D **87**, 123538 (2013).