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

# Fitting the Nonlinear Matter Bispectrum by the Halofit Approach

Ryuichi Takahashi<fnref n="1" />, Takahiro Nishimichi<fnref n="2" />,<fnref n="3" />, Toshiya Namikawa<fnref n="4" />, Atsushi Taruya<fnref n="2" />,<fnref n="3" />, Issha Kayo<fnref n="5" />, Ken Osato<fnref n="6" />, and Yosuke Kobayashi<fnref n="3" />

Received 2019 November 19; revised 2020 May 2; accepted 2020 May 4; published 2020 June 3

## Abstract

We provide a new fitting formula of the matter bispectrum in the nonlinear regime calibrated by high-resolution cosmological N-body simulations of 41 cold dark matter (wCDM, w = constant) models around the Planck 2015 best-fit parameters. As the parameterization in our fitting function is similar to that in Halofit, our fitting is named BiHalofit. The simulation volume is sufficiently large (>10 Gpc<fnref n="3" />) to cover almost all measurable triangle bispectrum configurations in the universe. The function is also calibrated using one-loop perturbation theory at large scales ($k < 0.3 h \text{ Mpc}^{-1}$). Our formula reproduced the matter bispectrum to within 10% (15%) accuracy in the Planck 2015 model at wavenumber $k < 3 (10) h \text{ Mpc}^{-1}$ and redshifts $z = 0-3$. The other 40 wCDM models obtained poorer fits, with accuracy approximating 20% at $k < 3 h \text{ Mpc}^{-1}$ and $z = 0-1.5$ (the deviation includes the 10%-level sample variance of the simulations). We also provide a fitting formula that corrects the baryonic effects such as radiative cooling and active galactic nucleus feedback, using the latest hydrodynamical simulation IllustrisTNG. We demonstrate that our new formula more accurately predicts the weak-lensing bispectrum than the existing fitting formulae. This formula will assist current and future weak-lensing surveys and cosmic microwave background lensing experiments. Numerical codes of the formula are available, written in Python (https://toshiyan.github.io/clpdoc/html/basic/basic.html#module-basic.bispec), C, and Fortran (http://cosmo.phys.hirosaki-u.ac.jp/takahasi/codes_e.htm).

Unified Astronomy Thesaurus concepts: N-body simulations (1083); Weak gravitational lensing (1797); Cosmic microwave background radiation (322); Large-scale structure of the universe (902); Cosmological parameters from large-scale structure (340)

## 1. Introduction

Observations of the cosmic microwave background (CMB) have revealed that the primordial density fluctuations are well described by a Gaussian field (Planck Collaboration 2019). The statistical property of a Gaussian field is fully described by the two-point correlation function or its Fourier transform, the power spectrum (PS). However, at late times, the density fluctuations become non-Gaussian through small-scale gravitational evolution. To fully characterize the statistical property of the non-Gaussian field and to access its cosmological information beyond the two-point (2pt) statistics, higher-order statistics are Required. The leading correction term of the commonly used PS is the bispectrum (BS), the Fourier transform of the three-point (3pt) correlation function.

The 3pt correlation function was first measured in the angular clustering of galaxies (Peebles & Groth 1975; Groth & Peebles 1977). Several groups later measured the 3pt statistics in redshift space using spectroscopic survey data (e.g., Jing & Börner 1998; Scoccimarro et al. 2001; Kayo et al. 2004). The use of 3pt statistics breaks a degeneracy between the galaxy bias and cosmological parameters (e.g., Fry & Gaztanaga 1993; Matarrese et al. 1997; Nishimichi et al. 2007). From the recent analyses for the Baryon Oscillation Spectroscopic Survey data,<fnref n="8" /> as a part of the Sloan Digital Sky Survey,<fnref n="9" /> the baryon

acoustic oscillation features were detected in the BS and 3pt correlation function (Gil-Marín et al. 2016; Slepian et al. 2017). The 3pt statistics contain valuable information complementary to the 2pt statistics and helped to tighten the constraints on the angular diameter distance to galaxies and the redshift space distortion.

Among various observables of large-scale structure, weak lensing can map a projected density field through the coherent distortion of background galaxies (e.g., Bartelmann & Schneider 2001). Current active weak-lensing surveys include the Subaru Hyper Suprime-Cam (HSC),<fnref n="10" /> the Dark Energy Survey (DES),<fnref n="11" /> and the Kilo-Degree Survey (KiDS).<fnref n="12" /> These surveys have placed strong constraints on the cosmological parameters such as the matter density $\Omega_m$ and the amplitude of density fluctuations $\sigma_8$ from the cosmic-shear two-point function (e.g., Abbott et al. 2018; van Uitert et al. 2018; Hamana et al. 2020; Hikage et al. 2019). In the 2020s, ground- and space-based missions such as the Large Synoptic Survey Telescope (LSST),<fnref n="13" /> Wide Field Infrared Survey Telescope (WFIRST),<fnref n="14" /> and Euclid<fnref n="15" /> will commence operations.

<footnote n="1">¹Faculty of Science and Technology, Hirosaki University, 3 Bunkyo-cho, Hirosaki, Aomori 036-8588, Japan</footnote>

<footnote n="2">²Center for Gravitational Physics, Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan</footnote>

<footnote n="3">³Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study (UTIAS), The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa-shi, Chiba, 277-8583, Japan</footnote>

<footnote n="4">⁴Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Cambridge CB3 0WA, UK</footnote>

<footnote n="5">⁵Department of Liberal Arts, Tokyo University of Technology, Ota-ku, Tokyo 144-8650, Japan</footnote>

<footnote n="6">⁶Institut d'Astrophysique de Paris, Sorbonne Université, CNRS, UMR 7095, 98bis boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="7">⁷National Astronomical Observatory of Japan, Mitaka, Tokyo, 181-8588, Japan</footnote>

<footnote n="8">⁸ https://www.sdss.org/surveys/boss/</footnote>

<footnote n="9">⁹ https://www.sdss.org/</footnote>

<footnote n="10">¹⁰ https://hsc.mtk.nao.ac.jp/ssp/</footnote>

<footnote n="11">¹¹ https://www.darkenergysurvey.org</footnote>

<footnote n="12">¹² http://kids.strw.leidenuniv.nl</footnote>

<footnote n="13">¹³ https://www.lsst.org/</footnote>

<footnote n="14">¹⁴ https://wfirst.gsfc.nasa.gov/</footnote>

<footnote n="15">¹⁵ https://www.euclid-ec.org/</footnote>
The weak-lensing BS contains additional information that complements the PS. Because it arises from the non-Gaussian properties, the BS is more sensitive to smaller-scale and lower-redshift structures than the PS. A joint analysis of both the PS and BS breaks parameter degeneracy and provides tighter constraints (e.g., Takada & Jain 2004; Kilbinger & Schneider 2005; Sefusatti et al. 2006; Munshi et al. 2011; Kayo & Takada 2013; Byun et al. 2017; Gatti et al. 2019). The BS can be comparable to or more powerful than the PS (Bergé et al. 2010; Sato & Nishimichi 2013; Coulton et al. 2019). Several groups have derived useful constraints from the three-point cosmic-shear statistics of real data (Bernardeau et al. 2002b; Jarvis et al. 2004; Semboloni et al. 2011b; Van Waerbeke et al. 2013; Fu et al. 2014; Simon et al. 2015). The higher-order moments of weak-lensing convergence also contain the non-Gaussian information (e.g., Petri et al. 2015). The DES will set the observational constraint from a joint analysis of the second- and third-order moments (Chang et al. 2018; Gatti et al. 2019).

CMB lensing is another promising cosmological probe of the density fluctuations at higher redshifts ($z \simeq 1-3$) than cosmic shear (e.g., Lewis & Challinor 2006). Recent CMB experiments have measured the lensing signals from temperature and polarization fluctuations (Planck Collaboration 2018a). The CMB lensing-potential PS provides rich cosmological information that complements the information in galaxy weak lensing (e.g., Planck Collaboration 2018a). The BS and higher-order spectra representing the non-Gaussian density fluctuations would be important in future CMB lensing observations. The non-Gaussianity slightly affects the lensing PS (Pratten & Lewis 2016), as well as the CMB temperature and polarization PSs (Lewis & Pratten 2016; Marozzi et al. 2018). It also contaminates the CMB lensing reconstruction (Böhm et al. 2016; Beck et al. 2018; Fabbian et al. 2019). On the other hand, the lensing BS can be measured as a useful signal in future CMB experiments (Namikawa 2016).

Against this background, an accurate model of nonlinear BS is highly demanded. A nonlinear model of the PS with a few percent accuracy up to $k = 10 h \text{ Mpc}^{-1}$ is also required to meet the statistical accuracy requirements of forthcoming weak-lensing surveys<fnref n="16" /> (Huterer & Takada 2005; Hearin et al. 2012). Scoccimarro & Couchman (2001) calibrated a fitting formula of BS in N-body simulations, which was later improved by Gil-Marín et al. (2012). However, the squeezed BS computed by these formulae is double (in the worst cases) that obtained in the latest numerical simulations (Fu et al. 2014; Coulton et al. 2019; Munshi et al. 2020; Namikawa et al. 2019). In this paper, we construct an improved fitting formula of the nonlinear matter BS calibrated in high-resolution cosmological N-body simulations of 41 wCDM models (where wCDM refers to cold dark matter and dark energy with a constant equation of state w). Mainly, we aim to construct the formula for the Planck 2015 ΛCDM model up to $k = 30 h \text{ Mpc}^{-1}$ in the redshift range $z = 0-10$, hoping that the formula has little dependence on cosmology. The other 40 wCDM models supplement the calibration at relatively low redshifts ($z = 0-1.5$). This allows us to explicitly examine the cosmological model dependence, which was not thoroughly done previously. We also include the calibration from one-loop perturbation theory at $k < 0.3 h \text{ Mpc}^{-1}$ in the $z = 0-10$ range, because the simulations have large sample variance at the largest scales.

To ensure an accurate calibration, we bin the simulation data and theoretical prediction into wavenumbers ($k_1, k_2, k_3$). We also consider the baryonic effects in a public hydrodynamic simulation package called the IllustrisTNG suite (Nelson et al. 2019).

The outline of this paper is as follows. Section 2 discusses the basics of matter BS and gives the previous and our own fitting formulae. Section 3 details our simulations. Section 4 describes the fitting procedures and presents the resulting fitting function (Figures 5–9). Section 5 discusses the baryonic effects on the BS using the IllustrisTNG data set. Section 6 compares the fitting formula predictions of the weak-lensing convergence BS with those of light-cone simulations. Section 7 discusses the systematics of cosmic-shear BS and CMB lensing BS. The main paper concludes with a summary in Section 8. Appendix A briefly discusses the halo model, and Appendixes B and C give the fitting formula and the baryonic correction, respectively.

## 2. Theory

### 2.1. Basics

The cosmological density contrast is usually described by its Fourier transform $\tilde{\delta}(k)$. The matter PS and BS are respectively defined as

$$ 
\begin{aligned}
P(k_1) \delta_D(k_1 + k_2) &= \langle \tilde{\delta}(k_1) \tilde{\delta}(k_2) \rangle, \\
B(k_1, k_2, k_3) \delta_D(k_1 + k_2 + k_3) &= \langle \tilde{\delta}(k_1) \tilde{\delta}(k_2) \tilde{\delta}(k_3) \rangle,
\end{aligned} 
\quad (1) $$

where $\delta_D$ is the Dirac delta function. Throughout this paper, we omit the redshift dependence in the arguments of functions because our discussion considers arbitrary redshifts.

At the tree level (i.e., the leading order in perturbation theory), the matter BS is given by the product of the linear matter PS, $P_L(k)$, as follows (e.g., Bernardeau et al. 2002a):

$$ B_{\text{tree}}(k_1, k_2, k_3) = 2F_2(k_1, k_2)P_L(k_1)P_L(k_2) + 2 \text{perm.} \quad (2) $$

Here the last term describes two permutations $(k_1, k_2) \rightarrow (k_2, k_3)$ and $(k_3, k_1)$, which are applied to the wavevectors in the first term. The $F_2$ kernel is

$$ F_2(k_1, k_2) = \frac{5}{7} + \frac{1}{2} \left( \frac{k_1}{k_2} + \frac{k_2}{k_1} \right) \mu_{12} + \frac{2}{7} \mu_{12}^2, \quad (3) $$

where $\mu_{12}$ is the cosine of the angle between $k_1$ and $k_2$, i.e., $\mu_{12} = k_1 \cdot k_2/(k_1 k_2)$.

To explore the nonlinear regime beyond the tree level, one usually relies on higher-order perturbation theories (e.g., Scoccimarro et al. 1998; Rampf & Wong 2012; Angulo et al. 2015; Hashimoto et al. 2017; Bose & Taruya 2018; Lazanu & Liguori 2018). However, these are reliable only up to the mildly nonlinear regime ($k \lesssim 0.2 h \text{ Mpc}^{-1}$). Another strategy adopts the analytical halo model (e.g., Cooray & Sheth 2002), which assumes that all matter is confined to halos. This model is valid over a wide range of scales and redshifts, but its current accuracy is approximately 30% (e.g., Lazanu et al. 2016; Bose et al. 2020). The last one is a fitting function calibrated in N-body simulations over various scales, epochs, and cosmological models.

### 2.2. Previous Fitting Formulae

Scoccimarro & Couchman (2001) (SC01) provided a fitting formula for the nonlinear BS. Their function is similar to the

<footnote n="16">¹⁶ To our knowledge, the required accuracy of the BS model for current and forthcoming surveys has not been estimated.</footnote>
tree-level formula (Equation (2)) but replaces the linear PS with a nonlinear model and modifies the $F_2$ kernel to enhance the BS amplitude at small scales. In the low-$k$ limit, their formula is consistent with the tree level. In the high-$k$ limit, the BS is proportional to $P(k_1)P(k_2) + P(k_2)P(k_3) + P(k_3)P(k_1)$ according to the hyperextended perturbation theory (Scoccolmarro & Frieman 1999). Their modified $F_2$ kernel contains six free parameters, which are fitted by their N-body results in four CDM models with $k < 3$ h Mpc$^{-1}$ and $0 \le z \le 1$. Later, Gil-Marín et al. (2012, hereafter GM12) increased the number of free parameters in $F_2$ to nine and recalibrated them from their N-body simulations in a single ΛCDM model over a relatively narrow range of wavenumbers ($k < 0.4$ h Mpc$^{-1}$) with $0 \le z \le 1.5$.

However, these formulae have several shortcomings. First, they are based on a nonlinear PS model such as Halofit (Smith et al. 2003; Takahashi et al. 2012), HMcode (Mead et al. 2015), or COSMIC EMULATOR (Lawrence et al. 2017). This PS model needs to be prepared by the user along with the BS formula. The discrepancies among these PS models are small but negligible, typically a few percent (e.g., Schneider et al. 2016); accordingly, they degrade the BS accuracy. Second, as indicated by Namikawa et al. (2019), these models overestimate the squeezed BS (i.e., the configuration of $k_1 \simeq k_2 \gg k_3$). Third, their fitting range of $k$ and $z$ is narrow. The current weak-lensing surveys measure the correlation function down to arcmin scales, requiring calibration up to $k = 10$ h Mpc$^{-1}$. In addition, as the CMB lensing probes the high-redshift structures ($z \simeq 1-3$), the calibration must extend at least to $z = 3$. Finally, these models do not consider the baryonic effects, which are important at $k \gtrsim 1$ h Mpc$^{-1}$.

## 2.3. Our Fitting Formula

Our fitting formula is based on the halo model and is similar to Halofit for the nonlinear PS (Smith et al. 2003). The halo model is popular for evaluating the multipoint statistics of nonlinear density fields (the halo model BS is detailed in Appendix A). Assuming that all particles are contained in halos, it decomposes the BS into three terms: one-, two-, and three-halo terms (hereafter denoted as 1h, 2h, and 3h, respectively). The 1h term describes the correlation in an individual halo, and the 2h (3h) term accounts for the correlation among two (three) different halos. The 1h and 3h terms dominate at small and large scales, respectively. Because the 2h term is subdominant in most of the triangle configurations (except in the squeezed case; see, e.g., Valageas & Nishimichi 2011; Valageas et al. 2012), it is dropped in our formulation and is absorbed by enhancing the 3h term at intermediate scales.

The fitting function consists of two terms,

$$B(k_1, k_2, k_3) = B_{1h}(k_1, k_2, k_3) + B_{3h}(k_1, k_2, k_3), \quad (4)$$

and approaches the tree-level formula in the low-$k$ limit. The function contains 52 free parameters to be fitted by our N-body data. The fitting function is explicitly given in Appendix B.

One may consider that 52 free parameters are many. However, given the huge number of triangle configurations ($\sim 5 \times 10^5$) for all wavenumbers, redshifts, and cosmological models in our calibration, the number of parameters is rather small. In addition, recalling that the revised Halofit (Takahashi et al. 2012) for the PS already contains 34 free parameters, SC01 and GM12 contain 40 and 43 parameters

in total, respectively (using the nonlinear PS from Halofit). Therefore, our parameters are only slightly more than in these previous models.

## 3. Numerical Results

The fitting formula was calibrated in cosmological dark matter N-body simulations. We used the N-body data set prepared by the DARK EMULATOR project (Nishimichi et al. 2019, hereafter N19). The N19 project has prepared 101 flat cosmological models (a fiducial ΛCDM and additional 100 wCDM models) in the range $z = 0-1.48$. The project aims to emulate several halo observables such as the halo-matter correlation function, the halo mass function, and the halo bias for ongoing weak-lensing surveys. The emulator will be publicly available soon.

### 3.1. Cosmological Models

We used the N19 simulations of 41 flat cosmological models.<fnref n="17" /> The fiducial ΛCDM model is consistent with the Planck 2015 best fit (Planck Collaboration 2016), with matter density $\Omega_m = 1 - \Omega_\Lambda = 0.3156$, baryon density $\Omega_b = 0.0492$, Hubble parameter $h = 0.6727$, spectral index $n_s = 0.9645$, and amplitude of matter density fluctuation on the scale of $8 h^{-1}$ Mpc $\sigma_8 = 0.831$.

The other wCDM models have six cosmological parameters: $\Omega_b h^2$, $\Omega_{cdm} h^2$, $\Omega_w$, $A_s$, $n_s$, and $w$. Here the dark energy equation of state $w$ is assumed to be constant, and $A_s$ is the amplitude of the primordial PS. These parameters are distributed around the fiducial model in the ranges $\pm 5\%$ for $\Omega_b h^2$ and $n_s$, $\pm 10\%$ for $\Omega_{cdm} h^2$, and $\pm 20\%$ for $\Omega_w$, $\ln A_s$, and $w$. In the N19 project, the cosmological parameters were sampled using a Latin Hypercube Design (e.g., Heitmann et al. 2009). The models were placed into five subsets, each containing 20 models. Figure 1 shows the distributions of $w$ and $\sigma_8$ versus $\Omega_m$ in the 41 models (the fiducial ΛCDM models and two subsets of N19) considered in the present study. The parameter range is wide enough for current and future weak-lensing surveys. In fact, the current constraint from the HSC (DES) cosmic-shear 2pt statistics alone is $S_8 = 0.795^{+0.043}_{-0.047}$ ($0.789^{+0.036}_{-0.038}$) in the flat wCDM model (Troxel et al. 2018; Hamana et al. 2020).

Although these simulations are dark-matter-only simulations, their initial condition accounts for the free-streaming damping by massive neutrinos. To compute the linear matter transfer function at the initial redshift of the simulations, N19 first computed the one at $z = 0$ with massive neutrinos and then multiplied it by the ratio of the linear growth factor between $z = 0$ and the target redshift, in which the scale-dependent growth due to neutrinos was neglected. The same procedure was done in this work. The neutrino density in all models was fixed at $\Omega_\nu h^2 = 6.4 \times 10^{-4}$, corresponding to a total mass 0.06 eV. This $\Omega_\nu$ is included in $\Omega_m$.

### 3.2. N-body Simulations

Our simulation settings are summarized in Table 1. To cover a wide range of length scales, we set four box sizes ($L = 4, 2, 1$, and $0.2 h^{-1}$ Gpc, where $L$ is the side length of the cubic box in the comoving scale). Note that the large simulation volume can

<footnote n="17">¹⁷ Unfortunately, the particle position data were lost for the rest (60) of the models owing to hard-disk trouble.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Distribution of cosmological parameters in the 41 models. The central red diamond is the Planck 2015 best-fit ΛCDM, while the others are the 40 wCDM models. The latter are divided into three groups of $S_8 \equiv \sigma_8(\Omega_m/0.3)^{0.5}$; $S_8 > 1.0$ (magenta circles), $S_8 < 0.6$ (cyan circles), and $0.6 < S_8 < 1.0$ (gray circles). These three groups are separated by the boundaries (dashed lines) at $S_8 = 0.6$ and 1.0.</figcaption></figure>

include almost all measurable triangle configurations of BS in the real universe. The large-volume simulations ($L = 4 h^{-1}$ Gpc) reduce the sample variance in the measured BS at small $k$, whereas the small-volume simulations ($L = 0.2 h^{-1}$ Gpc) reveal the asymptotic behavior at high $z$. Here the simulations with $L = 1$ and $2 h^{-1}$ Gpc at $z = 0-1.48$ are taken from N19, while the others are newly prepared in this work. The largest- and smallest-box simulations supplement the dynamic range covered by N19. The number of particles was set to $2048^3$ except for $L = 4 h^{-1}$ Gpc (where it was $4096^3$). The resulting particle Nyquist wavenumber is $k_{Ny} = \pi/l_p$, where $l_p (=n_p^{-1/3})$ is the mean interparticle separation at particle number density $n_p$. The $k_{Ny}$ values are listed in Table 1. The fiducial ΛCDM model has dozens of independent realizations, whereas each wCDM model has a single realization.

The initial matter PS was prepared by the public Boltzmann code CAMB (Lewis et al. 2000). The initial particle distribution was determined by the second-order Lagrangian perturbation theory (2LPT; Crocce et al. 2006; Nishimichi et al. 2009)<fnref n="18" /> at redshifts of $z_{in} = 31, 29, 59$, and 99 for $L = 4, 2, 1$, and $0.2 h^{-1}$ Gpc, respectively. The initial redshifts in the 40 wCDM models were changed because the initial amplitudes differed among the models. The initial redshift was determined by requiring the rms displacement to be 25% of the mean interparticle separation to achieve the optimal balance between the artificial force due to the grid pre-initial configuration and the transient due to the truncation of the LPT at the second order. The nonlinear gravitational evolution was followed using a tree-PM (particle mesh) code GADGET2 (Springel et al. 2001; Springel 2005). The number of

PM grid cells was $4096^3$ ($8192^3$ for $L = 4 h^{-1}$ Gpc). The gravitational softening length was set to 5% of the mean interparticle separation. The GADGET2 parameters (such as time step and force calculation parameters) were fine-tuned to determine the matter PS with percentage-level accuracy in N19 (see Section 3.4 of their paper). The particle snapshots were dumped at 14 redshifts ranging from $z = 0$ to 10 (see Table 1 for the exact redshift outputs).

To measure the density contrast, we assigned the N-body particles to the $1024^3$ regular grid cells in the box using the cloud-in-cell (CIC) interpolation with the interlacing scheme (e.g., Jing 2005; Sefusatti et al. 2016). The Fourier transform $\delta(k)$ of the density field was then obtained by fast Fourier transform (FFT).<fnref n="19" /> To explore smaller scales, we also employed the folding method (Jenkins et al. 1998), which folds the particle positions $x$ into a smaller box of side length $L/n$ by replacing $x$ with $x%(L/n)$ (where $a%b$ obtains a reminder of $a/b$). This procedure effectively increases the resolution by $n$ times. Here we set $n = 4$ and 10 at $L = 4$ and $1 h^{-1}$ Gpc, respectively. The minimum and maximum wavenumbers in the $1024^3$ cells were $k_{min}^L = 2\pi/L = 6.3 \times 10^{-3} h$ Mpc$^{-1}$ [$L/(h^{-1}$ Gpc)]$^{-1}$ and $k_{max}^L = 512 k_{min}^L = 3.2 h$ Mpc$^{-1}$ [$L/(h^{-1}$ Gpc)]$^{-1}$, respectively. The folding scheme simply enlarged both $k_{min}^L$ and $k_{max}^L$ by 4 or 10 times. The resultant $k_{max}^L$ values are given in Table 1.

### 3.3. PS Measurement for Accuracy Check

The numerical accuracy was checked by comparing the simulated matter PS with the results of a previous fitting formula. The PS estimator is given by

$$ \hat{P}(k) = \frac{1}{N_{\text{mode}}^{\text{PS}}} \sum_{|k'| \in k} |\delta(k')|^2, \qquad (5) $$

where the summation is performed over $k - \Delta k/2 < |k'| < k + \Delta k/2$ and $N_{\text{mode}}^{\text{PS}}$ is the number of modes in a fixed bin width ($\Delta \log_{10} k = 0.1$). Figure 2 plots the PS ratio of the simulation to the revised Halofit prediction (Smith et al. 2003; Takahashi et al. 2012). We here plot the average $P(k)$ and its $1\sigma$ error measured from the realizations. The results in the different boxes were nicely consistent. In larger simulation boxes, the measured PS was smaller than the Halofit prediction at large $k$ because of the lack of spatial resolution. The shot noise was not subtracted because it contributed less than 3% on the scales shown in Figure 2. The simulations agreed with the fitting formula within 5% for $k < 1 h$ Mpc$^{-1}$ at $z = 0-10$ and $k < 10 h$ Mpc$^{-1}$ at $z = 0-1.5$.

### 3.4. Bispectrum Measurement

The BS estimator is given by

$$ \hat{B}(k_1, k_2, k_3) = \frac{1}{N_{\text{triangle}}} \sum_{|k'_1| \in k_1} \sum_{|k'_2| \in k_2} \sum_{|k'_3| \in k_3} \times \delta(k'_1) \delta(k'_2) \delta(k'_3) \delta^{\text{K}}_{k'_1+k'_2+k'_3}, \qquad (6) $$

where the summation is performed over all modes in the bin, $|k'_i| \in k_i$ ($i = 1, 2, 3$), $N_{\text{triangle}}$ is the number of triangles, and $\delta^K$ is the Kronecker delta. Throughout this paper, the log-scale

<footnote n="18">¹⁸ The 2LPT reduces the error in the BS estimate caused by transients from the initial condition to below 2% at $z \le 1$ (McCullagh et al. 2016).</footnote>

<footnote n="19">¹⁹ FFTW (Fast Fourier Transform in the West) is available at http://www.fftw.org.</footnote>
<table><caption>**Table 1**
N-body Simulation Parameters</caption><thead><tr><th>Cosmological Model</th><th>Box Size (h<sup>-1</sup> Gpc)</th><th>Number of Particles</th><th>Number of Realizations</th><th>Particle Nyquist Wavenumber (h Mpc<sup>-1</sup>)</th><th>Maximum Wavenumber (h Mpc<sup>-1</sup>)</th><th>Output Redshifts</th></tr></thead><tbody><tr><td rowspan="4">Planck 2015 ΛCDM</td><td>4</td><td>4096<sup>3</sup></td><td>8</td><td>3.22</td><td>2.85</td><td rowspan="4">0, 0.55, 1.03, 1.48 and High Z</td></tr><tr><td>2</td><td>2048<sup>3</sup></td><td>15</td><td>3.22</td><td>1.42</td></tr><tr><td>1</td><td>2048<sup>3</sup></td><td>21</td><td>6.43</td><td>28.5</td><td rowspan="2">Low Z and High Z</td></tr><tr><td>0.2</td><td>2048<sup>3</sup></td><td>10</td><td>32.2</td><td>14.2</td></tr><tr><td rowspan="2">40 wCDM</td><td>2</td><td>2048<sup>3</sup></td><td>1</td><td>3.22</td><td>1.42</td><td rowspan="2">Low Z</td></tr><tr><td>1</td><td>2048<sup>3</sup></td><td>1</td><td>6.43</td><td>28.5</td><td>Low Z</td></tr></tbody></table>

**Note.** In the output redshifts column, Low Z covers 10 low redshifts ($z = 0, 0.15, 0.31, 0.42, 0.55, 0.69, 0.85, 1.03, 1.23$, and $1.48$), and High Z covers four high redshifts ($z = 2, 3, 5$, and $10$). The Low Z simulations with $L = 1$ and $2 h^{-1}$ Gpc are taken from Nishimichi et al. (2019); the others are newly prepared in this work.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** PS ratio of the simulation to the Halofit prediction in the Planck 2015 best-fit ΛCDM model. The results in simulation box sizes of $L = 4, 2, 1$, and $0.2 h^{-1}$ Gpc are indicated by blue, green, black, and gray symbols, respectively. The shot-noise contribution is less than 3%.</figcaption></figure>

bin width is constant ($\Delta \log_{10} k = 0.1$), unless otherwise stated. Equation (6) was calculated by the FFT-based quick estimator (e.g., Scoccimarro 2015). Using the identity $\delta_{k'_1+k'_2+k'_3}^K = N_{cell}^{-1} \sum_x e^{i(k'_1+k'_2+k'_3) \cdot x}$, Equation (6) reduces to

$$ \hat{B}(k_1, k_2, k_3) = \frac{1}{N_{\text{triangle}}} \frac{1}{N_{\text{cell}}} \sum_x \left[ \sum_{|k'_1| \in k_1} \delta(k'_1) e^{ik'_1 \cdot x} \times \sum_{|k'_2| \in k_2} \delta(k'_2) e^{ik'_2 \cdot x} \sum_{|k'_3| \in k_3} \delta(k'_3) e^{ik'_3 \cdot x} \right], \quad (7) $$

where $x$ is a discrete grid coordinate and $N_{cell} = 1024^3$ is the total number of cells. The summation over $k'_i$ is easily performed by FFT. Although Equation (7) can be quickly

computed, the FFTs in all $k'_i$ bins require large memory resources. This demand limits the grid resolution ($N_{cell}$).

The shot noise is measured as

$$ \hat{B}_{\text{sn}}(k_1, k_2, k_3) = \frac{1}{n_{\text{p}}} [\hat{P}(k_1) + \hat{P}(k_2) + \hat{P}(k_3)] - \frac{2}{n_{\text{p}}^2}, $$

where $n_{\text{p}}$ is the particle number density and $\hat{P}(k)$ is the PS estimator including the shot noise.

In the fiducial model, we calculated the average and standard deviation of BS from the realizations (the number of realizations is listed in Table 1). However, the results of the wCDM models have relatively large scatters because each wCDM model has only a single realization. Therefore, the fitting formula was the main calibration formula for the Planck 2015 model, while the other models supplementarily checked its dependence on the cosmological parameters.

# 4. Fitting Procedure and Results

## 4.1. Fitting to the N-body Results

This subsection presents our fitting procedure to the N-body results. The BS fitting function in Equation (4) contains 52 parameters. Arraying these parameters as $p = (p_1, p_2, \dots)$, the best fit is determined by the standard chi-squared analysis:

$$ \chi^2_{\text{sim}}(\mathbf{p}) = \sum_{c=1}^{41} \sum_{z=0}^{10} \sum_{L} \sum_{k_1, k_2, k_3}^{30 h \text{ Mpc}^{-1}} W_c W_z W_k \\ \times \left[ \frac{B^{\text{bin}}(k_1, k_2, k_3; \mathbf{p}) - B_{\text{sim}}(k_1, k_2, k_3)}{\Delta B_{\text{sim}}(k_1, k_2, k_3) + \epsilon(k_1, k_2, k_3)} \right]^2, \quad (8) $$

where the summation is performed for the 41 cosmological models (subscripted by $c$), all redshifts $z = 0-10$ (subscripted by $z$), all simulations in different box sizes $L$, and all triangles ($k_1, k_2, k_3$) up to $30 h \text{ Mpc}^{-1}$. Here $B^{\text{bin}}$ is the binned prediction of the fitting formula (given by Equation (9)), $B_{\text{sim}}$ is the simulation result, and $\Delta B_{\text{sim}}$ is the standard deviation estimated from the N-body realizations. As each of the 40 wCDM models has one realization and these cosmological parameters are basically similar,<fnref n="20" /> their relative standard deviations $\Delta \ln B_{\text{sim}}$ ($\equiv \Delta B_{\text{sim}}/B_{\text{sim}}$) are assumed equal to those of the Planck 2015. Although a small change in $B_{\text{sim}}$ may give a large impact on the resultant fitting formula, a change in $\Delta B_{\text{sim}}$ would not significantly change the result. We also include a “softening” term $\epsilon = 0.02 \times B_{\text{sim}}$ that reduces the influence of data points

<footnote n="20">²⁰ The differences of cosmological parameters among the 41 models are less than 20% (see Section 3.1).</footnote>
with very small $\Delta \ln B_{sim}$ ($\ll \epsilon/B_{sim} = 2\%$).$^{21}$ At large (small) scales where $\Delta \ln B_{sim}$ is larger (smaller) than 2%, the $\Delta B_{sim}$ ($\epsilon$) term dominates the denominator of Equation (8). The $\Delta B_{sim}$ term gives more weight to smaller-scale data (because $\Delta \ln B_{sim}$ is smaller), whereas the $\epsilon$ term gives an equal weight, irrespective of scale. The weight factors ($W$) were introduced to place greater importance on the lower-redshift data ($W_z$) because cosmic shear probes the low-$z$ ($\lesssim 0.5$) structures, and on larger-scale data ($W_k$) because the simulation results are reliable at least up to the particle Nyquist wavenumber (in Table 1) and the unaccounted baryonic effects can influence the small-scale results ($k \gtrsim 1$ h Mpc$^{-1}$). The fiducial cosmological model also received a high weighting ($W_c$).$^{22}$

The analysis included all triangles ($k_1, k_2, k_3$) satisfying the following three conditions:

1. Relative standard deviation below 10% (i.e., $\Delta \ln B_{sim} < 0.1$).

2. Shot-noise contribution below 3%.

3. If the deviation between the larger- and smaller-box simulation results exceeds 3% and the statistical error of $\Delta \ln B_{sim}$ is below 3%, we reject the larger-box result and use the smaller-box result only. In the larger-box simulation, the $B_{sim}$ at high $k$ is reduced by the lack of spatial resolution (see also Figure 2 in the PS case).

Conditions 1 and 2 exclude the data points at very small $k$ and large $k$, respectively. Condition 3 negligibly affects the data selection. Figure 3 plots the triangles ($k_1, k_2, k_3$) satisfying the above three conditions for $L = 1, 2, 4$ h$^{-1}$ Gpc at $z = 0.55$. In the range $0.1 \lesssim k/(h$ Mpc$^{-1}) \lesssim 2$, the simulation results for all box sizes were overlapping and the fit was reliable. As clarified in Figure 3, the simulations covered almost all triangles up to $k = 3$ h Mpc$^{-1}$. Note that most of the triangles were squeezed; the instances of equilateral and flattened cases were minor. Therefore, the fitting to squeezed cases is critically important. In the bottom panel, the discontinuity at $k_{max} \approx 3$ (0.8) h Mpc$^{-1}$ for $L = 1$ (4) h$^{-1}$ Gpc can be explained by the box-size change from $L$ to $L/10$ ($L/4$) when implementing the folding scheme (see also Section 3.2). The bottom panel is devoid of triangles in the lower right part, indicating that the calibration did not include very squeezed cases ($k_{max} \gg k_{min}$). These cases lie outside the maximum $k_{max}/k_{min}$ (512), which is determined by the number of FFT grids (1024$^3$). The folding method does not change this ratio (512). The number of independent triangular bins calibrated in the simulations of each cosmological model was approximately 950 at low $z$ ($z = 0, 0.55, 1$, and $1.48$) and 690 at high $z$ ($z = 2, 3$, and $5$), respectively.

Figure 4 plots the maximum and minimum wavenumbers in the calibration. The minimum $k$ of simulation is larger at higher $z$ because the relative error $\Delta \ln B_{sim}$ is larger. The maximum $k$ decreases at higher $z$ because the shot noise is not negligible at small scales. Note again that all the triangles in this $k$ range are not included in the calibration (e.g., very squeezed cases are missing; see also Figure 3).

<figure><img src="image_4.png" /><figcaption>Figure 3. Triangle configurations included in the calibration from the simulations with $L = 4, 2, 1$ h$^{-1}$ Gpc (blue, green, and black squares) and from the one-loop standard perturbation theory (SPT; orange diamonds). Here $k_{min}$, $k_{mid}$, and $k_{max}$ are the minimum, middle, and maximum side lengths of a triangle. The dashed red lines correspond to particular triangles: $k_{mid} = k_{max}$ is the squeezed, $k_{mid} = k_{max}/2$ is the flattened, and $k_{min} = k_{max}$ is the equilateral. The arrows in the top panel indicate the maximum wavenumbers in the calibration from the simulations (blue, green, and black) and from the SPT (orange). The thick (thin) arrows are with (without) the folding scheme. All the points satisfy conditions 1–3 in Section 4.1.</figcaption></figure>

For a fair comparison, the simulation results and the fitting formula predictions should be binned consistently because the BS is sensitive to the binning, especially at the squeezed limit (Sefusatti et al. 2010; Namikawa et al. 2019). Throughout this paper, the binned fitting was computed as

$$B^{\text{bin}}(k_1, k_2, k_3) = \frac{1}{N_{\text{triangle}}} \int_{|k'_1| \in k_1} d^3 k'_1 \int_{|k'_2| \in k_2} d^3 k'_2 \\ \times \int_{|k'_3| \in k_3} d^3 k'_3 \\ \times B(k'_1, k'_2, k'_3) \delta_D(k'_1 + k'_2 + k'_3), \quad (9)$$

<footnote>$^{21}$ Since the number of realizations are not large enough for estimating the variance accurately, some data points accidentally have very small $\Delta B_{sim}$.</footnote>

<footnote>$^{22}$ Accordingly, the weights were set to $W_z = 8, 3, 1$, and $0.3$ for $z \le 0.1, 0.1 < z \le 1, 1 < z \le 3$, and $z > 3$, respectively; $W_k = 3, 1$ and $0.3$ for $k_{max}/(h$ Mpc$^{-1}) \le 3.2, 3.2 < k_{max}/(h$ Mpc$^{-1}) \le 10$ and $k_{max}/(h$ Mpc$^{-1}) > 10$, respectively; and $W_c = 1$ ($8 \times 10^{-4}$) in the Planck 2015 model (otherwise). These values were chosen to achieve 10% accuracy of the fitting with $k < 3$ h Mpc$^{-1}$ at $z = 0-3$ in the Planck 2015.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Maximum and minimum wavenumbers included in the calibration as a function of 1 + z (in log scale). Symbols denote the simulations with L = 4, 2, 1, and 0.2 h⁻¹ Mpc (blue, green, black, and gray) and the one-loop SPT (orange).</figcaption></figure>

where B(k'_1, k'_2, k'_3) is the unbinned fitting and the number of triangles is

$$N_{\text{triangle}} = \int_{|k'_1| \in k_1} d^3 k'_1 \int_{|k'_2| \in k_2} d^3 k'_2 \int_{|k'_3| \in k_3} d^3 k'_3 \times \delta_{\text{D}}(k'_1 + k'_2 + k'_3). \quad (10)$$

Here kᵢ is the weighted mean wavenumber, defined as $k_i = \int_{|k'_i| \in k_i} d^3 k'_i k'_i \times \left[ \int_{|k'_i| \in k_i} d^3 k'_i \right]^{-1}$. The effect of the binning on BS is shown in Appendix B. Note that although the unbinned triangle (k'_1, k'_2, k'_3) satisfies the triangle condition (i.e., |k'_1 - k'_2| < k'_3 < k'_1 + k'_2), the bin center (k₁, k₂, k₃) may violate this condition.

We now comment on the effect of bin width Δk on the calibration result. A finer bin width reduces the binning uncertainty and improves the calibration, at the cost of increased sample variance (as ΔBsim ∝ Δk⁻³/² under the Gaussian approximation). Therefore, the appropriate Δk is not easily interpreted. Because binning smooths out the fine-k BS features over the bin width, the accuracy of our fitting formula may be degraded if the user adopts a finer Δk than ours. In Section 4.3, we will check the bin width dependence by comparing Δlog₁₀k = 0.1 and 0.05 in the Planck 2015 and confirm the agreement of the two results.

## 4.2. Fitting to Perturbation Theory

As the simulation result is noisy at large scales, the calibration on the linear to quasi-linear scales was also performed by perturbation theory. The same approach was adopted by Smith & Angulo (2019) for modeling the nonlinear matter PS. Here we applied one-loop standard perturbation theory (SPT), which includes the tree level and the next-to-leading-order terms (e.g., Scoccimarro 1997; Scoccimarro et al. 1998). The chi-square was defined analogously to Equation (8):

$$\chi^2_{\text{spt}}(\mathbf{p}) = W_{\text{spt}} \sum_{c=1}^{41} \sum_{z=0}^{10} \sum_{k_1, k_2, k_3}^{0.3 \ h \ \text{Mpc}^{-1}} W_z W_c \times \left[ \frac{B(k_1, k_2, k_3; \mathbf{p}) - B_{\text{spt}}(k_1, k_2, k_3)}{\Delta B(k_1, k_2, k_3) + \epsilon(k_1, k_2, k_3)} \right]^2, \quad (11)$$

where B_spt is the SPT prediction. Note that B and B_spt require no binning in this case. We set ΔB = 0.5 |B_spt - Btree| and ϵ = 0.01 × B_spt. We also set W_spt = 0.08 to bias the simulation calibration,<fnref n="23" /> W_c = 1 (3 × 10⁻⁵) in the Planck 2015 (otherwise), and W_z as prescribed in Section 4.1. All triangles (k₁, k₂, k₃) satisfying that B_spt agrees with Btree within 5% and up to 0.3 h Mpc⁻¹ were included, thus restricting the fitting to large scales. In the low-k limit, the ϵ term dominates the denominator of Equation (11). As k approaches 0.3 h Mpc⁻¹, the ΔB term (<0.025 × B_spt) dominates. We used the central bin values of (k₁, k₂, k₃), which were also used in the L = 4 h⁻¹ Gpc simulation with bin width Δlog₁₀k = 0.1. Figure 3 shows all triangles used in the SPT calibration at z = 0.55 (orange diamonds). The average number of triangles was 350 in each cosmological model at each redshift. Figure 4 shows the maximum and minimum k in the SPT calibration. The minimum k is 1.6 × 10⁻³ h Mpc⁻¹ for all the redshifts. The maximum k slightly increases from 0.18 to 0.28 h Mpc⁻¹ from low to high z, because the SPT approaches the tree level at higher z.

## 4.3. Results

The total χ² was computed as

$$\chi^2(\mathbf{p}) = \chi^2_{\text{sim}}(\mathbf{p}) + \chi^2_{\text{spt}}(\mathbf{p}). \quad (12)$$

The best-fitting parameters p were then numerically searched by minimizing χ². The resulting best-fit model is presented in Appendix B. The minimum was found by the downhill simplex routine (amoeba) in Numerical Recipes (Press et al. 2002).

Figure 5 plots the matter BSs computed by the tree-level formula, our fitting formula, SC01, and GM12, along with the simulation results of the Planck 2015 model with z = 0–2. From left to right, the four panels correspond to particular triangle configurations: equilateral (i.e., k₁ = k₂ = k₃), flattened (k₁ = 2k₂ = 2k₃), and two squeezed cases (k₁ = k₂ ≫ k₃ with k₃ = 0.045 and 0.45 h Mpc⁻¹). In this and the following figures, the simulation data points satisfy conditions 1–3 in Section 4.1. In the fitting formulae of SC01 and GM12, we applied the measured PS of the simulations to remove the inaccuracy of the PS appearing in these models. Figure 6 plots the ratios of SC01, GM12, and the simulation results to our fitting formula results. Clearly, our fitting formula agreed with the simulations over the tested scales, redshifts, and triangle shapes. In contrast, the previous formulae overpredicted the squeezed BS, as previously reported by Namikawa et al. (2019). The simulations performed in different box sizes were also consistent.

Figure 7 shows the ratios of the modeled and simulated BSs to the tree-level BS on quasi-nonlinear scales. On larger scales, both our simulations and fitting formula were consistent with the tree-level prediction. Meanwhile, the one-loop SPT slightly overpredicted the BS on quasi-nonlinear scales at low z (z < 1) (consistent with Figure 19 of Lazanu et al. 2016), but its inaccuracy improved at higher redshifts. In the flattened case, the SPT slightly suppressed the BS at k ~ 0.1 h Mpc⁻¹, and our model captures this trend. Some data points at k < 0.1 h Mpc⁻¹ were omitted because their relative error exceeded

<footnote n="23">²³ The resulting χ²sim was approximately 40 times larger than χ²spt.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Comparison of matter BSs obtained in *N*-body simulations and the fitting formula for the Planck 2015 best-fit ΛCDM model. The curves denote the theoretical models: our fit (BiHalofit; solid red), Gil-Marín et al. (2012) (GM12; long-dashed orange), Scoccimarro & Couchman (2001) (SC01; short-dashed brown), and the tree level (dotted purple). Symbols denote the simulation results in various box sizes: L = 4, 2, and 1 h⁻¹ Gpc (blue triangles, green triangles, and black circles) and L = 200 h⁻¹ Mpc (gray diamonds). The theoretical and simulation results are consistently binned (bin width Δ log₁₀ *k* = 0.1). Along the vertical axis, the BS is multiplied by ∝*k*ⁿ (*n* = 1, 2, or 3), as denoted above the top panels to clarify the presentation. Throughout this paper, the units of *k* and *B*(*k*₁, *k*₂, *k*₃) are *h* Mpc⁻¹ and (*h*⁻¹ Mpc)⁶, respectively.</figcaption></figure>

10%. The rms deviation between the formula and SPT was 0.96% for all the triangles in our sample of *z* = 0–10. Therefore, the accuracy reached the percent level at largest scales ($k < 0.3 \ h\ \text{Mpc}^{-1}$).

Figure 8 shows the BS ratios of the simulation results to our formula for *all* triangles satisfying conditions 1–3 in Section 4.1. There are approximately 1800 data points in each redshift. Our model agreed with the simulations within 10%
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Same as Figure 5, but relative to the red curves (BiHalofit).</figcaption></figure>

(15%) up to $k = 3$ (10) h Mpc⁻¹ for $z = 0–3$. At $z = 5$ (10), the agreement was 20% up to $k = 3$ (1) h Mpc⁻¹. The rms deviation was 2.7% (3.2%) and 3.7% (5.0%) up to $k = 3$ (10) h Mpc⁻¹ for $z = 0–3$ and $z = 0–10$, respectively. Moreover, the accuracy was independent of bin width, as confirmed by setting a narrower bin width ($\Delta\log_{10}k = 0.05$) in the same tests. The narrow bins yielded an rms deviation of 2.9% (3.4%) up to $k = 3$ (10) h Mpc⁻¹ at $z = 0–3$, quantitatively consistent with the above results. Therefore, the accuracy is approximately 3% at $k < 10$ h Mpc⁻¹ and $z = 0–3$ for most of the triangles (but it reaches 10%–15% in the worst cases).

Figure 9 plots the BS ratios of the simulations to our formula in the 40 wCDM models. In this case, as we prepared a single realization for each cosmological model, the BS measurements

had a relatively large scatter (typically 10%). All data points satisfied conditions 1–3 in Section 4.1. There are a huge number of data points (~5 × 10⁴) at each redshift. The rms deviation was 8.0% (11.2%) up to $k = 3$ (10) h Mpc⁻¹ for $z = 0–1.5$. The deviation includes the 10%-level sample variance of the simulations.

To further investigate the cosmological dependence of the accuracy, we divided the models into three groups with different ranges of $S_8 = \sigma_8(\Omega_m/0.3)^{0.5}$. The data points shown in Figure 9 are color-coded as described in the caption of Figure 1. Our formula agreed with the simulations within ~20% for $S_8 = 0.6–1.0$, but the agreement degraded outside this $S_8$ range because the fluctuation amplitude ($\sigma_8$) and the linear growth factor largely differed between these models and
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Similar to Figure 6, but relative to the tree-level results on quasi-linear scales. The dotted-dashed orange curves are the binned one-loop SPT predictions.</figcaption></figure>

the Planck 2015 model. As all cosmological models converged to the Einstein-de Sitter model at high z, the fits improved at higher redshifts.

In the 40 wCDM models, the rms deviation between the formula and the SPT is 1.3% at $k < 0.3 \ h \ \text{Mpc}^{-1}$ and $z = 0–10$. Therefore, our formula is well consistent with the SPT at largest scales.

## 5. Baryonic Effects

Our N-body simulations did not include the baryonic processes such as gas cooling, star formation, supernovae, and active galactic nucleus (AGN) feedbacks. Baryons are known to significantly affect the nonlinear PS at $k \gtrsim 1 \ h \ \text{Mpc}^{-1}$ (e.g., Semboloni et al. 2011a; van Daalen et al. 2011; Osato et al. 2015; Hellwing et al. 2016; Chisari et al. 2018, 2019). In this

section, the baryonic effects on the BS fitting formula are investigated in state-of-the-art hydrodynamic simulations using the IllustrisTNG data set<fnref n="24" /> (Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018, 2019; Pillepich et al. 2018; Springel et al. 2018). The simulations incorporate astrophysical processes in a subgrid model and thereby follow the galaxy formation and evolution processes. The IllustrisTNG project conducted three sets of simulations in different box sizes, with three mass resolutions in each box size. Here we used the highest-resolution simulation in the largest box (referred to as TNG300-1) of size $L = 205 \ h^{-1} \ \text{Mpc}$ (≈300 Mpc). This box contains 2500³ dark matter particles and the same number of baryon particles. The cosmological model of IllustrisTNG is based on the Planck 2015 best-fit ΛCDM (Planck

<footnote n="24">²⁴ http://www.tng-project.org</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** BS ratios of simulation results to BiHalofit for all triangles in the Planck 2015 model. The horizontal axis presents the maximum wavenumber of ($k_1, k_2, k_3$). Dozens of points are distributed along the vertical axis at each $k_{max}$. The horizontal positions for different $L$ are slightly offset to clarify the presentation.</figcaption></figure>

Collaboration 2016). The collaboration has released the particle positions and masses of dark matter and baryons (in the forms of gas, stars, and black holes) at $z = 0-20$. The IllustrisTNG team also performed dark-matter-only (dmo) runs. By comparing the simulations in the presence and absence of baryons, we can single out the impact of baryons on matter clustering.

To calculate the density contrast, we assigned the particle masses to 1024³ grid cells and measured the BS as described in Section 3.4. The bin width was set to $\Delta\log_{10}k = 0.1$. We calculated the BS ratio of the simulations with baryons ($B_b$) to the dmo run ($B_{dmo}$),

$$R_b(k_1, k_2, k_3) = \frac{B_b(k_1, k_2, k_3)}{B_{dmo}(k_1, k_2, k_3)} \quad (13)$$

We measured this ratio at 11 redshifts: $z = 0, 0.2, 0.4, 0.7, 1, 1.5, 2, 3, 5, 7, \text{ and } 10$.

Figure 10 plots the ratios in Equation (13) for three triangle configurations in the range of $z = 0-2$. To reduce the sample-variance scatter in the ratio at large scales, the simulations with and without baryons had the same seed in their initial conditions. The baryons suppressed the BS amplitude at $k \sim 10 h \, \text{Mpc}^{-1}$ by the AGN feedback but strongly enhanced the BS amplitude at high $k$ ($> 10 h \, \text{Mpc}^{-1}$) by the gas cooling. This trend is consistent with the PS (see also Figure 11). However, at intermediate scales ($k \sim 1-10 h \, \text{Mpc}^{-1}$) and low redshifts ($z < 1$), the baryons slightly enhanced the amplitude by $\sim 10\%$. To our knowledge, this small enhancement has not been commonly observed in PS.

Figure 11 plots the PS ratios with and without baryons, computed in TNG300-1. The circles (crosses) are the PSs of the total matter (dark matter component only) in the hydrodynamic
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Same as Figure 8, but showing the results of the 40 wCDM models with $z = 0–1.48$. The colors correspond to the cosmological models shown in Figure 1. The cyan, gray, and magenta points were obtained in models with different ranges of $S_8 \equiv \sigma_8(\Omega_m/0.3)^{0.5}$. The plus signs (crosses) are the simulation results of $L = 1 (2) h^{-1}$ Gpc. Note that every data point has an intrinsic scatter of ~10% because there is a single realization in each wCDM model.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Ratio of BS with baryons to BS without baryons, $R_b = B_b/B_{dmo}$, defined in Equation (13), measured from the TNG300-1. The filled circles are the total-matter (dark matter and baryons) BS divided by that from the dark-matter-only run. The red curves are our fit given in Appendix C. We comment that $R_b$ varies by ~10% among hydrodynamical simulations performed in different groups, because the baryonic feedback models differ (see the main text).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11. Similar to Figure 10, but ratio of P(k) measured from the TNG300-1. The filled circles (crosses) are the total-matter (dark matter) P(k) from the hydrodynamic run divided by that from the dark-matter-only run.</figcaption></figure>

run divided by that in the dmo run. At intermediate scales, the plots of the crosses were slightly enhanced, whereas those of the circles were not. The same feature is mentioned in Section 3 of Springel et al. (2018), consolidating that the enhancement source is the dark matter component. Moreover, the dark matter PS and the total-matter BS are enhanced at almost the same wavenumbers.

During the preparation of this paper, Foreman et al. (2019) posted an arXiv paper concerning the baryonic effects on BS measured in hydrodynamic simulations (including TNG300-1). They reported the same trend and clarified its cause. At late times ($z < 1$), the AGN feedback becomes less effective and the expelled gas reaccretes into a halo. Gas contraction then affects the dark matter distribution in the halo. As the BS is more sensitive to dark matter than the PS (see their Section 3.1.1), the enhancement at intermediate scales appears only in BS. By studying the baryonic effects on both PS and BS, one can discriminate among baryonic models (Semboloni et al. 2013; Foreman et al. 2019).

To incorporate the baryonic effect in our BS model, we constructed a fitting function of the ratio $R_b$ in Equation (13). The results are plotted as the solid red curves in Figure 10, and the functional form is given in Appendix C. This fitted the measurements within 7.3% (5.3%) for $k < 10 h \text{ Mpc}^{-1}$ at low (high) redshift, $z = 0-1$ (1.5–10). The rms deviation was 1.8% (2.9%) for $k < 30 h \text{ Mpc}^{-1}$ at $z = 0-1$ (1.5–10). In this data fitting, approximately 760 (8300) triangles existed at each redshift (over the full range $z = 0-10$). To include the baryonic effects, the user can simply multiply $R_b$ by the BS fitting formula. The same approach was adopted by Harnois-Déraps et al. (2015), who studied the baryonic effects on PS.

We comment that the BS ratio $R_b$ varies by approximately 10% among hydrodynamical simulations, because the baryonic feedback models differ. Foreman et al. (2019) measured the BS ratio $R_b$ for the equilateral case in four simulations: IllustrisTNG, Illustris (Vogelsberger et al. 2014), BAHAMAS (McCarthy et al. 2017), and EAGLE (Schaye et al. 2015). They

reported a 10%–20% variation in the results for $k > 1 h \text{ Mpc}^{-1}$ and $z = 0-3$. In Illustris and BAHAMAS, the small enhancement at intermediate scales ($k = 0.1-1 h \text{ Mpc}^{-1}$; see Figure 10) was absent, but the suppression at small scales ($k \gtrsim 10 h \text{ Mpc}^{-1}$) was amplified because these models implemented a stronger AGN feedback than IllustrisTNG. Therefore, the uncertainty in our fitting formula also hovered around 10%.

## 6. Comparison with Weak-lensing Simulations

Using the fitting formula of matter BS calibrated over wide ranges of wavenumbers and redshifts, we can predict the lensing observables by integrating along the line of sight. This section compares our theoretical prediction with the weak-lensing BS measured in ray-tracing simulations. We consider the convergence BS in two cases: CMB lensing (Section 6.1) and cosmic shear (Section 6.2).

The convergence field is a dimensionless matter density integrated along the line of sight toward the source. The convergence at angular position $\theta$ for a source distance $r_s$ is given by (e.g., Bartelmann & Schneider 2001)

$$ \kappa(\theta) = \int_{0}^{r_s} dr \, W(r, r_s) \delta(r\theta, r; z), \qquad (14) $$

with the weight function

$$ W(r, r_s) = \frac{3H_0^2 \Omega_m}{2c^2} \frac{r(r_s - r)}{a(r) r_s}, \qquad (15) $$

where $r(r_s)$ is the comoving distance (to the source) and $a(r)$ is the scale factor. The convergence BS is

$$ B_{\kappa}(\ell_1, \ell_2, \ell_3) = \int_{0}^{r_s} dr \frac{W^3(r, r_s)}{r^4} B\left(\frac{\ell_1}{r}, \frac{\ell_2}{r}, \frac{\ell_3}{r}; z\right), \qquad (16) $$

where $\ell_i (=k_i r)$ is the multipole moment and $B(k_1, k_2, k_3; z)$ is the matter BS at $z$. This formula was derived under the flat-sky and the Born approximations. When the source has a high redshift, the Born approximation is less accurate and must be adjusted by post-Born corrections (Pratten & Lewis 2016). These corrections are necessary only in CMB lensing (in cosmic shear, their contribution is $O(1%)$; see Figure 7 of Pratten & Lewis 2016).

For a source with a given redshift, the convergence BS is more sensitive to lower-$z$ structures than the convergence PS (see, e.g., Figure 4 of Takada & Jain 2002) because the matter BS (PS) evolves proportionally to the fourth (second) power of the linear growth factor in the linear regime. Therefore, the matter BS and PS can probe structures with different redshifts in a complimentary manner.

### 6.1. CMB Lensing

Namikawa et al. (2019) recently measured the convergence BS in full-sky light-cone simulations (Takahashi et al. 2017). Here we compare their measurements with the theoretical predictions. Takahashi et al. (2017) ran cosmological N-body simulations of the inhomogeneous mass distribution in the universe, from the present to the last scattering surface. Their cosmological model was consistent with the WMAP 9 yr result (Hinshaw et al. 2013). The authors also calculated the light-ray paths deflected by the intervening matter in a ray-tracing simulation, which tracks the trajectories of the light rays emitted from the observer’s position ($z = 0$) to the last
scattering surface (the ray-tracing scheme is detailed in Shirasaki et al. 2015). Their results included the post-Born effects. Takahashi et al. (2017) provided 108 full-sky convergence maps<sup>25</sup> based on the HEALPix pixelization with $N_{side} = 8192$ (4096), corresponding to a pixel size of 0.48 (0.96) arcmin (Górski et al. 2005). They confirmed that the convergence PS agrees with the theoretical CAMB prediction using the Halofit PS option (within 5% at $\ell \le 2000$ on the high-resolution maps with $N_{side} = 8192$).

Figure 12 plots the BS measurements obtained from the 108 maps with $N_{side} = 8192$ (Namikawa et al. 2019). The theoretical predictions were computed for the WMAP 9 yr cosmological model to be consistent with the simulations. Here the nonlinear PS for GM12, SC01, and the post-Born correction was computed by the revised Halofit. For a fair comparison, both the theoretical predictions and simulation results were binned with the same bin width ($\Delta\ell = 100$). The error bars were computed for the ideal full-sky measurement (i.e., the cosmic-variance limit) and scaled as $[(\text{survey area})/(4\pi)]^{-1/2}\Delta\ell^{-3/2}$, assuming Gaussian variance. Overall, our fitting formula better predicted the BS of CMB lensing than the previously proposed formulae. In the equilateral case, the analytical and simulated BS agreed within ~10% on most angular scales. The differences were within $0.2\sigma$ (bottom panels of Figure 12). In the flattened case, the ratio (middle panel) was far from unity because the BS approaches zero at $\ell \gtrsim 1000$. This discrepancy is approximately $0.2\sigma$ of the cosmic variance. In the squeezed and isosceles configurations, our fitting formula significantly reduced the discrepancy between the simulation result and the analytical prediction.

Although our fitting surely improved the prediction accuracy, noticeable discrepancies from the simulations were introduced by several sources. First, the finite thickness of the lens planes employed in the ray-tracing simulations may affect the simulations at $\ell < 200$ (the same effect on convergence PS is demonstrated in Figure 10 of Takahashi et al. 2017). Second, the flat-sky formula in Equation (16) is inaccurate at large angular scales (the accuracy of the flat-sky approximation in the cosmic-shear PS is detailed in Kilbinger et al. 2017; Kitching et al. 2017). For example, in the squeezed limit, the minimum multipole is fixed as $\ell_3 = 50$, but a larger $\ell_3$ can mitigate the discrepancy (Namikawa et al. 2019). Since reducing the finite thickness of lens planes requires more numerically expensive simulations, we will leave the detailed study for future work.

## 6.2. Cosmic Shear

Let us now consider the cosmic-shear signals in galaxy-shape measurements, which probe lower redshifts than CMB lensing. Sato et al. (2009) ran cosmological N-body simulations and subsequent ray-tracing simulations under the flat-sky approximation. Despite their small field of view ($5 \times 5$ deg$^2$), they acquired sufficiently many weak-lensing maps (1000) for an accurate BS measurement. Their cosmological model was consistent with the WMAP 3 yr result (Spergel et al. 2007).

Figure 13 plots the convergence BS at a source redshift of $z_s = 1$, measured from the 1000 maps by Kayo et al. (2013).

The theoretical and simulation results were consistently binned with $\Delta \log_{10} \ell = 0.13$. The simulation results were valid (within 5% error) up to $\ell \simeq 4000$, as confirmed by comparing the results with those from low- and high-resolution maps (Sato et al. 2009; Valageas et al. 2012). Overall, the plot shows a similar trend to the matter BS at $z = 0.55$ (see Figure 5). Our fitting formula well agreed with the simulation (within 10% level up to $\ell = 4000$). The deviation at small scales ($\ell \gtrsim 4000$) was attributed to the limited resolution of the simulation.

As the cosmological models in this and the previous Sections 6.1 and 6.2 differ from the Planck 2015, the agreement with the weak-lensing simulations provides a nontrivial validation of our formula for other cosmological models.

# 7. Discussion

## 7.1. Systematics in CMB Lensing

In CMB lensing measurements, the lensing map is reconstructed through mode-mixing of the CMB anisotropies induced by lensing (Hu & Okamoto 2002). Therefore, any other sources of mode-mixing can bias the lensing measurements and hence the BS of CMB lensing. Bias can be sourced from instrumentation factors such as masking, inhomogeneous noise, beam, and point sources (Hanson et al. 2009; Namikawa et al. 2013) and from extragalactic foregrounds such as the thermal Sunyaev–Zel’dovich effect, the cosmic infrared background (Osborne et al. 2014; van Engelen et al. 2014; Madhavacheril & Hill 2018), and its lensing (Mishra & Schaan 2019). Calibration uncertainties in the CMB map are also important sources of systematic error, because when the lensing reconstruction is performed by a quadratic estimator, the measured BS depends on the sixth power of the map-calibration uncertainties. In contrast, the lensing PS depends on the fourth power of the map-calibration uncertainties. Combining the BS and PS is expected to constrain the bias contributed by the instrumental uncertainties and astrophysical sources, because these sources affect the spectra in different ways. A joint analysis of the PS and BS is therefore crucial for a robust cosmological analysis in future CMB experiments.

## 7.2. Intrinsic Alignment

The intrinsic alignment (IA) of galaxies is a major systematic error in cosmic shear (reviewed by Joachimi et al. 2015; Troxel & Ishak 2015). A massive structure near the source galaxy exerts a tidal force that distorts the shapes and contaminates the lensing signal. Approximately 10% of the cosmic-shear BS is contaminated by this mechanism (Semboloni et al. 2008). Several authors have proposed methods for mitigating or removing the contamination from the signal (Shi et al. 2010; Troxel & Ishak 2012). The combined PS and BS can strongly constrain not only the cosmological parameters but also the IA.

## 7.3. Bispectrum Covariance

Thus far, we have not discussed the modeling of BS covariance, which is another important ingredient of cosmological likelihood analysis. The BS covariance of Gaussian fluctuations has a simple form given by the PS and the shot noise (Sefusatti et al. 2006). However, in the nonlinear regime, one should consider the non-Gaussian and supersample contributions (e.g., Takada & Hu 2013), which complicate the evaluation. In such cases, the covariance has been estimated

<footnote n="25"><sup>25</sup> These maps are available at http://cosmo.phys.hirosaki-u.ac.jp/takahasi/allskyINI.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Convergence BSs measured from simulation maps of CMB lensing. The black symbols were averaged from 108 full-sky maps (Takahashi et al. 2017; Namikawa et al. 2019). The error bars are the standard deviations scaled by [(survey area)/(4π)]⁻¹/². The solid red, dashed orange, and dotted brown curves are the theoretical predictions based on BiHalofit, GM12, and SC01, respectively. The middle panels plot the relative deviations from the red curves: ΔB/B ≡ B<sub>κ</sub>/B<sub>κ</sub><sup>BiHalofit</sup> − 1. In the bottom panels, these deviations are further divided by the relative standard deviation (σ/B).</figcaption></figure>

by perturbation theory (e.g., Sugiyama et al. 2019), the halo
model (e.g., Kayo et al. 2013; Rizzato et al. 2019), and an
ensemble of simulation mocks (e.g., Sato & Nishimichi 2013;
Chan & Blot 2017; Chan et al. 2018; Colavincenzo et al. 2019).
To estimate an unbiased inverse covariance, the last approach
should generate more mocks than a number of k-bins (e.g.,
Hartlap et al. 2007); consequently, the number of mocks can be
huge (>10<sup>-3</sup>). This topic is reserved for future work.

7.4. Emulator

Several groups are developing nonlinear PS emulators that
interpolate simulation results over a wide range of wavenum-
bers, redshifts, and cosmological models (Lawrence et al. 2017;
Garrison et al. 2018; DeRose et al. 2019; Knabenhans et al.
2019; Nishimichi et al. 2019). We expect that developing a
similar emulator for BS is much more formidable, for two
reasons. First, we measure a *binned* BS but require an *unbinned*
BS (recall that BS is sensitive to binning). Therefore, we cannot
simply interpolate the measured quantities. Second, BS
measurements have larger sample variances than PS measure-
ments, which demand many realizations in each cosmological
model. This is computationally expensive.

8. Conclusions

We have constructed a fitting formula of the matter BS
calibrated in high-resolution *N*-body simulations of 41 *w*CDM
models around the Planck 2015 best-fit ΛCDM model. The
calibration covers a wide range of wavenumbers (up to *k* =
30 *h* Mpc⁻¹) and redshifts (*z* = 0–10) for the Planck 2015

model. The 40 wCDM models supplement the calibration
at z = 0–1.5. We also performed a large-scale calibration
using perturbation theory for all the cosmological models
(at k < 0.3 h Mpc⁻¹ and z = 0–10). The calibration range is
summarized in Table 2. The simulation boxes are sufficiently
large (side length L = 1, 2, and 4 h⁻¹ Gpc) to cover almost all
triangles (k₁, k₂, k₃) measured in forthcoming weak-lensing
surveys and CMB lensing experiments. The accuracy was
within 10% (15%) up to k = 3 (10) h Mpc⁻¹ in the redshift
range z = 0–3 for the Planck 2015 model. The rms deviation
was 2.7% (3.2%) up to k = 3 (10) h Mpc⁻¹ for z = 0–3.
Therefore, the accuracy was approximately 3% for most of
the triangles and 10%–15% only for the worst cases. Mean-
while, the accuracy of the 40 wCDM models was around 20%
for k < 3 h Mpc⁻¹ and z = 0–1.5. In these models, a 10%
intrinsic scatter was introduced to the simulation data by the
single realization. The rms deviation was 8.0% (11.2%) up to
k = 3 (10) h Mpc⁻¹ for z = 0–1.5. The user can easily
incorporate the baryonic effects (calibrated using IllustrisTNG)
into the fitting formula. We also confirmed that the formula
reproduces the weak-lensing convergence BS measured in
light-cone simulations.

The σ<sub>8</sub> inferred from the Planck results is larger than that
estimated from cosmic shear and galaxy–galaxy lensing σ<sub>8</sub>
(σ<sub>8</sub> ~ 0.81 vs. ~0.77; e.g., MacCrann et al. 2015; Abbott et al.
2018; Planck Collaboration 2018b). Combining weak-lensing
PS and BS can tighten the σ<sub>8</sub> constraint by a factor of 1.6–3
(e.g., Takada & Jain 2004; Kayo & Takada 2013), providing
new clues for solving this controversy.
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Convergence BSs measured from 1000 simulation maps at the source redshift $z_s = 1$ (Sato et al. 2009; Kayo et al. 2013). Black circles are the measured averages, and the red curves were predicted by our model. The dashed orange and brown dotted curves were obtained by GM12 and SC01, respectively. The field of view of each map is 5 × 5 deg², and the error bars are the standard deviations scaled by [(survey are)/(25 deg²)]⁻¹/². The bottom panels plot the relative deviations from the red curves.</figcaption></figure>

<table><caption>**Table 2**
Calibration Range</caption><thead><tr><th>Cosmological Model</th><th>Minimum k (h Mpc<sup>-1</sup>)</th><th>Maximum k (h Mpc<sup>-1</sup>)</th><th>Redshift</th><th>Calibration</th></tr></thead><tbody><tr><td>Planck 2015 ΛCDM</td><td>1.6 × 10<sup>-3</sup></td><td>14–48</td><td>0–10</td><td>sim. and SPT</td></tr><tr><td rowspan="2">40 wCDM</td><td>1.6 × 10<sup>-3</sup></td><td>14–28</td><td>0–1.5</td><td>sim. and SPT</td></tr><tr><td>1.6 × 10<sup>-3</sup></td><td>0.18 – 0.28</td><td>2–10</td><td>SPT</td></tr></tbody></table>

**Note.** Calibration range of *k* and *z* in the simulations and the one-loop SPT. The maximum *k* slightly depends on redshift (see Figure 4). In the 40 wCDM models at *z* = 2–10, the calibration was done only by SPT.

We thank the IllustrisTNG team for making their simulations publicly available. This work was in part supported by Grant-in-Aid for Scientific Research from the Japan Society for the Promotion of Science (JSPS) (Nos. JP16H03977, JP17H01131, JP17K14273, JP19H00677, and JP19K14767) and MEXT Grant-in-Aid for Scientific Research on Innovative Areas (Nos. JP15H05893, JP15H05896, JP15H05889, JP18H04358, and JP20H04723). K.O. is supported by JSPS Overseas Research Fellowships. T.N. is supported by Japan Science and Technology Agency CREST JPMHCR1414. Numerical computations were in part carried out on Cray XC30 and XC50 at Centre for Computational Astrophysics, National Astronomical Observatory of Japan.

## Appendix A
Halo Model

The halo model, which assumes that all matter is confined in halos, is widely applied in nonlinear BS estimation (e.g., Cooray & Hu 2001; Cooray & Sheth 2002; Valageas & Nishimichi 2011;

Kayo et al. 2013; Yamamoto et al. 2017). The basic properties of a halo of mass *M* are characterized by the mass function $dn(M)/dM$, the spherical density profile $\rho(r; M)$, and the first- and second-order halo biases $b_{1,2}(M)$. This model decomposes the matter BS into three terms: one- (1h), two- (2h), and three-halo (3h) terms. The 1h and 3h terms dominate at small and large scales, respectively. The 2h term fills the gap between the 1h and 3h terms and only minimally contributes at intermediate scales, except at the squeezed limit. The BS is given by

$$B(k_1, k_2, k_3) = B_{1h}^{\mathrm{HM}}(k_1, k_2, k_3) + B_{2h}^{\mathrm{HM}}(k_1, k_2, k_3) + B_{3h}^{\mathrm{HM}}(k_1, k_2, k_3). \quad (A1)$$

The 1h term comes from the density profile of a single halo:

$$B_{1h}^{\mathrm{HM}}(k_1, k_2, k_3) = \int dM \frac{dn(M)}{dM} \left(\frac{M}{\bar{\rho}}\right)^3 \times u(k_1; M)u(k_2; M)u(k_3; M), \quad (A2)$$

where $\bar{\rho}$ is the cosmic mean density and $u(k; M)$ is the Fourier transform of the scaled density profile $\rho(r; M)/M$. The 2h term describes the correlation among two points in the same halo and a third point in another halo:

$$B_{2h}^{\mathrm{HM}}(k_1, k_2, k_3) = I_2^1(k_1, k_2)I_1^1(k_3)P_L(k_3) + 2 \mathrm{perm.}, \quad (A3)$$

with

$$I_2^1(k_1, k_2) = \int dM \frac{dn(M)}{dM} \left(\frac{M}{\bar{\rho}}\right)^2 b_1(M)u(k_1; M)u(k_2; M). \quad (A4)$$
The 3h term describes the spatial correlation among three different halos:

$$
\begin{align}
B_{3h}^{\text{HM}}(k_1, k_2, k_3) &= I_1^1(k_1)I_1^1(k_2)I_1^1(k_3)B_{\text{tree}}(k_1, k_2, k_3) \nonumber \\
&\quad + [I_1^1(k_1)I_1^1(k_2)I_1^2(k_3)P_L(k_1)P_L(k_2) + 2 \text{perm.}] \nonumber \\
&= 2\left[F_2(\mathbf{k}_1, \mathbf{k}_2) + \frac{I_1^2(k_3)}{2I_1^1(k_3)}\right]I_1^1(k_1)I_1^1(k_2)I_1^1(k_3) \nonumber \\
&\quad \times P_L(k_1)P_L(k_2) + 2 \text{perm.}, \tag{A5}
\end{align}
$$

with

$$
I_1^\beta(k) = \int dM \frac{dn(M)}{dM} \frac{M}{\bar{\rho}} b_\beta(M) u(k; M). \quad (A6)
$$

The 2h and 3h terms are proportional to $P_L$ and $(P_L)^2$, respectively.

**Appendix B**

**Fitting Formula**

Our fitting formula adopts the Halofit parameterization for nonlinear PS (Smith et al. 2003). The dimensionless linear PS is defined as $\Delta_L^2(k) = k^3 P_L(k)/(2\pi^2)$. The nonlinear scale $k_{\text{NL}}^{-1}$ is determined as

$$
\sigma^2(k_{\text{NL}}^{-1}) = 1, \quad \sigma^2(R) = \int d \ln k \, \Delta_L^2(k) \, e^{-k^2 R^2}. \quad (B1)
$$

The effective spectral index at $k_{\text{NL}}$ is defined as

$$
n_{\text{eff}} + 3 = - \frac{d \ln \sigma^2(R)}{d \ln R} \bigg|_{R=k_{\text{NL}}^{-1}} . \quad (B2)
$$

We also introduce a scaled wavenumber, $q_i = k_i/k_{\text{NL}}$ ($i=1, 2$, and 3). Note that the quantities $k_{\text{NL}}$ and $n_{\text{eff}}$ are evaluated at a given redshift. Identical parameters were defined in Smith et al. (2003).

The fitting function is the sum of the 1h and 3h terms:

$$
B(k_1, k_2, k_3) = B_{1h}(k_1, k_2, k_3) + B_{3h}(k_1, k_2, k_3). \quad (B3)
$$

The 1h term is

$$
B_{1h}(k_1, k_2, k_3) = \prod_{i=1}^{3} \left[ \frac{1}{a_n q_i^{\alpha_n} + b_n q_i^{\beta_n} \frac{1}{1 + (c_n q_i)^{-1}}} \right]. \quad (B4)
$$

Here $B_{1h}$ is assumed as the product of identical functions of $q_1$, $q_2$, and $q_3$. Similarly, the halo model $B_{1h}^{\text{HM}}$ given by Equation (A2) is the product of $u(k_i)$ terms. The 3h term is given by

$$
\begin{equation}
\begin{split}
B_{3h}(k_1, k_2, k_3) = {}& 2[F_2(\mathbf{k}_1, \mathbf{k}_2) + d_n q_3]I(k_1) \\
& \times I(k_2)I(k_3)P_E(k_1)P_E(k_2) + 2 \text{perm.},
\end{split}
\tag{B5}
\end{equation}
$$

with

$$
\begin{equation}
\begin{split}
P_E(k) = \frac{1 + f_n q^2}{1 + g_n q + h_n q^2} P_L(k) &+ \frac{1}{m_n q^{\mu_n} + n_n q^{\nu_n}} \\
& \times \frac{1}{1 + (p_n q)^{-3}}, \quad I(k) = \frac{1}{1 + e_n q}.
\end{split}
\tag{B6}
\end{equation}
$$

Here $P_E(k)$ defines the “enhanced” PS, obtained by adding a small-scale enhancement to the linear PS. The first (second) term of $P_E$ is similar to the 2h (1h) term of Halofit for the nonlinear PS. Similarly, $I(k)$ and $d_n q$ correspond to $I_1^1(k)$ and

$I_1^2(k)/[2I_1^1(k)]$ in the halo model, respectively. This 3h term approaches the tree level in the low-$k$ limit.

The 3h term $B_{3h}$ includes the 2h contribution $B_{2h}^{\text{HM}}$ in the halo model, as discussed below. As $B_{2h}^{\text{HM}}$ is proportional to $P_L$, several terms proportional to $P_L$ in $B_{3h}$ correspond to $B_{2h}^{\text{HM}}$. The enhanced PS can be decomposed into the linear PS and the small-scale enhancement: $P_E(k) \simeq P_L(k) + P_E^{\text{2nd}}(k)$ (where the prefactor of $P_L$ in Equation (B6) is ignored). The terms proportional to $P_L$ are given from Equation (B5) by

$$
\begin{equation}
\begin{split}
B_{3h}(k_1, k_2, k_3)|_{\propto P_L} \approx {}& 2\{F_2(\mathbf{k}_1, \mathbf{k}_3) + d_n q_2\} P_E^{\text{2nd}}(k_1) \\
& + (\mathbf{k}_1 \leftrightarrow \mathbf{k}_2)]I(k_1)I(k_2)I(k_3)P_L(k_3) + 2 \text{perm.}
\end{split}
\tag{B7}
\end{equation}
$$

Therefore, $2\{[F_2(\mathbf{k}_1, \mathbf{k}_3) + d_n q_2] P_E^{\text{2nd}}(k_1) + (\mathbf{k}_1 \leftrightarrow \mathbf{k}_2)]I(k_1)I(k_2)$ in the above equation corresponds to $I_2^1(k_1, k_2)$ in $B_{2h}^{\text{HM}}$. Equation (B7) enhances the squeezed $B_{3h}$ at intermediate scales (see also Figure C1).

The above fitting parameters ($a_n$, $b_n$,...) are polynomials in terms of $n_{\text{eff}}$ and $\log_{10} \sigma_8$, where $\sigma_8$ is the spherical overdensity at a radius of $8 h^{-1}$ Mpc at redshift $z$ (i.e., $\sigma_8(z=0)$ multiplied by the linear growth factor). The amplitude-determining parameters (i.e., $a_n$, $b_n$, $m_n$, and $n_n$) are functions of $\log_{10} \sigma_8$, whereas most of the other parameters are functions of $n_{\text{eff}}$. As $B(k_1, k_2, k_3)$ and $P(k)$ have dimensions of $[L^6]$ and $[L^3]$, respectively, $a_n$ and $b_n$ have dimensions of $[L^{-2}]$, $m_n$ and $n_n$ have dimensions of $[L^{-3}]$, and all other parameters are dimensionless. Here the length unit is chosen as $[L] = [h^{-1} \text{Mpc}]$.

The fitting parameters of the 1h term are given by

$$
\begin{align*}
\log_{10} a_n &= -2.167 - 2.944 \log_{10} \sigma_8 - 1.106 (\log_{10} \sigma_8)^2 \\
&\quad - 2.865 (\log_{10} \sigma_8)^3 - 0.310 r_1^{\gamma_n}, \\
\log_{10} b_n &= -3.428 - 2.681 \log_{10} \sigma_8 + 1.624 (\log_{10} \sigma_8)^2 \\
&\quad - 0.095 (\log_{10} \sigma_8)^3,
\end{align*}
$$

$$
\begin{align*}
\log_{10} b_n &= -3.428 - 2.681 \log_{10} \sigma_8 + 1.624 (\log_{10} \sigma_8)^2 \\
&\quad - 0.095 (\log_{10} \sigma_8)^3,
\end{align*}
$$

$$
\log_{10} c_n = 0.159 - 1.107 n_{\text{eff}},
$$

$$
\log_{10} \alpha_n = \min \left[ -4.348 - 3.006 n_{\text{eff}} - 0.5745 n_{\text{eff}}^2 + 10^{-0.9+0.2 n_{\text{eff}}} r_2^2, \log_{10} \left(1 - \frac{2}{3} n_s\right) \right],
$$

$$
\begin{align*}
\log_{10} \beta_n = {}& -1.731 - 2.845 n_{\text{eff}} - 1.4995 n_{\text{eff}}^2 \\
& - 0.2811 n_{\text{eff}}^3 + 0.007 r_2,
\end{align*}
$$

$$
\log_{10} \gamma_n = 0.182 + 0.570 n_{\text{eff}},
$$

(B8)

where $r_{1,2}$ are ratios of the minimum ($k_{\text{min}}$) and middle ($k_{\text{mid}}$) wavenumbers to the maximum ($k_{\text{max}}$) wavenumber of the triangle, respectively, given by

$$
r_1 = \frac{k_{\text{min}}}{k_{\text{max}}}, \quad r_2 = \frac{k_{\text{mid}} + k_{\text{min}} - k_{\text{max}}}{k_{\text{max}}}. \qquad (B9)
$$

The $r_{1,2}$ represents a “halo triaxiality” in the 1h term (Smith et al. 2006): $r_{1,2} \to 0$ in the squeezed case ($k_{\text{min}} \ll k_{\text{mid}} \simeq k_{\text{max}}$), $r_1 (r_2) \to 0.5 (0)$ in the flattened case ($k_{\text{min}} \simeq k_{\text{mid}} \simeq k_{\text{max}}/2$), and $r_{1,2} \to 1$ in the equilateral case ($k_{\text{min}} \simeq k_{\text{mid}} \simeq k_{\text{max}}$). These terms slightly enhance (suppress) the squeezed (equilateral) BS at $k \gtrsim 5 h \text{Mpc}^{-1}$. To ensure that the 1h term is smaller than the tree level in the low-$k$ limit, the maximum $\alpha_n$ was set to $\alpha_{n,\text{max}} = 1 - (2/3)n_s$ (where $n_s$ is the spectral index of the initial PS). The parameters
<figure><img src="image_15.png" /><figcaption>Figure C1. Contributions of the one-halo (1h) and three-halo (3h) terms to the total-matter BS in the fitting formula.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure C2. Effect of k-binning on the BS fitting formula. The black curve is the unbinned result, and the dashed red curve is the binned result with a bin width of $\Delta \log_{10} k = 0.1$.</figcaption></figure>

of the 3h term are given by

$$
\begin{align*}
\log_{10} f_n &= -10.533 - 16.838 n_{\text{eff}} - 9.3048 n_{\text{eff}}^2 - 1.8263 n_{\text{eff}}^3, \\
\log_{10} g_n &= 2.787 + 2.405 n_{\text{eff}} + 0.4577 n_{\text{eff}}^2, \\
\log_{10} h_n &= -1.118 - 0.394 n_{\text{eff}}, \\
\log_{10} m_n &= -2.605 - 2.434 \log_{10} \sigma_8 + 5.710(\log_{10} \sigma_8)^2, \\
\log_{10} n_n &= -4.468 - 3.080 \log_{10} \sigma_8 + 1.035(\log_{10} \sigma_8)^2, \\
\log_{10} \mu_n &= 15.312 + 22.977 n_{\text{eff}} + 10.9579 n_{\text{eff}}^2 + 1.6586 n_{\text{eff}}^3, \\
\log_{10} \nu_n &= 1.347 + 1.246 n_{\text{eff}} + 0.4525 n_{\text{eff}}^2, \\
\log_{10} p_n &= 0.071 - 0.433 n_{\text{eff}}, \\
\log_{10} d_n &= -0.483 + 0.892 \log_{10} \sigma_8 - 0.086 \Omega_m, \\
\log_{10} e_n &= -0.632 + 0.646 n_{\text{eff}}.
\end{align*}
\tag{B10}
$$

Here $\Omega_m$ is the matter density parameter at z. Note again that $a_n$ and $b_n$ ($m_n$ and $n_n$) have the units of $h^2 \text{Mpc}^{-2}$ ($h^3 \text{Mpc}^{-3}$) and the other parameters are dimensionless. As the calibration was performed from $z=0$ to 10, the formula should be switched to the tree level at $z>10$.

We checked that the above fitting parameters ($a_n$, $b_n$, ...) did not depend on other cosmological parameters as follows: we fitted the formula to the Planck 2015 model at each redshift and to each wCDM model at $z=0$: then, it turned out that the best-fit values of them ($a_n$, $b_n$, ...) mainly depend on two parameters

of $n_{\text{eff}}$ and $\sigma_8$ and did not correlate with the other parameters (i.e., cosmological parameters and redshift).

Figure C1 plots the separate contributions of $B_{1h}$ and $B_{3h}$ at $z=0.55$ in the Planck 2015 model. The results are unbinned. In the equilateral and flattened cases, the 1h (3h) term clearly dominated at small (large) scales. In the squeezed case, the 1h (3h) term dominated in the nonlinear (linear) regime of $k_3$. In all cases, the second term of $P_E$ enhances the 3h term at intermediate scales ($k \sim 1-10 h \text{Mpc}^{-1}$).

Figure C2 shows the binning effect on BS. The binning affected the squeezed BS, because the cosine term in the $F_2$ kernel is very sensitive to the squeezed triangle configuration (for details, see Section IIB of Namikawa et al. 2019).

## Appendix C

### Fitting the Ratio of the Bispectrum with Baryons to That without Baryons

This appendix fits the ratio of the BS with baryons to that without baryons. The ratio, defined as $R_b$ in Equation (13), was calibrated in the TNG300-1 simulation (Nelson et al. 2019). The analysis included all triangle configurations ($k_1$, $k_2$, $k_3$) satisfying the following two conditions: (a) the number of triangles in the bin exceeds $10^6$ to remove noisy data points, and (b) the shot-noise contribution is less than 3%. The fitting
range was $k = 0.03–30 h \, Mpc^{-1}$ and $z = 0–10$ (11 redshifts of $z = 0, 0.2, 0.4, 0.7, 1, 1.5, 2, 3, 5, 7, \text{ and } 10$).

For $z \le 5$, the BS ratio $R_b$ was given by

$$R_b(k_1, k_2, k_3) = \prod_{i=1}^{3} \left[ A_0 \exp\left\{-\left|\frac{x_i - \mu_0}{\sigma_0}\right|^{\alpha_0}\right\} - A_1 \exp\left\{-\left(\frac{x_i - \mu_1}{\sigma_1}\right)^2\right\} + \left\{\left(\frac{k_i}{k_*}\right)^{\alpha_2} + 1\right\}^{\beta_2} \right], \quad (C1)$$

where $x_i = \log_{10}[k_i/(h \, Mpc^{-1})]$. The fitting parameters are the following functions of the scale factor $a$:

$$\begin{aligned} A_0 &= 0.068(a - 0.5)^{0.47}\Theta(a - 0.5), \\ \mu_0 &= 0.018 a + 0.837 a^2, \\ \sigma_0 &= 0.881 \mu_0, \\ \alpha_0 &= 2.346, \\ A_1 &= 1.052(a - 0.2)^{1.41}\Theta(a - 0.2), \\ \mu_1 &= |0.172 + 3.048 a - 0.675 a^2|, \\ \sigma_1 &= (0.494 - 0.039 a) \mu_1, \\ k_* &= 29.90 - 38.73 a + 24.30 a^2, \\ \alpha_2 &= 2.25, \\ \alpha_2\beta_2 &= \frac{0.563}{(a/0.06)^{0.02} + 1}, \end{aligned} \quad (C2)$$

where $k_*$ has units of $h \, Mpc^{-1}$ and $\Theta(x)$ is the step function: $\Theta(x) = 1$ and $0$ for $x \ge 0$ and $x < 0$, respectively. The first term of Equation (C1) represents the small enhancement at intermediate scales ($k \approx 1–10 h \, Mpc^{-1}$) and low $z$ ($<1$), the second term is the depression at $k \approx 10 h \, Mpc^{-1}$, and the last term is the strong enhancement at high $k$ ($\gtrsim 10 h \, Mpc^{-1}$). The ratio $R_b$ approaches unity in the low-$k$ limit.

At higher redshifts ($z > 5$; $z = 7$ and $10$ in our data), $R_b = 1$ was a good approximation in our fitting range ($k < 30 h \, Mpc^{-1}$) because the effects of AGN feedback and star formation on BS were suppressed at high $z$.

## ORCID iDs

Ryuichi Takahashi https://orcid.org/0000-0001-6021-0147

Takahiro Nishimichi https://orcid.org/0000-0002-9664-0760

Atsushi Taruya https://orcid.org/0000-0002-4016-1955

Ken Osato https://orcid.org/0000-0002-7934-2569

## References

Abbott, T. M. C., Abdalla, F. B., Alarcon, A., et al. 2018, PhRvD, 98, 043526

Angulo, R. E., Foreman, S., Schmittfull, M., & Senatore, L. 2015, JCAP, 10, 039

Bartelmann, M., & Schneider, P. 2001, PhR, 340, 291

Beck, D., Fabbian, G., & Errard, J. 2018, PhRvD, 98, 043512

Bergé, J., Amara, A., & Réfrégier, A. 2010, ApJ, 712, 992

Bernardeau, F., Colombi, S., Gaztañaga, E., & Scoccimarro, R. 2002a, PhR, 367, 1

Bernardeau, F., Mellier, Y., & van Waerbeke, L. 2002b, A&A, 389, L28

Böhm, V., Schmittfull, M., & Sherwin, B. D. 2016, PhRvD, 94, 043519

Bose, B., Byun, J., Lacasa, F., Moradinezhad Dizgah, A., & Lombriser, L. 2020, JCAP, 02, 025

Bose, B., & Taruya, A. 2018, JCAP, 10, 019

Byun, J., Eggemeier, A., Regan, D., Seery, D., & Smith, R. E. 2017, MNRAS, 471, 1581

Chan, K. C., & Blot, L. 2017, PhRvD, 96, 023528

Chan, K. C., Moradinezhad Dizgah, A., & Noreña, J. 2018, PhRvD, 97, 043532

Chang, C., Pujol, A., Mawdsley, B., et al. 2018, MNRAS, 475, 3165

Chisari, N. E., Mead, A. J., Joudaki, S., et al. 2019, OJAp, 2, 4

Chisari, N. E., Richardson, M. L. A., Devriendt, J., et al. 2018, MNRAS, 480, 3962

Colavincenzo, M., Sefusatti, E., Monaco, P., et al. 2019, MNRAS, 482, 4883

Cooray, A., & Hu, W. 2001, ApJ, 548, 7

Cooray, A., & Sheth, R. 2002, PhR, 372, 1

Coulton, W. R., Liu, J., Madhavacheril, M. S., Böhm, V., & Spergel, D. N. 2019, JCAP, 05, 043

Crocce, M., Pueblas, S., & Scoccimarro, R. 2006, MNRAS, 373, 369

DeRose, J., Wechsler, R. H., Tinker, J. L., et al. 2019, ApJ, 875, 69

Fabbian, G., Lewis, A., & Beck, D. 2019, JCAP, 10, 057

Foreman, S., Coulton, W., Villaescusa-Navarro, F., & Barreira, A. 2019, arXiv:1910.03597

Fry, J. N., & Gaztanaga, E. 1993, ApJ, 413, 447

Fu, L., Kilbinger, M., Erben, T., et al. 2014, MNRAS, 441, 2725

Garrison, L. H., Eisenstein, D. J., Ferrer, D., et al. 2018, ApJS, 236, 43

Gatti, M., Chang, C., Friedrich, O., et al. 2019, arXiv:1911.05568

Gil-Marín, H., Percival, W. J., Cuesta, A. J., et al. 2016, MNRAS, 460, 4210

Gil-Marín, H., Wagner, C., Fragkoudi, F., Jimenez, R., & Verde, L. 2012, JCAP, 2, 047

Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

Groth, E. J., & Peebles, P. J. E. 1977, ApJ, 217, 385

Hamana, T., Shirasaki, M., Miyazaki, S., et al. 2020, PASJ, 72, 16

Hanson, D., Rocha, G., & Gorski, K. 2009, MNRAS, 400, 2169

Harnois-Déraps, J., van Waerbeke, L., Viola, M., & Heymans, C. 2015, MNRAS, 450, 1212

Hartlap, J., Simon, P., & Schneider, P. 2007, A&A, 464, 399

Hashimoto, I., Rasera, Y., & Taruya, A. 2017, PhRvD, 96, 043526

Hearin, A. P., Zentner, A. R., & Ma, Z. 2012, JCAP, 4, 034

Heitmann, K., Higdon, D., White, M., et al. 2009, ApJ, 705, 156

Hellwing, W. A., Schaller, M., Frenk, C. S., et al. 2016, MNRAS, 461, L11

Hikage, C., Oguri, M., Hamana, T., et al. 2019, PASJ, 71, 43

Hinshaw, G., Larson, D., Komatsu, E., et al. 2013, ApJS, 208, 19

Hu, W., & Okamoto, T. 2002, ApJ, 574, 566

Huterer, D., & Takada, M. 2005, ApH, 23, 369

Jarvis, M., Bernstein, G., & Jain, B. 2004, MNRAS, 352, 338

Jenkins, A., Frenk, C. S., Pearce, F. R., et al. 1998, ApJ, 499, 20

Jing, Y. P. 2005, ApJ, 620, 559

Jing, Y. P., & Börner, G. 1998, ApJ, 503, 37

Joachimi, B., Cacciato, M., Kitching, T. D., et al. 2015, SSRv, 193, 1

Kayo, I., Suto, Y., Nichol, R. C., et al. 2004, PASJ, 56, 415

Kayo, I., & Takada, M. 2013, arXiv:1306.4684

Kayo, I., Takada, M., & Jain, B. 2013, MNRAS, 429, 344

Kilbinger, M., Heymans, C., Asgari, M., et al. 2017, MNRAS, 472, 2126

Kilbinger, M., & Schneider, P. 2005, A&A, 442, 69

Kitching, T. D., Alsing, J., Heavens, A. F., et al. 2017, MNRAS, 469, 2737

Knabenhans, M., Stadel, J., Marelli, S., et al. 2019, MNRAS, 484, 5509

Lawrence, E., Heitmann, K., Kwan, J., et al. 2017, ApJ, 847, 50

Lazanu, A., Giannantonio, T., Schmittfull, M., & Shellard, E. P. S. 2016, PhRvD, 93, 083517

Lazanu, A., & Liguori, M. 2018, JCAP, 04, 055

Lewis, A., & Challinor, A. 2006, PhR, 429, 1

Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473

Lewis, A., & Pratten, G. 2016, JCAP, 12, 003

MacCrann, N., Zuntz, J., Bridle, S., Jain, B., & Becker, M. R. 2015, MNRAS, 451, 2877

Madhavacheril, M. S., & Hill, J. C. 2018, PhRvD, 98, 023534

Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113

Marozzi, G., Fanizza, G., Di Dio, E., & Durrer, R. 2018, PhRvD, 98, 023535

Matarrese, S., Verde, L., & Heavens, A. F. 1997, MNRAS, 290, 651

McCarthy, I. G., Schaye, J., Bird, S., & Le Brun, A. M. C. 2017, MNRAS, 465, 2936

McCullagh, N., Jeong, D., & Szalay, A. S. 2016, MNRAS, 455, 2945

Mead, A. J., Peacock, J. A., Heymans, C., Joudaki, S., & Heavens, A. F. 2015, MNRAS, 454, 1958

Mishra, N., & Schaan, E. 2019, PhRvD, 100, 123504

Munshi, D., Namikawa, T., Kitching, T. D., et al. 2020, MNRAS, 493, 3985

Munshi, D., Smidt, J., Heavens, A., Coles, P., & Cooray, A. 2011, MNRAS, 411, 2241

Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206

Namikawa, T. 2016, PhRvD, 93, 121301

Namikawa, T., Bose, B., Bouchet, F. R., Takahashi, R., & Taruya, A. 2019, PhRvD, 99, 063511
Namikawa, T., Hanson, D., & Takahashi, R. 2013, MNRAS, 431, 609

Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624

Nelson, D., Springel, V., Pillepich, A., et al. 2019, ComAC, 6, 2

Nishimichi, T., Kayo, I., Hikage, C., et al. 2007, PASJ, 59, 93

Nishimichi, T., Shirata, A., Taruya, A., et al. 2009, PASJ, 61, 321

Nishimichi, T., Takada, M., Takahashi, R., et al. 2019, ApJ, 884, 29

Osato, K., Shirasaki, M., & Yoshida, N. 2015, ApJ, 806, 186

Osborne, S. J., Hanson, D., & Dore, O. 2014, JCAP, 03, 024

Peebles, P. J. E., & Groth, E. J. 1975, ApJ, 196, 1

Petri, A., Liu, J., Haiman, Z., et al. 2015, PhRvD, 91, 103511

Pillepich, A., Nelson, D., Hernquist, L., et al. 2018, MNRAS, 475, 648

Planck Collaboration 2016, A&A, 594, A13

Planck Collaboration 2018a, arXiv:1807.06210

Planck Collaboration 2018b, arXiv:1807.06209

Planck Collaboration 2019, arXiv:1905.05697

Pratten, G., & Lewis, A. 2016, JCAP, 08, 047

Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2002, Numerical Recipes in C++ : the Art of Scientific Computing (Cambridge: Cambridge Univ. Press)

Rampf, C., & Wong, Y. Y. Y. 2012, JCAP, 06, 018

Rizzato, M., Benabed, K., Bernardeau, F., & Lacasa, F. 2019, MNRAS, 490, 4688

Sato, M., Hamana, T., Takahashi, R., et al. 2009, ApJ, 701, 945

Sato, M., & Nishimichi, T. 2013, PhRvD, 87, 123538

Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521

Schneider, A., Teyssier, R., Potter, D., et al. 2016, JCAP, 4, 047

Scoccimarro, R. 1997, ApJ, 487, 1

Scoccimarro, R. 2015, PhRvD, 92, 083532

Scoccimarro, R., Colombi, S., Fry, J. N., et al. 1998, ApJ, 496, 586

Scoccimarro, R., & Couchman, H. M. P. 2001, MNRAS, 325, 1312

Scoccimarro, R., Feldman, H. A., Fry, J. N., & Frieman, J. A. 2001, ApJ, 546, 652

Scoccimarro, R., & Frieman, J. A. 1999, ApJ, 520, 35

Sefusatti, E., Crocce, M., & Desjacques, V. 2010, MNRAS, 406, 1014

Sefusatti, E., Crocce, M., Pueblas, S., & Scoccimarro, R. 2006, PhRvD, 74, 023522

Sefusatti, E., Crocce, M., Scoccimarro, R., & Couchman, H. M. P. 2016, MNRAS, 460, 3624

Semboloni, E., Heymans, C., van Wærbeke, L., & Schneider, P. 2008, MNRAS, 388, 991

Semboloni, E., Hoekstra, H., & Schaye, J. 2013, MNRAS, 434, 148

Semboloni, E., Hoekstra, H., Schaye, J., van Daalen, M. P., & McCarthy, I. G. 2011a, MNRAS, 417, 2020

Semboloni, E., Schrabback, T., van Wærbeke, L., et al. 2011b, MNRAS, 410, 143

Shi, X., Joachimi, B., & Schneider, P. 2010, A&A, 523, A60

Shirasaki, M., Hamana, T., & Yoshida, N. 2015, MNRAS, 453, 3043

Simon, P., Semboloni, E., van Wærbeke, L., et al. 2015, MNRAS, 449, 1505

Slepian, Z., Eisenstein, D. J., Brownstein, J. R., et al. 2017, MNRAS, 469, 1738

Smith, R. E., & Angulo, R. E. 2019, MNRAS, 486, 1448

Smith, R. E., Peacock, J. A., Jenkins, A., et al. 2003, MNRAS, 341, 1311

Smith, R. E., Watts, P. I. R., & Sheth, R. K. 2006, MNRAS, 365, 214

Spergel, D. N., Bean, R., Doré, O., et al. 2007, ApJS, 170, 377

Springel, V. 2005, MNRAS, 364, 1105

Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676

Springel, V., Yoshida, N., & White, S. D. M. 2001, NewA, 6, 79

Sugiyama, N. S., Saito, S., Beutler, F., & Seo, H.-J. 2019, arXiv:1908.06234

Takada, M., & Hu, W. 2013, PhRvD, 87, 123504

Takada, M., & Jain, B. 2002, MNRAS, 337, 875

Takada, M., & Jain, B. 2004, MNRAS, 348, 897

Takahashi, R., Hamana, T., Shirasaki, M., et al. 2017, ApJ, 850, 24

Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M. 2012, ApJ, 761, 152

Troxel, M. A., & Ishak, M. 2012, MNRAS, 419, 1804

Troxel, M. A., & Ishak, M. 2015, PhR, 558, 1

Troxel, M. A., MacCrann, N., Zuntz, J., et al. 2018, PhRvD, 98, 043528

Valageas, P., & Nishimichi, T. 2011, A&A, 532, A4

Valageas, P., Sato, M., & Nishimichi, T. 2012, A&A, 541, A161

van Daalen, M. P., Schaye, J., Booth, C. M., & Dalla Vecchia, C. 2011, MNRAS, 415, 3649

van Engelen, A., Bhattacharya, S., Sehgal, N., et al. 2014, ApJ, 786, 13

van Uitert, E., Joachimi, B., Joudaki, S., et al. 2018, MNRAS, 476, 4662

Van Wærbeke, L., Benjamin, J., Erben, T., et al. 2013, MNRAS, 433, 3373

Vogelsberger, M., Genel, S., Springel, V., et al. 2014, MNRAS, 444, 1518

Yamamoto, K., Nan, Y., & Hikage, C. 2017, PhRvD, 95, 043528