A Measurement of the Cosmic Microwave Background $B$-Mode Polarization Power Spectrum at Sub-Degree Scales from 2 years of POLARBEAR Data

TL;DR

This paper presents an improved measurement of the CMB B-mode polarization power spectrum at sub-degree scales using two years of POLARBEAR data, providing evidence for lensing-induced B-modes and testing consistency with the DM model.

Contribution

The study combines two seasons of POLARBEAR data to significantly reduce uncertainties in the B-mode power spectrum measurement and tests its consistency with DM predictions.

Findings

Rejects the null hypothesis of no B-mode polarization at 3.1C3

Measures a lensing amplitude parameter A_L = 0.60^{+0.26}_{-0.24} (stat+inst+foreground+multi)

Provides evidence for lensing-induced B-modes in the CMB

Abstract

We report an improved measurement of the cosmic microwave background (CMB) $B$-mode polarization power spectrum with the POLARBEAR experiment at 150 GHz. By adding new data collected during the second season of observations (2013-2014) to re-analyzed data from the first season (2012-2013), we have reduced twofold the band-power uncertainties. The band powers are reported over angular multipoles $500 \leq \ell \leq 2100$, where the dominant $B$-mode signal is expected to be due to the gravitational lensing of $E$-modes. We reject the null hypothesis of no $B$-mode polarization at a confidence of 3.1$\sigma$ including both statistical and systematic uncertainties. We test the consistency of the measured $B$-modes with the $\Lambda$ Cold Dark Matter ($\Lambda$CDM) framework by fitting for a single lensing amplitude parameter $A_L$ relative to the Planck best-fit model prediction. We obtain $A_L = 0.60 ^{+0.26} _{-0.24} ({\rm stat}) ^{+0.00} _{-0.04}({\rm inst}) \pm 0.14 ({\rm foreground}) \pm 0.04 ({\rm multi})$, where $A_{L}=1$ is the fiducial $\Lambda$CDM value, and the details of the reported uncertainties are explained later in the manuscript.

Tables

Table 1. Table 1 Uncertainties in polarization angle.

Angle uncertainty (∘)	Global	Wafer-Averaged
Beam uncertainties	0.11	0.33
Relative gain uncertainties	0.11	0.28
Non-ideality of HWP	0.18	0.19
Circular polarization of Tau A	0.10	0.05
HWP angle uncertainties	0.31	0.04
Pixel pointing uncertainties	0.01	0.05
Bolometer time constant	0.01	0.04
Filtering effect	0.02	0.04
Total uncertainty (∘)	0.40	0.48

Table 2. Table 2 PTEs resulting from the null test framework.

Patch		Average of ${\chi }_{\mathrm{null}}\left(b\right)$			Extreme of ${\chi }_{\mathrm{null}}\left(b\right)$			Extreme of ${\chi }_{\mathrm{null}}^2$ by $C_{\ell }^{EB/BB}$			Extreme of ${\chi }_{\mathrm{null}}^2$ by test			Total ${\chi }_{\mathrm{null}}^2$
		Pipeline			Pipeline			Pipeline			Pipeline			Pipeline
		A	B		A	B		A	B		A	B		A	B
RA4.5		84.0%	100%		70.8%	49%		86.0%	54%		64.0%	73%		35.4%	47%
RA12		70.8%	109%		29.8%	19%		66.2%	52%		57.8%	54%		51.0%	52%
RA23		95.0%	191%		40.8%	13%		69.4%	11%		42.8%	16%		05.0%	05%

Table 3. Table 3 Sources of foreground power and their expected power in D ℓ B ​ B superscript subscript 𝐷 ℓ 𝐵 𝐵 D_{\ell}^{BB} .

Foreground	Expected power in $D_{\ell }^{BB}$ (${10}^{-4}$ $\mu$K2)
	$\ell =500\text{–}900$	900–1300	1300–1700	1700–2100
Galactic dust	$63.5\pm 123.3$	$53.6\pm 102.4$	$48.8\pm 91.4$	$46.9\pm 85.1$
Galactic synchrotron	$1.4\pm 2.1$	$1.2\pm 1.9$	$1.1\pm 1.7$	$1.0\pm 1.5$
Radio & dusty galaxies	$13.4\pm 5.5$	$26.8\pm 9.8$	$24.8\pm 15.3$	$67.4\pm 21.9$
Total	$78.3\pm 123.4$	$81.6\pm 102.9$	$74.7\pm 92.7$	$115.3\pm 87.9$

Table 4. Table 4 Reported Polarbear  band powers D ℓ B ​ B superscript subscript 𝐷 ℓ 𝐵 𝐵 D_{\ell}^{BB} = = ℓ ​ ( ℓ + 1 ) ​ C ℓ B ​ B / 2 ​ π ℓ ℓ 1 superscript subscript 𝐶 ℓ 𝐵 𝐵 2 𝜋 \ell\left(\ell+1\right)C_{\ell}^{BB}/2\pi  from pipeline A

Central $\ell$	$D_{\ell }^{BB}$ [$\mu$K2]
0700	${0.053}_{-0.026}^{+0.029}$
1100	${0.074}_{-0.054}^{+0.057}$
1500	${0.045}_{-0.109}^{+0.116}$
1900	${0.352}_{-0.239}^{+0.260}$

Table 5. Table 5 Summary of the reported Polarbear  systematic uncertainties.

Type	Source of systematics	Effect on $D_{\ell }^{BB}$ [${10}^{-4}$ $\mu$K2]	Effect on $A_{BB}$
Instrument (Sec. 4.4)	Gain drift	$8.5$	$0.009$
	Differential gain	$9.3$	$0.010$
	Differential beam size	$0.4$	$0.000$
	Differential beam ellipticity	$0.1$	$0.000$
	Differential & Boresight pointing	$5.7$	$0.008$
	Instrument & Relative polarization angle	$6.7$	$0.008$
	Electrical crosstalk	$2.5$	$0.003$
	Total	$33.3$	$0.037$
Astrophysical foreground	Galactic dust (Sec. 5.1)	$63.5\pm 123.3$	$0.071\pm 0.138$
	Galactic synchrotron (Sec. 5.1)	$1.4\pm 2.1$	$0.002\pm 0.002$
	Radio & Dusty galaxies (Sec. 5.2)	$13.4\pm 5.5$	$0.019\pm 0.005$
	sub total	$78.3\pm 123.4$	$0.092\pm 0.138$
Analysis999The residual of the ground pickup removal and $E$-to-$B$ leakage subtraction and the sensitivity of the transfer function to different cosmologies are evaluated using the same methodology of PB14. The systematic bias and uncertainty from “Analysis” is subtracted from the measured band powers and propagated in pipeline A.	Ground pickup removal	$0.5\pm 1.7$	$0.001\pm 0.002$
	$E$-to-$B$ leakage due to filter subtraction	$2.5$	$0.003$
	sub total	$2.9\pm 1.7$	$0.003\pm 0.002$
	Total	$81.2\pm 123.4$	$0.095\pm 0.138$
Multiplicative effect	Absolute gain uncertainty (Sec. 6)		$\pm 3.0\%$
	Beam uncertainty (Sec. 6)		$\pm 1.0\%$
	Polarization efficiency (Sec. 3.4)		$\pm 3.3\%$
	Transfer functionaafootnotemark:		$\pm 3.9\%$
	Total		$\pm 6.0\%$

Equations

$$-2\ln \mathcal{L}(X_{\mathrm{meas}}|X_{\mathrm{cosm}})=$$

$$\displaystyle-\ln[X_{\rm cosm}+x_{0}]\Big{]}^{2}$$

$$-\ln [\sigma (X_{\mathrm{meas}}+x_0)\sqrt{2\pi }],$$

Full text

††thanks:

Corresponding authors: Yuji Chinone, Julien Peloton.

      chinoney@berkeley.edu, j.peloton@sussex.ac.uk

1. Introduction

