# Flat to nonflat:
Calculating nonlinear power spectra of biased tracers for nonflat $\Lambda$CDM
model

Ryo Terasawa  0000-0002-1193-623X<EMAIL_ADDRESS>Kavli Institute for the
Physics and Mathematics of the Universe (WPI), The University of Tokyo
Institutes for Advanced Study (UTIAS), The University of Tokyo, Chiba
277-8583, Japan  Department of Physics, The University of Tokyo, Bunkyo, Tokyo
113-0031, Japan  Center for Data-Driven Discovery (CD3), Kavli IPMU (WPI),
UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan  Ryuichi
Takahashi  0000-0001-6021-0147 Faculty of Science and Technology, Hirosaki
University, 3 Bunkyo-cho, Hirosaki, Aomori 036-8561, Japan  Takahiro
Nishimichi  0000-0002-9664-0760 Department of Astrophysics and Atmospheric
Sciences, Faculty of Science, Kyoto Sangyo University, Motoyama, Kamigamo,
Kita-ku, Kyoto 603-8555, Japan  Center for Gravitational Physics and Quantum
Information, Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto
606-8502, Japan  Kavli Institute for the Physics and Mathematics of the
Universe (WPI), The University of Tokyo Institutes for Advanced Study (UTIAS),
The University of Tokyo, Chiba 277-8583, Japan  Masahiro Takada
0000-0002-5578-6472 Kavli Institute for the Physics and Mathematics of the
Universe (WPI), The University of Tokyo Institutes for Advanced Study (UTIAS),
The University of Tokyo, Chiba 277-8583, Japan  Center for Data-Driven
Discovery (CD3), Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa,
Chiba 277-8583, Japan

###### Abstract

The growth of large-scale structure, together with the geometrical information
of cosmic expansion history and cosmological distances, can be used to obtain
constraints on the spatial curvature of the universe that probes the early
universe physics, whereas modeling the nonlinear growth in a nonflat universe
is still challenging due to computational expense of simulations in a high-
dimensional cosmological parameter space. In this paper, we develop an
approximate method to compute the halo-matter and halo-auto power spectra for
nonflat $\Lambda$CDM model, from quantities representing the nonlinear
evolution of the corresponding flat $\Lambda$CDM model, based on the separate
universe (SU) method. By utilizing the fact that the growth response to long-
wavelength fluctuations (equivalently the curvature), $T_{\delta_{\rm b}}(k)$,
is approximated by the response to the Hubble parameter, $T_{h}(k)$, our
method allows one to estimate the nonlinear power spectra in a nonflat
universe efficiently from the power spectra of the flat universe. We use
$N$-body simulations to show that the estimator can provide the halo-matter
(halo-auto) power spectrum at $\sim 1\%$ ($\sim 2\%$ ) accuracy up to $k\simeq
3(1)\,h{\rm Mpc}^{-1}$ even for a model with large curvature $\Omega_{K}=\pm
0.1$. Using the estimator we can extend the prediction of the existing
emulators such as Dark Emulator to nonflat models without degrading their
accuracy. Since the response to long-wavelength fluctuations is also a key
quantity for estimating the super sample covariance (SSC), we discuss that the
approximate identity $T_{\delta_{\rm b}}(k)\approx T_{h}(k)$ can be used to
calculate the SSC terms analytically.

##  I Introduction

The spatial curvature of the universe ($\Omega_{K}$) is an important quantity
that characterizes the geometry of the universe and probes the physics of the
early universe including the mechanism of inflation. The curvature of the
universe influences various observables, including the spatial distribution of
galaxies, through both its geometric properties and its impact on the growth
of large-scale structures (LSS). The geometrical constraint, inferred from the
primary CMB anisotropy information of the Planck data (Planck Collaboration
_et al._ , 2020), is given as $\Omega_{K}=-0.044^{+0.018}_{-0.015}$ (68% CL,
Planck TT, TE, EE+lowE), implying a $2\sigma$ hint of the close geometry,
although most of the constraints are consistent with flat universe (e.g. Alam
_et al._ , 2017; Kumar  _et al._ , 2021; Vagnozzi  _et al._ , 2021a; Dhawan
_et al._ , 2021; Alam  _et al._ , 2021; Vagnozzi  _et al._ , 2021b; Liu  _et
al._ , 2022; DES Collaboration  _et al._ , 2022; Tristram  _et al._ , 2023).
This indicates the importance of constraining the curvature from LSS as an
independent probe. Furthermore, to obtain more precise constraints on the
curvature close to the amplitudes of primordial fluctuations,
$|\Omega_{K}|\sim 10^{-4}$, it is crucial to combine observations of these two
effects (e.g. Leonard  _et al._ , 2016). However, achieving this requires
development of a theoretical model that adequately incorporates mode coupling
effects. Mode coupling refers to the phenomenon where the non-linear nature of
gravity causes fluctuations of different wavelengths to interact with each
other. Constructing such a model is generally challenging. Despite its
inherent difficulty, incorporating small-scale information is vital for
obtaining stronger constraints on cosmological parameters. For flat
cosmological models in particular, fitting formulae or emulators based on the
predictions of $N$-body simulations have been developed (e.g. Peacock and
Dodds, 1996; Smith  _et al._ , 2003; Heitmann  _et al._ , 2010, 2009;
Takahashi  _et al._ , 2012; Heitmann  _et al._ , 2016; Lawrence  _et al._ ,
2017; Nishimichi  _et al._ , 2019; Reimberg  _et al._ , 2020; Mead  _et al._ ,
2021; Smith and Angulo, 2019; Euclid Collaboration  _et al._ , 2021). For
example, Dark Emulator (Nishimichi  _et al._ , 2019) predicts the basic
statistical quantities of dark matter halos (halos in what follows for
simplicity) such as their abundance as a function of mass and the halo power
spectrum down to the nonlinear scale. With these nonlinear predictions of the
halo statistics, together with the halo occupation distribution (hereafter,
HOD) description, Ref. Miyatake  _et al._ (2023) analyzed the two-point
correlation function of galaxy clustering, galaxy-galaxy lensing and cosmic
shear (so-called $3\times 2$ pt analysis), and obtained tighter constraints on
the cosmological parameters compared to the traditional analysis based on the
linear bias model (Sugiyama  _et al._ , 2023).

Theoretical models for the power spectrum in a nonflat universe are still in
the development stage compared to models in flat universe due to the
computational expense of covering a wider range in a multi-dimensional
cosmological parameter space. Because of this limitation of knowledge on
nonlinear scales, in the $3\times 2$ pt analysis of the DES Y3 (DES
Collaboration  _et al._ , 2022) data they used only the linear scales. In this
work, we provide a method to model the power spectrum down to nonlinear scales
so that one can gain information on the curvature from small scales.

The effect of curvature on structure formation can be identified with the very
long-wavelength density fluctuations in the $\Lambda$CDM model based on the
separate universe (SU) approach Baldauf  _et al._ (2011); Takada and Hu
(2013); Li  _et al._ (2014a, b); Wagner  _et al._ (2015a). In Terasawa  _et
al._ (2022), we developed a method to compute the nonlinear matter power
spectrum for nonflat cosmologies utilizing the SU approach. We utilized the
approximate identity that states the effect of these long-wavelength density
fluctuations on structure formation is well reproduced by the response to the
Hubble parameter $h$ (Schaan  _et al._ , 2014; Terasawa  _et al._ , 2022).
Using the response to long-wavelength fluctuations modeled in this way, we
have shown that the model of the matter power spectrum for flat
universes,whose accurate model calculation is already available from the
fitting formula or emulator, can be extended to calculate the nonlinear matter
power spectrum for a universe with non-zero curvature.

In this work, we apply the SU approach to predict the halo-matter and halo-
auto power spectra down to the nonlinear scales for nonflat cosmologies. We
will verify that the approximate identity holds for the halo-matter and halo-
auto power spectrum responses and construct the estimators of these power
spectra using the response to the Hubble parameter. The nonlinear prediction
of these spectra together with the HOD prescription provides a way to predict
the galaxy-galaxy lensing or galaxy clustering data down to small scales (e.g.
Miyatake  _et al._ , 2022a, 2023).

Another application of the modeling of the response to long-wavelength
fluctuations is to compute the super sample covariance (SSC) (Takada and Hu,
2013), which is the sample variance contribution caused by mode coupling with
long-wavelength fluctuations. SSC for the power spectrum can be computed using
its response to the super survey modes (e.g. Schaan  _et al._ , 2014). We also
provide a way to compute the response to the super survey modes using the
power spectrum response to $h$, which enables us to utilize a fitting formula
or an emulator validated only for flat geometry.

In the literature, there are analytical methods to compute the response, with
the perturbation theory (e.g. Baldauf  _et al._ , 2016) or halo model (e.g.
Krause and Eifler, 2017). Using these analytical methods, we can calculate the
response quickly. However, these models suffer from the limitation of scales
or inaccuracy. On the other hand, measuring the SSC from the scatter among an
ensemble of simulations (e.g. Bayer  _et al._ , 2023) is expected to provide
an accurate covariance matrix with a relatively high computational cost. The
size of the data vector is expected to become larger for future surveys,
especially the ones using multi-tracer, for which the number of simulations
needed could be $\mathcal{O}(10^{2})-\mathcal{O}(10^{3})$. In this paper, our
model is compared with simulations, and shown to be more accurate than the
analytical methods used in the literature, while the computational cost is
much cheaper than a simulation ensemble. Our method for computing the total
response would be implemented into the Core Cosmology Library (Chisari  _et
al._ , 2019).

This paper is organized as follows. In Sec. II we first review the SU approach
for the nonlinear matter power spectrum in Ref. Terasawa  _et al._ (2022) and
generalize it to predict the halo-matter or halo power spectrum. In Sec. III
we introduce our method to compute the total response to the super survey
modes using the response to $h$. In Sec. IV we describe details of $N$-body
simulations for flat and nonflat $\Lambda$CDM models used in this paper. In
Sec. V we present the main results of this paper and show a numerical
validation of our methods. Sec. VI is devoted to discussion and conclusion. In
Appendix A we describe the details of how to compute the power spectrum
response to the long-wavelength fluctuations. Throughout the paper, we assume
$\Lambda$CDM model and flat geometry for the fiducial cosmology. We also
assume that the SSC is calculated in this fiducial cosmology.

##  II SU approach for $P(k;\Omega_{K})$

###  II.1 SU approach for $P_{\rm mm}(k;\Omega_{K})$

In this section we briefly review the SU approach to predict the nonlinear
matter power spectrum $P_{\rm mm}(k)$ for a nonflat universe following Ref.
Terasawa  _et al._ (2022).

First, since the effect of the curvature or the long-wavelength density
fluctuations on structure formation appears only in the late universe,
throughout this paper we consider a model where structure formation in the
early universe is identical to that of the fiducial flat universe.
Specifically, we keep the parameters

$\displaystyle\left\\{\omega_{\rm c},\omega_{\rm b},A_{s},n_{s}\right\\},$ (1)

fixed, where $\omega_{\rm c}(\equiv\Omega_{\rm c}h^{2})$ and $\omega_{\rm
b}(\equiv\Omega_{\rm b}h^{2})$ are the physical density parameters of CDM and
baryon, respectively, and $A_{s}$ and $n_{s}$ are the amplitude (at the pivot
scale $k_{\rm pivot}=0.05\,{\rm Mpc}^{-1}$) and the spectral tilt of the power
spectrum of primordial curvature perturbations. Note that we fix the sum of
neutrino masses so that the early universe physics remains unchanged, and
treat its impact only through the linear transfer function of total matter
fluctuations at $z=0$ (see Refs. Nishimichi  _et al._ , 2019; Bayer  _et al._
, 2021, for more details). These models with the four fixed parameters (Eq. 1)
share the same linear power spectrum $P_{L}(k)$ at sufficiently high redshifts
$z_{i}\gg 1$. The remaining parameters that affect structure formation are
$\Omega_{K}$ and $h$ within the nonflat $\Lambda$CDM model. Note that the
density parameter for the cosmological constant, $\Omega_{\Lambda}$, is
automatically determined once the parameters above are all fixed.

In the SU approach, the effects of the background density modulation
$\delta_{\rm b}(t)$ in a flat universe are interpreted as the local effective
cosmology with modified background density

$\displaystyle\bar{\rho}_{\rm m}(t)=\bar{\rho}_{{\rm
m}f}(t)\left[1+\delta_{\rm b}(t)\right],$ (2)

and non-zero curvature corresponding to $\delta_{\rm b}$ (see below). Hence it
gives a mapping between nonflat and flat $\Lambda$CDM models. Hereafter we
call the two models $\Omega_{K}$-$\Lambda$CDM and $f\Lambda$CDM, respectively,
and we denote quantities in the $f\Lambda$CDM model by sub/superscript “$f$”.
We assume $\delta_{\rm b}(t)$ evolves according to the linear growth factor
$D_{f}(t)$ as $\delta_{\rm b}(t)\propto D_{f}(t)$.

The curvature and overdensity are related to each other via

$\displaystyle\Omega_{K}=-\frac{5\Omega_{\mathrm{m}}}{3}\frac{{\delta_{\rm
b}(t)}}{D_{f}(t)}.$ (3)

The remaining parameter, $h$, is mapped as

$\displaystyle h_{f}=h(1-\Omega_{K})^{1/2}.$ (4)

Note that the redshifts in these two cosmologies at a given cosmic time $t$
are related as

$\displaystyle(1+z_{f})[1+\delta_{\rm b}(z_{f})]^{1/3}=1+z.$ (5)

Inversely, we can treat the nonflat universe as an overdense/underdense region
in the corresponding flat universe, and the power spectrum in the target
nonflat universe can be approximated by Taylor expansion around the flat
universe as discussed in Ref. Terasawa  _et al._ (2022):

$\displaystyle{P}(k,z;\Omega_{K})$ $\displaystyle\simeq
P^{f}(k,z_{f};\delta_{\rm b})$
$\displaystyle\simeq\left.P^{f}(k,z_{f})\right|_{\delta_{\rm
b}=0}\left[1+\left.\frac{\partial\ln P^{f}(k,z_{f};\delta_{\rm
b})}{\partial\delta_{\rm b}}\right|_{{\rm G},\delta_{\rm b}=0}\delta_{\rm
b}\right]$ $\displaystyle\equiv\tilde{P}(k,z;\Omega_{K}),$ (6)

where $\delta_{\rm b}\equiv\delta_{\rm b}(z_{f})$. In the last equality on the
r.h.s., we have put the tilde symbol $\tilde{\hskip 10.00002pt}$ on
$\tilde{P}(k,z;\Omega_{K})$ to explicitly denote that $\tilde{P}$ is an
“estimator” of the nonlinear matter power spectrum for the
$\Omega_{K}$-$\Lambda$CDM model. We introduced subscript “G” to $\partial
P^{f}(k)/\partial\delta_{\rm b}$ to emphasize that it is the growth response
(Li  _et al._ , 2014a). The growth response can be computed as the difference
between the power spectra at a fixed comoving wavenumber $k$ in the two SU
cosmologies, omitting the dilation effect that originates from the difference
in the scale factors (Sec. III). Throughout the paper, the wavenumber $k$
refers to a given comoving wavenumber even if the corresponding physical
scales differ among the cosmologies. Note that the expansion of Eq. (6) is
applicable to not only the matter power spectrum but also the halo-matter and
halo-auto power spectra as we will show below.

