# Cosmology before noon with multiple galaxy populations

Haruki Ebina *  — and Martin White *  —

###### Abstract

Near-future facilities observing the high-redshift universe ($2<z<5$) will
have an opportunity to take advantage of “multi-tracer” cosmology by observing
multiple tracers of the matter density field: Lyman alpha emitters (LAE),
Lyman break galaxies (LBG), and CMB lensing $\kappa$. In this work we use
Fisher forecasts to investigate the effect of multi-tracers on next-generation
facilities. In agreement with previous work, we show that multiple tracers
improve constraints primarily from degeneracy breaking, instead of the
traditional intuition of sample variance cancellation. Then, we forecast that
for both BBN and CMB primary priors, the addition of lensing and LAEs onto a
LBG-only sample will gain 25% or more in many parameters, with the largest
gains being factor of $\sim 10$ improvement for $f_{\rm EDE}$. We include a
preliminary approach towards modelling the impact of radiative transfer (RT)
on forecasts involving LAEs by introducing a simplified model at linear theory
level. Our results, albeit preliminary, show that the while RT influences LAE-
only forecasts strongly, its effect on composite multi-tracer forecasts are
limited.

##  1 Introduction

The large-scale structure of the Universe arises from primordial fluctuations
in space-time in the very early Universe, acted upon by gravitational
instability in a dark-matter-dominated Universe over many Gyr [1, 2, 3, 4]. As
such it provides us with a window into fundamental physics and all of the
constituents, processes and parameters that influence the evolution of large-
scale structure even in subtle ways. In turn this has led to exciting
proposals for next-generation facilities that can make high precision
measurements of the quasi-linear modes that contain such a wealth of
information [5, 6, 7, 8, 9].

Much of the focus of these next generation facilities will be measurements of
large-scale structure at redshifts above $z\sim 1-2$ in order to take
advantage of several simultaneous, advantageous trends [10]. First, probing a
wide lever arm in time (redshift) is critical to disentangle evolutionary
influences from primordial ones and leads to rotated degeneracy directions
that tighten joint constraints. Second, the larger volume and smaller scale of
non-linearity both allow a larger lever arm in scale to break degeneracies
between cosmological and nuisance parameters. Third, the modes are in the
linear or quasi-linear regime and carry information from the early Universe
before it has been lost to non-linear evolution.

A fundamental limitation to how well we can probe fundamental physics and
cosmology comes from the influence of complex astrophysical processes that
affect the clustering of objects on the length scales that we observe. This
influence is usually described by a set of bias parameters [11] whose form is
dictated by the symmetries of our problem (e.g. rotational invariance, the
equivalence principle, …). The ability to break degeneracies between these
bias parameters and the fundamental parameters of interest is key to the
success of the program. Such degeneracy breaking is dramatically improved if
we have multiple tracers of the same underlying density fluctuations with
different (or no) biases, such as multiple galaxy populations (known as
“multi-tracer”; [12]) or gravitational lensing (of the CMB; [13, 14])).
Traditionally, multitracer has been considered in the context of ‘sample
variance cancellation’, which avoids the payment of sample variance by
measuring quantities, such as ratios of biases, at field level. However, for
our forecast we will not have high enough number densities for this to be
interesting, nor will we present forecasts for parameters that can be measured
at field level.

At $z>2$ there are two sets of galaxies that make natural targets for
spectroscopic surveys: Lyman break galaxies (LBGs) and Lyman alpha emitters
(LAEs). LBGs are massive, actively star-forming galaxies with a luminosity
approximately proportional to their stellar mass. They are abundant, well
studied and well understood [15, 16, 17]. By contrast LAEs are compact, metal-
poor, star-forming galaxies with SFR$\sim 1-10\ M_{\odot}\mathrm{yr}^{-1}$
that live in much lower mass halos and consequently are expected to have
significantly lower bias [18]. These two classes of targets naturally possess
different biases, making them ideal laboratories for the multi-tracer
approach. In addition we have lensing of the CMB, which provides a (projected)
tracer of the matter field itself (i.e. with no biases) which further serves
to break degeneracies.

In this paper we extend the Fisher formalism developed in ref. [5] to the case
of multiple galaxy populations (plus the CMB lensing already included in [5])
and forecast the improvement in constraining power that this provides. The
outline of the paper is as follows. In §2, we introduce and characterize the
three matter density tracers. In §3 we describe the theoretical framework used
for forecasting, including the modelling of biased tracers and the Fisher
information formalism. In §4 we then investigate the mechanism of multi-tracer
cosmology through forecasting of simplified galaxy tracers with varying
characteristics. §5 explains our approach to radiative transfer (RT), which
potentially poses concerns toward the use of LAEs as tracers. Finally in §6 we
present our forecasts for two experiment designs, Stage 4.5, motivated by DESI
II, and Stage V, and conclude in §7 with discussion.

##  2 Tracers of large-scale structure at high redshift

In this section we describe the relevant properties of the three tracers used
for this study; Lyman alpha emitters (LAEs), Lyman break galaxies (LBGs), and
gravitational lensing. For the galaxy tracers we determine the large-scale
bias, $b$, and comoving number density, $\mathrm{bar}{n}$, of the objects
based on previous astronomical surveys that have probed these tracers at the
redshifts of relevance. The distinction between these tracers is made
primarily by observational technique than physical properties, causing overlap
of these populations in many cases.

###  2.1 Lyman Alpha Emitters

Lyman alpha emitters (LAEs) are young, low-mass, actively star-forming
galaxies at high redshift [19, 18]. In the astronomical literature the
designation is primarily determined by the method of sample selection, and
LAEs are normally selected from narrow or intermediate band imaging surveys as
galaxies having large rest-frame equivalent width Ly$\alpha$ emission [20,
18]. While practical, this designation complicates the physical interpretation
of the population. Although the dominant physical origins of these emissions
(e.g. young massive stars, AGNs, and shock heating with minimal absorption by
dust) have been identified [18], a complete understanding of what determines
the high Ly$\alpha$ emissivity of LAEs has not yet been obtained. Although we
will largely follow the observational designation, from the perspective of
cosmology the primary interest will be in the ability to obtain a redshift for
the galaxy, which is more closely associated with line flux. There may also be
arguments for having detectable continuum emission in order to minimize
radiative transfer effects on the galaxy selection.

Due to the method of selection, the characterization of the LAE population for
the purpose of cosmological forecasts contains some subtleties and
difficulties. The majority of LAE surveys probe faint populations over
relatively small areas of sky, whereas the ideal objects for wide-field
spectroscopy tend to be brighter and not well-represented in such surveys.
Pure narrow-band selection is also not likely to be an efficient means of
finding galaxies over large cosmological volumes due to the limited comoving
depth each filter provides. This, along with varying observational definitions
of LAE across literature, means our estimates have an element of uncertainty.

The two key parameters for our forecasts will be the abundance and large-scale
clustering of LAEs, as determined by the 3D number density ($\mathrm{bar}{n})$
and large-scale bias ($b$). In order to determine these we make use of
published studies [21, 22, 23, 24, 25, 26, 27, 28, 29], which we discuss
further below. The final inputs are shown in Figure 1.

Figure 1: Left: 3D number density ($\mathrm{bar}{n}$) of LAEs found using
Schechter function fits to LF from previous studies [22, 23, 24]. The point
markers refer to the values reported in the original studies, while the lines
use interpolations of Schechter LF parameters. We see that although there is
rough consistency across the literature, a ‘factor of two’ difference is not
uncommon. We take the LF parameters from SC4K [22] as our fiducial model in
what follows. Right: Large scale bias $b$ of LAEs reported across the
literature with various redshift and line flux limits [21, 26, 30, 27, 28, 29,
25]. The solid lines show the fit of Eq. (2.4) for $f_{\rm lim}=3$, 5 and
$10\times 10^{-17}\,\text{erg}\,\text{cm}^{-2}\,\text{s}^{-1}$. The region
$b\geq 6$ is shaded, as populations with such high bias are not traditionally
used for cosmology.

LAE luminosity functions (LF) are frequently quoted in terms of the Ly$\alpha$
luminosity and fit with the Schechter functional form [31]. While the
uncertainties in the three parameters of the Schechter function are often
covariant, this information is not usually provided which makes assessing
uncertainties in $\mathrm{bar}{n}$ more difficult. Thus, we choose to compare
different surveys as an estimate of uncertainty. This has the advantage of
also including sample variance due to survey volume, though it unfortunately
mixes differences in observational selection into the uncertainty.

Assuming a Schechter LF the number density, $\mathrm{bar}{n}$, of objects
brighter than some minimum luminosity, $L_{\rm lim}$, can be written in terms
of three parameters ($\phi^{*}$, $L^{*}$ and $\alpha$) as

$\mathrm{bar}{n}(z)=\ln{10}\int^{\infty}_{\log_{10}L_{\rm lim}}\
\phi^{\star}\left(\frac{L}{L^{*}}\right)^{1+\alpha}e^{-L/L^{*}}d\log_{10}L$
(2.1)

where $L$ is the Ly$\alpha$ line luminosity. As it ties more directly to
observational ‘cost’, we will typically use the line flux ($f$) instead of
luminosity, the two are related by $L=4\pi f\,D_{L}^{2}(z)$, with $D_{L}$ the
luminosity distance.

In Figure 1, we compare number densities using various LF parameters from
different surveys [22, 24, 23]. While there is rough consistency, there are
frequently factors of $\sim 2$ between curves, suggesting our estimates of
$\mathrm{bar}{n}$ are likely uncertain ‘at the factor of two level’. We will
investigate the impact of this later in §6. We adopt the LF from the Slicing
COSMOS 4K (SC4K) survey [22] as our fiducial model. Since the LF-derived
number densities reflect the total number of galaxies at a given line flux
limit after correcting for incompleteness, we further include an estimated 50%
selection and redshift efficiency.

For our forecasts we will also need to understand the large-scale bias of
these populations, which is more complex than the number densities. Existing
surveys quote a range of different measures (see e.g. Fig. 3 of ref. [21] for
a recent compilation). The clustering strength in LAE surveys is often quoted
in terms of a correlation length, $r_{0}$, assuming a power-law (real-space)
correlation function:

$\xi_{gg}(r)=\left(\frac{r}{r_{0}}\right)^{\gamma}$ (2.2)

and the power-law index $\gamma$ is often held fixed (typically to $-1.8$) in
the fit. Given that astronomical LAE surveys typically have discrete narrow-
or medium-band selection over a small field, in practice $r_{0}$ is obtained
by fitting the angular correlation function, $w(\theta)$, to a power-law.
Ignoring redshift-space distortions (RSD) the projection of an isotropic 3D
correlation, $\xi(r)$, to 2D would give $w(\theta)\propto\theta^{\beta}$ with
$\beta=\gamma+1$. The presence of redshift-space distortions leads to
appreciable changes in this relation on scales comparable to the depth of the
shell in comoving distance, an effect that is usually neglected in the
literature.

Given the limited information available, we will convert the reported $r_{0}$
values to large-scale biases by working to linear order in bias, in which case

$b^{2}=\frac{\xi_{gg}}{\xi_{mm}}$ (2.3)

Again, any RSD contribution to the linear power-spectrum is ignored assuming
that the constraints are primarily from scales much smaller than the depth of
the survey. A subtlety involved here is that while we are using the linear
approximation for the bias, the power-law approximation for $\xi_{gg}$ assumes
nonlinear growth which is normally associated with non-linear, scale-dependent
bias. Although this poses an inconsistency in our approach, considering higher
order biases or re-investigating raw data to find the exact form of $\xi_{gg}$
goes well beyond the scope of this work and likely will not alter the
conclusions of this study significantly. It would be worthwhile to refine
these estimates when a new analysis of the clustering of LAEs over wide areas
becomes available.

This method for estimating the bias requires us to make a number of choices.
We evaluate $b$ by computing the ratio in Eq. 2.3 at a fixed scale, $r$. Since
the above estimation assumes that the power-law approximation of $\xi_{gg}$ is
valid we choose an angular scale of $\theta\sim 235^{\prime\prime}$. The
validity of the power-law approximation at this scale across surveys and
redshifts can be confirmed through inspection of the power-law fit to the
angular correlation in the literature [25, 21]. This scale corresponds to
$r\simeq 5\,h^{-1}\,\text{Mpc}$ at $z=3$, and $5\,h^{-1}$Mpc is large enough
to be in the quasi-linear regime for $z>2$ while being smaller than the
typical survey depth to reduce the impact of RSD. The numerator is estimated
from the power-law fit, while for the denominator we use the Hankel transform
of the HaloFit [32] power spectrum as computed by CLASS [33] which is
approximately a power-law on these scales.

We utilize clustering measurements from various surveys in the literature [21,
26, 30, 27, 28, 29, 25]. Note that we adopt $r_{0}$ values from [21] for all
surveys, as they have reported the values from the literature after
normalizing for some methodology differences. We estimate $b$ as described
above for a range of redshifts ($z$) and flux limits ($f_{\rm lim})$ and the
results are shown in Fig. 1. For convenience we fit these results with

$b(f_{\rm lim},z)=A(f_{\rm lim})(1+z)+B(f_{\rm lim})(1+z)^{2}$ (2.4)

where

$\displaystyle A(f_{\rm lim})$ $\displaystyle=0.457-1.755(\log_{10}f_{\rm
lim}+17)+0.720(\log_{10}f_{\rm lim}+17)^{2}$ (2.5) $\displaystyle B(f_{\rm
lim})$ $\displaystyle=0.012+0.318(\log_{10}f_{\rm
lim}+17)+0.043(\log_{10}f_{\rm lim}+17)^{2}$ (2.6)

are both quadratic functions of $\log_{10}{f_{\rm lim}}$, with the flux
measured in $\text{erg}\,\text{cm}^{-2}\,\text{s}^{-1}$ by standard
convention. The shape of this fit is motivated by a similar fit to the LBG
bias in ref. [17]. Since at these high redshifts the linear growth rate
$D(z)\propto(1+z)$ the $(1+z)$ term corresponds to co-evolution [34] while the
$(1+z)^{2}$ term is its first order correction. Ultimately, however, this is
simply a convenient fit to the observed results. The solid lines in Fig. 1
show the fit for $f_{\rm lim}=3$, 5 and $10\times
10^{-17}\,\text{erg}\,\text{cm}^{-2}\text{s}^{-1}$. Note that while we choose
to model $b$ as a function of $f_{\rm lim}$, which can be derived from $L_{\rm
lim}$, it is common for surveys to report the limiting luminosity in terms of
the median luminosity $L_{\rm med}$. The conversion from $L_{\rm med}$ to
$L_{\rm lim}$, assuming a Schechter luminosity function, is described in
Appendix B.

###  2.2 Lyman Break Galaxies

Figure 2: Left: Number density $\mathrm{bar}{n}$ of LBGs using Schechter
function fits to the UV LF [17], using fit parameters reported from preceding
studies [35, 36, 37]. Right: Large-scale bias $b$ of LBGs using a fit function
$b(m,z)$ in Eq. 2.9 [17]. The region $b\geq 6$ is shaded to indicate that
these high bias tracers are traditionally not used for cosmology.