The polarization of the cosmic microwave background (CMB) encodes broad cosmological information that is the focus of current and future generations of CMB experiments. The pattern of linear polarization separates into gradient-like $E$-mode and curl-like $B$-mode components. On one hand, $E$-mode polarization is generated by the same scalar density fluctuations in the ionized plasma before recombination that generate temperature anisotropies. In contrast, $B$-mode polarization is not generated by these scalar perturbations, and could be generated by either tensor perturbations (gravitational waves) from inflation or by conversion of $E$-modes to $B$-modes by gravitational lensing along the line of sight. On degree scales, where the inflation scenario predicts $B$-mode polarization of the CMB from primordial gravitational waves, no such signal has yet been detected (The Bicep2/Keck and Planck Collaborations: P. A. R. Ade et al., 2015; Planck Collaboration et al., 2016d; The Bicep2 and Keck Array Collaborations: P. A. R. Ade et al., 2016). In addition, gravitational lensing, which is not related to this primordial signal, induces a characteristic peak in the BB angular power spectrum at $\ell\sim 1000$. In this case, the primary $E$-mode spectrum is converted to $B$-modes by lensing from the large-scale distribution of matter. Several experiments have begun to measure this lensing $B$-mode signal, including Polarbear (The Polarbear Collaboration: P. A. R. Ade et al., 2014a), bicep2 (The Bicep2 Collaboration: P. A. R. Ade et al., 2014), Keck Array (The Bicep2 and Keck Array Collaborations: P. A. R. Ade et al., 2015), SPTpol (Keisler et al., 2015), cross-spectra between bicep2/Keck Array and Planck 2015 $B$-mode maps (The Bicep2/Keck and Planck Collaborations: P. A. R. Ade et al., 2015; The Bicep2 and Keck Array Collaborations: P. A. R. Ade et al., 2016), and actpol (Louis et al., 2016).

This paper reports results from the Polarbear telescope, which has a $3\farcm 5$ resolution, giving its best sensitivity through intermediate angular scales around the $\ell\sim 1000$ lensing peak. A previous Polarbear publication (The Polarbear Collaboration: P. A. R. Ade et al., 2014a), hereafter PB14, reported the first direct measurement of a non-zero $B$-mode signal, with modest significance. In this paper, we present an improved measurement of the angular power spectrum as measured by Polarbear in the same survey areas (a few degrees across for each sky patch), to greater depth. The procedure used to analyze the dataset shares many similarities with that described in PB14. In particular, the analysis is conducted blindly, with angular power spectra revealed only when the data pass a series of null tests and systematic error checks. We also incorporate many new features and improvements in the analysis as well as an improved discussion of the contribution of astrophysical foregrounds. We include new data taken between October, 2013 and April, 2014 (hereafter second season), but also a re-processing of the first season data (May, 2012 – June, 2013) in the light of the changes, which expands the total data volume with respect to PB14 by a factor of 61%. This means 62% of the total data is reprocessed and 38% of the total data is newly processed. We therefore refer the reader to descriptions in PB14 when necessary, emphasizing here only the additional steps.

In Section 2 we describe the two seasons of data collected by the Polarbear instrument. In Section 3 we discuss the calibration procedure, and the data analysis steps are outlined in Sections 4, 5 and 6. In Section 7 we show the power spectra results, and in Section 8 we draw our conclusions.

2. First and second season observations of the Polarbear instrument

Polarbear is a CMB experiment that has been observing from the 2.5 m Huan Tran Telescope since January 2012. The telescope is located at the James Ax Observatory at an elevation of 5,190 m in the Atacama Desert in Chile. The Polarbear receiver consists of an array of 1,274 transition edge sensor (TES) bolometers cooled to 0.3 K and observing the sky through lenslet-coupled double-slot dipole antennas. Within the array, the bolometers are grouped into 7 different wafers. More details on the receiver and telescope can be found in Arnold et al. (2012) and Kermish et al. (2012). A cold half-wave plate (HWP) is positioned on the sky side of cryogenic lenses in the receiver. While this HWP was stepped almost daily for the first half of the first season, it was not stepped during the rest of the first season or during the second season. We discuss the role of the HWP on the mitigation of instrumental systematic effects in Sec. 4.4.

The Polarbear observing strategy is described in PB14, and is summarized here. We observe three CMB fields, with each one visible for 6 to 8 hours per day. The three patches are centered at (RA, Dec)=(4h40m12s, $-45$∘00′), (11h53m0s, $-0$∘30′), and (23h1m48s, $-32$∘48′) which we call RA4.5, RA12 and RA23, respectively. We divide the observations into constant elevation scans (CESs) during which the telescope scans back and forth in azimuth at constant elevation. We call each sweep in azimuth a subscan. After approximately 15 minutes, when the patch has moved out of the field of view, we adjust the elevation and begin another CES.

The total observation time for the two seasons and the three CMB patches is 4,700 hours, corresponding to 33% of the total calendar time available for the two seasons. Of this time, 2,800 hours pass all the data quality checks described in Sec. 4.2 and are used to compute the power spectrum. As mentioned, this is an increase in data volume of 61% over PB14.

The process of extracting cosmological results from raw time-ordered data (TOD) can be summarized in three steps: calibration of the raw data, map-making, and power-spectrum estimation. In the following two sections, we describe these processes.

3. Calibration

Before creating maps of CMB polarization anisotropy, we characterize detector and telescope performance. This includes reconstructing the telescope pointing, measuring the beam, calibrating the detector time-ordered data (TOD), and determining the detectors’ polarization properties. We will describe each of these four steps in the following sections.

3.1. Pointing

To go from TOD to maps, we need to know where each detector was pointing as a function of time. We determine the pointing model using observations of bright extended and point-like millimeter sources selected from known source catalogs (Wrobel et al., 1998; Murphy et al., 2010) across a wide range in azimuth and elevation. The approach is similar to that of PB14, although the model had to be extended to handle the increased pointing data volume and sky coverage. Relative to PB14, we use three times more pointing data across a $\sim 19\%$ ($\sim 30\%$) larger range in azimuth (elevation). Five new parameters were added to the pointing model to enable the pointing to be reconstructed over a larger fraction of the sky. The first was a timing error, i.e., an offset in hour angle per period of time, caused by small synchronization errors between clocks and ephemerides. The next two, solar radiation flexure in azimuth and elevation, account for the small deformation of the telescope due to temperature gradients related to the Sun’s position. The final two, ambient temperature flexure in azimuth and elevation, are similar except with the deformation correlated with the ambient temperature (Matsuda, 2017). Neglecting the four new flexure terms would worsen the RMS pointing uncertainty by 8″.

The final pointing accuracy is similar to PB14, with an RMS scatter measured on known radio source positions of 27″ and 30″ for seasons 1 and 2 respectively.

3.2. Beam

We estimate the beam map and its effective window function, $B_{\ell}$, using dedicated observations of Jupiter, following the procedure outlined in PB14. However, for this analysis we conservatively discarded the observations of Jupiter taken when the atmospheric precipitable water vapor (PWV) exceeded a threshold of 4 mm (as was also done for CMB observations), despite no clear sign of detector saturation in the reported range of angular scales. As a consistency check we estimated the $B_{\ell}$ for each observation season separately to test potential variation over time and found no significant deviation. Consistent results have been derived using Saturn observations. We tested that the deconvolution of the bolometer time constants induced negligible change on the beam properties and window function. The main lobes of the beam are well approximated by a Gaussian core having $3\farcm 5\pm 0\farcm 1$ FWHM plus a diffraction tail asymptotically decaying as $1/\theta^{3}$ where $\theta$ is the radial coordinate of the beam profile (Hasselfield et al., 2013). The $1/\theta^{3}$ decay is a natural consequence of the presence of a Lyot stop in the optical chain as well as the finite detector size and difference in the spectral response of the detectors. In Figure 1 we show the mean radial profile of the beam obtained as an inverse noise weighted average of the beam profiles estimated for each single Jupiter observation in our dataset together with the best fit results of the beam profile model. The difference between the $B_{\ell}$ bandpowers derived from the best fit profile model (orange curve) and the reference one derived from the Jupiter maps (red curve) is less than $1.5\%$ in all the $\ell$ range considered in this work.

