Further author information: (Send correspondence to G.C.)
G.C.: E-mail<EMAIL_ADDRESS>
# CiaoCiao WFS: sensing phase discontinuities at the Extremely Large Telescope

Giulia Carlà  INAF, Osservatorio Astrofisico di Arcetri, Largo Enrico Fermi 5,
I-50125 Firenze, Italy  Lorenzo Busoni  INAF, Osservatorio Astrofisico di
Arcetri, Largo Enrico Fermi 5, I-50125 Firenze, Italy  Simone Esposito  INAF,
Osservatorio Astrofisico di Arcetri, Largo Enrico Fermi 5, I-50125 Firenze,
Italy  Guido Agapito  INAF, Osservatorio Astrofisico di Arcetri, Largo Enrico
Fermi 5, I-50125 Firenze, Italy  Ronald Holzlöhner  European Southern
Observatory (ESO), Karl-Schwarzschild-Str. 2, 85748 Garching, Germany

###### Abstract

The upcoming extremely large telescopes will have to deal with the so-called
“pupil fragmentation” effect: for what concerns the Extremely Large Telescope
(ELT), the presence of thick spider legs supporting the secondary mirror may
induce unseen phase discontinuities across the pupil sectors that could limit
the performance of the adaptive optics correction. In this context, we propose
a wavefront sensor (WFS), the CiaoCiao WFS, consisting in a rotational
shearing interferometer to sense phase differences between the pupil sectors.
In this work, we present the CiaoCiao WFS concept and the first analyses
carried out through numerical simulations. In particular, we analyze the
performance of such a wavefront sensor in case the phase discontinuities are
induced by low-wind effect during observations with the Multiconjugate
adaptive Optics Relay For ELT Observations (MORFEO).

###### keywords:

Pupil fragmentation, ELT, MORFEO, wavefront sensor, low-wind effect

##  1 INTRODUCTION

The upcoming 25-39-m class of extremely large telescopes (ELTs) is currently
under design and realization with the aim to provide astrophysical
observations in near-infrared wavelengths with unprecedented resolutions. Due
to their large sizes, the pupils of such telescopes are inevitably segmented
and partially shaded by thick support structures. This leads to the effect
known as “pupil fragmentation” that can severely degrade the performance of
the adaptive optics (AO) systems that will equip the ELTs to reach the
diffraction limit. For what concerns the 39-m Extremely Large Telescope (ELT)
[1], such an effect is generated by the presence of 50-cm-wide spider legs
supporting the secondary mirror that creates thick shadows in the pupil, thus
a fragmentation of the pupil in six disjoint segments whose separation is
larger than the typical Fried parameter. Due to this configuration, the
wavefront sensor (WFS) subapertures under the spiders are blind and thus the
reconstructed wavefront may present discontinuities between the pupil sectors.
If not sensed and controlled, phase discontinuities can severely degrade the
AO performance even during the best atmospheric turbulence conditions.

Low Wind Effect (LWE) can be one of the contributors to the discontinuities
between pupil sectors that has first been observed in 2015 at the Very Large
Telescope during the commissioning of SPHERE [2, 3]. The sharp steps that were
observed in the optical path difference (OPD) around the VLT spiders have then
been explained through the cooling of the spiders due to the radiative
transfer between the spiders and the sky background [4]: the air flowing near
the spider cools down and, in case of low-wind conditions ($<$ 5 m/s), is not
well mixed; this effect leads to temperature gradients that induce
discontinuities in the OPD across the spiders. Several strategies have been
proposed for a passive mitigation of LWE-induced phase discontinuities [5],
but no solutions foreseeing direct measurements have been defined so far.

The Multiconjugate adaptive Optics Relay For ELT Observations (MORFEO) [6],
that will deliver multiconjugate adaptive optics (MCAO) correction to the
first-light Multi-AO Imaging Camera for Deep Observations (MICADO) [7] of the
ELT, will be equipped with Shack Hartmann (SH) WFSs. Being slope sensors, SH
WFSs are not able to measure phase jumps across the spiders of the ELT since
the distance between the sectors is larger than the atmospheric coherence
length. Pupil fragmentation can then be a major issue for MORFEO performance.

In this context, we propose the CiaoCiao WFS concept [8], evolved from the
work presented in Ref. [9], to sense phase discontinuities across the pupil
sectors of the ELT: by taking advantage of the rotational symmetry of the ELT
pupil, we aim at using a rotational shearing interferometer [10] to make
interference between two pupil sectors separated by a spider in order to sense
phase differences between adjacent sectors. Though an extended analysis of the
CiaoCiao WFS is being carried out - from both a numerical and an experimental
point of view - and is going to be provided in a future publication, in this
work we want to focus on the first results of the numerical analysis
concerning the capability of the CiaoCiao WFS to measure LWE-induced phase
differences between the ELT pupil sectors.

In Section 2, we present the CiaoCiao WFS concept; in Section 3, we go through
the steps of the Python software that we setup to perform numerical
simulations in order to analyze the performance of the CiaoCiao WFS; in
Section 4, we show first results on CiaoCiao WFS performance, including the
measurements of phase differences induced by LWE on MORFEO residual wavefront.

##  2 The CiaoCiao WFS

We propose to use a rotational shearing interferometer as WFS - the CiaoCiao
WFS \- to measure phase differences across the spiders of the ELT pupil. The
basic principle is to have a setup (e.g. a configuration based on a Mach-
Zender interferometer as schematized in Fig.1) that simultaneously acquires
the image of the pupil and the same image rotated by an angle $\Delta\theta$,
and that makes interference between the two.

Figure 1: Sketch of a simplifed setup of the CiaoCiao WFS. A cube beam
splitter divides the incident wavefront into two beams and sends them into two
interferometer arms. Each arm includes a Dove prism: the first one (bottom)
rotates the wavefront with respect to the wavefront in the other
interferometer arm; the second one (top) is stationary and is used to
compensate the changes in the optical path. The rotated beam and the reference
one overlap again at a second beam splitter and the interferograms are
acquired on a camera.

Considering the electric fields of the two beams:

$\begin{split}E_{1}(\textbf{r},t)&=A_{1}(\textbf{r})\cos(\omega
t+\varphi_{1}(\textbf{r}))\\\
E_{2}(\textbf{r},t)&=A_{2}(\textbf{r})\cos(\omega
t+\varphi_{2}(\textbf{r}))\,,\end{split}$ (1)

the measured intensity is the time-averaged squared magnitude of the
interfered electric field:

$\begin{split}I(\textbf{r})&=\langle|E_{1}(\textbf{r},t)+E_{2}(\textbf{r},t)|^{2}\rangle_{T}\\\
&=\dfrac{I_{1}(\textbf{r})}{2}+\dfrac{I_{2}(\textbf{r})}{2}+\sqrt{I_{1}(\textbf{r})I_{2}(\textbf{r})}\cos(\Delta\varphi(\textbf{r}))\,,\end{split}$
(2)

where

$\begin{split}\Delta\varphi(\textbf{r})&=\varphi_{1}(\textbf{r})-\varphi_{2}(\textbf{r})\\\
&=\varphi_{1}(r,\theta)-\varphi_{1}(r,\theta-\Delta\theta)\,.\end{split}$ (3)

Thus, by inverting Eq. (2), we can estimate the phase difference between two
pupil points having the same radial coordinates and the azimuthal coordinates
differing by the rotational angle $\Delta\theta$ that is defined from the
setup of the shearing interferometer. If we average the measurements of
$\Delta\varphi$ on the overlapping sectors, we get the estimate of the phase
difference between them (Fig. 2):

$\overline{\Delta\varphi}=\dfrac{1}{S}\int_{0}^{R}dr\int_{0}^{\Delta\theta}d\theta\left[\varphi(r,\theta)-\varphi(r,\theta-\Delta\theta)\right]\,,$
(4)

where $S$ is the overlapped area and $R$ is the pupil radius. Considering a
sensing wavelength $\lambda$, it follows the measurement of the optical path
difference (OPD) as

$\mathrm{OPD}=\dfrac{\lambda}{2\pi}\overline{\Delta\varphi}\,.$ (5)

It is worth noting that if $\Delta\theta$ is small, the CiaoCiao WFS measures
phase differences across the spiders; on the other hand, if $\Delta\theta$ is
60∘, adjacent sectors are completely overlapped in the interferogram and then
the CiaoCiao WFS measures the differential average phase between the ELT pupil
sectors.

Figure 2: Left: image of the pupil; center: image of the pupil rotated by 20∘;
right: overlapping sectors after interference.

##  3 Software implementation

We have written a Python software to perform numerical simulations to test the
CiaoCiao WFS performance. The software is based on POPPY [11] that we use to
create an optical model as stack of pupil layers. In our model we can include
the ELT pupil aperture, the petaled M4, the turbulence, low order Zernike and
also load maps (e.g. AO residual turbulence, LWE) from .fits files.

The steps of a single simulation are:

  * •

setup two identical optical models and rotate one of them by $\Delta\theta$;

  * •

generate interferograms with proper phase shifts for phase reconstruction -
the software implements the four step phase shifting interferometry algorithm;

  * •

extract phase difference maps and convert them to OPD maps;

  * •

estimate the jumps across adjacent regions that have been overlapped - i.e.
average OPD on the overlapping sectors;

  * •

derive petal signals $\textbf{p}_{m}$ from the cumulative sum of the measured
jumps $\textbf{j}_{m}$:

$\begin{split}\textbf{p}_{m}&=-\left[j_{1},j_{1}+j_{2},j_{1}+j_{2}+j_{3},\cdots\right]\\\
&=-\left[\left(p_{6}-p_{1}\right),\left(p_{6}-p_{1}\right)+\left(p_{1}-p_{2}\right),\left(p_{6}-p_{1}\right)+\left(p_{1}-p_{2}\right)+\left(p_{2}-p_{3}\right),\cdots\right]\\\
&=\left[p_{1}-p_{6},p_{2}-p_{6},p_{3}-p_{6},p_{4}-p_{6},p_{5}-p_{6},0\right]\,,\end{split}$
(6)

that is, in the algorithm the petals are estimated with an offset given by the
signal on the 6-th sector and the correction applied based on these
measurements will introduce a global piston given by $p_{6}$. Of course, this
is not an issue for what concerns the AO performance.

##  4 Numerical simulations

We show first tests concerning the CiaoCiao WFS measurement of pure pistons
injected on the ELT pupil sectors; then, we focus on the effects produced by
the LWE on MORFEO residual wavefront and on the capability of the CiaoCiao WFS
to compensate for such effects.

###  4.1 Injection of pure pistons

We consider pure pistons applied on the ELT pupil sectors, as shown in Fig. 3
(left plot). We assume the CiaoCiao WFS sensing at $\lambda$=2.2 $\mu$m and
with a rotation angle of 20∘. Considering the ordering of the sectors shown in
Fig. 2, the injected pistons are:

$\textbf{p}=[0,100,-200,300,-400,500]\>\mathrm{nm}\,.$ (7)

The CiaoCiao WFS measures the following jumps across the sectors (right plot
of Fig. 3):

$\textbf{j}_{m}=[500,-100,300,-500,700,-900]\>\mathrm{nm}\,,$ (8)

that lead to the estimation of the following signals on the sectors (from Eq.
(6)):

$\textbf{p}_{m}=[-500,-400,-700,-200,-900,-0]\>\mathrm{nm}\,.$ (9)

Thus, the measurements from the CiaoCiao WFS are correct up to a global piston
of -500 nm; indeed, as shown in Section 3, they have an intrinsic offset given
by the signal on the 6-th sector. In this case, the correction from the
CiaoCiao WFS will leave a residual global piston, that clearly does not
represent an issue for the AO performance.

Figure 3: Left: OPD map showing the pistons injected on the ELT pupil sectors
and the expected jumps (in the white boxes). Right: jumps measured by the
CiaoCiao WFS from the interferogram of the two pupil images having a
differential rotation angle of 20∘.

In Fig. 4, we apply pistons up to 1000 nm on the ELT pupil sectors:

$\textbf{p}=[0,200,-400,600,-800,1000]\>\mathrm{nm}\,.$ (10)

Figure 4: Left: OPD map showing the pistons injected on the ELT pupil sectors
and the expected jumps (in the white boxes). Right: jumps measured by the
CiaoCiao WFS from the interferogram of the two pupil images having a
differential rotation angle of 20∘. The injected pistons here lead to jumps
exceeding the dynamic range [-$\lambda$/2; +$\lambda$/2].

Though the pistons are within the dynamic range of the interferometric
measurements, the jumps exceed such range and the CiaoCiao WFS measures:

$\textbf{j}_{m}=[1000,-200,600,-1000,-800,400]\>\mathrm{nm}\,,$ (11)

that is, the last two measurements are wrong by -/+$\lambda$. This leads to
the following petals estimation:

$\textbf{p}_{m}=[0,200,-400,600,1400,1000]\>\mathrm{nm}\,,$ (12)

where the 5-th petal measurement is not correct.

This simple test shows the issue given by the modulo $\lambda$ ambiguity if
carrying out the OPD measurements with monochromatic light. Considering a
sensing wavelength of 2.2 $\mu$m - i.e. the longest wavelength that we can
access in the framework of the ELT/MORFEO -, the dynamic range is limited to
$\pm$ 1100 nm and jumps larger than this will lead to wrong measurements from
the CiaoCiao WFS. In this context, we are investigating the possibility to
implement multi-wavelength reconstruction techniques.

###  4.2 MORFEO residuals including the compensation for LWE

In this section, we report first results of the analysis concerning the
measurement of LWE-induced differential phase on MORFEO residuals.

In Fig. 5, we show the OPD maps from LWE that we used and that have been
derived from the Computational Fluid Dynamics (CFD) simulations that are
reported in Ref. [12]. The simulations consist of 99 steps of 0.2 s each and
the wind speed is 0.5 m/s. As a first analysis, we extract the average OPD
produced by LWE on each sector and we show the variation over time in Fig. 6:
the average OPD on the ELT pupil sectors determined by LWE is characterized by
a slow temporal evolution and the differential average phase between adjacent
sectors is up to $\sim$1.2 $\mu$m.

Figure 5: OPD map from CFD simulations of LWE [12] (one instantaneous frame).
The wind is directed along the x axis and its speed is 0.5 m/s. Figure 6:
Average OPD on each ELT pupil sector as a function of time. Colors represent
each sector and the ordering is the same as the one of Fig. 2. S1 and S3 are
superimposed, as well as S4 and S6.

For what concerns the OPD maps of MORFEO residual turbulence, we got them from
PASSATA end-to-end simulations [13]. The simulations consist of 2 s at 500 Hz
and implement a typical MORFEO configuration - 3 NGSs on 55” equilateral
asterism, 6 LGSs on 45” asterism, control and reconstruction based on POLC-
MMSE algorithms - with atmospheric turbulence characterized by a seeing of
0.644 arcsec (i.e. median conditions). The OPD maps of MORFEO residual
turbulence are shown in Fig. 7 (left image) and present a flat wavefront on
which we measure an RMS of 255 nm.

If we include LWE as input together with atmospheric turbulence within PASSATA
simulations, we get the OPD maps shown in the right image of Fig. 7: MORFEO
residual wavefront appears now to be “petalled” and the residual RMS goes to
1744 nm, showing that the LWE can severely affect the operation of MORFEO.

Figure 7: Left: MORFEO residual wavefront from PASSATA simulations (one frame)
that include only atmospheric turbulence as input. Right: MORFEO residual
wavefront obtained from PASSATA simulations (one frame) when also LWE is
included.

We test the CiaoCiao WFS on MORFEO residual wavefront maps including LWE, in
order to verify if it is able to restore the RMS expected for the considered
MORFEO configuration: since the impact of LWE is to introduce differential
pistons between the ELT pupil sectors on MORFEO residual wavefront, the effect
could be indeed easily detected by the CiaoCiao WFS having a setup with a
rotation angle of 60∘. One issue is the dynamic range: from the left plot of
Fig. 7, we can measure jumps between adjacent sectors that exceed the range
available with a 2.2 $\mu$m sensing wavelength. Therefore, we simulated a
setup including a dual wavelength reconstruction [14] in H band ($\lambda_{1}$
= 1.5 $\mu$m, $\lambda_{2}$ = 1.6 $\mu$m) that provides a synthetic wavelength
of 24 $\mu$m. In Fig. 8, we show the interferogram (left) obtained with a
rotation angle of 60∘ and the results of the simulation (right): the CiaoCiao
WFS measures pistons up to $\sim$ 2.7 $\mu$m that, used to apply a correction
to MORFEO residual wavefront, bring the residual RMS back to 256 nm.

Figure 8: Left: Interferogram obtained with a rotation angle of 60∘ on MORFEO
residuals including LWE. Right: Table reporting the RMS without and with LWE
included in the simulations, the piston measurements from the CiaoCiao WFS and
the residual RMS after the correction is applied based on the CiaoCiao WFS
signals.

Considering the above analysis, it is worth noting the following aspects that
we aim to investigate:

  * •

the CiaoCiao WFS should not have stringent requirements on sky coverage for
what concerns the measurement of LWE-induced differential pistons. Indeed, as
shown in Fig. 6 the temporal evolution of LWE is slow and long integration
times could be used allowing for the sensing on faint stars;

  * •

since LWE appears to ”apply” six pistons on the ELT pupil sectors, we may not
need a detector with high spatial resolution. We then intend to test different
spatial samplings;

  * •

a K-band WFS does not provide the dynamic range required by LWE-induced phase
differences. We did first tests with dual-wavelength reconstruction and we
want to investigate further multi-wavelength techniques.

##  5 Conclusions

We presented the CiaoCiao WFS, a rotational shearing interferometer to sense
phase discontinuities at the ELT. Thanks to the rotational symmetry of the ELT
pupil, it could be a simple solution to tackle the issue of pupil
fragmentation that, if not controlled, will severely affect the performance of
the adaptive optics systems. Depending on the rotation angle of the shearing
interferometer, the CiaoCiao WFS can sense either phase discontinuities at the
edge of the spiders or the differential average phase between the adjacent
sectors. We focused the analysis on the performance of the CiaoCiao WFS when
dealing with LWE. From PASSATA simulations of MORFEO residuals, which include
LWE as input together with atmospheric turbulence, LWE shows to be critical
for the operation of MORFEO as it introduces differential pistons on the
residual wavefront that determine an RMS of 1744 nm. Our first tests of the
CiaoCiao WFS performance when measuring the LWE-induced differential phase
show promising results: the correction of the piston signals reconstructed by
the CiaoCiao WFS reduces the residual RMS of MORFEO to 256 nm.

