A Measurement of the Degree-scale CMB B-mode Angular Power Spectrum with Polarbear

Polarbear is a CMB experiment installed on the 2.5 m aperture Huan Tran Telescope in 2012 January. The telescope is located at the James Ax Observatory at an elevation of 5190 m in the Atacama Desert in Chile. The Polarbear receiver consists of seven wafers containing a total of 1274 transition-edge sensor bolometers cooled to 0.3 K observing the sky through lenslet-coupled double-slot dipole antennas. The instantaneous field of view of the focal plane is approximately 24 in diameter. More details on the receiver and telescope can be found in Arnold et al. (2012) and Kermish et al. (2012).

In 2014 February, we installed a continuously rotating HWP at the focus of the primary mirror to mitigate low-frequency noise. The HWP consists of an antireflective-coated birefringent single crystal sapphire plate rotating at two revolutions per second. More details of the installation of the HWP in the telescope optical path and its performance in a series of test observations can be found in Takakura et al. (2017), hereafter T17.

The Polarbear observations used in this analysis target a 670 effective square degree patch of sky centered on (R.A., decl.) = (+0^h12^m0^s, −59°18'). This patch has significant overlap with the area mapped by South Pole experiments including BICEP, the Keck Array, and SPTpol. The patch coverage is shown in Figure 1. Regular science observations with the HWP began on 2014 July 25 and continued until 2016 December 30. The data set is broken up into three seasons, as described in Table 1.

Zoom In Zoom Out Reset image size

The scan strategy consists of three sets of scans repeated every sidereal day. We scan for a 4 hr 45 minute block as the patch rises above the horizon, a 3 hr 53 minute block at high elevation as the patch transits, and a 4 hr 45 minute block as the patch sets. The rising and setting scan occur at a boresight elevation of 30° and 352 respectively. The elevation of the high-elevation scan is stepped through ten offsets from 455 to 655 to provide even coverage of the patch with one offset chosen per day. Each observation block is broken into hour-long scans referred to as "constant elevation scans" (CESs), after which the detectors are retuned. More details on the calibration procedures are given in Section 3. The telescope is repointed between four hour long blocks but not between hour long CESs.

During scans, the telescope is scanned at a constant velocity of 04 s⁻¹ with a throw of 20° and 35° on the sky for the low-elevation and high-elevation scans, respectively. This results in approximately 70 subscans (defined as one left-going or right-going motion of the telescope) in each hour long CES.

There are 429 days that have data included in the final data set. As a result of the ten high-elevation scan offsets, a complete map of the field is produced once every ten days. Several of these ten-day data subsets have been combined in post-processing due to low yield or incomplete coverage after data selection. This results in 38 approximately even splits of the data set. These splits form the basis for the cross-spectra used in the power spectrum estimation.

Data are discarded from two blocks of time due to mechanical problems with the weatherproof enclosure of the HWP from 2014 October 17 to 2014 November 7 and the eruption of a nearby volcano on 2015 October 30. Additionally, the thermal source used for relative gain calibration was replaced multiple times, due to mechanical problems. The gain calibration is described in Section 3.4.

When the science patch is not available, the cryogenic refrigeration system is recycled, and we perform calibration measurements. These measurements are described in further detail in Section 3. For most of the data set, the fridge was cycled every 24 sidereal hours; however, starting in October of 2015, the fridge cycle was done every 48 sidereal hours in order to provide time to scan the Northern science patch referred to as RA12 in PB14 and PB17. Those data are not included in this analysis.

We reconstruct the HWP angle and interpolate this to match the bolometer readout time stamps. We see a jitter in the reconstructed HWP angle, on the order of 10⁻⁷ rad² Hz⁻¹, due to a combination of physical angle jitter and reconstruction error, which is expected to be subdominant at approximately 10⁻⁹ rad² Hz⁻¹. We simulate the impact of this jitter in Section 5.2 and find the impact on our results to be small.

The telescope pointing reconstruction is performed in a manner very similar to those in PB14 and PB17. The telescope is used for dedicated raster observations of bright point sources selected from PCCS and ATCA catalogs (Murphy et al. 2010; Planck Collaboration et al. 2014) prior to each science observation. The source selection has been modified from PB17 to better match the azimuth and elevation ranges used during the science scans. The observed position computed based on the telescope's azimuth and elevation encoder values are compared to the expected positions. The resulting azimuth and elevation offsets ΔAz(Az, El), ΔEl(Az, El) are fit to an eight-parameter model including terms for uneven heating of the telescope due to insolation (Matsuda 2017). The fiducial model uses the solar radiation as an input parameter—in contrast with PB17, which used ambient temperature. We find a root-mean-square (rms) pointing model residual of 50''.

We also construct pointing solutions that include the Crab Nebula (Tau A) and Jupiter scans performed for polarization angle and beam calibrations in the same way. However, these data are not used in the fiducial boresight pointing solution, as they are observed in a significantly different range of azimuth and elevation from the science observations. We perform this calibration with several different subsets of pointing observations and parameter combinations. In Section 5.2, we show that the difference between these pointing solutions is negligible for the ℓ range considered in this paper.

Following results reported in PB14 and PB17, the detector beam offsets are derived from array raster scans over Jupiter. We find that the reconstructed offsets are consistent with previous results at the level of several tens of arcseconds. We explicitly cut detectors whose mean fitted beam offset differs from PB17 by more than one arcminute. This cut has a negligible impact on the overall efficiency of our data selection.

Following PB14 and PB17, the instrumental beam and the associated window function B_ℓ are measured using dedicated raster observations of Jupiter. We take Jupiter observations with the HWP rotating nominally at 2 Hz at a scan speed of 02 s⁻¹ on the sky. The HWP synchronous structure produced by the instrumental polarization of the primary and differential transmission and emissivity of the HWP is subtracted by masking off a 25' disk in the map domain centered on Jupiter and fitting a first-order polynomial to the time-dependent amplitude of each HWP harmonic. After the HWP synchronous structure is subtracted the TOD are projected into 05 pixels on the sky. The beam window function is taken to be the average of the azimuthally averaged two-dimensional Fourier transform after dividing out the Jupiter disk for each observation and correcting for the pixel window function analytically. The reconstructed beam window function is shown in Figure 2.

Zoom In Zoom Out Reset image size

We find a beam window function that is similar to, but marginally wider than, PB17 with a median Gaussian 36 full width at half maximum. This difference may be due to the telescope being imperfectly focused because the addition of the sapphire HWP lengthens the optical path between the primary and the secondary mirrors, or it may be due to diffraction off the HWP structure itself. We see no weather dependence or variation between seasons in the fitted beamwidth.

Boresight and detector statistical pointing error adds an additional suppression of high-ℓ power in the maps. The pointing model predicts residual pointing errors of 50'' rms. In Section 4.6, we fit an ℓ dependence to the overall gain amplitude, corresponding to the widening of the effective beam due to pointing error, and find no statistical preference for an additional widening of the effective beam due to pointing error. As a result, we convolve the beam function with a Gaussian corresponding to the best-fit pointing model error from the E-mode spectrum and include the effective beamwidth uncertainty in our multiplicative error estimate.

An additional source of systematic uncertainty in the beam comes from detector crosstalk. We measure the detector beam window function using temperature data but compute polarization power spectra. As described in Crowley et al. (2018), crosstalk acts differently in temperature, and polarization also acts differently in the presence of a continuous HWP, resulting in an effective beam mismatch. This effect is described in detail in Section 5.2. We find it to be subdominant to the statistical uncertainty for all spectra.

We estimate the polarization response to Jupiter and decompose this into scalar, dipole, and quadrupole terms. We subtract the scalar term as temperature-to-polarization leakage; this is described in Section 4.4. In Section 5.2, we show the contamination due to higher-order terms is negligible.

Following PB14 and PB17, we use a three-step gain calibration to turn our TOD into physical temperature units with an additional correction for polarization efficiency. We measure a time-dependent gain calibration between detectors using a chopped thermal source before and after each set of four CESs. The gain is linearly interpolated between these measurements. We then calibrate these measurements to temperature on the sky using the same Jupiter scans conducted to measure the beam window function. Finally, we construct a CMB map and scale the overall amplitude of the polarization E-mode fluctuations to match the Planck 2018 data. Performing the absolute gain calibration with the E-mode spectrum jointly constrains the overall gain and polarization efficiency.

At several points during the observations, the thermal source used in the relative calibration was replaced for mechanical reasons. The conversion from the thermal source amplitude to temperature on the sky is performed separately for each source. PB14 describes a correction in this analysis due to the position of a cold HWP in the receiver. The cold HWP was stepped almost daily during the first half of the first season of observations, but it was never moved during the three seasons comprising this data set. We simulate the systematic error introduced by uncertainty on the measured gain acting on the CMB signal and find the contamination in all polarization spectra to be negligible. This is described in Section 5.2.

The detector time constants are measured by sweeping the frequency of the chopped thermal calibration source from 4 to 44 Hz. The source amplitude versus frequency is fit to a single pole transfer function. These time constants are interpolated linearly between the calibrations and deconvolved in the subsequent TOD processing.

The detector polarization angles and efficiencies are derived from dedicated raster scans of Tau A. The raster scan data is taken at 02 s⁻¹ scan speed on the sky with the HWP rotating nominally. We deconvolve the detector time constants measured from a chopped thermal source calibration immediately before and after the Tau A raster scan. This is necessary because a phase lag between the detector TOD and the reconstructed HWP angle appears as a polarization angle error. After correcting for detector time constants, we find no significant correlation between the measured polarization angle and the precipitable water vapor (PWV) measured at the nearby Atacama Pathfinder Experiment⁴¹ (APEX) site. Without time-constant deconvolution, such a correlation is produced by the dependence of the time constants on the atmospheric loading and therefore the local PWV. The demodulated polarization TOD described in Section 4.1 is fit to a beam-convolved polarization map of Tau A from a measurement taken with the IRAM 30 m telescope (Aumont et al. 2010). We find statistically consistent detector polarization angles in focal plane coordinates compared to PB17 and no significant variation between seasons. We find the difference in measured polarization angle between detectors in a pair to be 905 ± 12. This is consistent with the design value of 90°. In Section 5.2, we simulate the expected systematic contamination due to detector polarization angle uncertainty at this level and find the effect to be negligible.

The measured polarized flux from these scans is used to estimate the polarization efficiency of the telescope and receiver. The polarization efficiency is degraded in three ways in addition to the non-ideality of the cold HWP described in PB14 and PB17.

These three effects are the non-ideality of the rotating HWP, the breaking of the Mizuguchi–Dragone (MD) condition (Mizugutch et al. 1976; Dragone 1978), and nonzero bolometer time constants. The modulation efficiency of the warm HWP is estimated from coherent source lab measurements and design detector bandpasses (Arnold et al. 2012). The polarization efficiency term due to the MD breaking is estimated from a GRASP⁴² physical optics simulation (Matsuda et al. 2018). The detector time constant acts as a time domain low-pass filter on the TOD, which has a response less than unity at the polarization modulation frequency. We find the measured values from the Tau A calibration to be consistent with the predictions—but with a larger statistical error, as expected. The non-ideality of the rotating HWP, cold HWP, and MD breaking polarization efficiency terms are intrinsic to the detectors and telescope geometry, and are corrected by rescaling the TOD and noise weights with the predicted values. The average contributions to the overall polarization efficiency from the rotating HWP, cold HWP, and MD breaking are ε ≈ 0.98, 0.98, and 0.93, respectively. The polarization efficiency due to the detector time constant depends on the numerical value of the time constant, which in turn depends on the detector bias point and optical loading. As a result, this polarization efficiency term is corrected by deconvolving the measured time constant from the chopped thermal source before and after each set of four CESs. The average polarization efficiency contribution from the detector time constants is ε ≈ 0.98.

We self-calibrate the overall polarization angle by fitting ${C}_{{\ell }}^{\mathrm{EB}}=0$ following Keating et al. (2013). The overall polarization efficiency is degenerate with the absolute gain of the polarization maps and is set by matching the Polarbear E-mode spectrum to the Planck 143 GHz E-mode data on our patch. This is described in more detail in Section 4.6.

This analysis uses a pre-processing and demodulation procedure similar to that in T17 and Takakura et al. (2019), hereafter T19. A brief overview is provided here for completeness. For simplicity, terms relating to detector nonlinearity and instrumental polarization are neglected here. The raw bolometer TOD in the telescope frame can be modeled following

$$\begin{eqnarray}\begin{array}{rcl}{d}_{m}(t) & = & I(t)+\varepsilon \mathrm{Re}\{[Q(t)+{iU}(t)]\,{e}^{-i(4\chi +2{\theta }_{\det }})\}\\ & & +\ A(\chi ,t)+{{ \mathcal N }}_{m},\end{array}\end{eqnarray} \tag{ 1 }$$

where χ ≈ ωt is the HWP angle, e^−i4χ ≡ m(χ) describes the polarization modulation of the HWP, ε is a polarization efficiency, θ_det is the detector angle with respect to instrument coordinates, ${{ \mathcal N }}_{m}$ is the detector white noise, and A(χ, t) is a slowly varying HWP-synchronous structure in the TOD.

Prior to demodulation, we subtract the HWP-synchronous structure A(χ, t) using an iterative method similar to Johnson et al. (2007). The HWP angle is reconstructed from the encoder. The HWP synchronous structure is estimated following Kusaka et al. (2014) and is decomposed into harmonics. Gaps in the TOD are filled with the HWP synchronous structure. At each harmonic n ∈ {1, 2, ... 7}, the TOD are bandpass filtered with a bandwidth of 1 Hz and demodulated to form a time-dependent amplitude for that harmonic of the HWP synchronous structure. A linear drift is fit to the amplitude of each harmonic. The sum over n of these templates is subtracted from each TOD, and the process is iterated again. After the HWP structure has been subtracted, the polarization signal is reconstructed by multiplying the data by twice the conjugate of the modulation function 2m*(χ) and low-pass filtering to recover the polarization signal,

$$\begin{eqnarray}&&{d}_{d}(t)=\varepsilon [Q(t)+{iU}(t)]{e}^{-2i{\theta }_{\det }}+{{ \mathcal N }}_{d}.\end{eqnarray} \tag{ 2 }$$

The TOD sampled at 8 Hz are used as the input to all subsequent analysis pipeline steps. Due to the high sample rate of the raw TOD and the number of Fourier Transform operations, including the demodulation in our main simulations would be computationally prohibitive. The window function associated with the 8 Hz sample rate is negligible for our ℓ range. The detector time constant deconvolution is performed on the demodulated polarization TOD.

Note that the noise in the demodulated TOD ${{ \mathcal N }}_{d}$ is complex and both the real and imaginary components have twice the variance of the modulated detector white noise ${{ \mathcal N }}_{m}$ due to the factor of two necessary to recover Q and U correctly. There is no intrinsic difference in the polarization sensitivity compared to the pair differencing case. The noise ${{ \mathcal N }}_{d}$ can be well-described as white noise plus a single low-frequency component common to all detectors; this is described in Section 4.3. The demodulation algorithm assumes perfect separation of intensity and polarization signals in time domain frequency—in other words, the HWP frequency is much higher than the scan speed divided by the beam size. In Section 5.2, we quantify the impact of imperfect separation in the real data and find the effect to be negligible. These demodulated data are used as the input to the subsequent data characterization and mapmaking pipelines.

We record the low-pass filtered HWP structure subtracted TOD (the 0f or intensity component) and the demodulation at twice the HWP frequency (the 2f component) for data characterization and selection.

The data selection framework used in this analysis consists of several rounds of increasingly selective criteria to characterize low-frequency noise and mitigate possible systematic contamination. The stages can be roughly described as cutting out glitches and effects well-localized in time, cutting data based on individual detector noise properties over the course of an observation, cutting data based on array noise properties, and cutting data based on map domain noise.

In all stages, the fundamental unit of data considered is the detector subscan. Data where the telescope is accelerating (turnarounds) are rejected completely. A table of efficiencies is shown in Table 2. The primary contributions to the end-to-end efficiency are the fraction of time the patch is available for observations, instrument downtime due to maintenance and weather, and overall focal plane yield due to nonfunctioning detectors and glitch cuts. A total of 2985 CES observations are used in the final mapmaking from the full observation period from 2014 July 25 until 2016 December 30.

In the first stage, a set of time domain glitch criteria similar to those used in PB14 are applied to the pre-demodulated timestreams in order to find and remove high-frequency features. The TOD are convolved with a kernel designed to pick up sharp temporal spikes in the data while nulling the HWP synchronous structure at multiples of 2 Hz. Subscans where the maximum deviation of the convolved TOD is greater than ten times the standard deviation are discarded. This operation is performed on the full sample rate data to improve sensitivity to fast glitches. The post-demodulation intensity, 2f, and polarization TOD are convolved with an additional series of kernels sensitive to fast glitches in the TOD, and subscans are cut following the same criteria.

After these flagging steps, an analysis is performed to remove the low-frequency array common mode glitches in the polarization data that were shown in T19 to be correlated with polarized emission from clouds. We begin by performing a principal component analysis (PCA) on the detector intensity TOD. Working detectors are selected by removing channels whose intensity timestreams are not dominated by the largest eigenmode of the detector covariance matrix, due to atmospheric fluctuations. The Q + iU timestreams from each detector are rotated into the instrument frame via multiplication with ${e}^{i2{\theta }_{\det }}$. These TOD are coadded to form a full focal plane common mode signal that is rotated again into minimum and maximum variance eigenmodes of the Q × U covariance matrix. Subscans where the ratio of the standard deviations of the two modes is greater than ten are cut. This corresponds to the TOD cloud detection criteria used in T19.

At this stage channels that have low signal-to-noise or bad fits to the chopped thermal calibration source, channels that have no position as measured by the beam calibration, and channels with no polarization angle measured from Tau A are cut.

In the second stage, each demodulated timestream is processed through the first portion of the mapmaking TOD filters described in Section 4.4. The TOD are filtered by a first-order subscan polynomial, ground-fixed structure is subtracted by binning the TOD in azimuth and subtracting the resulting template, and monopole temperature-to-polarization leakage is removed. Turnarounds and subscans flagged by the glitch cuts are filled in with white noise matched to the rms of the surrounding samples. The power spectral density (PSD) of the TOD is then measured and fit to a model consisting of white noise and a 1/f² low-frequency noise term. Detectors with anomalously high white noise floors, high low-frequency noise, or poor fits to the model are discarded. For most individual bolometers, the low-frequency noise is undetectably small after the subscan polynomial filter. The action of the high-pass subscan polynomial filter imposes a floor of several times the scan frequency f_scan ≈ 10 mHz on the knee frequency that can be resolved. Fourier domain lines in the PSD are noted and notched out in subsequent data processing. This notch filter is described in Section 4.4.

In the third stage, the common mode timestream is recomputed using the inverse variance individual detector weights, the additional data cuts, and Fourier notch filters defined in the previous stage. We fit the common mode PSD to a white noise and 1/f^α term for instrument frame Q and U separately. White noise is largely uncorrelated between bolometers, whereas the low-frequency component is coherent across detectors—meaning that averaging detectors results in a higher knee frequency. The power-law index α is allowed to float because the higher knee frequency provides a sufficient lever arm to resolve the low-frequency exponent. Observations where the common mode knee frequency f_knee in telescope frame Q (U) is greater than 150 mHz (100 mHz) or where the common mode noise floor is anomalously high are rejected. We find median knee frequencies of 45 mHz (24 mHz) in the common mode Q (U) PSDs and a median polarization white noise floor of $46\,\mu {\rm{K}}\sqrt{s}$ for Q and U individually in CMB temperature units. This corresponds to an array noise equivalent temperature of ${\mathrm{NET}}_{\mathrm{array}}=23\,\mu {\rm{K}}\sqrt{{\rm{s}}}$, which is comparable to PB14 and PB17. The knee frequency distributions differ because common mode Q and U noise are produced by different physical mechanisms. The Q low-frequency noise can be produced by temperature noise brought into the polarization TOD by temperature-to-polarization leakage subtraction or by spurious gain drifts acting on the HWP 4f structure. In contrast, U is out of phase with the 4f HWP structure and is approximately out of phase with the temperature-to-polarization leakage, meaning that low-frequency noise requires a phase drift of the HWP structure in the TOD produced by effects such as detector time constant drift. It is worth noting that we do not see a correlation between the array knee frequencies and the amplitude of the atmospheric fluctuations in the temperature data.

In the fourth stage, the detector weights computed in the second stage and the common mode PSD defined in the third stage are used to create individual observation maps, following the standard TOD filtering outlined in Section 4.4, with the modification that the telescope frame polarization is treated as a scalar field and is not rotated into R.A. and decl. coordinates. The maps are downsampled to degree pixels and a map χ² is computed by comparing the fluctuation in the data to the expectation from the detector noise weights. Maps with an anomalously high χ² are rejected as a check for low-frequency pathologies that are not readily visible in the common mode PSDs. In practice, this cut strongly overlaps with the common mode PSD criteria.

We define two noise models for use in our simulation pipeline, and show good overall agreement between the two. The noise in the data at low frequencies in the time domain can be described by a single common mode 1/f^α component, meaning higher-order modes in the TOD covariance are negligible. We run the end-to-end analysis pipeline in two configurations: one using a TOD noise model and the other using random sign coaddition of individual CES maps to generate matched "signflip" noise realizations.

In the TOD noise model, we generate white noise on a per-detector basis, assuming no correlations between bolometers and the noise weights derived in Section 4.2. The common mode low-frequency noise is synthesized in telescope coordinates for Q and U separately, and is matched to the amplitude and power-law index α fit from the real data. This mode is then rotated into the polarization angle of each detector and added to the uncorrelated white noise.

In the "signflip" configuration, the noise realizations are generated by randomly assigning a + 1 or −1 factor to each map during the coaddition of CESs. This creates noise realizations with the exact power spectrum and correlation structure of the real data, assuming that the noise fluctuations between CESs are uncorrelated. This assumption can be broken by coherent drift in the amplitude of the ground synchronous structure, as described in Section 5.2. We find this effect to be negligible.

We find good agreement between the "signflip" noise realizations and the TOD noise realizations using our cross-spectrum estimator described in Section 4.5. The noise bias derived from both pipelines is shown in Figure 3. The full coadd and all null test splits described in Section 5.1 agree well, except for one null spectrum. In the "top versus bottom bolometers" case that explicitly splits paired detectors, we observe excess variance and an anticorrelation in map space between the two halves that cannot be reproduced by the simple TOD noise model. Temperature noise aliasing into the polarization frequencies in the TOD can create a similar noise anticorrelation. This is naturally accounted for in the "signflip" noise realizations.

Zoom In Zoom Out Reset image size

We use the random sign coaddition pipeline to generate the noise realizations used in our fiducial null tests and error bar estimation.

In this analysis, we use a MASTER "filter and bin" mapmaker (Hivon et al. 2002) where the TOD is filtered to suppress low-frequency noise and projected onto the sky using inverse-variance noise weights.

The mapmaking pipeline takes the demodulated data described in Section 4.1 and applies an additional stack of time domain filters. The ordering of the filtering is set to avoid bias in the fitted temperature leakage coefficients by spurious modes in the temperature and polarization TOD.

The first filtering operation works in the Fourier domain and low-pass filters signal above 1.2 Hz. A set of notch filters are applied to Fourier domain glitches flagged in the data selection pipeline. A notch width of 10 mHz is used on every bolometer on the focal plane when a high-significance glitch is seen in more than 50 detectors for a given CES. An identical set of notch filters are applied to the simulation data.

The second filtering operation removes a second-order polynomial from each bolometer TOD for the whole CES. This mode is expected to be dominated by thermal fluctuations of the focal plane and cryogenics.

The third filter subtracts a ground-fixed template in I, Q, and U. The filter is constructed by averaging the telescope frame TOD in 144 azimuth bins and subtracting the resulting template from the TOD. This operation is performed independently for each bolometer and CES. Unlike PB14 and PB17, the same template is used for both left-going and right-going subscans. We expect that the ground synchronous structure in our data is due to telescope sidelobes far from the main beam interacting with the surrounding terrain. We conservatively assume no correlation between CESs, and subtract a unique ground-fixed template for each CES even though the ground fixed signal that we observe is generally stable between CESs. In Section 5.2, we place an upper limit on the error introduced by possible time variability in the ground-synchronous structure within each CES. We test varying the bin size and subtracting out a smoothed version of the template, and find no significant difference in the final measured power spectrum.

Once these modes have been projected out of the data, we perform a PCA similar to that in T17 to remove temperature-to-polarization leakage due to detector nonlinearity and instrumental polarization from the off-axis telescope design. We see a weak frequency dependence in the temperature leakage coefficients below the telescope scan frequency, and Fourier domain glitches in the TOD at high frequencies. As a result, to estimate the leakage coefficients, we form a copy of the TOD, apply a low-pass filter at 400 mHz, and subtract a first-order polynomial from each subscan. We then compute a 3 × 3 covariance matrix between the I, Q, and U TOD and average this between subscans. The leakage coefficients are determined using the PCA described in T17 with this covariance matrix. We find the estimated leakage coefficients to be stable over the course of our observations. The temperature leakage is subtracted from the original polarization TOD without the subscan polynomial and low-pass filters. Since the leakage coefficient determination is heavily dominated by atmospheric fluctuations, we do not expect this process to be significantly biased by cosmological signal. In Section 5.2, we simulate the error expected in the B-mode power spectrum due to multiplicative detector nonlinearity and find the effect to be negligible. We do not include this leakage or the PCA filter in our main simulation pipeline.

Following this, a first-order polynomial is subtracted from the Q and U TOD for each subscan in polarization to mitigate low-frequency noise. No further processing is applied to the I TOD, as these data are not used in the subsequent analyses.

Finally, a filter is applied to the data to suppress the low-frequency mode seen in all detectors. A low-pass filtered version of the Q and U array common mode is subtracted from each bolometer TOD in the instrument frame. The spectral shape of the low-pass filter is the inverse of the fit PSD of the stacked timestream derived in the data selection pipeline. The same filter is applied to simulated data. As described in Section 4.3, we do not see significant low-frequency noise beyond a single focal plane common mode.

We project the TOD into 8' pixels on the sky using a Lambert Cylindrical equal area projection. The resulting maps for the real data and a sample noise realization are shown in Figure 4. We achieve an effective map depth of 32 μK arcmin after correcting for the beam and TOD filtering.

Zoom In Zoom Out Reset image size

The power spectrum estimation pipeline closely follows PB14 and "Pipeline A" in PB17, with several minor changes to improve numerical accuracy at low ℓ. A brief overview is provided for completeness.

The data set is grouped into 38 bundles of approximately uniform weight and sky coverage. We form power pseudo-spectra by taking cross-spectra between these bundles to remove the noise bias,

$$\begin{eqnarray}&&{\tilde{C}}^{{XY}}=\displaystyle \frac{1}{{\sum }_{i\ne j}{w}_{i}{w}_{j}}\sum _{i\ne j}{w}_{i}{w}_{j}{{\boldsymbol{m}}}_{i}^{X}{{\boldsymbol{m}}}_{j}^{Y* },\end{eqnarray} \tag{ 3 }$$

where ${\boldsymbol{m}}$ and w are the apodized Fourier transform and weight of each map bundle, respectively. The pseudo-spectrum is averaged into bins of Δℓ = 2. Following PB14, we use a pure B-mode estimator based on Smith (2006) and Smith & Zaldarriaga (2007). The apodization mask used in the pseudo-spectrum is significantly more aggressively smoothed than PB17. We apply an 8° cosine square edge taper and an 8° Hamming window to the pixel weight map. This improves the numerical stability of the reconstructed power spectrum at ℓ ≤ 100. We do not mask point sources in the power spectrum estimation. As described in Section 6.4, we do not see evidence for bright polarized point sources in higher-resolution versions of our maps, and we expect the contamination from unresolved polarized point sources to be negligible for this ℓ range.

The noise pseudo-spectrum ${\tilde{N}}^{{XY}}$ is taken to be the pseudo-spectrum of the sum of the apodized map bundles minus the cross pseudo-spectrum ${\tilde{C}}^{{XY}}$.

The pseudo-spectrum is taken to be a linear function of the true underlying power spectrum on the sky

$$\begin{eqnarray}&&{\tilde{C}}_{{\ell }}=\sum _{{{\ell }}^{{\prime} }}{{\boldsymbol{K}}}_{{\ell }{{\ell }}^{{\prime} }}{C}_{{{\ell }}^{{\prime} }},\end{eqnarray} \tag{ 4 }$$

where the transformation is given by

$$\begin{eqnarray}&&{{\boldsymbol{K}}}_{{\ell }{{\ell }}^{{\prime} }}={{\boldsymbol{M}}}_{{\ell }{{\ell }}^{{\prime} }}{F}_{{{\ell }}^{{\prime} }}{B}_{{{\ell }}^{{\prime} }}^{2}.\end{eqnarray} \tag{ 5 }$$

The mode-mixing matrix ${{\boldsymbol{M}}}_{{\ell }{{\ell }}^{{\prime} }}$ describes the mixing of ℓ modes due to finite sky coverage and is estimated directly by computing the pseudo-spectrum of narrowband noise realizations. Here, B_ℓ is the beam window function defined in Section 3.3. The filter transfer function is found via an iterative procedure, following PB14:

$$\begin{eqnarray}&&{F}_{{\ell }}^{n}={F}_{{\ell }}^{n-1}+\displaystyle \frac{{\tilde{C}}_{{\ell }}-{\sum }_{{{\ell }}^{{\prime} }}{{\boldsymbol{M}}}_{{\ell }{{\ell }}^{{\prime} }}{F}_{{{\ell }}^{{\prime} }}^{n-1}{C}_{{{\ell }}^{{\prime} }}{B}_{{{\ell }}^{{\prime} }}^{2}}{{C}_{{\ell }}{B}_{{\ell }}^{2}},{F}_{{\ell }}^{0}=1.\end{eqnarray} \tag{ 6 }$$

In contrast to PB17, we cut the iterative series off at n = 3 to avoid overfitting fluctuations in the simulated pseudo-spectra. This results in a negligible bias on the reconstructed power spectrum. The filter transfer function is calculated using simulated ΛCDM EE-only and BB-only skies drawn from the Planck 2018 baseline TT, EE, TE + lowE + lensing ΛCDM cosmology (Planck Collaboration et al. 2018b) with 1' pixels. We verify that the numerical value of the filter transfer function does not depend on the underlying cosmology used in the simulations.

We estimate the power spectrum in coarser bins of width Δℓ = 50. The binned estimate for the true power spectrum can be written using binning and interpolation operators ${\boldsymbol{P}}$ and ${\boldsymbol{Q}}$,

$$\begin{eqnarray}&&{\hat{C}}_{b}=\sum _{{b}^{{\prime} }{\ell }}{{\boldsymbol{K}}}_{{{bb}}^{{\prime} }}^{-1}{{\boldsymbol{P}}}_{{b}^{{\prime} }{\ell }}{\tilde{C}}_{{\ell }},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{{\boldsymbol{K}}}_{{{bb}}^{{\prime} }}=\sum _{{\ell }\ {{\ell }}^{{\prime} }}{{\boldsymbol{P}}}_{b{\ell }}{{\boldsymbol{M}}}_{{\ell }{{\ell }}^{{\prime} }}{F}_{{{\ell }}^{{\prime} }}{B}_{{{\ell }}^{{\prime} }}^{2}{{\boldsymbol{Q}}}_{{{\ell }}^{{\prime} }{b}^{{\prime} }}.\end{eqnarray} \tag{ 8 }$$

The dependence of the binned spectrum on the underlying spectrum is given by the band power window functions ${{\boldsymbol{w}}}_{b}{\ell }$,

$$\begin{eqnarray}&&{\hat{C}}_{b}=\sum _{{\ell }}{{\boldsymbol{w}}}_{b{\ell }}{C}_{{\ell }},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{{\boldsymbol{w}}}_{b{\ell }}=\sum _{{b}^{{\prime} }{{\ell }}^{{\prime} }}{{\boldsymbol{K}}}_{{{bb}}^{{\prime} }}^{-1}{{\boldsymbol{P}}}_{{b}^{{\prime} }{{\ell }}^{{\prime} }}{{\boldsymbol{K}}}_{{{\ell }}^{{\prime} }{\ell }}.\end{eqnarray} \tag{ 10 }$$

These band power window functions are shown in Figure 5. The shape of the lowest bandpower is due to sensitivity degradation at low ℓ.

Zoom In Zoom Out Reset image size

We have checked that the flat-sky approximation does not introduce a significant bias in the measured power spectra.

The statistical uncertainty on the reconstructed power spectrum is taken to be the standard deviation of the reconstructed power spectrum of Monte Carlo (MC) simulations containing filtered ΛCDM sky realizations and "signflip" noise realizations. We compute analytic uncertainties in the same way as PB14, and find those results to be consistent with our MC uncertainties.

For all spectra, we follow the D_ℓ ≡ ℓ(ℓ + 1) C_ℓ/(2π) convention unless otherwise noted.

4.5.1. E-mode and B-mode Mixing

The timestream filtering mixes E-mode power into the B-mode spectrum. This is subtracted in the pseudo-spectra following PB14,

$$\begin{eqnarray}&&{\tilde{C}}_{{\ell }}^{E\to B}=\displaystyle \frac{{F}_{{\ell }}^{E\to B}}{{F}_{{\ell }}^{E\to E}}{\tilde{C}}_{{\ell }}^{E}.\end{eqnarray} \tag{ 11 }$$

This reconstructs the correct central value of ${C}_{{\ell }}^{\mathrm{BB}}$; however, it does not remove the excess variance in the B-mode spectrum. We find that the level of B-mode power introduced by timestream filtering is comparable to the expected lensing B-mode signal in our lowest band power. Because the fluctuations due to leaked E-modes are below the fluctuations of the noise bias, the contribution to the statistical uncertainties is small. We nonetheless account for this effect in our statistical uncertainty by subtracting the measured E-mode spectrum in each realization in our simulations. We do not perform any matrix-based separation of E- and B-modes as described in BICEP2 Collaboration et al. (2016).

4.5.2. Quantifying Low-ℓ Statistical Performance

To accurately describe the low-ℓ statistical performance of Polarbear, it is necessary to account for the effect of the TOD filtering and the beam window function in the measured noise bias. We do this following BK15 by referring the noise pseudo-spectrum ${\tilde{N}}_{{\ell }}$ to a true power spectrum on the sky N_b following Equation (7). We write the number of degrees of freedom per band power as an ℓ-dependent effective sky area. These two quantities represent the noise contribution to the power spectrum statistical uncertainties. Figure 6 and Table 3 show the results of this analysis. It should be noted that this is not exactly equivalent to an analytic estimate of the statistical errors, because it does not account for terms arising from the CMB signal variance.

Zoom In Zoom Out Reset image size

Table 3.  Noise Bias and Effective Degrees of Freedom

Nominal Band Definition (ℓ)	Band Power Centroid (ℓ)	Noise Bias Nb (10−4 μK2)	Effective Noise (fsky)
50 < ℓ ≤ 100	81.5	1.665	0.0089
100 < ℓ ≤ 150	125.6	1.006	0.0130
150 < ℓ ≤ 200	175.1	0.875	0.0125
200 < ℓ ≤ 250	224.8	0.872	0.0147
250 < ℓ ≤ 300	274.4	0.872	0.0149
300 < ℓ ≤ 350	324.8	0.891	0.0181
350 < ℓ ≤ 400	374.9	0.921	0.0152
400 < ℓ ≤ 450	425.0	0.943	0.0136
450 < ℓ ≤ 500	474.8	0.963	0.0137
500 < ℓ ≤ 550	524.7	0.998	0.0130
550 < ℓ ≤ 600	574.7	1.030	0.0137

Note. Data points for the curves shown in Figure 6. These data use the "signflip" noise model and the fiducial internal cross power spectrum estimator. These numbers are written in C_ℓ units and do not include the factor of ℓ(ℓ + 1)/2π used in the rest of this work.

Download table as:  ASCIITypeset image

We additionally compute an ℓ knee that can be compared to the numbers presented in The Simons Observatory Collaboration et al. (2019). We fit our estimated power spectrum uncertainties following the same formalism and find a knee of ℓ = 90. This ℓ knee can also be computed by incorporating the effective sky area degradation into the noise bias, ${N}_{b}/\sqrt{{f}_{\mathrm{sky}}}$, using the numbers in Table 3 giving consistent results.

We perform an overall gain calibration by cross-correlating our data to the Planck satellite. We use the Planck 2018 PR3 143 GHz full mission maps⁴³ and process them through our filtering and mapmaking pipeline. We compute debiased spectra for both the Polarbear internal cross-spectra using the fiducial power spectrum estimate and the fully coadded Polarbear maps correlated with the scanned Planck maps using the full coadd power spectrum estimator described in Appendix. We fit an overall gain calibration factor based on the ratio of the E-mode spectra,

$$\begin{eqnarray}&&{\hat{g}}_{b}={\hat{C}}_{b}^{\mathrm{EE},\mathrm{PB}}/\,{\hat{C}}_{b}^{\mathrm{EE},\mathrm{PB}\times \mathrm{Planck}}.\end{eqnarray} \tag{ 12 }$$

The uncertainty on this ratio is computed via MC simulation, holding the underlying sky realization fixed. We use 96 realizations of the Planck FFP10 noise model (Planck Collaboration et al. 2018a) to approximate the Planck map noise. We fit this calibration to an ℓ-dependent gain model, accounting for a smearing of the beam profile following:

$$\begin{eqnarray}&&g({\ell })={g}_{0}\,\exp \left[-\displaystyle \frac{{\ell }({\ell }\ +1)}{2}\,{\sigma }^{2}\right],\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\hat{g}}_{b}=\sum _{{\ell }}{{\boldsymbol{w}}}_{b{\ell }}{g}_{{\ell }}.\end{eqnarray} \tag{ 14 }$$

