# Probing cosmic strings by reconstructing polarization rotation of the cosmic
microwave background

Weichen Winston Yin  Liang Dai  and Simone Ferraro

###### Abstract

Cosmic birefringence—the rotation of the polarization of the cosmic microwave
background (CMB) photons as they travel to the Earth—is a smoking gun for
axion-like particles (ALPs) that interact with the photon. It has recently
been suggested that topological defects in the ALP field called cosmic strings
can cause polarization rotation in quantized amounts that are proportional to
the electromagnetic chiral anomaly coefficient $\mathcal{A}$, which holds
direct information about physics at very high energies. In this work, we study
the detectability of a random network of cosmic strings through estimating
rotation using quadratic estimators (QEs). We show that the QE derived from
the maximum likelihood estimator is equivalent to the recently proposed
global-minimum-variance QE; the classic Hu-Okamoto QE equals the global-
minimum-variance QE under special conditions, but is otherwise still nearly
globally optimal. We calculate the sensitivity of QEs to cosmic birefringence
from string networks, for the Planck satellite mission, as well as for third-
and fourth-generation ground-based CMB experiments. Using published binned
rotation power spectrum derived from the Planck 2015 polarization data, we
obtain a constraint $\mathcal{A}^{2}\,\xi_{0}<0.93$ at the 95% confidence
level, where $\xi_{0}$ is the total length of strings in units of the Hubble
scale per Hubble volume, for a phenomenological but reasonable string network
model describing a continuous distribution of string sizes. Looking forward,
experiments such as the Simons Observatory and CMB-S4 will either discover or
falsify the existence of an ALP string network for the theoretically plausible
parameter space $\mathcal{A}^{2}\,\xi_{0}\gtrsim 0.01$.

##  1 Introduction

The intensity and polarization anistropies of the cosmic microwave background
(CMB) are the most stringently tested physical observables in cosmology. At
the leading order, the anisotropies are Gaussian random quantities and follow
isotropic power spectra that are precisely predictable from physical models of
the primordial Universe. Nearly originated from the cosmic horizon, the CMB
photons traverse the longest distances in the Universe without major
absorption or scattering.

The above two reasons render the CMB a sensitive probe of cosmic
birefringence—the coherent rotation of the polarization planes as CMB photons
travel along the line of sight. This effect can be mediated by an ultralight
axion-like particles (ALPs, or simply “axions” in this work) that couples to
the photon via a Chern-Simons coupling [1, 2, 3]. Such ALPs are a ubiquitous
outcome of string theory compactifications [4]; they may account for the dark
matter [5], or hold clues to the mystery of the Dark Energy [6, 7]. Recently,
Ref. [8] pointed out that topological defects in the ALP field, termed cosmic
strings, can imprint birefringence patterns in the CMB in a quantized manner.
Remarkably, the observable rotation angle is insensitive to the ALP-photon
coupling constant, which is unconstrained and subject to renormalization. In
fact, the rotation angle varies in quantized amounts as the line of sight
crosses individual strings, which is determined by the electromagnetic chiral
anomaly coefficient $\mathcal{A}$, free from renormalization [9]. Not only
does $\mathcal{A}$ offer insight into the smallest unit of electric charge
[10], it is also useful for distinguishing between beyond-the-Standard-Model
theories, whose matter content and charge assignment at very high energies
lead to predictions for the value of $\mathcal{A}$ [8].

Recently, the study of CMB birefringence has heated up after tentative
evidence for birefringence that is highly coherent across the sky was reported
in a refined analysis of the Planck 2018 polarization data [11], with
implications for ALP searches [12]. While an isotropic polarization rotation
would be evidence for symmetry-breaking phenomena beyond the Standard Model,
it is the power spectrum of direction-dependent rotations that will inform us
about the spatial structure of birefringence sources.

While the anisotropic rotation induced by the strings is non-Gaussian by
nature, in this work we study detectability using quadratic estimators which
are optimal for Gaussian random fields. In this method, one seeks statistical
correlation between Fourier modes of the anisotropies with different Fourier
wave vectors. This method is appropriate for the scenario of a large number of
strings within the observable Universe, but should also offer decent detection
statistics in the case of relatively few strings. This simplification may be
adequate if one aims for a first statistical detection of rotation, before
characterization of the non-Gaussian behaviors of the rotation pattern is
carried out.

Quadratic estimators for the purpose of estimating direction-dependent
birefringence were first investigated in Refs. [13, 14] in the full-sky
formalism and in Ref. [15] in the flat-sky approximation. This method was used
to derive bounds on the rotation power spectrum in the WMAP data [16]. In this
work, we compare a few differently constructed quadratic estimators in general
as well as in the context of cosmic birefringence. We study the detectability
of direction-dependent birefringence as a result of a network of ALP cosmic
strings [8, 17], for CMB anisotropy data collected by the Planck mission as
well as forthcoming data to be collected by third and fourth generation
ground-based CMB experiments. Instead of using the $EB$ and/or $TB$ estimators
alone (as in Refs. [15] and [8]), we maximize information gain by including
all quadratic combinations of CMB anisotropies. We give signal-to-noise ratio
(SNR) forecasts for measuring the parameter combination
$\mathcal{A}^{2}\,\xi_{0}$ through these CMB experiments, where $\xi_{0}$ is
the total length of strings per Hubble volume in units of the Hubble length,
as described in Section 7.

In the absence of cosmic birefringence, the intervening large-scale structure
of the Universe induces cross-correlation between Fourier modes of the CMB
anisotropies via gravitational lensing. Nevertheless, we will show that
estimation of polarization rotation is nearly orthogonal to the estimation of
the lensing potential, as was also discussed in Ref. [14]. Therefore, for the
purpose of detecting polarization rotation, it is a good approximation to
treat the lensed anistropies as Gaussian random fields with statistical
isotropy and to completely neglect mode correlations induced by lensing.

We study three ways to construct quadratic estimators for the rotation field
using primary anisotropies. The first one is the famous Hu-Okamoto (HO)
estimator [18], which has been widely employed in weak lensing reconstruction
of the CMB. In the second way, the so-called global-minimum-variance (GMV)
estimator is constructed following the recent work of Ref. [19], which
improves upon the HO estimator. For the third one, we derive an optimal
quadratic estimator following the maximization of the full likelihood function
for the anisotropy observables.

We explicitly show that the third quadratic estimator, obtained as an
approximation to the maximal likelihood estimator, is equivalent to the GMV
estimator, as one would anticipate based on the definition of the GMV
estimator. Although the GMV estimator is in general distinct from the HO
estimator and is mathematically guaranteed to be at least as sensitive as the
latter [19], we show that the two are equivalent when rotation estimators are
constructed only from the $E$\- and $B$-mode polarization anisotropies, which
includes the dominant $EB$ estimators. We phrase a general sufficient
condition for such equivalence. Furthermore, we numerically find that if the
temperature anisotropies are also included, then the difference in the
reconstruction noise between the GMV and HO estimators for the rotation field
is insignificant at both current and forthcoming CMB experiments, with less
than 0.5% improvement from GMV across the range of angular scales.

We find that information about the plausible cosmic string network is
dominated by rotation reconstruction at large scales $L\lesssim 100$, and that
data from the Planck mission should already have constraining power for a
significant part of the theoretically favored parameter space,
$\mathcal{A}^{2}\,\xi_{0}\gtrsim\mathcal{O}(1)$. For the entire theoretically
plausible parameter range $\mathcal{A}^{2}\,\xi_{0}\gtrsim 0.01$, we find that
future CMB experiments will have the definitive statistical power to either
substantiate or falsify the existence of cosmic strings from ALPs that couple
to electromagnetism.

Using published birefringence power spectrum extracted from Planck 2015 data
[20], we place upper bounds on the parameter $\mathcal{A}^{2}\,\xi_{0}$. For
two phenomenological string network models developed in [17], we find that
$\mathcal{A}^{2}\,\xi_{0}<8.0$ if the Universe is populated by string loops of
the same but arbitrary size, and $\mathcal{A}^{2}\,\xi_{0}<0.93$ if, more
realistically, the string loops have a continuous distribution of sizes, both
at the 95% confidence level. These results are consistent with the absence of
a cosmic string network.

This paper is organized as follows. Section 2 reviews the construction of the
HO and GMV estimators. In Section 3, we derive a quadratic estimator from the
maximum likelihood method. In Section 4, we show that the maximum likelihood
quadratic estimator is equivalent to the GMV estimator; in addition, HO and
GMV estimators constructed from only the polarization anisotropies are
equivalent, and a general sufficient condition for such equivalence is stated.
In Section 6, we numerically compute the reconstruction noise spectra for the
HO and GMV estimators, respectively, showing that the difference is
insignificant in practice. In Section 7, we compute the overall signal-to-
noise ratio for inferring the existence of a cosmic string network. We show in
Section 5 that the bias induced by lensing on estimating rotation is at most
second order in the lensing potential, justifying the approximation of lensed
CMB anisotropies as Gaussian random fields. Finally in Section 8,
birefringence power spectra extracted from Planck 2015 is used to constrain
parameters in string network models. Concluding remarks are given in Section
9.

##  2 Quadratic estimators

As the CMB photons travel from the last scattering surface, they follow
deflected trajectories due to gravitational deflection by matter overdensity
and underdensity. Along the perturbed trajectories, the axion-photon coupling
causes a rotation of the plane of polarization. At the leading order, it
suffices to compute the rotation effect along the unperturbed line of sight.
The absence of interference between these two effects at the leading order is
discussed and justified in Section 5. As we focus on the rotation effect, we
start with the lensed but unrotated CMB anisotropies, and consider the effect
of rotation on them.

###  2.1 Notation

Throughout this paper, we assume a flat-sky approximation, in which the
observables are functions of
$\bm{\mathbf{\hat{n}}}=n_{x}\,\bm{\mathbf{\hat{x}}}+n_{y}\,\bm{\mathbf{\hat{y}}}$.
We use the tilded notation to indicate lensed but unrotated quantities: the
temperature anisotropies $\tilde{T}(\bm{\mathbf{\hat{n}}})$, and the
polarization anisotropies $\tilde{E}(\bm{\mathbf{\hat{n}}})$ and
$\tilde{B}(\bm{\mathbf{\hat{n}}})$. Their power spectra in Fourier space
$\tilde{C}^{XY}_{l}$ are defined through:

$\expectationvalue{\tilde{X}(\bm{\mathbf{l}}_{1})\tilde{Y}^{*}(\bm{\mathbf{l}}_{2})}=(2\pi)^{2}\,\delta_{D}(\bm{\mathbf{l}}_{1}-\bm{\mathbf{l}}_{2})\,\tilde{C}^{XY}_{l_{1}}.$
(2.1)

We parameterize cosmic birefringence using an anisotropic rotation angle for
the polarization $\alpha(\bm{\mathbf{\hat{n}}})$. The lensed and rotated
quantities are denoted by $T(\bm{\mathbf{\hat{n}}})$,
$E(\bm{\mathbf{\hat{n}}})$, and $B(\bm{\mathbf{\hat{n}}})$ (birefringence
leaves the temperature anisotropies unchanged). The cross correlations are
consequently modified to

$\expectationvalue{X(\bm{\mathbf{l}}_{1})Y^{*}(\bm{\mathbf{l}}_{2})}=\expectationvalue{\tilde{X}(\bm{\mathbf{l}}_{1})\tilde{Y}^{*}(\bm{\mathbf{l}}_{2})}+f_{XY}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})\,\alpha(\bm{\mathbf{l}}_{1}-\bm{\mathbf{l}}_{2})+O(\alpha^{2}).$
(2.2)

Based on the latest constraints from the CMB data [21] and predictions of
theoretical models [17], we expect $\alpha(\bm{\mathbf{\hat{n}}})$ to be small
angles (in radians). Therefore, it suffices to expand the expression up to the
linear order.

The observed anisotropies (denoted with a subscript $o$) include both the
cosmological signals and noise (denoted with a subscript $n$):

$X_{o}(\bm{\mathbf{l}})=X(\bm{\mathbf{l}})+X_{n}(\bm{\mathbf{l}}).$ (2.3)

The Fourier space power spectra can be separately defined for the observed
anisotropies and for the noise, respectively:

$\displaystyle\expectationvalue{X_{o}(\bm{\mathbf{l}}_{1})Y_{o}^{*}(\bm{\mathbf{l}}_{2})}$
$\displaystyle=(2\pi)^{2}\,\delta_{D}(\bm{\mathbf{l}}_{1}-\bm{\mathbf{l}}_{2})\,C^{XY}_{l_{1}}$
(2.4)
$\displaystyle\expectationvalue{X_{n}(\bm{\mathbf{l}}_{1})Y_{n}^{*}(\bm{\mathbf{l}}_{2})}$
$\displaystyle=(2\pi)^{2}\,\delta_{D}(\bm{\mathbf{l}}_{1}-\bm{\mathbf{l}}_{2})\,C^{XY,n}_{l_{1}}.$
(2.5)

In this paper, we concern the variance of rotation estimators under the null
hypothesis, _i.e._ $\alpha(\bm{\mathbf{\hat{n}}})=0$. Since cosmic signals and
noise are uncorrelated, we have

$C^{XY}_{l}=\tilde{C}^{XY}_{l}+C^{XY,n}_{l}.$ (2.6)

###  2.2 Rotation of polarization

