# Toward a more stringent test of gravity with redshift space power spectrum:
simultaneous probe of growth and amplitude of large-scale structure

Yong-Seon Song1, Yi Zheng2, Atsushi Taruya3,4<EMAIL_ADDRESS>1 Korea
Astronomy and Space Science Institute, Daejeon 34055, Republic of Korea
2 School of Physics and Astronomy, Sun Yat-sen University, 2 Daxue Road,
Tangjia, Zhuhai, 519082, China
3 Center for Gravitational Physics, Yukawa Institute for Theoretical Physics,
Kyoto University, Kyoto 606-8502, Japan,
4 Kavli Institute for the Physics and Mathematics of the Universe, Todai
Institutes for Advanced Study, the University of Tokyo, Kashiwa, Chiba
277-8583, Japan (Kavli IPMU, WPI)

###### Abstract

Redshift-space distortions (RSD) offers an exciting opportunity to test the
gravity on cosmological scales. In the presence of galaxy bias, however, the
RSD measurement at large scales, where the linear theory prediction is safely
applied, is known to exhibit a degeneracy between the parameters of structure
growth $f$ and fluctuation amplitude $\sigma_{8}$, and one can only constrain
the parameters in the form of $f\sigma_{8}$. In order to disentangle this
degeneracy, in this paper, we go beyond the linear theory, and consider the
model of RSD applicable to a weakly nonlinear regime. Based on the Fisher
matrix analysis, we show explicitly that the degeneracy of the parameter
$f\sigma_{8}$ can be broken, and $\sigma_{8}$ is separately estimated in the
presence of galaxy bias. Performing further the Markov chain Monte Carlo
analysis, we verify that our model correctly reproduces the fiducial values of
$f\sigma_{8}$ and $\sigma_{8}$, with the statistical errors consistent with
those estimated from the Fisher matrix analysis. We show that upcoming galaxy
survey of the stage-IV class can unambiguously determine $\sigma_{8}$ at the
precision down to $\lesssim 10$% at higher redshifts even if we restrict the
accessible scales to $k\lesssim 0.16\,h$ Mpc-1.

Large-scale structure formation

###### pacs:

98.80.-k;04.50.Kd;98.65.Dx

††preprint: YITP-21-09

##  I Introduction

Since its discovery two decades ago Riess et al. (1998); Perlmutter et al.
(1999), the origin of late-time cosmic acceleration has remained puzzled.
While the flat Lambda cold dark matter ($\Lambda$CDM) model is currently the
best as concordant cosmological model to describe both the cosmic expansion
and structure formation, consistent with observations of cosmic microwave
background and large-scale structure, the tension with the cosmological
parameters determined at the local universe has been recently highlighted,
suggesting a need of new physics beyond $\Lambda$CDM model Riess et al.
(2016); Wong et al. (2020); Jedamzik et al. (2020).

Theoretically, the origin of cosmic acceleration can be explained by either
the presence of a mysterious energy component called dark energy or a long-
distance modification of gravity, referred to as modified gravity Lue (2006);
Frieman et al. (2008); Li et al. (2011); Clifton et al. (2012); Joyce et al.
(2016); Koyama (2016); Nojiri et al. (2017); Arun et al. (2017); Wang et al.
(2017); Brax (2018); Ishak (2019). In order to realize the accelerated cosmic
expansion, the former introduces a negative pressure support, and the latter
changes a law of gravitational physics on large scales. Observationally
discriminating between two scenarios therefore requires a simultaneous
measurements of the cosmic expansion and growth of structure.

Among various cosmological probes, the redshift-space galaxy clustering offers
a sensible probe of both the cosmic expansion and growth of structure. At
large scales, the baryon acoustic oscillations (BAO) imprinted on the
clustering pattern of galaxies appears statistically manifest, and it can be
used for a standard ruler to determine the angular diameter distance ($D_{A}$)
and Hubble parameters $H$ at high redshifts through the Alcock-Paczynski test
Alcock and Paczynski (1979). Further, the observed galaxy distribution via
spectroscopic surveys is statistically anisotropic due to the peculiar
velocity of galaxies by Doppler effect, referred to as the redshift-space
distortions (RSD). In linear theory, the strength of anisotropies is solely
characterized by the growth rate $f$, defined by $f=d\ln D_{+}/d\ln a$, with
$D_{+}$ and $a$ being linear growth factor and scale factor of the Universe,
respectively. Thus, through a precision measurement of redshift-space galaxy
power spectrum or correlation function, one can in principle obtain
simultaneously the information on the three parameters (i.e., $D_{A}$, $H$,
and $f$).

However, in the presence of the galaxy bias, the situation becomes bit
complicated. To be precise, consider the large scales where the linear theory
is safely applied. Then, the galaxy power spectrum is generally modeled as

$P^{(S)}_{g}(k,\mu)=\left[b_{1}\sigma_{8}(z)+f\sigma_{8}(z)\mu^{2}\right]^{2}P_{\rm
m}(k)\,,$ (1)

where the variable $\mu$ is the directional cosine between wavevector and
line-of-sight direction. The function $P_{\rm m}$ is the matter power spectrum
normalized at the present day. The quantity $\sigma_{8}(z)$ is the root-mean-
square mass fluctuation in spheres with radius $8\,h^{-1}\,{\rm Mpc}$ at
redshift $z$, and is recast in linear theory as
$\sigma_{8}(z)=\sigma_{8}(0)D_{+}(z)$ with $D_{+}$ being normalized to unity
at $z=0$. The $b_{1}$ is the linear bias parameter. Note that taking further
the Alcock-Paczynski effect into account, the projected wavenumbers
perpendicular and parallel to the line-of-sight direction, $k_{\perp}$ and
$k_{\parallel}=\mu\,k$, are respectively replaced with $D_{A}/D_{A,{\rm
fid}}k_{\perp}$ and $(H/H_{\rm fid})^{-1}k_{\parallel}$, and the power
spectrum given above is further multiplied by $(H/H_{\rm fid})(D_{A}/D_{A,{\rm
fid}})^{-2}$, where the quantities with subscript indicate those estimated in
a fiducial cosmological model.

The structure of Eq. (1) indicates that the parameter $\sigma_{8}$ is
degenerated with growth rate $f$ and bias $b_{1}$, and one can only determine
the combinations of the parameters, i.e., $b_{1}\sigma_{8}$ and $f\sigma_{8}$,
through the observed power spectrum. In other words, unless the bias parameter
$b_{1}$ is known a priori, one cannot break the degeneracy between growth rate
$f$ and $\sigma_{8}$. Note that the Alcock-Paczynski effect induces
distinctive anisotropies in the measured power spectrum, and making use of BAO
features, one can separately determine $D_{A}/D_{A,{\rm fid}}$ and $H/H_{\rm
fid}$ without any degeneracy with $\sigma_{8}$.

Toward a solid test of gravity and cosmic acceleration, the degeneracy between
$f$ and $\sigma_{8}$ has to be broken. To do this, one simple approach is to
combine the power spectrum with other cosmological probe. In this respect, the
use of galaxy bispectrum would be obviously important, and this can also
provide an additional cosmological information, further tightening the
cosmological constraints Song et al. (2015); Gil-Marín et al. (2015). Another
approach is to stick to the power spectrum, and to use the small-scale
information beyond the linear regime. Recalling that Eq. (1) is valid only at
large scales, if we go to nonlinear regime, there appear corrections involving
the parameters $\sigma_{8}$ and $f$ but with a different combination. In fact,
using the perturbation theory calculation, we can identify the directional-
dependent terms proportional to $f^{2}\sigma_{8}^{4}$ or $f^{3}\sigma_{8}^{4}$
in the matter power spectrum Taruya et al. (2010, 2013). Thus, provided an
accurate theoretical model, accessing (weakly) nonlinear scales would give a
way to break the degeneracy between $f$ and $\sigma_{8}$. A price to pay is,
however, the new degrees of freedom to characterize the galaxy bias at
nonlinear scales. That is, beyond linear scales, we need to introduce several
bias parameters describing the nonlinear modification to power spectrum, and
one has to marginalize them to determine parameters sensitive to the
cosmology. It is thus not trivial at all that the growth rate $f$ is uniquely
determined without any degeneracy.

In this paper, based on a nonlinear theoretical model of the redshift-space
galaxy power spectrum, we explicitly demonstrate that the degeneracy between
$f$ and $\sigma_{8}$ is broken at weakly nonlinear scales. The model
specifically considered here is called the hybrid RSD model that has been
developed by Ref. Zheng and Song (2016); Song et al. (2018); Zheng et al.
(2019). This is a perturbation theory based model Taruya et al. (2010, 2013),
but partly including the correction terms calibrated by $N$-body simulations,
the accuracy is improved and the applicable range is extended Song et al.
(2018). In this paper, incorporating further the nonlinear bias prescription
into the hybrid RSD model, we will investigate how the degeneracy between $f$
and $\sigma_{8}$ is broken, and quantify the expected constraints on
parameters $\sigma_{8}$, $f\sigma_{8}$ as well as $D_{A}$ and $H$ for a
representative galaxy survey (Dark Energy Spectroscopy Instrument, DESI).

Related to the present study, one may comment on the full-shape analysis that
recently performed using SDSS BOSS galaxies Ivanov et al. (2020); D’Amico et
al. (2020). In this study, based on the effective field theory of large-scale
structure, the input linear power spectrum is allowed to vary, enabling us to
directly constrain each cosmological parameters. Going beyond linear regime,
not only tightening constraints but also breaking parameter degeneracy is
shown to be manifest Nishimichi et al. (2020), assuming the underlying theory
of gravity. The spirit of this approach is close to the present paper, but we
here stick to the consistent test of gravity, and take $f$ or $f\sigma_{8}$ to
be a free parameter, independent of other cosmological parameters.

This paper is organized as follows. Sec. II presents a model of the redshift-
space galaxy power spectrum applicable to the weakly nonlinear regime. Based
on this model as a theoretical template, in Sec. III, we use the Fisher matrix
formalism to quantitatively investigate how well one can break the degeneracy
of the parameters $f\sigma_{8}$. Sec. IV examines the Markov chain Monte Carlo
analysis to estimate the cosmological parameters in the hal catalog. We verify
that our model of RSD faithfully reproduces the fiducial parameters in the
$N$-body simulations, consistent with the Fisher forecast results. Finally,
Sec. V is devoted to conclusion and discussion.

##  II Hybrid RSD model for galaxy clustering

In this section, we present the model of redshift-space galaxy power spectrum
beyond the linear regime. After briefly reviewing the hybrid RSD model for
matter fluctuations in Sec. II.1, we extend it to incorporate the nonlinear
galaxy bias in Sec. II.2.

###  II.1 Hybrid RSD model for matter power spectrum

In order to model the observed galaxy power spectrum, one critical ingredient
is the redshift-space distortions (RSD). In principle, all the effects of RSD
is accounted for by the simple relation between real-space ($r$) and redshift-
space ($s$) positions:

$\mbox{\boldmath$s$}=\mbox{\boldmath$r$}+\frac{\mbox{\boldmath$v$}\cdot\hat{z}}{aH}\hat{z},$
(2)

where the quantities $v$, $a$ and $H$ respectively denote the physical
peculiar velocity, the scale factor of the Universe, and the Hubble parameter.
In this paper, we will take the distant-observer limit, and choose the
$z$-direction as the line-of-sight direction. With the mapping relation given
above, the power spectrum in redshift space is generally expressed in terms of
the real-space quantities (e.g., Ref. Taruya et al. (2010)):

$P^{\rm(S)}(k,\mu)=\int
d^{3}\mbox{\boldmath$x$}\,e^{i\,\mbox{\boldmath$k$}\cdot\mbox{\boldmath$x$}}\bigl{\langle}e^{j_{1}A_{1}}A_{2}A_{3}\bigr{\rangle}\,,$
(3)

where the variable $j_{1}$ and functions $A_{i}$ are defined as follows:

$\displaystyle j_{1}=-i\,k\mu,$ $\displaystyle
A_{1}=u_{z}(\mbox{\boldmath$r$})-u_{z}(\mbox{\boldmath$r$}^{\prime}),$
$\displaystyle
A_{2}=\delta(\mbox{\boldmath$r$})+\,\nabla_{z}u_{z}(\mbox{\boldmath$r$}),$
$\displaystyle
A_{3}=\delta(\mbox{\boldmath$r$}^{\prime})+\,\nabla_{z}u_{z}(\mbox{\boldmath$r$}^{\prime}).$

The quantities $x$ and $u$ are given respectively by
$\mbox{\boldmath$x$}=\mbox{\boldmath$r$}-\mbox{\boldmath$r$}^{\prime}$ and
$\mbox{\boldmath$u$}\equiv-\mbox{\boldmath$v$}/(aH)$. The function $u_{z}$ is
the line-of-sight component of $u$.

Based on Eq. (3), a rigorous calculation of the redshift-space power spectrum
generally requires a non-perturbative treatment. This is true even if the
density and velocity follow the Gaussian statistics. We thus employ the
perturbative treatment and derive the expression relevant to the weakly
nonlinear regime. To do this, one proposition made in Refs. Taruya et al.
(2010); Zheng et al. (2015) is that a part of the zero-lag correlation in the
exponent is kept as a non-perturbative contribution, while rest of the terms
is Taylor-expanded. The resultant expression for the power spectrum, relevant
at the next-to-next-to-leading order becomes Taruya et al. (2010, 2013); Zheng
et al. (2015)

$\displaystyle P^{\rm(S)}(k,\mu)$ $\displaystyle=D^{\rm
FoG}(k\mu\sigma_{z})\,\Bigl{[}P_{\delta\delta}(k)+2\mu^{2}P_{\delta\Theta}(k)+\mu^{4}P_{\Theta\Theta}(k)$
$\displaystyle\qquad+A(k,\mu)+B(k,\mu)+T(k,\mu)+F(k,\mu)\Bigr{]},$

where the quantity $\Theta$ is the velocity-divergence field,
$\Theta\equiv-\nabla\cdot\mbox{\boldmath$v$}/(aH)=\nabla\cdot\mbox{\boldmath$u$}$.
The spectra, $P_{\delta\delta}$, $P_{\delta\Theta}$, and $P_{\Theta\Theta}$
are respectively the auto-power spectrum of density, velocity-divergence
fields, and their cross-power spectrum. The first line in the bracket is
originated from the term $\langle A_{2}A_{3}\rangle_{c}$, and is obtained
assuming the irrotational flow. At the linear order, it is reduced to the
squashing Kaiser term, i.e., Eq. (1). The rest of the terms in the bracket are
the higher-order corrections characterizing the nonlinear correlations between
the density and velocity fields, defined by

$\displaystyle A(k,\mu)$ $\displaystyle=$ $\displaystyle j_{1}\,\int
d^{3}\mbox{\boldmath$x$}\,\,e^{i\mbox{\boldmath$k$}\cdot\mbox{\boldmath$x$}}\,\,\langle
A_{1}A_{2}A_{3}\rangle_{c},$ $\displaystyle B(k,\mu)$ $\displaystyle=$
$\displaystyle j_{1}^{2}\,\int
d^{3}\mbox{\boldmath$x$}\,\,e^{i\mbox{\boldmath$k$}\cdot\mbox{\boldmath$x$}}\,\,\langle
A_{1}A_{2}\rangle_{c}\,\langle A_{1}A_{3}\rangle_{c},$ $\displaystyle
T(k,\mu)$ $\displaystyle=$ $\displaystyle\frac{1}{2}j_{1}^{2}\,\int
d^{3}\mbox{\boldmath$x$}\,\,e^{i\mbox{\boldmath$k$}\cdot\mbox{\boldmath$x$}}\,\,\langle
A_{1}^{2}A_{2}A_{3}\rangle_{c},$ $\displaystyle F(k,\mu)$ $\displaystyle=$
$\displaystyle-j_{1}^{2}\,\int
d^{3}\mbox{\boldmath$x$}\,\,e^{i\mbox{\boldmath$k$}\cdot\mbox{\boldmath$x$}}\,\,\langle
u_{z}u_{z}^{\prime}\rangle_{c}\langle A_{2}A_{3}\rangle_{c}.$

Note that the factorized term $D^{\rm FoG}$ in Eq. (II.1) represents a non-
perturbative contribution coming from the zero-lag correlation of velocity
fields, and it plays a role to suppress the overall amplitude at small scales.
The explicit functional form will be later specified (see Sec. III).

To explicitly compute Eq. (II.1) in the case of the matter fluctuations, in
Ref. Song et al. (2018), we have used the re-summed perturbation theory
treatment by the multi-point propagator expansion Taruya et al. (2012, 2013);
Bernardeau et al. (2008). To account further for a small but non-negligible
flaw in the perturbative calculations, we have added the corrections
calibrated by the $N$-body simulations, which enables us to predict the
redshift-space matter power spectrum at the precision of $\lesssim 1$% down to
$k\simeq 0.18\,h\,{\rm Mpc}^{-1}$. In Appendix A, the description of our
hybrid model is presented for the matter power spectrum, and the explicit
dependence on the quantities $\sigma_{8}(z)$ and $f(z)$ are shown.

Although the expressions of the higher-order contributions are rather
intricate, these terms involve the contributions expressed in the polynomial
form of $f^{m}\sigma^{n}$, with the power-law indices $m$ and $n$ running
respectively from $1$ to $4$ and $4$ to $6$ for the perturbative calculations
at next-to-next-to-leading order. Besides, beyond the linear order, even the
auto-power spectra $P_{\delta\delta}$ and $P_{\Theta\Theta}$ do not simply
scale as $\sigma_{8}^{2}$ and $(f\,\sigma_{8})^{2}$, respectively. Thus, if
one gets access to the weakly nonlinear regime where the higher-order terms
play a role, we expect that the degeneracy of the parameter $f\sigma_{8}$ at
Eq. (1) is broken, and $f$ and $\sigma_{8}$ can be separately determined.

Figure 1: Left: Logarithmic derivative of the real-space halo power spectrum
with respect to the linear growth factor $G_{\delta}$ (black) and linear bias
parameter $b_{1}$ (blue), plotted as function of the wavenumber $k$. In each
panel, the solid lines represent the results for hybrid model,
$P_{\delta\delta}^{\rm hybrid}$, while dashed and dotted lines are for the
linear theory prediction, $P_{\delta\delta}^{\rm lin}$. Right: logarithmic
derivative of the redshift-space halo power spectra with respect to
$G_{\delta}$, $b_{1}$ and $G_{\Theta}$, respectively shown from top to bottom
panels. Both of the Kaiser (linear theory) and hybrid models are shown as blue
and red curves, specifically fixing the directional cosine to $\mu=0.0$
(dotted), $0.3$ (short-dashed), $0.6$ (long-dashed) and $0.9$ (solid).

###  II.2 Modeling galaxy bias

On top of the hybrid RSD model in previous subsection, there is one more step
toward a practical application to the observed galaxy power spectrum. Since
the galaxy distribution is a biased tracer of matter distribution, accounting
for the galaxy bias is another crucial task. In order to do this, one may
recall that the expression given at Eq. (II.1) is fairly general. Replacing
the matter density field $\delta$ with the galaxy/halo density field
$\delta_{h}$, Eq. (II.1) can be applied to the observed power spectrum, and
hence our hybrid model of matter power spectrum is used as a building block to
compute accurately the redshift-space galaxy power spectrum.

For our interest at the weakly nonlinear regime, a perturbative description of
the tracer field is valid, and $\delta_{h}$ is expanded in powers of the
matter density field $\delta$, including the non-local contributions. Here, we
adopt the prescription proposed by Ref. McDonald and Roy (2009); Saito et al.
(2014), which has been applied to the SDSS BOSS (e.g., Gil-Marín et al.
(2015)) and eBOSS galaxies de Mattia et al. (2020). Apart from the stochastic
terms, this is a general perturbative expansion valid at the next-to-leading
order. While the hybrid model for the matter power spectrum in Sec. II.1 and
Appendix A includes the corrections valid at next-to-next-to-leading order, we
shall below consider the scales where the next-to-leading order is important,
but the next-to-next-to-leading order is still subdominant. A fully consistent
treatment including the bias at next-to-next-to-leading order will be
discussed in our future work.

Based on Refs. McDonald and Roy (2009); Saito et al. (2014), the auto-power
spectrum of the tracer density field,
$P_{\tilde{\delta}_{h}\tilde{\delta}_{h}}$, is expressed as follows:

$\displaystyle P_{\tilde{\delta}_{h}\tilde{\delta}_{h}}(k)$ $\displaystyle=$
$\displaystyle P_{\delta_{h}\delta_{h}}(k)+P_{\epsilon\epsilon}\,,$ (4)

where the term $P_{\epsilon\epsilon}$ is the stochastic contribution mainly
characterizing the shot noise, which is usually constant, and is estimated,
assuming the Poisson noise, from the number density of tracer field. The first
term, $P_{\delta_{h}\delta_{h}}$, is the deterministic part, and is given in
the following form:

$\displaystyle P_{\delta_{h}\delta_{h}}(k)$
$\displaystyle=b_{1}^{2}P_{\delta\delta}(k)+2b_{1}b_{2}P_{b2,\delta}(k)+2b_{2}b_{s2}P_{bs2,\delta}(k)$
$\displaystyle+2b_{1}b_{3\rm{nl}}\sigma_{3}^{2}(k)P^{\rm{L}}_{\rm
m}(k)+b_{2}^{2}P_{b22}(k)$
$\displaystyle+2b_{2}b_{s2}P_{b2s2}(k)+b_{s2}^{2}P_{b22}(k).$ (5)

On the other hand, the cross-power spectrum, $P_{\delta_{h}\theta_{h}}$,
contains only the deterministic contribution, whose expression is given by

$\displaystyle P_{\delta_{h}\Theta_{h}}(k)$
$\displaystyle=b_{1}P_{\delta\Theta}(k)+b_{2}P_{b2,\Theta}(k)+b_{s2}P_{bs2,\Theta}(k)$
$\displaystyle+b_{3\rm{nl}}\sigma_{3}^{2}(k)P^{\rm{L}}_{\rm m}(k)\,.$ (6)

Note that in the absence of velocity bias, the auto-power spectrum of
galaxy/halo velocity field, $P_{\Theta_{h}\Theta_{h}}$, is identical to
$P_{\Theta\Theta}$. In Eqs. (4) and (6), while the first terms,
$b_{1}^{2}P_{\delta\delta}$ and $b_{1}P_{\delta\Theta}$, as well as
$P_{\Theta_{h}\Theta_{h}}$, are the leading-order bias contribution, and can
be computed with the hybrid RSD model for the matter power spectrum, the terms
involving parameters $b_{2}$, $b_{s2}$, and $b_{\rm 3nl}$ are the higher-order
bias terms. The explicit expressions for their scale-dependent functions,
$P_{b2,\delta}$, $P_{bs2,\delta}$, $\sigma_{3}^{2}$, $P_{b22}$, …, can be
found in Ref. McDonald and Roy (2009); Saito et al. (2014).

Regarding the $A$, $B$, $T$ and $F$ terms in Eq. (II.1), they are all regarded
as the higher-order corrections. We only apply the linear bias $b_{1}$ to
these terms Zheng et al. (2019),

$\displaystyle A_{h}(k,\mu)$ $\displaystyle=$ $\displaystyle
b_{1}^{3}A(k,\mu,f/b_{1})\,,$ (7) $\displaystyle B_{h}(k,\mu)$
$\displaystyle=$ $\displaystyle b_{1}^{4}B(k,\mu,f/b_{1})\,,$ (8)
$\displaystyle F_{h}(k,\mu)$ $\displaystyle=$ $\displaystyle
b_{1}^{4}F(k,\mu,f/b_{1})\,,$ (9) $\displaystyle T_{h}(k,\mu)$
$\displaystyle=$ $\displaystyle b_{1}^{4}T(k,\mu,f/b_{1})\,,$ (10)