We find an overall gain calibration factor of 1.08 ± 0.04 in amplitude and an effective beam smearing of σ² = 1.14 ± 5.58 arcmin². We convolve the beam with the best-fit value for the pointing model error and treat the uncertainties in g₀ and σ² as a calibration error term in our final results. We compute an alternate absolute gain calibration fitting the Polarbear spectrum to the best-fit Planck ΛCDM theory spectrum and find agreement with the fiducial calibration at the percent level across our ℓ range.

After applying this absolute calibration, we compare the E-mode spectrum to Planck and form a null spectrum:

$$\begin{eqnarray}&&{\hat{C}}_{b,\mathrm{null}}={\hat{C}}_{b}^{\mathrm{EE},\mathrm{PB}}-{\hat{C}}_{b}^{\mathrm{EE},\mathrm{PB}\times \mathrm{Planck}}.\end{eqnarray} \tag{ 15 }$$

The uncertainty on this null spectrum is computed via MC. We find this null spectrum to be consistent with zero with χ²/ν = 7.0/9, corresponding to a probability to exceed (PTE) of 64%. When we compare our measured E-mode spectrum to the best-fit ΛCDM theory, we observe a marginally significant discrepancy in our highest two ℓ bins. This appears to be due to an anisotropic feature seen in the two-dimensional power spectrum at approximately the size scale and orientation of the detector wafers. These fluctuations have no significant counterpart in the null tests described in Section 5.1. These fluctuations do not depend on any of the operations in the TOD filtering and mapmaking pipeline that have characteristic scales on the sky or frequencies in the time domain. The overall gain and beam calibration does not significantly shift if these two bins are removed from the analysis. We find that there is no significant shift in the B-mode spectrum when these regions of the Fourier plane are masked in the pseudo-spectrum estimation. We show our measured E-mode spectrum as well as the residuals compared to Planck and the ΛCDM theory in Figure 7.