Let $\tilde{Q}(\bm{\mathbf{\hat{n}}})\pm i\tilde{U}(\bm{\mathbf{\hat{n}}})$ be
the Stokes’ parameters for the lensed but unrotated polarizations. For a line-
of-sight specific polarization rotation $\alpha(\bm{\mathbf{n}})$, the
perturbed Stokes fields are given by:

$Q(\bm{\mathbf{\hat{n}}})\pm
iU(\bm{\mathbf{\hat{n}}})=\left[\tilde{Q}(\bm{\mathbf{\hat{n}}})\pm
i\tilde{U}(\bm{\mathbf{\hat{n}}})\right]\,e^{\pm
2\,i\,\alpha(\bm{\mathbf{\hat{n}}})},$ (2.7)

The $E$\- and $B$-modes can be derived from the Stokes’ parameters via [22]

$E(\bm{\mathbf{l}})\pm
iB(\bm{\mathbf{l}})=\int\mathrm{d}^{2}\bm{\mathbf{\hat{n}}}\,e^{-i\bm{\mathbf{l}}\cdot\bm{\mathbf{\hat{n}}}}\,e^{\mp
2\,i\,\varphi_{\bm{\mathbf{l}}}}\,[Q(\bm{\mathbf{\hat{n}}})\pm
iU(\bm{\mathbf{\hat{n}}})],$ (2.8)

where $\varphi_{\bm{\mathbf{l}}}$ is the angle between the vector
${\bm{\mathbf{l}}}$ and the first coordinate axis on the sky. A similar
expression can be written for the lensed but unrotated anisotropies.

We have assumed that all primary CMB anisotropies originate from the last
scattering surface, neglecting the possibility that the effects of the ALP
field may be different for polarizations produced at recombination and at
reionization. Future CMB satellites such as LiteBIRD will be able to detect
any difference in the birefringence angle between the two epochs to within
$0.05^{\circ}$ [23]. For the experiments we consider, this is much smaller
than the reconstruction noise at all but the lowest few $L$ modes (see Figure
1). Due to the limited number of modes, reionization polarization signals are
not expected to affect the overall detectability for even the lowest $L$
modes, regardless of the correct treatment of the rotation effect up to the
epoch of reionization. Having the small-scale primary anisotropies available
is sufficient for reconstructing the rotation on large angular scales.

In Table 1, we summarize the response functions (introduced in Eq. (2.2)) due
to the action of rotation.

$XY$ | $f_{XY}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})$
---|---
$TT$ | $0$
$TE$ | $-2\tilde{C}_{l_{1}}^{TE}\sin 2\varphi_{12}$
$TB$ | $2\tilde{C}_{l_{1}}^{TE}\cos 2\varphi_{12}$
$EE$ | $-2\left(\tilde{C}_{l_{1}}^{EE}-\tilde{C}_{l_{2}}^{EE}\right)\sin 2\varphi_{12}$
$EB$ | $2\left(\tilde{C}_{l_{1}}^{EE}-\tilde{C}_{l_{2}}^{BB}\right)\cos 2\varphi_{12}$
$BB$ | $-2\left(\tilde{C}_{l_{1}}^{BB}-\tilde{C}_{l_{2}}^{BB}\right)\sin 2\varphi_{12}$
Table 1: Response functions under the action of line-of-sight dependent
rotation of polarization. Here we have defined the angle
$\varphi_{12}=\varphi_{\mathbf{l}_{1}}-\varphi_{\mathbf{l}_{2}}$.

###  2.3 Hu-Okamoto estimator

Following the treatment of Hu & Okamoto for estimating weak lensing from the
CMB anisotropies [18], quadratic estimators can be constructed from the
observed quantities [15]:

$\hat{\alpha}_{XY}(\bm{\mathbf{L}})=\int\frac{\mathrm{d}^{2}\bm{\mathbf{l}}_{1}\,\mathrm{d}^{2}\bm{\mathbf{l}}_{2}}{(2\pi)^{2}}\,\delta_{D}(\bm{\mathbf{l}}_{1}-\bm{\mathbf{l}}_{2}-\bm{\mathbf{L}})\,X_{o}(\bm{\mathbf{l}}_{1})\,Y_{o}^{*}(\bm{\mathbf{l}}_{2})\,F_{XY}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2}),$
(2.9)

where $F_{XY}$ is a yet undetermined real-valued function,

and we have assumed $\bm{\mathbf{L}}\neq 0$. The $\bm{\mathbf{L}}=0$ case
corresponds to a constant rotation of the polarization on the whole sky. See
Refs. [24, 25, 26, 27, 28, 29, 11] for recent estimates. Henceforth, we will
write $\bm{\mathbf{L}}=\bm{\mathbf{l}}_{1}-\bm{\mathbf{l}}_{2}$. Following
Ref. [19], we also introduce the shorthand notation:

$\int\frac{\mathrm{d}^{2}\bm{\mathbf{l}}_{1}\,\mathrm{d}^{2}\bm{\mathbf{l}}_{2}}{(2\pi)^{2}}\,\delta_{D}(\bm{\mathbf{l}}_{1}-\bm{\mathbf{l}}_{2}-\bm{\mathbf{L}})\to\int_{\bm{\mathbf{L}}=\mathbf{l}_{1}-\mathbf{l}_{2}}.$
(2.10)

To find the coefficients $F_{XY}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})$,
three conditions need to be imposed on $\hat{\alpha}_{XY}(\bm{\mathbf{L}})$.
First, the reality condition
$\hat{\alpha}_{XY}^{*}(\bm{\mathbf{L}})=\hat{\alpha}_{XY}(-\bm{\mathbf{L}})$
is needed because the rotation field to be estimated is real-valued. This
requires
$F_{XY}(-\bm{\mathbf{l}}_{1},-\bm{\mathbf{l}}_{2})=F_{XY}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})$.
Second, the estimator is required to be unbiased,
$\expectationvalue{\hat{\alpha}_{XY}(\bm{\mathbf{L}})}=\alpha(\bm{\mathbf{L}})$,
which leads to equation

$1=\int_{\bm{\mathbf{L}}=\mathbf{l}_{1}-\mathbf{l}_{2}}f_{XY}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})\,F_{XY}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2}).$
(2.11)

Finally, subject to the above two conditions, we seek to minimize the variance
of the estimator via the following condition:

$\frac{\delta}{\delta
F_{XY}}\left[\expectationvalue{\hat{\alpha}_{XY}(\bm{\mathbf{L}})\,\hat{\alpha}_{XY}^{*}(\bm{\mathbf{L}}^{\prime})}-\expectationvalue{\hat{\alpha}_{XY}(\bm{\mathbf{L}})}\expectationvalue{\hat{\alpha}_{XY}^{*}(\bm{\mathbf{L}}^{\prime})}\right]=0.$
(2.12)

The solution to this minimization problem reads:

$\displaystyle F_{XY}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})$
$\displaystyle=\lambda_{XY}(\bm{\mathbf{L}})\,\frac{f_{XY}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})}{(1+\delta_{XY})\,C_{l_{1}}^{XX}\,C_{l_{2}}^{YY}}\quad\text{for
$XY\neq TE$},$ (2.13) $\displaystyle
F_{TE}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})$
$\displaystyle=\lambda_{TE}(\bm{\mathbf{L}})\,\frac{C^{EE}_{l_{1}}\,C^{TT}_{l_{2}}\,f_{TE}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})-C^{TE}_{l_{1}}\,C^{TE}_{l_{2}}\,f_{TE}(\bm{\mathbf{l}}_{2},\bm{\mathbf{l}}_{1})}{C^{TT}_{l_{1}}\,C^{EE}_{l_{2}}\,C^{EE}_{l_{1}}\,C^{TT}_{l_{2}}-(C^{TE}_{l_{1}}\,C^{TE}_{l_{2}})^{2}},$
(2.14)

where we define the normalization constants

$\displaystyle[\lambda_{XY}(\bm{\mathbf{L}})]^{-1}$
$\displaystyle=\int_{\bm{\mathbf{L}}=\mathbf{l}_{1}-\mathbf{l}_{2}}\frac{[f_{XY}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})]^{2}}{(1+\delta_{XY})\,C^{XX}_{l_{1}}\,C^{YY}_{l_{2}}}\quad\text{for
$XY\neq TE$},$ (2.15) $\displaystyle[\lambda_{TE}(\bm{\mathbf{L}})]^{-1}$
$\displaystyle=\int_{\bm{\mathbf{L}}=\mathbf{l}_{1}-\mathbf{l}_{2}}f_{TE}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})\,\frac{C^{EE}_{l_{1}}\,C^{TT}_{l_{2}}\,f_{TE}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})-C^{TE}_{l_{1}}\,C^{TE}_{l_{2}}\,f_{TE}(\bm{\mathbf{l}}_{2},\bm{\mathbf{l}}_{1})}{C^{TT}_{l_{1}}\,C^{EE}_{l_{2}}\,C^{EE}_{l_{1}}\,C^{TT}_{l_{2}}-(C^{TE}_{l_{1}}\,C^{TE}_{l_{2}})^{2}}.$
(2.16)

The covariance between two estimators is

$\expectationvalue{\hat{\alpha}_{XY}(\bm{\mathbf{L}})\hat{\alpha}_{X^{\prime}Y^{\prime}}^{*}(\bm{\mathbf{L}}^{\prime})}-\expectationvalue{\hat{\alpha}_{XY}(\bm{\mathbf{L}})}\expectationvalue{\hat{\alpha}_{X^{\prime}Y^{\prime}}^{*}(\bm{\mathbf{L}}^{\prime})}=(2\pi)^{2}\delta(\bm{\mathbf{L}}-\bm{\mathbf{L}}^{\prime})N_{XYX^{\prime}Y^{\prime}}(\bm{\mathbf{L}}),$
(2.17)

where

$N_{XYX^{\prime}Y^{\prime}}(\bm{\mathbf{L}})=\int_{\bm{\mathbf{L}}=\bm{\mathbf{l}}_{1}-\bm{\mathbf{l}}_{2}}F_{XY}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})\left[F_{X^{\prime}Y^{\prime}}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})C^{XX^{\prime}}_{l_{1}}C^{YY^{\prime}}_{l_{2}}+F_{X^{\prime}Y^{\prime}}(\bm{\mathbf{l}}_{2},\bm{\mathbf{l}}_{1})C^{XY^{\prime}}_{l_{1}}C^{YX^{\prime}}_{l_{2}}\right].$
(2.18)

The estimator constructed for any single pair $XY$ has a variance

$\expectationvalue{\hat{\alpha}_{XY}(\bm{\mathbf{L}})\,\hat{\alpha}_{XY}^{*}(\bm{\mathbf{L}}^{\prime})}-\expectationvalue{\hat{\alpha}_{XY}(\bm{\mathbf{L}})}\expectationvalue{\hat{\alpha}_{XY}^{*}(\bm{\mathbf{L}}^{\prime})}=(2\pi)^{2}\,\delta_{D}(\bm{\mathbf{L}}-\bm{\mathbf{L}}^{\prime})\,N_{XY}(\bm{\mathbf{L}}),$
(2.19)

where it can be shown that
$N_{XY}(\bm{\mathbf{L}})=N_{XYXY}=\lambda_{XY}(\bm{\mathbf{L}})$. We note that
$N_{XY}(\bm{\mathbf{L}})$ is only dependent on $L=|{\bm{\mathbf{L}}}|$ due to
statistical invariance under rotation of the sky. From here on, we will write
$N_{XY}(L)$ to reflect this fact.

The Hu-Okamoto (HO) estimator is constructed as the unique linear combination
of the estimators $\hat{\alpha}_{XY}(\bm{\mathbf{L}})$ corresponding to five
pairs, $TE$, $TB$, $EE$, $EB$, and $BB$ (the $TT$ estimator does not exist, as
pointwise polarization rotations have no effect on the temperature
anisotropies), such that the linear combination remains unbiased and has a
minimized variance:

$\hat{\alpha}_{HO}(\bm{\mathbf{L}})=\sum_{I}\,w_{I}(L)\,\hat{\alpha}_{I}(\bm{\mathbf{L}}),$
(2.20)

where the index $I$ ranges over unique pairs of observables $XY$. Defining the
symmetric matrix $\bm{\mathbf{N}}(L)$ with elements
$[\bm{\mathbf{N}}]_{IJ}(L)=N_{IJ}(L)$ as defined in Eq. (2.18), the weights
can be computed from

$w_{I}(L)=\frac{\sum_{J}[\bm{\mathbf{N}}^{-1}(L)]_{IJ}}{\sum_{JK}[\bm{\mathbf{N}}^{-1}(L)]_{JK}}.$
(2.21)

Under the null hypothesis that there is no polarization rotation, the HO
estimator has a variance

$\expectationvalue{\hat{\alpha}_{HO}(\bm{\mathbf{L}})\,\hat{\alpha}_{HO}^{*}(\bm{\mathbf{L}}^{\prime})}-\expectationvalue{\hat{\alpha}_{HO}(\bm{\mathbf{L}})}\expectationvalue{\hat{\alpha}_{HO}^{*}(\bm{\mathbf{L}}^{\prime})}=(2\pi)^{2}\,\delta_{D}(\bm{\mathbf{L}}-\bm{\mathbf{L}}^{\prime})\,N_{HO}(\bm{\mathbf{L}}),$
(2.22)

where

$N_{HO}(\bm{\mathbf{L}})=N_{HO}(L)=\left(\sum_{IJ}[\bm{\mathbf{N}}^{-1}(L)]_{IJ}\right)^{-1}.$
(2.23)

