# Parameterization of Stochasticity in Galaxy Clustering and Reconstruction of
Tomographic Matter Clustering

Shuren Zhou 1,2<EMAIL_ADDRESS>Pengjie Zhang 1,3,2<EMAIL_ADDRESS>1 Department of Astronomy, School of Physics and Astronomy, Shanghai Jiao Tong
University, Shanghai, China
2 Key Laboratory for Particle Astrophysics and Cosmology (MOE)/Shanghai Key
Laboratory for Particle Physics and Cosmology, China
3 Division of Astronomy and Astrophysics, Tsung-Dao Lee Institute, Shanghai
Jiao Tong University, Shanghai, China

###### Abstract

The stochasticity in galaxy clustering, the mismatch between galaxy and
underlying matter distribution, suppresses the matter clustering amplitude
reconstructed by the combination of galaxy auto-correlation and galaxy-galaxy
lensing cross-correlation. In this work, we solve the stochasticity
systematics by parameterizing the cross correlation coefficient $r(k)$ between
galaxy and matter. We investigate the performance of 12 kinds of
parameterization schemes, against the cosmoDC2 $\&$ TNG300-1 galaxy samples
over a wide range of redshift and flux cut. The 2-parameter fits are found to
describe the stochasticity up to $k_{\rm max}=0.9\,{\rm Mpc^{-1}}h$, while the
best performing quadratic scheme $r^{2}_{s}(k)=1+c_{1}k+c_{2}k^{2}$ reaches
better than $1\%$ accuracy for both the direct ${r}^{2}_{s}(k)$ fit and
reconstructing matter clustering. Then, we apply the accurate quadratic scheme
to forecast the tomographic matter clustering reconstruction by the
combination DESI-like LRG $\times$ CSST-like cosmic shear. Depending on
assumption of stochasticity, we find that the neglect of a serious
stochasticity would result in significant systematic bias in both the
reconstruction and the inferred cosmological parameters, even if we adopt
scale cut $k_{\rm max}=0.1\,{\rm Mpc^{-1}}h$. We demonstrate the necessity of
including stochasticity in reconstruction, and forecast that the
reconstruction alone enables a $S_{8}$ constraint at about $1.5\%$ precision,
free from galaxy bias and stochasticity. We will validate our method for DESI
spectroscopic survey, and the analysis is expected to be complementary to DESI
cosmological constraint by BAO and RSD.

##  I Introduction

Galaxy clustering contains a wealth of cosmological information, and serves as
an important probe to measure the large scale structure in the ongoing
spectroscopic galaxy surveys such as DESI [1, 2, 3] and LSST [4]. The upcoming
million spectroscopic objects offer precision but biased matter clustering
measurement. Meanwhile, other cosmological observables, such as galaxy lensing
delivered from Euclid satellite mission [5] and the future CSST cosmic shear
survey [6, 7, 8], is complementary to the galaxy clustering. The gravitational
lensing is the projected distortion effect of the background object due to the
deflection of arrival photons by the foreground matter clustering, providing
the unbiased but integrated matter clustering amplitude. The joint analysis of
both two, specifically, galaxy auto-power spectrum ${C}^{gg}_{\ell}$ and
cross-power spectrum with lensing convergence ${C}^{g\kappa}_{\ell}$, can
break the degeneracy between galaxy bias and matter clustering amplitude,
obtaining the tomographic matter clustering power spectrum [9, 10, 11, 12].

The galaxy deterministic bias can be cancelled model-independently combining
different data set (e.g. the gravitational potential decay rate measurement
[13, 14] and $E_{G}$ statistic [15]). This would motivate a simple estimator
for the matter power spectrum,

$\displaystyle\hat{P}_{mm}(k)=\hat{P}_{gm}^{2}(k)/\hat{P}_{gg}(k)\quad.$ (1)

Here $\hat{P}_{gg}$ and $\hat{P}_{gm}$ are the galaxy-galaxy and galaxy-matter
3D power spectrum, realized by the angular power spectrum $\hat{C}^{gg}$ and
$\hat{C}^{g\kappa}$ with narrow redshift width. It seems a promising estimate
of matter clustering, since any galaxy bias $b(k)$ spontaneously disappears if
the dependence of galaxy on matter distribution is deterministic. However,
apart from the deterministic galaxy bias, there are stochastic components in
the galaxy-matter relation, arising from the intrinsic discrete property and
the complex physical processes in halo collapse and galaxy formation [16, 17,
18, 19]. Therefore, the ensemble average of the estimator is biased by a
multiplicative $r^{2}$,

$\displaystyle\langle\hat{P}_{mm}(k)\rangle=r^{2}(k)P_{mm}(k)\quad.$ (2)

Here $P_{mm}$ is the theoretical expectation, and $r^{2}$ is the cross
correlation coefficient between galaxy number counts overdensity
$\delta_{g}({\bf k})$ and underlying matter overdensity $\delta_{m}({\bf k})$.
The function $r^{2}(k)$ quantifies the efficiency of galaxy tracing the matter
distribution, and we always have $r(k)<1$ for stochastic tracers. The neglect
of the impact of stochasticity in galaxy clustering will suppress the
reconstructed matter power spectrum by $r^{2}(k)$, equivalently suppressing
the amplitude by $r(k)$, depending on specific scale $k$ we interest in. The
overall suppression also propagates to the cosmological parameters constraint
and results in potentially systematics. In order to leverage the unbiased
clustering information in non-linear scale, which is contaminated by
stochasticity, we expect an accurate $r(k)$ parameterization to describe the
suppression in estimator Eq. (1).

Two pathways have been explored to alleviate the stochasticity induced
systematics. One pathway is to realize the reduced stochasticity by some
proper weights. The weighting schemes are adopted by weighing $\delta_{g}$
with galaxy bias [17], function of host halo mass [18, 20, 21, 22] or some
principal components [23, 24, 25]. Especially, the principal component
analysis (PCA) suggests that the principal eigenmode of galaxy clustering
traces the matter clustering, while only two eigenmodes significantly deviate
from the shotnoise expectation [24, 25]. It implies that the effective degree
of freedom (${\rm DoF}$) in stochasticity is finite and approximately ${\rm
DoF}\simeq 2$. Another pathway is to parameterize the stochasticity with some
phenomenological parameters. In the perspective of effective field theory of
large scale structure (EFTofLSS), the stochasticity contribution to the galaxy
auto-/cross-power spectrum is expanded in even power of Fourier mode $k$ [16,
26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37]. In the widely adopted scheme,
the expansion series in real space is truncated at the first two terms, where
the first constant term characterizes the deviation of Poisson shot noise and
the second $k^{2}$-scaled term describes the scale dependence. In the redshift
space, an additional free parameter is required to describe the anisotropic
directional dependence [35, 37]. These phenomenological descriptions are fully
investigated in the content of galaxy bias expansion scheme, but it is
difficult to ascertain the impact of stochasticity because of the mixture with
higher order terms omitted in expansion truncation [28, 29, 31]. Moreover,
perturbation description only applies to the scale $k\lesssim 0.3\,{\rm
Mpc^{-1}}h$ due to the strong non-linearity [38, 28, 29, 30], losing modeling
ability beyond linear region.

In this work, we aim to seek a parameterization for the cross correlation
coefficient $r(k)$, following the perspective of EFTofLSS. In the following
context, we investigate various parameterization forms and test their
performance against the galaxy samples in simulations. The simple quadratic
scheme $r_{s}^{2}=1+c_{1}k+c_{2}k^{2}$ is identified as the most promising
parameterization. Then, we demonstrate the reconstruction of the matter power
spectrum combining galaxy clustering and galaxy lensing, together with the
best performing quadratic $r(k)$ scheme. Our method is suited for the
tomographic clustering analysis in DESI Luminous Red Galaxies (LRG) samples.
DESI was designed to conduct a survey covering $14000\,{\rm deg^{2}}$
footprint, obtaining about $40$ millions of spectroscopic targets ranging from
redshift 0 to $3.5$. These precision measurements will enable the constraint
of $f\sigma_{8}$ to an aggregate precision $0.95\%$ by combining the baryonic
acoustic oscillations (BAO) and redshift space distortion (RSD) [2], where $f$
is the growth and $\sigma_{8}$ is the amplitude of matter fluctuation in
spheres of $8\,h^{-1}{\rm Mpc}$. Our analysis enables the direct constraint on
fluctuation amplitude $\sigma_{8}$, which could serve as a strong supplement
for DESI’s fiducial analysis on structure growth. It is also an immediate
comparison with the results from cosmic microwave background (CMB) constraint
[39], to help to clarify whether the anomaly in low redshift observations
compared to CMB is attributed to new physics or merely systematic
contamination [40, 10, 11].

The paper is organized as followed. In Section II, we introduce the
theoretical framework under the perturbation expansion perspective to
illustrate the impact of stochasticity on the galaxy power spectrum, and then
introduce our parameterization schemes. Readers with limited interest in the
background can skip Section II and refer to Table 1 directly. In Section III,
we introduce the methodology to test the performance of the parameterization
and the methodology combining the galaxy clustering and galaxy lensing. In
Section IV, we present results of the performance test. In Section V, we
present the fisher forecast for the combination of DESI-like spectroscopic
tracers and CSST-like cosmic shear.

##  II Parameterization Scheme

Table 1:  Parameterization expressions of $r^{2}(k)$ and $r^{2}_{s}(k)$, defined in Eq. (4) and Eq. (5). Notation $r^{2}_{(s)}$ refers to the parameterization that are applied to both $r^{2}$ and $r^{2}_{s}$. $P_{mm}$ and $P_{gg,s}$ are the matter power spectrum and the galaxy power spectrum subtracted the Poisson expectation. $\Delta^{2}(k)\equiv k^{3}P(k)/(2\pi^{2})$ is the corresponding variance.  | Label | Formula | Free Parameters |
---|---|---|---|---
| Expan-mm-1 | $r^{2}_{(s)}=\left[1+\alpha\,/\left(\bar{n}_{g}P_{mm}\right)\right]^{-1}$ | $\;\alpha\;$ |
| Expan-mm-2 | $r^{2}_{(s)}=\left[1+\alpha\,(1+\beta k^{2})/\left(\bar{n}_{g}P_{mm}\right)\right]^{-1}$ | $\;\alpha,\;\beta\;$ |
| Expan-mm-3 | $r^{2}=\left[1+\alpha\,(1+\beta k^{2}+\gamma k^{4})/\left(\bar{n}_{g}P_{mm}\right)\right]^{-1}$ | $\;\alpha,\;\beta,\;\gamma\;$ |
| Expan-gg-1 | $r^{2}=\left[1+\alpha\,/\left(\bar{n}_{g}P_{gg,s}\right)\right]^{-1}$ | $\;\alpha\;$ |
| Expan-gg-2 | $r^{2}=\left[1+\alpha\,(1+\beta k^{2})/\left(\bar{n}_{g}P_{gg,s}\right)\right]^{-1}$ | $\;\alpha,\;\beta\;$ |
| Expan-gg-3 | $r^{2}=\left[1+\alpha\,(1+\beta k^{2}+\gamma k^{4})/\left(\bar{n}_{g}P_{gg,s}\right)\right]^{-1}$ | $\;\alpha,\;\beta,\;\gamma\;$ |
| Expan-new-2 | $r^{2}_{(s)}=(1+\beta k^{2})\left[1+\alpha/\left(\bar{n}_{g}P_{mm}\right)\right]^{-1}$ | $\;\alpha,\;\beta\;$ |
| Zheng2013-mm | $r^{2}_{(s)}=\left(1+\beta\Delta^{2}_{mm}\right)^{-1}\left[1+\alpha/\left(\bar{n}_{g}P_{mm}\right)\right]^{-1}$ | $\;\alpha,\;\beta\;$ |
| Zheng2013-gg | $r^{2}=\left(1+\beta\Delta^{2}_{gg,s}\right)^{-1}\left[1+\alpha/\left(\bar{n}_{g}P_{gg,s}\right)\right]^{-1}$ | $\;\alpha,\;\beta\;$ |
| Q-bias | $r_{s}^{2}=\left(1+c_{2}k^{2}\right)/\left(1+c_{1}k\right)$ | $\;c_{1},\;c_{2}\;$ |
| Quadratic | $r_{s}^{2}=1+c_{1}k+c_{2}k^{2}$ | $\;c_{1},\;c_{2}\;$ |
| Quadratic-1 | $r_{s}^{2}=1+c_{2}k^{2}$ | $\;c_{2}\;$ |

In Fourier space, the galaxy overdensity $\delta_{g}(\rm x)$ can be decomposed
into the deterministic and the stochastic part

$\delta_{g}({\bf k})=b(k)\delta_{m}({\bf k})+\delta_{\mathcal{S}}({\bf k})$
(3)

Here the scale dependent galaxy bias is defined as $b(k)\equiv
P_{gm}(k)/P_{mm}(k)$, therefore $\langle\delta_{m}({\bf
k})\delta_{\mathcal{S}}^{*}({\bf k})\rangle=0$. In the galaxy auto power
spectrum, the stochasticity induces an additionally inconsistent term, that
$P_{gg}(k)=b^{2}(k)P_{mm}(k)+P_{\mathcal{S}}(k)$. We can quantify the
efficiency of the galaxy tracing the matter field by the cross correlation
coefficient