The median beam ellipticity measured across the array is 5% and the median difference in ellipticity measured for two channels in each focal plane pixel is 1.6%, with a subdegree difference in the ellipticity orientation. In the systematics simulations of Sec. 4.4, we use the full distributions of those two quantities to quantify the potential biases due to beam effects on the $B$-mode power spectrum.

The Jupiter-measured beam is then symmetrized radially and convolved with a pointing error term for each of the three fields to calculate the per-field beam as a function of angular multipole that will be used in the power spectrum analysis. The profile of the symmetrized beam is consistent with the radial profile data as shown in Figure 1. We refer the reader to PB14 for more details on the process.

The pointing error calculation has changed slightly from PB14, although the underlying model is unchanged. We assume the real instrumental pointing is Gaussian distributed around the reconstructed pointing with a standard deviation referred to as the pointing jitter. In PB14, the pointing jitter and uncertainty on the jitter were calculated by fitting a beam width to subsets of the data to obtain a probability density function for the real jitter. However, further analysis showed that the underlying CMB anisotropies at each source’s location are a significant bias in the fitted beam width, an effect not captured in the previous analysis. Like PB14, each point source is fitted using free parameters for source position, amplitude and beam width. The uncertainty on the parameters is then estimated from the scatter of the fitted parameters for point sources simulated using the central fitted values against a $\Lambda$CDM CMB temperature background. The pointing jitter and its uncertainty is then backed out of the fitted parameters and their uncertainties obtained from simulations. Using this algorithm, we find pointing jitter values of $24\farcs 5\pm 2\farcs 7$, $24\farcs 5\pm 5\farcs 9$, and $57\farcs 1\pm 6\farcs 1$ for RA4.5, RA12, and RA23, respectively. These jitter values are consistent with those found in the PB14 analysis, but with significantly reduced error bars due to the correct treatment to the underlying CMB fluctuations.

3.3. Detector gains

We calibrate the TOD to physical Rayleigh-Jeans Kelvin temperature units (K${}_{\text{RJ}}$) in a multi-step process, following the methods of PB14. First we determine the relative calibration between detectors over time using a combination of an internal thermal source and Saturn observations. With this relative calibration in hand, we can turn the TOD into a temperature map. Finally, we determine the absolute calibration by looking at the temperature anisotropy CMB power spectrum as discussed in Sec. 6. In this section, we will focus on the relative calibration process.

In brief, the relative calibration of the TOD proceeds in three steps. The first step combines an internal thermal calibration source (the stimulator) with the flux from Saturn. We discuss the complications introduced by Saturn’s rings in Appendix A. Next, we correct for the effect due to the polarized emission from the stimulator, which is then rotated by the HWP. We thus define the HWP angle-dependent template by considering observations of a variety of astrophysical sources at different HWP rotation angles. In addition to the HWP position, the polarized response also depends on the status of the stimulator (a few hardware changes were made over the observing period) and potentially a long-term temporal drift. As a consequence we generate a new polarization template for each of six epochs (four for season 1 and two for season 2) corresponding to changes in the stimulator hardware. Finally, we combine the stimulator and Saturn observations to calculate the relative gain, which is the conversion from electrical current into K${}_{\text{RJ}}$ units, for each detector and observation.

The median variation of the relative gain calibration for all detectors between two consecutive stimulator measurements is 0.5%; we use this information for the systematic error estimate to quantify the impact of the drift of the gains on the $B$-mode power spectrum (see Sec. 4.4).

3.4. Polarization angle

As in PB14, we determine the relative detector polarization angles using daily observations of Tau A (the Crab nebula). We then further improve the accuracy in the global polarization angle by nulling the $C_{\ell}^{EB}$ cross-spectrum (see Sec. 6), assumed to be zero. Here we summarize this approach, with more details available in PB14. We also tabulate the achieved polarization angle uncertainties from this analysis, which are important to the systematics budget for the $C_{\ell}^{BB}$ power spectrum (see Sec. 4.4).

We estimate polarization angles and efficiencies by observing Tau A daily for roughly 30 minutes and fitting each detector’s TOD from the two seasons of data to a reference Tau A map from the IRAM111Institut de Radioastronomie Millimétrique 30 m telescope observation (Aumont et al., 2010) using the Polarbear Jupiter-based beam and known pointing information. The individual-pixel angular uncertainty in each wafer is estimated to be 1$\fdg$2 using the measured angular dispersions in each wafer. In addition, the total systematic uncertainties for the array-averaged and wafer-averaged angles are 0$\fdg$40 and 0$\fdg$48. A detailed breakdown can be found in Table 1. By combining the two seasons of data, the systematic uncertainties from the beam and relative gain uncertainties causing temperature-to-polarization leakage are improved compared to PB14. However, the polarization angle difference between the seasons described below shows that the HWP angle uncertainty is now more significant than what was found in PB14. These changes make the total polarization angle uncertainty at the same level to that of PB14.

The array-averaged polarization efficiency was 97.2%, consistent with the expected value of 97.6% from HWP and anti-reflection coatings. The systematic uncertainty for the polarization efficiency is 1.6%, which introduces a 3.3% multiplicative uncertainty in the CMB $B$-mode power spectrum.

We also checked the stability of the array-averaged polarization angles over time. Figure 2 shows a histogram of the measured array-averaged Tau A polarization angles. Each angle is calculated by computing $(1/2)\arctan\left(\sum{U_{j}}/\sum{Q_{j}}\right)$ across map pixels within 10′ of the center of Tau A, where $Q_{j}$ and $U_{j}$ are Stokes $Q$ and $U$ in the map pixel $j$. Two effects are apparent: (1) the average angle shifted by 0$\fdg$9 between the two seasons; (2) the scatter was larger in season 1. We believe that both effects are due to the uncertainty of the HWP position. The polarization angle calculation is based on the commanded HWP positions (quantized at 11$\fdg$25). Any offset, $\Theta_{\rm err}$, between the commanded and actual HWP position will shift the calculated polarization angle by 2$\Theta_{\rm err}$. The first season data was taken at a number of HWP positions, (where multiple HWP offsets lead to the increased scatter in the polarization angle), while the second season data was taken at single HWP position (where the single HWP offset shows up as an polarization angle offset between first and second seasons). We can reconstruct the distribution of HWP offsets that would lead to the observed polarization angle variations, and find the typical magnitude to be 0$\fdg$28 (a polarization angle uncertainty of 0$\fdg$56). After unblinding the $C_{\ell}^{EB}$ spectrum, we found that season-by-season $C_{\ell}^{EB}$ nulling led to a consistent shift in the polarization angle. The two-season combined Tau A polarization angle was $150\fdg 4\pm 0\fdg 2({\rm stat})\pm 0\fdg 8({\rm sys})$ in equatorial coordinates using the $C_{\ell}^{EB}$-derived polarization angle described in Sec. 6.

4. Data analysis

As mentioned, the two-season dataset contains about 61% more data than the PB14 dataset. This increase corresponds to new data from the second season of observation, as well as a re-analysis of the data from the first season with the improvements made in calibration and data selection, as detailed in other sections. The main update to the data analysis with respect to PB14 has been the implementation of a second and complementary pipeline to analyze the dataset, for more robust results through consistency checks and improved systematic error control. While the calibration of the TOD (see Sec. 3) is shared between the two pipelines, they perform map-making and the power-spectrum estimation differently.

4.1. Map-making and power-spectra estimation

We adopt two independent and algorithmically different pipelines in the analysis: the first pipeline (hereafter “pipeline A”) is based on the MASTER method (Hivon et al., 2002) and was described in PB14, and the second pipeline (hereafter “pipeline B”) is based on the work described in Poletti et al. (2016). The main steps performed by the two pipelines to estimate sky maps, angular cross-spectra, and null cross-spectra are described in Fig. 3.