We note that the $5\times 5$ symmetric matrix $\bm{\mathbf{N}}(L)$ is block
diagonal; it is comprised of the $1\times 1$ block for $BB$, the $2\times 2$
block for $TE$ and $EE$, and the $2\times 2$ block for $TB$ and $BB$. This
greatly simplifies the inversion of the matrix.

###  2.4 Global-minimum-variance estimator

Recently in the context of weak lensing reconstruction, Ref. [19] proposes
another quadratic estimator called the global-minimum-variance (GMV)
estimator, which is directly constructed from a minimum-variance linear
combination of pairs of all observables. This idea can be straightforwardly
applied to the case of cosmic birefringence, leading to the following
quadratic estimator

$\hat{\alpha}_{GMV}(\bm{\mathbf{L}})=\int_{\bm{\mathbf{L}}=\mathbf{l}_{1}-\mathbf{l}_{2}}X^{i}_{o}(\bm{\mathbf{l}}_{1})\,X^{j*}_{o}(\bm{\mathbf{l}}_{2})\,\Xi_{ij}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2}),$
(2.24)

where Einstein summation is assumed over $i$ and $j$, and $X^{i}=T$, $E$, or
$B$. In contrast, the HO estimator, discussed in the previous subsection and
historically widely used, is constructed as the minimum-variance linear
combination of estimators which are in turn constructed from pairs of
observables that correspond to different Fourier vectors. The work of Ref.
[19] shows that the reconstruction noise (variance) of the HO estimator is
strictly no less than that of the GMV estimator. Hence, the HO estimator is a
sub-optimal quadratic estimator. The two-step process to construct the HO
estimator compromises the optimality. In fact, the GMV estimator achieves the
least variance among all unbiased quadratic estimators.

The derivation needed to specify the GMV estimator closely parallels that in
the previous subsection. Here, we summarize the relevant results following the
notation of Ref. [19]. To this end, we introduce bold face symbols,
$\bm{\mathbf{f}}$, $\bm{\mathbf{C}}_{l}$, and $\bm{\mathbf{\Xi}}$, to denote
matrices whose matrix elements are given by $f_{ij}$, $C_{l}^{ij}$, and
$\Xi_{ij}$, respectively.

First of all, the GMV estimator can be made unbiased if we require

$1=\int_{\bm{\mathbf{L}}=\mathbf{l}_{1}-\mathbf{l}_{2}}f_{ij}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})\,\Xi_{ij}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2}).$
(2.25)

The variance of $\hat{\alpha}_{GMV}$ can be defined analogously to Eq. (2.22).
Under the condition of unbiasedness, the weights that minimize this variance
are

$\bm{\mathbf{\Xi}}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})=\frac{1}{2}\,N_{GMV}(L)\,\bm{\mathbf{C}}_{l_{1}}^{-1}\,\bm{\mathbf{f}}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})\,\bm{\mathbf{C}}_{l_{2}}^{-1},$
(2.26)

where the reconstruction noise spectrum $N_{GMV}(L)$ is given by the integral
of a matrix trace

$[N_{GMV}(L)]^{-1}=\frac{1}{2}\,\int_{\bm{\mathbf{L}}=\mathbf{l}_{1}-\mathbf{l}_{2}}{\rm
Tr}\left[\bm{\mathbf{C}}_{l_{1}}^{-1}\,\bm{\mathbf{f}}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})\,\bm{\mathbf{C}}_{l_{2}}^{-1}\,\bm{\mathbf{f}}(\bm{\mathbf{l}}_{2},\bm{\mathbf{l}}_{1})\right].$
(2.27)

Using the symmetry
$f_{YX}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})=f_{XY}(\bm{\mathbf{l}}_{2},\bm{\mathbf{l}}_{1})$,
we can rewrite
$\bm{\mathbf{f}}(\bm{\mathbf{l}}_{2},\bm{\mathbf{l}}_{1})=\bm{\mathbf{f}}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})^{T}$.

The precise gap between the GMV estimator and the HO estimator is more
manifest if we recast Eq. (2.20)) into the form

$\hat{\alpha}_{HO}(\bm{\mathbf{L}})=\sum_{i,j}\,\int_{\bm{\mathbf{L}}=\mathbf{l}_{1}-\mathbf{l}_{2}}X^{i}(\bm{\mathbf{l}}_{1})\,X^{j*}(\bm{\mathbf{l}}_{2})\,w_{ij}(|\bm{\mathbf{l}}_{1}-\bm{\mathbf{l}}_{2}|)\,F_{ij}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2}).$
(2.28)

A comparison with Eq. (2.24) suggests that $\hat{\alpha}_{HO}$ would be made
to match $\hat{\alpha}_{GMV}$ precisely if we allowed $w_{ij}$ to depend
arbitrarily on $\bm{\mathbf{l}}_{1}$ and $\bm{\mathbf{l}}_{2}$, i.e.
$w_{ij}=w_{ij}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})$, rather than on the
combination $L=|\bm{\mathbf{l}}_{1}-\bm{\mathbf{l}}_{2}|$. In that case, we
would have simply treated the product of
$w_{ij}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})$ and
$F_{ij}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})$ as a single coefficient
$\Xi_{ij}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})$, and hence followed the
logic of constructing the GMV estimator by minimizing the variance.

Nevertheless, the conceptual distinction between the GMV and HO estimators
does not lead to a drastic difference in the estimation error, as we
demonstrate both analytically and numerically in Section 6.

##  3 Maximum likelihood estimator

In this section, we sketch another approach to estimating birefringence
rotation in the CMB observables from the perspective of maximizing the
likelihood function. This approach is inspired by the work of Ref. [30] on
constructing the maximum likelihood estimator for weak gravitational lensing.
By modifying the maximum likelihood estimator, a quadratic estimator can be
derived, which we show is precisely the GMV estimator. The connection between
the maximum likelihood estimator and the quadratic estimator in terms of
Fisher matrices in the context of weak gravitational lensing has been noted in
Ref. [31].

Let $\tilde{\bm{\mathbf{d}}}$ denote a column vector that collectively
contains all the available CMB observables,
$\tilde{T}(\bm{\mathbf{l}}),\tilde{E}(\bm{\mathbf{l}}),\tilde{B}(\bm{\mathbf{l}})$,
for all wave vectors $\bm{\mathbf{l}}$, free of polarization rotation or
noise. The vector $\tilde{\bm{\mathbf{d}}}$ has a covariance matrix:

$\tilde{\bm{\mathbf{C}}}=\expectationvalue{\tilde{\bm{\mathbf{d}}}\,\tilde{\bm{\mathbf{d}}}^{T}},$
(3.1)

where random realizations of CMB anisotropies are averaged over. The
covariance matrix $\tilde{\bm{\mathbf{C}}}$ in principle follows from the
theoretical prediction that the CMB anisotropies of cosmological origin are
Gaussian random variables.

While the Fourier components $\tilde{T}(\bm{\mathbf{l}})$,
$\tilde{E}(\bm{\mathbf{l}})$ and $\tilde{B}(\bm{\mathbf{l}})$ are complex-
valued, in principle a set of real-valued variables can be instead used to
build the vector $\tilde{\bm{\mathbf{d}}}$. Correspondingly, the covariance
matrix is real symmetric and positive definite. For this reason, we shall have
in mind the machinery of linear algebra for real variables, and apply it to
the following derivations.

We note that there are no primary $B$ modes in the absence of primordial
tensor perturbations; [32, 33, 34] the cosmic $B$-mode polarization detected
in the CMB results from weak lensing by the large-scale structure (see Ref.
[35] for a review), and for this reason the CMB anisotropies acquire mode
couplings. Nevertheless, for our purpose it suffices to treat the $B$ modes in
the likelihood function as if they were genuinely Gaussian random primary
anisotropies, as long as the corresponding (isotropic) power spectrum is used
in constructing $\tilde{\bm{\mathbf{C}}}$.

The log likelihood of the vector $\tilde{\bm{\mathbf{d}}}$ given the
covariance matrix $\tilde{\bm{\mathbf{C}}}$ is

$\ln\mathcal{L}[\tilde{\bm{\mathbf{d}}}|\tilde{\bm{\mathbf{C}}}]=-\frac{1}{2}\,\tilde{\bm{\mathbf{d}}}^{T}\,\tilde{\bm{\mathbf{C}}}^{-1}\,\tilde{\bm{\mathbf{d}}}-\frac{1}{2}\,\ln\det\tilde{\bm{\mathbf{C}}}+\mathrm{const}.$
(3.2)

In the presence of rotation $\alpha$ and noise $\bm{\mathbf{n}}$, the measured
CMB observables can be written as

$\bm{\mathbf{d}}=\left[\bm{\mathbf{M}}(\alpha)\right]^{-1}\,\tilde{\bm{\mathbf{d}}}+\bm{\mathbf{n}},$
(3.3)

where $\bm{\mathbf{M}}(\alpha)$ is an orthogonal matrix (satisfying
$[\bm{\mathbf{M}}(\alpha)]^{-1}=[\bm{\mathbf{M}}(\alpha)]^{T}$) that describes
the action of pointwise polarization rotation on the anisotropies. Assuming
that the noise $\bm{\mathbf{n}}$ has a covariance matrix
$\bm{\mathbf{C}}_{n}$, the joint log likelihood of $\bm{\mathbf{d}}$ and
$\bm{\mathbf{n}}$ is now

$\ln\mathcal{L}[\bm{\mathbf{d}},\bm{\mathbf{n}}|\tilde{\bm{\mathbf{C}}},\alpha,\bm{\mathbf{C}}_{n}]=-\frac{1}{2}\,\bm{\mathbf{n}}^{T}\,\bm{\mathbf{C}}_{n}^{-1}\,\bm{\mathbf{n}}-\frac{1}{2}\,(\bm{\mathbf{d}}-\bm{\mathbf{n}})^{T}\,[\bm{\mathbf{M}}(\alpha)]^{T}\,\bm{\mathbf{\tilde{C}}}^{-1}\,\bm{\mathbf{M}}(\alpha)\,(\bm{\mathbf{d}}-\bm{\mathbf{n}}),$
(3.4)

where we have dropped irrelevant additive constant terms. For a fixed
$\bm{\mathbf{d}}$, Eq. (3.4) is maximized at

$\bm{\mathbf{n}}=\bm{\mathbf{n}}_{\mathrm{max}}=\left(\bm{\mathbf{C}}_{n}^{-1}+\bm{\mathbf{M}}^{T}\,\bm{\mathbf{\tilde{C}}}^{-1}\,\bm{\mathbf{M}}\right)^{-1}\,\bm{\mathbf{M}}^{T}\,\bm{\mathbf{\tilde{C}}}\,\bm{\mathbf{M}}\,\bm{\mathbf{d}},$
(3.5)

and the maximum value of the likelihood is

$\ln\mathcal{L}[\bm{\mathbf{d}}|\tilde{\bm{\mathbf{C}}},\alpha,\bm{\mathbf{n}}=\bm{\mathbf{n}}_{\mathrm{max}}]=-\frac{1}{2}\,\bm{\mathbf{d}}^{T}\,(\bm{\mathbf{C}}_{n}+\bm{\mathbf{M}}^{T}\,\tilde{\bm{\mathbf{C}}}\,\bm{\mathbf{M}})^{-1}\,\bm{\mathbf{d}},$
(3.6)

which we shall simply call $\ln\mathcal{L}$. We can now estimate the rotation
field by maximising this log likelihood with respect to $\alpha$. Anticipating
$\alpha$ to be numerically small in radians [17], we consider an approximation
of $\ln\mathcal{L}$ as a second-order polynomial in $\alpha$:

$\ln\mathcal{L}=\mathrm{const}+\bm{\mathbf{d}}^{T}\,\bm{\mathbf{K}}_{1}(\alpha)\,\bm{\mathbf{d}}-\frac{1}{2}\,\bm{\mathbf{d}}^{T}\,\bm{\mathbf{K}}_{2}(\alpha)\,\bm{\mathbf{d}}+O(\alpha^{3}).$
(3.7)

The matrices $\bm{\mathbf{K}}_{1}(\alpha)$ and $\bm{\mathbf{K}}_{2}(\alpha)$
have linear and quadratic dependence on $\alpha$, respectively. To find these
matrices, we first Taylor-expand $\bm{\mathbf{M}}(\alpha)$:

$\bm{\mathbf{M}}(\alpha)=\bm{\mathbf{I}}-\bm{\mathbf{T}}(\alpha)+O(\alpha^{2}),$
(3.8)

where $\bm{\mathbf{T}}(\alpha)$ is an antisymmetric matrix with linear
dependence on $\alpha$. Substituting this expansion into Eq. (3.6) and
defining $\bm{\mathbf{C}}=\tilde{\bm{\mathbf{C}}}+\bm{\mathbf{C}}_{n}$, we
obtain