where all the subscripts $h$ at left-hand-side represent the higher-order
terms of the tracer fields. The accuracy of this treatment and its impacts on
the RSD model was discussed in Ref. Zheng et al. (2019).

Figure 2: Left: the expected errors on $f\sigma_{8}$ for DESI, with
$\sigma_{8}$ marginalized over. The result derived from the Fisher matrix,
depicted as the error bar around the fiducial value, are plotted against
redshift. The upper and lower panels respectively show the estimated
constraints on $f\sigma_{8}$ derived from the LRG and ELG type galaxies.
Right: Same as in the left panel, but the parameter $\sigma_{8}$ is fixed, and
is not treated as a free parameter.

Incorporating the bias prescription given above into Eq. (II.1), our hybrid
RSD model in Sec. II.1 and Appendix A enables us to predict the redshift-space
galaxy/halo power spectrum beyond the linear regime. To see the behavior of
the model beyond the linear bias prescription, left panel of Fig. 1 shows the
logarithmic response of the real-space power spectrum to the quantities
$b_{1}$ and $G_{\delta}$, depicted as respectively blue and black lines. In
linear theory, there is no difference between these responses, which exactly
give $2$ (dashed and dotted). On the other hand, the hybrid model involving
the galaxy bias expansion exhibits a wiggle feature for the response to
$G_{\delta}$, showing a small but visible deviation from $2$. This is caused
by non–linear smoothing around BAO peaks. On the other hand, the response to
$b_{1}$ still lies at $2$ even in the hybrid model, indicating that the
dependence of $G_{\delta}$ or $\sigma_{8}$ and $b_{1}$ becomes distinguishable
beyond the linear regime. The selected fiducial values of galaxy biases are
presented in Table 1. The fiducial values for coherent bias $b_{1}$ are taken
from DESI prediction DESI Collaboration et al. (2016), and the redshift
dependence of $b_{2}$ biases are fitting results from halo catalogue.

To feel some flavors on how the degeneracy is broken, left panel of Fig. 1
presents the response of the redshift-space power spectrum to the variation of
$G_{\delta}$ and $G_{\Theta}$. The quantity $G_{\delta}$ is the linear growth
factor normalized at primordial epoch as $G_{\delta}(z_{i})(1+z_{i})=1$, and
is related to $\sigma_{8}(z)$ through
$\sigma_{8}(z)=[G_{\delta}(z)/G_{\delta}(0)]\,\sigma_{8}(0)$. In left panel of
Fig. 1, the response of the redshift-space galaxy/halo power spectrum to the
variation of $G_{\delta}$ is shown for both linear theory prediction and
hybrid model, depicted as black dotted and solid curves respectively. The
function $G_{\Theta}$ is related to the linear growth rate $f$ through
$G_{\Theta}=f\,G_{\delta}$. On the other hand, right panel of Fig. 1 plots the
response of the 2D redshift-space power spectrum to the variations of
$G_{\delta}$ and $G_{\Theta}$ at top and bottom panels respectively,
specifically fixing the directional cosine $\mu$ to $0.0$, $0.3$, $0.6$ and
$0.9$. The distinct behavior of the power spectrum is observed by varying
$\mu$, and the results of hybrid model start to deviate from linear theory
predictions, indicating that the parameter degeneracy between $G_{\delta}$ and
$G_{\Theta}$ is broken, and so is the case for $\sigma_{8}$ and $f\sigma_{8}$.

##  III Forecast constraints on $\sigma_{8}$ and $f\sigma_{8}$

$z$ | $n_{g}^{\rm LRG}$ | $n_{g}^{\rm ELG}$ | $V$ | $b_{1}^{\rm LRG}$ | $b_{1}^{\rm ELG}$ | $b_{2}^{\rm LRG}$ | $b_{2}^{\rm ELG}$
---|---|---|---|---|---|---|---
0.4–0.6 | $4.9\times 10^{-4}$ | $1.6\times 10^{-4}$ | 3.5 | 2.22 | 1.10 | 0.10 | 0.10
0.6–0.8 | $9.9\times 10^{-4}$ | $1.4\times 10^{-3}$ | 5.4 | 2.45 | 1.21 | 0.67 | 0.67
0.8–1.0 | $3.9\times 10^{-4}$ | $1.7\times 10^{-3}$ | 7.0 | 2.69 | 1.33 | 1.40 | 1.40
1.0–1.2 | $2.4\times 10^{-5}$ | $9.8\times 10^{-4}$ | 8.4 | 2.94 | 1.45 | 2.48 | 2.48
1.2–1.4 | — | $5.8\times 10^{-4}$ | 9.4 | — | 1.57 | — | 3.91
1.4–1.6 | — | $2.3\times 10^{-4}$ | 10. | — | 1.70 | — | 5.68
Table 1:  The expected number density of LRG and ELG galaxies $n_{g}^{\rm
LRG}[h^{3}{\rm Mpc}^{-3}]$ and $n_{g}^{\rm ELG}[h^{3}{\rm Mpc}^{-3}]$, and
survey volume $V[h^{-3}\,{\rm Gpc}^{3}]$ at each redshift bin used in the
Fisher matrix analysis. These specific values are taken from those assumed in
the DESI experiment DESI Collaboration et al. (2016). The fiducial values of
bias for LRG and ELG are presented as well, calculated by the formulas
$b_{1}^{\rm LRG}(z)D(z)=1.7$ and $b_{1}^{\rm ELG}(z)D(z)=0.84$ DESI
Collaboration et al. (2016), with $D$ being the linear growth factor
normalized to unity at present time, i.e., $D(z)=G_{\delta}(z)/G_{\delta}(0)$.
Concerned with $b_{2}^{\rm LRG}$ and $b_{2}^{\rm ELG}$, we exploits the halo
fit results. To compute these parameters, Planck $\Lambda$CDM model Planck
Collaboration et al. (2015) is adopted as the fiducial cosmology.

Figure 3: Left: The fractional errors on $D_{A}$ from DESI LRG (upper) and ELG
(lower) samples. The results in units of percentage are plotted as function of
redshift. Solid and dashed curves represent the cases with $\sigma_{8}$
marginalized over and fixed, respectively. Right: Same as in left panel, but
the fractional errors on $H^{-1}$ is shown.

In this section, based on the model described in Sec. II as a theoretical
template, we demonstrate explicitly how well the degeneracy of the parameter
$f\sigma_{8}$ can be broken, adopting specifically the Dark Energy
Spectroscopy Instrument (DESI) DESI Collaboration et al. (2016), which is a
representative galaxy redshift survey of the so-called stage-IV class Albrecht
et al. (2006), dedicated for a precision measurement of BAO and RSD utilizing
the emission-line galaxy (ELG) and luminous red galaxy (LRG) samples (see
Table 1 for parameter specification, which is discussed below).

###  III.1 Fisher matrix formalism

To estimate quantitatively the size of the expected errors, we use the Fisher
matrix formalism. Regarding the model in Sec. II as an observed power
spectrum, the Fisher matrix is evaluated with

$\displaystyle F_{\alpha\beta}=\sum_{k,\mu}\frac{\partial P^{\rm
obs}(k,\mu)}{\partial x_{\alpha}}C_{PP}^{-1}\frac{\partial P^{\rm
obs}(k,\mu)}{\partial x_{\beta}}$ (11)

where the derivative of the power spectrum is taken with respect to the
parameter to estimate, $x_{\alpha}$. The observed power spectrum $P^{\rm obs}$
is related to $P^{\rm(S)}$ given at Sec. II through the Alcock-Paczynski
effect (Ref. Alcock and Paczynski (1979), see also Sec. I):

$\displaystyle P^{\rm obs}(k,\mu)=\Bigl{(}\frac{H}{H_{\rm
fid}}\Bigr{)}\Bigl{(}\frac{D_{A}}{D_{A,{\rm
fid}}}\Bigr{)}^{-2}\,P^{\rm(S)}(q,\nu)\,;\,\,$ (12) $\displaystyle\quad
q=\Bigl{\\{}\Bigl{(}\frac{D_{A}}{D_{A,{\rm
fid}}}\Bigr{)}^{2}(1-\mu^{2})+\Bigl{(}\frac{H}{H_{\rm
fid}}\Bigr{)}^{-2}\mu^{2}\Bigr{\\}}^{1/2}\,k,$
$\displaystyle\quad\nu=\frac{k}{q}\,\Bigl{(}\frac{H}{H_{\rm
fid}}\Bigr{)}^{-1}\mu.$

In Eq. (11), the quantity $C_{PP}$ describes the error covariance of the
measured power spectra, whose dominant contributions are the shot noise
arising from the discreteness of galaxy distribution, and the cosmic variance
due to the limited number of Fourier modes for a finite-volume survey. While
the non-Gaussian contribution, leading to the non-zero off-diagonal
components, is known to play an important role beyond the linear regime, we
will work with the linear Gaussian covariance given below,

$\displaystyle C_{\rm
PP}=\frac{1}{N_{P}}\left[P^{\rm(S)}(k,\mu)+\frac{1}{n_{g}^{X}}\right]^{2}\,,$
(13)

where the quantity $n_{g}^{X}$ denotes the number density of galaxies for a
specific galaxy type $X$, i.e., $X=$LRG or ELG for DESI, whose values are
summarized in Table 1. The quantity $N_{P}$ represents the number of available
Fourier modes, which is estimated to be

$\displaystyle N_{P}=\frac{V_{\rm survey}}{2(2\pi)^{2}}k^{2}\Delta
k\Delta\mu\,$ (14)

with $\Delta k$ and $\Delta\mu$ being the bin width of the power spectrum data
given in the $(k,\mu)$ plane, with $\Delta k=0.01h{\rm Mpc^{-1}}$ and
$\Delta\mu=0.1$. Here, the quantity $V_{\rm survey}$ is the survey volume,
whose values is also listed in Table 1.

Adapting the Gaussian covariance at Eq. (13), the forecast results of our
Fisher matrix analysis would be optimistic. Nevertheless, at the scales where
the perturbative corrections of the next-to-next-to-leading order is still
sub-dominant, the impact of non-Gaussian covariance would be small. The
forecast results presented below can thus give an important guideline toward a
more quantitative analysis. Indeed, to compute the Fisher matrix at Eq. (11)
below, we will conservatively fix the maximum wavenumber to $k_{\rm
max}=0.16\,h$ Mpc-1, at which the non-linearity is still mild at the range of
redshifts, $0.4\leq z\leq 1.6$.

To evaluate Eq. (11), we will use the power spectrum of $P^{\rm(S)}(k,\mu)$ in
Eq. II.1. As the free parameters to estimate, we consider, at each redshift
slice,
$x_{\alpha}=\\{G_{\delta},\,G_{\Theta},\,D_{A},\,H,\,b_{1},b_{2},\,\sigma_{z}\\}$,
where the former two are simply proportional to $\sigma_{8}(z)$ and
$f\,\sigma_{8}(z)$. The parameter $\sigma_{z}$ quantifies the non-perturbative
suppression of the power spectrum by the so-called Fingers-of-God damping that
appears in the function $D_{\rm FoG}$ at Eq. (II.1). We shall below adopt the
Gaussian form of the damping function:

$\displaystyle D_{\rm FoG}(x)=\exp(-x^{2}).$ (15)

where the fiducial $\sigma_{p}=3.3h{\rm Mpc^{-1}}$ computed using linear
theory only. Then, the derivative of the power spectrum $P^{\rm obs}$, given
at Eq. (12), is taken with respect to the parameters, $x_{\alpha}$. Here, the
higher-order bias parameters, $b_{s2}$ and $b_{3\rm{nl}}$, are fixed, adopting
the following relations:

$\displaystyle b_{s2}=-\frac{4}{7}(b_{1}-1)\,,$ $\displaystyle
b_{3\rm{nl}}=\frac{32}{315}(b_{1}-1)\,,$

which are derived assuming the Lagrangian local bias, and are shown to be a
good approximation for halos and galaxies residing at the halo center (e.g.,
Baldauf et al. (2012); Chan et al. (2012); Saito et al. (2014)). The fiducial
values of the linear bias parameters for ELG and LRG are presented in Table 1.