Furthermore, we define a normalized growth response $T_{\delta_{\rm b}}^{\rm
mm}(k)$ as

$\displaystyle T_{\delta_{\rm b}}^{\rm mm}(k)\equiv\left[2\frac{\partial\ln
D(\delta_{\rm b})}{\partial\delta_{\rm b}}\right]^{-1}\left.\frac{\partial\ln
P_{\rm mm}(k;\delta_{\rm b})}{\partial\delta_{\rm b}}\right|_{{\rm
G},\delta_{\rm b}=0}.$ (7)

The normalized response has an asymptotic behavior of $T_{\delta_{\rm b}}^{\rm
mm}\rightarrow 1$ at the linear limit $k\rightarrow 0$, because $P_{\rm
mm}(k,z)\propto D(z)^{2}P_{L}(k,z_{i})$ in such linear regime. The linear
limit of the matter power spectrum growth response is given in Refs. Baldauf
_et al._ (2011); Takada and Hu (2013) as

$\displaystyle 2\frac{\partial\ln D(\delta_{\rm b})}{\partial\delta_{\rm
b}}\approx\frac{26}{21}.$ (8)

Since we have a prediction for the linear power spectrum for the
$\Omega_{K}$-$\Lambda$CDM model, we only need to expand the nonlinear
correction $B(k,z)$ defined as $P(k,z)\equiv P_{L}(k,z)B(k,z)$. Using the
normalized response $T_{\delta_{\rm b}}(k)$, we can write the expansion as

$\displaystyle\tilde{P}(k;\Omega_{K})=\left(\frac{D(z)}{D_{f}(z_{f})}\right)^{2}P^{f}(k)\left[1+\frac{26}{21}(T_{\delta_{\rm
b}}(k)-1)\delta_{\rm
b}(z_{f})+\left(\frac{G_{2}(k)}{2}-\frac{1501}{1323}-\left(\frac{26}{21}\right)^{2}(T_{\delta_{\rm
b}}(k)-1)\right)\delta_{\rm b}(z_{f})^{2}\right],$ (9)

where

$\displaystyle
G_{2}(k)\equiv\frac{1}{P^{f}(k)}\left.\frac{\partial^{2}P(k)}{\partial\delta_{\rm
b}^{2}}\right|_{{\rm G},\delta_{\rm b}=0}$ (10)

and $1501/1323$ is the linear limit of $G_{2}(k)/2$ for the matter power
spectrum (Wagner  _et al._ , 2015b), derived in an Einstein-de Sitter
cosmology. We explicitly wrote down the expansion up to the second order in
$\delta_{\rm b}$ and show below that the second order term is negligible for
moderate values of $|\Omega_{K}|$ that are consistent with the current bounds,
$|\Omega_{K}|\lesssim\mathcal{O}(10^{-1})$. For the matter power spectrum, the
normalized growth response $T_{\delta_{\rm b}}(k)$ ranges about $[0.3,1.6]$ at
redshift $z\simeq[0,1.5]$ in the range of $k\simeq[10^{-2},6]\,h{\rm
Mpc}^{-1}$ (see Fig. 2 of Ref. (Terasawa  _et al._ , 2022)) and $G_{2}(k)/2$
ranges about $[0.5,2]$ at $z=0$ in the range of $k\simeq[10^{-2},2]\,h{\rm
Mpc}^{-1}$ (Fig. 2 of (Wagner  _et al._ , 2015a)). Since the prefactors of the
first and second order terms before $\delta_{\rm b}$ and $\delta_{\rm b}^{2}$
are at most $\mathcal{O}(10^{-1})$, the contributions from the first and
second order terms are at percent and sub-percent level, respectively, even
for a large curvature case of $|\Omega_{K}|\sim\mathcal{O}(10^{-1})$
(corresponding to $\delta_{\rm b}(t)\sim\mathcal{O}(10^{-1})$). As we will
show below, the second or higher order terms are also negligible for the halo-
matter or halo-auto power spectrum. Hence we consider the expansion up to the
first order of $\delta_{\rm b}$ throughout this paper.

We also define the normalized growth response to $h$ within the flat model as:

$\displaystyle T_{h}^{\rm mm}(k)\equiv\left[2\frac{\partial\ln D}{\partial
h}\right]^{-1}\frac{\partial\ln P_{\rm mm}(k;h)}{\partial h}.$ (11)

As shown numerically in Refs. (Li  _et al._ , 2014b; Terasawa  _et al._ ,
2022) (also as will be shown in Fig. 1), these two responses agree well even
in the nonlinear regime: $T_{\delta_{\rm b}}^{\rm mm}(k)\simeq T_{h}^{\rm
mm}(k)$. Finally, using this approximate identity $T_{\delta_{\rm b}}(k)\simeq
T_{h}(k)$ and ignoring the terms of second or higher order in $\delta_{\rm
b}(t)$, we can approximate the estimator as

$\displaystyle\tilde{P}(k,z;\Omega_{K})$
$\displaystyle=\left(\frac{D(z)}{D_{f}(z_{f})}\right)^{2}P^{f}(k,z_{f})$
$\displaystyle\times\left[1+\frac{26}{21}(T_{h}(k)-1)\delta_{\rm
b}(z_{f})\right].$ (12)

By construction, the estimator reproduces the linear prediction at
$k\rightarrow 0$, where the first order term vanishes because $T_{h}^{\rm
mm}(k\rightarrow 0)=1$. The novel feature of the above estimator (Eq. 12) is
that it allows one to compute the nonlinear matter power spectrum for nonflat
universe from the quantities in flat universe. With this estimator we can
extend emulators available in the community, which are applicable only for
flat universes, to predict the power spectrum in a nonflat universe.

###  II.2 SU approach for halo power spectra

Figure 1: The responses of the power spectra to $\delta_{\rm b}$ and $h$
measured from the $N$-body simulations (see Section IV for details of the
simulations). The circle, triangle, and square symbols show the responses of
the matter power spectrum ($T^{\rm mm}(k)$), halo-matter power spectrum
($T^{\rm hm}(k)$), and halo power spectrum ($T^{\rm hh}(k)$), respectively.
For $T^{\rm hh}(k)$, we only plot the scales where halo power spectra after
subtracting the shot noise have positive values. We slightly shift the symbols
of $T_{h}(k)$ along the $x$-axis for illustration.

In this paper, we extend the SU approach mentioned above to the halo-matter
and halo-auto power spectra. Similarly to the above discussion, we define
$T_{\delta_{\rm b}}^{\rm XY}(k)$ and $T_{h}^{\rm XY}(k)$ as

$\displaystyle T_{\delta_{\rm b}}^{\rm XY}(k)$
$\displaystyle\equiv\left[2\frac{\partial\ln D(\delta_{\rm
b})}{\partial\delta_{\rm b}}\right]^{-1}\left.\frac{\partial\ln P_{\rm
XY}(k;\delta_{\rm b})}{\partial\delta_{\rm b}}\right|_{{\rm G},\delta_{\rm
b}=0},$ $\displaystyle T_{h}^{\rm XY}(k)$
$\displaystyle\equiv\left[2\frac{\partial\ln D}{\partial
h}\right]^{-1}\frac{\partial\ln P_{\rm XY}(k;h)}{\partial h},$ (13)

where $\rm{XY}=\\{\rm mm,hm,hh\\}$. When considering the response of $P_{\rm
hm}(k)$ or $P_{\rm hh}(k)$, we perform the derivative keeping the comoving
halo number density fixed. As we will specify the latter in Sec. A.3, these
responses are not the same as the growth response in the SU approach where the
derivative is performed keeping the halo mass threshold fixed.

As a highlight of our approach, in Fig. 1 we show the responses of power
spectra of matter-matter, halo-matter, and halo-halo to $\delta_{\rm b}$ or
$h$, $T^{\rm mm}$, $T^{\rm hm}$ and $T^{\rm hh}$, which are computed using Eq.
(13). At the large-scale limit $k\rightarrow 0$, the matter-matter response
$T^{\rm mm}_{\delta_{\rm b}}=1$ and $T^{\rm mm}_{h}=1$ by definition. Given
that the abundance-matched halo samples correspond to the same initial density
peaks, which is true as long as mergers do not severely affect this
correspondence, the clustering amplitudes of the abundance-matched halos on
large scales barely change. This is why $T^{\rm hh}(k\rightarrow 0)\simeq 0$
and $T^{\rm hm}$ is smaller than $T^{\rm mm}$. The figure also shows that the
responses for “mm”, “hm” and “hh” have quite different $k$-dependence. Hence
accurate calibration is important to capture the scale dependence. The
responses for $\delta_{\rm b}$ and $h$ are in good agreement with each other,
for all the “mm”, “hm” and “hh” power spectra, which validates our method to
approximate the response to $\delta_{\rm b}$ by the response to $h$.

Finally, using the approximate identity $T_{\delta_{\rm b}}(k)\simeq
T_{h}(k)$, we can obtain the estimator as

$\displaystyle\tilde{P}_{\rm XY}(k,z;\Omega_{K})$
$\displaystyle=\left(\frac{D(z)}{D_{f}(z_{f})}\right)^{2}P_{\rm
XY}^{f}(k,z_{f})$ $\displaystyle\times\left[1+\frac{26}{21}(T^{\rm
XY}_{h}(k)-1)\delta_{\rm b}(z_{f})\right].$ (14)

##  III Total response estimator using flat $\Lambda$CDM model

In this section, we discuss that our approximation using the response of the
power spectrum to the Hubble parameter can be used to calibrate the SSC of a
cosmological observable. Here the SSC is the sampling variance error of the
observable in a finite-volume survey, which arises from the mode coupling of
density fluctuations in the survey with density fluctuations on scales greater
than the survey window, i.e. super-survey modes (Takada and Hu, 2013; Li  _et
al._ , 2014a; Takahashi  _et al._ , 2019). For instance, the SSC gives a
dominant source of the sample variance on scales larger than the scales where
the shot noise dominates, e.g., for the two-point correlation function of
cosmic shear, which is given by the weighted line-of-sight integration of the
matter power spectrum.

SSC for the X- and Y-observables in the $i$-th and $j$-th bins, respectively,
is generally expressed, as proposed by Ref. Takada and Hu (2013), as

$\displaystyle C^{({\rm XY}){\rm SSC}}_{ij}\approx\sigma_{\rm
b}^{2}\frac{\partial\mathcal{O}_{{\rm X}i}}{\partial\delta_{\rm
b}}\frac{\partial\mathcal{O}_{{\rm Y}j}}{\partial\delta_{\rm b}},$ (15)

where $\sigma_{\rm b}^{2}$ is the variance in the linear density filtered by a
window function assuming a sufficiently large survey volume in that the super-
survey modes ($\delta_{\rm b}$) are in the linear regime, and $\partial{\cal
O}/\partial\delta_{\rm b}$ is the total response of the observable to
$\delta_{\rm b}$ including both the growth and dilation responses (see below).
The observables can be any statistical cosmological quantities such as power
spectrum, bispectrum, and cluster mass function. In this paper, we consider,
as $\mathcal{O}_{\rm X}$ and/or $\mathcal{O}_{\rm Y}$, either of the matter,
halo-matter, or halo-halo power spectum and the $i$-th index in the above
equation corresponds to the $i$-th $k$ bin. Otherwise, the above equation is a
general expression of SSC. We also note that, once the total response for the
3D observable (e.g. the matter power spectrum) is given, the SSC to the
corresponding angular observable (e.g., cosmic power spectrum) can be obtained
by a weighted line-of-sight integral of the 3D SSC term (Takahashi  _et al._ ,
2019).

In the following, we provide an approximate way to compute the total response
using the power spectrum response to $h$.

###  III.1 Preliminary

First, we summarize analytical methods of computing the total response
proposed in the literature.

$\bullet$Perturbation theory – The total response can be computed using the
perturbation theory (PT) Baldauf  _et al._ (2016) as

$\displaystyle\frac{\partial\ln P_{\rm mm,PT}(k)}{\partial\delta_{\rm b}}$
$\displaystyle=\frac{68}{21}-\frac{1}{3}\frac{\partial\ln
k^{3}P_{L}(k)}{\partial\ln k},$ $\displaystyle\frac{\partial\ln P_{\rm
gm,PT}(k)}{\partial\delta_{\rm b}}$
$\displaystyle=\frac{68}{21}+\frac{b_{g,2}}{b_{g}}-b_{g}-\frac{1}{3}\frac{\partial\ln
k^{3}P_{L}(k)}{\partial\ln k},$ $\displaystyle\frac{\partial\ln P_{\rm
gg,PT}(k)}{\partial\delta_{\rm b}}$
$\displaystyle=\frac{68}{21}+2\frac{b_{g,2}}{b_{g}}-2b_{g}-\frac{1}{3}\frac{\partial\ln
k^{3}P_{L}(k)}{\partial\ln k},$ (16)

where $b_{g}$ and $b_{g,2}$ denote the first- and second-order Eulerian galaxy
biases. Using the halo occupation distribution (HOD) prescription (e.g. Zheng
_et al._ , 2005), $b_{g}$ is given as

$\displaystyle
b_{g}=\frac{1}{\bar{n}_{g}}\int~{}\mathrm{d}M\frac{\mathrm{d}n}{\mathrm{d}M}(M)b_{{\rm
h},1}(M)\langle N_{\mathrm{g}}(M)\rangle,$ (17)

with the mean number density of galaxies

$\displaystyle\bar{n}_{g}=\int~{}\mathrm{d}M\frac{\mathrm{d}n}{\mathrm{d}M}(M)\langle
N_{\mathrm{g}}(M)\rangle,$ (18)

where $\langle N_{\mathrm{g}}(M)\rangle$ is the mean number of galaxies in a
halo, $\mathrm{d}n/\mathrm{d}M$ is the halo mass function in the mass range
$[M,M+\mathrm{d}M]$, and $b_{{\rm h},1}(M)$ is the linear bias of halos of
mass $M$. Throughout the paper, we use the fitting function by Ref. Tinker
_et al._ (2010) for the linear halo bias $b_{{\rm h},1}(M)$.

Similarly, $b_{g,2}$ is given as

$\displaystyle
b_{g,2}=\frac{1}{\bar{n}_{g}}\int\\!\mathrm{d}M\frac{\mathrm{d}n}{\mathrm{d}M}(M)b_{{\rm
h},2}(M)\langle N_{\mathrm{g}}(M)\rangle,$ (19)

where $b_{{\rm h},2}(M)$ is the 2nd-order halo bias of halos of mass $M$.
Throughout this paper, we use the fitting formula of $b_{{\rm h},2}$ proposed
in Ref. Hoffmann  _et al._ (2017), which gives $b_{{\rm h},2}(M)$ in terms of
the linear halo bias, i.e. $b_{{\rm h},2}=b_{{\rm h},2}(b_{{\rm h},1})$.

