# Extracting bias using the cross-bispectrum: An EoR and 21 cm-[CII]-[CII]
case study

Angus Beane11affiliation: Department of Physics & Astronomy, University of
Pennsylvania, 209 South 33rd Street, Philadelphia, PA 19104, USA & Adam
Lidz11affiliation: Department of Physics & Astronomy, University of
Pennsylvania, 209 South 33rd Street, Philadelphia, PA 19104, USA
<EMAIL_ADDRESS>
###### Abstract

The amplitude of redshifted 21 cm fluctuations during the Epoch of
Reionization (EoR) is expected to show a distinctive “rise and fall” behavior
with decreasing redshift as reionization proceeds. On large scales ($k\lesssim
0.1$ Mpc-1) this can mostly be characterized by evolution in the product of
the mean 21 cm brightness temperature and a bias factor, $\langle
T_{21}\rangle b_{21}(z)$. This quantity evolves in a distinctive way that can
help in determining the average ionization history of the intergalactic medium
(IGM) from upcoming 21 cm fluctuation data sets. Here we consider extracting
$\langle T_{21}\rangle b_{21}(z)$ using a combination of future redshifted 21
cm and [CII] line-intensity mapping data sets. Our method exploits the
dependence of the 21 cm-[CII]-[CII] cross-bispectrum on the shape of triangle
configurations in Fourier space. This allows one to determine $\langle
T_{21}\rangle b_{21}(z)$ yet, importantly, is less sensitive to foreground
contamination than the 21 cm auto-spectrum, and so can provide a valuable
cross-check. We compare the results of simulated bispectra with second-order
perturbation theory: on the largest scales well-probed by our simulations
($k\sim 0.05\,\text{Mpc}^{-1}$), the perturbative estimate of $\langle
T_{21}\rangle b_{21}$ matches the true value to within $10\%$ for $\langle
x_{i}\rangle\lesssim 0.8$. The perturbative formula is most accurate early in
the EoR. We consider the 21 cm auto-bispectrum and show that this statistic
may also be used to extract the 21 cm bias factor. Finally, we discuss the
survey requirements for measuring the cross-bispectrum. Although we focus on
the 21 cm-[CII]-[CII] bispectrum during reionization, our method may be of
broader interest and can be applied to any two fields throughout cosmic
history.

###### Subject headings:

cosmology: theory – intergalactic medium – large scale structure of universe

tablenum SIXtablenum

##  1\. Introduction

Observations of the redshifted 21 cm line promise to reveal the timing and
spatial structure of the Epoch of Reionization (EoR) in the near future. This
will help determine the formation time and properties of the first stars,
galaxies, and accreting black holes, and the nature of large-scale structure
at high redshift (Loeb & Furlanetto, 2013). In order to understand the full
astrophysical implications of the upcoming data, a range of challenges must
first be overcome, both to mitigate systematic effects from foreground
contamination and instrumental artifacts, and to faithfully model and
interpret the measurements.

One goal of redshifted 21 cm surveys is to measure the redshift evolution of
the power spectrum of 21 cm fluctuations. On large scales, the amplitude of
the 21 cm power spectrum is expected to rise and fall with decreasing redshift
as reionization proceeds (e.g. Lidz et al. 2008). This can in turn be used to
extract information about the volume-averaged ionization fraction and its
redshift evolution. On scales larger than the size of the ionized bubbles,
this redshift evolution should mostly be characterized by a bias factor,
$b_{21}(z)$, relating the 21 cm fluctuations to fluctuations in the underlying
matter density field on large scales (Furlanetto et al., 2004), and by the
spatial average 21 cm brightness temperature $\langle T_{21}\rangle(z)$. The
bias factor description, if accurate, has the important virtue of providing a
model-independent characterization of the 21 cm fluctuation measurements. The
measured bias factors can then be compared with simulations to extract
information about the ionizing sources and the reionization history.

In conjunction with the 21 cm experiments, a number of efforts are underway to
carry out line-intensity mapping surveys of the EoR in various other emission
lines, such as [CII] (Crites et al., 2014), CO (Chung et al., 2017), and
Ly$\,\alpha$ (Doré et al., 2016). In these line-intensity mapping
observations, one measures the spatial fluctuations in the collective emission
from many individually unresolved sources (Kovetz et al., 2017). Like the 21
cm experiments, these observations span large regions on the sky and provide
accurate redshift information, enabling cross-correlation measurements with
the 21 cm data sets (Lidz et al., 2011; Gong et al., 2012, 2011; Silva et al.,
2015). Surveys in lines such as [CII] and CO should trace large-scale
structure in the galaxy distribution and complement the redshifted 21 cm
experiments, which probe the IGM (Lidz et al., 2009). (Ly$\,\alpha$ provides a
useful intermediate case, reflecting fluctuations in both the source
distribution and that of neutral hydrogen in the IGM, e.g. Pullen et al.
2014.) In addition, the line-intensity 21 cm cross-correlation is less
susceptible to foreground contamination and other systematic effects than the
auto-spectrum extracted from either data set alone (Furlanetto & Lidz, 2007;
Lidz et al., 2009, 2011). The advantage of the cross-spectrum is that residual
21 cm foregrounds from galactic synchrotron, for example, may produce a strong
bias in the auto-spectrum yet do not correlate on average with other tracers
of high redshift structure.

Here we consider an approach for extracting $\langle T_{21}\rangle b_{21}(z)$
from a combination of line-intensity mapping and redshifted 21 cm data. Our
aim is to develop a way of extracting this key quantity from upcoming data
that is less susceptible to systematic contamination than the usual 21 cm
auto-spectrum measurements. In this work we focus on the [CII] line, since
surveys are already underway to detect this line from the EoR (Crites et al.,
2014; Lagache, 2017), but related measurements could be carried out using
other lines.

The method developed here exploits the fact that the growth of large-scale
structure under gravity generates non-Gaussianity, and that the resulting
matter bispectrum has a distinctive dependence on triangle configuration. In
the context of galaxy surveys, this has long been recognized and used to
constrain galaxy bias (e.g. Fry 1994; Matarrese et al. 1997; Scoccimarro et
al. 2001; Verde et al. 2002). In principle, a closely related measurement may
be used to constrain 21 cm biasing as recently discussed by Hoffmann et al.
(2018) in the context of the three-point correlation of the 21 cm field in
configuration space. This recent work follows a number of papers on the
general theme of using the 21 cm bispectrum to complement 21 cm power spectrum
measurements during reionization (Bharadwaj & Pandey, 2005; Shimabukuro et
al., 2016, 2017; Majumdar et al., 2018). In addition to their utility in
constraining 21 cm biasing, higher-order statistics are required to access the
full information content of the highly non-Gaussian 21 cm signal expected
during reionization.

Here we consider a variant of the usual approach for constraining biasing: we
propose to extract $\langle T_{21}\rangle b_{21}(z)$ from a cross-bispectrum
statistic, specifically the 21 cm-[CII]-[CII] bispectrum (involving one 21 cm
field and two [CII] fields). This statistic has the virtue that it is less
sensitive to foreground contamination than the 21 cm auto-spectrum and the 21
cm auto-bispectrum (for similar reasons to the two-point cross-correlation) ,
and can potentially provide a powerful cross-check of inferences from the 21
cm auto-spectrum. We also consider extracting $\langle T_{21}\rangle
b_{21}(z)$ from the 21 cm auto-bispectrum, and we find that this provides
another route for constraining 21 cm biasing, in broad agreement with earlier
work from Hoffmann et al. (2018). In contrast to this previous study, we work
in Fourier space since this is a more natural basis for the interferometric 21
cm measurements. The 21 cm auto-bispectrum is, however, more sensitive to
foreground contamination than the cross-bispectrum advocated here. On the
other hand, it is likely a better probe of reionization’s early phases since
[CII] emission may be dim at early times, i.e. before metal enrichment is well
underway.

The outline of this paper is as follows. In § 2 we present our approach for
measuring the redshift evolution of the 21 cm bias factor. In § 3 we discuss
the reionization simulations used to develop and test our analysis technique
and the algorithm used to measure the bispectrum from simulations. We then
discuss the simulated cross-bispectrum measurements in § 4, and consistency
tests of the perturbative framework in § 5. In § 6 we discuss the prospects
for measuring the cross-bispectra from future surveys. We summarize our
results in § 7. We assume a $\Lambda$CDM cosmology, parameterized by
$(\Omega_{m},\Omega_{\Lambda},\Omega_{b},h,\sigma_{8},n_{s})=(0.27,0.73,0.046,0.7,0.8,1)$,
in broad consistency with recent Planck measurements (Planck Collaboration et
al., 2016a)

##  2\. Constraining 21 cm Bias with the Cross-Bispectrum

Here we explain our technique, describing how the 21 cm bias may be extracted
from the 21 cm-[CII]-[CII] cross-bispectrum. Consider first the 21 cm
brightness temperature field, which describes the brightness temperature
contrast between a neutral hydrogen cloud and the cosmic microwave background:

$T_{21}({\mathchoice{{\hbox{\BFd x}}}{{\hbox{\BFt x}}}{{\hbox{\BFs
x}}}{{\hbox{\BFss x}}}})=T_{0}X_{\text{HI}}({\mathchoice{{\hbox{\BFd
x}}}{{\hbox{\BFt x}}}{{\hbox{\BFs x}}}{{\hbox{\BFss
x}}}})\left[1+\delta_{\rho}({\mathchoice{{\hbox{\BFd x}}}{{\hbox{\BFt
x}}}{{\hbox{\BFs x}}}{{\hbox{\BFss x}}}})\right],$ (1)

Here $T_{0}=28.1745\,\text{mK}\left[(1+z)/10\right]^{1/2}$ (e.g. Zaldarriaga
et al. 2004), $X_{\text{HI}}({\mathchoice{{\hbox{\BFd x}}}{{\hbox{\BFt
x}}}{{\hbox{\BFs x}}}{{\hbox{\BFss x}}}})$ is the neutral hydrogen fraction
(at position x), and $\delta_{\rho}({\mathchoice{{\hbox{\BFd x}}}{{\hbox{\BFt
x}}}{{\hbox{\BFs x}}}{{\hbox{\BFss x}}}})$ is the gas density contrast, which
is assumed to follow the overall matter density field on the large scales of
interest. We can further define the fractional 21 cm brightness temperature
fluctuation, $\delta_{21}({\mathchoice{{\hbox{\BFd x}}}{{\hbox{\BFt
x}}}{{\hbox{\BFs x}}}{{\hbox{\BFss x}}}})=(T_{21}({\mathchoice{{\hbox{\BFd
x}}}{{\hbox{\BFt x}}}{{\hbox{\BFs x}}}{{\hbox{\BFss x}}}})-\langle
T_{21}\rangle)/\langle T_{21}\rangle$. Since $\langle T_{21}\rangle$ is not
directly observable in an interferometric measurement, we will ultimately work
with the statistics of $T_{21}({\mathchoice{{\hbox{\BFd x}}}{{\hbox{\BFt
x}}}{{\hbox{\BFs x}}}{{\hbox{\BFss x}}}})$, or more precisely its Fourier
counterpart, $T_{21}({\mathchoice{{\hbox{\BFd k}}}{{\hbox{\BFt
k}}}{{\hbox{\BFs k}}}{{\hbox{\BFss k}}}})$.111Our Fourier convention is:
$T_{21}({\mathchoice{{\hbox{\BFd k}}}{{\hbox{\BFt k}}}{{\hbox{\BFs
k}}}{{\hbox{\BFss k}}}})=\int\text{d}^{3}x\,T_{21}({\mathchoice{{\hbox{\BFd
x}}}{{\hbox{\BFt x}}}{{\hbox{\BFs x}}}{{\hbox{\BFss
x}}}})e^{i{\mathchoice{{\hbox{\BFd k}}}{{\hbox{\BFt k}}}{{\hbox{\BFs
k}}}{{\hbox{\BFss k}}}}\cdot{\mathchoice{{\hbox{\BFd x}}}{{\hbox{\BFt
x}}}{{\hbox{\BFs x}}}{{\hbox{\BFss x}}}}}$ and
$T_{21}({\mathchoice{{\hbox{\BFd x}}}{{\hbox{\BFt x}}}{{\hbox{\BFs
x}}}{{\hbox{\BFss
x}}}})=\int\frac{\text{d}^{3}k}{(2\pi)^{3}}\,T_{21}({\mathchoice{{\hbox{\BFd
k}}}{{\hbox{\BFt k}}}{{\hbox{\BFs k}}}{{\hbox{\BFss
k}}}})e^{-i{\mathchoice{{\hbox{\BFd k}}}{{\hbox{\BFt k}}}{{\hbox{\BFs
k}}}{{\hbox{\BFss k}}}}\cdot{\mathchoice{{\hbox{\BFd x}}}{{\hbox{\BFt
x}}}{{\hbox{\BFs x}}}{{\hbox{\BFss x}}}}}$. However, it is cleaner to work
with dimensionless fluctuation fields, such as
$\delta_{21}({\mathchoice{{\hbox{\BFd x}}}{{\hbox{\BFt x}}}{{\hbox{\BFs
x}}}{{\hbox{\BFss x}}}})$, in describing biasing relations. Eq. 1 ignores the
impact of peculiar velocities and spin temperature fluctuations, but these
should be good approximations during most of the EoR (e.g. Jensen et al. 2013;
Pritchard & Furlanetto 2007).