Lyman break galaxies (LBGs) [17, 15, 16] are a set of massive, actively star-
forming galaxies at high redshift. The population is typically selected
through a ‘dropout’ technique in which the selection targets the distinct drop
in flux density $F_{\nu}$ that occurs at 912Å (the Lyman break) due to
absorption by stellar atmospheres and interstellar medium on top of an
otherwise flat spectrum. For the redshifts of our interest, the spectra
bluewards of 1216Å is further suppressed due to Lyman-series blanketing [38].
Importantly, selection for this drop yields physically similar galaxies across
redshift. Each selection sample is typically named after the filter that is
used to detect this drop. For instance, a $u$-dropout sample will have a lack
of flux in the $u$ band compared to significant flux in redder bands,
corresponding to galaxies at $z\approx 3$. The typical width for each sample
is $\Delta z\approx 0.7$. Similar to LAEs, the number density
$\mathrm{bar}{n}$ and large scale bias $b$ are our key parameters, allowing us
to use a nearly identical approach for parameter estimation. The results are
summarized in Figure 2.

The Schechter formalism is once again useful to predict the number density.
Here, we will be using the UV ($\approx 1500$Å rest-frame) LF with absolute
magnitude $M_{\rm UV}$

$dn(M_{\rm
UV})=\left(\frac{\ln{10}}{2.5}\right)\phi^{*}10^{-0.4(1+\alpha)(M_{\rm
UV}-M_{\rm UV}^{*})}\exp{-10^{-0.4(M_{\rm UV}-M_{\rm UV}^{*})}}dM_{\rm UV}$
(2.7)

with the LF parameters from Table 3 of ref. [17] (in turn based on refs. [35,
36, 37]). The absolute magnitude can be converted to apparent magnitude, $m$,
using the appropriate (near-zero) $k$-correction [17]

$M_{\rm
UV}=m-5\log_{10}\left(\frac{D_{L}(z)}{10\text{pc}}\right)+2.5\log_{10}(1+z)+\underbrace{m_{\rm
UV}-m}_{\approx 0}$ (2.8)

We will use this to model (apparent) magnitude limited surveys.

The LBG biases are calculated using the same power-law approximation of
correlation functions as in LAEs. Ref. [17] fit a function $b(m,z)$ to results
from the literature

$b(m,z)=A(m)(1+z)+B(m)(1+z)^{2}$ (2.9)

where $A(m)=-0.98(m-25)+0.11$ and $B(m)=0.12(m-25)+0.17$. Again, while this
fitting model has some physical basis on the co-evolution approximation of
galaxies, it is ultimately a useful fit to the data.

$z$ | 2.2 | 2.5 | 3.1 | 3.9 | 4.7 | 5.4
---|---|---|---|---|---|---
LAE $f_{\rm lim}^{*}$ | 13 | 11 | 5.6 | 5.1 | 5.4 | 6.9
$z$ | 2 | – | 3 | 3.8 | 4.9 | 5.9
LBG $m^{*}$ | 24.2 | – | 24.7 | 25.4 | 25.5 | 25.8
Table 1: Characteristic luminosity of Schechter function fits to the LAE LF
[22] and LBG LFs [17] (in turn based on [35, 36, 37]), converted to line flux
(in units of $10^{-17}\,\text{erg}\,\text{cm}^{-2}\,\text{s}^{-1}$) and
apparent magnitude, respectively. This can be used to gauge what the optimal
limiting luminosity is for experiments, as the LF will transition from an
exponential to a power-law fainter than these characteristic luminosities and
provide diminishing returns in galaxy counts.

The summary of the characteristic luminosities for LAEs and LBGs, shown in
Table 1 provide us with an intuitive limit on where our gains are maximized
when pushing to fainter luminosities. Starting at high luminosities, the
number densities increase exponentially as we push fainter until we reach the
“knee” at the characteristic luminosity ($L^{\star}$), where the increase
turns to a power-law. Near-term surveys will be able to reach such regimes.

Experiment | $2.4<z<3.2$ | $3.2<z<4.4$ | $f_{\rm sky}$
---|---|---|---
LAE $10^{17}f_{\rm lim}$ | LBG $m$ | LAE $10^{17}f_{\rm lim}$ | LBG $m$
Stage 4.5 | 10 | 24.5 | – | – | 0.1212
Stage V Plan A | 10 | 24.5 | 10 | 24.5 | 0.4
Stage V Plan B | 10 | 24.5 | 5 | 24.5 | 0.3211
Stage V Plan C | 10 | 24.5 | 10 | 25.0 | 0.2796
Table 2: Experiment designs considered for this work. The Stage 4.5 design is
strongly motivated by the proposed DESI II design.

For this work, we will consider a ‘Stage 4.5’ design motivated by the DESI II
design with $2.4<z<3.2$, and a ‘Stage V’ design with $2.4<z<4.4$. Both designs
are constrained by the limiting luminosities for LAEs and LBGs. While Stage
4.5, with a closer timeline, is specified in luminosity limits, we will
consider multiple scenarios of Stage V luminosity limits to perform a crude
optimization. The optimization will be done by assuming both a fixed total
integration time for the experiment and a fixed individual integration time
per galaxy, regardless of flux, magnitude, or galaxy type. This is equivalent
to fixing the total number of galaxies observed. The experimental designs are
shown in Table 2. A comparison against Table 1 shows that for a Stage 4.5
survey with $f_{\rm lim}=10\times 10^{-17}$, $m=24.5$, and $z\lesssim 3$, the
luminosity limits are near the characteristic luminosities with the exception
of LAEs at $z=3.1$.

A more quantitative measure of the strength of our signal is where in
$k$-space we become shot noise limited in the monopole, as shown in Figure 3.
For the Stage 4.5 survey, the LAE signals become shot noise limited at
$k\simeq 0.15\,h\,\mathrm{Mpc}^{-1}$ and LBGs at $k\simeq
0.3-0.4\,h\,\mathrm{Mpc}^{-1}$. Furthermore, if we extend these luminosity
limits to $z=4$, we are totally shot noise dominated for LAEs. Notice that the
magnitude limits on the right axis do not go to 26 for $z=2.5$ and 3.0. This
is because the bias-corrected shot noise
$[\mathrm{bar}{n}(b^{2}+2bf/3+f^{2}/5)]^{-1}$ increases at fainter magnitudes
due to the decrease in bias over-powering the increase in number density as we
push fainter than the corresponding characteristic magnitudes of the LBG
Schechter function fit.

Figure 3: Monopole of the signal ($P-P_{\rm stoch}$) against the shot noise at
various luminosity limits (right axis) for LAEs (top) and LBGs (bottom). The
$k$ where signal equals shot noise reflects where we enter a shot noise
dominated regime ($\bar{n}P$=1). The bias-correction simply refers to dividing
out the monopole Kaiser factor, $b^{2}+2bf/3+f^{2}/5$, so that the large-scale
clustering is independent of magnitude. The spread in the curves at high-$k$
is due to scale-dependent bias, which we model as described in §3. The dotted
lines illustrate where the fiducial Stage 4.5 survey become shot noise
dominated.

###  2.3 CMB Lensing

We will also include CMB lensing $\kappa$ as a third tracer of the matter
density. The aspect of lensing that makes it a particularly interesting
addition to the set of tracers is that it introduces fewer nuisance parameters
than a galaxy tracer; it is an unbiased tracer of matter, with no shot noise.
The high redshifts of near-term experiments further motivate the use of
lensing, as the CMB lensing kernel peaks near these redshifts. On the other
hand, we will not be including cosmic shear as a tracer, as this will require
galaxy sources at extremely high redshift.

We consider Simons Observatory (SO) [39] with $f_{\rm sky}=0.4$ as the CMB
experiment for our forecasts and use the noise $N^{\kappa}_{\ell}$ obtained
from the SO noise calculator as in ref. [5]. The calculator uses an iterated
minimum variance (MV) quadratic estimator and MV internal linear combination
of both CMB temperature and polarization data for calculation, and includes
atmospheric effects [5].

##  3 Theory and formalism

In this section we outline the theoretical models we employ and the
forecasting formalism we adopt, which is based on the Fisher matrix and is
substantially similar to that used111https://github.com/NoahSailer/FishLSS in
ref. [5]. We take $\Lambda$CDM as our fiducial cosmology with parameters given
in Table 3.

###  3.1 Model

As in ref. [5], our theoretical model is based upon 1-loop perturbation theory
with a general, quadratic bias model. This model has been shown to be accurate
(at the relevant scales) well beyond the demands of current and next-
generation surveys [40] while including scale-dependent bias and mode-coupling
due to quasi-linear evolution that improve the reliability of the forecasts.
It has been applied to earlier surveys in ref. [41, 42, 43] and is currently
being evaluated for use in the DESI and Euclid surveys [44]. For further
details of the model see [45, 46, 40].

Within this formalism, the auto- or cross-spectra for our tracers can be
written as a linear combination of different ‘basis’ spectra multiplied by
bias coefficients, plus contributions from counterterms and stochastic terms:

$P(k,\mu)=\sum_{X,Y}b_{X}b_{Y}P_{XY}(k,\mu)+P_{c.t.}(k,\mu)+P_{\rm
stoch}(k,\mu)$ (3.1)

where
$X,Y\in\\{m,\delta_{L},\delta_{L}^{2},s_{L}^{2},\nabla^{2}\delta_{L},\cdots\\}$
label the component spectra, $P_{c.t.}$ are counterterms arising from small-
scale physics not directly included in the model, and $P_{\rm stoch}$ are
stochastic terms such as shot noise and fingers of god. Figure 4 shows the
contribution of each component in a sample case of constant biases, $b_{a}=2$
and $b_{b}=5$, constant number density, $\bar{n}=(100b^{2})^{-1}$, and no
correlation between the shot noise in the two samples (these assumptions
correspond to $\tilde{N}=100$ and $f_{\rm over}=0$, as introduced later in
this section). Note how the scale and angle dependence of each component
spectra differ, allowing us to potentially break degeneracies between
parameters.

Figure 4: The contributions to the power spectra from each component in Eq.
(3.1), where $b_{a}=2$, $b_{b}=5$, $z=3$, and $P_{\rm PT}$ refers to the sum
of quadratic bias component spectra $\sum_{X,Y}b_{X}b_{Y}P_{XY}(k,\mu)$. The
columns refer to the monopole, quadrupole, and hexadecapole of each component.
Red lines show the $aa$ auto-spectra, green lines the $ab$ cross-spectra, and
blue lines the $bb$ auto-spectra. The shaded region reflects the $k_{\rm nl}$
limit at $z=3$.

###  3.2 Parameters

The parameters that we use can be categorized into ‘global’ and ‘local’
parameters. The global parameters are those that are common to all redshifts,
such as cosmological parameters. The local parameters are those that are
assigned for each redshift that we consider. The perturbation theory
parameters that we describe here ($b_{i}$, $N_{i}$, $\alpha_{i}$) are examples
of local parameters. For a complete list of parameters, their categorization
(global or local), and their fiducial values, refer to Table 3.

The bias parameters $b_{X}$ contribute to the spectra are organized in terms
of the dependence of the galaxy overdensity field $\delta_{g}$ on the matter
field $\delta$ and shear $s_{ij}$. Up to one-loop, the bias terms that we
include in our model are

$\delta_{g}=b_{1}\delta+\frac{b_{2}}{2}\delta^{2}+b_{s}s^{2}$ (3.2)

where $s^{2}=(s_{ij})^{2}$ [45]. The cubic bias terms are nearly degenerate
with our counter- and stochastic terms and so do not need to be included
explicitly. This expression corresponds to the parametrization used in
Eulerian perturbation thoery (EPT). Each Eulerian bias can be related to the
Lagrangian biases $b_{X}^{L}$ as [11, 45]

$b_{1}=b^{L}_{1}+1,\quad b_{2}=b^{L}_{2}+\frac{8}{21}b^{L}_{1},\quad
b_{s}=b^{L}_{s}-\frac{2}{7}b^{L}_{1}\quad.$ (3.3)

The fiducial values for each parameter is dependent on that of $b_{1}$, which
we have estimated for both LAEs and LBGs in §2. For $b_{2}$, we use the Sheth-
Tormen mass function [47] and peak-background split model, for which [48, 49]

$\displaystyle b^{L}_{1}$
$\displaystyle=\frac{1}{\delta_{c}}\left[a\nu^{2}-1+\frac{2p}{1+(a\nu^{2})^{p}}\right]$
(3.4) $\displaystyle b^{L}_{2}$
$\displaystyle=\frac{1}{\delta_{c}^{2}}\left[a^{2}\nu^{4}-3a\nu^{2}+\frac{2p(2a\nu^{2}+2p-1)}{1+(a\nu^{2})^{p}}\right]$
(3.5)

with $\delta_{c}=1.686$, $a=0.707$ and $p=0.3$. This is shown to agree well
with N-body simulations for $b_{2}$ [50] and thus motivates a reasonable
fiducial value for our samples. Furthermore, we have confirmed that our
results are qualitatively invariant to our choice of $b_{2}$. For $b_{s}$, we
use the coevolution approximation, which assumes $b^{L}_{s}=0$. While this
does not agree as closely with simulations as our choice for $b_{2}$,
simulations suggest that for of our samples we expect $|b^{L}_{s}|\lesssim 1$
[50], and the differences have little impact on our forecasts.

Our model for the counterterm consists of parameters $\alpha_{2n,a}$, such
that

$P_{a}\supset\alpha_{2n,a}\,\mu^{2n}(k/k_{\star})^{2}P_{cb}$ (3.6)

for each auto-spectrum, with $k_{\star}=1\,h\,\text{Mpc}^{-1}$. We inherit the
fiducial values of galaxy $\alpha$’s from ref. [5], where $\alpha_{0}$ values
are determined using fits towards Halofit simulations and higher order
parameters are set to 0 for convenience. For $\alpha_{0,\kappa}$, we obtain
fiducial values by simply fitting $P_{HF}(z)=P_{\rm
lin}(1+\alpha_{0,\kappa}k^{2})$. Cross-spectra counterterm contributions are
fully determined by those of the auto-spectra, thus we do not have degrees of
freedom for $\alpha_{ab,a\neq b}$, for example. A detailed discussion of
cross-spectrum counterterms is provided in Appendix C.

For the stochastic terms, we use a set of parameters $N_{2n,ab}$, such that

$P_{ab}\supset N_{2n,ab}(k\mu)^{2n}$ (3.7)

Our fiducial values of the parameters are survey dependent for $n=0$ or 1, and
zero for higher orders. Note that $n=0$ is the “shot noise term”, and hence
our fiducial value is $N_{0,aa}=1/n_{a}$ for auto-spectra. Extending this
model to cross-spectra, we let $N_{0,ab}=f_{\rm over}/\sqrt{n_{a}n_{b}}$
fiducially for $a\neq b$, where $f_{\rm over}\in[0,1]$ is a measure of how
much correlation we have between the tracers’ shot noise (e.g. to what degree
the populations occupy the same halos). As the fiducial value for $N_{2}$, we
set $N_{2,ab}=-N_{0,ab}\sigma_{v}^{2}$, with $\sigma_{v}=(1+z)(100{\rm
km/s})/H(z)$ being the typical velocity dispersion in comoving units [5].
Notice that the stochastic parameters are the only new degrees of freedom
introduced when including the cross-spectra, motivating the use of all
available spectra.

