# Maximum likelihood kinetic Sunyaev Zel’dovich velocity reconstruction

Dagoberto Contreras1 Fiona McCarthy2,3,4 Matthew C. Johnson1,3 1Department of
Physics and Astronomy, York University, Toronto, Ontario, M3J 1P3, Canada
2Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue,
New York, NY 10010 USA  3Perimeter Institute for Theoretical Physics,
Waterloo, Ontario N2L 2Y5, Canada  4Department of Physics and Astronomy,
University of Waterloo, Waterloo, Ontario, N2L 3G1, Canada

###### Abstract

The kinetic Sunyaev Zel’dovich (kSZ) effect, cosmic microwave background (CMB)
temperature anisotropies induced by the scattering of CMB photons from free
electrons, will be measured by near-term CMB experiments at high significance.
By combining CMB temperature anisotropies with a tracer of structure, such as
a galaxy redshift survey, previous literature introduced a number of
techniques to reconstruct the radial velocity field. This reconstructed radial
velocity field encapsulates the majority of the cosmological information
contained in the kSZ temperature anisotropies, and can provide powerful new
tests of the standard cosmological model and theories beyond it. In this
paper, we introduce a new estimator for the radial velocity field based on a
coarse-grained maximum likelihood fit for the kSZ component of the temperature
anisotropies, given a tracer of the optical depth. We demonstrate that this
maximum likelihood estimator yields a higher fidelity reconstruction than
existing quadratic estimators in the low-noise and high-resolution regime
targeted by upcoming CMB experiments. We describe implementations of the
maximum likelihood estimator in both harmonic space and map space, using
either direct measurements of the optical depth field or a galaxy survey as a
tracer. We comment briefly on the impact of biases introduced by imperfect
reconstruction of the optical depth field.

##  I Introduction

Future cosmic microwave background (CMB) experiments such as Simons
Observatory Aguirre et al. (2019), CMB-S4 Abazajian et al. (2016), and CMB-HD
Sehgal et al. (2019) will probe the high-resolution, low-noise frontier where
secondary CMB anisotropies dominate. Secondary anisotropies arise due to the
electromagnetic (Sunyaev Zel’dovich effect) or gravitational (lensing)
scattering of CMB photons from structure, and encode a significant amount of
astrophysical and cosmological information. How to access this information is
an active area of investigation. In this paper, we focus on the kinetic
Sunyaev Zel’dovich (kSZ) effect, blackbody temperature anisotropies induced by
the scattering of CMB photons from free electrons in bulk motion along the
line of sight Sunyaev and Zeldovich (1980). We introduce a new estimator for
the radial velocity field based on a coarse-grained maximum likelihood
solution constructed from the CMB temperature anisotropies and a tracer of the
optical depth.

Previous literature introduced a set of quadratic estimators (QEs) for the
radial velocity field Deutsch et al. (2018a); Smith et al. (2018); Cayuso et
al. (2021) based on the CMB temperature anisotropies and a galaxy survey as a
tracer of the optical depth111More precisely, one reconstructs the remote
dipole field: the CMB dipole seen in the local rest frame of the scatterer,
projected along the line of sight. The largest contribution to the remote
dipole field is given by the radial peculiar velocity, but there are
additional ’primordial’ contributions; see Refs. Terrana et al. (2017);
Deutsch et al. (2018a) for discussion.. Fisher forecasts Deutsch et al.
(2018a); Smith et al. (2018); Cayuso et al. (2021), the analysis of N-body
simulations Cayuso et al. (2018); Giri and Smith (2020), and the analysis of
Gaussian mocks containing various foregrounds and systematics Cayuso et al.
(2021) indicate that these QEs will yield high fidelity reconstructions of the
radial velocity field with near-term CMB and galaxy surveys. The QE
reconstructs the radial velocity field from statistical anisotropies in the
cross-power spectrum of the CMB temperature and a redshift-binned galaxy
survey. The estimator relies on the assumption that all the fields involved
are Gaussian, and is based on the theoretical ensemble average spectra, which
are used to find the appropriate weights for the estimators to be unbiased and
have the minimum variance. However, these properties of the QE lead to
significant signal-to-noise penalties and biases in the high signal-to-noise
regime. For example, the non-Gaussian nature of the optical depth and galaxy
density on small scales leads to the so-called $N^{(3/2)}$ bias which can
exceed the Gaussian reconstruction noise in the high signal-to-noise regime
Giri and Smith (2020). For sufficiently low noise and high resolution CMB
experiments, and in the absence of foregrounds, the dominant contribution to
the observed temperature anisotropies on small angular scales will be the kSZ
effect itself. In this limit, the QE is limited by sample variance, since it
relies on searching for statistical anisotropies beyond those expected in the
ensemble average. However, because one is able to accurately measure the
detailed kSZ temperature anisotropies in this limit, it should be possible to
do better by using the properties of our realization rather than the ensemble
statistics (e.g. to use a form a sample variance cancellation).

To overcome these two shortcomings of the QE, we introduce a new estimator for
the radial velocity field based on a coarse-grained maximum likelihood
condition. This maximum likelihood (MaxL) estimator does not rely on the
assumption that the underlying fields are Gaussian (hence we do not expect the
analogue of the $N^{(3/2)}$ bias), and relies on the realization of the
underlying fields rather than their statistical properties (and hence should
incorporate a form of sample variance cancellation). We therefore expect the
MaxL estimator to have advantages over the QE in the high signal to noise
regime, where the kSZ effect is signal dominated (e.g. over the primary CMB,
instrumental noise, and foregrounds) in the CMB on small angular scales.
Because it relies on the properties of a realization, forecasting with the
MaxL estimator requires simulations, unlike the QE, whose statistical
properties can be calculated analytically. We create properly correlated
random Gaussian realizations of the underlying fields and make direct
comparisons between the reconstruction of the QE as implemented previously
Cayuso et al. (2021) and the MaxL estimator. Because the realizations are
purely Gaussian, we are unable to demonstrate the absence of a non-Gaussian
bias for the MaxL estimator, and therefore it will be important in future work
to compare the results of applying the MaxL estimator and QE to N-body
simulations.

A key ingredient of the MaxL estimator is an estimate of the optical depth
field. In the absence of a direct measurement (e.g. by using a large number of
Fast Radio Bursts Madhavacheril et al. (2019)), we must infer the optical
depth from a tracer of large scale structure, such as a galaxy survey. This
inference will be imperfect, leading to a multiplicative optical depth bias on
the reconstructed velocities Battaglia (2016); Smith et al. (2018);
Madhavacheril et al. (2019); Hotinli et al. (2021). Further multiplicative
biases will be introduced by redshift errors in the galaxy survey and various
systematics in the galaxy and CMB surveys. These multiplicative biases on the
QE are described in detail in Ref. Cayuso et al. (2021), and in this paper we
describe their nature using a simplistic estimator for the optical depth
field. Our implementation of the MaxL estimator can likely be improved upon
through the use of more accurate methods to estimate the optical depth field,
e.g. using machine learning techniques trained on accurate hydrodynamical
simulations (see e.g. Thiele et al. (2020)).

The plan of the paper is as follows. In Sec. II we provide some intuition and
motivation for the MaxL estimator. In Sec. III we review our model for the kSZ
component of the temperature anisotropies, which is largely a review of Ref.
Cayuso et al. (2021). In Sec. IV, we derive an implementation of the maximum
likelihood estimator (given the optical depth field) for the radial velocity
field in both harmonic space and map space. In Sec. V we compare the
performance of the MaxL estimator and QE on a set of correlated Gaussian mock
datasets. In Sec. VI we discuss inference of the optical depth field from a
galaxy redshift survey and the associated optical depth bias on the
reconstructed velocity; we conclude in Sec. VII.

####  Some notation

Throughout the paper, we refer to objects defined on a sphere in both in
“harmonic” and “map” space. We label the multipole moments of fields in
harmonic space by the pairs $(\ell,m),(\ell^{\prime},m^{\prime})\ldots$ and of
the pixels in map space by the letters $i,j,\ldots$. For coarse-grained
quantities, we use capital letters $(L,M),(L^{\prime},M^{\prime}),\ldots$ or
$I,J\ldots$. We denote operators in boldface, such that $\tau$ is a field but
$\boldsymbol{\tau}$ is an operator.

##  II Intuition and motivation

In this section we give a heuristic description of the QE that recovers the
radial velocity field from a temperature and optical depth measurement, and
motivate a new estimator that can avoid some of the QE reconstruction noise in
the high-signal-to-noise regime of velocity reconstruction.

The kSZ contribution to the CMB temperature $\Theta$ is of the form $\tau v$,
where $v$ is the radial velocity field and $\tau$ is the optical depth field,
such that

$\displaystyle\Theta=$ $\displaystyle\ \Theta^{\rm pCMB}+\Theta^{\rm kSZ},$
$\displaystyle\sim$ $\displaystyle\ \Theta^{\rm pCMB}+\tau v,$ (1)

where $\Theta^{\rm pCMB}$ describes all contributions to the temperature
anisotropies that are not kSZ (this includes primary CMB contributions as well
as other foreground contributions and instrumental noise). Here, we employ a
toy model where kSZ is sourced by a two dimensional field; later we include
line of sight effects. The goal is to use a measurement of $\Theta$ and $\tau$
to recover $v$. The QE recovers the velocity from this combination and their
statistics:

$\hat{v}^{\rm QE}\sim\frac{1}{{A}}\Theta\ \tau,$ (2)

where $A$ is a theoretically calculated filter that can be thought of as
$\sim\left<\tau^{2}\right>$; it is calculated from the ensemble statistics of
$\tau$. The QE differs from the true velocity by

$\displaystyle\hat{v}^{\rm QE}-v\sim$
$\displaystyle\frac{1}{\left<\tau^{2}\right>}\Theta\ \tau-v$
$\displaystyle\sim$ $\displaystyle\
v\left(\frac{\tau^{2}}{\left<\tau^{2}\right>}-1\right)+\frac{\Theta^{\rm
pCMB}}{\left<\tau^{2}\right>}.$ (3)

In particular, there is a contribution to the residual sourced by the
difference of the realization of $\tau^{2}$ from its theoretically expected
variance. This contribution persists even when $\Theta^{\rm pCMB}$ is
vanishingly small.

However, in the limit where the $\Theta^{\rm kSZ}$ contributions dominate over
the non-kSZ contributions, we can avoid this contribution to the
reconstruction noise by neglecting $\Theta^{\rm{pCMB}}$ and simply “inverting”
Eq. (1):

$\hat{v}\sim\frac{\Theta}{\tau}.$ (4)

In this case, the “reconstruction noise” only has contributions from the non-
kSZ contributions to the temperature:

$\hat{v}-v\sim\frac{\Theta^{\rm pCMB}}{\tau}.$ (5)

This allows us to make a lower-noise estimate of $v$ by utilizing angular
scales where the kSZ temperature dominates over the other contributions to the
CMB temperature, similar to what is done in the gradient-inversion approach to
CMB lensing reconstruction Hadzhiyska et al. (2019).

In the remainder of this Paper, we develop this idea more rigorously by
deriving a new estimator from the maximization of a likelihood, and compare
explicitly the performance of the “MaxL” estimator with that of the QE. Note
that further complication arises from the fact that the kSZ temperature is not
truly sourced by the product of one velocity and one electron density field,
as in Eq. (1), but by the integral of this product along the line of sight. We
will deal with this explicitly by “redshift-binning” the velocity and optical
depth fields, with the end-goal being to reconstruct the redshift-binned
velocity field.

##  III The kinetic Sunyaev Zel’dovich effect

The observed kSZ temperature anisotropies in the direction $\hat{n}$ are given
by the line-of-sight integral over radial comoving distance $\chi$

$\Theta^{\rm kSZ}(\hat{n})=\int_{0}^{\chi_{\rm max}}d\chi\
\dot{\tau}(\hat{n},\chi)v(\hat{n},\chi).$ (6)

Here, $v(\hat{n},\chi)$ is the radial velocity field and the differential
optical depth $\dot{\tau}(\hat{n},\chi)$ is

$\dot{\tau}(\hat{n},\chi)=-\sigma_{T}a(\chi)\bar{n}_{e}(\chi)\left(1+\delta_{e}(\hat{n},\chi)\right),$
(7)

with $\sigma_{T}$ the Thompson cross section, $a(\chi)$ the scale factor,
$\bar{n}_{e}(\chi)$ the cosmological mean electron density, and
$\delta_{e}(\hat{n},\chi)$ the electron overdensity field. There are
contributions to the kSZ signal from during and after the epoch of
reionization. We are primarily concerned with the late-time kSZ effect in this
paper222See Ref. Hotinli and Johnson (2020) for a discussion of kSZ velocity
reconstruction during the epoch of reionization., and so in Eq. (6) we
introduce a cutoff $\chi_{\rm max}$ on the line of sight integral at some time
after reionization (which also justifies our neglect of the optical depth
factor $e^{-\tau}$ in this equation).

Following Ref. Cayuso et al. (2021), we decompose fields along the radial
direction in the basis of Haar wavelets:

$v(\hat{n},\chi)=\sum_{s=0}^{\infty}v^{s}(\hat{n})h^{s}(\chi),$ (8)

where $h^{s}$ are the Haar wavelet basis functions; these form a complete
basis on the interval $0\leq\chi\leq\chi_{\rm max}$. The Haar expansion has
the additional useful property that if we truncate at a fixed index $N$, the
sum in the Haar basis is simply related to a sum over radial tophat bin-
averages over bins of equal comoving size $\Delta\chi=\chi_{\rm max}/N$:

$v(\hat{n},\chi)=\sum_{\alpha=0}^{N-1}v^{\alpha}(\hat{n})\Pi^{\alpha}(\chi)+\sum_{s=N}^{\infty}v^{s}(\hat{n})h^{s}(\chi),$
(9)

where $\Pi^{\alpha}$ are radial top-hat bins. We will reserve Greek
superscripts $\alpha,\beta,\gamma,\ldots$ to denote the tophat bin-averages.
The bin averages $v^{\alpha}$ give a coarse-grained description of the field
on the past light cone, and the remaining sum over Haar moments gives the
residual fine-grained information. We can therefore write the line of sight
integral in Eq. (6) as the sum:

$\Theta^{\rm
kSZ}(\hat{n})=\sum_{s=0}^{\infty}\tau^{s}(\hat{n})v^{s}(\hat{n})=\sum_{\alpha=0}^{N-1}\tau^{\alpha}(\hat{n})v^{\alpha}(\hat{n})+\sum_{s=N}^{\infty}\tau^{s}(\hat{n})v^{s}(\hat{n}),$
(10)

where we have defined $\tau^{\alpha}(\hat{n})\equiv\Delta\chi\
\dot{\tau}^{\alpha}(\hat{n})$. Below, we approximate the continuum limit by
truncating at fixed $N$. For the contribution to kSZ from the redshift range
$0<z<5$ most of the kSZ autopower is captured for $N=512$ as shown in Ref.
Cayuso et al. (2021). However, in practice we will have to truncate the sum at
a somewhat lower value of $N$ for the application of the estimators derived
below.

It will be useful to think of the optical depth as a set of $N$ operators that
map the bin-averaged radial velocity field to the observed temperature
anisotropies, in which case we can write the sky-basis-independent equation
truncated at fixed $N$:

$\Theta^{\rm kSZ}\simeq\sum_{\alpha=0}^{N-1}\boldsymbol{\tau}^{\alpha}\cdot
v^{\alpha}.$ (11)

Here, just as any field defined on a sphere, $\Theta^{\rm kSZ}$ can be
represented either in map space or harmonic space. In map space, $\Theta_{i}$
has dimension equal to the number of pixels $N_{\rm pix}$ at some resolution
(e.g. using the Healpix pixelization scheme Górski et al. (2005));
$v_{i}^{\alpha}$ is a set of $N$ fields each with dimension equal to
$N_{\rm{pix}}$; $\boldsymbol{\tau}_{ij}^{\alpha}$ is a set of $N$ operators
each with dimension equal to $N_{\rm{pix}}\times N_{\rm{pix}}$. In harmonic
space, $\Theta_{\ell m}$ has dimension equal to the number of modes required
to describe the field at the appropriate resolution (which should in principle
be equal to $N_{\rm{pix}}$, but in practice is somewhat larger); we denote
this number as $N_{\rm harm}$. In Appendix A, we present the form of
$\boldsymbol{\tau}$ in each representation, as well as a mixed harmonic-map
space representation.

##  IV Maximum likelihood estimator

The dominant blackbody components of the observed CMB temperature anisotropies
are the lensed primary CMB and the late-time and reionization kSZ. Neglecting
the reionization kSZ, the blackbody component of the temperature anisotropy is
given by:

$\Theta=\Theta^{\rm pCMB}+\sum_{s=0}^{\infty}{\tau}^{s}\cdot v^{s}.$ (12)

Here $\Theta^{\rm pCMB}$ are the primary CMB anisotropies in map or harmonic
space, which we assume purely Gaussian with covariance
$\boldsymbol{C}^{\mathrm{pCMB}}$. The likelihood function $\mathcal{L}$ for
the primary CMB is then

$\ln\mathcal{L}=(\Theta-\sum_{s=0}^{\infty}\boldsymbol{\tau}^{s}\cdot
v^{s})^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}(\Theta-\sum_{s=0}^{\infty}\boldsymbol{\tau}^{s}\cdot
v^{s})+\ln{\rm det}\ \boldsymbol{C}^{\mathrm{pCMB}}$ (13)

Extremizing $\mathcal{L}$ with respect to $(v^{q})^{\dagger}$:

$\displaystyle\frac{\delta}{\delta(v^{q})^{\dagger}}\ln\mathcal{L}=$
$\displaystyle(\boldsymbol{\tau}^{q})^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}$
$\displaystyle(\Theta-\sum_{s=0}^{\infty}\boldsymbol{\tau}^{s}\cdot v^{s})=0$
$\displaystyle\implies\ \ $
$\displaystyle(\boldsymbol{\tau}^{q})^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\Theta$
$\displaystyle=\sum_{s=0}^{\infty}\left[(\boldsymbol{\tau}^{q})^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}^{s}\right]\cdot
v^{s}$ (14)
$\displaystyle(\boldsymbol{\tau}^{\alpha})^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\Theta$
$\displaystyle\simeq\sum_{\beta=0}^{N-1}\left[(\boldsymbol{\tau}^{\alpha})^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}^{\beta}\right]\cdot
v^{\beta}.$

In the last line we truncate the radial sum at $N$ and work with the bin
coefficients (denoted by the Greek indices). Note that $\Theta$ and
$\boldsymbol{\tau}$ are fixed, and our goal is to find the set of $v^{\alpha}$
that maximize the likelihood. Organizing both $v^{\beta}$ and
$(\boldsymbol{\tau}^{\alpha})^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\Theta$
into $\left(N_{\rm pix}\times N\right)$-dimensional vectors containing both
the bin and pixel information, and organizing
$\left[(\boldsymbol{\tau}^{\alpha})^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}^{\beta}\right]$
into an $(N_{\rm pix}N)\times(N_{\rm pix}N)$-dimensional operator, Eq. (IV) is
a linear equation for the velocity field:

$\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\Theta\right]=\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}\right]\cdot
v.$ (15)

This equation will have a solution when the operator
$\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}$
has a well-defined inverse. In general, this isn’t the case, since there will
be many highly oscillatory functions $v$ that nearly satisfy the maximum
likelihood condition and pixels at which the operator (whose eigenbasis is in
map space) is nearly zero (indicating small eigenvalues). However, we can
first perform the matrix products in
$\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}$
at full resolution, and then coarse-grain this operator so that it is
invertible. The maximum-likelihood condition is then

$\displaystyle
W\left(\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\Theta\right]\right)=W\left(\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}\right]\right)\cdot
W\left(v\right)$ $\displaystyle\implies W\left(v\right)=$
$\displaystyle\left[W\left(\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}\right]\right)\right]^{-1}W\left(\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\Theta\right]\right)$
(16)

where $W(\cdot)$ is the coarse-graining procedure; we demonstrate explicitly
below that it is appropriate to separately coarse-grain
$\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}\right]$
and $v$. The proposal is that with Eq. (16), we can reconstruct the coarse-
grained velocity $W(v)$. We now separately examine if there is a useful
solution to the maximum likelihood condition in harmonic or map space for the
velocity field.

###  IV.1 Harmonic space estimator

Before coarse graining, the maximum likelihood condition in harmonic space is

$\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\Theta\right]_{\ell
m}^{\alpha}=\sum_{\ell^{\prime}m^{\prime};\
\beta}\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}\right]^{\alpha\beta}_{\ell
m;\ell^{\prime}m^{\prime}}v^{\beta}_{\ell^{\prime}m^{\prime}}.$ (17)

In Appendix A we evaluate the components of the operators and vectors in Eq.
(17) explicitly, obtaining:

$\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}\right]^{\alpha\beta}_{\ell
m;\ell^{\prime}m^{\prime}}=\int d\Omega\
\tau^{\alpha}(\hat{n})\bar{\tau}^{\beta}(\hat{n})Y^{*}_{\ell
m}(\hat{n})Y_{\ell^{\prime}m^{\prime}}(\hat{n}),$ (18)

and

$\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\Theta\right]_{\ell
m}^{\beta}=\int d\Omega\ \tau^{\beta}(\hat{n})\bar{\Theta}(\hat{n})Y^{*}_{\ell
m}(\hat{n}),$ (19)

where the bar denotes inverse-variance filtering using
$C_{\ell}^{\mathrm{pCMB}}$:

$\bar{A}(\hat{n})\equiv\sum_{\ell m}\frac{A_{\ell
m}}{C_{\ell}^{\mathrm{pCMB}}}Y_{\ell m}(\hat{n}).$ (20)

The coarse-graining operation required to make the operator
$\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}\right]_{\ell
m;\ell^{\prime}m^{\prime}}^{\alpha\beta}$ invertible is to cut off the sum at
a maximum multipole $L_{\mathrm{max}}$, which will be the maximum multipole of
the $v$ field that we recover; for more details, see Appendix B. After this
coarse-graining, the maximum-likelihood solution for the velocity is:

$\hat{v}_{LM}^{\alpha}=\sum_{L^{\prime}M^{\prime}}^{L_{\rm{max}}}\sum_{\beta=0}^{N-1}\left([\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{-1}\right)^{\alpha\beta}_{LM;L^{\prime}M^{\prime}}(\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\Theta)_{L^{\prime}M^{\prime}}^{\beta},$
(21)

where the $L_{\mathrm{max}}$ cut-off should be done after the calculation of
$\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}$
but before its inversion. The maximum values $N$ and $L_{\rm{max}}$ must be
empirically determined from the numerical inversion of
$\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}\right]_{\ell
m;\ell^{\prime}m^{\prime}}^{\alpha\beta}$ for each realization. A crude
approach to finding an approximate solution to the maximum likelihood
condition Eq. (17) is to substitute with the ensemble-average operator over
realizations of the optical depth, which is explicitly invertible. This
approach is presented in Appendix C, where we demonstrate that it yields an
expression equivalent to the QE in the low signal-to-noise regime. We note
that Eq. (21) is nearly identical to an estimator that was successfully used
to reconstruct our peculiar velocity with respect the CMB rest frame using the
thermal Sunyaev Zel’dovich effect Akrami et al. (2020).

We now compute the mean and variance of the MaxL estimator. The estimator Eq.
(21) is based on a fixed realization of the optical depth, so the statistics
of the estimator should be computed over an ensemble of realizations of the
primary CMB which we denote by $\langle\rangle_{\Theta}$. In this ensemble,
the mean of the estimator is:

$\langle\hat{v}_{LM}^{\alpha}\rangle_{\Theta}=\sum_{L^{\prime}M^{\prime}}^{L^{\prime}_{\rm{max}}}\sum_{\beta=0}^{N-1}\left([\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{-1}\right)^{\alpha\beta}_{LM;L^{\prime}M^{\prime}}\langle\left[\boldsymbol{\tau}^{\dagger}(C^{\mathrm{pCMB}})^{-1}\Theta\right]_{L^{\prime}M^{\prime}}^{\beta}\rangle_{\Theta}.$
(22)

We find that

$\langle\hat{v}_{LM}^{\alpha}\rangle_{\Theta}=v_{LM}^{\alpha}+\beta_{LM}^{\alpha},$
(23)

where $\beta$ is a “coarse-graining” bias

$\beta_{LM}^{\alpha}\equiv\sum_{L^{\prime}M^{\prime}}^{L_{\rm{max}}}\sum_{\gamma=0}^{N-1}\sum_{c=N}^{\infty}\sum_{L^{\prime\prime}M^{\prime\prime}}^{L_{\rm{max}}}\left([\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{-1}\right)^{\alpha\gamma}_{LM;L^{\prime}M^{\prime}}[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{\gamma
c}_{L^{\prime}M^{\prime};L^{\prime\prime}M^{\prime\prime}}v_{L^{\prime\prime}M^{\prime\prime}}^{c}.$
(24)

Details of the calculation that leads to (23) are in Appendix D.1. We also
calculate the variance of the estimator, and find that

$\langle\hat{v}^{\alpha}_{LM}\hat{v}^{\beta}_{L^{\prime}M^{\prime}}\rangle_{\Theta}=(C^{vv})_{LM;L^{\prime}M^{\prime}}^{\alpha\beta}+N_{LM;L^{\prime}M^{\prime}}^{\alpha\beta},$
(25)

with

$(C^{vv})_{LM;L^{\prime}M^{\prime}}^{\alpha\beta}=(v_{LM}^{\alpha}+\beta_{LM}^{\alpha})(v_{L^{\prime}M^{\prime}}^{\beta}+\beta_{L^{\prime}M^{\prime}}^{\beta}),$
(26)

and

$N_{LM;L^{\prime}M^{\prime}}^{\alpha\beta}=\left([\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{-1}\right)^{\alpha\beta}_{LM;L^{\prime}M^{\prime}}.$
(27)

Unlike the QE Deutsch et al. (2018a); Cayuso et al. (2021), where one must
assume that the underlying fields are Gaussian, it is not necessary to make
any assumptions about the statistics of the optical depth or velocity fields.
Therefore, we do not expect additional contributions to the estimator variance
beyond what is presented here. This is in contrast to the QE, where extra
contributions to the estimator variance arise from the non-Gaussian nature of
the optical depth and velocities Giri and Smith (2020). We view this as a
significant advantage of the MaxL estimator, since these non-Gaussian
contributions are difficult to compute and remove.

Evaluating the harmonic space MaxL estimator in Eq. (21) requires costly, and
potentially numerically inaccurate, forward and reverse spherical harmonic
transforms. Motivated by this, in the next section we construct a map space
estimator. Ultimately, the implementation of the map space estimator is
simpler and more accurate, and will be the method used for the forecasts
presented below. However, there may be contexts in which the harmonic space
estimator is superior, as we comment on below.

###  IV.2 Map space estimator

We now derive an alternative coarse graining procedure and obtain a MaxL
estimator defined in terms of fields in map space. In map space, (as
demonstrated in Appendix A) the maximum likelihood condition can be written
explicitly as

$\tau^{\alpha}(\hat{n})\bar{\Theta}(\hat{n})=\sum_{\beta=0}^{N-1}\tau^{\alpha}(\hat{n})\bar{\tau}^{\beta}(\hat{n})v^{\beta}(\hat{n}),$
(28)

where the bar denotes inverse-variance filtering by
$C_{\ell}^{\mathrm{pCMB}}$:

$\bar{A}(\hat{n})\equiv\sum_{\ell m}\frac{A_{\ell
m}}{C_{\ell}^{\mathrm{pCMB}}}Y_{\ell m}(\hat{n}).$ (29)

To coarse-grain, we define a set of smoothing kernels $W_{I}(\hat{n})$ with
compact support in regions of the map labeled by $I=1,2,\ldots I_{\rm max}$
such that

$\int d\Omega\
W_{I}(\hat{n})W_{J}(\hat{n})=\frac{\delta_{IJ}}{\Delta\theta^{2}}$ (30)

where we have assumed each kernel is non-zero over the same solid angle
$\Delta\theta^{2}$ for all $I$. A simple choice for $W_{I}(\hat{n})$
satisfying this criterion, utilized below, is a simple average over equal-area
pixels. Note that this is only one choice of definition for $W_{I}(\hat{n})$
among many possibilities. Multiplying both sides of Eq. (28) by
$W_{I}(\hat{n})$ and integrating over the map, we obtain

$\int d\Omega\
W_{I}(\hat{n})\tau^{\alpha}(\hat{n})\bar{\Theta}(\hat{n})=\sum_{\beta=0}^{N-1}\int
d\Omega\
W_{I}(\hat{n})\tau^{\alpha}(\hat{n})\bar{\tau}^{\beta}(\hat{n})v^{\beta}(\hat{n}).$
(31)

Expanding the velocity in terms of an average velocity $v^{\beta}_{I}$ over
the region of the map labeled by $I$ and a residual $\Delta
v^{\beta}_{I}(\hat{n})$ in the same region, we have

$\int d\Omega\
W_{I}(\hat{n})\tau^{\alpha}(\hat{n})\bar{\Theta}(\hat{n})=\sum_{\beta=0}^{N-1}v^{\beta}_{I}\int
d\Omega\
W_{I}(\hat{n})\tau^{\alpha}(\hat{n})\bar{\tau}^{\beta}(\hat{n})+\sum_{\beta=0}^{N-1}\int
d\Omega\ W_{I}(\hat{n})\tau^{\alpha}(\hat{n})\bar{\tau}^{\beta}(\hat{n})\Delta
v^{\beta}_{I}(\hat{n}).$ (32)

Defining the operators

$\boldsymbol{\mathcal{W}}_{I;\alpha\beta}\equiv\int d\Omega\
W_{I}(\hat{n})\tau^{\alpha}(\hat{n})\bar{\tau}^{\beta}(\hat{n}),$ (33)

we have the following estimator for the average velocity in region $I$:

$\hat{v}_{I}^{\alpha}=\sum_{\beta=0}^{N-1}\boldsymbol{\mathcal{W}}_{I;\alpha\beta}{}^{-1}\int
d\Omega\ W_{I}(\hat{n})\tau^{\beta}(\hat{n})\bar{\Theta}(\hat{n}).$ (34)

In Appendix C, we demonstrate that the QE in the low signal-to-noise regime is
recovered upon replacing
$\boldsymbol{\mathcal{W}}_{I;\alpha\beta}\rightarrow\langle\boldsymbol{\mathcal{W}}_{I;\alpha\beta}\rangle_{\tau}$.

In our implementation of the map space estimator, we employ the Healpix equal-
area pixelization scheme Górski et al. (2005). If we define the input fields
at a native Healpix resolution with pixels indexed by $i=1,2,\ldots N_{\rm
pix}^{\rm in}$ and the regions $I$ are defined by pixels at a degraded Healpix
resolution $N_{\rm pix}^{\rm out}$, then we can write a discrete version of
the MaxL estimator:

$\hat{v}_{I}^{\alpha}=\sum_{\beta=0}^{N-1}\left[\sum_{i=1}^{N_{\rm pix}^{\rm
in}}W_{Ii}\
\tau^{\alpha}_{i}\bar{\tau}^{\beta}_{i}\right]^{-1}\left[\sum_{i=1}^{N_{\rm
pix}^{\rm in}}\ W_{Ii}\ \tau^{\beta}_{i}\bar{\Theta}_{i}\right].$ (35)

Note that $W_{Ii}$ as defined here is simply the Healpix ud_grade function –
an operator that averages pixels at a fine resolution labeled by $i$ that fit
within a coarse pixel labeled by $I$ (see Ref. Górski et al. (2005)).

Taking an ensemble average over realizations of the primary CMB, the estimator
mean is

$\langle\hat{v}_{I}^{\alpha}\rangle_{\Theta}=v_{I}^{\alpha}+\beta_{I}^{\alpha}$
(36)

where the additive bias $\beta_{I}^{\alpha}$ arises from coarse graining:

$\displaystyle\beta_{I}^{\alpha}$ $\displaystyle=$
$\displaystyle\sum_{\beta;\gamma=0}^{N-1}\boldsymbol{\mathcal{W}}_{I;\alpha\beta}{}^{-1}\int
d\Omega\ W_{I}(\hat{n})\tau^{\beta}(\hat{n})\bar{\tau}^{\gamma}(\hat{n})\Delta
v^{\gamma}_{I}(\hat{n})$ (37) $\displaystyle+$
$\displaystyle\sum_{\beta=0}^{N-1}\sum_{c=N}^{\infty}\boldsymbol{\mathcal{W}}_{I;\alpha\beta}{}^{-1}\int
d\Omega\
W_{I}(\hat{n})\tau^{\beta}(\hat{n})\bar{\tau}^{c}(\hat{n})v^{c}(\hat{n}).$

The first term is the bias from the pixel-averaging operation and the second
term is the bias from coarse graining along the line of sight. Moving to the
estimator variance:

$\displaystyle\langle\hat{v}_{I}^{\alpha}\hat{v}_{J}^{\beta}\rangle_{\Theta}$
$\displaystyle=$
$\displaystyle\left(v_{I}^{\alpha}+\beta_{I}^{\alpha}\right)\left(v_{J}^{\beta}+\beta_{J}^{\beta}\right)+\frac{1}{\Delta\theta^{2}}\boldsymbol{\mathcal{W}}_{I;\alpha\beta}{}^{-1}\delta_{IJ}$
(38)

Just as for the harmonic space estimator, the estimator variance includes an
additive bias due to coarse graining as well as a reconstruction noise. Note
that the reconstruction noise is local, and does not have any pixel-pixel
covariance (although there is bin-bin covariance). A derivation of the
estimator mean and variance can be found in Appendix D.2.

The map space estimator is local, only depending on the fields in the region
supported by $W_{I}(\hat{n})$. This allows us to construct $I_{\rm max}$
independent $N\times N$ operators $\boldsymbol{\mathcal{W}}_{I;\alpha\beta}$,
which are more easily (and accurately) inverted than the single $(N_{\rm
harm}^{2}N)\times(N_{\rm harm}^{2}N)$ operator needed for the harmonic space
MaxL estimator. The local nature of the estimator also makes it a convenient
choice for masked datasets. Indeed, one could generalize the discussion above
to arbitrary complex pixelizations, suited to a particular procedure for
masking, foreground removal, or anisotropic noise.

##  V Implementation and forecasts

In this section, we describe an implementation of the MaxL estimator and
compare its performance to the existing QE introduced in Ref. Deutsch et al.
(2018a). Unlike the QE, whose performance can be forecast from theoretical
spectra, the performance of the MaxL estimator depends on the realization. We
therefore must generate ensembles of simulations to compare the performance of
the MaxL estimator to the QE. In this paper, we perform reconstructions on
ensembles of properly correlated Gaussian simulations generated using the
algorithm of Cayuso et al. (2021); an analysis of N-body simulations is an
important direction for future work. Specifically, we first calculate the
theoretical power spectra of the lensed primary CMB $C_{\ell}^{pCMB}$; of the
velocity field $C_{\ell}^{vv}$; and of the electron density field
$C_{\ell}^{\tau\tau}$. We neglect the cross-power spectrum of electrons and
velocity $C_{\ell}^{\tau v}$ since the estimators utilize the optical depth
and velocities on very different scales. We use CAMB Lewis et al. (2000) to
calculate quantities such as the primary CMB and linear matter power spectrum,
with the non-linear regimes modelled with a halo model Cayuso et al. (2021)
(see e.g. Ref. Cooray and Sheth (2002) for a review of the halo model). The
electrons are modelled with the density profiles of Battaglia (2016). These
power spectra are then used to create appropriately correlated Gaussian
simulations of the binned velocity $v$ and $\tau$, as well as the primary CMB.
To model the non-Gaussian kSZ map, we simply multiply and sum the simulated
$v$ and $\tau$ maps: $\Theta^{\rm{kSZ}}=\sum v^{\alpha}\tau^{\alpha}$. We use
32 redshift bins between $z=0.2$ and $z=5$, and employ a binning scheme such
that each bin has equal comoving width. Further details of this implementation
are given in Sec. III of Cayuso et al. (2021). We consider the idealized
scenario where all non-blackbody foregrounds are absent and where we neglect
the non-Gaussianities induced by lensing of the primary CMB. A discussion of
these effects in the context of the QE can be found in Ref. Cayuso et al.
(2021).

In Fig. 1 we show the ensemble averaged power spectra of the (lensed) primary
CMB and of the kSZ effect. Note that we only include power from the 32
redshift bins between $0.2<z<5$ used in our analysis, while Ref. Cayuso et al.
(2021) demonstrated that approximately $512$ radial bins are necessary to
capture the majority of the kSZ power in this redshift range. We continue with
this model in order to be self-consistent, as using a larger number of bins
would be too computationally expensive. We do not expect this approximation to
affect our conclusions.

Figure 1: The primary CMB power spectrum and the kSZ power spectrum using 32
radial bins of equal comoving size between $0.2<z<5$. At large-$\ell$, the kSZ
power dominates over the CMB power. It is only once we access these angular
scales that we expect the MaxL estimator to show improvements over the QE.

Throughout this section, we consider the case where we have access to a direct
measurement of the optical depth. In practice, we will not have access to such
a map: one must use a tracer of the electrons such as a galaxy survey.
Connecting such a tracer to the optical depth necessitates further modelling,
in particular through the cross power spectrum $C_{\ell}^{\tau g}$. Incorrect
modelling of $C_{\ell}^{\tau g}$ leads to a bias on the reconstruction known
as the optical depth bias; we comment on this further in Section VI.
Additionally, we first consider the case where $C_{\ell}^{\tau\tau}$ is known
to arbitrarily high resolution, and don’t consider realistic effects such as
the shot noise of a galaxy survey or redshift error on a photometric survey.

We have validated two pipelines, based the map and harmonic space MaxL
estimators presented above. We find that the map space estimator has a number
of computational advantages. The harmonic space estimator involves several
forward and reverse spherical harmonic transforms, which in the Healpix
pixelization are not information-preserving. This loss of information leads to
numerical inaccuracies in the reconstruction which can be overcome by using
exact, but less efficient, pixelization schemes. In addition, the operator
inversion required in the velocity estimator is costly when incorporating the
full multipole and bin covariance; this places computational limitations on
the degree of coarse graining that can be considered. The map space estimator
overcomes both of these challenges, and so we focus on this implementation in
what follows.

In our first illustration of the MaxL estimator, we completely remove the
primary CMB and perform reconstruction on kSZ maps alone (we add a negligible
amount of white noise to ensure that we are not dividing by zero when we
filter). We use the equal-area Healpix pixelization Górski et al. (2005),
where the resolution is set by the parameter $N_{\rm side}$, related to the
number of pixels in the map by $N_{\rm pix}=12N_{\rm side}^{2}$. We apply both
the MaxL estimator and the QE to the kSZ-only and optical depth maps for a
single realization at a native resolution of $N_{\rm side}^{\rm in}=2048$.
Recalling the discussion in Sec. II, we expect the QE to be limited by
reconstruction noise, while the MaxL estimator is formally noiseless if the
maximum likelihood condition can be inverted without coarse graining. Some
degree of coarse graining is necessary for a well-defined solution to the
maximum likelihood condition, and we therefore expect the MaxL estimator to be
limited by the additive coarse graining bias as described in the previous
section.

The power spectrum of the true velocity and the power spectrum of the
reconstruction residuals (defined by $\hat{v}-v$) are shown in Fig. 2. For the
QE, we perform the reconstruction at a resolution of $N_{\rm side}^{\rm
out}=128$; for the MaxL estimator, we perform the reconstruction at a variety
of resolutions $N_{\rm side}^{\rm out}=16,64,128$. For low-resolution output
maps ($N_{\rm side}^{\rm out}=16$) the MaxL estimator in fact performs worse
than the QE due to a large coarse-graining bias $\beta$, but as we increase
$N_{\rm side}^{\rm out}$ the MaxL estimator becomes superior as expected. We
have explicitly checked that the residual agrees with $\beta$ as calculated
directly from (37).

Figure 2: We show the true velocity power spectrum as well as the power of the
residuals $v-\hat{v}^{\rm{QE}}$ and $v-\hat{v}^{\rm{MaxL}}$. On the left is a
low-$z$ bin and on the right is a high-$z$ bin. We show $\hat{v}^{\rm{MaxL}}$
for increasing resolution $N_{\rm side}^{\rm out}$. At low resolution, the
“coarse-graining” bias $\beta$ is significant and $\hat{v}^{\rm{MaxL}}$
performs worse than $\hat{v}^{\rm{QE}}$, but as we increase $N_{\rm side}^{\rm
out}$ this bias decreases and we get a better reconstruction than
$\hat{v}^{\rm{QE}}$.

When we include the primary CMB in the temperature maps, we do not expect the
improvement of the MaxL estimator over the QE to be as dramatic. In this case,
the MaxL estimator will have both reconstruction noise and an additive coarse-
graining bias (see Eq. (38)). We only expect the MaxL estimator to improve
upon the QE when the resolution of the input maps is high enough to probe the
regime where the amplitude of the kSZ anisotropies are comparable to those of
the primary CMB. To test this expected resolution-dependence, we produce input
maps at $N_{\rm side}^{\rm in}=1024,\ 2048$. In Fig. 3 we compare the velocity
reconstruction using the quadratic and MaxL estimators for $N_{\rm side}^{\rm
out}=16,64,128$ using maps with $N_{\rm side}^{\rm in}=2048$. In Fig. 4 we
plot the ratio of the total signal to noise for the MaxL estimator and the QE;
the signal to noise is defined as

$S/N_{L}=\sqrt{\mathrm{Tr}\left[C_{L}^{v^{\alpha}v^{\beta}}\left(N_{L}^{\hat{v}^{\alpha}\hat{v}^{\beta}}{}\right)^{-1}\right]}$
(39)

where $C_{L}^{v^{\alpha}v^{\beta}}$ is the (true) signal matrix and we use the
power of the residual as a proxy for the noise
$N_{L}^{\hat{v}^{\alpha}\hat{v}^{\beta}}$.

We see that the signal to noise for the MaxL estimator is comparable to that
of the QE for $N_{\rm side}^{\rm in}=1024$. Increasing to $N_{\rm side}^{\rm
in}=2048$, the MaxL signal to noise per mode is approximately 25% higher than
for the QE on small angular scales, and significantly higher on large angular
scales. This demonstrates the improvement in performance expected once angular
scales on which kSZ power is significant compared to the primary CMB are
included.

Figure 3: We show the true velocity power spectrum as well as the power of the
residuals $v-\hat{v}^{\rm QE}$ and $v-\hat{v}^{\rm MaxL}$. On the left is a
low-$z$ bin and on the right is a high-$z$ bin. We show $\hat{v}^{\rm MaxL}$
for an increasing output resolution. Here we see less of an improvement over
the QE than we did in Fig. 2 where we did not include the primary CMB. The
MaxL residual power is $~{}75\%$ of the QE residual power on small angular
scales, and significantly better on large angular scales.  Figure 4: The ratio
of the total signal to noise per mode of the MaxL and quadratic estimators. We
find that at $N_{\rm side}^{\rm in}$ of 1024, the estimators perform
similarly; however, as we go to $N_{\rm{side}}^{\rm{in}}$ of 2048 and higher,
the MaxL estimator has higher signal to noise than the QE. In both cases we
have chosen a different $N_{\rm side}^{\rm out}$ at which the MaxL performs
the best: if $N_{\rm side}^{\rm out}$ is too low, the coarse-graining bias is
too large, and if it is to high, the operator $\mathcal{W}$ becomes un-
invertible. These $N_{\rm side}^{\rm out}$ corresponded to $N_{\rm side}^{\rm
out}=\\{16,64\\}$ for $N_{\rm side}^{\rm out}=\\{1024,2048\\}$ respectively.

From these results we see that with perfect knowledge of the optical depth
field and including only the primary CMB and kSZ temperature anisotropies, the
MaxL estimator does yield a higher fidelity reconstruction than the QE for
sufficiently high resolution data. In the next section, we take one step
towards increasing realism and compare the performance of the MaxL and
quadratic estimators when the optical depth is inferred from a galaxy survey.

##  VI Estimating the optical depth

In the absence of a direct measure of the differential optical depth, we must
estimate it from a tracer in order to implement the maximum likelihood
estimators introduced above. A simple estimator is:

$\hat{\tau}_{\ell m}^{\alpha}=\sum_{\beta,\gamma=0}^{N-1}\left(C^{\tau
g}\right)_{\ell}^{\alpha\beta}\left(\left[C^{gg}\right]^{-1}\right)_{\ell}^{\beta\gamma}g_{\ell
m}^{\gamma},$ (40)

where $g_{\ell m}^{\gamma}$ is the redshift-binned galaxy density,
$\left(C^{gg}\right)_{\ell}^{\beta\gamma}$ is the galaxy power spectrum, and
$\left(C^{\tau g}\right)_{\ell}^{\alpha\beta}$ is the cross-power with the
optical depth. The ensemble mean of this estimator is zero while the variance
is:

$\displaystyle\langle\hat{\tau}_{\ell m}^{\alpha}(\hat{\tau}_{\ell
m}^{\beta})^{\dagger}\rangle$ $\displaystyle=$
$\displaystyle\sum_{\gamma,\delta=0}^{N-1}\left(C^{\tau
g}\right)_{\ell}^{\alpha\gamma}\left(\left[C^{gg}\right]^{-1}\right)_{\ell}^{\gamma\delta}\left(C^{g\tau}\right)_{\ell}^{\delta\beta}.$
(41)

For correlated Gaussian fields Eq. (40) will yield an optical depth that has
the expected statistics. However, it is likely that there are superior ways to
faithfully estimate the optical depth from the galaxy density. We leave this
question to future work. Note that the estimator variance reduces to
$\left(C^{\tau\tau}\right)_{\ell}^{\alpha\beta}$ in the limit where the actual
realization of the optical depth is related to the realization of the galaxy
density by $\tau_{\ell m}^{\alpha}=\sum b_{\ell}^{\alpha\gamma}g_{\ell
m}^{\gamma}$; that is, when the optical depth in a redshift bin can be
obtained by isotropically filtering a linear combination of the tracers in
different bins.

An implementation of the MaxL estimator in this scenario involves two steps:
first one estimates the optical depth field, then one uses in the MaxL
estimator to obtain the reconstructed velocity. To compare the MaxL and
quadratic estimators we must now include a galaxy field in our correlated
Gaussian simulations. Following in the implementation in Ref. Cayuso et al.
(2021), we generate galaxy auto- and cross-spectra for an LSST-like galaxy
survey. We refer the reader to Ref. Cayuso et al. (2021) for complete details,
but in brief: the galaxy spectra are generated using a halo model
incorporating a realistic halo occupation distribution (HOD), with number
density of galaxies in the survey assumed to be

$n(z)=\frac{n_{g}}{2z_{0}}\left(\frac{z}{z_{0}}\right)^{2}\exp\left(-\frac{z}{z_{0}}\right),$
(42)

with $z_{0}=0.3$ and $n_{g}=40\ {\rm arcmin}^{-2}$. The HOD and number density
of objects determine the galaxy bias and galaxy shot noise in
$\mathbf{C}_{\ell}^{gg}$; the halo model is used to determine a self-
consistent $\mathbf{C}_{\ell}^{g\tau}$. Using these spectra as inputs, we
generate properly correlated Gaussian maps of all fields at a Healpix
resolution of $N_{\rm side}^{\rm in}=2048$. Redshift errors associated with a
photometric galaxy survey such as LSST lead to bin-bin correlations in the
galaxy spectra on both large and small scales; we defer the inclusion of this
effect to future work.

Using the Gaussian simulations described above, we implement the QE on the
galaxy density and temperature (kSZ and primary CMB) maps using the theory
spectra for $\mathbf{C}_{\ell}^{gg}$ and $\mathbf{C}_{\ell}^{g\tau}$ as
inputs. We then use Eq. (40) to estimate the optical depth from the galaxy
density and apply the MaxL estimator to the estimated optical depth and
temperature map. The reconstruction residuals $\hat{v}-v$ for the MaxL and
quadratic estimators are shown in Fig. 5. At low-redshfit (left panel), a
high-fidelity reconstruction is possible, albeit at a lower signal-to-noise
than what is possible using the optical depth itself as a tracer. At high
redshift all galaxy density measurements are shot-noise dominated, as the
galaxy density is very low. As such, it is not possible to reconstruct the
high-$z$ velocity; this is evident on the right hand side of Fig. 5. At both
low- and high-redshift, the MaxL estimator performs at least as well as the
QE, in this regime ($N^{\rm{in}}_{\rm{side}}=2048$). Given that an LSST-like
experiment is shot-noise dominated on small angular scales, the additional
information leveraged by the MaxL estimator from the optical depth on small
scales is significantly degraded. Because of this, it may be difficult to
identify an implementation of the MaxL estimator that would outperform the QE
for near-term experiments.

Figure 5: In this Figure we show the true velocity power spectrum $C_{L}^{vv}$
and the power of the residuals $v-\hat{v}^{QE}$ and $v-\hat{v}^{MaxL}$, using
galaxies to estimate the optical depth field. We do this for the MaxL by using
the procedure defined in (40). For the QE we use the simulated galaxy maps
directly and modify the theoretical weights in the estimator to ensure an
unbiased output.

When a tracer is used in the MaxL estimator, we must confront the inevitable
modelling uncertainty involved in estimating the optical depth field. In the
estimator Eq. (40), this modelling uncertainty arises through the theory
spectra $\left(C^{\tau g}\right)_{\ell}^{\alpha\beta}$ and
$\left(C^{gg}\right)_{\ell}^{\alpha\beta}$ that must be provided as inputs. As
we demonstrate explicitly in Appendix E, the consequence is an overall
multiplicative bias on the estimator mean, referred to in previous literature
(e.g. Battaglia (2016); Smith et al. (2018); Madhavacheril et al. (2019);
Hotinli et al. (2021)) as the optical depth bias. Neglecting the coarse-
graining bias, for the map space estimator we have:

$\langle\hat{v}_{I}^{\alpha}\rangle=\sum_{\gamma=0}^{N-1}R^{\alpha\gamma}_{I}v_{I}^{\gamma},$
(43)

where the optical depth bias matrix is

$R^{\alpha\gamma}_{I}\equiv\sum_{\beta=0}^{N-1}\left[\int d\Omega\
W_{I}(\hat{n})\hat{\tau}^{\beta}(\hat{n})\hat{\bar{\tau}}^{\alpha}(\hat{n})\right]^{-1}\int
d\Omega\
W_{I}(\hat{n})\hat{\tau}^{\beta}(\hat{n})\bar{\tau}^{\gamma}(\hat{n}).$ (44)

A hat on $\tau$ implies that it has been reconstructed from a tracer, while
the absence of a hat indicates the actual optical depth. In complete
generality, the optical depth bias matrix contains $N^{2}\times N_{\rm pix}$
free parameters. These free parameters are completely degenerate with the
$N\times N_{\rm pix}$ independent degrees of freedom in the coarse grained
velocity. Fortunately, as argued in Appendix E, we expect that the optical
depth bias does not depend on the pixel index $I$. Roughly, this follows from
the assumption that the optical depth field is statistically isotropic.
Furthermore, we expect that the dominant contribution to the off-diagonal
$\alpha\neq\gamma$ components of the optical depth matrix are expected to
arise from redshift errors in the tracer used to reconstruct the optical
depth. This implies that there are only $N$ free parameters in the optical
depth bias matrix, so long as the redshift errors of the tracer can be well-
characterized. This is consistent with the QE, where there are $N$ free
parameters expected in the optical depth bias Cayuso et al. (2021). A complete
investigation of the optical depth bias will be the focus of future work.

##  VII Conclusions

In this paper, we have introduced a new class of estimators for the radial
velocity field (more generally, the remote dipole field) from the CMB
temperature anisotropies induced by the kinetic Sunyaev Zel’dovich effect and
a tracer of the optical depth field. These MaxL estimators are based on a
coarse-grained maximum likelihood condition, and we have introduced two
coarse-graining schemes in harmonic and map space that can be used to solve
for the maximum likelihood radial velocity field. We derived the expected
properties of the MaxL estimator, and compared the performance of the MaxL
estimator to the existing quadratic estimator using the correlated Gaussian
simulation framework of Ref. Cayuso et al. (2021).

The MaxL estimator introduced here has several advantages over the existing QE
for velocity reconstruction. As we have demonstrated, the MaxL estimator can
in principle yield a higher fidelity reconstruction than the QE in the high-
resolution, low-noise frontier of future CMB experiments. The MaxL estimator
is not derived based on the assumption that the optical depth field is
Gaussian, and it therefore does not receive contributions from the non-
Gaussian bias terms present in the QE Giri and Smith (2020). The map space
version of the MaxL estimator is numerically stable and simple to implement.
It is absolutely local, making it a useful and versatile tool to apply to
datasets containing masks and spatially varying foregrounds. The modular
nature of the MaxL estimator makes it amenable to improvements in techniques
used to estimate the optical depth field or improved coarse-graining schemes.
Furthermore, applying both the quadratic and MaxL estimators to future
datasets will be a useful diagnostic of various biases and systematics in each
analysis.

There are a variety of directions to take in future work. An important first
step is a comparison of the quadratic and MaxL estimators using N-body data,
where some of the advantages of the MaxL estimator described above can be
confirmed in a more realistic setting. Another line of inquiry is to determine
better methods for estimating the optical depth field from a galaxy redshift
survey, or other tracers of large scale structure such as the cosmic infrared
background McCarthy and Johnson (2020). For example, machine learning
techniques trained on simulations could capture non-Gaussian information to
better model the tracer-optical depth correlations; such a framework could
incorporate multiple tracers into an overall estimate of the optical depth
field over a wide range of redshifts. It will also be important to develop an
analysis framework incorporating a mask and spatially varying foreground
mitigation techniques as exists for the QE Cayuso et al. (2021). Finally, one
can construct an analogous set of MaxL estimators based on other CMB
secondaries; for example, transverse velocity reconstruction from the moving
lens effect Hotinli et al. (2019) or remote quadrupole reconstruction from the
polarized Sunyaev Zel’dovich effect Deutsch et al. (2018b, a). Improved
techniques to produce unbiased high-fidelity reconstruction of cosmological
fields will yield significant scientific dividends for future experiments,
maximizing the ability of these datasets to constrain fundamental properties
of our Universe.

###### Acknowledgements.

We thank R. Bloch, C. Hill, S. Hotinli, M. Madhavacheril, J. Mertens, M.
Münchmeyer, and B. Sherwin for useful conversations and input at various
stages of this project. MCJ is supported by the National Science and
Engineering Research Council through a Discovery grant. 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.
Some of the results in this paper have been derived using the HEALPix package
Górski et al. (2005).

## References

  * Aguirre et al. (2019) J. Aguirre et al. (Simons Observatory), JCAP 1902, 056 (2019), eprint 1808.07445.
  * Abazajian et al. (2016) K. N. Abazajian et al. (CMB-S4) (2016), eprint 1610.02743.
  * Sehgal et al. (2019) N. Sehgal et al. (2019), eprint 1906.10134.
  * Sunyaev and Zeldovich (1980) R. A. Sunyaev and I. B. Zeldovich, MNRAS 190, 413 (1980).
  * Deutsch et al. (2018a) A.-S. Deutsch, E. Dimastrogiovanni, M. C. Johnson, M. Münchmeyer, and A. Terrana, Phys. Rev. D98, 123501 (2018a), eprint 1707.08129.
  * Smith et al. (2018) K. M. Smith, M. S. Madhavacheril, M. Münchmeyer, S. Ferraro, U. Giri, and M. C. Johnson (2018), eprint 1810.13423.
  * Cayuso et al. (2021) J. Cayuso, R. Bloch, S. C. Hotinli, M. C. Johnson, and F. McCarthy (2021), eprint 2111.11526.
  * Terrana et al. (2017) A. Terrana, M.-J. Harris, and M. C. Johnson, Journal of Cosmology and Astroparticle Physics 2017, 040 (2017), eprint 1610.06919.
  * Cayuso et al. (2018) J. I. Cayuso, M. C. Johnson, and J. B. Mertens, Phys. Rev. D 98, 063502 (2018), URL https://link.aps.org/doi/10.1103/PhysRevD.98.063502.
  * Giri and Smith (2020) U. Giri and K. M. Smith (2020), eprint 2010.07193.
  * Madhavacheril et al. (2019) M. S. Madhavacheril, N. Battaglia, K. M. Smith, and J. L. Sievers (2019), eprint 1901.02418.
  * Battaglia (2016) N. Battaglia, JCAP 1608, 058 (2016), eprint 1607.02442.
  * Hotinli et al. (2021) S. C. Hotinli, K. M. Smith, M. S. Madhavacheril, and M. Kamionkowski (2021), eprint 2108.02207.
  * Thiele et al. (2020) L. Thiele, F. Villaescusa-Navarro, D. N. Spergel, D. Nelson, and A. Pillepich, Astrophys. J. 902, 129 (2020), eprint 2007.07267.
  * Hadzhiyska et al. (2019) B. Hadzhiyska, B. D. Sherwin, M. Madhavacheril, and S. Ferraro, Phys. Rev. D 100, 023547 (2019), eprint 1905.04217.
  * Hotinli and Johnson (2020) S. C. Hotinli and M. C. Johnson (2020), eprint 2012.09851.
  * Górski et al. (2005) K. M. Górski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and M. Bartelmann, ApJ 622, 759 (2005), eprint astro-ph/0409513.
  * Akrami et al. (2020) Y. Akrami et al. (Planck), Astron. Astrophys. 644, A100 (2020), eprint 2003.12646.
  * Lewis et al. (2000) A. Lewis, A. Challinor, and A. Lasenby, ApJ 538, 473 (2000), eprint astro-ph/9911177.
  * Cooray and Sheth (2002) A. Cooray and R. K. Sheth, Phys. Rept. 372, 1 (2002), eprint astro-ph/0206508.
  * Battaglia (2016) N. Battaglia, J. Cosmology Astropart. Phys 2016, 058 (2016), eprint 1607.02442.
  * McCarthy and Johnson (2020) F. McCarthy and M. C. Johnson, Phys. Rev. D 102, 043520 (2020), eprint 1907.06678.
  * 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, Phys. Rev. Lett. 123, 061301 (2019), eprint 1812.03167.
  * Deutsch et al. (2018b) A.-S. Deutsch, M. C. Johnson, M. Münchmeyer, and A. Terrana, Journal of Cosmology and Astroparticle Physics 2018, 034 (2018b), eprint 1705.08907v3.

##  Appendix A Harmonic and map space representations

In this paper, we work both in harmonic space and in map space (also referred
to as “pixel” or “real” space). In each case, we are considering fields in a
different basis: in map space, we describe fields by associating a value to
each pixel on a sphere, while in harmonic space we describe fields by
expanding in the basis of spherical harmonics. In this Appendix, we briefly
recall how to transform from one basis to another, considering both the fields
themselves, and the operators on these fields.

###  A.1 Fields

A 2-dimensional field $A(\hat{n})$ defined on a sphere can be represented in
the map or harmonic basis:

$a(\hat{n})=\sum_{i}a_{i}P_{i}(\hat{n})=\sum_{\ell m}a_{\ell m}Y_{\ell
m}(\hat{n})$ (45)

where $a_{i}$ is the value of the field at each pixel $P_{i}(\hat{n})$ and
$a_{\ell m}$ are the spherical harmonic coefficients of the field decomposed
onto the spherical harmonic basis $Y_{\ell m}(\hat{n})$. In the map basis,
each pixel can be denoted $i$ and the value $a_{i}$ is given by

$a_{i}=\int d\Omega\ a(\hat{n})P_{i}(\hat{n})$ (46)

where $P_{i}$ is an element of the discrete pixelization basis. In the
harmonic basis, fields can be represented as a sum over the spherical harmonic
basis $Y_{\ell m}$ with the $a_{\ell m}$’s given by

$a_{\ell m}=\int d\Omega\ a(\hat{n})Y_{\ell m}^{*}(\hat{n}).$ (47)

In particular, we can go between map-space and harmonic-space according to

$\displaystyle a_{\ell m}=$ $\displaystyle\int
d\Omega\left(\sum_{i}a_{i}P_{i}(\hat{n})\right)Y_{\ell
m}^{*}(\hat{n})=\sum_{i}a_{i}Y_{\ell m}^{*}(\hat{n}_{i});$ (48) $\displaystyle
a_{i}=$ $\displaystyle\int d\Omega\sum_{\ell m}a_{\ell m}Y_{\ell
m}(\hat{n})P_{i}(\hat{n})=\sum_{\ell m}a_{\ell m}Y_{\ell m}(\hat{n}_{i}),$
(49)

where $Y_{\ell m}(\hat{n}_{i})\equiv\int d\Omega P_{i}(\hat{n})Y_{\ell
m}(\hat{n})$, e.g. the spherical harmonic weighted over the pixel centered on
$\hat{n}_{i}$.

###  A.2 Operators

We can decompose operators $\mathbf{O}$ on either basis. Thus the expression

$\boldsymbol{O}b=c$ (50)

can be written as

$\sum_{j}\boldsymbol{O}_{ij}b_{j}=c_{i},$ (51)

or

$\sum_{LM}\boldsymbol{O}_{\ell m;LM}b_{LM}=c_{\ell m}.$ (52)

We can also consider a mixed basis

$\sum_{LM}\boldsymbol{O}_{i;LM}b_{LM}=c_{i}.$ (53)

We also have the relation

$\displaystyle\boldsymbol{O}_{\ell m;LM}b_{LM}=$
$\displaystyle\sum_{i}\boldsymbol{O}_{ij}b_{j}Y_{\ell m}^{*}(\hat{n}_{i})$
(54) $\displaystyle=$
$\displaystyle\sum_{i}\boldsymbol{O}_{ij}\sum_{\ell^{\prime}m^{\prime}}b_{\ell^{\prime}m^{\prime}}Y_{\ell^{\prime}m^{\prime}}(\hat{n}_{j})Y_{\ell
m}^{*}(\hat{n}_{i}).$ (55)

###  A.3 The maximum likelihood condition in map and harmonic space

The maximum likelihood condition is

$\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\Theta\right]=\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}\right]\cdot
v.$ (56)

To write this in harmonic or map space, we need to form both the operator
$\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}$
and the vector
$\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\Theta$ in
map and harmonic space. Let us first consider the operator
$\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}$.