$\displaystyle r^{2}(k)$ $\displaystyle\equiv\;$
$\displaystyle\frac{P_{gm}^{2}}{P_{mm}P_{gg}}=\frac{1}{1+P_{\mathcal{S}}\left(b^{2}\,P_{mm}\right)^{-1}}\quad.$
(4)

Another widely adopted convention is

$\displaystyle r_{s}^{2}(k)\equiv\frac{P_{gm}^{2}}{P_{mm}\,P_{gg,s}}\quad,$
(5)

where the Poisson noise expectation is subtracted from the galaxy power
spectrum, $P_{gg,s}(k)=P_{gg}(k)-1/\bar{n}_{g}$. In the following context, we
utilize both $r^{2}$ and $r^{2}_{s}$ to describe stochasticity since they are
equivalent in practical applications.

In the perspective of EFTofLSS, a complete description of the general bias
expansion requires the stochastic variables due to the coupling of the small
scale modes, statistic properties of which in general depend on the local
observables [41, 42, 16, 38]. Up to the third order, the analogous expansion
of the stochastic components in Eq. (3) is

$\displaystyle\delta_{\mathcal{S}}({\bf x})$ $\displaystyle=\,$
$\displaystyle\varepsilon_{0}({\bf x})+\varepsilon_{1}({\bf x})\delta_{m}({\bf
x})+\varepsilon_{2}({\bf x})\delta_{m}^{2}({\bf x})+\varepsilon_{s^{2}}({\bf
x})s^{2}({\bf x})$ (6) $\displaystyle+\,$ $\displaystyle
b_{s^{2}}\left[s^{2}_{\mathcal{S}}({\bf x})-\langle s^{2}_{\mathcal{S}}({\bf
x})\rangle\right]+b_{s^{3}}\left[s^{3}_{\mathcal{S}}({\bf x})-\langle
s^{3}_{\mathcal{S}}({\bf x})\rangle\right]$ $\displaystyle+\,$ $\displaystyle
b_{\delta s^{2}}\,\delta_{m}s^{2}_{\mathcal{S}}({\bf
x})+\,b_{st}\left[st_{\mathcal{S}}({\bf x})-\langle st_{\mathcal{S}}({\bf
x})\rangle\right]$ $\displaystyle+\,$ $\displaystyle
b_{\psi}\left[\psi_{\mathcal{S}}({\bf x})-\langle\psi_{\mathcal{S}}({\bf
x})\rangle\right]\;+\,\mathcal{O}(\varepsilon^{4})\quad,$

where the tidal field and other higher order terms are defined as [41]

$\displaystyle s_{ij}({\bf x})$ $\displaystyle=$
$\displaystyle\left(\nabla_{i}\nabla_{j}\nabla^{-2}-\frac{1}{3}\delta_{ij}^{K}\right)\delta({\bf
x})$ $\displaystyle t_{ij}({\bf x})$ $\displaystyle=$
$\displaystyle\left(\nabla_{i}\nabla_{j}\nabla^{-2}-\frac{1}{3}\delta_{ij}^{K}\right)\left[\theta({\bf
x})-\delta({\bf x})\right]$ $\displaystyle\psi({\bf x})$ $\displaystyle=$
$\displaystyle\theta({\bf x})-\delta({\bf x})-\frac{2}{7}s^{2}({\bf
x})+\frac{4}{21}\delta^{2}({\bf x})\quad.$ (7)

Here the stochastic variables
$\varepsilon_{0},\varepsilon_{1},\varepsilon_{2}\,\cdots$ are uncorrelated
with the large scale perturbations, completely described by their own joint
moments. The multiplicative constants $b_{s^{2}},\,b_{s^{3}},\,b_{\delta
s^{2}},\,b_{st}$ and $b_{\psi}$ are higher order bias, and the corresponding
scalars are contracted by $s^{2}=s_{ij}s_{ji}$, $s^{3}=s_{ij}s_{jk}s_{ki}$ and
$st=s_{ij}t_{ji}$. In Eq. (6), different from the usual galaxy bias expansion,
we have to further decompose the tidal field
$s^{2}=s^{2}_{\mathcal{D}}+s^{2}_{\mathcal{S}}$, where $s^{2}_{\mathcal{D}}$
is completely correlated with matter field and resummed to the
$b(k)\delta_{m}({\bf k})$ term, while $s^{2}_{\mathcal{S}}$ serves as one of
the stochastic components with $\langle
s^{2}_{\mathcal{S}}\,\delta_{m}\rangle=0$. The same arguments are applied to
components $s^{3}$, $\delta_{m}s^{2}$, $st$ and $\psi$. Taking the Fourier
transform of Eq. (6), we obtain the stochasticity power spectrum,

$\displaystyle P_{\mathcal{S}}(k)$ $\displaystyle=\;$ $\displaystyle
P_{\varepsilon_{0}\varepsilon_{0}}(k)+\int\frac{d^{3}{\bf
q}}{(2\pi)^{3}}\,P_{\varepsilon_{1}\varepsilon_{1}}(k)P_{mm}(|{\bf q-k}|)$ (8)
$\displaystyle+\,$
$\displaystyle\frac{13}{9}\sigma^{2}_{m}P_{\varepsilon_{0}\varepsilon_{2}}(k)+b_{s^{2}}^{2}P_{s^{2},\mathcal{S}}(k)+\mathcal{O}(\varepsilon^{5})\quad,$

where the RMS of the fluctuation of the matter field is defined by
$\sigma^{2}_{m}\equiv\int\,P_{mm}(q)q^{2}dq\,/(2\pi^{2})$.

Eq. (6) shows that, although $\delta_{\mathcal{S}}({\bf x})$ does not
contribute to cross power spectrum, it still contributes in the higher order
statistics in a complex manner. We can classify these terms in Eq. (6) into
three types. The first type is the isolated stochasticity parameter
$\varepsilon_{0}$, which would contribute common shot noise
$P_{\varepsilon_{0}\varepsilon_{0}}\rightarrow 1/\bar{n}_{g}$ for Poisson
sampling. The second type is the coupling term between the stochastic
variables $\varepsilon_{1},\varepsilon_{2},\varepsilon_{s^{2}}$ and the large
scale fluctuation field, reflecting the coevolution of gravity and galaxy
stochasticity [16]. The coevolution also appears in the coupling terms in
power spectrum Eq. (8), in the form of four points correlation or higher order
correlation. The third type is the residual decorrelation of the underlying
large scale fluctuation field with respect to the matter field, and these
terms do not couple the stochastic variables in the power spectrum.

The parameterization of $r^{2}(k)$ is equivalent to parameterize the
stochasticity power spectrum Eq. (8). Though some stochasticity effects have
been investigated analytically, i.e. Ref. [19] simplify the dark matter halo
as hard sphere to derive the halo exclusion effect, it is still unable to
model galaxy stochasticity in a comprehensive way due to the complex nature of
galaxy formation and evolution. In order to construct an accurate
parameterization for $r^{2}(k)$, we draw inspiration from the usual
methodology in galaxy bias expansion, where the leading order stochasticity
power spectrum $P_{\varepsilon_{0}\varepsilon_{0}}(k)$ is expanded as a series
of $k$ and keep only the first two terms [16, 26, 27, 28, 29, 30, 31, 37].
Therefore, the stochasticity power spectrum is parameterized as

$P_{\mathcal{S}}(k)=\frac{\alpha}{\bar{n}_{g}}(1+\beta k^{2}+\gamma
k^{4})+\mathcal{O}(k^{6})\quad.$ (9)

Here the terms in odd power of $k$ vanish due to the isotropy of the
stochasticity field. The multiplication constant $\alpha/\bar{n}_{g}$ accounts
for the Poisson-like noise, and the $k$ series in bracket describe the scale
dependence of stochasticity. As mentioned above, the previous PCA
investigations have strong implication that the effective ${\rm DoF}\simeq 2$,
namely two free parameters are sufficient to characterize stochasticity.
Therefore, we propose the first parameterization formula in article,

$r^{2}(k)=\;\frac{1}{1+\alpha\left(1+\beta
k^{2}\right)/\left(\bar{n}_{g}P_{gg,s}(k)\right)}\quad,$ (10)

where $(\alpha,\beta)$ are free parameters and observables $P_{gg,s}$ replace
$b^{2}P_{mm}$ compared to definition Eq. (4). Additionally, we keep $\alpha$
as the only non-zero free parameter to investigate the ability of 1-parameter
fit describing stochasticity, and also test the 3-parameter fit to investigate
the robustness of the assumption ${\rm DoF}\simeq 2$.

We can replace the measurable $P_{gg}$ with theoretical $P_{mm}$, which is
determined simultaneously during the matter clustering modeling and
cosmological parameters constraint (see Expan-mm-2 in Table 1), and now the
parameterization is free from uncertainty in $P_{gg}$ measurement. Compared to
Eq. (10), $\alpha$ is roughly rescaled by the square of linear galaxy bias,
and the quadratic term absorbs the scale dependence of galaxy bias to some
extent. In large scale limit, the absorption of the scale dependence by
quadratic term is expected, since the leading order scale dependence arising
from $\nabla^{2}\delta$ contributes to $b(k)$ scales as $k^{2}$ [41]. Another
straightforward thinking is whether the direct expansion of $r^{2}(k)$ with
respect to $k$ is better than the expansion of $P_{\mathcal{S}}$. So we expand
the scale dependence of $r^{2}$ as $k$ series and keep the first correction
term, then obtain a new parameterization form (see Expan-new-2 in Table 1).

Another class of parameterization schemes inspired by the modeling of galaxy
peculiar velocity in Refs. [43, 44], in which the ratio between the velocity
divergence component $\theta$ and the matter overdensity $\delta_{m}$ is
statistically compacted as a window function $W(k)$, proportional to their
cross correlation coefficient $r_{\theta\delta}$. $W(k)$ describes both the
non-linear evolution and the stochastic component between two fields
$(\theta,\delta_{m})$, similar to $(\delta_{g},\delta_{m})$ relation in our
case. Further its fitting formula motivated by perturbation theory achieves
accuracy within $2\%$ in mildly non-linear region [44]. With these
considerations, we adopt the proposed parameterized form of $W(k)$ for the
$r^{2}(k)$ fit, and multiply the Poisson-like term in Eq. (10) to mimic the
suppression due to the intrinsic discrete noise.

$r^{2}(k)\;=\;\frac{1}{1+\beta\Delta_{mm}^{2}}\,\frac{1}{1+\alpha/\left(\bar{n}_{g}P_{mm}\right)}\quad.$
(11)

Here $\Delta_{mm}^{2}(k)\equiv k^{3}P_{mm}(k)/(2\pi^{2})$ is the matter power
spectrum variance, and $(\alpha,\beta)$ are the free parameters. Similarly, we
can replace $P_{mm}$ with $P_{gg}$ and obtain a new parameterization.

The most straightforward parameterization scheme is to resum the $k$\- $\&$
$k^{2}$\- dependence of $r_{s}^{2}$ to quadratic expansion
$r_{s}^{2}=1+c_{1}k+c_{2}k^{2}$, where $r_{s}\rightarrow 1$ as $k\rightarrow
0$. Physically, we expected all the stochasticity are localized within scale
$R_{*}$, where $R_{*}$ is the typical scale of galaxy formation [16], because
galaxy distribution in very large scale is only determined by the large scale
fluctuation traced by matter field. However, it is not rigorously true,
because the constant discrete noise deviates from $1/\bar{n}_{g}$ by the sub-
Poisson suppression or super-Poisson enhancement. Nevertheless, in practical
cases, it does not prohibit the excellent performance of the formula such as
$r_{s}^{2}=1+c_{1}k+c_{2}k^{2}$. On the one hand, we expect deviation of
$r_{s}$ from $1$ is negligible, since the clustering signal dominates in large
scale and cosmic variance is significant. On the other hand, in realistic
case, the formula allows $r_{s}^{2}<1$ at the largest scale by the adjustment
of the $c_{1},c_{2}$ terms. Therefore, we expect the quadratic formula is
accurate if we set a small enough $k_{\rm max}$ during fitting. We also
further reduce the number of free parameters by keeping only the quadratic
terms (see Quadratic-1 in Table 1). Besides, based on the analysis of mock
galaxy catalogs, the non-linear scale dependence of galaxy bias is proposed to
be parameterized as $\left(1+c_{2}k^{2}\right)/\left(1+c_{1}k\right)$, known
as ’Q-bias’ prescription [45, 46, 47, 48]. We migrate it to our stochasticity
description.

The parameterization schemes are summarized in Table 1. In the following
context, we use the label shown in the first column of the table to refer to
these parameterization schemes. There are two distinctive types of
parameterization for $r^{2}$ and $r^{2}_{s}$ when applying to realistic galaxy
survey. One type of parameterization includes $P_{gg,s}$ that is directly
measurable in the realistic data. While another type of parameterization
includes $P_{mm}$ that is simultaneously determined in the theoretical
framework fitting. We omit this distinction below for the purpose of testing.

Table 2:  Galaxy samples: cosmoDC2 .