Figure 4: The expected error on $\sigma_{8}$ for the DESI LRG (upper) and ELG
(lower) type samples.

Figure 5: Left: The measured and estimated cosmological constraints on
$\sigma_{8}$ and $f\sigma_{8}$ at 68$\%$ and 95 $\%$ confidence levels are
presented with the blue filled and the black unfilled contours respectively.
The true values are presented with a black dot in the middle of contours.
Those constraints are computed after fully marginalizing other parameters such
as ($D_{A},H^{-1},\sigma_{p},b1,b2$). Right: The same presentation as for the
left panel, but for $D_{A}$ and $H^{-1}$. Those constraints are computed after
fully marginalizing other parameters such as
($\sigma_{8},f\sigma_{8},\sigma_{p},b_{1},b_{2}$).

###  III.2 Results

With the setup described in Sec. III.1, the one-dimensional marginalized
errors on the parameters $f\sigma_{8}$ and $\sigma_{8}$, as well as the
geometric distances, are obtained from the inverse of Fisher matrix, $\Delta
x_{\alpha}=\sqrt{F^{-1}_{\alpha\alpha}}$.

Let us first see how the expected error on the parameter $f\sigma_{8}$ is
altered when the parameter $\sigma_{8}$ or $G_{\delta}$ is taken to be free,
and is marginalized. Fig. 2 plots the results with $\sigma_{8}$ marginalized
(left) and fixed (right), the latter of which is obtained from the inverse of
the sub-matrix $F_{\alpha\beta}$ subtracting the component of $G_{\delta}$. In
each panel, the one-dimensional error on $f\sigma_{8}$ is shown at each
redshift slice for the LRG (top) and ELG (bottom) samples. Remarkably, we
found that the size of the error on $f\sigma_{8}$ almost remains unchanged
irrespective of the treatment of the parameter $\sigma_{8}$, leading to the
constraint at the level of a few percent in both ELG and LRG samples, which is
expected from the stage-IV class surveys.

To elucidate further the impact of marginalizing the $\sigma_{8}$ parameter,
we also plot in Fig. 3 the constraints on the geometric distances. Here, the
fractional errors on the angular diameter distances and Hubble parameters ,
$(D_{A}-D_{A,{\rm fid}})/D_{A,{\rm fid}}$ and $(H-H_{\rm fid})/H_{\rm fid}$,
are respectively shown in left and right panels as function of the redshift.
In each panel, the results with $\sigma_{8}$ marginalized and fixed are
depicted as solid and dashed lines. Again, no notable difference is found
between the two cases, and even the results marginalizing $\sigma_{8}$ reach
at 1 % precision for the constraints on both $D_{A}$ and $H$. In particular,
with the ELG samples, the constraint is improved as increasing the redshift
out to $z\sim 1.5$, with the statistical error down to a sub-percent level.

In Fig. 4, we plot the statistical errors on $\sigma_{8}$ that is separately
determined as a free parameter together with $f\,\sigma_{8}$. As opposed to
the tight constraint on $f\sigma_{8}$, the $\sigma_{8}$ appears poorly
constrained. Considering the fact that the nonlinear corrections in the power
spectrum template that can break the parameter degeneracy are small, this
result is reasonable. Rather, accessing the scales to the weakly nonlinear
regime, the number of available Fourier modes gets increasing, and the
constraints on $f\sigma_{8}$, as well as geometric distances are improved, as
we have seen in Figs. 2 and 3. Nevertheless, using the ELG, the expected error
on $\sigma_{8}$ becomes also improved with redshifts, achieving $10$%
precision at $z=1.5$. Since constraining the growth history of structure at
higher redshifts particularly helps testing and constraining the gravity as
well as the cosmic acceleration, the simultaneous determination of
$f\sigma_{8}$ and $\sigma_{8}$ is beneficial. Combining the power spectrum
with bispectrum, the constraining power will be further improved.

##  IV Breaking the $f\,\sigma_{8}$ degeneracy in $N$-body simulations

Given the survey setup and the theoretical template of power spectrum, the
Fisher matrix analysis in previous section predicts the statistical errors on
each free parameter and their parameter degeneracy, but cannot tell how the
best-fit parameters are accurately determined. In this section, to check the
validity of the forecast results in Sec. III as well as to test the hybrid RSD
model of power spectrum in Sec. II, we here examine the parameter estimation
analysis, based on the Markov chain Monte Carlo (MCMC) technique.

For this purpose, the cosmological $N$-body simulation is carried out by the
publicly available code, GADGET-2 Springel et al. (2001), and we ran in total
$100$ simulations, with $1024^{3}$ particles in comoving periodic cubes of the
side length $1,890\,h^{-1}$ Mpc. Using the output results at $z=0.9$, the halo
catalog was created using the halo finder, ROCKSTAR Behroozi et al. (2013),
with a halo mass range of
$10^{13}\,h^{-1}\,M_{\odot}\,<\,M_{h}\,<\,10^{13.5}\,h^{-1}\,M_{\odot}$, and
the halo number density $n_{h}=1.9\times 10^{-4}\,h^{3}{\rm Mpc}^{-3}$. The
resultant volume of the halo catalog is $6.8\,h^{-3}\,{\rm Gpc}^{3}$, roughly
corresponding to the survey volume of DESI at the redshift slice of $z=0.8--
1.0$. The number density of our halo sample is smaller than those expected
from DESI LRG and ELG samples (see Table 1), but this does not affect the
estimation of statistical errors, as we consider the scales where the shot
noise is subdominant.

The created halo catalog at $z=0.9$ is then used to measure the power spectrum
in redshift space, which is compared with the hybrid RSD model to estimate the
parameters
$(G_{\delta},\,G_{\Theta},\,D_{A},\,H,\,b_{1},\,b_{2},\,\sigma_{z})$,
similarly to what has been done in Sec. III, assuming the Lagrangian local
bias. Adopting the same maximum wavenumber as used in Sec. III, i.e., $k_{\rm
max}=0.16\,h$ Mpc-1, the MCMC results marginalized over other parameters are
shown in Fig. 5, where we plot the two-dimensional error contours on
$(\sigma_{8},\,f\sigma_{8})$ (left) and $(D_{A},\,H^{-1})$ (right). The blue
shaded regions represent the $68$% (dark) and $95$% (light) credible regions
obtained from the MCMC analysis, which exhibit roughly the

In Fig. 5, we also plot the forecast results for the Fisher matrix analysis.
The solid and dashed contours centered at the fiducial parameters, indicated
by the black filled circles, are respectively the the $68$% and $95$% credible
regions for the forecast errors. The MCMC results reproduce well the forecast
results of the Fisher matrix analysis, and the $68$% credible regions
consistently include the fiducial parameters. Thus, with the hybrid RSD model
presented in Sec. II, unbiased parameter estimation is shown to be possible,
with the degeneracy between the parameters $f$ and $\sigma_{8}$ broken.
Although the statistical error of $\sigma_{8}$ is large, this is the first
demonstration that the simultaneous determination of $f$ and $\sigma_{8}$ is
possible only with the power spectrum at weakly nonlinear scales.

##  V Conclusion

Redshift-space distortions (RSD) that appear in the observed galaxy
distributions via spectroscopic surveys offer an important clue to test the
gravity on cosmological scales. Combining the measurement of baryon acoustic
osculations (BAO), RSD can be also used to clarify the nature of cosmic
acceleration. Toward an unambiguous estimation of cosmological parameters, a
crucial issue is not only to improve the precision of RSD measurement, but
also to exploit the method to disentangle the parameter degeneracy inherent in
the observables.

In this paper, on the basis of an accurate template for the galaxy/halo
redshift-space power spectrum, we get access to the weakly nonlinear regime,
and showed that the parameter degeneracy inherent in the linear-theory power
spectrum can be broken. To be precise, in linear theory, the linear growth
rate $f$ and the fluctuation amplitude $\sigma_{8}$ appear in the form of
$f\sigma_{8}$ [Eq. (1)], and this degeneracy cannot be broken unless the
galaxy bias parameter is a priori known or is accurately determined. In order
to break the degeneracy, one way is to go beyond the linear theory. Here, we
use the hybrid RSD model of power spectrum that has been developed in our
previous papers. Based on the perturbation theory calculations, the model
incorporates the higher-order corrections calibrated with $N$-body simulations
into the power spectrum expressions, and the accuracy of predictions is
improved. Taking further the galaxy bias into account, the hybrid model
enables us to get access to the observed galaxy power spectrum at the weakly
nonlinear regime.

Employing the Fisher matrix analysis, we show explicitly that the degeneracy
of the parameter $f\sigma_{8}$ can be broken, and $\sigma_{8}$ is separately
estimated in the presence of galaxy bias. The statistical errors on
$f\sigma_{8}$ as well as the geometric distances $D_{A}$ and $H$, determined
from the BAO via the Alcock-Paczynski effect, are found to remain unchanged,
irrespective of whether we treat $\sigma_{8}$ or the growth factor
$D_{\delta}$ as a free parameter to marginalize or not. As a result, we have
shown that the Dark Energy Survey Instrument, as a representative stage-IV
class galaxy survey, can unambiguously determine $\sigma_{8}$ at the precision
of $\sim 10$% at higher redshifts even if we restrict the accessible scales to
$k\lesssim 0.16\,h$ Mpc-1. Further, performing the Markov chain Monte Carlo
analysis, we explicitly demonstrate that with the hybrid RSD model, the
parameters $f\sigma_{8}$ and $\sigma_{8}$ are simultaneously estimated, and
their fiducial values can be properly recovered, with the statistical errors
fully consistent with the forecast results of the Fisher matrix analysis.

While the analysis in the present paper gives a first explicit demonstration
on how well the parameter degeneracy in the measurement of RSD can be broken,
the ability to achieve this heavily relies on the theoretical template of the
observed power spectrum. Toward a further improvement of the cosmological
constraints in an unbiased way, one needs to develop a model that can get
access to smaller scales. Though the present paper considered the perturbation
theory based model aided by the $N$-body simulations, a simulation based model
such as the so-called emulator would be certainly powerful (e.g. Kobayashi et
al. (2020)). With such a model, the accessible range of wavenumber becomes
broader, and the simultaneous constraints on $f\sigma_{8}$ and $\sigma_{8}$,
as well as the geometric distances, will be improved in a greater precision.
The discussion along the direction of this should be done in the near future.

## Acknowledgments

We would like to thank Takahiro Nishimichi for useful discussions and
comments. Numerical calculations were performed by using a high performance
computing cluster in the Korea Astronomy and Space Science Institute. YZ
acknowledges the support from the Guangdong Basic and Applied Basic Research
Foundation No.2019A1515111098. AT acknowledges the support from MEXT/JSPS
KAKENHI Grant No. JP16H03977, JP17H06359 and JP20H05861. AT was also supported
by JST AIP Acceleration Research Grant NO. JP20317829, Japan.

##  Appendix A Description of hybrid RSD model for matter power spectrum

The evolution of the density and peculiar velocity fields is traditionally
parameterized by $\sigma_{8}$ and $f\sigma_{8}$. Both are related to the
growth functions of $G_{\delta}$ and $G_{\Theta}$ written as below,

$\displaystyle\sigma_{8}^{2}(z)=\frac{G_{\delta}^{2}(z)}{2\pi^{2}}\int
W_{8}^{2}(k)k^{2}P^{i}_{\delta\delta}(k)dk\,,$ (16)
$\displaystyle(f\sigma_{8})^{2}(z)=\frac{G_{\Theta}^{2}(z)}{2\pi^{2}}\int
W_{8}^{2}(k)k^{2}P^{i}_{\delta\delta}(k)dk\,.$ (17)

The window function $W_{8}$ is given by,

$\displaystyle W_{8}=\frac{3j_{1}(kR_{8})}{kR_{8}}\,,$ (18)

where $j_{1}$ is the first-order spherical Bessel function, and
$R_{8}=8h^{-1}{\rm Mpc}$. We intend to simultaneously probe $\sigma_{8}$ and
$f\sigma_{8}$ by exploiting the $G_{\delta}$ and $G_{\Theta}$ parameters.