Several aspects are currently under investigation and will be the object of
future works: among these, multi-wavelength reconstruction techniques to
extend the dynamic range of the WFS, integration time and spatial sampling.

## References

  * [1] R. Tamai et al.  “The ESO’s ELT construction progress”  In _Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series_ 11445, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 2020, pp. 114451E  DOI: 10.1117/12.2562828
  * [2] J.-F. Sauvage et al.  “Low Wind Effect, the main limitation of the SPHERE instrument”  In _Adaptive Optics for Extremely Large Telescopes IV (AO4ELT4)_ , 2015, pp. E9
  * [3] J.-F. Sauvage et al.  “Tackling down the low wind effect on SPHERE instrument”  In _Adaptive Optics Systems V_ 9909, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 2016, pp. 990916  DOI: 10.1117/12.2232459
  * [4] R. Holzlöhner, S. Kimeswenger, W. Kausch and S. Noll  “Bolometric night sky temperature and subcooling of telescope structures”  In _A &A_ 645, 2021, pp. A32  DOI: 10.1051/0004-6361/202038853
  * [5] J. Milli et al.  “Low wind effect on VLT/SPHERE: impact, mitigation strategy, and results”  In _Adaptive Optics Systems VI_ 10703, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 2018, pp. 107032A  DOI: 10.1117/12.2311499
  * [6] L. Busoni et al.  “MORFEO enters final design phase”, pp. These proceedings
  * [7] R. Davies et al.  “MICADO: The Multi-Adaptive Optics Camera for Deep Observations”  In _The Messenger_ 182, 2021, pp. 17–21  DOI: 10.18727/0722-6691/5217
  * [8] L. Busoni et al.  “A rotational shearing interferometer to sense across-the-spider phase discontinuity at ELT”, pp. Wavefront Sensing in the VLT/ELT era VII
  * [9] R. Holzlöhner, A. Kellerer and S. Lèvêque  “Simulation Environment for Segment and Petal Phasing of Large Telescopes”  In _Video Memorie della Societa Astronomica Italiana_ 1.1, 2019, pp. 32  DOI: 10.36116/WFS2019.32
  * [10] D. Malacara  “Optical shop testing”  John Wiley & Sons, 2007
  * [11] “POPPY documentation”, https://poppy-optics.readthedocs.io/en/latest/
  * [12] R. Martins, C. Verinaud and C. Kleinclaus  “Transient wavefront error from cooled air downwind of telescope spiders”  In _Modeling, Systems Engineering, and Project Management for Astronomy X_ 12187, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 2022, pp. 1218712  DOI: 10.1117/12.2627683
  * [13] G. Agapito, A. Puglisi and S. Esposito  “PASSATA: object oriented numerical simulation software for adaptive optics”  In _Adaptive Optics Systems V_ 9909, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 2016, pp. 99097E  DOI: 10.1117/12.2233963
  * [14] K. Houairi and F. Cassaing  “Two-wavelength interferometry: extended range and accurate optical path difference analytical estimator”  In _Journal of the Optical Society of America A_ 26.12, 2009, pp. 2503–2511  DOI: 10.1364/JOSAA.26.002503


# Cosmic birefringence from Planck data release 4

P. Diego-Palazuelos<EMAIL_ADDRESS>Instituto de Física de Cantabria
(CSIC-Universidad de Cantabria), Avda. de los Castros s/n, E-39005 Santander,
Spain  Dpto. de Física Moderna, Universidad de Cantabria, Avda. de los Castros
s/n, E-39005 Santander, Spain  J. R. Eskilt<EMAIL_ADDRESS>Institute
of Theoretical Astrophysics, University of Oslo, P.O. Box 1029 Blindern,
N-0315 Oslo, Norway  Y. Minami  Research Center for Nuclear Physics, Osaka
University, Ibaraki, Osaka, 567-0047 Japan  M. Tristram  Université Paris-
Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France  R. M. Sullivan  Department of
Physics & Astronomy, University of British Columbia, 6224 Agricultural Road,
Vancouver, British Columbia, Canada  A. J. Banday  Université de Toulouse,
UPS-OMP, IRAP, F-31028 Toulouse cedex 4, France  CNRS, IRAP, 9 Av. colonel
Roche, BP 44346, F-31028 Toulouse cedex 4, France  R. B. Barreiro  Instituto
de Física de Cantabria (CSIC-Universidad de Cantabria), Avda. de los Castros
s/n, E-39005 Santander, Spain  H. K. Eriksen  Institute of Theoretical
Astrophysics, University of Oslo, P.O. Box 1029 Blindern, N-0315 Oslo, Norway
K. M. Górski  Jet Propulsion Laboratory, California Institute of Technology,
4800 Oak Grove Drive, Pasadena, California, U.S.A.  Warsaw University
Observatory, Aleje Ujazdowskie 4, 00-478 Warszawa, Poland  R. Keskitalo
Computational Cosmology Center, Lawrence Berkeley National Laboratory,
Berkeley, California, U.S.A.  Space Sciences Laboratory, University of
California, Berkeley, California, U.S.A.  E. Komatsu  Max Planck Institute for
Astrophysics, Karl-Schwarzschild-Str. 1, D-85748 Garching, Germany  Kavli
Institute for the Physics and Mathematics of the Universe (Kavli IPMU, WPI),
Todai Institutes for Advanced Study, The University of Tokyo, Kashiwa
277-8583, Japan  E. Martínez-González  Instituto de Física de Cantabria (CSIC-
Universidad de Cantabria), Avda. de los Castros s/n, E-39005 Santander, Spain
D. Scott  Department of Physics & Astronomy, University of British Columbia,
6224 Agricultural Road, Vancouver, British Columbia, Canada  P. Vielva
Instituto de Física de Cantabria (CSIC-Universidad de Cantabria), Avda. de los
Castros s/n, E-39005 Santander, Spain  I. K. Wehus  Institute of Theoretical
Astrophysics, University of Oslo, P.O. Box 1029 Blindern, N-0315 Oslo, Norway

###### Abstract

We search for the signature of parity-violating physics in the cosmic
microwave background, called cosmic birefringence, using the Planck data
release 4. We initially find a birefringence angle of $\beta=0.30\pm 0.11$
(68% C.L.) for nearly full-sky data. The values of $\beta$ decrease as we
enlarge the Galactic mask, which can be interpreted as the effect of polarized
foreground emission. Two independent ways to model this effect are used to
mitigate the systematic impact on $\beta$ for different sky fractions. We
choose not to assign cosmological significance to the measured value of
$\beta$ until we improve our knowledge of the foreground polarization.

##  I Introduction

Dark matter and dark energy in the Universe [1] may be a parity-violating
pseudoscalar field, $\phi$, which changes sign under inversion of spatial
coordinates [2, 3]. This field can couple to the electromagnetic tensor
$F_{\mu\nu}$ and its dual tensor $\tilde{F}^{\mu\nu}$ via a Chern-Simons term
in the Lagrangian density, $\frac{1}{4}g_{\phi\gamma}\phi
F_{\mu\nu}\tilde{F}^{\mu\nu}$ [4, 5], which makes the phase velocities of
right- and left-handed states of photons different; thus, the plane of linear
polarization rotates clockwise on the sky by an angle
$\beta=-\frac{1}{2}g_{\phi\gamma}\int dt~{}\partial{\phi}/\partial t$, where
$g_{\phi\gamma}$ is the coupling constant [6, 7, 8, 9]. The space filled with
$\phi$ therefore behaves as if it were a birefringent material. For this
reason, such an effect is often called “cosmic birefringence.”

Linear polarization of the cosmic microwave background (CMB) photons is
sensitive to $\beta$ [10]. The polarization pattern on the sky can be
decomposed into parity-even $E$ modes and parity-odd $B$ modes [11, 12]. The
correlation functions of polarization fields (or the power spectra $C_{\ell}$
in spherical harmonics space with angular wavenumber $\ell$), contain two
parity-even $EE$ and $BB$ auto-spectra, and one parity-odd $EB$ cross-
spectrum. Cosmic birefringence then yields
$C_{\ell}^{EB,\mathrm{o}}=\frac{1}{2}\sin(4\beta)(C_{\ell}^{EE}-C_{\ell}^{BB})$,
even when the intrinsic polarization contains no $EB$ [10, 13, 14, 15, 16].
Here, the superscript “o” denotes the observed value, whereas the $EE$ and
$BB$ power spectra on the right side are the ones before undergoing the cosmic
birefringence.

Recently, a weak signal of $\beta=0.35^{\circ}\pm 0.14^{\circ}$ (68% C.L.) was
reported with a statistical significance of $2.4\,\sigma$ [17], using analysis
of $C_{\ell}^{EB}$ data from the European Space Agency Planck mission high-
frequency instrument (HFI) public release 3 (PR3) [18]. Refs. [19, 20, 21, 22,
23, 24, 25, 26] discuss possible cosmological implications of this particular
measurement, and Refs. [27, 28, 29, 30, 31, 32, 33] give previous constraints
on $\beta$. In this paper, we apply the analysis method of Ref. [17] to the
Planck public release 4 (PR4), the so-called “NPIPE” reprocessing of the
Planck data [34], to further investigate the origin and robustness of this
signal with respect to the known systematics in the Planck data sets and data
cuts, as well as to the $EB$ correlation intrinsic to the polarized Galactic
emission. Throughout this paper, we quote uncertainties at the 68 % confidence
level (C.L.).