Zoom In Zoom Out Reset image size

After applying the overall gain calibration, we self-calibrate the overall instrument polarization angle following Keating et al. (2013). We apply an overall polarization angle correction Δψ, such that the measured ${C}_{{\ell }}^{\mathrm{EB}}$ power is minimized:

$$\begin{eqnarray}&&-2\mathrm{ln}({ \mathcal L })=\sum _{b}{\left[\displaystyle \frac{{\hat{C}}_{b}^{\mathrm{EB}}-(1/2)\sin (4{\rm{\Delta }}\psi ){\sum }_{{\ell }}{{\boldsymbol{w}}}_{b{\ell }}{C}_{{\ell }}^{\mathrm{EE}}}{{\rm{\Delta }}{\hat{C}}_{b}^{\mathrm{EB}}}\right]}^{2}.\end{eqnarray} \tag{ 16 }$$

Here, ${\rm{\Delta }}{\hat{C}}_{b}^{\mathrm{EB}}$ is the uncertainty in the reconstructed EB spectrum from the MC simulations. We find an overall calibration Δψ = −061 ± 022. After applying this calibration, we find that ${\hat{C}}_{b}^{\mathrm{EB}}$ is consistent with zero, with a total χ² PTE of 77%. Our absolute angle calibration is compatible with PB17, which reported Δψ _PB17 = −079 ± 016. This calibration is shown in Figure 8. It should be noted that this calibration is independent of PB17.

