# Joint analysis of BICEP2/keck array and Planck data

A. Ade, N. Aghanim, Z. Ahmed, R. Aikin, D. Alexander, M. Arnaud, J.

Aumont, C. Baccigalupi, A. Banday, D. Barkats, et al.

► To cite this version:

A. Ade, N. Aghanim, Z. Ahmed, R. Aikin, D. Alexander, et al. Joint analysis of BICEP2/keck array and Planck data. Physical Review Letters, 2015, 114 (10), pp.101301. 10.1103/Phys-RevLett.114.101301. hal-01593518

HAL Id: hal-01593518

https://hal.science/hal-01593518

Submitted on 12 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Joint Analysis of BICEP2/Keck Array and Planck Data

P. A. R. Ade et al.*

(BICEP2/Keck and Planck Collaborations)

(Received 21 January 2015; published 9 March 2015)

We report the results of a joint analysis of data from BICEP2/Keck Array and Planck. BICEP2 and Keck Array have observed the same approximately 400 deg² patch of sky centered on RA 0 h, Dec. -57.5°. The combined maps reach a depth of 57 nK deg in Stokes Q and U in a band centered at 150 GHz. Planck has observed the full sky in polarization at seven frequencies from 30 to 353 GHz, but much less deeply in any given region (1.2 µK deg in Q and U at 143 GHz). We detect 150 × 353 cross-correlation in B modes at high significance. We fit the single- and cross-frequency power spectra at frequencies ≥ 150 GHz to a lensed-ΛCDM model that includes dust and a possible contribution from inflationary gravitational waves (as parametrized by the tensor-to-scalar ratio r), using a prior on the frequency spectral behavior of polarized dust emission from previous Planck analysis of other regions of the sky. We find strong evidence for dust and no statistically significant evidence for tensor modes. We probe various model variations and extensions, including adding a synchrotron component in combination with lower frequency data, and find that these make little difference to the r constraint. Finally, we present an alternative analysis which is similar to a map-based cleaning of the dust contribution, and show that this gives similar constraints. The final result is expressed as a likelihood curve for r, and yields an upper limit $r_{0.05} < 0.12$ at 95% confidence. Marginalizing over dust and r, lensing B modes are detected at $7.0\sigma$ significance.

DOI: 10.1103/PhysRevLett.114.101301

PACS numbers: 98.70.Vc, 04.80.Nn, 95.85.Bh, 98.80.Es

## I. INTRODUCTION

The cosmic microwave background (CMB) [1] is an essential source of information about all epochs of the Universe. In the past several decades, characterization of the temperature and polarization anisotropies of the CMB has helped to establish the standard cosmological model (ΛCDM) and to measure its parameters to high precision (see, for example, Refs. [2,3]).

An extension to the standard big bang model, inflation, postulates a short period of exponential expansion in the very early Universe, naturally setting the initial conditions required by ΛCDM, as well as solving a number of additional problems in standard cosmology. Inflation's basic predictions regarding the Universe's large-scale geometry and structure have been borne out by cosmological measurements to date (see Ref. [4] for a review). Inflation makes an additional prediction, the existence of a background of gravitational waves, or tensor mode perturbations [5-8]. At the recombination epoch, the inflationary gravitational waves (IGW) contribute to the anisotropy of the CMB in both total intensity and linear polarization. The amplitude of tensors is conventionally parametrized by *r*, the tensor-to-scalar ratio at a fiducial scale. Theoretical predictions of the value of *r* cover a very wide range. Conversely, a measurement of *r* can discriminate between models of inflation.

Tensor modes produce a small increment in the temperature anisotropy power spectrum over the standard ΛCDM scalar perturbations at multipoles $l \lesssim 60$; measuring this increment requires the large sky coverage traditionally achieved by space-based experiments, and an understanding of the other cosmological parameters. The effects of tensor perturbations on B-mode polarization is less ambiguous than on temperature or E-mode polarization over the range $l \lesssim 150$. The B-mode polarization signal produced by scalar perturbations is very small and is dominated by the weak lensing of E-mode polarization on small angular scales, making the detection of an IGW contribution possible [9-12].

Planck [13] was the third generation CMB space mission, which mapped the full sky in polarization in seven bands centered at frequencies from 30 to 353 GHz to a resolution of 33 to 5 arc min [14,15]. The Planck Collaboration has published the best limit to date on tensor modes using CMB data alone [3]: $r_{0.002} < 0.11$ (at 95% confidence) using a combination of Planck, SPT, and ACT temperature data, plus WMAP polarization, although the Planck *r* limit is model dependent, with running of the scalar spectral index or additional relativistic degrees of freedom being well-known degeneracies which allow larger values of *r*.

Interstellar dust grains produce thermal emission, the brightness of which increases rapidly from the 100-150 GHz frequencies favored for CMB observations, becoming dominant at $\ge 350$ GHz even at high galactic latitude. The dust grains align with the Galactic magnetic

<footnote>*Full author list given at the end of the article. Please cite as "BICEP2/Keck and Planck Collaborations."</footnote>
field to produce emission with a degree of linear polariza-
zation [16]. The observed degree of polarization depends
on the structure of the Galactic magnetic field along the line
of sight, as well as the properties of the dust grains (see, for
example, Refs. [17,18]). This polarized dust emission
results in both *E* mode and *B* mode, and acts as a potential
contaminant to a measurement of *r*. Galactic dust polari-
zation was detected by Archeops [19] at 353 GHz and by
WMAP [2,20] at 90 GHz.

BICEP2 was a specialized, low angular resolution
experiment, which operated from the South Pole from
2010 to 2012, concentrating 150 GHz sensitivity compa-
rable to Planck on a roughly 1% patch of sky at high
Galactic latitude [21]. The BICEP2 Collaboration pub-
lished a highly significant detection of B-mode polarization
in excess of the r = 0 lensed-ΛCDM expectation over the
range 30 < l < 150 in Ref. [22] (hereafter BK-I). Modest
evidence against a thermal Galactic dust component domi-
nating the observed signal was presented based on the cross
spectrum against 100 GHz maps from the previous BICEP1
experiment. The detected B-mode level was higher than
that projected by several existing dust models [23,24]
although these did not claim any high degree of reliability.

The *Planck* survey released information on the structure
of the dust polarization sky at intermediate latitudes [25],
and the frequency dependence of the polarized dust
emission at frequencies relevant to CMB studies [26].
Other papers argued that the BICEP2 region is significantly
contaminated by dust [27,28]. Finally, *Planck* released
information on dust polarization at high latitude [29]
(hereafter PIP-XXX), and in particular examined a field
centered on the BICEP2 region (but somewhat larger than
it) finding a level of polarized dust emission at 353 GHz
sufficient to explain the 150 GHz excess observed by
BICEP2, although with relatively low signal-to-noise.

*Keck Array* is a system of BICEP2-like receivers also
located at the South Pole. During the 2012 and 2013
seasons *Keck Array* observed the same field as BICEP2 in
the same 150 GHz frequency band [30] (hereafter BK-V).
Combining the BICEP2 and *Keck Array* maps yields *Q* and
*U* maps with rms noise of 57 nK in nominal 1 deg² pixels
—by far the deepest made to date.