As we will see in Figs. 6 and 7, the perturbation theory well describes the
response on large scale, but it fails to reproduce the small-scale behavior,
where perturbation theory itself starts to break down.

$\bullet$Halo Model – The halo model (e.g. Sheth and Tormen, 2002) is commonly
used to calculate the response of the matter power spectrum in the literature
(Takada and Hu, 2013; Li  _et al._ , 2014a; Chiang  _et al._ , 2014). Ref.
(Krause and Eifler, 2017) derived the responses for the multiprobe power
spectrum, denoted as $P_{\rm XY}$, based on the halo model. Following the
notations in Ref. (Krause and Eifler, 2017), the responses are written as

$\displaystyle\frac{\partial P_{\rm XY}(k)}{\partial\delta_{\rm b}}$
$\displaystyle=\left(\frac{47}{21}+\frac{I_{\rm X}^{2}(k)}{I_{\rm
X}^{1}(k)}+\frac{I_{\rm Y}^{2}(k)}{I_{\rm
Y}^{1}(k)}-\frac{1}{3}\frac{\partial\ln P_{L}(k)}{\partial\ln k}\right)$
$\displaystyle\times I_{\rm X}^{1}(k)I_{\rm Y}^{1}(k)P_{L}(k)$
$\displaystyle+I_{\rm XY}^{1}(k,k)-\left[b_{\rm X,X=g}+b_{\rm
Y,Y=g}\right]P_{\rm XY}(k),$ (20)

where $b_{{\rm X,X}={\rm g}}=b_{\rm g}$ for ${\rm X}={\rm g}$ and otherwise
$b_{{\rm X,X}={\rm g}}=0$ and the functions $I_{\rm X}^{\alpha}$ and $I_{\rm
XY}^{\alpha}$ are defined as

$\displaystyle I_{\rm X}^{\alpha}(k)$
$\displaystyle=\int\\!\mathrm{d}M\frac{\mathrm{d}n}{\mathrm{d}M}(M)b_{\rm{h},\alpha}(M)\tilde{u}_{\rm
X}(k;M),$ $\displaystyle I_{\rm XY}^{\alpha}(k,k^{\prime})$
$\displaystyle=\int\\!\mathrm{d}M\frac{\mathrm{d}n}{\mathrm{d}M}(M)b_{\rm{h},\alpha}(M)\tilde{u}_{\rm
X}(k;M)\tilde{u}_{\rm Y}(k^{\prime};M),$ (21)

where $\tilde{u}_{\rm X}(k;M)$ is the Fourier transform of the radial profile
of tracers X in host halos, multiplied by the number density normalization and
the HOD function for galaxies following the notations in Ref. Krause and
Eifler (2017). Note that $\tilde{u}_{\rm X}(k;M)$ has a dimension of volume.

We stress that we included the terms proportional to $I_{\rm X}^{2}(k)$ in the
parentheses of Eq. (20), which were missing in Ref. (Krause and Eifler, 2017).
These terms depend on the 2nd-order halo bias, $b_{{\rm h},2}$, which arises
from the response of the linear halo bias $b_{{\rm h},1}$ to $\delta_{\rm b}$.
For the matter power spectrum (for the case ${\rm X}={\rm Y}={\rm m}$), the
terms proportional to $I_{\rm m}^{2}(k)$ are negligible due to the halo model
consistency relation (Wagner  _et al._ , 2015b). Since
$I_{g}^{1}(k)\rightarrow b_{g}$ and $I_{g}^{2}(k)\rightarrow b_{g,2}$ at the
$k\rightarrow 0$ limit, $I_{g}^{2}(k)/I_{g}^{1}(k)\rightarrow b_{g,2}/b_{g}$
at the limit, which reproduces the PT response (Eq. 16). Thus we believe that
the above response formula (Eq. 20) is more accurate in the sense that it
includes the PT theory at the limit of $k\rightarrow 0$.

We will below assess the accuracy and limitation of these analytical formulae
of the power spectrum responses by comparing the model predictions with the
simulation results.

###  III.2 $h$ response method

The long-wavelength modes (super-survey modes) whose wavelengths are larger
than the survey volume/simulation box affect the growth of LSS via mode
coupling. The effects of the long-wavelength modes can be considered as the
background density (mean density) modulation, $\delta_{\rm b}$, which in turn
can be interpreted as the local effective cosmology with non-zero curvature
corresponding to $\delta_{\rm b}$ (the SU approach (e.g. Baldauf  _et al._ ,
2011; Takada and Hu, 2013; Li  _et al._ , 2014a, b; Wagner  _et al._ ,
2015a)).

Different large-scale structure tracers are measured with respect to either
the global or local mean density of the tracers, where the local mean is the
average density of the tracers in a finite-volume survey. Depending on this
difference, the response of power spectrum of X and Y tracers, $P_{\rm XY}$,
to the super survey modes can be decomposed into three contributions (e.g. Li
_et al._ , 2014a; Takahashi  _et al._ , 2019):

$\displaystyle\left.\frac{\partial\ln P_{\rm XY}(k;\delta_{\rm
b})}{\partial\delta_{\rm b}}\right|_{\rm total}$
$\displaystyle=n+\left.\frac{\partial\ln P_{\rm XY}(k;\delta_{\rm
b})}{\partial\delta_{\rm b}}\right|_{\rm G}$
$\displaystyle-\frac{1}{3}\frac{\partial\ln k^{3}P_{\rm XY}(k)}{\partial\ln
k}.$ (22)

The first term on the r.h.s. accounts for the change in the mean density of
the tracers by $\delta_{\rm b}$, used in the definition of the density
fluctuation field; for example, for the matter field (XY=mm),
$\delta_{\mathrm{m}}=\rho_{\mathrm{m}}/\bar{\rho}_{\mathrm{m}}-1$. The cases
of $n=2,1$ and $0$ correspond to “${\rm XY}$”$=$“${\rm mm}$”, “${\rm gm}$” and
“gg”, respectively, which are relevant to cosmic shear, galaxy-galaxy weak
lensing and galaxy-galaxy clustering, respectively (Takahashi  _et al._ ,
2019). The second term is the growth response, which describes the fractional
change in the power spectrum amplitude by the presence of $\delta_{\rm b}$.
The last term is the dilation response, which originates from the change in
the physical scale corresponding to a given comoving scale due to the change
in cosmic expansion.

Using this decomposition, we propose the matter power spectrum response
estimator as

$\displaystyle\left.\frac{\partial\ln P_{\rm mm}(k;\delta_{\rm
b})}{\partial\delta_{\rm b}}\right|_{\rm total}$
$\displaystyle=2+\frac{26}{21}T_{h}^{\rm mm}(k)-\frac{1}{3}\frac{\partial\ln
k^{3}P_{\rm mm}(k)}{\partial\ln k},$

where we have used the fact that the normalized growth response to
$\delta_{\rm b}$ is approximated by that to $h$:

$\displaystyle\left.\frac{\partial\ln P_{\rm mm}(k;\delta_{\rm
b})}{\partial\delta_{\rm b}}\right|_{G}=\frac{26}{21}T_{\delta_{\rm b}}^{\rm
mm}(k)\simeq\frac{26}{21}T_{h}^{\rm mm}(k).$ (24)

Similarly, as we will show in Appendix A.3, we can use the approximate
identity of the responses of halo-matter and halo-auto power spectra to
predict the growth responses of galaxy-matter and galaxy-auto power spectra
with the HOD description:

$\displaystyle\left.\frac{\partial P_{\rm gm}(k)}{\partial\delta_{\rm
b}}\right|_{\rm G}$
$\displaystyle\simeq\frac{1}{\bar{n}_{g}}\int\mathrm{d}M\left\\{\dfrac{\mathrm{d}\langle\tilde{U}_{g}(k;M)\rangle}{\mathrm{d}M}n(>M)\right.\left[b_{\mathrm{h},1}^{L}(>M)P_{\rm
hm}(k;M)+\frac{26}{21}T_{h}^{\rm hm}(k;>M)P_{\rm hm}(k;>M)\right]$
$\displaystyle~{}~{}~{}~{}~{}\left.+\frac{1}{3}P_{\mathrm{hm}}(k;>M)\frac{\partial}{\partial\ln
k}\frac{\mathrm{d}\ln\langle\tilde{U}_{g}(k;M)\rangle}{\mathrm{d}M}\right\\}-b_{g}^{L}P_{\rm
gm}(k),$ (25) $\displaystyle\left.\frac{\partial P_{\rm
gg}(k)}{\partial\delta_{\rm b}}\right|_{\rm G}$
$\displaystyle\simeq\frac{1}{\bar{n}_{g}^{2}}\iint\mathrm{d}M\mathrm{d}M^{\prime}\left\\{\frac{\mathrm{d}\langle\tilde{U}_{g}(k;M)\rangle}{\mathrm{d}M}\frac{\mathrm{d}\tilde{U}_{g}(k;M^{\prime})}{\mathrm{d}M^{\prime}}\right.n(>M)n(>M^{\prime})\left[2b_{\mathrm{h},1}^{L}(>M^{\prime})P_{\rm
hh}(k;>M,M^{\prime})\right.$
$\displaystyle~{}~{}~{}~{}~{}+\left.\left.\left(\frac{26}{21}\right)T_{h}^{\rm
hh}(k;>M,>M^{\prime})P_{\rm hh}(k;>M,>M^{\prime})+\frac{2}{3}{P_{\rm
hh}(k;>M,>M^{\prime})}\frac{\partial}{\partial\ln
k}\frac{\mathrm{d}\ln\langle\tilde{U}_{g}(k;M)\rangle}{\mathrm{d}M}\right]\right\\}$
$\displaystyle~{}~{}~{}~{}~{}+\frac{1}{\bar{n}_{g}^{2}}\int\mathrm{d}M\frac{\mathrm{d}n}{\mathrm{d}M}(M)\left[b_{\mathrm{h},1}^{L}(M)\langle\tilde{U}_{g}^{2}(k;M)\rangle+\frac{1}{3}\frac{\partial\langle\tilde{U}_{g}^{2}(k;M)\rangle}{\partial\ln
k}\right]-2b_{g}^{L}P_{\rm gg}(k),$ (26)

where the super-script $L$ denotes the Lagrangian bias and
$\tilde{U}_{g}(k;M)$ describes the quantity arising from the galaxy-halo
connection (see Appendix A for the definition).

Since the $h$ response is defined within flat cosmologies, we can predict the
power spectrum response from flat cosmology predictions or simulations without
performing a pair of nonflat universe simulations following the SU framework.
As we will show below, using the simulation-based emulator we can predict the
response accurately down to the nonlinear scale.

Finally, to ensure the estimator is correct at large scales, we stitch the
estimator with perturbation theory prediction (Baldauf  _et al._ , 2016)
described above in Eq. (16). Specifically, we smoothly stitch the estimator
with perturbation theory prediction as

$\displaystyle\frac{\partial\ln P_{\rm XY}(k)}{\partial\delta_{\rm b}}$
$\displaystyle=\frac{\partial\ln P_{\rm XY,PT}(k)}{\partial\delta_{\rm
b}}e^{-(k/k_{\rm switch})}$ $\displaystyle+\frac{\partial\ln P_{\rm
XY,est}(k)}{\partial\delta_{\rm b}}\left[1-e^{-(k/k_{\rm switch)}}\right].$
(27)

Throughout the paper, we adopt the switching scale $k_{\rm
switch}=0.08~{}h~{}{\rm Mpc}^{-1}$ for XY = mm, gm and gg.

##  IV Simulations

Table 1: Details of $N$-body simulations for different cosmological models. The columns “$\Omega_{K}$” and “$h$” give their values of the curvature parameter and Hubble parameter that are employed in each simulation, while we fix other cosmological parameters $\\{\omega_{\rm c},\omega_{\rm b},A_{s},n_{s}\\}$, which are needed to specify the linear power spectrum for the initial conditions, to the values for the fiducial Planck cosmology (see text for details). $\Omega_{\rm m}$ and $\Omega_{\Lambda}$ are specified by a given set of $\Omega_{K}$ and $h$, because we keep $\Omega_{\rm m}h^{2}$ fixed and $\Omega_{\Lambda}=1-\Omega_{\rm m}-\Omega_{\rm K}$. The column “$N_{\rm real}$” denotes the number of realizations, with different initial seeds, used for each model. We employ the “paired-and-fixed” method in Angulo and Pontzen (2016) to reduce the sample variance effect in small $k$ bins for the power spectrum measurement: it uses the paired (2) simulations by design (see text for details). The column “redshift ($z$)” gives the redshifts of simulation outputs: for $\Omega_{K}$-$\Lambda$CDM model, we properly choose the redshifts corresponding to the same cosmic time for each of redshifts, $z=\\{0.0,0.549,1.025,1.476\\}$ in the “fiducial” model in the SU approach (see around Eq. 5 in Section II). All the simulations are done in the fixed comoving box size without $h$ in its units, i.e. $L\simeq 2.97\,{\rm Gpc}$ (corresponding to $2\,h_{f}^{-1}{\rm Gpc}$ for the fiducial model) and with the same particle number, i.e. $N_{p}=2048^{3}$. Name | $\Omega_{K}$ | $h$ | $N_{\rm real}$ | redshift ($z$)
---|---|---|---|---
flat (fiducial) | 0 | $0.6727$ | 2 | $\\{0.0,0.549,1.025,1.476\\}$
$\Omega_{K}$-$\Lambda$CDM1 | 0.00663 | 0.6749 | 8 | $\\{-0.0033,0.544,1.018,1.467\\}$
| $-0.00672$ | 0.6705 | 8 | $\\{0.0033,0.554,1.031,1.484\\}$
$\Omega_{K}$-$\Lambda$CDM2 | 0.1 | 0.7091 | 2 | $\\{-0.059,0.482,0.955,1.405\\}$
| $-0.1$ | 0.6414 | 2 | $\\{0.043,0.600,1.079,1.531\\}$
$h$-$\Lambda$CDM | 0 | 0.6927 | 8 | $\\{0.0,0.549,1.025,1.476\\}$
| 0 | 0.6527 | 8 | $\\{0.0,0.549,1.025,1.476\\}$

###  IV.1 $N$-body simulations

In this section we give a brief summary of the simulations used in this paper.
Our simulations follow the method in Nishimichi  _et al._ (2019).