The pipeline-A and pipeline-B map-making algorithms differ slightly in their filtering of the calibrated data and substantially in the way they project the filtered TOD into maps. Both perform high-pass and azimuthal filters to remove atmospheric noise and ground pickup, respectively. Pipeline A performs these two operations sequentially while pipeline B performs them simultaneously.222This is desirable in general as the templates that need to be filtered out are not always orthogonal from the outset. However, this step makes the map-making process more complex than filtering the non-orthogonal templates one after another. In addition, pipeline A applies a low-pass filter prior to map-making. The pipeline-A map-making algorithm projects the time-domain data into maps without accounting for the filtering-induced power suppression, but accounts for this during power-spectrum estimation. Combined with the sequential filtering, this allows a substantial speedup compared to the other pipeline. Pipeline B accounts for the filtering of the TOD and produces unbiased maps. The two pipelines also use different pixelization schemes. Pipeline A projects the time-domain data onto flat sky maps using a cylindrical equal-area projection, with map pixels of width $2\arcmin$. Pipeline B projects the time-domain data onto curved-sky maps using the HEALPix333http://healpix.jpl.nasa.gov/ pixelization ($N_{\rm side}=2048$, i.e., map pixel width $\sim 1\farcm 7$; Górski et al. (2005)). Given the relatively small fraction of sky observed (a few degrees across for each patch), and the resolution of the telescope ($3\farcm 5$), both approaches remain equivalent. In both pipelines, the data are combined into chunks which are later cross-correlated to avoid noise biases in the estimated power spectra. After the data selection (see Sec. 4.2), pipeline A splits the two-season datasets into 214, 212, and 269 daily maps for RA4.5, RA12, and RA23, respectively; pipeline B splits the datasets into eight chunks of data for each patch corresponding to roughly 1.5 cumulative months of data each. The full-season temperature and Stokes $Q$ and $U$ maps of RA23 produced by the two pipelines are plotted in Figure 4. The resulting polarization white-noise levels reach 7, 6, and 5 $\mu$K-arcmin (10, 7, and 6 $\mu$K-arcmin with the beam and filter transfer function divided out) for RA4.5, RA12, and RA23, respectively, for pipeline A.

In both pipelines the power-spectrum estimators are based on the pure pseudo-$C_{\ell}$ technique (Smith, 2006) although they differ in a few aspects, such as the computation of the mode-mode mixing matrix, and the handling of the off-diagonal elements controlling the level of $E$-to-$B$ leakage. For pipeline A, the formalism was described in PB14. For pipeline B, we use the software X2pure (Grain et al., 2009, 2012; Ferté et al., 2013). As described in Poletti et al. (2016), we perform a map-domain removal of the noisiest modes in order to control large-scale noise. Figure 4 shows the effective maps after performing this mode removal. We estimate the statistical uncertainty on the $C_{\ell}^{EB}$ and $C_{\ell}^{BB}$ spectra from 500 signal and noise Monte Carlo simulations (100 for pipeline B) as described in Sec. 7. For pipeline A, thanks to a realistic uncertainty estimate as well as overall larger dataset and improved calibration, this release achieves almost twice the sensitivity of PB14. As a result of applying different treatments to the data, the estimated statistical uncertainties from the two pipelines are different (see Poletti et al. (2016) for a complete comparison of methods), being on average 20% higher for pipeline B. Part of this difference is also attributed to the difference of resulting sky area used in the two analyses.

In PB14, Polarbear implemented a blind analysis method in our measurement of the $B$-mode polarization, and we have adopted the same procedure for the two-season analysis. In the following, we describe the data selection, analysis validation by a null-test framework, instrumental systematic-error estimation, and foreground estimation including biases from contamination due to point-like extra-Galactic and polarized Galactic foregrounds.

4.2. Data selection

In this section, we describe the data quality cuts. As mentioned in Sec. 2, 60% of the total CMB observation data (2,800 out of 4,700 hours) ends up being included in the science analysis of this work. We do not select the other 40% due to some combination of bad weather, incomplete observations, or hardware glitches. Note that the thresholds for the latter were validated and finalized while running the suite of null tests described in the next section, and before unblinding.

The data quality checks proceed in three stages. First, we require a successful measurement of the detector gains for each CES (see Sec. 3.3 for details on how detector gains are estimated). Second, we discard data based on conditions such as the level of PWV in the atmosphere, the angle between the observed patch and the Sun or the Moon, as well as measures of the quality of the data, such as the bolometer yield being too low, problematic scan length, bad azimuth encoder data, or the array median gain changing too rapidly, for example. This step can remove data from an entire CES or some of the detector channels within a CES, leaving $\sim$ 10,000 and $\sim$4,000 CES for season 1 and season 2 respectively. Note that, thanks to the improved calibration and slightly different data selection, the definition of season 1 here is not strictly equivalent to that used in PB14; the change increased by approximately 15% of the observation-hours that went into the maps. In the end, 2,800 observation-hours pass both the first and second data quality checks. The volume of data from the second season is smaller than the one from first season mainly because the period of observing time is smaller (13 months for first season and 6 months for second season). The third step goes further and consists of defining subscans (a constant-velocity segment of a timestream) and among them, identifying contaminated subscans (such as by finding glitches in them). As in PB14, we discard all data obtained while the telescope is accelerating.

4.3. Analysis validation: null tests

We perform a suite of null tests to evaluate the calibration, data-selection criteria and filtering methods, and to test for unknown systematic errors, before unblinding the data. This task consists of iteratively running the null-test framework described in PB14, with the addition of three new data splits, until a set of predefined criteria (described below) are passed. Due to the non-negligible extra computational cost of pipeline B with respect to pipeline A,444A single map-making run of pipeline B uses roughly as much computation time as one full run of the null-test framework of pipeline A. the full set of null tests (see Sec. 4.3.1) is performed only by pipeline A as detailed in PB14. We primarily used these null tests by pipeline A to define the data selection criteria. Once the science dataset is defined, pipeline B performs a subset of the data splits probed by pipeline A. The way pipeline B computes null-spectra differs slightly from pipeline A. It first splits the dataset in two according to the null-test criterion; second, it computes unbiased maps of 4 (or 8) disjoint subsets of each split; null-maps are then computed by taking the difference of maps that belong to different splits; and finally it computes and co-adds the cross-spectra of the null maps to form a null spectrum. Nevertheless, pipeline B performs further checks allowed by the production of null maps in addition to null spectra. These include the visual inspection of null maps, the study of the individual cross spectra of null maps, and their distribution (the same data can contribute to several cross-spectra). More specifically, if the band powers of a null spectrum have a probability-to-exceed (PTE, see also Sec. 4.3.2) less than 5%, we study the cross-spectra of which it is the average and we make sure that the low PTE is not determined by the cross-spectra involving the same particular portion of the data set. None of those extra checks found signs of contamination or inconsistencies in the dataset. We emphasize that this large variety of the validation tests is made possible by combining the complementary strengths of the two pipelines.

4.3.1 Data splits

The null tests are performed for several interesting splits of the data, chosen to be sensitive to various sources of systematic contamination or miscalibration. In addition to the null tests we performed in PB14 (“first half versus second half of the dataset”, “rising versus setting”, “high elevation versus low elevation”, “High gain versus low gain”, “Good versus bad weather”, “pixel type”,555Each detector wafer has two different pixel polarization angles. “left versus right side of the focal plane”, “left- versus right-going subscan” and “Moon distance”), we introduced the following three new tests to check the difference between seasons and our possible concerns:

• •

“First season versus second season”: probing seasonal variation on year-long time-scales. This test is sensitive to systematic changes in the calibration, beams, telescope, and detectors.

• •

“Sun distance”: checking for residual contamination after setting the Sun-proximity threshold for an observation to be considered for analysis.

• •

“Sun above the horizon versus sun below the horizon”: checking for contamination from the far sidelobe of the beam, and systematic changes of the pointing due to the small deformation of the telescope by solar heating.