##  II The NPIPE Data

Details of the NPIPE reprocessing of the Planck data are reported in Ref.
[34]. Here, we briefly summarize the key aspects relevant to our analysis. The
NPIPE pipeline processes raw, uncalibrated detector data into polarized
frequency and detector-set maps in the HEALPix format [35]. It is the only
Planck pipeline that is designed to work with both low-frequency instrument
(LFI) [36] and HFI [18] data.

The NPIPE reprocessing achieved smaller noise by (1) including more data
acquired during each 4-minute repointing maneuver (which makes up 9 % of the
mission) between each of the 45-minute scans; and (2) better modeling the data
via a short (167 ms) baseline offset model for noise, suppressing degree-scale
noise residuals. HFI PR3 data were fitted with pointing period (between 35 and
75 min long) offsets. A multi-frequency polarization model used in calibration
greatly reduces large-scale polarization uncertainty but introduces a pipeline
transfer-function that suppresses CMB polarization power at $\ell<20$.
Finally, a second-order analog-to-digital conversion nonlinearity (ADCNL)
model is used.

The net effect of these differences is a scale-dependent reduction in the
total uncertainty of $EE$ and $BB$: (1) About 50 % lower noise power spectrum
($N_{\ell}$) at $\ell\sim 10$; (2) 20–30 % lower $N_{\ell}$ at $\ell\sim 100$;
and (3) 10–20 % lower $N_{\ell}$ at $\ell\sim 1000$ (this applies to
temperature as well).

Much like the HFI PR3, NPIPE fits and corrects for gain fluctuations, ADCNL,
bolometric transfer-function residuals and bandpass mismatch by fitting time-
domain templates while solving for the polarized map. The release is
accompanied by 600 total (signal, noise and systematics) simulations, each
with realistic beam and pipeline transfer-functions and improved noise
consistency.

##  III Analysis method

Our analysis pipeline is based upon Refs. [37, 38, 39], built independently by
four groups (JRE, PDP, YM, MT): JRE, YM and MT follow the original
implementation but with different $C_{\ell}$ estimation codes [40, 41, 42],
while PDP uses an alternative implementation based on the small-angle
approximation [43].

If we relied only on $C_{\ell}^{EB}$ of the CMB, it would not be possible to
distinguish between $\beta$ and miscalibration of the instrumental
polarization angle, $\alpha$; thus, $\beta$ and $\alpha$ would be degenerate
[44, 45, 46, 47]. However, since $\beta$ is proportional to the path length of
photons, the polarized emission from our Galaxy is only negligibly affected by
$\beta$. We can use this property to break degeneracy between $\beta$ and
$\alpha$ [37]. Specifically, $C_{\ell}^{EB}$ of the CMB yields the sum
$\alpha+\beta$, whereas foreground emission yields $\alpha$ in the absence of
$EB$ intrinsic to the foreground.

We use four polarization-sensitive channels of the HFI at central frequencies
of $\nu=(100,143,217,353)$ GHz. We split the detector sets on the focal plane
into A and B sets following the definition given in Ref. [34]. There are eight
miscalibration angles, $\alpha_{i}$, with
$i=100\textrm{A},\dots,353\textrm{B}$. We then exclude the auto-power spectra
of the same maps, e.g., $100\textrm{A}\times 100\textrm{A}$, to avoid
contamination of possible correlated noise. There are 28 unique pairs of A/B
sets for $C_{\ell}^{EE}$ and $C_{\ell}^{BB}$, and 56 unique pairs for
$C_{\ell}^{EB}$.

We use a set of four Galactic masks removing respectively 5, 10, 20, and 30 %
of the sky. Each mask is constructed by thresholding the 353-GHz polarization
and total intensity maps smoothed with a Gaussian with a full-width-at-half-
maximum (FWHM) of $10^{\circ}$, to avoid polarized foreground residuals and
potential intensity-to-polarization ($I\to P$) leakage residuals.

We also exclude pixels in which the carbon monoxide (CO) line is bright. While
CO is not polarized, the mismatch of detector bandpasses creates a spurious
polarization signal via $I\to P$ leakage. We exclude pixels where CO is
brighter than 45 KRJ km s-1. The Galactic and CO masks are apodized with a
$1^{\circ}$ FWHM Gaussian taper. While the CO strength varies over frequencies
and there is no CO at 143 GHz, we use a common CO mask for all frequencies to
simplify the analysis because the CO mask removes less than 5 % of the sky,
essentially on the Galactic plane where we expect the foreground model to fail
anyway. Finally, a common mask for point sources is constructed from the union
of the Planck point-source masks at 100, 143, 217, and 353 GHz.

The effective sky fraction in the covariance matrix of $C_{\ell}$ is given by
$f_{\mathrm{sky}}=N_{\mathrm{pix}}^{-1}(\sum_{j=1}^{N_{\mathrm{pix}}}w_{j}^{2})^{2}/(\sum_{j=1}^{N_{\mathrm{pix}}}w_{j}^{4})$,
where $w_{i}$ is the value of the (noninteger) apodized mask [48, 49] and
$N_{\mathrm{pix}}$ is the number of pixels. Our combined masks yield
$f_{\mathrm{sky}}=0.93$, $0.90$, $0.85$, $0.75$, and $0.63$ for 0, 5, 10, 20,
and 30 % Galactic masks, respectively.

We compute $C_{\ell}$ using PolSpice (JRE)
[40]111http://www2.iap.fr/users/hivon/software/PolSpice/, NaMaster (PDP and
YM) [41]222https://github.com/LSSTDESC/NaMaster, and Xpol (MT)
[42]333https://gitlab.in2p3.fr/tristram/Xpol. The $EB$ power spectra at low
($\ell\lesssim 10$) and high multipoles are sensitive to cosmic birefringence
up to the epochs of reionization (a redshift of $z\approx 10$) and decoupling
($z\approx 1090$), respectively [50, 51]. We follow previous work [29, 17] and
focus on the high-$\ell$ data, and bin $C_{\ell}$ from
$\ell_{\mathrm{min}}=51$ to $\ell_{\mathrm{max}}=1490$ with a spacing of
$\Delta\ell=20$. The number of bins is $N_{\mathrm{bins}}=72$.

JRE, YM and MT use a Markov chain Monte Carlo to evaluate $-2\ln
L=\sum_{b=1}^{N_{\mathrm{bins}}}(\vec{v}^{\sf
T}_{b}\mathbf{M}^{-1}_{b}\vec{v}_{b}+\ln|\mathbf{M}_{b}|)$, and obtain the
posterior distributions of $\beta$ and $\alpha_{i}$, while PDP minimizes
$-2\ln L$ analytically in the small-angle approximation [43]. Here, $b$ is the
index for bins, and
$\vec{v}_{b}\equiv\mathbf{A}\vec{C}_{b}^{\mathrm{o}}-\mathbf{B}\vec{C}_{b}^{\mathrm{CMB,th}}$
with
$\vec{C}_{b}^{\mathrm{o}}=(C_{b}^{E_{i}E_{j},\mathrm{o}}~{}C_{b}^{B_{i}B_{j},\mathrm{o}}~{}C_{b}^{E_{i}B_{j},\mathrm{o}})^{\sf
T}$ and $\vec{C}_{b}^{\mathrm{CMB,th}}$ being the beam-smoothed and binned
theoretical $\Lambda$ Cold Dark Matter ($\Lambda$CDM) best-fitting CMB spectra
[52]. The covariance matrix, $\mathbf{M}_{b}$, is binned from
$\mathbf{M}_{\ell}=\mathbf{A}\mathrm{cov}(\vec{C}_{\ell}^{\mathrm{o}},\vec{C}_{\ell}^{\mathrm{o}}{}^{\sf
T})\mathbf{A}^{\sf T}$, and we divide $\mathbf{M}_{b}$ by $f_{\mathrm{sky}}$
to account for the mask. Both $\mathbf{A}=\mathbf{A}(\alpha_{i},\alpha_{j})$
and $\mathbf{B}=\mathbf{B}(\alpha_{i},\alpha_{j},\beta)$ are block-diagonal
matrices defined in Refs. [39, 17].

##  IV The NPIPE Simulation

Not only $\alpha_{i}$, but also other instrumental systematics can create
spurious $TB$ and $EB$ correlations [45, 53]. To quantify the impact of the
known systematics on $\beta$, we use realistic simulations of the NPIPE
processing, which include beam systematics, gain calibration and bandpass
mismatches, ADCNL and the transfer-function correction, among others [34].

We use the A/B sets of CMB$+$noise$+$foreground realizations, with the
foreground being the Commander sky model [54]. To isolate the systematics
coupled to the CMB, we remove the beam-smoothed Commander maps from the
simulations. The foreground-removed map still contains systematics coupled to
the foreground, such as spurious polarization from the bandpass mismatch.

Without foregrounds we can determine either $\alpha_{i}$ or $\beta$, assuming
that the other vanishes. We find $\beta=-0.009^{\circ}\pm 0.003^{\circ}$ for
$f_{\mathrm{sky}}=0.93$, where the uncertainty is that of the mean of 100
realizations. Therefore the impact of the known systematics on $\beta$ is
negligible compared to the statistical uncertainty of the measurement,
$0.11^{\circ}$, reported in Table 1. From the simulations we also find
$\alpha_{i}=0.188^{\circ}\pm 0.009^{\circ}$, $-0.305^{\circ}\pm
0.007^{\circ}$, $0.047^{\circ}\pm 0.006^{\circ}$, $0.039^{\circ}\pm
0.005^{\circ}$, $-0.063^{\circ}\pm 0.008^{\circ}$, $0.020^{\circ}\pm
0.008^{\circ}$, $0.01^{\circ}\pm 0.03^{\circ}$, and $-0.06^{\circ}\pm
0.04^{\circ}$ for $i=100\textrm{A},\dots,353\textrm{B}$. These values do not
need to agree with those of the data because the (unknown) miscalibration of
the instrumental polarization angles and the cross-polarization response of
beams are not included in the simulations.