We use Gadget-2 (Springel, 2005) to carry out $N$-body simulation for a given
cosmological model. The initial conditions are set up at redshift $z_{i}=29$
using the second-order Lagrangian perturbation theory (Scoccimarro, 1998;
Crocce and Scoccimarro, 2006) implemented by Nishimichi  _et al._ (2009) and
then parallelized in Valageas and Nishimichi (2011). We use the public code
CAMB (Lewis  _et al._ , 2000) to compute the transfer function for a given
model, which is used to compute the input linear power spectrum. For all
simulations in this paper, we use the same simulation box size in Gpc (i.e.
without $h$ in the units) and the same number of particles:
$L=2\,h_{f}^{-1}\mathrm{Gpc}\simeq 2.97\,{\rm Gpc}$ (without $h$ in units) and
$N_{p}=2048^{3}$, which correspond to the particle Nyquist wavenumber,
$k=3.2~{}h_{f}\rm{Mpc}^{-1}$. In the following we will show the results at
wavenumbers smaller than this Nyquist wavenumber. In this paper we use
simulations for 4 different cosmological models, denoted as “fiducial” flat
$\Lambda$CDM, “$\Omega_{K}$-$\Lambda$CDM1”, “$\Omega_{K}$-$\Lambda$CDM2”, and
“$h$-$\Lambda$CDM” models, respectively, as given in Table 1. Here the
cosmological parameters for the “fiducial” model are chosen to be consistent
with those for the Planck 2015 best-fit cosmology Planck Collaboration  _et
al._ (2016). The cosmological parameters for each of the nonflat cosmological
models are chosen so that it has the fiducial $\Lambda$CDM model as the
corresponding flat $\Lambda$CDM model in the SU approach. We use paired
simulations for “$\Omega_{K}$-$\Lambda$CDM1” model to compute the power
spectrum response with respect to $\delta_{\rm b}$ ($T_{\delta_{\rm b}}$),
where the curvature parameters are specified by $\delta_{\rm b}=\pm 0.01$ at
$z_{f}=0$. The “$h$-$\Lambda$CDM” model is for computing the response with
respect to $h$ ($T_{h}$): here, we chose a step size of $\delta h=\pm 0.02$
for the numerical derivative. We also use the simulations for nonflat
$\Lambda$CDM models with $\Omega_{K}=\pm 0.1$, named as
“$\Omega_{K}$-$\Lambda$CDM2”, to assess how our method can approximate the
halo-matter and halo-auto power spectrum for nonflat models.

Table 1 gives the values of $\Omega_{K}$ and $h$, and we use the fixed values
of other cosmological parameters, given as $(\omega_{\rm c},\omega_{\rm
b},A_{s},n_{s})=(0.1198,0.02225,2.2065\times 10^{-9},0.9645)$, which specify
the transfer function and the primordial power spectrum, or equivalently the
linear matter power spectrum. Note that we also include the effect of massive
neutrinos on the linear matter power spectrum, assuming
$\Omega_{\nu}h^{2}=0.00064$ corresponding to $m_{\nu,{\rm tot}}=0.06\,{\rm
eV}$, the lower limit inferred from the oscillation experiments (see Ref.
Nishimichi  _et al._ , 2019, for details). Hence the physical density
parameter of total matter is $\Omega_{\rm m}h^{2}=\omega_{\rm c}+\omega_{\rm
b}+\omega_{\nu}$. Note that $\Omega_{\rm m}$ and $\Omega_{\Lambda}$ are
specified by a given set of the parameters for each model: $\Omega_{\rm
m}=\Omega_{{\rm m}f}h_{f}^{2}/h^{2}$ and $\Omega_{\Lambda}=1-\Omega_{\rm
m}-\Omega_{K})$. For each model, we use the outputs at 4 redshifts,
$z_{f}\simeq 0,0.55,1.03$ and $1.48$. Since the “fiducial” flat $\Lambda$CDM
model is the flat model in the SU method, each redshift for the fiducial flat
model corresponds to a slightly different redshift in each nonflat model,
which is computed from Eq. (5).

Furthermore, we use simulations that are run using the “paired-and-fixed”
method in Angulo and Pontzen (2016), where the initial density field in each
Fourier mode is generated from the fixed amplitude of the power spectrum
$\sqrt{P(k)}$ and the paired simulations with reverse phases, i.e.
$\delta_{{\bf k}}$ and $-\delta_{{\bf k}}$. The mean power spectrum of the
paired runs fairly well reproduces the ensemble average of many realizations
even in the nonlinear regime Angulo and Pontzen (2016); Villaescusa-Navarro
_et al._ (2018). The paired-and-fixed simulations allow us to significantly
reduce the sample variance in the power spectrum estimation. “2 (8)” on the
column $N_{\rm real}$ denotes one (four) pair(s) of the paired-and-fixed
simulations. For all simulations, halos are identified using Rockstar
(Behroozi  _et al._ , 2013). We adopt $M\equiv
M_{\Delta}=4\pi/3(R_{\Delta})^{3}(\Delta\bar{\rho}_{\rm m0})$ with
$\Delta=200$ throughout the paper. Note that when measuring the total response
in the $\Omega_{K}$-$\Lambda$CDM1 simulations, we use the spherical
overdensity $\Delta_{W}=200/(1+\delta_{\rm b}(t))$ so that halos are
identified using the same physical overdensity as in the corresponding global
universe (e.g. Li  _et al._ , 2016). Note that all the $N$-body simulations
for different cosmological models have a fixed mass resolution because they
share the same $\Omega_{\rm m}h^{2}$ and comoving volume.

After we identified halo candidates, we determine whether they are central or
satellite halos. When the separation of two different halos (between their
centers) is closer than $R_{\rm 200m}$ of the more massive one, we mark the
less massive one as a satellite halo. In the following, we use only central
halos.

###  IV.2 Mock catalogs of galaxies in $\Omega_{K}$-$\Lambda$CDM simulations

Table 2: The parameters of halo occupation distribution (HOD, Eq. 28) used for making the galaxy mock catalogs in this paper. $\log M_{\rm min}[M_{\odot}/h_{f}]$ | $13.94$
---|---
$\sigma_{\ln{M}}$ | $0.5$
$\log M_{0}[M_{\odot}/h_{f}]$ | $13.72$
$\log M_{1}[M_{\odot}/h_{f}]$ | $14.46$
$\alpha_{\rm sat}$ | $1.192$

To validate our method to compute the galaxy clustering observables such as
galaxy-matter and galaxy-auto power spectra in nonflat cosmology and their
total response, we build the galaxy mock catalogs using the halo catalogs of
$\Omega_{K}$-$\Lambda$CDM simulations. We use the halo catalogs of
$\Omega_{K}$-$\Lambda$CDM1 simulations with the halo mass defined using
$\Delta_{W}=200/(1+\delta_{\rm b}(t))$ to measure the total response, whereas
we use the halo catalogs of $\Omega_{K}$-$\Lambda$CDM2 simulations with the
halo mass defined using $\Delta=200$ to measure the galaxy clustering
observables in the nonflat universe. We assume that the galaxy profile around
halos and the mean HOD depend only on halo mass and are invariant to
cosmology, and populate galaxies into halos of each realization assuming the
same HOD in each $\Omega_{K}$-$\Lambda$CDM simulations. Specifically, we adopt
the following central and satellite HODs (Zheng  _et al._ , 2005):

$\displaystyle\langle N_{c}(M)\rangle$ $\displaystyle=\frac{1}{2}\left[1+{\rm
erf}\left(\frac{\log M-\log M_{\rm min}}{\sigma_{\ln M}}\right)\right],$
$\displaystyle\langle N_{S}(M)\rangle$
$\displaystyle=\Theta(M-M_{0})\left(\frac{M-M_{0}}{M_{1}}\right)^{\alpha_{\rm
sat}},$ (28)

where $\Theta$ is the step function. The total galaxy occupation is written as

$\displaystyle\langle N_{g}(M)\rangle=\langle N_{c}(M)\rangle\left[1+\langle
N_{s}(M)\rangle\right].$ (29)

We populate the central galaxies into the center of halos according to the
Bernoulli distribution with mean $\langle N_{c}(M)\rangle$. For satellite
galaxies, we populate them only into host halos each of which already hosts a
central galaxy, assuming they obey the Poisson distribution with mean $\langle
N_{S}(M)\rangle$. We assume that the satellite galaxy density profile
$u_{S}(r,M)$ follows the NFW profile (Navarro  _et al._ , 1996) with mass-
concentration relation $c(M)$ in Ref. Diemer and Kravtsov (2015). We employ
the HOD parameters as denoted in Table 2. These values correspond to the
fiducial vaules of SDSS “CMASS1” sample in Ref. Miyatake  _et al._ (2022b)
except for $\sigma_{\log{M}}$, for which we use $\sigma_{\log{M}}=0.5$ instead
of the fiducial value $\sigma_{\log{M}}=0.7919$.

As noted in the previous section, each pair of the SU simulation uses the same
initial seeds to reduce the sample variance. However, as we populated galaxies
randomly, the reduction of the sample variance is partially ruined. Hence, we
marginalized over 10 HOD seeds for each realization of $N$-body simulation to
obtain the converged results.

Figure 2: Growth responses of halo-matter and halo-auto power spectra with
respect to $\delta_{\rm b}$ ($h$), $T_{\delta_{\rm b}}^{\rm hm}(k)$
($T_{h}^{\rm hm}(k)$) and $T_{\delta_{\rm b}}^{\rm hh}(k)$ ($T_{h}^{\rm
hh}(k)$), at the four redshifts as denoted by the legend in each panel. We use
8 paired simulations for $\mbox{$\Omega_{K}$-$\Lambda$CDM}1$ and
$h$-$\Lambda$CDM models in Table 1 to compute these responses. The symbols
denote the mean of $T_{\delta_{\rm b}}$ or $T_{h}$ in each $k$ bin, and the
error bars (although not visible in some $k$ bins) denote the statistical
errors for the simulation box with side length $L=2\,h_{f}^{-1}{\rm Gpc}$,
which are estimated from the standard deviations among the 8 paired
simulations. Note that the range of $y$-axis is different in different panels.
For $T^{\rm hh}(k)$, we only plot the scales where halo power spectra after
subtracting the shot noise have positive values. We slightly shift the symbols
of $T_{h}(k)$ along the $x$-axis for illustration.

###  IV.3 Measurements of power spectra and the responses

To calculate the power spectrum from each simulation output, we assign the
$N$-body particles, halos, or galaxies on $2048^{3}$ grids using the cloud-in-
cells (CIC) method (Hockney and Eastwood, 1981) to obtain the density fields
of matter, halo or galaxy. After performing the Fourier transform, we correct
for the window function of CIC following the method described in Ref. Jing
(2005). We will show the results at wavenumbers smaller than the Nyquist
frequency, $k=3.2~{}h_{f}{\mathrm{Mpc}^{-1}}$. Furthermore, to evaluate the
power spectrum at $k\geq 1.6~{}h_{f}{\mathrm{Mpc}^{-1}}$ accurately, we fold
the particle positions into a smaller box by replacing
$\mathbb{x}\rightarrow\mathbb{x}\%(L/2)$, where the operation $a\%b$ stands
for the remainder of the division of $a$ by $b$. This procedure leads to
effectively $2$ times higher resolution. For the halo-auto and galaxy-auto
power spectrum, we subtract the shot noise from the measurements, where we
simply assume the shot noise to be given by the number density of the tracers.

Since we use the fixed box size and the same particle number, we use the same
$k$ binning to estimate the average of $|\delta_{\bf k}|^{2}$ in each $k$ bin
to estimate the band power. We then use the two-side numerical derivative
method to compute the power spectrum responses. We compute the growth response
as

$\displaystyle\left.\frac{\partial\ln P(k,z)}{\partial\delta_{\rm
b}}\right|_{{\rm G},\rm sim}$ $\displaystyle=\frac{\ln
P_{+}(k;\delta_{\rm{b}0}=+\epsilon)-\ln
P_{-}(k;\delta_{\rm{b}0}=-\epsilon)}{2\epsilon\times D(z)}$ (30)

and compute the total response as

$\displaystyle\left.\frac{\partial P(k,z)}{\partial\delta_{\rm b}}\right|_{\rm
total,sim}=(n-1)P(k,\delta_{\rm b0}=0)$
$\displaystyle+\frac{P_{+}((1-\delta_{\rm
b}/3)k;\delta_{\rm{b}0}=+\epsilon)-P_{-}((1+\delta_{\rm
b}/3)k;\delta_{\rm{b}0}=-\epsilon)}{2\epsilon\times D(z)},$ (31)

where $\delta_{\rm{b}0}\equiv\delta_{\rm b}(z_{f}=0)$ and $\epsilon=0.01$ for
the $\Omega_{K}$-$\Lambda$CDM1 simulations. To reduce statistical
stochasticity (or sample variance), we employ the same initial seeds as those
for the “fiducial” model. The column “$N_{\rm real}$” in Table 1 denotes the
number of realizations for paired simulations, where each pair uses the same
initial seeds. For $\Omega_{K}$-$\Lambda$CDM1 and $h$-$\Lambda$CDM models, we
run 8 paired simulations to estimate the statistical scatters.

##  V Results

###  V.1 Growth response

Figure 3: Comparison of the simulation result for $T_{h}^{\rm hm}(k)$,
$T_{h}^{\rm hh}(k)$ at $z_{f}=0.55$ with those computed using Dark Emulator
(Nishimichi  _et al._ , 2019), $T_{h}^{\rm hm,emu}(k)$, $T_{h}^{\rm
hh,emu}(k)$.

First, we study the approximate identity of the growth response
$T_{\delta_{\rm b}}(k)\approx T_{h}(k)$ with our $N$-body simulations.

In Fig. 2, we compare the growth response of halo-matter power spectrum to
$\delta_{\rm b}$ and $h$, for the abundance-matched halo samples, at the four
redshifts as in Table 1. These responses are calculated from $N$-body
simulations for $\Omega_{K}$-$\Lambda$CDM1 and $h$-$\Lambda$CDM models. We can
see that the approximate identity $T_{\delta_{\rm b}}^{\rm hm}(k)\approx
T_{h}(k)$ holds for all four redshifts and down to the nonlinear scale. Note
that halos are selected using the abundance matching method, where the same
number of halos, selected in the ascending order of masses, are identified in
the two simulations when computing the responses.

In Fig. 2, we also compare the growth response of the halo-auto power
spectrum. Again it is clear that the approximate identities $T_{\delta_{\rm
b}}^{\rm hh}(k)\approx T_{h}^{\rm hh}(k)$ hold for all four redshifts and down
to mildly nonlinear scales. Since the pairs of SU simulations
(“$\Omega_{K}$-$\Lambda$CDM1”) share the same Gaussian initial condition, the
abundance-matched halo samples correspond to the same initial density peaks
unless mergers severely affect this correspondence. Hence the pairs of SU
simulations have the similar clustering amplitudes of the abundance-matched
halos on large scales, which lead to $T^{\rm hh}(k\rightarrow 0)\simeq 0$. The
large change in the response $T^{\rm hh}(k)$ around $k\sim
1~{}h_{f}~{}\rm{Mpc}^{-1}$ can be attributed to the exclusion effects of the
halos, whereas the scales smaller than the exclusion scales does not
contribute to the galaxy-auto power spectrum.

We showed the results for only one case of the number density selected halo
sample, but we confirmed the approximate identity also holds for other number
densities.