Our $k_{\rm max}$ is set by the non-linear scale. Unless explicitly stated, we
will use $k_{\rm max}(z)=k_{\rm nl}(z)=\Sigma^{-1}(z)$ where $\Sigma$ is the
rms displacement in the Zel’dovich approximation. We take $k_{\rm
min}=0.003\,h\,\text{Mpc}^{-1}$, motivated by the longest mode that fit in our
surveys, i.e. $k_{\rm min}\sim 2\pi/V^{1/3}$. Our results are in general not
very sensitive to $k_{\rm min}$.

###  3.3 2D power spectrum

To include CMB lensing in our analysis, we must also have a method to
calculate the 2D power-spectra, where we can cross-examine the LSS and CMB
data on equal footing. The 2D power-spectra can be calculated from the 3D
power-spectra using the Limber approximation,

$C_{\ell}^{ab}=\int_{0}^{\chi_{*}}d\chi\frac{W^{a}(\chi)W^{b}(\chi)}{\chi^{2}}P_{ab}\left(k_{\perp}=\frac{\ell+1/2}{\chi},k_{\parallel}=0\right)$
(3.8)

where the $a$ and $b$ are available tracers, either $\kappa$ or galaxies.
$\chi_{*}$ is the comoving distance to the surface of last scattering and
$W$’s are the lensing kernels

$W^{\kappa}(\chi)=\frac{3}{2}\Omega_{m}H_{0}^{2}(1+z)\frac{\chi(\chi_{*}-\chi)}{\chi_{*}}\quad\text{and}\quad
W^{g_{a}}(\chi)\propto H(z)\frac{dN_{a}}{dz}$ (3.9)

where $W^{g}$ is normalized such that $\int_{\chi}W^{g}(\chi)=1$. Note that we
take the lowest order correction to the Limber approximation
$\ell\rightarrow\ell+1/2$, improving our precision to $\mathcal{O}(\ell^{-2})$
[51, 5]. We do not marginalize over higher order stochastic and counter-term
parameters for the 2D power-spectra, as the Limber approximation limits us to
$\mu=0$. For the lensing-lensing auto-spectrum, the $\ell_{\rm max}$ is
limited not by $k_{\rm nl}$, but by a stricter constraint from baryonic
feedback, which limits us to $\ell<500$. For the rest of the spectra, we
translate our $k_{\rm max}$ constraints into $\ell_{\rm max}$.

###  3.4 Fisher formalism

The Fisher matrix formalism [52] allows us to forecast the parameter
constraints on an assumed underlying model, given a experimental
configuration. It is widely used to understand the impact of survey design
choices and parameter degeneracies on the inference of cosmological parameters
[53, 54, 55, 56, 5, 57]. The Fisher matrix can be calculated as

$F_{ij}=\sum_{X,Y}\int\frac{k^{2}\,dk\,d\mu}{2\pi^{2}}\ \frac{\partial
P_{X}}{\partial\theta_{i}}\mathbf{C}^{-1}_{XY}\frac{\partial
P_{Y}}{\partial\theta_{j}}$ (3.10)

where $X,Y$ label available auto- and cross-spectra, $\theta_{i}$ are
parameters, $\mathbf{C}^{-1}_{XY}=(\mathbf{C}^{-1})_{XY}$ and the covariance
matrix $\mathbf{C}$ is

$\mathbf{C}_{abcd}=P_{ac}P_{bd}+P_{ad}P_{bc}$ (3.11)

where the subscripts represent the samples used to calculate spectra. Note
that the power-spectra ($P_{X}$, $P_{ab}$) above include the counter and
stochastic terms. We compute these matrices using the FishLSS code, which was
originally developed for single-tracer forecasts [5] and extended to
multitracer forecasts for this work. In the development process we have
confirmed that the new version of the code reproduces previous results in
single-tracer cases. Without full consideration of multi-tracer effects, one
would calculate the Fisher matrix with $X,\,Y$ summing over only the auto-
spectra and thus lose information gained from the cross-spectra.

For a general galaxy-galaxy-lensing sample, the formalism is modified
slightly. While the Fisher matrix earlier still holds, one must compute the
matrix using the angular power-spectra $C_{\ell}^{X}$ as well. One obtains the
full set of information by combining the Fisher matrices from the 3D power-
spectra and 2D power-spectra. By restricting the full-shape signal to
$k_{\parallel}>10^{-3}$, there will be no covariance between the full-shape
and the angular spectra that probe $k_{\parallel}\sim 0$. This will allow for
the Fisher matrices to simply add together [5, 52].

The covariance matrix for a galaxy-galaxy-lensing cross-analysis is