$\displaystyle\bm{\mathbf{K}}_{1}$
$\displaystyle=\bm{\mathbf{C}}^{-1}\,\bm{\mathbf{T}}\,\tilde{\bm{\mathbf{C}}}\,\bm{\mathbf{C}}^{-1},$
(3.9) $\displaystyle\bm{\mathbf{K}}_{2}$
$\displaystyle=\bm{\mathbf{C}}^{-1}\,\bm{\mathbf{T}}^{T}\,\tilde{\bm{\mathbf{C}}}\,\bm{\mathbf{C}}^{-1}\,\tilde{\bm{\mathbf{C}}}\,\bm{\mathbf{T}}\,\bm{\mathbf{C}}^{-1}+\bm{\mathbf{C}}^{-1}\,\tilde{\bm{\mathbf{C}}}\,\bm{\mathbf{T}}^{T}\,\bm{\mathbf{C}}^{-1}\,\bm{\mathbf{T}}\,\tilde{\bm{\mathbf{C}}}\,\bm{\mathbf{C}}^{-1}-2\,\bm{\mathbf{C}}^{-1}\,\bm{\mathbf{T}}^{T}\,\tilde{\bm{\mathbf{C}}}\,\bm{\mathbf{C}}^{-1}\,\bm{\mathbf{T}}\,\tilde{\bm{\mathbf{C}}}\,\bm{\mathbf{C}}^{-1}$
$\displaystyle\quad+\bm{\mathbf{C}}^{-1}\,\bm{\mathbf{T}}^{T}\,\bm{\mathbf{T}}\,\tilde{\bm{\mathbf{C}}}\,\bm{\mathbf{C}}^{-1}-\bm{\mathbf{C}}^{-1}\,\bm{\mathbf{T}}^{T}\,\tilde{\bm{\mathbf{C}}}\,\bm{\mathbf{T}}\,\bm{\mathbf{C}}^{-1}.$
(3.10)

We can make the dependence on $\alpha$ explicit by writing

$\displaystyle\bm{\mathbf{T}}(\alpha)$
$\displaystyle=\sum_{\bm{\mathbf{L}}}\,\bm{\mathbf{T}}_{\bm{\mathbf{L}}}\,\alpha_{\bm{\mathbf{L}}},$
(3.11)
$\displaystyle\bm{\mathbf{d}}^{T}\,\bm{\mathbf{K}}_{1}(\alpha)\,\bm{\mathbf{d}}$
$\displaystyle=\sum_{\bm{\mathbf{L}}}\,\mathcal{K}_{1,\bm{\mathbf{L}}}\,\alpha_{\bm{\mathbf{L}}},$
(3.12)
$\displaystyle\bm{\mathbf{d}}^{T}\,\bm{\mathbf{K}}_{2}(\alpha)\,\bm{\mathbf{d}}$
$\displaystyle=\sum_{\bm{\mathbf{L}}}\,\sum_{\bm{\mathbf{L}}^{\prime}}\,\mathcal{K}_{2,\bm{\mathbf{L}}\bm{\mathbf{L}}^{\prime}}\,\alpha_{\bm{\mathbf{L}}}^{*}\,\alpha_{\bm{\mathbf{L}}^{\prime}},$
(3.13)

where the indices $\bm{\mathbf{L}}$ and $\bm{\mathbf{L}}^{\prime}$ label the
components of the anisotropic rotation field on the sky (i.e. its various
Fourier modes $\alpha(\bm{\mathbf{L}})$). With the cubic and higher order
terms neglected, Eq. (3.6) takes a quadratic form:

$\ln\mathcal{L}=\mathrm{const}+\sum_{\bm{\mathbf{L}}}\,\mathcal{K}_{1,\bm{\mathbf{L}}}\,\alpha_{\bm{\mathbf{L}}}-\frac{1}{2}\,\sum_{\bm{\mathbf{L}}}\,\sum_{\bm{\mathbf{L}}^{\prime}}\,\mathcal{K}_{2,\bm{\mathbf{L}}\bm{\mathbf{L}}^{\prime}}\,\alpha_{\bm{\mathbf{L}}}^{*}\,\alpha_{\bm{\mathbf{L}}^{\prime}}.$
(3.14)

This log likelihood is maximized for

$\alpha_{\bm{\mathbf{L}}}=\sum_{\bm{\mathbf{L}}^{\prime}}\,[\mathcal{K}_{2}^{-1}]_{\bm{\mathbf{L}}\bm{\mathbf{L}}^{\prime}}\,\mathcal{K}_{1,\bm{\mathbf{L}}^{\prime}},$
(3.15)

where the matrix inverse of $\mathcal{K}_{2}$ is taken with respect to the
indices $\bm{\mathbf{L}}$ and $\bm{\mathbf{L}}^{\prime}$. Note that the
matrices $\mathcal{K}_{1,\bm{\mathbf{L}}}$ and
$\mathcal{K}_{2,\bm{\mathbf{L}}\bm{\mathbf{L}}^{\prime}}$ both depend
quadratically on the specific realization of the measured CMB observables
$\bm{\mathbf{d}}$.

##  4 Relation between various quadratic estimators

###  4.1 Equivalence to GMV estimator

To make contact with quadratic estimators, we consider replacing the matrix
$\mathcal{K}_{2}$ with
$\mathrm{bar}{\mathcal{K}}_{2}=\expectationvalue{\mathcal{K}_{2}}$, which is
the ensemble expectation after averaging over random realizations of
$\bm{\mathbf{d}}$. The matrix $\mathcal{K}_{2}$ that correspond to any
specific realization of $\bm{\mathbf{d}}$ fluctuates around
$\mathrm{bar}{\mathcal{K}}_{2}$, but the fluctuation is usually small in the
limit of a large number of Fourier modes. Without having any statistically
preferred direction on the sky, the matrix $\mathrm{bar}{\mathcal{K}}_{2}$ is
shown below to be diagonal in $\bm{\mathbf{L}}$ space. We therefore define the
reconstructed rotation field as follows, by modifying Eq. (3.15):

$\alpha_{\mathrm{rec},{\bm{\mathbf{L}}}}=\sum_{\bm{\mathbf{L}}^{\prime}}\,[\mathrm{bar}{\mathcal{K}}_{2}^{-1}]_{\bm{\mathbf{L}}\bm{\mathbf{L}}^{\prime}}\,\mathcal{K}_{1,\bm{\mathbf{L}}^{\prime}}.$
(4.1)

Since the only dependence on the data $\bm{\mathbf{d}}$ is in
$\mathcal{K}_{1,\bm{\mathbf{L}}^{\prime}}$ as a bilinear form in
$\bm{\mathbf{d}}$, $\alpha_{\mathrm{rec},\bm{\mathbf{L}}}$ is a quadratic
estimator. After some algebra and using Eq. (3.1), we derive

$\mathrm{bar}{\mathcal{K}}_{2,\bm{\mathbf{L}}\bm{\mathbf{L}}^{\prime}}={\rm
Tr}\left[\tilde{\bm{\mathbf{C}}}\,\bm{\mathbf{C}}^{-1}\,\tilde{\bm{\mathbf{C}}}\,\bm{\mathbf{T}}^{T}_{-\bm{\mathbf{L}}}\,\bm{\mathbf{C}}^{-1}\,\bm{\mathbf{T}}_{\bm{\mathbf{L}}^{\prime}}\right]-{\rm
Tr}\left[\tilde{\bm{\mathbf{C}}}\,\bm{\mathbf{C}}^{-1}\,\bm{\mathbf{T}}^{T}_{-\bm{\mathbf{L}}}\,\tilde{\bm{\mathbf{C}}}\,\bm{\mathbf{C}}^{-1}\,\bm{\mathbf{T}}_{\bm{\mathbf{L}}^{\prime}}\right].$
(4.2)

Finally, we show that $\alpha_{\mathrm{rec},{\bm{\mathbf{L}}}}$ coincides with
the GMV estimator $\hat{\alpha}_{GMV}(L)$ constructed in Section 2.4. To this
end, we first recast some of the equations in the preceding sections in terms
of matrices. Eq. (2.1) becomes Eq. (3.1), while Eq. (2.2) is replaced by

$\expectationvalue{[\bm{\mathbf{M}}(\alpha)]^{T}\,\tilde{\bm{\mathbf{d}}}\,\tilde{\bm{\mathbf{d}}}^{T}\,\bm{\mathbf{M}}(\alpha)}=\expectationvalue{\tilde{\bm{\mathbf{d}}}\,\tilde{\bm{\mathbf{d}}}^{T}}+\sum_{\bm{\mathbf{L}}}\,\bm{\mathbf{f}}_{\bm{\mathbf{L}}}\,\alpha_{\bm{\mathbf{L}}}+O(\alpha^{2}).$
(4.3)

In other words, an arbitrary matrix element of
$\bm{\mathbf{f}}_{\bm{\mathbf{L}}}$ has the form
$f_{ij}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})\,(2\pi)^{2}\,\delta_{D}(\bm{\mathbf{l}}_{1}-\bm{\mathbf{l}}_{2}-\bm{\mathbf{L}})$.
Using Eq. (3.8) and Eq. (3.11), we find

$\bm{\mathbf{f}}_{\bm{\mathbf{L}}}=\bm{\mathbf{T}}_{\bm{\mathbf{L}}}\,\tilde{\bm{\mathbf{C}}}+\tilde{\bm{\mathbf{C}}}\,\bm{\mathbf{T}}^{T}_{\bm{\mathbf{L}}}.$
(4.4)

On the other hand, an arbitrary element of $\bm{\mathbf{C}}$, denoted
$C^{i\bm{\mathbf{l}}_{1},j\bm{\mathbf{l}}_{2}}$, is equal to
$C^{ij}_{l_{1}}\,(2\pi)^{2}\,\delta_{D}(\bm{\mathbf{l}}_{1}-\bm{\mathbf{l}}_{2})$.
With these in mind, we can rewrite Eq. (2.27) as

$\displaystyle(2\pi)^{2}\delta_{D}(\bm{\mathbf{L}}-\bm{\mathbf{L}}^{\prime})N_{GMV}(L)^{-1}$
$\displaystyle=\frac{1}{2}\,{\rm
Tr}\left[\bm{\mathbf{C}}^{-1}\,\bm{\mathbf{f}}_{-\bm{\mathbf{L}}}\,\bm{\mathbf{C}}^{-1}\,\bm{\mathbf{f}}_{\bm{\mathbf{L}}^{\prime}}\right]$
$\displaystyle=\frac{1}{2}\,{\rm
Tr}\left[\bm{\mathbf{C}}^{-1}\,\bm{\mathbf{T}}_{-\bm{\mathbf{L}}}\,\tilde{\bm{\mathbf{C}}}\,\bm{\mathbf{C}}^{-1}\,\bm{\mathbf{T}}_{\bm{\mathbf{L}}^{\prime}}\,\tilde{\bm{\mathbf{C}}}+\bm{\mathbf{C}}^{-1}\,\tilde{\bm{\mathbf{C}}}\,\bm{\mathbf{T}}^{T}_{-\bm{\mathbf{L}}}\,\bm{\mathbf{C}}^{-1}\,\bm{\mathbf{T}}_{\bm{\mathbf{L}}^{\prime}}\,\tilde{\bm{\mathbf{C}}}\right.$
$\displaystyle\quad+\left.\bm{\mathbf{C}}^{-1}\,\tilde{\bm{\mathbf{C}}}\,\bm{\mathbf{T}}^{T}_{-\bm{\mathbf{L}}}\,\bm{\mathbf{C}}^{-1}\,\bm{\mathbf{T}}_{\bm{\mathbf{L}}^{\prime}}\,\tilde{\bm{\mathbf{C}}}+\bm{\mathbf{C}}^{-1}\,\tilde{\bm{\mathbf{C}}}\,\bm{\mathbf{T}}^{T}_{-\bm{\mathbf{L}}}\,\bm{\mathbf{C}}^{-1}\,\tilde{\bm{\mathbf{C}}}\,\bm{\mathbf{T}}^{T}_{\bm{\mathbf{L}}^{\prime}}\right]$
$\displaystyle=\mathrm{bar}{\mathcal{K}}_{2,\bm{\mathbf{L}}\bm{\mathbf{L}}^{\prime}},$
(4.5)

where in the last line we have used the cyclic property of the trace and the
antisymmetry of $\bm{\mathbf{T}}_{L}$ to match Eq. (4.2). We note that
$\mathrm{bar}{\mathcal{K}}_{2,\bm{\mathbf{L}}\bm{\mathbf{L}}^{\prime}}$ is
indeed diagonal in $\bm{\mathbf{L}}$ space, as expected from statistical
isotropy. Finally, Eq. (2.24) and Eq. (2.26) can be rewritten as

$\displaystyle\hat{\alpha}_{GMV}(\bm{\mathbf{L}})$
$\displaystyle=\frac{1}{2}\,N_{GMV}(L)^{-1}\,\bm{\mathbf{d}}^{T}\,\bm{\mathbf{C}}^{-1}\,\bm{\mathbf{f}}_{\bm{\mathbf{L}}}\,\bm{\mathbf{C}}^{-1}\,\bm{\mathbf{d}}$
(4.6)
$\displaystyle=\frac{1}{2}\,\sum_{\bm{\mathbf{L}}^{\prime}}\,\mathrm{bar}{\mathcal{K}}_{2,\bm{\mathbf{L}}\bm{\mathbf{L}}^{\prime}}\,\bm{\mathbf{d}}^{T}\,\bm{\mathbf{C}}^{-1}\,\left(\bm{\mathbf{T}}_{\bm{\mathbf{L}}^{\prime}}\,\tilde{\bm{\mathbf{C}}}+\tilde{\bm{\mathbf{C}}}\,\bm{\mathbf{T}}^{T}_{\bm{\mathbf{L}}^{\prime}}\right)\,\bm{\mathbf{C}}^{-1}\,\bm{\mathbf{d}}$
(4.7)
$\displaystyle=\sum_{\bm{\mathbf{L}}^{\prime}}\,\mathrm{bar}{\mathcal{K}}_{2,\bm{\mathbf{L}}\bm{\mathbf{L}}^{\prime}}\,\bm{\mathbf{d}}^{T}\,\bm{\mathbf{C}}^{-1}\,\bm{\mathbf{T}}_{\bm{\mathbf{L}}^{\prime}}\,\tilde{\bm{\mathbf{C}}}\,\bm{\mathbf{C}}^{-1}\,\bm{\mathbf{d}}$
(4.8)
$\displaystyle=\sum_{\bm{\mathbf{L}}^{\prime}}\,\mathrm{bar}{\mathcal{K}}_{2,\bm{\mathbf{L}}\bm{\mathbf{L}}^{\prime}}\,\mathcal{K}_{1\bm{\mathbf{L}}^{\prime}}=\hat{\alpha}_{\mathrm{rec},\bm{\mathbf{L}}}.$
(4.9)