Similarly, we can consider the spatial fluctuations in the [CII] emission
field, as probed by upcoming [CII] line-intensity mapping experiments. The
[CII] emission is generally characterized by the three-dimensional specific
intensity field, denoted here by $I_{\text{CII}}({\mathchoice{{\hbox{\BFd
x}}}{{\hbox{\BFt x}}}{{\hbox{\BFs x}}}{{\hbox{\BFss x}}}})$ (Lidz et al.,
2011; Kovetz et al., 2017), or its Fourier partner,
$I_{\text{CII}}({\mathchoice{{\hbox{\BFd k}}}{{\hbox{\BFt k}}}{{\hbox{\BFs
k}}}{{\hbox{\BFss k}}}})$. The fractional fluctuations in the [CII] specific
intensity field are $\delta_{\text{CII}}({\mathchoice{{\hbox{\BFd
x}}}{{\hbox{\BFt x}}}{{\hbox{\BFs x}}}{{\hbox{\BFss
x}}}})=(I_{\text{CII}}({\mathchoice{{\hbox{\BFd x}}}{{\hbox{\BFt
x}}}{{\hbox{\BFs x}}}{{\hbox{\BFss x}}}})-\langle
I_{\text{CII}}\rangle)/\langle I_{\text{CII}}\rangle$.

As we will test subsequently using simulations of the EoR, we assume that on
large scales each of these fluctuation fields is a deterministic and local
function of the underlying matter density distribution (Scherrer & Weinberg,
1998; Dekel & Lahav, 1999). For example, we suppose that
$\delta_{21}({\mathchoice{{\hbox{\BFd x}}}{{\hbox{\BFt x}}}{{\hbox{\BFs
x}}}{{\hbox{\BFss x}}}})=f\left[\delta_{\rho}({\mathchoice{{\hbox{\BFd
x}}}{{\hbox{\BFt x}}}{{\hbox{\BFs x}}}{{\hbox{\BFss x}}}})\right]$, where $f$
is the biasing function that specifies the relation between the 21 cm and
matter density fluctuations. It is not obvious that the local biasing
assumption is adequate for describing 21 cm fluctuations during reionization,
since distant sources may impact the ionization state of the gas – especially
at the late stages of reionization. However, as we will quantify, local
biasing may nevertheless be a good approximation if we consider sufficiently
large scales. Considering such scales, we Taylor expand the biasing function
for small $\delta_{\rho}({\mathchoice{{\hbox{\BFd x}}}{{\hbox{\BFt
x}}}{{\hbox{\BFs x}}}{{\hbox{\BFss x}}}})$, keeping terms to second-order in
the density fluctuations:

$\delta_{21}({\mathchoice{{\hbox{\BFd x}}}{{\hbox{\BFt x}}}{{\hbox{\BFs
x}}}{{\hbox{\BFss x}}}})=b_{21}\delta_{\rho}({\mathchoice{{\hbox{\BFd
x}}}{{\hbox{\BFt x}}}{{\hbox{\BFs x}}}{{\hbox{\BFss
x}}}})+\frac{b_{21}^{(2)}}{2}\delta^{2}_{\rho}({\mathchoice{{\hbox{\BFd
x}}}{{\hbox{\BFt x}}}{{\hbox{\BFs x}}}{{\hbox{\BFss x}}}}).$ (2)

This equation defines the first-order and second-order 21 cm bias factors,
$b_{21}$ and $b_{21}^{(2)}$. (For brevity of notation, we do not include a
super-script (1) in our notation for the first-order bias factors.) The
Fourier space counterpart of Eq. 2 is,

$\delta_{21}({\mathchoice{{\hbox{\BFd k}}}{{\hbox{\BFt k}}}{{\hbox{\BFs
k}}}{{\hbox{\BFss k}}}})=b_{21}\delta_{\rho}({\mathchoice{{\hbox{\BFd
k}}}{{\hbox{\BFt k}}}{{\hbox{\BFs k}}}{{\hbox{\BFss
k}}}})+\frac{b_{21}^{(2)}}{2}\int\frac{d^{3}q}{(2\pi)^{3}}\delta_{\rho}({\mathchoice{{\hbox{\BFd
q}}}{{\hbox{\BFt q}}}{{\hbox{\BFs q}}}{{\hbox{\BFss
q}}}})\delta_{\rho}({\mathchoice{{\hbox{\BFd k}}}{{\hbox{\BFt k}}}{{\hbox{\BFs
k}}}{{\hbox{\BFss k}}}}-{\mathchoice{{\hbox{\BFd q}}}{{\hbox{\BFt
q}}}{{\hbox{\BFs q}}}{{\hbox{\BFss q}}}}),$ (3)

i.e. the $\delta^{2}_{\rho}({\mathchoice{{\hbox{\BFd x}}}{{\hbox{\BFt
x}}}{{\hbox{\BFs x}}}{{\hbox{\BFss x}}}})$ term in configuration space becomes
a convolution in Fourier space. In practice, since $\langle T_{21}\rangle$ is
not directly observable in an interferometric measurement, we generally
consider the bias factors multiplied by $\langle T_{21}\rangle$. Under the
assumption of a local and deterministic biasing relation for the [CII]
specific intensity field, we define analogous bias factors $b_{\text{CII}}$
and $b_{\text{CII}}^{(2)}$.

As large-scale structure grows under gravity, the matter density field will
develop non-Gaussianity. At second-order in perturbation theory, the density
fluctuations $\delta_{\rho}({\mathchoice{{\hbox{\BFd k}}}{{\hbox{\BFt
k}}}{{\hbox{\BFs k}}}{{\hbox{\BFss k}}}})$ follow (e.g. Bernardeau et al.
2002):

$\delta_{\rho}(\bm{k})=\delta_{\rho}^{(1)}({\mathchoice{{\hbox{\BFd
k}}}{{\hbox{\BFt k}}}{{\hbox{\BFs k}}}{{\hbox{\BFss
k}}}})+\int\frac{\text{d}^{3}q}{(2\pi)^{3}}F_{2}(\bm{q},\bm{k}-\bm{q})\delta_{\rho}^{(1)}(\bm{q})\delta_{\rho}^{(1)}(\bm{k}-\bm{q}),$
(4)

where $\delta_{\rho}^{(1)}$ denotes the first-order fluctuation field and
$F_{2}({\mathchoice{{\hbox{\BFd k}}}{{\hbox{\BFt k}}}{{\hbox{\BFs
k}}}{{\hbox{\BFss k}}}}_{1},{\mathchoice{{\hbox{\BFd k}}}{{\hbox{\BFt
k}}}{{\hbox{\BFs k}}}{{\hbox{\BFss k}}}}_{2})$ describes mode-coupling from
non-linear growth. In a flat $\Lambda$CDM universe with $\Omega_{m}(z)\approx
1$, appropriate for high redshifts, the mode-coupling kernel is given by:

$F_{2}(\bm{q}_{1},\bm{q}_{2})=\frac{5}{7}+\frac{\bm{q_{1}}\cdot\bm{q}_{2}}{2q_{1}q_{2}}\left(\frac{q_{1}}{q_{2}}+\frac{q_{2}}{q_{1}}\right)+\frac{2}{7}\left(\frac{\bm{q}_{1}\cdot\bm{q}_{2}}{q_{1}q_{2}}\right)^{2}.$
(5)

This mode-coupling sources the matter density bispectrum, which is defined by:

$\begin{split}\langle\delta_{\rho}(\bm{k}_{1})\delta_{\rho}(\bm{k}_{2})\delta_{\rho}(\bm{k}_{3})\rangle&\equiv(2\pi)^{3}\delta_{D}(\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3})\\\
&\times
B_{\delta,\delta,\delta}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})\text{.}\end{split}$
(6)

Using the result of second-order perturbation theory, Eqs. 4 and 5 gives the
bispectrum to lowest non-vanishing order as:

$\begin{split}B^{(0)}_{\delta,\delta,\delta}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})=&2F_{2}(\bm{k}_{1},\bm{k}_{2})P_{\delta,\delta}^{\text{lin}}(k_{1})P_{\delta,\delta}^{\text{lin}}(k_{2})\\\
&+\text{2 perm.},\end{split}$ (7)

where $P^{\text{lin}}_{\delta,\delta}$ is the linear matter auto-spectrum and
2 perm. indicates permutations of the products of power spectra that enter,
with each permutation picking out a product of two of the three wavevectors
involved. We make a distinction here between the linear and nonlinear matter
power spectra, but there is only a $\sim 5\%$ difference between the two for
the smallest scales ($k\sim 0.4\,\text{Mpc}^{-1}$) and lowest redshift ($z\sim
6$) we consider. Since $B\propto P^{2}$, it is also convenient to define the
reduced bispectrum. In the case of the matter density field,

$Q_{\delta,\delta,\delta}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})=\frac{B_{\delta,\delta,\delta}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})}{P_{\delta,\delta}^{\text{nl}}(k_{1})P_{\delta,\delta}^{\text{nl}}(k_{2})+\text{2
perm.}}.$ (8)

It is common to specify $Q$ by the magnitude of wavevectors,
$k_{1}=|\bm{k}_{1}|$, and $k_{2}=|\bm{k}_{2}|$, and the angle between
$\bm{k}_{1}$ and $\bm{k}_{2}$, $\theta_{12}$.222We will use a different
Fourier space characterization when we actually extract $b_{21}$ from our
simulations, see § 3.4 The reduced matter bispectrum has a distinctive form:
first, gravitational mode-coupling (Eqs. 4-8) enhances the small-scale power
spectrum in large-scale overdense regions compared to that in underdense
regions. Second, the reduced bispectrum is smaller for triangles that are
close to isosceles with $\theta_{12}\sim\pi/2$, than for triangles that are
nearly co-linear with $\theta_{12}\sim 0,\pi$. This is a consequence of the
filamentary nature of large-scale structure, and leads to a characteristic
“U”-shaped dependence of $Q_{\delta,\delta,\delta}$ on $\theta_{12}$ (e.g.
Bernardeau et al. 2002).

Turning to biased-tracers, such as the 21 cm or [CII] fluctuations, it is
clear that (at lowest non-vanishing order) the bispectra of these fields
receive contributions both from non-Gaussianity in the matter distribution
(described by Eq. 8) and owing to non-linearities in the biasing relation
(Eqs. 2-3). Here we consider the cross-bispectrum between the 21 cm field and
two copies of the [CII] fluctuation field. This quantity is less sensitive to
foreground contamination – and to unshared systematics – than the auto-
bispectra of these fields. This is because the foregrounds in the two surveys
should be mostly uncorrelated (asides for common foregrounds at the widely
separated observing frequencies of the two experiments), and so should not
contribute to the ensemble-averaged cross-bispectrum.333Any correlations
between the signal in one survey and the foregrounds in the other could also
produce a small ensemble-averaged cross-bispectrum. For example, the high
redshift [CII] emitting galaxies – or other correlated sources at the redshift
of these galaxies – may emit synchrotron radiation which constitutes a (very
small) “foreground” for the 21 cm survey. Foreground cleaning or avoidance is
still important because residual foregrounds will increase the variance of a
cross-bispectrum estimate, but the cross-bispectrum has the important virtue
that unshared foregrounds will not produce a spurious signal on average. The
cross-power spectrum is also less sensitive to foreground contamination;
however on large scales, the cross-spectrum is proportional to the overall
product $b_{21}\langle T_{21}\rangle b_{\text{CII}}\langle
I_{\text{CII}}\rangle$. That is, the cross-power spectrum depends on both the
21 cm and [CII] biasing relations, as well as the mean intensity in each line.
On the other hand, we will show that the shape of the cross-bispectrum – when
suitably defined – is sensitive only to $\langle T_{21}\rangle b_{21}$ (and
independent of the [CII] biasing and mean intensity).