$\mathsf{FluxLimitLabel}$ | Flux-limit of $ugrizY$ bands | $n_{g}\;[\,{\rm Mpc^{-3}h^{3}}\,]$ |
---|---|---|---
z=0.15 | z=0.5 | z=1.0 | z=1.5 |
1 | ( 28.1, 29.4, 29.5, 28.8, 28.1, 26.9 ) | 0.0845 | 0.0592 | 0.0582 | 0.0332 |
2 | ( 27.1, 28.4, 28.5, 27.8, 27.1, 25.9 ) | 0.0823 | 0.0494 | 0.0458 | 0.0240 |
3 | ( 26.1, 27.4, 27.5, 26.8, 26.1, 24.9 ) | 0.0793 | 0.0394 | 0.0344 | 0.0167 |
4 | ( 25.1, 26.4, 26.5, 25.8, 25.1, 23.9 ) | 0.0749 | 0.0300 | 0.0256 | 0.0111 |
5 | ( 24.1, 25.4, 25.5, 24.8, 24.1, 22.9 ) | 0.0686 | 0.0225 | 0.0196 | 0.0067 |
6 | ( 23.1, 24.4, 24.5, 23.8, 23.1, 21.9 ) | 0.0609 | 0.0163 | 0.0140 | 0.0033 |
7 | ( 22.1, 23.4, 23.5, 22.8, 22.1, 20.9 ) | 0.0520 | 0.0115 | 0.0090 | 0.0014 |
8 | ( 21.1, 22.4, 22.5, 21.8, 21.1, 19.9 ) | 0.0428 | 0.0076 | 0.0057 | 0.0005 |

Table 3:  Galaxy samples: TNG300-1 .

$\mathsf{FluxLimitLabel}$ | Flux-limit of $griz$ bands | $n_{g}\;[\,{\rm Mpc^{-3}h^{3}}\,]$ |
---|---|---|---
z=0.15 | z=0.5 | z=1.0 | z=1.5 |
1 | ( 29.4, 29.5, 28.8, 28.1) | 0.1353 | 0.0915 | 0.0698 | 0.0550 |
2 | ( 28.4, 28.5, 27.8, 27.1) | 0.1221 | 0.0794 | 0.0587 | 0.0470 |
3 | ( 27.4, 27.5, 26.8, 26.1) | 0.1129 | 0.0700 | 0.0497 | 0.0377 |
4 | ( 26.4, 26.5, 25.8, 25.1) | 0.1032 | 0.0611 | 0.0422 | 0.0298 |
5 | ( 25.4, 25.5, 24.8, 24.1) | 0.0914 | 0.0525 | 0.0352 | 0.0244 |
6 | ( 24.4, 24.5, 23.8, 23.1) | 0.0815 | 0.0439 | 0.0295 | 0.0198 |
7 | ( 23.4, 23.5, 22.8, 22.1) | 0.0707 | 0.0379 | 0.0239 | 0.0163 |
8 | ( 22.4, 22.5, 21.8, 21.1) | 0.0608 | 0.0316 | 0.0195 | 0.0125 |

##  III Methodology

We present the accuracy analysis of the parameterization schemes in two ways:
(1) Quantify absolute difference between the measurable $\hat{r}$ and the fit
values $r$ by quantity ${\bf Q}$. (2) Investigate the impact on the matter
clustering amplitude $A$ by the parameterization schemes, where the accuracy
and precision of $A$ is quantified by the systematic bias $\delta_{A}$ and
statistic uncertainty $\sigma_{A}$.

We measure $\hat{r}^{2}$ in the simulations, and then fit with the
parameterized formula $r^{2}$ by minimizing $\chi^{2}$.

$\displaystyle\chi^{2}=\sum_{k=k_{\rm min}}^{k_{\rm
max}}{\frac{\left[\hat{r}^{2}(k)-r^{2}(k)\right]^{2}}{\sigma^{2}(k)}}\quad,$
(12)

where $\sigma^{2}$ is approximated as the Gaussian variance of the measured
$r^{2}$ or $r^{2}_{s}$. In principle, both the matter and galaxy overdensity
field are non-Gaussian, especially in the small scale region where the strong
gravitational clustering leads to non-linear evolution. The Gaussian variance
only accounts for part of the fluctuation in overdensity field, and the
underestimation of Gaussian variance increase severely as deep into non-linear
scale. So $\chi^{2}$ defined by Eq. (12) does not accurately characterize the
statistic significance of the parameterization fitting with measurement.
However, it still serves as a practical criterion to fit the data and find the
bestfit parameters. To avoid the computational challenges of estimating large
amount of covariance of test samples, we utilize the absolute difference of
$r_{(s)}$ to quantify the goodness of fit, rather than the $\chi^{2}$
statistic. We define the absolute deviation between the measured $\hat{r}$ and
$r$ over the fit range $(k_{\rm min},\,k_{\rm max})$

$\displaystyle{\bf
Q}\;=\;\sqrt{\frac{\sum_{i}{N_{i}\,\left(\hat{r}_{i}-r_{i}\right)^{2}}}{\sum_{i}N_{i}}}\quad,$
(13)

here $N_{i}\propto k^{2}_{i}\Delta k_{i}$ is the number of independent $k$
modes in $i$-th $k$ bin, and ${\bf Q}$ is the RMS over all the fitting bin
weighted by the number of $k$ modes in each bin.

In the cosmological applications, we are interested in the impact on the
accuracy of the reconstructed matter power spectrum, rather than the goodness
of fit. Specifically, we focus on two questions: how much the bias in the
constrained total amplitude of the matter power spectrum (or equivalently
$\sigma_{8}$ in linear scale) is and how much the constraint precision should
be sacrificed using these parameterizations. As an illustration, we forecast
these impacts by combining the projected observable $\hat{C}^{g\kappa}_{\ell}$
and $\hat{C}^{gg}(k)$, where the galaxy is the spectroscopic tracer that
targeted by DESI [1, 2], and the lensing convergence $\kappa$ could be
realized by the future CSST cosmic shear survey [6, 7, 8].

Under the Limber’s approximation, the measured angular power spectrums of
galaxy clustering and galaxy lensing are given by

$\displaystyle\hat{C}^{g\kappa}_{\ell}$ $\displaystyle=$
$\displaystyle\int{d\chi\,\chi^{-2}\,n_{g}(\chi)W_{L}(\chi)\hat{P}_{gm}\left(\frac{\ell}{\chi},\chi\right)}$
(14) $\displaystyle\hat{C}^{gg}_{\ell}$ $\displaystyle=$
$\displaystyle\int{d\chi\,\chi^{-2}\,n_{g}(\chi)n_{g}(\chi)\hat{P}_{gg}\left(\frac{\ell}{\chi},\chi\right)}\quad.$

Here $W_{L}(\chi)$ is the lensing kernel. For simplicity, we assume the both
the redshift distributions of the foreground galaxy are narrow enough that
$n_{g}(z)\rightarrow\delta^{D}(z-z_{g})$, then we have

$\displaystyle\hat{C}^{g\kappa}_{\ell}$ $\displaystyle=$
$\displaystyle\chi_{g}^{-2}\,W_{L}(\chi_{g})\,\hat{P}_{gm}\left(\frac{\ell}{\chi_{g}}\right)\quad,$
(15) $\displaystyle\hat{C}^{gg}_{\ell}$ $\displaystyle=$
$\displaystyle\chi_{g}^{-2}\,\Delta\chi_{g}^{-1}\,\hat{P}_{gg}\left(\frac{\ell}{\chi_{g}}\right)\quad.$
(16)

Here $\Delta\chi_{g}$ is the comoving radical distance width of the galaxy
bin. The reconstructed power spectrum $\hat{P}(k)$ differs from the direct
measurable $\hat{C}_{\ell}$ with a geometric factor. We can rewrite
$\hat{P}_{gg}$ in Eq. (16) in the form of cross power spectrum.

$\displaystyle C^{g\kappa}_{\ell}$ $\displaystyle=$
$\displaystyle\chi_{g}^{-2}\,W_{L}(\chi_{g})\,r\,\sqrt{\hat{P}_{gg}\left(\frac{\ell}{\chi_{g}}\right)P_{mm}\left(\frac{\ell}{\chi_{g}}\right)}$
(17)

Here $r$ is the model of cross correlation. If we subtract the Poisson
expectation in the galaxy power spectrum, a replacement should be made
$(r\,,\hat{P}_{gg})\rightarrow(r_{s}\,,\hat{P}_{gg,s})$. Since the fluctuation
in cosmic shear is the dominant noise source in the realistic measurements, we
can neglect the uncertainty in $\hat{P}_{gg}$ and consider only the noise
contribution from $\hat{C}^{g\kappa}_{\ell}$. Therefore we can treat
$C^{g\kappa}_{\ell}$ in Eq. (17) as a theoretical template for the observable
of Eq. (15). The matter power spectrum $P_{mm}$ can be estimated by minimizing

$\displaystyle\chi^{2}_{A}=\left[\hat{C}^{g\kappa}_{\ell}-C^{g\kappa}_{\ell}\right]^{T}{\bf
Cov}^{-1}\left[\hat{C}^{g\kappa}_{\ell}-C^{g\kappa}_{\ell}\right]\;.$ (18)

Further, we only consider the Gaussian uncertainty, which leads to null off-
diagonal elements of covariance ${\bf Cov}$ and only its diagonal elements
contribute. Under these assumptions, the loss function Eq. (18) combined with
Eq. (15) and Eq. (17) can be simplified as

$\displaystyle\chi_{A}^{2}=\left[\frac{W_{L}(\chi_{g})}{\chi_{g}^{2}}\right]^{2}\sum_{\ell}\frac{\hat{P}_{gg}P_{mm}}{\sigma_{\ell}^{2}}\left(\hat{r}-A\,r\right)^{2}$
(19)

where we make replacement $P_{mm}\rightarrow A^{2}P_{mm}$ to fix $P_{mm}$ as
the fiducial true value. The free rescaling factor $A$ marginalizes the
overall amplitude of the measured matter field, $A=1$ meaning the unbiased
measurement. $\hat{r}\equiv\hat{P}_{gm}/\sqrt{\hat{P}_{gg}P_{mm}}$ is the
fiducial true cross correlation coefficient by definition. Now all the gained
cosmological information are contained in the amplitude $A$, of which the
systematic bias is characterized by $\delta_{A}=A-1$ and the statistic
uncertainty is $\sigma_{A}$. For a stochasticity parameterization scheme
given, the matter clustering amplitude as well as the nuisance stochasticity
parameters $\lambda_{\mu}=\left(\,A,\alpha,\beta,\cdots\right)$ are given by
the bestfit values of Eq. (19), then the systematic bias $\delta_{A}$ is also
determined. The statistic uncertainty $\sigma_{A}$ is given by the Fisher
analysis [49],

$\displaystyle\sigma^{2}_{\lambda_{\mu}}\;=\;\left(\frac{1}{2}\frac{\partial^{2}\chi_{A}^{2}}{\partial\lambda_{\alpha}\partial\lambda_{\beta}}\biggl{|}_{\rm
bestfit}\right)^{-1}\biggl{|}_{\mu\mu}$ (20)

We have assumed the Gaussian variance for measurement in Eq. (19), which is
given by

$\displaystyle\sigma^{2}_{\ell}$ $\displaystyle=$
$\displaystyle\frac{1}{f_{\rm
sky}(2\ell+1)\Delta\ell}\left[\hat{C}^{gg}_{\ell}\hat{C}^{\kappa\kappa}_{\ell}+(C^{g\kappa}_{\ell})^{2}\right]\quad.$
(21)

Specifically,

$\displaystyle\hat{C}^{gg}_{\ell}$ $\displaystyle=$
$\displaystyle{C}^{gg}_{\ell}+\frac{1}{\Sigma_{g}}\;,$
$\displaystyle\hat{C}^{\kappa\kappa}_{\ell}$ $\displaystyle=$
$\displaystyle{C}^{\kappa\kappa}_{\ell}+\frac{\sigma_{\gamma}^{2}}{\Sigma_{\gamma}}+N^{\gamma}_{\rm
add}\;,$ (22)

where, $\Sigma$ is the surface number density of tracers, $\sigma_{\gamma}$ is
the shape noise and $N_{\rm add}^{\gamma}$ is the additive noise term, all of
which vary with survey details.

If the impact of $r^{2}$ is neglected, namely setting a model $r=1$, the
systematic bias and the statistic uncertainty can be analytically expressed as

$\displaystyle\delta_{A}$ $\displaystyle\,=\,$
$\displaystyle\frac{\sum_{\ell}{\sigma_{\ell}^{-2}\,\hat{P}_{gg}P_{mm}\hat{r}}}{\sum_{\ell}{\sigma_{\ell}^{-2}\,\hat{P}_{gg}P_{mm}}}-1$
(23) $\displaystyle\sigma_{A}$ $\displaystyle\,=\,$
$\displaystyle\left[\frac{W_{L}(\chi_{g})}{\chi_{g}^{2}}\right]^{-1}\left(\sum_{\ell}\frac{\hat{P}_{gg}P_{mm}}{\sigma_{\ell}^{2}}\right)^{-1/2}$