In the hybrid approach of modelling the RSD effect, the theoretical spectra
$\bar{P}^{\rm th}_{XY}(k,z)$ $(X,Y=\delta\,\,\mbox{or}\,\,\Theta)$ are
computed by the RegPT treatment Taruya et al. (2012), in which all the
statistical quantities including power spectrum are expanded in terms of the
multi-point propagators up to the two–loop order as,

$\displaystyle\bar{P}_{XY}(k,z)=\bar{\Gamma}_{X}^{(1)}(k,z)\bar{\Gamma}_{Y}^{(1)}(k,z)\bar{P}^{i}(k)$
$\displaystyle\quad+2\int\frac{d^{3}\vec{q}}{(2\pi)^{3}}\bar{\Gamma}_{X}^{(2)}(\vec{q},\vec{k}-\vec{q},z)\bar{\Gamma}_{Y}^{(2)}(\vec{q},\vec{k}-\vec{q},z)$
$\displaystyle\quad\quad\quad\times\bar{P}^{i}(q)\bar{P}^{i}(|\vec{k}-\vec{q}|)$
$\displaystyle\quad+6\int\frac{d^{3}\vec{p}d^{3}\vec{q}}{(2\pi)^{6}}$
$\displaystyle\quad\quad\quad\times\bar{\Gamma}_{X}^{(3)}(\vec{p},\vec{q},\vec{k}-\vec{p}-\vec{q},z)\bar{\Gamma}_{Y}^{(3)}(\vec{p},\vec{q},\vec{k}-\vec{p}-\vec{q},z)$
$\displaystyle\quad\quad\quad\times\bar{P}^{i}(p)\bar{P}^{i}(q)\bar{P}^{i}(|\vec{k}-\vec{p}-\vec{q}|).$
(19)

Here $\bar{P}^{i}$ is the initial power spectrum, and $\Gamma_{X}^{(n)}$ is
the $(n+1)$-point propagator. In RegPT treatment, the propagators are
constructed with the standard PT calculations. While the standard PT is
usually applied to a limited range of wavenumbers, incorporating the result of
a partial resummation in the high-$k$ limit, a regularized prediction of the
propagators could be applicable to a larger $k$.

Let us first see the two-point propagator $\bar{\Gamma}_{X}^{(1)}$. The
expression relevant at the two-loop order is summarized as

$\displaystyle\bar{\Gamma}_{X}^{(1)}(k,z)={\rm
exp}\left(-\bar{G}_{\delta}^{2}\,\bar{\gamma}\right)\sum_{n}\bar{G}_{X}\bar{G}_{\delta}^{n-1}\bar{\cal
C}^{(1)}_{n}(\bar{\gamma}).$ (20)

Here, $\bar{\gamma}$ is defined by $\bar{\gamma}=k^{2}\bar{\sigma}^{2}_{\rm
d}/2$ with $\bar{\sigma}_{\rm d}^{2}$ being the dispersion of displacement
field. The $\sigma_{\rm d}$ is computed with the initial power spectrum
through 111Note that choice of the upper bound of the integral has been
specified inTaruya et al. (2012) in somewhat phenomenological way, and there
might be a possible uncertainty. However, it can be absorbed into $\bar{\cal
O}$ in our prescription at least perturbatively. $\bar{\sigma}_{\rm
d}^{2}=\int_{0}^{k/2}(dq/6\pi^{2})\bar{P}^{i}(q)$. The coefficients $\bar{\cal
C}^{(n)}$ in Eq. (20) are expressed in terms of the standard PT results, and
including the theoretical uncertainties, they are given by

$\displaystyle\bar{\cal C}^{(1)}_{1}(\bar{\gamma})$ $\displaystyle=$
$\displaystyle 1,$ $\displaystyle\bar{\cal C}^{(1)}_{3}(\bar{\gamma})$
$\displaystyle=$ $\displaystyle\bar{\gamma}+\bar{\Gamma}_{X,{\rm
1-loop}}^{(1)}(k),$ $\displaystyle\bar{\cal C}^{(1)}_{5}(\bar{\gamma})$
$\displaystyle=$
$\displaystyle\bar{\gamma}^{2}/2+\bar{\gamma}\bar{\Gamma}_{X,{\rm
1-loop}}^{(1)}(k)+\bar{\Gamma}_{X,{\rm 2-loop}}^{(1)}(k)+\bar{\cal
O}^{(1)}_{X,5},$ $\displaystyle\bar{\cal C}^{(1)}_{n}(\bar{\gamma})$
$\displaystyle=$ $\displaystyle\bar{\cal O}^{(1)}_{X,n}\,,$ (21)

and $\bar{\cal C}^{(1)}_{n}=0$ for even numbers of $n$. The $\bar{G}_{X}$
denotes the density ($X=\delta$) and velocity ($X=\Theta$) growth functions
for the fiducial cosmology at the redshift $z$. These growth functions are the
key quantities to be estimated from observations, and are related to the
linear growth factor $D_{+}$ and linear growth rate $f$ through
$G_{\delta}=D_{+}$ and $G_{\theta}=f\,D_{+}$. The function $\Gamma_{X,{\rm
n-loop}}^{(p)}$ represents the standard PT $(p+1)$-point propagator at the
$n$-loop order, whose explicit expression is given in Taruya et al. (2012);
Bernardeau et al. (2014). The quantity $\bar{\cal O}^{(1)}_{X,n}$
characterizes the uncertainties or systematics in PT, which will be later
calibrated with $N$-body simulations. We assume that the uncertainties arise
not only from the higher-order (three-loop) but also from the two-loop order,
partly due to the UV sensitive behavior of the single-stream PT calculation.

Similarly, the expression of the three-point propagator,
$\bar{\Gamma}_{X}^{(2)}(k,z)$ is given by

$\displaystyle\bar{\Gamma}_{X}^{(2)}(k,z)$ $\displaystyle=$ $\displaystyle{\rm
exp}\left(-\bar{G}_{\delta}^{2}\,\bar{\gamma}\right)\sum_{n}\bar{G}_{X}\,\bar{G}_{\delta}^{n-1}\bar{\cal
C}^{(2)}_{n}\,\,;$ $\displaystyle\bar{\cal C}^{(2)}_{2}(\bar{\gamma})$
$\displaystyle=$ $\displaystyle\bar{F}_{X}^{(2)}(\vec{q},\vec{k}-\vec{q}),$
$\displaystyle\bar{\cal C}^{(2)}_{4}(\bar{\gamma})$ $\displaystyle=$
$\displaystyle\frac{\bar{\gamma}}{2}\bar{F}_{X}^{(2)}(\vec{q},\vec{k}-\vec{q})+\bar{\Gamma}_{X,{\rm
1-loop}}^{(2)}(\vec{q},\vec{k}-\vec{q})+\bar{\cal O}^{(2)}_{X,4}\,,$
$\displaystyle\bar{\cal C}^{(2)}_{n}(\bar{\gamma})$ $\displaystyle=$
$\displaystyle\bar{\cal O}^{(2)}_{X,n}\,,$ (22)

and $\bar{\cal C}^{(2)}_{n}=0$ for odd number of $n$. Also, the expression of
the four-point propagator, $\bar{\Gamma}_{X}^{(3)}(k,z)$, relevant at the two-
loop order, is

$\displaystyle\bar{\Gamma}_{X}^{(3)}(k,z)$ $\displaystyle=$ $\displaystyle{\rm
exp}\left(-\bar{G}_{\delta}^{2}\,\gamma\right)\sum_{n}\bar{G}_{X}\,\bar{G}_{\delta}^{n-1}\bar{\cal
C}^{(3)}_{n}\,\,;$ (23) $\displaystyle\bar{\cal C}^{(3)}_{3}(\bar{\gamma})$
$\displaystyle=$
$\displaystyle\bar{F}_{X}^{(3)}(\vec{p},\vec{q},\vec{k}-\vec{p}-\vec{q})+\bar{\cal
O}^{(3)}_{X,3},$ (24) $\displaystyle\bar{\cal C}^{(3)}_{n}(\bar{\gamma})$
$\displaystyle=$ $\displaystyle\bar{\cal O}^{(3)}_{X,n}.$ (25)

Note that $\bar{\cal O}^{(2)}$ and $\bar{\cal O}^{(3)}$ represent the possible
uncertainties.

While the prescription given above is supposed to give an accurate theoretical
prediction at higher redshifts and larger scales, Ref. Song et al. (2018)
reported that the RegPT prediction of $\bar{P}_{XY}$ at a low redshift (to be
precise $z=0.5$) exhibits a small deviation from $N$-body simulations at
$k>0.1\,h\,{\rm Mpc}^{-1}$ . Although this might be partly ascribed to the
systematics in the $N$-body simulations, a lack of higher-order terms as well
as a small systematics in the PT calculations is known to sensitively affect
the high-$k$ prediction. Here, we characterize the difference between the
measured and predicted power spectra by $\bar{P}_{XY}^{\rm res}$. We then
divide the power spectrum into two pieces:

$\displaystyle\bar{P}_{XY}(k,z)=\bar{P}^{\rm th}_{XY}(k,z)+\bar{P}^{\rm
res}_{XY}(k,z),$ (26)

where $\bar{P}^{\rm th}_{XY}$ represents the PT prediction with $\bar{\cal
O}^{(m)}_{X,n}\rightarrow 0$. Collecting all the uncertainties introduced in
the multi-point propagators, the residual power spectrum $\bar{P}_{XY}^{\rm
res}$ is schematically expressed as

$\displaystyle\bar{P}^{\rm
res}_{XY}=\bar{G}_{X}\bar{G}_{Y}\bar{G}_{\delta}^{4}\,\Biggl{\\{}\left[{\cal
O}^{(1)}_{Y,5}+{\rm higher}\right]\bar{P}^{i}$ $\displaystyle+\left[\bar{\cal
O}^{(1)}_{X,5}+{\rm higher}\right]\bar{P}^{i}+\int\left[\bar{\cal
O}^{(2)}_{Y,4}\bar{F}_{Y}^{(2)}+{\rm higher}\right]\bar{P}^{i}\bar{P}^{i}$
$\displaystyle+\int\left[\bar{\cal O}^{(2)}_{X,4}\bar{F}_{X}^{(2)}+{\rm
higher}\right]\bar{P}^{i}\bar{P}^{i}$ $\displaystyle+\int\int\left[\bar{\cal
O}^{(3)}_{Y,3}\bar{F}_{Y}^{(3)}+{\rm
higher}\right]\bar{P}^{i}\bar{P}^{i}\bar{P}^{i}$
$\displaystyle+\bar{\int}\int\left[\bar{\cal
O}^{(3)}_{X,3}\bar{F}_{X}^{(3)}+{\rm
higher}\right]\bar{P}^{i}\bar{P}^{i}\bar{P}^{i}\Biggr{\\}}.$ (27)

Here, the uncertainty $\bar{\cal O}^{(m)}_{X,n}$ is assumed to be small, and
to be perturbatively treated. The expression implies that apart from a
detailed scale-dependent behavior, time dependence is characterized by
$G_{X}G_{Y}G_{\delta}^{4}$. Thus, once we calibrate the $\bar{P}_{XY}^{\rm
res}$ at a given redshift, we may use it for the prediction at another
redshift by simply rescaling the calibrated residuals. Furthermore, for the
cosmological models close to the fiducial model, the scale dependence of the
higher-order PT corrections is generally insensitive to the cosmology, and we
may also apply the calibrated $\bar{P}_{XY}^{\rm res}$ to other cosmological
models.

Next we introduce the way to calculate the dark matter higher order terms.
They are incorporated with Eq. (7-10) to calculate the halo higher order
terms. To begin with, let us consider the $A$ term. From the explicit form,
the $A$ term is divided into six pieces. Here, we specifically write down the
expressions in fiducial cosmological model:

$\displaystyle\bar{A}(k,\mu)$ $\displaystyle=$ $\displaystyle j_{1}\,\int
d^{3}\mbox{\boldmath$x$}\,\,e^{i\mbox{\boldmath$k$}\cdot\mbox{\boldmath$x$}}\,\,\langle
A_{1}A_{2}A_{3}\rangle_{c}$ (28) $\displaystyle=$
$\displaystyle\sum_{n=1}^{6}\bar{\cal A}_{n}$

