# Reconstructing patchy helium reionization using the cosmic microwave
background and large-scale structure

Mesut Çalışkan<EMAIL_ADDRESS>William H. Miller III Department of Physics
and Astronomy, 3400 N. Charles St., Baltimore, Maryland 21218, USA  Neha Anil
Kumar<EMAIL_ADDRESS>William H. Miller III Department of Physics and
Astronomy, 3400 N. Charles St., Baltimore, Maryland 21218, USA  Selim C.
Hotinli<EMAIL_ADDRESS>William H. Miller III Department of
Physics and Astronomy, 3400 N. Charles St., Baltimore, Maryland 21218, USA
Perimeter Institute for Theoretical Physics, 31 Caroline St N, Waterloo, ON
N2L 2Y5, Canada  Marc Kamionkowski  William H. Miller III Department of
Physics and Astronomy, 3400 N. Charles St., Baltimore, Maryland 21218, USA

(August 28, 2024)

###### Abstract

The intergalactic helium became fully ionized by the end of cosmic noon
($z\sim 2$). Similarly to the reionization of hydrogen, helium reionization is
expected to be patchy, driven by luminous quasars that ionize the
intergalactic gas in their surrounding environment. Probing the morphology of
ionized electrons during this epoch can provide crucial information about
early structure formation, including the clustering and luminosities of
quasars, the accretion rates, variability, and lifetimes of active galactic
nuclei, as well as the growth and evolution of supermassive black holes. In
this study, we present how measurements of the cosmic microwave background
(CMB) can be used to reconstruct the optical-depth fluctuations resulting from
patchy helium reionization. As helium reionization occurred at lower
redshifts, upcoming probes of large-scale structure surveys will present a
significant opportunity to enhance the prospects of probing this epoch by
their combined analysis with the CMB. Using a joint information-matrix
analysis of hydrogen and helium reionization, we show that near-future galaxy
and CMB surveys will have enough statistical power to detect optical-depth
fluctuations due to doubly-ionized helium, providing a way of measuring the
redshift and duration of helium reionization to high significance. We also
show that modeling uncertainties in helium reionization can impact the
measurement precision of parameters characterizing hydrogen reionization.

###### Contents

  1. I Introduction
  2. II Helium reionization
    1. II.1 modeling the inhomogeneous optical depth
    2. II.2 Varying the model parameters
    3. II.3 Angular power spectrum of the optical depth
    4. II.4 Angular cross-power between the galaxy and ionized electron fluctuations
    5. II.5 Effects on the CMB
  3. III Patchy optical depth estimator
  4. IV Forecasts
    1. IV.1 Experiments
    2. IV.2 Detection SNR of the patchy optical depth
    3. IV.3 Probing helium reionization with the patchy optical depth and LSS
  5. V Discussion
  6. A Bubble model and the binned optical-depth spectra
  7. B Optical depth reconstruction noise
  8. C Improvements beyond ILC-cleaned foregrounds

##  I Introduction

Throughout the cosmic history of structure formation, the intergalactic medium
undergoes two significant transitions: the reionization epochs of hydrogen and
helium.111In this work we refer to the reionization of the second electron of
helium as “helium reionization.” While the epoch of hydrogen reionization has
attracted much attention from the cosmology and astrophysics communities [see,
e.g., 1, 2, for reviews], the focus on the reionization of helium has been
more limited [see, however, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32] despite its
observational accessibility. Since the second electron in helium has an
ionization energy of $\sim 54{\,\rm eV}$, photons expected to be produced by
stars do not doubly ionize the intergalactic helium efficiently at
cosmological scales. Diffuse helium remains singly ionized until populations
of active galactic nuclei and quasars start emitting sufficiently energetic
photons. Consequently, the epoch of helium reionization is expected to take
place at lower redshifts ($z\sim 3$), a time which is accessible through
observations of the ionizing sources. The morphology of this epoch is closely
linked to quasar number counts, clustering and luminosities [33, 34, 35, 36,
37], the formation of active galactic nuclei, their variability and lifetimes
[38, 39, 40], as well as early formation of supermassive black holes [41].
There is now also growing evidence that this epoch will be observationally
accessible from joint analyses of the CMB with other cosmological surveys
[e.g., 31, 32].

Here, we study the prospects of detecting and characterizing helium
reionization by extracting the inhomogeneous reionization signal from joint
analyses of upcoming cosmic microwave background (CMB) maps and galaxy
surveys. Similar to the epoch of hydrogen reionization, the reionization of
helium is anticipated to be anisotropic (or ‘patchy’) due to the formation and
growth of ionized bubbles around luminous sources. Unlike the reionization of
hydrogen, which is sourced by abundant hot stars, the reionization of helium
follows more dispersed sources such as quasars. This suggests that the
morphology of helium reionization (in addition to its redshift) should be
different than that of the hydrogen, holding valuable information about the
distribution and abundance of quasars (see, for example, Ref. [7]). The
patchiness from both epochs would lead to varying optical depth to
reionization across the sky, altering the statistics of the observed CMB maps
in multiple ways. These include the screening of the surface of the last
scattering, generation of new polarization via Thomson scattering from
reionization bubbles, and the kinetic Sunyaev Zel’dovich effect [42].

In this work, we use the statistical technique developed in Refs. [43, 44, 45,
46, 47], which was previously applied to extract the hydrogen reionization
signal from the CMB, together with a novel application of cross-correlating
the patchy optical depth with other surveys of large-scale structure (LSS),
such as galaxies. Given the increasing precision of the upcoming CMB
experiments, we fold in the epoch of helium reionization and calculate new
detection probabilities, forecasting the future ability to characterize each
of the independent epochs. We construct a quadratic estimator for the modes of
the patchy optical-depth field, which separates the patchy reionization signal
from CMB in the form of a noisy optical-depth map. In our forecasts, we cross-
correlate this estimator with tomographic measurements of the galaxy over-
density field while also considering its use as a standalone probe.

Unlike the optical depth reconstruction for hydrogen reionization [e.g., 44,
45], the potential science returns from probing patchy helium reionization
optical depth have not yet been extensively studied in the context of CMB. In
this study, we bridge this gap by focusing our forecasts on two upcoming CMB
experiments, namely, the stage-4 telescope CMB-S4 [48, 49, 50] and its stage-5
counterpart CMB-HD [51, 52, 53]. These upcoming experiments, alongside galaxy
surveys such as Vera Rubin Observatory (LSST) [54, 55] and MegaMapper [56],
stand to transform our understanding of cosmology with an influx of high-
precision data in the upcoming years. Among the new windows of opportunity
that are being opened by these experiments, the prospects of using the CMB as
a cosmological ‘back-light’ to probe LSS from signals sourced by interactions
between CMB photons and the intervening cosmological structures particularly
motivates this work.

These programs include the reconstruction of the lensing potential [see, e.g.,
57, for a review], radial and transverse velocity [58, 59, 60, 61, 62, 63, 64,
65, 66, 67], and quadrupole fields [66, 68, 69, 70, 71], as well as the patchy
optical depth, which we focus on in this paper. While the prospects of joint
analyses between reconstructed optical depth and tracers of LSS such as 21-cm
intensity mapping have been studied previously in literature [72, 73, 74] in
the context of hydrogen reionization, the use of these methods for probing
helium reionization has not been explored. Given that helium reionization
takes place at lower redshifts—a range well-covered by upcoming galaxy surveys
[54, 75], which will detect high numbers of galaxies—it stands to benefit
significantly from the statistical power provided by cross-correlations with
the galaxy distribution, in contrast to hydrogen reionization.

Little is known about the precise characteristics of helium reionization, such
as its morphology, timing, and duration. While surveys measuring helium and
hydrogen Lyman-$\alpha$ forests can, in principle, help characterize this
epoch, the analyses are subject to a range of systematic and astrophysical
uncertainties, such as those affecting the inferred flux levels of the
Lyman-$\alpha$ forests [76, 77, 78, 79], for example.

Specifically, measurements of hydrogen Lyman-$\alpha$ serve only as an
indirect probe of helium reionization, inferring insights through the thermal
history of the intergalactic medium in a model-dependent manner, and are
subject to systematic uncertainties. On the other hand, these measurements are
also significantly obscured by intervening Lyman-limit systems at lower
redshifts. Although further simulations and analyses will be needed to
understand whether these issues impose a significant limit on the prospects of
probing helium reionization from Lyman-$\alpha$ forests, it is nevertheless
suggestive that additional probes will be valuable for the unambiguous
characterization of this epoch.

In addition to enhancing our understanding of the astrophysics of quasars and
active galactic nuclei, probing helium reionization offers additional benefits
for cosmological inference. For example, the determination of the total change
in the free-electron fraction during this epoch serves as an indicator of the
primordial helium abundance $Y_{p}$. In turn, improved measurements of $Y_{p}$
can significantly enhance our understanding of weak interaction rates, neutron
lifetime, and Big Bang nucleosynthesis, as discussed in Refs. [80, 81]. Such
advancements provide valuable insights into the intricate details of our
cosmological history. As the effective number of degrees of freedom of light
relic species $N_{\rm eff}$ and $Y_{p}$ affect the small-scale CMB damping
tail in similar ways [82], improving the $Y_{p}$ measurement also improves
measurements of $N_{\rm eff}$, a driving science goal of the upcoming CMB
experiments.222This improvement is achieved, for instance, by breaking the
degeneracy between $Y_{p}$ and $N_{\rm eff}$ [83, 84, 85].

This paper is organized as follows. In Sec. II, we introduce our model of
patchy reionization characterizing the redshift evolution of the ionized
regions during both helium and hydrogen reionization. We introduce the primary
model of the integrated, auto-correlation signal $C_{\ell}^{\tau\tau}$,
accounting for optical-depth fluctuations sourced by both epochs. Furthermore,
we also present the central model for the cross-correlation of the patchy-
reionization optical-depth field with large-scale, tomographic galaxy survey
data. We conclude this modeling section by deriving the modifications to CMB
spectra arising from patchy reionization. We present the estimator for the
patchy optical-depth reconstruction in Sec. III. With the signal and noise
established, we finally show our results from the forecasts in Sec. IV. We
introduce the assumed experiment specifications and the resulting measurement
noise for near-future CMB experiments and galaxy surveys. We then forecast the
ability of these future surveys to probe optical-depth fluctuations sourced by
helium reionization via computation of the measurement SNR. We also forecast
the measurability of the morphology of this epoch by presenting results from
our information-matrix campaign. These results include the anticipated
fractional errors on the model parameters and an exploration of how our chosen
fiducial parameters impact the measurability of both helium and hydrogen
reionization. We conclude with a discussion in Sec. V. Appendix A contains
more details on the derivation of the binned auto-power spectrum of the
optical-depth field, as well as its cross-correlation with the galaxy density
field. Appendix B contains supplementary information on the derivation of the
noise in the CMB-reconstructed optical-depth field. Finally, Appendix C
contains selected forecast results for experiment configurations assuming a
more optimistic foreground cleaning methodology.

##  II Helium reionization

In this section, we detail the modeling of the inhomogeneous optical depth,
calculate the pertinent angular power spectra, and examine the observable
impacts of patchy reionization on the CMB. Section II.1 starts with a
definition of the optical depth of photons scattering off free electrons
originating from hydrogen and helium reionization. We then model the mean
ionization fraction and describe the assumed size distribution of ionized
regions. Section II.2 explores the influence of varying helium reionization
parameters. In Sec. II.3, we compute the angular power spectrum of the average
optical depth for both hydrogen and helium, incorporating the cross-power term
between the two. Section II.4 is dedicated to deriving the angular cross-power
between galaxy and ionized electron fluctuations and presents a comparative
analysis of all the binned power spectra. The section concludes with a
comprehensive discussion of the effects of patchy reionization on the CMB in
Section II.5.

###  II.1 modeling the inhomogeneous optical depth

The optical depth of photons scattering off of free electrons sourced by
reionization, out to redshift $z$, is

$\tau_{\theta}(\hat{\mathbf{n}},z)=\sigma_{T}n_{\theta,0}\int_{0}^{z}\frac{\differential
z^{\prime}(1+z^{\prime})^{2}}{H(z^{\prime})}x_{e}^{\theta}(\hat{\mathbf{n}},z^{\prime})\,,$
(1)

where we use the symbol $\theta\in\\{\rm H,He\\}$ to distinguish between the
contributions from reionization of hydrogen and helium, respectively. Here,
$\sigma_{T}$ is the Thomson scattering cross section, $H(z)$ is the Hubble
parameter at $z$, $x_{e}(\hat{\mathbf{n}},z)$ is the ionization fraction at
redshift $z$ in the line-of-sight direction $\hat{\mathbf{n}}$. We take
$n_{\theta,0}=\\{n_{p,0},f_{\rm He}n_{p,0}\\}$, where $n_{p,0}$ is the number
density of protons today, and $f_{\rm He}$ is the fraction of helium to
hydrogen atoms satisfying $f_{\rm He}\simeq 0.08$.

For each of the epochs, we model the mean ionization fraction with a
hyperbolic tangent

$\bar{x}^{\theta}_{e}(z)=\frac{1}{2}\left[1-\tanh\left(\frac{y(z)-y_{\rm
re}^{\theta}}{\Delta_{y}^{\theta}}\right)\right]\,,$ (2)

where $y(z)=(1+z)^{3/2}$, and $y_{\rm re}$ and $\Delta_{y}$ are free
parameters of the reionization model which roughly specify the time and
duration of the reionization process. Although simple, this functional form
has been used by various studies focusing on patchy reionization [e.g., 86,
87, 88] and enables us to explore the impact of various reionization times and
durations (i.e., when the epoch at which the reionization starts and how fast
the media transitions from nonionized to fully reionized) on the detection
signal-to-noise ratio (SNR) and parameter estimation, as we explain in later
sections.

For hydrogen reionization, we choose the fiducial values $\\{y_{\rm re}^{\rm
H},\Delta^{\rm H}_{y}\\}=\\{27.0,7.0\\}$, corresponding to a central redshift
of $z\sim 8$ and a duration spanning $11\gtrsim z\gtrsim 5$. This fiducial
choice is not only physically motivated by the approximate redshifts of star-
formation and measurements of the Lyman-$\alpha$ forests [1, 2], but is also
consistent with recent CMB measurements of the optical depth to recombination
setting $\bar{\tau}\approx 0.060$ [89, 90]. Similarly, for the helium
counterpart, we choose $\\{y_{\rm re}^{\rm He},\Delta^{\rm
He}_{y}\\}=\\{8.0,3.1\\}$, corresponding to a central redshift of $z\sim 3$
and a duration spanning $5\gtrsim z\gtrsim 1$. This choice is primarily
motivated by the expected redshifts of galaxy and quasar formation and its
impact on helium reionization [91]. Note that the relatively low abundance of
helium allows more freedom in the characterization of $\bar{x}_{e}^{\rm He}$,
with minimal impact on the average optical depth to reionization.

Figure 1: Upper panel: mean ionization fraction $\bar{x}_{e}(z)$ for hydrogen
(H) and helium (He) as a function of redshift $z$. Lower panel: spherical
bubble radius distribution $R\,P(R)$ as a function of bubble radius $R$. The
patchy reionization model is based on Eqs. (2) and (3), and the fiducial model
parameters describing the ionization fraction functions and the bubble radius
distribution are listed in Tab. 1 for each of the epochs.

Consistent with existing analysis and modeling of hydrogen reionization [e.g.,
92, 93, 94, 95, 45, 96], we represent ionized regions as spherical
‘bubble’-like volumes centered at the ionizing source. To allow for bubbles of
various sizes, we assume that the bubbles have radius $R$, distributed
according to a log-normal distribution as follows:

$P(R)=\frac{1}{R}\frac{1}{\sqrt{2\pi\sigma_{\ln
R}^{2}}}e^{-[\ln(R/\bar{R})]^{2}/(2\sigma_{\ln{R}}^{2})}\,,$ (3)

where $\bar{R}$ is the characteristic size of the bubbles, and
$\sigma_{\ln{R}}$ is the width of the distribution. In addition to this, we
assume that the number density of reionization bubbles is a biased tracer of
matter fluctuations on large scales, and thus we model the spatial
distribution of bubbles with a bubble bias $b^{\theta}$. Following Ref. [45],
we set $\\{\bar{R}^{\rm H},\sigma^{\rm H}_{\ln R}\\}=\\{5\
\rm{Mpc},\,\ln(2)\\}$ for hydrogen reionization. On the other hand, since
helium reionization is sourced by luminous quasars and active galactic nuclei
emitting hard photons, the typical size of bubbles during the reionization of
helium is anticipated to be larger than bubbles forming during the
reionization of hydrogen [e.g., 7]. Therefore, we choose $\\{\bar{R}^{\rm
He},\sigma^{\rm He}_{\ln R}\\}=\\{15\ \rm{Mpc},\,\ln(2)\\}$ for the helium
reionization.

The upper panel of Fig. 1 shows the mean ionization fraction as a function of
the redshift $z$ for hydrogen (helium) based on these fiducial parameters in
solid black (dashed red) line, while the lower panel represents the spherical
bubble radius distribution $RP(R)$ as a function of radius $R$. As seen,
helium reionization starts much later ($z\sim 5$) compared to that of hydrogen
and goes through a period of maximal patchiness around $z\sim 3$. As the lower
panel indicates, there is a lower probability for helium bubbles to be small
($R\lesssim 5\ \rm{Mpc}$). Although most bubbles are with $R\sim 15\
\rm{Mpc}$, the helium bubbles can be much larger than hydrogen bubbles (up to
$R\lesssim 100\ \rm{Mpc}$).

###  II.2 Varying the model parameters

Our choice of fiducial helium reionization parameters $\\{y_{\rm re}^{\rm
He},\Delta_{y}^{\rm He},\bar{R}^{\rm He},\sigma_{\ln R}^{\rm
He}\\}=\\{8.0,3.1,15.0\ {\rm{Mpc}},\ln{(2.0)}\\}$ is informed by the expected
time and duration of helium reionization, as well as the impact of the
difference in ionizing sources on the bubble size distribution [e.g., 7].
Nonetheless, there is great uncertainty as to what these parameters actually
are [76, 77, 78, 79]. The simplicity of our modeling enables us to explore a
wide range of possible reionization configurations by scanning the parameter
space of helium reionization parameters in the relevant ranges.