The 12 null tests are used to analyze the dataset, and the correlations between tests are taken into account in the analysis by also running the same suite of null tests on Monte Carlo (MC) simulations as described below. For the 3 sky patches, pipeline A runs then a total of 36 null tests, while pipeline B focuses on a subset of 11 null tests.

4.3.2 Analysis

In the null spectra, for each band power bin $b$, we calculate the statistic $\chi_{\rm null}(b)\equiv\hat{C}_{b}^{\rm null}/\sigma_{b}$, where $\sigma_{b}$ is an MC-based estimate of the corresponding standard deviation, and its square $\chi_{\rm null}^{2}(b)$. $\chi_{\rm null}(b)$ is sensitive to systematic biases in the null spectra, while $\chi_{\rm null}^{2}(b)$ is more sensitive to outlier bins. To probe for systematic contamination affecting a particular power spectrum or null-test data split, we calculate the sum of $\chi_{\rm null}^{2}(b)$ over $500<b<2100$ by spectrum (“$\chi^{2}_{\rm null}$ by spectrum”), and the sum for both spectra for a specific test (“$\chi^{2}_{\rm null}$ by test”). We consider both $C_{\ell}^{EB}$ and $C_{\ell}^{BB}$ in the null tests in order to investigate sources of spurious $B$-mode signals. We require that each of these sets of PTEs are consistent with a uniform distribution, as evaluated using a Kolmogorov-Smirnov test, requiring a $p$-value (probability of seeing deviation from uniformity greater than that which is observed given the hypothesis of uniformity) to be equal to or greater than 5%. These distributions are consistent with the uniform distribution and Figure 5 shows the PTE distribution of the $\chi^{2}_{\rm null}$ for the two pipelines.

We create test statistics based on these quantities to search for different manifestations of systematic contamination. The five test statistics are PTEs from (1) the average value of $\chi_{\rm null}$; the extreme value of $\chi^{2}_{\rm null}$ (2) by bin, (3) by spectrum, and (4) by test; (5) the total $\chi^{2}_{\rm null}$ summed over the twelve null tests. In each case, the result from the data is compared to the result from simulations, and PTEs are calculated. Finally, we combine each of the test statistics, and calculate the final PTE, requiring it to be equal to or greater than 5%. Table 2 shows the PTE values for the final configuration.

Comparing the most significant outlier from the five test statistics with that from signal and noise simulations, pipeline A (pipeline B) gets PTEs of 71.8% (77%), 65.2% (16%), and 16.6% (13%) for RA4.5, RA12, and RA23 respectively. We therefore achieve the requirements described above, finding no evidence for systematics or miscalibration in the Polarbear dataset used for the analysis and in the analysis process itself.

4.4. Analysis validation: simulations of instrumental effects

We describe in this section the signal-only simulations used to determine the effect of uncertainties in the instrument model on the power spectrum listed in Table 5 (sub-section “Instrument”). We investigate nine systematic instrumental effects: uncertainty in instrument polarization angle; uncertainty in relative pixel polarization angles; uncertainty in instrument boresight pointing model; differential pointing between the two detectors in a pixel; the drift of the gains between two consecutive thermal source calibrator measurements; relative gain calibration uncertainty between the two detectors in a pixel; crosstalk in the multiplexed readout; differential beam size; and differential beam ellipticity.

The pipeline used to analyze those systematic effects was described in PB14, while here we describe two improvements: (1) the systematics pipeline follows exactly all the data analysis pipeline steps, and (2) all instrumental systematics studied are now included within the systematics pipeline. The first modification makes the systematic error study more comprehensive than in PB14, i.e., the simulations performed here also include time-domain filtering and the cross-correlation of submaps. The second modification allows us to have a common framework for all effects, and therefore the same metric to analyze the results (i.e., we now include beam effects and electrical crosstalk, which were analyzed separately in PB14). While the procedure was developed for both data analysis pipelines, we report here only the results from pipeline A due to computing constraints.

A major change in the hardware configuration with respect to PB14 is the position of the HWP. As mentioned in Sec. 3.4, we recorded the data beginning from the middle of the first season without stepping the HWP, in order to decrease the scatter in the polarization angle of the detector. We found that this results in less mitigation of some instrumental effects such as those related to instrument and pixel polarization angle, drift of the gains, or crosstalk in the multiplexed readout. We note, however, that pointing-related effects remain stable with respect to PB14, and we reduce the uncertainty coming from beam-related effects mainly thanks to an improved characterization.

All nine effects and their combination were found to produce spurious BB power well below the statistical uncertainty in the measurement of $C_{\ell}^{BB}$, as also shown in Fig. 6. The individual effects are combined linearly to give the total contamination.

5. Foregrounds

Polarized Galactic and extra-Galactic foregrounds are a potential contaminant to the CMB, and a particular concern for the very faint $B$-modes. As described in PB14, there are four foregrounds sources of interest: the first two are polarized Galactic dust and synchrotron emission, dominating at large angular scales (down to a few arcminutes) at intermediate and high Galactic latitudes; the second two are emission from polarized radio and dusty galaxies, on scales of a few arcminutes scale and smaller. In this section, we describe how we estimate the band power contribution from each of these sources. Note that the methods employed to estimate the polarized Galactic dust and synchrotron contamination is different from the one described in PB14, while the method used to estimate contamination from radio and dusty galaxies remains similar. Table 3 reports the mean values and 68.3% confidence intervals obtained for the band powers from the combination of the three patches for both components, and Figure 7 shows the estimated foreground contributions to $\ell\left(\ell+1\right)C_{\ell}^{BB}/2\pi$ at a 68.3% confidence level. As can be seen, the total foreground contribution is estimated to be small in any band power, though not completely negligible.

Note that the methods employed to estimate the polarized Galactic dust and synchrotron contamination is different from the one described in PB14, while the method used to estimate contamination from radio and dusty galaxies remains similar. Indeed, in PB14 the estimates of diffuse polarized foregrounds in the observed fields was based on models built upon the public data from wmap and Planck at the time for the dust (Bennett et al., 2013; Planck Collaboration et al., 2014), as well as the QUIET levels reported for synchrotron (QUIET Collaboration et al., 2012). Between the PB14 and this paper, the Planck 2015 polarization sky maps, including in particular the 353 GHz as the main tracer of the polarized Galactic dust emission, have become available (Planck Collaboration et al., 2016a).666https://www.cosmos.esa.int/web/planck/pla We therefore use this new data as the main reference for both the dust and synchrotron estimates reported below. It is relevant to notice also that, while carrying out the present analysis, we discovered a conversion error in the evaluation of the contamination from synchrotron emission in PB14, which however does not change significantly the conclusions there, and that has been reported in a separate erratum.

5.1. Polarized Galactic dust & synchrotron

The polarized Galactic foregrounds, synchrotron and thermal dust, dominate at frequencies smaller and larger than 70 GHz, respectively (Planck Collaboration et al., 2016c); at intermediate and high Galactic latitudes, and at the degree angular scale, the frequency of foreground minimum seems to vary substantially, almost equally distributed between 60 and 90 GHz in the regions where both foregrounds are detected with high significance (Krachmalnicoff et al., 2016). These Galactic foregrounds substantially contaminate CMB $B$-mode measurements at all frequencies on large angular scales, even in the cleanest regions of the sky, as shown in several recent studies (Planck Collaboration Int. XXX, 2016; Krachmalnicoff et al., 2016). Nevertheless, Galactic foregrounds are expected to be sub-dominant with respect to the lensing $B$-modes at arcminute scales, and our purpose here is to assess their relevance and provide upper limits on their contribution to our observations.