$(\mathbf{C}_{\ell})_{XY}=\left\\{\begin{array}[]{ll}2(C_{\ell}^{\kappa\kappa}+N_{\ell}^{\kappa\kappa})^{2}/f_{\rm
sky}^{\rm CMB}&X=Y=\kappa\kappa\\\ (C_{\ell}^{\kappa g_{a,i}}C_{\ell}^{\kappa
g_{b,j}}+\delta_{ij}(C_{\ell}^{\kappa\kappa}+N_{\ell}^{\kappa\kappa})C_{\ell}^{g_{a,i}g_{b,i}})/f_{\rm
sky}^{\cap}&X=\kappa g_{a,i},Y=\kappa g_{b,j}\\\
\delta_{ij}(C_{\ell}^{g_{a,i}g_{c,i}}C_{\ell}^{g_{b,i}g_{d,i}}+C_{\ell}^{g_{a,i}g_{d,i}}C_{\ell}^{g_{b,i}g_{d,i}})/f_{\rm
sky}^{\rm LSS}&X=g_{a,i}g_{b,i},Y=g_{c,j}g_{d,j}\\\
2(C_{\ell}^{\kappa\kappa}+N_{\ell}^{\kappa\kappa})C_{\ell}^{\kappa
g_{a,i}}/f_{\rm sky}^{\rm CMB}&X=\kappa\kappa,Y=\kappa g_{a,i}\\\
2C_{\ell}^{\kappa g_{a,i}}C_{\ell}^{\kappa g_{b,i}}f_{\rm sky}^{\cap}/(f_{\rm
sky}^{\rm CMB}f_{\rm sky}^{\rm LSS})&X=\kappa\kappa,Y=g_{a,i}g_{b,i}\\\
\delta_{ij}(C_{\ell}^{\kappa
g_{b,i}}C_{\ell}^{g_{a,i}g_{c,i}}+C_{\ell}^{\kappa
g_{c,i}}C_{\ell}^{g_{a,i}g_{b,i}})/f_{\rm sky}^{\rm LSS}&X=\kappa
g_{a,i},Y=g_{b,j}g_{c,j}\\\ \end{array}\right.$ (3.12)

and the Fisher matrix is

$F^{\rm 2D}_{ij}=\sum_{\ell=\ell_{\rm min}}^{\ell_{\rm
max}}(2\ell+1)\sum_{XY}\frac{\partial
C_{\ell}^{X}}{\partial\theta_{i}}(\mathbf{C}_{\ell})^{-1}_{XY}\frac{\partial
C_{\ell}^{Y}}{\partial\theta_{j}}$ (3.13)

where $f_{\rm sky}^{\cap}$ is the overlapping sky area between LSS and CMB
surveys, and $X,Y$ sum over the available spectrum $\\{\kappa\kappa,\kappa
g_{a,i},g_{a,i}g_{b,i}\\}$ for all available galaxy tracers $a,b$ and redshift
bins $i$ [5]. For the purposes of our forecast we assume maximal overlap
between LSS and CMB experiments, i.e. $f_{\rm sky}^{\cap}=\min(f_{\rm
sky}^{\rm CMB},f_{\rm sky}^{\rm LSS})$. Note that the cross-spectra between
tracers at different redshift bins is zero, so we do not sum over them.

The simple treatment of priors is a strength of the Fisher formalism, with the
inverse covariance of the priors simply treated as an added information
matrix. We will be using conservative priors on the nuisance parameters to
regularize forecast behavior. The set of priors used is shown in Table 4. Most
of the nuisance parameter priors reflect the fidelity we have on constraining
these parameters, by applying uncertainties proportional to their fiducial
value. This cannot be done for parameters with a fiducial value of 0, so we
describe this briefly. For $N_{ab,a\neq b}$ and $N_{2,ab,a\neq b}$, we apply
the same errors as those of the auto-spectra, assuming there is full
correlation between the shot noise of two galaxy populations (e.g.
$N_{ab,a\neq b}=1/\sqrt{n_{a}n_{b}}$). The other terms with fiducial value
have no prior, due to the lack of conservative expectation for their values.
We have confirmed that our conservative priors have little effect on our final
results. We will be considering both BBN and Planck priors for cosmological
parameters.

###  3.5 Benefits of multi-tracer cosmology

The use of multiple tracers is closely associated with the idea of ‘sample
variance cancellation’ [12, 58]. For two quantities, e.g. a density and
velocity field, some relationships can be measured at the field level and thus
do not pay any sample variance penalty. For example, within linear theory and
in general relativity, the density contrast $\delta(\mathbf{k})$ is related to
the velocity divergence, $\theta(\mathbf{k})$ mode-by-mode. A ratio of
measured modes can thus measure $f$ with no sample variance penalty.
Similarly, for linear theory and scale-independent, deterministic bias the
ratio of $\delta_{i}(\mathbf{k})$ for two galaxy samples measures the ratio of
their biases in a way that is not limited by sample variance. Colloquially we
say that sample variance has ‘canceled’, though it would be more correct to
state that sample variance only enters when comparing a particular realization
to the ensemble average. While in general, cosmology dependence enters
quantities (such as the linear power spectrum) that are expectation values,
for some situations (such as the scale-dependent bias introduced by primordial
non-Gaussianity) this is not the case and ‘cancellation’ can lead to
tremendous gains in precision.

Sample variance cancellation can be recast as a fitting problem, in which case
it becomes the statement that when including cross-correlations in addition to
auto-correlations noise that is highly correlated between probes does not get
double counted [55]. In such situations some eigenvectors ‘pay’ the sample
variance penalty while others do not.

Past applications of sample variance cancellation [59] have not found drastic
gains in constraining power from easily achievable number densities and we are
unlikely to have high-enough number densities for this either. However,
another use of multiple tracers is to break degeneracies between parameters.
Including multiple tracers allow us to distinguish between parameter
dependencies of the data, breaking parameter degeneracies that limit our
precision. Recent studies indicate that in a general case, it is this
degeneracy breaking that drives the improvements in constraints, rather than
sample variance cancellation [60, 61].

Parameter | Definition | Fiducial value
---|---|---
$\log A_{s}$ | Primordial amplitude | $2.10732\times 10^{-9}$
$n_{s}$ | Spectral index | 0.96824
$h$ | Hubble parameter: $H_{0}/$100 km/s/Mpc | 0.6770
$\tau$ | Optical depth to reionization | 0.0568
$\omega_{c}$ | fractional dark matter density: $\Omega_{c}h^{2}$ | 0.11923
$\omega_{b}$ | fractional baryonic matter density: $\Omega_{b}h^{2}$ | 0.02247
$b$ | First order bias $\delta_{g}\supset b_{g}\delta_{m}$ |  $b_{\rm LAE}$: Eq. 2.4
|  |  $b_{\rm LBG}$: Eq. 2.9 [17]
$b_{2}$ | Second order bias | Eq. 3.5
$b_{s}$ | Shear bias | $-\frac{2}{7}(b_{g}-1)$
$\alpha_{0}$ | Lowest order counter term $P_{aa}\supset\alpha_{0,a}k^{2}P_{m}$ | $1.22+0.24b_{a}^{2}(z-5.96)$
$\alpha_{2n}$ | Higher order counter terms $P_{aa}\supset\alpha_{2n,a}\mu^{2n}k^{2}P_{m}$ | $0$
$N_{ab}$ | Shot noise term $P_{ab}\supset N_{ab}$ | $N_{aa}=1/\bar{n_{a}}$
|  | $N_{ab,a\neq b}=0$
$N_{2,ab}$ | FOG-like effect $P_{ab}\supset N_{2,ab}(k\mu)^{2}$ | $-N_{ab}\sigma_{v}^{2}$
$N_{4,ab}$ | FOG-like effect $P_{ab}\supset N_{4,ab}(k\mu)^{4}$ | 0
Table 3: The parameters for our fiducial cosmological model, $\Lambda$CDM with massive neutrinos. The cosmological parameters are placed on the top and model parameters are on the bottom. The model parameters are individually assigned either to each tracer or each spectra. Furthermore, we allow model parameters to evolve over redshift, which gives us new degrees of freedom for each additional $z$ bin.  Parameter | Uncertainty | Parameter | Uncertainty | Parameter | Uncertainty
---|---|---|---|---|---
$b$ | 100% | $N_{aa}$ | 100% | $\alpha_{0}$ | 100%
$b_{2}$ | 100% | $N_{ab,a\neq b}$ | $\sqrt{N_{aa}N_{bb}}$ | $\alpha_{2}$ | $\infty$
$b_{s}$ | 100% | $N_{2,aa}$ | 300% | $\alpha_{4}$ | $\infty$
$N_{4}$ | $\infty$ | $N_{2,ab,a\neq b}$ | $3\sqrt{N_{2,aa}N_{2,bb}}$ |  |
Table 4: Prior uncertainties for nuisance parameters. For $\Lambda$CDM
parameters, we use BBN or CMB primary priors appropriately.

##  4 Multitracer

Before we turn to forecasts for particular models let us investigate some test
samples with semi-realistic parameters to gain a better understanding of the
multitracer approach. To do this, we make forecasts for the amplitude of the
linear theory power spectrum, $\sigma_{8}$, holding fixed the other
$\Lambda$CDM parameters $\\{h,\tau,n_{s},\omega_{c},\omega_{b}\\}$ at some
fiducial values. We allow the biases, counterterms and stochastic terms for
each sample to vary.

Throughout this section we use the notation $a$ and $b$ to refer to the auto-
spectrum analyses of each tracer, “auto” or $aa+bb$ refers to using two tracer
populations, but with only the galaxy auto-spectra considered (i.e. incomplete
multitracer) and “full” or $a\times b$ refers to the complete multitracer
approach including the galaxy-galaxy cross-spectrum.

###  4.1 Linear

We will start by investigating multitracer cosmology in the Fisher formalism
using linear theory. In linear theory, with the exception of shot noise, the
power spectra of each galaxy population is proportional to the linear matter
power spectrum

$P_{ab}(k,\mu)=(b_{a}+f\mu^{2})(b_{b}+f\mu^{2})P_{L}(k)+N_{ab}$ (4.1)

where $b$ is the (deterministic, scale-independent) linear bias, $f$ is the
scale-independent growth factor, $P_{L}$ is the linear theory matter power
spectrum, $N_{ab}$ is the shot noise term, and $\mu$ is the cosine of the
angle to the line of sight. The relevant derivatives for linear theory power
spectrum are expressed as

$\begin{split}\frac{\partial
P_{ab}}{\partial\log\sigma_{8}}&=2(b_{a}+f\mu^{2})(b_{b}+f\mu^{2})P_{L}\\\
\frac{\partial P_{ab}}{\partial
b_{c}}&=[\delta_{ac}(b_{b}+f\mu^{2})+\delta_{bc}(b_{a}+f\mu^{2})]P_{L}\\\
\frac{\partial P_{ab}}{\partial N_{cd}}&=\delta_{ab,cd}\end{split}$ (4.2)

We use this method first as a check on our code (since the calculation can be
done analytically) and second to illustrate some of the features of
multitracer cosmology.

The Fisher information is a function of both the amount of data and the
degeneracy between the derivatives of the data with respect to the parameters.
In multitracer cosmology, we hope to see an increase in constraints that goes
beyond the simple increase in data quantity, which must come from breaking of
degeneracy. However, in linear theory, the lack of parameter dependencies
beyond the $\sigma_{8}$-$b$ degeneracy suggests that this will not occur.

At the minimum, we expect the low-bias tracer to improve the constraints by a
significant margin, as low-bias tracers have an inherent advantage in
constraining structure growth compared to high-bias tracers. Even if there is
no synergy between the measurement, matching the low-bias constraints will
significantly reduce the error from having a single, high-bias tracer. This
can easily be explained in the following sense. In linear theory, the auto-
spectrum of a biased tracer is given by

$P(k,\mu)=(b+f\mu^{2})^{2}P_{L}(k)+N=(b\sigma_{8}+f\sigma_{8}\mu^{2})^{2}\hat{P_{L}}(k)+N$
(4.3)

where $\hat{P_{L}}$ is the normalized linear power-spectrum,
$P_{L}=\sigma_{8}^{2}\hat{P_{L}}$. As $b$ is a nuisance parameter, one must
probe for the $\mu$-dependence of the spectrum to obtain $\sigma_{8}$. With
the error proportional to the overall spectrum mode-by-mode ($\delta P\propto
P$), one loses constraining power on $\sigma_{8}$ with greater $b$. To
maximize SNR, it’s clear that we must minimize $b$ and maximize $\bar{n}$,
both of which are satisfied by choosing a low-bias sample. In the limit that
we go to $b=0$, we are measuring the velocity divergence $\theta$. In this
case, our precision of the amplitude of $P$ is sample variance limited such
that $\delta A/A=\sqrt{2/N_{\rm modes}}$, making
$\delta\sigma_{8}/\sigma_{8}=(1/2)\sqrt{2/N_{\rm modes}}$.

Before proceeding further, let us define an effective shot noise $\tilde{N}$,
as

$b^{2}\tilde{N}\equiv\frac{1}{\bar{n}}$ (4.4)

We will use $\tilde{N}$ to explore different noise regimes throughout this
section. For reference, at $z=3$, LAEs with $f_{\rm lim}$ of 3, 5, 10$\times
10^{-17}\,\text{erg}\,\text{cm}^{-2}\,\text{s}^{-1}$ have $\tilde{N}\approx$
100, 170, 410, and LBGs with $m=$ 24, 24.5, 25 have $\tilde{N}\approx$ 570,
230, 140, respectively, in units of $h^{3}\,\text{Mpc}^{-3}$.

Let us consider samples in linear theory, with constant biases,
$b_{a}=\\{1.5,2,3,4\\}$, $b_{b}=5$, and a linear regime $k_{\rm
max}=0.1\,h\,\text{Mpc}^{-1}$. A subset of our constraint calculations are
shown in Figure 5. As we predict there is little to no improvement from adding
the cross-spectrum, as shown in Figure 6, indicating that we are losing little
to degeneracy.

In the same figure, we detect that the ‘auto’ ($aa+bb$) results are
significantly better than the single low-bias result. This indicates that a
naive error propagation of $1/\sigma^{2}=1/\sigma_{1}^{2}+1/\sigma_{2}^{2}$ is
not appropriate, as this predicts the low-bias and ‘auto’ constraints to be
similar for most of our samples. The reason is that we must include the
correlation between measurements (e.g. the velocity part of the power spectra
is the same). Including the correlation, $r$, the variance $\sigma^{2}$
becomes

$\sigma^{2}=\frac{(1-r^{2})\sigma_{1}^{2}\sigma_{2}^{2}}{\sigma_{1}^{2}-2r\sigma_{1}\sigma_{2}+\sigma_{2}^{2}}$
(4.5)

for the composite measurement, with $\sigma_{i}$ referring to the uncertainty
of the $i$’th measurement. Refer to Appendix A for detailed calculations.

This result explains the two main trends that we see. The first trend, the
decreasing improvement as one takes $b_{a}\rightarrow b_{b}$, can be
illustrated by the case with redundant measurements ($\sigma_{1}=\sigma_{2}$,
$r\rightarrow 1$), which yields $\sigma^{2}=\sigma_{1}^{2}=\sigma_{2}^{2}$.
The second trend, the increasing improvement as one goes to lower noise,
follows the limit of non-redundant measurements with high correlation
($\sigma_{1}\neq\sigma_{2}$, $\absolutevalue{r}\rightarrow 1$), which yields
$\sigma^{2}\rightarrow 0$.

One can directly compare the $\delta\sigma_{8}/\sigma_{8}$ results with
improvements in constraints for the linear bias ratio $b_{a}/b_{b}$, where
sample variance cancellation is present [12]. A direct comparison is made in
Figure 6. The lack of improvement in $\delta\log\sigma_{8}$ by including the
cross-correlation, as compared to the significant improvement in
$\delta(b_{a}/b_{b})$, illustrates that sample variance cancellation is not
relevant for constraining $\sigma_{8}$.

Figure 5: The relative constraints of $\sigma_{8}$ computed in linear theory,
with constant linear biases and number densities.
$\tilde{N}=(b^{2}\bar{n})^{-1}$ is the effective shot noise. $a$ and $b$ refer
to the one-tracer constraints from each tracer, $aa+bb$ refer to the auto-
spectra (incomplete multitracer) constraints, and $a\times b$ refer to the
complete multitracer constraint.  Figure 6: Left: Showing improvement in
relative constraints of $\sigma_{8}$, when including the high bias tracer to
the low bias (auto/$aa$) and including cross-spectrum to auto-spectra
(full/auto). Right: Showing improvement in relative constraints of
$b_{a}/b_{b}$ when including cross-spectrum to auto-spectra (full/auto). The
significant improvement in $\delta(b_{a}/b_{b})$ compared to
$\delta\log\sigma_{8}$ show the power of sample variance cancellation where
present.

###  4.2 Nonlinear

Now let us study the nonlinear scenario, with the fully general power spectrum
(Eq. 3.1). The existence of many more parameters allows for the existence of
significant parameter degeneracies, and hence the potential for significant
degeneracy breaking from multitracer techniques. We test on samples of
identical linear biases and noise to the linear regime examples, but with
further analysis including CMB lensing as well.

A subset of the results are shown in Figure 7. For some samples there is
visible improvement of constraints from adding the $ab$ cross-spectrum,
contrary to the marginal improvement seen in the linear regime. This indicates
that we have degeneracy breaking at hand, as the test case in the linear
regime has shown that the increase in information from the additional spectrum
only improves our constraints marginally. The average improvement in
constraining power by adding the cross-spectrum (auto$\to$full), which can be
interpreted as the average degeneracy breaking, is 17%. This marks a 22%
improvement from the low-bias tracer $a$ and a 68% improvement from the high-
bias tracer $b$. This effect of degeneracy breaking can vary significantly
with the noise regime that we are working in. For instance, at $\tilde{N}=10$
the improvement from auto$\to$full is 30% on average, while at $\tilde{N}=500$
it is only 5%.

One can hope to further improve these results by reducing the number of added
free parameters. A common practice in the field is to neglect the cross-
spectrum stochastic terms $N_{2n,ab}$, by assuming that the two galaxy tracers
have completely disjoint inhabitance of halos. While past work has shown that
the inclusion of these parameters do not crucially change the forecast
behavior [61], it will be worthwhile to investigate this further by
considering small-scale cross-clustering of LBGs and LAEs at the redshifts and
luminosities of relevance to future cosmology surveys.

From the perspective of added parameter space, the lack of new parameters
other than $\alpha_{0,\kappa}$ makes CMB lensing an ideal candidate for cross-
analysis. Our results, also shown in Figure 7, indicate that we have 11%
improvement on average when adding lensing to our full multitracer case. When
adding to each one-tracer case, we see a 6% improvement and 19% improvement
for tracers $a$ and $b$, respectively. Although the increments are
significantly less than that of the galaxy-galaxy multitracer analysis, it is
worthwhile to note that there is still improvement from lensing after
including two galaxy tracers. We may attribute the little improvement from
lensing to the limitations placed on the number of modes measured in a 2D
power-spectrum, which is significantly less than that of the 3D spectrum, and
the lack of cosmological parameters in this “fixed-shape” approach.

Figure 7: The relative constraints of $\sigma_{8}$ in the nonlinear regime,
with various noise and $b_{a}$ values. The constraints with smaller caps show
the respective constraints including lensing.

##  5 Radiative transfer

Radiative transfer (RT) is a known, yet unresolved potential cause of concern
for the use of LAEs as matter density tracers. Ref. [62] suggests that RT can
significantly vary the number of galaxies in catalog depending on the local
density and LOS velocity divergence. The latter effect is degenerate with the
$\mu$ dependence used to constrain $f\sigma_{8}$ in RSD studies (refer to §4
for a brief discussion in linear theory). A more recent work [63] argues that
the limited resolution of early simulations are the source of such a large
effect. Simulations with higher resolution show that regions of Ly$\alpha$
emission tend to have higher density, such that the photons escaping the
galaxy have shifted further away from resonance and thus are less susceptible
to scattering by the surrounding environment. While direct observations of RT
effect are currently unavailable, the impact of Hi absorption on LAE
clustering has been detected [64].

For an analysis using only LAEs, a large RT contribution to the selection of
galaxies would lead to very large errors on $\sigma_{8}$. However, with a
complex model that can potentially break degeneracies through multi-tracers,
these samples may prove to be useful. In order to make a preliminary
investigation of whether RT-affected tracers can still be valuable in multi-
tracers, we take the following approach. We imagine that the RT effect is
small enough that it can be included by taking

$P(k,\mu)\to P(k,\mu)+2\delta
f\left(b\mu^{2}+f\mu^{4}\right)P_{L}+\mathcal{O}(\delta f^{2})$ (5.1)

while leaving the higher loop terms unchanged. In principle there are
corrections to all of the higher order biases and the 1-loop terms as well,
but we are assuming that these terms are sub-dominant and multiplied by a
“small” RT effect and thus they can be neglected. A more complete and
consistent calculation will be attempted in a future publication.

##  6 Results

In this section we outline our main forecast results for two experiment
designs: ‘Stage 4.5’, motivated by DESI II, and ‘Stage V’. As discussed in §2
we will use one luminosity limit ($f_{\rm lim}=10\times 10^{-17}$ for LAEs and
$m=24.5$ for LBGs) as the fiducial design for Stage 4.5 and consider multiple
scenarios of Stage V luminosity limits for a preliminary optimization
constrained on the total number of galaxies observed. The experimental designs
are shown in Table 2.

In this section. we will consider gains from any additional tracer as “multi-
tracer improvement”, e.g. improvements from LBG-only analysis to an
LAE$\times$LBG or $\kappa\times$LBG analysis. This is unlike §4 where the
focus was on the improvement gained through the inclusion of cross-spectra,
and hence now the appropriate cross-spectra are always included when multiple
tracers are in use. It should not be a surprise, however, that an inclusion of
additional tracer improves the constraints with all else unchanged. To
distinguish the gains from multitracer from the simple increment in data, we
will be adjusting the number density of the LBG data available when including
LAEs. The constraint here will be the total number of galaxies observed, as we
have done when deciding the optimization for Stage V. As we have shown in
Figure 4 that LBGs are less shot-noise dominated with existing luminosity cuts
and the number density of LAEs are several factors smaller than LBGs, we will
observe LAEs at maximum efficiency and decrease the LBG observation efficiency
necessarily, to 0.61 (note that these efficiencies are multiplied by the 50%
success rate for redshift detection that is always in place)222We have
confirmed that the lack of inclusion of this LBG efficiency decrease affects
the improvement results by $<10\%$. A similar treatment cannot be made for
lensing as they are provided by CMB surveys separate from LSS experiments.
Hence, we will compare the galaxy-only and $\kappa\times$galaxy with the same
galaxy samples and for this reason the LAE and lensing multitracer
improvements should not be regarded on the same basis.

For all forecasts using full shape information, we include all relevant
$\Lambda$CDM and nuisance parameters (Table 3) and treat extra parameters of
interest as one parameter extensions. For cosmological priors, we will use
conservative BBN priors of $\sigma(\omega_{b})=0.0005$ or CMB primary priors
from Planck. A summary of the multi-tracers improvement in the Stage 4.5
design are provided in Figure 8.

The blue curve in the figure shows a case for when LAE number densities are
doubled, in the BBN prior scenario. This is to address the concern raised in
§2.1 that offsets of factor 2 on LAE number density estimates are not
unlikely. While this increases the improvement $\sim 10\%$ for most
parameters, there are no qualitative changes to the final results.

Figure 8: A summary of parameter constraint improvements from using
multitracers in the Stage 4.5 design. The ‘ST’ and ‘MT’ on the y-axis refer to
single tracer and multitracers, where the LBG-only sample is considered the
base single tracer sample. The type of marker depicts what tracer was added on
to the sample. The shown improvements for $f_{\rm EDE}$, $A_{\rm lin}$,
$A_{\rm log}$ are the average improvements across the parameter space probed.

###  6.1 Growth rate

In this section, we measure structure growth by measuring $\sigma_{8}$
constraints. While $\sigma_{8}$ is a well-motivated variable to pursue in
multi-tracers both in terms of existing tensions in the field and linear
theory arguments provided in §4, it is also a variable where we expect
significant impact from RT effects as mentioned in §5. It is possible that
multi-tracers break the RT-induced degeneracies enough such that inclusion of
LAEs significantly improve $\sigma_{8}$ constraints.

We use two approaches in this section. First we have our full-shape approach,
which marginalizes over all of the $\Lambda$CDM parameters. Another approach
is to use the ‘fixed-shape’ approach of §4, which takes a model-agonistic
approach and only marginalizes over amplitude. The latter approach
marginalizes over a subset of the former approach’s parameters and thus always
yields tighter constraints. In both approaches, the lensing auto-spectrum is
not included. This ensures that the $\sigma_{8}(z)$ constraints are local,
using only the data collected in each redshift bin up to the shared BBN prior.

The results for Stage 4.5 are shown in Figure 9. We find that the both
approaches gain $\approx 10\%$ from the introduction of LAEs. This is less
than the improvements from including lensing, which are $\approx 20-30\%$.
These improvements in total amount to $\approx 25-35\%$. Notice that unlike in
§4, the low bias tracer LAE do not have tighter $\sigma_{8}$ constraints than
the high bias tracer LBG. This is due to LAEs being more shot noise dominated
than LBGs at the luminosity limits of Stage 4.5. This increases the importance
of lensing as well. Multi-tracer benefits from LAEs diminishes significantly
when RT is present with the constraints weakening by approximately 5%. This
turns to an improvement when lensing is included
($\kappa\times$LBG$\to\kappa\times$LAE$\times$LBG), but the improvements
remain less than 5%.

Figure 9: Relative constraints on $\sigma_{8}(z)$ from Stage 4.5 in the
$\Lambda$CDM-derived (left) and fixed shape (right) methods. The smaller cap
sizes correspond to the constraints including lensing. We see that the
inclusion of RT effects diminish majority of the improvement from including
LAEs.

Figure 10 shows the Stage V constraints for all 3 scenarios, with both galaxy
tracers included without RT. While there are sizable differences between
scenarios, they are not drastic, especially with the inclusion of lensing.
This indicates that $\sigma_{8}$ constraints are robust to experiment designs.

Figure 10: Relative constraints on $\sigma_{8}(z)$ from Stage V surveys in the
$\Lambda$CDM-derived (left) and fixed shape (right) methods, with both LAEs
and LBGs, without RT. The smaller cap sizes correspond to the constraints
including CMB lensing. All error bars assume inclusion of all multi-tracer
spectra.

###  6.2 Standard model

Multi-tracers will further assist next-generation LSS experiments to provide
constraints similar to CMB experiments alone. Figure 11 shows a summary of
forecasts. It does not include $\omega_{b}$ as it shows no influence from
either multi-tracers nor RT, likely due to the prior dominating the
constraining power. Forecasts indicate that with a BBN prior, multi-tracers
through LAEs will improve the constraints on all other cosmological parameters
by $\lesssim 10\%$. Inclusion of lensing generally yields larger improvements,
with the improvement varying from $\sim 20\%$ for $h$ to $\sim 85\%$ for
$\Omega_{k}$. Overall, the improvements amount to $25-45\%$, with the
exception of $\Omega_{k}$ at 89%.

With Planck primary priors, there are varying multi-tracer effects. Multi-
tracers with LAEs on a LBG-only sample will have even less improvement
compared to when using BBN priors, at $\sim 5\%$, but when lensing is included
(i.e. LBG$\times\kappa\to\text{LAE}\times\text{LBG}\times\kappa$) the
improvements increase, marking $\sim 35\%$ for $h$ and $\omega_{c}$.
Similarly, the improvements from adding lensing vary across parameters,
marking a wide range between $\sim 20\%$ to $\sim 65\%$. $h$ and $\omega_{c}$
are again most affected, with constraints improving $\sim 60\%$ when adding
lensing to the LAE$\times$LBG sample. Refer to Figure 8 for specific parameter
constraint improvements.

Inclusion of RT generally has little effect on the cosmological constraints,
except for $\log(A_{s})$ that experience significant effects for forecasts
without $\kappa$ due to reasons similar to $\sigma_{8}$. The final constraints
with inclusion of all tracers experience sub-percent effect from RT, however.
Our results also indicate that Stage V designs provide constraints
significantly stronger than that of Stage 4.5, regardless of the luminosity
limits adopted.

Figure 11: $\Lambda$CDM and extension results. The ‘4.5+BBN’, ‘4.5+Planck’ and
‘SV+BBN’ each refer to Stage 4.5 with a BBN prior, Stage 4.5 with a Planck
prior, and Stage V with a BBN prior, respectively. The massive neutrino
constraints include DESI BAO and lensing auto-spectra for all results. The
magenta errorbars at $z_{*}$ are the constraints given by Planck [65].

###  6.3 Distances

One of the most important constraints from LSS surveys are the distance
measurements made through the BAO peaks in the power-spectra. Previous work
have forecasted that near-term surveys such as DESI, Euclid and HIRAX can
measure angular diameter distance and Hubble parameter to sub-percent and
percent precision, respectively, up to $z\sim 2$, and future experiments to
sub-percent precision up to $z\sim 5$ [5]. We will follow the modelling
methods used in ref. [5], marginalizing over linear bias and 15 polynomial
coefficients defined as

$P_{\rm obs}=\alpha_{\parallel}^{-1}\alpha_{\perp}^{-2}P_{\rm
recon}+\sum_{n=0}^{4}\sum_{m=0}^{2}c_{nm}k^{n}\mu^{2m}$ (6.1)

where $\alpha_{\parallel}$ and $\alpha_{\perp}$ are A-P parameters and $P_{\rm
recon}$ is the post-reconstruction power-spectrum calculated within the
Zeldovich approximation using the “RecSym” convention [66, 67]. We refer the
reader to ref. [5] for further details about the model.

Our results, shown in Figure 12, show that Stage 4.5 will reach sub-percent
precision in both distance measures in its entire redshift range $2.4<z<3.2$,
with the aid of multi-tracers improving constraints by $\sim 10-15\%$. More
significantly, these improvements are uninfluenced by the presence of RT. This
is not surprising, given that BAOs are a well-localized feature with a known
shape and the modelling process already marginalizes over any k dependence for
each spectrum beyond $P_{\rm recon}$.

Figure 12: Stage 4.5 constraints on A-P parameters ($\alpha_{\perp}$,
$\alpha_{\parallel}$), which can be interpreted as relative constraints on BAO
parameters ($D_{A}(z)/r_{d}$, $r_{d}H(z)$), respectively. The inclusion of RT
clearly have minimal effect for BAO constraints.

Forecasts for Stage V surveys, shown in Figure 13, indicate that future survey
of such kind can improve and further extend these constraints up to $z=4.4$.
The comparable results between the different designs of Stage V indicate that
such constraints are robust to the survey optimization, although fractional
improvements can be gained at some redshifts in price for constraining power
at other redshifts.

Figure 13: Stage V constraints for A-P parameters ($\alpha_{\perp}$,
$\alpha_{\parallel}$), which can be interpreted as relative constraints on
$D_{A}(z)/r_{d}$ and $r_{d}H(z)$, respectively. Each constraint assumes the
inclusion of all multi-tracer spectra.

###  6.4 Neutrino mass and light relics

Variation of neutrino mass is a reasonable extension to the $\Lambda$CDM
model. Measurements have yet to determine whether the neutrino mass hierarchy
is normal ($m_{\nu_{2}}<m_{\nu_{3}}$, $\Sigma m_{\nu_{i}}\gtrsim 58$meV) or
inverted ($m_{\nu_{2}}>m_{\nu_{3}}$, $\Sigma m_{\nu_{i}}\gtrsim 105$meV) [68].
Sum of neutrino masses can be constrained cosmologically by measuring its
effect on power-spectra. As neutrinos transition from relativistic to non-
relativistic species at scale factor $a_{\rm nr}$, neutrinos only cluster at
large scales ($k\ll k_{\rm nr}=a_{\rm nr}H(a_{\rm nr})$), and damp the power-
spectra on smaller scales [5].

For neutrino mass ($M_{\nu}=\Sigma m_{\nu_{i}}$) forecasts, we include DESI
BAO and SO lensing auto-spectrum for all results. The DESI BAO data included
is from the Bright Galaxy Survey (BGS; [69, 70]), Luminous Red Galaxies (LRG;
[71, 70]), and Emission Line Galaxies (ELG; [72]), post-reconstruction. We do
not include the Lyman-alpha data due to its redshift overlap with next-
generation facility galaxies. The results shown in Figure 11, show that Stage
4.5 will be able to constrain the total neutrino mass to $157$meV ($28$meV)
with BBN (Planck) priors. While improvements from individual tracers are
$\lesssim 10\%$ for BBN priors and $\lesssim 4\%$ for Planck priors, the total
improvement from both LAEs and lensing combined are double that, at 20% and
11% for BBN and Planck priors, respectively.

Ref. [73, 74, 5, 75] notes that improvement of $\tau$ measurements will
significantly tighten $M_{\nu}$ measurements by removing the $\tau$-$A_{s}$
degeneracy in primary CMB. Ref. [76] has shown that this leads to immense
advantages when considering the $\tau$-constraints provided by near-future CMB
surveys. Multiple galaxy tracers allow us to further take advantage of this.
Figure 14 shows the effect, by showing forecasts for $\sigma$s towards normal
hierarchy detection using external $\tau$ constraints from current and near-
future sources. In particular, LiteBIRD [77] will be able to provide cosmic
variance limited constraints on $\tau$ for CMB experiments at
$\sigma(\tau)=0.002$ [78, 79, 80]. At the tightest $\tau$ constraints from
21cm surveys [75, 73, 81, 82], the inclusion of LAEs will improve the
measurements by 0.7$\sigma$ and 0.3$\sigma$ for Stage 4.5 and Stage V,
respectively, allowing for the 5.5$\sigma$ detection of the normal mass
hierarchy.

Figure 14: Fidelity of detection for normal neutrino mass hierarchy given
$\tau$ constraints provided by external facilities. Each contour represents
different information used for forecast. The red contour is only the default
set of data, from Planck+SO CMB prior, DESI BAO measurements, and SO lensing
auto-spectrum. The other contours include Stage 4.5 of Stage V galaxies and
appropriate galaxy-lensing cross-correlations in addition to the default set.
The labels show external constraints on $\tau$ that can be obtained currently
or in the near-future. $\sigma(\tau)$ in each case are 0.0071 [83, 5, 84],
0.0058 [5, 75], 0.002 [79, 80], and 0.0009 [75], for Planck, Planck+SO+DESI
(PSD), LiteBIRD, and 21cm surveys, respectively. In particular, the constraint
provided by LiteBIRD is the cosmic variance limit [78] for CMB experiments.

The light relic forecasts on $N_{\rm eff}$ show that while near future
constraints from Stage 4.5 and Stage V will be comparable to or tighter than
the Planck constraints. The assistance of multi-tracers are more significant
with a BBN prior, with the total improvement at $38\%$. With a Planck prior,
the help of both tracers decrease, marking a total improvement of 25%.
Finally, we see that the choice of Stage V design has minimal effect on the
final constraints, providing errors tighter than that of Planck using only
galaxy tracers and a BBN prior.

###  6.5 Primordial features

Current constraints of inflation, primarily by CMB observations, are
consistent with a simple, single scalar field slow-roll inflation. The (near)
flatness of the inflationary potential leads to a simple explanation of the
nearly scale-invariant primordial power-spectrum. When formulating a model
from a more fundamental basis, however, one often finds deviation from said
scale-invariance [5, 85, 86]. A canonical subclass of such models approximate
deviation from scale-invariance over a narrow range of scales, using
oscillations in $k$ or $\log{k}$. For this work, we model this as oscillations
in the linear power-spectrum

$P_{m}(k)\rightarrow[1+A_{\rm lin}\sin{(\omega_{\rm lin}k+\phi_{\rm
lin})}+A_{\rm log}\sin{(\omega_{\rm log}\log{k/k_{*}}+\phi_{\rm
log})}]P_{m}(k)$ (6.2)

where $k_{*}=0.05\,h^{-1}\,\text{Mpc}$ is the pivot scale [86, 87, 5]. We fix
$\phi_{\rm lin}$ and $\phi_{\rm log}$ to $\pi/2$ for all forecasts and
characterize the change in constraining power as we vary $\omega_{\rm lin}$
and $\omega_{\rm log}$.

Our results, shown in Figures 15 indicate that use of multiple galaxy tracers
benefits our sensitivity by 8% for $A_{\rm lin}$ and 16% for $A_{\rm log}$
with BBN priors, and less than 2% for both with Planck priors. Furthermore,
the inclusion of RT minimally affect the constraints, as increasing
uncertainty in $\mu$ dependence does not directly alter our ability to detect
$k$ dependence. We find that these constraints are robust against crude
optimization strategies for Stage V.

Figure 15: Forecast results for $A_{\rm lin}$ and $A_{\rm log}$ for a BBN-
prior. The loss in $\delta A_{\rm lin}$ near $\omega_{\rm lin}=100$ with the
BBN prior is due to the degeneracy with BAO [5], but this degeneracy is
removed when using a Planck prior.

###  6.6 Early dark energy

While current measurements of cosmology are consistent with that of the
cosmological constant $\Lambda$ as the source of dark energy, dynamical dark
energy is an active area of research. Early dark energy (EDE) is among these
models, with $\Omega_{\rm DE}$ rapidly dropping at high redshift. The specific
model that we probe has the potential $V(\phi)\propto[1-\cos{(\phi/f})]^{3}$,
where $\phi$ is a light scalar field [88, 5]. The field is parameterized by
redshift $z_{c}$, where $f_{\rm EDE}\equiv\rho_{\rm EDE}/\rho_{\rm tot}$ is
maximum, and the initial value of the field $\theta_{i}\equiv\phi_{i}/f$. This
model, as well as other major EDE models, invokes a phase shift in BAO when
$z_{c}\gtrsim z_{*}$, where $z_{*}\sim 1100$ is recombination. Otherwise, when
$z_{c}\lesssim z_{*}$, EDE damps the power-spectrum equally across scales for
$k\gtrsim 0.01\,h\,\text{Mpc}^{-1}$ [5]. This can be understood by simply
considering the derivation for sound horizon $r_{s}$ and angular diameter
distance $D_{A}(z)$

$r_{s}=\int_{z_{*}}^{\infty}\frac{c_{s}(z)}{H(z)}dz\quad\quad
D_{A}(z)=\int_{0}^{z}\frac{dz}{H(z)}$ (6.3)

When $z_{c}\gtrsim z_{*}$, the BAO scale is impacted by influence on $r_{s}$,
while for $z_{c}\lesssim z_{*}$, there is no influence on either parameters
(as long as $z<z_{c}$, which is always the case for realistic models and any
quasi-near-time surveys). The phase shifts on BAO, or more specifically
modifications of the Hubble parameter measurements from the CMB, has the
potential to relieve the Hubble tension. For the forecasts in this section, we
fix $\theta_{i}=2.83$ and probe $1.5<\log_{10}(z_{c})<6.5$ using CLASS_EDE
[89].

The results shown in Figure 16 indicate that multi-tracers play a significant
role in EDE constraints. While the improvements through LAEs alone remain
modest at $\lesssim 10\%$ for both priors, lensing improves the constraints by
$91\%$ with a BBN prior and $\sim 64\%$ with a Planck prior, totalling
together to 97% and 85%, respectively.

We also observe that composite constraint is affected marginally by the
presence of RT, while for LAE-only constraints this is not always true. Ref.
[5] makes note that LSS-only constraints (without RT) are worsened at low
$z_{c}$ ($z_{c}\lesssim z_{*}$) and attribute this to the degeneracy between
the damping of power-spectra at $k\gtrsim 0.01\,h\,\text{Mpc}^{-1}$, $A_{s}$
and $b$. We recognize that for LAE-only constraints in the Planck prior case,
RT significantly diminishes our constraining power further at this regime.
This is entirely expected, as the variational $f$ at linear-theory level
introduces an additional degree of freedom on the existing amplitude
degeneracy. Notably, the inclusion of lensing nearly recovers our constraining
power lost from RT. This again suggests that the inclusion of lensing aids the
breaking of degeneracy introduced by RT.

Figure 16: Forecast results for early dark energy. RT induces a loss of
constraining power at with a Planck prior due to the additional amplitude
degeneracy.

##  7 Conclusions

Next generation experiments in large-scale structure aim to probe the high
redshift universe to constrain the $\Lambda$CDM model and its extensions.
Modelling of large-scale structure using perturbation theory introduces
nuisance parameters (Table 3) such as bias $b_{X}$, counter-terms
$\alpha_{2n}$, and stochastic terms $N_{2n}$ in order to make predictions in a
controlled manner. A major limitation to cosmological constraints are the
degeneracy that these parameters introduce in the data, as well as inherent
degeneracies between cosmological parameters. Tracing multiple different
tracers at the same redshift, called the “multi-tracer” technique, may allow
us to break these degeneracies by using the different dependence of tracers to
each parameter. Next generation facilities may choose to do this with Lyman
alpha emitters (LAEs), Lyman break galaxies (LBGs), and CMB lensing $\kappa$.
With different large-scale biases (Figures 1, 2, or the lack thereof), these
tracers will be well-suited for multi-tracer cosmology. In this work, we have
compared forecasted cosmological constraints from single tracers, an
“incomplete” multi-tracer technique (where one uses multiple galaxy tracers
but only includes the two galaxy auto-spectra as signals neglecting the
galaxy-galaxy cross-spectrum) and one which includes all available spectra.

We start with a generic investigation of the mechanism of multi-tracer
cosmology, using $\sigma_{8}$ forecasts for a set of simplified tracers. We
find that degeneracy breaking indeed drives a majority of the improvement in
constraining power, in agreement with recent studies [60, 61].

We then apply the multi-tracer technique for Fisher forecasts of next
generation surveys, namely “Stage 4.5”, inspired by DESI II, and “Stage V”,
using both BBN and Planck primary priors. We find that the multi-tracer
technique has varying effect on constraining power, as shown in Figure 8. In
general, multi-tracer with LAEs on a LBG-only sample do not improve the
constraints drastically with the improvements $\lesssim 15\%$ for most
parameters. Interestingly, when using LAEs as an addition onto a
$\kappa\times$LBG sample, there are parameters that improve considerably,
namely $h$ and $\omega_{c}$ with a Planck prior, at $\sim 35\%$, and $f_{\rm
EDE}$ at $\sim 60\%$. CMB lensing proves to be a better contribution to multi-
tracers than LAEs in many cases, although the improvements cannot be compared
on equal footing as they are adopted without any cost, while observing LAEs
cost observation time that can otherwise be used for LBGs. With lensing, many
of cosmological parameters being improved $\gtrsim 20\%$, with a factor of
$\gtrsim 10$ improvement for $f_{\rm EDE}$ being the most drastic.
Importantly, we find that these improvements aggregate well in general (e.g.,
for $\delta\sigma_{8}$ the $\sim 10\%$ and $\sim 20\%$ improvements from LAEs
and lensing, respectively, lead to $\sim 30\%$ improvements when combined),
proving that we are not in a regime where there are enough tracers that
additional tracers become insignificant.

Where they overlap our results agree well with the recent study of ref. [61],
although the results are presented differently. Ref. [61] quotes $\sim 50\%$
improvement in constraining power from multi-tracers with two galaxy samples,
using $\sigma_{ST}/\sigma_{MT}-1$ as their metric. We use
$1-\sigma_{MT}/\sigma_{ST}$ which leads to smaller quoted improvements, e.g. a
50% improvement in ref. [61] corresponds to a 33% improvement in our metric,
bridging a significant portion of the apparent gap between the final results.
Remaining differences in results can be attributed to differences in the
galaxy samples used, such as the difference between galaxy biases. A similar
comment can be made for differences between this work and ref. [90], which
show improvements ranging from 5% to 33% from ELGs and LRGs. In both studies,
there is a factor $\sim 2$ difference between linear biases, whereas LAE and
LBG biases differ only by 30% at most (for Stage 4.5) in this work.

We also perform a crude optimization of Stage V survey designs by changing
luminosity limits for both LAEs and LBGs while keeping the total integration
time fixed. We find that most Stage V constraints are robust to this
optimization, leaving room for future optimization involving specific targets
and observation techniques. A detailed study of survey optimization, including
galaxy mix as a function of redshift, time-to-redshift and redshift efficiency
factors, is left for future work.

A possible source of concern for the use of LAEs as cosmological probes is the
impact of radiative transfer (RT), specifically the absorption and scattering
of Ly$\alpha$ emission lines by neutral hydrogen near the LAE. By assuming
that RT effects are small, we attempt a preliminary investigation of the
impact of RT by including its effects at the linear theory level, ignoring
higher order contributions. Forecasts using this model show that $\sigma_{8}$
loses nearly half of its multi-tracer gains by LAEs, and LAE-only forecasts
experience significant loss of constraining power for $\sigma_{8}$, $A_{s}$,
neutrino mass, and early dark energy with Planck priors at low $z_{c}$ due to
the additional amplitude degeneracy. On the contrary, however, we find that
other forecasts are fairly unaffected, and all constraints experience little
effect for the final, three-tracer forecast. Of course, this treatment of RT
is merely preliminary and a more complete calculation will be attempted in the
future.

## Acknowledgments

We thank Arjun Dey, Shi-Fan Chen, Simone Ferraro, Noah Sailer, and Ruiyang
Zhao for useful discussions throughout the course of this work. H.E. and M.W.
are supported by the DOE. This research has made use of NASA’s Astrophysics
Data System and the arXiv preprint server. This research is supported by the
Director, Office of Science, Office of High Energy Physics of the U.S.
Department of Energy under Contract No. DE-AC02-05CH11231, and by the National
Energy Research Scientific Computing Center, a DOE Office of Science User
Facility under the same contract.

##  Appendix A Variance of measurements with correlation

In this section we walk through the calculation of the composite variance from
two correlated measurements, leading to Eq. 4.5. Suppose that two measurements
$(x_{1},x_{2})$ are made drawing from a Gaussian distribution of mean $\mu$
and covariance $C$. The covariance matrix will be of the form

$C=\begin{pmatrix}\sigma_{1}^{2}&r\sigma_{1}\sigma_{2}\\\
r\sigma_{1}\sigma_{2}&\sigma_{2}^{2}\end{pmatrix}$ (A.1)

where $r\in[-1,1]$ is the correlation coefficient. The probability of
measuring $(x_{1},x_{2})$, then, is

$P(x_{1},x_{2})\propto\exp{(\vec{x}-\vec{\mu})^{T}C^{-1}(\vec{x}-\vec{\mu})/2}$
(A.2)

where $\vec{x}=(x_{1},x_{2})$ and $\vec{\mu}=(\mu,\mu)$. Maximizing the
likelihood of observation provides the best estimate of $\mu$, with variance
$\sigma^{2}$

$\mu=\frac{\vec{1}^{T}C^{-1}\vec{x}}{\vec{1}^{T}C^{-1}\vec{1}},\quad\quad\sigma^{2}=\frac{(1-r^{2})\sigma_{1}^{2}\sigma_{2}^{2}}{\sigma_{1}^{2}-2r\sigma_{1}\sigma_{2}+\sigma_{2}^{2}}$
(A.3)

where $\vec{1}=(1,1)$. Taking some limits with these expressions reveal
important insights. In the fully uncorrelated case $r=0$, one recovers the
standard weighted average
$\mu=(x_{1}/\sigma_{1}^{2}+x_{2}/\sigma_{2}^{2})(1/\sigma_{1}^{2}+1/\sigma_{2}^{2})^{-1}$,
$\sigma^{2}=(1/\sigma_{1}^{2}+1/\sigma_{2}^{2})^{-1}$. When two measurements
are redundant, with $\sigma_{1}=\sigma_{2}$ and $r\to 1$, one obtains
$\sigma^{2}=\sigma_{1}^{2}$, as expected. The best case is where there is near
full correlation $\absolutevalue{r}\to$, with non-redundant information
$\sigma_{1}\neq\sigma_{2}$. Then, one obtains $\sigma^{2}\to 0$.

##  Appendix B Conventions of limiting luminosity in literature

Literature from the astronomy community often use two different conventions to
make note of their luminosity function integration range. In this section we
show the conversion between the two conventions.

The first approach is to simply use the lower limit of the integral, e.g.

$\bar{n}=\int_{L_{\rm lim}}^{\infty}\Phi(L^{\prime})dL^{\prime}$ (B.1)

Another is to use the median, e.g. $L_{\rm med}$, such that

$\int_{L_{\rm
med}}^{\infty}\Phi(L^{\prime})dL^{\prime}=\frac{1}{2}\int_{L_{\rm
lim}}^{\infty}\Phi(L^{\prime})dL^{\prime}$ (B.2)

Note that in some authors place an explicit, finite upper limit on the
integrals depending on the context. These definitions are extended for
limiting line fluxes or magnitudes analogously. It is clear that with the
luminosity function (LF) parameters, one can convert one to the other easily.
For the regimes of our interest, the differences are typically $\sim 0.1-0.2$
in log-space.

##  Appendix C Counterterms

The theory of Lagrangian EFT has been developed in refs. [91, 92, 93, 94, 95,
96, 97, 98, 99, 100, 101, 102, 49, 103, 104, 105, 106]. However in those works
the cross-spectra between different tracers was not highlighted. In this
appendix we briefly recap the key points of Lagrangian EFT, with a focus on
the counterterms for the cross-spectra which have not appeared explicitly
before. This will also allow us to make contact with the Eulerian approach of
ref. [61], which has the same counterterm structure.

We will perform calculations at both the density field-level and the power-
spectrum-level. The latter is a more precise follow-through of the
calculation, as the expectation values are taken explicitly, but we will see
that the former will yield the same result with a simpler interpretation.

We will start with the field-level interpretation. In Lagrangian EFT, the
density field of a tracer is given by

$\delta_{a}+(2\pi)^{3}\delta_{D}(\textbf{k})=\int_{\textbf{q}}F[\delta(q)]\,e^{i\textbf{k}\cdot(\textbf{q}+\Psi(\textbf{q}))}$
(C.1)

where $\Psi$ is the displacement field and the functional $F$ describes the
tracer’s dependence on underlying fields using bias. For this section, we
limit ourselves to linear order, $F=1+b_{1}\delta$, with $b_{1}$ the
Lagrangian bias. Including higher-order terms ($\delta^{2}$, $s^{2}$, etc.)
will not add to this section as their contributions to the counterterm will
exceed one-loop order. For similar reasons we will limit ourselves to the
linear solution of $\Psi$.

To obtain the counterterms, we begin by dividing the displacements into long-
wavelength components and a short-wavelength component. The former are treated
perturbatively, while the latter can only be constrained by symmetry arguments
and are included as a derivative expansion. Some of these terms correlate with
the long-wavelength contributions while some do not. We refer to the former as
counterterms and the latter as stochastic terms. These terms allow for the
final results to be independent of the cutoff scale $\Lambda$ used in defining
the perturbative terms.

To linear order the (long-wavelength) displacement field $\Psi^{L}$ is simply
$(i\vec{k}/k^{2})\,\delta_{\rm lin}(k)$ [91]. In redshift space we apply an
additional displacement along the line of sight [107]. Since
$\dot{\Psi}^{L}=f\Psi^{L}$, in Hubble units, for the linear displacement the
transition to redshift space is accomplished by multiplying the line-of-sight
component of $\Psi^{L}$ by $1+f$. Then, in the absence of counterterms (i.e.
absence of small-scale contributions), we obtain the Kaiser power spectrum by
expanding the exponent $e^{i\textbf{k}\cdot\Psi}$ as

$e^{i\textbf{k}\cdot\Psi(s)}\simeq
1+i\textbf{k}\cdot\Psi(s)=1+i\textbf{k}\cdot[\Psi(q)+f\hat{n}(\hat{n}\cdot\Psi(q))]$
(C.2)

where $q$ and $s$ are the real-space and redshift-space coordinates,
respectively. Substituting this to Eq. C.1 and using
$\mathbf{k}\cdot\hat{n}=k\mu$ yields
$\delta_{a}=([1+b_{1}^{(a)}]+f\mu^{2})\delta_{\rm lin}$, where $1+b_{1}$ is
equivalent to linear bias in Eulerian PT.

Now let us divide our displacement field into
$\Psi=\Psi^{L}+c_{0}(t)\nabla\delta_{\rm lin}+\epsilon_{i}+\cdots$. The term
$c_{0}\nabla\delta_{\rm lin}$ is the first-order counterterm (i.e. due to
small scale dynamics) contribution to $\Psi$ allowed by symmetry, whereas
$\epsilon$ is the stochastic noise. As the time-dependence of $c_{0}$ and
$\epsilon$ are unknown, we gain additional degrees of freedom for
$\dot{\Psi}$. Let us define this as $c_{2}$ and $\dot{\epsilon}$ such that

$\displaystyle\Psi_{i}$ $\displaystyle=\Psi_{i}^{L}+c_{0}\nabla\delta_{\rm
lin}+\epsilon_{i}$ (C.3) $\displaystyle\dot{\Psi}_{i}$
$\displaystyle=\dot{\Psi}_{i}^{L}+c_{2}\nabla\delta_{\rm
lin}+\dot{\epsilon}_{i}$ (C.4)

Performing a direct expansion for small-scale contributions as we did for
linear theory calculations, we find

$\displaystyle\delta_{a}$
$\displaystyle\supset\int_{\textbf{q}}e^{i\textbf{k}\cdot\textbf{q}}\sum_{n}\frac{1}{n!}[ik_{i}(\Psi_{i}+\dot{\Psi}_{\hat{n},i})]^{n}$
(C.5)
$\displaystyle\supset\int_{\textbf{q}}e^{i\textbf{k}\cdot\textbf{q}}\left(-c_{a,0}-c_{a,2}\mu^{2}+\frac{1}{2}(1+f\mu^{2})((\hat{k}_{i}\epsilon_{i})^{2}+(\hat{n}_{i}\dot{\epsilon}_{i})^{2}\mu^{2})\right)k^{2}\delta_{\rm
lin}$ (C.6)

where we use
$\dot{\Psi}^{L}_{\hat{n},i}=\hat{n}_{i}\hat{n}_{j}\dot{\Psi}^{L}_{j}$ for
shorthand and the subscripts $a$ on coefficients allow for tracer-dependent
small-scale dynamics (whereas stochastic contributions $\epsilon$ do not have
tracer-dependence). Note that we have neglected terms that go as even powers
of $\delta_{\rm lin}$ or odd powers of $\epsilon,\,\dot{\epsilon}$ as they
vanish when taking expectation values at the power-spectrum level.

To lowest order the counterterms for the autospectrum of tracer $a$ are thus
$k^{2}\,P_{\rm lin}$ times

$2\left([1+b_{1}^{(a)}]+f\mu^{2}\right)\left(-\tilde{c}_{a,0}+-\tilde{c}_{a,2}\mu^{2}-\tilde{c}_{4}\mu^{4}\right)$
(C.7)

while those for the cross-spectrum between tracers $a$ and $b$ are
$k^{2}\,P_{\rm lin}$ times

$\left([1+b_{1}^{(a)}]+f\mu^{2}\right)\left(-\tilde{c}_{b,0}+-\tilde{c}_{b,2}\mu^{2}-\tilde{c}_{4}\mu^{4}\right)+\left([1+b_{1}^{(b)}]+f\mu^{2}\right)\left(-\tilde{c}_{a,0}+-\tilde{c}_{a,2}\mu^{2}-\tilde{c}_{4}\mu^{4}\right)$
(C.8)

where we have redefined coefficients for simplicity. The cross-spectra thus
introduce no new free parameters beyond those already needed for the
counterterms of the auto-spectra, and this increases the utility of the cross-
spectrum in breaking degeneracies. Note that although $\tilde{c}_{4}$ is a
shared contribution across tracers, in this work we will allow them to be
tracer-dependent, in order to be conservative. We expect that this will change
the results minimally.

Now, we will perform the calculation at the power-spectrum-level to see that
our results above are indeed correct. To do this, we will follow the moment
expansion approach [106, 108, 109]. In this approach, one uses the fact that
the redshift space power-spectrum is a specific case of the moment generating
function and expands it such that

$\displaystyle\frac{k^{3}}{2\pi^{2}}P_{s}(\textbf{k})=\tilde{M}(\textbf{J}=\textbf{k},\textbf{k})$
$\displaystyle=\frac{k^{3}}{2\pi^{2}}\int
d^{3}re^{i\textbf{k}\cdot\textbf{r}}\expectationvalue{(1+\delta_{g}(\textbf{x}_{1}))(1+\delta_{g}(\textbf{x}_{2}))e^{i\textbf{k}\cdot\Delta\textbf{u}}}_{\mathbf{x}_{1}-\mathbf{x}_{2}=\mathbf{r}}$
(C.9)
$\displaystyle=\frac{k^{3}}{2\pi^{2}}\sum_{n=0}^{\infty}\frac{i^{n}}{n!}k_{i_{1}}\dots
k_{i_{n}}\tilde{\Xi}_{i_{1}\dots i_{n}}^{(n)}(\textbf{k})$ (C.10)

where $\mathbf{u}=\hat{n}(\hat{n}\cdot\mathbf{v})/\mathcal{H}$ is the LOS
velocity,
$\Delta\mathbf{u}=\mathbf{u}(\mathbf{x}_{1})-\mathbf{u}(\mathbf{x}_{2})$ their
difference, and $\tilde{\Xi}^{(n)}$ are Fourier transforms of the velocity
moments

$\Xi_{i_{1}\dots
i_{n}}^{(n)}=\expectationvalue{(1+\delta_{1})(1+\delta_{2})\Delta\textbf{u}_{i_{1}}\dots\Delta\textbf{u}_{i_{n}}}\quad.$
(C.11)

Notice that here, the first moment is simply the real-space power-spectrum and
the higher moment contributions provide information about redshift space
distortions (RSD). As RSD depend only on the line-of-sight component of
velocities, the only non-vanishing contributions to $P_{s}$ come from the LOS
component of $\tilde{\Xi}$. With this in mind, we rewrite the equation above
as

$P_{s}(\textbf{k})=\sum_{n=0}^{\infty}\frac{i^{n}}{n!}(k\mu)^{n}\sum_{\ell=0}^{n}\tilde{\Xi}_{\ell}^{(n)}\mathcal{L}_{\ell}(\mu)$
(C.12)

where $\tilde{\Xi}_{\ell}^{(n)}$ are Legendre moments of the velocity moments.

Formally, counterterm contributions arise from UV-divergences that arise from
contractions of operators at small distances. As our theory knows only of
large-scale dynamics, the only constraints we have on these contributions are
from symmetry. Those that appear in the final expression are

$\displaystyle\expectationvalue{\Psi_{i}\Psi_{j}}$
$\displaystyle\supset\beta_{1}\delta_{ij}+\beta_{2}\delta_{ij}\delta+\beta_{3}\hat{k}_{i}\hat{k}_{j}\delta$
(C.13) $\displaystyle\expectationvalue{\delta\Psi_{i}}$
$\displaystyle\supset\beta_{4}ik_{i}\delta$ (C.14)

and similarly for contact terms with derivatives, e.g.
$\Psi_{i}\dot{\Psi}_{j}\supset\dot{\beta}_{1}\delta_{ij}+\dot{\beta}_{2}\delta_{ij}\delta+\dot{\beta}_{3}\hat{k}_{i}\hat{k}_{j}\delta$.

Using the above contributions, one can compute the contributions to each
velocity moment multipole. For example, we have

$\displaystyle\dddot{W}_{ijk}$
$\displaystyle=\expectationvalue{\dot{\Delta}_{i}\dot{\Delta}_{j}\dot{\Delta}_{k}}$
(C.15)
$\displaystyle=\underbrace{\expectationvalue{\dot{\Psi}_{i}(q_{1})\dot{\Psi}_{j}(q_{1})\dot{\Psi}_{k}(q_{1})-\left(q_{1}\to
q_{2}\right)}}_{\propto\expectationvalue{\delta}=0\,\text{for
c.t.}}+\expectationvalue{\dot{\Psi}_{\\{i}(q_{1})\dot{\Psi}_{j}(q_{2})\dot{\Psi}_{k\\}}(q_{2})-(q_{1}\leftrightarrow
q_{2}))}$ (C.16)
$\displaystyle\supset\expectationvalue{\frac{2ik_{\\{i}}{k^{2}}\delta_{1}\left(\ddot{\beta}_{1}\delta_{jk\\}}+\ddot{\beta}_{2}\delta_{jk\\}}\delta_{2}+\ddot{\beta}_{3}\hat{k}_{j}\hat{k}_{k\\}}\delta_{2}\right)}$
(C.17)
$\displaystyle=\frac{2i}{k}P_{L}f\left(\ddot{\beta}_{2}\hat{k}_{\\{i}\delta_{jk\\}}+3\ddot{\beta}_{3}\hat{k}_{i}\hat{k}_{j}\hat{k}_{k}\right)$
(C.18)