Zoom In Zoom Out Reset image size

The self-calibrated polarization angle correction is applied in the map domain. The absolute gain computed with the EE spectrum differs slightly from the MC simulations used to establish the band power errors and covariances. To account for this in the final power spectrum uncertainties, we assume that the total variance in each band power is the sum of the signal variance and the noise variance as derived from signal-only and noise-only simulations, respectively. The noise variance is rescaled to the best-fit absolute gain calibration.

We quantify a calibration uncertainty in the measured power spectrum using the uncertainty in the overall gain calibration g₀, pointing model error σ², polarization angle self calibration Δψ, and statistical uncertainty on the beam window function ΔB_ℓ. The correlation between g₀ and σ² is modeled as a simple Gaussian covariance. The three effects are added in quadrature to form a total calibration error. This represents an upper limit, as the overall gain calibration and beam window function are degenerate. The numerical values of these calibration errors in the estimated power spectrum are shown in Section 6.1.

We perform a set of null tests to establish the internal consistency of the data set and search for possible systematic contamination in the final power spectra. In general, it is not possible to construct difference maps between halves of the data with zero signal due to anisotropic scanning and filtering effects. This effect becomes particularly pronounced at low ℓ. As a result, we follow the formalism developed originally by the QUIET collaboration (Bischoff 2010) and used in PB14 and PB17. We include the filtering and mode mixing explicitly in the construction of the null spectrum