For thermal dust, we adopt the following procedure. We take the publicly available Planck 2015 sky map at 353 GHz as a tracer of polarized emission from thermal dust. To avoid noise bias in the computation of spectra, we calculate cross-spectra using half-mission (HM) splittings using the X2pure power spectrum estimator. We evaluate statistical errors by means of white noise MC simulations, using the pixel-pixel noise covariance matrices of the input maps. The low signal-to-noise ratio of Planck 2015 polarization maps at high Galactic latitudes prevents the estimation of $B$-mode spectra directly in the regions corresponding to our patches, because of their reduced size and the noise level at the small angular scales we consider. Therefore, in order to provide an upper limit on the amplitude of thermal dust, we compute power spectra on larger circular regions, with $10^{\circ}$ radius, centered on our sky patches. We extrapolate the measured amplitude of the spectra at $\ell\simeq 80$ (a multipole bin between 60 and 99) in these regions to higher multipoles by applying the power law scaling ${\cal D}_{l}\propto\ell^{\alpha_{d}}$ with $\alpha_{d}=-0.42\pm 0.02$ (Planck Collaboration Int. XXX, 2016). We then scale the Planck 2015 measurements to the Polarbear frequency assuming a modified blackbody spectral dependence for the thermal dust, with temperature $T_{d}\simeq 19.6$ K and $\beta_{d}\simeq 1.59\pm 0.14$ (Planck Collaboration Int. XXII, 2015). For computing the extrapolation in frequency we take into account both the Planck 2015 and the Polarbear frequency band passes, and for the extrapolation in $\ell$ we consider the actual Polarbear band power window function. To account for the fact that we are measuring the foreground amplitude on larger regions, we include in the error budget the Gaussian approximation of the signal sample variance evaluated for the actual sky area of the Polarbear patches.

A similar procedure has been used to estimate the amplitude of polarized synchrotron emission. We computed power spectra for synchrotron $B$-modes by cross-correlating the Planck 2015-Low Frequency Instrument 30 GHz map at the effective frequency of 28.4 GHz (Planck Collaboration et al., 2016a) with the wmap-K map at a frequency of 22.8 GHz (Bennett et al., 2013). Spectra are computed on the same circular regions of $10^{\circ}$ radius. The measured amplitude at $\ell\simeq 80$ is rescaled in frequency considering a power-law frequency dependence with $\beta_{s}=-3.12\pm 0.02$ (Fuskeland et al., 2014) and in multipole considering ${\cal D}_{l}\propto\ell^{\alpha_{s}}$ with $\alpha_{s}=-0.31\pm 0.13$ (Planck Collaboration et al., 2016c).

We note that for both polarized Galactic dust and synchrotron, the limiting factor in our calculation of the foreground mean values and uncertainties comes from the Planck 2015 polarization noise level at such small angular scales. The level of contamination reported in this paper is therefore higher than our previous PB14 estimates, which was based on available models at the time. In particular, the dominant contribution from Galactic dust, represented by the top of the orange area in Figure 7, raises by a factor of about 5 with respect to PB14.

5.2. Radio & dusty galaxies

We estimate the radio and dusty galaxy power by drawing 10,000 realizations from the distributions described below. In both cases, we take recent measurements of the temperature power and then convert to polarization power using an estimate of the mean square polarization fraction. For radio galaxies, the temperature power estimate comes from the De Zotti et al. (2005) model at 150 GHz, which is scaled to the Polarbear effective frequency according to the measured spectral index for the radio sources of $\alpha_{\rm rg}=-0.90\pm 0.20$ from George et al. (2015) (hereafter G14). We ignore shot noise due to finite sky area, but assume a 10% calibration uncertainty on our nominal masking threshold of 25 mJy (i.e., the modeled threshold is $25\pm 2.5$ mJy). We take the recent estimate from SPTpol of the mean square polarization fraction for synchrotron sources of $\left<p^{2}\right>=(1.42\pm 0.15)\times 10^{-3}$ (T. Crawford, private communication).777https://cmb-s4.org/CMB-S4workshops/index.php/File:Sptpol_ptsrc_polfrac_500d.pdf These numbers are consistent with a recent analysis of brighter synchrotron sources in the 143 GHz Planck 2015 maps which found $\sqrt{\left<p^{2}\right>}=0.043\pm 0.0018$, corresponding to $\left<p^{2}\right>=(1.85\pm 0.15)\times 10^{-3}$ (Bonavera et al., 2017).

We use direct observations of the dusty galaxy power in temperature from G14. These include a Poisson term with $D_{3000}^{P}=9.16\pm 0.36$ $\mu$K2 at the spt frequency (150 GHz) and a spectral index of $\alpha_{P}=3.267\pm 0.077$, and a clustered term of the form $\ell^{0.8}$ normalized to $D_{3000}^{C}=3.46\pm 0.54$ $\mu$K2 with a spectral index of $\alpha_{C}=4.27\pm 0.2$. Data on the polarization fraction of these dusty galaxies are still poor, but they are expected to have lower polarization fractions than synchrotron sources. We conservatively draw the mean square polarization fraction $\left<p^{2}\right>$ from a uniform distribution between $1\times 10^{-4}$ and $1.57\times 10^{-3}$ (i.e., 1% polarized to the +1 sigma limit from the G14 synchrotron measurement).

For each realization, we multiply the inferred dusty and radio galaxy spectra by the appropriate window functions to directly compare to the measured band powers.

6. Calibration using CMB spectra

The absolute gain calibration is performed differently by the two pipelines. Pipeline A combines the estimate of $C_{\ell}^{TT}$ from each patch into one single estimate according to their statistical and beam uncertainties, and then estimates the absolute gain by fitting the patch-combined $C_{\ell}^{TT}$ to a theoretical Planck 2015 $\Lambda$CDM spectrum.888All references to the Planck 2015 $\Lambda$CDM model in this work refer to the best-fit values for the base$\_$plikHM$\_$TT$\_$lowTEB$\_$lensing configuration. Finally the correction is applied to the individual maps. Pipeline B estimates the absolute gain for each patch by first cross-correlating the Polarbear temperature maps with the foreground-cleaned Planck 2015 temperature map produced by the Spectral Matching Independent Component Analysis (SMICA) method  (Planck Collaboration et al., 2016b). These cross-spectra are then compared to the corresponding Polarbear auto-spectra, before combining the estimates into a single absolute gain factor.

Both methods give consistent results, however the second method gives a higher uncertainty due to the propagation of the noise in Planck 2015 data: 3.0% uncertainty on the absolute gain factor for pipeline A, and 4.3% uncertainty for pipeline B. We note, however, that the method used by pipeline B has the advantage of being more cosmology independent. These calibration factors are then applied to all spectra by the respective pipelines to produce gain-calibrated spectra. The beam uncertainty is estimated from uncertainty in the point-source-derived beam-smoothing correction (Sec. 3.2), and the variation in that correction across each field. We shift the simulated beam by $1\sigma$, and find this beam shift leads to a $\pm 1.0\%$ change in the best-fit lensing $B$-mode amplitude. We include this $1.0\%$ uncertainty as part of the multiplicative error budget. A complete breakdown of the multiplicative uncertainties can be found in Table 5.

The global instrument polarization angle correction $\Delta\psi$ is obtained by fitting the patch-combined $C_{\ell}^{EB}$ to $2\Delta\psi C_{\ell}^{EE}$ (Keating et al., 2013) for both pipelines. The best-fit value and statistical uncertainty in the global $C_{\ell}^{EB}$-derived instrument polarization angle correction is $-0\fdg 79\pm 0\fdg 16$ ($-0\fdg 67\pm 0\fdg 17$) for pipeline A (pipeline B). Combining this value with the relative shift from the Tau A-derived angle (see Sec. 3.4), this result is consistent with the results obtained previously in PB14. Finally, pipeline A applies the polarization angle correction to the individual maps and pipeline B to the power spectra.

7. $B$-mode power spectrum results

The final analysis procedure follows closely the final first season analysis described in PB14. For each patch, and each spectrum, we first form the covariance matrices using our set of 500 MC simulations (100 MC for pipeline B), and we then add the beam and gain calibration uncertainties as described later in this section. We note that the use of MC simulations to construct the band power covariance matrices gives more realistic estimates than the analytical expressions used in PB14. Finally for both pipelines, a single estimate of the power spectra from the three patches is created using the measured band powers and their covariance matrices.