Using the velocity moment expressions [45], we see that this contributes to
the velocity moment multipoles

$-i\tilde{\Xi}^{(3)}_{1,ct}=\frac{6}{k}P_{L}f\left(\ddot{\beta}_{2}+\frac{3}{5}\ddot{\beta}_{3}\right)\quad,\quad-i\tilde{\Xi}^{(3)}_{3,ct}=\frac{12}{5k}P_{L}f\ddot{\beta}_{3}$
(C.19)

which in turn contributes

$P_{s.ct}\supset
k^{2}P_{L}(\ddot{\beta}_{2}f\mu^{4}+\ddot{\beta}_{3}f\mu^{6})$ (C.20)

Note that for $n=3$ we use the Legendre basis instead of the polynomial basis
shown in the Appendix of [45]. Performing a similar computation for all terms,
one obtains

$P_{s,ct}=k^{2}P_{L}(\beta_{2}+\beta_{3}-2b\beta_{4}+(2\dot{\beta}_{2}+2\dot{\beta}_{3}-2b\dot{\beta}_{4}+\ddot{\beta}_{2})\mu^{2}+\ddot{\beta}_{3}\mu^{4})(1+b+f\mu^{2})$
(C.21)

in agreement with the field-level calculation above.

## References

  * [1] J. A. Peacock, Cosmological Physics.  Jan., 1999.
  * [2] S. Dodelson, Modern cosmology.  2003\.
  * [3] S. Dodelson and F. Schmidt, Modern Cosmology.  Elsevier Science, 2020.
  * [4] D. Baumann, Cosmology.  Cambridge University Press, 2022.
  * [5] N. Sailer, E. Castorina, S. Ferraro, and M. White, Cosmology at high redshift — a probe of fundamental physics, JCAP 12 (2021), no. 12 049, [arXiv:2106.09713].
  * [6] R. Ellis and K. Dawson, SpecTel: A 10-12 meter class Spectroscopic Survey Telescope, in Bull. Am. Astron. Soc. , vol. 51, p. 45, Sept., 2019.  arXiv:1907.06797.
  * [7] K. Dawson, J. Frieman, K. Heitmann, et al., Cosmic Visions Dark Energy: Small Projects Portfolio, arXiv:1802.07216.
  * [8] A. Blanchard et al., Gravitation And the Universe from large Scale-Structures: The GAUSS mission concept. Mapping the cosmic web up to the reionization era, arXiv:2102.03931.
  * [9] D. J. Schlegel, S. Ferraro, G. Aldering, et al., A Spectroscopic Road Map for Cosmic Frontier: DESI, DESI-II, Stage-5, arXiv e-prints (Sept., 2022) arXiv:2209.03585, [arXiv:2209.03585].
  * [10] S. Ferraro, N. Sailer, A. Slosar, and M. White, Snowmass2021 Cosmic Frontier White Paper: Cosmology and Fundamental Physics from the three-dimensional Large Scale Structure, arXiv e-prints (Mar., 2022) arXiv:2203.07506, [arXiv:2203.07506].
  * [11] V. Desjacques, D. Jeong, and F. Schmidt, Large-scale galaxy bias, Phys. Rep. 733 (Feb., 2018) 1–193, [arXiv:1611.09787].
  * [12] P. McDonald and U. Seljak, How to evade the sample variance limit on measurements of redshift-space distortions, Journal of Cosmology and Astro-Particle Physics 2009 (Oct., 2009) 007, [arXiv:0810.0323].
  * [13] A. Lewis and A. Challinor, Weak gravitational lensing of the CMB, Phys. Rep. 429 (June, 2006) 1–65, [astro-ph/0601594].
  * [14] D. Hanson, A. Challinor, and A. Lewis, Weak lensing of the CMB, General Relativity and Gravitation 42 (Sept., 2010) 2197–2218, [arXiv:0911.0612].
  * [15] M. Giavalisco, Lyman-Break Galaxies, Ann. Rev. Astron. & Astrophys.  40 (2002) 579–641.
  * [16] A. E. Shapley, Physical Properties of Galaxies from z = 2-4, Ann. Rev. Astron. & Astrophys.  49 (Sept., 2011) 525–580, [arXiv:1107.5060].
  * [17] M. J. Wilson and M. White, Cosmology with dropout selection: straw-man surveys & CMB lensing, Journal of Cosmology and Astro-Particle Physics 2019 (Oct., 2019) 015, [arXiv:1904.13378].
  * [18] M. Ouchi, Y. Ono, and T. Shibuya, Observations of the Lyman-$\alpha$ Universe, Ann. Rev. Astron. & Astrophys.  58 (Aug., 2020) 617–659, [arXiv:2012.07960].
  * [19] M. Dijkstra, Ly$\alpha$ emitting galaxies as a probe of reionisation, Publications of the Astronomical Society of Australia 31 (2014) e040.
  * [20] T. Garel, J. Blaizot, B. Guiderdoni, et al., The UV, Lyman $\alpha$, and dark matter halo properties of high-redshift galaxies, Mon. Not. R. Astron. Soc.  450 (June, 2015) 1279–1294, [arXiv:1503.06635].
  * [21] A. A. Khostovan, D. Sobral, B. Mobasher, et al., The clustering of typical ly$\alpha$ emitters from z $\sim$ 2.5–6: host halo masses depend on ly$\alpha$ and UV luminosities, Monthly Notices of the Royal Astronomical Society 489 (aug, 2019) 555–573.
  * [22] D. Sobral, S. Santos, J. Matthee, et al., Slicing COSMOS with SC4k: the evolution of typical ly$\alpha$ emitters and the ly$\alpha$ escape fraction from z $\sim$ 2 to 6, Monthly Notices of the Royal Astronomical Society 476 (feb, 2018) 4725–4752.
  * [23] A. Konno, M. Ouchi, K. Nakajima, et al., Bright and faint ends of ly$\alpha$ luminosity functions at z = 2 determined by the subaru survey: Implications for agns, magnification bias, and ism h i evolution, The Astrophysical Journal 823 (may, 2016) 20.
  * [24] M. Ouchi, K. Shimasaku, M. Akiyama, et al., The Subaru/XMM-Newton Deep Survey (SXDS). IV. Evolution of Ly$\alpha$ Emitters from z = 3.1 to 5.7 in the 1 deg2 Field: Luminosity Functions and AGN, Astrophys. J. Suppl.  176 (June, 2008) 301–330, [arXiv:0707.3161].
  * [25] L. Guaita, E. Gawiser, N. Padilla, et al., Ly$\alpha$-emitting galaxies at z = 2.1 in ecdf-s: Building blocks of typical present-day galaxies?*, The Astrophysical Journal 714 (apr, 2010) 255.
  * [26] R. Bielby, P. Tummuangpak, T. Shanks, et al., The vlt lbg redshift survey \- v. characterising the z = 3.1 lyman alpha emitter population, Monthly Notices of the Royal Astronomical Society 456 (01, 2015).
  * [27] Y. Shioya, Y. Taniguchi, S. S. Sasaki, et al., Photometric properties of ly$\alpha$ emitters at z $\approx$ 4.86 in the cosmos 2 square degree field*, The Astrophysical Journal 696 (apr, 2009) 546.
  * [28] C.-N. Hao, J.-S. Huang, X. Xia, et al., A deep ly$\alpha$ survey in ecdf-s and cosmos. i. general properties of ly$\alpha$ emitters at z $\sim$ 2, The Astrophysical Journal 864 (sep, 2018) 145.
  * [29] E. Gawiser, H. Francke, K. Lai, et al., Ly$\alpha$-emitting galaxies at z = 3.1: L* progenitors experiencing rapid star formation*, The Astrophysical Journal 671 (dec, 2007) 278.
  * [30] M. Ouchi, K. Shimasaku, H. Furusawa, et al., Statistics of 207 ly$\alpha$ emitters at a redshift near 7: Constraints on reionization and galaxy formation models*, The Astrophysical Journal 723 (oct, 2010) 869.
  * [31] P. Schechter, An analytic expression for the luminosity function for galaxies., Astrophys. J.  203 (Jan., 1976) 297–306.
  * [32] R. Takahashi, M. Sato, T. Nishimichi, et al., Revising the Halofit Model for the Nonlinear Matter Power Spectrum, Astrophys. J.  761 (Dec., 2012) 152, [arXiv:1208.2701].
  * [33] D. Blas, J. Lesgourgues, and T. Tram, The Cosmic Linear Anisotropy Solving System (CLASS). Part II: Approximation schemes, Journal of Cosmology and Astro-Particle Physics 7 (July, 2011) 034, [arXiv:1104.2933].
  * [34] M. Tegmark and P. J. E. Peebles, The Time Evolution of Bias, Astrophys. J. Lett.  500 (June, 1998) L79–L82, [astro-ph/9804067].
  * [35] N. A. Reddy, C. C. Steidel, M. Pettini, et al., Multiwavelength Constraints on the Cosmic Star Formation History from Spectroscopy: The Rest-Frame Ultraviolet, H$\alpha$, and Infrared Luminosity Functions at Redshifts $1.9<z<3.4$, Astrophys. J. Suppl.  175 (Mar., 2008) 48–85, [arXiv:0706.4091].
  * [36] M. A. Malkan, D. P. Cohen, M. Maruyama, et al., Lyman-break galaxies at z $\sim$ 3 in the subaru deep field: Luminosity function, clustering, and [o iii] emission, The Astrophysical Journal 850 (nov, 2017) 5, [arXiv:1711.04787].
  * [37] Y. Ono, M. Ouchi, Y. Harikane, et al., Great optically luminous dropout research using subaru HSC (GOLDRUSH). i. UV luminosity functions at z $\sim$ 4–7 derived with the half-million dropouts on the 100 deg2 sky†, Publications of the Astronomical Society of Japan 70 (nov, 2017) [arXiv:1704.06004].
  * [38] P. Madau, Radiative Transfer in a Clumpy Universe: The Colors of High-Redshift Galaxies, Astrophys. J.  441 (Mar., 1995) 18.
  * [39] N. Galitzki, A. Ali, K. S. Arnold, et al., The Simons Observatory: instrument overview, in Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy IX, vol. 10708 of Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, p. 1070804, July, 2018.  arXiv:1808.04493.
  * [40] S.-F. Chen, Z. Vlah, E. Castorina, and M. White, Redshift-Space Distortions in Lagrangian Perturbation Theory, arXiv e-prints (Dec., 2020) arXiv:2012.04636, [arXiv:2012.04636].
  * [41] M. White, R. Zhou, J. DeRose, et al., Cosmological constraints from the tomographic cross-correlation of DESI Luminous Red Galaxies and Planck CMB lensing, Journal of Cosmology and Astro-Particle Physics 2022 (Feb., 2022) 007, [arXiv:2111.09898].
  * [42] S.-F. Chen, Z. Vlah, and M. White, A new analysis of galaxy 2-point functions in the BOSS survey, including full-shape information and post-reconstruction BAO, Journal of Cosmology and Astro-Particle Physics 2022 (Feb., 2022) 008, [arXiv:2110.05530].
  * [43] S.-F. Chen, M. White, J. DeRose, and N. Kokron, Cosmological analysis of three-dimensional BOSS galaxy clustering and Planck CMB lensing cross correlations via Lagrangian perturbation theory, Journal of Cosmology and Astro-Particle Physics 2022 (July, 2022) 041, [arXiv:2204.10392].
  * [44] M. Maus, “An analysis of parameter compression and full-modeling techniques on desi abacus mocks with velocileptors.”.
  * [45] S.-F. Chen, Z. Vlah, and M. White, Consistent Modeling of Velocity Statistics and Redshift-Space Distortions in One-Loop Perturbation Theory, JCAP 07 (2020) 062, [arXiv:2005.00523].
  * [46] S.-F. Chen, Z. Vlah, and M. White, Modeling features in the redshift-space halo power spectrum with perturbation theory, arXiv e-prints (July, 2020) arXiv:2007.00704, [arXiv:2007.00704].
  * [47] R. K. Sheth and G. Tormen, Large-scale bias and the peak background split, Mon. Not. R. Astron. Soc.  308 (Sept., 1999) 119–126, [astro-ph/9901122].
  * [48] S. Cole and N. Kaiser, Biased clustering in the cold dark matter cosmogony, Mon. Not. R. Astron. Soc.  237 (Apr., 1989) 1127–1146.
  * [49] M. White, The Zel’dovich approximation, Mon. Not. R. Astron. Soc.  439 (Apr., 2014) 3630–3640, [arXiv:1401.5466].
  * [50] M. M. Abidi and T. Baldauf, Cubic halo bias in Eulerian and Lagrangian space, Journal of Cosmology and Astro-Particle Physics 2018 (July, 2018) 029, [arXiv:1802.07622].
  * [51] M. Loverde and N. Afshordi, Extended Limber approximation, Phys. Rev. D  78 (Dec., 2008) 123506, [arXiv:0809.5112].
  * [52] B. R. Frieden, Physics from Fisher Information: A Unification.  Cambridge University Press, 1998.
  * [53] M. Tegmark, Measuring cosmological parameters with galaxy surveys, Physical Review Letters 79 (nov, 1997) 3806–3809.
  * [54] M. S. Vogeley and A. S. Szalay, Eigenmode Analysis of Galaxy Redshift Surveys. I. Theory and Methods, Astrophys. J.  465 (July, 1996) 34, [astro-ph/9601185].
  * [55] M. White, Y.-S. Song, and W. J. Percival, Forecasting cosmological constraints from redshift surveys, Mon. Not. R. Astron. Soc.  397 (Aug., 2009) 1348–1354, [arXiv:0810.1518].
  * [56] S.-F. Chen, E. Castorina, M. White, and A. Slosar, Synergies between radio, optical and microwave observations at high redshift, Journal of Cosmology and Astro-Particle Physics 2019 (July, 2019) 023, [arXiv:1810.00911].
  * [57] Euclid Collaboration, A. Blanchard, S. Camera, et al., Euclid preparation. VII. Forecast validation for Euclid cosmological probes, Astron. Astrophys.  642 (Oct., 2020) A191, [arXiv:1910.09273].
  * [58] G. M. Bernstein and Y.-C. Cai, Cosmology without cosmic variance, Mon. Not. R. Astron. Soc.  416 (Oct., 2011) 3009–3016, [arXiv:1104.3862].
  * [59] C. Blake, I. K. Baldry, J. Bland-Hawthorn, et al., Galaxy And Mass Assembly (GAMA): improved cosmic growth measurements using multiple tracers of large-scale structure, Mon. Not. R. Astron. Soc.  436 (Dec., 2013) 3089–3105, [arXiv:1309.5556].
  * [60] T. Mergulhã o, H. Rubira, R. Voivodic, and L. R. Abramo, The effective field theory of large-scale structure and multi-tracer, Journal of Cosmology and Astroparticle Physics 2022 (apr, 2022) 021.
  * [61] T. Mergulhão, H. Rubira, and R. Voivodic, The effective field theory of large-scale structure and multi-tracer ii: redshift space and realistic tracers, 2023.
  * [62] Z. Zheng, R. Cen, H. Trac, and J. Miralda-Escudé, Radiative Transfer Modeling of Ly$\alpha$ Emitters. II. New Effects on Galaxy Clustering, Astrophys. J.  726 (Jan., 2011) 38, [arXiv:1003.4990].
  * [63] C. Behrens, C. Byrohl, S. Saito, and J. C. Niemeyer, The impact of Lyman-$\alpha$ radiative transfer on large-scale clustering in the Illustris simulation, Astron. Astrophys.  614 (June, 2018) A31, [arXiv:1710.06171].
  * [64] R. Momose, K. Shimasaku, K. Nagamine, et al., Catch Me if You Can: Biased Distribution of Ly$\alpha$-emitting Galaxies according to the Viewing Direction, Astrophys. J. Lett.  912 (May, 2021) L24, [arXiv:2104.10580].
  * [65] Planck Collaboration, N. Aghanim, Y. Akrami, et al., Planck 2018 results. VI. Cosmological parameters, arXiv e-prints (July, 2018) arXiv:1807.06209, [arXiv:1807.06209].
  * [66] M. White, Reconstruction within the Zeldovich approximation, Mon. Not. R. Astron. Soc.  450 (July, 2015) 3822–3828, [arXiv:1504.03677].
  * [67] S.-F. Chen, Z. Vlah, and M. White, The reconstructed power spectrum in the Zeldovich approximation, Journal of Cosmology and Astro-Particle Physics 2019 (Sept., 2019) 017, [arXiv:1907.00043].
  * [68] K. N. Abazajian, P. Adshead, Z. Ahmed, et al., CMB-S4 Science Book, First Edition, ArXiv e-prints (Oct., 2016) [arXiv:1610.02743].
  * [69] C. Hahn, M. J. Wilson, O. Ruiz-Macias, et al., The desi bright galaxy survey: Final target selection, design, and validation, The Astronomical Journal 165 (May, 2023) 253.
  * [70] DESI Collaboration, A. Aghamousa, J. Aguilar, et al., The DESI Experiment Part I: Science,Targeting, and Survey Design, ArXiv e-prints (Oct., 2016) [arXiv:1611.00036].
  * [71] R. Zhou, S. Ferraro, M. White, et al., Desi luminous red galaxy samples for cross-correlations, Journal of Cosmology and Astroparticle Physics 2023 (Nov., 2023) 097.
  * [72] N. Mostek, A. L. Coil, M. Cooper, et al., The deep2 galaxy redshift survey: Clustering dependence on galaxy stellar mass and star formation rate atz $\sim$ 1, The Astrophysical Journal 767 (Mar., 2013) 89.
  * [73] A. Liu, J. R. Pritchard, R. Allison, et al., Eliminating the optical depth nuisance from the CMB with 21 cm cosmology, Phys. Rev. D 93 (2016), no. 4 043013, [arXiv:1509.08463].
  * [74] A. Chudaykin and M. M. Ivanov, Measuring neutrino masses with large-scale structure: Euclid forecast with controlled theoretical error, Journal of Cosmology and Astro-Particle Physics 2019 (Nov, 2019) 034, [arXiv:1907.06666].
  * [75] N. Sailer, S.-F. Chen, and M. White, Optical depth to reionization from perturbative 21 cm clustering, JCAP 10 (2022) 007, [arXiv:2205.11504].
  * [76] N. Sailer, S. Ferraro, and M. White, “Measuring neutrino masses from the cosmic frontier.”.
  * [77] H. Sugai, P. A. R. Ade, Y. Akiba, et al., Updated design of the cmb polarization experiment satellite litebird, Journal of Low Temperature Physics 199 (Jan, 2020) 1107–1117.
  * [78] D. J. Watts, B. Wang, A. Ali, et al., A projected estimate of the reionization optical depth using the class experiment’s sample variance limited e-mode measurement, The Astrophysical Journal 863 (Aug., 2018) 121.
  * [79] E. Allys, K. Arnold, J. Aumont, et al., Probing cosmic inflation with the litebird cosmic microwave background polarization survey, Progress of Theoretical and Experimental Physics 2023 (Nov., 2022).
  * [80] E. D. Valentino, T. Brinckmann, M. Gerbino, et al., Exploring cosmic origins with core: Cosmological parameters, Journal of Cosmology and Astroparticle Physics 2018 (Apr., 2018) 017–017.
  * [81] W. Qin, K. Schutz, A. Smith, et al., An Effective Bias Expansion for 21 cm Cosmology in Redshift Space, arXiv:2205.06270.
  * [82] G. Shmueli, D. Sarkar, and E. D. Kovetz, Mitigating the optical depth degeneracy in the cosmological measurement of neutrino masses using 21-cm observations, 2023.
  * [83] Planck Collaboration, Y. Akrami, F. Arroja, et al., Planck 2018 results. I. Overview and the cosmological legacy of Planck, arXiv e-prints (July, 2018) arXiv:1807.06205, [arXiv:1807.06205].
  * [84] Pagano, L., Delouis, J.-M., Mottet, S., et al., Reionization optical depth determination from planck hfi data with ten percent accuracy, A&A 635 (2020) A99.
  * [85] A. Slosar, R. Mandelbaum, and D. Eisenstein, Dark Energy and Modified Gravity, Bull. Am. Astron. Soc.  51 (May, 2019) 97, [arXiv:1903.12016].
  * [86] F. Beutler, M. Biagetti, D. Green, et al., Primordial features from linear to nonlinear scales, Physical Review Research 1 (Dec., 2019) 033209, [arXiv:1906.08758].
  * [87] M. Ballardini, F. Finelli, C. Fedeli, and L. Moscardini, Probing primordial features with future galaxy surveys, JCAP 10 (2016) 041, [arXiv:1606.03747]. [Erratum: JCAP 04, E01 (2018)].
  * [88] V. Poulin, T. L. Smith, T. Karwal, and M. Kamionkowski, Early dark energy can resolve the hubble tension, Physical Review Letters 122 (Jun, 2019).
  * [89] J. C. Hill, E. McDonough, M. W. Toomey, and S. Alexander, Early Dark Energy Does Not Restore Cosmological Concordance, arXiv e-prints (Mar., 2020) arXiv:2003.07355, [arXiv:2003.07355].
  * [90] R. Zhao, X. Mu, R. Gsponer, et al., A multi-tracer analysis for the eboss galaxy sample based on the effective field theory of large-scale structure, 2023\.
  * [91] Y. B. Zel’dovich, Gravitational instability: An approximate theory for large density perturbations., Astron. Astrophys.  5 (Mar., 1970) 84–89.
  * [92] T. Buchert, A class of solutions in Newtonian cosmology and the pancake theory, Astron. Astrophys.  223 (Oct., 1989) 9–24.
  * [93] F. Moutarde, J.-M. Alimi, F. R. Bouchet, et al., Precollapse scale invariance in gravitational instability, Astrophys. J.  382 (Dec., 1991) 377–381.
  * [94] E. Hivon, F. R. Bouchet, S. Colombi, and R. Juszkiewicz, Redshift distortions of clustering: a Lagrangian approach., Astron. Astrophys.  298 (June, 1995) 643, [astro-ph/9407049].
  * [95] A. N. Taylor and A. J. S. Hamilton, Non-linear cosmological power spectra in real and redshift space, Mon. Not. R. Astron. Soc.  282 (Oct., 1996) 767–778, [astro-ph/9604020].
  * [96] V. Zheligovsky and U. Frisch, Time-analyticity of Lagrangian particle trajectories in ideal fluid flow, Journal of Fluid Mechanics 749 (June, 2014) 404–430, [arXiv:1312.6320].
  * [97] T. Matsubara, Recursive solutions of Lagrangian perturbation theory, Phys. Rev. D  92 (July, 2015) 023534, [arXiv:1505.01481].
  * [98] M. McQuinn and M. White, Cosmological perturbation theory in 1+1 dimensions, Journal of Cosmology and Astro-Particle Physics 1 (Jan., 2016) 043, [arXiv:1502.07389].
  * [99] Z. Vlah, E. Castorina, and M. White, The Gaussian streaming model and convolution Lagrangian effective field theory, Journal of Cosmology and Astro-Particle Physics 12 (Dec., 2016) 007, [arXiv:1609.02908].
  * [100] T. Matsubara, Resumming cosmological perturbations via the Lagrangian picture: One-loop results in real space and in redshift space, Phys. Rev. D  77 (Mar., 2008) 063530, [arXiv:0711.2521].
  * [101] T. Matsubara, Nonlinear perturbation theory with halo bias and redshift-space distortions via the Lagrangian picture, Phys. Rev. D  78 (Oct., 2008) 083519, [arXiv:0807.1733].
  * [102] J. Carlson, B. Reid, and M. White, Convolution Lagrangian perturbation theory for biased tracers, Mon. Not. R. Astron. Soc.  429 (Feb., 2013) 1674–1685, [arXiv:1209.0780].
  * [103] R. A. Porto, L. Senatore, and M. Zaldarriaga, The Lagrangian-space Effective Field Theory of large scale structures, Journal of Cosmology and Astro-Particle Physics 5 (May, 2014) 022, [arXiv:1311.2168].
  * [104] Z. Vlah, M. White, and A. Aviles, A Lagrangian effective field theory, Journal of Cosmology and Astro-Particle Physics 9 (Sept., 2015) 014, [arXiv:1506.05264].
  * [105] Z. Vlah, U. Seljak, M. Yat Chu, and Y. Feng, Perturbation theory, effective field theory, and oscillations in the power spectrum, Journal of Cosmology and Astro-Particle Physics 2016 (Mar, 2016) 057, [arXiv:1509.02120].
  * [106] Z. Vlah and M. White, Exploring redshift-space distortions in large-scale structure, Journal of Cosmology and Astro-Particle Physics 2019 (Mar, 2019) 007, [arXiv:1812.02775].
  * [107] N. Kaiser, Clustering in real space and in redshift space, Mon. Not. R. Astron. Soc.  227 (July, 1987) 1–21.
  * [108] U. Seljak and P. McDonald, Distribution function approach to redshift space distortions, Journal of Cosmology and Astro-Particle Physics 11 (Nov., 2011) 039, [arXiv:1109.1888].
  * [109] Z. Vlah, U. Seljak, T. Okumura, and V. Desjacques, Distribution function approach to redshift space distortions. Part V: perturbation theory applied to dark matter halos, Journal of Cosmology and Astro-Particle Physics 10 (Oct., 2013) 053, [arXiv:1308.6294].