Figure 1:  The measured $\hat{r}^{2}$ and their fitting results $r^{2}$ by Expan-mm-2 and Zheng2013-mm on cosmoDC2 galaxies ($\mathsf{FluxLimitLabel}=3$). Four groups of rows show the results at $z=0.15,0.5,1.0$ and $1.5$. The first three columns show the results of 3 selected galaxy $g-r$ color bins (total 6 bins), and the last column shows the result of flux limited full samples. Thick Rows: The solid line is $\hat{r}^{2}$ measured in simulation, and the dash or dot lines are the result of fitting formula Expan-mm-2 fitted with $k_{\rm max}=0.3,0.5,0.7$ and $0.9\,{\rm Mpc^{-1}}h$. Narrow Rows: They show the ratio between $r^{2}$ and $\hat{r}^{2}$, and the blue lines are fitted with Expan-mm-2 while the red lines are fitted with Zheng2013-mm, where horizontal lines mark the region of $\pm 1\%$. Notice that $\chi^{2}/N$ do not serve as a goodness of fit criterion here.  Figure 2:  Same plot as Fig. 1 for cosmoDC2 galaxies ($\mathsf{FluxLimitLabel}=3$), but here we present results for $r^{2}_{s}$. The blue lines are fitted with Quadratic while the red lines are fitted with $\texttt{Expan-mm-2}\;(r_{s}^{2})$.  Table 4:  Summarized Index for Stochasticity Parameterization. Two values provided for each index shown in form $A\,|\,B$, meaning the values measured in cosmoDC2 ($A$) or TNG300-1 ($B$) galaxy samples. The range of measurement is $0.03<k<0.9\;{\rm Mpc^{-1}h}$. The bar over the quantity means the value averaging over the all the fiducial samples of 4 redshift snapshots, 8 $\mathsf{FluxLimitLabel}$, 15$\,|\,$6 galaxy color and 6 bins for a color bin.  |  | $\bar{{\bf Q}}\times 100$ | worst ${\bf Q}\times 100$ | $\bar{\delta}_{A}$ | $\bar{\sigma}_{A}$ | worst $\delta_{A}$ |
---|---|---|---|---|---|---|---
| Expan-mm-1 | $1.92\,|\,1.41$ | $4.65\,|\,3.92$ | $0.095\,|\,0.044$ | $0.040\,|\,0.031$ | $0.327\,|\,0.112$ |
| Expan-mm-2 | $0.83\,|\,0.71$ | $1.45\,|\,1.27$ | $0.019\,|\,0.013$ | $0.043\,|\,0.037$ | $0.133\,|\,0.048$ |
| Expan-mm-3 | $0.77\,|\,0.65$ | $1.05\,|\,0.95$ | $0.011\,|\,0.010$ | $0.049\,|\,0.043$ | $0.071\,|\,0.029$ |
| Expan-gg-1 | $1.07\,|\,0.86$ | $2.62\,|\,1.63$ | $0.048\,|\,0.016$ | $0.035\,|\,0.028$ | $0.301\,|\,0.042$ |
| Expan-gg-2 | $0.87\,|\,0.71$ | $1.05\,|\,1.03$ | $0.013\,|\,0.012$ | $0.040\,|\,0.035$ | $0.094\,|\,0.026$ |
| Expan-gg-3 | $0.87\,|\,0.69$ | $1.06\,|\,1.04$ | $0.010\,|\,0.010$ | $0.047\,|\,0.041$ | $0.050\,|\,0.021$ |
| Expan-new-2 | $0.92\,|\,0.72$ | $1.65\,|\,1.20$ | $0.025\,|\,0.011$ | $0.043\,|\,0.037$ | $0.103\,|\,0.031$ |
| Zheng2013-mm | $0.72\,|\,0.70$ | $1.59\,|\,2.08$ | $0.015\,|\,0.009$ | $0.042\,|\,0.037$ | $0.116\,|\,0.033$ |
| Zheng2013-gg | $0.79\,|\,0.72$ | $1.65\,|\,1.80$ | $0.012\,|\,0.007$ | $0.042\,|\,0.038$ | $0.106\,|\,0.016$ |
| Expan-mm-1 ($r_{s}^{2}$) | $1.03\,|\,0.85$ | $3.38\,|\,2.00$ | $0.020\,|\,0.021$ | $0.030\,|\,0.027$ | $0.039\,|\,0.044$ |
| Expan-mm-2 ($r_{s}^{2}$) | $0.65\,|\,0.51$ | $3.14\,|\,0.87$ | $0.011\,|\,0.014$ | $0.038\,|\,0.034$ | $0.023\,|\,0.030$ |
| Expan-new-2 ($r_{s}^{2}$) | $0.73\,|\,0.57$ | $3.11\,|\,0.98$ | $0.010\,|\,0.012$ | $0.040\,|\,0.036$ | $0.019\,|\,0.026$ |
| Zheng2013-mm ($r_{s}^{2}$) | $1.07\,|\,0.93$ | $3.51\,|\,1.95$ | $0.008\,|\,0.011$ | $0.039\,|\,0.035$ | $0.019\,|\,0.030$ |
| Q-bias | $0.55\,|\,0.26$ | $3.15\,|\,0.86$ | $0.018\,|\,0.012$ | $0.050\,|\,0.045$ | $0.051\,|\,0.034$ |
| Quadratic | $0.35\,|\,0.19$ | $3.15\,|\,0.42$ | $0.006\,|\,0.006$ | $0.037\,|\,0.037$ | $0.024\,|\,0.018$ |
| Quadratic-1 | $5.41\,|\,3.62$ | $14.86\,|\,12.37$ | $0.088\,|\,0.070$ | $0.016\,|\,0.015$ | $0.185\,|\,0.166$ |
Figure 3:  The matter clustering amplitude $A$ measured in cosmoDC2 galaxies.
We show the bias $\delta_{A}=A-1$ and corresponding uncertainty $\sigma_{A}$
as function of maximum fit scale $k_{\rm max}$. Four rows show the results at
redshift $z=0.15,0.5,1.0$ and $1.5$ respectively. The first column shows the
results setting $r_{s}=1$, and the rest of 7 columns show the results for
$r^{2}_{s}$ parameterizations. In each subfigure, different colors denote the
full flux limited samples with different $\mathsf{FluxLimitLabel}$, and the
dash/dot lines mark the region of $\pm 2\%\,/\,\pm 5\%$.  Figure 4:  Same plot
as Fig. 3 for the full flux limited samples with different
$\mathsf{FluxLimitLabel}$ in cosmoDC2 galaxies. Here we present results for
$r^{2}$ parameterizations.  Figure 5:  Scatter plot of relation between
$\bar{\delta}_{A}$ and $\bar{\sigma}_{A}$ of fourth and fifth columns in Table
4. Black/Red markers denote the results of cosmoDC2/TNG300-1 galaxies. The
green solid line marks $\bar{\sigma}_{A}=\bar{\delta}_{A}$, and the dash lines
mark $\bar{\sigma}_{A}=\bar{\delta}_{A}\pm 1\%$ and
$\bar{\sigma}_{A}=\bar{\delta}_{A}\pm 2\%$.

##  IV  Performance of Proposed Parameterizations

In this section, we present the results of testing the fitting accuracy and
reconstruction ability of the parameterizations. In Section IV.1, we present
the direct $r^{2}$ fitting accuracy for our fiducial galaxy samples. In
Section IV.2, we further analyze the impact on the reconstructed amplitude
$A$, based on the fiducial galaxy samples and assuming the cosmic variance
limit for lensing measurement.

The fiducial galaxy samples for testing stochasticity consist of two sets of
galaxy mocks, cosmoDC2 galaxy mock [50] and TNG300-1 galaxy samples [51, 52,
53, 54, 55], combined with $1200\,h^{-1}{\rm Mpc}$ dark matter simulation. The
details of fiducial galaxy samples are listed in Table 3 $\&$ 3. The galaxy
samples are defined by the redshift, flux and the galaxy color. We apply
various magnitude cut of $ugrizY$/$griz$ bands to the cosmoDC2/TNG300-1
galaxies. These flux limited samples are labeled by tag
$\mathsf{FluxLimitLabel}$. Notice that $\mathsf{FluxLimitLabel}=3$ corresponds
to the expected depth of coadded images from LSST [4]. In order to test our
parameterization with various stochasticity, we further separate the full flux
limited galaxy sample into multi-bins according to the galaxy color. For
cosmoDC2 galaxy, we utilize all the possible color by $ugrizY$ bands
combination ($u-g$, $g-r$ and so on), obtaining totally $15$ kinds of galaxy
colors, and for TNG300-1 there is $6$ kinds of galaxy colors. The number of
galaxy color bins is $N_{\rm bins}=6$, and the number of galaxies in each bin
is identical. We also get the consistent results with $N_{\rm bins}=12$, but
we do not show details here to avoid redundancy.

###  IV.1 Goodness of Fit

In Fig. 1, we show the measured $\hat{r}^{2}$ on cosmoDC2 galaxy samples and
their fitting results $r^{2}$ by schemes Expan-mm-2 and Zheng2013-mm. As the
fit range $k_{\rm max}\lesssim 0.5\;{\rm Mpc^{-1}h}$, the accuracy generally
reaches $\lesssim 1\%$. When $k_{\rm max}$ extends to $0.9\;{\rm Mpc}^{-1}h$,
it still maintains the accuracy level at $3\%$ conservatively, and for most of
the samples shown in figure, we still have accuracy $\lesssim 1\%$. In Fig. 2,
we show the same results for $r^{2}_{s}$ by Quadratic and Expan-mm-2. We can
see a tighter fitting for $r_{s}^{2}$, especially for the full sample, the
Quadratic accuracy is strictly at the level of $<1\%$ at $k_{\rm
max}<0.9{\,\rm Mpc^{-1}}h$ for samples at redshift $z>0.5$. A visual advantage
of $r_{s}^{2}$ fit results is that, our parameterizations are able to
characterize the upwarp around $k\sim 1{\,\rm Mpc^{-1}}h$ due to the halo
exclusion effect, for instance Quadratic shown in Fig. 2. The terms such as
$k^{2}$ are unfavorable in direct $P_{\mathcal{S}}$ fitting, due to the
divergence in Fourier transform, but fortunately this property is able to
capture the sub-Poisson behavior appearing in the non-linear scale, enhancing
the flexibility of $r_{s}^{2}$ fitting in turn. Similar results are obtained
for other parameterization schemes.

We utilize the mean ${\bf\bar{Q}}$ values to summarize the performance of the
proposed parameterizations against the above galaxy samples. Here,
${\bf\bar{Q}}$ means averaging ${\bf Q}$ values over the all the galaxy
samples of 4 redshift snapshots, 8 $\mathsf{FluxLimitLabel}$, 15 (6) galaxy
color and 6 bins for a color bin. The range of k is fixed as $0.03\;{\rm
Mpc^{-1}}h<k<0.9\;{\rm Mpc^{-1}}h$, where the larger scale modes $k\lesssim
0.03{\,\rm Mpc^{-1}}h$ are affected greatly by cosmic variance and the
definition of $\bf Q$ does not account for the fluctuation in measurement. We
list the results in first and second columns in Table 4. Based on
${\bf\bar{Q}}$ comparison, we can derive two immediate conclusions:

  * $\bf(1)$

Two free parameters in parameterization are able to describe stochasticity to
small scale $k_{\rm max}=0.9\;{\rm Mpc^{-1}}h$ in $\sim 1\%$ accuracy,
verifying the implication by the previously found ${\rm DOF}\simeq 2$
eigenmodes in simulated galaxy clustering.

  * $\bf(2)$

The Quadratic scheme is the most promising parameterization with the best
$<0.5\%$ accuracy in $r_{s}^{2}$ fitting among all 12 proposed
parameterization schemes, quantified by the summarized performance over all
the samples.

Specifically, the inclusion of more free parameters do benefit the
$r^{2}_{(s)}$ fitting. For instance, when we replace the scheme Expan-mm-1 by
Expan-mm-2, the ${\bf\bar{Q}}$ value decreases $1.92\%\rightarrow 0.83\%$ for
cosmoDC2 galaxies and $1.41\%\rightarrow 0.71\%$ for TNG300-1 galaxies, almost
twice as much. However, the ${\bf\bar{Q}}$ decrement when replacing Expan-mm-2
by Expan-mm-3 is not as much as for Expan-mm-1 $\rightarrow$ Expan-mm-2,
producing little improvement especially for Expan-gg-2 $\rightarrow$ Expan-
gg-3. For all $2$-parameter $r^{2}$ parameterizations, we have
${\bf\bar{Q}}\sim 0.8\%$, and the worst values appears among all the galaxy
samples is ${\bf Q}<2\%$. For all $2$-parameter $r^{2}_{s}$ parameterizations,
we have similar value ${\bf\bar{Q}}\lesssim 1\%$ for all samples, but it seems
that worst ${\bf Q}\sim 3.5\%$ of cosmoDC2 samples is larger than that for
$r^{2}$. It arises from the intrinsic true that $r_{s}>r$, and the deviation
$\delta r$ is estimated by ${\bf Q}(r)\sim\delta r$. Therefore we expect that
${\bf Q}(r_{s})>{\bf Q}(r)$ in general, and the larger ${\bf Q}(r_{s})$ value
compared to ${\bf Q}(r)$ do not confirm $r^{2}$ fits outperforms $r^{2}_{s}$
fits. Even so, some $r^{2}_{s}$ parameterizations still achieve a smaller
${\bf\bar{Q}}$ value compared to $r^{2}$ parameterizations, especially
Quadratic, of which the ${\bf\bar{Q}}$ is smaller than any other
parameterization even those with $3$ parameters. It suggests that Quadratic is
a promising candidate to characterize the stochasticity.

