Improved Upper Limit on Degree-scale CMB B-mode Polarization Power from the 670 Square-degree POLARBEAR Survey

We report an improved measurement of the degree-scale cosmic microwave background B-mode angular-power spectrum over 670 deg2 sky area at 150 GHz with Polarbear. In the original analysis of the data, errors in the angle measurement of the continuously rotating half-wave plate, a polarization modulator, caused significant data loss. By introducing an angle-correction algorithm, the data volume is increased by a factor of 1.8. We report a new analysis using the larger data set. We find the measured B-mode spectrum is consistent with the ΛCDM model with Galactic dust foregrounds. We estimate the contamination of the foreground by cross-correlating our data and Planck 143, 217, and 353 GHz measurements, where its spectrum is modeled as a power law in angular scale and a modified blackbody in frequency. We place an upper limit on the tensor-to-scalar ratio r < 0.33 at 95% confidence level after marginalizing over the foreground parameters.


INTRODUCTION
Anisotropies in the cosmic microwave background (CMB) bring us fundamental information about our universe.If detected, degree-scale B-mode polarization, the parity-odd component of the linear polarization anisotropies, is a footprint of the primordial gravitational waves generated during the cosmic inflation era.By measuring the amplitude of the B-modes, we can determine the tensor-to-scalar ratio r and test the physical mechanisms of the inflation.
The Polarbear experiment is a ground-based experiment in the Atacama desert in Chile.It consists of the 2.5 m aperture Huan Tran Telescope with 1274 transition-edge-sensor bolometers sensitive to the Corresponding author: S. Takakura satoru.takakura-1@colorado.edu 150 GHz band (Kermish et al. 2012;Arnold et al. 2012).In 2014, we installed a continuously rotating half-wave plate (HWP) at the prime focus (Takakura et al. 2017).The HWP modulates incoming linear polarization signals and therefore reduces low-frequency noise due to both the atmosphere and the instrument.
In Polarbear Collaboration (2020, hereafter PB20), we reported a measurement of the degree-scale B-mode angular-power spectrum using data from three years of observations from 2014 to 2016.The HWP modulation results in a relatively low knee in the noise spectrum at knee = 90, where the contribution of the low-frequency noise to the power spectrum uncertainty becomes comparable to that of detector white noise.We place an upper limit of r < 0.90 at the 95% confidence level.
In PB20, however, we used only 29.2% of data after eliminating data from detectors that failed to tune correctly or that have glitches due to various disturbances.Data containing glitches can potentially be made available by improvements to the analysis process.We find that most of the glitches in the detector polarization timestream come from an angle error of the HWP.Here, the angle error is the offset of the measured angle from the real angle, which occasionally occurs due to electrical noise within the encoder circuit.By improving the The main improvement in this study is in the processing of the encoder data of the HWP angle.In this section, we explain details of the HWP angle error: how the angle error causes a glitch on the detector signal, how the angle error is caused, and how we have improved the correction of the angle error.
Except for the improved correction of the HWP encoder data, we follow the data processing presented in PB20 and references therein.