The $I\to P$ leakage gives $C_{\ell}^{EB}\propto C_{\ell}^{TT}$, whereas the
cross-polarization effect gives $C_{\ell}^{EB}\propto C_{\ell}^{EE}$. We find
that $C_{\ell}^{EB}$ of 100A and 100B resembles $C_{\ell}^{EE}$ of the CMB;
thus, $\alpha_{100\textrm{A}}$ and $\alpha_{100\textrm{B}}$ detected in the
simulations are due to the cross-polarization effect. Because the impact on
$\beta$ is small, we do not correct the values of $\alpha_{i}$ for these
systematics.

Table 1: Cosmic birefringence and miscalibration angles in units of degrees with $1\,\sigma~{}(68\,\%)$ uncertainties after correcting for the foreground $EB$ correlation. The foreground $EB$ amplitudes, $A_{\ell}$, in four multipole bins are also shown.  $f_{\mathrm{sky}}$ | 0.93 | 0.90 | 0.85 | 0.75 | 0.63
---|---|---|---|---|---
$\beta$ | $\phantom{-}0.36\pm 0.11$ | $\phantom{-}0.26\pm 0.14$ | $\phantom{-}0.14\pm 0.17$ | $\phantom{-}0.10\pm 0.21$ | $\phantom{-}0.29\pm 0.28$
$\alpha_{100\textrm{A}}$ | $-0.32\pm 0.13$ | $-0.17\pm 0.16$ | $-0.07\pm 0.19$ | $-0.01\pm 0.23$ | $-0.21\pm 0.29$
$\alpha_{100\textrm{B}}$ | $-0.43\pm 0.13$ | $-0.32\pm 0.16$ | $-0.20\pm 0.19$ | $-0.14\pm 0.22$ | $-0.28\pm 0.29$
$\alpha_{143\textrm{A}}$ | $\phantom{-}0.03\pm 0.11$ | $\phantom{-}0.13\pm 0.14$ | $\phantom{-}0.29\pm 0.18$ | $\phantom{-}0.40\pm 0.21$ | $\phantom{-}0.22\pm 0.28$
$\alpha_{143\textrm{B}}$ | $\phantom{-}0.15\pm 0.11$ | $\phantom{-}0.25\pm 0.14$ | $\phantom{-}0.37\pm 0.18$ | $\phantom{-}0.39\pm 0.22$ | $\phantom{-}0.21\pm 0.28$
$\alpha_{217\textrm{A}}$ | $-0.06\pm 0.11$ | $\phantom{-}0.10\pm 0.14$ | $\phantom{-}0.22\pm 0.17$ | $\phantom{-}0.21\pm 0.21$ | $\phantom{-}0.02\pm 0.28$
$\alpha_{217\textrm{B}}$ | $-0.07\pm 0.11$ | $\phantom{-}0.07\pm 0.14$ | $\phantom{-}0.23\pm 0.17$ | $\phantom{-}0.23\pm 0.21$ | $0.003\pm 0.28$
$\alpha_{353\textrm{A}}$ | $-0.19\pm 0.10$ | $-0.08\pm 0.13$ | $\phantom{-}0.12\pm 0.17$ | $\phantom{-}0.03\pm 0.21$ | $-0.09\pm 0.28$
$\alpha_{353\textrm{B}}$ | $-0.23\pm 0.11$ | $-0.10\pm 0.13$ | $\phantom{-}0.10\pm 0.17$ | $\phantom{-}0.02\pm 0.21$ | $-0.02\pm 0.29$
$10^{2}A_{51-130}$ | $2.5^{+1.6}_{-1.4}$ | $5.7^{+2.5}_{-2.4}$ | $3.4^{+1.8}_{-1.7}$ | $18.8^{+6.0}_{-6.1}$ | $14.1^{+3.8}_{-3.7}$
$10^{2}A_{131-210}$ | $0.8^{+1.2}_{-0.6}$ | $4.3^{+5.3}_{-3.1}$ | $9.8^{+4.2}_{-4.0}$ | $4.3^{+3.4}_{-2.8}$ | $2.6^{+2.9}_{-1.8}$
$10^{2}A_{211-510}$ | $1.5^{+2.4}_{-1.1}$ | $7.3^{+6.0}_{-4.7}$ | $5.1^{+4.9}_{-3.4}$ | $1.6^{+2.2}_{-1.2}$ | $3.1^{+3.2}_{-2.1}$
$10^{2}A_{511-1490}$ | $6.2^{+5.7}_{-4.1}$ | $4.2^{+4.2}_{-2.9}$ | $6.2^{+6.7}_{-4.3}$ | $4.9^{+5.3}_{-3.4}$ | $5.8^{+5.3}_{-3.8}$

##  V Results

Figure 1:  Constraints on $\beta$ for various values of $f_{\mathrm{sky}}$
with and without accounting for the foreground $EB$ correlations. For the
former, the dashed and dotted lines show corrections using the filament model
(Eq. (2)) and the Commander sky model, respectively. For the latter, the
results of four pipelines (JRE, PDP, YM, MT) are shown.

We show the measured values of $\beta$ in Fig. 1. The results of the four
pipelines agree, and we quote the results of the JRE pipeline throughout this
paper, unless otherwise noted. When we use all the data with
$f_{\mathrm{sky}}=0.93$, we find $\beta=0.30^{\circ}\pm 0.11^{\circ}$. When we
remove one frequency channel at a time, we find $\beta=0.34^{\circ}\pm
0.12^{\circ}$ (without 100 GHz), $0.38^{\circ}\pm 0.12^{\circ}$ (143 GHz),
$0.27^{\circ}\pm 0.13^{\circ}$ (217 GHz), and $0.08^{\circ}\pm 0.21^{\circ}$
(353 GHz). The cross-power spectra with 353 GHz make the largest contribution.

We observe a decreasing $\beta$ for smaller $f_{\mathrm{sky}}$. This shift of
$\beta$ goes along the degeneracy line of
$\alpha_{i}+\beta=\mathrm{constant}$. A possible explanation for this trend is
an $EB$ correlation intrinsic to the polarized dust emission [37, 38, 17, 55].
When the foreground $EB$ power spectrum $C_{\ell}^{EB,\mathrm{FG}}$ exists, we
have the relation
$C_{\ell}^{EB,\mathrm{FG,o}}=\frac{1}{2}\sin(4\alpha)(C_{\ell}^{EE,\mathrm{FG}}-C_{\ell}^{BB,\mathrm{FG}})+C_{\ell}^{EB,\mathrm{FG}}\cos(4\alpha)$
for a single channel [56]. We can formally rewrite this as

$C_{\ell}^{EB,\mathrm{FG,o}}=\sqrt{J_{\ell}^{2}+\left(C_{\ell}^{EB,\mathrm{FG}}\right)^{2}}\sin(4\alpha+4\gamma_{\ell})\,,$
(1)

where $J_{\ell}\equiv(C_{\ell}^{EE,\mathrm{FG}}-C_{\ell}^{BB,\mathrm{FG}})/2$
and $\tan(4\gamma_{\ell})\equiv C_{\ell}^{EB,\mathrm{FG}}/J_{\ell}$. Here,
$\gamma_{\ell}$ is an effective angle for the foreground $EB$. For
$C_{\ell}^{EB,\mathrm{FG}}\propto J_{\ell}$, $\gamma_{\ell}$ becomes
independent of $\ell$, $\gamma_{\ell}=\gamma$, and is degenerate with
$\alpha$. In this limit and $|\gamma|\ll 1$, the foreground does not yield
$\alpha$ but $\alpha+\gamma$, and we measure $\beta-\gamma$ [37]. The sum,
$\alpha+\beta$, is not affected.

As both positive $TE$ and $TB$ correlations are found for polarized dust
emission in the Planck data when the Galactic plane is masked [57, 58], we
expect $C_{\ell}^{EB,\mathrm{dust}}>0$ for the same mask [59]. This is also
confirmed by an independent analysis using the distribution of filaments of
neutral hydrogen atoms [55]. From the ratio of $C_{\ell}^{TB}$ and
$C_{\ell}^{TE}$ at 353 GHz, $\gamma_{\ell}\approx\gamma>0$ is inferred for
$\ell\lesssim 500$ when the Galactic plane is masked with
$f_{\mathrm{sky}}\approx 0.7$ [55]. On the other hand, $\gamma_{\ell}$
oscillates around zero for nearly full-sky data. Therefore, we expect the
observed $\beta$ (which measures some average of $\beta-\gamma_{\ell}$ over
$\ell$) to decrease as we enlarge the Galactic mask, as observed.