We note that the above derivation of the connection between the maximum
likelihood estimator and GMV estimator assumes an antisymmetric
$\bm{\mathbf{T}}$ matrix (or more generally an anti-Hermitian matrix). In more
familiar notation, if

$\begin{pmatrix}\delta T(\bm{\mathbf{l}})\\\ \delta E(\bm{\mathbf{l}})\\\
\delta
B(\bm{\mathbf{l}})\end{pmatrix}=\int\frac{\mathrm{d}^{2}\bm{\mathbf{l}}^{\prime}}{(2\pi)^{2}}\bm{\mathbf{W}}(\bm{\mathbf{l}},\bm{\mathbf{l}}^{\prime})\begin{pmatrix}\tilde{T}(\bm{\mathbf{l}}^{\prime})\\\
\tilde{E}(\bm{\mathbf{l}}^{\prime})\\\
\tilde{B}(\bm{\mathbf{l}}^{\prime})\end{pmatrix},$ (4.10)

for some $3\times 3$ matrix $\bm{\mathbf{W}}$ that depends on
$\bm{\mathbf{l}}$ and $\bm{\mathbf{l}}^{\prime}$, then we require
$\bm{\mathbf{W}}^{T}(\bm{\mathbf{l}}^{\prime},\bm{\mathbf{l}})=-\bm{\mathbf{W}}(\bm{\mathbf{l}},\bm{\mathbf{l}}^{\prime})$
[or more generally
$\bm{\mathbf{W}}^{T}(\bm{\mathbf{l}}^{\prime},\bm{\mathbf{l}})=-\bm{\mathbf{W}}^{*}(\bm{\mathbf{l}},\bm{\mathbf{l}}^{\prime})$].
While this is true for polarization rotations, it is not true for
gravitational lensing (see Eq. (4) of [18]).

###  4.2 Equivalence between HO and GMV estimators

Having summarized the quadratic estimators, we now quantify their performance
using the power spectrum of reconstruction noise, namely the variance of
$\hat{\alpha}(\bm{\mathbf{L}})$ as a function of $L=|\bm{\mathbf{L}}|$. We
will perform calculations under the null hypothesis as compelling evidence for
anisotropic CMB birefringence has not yet been found. Since we have already
shown in Section 3 that the simple modification of the maximum likelihood
estimator is equivalent to the GMV estimator, it remains to compare the GMV
and HO estimators. We show that, if the two are constructed using only the
polarization but not the temperature anisotropies, then they have the same
reconstruction noise power spectrum. When the temperature anisotropies are
included in estimating rotation, we find numerically in Section 6 that the
difference in the reconstruction noise spectra between the GMV and HO
estimators is insignificant.

From Eq. (2.14) and Eq. (2.16), we can see that the presence of a non-
vanishing cross-variable spectrum, $C^{TE}_{l}$, makes formulas more
complicated (note that $C^{EB}_{l}$ vanishes due to statistical invariance
under parity). This complication is avoided if temperature anisotropies
$T(\bm{\mathbf{l}})$ are entirely excluded from the construction of the GMV
and HO estimators. When only the $E$ and $B$ polarization anisotropies are
considered, it is possible to show that the GMV and HO estimators are
equivalent. In fact, this is a special case of the general statement that the
GMV and HO estimators are identical whenever $C^{XY}_{l}=0$ for any $X\neq Y$,
which we now show.

Assuming $C^{XY}_{l}=0$ for $X\neq Y$, Eq. (2.18) implies that
$\bm{\mathbf{N}}$ is a diagonal matrix, _i.e._ there is no covariance between
different single-pair estimators. Using Eq. (2.9), Eq. (2.13), and Eq. (2.20),
the HO estimator becomes

$\hat{\alpha}_{HO}(\bm{\mathbf{L}})=N_{HO}(L)\sum_{(X,Y)}\int_{\bm{\mathbf{L}}=\bm{\mathbf{l}}_{1}-\bm{\mathbf{l}}_{2}}X_{o}(\bm{\mathbf{l}}_{1})\,Y^{*}_{o}(\bm{\mathbf{l}}_{2})\,\frac{f_{XY}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})}{(1+\delta_{XY})\,C^{XX}_{l_{1}}\,C^{YY}_{l_{2}}},$
(4.11)

where the sum is over unique pairs of observables $XY$, and

$N_{HO}(L)=\left(\sum_{(X,Y)}N^{-1}_{XY}(L)\right)^{-1}.$ (4.12)

On the other hand, Eq. (2.24) can be rewritten as

$\hat{\alpha}_{GMV}=N_{GMV}(L)\sum_{(X,Y)}\int_{\bm{\mathbf{L}}=\bm{\mathbf{l}}_{1}-\bm{\mathbf{l}}_{2}}X_{o}(\bm{\mathbf{l}}_{1})\,Y^{*}_{o}(\bm{\mathbf{l}}_{2})\,\frac{f_{XY}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})}{(1+\delta_{XY})\,C^{XX}_{l_{1}}\,C^{YY}_{l_{2}}}.$
(4.13)

The equality of the reconstruction noise spectra $N_{HO}(L)=N_{GMV}(L)$ can be
straightforwardly derived from Eq. (2.27) and the symmetry property of
$\bm{\mathbf{f}}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})$. Therefore, we
conclude that $\hat{\alpha}_{HO}(L)=\hat{\alpha}_{GMV}(L)$ if $C^{XY}_{l}=0$
for $X\neq Y$.

In the special case where only polarization anisotropies $E$ and $B$ are
considered, the reconstruction noise spectrum has the form

$N_{HO}(L)=N_{GMV}(L)=\left(\int_{\bm{\mathbf{L}}=\mathbf{l}_{1}-\mathbf{l}_{2}}\left[\frac{f_{EE}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})^{2}}{2\,C_{l_{1}}^{EE}C_{l_{2}}^{EE}}+\frac{f_{EB}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})^{2}}{C_{l_{1}}^{EE}C_{l_{2}}^{BB}}+\frac{f_{BB}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})^{2}}{2\,C_{l_{1}}^{BB}C_{l_{2}}^{BB}}\right]\right)^{-1}.$
(4.14)

##  5 Simultaneous estimation of lensing and rotation

We have carried out the analysis under the assumption that the power spectra
of lensed (but unrotated) CMB observables are known. The estimator for the
rotation $\hat{\alpha}$ is constructed from these lensed power spectra through
the weights $F_{XY}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})$ or
$\Xi_{ij}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})$. However, one might worry
that the estimate for lensing is contaminated by the presence of anisotropic
rotation, and conversely the estimate for rotation might be contaminated by
lensing. We show in this section that any dependence of quadratic estimators
$\expectationvalue{\hat{\alpha}_{XY}}$ on the lensing potential $\phi$ is of
$\mathcal{O}(\phi^{2})$, and likewise any dependence of quadratic estimators
for the lensing potential $\expectationvalue{\hat{\phi}_{XY}}$ on the
polarization rotation is of order $\mathcal{O}(\alpha^{2})$. In other words,
at linear order $\alpha$ and $\phi$ can be simultaneously estimated from the
CMB anisotropies without bias [14]. However, the second-order contamination of
the rotation power spectrum due to lensing will eventually need to be
corrected for in future experiments (see Fig. 2 in Ref. [36]).

We first rewrite Eq. (2.2) to include the linear response due to the lensing
potential:

$\expectationvalue{X(\bm{\mathbf{l}}_{1})Y^{*}(\bm{\mathbf{l}}_{2})}=\expectationvalue{\mathrm{bar}X(\bm{\mathbf{l}}_{1})\mathrm{bar}Y^{*}(\bm{\mathbf{l}}_{2})}+f_{XY}^{\phi}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})\phi(\bm{\mathbf{l}}_{1}-\bm{\mathbf{l}}_{2})+f_{XY}^{\alpha}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})\alpha(\bm{\mathbf{l}}_{1}-\bm{\mathbf{l}}_{2})+O(\alpha^{2},\phi^{2},\alpha\phi),$
(5.1)

where $\mathrm{bar}X,\mathrm{bar}Y$ denote the unlensed and unrotated
anisotropies. We now write down the “naive” minimum-variance estimator for
$\alpha$ and $\phi$, respectively, from the observed quantities (with
subscript $o$), pretending that only one effect takes place at a time:

$\displaystyle\hat{\alpha}_{XY}$
$\displaystyle=\int_{\bm{\mathbf{L}}=\bm{\mathbf{l}}_{1}-\bm{\mathbf{l}}_{2}}X_{o}(\bm{\mathbf{l}}_{1})Y_{o}^{*}(\bm{\mathbf{l}}_{2})F^{\alpha}_{XY}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2}),$
(5.2) $\displaystyle\hat{\phi}_{XY}$
$\displaystyle=\int_{\bm{\mathbf{L}}=\bm{\mathbf{l}}_{1}-\bm{\mathbf{l}}_{2}}X_{o}(\bm{\mathbf{l}}_{1})Y_{o}^{*}(\bm{\mathbf{l}}_{2})F^{\phi}_{XY}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2}).$
(5.3)

Here the weights $F^{\alpha}_{XY}$ are the same as in Eq. (2.13) and Eq.
(2.14). The weights $F^{\phi}_{XY}$ are defined in Ref. [19].

For $\bm{\mathbf{L}}\neq 0$, the rotation estimator has an expectation value

$\displaystyle\expectationvalue{\hat{\alpha}_{XY}(\bm{\mathbf{L}})}$
$\displaystyle=\int_{\bm{\mathbf{L}}=\bm{\mathbf{l}}_{1}-\bm{\mathbf{l}}_{2}}\expectationvalue{\mathrm{bar}X(\bm{\mathbf{l}}_{1})\mathrm{bar}Y^{*}(\bm{\mathbf{l}}_{2})}\,F^{\alpha}_{XY}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})+\int_{\bm{\mathbf{L}}=\bm{\mathbf{l}}_{1}-\bm{\mathbf{l}}_{2}}\expectationvalue{\mathrm{bar}X_{n}(\bm{\mathbf{l}}_{1})\mathrm{bar}Y_{n}^{*}(\bm{\mathbf{l}}_{2})}\,F^{\alpha}_{XY}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})$
$\displaystyle\quad+\alpha(\bm{\mathbf{L}})\int_{\bm{\mathbf{L}}=\bm{\mathbf{l}}_{1}-\bm{\mathbf{l}}_{2}}f^{\alpha}_{XY}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})\,F^{\alpha}_{XY}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})+\phi(\bm{\mathbf{L}})\int_{\bm{\mathbf{L}}=\bm{\mathbf{l}}_{1}-\bm{\mathbf{l}}_{2}}f^{\phi}_{XY}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})\,F^{\alpha}_{XY}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2}),$
(5.4)

where $f^{\alpha}_{XY}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})$ are found in
Table 1, and $f^{\phi}_{XY}$ are found in Table 1 of Ref. [19]. The first two
terms vanish because $\bm{\mathbf{l}}_{1}\neq\bm{\mathbf{l}}_{2}$. The third
term reduces to $\alpha(\bm{\mathbf{L}})$ due to the unbiasedness of the
estimator (Eq. (2.11)). It remains to see whether the last term vanishes, or
in other words, whether the naive estimator for $\alpha$ is contaminated by
lensing.

Here, we show explicitly the proof for $XY=EB$, but the result holds
analogously for all pairs $XY$.

$\displaystyle\quad\int_{\bm{\mathbf{L}}=\bm{\mathbf{l}}_{1}-\bm{\mathbf{l}}_{2}}f^{\phi}_{EB}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})\,F^{\alpha}_{EB}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})$
(5.5)
$\displaystyle=\int_{\bm{\mathbf{L}}=\bm{\mathbf{l}}_{1}-\bm{\mathbf{l}}_{2}}\frac{f^{\phi}_{EB}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})\,f^{\alpha}_{XY}(\bm{\mathbf{l}}_{1},\bm{\mathbf{l}}_{2})}{C^{EB}_{l_{1}}C^{EB}_{l_{2}}}$
(5.6)
$\displaystyle=\int\frac{\mathrm{d}^{2}\bm{\mathbf{l}}_{1}\,\mathrm{d}^{2}\bm{\mathbf{l}}_{2}}{(2\pi)^{2}}\delta_{D}(l_{1x}-l_{2x}-L_{x})\,\delta_{D}(l_{1y}-l_{2y}-L_{y})\,\left(C^{EB}_{l_{1}}\,C^{EB}_{l_{2}}\right)^{-1}$
(5.7) $\displaystyle\quad\times\left(\tilde{C}^{E\nabla
E}_{l_{1}}\,\bm{\mathbf{L}}\cdot\bm{\mathbf{l}}_{1}+\tilde{C}^{B\nabla
B}_{l_{2}}\,\bm{\mathbf{L}}\cdot\bm{\mathbf{l}}_{2}\right)\,\sin
2\varphi_{12}\times
2\,\left(\tilde{C}^{BB}_{l_{1}}-\tilde{C}^{EE}_{l_{2}}\right)\,\cos
2\varphi_{12}.$ (5.8)