Glitches Due to the Angle Error of the HWP
Glitches are spurious signals in detector timestreams.They have several causes: transient physical events such as cosmic-ray hits, atmospheric noise in bad weather conditions, electrical noise pickup, and unexpected data drop in the readout system.Thus, glitches have various timescales and shapes.To drop all kinds of glitches, we apply several filters in our analysis process.
We apply glitch detection for three types of timestreams: full-sampling timestreams for each detector, demodulated and downsampled timestreams for each detector, and timestreams averaged among all detectors.In the first step, we detect short-timescale glitches such as cosmic ray hits.In the next step, we focus on glitches below 4 Hz after demodulation,1 which contaminate our science signal.Electrical pickup and bad weather data are flagged.Finally, we catch faint but correlated glitches.Polarized bursts due to clouds (Takakura et al. 2019) are detected here.The commonmode glitch detection has the largest impact on the data selection because it has the highest sensitivity and affects all detectors.
The angle error of the HWP is another source of correlated glitches.In observations with the rotating HWP, the detector signal, d(t), is modeled as (Takakura et al. 2017) The unpolarized Stokes component, I(t), is not modulated, while the linear polarization components, Q(t) and U (t), are modulated by the angle of HWP, θ(t) = ωt, where ω/2π = 2.0 Hz.N (t) is detector noise.The last term is instrumental signals called HWP synchronous signals (HWPSSs), which are classified by the order of the harmonic n.We obtain polarization components by demodulating d(t) as where F LP and F BP are low-pass and bandpass filters used to select signals around the modulation frequency 4ω.Here, we use the measured angle of the HWP θ (t), which is reconstructed from the HWP encoder data.If the measured angle has an error from the actual angle as ∆θ(t) = θ (t) − θ(t), the demodulated signal becomes where N Q (t) and N U (t) are demodulated detector noise, A 4 (t) is the fourth harmonic of the HWPSS, and A 4 is its average.Here, we assume that ∆θ(t) 1 and A 4 is much larger than Q(t) and U (t).In the case of Polarbear, instrumental polarization due to the primary mirror produces A 4 ∼ 0.1 K uniformly among all detectors (Takakura et al. 2017).Therefore, the last term of Equation (3) becomes a source of correlated glitches.The blue line shows the interpolated angle using a timestamp.The wavy pattern on it is the actual rotation jitter of the HWP.In the third rotation, it jumps by 0. • 5 due to an error count and is reset by the indexing signal.The orange line shows the correction in PB20, which linearly interpolates the jump.The green points show the new method, in which we fix the error at the encoder count level and then apply timestamp interpolation.The green line is the result, which keeps the actual HWP angle jitter.The bottom panel shows the detector signal of Stokes U component in the instrumental coordinates.PB20 data, the orange line, show a glitch, which is consistent with the expectation by Equation (3), the black line.In the new data, the green line, the HWPSS is subtracted correctly, and the glitch disappears.

Example of Data with the Angle Error
Figure 1 shows an example of HWP encoder data causing the angle error ∆θ(t) and the resulting glitch in the detector timestream.To understand this, it is necessary to explain how the encoder works.
We measure the angle of the HWP using an encoder plate with 360 precisely machined teeth and one index hole.Optocouplers regularly sample whether the gate is open or closed at 40 kHz.Synchronizing with detector sampling at 191 Hz, we store the count of the edges and the timing of the last edge.Since the timings of the edge and detector sampling are asynchronous, the raw encoder count (the blue points) has a quantization noise of 0. • 5. We fix this quantization error by interpolating the middle angle at the detector sampling between edges using the timestamp information (the blue line).The statistical uncertainty of the angle (the width of the blue line) comes from the quantization error on the timestamp of 3×10 −6 rad √ s.This uncertainty causes an angle error noise of 1.4 µK √ s, which is smaller than the detector array sensitivity.On the other hand, the HWP rotation is not perfectly continuous and contains a jitter of δθ(t) = θ(t) − ωt.The wavy fluctuation of the blue line is the jitter, which has a weak resonance around 8 Hz probably due to the combination of the spring constant of the system and feedback parameters for the servo motor driver.This actual jitter does not introduce noise on the demodulated signal if the encoder measures the jitter as θ (t) = θ(t) = ωt + δθ(t) and ∆θ(t) = 0.
The problem in our HWP encoder is that it occasionally has erroneous counts due to electrical noise (the jump in the blue line).We detect these bad counts by comparing two optocouplers.In PB20, we dropped encoder samples from the bad counts to the next index signal and interpolated linearly (the orange line).This interpolation nicely tracks the continuous rotation, but not the jitter, i.e., θ (t) = ωt.The angle error in the interpolated samples becomes ∆θ(t) = δθ(t).As shown in the bottom panel of Figure 1, this angle error causes the glitches (the orange points) by the last term of Equation (3) (the black line).Here, the amplitude of the glitch is ∼ 1 mK, which is comparable to the instantaneous sensitivity of a single detector.This means that we cannot detect this glitch in a single detector analysis.However, this glitch is correlated among detectors, but the white noise is not.Therefore, averaging hundreds of detectors improves the significance dramatically.
We improved the correction method to directly decrement the counter in the offline analysis (the green points).Then, we apply quantization noise reduction as other normal data (the green line).In this method, the fixed angle keeps information of the actual jitter as θ (t) = ωt + δθ(t), and the angle error ∆θ(t) becomes zero.Therefore, we can clean the glitches as the green points in the bottom panel of Figure 1.
We reprocessed all data with the new method and successfully cleaned this type of glitch as described in Section 3.
To mitigate the risk of this type of error in HWP data acquisition in future experiments such as the Simons Array experiment (Suzuki et al. 2016), there are two solutions.One simple way is making the hardware for the encoder data acquisition robust.In PB-2a, the first receiver of the Simons Array, we use a commercial encoder (Hill et al. 2016).The other solution is making the rotation stable, which makes the data more robust for offline data correction even if there were some faulty data.The cryogenic superconducting bearing technique developed for PB-2b, the second receiver of the Simons Array, is promising to achieve very stable rotation (Hill et al. 2020).