We estimate $\gamma=r^{EB,\mathrm{dust}}\sqrt{\xi}/[2(1-\xi)]$ [37], where
$r^{EB,\mathrm{dust}}$ is the dust $EB$ cross-correlation coefficient and
$\xi\equiv C_{\ell}^{BB,\mathrm{dust}}/C_{\ell}^{EE,\mathrm{dust}}\approx 0.5$
[60]. Ref. [55] suggests that
$r^{EB,\mathrm{dust}}=4|r^{TB,\mathrm{dust}}|\psi^{\mathrm{dust}}$ with
$r^{TB,\mathrm{dust}}\approx 0.05$ [57] and $\psi^{\mathrm{dust}}\approx
5^{\circ}$ as a model for the signed upper bound on $r^{EB,\mathrm{dust}}$
with $f_{\mathrm{sky}}\approx 0.7$. This gives $\gamma\lesssim 0.7^{\circ}$.
The observed shift in $\beta$ from nearly full-sky data to
$f_{\mathrm{sky}}\approx 0.7$ is $0.5^{\circ}$, consistent with the estimate.

##  VI Modeling the impact of the foreground $EB$ correlation

We apply two independent approaches to assess the foreground impact. The first
model for $C_{\ell}^{EB,\mathrm{dust}}$ is based on filaments of hydrogen
clouds producing the thermal dust emission and polarization [59, 55]. We use
this model because it is the only physical model for the foreground $EB$
available today. When the filaments and the magnetic field lines are perfectly
aligned, the model produces a positive $TE$ but no $TB$ or $EB$ correlations.
When they misalign by a small angle $\psi^{\mathrm{dust}}$, $TB$ and $EB$
correlations emerge with the same sign. We thus use the ansatz
$C_{\ell}^{EB,\mathrm{dust}}=A_{\ell}C_{\ell}^{EE,\mathrm{dust}}\sin(4\psi_{\ell}^{\mathrm{dust}})$
with $0\leq A_{\ell}\ll 1$ and
$\psi_{\ell}^{\mathrm{dust}}=\frac{1}{2}\arctan(C_{\ell}^{TB,\mathrm{dust}}/C_{\ell}^{TE,\mathrm{dust}})$.
In the model of Ref. [55] $A_{\ell}\to|r_{\ell}^{TB,\mathrm{dust}}|$ and
$C_{\ell}^{EE,\mathrm{dust}}\to\sqrt{C_{\ell}^{EE,\mathrm{dust}}C_{\ell}^{BB,\mathrm{dust}}}$;
however, we treat $A_{\ell}$ as a free parameter here, and use
$C_{\ell}^{EE,\mathrm{dust}}$ because the data suggest
$C_{\ell}^{EE,\mathrm{dust}}\propto C_{\ell}^{BB,\mathrm{dust}}$ and taking
the square-root of noisy data is not numerically stable. When the angles are
small, we obtain

$\gamma_{\ell}\simeq\frac{A_{\ell}C_{\ell}^{EE,\mathrm{dust}}}{C_{\ell}^{EE,\mathrm{dust}}-C_{\ell}^{BB,\mathrm{dust}}}\frac{C_{\ell}^{TB,\mathrm{dust}}}{C_{\ell}^{TE,\mathrm{dust}}}\,.$
(2)

We use this model to account for the possible impact of the dust $EB$
correlation.

We specifically modify $\mathbf{A}$ and $\mathbf{B}$ in $-2\ln L$ as
$\mathbf{A}=\left(\begin{array}[]{cc}-\vec{\Lambda}^{\sf
T}_{\ell}\mathbf{\Lambda}_{\ell}^{-1},1\end{array}\right)$ and

$\mathbf{B}=\left(\begin{array}[]{cc}\vec{R}^{\sf
T}(\alpha_{i}+\beta,\alpha_{j}+\beta)-\vec{\Lambda}^{\sf
T}_{\ell}\mathbf{\Lambda}_{\ell}^{-1}\mathbf{R}(\alpha_{i}+\beta,\alpha_{j}+\beta)\end{array}\right)\,,$
(3)

where $\mathbf{R}$ and $\vec{R}$ are defined in equations 8 and 9 of Ref.
[39], respectively. We define the new matrices,
${\bf\Lambda}_{\ell}=\mathbf{R}(\alpha_{i},\alpha_{j})+\mathbf{D}(\alpha_{i},\alpha_{j})\mathbf{F}_{\ell}$
and $\vec{\Lambda}_{\ell}^{\sf T}=\vec{R}^{\sf
T}(\alpha_{i},\alpha_{j})+\vec{D}^{\sf
T}(\alpha_{i},\alpha_{j})\mathbf{F}_{\ell}$, where

$\displaystyle\bf{D}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}-\cos(2\alpha_{i})\sin(2\alpha_{j})&-\sin(2\alpha_{i})\cos(2\alpha_{j})\\\
\sin(2\alpha_{i})\cos(2\alpha_{j})&\cos(2\alpha_{i})\sin(2\alpha_{j})\end{array}\right)\,,$
(6) $\vec{D}=\left(\begin{array}[]{c}\cos(2\alpha_{i})\cos(2\alpha_{j})\\\
-\sin(2\alpha_{i})\sin(2\alpha_{j})\end{array}\right)\,,\quad\mathbf{F}_{\ell}=\left(\begin{array}[]{cc}2\gamma_{\ell,j}&-2\gamma_{\ell,j}\\\
2\gamma_{\ell,i}&-2\gamma_{\ell,i}\end{array}\right)\,.$ (7)

We retain the possibility that $\gamma_{\ell,i}$ may depend on the $i$th
frequency channel.

We compute the ratios of $C_{\ell}^{XY,\mathrm{dust}}$ in Eq. (2) using the
353 GHz data. To reduce the scatter, we smooth $C_{\ell}^{XY}$ by applying a
one-dimensional Gaussian filter before computing the ratios. Since the dust
spectral energy distribution (SED) cancels, these ratios can be used to fit
the cross-power spectra of all the other frequency channels. In this model, we
expect $A_{\ell}$, hence $\gamma_{\ell}$, to be independent of frequencies. We
confirm that allowing for frequency-dependent $A_{\ell}$ yields similar
results; thus, we use the same $A_{\ell}$ for all the frequency combinations.
To account for possible dependence on $\ell$, we split $A_{\ell}$ into four
bins ($51\leq\ell\leq 130$, $131\leq\ell\leq 210$, $211\leq\ell\leq 510$ and
$511\leq\ell\leq 1490$).

In Fig. 1, we show $\beta$ from the simultaneous fit of $\alpha_{i}$, $\beta$,
and $\gamma_{\ell}$. The model brings $\beta$ back to positive values for all
$f_{\mathrm{sky}}$, confirming that the declining $\beta$ is caused by
$C_{\ell}^{EB,\mathrm{dust}}$. We report the numerical values in Table 1.

We also try another, completely different approach to assessing the foreground
impact using the NPIPE simulation of the Commander sky model [54] based on
power-law synchrotron and one-component modified blackbody (MBB) dust SED.
When the angular resolution of the foreground model is higher than that of the
target frequency channel, the foreground component is smoothed to match the
QuickPol [53] beam specific to the NPIPE data set. To avoid divergence in the
deconvolution, the Gaussian beam of $5^{\prime}$ FWHM present in Commander’s
dust component is maintained at 217 and 353 GHz. We measure
$C_{\ell}^{EB,\mathrm{FG}}$ from the simulations and include it in the PDP
pipeline with a free amplitude ${\cal D}$, which is fit simultaneously with
$\beta$ and $\alpha_{i}$ [61].

The measured $C_{\ell}^{EB,\mathrm{FG}}$ does not represent a signal-dominated
template for the foreground $EB$, but is dominated by the statistical
fluctuation with the variance proportional to
$C_{\ell}^{EE,\mathrm{FG}}C_{\ell}^{BB,\mathrm{FG}}/(2\ell+1)$. Thus, using
this in the fit requires great care; since the Commander sky model was taken
from the data, $C_{\ell}^{EB,\mathrm{FG}}$ fluctuates in the same way as the
data, resulting in surprisingly tight constraints on the parameters. We find
amplitude parameters of ${\cal D}=1.03\pm 0.01$, $1.02\pm 0.01$, $1.03\pm
0.01$, $1.01\pm 0.04$, and $1.13\pm 0.05$ in descending order of
$f_{\mathrm{sky}}$. This does not mean that we detect the foreground $EB$ with
high statistical significance, but merely shows that the data are consistent
with $C_{\ell}^{EB,\mathrm{FG}}$ taken from the simulation on a mode by mode
basis.

In Fig. 1, we show that the behavior of $\beta$ is qualitatively similar to
that from the filament model. The values are $\beta=0.16^{\circ}\pm
0.05^{\circ}$, $0.15^{\circ}\pm 0.05^{\circ}$, $0.19^{\circ}\pm 0.06^{\circ}$,
$0.12^{\circ}\pm 0.13^{\circ}$, and $0.34^{\circ}\pm 0.19^{\circ}$ in
descending order of $f_{\mathrm{sky}}$, which agree well with those of the
filament model except for $f_{\mathrm{sky}}=0.93$. This discrepancy could be
due to the complexity of the foreground emission near the Galactic plane that
is not captured by the power-law synchrotron and MBB dust SED in the Commander
model. The derived uncertainties in $\beta$ are smaller because
$C_{\ell}^{EB,\mathrm{FG}}$ is correlated well with the data, including the
statistical fluctuations, hence reducing (somewhat artificially) the
covariance in the likelihood.