Figure 2: Left panels: average ionization fraction $\bar{x}_{e}(z)$ as a
function of redshift $z$ for various model parameters characterizing the
patchy helium (He) reionization. In the upper panel, $y_{\rm re}^{\rm He}=8.0$
while $\Delta_{y}^{\rm He}$ is varied. In the lower panel, the former is
varied while the latter is fixed at $\Delta_{y}^{\rm He}=3.1$. Right panels:
spherical bubble radius distribution $R\,P(R)$ as a function of bubble radius
$R$ for various model parameters. Similar to the left panels, $\sigma_{\ln
R}^{\rm He}$ ($\bar{R}^{\rm He}$) is fixed in the upper (lower) panel while
the other parameter is varied. In each case, the red line shows the results
for the set of fiducial parameters.

The left panels in Fig. 2 show the average ionization fraction $\bar{x}^{\rm
He}_{e}(z)$ for various model parameters. In the upper left panel, $y_{\rm
re}^{\rm He}$ is fixed at the fiducial value of $8.0$ while $\Delta_{y}^{\rm
He}$ is varied between $1.0$ and $6.0$ in blue and green, respectively. On the
other hand, in the bottom left panel, $\Delta_{y}^{\rm He}$ is fixed at the
fiducial value of $3.1$ while $y_{\rm re}^{\rm He}$ is varied between $3.0$
and $10.0$ in orange and black, respectively. Given the two plots on the left
of Fig. 2, we can understand the effects of changing $\Delta_{y}^{\rm He}$ and
$y_{\rm re}^{\rm He}$ on the redshift evolution of $\bar{x}^{\rm He}_{e}(z)$.
The upper left panel indicates that the parameter $\Delta_{y}^{\rm He}$
determines how smoothly the transition to reionization occurs. In other words,
the effect of lowering (increasing) $\Delta_{y}^{\rm He}$ is roughly speeding
up (slowing down) the rate at which helium ionizes as a function of redshift.
Note that changing $\Delta_{y}^{\rm He}$, with $y_{\rm re}^{\rm He}$ fixed,
does not affect the mean reionization redshift [redshift at which
$\bar{x}_{e}^{\rm He}(z)\approx 0.5$]. On the contrary, the bottom left panel
suggests that the parameter $y_{\rm re}^{\rm He}$ alters the midpoint’s
location, essentially shifting the mean redshift of reionization. Increasing
(decreasing) $y_{\rm re}^{\rm He}$ shifts the entire function to the right
(left) in redshift, thus changing the time at which reionization occurs. These
minimum and maximum ranges are determined to account for earlier/ later helium
reionization epochs as well as faster/ slower transitions to fully reionized
helium, and represent the limits of the ranges we consider in our forecasts.

The right panels in Fig. 2 show the spherical bubble radius distribution for
various parameter values. In the upper right panel, the width of the
distribution is fixed, while the characteristic bubble size $\bar{R}$ is
varied between $5.0\ \rm{Mpc}$ and $20.0\ \rm{Mpc}$ in blue and green,
respectively. In the lower right panel, the characteristic size is set to the
fiducial value while the width of the distribution is changed between
$\ln{1.5}$ and $\ln{4.0}$ in orange and black, respectively. These limits are
chosen to account for much larger (smaller) helium reionization bubbles than
expected, account for uncertainty in the width of the distribution, and
represent the limits of the ranges we consider in our forecasts. In each of
the subplots, the red line represents the results for the set of fiducial
parameters previously outlined for ease of comparison.

###  II.3 Angular power spectrum of the optical depth

Given that the observed optical depth $\tau$ will get contributions from both
the epochs of reionization, we use the derivations in Refs. [45, 93] to write
the angular power spectrum of the optical depth as a line-of-sight integral as
follows:

$C_{\ell}^{\tau\tau}=\underbrace{\int_{0}^{\chi_{*}}\differential\chi\frac{\sigma_{T}^{2}n_{p,0}^{2}}{a^{4}\chi^{2}}P^{\rm
H}_{\Delta_{e}\Delta_{e}}(\chi,k)}_{\mathrm{H-only}}+\underbrace{\int_{0}^{\chi_{*}}\differential\chi\,f_{\rm
He}^{2}\frac{\sigma_{T}^{2}n_{p,0}^{2}}{a^{4}\chi^{2}}P^{\rm
He}_{\Delta_{e}\Delta_{e}}(\chi,k)}_{\mathrm{He-
only}}+\underbrace{2\int_{0}^{\chi_{*}}\differential\chi\,f_{\rm
He}\frac{\sigma_{T}^{2}n_{p,0}^{2}}{a^{4}\chi^{2}}\bar{x}_{e}^{\rm
He}\bar{x}_{e}^{\rm H}P(\chi,k)}_{\mathrm{H-He\ cross}}\,,$ (4)

where $a$ is the cosmic scale factor, $\chi_{*}$ is the comoving distance to
recombination, $k$ is the Fourier wave number, $P(\chi,k)$ is the nonlinear
matter power spectrum, and we use the Limber approximation [97], setting
$k=\ell/\chi$. Here $P_{\Delta_{e}\Delta_{e}}^{\rm H}$ and
$P_{\Delta_{e}\Delta_{e}}^{\rm He}$ are the power spectra of hydrogen and
helium ionization fluctuations, respectively, which we define next. Note that
we arrive at this model by categorizing the free-electrons sourced by helium
reionization separately from those sourced by hydrogen, with the third term in
the above equation accounting for correlations in the distributions of the two
separate ‘populations.’ Therefore, the helium-fraction factor appropriately
weights each of the helium-dependent terms, in order to account for its
relatively low abundance.

Table 1: The set of reionization model parameters for hydrogen and helium and the galaxy bias factor. The middle columns present the set of fiducial values while the right-side columns show the range of values used in the forecasts described in Sec. IV. Parameter | Fiducial Value | Range of Values
---|---|---
|  Hydrogen | Helium |  Hydrogen | Helium
$y_{\rm re}^{\theta}$ | $27.0$ | $8.0$ | Fixed | $3.0\leq y_{\rm re}^{\rm He}\leq 10.0$
$\Delta_{y}^{\theta}$ | $7.0$ | $3.1$ | Fixed | $1.0\leq\Delta_{y}^{\rm He}\leq 6.0$
$\bar{R}^{\theta}$ | $5\ \rm{Mpc}$ | $15\ \rm{Mpc}\ \ $ | Fixed | $5.0\ \textrm{Mpc}\leq\bar{R}^{\rm He}\leq 20.0\ \textrm{Mpc}$
$\sigma_{\ln{R}}^{\theta}$ | $\ln 2.0$ | $\ln 2.0$ | Fixed | $\ln{1.5}\leq\sigma^{\rm He}_{\ln{R}}\leq\ln{4.0}$
Bubble bias $b^{\theta}$ | 6.0 | 6.0 | Fixed | Fixed
Galaxy bias factor $b_{g,0}$ | 0.95 | Fixed

To model the power-spectra $P_{\Delta_{e}\Delta_{e}}^{\theta}$ of ionization
fluctuations, we assume that the basic morphology of reionization remains the
same across the two epochs, and therefore, the functional form for the spectra
are identical:

$\begin{split}P_{\Delta_{e}\Delta_{e}}^{\theta}(\chi,k)&=P^{1b,\,\theta}_{\Delta_{e}\Delta_{e}}(\chi,k)+P^{2b,\,\theta}_{\Delta_{e}\Delta_{e}}(\chi,k)\,,\\\
\end{split}$ (5)

where, the 2-bubble $P^{2b,\,\theta}_{\Delta_{e}\Delta_{e}}(\chi,k)$ and
1-bubble $P^{1b,\,\theta}_{\Delta_{e}\Delta_{e}}(\chi,k)$ contributions are

$\begin{split}P^{2b,\,\theta}_{\Delta_{e}\Delta_{e}}(\chi,k)&=[(1-\bar{x}_{e}^{\theta})\ln(1-\bar{x}_{e}^{\theta}){{b^{\theta}}}\,I^{\theta}(k)-\bar{x}_{e}^{\theta}]^{2}P(\chi,k)\,,\\\
P^{1b,\,\theta}_{\Delta_{e}\Delta_{e}}(\chi,k)&=\bar{x}_{e}^{\theta}(1-\bar{x}_{e}^{\theta})[F^{\theta}(k)+G^{\theta}(k)]\,.\end{split}$
(6)

Here, the average ionization fraction $\bar{x}_{e}^{\theta}(z)$ varies from
zero to unity according to Eq. 2 (cf. Fig. 1), $b^{\theta}$ is the bubble bias
and is set to $6.0$ for both helium and hydrogen for simplicity, and terms
$\bar{x}_{e}^{\theta}$, $I^{\theta}(k)$, $G^{\theta}(k)$, $F^{\theta}(k)$ each
depend on the details of the reionization: the time and duration $\\{y_{\rm
re}^{\theta},\Delta_{y}^{\theta}\\}$ as well as the bubble size distribution
$\\{\bar{R}^{\theta},\sigma^{\theta}_{\ln R}\\}$. These spectra are modeled
following Ref. [45], details of which can be found in Appendix A.333The final
term in Eq. (4) comes from the free-electron density fluctuations which we
assume follow matter fluctuations, often referred to as the ‘late-time’
contribution to the electron fluctuations. The power-spectrum of the free-
electron fluctuations at the limit $z\rightarrow 0$ satisfies
$P_{ee}(\chi,k)=(1+f_{\rm He})^{2}n_{p,0}^{2}P(k)$. The set of fiducial
parameters used in modeling the angular power spectrum and the respective
range of values used in the forecasts of Sec. IV can be seen in Tab. 1. In
combination, there are 11 parameters (5 for the description of helium
reionization, 5 for hydrogen, and the galaxy bias factor $b_{g,0}$).

To understand the relative weight of each of the separate terms in Eq. (4),
Fig. 3 plots each contribution as a function of multipole $\ell$ for the
chosen set of fiducial parameters presented in Tab. 1. Relative suppression of
He-only and H-He cross terms is primarily sourced by the factor $f_{\rm He}$,
which manifests a significantly lower abundance of helium compared to
hydrogen. While the He-only term is slightly higher than the H-He cross term
at lower $\ell$, the latter dominates by about an order of magnitude at
smaller scales due to comparably greater hydrogen fraction.

Figure 3: Comparison of the angular power spectrum of optical-depth
fluctuations $C_{\ell}^{\tau\tau}$ as described in Eq. (4), where H-only
(solid black), He-only (dashed red), and the H-He cross terms (dot-dashed
blue) are given by the first, second, and third term of the equation,
respectively.  Figure 4: Angular power spectrum of the optical-depth
fluctuations $C_{\ell}^{\tau\tau}$ for various model parameters characterizing
patchy helium (He) reionization for the He-only (H-He cross) terms of Eq. (4)
in the left and middle (right) panels. In the left and right panels, the
bubble radius distribution parameters are fixed to the fiducial values while
the ionization fraction function parameters are varied. In the middle panel,
the ionization fraction function parameters are fixed to the fiducial values
while the bubble radius distribution parameters are varied. In each case, the
red line shows the results for the set of fiducial parameters.

To elucidate the influence of specific fiducial helium reionization parameters
on the helium-related terms in Eq. (4), Fig. 4 presents the relevant
components of $C_{\ell}^{\tau\tau}$ under different reionization scenarios.
The figure’s three panels highlight the maximum deviations in the $\tau$ auto-
correlation permissible within the parameter ranges specified in Tab. 1.

The left-most panel displays the minimum (orange) and maximum (black) possible
‘He-only’ contributions considered, holding the bubble parameters
$\bar{R}^{\rm He}$ and $\sigma_{\ln R}^{\rm He}$ constant while varying the
ionization fraction parameters. Similarly, the middle panel shows the impact
of the assumed He bubble size distribution by displaying the He-only
contributions arising from a $P(R)$ strongly preferring smaller bubbles
(orange) and one that more loosely prefers larger bubbles (black). This is
achieved by holding the ionization parameters $y_{\rm re}^{\rm He}$ and
$\Delta_{y}^{\rm He}$ constant while varying the bubble parameters to their
extremes. Finally, the third panel displays the possible contributions from
the ‘H-He Cross’ term arising from the widest variations in $y_{\rm re}^{\rm
He}$ and $\Delta_{y}^{\rm He}$ allowed by Tab. 1, similar to the left-most
plot. In all three cases, the solid line in red displays the expectation from
the assumed fiducial parameters, and the parameter variations considered
correspond to the differing ionization fraction functions and bubble
distributions plotted in Fig. 2. These extremal $C_{\ell}^{\tau\tau}$
encompass the variations considered in our forecasts in Sec. IV.

###  II.4 Angular cross-power between the galaxy and ionized electron
fluctuations

Looking at Eq. (4) and Fig. 3, one might naturally expect the signal from
helium reionization to be washed-out by the (relatively) larger hydrogen
counterpart in the $\tau$ auto-correlation measurement. However, noting that
the ionized helium bubbles are expected to form around AGN/quasars, which
trace the distribution of galaxies on large scales, we can amplify the helium
signal by cross-correlating optical depth measurements with galaxy survey
data. This cross-correlation takes the following functional form:

$C_{\ell}^{\tau
g}=\underbrace{\int_{0}^{\chi_{*}}\differential\chi\frac{\sigma_{T}n_{p,0}}{a^{2}\chi^{2}}P^{\rm
H}_{\Delta_{e}g}(\chi,k)}_{\mathrm{H-only}}+\underbrace{\int_{0}^{\chi_{*}}\differential\chi\frac{\sigma_{T}f_{\rm
He}n_{p,0}}{a^{2}\chi^{2}}P^{\rm He}_{\Delta_{e}g}(\chi,k)}_{\mathrm{He-
only}}\,.$ (7)

Here,

$P_{\Delta_{e}g}^{\theta}(\chi,k)=[\bar{x}^{\theta}_{e}-(1-\bar{x}^{\theta}_{e})\ln(1-\bar{x}^{\theta}_{e})b^{\theta}\,I^{\theta}(k)]b_{g}(z)P(\chi,k)\,,$
(8)

where $\bar{x}_{e}^{\theta}(z)$ similarly varying from zero to unity,
$b_{g}(z)=b_{g,0}(1+z)$ is the galaxy bias, and we use the Limber
approximation as in Eq. (4). We show the cross-correlation signal in the
absence of this approximation in Appendix A. Note that, unlike
$C_{\ell}^{\tau\tau}$, $C_{\ell}^{\tau g}$ consists only of a 2-bubble-like
term. The resulting effect is that $C_{\ell}^{\tau g}$ is a (relatively)
large-scale probe. While it can appreciably boost the signal, it may not be
able to carry sufficient information regarding the small-scale features of He
reionization, as we will discuss in Sec. IV.

Figure 5: Comparison of the angular cross-power spectrum between the galaxy
and the optical-depth field $C_{\ell}^{\tau g}$ based on our fiducial hydrogen
(H) and helium (He) reionization models, as described in Eq. 7 and 8. The
black, solid curve corresponds to the contribution from the ‘H-only’ term to
the total $C_{\ell}^{\tau g}$ signal [see Eq. (7)]. Similarly, the red, dashed
curve corresponds to the contribution from the ‘He-only’ term.

Figure 5 depicts the relative contribution from each of the terms in Eq. (7).
Although it appears as if, once again, the H term is far more dominant than
the He contribution, it is important to note that at the redshifts relevant to
galaxy observations ($z\lesssim 5$) the power in the ‘H-only’ term is coming
from the correlation between the matter and galaxy distribution on large
scales. That is, this signal varies smoothly at late-times and does not
contain information on the ‘patchy’ morphology of hydrogen reionization. On
the other hand, although the helium signal is subdominant in absolute terms,
it is actually this signal that causes significant redshift evolution of the
$\tau-g$ signal at late times. Therefore, the binned SNR and information-
matrix analysis performed in this work, that accounts for the tomographic
measurement of galaxy distributions, shows that most of the constraining power
in the measurement of parameters characterizing helium reionization is sourced
by this cross-correlation.

Finally, Fig. 6 displays the contributions of the He and H-He cross terms from
different redshift bins to the line-of-sight integrated optical depth signal
$C_{\ell}^{\tau\tau}$ and the cross-correlated signal $C_{\ell}^{\tau g}$.
These contributions are computed for 8 equal-width redshift bins between
$z=0.2$ and $z=5.0$. The contribution from each independent bin $[z_{\rm
low},z_{\rm high}]$ calculated by performing the line-of-sight integrations
[Eqs. (4) and (7)] from lower bound $\chi(z_{\rm low})$ to upper bound
$\chi(z_{\rm high})$.444A detailed description of the calculation of the
redshift binned power spectra, $C_{\ell,\alpha\alpha}^{\tau\tau}$ and
$C_{\ell,\alpha\alpha}^{\tau g}$, can be found in Appendix A.

The left-most plot of Fig. 6 shows the binned contributions to the ‘He-only’
term in the model for $C_{\ell}^{\tau\tau}$. As expected, the signal is
dominated by redshift bins spanning $5\gtrsim z\gtrsim 1.5$, which correspond
to the ‘patchy epoch’ of helium reionization as set by the fiducial parameters
(cf. Fig. 1). Note also that while the patchy epoch dominates the signal at
relatively large scales ($\ell\lesssim 800$), the angular power spectrum is
dominated by lower redshift bins at smaller scales ($\ell\gtrsim 800$). The
middle panel displays the contributions from the same set of bins to the ‘H-He
cross’ term. In contrast to the ‘He-only’ contributions, this plot indicates
that bins at later redshifts contribute more to the cross signal independent
of the scale. This behaviour can be attributed to the fact that our fiducial
model for reionization predicts that the patchy-hydrogen and -helium power
spectra approximate to the nonlinear matter power spectrum $P(\chi,k)$ at late
times. In the cross term, at the displayed redshifts, $\bar{x}_{e}^{\rm
H}\approx 1$ and therefore we integrate $P(\chi,k)$ modulated by
$\bar{x}_{e}^{\rm He}$ over increasingly larger redshift bins (in
$\chi$-space) at later times, explaining the increasing power with decreasing
redshift. Similar behaviour can be seen in the right-most plot as well, which
depicts the binned contributions to the ‘He-only’ term in the $C_{\ell}^{\tau
g}$ signal. This is once again explained by $\bar{x}_{e}^{\rm He}\rightarrow
1$, modulating the galaxy power spectrum $b_{g}P(\chi,k)$, integrated over
increasingly large redshift bins at later times. The three plots clearly
depict the power of the cross-correlation in absolute terms as well, with the
helium contribution to $C_{\ell}^{\tau g}$ being $\sim 4$ orders of magnitude
larger than both the helium-related terms in $C_{\ell}^{\tau\tau}$. Therefore,
the cross correlation of the optical-depth field with the large-scale galaxy
distribution will be an integral tool in the characterization of helium
reionization.