DATA SELECTION
We applied the same PB20 data selection method to the data processed with the new encoder correction.The results of data volume and data selection efficiencies are summarized in Tables 1 and 2. The final volume of data available increased by 81.6% from PB20.
The total observation time has slightly increased from PB20 by 3.1%.We use observations from 2014 July 25 until 2016 December 30.We also recovered some missing observations by reconstructing the database of all observations.
The detectors and their calibrations (pointing offset, relative gain, and relative polarization angle) used in this analysis are identical to those in PB20.Thus, we have no increase in detector yields.
The main data increase comes from improvements in data selection efficiency.To explain how the encoder correction affects the data selection, we briefly explain our data selection procedure.
The data selection is done by flagging bad data whose data quality exceeds some threshold.We use various types of data qualities evaluated using detector data, calibration data, housekeeping data, and external data.Here, some data quality evaluations require data selection based on other low-level data qualities.For example, bad sections of detector timestreams with glitches are masked in the evaluation of the power spectrum density (PSD) to prevent spurious contamination.
The data quality items and the requirements in the new analysis are almost the same as in PB20.One property added is the fraction of remaining subscans.Our one-hour observation consists of about 70 scans left and right at constant elevation.A subscan is each one-way scan.The minimum unit of our data selection flagging is each subscan.Sometimes, most of the subscans are flagged, which causes a problem in the polynomial filter to detrend the baseline drift of the detector.In PB20, such data are removed in the timestream filtering step on the fly, and flagged as one of the reasons of failure in the evaluation of PSD.In the new analysis, it is explicitly captured in the data selection.
Even with the same criteria for the data selection, the efficiencies improve for good-quality data.The glitches due to the encoder error are detected in the commonmode glitch stage.Here, to test the Gaussianity of the averaged detector timestream d t , we compute its kurto- − 3 and require −1.5 < K < 10 as shown in Figure 2. Thanks to the encoder correction, the glitches in the data have disappeared, and thus the efficiency of the step has significantly improved from 60.6% to 92.3%, which means an increase in the remaining data by 52.4%.
In addition, the selection efficiency based on the remaining subscan fraction has also improved by 22.6%.Sometimes, the encoder error occurs so often that most of subscans in a one-hour observation are flagged, then we have to drop the observation entirely.As the former selection efficiency based on glitches increases, the efficiency of this following data selection also improves.

ANALYSIS PIPELINE
Here, we analyze the new data set and measure the angular-power spectrum of the CMB B-mode at the multipole range of 50 < < 600. Figure 3 shows the overview of the analysis pipeline.All the details are described in PB20.Here, we only report updates on the data validation, overall calibrations, and measured power spectra.