As discussed in Ref. Terasawa  _et al._ (2022), we expect the responses
$T_{\delta_{\rm b}}(k)$ and $T_{h}(k)$ agree when the power spectrum is a
functional of the amplitude of the linear power spectrum. If we assume the
universal halo mass function, which depends on cosmology only through
$P_{L}(k)$, the corresponding halo bias is also determined only by $P_{L}(k)$.
Hence, we can expect the responses $T(k)$ to the different parameters that
leave the shape of the linear power spectrum unchanged agree with each other
for the linear or quasi-nonlinear scales, where bias expansion of the density
field is valid. Especially, on large scales, where $P_{\rm
hm}(k)=b_{1}P_{L}(k)$ ($P_{\rm hh}(k)=b_{1}^{2}P_{L}(k)$) holds, the response
$T(k)$ is expected to be constant. On the other hand, the responses in the
nonlinear regime could be different due to the difference in the growth
history of the structure (Terasawa  _et al._ , 2022), which leads to a change
in the concentration of halos (Li  _et al._ , 2014b). For $T_{\delta_{\rm
b}}(k)$ and $T_{h}(k)$, we found out they agree with each other even in the
nonlinear scale by numerical simulations.

In Fig. 3 we assess the accuracy of the responses of halo-matter and halo-auto
power spectra to $h$ calculated by Dark Emulator. The emulator can predict the
response accurately down to the nonlinear scale. For the halo-auto power
spectrum, it turned out that Dark Emulator does not predict the $h$ response
correctly, and hence we also show the response to $A_{s}$ and $z$ calculated
by Dark Emulator, $T_{\ln A_{s}}^{\rm hh,emu}(k)$ and $T_{z}^{\rm hh,emu}(k)$.
These responses approximate the response to $h$ well. Since the halo-auto
power spectrum only contains the 2-halo term, it is less affected by a change
in the concentration of halos and we can expect that the responses to $A_{s}$,
$z$, and $h$ have similar features.

Figure 4: An assessment of the accuracy of our method (Eq. 14) for predicting
the halo-matter power spectrum, $P_{\rm hm}(k)$ (left panel), and the halo-
halo power spectrum, $P_{\rm hh}(k)$ (right), for nonflat $\Lambda$CDM models.
The different symbols in each panel denote $P(k)$, directly estimated from
$N$-body simulations for nonflat models with $\Omega_{K}=\pm 0.1$
($\mbox{$\Omega_{K}$-$\Lambda$CDM}{2}$ models in Table 1), at $z_{f}=0.55$,
while the lines denote the results from our method, $\tilde{P}(k)$ in Eq. (14)
Note that we used the simulation results for $P^{f}(k,z_{f})$ and
$T_{h}(k,z_{f})$ in Eq. (14). For the simulation results we used the “paired-
and-fixed” method in Angulo and Pontzen (2016) to reduce the stochasticity,
and we considered the halo sample with number density $n_{\rm
h}=10^{-3}~{}h_{f}^{3}~{}{\rm Mpc}^{-3}$. We only plot the range of scales
where the halo power spectra after subtracting the shot noise have positive
values. For comparison, we also show the simulation result for the flat
fiducial simulation by triangle symbols. The lower plot in each panel shows
the ratio between the simulation result and our method. The horizontal solid
and dashed lines denote $\pm 2,\pm 1\%$ fractional accuracy, respectively.

Figure 5: An assessment of the accuracy of our method (Eq. 14) for predicting
the galaxy-matter power spectrum, $P_{\rm gm}(k)$ (left panel), and the
galaxy-auto power spectrum, $P_{\rm gg}(k)$ (right), for nonflat $\Lambda$CDM
models. The different symbols in each panel denote $P(k)$, directly estimated
from $N$-body simulations for flat and nonflat models with $\Omega_{K}=\pm
0.1$, at $z_{f}=0.55$, while the lines denote the prediction using our method,
$\tilde{P}(k)$ in Eq. (14) together with Dark Emulator. Error bars are
estimated from the standard deviations among the 20 mock catalogs. The lower
plot in each panel shows the ratio between the simulation result and our
method. The horizontal solid and dashed lines denote $\pm 5,\pm 2\%$
fractional accuracy, respectively.

###  V.2 SU approach for $P(k;\Omega_{K})$

To assess the performance of our estimator for the halo-matter and halo-auto
power spectra for nonflat universe $\tilde{P}(k;\Omega_{K})$ (Eq. 14), we
compare it with the power spectra measured from $N$-body simulations with
$\Omega_{K}=\pm 0.1$ ($\Omega_{K}$-$\Lambda$CDM2 model).

In Fig. 4, the data points show the power spectrum measured from simulations,
for $f\Lambda$CDM model and $\Omega_{K}$-$\Lambda$CDM2 models. The curves show
the predictions of Eq. 14, where we used $P^{f}(k,z_{f})$ and $T_{h}(k,z_{f})$
measured from simulations of $f\Lambda$CDM and $h$-$\Lambda$CDM models,
respectively. For the halo-matter power spectrum (left panel), the estimator
has $\sim 1\%$ accuracy even for such large curvature, $\Omega_{K}=+0.1,-0.1$
(corresponding to $\delta_{\rm b}(z_{f}=0.55)=-0.12,+0.10$). For the halo-auto
power spectrum (right panel), the estimator has $\sim 2\%$ accuracy up to
$k\simeq 1~{}h_{f}\rm{Mpc}^{-1}$. The relatively large deviation at $k\gtrsim
1~{}h_{f}\rm{Mpc}^{-1}$ is due to both large deviation from $P^{f}(k,z_{f})$
and the inaccuracy of approximation, $T_{\delta_{\rm b}}(k)\approx T_{h}(k)$.
The smaller fractional change in the amplitudes of $P_{\rm hm}$ and $P_{\rm
hh}$ by the non-zero $\Omega_{K}$ than the change in $P_{\rm mm}$ (see Fig. 5
in Ref. Terasawa  _et al._ , 2022) is ascribed to the result of Fig. 1 (see
the discussion around the figure).

Further, we tested the accuracy of the prediction for the galaxy-matter and
galaxy-auto power spectra by Dark Emulator, using the estimator in Eq. (14).
We calculated these spectra as the weighted integral of $P_{\rm hm}(k;>M)$ or
$P_{\rm hh}(k;>M)$ as in Eq. 37. For the nonflat cosmology, we use the Eq. 14
to estimate $P_{\rm hm}(k;>M)$ and $P_{\rm hh}(k;>M)$, using $P^{f}(k,z_{f})$
and $T_{h}(k,z_{f})$ predicted by the emulator. We used the halo mass function
predicted by the emulator for the flat model, while we used the model in Ref
(Tinker  _et al._ , 2008) for the nonflat model because the emulator is not
trained for the nonflat model. On the large scale, we stitched the predictions
with those of the linear theory, similarly as Eq. 27.

In Fig. 5, we can see the predictions of the galaxy-matter power spectrum by
the emulator, $P^{\rm emu}(k)$ have $\sim 5\%$ accuracy. Compared to the
prediction for the flat universe, we can extend the prediction to the non-zero
curvature without significant degradation.

For the galaxy-auto power spectrum, due to the poor accuracy of calculating
the $h$ response by the emulator (Fig. 3), we use $T_{\ln A_{s}}^{\rm hh}(k)$
instead. The accuracy of predicting the galaxy-auto power spectra in nonflat
universe is degraded compared to that for the flat universe, but the estimator
still can predict the spectra with $\sim 10\%$ accuracy.

Figure 6: Total response of matter power spectrum at $z_{f}=0.55$ measured
from the simulations and those calculated using our method (Eq.
LABEL:eq:p_mm_total; labeled as “This work”), perturbation theory (Eq. 16) and
halo model (Eq. 20). The error bars (although not visible) are estimated from
the standard deviations among the 8 paired simulations.

Figure 7: Total response of power spectrum at $z_{f}=0.55$ measured from the
simulations and those calculated using our method (Eqs. 22, 25, and 26);
labeled as “This work”), perturbation theory (Eq. 16) and halo model. The
dotted and dashed-dotted lines denote the halo model (+HOD) prediction with or
without $I_{g}^{2}(k)$ term in Eq. (20). Left: the response of galaxy-matter
power spectrum. Right: the response of galaxy-auto power spectrum. The error
bars are estimated from the standard deviations among the 8 paired
simulations.

###  V.3 Total response

We tested the accuracy of our estimator for the power spectrum total response
to the super-survey modes compared with the measurement from the mock galaxy
catalog and $N$-body results.

In Fig. 6 we show the accuracy of the matter power spectrum total-response at
$z_{f}=0.55$ calculated by our method (Eq. LABEL:eq:p_mm_total) against the
$N$-body simulations, along with other theoretical predictions. We use Halofit
in Smith et al. (Smith+03) (Smith  _et al._ , 2003) and Takahashi et al.
(Takahashi+12) (Takahashi  _et al._ , 2012) to compute $T_{h}^{\rm mm}(k)$ and
$P_{\rm mm}(k)$ in Eq. (LABEL:eq:p_mm_total), respectively. This different
choice of Halofit is because Smith+03 predicts $T_{h}^{\rm mm}(k)$ better than
Takahashi+12 (See the Fig. 3 in (Terasawa  _et al._ , 2022)), while
Takahashi+12 is known to be more accurate for $P_{\rm mm}(k)$ especially on
small scales. The perturbation theory prediction agrees with the $N$-body
results on large scales, but it begins to deviate around $k\gtrsim
0.1\,h_{f}{\rm Mpc}^{-1}$. The halo model reproduces the behavior of the
response out to the nonlinear scale, but it underestimates the response by
$\sim 20\%$ at $k\sim 1\,h{\rm Mpc}^{-1}$. Compared with these two analytical
models which are commonly used in the literature, our method predicts the
response better over a wide range of scales.

In Fig. 7 we show the accuracy of the galaxy-matter and galaxy-auto power
spectrum total response at $z_{f}=0.55$ calculated by our method (Eqs. 22, 25,
and 26.) against the $N$-body simulations, along with other theoretical
predictions. We use Dark Emulator to compute $P(k)$, $T_{h}(k)$, and mass
function in Eqs. (22), (25), and (26). For the galaxy-auto power spectrum, as
discussed in the Sec. V.1, since Dark Emulator does not predict the $h$
response of the halo-auto power spectrum correctly, we used the response to
$A_{s}$ calculated by Dark Emulator instead. As with the matter power spectrum
response, our method can predict the galaxy-matter and galaxy-auto power
spectrum responses more accurately compared to the analytical ones in the
literature.

Lastly, we mention the linear limit of the halo model predictions. In Figs. 6
and 7, as we discussed in Sec. III.1, we can see the predictions of the halo
model converge to those of the perturbation theory at large scale. We also
show the impact of ignoring the $I_{\rm g}^{2}(k)$ term in the halo model
prediction. It is clear that ignoring the term results in an inaccurate
prediction on linear and quasi-nonlinear scales.

##  VI conclusion

In this paper, we have developed an approximate method to model the halo-
matter and halo-auto power spectra for the nonflat $\Lambda$CDM model, from
quantities representing the nonlinear evolution of the corresponding flat
$\Lambda$CDM model, based on the SU method. The key points to build the
estimator are the correspondence between the nonflat and flat universes
through the SU picture and the equivalence of the growth responses to long-
wavelength modes and the Hubble parameter. This work is a sequel of our
previous research Terasawa  _et al._ (2022), in which we proposed the
approximate method for predicting the nonlinear matter power spectrum for
nonflat $\Lambda$CDM model.

The estimator of the halo-matter (halo-auto) power spectrum has $\sim 1\%$
($\sim 2\%$) accuracy even for large curvature model with $\Omega_{K}=\pm
0.1$. Using the estimator we can extend the existing emulators to predict the
nonlinear power spectra for nonflat universe without degrading its accuracy.
In particular, we showed we can extend the Dark Emulator to predict the
galaxy-matter and galaxy-auto power spectra for a non-zero curvature model.

The response $T_{\delta_{\rm b}}(k)$ is also a key quantity for estimating
SSC. We utilized the approximate identity $T_{\delta_{\rm b}}(k)\approx
T_{h}(k)$ and proposed the calculation of the total response (SSC terms) using
Halofit or Dark Emulator. We showed that our method can predict the total
response at an accuracy better than the analytical methods used in the
literature such as the perturbation theory and the halo model, thanks to
capturing the nonlinear response through $T_{h}(k)$. Our method for computing
the total response would be implemented into the Core Cosmology Library
(Chisari  _et al._ , 2019).

Although we assume $\Lambda$CDM model in this paper, we expect our approach to
approximate the response to $\delta_{\rm b}$ is also applicable to $w$CDM
model, where dark energy is not the cosmological constant. Once validated with
numerical simulations, we can ease the computational cost of the SU
simulations by substituting it with the response to another parameter which
can be evaluated fast. In this paper, we focused on the responses to
$\delta_{\rm b}$ and $h$, but we can think of other parameters that share the
almost same response in wider cosmological parameter space, which can be used
to extend the simulation-based theory predictions and to ease the
computational cost for building such models.

###### Acknowledgements.

We would like to thank Kaz Akitsu, Yosuke Kobayashi, Yue Nan, Tim Eifler, and
Elisabeth Krause for their useful and stimulating discussion. This work was
supported in part by World Premier International Research Center Initiative
(WPI Initiative), MEXT, Japan, JSPS KAKENHI Grant Numbers JP23KJ0747,
JP22H00130, JP20H05850, JP20H05855, JP20H05861, JP20H04723, JP19H00677,
JP21H01081, JP22K03634, Basic Research Grant (Super AI) of Institute for AI
and Beyond of the University of Tokyo, and Japan Science and Technology Agency
(JST) AIP Acceleration Research Grant Number JP20317829. Numerical
computations were carried out on Cray XC50 at Center for Computational
Astrophysics, National Astronomical Observatory of Japan.

##  Appendix A Calculation methods of the response functions of halo and
galaxy power spectra

###  A.1 The response of galaxy-matter power spectrum

We define the cumulative number density of halos as

$\displaystyle n(>M_{\rm th};\delta_{\rm b})\equiv\int_{\rm
M_{th}}^{\infty}\frac{\mathrm{d}M}{M}n_{\ln M}(M;\delta_{\rm b}),$ (32)

where $n_{\ln M}\equiv\mathrm{d}n/\mathrm{d}\ln M$. We define the following
dimensionless quantities:

$\displaystyle\mathscr{P}_{\mathrm{hm}}(k;>M_{\mathrm{th}})\equiv
n(>M_{\mathrm{th}})P_{\mathrm{hm}}(k;>M_{\mathrm{th}}),$ (33)

for mass-threshold samples (above mass $M_{\mathrm{th}}$) and

$\displaystyle\mathscr{P}_{\mathrm{hm}}(k;M_{-},M_{+})\equiv
n(M_{-},M_{+})P_{\mathrm{hm}}(k;M_{-},M_{+}),$ (34)

for mass bin samples in the mass range of $[M_{-},M_{+}]$, respectively. Using
these quantities, all the equations for $P_{\mathrm{gm}}$ and its response can
be simplified. The idea behind this is as follows. First, mass-threshold halo
samples are easier to analyze in simulations. Once the halo-matter power
spectrum for various different mass thresholds is known, one can convert them
to that for mass-bin samples by taking a derivative w.r.t. the mass threshold.
Numerically, this can be done by taking the finite difference:

$P_{\mathrm{hm}}(k;M)\simeq\dfrac{n(>M_{-})P_{\mathrm{hm}}(k;>M_{-})-n(>M_{+})P_{\mathrm{hm}}(k;>M_{+})}{n(>M_{-})-n(>M_{+})},$
(35)