Specifically, consider the cross-bispectrum between a single 21 cm field and
two [CII] emission fields:

$\begin{split}\langle
T_{21}(\bm{k}_{1})I_{\text{CII}}(\bm{k}_{2})I_{\text{CII}}(\bm{k}_{3})\rangle\equiv&\,(2\pi)^{3}\delta_{D}(\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3})\\\
&\times
B_{21,\text{CII},\text{CII}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3}).\end{split}$
(9)

Note that we have now switched from considering $\delta_{21}$ and
$\delta_{\text{CII}}$ to $T_{21}$ and $I_{\text{CII}}$, since only the latter
quantities are directly observable from interferometric and line-intensity
mapping measurements. Apart from the overall $\bm{k}=0$ mode,
$T_{21}(\bm{k})=\langle T_{21}\rangle\delta_{21}(\bm{k})$, and
$I_{\text{CII}}(\bm{k})=\langle
I_{\text{CII}}\rangle\delta_{\text{CII}}(\bm{k})$.

We can define a reduced cross-bispectrum that is formed from Eq. 9 and the
cross-power spectrum between the 21 cm and [CII] fields,
$P_{21,\text{CII}}(k_{i})$, at each of the three wavevectors $k_{i}$ with
$i=1,2,3$. Specifically,

$\hat{Q}_{21,\text{CII},\text{CII}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})=\frac{B_{21,\text{CII},\text{CII}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})}{P_{21,\text{CII}}(k_{1})P_{21,\text{CII}}(k_{2})+\text{2
perm.}}$ (10)

Note that we deliberately construct the reduced cross-bispectrum using only
the cross-power spectrum in the denominator of Eq. 10: explicitly, “2 perm.”
is specified by
$P_{21,\text{CII}}(k_{2})P_{21,\text{CII}}(k_{3})+P_{21,\text{CII}}(k_{3})P_{21,\text{CII}}(k_{1})$
and does not involve the auto-power spectrum of either field. Similar
definitions hold for $\hat{Q}_{\text{CII},21,\text{CII}}$ and
$\hat{Q}_{\text{CII},\text{CII},21}$; these merely rearrange which wavevector
is attached to the 21 cm field. Note that with this definition, using ${\rm
mK}$ units for $T_{21}$ and ${\rm Jy/str}$ units for $I_{\text{CII}}$, the
cross-bispectrum $B_{21,\text{CII},\text{CII}}$ has units of ${\rm mK}({\rm
Jy/str})^{2}({\rm Mpc})^{6}$, while the reduced cross-bispectrum
$\hat{Q}_{21,\text{CII},\text{CII}}$ has units of ${\rm mK}^{-1}$.

$\displaystyle\hat{Q}_{21,\text{CII},\text{CII}}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})=\frac{Q_{\delta,\delta,\delta}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})}{\langle
T_{21}\rangle b_{21}}$ $\displaystyle+\frac{b_{21}^{(2)}}{2\langle
T_{21}\rangle b_{21}^{2}}L(k_{1},k_{2})$
$\displaystyle+\frac{b_{\text{CII}}^{(2)}}{2\langle T_{21}\rangle
b_{21}b_{\text{CII}}}\big{[}L(k_{2},k_{3})+L(k_{3},k_{1})\big{]}\text{,}$

where we have defined

$L(k_{i},k_{j})=\frac{P_{1,2}(k_{i})P_{1,2}(k_{j})}{P_{1,2}(k_{1})P_{1,2}(k_{2})+\text{2
perm.}}\text{.}$ (11)

As will become clear shortly, these $L$ terms are inconvenient, but note that

$L(k_{1},k_{2})+L(k_{2},k_{3})+L(k_{3},k_{1})=1\text{.}$ (12)

Therefore, summing permutations of Eq. 2 and using Eq. 12 rids us of the $L$
terms:

$Q_{21,\text{CII},\text{CII}}\equiv\frac{1}{3}\left(\hat{Q}_{21,\text{CII},\text{CII}}+\hat{Q}_{\text{CII},21,\text{CII}}+\hat{Q}_{\text{CII},\text{CII},21}\right)\text{.}$
(13)

Throughout we will only work with $Q_{1,2,2}$, and so these permutations are
implicit in what follows. Note that the permutations in Eq. 13 imply that
swapping any two of the wavevector arguments results in the same value of $Q$.
Therefore, we will also implicitly enforce $k_{1}\geq k_{2}\geq k_{3}$ where
relevant, to avoid double counting.

With these permutations in mind, we now have the following formula:

$\begin{split}Q_{21,\text{CII},\text{CII}}=&\frac{Q_{\delta,\delta,\delta}}{\langle
T_{21}\rangle b_{21}}+\frac{1}{6}\frac{b_{21}^{(2)}}{\langle T_{21}\rangle
b_{21}^{2}}\\\ &+\frac{1}{3}\frac{b_{\text{CII}}^{(2)}}{\langle T_{21}\rangle
b_{21}b_{\text{CII}}}\text{.}\end{split}$ (14)

As we have described, this formula is valid at second-order in perturbation
theory. In our actual analysis, we will use the simulated
$Q_{\delta,\delta,\delta}$ instead of the perturbative
$Q_{\delta,\delta,\delta}$. This is mostly to minimize the impact of sample
variance and is discussed further in § 4.

Our analysis strategy is now clear, and analogous to methods developed
previously to constrain galaxy bias (Fry, 1994; Matarrese et al., 1997;
Scoccimarro et al., 2001; Verde et al., 2002). The first term in
$Q_{21,\text{CII},\text{CII}}$ depends on the reduced density bispectrum, with
its distinctive “U”-shaped dependence on $\theta_{12}$, and the 21 cm linear
bias factor. In this term, only the 21 cm bias enters and so this term is
entirely independent of the statistical properties of the [CII] emission
fluctuations. In contrast, the other terms depend on the first and second-
order bias factors of both the 21 cm and [CII] fields, and so these factors
are harder to interpret, but produce only a constant offset. The rationale for
taking the permuted version (Eq. 13) is now clarified, since as a result the
[CII] bias factors enter only as an overall constant, independent of
triangular configuration. In summary, by measuring the triangular shape
dependence of $Q_{21,\text{CII},\text{CII}}$, we should be able to extract
$\langle T_{21}\rangle b_{21}$.

In testing the perturbative framework described above, we will consider
various additional bispectra. At lowest non-vanishing order in perturbation
theory, these all have the general form:

$Q_{21,\text{X},\text{X}}^{(0)}=\frac{Q_{\delta,\delta,\delta}}{\langle
T_{21}\rangle b_{21}}+C_{21,\text{X},\text{X}},$ (15)

where the various values of $C$ depend on which field “X” is being combined
with the 21 cm fluctuations. Specifically,

$C_{21,21,21}=\frac{1}{2}\frac{b_{21}^{(2)}}{\langle T_{21}\rangle
b_{21}^{2}},$ (16)
$C_{21,\delta,\delta}=\frac{1}{6}\frac{b_{21}^{(2)}}{\langle T_{21}\rangle
b_{21}^{2}},$ (17)

and

$C_{21,\text{CII},\text{CII}}=\frac{1}{6}\frac{b_{21}^{(2)}}{\langle
T_{21}\rangle b_{21}^{2}}+\frac{1}{3}\frac{b_{\text{CII}}^{(2)}}{\langle
T_{21}\rangle b_{21}b_{\text{CII}}}\text{.}$ (18)

##  3\. Reionization Simulations and Methodology

In order to characterize the cross-bispectrum signal at different stages of
reionization, and to test the accuracy of the perturbative formulas from the
previous section, we turn to semi-numerical simulations of the EoR (Zahn et
al., 2006; Mesinger & Furlanetto, 2007), based on the excursion-set (Bond et
al., 1991) model of reionization (Furlanetto et al., 2004). Specifically, we
use the publicly available 21cmFAST code v1.12 (Mesinger et al., 2011).

Our simulations are unlikely to provide fully accurate models of the small-
scale 21 cm signal. In the context of this paper, the small-scale modes mostly
set the overall values of the bias factors ($\langle T_{21}\rangle b_{21}$) at
different redshifts and neutral fractions. The small-scale modes should
therefore impact the precise evolution of the bias factors with redshift, and
not the overall framework we propose here.

###  3.1. L-PICOLA

By default, 21cmFAST uses the Zel’dovich approximation (ZA) to generate the
density field (Zel’dovich, 1970). This is appropriate for capturing the large-
scale two-point statistics of the 21 cm field at high redshifts, but it is
inadequate for modeling three-point statistics. The ZA is known to
underestimate the density bispectrum (Scoccimarro, 1997; Leclercq et al.,
2013), and we have verified that it provides a poor description of the density
bispectrum at the scales and redshifts of interest for our study. We thus turn
to the publicly available L-PICOLA v1.3 code (Howlett et al., 2015; Tassev et
al., 2013), which is a hybrid between a 2nd order Lagrangian perturbation
theory (2LPT) and a particle mesh (PM) code. L-PICOLA is a parallel version of
the COLA method, and accurately reproduces the two and three-point statistics
of interest of our study, yet is less computationally expensive than full
GADGET-2 (Springel, 2005) runs.

In order to accurately estimate the bispectra starting at $z\sim 10.5$, two
potential issues are: transients associated with the initial conditions
(Crocce et al., 2006) and discreteness noise from the finite number of
particles in the simulation. These issues are more of a concern at the high
redshifts of interest for our study (than near $z=0$), since there is less
time for transients to relax and because the density fluctuations are small at
early times. The transients are suppressed in part by adopting 2LPT initial
conditions (as opposed to Zel’dovich initial conditions) and by starting the
simulation early. We guard against shot-noise by simulating a fairly high
particle density. Specifically, we simulate $1024^{3}$ particles in a
simulation box with a co-moving side length of $L_{\rm box}=800\,\text{Mpc}$.
We chose an extremely high starting redshift of $z_{i}=5000$.444The starting
redshift is artificially high; the initial conditions are set by rescaling
from the $z=0$ matter power spectrum assuming that the universe is dominated
by pressure-free matter at the initial redshift (even though this redshift is
formally before matter-radiation equality.) The subsequent evolution is
tracked by L-PICOLA and this returns the correct linear evolution (under the
approximation that the universe is entirely composed of pressure-free matter)
after recombination. The early start was just meant to guard against initial
transients and to ensure the correct density statistics at $z\sim 6-11$. The
code takes 300 times steps between $z_{i}$ and $z=10.49$.

In order to beat-down noise in our bispectrum estimates on large scales, we
run ten different realizations of the L-PICOLA simulations, resulting in a
total simulated volume of $5.12\,\text{Gpc}^{3}$. The boxsize was chosen so
that the fundamental mode of the simulation box is much larger than the scales
of interest for our study. Specifically, the largest scale considered is
$k=0.05\,\text{Mpc}^{-1}$ and the fundamental mode of the simulation box is
$k_{f}=0.00785\,\text{Mpc}^{-1}$. We save simulation snapshots at redshifts of
$z=10.49$, $9.41$, $8.43$, $7.88$, $7.37$, $7.05$, $6.73$, $6.43$, and $6.00$,
and take ten time steps between each output. For each redshift snapshot, we
use nearest grid-point interpolation to estimate the density field on a
$512^{3}$ Cartestian grid from the particle positions.

###  3.2. 21cmFAST

The $512^{3}$ particle, $L_{\rm box}=800\,\text{Mpc}$, gridded density fields
from the L-PICOLA runs are then passed to the 21cmFAST code to generate
simulated ionization and 21 cm fields at each redshift. Specifically, the
excursion set approach is used to identify ionized regions (Mesinger et al.,
2011) using the following simulation parameters: the ionizing efficiency is
set to $\zeta=10$, the minimum virial temperature of galaxy hosting halos is
$T_{\text{vir}}=10^{4}$ K; the largest smoothing scale – for generating the
ionization field using the excursion-set methodology – is $R_{\text{max}}=30$
co-moving Mpc. (The latter quantity is sometimes loosely termed the “mean-free
path”, $R_{\text{mfp}}$, in the literature.)