$$\begin{eqnarray}&&{\hat{C}}_{{\ell }}^{\mathrm{null}}={\hat{C}}_{{\ell }}^{{\rm{A}}}+{\hat{C}}_{{\ell }}^{{\rm{B}}}-2{\hat{C}}_{{\ell }}^{\mathrm{AB}},\end{eqnarray} \tag{ 17 }$$

where ${\hat{C}}_{{\ell }}^{{\rm{A}}}$ and ${\hat{C}}_{{\ell }}^{{\rm{B}}}$ are the debiased autospectrum for each half of the split. Likewise, ${\hat{C}}_{{\ell }}^{\mathrm{AB}}$ is the cross-spectrum between the splits. The spectra are computed using the same cross-spectrum formalism and map bundles used in the main pipeline.

The filter transfer function and mode-mixing matrix are computed in the same way as the fiducial power spectrum pipeline, using 92 EE-only and 92 BB-only Planck 2018 ΛCDM input maps for each test. Only map regions present in both halves of the split are used to form the pseudo-spectrum and mode-mixing matrix. The null spectra are computed using the same ℓ binning as the final power spectrum. The EE, EB, and BB null spectra are compared to 192 EE + BB signal and noise simulations. We find that this number of simulations is adequate for percent-level statistical uncertainties on the null spectrum PTE values and filter transfer functions across our full ℓ range.

For our fiducial null test statistics, we use noise realizations generated with the "signflip" pipeline. There is no significant difference from the PTE values computed with the TOD noise model, with the exception of the "top versus bottom" null test that explicitly separates detector pairs. This is due to the presence of an additional anticorrelated noise term when detector pairs are separated that is not included in the TOD noise model described in Section 4.3.

5.1.1. Null Test Data Splits

5.1.2. Null Test Statistics and Analysis

For each bin in each null spectrum, we compute the statistic ${\chi }_{\mathrm{null}}\equiv {\hat{C}}_{b}^{\mathrm{null}}/\sigma ({\hat{C}}_{b}^{\mathrm{null}})$, where $\sigma ({\hat{C}}_{b}^{\mathrm{null}})$ is the standard deviation of the MC null spectra. We use both χ_null and ${\chi }_{\mathrm{null}}^{2}$ because the former is more sensitive to systematic biases and the latter is more sensitive to outliers. Figure 9 shows the χ_null and ${\chi }_{\mathrm{null}}^{2}$ distributions compared to the expectation from MC simulations. Table 4 shows the PTE of the total ${\chi }_{\mathrm{null}}^{2}$ summed over ℓ bins for each test and summed over tests for each bin.

Zoom In Zoom Out Reset image size

Table 4.  Null Test PTE Values

	Total EEχ2 PTE	Total EBχ2 PTE	Total BBχ2 PTE
Null test summed over ℓ bins
First half versus second half	86.5%	43.2%	79.7%
Rising versus middle and setting	85.4%	70.8%	6.8%
Middle versus rising and setting	63.0%	63.0%	39.1%
Setting versus rising and middle	50.5%	53.1%	19.3%
Left-going versus right-going subscans	36.5%	26.0%	6.2%
High-gain versus low-gain CESs	45.3%	83.3%	3.1%
High PWV versus low PWV	50.5%	33.9%	28.6%
Common mode Q knee frequency	57.8%	36.5%	26.6%
Common mode U knee frequency	91.7%	54.7%	79.7%
Mean temperature leakage by bolometer	78.6%	54.2%	56.8%
2f amplitude by bolometer	5.7%	32.8%	83.9%
4f amplitude by bolometer	35.4%	18.2%	33.3%
Q versus U pixels	72.4%	64.6%	28.1%
Sun above or below the horizon	76.6%	91.1%	32.2%
Moon above or below the horizon	76.0%	77.6%	57.3%
Top half versus bottom half	79.7%	35.4%	27.1%
Left half versus right half	53.6%	33.9%	76.0%
Top versus bottom bolometers	36.5%	55.7%	35.4%
ℓ bin summed over null tests
50 ≤ ℓ ≤ 100	31.8%	58.3%	44.8%
100 < ℓ ≤ 150	64.1%	14.1%	24.5%
150 < ℓ ≤ 200	61.5%	46.9%	96.9%
200 < ℓ ≤ 250	71.4%	74.0%	28.6%
250 < ℓ ≤ 300	83.9%	7.3%	26.6%
300 < ℓ ≤ 350	50.5%	92.7%	6.8%
350 < ℓ ≤ 400	64.1%	97.9%	92.2%
400 < ℓ ≤ 450	44.3%	84.4%	5.2%
450 < ℓ ≤ 500	96.9%	63.5%	3.1%
500 < ℓ ≤ 550	68.8%	84.5%	49.0%
550 < ℓ ≤ 600	49.5%	16.1%	49.5%