Note again that the barred quantities indicate those computed/measured in
fiducial cosmological model. The explicit form of $\bar{\cal A}_{n}$ is given
below:

$\displaystyle\bar{\cal A}_{1}$ $\displaystyle=$ $\displaystyle 2j_{1}\,\int
d^{3}\mbox{\boldmath$x$}\,\,e^{i\mbox{\boldmath$k$}\cdot\mbox{\boldmath$x$}}\,\,\langle
u_{z}(\mbox{\boldmath$r$})\delta(\mbox{\boldmath$r$})\delta(\mbox{\boldmath$r$}^{\prime})\rangle_{c},$
(29) $\displaystyle\bar{\cal A}_{2}$ $\displaystyle=$ $\displaystyle
j_{1}\,\int
d^{3}\mbox{\boldmath$x$}\,\,e^{i\mbox{\boldmath$k$}\cdot\mbox{\boldmath$x$}}\,\,\langle
u_{z}(\mbox{\boldmath$r$})\delta(\mbox{\boldmath$r$})\,\nabla_{z}u_{z}(\mbox{\boldmath$r$}^{\prime})\rangle_{c},$
(30) $\displaystyle\bar{\cal A}_{3}$ $\displaystyle=$ $\displaystyle
j_{1}\,\int
d^{3}\mbox{\boldmath$x$}\,\,e^{i\mbox{\boldmath$k$}\cdot\mbox{\boldmath$x$}}\,\,\langle
u_{z}(\mbox{\boldmath$r$})\,\nabla_{z}u_{z}(\mbox{\boldmath$r$})\delta(\mbox{\boldmath$r$}^{\prime})\rangle_{c},$
(31) $\displaystyle\bar{\cal A}_{4}$ $\displaystyle=$ $\displaystyle
2j_{1}\,\int
d^{3}\mbox{\boldmath$x$}\,\,e^{i\mbox{\boldmath$k$}\cdot\mbox{\boldmath$x$}}\,\,\langle
u_{z}(\mbox{\boldmath$r$})\,\nabla_{z}u_{z}(\mbox{\boldmath$r$})\,\nabla_{z}u_{z}(\mbox{\boldmath$r$}^{\prime})\rangle_{c},$
(32) $\displaystyle\bar{\cal A}_{5}$ $\displaystyle=$ $\displaystyle
j_{1}\,\int
d^{3}\mbox{\boldmath$x$}\,\,e^{i\mbox{\boldmath$k$}\cdot\mbox{\boldmath$x$}}\,\,\langle-\delta(\mbox{\boldmath$r$})u_{z}(\mbox{\boldmath$r$}^{\prime})\,\nabla_{z}u_{z}(\mbox{\boldmath$r$}^{\prime})\rangle_{c},$
(33) $\displaystyle\bar{\cal A}_{6}$ $\displaystyle=$ $\displaystyle
j_{1}\,\int
d^{3}\mbox{\boldmath$x$}\,\,e^{i\mbox{\boldmath$k$}\cdot\mbox{\boldmath$x$}}\,\,\langle-\,\nabla_{z}u_{z}(\mbox{\boldmath$r$})u_{z}(\mbox{\boldmath$r$}^{\prime})\delta(\mbox{\boldmath$r$}^{\prime})\rangle_{c}.$
(34)

These terms are measured from $N$-body simulations according to Ref. Zheng and
Song (2016). To apply the measured results to the prediction in other
cosmological models, we assume the scaling ansatz, as similarly adopted in the
prediction of power spectrum $P_{\rm XY}$. That is, assuming that the scale-
dependence of each term is insensitive to the cosmology, the prediction of
each term is made by simply rescaling the measured results. The proposition
made here is that the time-dependence of each term is approximately determined
by the leading-order growth factor dependence of $u_{z}$ and $\delta$. Then,
$A$ term in general cosmological model is expressed as

$\displaystyle A(k,\mu)=\sum_{n=1}^{6}{\cal A}_{n}$ (35)
$\displaystyle=\left(G_{\delta}/\bar{G}_{\delta}\right)^{2}\left(G_{\Theta}/\bar{G}_{\Theta}\right)\bar{\cal
A}_{1}+\left(G_{\delta}/\bar{G}_{\delta}\right)\left(G_{\Theta}/\bar{G}_{\Theta}\right)^{2}\bar{\cal
A}_{2}$
$\displaystyle+\left(G_{\delta}/\bar{G}_{\delta}\right)\left(G_{\Theta}/\bar{G}_{\Theta}\right)^{2}\bar{\cal
A}_{3}+\left(G_{\Theta}/\bar{G}_{\Theta}\right)^{3}\bar{\cal A}_{4}$
$\displaystyle+\left(G_{\delta}/\bar{G}_{\delta}\right)\left(G_{\Theta}/\bar{G}_{\Theta}\right)^{2}\bar{\cal
A}_{5}+\left(G_{\delta}/\bar{G}_{\delta}\right)\left(G_{\Theta}/\bar{G}_{\Theta}\right)^{2}\bar{\cal
A}_{6}$

We then apply the same strategy to other higher-order corrections, $B$, $F$
and $T$. Dividing these corrections into several pieces, the scaling ansatz
leads to the following predictions:

$\displaystyle B(k,\mu)=\sum_{n=1}^{4}{\cal B}_{n}$ (36)
$\displaystyle=\left(G_{\delta}/\bar{G}_{\delta}\right)^{2}\left(G_{\Theta}/\bar{G}_{\Theta}\right)^{2}\bar{\cal
B}_{1}+\left(G_{\delta}/\bar{G}_{\delta}\right)\left(G_{\Theta}/\bar{G}_{\Theta}\right)^{3}\bar{\cal
B}_{2}$
$\displaystyle+\left(G_{\delta}/\bar{G}_{\delta}\right)\left(G_{\Theta}/\bar{G}_{\Theta}\right)^{3}\bar{\cal
B}_{3}+\left(G_{\Theta}/\bar{G}_{\Theta}\right)^{4}\bar{\cal B}_{4}$
$\displaystyle F(k,\mu)=\sum_{n=1}^{3}{\cal F}_{n}$ (37)
$\displaystyle=\left(G_{\delta}/\bar{G}_{\delta}\right)^{2}\left(G_{\Theta}/\bar{G}_{\Theta}\right)^{2}\bar{\cal
F}_{1}+\left(G_{\delta}/\bar{G}_{\delta}\right)\left(G_{\Theta}/\bar{G}_{\Theta}\right)^{3}\bar{\cal
F}_{2}$ $\displaystyle+\left(G_{\Theta}/\bar{G}_{\Theta}\right)^{4}\bar{\cal
F}_{3}$ $\displaystyle T(k,\mu)=\sum_{n=1}^{7}{\cal T}_{n}$ (38)
$\displaystyle=\left(G_{\delta}/\bar{G}_{\delta}\right)^{2}\left(G_{\Theta}/\bar{G}_{\Theta}\right)^{2}\bar{\cal
T}_{1}+\left(G_{\delta}/\bar{G}_{\delta}\right)\left(G_{\Theta}/\bar{G}_{\Theta}\right)^{3}\bar{\cal
T}_{2}$
$\displaystyle+\left(G_{\delta}/\bar{G}_{\delta}\right)\left(G_{\Theta}/\bar{G}_{\Theta}\right)^{3}\bar{\cal
T}_{3}+\left(G_{\Theta}/\bar{G}_{\Theta}\right)^{4}\bar{\cal T}_{4}$
$\displaystyle+\left(G_{\delta}/\bar{G}_{\delta}\right)^{2}\left(G_{\Theta}/\bar{G}_{\Theta}\right)^{2}\bar{\cal
T}_{5}+\left(G_{\delta}/\bar{G}_{\delta}\right)\left(G_{\Theta}/\bar{G}_{\Theta}\right)^{3}\bar{\cal
T}_{6}$ $\displaystyle+\left(G_{\Theta}/\bar{G}_{\Theta}\right)^{4}\bar{\cal
T}_{7}$

Here, the quantities $\bar{\cal B}_{n}$, $\bar{\cal F}_{n}$, and $\bar{\cal
T}_{n}$ are measured in the fiducial cosmological model. Both $\sigma_{8}$ and
$f\sigma_{8}$ are estimated from the measured growth functions of $G_{\delta}$
and $G_{\Theta}$.