The quantities $\tilde{C}^{X\nabla Y}_{l}$ are the lensed gradient spectrum,
defined in [37, 38]. Under a change of variables: $l_{iy}\to-l_{iy}$ for
$i=1,2$, we have $\varphi_{12}\to-\varphi_{12}$ and

$\displaystyle\delta_{D}(l_{1y}-l_{2y}-L_{y})\,\bm{\mathbf{L}}\cdot\bm{\mathbf{l}}_{i}$
$\displaystyle\to\delta_{D}(-l_{1y}+l_{2y}-L_{y})\,(L_{x}l_{ix}-L_{y}l_{iy})$
(5.9)
$\displaystyle=\delta_{D}(-l_{1y}+l_{2y}-L_{y})\,[L_{x}l_{ix}+(l_{1y}-l_{2y})\,l_{iy}]$
(5.10)
$\displaystyle=\delta_{D}(l_{1y}-l_{2y}-L_{y})\,[L_{x}l_{ix}+(l_{1y}-l_{2y})\,l_{iy}]$
(5.11)
$\displaystyle=\delta_{D}(l_{1y}-l_{2y}-L_{y})\,\bm{\mathbf{L}}\cdot\bm{\mathbf{l}}_{i}.$
(5.12)

The whole integrand is odd under this change of variables, and hence the
integral vanishes.

By putting together Table 1, and Table 1 of Ref. [19], it is easy to see that
the corresponding integrals for other pairs $XY$ also vanish. We conclude that

$\expectationvalue{\hat{\alpha}_{XY}(\bm{\mathbf{L}})}=\alpha(\bm{\mathbf{L}})+O(\alpha^{2},\phi^{2},\alpha\phi),$
(5.13)

where $\hat{\alpha}_{XY}$ here is constructed from unlensed and unrotated
spectra. A similar statement can be made about
$\hat{\phi}_{XY}(\bm{\mathbf{L}})$ by the same parity argument. Ref. [14]
reached the same conclusion based on the parity argument, for the full-sky
version of these estimators.

##  6 Numerical forecasts

We compute the reconstruction noise spectra for the HO and GMV estimators
using polarization and temperature anisotropies. The lensed CMB power spectra
($\tilde{C}^{XY}_{l}$) are generated by CAMB [39] up to
$l_{\mathrm{max}}=3000$ [40], with concordance cosmological parameters
$H_{0}=67.4\,\mathrm{km\,s^{-1}\,s^{-1}\,Mpc}$, $\Omega_{b}h^{2}=0.0224$,
$\Omega_{c}h^{2}=0.120$, $m_{\nu}=0.06\,\mathrm{eV}$, $\Omega_{K}=0$, and
$\tau=0.054$.

We assume that the noise is Gaussian, with an isotropic power spectrum in the
form [18]

$C^{XY,n}_{l}=\delta_{XY}\,\Delta_{X}^{2}\,e^{(l^{2}\,\Theta^{2})/(8\ln 2)},$
(6.1)

where $\Theta$ is the full width at half maximum (FWHM) of the Gaussian beam
profile, and $\Delta_{X}=\Delta_{T}$ and $\Delta_{P}$ for temperature ($X=T$)
and polarization ($X=E,\,B$), respectively. We assume a fully polarized
detector, with $\sqrt{2}\,\Delta_{T}=\Delta_{P}$ [18].

We perform calculations for three sets of detector parameters $\Delta_{T}$ and
$\Theta$, which are appropriate for a range of CMB experiments: Planck 2018
[41], Simons Observatory [42], and a Stage-IV CMB experiment (CMB-S4) [43].
See Ref. [44] for details on the the detection of birefringence using the
latter two experiments. In all three experiments, we include primary
anisotropy modes in the range $l_{\rm min}\leqslant l\leqslant l_{\rm max}$.
Refer to Table 2 for the parameter values we use.

| $\Delta T\,[\mathrm{\mu K\,arcmin}]$ | $\Theta\,[\mathrm{arcmin}]$ | $l_{\mathrm{min}}$ | $l_{\mathrm{max}}$ | $f_{\mathrm{sky}}$
---|---|---|---|---|---
Planck SMICA | 45.0 | 5.0 | 5 | 3000 | 0.7
Simons Observatory | 7.0 | 1.4 | 30 | 3000 | 0.4
CMB-S4 | 1.0 | 1.4 | 30 | 3000 | 0.4
Table 2: Benchmark experimental parameters considered in this work [45].

$\displaystyle{10^{-3}}$$\displaystyle{10^{0}}$$\displaystyle{10^{3}}$$\displaystyle{10^{6}}$$\displaystyle[L(L+1)N(L)/2\pi]^{1/2}$(deg)Hu-
OkamotoEB$\displaystyle{10^{0}}$$\displaystyle{10^{1}}$$\displaystyle{10^{2}}$$\displaystyle{10^{3}}$$\displaystyle
L$$\displaystyle{10^{-3}}$$\displaystyle{10^{0}}$$\displaystyle{10^{3}}$$\displaystyle{10^{6}}$$\displaystyle[L(L+1)N(L)/2\pi]^{1/2}$(deg)TB$\displaystyle{10^{0}}$$\displaystyle{10^{1}}$$\displaystyle{10^{2}}$$\displaystyle{10^{3}}$$\displaystyle
L$TEPlanck SMICASimonsCMB-S4

Figure 1: Reconstruction noise spectra for the HO estimator computed for the
three sets of experimental parameters considered in this work. We plot the
reconstruction noise spectra for the combined HO estimator, as well as those
for using only one pair of CMB observables at a time with the lowest
reconstruction noise. That is, we selectively present the results for three
such pairs that make the greatest contribution to the overall detectability of
the Hu-Okamoto estimator.

Mathematically, the reconstruction noise of the GMV estimator $N_{GMV}(L)$ is
strictly smaller than or equal to that for the HO estimator $N_{HO}(L)$.
Numerically, we find that the fractional difference between the two is no
greater than 0.5% for $L=1-10^{3}$. For this reason, we have not included
separate curves for $N_{GMV}(L)$ in Figure 1, as the difference is too small
to be visually discernible (they would overlap with the curves in the top-left
plot for the Hu-Okamoto estimator).

The goal of estimating spatially dependent polarization rotation is to detect
any new physics that might possibly cause this phenomenon, such as a
pseudoscalar axion-like field that couples to the electromagnetic field. In
this section, we study the constraining power of the quadratic estimators on
the scenario that the rotation results from a network of cosmic strings. As
explained in Ref. [8], these are topological remnants reminiscent of a phase
transition in the pseudoscalar field in the primordial Universe.

Since the cosmic strings that comprise the string network have random
configurations, a key theoretical prediction is the power spectrum for the
rotation field $\alpha$:

$\expectationvalue{\alpha(\bm{\mathbf{L}})\alpha(\bm{\mathbf{L}}^{\prime})}=(2\pi)^{2}\,\delta_{D}(\bm{\mathbf{L}}-\bm{\mathbf{L}}^{\prime})\,K\,P(L),$
(6.2)

where the notation for ensemble averaging $\langle\cdots\rangle$ is over all
string network realizations. Here $P(L)$ is a fiducial functional form for the
power spectrum, and $K$ is a normalization constant that we would like to
detect or constrain. In the presence of strings, a general quadratic estimator
$\hat{\alpha}(\bm{\mathbf{L}})$ has a covariance

$\expectationvalue{\hat{\alpha}(\bm{\mathbf{L}})\,\hat{\alpha}(\bm{\mathbf{L}}^{\prime})}=(2\pi)^{2}\,\delta_{D}(\bm{\mathbf{L}}-\bm{\mathbf{L}}^{\prime})\,\left(N(L)+K\,P(L)\right),$
(6.3)

where $N(L)$ is the null-hypothesis reconstruction noise power spectrum as we
have discussed in the previous sections. In practical analysis, one works with
a discretized set of wave vectors $\bm{\mathbf{L}}$, and we may write

$\expectationvalue{\hat{\alpha}_{\bm{\mathbf{L}}}\hat{\alpha}_{\bm{\mathbf{L}}^{\prime}}}=A\,\delta_{\bm{\mathbf{L}},\bm{\mathbf{L}}^{\prime}}\,\left(N_{L}+K\,P_{L}\right),$
(6.4)

where $A$ is the total area of the sky that is surveyed.

An unbiased estimator for the normalization $K$ may then be constructed from
each $\hat{\alpha}(\bm{\mathbf{L}})$ by

$\hat{K}_{\bm{\mathbf{L}}}=P_{L}^{-1}(A^{-1}|\hat{\alpha}_{\bm{\mathbf{L}}}|^{2}-N_{L}).$
(6.5)

To take full advantage of the information, these estimators can then be
optimally combined to make a single estimator:

$\hat{K}=\sum_{\bm{\mathbf{L}}}\,w_{\bm{\mathbf{L}}}\,\hat{K}_{\bm{\mathbf{L}}},$
(6.6)

The weights $w_{\bm{\mathbf{L}}}$ are chosen to minimize the variance of
$\hat{K}$ in the absence of polarization rotation (i.e. $K=0$):

$w_{\bm{\mathbf{L}}}=\frac{\sum_{\bm{\mathbf{L}}^{\prime}}[\mathcal{C}_{K}^{-1}]_{\bm{\mathbf{L}},\bm{\mathbf{L}}^{\prime}}}{\sum_{\bm{\mathbf{L}}^{\prime},\bm{\mathbf{L}}^{\prime\prime}}[\mathcal{C}_{K}^{-1}]_{\bm{\mathbf{L}}^{\prime},\bm{\mathbf{L}}^{\prime\prime}}},$
(6.7)

where $\mathcal{C}_{K}$ is the covariance matrix of $K_{\bm{\mathbf{L}}}$’s
for the various $\bm{\mathbf{L}}$’s. The variance of $\hat{K}$ is then

$\langle\hat{K}^{2}\rangle=\left[\sum_{\bm{\mathbf{L}},\bm{\mathbf{L}}^{\prime}}[\mathcal{C}_{K}^{-1}]_{\bm{\mathbf{L}},\bm{\mathbf{L}}^{\prime}}\right]^{-1}.$
(6.8)

In fact, the covariance matrix is diagonal:

$(\mathcal{C}_{K})_{\bm{\mathbf{L}},\bm{\mathbf{L}}^{\prime}}=\delta_{\bm{\mathbf{L}},\bm{\mathbf{L}}^{\prime}}\langle\hat{K}_{\bm{\mathbf{L}}}^{2}\rangle,$
(6.9)

which is a consequence of statistical isotropy of the CMB anisotropies, namely
that no direction on the sky has special statistics. While statistical
isotropy holds for the primordial anisotropies, weak lensing by a specific
realization of the foreground large-scale structure induces couplings between
different Fourier modes [18]. For the purpose of estimating the signal-to-
noise ratio of $\hat{K}$, we neglect the mode couplings due to lensing. This
omission is justified, as we will show in Section 5 that any dependence of the
rotation quadratic estimators on the lensing potential $\phi$ is of order
$\mathcal{O}(\phi^{2})$; when random realizations of foreground lens are
averaged over, statistical isotropy is restored.

After some calculation, we have

$\langle\hat{K}_{\bm{\mathbf{L}}}^{2}\rangle=\left(\frac{N_{L}}{P_{L}}\right)^{2}.$
(6.10)

The variance of $\hat{K}$ under the null hypothesis is given by

$\langle\hat{K}^{2}\rangle^{-1}=\sum_{\bm{\mathbf{L}}}\left(\frac{P_{L}}{N_{L}}\right)^{2}.$
(6.11)

The square of the signal-to-noise ratio is

$\frac{K^{2}}{\langle\hat{K}^{2}\rangle}=\sum_{\bm{\mathbf{L}}}\left(\frac{KP_{L}}{N_{L}}\right)^{2}=\sum_{L}f_{\mathrm{sky}}(2L+1)\left(\frac{KP_{L}}{N_{L}}\right)^{2},$
(6.12)

where $f_{\mathrm{sky}}$ is the fraction of the full sky that is surveyed, and
$2L+1$ is the number of discrete $\bm{\mathbf{L}}$ vectors at a fixed
magnitude $L$.

##  7 String network models

As an example, we consider phenomenological models of string networks
developed in Ref. [17]. The theoretical framework is as follows. The Universe
is populated by many circular string loops oriented in random directions,
arising from a pseudoscalar field with anomaly coefficient $\mathcal{A}$. At
any given conformal time $\tau$, the number of string loops with radius
between $r$ and $r+\mathrm{d}r$ is given by the distribution function
$\nu(r,\tau)\,\mathrm{d}r$, specified by the model.

In order to compare models with different distributions $\nu$, a dimensionless
quantity is defined:

$\xi_{\infty}(\tau)=\frac{1}{a(\tau)^{2}H(\tau)^{2}}\int_{0}^{\infty}\mathrm{d}r\,2\pi
r\,\nu(r,\tau).$ (7.1)