Note. PTE values for the total χ² of each null spectrum summed over ℓ bins and each ℓ bin summed over null spectra. None of the null spectra indicate significant problems. PTE values are computed directly from the 192 signal+noise simulations and are therefore quantized at the 0.5% level.

Download table as:  ASCIITypeset image

To probe for systematic contamination and consistency of the noise model, we compute five statistics on the χ_null values. These statistics were determined before computing spectra with the real data.

• 1.  
  "Average χ overall:" the mean value of χ_null for all tests and bins
• 2.  
  "Most extreme χ² by bin:" the most extreme ${\chi }_{\mathrm{null}}^{2}$ when summing spectra over tests
• 3.  
  "Most extreme χ² by test:" the most extreme ${\chi }_{\mathrm{null}}^{2}$ when summing spectra over bins
• 4.  
  "Most extreme χ² overall:" the most extreme ${\chi }_{\mathrm{null}}^{2}$ for all bins and tests
• 5.  
  "Total χ² overall:" the sum of ${\chi }_{\mathrm{null}}^{2}$ for all spectra

For each statistic, we compute a PTE by comparing the real data to same statistic computed with the MC realizations. Directly comparing the data to the simulations accounts for any correlations that exist between bins and tests in the computation of the PTE values. We define an additional statistic P_low, which is the lowest of the five PTE values. We require that the PTE of P_low be greater than 5%. The numerical value for these statistics can be seen in Table 5. All spectra (EE, EB, and BB) pass these criteria.

Table 5.  Null Test PTE Values

Null Statistic	EE PTE	EB PTE	BB PTE
Average χ overall	73.4%	80.7%	9.9%
Most extreme χ2 by bin	96.9%	43.8%	31.7%
Most extreme χ2 by test	70.3%	98.4%	57.3%
Most extreme χ2 overall	48.4%	84.9%	66.1%
Total χ2 overall	90.6%	78.7%	12.5%
Lowest statistic	85.4%	86.4%	33.9%
KS test on all bins	10.1%	60.4%	15.9%
KS test on all spectra	6.4%	31.8%	10.5%
KS test overall	35.9%	27.5%	14.6%

Note. PTE values for each of the high-level null test statistics.

Download table as:  ASCIITypeset image

In addition, we require that the PTE of the ${\chi }_{\mathrm{null}}^{2}$ values by test, by bin, and overall be consistent with a uniform distribution, using a Kolmogorov–Smirnov (KS) test to check for systematic mismatch between the real and simulated uncertainties. Figure 10 shows the PTE distribution for all bins and tests. We find the PTE distributions to be consistent with uniform.

Zoom In Zoom Out Reset image size

In addition to the null tests, we simulate several known sources of systematic error. We find upper bounds on the contamination introduced by these effects to be subdominant to the statistical error in all spectra.

The systematic error estimate uses a modified version of the simulation pipeline developed for the null tests and the power spectrum estimation. For most systematic effects, a signal-only ΛCDM sky is scanned to form simulated TOD. Next, these TOD are distorted following a model of the given effect, and then filtered and projected onto the sky using the fiducial mapmaker. We compute pseudo-spectra from these distorted maps and refer them to the underlying sky using the same mode coupling and filter transfer functions used for the fiducial power spectrum estimate. Since several systematics are suppressed by the overall gain and polarization angle calibration, we perform these overall calibrations in each simulation. We do not model the ℓ-dependent gain model applied to the real data, and only fit the overall gain g₀. In parallel, the same mapmaking and power spectrum estimation is performed without distorting the TOD to form a reference simulation. The systematic error is taken to be the absolute difference between the contaminated and reference power spectra.

To form an overall systematic error estimate, we linearly add the power spectrum contamination from each set of simulations. In practice, we expect that each source of error will be largely uncorrelated with the others, meaning this is a conservative upper limit. The total systematic error estimate as well as the contributions from groups of effects in EE and BB are shown in Figure 11. We find the expected systematic contamination in the EB spectrum to be likewise subdominant to our statistical errors.

Zoom In Zoom Out Reset image size

5.2.1. Gain Miscalibration, Time Constant Drift, Detector Nonlinearity

5.2.2. Polarization Angle Error

5.2.3. Boresight Pointing Error

5.2.4. Crosstalk

We observe an electrical coupling between detector readouts on the same cables. We assume that this crosstalk is linear and constant through the entire data set. We estimate the amplitude of this effect using the observations of Jupiter used for the gain and beam calibration. Crosstalk appears in these data as an apparent negative copy of the beam shape several tens of arcminutes away from the main beam. We estimate the amplitude of these images using a matched filter and construct a matrix representing the coupling of signal in detector i to observed signal in detector j,

$$\begin{eqnarray}&&{d}_{j,\mathrm{observed}}={d}_{j,\mathrm{real}}+\sum _{i\ne j}{{\boldsymbol{L}}}_{{ij}}{d}_{i,\mathrm{real}}.\end{eqnarray} \tag{ 18 }$$

We find the matrix ${\boldsymbol{L}}$ to be sparse, and do not see significant crosstalk for detector readout on different cables. The median nonzero off-diagonal element of the matrix ${\boldsymbol{L}}$ is 1%, and the median row sum representing the ratio of crosstalked power to power in the main beam is approximately 4%.

We estimate this systematic error as the sum of two effects. We first quantify error introduced in the beam calibration. This is due to the fact that crosstalk is strongly suppressed in polarization because of the different polarization angles of each detector (Crowley et al. 2018) in HWP experiments. As a result, the temperature and polarization see different effective beam profiles. The second effect is the direct distortion of the polarization signal, which we simulate by injecting crosstalk at the TOD level in signal-only simulations. We find the beam miscalibration effect to be the dominant systematic in our E-mode spectrum. It is possible to strongly suppress crosstalk by inverting the mixing matrix ${\boldsymbol{L}}$ in the data processing pipeline as done in Henning et al. (2018). We do not take this step, however, as the expected contamination is already below our statistical error.

5.2.5. Ground Pickup

5.2.6. HWP Signal Aliasing

5.2.7. HWP Imperfections

5.2.8. Temperature-to-polarization Leakage

5.2.9. Cross Polarization from MD Breaking

Our estimated B-mode power spectrum is shown in Figure 12. We observe a modest excess above the ΛCDM lensing expectation in our lowest two ℓ bins.

Zoom In Zoom Out Reset image size

The band power uncertainties from all terms (statistical, systematic, and calibration) are shown in Table 6. We find the statistical error to be dominant in all ℓ bins.

Table 6.  Band Powers and Uncertainties

Nominal Band Definition	${D}_{{\ell }}^{\mathrm{BB}}$ (μK2)	${D}_{{\ell }}^{\mathrm{BB}}$ stat error (μK2)	${D}_{{\ell }}^{\mathrm{BB}}$ syst error (μK2)	${D}_{{\ell }}^{\mathrm{BB}}$ cal error (μK2)
50 < ℓ ≤ 100	0.0390	0.0268	0.0040	0.0027
100 < ℓ ≤ 150	0.0449	0.0288	0.0010	0.0023
150 < ℓ ≤ 200	0.0194	0.0415	0.0012	0.0009
200 < ℓ ≤ 250	0.0345	0.0559	0.0007	0.0014
250 < ℓ ≤ 300	−0.0566	0.0747	0.0015	0.0019
300 < ℓ ≤ 350	0.0471	0.0910	0.0022	0.0019
350 < ℓ ≤ 400	0.3731	0.1236	0.0022	0.0168
400 < ℓ ≤ 450	0.0503	0.1641	0.0021	0.0031
450 < ℓ ≤ 500	−0.0189	0.1907	0.0029	0.0017
500 < ℓ ≤ 550	0.0143	0.2377	0.0019	0.0015
550 < ℓ ≤ 600	−0.1037	0.2765	0.0049	0.0126

Notes. Measured bandpowers as well as the statistical, combined systematic, and calibration uncertainty error estimates. The calibration and systematic error estimates are described in Sections 4.6 and 5.2, respectively.

Download table as:  ASCIITypeset image

We compute an estimate of the overall amplitude of our observed B-mode signal relative to previous measurements. For the purposes of this calculation, we assume that the underlying sky consists of a lensing CMB component corresponding to the Planck 2018 ΛCDM lensing B-mode spectrum and a foreground component modeled by a power law ${D}_{{\ell },\mathrm{dust}}=9\,\times {10}^{-3}{({\ell }/80)}^{-0.6}\,\mu {{\rm{K}}}^{2}$ taken from the BK15 spectral decomposition at 150 GHz. We find a reduced χ² of 11.6/11 compared to this model indicating good agreement. Naively fitting for an overall B-mode amplitude to rescale this template, we find A_BB = 1.8 ± 0.8, disfavoring the null BB hypothesis at 2.2σ. This estimate neglects the slightly non-Gaussian shape of the band power distributions. However, that is accounted for in our cosmological parameter constraints, as shown in subsequent sections.

We additionally compute a naive uncertainty on r from our ΛCDM lensing-only (r = 0) simulations. We find σ(r) = 0.34, neglecting the non-Gaussian shape of the likelihood.

We compute cross-spectra with three Planck 2018 High Frequency Instrument (HFI) maps to quantify the contribution of galactic dust to our B-mode autospectrum. We form ten unique auto- and cross-spectra between four maps:

• 1.  
  Polarbear 150 GHz map;
• 2.  
  Planck 143 GHz frequency map;
• 3.  
  Planck 217 GHz frequency map;
• 4.  
  Planck 353 GHz frequency map.

We process the Planck PR3 full mission frequency maps through the Polarbear observing pipeline to create Planck maps as seen by Polarbear. Additionally, we process 96 noise realizations from the Planck FFP10 noise simulations for each frequency, to establish the Planck noise bias. The Polarbear noise bias is estimated using 192 "signflip" noise realizations.

The spectra involving the Planck maps are formed using the alternate full autospectrum pipeline described in the Appendix, where we simply compute the power spectrum of the fully coadded map rather than the cross-spectrum of the ten-day data subsets. We compute the measured debiased power spectrum ${\hat{C}}_{b}$ for display purposes by subtracting the mean autospectrum of the FFP10 noise simulations from the autospectrum of the real map. This is done to avoid relying on same-frequency cross-spectra in the Planck 2018 PR3 maps, which are known to be contaminated by systematics at the lowest ℓ values (Planck Collaboration et al. 2018a).

The Polarbear autospectrum and noise bias are computed using the fiducial cross-spectrum pipeline. Using the autospectrum of the fully coadded map in place of the internal cross-spectrum formalism does not significantly shift the final parameter estimate. However, we use the fiducial power spectrum pipeline for consistency with previous results, as well as to bolster robustness against noise-bias misestimation.

The six cross-spectra between frequency bands are formed directly from the cross-spectra of the fully coadded maps, assuming that systematics and noise are uncorrelated between frequencies and experiments.

We use a quasi-analytic approach to estimate uncertainty in the cross- and autospectra, motivated by the computational cost of running the number of simulations necessary to directly estimate the full covariance matrix. We assume that every spectrum follows a reduced χ-squared distribution, with the same effective number of degrees of freedom per band power (denoted ν_b) determined by the patch geometry and TOD filtering. We estimate the number of degrees of freedom from the fractional uncertainty in the autospectrum of the scanned noise realizations,

$$\begin{eqnarray}&&{\nu }_{b}=2{\left(\displaystyle \frac{{\hat{C}}_{b}}{\sigma ({\hat{C}}_{b})}\right)}^{2}.\end{eqnarray} \tag{ 19 }$$

For the autospectra, we estimate the uncertainty of the measured power spectrum ${\hat{C}}_{b}$ using the best-fit signal power spectrum C_b and mean noise bias N_b,

$$\begin{eqnarray}&&{\rm{\Delta }}{\hat{C}}_{b}=\sqrt{\displaystyle \frac{2}{{\nu }_{b}}}\,({C}_{b}+{N}_{b}).\end{eqnarray} \tag{ 20 }$$

Equivalently, the uncertainty in the measured cross-spectrum can be written in the form:

$$\begin{eqnarray}&&{\rm{\Delta }}{\hat{C}}_{b}^{\mathrm{AB}}=\sqrt{\displaystyle \frac{1}{{\nu }_{b}}\left({C}_{b,{\rm{A}}}+{N}_{b,{\rm{A}}})({C}_{b,{\rm{B}}}+{N}_{b,{\rm{B}}})+{C}_{b,\mathrm{AB}}^{\,2}\right)}.\end{eqnarray} \tag{ 21 }$$

It is important to note that these uncertainties are plotted for visualization purposes and do not directly enter the likelihood model. We perform an end-to-end validation of the pipeline by simulating all cross-spectra using 12 CMB realizations for each value of r ∈ {0.0, 0.1, 0.2} and a single PySM dust realization (Thorne et al. 2017) scaled to match the measured dust emission reported in BK15.

The fixed ν_b model is an approximation to the real data uncertainties, because signal and noise will generally have different ν_b values. This is due to two significant effects: the approximately isotropic Planck noise in our patch has systematically fewer degrees of freedom per band power than the anisotropic Polarbear noise, and the $E\to B$ leakage subtraction in Equation (11) suppresses the effective number of degrees of freedom at low ℓ. We combine the Polarbear- and Planck-derived ν_b by taking the geometric mean. In Section 6.5, we discuss this choice and define a goodness-of-fit metric for our likelihood model that is sensitive to systematic misestimation of the auto- and cross-spectrum uncertainties between the simulations and the real data.

We compute an upper limit on polarized B-mode Galactic synchrotron emission in our field at 150 GHz.

The power spectrum estimation formalism follows the formalism used to compute the cross-spectra with high-frequency Planck data. The autospectrum of the Planck Low Frequency Instrument (LFI) 30 GHz full missing map (Planck Collaboration et al. 2018c) as seen by Polarbear is computed and the noise bias is subtracted using the Planck FFP10 noise realizations. We validate our pipeline using a set of CMB simulations with r = 0.0 added to a PySM map of Galactic synchrotron emission as a fiducial signal model for the 30 GHz channel. Our simulations assume that dust foregrounds are negligible at 30 GHz and that synchrotron foregrounds are negligible at 150 GHz. We use the same fixed ν_b approximation as the high-frequency case using the signal cross-spectrum C_b,AB computed from the PySM models. Due to the Planck beam size at 30 GHz, we only consider the first five band powers corresponding to ℓ ≤ 300.

We model the synchrotron emission as a power law in ℓ and frequency. Specifically, we assume that the power spectrum takes the form

$$\begin{eqnarray}&&{D}_{{\ell },\mathrm{sync}}={A}_{\mathrm{sync}}{\left(\displaystyle \frac{{\ell }}{{{\ell }}_{0}}\right)}^{{\alpha }_{\mathrm{sync}}}{\left(\displaystyle \frac{\nu }{{\nu }_{0}}\right)}^{2{\beta }_{\mathrm{sync}}},\end{eqnarray} \tag{ 22 }$$

in brightness units where the 2 is due to the fact that D_ℓ is quadratic in signal amplitude. We assumed fixed values for the power-law indices α_sync = −1.18 taken from the highest Galactic latitudes in a recent measurement by S-PASS in conjunction with WMAP and Planck data (Krachmalnicoff et al. 2018). We assume a value of β_sync = −3.2 from the same analysis, which is consistent with the prior used in BK15. Following previous work, we choose pivots ℓ₀ = 80 and ν₀ = 23 GHz. We construct a one-parameter likelihood following Hamimeche & Lewis (2008), hereafter HL08, relating the amplitude A_sync to the autospectrum of the Planck 30 GHz map and assuming a fixed ΛCDM CMB lensing component with r = 0. We integrate over the Planck 30 GHz average bandpass from the LFI reduced instrument model available from the Planck Legacy Archive and the Polarbear design bandpass.

We find the likelihood of A_sync peaks at zero, with a 95% upper limit of A_{sync,23 GHz} ≤ 12.5 μK². Using the fixed prior on α_sync and β_sync, this amplitude corresponds to a 95% upper limit on the synchrotron contamination at our frequency band and ℓ = 80 of 2.7 × 10⁻⁴ μK². We do not include synchrotron contamination in our fiducial r likelihood model, as it is deeply subdominant to dust foregrounds at 150 GHz.

We find the cross-spectrum between Planck 30 GHz and Polarbear to be consistent with null. We do not directly use this spectrum in our estimate of synchrotron contamination, as this adds no meaningful constraining power on the synchrotron spectral index.

All 12 spectra (10 with Planck HFI and 2 with Planck LFI) from the real data are shown in Figure 13. The HFI spectra are used as the primary input to our foreground and cosmological parameter constraints described in Section 6.5.

Zoom In Zoom Out Reset image size

Polarized point sources contribute to the level of B-mode power in our maps as a Poisson noise term. At frequencies ≲150 GHz, the source number counts are expected to be dominated by blazars and flat spectrum quasars. ACTPol has reported 26 detections of radio sources during the first two seasons on a comparable fraction of the sky at the same frequency (Datta et al. 2019). Fourteen of these sources showed linear polarization at more than 3σ significance. The mean polarization fraction of these sources is 0.028 ± 0.005, in agreement with previous studies in Puglisi et al. (2018) and Bonavera et al. (2017).

We search our maps for statistically significant polarized point sources. We first apply a matched spatial filter to a high-resolution version of our maps, similar to Vieira et al. (2010) and Marriage et al. (2011). We detect 19 point sources in intensity above 5σ, but do not detect any sources in polarization at the same significance. This null result is consistent with forecasts at our frequency (Tucci et al. 2011). We therefore consider the emission from polarized sources below our detection flux.

We estimate the level of B-mode power in our maps from undetected radio sources using the PS4C forecasting tools presented in Puglisi et al. (2018). We convert our sensitivity into a 5σ equivalent detection flux of 60 mJy in polarization. The upper limit on the contamination at ℓ = 80 is estimated to be D_ℓ = 1.6 × 10⁻⁴ μK². This is approximately equivalent to Δr ≲ 0.002. We therefore neglect the point-source contamination in our parameter constraints.

We fit the B-mode power spectrum from our data and cross-correlation with Planck high-frequency data to a CMB and single dust component, as done by The BICEP2/Keck & Planck Collaborations et al. (2015), using a likelihood model similar to HL08 and Cardoso et al. (2008). For each bandpower, we write a n_freq × n_freq matrix of the measured cross-spectra ${\hat{{\boldsymbol{C}}}}_{b}$, where n_freq = 4 is the number of frequency channels considered. The diagonal elements of this matrix contain the noise bias. For the Planck channels, this is simply a consequence of computing the autospectrum of the full mission frequency maps. For Polarbear, this is obtained by adding the noise bias to the fiducial B-mode power spectrum described in Section 4.5. This makes the estimation of r more robust to misestimations of the Polarbear noise bias than simply using the autospectrum of the full coadd. We simulate artificially changing the Polarbear noise bias by 10% and find a small shift in the uncertainty but no significant shift in the reconstructed r value.