## References

  * Riess et al. (1998) A. G. Riess, A. V. Filippenko, P. Challis, A. Clocchiatti, A. Diercks, P. M. Garnavich, R. L. Gilliland, C. J. Hogan, S. Jha, R. P. Kirshner, et al., AJ 116, 1009 (1998), eprint astro-ph/9805201.
  * Perlmutter et al. (1999) S. Perlmutter, G. Aldering, G. Goldhaber, R. A. Knop, P. Nugent, P. G. Castro, S. Deustua, S. Fabbro, A. Goobar, D. E. Groom, et al., Astrophys. J. 517, 565 (1999), eprint astro-ph/9812133.
  * Riess et al. (2016) A. G. Riess, L. M. Macri, S. L. Hoffmann, D. Scolnic, S. Casertano, A. V. Filippenko, B. E. Tucker, M. J. Reid, D. O. Jones, J. M. Silverman, et al., Astrophys. J. 826, 56 (2016), eprint 1604.01424.
  * Wong et al. (2020) K. C. Wong, S. H. Suyu, G. C. F. Chen, C. E. Rusu, M. Millon, D. Sluse, V. Bonvin, C. D. Fassnacht, S. Taubenberger, M. W. Auger, et al., MNRAS 498, 1420 (2020), eprint 1907.04869.
  * Jedamzik et al. (2020) K. Jedamzik, L. Pogosian, and G.-B. Zhao, arXiv e-prints arXiv:2010.04158 (2020), eprint 2010.04158.
  * Lue (2006) A. Lue, Physics reports 423, 1 (2006), eprint astro-ph/0510068.
  * Frieman et al. (2008) J. A. Frieman, M. S. Turner, and D. Huterer, Ann. Rev. Astron. Astrophys. 46, 385 (2008), eprint 0803.0982.
  * Li et al. (2011) M. Li, X.-D. Li, S. Wang, and Y. Wang, Communications in Theoretical Physics 56, 525 (2011), eprint 1103.5870.
  * Clifton et al. (2012) T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Physics reports 513, 1 (2012), eprint 1106.2476.
  * Joyce et al. (2016) A. Joyce, L. Lombriser, and F. Schmidt, Annual Review of Nuclear and Particle Science 66, 95 (2016), eprint 1601.06133.
  * Koyama (2016) K. Koyama, Reports on Progress in Physics 79, 046902 (2016), eprint 1504.04623.
  * Nojiri et al. (2017) S. Nojiri, S. D. Odintsov, and V. K. Oikonomou, Physics reports 692, 1 (2017), eprint 1705.11098.
  * Arun et al. (2017) K. Arun, S. B. Gudennavar, and C. Sivaram, Advances in Space Research 60, 166 (2017), eprint 1704.06155.
  * Wang et al. (2017) S. Wang, Y. Wang, and M. Li, Physics reports 696, 1 (2017), eprint 1612.00345.
  * Brax (2018) P. Brax, Reports on Progress in Physics 81, 016902 (2018).
  * Ishak (2019) M. Ishak, Living Reviews in Relativity 22, 1 (2019), eprint 1806.10122.
  * Alcock and Paczynski (1979) C. Alcock and B. Paczynski, Nature (London) 281, 358 (1979).
  * Song et al. (2015) Y.-S. Song, A. Taruya, and A. Oka, JCAP 2015, 007 (2015), eprint 1502.03099.
  * Gil-Marín et al. (2015) H. Gil-Marín, J. Noreña, L. Verde, W. J. Percival, C. Wagner, M. Manera, and D. P. Schneider, MNRAS 451, 539 (2015), eprint 1407.5668.
  * Taruya et al. (2010) A. Taruya, T. Nishimichi, and S. Saito, Phys. Rev. D 82, 063522 (2010), eprint 1006.0699.
  * Taruya et al. (2013) A. Taruya, T. Nishimichi, and F. Bernardeau, Phys. Rev. D 87, 083509 (2013), eprint 1301.3624.
  * Zheng and Song (2016) Y. Zheng and Y.-S. Song, JCAP 8, 050 (2016), eprint 1603.00101.
  * Song et al. (2018) Y.-S. Song, Y. Zheng, A. Taruya, and M. Oh, ArXiv e-prints (2018), eprint 1801.04950.
  * Zheng et al. (2019) Y. Zheng, Y.-S. Song, and M. Oh, JCAP 2019, 013 (2019), eprint 1807.08115.
  * Ivanov et al. (2020) M. M. Ivanov, M. Simonović, and M. Zaldarriaga, JCAP 05, 042 (2020), eprint 1909.05277.
  * D’Amico et al. (2020) G. D’Amico, J. Gleyzes, N. Kokron, K. Markovic, L. Senatore, P. Zhang, F. Beutler, and H. Gil-Marín, JCAP 05, 005 (2020), eprint 1909.05271.
  * Nishimichi et al. (2020) T. Nishimichi, G. D’Amico, M. M. Ivanov, L. Senatore, M. Simonović, M. Takada, M. Zaldarriaga, and P. Zhang (2020), eprint 2003.08277.
  * Zheng et al. (2015) Y. Zheng, P. Zhang, and Y. Jing, Phys. Rev. D 91, 043523 (2015), eprint 1409.6809.
  * Taruya et al. (2012) A. Taruya, F. Bernardeau, T. Nishimichi, and S. Codis, Phys. Rev. D 86, 103528 (2012), eprint 1208.1191.
  * Bernardeau et al. (2008) F. Bernardeau, M. Crocce, and R. Scoccimarro, Phys. Rev. D 78, 103521 (2008), eprint 0806.2334.
  * McDonald and Roy (2009) P. McDonald and A. Roy, JCAP 8, 020 (2009), eprint 0902.0991.
  * Saito et al. (2014) S. Saito, T. Baldauf, Z. Vlah, U. Seljak, T. Okumura, and P. McDonald, Phys. Rev. D 90, 123522 (2014), eprint 1405.1447.
  * de Mattia et al. (2020) A. de Mattia, V. Ruhlmann-Kleider, A. Raichoor, A. J. Ross, A. Tamone, C. Zhao, S. Alam, S. Avila, E. Burtin, J. Bautista, et al., MNRAS (2020), eprint 2007.09008.
  * DESI Collaboration et al. (2016) DESI Collaboration, A. Aghamousa, J. Aguilar, S. Ahlen, S. Alam, L. E. Allen, C. Allende Prieto, J. Annis, S. Bailey, C. Balland, et al., ArXiv e-prints (2016), eprint 1611.00036.
  * Planck Collaboration et al. (2015) 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, et al., ArXiv e-prints (2015), eprint 1502.01589.
  * Albrecht et al. (2006) A. Albrecht, G. Bernstein, R. Cahn, W. L. Freedman, J. Hewitt, W. Hu, J. Huth, M. Kamionkowski, E. W. Kolb, L. Knox, et al., arXiv e-prints astro-ph/0609591 (2006), eprint astro-ph/0609591.
  * Baldauf et al. (2012) T. Baldauf, U. Seljak, V. Desjacques, and P. McDonald, Phys. Rev. D 86, 083540 (2012), eprint 1201.4827.
  * Chan et al. (2012) K. C. Chan, R. Scoccimarro, and R. K. Sheth, Phys. Rev. D 85, 083509 (2012), eprint 1201.3614.
  * Springel et al. (2001) V. Springel, N. Yoshida, and S. D. M. White, New Astronomy 6, 79 (2001), eprint arXiv:astro-ph/0003162.
  * Behroozi et al. (2013) P. S. Behroozi, R. H. Wechsler, and H.-Y. Wu, Astrophys. J. 762, 109 (2013), eprint 1110.4372.
  * Kobayashi et al. (2020) Y. Kobayashi, T. Nishimichi, M. Takada, R. Takahashi, and K. Osato, Phys. Rev. D 102, 063504 (2020), eprint 2005.06122.
  * Bernardeau et al. (2014) F. Bernardeau, A. Taruya, and T. Nishimichi, Phys. Rev. D 89, 023502 (2014), eprint 1211.1571.


# Direct determination of the excitation energy of quasi-stable isomer 180mTa

D.A. Nesterenko  K. Blaum  P. Delahaye  S. Eliseev  T. Eronen  P. Filianin  Z.
Ge  M. Hukkanen  A. Kankainen  Yu.N. Novikov  A.V. Popov  A. Raggio  M.
Stryjczyk  V. Virtanen  University of Jyväskylä, Accelerator Laboratory,
Department of Physics, University of Jyväskylä, Finland  Max-Planck Institut
f. Kernphysik, Heidelberg, Germany  GANIL, CEA/DSM-CNRS/IN2P3, Bd Henri
Becquerel, 14000 Caen, France  GSI Helmholtzzentrum f. Schwerionenforschung,
Darmstadt, Germany  NRC “Kurchatov Institute” - Petersburg Nuclear Physics
Institute, Gatchina, Russia

###### Abstract

180mTa is a naturally abundant quasi-stable nuclide and the longest-lived
nuclear isomer known to date. It is of interest for, among others, the search
for dark matter, for the development of a gamma laser and for astrophysics. So
far, its excitation energy has not been measured directly but has been based
on an evaluation of available nuclear reaction data. We have determined the
excitation energy of this isomer with high accuracy using the Penning-trap
mass spectrometer JYFLTRAP. The determined mass difference between the ground
and isomeric states of 180Ta yields an excitation energy of 76.79(55) keV for
180mTa. This is the first direct measurement of the excitation energy and
provides a better accuracy than the previous evaluation value, 75.3(14) keV.

###### keywords:

Penning trap , high-precision mass spectrometry , nuclear isomers

††journal: PRC

##  1 Introduction

Nuclear isomers were discovered around 100 years ago [1]. Today, around 1938
excited isomeric ($T_{1/2}\geq 100$ ns) states are known [2]. 180mTa isomer
with spin-parity $J^{\pi}=9^{-}$ and excitation energy of around 75 keV is a
unique nuclear isomer. The half-life of 180mTa is longest of all isomeric
states, longer than $4.5\times 10^{16}$ y [2]. The isomer is much longer-
living than the $J^{\pi}=1^{+}$ ground state of 180Ta, which has a half-life
of 8.1 h. As a consequence, 180Ta is present in nature in its isomeric state,
making a 0.01201(8)% fraction of natural tantalum.

Due to the large spin difference between the $1^{+}$ ground state and the
$9^{-}$ isomer, the internal transition to the ground state is not possible.
The ground state decays via $\beta^{-}$ decay to 180W and via electron capture
($EC$) decay to 180Hf (see Fig. 1). For the 180mTa isomer, the beta decay
branches are non-existent or extremely small as there are no $8^{-}$, $9^{-}$,
$10^{-}$ states within the $Q_{\beta}-/Q_{EC}$ window.

180mTa has been studied in detail from many perspectives. The synthesis of
180mTa, the rarest stable isotope in nature, is still not well understood.
Thermal excitation and de-excitations are also affected by nuclear structure
of 180mTa, which is a prolate nucleus with the projection of the spin on the
nuclear deformation axis $K=9$. The bands built on the $K^{\pi}=9^{-}$ isomer
and $K^{\pi}=1^{+}$ ground state have been studied in detail [3, 4, 5, 6]. At
high stellar temperatures, the isomer can be depopulated via thermal
excitation to higher-lying states, which eventually de-excite to the short-
lived ground state. The production of 180mTa via the slow neutron-capture
process ($s$-process) as well as $p$\- and $\nu$-processes has been studied
[7, 8]. It has been proposed as a good candidate for the development of a
gamma laser [3]. 180mTa has also potential to be used in experiments devoted
to the search for dark matter [9].

So far, a direct decay of 180mTa has never been observed. Its excitation
energy has been evaluated based on various nuclear reactions, resulting in
values 75.3(14) keV [2] and 77.2(12) keV [10]. There are several examples in
the literature, where an evaluated value has turned out to be in a strong
disagreement with the result achieved by direct Penning-trap mass measurements
(for example, a $10\sigma$ discrepancy was found for 102Pd–102Ru in [11]).
Thus, a direct mass measurement is of outermost importance. In this work, we
have used Penning-trap mass spectrometry to reliably measure the key
spectroscopic parameter of the isomer – its excitation energy – as a mass
difference between the isomeric and ground state of 180Ta.

Figure 1: (colour online) Decay scheme of 180Ta. Energy levels in 180Ta and
$Q$ values [2] are given in keV. Due to the high spin difference, a
depopulation of the naturally abundant 76.79 keV state in 180mTa can only
occur by a photoexcitation into resonance intermediate states which have a
decay branch into the ground state. The ground state, in turn, decays either
by $\beta^{-}$ decay to 180W or by electron capture ($EC$) to 180Hf.

##  2 Experimental method and results

The measurement was performed with singly charged 180Ta+ and 180mTa+ ions with
the Penning-trap mass spectrometer JYFLTRAP [12] at the Ion Guide Isotope
Separator On-Line (IGISOL) facility [13]. Both 180Ta+ and 180mTa+ ions were
simultaneously produced using nuclear reactions with a 40 MeV proton beam
impinging into a ${}^{\text{nat}}$Ta target with a thickness of about 5
mg/cm2. The experimental reaction cross section for a production of 180Ta at
this proton energy is a few hundred mb [14], and a similar cross section is
expected for the production of 180mTa according to TALYS simulations [15]. The
reaction products were thermalized in the IGISOL gas cell filled with helium
gas at a pressure of about 150 mbar. A large fraction of ions ends up as
singly-charged. The ions are extracted from the IGISOL gas cell and guided via
a sextupole ion guide [16] to a high-vacuum region. Then they were accelerated
to 30 keV energy with subsequent mass separation by a $55^{\circ}$ dipole
magnet. The continuous ion beam with the selected mass number of $A=180$ was
injected into a gas-filled radio-frequency quadrupole [17], where it is
converted into compact ion bunches. Finally, the ion bunches are guided to the
JYFLTRAP double Penning trap mass spectrometer.

Figure 2: (colour online) Projection of the cyclotron phase image (blue bins)
and the trap center (green bins) onto the position-sensitive MCP detector for
the ${}^{180}\text{Ta}^{+}$ and ${}^{180m}\text{Ta}^{+}$ ions in a single
cyclotron frequency measurement with the PI-ICR method. The phase accumulation
time $t_{acc}$ was about of 1125 ms.

In the mass spectrometer the ions were cooled, centered and purified in the
first preparation Penning trap by means of the mass-selective buffer gas
cooling technique [18]. The process requires approximately 300 ms. Since the
mass difference between the isobaric contaminant ${}^{180}\text{W}^{+}$ and
the ion of interest ${}^{180}\text{Ta}^{+}$ is only about 700 keV/c2,
corresponding to about 2.5 Hz in cyclotron frequency difference, the buffer
gas cooling technique is unable to resolve them. To get rid of
${}^{180}\text{W}^{+}$ ions, after cooling in the preparation trap the ions
were transported to the second (measurement) trap where a 200-ms dipolar pulse
at the corresponding reduced cyclotron frequency $\nu_{+}$ was applied to
excite these ions to large cyclotron radii, larger then the radius of the
diaphragm between the traps. Then the ions were sent back to the first trap,
and only ions of interest that remain cooled at the trap center passed the
diaphragm and were subject to the final cooling that takes 160 ms. Finally,
only ${}^{180}\text{Ta}^{+}$ ions, both the ground and the isomeric state,
were transported to the measurement trap for the final determination of their
mass. Tantalum is a highly reactive chemical element and easily forms
${}^{180}\mathrm{T}^{16}\mathrm{O}^{+}$ compound in the measurement trap by
interactions of the tantalum ions with residual gas particles. To get rid of
this constantly forming tantalum oxide contamination, a continuous dipolar
excitation at the corresponding $\nu_{+}$ frequency was constantly applied
over the entire duration of the measurement. Since the masses and thus the
frequencies of ${}^{180}\mathrm{T}^{16}\mathrm{O}^{+}$ and
${}^{180}\mathrm{Ta}^{+}$ ions differ significantly, the continuous cleaning
excitation pulse does not affect the actual mass measurement of
${}^{180}\text{Ta}^{+}$ states.