Figure 1.— The ionization history and the redshift evolution of the 21 cm bias
factor. Top panel: The volume-averaged neutral hydrogen fraction, $\langle
x_{\text{HI}}\rangle$, as a function of redshift from our 21cmFAST simulation.
Bottom panel: The corresponding evolution of the 21 cm bias factor, $\langle
T_{21}\rangle b_{21}$, which rises and falls with decreasing redshift and
peaks (in amplitude) around reionization’s midpoint, when $\langle
x_{i}\rangle\sim 0.5$ – or slightly thereafter – of the IGM volume is
reionized. We compute $\langle T_{21}\rangle b_{21}$ from the 21 cm-density
cross-spectrum at various $k$. At high redshifts (early in reionization)
$\langle T_{21}\rangle b_{21}$ is insensitive to the precise wavenumber
considered, but differences appear as reionization proceeds. This provides one
indication of the redshifts and scales where the linear biasing description
becomes imperfect.

The Top panel of Fig. 1 shows the redshift evolution of the volume-averaged
neutral fraction, $\langle x_{\text{HI}}\rangle$, in our 21cmFAST model. In
this model, reionization completes shortly after $z\sim 6$, while $50\%$ of
the volume is ionized at $z\sim 7.5$, and $20\%$ is ionized by $z\sim 10$;
this is broadly consistent with current constraints on the EoR (e.g. Becker et
al. (2001); Fan et al. (2006); Venemans et al. (2013); Becker et al. (2015);
Planck Collaboration et al. (2016b). ) The Bottom panel shows the
corresponding evolution of the 21 cm bias factor with redshift, $\langle
T_{21}\rangle b_{21}$ (estimated from the 21 cm density cross-power spectrum
as we detail subsequently) calculated at each of
$k=0.05\,\text{Mpc}^{-1},0.1\,\text{Mpc}^{-1}$, and $0.2\,\text{Mpc}^{-1}$.
Fig. 1 exemplifies the distinctive “rise and fall” evolution in the magnitude
of $\langle T_{21}\rangle b_{21}$: a key goal of upcoming 21 cm surveys is to
measure this evolution, and to use it to extract the corresponding ionization
history. Reiterating, the motivation of our paper is to develop an approach
for extracting this evolution which is less susceptible to foreground
contamination than the 21 cm auto spectrum. At early times, the bias factor
inferred is identical for all three spatial scales considered. However, as
reionization proceeds the bias factor starts to depend on wavenumber: this
occurs because the ionized regions grow in size during reionization, and the
linear biasing description breaks down on progressively larger scales. Hence,
even on the rather large scales spanned by our simulation box we expect linear
biasing to be an imperfect description during the late stages of reionization
(as we will quantify further subsequently). In practice, it will probably be
challenging to access still larger scales from future surveys (where linear
biasing might be a better description), since foreground avoidance or cleaning
limit the prospects for extracting very large scale modes. (See, e.g., Fig. 1
of Dillon et al. 2014 and the discussion therein.)

Lastly, note that the bias factor in the Bottom panel of Fig. 1 is negative:
this indicates that large-scale overdense regions are dimmer in 21 cm emission
than typical regions. The negative sign results because the ionizing sources
form first in large-scale overdense regions and so these regions are ionized
early. It is worth remarking that the 21 cm auto-power spectrum is sensitive
only to $\left[\langle T_{21}\rangle b_{21}\right]^{2}$ and hence can not
distinguish the overall sign here. On the other hand, the 21 cm-[CII] cross-
spectrum, the 21 cm auto-bispectrum, and the 21 cm-[CII]-[CII] cross-
bispectrum each depend on both the magnitude and the sign of $\langle
T_{21}\rangle b_{21}$, as has also been shown by Majumdar et al. (2018);
Hoffmann et al. (2018).

Figure 2.— Linear biasing of the 21 cm fluctuations and its range of validity.
Top panel: Estimates of the scale-dependent bias factor, $\langle
T_{21}\rangle b_{21}(k)$, from our reionization simulation for several
different redshifts and ionization fractions ($\langle x_{i}\rangle$), as
described in the legend. The solid lines show estimates from the 21 cm auto-
power spectrum, while the dashed lines show the bias factor inferred from the
21 cm-density cross-power spectrum. Bottom panel: The cross-correlation
coefficient between the 21 cm and matter density fields, $r_{21,\delta}(k)$.
At $k\lesssim 0.1\,\text{Mpc}^{-1}$, the bias factor is close to scale-
independent, and $r_{21,\delta}(k)\sim-1$.

We can carry out further tests of the linear biasing approximation using
additional two-point statistical measures. As already mentioned, we expect
this approximation to break down on sufficiently small scales, especially
during the later stages of reionization. On scales smaller than the size of
the ionized regions, there are order unity fluctuations in the ionization
fraction, as some areas in the IGM are nearly completely ionized while others
remain almost completely neutral. In addition, the ionization and density
field de-correlate on scales smaller than the size of the ionized bubbles
(e.g. Furlanetto et al. 2004; Zahn et al. 2006). Hence on these scales, the 21
cm and density fluctuation fields do not share the same phases and an
expansion along the lines of Eq. 2 is inadequate. To explore this we estimated
the bias factor from the cross-power spectrum of the 21 cm and density fields,
$b_{21,x}(k)=P_{21,\delta}(k)/P_{\delta,\delta}(k)$, and from the 21 cm auto-
power spectrum
$b_{21,a}(k)=\left[P_{21,21}(k)/P_{\delta,\delta}(k)\right]^{1/2}$ (Top panel
of Fig. 2). Note that these two estimates should agree in the linear-biasing
regime, but may depart from each other on smaller scales (at higher $k$).
Furthermore, we calculate the cross-correlation coefficient between the 21 cm
and density fields,
$r_{21,\delta}(k)=P_{21,\delta}(k)/\sqrt{P_{21,21}(k)P_{\delta,\delta}(k)}$
(Bottom panel). The cross-correlation coefficient quantifies how well the
phases of the two fields track each other.

At $k\lesssim 0.1\,\text{Mpc}^{-1}$, for the redshifts shown: the correlation-
coefficient between the two fields is $r_{21,\delta}\sim-1$, and the bias
factor is fairly scale-independent. The auto/cross-power spectrum approaches
for calculating the bias factor also agree with each other to within $\sim
20\%$ for $k\sim 0.1\,\text{Mpc}^{-1}$, and to within $\sim 5\%$ for $k\sim
0.05\,\text{Mpc}^{-1}$ at $\langle x_{i}\rangle=0.8$ ($z=6.43$). The linear
biasing description appears best in the middle of reionization. Towards the
end of reionization, as illustrated by the blue $\langle x_{i}\rangle=0.80$
curves, the bubbles are large enough that $k\sim 0.1\,\text{Mpc}^{-1}$ is only
marginally in the linear biasing regime. On the other hand, in the earliest
phases of reionization ($\langle x_{i}\rangle\sim 0.1$), the correlation-
coefficient departs significantly from $r_{21,\delta}\sim-1$: this likely
results because the sign of $r_{21,\delta}$ is initially positive and
gradually reverses (Lidz et al., 2008). This occurs because the large-scale
overdense regions prior to reionization contain more neutral hydrogen and are
brighter in 21 cm emission than underdense areas, while the large-scale
overdense regions ionize first and subsequently become dimmer in 21 cm. In our
model large-scale overdense and underdense regions are nearly at the same
brightness temperature at $\langle x_{i}\rangle\sim 0.1$ while the reversal is
almost complete by $\langle x_{i}\rangle=0.22$ at $z=9.41$ (the earliest case
shown in Fig. 2). We do not expect our perturbative formula to be accurate in
the brief earlier phase near $\langle x_{i}\rangle\sim 0.1$, where
$r_{21,\delta}$ changes sign rapidly. Note also that we neglect spin
temperature fluctuations, which should impact our predictions at early times.
In general, however, our results suggest that linear biasing is a fairly good
description at $k\lesssim 0.1\,\text{Mpc}^{-1}$ during the bulk of
reionization, further motivating us to test our cross-bispectrum method.

###  3.3. [CII] Intensity Field

In order to model the [CII] intensity field, our simulations would ideally
resolve halos with masses down to the atomic cooling limit ($\sim
10^{8}\,M_{\sun}$). In this case, we could model the high redshift [CII]
intensity field by assuming correlations between halo properties and [CII]
luminosity. However, it is difficult to resolve such small mass halos while
capturing a large enough volume to accurately model the bispectrum. For
example, the particle mass in our simulations is $\sim 10^{10}M_{\sun}$ and so
we don’t resolve many of the small mass halos which may host [CII] emitting
galaxies. Consequently, we defer a full treatment of [CII] emission to future
work.

Note, however, that the perturbative formulae (Eq. 14 and Eq. 15) predict –
apart from the overall constant term – that the $Q_{21,\text{CII},\text{CII}}$
bispectrum should be identical to that of the $Q_{21,\delta,\delta}$
bispectrum. That is, both of these cross-bispectra are expected to follow the
$Q_{\delta,\delta,\delta}/(\langle T_{21}\rangle b_{21})$ form. Consequently,
we study the statistics of $Q_{21,\delta,\delta}$ here rather than
$Q_{21,\text{CII},\text{CII}}$ since the latter is more difficult to model.
The main shortcoming of this approach is that it neglects the additional non-
linearities in the [CII] biasing relation, and these effects may reduce the
agreement with the perturbative formula. In estimating the signal-to-noise
ratio at which $Q_{21,\text{CII},\text{CII}}$ may be measured in future
surveys (§ 6), we assume a linear biasing model for [CII] emission. This
allows us to predict the [CII] auto-spectrum and [CII]-21 cm cross-spectrum
(as required for the signal-to-noise ratio calculation) from the statistics of
the simulated density field. It may be possible to construct a sub-grid model
to include unresolved small halos with the correct statistical properties
(e.g. McQuinn et al. (2007)). However, the approach adopted in that work would
need to be extended to capture the correct quadratic halo bias and is beyond
the scope of this paper.

###  3.4. Estimating the Bispectrum

In order to estimate the bispectrum from our 21cmFAST simulation outputs, we
follow the method described in Smith et al. (2008). Throughout we consider two
characterizations of the wavevector arguments. In both cases we calculate the
average bispectrum in linearly-spaced wavevector bins.

Our first characterization is for illustration purposes only. Here we fix the
magnitude of $|{\mathchoice{{\hbox{\BFd k}}}{{\hbox{\BFt k}}}{{\hbox{\BFs
k}}}{{\hbox{\BFss k}}}}_{1}|$ and $|{\mathchoice{{\hbox{\BFd k}}}{{\hbox{\BFt
k}}}{{\hbox{\BFs k}}}{{\hbox{\BFss k}}}}_{2}|$ while varying the angle
$\theta_{12}$ between the two wavevectors. The third side of the triangle,
${\mathchoice{{\hbox{\BFd k}}}{{\hbox{\BFt k}}}{{\hbox{\BFs k}}}{{\hbox{\BFss
k}}}}_{3}$, is determined by these three parameters and the requirement of
translation invariance (${\mathchoice{{\hbox{\BFd k}}}{{\hbox{\BFt
k}}}{{\hbox{\BFs k}}}{{\hbox{\BFss k}}}}_{1}+{\mathchoice{{\hbox{\BFd
k}}}{{\hbox{\BFt k}}}{{\hbox{\BFs k}}}{{\hbox{\BFss
k}}}}_{2}+{\mathchoice{{\hbox{\BFd k}}}{{\hbox{\BFt k}}}{{\hbox{\BFs
k}}}{{\hbox{\BFss k}}}}_{3}=\bm{0}$.) We consider spherical bins around each
wavevector of thickness $\Delta k=k_{f}$, where
$k_{f}=2\pi/L_{\text{box}}=0.00785\,\text{Mpc}^{-1}$ is the fundamental mode
of our simulation box.