Intuitively, if all the string loops at conformal time $\tau$ were rearranged
into a minimal number of Hubble-scale loops while conserving the energy
density of the string network, then $\xi_{\infty}(\tau)$ would be equal to the
number of such Hubble-scale loops per Hubble volume. Equivalently,
$\xi_{\infty}(\tau)$ is the total length of strings in units of the Hubble
scale per Hubble volume at conformal time $\tau$. Both models considered here
are _in scaling_. That is to say, the string network’s energy density scales
like the dominant energy density in the Universe, and
$\xi_{\infty}(\tau)\equiv\xi_{0}$ is independent of $\tau$.111Whether cosmic
string networks exhibit this in-scaling property is a matter of debate, with
some studies concluding the affirmative and others finding a logarithmic
growth in $\xi_{\infty}$. See Ref. [17] for a brief discussion and references
on this subject. The model parameter $\xi_{0}$ is also referred to as the
effective number of strings per Hubble volume.

The differences and additional parameters of these models are briefly stated
below.

  * •

Model I: All string loops have the same radius $\zeta_{0}$ (in units of the
Hubble scale at each redshift).

  * •

Model II: A fraction $f_{\mathrm{sub}}$ of the string loops are at the sub-
Hubble scale, with radii following a continuous, logarithmically flat
distribution between $\zeta_{\mathrm{min}}=0.1$ and $\zeta_{\mathrm{max}}=1$.
The remaining $(1-f_{\mathrm{sub}})$ fraction of the string loops have the
same radius $\zeta_{\mathrm{max}}=1$.

For both models, the predicted power spectrum of the rotation field due to the
string network is proportional to the product $\mathcal{A}^{2}\,\xi_{0}$:

$\expectationvalue{\alpha_{\bm{\mathbf{L}}}\,\alpha_{\bm{\mathbf{L}}^{\prime}}}=A\,\delta_{\bm{\mathbf{L}},\bm{\mathbf{L}}^{\prime}}\,\mathcal{A}^{2}\,\xi_{0}\,P_{L}(\cdots),$
(7.2)

where the fiducial power spectrum shape $P_{L}(\cdots)$ does not depend on
either $\mathcal{A}$ or $\xi_{0}$, but only on $\zeta_{0}$ for Model I or
$f_{\mathrm{sub}}$ for Model II. The formalism we have laid down above is
therefore applicable to constraining the overall normalization parameter
$K=\mathcal{A}^{2}\,\xi_{0}$, given the value of the remaining model
parameter.

The overall SNR as in Eq. (6.12) accumulates as we include more discrete
$\bm{\mathbf{L}}$ modes. In Figure 2 and Figure 3, we plot the SNR of the
$\hat{K}$ estimator as a function of $L_{\mathrm{max}}$, the maximum magnitude
of $\bm{\mathbf{L}}$ included into the analysis, for the three experimental
setups we consider in Section 6 and the two models stated above. For the large
scale cutoff, we set $L_{\mathrm{min}}=2$ for Planck SMICA and
$L_{\mathrm{min}}=30$ for Simons Observatory and CMB-S4. The survey sky
coverage $f_{\mathrm{sky}}$ used for each experimental setup can be found in
Table 2. The estimator $\hat{K}$ is defined from the HO estimator through Eq.
(6.5), and there would be a negligible effect on the SNR curves if the GMV
estimator were used instead of the HO estimator.

$\displaystyle 2$$\displaystyle 10$$\displaystyle 100$$\displaystyle
L_{\mathrm{max}}$$\displaystyle{10^{0}}$$\displaystyle{10^{1}}$$\displaystyle{10^{2}}$$\displaystyle{10^{3}}$Signal-
to-noise ratioPlanck SMICA$\displaystyle 30$$\displaystyle 100$$\displaystyle
L_{\mathrm{max}}$Simons Observatory$\displaystyle 30$$\displaystyle
100$$\displaystyle L_{\mathrm{max}}$CMB-S4

Figure 2: Signal-to-noise ratio (SNR) for detecting the polarization rotation
power spectrum normalization $K=\mathcal{A}^{2}\,\xi_{0}$, as a function of
$L_{\mathrm{max}}$ for string network Model I. The solid curves correspond to
$\mathcal{A}=1$ and $\xi_{0}=10$, and the dashed curves correspond to
$\mathcal{A}=1/3$ and $\xi_{0}=100$. The blue, orange, and green curves are
computed for three different string loop radii
$\zeta_{0}=1,10^{-0.5},10^{-1}$, respectively. The SNR saturates at $L_{\rm
max}\gtrsim 30$ for Planck SMICA with $L_{\mathrm{min}}=2$. For the ground-
based Simons Observatory and CMB-S4 with $L_{\mathrm{min}}=30$, the SNR
saturates at $L_{\rm max}\gtrsim 100$.

$\displaystyle 2$$\displaystyle 10$$\displaystyle 100$$\displaystyle
L_{\mathrm{max}}$$\displaystyle{10^{1}}$$\displaystyle{10^{2}}$$\displaystyle{10^{3}}$Signal-
to-noise ratioPlanck SMICA$\displaystyle 30$$\displaystyle 100$$\displaystyle
L_{\mathrm{max}}$Simons Observatory$\displaystyle 30$$\displaystyle
100$$\displaystyle L_{\mathrm{max}}$CMB-S4

Figure 3: Signal-to-noise ratio (SNR) for detecting the polarization rotation
power spectrum normalization $K=\mathcal{A}^{2}\,\xi_{0}$, as a function of
$L_{\mathrm{max}}$ for string network Model II. Assuming
$\zeta_{\mathrm{min}}=0.1$ and $\zeta_{\mathrm{max}}=1$, the blue, orange, and
green curves are computed for three different fractions of sub-Hubble scale
strings $f_{\mathrm{sub}}=0.2,0.6,0.9$, respectively. The SNR saturates at
$L_{\rm max}\gtrsim 20$ for Planck SMICA with $L_{\mathrm{min}}=2$. For the
ground-based Simons Observatory and CMB-S4 with $L_{\mathrm{min}}=30$, the SNR
saturates at $L_{\rm max}\gtrsim 100$.

As shown in Figure 2 and Figure 3, for all three experimental setups and both
models considered, the SNR accumulates as $L_{\mathrm{max}}$ increases until
it is saturated above $L_{\mathrm{max}}=100$. This is a result of the string
network producing most of the correlated rotations on very large angular
scales [8]. Therefore, it is sufficient to consider quadratic estimators
$\hat{\alpha}(\bm{\mathbf{L}})$ for $|\bm{\mathbf{L}}|<100$ for the purpose of
constraining cosmic strings signals.

For the Planck SMICA experiment, we have chosen $L_{\mathrm{min}}=2$ on
account of the analysis of Ref. [20]. An alternative choice of
$L_{\mathrm{min}}=5$ results in the reduction of the SNR by a little less than
a factor of two.

For the more simplistic model of uniform string sizes (Model I), Figure 2
shows that the SNR depends strongly on the value of the string loop radius
$\zeta_{0}$. As an improvement upon Model I, Model II exhibits a much weaker
dependence of the SNR on the model parameter $f_{\mathrm{sub}}$.

According to Figure 3, we note that the Planck SMICA data should have already
been sufficient for placing tight constraints on the more realistic string
model with a continuous distribution of string loop sizes (Model II). In
Section 8, we obtain such a constraint using the rotation power spectrum
estimated from the Planck 2015 data by Ref. [20].

Fixing $f_{\mathrm{sub}}=0.6$ and $L_{\mathrm{max}}=100$, we plot in Figure 4
the SNR contours in the $(\mathcal{A},\xi_{0})$ parameter plane.

$\displaystyle{10^{-2}}$$\displaystyle{10^{-1}}$$\displaystyle{10^{0}}$$\displaystyle{10^{1}}$$\displaystyle\mathcal{A}$$\displaystyle{10^{-1}}$$\displaystyle{10^{0}}$$\displaystyle{10^{1}}$$\displaystyle{10^{2}}$$\displaystyle{10^{3}}$$\displaystyle\xi_{0}$Signal-
to-noise ratio ($\displaystyle f_{\mathrm{sub}}=0.6$, $\displaystyle
L_{\mathrm{max}}=100$)

Figure 4: Signal-to-noise ratio of the $\hat{K}$ estimator in the
$(\mathcal{A},\xi_{0})$ parameter space, assuming $f_{\mathrm{sub}}=0.6$ and
$L_{\mathrm{max}}=100$. Solid lines mark the $3\sigma$ levels, while dashed
lines mark the $30\sigma$ levels. The blue, orange, and green levels refer to
the Planck SMICA, Simons Observatory, and CMB-S4 experiments, respectively.
The theoretically favored region of the parameter space, $0.1<\mathcal{A}<1$
and $1<\xi_{0}<100$, is indicated as the grey shaded area.

##  8 Planck 2015 bounds

In Ref. [20], Planck 2015 data has been used to extract a (binned)
birefringence power spectrum for $1\leq L\leq 730$. For both models studied in
Ref. [17] (stated in Section 7), we obtain posteriors for the relevant model
parameters using the Markov Chain Monte Carlo (MCMC) method, with particular
emphasis on the power spectrum normalization $\mathcal{A}^{2}\,\xi_{0}$.
Assuming $\mathcal{A}$ is an $\mathcal{O}(1)$ number, this allows us to place
a constraint on the number of cosmic strings in the Universe.

The theoretical birefringence power spectra predicted by the two models are
computed by direct integration according to Ref. [17]. To compare with the
model predictions, we only include the birefringence power spectrum for $2\leq
L\leq 305$, where the bin sizes are $\Delta L=1$ for $2\leq L\leq 30$ and
$\Delta L=50$ for $31\leq L\leq 330$. The dipole mode ($L=1$) has been
excluded from this comparison, as it is most susceptible to contamination by
residual foreground that has not been accounted for [20]. The high-$L$ modes
have also been excluded, as we anticipate little to no constraining power from
their inclusion (see Figure 2 and Figure 3).

Figure 5: 2D joint distribution and 1D marginalized distributions for the
power spectrum normalization $\mathcal{A}^{2}\,\xi_{0}$ and the string loop
radius $\zeta_{0}$ for string network Model I using Planck 2015 data from Ref.
[20].

We present in Figure 5 the 2D joint distribution and 1D marginalized
distributions for the power spectrum normalization $\mathcal{A}^{2}\,\xi_{0}$
and string loop radius $\zeta_{0}$ for Model I, obtained through MCMC
sampling. Although the marginalized distribution of $\mathcal{A}^{2}\,\xi_{0}$
is highly peaked near 0, it also has a wide tail, leading to poorer
constraining power of this model. Setting a flat prior in $\zeta_{0}$ for
$0.1<\zeta_{0}<1$, we report an upper bound of $\mathcal{A}^{2}\,\xi_{0}<8.0$
at the 95% confidence level.

Figure 6: 2D joint distribution and 1D marginalized distributions for the
power spectrum normalization $\mathcal{A}^{2}\,\xi_{0}$ and the fraction
$f_{\mathrm{sub}}$ of sub-Hubble scale string loops for Model II using data
from Ref. [20].

We present in Figure 6 the results of MCMC sampling for the parameters in the
more realistic Model II. Setting $\zeta_{\mathrm{min}}=0.1$ and
$\zeta_{\mathrm{max}}=1$ as in Figure 3, we report an upper bound of
$\mathcal{A}^{2}\,\xi_{0}<0.93$ at the 95% confidence level.

Both Figure 5 and Figure 6 show that the data is consistent with the absence
of a network of cosmic string loops. The model parameter $\zeta_{0}$ or
$f_{\mathrm{sub}}$ is poorly constrained by the data due to the absence of a
signal.

##  9 Conclusion

We have summarized three different constructions of the quadratic estimator
for direction dependent rotation of CMB polarization—the HO estimator [18],
the GMV estimator [19], and one that is derived from the maximum likelihood
estimator. We have shown that the latter two are equivalent in the context of
cosmic birefringence, and the first two are equivalent if temperature
anisotropies are not used (a general sufficient condition for such equivalence
is also stated). While the GMV estimator improves upon the HO estimator when
both temperature and polarization anisotropies are taken into account, the
improvement is numerically computed to be below 0.5% on all relevant angular
scales.

Although we have not seen a substantial reduction in the reconstruction noise
by replacing the HO estimator with the GMV estimator, it might still be
preferable to adopt the former for two reasons. First, the GMV estimator (as
well as its reconstruction noise spectrum) is computed from a single integral
(see Eq. (2.24) and Eq. (2.27)), as opposed to one integral for every pair of
anisotropy observables. Second, as pointed out in Ref. [19], the weights
$\bm{\mathbf{\Xi}}$ used in the GMV estimator are separable into sums of
products of functions of $\bm{\mathbf{l}}_{1}$ and functions of
$\bm{\mathbf{l}}_{2}$, allowing efficient computation via Fast Fourier
Transform (FFT) without any further approximations (such as by setting
$\tilde{C}_{l}^{TE}=0$ in some formulae).

We present the reconstruction noise spectra of the HO estimator for the Planck
mission, the Simons Observatory, as well as CMB-S4. This is compared to the
reconstruction noise spectra of the single-pair $EB$ and $TB$ estimators [15].
See Figure 1.

We have made SNR forecasts for phenomenological string network models (see
Figure 2 and Figure 3). It is found that contributions to the SNR are
dominated by the reconstruction of modes on large angular scales: $L\lesssim
30$ for Planck SMICA, and $L\lesssim 100$ for the Simons Observatory and
CMB-S4. When a continuous distribution of string loop sizes is accounted for,
detectability does not strongly depend on the fraction of sub-horizon string
loop sizes. In Figure 4, it can be seen that existing Planck data should
provide tight constraints in a significant portion of the theoretically
plausible parameter space at $\mathcal{A}^{2}\,\xi_{0}\gtrsim 1$.