Figure 6: Comparison of the binned angular power spectrum of the optical depth
for He-only and H-He cross components (left and middle panels, respectively)
and the angular cross-power between the galaxy and the optical-depth
fluctuations for the He-only component (right panel). Different colors
correspond to the different redshift bins. The limits of the colors in the
color map are consistent with the limits of the redshift bins used.

Calculating the power spectra described above can be computationally
exhaustive, especially if a wide range of parameter values and experiment
configurations are considered. To improve the efficiency of the computations,
we utilize the Python package numba [98, 99], which is an open-source just-in-
time compiler that translates a subset of Python and NumPy into fast machine
code. This shortens the time required to evaluate the power spectra by up to
two orders of magnitude, allowing us to produce high-resolution binned spectra
of the modeled signal $C_{\ell}^{\tau\tau}$ and $C_{\ell}^{\tau g}$ for
varying fiducial parameter scenarios.

###  II.5 Effects on the CMB

Having modelled the integrated optical-depth signal in terms of reionization
parameters, we turn our attention to calculating the observable effects of
patchy reionization on the CMB temperature and polarization fluctuations. Most
generally, the patchy reionization leads to three effects on the CMB:

  * •

Screening of the CMB photons, where photon intensity gets multiplied by a
factor $e^{-\tau(\hat{\mathbf{n}})}$.

  * •

Generation of new polarization due to Thomson scattering of CMB photons off
electrons in ionized bubbles.

  * •

The kinetic Sunyaev Zel’dovich (kSZ) effect generated from the radial motion
of reionization bubbles relative to the observer.

In this work, we focus on the first two effects. The kSZ effect, which
requires careful modeling of small-scale electron fluctuations, also contains
valuable cosmological information [see, e.g., 100, 101], which we study in an
upcoming work.

The next step, therefore, is to propagate the effect of Thomson scattering
from reionization on to the CMB polarization and temperature maps. We start by
expressing the observed polarization signal along the line-of-sight as
follows:

$(Q\pm
iU)(\hat{\mathbf{n}})=\int_{0}^{\infty}\differential\chi\frac{\mathrm{d}\tau}{\mathrm{d}\chi}e^{-\tau(\chi\hat{\mathbf{n}})}S^{\pm}_{\rm
pol}(\chi\hat{\mathbf{n}})\,,$ (9)

where $Q$ and $U$ are the standard Stokes parameters, $\chi\equiv\chi(z)$ is
the comoving distance to redshift $z$, $S^{\pm}_{\rm
pol}(\chi\hat{\mathbf{n}})=-\sqrt{6}/10\sum_{m}\,{}_{\pm
2}Y_{2m}(\hat{\mathbf{n}})a^{T}_{2m}(\chi\hat{\mathbf{n}})$ is the local
temperature quadrupole that the electron sees,
$a_{2m}^{T}(\chi\hat{\mathbf{n}})$ are the temperature quadrupole moments, and
$\,{}_{\pm 2}Y_{2m}(\hat{\mathbf{n}})$ are the spin $\pm 2$ spherical
harmonics. Similar to Ref. [45], we parameterize the effect of reionization on
CMB polarization anisotropies as

$\begin{split}(Q\pm iU)(\hat{\mathbf{n}})=&(Q\pm
iU)_{0}(\hat{\mathbf{n}})+\sigma_{T}n_{p,0}\int\frac{\differential\chi}{a^{2}}x_{e}(\chi\hat{\mathbf{n}})(Q\pm
iU)_{1}(\chi\hat{\mathbf{n}})\,,\end{split}$ (10)

where $(Q\pm iU)_{0}(\hat{\mathbf{n}})$ is the polarization signal in the
absence of these effects, $x_{e}(\chi\hat{\mathbf{n}})$ is the spatial
fluctuations in the total ionized electron fraction introduced above, and

$(Q\pm iU)_{1}(\chi\hat{\mathbf{n}})=e^{\tau(\chi)}S^{\pm}_{\rm
pol}(\chi\hat{\mathbf{n}})-\int_{\chi}^{\infty}\differential\chi^{\prime}\frac{\differential\tau}{\differential\chi^{\prime}}e^{-\tau(\chi)}S_{\rm
pol}^{\pm}(\chi\hat{\mathbf{n}})\,,$ (11)

is the modulation to the polarization sourced by patchy reionization. Here,
$\tau(\chi)$ is the spatially-averaged mean optical depth at comoving distance
$\chi$. Throughout, we represent the CMB signals in redshift-binned form,
introducing $N$ redshift bins covering the patchy epoch of reionization. The
polarization signal takes the form

$(Q\pm iU)(\hat{\mathbf{n}})=(Q\pm
iU)_{0}(\hat{\mathbf{n}})+\sum\limits^{N}_{\alpha}\Delta\tau^{\alpha}(\hat{\mathbf{n}})(Q\pm
iU)_{1}(\chi_{\alpha}\hat{\mathbf{n}})\,,$ (12)

where $\Delta\tau^{\alpha}$ is the redshift-bin–averaged optical-depth
fluctuation. The effect of reionization on the CMB temperature anisotropy can
also be written in a form similar to polarization, satisfying

$T(\hat{\mathbf{n}})=T_{0}(\hat{\mathbf{n}})+\sum\limits^{N}_{\alpha}\Delta\tau^{\alpha}(\hat{\mathbf{n}})T_{1}(\chi_{\alpha}\hat{\mathbf{n}})\,,$
(13)

with

$T(\hat{\mathbf{n}})=\int_{0}^{\infty}\differential\chi
S_{T}(\chi\hat{\mathbf{n}})\,,$ (14)

where $S_{T}(\chi\hat{\mathbf{n}},\tau)$ is a function of local quantities,
dependent on the optical depth, albeit in a more complicated way due to
multiple contributing effects. The modulation
$T^{\alpha}_{1}(\chi\hat{\mathbf{n}})$ in this case can be expressed with a
functional integral of the form

$T_{1}^{\alpha}(\chi\hat{\mathbf{n}})=\int_{\chi_{\alpha}}^{\infty}\frac{\delta
T(\chi\hat{\mathbf{n}})}{\delta\tau(\chi)}\,,$ (15)

and can be found in Ref. [45], for example.

Since $Q$ and $U$ are not coordinate invariant, it is more convenient to re-
cast these maps in terms of the scalar $E$ and $B$ fields using spin-raising
and -lowering operators. The harmonic coefficients of the CMB $E$\- and
$B$-mode maps satisfy

$\begin{split}a^{E}_{\ell m}&\equiv\frac{1}{2}(\,_{+2}a_{\ell
m}+\,_{-2}a_{\ell m})\,,\\\ a^{B}_{\ell m}&\equiv\frac{1}{2i}(\,_{+2}a_{\ell
m}-\,_{-2}a_{\ell m})\,,\end{split}$ (16)

where

$\,{}_{\pm 2}a_{\ell m}=\int\differential^{2}\hat{\mathbf{n}}(Q\pm
iU)(\hat{\mathbf{n}})\,_{\pm 2}Y_{\ell m}(\hat{\mathbf{n}}).$ (17)

The harmonic coefficients for the $E$-mode polarization, $a_{\ell m}^{E_{0}}$
and $a_{\ell m}^{E_{1}^{\alpha}}$, then satisfy

$\begin{split}a_{\ell m}^{E_{0}}&=\frac{1}{2}(\,_{+2}a^{E_{0}}_{\ell
m}+\,_{-2}a^{E_{0}}_{\ell m})\,,\\\ a_{\ell
m}^{E_{1}^{\alpha}}&=\frac{1}{2}(\,_{+2}a^{E_{1}^{\alpha}}_{\ell
m}+\,_{-2}a^{E_{1}^{\alpha}}_{\ell m})\,,\end{split}$ (18)

where $a^{E_{0}}_{\ell m}=\int\differential^{2}\hat{\mathbf{n}}(Q\pm
iU)_{0}(\hat{\mathbf{n}})\,_{\pm 2}Y_{\ell m}(\hat{\mathbf{n}})$ for the
homogeneous contribution and $a^{E_{1}^{\alpha}}_{\ell
m}=\int\differential^{2}\hat{\mathbf{n}}(Q\pm
iU)_{1}(\chi_{\alpha}\hat{\mathbf{n}})\,_{\pm 2}Y_{\ell m}(\hat{\mathbf{n}})$
for the contribution from patchy reionization. Similarly, the temperature
harmonic coefficients from patchy reionization satisfy

$T_{1}(\chi_{\alpha}\hat{\mathbf{n}})=\sum\limits_{\ell m}a_{\ell
m}^{T_{1}(\chi_{\alpha})}Y_{\ell m}(\hat{\mathbf{n}})\,,$ (19)

and

$T_{0}(\hat{\mathbf{n}})=\sum\limits_{\ell m}a_{\ell m}^{T_{0}}Y_{\ell
m}(\hat{\mathbf{n}})\,.$ (20)

In order to calculate the effect of patchy reionization on the observed CMB
spectra, which we use to define optical-depth estimators in the next section,
it is sufficient to calculate the angular cross-correlation functions from
these coefficients, which are given by

$C_{\ell}^{E_{0}E_{1}^{\alpha}}\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}=\langle
a^{E_{0}}_{\ell m}a^{E_{1}^{\alpha}}_{\ell^{\prime}m^{\prime}}\rangle\,,$ (21)
$C_{\ell}^{T_{0}E_{1}^{\alpha}}\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}=\langle
a^{T_{0}}_{\ell m}a^{E_{1}^{\alpha}}_{\ell^{\prime}m^{\prime}}\rangle\,,$ (22)
$C_{\ell}^{T_{0}T_{1}^{\alpha}}\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}=\langle
a^{T_{0}}_{\ell m}a^{T_{1}^{\alpha}}_{\ell^{\prime}m^{\prime}}\rangle\,,$ (23)

where $\delta_{ij}$ is the Dirac delta function. The specific form of the
coefficients $E_{1}^{\alpha}$, $T_{1}^{\alpha}$ are given in Ref. [45].

##  III Patchy optical depth estimator

The three effects described in Sec. II.5 introduce a statistical anisotropy in
the small-scale CMB temperature and polarization, which can be used to
reconstruct the underlying fluctuations of the optical depth $\Delta\tau$. At
leading order in $\Delta\tau$, the cross-correlation of two CMB fields can be
written as

$\begin{split}\langle a^{X}_{\ell
m}a^{Y}_{\ell^{\prime}m^{\prime}}\rangle=&(-1)^{m}C_{\ell}^{X_{0}Y_{0}}\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}+\sum\limits^{N}_{\alpha}\\!\sum\limits_{LML^{\prime}M^{\prime}}\\!\\!\Delta\tau^{\alpha}_{L^{\prime}M^{\prime}}\langle
a^{X_{0}}_{LM}a^{Y_{1}^{\alpha}}_{L^{\prime}M^{\prime}}\rangle\int\differential^{2}\hat{\mathbf{n}}Y_{\ell
m}^{*}(\hat{\mathbf{n}})Y_{LM}(\hat{\mathbf{n}})Y_{L^{\prime}M^{\prime}}(\hat{\mathbf{n}})\\\
=&(-1)^{m}C_{\ell}^{X_{0}Y_{0}}\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}+\sum\limits^{N}_{\alpha}\sum\limits_{ML}\Delta\tau^{\alpha}_{LM}\Gamma_{\ell\ell^{\prime}L}^{X_{0}Y^{\alpha}_{1}}\begin{pmatrix}\ell&\ell^{\prime}&L\\\
m&m^{\prime}&M\end{pmatrix}\,,\end{split}$ (24)

where $\Delta\tau^{\alpha}$ is the redshift-bin averaged optical-depth
fluctuations, for which a biased minimum variance estimator can be written as

$\begin{split}&\widehat{\Delta\tau_{\ell
m}^{\alpha}}\,{\rm(biased)}=A_{\ell,\alpha\alpha}^{XY}(-1)^{m}\\!\\!\\!\\!\sum\limits_{LL^{\prime}MM^{\prime}}\\!\\!\begin{pmatrix}L&L^{\prime}&\ell\\\
\\!M&\\!M^{\prime}&m\end{pmatrix}\Gamma^{X_{0}Y^{\alpha}_{1}}_{LL^{\prime}\ell}\frac{a^{X}_{LM}a^{Y}_{L^{\prime}M^{\prime}}}{C_{L}^{XX,\rm
obs}C_{L^{\prime}}^{YY,\rm obs}}\,,\end{split}$ (25)

where

$A^{XY}_{\ell,\alpha\beta}=\left[\frac{1}{2\ell+1}\sum\limits_{\ell^{\prime}L}\frac{\Gamma^{X_{0}Y^{\alpha}}_{LL^{\prime}\ell}(\Gamma^{X_{0}Y^{\beta}_{1}}_{LL^{\prime}\ell})^{*}}{C_{L}^{XX,\rm
obs}C_{L^{\prime}}^{YY,\rm obs}}\right]^{-1}\,.$ (26)

The above estimator has a biased optical-depth reconstruction noise given by

${N}_{\ell,\alpha\beta}^{\tau\tau;XY}\,(\rm
biased)=\frac{A^{XY}_{\ell,\alpha\alpha}A_{\ell,\beta\beta}^{XY}}{A^{XY}_{\ell,\alpha\beta}}\,.$
(27)

We can also define an unbiased quadratic estimator for the optical depth as

$\widehat{\Delta\tau_{\ell
m}^{\alpha}}=(\bm{R}^{-1})_{\alpha\gamma}\widehat{\Delta\tau_{\ell
m}^{\gamma}}{{\,(\rm biased)}}\,,$ (28)

which satisfies $\langle\widehat{\Delta\tau_{\ell
m}^{\alpha}}\rangle=\Delta\tau_{\ell m}^{\alpha}$, where $\bm{R}$ is a
rotation matrix that de-biases the reconstructed optical depth, whose elements
can be found to satisfy

$R_{\alpha\beta}^{XY}=\frac{A^{XY}_{\ell,\alpha\alpha}}{A^{XY}_{\ell,\alpha\beta}}\,.$
(29)

The reconstruction noise for the unbiased estimator can then be calculated as

$N_{\ell,\alpha\beta}^{\tau\tau;XY}=(\bm{R}^{-1})_{\alpha\gamma}^{XY}(\bm{R}^{-1})_{\beta\delta}^{XY}\,\,\tilde{N}_{\ell,\gamma\delta}^{\tau\tau;XY}\,\,\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}\,,$
(30)

and we define a minimum-variance reconstruction noise by combining these
estimators as

${{N}_{\ell,\alpha\beta}^{\tau\tau;\rm
MV}}=(\sum_{XY}({\bm{N}_{\ell}^{\tau\tau;XY}})^{-1})^{-1}_{\alpha\beta}\,,$
(31)

where $X,Y\in\\{T,E,B\\}$ and $\alpha,\beta$ span the redshift-bin indices.
Note that we omit correlations between different redshift bins when
calculating the signal from Eq. (24), as we find these terms to be small and
inconsequential to our minimum-variance reconstruction noise. In close analogy
with CMB lensing reconstruction, the covariance of the reconstruction noise
spectra also includes off-diagonal cross-correlation terms between different
estimators. We surmise that taking into account these off-diagonal covariance
terms may increase the minimum-variance reconstruction noise up to around a
factor of $\mathcal{O}(0.1)$, similar to the case of lensing reconstruction.
We leave studying the impact of off-diagonal contributions to the lensing
reconstruction noise covariance to future work [see, e.g., 84, 102, for
calculations corresponding to lensing quadratic estimator]. We show our
calculations of different optical-depth estimators in Appendix B.

##  IV Forecasts

Building on the target optical depth signal outlined in Sec. II and the
measurement methodology and noise considerations for reconstructing this field
presented in Sec. III, we now advance to forecasting measurement accuracy for
forthcoming CMB and LSS surveys. These forecasts are based on the standard
Planck 2018 [90] $\Lambda\textrm{CDM}$ cosmology, employing the six parameters
specified in Table 2. Sec. IV.1 provides an overview of experimental setups,
including beam characteristics and noise assessments. In Sec. IV.2, we analyze
the auto-correlation SNR of the reconstructed patchy optical depth and the
optical-depth-galaxy cross-correlation SNR (Eqs. (33) and (34)), presenting
results for various experimental scenarios and exploring the impact of He
reionization parameter variations. Finally, Sec. IV.3 delves into the
measurability of H and He reionization parameters, assessing potential
parameter degeneracies and the influence of varying He reionization
parameters.

###  IV.1 Experiments

Table 2: Fiducial cosmological parameters for the 6-parameter $\Lambda$CDM model considered in our calculations throughout this paper matching Planck 2018 [90]. Parameter | Fiducial Value
---|---
Cold dark matter density ($\Omega_{c}h^{2}$) | 0.120
Baryon density ($\Omega_{b}h^{2}$) | 0.022
Angle subtended by acoustic scale ($\theta_{s}$) | 0.010410
Optical depth to recombination $(\tau)$ | 0.060
Primordial scalar fluctuation amplitude ($A_{s}$) | $2.196\\!\times\\!10^{-9}$
Primordial scalar fluctuation slope ($n_{s}$) | 0.965
Table 3: Inputs to ILC noise for the baseline CMB configurations: The beam and temperature noise RMS parameters are chosen to roughly match CMB-S4 [48] and CMB-HD [52]. We model the CMB noise as shown in Eq. (32). In both cases, we account for the degradation due to Earth’s atmosphere by including the CMB red-noise with $\ell_{\rm knee}=100$ and $\alpha_{\rm knee}=-3$. The polarization noise satisfies $\Delta_{E}=\Delta_{B}=\sqrt{2}\Delta_{T}$. | Beam FWHM | Noise RMS
---|---|---
|  | ($\mu$K-arcmin)
| S4 | HD | S4 | HD
39 GHz | $5.1^{\prime}$ | $36.3^{\prime\prime}$ | 12.4 | 3.4
93 GHz | $2.2^{\prime}$ | $15.3^{\prime\prime}$ | 2.0 | 0.6
145 GHz | $1.4^{\prime}$ | $10.0^{\prime\prime}$ | 2.0 | 0.6
225 GHz | $1.0^{\prime}$ | $6.6^{\prime\prime}$ | 6.9 | 1.9
280 GHz | $0.9^{\prime}$ | $5.4^{\prime\prime}$ | 16.7 | 4.6