Our second characterization sets a minimum ($k_{\text{min}}$) and a maximum
($k_{\text{max}}$) magnitude for each wavevector argument. That is, we enforce
$k_{\text{min}}<|{\mathchoice{{\hbox{\BFd k}}}{{\hbox{\BFt k}}}{{\hbox{\BFs
k}}}{{\hbox{\BFss k}}}}_{i}|<k_{\text{max}}$ for each wavevector
${\mathchoice{{\hbox{\BFd k}}}{{\hbox{\BFt k}}}{{\hbox{\BFs k}}}{{\hbox{\BFss
k}}}}_{i}$, along with the usual constraint that ${\mathchoice{{\hbox{\BFd
k}}}{{\hbox{\BFt k}}}{{\hbox{\BFs k}}}{{\hbox{\BFss
k}}}}_{1}+{\mathchoice{{\hbox{\BFd k}}}{{\hbox{\BFt k}}}{{\hbox{\BFs
k}}}{{\hbox{\BFss k}}}}_{2}+{\mathchoice{{\hbox{\BFd k}}}{{\hbox{\BFt
k}}}{{\hbox{\BFs k}}}{{\hbox{\BFss k}}}}_{3}=\bm{0}$. Owing to homogeneity and
isotropy, all bispectra have three arguments: $|{\mathchoice{{\hbox{\BFd
k}}}{{\hbox{\BFt k}}}{{\hbox{\BFs k}}}{{\hbox{\BFss k}}}}_{1}|$,
$|{\mathchoice{{\hbox{\BFd k}}}{{\hbox{\BFt k}}}{{\hbox{\BFs k}}}{{\hbox{\BFss
k}}}}_{2}|$, $|{\mathchoice{{\hbox{\BFd k}}}{{\hbox{\BFt k}}}{{\hbox{\BFs
k}}}{{\hbox{\BFss k}}}}_{3}|$. Our bins are cubes of side length $4\,k_{f}$ in
Fourier space. This characterization captures nearly all the available modes.
As discussed in § 2, we also enforce that $k_{1}\geq k_{2}\geq k_{3}$. In
estimating bispectra, we loop through all of the simulated Fourier modes that
lie within this wavevector range. We determine $B_{21,\delta,\delta}$ and
$P_{21,\delta}$ from the same set of Fourier modes to reduce sample variance
in our estimates of $Q_{21,\delta,\delta}$. In general, one would like to make
$k_{\text{min}}$ and $k_{\text{max}}$ as small and large as possible,
respectively. However, our estimates at small $k$ are noisy and the
perturbative formula becomes less accurate at high $k$. We chose values of
$k_{\text{min}}$ and $k_{\text{max}}$ throughout that give the best results,
but also show results for different values of $k_{\text{min}}$ and
$k_{\text{max}}$ in Appendix A.

##  4\. Results

Figure 3.— Reduced bispectrum of the density field. Here we fix
$k_{2}=2k_{1}=0.2\,\text{Mpc}^{-1}$ in the Upper panel and
$k_{2}=2k_{1}=0.4\,\text{Mpc}^{-1}$ in the Bottom panel. We plot the results
as a function of $\theta_{12}$, the angle between the two wavevectors. The
solid line is the prediction of second-order perturbation theory, as given in
Eq. 5. Each colored point shows the average reduced bispectrum from ten
different simulation realizations. The error bars give estimates of the error
on the average, as determined from the scatter across the ten simulation
realizations. Note that the poorer sampling for nearly collinear triangle
configurations leads to noisier estimates near $\theta_{12}/\pi\approx
0\text{,}1$. The dashed line is the perturbative $Q_{\delta,\delta,\delta}$
but with a constant added to minimize the $\chi^{2}$ difference with the
$z=7.37$ measurements. This indicates that, for the most part, the
perturbative $Q_{\delta,\delta,\delta}$ differs from the simulation results
only by an overall constant.

First we consider the reduced auto-bispectrum of the simulated matter density
field, $Q_{\delta,\delta,\delta}(\theta_{12})$, at the redshifts of interest
in Fig. 3. As discussed in § 3.1, we average over ten independent simulation
realizations. We show the simulated reduced density bispectrum for triangles
with $k_{1}=0.1\,\text{Mpc}^{-1}$ and $k_{2}=2k_{1}$ and a case where the
wavenumbers are twice as large. In the larger-scale case, the simulated
bispectra are somewhat noisy while the smaller-scale case provides a smoother
bispectrum estimate, because of the larger number of sampled modes at smaller
scales. The simulation measurements are generally similar to the predictions
of second order perturbation theory (Eq. 7) for the large scales shown and the
high redshifts considered, although there are noticable differences. The
departures from the predictions of second-order perturbation theory evident in
the figure reflect non-linearities that are not fully captured at second-order
in $\delta$. In order to partly account for these non-linear density
fluctuations, we subsequently use the simulated density auto-bispectrum,
$Q_{\delta,\delta,\delta}$, as input to the perturbative formulae (Eqs. 14 and
15). Throughout, we refer to this prediction loosely as the “second-order
perturbation theory prediction” although strictly speaking it accounts for the
impact of additional non-linearities in the form of
$Q_{\delta,\delta,\delta}$. This also serves to reduce the impact of sample
variance in our comparisons between simulations and perturbation theory
estimates.

That being said, the departures from non-linearity manifest mainly as an
overall constant. We illustrate this by showing in Fig. 3 (dashed line) the
perturbative $Q_{\delta,\delta,\delta}$ but with a constant added to minimize
the $\chi^{2}$ statistic at $z=7.37$. The agreement between the dashed line
and the points at $k_{2}=2k_{1}=0.4\,\text{Mpc}^{-1}$ indicates that the non-
linearities will propagate into $C_{21,\text{X},\text{X}}$ but not into
$\langle T_{21}\rangle b_{21}$. While this is not the case for the noisier
$k_{2}=2k_{1}=0.2\,\text{Mpc}^{-1}$ panel, we nonetheless proceed with these
caveats in mind.

Figure 4.— Evolution of various reduced bispectra with redshift/ionization
fraction. Left: the reduced 21 cm auto-bispectrum. Right: the reduced 21 cm-
density-density cross-bispectrum. In each panel, the solid lines show the best
fit models of the form $Q_{\delta,\delta,\delta}/\langle T_{21}\rangle
b_{21}+C$, using the simulated values of $Q_{\delta,\delta,\delta}$. We
restrict ourself to modes that satisfy $k_{1}=0.1\,\text{Mpc}^{-1}$ and
$k_{2}=0.2\,\text{Mpc}^{-1}$, varying the angle $\theta_{12}$ between the two
wavevectors. Note that the best fit values of $\langle T_{21}\rangle b_{21}$
and $C$ here are obtained by fitting only to this restricted range of triangle
configurations and differ somewhat from the best fits obtained subsequently,
which match to a broader range of triangles. The perturbative formula provides
a better match to $Q_{21,\delta,\delta}$ than $Q_{21,21,21}$ as we discuss in
the text. Note that since we use the simulated $Q_{\delta,\delta,\delta}$ in
our fits, the curves share some of the noise in our simulation estimates.

We next consider the reduced 21 cm auto-bispectrum as a function of ionization
fraction, along with the 21 cm-$\delta$-$\delta$ cross-bispectrum.555The 21
cm-[CII]-[CII] cross-bispectrum is ultimately the quantity of interest, but as
discussed in § 3.3, we are unable to capture small mass halos and so we use
the 21 cm-$\delta$-$\delta$ bispectrum as a proxy for the 21 cm-[CII]-[CII]
bispectrum. If the second-order predictions in § 2 are correct, then we expect
that all three of these bispectra differ only by a redshift-dependent
constant. Our goals here are two-fold: first, we want to ascertain whether
these quantities evolve strongly with $\langle x_{i}\rangle$, as expected if
their evolution is driven, in part, by the “rise and fall” behavior
illustrated in Fig. 1. Second, we aim to see if the simulated bispectra follow
the perturbative expectation that
$Q_{21,\text{X},\text{X}}=Q_{\delta,\delta,\delta}/\langle T_{21}\rangle
b_{21}+C_{21,\text{X},\text{X}}$.

Fig. 4 shows that the various bispectra evolve strongly with ionization
fraction. This is encouraging for our ultimate goal of using these
measurements to extract information about the ionization history. The overall
value and evolution of the reduced bispectra vary somewhat depending on which
bispectrum is calculated, but this is expected since the values of
$C_{21,21,21}$ and $C_{21,\delta,\delta}$ should differ by a factor of 3.
Nevertheless, the overall shape of the two bispectra are similar, and the
bispectra share the same ordering – from smallest to largest $Q$ – with
increasing ionization fraction.

In order to test the perturbative formula, we determine the best fit $\langle
T_{21}\rangle b_{21}$ and $C_{21,\text{X},\text{X}}$ for
$Q_{21,\text{X},\text{X}}$. To carry out these fits we use the variance
estimated from the scatter across our ten simulation realizations. We neglect
the covariance between different wavevector bins, since our estimates of the
off-diagonal elements of the covariance matrix are too noisy to be reliable.
We expect this to have relatively little impact on the best fit parameters,
but it does prevent us from assessing the overall goodness of fit of the
perturbative formula since the effective number of degrees of freedom are
uncertain. The inability to measure covariance also prevents us from
accurately computing error bars on the inferred $\langle T_{21}\rangle b_{21}$
and $C_{21,\text{X},\text{X}}$. We defer this calculation to future studies.

##  5\. Checking the Perturbative Formulas

Figure 5.— Consistency between different estimates of $\langle T_{21}\rangle
b_{21}$. The black line show the linear 21 cm bias factor, $\langle
T_{21}\rangle b_{21}$, inferred from the cross-power spectrum, $\langle
T_{21}\rangle b_{21}=P_{21,\delta}/P_{\delta,\delta}$ at
$k=0.05\,\text{Mpc}^{-1}$, as a function of $\langle x_{i}\rangle$. These are
compared with the linear bias inferred from $Q_{21,21,21}$ (red) and
$Q_{21,\delta,\delta}$ (blue).

An interesting test of the perturbative framework is to compare the values of
$\langle T_{21}\rangle b_{21}$ and $C_{21,\text{X},\text{X}}$ inferred from a
broad range of simulated two and three-point statistics. First, we compare the
values of $\langle T_{21}\rangle b_{21}$ extracted from our fits to
$Q_{21,\text{X},\text{X}}$ with those from the 21 cm-density cross-power
spectrum. The results of this comparison are shown for a range of ionization
fractions in Fig. 5. In general, all three approaches return similar values of
$\langle T_{21}\rangle b_{21}$ and share the same qualitative evolution with
ionization fraction. However, the $\langle T_{21}\rangle b_{21}$ estimates
from the cross-power spectrum and the cross-bispectrum agree better with each
other than with the auto-bispectrum inferences. This likely results because
the 21 cm auto-bispectrum involves three significantly non-Gaussian fields (at
least on scales smaller than the size of the ionized regions) and so our
perturbative expansion is presumably less accurate for the auto-bispectrum. We
find that by using $Q_{21,\delta,\delta}$, we are able to receover $\langle
T_{21}\rangle b_{21}$ to within $5\%$ for $\langle x_{i}\rangle<0.5$, and
$10\%$ for $\langle x_{i}\rangle>0.5$, except at the lowest redshift we tested
($z=6$). On the other hand, $Q_{21,21,21}$ is only able to recover $\langle
T_{21}\rangle b_{21}$ to within $\sim 20\%$ accuracy for all of the EoR,
although the precise number quoted here may be subject to noise in our
$Q_{21,21,21}$ estimates.

The agreement between the cross-power spectrum and cross-bispectrum estimates
is best at $\langle x_{i}\rangle\lesssim 0.6$, at which point they separate a
bit, before briefly coming together near $\langle x_{i}\rangle\sim 0.8$; they
then split-off slightly again towards the tail-end of reionization. The
results in Fig. 5 adopt triangles with wavevectors between
$k_{\text{min}}=0.05\,\text{Mpc}^{-1}$ and
$k_{\text{max}}=0.4\,\text{Mpc}^{-1}$. In Appendix A (Fig. 10) we investigate
other choices for the range of wavevectors used in our fitting procedure.
There we show that the values of $\langle T_{21}\rangle b_{21}$ inferred from
$Q_{21,\delta,\delta}$ are insensitive to the precise choices of
$k_{\text{min}}$ and $k_{\text{max}}$, while the results from $Q_{21,21,21}$
are noisier and show more sensitivity to the choice of scales.

Figure 6.—  Evolution of the constant offset terms, $C_{21,X,X}$ with
ionization fraction. Note that it appears that $C_{21,21,21}\sim
C_{21,\delta,\delta}$ for $\langle x_{i}\rangle\sim 0.7$. This occurs,
however, because both constant offset terms are very close to zero at this
ionization fraction. See Fig. 7.

We can also test whether the constant offset terms agree with the formulas
from perturbation theory (Eqs. 16-18). Before considering this test, we note
that the behavior of the constant offset term with redshift may offer an
additional handle on the ionization history since e.g. $C_{21,21,21}$ is a
measurable quantity. Indeed, Fig. 6 shows that the constant offsets evolve
fairly strongly with redshift/ionization fraction, and that the different
offset terms evolve in qualitatively similar ways with average ionization
fraction. However, note that $C_{21,\text{CII},\text{CII}}$, which is most
robust to foreground contamination, depends on $b_{\text{CII}}$ and
$b_{\text{CII}}^{(2)}$; that is, this quantity depends on both the [CII] and
the 21 cm biasing and is therefore more challenging to interpret.