In this Letter, we take cross spectra between the joint
BICEP2/*Keck* maps and all the polarized bands of *Planck*.
The structure is as follows. In Sec. II we describe the
preparation of the input maps, the expectations for dust, and
the power spectrum results. In Sec. III the main multi-
frequency cross-spectrum likelihood method is introduced
and applied to the data, and a number of variations from the
selected fiducial analysis are explored. Section IV describes
validation tests using simulations as well as an alternate
likelihood. In Sec. V we investigate whether there could be
decorrelation between the *Planck* and BICEP2/*Keck* maps
due to the astrophysics of dust and/or instrumental effects.
Finally we conclude in Sec. VI.

## II. MAPS TO POWER SPECTRA

A. Maps and preparation

We primarily use the BICEP2/*Keck* combined maps, as
described in BK-V. We also use the BICEP2-only and
*Keck*-only maps as a cross-check. The *Planck* maps used
for cross-correlation with BICEP2/*Keck* are the full-
mission polarized maps from the PR2 *Planck* science
release [31,32], a subset of which was presented in
PIP-XXX. We compute *Planck* single-frequency spectra
as the cross-power spectra of two data-split maps, in which
the data are split into two subsets with independent noise.
We consider three data split maps: (i) detector-set maps,
where the detectors at a given frequency are divided into
two groups, (ii) yearly maps, where the data from the first
and second years of observations are used for the two maps,
and (iii) half-ring maps, where the data from each pointing
period is divided in halves. To evaluate uncertainties due to
*Planck* instrumental noise, we use 500 noise simulations of
each map; these are the standard set of time-ordered data
noise simulations projected into sky maps (the full focal
plane 8 (FFP8) simulations defined in Ref. [33]).

While the *Planck* maps are filtered only by the instru-
ment beam (the effective beam defined in Refs. [34] and
[35]), the BICEP2/*Keck* maps are in addition filtered due
to the observation strategy and analysis process. In par-
ticular, large angular scales are suppressed anisotropically
in the BICEP2/*Keck* mapmaking process to avoid atmos-
pheric and ground-fixed contamination; this suppression is
corrected in the power spectrum estimate. In order to
facilitate comparison, we therefore prepare “*Planck* as
seen by BICEP2/*Keck*” maps. In the first step we use
the ANAFAST, ALTERALM, and SYNFAST routines from the
HEALPix [36] package [37] to resmooth the *Planck* maps
with the BICEP2/*Keck* beam profile, assuming azimuthal
symmetry of the beam. The coordinate rotation from
Galactic to celestial coordinates of the *T*, *Q*, and *U* maps
is performed using the ALTERALM routine in the HEALPix
package. The sign of the Stokes *U* map is flipped to convert
from the HEALPix to the IAU polarization convention. Next
we pass these through the “observing” matrix *R*, described
in Sect. VI. B of BK-I, to produce maps that include the
filtering of modes occurring in the data processing pipeline
(including polynomial filtering and scan-synchronous
template removal, plus deprojection of beam systematics).

Figure 1 shows the resmoothed *Planck* 353 GHz *T*, *Q*, and *U* maps before and after filtering. In both cases the BICEP2/*Keck* inverse variance apodization mask has been applied. This figure emphasizes the need to account for the filtering before any comparison of maps is attempted, either qualitative or quantitative.

B. Expected spatial and frequency spectra of dust

Before examining the power spectra it is useful to review
expectations for the spatial and frequency spectra of dust.
<figure><img src="image_2.png" /><figcaption>FIG. 1 (color). Planck 353 GHz *T*, *Q*, and *U* maps before (left) and after (right) the application of BICEP2/Keck filtering. In both cases the maps have been multiplied by the BICEP2/Keck apodization mask. The Planck maps are presmoothed to the BICEP2/Keck beam profile and have the mean value subtracted. The filtering, in particular the third order polynominal subtraction to suppress atmospheric pickup, removes large-angular scale signal along the BICEP2/Keck scanning direction (parallel to the right ascension direction in the maps here).</figcaption></figure>

Figure 2 of PIP-XXX shows that the dust *BB* (and *EE*) angular power spectra are well fit by a simple power law $D_{\ell} \propto \ell^{-0.42}$, where $D_{\ell} = C_{\ell} \ell(\ell + 1)/2\pi$, when averaging over large regions of sky outside of the Galactic plane. Section V B of the same paper states that there is no evidence for departure from this behavior for 1% sky patches, although the signal-to-noise ratio is low for some regions. Presumably we expect greater fluctuation from the mean behavior than would be expected for a Gaussian random field.

The spectral energy distribution (SED) of dust polarization was measured in Ref. [26] for 400 patches with 10° radius at intermediate Galactic latitudes. Figure 10 of this reference shows empirically that the mean polarized dust SED is described by a simple modified blackbody spectrum with $T_d = 19.6$ K and $\beta_d = 1.59 \pm 0.17$ to within an accuracy of a few percent over the frequency range 100–353 GHz. Within this frequency range variations in the two parameters are highly degenerate and the choice is made to hold $T_d$ fixed at the value obtained from a fit to the SED. The power spectrum estimation proceeds exactly as in

BK-I, including the matrix based purification operation to prevent *E* to *B* mixing. Figure 2 shows the results for BICEP2/Keck and Planck 353 GHz for *TT*, *TE*, *EE*, and *BB*. In all cases the error bars are the standard deviations of lensed-ΛCDM + noise simulations [38] and hence contain no sample variance on any other component. The results in the left column are autospectra, identical to those given in BK-I and BK-V—these spectra are consistent with lensed ΛCDM + noise except for the excess in *BB* for $\ell < 200$.

The right column of Fig. 2 shows cross spectra between two halves of the Planck 353 GHz data set, with three different splits shown. The Planck Collaboration prefers the use of cross spectra even at a single frequency to gain additional immunity to systematics and to avoid the need to noise debias autospectra. The *TT* spectrum is higher than ΛCDM around $\ell = 200$—presumably due to a dust contribution. The *EE* and *BB* spectra are noisy, but both appear to show an excess over ΛCDM for $\ell < 150$—again presumably due to dust. We note that these spectra do not appear to follow the power-law expectation mentioned in Sec. II B, but we emphasize that the error bars
<figure><img src="image_3.png" /></figure>

FIG. 2 (color). Single- and cross-frequency spectra between BICEP2/Keck maps at 150 GHz and Planck maps at 353 GHz. The red curves show the lensed-ΛCDM expectations. The left column shows single-frequency spectra of the BICEP2, *Keck Array* and combined BICEP2/Keck maps. The BICEP2 spectra are identical to those in BK-I, while the *Keck Array* and combined are as given in BK-V. The center column shows cross-frequency spectra between BICEP2/Keck maps and Planck 353 GHz maps. The right column shows Planck 353 GHz data-split cross spectra. In all cases the error bars are the standard deviations of lensed-ΛCDM + noise simulations and hence contain no sample variance on any other component. For *EE* and *BB* the χ² and χ (sum of deviations) versus lensed ΛCDM for the nine band powers shown is marked at upper and lower left (for the combined BICEP2/Keck points and DS1 × DS2, respectively). In the bottom row (for *BB*) the center and right panels have a scaling applied such that signal from dust with the fiducial frequency spectrum would produce signal with the same apparent amplitude as in the 150 GHz panel on the left (as indicated by the right-side y axes). We see from the significant excess apparent in the bottom center panel that a substantial amount of the signal detected at 150 GHz by BICEP2 and *Keck Array* indeed appears to be due to dust.

contain no sample variance on any dust component (Gaussian or otherwise).

The center column of Fig. 2 shows cross spectra between BICEP2/Keck and Planck maps. For TE one can use the T modes from BICEP2 and the E modes from Planck or vice versa and both options are shown. Since the T modes are very similar between the two experiments, these TE spectra look similar to the single-experiment TE spectrum which shares the E modes. The EE and BB cross spectra are the most interesting—there appears to be a highly significant detection of correlated B-mode power between 150 and

353 GHz, with the pattern being much brighter at 353, consistent with the expectation from dust. We also see hints of detection in the *EE* spectrum—while dust E modes are subdominant to the cosmological signal at 150 GHz, the weak dust contribution enhances the BK150 × P353 cross spectrum at $\ell \approx 100$.

The polarized dust SED model mentioned in Sec. II B implies that dust emission is approximately 25 times brighter in the Planck 353 GHz band than it is in the BICEP2/Keck 150 GHz band (integrating appropriately over the instrumental bandpasses). The expectation for a
<figure><img src="image_4.png" /><figcaption>FIG. 3 (color). *EE* (left column) and *BB* (right column) cross spectra between BICEP2/Keck maps and all of the polarized frequencies of Planck. In all cases the quantity plotted is $\ell(\ell + 1)C_l/2\pi$ in units of $\mu K_{\text{CMB}}^2$, and the red curves show the lensed-$\Lambda$CDM expectations. The error bars are the standard deviations of lensed-$\Lambda$CDM + noise simulations and hence contain no sample variance on any other component. Also note that the y-axis scales differ from panel to panel in the right column. The $\chi^2$ and $\chi$ (sum of deviations) versus lensed $\Lambda$CDM for the five band powers shown is marked at upper left. There are no additional strong detections of deviation from lensed $\Lambda$CDM over those already shown in Fig. 2 although BK150 $\times$ P217 shows some evidence of excess.</figcaption></figure>

dust-dominated spectrum is thus that the BK150 $\times$ P353 cross spectrum should have an amplitude 25 times that of BK150 $\times$ BK150, and P353 $\times$ P353 should be 25 times higher again. The y-axis scaling in the bottom row of Fig. 2 has been adjusted so that a dust signal obeying this rule will have equal apparent amplitude in each panel. We see that a substantial amount of the BK150 $\times$ BK150 signal indeed appears to be due to dust.

To make a rough estimate of the significance of deviation from lensed $\Lambda$CDM, we calculate $\chi^2$ and $\chi$ (sum of deviations) for each of the *EE* and *BB* spectra and show these in Fig. 2. For the nine band powers used the expectation value (standard deviation) for $\chi^2$ and $\chi$ are 9 (4.2) and 0 (3), respectively. We see that BK150 $\times$ BK150 and BK150 $\times$ P353 are highly significant in *BB*, while P353 $\times$ P353 has modest significance in both *EE* and *BB*.

Figure 3 shows *EE* and *BB* cross spectra between BICEP2/Keck and all of the polarized frequencies of Planck (also including the BICEP2/Keck autospectra). For the five band powers shown the expectation value (standard deviation) for $\chi^2$ and $\chi$ are 5 (3.1) and 0 (2.2), respectively. As already noted, the BK150 $\times$ BK150 and BK150 $\times$ P353 *BB* spectra show highly significant excesses. Additionally, there is evidence for excess *BB* in BK150 $\times$ P217 spectrum, and for excess *EE* in BK150 $\times$ P353. The other spectra in Fig. 3 show no strong evidence for excess, although we note that only one of the $\chi$ values is negative. There is weak evidence for excess in the BK150 $\times$ P70 *BB* spectrum but none in BK150 $\times$ P30 so this is presumably just a noise fluctuation.

There are a large number of additional Planck-only spectra, which are not plotted here. The noise on these is large and all are consistent with $\Lambda$CDM, with the possible exception of P217 $\times$ P353, where modest evidence for an excess is seen in both *EE* and *BB* (see, e.g., Fig. 10 of PIP-XXX).

## D. Consistency of BICEP2 and Keck Array spectra

The *BB* autospectra for BICEP2 and Keck Array in the lower left panel of Fig. 2 appear to differ by more than might be expected, given that the BICEP2 and Keck maps cover almost exactly the same region of sky. However, the error bars in this figure are the standard deviations of lensed-$\Lambda$CDM + noise simulations; while the signal is largely common between the two experiments the noise is not, and the signal-noise cross terms produce substantial additional fluctuation of the difference. The correct way to quantify this is to compare the difference of the real data to the pairwise differences of simulations, using common input skies that have power similar to that observed in the real data. This was done in Sec. 8 of BK-V and the BICEP2 and Keck maps were shown to be statistically compatible. In an analogous manner we can also ask if the B150 $\times$ P353 and K150 $\times$ P353 *BB* cross spectra shown in the bottom middle panel of Fig. 2 are compatible. Figure 4 shows the results. We calculate the $\chi^2$ and $\chi$ statistics on these
<figure><img src="image_5.png" /></figure>

FIG. 4. Differences of B150 × P353 and K150 × P353 BB cross spectra. The error bars are the standard deviations of the pairwise differences of signal + noise simulations that share common input skies. The probability to exceed the observed values of $\chi^2$ and $\chi$ statistics, as evaluated against the simulations, is quoted for band-power ranges 1–5 (20 < $\ell$ < 200) and 1–9 (20 < $\ell$ < 330). There is no evidence that these spectra are statistically incompatible.

difference spectra and compare to the simulated distributions exactly as in BK-V. The probability to exceed (PTE) the observed values is given in the figure for band powers 1–5 (20 < $\ell$ < 200) and 1–9 (20 < $\ell$ < 330). There is no evidence that these spectra are statistically incompatible.

### E. Alternative power spectrum estimation

We check the reliability of the power spectrum estimation with an alternative pipeline. The filtered and purified Planck and BICEP2 maps used to make the spectra shown in Fig. 2 are transformed back into the HEALPIX pixelization using cubic spline interpolation. The B-mode cross power is then computed with the XPOL [39] and PURECL [40] estimators. Figure 5 shows the difference between these alternative band powers and the standard band powers for the B150 × P353 BB cross spectrum. As in Fig. 4 the error bars are the standard deviations of pairwise differences of simulations, which share common input skies and have power similar to that observed in the real data. The agreement is not expected to be exact due to the differing band-power window functions, but the differences of the real band powers are consistent with those of the simulations.

## III. LIKELIHOOD ANALYSIS

### A. Algorithm

While it is conventional in plots like Fig. 2 to present band powers with symmetric error bars, it is important to appreciate that this is an approximation. The likelihood of an observed band power for a given model expectation value is generally an asymmetric function, which can be

<figure><img src="image_6.png" /></figure>

FIG. 5 (color). Differences of B150 × P353 BB cross spectra from the standard power spectrum estimator and alternate estimators. The error bars are the standard deviations of the pairwise differences of signal+noise simulations that share common input skies. The probability to exceed the observed values of $\chi^2$ and $\chi$ statistics, as evaluated against the simulations, is quoted for band-power ranges 1–5 (20 < $\ell$ < 200) and 1–9 (20 < $\ell$ < 330). We see that the differences of the real spectra are consistent with the differences of the simulations.

computed given knowledge of the noise level(s). To compute the joint likelihood of an ensemble of measured band-power values it is, of course, necessary to consider their full covariance—this is especially important when using spectra taken at different frequencies on the same field, where the signal covariance can be very strong.

We compute the band-power covariance using full simulations of signal-cross-signal, noise-cross-noise, and signal-cross-noise. From these, we can construct the covariance matrix for a general model containing multiple signal components with any desired set of SEDs. When we do this we deliberately exclude terms whose expectation value is zero, in order to reduce noise in the resulting matrix due to the limited number of simulated realizations.

To compute the joint likelihood of the data for any given proposed model we use the Hamimeche-Lewis [41] approximation (HL; see Sec. 9.1 of Ref. [42] for mentionation details). Here we extend the method to deal with single- and cross-frequency spectra, and the covariances thereof, in an analogous manner to the treatment of, for example, *TT*, *TE*, and *EE* in the standard HL method. The HL formulation requires that the band-power covariance matrix be determined for only a single “fiducial model.” We compute multidimensional grids of models explicitly and/or use COSMOC [43] to sample the parameter space.

### B. Fiducial analysis

As an extension of the simplest lensed-ΛCDM paradigm, we initially consider a two component model of IGW with amplitude $r$, plus dust with amplitude $A_d$ (specified at
<figure><img src="image_7.png" /></figure>

FIG. 6 (color). Likelihood results from a basic lensed-ΛCDM + r + dust model, fitting BB autospectra and cross spectra taken between maps at 150 GHz, 217, and 353 GHz. The 217 and 353 GHz maps come from Planck. The primary results (heavy black) use the 150 GHz combined maps from BICEP2/Keck. Alternate curves (light blue and red) show how the results vary when the BICEP2 and Keck Array only maps are used. In all cases a Gaussian prior is placed on the dust frequency spectrum parameter $\beta_d = 1.59 \pm 0.11$. In the right panel the two-dimensional contours enclose 68% and 95% of the total likelihood.

353 GHz and $\ell = 80$). [Here we assume that the spectral index of the tensor modes ($n_t$) is zero, and a scalar pivot scale of 0.05 Mpc$^{-1}$; all values of r quoted in this Letter are $r_{0.05}$ unless noted otherwise.] Figure 6 shows the results of fitting such a model to BB band powers taken between BICEP2/Keck and the 217 and 353 GHz bands of Planck, using band powers 1–5 ($20 < \ell < 200$). For the Planck single-frequency case, the cross spectrum of detector sets (DS1 × DS2) is used, following PIP-XXX. The dust is modeled as a power law $D_{\ell} \propto \ell^{-0.42}$, with free amplitude $A_d$ and scaling with frequency according to the modified blackbody model.

As discussed in Sec. II B the simple modified blackbody model is shown empirically in Ref. [26] to describe the mean polarized dust SED at mid-Galactic latitudes to an accuracy of a few percent over the frequency range 100–353 GHz, with variation of the $\beta_d$ parameter being sufficient to characterize the patch-to-patch variation. Since it is not possible to constrain $\beta_d$ using the BICEP2/Keck and Planck cross-spectral band powers alone a tight Gaussian prior $\beta_d = 1.59 \pm 0.11$ is imposed, the uncertainty being scaled from the observed patch-to-patch variation at intermediate Galactic latitudes in Ref. [26], as explained in PIP-XXX. This prior assumes that the SED of dust polarization at intermediate latitudes [26] applies to the high latitude BICEP2/Keck field. From dust astrophysics, we expect variations of the dust SED in intensity and polarization to be correlated [18]. We thus tested our assumption by measuring the $\beta_d$ of the dust total intensity in the BICEP2/Keck field using the template fitting analysis described in Ref. [44], and find the same value.

In Fig. 6 we see that the BICEP2 data produce an r likelihood that peaks higher than that for the Keck Array data. This is because for $\ell < 120$ the autospectrum $B_{150} \times B_{150}$ is higher than for $K_{150} \times K_{150}$, while the cross spectrum $B_{150} \times P_{353}$ is lower than $K_{150} \times P_{353}$ (see

Fig. 2). However, recall that both pairs of spectra $B_{150} \times B_{150}/K_{150} \times K_{150}$ and $B_{150} \times P_{353}/K_{150} \times P_{353}$ have been shown to be consistent within noise fluctuation (see Sec. II D). In Sec. IVA these likelihood results are also found to be compatible. Given the consistency between the two experiments, the combined result gives the best available measurement of the sky.

The combined curves (BK + P) in the left and center panels of Fig. 6 yield the following results: $r = 0.048_{-0.032}^{+0.035}$, $r < 0.12$ at 95% confidence, and $A_d = 3.3_{-0.8}^{+0.9}$. For r the zero-to-peak likelihood ratio is 0.38. Taking $\frac{1}{2}[1-f(-2 \log L_0/L_{\text{peak}})]$, where f is the $\chi^2$ cdf (for one degree of freedom), we estimate that the probability to get a number smaller than this is 8% if in fact $r = 0$. For $A_d$ the zero-to-peak ratio is $1.8 \times 10^{-6}$ corresponding to a smaller-than probability of $1.4 \times 10^{-7}$, and a $5.1\sigma$ detection of dust power.

The maximum likelihood model on the grid has parameters $r = 0.05$, $A_d = 3.30 \ \mu\text{K}^2$ (and $\beta_d = 1.6$). Computing the band-power covariance matrix for this model, we obtain a $\chi^2$ of 40.8. Using 28 degrees of freedom—5 band powers times 6 spectra, minus 2 fit parameters (since $\beta_d$ is not really free)—gives a PTE of 0.06. The largest contributions to $\chi^2$ come from the $P_{353} \times P_{353}$ spectrum shown in the lower right panel of Fig. 2.

### C. Variations from the fiducial data set and model

We now investigate a number of variations from the fiducial analysis to see what difference these make to the constraint on r.

*Choice of Planck single-frequency spectra:* Switching the Planck single-frequency spectra to use one of the alternative data splits (yearly or half-ring instead of detector set) makes little difference (see Fig. 7).

*Using only 150 and 353 GHz:* Dropping the spectra involving 217 GHz from consideration also has little effect (see Fig. 7).
<figure><img src="image_8.png" /><figcaption>FIG. 7 (color). Likelihood results when varying the data sets used and the model priors—see Sec. III C for details.</figcaption></figure>

**Using only BK150 × BK150 and BK150 × P353:** Also excluding the 353 GHz single-frequency spectrum from consideration makes little difference. The statistical weight of the BK150 × BK150 and BK150 × P353 spectra dominate (see Fig. 7).

**Extending the band-power range:** Going back to the base data set and extending the range of band powers considered to 1–9 (corresponding to $20 < \ell < 330$) makes very little difference—the dominant statistical weight is with the lower band powers (see Fig. 7).

**Including EE spectra:** We can also include in the fits the EE spectra shown in Fig. 3. PIP-XXX (Figs. 5 and A.3) shows that the level of EE from Galactic dust is on average around twice the level of BB. However, there are substantial variations in this ratio from sky patch to sky patch. Setting $EE/BB = 2$ we find that the constraint on $A_d$ narrows, while the $r$ constraint changes little; this latter result is also shown in Fig. 7. The maximum likelihood model on the grid is unchanged and its $\chi^2$ PTE is acceptable.

**Relaxing the $\beta_d$ prior:** Relaxing the prior on the dust spectral index to $\beta_d = 1.59 \pm 0.33$ pushes the peak of the $r$ constraint up (see Fig. 7). However, it is not clear if this looser prior is self-consistent; if the frequency spectral index varied significantly across the sky it would invalidate cross-spectral analysis, but there is strong evidence against such variation at high latitude, as explained in Sec. VA. Nevertheless, it is important to appreciate that the $r$ constraint curves shown in Fig. 6 shift left (right) when assuming a lower (higher) value of $\beta_d$. For $\beta_d = 1.3 \pm 0.11$ the peak is at $r = 0.021$ and for $\beta_d = 1.9 \pm 0.11$ the peak is at $r = 0.073$.

**Varying the dust power spectrum shape:** In the fiducial analysis the dust spatial power spectrum is assumed to be a power law with $D_\ell \propto \ell^{-0.42}$. Marginalizing over spectral indices in the range −0.8 to 0 we find little change in the $r$ constraint (see also Sec. IV B for an alternate relaxation of

the assumptions regarding the spatial properties of the dust pattern).

**Using Gaussian determinant likelihood:** The fiducial analysis uses the HL likelihood approximation, as described in Sec. III A. An alternative is to recompute the covariance matrix **C** at each point in parameter space and take $L = \det(\mathbf{C})^{-1/2}\exp[-(\mathbf{d}^T\mathbf{C}^{-1}\mathbf{d})/2]$, where **d** is the deviation of the observed band powers from the model expectation values. This results in an $r$ constraint which peaks slightly lower, as shown in Fig. 7. Running both methods on the simulated realizations described in Sec. IVA, indicates that such a difference is not unexpected and that there may be a small systematic downward bias in the Gaussian determinant method.

**Varying the HL fiducial model:** As mentioned in Sec. III A the HL likelihood formulation requires that the expectation values and band-power covariance matrix be provided for a single “fiducial model” (not to be confused with the “fiducial analysis” of Sec. III B). The results from HL are supposed to be rather insensitive to the choice of this model, although preferably it should be close to reality. Normally, we use the lensed-$\Lambda$CDM + dust simulations described in Sec. IVA below. Switching this to lensed $\Lambda$CDM + $r = 0.2$ produces no change on average in the simulations, although it does cause any given realization to shift slightly—the change for the real data case is shown in Fig. 7.

**Adding synchrotron:** BK-I took the WMAP K-band (23 GHz) map, extrapolated it to 150 GHz according to $\nu^{-3.3}$ (mean value within the BICEP2 field of the MCMC “Model f” spectral index map provided by WMAP [2]), and found a negligible predicted contribution ($r_{sync,150} = 0.0008 \pm 0.0041$). Figure 3 does not offer strong motivation to reexamine this finding—the only significant detections of correlated BB power are in the BK150 × P353 and, to a lesser extent, BK150 × P217 spectra. However, here we proceed to a fit including all the polarized bands of Planck (as shown in Fig. 3) and adding a synchrotron component to the base lensed-$\Lambda$CDM + noise + $r$ + dust model. We take synchrotron to have a power law spectrum $D_\ell \propto \ell^{-0.6}$ [23], with free amplitude $A_{sync}$, where $A_{sync}$ is the amplitude at $\ell = 80$ and at 150 GHz, and scaling with frequency according to $\nu^{-3.3}$. In such a scenario we can vary the degree of correlation that is assumed between the dust and synchrotron sky patterns. Figure 8 shows results for the uncorrelated and fully correlated cases. Marginalizing over $r$ and $A_d$ we find $A_{sync} < 0.0003 \ \mu\text{K}^2$ at 95% confidence for the uncorrelated case, and many times smaller for the correlated. This last is because once one has a detection of dust it effectively becomes a template for the synchrotron. This synchrotron limit is driven by the Planck 30 GHz band—we obtain almost identical results when adding only this band, and a much softer limit when not including it. If we instead assume synchrotron scaling of $\nu^{-3.0}$ the limit on $A_{sync}$ is
<figure><img src="image_9.png" /><figcaption>FIG. 8 (color). Likelihood results for a fit when adding the lower frequency bands of Planck, and extending the model to include a synchrotron component. The results for two different assumed degrees of correlation between the dust and synchrotron sky patterns are compared to those for the comparable model without synchrotron (see text for details).</figcaption></figure>

approximately doubled for the uncorrelated case and reduced for the correlated. (Because the DS1 × DS2 data split is not available for the *Planck* LFI bands we switch to Y1 × Y2 for this variant analysis, and so we compare to this case in Fig. 8 rather than the usual fiducial case.)

**Varying lensing amplitude:** In the fiducial analysis the amplitude of the lensing effect is held fixed at the ΛCDM expectation ($A_L = 1$). Using their own and other data, the *Planck* Collaboration quote a limit on the amplitude of the lensing effect versus the ΛCDM expectation of $A_L = 0.99 \pm 0.05$ [3]. Allowing $A_L$ to float freely, and using all nine band powers, we obtain the results shown in Fig. 9—there is only weak degeneracy between $A_L$ and both $r$ and $A_d$. Marginalizing over $r$ and $A_d$ we find $A_L = 1.13 \pm 0.18$ with a likelihood ratio between zero and peak of $3 \times 10^{-11}$. Using the expression given in Sec. III B this corresponds to a smaller-than probability of $2 \times 10^{-12}$, equivalent to a $7.0\sigma$ detection of lensing in the *BB* spectrum. We note this is the most significant to-date direct measurement of lensing in B-mode polarization.

<figure><img src="image_10.png" /><figcaption>FIG. 9 (color). Likelihood results for a fit allowing the lensing scale factor $A_L$ to float freely and using all nine band powers. Marginalizing over $r$ and $A_d$, we find that $A_L = 1.13 \pm 0.18$ and $A_L = 0$ is ruled out with $7.0\sigma$ significance.</figcaption></figure>

# IV. LIKELIHOOD VALIDATION

## A. Validation with simulations

We run the algorithm used in Sec. III B on ensembles of simulated realizations to check its performance. We first consider a model where $r = 0$ and $A_d = 3.6 \ \mu\text{K}^2$, this latter being close to the value favored by the data in a dust-only scenario [45]. We generate Gaussian random realizations using the fiducial spatial power law $D_\ell \propto \ell^{-0.42}$, scale these to the various frequency bands using the modified black-body law with $T_d = 19.6$ K and $\beta_d = 1.59$, and add to the usual realizations of lensed-ΛCDM + noise. Figure 10 shows some of the resulting $r$ and $A_d$ constraint curves, with the result for the real data from Fig. 6 overplotted. As expected, approximately 50% of the $r$ likelihoods peak above zero. The median 95% upper limit is $r < 0.075$. We find that 8% of the realizations have a ratio $L_0/L_{\text{peak}}$ less than the 0.38 observed in the real data, in agreement with the estimate in Sec. III B. Running these dust-only realizations for BICEP2 only and *Keck Array* only, we find that the shift in the maximum likelihood value of $r$ seen in the real data in Fig. 6 is exceeded in about 10% of the simulations.

The above simulations assume that the dust component follows on average the fiducial $D_\ell \propto \ell^{-0.42}$ spatial power law, and fluctuates around it in a Gaussian manner. To obtain sample dust sky patterns that may deviate from this behavior in a way which better reflects reality, we take the prelaunch version of the *Planck* Sky Model (PSM; version 1.7.8 run in “simulation” mode) [24] evaluated in the *Planck* 353 GHz band and pull out the same 352 $|b| > 35^\circ$ partially overlapping regions used in PIP-XXX. We then scale these to the other bands and proceed as before. Some of the regions have dust power orders of magnitude higher than the real data and we cut them out (selecting 139 regions with peak $A_d < 20 \ \mu\text{K}^2$). Figure 11 presents the results. The $r$ likelihoods will broaden as the level of $A_d$ increases, and we should therefore not be surprised if the fraction of realizations peaking at a value higher than the real data is increased compared to the simulations with

<figure><img src="image_11.png" /><figcaption>FIG. 10. Likelihoods for $r$ and $A_d$, using BICEP2/Keck and *Planck*, as plotted in Fig. 6, overplotted on constraints obtained from realizations of a lensed-ΛCDM + noise + dust model with dust power similar to that favored by the real data ($A_d = 3.6 \ \mu\text{K}^2$). Half of the $r$ curves peak at zero as expected.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIG. 11. Constraints obtained when adding dust realizations from the Planck Sky Model version 1.7.8 to the base lensed-ΛCDM + noise simulations. (Curves for 139 regions with peak $A_d < 20~\mu\text{K}^2$ are plotted.) We see that the results for $r$ are unbiased in the presence of dust realizations which do not necessarily follow the $\ell^{-0.42}$ power law or have Gaussian fluctuations about it.</figcaption></figure>

mean $A_d = 3.6~\mu\text{K}^2$. However, we still expect that on average 50% will peak above zero and approximately 8% will have an $L_0/L_{\text{peak}}$ ratio less than the 0.38 observed in the real data. In fact we find 57% and 7%, respectively, consistent with the expected values. There is one realization which has a nominal (false) detection of nonzero $r$ of $3.3\sigma$, although this turns out to also have one of the lowest $L_0/L_{\text{peak}}$ ratios in the Gaussian simulations shown in Fig. 10 (with which it shares the CMB and noise components), so this is apparently just a relatively unlikely fluctuation.

## B. Subtraction of scaled spectra

As previously mentioned, the modified blackbody model predicts that dust emission is 4% as bright in the BICEP2 band as it is in the Planck 353 GHz band. Therefore, taking the autospectra and cross spectra of the combined BICEP2/Keck maps and the Planck 353 GHz maps, as shown in the bottom row of Fig. 2, and evaluating $(\text{BK} \times \text{BK} - \alpha \text{BK} \times \text{P})/(1-\alpha)$, at $\alpha = \alpha_{\text{fid}}$ cleans out the dust contribution (where $\alpha_{\text{fid}} = 0.04$). The upper panel of Fig. 12 shows the result.

As an alternative to the full likelihood analysis presented in Sec. III B, we can instead work with the differenced spectra from above, a method we denote the “cleaning” approach. If $\alpha_{\text{fid}}$ were the true value, the expectation value of this combination over CMB and noise would have no dust contribution. However, dust would still contribute to its variance, but only through its 2-point function. In practice, we do not know $\alpha$ perfectly, and this uncertainty needs to be accounted for in a likelihood constructed from the differenced spectra. Our approach is to treat the differenced spectra as a form of data compression, and to compute the expectation value as a function of $r$, $A_d$, and $\beta_d$ at each point in parameter space [the dust dependence enters for $\alpha(\beta_d) \neq \alpha_{\text{fid}}$]. We use the method of Ref. [41], with a fiducial covariance matrix, to build a likelihood for

<figure><img src="image_13.png" /><figcaption>FIG. 12 (color). (Upper) BB spectrum of the BICEP2/Keck maps before and after subtraction of the dust contribution, estimated from the cross spectrum with Planck 353 GHz. The error bars are the standard deviations of simulations, which, in the latter case, have been scaled and combined in the same way. The inner error bars are from lensed-ΛCDM + noise simulations as in the previous plots, while the outer error bars are from the lensed-ΛCDM + noise + dust simulations. The red curve shows the lensed-ΛCDM expectation. (Lower) Constraint on $r$ derived from the cleaned spectrum compared to the fiducial analysis shown in Fig. 6.</figcaption></figure>

the difference spectra, and marginalize over $A_d$ and $\beta_d$, and hence $\alpha$, adopting the prior $\beta_d = 1.59 \pm 0.11$. This alternative likelihood has the advantage of being less sensitive to non-Gaussianity of the dust, since only the 2-point function of the dust affects the covariance of the differenced spectra close to $\alpha_{\text{fid}}$, while the full analysis may, in principle, be affected by the non-Gaussianity of the dust through 4-point contributions to power spectra covariances. This cleaning approach does, however, ignore the (small amount of) additional information available at other frequencies. The lower panel of Fig. 12 compares the result to the fiducial analysis with the full multispectra likelihood. It is clear from the widths of the likelihood curves that compressing the spectra to form the cleaned difference results in very little loss of information on $r$. The difference in peak values arises from the different data
treatments and is consistent with the scatter seen across
simulations. Finally, we note that one could also form a
combination (BK × BK − 2αBK × P + α²P × P)/(1 − α)²
in which dust does not enter at all for α = α<sub>fid</sub>. However,
the variance of this combination of spectra is large due to
the *Planck* noise levels, and likelihoods built from this
combination are considerably less constraining.

V. POSSIBLE CAUSES OF DECORRELATION

Any systematic error that suppresses the BK150 × P353 cross-frequency spectrum with respect to the BK150 × BK150 and P353 × P353 single-frequency spectra would cause a systematic upward bias on the *r* constraint. Here we investigate a couple of possibilities.

A. Spatially varying dust frequency spectrum

If the frequency dependence of polarized dust emission
varied from place to place on the sky, it would cause the
150 and 353 GHz dust sky patterns to decorrelate and
suppress the BK150 × P353 cross-frequency spectrum
relative to the single-frequency spectra. The assumption
made so far in this Letter is that such decorrelation is
negligible. In fact, PIP-XXX implicitly tests for such
variation in their Fig. 6, where the *Planck* single- and
cross-frequency spectra are compared to the modified
blackbody model (with the cross-frequency spectra plotted
at the geometric mean of their respective frequencies). This
plot is for an average over a large region of low foreground
sky (24%); however, note that if there were spatial variation
of the spectral behavior anywhere in this region it would
cause suppression of the cross-frequency spectra with
respect to the single-frequency spectra.

PIP-XXX also tests explicitly for evidence of decorre-
lation of the dust pattern across frequencies. Their Fig. E.1
shows the results for large and small sky patches. The
signal-to-noise ratio is low in clean regions, but no evidence
of decorrelation is found.

As a further check, we artificially suppress the amplitude
of the BK150 × P353 spectra in the Gaussian dust-only
simulations (see Sec. IVA) by a conservative 10% (PIP-
XXX sets a 7% upper limit). We find that the maximum
likelihood value for *r* shifts up by an average of 0.018,
while *A*<sub>d</sub> shifts down by an average of 0.43 µK<sup>2</sup>, with the
size of the shift proportional to the magnitude of the dust
power in each given realization. This behavior is readily
understandable—since the BK150 × BK150 and BK150 ×
P353 spectra dominate the statistical weight, a decrease of
the latter is interpreted as a reduction in dust power, which
is compensated by an increase in *r*. The bias on *r* will be
linearly related to the assumed decorrelation factor.

B. Calibration, analysis, etc.

Figure 3 shows that the *EE* spectrum BK150 × BK150 is
extremely similar to that for BK150 × P143. We can

compare such spectra to set limits on possible decorrelation
between the BICEP2/Keck and Planck maps arising from
any instrumental or analysis related effect, including differ-
ential pointing, polarization angle mischaracterization, etc.
Taking the ratio of BK150 × P143 to the geometric mean of
BK150 × BK150 and P143H1 × P143H2, we find that for
TT the decorrelation is approximately 0.1%. For EE the
signal-to-noise ratio is lower, but decorrelation is limited to
below 2%, and consistent with zero when compared to the
fluctuation of signal + noise simulations.

VI. CONCLUSIONS

BK-I reported a highly significant detection of *B*-mode
polarization, at 150 GHz, in excess of the lensed-ΛCDM
expectation over the range 30 < *ℓ* < 150. This excess has
been confirmed by additional data on the same field from
the successor experiment *Keck Array*. PIP-XXX found that
the level of dust power in a field centered on the
BICEP2/Keck region (but somewhat larger than it) is of
the same magnitude as the reported excess, but noted that,
“the present uncertainties are large,” and that a joint
analysis was required.

In this Letter we have performed this joint analysis, using
the combined BICEP2/Keck maps. Cross-correlating
these maps against all of the polarized frequency bands
of Planck we find a highly significant B-mode detection
only in the cross spectrum with 353 GHz. We emphasize
that this 150 × 353 GHz cross spectrum has a much higher
signal-to-noise ratio than the 353 GHz single-frequency
spectrum that PIP-XXX analyzed.

We have analyzed the data using a multifrequency,
multicomponent fit. In this fit it is necessary to impose a
prior on the variation of the brightness of the polarized dust
emission with observing frequency, since the available data
are unable to constrain this alone, due to the relatively low
signal-to-noise ratio in B-mode polarization at 353 GHz.
However, based on the available information from Planck
on the frequency dependence of polarized dust emission
across the mid- and high-Galactic latitude sky, and the
patch-to-patch stability thereof, this prior appears to be
justified and conservative.

We have shown that the final constraint on the tensor-to-
scalar ratio *r* is very stable when varying the frequency
bands used, as well as the model priors. The result does
differ when using the BICEP2 and *Keck Array* data alone
rather than in combination, but the difference is compatible
with noise fluctuation. Expanding the model to include
synchrotron emission, while also including lower *Planck*
frequencies, does not change the result.

Allowing the amplitude of lensing to be free, we obtain
*A*<sub>L</sub> = 1.13 ± 0.18, with a significance of detection of 7.0σ.
This is the most significant direct detection to date of
lensing in *B*-mode polarization, even compared to experi-
ments with higher angular resolution. The POLARBEAR
experiment has reported a detection of *B*-mode lensing
on smaller angular scales ($500 < l' < 2100$), rejecting the $A_L = 0$ hypothesis at 97.2% confidence [46]. Additionally, ACT [47] and SPT [48] have reported lensing detections in polarization in cross-correlation with some other tracer of the dark matter distribution on the sky.

We have validated the main likelihood analysis on simulations of a dust-only model and performed a simple subtraction of scaled spectra, which approximates a map-based dust cleaning (obtaining an *r* constraint curve that peaks somewhat lower). Finally, we investigated the possibility of astrophysical or instrumental decorrelation of the sky patterns between experiments or frequencies and find no evidence for relevant bias.

The final result is expressed as a likelihood curve for *r*, and yields an upper limit *r* < 0.12 at 95% confidence. The median limit in the lensed-ΛCDM + noise + dust simulations is *r* < 0.075. It is interesting to compare this latter to dust-free simulations using only BICEP2/Keck where the median limit is *r* < 0.03—the difference represents the limitation due to noise in the Planck maps, when marginalizing over dust. The *r* constraint curve peaks at *r* = 0.05 but disfavors zero only by a factor of 2.5. This is expected by chance 8% of the time, as confirmed in simulations of a dust-only model. We emphasize that this significance is too low to be interpreted as a detection of primordial *B* modes. Transforming the Planck temperature-only 95% confidence limit of *r*<sub>0.002</sub> < 0.11 [3] to the pivot scale used in this Letter yields *r*<sub>0.05</sub> < 0.12, compatible with the present result.

A COSMOMC module containing the band powers for all cross spectra between the combined BICEP2/Keck maps and all of the frequencies of Planck is available for download in Ref. [49].

In order to further constrain or detect IGW, additional data are required. The Planck Collaboration may be able to make progress alone using the large angular scale “reionization bump,” if systematics can be appropriately controlled [50]. To take small patch “recombination bump” studies of type pursued here to the next level, data with signal-to-noise comparable to that achieved by BICEP2/Keck at 150 GHz are required at more than one frequency. Figure 13 summarizes the situation. The BICEP2/Keck noise is much lower in the BICEP2/Keck field than the Planck noise. However, since dust emission is dramatically brighter at 353 GHz, it is detected in the cross spectrum between BICEP2/Keck and Planck 353 GHz. Synchrotron is not detected and the crossover frequency with dust is ≲100 GHz. Planck’s PR2 data release [51] shows that for the cleanest 73% of the sky, at 40 arc min scales, the polarized foreground minimum is at ~80–90 GHz. During the 2014 season, two of the Keck Array receivers observed in the 95 GHz band and these data are under active analysis. BICEP3 will add substantial additional sensitivity at 95 GHz in the 2015, and especially 2016, seasons. Meanwhile, many other ground-based and

<figure><img src="image_14.png" /></figure>

FIG. 13 (color). Expectation values, and uncertainties thereon, for the $\ell \sim 80$ BB band power in the BICEP2/Keck field. The green and magenta lines correspond to the expected signal power of lensed ΛCDM and $r = 0.05$. Since CMB units are used, the levels corresponding to these are flat with frequency. The gray band shows the best fit dust model (see Sec. III B) and the blue shaded region shows the allowed region for synchrotron (see Sec. III C). The BICEP2/Keck noise uncertainty is shown as a single starred point, and the noise uncertainties of the Planck single-frequency spectra evaluated in the BICEP2/Keck field are shown in red. The blue points show the noise uncertainty of the cross spectra taken between BICEP2/Keck and, from left to right, Planck 30, 44, 70, 100, 143, 217, and 353 GHz, and plotted at horizontal positions such that they can be compared vertically with the dust and sync curves.

suborbital experiments are making measurements at a variety of frequencies and sky coverage fractions.

ACKNOWLEDGMENTS

BICEP2 was supported by the U.S. National Science Foundation under Grants No. ANT-0742818 and No. ANT-1044978 (Caltech and Harvard) and No. ANT-0742592 and No. ANT-1110087 (Chicago and Minnesota). The development of antenna-coupled detector technology was supported by the JPL Research and Technology Development Fund and Grants No. 06-ARPA206-0040 and No. 10-SAT10-0017 from the NASA APRA and SAT programs. The Keck Array project was supported by the National Science Foundation under Grants No. ANT-1145172 (Harvard), No. ANT-1145143 (Minnesota) and No. ANT-1145248 (Stanford), and from the Keck Foundation (Caltech). We thank the staff of the U.S. Antarctic Program and in particular the South Pole Station without whose help this research would not have been possible. The Planck Collaboration acknowledges the support of the following: ESA; CNES, and CNRS/INSU-IN2P3-INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MINECO,
JA and RES (Spain); Tekes, AoF, and CSC (Finland); DLR and MPG (Germany); CSA (Canada); DTU Space (Denmark); SER/SSO (Switzerland); RCN (Norway); SFI (Ireland); FCT/MCTES (Portugal); ERC and PRACE (EU). A description of the Planck Collaboration and a list of its members, indicating which technical or scientific activities they have been involved in, can be found in Ref. [52].

[1] A. A. Penzias and R. W. Wilson, A measurement of excess antenna temperature at 4080 Mc/s, Astrophys. J. **142**, 419 (1965).

[2] C. L. Bennett, D. Larson, J. L. Weiland, N. Jarosik, G. Hinshaw, N. Odegard, K. M. Smith, R. S. Hill, B. Gold, M. Halpern, E. Komatsu, M. R. Nolta, L. Page, D. N. Spergel, E. Wollack, J. Dunkley, A. Kogut, M. Limon, S. S. Meyer, G. S. Tucker, and E. L. Wright, Nine-year Wilkinson Microwave Anisotropy Probe (WMAP) observations: Final maps and results, Astrophys. J. Suppl. Ser. **208**, 20 (2013).

[3] Planck Collaboration XVI, Cosmological parameters, Astron. Astrophys. **571**, A16 (2014).

[4] Planck Collaboration XXII, Constraints on inflation, Astron. Astrophys. **571**, A22 (2014).

[5] A. A. Starobinsky, Zh. Eksp. Teor. Fiz. Pisma Red. **30**, 719 (1979).

[6] V. A. Rubakov, M. V. Sazhin, and A. V. Veryaskin, Graviton creation in the inflationary universe and the grand unification scale, Phys. Lett. B **115**, 189 (1982).

[7] R. Fabbri and M. D. Pollock, The effect of primordially produced gravitons upon the anisotropy of the cosmological microwave background radiation, Phys. Lett. B **125**, 445 (1983).

[8] L. F. Abbott and M. B. Wise, Constraints on generalized inflationary cosmologies, Nucl. Phys. B**244**, 541 (1984).

[9] A. G. Polnarev, Sov. Astron. **29**, 607 (1985).

[10] U. Seljak and M. Zaldarriaga, Signature of Gravity Waves in Polarization of the Microwave Background, Phys. Rev. Lett. **78**, 2054 (1997).

[11] U. Seljak, Measuring polarization in cosmic microwave background, Astrophys. J. **482**, 6 (1997).

[12] M. Kamionkowski, A. Kosowsky, and A. Stebbins, A Probe of Primordial Gravity Waves and Vorticity, Phys. Rev. Lett. **78**, 2058 (1997).

[13] Planck (http://www.esa.int/planck) is a project of the European Space Agency (ESA) with instruments provided by two scientific consortia funded by ESA member states (in particular the lead countries, France and Italy) with contributions from NASA (USA), and telescope reflectors provided in a collaboration between ESA and a scientific consortium led and funded by Denmark.

[14] Planck Collaboration 2011A, The Planck mission, Astron. Astrophys. **536**, A1 (2011).

[15] Planck Collaboration I, Overview of products and scientific results, Astron. Astrophys. **571**, A1 (2014).

[16] R. H. Hildebrand, J. L. Dotson, C. D. Dowell, D. A. Schleuning, and J. E. Vaillancourt, The far-infrared polarization spectrum: First results and analysis, Astrophys. J. **516**, 834 (1999).

[17] B. T. Draine, Astrophysics of dust in cold clouds, in *The Cold Universe, Saas-Fee Advanced Course 32*, Lecture Notes 2002 of the Swiss Society for Astronomy and Astrophysics (SSAA), edited by A. W. Blain, F. Combes, B. T. Draine, D. Pfenniger, and Y. Revaz (Springer-Verlag, Berlin, 2004), p. 213, astro-ph/0304488.

[18] P. G. Martin, On predicting the polarization of low frequency emission by diffuse interstellar dust, in *EAS Publications Series*, edited by M.-A. Miville-Deschênes and F. Boulanger (EAS Publications, Paris, 2007), Vol. 23, pp. 165–188.

[19] A. Benoît, P. Ade, A. Amblard, R. Ansari, É. Aubourg, S. Bargot, J. G. Bartlett, J.-P. Bernard, R. S. Bhatia, A. Blanchard, J. J. Bock, A. Boscaleri, F. R. Bouchet, A. Bourrachot, P. Camus, F. Couchot, P. de Bernardis, J. Delabrouille, F.-X. Désert, O. Doré et al., First detection of polarization of the submillimetre diffuse galactic dust emission by Archeops, Astron. Astrophys. **424**, 571 (2004).

[20] L. Page, G. Hinshaw, E. Komatsu, M. R. Nolta, D. N. Spergel, C. L. Bennett, C. Barnes, R. Bean, O. Doré, J. Dunkley, M. Halpern, R. S. Hill, N. Jarosik, A. Kogut, M. Limon, S. S. Meyer, N. Odegard, H. V. Peiris, G. S. Tucker, L. Verde et al., Three year Wilkinson Microwave Anisotropy Probe (WMAP) observations: Polarization analysis, Astrophys. J. Suppl. Ser. **170**, 335 (2007).

[21] BICEP2 Collaboration II, Experiment and three-year data set, Astrophys. J. **792**, 62 (2014).

[22] BICEP2 Collaboration I, Detection of B-Mode Polarization at Degree Angular Scales, Phys. Rev. Lett. **112**, 241101 (2014).

[23] J. Dunkley, A. Amblard, C. Baccigalupi, M. Betoule, D. Chuss, A. Cooray, J. Delabrouille, C. Dickinson, G. Dobler, J. Dotson, H. K. Eriksen, D. Finkbeiner, D. Fixsen, P. Fosalba, A. Fraisse, C. Hirata, A. Kogut, J. Kristiansen, C. Lawrence, A. M. Magalhães et al., Prospects for polarized foreground removal, AIP Conf. Proc. **1141**, 222 (2009).

[24] J. Delabrouille, M. Betoule, J.-B. Melin, M.-A. Miville-Deschênes, J. Gonzalez-Nuevo, M. Le Jeune, G. Castex, G. de Zotti, S. Basak, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, J.-P. Bernard, F. R. Bouchet, D. L. Clements, A. da Silva, C. Dickinson, F. Dodu, K. Dolag et al., The pre-launch Planck Sky Model: A model of sky emission at submillimetre to centimetre wavelengths, Astron. Astrophys. **553**, A96 (2013).

[25] Planck Collaboration Int. XIX, An overview of the polarized thermal emission from Galactic dust, arXiv:1405.0871v1.

[26] Planck Collaboration Int. XXII, Frequency dependence of thermal emission from Galactic dust in intensity and polarization, arXiv:1405.0874.

[27] R. Flauger, J. C. Hill, and D. N. Spergel, Toward an understanding of foreground emission in the BICEP2 region, J. Cosmol. Astropart. Phys. **8** (2014) 039.

[28] M. J. Mortonson and U. Seljak, A joint analysis of Planck and BICEP2 B modes including dust polarization uncertainty, J. Cosmol. Astropart. Phys. **10** (2014) 035.

[29] Planck Collaboration Int. XXX, The angular power spectrum of polarized dust emission at intermediate and high Galactic latitudes, arXiv:1409.5738.
[30] *Keck Array* and BICEP2 Collaborations V, Measurements of B-mode polarization at degree angular scales and 150 GHz by the Keck array, arXiv:1502.00643.

[31] http://archives.esac.esa.int/pla2.

[32] Planck Collaboration 2015-I, Overview of products and scientific results, arXiv:1502.01582.

[33] Planck Collaboration 2015-XII, Simulations (to be published).

[34] Planck Collaboration IV, Low Frequency Instrument beams and window functions, Astron. Astrophys. 571, A4 (2014).

[35] Planck Collaboration VII, HFI time response and beams, Astron. Astrophys. 571, A7 (2014).

[36] http://healpix.sourceforge.net/.

[37] K. M. Górski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and M. Bartelmann, HEALPix: A framework for high-resolution discretization and fast analysis of data distributed on the sphere, Astrophys. J. 622, 759 (2005).

[38] With parameters taken from Planck [3].

[39] M. Tristram, J. F. Macías-Pérez, C. Renault, and D. Santos, Xspect, estimation of the angular power spectrum by computing cross-power spectra with analytical error bars, Mon. Not. R. Astron. Soc. 358, 833 (2005).

[40] M. Preece, Analysis of Cosmic Microwave Background Polarisation, Ph.D. thesis, The University of Manchester, U.K., 2011.

[41] S. Hamimeche and A. Lewis, Likelihood analysis of CMB temperature and polarization power spectra, Phys. Rev. D 77, 103013 (2008).

[42] D. Barkats, R. Aikin, C. Bischoff, I. Buder, J. P. Kaufman, B. G. Keating, J. M. Kovac, M. Su, P. A. R. Ade, J. O. Battle, E. M. Bierman, J. J. Bock, H. C. Chiang, C. D. Dowell, L. Duband, J. Filippini, E. F. Hivon, W. L. Holzapfel, V. V. Hristov, W. C. Jones et al., Degree-scale cosmic microwave background polarization measurements from three years of BICEP1 data, Astrophys. J. 783, 67 (2014).

[43] A. Lewis and S. Bridle, Cosmological parameters from CMB and other data: A Monte Carlo approach, Phys. Rev. D 66, 103511 (2002).

[44] Planck Collaboration Int. XVII, Emission of dust in the diffuse interstellar medium from the far-infrared to microwave frequencies, Astron. Astrophys. 566, A55 (2014).

[45] Note that this is the number evaluated at 353 GHz exactly—the equivalent number as integrated over the Planck 353 GHz passband is 4.5 $\mu$K$^2$ and the mask used in PIP-XXX is somewhat different (larger) than the BICEP2/Keck mask used here.

[46] POLARBEAR Collaboration, A measurement of the cosmic microwave background B-mode polarization power spectrum at sub-degree scales with POLARBEAR, Astrophys. J. 794, 171 (2014).

[47] A. van Engelen, B. D. Sherwin, N. Sehgal, G. E. Addison, R. Allison, N. Battaglia, F. de Bernardis, E. Calabrese, K. Coughlin, D. Crichton, J. R. Bond, R. Datta, R. Dunner, J. Dunkley, E. Grace, M. Gralla, A. Hajian, M. Hasselfield, S. Henderson, J. C. Hill et al., The Atacama Cosmology Telescope: Lensing of CMB temperature and polarization derived from cosmic infrared background cross-correlation, arXiv:1412.0626.

[48] D. Hanson, S. Hoover, A. Crites, P. A. R. Ade, K. A. Aird, J. E. Austermann, J. A. Beall, A. N. Bender, B. A. Benson, L. E. Bleem, J. J. Bock, J. E. Carlstrom, C. L. Chang, H. C. Chiang, H.-M. Cho, A. Conley, T. M. Crawford, T. de Haan, M. A. Dobbs, W. Everett et al., Detection of B-Mode Polarization in the Cosmic Microwave Background with Data from the South Pole Telescope, Phys. Rev. Lett. 111, 141301 (2013).

[49] http://bicepkeck.org.

[50] G. Efstathiou and S. Gratton, B-mode detection with an extended Planck mission, J. Cosmol. Astropart. Phys. 6 (2009) 011.

[51] Planck Collaboration 2015-X, Diffuse component separation: Foreground maps, arXiv:1502.01588.

[52] http://www.cosmos.esa.int/web/planck/planck-collaboration.

P. A. R. Ade,¹ N. Aghanim,² Z. Ahmed,³ R. W. Aikin,⁴ K. D. Alexander,⁵ M. Arnaud,⁶ J. Aumont,² C. Baccigalupi,⁷ A. J. Banday,⁸,⁹ D. Barkats,¹⁰ R. B. Barreiro,¹¹ J. G. Bartlett,¹²,¹³ N. Bartolo,¹⁴,¹⁵ E. Battaner,¹⁶,¹⁷ K. Benabed,¹⁸,¹⁹ A. Benoît,²⁰ A. Benoit-Lévy,²¹,¹⁸,¹⁹ S. J. Benton,²² J.-P. Bernard,⁸,⁹ M. Bersanelli,²³,²⁴ P. Bielewicz,⁸,⁹,⁷ C. A. Bischoff,⁵ J. J. Bock,¹³,⁴ A. Bonaldi,²⁵ L. Bonavera,¹¹ J. R. Bond,²⁶ J. Borrill,²⁷,²⁸ F. R. Bouchet,¹⁸,²⁹ F. Boulanger,² J. A. Brevik,⁴ M. Bucher,¹² I. Buder,⁵ E. Bullock,³⁰ C. Burigana,³¹,³²,³³ R. C. Butler,³¹ V. Buza,⁵ E. Calabrese,³⁴ J.-F. Cardoso,³⁵,¹²,¹⁸ A. Catalano,³⁶,³⁷ A. Challinor,³⁸,³⁹,⁴⁰ R.-R. Chary,⁴¹ H. C. Chiang,⁴²,⁴³ P. R. Christensen,⁴⁴,⁴⁵ L. P. L. Colombo,⁴⁶,¹³ C. Combet,³⁶ J. Connors,⁵ F. Couchot,⁴⁷ A. Coulais,³⁷,† B. P. Crill,¹³,⁴ A. Curto,⁴⁸,¹¹ F. Cuttaia,³¹ L. Danese,⁷ R. D. Davies,²⁵ R. J. Davis,²⁵ P. de Bernardis,⁴⁹ A. de Rosa,³¹ G. de Zotti,⁵⁰,⁷ J. Delabrouille,¹² J.-M. Delouis,¹⁸,¹⁹ F.-X. Désert,⁵¹ C. Dickinson,²⁵ J. M. Diego,¹¹ H. Dole,²,⁵² S. Donzelli,²⁴ O. Doré,¹³,⁴ M. Douspis,² C. D. Dowell,¹³ L. Duband,⁵³ A. Ducout,¹⁸,⁵⁴ J. Dunkley,³⁴ X. Dupac,⁵⁵ C. Dvorkin,⁵ G. Efstathiou,³⁸ F. Elsner,²¹,¹⁸,¹⁹ T. A. Enßlin,⁵⁶ H. K. Eriksen,⁵⁷ E. Falgarone,³⁷ J. P. Filippini,⁴,⁵⁸ F. Finelli,³¹,³³ S. Fliescher,⁵⁹ O. Forni,⁸,⁹ M. Frailis,⁶⁰ A. A. Fraisse,⁴² E. Franceschi,³¹ A. Frejsel,⁴⁴ S. Galeotta,⁶⁰ S. Galli,¹⁸ K. Ganga,¹² T. Ghosh,² M. Giard,⁸,⁹ E. Gjerløw,⁵⁷ S. R. Golwala,⁴ J. González-Nuevo,¹¹,⁷ K. M. Górski,¹³,⁶¹ S. Gratton,³⁹,³⁸ A. Gregorio,⁶²,⁶⁰,⁶³ A. Gruppuso,³¹ J. E. Gudmundsson,⁴² M. Halpern,⁶⁴ F. K. Hansen,⁵⁷ D. Hanson,⁶⁵,¹³,²⁶ D. L. Harrison,³⁸,³⁹ M. Hasselfield,⁶⁴ G. Helou,⁴ S. Henrot-Versillé,⁴⁷ D. Herranz,¹¹ S. R. Hildebrandt,¹³,⁴ G. C. Hilton,⁶⁶ E. Hivon,¹⁸,¹⁹ M. Hobson,⁴⁸ W. A. Holmes,¹³ W. Hovest,⁵⁶ V. V. Hristov,⁴
K. M. Huffenberger,⁶⁷ H. Hui,⁴ G. Hurier,² K. D. Irwin,³,⁶,⁸,⁶⁶ A. H. Jaffe,⁵⁴ T. R. Jaffe,⁸,⁹ J. Jewell,¹³ W. C. Jones,⁴² M. Juvela,⁶⁹ A. Karakci,¹² K. S. Karkare,⁵ J. P. Kaufman,⁷⁰ B. G. Keating,⁷⁰ S. Kefeli,⁴ E. Keihänen,⁶⁹ S. A. Kernasovskiy,³ R. Keskitalo,²⁷ T. S. Kisner,⁷¹ R. Kneissl,⁷²,⁷³ J. Knoche,⁵⁶ L. Knox,⁷⁴ J. M. Kovac,⁵ N. Krachmalnicoff,²³ M. Kunz,⁷⁵,²,⁷⁶ C. L. Kuo,³,⁶⁸ H. Kurki-Suonio,⁶⁹,⁷⁷ G. Lagache,⁷⁸,² A. Lähteenmäki,⁷⁹,⁷⁷ J.-M. Lamarre,³⁷ A. Lasenby,⁴⁸,³⁹ M. Lattanzi,³² C. R. Lawrence,¹³ E. M. Leitch,⁸⁰ R. Leonardi,⁵⁵ F. Levrier,³⁷ A. Lewis,⁸¹ M. Liguori,¹⁴,¹⁵ P. B. Lilje,⁵⁷ M. Linden-Vørnle,⁸² M. López-Caniego,⁵⁵,¹¹ P. M. Lubin,⁸³ M. Lueker,⁴ J. F. Macías-Pérez,³⁶ B. Maffei,²⁵ D. Maino,²³,²⁴ N. Mandolesi,³¹,³² A. Mangilli,²,⁴⁷ M. Maris,⁶⁰ P. G. Martin,²⁶ E. Martínez-González,¹¹ S. Masi,⁴⁹ P. Mason,⁴ S. Matarrese,¹⁴,¹⁵,⁸⁴ K. G. Megerian,¹³ P. R. Meinhold,⁸³ A. Melchiorri,⁴⁹,⁸⁵ L. Mendes,⁵⁵ A. Mennella,²³,²⁴ M. Migliaccio,³⁸,³⁹ S. Mitra,⁸⁶,¹³ M.-A. Miville-Deschênes,²,²⁶ A. Moneti,¹⁸ L. Montier,⁸,⁹ G. Morgante,³¹ D. Mortlock,⁵⁴ A. Moss,⁸⁷ D. Munshi,¹ J. A. Murphy,⁸⁸ P. Naselsky,⁴⁴,⁴⁵ F. Nati,⁴² P. Natoli,³²,⁸⁹,³¹ C. B. Netterfield,⁹⁰ H. T. Nguyen,¹³ H. U. Nørgaard-Nielsen,⁸² F. Noviello,²⁵ D. Novikov,⁹¹ I. Novikov,⁴⁴,⁹¹ R. O’Brien,¹³ R. W. Ogburn IV,³,⁶⁸ A. Orlando,⁷⁰ L. Pagano,⁴⁹,⁸⁵ F. Pajot,² R. Paladini,⁴¹ D. Paoletti,³¹,³³ B. Partridge,⁹² F. Pasian,⁶⁰ G. Patanchon,¹² T. J. Pearson,⁴,⁴¹ O. Perdereau,⁴⁷ L. Perotto,³⁶ V. Pettorino,⁹³ F. Piacentini,⁴⁹ M. Piat,¹² D. Pietrobon,¹³ S. Plaszczynski,⁴⁷ E. Pointecouteau,⁸,⁹ G. Polenta,⁸⁹,⁹⁴ N. Ponthieu,²,⁵¹ G. W. Pratt,⁶ S. Prunet,¹⁸,¹⁹,‡ C. Pryke,⁵⁹,³⁰ J.-L. Puget,² J. P. Rachen,⁹⁵,⁵⁶ W. T. Reach,⁹⁶ R. Rebolo,⁹⁷,⁹⁸,⁹⁹ M. Reinecke,⁵⁶ M. Remazeilles,²⁵,²¹,²² C. Renault,³⁶ A. Renzi,¹⁰⁰,¹⁰¹ S. Richter,⁵ I. Ristorcelli,⁸,⁹ G. Rocha,¹³,⁴ M. Rossetti,²³,²⁴ G. Roudier,¹²,³⁷,¹³ M. Rowan-Robinson,⁵⁴ J. A. Rubiño-Martín,⁹⁷,⁹⁹ B. Rusholme,⁴¹ M. Sandri,³¹ D. Santos,³⁶ M. Savelainen,⁶⁹,⁷⁷ G. Savini,¹⁰² R. Schwarz,⁵⁹ D. Scott,⁶⁴ M. D. Seiffert,¹³,⁴ C. D. Sheehy,⁵⁹,¹⁰³ L. D. Spencer,¹ Z. K. Staniszewski,⁴,¹³ V. Stolyarov,⁴⁸,³⁹,¹⁰⁴ R. Sudiwala,¹ R. Sunyaev,⁵⁶,¹⁰⁵ D. Sutton,³⁸,³⁹ A.-S. Suur-Uski,⁶⁹,⁷⁷ J.-F. Sygnet,¹⁸ J. A. Tauber,¹⁰⁶ G. P. Teply,⁴ L. Terenzi,¹⁰⁷,³¹ K. L. Thompson,³ L. Toffolatti,¹⁰⁸,¹¹,³¹ J. E. Tolan,³ M. Tomasi,²³,²⁴ M. Tristram,⁴⁷ M. Tucci,⁷⁵ A. D. Turner,¹³,⁸⁰ L. Valenziano,³¹ J. Valiviita,⁶⁹,⁷⁷ B. Van Tent,¹⁰⁹ L. Vibert,² P. Vielva,¹¹ A. G. Vieregg,¹⁰³,¹¹⁰ F. Villa,³¹ L. A. Wade,¹³ B. D. Wandelt,¹⁸,¹⁹,⁵⁸ R. Watson,²⁵ A. C. Weber,¹³ I. K. Wehus,¹³ M. White,¹¹¹ S. D. M. White,⁵⁶ J. Willmert,⁵⁹ C. L. Wong,⁵ K. W. Yoon,³,⁶⁸ D. Yvon,¹¹² A. Zacchei,⁶⁰ and A. Zonca⁸³

(BICEP2/Keck and Planck Collaborations)

¹School of Physics and Astronomy, Cardiff University, Queens Buildings, The Parade, Cardiff, CF24 3AA, United Kingdom

²Institut d'Astrophysique Spatiale, CNRS (UMR8617) Université Paris-Sud 11, Bâtiment 121, Orsay, France

³Department of Physics, Stanford University, Stanford, California 94305, USA

⁴California Institute of Technology, Pasadena, California, USA

⁵Harvard-Smithsonian Center for Astrophysics, 60 Garden Street MS 42, Cambridge, Massachusetts 02138, USA

⁶Laboratoire AIM, IRFU/Service d'Astrophysique-CEA/DSM-CNRS-Université Paris Diderot, Bâtiment 709, CEA-Saclay, F-91191 Gif-sur-Yvette Cedex, France

⁷SISSA, Astrophysics Sector, via Bonomea 265, 34136, Trieste, Italy

⁸Université de Toulouse, UPS-OMP, IRAP, F-31028 Toulouse cedex 4, France

⁹CNRS, IRAP, 9 Avenue colonel Roche, BP 44346, F-31028 Toulouse cedex 4, France

¹⁰Joint ALMA Observatory, Vitacura, Santiago, Chile

¹¹Instituto de Física de Cantabria (CSIC-Universidad de Cantabria), Avenida de los Castros s/n, Santander, Spain

¹²APC, AstroParticule et Cosmologie, Université Paris Diderot, CNRS/IN2P3, CEA/Irfu, Observatoire de Paris, Sorbonne Paris Cité, 10, rue Alice Domon et Léonie Duquet, 75205 Paris Cedex 13, France

¹³Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, California, USA

¹⁴Dipartimento di Fisica e Astronomia G. Galilei, Università degli Studi di Padova, via Marzolo 8, 35131 Padova, Italy

¹⁵Istituto Nazionale di Fisica Nucleare, Sezione di Padova, via Marzolo 8, I-35131 Padova, Italy

¹⁶University of Granada, Departamento de Física Teórica y del Cosmos, Facultad de Ciencias, Granada, Spain

¹⁷University of Granada, Instituto Carlos I de Física Teórica y Computacional, Granada, Spain

¹⁸Institut d'Astrophysique de Paris, CNRS (UMR7095), 98 bis Boulevard Arago, F-75014, Paris, France

¹⁹UPMC Université de Paris 06, UMR7095, 98 bis Boulevard Arago, F-75014, Paris, France

²⁰Institut Néel, CNRS, Université Joseph Fourier Grenoble I, 25 rue des Martyrs, Grenoble, France

²¹Department of Physics and Astronomy, University College London, London WC1E 6BT, United Kingdom

²²Department of Physics, University of Toronto, Toronto, Ontario, M5S 1A7, Canada

²³Dipartimento di Fisica, Università degli Studi di Milano, Via Celoria, 16, Milano, Italy

²⁴INAF/IASF Milano, Via E. Bassini 15, Milano, Italy

²⁵Jodrell Bank Centre for Astrophysics, Alan Turing Building, School of Physics and Astronomy, The University of Manchester, Oxford Road, Manchester, M13 9PL, United Kingdom

²⁶CITA, University of Toronto, 60 St. George Street, Toronto, Ontario M5S 3H8, Canada

²⁷Computational Cosmology Center, Lawrence Berkeley National Laboratory, Berkeley, California, USA
<sup>28</sup>Space Sciences Laboratory, University of California, Berkeley, California, USA

<sup>29</sup>Sorbonne Université-UPMC, UMR7095, Institut d'Astrophysique de Paris, 98 bis Boulevard Arago, F-75014, Paris, France

<sup>30</sup>Minnesota Institute for Astrophysics, University of Minnesota, Minneapolis, Minnesota 55455, USA

<sup>31</sup>INAF/IASF Bologna, Via Gobetti 101, Bologna, Italy

<sup>32</sup>Dipartimento di Fisica e Scienze della Terra, Università di Ferrara, Via Saragat 1, 44122 Ferrara, Italy

<sup>33</sup>INFN, Sezione di Bologna, Via Irnerio 46, I-40126, Bologna, Italy

<sup>34</sup>Sub-Department of Astrophysics, University of Oxford, Keble Road, Oxford OX1 3RH, United Kingdom

<sup>35</sup>Laboratoire Traitement et Communication de l'Information, CNRS (UMR 5141) and Télécom ParisTech, 46 rue Barrault F-75634 Paris Cedex 13, France

<sup>36</sup>Laboratoire de Physique Subatomique et Cosmologie, Université Grenoble-Alpes, CNRS/IN2P3, 53, rue des Martyrs, 38026 Grenoble Cedex, France

<sup>37</sup>LERMA, CNRS, Observatoire de Paris, 61 Avenue de l'Observatoire, Paris, France

<sup>38</sup>Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, United Kingdom

<sup>39</sup>Kavli Institute for Cosmology, Madingley Road, Cambridge, CB3 0HA, United Kingdom

<sup>40</sup>Centre for Theoretical Cosmology, DAMTP, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, United Kingdom

<sup>41</sup>Infrared Processing and Analysis Center, California Institute of Technology, Pasadena, California 91125, USA

<sup>42</sup>Department of Physics, Princeton University, Princeton, New Jersey, USA

<sup>43</sup>Astrophysics & Cosmology Research Unit, School of Mathematics, Statistics & Computer Science, University of KwaZulu-Natal, Westville Campus, Private Bag X54001, Durban 4000, South Africa

<sup>44</sup>Niels Bohr Institute, Blegdamsvej 17, Copenhagen, Denmark

<sup>45</sup>Discovery Center, Niels Bohr Institute, Blegdamsvej 17, Copenhagen, Denmark

<sup>46</sup>Department of Physics and Astronomy, Dana and David Dornsife College of Letter, Arts and Sciences, University of Southern California, Los Angeles, California 90089, USA

<sup>47</sup>LAL, Université Paris-Sud, CNRS/IN2P3, Orsay, France

<sup>48</sup>Astrophysics Group, Cavendish Laboratory, University of Cambridge, J J Thomson Avenue, Cambridge CB3 0HE, United Kingdom

<sup>49</sup>Dipartimento di Fisica, Università La Sapienza, Piazzale Aldo Moro 2, Roma, Italy

<sup>50</sup>INAF-Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, Padova, Italy

<sup>51</sup>IPAG: Institut de Planétologie et d'Astrophysique de Grenoble, Université Grenoble Alpes, IPAG, F-38000 Grenoble, France, CNRS, IPAG, F-38000 Grenoble, France

<sup>52</sup>Institut Universitaire de France, 103, bd Saint-Michel, 75005, Paris, France

<sup>53</sup>Service des Basses Températures, Commissariat à l'Energie Atomique, 38054 Grenoble, France

<sup>54</sup>Imperial College London, Astrophysics group, Blackett Laboratory, Prince Consort Road, London, SW7 2AZ, United Kingdom

<sup>55</sup>European Space Agency, ESAC, Planck Science Office, Camino bajo del Castillo, s/n, Urbanización Villafranca del Castillo, Villanueva de la Cañada, Madrid, Spain

<sup>56</sup>Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Strasse 1, 85741 Garching, Germany

<sup>57</sup>Institute of Theoretical Astrophysics, University of Oslo, Blindern, Oslo, Norway

<sup>58</sup>Department of Physics, University of Illinois at Urbana-Champaign, 1110 West Green Street, Urbana, Illinois, USA

<sup>59</sup>School of Physics and Astronomy, University of Minnesota, Minneapolis, Minnesota 55455, USA

<sup>60</sup>INAF-Osservatorio Astronomico di Trieste, Via G.B. Tiepolo 11, Trieste, Italy

<sup>61</sup>Warsaw University Observatory, Aleje Ujazdowskie 4, 00-478 Warszawa, Poland

<sup>62</sup>Dipartimento di Fisica, Università degli Studi di Trieste, via Alfonso Valerio 2, Trieste, Italy

<sup>63</sup>INFN/National Institute for Nuclear Physics, Via Valerio 2, I-34127 Trieste, Italy

<sup>64</sup>Department of Physics & Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, British Columbia, Canada

<sup>65</sup>McGill Physics, Ernest Rutherford Physics Building, McGill University, 3600 rue University, Montréal, Quebec, H3A 2T8, Canada

<sup>66</sup>National Institute of Standards and Technology, Boulder, Colorado 80305, USA

<sup>67</sup>Department of Physics, Florida State University, Keen Physics Building, 77 Chieftan Way, Tallahassee, Florida, USA

<sup>68</sup>Kavli Institute for Particle Astrophysics and Cosmology, SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, California 94025, USA

<sup>69</sup>Department of Physics, Gustaf Hällströmin katu 2a, University of Helsinki, Helsinki, Finland

<sup>70</sup>Department of Physics, University of California at San Diego, La Jolla, California 92093, USA

<sup>71</sup>Lawrence Berkeley National Laboratory, Berkeley, California, USA

<sup>72</sup>European Southern Observatory, ESO Vitacura, Alonso de Cordova 3107, Vitacura, Casilla 19001, Santiago, Chile

<sup>73</sup>Atacama Large Millimeter/submillimeter Array, ALMA Santiago Central Offices, Alonso de Cordova 3107, Vitacura, Casilla 763 0355, Santiago, Chile

<sup>74</sup>Department of Physics, University of California, One Shields Avenue, Davis, California, USA

<sup>75</sup>Département de Physique Théorique, Université de Genève, 24, Quai E. Ansermet, 1211 Genève 4, Switzerland

<sup>76</sup>African Institute for Mathematical Sciences, 6-8 Melrose Road, Muizenberg, Cape Town, South Africa

<sup>77</sup>Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of Helsinki, Helsinki, Finland

<sup>78</sup>Aix Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, 13388, Marseille, France
<sup>79</sup>Aalto University Metsähovi Radio Observatory and Department of Radio Science and Engineering, P.O. Box 13000, FI-00076 AALTO, Finland

<sup>80</sup>University of Chicago, Chicago, Illinois 60637, USA

<sup>81</sup>Department of Physics and Astronomy, University of Sussex, Brighton BN1 9QH, United Kingdom

<sup>82</sup>DTU Space, National Space Institute, Technical University of Denmark, Elektrovej 327, DK-2800 Kongens Lyngby, Denmark

<sup>83</sup>Department of Physics, University of California, Santa Barbara, California, USA

<sup>84</sup>Gran Sasso Science Institute, INFN, viale F. Crispi 7, 67100 L'Aquila, Italy

<sup>85</sup>INFN, Sezione di Roma 1, Università di Roma Sapienza, Piazzale Aldo Moro 2, 00185, Roma, Italy

<sup>86</sup>IUCAA, Post Bag 4, Ganeshkhind, Pune University Campus, Pune 411 007, India

<sup>87</sup>School of Physics and Astronomy, University of Nottingham, Nottingham NG7 2RD, United Kingdom

<sup>88</sup>National University of Ireland, Department of Experimental Physics, Maynooth, County Kildare, Ireland

<sup>89</sup>Agenzia Spaziale Italiana Science Data Center, Via del Politecnico snc, 00133, Roma, Italy

<sup>90</sup>Department of Astronomy and Astrophysics, University of Toronto, 50 Saint George Street, Toronto, Ontario, Canada

<sup>91</sup>Lebedev Physical Institute of the Russian Academy of Sciences, Astro Space Centre, 84/32 Profsoyuznaya st., Moscow, GSP-7, 117997, Russia

<sup>92</sup>Haverford College Astronomy Department, 370 Lancaster Avenue, Haverford, Pennsylvania, USA

<sup>93</sup>HGSFP and University of Heidelberg, Theoretical Physics Department, Philosophenweg 16, 69120, Heidelberg, Germany

<sup>94</sup>INAF-Osservatorio Astronomico di Roma, via di Frascati 33, Monte Porzio Catone, Italy

<sup>95</sup>Department of Astrophysics/IMAPP, Radboud University Nijmegen, P.O. Box 9010, 6500 GL Nijmegen, The Netherlands

<sup>96</sup>Universities Space Research Association, Stratospheric Observatory for Infrared Astronomy, MS 232-11, Moffett Field, California 94035, USA

<sup>97</sup>Instituto de Astrofísica de Canarias, C/Vía Láctea s/n, La Laguna, Tenerife, Spain

<sup>98</sup>Consejo Superior de Investigaciones Científicas (CSIC), Madrid, Spain

<sup>99</sup>Departamento Astrofísica, Universidad de La Laguna (ULL), E-38206 La Laguna, Tenerife, Spain

<sup>100</sup>Dipartimento di Matematica, Università di Roma Tor Vergata, Via della Ricerca Scientifica, 1, Roma, Italy

<sup>101</sup>INFN, Sezione di Roma 2, Università di Roma Tor Vergata, Via della Ricerca Scientifica, 1, Roma, Italy

<sup>102</sup>Optical Science Laboratory, University College London, Gower Street, London, United Kingdom

<sup>103</sup>Kavli Institute for Cosmological Physics, University of Chicago, Chicago, Illinois 60637, USA

<sup>104</sup>Special Astrophysical Observatory, Russian Academy of Sciences, Nizhniy Arkhyz, Zelenchukskiy region, Karachai-Cherkessian Republic, 369167, Russia

<sup>105</sup>Space Research Institute (IKI), Russian Academy of Sciences, Profsoyuznaya Street, 84/32, Moscow, 117997, Russia

<sup>106</sup>European Space Agency, ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands

<sup>107</sup>Facoltà di Ingegneria, Università degli Studi e-Campus, Via Isimbardi 10, Novedrate (CO), 22060, Italy

<sup>108</sup>Departamento de Física, Universidad de Oviedo, Avda. Calvo Sotelo s/n, Oviedo, Spain

<sup>109</sup>Laboratoire de Physique Théorique, Université Paris-Sud 11 & CNRS, Bâtiment 210, 91405 Orsay, France

<sup>110</sup>Department of Physics, Enrico Fermi Institute, University of Chicago, Chicago, Illinois 60637, USA

<sup>111</sup>Department of Physics, University of California, Berkeley, California, USA

<sup>112</sup>DSM/Irfu/SPP, CEA-Saclay, F-91191 Gif-sur-Yvette Cedex, France

<footnote>†Brendan.P.Crill@jpl.nasa.gov</footnote>

<footnote>‡pryke@physics.umn.edu</footnote>