Besides the fiducial color-bins samples, we further separate the full flux
limited galaxy samples into multi-bins by stellar mass and host-halo mass. In
the stellar mass classified galaxy samples, the parameterization performance
is consistent with the fiducial color-bins case. The absolute values of the
measured $\bf\bar{Q}$ shows little difference compared to the color-bins
samples, and the Quadratic scheme is still the best performing
parameterization with ${\bf\bar{Q}}=0.39\%\,\&\,0.17\%$ for
cosmoDC2$\,\&\,$TNG300-1 samples, closed to the fiducial case. However, the
situation changes if we adopt host-halo mass as the classification criteria,
since it separates the galaxies into satellite-rich and satellite-poor galaxy
populations, as well as into the halo exclusion-effective and exclusion-
ineffective galaxy populations. As it has derived in [19, 31], the large scale
modes can not identify the galaxies in a same host-halo, resulting in the
host-halo dominating shot noise in large scale, sourcing the sup-Poisson
enhancement. So the satellite-poor bin is affected less by the enhancement,
while the satellite-rich bin is enhanced greatly. Apart from the enhancement,
the halo exclusion effect is more and more significant as the host-halo mass
increase, and it sources the sub-Poisson suppression. Within our performance
tests on cosmoDC2$\,\&\,$TNG300-1, we find the sup-Poisson enhancement is
dominant in the first 5 less-massive bins, and the there is a trend of
increasing $r_{s}^{2}$ suppression as host-halo mass increases. For these 5
bins, a particular feature is the significant suppression $r_{s}^{2}\lesssim
0.95$ even just at $k=0.05{\,\rm Mpc^{-1}}h$, distinguished from the color-
bins as well as the stellar mass-bins samples. While for the most massive bin,
the halo exclusion overcomes the sup-Poisson effect and alleviates the
suppression on $r_{s}^{2}$, allowing $r_{s}^{2}\sim 1$ for cosmoDC2 galaxies.
Despite significant distinction in these host-halo mass bin samples, our 2-
$\&$ 3-parameter parameterizations still work well within range
$0.03<k<0.9{\,\rm Mpc^{-1}}h$. In contrast, the 1-parameter fits such as
Expan-mm-1 reach ${\bf\bar{Q}}=3.7\%\,\&\,3.6\%$, revealing that one degree of
freedom is not flexible enough to parameterize the stochasticity. For the
2-parameter fitting, ${\bf\bar{Q}}$ deviates slightly from the results in
Table 4, and the worst $\bf Q$ appearing in $r^{2}$ parameterizations is
limited to ${\bf Q}<3\%$. Especially, Quadratic still provides the best
performance with ${\bf\bar{Q}}=0.53\%\,\&\,0.61\%$. To sum up, all the tests
on various samples reinforce two conclusions above, reiterating the suggestion
that the simple Quadratic is a promising choice for stochasticity
parameterization

###  IV.2 Impact on Matter Clustering Amplitude

We utilize $\delta_{A}$ and $\sigma_{A}$ to quantify the ability of
parameterizations alleviating the stochasticity impact. Solely for the purpose
of comparison, we assume the distribution of the lensed galaxy sources is an
infinite thin plane, locating at $z_{s}=3$, and the redshift bin width of the
foreground clustering galaxies is $\Delta z=0.1$, with $f_{\rm sky}=1$ for
both sets. For simplicity, we assume cosmic variance limit for the cosmic
shear measurements, equivalently neglecting shotnoise in $\kappa$ field. The
neglected complexities such as the redshift distribution of source galaxies
and shape measurement noise will be taken into account later in forecasting
the DESI-CSST performance. The galaxy clustering shotnoise measured in
simulation is kept in variance, since it is an important ingredient of
stochasticity.

To make an intuitive comparison, we demonstrate $(\delta_{A},\sigma_{A})$
measurement for the full flux limited samples in cosmoDC2 galaxies in Fig. 3
$\&$ 4. We also present the results neglecting the impact of stochasticity,
equivalently setting $r_{s}=1$, as the baseline of the systematic bias and the
statistic uncertainty. In the context of amplitude $A$ quantification, the
conclusions are same as that based on the ${\bf\bar{Q}}$. Some details are
listed below.

  * •

The statistic uncertainty decrease as redshift increase, as it is shown in the
first row of Fig. 3, resulted from the larger cosmic volume therefore higher
surface number density.

  * •

At redshift $z=1.0$ and $z=1.5$, the systematic bias of luminous populations
is softened, counteracted by the enhanced sub-Poisson feature for the massive
host halo population.

  * •

When an appropriate stochasticity parameterization is included in modeling,
$\delta_{A}$ is alleviated while $\sigma_{A}$ is enlarged in different levels.
The 2-parameter $r^{2}_{s}$ fits present a maximum benefit on constraining
matter amplitude $A$, while the 1-parameter and 3-parameter fits do not
achieve a competitive gain.

  * •

For 1-parameter fits, most of the galaxy samples retain bias at the level
$\delta_{A}\gtrsim 2\%$ at $k_{\rm max}=0.9{\,\rm Mpc^{-1}}h$. The Quadratic-1
appears as the worst case, only few percent systematic bias is corrected, and
there is still $\delta_{A}\sim 10\%$ at $z=0.5$. While for parameterizations
such as Expan-gg-1 and Expan-mm-1, the systematic bias is controlled at the
level $\delta_{A}\lesssim 5\%$ at $k_{\rm max}\simeq\,1{\,\rm Mpc^{-1}}h$. The
distinction is expected in the Eq. (8), which reveal that the stochasticity
power spectrum is correlated with large scale clustering in higher order
correlation terms, and the parameterizations such as Expan-gg-1 and Expan-mm-1
make use of the galaxy/matter clustering information while a naive low $k$
expansion Quadratic-1 makes flat assumption of stochasticity.

  * •

For 2-parameter fits, the reconstruction accuracy is almost
$\delta_{A}\lesssim 2\%$ at $k_{\rm max}=0.9{\,\rm Mpc^{-1}}h$ for $r^{2}$,
and is generally $\delta_{A}<2\%$ for $r^{2}_{s}$. However, for the luminous
galaxy populations at $z=1.5$, a deviation $2\%\lesssim\delta_{A}\lesssim 5\%$
appears in $r^{2}$ parameterizations including $P_{mm}$ (i.e. Expan-mm-2) but
not in those including $P_{gg,s}$ (i.e. Expan-gg-2). We suppose the poor
performance is caused by the strong halo exclusion in these massive host-halo
populations, which contributes an enhanced scale dependent bias in non-linear
region. The prominent non-linear bias is naturally included in the $P_{gg,s}$,
but not related with $P_{mm}$. Fortunately, this issue does not affect the
$r_{s}^{2}$ fits, because the overall halo exclusion always emerges as an
upwarp in $\hat{r}_{s}^{2}$, which is well characterized by our $r^{2}_{s}$
parameterizations.

  * •

For 3-parameter fits, $\delta_{A}$ is not corrected significantly compared to
2-parameter fits while $\sigma_{A}$ is magnified. For instance, when we
replace parameterization Expan-mm-2 by Expan-mm-3 for the worst performing
luminous galaxy populations, the correction of $\delta_{A}$ is improved by
several percent, but $\sigma_{A}$ is amplified by about $20\%$.

The above discussion is based on cosmoDC2 samples, and the similar results are
also obtained in the TNG300-1 galaxy samples.

In Table 4, we also show the mean value $\bar{\delta}_{A}$ and
$\bar{\sigma}_{A}$ for the fiducial samples. Here the average value is defined
by RMS, $\bar{\delta}_{A}=\sqrt{\sum_{i}{\left(\delta_{A,i}\right)^{2}}}$ and
$\bar{\sigma}_{A}=\sqrt{\sum_{i}{\sigma_{A,i}^{2}}}$, where the summation
index $i$ runs over all the galaxy samples. Now $\delta_{A}$ serves as a fair
criterion to quantify the parameterization ability reconstructing the matter
amplitude $A$ among different parameterizations. These summarized indexes
reinforce above conclusions that the simple Quadratic outperforms other
parameterizations. Specifically, the recommended 2-parameter fits reach about
$1\%\sim 2\%$ accuracy of $\delta_{A}$, and approximate $3.5\%\sim 4\%$
precision of $\sigma_{A}$ for two sets of galaxy samples. Remarkably,
Quadratic scheme outperforms other schemes with a minimal $\delta_{A}=0.006$.
What’s more, we present the $\delta_{A}$-$\sigma_{A}$ relation scatter plot in
Fig. 5 for a visual comparison. The scatters distribute around the diagonal,
exhibiting a negative correlation, and it reflects the balance of benefit in
systematic bias correction, better reducing $\delta_{A}$ while worse
amplifying $\sigma_{A}$. We exclude the 1-parameter $r^{2}$ fits on the right
side in Fig. 5 because of the large systematic bias, and exclude the
3-parameter $r^{2}$ fits and Q-bias scheme on the top side because of the
large uncertainty. The rest 2-parameter $r^{2}\,\&\,r^{2}_{s}$ schemes and
$\texttt{Expan-mm-1}\;(r^{2}_{s})$ scheme are applicable, in particular the
Quadratic scheme.

Apart from the intrinsic distinction in galaxy samples such as redshift, flux
limit and mass, though not presented specifically, the simulation details also
play an important role in stochasticity. Two galaxy catalogs utilized here
exhibit obviously different behaviors. For instance, there is relatively
subtle stochasticity for TNG300-1 samples compared to cosmoDC2, though with
completely same selection criteria. Conclusions derived from two individual
sets are consistent, but the performance details are different.

Table 5:  Survey Parameters222 Totally 7 redshift bins for cross correlation analysis. Here, $\Sigma_{g}$ /$n_{g}$ are the surface/volume number density of DESI LRG in each spec-z bin [2], and $\Sigma_{\gamma}$ is the surface number density of CSST cosmic shear galaxy in each phot-z bin.  DESI LRG | CSST cosmic shear
---|---
$z$ | $\Sigma_{g}\,[{\rm deg^{-2}}]$ | $n_{g}\,[{\rm Mpc^{-3}h^{3}}]$ | $z_{p}$ | $\Sigma_{\gamma}\,[{\rm deg^{-2}}]$
$[0.4,\,0.5]$ | $47.5$ | $4.38\times 10^{-4}$ | $[0.7,\,3]$ | 68208
$[0.5,\,0.6]$ | $65.6$ | $4.49\times 10^{-4}$ | $[0.8,\,3]$ | 60462
$[0.6,\,0.7]$ | $80.0$ | $4.37\times 10^{-4}$ | $[0.9,\,3]$ | 52983
$[0.7,\,0.8]$ | $93.2$ | $4.25\times 10^{-4}$ | $[1.0,\,3]$ | 45959
$[0.8,\,0.9]$ | $99.3$ | $3.93\times 10^{-4}$ | $[1.2,\,3]$ | 33693
$[0.9,\,1.0]$ | $63.7$ | $2.24\times 10^{-4}$ | $[1.3,\,3]$ | 28528
$[1.0,\,1.1]$ | $28.3$ | $0.91\times 10^{-4}$ | $[1.3,\,3]$ | 28528

Figure 6:  Forecasted stochasticity impact on the matter amplitude uncertainty
$\sigma_{A}$ for DESI-like galaxy $\times$ CSST-like cosmic shear. In both top
and bottom panels, 7 columns show the results at different redshift bins. The
top subfigures show systematic bias $\delta_{A}$ by assuming $r_{s}=1$ (black
line) and the fiducial case $A=1$ with $r_{s}^{2}$ parameterized by Quadratic
(cyan line), and the $1\sigma$ regions are marked. The bottom subfigures show
statistic uncertainty $\sigma_{A}$ under various parameterizations compared to
the case $r_{s}=1$ (black line). The cosmological parameters are fixed as dark
matter simulation described in Appendix A.  Figure 7:  Forecasted statistic
uncertainty $\sigma_{A}$ as function of redshift bin $z$ for DESI-like galaxy
$\times$ CSST-like cosmic shear. The result shown here is fixed with
$\ell_{\rm max}\chi_{g}=0.9\,h^{-1}{\rm Mpc}$.  Figure 8:  Forecasted
statistic uncertainty of cosmological parameter $S_{8}$ as function of
stochasticity parameters $(c_{1},c_{2})$. The blank region in left bottom is
null value because of unphysical $r^{2}<0$ for these parameters.

##  V Forecast for DESI-like LRG

The stochasticity parameterization alleviates the systematic bias to a large
extent, but including nuisance parameters reduces the power of cosmological
parameters constraint. Here we forecast the specific impact on the combination
of DESI-like galaxy and CSST-like cosmic shear. We adopt Quadratic scheme as
the fiducial stochasticity behavior for DESI-like galaxies, because it has
been shown to present the best performance among all candidates.

The forecasted surface number density of DESI-like galaxies and the estimation
of CSST-like cosmic shear survey are summarized in Table. 5. For the DESI
spectroscopic tracer, we only utilize LRG in our analysis. Though the previous
galaxy samples selection differs from DESI LRG selecting strategy [56], it is
supposed to affect little on our forecast. On the one hand, DESI LRG sample is
selected using both magnitude and color, approximately being flux limited
samples for which we have performed various tests. On the other hand, our
parameterizations have been shown to be powerful enough to characterize
various kinds of stochasticity, so it is expected to be applicable to
realistic DESI LRG. The validation and the forecast for the DESI ELG samples,
the major targets with complex identification [57], is remained to future
investigation. In this work, LRG samples are divided into $7$ redshift bins
covering redshift from $0.4$ to $1.1$ with $\Delta z=0.1$. We adopt the galaxy
bias $b(k)\equiv b_{L}=1.7/D(z)$, where $D(z)$ is the linear growth factor
[2].