where $M_{\pm}=M\pm\epsilon_{M}$ with some small $\epsilon_{M}$ compared to
$M_{\pm}$. In the above, note that the halo-matter power spectrum is a halo
number-weighted quantity, and thus $n$ must be considered appropriately. This
equation is actually implemented in Dark Emulator. Now, using
$\mathscr{P}_{\mathrm{hm}}$, this relation is simplified as

$\displaystyle\mathscr{P}_{\mathrm{hm}}(k;M_{-},M_{+})=\mathscr{P}_{\mathrm{hm}}(k;>M_{-})-\mathscr{P}_{\mathrm{hm}}(k;>M_{+}).$
(36)

We can use a similar trick for the galaxy-matter power spectrum. Assuming a
halo occupation distribution model, we can have

$\displaystyle
P_{\mathrm{gm}}(k)=\frac{1}{\bar{n}_{\mathrm{g}}}\int\mathrm{d}M\langle\tilde{U}_{g}(k;M)\rangle\dfrac{\mathrm{d}n}{\mathrm{d}M}(M)P_{\mathrm{hm}}(k;M),$
(37)

where $\tilde{U}_{g}(k;M)$ describes the Fourier transform of the average
radial profile of galaxy number density in host halos with mass $M$, which can
be computed by taking the product of the halo occupation distribution
$N_{\mathrm{g}}(M)$ and the radial profile of each galaxy.

This can be evaluated by substituting Eq. (35). Now, multiplying both sides of
Eq. (37) by $\bar{n}_{\mathrm{g}}$ and defining
$\mathscr{P}_{\mathrm{gm}}=\bar{n}_{\mathrm{g}}P_{\mathrm{gm}}$, one can
compute

$\displaystyle\mathscr{P}_{\mathrm{gm}}(k)$ $\displaystyle=$
$\displaystyle\int\mathrm{d}M\langle\tilde{U}_{g}(k;M)\rangle\dfrac{\mathrm{d}}{\mathrm{d}M}\Bigl{(}-\mathscr{P}_{\mathrm{hm}}(k;>M)\Bigr{)}$
(38) $\displaystyle=$
$\displaystyle-\langle\tilde{U}_{g}(k;M)\rangle\mathscr{P}_{\mathrm{hm}}(k;>M)\Bigr{|}_{0}^{\infty}+\int\mathrm{d}M\dfrac{\mathrm{d}\langle\tilde{U}_{g}(k;M)\rangle}{\mathrm{d}M}\mathscr{P}_{\mathrm{hm}}(k;>M),$
$\displaystyle=$
$\displaystyle\int\mathrm{d}M\dfrac{\mathrm{d}\langle\tilde{U}_{g}(k;M)\rangle}{\mathrm{d}M}\mathscr{P}_{\mathrm{hm}}(k;>M).$

In the above, we have performed an integration by parts to obtain the second
line and then use the fact that $N_{\mathrm{g}}(M)$ tends to zero at the low
mass end and $\mathscr{P}_{\mathrm{hm}}(M)$ tends to zero at the high mass end
to reach the final line. This is how we can avoid estimating the power
spectrum for the mass-bin samples and directly compute the galaxy-matter power
spectrum from the mass-threshold halo samples.

The responses are also related by a similar equation as

$\displaystyle\left.\frac{\partial\mathscr{P}_{\mathrm{gm}}(k)}{\partial\delta_{\mathrm{b}}}\right|_{\rm
G}=\int\mathrm{d}M\left[\frac{\mathrm{d}\langle\tilde{U}_{g}(k;M)\rangle}{\mathrm{d}M}\left.\frac{\partial\mathscr{P}_{\mathrm{hm}}(k;>M)}{\partial\delta_{\mathrm{b}}}\right|_{{\rm
G},M}+\frac{\partial}{\partial\delta_{\mathrm{b}}}\left.\frac{\mathrm{d}\langle\tilde{U}_{g}(k;M)\rangle}{\mathrm{d}M}\right|_{\rm
G}\mathscr{P}_{\mathrm{hm}}(k;>M)\right],$ (39)

where the subscript $M$ denotes the derivative with fixed mass threshold.
Assuming that the cosmological dependence is fully encoded in
$\mathscr{P}_{\mathrm{gm}}$ and that the galaxy profile around halos in
physical scale and the mean HOD are given irrespective of the cosmological
model, which is partly validated in Ref (Voivodic and Barreira, 2021), one can
trivially derive

$\displaystyle\frac{\partial}{\partial\delta_{\mathrm{b}}}\left.\frac{\mathrm{d}\langle\tilde{U}_{g}(k;M)\rangle}{\mathrm{d}M}\right|_{\rm
G}=\frac{1}{3}\frac{\partial}{\partial\ln
k}\frac{\mathrm{d}\langle\tilde{U}_{g}(k;M)\rangle}{\mathrm{d}M}.$ (40)

Finally, one can convert this to the standard response function by

$\displaystyle\left.\dfrac{\partial
P_{\mathrm{gm}}(k)}{\partial\delta_{\mathrm{b}}}\right|_{\rm
G}=\frac{1}{\bar{n}_{\mathrm{g}}}\left.\dfrac{\partial\mathscr{P}_{\mathrm{gm}}(k)}{\partial\delta_{\mathrm{b}}}\right|_{\rm
G}-b_{\mathrm{g}}^{\mathrm{L}}P_{\mathrm{gm}}(k),$ (41)

where $b_{g}^{L}$ is the Lagrangian galaxy bias, related to the Eulerian bias
as $b_{g}=b_{g}^{L}+1$. We can give a physical interpretation of the bias
factor as follows. The Lagrangian galaxy bias, $b_{g}^{L}$, describes the
change in the number density of galaxies at a fixed Lagrangian, or,
equivalently, local comoving volume in the SU picture, and the difference
between the Eulerian and the Lagrangian bias, $b_{g}-b_{g}^{L}=1$ represents
the contribution from the change in the physical volume, i.e., equivalent to
the dilation effect. Note that the second term in Eq. 41, $-b_{g}^{L}P_{\rm
gm}(k)$ comes from the response of $\bar{n}_{g}$ and since the growth response
should be evaluated at a fixed comoving volume, the response of $\bar{n}_{g}$
gives $b_{g}^{L}$ instead of $b_{g}$.

###  A.2 The response of galaxy auto-power spectrum

Similarly, we can also evaluate the galaxy-auto power spectrum using mass-
threshold quantities. First, the halo-auto power spectrum of the mass-bin
samples can be calculated by taking the finite difference:

$\displaystyle P_{\rm hh}(k;M_{1},M_{2})=\frac{\frac{\partial^{2}}{\partial
M\partial M^{\prime}}\left.\left[n(>M)n(>M^{\prime})P_{\rm
hh}(k;>M,>M^{\prime})\right]\right|_{M=M_{1},M^{\prime}=M_{2}}}{\frac{\mathrm{d}n}{\mathrm{d}M}(M_{1})\frac{\mathrm{d}n}{\mathrm{d}M}(M_{2})}.$
(42)

Following the halo model description, we split the galaxy power spectrum into
the 2-halo and 1-halo terms as $P_{\rm gg}(k)=P_{\rm gg}^{2h}(k)+P_{\rm
gg}^{1h}(k)$. Using Eq. 42, we can write the 2-halo term of the galaxy power
spectrum as

$\displaystyle P_{\rm gg}^{2h}(k)$
$\displaystyle=\frac{1}{\bar{n}_{g}^{2}}\int\mathrm{d}M\frac{\mathrm{d}n}{\mathrm{d}M}(M)\langle\tilde{U}_{g}(k;M)\rangle\int\mathrm{d}M^{\prime}\frac{\mathrm{d}n}{\mathrm{d}M^{\prime}}(M^{\prime})\langle\tilde{U}_{g}(k;M^{\prime})\rangle
P_{\rm hh}(k;M,M^{\prime})$
$\displaystyle=\frac{1}{\bar{n}_{g}^{2}}\int\mathrm{d}M\frac{\mathrm{d}\langle\tilde{U}_{g}(k;M)\rangle}{\mathrm{d}M}n(>M)\int\mathrm{d}M^{\prime}\frac{\mathrm{d}\langle\tilde{U}_{g}(k;M^{\prime})\rangle}{\mathrm{d}M^{\prime}}n(>M^{\prime})P_{\rm
hh}(k;>M,>M^{\prime})$
$\displaystyle=\frac{1}{\bar{n}_{g}^{2}}\iint\mathrm{d}M\mathrm{d}M^{\prime}\frac{\mathrm{d}\langle\tilde{U}_{g}(k;M)\rangle}{\mathrm{d}M}\frac{\mathrm{d}\tilde{U}_{g}(k;M^{\prime})}{\mathrm{d}M^{\prime}}\mathscr{P}_{\rm
hh}(k;>M,>M^{\prime}),$ (43)

where we did integral by parts and ignored the surface terms, and we define

$\displaystyle\mathscr{P}_{\rm hh}(k;>M_{\rm th},>M^{\prime}_{\rm th})$
$\displaystyle\equiv n(>M_{\rm th})n(>M^{\prime}_{\rm th})P_{\rm hh}(k;>M_{\rm
th},>M^{\prime}_{\rm th}).$ (44)

The SU growth responses of the 2-halo term of the galaxy-auto power spectrum
can be calculated as

$\displaystyle\left.\frac{\partial P_{\rm gg}^{2h}(k)}{\partial\delta_{\rm
b}}\right|_{\rm G}$
$\displaystyle=\frac{1}{\bar{n}_{g}^{2}}\iint\mathrm{d}M\mathrm{d}M^{\prime}\frac{\mathrm{d}\tilde{U}_{g}(k;M^{\prime})}{\mathrm{d}M^{\prime}}\left[\frac{\mathrm{d}\langle\tilde{U}_{g}(k;M)\rangle}{\mathrm{d}M}\left.\frac{\partial\mathscr{P}_{\rm
hh}(k;>M,>M^{\prime})}{\partial\delta_{\rm b}}\right|_{{\rm
G},M,M^{\prime}}\right.$
$\displaystyle~{}~{}~{}~{}~{}+\left.\left.2\frac{\partial}{\partial\delta_{\rm
b}}\frac{\mathrm{d}\langle\tilde{U}_{g}(k;M)\rangle}{\mathrm{d}M}\right|_{\rm
G}{\mathscr{P}_{\rm hh}(k;>M,>M^{\prime})}\right]-2b_{g}^{L}P_{\rm
gg}^{2h}(k).$ (45)

For the 1-halo term, we can compute the SU growth response as

$\displaystyle\left.\frac{\partial P_{\rm gg}^{1h}(k)}{\partial\delta_{\rm
b}}\right|_{\rm G}$
$\displaystyle=\frac{1}{\bar{n}_{g}^{2}}\int\mathrm{d}M\frac{\mathrm{d}n}{\mathrm{d}M}(M)\left[b_{\mathrm{h},1}^{L}(M)\langle\tilde{U}_{g}^{2}(k;M)\rangle+\left.\frac{\partial\langle\tilde{U}_{g}^{2}(k;M)\rangle}{\partial\delta_{\rm
b}}\right|_{\rm G}\right]-2b_{g}^{L}P_{\rm gg}^{1h}(k).$ (46)

Note that the response of the 1-halo term in our approach is identical to that
of the halo model description (e.g. Krause and Eifler, 2017; Nicola  _et al._
, 2020).

###  A.3 The response calibration of halo and galaxy power spectra using the
abundance matching method

In this subsection, we describe the relation between the SU growth response
and the number density fixed or abundance matched response which we clarify
later. We denote the derivative performed keeping the comoving halo number
density fixed as “AM” response, which stands for Abundance Matching method.
When evaluating the AM response, we change the threshold $M_{\rm
th}(\delta_{\rm b})$ so that the cumulative halo number density in the
comoving volume is kept fixed when varying $\delta_{\rm b}$ (Li  _et al._ ,
2016):

$\displaystyle\frac{\mathrm{d}n(>M_{\rm th};\delta_{\rm
b})}{\mathrm{d}\delta_{\rm b}}=0.$ (47)

We define the mass threshold shift $s(M)$ to keep the number density as

$\displaystyle\mathrm{d}\ln M_{\rm th}\equiv s(M_{\rm
th})\mathrm{d}\delta_{\rm b}.$ (48)

We can calculate Eq. (47) as

$\displaystyle\frac{\mathrm{d}n(>M_{\rm th};\delta_{\rm
b})}{\mathrm{d}\delta_{\rm b}}$
$\displaystyle=\frac{\mathrm{d}}{\mathrm{d}\delta_{\rm b}}\int_{\rm
M_{th}(\delta_{\rm b})}^{\infty}\frac{\mathrm{d}M}{M}n_{\ln M}(M;\delta_{\rm
b})$ $\displaystyle=\int_{\rm
M_{th}}^{\infty}\frac{\mathrm{d}M}{M}\frac{\partial n_{\ln M}(M;\delta_{\rm
b})}{\partial\delta_{\rm b}}+\frac{\partial}{\partial\delta_{\rm
b}}\left[F_{1}(M=\infty)-F_{1}(M=M_{\rm th})\right]$ $\displaystyle=\int_{\rm
M_{th}}^{\infty}\frac{\mathrm{d}M}{M}\frac{\partial n_{\ln M}(M;\delta_{\rm
b})}{\partial\delta_{\rm b}}-n_{\ln M}(M_{\rm th})s(M_{\rm th})$
$\displaystyle=0,$ (49)

where we defined $F_{1}(M)$ as the function satisfying $\partial
F_{1}(M)/\partial\ln M=n_{\ln M}(M)$ and we can calculate its response as

$\displaystyle\frac{\partial F_{1}(M)}{\partial\delta_{\rm b}}$
$\displaystyle=\frac{\partial F_{1}(M)}{\partial\ln M}\frac{\partial\ln
M}{\partial\delta_{\rm b}}$ $\displaystyle=n_{\ln M}(M)s(M).$ (50)

Hence, we have

$\displaystyle\int_{\rm M_{th}}^{\infty}\frac{\mathrm{d}M}{M}\frac{\partial
n_{\ln M}(M;\delta_{\rm b})}{\partial\delta_{\rm b}}=n_{\ln M}(M_{\rm
th})s(M_{\rm th}).$ (51)

Using this, the Lagrangian halo bias above the mass threshold is given as (Li
_et al._ , 2016)

$\displaystyle b_{\mathrm{h},1}^{L}(>M_{\rm th})$
$\displaystyle\equiv\frac{1}{n(>M_{\rm th})}\int_{M_{\rm
th}}^{\infty}\frac{\mathrm{d}M}{M}b_{\mathrm{h},1}^{L}n_{\ln M}$
$\displaystyle=\frac{1}{n(>M_{\rm th})}\int_{M_{\rm
th}}^{\infty}\frac{\mathrm{d}M}{M}\frac{\partial\ln n_{\ln
M}}{\partial\delta_{\rm b}}n_{\ln M}$ $\displaystyle=\frac{n_{\ln M}(M_{\rm
th})s(M_{\rm th})}{n(>M_{\rm th})}.$ (52)