In Penning trap mass spectrometry, the mass $m$ of an ion is determined by
measuring the ion’s cyclotron frequency $\nu_{c}$

$\nu_{c}=\frac{1}{2\pi}\frac{q}{m}B,$ (1)

where $q/m$ is the ion’s charge-to-mass ratio and $B$ is the magnetic field
strength. The cyclotron frequency is measured with the phase-imaging ion-
cyclotron-resonance (PI-ICR) technique that is described in detail in [19, 20,
21]. Here only the main concept and the details specific to the presented
measurement are given. The PI-ICR method is based on the observation of a
phase evolution of the magnetron and cyclotron ion motions with frequencies
$\nu_{-}$ and $\nu_{+}$, respectively, by projecting the ion’s position in the
trap onto a position-sensitive detector. The projection of the position of the
ion that performs pure magnetron and cyclotron motions are called $\nu_{-}$
and $\nu_{+}$ phase images (spots), respectively. The positions of the phase
images are described with polar angles $\alpha_{-}$ and $\alpha_{+}$,
respectively, with respect to the trap center. The cyclotron frequency
$\nu_{c}$ calculated as a sum of two radial frequencies [22] is determined
from the angle $\alpha_{c}=\alpha_{+}-\alpha_{-}$ between these two phase
images as

$\nu_{c}=\nu_{-}+\nu_{+}=\frac{\alpha_{c}+2\pi n}{2\pi t_{\mathrm{acc}}},$ (2)

where $n$ is the full number of revolutions during the phase accumulation time
$t_{\text{acc}}$. For $\nu_{c}$ frequency measurement of a certain ion species
the phase spots and the center spot were alternately accumulated with a period
of about 3 minutes. In turn, measurements between the two ion species were
alternated after 4 such rounds, i.e with a period of about 12 minutes.

Figure 3: (colour online) Cyclotron frequency ratios
$R=\nu_{c}(^{180}$Ta${}^{+})/\nu_{c}(^{180m}$Ta${}^{+})$ determined for each
2-hour measurement period. The grey band represents the total 1$\sigma$
uncertainty of the weighted mean frequency ratio $\bar{R}=1.0000004581(33)$.

The phase accumulation time $t_{\text{acc}}$ of 950 ms has been used to
resolve the ground and isomeric state in 180Ta. As a crosscheck, a part of the
data was acquired with $t_{\text{acc}}=1125$ ms (Fig. 2) to ensure that the
cyclotron spots on the detector are not contaminated by any possible isobaric
impurities. Since both states are simultaneously produced, have similar
yields, and cannot be efficiently separated prior to the actual PI-ICR
measurement, for minimization of the cyclotron frequency shift due to an ion-
ion interaction between these different ion species, the measurement was
carried out with a low count rate of approximately one detected ion per five
bunches. The measurement settings for both the ground and isomeric state are
the same except for the delay and frequency of the conversion pulse, which are
adjusted to minimize angle $\alpha_{c}$ between the magnetron and cyclotron
phases for each state. The average angle $\alpha_{c}$ in the measurements did
not exceed a few degrees to minimize the systematic shifts due to the
distortion of the phase projections.

Due to the intentionally low count rate, a couple of hours were required to
collect reasonable statistics for an individual $\nu_{c}$ determination. Thus,
the total measurement time was divided into 23 approximately 2-hour periods.
For each period the $\nu_{c}$ frequency for both ion species was determined
and a single cyclotron frequency ratio
$R=\nu_{c}(^{180}\mathrm{Ta}^{+})/\nu_{c}(^{180\text{m}}\mathrm{Ta}^{+})$ was
calculated, as shown in Fig. 3. The final cyclotron-frequency ratio
$\bar{R}=1.0000004581(33)$ was obtained as the weighted mean of these single
ratios with the uncertainty taken as the larger of the internal and external
statistical uncertainties [23]. The associated Birge ratio was about a unity.
The systematic uncertainties specific to the PI-ICR method are discussed in
[24]. The mass-dependent systematic effects for mass doublets are negligible
[25], thus the statistical uncertainty is dominant in the final uncertainty.

The excitation energy $E^{*}$ of 180mTa was calculated from the cyclotron
frequency ratio as

$E^{*}=\left(M(^{180m}\text{Ta})-M(^{180}\text{Ta})\right)c^{2}=\\\
(\bar{R}-1)\left(M(^{180}\text{Ta})-m_{e}\right)c^{2},$ (3)

where $M$(180Ta) and $M$(180mTa) are the atomic masses of 180Ta and 180mTa,
respectively, and $m_{e}$ is the electron mass. The uncertainty of the
tantalum mass value $\delta M(^{180}\text{Ta})=2.1$ keV/$c^{2}$ [26] does not
affect the precision since the first term $(\bar{R}-1)<10^{-5}$. For the same
reason, the binding energy of the valence electron can be neglected.
Ultimately, the determined excitation energy of ${}^{180\text{m}}$Ta is
$E^{*}=76.79(55)$ keV.

##  3 Conclusion

${}^{180\text{m}}$Ta is an unusual isomeric state. It is the only excited
nuclear state present in nature and the longest-lived isomer. It is also a
state of potential interest in various fields of physics. However, its main
characteristics, namely the excitation energy, had not been measured directly.
In this work, we have measured the cyclotron frequency ratio between the
${}^{180}\text{Ta}^{+}$ ground-state and isomeric-state ions in a single
experiment at the Penning-trap facility JYFLTRAP. With the obtained cyclotron
frequency ratio, $R=1.0000004581(33)$, we calculated the excitation energy of
the longest-lived isomer ${}^{180m}\text{Ta}$ to be $E^{*}=76.79(55)$ keV. Our
new value is a factor of two more precise and in a good agreement with the
previously known literature values [2, 10] evaluated on the basis of various
nuclear reactions with ${}^{180}\text{Ta}$. The new accurate value improves
the precision of the excitation energy and related reaction $Q$ values
relevant for stellar nucleosynthesis [8, 3], driven gamma emission in 180Ta
and gamma laser development [27], and to a dark matter search [28].

## Acknowledgements

The experiment was performed in 2021. This work has been supported by the
European Union’s Horizon 2020 research and innovation programme under grant
agreements No 771036 (ERC CoG MAIDEN) and No 861198–LISA–H2020-MSCA-ITN-2019
and by the Academy of Finland under projects No 295207 and 327629.

## References

  * Walker and Podolyák [2020] P. Walker and Z. Podolyák, Physica Scripta 95, 044004 (2020).
  * Kondev  _et al._ [2021] F. Kondev, M. Wang, W. Huang, S. Naimi, and G. Audi, Chinese Physics C 45, 030001 (2021).
  * Walker  _et al._ [2001] P. M. Walker, G. D. Dracoulis, and J. J. Carroll, Phys. Rev. C 64, 061302 (2001).
  * Dracoulis  _et al._ [1998] G. D. Dracoulis, S. M. Mullins, A. P. Byrne, F. G. Kondev, T. Kibédi, S. Bayer, G. J. Lane, T. R. McGoram, and P. M. Davidson, Phys. Rev. C 58, 1444 (1998).
  * Saitoh  _et al._ [1999] T. Saitoh, N. Hashimoto, G. Sletten, R. Bark, S. Törmänen, M. Bergström, K. Furuno, K. Furutaka, G. Hagemann, T. Hayakawa, T. Komatsubara, A. Maj, S. Mitarai, M. Oshima, J. Sampson, T. Shizuma, and P. Varmette, Nuclear Physics A 660, 121 (1999).
  * Dracoulis  _et al._ [2000] G. D. Dracoulis, T. Kibédi, A. P. Byrne, R. A. Bark, and A. M. Baxter, Phys. Rev. C 62, 037301 (2000).
  * Käppeler  _et al._ [2011] F. Käppeler, R. Gallino, S. Bisterzo, and W. Aoki, Rev. Mod. Phys. 83, 157 (2011).
  * Mohr [2009] P. Mohr, PoS NIC X, 081 (2009).
  * Pospelov  _et al._ [2020] M. Pospelov, S. Rajendran, and H. Ramani, Phys. Rev. D 101, 055001 (2020).
  * McCutchan [2015] E. McCutchan, Nuclear Data Sheets 126, 151 (2015).
  * Nesterenko  _et al._ [2019] D. Nesterenko, L. Canete, T. Eronen, A. Jokinen, A. Kankainen, Y. Novikov, S. Rinta-Antila, A. de Roubin, and M. Vilen, International Journal of Mass Spectrometry 435, 204 (2019).
  * Eronen  _et al._ [2012] T. Eronen _et al._ , Eur. Phys. J. A 48, 46 (2012).
  * Moore  _et al._ [2013] I. Moore _et al._ , Nucl. Instrum. Methods Phys. Res. B 317, 208 (2013).
  * Uddin  _et al._ [2004] M. Uddin, M. Hagiwara, N. Kawata, T. Itoga, N. Hirabayashi, M. Baba, F. Tarkanyi, F. Ditroi, and J. Csikai, Journal of Nuclear Science and Technology 41, 160 (2004).
  * Koning  _et al._ [2019] A. Koning, D. Rochman, J.-C. Sublet, N. Dzysiuk, M. Fleming, and S. van der Marck, Nuclear Data Sheets 155, 1 (2019), special Issue on Nuclear Reaction Data.
  * Karvonen  _et al._ [2008] P. Karvonen, I. Moore, T. Sonoda, T. Kessler, H. Penttilä, K. Peräjärvi, P. Ronkanen, and J. Äystö, Nucl. Instrum. Meth. Phys. Res. B 266, 4794 (2008).
  * Nieminen  _et al._ [2001] A. Nieminen, J. Huikari, A. Jokinen, J. Äystö, P. Campbell, and E. Cochrane, Nucl. Instrum. Meth. Phys. Res. A 469, 244 (2001).
  * Savard  _et al._ [1991] G. Savard, S. Becker, G. Bollen, H. J. Kluge, R. B. Moore, T. Otto, L. Schweikhard, H. Stolzenberg, and U. Wiess, Phys. Lett. A 158, 247 (1991).
  * Eliseev  _et al._ [2013] S. Eliseev, K. Blaum, M. Block, C. Droese, M. Goncharov, E. Minaya Ramirez, D. A. Nesterenko, Y. N. Novikov, and L. Schweikhard, Phys. Rev. Lett. 110, 082501 (2013).
  * Eliseev  _et al._ [2014] S. Eliseev _et al._ , Applied Physics B 114, 107 (2014).
  * Nesterenko  _et al._ [2018] D. A. Nesterenko _et al._ , Eur. Phys. J. A 54, 154 (2018).
  * Gabrielse [2009] G. Gabrielse, Phys. Rev. Lett. 102, 172501 (2009).
  * Birge [1932] R. T. Birge, Phys. Rev. 40, 207 (1932).
  * Nesterenko  _et al._ [2021] D. Nesterenko, T. Eronen, Z. Ge, A. Kankainen, and M. Vilen, The European Physical Journal A 57, 302 (2021).
  * Roux  _et al._ [2013] C. Roux _et al._ , The European Physical Journal D 67, 146 (2013).
  * Wang  _et al._ [2021] M. Wang, W. Huang, F. Kondev, G. Audi, and S. Naimi,  45, 030003 (2021).
  * Carroll  _et al._ [2001] J. Carroll, S. Karamian, L. Rivlin, and A. Zadernovsky, Hyperfine Interactions 135, 3 (2001).
  * Lehnert  _et al._ [2020] B. Lehnert, H. Ramani, M. Hult, G. Lutter, M. Pospelov, S. Rajendran, and K. Zuber, Phys. Rev. Lett. 124, 181802 (2020).