Calibrations Based on the Power Spectrum
As described in PB20, we perform calibrations of the instruments in two steps.In the first step, we use calibration measurements, e.g., the thermal source calibrator and observations of planets and other bright sources, including Tau A. We calibrate the pointing model of the telescope, the pointing offsets of detectors, the effective beam function, relative gain variations, detector time constants, relative polarization angles, and polarization efficiencies.All these calibrations are the same as in PB20.In the second step, we determine overall calibration parameters using the measured CMB angularpower spectra.We calibrate the absolute gain and the beam uncertainty using the E-mode autospectrum, and the absolute polarization angle using the EB crossspectrum.
The calibration error on power spectra due to the overall gain g 0 and beam uncertainty σ2 is modeled as where the exponential term represents Gaussian smearing due to the pointing jitter.Here, we treat σ 2 as a single parameter and allow it to be negative, which helps to remove a potential systematic error in our beam calibration.We estimate g 0 and σ 2 using the Planck 2018 PR3 143 GHz full mission maps.We compute Polarbear-auto, Planck-auto, and Polarbear-cross-Planck E-mode binned power spectrum estimates, The Polarbear autospectrum is calculated by crosscorrelating 38 submaps grouped every 10 days, thus free from the noise bias.The uncertainties of the spectra are estimated using the noise simulations.
We fit Equation (5) varying g 0 , σ 2 , and DEE b .We find the overall gain calibration factor of g 0 = 1.106 ± 0.021 and the beam uncertainty factor of σ 2 = 0.99 ±  3.00 arcmin 2 .These are in good agreement with PB20, as expected.
We calibrate the absolute polarization angle using the observed EB cross-spectrum.We assume that the original EB correlation is null, but the observed one has finite values from the leakage of E-modes due to the angle error ψ as (Minami et al. 2019) Here, we use the observed E-mode spectra instead of the theoretical one used in PB20.It makes the result independent of the absolute gain calibration and the fiducial cosmological parameters.The uncertainty of the spectrum ∆ DEB b,PB is estimated using the quasi-analytic method (PB20) as where N EE b,PB is the noise bias of Polarbear autospectrum estimated using the noise simulations, and ν EE b,PB is the degrees of freedom, estimated as Here we assume that the noise bias and degrees of freedom for B-modes are similar to E-modes.

Data Validations
As described in Section 2, the main difference of the new data set from PB20 is the correction of the glitch due to the HWP angle error.It should not inject any systematic uncertainties.Therefore, we use the same estimates of the systematic uncertainties done in PB20.However, thanks to the increase of statistics, unknown systematics may become noticeable.Thus, we perform the null tests in PB20 with the new data set.
The null tests are performed for the same 18 splits as PB20.The null spectra are compared to 500 noise-only simulations generated using the "signflip" method.We estimate the uncertainty of the null spectra from the noise simulations and evaluate a χ value for each spectrum, each null split, and each bin.We compute the sum of χ 2 over bins or over null splits.Then we evaluate the probability to exceed (PTE) values by counting the fraction of the noise simulations whose total χ 2 exceeds the value of real data.Table 3 shows the results.Next, we compute five representative statistics: (1) the average of χ among all bins and all null splits, (2) the most extreme total χ 2 by bin summed over null splits, (3) the most extreme total χ 2 by test summed over bins, (4) the most extreme χ 2 among all bins and all null splits, and (5) the total χ 2 summed over bins and null splits.The PTEs are computed by comparing the statistics from real data with those evaluated from every realization of the noise simulations.Finally, we choose the lowest PTE of the five statistics and evaluate its  PTE again by comparing with the lowest PTEs from the noise simulations.Table 4 shows the results.In the EE spectrum, for example, the lowest PTE is 3.6%, but in 14.0% of the noise simulations, the lowest PTE becomes lower than 3.6%.Thus, we obtain the final PTE of 14.0%.In addition, we test the consistency of the distribution of PTE estimates with a uniform distribu-  4. Finally, we check the power spectra except for B-mode autocorrelation.As absolute calibrations in Section 4.1, we use Planck 143 GHz maps as reference.Figure 4 shows the power spectra from Polarbear and Planck maps.We take the inverse-variance weighted average3 of the three or four auto-and cross spectra and compute the total χ 2 for each spectrum compared to the averaged one.Table 5 shows the PTE values of the total χ 2 as well as that of the overall total χ 2 .We find that all these spectra are consistent between Polarbear and Planck.Here, we apply the overall calibrations, which use EE and EB spectra as Section 4.1.The consistencies of T T , T E, and T B spectra support the robustness of our calibrations and analysis methods.