In the case that the underlying fields are isotropic, Gaussian, and measured on the full sky, ${\hat{{\boldsymbol{C}}}}_{b}$ follows a Wishart distribution with ${\sum }_{{\ell }}{{\boldsymbol{P}}}_{b{\ell }}(2{\ell }\ +1)$ degrees of freedom per band power, where ${{\boldsymbol{P}}}_{b}{\ell }$ is the binning operator defined in Section 4.5. We assume that our ${\hat{{\boldsymbol{C}}}}_{b}$ follows the same distribution, but with an effective number of degrees of freedom ν_b estimated using Equation (19) for our partial sky area. This is an approximation because the effective number of degrees of freedom for the Planck and the Polarbear maps differ somewhat, due to the anisotropy of the Polarbear noise. The choice between these sets of ν_b values has little influence on the results; in both cases, the simulations show that there is no bias on the estimation of r and that the width of the marginal posterior on r is compatible with the dispersion of the best-fit values. We choose to use the geometric mean of the two ν_b estimates for the analysis that we describe in the rest of this section.

Under these assumptions, the likelihood ${ \mathcal L }$ of a true spectrum ${{\boldsymbol{C}}}_{b}$ given measured ${\hat{{\boldsymbol{C}}}}_{b}$ is given by

$$\begin{eqnarray}&&-2\mathrm{ln}{ \mathcal L }=\sum _{b}{\nu }_{b}\left\{\mathrm{Tr}[{\hat{{\boldsymbol{C}}}}_{b}{{\boldsymbol{C}}}_{b}^{-1}]-\mathrm{ln}| {\hat{{\boldsymbol{C}}}}_{b}{{\boldsymbol{C}}}_{b}^{-1}| -{n}_{\mathrm{freq}}\right\},\end{eqnarray} \tag{ 23 }$$

up to a constant offset. We model the underlying ${{\boldsymbol{C}}}_{b}$ as a sum of CMB, dust, and noise components,

$$\begin{eqnarray}&&{{\boldsymbol{C}}}_{b}={{\boldsymbol{C}}}_{b}^{\mathrm{CMB}}+{{\boldsymbol{C}}}_{b}^{\mathrm{dust}}+{{\boldsymbol{N}}}_{b}.\end{eqnarray} \tag{ 24 }$$

Every element of the CMB cross-spectrum matrix ${{\boldsymbol{C}}}_{b}^{\mathrm{CMB}}$ is equal within a bandpower, because all spectra are computed in CMB temperature units. The CMB receives contributions from lensing and tensor B-mode signals,

$$\begin{eqnarray}&&{{\boldsymbol{C}}}_{b}^{\mathrm{CMB}}=r{{\boldsymbol{C}}}_{b}^{\mathrm{tens}}+{A}_{\mathrm{lens}}{{\boldsymbol{C}}}_{b}^{\mathrm{lens}}.\end{eqnarray} \tag{ 25 }$$

Here, r is the tensor-to-scalar ratio, A_lens is the normalized amplitude of the ΛCDM lensing signal, ${{\boldsymbol{C}}}_{b}^{\mathrm{tens}}$ is the binned tensor B-mode signal, and ${{\boldsymbol{C}}}_{b}^{\mathrm{lens}}$ is the binned lensing signal. The dust component is treated as a power law in ℓ and a modified blackbody in frequency, following Planck Collaboration Int. XXX. (2016), Planck Collaboration et al. (2018e), and subsequent work. We define a vector f(β_dust, T_dust) that represents the dust emission for each frequency bandpass converted into CMB temperature units, where β_dust and T_dust are the spectral index and temperature of the modified blackbody, respectively. The dust component of ${{\boldsymbol{C}}}_{b}$ for frequencies i, j can be written as

$$\begin{eqnarray}&&{{\boldsymbol{C}}}_{b,{ij}}^{\mathrm{dust}}={A}_{\mathrm{dust}}({f}_{i}{f}_{j}){\left(\displaystyle \frac{{\ell }}{{{\ell }}_{0}}\right)}^{{\alpha }_{\mathrm{dust}}},\end{eqnarray} \tag{ 26 }$$

where A_dust is the amplitude of the dust signal and α_dust is the power-law index in ℓ. We assume a pivot value of ℓ₀ = 80 and normalize f such that A_dust represents the dust emission at 353 GHz. It should be noted that dust foregrounds are not expected to be Gaussian nor to follow a power law in ℓ. The physics of interstellar dust may result in a polarized frequency scaling significantly different from a single modified blackbody, a more complex ℓ dependence, or spatial decorrelation between high frequencies and CMB channels. However, we do not have the sensitivity to meaningfully constrain more complex foreground models with this data set, and therefore only consider the fiducial case.

For each frequency channel, we integrate over the instrument bandpass. The Polarbear channel uses the design bandpasses from Arnold et al. (2012) and the Planck channels use the HFI-reduced instrument model.

The noise component ${{\boldsymbol{N}}}_{b}$ is entirely diagonal, as the noise between frequency bands and experiments is expected to be uncorrelated. The value of the diagonal elements is taken to be the power spectrum of the "signflip" (FFP10) noise realizations for the Polarbear (Planck) frequency channels using the cross-spectrum (autospectrum) pipelines.

The model contains six parameters: r, A_lens, α_dust, β_dust, T_dust, and A_dust. We fix the values of A_lens = 1 and T_dust = 19.6, and allow the other four to float. We apply Gaussian priors on α_dust = −0.58 ± 0.21 and β_dust = 1.59 ± 0.11, respectively, based on the marginal posterior from BK15 and the estimate of patch-to-patch variation found by Planck Collaboration Int. XXX. (2016). In our fiducial likelihood, we impose the prior r ≥ 0.

We define a goodness-of-fit criterion following HL08:

$$\begin{eqnarray}&&{\chi }_{\mathrm{eff}}^{2}=-2\mathrm{ln}{ \mathcal L },\end{eqnarray} \tag{ 27 }$$

and analogously find that, in the limit that ν_b ≫ n_freq and that the number of fit parameters is negligible compared to the total number of bins across all spectra, ${\chi }_{\mathrm{eff}}^{2}$ has expectation value and variance

$$\begin{eqnarray}&&\langle {\chi }_{\mathrm{eff}}^{2}\rangle \simeq {n}_{\mathrm{bins}}\displaystyle \frac{{n}_{\mathrm{freq}}({n}_{\mathrm{freq}}+1)}{2}\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&\mathrm{var}({\chi }_{\mathrm{eff}}^{2})\simeq 2\langle {\chi }_{\mathrm{eff}}^{2}\rangle ,\end{eqnarray} \tag{ 29 }$$

where n_bins is the number of ℓ bins in each spectrum.⁴⁴ This is consistent with a χ² distribution with a number of degrees of freedom equal to the total number of bins across all unique spectra in $\hat{{\boldsymbol{C}}}$, in our case 110. Figure 14 shows that the simulated ${\chi }_{\mathrm{eff}}^{2}$ values are consistent with expectations. The value of ${\chi }_{\mathrm{eff}}^{2}$ from the real data lies in the middle of the simulations with PTE = 70%.

Zoom In Zoom Out Reset image size

Our parameter constraints are shown in Figure 15. Our posteriors on α_dust and β_dust are dominated by the input prior, meaning the data have little additional constraining power on these parameters. The prior on α_dust is not critical for the estimation of r, since the dust power in the bins and spectra with the most constraining power on r is largely set by A_dust and β_dust. We include this prior because it results in a more efficient exploration of the parameter space by the Markov chain Monte Carlo (MCMC). Using a prior on β_dust is more important for our r constraint, since it is necessary for the MCMC chains to converge reliably. This prior improves our upper limit on r by ∼30% in simulations. We find evidence for dust B-modes in our patch with a best-fit value for A_dust of 3.6 μK². This is marginally lower than the value reported in BK15. The expected difference between the two results is nontrivial to compute, given the partial overlap in data sets and different observation strategies. We exclude zero dust foregrounds with 99% confidence. We find a 95% confidence level upper limit on the primordial tensor-to-scalar ratio of r < 0.90, after marginalizing over foreground parameters.

Zoom In Zoom Out Reset image size

We estimate the impact of our instrumental systematic and calibration errors on the final r constraint as follows. We add the upper bounds on systematic contamination reported in Table 6 to a reference theoretical power spectrum containing CMB (r = 0), dust, and noise. We analyze these spectra following the real data and find Δr_syst = 0.02. We additionally generate random multiplicative calibration errors according to the levels reported in Table 6 and run the component separation. We find that the bias on the best-fit value of r is smaller than Δr_cal ≤ 0.01 and the effect on the estimated uncertainty is less than 10%. We therefore neglect multiplicative calibration effects and only consider additive systematic errors. We do not attempt to quantify the impact of Planck systematics in our results.

Finally, for a sensitivity comparison to other experiments, we remove the prior on r and provide a best-fit estimate as the maximum of the marginal posterior. The statistical uncertainty is estimated as the most narrow interval that contains 68% of the integral of the distribution. We find $r={0.26}_{-0.30}^{+0.32}\ (\mathrm{stat}.)\,\pm 0.02\ (\mathrm{syst}.)$, where the statistical error refers to the narrowest interval containing 68% of the probability weight.