We model atmospheric and instrumental noise contributions to the CMB
temperature and polarization as

$N_{\ell}^{\rm TT}=\Delta_{T}^{2}\exp\left(\frac{\ell(\ell+1)\theta^{2}_{\rm
FWHM}}{8\ln 2}\right)\left[1+(\ell_{\rm knee}/\ell)^{\alpha_{\rm
knee}}\right]\,,\\\ $ (32)

assuming $N_{\ell}^{EE}=N_{\ell}^{BB}=2N_{\ell}^{\rm TT}$, where $\Delta_{T}$
is the detector RMS noise, and $\theta_{\rm FWHM}$ is the Gaussian beam full
width at half maximum following internal-linear-combination (ILC) cleaning
[103]. The term inside the square-brackets in Eq. (32) corresponds to the
‘red’ noise due to Earth’s atmosphere, parameterised by the parameters
$\ell_{\rm knee}$ and $\alpha_{\rm knee}$. We define our choices for these
parameters in Table 3, which we set to match the ongoing and upcoming CMB
surveys—CMB-S4 [104, 105, 106] and CMB-HD [51, 52].

To account for residual contamination from foregrounds to the CMB, we consider
Poisson and clustered cosmic infrared background (CIB), as well as the thermal
Sunyaev-Zel’dovich (tSZ) foreground, which we calculate following Refs. [107,
108]. We omit the cross-correlation between tSZ and CIB. We include radio
point sources following Ref. [109]. In addition to frequency-dependent
contributions, we also consider black-body, late-time, and reionization kSZ
and calculate the lensed CMB black-body using CAMB [110]. In what follows, we
perform forecasts for both ILC-cleaned CMB and ‘black-body’ CMB, where for the
latter, we omit contributions from frequency-dependent foregrounds but include
residual white noise after ILC cleaning.