For the cosmic shear tracers, CSST photometric galaxy survey is expected to
reach a total galaxy surface number density $\Sigma_{g}=28\,{\rm arcmin^{-2}}$
with redshift scatter $\sigma_{z}=0.05$. In our analysis, the redshift
distribution is adopted as $n(z)\propto z^{2}e^{-z/z_{*}}$ covering redshift
$0<z<4$, where $z_{*}=0.35$. The photometric redshift bins are selected to
avoid the overlapping of spectroscopic galaxy and the cosmic shear tracers, as
well as minimize the shotnoise of $\kappa$ field. We assume the Gaussian
variance for CSST-like cosmic shear, and adopt the estimated shape noise
$\sigma_{\gamma}=0.2$ and additive noise $N_{\rm add}^{\gamma}=10^{-9}$ [7].
The overlapping fraction of sky with DESI-like galaxies is $f_{\rm sky}=0.3$.
To avoid complexity, we have neglected the impacts of magnification bias,
photometric error and other systematic uncertainties.

For illustrative purpose, we simply make an aggressive guess of stochasticity
with Quadratic parameters $c_{1}=-1.57$, $c_{2}=0.74$ for all 7 galaxy bins
(called Scenario-A below), and treat all $2\times 7$ nuisance stochasticity
parameters as free parameters when applying fisher forecast. We do not show
the results for a scheme such as Expan-gg-2, because it is equivalent with
Expan-mm-2 under the linear bias assumption. We show the results under
Scenario-A stochasticity in Fig. 6. If we neglect the impact of stochasticity
and adopt $r_{s}=1$, the systematic bias $\delta_{A}$ would exceed the
statistic uncertainty at $k_{\rm max}\equiv\ell_{\rm max}\,\chi_{g}\simeq
0.1{\,\rm Mpc^{-1}}h$, where $\delta_{A}\simeq\sigma_{A}\simeq 0.05$. As
$k_{\rm max}$ extends to $0.9{\,\rm Mpc^{-1}}h$, as shown in Fig. 7, the
uncertainty of $r_{s}=1$ is reduced to $\sigma_{A}\simeq 0.01$ while the
underestimation is serious $\delta_{A}=0.3\sim 0.4$. When we adopt
parameterization Quadratic, the matter clustering amplitude is estimated
unbiasedly with precision $\sigma_{A}\simeq 0.04$ for the first 5 bins and
$\sigma_{A}\simeq 0.045,\;0.059$ for the latter 2 bins. In the bottom
subfigures of each column in Fig. 6, we further compare the abilities of other
parameterization schemes. Among these various schemes, the values of
$\sigma_{A}$ bifurcate into two branches as $\ell_{\rm max}$ increase. One
branch has larger $\sigma_{A}$ for $r^{2}$ parameterizations, while another
branch has smaller $\sigma_{A}$ for $r^{2}_{s}$ parameterizations. The reason
for these distinctive behaviors of two kinds of parameterization schemes is
the flattened effect of the discrete noise in galaxy power spectrum. If we do
not subtract the Poisson noise expectation when $\bar{n}_{g}P_{gg}\ll 1$, the
curve of $r^{2}$ is suppressed to a nearly flat shape, that $r_{s}>r$. Thus
the modeling $\hat{A}=rA$ is less sensitive compared to the modeling
$\hat{A}_{s}=r_{s}A$, quantitatively,
$|\frac{\partial\hat{A}_{s}}{\partial\lambda}|>|\frac{\partial\hat{A}}{\partial\lambda}|$,
where $\lambda$ is the cosmological parameter we interest in. It leads to a
larger uncertainty of $A$ in $r^{2}$ modeling. Consistent with the argument,
the bifurcation is significant for the last redshift bin, of which the number
density is lowest, reaching $\bar{n}_{g}P_{gg,s}=0.1$ at $\ell_{\rm max}\simeq
1500$. For the relative higher number density bin, where
$\bar{n}_{g}P_{gg,s}=0.1$ at $\ell_{\rm max}\simeq 3200\sim 5000$, the
bifurcation is mild.

The underestimation in matter clustering amplitude would propagate to
cosmological constraint and result into a substantial effect, sourcing a
potential systematic bias in cosmological parameters. Different experiments
have reported the tendency that the lower redshift measurements prefer lower
matter clustering amplitude compared with the CMB constraint in the context of
$\Lambda{\rm CDM}$ [10, 11, 58, 12]. For the cosmic shear measurement, this
anomaly is usually summarized as the structure growth parameter
$S_{8}\equiv\sigma_{8}(\Omega_{m}/0.3)^{0.5}$. Since the value of $S_{8}$ (or
$\sigma_{8}$) is degenerated with $r$, an appropriate $r^{2}$ modeling
excludes the potential systematic bias thus obtain an unbiased constraint. We
demonstrate the stochasticity impact on cosmological constraint by varying the
parameters $\Omega_{m}$ and $\sigma_{8}$, and present the constraint
performance for $S_{8}$ together with/without stochasticity parameterization.
We find:

  * •

In the Scenario-A, an aggressive guess of serious stochasticity, our fiducial
results with $k_{\rm max}=0.9{\,\rm Mpc^{-1}}h$ are $\sigma_{8}=0.83\pm 0.03$
and $\Omega_{m}=0.268\pm 0.012$, with precision $3.5\%$ and $4.5\%$
respectively. The structure growth is $S_{8}=0.784\pm 0.012$, with precision
$1.5\%$.

  * •

In the Scenario-A, the neglect of stochasticity imports significant bias on
cosmological parameter. Setting $r_{s}=1$ and even at the linear scale $k_{\rm
max}=0.1{\,\rm Mpc^{-1}}h$, the inferred $S_{8}=0.746\pm 0.012$ deviates from
the fiducial values $S_{8}$ on the level of $3\sigma$. The deviation is
moderate if we directly constrain $\sigma_{8}$, which is only about
$0.3\sigma$ deviated from fiducial value at $0.1{\,\rm Mpc^{-1}}h\leq k_{\rm
max}\leq 0.3{\,\rm Mpc^{-1}}h$. But the $\sigma_{8}$ precision is limited to
$>3\%$ at $k_{\rm max}<0.3{\,\rm Mpc^{-1}}h$ due to the degeneration of
$\sigma_{8}$ and $\Omega_{m}$, and its systematic bias reaches about $3\sigma$
again at $k_{\rm max}\gtrsim 0.4{\,\rm Mpc^{-1}}h$.

  * •

In the Scenario-B, a mild stochasticity hypothesis with
$c_{1}=-0.63,c_{2}=0.27$, the fiducial results are $\sigma_{8}=0.83\pm 0.03$,
$\Omega_{m}=0.268\pm 0.012$ and $S_{8}=0.784\pm 0.013$, with constraint
precision $3.4\%,\,4.5\%$ and $1.7\%$ respectively. The precision varies
slightly compared with the Scenario-A, though we adopt quite different
$r_{s}^{2}$ value.

  * •

In the Scenario-B, the deviation is mitigated to a large extent compared to
the Scenario-A, but when we set $r_{s}=1$, the true value of $S_{8}$ is still
outside the $2\sigma$ region of inferred values at $k_{\rm max}=0.3{\,\rm
Mpc^{-1}}h$. As $k_{\rm max}$ increases, the deviation is more and more
serious, obtaining $S_{8}=0.736\pm 0.007$ at $k_{\rm max}=0.5{\,\rm
Mpc^{-1}}h$, a $7\sigma$ deviation. While for the $\sigma_{8}$ constraint, the
systematic bias is not significant compared with statistic fluctuation, but we
still have $\sigma_{8}=0.79\pm 0.02$ at $k_{\rm max}=0.5{\,\rm Mpc^{-1}}h$, a
$2\sigma$ deviation.

  * •

In Fig. 8, we present the statistic uncertainty for the fiducial forecasted
$S_{8}$ by varying the stochasticity parameters. In the parameter space, the
variation of $(c_{1},c_{2})$ introduces distinguished behavior of
stochasticity, as well as suppression of matter clustering amplitude in
different levels. But the statistic uncertainty for the final constrained
$S_{8}$ only varies in a narrow range $0.012\lesssim\sigma_{S_{8}}\lesssim
0.014$. We suppose the uncertainty is limited by the large number of
constraint parameters (total 14 nuisance parameters plus 2 cosmological
parameters). Thus there are tight bounds for $\sigma_{S_{8}}$ with various
$(c_{1},c_{2})$, rather than a strong dependence on the specific shape of
$r_{s}^{2}$. Similar tight bounds arise for both $\Omega_{m}$ and $\sigma_{8}$
uncertainty. Consequently, the uncertainty of the inferred parameters is
controlled in a limited range, not matter what kinds of stochasticity.

We aim to provide an unbiased cosmological parameters prediction, but one may
wonder that an accurate stochasticity model seems to pay a high price in
constraining power. We suppose this is not an issue for the following reasons.
Firstly, the vast parameter set (14 nuisance parameters for 7 redshift bins)
is not necessary in realistic case, since $\Delta z$ is so small that we
expect the clustering as well as stochasticity slowly evolves relative to
adjacent redshift bins. Thus some simple continuity hypothesis is enabled and
the number of nuisance parameters is reduced. Secondly, with the assistance of
calibration in the targeted mock, some proper priors of stochasticity
parameters are available [31], then reducing the uncertainty. Finally, the
stochasticity components are considerable compared to clustering signal in
current spectroscopic survey, such as Ref. [48, 59]. And stochasticity
modeling is also important to obtain a consistent result in the future multi-
probes applicants [28, 30]. Thus the inclusion of stochasticity nuisance
parameters is inevitable, and the key, also our future work, is how to further
reduce the degree of freedom in stochasticity parameters for tomographic
analysis.

##  VI Conclusion and Discussions

In this work, we addressed two leading questions: how to parameterize
stochasticity in galaxy clustering, and how to reconstruct the unbiased
tomographic matter clustering. In order to seek an appropriate
phenomenological description of stochasticity, we investigate 12 kinds of
parameterization schemes for the cross correlation coefficients between galaxy
overdensity $\delta_{g}({\bf k})$ and the underlying matter overdensity
$\delta_{m}({\bf k})$. We test the performance against the galaxy samples
selected from cosmoDC2 and TNG300-1 simulations, over a wide range of
redshift, flux and various galaxy colors. We quantify the fitting performance
by the absolute difference of $r_{(s)}(k)$ values and the accuracy $\&$
precision of reconstructed matter clustering amplitude. Then we choose the
best performing quadratic scheme as the fiducial stochasticity, and present
the fisher forecast of the matter clustering amplitude reconstruction as well
as the cosmological parameters constraint, assuming the combination DESI-like
LRG and CSST-like cosmic shear. Our main conclusions are as following:

  1. (i)

We verify the suggestion the effective ${\rm DoF}\simeq 2$ in galaxy
stochasticity, and that 2-parameter scheme provides the maximum gain in the
stochasticity description compared to 1-parameter or 3-parameter scheme.
Against two sets of galaxy samples, we have shown that 2-parameter fits are
able to reach about $1\%$ accuracy for $r_{(s)}^{2}$, and obtain general
$0.6\%\sim 2\%$ accuracy for the reconstructed matter clustering amplitude at
$k_{\rm max}=0.9{\,\rm Mpc^{-1}}h$.

  2. (ii)

The quadratic scheme $r^{2}_{s}=1+c_{1}k+c_{2}k^{2}$ is the most promising
parameterization among all the schemes, performing the best direct fit for
$\hat{r}_{s}^{2}$ and reconstructing the matter clustering amplitude with
better than $1\%$ accuracy.

  3. (iii)

For the combination DESI-like LRG $\times$ CSST-like cosmic shear, we forecast
that there is a suppression in the reconstructed tomographic matter clustering
amplitude when neglecting the stochasticity. The systematic bias sourced by
the suppression is comparable with statistic uncertainty at linear scale
$k_{\rm max}=0.1{\,\rm Mpc^{-1}}h$ if the stochasticity contamination is
serious. The serious stochasticity also biases the inferred $S_{8}$ at
$3\sigma$ level even at $k_{\rm max}=0.1{\,\rm Mpc^{-1}}h$.

  4. (iv)

Together with the quadratic stochasticity parameterization, the combination
DESI-like LRG $\times$ CSST cosmic shear is expected to achieve $S_{8}$
constraint at $1.5\%$ precision, free from the stochasticity systematic. The
uncertainty is likely to be further reduced with stochasticity parameter
number reduced.

We caution that this work is still theoretical. Its applicability in real data
requires validation in targeted mocks, which must be done case to case.