Results of B-mode Power Spectrum Estimates
The results of B-mode power spectrum measurements with the new data set are summarized in Table 6 and shown in Figure 5. Again, we estimate the statistical uncertainties using the quasi-analytic method with 500 noise-only simulations based on the "signflip" method.The estimates of systematic uncertainties are same with PB20 because we use the same period of observations and the same instrumental calibrations.The overall calibration uncertainties are multiplicative error due to the overall gain and beam calibrations in Section 4.1.
We compare our B-mode measurements with a model of B-mode signal based on the Planck 2018 ΛCDM lensing B-mode spectrum and a foreground component modeled by BICEP2/Keck Collaboration & Planck Collaboration (2015).Direct comparison gives the reduced χ 2 of 9.34/11, indicating good agreement.Fitting an overall amplitude rescaling this template, we ) )  are uncertainties on the band powers due to noise statistics, instrumental systematics, and overall calibration uncertainties, respectively.We assume the same estimates of ∆ DBB b,syst as PB20.
find A BB = 0.59 +0.46  −0.31 , and the null hypothesis is disfavored at 1.4 σ.Note that this significance is lower than the 2.2 σ in PB20, but the uncertainty on A BB has improved from 0.8 thanks to the increase of data volume.See Section 6 for more detailed discussions about the consistency with PB20.
In Section 5 we perform more detailed parameter estimations combining Planck 2018 PR3 observations at 143, 273, and 353 GHz.As we do for the 143 GHz map used for the absolute gain calibration in Section 4.1, we emulate Polarbear observations on the Planck maps, and process map-making as for Polarbear data.We stack maps from all emulated observations by frequency, and compute auto-and cross spectra from these Planck three-frequency maps as well as from the Polarbear map.Here, auto-spectra from the same Planck observation contain the noise bias, which is estimated using the 96 Planck noise simulation maps for each frequency.

PARAMETER ESTIMATION
In this section we briefly explain our likelihood, the constraints on the tensor-to-scalar ratio r, and the foreground contamination in our BB spectrum.We follow the same procedure as PB20, which we summarize briefly here.Our estimation uses BB signal and noise spectra of four auto-and six cross spectra from Polarbear 150 GHz, Planck 143 GHz, Planck 273 GHz, and Planck 353 GHz.We arrange these 10 measured signal spectra in such a way that each bandpower is a matrix Db,ij , where i and j are one of the four observations.The diagonal block of this matrix contains four autospectrum values, and the off-diagonal block contains six cross-spectrum values.The noise spectra are also arranged as N b,ij .Since noise between frequency bands and experiments is expected to be uncorrelated, N b,ij is diagonal.
The signal spectra D b,ij can be modeled as a sum of CMB and foregrounds as where D fg b,ij is the total foreground component, which depends on frequencies.The CMB component D CMB b is the same among the frequencies and has contributions from primordial B-modes and from the weak lensing of the CMB.Thus it can be modeled as where r is the tensor-to-scalar ratio, and A lens is the normalized amplitude of the ΛCDM lensing signal.D tens b and D lens b are the binned tensor and lensing signals, respectively.The tensor spectrum at r = 1 is computed using CAMB (Lewis et al. 2000) with a tensor power amplitude of 2.46 × 10 −9 and a spectral index of zero.As shown in PB20, the upper limit of contamination from the polarized sources and synchrotron emission is subdominant to the dust component at 150 GHz.So we assume that our measured spectra Db,ij are contaminated by foregrounds mostly from Galactic dust, D fg b,ij ≈ D dust b,ij , where D dust b,ij is the dust component, which we consider a power law in and a modified blackbody in frequencies i, j as defined in Adam et al. (2016) and Akrami et al. (2020), Here w b is the bandpass window function.A dust is the amplitude of the dust component, and α dust is the power-law index in .We consider a pivot value of 0 = 80.The f i is the dust emission at each frequency bandpass in CMB temperature units defined as f (β dust , T dust ), where β dust is the spectral index, and T dust is the temperature of the modified blackbody.The f i is normalized such that A dust corresponds to the dust emission at 353 GHz.

Likelihood
Under the assumption according to Hamimeche & Lewis (2008) that the measured Dtot b,ij = Db,ij + N b,ij follows a Wishart distribution (Wishart 1928) with an effective number of degrees of freedom ν b , we define our likelihood L of the true spectrum (12) The effective number of degrees of freedom ν b is estimated from the standard deviation of the spectrum of the noise realizations.The standard deviation in autoand cross-spectrum is determined following PB20 as For our estimation, we use the geometrical mean of ν b of Polarbear and Planck.We sample our likelihood using EMCEE (Foreman-Mackey et al. 2013) for the parameter estimation.Our model contains four free parameters, r, A dust , α dust , and β dust .We fixed the values of A lens = 1 and T dust = 19.6K (Planck Collaboration 2015).For α dust and β dust , we considered Gaussian priors −0.58 ± 0.21 and 1.59 ± 0.11, respectively (BICEP2/Keck Collaboration & Planck Collaboration 2015; Adam et al. 2016).