The birefringence power spectrum extracted from Planck 2015 data [20] has
allowed us to place a constraint on $\mathcal{A}^{2}\,\xi_{0}$ using either of
the string network models. Model I, which postulates circular string loops of
identical radii, leads to an upper bound of $\mathcal{A}^{2}\,\xi_{0}<8.0$ at
the 95% confidence level, while the more realistic Model II, which postulates
circular string loops of a continuous distribution of radii, leads to a
stronger upper bound of $\mathcal{A}^{2}\,\xi_{0}<0.93$ at the 95% confidence
level. Model I leads to a weaker constraint because much of the parameter
space describes a Universe populated by sub-Hubble scale string loops
($\zeta_{0}<1$), and the signal in the power spectrum is suppressed by an
overall factor $\zeta_{0}$ [17]. On the other hand, most of the parameter
space in Model II describes a Universe with many Hubble-scale string loops,
leading to the prediction of a stronger signal. Simulations have shown that
less than 20% of the axion string loops are below the Hubble scale, while the
remaining strings wrap the simulation box [46], interpreted as string loops at
the Hubble scale in Model II. It is for this reason that we view Model II,
along with its bound on $\mathcal{A}^{2}\,\xi_{0}$, as a more accurate
description of the phenomenology of cosmic string networks. In the future,
more detailed numerical simulations of cosmic string network formation should
provide more accurate models of the rotation field. The analysis presented in
this work can also be performed for other models of cosmic birefringence, such
as ALP domain wall scenarios presented in Ref. [47].

Third- and fourth-generation CMB experiments from the ground are expected to
improve upon the SNR of the Planck mission by another one and two orders of
magnitude, respectively, and will be able to either discover or rule out such
cosmic string networks for most of the theoretically favored region of the
parameter space (see Figure 4).

We have justified the simplifying assumption that for estimating polarization
rotation the lensed temperature and polarization anisotropies may be taken to
be Gaussian random fields, by showing that the lensing potential $\phi$, which
introduces mode couplings to the anisotropies between different Fourier modes,
cause no bias at $\mathcal{O}(\phi)$ in estimating rotation.

While the quadratic estimators studied in this work can statistically detect a
network of strings, it is sensitive to the combination
$\mathcal{A}^{2}\,\xi_{0}$. The shape of the rotation power spectrum is
expected to be nearly universal for $L\lesssim 50$, whether there are many
small string loops or fewer larger loops [17]. Breaking the degeneracy between
$\mathcal{A}$ and $\xi_{0}$ is crucial if one wishes to measure the axion-
photon anomaly coefficient $\mathcal{A}$, whose value will have profound
implications for particle physics at high energies as explained in Section 1.
This is in principle possible by measuring the discontinuity in the rotation
angle across individual strings [8]. Non-Gaussian statistics of the rotation
field are expected to hold this information, which quadratic estimators are
ill-equipped to probe. Moreover, linear stringy features imply phase
correlations between Fourier modes of the rotation field, the confirmation of
which will provide decisive evidence for topological defects as the origin of
the ALP field. Extracting this information calls for further development of
non-Gaussian statistical tests for polarization rotations in the CMB, or the
employment of neural-network-based reconstruction techniques [48].

## Acknowledgments

We thank Prateek Agrawal, Anson Hook and Junwu Huang for suggesting this
research, and for their valuable input during both the early stage and the
completion of this work. LD acknowledges the research grant support from the
Alfred P. Sloan Foundation (award number FG-2021-16495). SF is supported by
the Physics Division of Lawrence Berkeley National Laboratory.

## References

  * [1] S. M. Carroll, _Quintessence and the rest of the world: Suppressing long-range interactions_ , _Phys. Rev. Lett._ 81 (1998) 3067.
  * [2] S. M. Carroll, G. B. Field and R. Jackiw, _Limits on a lorentz- and parity-violating modification of electrodynamics_ , _Phys. Rev. D_ 41 (1990) 1231.
  * [3] D. Harari and P. Sikivie, _Effects of a Nambu-Goldstone boson on the polarization of radio galaxies and the cosmic microwave background_ , _Physics Letters B_ 289 (1992) 67.
  * [4] P. Svrcek and E. Witten, _Axions in string theory_ , _Journal of High Energy Physics_ 2006 (2006) 051 [hep-th/0605206].
  * [5] L. Hui, J. P. Ostriker, S. Tremaine and E. Witten, _Ultralight scalars as cosmological dark matter_ , _Physical Review D_ 95 (2017) .
  * [6] M. Kamionkowski, J. Pradler and D. G. E. Walker, _Dark Energy from the String Axiverse_ , _Phys. Rev. Lett._ 113 (2014) 251302 [1409.0549].
  * [7] R. Emami, D. Grin, J. Pradler, A. Raccanelli and M. Kamionkowski, _Cosmological tests of an axiverse-inspired quintessence field_ , _Phys. Rev. D_ 93 (2016) 123005 [1603.04851].
  * [8] P. Agrawal, A. Hook and J. Huang, _A CMB Millikan experiment with cosmic axiverse strings_ , _JHEP_ 07 (2020) 138 [1912.02823].
  * [9] G. Hooft, _Naturalness, chiral symmetry, and spontaneous chiral symmetry breaking_ , in _Recent Developments in Gauge Theories_ , G. Hooft, C. Itzykson, A. Jaffe, H. Lehmann, P. K. Mitter, I. M. Singer et al., eds., (Boston, MA), pp. 135–157, Springer US, (1980), DOI.
  * [10] E. Witten, _Dyons of Charge e theta/2 pi_ , _Phys. Lett. B_ 86 (1979) 283.
  * [11] Y. Minami and E. Komatsu, _New extraction of the cosmic birefringence from the planck 2018 polarization data_ , _Physical Review Letters_ 125 (2020) .
  * [12] T. Fujita, K. Murai, H. Nakatsuka and S. Tsujikawa, _Detection of isotropic cosmic birefringence and its implications for axionlike particles including dark energy_ , _Physical Review D_ 103 (2021) .
  * [13] M. Kamionkowski, _How to derotate the cosmic microwave background polarization_ , _Physical Review Letters_ 102 (2009) .
  * [14] V. Gluscevic, M. Kamionkowski and A. Cooray, _Derotation of the cosmic microwave background polarization: Full-sky formalism_ , _Phys. Rev. D_ 80 (2009) 023510 [0905.1687].
  * [15] A. P. S. Yadav, R. Biswas, M. Su and M. Zaldarriaga, _Constraining a spatially dependent rotation of the cosmic microwave background polarization_ , _Physical Review D_ 79 (2009) .
  * [16] V. Gluscevic, D. Hanson, M. Kamionkowski and C. M. Hirata, _First CMB constraints on direction-dependent cosmological birefringence from WMAP-7_ , _Phys. Rev. D_ 86 (2012) 103529 [1206.5546].
  * [17] M. Jain, A. J. Long and M. A. Amin, _Cmb birefringence from ultralight-axion string networks_ , _Journal of Cosmology and Astroparticle Physics_ 2021 (2021) 055.
  * [18] W. Hu and T. Okamoto, _Mass reconstruction with cosmic microwave background polarization_ , _The Astrophysical Journal_ 574 (2002) 566–574.
  * [19] A. S. Maniyar, Y. Ali-Haïmoud, J. Carron, A. Lewis and M. S. Madhavacheril, _Quadratic estimators for cmb weak lensing_ , _Physical Review D_ 103 (2021) .
  * [20] D. Contreras, P. Boubel and D. Scott, _Constraints on direction-dependent cosmic birefringence fromplanckpolarization data_ , _Journal of Cosmology and Astroparticle Physics_ 2017 (2017) 046–046.
  * [21] D. Contreras, P. Boubel and D. Scott, _Constraints on direction-dependent cosmic birefringence fromplanckpolarization data_ , _Journal of Cosmology and Astroparticle Physics_ 2017 (2017) 046–046.
  * [22] U. Seljak and M. Zaldarriaga, _Polarization of microwave background: Statistical and physical properties_ , in _33rd Rencontres de Moriond: Fundamental Parameters in Cosmology_ , pp. 107–116, 1998, astro-ph/9805010.
  * [23] B. D. Sherwin and T. Namikawa, _Cosmic birefringence tomography and calibration-independence with reionization signals in the cmb_ , 2021.
  * [24] G. Hinshaw, D. Larson, E. Komatsu, D. N. Spergel, C. L. Bennett, J. Dunkley et al., _Nine-year wilkinson microwave anisotropy probe ( wmap ) observations: Cosmological parameter results_ , _The Astrophysical Journal Supplement Series_ 208 (2013) 19.
  * [25] N. Aghanim, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday et al., _Planckintermediate results_ , _Astronomy & Astrophysics_ 596 (2016) A110.
  * [26] S. Adachi, M. A. O. Aguilar Faúndez, K. Arnold, C. Baccigalupi, D. Barron, D. Beck et al., _A measurement of the degree-scale cmb b-mode angular power spectrum with polarbear_ , _The Astrophysical Journal_ 897 (2020) 55.
  * [27] F. Bianchini, W. L. K. Wu, P. A. R. Ade, A. J. Anderson, J. E. Austermann, J. S. Avva et al., _Searching for anisotropic cosmic birefringence with polarization data from sptpol_ , _Physical Review D_ 102 (2020) .
  * [28] T. Namikawa, Y. Guan, O. Darwish, B. D. Sherwin, S. Aiola, N. Battaglia et al., _Atacama cosmology telescope: Constraints on cosmic birefringence_ , _Physical Review D_ 101 (2020) .
  * [29] S. K. Choi, M. Hasselfield, S.-P. P. Ho, B. Koopman, M. Lungu, M. H. Abitbol et al., _The atacama cosmology telescope: a measurement of the cosmic microwave background power spectra at 98 and 150 ghz_ , _Journal of Cosmology and Astroparticle Physics_ 2020 (2020) 045–045.
  * [30] C. M. Hirata and U. Seljak, _Analyzing weak lensing of the cosmic microwave background using the likelihood function_ , _Phys. Rev. D_ 67 (2003) 043001 [astro-ph/0209489].
  * [31] C. M. Hirata and U. Seljak, _Reconstruction of lensing from the cosmic microwave background polarization_ , _Physical Review D_ 68 (2003) .
  * [32] U. Seljak and M. Zaldarriaga, _Signature of gravity waves in the polarization of the microwave background_ , _Phys. Rev. Lett._ 78 (1997) 2054.
  * [33] M. Zaldarriaga and U. Seljak, _All-sky analysis of polarization in the microwave background_ , _Physical Review D_ 55 (1997) 1830–1840.
  * [34] M. Kamionkowski, A. Kosowsky and A. Stebbins, _A probe of primordial gravity waves and vorticity_ , _Phys. Rev. Lett._ 78 (1997) 2058.
  * [35] A. Lewis and A. Challinor, _Weak gravitational lensing of the cmb_ , _Physics Reports_ 429 (2006) 1–65.
  * [36] T. Namikawa, _Testing parity-violating physics from cosmic rotation power reconstruction_ , _Physical Review D_ 95 (2017) .
  * [37] A. Lewis, A. Challinor and D. Hanson, _The shape of the cmb lensing bispectrum_ , _Journal of Cosmology and Astroparticle Physics_ 2011 (2011) 018–018.
  * [38] G. Fabbian, A. Lewis and D. Beck, _Cmb lensing reconstruction biases in cross-correlation with large-scale structure probes_ , _Journal of Cosmology and Astroparticle Physics_ 2019 (2019) 057–057.
  * [39] A. Lewis and A. Challinor, _CAMB: Code for Anisotropies in the Microwave Background_ , Feb., 2011.
  * [40] Planck Collaboration, Aghanim, N., Akrami, Y., Ashdown, M., Aumont, J., Baccigalupi, C. et al., _Planck 2018 results - vi. cosmological parameters_ , _A &A_ 641 (2020) A6.
  * [41] Planck Collaboration, Akrami, Y., Ashdown, M., Aumont, J., Baccigalupi, C., Ballardini, M. et al., _Planck 2018 results - iv. diffuse component separation_ , _A &A_ 641 (2020) A4.
  * [42] P. Ade, J. Aguirre, Z. Ahmed, S. Aiola, A. Ali, D. Alonso et al., _The simons observatory: science goals and forecasts_ , _Journal of Cosmology and Astroparticle Physics_ 2019 (2019) 056–056.
  * [43] K. N. Abazajian, P. Adshead, Z. Ahmed, S. W. Allen, D. Alonso, K. S. Arnold et al., _Cmb-s4 science book, first edition_ , 2016.
  * [44] L. Pogosian, M. Shimon, M. Mewes and B. Keating, _Future cmb constraints on cosmic birefringence and implications for fundamental physics_ , _Physical Review D_ 100 (2019) .
  * [45] S. Ferraro and J. C. Hill, _Bias to cmb lensing reconstruction from temperature anisotropies due to large-scale galaxy motions_ , _Physical Review D_ 97 (2018) .
  * [46] M. Gorghetto, E. Hardy and G. Villadoro, _Axions from strings: the attractive solution_ , _Journal of High Energy Physics_ 2018 (2018) .
  * [47] F. Takahashi and W. Yin, _Kilobyte cosmic birefringence from alp domain walls_ , _Journal of Cosmology and Astroparticle Physics_ 2021 (2021) 007.
  * [48] E. Guzman and J. Meyers, _Reconstructing cosmic polarization rotation with resunet-cmb_ , 2021.