For the galaxy density, we consider a galaxy catalogue with specifications
anticipated to match the photometric LSST survey [111]. We approximate the
galaxy density of the “gold” sample with
$n_{\text{gal}}(z)=n_{0}[({z}/{z_{0}}]^{2}\exp(-z/z_{0})/{2z_{0}}$, where
$n_{0}=40\leavevmode\nobreak\ \text{arcmin}^{-2}$, $z_{0}=0.3$, and we take
the galaxy bias as $b_{g}(z)=0.95(1+z)$. While accounting for the standard
anticipated photometric redshift (photo-$z$) error $\sigma_{z}=0.03(1+z)$, it
is noteworthy that the galaxy redshift bin sizes chosen in this analysis—eight
equal-redshift bins within the range $z\in[0.2,5.0]$—are considerably larger
than the photometric redshift errors at all redshifts. As a result, we
anticipate that photo-$z$ errors will not significantly impact our results.
While this study does not include it, a future spectroscopic survey such as
MegaMapper could, in principle, enable much finer redshift binning. This would
enhance the statistical significance of the joint analysis with the CMB data.
We leave forecasting for MegaMapper to future work.

###  IV.2 Detection SNR of the patchy optical depth

To evaluate the statistical power of upcoming CMB and galaxy surveys in
detecting optical-depth fluctuations caused by the ionized second electron of
helium, we initially investigate the SNR of the reconstructed patchy optical
depth. We define the auto-correlation SNR as

${\rm SNR}^{2}\leavevmode\nobreak\ ({\rm
auto})=\sum_{\ell\ell^{\prime}}\sum_{\alpha\beta\gamma\delta}C_{\ell}^{{{\tau}}_{\alpha}{{\tau}}_{\beta}}\textbf{cov}^{-1}\left(\tilde{C}^{{\tau}_{\alpha}{\tau}_{\beta}}_{\ell},\tilde{C}^{{\tau}_{\gamma}{\tau}_{\delta}}_{\ell^{\prime}}\right)C_{\ell^{\prime}}^{{{\tau}}_{\gamma}{{\tau}}_{\delta}}\,,$
(33)

where the index $\alpha$ labels a redshift bin, such that
$C_{\ell}^{{\tau}_{\alpha}{\tau}_{\beta}}$ represents the correlation of
$\tau$ across two redshift bins centered at $z_{\alpha}$ and $z_{\beta}$. Note
that spectra with tilde represent the observed spectra
$\tilde{C}_{\ell}^{{\tau}_{\alpha}{\tau}_{\beta}}\equiv
C_{\ell}^{{\tau}_{\alpha}{\tau}_{\beta}}+N^{\tau\tau;\rm
MV}_{\ell,\alpha\beta}$, including the optical-depth reconstruction noise
$N^{\tau\tau;\rm MV}_{\ell,\alpha\beta}$ defined in Eq. (31), and the optical-
depth signal $C_{\ell}^{\tau_{\alpha}\tau_{\beta}}$, which we will assume to
be diagonal in redshift bins.

Using the same notational conventions, we define the reconstructed optical-
depth–galaxy cross-correlation SNR as

${\rm SNR}^{2}\leavevmode\nobreak\ (\rm
cross)=\sum_{\ell\ell^{\prime}}\sum_{\alpha\beta\gamma\delta}C_{\ell}^{{\tau}_{\alpha}g_{\beta}}\textbf{cov}^{-1}\left(\tilde{C}^{{\tau}_{\alpha}g_{\beta}}_{\ell},\tilde{C}^{{\tau}_{\gamma}g_{\delta}}_{\ell^{\prime}}\right)C_{\ell^{\prime}}^{{\tau}_{\gamma}g_{\delta}}\,.$
(34)

In both the SNR equations, the covariance satisfies

$\textbf{cov}\left(\tilde{C}^{XY}_{\ell},\tilde{C}^{WZ}_{\ell^{\prime}}\right)=\frac{\delta_{\ell\ell^{\prime}}}{2\ell+1}f_{\mathrm{sky}}^{-1}\left(\tilde{C}_{\ell}^{XW}\tilde{C}_{\ell}^{YZ}+\tilde{C}_{\ell}^{XZ}\tilde{C}_{\ell}^{YW}\right)\,,$
(35)

where $f_{\rm sky}$ is the survey coverage sky fraction.

We show the SNR forecasts for a range of experimental specifications in Figs.
7 and 8. Figure 7 shows the detection SNR as a function of $L_{\rm max}$ for
experiment specifications corresponding to CMB-S4 (green colored lines
labelled S4) and CMB-HD (blue colored lines labelled S5). The dashed, dot-
dashed, and solid lines correspond to SNR obtained from the the auto-
correlation of the reconstructed optical-depth field $\langle\tau\tau\rangle$
[cf. Eq. (33)], cross correlation of reconstructed optical-depth field with
LSST $\langle\tau g\rangle$ [cf. Eq. (34)], and the total of the two,
respectively. The left panel displays the resulting SNR from considering only
the black-body contributions to the CMB (in addition to white noise). In
contrast, the right panel corresponds to using ILC-cleaned CMB spectra
including all foregrounds as described in Sec. IV.1. Throughout, we take
$f_{\rm sky}=1$ due to the uncertainty in the sky coverage of upcoming
experiments. The detection SNR can be straight-forwardly scaled by
$\sqrt{f_{\rm sky}}$ based on the coverage of a given survey. For CMB-S4, the
anticipated sky coverage fraction is $\sim 0.4$. For the joint sky area
surveyed by LSST and CMB-S4, $f_{\rm sky}$ is also expected to be $\sim 0.4$,
while CMB-HD is conceived to cover over half of the sky [52].

Figure 7: Detection SNR as a function of maximum multiple $L_{\rm max}$
considered, depicting the statistical power of CMB and LSS experiments for
detecting optical-depth fluctuations from the doubly-ionized helium. The
survey sky fraction $f_{\rm sky}$ is normalized to unity. The left panel
corresponds to considering only the black-body contribution to the CMB. The
right panel corresponds to using ILC-cleaned CMB spectra, including residual
frequency-dependent foregrounds as described in Sec. IV.1. The dotted lines
labeled $\langle\tau\tau\rangle$ correspond to the detection SNR from auto-
correlation of the reconstructed $\tau$ [including both the He-only and H-He
‘cross’ terms in Eq. (4)]. The dot-dashed lines labeled $\langle\tau g\rangle$
correspond to detection SNR from cross-correlating the reconstructed optical
depth and observed galaxy number density with survey specifications matching
LSST [He-only term in Eq. (7)]. Solid lines correspond to the combination of
both signals. The green (labelled ‘S4’) and blue (labelled ‘S5’) correspond to
the assumption of a CMB-S4-like and CMB-HD-like CMB experiment, respectively.
We find that near-future CMB and LSS experiments will have enough statistical
power to detect optical-depth fluctuations from reionized helium, with
detection significance reaching $\gtrsim 3\sigma$ ($\gtrsim 10\sigma$) from
optical-depth reconstruction up to scales $L_{\rm max}\sim 3000$ for CMB-S4
(CMB-HD) together with LSST, or $L_{\rm max}\sim 4000$ if $f_{\rm sky}=0.4$.
For most $L_{\rm max}$ values, the SNR is dominated by the cross-correlation
of the reconstructed optical depth and galaxy field, whereas the reconstructed
optical depth can be detected as a standalone probe with a CMB-HD-like survey
after ILC cleaning.  Figure 8: The dependence of the detection SNR of optical-
depth fluctuations from helium reionization on the CMB white-noise RMS
$\Delta_{T}$ in ${\rm\mu K}$-arcmin. Here, $L_{\rm max}$ is fixed to $6000$.
The left (right) panel corresponds to black-body (ILC-cleaned) CMB spectra.
The blue (green) line corresponds to a $1.4^{\prime}$ ($15^{\prime\prime}$)
CMB beam representing a CMB-S4-like (CMB-HD-like) telescope. Here, $f_{\rm
sky}$ is normalized to unity. Detection SNR can be appropriately scaled by
$\sqrt{f_{\rm sky}}$ based on the sky-coverage fraction of a given survey.
Note that while improving the beam quality alone does not significantly
enhance the SNR at a given noise level, it is important to consider that the
white noise RMS $\Delta_{T}$ for CMB-S4 (CMB-HD) is approximately $1\mu{\rm
K}$-arcmin ($0.2\mu{\rm K}$-arcmin).

We find that a joint analysis of (near future) CMB-S4 and LSST will have
enough statistical power to make a high-significance ($\gtrsim 3\sigma$)
detection of optical depth fluctuations sourced by the second electron in
helium, once the field is reconstructed up to scales $L_{\rm max}\sim 3000$
(or $L_{\rm max}\sim 4000$, if $f_{\rm sky}=0.4$), following standard ILC
cleaning. An S5 experiment like CMB-HD, moreover, could reach beyond $\gtrsim
10\sigma$ detection significance for the same $L_{\rm max}$ choices. If the
reconstructed optical depth could be utilized up to smaller scales, $L_{\rm
max}\simeq 10000$ for example, detection SNR for both survey combinations can
be boosted by more than an order of magnitude.555The feasibility of this
improvement depends on the modeling uncertainties of the nonlinear
contributions to galaxy and reconstructed optical depth.

The statistical power of these measurements depends significantly on the
residual foregrounds after standard ILC cleaning, as can be seen by comparing
the ILC-cleaned results (right panel of Fig. 7) to the detection SNR forecasts
assuming black-body CMB (left panel of Fig. 7), which are more optimistic. We
also find, for example, if frequency-dependent-foregrounds could be removed
more substantially, the optical-depth fluctuations can be detected from auto-
correlations with a CMB-HD-like survey at a comparable SNR to that of the
$\tau$-galaxy cross-correlation. Recent improvements of the ILC-cleaning
methods, such as the needlet ILC techniques discussed in Ref. [112, and
references therein], suggest there is room for improvement. Here, we show both
scenarios to demonstrate the lower and upper limits of the statistical power
of our analysis, given our survey specifications.

Figure 8 further demonstrates the dependence of the detection SNR on the CMB
RMS noise $\Delta_{T}$ in $\rm\mu K$-arcmin for a fixed choice of $L_{\rm
max}$. Here, our S4 (S5) survey corresponds to a $1.4^{\prime}$
($15^{\prime\prime}$) beam.

Figure 9:  Helium reionization optical-depth detection SNR for a range of
values of helium-reionization parameters $y_{\rm re}^{\rm He}$ and
$\Delta_{y}^{\rm He}$. Here, $y_{\rm re}^{\rm He}$ sets the mean-redshift of
reionization while $\Delta_{y}^{\rm He}$ determines the duration (cf. Fig. 2).
The bubble distribution parameters are fixed to their fiducial values
indicated in Tab. 1, and the sky-coverage faction $f_{\rm sky}$ is normalized
to unity. The top, middle, and bottom rows set the maximum multipole of the
reconstructed optical depth to 4095, 6095 and 9995, respectively. The left
(right) column corresponds to a CMB-HD-like (CMB-S4-like) survey specification
after ILC cleaning. We find that longer duration (slower transition) combined
with later reionization (patchy epoch coinciding with galaxies) leads to
higher SNRs. In principle, the SNR here can be used to put bounds on helium-
reionization parameters. In case the optical-depth signal from helium is not
detected above $\textrm{SNR}\sim 5$ with CMB-S4 for $L_{\rm max}=4095$ (or
$\sim 3\sigma$ if $f_{\rm sky}=0.4$), for example, this could suggest helium
reionization might have happened more rapidly, i.e. smaller $\Delta_{y}^{\rm
He}$, or started at a an earlier redshift, i.e. lower $y_{\rm re}^{\rm He}$.
Figure 10:  The maximum CMB RMS noise level, $\Delta_{T}^{(\rm max)}$, that
warrants an $\textrm{SNR (cross)}\geq 3$ detection for the optical-depth
signal, for a range of ionization fraction function parameter values. Here,
$L_{\rm max}$ is assumed to be $6000$, and the bubble size distribution
parameters are fixed to their fiducial values indicated in Tab. 1. The right
(left) panel corresponds to a CMB-S4-like (CMB-HD-like) beam. Unlike previous
detection-SNR-related figures, here we set the sky fraction $f_{\rm sky}=0.4$,
matching anticipated sky coverage from combinations of CMB-S4 and CMB-HD with
LSST. Note the significant variation of $\Delta_{T}^{(\rm max)}$ values with
varying $y_{\rm re}^{\rm He}$ and $\Delta_{y}^{\rm He}$ values. This figure is
similar to Fig. 9 in demonstrating that an unambiguous detection of optical
depth at a certain significance can in principle be used to pose bounds on the
helium-reionization parameter space.  Figure 11:  Similar to Fig. 9, except
bubble size distribution parameters $\bar{R}^{\rm He}$ and $\sigma^{\rm
He}_{\ln{R}}$ are varied while the ionization fraction parameters are set to
their fiducial values highlighted in Tab. 1. Unlike $y_{\rm re}^{\rm He}$ and
$\Delta_{y}^{\rm He}$ (cf. Fig. 9), a variation in the bubble size
distribution parameters has little effect on the detection SNR. This is
expected, as the detection SNR is dominated by the $\langle\tau
g\rangle$-cross signal, which is only weakly sensitive to the bubble size
distributions considered in our study.

While the results shown in Figs. 7 and 8 are based on the fiducial hydrogen
and helium reionization parameters outlined in Tab. 1, we show the dependence
of the optical-depth detection SNR on varying helium reionization parameters
in Figs. 9 and 10. Figure 9 shows the detection SNR as a function varying
helium reionization parameters $y_{\rm re}^{\rm He}$ and $\Delta_{y}^{\rm
He}$. Here, the bubble size distribution parameters are fixed at their
fiducial values, and $f_{\rm sky}$ is normalized to unity. Similar to before,
the detection SNR can be scaled by $\sqrt{f_{\rm sky}}$ depending the sky
coverage fraction of a given survey. The top, middle and bottom panels
correspond to $L_{\rm max}$ choices of 4095, 6095, and 9995, respectively. The
left (right) column corresponds to CMB-HD (CMB-S4) like survey specifications.
In all cases, we find the detection SNR depends significantly on the helium
reionization parameter values, indicative of the sensitivity of the total
signal to characteristics of helium reionization. In addition to the trivial
increase of SNR with increasing $L_{\rm max}$, the inclination of the fixed-
SNR contours indicates that detection prospects improve with decreasing
$y_{\rm re}^{\rm He}$ and increasing $\Delta_{y}^{\rm He}$. In other words,
the detection SNR for helium reionization is higher when helium reionization
occurs at a relatively later redshift over a longer duration, with highest SNR
achieved for $\bar{x}_{e}^{\rm He}$ modelled according to the orange line in
Fig. 2. An increased duration corresponds to an extended period of
‘patchiness’, leading to larger effect of the optical depth fluctuations on
the integrated signal measurement. Furthermore, $\langle\tau g\rangle$ is the
dominant contribution to the detection SNR (see Figs. 7 and 8), with
tomographic galaxy measurements improving in later redshift bins. As a result,
a later epoch of reionization allows for a stronger cross-correlation between
the optical-depth and the galaxy-density field leading to a higher detection
SNR.

Similarly, Fig. 10 shows the highest RMS CMB noise level
$\Delta_{T}^{\rm(max)}$ resulting in an $\textrm{SNR}\geq 3$ detection for the
reconstructed optical depth for a $1.4^{\prime}$ ($15^{\prime\prime}$) beam on
the left (right) panel. Note that, unlike the previous detection SNR figures,
we set $f_{\rm sky}=0.4$ here. Again, we find that varying helium reionization
parameters leads to significant variations of $\Delta_{T}^{\rm(max)}$,
suggesting a detection of the optical-depth fluctuations at a given SNR and
experimental noise can in principle be translated to upper or lower limits on
the helium reionization parameter space.

Figure 11 shows the same as Fig. 9 except bubble size distribution parameters
$\bar{R}^{\rm He}$ and $\sigma^{\rm He}_{\ln{R}}$ are varied while the
ionization fraction parameters are set to their fiducial values highlighted in
Tab. 1. We find that variation in the bubble size distribution has little
effect on the detection SNR compared to the previous case (cf. Fig. 9). This
is expected, as the detection SNR is dominated by the $\tau\ \textrm{--}\ g$
cross signal, which is mainly sensitive to scales larger than that of bubble
sizes considered in this study. We omit the version of Fig. 10 where the
bubble size distribution parameters are varied since, similar to before,
variation of these parameters does not greatly change the results.

The detection SNR results shown in this section demonstrate the total
statistical power of CMB and galaxy surveys for detecting optical-depth
fluctuations from the ionized second electron in helium, including
contributions from electron fluctuations during late times after the end of
helium reionization. As a result, they do not directly correspond to the
detection prospects of helium reionization parameters, which we investigate in
the next section with information-matrix forecasts.

###  IV.3 Probing helium reionization with the patchy optical depth and LSS

To assess future prospects of characterizing helium reionization via the
detection of reionization model parameters from CMB and LSS (galaxies), we
define an ensemble-information matrix as

$\mathcal{F}_{ik}=\sum\limits_{XYWZ}\sum\limits_{\ell=L_{\rm min}}^{L_{\rm
max}}\\!\\!\frac{\partial
C_{\ell}^{XY}}{\partial{\pi_{i}}}\textbf{cov}^{-1}\left(\tilde{C}^{XW}_{\ell},\tilde{C}^{YZ}_{\ell^{\prime}}\right)\frac{\partial
C_{\ell^{\prime}}^{W{Z}}}{\partial{\pi_{k}}}\,,$ (36)

where
$XY,ZW\in\\{{\tau}_{\alpha}{\tau}_{\beta},{\tau}_{\alpha}g_{\beta},g_{\alpha}g_{\beta}\\}$
for every unique pair of redshift bins centered at $z_{\alpha}$ and
$z_{\beta}$. Here, $\partial C_{\ell}^{X}/\partial\pi_{i}$ represents the
derivative of the signal with respect to the parameter $\pi_{i}$. Our
parameter array includes five helium reionization parameters $\\{y_{\rm
re}^{\rm He},\Delta_{\rm re}^{\rm He},b_{\rm He},\sigma_{\rm\ln R}^{\rm
He},\bar{R}^{\rm He}\\}$, as well as another five parameters characterising
hydrogen reionization $\\{y_{\rm re}^{\rm H},\Delta_{\rm re}^{\rm H},b_{\rm
H},\sigma_{\rm\ln R}^{\rm H},\bar{R}^{\rm H}\\}$ and the galaxy bias $b_{g}$,
as defined in Sec. II.

Figure 12: Fractional errors of helium (right panel) and hydrogen (left panel)
reionization parameters. The $y$-axes are fractional errors defined as
$\Delta(\pi_{i})/\pi_{i}$, where $\pi=\\{\bar{R}^{\theta},\sigma_{\ln
R}^{\theta},\Delta_{y}^{\theta},y_{\rm re}^{\theta},b^{\theta}\\}$ with
$\theta=\\{\rm He,H\\}$ for helium and hydrogen reionization parameters,
respectively. The $x$-axes correspond to the highest multipole of the
reconstructed optical depth $L_{\rm max}$ taken in the information-matrix
calculation. Here, we consider a CMB-S4 like survey after ILC cleaning,
together with LSST. The gray-shaded region corresponds to fractional errors
larger than unity, below which parameters could be detected. We find that
while the optical-depth reconstruction will provide high-accuracy measurements
of all hydrogen-reionization-model parameters we consider, measuring helium
reionization may be more difficult. For our fiducial helium-reionization model
choices, only the redshift of the patchy epoch of helium reionization may be
measured to high precision with CMB-S4 and LSST. We note, however, that these
results depend on our modeling choices and if helium reionization occurs at
lower redshift or lasts longer, detection prospects could be more optimistic
as shown in Fig. 14.

Our data array consists of reconstructed optical-depth fluctuations in 16
redshift bins within the range $z\in[0.2,20]$—8 bins distributed within
$z\in[0.2,5.0]$ and another 8 bins within $z\in[5.0,20]$, with equal spacing
in redshift. We include the galaxy field only for the 8 lower-redshift bins
with otherwise identical redshift spacing. It is important to note that
although our information matrix construction appropriately assumes that the
signals at different $\ell$-modes are independent, it does not assume that the
observed field is uncorrelated across different redshift bins. As stated
earlier, we do assume that the optical depth field and the galaxy over-density
field show minimal correlation across redshift bins, the non-zero correlation
is purely sourced by the unbiased ${\tau}$-reconstruction noise
${{N}_{\ell,\alpha\beta}^{\tau\tau}}$ derived in Sec. III.

With the previously defined experiment configurations, and the above
assumptions on redshift-binning and signal correlations, we finally construct
an 11-dimensional information matrix via Eq. (36) to estimate our constraining
power on the fiducial reionization parameters (see Tab. 1). Note that the
inclusion of both H and He reionization signals in this single information
matrix allows for forecasts that account for possible parameter degeneracies
between the two epochs. Furthermore, as we will show below, an exploration of
the helium reionization parameter space under this information-matrix
formalism will allow for predictions on the impact of the helium signal on our
ability to constrain standard hydrogen reionization parameters. Note that we
only perform the following analysis with the inclusion of the galaxy cross-
correlation i.e., we do not provide parameter constraints using the optical
depth field alone because the unbiased reconstruction noise
${{N}_{\ell,\alpha\beta}^{\tau\tau}}$ is too large to constrain all 10
reionization parameters using the $C_{\ell}^{\tau\tau}$ signal alone.

The fractional errors on the fiducial helium and hydrogen reionization
parameters can be found in Figs. 12 and 13. Figure 12 shows the fractional
errors from combination of CMB-S4 and LSST, whereas Fig. 13 displays the same
results, this time assuming an experimental configuration corresponding to
CMB-HD and LSST. Both sets of results assume that the optical-depth signal is
obtained from the CMB using ILC cleaning and the constraints are displayed as
a function of the maximum reconstructed multipole $L_{\rm max}$ in Eq. (36).
In each case, the right (left) panel corresponds to fractional errors on
hydrogen (helium) reionization parameters.

Figure 13:  Similar to Fig. 12, but in this case, we consider a CMB-HD-like
experiment with ILC cleaning. We find that, for sufficiently high values of
$L_{\rm max}$, both of the helium reionization parameters characterising the
redshift and duration of the patchy epoch can be measured to high accuracy. We
also note that the measurement accuracy for hydrogen reionization parameters
increases significantly compared to CMB-S4, suggesting future CMB surveys such
as CMB-HD, when combined with LSST, will play a significant role in
characterizing both hydrogen and helium reionization to high precision.
Figure 14: Dependence of fractional errors on varying values of helium
reionization parameters $y_{\rm re}^{\rm He}$ and $\Delta_{y}^{\rm He}$,
similar to Fig. 9, for a CMB-S4 like survey together with LSST. Here, we set
$L_{\rm max}=6000$ and assume the CMB spectra to be ILC cleaned. The lower two
panels in the right column correspond to fractional errors on helium
ionization-fraction parameters. The remainder of the panels correspond to
parameters modeling hydrogen reionization. We find that for all parameters,
fractional errors depend on the fiducial choice of $y_{\rm re}^{\rm He}$ and
$\Delta_{y}^{\rm He}$. This is particularly important for probing hydrogen
reionization, for which the constraints on parameters such as $\Delta_{y}^{\rm
H}$ and $y_{\rm re}^{\rm H}$ vary by up to a factor of $\mathcal{O}(5)$ within
the range of the parameter space we consider here, with errors becoming larger
with shorter duration of helium reionization (lower $\Delta_{y}^{\rm He}$) and
with higher values of $y_{\rm re}^{\rm He}$. This is due to the degeneracy
between helium and hydrogen reionization parameters in case helium
reionization occurs before hydrogen reionization ends completely. Our results
suggest that reionization of helium could, in principle, play a limiting role
on the non-ambiguous characterisation of hydrogen reionization from optical-
depth reconstruction.  Figure 15: Similar to Fig. 14, except the results are
for a CMB-HD-like survey together with LSST. As seen, the measurability of
both the He ionization-fraction parameters improve compared to case with
CMB-S4, and both $y_{\rm re}^{\rm He}$ and $\Delta_{y}^{\rm He}$ can be
constrained across most of the grid. This is followed with an additional
improvement to the measurement of H reionization parameters compared to the
previous case. These results indicate that the upcoming CMB and galaxy surveys
will have the statistical power to probe the midpoint of helium reionization
in redshift and the duration of the patchy epoch (parametrized by $y_{\rm
re}^{\rm He}$ and $\Delta_{y}^{\rm He}$, respectively) to a sufficient
precision via optical-depth reconstruction cross-correlated with the
statistically powerful large-scale structure data sets.  Figure 16:
Dependence of the measurability of He bubble size distribution parameters,
$\bar{R}^{\rm He}$ and $\sigma^{\rm He}_{\ln{R}}$, and the Helium bubble bias
$b^{\rm He}$ on varying $\bar{R}^{\rm He}$ and $\sigma^{\rm He}_{\ln{R}}$ for
CMB-HD-like survey specifications after ILC cleaning together with LSST. For
simplicity, we only show the results for $L_{\rm max}=6000$. Colored regions
represent where the fractional errors $\Delta\pi_{i}/\pi_{i}<1$, while the
regions where $\Delta\pi_{i}/\pi_{i}\geq 1$ are colored in black. Results
indicate $\bar{R}^{\rm He}$, $\sigma^{\rm He}_{\ln{R}}$, and $b^{\rm He}$ can
be constrained if the bubble sizes are small enough ($\bar{R}^{\rm He}\lesssim
6\textrm{--}9$ Mpc) and the distribution width is low $\sigma^{\rm
He}_{\ln{R}}\lesssim\ln{1.8}$. Note that the range of $\bar{R}^{\rm He}$ and
$\sigma_{\ln{R}}^{\rm He}$ displayed here is narrower than that in previous
plots to enhance the visibility of the measurable parameter space.

In agreement with forecasts in the literature [see, e.g., 113, 96, 114, 115,
116, 47, 117, 118, 119], we find that optical depth reconstruction, via both
the assumed CMB experiment specifications, will allow for high-precision
measurements of all hydrogen reionization parameters, further demonstrating
the potential value of this program in the near future. In contrast, the
information-matrix results indicate that characterizing the helium epoch will
be more cumbersome. For the combination of CMB-S4 and LSST, we find that
fractional errors reach values below unity only for $y_{\rm re}^{\rm He}$,
indicating that only the redshift of the patchy epoch (midpoint) of helium
reionization may be constrained.

Nonetheless, for the combination of CMB-HD and LSST, we see an improvement,
with below-unity fractional errors achieved for both $y_{\rm re}^{\rm He}$ and
$\Delta_{\rm re}^{\rm He}$ at sufficiently high $L_{\rm max}$. This is
followed with an additional improvement to the measurement of hydrogen
reionization parameters compared to the previous case. This implies that
upcoming CMB and galaxy surveys will have the statistical power to constrain
the midpoint of helium reionization in redshift and the duration of the patchy
epoch (parametrized by parametrized by $y_{\rm re}^{\rm He}$ and $\Delta_{\rm
re}^{\rm He}$, respectively) to a sufficient precision via optical-depth
reconstruction cross-correlated with the statistically powerful galaxy survey
data sets.

Advances in ILC-cleaning methods can, in principle, further improve these
constraints, as shown in Figs. 17 and 18, resulting in, for example, multiple
helium reionization parameters being detected for sufficiently large $L_{\rm
max}$ values even with CMB-S4. In all the experiment configurations
considered, under the fiducial parameter values (Tab. 1), the helium bubble-
size-distribution parameters $\\{\bar{R}^{\rm He},\sigma_{\ln R}^{\rm He}\\}$
and bubble bias $b^{\rm He}$ remain unmeasurable
($\Delta\pi_{i}/\pi_{i}>1.0$). This is because the constraining power during
the epoch of He reionization is primarily sourced by the cross-correlation
$\langle\tau g\rangle$. Our current model for this correlation only accounts
for galaxies as a large-scale tracer of bubbles i.e.; it only accounts for
correlations of the 2-bubble term of galaxy distributions and optical depth
fluctuations, where the bubble bias is degenerate with fluctuations in the
bubble size distribution [see Eq. 8]. Therefore, relying solely on
$P_{\Delta_{e}\Delta_{e}}^{1b,{\rm He}}$ for constraints on the bubble
parameters does not result in enough power to characterize the morphology of
ionized regions during this epoch.

Similarly to our SNR analysis, we once again vary the assumed set of fiducial
He reionization parameters to asses the impact of the morphology of this epoch
on the measurement prospects of not only the helium parameters but also their
hydrogen counterparts. We first consider the resulting parameter constraints
under variations of the redshift dependence/evolution of the helium
reionization epoch. We run the information-matrix formalism, varying the
assumed shape of $\bar{x}_{e}^{\rm He}(z)$ with $3.0\leq y_{\rm re}^{\rm
He}\leq 10.0$ and $1.0\leq\Delta_{y}^{\rm He}\leq 6.0$ (as summarized in Tab.
1), holding all other parameters fixed at their fiducial values.666Note that
the model-extremes of $\bar{x}_{e}^{\rm He}$ [$P(R)^{\rm He}$] considered in
this analysis correspond to the two variations plotted in the top [bottom]
panel of Fig. 2, over the assumed fiducial model in red. The results of these
parameter variations, assuming ILC-cleaned CMB observations, are summarized in
Figs. 14 and 15, with the former (latter) displaying the relevant results for
the experiment configuration corresponding to CMB-S4 (CMB-HD) and LSST. In
each of the figures, the top row presents the effects of varying $y_{\rm
re}^{\rm He}$ and $\Delta_{y}^{\rm He}$ on the fractional error of the
hydrogen bubble parameters $\bar{R}^{\rm H}$ and $\sigma_{\ln R}^{\rm H}$.
Similarly, the middle (bottom) panel displays the effects of this variation on
the duration (mean-reionization redshift) of each of the epochs, characterized
by $\Delta_{y}^{\rm H}$ and $\Delta_{y}^{\rm He}$ ($y_{\rm re}^{\rm H}$ and
$y_{\rm re}^{\rm He}$). We do not display the fractional error contour plots
for helium bubble parameters $\bar{R}^{\rm He}$, $\sigma_{\ln R}^{\rm He}$,
and the bubble bias $b^{\rm He}$, since these parameters are not measurable
($\Delta\pi_{i}/\pi_{i}>1.0$) at any point in the grid.

Across both the considered experiment configurations, the results indicate
that if helium reionization occurs at lower redshifts (lower $y_{\rm re}^{\rm
He}$ values), or lasts a longer time (larger $\Delta_{y}^{\rm He}$ values)
fractional errors on $y_{\rm re}^{\rm He}$ improve significantly, reducing by
around an order of magnitude within the span of the parameter space we
consider. The behaviour of fractional errors in $\Delta_{y}^{\rm He}$ is more
complicated, with information-matrix results indicating improved constraints
for mid-range $y_{\rm re}\sim 5.0$ with shorter duration of transition (lower
$\Delta_{y}^{\rm He}$).

Furthermore, it is clear that the measurability of some hydrogen reionization
parameters is also visibly impacted, up to a factor of $\mathcal{O}(5)$, by
variations in $y_{\rm re}^{\rm He}$ and $\Delta_{y}^{\rm He}$. The parameters
characterizing $\bar{x}_{e}^{\rm H}$ see decreasing fractional errors with
increasing duration of the helium epoch and decreasing mean redshift, whereas
the hydrogen bubble parameters see minimal change across the grid. This
behaviour can be attributed to an enhanced degeneracy between the helium and
hydrogen reionization parameters if helium reionization occurs at higher
redshifts, where the cross-correlation $\langle\tau g\rangle\rightarrow 0$ ,
before hydrogen reionization comes to a complete end. Therefore, Figs. 14 and
15 suggest that characterization of hydrogen reionization from the
reconstructed optical depth can, in principle, be biased by helium
reionization.

Our analysis indicates that varying the He bubble size distribution parameters
$\bar{R}^{\rm He}$ and $\sigma^{\rm He}_{\ln{R}}$ have minimal $(\lesssim
10\%)$ impact on the measurability of hydrogen reionization parameters.
Similarly, the measurability of the helium ionization fraction parameters
$y_{\rm re}^{\rm He}$ and $\Delta_{y}^{\rm He}$ also do not depend heavily on
varying He bubble parameters. Nonetheless, we find that varying the He bubble
size distribution have significant impact on the measurability of He bubble
parameters $\bar{R}^{\rm He}$ and $\sigma^{\rm He}_{\ln{R}}$ as well as the He
bubble bias $b^{\rm He}$.

To quantify the extent of this impact, we vary helium bubble distribution
$P(R)^{\rm He}$, with $5.0\ \textrm{Mpc}\leq\bar{R}^{\rm He}\leq 20.0\
\textrm{Mpc}$ and $\ln 1.5\leq\sigma_{\ln R}^{\rm He}\leq\ln 4.0$ (as
summarized in Tab. 1), holding all other parameters fixed at their fiducial
values6. Figure 16 shows the fractional errors for $\bar{R}^{\rm He}$ (left),
$\sigma_{\ln R}^{\rm H}$ (middle) and helium bubble bias $b^{\rm He}$ (right).
Here, we show the results for $L_{\rm max}=6000$ and assume that the CMB
spectra are ILC cleaned.

Despite the concentrated effects of bubble-parameter variations, all three
contour plots agree on an important take-away—although initial information
matrix results on the physically motivated, fiducial parameter space may
indicate that some helium parameters may never be measurable, results can
change starkly with changes in the assumed bubble distribution. Specifically,
the three previously unmeasurable ($\Delta\pi_{i}/\pi_{i}>1.0$) parameters
$\bar{R}^{\rm He}$, $\sigma_{\ln R}^{\rm He}$, and $b^{\rm He}$ can be
measured ($\Delta\pi_{i}/\pi_{i}<1.0$) if the regions of ionized helium are
smaller ($\bar{R}^{\rm He}\lesssim 6$ Mpc). The conditions are slightly looser
for $\sigma_{\ln R}^{\rm He}$ in specific, which remains measurable up to
$\bar{R}^{\rm He}\lesssim 9$ Mpc, as long as the distribution of bubbles is
sharply peaked at the mean value. This behaviour can be attributed to the fact
that the ionized bubble distribution strongly effects the small-scale 1-bubble
term of $P^{\rm He}_{\Delta_{e}\Delta_{e}}$ appearing in our model for
$C_{\ell}^{\tau\tau}$. Smaller bubble size boosts the power in the $P^{1b,{\rm
He}}_{\Delta_{e}\Delta_{e}}$ making the model more sensitive to changes in the
$P(R)^{\rm He}$. The effect of smaller characteristic bubble sizes can be
seen, for example, in the middle panel (orange line) of Fig. 4. This manifests
in our information-matrix analysis in the form of increasingly tight
constraints on $\bar{R}^{\rm He}$, $\sigma_{\ln R}^{\rm He}$, and $b^{\rm
He}$.

##  V Discussion

Over the course of cosmic history, the intergalactic medium experienced two
pivotal transitions: the reionization of hydrogen and helium. While hydrogen
reionization, primarily driven by photons from hot stars, has been extensively
studied, helium reionization, occurring at lower redshifts around $z\sim 3$,
has garnered less attention. Helium reionization is anticipated to be
anisotropic, like that of hydrogen, characterized by the formation and growth
of ionized bubbles around luminous sources [3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31,
32]. However, the reionization of helium is driven by energetic photons from
active galactic nuclei and quasars, since the second electron in helium
requires $\sim 54{\rm eV}$ of energy to be ionized. As a result, the
morphology of helium reionization can provide valuable insight into quasar
abundance, clustering, luminosities, variability, and lifetimes, and the early
formation of supermassive black holes [33, 34, 35, 36, 37, 38, 39, 40].

In this study, we have evaluated the potential of forthcoming surveys to probe
signatures of helium reionization through the reconstruction of the optical-
depth field $\tau$. We repurpose existing statistical tools of optical depth
reconstruction, previously developed for probing hydrogen reionization in
Refs. [43, 44, 45, 46, 47], and combine it with the statistical power of
galaxy surveys, taking advantage of the late-time occurrence of helium
reionization to forecast that the characterization of the morphology of this
epoch is achievable in the near future.

We characterize the epoch of helium reionization by modeling the evolution of
the mean ionization fraction of helium $\bar{x}_{e}^{\rm He}(z)$ and the size
distribution of ‘bubble’-like regions of ionized helium. The modeling choices
are an extension of the functional forms previously used to characterize the
epoch of hydrogen reionization (as seen in Refs. [45, 46]), with an updated
parameter space to account for the later occurrence of the helium epoch and
the difference in the ionizing sources. Based on this morphology, we model the
expected optical-depth-fluctuation power spectrum $C_{\ell}^{\tau\tau}$,
critically accounting for both helium and hydrogen reionization (as an
extension of the model presented in Ref. [45]). For the first time, we also
present a model for the cross-correlation of the $\tau$-field with the galaxy
density field that accounts for both of these separate epochs. Additionally,
we present a derivation for the expected noise in the reconstruction of the
optical-depth field, assuming that optical-depth measurements are derived from
their effects on the CMB temperature and polarization maps.

With the noise and signal defined, we present forecasts on our ability to
measure the effects of optical-depth fluctuations sourced by helium
reionization in the form of measurement SNR. To estimate how well the
conjunction of future CMB experiments and galaxy surveys may be able to
characterize the reionization morphology—including the time, duration, and
patchiness of the epoch—we perform a comprehensive information-matrix analysis
that incorporates the parameter space characterizing both the helium and
hydrogen reionization. Thus, we also account for any possible modeling
degeneracies that may arise due to the similarities in the signal from the two
separate epochs.

In our forecasts, we primarily consider two sets of experimental
configurations, both of which approximate the specifications of either CMB-S4
[48] or CMB-HD [52]. In each case, we assume that the CMB-reconstructed
optical depth is cross-correlated with the tomographic galaxy survey data
obtained from the LSST [55]. We perform forecasts for varying choices of
$L_{\rm max}$ (maximum multipole of $\tau$-reconstruction) and different
assumptions on the CMB foreground. The conservative forecasts assume that the
optical-depth field is obtained from the CMB using ILC cleaning, while the
more optimistic ones assume that the CMB dataset is reduced only to its black-
body contributions. Note also that we omitted the forecasting for MegaMapper
[56], which is planned to be the spectroscopic follow-up to LSST. Such a
survey will measure galaxy redshifts to high precision, allowing the use of
the full three-dimensional galaxy distribution much more thoroughly. We leave
the inclusion of this survey in forecasting to future work.

Through this analysis, we find that the characterization of the epoch of
helium reionization is primarily driven by the power in the cross-correlation
$\langle\tau g\rangle$, irrespective of the assumed experiment configuration.
In specific, while combination of CMB-S4 and LSST will have sufficient
statistical power to detect optical-depth fluctuations sourced by the epoch of
helium reionization, connecting these detections to an unambiguous measurement
of helium reionization parameters may be more difficult. Nevertheless, our
information-matrix forecasts (on fiducial parameters) indicate that this
experiment configuration will likely be able to constrain the mid-point/mean-
redshift of reionization, with constraints on the duration being achievable
only for more optimistic assumptions on the cleaning of the CMB maps and
higher $L_{\rm max}$ of $\tau$-reconstruction.

On the other hand, forecasts on the experiment configuration corresponding to
CMB-HD and LSST indicate significant improvements, with the detection SNR for
$\tau$-fluctuations sourced by helium increasing above $\mathcal{O}(10)$ for
reasonable choices of foreground cleaning and $L_{\rm max}$. This improvement
is also manifested in the parameter measurements forecasted via the
information-matrix formalism. Our analysis suggests that the combination of
CMB-HD and LSST will be able to probe the redshift evolution of helium
reionization to high fidelity, pinning down the time and duration of the
epoch. In both the considered experiment configurations, unfortunately, the
patchy bubble parameters and bubble bias remain elusive under the fiducial
assumptions on bubble size distributions.

Considering the limited constraints on the epoch of helium reionization and
the subsequent flexibility in modeling, we strive to accommodate the effects
of our chosen fiducial parameters. This is achieved by conducting (SNR and
information-matrix) forecasts over a broad parameter space, thereby
encompassing a diverse range of morphologies characterizing this epoch. We
vary the redshift evolution (both time and duration) of reionization and
bubble size distribution to explore the parameter space where higher
measurement SNRs and improved parameter constraints are likely. Our analysis
indicates that, for both experimental configurations, the measurement SNR of
$\tau$-fluctuations sourced by helium reionization improves significantly when
the reionization happens later, lasts longer, and the ionized regions are
smaller. In fact, for both experiment configurations, fractional errors on the
mean redshift of reionization see an improvement by up to an order of
magnitude as the epoch of reionization is pushed to lower redshifts and longer
durations. Furthermore, we find that the probing the bubble size distribution
becomes accessible for CMB-HD and LSST when the average size of ionized
regions is small (radius of $\lesssim 7$ Mpc).

Moreover, since the modeling of the signal and the resulting information-
matrix formalism accounts for both the hydrogen- and helium-reionization
parameter space, we explore the effects of different helium reionization
morphologies on our ability to characterize its hydrogen counterpart. Across
both the experiment configurations, the fractional errors in the hydrogen
reionization parameter space are affected the most by the changes in the
redshift evolution of the helium reionization epoch. In fact, fractional
errors from our information-matrix analyses indicate that the characterization
of the hydrogen epoch is easiest (and likely least biased) when the helium
reionization epoch occurs at later times for longer durations. This behaviour
can be attributed to the alleviation of degeneracies across the two similarly
modelled epochs, in models where helium reionization begins after hydrogen
reionization is over. Furthermore, the galaxy dataset can play a stronger role
in the alleviation of this degeneracy if the patchiness of helium reionization
happens at lower redshifts (where the galaxies are more abundant) and lasts
for a longer time (longer duration of $\langle\tau g\rangle$ correlation
sourced by helium reionization). As a result, this work suggests that the
inclusion of helium reionization in future measurements of the hydrogen
reionization epoch, especially when using line-of-sight integrated signals,
will be integral to ensure unbiased results on the hydrogen parameter space.

It is vital for us to address the fact that the results presented in this
paper are dependent on our simplified model of helium reionization. This
choice was made not only to mirror existing forecasts on probing hydrogen
reionization with similar techniques [e.g., 113, 96, 114, 115, 116, 47, 117,
118, 119], but also to model the epoch of helium reionization without
including too many astrophysical priors. In practice, both the bubble-size
distribution and its mean will have redshift dependence set by the abundance
of quasars (or active galactic nuclei) and their efficiency in ionizing their
environment [e.g., 15]. Given a distribution of sources, the mean ionization
fraction can, in principle, be predicted from the typical expansion speed of
ionizing bubbles. This speed may vary based on the intergalactic environment,
suggesting that it is not strictly a free parameter. Such improvements in
helium (and hydrogen) reionization model building can lead to different SNR
values as suggested by our analysis varying the fiducial model parameters.
Furthermore, given the similarity of the modeling of the two epochs, the
forecasts presented here are heavily impacted by parameter degeneracies across
the two, i.e., forecasts might improve if we model the epoch of helium
reionization using more priors from astrophysics, specifically given that this
epoch is visibly accessible. We leave improvements to the modeling of helium
reionization to our upcoming works on the subject. Our choices here provide a
general assessment on the statistical power accessible with CMB and galaxy
surveys and will guide these future efforts.

###### Acknowledgements.

We are grateful to Simone Ferraro, Gil Holder, Amalia Madden, and Kendrick
Smith for useful discussions and collaboration. M.Ç., N.A.K, and M.K.
acknowledge support from the John Templeton Foundation through Grant No.
62840. M.Ç. is also supported by NSF Grants No. AST-2006538, PHY-2207502,
PHY-090003 and PHY-20043, NASA Grants No. 20-LPS20-0011 and 21-ATP21-0010, and
by Johns Hopkins University through the Rowland Research Fellowship. M.K. was
supported by NSF Grant No. 2112699 and the Simons Foundation. S.C.H. was
supported by the P. J. E. Peebles Fellowship at Perimeter Institute for
Theoretical Physics and the Horizon Fellowship from Johns Hopkins University.
This research was supported in part by Perimeter Institute for Theoretical
Physics. Research at Perimeter Institute is supported by the Government of
Canada through the Department of Innovation, Science and Economic Development
Canada and by the Province of Ontario through the Ministry of Research,
Innovation and Science. This work was in part carried out at the Advanced
Research Computing at Hopkins (ARCH) core facility (rockfish.jhu.edu), which
is supported by the NSF Grant No. OAC1920103. This work was performed in part
at Aspen Center for Physics, which is supported by NSF Grant No. PHY-2210452.

##  Appendix A Bubble model and the binned optical-depth spectra

The optical depth from a range of redshifts within $z\in[z_{\rm
min}^{\alpha},z_{\rm max}^{\alpha}]$ in the sky satisfies

$\tau_{X}^{\alpha}(\hat{\mathbf{n}},z)=\sigma_{T}n_{p,0}\int_{z^{\alpha}_{\rm
min}}^{z^{\alpha}_{\rm max}}\frac{\differential
z^{\prime}(1+z^{\prime})^{2}}{H(z^{\prime})}x_{e}^{X}(\hat{\mathbf{n}},z^{\prime})\,,$
(37)

where $\alpha$ indicates the redshift bin, and the cross-power spectra can be
defined as

$\begin{split}C_{\ell}^{\tau_{\alpha}\tau_{\beta}}=&\sigma_{T}^{2}n_{\rm
p,0}^{2}\int_{\chi_{\rm min}^{\alpha}}^{\chi_{\rm
max}^{\alpha}}\\!\\!\differential\chi\\!\\!\int_{\chi_{\rm
min}^{\beta}}^{\chi_{\rm
max}^{\beta}}\differential\chi^{\prime}\int\differential
k\,k^{2}(2/\pi)\frac{j_{\ell}(k\chi)}{a(\chi)^{2}}\frac{j_{\ell}(k\chi^{\prime})}{a(\chi^{\prime})^{2}}\sqrt{P^{\rm
H}_{\Delta_{e}\Delta_{e}}(k,\chi)}\sqrt{P^{\rm
H}_{\Delta_{e}\Delta_{e}}(k,\chi^{\prime})}\\\ &+\sigma_{T}^{2}f_{\rm
He}^{2}n_{\rm p,0}^{2}\int_{\chi_{\rm min}^{\alpha}}^{\chi_{\rm
max}^{\alpha}}\\!\\!\differential\chi\\!\\!\int_{\chi_{\rm
min}^{\beta}}^{\chi_{\rm
max}^{\beta}}\differential\chi^{\prime}\int\differential
k\,k^{2}(2/\pi)\frac{j_{\ell}(k\chi)}{a(\chi)^{2}}\frac{j_{\ell}(k\chi^{\prime})}{a(\chi^{\prime})^{2}}\sqrt{P^{\rm
He}_{\Delta_{e}\Delta_{e}}(k,\chi)}\sqrt{P^{\rm
He}_{\Delta_{e}\Delta_{e}}(k,\chi^{\prime})}\\\ &+2\sigma_{T}^{2}f_{\rm
He}n_{\rm p,0}^{2}\int_{\chi_{\rm min}^{\alpha}}^{\chi_{\rm
max}^{\alpha}}\\!\\!\differential\chi\\!\\!\int_{\chi_{\rm
min}^{\beta}}^{\chi_{\rm
max}^{\beta}}\differential\chi^{\prime}\int\differential
k\,k^{2}(2/\pi)\frac{j_{\ell}(k\chi)}{a(\chi)^{2}}\frac{j_{\ell}(k\chi^{\prime})}{a(\chi^{\prime})^{2}}\sqrt{\bar{x}_{e}^{\rm
He-H}(\chi)P(k,\chi)}\sqrt{\bar{x}_{e}^{\rm
He-H}(\chi^{\prime})P(k,\chi^{\prime})}\,,\\\ \end{split}$ (38)

with $\bar{x}_{e}^{\rm He-H}(\chi)\equiv\bar{x}_{e}^{\rm
H}(\chi)\bar{x}_{e}^{\rm He}(\chi)$. For the diagonal redshift bins, Eq. (38)
simplifies to

${C_{\ell}^{\tau_{\alpha}\tau_{\alpha}}=\int_{\chi_{\rm
min}^{\alpha}}^{\chi_{\rm
max}^{\alpha}}\\!\\!\\!\differential\chi\\!\frac{\sigma_{T}^{2}n_{p,0}^{2}}{a^{4}\chi^{2}}P^{\rm
H}_{\Delta_{e}\Delta_{e}}(\chi,k)+\int_{\chi_{\rm min}^{\alpha}}^{\chi_{\rm
max}^{\alpha}}\\!\\!\\!\differential\chi\\!\frac{\sigma_{T}^{2}f_{\rm
He}^{2}n_{p,0}^{2}}{a^{4}\chi^{2}}P^{\rm
He}_{\Delta_{e}\Delta_{e}}(\chi,k)+2\int_{\chi_{\rm min}^{\alpha}}^{\chi_{\rm
max}^{\alpha}}\\!\\!\\!\differential\chi\bar{x}_{e}^{\rm He}\bar{x}_{e}^{\rm
H}\frac{\sigma_{T}^{2}f_{\rm He}n_{p,0}^{2}}{a^{4}\chi^{2}}P(\chi,k)\,,}$ (39)

where $k=\ell/\chi$, and we omitted showing the redshift dependence of mean
ionization fractions for brevity.

The cross power between the galaxy and the ionized electron fluctuations
satisfies

$\begin{split}C_{\ell}^{\tau_{\alpha}g_{\beta}}&=\sigma_{T}n_{\rm
p,0}\int_{\chi_{\rm min}^{\alpha}}^{\chi_{\rm
max}^{\alpha}}\\!\\!\differential\chi\\!\\!\int_{\chi_{\rm
min}^{\beta}}^{\chi_{\rm
max}^{\beta}}\differential\chi^{\prime}\int\differential
k\,k^{2}(2/\pi)\frac{j_{\ell}(k\chi)}{a(z)^{2}}{j_{\ell}(k\chi^{\prime})}P_{\Delta
x_{e}\,g}^{\rm H}(\chi,\chi^{\prime},k)\\\ &+\sigma_{T}f_{\rm He}n_{\rm
p,0}\int_{\chi_{\rm min}^{\alpha}}^{\chi_{\rm
max}^{\alpha}}\\!\\!\differential\chi\\!\\!\int_{\chi_{\rm
min}^{\beta}}^{\chi_{\rm
max}^{\beta}}\differential\chi^{\prime}\int\differential
k\,k^{2}(2/\pi)\frac{j_{\ell}(k\chi)}{a(z)^{2}}{j_{\ell}(k\chi^{\prime})}P_{\Delta
x_{e}\,g}^{\rm He}(\chi,\chi^{\prime},k)\,,\end{split}$ (40)

where we defined

$P_{\Delta{x_{e}}g}(\chi,\chi^{\prime},k)=[\bar{x}_{e}-(1-\bar{x}_{e})\ln(1-\bar{x}_{e})b\,I(k)]b_{g}\,\sqrt{P(\chi,k)}\sqrt{P(\chi^{\prime},k)}\,.$
(41)

Similarly for the diagonal redshift bins, the bin-bin auto-spectrum of the
cross correlation reduces to the simpler form as follows:

${C_{\ell}^{\tau_{\alpha}g_{\alpha}}=\int_{\chi_{\rm min}^{\alpha}}^{\chi_{\rm
max}^{\alpha}}\\!\\!\\!\differential\chi\\!\frac{\sigma_{T}n_{p,0}}{a^{2}\chi^{2}}P^{\rm
H}_{\Delta_{e}g}(\chi,k)+\int_{\chi_{\rm min}^{\alpha}}^{\chi_{\rm
max}^{\alpha}}\\!\\!\\!\differential\chi\\!\frac{\sigma_{T}f_{\rm
He}n_{p,0}}{a^{2}\chi^{2}}P^{\rm He}_{\Delta_{e}g}(\chi,k)}\,,$ (42)

with terms $P^{\rm H/He}_{\Delta x_{e}\,g}(\chi,k)$ defined in Eq. (8).

##  Appendix B Optical depth reconstruction noise

The modulations of temperature and polarization due to the effects described
in Sec. II.5 allows defining multiple estimators for the optical depth from a
combination of temperature $T$ and polarization $E$ and $B$ fields. The
corresponding $\Gamma$ couplings for these combinations are [45]

$\Gamma_{\ell_{1}\ell_{2}\ell}^{TT(\chi_{\alpha})}=(C_{\ell_{1}}^{T_{0}T_{1}^{\alpha}}+C_{\ell_{2}}^{T_{0}T_{1}^{\alpha}})J_{000}^{\ell_{1}\ell_{2}\ell}\,,$
(43)
$\Gamma_{\ell_{1}\ell_{2}\ell}^{EE(\chi_{\alpha})}=\frac{1}{2}(C_{\ell_{1}}^{E_{0}E_{1}^{\alpha}}+C_{\ell_{2}}^{E_{0}E_{1}^{\alpha}})(J_{-220}^{\ell_{1}\ell_{2}\ell}+J_{2-20}^{\ell_{1}\ell_{2}\ell})\,,$
(44)
$\begin{split}\Gamma_{\ell_{1}\ell_{2}\ell}^{TE(\chi_{\alpha})}=\frac{C_{\ell_{1}}^{T_{0}E_{1}^{\alpha}}}{2}(J_{-220}^{\ell_{1}\ell_{2}\ell}+J_{2-20}^{\ell_{1}\ell_{2}\ell})+C_{\ell}^{T_{1}E_{0}^{\alpha}}J_{000}^{\ell_{1}\ell_{2}\ell}\,,\end{split}$
(45)
$\Gamma_{\ell_{1}\ell_{2}\ell}^{TB(\chi_{\alpha})}=\frac{C_{\ell_{1}}^{T_{0}E_{1}^{\alpha}}}{2i}(J_{-220}^{\ell_{1}\ell_{2}\ell}-J_{2-20}^{\ell_{1}\ell_{2}\ell})\,,$
(46)

and

$\Gamma_{\ell_{1}\ell_{2}\ell}^{EB(\chi_{\alpha})}=\frac{C_{\ell_{1}}^{E_{0}E_{1}^{\alpha}}}{2i}(J_{-220}^{\ell_{1}\ell_{2}\ell}-J_{2-20}^{\ell_{1}\ell_{2}\ell})\,,$
(47)

where

$J_{m_{1}m_{2}m}^{\ell_{1}\ell_{2}\ell}=\sqrt{\frac{(2\ell_{1}+1)(2\ell_{2}+1)(2\ell+1)}{4\pi}}\begin{pmatrix}\ell_{1}&\ell_{2}&\ell\\\
m_{1}&m_{2}&m\end{pmatrix}\,.$ (48)

The (biased) noise terms can be found as

$\begin{split}\frac{1}{\tilde{N}_{\ell,\alpha\beta}^{\tau\tau;EB}}={\pi}\\!\\!\int\\!\\!\differential\theta\Big{[}\zeta^{E_{0}E_{1}(\chi_{\alpha},\chi_{\beta})}_{22}(\theta)\zeta^{B}_{22}-\zeta^{E_{0}E_{1}(\chi_{\alpha},\chi_{\beta})}_{2-2}(\theta)\zeta^{B}_{2-2}(\theta)\Big{]}d_{00}^{\ell}(\theta)\,,\end{split}$
(49)

where

$\zeta^{E_{0}E_{1}(\chi_{\alpha})}_{mm^{\prime}}(\theta)=\sum\limits_{\ell}\frac{(2\ell+1)}{4\pi}\frac{C_{\ell}^{E_{0}E_{1}^{\alpha}}C_{\ell}^{E_{0}E_{1}^{\beta}}}{C_{\ell}^{EE}+N_{\ell}^{EE}}d_{mm^{\prime}}^{\ell}(\theta)\,,$
(50)
$\zeta^{B}_{mm^{\prime}}(\theta)=\sum\limits_{\ell}\frac{(2\ell+1)}{4\pi}\frac{1}{C_{\ell}^{BB}+N_{\ell}^{BB}}d_{mm^{\prime}}^{\ell}(\theta)\,,$
(51)

and

$\frac{1}{\tilde{N}_{\ell,\alpha\beta}^{\tau\tau;TT}}={4\pi}\\!\\!\int\\!\\!\differential\theta\\!\left[\zeta^{T_{0}T_{1}(\chi_{\alpha},\chi_{\beta})}_{00}(\theta)\zeta^{T}_{00}+\zeta^{T_{0}T_{1}^{\alpha}}_{00,1}(\theta)\zeta^{T_{0}T_{1}^{\beta}}_{00,1}(\theta)\right]d_{00}^{\ell}(\theta)\,,$
(52)

where

$\zeta^{T_{0}T_{1}(\chi_{\alpha},\chi_{\beta})}_{00}(\theta)=\sum\limits_{\ell}\frac{(2\ell+1)}{4\pi}\frac{C_{\ell}^{T_{0}T_{1}^{\alpha}}C_{\ell}^{T_{0}T_{1}^{\beta}}}{C_{\ell}^{TT}+N_{\ell}^{TT}}d_{00}^{\ell}(\theta)\,,$
(53)
$\zeta^{T}_{00}(\theta)=\sum\limits_{\ell}\frac{(2\ell+1)}{4\pi}\frac{1}{C_{\ell}^{TT}+N_{\ell}^{TT}}d_{00}^{\ell}(\theta)\,,$
(54)
$\zeta^{T_{0}T_{1}(\chi_{\alpha})}_{00,1}(\theta)=\sum\limits_{\ell}\frac{(2\ell+1)}{4\pi}\frac{C_{\ell}^{T_{0}T_{1}^{\alpha}}}{C_{\ell}^{TT}+N_{\ell}^{TT}}d_{00}^{\ell}(\theta)\,,$
(55)

and

$\begin{split}\frac{1}{\tilde{N}_{\ell,\alpha\alpha}^{\tau\tau;EE}}=2{\pi}\\!\\!\int\\!\\!\differential\theta\Big{[}\zeta^{E_{0}E_{1}(\chi_{\alpha},\chi_{\beta})}_{22}(\theta)\zeta^{E}_{22}&\\!+\\!\zeta^{E_{0}E_{1}(\chi_{\alpha},\chi_{\beta})}_{2-2}(\theta)\zeta^{E}_{2-2}(\theta)\\\
&\\!+\\!\zeta^{E_{0}E_{1}(\chi_{\alpha})}_{22,1}(\theta)\zeta^{E_{0}E_{1}(\chi_{\beta})}_{22,1}(\theta)\\!+\\!\zeta^{E_{0}E_{1}(\chi_{\alpha})}_{2-2,1}(\theta)\zeta^{E_{0}E_{1}(\chi_{\beta})}_{2-2,1}(\theta)\Big{]}d_{00}^{\ell}(\theta)\,,\end{split}$
(56)

where

$\zeta^{E_{0}E_{1}(\chi_{\alpha},\chi_{\beta})}_{mm^{\prime}}(\theta)=\sum\limits_{\ell}\frac{(2\ell+1)}{4\pi}\frac{C_{\ell}^{E_{0}E_{1}^{\alpha}}C_{\ell}^{E_{0}E_{1}^{\beta}}}{C_{\ell}^{EE}+N_{\ell}^{EE}}d_{mm^{\prime}}^{\ell}(\theta)\,,$
(57)
$\zeta^{E}_{mm^{\prime}}(\theta)=\sum\limits_{\ell}\frac{(2\ell+1)}{4\pi}\frac{1}{C_{\ell}^{EE}+N_{\ell}^{EE}}d_{mm^{\prime}}^{\ell}(\theta)\,,$
(58)
$\zeta^{E_{0}E_{1}(\chi_{\alpha})}_{mm^{\prime},1}(\theta)=\sum\limits_{\ell}\frac{(2\ell+1)}{4\pi}\frac{C_{\ell}^{E_{0}E_{1}^{\alpha}}}{C_{\ell}^{EE}+N_{\ell}^{EE}}d_{mm^{\prime}}^{\ell}(\theta)\,,$
(59)

and

$\frac{1}{\tilde{N}_{\ell,\alpha\beta}^{\tau\tau;TB}}={\pi}\\!\\!\int\\!\\!\differential\theta\left[\zeta^{T_{0}E_{1}(\chi_{\alpha},\chi_{\beta})}_{22}(\theta)\zeta^{B}_{22}{-}\zeta^{T_{0}E_{1}(\chi_{\alpha},\chi_{\beta})}_{2-2}(\theta)\zeta^{B}_{2-2}(\theta)\right]d_{00}^{\ell}(\theta)\,,$
(60)

where

$\zeta^{T_{0}E_{1}(\chi_{\alpha},\chi_{\beta})}_{mm^{\prime}}(\theta)=\sum\limits_{\ell_{1}}\frac{(2\ell_{1}+1)}{4\pi}\frac{C_{\ell_{1}}^{T_{0}E_{1}^{\alpha}}C_{\ell_{1}}^{T_{0}E_{1}^{\beta}}}{C_{\ell_{1}}^{TT}+N_{\ell_{1}}^{TT}}d_{mm^{\prime}}^{\ell_{1}}(\theta)\,,$
(61)

and finally

$\begin{split}\frac{1}{\tilde{N}_{\ell,\alpha\beta}^{\tau\tau;TE}}={\pi}\\!\\!\int\\!\\!\differential\theta\big{[}\zeta^{T_{0}E_{1}(\chi_{\alpha},\chi_{\beta})}_{22}(\theta)\zeta^{E}_{22}&-\zeta^{T_{0}E_{1}(\chi_{\alpha},\chi_{\beta})}_{2-2}(\theta)\zeta^{E}_{2-2}(\theta)+\zeta^{T_{0}E_{1}(\chi_{\alpha})}_{-20,1}(\theta)\zeta^{T_{0}E_{1}(\chi_{\beta})}_{20,1}(\theta)\\\
&+\zeta^{T_{0}E_{1}(\chi_{\alpha})}_{20,1}(\theta)\zeta^{T_{0}E_{1}(\chi_{\beta})}_{-20,1}(\theta)+2\zeta^{T_{0}E_{1}(\chi_{\alpha},\chi_{\beta})}_{00}(\theta)\zeta^{E}_{00}(\theta)\big{]}d_{00}^{\ell}(\theta)\,,\end{split}$
(62)

where

$\zeta^{T_{0}E_{1}(\chi_{\alpha})}_{mm^{\prime},1}(\theta)=\sum\limits_{\ell}\frac{(2\ell+1)}{4\pi}\frac{C_{\ell}^{T_{0}E_{1}^{\alpha}}}{C_{\ell}^{TT}+N_{\ell}^{TT}}d_{mm^{\prime}}^{\ell}(\theta)\,.$
(63)

Throughout, we assume no tensor fluctuations and set
$({\tilde{N}_{\ell,\alpha\beta}^{\tau\tau;BB}})^{-1}\rightarrow 0$. These
noise terms can be rotated using Eq. (29) to get the noise terms on the
unbiased optical-depth and the minimum-variance noise, which can be calculated
using Eq. (31).

##  Appendix C Improvements beyond ILC-cleaned foregrounds

Our SNR and reionization parameter forecasts in Sec. IV taken into account the
residual foregrounds in the CMB maps after standard ILC cleaning. However,
there is growing evidence that the foreground removal techniques will be
improved in the near future [e.g., 112, 120, 121, 122, 123, 124] with more
advanced methods like constrained-ILC [123] and needlet-based ILC cleaning
[124]. Better removal of frequency-dependent foregrounds improves the
prospects of detecting signatures of helium reionization significantly, as can
be seen from Figs. 7 and 8. To highlight the reach of our method in case
frequency-dependent foregrounds can be removed perfectly from the data, we
show the SNR and parameter forecasts considering only the black-body
contributions to the CMB in Figs. 17–20.

In Figs. 17 and 18, we show the fractional parameter measurement errors
assuming black-body CMB spectra. The constraints on the helium reionization
parameters are observed to improve. Most notably, our findings indicate that
the helium reionization parameters, characterizing the time ($y_{\rm re}^{\rm
He}$) and duration ($\Delta_{y}^{\rm He}$) of this epoch, become detectable
with high significance for $L_{\rm max}>3000$.

Figures 19 and 20 illustrate the fractional errors as functions of the varying
values of $y_{\rm re}^{\rm He}$ and $\Delta_{y}^{\rm He}$, assuming black-body
CMB-S4 and CMB-HD scenarios, respectively. We observe that the fractional
errors improve by a factor of approximately 3 to 5 across all parameters.
Compared to the ILC-cleaned CMB results shown in Figs.12 and13, the dependency
on the varying helium reionization parameter values is similar yet more
pronounced.

Figure 17: Similar to Fig. 12, but assuming black-body CMB spectra instead of
ILC cleaning. Figure 18: Similar to Fig. 13, but assuming black-body CMB
spectra instead of ILC cleaning. Figure 19: Similar to Fig. 14, but assuming
black-body CMB spectra instead of ILC cleaning.

Figure 20: Similar to Fig. 15, but assuming black-body CMB spectra instead of
ILC cleaning.

## References

  * Barkana and Loeb [2001] R. Barkana and A. Loeb, In the beginning: The First sources of light and the reionization of the Universe, Phys. Rept. 349, 125 (2001), arXiv:astro-ph/0010468 .
  * Pritchard and Loeb [2012] J. R. Pritchard and A. Loeb, 21-cm cosmology, Rept. Prog. Phys. 75, 086901 (2012), arXiv:1109.6012 [astro-ph.CO] .
  * Wyithe and Loeb [2003] J. S. B. Wyithe and A. Loeb, Reionization of hydrogen and helium by early stars and quasars, Astrophys. J. 586, 693 (2003), arXiv:astro-ph/0209056 .
  * Madau and Haardt [2015] P. Madau and F. Haardt, Cosmic Reionization after Planck: Could Quasars Do It All?, Astrophys. J. Lett. 813, L8 (2015), arXiv:1507.07678 [astro-ph.CO] .
  * McQuinn [2012] M. McQuinn, Constraints on X-ray Emissions from the Reionization Era, Mon. Not. Roy. Astron. Soc. 426, 1349 (2012), arXiv:1206.1335 [astro-ph.CO] .
  * Worseck  _et al._ [2016] G. Worseck, J. X. Prochaska, J. F. Hennawi, and M. McQuinn, Early and Extended Helium Reionization Over More Than 600 Million Years of Cosmic Time, Astrophys. J. 825, 144 (2016), arXiv:1405.7405 [astro-ph.CO] .
  * Furlanetto and Oh [2008a] S. Furlanetto and S. P. Oh, The History and Morphology of Helium Reionization, Astrophys. J. 681, 1 (2008a), arXiv:0711.1542 [astro-ph] .
  * Worseck  _et al._ [2011] G. Worseck, J. X. Prochaska, M. McQuinn, A. Dall’Aglio, C. Fechner, J. F. Hennawi, D. Reimers, P. Richter, and L. Wisotzki, The End of Helium Reionization at z~2.7 Inferred from Cosmic Variance in HST/COS HeII Lyman Alpha Absorption Spectra, Astrophys. J. Lett. 733, L24 (2011), arXiv:1103.5752 [astro-ph.CO] .
  * Sokasian  _et al._ [2002] A. Sokasian, T. Abel, and L. E. Hernquist, The epoch of helium reionization, Mon. Not. Roy. Astron. Soc. 332, 601 (2002), arXiv:astro-ph/0112297 .
  * Compostella  _et al._ [2013] M. Compostella, S. Cantalupo, and C. Porciani, The imprint of inhomogeneous HeII reionization on the HI and HeII Ly-alpha forest, Mon. Not. Roy. Astron. Soc. 435, 3169 (2013), arXiv:1306.5745 [astro-ph.CO] .
  * Oh  _et al._ [2001] S. P. Oh, Z. Haiman, and M. J. Rees, HeII recombination lines from the first luminous objects, Astrophys. J. 553, 73 (2001), arXiv:astro-ph/0007351 .
  * Furlanetto and Oh [2008b] S. Furlanetto and S. P. Oh, Inhomogeneous Helium Reionization and the Equation of State of the Intergalactic Medium, Astrophys. J. 682, 14 (2008b), arXiv:0711.0751 [astro-ph] .
  * Furlanetto [2009] S. Furlanetto, Fluctuations in the Ionizing Background During and After Helium Reionization, Astrophys. J. 703, 702 (2009), arXiv:0812.3411 [astro-ph] .
  * Dixon and Furlanetto [2009] K. L. Dixon and S. R. Furlanetto, The evolution of the helium-ionizing background at z ~ 2-3, Astrophys. J. 706, 970 (2009), arXiv:0906.4116 [astro-ph.CO] .
  * La Plante  _et al._ [2017] P. La Plante, H. Trac, R. Croft, and R. Cen, Helium Reionization Simulations. II. Signatures of Quasar Activity on the IGM, Astrophys. J. 841, 87 (2017), arXiv:1610.02047 [astro-ph.CO] .
  * Caleb  _et al._ [2019] M. Caleb, C. Flynn, and B. Stappers, Constraining the era of helium reionization using fast radio bursts, Mon. Not. Roy. Astron. Soc. 485, 2281 (2019), arXiv:1902.06981 [astro-ph.HE] .
  * Linder [2020] E. V. Linder, Detecting Helium Reionization with Fast Radio Bursts, Phys. Rev. D 101, 103019 (2020), arXiv:2001.11517 [astro-ph.CO] .
  * Meiksin and Tittley [2012] A. Meiksin and E. R. Tittley, The impact of helium reionization on the structure of the intergalactic medium, Mon. Not. Roy. Astron. Soc. 423, 7 (2012), arXiv:1109.5037 [astro-ph.CO] .
  * Compostella  _et al._ [2014] M. Compostella, S. Cantalupo, and C. Porciani, AGN-driven helium reionization and the incidence of extended He III regions at redshift $z>3$, Mon. Not. Roy. Astron. Soc. 445, 4186 (2014), arXiv:1407.1316 [astro-ph.CO] .
  * Eide  _et al._ [2020] M. B. Eide, B. Ciardi, L. Graziani, P. Busch, Y. Feng, and T. Di Matteo, Large scale simulations of H and He reionization and heating driven by stars and more energetic sources, Mon. Not. Roy. Astron. Soc. 498, 6083 (2020), arXiv:2009.06631 [astro-ph.CO] .
  * Upton Sanderbeck and Bird [2020] P. Upton Sanderbeck and S. Bird, Inhomogeneous He ii reionization in hydrodynamic simulations, Mon. Not. Roy. Astron. Soc. 496, 4372 (2020), arXiv:2002.05733 [astro-ph.CO] .
  * Bhattacharya  _et al._ [2021] M. Bhattacharya, P. Kumar, and E. V. Linder, Fast Radio Burst Dispersion Measure Distribution as a Probe of Helium Reionization, Phys. Rev. D 103, 103526 (2021), arXiv:2010.14530 [astro-ph.CO] .
  * Villasenor  _et al._ [2022] B. Villasenor, B. Robertson, P. Madau, and E. Schneider, Inferring the Thermal History of the Intergalactic Medium from the Properties of the Hydrogen and Helium Ly$\alpha$ Forest, Astrophys. J. 933, 59 (2022), arXiv:2111.00019 [astro-ph.CO] .
  * Meiksin  _et al._ [2010] A. Meiksin, E. R. Tittley, and C. K. Brown, Helium reionization and the thermal proximity effect, Mon. Not. Roy. Astron. Soc. 401, 77 (2010), arXiv:1004.4765 [astro-ph.CO] .
  * Gotberg  _et al._ [2020] Y. Gotberg, S. E. de Mink, M. McQuinn, E. Zapartas, J. H. Groh, and C. Norman, Contribution from stars stripped in binaries to cosmic reionization of hydrogen and helium, Astron. Astrophys. 634, A134 (2020), arXiv:1911.00543 [astro-ph.GA] .
  * La Plante and Trac [2016] P. La Plante and H. Trac, Helium Reionization Simulations. I. Modeling Quasars as Radiation Sources, Astrophys. J. 828, 90 (2016), arXiv:1507.03021 [astro-ph.CO] .
  * Syphers  _et al._ [2011] D. Syphers, S. F. Anderson, W. Zheng, B. Smith, M. Pieri, G. A. Kriss, A. Meiksin, D. P. Schneider, J. M. Shull, and D. G. York, He II Ly\beta Gunn-Peterson Absorption: New HST Observations, and Theoretical Expectations, Astrophys. J. 742, 99 (2011), arXiv:1108.4727 [astro-ph.CO] .
  * Dixon  _et al._ [2014] K. L. Dixon, S. R. Furlanetto, and A. Mesinger, Semi-numeric simulations of helium reionization and the fluctuating radiation background, Mon. Not. Roy. Astron. Soc. 440, 987 (2014), arXiv:1306.1255 [astro-ph.CO] .
  * La Plante  _et al._ [2018] P. La Plante, H. Trac, R. Croft, and R. Cen, Helium Reionization Simulations. III. The Helium Ly$\alpha$ Forest, Astrophys. J. 868, 106 (2018), arXiv:1710.03286 [astro-ph.CO] .
  * Lau  _et al._ [2021] A. W. K. Lau, A. Mitra, M. Shafiee, and G. Smoot, Constraining HeII reionization detection uncertainties via fast radio bursts, New Astron. 89, 101627 (2021), arXiv:2006.11072 [astro-ph.HE] .
  * Hotinli  _et al._ [2023a] S. C. Hotinli, S. Ferraro, G. P. Holder, M. C. Johnson, M. Kamionkowski, and P. La Plante, Probing helium reionization with kinetic Sunyaev-Zel’dovich tomography, Phys. Rev. D 107, 103517 (2023a), arXiv:2207.07660 [astro-ph.CO] .
  * Hotinli [2023] S. C. Hotinli, Cosmological probes of helium reionization, Phys. Rev. D 108, 043528 (2023), arXiv:2212.08004 [astro-ph.CO] .
  * Masters  _et al._ [2012] D. Masters, P. Capak, M. Salvato, F. Civano, B. Mobasher, B. Siana, G. Hasinger, C. D. Impey, T. Nagao, J. R. Trump, H. Ikeda, M. Elvis, and N. Scoville, Evolution of the Quasar Luminosity Function over 3 < z < 5 in the COSMOS Survey Field, ApJ 755, 169 (2012), arXiv:1207.2154 [astro-ph.CO] .
  * Ross  _et al._ [2013] N. P. Ross, I. D. McGreer, M. White, G. T. Richards, A. D. Myers, N. Palanque-Delabrouille, M. A. Strauss, S. F. Anderson, Y. Shen, W. N. Brandt, C. Yèche, M. E. C. Swanson, É. Aubourg, S. Bailey, D. Bizyaev, J. Bovy, H. Brewington, J. Brinkmann, C. DeGraf, T. Di Matteo, G. Ebelke, X. Fan, J. Ge, E. Malanushenko, V. Malanushenko, R. Mandelbaum, C. Maraston, D. Muna, D. Oravetz, K. Pan, I. Pâris, P. Petitjean, K. Schawinski, D. J. Schlegel, D. P. Schneider, J. D. Silverman, A. Simmons, S. Snedden, A. Streblyanska, N. Suzuki, D. H. Weinberg, and D. York, The SDSS-III Baryon Oscillation Spectroscopic Survey: The Quasar Luminosity Function from Data Release Nine, ApJ 773, 14 (2013), arXiv:1210.6389 [astro-ph.CO] .
  * McGreer  _et al._ [2013] I. D. McGreer, L. Jiang, X. Fan, G. T. Richards, M. A. Strauss, N. P. Ross, M. White, Y. Shen, D. P. Schneider, A. D. Myers, W. N. Brandt, C. DeGraf, E. Glikman, J. Ge, and A. Streblyanska, The z = 5 Quasar Luminosity Function from SDSS Stripe 82, ApJ 768, 105 (2013), arXiv:1212.4493 [astro-ph.CO] .
  * McGreer  _et al._ [2018] I. D. McGreer, X. Fan, L. Jiang, and Z. Cai, The Faint End of the z = 5 Quasar Luminosity Function from the CFHTLS, Astron. J. 155, 131 (2018), arXiv:1710.09390 [astro-ph.GA] .
  * Pan  _et al._ [2022] Z. Pan, L. Jiang, X. Fan, J. Wu, and J. Yang, Quasar UV Luminosity Function at 3.5 < z < 5.0 from SDSS Deep Imaging Data, ApJ 928, 172 (2022), arXiv:2112.07801 [astro-ph.GA] .
  * Shen and Ho [2014] Y. Shen and L. C. Ho, The diversity of quasars unified by accretion and orientation, Nature 513, 210 (2014), arXiv:1409.2887 [astro-ph.GA] .
  * Hopkins  _et al._ [2007] P. F. Hopkins, A. Lidz, L. Hernquist, A. L. Coil, A. D. Myers, T. J. Cox, and D. N. Spergel, The Co-Formation of Spheroids and Quasars Traced in their Clustering, Astrophys. J. 662, 110 (2007), arXiv:astro-ph/0611792 .
  * Schmidt  _et al._ [2017] T. M. Schmidt, G. Worseck, J. F. Hennawi, J. X. Prochaska, and N. H. M. Crighton, Statistical Detection of the He II Transverse Proximity Effect: Evidence for Sustained Quasar Activity for >25 Million Years, ApJ 847, 81 (2017), arXiv:1701.08769 [astro-ph.GA] .
  * Inayoshi  _et al._ [2020] K. Inayoshi, E. Visbal, and Z. Haiman, The Assembly of the First Massive Black Holes, Ann. Rev. Astron. Astrophys. 58, 27 (2020), arXiv:1911.05791 [astro-ph.GA] .
  * Sunyaev and Zeldovich [1972] R. A. Sunyaev and Y. B. Zeldovich, The Observations of Relic Radiation as a Test of the Nature of X-Ray Radiation from the Clusters of Galaxies, Comments on Astrophysics and Space Physics 4, 173 (1972).
  * Dore  _et al._ [2007] O. Dore, G. Holder, M. Alvarez, I. T. Iliev, G. Mellema, U.-L. Pen, and P. R. Shapiro, The Signature of Patchy Reionization in the Polarization Anisotropy of the CMB, Phys. Rev. D 76, 043002 (2007), arXiv:astro-ph/0701784 .
  * Holder  _et al._ [2007] G. P. Holder, I. T. Iliev, and G. Mellema, Reconstructing the Thomson Optical Depth due to Patchy Reionization with 21 cm Fluctuation Maps, ApJ 663, L1 (2007), arXiv:astro-ph/0609689 [astro-ph] .
  * Dvorkin and Smith [2009] C. Dvorkin and K. M. Smith, Reconstructing Patchy Reionization from the Cosmic Microwave Background, Phys. Rev. D 79, 043003 (2009), arXiv:0812.1566 [astro-ph] .
  * Dvorkin  _et al._ [2009] C. Dvorkin, W. Hu, and K. M. Smith, B-mode CMB Polarization from Patchy Screening during Reionization, Phys. Rev. D 79, 107302 (2009), arXiv:0902.4413 [astro-ph.CO] .
  * Feng and Holder [2019] C. Feng and G. Holder, Searching for patchy reionization from cosmic microwave background with hybrid quadratic estimators, Phys. Rev. D 99, 123502 (2019), arXiv:1808.01592 [astro-ph.CO] .
  * Abazajian  _et al._ [2016a] K. N. Abazajian _et al._ (CMB-S4), CMB-S4 Science Book, First Edition, (2016a), arXiv:1610.02743 [astro-ph.CO] .
  * Alvarez  _et al._ [2021] M. A. Alvarez, S. Ferraro, J. C. Hill, R. Hložek, and M. Ikape, Mitigating the optical depth degeneracy using the kinematic Sunyaev-Zel’dovich effect with CMB-S4, Phys. Rev. D 103, 063518 (2021), arXiv:2006.06594 [astro-ph.CO] .
  * Abazajian  _et al._ [2022] K. Abazajian _et al._ , Snowmass 2021 CMB-S4 White Paper, in  _2022 Snowmass Summer Study_ (2022) arXiv:2203.08024 [astro-ph.CO] .
  * Sehgal  _et al._ [2019] N. Sehgal _et al._ , CMB-HD: An Ultra-Deep, High-Resolution Millimeter-Wave Survey Over Half the Sky, (2019), arXiv:1906.10134 [astro-ph.CO] .
  * Aiola  _et al._ [2022] S. Aiola _et al._ (CMB-HD), Snowmass2021 CMB-HD White Paper, (2022), arXiv:2203.05728 [astro-ph.CO] .
  * MacInnis  _et al._ [2023] A. MacInnis, N. Sehgal, and M. Rothermel, Cosmological Parameter Forecasts for a CMB-HD Survey, (2023), arXiv:2309.03021 [astro-ph.CO] .
  * Ivezić  _et al._ [2019] v. Ivezić _et al._ (LSST), LSST: from Science Drivers to Reference Design and Anticipated Data Products, Astrophys. J. 873, 111 (2019), arXiv:0805.2366 [astro-ph] .
  * The LSST Dark Energy Science Collaboration  _et al._ [2018] The LSST Dark Energy Science Collaboration, R. Mandelbaum, T. Eifler, R. Hložek, T. Collett, E. Gawiser, and o. Scolnic, The LSST Dark Energy Science Collaboration (DESC) Science Requirements Document, arXiv e-prints , arXiv:1809.01669 (2018), arXiv:1809.01669 [astro-ph.CO] .
  * Schlegel  _et al._ [2019] D. J. Schlegel _et al._ , Astro2020 APC White Paper: The MegaMapper: a z $>$ 2 Spectroscopic Instrument for the Study of Inflation and Dark Energy, (2019), arXiv:1907.11171 [astro-ph.IM] .
  * Lewis and Challinor [2006] A. Lewis and A. Challinor, Weak gravitational lensing of the cmb, Phys. Rept. 429, 1 (2006), arXiv:astro-ph/0601594 [astro-ph] .
  * Cayuso  _et al._ [2023] J. Cayuso, R. Bloch, S. C. Hotinli, M. C. Johnson, and F. McCarthy, Velocity reconstruction with the cosmic microwave background and galaxy surveys, JCAP 02, 051, arXiv:2111.11526 [astro-ph.CO] .
  * Smith  _et al._ [2018] K. M. Smith, M. S. Madhavacheril, M. Münchmeyer, S. Ferraro, U. Giri, and M. C. Johnson, KSZ tomography and the bispectrum, (2018), arXiv:1810.13423 [astro-ph.CO] .
  * Hotinli  _et al._ [2019] S. C. Hotinli, J. Meyers, N. Dalal, A. H. Jaffe, M. C. Johnson, J. B. Mertens, M. Münchmeyer, K. M. Smith, and A. van Engelen, Transverse Velocities with the Moving Lens Effect, Phys. Rev. Lett. 123, 061301 (2019), arXiv:1812.03167 [astro-ph.CO] .
  * Anil Kumar  _et al._ [2022] N. Anil Kumar, G. Sato-Polito, M. Kamionkowski, and S. C. Hotinli, Primordial trispectrum from kinetic Sunyaev-Zel’dovich tomography, Phys. Rev. D 106, 063533 (2022), arXiv:2205.03423 [astro-ph.CO] .
  * Kumar  _et al._ [2023] N. A. Kumar, S. C. Hotinli, and M. Kamionkowski, Uncorrelated compensated isocurvature perturbations from kinetic Sunyaev-Zeldovich tomography, Phys. Rev. D 107, 043504 (2023), arXiv:2208.02829 [astro-ph.CO] .
  * Hotinli  _et al._ [2021] S. C. Hotinli, K. M. Smith, M. S. Madhavacheril, and M. Kamionkowski, Cosmology with the moving lens effect, Phys. Rev. D 104, 083529 (2021), arXiv:2108.02207 [astro-ph.CO] .
  * Hotinli and Johnson [2022] S. C. Hotinli and M. C. Johnson, Reconstructing large scales at cosmic dawn, Phys. Rev. D 105, 063522 (2022), arXiv:2012.09851 [astro-ph.CO] .
  * Münchmeyer  _et al._ [2019] M. Münchmeyer, M. S. Madhavacheril, S. Ferraro, M. C. Johnson, and K. M. Smith, Constraining local non-Gaussianities with kinetic Sunyaev-Zel’dovich tomography, Phys. Rev. D 100, 083508 (2019), arXiv:1810.13424 [astro-ph.CO] .
  * Deutsch  _et al._ [2018a] A.-S. Deutsch, E. Dimastrogiovanni, M. C. Johnson, M. Münchmeyer, and A. Terrana, Reconstruction of the remote dipole and quadrupole fields from the kinetic Sunyaev Zel’dovich and polarized Sunyaev Zel’dovich effects, Phys. Rev. D 98, 123501 (2018a), arXiv:1707.08129 [astro-ph.CO] .
  * Hotinli  _et al._ [2023b] S. C. Hotinli, E. Pierpaoli, S. Ferraro, and K. Smith, Transverse velocities and matter gradient correlations: A new signal and a new challenge to moving-lens analyses, Phys. Rev. D 108, 083508 (2023b), arXiv:2305.15462 [astro-ph.CO] .
  * Hotinli  _et al._ [2022a] S. C. Hotinli, G. P. Holder, M. C. Johnson, and M. Kamionkowski, Cosmology from the kinetic polarized Sunyaev Zel’dovich effect, JCAP 10, 026, arXiv:2204.12503 [astro-ph.CO] .
  * Lee  _et al._ [2022] N. Lee, S. C. Hotinli, and M. Kamionkowski, Probing cosmic birefringence with polarized Sunyaev-Zel’dovich tomography, Phys. Rev. D 106, 083518 (2022), arXiv:2207.05687 [astro-ph.CO] .
  * Namikawa and Obata [2023] T. Namikawa and I. Obata, Cosmic birefringence tomography with polarized Sunyaev-Zel’dovich effect, Phys. Rev. D 108, 083510 (2023), arXiv:2306.08875 [astro-ph.CO] .
  * Deutsch  _et al._ [2018b] A.-S. Deutsch, M. C. Johnson, M. Münchmeyer, and A. Terrana, Polarized Sunyaev Zel’dovich tomography, JCAP 04, 034, arXiv:1705.08907 [astro-ph.CO] .
  * Meerburg  _et al._ [2013] P. D. Meerburg, C. Dvorkin, and D. N. Spergel, Probing Patchy Reionization through $\tau$-21 cm Correlation Statistics, Astrophys. J. 779, 124 (2013), arXiv:1303.3887 [astro-ph.CO] .
  * Roy  _et al._ [2020] A. Roy, A. Lapi, D. Spergel, S. Basak, and C. Baccigalupi, Detectability of the $\tau_{es}$$-$21cm cross-correlation: a tomographic probe of patchy reionization, JCAP 03, 062, arXiv:1904.02637 [astro-ph.CO] .
  * Orlando  _et al._ [2023] G. Orlando, T. Flöss, P. D. Meerburg, and J. Silk, Local non-Gaussianities from cross-correlations between the CMB and 21-cm, (2023), arXiv:2307.15046 [astro-ph.CO] .
  * Sailer  _et al._ [2021] N. Sailer, E. Castorina, S. Ferraro, and M. White, Cosmology at high redshift — a probe of fundamental physics, JCAP 12 (12), 049, arXiv:2106.09713 [astro-ph.CO] .
  * Becker  _et al._ [2011] G. D. Becker, J. S. Bolton, M. G. Haehnelt, and W. L. W. Sargent, Detection of extended He II reionization in the temperature evolution of the intergalactic medium, MNRAS 410, 1096 (2011), arXiv:1008.2622 [astro-ph.CO] .