Likewise, the AM response of $\mathscr{P}_{\rm hm}(k;>M_{\rm th})$ can be
calculated as

$\displaystyle\left.\frac{\mathrm{d}\mathscr{P}_{\rm hm}(k;>M_{\rm
th})}{\mathrm{d}\delta_{\rm b}}\right|_{\rm G}$
$\displaystyle=\frac{\mathrm{d}}{\mathrm{d}\delta_{\rm b}}\int_{\rm
M_{th}}^{\infty}\mathrm{d}\ln Mn_{\ln M}(M)P_{\rm hm}(k;M)$
$\displaystyle=\int_{\rm M_{th}}^{\infty}\mathrm{d}\ln
M\frac{\partial~{}}{\partial\delta_{\rm b}}\left[n_{\ln M}(M)P_{\rm
hm}(k;M)\right]+\frac{\mathrm{d}}{\mathrm{d}\delta_{\rm
b}}\left[F_{2}(M=\infty)-F_{2}(M=M_{\rm th})\right]$
$\displaystyle=\left.\frac{\partial\mathscr{P}_{\rm hm}(k;>M_{\rm
th})}{\partial\delta_{\rm b}}\right|_{{\rm G},M_{\rm th}}-n_{\ln M}(M_{\rm
th})P_{\rm hm}(k;M_{\rm th})s(M_{\rm th})$
$\displaystyle=\left.\frac{\partial\mathscr{P}_{\rm hm}(k;>M_{\rm
th})}{\partial\delta_{\rm b}}\right|_{{\rm G},M_{\rm
th}}-b_{\mathrm{h},1}^{L}(>M_{\rm th})n(>M_{\rm th})P_{\rm hm}(k;M_{\rm th}),$
(53)

where $F_{2}(M)$ is a function satisfying $\partial F_{2}(M)/\partial\ln
M=n_{\ln M}(M)P_{\rm hm}(k;M)$. On the other hand, the l.h.s. of Eq. (53) can
be written as

$\displaystyle\left.\frac{\mathrm{d}\mathscr{P}_{\rm hm}(k;>M_{\rm
th})}{\mathrm{d}\delta_{\rm b}}\right|_{\rm G}$ $\displaystyle=n(>M_{\rm
th})\left.\frac{\partial P_{\rm hm}(k;>M_{\rm th})}{\partial\delta_{\rm
b}}\right|_{\rm G}$ $\displaystyle=n(>M_{\rm
th})\left(\frac{26}{21}\right)P_{\rm hm}(k;>M_{\rm th})\,T_{\delta_{\rm
b}}^{\rm hm}(k;>M_{\rm th})$ $\displaystyle\simeq n(>M_{\rm
th})\left(\frac{26}{21}\right)P_{\rm hm}(k;>M_{\rm th})\,T_{h}^{\rm
hm}(k;>M_{\rm th}),$ (54)

where we used the approximate relation $T_{h}^{\rm hm}(k;>M_{\rm th})\approx
T_{\delta_{\rm b}}^{\rm hm}(k;>M_{\rm th})$ in the last equation. Hence, we
can calculate the $P_{\rm gm}(k)$ SU-growth response using $T_{h}(k;>M_{\rm
th})$ as

$\displaystyle\left.\frac{\partial P_{\rm gm}(k)}{\partial\delta_{\rm
b}}\right|_{\rm G}$
$\displaystyle=\frac{1}{\bar{n}_{g}}\left.\frac{\partial\mathscr{P}_{\mathrm{gm}}(k)}{\partial\delta_{\rm
b}}\right|_{\rm G}-b_{g}^{L}P_{\rm gm}(k)$
$\displaystyle=\frac{1}{\bar{n}_{g}}\int\mathrm{d}M\left\\{\dfrac{\mathrm{d}\langle\tilde{U}_{g}(k;M)\rangle}{\mathrm{d}M}\left.\frac{\partial\mathscr{P}_{\mathrm{hm}}(k;>M)}{\partial\delta_{\rm
b}}\right|_{{\rm
G},M}+\frac{1}{3}\mathscr{P}_{\mathrm{hm}}(k;>M)\frac{\partial}{\partial\ln
k}\frac{\mathrm{d}\langle\tilde{U}_{g}(k;M)\rangle}{\mathrm{d}M}\right\\}-b_{g}^{L}P_{\rm
gm}(k)$
$\displaystyle=\frac{1}{\bar{n}_{g}}\int\mathrm{d}M\left\\{\dfrac{\mathrm{d}\langle\tilde{U}_{g}(k;M)\rangle}{\mathrm{d}M}n(>M)\left[b_{\mathrm{h},1}^{L}(>M)P_{\rm
hm}(k;M)+\left(\frac{26}{21}\right)T_{h}^{\rm hm}(k;>M)P_{\rm
hm}(k;>M)\right]\right.$
$\displaystyle~{}~{}~{}~{}~{}+\left.\frac{1}{3}P_{\mathrm{hm}}(k;>M)\frac{\partial}{\partial\ln
k}\frac{\mathrm{d}\ln\langle\tilde{U}_{g}(k;M)\rangle}{\mathrm{d}M}\right\\}-b_{g}^{L}P_{\rm
gm}(k).$ (55)

The AM response of $\mathscr{P}_{\rm hh}(k;>M_{\rm th},>M^{\prime}_{\rm th})$
is calculated as

$\displaystyle\left.\frac{\mathrm{d}\mathscr{P}_{\rm
hh}}{\mathrm{d}\delta_{\rm b}}\right|_{\rm G}$
$\displaystyle=\left.\frac{\mathrm{d}\mathscr{P}_{\rm hh}(k;>M_{\rm
th},>M^{\prime}_{\rm th})}{\mathrm{d}\delta_{\rm b}}\right|_{{\rm G},M_{\rm
th},M^{\prime}_{\rm th}}+\int_{M_{\rm th}}\mathrm{d}\ln Mn_{\ln
M}(M)\left\\{\frac{\partial}{\partial\delta_{\rm
b}}\left[F_{3}(M^{\prime}=\infty)-F_{3}(M^{\prime}=M^{\prime}_{\rm
th})\right]\right\\}$ $\displaystyle~{}~{}~{}~{}~{}+\int_{M^{\prime}_{\rm
th}}\mathrm{d}\ln M^{\prime}n_{\ln
M}(M^{\prime})\left\\{\frac{\partial}{\partial\delta_{\rm
b}}\left[F_{3}(M=\infty)-F_{3}(M=M_{\rm th})\right]\right\\}$
$\displaystyle=\left.\frac{\mathrm{d}\mathscr{P}_{\rm hh}(k;>M_{\rm
th},>M^{\prime}_{\rm th})}{\mathrm{d}\delta_{\rm b}}\right|_{{\rm G},M_{\rm
th},M^{\prime}_{\rm th}}-b_{\mathrm{h},1}^{L}(>M^{\prime}_{\rm
th})n(>M^{\prime}_{\rm th})\int_{M_{\rm th}}\mathrm{d}\ln Mn_{\ln M}(M)P_{\rm
hh}(k;M,M^{\prime}_{\rm th})$
$\displaystyle~{}~{}~{}~{}~{}-b_{\mathrm{h},1}^{L}(>M_{\rm th})n(>M_{\rm
th})\int_{M^{\prime}_{\rm th}}\mathrm{d}\ln M^{\prime}n_{\ln
M}(M^{\prime})P_{\rm hh}(k;M^{\prime},M_{\rm th})$
$\displaystyle=\left.\frac{\mathrm{d}\mathscr{P}_{\rm hh}(k;>M_{\rm
th},>M^{\prime}_{\rm th})}{\mathrm{d}\delta_{\rm b}}\right|_{{\rm G},M_{\rm
th},M^{\prime}_{\rm th}}-n(>M_{\rm th})n(>M^{\prime}_{\rm th})$
$\displaystyle~{}~{}~{}~{}~{}\times\left[b_{\mathrm{h},1}^{L}(>M^{\prime}_{\rm
th})P_{\rm hh}(k;>M_{\rm th},M^{\prime}_{\rm th})+b_{\mathrm{h},1}^{L}(>M_{\rm
th})P_{\rm hh}(k;>M^{\prime}_{\rm th},M_{\rm th})\right],$ (56)

where $F_{3}(M)$ satisfying $\partial F_{3}(M)/\partial\ln M=n_{\ln
M}(M)P_{\rm hh}(k;M,M^{\prime})$. Since the l.h.s. of Eq. (56) can be written
using the approximate relation $T_{\delta_{\rm b}}^{\rm hh}(k;>M_{\rm
th},>M^{\prime}_{\rm th})\approx T_{h}^{\rm hh}(k;>M_{\rm th},>M^{\prime}_{\rm
th})$ as

$\displaystyle\left.\frac{\mathrm{d}\mathscr{P}_{\rm
hh}}{\mathrm{d}\delta_{\rm b}}\right|_{\rm G}$ $\displaystyle=n(>M_{\rm
th})n(>M^{\prime}_{\rm th})\left.\frac{\mathrm{d}P_{\rm hh}(k;>M_{\rm
th},>M^{\prime}_{\rm th})}{\mathrm{d}\delta_{\rm b}}\right|_{\rm G}$
$\displaystyle=n(>M_{\rm th})n(>M^{\prime}_{\rm
th})\left(\frac{26}{21}\right)T_{\delta_{\rm b}}^{\rm hh}(k;>M_{\rm
th},>M^{\prime}_{\rm th})P_{\rm hh}(k;>M_{\rm th},>M^{\prime}_{\rm th})$
$\displaystyle\simeq n(>M_{\rm th})n(>M^{\prime}_{\rm
th})\left(\frac{26}{21}\right)T_{h}^{\rm hh}(k;>M_{\rm th},>M^{\prime}_{\rm
th})P_{\rm hh}(k;>M_{\rm th},>M^{\prime}_{\rm th}).$ (57)

The SU-growth response of $P_{\rm gg}^{2h}$ can be related to $T_{h}^{\rm
hh}(k;>M_{\rm th},>M^{\prime}_{\rm th})$ as

$\displaystyle\left.\frac{\partial P_{\rm gg}^{2h}(k)}{\partial\delta_{\rm
b}}\right|_{\rm G}$
$\displaystyle=\frac{1}{\bar{n}_{g}^{2}}\iint\mathrm{d}M\mathrm{d}M^{\prime}\frac{\mathrm{d}\tilde{U}_{g}(k;M^{\prime})}{\mathrm{d}M^{\prime}}\left[\frac{\mathrm{d}\langle\tilde{U}_{g}(k;M)\rangle}{\mathrm{d}M}\left.\frac{\partial\mathscr{P}_{\rm
hh}(k;>M,>M^{\prime})}{\partial\delta_{\rm b}}\right|_{{\rm
G},M,M^{\prime}}\right.$
$\displaystyle~{}~{}~{}~{}~{}\left.+\frac{2}{3}{\mathscr{P}_{\rm
hh}(k;>M,>M^{\prime})}\frac{\partial}{\partial\ln
k}\frac{\mathrm{d}\langle\tilde{U}_{g}(k;M)\rangle}{\mathrm{d}M}\right]-2b_{g}^{L}P_{\rm
gg}^{2h}(k)$
$\displaystyle\simeq\frac{1}{\bar{n}_{g}^{2}}\iint\mathrm{d}M\mathrm{d}M^{\prime}\left\\{\frac{\mathrm{d}\langle\tilde{U}_{g}(k;M)\rangle}{\mathrm{d}M}\frac{\mathrm{d}\tilde{U}_{g}(k;M^{\prime})}{\mathrm{d}M^{\prime}}\right.n(>M)n(>M^{\prime})\left[2b_{\mathrm{h},1}^{L}(>M^{\prime})P_{\rm
hh}(k;>M,M^{\prime})\right.$
$\displaystyle~{}~{}~{}~{}~{}+\left.\left.\left(\frac{26}{21}\right)T_{h}^{\rm
hh}(k;>M,>M^{\prime})P_{\rm hh}(k;>M,>M^{\prime})+\frac{2}{3}{P_{\rm
hh}(k;>M,>M^{\prime})}\frac{\partial}{\partial\ln
k}\frac{\mathrm{d}\ln\langle\tilde{U}_{g}(k;M)\rangle}{\mathrm{d}M}\right]\right\\}$
$\displaystyle~{}~{}~{}~{}~{}-2b_{g}^{L}P_{\rm gg}^{2h}(k).$ (58)

Using Eqs. (46) and (58), we can compute the SU growth response for the
galaxy-auto power spectrum as Eq. (26).