The band power covariance matrix accounts for sample variance and instrumental noise variance. These contributions are estimated by means of 500 MC simulations which use the same pointing and detector weighting as the data, and take into account signal from a beam-convolved realization of a Planck 2015 $\Lambda$CDM power spectrum that includes the effect of gravitational lensing. As in PB14, we model the noise in the simulation as white, and we add random spectrally flat noise to the TOD of each detector variance equivalent to that measured from the detectors. For each detector, the noise variance is estimated from the average of the time domain power spectral densities over 1–3 Hz, corresponding approximately to an $\ell$ range of 500–1500, and the filtering of the TOD is designed to suppress the correlated part of the noise in this frequency range. This choice of using white noise, as opposed to a model of correlated noise, has been validated with the null-test framework (Sec. 4.3).

Finally, we combine the per-patch power spectra into a single patch-combined spectrum according to the band power covariance matrices. The diagonal elements of the patch-combined covariance matrix are used as the pipeline-B error bars shown in Figure 9. Note that the multiplicative uncertainty (due to e.g., the absolute calibration uncertainty), and systematic uncertainties are not shown in Figure 9, but can be found in Table 5. For the released $B$-mode band powers from pipeline A, we go a step further and account for the slightly non-Gaussian shape of the band power likelihoods using an offset log-normal distribution (Bond et al., 2000):

[TABLE]

where $X_{\rm meas}$ is the measure of a parameter of interest $X$, and $X_{\rm cosm}$ a fiducial value for this parameter (the $B$-mode band powers in this case). We fit for the values of $\sigma$ and $x_{0}$ in each band using simulations with 0%, 69% and 100% of the $B$-mode power predicted by the Planck 2015 $\Lambda$CDM model. The intermediate value of 69% was chosen as being close to the observed $B$-mode power. The values are reported in Table 4 along with the band powers. We use the same formalism to empirically model the likelihood for the $B$-mode amplitude parameters $A_{BB}$ and $A_{L}$ described below. We find the average of the upper and lower error bar of the band powers (the $B$-mode amplitude) is few $\%$ smaller or larger ($\sim 7\%$ smaller) than that from using the Gaussian distribution. The shape of the likelihood is positively skewed and the upper error bar is $\lesssim 10\%$ ($\sim 7\%$) larger than the lower error bar.

The combined $C_{\ell}^{TT}$, $C_{\ell}^{EE}$, $C_{\ell}^{TE}$, $C_{\ell}^{TB}$, and $C_{\ell}^{EB}$ spectra after the global calibration of absolute gain and polarization angle are plotted in Figure 8. We find that the patch-combined and individual patch spectra are consistent with the $\Lambda$CDM model. The patch-combined spectra from pipeline A (pipeline B) have a PTE with respect to the $\Lambda$CDM cosmology of 32% (28%), 62%(26%), 92%(67%), 79%(29%), and 60%(75%) for $C_{\ell}^{TT}$, $C_{\ell}^{EE}$, $C_{\ell}^{TE}$, $C_{\ell}^{TB}$, and $C_{\ell}^{EB}$, respectively. The spectra from both pipelines are in good agreement.

The difference of the individual-patch $C_{\ell}^{BB}$ spectra from pipeline A and pipeline B has a PTE with respect to a null spectrum of 23%, where the variance of the null spectrum is derived from an analytical estimate. We remind the reader that the two pipelines treat noisy modes differently (see Poletti et al. (2016) for a detailed discussion), and we therefore expect a larger scatter between results in the noise dominated regime ($\ell>1300$).

The combined $C_{\ell}^{BB}$ spectra are shown in Figure 9, indicating that the spectra from both pipelines are in good agreement. The PTE of these band powers with respect to the Planck 2015 $\Lambda$CDM spectrum are 55% (pipeline A) and 41% (pipeline B). Since the results using pipeline A has satisfied a larger set of the null tests, and the instrumental systematic error analysis was fully performed only on pipeline A, we adopt its power spectra results (tabulated in Table 4) as the reference (official) Polarbear results for this release. The results in Table 4 are the data that should be quoted or used for any further cosmological analysis.

First we consider the significance with which these data rule out the null hypothesis of no $B$-modes. After setting the sample variance to zero, we fit for the amplitude parameter, finding $A_{BB}=0.75^{+0.21}_{-0.20}({\rm stat})^{+0.00}_{-0.04}({\rm inst})$, where in this expression, “stat” refers to the expected statistical fluctuation of the measurement ($X_{\rm meas}$ in Equation 1 with $X=A_{BB}$) evaluated with the likelihood of no lensing signal ($X_{\rm cosm}=0$), and “inst” to the systematic uncertainty associated with possible biases from the instrument. To calculate the upper bound on the additive uncertainties from instrumental systematic errors, we linearly add, in each band, the upper bound of the band powers of all the instrument systematic effects discussed in Sec. 4.4. This produces a lower $A_{BB}$, and sets the lower bound of the additive uncertainty. To be conservative, we evaluate the detection significance by subtracting this systematic uncertainty of 0.04 from 0.75. The likelihood ratio is ${\cal{L}}(X_{\rm meas}=0.71|X_{\rm cosm}=0)/{\rm max}_{X_{\rm cosm}}{\cal{L}}(X_{\rm meas}=0.71|X_{\rm cosm})=9.0\times 10^{-3}$, corresponding to a 3.1$\sigma$ rejection of the no $B$-mode null hypothesis.

Next we fit the band powers to the Planck 2015 $\Lambda$CDM cosmological model with a single $B$-mode amplitude parameter, $A_{BB}$. We find $A_{BB}=0.69^{+0.26}_{-0.25}({\rm stat})^{+0.00}_{-0.04}({\rm inst})\pm 0.04({\rm multi})$, where in this expression, “stat” refers to the 68.3% confidence interval of the estimated quantity ($X_{\rm cosm}$ in Equation 1 with $X=A_{BB}$) given our observation with non-zero lensing signal ($X_{\rm meas}=0.69$), and “multi” to multiplicative calibration uncertainties. The likelihood of the $B$-mode amplitude is shown in Figure 10. The shift in the estimated amplitude with respect to the null result is due to different field and bandpower weighting, which uses the covariance matrix with the Planck 2015 $\Lambda$CDM model, or $A_{BB}=1$.

Finally, we fit for the amplitude parameter $A_{BB}$ after subtracting the foreground terms from Sec. 5, which we denote finally as $A_{L}$. We find $A_{L}=0.60^{+0.26}_{-0.24}({\rm stat})^{+0.00}_{-0.04}({\rm inst})\pm 0.14({\rm foreground})\pm 0.04({\rm multi})$, where “foreground” refers to the total foreground uncertainty. This amplitude can be interpreted as the measured amplitude of the lensing $B$-modes. Table 5 summarizes all the systematic uncertainties in the measurement of $A_{L}$. Although the mean value shifts from the previous measurement in PB14, the change in the mean value is consistent with statistical fluctuation even accounting for the fact that data have overlap between this result and PB14. Compared to PB14, the uncertainty on the lensing amplitude has been reduced by a factor of two thanks to the realistic uncertainty estimate as well as the overall larger dataset and improved calibration.

8. Conclusions

We have measured the CMB $B$-mode angular power spectrum, $C_{\ell}^{BB}$, over the multipole range $500<\ell<2100$ from a blind analysis of data from the first two seasons of Polarbear observations. We have doubled the sensitivity on the lensing amplitude compared to the first-season result (PB14) and rejected the null hypothesis of no $B$-mode polarization with 3.1$\sigma$ confidence. After subtracting the new estimated foreground contamination based on the Planck 2015 data which were unavailable in PB14, we found the amplitude of the lensing signal to be $A_{L}=0.60^{+0.26}_{-0.24}({\rm stat})^{+0.00}_{-0.04}({\rm inst})\pm 0.14({\rm foreground})\pm 0.04({\rm multi})$.