The galaxy-matter cross power spectrum is also available by the cross
correlation with CMB lensing. The CMB lensing is much cleaner than the cosmic
shear estimated from the galaxy ellipticity, and the photometric redshift
distribution of lensed galaxies introduces large uncertainty while the CMB
redshift is accurately known. Thus we expect the combination with galaxy
clustering $\times$ CMB lensing achieves a tighter tomographic matter
clustering constraint [9]. Apart from the direct reconstruction of tomographic
matter clustering from galaxy/CMB lensing, there are wide applications for an
accurate $r^{2}$ parameterization. For example, it could alleviate the
systematic bias due to stochastic components in the estimator of $D_{G}$
statistic [60, 61, 62, 12]. In the condition of non-negligible stochasticity,
the theoretical expectation of estimator $\hat{D}_{G}$ (refer to Ref. [60] for
detailed expression) can be rewritten to correct the stochasticity
suppression,
$\langle\hat{D}_{G}(\ell)/r\left(\frac{\ell}{\chi}\right)\rangle=D(z)$, where
$D(z)$ is the linear growth function.

We plan to validate our methods with the simulations designed for DESI
spectroscopic redshift surveys, in particular the ELG mocks. After the
validation, we will apply our methods to the DESI spectroscopic sources to
reconstruct the matter power spectrum. Meanwhile, our analysis simultaneously
provides the stochasticity parameters for the galaxy samples, therefore we can
determine whether the stochasticity is a serious systematic bias in the joint
analysis of galaxy auto- and cross-power spectrum or not. It would serve as an
important supplementation for systematic diagnosis in the future DESI multi-
probes applications. We expect our reconstructions will be complementary to
DESI constraint of dark energy and modified gravity by BAO and RSD
measurement.

###### Acknowledgements.

We would like to thank Yu Yu’s help in the access and usage of simulations,
and thank Hong-Ming Zhu for kind suggestions. This work is supported by the
National Key R&D Program of China
(2023YFA1607800,2023YFA1607801,2020YFC2201602), the National Science
Foundation of China (11621303), CMS-CSST-2021-A02, and the Fundamental
Research Funds for the Central Universities. This work made use of the Gravity
Supercomputer at the Department of Astronomy, Shanghai Jiao Tong University.
We acknowledge the use of code package $\mathtt{nbodykit}$ [63], CAMB [64] and
CCL [65].

Figure 9:  Scatter plots of parameters distribution in 4 redshift snapshots
based on cosmoDC2 samples. The top panels show parameters $(\alpha,\beta)$ in
Expan-mm-2, and the bottom panels show parameters $(c_{1},c_{2})$ in
Quadratic, with $k_{\rm max}=0.9{\,\rm Mpc^{-1}}h$. In each subfigure, all
various galaxy color-bins with same $\mathsf{FluxLimitLabel}$ are marked with
same color.  Figure 10:  Comparison of parameters $(c_{1},c_{2})$ in Quadratic
between cosmoDC2 and TNG300-1 galaxy samples. In each subfigure, all different
$\mathsf{FluxLimitLabel}$ and different galaxy color-bins in same redshift are
marked with same color, with $k_{\rm max}=0.9{\,\rm Mpc^{-1}}h$.

##  Appendix A Galaxy Samples

The performance testing for the parameterizations in Table 1 is based on the
fiducial color-bins samples, where the specific definitions of various
galaxies refer to Table. 3 and Table. 3. The fiducial galaxy samples consist
of two sets of galaxy samples, cosmoDC2 mock [50] and TNG300-1 galaxy samples
from the IllustrisTNG project [51, 52, 53, 54, 55]. The cosmoDC2 is a
synthetic galaxy catalog generated for LSST survey [50], and it provides a
lightcone covering $440\,{\rm deg}^{2}$ of sky area, with $\sim 2.6$ billions
of galaxies distributed to redshift $z=3$. The magnitude $ugrizY$ photometirc
bands and the host halo are provided for the galaxies. The hydrodynamic
simulation TNG300-1 is a publicly available simulation from the IllustrisTNG
project [51, 52, 53, 54, 55]. IllustrisTNG is a suite of large volume,
cosmological, gravo-magnetohydrodynamical simulations run with the moving-mesh
code ${\tt APEPO}$ [66]. IllustrisTNG implements a set of physical processes
to model the galaxy formation process including gas cooling, star formation
and evolution, supernova feedback [67], and AGN feedback [68]. This project
includes TNG50, TNG100 and TNG300, with boxsize of $50$ Mpc, $100$ Mpc and
$300$ Mpc respectively. We choose TNG300-1 for our analysis for its largest
amount of galaxy samples. We identify all luminous subhalos in the TNG300-1 as
galaxies.

However, the accessed cosmoDC2 products do not included the associated 3D
matter field, which is based on the Outer Rim simulation [69]. Regarding the
TNG300-1 simulation, the products we assess are complete but the boxisize
$L\,=\,205\,{\rm Mpc/h}$ is limited for our analysis, which lacks the
statistic for linear scale $k\lesssim 0.1{\rm\,Mpc^{-1}}h$ and induces large
fluctuation in $r^{2}$. To solve these problems, we combine the galaxies in
cosmoDC2$\;\&\;$TNG300-1 with dark matter CosmicGrowth simulation [70]. The
CosmicGrowth simulation adopts a flat WMAP cosmology [71], of which the
boxsize is $L\,=\,1200\,{\rm Mpc/h}$ and particle number is $3072^{3}$. The
cosmological parameters adopt
$\Omega_{b}=0.0445\,,\;\Omega_{c}=0.2235\,,\;h=0.71$ and $\sigma_{8}=0.83$. We
use the galaxy-host halo relation of the galaxy in cosmoDC2$\;\&\;$TNG300-1 to
assign galaxies to the CosmicGrowth simulation halos, preserving the
conditional probability of galaxy contents given hosthalo mass [50]. Because
of the distinction of these simulations resulted from different cosmological
parameters and the simulation details, we do an abundance matching between
cosmoDC2$\;\&\;$TNG300-1 and CosmicGrowth halos, and then assign galaxies in
cosmoDC2$\;\&\;$TNG300-1 to CosmicGrowth halos.

##  Appendix B Gaussian Variance

The direct fitting of the measured $r^{2}$ in the simulations to our
parameterized formula is through minimizing the Eq. (12), where the covariance
is assumed to be Gaussian. Start with the covariance of power spectrum of
Gaussian fields $i,\,j,\,p,\,q$,

$\displaystyle{\rm\bf
Cov}\left[P_{ij}(k),P_{pq}(k^{\prime})\right]=\frac{\delta_{k,k^{\prime}}^{K}}{N_{\bf
k}}\left(P_{ip}P_{jq}+P_{iq}P_{jp}\right)\;,$ (24)

it is straightforward to derive the Gaussian variance of $r^{2}$ and
$r^{2}_{s}$,

$\displaystyle\sigma^{2}|_{r^{2}}$ $\displaystyle=$
$\displaystyle\frac{4}{N_{k}}\left(r^{2}\right)^{2}\left(\frac{1}{r^{2}}+r^{2}-2\right)\quad,$
(25) $\displaystyle\sigma^{2}|_{r^{2}_{s}}$ $\displaystyle=$
$\displaystyle\frac{4}{N_{k}}\left(r_{s}^{2}\right)^{2}\left[\frac{1}{r_{s}^{2}}+r_{s}^{2}-2\right.$
$\displaystyle+\left.\frac{1}{\bar{n}_{g}P_{gg,s}}\left(\frac{1}{r_{s}^{2}}-1\right)+\frac{1}{2\left(\bar{n}_{g}P_{gg,s}\right)^{2}}\right]\quad.$

##  Appendix C Parameter Distribution

In Fig. 9, we demonstrate the parameters distribution for two schemes Expan-
mm-2 and Quadratic. As shown in the figure, stochasticity parameters exhibit
certain relevance with respect to the galaxy properties, and this is a natural
outcome. On the one hand, the stochasticity traces the complex nature of
galaxy distribution, thus they should exhibit regular evolution to varying
degree in parametric space. On the other hand, as it is indicated by Eq. (8),
the stochasticity is coupled with large scale structure fluctuation in higher
order correlation. The coevolution effect also promotes the stochasticity
parameters to behave regularly. For the $(\alpha,\beta)$ parameters in Expan-
mm-2, the bins in different $\mathsf{FluxLimitLabel}$, equivalently different
$\bar{n}_{g}$, are separated into different tracks. It can be understood by
the form of Expan-mm-2, where the parameter $\alpha$ is degenerated with
$\bar{n}_{g}$ (as well as galaxy bias in fact), thus any variation in
stochasticity corresponding to number count will propagate to $\alpha$
directly though $\bar{n}_{g}$. It turns out as that different galaxy samples
populations are offset according to $\mathsf{FluxLimitLabel}$, modulate by
$\alpha$. Such obvious separation does not appear for $(c_{1},c_{2})$ in
Quadratic, as we expected.

We have emphasized the variation of stochasticity behavior for different
galaxy definitions and simulation details. We demonstrate the latter
intuitively by comparing the parameters distribution of two sets of simulation
samples in Fig. 10 The global features for two sets are similar, but the
details in each redshift is different to a large extent.