First, it will be useful to write the elements of $\boldsymbol{\tau}$ and
$\boldsymbol{C}^{\mathrm{\mathrm{pCMB}}}$. In map space, $\boldsymbol{\tau}$
is diagonal:

$\boldsymbol{\tau}^{\alpha}_{ij}=\tau^{\alpha}_{i}\delta_{ij}$ (57)

where $\tau^{\alpha}_{i}$ are the components of the $\tau^{\alpha}$ field in
the position basis and $\delta_{ij}$ is the Kronecker delta, the identity on
map space. Using the basis-transformation Eq. (55) for a diagonal operator on
map space gives

$\boldsymbol{\tau}^{\alpha}_{\ell
m;LM}=\sum_{\ell^{\prime}m^{\prime}}\tau^{\alpha}_{\ell^{\prime}m^{\prime}}\int
d\Omega\ Y_{\ell^{\prime}m^{\prime}}(\hat{n})Y_{LM}(\hat{n})Y_{\ell
m}^{*}(\hat{n}).$ (58)

We can also write a mixed-basis expression

$\boldsymbol{\tau}^{\alpha}_{i;LM}=\tau^{\alpha}_{i}Y_{LM}(\hat{n}_{i}).$ (59)

It will be most helpful to consider $\boldsymbol{C}^{\mathrm{\mathrm{pCMB}}}$
in harmonic space, where it is diagonal:

$\boldsymbol{C}^{\mathrm{\mathrm{pCMB}}}_{\ell
m;\ell^{\prime}m^{\prime}}=C_{\ell}^{\mathrm{\mathrm{pCMB}}}\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}},$
(60)

where $C_{\ell}^{\mathrm{\mathrm{pCMB}}}$ is the power spectrum of the primary
CMB. Then we can evaluate

$\displaystyle\left[(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}\right]^{\alpha}_{\ell^{\prime\prime}m^{\prime\prime}LM}$
$\displaystyle=$ $\displaystyle\sum_{\ell
m}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}_{\ell^{\prime\prime}m^{\prime\prime}\ell
m}\boldsymbol{\tau}^{\alpha}_{\ell m;LM}$ $\displaystyle=$
$\displaystyle\sum_{\ell m;\ \ell^{\prime}m^{\prime}}\left[\int d\Omega\
Y_{\ell^{\prime}m^{\prime}}Y_{LM}Y^{*}_{\ell
m}\right]\frac{\tau^{\alpha}_{\ell^{\prime}m^{\prime}}}{C_{\ell}^{\mathrm{pCMB}}}\delta_{\ell\ell^{\prime\prime}}\delta_{mm^{\prime\prime}}$
$\displaystyle=$
$\displaystyle\sum_{\ell^{\prime}m^{\prime}}(-1)^{m^{\prime\prime}}\sqrt{\frac{(2\ell^{\prime}+1)(2L+1)(2\ell^{\prime\prime}+1)}{4\pi}}\begin{pmatrix}L&\ell^{\prime}&\ell^{\prime\prime}\\\
0&0&0\end{pmatrix}\begin{pmatrix}L&\ell^{\prime}&\ell^{\prime\prime}\\\
M&m^{\prime}&-m^{\prime\prime}\end{pmatrix}\frac{\tau^{\alpha}_{\ell^{\prime}m^{\prime}}}{C_{\ell^{\prime\prime}}^{\mathrm{pCMB}}}$

where $\begin{pmatrix}L&\ell^{\prime}&\ell^{\prime\prime}\\\
M&m^{\prime}&-m^{\prime\prime}\end{pmatrix}$ is the Wigner $3J$-symbol.
Without making further assumptions, this is as far as we can go. However, we
will mainly be interested in the regime where $L\ll\ell^{\prime},\
\ell^{\prime\prime}$. The triangle inequality restricts $\ell$ and
$\ell^{\prime}$ to the range $|\ell^{\prime\prime}-\ell^{\prime}|\leq
L\leq\ell^{\prime\prime}+\ell^{\prime}$. So, as long as
$C_{\ell^{\prime\prime}}^{\mathrm{pCMB}}$ is not a strongly varying function,
we can approximate:

$C_{\ell^{\prime\prime}}^{\mathrm{pCMB}}\simeq C_{\ell^{\prime\prime}\pm
L}^{\mathrm{pCMB}}\simeq C_{\ell^{\prime}}^{\mathrm{pCMB}}.$ (62)

Under this approximation, we have:

$\displaystyle\sum_{\ell
m}\left[(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\right]_{\ell^{\prime\prime}m^{\prime\prime};\ell
m}\boldsymbol{\tau}^{\alpha}_{\ell m;LM}$ $\displaystyle=$ $\displaystyle\int
d\Omega\
\left[\sum_{\ell^{\prime}m^{\prime}}\frac{\tau^{\alpha}_{\ell^{\prime}m^{\prime}}}{C_{\ell^{\prime}}^{\mathrm{pCMB}}}Y_{\ell^{\prime}m^{\prime}}\right]Y_{LM}Y^{*}_{\ell^{\prime\prime}m^{\prime\prime}},$
(63)

which has a nice mixed pixel-space harmonic-space representation:

$\displaystyle\left[(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\right]_{ij}\boldsymbol{\tau}^{\alpha}_{j;LM}$
$\displaystyle=$ $\displaystyle\bar{\tau}^{\alpha}_{i}Y_{LM}(\hat{n}_{i})$
(64)

where $\bar{\tau}^{\alpha}_{i}$ is the inverse-variance filtered field

$\bar{\tau}^{\alpha}_{i}\equiv\sum_{\ell^{\prime}m^{\prime}}\frac{\tau^{\alpha}_{\ell^{\prime}m^{\prime}}}{C_{\ell^{\prime}}^{\mathrm{pCMB}}}Y_{\ell^{\prime}m^{\prime}}(\hat{n}_{i}).$
(65)

With these results, we can now form the operator
$\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}\right]$.
It is most convenient to do this using the mixed representation:

$\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}\right]^{\alpha\beta}_{L^{\prime}M^{\prime};LM}=\int
d\Omega\
\tau^{\alpha}(\hat{n})\bar{\tau}^{\beta}(\hat{n})Y_{L^{\prime}M^{\prime}}(\hat{n})^{*}Y_{LM}(\hat{n}).$
(66)

Using an analogous set of arguments, we have:

$\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\Theta\right]_{L^{\prime}M^{\prime}}^{\beta}=\int
d\Omega\
\tau^{\beta}(\hat{n})\bar{\Theta}(\hat{n})Y_{L^{\prime}M^{\prime}}^{*}(\hat{n})$
(67)

where, again, the bar denotes inverse-variance filtering:

$\bar{\Theta}_{i}\equiv\sum_{\ell^{\prime}m^{\prime}}\frac{\Theta_{\ell^{\prime}m^{\prime}}}{C_{\ell^{\prime}}^{\mathrm{pCMB}}}Y_{\ell^{\prime}m^{\prime}}(\hat{n}_{i}).$
(68)

The maximum likelihood condition is

$\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\Theta\right]_{L^{\prime}M^{\prime}}^{\alpha}=\sum_{LM;\
\beta}\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}\right]^{\alpha\beta}_{L^{\prime}M^{\prime};LM}v^{\beta}_{LM}.$
(69)

We can also rotate into map space by first writing

$\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}\right]^{\alpha\beta}_{i;LM}=\tau^{\alpha}_{i}\bar{\tau}^{\beta}_{i}\
Y_{LM}(\hat{n}_{i})$ (70)