The perturbative prediction (Eqs. 16 and 17) is that the constant offsets are
related by $C_{21,21,21}=3C_{21,\delta,\delta}$. Since we extracted each of
these terms independently, we are able to test this statement without
extracting $\langle T_{21}\rangle b_{21}^{(2)}$. We plot the ratio
$C_{21,21,21}/3C_{21,\delta,\delta}$ in Fig. 7. It is encouraging that for the
much of the reionization process this ratio lies within $10-20\%$ of unity.
The departures shown in the figure arise: i) at very early times (near
$\langle x_{i}\rangle\sim 0.1$) when the correlation coefficient reverses sign
(see §3.2), ii) near $\langle x_{i}\rangle\sim 0.7$ when each of the constant
offset terms happen to be close to zero and sample variance likely impacts our
estimate of this ratio, and iii) near the end of reionization when the ionized
regions are large and the perturbative framework is expected to break down.

Figure 7.—  Consistency check on the constant offset terms. If second-order
perturbation theory is accurate, we expect
$C_{21,21,21}/3C_{21,\delta,\delta}=1$. We find that this ratio is within
$20\%$ of unity during the bulk of reionization. The dip at $\langle
x_{i}\rangle\sim 0.7$ results because
$C_{21,21,21}\text{,}C_{21,\delta,\delta}\sim 0$: at this stage of
reionization, more simulation realizations would be required to accurately
estimate the ratio of these small values.

For completeness, we note that with more detailed [CII] emission models, one
could carry out further tests of the perturbative framework. Specifically, the
perturbative formulas demand specific relationships between the constant
offset terms for the 21 cm-[CII]-[CII] cross-bispectrum, the 21
cm-$\delta$-$\delta$ cross-bispectrum, the 21 cm auto-bispectrum, and the
[CII] auto-bispectrum. Using Eqs. 15 \- 18, and analogous formulas for the
[CII] intensity field, we have that,

$\frac{3}{2}\frac{C_{21,\text{CII},\text{CII}}-C_{21,\delta,\delta}}{C_{\text{CII},\text{CII},\text{CII}}}=\frac{b_{\text{CII}}}{b_{21}}$
(19)

and, similarly,

$\frac{3}{2}\frac{C_{21,\text{CII},\text{CII}}-\frac{1}{3}C_{21,21,21}}{C_{\text{CII},\text{CII},\text{CII}}}=\frac{b_{\text{CII}}}{b_{21}}\text{.}$
(20)

We comment in passing that another possibility is to extract $\langle
I_{\text{CII}}\rangle b_{\text{CII}}$ from the analogous cross-bispectrum,
$Q_{\text{CII},21,21}$. If line-intensity mapping experiments are performed
using additional emission lines such as Ly-$\alpha$ or using rotational
transitions from CO molecules, the cross-bispectra methodology may allow
further inferences and consistency checks.

##  6\. Detectability

We now discuss the prospects of measuring the 21 cm-[CII]-[CII] cross-
bispectrum from future surveys. Here our aim is to obtain rough estimates of
the joint sky coverage and noise power spectra required for upcoming [CII] and
21 cm surveys to detect this signal.

###  6.1. [CII] Statistics

As discussed in § 3.3, our current modeling assumes that the reduced cross-
bispectrum $Q_{21,\text{CII},\text{CII}}$ tracks $Q_{21,\delta,\delta}$ apart
for an overall constant. In order to calculate the variance of
$Q_{21,\text{CII},\text{CII}}$, we require a model for the [CII] emission
power spectrum and the cross-power spectrum with the 21 cm field. For this
purpose, we assume a linear biasing model for the [CII] emission fluctuations:

$\begin{split}P_{\text{CII},\text{CII}}&=\big{(}b_{\text{CII}}\langle
I_{\text{CII}}\rangle\big{)}^{2}P_{\delta,\delta}\\\
P_{21,\text{CII}}&=b_{\text{CII}}\langle I_{\text{CII}}\rangle
P_{21,\delta}\text{.}\end{split}$ (21)

Here $b_{\text{CII}}$ is the luminosity-weighted bias of the halos that host
[CII]-emitting galaxies. For simplicity, we suppose that $b_{\text{CII}}=3$ at
all redshifts of interest; our results are relatively insensitive to this
choice. We follow Lidz & Taylor (2016) (which is based on Lidz et al. (2011);
Pullen et al. (2014)) in calculating $\langle I_{\text{CII}}\rangle$.
Specifically, we assume a one-to-one mapping between star formation rate, SFR,
and [CII] luminosity. The average specific intensity of the [CII] emission is:

$\langle
I_{\text{CII}}\rangle(z)=\frac{\epsilon_{\text{CII}}}{4\pi\nu_{\text{rest}\text{,}\text{CII}}}\frac{c}{H(z)}\text{,}$
(22)

where $\epsilon_{\text{CII}}$ denotes the average co-moving emissivity
(luminosity density) in the line, and $\nu_{\text{rest}\text{,}\text{CII}}$ is
the rest frame frequency of the [CII] line. We have further assumed a delta
function [CII] line emission profile.

The one-to-one relation between [CII] luminosity and star-formation rate (SFR)
we adopt is:

$L=L_{0}\frac{\text{SFR}}{1M_{\sun}\text{yr}^{-1}}\text{,}$ (23)

with $L_{0}=6\times 10^{6}L_{\sun}$ (Righi et al., 2008; Visbal & Loeb, 2010).
We neglect any redshift evolution in $L_{0}$ (see e.g. Lidz & Taylor 2016 for
a discussion). We assume that the abundance of galaxies as a function of their
SFR obeys a Schechter function form (Schechter, 1976):

$\phi(\text{SFR})\,\text{d}\text{SFR}=\phi_{*}\left(\frac{\text{SFR}}{\text{SFR}_{*}}\right)^{\alpha}\exp\left[-\frac{\text{SFR}}{\text{SFR}_{*}}\right]\frac{\text{d}\text{SFR}}{\text{SFR}_{*}}\text{.}$
(24)

Here $\alpha$ denotes the faint-end slope, while $\phi_{*}$ and
$\text{SFR}_{*}$ are characteristic number densities and SFRs. The co-moving
[CII] emissivity is then:

$\epsilon_{\text{CII}}=\phi_{*}L_{0}\frac{\text{SFR}_{*}}{1\,M_{\sun}\text{yr}^{-1}}\Gamma(2+\alpha)\text{.}$
(25)

We use the Schechter SFR function parameters in Smit et al. (2012), determined
using dust-corrected ultraviolet luminosity functions at redshifts of
$z=2.2\text{, }3.8\text{, }5.0\text{, }5.9\text{, }6.8$. We find that the
following functional form provides a reasonable approximation to the redshift
evolution of the Schechter parameters:

$\begin{split}\alpha&=-1.96\\\
\text{SFR}_{*}(z)&=134.39\,z^{-1.18}\,$\mathrm{M_{\sun}}\text{\,}{\mathrm{yr}}^{-1}$\\\
\phi_{*}(z)&=0.00358\,z^{-0.85}\,${\mathrm{Mpc}}^{-3}$\text{.}\end{split}$
(26)

This procedure allows us to compute $\langle I_{\text{CII}}\rangle$ at each
relevant redshift. We find: $\langle I_{\text{CII}}\rangle=3.1\times
10^{2}\,\text{Jy}/\text{str}\text{,}\,7.1\times
10^{2}\,\text{Jy}/\text{str}\text{,}\,1.1\times 10^{3}\,\text{Jy}/\text{str}$
for $z=9.41\text{,}\,7.37\text{,}\,6.43$, respectively.

###  6.2. SNR Formula

In order to estimate the signal-to-noise ratio (SNR), we require a formula for
the cross-bispectrum variance, $\mathrm{Var}[Q_{21,\text{CII},\text{CII}}]$.
In calculating the variance, we suppose that the [CII] and 21 cm fluctuation
fields obey Gaussian statistics. The variance of the cross-bispectrum (for
$k_{1}\neq k_{2}\neq k_{3}$) is (see e.g. Greig et al. 2013; Joachimi et al.
2009 for related calculations):

$\begin{split}&\mathrm{Var}[B_{21,\text{CII},\text{CII}}+\text{2 perm.}]=\\\
&\frac{V_{s}}{N_{t}}\big{(}P_{\text{tot,21}}(k_{1})P_{\text{tot,CII}}(k_{2})P_{\text{tot,CII}}(k_{3})+\text{2
perm.}\big{)}\text{,}\end{split}$ (27)

where $P_{\text{tot}}=P+N$ is the total signal plus detector noise auto-power
spectrum for each of the 21 cm and [CII] emission fields. The noise power
spectrum for the 21 cm survey is subsequently denoted by $N_{21}$ and the
[CII] noise power by $N_{\text{CII}}$. We neglect contributions to the
variance from residual foregrounds. Later we will add in the permutations on
$k_{1}$, $k_{2}$, and $k_{3}$. Note that if any of $k_{1}$, $k_{2}$, or
$k_{3}$ are equal, there are extra terms (which we write down explicitly in
Appendix B). These extra terms are what give rise to the symmetry factor,
$s_{B}$, in Greig et al. (2013) and related works.

Now, we use propagation of errors to recover $\mathrm{Var}[Q]$,

$\begin{split}\mathrm{Var}[Q]&=\left(\frac{\partial Q}{\partial
B}\right)^{2}\mathrm{Var}[B]\\\ &+\sum_{i=1}^{3}\left(\frac{\partial
Q}{\partial P(k_{i})}\right)^{2}\mathrm{Var}[P(k_{i})]\\\
&+\sum_{j=1}^{3}\frac{\partial Q}{\partial P(k_{j})}\frac{\partial Q}{\partial
B}\mathrm{Cov}[B,P(k_{j})]\text{,}\end{split}$ (28)

where the summations $i$ and $j$ span the three wavevector arguments, $k_{1}$,
$k_{2}$, and $k_{3}$. Under the Gaussian approximation adopted here,
$\mathrm{Cov}[B,P(k_{i})]$ vanishes because it is a five-point function with
vanishing unconnected pieces. Furthermore, we suppose that the term involving
$\mathrm{Var}[P(k_{i})]$ is sub-dominant and neglect it in what follows.
Finally, accounting for wavevector permutations, we arrive at our formula for
the variance on $Q$:

$\mathrm{Var}[Q]=\frac{V_{s}}{N_{t}}\frac{\left(P_{\text{tot,21}}(k_{1})P_{\text{tot,CII}}(k_{2})P_{\text{tot,CII}}(k_{3})+\text{2
perm.}\right)}{9\left(P_{21,\text{CII}}(k_{1})P_{21,\text{CII}}(k_{2})+\text{2
perm.}\right)^{2}}\text{,}$ (29)

where $V_{s}$ is the survey volume, and $N_{t}$ is the number of triangles,
which we will derive shortly. Note that the $1/9$ factor comes from the $1/3$
in our definition of $Q_{21,\text{CII},\text{CII}}$ in Eq. 13

In order to estimate the expected error bars on the parameters of interest,
specifically $\langle T_{21}\rangle b_{21}$ and
$C_{21,\text{CII},\text{CII}}$, we use the Fisher matrix formalism. The Fisher
matrix is calculated as (e.g. Greig et al. 2013):

$F_{ij}=\sum_{k_{1}\leq k_{2}\leq
k_{3}}\frac{1}{\mathrm{Var}[Q](A,N_{21},N_{\text{CII}})}\frac{\partial
Q}{\partial\theta_{i}}\frac{\partial Q}{\partial\theta_{j}}\text{,}$ (30)

where $\theta_{i}$ denotes our two model parameters ($\langle T_{21}\rangle
b_{21}$ and $C_{21,\text{CII},\text{CII}}$.)

###  6.3. Survey Parameters

We adopt a simple, yet flexible, description for the upcoming 21 cm and [CII]
surveys. First, we assume that the noise power spectrum for each survey is
well-approximated by isotropic, white-noise. In detail, this is an imperfect
approximation, especially for the interferometric 21 cm observations in which
case the noise power is a strong function of $k_{\perp}$; in reality, the 21
cm noise power is sensitive to the precise distribution of interferometric
baselines. Our estimates here should be refined in future work. In our
fiducial model, we match the amplitude of the white-noise power spectrum in
each of the 21 cm and [CII] fields to their signal power at
$k=0.05\,\text{Mpc}^{-1}$, i.e. we set $N_{21}=P_{21,21}$ at this wavenumber,
and similarly for $N_{\text{CII}}$. This amounts to assuming that larger
spatial scales are sample-variance dominated, while smaller-scales are limited
by detector noise. For the survey volume, we assume that a common survey area,
$A_{\text{survey}}$, is shared between the 21 cm and [CII] surveys, and adopt
a redshift bin width of $\Delta z=0.3$. Our fiducial model takes
$A_{\text{survey}}=50\,\text{deg}^{2}$.