The data were analyzed by two independent pipelines giving results in agreement. These results are supported by an extensive suite of null tests in which 12 divisions of data were used to finalize the data selection and by an estimate of systematic errors from 9 sources of instrumental contamination using a detailed instrument model. We found that all the systematic uncertainties were small compared to the statistical uncertainty in the measurement. To motivate comprehensive evaluation of the dataset and prevent observer bias in data selection and analysis, the selection criteria and systematic errors were determined before the $B$-mode power spectra themselves were examined by the two independent analysis pipelines. In Figure 11, we present our measurement alongside recent measurements by SPTpol (Keisler et al., 2015), actpol (Louis et al., 2016) and Keck Array (The Bicep2 and Keck Array Collaborations: P. A. R. Ade et al., 2015).

We reported previously evidence for gravitational lensing of the polarized CMB in Polarbear data (The Polarbear Collaboration: P. A. R. Ade et al., 2014b) on the same sky area as studied in this work using the PB14 dataset. An updated lensing analysis using the first two seasons of Polarbear observations described in this paper will be presented soon in a separate publication. After the two seasons reported here, we installed a continuously rotating half-wave plate in May 2014 (Takakura et al., 2017), and we are observing a larger patch of the sky ($\sim 700$ $\deg^{2}$) since then, targeting inflationary $B$ modes at degree angular scales. Results from the analysis of this dataset will be reported in future publications.

Foreground contamination is one of the main limiting factors in our measurement of $B$ modes. Neither wmap nor Planck 2015 has enough sensitivity at small angular scales to sufficiently constrain the polarized synchrotron and dust amplitude in the Polarbear observations reported here. Dedicated multi-frequency observations are therefore needed to obtain better constraints. New receivers (Polarbear-2) at $95\,{\rm GHz}$, $150\,{\rm GHz}$, $220\,{\rm GHz}$ and $270\,{\rm GHz}$ with the sensitivity to reach inflationary tensor-to-scalar ratio $\sigma(r)<0.01$ are under development (Inoue et al., 2016) and are being implemented on the Simons Array telescopes (Suzuki et al., 2016). Such improvements will enable the Simons Array to enter into an era of precise CMB polarization measurements, improving our understanding of the early and late Universe physics.

JP acknowledges support from the Science and Technology Facilities Council [grant number ST/L000652/1] and from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement No. [616170]. The Polarbear project is funded by the National Science Foundation under Grants No. AST-0618398 and No. AST-1212230. The James Ax Observatory operates in the Parque Astronómico Atacama in Northern Chile under the auspices of the Comisión Nacional de Investigación Científica y Tecnológica de Chile (CONICYT). This research used resources of the Central Computing System, owned and operated by the Computing Research Center at KEK, the HPCI system (Project ID:hp150132), and the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. In Japan, this work was supported by MEXT KAKENHI Grant Numbers JP15H05891, 21111002, JSPS KAKENHI Grant Numbers JP26220709, JP24111715, JP26800125. This work was supported by World Premier International Research Center Initiative (WPI), MEXT, Japan. This work was supported by JSPS Core-to-Core program. In Italy, this work was supported by the RADIOFOREGROUNDS grant of the European Union’s Horizon 2020 research and innovation programme (COMPET-05-2015, grant agreement number 687312) as well as by the INDARK INFN Initiative. The McGill authors acknowledge funding from the Natural Sciences and Engineering Research Council of Canada and the Canadian Institute for Advanced Research. GF acknowledges support from the CNES postdoctoral program. DB acknowledges support from NSF grant AST-1501422. DB acknowledges support from FONDECYT postdoctoral grant number 3150504. CF acknowledges support from grants HST-AR-13886.001-A, IGPP LANL 368641, NASA NNX16AF39G and Ax Foundation for Cosmology at UC Irvine. CR acknowledges support from a Australian Research Council’s Future Fellowship (FT150100074). MA acknowledges support from CONICYT UC Berkeley-Chile Seed Grant (CLAS fund) Number 77047, Fondecyt project 1130777, DFI postgraduate scholarship program and DFI Postgraduate Competitive Fund for Support in the Attendance to Scientific Events. RD acknowledges support from CONICYT grants FONDECYT 1141113, QUIMAL 160009 and Anillo ACT-1417. ST was supported by Grant-in-Aid for JSPS Research Fellow.

Appendix A Detector gains: modeling Saturn’s rings

For the PB14 data release, we used only observations of Saturn taken between June and September 2012 (4 months of data), and we assumed the temperature of Saturn to be constant, with T${}_{\text{sat}}=148$ K${}_{\text{RJ}}$. However, the analysis of Saturn observations over a longer period (May 2012 – April 2014) exhibited variations up to 10% over time (Peloton, 2015). We investigated several possible explanations, such as a misestimation of the apparent diameter of the planet, dependency in the elevation, miscorrection of the atmosphere contribution, stability of the stimulator, change in the main beam parameters, or even changes in the sidelobes over time, but none gave an explanation for the shift of approximately $10\%$ in between the two seasons. This shift was however explained by taking into account the influence of Saturn’s rings on the planet flux. Even if we cannot resolve the rings,101010Given the angular size of Saturn with respect to our beam size, we cannot resolve the details of the planet. The variations in the atmosphere of Saturn (pole/equator differences, clouds, seasonal variations, etc.) are also not taken into account here. their inclination with respect to the plane of observation from the Earth produces a variation over time in the measured microwave brightness temperature, as described and shown in, e.g., Dunn et al. (2002), Weiland et al. (2011), and Hasselfield et al. (2013). The contribution of the rings to the total temperature brightness of Saturn is twofold: they obscure and therefore reduce the flux from the main disk, but they also radiate at a lower temperature and contribute to the total signal through a mix of scattering and thermal re-emission of planetary radiation (at 150 GHz, the thermal emission dominates over the scattering). The change in opening angle (i.e., the angle made by the plane of the rings and the line of sight of the observer) allows us to break the degeneracy between the contribution of the disk and that of the rings. We follow the two-parameter empirical model proposed by wmap (Weiland et al., 2011), by modeling the total brightness temperature of Saturn as a contribution from the disk $T_{\text{disk}}$ with seven extended surrounding rings described by one global effective temperature $T_{\text{ring}}$.111111$T_{\text{disk}}$ and $T_{\text{ring}}$ temperatures are assumed time invariant. The individual contributions from those two parameters to the total brightness is given by the orientation of the system at the moment of the observation.

We use 132 observations of Saturn between May 2012 and April 2014, selected for their high quality of data, and spanning ring opening angles from $\sim 12^{\circ}$ to $\sim 23^{\circ}$. The observations are combined into periods of 15–30 days, taking the uncertainty in each combined measurement to be the error in the mean of the contributing data. The model fit parameters and uncertainties ($\chi^{2}\sim 2$ per degree of freedom) are $T_{\text{disk}}=126.8\pm 2.7$ K${}_{\text{RJ}}$ and $T_{\text{ring}}=16.3\pm 6.0$ K${}_{\text{RJ}}$ (after applying the absolute calibration derived from the CMB temperature spectrum as described in Sec. 6). The uncertainties include a statistical contribution and a systematic contribution from the wafer-to-wafer variations. act collaboration provides complementary measurements at 148 GHz for ring opening angles from $-2^{\circ}$ to $12^{\circ}$ in Hasselfield et al. (2013). They found $T_{\text{disk}}=133.8\pm 3.2$ K${}_{\text{RJ}}$ and $T_{\text{ring}}=17.7\pm 2.2$ K${}_{\text{RJ}}$, which is consistent with our values. On the other hand Planck 2015 collaboration reported recently a measurement of Saturn’s disk temperature rescaled at 147 GHz (Planck Collaboration et al., 2016e). They found a temperature of the main disk $T_{\text{disk}}=145.7\pm 1.1$ K${}_{\text{RJ}}$ higher than the Polarbear and act values. However, Figure 10 in Planck Collaboration et al. (2016e) seems to suggest that the ring temperature between Planck 2015 and act, and therefore Polarbear, are in agreement.