Constraints on Parameters
Marginalized 68% and 95% parameter constraint contours and the posteriors are shown in Figure 6.Similar to PB20, the posteriors of α dust and β dust are dominated by the priors because these parameters are much less sensitive to the Polarbear data.Our estimate excludes the zero dust foregrounds with 99% confidence and shows evidence of dust B-modes with amplitude A dust = 4.0 +1.7  −1.6 .The parameter A dust is mostly constrained by Planck data.The 10% increase in the bestfit value compared to PB20 is due to the degeneracy with the parameter r.We find a 68% confidence level maximum likelihood value of r = −0.04+0.18 −0.15 (stat) ± 0.03(sys).We report the improved 95% confidence upper limit of r < 0.33 after marginalizing over foreground parameters, requiring r and A dust to be positive posteriori.The new addition of data tightens the constraint on r by a factor of 2.7 compared to PB20.We validate the constraint in Section 6.

Goodness of Fit
As a measure of the goodness of fit of D b,ij to Db,ij , we can define an effective chi-square following Hamimeche & Lewis (2008), Here, χ 2 eff is zero if D b,ij = Db,ij .For ν b n freq and in the limit of a negligible number of fit parameters compared to the total number of bins across all spectra, the expectation value and variance of the effective chi-square under the Wishart distribution is given by where n bins is the number of multipole bins of the spectra.In Figure 7 we show our maximum likelihood χ 2 eff = 122.98 with the red line.The solid vertical black line shows the expected value, and the gray shaded area shows the variance.The χ 2 eff of our data is consistent with the simulated χ 2 eff , and it lies within the expected variance under the Wishart distribution.In Figure 8 we show the normalized difference between the measured cross spectra and the best-fit CMB+foregrounds model.

COMPARISON WITH THE RESULTS OF PB20
Here, we compare the new results with those of PB20 and discuss their consistency.Since both data sets pass the null tests, the results are statistically valid.However, the new data set significantly overlapsthe PB20 data set, and the change should be the result of the additional data recovered in Section 3.
In this section, we evaluate the probability of the consistency using a Monte Carlo (MC) simulation as follows.First, we evaluate the effective data increase between PB20 and the new data set.Next, we simulate expected shifts of measurements from PB20 due to the increase of statistics.Finally, we evaluate the probability by comparing the actual new measurement with the MC simulations.First, we evaluate the effective data increase by comparing the noise power spectrum, N b .The estimate of fractional data increase of 81.6% in Section 3 assumes uniform weights among detectors and observations.In practice, however, each observation has a different noise performance depending on the weather and other instrumental conditions.In the map-making and coadding steps, we average data using the inverse-variance weight-ing.If the recovered data have a better noise performance than the PB20 data set, the effective data increase can exceed the naive estimation.
In the analytic power spectrum uncertainty estimation known as the Knox formula (Knox 1995), the inverse of the noise power spectrum is proportional to the sensitivity-weighted data volume per solid angle as where n det is the total number of detectors, t obs is the total observation time, N NET is the instantaneous sensitivity of a single detector, 4πf sky is the solid angle of the observing patch, and the exponential term is the Gaussian beam smearing effect.We generalize this formula and apply it to our binned noise power spectrum, 2π w b N .Since the observing patch and the scan strategy are the same, the patch size f sky and the mode mixing matrix w b are the same between PB20 and the new data set.Therefore, the ratio of N b between PB20 and the new one corresponds to the ratio of the total data volume, n det t obs , weighted by N −2 NET .We obtain +106.5% for the lowest bin and +95.0% on average.

Simulation of the Power Spectrum with Additional Data
Next, we perform an MC simulation to compute the expected measurements of power spectra by increasing statistics from PB20.Here, we do not simulate detector timestreams or maps, but directly simulate power spectra assuming their statistical behavior.We simulate a set of correlated first and second power spectrum measurements so that their statistical uncertainties are equal to PB20 and the new measurements.We select simulations whose first measurement is the same as PB20.Then, the second measurements of these simulations are the statistical expectations of the new measurement.The details of the calculation are explained in Appendix A.
Figure 9 shows the results from 100,000 simulations for Polarbear autospectrum and cross spectra with Planck 143, 217, and 353 GHz.Here, the MC simulations of the second measurement have shifts from PB20, even though we require the first measurement to be PB20.This is because the second measurement tends to approach the true signal we assume, which is the lensing B-mode with Galactic dust foreground and r = 0.