Let us comment on another limitation of this approach. The existence of
$\alpha_{i}$ inevitably yields a spurious $C_{\ell}^{EB,\mathrm{FG}}$ in the
Commander map. If $\alpha_{i}$ varies over frequency channels $i$, the
foreground SED of polarization fields in the sky assumed for Commander is no
longer adequate for describing the data, further yielding a spurious $EB$.
Therefore we must be careful when interpreting $C_{\ell}^{EB,\mathrm{FG}}$ in
the Commander map. The modeling approach (Eq. (2)) does not suffer from this
issue, but the model itself may be limited. These two approaches are
complementary, and it is encouraging that they yield similar results. We adopt
$\beta$ from the filament model as our baseline results.

##  VII Conclusions

In this paper, we searched for the signal of cosmic birefringence in the
Planck PR4 data. First, we applied the methodology of Ref. [17], finding
$\beta=0.30^{\circ}\pm 0.11^{\circ}$ for nearly full-sky data. This agrees
with, and is more precise than, the previous estimate from the PR3 data set
[17].

We then expanded the methodology and explored the dependence of $\beta$ on
Galactic masks. We found a trend of decreasing $\beta$ for smaller
$f_{\mathrm{sky}}$, which can be understood as the effect of
$C_{\ell}^{EB,\mathrm{dust}}$ [37, 38, 17, 55]. Accounting for this effect in
two independent ways, we found that $\beta$ was positive for all
$f_{\mathrm{sky}}$. Which $f_{\mathrm{sky}}$ value should we choose for the
most robust determination of $\beta$? We cannot provide a definitive answer
until we improve our understanding of the foreground $EB$.

If confirmed as a cosmological signal, this would provide evidence for physics
beyond the standard model of elementary particles and fields, with profound
implications for fundamental physics. To make progress, we must search for
$\beta$ in independent data sets of on-going [62, 63, 64, 65, 66, 67] and
future [68, 69, 70, 71, 72] experiments. Since the largest statistical
significance is seen for nearly full-sky data, a full-sky mission such as
LiteBIRD [72] will play an important role.

An overarching goal of many planned CMB experiments is to find the signature
of primordial gravitational waves in $B$-mode polarization [73], which drives
requirements for experimental design. The $EB$ science may provide additional
requirements for, e.g., cross-polarization coupling and beam systematics.
Equally important is the need for high-fidelity end-to-end simulations, with
the $EB$ science in mind.

###### Acknowledgements.

We thank D. J. Watts, S. Clark and J. C. Hill for useful discussions and B.
Partridge for comments on the draft. Planck is a project of the European Space
Agency (ESA) with instruments provided by two scientific consortia funded by
ESA member states and led by Principal Investigators from France and Italy,
telescope reflectors provided through a collaboration between ESA and a
scientific consortium led and funded by Denmark, and additional contributions
from NASA (USA). Some of the results in this paper have been derived using the
HEALPix package. This research used resources of the National Energy Research
Scientific Computing Center (NERSC), a U.S. Department of Energy Office of
Science User Facility operated under Contract No. DE-AC02-05CH11231. Part of
the research was carried out at the Jet Propulsion Laboratory, California
Institute of Technology, under a contract with the National Aeronautics and
Space Administration (80NM0018D0004). JRE, HKE and IW acknowledge funding from
the European Research Council (ERC) under the Horizon 2020 Research and
Innovation Programme (Grant agreement No. 819478). The work of YM was
supported in part by the Japan Society for the Promotion of Science (JSPS)
KAKENHI, Grants No. JP20K14497. RS and DS acknowledge the support of the
Natural Sciences and Engineering Research Council of Canada. PDP acknowledges
financial support from the Formación del Profesorado Universitario (FPU)
Programme of the Spanish Ministerio de Ciencia, Innovación y Universidades,
and thanks the Max Planck Institute for Astrophysics for hospitality when this
work was being finalized. PDP, BB, EMG and PV thank the Spanish Agencia
Estatal de Investigación (AEI, MICIU) for the financial support provided under
the projects with references PID2019-110610RB-C21, ESP2017-83921-C2-1-R and
AYA2017-90675-REDC, co-funded with EU FEDER funds, and acknowledge supports
from Universidad de Cantabria and Consejería de Universidades, Igualdad,
Cultura y Deporte del Gobierno de Cantabria via the “Instrumentación y ciencia
de datos para sondear la naturaleza del universo” project, as well as from
Unidad de Excelencia María de Maeztu (MDM-2017-0765). The work of EK was
supported in part by JSPS KAKENHI Grant No. JP20H05850 and JP20H05859, and the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy - EXC-2094 - 390783311. The Kavli IPMU is
supported by World Premier International Research Center Initiative (WPI),
MEXT, Japan.