The number of triangles (with wavevectors within a given range) that fit into
the overlap region of the two surveys may be calculated as:

$N_{t}=\frac{V_{B}}{V_{\text{fund}}^{2}}\text{.}$ (31)

Here $V_{B}$ is the 6D Fourier volume satisfying
$\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3}=\bm{0}$ of the 9D Fourier space spanned by
$(\bm{k}_{1}\text{,}\bm{k}_{2}\text{,}\bm{k}_{3})$. More specifically,

$V_{B}\equiv\int_{k_{1}}\text{d}^{3}q_{1}\int_{k_{2}}\text{d}^{3}q_{2}\int_{k_{3}}\text{d}^{3}q_{3}\,\delta_{D}(\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3})\text{,}$
(32)

where $\int_{k_{i}}$ refers to the integral over a wavevector bin centered
around $|\bm{q_{i}}|=|\bm{k_{i}}|$. Assuming that $k_{1}\geq k_{2}\geq k_{3}$,
one can show that (e.g. Greig et al. (2013)),

$V_{B}\simeq 8\pi^{2}k_{1}k_{2}k_{3}\,(\Delta
k)^{3}\,\theta(k_{1},k_{2},k_{3})\text{,}$ (33)

where,

$\theta(k_{1},k_{2},k_{3})=\begin{cases}\frac{1}{2}&\text{if
$k_{i}=k_{j}+k_{k}$}\\\ 1&\text{if $k_{i}\neq k_{j}+k_{k}$,}\end{cases}$ (34)

is a symmetry factor ensuring that equivalent modes are not double counted. We
adopt a bin size of $\Delta k=0.03\,\text{Mpc}^{-1}$. The quantity
$V_{\text{fund}}$ in Eq. 31 denotes the $k$-space volume spanned by the
fundamental modes of the joint survey regions. For this, we assume that the
joint survey area is a square, giving:

$V_{\text{fund}}=\frac{(2\pi)^{3}}{L_{\perp}^{2}L_{\parallel}}$ (35)

where $L_{\perp}$ and $L_{\parallel}$ are the co-moving extent of the survey
volume perpendicular, parallel to the line of sight. We calculate $L_{\perp}$
at the central redshift of each redshift bin.

We then proceed to calculate the expected SNR on $\langle T_{21}\rangle
b_{21}$, marginalizing over $C_{21,21,21}$, using the inverse Fisher matrix
(Eq. 30). We perform these calculations at each of three different redshifts,
while varying the survey parameters, $N_{\text{CII}}$, $A_{\text{survey}}$,
and $N_{21}$ around our fiducial values to explore how our forecasts depend on
these inputs.

Figure 8.— Estimated SNR of $Q_{21,\text{CII},\text{CII}}$ as a function of
(Left-hand panel) $N_{\text{CII}}$, (Middle panel) $A_{\text{survey}}$, and
(Right-hand panel) $N_{\text{21}}$. When not varied, we assumed that
$N_{\text{CII}}=P_{\text{CII},\text{CII}}$ at $k=0.05\,\text{Mpc}^{-1}$,
$A_{\text{survey}}=50\,\text{deg}^{2}$ and $N_{21}=P_{21,21}$ at
$k=0.05\,\text{Mpc}^{-1}$. Fiducial values are given by vertical dashed lines
in each plot – for $N_{21}$ and $N_{\text{CII}}$ these are redshift dependent.

The results of these calculations are shown in Fig. 8. These estimates are
generally encouraging: for example, in order to achieve a $10-\sigma$
detection (i.e. $\text{SNR}=10$) of the 21 cm-[CII]-[CII] bispectrum near the
mid-point of reionization (at $z=7.37$ in this model) a $50\,\text{deg}^{2}$
survey with $N_{\text{CII}}\lesssim 1.6\times
10^{8}\text{\,}(\text{Jy}/\text{sr})^{2}\,\text{Mpc}^{3}$ is required. This is
more demanding than the noise expected for the “Stage II” [CII] survey
described in Silva et al. (2015) which has a noise power spectrum of
$N_{\text{CII}}\sim 2.5\times
10^{9}\,(\text{Jy}/\text{sr})^{2}\,\text{Mpc}^{3}$ at $z=7$ (Lidz & Taylor,
2016). Nonetheless, this improved sensitivity may be achievable with future
improvements in detector technology, for example. As quantified in the figure,
the requirements are somewhat less stringent at $z=6.43$ and somewhat more so
at $z=9.41$. Note that our model assumes that the $L_{\text{CII}}$-SFR
relation is independent of redshift. This assumption is most suspect for this
highest redshift bin, since the galaxies at this early time may have low
metallicities and little [CII] emission. In any case, none of the currently
planned [CII] emission surveys target such high redshifts (Kovetz et al.,
2017). The requirements on the 21 cm thermal noise appear less severe, with
the SNR($\langle T_{21}\rangle b_{21}$) saturating for $N_{21}\lesssim
10^{5}\,\text{mK}^{2}\,\text{Mpc}^{3}$. The results are less sensitive to the
21 cm noise power spectrum than the [CII] noise both because our cross-
bispectrum involves only one 21 cm field (yet two [CII] fields) and since the
21 cm power spectrum becomes sample-variance limited at the wavenumbers of
interest. The latter fact explains the saturation in the SNR with decreasing
$N_{21}$. The required $N_{21}$ seems feasible since HERA-350 will image some
large-scale modes (DeBoer et al., 2017), implying that sample-variance limited
sensitivity will be achieved on large scales.

Another approach for improving the SNR($\langle T_{21}\rangle b_{21}$) is to
increase the joint sky-coverage. For example, one possibility is to measure
the 21 cm-Ly-$\alpha$-Ly-$\alpha$ bispectrum (rather than the 21
cm-[CII]-[CII] bispectrum considered here). This may be feasible with the
planned all-sky survey SPHEREX (Doré et al., 2018), in which case one could
match the entire $\sim 1440\,\text{deg}^{2}$ coverage of HERA-350 (DeBoer et
al., 2017). The figure illustrates that our forecasts improve dramatically for
wider survey areas.

To further explore the potential reach of our method, we consider more
futuristic surveys and plot $1$ and $2-\sigma$ contours in the $\langle
T_{21}\rangle b_{21}$-$C_{21,\text{CII},\text{CII}}$ plane at different
redshifts in Fig. 9. In this case, we set the joint survey area to be
$A=1000\,\text{deg}^{2}$ and we match $N_{21}$ and $N_{\text{CII}}$ to
$P_{21}$ and $P_{\text{CII}}$ at $k=0.1\,\text{Mpc}^{-1}$. This corresponds to
noise values of $N_{21}=1.7\times 10^{5}\text{,}\,2.6\times
10^{5}\text{,}\,7.6\times 10^{5}\,\text{mK}^{2}\,\text{Mpc}^{3}$, and
$N_{\text{CII}}=3.4\times 10^{9}\text{,}\,1.1\times 10^{9}\text{,}\,1.4\times
10^{8}\,(\text{Jy}/\text{sr})^{2}\,\text{Mpc}^{3}$ at
$z=6.43\text{,}\,7.37\text{,}\,9.41$. This survey should indeed enable a sharp
test of the “rise and fall” behavior from the cross-bispectrum. The contours
also help illustrate the level of degeneracy between $\langle T_{21}\rangle
b_{21}$ and $C_{21,\text{CII},\text{CII}}$; these parameters are less
covariant during the middle of reionization (when $\langle T_{21}\rangle
b_{21}$ is large). This results because $Q_{21,\text{CII},\text{CII}}$ has a
flatter dependence on $\theta_{12}$ at this stage of reionization, which makes
it easier to estimate $C_{21,\text{CII},\text{CII}}$.

Figure 9.— Forecasted constraints in the $\langle T_{21}\rangle
b_{21}$-$C_{21,\text{CII},\text{CII}}$ plane at different redshifts. The
solid/dashed ellipses show the expected $1-\sigma$/$2-\sigma$ contours at
different redshifts, with each contour centered on the values of $\langle
T_{21}\rangle b_{21}$ and $C_{21,\text{CII},\text{CII}}$ in our simulated
models. These forecasts assume an ambitious, futuristic survey, with a sky
coverage of $A=1000\,\text{deg}^{2}$ and low detector noise with $N_{21}$ and
$N_{\text{CII}}$ matching $P_{21}$ and $P_{\text{CII}}$ at the scale of
$k=0.1\,\text{Mpc}^{-1}$ (see text).

##  7\. Conclusions

Our aim in this work has been to help circumvent two of the primary challenges
associated with 21 cm fluctuation measurements of the EoR. The two concerns
are that: residual foregrounds may produce errors in the inferred 21 cm auto-
spectrum (or lead to a spurious detection), and that the 21 cm auto-spectrum
is itself difficult to interpret given our imperfect models of the EoR. The
first issue may, in part, be side-stepped by cross-correlating 21 cm
fluctuation measurements with an additional tracer of the high redshift
universe, such as line-intensity mapping data cubes in the [CII] emission
line. The second concern is, in part, alleviated by identifying statistics
that are amenable to an analytic treatment.

Here, we proposed that the 21 cm-[CII]-[CII] cross-bispectrum may help in both
regards. First, residual foregrounds only impact the variance of the cross-
bispectrum and not the ensemble average (modulo foregrounds that are common to
the two surveys.) Second, we showed that the configuration dependence of the
21 cm-$\delta$-$\delta$ cross-bispectrum (a proxy for the 21 cm-[CII]-[CII]
cross-bispectrum which is more challenging to simulate) is fairly well
described by second-order perturbation theory, especially provided one
considers large spatial scales before the ionization fraction reaches $\langle
x_{i}\rangle\sim 0.7$. Further, this may be used to extract the redshift
evolution of $\langle T_{21}\rangle b_{21}$ without resorting to reionization
simulations. While such simulations are nevertheless required to understand
the full implications of the $\langle T_{21}\rangle b_{21}$ measurements for
the reionization history and the properties of the ionizing sources, we
believe it is still valuable to consider upcoming 21 cm measurements in the
context of cosmological perturbation theory. For example, the redshift
evolution of $\langle T_{21}\rangle b_{21}$ inferred (on large scales) from
the 21 cm auto-power spectrum, the 21 cm auto-bispectrum, and the 21
cm-[CII]-[CII] cross-bispectrum should each show the “rise and fall”
signature. More generally, demonstrating the consistency of inferences from
two and three-point statistics should help in establishing the robustness of
initial 21 cm detections. Unsurprisingly, the perturbative description is
imperfect and future simulation efforts may help in exploiting information on
small spatial scales, especially during the late stages of reionization.

We made rough estimates of the survey specifications required to detect the 21
cm-[CII]-[CII] cross-bispectrum. Assuming a joint survey area of
$A_{\text{survey}}=50\,\text{deg}^{2}$, we found that a detection generally
requires a more sensitive [CII] survey than the Stage-II [CII] emission survey
described in Silva et al. (2015). It would be interesting to refine our
analysis by considering a more detailed treatment of the noise in each survey,
by quantifying the impact of residual foregrounds on the cross-bispectrum
variance, and by exploring to what extent foregrounds may be shared by the two
surveys. Since the perturbative description is accurate on large spatial
scales, this approach is most valuable if measurements can be made robustly at
small wavenumbers ($k\lesssim 0.05-0.1\,\text{Mpc}^{-1}$).

In terms of modeling, we showed that the L-PICOLA code may be used in
conjunction with 21cmFAST (replacing the Zel’dovich generated density field
used in 21cmFAST by default) to improve the treatment of three-point
statistics in EoR calculations, at only modest additional computational
expense. In the future, it would be interesting to include a more detailed
model for the [CII]-emission fluctuations: it is challenging to capture the
large spatial scales of interest for our bispectrum calculations, while
simultaneously resolving the small mass halos hosting [CII]-emitting galaxies.
We considered only a single representative reionization model: we expect the
perturbative description to work better if the ionized regions are smaller
than in this model, while it will perform less well in scenarios with larger
bubbles. While we considered only the 21 cm and [CII] emission lines in this
work, the same methodology may be applied to other emission lines and it will
be interesting to consider these prospects. In summary, higher-order
statistics should help in extracting key information about the EoR from
upcoming surveys.