and

$\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\Theta\right]_{i}^{\alpha}=\tau^{\alpha}_{i}\bar{\Theta}_{i}.$
(71)

The maximum likelihood condition is then:

$\tau^{\alpha}_{i}\bar{\Theta}_{i}=\sum_{LM;\
\beta}\tau^{\alpha}_{i}\bar{\tau}^{\beta}_{i}v^{\beta}_{LM}\
Y_{LM}(\hat{n}_{i})=\sum_{\beta=0}^{N-1}\tau^{\alpha}_{i}\bar{\tau}^{\beta}_{i}v^{\beta}_{i}.$
(72)

##  Appendix B Coarse-graining

To determine if there is a solution to the maximum likelihood condition Eq.
(56), we must determine if the operator
$\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}$
is invertible for a particular coarse graining procedure. Note that this
operator depends on the actual realization of the optical depth anisotropies,
so there is some danger that an inverse exists for some realizations but not
others. Nevertheless, we can get some idea of what to expect by examining the
statistical properties of this operator over an ensemble of realizations of
the optical depth anisotropies. For example, if we focus on the harmonic
basis, the mean of the operator is:

$\displaystyle\langle\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}\right]^{\alpha\beta}_{L^{\prime}M^{\prime};LM}\rangle_{\tau}$
$\displaystyle=$
$\displaystyle\sum_{\ell,\ell^{\prime}}\frac{(2\ell+1)(2\ell^{\prime}+1)}{4\pi}\left(\begin{array}[]{ccc}\ell&\ell^{\prime}&L\\\
0&0&0\end{array}\right)^{2}\frac{C_{\ell^{\prime}}^{\tau^{\alpha}\tau^{\beta}}}{C_{\ell}^{\mathrm{pCMB}}}\delta_{LL^{\prime}}\delta_{MM^{\prime}}$
(75)

where the $\langle\ldots\rangle_{\tau}$ denotes an ensemble average over the
optical depth anisotropies. On average, we see that the operator is diagonal
in harmonic space. If the sum is dominated by large-$\ell$ (which we expect it
to be) then it will also be nearly diagonal in bin-space. This is because the
bin-bin correlations in $C_{\ell}^{\tau^{a}\tau^{b}}$ become vanishingly small
at high-$\ell$ as the scales probed by the angular projection become much
smaller than the extent of the bins. Finding the variance of the operator:

$\displaystyle\langle\left(\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}\right]^{\alpha\beta}_{L^{\prime}M^{\prime};LM}\right)^{2}\rangle$
$\displaystyle\simeq$
$\displaystyle\langle\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}\right]^{\alpha\beta}_{L^{\prime}M^{\prime};LM}\rangle^{2}$
(80) $\displaystyle+$
$\displaystyle\left(\sum_{\ell}\frac{2\ell+1}{4\pi}\left[\frac{C_{\ell}^{\tau^{\alpha}\tau^{\alpha}}C_{\ell}^{\tau^{\beta}\tau^{\beta}}}{(C_{\ell}^{\mathrm{pCMB}})^{2}}+\frac{C_{\ell}^{\tau^{\alpha}\tau^{\beta}}C_{\ell}^{\tau^{\alpha}\tau^{\beta}}}{(C_{\ell}^{\mathrm{pCMB}})^{2}}\right]\right)\sum_{\ell^{\prime\prime}}(2\ell^{\prime\prime}+1)(2L+1)(2L^{\prime}+1)$
$\displaystyle\times$
$\displaystyle\left(\begin{array}[]{ccc}L&L^{\prime}&\ell^{\prime\prime}\\\
0&0&0\end{array}\right)^{2}\left(\begin{array}[]{ccc}L&L^{\prime}&\ell^{\prime\prime}\\\
M&M^{\prime}&M-M^{\prime}\end{array}\right)^{2}$

where we have assumed that the optical depth anisotropies are Gaussian and
used the limit where $L,L^{\prime}$ are small. Notably, the variance is not
diagonal in the harmonic coefficients, and becomes increasingly non-diagonal
at higher $L$ and $L^{\prime}$. It is also non-diagonal in bin space.
Therefore, in a particular realization we can expect that the operator
$\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}$
becomes highly non-diagonal on small angular scales. This implies that the
condition number could be quite large at full resolution, and there may be no
well-defined inverse. Our strategy is therefore to coarse-grain the maximum
likelihood condition by truncating at some empirically determined $L_{\rm
max}$ and $N$ before solving for the velocity field.

##  Appendix C Connecting to the quadratic estimator

###  C.1 Harmonic space estimator

In the low signal-to-noise regime, we can obtain the quadratic estimator from
the maximum likelihood condition Eq. (21) by taking an ensemble average over
the optical depth $\langle\rangle_{\tau}$:

$\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\Theta\right]_{L^{\prime}M^{\prime}}^{\alpha}=\sum_{LM;\
\beta}\langle\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}\right]^{\alpha\beta}_{L^{\prime}M^{\prime};LM}\rangle_{\tau}\
v^{\beta}_{LM}.$ (81)

The ensemble average is given by Eq. (75), which using the notation of Ref.
Cayuso et al. (2021) can be recognized as:

$\langle\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}\right]^{\alpha\beta}_{L^{\prime}M^{\prime};LM}\rangle_{\tau}=\frac{1}{2L+1}\sum_{\ell,\ell^{\prime}}G^{v^{\beta}\tau^{\beta}}_{\ell\ell^{\prime}L}f^{v^{\beta}\tau^{\alpha}}_{\ell
L\ell^{\prime}}\delta_{LL^{\prime}}\delta_{MM^{\prime}}$ (82)

with

$G^{v^{\beta}\tau^{\beta}}_{\ell\ell^{\prime}L}=\frac{f_{\ell
L\ell^{\prime}}^{v^{\beta}\tau^{\beta}}}{C_{\ell}^{\mathrm{pCMB}}C_{\ell^{\prime}}^{\tau^{\beta}\tau^{\beta}}};\
\ \ f_{\ell
L\ell^{\prime}}^{v^{\beta}\tau^{\alpha}}=\sqrt{\frac{(2\ell+1)(2\ell^{\prime}+1)(2L+1)}{4\pi}}\left(\begin{array}[]{ccc}\ell&\ell^{\prime}&L\\\
0&0&0\end{array}\right)C_{\ell^{\prime}}^{\tau^{\beta}\tau^{\alpha}}.$ (83)

The left hand side of the maximum likelihood condition can be expanded to:

$\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\Theta\right]_{LM}^{\alpha}=\sum_{\ell
m;\ell^{\prime}m^{\prime}}(-1)^{M}\left(\begin{array}[]{ccc}\ell&\ell^{\prime}&L\\\
m&m^{\prime}&-M\end{array}\right)G^{v^{\alpha}\tau^{\alpha}}_{\ell\ell^{\prime}L}\Theta_{\ell
m}\tau^{\alpha}_{\ell^{\prime}m^{\prime}}.$ (84)

For a completely diagonal
$C_{\ell^{\prime}}^{\tau^{\beta}\tau^{\alpha}}=C_{\ell^{\prime}}^{\tau^{\alpha}\tau^{\alpha}}\delta_{\alpha\beta}$,
we solve the ensemble-averaged maximum likelihood condition to arrive at the
traditional form for the quadratic estimator appearing in Refs. Deutsch et al.
(2018a); Cayuso et al. (2021):

$v^{\alpha}_{LM}=A_{L}^{v^{\alpha}}\sum_{\ell
m;\ell^{\prime}m^{\prime}}(-1)^{M}\left(\begin{array}[]{ccc}\ell&\ell^{\prime}&L\\\
m&m^{\prime}&-M\end{array}\right)G^{v^{\alpha}\tau^{\alpha}}_{\ell\ell^{\prime}L}\Theta_{\ell
m}\tau^{\alpha}_{\ell^{\prime}m^{\prime}}$ (85)

where

$A_{L}^{v^{\alpha}}\equiv\left[\frac{1}{2L+1}\sum_{\ell,\ell^{\prime}}G^{v^{\alpha}\tau^{\alpha}}_{\ell\ell^{\prime}L}f^{v^{\alpha}\tau^{\alpha}}_{\ell
L\ell^{\prime}}\right]^{-1}.$ (86)

There are, however, a few crucial differences with the quadratic estimator
presented in Refs. Deutsch et al. (2018a); Cayuso et al. (2021). In the
traditional quadratic estimator, the weights needed for an unbiased minimum
variance reconstruction require that everywhere above we replace
$C_{\ell}^{\mathrm{pCMB}}\rightarrow
C_{\ell}^{\mathrm{pCMB}}+C_{\ell}^{\rm{kSZ}}$. In the limit where
$C_{\ell}^{\mathrm{pCMB}}\ll C_{\ell}^{\rm{kSZ}}$, as occurs for a
sufficiently high-resolution low-noise CMB experiment, this is a non-
negligible difference – as expected, since in this case the temperature
anisotropies on small scales are very non-Gaussian, and far from the limit of
Gaussian fields used to derive the quadratic estimator. Therefore, we only
recover the traditional quadratic estimator in the limit where
$C_{\ell}^{\mathrm{pCMB}}\gg C_{\ell}^{\rm{kSZ}}$. This observation supports
the notion that the maximum likelihood estimator will be superior in the high
signal-to-noise regime. Another difference arises in the case where
$C_{\ell}^{\tau^{\beta}\tau^{\alpha}}$ is not diagonal. In this case, to find
a solution to the ensemble averaged maximum likelihood condition it is
necessary to find the inverse of $C_{\ell}^{\tau^{\beta}\tau^{\alpha}}$ (if it
exists). The resulting solution is different than the approach presented in
Ref. Cayuso et al. (2021), which first constructs a set of biased estimators
which are then de-biased by forming the appropriate linear combinations. The
relevant differences would be interesting to explore further in future work.

###  C.2 Map space estimator

Just as for the harmonic space estimator, we can derive an ensemble average
map space estimator. Starting from Eq. (32):

$\displaystyle\int d\Omega\
W_{I}(\hat{n})\tau^{\alpha}(\hat{n})\bar{\Theta}(\hat{n})$ $\displaystyle=$
$\displaystyle\sum_{\beta=0}^{N-1}v^{\beta}_{I}\ \int d\Omega\
W_{I}(\hat{n})\langle\tau^{\alpha}(\hat{n})\bar{\tau}^{\beta}(\hat{n})\rangle_{\tau}$
(87) $\displaystyle+$ $\displaystyle\sum_{\beta=0}^{N-1}\int d\Omega\
W_{I}(\hat{n})\langle\tau^{\alpha}(\hat{n})\bar{\tau}^{\beta}(\hat{n})\rangle_{\tau}\Delta
v^{\beta}_{I}(\hat{n})$ (88) $\displaystyle=$
$\displaystyle\sum_{\beta=0}^{N-1}v^{\beta}_{I}\
\langle\tau^{\alpha}(0)\bar{\tau}^{\beta}(0)\rangle_{\tau}\ \int d\Omega\
\left[W_{I}(\hat{n})+\int d\Omega\ W_{I}(\hat{n})\Delta
v^{\beta}_{I}(\hat{n})\right]$ (89) $\displaystyle=$
$\displaystyle\sum_{\beta=0}^{N-1}v^{\beta}_{I}\
\langle\tau^{\alpha}(0)\bar{\tau}^{\beta}(0)\rangle_{\tau}$ (90)

where we’ve assumed that $\tau^{\alpha}$ is statistically isotropic, $\int
d\Omega\ W_{I}(\hat{n})=1$, and $\int d\Omega\ W_{I}(\hat{n})\Delta
v^{\beta}_{I}(\hat{n})=0$ by definition. This suggests the following map space
quadratic estimator:

$\hat{v}^{\alpha}_{I}=\sum_{\beta=0}^{N-1}\left[\langle\tau^{\beta}(0)\bar{\tau}^{\alpha}(0)\rangle_{\tau}\right]^{-1}\int
d\Omega\ W_{I}(\hat{n})\tau^{\alpha}(\hat{n})\bar{\Theta}(\hat{n}).$ (91)

For a diagonal variance
$\langle\tau^{\beta}(0)\bar{\tau}^{\alpha}(0)\rangle_{\tau}=\langle\tau^{\alpha}(0)\bar{\tau}^{\alpha}(0)\rangle_{\tau}\delta_{\alpha\beta}$,
we can identify

$\langle\tau^{\alpha}(0)\bar{\tau}^{\alpha}(0)\rangle_{\tau}=\frac{1}{A_{L=0}^{v^{\alpha}}}$
(92)

as defined above in Eq. (86). Since Eq. (86) is nearly independent of $L$
Deutsch et al. (2018a); Cayuso et al. (2021), one can simply show that Eq.
(91) is equivalent to Eq. (85) at low $L$, and also equivalent to the map
based quadratic estimator presented in Refs. Deutsch et al. (2018a); Cayuso et
al. (2021) (with the same caveats as for the harmonic space estimator). In the
low signal-to-noise regime, Eq. (91) may be a simpler implementation of the
quadratic estimator used in previous work. We leave further investigation of
this point to future work.

##  Appendix D Mean and variance of the MaxL estimator

###  D.1 Harmonic space estimator

The mean of the harmonic space estimator is

$\langle\hat{v}_{LM}^{\alpha}\rangle_{\Theta}=\sum_{L^{\prime}M^{\prime}}^{L^{\prime}_{\rm{max}}}\sum_{\beta=0}^{N-1}\left([\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{-1}\right)^{\alpha\beta}_{LM;L^{\prime}M^{\prime}}\langle\left[\boldsymbol{\tau}^{\dagger}(C^{\mathrm{pCMB}})^{-1}\Theta\right]_{L^{\prime}M^{\prime}}^{\beta}\rangle_{\Theta}.$
(93)

To calculate this, we must evaluate:

$\displaystyle\langle\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\Theta\right]_{L^{\prime}M^{\prime}}^{\alpha}\rangle_{\Theta}$
$\displaystyle=$ $\displaystyle\int d\Omega\ \tau^{\alpha}(\hat{n})\
\langle\bar{\Theta}(\hat{n})\rangle_{\Theta}\
Y_{L^{\prime}M^{\prime}}^{*}(\hat{n}).$ (94)