## References

  * Abareshi  _et al._ [2022] B. Abareshi, J. Aguilar, S. Ahlen, S. Alam, D. M. Alexander, R. Alfarsy, L. Allen, C. A. Prieto, O. Alves, J. Ameel, _et al._ , The Astronomical Journal 164, 207 (2022).
  * Adame  _et al._ [2023] A. Adame, J. Aguilar, S. Ahlen, S. Alam, G. Aldering, D. Alexander, R. Alfarsy, C. A. Prieto, M. Alvarez, O. Alves, _et al._ , arXiv preprint arXiv:2306.06307 (2023).
  * Collaboration  _et al._ [2023] D. Collaboration, A. Adame, J. Aguilar, S. Ahlen, S. Alam, G. Aldering, D. Alexander, R. Alfarsy, C. A. Prieto, M. Alvarez, _et al._ , arXiv preprint arXiv:2306.06308 (2023).
  * Ivezić  _et al._ [2019] Ž. Ivezić, S. M. Kahn, J. A. Tyson, B. Abel, E. Acosta, R. Allsman, D. Alonso, Y. AlSayyad, S. F. Anderson, J. Andrew, _et al._ , The Astrophysical Journal 873, 111 (2019).
  * Scaramella  _et al._ [2022] R. Scaramella, J. Amiaux, Y. Mellier, C. Burigana, C. Carvalho, J.-C. Cuillandre, A. Da Silva, A. Derosa, J. Dinis, E. Maiorano, _et al._ , Astronomy & Astrophysics 662, A112 (2022).
  * Cao  _et al._ [2018] Y. Cao, Y. Gong, X.-M. Meng, C. K. Xu, X. Chen, Q. Guo, R. Li, D. Liu, Y. Xue, L. Cao, _et al._ , Monthly Notices of the Royal Astronomical Society 480, 2178 (2018).
  * Gong  _et al._ [2019] Y. Gong, X. Liu, Y. Cao, X. Chen, Z. Fan, R. Li, X.-D. Li, Z. Li, X. Zhang, and H. Zhan, The Astrophysical Journal 883, 203 (2019).
  * Lin  _et al._ [2022] H. Lin, Y. Gong, X. Chen, K. C. Chan, Z. Fan, and H. Zhan, Monthly Notices of the Royal Astronomical Society 515, 5743 (2022).
  * Peacock and Bilicki [2018] J. Peacock and M. Bilicki, Monthly Notices of the Royal Astronomical Society 481, 1133 (2018).
  * Krolewski  _et al._ [2021] A. Krolewski, S. Ferraro, and M. White, Journal of Cosmology and Astroparticle Physics 2021 (12), 028.
  * White  _et al._ [2022] M. White, R. Zhou, J. DeRose, S. Ferraro, S.-F. Chen, N. Kokron, S. Bailey, D. Brooks, J. Garcia-Bellido, J. Guy, _et al._ , Journal of Cosmology and Astroparticle Physics 2022 (02), 007.
  * Marques  _et al._ [2024] G. Marques, M. Madhavacheril, O. Darwish, S. Shaikh, M. Aguena, O. Alves, S. Avila, D. Bacon, E. Baxter, K. Bechtol, _et al._ , Journal of Cosmology and Astroparticle Physics 2024 (01), 033.
  * Zhang [2006] P. Zhang, Astrophys. J. 647, 55 (2006), arXiv:astro-ph/0512422 .
  * Dong  _et al._ [2022] F. Dong, P. Zhang, Z. Sun, and C. Park, Astrophys. J. 938, 72 (2022), arXiv:2206.04917 [astro-ph.CO] .
  * Zhang  _et al._ [2007] P. Zhang, M. Liguori, R. Bean, and S. Dodelson, Phys. Rev. Lett. 99, 141302 (2007), arXiv:0704.1932 [astro-ph] .
  * Desjacques  _et al._ [2018] V. Desjacques, D. Jeong, and F. Schmidt, Physics reports 733, 1 (2018).
  * Bonoli and Pen [2009] S. Bonoli and U.-L. Pen, Monthly Notices of the Royal Astronomical Society 396, 1610 (2009).
  * Seljak  _et al._ [2009] U. Seljak, N. Hamaus, and V. Desjacques, Physical Review Letters 103, 091303 (2009).
  * Baldauf  _et al._ [2013] T. Baldauf, U. Seljak, R. E. Smith, N. Hamaus, and V. Desjacques, Physical Review D 88, 083507 (2013).
  * Hamaus  _et al._ [2010] N. Hamaus, U. Seljak, V. Desjacques, R. E. Smith, and T. Baldauf, Physical Review D 82, 043515 (2010).
  * Cai  _et al._ [2011] Y.-C. Cai, G. Bernstein, and R. K. Sheth, Monthly Notices of the Royal Astronomical Society 412, 995 (2011).
  * Liu  _et al._ [2021] Y. Liu, Y. Yu, and B. Li, The Astrophysical Journal Supplement Series 254, 4 (2021).
  * Tegmark and Bromley [1999] M. Tegmark and B. C. Bromley, The Astrophysical Journal 518, L69–L72 (1999).
  * Hamaus  _et al._ [2012] N. Hamaus, U. Seljak, and V. Desjacques, Physical Review D 86, 103513 (2012).
  * Zhou  _et al._ [2023a] S. Zhou, P. Zhang, and Z. Chen, Monthly Notices of the Royal Astronomical Society 523, 5789 (2023a).
  * Cabass and Schmidt [2020] G. Cabass and F. Schmidt, Journal of Cosmology and Astroparticle Physics 2020 (07), 051.
  * Eggemeier  _et al._ [2019] A. Eggemeier, R. Scoccimarro, and R. E. Smith, Physical Review D 99, 123514 (2019).
  * Alexander  _et al._ [2020] E. Alexander, C. Martin, P. Andrea, _et al._ , Phys. Rev. D 102 (2020).
  * Eggemeier  _et al._ [2021] A. Eggemeier, R. Scoccimarro, R. E. Smith, M. Crocce, A. Pezzotta, and A. G. Sánchez, Physical Review D 103, 123550 (2021).
  * Pezzotta  _et al._ [2021] A. Pezzotta, M. Crocce, A. Eggemeier, A. G. Sánchez, and R. Scoccimarro, Physical Review D 104, 043531 (2021).
  * Kokron  _et al._ [2022] N. Kokron, J. DeRose, S.-F. Chen, M. White, and R. H. Wechsler, Monthly Notices of the Royal Astronomical Society 514, 2198 (2022).
  * Schmidt [2016] F. Schmidt, Physical Review D 93, 063512 (2016).
  * Perko  _et al._ [2016] A. Perko, L. Senatore, E. Jennings, and R. H. Wechsler, arXiv preprint arXiv:1610.09321 (2016).
  * Chudaykin  _et al._ [2020] A. Chudaykin, M. M. Ivanov, O. H. Philcox, and M. Simonović, Physical Review D 102, 063533 (2020).
  * Ashley  _et al._ [2016] P. Ashley, S. Leonardo, J. Elise, _et al._ , arXiv preprint arXiv:1610.09321 (2016).
  * Moretti  _et al._ [2023] C. Moretti, M. Tsedrik, P. Carrilho, and A. Pourtsidou, arXiv preprint arXiv:2306.09275 (2023).
  * Mergulhão  _et al._ [2023] T. Mergulhão, H. Rubira, and R. Voivodic, arXiv preprint arXiv:2306.05474 (2023).
  * Schmittfull  _et al._ [2019] M. Schmittfull, M. Simonović, V. Assassi, and M. Zaldarriaga, Physical Review D 100, 043514 (2019).
  * Aghanim  _et al._ [2020] N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. Banday, R. Barreiro, N. Bartolo, S. Basak, _et al._ , Astronomy & Astrophysics 641, A6 (2020).
  * Tutusaus  _et al._ [2023] I. Tutusaus, C. Bonvin, and N. Grimm, arXiv preprint arXiv:2312.06434 (2023).
  * McDonald and Roy [2009] P. McDonald and A. Roy, Journal of Cosmology and Astroparticle Physics 2009 (08), 020.
  * Mirbabayi  _et al._ [2015] M. Mirbabayi, F. Schmidt, and M. Zaldarriaga, Journal of Cosmology and Astroparticle Physics 2015 (07), 030.
  * Zhang  _et al._ [2013] P. Zhang, J. Pan, and Y. Zheng, Physical Review D 87, 063526 (2013).
  * Zheng  _et al._ [2013] Y. Zheng, P. Zhang, Y. Jing, W. Lin, and J. Pan, Physical Review D 88, 103510 (2013).
  * Smith  _et al._ [2007] R. E. Smith, R. Scoccimarro, and R. K. Sheth, Physical Review D 75, 10.1103/physrevd.75.063512 (2007).
  * Cole  _et al._ [2005] S. Cole, W. J. Percival, J. A. Peacock, P. Norberg, C. M. Baugh, C. S. Frenk, I. Baldry, J. Bland-Hawthorn, T. Bridges, R. Cannon, M. Colless, C. Collins, W. Couch, N. J. G. Cross, G. Dalton, V. R. Eke, R. De Propris, S. P. Driver, G. Efstathiou, R. S. Ellis, K. Glazebrook, C. Jackson, A. Jenkins, O. Lahav, I. Lewis, S. Lumsden, S. Maddox, D. Madgwick, B. A. Peterson, W. Sutherland, and K. Taylor, Monthly Notices of the Royal Astronomical Society 362, 505–534 (2005).
  * Song  _et al._ [2015] Y.-S. Song, A. Taruya, E. Linder, K. Koyama, C. G. Sabiu, G.-B. Zhao, F. Bernardeau, T. Nishimichi, and T. Okumura, Physical Review D 92, 043522 (2015).
  * Bose  _et al._ [2023] B. Bose, P. Carrilho, M. Marinucci, C. Moretti, M. Pietroni, E. Carella, L. Piga, B. Wright, F. Vernizzi, C. Carbone, _et al._ , arXiv preprint arXiv:2311.13529 (2023).
  * Schmittfull and Seljak [2018] M. Schmittfull and U. Seljak, Physical Review D 97, 123540 (2018).
  * Korytov  _et al._ [2019] D. Korytov, A. Hearin, E. Kovacs, P. Larsen, E. Rangel, J. Hollowed, A. J. Benson, K. Heitmann, Y.-Y. Mao, A. Bahmanyar, _et al._ , The Astrophysical Journal Supplement Series 245, 26 (2019).
  * Springel  _et al._ [2017] V. Springel, R. Pakmor, A. Pillepich, R. Weinberger, D. Nelson, L. Hernquist, M. Vogelsberger, S. Genel, P. Torrey, F. Marinacci, and J. Naiman, Monthly Notices of the Royal Astronomical Society 475, 676 (2017).
  * Nelson  _et al._ [2017] D. Nelson, A. Pillepich, V. Springel, R. Weinberger, L. Hernquist, R. Pakmor, S. Genel, P. Torrey, M. Vogelsberger, G. Kauffmann, F. Marinacci, and J. Naiman, Monthly Notices of the Royal Astronomical Society 475, 624 (2017).
  * Pillepich  _et al._ [2017] A. Pillepich, D. Nelson, L. Hernquist, V. Springel, R. Pakmor, P. Torrey, R. Weinberger, S. Genel, J. P. Naiman, F. Marinacci, and M. Vogelsberger, Monthly Notices of the Royal Astronomical Society 475, 648 (2017).
  * Naiman  _et al._ [2018] J. P. Naiman, A. Pillepich, V. Springel, E. Ramirez-Ruiz, P. Torrey, M. Vogelsberger, R. Pakmor, D. Nelson, F. Marinacci, L. Hernquist, R. Weinberger, and S. Genel, Monthly Notices of the Royal Astronomical Society 477, 1206 (2018).
  * Marinacci  _et al._ [2018] F. Marinacci, M. Vogelsberger, R. Pakmor, P. Torrey, V. Springel, L. Hernquist, D. Nelson, R. Weinberger, A. Pillepich, J. Naiman, and S. Genel, Monthly Notices of the Royal Astronomical Society 10.1093/mnras/sty2206 (2018).
  * Zhou  _et al._ [2023b] R. Zhou, B. Dey, J. A. Newman, D. J. Eisenstein, K. Dawson, S. Bailey, A. Berti, J. Guy, T.-W. Lan, H. Zou, J. Aguilar, S. Ahlen, S. Alam, D. Brooks, A. de la Macorra, A. Dey, G. Dhungana, K. Fanning, A. Font-Ribera, S. G. A. Gontcho, K. Honscheid, M. Ishak, T. Kisner, A. Kovács, A. Kremin, M. Landriau, M. E. Levi, C. Magneville, M. Manera, P. Martini, A. M. Meisner, R. Miquel, J. Moustakas, A. D. Myers, J. Nie, N. Palanque-Delabrouille, W. J. Percival, C. Poppett, F. Prada, A. Raichoor, A. J. Ross, E. Schlafly, D. Schlegel, M. Schubnell, G. Tarlé, B. A. Weaver, R. H. Wechsler, C. Yéche, and Z. Zhou, The Astronomical Journal 165, 58 (2023b).
  * Raichoor  _et al._ [2023] A. Raichoor, J. Moustakas, J. A. Newman, T. Karim, S. Ahlen, S. Alam, S. Bailey, D. Brooks, K. Dawson, A. de la Macorra, A. de Mattia, A. Dey, B. Dey, G. Dhungana, S. Eftekharzadeh, D. J. Eisenstein, K. Fanning, A. Font-Ribera, J. García-Bellido, E. Gaztañaga, S. G. A Gontcho, J. Guy, K. Honscheid, M. Ishak, R. Kehoe, T. Kisner, A. Kremin, T.-W. Lan, M. Landriau, L. Le Guillou, M. E. Levi, C. Magneville, M. Manera, P. Martini, A. M. Meisner, A. D. Myers, J. Nie, N. Palanque-Delabrouille, W. J. Percival, C. Poppett, F. Prada, A. J. Ross, V. Ruhlmann-Kleider, C. G. Sabiu, E. F. Schlafly, D. Schlegel, G. Tarlé, B. A. Weaver, C. Yèche, R. Zhou, Z. Zhou, and H. Zou, The Astronomical Journal 165, 126 (2023).
  * Abbott  _et al._ [2022] T. Abbott, M. Aguena, A. Alarcon, S. Allam, O. Alves, A. Amon, F. Andrade-Oliveira, J. Annis, S. Avila, D. Bacon, _et al._ , Physical Review D 105, 023520 (2022).
  * Pezzotta  _et al._ [2023] A. Pezzotta, C. Moretti, M. Zennaro, A. M. Dizgah, M. Crocce, E. Sefusatti, I. Ferrero, K. Pardede, A. Eggemeier, A. Barreira, _et al._ , arXiv preprint arXiv:2312.00679 (2023).
  * Giannantonio  _et al._ [2016] T. Giannantonio, P. Fosalba, R. Cawthon, Y. Omori, M. Crocce, F. Elsner, B. Leistedt, S. Dodelson, A. Benoit-Lévy, E. Gaztanaga, _et al._ , Monthly Notices of the Royal Astronomical Society 456, 3213 (2016).
  * Omori  _et al._ [2019] Y. Omori, T. Giannantonio, A. Porredon, E. Baxter, C. Chang, M. Crocce, P. Fosalba, A. Alarcon, N. Banik, J. Blazek, _et al._ , Physical Review D 100, 043501 (2019).
  * Marques and Bernui [2020] G. A. Marques and A. Bernui, Journal of Cosmology and Astroparticle Physics 2020 (05), 052.
  * Hand  _et al._ [2018] N. Hand, Y. Feng, F. Beutler, Y. Li, C. Modi, U. Seljak, and Z. Slepian, The Astronomical Journal 156, 160 (2018).
  * Lewis and Challinor [2011] A. Lewis and A. Challinor, CAMB: Code for Anisotropies in the Microwave Background, Astrophysics Source Code Library, record ascl:1102.026 (2011), ascl:1102.026 .
  * Chisari  _et al._ [2019] N. E. Chisari, D. Alonso, E. Krause, C. D. Leonard, P. Bull, J. Neveu, A. 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, and E. L. Wagoner, The Astrophysical Journal Supplement Series 242, 2 (2019).
  * Springel [2010] V. Springel, Monthly Notices of the Royal Astronomical Society 401, 791 (2010).
  * Pillepich  _et al._ [2018] A. Pillepich, V. Springel, D. Nelson, S. Genel, J. Naiman, R. Pakmor, L. Hernquist, P. Torrey, M. Vogelsberger, R. Weinberger, _et al._ , Monthly Notices of the Royal Astronomical Society 473, 4077 (2018).
  * Weinberger  _et al._ [2016] R. Weinberger, V. Springel, L. Hernquist, A. Pillepich, F. Marinacci, R. Pakmor, D. Nelson, S. Genel, M. Vogelsberger, J. Naiman, _et al._ , Monthly Notices of the Royal Astronomical Society 465, 3291 (2016).
  * Heitmann  _et al._ [2019] K. Heitmann, H. Finkel, A. Pope, V. Morozov, N. Frontiere, S. Habib, E. Rangel, T. Uram, D. Korytov, H. Child, S. Flender, J. Insley, and S. Rizzi, The Astrophysical Journal Supplement Series 245, 16 (2019).
  * Jing [2019] Y. Jing, Science China Physics, Mechanics & Astronomy 62, 1 (2019).
  * Hinshaw  _et al._ [2013] G. Hinshaw, D. Larson, E. Komatsu, D. N. Spergel, C. Bennett, J. Dunkley, M. Nolta, M. Halpern, R. Hill, N. Odegard, _et al._ , The Astrophysical Journal Supplement Series 208, 19 (2013).