Consistency Probabilities
We compare the new result in this work with these MC simulations and evaluate PTEs for each spectrum and each band.The results are shown in Table 7.We find no extremely low PTEs in the Polarbear autospectrum and cross-spectrum with Planck 143 GHz.
In cross spectra with Planck 217 and 353 GHz, we find some low PTEs.Table 8 also shows the total χ 2 of 11bins and its PTE.We find again that the cross-spectrum with Planck 353 GHz has a low PTE.The overall χ 2 of all 44 bins is 63.6, and its PTE is 3.0%.
We also perform the KS test to compare the distribution of PTEs with a uniform distribution.The results for each spectrum are shown in Table 8.The overall KS test PTE from all 44 PTEs in Table 7 is 44.0%.In this case, we do not find any low PTEs.
Although the overall χ 2 PTE is low, it is due to the two highest points of the cross-spectrum with Planck 353 GHz.The positive shift at 550 < < 600 may indicate the contamination of polarized dusty star-forming galaxies (Bonavera et al. 2017;Lagache et al. 2020), although this would not explain the negative shift at 500 < < 550.If we remove the highest point, the total χ 2 of the spectrum becomes 15.7, and its PTE improves to 11.0%.The overall χ 2 PTE also improves to 12.4%.The impacts of the two lowest PTE bins on the tensor-to-scalar ratio r are smaller than 2 × 10 −4 .Since all PTEs for Polarbear autospectrum are reasonable, we conclude that the PB20 result and the new result are consistent in terms of the r measurements.Future measurements with more statistics will give a conclusive understanding of the cross-spectrum with Planck 353 GHz.
Finally, we apply this consistency check to the evaluation of the tensor-to-scalar ratio r.Since lower bins are more sensitive to r, the overall PTE changes.We compare the estimate of r from the new measurement with those from the MC simulations.Note again that MC simulations assume a fiducial true signal with r = 0.Only for consistency checking, we use a naive estimation compared to that performed in Section 5. We fit only r by fixing the foreground parameters.We obtain the distribution of the best-fit r from MC simulations as r = 0.08 +0.13  −0.12 .Here, the bias comes from the contribution of the first measurement fixed to the result of PB20.By comparing the actual best-fit r = −0.04 from the new data set, we obtain the PTE of 32.9%.

CONCLUSION
We perform an improved analysis of the PB20 data, the three seasons of Polarbear observations on a 670 deg 2 patch of the sky.We successfully recover 80% more data by improving the angle error correction of the rotating HWP.
By processing the data using the same analysis pipeline as was used in PB20, we measure the CMB B-mode power spectrum at the multipole range 50 < < 600.We find no excess signal beyond that expected from the combination of a ΛCDM model and a Galactic dust foreground.We place an upper limit on the tensor-to-scalar ratio r < 0.33 at a 95% confidence level.The change in the Polarbear B-mode power spectrum from PB20 is consistent with statistical expectations due to the additional data.
Our result demonstrates the possibility of measuring degree-scale CMB anisotropies with a ground-based telescope located in the Atacama desert of Chile.The rotating HWP is a key technique for separating contamination from atmospheric fluctuations and achieving a good noise performance at low frequencies.The low-frequency noise performance enabled by the continuously rotating HWP is one of reasons why some future experiments, including the Simons Array experiment (Suzuki et al. 2016) and the Small Aperture Telescope of the Simons Observatory experiment (Ade et al. 2019), plan to employ rotating HWPs.As we show in this paper, accurate angle encoding is a key to the success of this methodology.

Figure 1 .
Figure1.Example of encoder data with a wrong count and corresponding detector signal.The horizontal axis is the rotation count of the HWP, which takes 0.5 s per rotation.The top panel shows HWP encoder data subtracted by the regular drift, ωt.The blue points are raw encoder counts sampled at 191 Hz.The uniform pattern is quantization noise due to the resolution of the encoder plate by 0. • 5.The blue line shows the interpolated angle using a timestamp.The wavy pattern on it is the actual rotation jitter of the HWP.In the third rotation, it jumps by 0. • 5 due to an error count and is reset by the indexing signal.The orange line shows the correction in PB20, which linearly interpolates the jump.The green points show the new method, in which we fix the error at the encoder count level and then apply timestamp interpolation.The green line is the result, which keeps the actual HWP angle jitter.The bottom panel shows the detector signal of Stokes U component in the instrumental coordinates.PB20 data, the orange line, show a glitch, which is consistent with the expectation by Equation (3), the black line.In the new data, the green line, the HWPSS is subtracted correctly, and the glitch disappears.