## Acknowledgements

We thank the anonymous referee for providing helpful comments. We would like
to thank Matthew McQuinn for helpful comments on a draft of this paper. AB
would like to thank Congzhou M. Sha for helpful comments on code used in this
work, and Todd Phillips for helpful discussions. AB was supported by the Roy &
Diana Vagelos Program in the Molecular Life Sciences and the Roy & Diana
Vagelos Challenge Award.

## References

  * Becker et al. (2015) Becker, G. D., Bolton, J. S., & Lidz, A. 2015, PASA, 32, e045
  * Becker et al. (2001) Becker, R. H., Fan, X., White, R. L., et al. 2001, AJ, 122, 2850
  * Bernardeau et al. (2002) Bernardeau, F., Colombi, S., Gaztanaga, E., & Scoccimarro, R. 2002, Phys. Rept., 367, 1
  * Bharadwaj & Pandey (2005) Bharadwaj, S., & Pandey, S. K. 2005, Mon. Not. Roy. Astron. Soc., 358, 968
  * Bond et al. (1991) Bond, J. R., Cole, S., Efstathiou, G., & Kaiser, N. 1991, Astrophys. J., 379, 440
  * Chung et al. (2017) Chung, D. T., Li, T. Y., Viero, M. P., Church, S. E., & Wechsler, R. H. 2017, Astrophys. J., 846, 60
  * Crites et al. (2014) Crites, A. T., Bock, J. J., Bradford, C. M., et al. 2014, in Proc. SPIE, Vol. 9153, Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy VII, 91531W
  * Crocce et al. (2006) Crocce, M., Pueblas, S., & Scoccimarro, R. 2006, MNRAS, 373, 369
  * DeBoer et al. (2017) DeBoer, D. R., Parsons, A. R., Aguirre, J. E., et al. 2017, PASP, 129, 045001
  * Dekel & Lahav (1999) Dekel, A., & Lahav, O. 1999, Astrophys. J., 520, 24
  * Dillon et al. (2014) Dillon, J. S., et al. 2014, Phys. Rev., D89, 023002
  * Doré et al. (2016) Doré, O., Werner, M. W., Ashby, M., et al. 2016, ArXiv e-prints, arXiv:1606.07039
  * Doré et al. (2018) Doré, O., Werner, M. W., Ashby, M. L. N., et al. 2018, ArXiv e-prints, arXiv:1805.05489
  * Fan et al. (2006) Fan, X., Strauss, M. A., Richards, G. T., et al. 2006, AJ, 131, 1203
  * Fry (1994) Fry, J. N. 1994, Phys. Rev. Lett., 73, 215
  * Furlanetto & Lidz (2007) Furlanetto, S., & Lidz, A. 2007, Astrophys. J., 660, 1030
  * Furlanetto et al. (2004) Furlanetto, S., Zaldarriaga, M., & Hernquist, L. 2004, Astrophys. J., 613, 1
  * Furlanetto et al. (2004) Furlanetto, S. R., Zaldarriaga, M., & Hernquist, L. 2004, ApJ, 613, 1
  * Gong et al. (2012) Gong, Y., Cooray, A., Silva, M., et al. 2012, ApJ, 745, 49
  * Gong et al. (2011) Gong, Y., Cooray, A., Silva, M. B., Santos, M. G., & Lubin, P. 2011, ApJ, 728, L46
  * Greig et al. (2013) Greig, B., Komatsu, E., & Wyithe, J. S. B. 2013, MNRAS, 431, 1777
  * Hoffmann et al. (2018) Hoffmann, K., Mao, Y., Mo, H., & Wandelt, B. D. 2018, ArXiv e-prints, arXiv:1802.02578
  * Howlett et al. (2015) Howlett, C., Manera, M., & Percival, W. J. 2015, Astronomy and Computing, 12, 109
  * Jensen et al. (2013) Jensen, H., et al. 2013, Mon. Not. Roy. Astron. Soc., 435, 460
  * Joachimi et al. (2009) Joachimi, B., Shi, X., & Schneider, P. 2009, A&A, 508, 1193
  * Kovetz et al. (2017) Kovetz, E. D., Viero, M. P., Lidz, A., et al. 2017, ArXiv e-prints, arXiv:1709.09066
  * Lagache (2017) Lagache, G. 2017, IAU Symp., 333, 228
  * Leclercq et al. (2013) Leclercq, F., Jasche, J., Gil-Marín, H., & Wandelt, B. 2013, J. Cosmology Astropart. Phys., 11, 048
  * Lidz et al. (2011) Lidz, A., Furlanetto, S. R., Oh, S. P., et al. 2011, ApJ, 741, 70
  * Lidz & Taylor (2016) Lidz, A., & Taylor, J. 2016, ApJ, 825, 143
  * Lidz et al. (2009) Lidz, A., Zahn, O., Furlanetto, S., et al. 2009, Astrophys. J., 690, 252
  * Lidz et al. (2008) Lidz, A., Zahn, O., McQuinn, M., Zaldarriaga, M., & Hernquist, L. 2008, ApJ, 680, 962
  * Loeb & Furlanetto (2013) Loeb, A., & Furlanetto, S. R. 2013, The First Galaxies in the Universe (Princeton University Press)
  * Majumdar et al. (2018) Majumdar, S., Pritchard, J. R., Mondal, R., et al. 2018, MNRAS, 476, 4007
  * Matarrese et al. (1997) Matarrese, S., Verde, L., & Heavens, A. F. 1997, Mon. Not. Roy. Astron. Soc., 290, 651
  * McQuinn et al. (2007) McQuinn, M., Lidz, A., Zahn, O., et al. 2007, MNRAS, 377, 1043
  * Mesinger & Furlanetto (2007) Mesinger, A., & Furlanetto, S. 2007, Astrophys. J., 669, 663
  * Mesinger et al. (2011) Mesinger, A., Furlanetto, S., & Cen, R. 2011, MNRAS, 411, 955
  * Planck Collaboration et al. (2016a) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016a, A&A, 594, A13
  * Planck Collaboration et al. (2016b) Planck Collaboration, Adam, R., Aghanim, N., et al. 2016b, A&A, 596, A108
  * Pritchard & Furlanetto (2007) Pritchard, J. R., & Furlanetto, S. R. 2007, Mon. Not. Roy. Astron. Soc., 376, 1680
  * Pullen et al. (2014) Pullen, A. R., Dore, O., & Bock, J. 2014, Astrophys. J., 786, 111
  * Righi et al. (2008) Righi, M., Hernández-Monteagudo, C., & Sunyaev, R. A. 2008, A&A, 489, 489
  * Schechter (1976) Schechter, P. 1976, ApJ, 203, 297
  * Scherrer & Weinberg (1998) Scherrer, R. J., & Weinberg, D. H. 1998, Astrophys. J., 504, 607
  * Scoccimarro (1997) Scoccimarro, R. 1997, ApJ, 487, 1
  * Scoccimarro et al. (2001) Scoccimarro, R., Feldman, H. A., Fry, J. N., & Frieman, J. A. 2001, Astrophys. J., 546, 652
  * Shimabukuro et al. (2016) Shimabukuro, H., Yoshiura, S., Takahashi, K., Yokoyama, S., & Ichiki, K. 2016, Mon. Not. Roy. Astron. Soc., 458, 3003
  * Shimabukuro et al. (2017) —. 2017, Mon. Not. Roy. Astron. Soc., 468, 1542
  * Silva et al. (2015) Silva, M., Santos, M. G., Cooray, A., & Gong, Y. 2015, ApJ, 806, 209
  * Smit et al. (2012) Smit, R., Bouwens, R. J., Franx, M., et al. 2012, ApJ, 756, 14
  * Smith et al. (2008) Smith, R. E., Sheth, R. K., & Scoccimarro, R. 2008, Phys. Rev., D78, 023523
  * Springel (2005) Springel, V. 2005, Mon. Not. Roy. Astron. Soc., 364, 1105
  * Tassev et al. (2013) Tassev, S., Zaldarriaga, M., & Eisenstein, D. J. 2013, J. Cosmology Astropart. Phys., 6, 036
  * Venemans et al. (2013) Venemans, B. P., Findlay, J. R., Sutherland, W. J., et al. 2013, ApJ, 779, 24
  * Verde et al. (2002) Verde, L., et al. 2002, Mon. Not. Roy. Astron. Soc., 335, 432
  * Visbal & Loeb (2010) Visbal, E., & Loeb, A. 2010, J. Cosmology Astropart. Phys., 11, 016
  * Zahn et al. (2006) Zahn, O., Lidz, A., McQuinn, M., et al. 2006, Astrophys. J., 654, 12
  * Zaldarriaga et al. (2004) Zaldarriaga, M., Furlanetto, S. R., & Hernquist, L. 2004, Astrophys. J., 608, 622
  * Zel’dovich (1970) Zel’dovich, Y. B. 1970, A&A, 5, 84

## Appendix A: 21cm Bias Results

Here we show how our $\langle T_{21}\rangle b_{21}$ results vary for different
choices of $k_{\text{min}}$ and $k_{\text{max}}$ (Fig. 10). As discussed in §
3.4, the results are noisy for small $k_{\text{min}}$, while the perturbative
description breaks down when $k_{\text{max}}$ is too large.

Figure 10.— Consistency between different estimates of $\langle T_{21}\rangle
b_{21}$. The black curves show the linear 21 cm bias factor, $\langle
T_{21}\rangle b_{21}$, inferred from the cross-power spectrum, $\langle
T_{21}\rangle b_{21}=P_{21,\delta}/P_{\delta,\delta}$ on large scales, as a
function of $\langle x_{i}\rangle$. These are compared with the linear bias
inferred from $Q_{21,21,21}$ (red) and $Q_{21,\delta,\delta}$ (blue). Left,
Center: Extracted $\langle T_{21}\rangle b_{21}$ for two different values of
$k_{\text{min}}$ and $k_{\text{max}}$. Right: The plot from Fig. 5 is
reproduced here for convenience. Note that the values of $\langle
T_{21}\rangle b_{21}$ inferred from the 21 cm auto-bispectrum at low $\langle
x_{i}\rangle$ show a marked difference from the cross-power spectrum and
cross-bispectrum inferences in the Central panel. We believe this relates to
the early-phase transition in the sign of the 21 cm-density cross-correlation
coefficient (see § 3.2)

## Appendix B: Cross-Bispectrum Variance

Here we write down a formula for the cross-bispectrum variance in the case
that some of the $k$-modes have equal magnitude. See e.g. Joachimi et al.
(2009) for more details. Recall that we have enforced $k_{1}\geq k_{2}\geq
k_{3}$.

$\mathrm{Var}[B_{21,\text{CII},\text{CII}}+\text{2
perm.}]=\begin{cases}P_{\text{tot,21}}(k_{1})P_{\text{tot,CII}}(k_{2})P_{\text{tot,CII}}(k_{3})+\text{2
perm.}&\text{if $k_{1}\neq k_{2}\neq k_{3}$}\\\ \\\\[12.0pt]
\begin{split}&2P_{\text{tot},21}(k_{3})P_{\text{tot},\text{CII}}(k_{1})^{2}\\\
&+2P_{\text{tot},21}(k_{1})P_{\text{tot},\text{CII}}(k_{3})P_{\text{tot},21}(k_{1})\\\
&+2P_{21,\text{CII}}(k_{1})^{2}P_{\text{tot},\text{CII}}(k_{3})\end{split}&\text{if
$k_{1}=k_{2}\neq k_{3}$}\\\ \\\\[12.0pt]
\begin{split}&2P_{\text{tot},21}(k_{1})P_{\text{tot},\text{CII}}(k_{2})^{2}\\\
&+2P_{\text{tot},21}(k_{2})P_{\text{tot},\text{CII}}(k_{1})P_{\text{tot},21}(k_{2})\\\
&+2P_{21,\text{CII}}(k_{2})^{2}P_{\text{tot},\text{CII}}(k_{1})\end{split}&\text{if
$k_{1}\neq k_{2}=k_{3}$}\\\ \\\\[12.0pt]
\begin{split}&6P_{\text{tot,21}}(k_{1})P_{\text{tot,CII}}(k_{1})P_{\text{tot,CII}}(k_{1})\\\
&+12P_{21,\text{CII}}(k_{1})^{2}P_{\text{tot,CII}}(k_{1})\end{split}&\text{if
$k_{1}=k_{2}=k_{3}$}\\\ \end{cases}$ (36)