The observed temperature anisotropies consist of the zero-mean primary CMB and
the kSZ:

$\displaystyle\langle\bar{\Theta}(\hat{n})\rangle_{\Theta}=\langle\bar{\Theta}^{\mathrm{pCMB}}(\hat{n})\rangle_{\Theta}+\bar{\Theta}^{\rm{kSZ}}(\hat{n})=\bar{\Theta}^{\rm{kSZ}}(\hat{n})$
(95)

Expanding the kSZ anisotropies into coarse- and fine-grained components, we
have

$\displaystyle\langle\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\Theta\right]_{L^{\prime}M^{\prime}}^{\alpha}\rangle_{\Theta}$
$\displaystyle=$ $\displaystyle\int d\Omega\ \tau^{\alpha}(\hat{n})\
\left[\sum_{\beta=0}^{N-1}\sum_{LM}^{L_{\rm{max}}}\bar{\tau}^{\beta}(\hat{n})v^{\beta}_{LM}Y_{LM}(\hat{n})\right.$
$\displaystyle\left.+\sum_{b=N}^{\infty}\sum_{LM}^{L_{\rm{max}}}\bar{\tau}^{b}(\hat{n})v^{b}_{LM}Y_{LM}(\hat{n})+\sum_{b=0}^{\infty}\sum_{LM=L_{max}}^{\infty}\bar{\tau}^{b}(\hat{n})v^{b}_{LM}Y_{LM}(\hat{n})\right]\
Y_{L^{\prime}M^{\prime}}^{*}(\hat{n})$ $\displaystyle\simeq$
$\displaystyle\sum_{\beta=0}^{N-1}\sum_{LM}^{L_{\rm{max}}}[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{\alpha\beta}_{L^{\prime}M^{\prime};LM}v_{LM}^{\beta}+\sum_{b=N}^{\infty}\sum_{LM}^{L_{\rm{max}}}[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{\alpha
b}_{L^{\prime}M^{\prime};LM}v_{LM}^{b}$

where we assumed that the dominant contributions from the velocity field are
at low $L$. We therefore have

$\langle\hat{v}_{LM}^{\alpha}\rangle_{\Theta}=v_{LM}^{\alpha}+\beta_{LM}^{\alpha}$
(96)

where

$\beta_{LM}^{\alpha}\equiv\sum_{L^{\prime}M^{\prime}}^{L_{\rm{max}}}\sum_{\gamma=0}^{N-1}\sum_{c=N}^{\infty}\sum_{L^{\prime\prime}M^{\prime\prime}}^{L_{\rm{max}}}\left([\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{-1}\right)^{\alpha\gamma}_{LM;L^{\prime}M^{\prime}}[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{\gamma
c}_{L^{\prime}M^{\prime};L^{\prime\prime}M^{\prime\prime}}v_{L^{\prime\prime}M^{\prime\prime}}^{c}.$
(97)

So, the estimator is unbiased up to an additive bias $\beta_{LM}^{\alpha}$,
which depends on the fine-grained radial velocity modes.

Moving to the variance:

$\displaystyle\langle\hat{v}^{\alpha}_{LM}\hat{v}^{\beta}_{L^{\prime}M^{\prime}}\rangle_{\Theta}$
$\displaystyle=$ $\displaystyle\sum_{L_{a}M_{a};L_{b}M_{b}}^{L_{\rm
max}}\sum_{\gamma;\delta=0}^{N-1}\left([\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{-1}\right)^{\alpha\gamma}_{LM;L_{a}M_{a}}\left([\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{-1}\right)^{\beta\delta}_{L^{\prime}M^{\prime};L_{b}M_{b}}$
(98) $\displaystyle\times$ $\displaystyle\int\int d\Omega d\Omega^{\prime}\
\tau^{\gamma}(\hat{n})\tau^{\delta}(\hat{n}^{\prime})\
\langle\bar{\Theta}(\hat{n})\bar{\Theta}(\hat{n}^{\prime})\rangle_{\Theta}\
Y_{L_{a}M_{a}}(\hat{n})Y_{L_{b}M_{b}}(\hat{n}^{\prime}).$

There is a contribution to the temperature correlation function from the
primary CMB and from the kSZ effect:

$\langle\bar{\Theta}(\hat{n})\bar{\Theta}(\hat{n}^{\prime})\rangle_{\Theta}=\sum_{\ell
m}\frac{1}{C_{\ell}^{\mathrm{pCMB}}}Y_{\ell m}(\hat{n})Y_{\ell
m}(\hat{n}^{\prime})+\bar{\Theta}^{\rm{kSZ}}(\hat{n})\bar{\Theta}^{\rm{kSZ}}(\hat{n}^{\prime}).$
(99)

The first term will give rise to the Gaussian estimator noise and the second
will describe our signal and biases arising due to our coarse-graining of the
maximum likelihood condition:

$\langle\hat{v}^{\alpha}_{LM}\hat{v}^{\beta}_{L^{\prime}M^{\prime}}\rangle_{\Theta}=(C^{vv})_{LM;L^{\prime}M^{\prime}}^{\alpha\beta}+N_{LM;L^{\prime}M^{\prime}}^{\alpha\beta}.$
(100)

The reconstruction noise $N_{LM;L^{\prime}M^{\prime}}^{\alpha\beta}$ arises
due to the contribution from the primary CMB:

$\displaystyle N_{LM;L^{\prime}M^{\prime}}^{\alpha\beta}$ $\displaystyle=$
$\displaystyle\sum_{L_{a}M_{a};L_{b}M_{b};\ell
m}\sum_{\gamma;\delta}\left([\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{-1}\right)^{\alpha\gamma}_{LM;L_{a}M_{a}}\left([\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{-1}\right)^{\beta\delta}_{L^{\prime}M^{\prime};L_{b}M_{b}}$
$\displaystyle\times$ $\displaystyle\int\int d\Omega d\Omega^{\prime}\
\tau^{\gamma}(\hat{n})\tau^{\delta}(\hat{n}^{\prime})\frac{1}{C_{\ell}^{\mathrm{pCMB}}}\
Y_{\ell m}(\hat{n})Y_{\ell
m}(\hat{n}^{\prime})Y_{L_{a}M_{a}}(\hat{n})Y_{L_{b}M_{b}}(\hat{n}^{\prime})$
$\displaystyle\simeq$
$\displaystyle\sum_{L_{a}M_{a};L_{b}M_{b}}\sum_{\gamma;\delta}\left([\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{-1}\right)^{\alpha\gamma}_{LM;L_{a}M_{a}}\left([\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{-1}\right)^{\beta\delta}_{L^{\prime}M^{\prime};L_{b}M_{b}}$
$\displaystyle\times$ $\displaystyle\int d\Omega\
\tau^{\gamma}(\hat{n})\bar{\tau}^{\delta}(\hat{n})\
Y_{L_{a}M_{a}}(\hat{n})Y_{L_{b}M_{b}}(\hat{n})$ (102) $\displaystyle=$
$\displaystyle\sum_{L_{a}M_{a};L_{b}M_{b}}\sum_{\gamma;\delta}\left([\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{-1}\right)^{\alpha\gamma}_{LM;L_{a}M_{a}}\left([\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{-1}\right)^{\beta\delta}_{L^{\prime}M^{\prime};L_{b}M_{b}}$
(103) $\displaystyle\times$
$\displaystyle\left(\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}\right)^{\gamma\delta}_{L_{b}M_{b};L_{a}M_{a}}$
(104) $\displaystyle=$
$\displaystyle\left([\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{-1}\right)^{\alpha\beta}_{LM;L^{\prime}M^{\prime}}$
(105)

where in going from the first to second lines we have assumed that
$C_{\ell}^{\mathrm{pCMB}}$ does not vary too quickly on scales of order
$\Delta\ell=L$ and used the identity

$\int d^{2}n\
Y_{\ell_{a}m_{a}}Y_{\ell_{b}m_{b}}Y_{\ell_{c}m_{c}}Y_{\ell_{d}m_{d}}=\sum_{\ell
m}\left(\int d^{2}n\ Y_{\ell_{a}m_{a}}Y_{\ell_{b}m_{b}}Y_{\ell
m}\right)\left(\int d^{2}n^{\prime}\ Y_{\ell_{c}m_{c}}Y_{\ell_{d}m_{d}}Y_{\ell
m}\right).$ (106)

For the signal and bias due to ‘fine modes’ we have:

$\displaystyle(C^{vv})_{LM;L^{\prime}M^{\prime}}^{\alpha\beta}$
$\displaystyle=$
$\displaystyle\sum_{L_{a}M_{a};L_{b}M_{b}}\sum_{\gamma;\delta}\left([\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{-1}\right)^{\alpha\gamma}_{LM;L_{a}M_{a}}\left([\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{-1}\right)^{\beta\delta}_{L^{\prime}M^{\prime};L_{b}M_{b}}$
$\displaystyle\times$ $\displaystyle\int\int d\Omega d\Omega^{\prime}\
\tau^{\gamma}(\hat{n})\tau^{\delta}(\hat{n}^{\prime})\
\bar{\Theta}^{\rm{kSZ}}(\hat{n})\bar{\Theta}^{\rm{kSZ}}(\hat{n}^{\prime})\
Y_{L_{a}M_{a}}(\hat{n})Y_{L_{b}M_{b}}(\hat{n}^{\prime})$ $\displaystyle\simeq$
$\displaystyle\sum_{L_{a}M_{a};L_{b}M_{b}}\sum_{\gamma;\delta;e;f}\left([\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{-1}\right)^{\alpha\gamma}_{LM;L_{a}M_{a}}\left([\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{-1}\right)^{\beta\delta}_{L^{\prime}M^{\prime};L_{b}M_{b}}$
$\displaystyle\times$ $\displaystyle\int\int d\Omega d\Omega^{\prime}\
\tau^{\gamma}(\hat{n})\tau^{\delta}(\hat{n}^{\prime})\bar{\tau}^{e}(\hat{n})v^{e}(\hat{n})\bar{\tau}^{f}(\hat{n}^{\prime})v^{f}(\hat{n}^{\prime})\
Y_{L_{a}M_{a}}(\hat{n})Y_{L_{b}M_{b}}(\hat{n}^{\prime})$ $\displaystyle=$
$\displaystyle\sum_{L_{a}M_{a};L_{b}M_{b}}\sum_{\gamma;\delta;e;f}\left([\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{-1}\right)^{\alpha\gamma}_{LM;L_{a}M_{a}}\left([\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{-1}\right)^{\beta\delta}_{L^{\prime}M^{\prime};L_{b}M_{b}}$
$\displaystyle\times$
$\displaystyle\sum_{L_{c}M_{c};L_{d}M_{d}}v_{L_{c}M_{c}}^{e}v_{L_{d}M_{d}}^{f}\int
d\Omega\
\tau^{\gamma}(\hat{n})\bar{\tau}^{e}(\hat{n})Y_{L_{a}M_{a}}(\hat{n})Y_{L_{c}M_{c}}(\hat{n})\int
d\Omega^{\prime}\tau^{\delta}(\hat{n}^{\prime})\bar{\tau}^{f}(\hat{n}^{\prime})Y_{L_{b}M_{b}}(\hat{n}^{\prime})Y_{L_{d}M_{d}}(\hat{n}^{\prime})$
$\displaystyle=$
$\displaystyle\sum_{L_{a}M_{a};L_{b}M_{b}}\sum_{L_{c}M_{c};L_{d}M_{d}}\sum_{\gamma;\delta;e;f}v_{L_{c}M_{c}}^{e}v_{L_{d}M_{d}}^{f}\left([\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{-1}\right)^{\alpha\gamma}_{LM;L_{a}M_{a}}\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}\right]^{\gamma
e}_{L_{a}M_{a};L_{c}M_{c}}$ $\displaystyle\times$
$\displaystyle\left([\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{-1}\right)^{\beta\delta}_{L^{\prime}M^{\prime};L_{b}M_{b}}\left[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}\right]^{\delta
f}_{L_{b}M_{b};L_{d}M_{d}}$ $\displaystyle=$
$\displaystyle(v_{LM}^{\alpha}+\beta_{LM}^{\alpha})(v_{L^{\prime}M^{\prime}}^{\beta}+\beta_{L^{\prime}M^{\prime}}^{\beta}).$
(108)

The main assumption made in the derivation of the estimator mean and variance
was that $C_{\ell}^{\mathrm{pCMB}}$ does not vary too quickly on scales of
order $\Delta\ell=L$, which can always be satisfied for sufficiently small
$L$. For modest values of $L\sim 100$, we expect this to be a reasonable
assumption at large $\ell$ where the primary CMB is smooth and dominated by
lensing, instrumental noise, or reionization kSZ (depending on the
instrument).

###  D.2 Map space estimator

We now compute the mean and variance of the map space estimator. Starting with
the mean:

$\displaystyle\langle\hat{v}_{I}^{\alpha}\rangle_{\Theta}$ $\displaystyle=$
$\displaystyle\sum_{\beta=0}^{N-1}\boldsymbol{\mathcal{W}}_{I;\alpha\beta}{}^{-1}\int
d\Omega\ W_{I}(\hat{n})\tau^{\beta}(\hat{n})\bar{\Theta}^{\rm kSZ}(\hat{n}),$
(109) $\displaystyle=$
$\displaystyle\sum_{\beta=0}^{N-1}\boldsymbol{\mathcal{W}}_{I;\alpha\beta}{}^{-1}\int
d\Omega\
W_{I}(\hat{n})\tau^{\beta}(\hat{n})\sum_{c=0}^{\infty}\bar{\tau}^{c}(\hat{n})\left(v^{c}_{I}+\Delta
v^{c}_{I}(\hat{n})\right)$ (110) $\displaystyle=$
$\displaystyle\sum_{\beta;\gamma=0}^{N-1}v^{\gamma}_{I}\boldsymbol{\mathcal{W}}_{I;\alpha\beta}{}^{-1}\
\boldsymbol{\mathcal{W}}_{I;\gamma\beta}+\sum_{\beta;\gamma=0}^{N-1}\boldsymbol{\mathcal{W}}_{I;\alpha\beta}{}^{-1}\int
d\Omega\ W_{I}(\hat{n})\tau^{\beta}(\hat{n})\bar{\tau}^{\gamma}(\hat{n})\Delta
v_{I}^{\gamma}(\hat{n})$ (111) $\displaystyle+$
$\displaystyle\sum_{\beta=0}^{N-1}\sum_{c=N}^{\infty}\boldsymbol{\mathcal{W}}_{I;\alpha\beta}{}^{-1}\int
d\Omega\
W_{I}(\hat{n})\tau^{\beta}(\hat{n})\bar{\tau}^{c}(\hat{n})v^{c}(\hat{n})$
(112) $\displaystyle=$ $\displaystyle v_{I}^{\alpha}+\beta_{I}^{\alpha}$ (113)

where we have defined

$\displaystyle\beta_{I}^{\alpha}$ $\displaystyle=$
$\displaystyle\sum_{\beta;\gamma=0}^{N-1}\boldsymbol{\mathcal{W}}_{I;\alpha\beta}{}^{-1}\int
d\Omega\ W_{I}(\hat{n})\tau^{\beta}(\hat{n})\bar{\tau}^{\gamma}(\hat{n})\Delta
v^{\gamma}_{I}(\hat{n})$ (114) $\displaystyle+$
$\displaystyle\sum_{\beta=0}^{N-1}\sum_{c=N}^{\infty}\boldsymbol{\mathcal{W}}_{I;\alpha\beta}{}^{-1}\int
d\Omega\
W_{I}(\hat{n})\tau^{\beta}(\hat{n})\bar{\tau}^{c}(\hat{n})v^{c}(\hat{n}).$

Moving to the estimator variance, we have:

$\displaystyle\langle\hat{v}_{I}^{\alpha}\hat{v}_{J}^{\beta}\rangle_{\Theta}$
$\displaystyle=$
$\displaystyle\sum_{\gamma;\delta=0}^{N-1}\boldsymbol{\mathcal{W}}_{I;\alpha\gamma}{}^{-1}\boldsymbol{\mathcal{W}}_{J;\beta\delta}{}^{-1}\int
d\Omega d\Omega^{\prime}\
W_{I}(\hat{n})W_{J}(\hat{n}^{\prime})\tau^{\gamma}(\hat{n})\tau^{\delta}(\hat{n}^{\prime})\langle\bar{\Theta}(\hat{n})\bar{\Theta}(\hat{n}^{\prime})\rangle_{\Theta}$
(115) $\displaystyle=$
$\displaystyle\sum_{\gamma;\delta=0}^{N-1}\boldsymbol{\mathcal{W}}_{I;\alpha\gamma}{}^{-1}\boldsymbol{\mathcal{W}}_{J;\beta\delta}{}^{-1}\left[\int
d\Omega\
W_{I}(\hat{n})\tau^{\gamma}(\hat{n})\bar{\Theta}^{\rm{kSZ}}(\hat{n})\int
d\Omega^{\prime}\
W_{J}(\hat{n}^{\prime})\tau^{\delta}(\hat{n}^{\prime})\bar{\Theta}^{\rm{kSZ}}(\hat{n}^{\prime})\right.$
$\displaystyle+$ $\displaystyle\left.\sum_{\ell m}\int d\Omega\
W_{I}(\hat{n})\tau^{\gamma}(\hat{n})Y_{\ell m}(\hat{n})\int d\Omega^{\prime}\
W_{J}(\hat{n}^{\prime})\frac{\tau^{\delta}(\hat{n}^{\prime})}{C_{\ell}^{\mathrm{pCMB}}}Y_{\ell
m}(\hat{n}^{\prime})\right]$ $\displaystyle\simeq$
$\displaystyle\left(v_{I}^{\alpha}+\beta_{I}^{\alpha}\right)\left(v_{J}^{\beta}+\beta_{J}^{\beta}\right)+\sum_{\gamma;\delta=0}^{N-1}\boldsymbol{\mathcal{W}}_{I;\alpha\gamma}{}^{-1}\boldsymbol{\mathcal{W}}_{J;\beta\delta}{}^{-1}\int
d\Omega\
W_{I}(\hat{n})W_{J}(\hat{n})\tau^{\gamma}(\hat{n})\bar{\tau}^{\delta}(\hat{n})$
$\displaystyle=$
$\displaystyle\left(v_{I}^{\alpha}+\beta_{I}^{\alpha}\right)\left(v_{J}^{\beta}+\beta_{J}^{\beta}\right)+\frac{1}{\Delta\theta^{2}}\boldsymbol{\mathcal{W}}_{I;\alpha\beta}{}^{-1}\delta_{IJ}$

where in going from the second to the third line, we have made the
approximation that $C_{\ell}^{\mathrm{pCMB}}$ is not too rapidly varying with
scale, as we did in the derivation of the harmonic space estimator variance.

##  Appendix E Optical depth bias

When we only have access to an estimate $\hat{\tau}$ of the differential
optical depth, the harmonic space estimator is given by

$\hat{v}_{LM}^{\alpha}=\sum_{L^{\prime}M^{\prime}}^{L^{\prime}_{\rm{max}}}\sum_{\beta=0}^{N-1}\left([\boldsymbol{\hat{\tau}}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\hat{\tau}}]^{-1}\right)^{\alpha\beta}_{LM;L^{\prime}M^{\prime}}(\boldsymbol{\hat{\tau}}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\Theta)_{L^{\prime}M^{\prime}}^{\beta}$
(116)

Computing the estimator mean, we have:

$\displaystyle\langle\hat{v}_{LM}^{\alpha}\rangle_{\Theta}$ $\displaystyle=$
$\displaystyle\sum_{L^{\prime}M^{\prime}}^{L^{\prime}_{\rm{max}}}\sum_{\beta=0}^{N-1}\left([\boldsymbol{\hat{\tau}}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\hat{\tau}}]^{-1}\right)^{\alpha\beta}_{LM;L^{\prime}M^{\prime}}\langle(\boldsymbol{\hat{\tau}}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\Theta)_{L^{\prime}M^{\prime}}^{\beta}\rangle_{\Theta}$
(117) $\displaystyle=$
$\displaystyle\left[\sum_{L^{\prime\prime}M^{\prime\prime}}^{L_{\rm{max}}}\sum_{\gamma=0}^{N-1}R^{\alpha\gamma}_{LM;L^{\prime\prime}M^{\prime\prime}}v_{L^{\prime\prime}M^{\prime\prime}}^{\gamma}\right]+\hat{\beta}^{\alpha}_{LM}$
(118)

where

$R^{\alpha\gamma}_{LM;L^{\prime\prime}M^{\prime\prime}}=\sum_{L^{\prime}M^{\prime}}^{L^{\prime}_{\rm{max}}}\sum_{\beta=0}^{N-1}\left([\boldsymbol{\hat{\tau}}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\hat{\tau}}]^{-1}\right)^{\alpha\beta}_{LM;L^{\prime}M^{\prime}}[\boldsymbol{\hat{\tau}}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{\beta\gamma}_{L^{\prime}M^{\prime};L^{\prime\prime}M^{\prime\prime}}$
(119)

and

$\hat{\beta}_{LM}^{\alpha}\equiv\sum_{L^{\prime}M^{\prime}}^{L_{\rm{max}}}\sum_{\gamma=0}^{N-1}\sum_{c=N}^{\infty}\sum_{L^{\prime\prime}M^{\prime\prime}}^{L_{\rm{max}}}\left([\boldsymbol{\hat{\tau}}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\hat{\tau}}]^{-1}\right)^{\alpha\gamma}_{LM;L^{\prime}M^{\prime}}[\boldsymbol{\hat{\tau}}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{\gamma
c}_{L^{\prime}M^{\prime};L^{\prime\prime}M^{\prime\prime}}v_{L^{\prime\prime}M^{\prime\prime}}^{c}.$
(120)

When $\hat{\tau}=\tau$, this is simply the identity matrix, and we recover an
unbiased estimator mean. In reality, an imperfect estimate of the differential
optical depth will result in the mixing of information and therefore a bias.
It is important to characterize the magnitude and properties of this bias to
obtain useful cosmological information from the reconstruction.

To explore the optical depth bias a bit more, suppose we have a pretty good
reconstruction of the differential optical depth, such that

$\hat{\tau}^{\alpha}(\hat{n})=\tau^{\alpha}(\hat{n})+n_{\tau}^{\alpha}(\hat{n});\
\ \ n_{\tau}^{\alpha}(\hat{n})\ll\tau^{\alpha}(\hat{n}).$ (121)

We can approximate

$[\boldsymbol{\hat{\tau}}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\hat{\tau}}]^{-1}\simeq\left[1-[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{-1}\cdot[\boldsymbol{n_{\tau}}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{n_{\tau}}]\right]\cdot[\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\tau]^{-1}$
(122)

and therefore

$R^{\alpha\gamma}_{LM;L^{\prime\prime}M^{\prime\prime}}\simeq\delta_{\alpha\gamma}\delta_{LM;L^{\prime\prime}M^{\prime\prime}}-\sum_{L^{\prime}M^{\prime}}^{L^{\prime}_{\mathrm{max}}}\sum_{\beta=0}^{N-1}\left([\boldsymbol{\tau}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{\tau}]^{-1}\right)^{\alpha\beta}_{LM;L^{\prime}M^{\prime}}[\boldsymbol{n_{\tau}}^{\dagger}(\boldsymbol{C}^{\mathrm{pCMB}})^{-1}\boldsymbol{n_{\tau}}]^{\beta\gamma}_{L^{\prime}M^{\prime};L^{\prime\prime}M^{\prime\prime}}.$
(123)

Taking an ensemble average over $\tau$ and $n_{\tau}$, and using Eq. (75), we
have

$\langle
R^{\alpha\gamma}_{LM;L^{\prime\prime}M^{\prime\prime}}\rangle\simeq\delta_{\alpha\gamma}\delta_{LM;L^{\prime\prime}M^{\prime\prime}}-\delta_{LM;L^{\prime\prime}M^{\prime\prime}}\sum_{\beta}^{\beta_{\mathrm{max}}}\langle\tau^{\alpha}(0)\bar{\tau}^{\beta}(0)\rangle^{-1}\langle
n_{\tau}^{\beta}(0)\bar{n}_{\tau}^{\gamma}(0)\rangle$ (124)

where $\langle\tau^{\alpha}(0)\bar{\tau}^{\beta}(0)\rangle^{-1}$ is the bin-
space inverse of Eq. (75) and $\langle
n_{\tau}^{\beta}(0)\bar{n}_{\tau}^{\gamma}(0)\rangle$ is

$\langle
n_{\tau}^{\beta}(0)\bar{n}_{\tau}^{\gamma}(0)\rangle=\sum_{\ell}\frac{2\ell+1}{4\pi}\frac{C_{\ell}^{n_{\tau}^{\beta}n_{\tau}^{\gamma}}}{C_{\ell}^{\mathrm{pCMB}}}.$
(125)

In particular, note that on average in this limit the optical depth bias is
independent of scale $L^{\prime\prime}M^{\prime\prime}$, and small when the
differential optical depth is reconstructed well.

For the map space estimator, we have:

$\hat{v}_{I}^{\alpha}=\sum_{\beta=0}^{N-1}\left[\int d\Omega\
W_{I}(\hat{n})\hat{\tau}^{\beta}(\hat{n})\hat{\bar{\tau}}^{\alpha}(\hat{n})\right]^{-1}\int
d\Omega\ W_{I}(\hat{n})\hat{\tau}^{\beta}(\hat{n})\bar{\Theta}(\hat{n}).$
(126)

Computing the estimator mean:

$\langle\hat{v}_{I}^{\alpha}\rangle=\sum_{\gamma=0}^{N-1}R^{\alpha\gamma}_{I}v_{I}^{\gamma}+\hat{\beta}_{I}^{\alpha}$
(127)

where

$R^{\alpha\gamma}_{I}\equiv\sum_{\beta=0}^{N-1}\left[\int d\Omega\
W_{I}(\hat{n})\hat{\tau}^{\beta}(\hat{n})\hat{\bar{\tau}}^{\alpha}(\hat{n})\right]^{-1}\int
d\Omega\
W_{I}(\hat{n})\hat{\tau}^{\beta}(\hat{n})\bar{\tau}^{\gamma}(\hat{n})$ (128)

and

$\displaystyle\hat{\beta}_{I}^{\alpha}$ $\displaystyle=$
$\displaystyle\sum_{\beta;\gamma=0}^{N-1}\left[\int d\Omega\
W_{I}(\hat{n})\hat{\tau}^{\beta}(\hat{n})\hat{\bar{\tau}}^{\alpha}(\hat{n})\right]^{-1}\int
d\Omega\
W_{I}(\hat{n})\hat{\tau}^{\beta}(\hat{n})\bar{\tau}^{\gamma}(\hat{n})\Delta
v^{\gamma}_{I}(\hat{n})$ (129) $\displaystyle+$
$\displaystyle\sum_{\beta=0}^{N-1}\sum_{c=N}^{\infty}\left[\int d\Omega\
W_{I}(\hat{n})\hat{\tau}^{\beta}(\hat{n})\hat{\bar{\tau}}^{\alpha}(\hat{n})\right]^{-1}\int
d\Omega\
W_{I}(\hat{n})\tau^{\beta}(\hat{n})\bar{\tau}^{c}(\hat{n})v^{c}(\hat{n}).$

Assuming a well-modelled optical depth, as in Eq. (121), a similar derivation
as for the harmonic estimator yields:

$\langle
R^{\alpha\gamma}_{I}\rangle\simeq\delta_{\alpha\gamma}-\sum_{\beta=0}^{N-1}\langle\tau^{\alpha}(0)\bar{\tau}^{\beta}(0)\rangle^{-1}\langle
n_{\tau}^{\beta}(0)\bar{n}_{\tau}^{\gamma}(0)\rangle.$ (130)

So, the optical depth bias in this limit is independent of the pixel index
$I$.