Figure 2 .
Figure 2. Distribution of kurtosis computed for the common mode of the detector timestream.The blue line shows PB20 data, and the red line shows new one with the encoder correction.The dashed lines are upper and lower thresholds for the data selection.When the timestream follows the normal distribution, the kurtosis becomes 0. Glitches make the kurtosis > 0 as in the PB20 data.Kurtosis < 0 is due to residual low-frequency noise.
2 DEE b,PB , DEE b,P143 , and DEE b,PB×P143 .The calibration errors in the Polarbear maps modify the spectra as true E-mode spectra in the Polarbear patch, ∆ DEE b is the statistical uncertainty for each spectrum, and ĝb = w b g is the binned calibration factor weighted by the bandpass window function w b .Here, we emulate Polarbear observations on the Planck maps and process the same map-making as Polarbear data.The resulting Planck map should contain the same signal DEE b as Polarbear.Noise bias in the Planck autospectrum is estimated and subtracted using 96 realizations of Planck noise simulation maps.

Figure 3 .
Figure3.Schematic view of the steps in our power spectrum estimation pipeline.Each box represents an analysis procedure, input data, or its product.The arrows show the flow of the pipeline, where the dotted arrows show the flow for the null test especially.All the details of the pipeline are described in PB20.The main change in this work is the update of the input Polarbear data to the new data set, as described in Section 2. The other procedures of the pipeline are almost the same as in PB20, except for the minor change in the overall calibration, as described in the text.

Figure 4 .Figure 5 .
Figure 4. Binned estimates for the angular-power spectrum for T T , T E, EE, T B, and EB spectra from the top to the bottom panel.The red and blue points show the spectra from Polarbear maps and Planck 143 GHz maps.The orange, green, and cyan points show the spectra from different combinations of cross correlation between Polarbear and Planck.The error bars include the sample variance.Gray lines show the theoretical estimates with our fiducial cosmological parameters, and the black points are their binned averages.

Figure 6 .
Figure6.Marginal posteriors for the four parameters r, A dust , α dust , and β dust .We compare our previous estimate (red) with the reanalysis estimate (blue).

6. 1 .Figure 7 .Figure 8 .
Figure 7. Effective chi-squared of the best fit as the red vertical line.The black line is the expectation value of the Wishart distribution, Equation (15), and the shaded area within the two dotted black verticals is the variance of the Wishart distribution, Equation (16).The histogram in red shows the effective chi-square distribution for a set of 96 simulations.

Figure 9 .
Figure 9. MC simulations of the second measurement from PB20 with additional data.The Polarbear autospectrum and cross spectra with Planck 143, 217, and 353 GHz are shown from left to right.The orange points show the median value and 1σ range of the MC simulations.The blue and red points are the central values of DBB b in PB20 and in this work, respectively.The gray line shows the signal spectrum assumed for the MC simulations.

Table 1 .
Data Volume

Table 2 .
Data Selection Efficiency Note-Since there are some updates in data selection, we cannot directly compare data selections in this work and PB20.Here, the PB20 data selection is reproduced applying the new data selection to PB20 intermediate data just for comparison, thus the cumulative efficiency is different from Table1.

Table 3 .
Null Test Total χ 2 PTE Values

Table 4 .
Null Test PTE Values

Table 5 .
PTE Values of the Planck Consistency Check Kolmogorov-Smirnov (KS) test.We test distributions of PTEs in Table 3, as well as PTEs for individual χ values.The results are shown in Table 4.We require that the PTE values of the lowest statistic and KS tests must be greater than 5%.All spectra (EE, EB, and BB) pass these criteria, as shown in Table

Table 6 .
B-mode Band Powers and Uncertainties b is the binned angular-power spectrum estimates from the Polarbear maps.∆ DBB b,stat , ∆ DBB b,syst , and ∆ DBB b,cal

Table 8 .
Comparison PTE for Each Spectrum PB × PB 143 × PB 217 × PB 353 × PB The number in parentheses shows the total χ 2 value.