## References

  * Planck Collaboration  _et al._ (2020) Planck Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, _et al._ , Astronomy and Astrophysics 641, A6 (2020), arXiv:1807.06209 [astro-ph.CO] .
  * Alam  _et al._ (2017) S. Alam, M. Ata, S. Bailey, F. Beutler, D. Bizyaev, J. A. Blazek, _et al._ , Monthly Notices of the Royal Astronomical Society 470, 2617 (2017), arXiv:1607.03155 [astro-ph.CO] .
  * Kumar  _et al._ (2021) D. Kumar, D. Jain, S. Mahajan, A. Mukherjee, and N. Rani, Phys. Rev. D 103, 063511 (2021), arXiv:2002.06354 [astro-ph.CO] .
  * Vagnozzi  _et al._ (2021a) S. Vagnozzi, A. Loeb, and M. Moresco, The Astrophysical Journal 908, 84 (2021a), arXiv:2011.11645 [astro-ph.CO] .
  * Dhawan  _et al._ (2021) S. Dhawan, J. Alsing, and S. Vagnozzi, Monthly Notices of the Royal Astronomical Society 506, L1 (2021), arXiv:2104.02485 [astro-ph.CO] .
  * Alam  _et al._ (2021) S. Alam, M. Aubert, S. Avila, C. Balland, J. E. Bautista, M. A. Bershady, _et al._ , Phys. Rev. D 103, 083533 (2021), arXiv:2007.08991 [astro-ph.CO] .
  * Vagnozzi  _et al._ (2021b) S. Vagnozzi, E. Di Valentino, S. Gariazzo, A. Melchiorri, O. Mena, and J. Silk, Physics of the Dark Universe 33, 100851 (2021b), arXiv:2010.02230 [astro-ph.CO] .
  * Liu  _et al._ (2022) T. Liu, S. Cao, M. Biesiada, and S. Geng, The Astrophysical Journal 939, 37 (2022), arXiv:2204.07365 [astro-ph.CO] .
  * DES Collaboration  _et al._ (2022) DES Collaboration, T. M. C. Abbott, M. Aguena, A. Alarcon, O. Alves, A. Amon, _et al._ , arXiv e-prints , arXiv:2207.05766 (2022), arXiv:2207.05766 [astro-ph.CO] .
  * Tristram  _et al._ (2023) M. Tristram, A. J. Banday, M. Douspis, X. Garrido, K. M. Górski, S. Henrot-Versillé, S. Ilić, R. Keskitalo, G. Lagache, C. R. Lawrence, B. Partridge, and D. Scott, arXiv e-prints , arXiv:2309.10034 (2023), arXiv:2309.10034 [astro-ph.CO] .
  * Leonard  _et al._ (2016) C. D. Leonard, P. Bull, and R. Allison, Phys. Rev. D 94, 023502 (2016), arXiv:1604.01410 [astro-ph.CO] .
  * Peacock and Dodds (1996) J. A. Peacock and S. J. Dodds, Monthly Notices of the Royal Astronomical Society 280, L19 (1996), arXiv:astro-ph/9603031 [astro-ph] .
  * Smith  _et al._ (2003) R. E. Smith, J. A. Peacock, A. Jenkins, S. D. M. White, C. S. Frenk, F. R. Pearce, P. A. Thomas, G. Efstathiou, and H. M. P. Couchman, Monthly Notices of the Royal Astronomical Society 341, 1311 (2003), astro-ph/0207664 .
  * Heitmann  _et al._ (2010) K. Heitmann, M. White, C. Wagner, S. Habib, and D. Higdon, The Astrophysical Journal 715, 104 (2010), arXiv:0812.1052 .
  * Heitmann  _et al._ (2009) K. Heitmann, D. Higdon, M. White, S. Habib, B. J. Williams, E. Lawrence, and C. Wagner, The Astrophysical Journal 705, 156 (2009), arXiv:0902.0429 [astro-ph.CO] .
  * Takahashi  _et al._ (2012) R. Takahashi, M. Sato, T. Nishimichi, A. Taruya, and M. Oguri, The Astrophysical Journal 761, 152 (2012), arXiv:1208.2701 .
  * Heitmann  _et al._ (2016) K. Heitmann, D. Bingham, E. Lawrence, S. Bergner, S. Habib, D. Higdon, A. Pope, R. Biswas, H. Finkel, N. Frontiere, and S. Bhattacharya, The Astrophysical Journal 820, 108 (2016), arXiv:1508.02654 .
  * Lawrence  _et al._ (2017) E. Lawrence, K. Heitmann, J. Kwan, A. Upadhye, D. Bingham, S. Habib, D. Higdon, A. Pope, H. Finkel, and N. Frontiere, The Astrophysical Journal 847, 50 (2017), arXiv:1705.03388 .
  * Nishimichi  _et al._ (2019) T. Nishimichi, M. Takada, R. Takahashi, K. Osato, M. Shirasaki, T. Oogi, H. Miyatake, M. Oguri, R. Murata, Y. Kobayashi, and N. Yoshida, The Astrophysical Journal 884, 29 (2019), arXiv:1981.09504 [astro-ph.CO] .
  * Reimberg  _et al._ (2020) P. Reimberg, F. Bernardeau, T. Nishimichi, and M. Rizzato, Monthly Notices of the Royal Astronomical Society 492, 5226 (2020), arXiv:1811.02976 [astro-ph.CO] .
  * Mead  _et al._ (2021) A. J. Mead, S. Brieden, T. Tröster, and C. Heymans, Monthly Notices of the Royal Astronomical Society 502, 1401 (2021), arXiv:2009.01858 [astro-ph.CO] .
  * Smith and Angulo (2019) R. E. Smith and R. E. Angulo, Monthly Notices of the Royal Astronomical Society 486, 1448 (2019), arXiv:1807.00040 [astro-ph.CO] .
  * Euclid Collaboration  _et al._ (2021) Euclid Collaboration, M. Knabenhans, J. Stadel, D. Potter, J. Dakin, S. Hannestad, _et al._ , Monthly Notices of the Royal Astronomical Society 505, 2840 (2021), arXiv:2010.11288 [astro-ph.CO] .
  * Miyatake  _et al._ (2023) H. Miyatake, S. Sugiyama, M. Takada, T. Nishimichi, X. Li, M. Shirasaki, S. More, Y. Kobayashi, A. J. Nishizawa, M. M. Rau, T. Zhang, R. Takahashi, R. Dalal, R. Mandelbaum, M. A. Strauss, T. Hamana, M. Oguri, K. Osato, W. Luo, A. Kannawadi, B.-C. Hsieh, R. Armstrong, Y. Komiyama, R. H. Lupton, N. B. Lust, L. A. MacArthur, S. Miyazaki, H. Murayama, Y. Okura, P. A. Price, T. Sunayama, P. J. Tait, M. Tanaka, and S.-Y. Wang, arXiv e-prints , arXiv:2304.00704 (2023), arXiv:2304.00704 [astro-ph.CO] .
  * Sugiyama  _et al._ (2023) S. Sugiyama, H. Miyatake, S. More, X. Li, M. Shirasaki, M. Takada, Y. Kobayashi, R. Takahashi, T. Nishimichi, A. J. Nishizawa, M. M. Rau, T. Zhang, R. Dalal, R. Mandelbaum, M. A. Strauss, T. Hamana, M. Oguri, K. Osato, A. Kannawadi, R. Armstrong, Y. Komiyama, R. H. Lupton, N. B. Lust, S. Miyazaki, H. Murayama, Y. Okura, P. A. Price, P. J. Tait, M. Tanaka, and S.-Y. Wang, arXiv e-prints , arXiv:2304.00705 (2023), arXiv:2304.00705 [astro-ph.CO] .
  * Baldauf  _et al._ (2011) T. Baldauf, U. Seljak, L. Senatore, and M. Zaldarriaga, Journal of Cosmology and Astroparticle Physics 2011, 031 (2011), arXiv:1106.5507 [astro-ph.CO] .
  * Takada and Hu (2013) M. Takada and W. Hu, Phys. Rev. D 87, 123504 (2013), arXiv:1302.6994 [astro-ph.CO] .
  * Li  _et al._ (2014a) Y. Li, W. Hu, and M. Takada, Phys. Rev. D 89, 083519 (2014a), arXiv:1401.0385 [astro-ph.CO] .
  * Li  _et al._ (2014b) Y. Li, W. Hu, and M. Takada, Phys. Rev. D 90, 103530 (2014b), arXiv:1408.1081 [astro-ph.CO] .
  * Wagner  _et al._ (2015a) C. Wagner, F. Schmidt, C. T. Chiang, and E. Komatsu, Monthly Notices of the Royal Astronomical Society 448, L11 (2015a), arXiv:1409.6294 [astro-ph.CO] .
  * Terasawa  _et al._ (2022) R. Terasawa, R. Takahashi, T. Nishimichi, and M. Takada, Phys. Rev. D 106, 083504 (2022), arXiv:2205.10339 [astro-ph.CO] .
  * Schaan  _et al._ (2014) E. Schaan, M. Takada, and D. N. Spergel, Phys. Rev. D 90, 123523 (2014), arXiv:1406.3330 [astro-ph.CO] .
  * Miyatake  _et al._ (2022a) H. Miyatake, S. Sugiyama, M. Takada, T. Nishimichi, M. Shirasaki, Y. Kobayashi, R. Mandelbaum, S. More, M. Oguri, K. Osato, Y. Park, R. Takahashi, J. Coupon, C. Hikage, B.-C. Hsieh, Y. Komiyama, A. Leauthaud, X. Li, W. Luo, R. H. Lupton, S. Miyazaki, H. Murayama, A. J. Nishizawa, P. A. Price, M. Simet, J. S. Speagle, M. A. Strauss, M. Tanaka, and N. Yoshida, Phys. Rev. D 106, 083520 (2022a), arXiv:2111.02419 [astro-ph.CO] .
  * Baldauf  _et al._ (2016) T. Baldauf, U. Seljak, L. Senatore, and M. Zaldarriaga, Journal of Cosmology and Astroparticle Physics 2016, 007 (2016), arXiv:1511.01465 [astro-ph.CO] .
  * Krause and Eifler (2017) E. Krause and T. Eifler, Monthly Notices of the Royal Astronomical Society 470, 2100 (2017), arXiv:1601.05779 [astro-ph.CO] .
  * Bayer  _et al._ (2023) A. E. Bayer, J. Liu, R. Terasawa, A. Barreira, Y. Zhong, and Y. Feng, Phys. Rev. D 108, 043521 (2023), arXiv:2210.15647 [astro-ph.CO] .
  * Chisari  _et al._ (2019) N. E. Chisari, D. Alonso, E. Krause, C. D. Leonard, P. Bull, J. Neveu, A. S. Villarreal, S. Singh, T. McClintock, J. Ellison, Z. Du, J. Zuntz, A. Mead, S. Joudaki, C. S. Lorenz, T. Tröster, J. Sanchez, F. Lanusse, M. Ishak, R. Hlozek, J. Blazek, J.-E. Campagne, H. Almoubayyed, T. Eifler, M. Kirby, D. Kirkby, S. Plaszczynski, A. Slosar, M. Vrastil, E. L. Wagoner, and LSST Dark Energy Science Collaboration, The Astrophysical Journal Supplement Series 242, 2 (2019), arXiv:1812.05995 [astro-ph.CO] .
  * Bayer  _et al._ (2021) A. E. Bayer, A. Banerjee, and U. Seljak, arXiv e-prints , arXiv:2108.04215 (2021), arXiv:2108.04215 [astro-ph.CO] .
  * Wagner  _et al._ (2015b) C. Wagner, F. Schmidt, C.-T. Chiang, and E. Komatsu, Journal of Cosmology and Astroparticle Physics 2015, 042 (2015b), arXiv:1503.03487 [astro-ph.CO] .
  * Takahashi  _et al._ (2019) R. Takahashi, T. Nishimichi, M. Takada, M. Shirasaki, and K. Shiroyama, Monthly Notices of the Royal Astronomical Society 482, 4253 (2019), arXiv:1805.11629 [astro-ph.CO] .
  * Zheng  _et al._ (2005) Z. Zheng, A. A. Berlind, D. H. Weinberg, A. J. Benson, C. M. Baugh, S. Cole, R. Davé, C. S. Frenk, N. Katz, and C. G. Lacey, The Astrophysical Journal 633, 791 (2005), astro-ph/0408564 .
  * Tinker  _et al._ (2010) J. L. Tinker, B. E. Robertson, A. V. Kravtsov, A. Klypin, M. S. Warren, G. Yepes, and S. Gottlöber, The Astrophysical Journal 724, 878 (2010), arXiv:1001.3162 .
  * Hoffmann  _et al._ (2017) K. Hoffmann, J. Bel, and E. Gaztañaga, Monthly Notices of the Royal Astronomical Society 465, 2225 (2017), arXiv:1607.01024 [astro-ph.CO] .
  * Sheth and Tormen (2002) R. K. Sheth and G. Tormen, Monthly Notices of the Royal Astronomical Society 329, 61 (2002), arXiv:astro-ph/0105113 .
  * Chiang  _et al._ (2014) C.-T. Chiang, C. Wagner, F. Schmidt, and E. Komatsu, Journal of Cosmology and Astroparticle Physics 2014, 048 (2014), arXiv:1403.3411 [astro-ph.CO] .
  * Angulo and Pontzen (2016) R. E. Angulo and A. Pontzen, Monthly Notices of the Royal Astronomical Society 462, L1 (2016), arXiv:1603.05253 [astro-ph.CO] .
  * Springel (2005) V. Springel, Monthly Notices of the Royal Astronomical Society 364, 1105 (2005).
  * Scoccimarro (1998) R. Scoccimarro, Monthly Notices of the Royal Astronomical Society 299, 1097 (1998), arXiv:astro-ph/9711187 .
  * Crocce and Scoccimarro (2006) M. Crocce and R. Scoccimarro, Phys. Rev. D 73, 063519 (2006), arXiv:astro-ph/0509418 .
  * Nishimichi  _et al._ (2009) T. Nishimichi, A. Shirata, A. Taruya, K. Yahata, S. Saito, Y. Suto, R. Takahashi, N. Yoshida, T. Matsubara, N. Sugiyama, I. Kayo, Y. Jing, and K. Yoshikawa, Publ. Astron. Soc. Japan 61, 321 (2009), arXiv:0810.0813 .
  * Valageas and Nishimichi (2011) P. Valageas and T. Nishimichi, Astronomy and Astrophysics 527, A87 (2011), arXiv:1009.0597 [astro-ph.CO] .
  * Lewis  _et al._ (2000) A. Lewis, A. Challinor, and A. Lasenby, The Astrophysical Journal 538, 473 (2000).
  * Planck Collaboration  _et al._ (2016) Planck Collaboration, P. A. R. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro, J. G. Bartlett, and et al., Astronomy and Astrophysics 594, A13 (2016), arXiv:1502.01589 .
  * Villaescusa-Navarro  _et al._ (2018) F. Villaescusa-Navarro, S. Naess, S. Genel, A. Pontzen, B. Wandelt, L. Anderson, A. Font-Ribera, N. Battaglia, and D. N. Spergel, The Astrophysical Journal 867, 137 (2018), arXiv:1806.01871 [astro-ph.CO] .
  * Behroozi  _et al._ (2013) P. S. Behroozi, R. H. Wechsler, and H.-Y. Wu, The Astrophysical Journal 762, 109 (2013), arXiv:1110.4372 [astro-ph.CO] .
  * Li  _et al._ (2016) Y. Li, W. Hu, and M. Takada, Phys. Rev. D 93, 063507 (2016), arXiv:1511.01454 .
  * Navarro  _et al._ (1996) J. F. Navarro, C. S. Frenk, and S. D. M. White, The Astrophysical Journal 462, 563 (1996), astro-ph/9508025 .
  * Diemer and Kravtsov (2015) B. Diemer and A. V. Kravtsov, The Astrophysical Journal 799, 108 (2015), arXiv:1407.4730 .
  * Miyatake  _et al._ (2022b) H. Miyatake, Y. Kobayashi, M. Takada, T. Nishimichi, M. Shirasaki, S. Sugiyama, R. Takahashi, K. Osato, S. More, and Y. Park, Phys. Rev. D 106, 083519 (2022b), arXiv:2101.00113 [astro-ph.CO] .
  * Hockney and Eastwood (1981) R. W. Hockney and J. W. Eastwood,  _Computer Simulation Using Particles, New York: McGraw-Hill, 1981_ (New York: McGraw-Hill, 1981).
  * Jing (2005) Y. P. Jing, The Astrophysical Journal 620, 559 (2005), arXiv:astro-ph/0409240 .
  * Tinker  _et al._ (2008) J. Tinker, A. V. Kravtsov, A. Klypin, K. Abazajian, M. Warren, G. Yepes, S. Gottlöber, and D. E. Holz, The Astrophysical Journal 688, 709 (2008), arXiv:0803.2706 .
  * Voivodic and Barreira (2021) R. Voivodic and A. Barreira, Journal of Cosmology and Astroparticle Physics 2021, 069 (2021), arXiv:2012.04637 [astro-ph.CO] .
  * Nicola  _et al._ (2020) A. Nicola, D. Alonso, J. Sánchez, A. Slosar, H. Awan, A. Broussard, J. Dunkley, E. Gawiser, Z. Gomes, R. Mandelbaum, H. Miyatake, J. A. Newman, I. Sevilla-Noarbe, S. Skinner, and E. L. Wagoner, Journal of Cosmology and Astroparticle Physics 2020, 044 (2020), arXiv:1912.08209 [astro-ph.CO] .