## References

  * Weinberg [2008] S. Weinberg,  _Cosmology_ (Oxford University Press, 2008).
  * Marsh [2016] D. J. E. Marsh, Phys. Rept. 643, 1 (2016), arXiv:1510.07633 [astro-ph.CO] .
  * Ferreira [2021] E. G. M. Ferreira, Astron. Astrophys. Rev. 29, 7 (2021), arXiv:2005.03254 [astro-ph.CO] .
  * Ni [1977] W.-T. Ni, Phys. Rev. Lett. 38, 301 (1977).
  * Turner and Widrow [1988] M. S. Turner and L. M. Widrow, Phys. Rev. D 37, 2743 (1988).
  * Carroll  _et al._ [1990] S. M. Carroll, G. B. Field, and R. Jackiw, Phys. Rev. D 41, 1231 (1990).
  * Carroll and Field [1991] S. M. Carroll and G. B. Field, Phys. Rev. D 43, 3789 (1991).
  * Harari and Sikivie [1992] D. Harari and P. Sikivie, Phys. Lett. B 289, 67 (1992).
  * Carroll [1998] S. M. Carroll, Phys. Rev. Lett. 81, 3067 (1998), arXiv:astro-ph/9806099 [astro-ph] .
  * Lue  _et al._ [1999] A. Lue, L.-M. Wang, and M. Kamionkowski, Phys. Rev. Lett. 83, 1506 (1999), arXiv:astro-ph/9812088 [astro-ph] .
  * Seljak and Zaldarriaga [1997] U. Seljak and M. Zaldarriaga, Phys. Rev. Lett. 78, 2054 (1997), arXiv:astro-ph/9609169 [astro-ph] .
  * Kamionkowski  _et al._ [1997] M. Kamionkowski, A. Kosowsky, and A. Stebbins, Phys. Rev. Lett. 78, 2058 (1997), arXiv:astro-ph/9609132 .
  * Feng  _et al._ [2005] B. Feng, H. Li, M. Li, and X. Zhang, Phys. Lett. B620, 27 (2005), arXiv:hep-ph/0406269 [hep-ph] .
  * Feng  _et al._ [2006] B. Feng, M. Li, J.-Q. Xia, X. Chen, and X. Zhang, Phys. Rev. Lett. 96, 221302 (2006), arXiv:astro-ph/0601095 [astro-ph] .
  * Liu  _et al._ [2006] G.-C. Liu, S. Lee, and K.-W. Ng, Phys. Rev. Lett. 97, 161303 (2006), arXiv:astro-ph/0606248 [astro-ph] .
  * Finelli and Galaverni [2009] F. Finelli and M. Galaverni, Phys. Rev. D 79, 063002 (2009), arXiv:0802.4210 [astro-ph] .
  * Minami and Komatsu [2020a] Y. Minami and E. Komatsu, Phys. Rev. Lett. 125, 221301 (2020a), arXiv:2011.11254 [astro-ph.CO] .
  * Planck Collaboration III [2020] Planck Collaboration III, Astron. Astrophys. 641, A3 (2020), arXiv:1807.06207 [astro-ph.CO] .
  * Fujita  _et al._ [2021a] T. Fujita, Y. Minami, K. Murai, and H. Nakatsuka, Phys. Rev. D 103, 063508 (2021a), arXiv:2008.02473 [astro-ph.CO] .
  * Fujita  _et al._ [2021b] T. Fujita, K. Murai, H. Nakatsuka, and S. Tsujikawa, Phys. Rev. D 103, 043509 (2021b), arXiv:2011.11894 [astro-ph.CO] .
  * Takahashi and Yin [2021] F. Takahashi and W. Yin, J. Cosmol. Astropart. Phys. 04 (2021), 007, arXiv:2012.11576 [hep-ph] .
  * Mehta  _et al._ [2021] V. M. Mehta, M. Demirtas, C. Long, D. J. E. Marsh, L. McAllister, and M. J. Stott, J. Cosmol. Astropart. Phys. 07 (2021), 033, arXiv:2103.06812 [hep-th] .
  * Nakagawa  _et al._ [2021] S. Nakagawa, F. Takahashi, and M. Yamada, Phys. Rev. Lett. 127, 181103 (2021), arXiv:2103.08153 [hep-ph] .
  * Alvey and Escudero Abenza [2021] J. Alvey and M. Escudero Abenza, Phys. Lett. B 823, 136752 (2021), arXiv:2106.04226 [hep-th] .
  * Choi  _et al._ [2021] G. Choi, W. Lin, L. Visinelli, and T. T. Yanagida, Phys. Rev. D 104, L101302 (2021), arXiv:2106.12602 [hep-ph] .
  * Obata [2021] I. Obata, arXiv e-prints (2021), arXiv:2108.02150 [astro-ph.CO] .
  * Hinshaw  _et al._ [2013] G. Hinshaw _et al._ (WMAP Collaboration), Astrophys. J. Suppl. 208, 19 (2013), arXiv:1212.5226 [astro-ph.CO] .
  * Kaufman  _et al._ [2016] J. P. Kaufman, B. G. Keating, and B. R. Johnson, Mon. Not. Roy. Astron. Soc. 455, 1981 (2016), arXiv:1409.8242 [astro-ph.CO] .
  * Planck Collaboration Int. XLIX [2016] Planck Collaboration Int. XLIX, Astron. Astrophys. 596, A110 (2016), arXiv:1605.08633 [astro-ph.CO] .
  * Adachi  _et al._ [2020a] S. Adachi _et al._ (POLARBEAR Collaboration), Astrophys. J. 897, 55 (2020a), arXiv:1910.02608 [astro-ph.CO] .
  * Bianchini  _et al._ [2020] F. Bianchini _et al._ (SPT Collaboration), Phys. Rev. D 102, 083504 (2020), arXiv:2006.08061 [astro-ph.CO] .
  * Namikawa  _et al._ [2020] T. Namikawa _et al._ (ACT Collaboration), Phys. Rev. D 101, 083527 (2020), arXiv:2001.10465 [astro-ph.CO] .
  * Choi  _et al._ [2020] S. K. Choi _et al._ (ACT Collaboration), J. Cosmol. Astropart. Phys. 12 (2020), 045, arXiv:2007.07289 [astro-ph.CO] .
  * Planck Collaboration Int. LVII [2020] Planck Collaboration Int. LVII, Astron. Astrophys. 643, A42 (2020), arXiv:2007.04997 [astro-ph.CO] .
  * Gorski  _et al._ [2005] K. M. Gorski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and M. Bartelman, Astrophys. J. 622, 759 (2005), arXiv:astro-ph/0409513 [astro-ph] .
  * Planck Collaboration II [2020] Planck Collaboration II, Astron. Astrophys. 641, A2 (2020), arXiv:1807.06206 [astro-ph.CO] .
  * Minami  _et al._ [2019] Y. Minami, H. Ochi, K. Ichiki, N. Katayama, E. Komatsu, and T. Matsumura, PTEP 2019, 083E02 (2019), arXiv:1904.12440 [astro-ph.CO] .
  * Minami [2020] Y. Minami, PTEP 2020, 063E01 (2020), arXiv:2002.03572 [astro-ph.CO] .
  * Minami and Komatsu [2020b] Y. Minami and E. Komatsu, PTEP 2020, 103E02 (2020b), arXiv:2006.15982 [astro-ph.CO] .
  * Chon  _et al._ [2004] G. Chon, A. Challinor, S. Prunet, E. Hivon, and I. Szapudi, Mon. Not. Roy. Astron. Soc. 350, 914 (2004), arXiv:astro-ph/0303414 .
  * Alonso  _et al._ [2019] D. Alonso, J. Sanchez, and A. Slosar (LSST Dark Energy Science), Mon. Not. Roy. Astron. Soc. 484, 4127 (2019), arXiv:1809.09603 [astro-ph.CO] .
  * Tristram  _et al._ [2005] M. Tristram, J. F. Macias-Perez, C. Renault, and D. Santos, Mon. Not. Roy. Astron. Soc. 358, 833 (2005), arXiv:astro-ph/0405575 .
  * de la Hoz  _et al._ [2021] E. de la Hoz, P. Diego-Palazuelos, E. Martínez-González, P. Vielva, R. B. Barreiro, and J. D. Bilbao-Ahedo, arXiv e-prints (2021), arXiv:2110.14328 [astro-ph.CO] .
  * Wu  _et al._ [2009] E. Y. S. Wu _et al._ (QUaD Collaboration), Phys. Rev. Lett. 102, 161302 (2009), arXiv:0811.0618 [astro-ph] .
  * Miller  _et al._ [2009] N. J. Miller, M. Shimon, and B. G. Keating, Phys. Rev. D 79, 103002 (2009), arXiv:0903.1116 [astro-ph.CO] .
  * Komatsu  _et al._ [2011] E. Komatsu _et al._ (WMAP Collaboration), Astrophys. J. Suppl. 192, 18 (2011), arXiv:1001.4538 [astro-ph.CO] .
  * Keating  _et al._ [2012] B. Keating, M. Shimon, and A. Yadav, Astrophys. J. 762, L23 (2012), arXiv:1211.5734 [astro-ph.CO] .
  * Hivon  _et al._ [2002] E. Hivon, K. M. Gorski, C. B. Netterfield, B. P. Crill, S. Prunet, and F. Hansen, The Astrophysical Journal 567, 2 (2002).
  * Challinor and Chon [2005] A. Challinor and G. Chon, Mon. Not. Roy. Astron. Soc. 360, 509 (2005), arXiv:astro-ph/0410097 .
  * Komatsu  _et al._ [2009] E. Komatsu _et al._ (WMAP Collaboration), Astrophys. J. Suppl. 180, 330 (2009), arXiv:0803.0547 [astro-ph] .
  * Sherwin and Namikawa [2021] B. D. Sherwin and T. Namikawa, arXiv e-prints (2021), arXiv:2108.09287 [astro-ph.CO] .
  * Planck Collaboration VI [2020] Planck Collaboration VI, Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Hivon  _et al._ [2017] E. Hivon, S. Mottet, and N. Ponthieu, Astron. Astrophys. 598, A25 (2017), arXiv:1608.08833 [astro-ph.CO] .
  * Planck Collaboration IV [2020] Planck Collaboration IV, Astron. Astrophys. 641, A4 (2020), arXiv:1807.06208 [astro-ph.CO] .
  * Clark  _et al._ [2021] S. E. Clark, C.-G. Kim, J. C. Hill, and B. S. Hensley, Astrophys. J. 919, 53 (2021), arXiv:2105.00120 [astro-ph.GA] .
  * Abitbol  _et al._ [2016] M. H. Abitbol, J. C. Hill, and B. R. Johnson, Mon. Not. Roy. Astron. Soc. 457, 1796 (2016), arXiv:1512.06834 [astro-ph.CO] .
  * Planck Collaboration XI [2020] Planck Collaboration XI, Astron. Astrophys. 641, A11 (2020), arXiv:1801.04945 [astro-ph.GA] .
  * Weiland  _et al._ [2020] J. L. Weiland, G. E. Addison, C. L. Bennett, M. Halpern, and G. Hinshaw, Astrophys. J. 893, 119 (2020), arXiv:1907.02486 [astro-ph.CO] .
  * Huffenberger  _et al._ [2020] K. M. Huffenberger, A. Rotti, and D. C. Collins, Astrophys. J. 899, 31 (2020), arXiv:1906.10052 [astro-ph.CO] .
  * Planck Collaboration Int. XXX [2016] Planck Collaboration Int. XXX, Astron. Astrophys. 586, A133 (2016), arXiv:1409.5738 [astro-ph.CO] .
  * Diego-Palazuelos  _et al._ [2021] P. Diego-Palazuelos _et al._ , Manuscript in preparation (2021).
  * Xu  _et al._ [2019] Z. Xu _et al._ , Astrophys. J. 891, 134 (2019), arXiv:1911.04499 [astro-ph.IM] .
  * Adachi  _et al._ [2020b] S. Adachi _et al._ (POLARBEAR Collaboration), Astrophys. J. 904, 65 (2020b), arXiv:2005.06168 [astro-ph.CO] .
  * Aiola  _et al._ [2020] S. Aiola _et al._ (ACT Collaboration), J. Cosmol. Astropart. Phys. 12 (2020), 047, arXiv:2007.07288 [astro-ph.CO] .
  * Dutcher  _et al._ [2021] D. Dutcher _et al._ (SPT-3G Collaboration), Phys. Rev. D 104, 022003 (2021), arXiv:2101.01684 [astro-ph.CO] .
  * Ade  _et al._ [2021a] P. A. R. Ade _et al._ (BICEP/Keck Collaboration), Phys. Rev. Lett. 127, 151301 (2021a), arXiv:2110.00483 [astro-ph.CO] .
  * Ade  _et al._ [2021b] P. A. R. Ade _et al._ (SPIDER Collaboration), arXiv e-prints (2021b), arXiv:2103.13334 [astro-ph.CO] .
  * Westbrook  _et al._ [2018] B. Westbrook _et al._ (POLARBEAR-2 and Simons Array Collaborations), Journal of Low Temperature Physics 193 (2018).
  * Hui  _et al._ [2018] H. Hui _et al._ (BICEP Array Collaboration), Proc. SPIE Int. Soc. Opt. Eng. 10708, 1070807 (2018), arXiv:1808.00568 [astro-ph.IM] .
  * Ade  _et al._ [2019] P. Ade _et al._ (Simons Observatory Collaboration), J. Cosmol. Astropart. Phys. 02 (2019), 056, arXiv:1808.07445 [astro-ph.CO] .
  * Carlstrom  _et al._ [2019] J. Carlstrom _et al._ , Bulletin of the AAS 51 (2019), https://baas.aas.org/pub/2020n7i209.
  * Hazumi  _et al._ [2019] M. Hazumi _et al._ (LiteBIRD Collaboration), J. Low Temp. Phys. 194, 443 (2019).
  * Kamionkowski and Kovetz [2016] M. Kamionkowski and E. D. Kovetz, Ann. Rev. Astron. Astrophys. 54, 227 (2016), arXiv:1510.06042 [astro-ph.CO] .