Issue 
A&A
Volume 596, December 2016



Article Number  A107  
Number of page(s)  52  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201628890  
Published online  12 December 2016 
Planck intermediate results
XLVI. Reduction of largescale systematic effects in HFI polarization maps and estimation of the reionization optical depth
^{1} APC, AstroParticule et Cosmologie, Université Paris Diderot, CNRS/IN2P3, CEA/lrfu, Observatoire de Paris, Sorbonne Paris Cité, 10 rue Alice Domon et Léonie Duquet, 75205 Paris Cedex 13, France
^{2} African Institute for Mathematical Sciences, 68 Melrose Road, Muizenberg, 7945 Cape Town, South Africa
^{3} Agenzia Spaziale Italiana Science Data Center, via del Politecnico snc, 00133 Roma, Italy
^{4} Aix Marseille Université, CNRS, LAM (Laboratoire d’Astrophysique de Marseille) UMR 7326, 13388 Marseille, France
^{5} Aix Marseille Université, Centre de Physique Théorique, 163 Avenue de Luminy, 13288 Marseille, France
^{6} Astrophysics Group, Cavendish Laboratory, University of Cambridge, J J Thomson Avenue, Cambridge CB3 0HE, UK
^{7} Astrophysics & Cosmology Research Unit, School of Mathematics, Statistics & Computer Science, University of KwaZuluNatal, Westville Campus, Private Bag X54001, 4000 Durban, South Africa
^{8} CITA, University of Toronto, 60 St. George St., Toronto, ON M5S 3H8, Canada
^{9} CNRS, IRAP, 9 Av. colonel Roche, BP 44346, 31028 Toulouse Cedex 4, France
^{10} California Institute of Technology, Pasadena, California, CA 91109 USA
^{11} Centre for Theoretical Cosmology, DAMTP, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK
^{12} Computational Cosmology Center, Lawrence Berkeley National Laboratory, Berkeley, California, CA 94720, USA
^{13} DTU Space, National Space Institute, Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark
^{14} Département de Physique Théorique, Université de Genève, 24 Quai E. Ansermet, 1211 Genève 4, Switzerland
^{15} Departamento de Astrofísica, Universidad de La Laguna (ULL), 38206 La Laguna, Tenerife, Spain
^{16} Departamento de Física, Universidad de Oviedo, Avda. Calvo Sotelo s/n, 33007 Oviedo, Spain
^{17} Department of Astrophysics/IMAPP, Radboud University Nijmegen, PO Box 9010, 6500 GL Nijmegen, The Netherlands
^{18} Department of Physics & Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, British Columbia, Canada
^{19} Department of Physics and Astronomy, Dana and David Dornsife College of Letter, Arts and Sciences, University of Southern California, Los Angeles, CA 90089, USA
^{20} Department of Physics and Astronomy, University College London, London WC1E 6BT, UK
^{21} Department of Physics and Astronomy, University of Sussex, Brighton BN1 9QH, UK
^{22} Department of Physics, Gustaf Hällströmin katu 2a, University of Helsinki, 00560 Helsinki, Finland
^{23} Department of Physics, Princeton University, Princeton, New Jersey, NJ 08544, USA
^{24} Department of Physics, University of California, Berkeley, California, CA 94720, USA
^{25} Department of Physics, University of California, One Shields Avenue, Davis, California, CA 93106, USA
^{26} Department of Physics, University of California, Santa Barbara, California, CA 93106, USA
^{27} Department of Physics, University of Illinois at UrbanaChampaign, 1110 West Green Street, Urbana, Illinois, USA
^{28} Dipartimento di Fisica e Astronomia G. Galilei, Università degli Studi di Padova, via Marzolo 8, 35131 Padova, Italy
^{29} Dipartimento di Fisica e Astronomia, Alma Mater Studiorum, Università degli Studi di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy
^{30} Dipartimento di Fisica e Scienze della Terra, Università di Ferrara, via Saragat 1, 44122 Ferrara, Italy
^{31} Dipartimento di Fisica, Università La Sapienza, P. le A. Moro 2, 00133 Roma, Italy
^{32} Dipartimento di Fisica, Università degli Studi di Milano, via Celoria 16, 20133 Milano, Italy
^{33} Dipartimento di Fisica, Università di Roma Tor Vergata, via della Ricerca Scientifica 1, 00133 Roma, Italy
^{34} Dipartimento di Matematica, Università di Roma Tor Vergata, via della Ricerca Scientifica 1, 00133 Roma, Italy
^{35} Discovery Center, Niels Bohr Institute, Copenhagen University, Blegdamsvej 17, 1165 Copenhagen, Denmark
^{36} European Space Agency, ESAC, Planck Science Office, Camino bajo del Castillo, s/n, Urbanización Villafranca del Castillo, 28691 Villanueva de la Cañada, Madrid, Spain
^{37} European Space Agency, ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands
^{38} Gran Sasso Science Institute, INFN, viale F. Crispi 7, 67100 L’ Aquila, Italy
^{39} HGSFP and University of Heidelberg, Theoretical Physics Department, Philosophenweg 16, 69120 Heidelberg, Germany
^{40} Haverford College Astronomy Department, 370 Lancaster Avenue, Haverford, Pennsylvania, PA 19041, USA
^{41} Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of Helsinki, 00560 Helsinki, Finland
^{42} INAF–Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, 35131 Padova, Italy
^{43} INAF–Osservatorio Astronomico di Roma, via di Frascati 33, 00040 Monte Porzio Catone, Italy
^{44} INAF–Osservatorio Astronomico di Trieste, via G.B. Tiepolo 11, 34127 Trieste, Italy
^{45} INAF/IASF Bologna, via Gobetti 101, 40127 Bologna, Italy
^{46} INAF/IASF Milano, via E. Bassini 15, 20133 Milano, Italy
^{47} INFN–CNAF, viale Berti Pichat 6/2, 40127 Bologna, Italy
^{48} INFN, Sezione di Bologna, viale Berti Pichat 6/2, 40127 Bologna, Italy
^{49} INFN, Sezione di Ferrara, via Saragat 1, 44122 Ferrara, Italy
^{50} INFN, Sezione di Roma 1, Università di Roma Sapienza, Piazzale Aldo Moro 2, 00185 Roma, Italy
^{51} INFN, Sezione di Roma 2, Università di Roma Tor Vergata, via della Ricerca Scientifica, 1, 00185 Roma, Italy
^{52} Imperial College London, Astrophysics group, Blackett Laboratory, Prince Consort Road, London, SW7 2AZ, UK
^{53} Institut d’Astrophysique Spatiale, CNRS, Univ. ParisSud, Université ParisSaclay, Bât. 121, 91405 Orsay Cedex, France
^{54} Institut d’Astrophysique de Paris, CNRS (UMR 7095), 98bis Boulevard Arago, 75014 Paris, France
^{55} Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
^{56} Institute of Theoretical Astrophysics, University of Oslo, Blindern, 0371 Oslo, Norway
^{57} Instituto de Astrofísica de Canarias, C/Vía Láctea s/n, 38205 La Laguna, Tenerife, Spain
^{58} Instituto de Física de Cantabria (CSICUniversidad de Cantabria), Avda. de los Castros s/n, 39005 Santander, Spain
^{59} Istituto Nazionale di Fisica Nucleare, Sezione di Padova, via Marzolo 8, 35131 Padova, Italy
^{60} Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, California, CA 91109, USA
^{61} Jodrell Bank Centre for Astrophysics, Alan Turing Building, School of Physics and Astronomy, The University of Manchester, Oxford Road, Manchester, M13 9PL, UK
^{62} Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA
^{63} Kavli Institute for Cosmology Cambridge, Madingley Road, Cambridge, CB3 0HA, UK
^{64} LAL, Université ParisSud, CNRS/IN2P3, 91405 Orsay, France
^{65} LERMA, CNRS, Observatoire de Paris, 61 Avenue de l’Observatoire, 75000 Paris, France
^{66} Laboratoire Traitement et Communication de l’Information, CNRS (UMR 5141) and Télécom ParisTech, 46 rue Barrault, 75634 Paris Cedex 13, France
^{67} Laboratoire de Physique Subatomique et Cosmologie, Université GrenobleAlpes, CNRS/IN2P3, 53 rue des Martyrs, 38026 Grenoble Cedex, France
^{68} Laboratoire de Physique Théorique, Université ParisSud 11 & CNRS, Bâtiment 210, 91405 Orsay, France
^{69} Lawrence Berkeley National Laboratory, Berkeley, California, USA
^{70} MaxPlanckInstitut für Astrophysik, KarlSchwarzschildStr. 1, 85741 Garching, Germany
^{71} Mullard Space Science Laboratory, University College London, Surrey RH5 6NT, UK
^{72} Nicolaus Copernicus Astronomical Center, Bartycka 18, 00716 Warsaw, Poland
^{73} Niels Bohr Institute, Copenhagen University, Blegdamsvej 17, 1165 Copenhagen, Denmark
^{74} Nordita (Nordic Institute for Theoretical Physics), Roslagstullsbacken 23, 106 91 Stockholm, Sweden
^{75} SISSA, Astrophysics Sector, via Bonomea 265, 34136 Trieste, Italy
^{76} School of Chemistry and Physics, University ofKwaZuluNatal, Westville Campus, Private Bag X54001, 4000 Durban, South Africa
^{77} School of Physics and Astronomy, University of Nottingham, Nottingham NG7 2RD, UK
^{78} Simon Fraser University, Department of Physics, 8888 University Drive, Burnaby BC, Canada
^{79} Sorbonne UniversitéUPMC, UMR 7095, Institut d’Astrophysique de Paris, 98bis Boulevard Arago, 75014 Paris, France
^{80} Space Research Institute (IKI), Russian Academy of Sciences, Profsoyuznaya Str, 84/32, 117997 Moscow, Russia
^{81} Space Sciences Laboratory, University of California, Berkeley, California, CA 94720, USA
^{82} SubDepartment of Astrophysics, University of Oxford, Keble Road, Oxford OX1 3RH, UK
^{83} The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, AlbaNova, 106 91 Stockholm, Sweden
^{84} UPMC Univ Paris 06, UMR7095, 98bis Boulevard Arago, 75014 Paris, France
^{85} Université de Toulouse, UPSOMP, IRAP, 31028 Toulouse Cedex 4, France
^{86} University of Granada, Departamento de Física Teórica y del Cosmos, Facultad de Ciencias, 18071 Granada, Spain
^{87} Warsaw University Observatory, Aleje Ujazdowskie 4, 00478 Warszawa, Poland
^{⋆}
Corresponding author: J. L. Puget, email: jeanloup.puget@ias.upsud.fr
Received: 10 May 2016
Accepted: 14 July 2016
This paper describes the identification, modelling, and removal of previously unexplained systematic effects in the polarization data of the Planck High Frequency Instrument (HFI) on large angular scales, including new mapmaking and calibration procedures, new and more complete endtoend simulations, and a set of robust internal consistency checks on the resulting maps. These maps, at 100, 143, 217, and 353 GHz, are early versions of those that will be released in final form later in 2016. The improvements allow us to determine the cosmic reionization optical depth τ using, for the first time, the lowmultipole EE data from HFI, reducing significantly the central value and uncertainty, and hence the upper limit. Two different likelihood procedures are used to constrain τ from two estimators of the CMB E and Bmode angular power spectra at 100 and 143 GHz, after debiasing the spectra from a small remaining systematic contamination. These all give fully consistent results. A further consistency test is performed using crosscorrelations derived from the Low Frequency Instrument maps of the Planck 2015 data release and the new HFI data. For this purpose, endtoend analyses of systematic effects from the two instruments are used to demonstrate the near independence of their dominant systematic error residuals. The tightest result comes from the HFIbased τ posterior distribution using the maximum likelihood power spectrum estimator from EE data only, giving a value 0.055 ± 0.009. In a companion paper these results are discussed in the context of the bestfit PlanckΛCDM cosmological model and recent models of reionization.
Key words: cosmology: observations / dark ages, reionization, first stars / cosmic background radiation / space vehicles: instruments / instrumentation: detectors
© ESO, 2016
1. Introduction
The Emode polarization signal of the cosmic microwave background (CMB) at multipoles less than 15 is sensitive to the value of the Thomson scattering optical depth τ. In polarization at large angular scales, the extra signal generated by reionization dominates over the signal from recombination. Cosmic microwave background polarization measurements thus provide an important constraint on models of early galaxy evolution and star formation, providing the integrated optical depth of the entire history of reionization, which is complementary information to the lower limit on the redshift of full reionization provided by Lymanα absorption in the spectra of high redshift objects (see Dijkstra 2014; Mashian et al. 2016; Robertson et al. 2015; Bouwens et al. 2015; Zitrin et al. 2015, for recent results and reviews). The reionization parameter τ is difficult to constrain with CMB temperature measurements alone as the TT power spectrum depends on the combination A_{s}e^{− 2τ}, and τ is degenerate with A_{s}, the amplitude of the initial cosmological scalar perturbations.
The Bmode polarization signal at low multipoles is created by the scattering of primordial tensor anisotropies in the CMB by reionized matter. This signal scales roughly as τ^{2}. The value of τ is thus also important for experiments constraining the tensortoscalar ratio, r, using the reionization peak of the B modes.
Owing to the large angular scale of the signals, τ has been measured only in fullsky measurements made from space. Hinshaw et al. (2013) report τ = 0.089 ± 0.014 from the WMAP9 analysis. Planck Collaboration XI (2016) used Planck^{1} polarized 353GHz data to clean the WMAP Ka, Q, and V maps of polarized dust emission, and WMAP Kband as a template to remove polarized synchrotron emission, lowering τ from WMAP data alone by about 1σ to τ = 0.075 ± 0.013. Planck Collaboration XIII (2016) used 70GHz polarization data at low multipoles, which were assumed to be nearly noiselimited (see Sect. 4 below for new estimates of systematic effects) and were cleaned using 353 GHz and 30 GHz to remove dust and synchrotron emission. This lowmultipole likelihood alone gave τ = 0.067 ± 0.023. In combination with the Planck highmultipole TT likelihood, it gave τ = 0.078 ± 0.019. The combination PlanckTT+lensing gave a comparable value of τ = 0.070 ± 0.024.
The HFI 100 and 143GHz lowmultipole polarization data combined with 353GHz to remove the Galatic dust foreground could provide tighter constraints on τ, reducing uncertainties by nearly a factor of about two with respect to the results quoted above; however, systematic errors remaining in the HFI polarized data at large angular scales in the 2015 Planck release led the Planck collaboration to delay their use. Since the time when the data for the 2015 release were frozen, we have made substantial improvements in the characterization and removal of HFI systematic errors, which we describe in this paper. The maps, simulations, and sometimes computer codes used here are referred to as “pre2016”, and are nearly identical to those in the 2016 data release to come. We note that the HFI submillimetre channels at 545 and 857 GHz are not polarization sensitive and their improvement is not discussed in this paper.
The paper is structured in the following way. Section 2 analyses systematic errors affecting the HFI timeordered information processing, describes tests of these errors, and describes new versions of polarized mapmaking and calibration. Section 3 discusses global tests of HFI polarization at the power spectrum level. Section 4 discusses quality tests of Low Frequency Instrument (LFI) polarization data. Section 5 describes component separation and crossspectra and their suitability for the measurement of τ. Section 6 describes the likelihood analysis of the τ parameter. Finally Sect. 7 discusses the implications of our new τ value.
These data described here are also used in an accompanying paper (Planck Collaboration Int. XLVII 2016) to constrain different models for the reionization history of the Universe.
2. Improvements of HFI mapmaking and calibration for low multipoles
2.1. Overview
Analysis of HFI polarization systematic effects is complex and technical in nature. This paper, together with Planck Collaboration VII (2016), completes the characterization of the HFI timeordered information (TOI) products, extending analysis of the polarization power spectra to the lowest multipoles.
The HFI TOI used in this paper are described in the 2015 data release (Planck Collaboration VII 2016), and a calibrated version is publicly available in the Planck Legacy Archive^{2}. The spacecraft pointing solution (Planck Collaboration I 2016) and all instrumental parameters, including the focal plane geometry, crosspolar leakage parameters, and polarization angle parameters, are identical to those used in the 2015 release.
HFI systematic errors have been discussed extensively in earlier papers (Planck Collaboration X 2014; Planck Collaboration VII 2016) and were shown to be under control on small angular scales. Large angular scales, however, were still affected by notfullyunderstood systematic effects.
The reduction of systematic errors in HFI polarization measurements at large angular scales reported in this paper is the result of applying a new mapmaking code, SRoll. SRoll is a polarized destriping algorithm, similar to the 2015 HFI mapmaking pipeline (Planck Collaboration VIII 2016), but it solves simultaneously for absolute calibration, leakage of temperature into polarization, and various other sources of systematic error. Systematic residuals are removed via template fitting, using HEALPixbinned rings (HPR) to compress the timeordered information. Each HPR consists of data from a given stable pointing period (ring) binned into HEALPix pixels (Górski et al. 2005). A version of this code with minor additions will be used for the final Planck data release in 2016. In its present form, SRoll reduces the systematic errors in the HFI EE polarization spectra at multipoles ℓ ≥ 2 to levels close to the nominal instrumental noise. This allows a measurement of the reionization optical depth τ from the HFI alone, which was not possible at the time of the 2015 Planck data release.
The rest of this section discusses each of the specific improvements. Section 2.2 introduces the endtoend simulations used in this paper. Section 2.3 analyses the detector system noise, which sets a fundamental limit to the measurement of the polarization spectra and to τ with HFI data. (Appendix A summarizes a number of effects removed in the TOI processing that make a negligible contribution to the systematic error budget of polarization measurements at large angular scales.) Section 2.4 describes removal of the zodiacal light, and far sidelobes. Section 2.5 characterizes the systematic errors induced by nonlinearities in the analoguetodigital converters (ADCs) used in the bolometer readout electronics. These nonlinearities, if uncorrected, would be the dominant source of systematic errors in polarization at large angular scales. Section 2.6 discusses all other systematic errors not corrected after the TOI processing, and mitigated by SRoll.
2.2. Endtoend simulations
Endtoend (E2E) HFI simulations, introduced for the first time in 2015 (Planck Collaboration XI 2016), include an accurate representation of the response of the instrument to sky signals, as well as known instrumental systematic errors at the TOI or ring level, all of which are propagated through the entire analysis pipeline to determine their effect on the final maps. The version used in this paper runs the pre2016 Sroll mapmaking code, described in detail in Sect. B.3. Simulations built with this pre2016 pipeline are referred to as HFI focal plane simulations (HFPSs). Three sets of HFPSs are used in this paper: HFPS1 contains 83 realizations; HFPS2 and HFPS3 contain 100 realizations each, which can be used either separately or together.
2.3. Detector noise
Detector noise is determined through a multistep process, starting with firstorder correction of the TOI for ADC nonlinearity, demodulation, deglitching, 4K line removal, and timeresponse deconvolution. The sky signal is then estimated (using the redundancy of multiple scans of the same sky during each stable pointing period) and removed, leaving the noise, which can then be characterized (see Planck Collaboration VII 2016, for details).
Fig. 1 Mean power spectra of the signalsubtracted, timeordered data from Survey 2 for each polarizationsensitive HFI frequency channel. The spectra are normalized at 0.25 Hz. Blue, green, red, and cyan represent 100, 143, 217, and 353 GHz, respectively. The vertical dashed line marks the spacecraft spin frequency. The sharp spikes at high frequencies are the socalled 4K cooler lines. These noise spectra are built before the time transfer function deconvolution. 
Within a frequency band, the noise power spectrum of each bolometer has a similar shape below 5 Hz. Figure 1 shows the mean noise spectrum for each of the polarizationsensitive HFI channels, normalized to unity at 0.25 Hz (we are interested here in the shapes of the noise spectra, not their absolute amplitude).
Fig. 2 Noise crosspower spectra of the 143GHz bolometers, with the unpolarized spiderweb bolometers (SWBs) in red and the polarizationsensitive bolometers (PSBs) in blue. The lowlevel correlated white noise component of the PSB noise is associated with common glitches below the detection threshold. Autospectra are shown in black. The uncorrelated noise is in green. 
Detector noise can be divided into three components (see Fig. 2), with spectra varying approximately as f^{0} (i.e., white), f^{1}, and f^{2}. The f^{0} and f^{1} components are uncorrelated between detectors; however, the f^{2} component is correlated between detectors. This component dominates below 10^{2} Hz, and thus below the spin frequency. Noise levels at frequencies around 0.002 Hz (much below the spin frequency) changed during the mission, reflecting variations of the Galactic cosmic ray flux with solar modulation and their effect on the temperature of the bolometer plate. These effects are corrected using a template built from the signal of the dark bolometers smoothed on minute timescales (Planck Collaboration VI 2014). The f^{2} component is different for spiderweb bolometers (SWBs, red curves) and polarizationsensitive bolometers (PSBs, blue curves). Common glitches seen in the two silicon wafers of a PSB produce correlated noise with a knee frequency of about 0.1 Hz (Planck Collaboration X 2014). Above 0.2 Hz, the correlated noise is at a level of 1% of the total detector noise; at twice the spin frequency (multipole ℓ = 2 in the maps), it contributes only 10% of the total noise.
After removal of the correlated part, an uncorrelated f^{1} component remains (green curve in Fig. 2) that has f_{knee} ≈ 0.2 Hz for all HFI detectors (see, e.g., Fig. 13 of Planck Collaboration X 2014). Above this frequency the noise is predominantly white, with amplitude in good agreement with groundbased measurements. The 10MΩ resistor and the capacitor in the focal plane, read out through the same electronics as the bolometers, show only 1 /f noise below 3 mHz, showing that this additional 1 /f component is not generated in the readout. The 1 /f component of the noise is discussed in Sect. A.1, and is well approximated by a Gaussian statistical distribution, as seen in Fig. 3. The low level (<10^{3}) nonGaussian wings which are cut by the 3.3σ clipping of the deglitching algorithm, are not seen on 100ring averages, which suggest that they are caused by rare events.
Fig. 3 Histogram of the noise between 0.018 and 0.062 Hz (frequencies at which it is dominated by the uncorrelated 1 /f noise) for detector 1431a in blue, together with the bestfit Gaussian distribution in red. 
1 /f noise around 0.06 Hz is largely uncorrelated, and is constant during the mission within uncertainties. This confirms that undetected glitches do not contribute significantly, apart from a small contribution near the spin frequency. No intrinsic 1 /f noise was detected in individual bolometer ground tests, but a 1 /f component was seen in HFI focal plane unit ground calibration measurements below 0.1 Hz (see, e.g., Fig. 10 of Lamarre et al. 2010), comparable to those seen in flight. Nevertheless a thermal fluctuation residual cannot be excluded.
In summary, the correlated component (f^{2}) between all bolometers is mostly removed by the dark bolometer baseline removal. The weak correlated white noise (1%) between pairs of bolometers within a PSB, caused by undetected glitches, is not removed by this baseline removal procedure. The oneminute smoothing of the dark bolometer timeline used to remove this common thermal mode will leave a small correlated residual in the frequencies just above the spin frequency. Furthermore the ADC nonlinearities discussed later will affect differently each bolometer and introduce some uncorrelated residuals which might account for at least part of this 1 /f component.
In conclusion, the component of TOI noise that is uncorrelated between detectors (except near the spin frequency) can be modelled with two Gaussian components, white and 1 /f. These are taken from version 8 of the FFP8 noise model Planck Collaboration XII (2016), adjusted to the power spectrum observed in the halfring null tests. The TOI noise, and its propagation to maps and power spectra, is taken in this paper as the fundamental limitation of HFI, and the maps built with FFP8 noise give the reference for maps and power spectra in simulations.
2.4. Sky components removed in the ring making
2.4.1. Removal of zodiacal dust emission
Thermal emission from interplanetary dust – the zodiacal emission – varies not only with frequency and direction, but also with time as Planck moves through the solar system. The size of the effect on temperature is small, as estimated in Planck Collaboration VIII (2016), and it was not removed from the TOI released in 2015. Rather, we reconstruct the emission using the COBE Zodiacal Model (Kelsall et al. 1998) and the zodiacal emissivity parameters found in Planck Collaboration VIII (2016), and subtract it from the TOI of each detector prior to mapmaking.
The zodiacal emission is not expected to be intrinsically polarized; however, HFI measures polarization by differencing the signals from quadruplets of PSBs. Since each detector has a slightly different spectral response, any component that has a different spectrum from the blackbody spectrum of the primordial CMB anisotropies, such as zodiacal emission, can introduce leakage of temperature into polarization. By subtracting a model of zodiacal emission from the TOI of each detector, this source of leakage is strongly suppressed.
Differences between Q and U maps made with and without removal of zodiacal emission show the size of the effect. At 100, 143, and 217 GHz, the zodiacal correction is less than about 100 nK_{CMB}. At 353 GHz, it is an order of magnitude larger. If our correction was so poor that it left 25% of the emission, temperature leakage would contribute errors of order 25 nK_{CMB} to the polarization maps at 100, 143, and 217 GHz, and perhaps 100 nK_{CMB} at 353 GHz. The effect on τ is negligible.
Estimates of the relative impact of the FSLs per bolometer within a frequency band.
2.4.2. Far sidelobes
HFI far sidelobes (FSLs) are discussed in detail in Planck Collaboration VII (2016). They are dominated by radiation from the feedhorns spilling over the edge of the secondary mirror, and by radiation reflected by the secondary mirror spilling over the edge of primary mirrorr or the main baffle^{3}. The effects of diffraction by edges other than the mouths of the feedhorns themselves is important only at 100 GHz, where diffraction can affect very low multipoles (see below). FSL variations inside a frequency band are negligible at HFI frequencies. The Planck optical system is modelled using the GRASP software^{4}. Higherorder effects with more than one reflection and diffraction by sharp edges can be computed by GRASP, but the complexity and computational time is a strong function of the order number. We have computed the firstorder FSLs for each bolometer, and then checked for a few representative bolometers that the addition of the next seven orders gives small corrections to the firstorder computations.
Convolution of the sky maps with the FSL model predicts small contributions to the maps. In addition, the contributions from FSLs close to the spin axis are reduced by destriping during the mapmaking process (a contribution of FSLs on the spinaxis is completely removed by the destriper). The Galactic contributions though the FSLs are similar at 100, 143, and 217 GHz since the decrease in the FSL amplitude and the increase in Galactic emission at higher frequencies roughly compensate each other. The solar dipole FSL contributions decrease with frequency, as summarized in Table 1. The FSLs produce a direct dipole calibration shift by increasing the effective beam etendue, though the shift is reduced by the destriping. The dispersion between detectors is caused by the variation in the main spillover (the rays that miss the secondary mirror), which is dominated by the position of the horns in the focal plane: the bolometers further away from the symmetry axis experience a larger spillover. This effect can be used as an additional test of the fidelity of the GRASP calculations used to model the FSLs.
The impact of FSLs on the HFI maps depends on the scanning strategy, and so must be corrected in either the TOI processing or in building the HPRs. A FSL model has been computed from the firstorder GRASP calculations. Higherorder effects are absorbed, together with other residuals, into the empirical complex transfer functions discussed in Sect. 2.6.2. The parameters of the transfer functions are then determined from the data.
Fig. 4 Autopower spectra, showing the level of the simulated FSL projected on the maps predicted using the GRASP model. At ℓ < 10, the FSL signal at all frequencies is at least two orders of magnitude smaller than the cosmological FEE signal. 
Angular power spectra^{5} from HFI pre2016 E2E simulations of the FSLs are displayed in Fig. 4. At all HFI frequencies, the FSL effects are smaller than the fiducial FEE power spectrum by two orders of magnitude, where we define Fxx (FTT, FEE, FBB) to be the base ΛCDM power spectra from bestfit Planck 2015 cosmological parameters (Planck Collaboration XIII 2016) τ = 0.066 and upper limit r = 0.11.
The FSL signals are removed using GRASP firstorder predictions, leaving residuals due to higher orders and uncertainties. These residuals are smaller than the effect displayed in Fig. 4.
Table 1 provides estimates of the effects of FSLs on the relative intercalibration with respect to the average of all detectors within the same frequency band. As discussed in Sect. 2.4.2, the effect of the FSLs depends on the scanning strategy, and must be removed at the HPR level. The effect on the dipole requires propagation of the FSLs through the HFI pre2016 E2E simulations to take into account the filtering by the destriper. These numbers can be compared directly with the main beam dipole amplitudes measured from each bolometer. At 100 and 143 GHz, the rms dispersions of the relative dipole calibration measured from individual bolometers are 5 × 10^{6} and 9 × 10^{6}, respectively (see Fig. 13). Comparing with the rms variation of the FSL contribution to the dipole calibration listed in Table 1 (1 × 10^{4} and 2 × 10^{5}) shows that corrections using the GRASP model are accurate to better than 5 and 2%, respectively, for 100 and 143 GHz. This is not surprising since the first order (spillover) dominates, as discussed earlier in this section.
The asymmetries of the FSLs are mainly caused by the secondary mirror spillover, which depends on the position of the detector in its line parallel to the scan direction. Smaller asymmetries are generated by the small higherorder effects in the FSL GRASP calculations discussed above. Uncertainties in the very long time constants can also leave small transfer function residuals. These can be tested on the data redundancies, and when detected, corrected by an empirical complex (phase and amplitude) transfer function (see Sect. 2.6.2) at low frequencies.
2.5. ADC nonlinearity systematic effects
Nonlinearity in the ADCs^{6} in the bolometer readouts (Planck Collaboration VI 2014; Planck Collaboration VII 2016) introduces systematic errors in the data.
Fig. 5 Relationship between input and output, for one spare ADC. The plot shows the difference between the measured digitized output signal level and the one with a perfectly linear ADC, as a function of the output level, over a signal range appropriate for the sky signals. Thus, on the one hand, a perfectly linear device would be a horizontal line at zero; in a real device such as shown here, on the other hand, the relationship between input and output is complicated and nonlinear everywhere, especially near the middle of the range around 0 ADU. 
Figure 5, for example, shows the deviations from linearity measured in a flight spare ADC. A perfectly linear device, in contrast, would lead to a horizontal line at zero. The strongest nonlinearity in Fig. 5 at zero ADU, lies in the middle of the ADC range, and drives the main ADC nonlinearity effect. ADC deviations from linearity are the dominant source of systematic errors in HFI lowℓ polarization data. They create a firstorder variable gain for the readout electronics of each detector, as well as higherorder effects as described below.
Fig. 6 Fastsampled signal from bolometer 1431a. Left: eighty samples from early in the mission, corresponding to one cycle of the 90Hz squarewave modulation of the bias current across the bolometer. Because of the squarewave modulation, the first and last sets of 40 samples are nearly mirror images of one another across the xaxis. The units are raw analoguetodigital units (ADU). Middle: normalized histogram of the fastsample signal values for day 91 of the mission, in ADU. The signal spread is dominated by the combination of the noise and the squarewave modulation of the bias current across the bolometer (left panel), with additional contributions from the CMB solar dipole. Right: normalized histograms for each day, starting with day 91, stacked left to right for the entire mission. Histogram values are given by colour, as indicated. The obvious symmetric trends during the mission are caused by drifts in temperature of the bolometer plate. Isolated days with large deviations from zero, seen as narrow black vertical lines, are due to solar flares. 
In the 2013 data release, the effects of ADC nonlinearities were partially corrected (during mapmaking) by application of a timevariable gain for each detector (i.e., calibration to astrophysics units; Planck Collaboration VI 2014), with residual effects at levels of a few percent or less in the TOI. This simplified approach to ADC nonlinearity was adequate for the analysis of temperature anisotropies; however, it leaves residual effects in the polarization maps at the lowest multipoles at levels about 30 times higher than the noise.
In 2013, a proper correction that would take into account the detailed characteristics of the ADCs themselves at the TOI processing level was not possible because the prelaunch measurements of the ADCs had inadequate precision (Planck Collaboration VII 2016). Accordingly, during the extended Planck observations with LFI, after the ^{3}He for the 0.1K cooler ran out (February 2012 to August 2013), we conducted an inflight measurement campaign to characterize the HFI ADC nonlinearities more accurately. The resulting model of ADC nonlinearity was used to correct the 2015 HFI data (in TOI processing). Residual effects were much reduced, to a level of 0.2−0.3% in the TOI, and had negligible effect on the 2015 polarization maps at multipoles ℓ ≳ 30. No additional correction of the type made in 2013 was performed because it did not bring clear improvement. (This was subsequently understood to be a consequence of degeneracies with other corrections.) However, a better treatment of ADC nonlinearities is required to achieve noiselimited polarization maps at multipoles less than 30.
The results in this paper are based on an improved treatment of ADC nonlinearities, made possible by the new mapmaking code (Sroll), which simultaneously solves for temperaturetopolarization leakage and residual gain variations from ADC nonlinearities. Specifically, the 2015 correction in the TOI (using the model of ADC nonlinearity derived from the extendedmission measurements) is followed by a step in Sroll that calculates a multiplicative gain correction of the residuals left by the first step.
In this section, we explain why application of a timevariable gain worked reasonably well for temperature in 2013, why a model of ADC nonlinearity worked well enough for highℓ polarization in 2015, and why a global determination of leakage levels, together with the residual gain variations induced by the ADC nonlinearity, works better yet, and can be used in this paper for lowℓ polarization. In addition, we show that a higherorder (but nonnegligible) of ADC nonlinearity acting on the CMB dipole also has to be taken into account.
To begin, consider how the signal levels at the inputs of the ADCs change with time. Figure 6 shows the fastsampled signal from a single bolometer (1431a) for three different time periods. The left panel shows the 80 fast samples in one cycle of the 90Hz squarewave modulation of the bias current across the bolometer. Planck Collaboration VII (2016) and Fig. 13 of Lamarre et al. (2010) present details and an explanation of the shape, which varies from bolometer to bolometer and also changes slowly throughout the mission. A squarewave compensation voltage is subtracted from the signal at the input of the readout electronics to bring both modulations close to zero, in order to limit nonlinearity effects in the analogue amplification stages. The middle panel of Fig. 6 shows a normalized histogram of all the fastsample signal values for the first day of observations, day 91^{7}. The signal spread is dominated by the noise, combined with the squarewave modulation of the bias current across the bolometer (left panel), with additional contributions from the CMB solar dipole. The right panel shows daily normalized histograms of the fast samples of detector 1431a for each day in the mission, starting with day 91, stacked left to right, with the histogram values colourcoded, as indicated. The two modulation states of the signal and their evolution on the ADC are clearly seen as positive and negative bands. The large, longterm, symmetric trends during the mission are caused mainly by the slow temperature drift in the bolometer plate. Additional asymmetric drifts are due to longterm variations in the readout electronics.
Because the spread of the signal is broader than the sky signal by an order of magnitude (comparing Figs. 6 and 7), it combines the various discontinuities shown in Fig. 5 into a small but complex relationship between the signal and the power on the detectors. Nevertheless, approximations are possible. In 2013, as mentioned at the beginning of this section, this effective gain was calculated for every pointing period. The TOI were then corrected with these gains smoothed by a boxcar average over 50 pointing periods. This corrected the main gain effects of the ADC nonlinearity as the signal level drifts slowly throughout the mission. In addition to the gain variation the higherorder nonlineariries distort significantly the shape of the dipole.
The shape of the input signal due to the 90Hz modulation discussed above, and an estimate of how it changed throughout the HFI lifetime, were established for each bolometer (see Figs. 2−4 of Planck Collaboration VII 2016). This permitted an approximate empirical reconstruction of the input signal at the level of the TOI.
The effects of ADC nonlinearity on the dipole signal amplitude give an excellent measure of the gain, which can be applied linearly to all signals. We note that the distortion of the shape of the largescale CMB anisotropies is negligible, but not the dipole distortion, which leaves nonnegligible, additive, largeangularscale residuals. Figure 7 shows a simulation that contains only the dipole signal and the thermal baseline drift throughout the mission, for four representative bolometers.
Fig. 7 Simulation of the modulated noiseless sky signal per ring shows the variable amplitude of the dipole and thermal drift, which dominates the signal, as a function of time (expressed in stable pointing period number) during the entire mission for four representative bolometers. The different signs of the drifts derivative depends on the level of the compensation of the modulation, which can bring one state of the modulation on either side of the middle of the ADC. 
The amplitude of the dipole signal in a given ring changes with the offset between the spin axis of Planck and the axis of the solar dipole. When the two are nearly aligned, the amplitude on a ring is small. When the two are far apart, the amplitude can reach 5−30 ADU units or so (noticeably less at 353 GHz).
In 2015, the model of ADC nonlinearities developed during the warm mission was applied to the TOI data (Planck Collaboration VII 2016). The correction takes into account the shape of the bolometer modulation (left panel of Fig. 6), and corrects for the ADC nonlinearity induced by longterm drifts, resulting in residual effects in the maps an order of magnitude smaller than the 2013 correction, and good enough to be usable for highℓ polarization. However, the detailed shape of the dipole signal after passage through the bolometer readout circuits, and at the input of the ADC, is still not taken into account, leaving residuals in the data that are too large for lowℓ polarization .
Figure 19 of Planck Collaboration VII (2016) shows the angular power spectra of difference maps, made with and without the ADC nonlinearity correction, for individual bolometers at 100, 143, 217, and 353 GHz. Before the ADC nonlinearity correction, errors caused by nonlinearities scale approximately as ℓ^{2} below multipoles ℓ ≲ 100. After correction, the power spectra are reduced in amplitude by factors of between 10 and 100 at low multipoles.
Clearly the ADC nonlinearity correction performed in the 2015 analysis is a big improvement over the variablegain adjustment used in the 2013 results. Nevertheless, the residual errors from the nonlinearity correction alone are still too large for accurate polarization measurements at low multipoles (Planck Collaboration VIII 2016). However, correction of these residual errors with a residual gain adjustment is now possible in the Sroll mapmaking, which solves simultaneously for temperaturetopolarization leakage levels and gain variations, something that was not possible in the previous mapmaking algorithm. The implementation of this approach is described in the following subsection.
2.6. Correction of temperaturetopolarization leakage
Once the HPRs are built with far sidelobes and zodiacal emission removed, 6% of the HPRs are discarded using the criteria described in Planck Collaboration VII (2016). Using a generalized destriper, which takes advantage of the redundancies in the scanning strategy, we solve selfconsistently for all temperaturetopolarization leakage terms:

ADC nonlinearityinduced gain variation;

additional empirical complex transfer functions;

calibration factors;

bandpass mismatch coefficients associated with foregrounds.
The destriper baseline solution avoids regions of the sky with strong gradients that could bias the baselines. We construct one sky mask per frequency using a threshold in temperature (Sect. B.2). The foreground templates for the bandpass mismatch leakage are the Commander dust, CO, and freefree maps (Planck Collaboration IX 2016). The detailed description of the method, performance from simulations, and tests are to be found in Appendix B.
2.6.1. ADCinduced gain correction
A less biased measure of gain changes is given in Fig. 8, which shows differences in dipole amplitudes measured on the same ring one year apart.
Fig. 8 Dipole amplitude difference on the rings observed twice, one year apart, for each of the 143GHz polarizationsensitive bolometers. This detects the timedependent response associated with the excursions of the signal on the ADC, illustrated in Fig. 7. The blue curve shows the dipole differences in units of μK after ADC correction in the TOI processing, with no other processing. The red curve shows the measured dipole amplitude solved by Sroll, demonstrating the reliability of the model, which can then be applied to small signals. 
Application of the ADC nonlinearity model to the TOI in the 2015 data release reduces the rms dipole amplitude variations from around 1% (no correction) to 2 × 10^{3} (blue curve). The red curve shows how well SRoll is able to solve for these residual gain adjustments, which are then applied. Algorithmic details and accuracies are described in Appendix B.
Distortions of the dipole shape caused by ADC nonlinearities within each modulation state, however, are still not accounted for. Because the drift of the signal on the ADC is large between the two halves of the mission, halfmission null tests both in the data and in simulations give an excellent test of the quality of the ADC nonlinearity corrections and of this residual dipole distortion. These are shown for 100 GHz in Fig. B.14 (see Sect. B.4.2 for details) and their associated EE power spectra in Fig. 9.
Fig. 9 Power spectra of the 100GHz halfmission nulltest maps shown in Fig. B.14, which are dominated by ADC nonlinearity effects. We compare six simulations at 100 GHz drawn from the Sroll uncertainties (red lines), and the average for the full frequency 100GHz data (blue line). In the top panel only the TOI nonlinearity correction is applied. In the bottom panel, both the TOI correction and the timevarying gain correction are applied. This leaves only the dipole distortion. For both panels the green lines are for simulations without the ADC nonlinearity effect (only noise and other systematic residuals). 
These spectra are calculated for three cases: (i) simulated data from six realizations corrected with the ADC TOI correction alone (top panel, red lines); (ii) simulated data corrected with the additional correction of Sroll residual gain variation (bottom panel, red lines); and (iii) simulated data without any ADC nonlinearity (i.e., the “ideal” case; both panels, green lines). In both panels, blue lines correspond to the data themselves processed in the same way.
The difference between the red lines in the top and bottom panels of Fig. 9 shows the efficiency of the Sroll residual gain variation correction. The difference between the red and green lines in the bottom panel shows the level of the residuals after this correction, which is dominated by dipole distortion. At ℓ ≥ 4, the red and green lines are at the same level within the noise plus small contributions from other systematic residuals. They all stay below 1 × 10^{3}μ. The only significant deviations are for the quadrupole and octopole terms. Although this is the largest systematic effect left uncorrected in the HFI polarization maps, it is rather weak (≤ 20%) compared to the expected EE signal. This effect has been wellsimulated, and the simulations will be used to remove it from the EE power spectra in the science analysis.
2.6.2. Empirical complex transfer function
In the TOI processing, the time response of the bolometers and associated readout electronics is modelled as a Fourier filter, which is determined from observations of Saturn, Jupiter, and stacked glitches. The data are deconvolved from this response function (Planck Collaboration VII 2016). In addition, an inscan, phaseshifted dipole caused by very long time constants (VLTC) is subtracted from the timeline, a necessary step for the convergence of the orbital dipole calibration, as described in Planck Collaboration VIII (2016). There are unavoidable uncertainties in the transfer function used to correct the timelines. We fit an additional empirical complex transfer function in the mapmaking, taking advantage of the redundancies in the data to capture any residuals from the abovementioned bolometer/electronics time response deconvolution. The fit also captures residuals from FSLs, which shift dipoles and low frequency signals (inscan and crossscan).
The empirical transfer function is composed of four complex quantities associated with four bands of spin frequency harmonics ([0.017 Hz], [0.033−0.050 Hz], [0.067−0.167 Hz], and [0.133−0.250 Hz]), which are adjusted to minimize the residuals in the global mapmaking (see Appendix B.4.3). Including higher spin frequency harmonics provides only negligible improvements to the maps. In the four spinharmonic frequency bands, the phase shifts of the transfer functions are easily fitted because they are not degenerate with the sky signal. However, for 100−217 GHz, the CMB+foreground signals are too weak to fix the amplitudes of the transfer functions accurately and no amplitude correction is applied.
Fig. 10 Bestfit solutions for the real and imaginary parts of the empirical additional transfer function as a function of frequency, for the 353GHz bolometers. 
Figure 10 shows this transfer function for the eight PSBs at 353 GHz, where the strong Galactic dust signal allows an accurate determination. The real part of the function measures the asymmetry between the crossscan and the inscan residuals. The imaginary part measures the shift along the scan. At 353 GHz, the imaginary part is almost negligible.
To estimate the accuracy of the empirical transfer function, we use oddminuseven Survey map differences, which are sensitive to phase shifts at low harmonics of the spin frequency. We compute a pattern map associated with a phase shift of signal. The correlation of the data with this pattern gives the residual error left in the signal after correction with the empirical transfer function. These relative errors on the signal are shown in Fig. 11. Comparing the residual errors at the four lowest sets of multipoles to those at higher multipoles, we see only upper limits below 10^{3} at 100 and 143 GHz, and three times lower at 217 and 353 GHz. This clearly demonstrates that the additional complex transfer function works well to correct phase shift residuals at the low multipoles that have been fitted. The mapmaking does not include corrections for temporal frequencies higher than 0.250 Hz, corresponding roughly to ℓ > 15, and the odd−even Survey difference test still detects some shifts in the data.
Fig. 11 Ratio of the fitted data to simulated patterns detecting the residual imaginary part of the empirical transfer function, measured in odd minus even Survey difference maps averaged for sets of harmonics. The transfer function correction has been applied only over the four first sets of harmonic ranges (ℓ < 15); higher harmonics have not been corrected by the empirical transfer function. 
Transfer function residuals also induce leakage of the Solar dipole into the orbital dipole. This leakage affects calibration differently in odd and even surveys. The solar dipole residual amplitudes per detector with respect to the average per frequency are displayed in Fig. 12.
Fig. 12 Residual solar dipole amplitude for each bolometer, by Survey. The average dipole at each frequency is subtracted. For 100 and 143 GHz (top panels), the variations are compatible with the relative calibration uncertainty of 10^{4}. At 353 GHz, the scale is expanded by a factor of five, and all detectors show an obvious odd/even pattern, which is marginally apparent at 217 GHz. 
The residual amplitude provides a strong test of the improvement provided by the transfer function correction in reducing the leakage between dipoles and gain differences between odd and even Surveys. At CMB frequencies (100 and 143 GHz), this figure does not show any systematic odd/even Survey behaviour at the level of 0.2 μK. This translates into an upper limit on dipoleleakageinduced miscalibration better than 0.01% for each bolometer. Nevertheless the odd−even differences of dipole amplitudes at 353 GHz are apparent for all bolometers, with an amplitude up to ± 1.5 μK or approximately 0.1% in odd−even miscalibration. The empirical real part of the transfer function cannot be determined for the dipole. We note that Survey 5 is affected by residuals from solar flares and endoflife tests; the last part of Survey 5 will be removed from the 2016 data release.
We propagate the residual uncertainties of the empirical transfer function (as determined from the fitting procedure; Appendix B.4.3) to the maps and power spectra using the E2E simulations. The autopower spectra (Fig. B.16) show that over the reionization peak the residuals are at a very low level: D(ℓ) < 10^{4}μK^{2} for CMB channels^{8}; and D(ℓ) ≈ 10^{3}μK^{2} for 353 GHz, except at ℓ = 2, where it reaches 10^{2}μK^{2}. This demonstrates that the residuals of these systematic effects have a negligible effect on τ measurements, except possibly at ℓ = 2.
2.6.3. Interdetector calibration of the polarizationsensitive bolometers
Calibration mismatches between detectors produce leakage of temperature to polarization. The interdetector relative calibration of the PSBs within a frequency band can be tested on singledetector, temperatureonly maps (the polarized signal of the fullfrequency map is subtracted from the detector TOI). The bestfit solar dipole (determined from the 100 and 143GHz fullfrequency maps) is removed. The residual dipoles in the maps measure the relative calibration of each detector with respect to the average over the frequency band. The results are shown in Fig. 13.
Fig. 13 Relative calibration measured by the dipole amplitude for each polarizationsensitive detector with respect to the mean dipole across a frequency band. 
The low level of variations constrain the residual calibration mismatch between detectors that could lead to leakage of temperature to polarization. The relative calibration factors averaged over the full mission show rms dispersions of ± 5 × 10^{6} for 100 GHz, 9 × 10^{6} for 143 GHz, 2 × 10^{5} for 217 GHz, and 2 × 10^{4} for 353 GHz. At 100 GHz, 143 GHz, and 217 GHz, these dispersions are very small, and unprecedented for a CMB experiment. Absolute calibrations of band averages are discussed in Sect. 5.1. At CMB frequencies, results are still affected by gain errors between bands and by residual gain variations over time. At these frequencies, the gain mismatch is consistent with statistical errors (Tristram et al. 2011), therefore, it is not possible to improve the gain mismatch any further at these frequencies. At 353 GHz, the gain mismatch is larger than the statistical errors. The worst outliers (bolometers 3536a and 3536b) also show large odd−even discrepancies in Fig. 11. Therefore there is hope for improving the relative calibration.
2.6.4. Tests of bandpass mismatch leakage coefficients on gain: comparisons with ground tests
The SRoll mapmaking procedure (Appendix B) solves for temperaturetopolarization leakage resulting from the different response that each bolometer has to a foreground with an SED different from that of the CMB anisotropies. The solved bandpass mismatch coefficient associated with thermal dust emission is compared in Fig. 14 to that expected from prelaunch measurements of the detector bandpass (Planck Collaboration IX 2014). The statistical uncertainties on the ground measurements are dominated by systematic effects in the measurements. For all bolometers except 1433b, the two are consistent to within the error bars estimated using simulations. The smaller error bars are the sky determinations by SRoll and are in closer agreement with ground measurements than the conservative estimates of the systematic errors on the ground measurements would predict. The only exception is bolometer 3b, for which the ground and sky measurements differ by nearly 3 × the more accurate sky uncertainty.
Fig. 14 Comparison of the response of each detector to dust as measured on the ground (blue) and as solved by SRoll (red). 
2.6.5. Summary of improvements
The generalized destriper solution, solving simultaneously for bandpassmismatch leakage, intercalibration errors, and ADCinduced gain variations and dipole distortions, has been shown to be necessary to achieve a nearly complete correction of the ADC nonlinearities. This leads to much improved maps at low multipoles compared to previous releases. In Sect. 6.2 and Appendix B we will demonstrate that SRoll mapmaking does not filter or affect the CMB signal itself.
At 100, 143, and 217 GHz, we are now close to being noiselimited on all angular scales, with small remaining systematic errors due to the empirical ADC corrections at the mapmaking level. These separate corrections should be integrated in the TOI/HPR processing for better correction. At 353 GHz, and to a lesser extent 217 GHz, we observe residual systematic calibration effects, as seen in Figs. 12 and 13, but we will show in the following sections that the residuals are small enough to have negligible effect on the determination of τ. The origin of this effect is not fully understood; correction algorithms are in development.
3. Consistency tests of the HFI polarization maps
As described above and in Appendix B, the SRoll mapmaking algorithm corrects simultaneously for several sources of temperaturetopolarization leakage that were not previously corrected.
Section 3.2 gives the results of null tests that show how systematic effects at large angular scales are very significantly reduced compared to those in the HFI polarization maps from the 2015 Planck data release. As a test of the accuracy of this process, the results are also compared with the HFI pre2016 E2E simulations. We begin in the next subsection (Sect. 3.1) by showing that detection of the a posteriori crosscorrelations of the final maps with leakage templates cannot work because of the degeneracy with the dipole distortion.
3.1. Temperaturetopolarization leakage
As discussed in Sect. 2, any bandpass and calibration mismatch between bolometers induces temperaturetopolarization leakage, and hence spurious polarization signals. Each leakage pattern on the sky for each bolometer is fully determined by the scanning strategy, along with a set of leakage coefficients, and the temperature maps involved (dipole and foregrounds). The SRoll approach improves on the 2013 and 2015 HFI mapmaking pipeline by correcting all temperaturetopolarization mismatches to levels where they are negligible. Detector intercalibration has been much improved, as shown in Sect. 2.6.3. Similarly the residual bandpass leakage (mainly due to dust and CO) in the Q and U HFI pre2016 maps is also greatly reduced (Appendix B.4.1).
In the 2015 data release, we used leakage template fitting (Planck Collaboration VIII 2016) to check a posteriori the level of temperaturetopolarization leakage residuals, although this leakage was not removed from the maps. From Sect. 2.6 and the appendices, we expect temperaturetopolarization leakage to be very much reduced for the HFI pre2016 data set. Contrary to expectations, however, the templatefitting test (expanded to account for synchrotron polarized emission) on the HFI pre2016 release used in this paper still reveals significant leakage. Suspecting that the problem is residual dipole distortion induced by ADC nonlinearity, which is not yet removed, we performed leakage tests on simulations both with and without the effect. This is done by forcing the gain to be constant within a pointing period, which removes only the dipole distortion within the ADC nonlinearity effects.
Figure 15 shows that autospectra exhibit a leakage in the pre2016 100GHz data (green line), comparable to the simulations with (blue line) full ADC nonlinearity effects. The red line, with the ADC dipole distortion effect removed, is lower by a factor of 5 or more.
Fig. 15 Simulation of the templatefitting tests for temperaturetopolarization leakage in the 100GHz maps, with (blue) and without (red) ADCnonlinearity dipole distortion. The green curve shows the result of the leakage fit in the HFI pre2016 data, which is comparable to the simulation with the ADC nonlinearities. The red curve is lower by a factor of 5 or more, showing that dipole distortions due to ADC nonlinearity are a significant contributor to the leakage. The black dashed line corresponds to the fiducial model FEE with τ = 0.066. 
This demonstrates that the detected leakage at 100 GHz contains a potentially significant, and possibly dominant, spurious detection due to a degeneracy between the dipole distortion and the leakage templates. This is in line with the low level of the leakage and the level of the dipole distortion discussed in Sect. 2.5. This can only be confirmed through E2E simulations and comparison with data residuals from an appropriate null test. The remaining systematic effects are at a level of 8 × 10^{4}μ for ℓ > 3. The quadrupole systematic term is an order of magnitude larger. As discussed before, quadrupole and octopole systematic effects dominate in Fig. 9. The residuals of the ADC nonlinearity dipole distortion at ℓ > 3 are small with respect to the FEE signal.
Figure 16 shows results for 100, 143, and 217GHz data. All spectra show a rise similar to the one seen for the blue line in Fig. 15, which simulates the ADC nonlinearity dipole distortion at 100 GHz. The higher frequencies also exhibit negligible temperaturetopolarization leakage, but have comparable ADC nonlinearity dipole distortions. This is a strong indication that, at these frequencies, the dominant residual systematic is also due to ADC nonlinearity acting on the dipole, which is not removed from the maps, in agreement with Fig. 17. This effect has to be accounted for in the likelihood before science results can be extracted from the pre2016 maps. The temperaturetopolarization leakage is at a level of 10^{3}μK^{2} or lower at ℓ > 4. The dominant systematic at 353 GHz, in contrast, is calibration uncertainty (see Sect. 2.6.3).
Fig. 16 EE auto and crossspectra of the global fit test of the temperaturetopolarization leakage, for 100, 143, and 217 GHz. Levels are similar to that for 100 GHz, for which Fig. 15 shows that this is dominated by the ADCnonlinearity dipole distortion. The dashed black line corresponds to the fiducial model FEE with τ = 0.066. 
3.2. Null tests
In this section we use null tests on power spectra, together with our understanding of the simulated systematic effects discussed in Appendix B, to demonstrate that we have identified all detectable systematic effects. In Fig. 17, the green lines show the total of noise and systematic effects. Noise (red lines) dominates at multipoles ℓ > 10 for all frequencies. For the CMB channels (100, 143, and 217 GHz) the main systematic is the ADCnonlinearity dipole distortion (purple lines), which is not removed in the mapmaking; for this systematic, it is not the residual but the effect itself that remains in the maps. The residuals from other effects of ADC nonlinearities are smaller (dark blue lines), and the bandpass leakage residual (light blue) is smaller still. Frequencyband intercalibration residuals (orange) are even lower, as are other systematic effects that are not shown.
For the 353GHz channel, which is used only to clean foreground dust, the simulated systematic effects at low multipoles are all negligible with respect to the FEE model level scaled by the foreground correction coefficients for 100 and 143 GHz (dashed and dotted black). Although the intersurvey calibration difference shown in Fig. 12 is greater for 353 GHz than for the CMB channels, this does not affect the present results.
We perform null tests on four different data splits. The first is the odd−even survey differences. These map differences test residuals associated with the scanning direction and far sidelobes, short and long time constants in the bolometers, and beam asymmetry. The only clear evidence for any residual systematic effects associated with scanning direction is the odd−even solar dipole amplitude oscillation displayed in Fig.12. This is above the noise level only at 353 GHz (Sects. 2.6.2 and 2.6.3). Other scanningdirection systematic effects are almost entirely eliminated by the SRoll mapmaking algorithm, and are not discussed further in this section.
The other three data splits used for null tests are by “detset” (Planck Collaboration VII 2016), by the first and second halves of each stable pointing period (Planck Collaboration VII 2016), and by the first and second halves of the mission (Planck Collaboration VIII 2016). In each case, we compute C_{ℓ} spectra of the difference between the maps made with the two subsets of the data. The spectra are rescaled to the full data set. The resulting spectra contain noise and systematic errors, and can be compared with the FFP8 simulations (Planck Collaboration XII 2016) of the Planck sky signal and TOI noise only.
Detset differences are sensitive to errors in detector parameters, including polarization angle, crosspolar leakage, detectormismatch leakage, far sidelobes, and time response.
Halfring differences are sensitive to detector noise at levels predicted by the analysis of TOI (Sect. 2.3). All systematic errors cancel that are associated with different detector properties, drifts on timescales longer than 30 min, and leakage patterns associated with the scanning strategy^{9}.
Fig. 17 Residual EE autopower spectra of systematic effects from the HFI pre2016 E2E simulations computed on 50% of the sky (colours specified in the top left panel apply to all panels). The purple line (ADC NL total residual) shows the sum of all effects associated with ADC nonlinearity. The dark blue line (ADC NL no distortion) shows the level without the dominant dipole distortion. The plots show also the FEE model (black curves). The 100GHz and 143GHz model scaled to 353 GHz with a dust SED is shown as dashed and dotted lines, respectively. 
Halfmission differences are sensitive to longterm time drifts, especially those related to the position of the signal on the ADC and related changes of the 4K lines affecting the modulated signal.
Fig. 18 EE autopowerspectra of detectorset, halfmission, and halfring difference maps for 100, 143, and 217 GHz. C_{ℓ} rather than D_{ℓ} is plotted here to emphasize low multipoles. Results from the Planck 2015 data release are on the left; results from the pre2016 maps described in this paper are on the right. Colourcoding is the same for all frequencies. We also show for reference an average of FFP8 simulations boosted by 20% to fit the halfring null tests. The black curves show the FEE model. 
Fig. 19 EE and BB crosspower spectra of the residual effect computed from null tests between halfmission maps based on full mission minimization (red curve) or independent minimization for each half mission (blue curve). This second approach clearly shows a systematic effect. The sum of all systematic effects, dominated by the ADCnonlinearity dipole distortion shown in green in Fig. 17, is at the same level as the simulated null test with independent minimization. The FFP8 power spectrum is given again for reference. 
Fig. 20 EE and BB spectra of the 2015 maps and the pre2016 maps used in this work at 100, 143, 217, and 353 GHz. The crossspectra of detset and halfmission maps and the autospectra of the detset and halfmission difference maps are shown. The maps are masked so that 43% of the sky is used. 
The pre2016 E2E simulations can be used to distinguish whether systematic effects are comparable to the level of the base ΛCDM EE spectrum and could therefore affect the τ determination, or whether the systematic effects are negligible or accurately corrected by SRoll.
Figure 18 compares EE autospectra of detset, halfmission, and halfring nulltest difference maps for the pre2016 data and the 2015 release data. The halfring null test results (blue lines) agree with FFP8 as expected. For detset (red lines) and half mission (green lines) null tests, the 2015 data show large excesses over the FFP8 simulation up to ℓ = 100. In contrast, the HFI pre2016 data detset differences for 100 and 143 GHz are in good agreement with the FFP8 reference simulation. This is no surprise as systematics detected by this test have been shown to be small. At 217 GHz, the detset test is not yet at the level of FFP8.
For the halfmission null tests, the analysis of systematic effects shows that the ADCnonlinearity dipole distortion, which has not been removed, dominates, and should leave an observable excess at low multipoles in this test, which is not however seen. Thus this null test does not agree with the systematic analysis. A possible explanation is that the destriping is done on the full mission and applied to the two halves of the mission. The correlation thereby introduced in the two halves could lead to an underestimate of the residual seen in the null test. We checked this by constructing a set of maps in which the destriping is done independently for each halfmission. Figure 19 shows the results of this check. The independent destriper for the two halfmissions (blue lines) shows a systematic effect at all frequencies in the halfmission null test that is not seen for the full mission minimization (red lines). The separated minimization can also be compared to the sum of all simulated effects (green lines); it shows the expected behaviour for the 100 to 217GHz bands. We conclude that the uncorrected ADCnonlinearity dipole distortion accounts for most of the systematic detected by this new halfmission null test. This, of course, does not imply that the maps should be constructed using the separate halfmission minimizations because these use fewer redundancies than the full mission ones, and are therefore less powerful.
At 353 GHz, the null test is below the sum of all systematic effects between multipoles 2 and 50, and is not catching all systematic effects. The calibration and transfer functions for 353 GHz are not captured by the halfmission null test; this can plausibly explain the difference.
We have shown that the destriping of the two halfmissions should be done independently in order to detect ADC nonlinearity residuals. The 2016 data release will therefore include halfmission SRoll fits of parameters.
Figure 20 shows EE and BB spectra of the difference maps and crossspectra (signal) over 43% of high latitude sky for both the 2015 HFI maps and for the HFI pre2016 polarization maps described in this paper, and shows the relative improvements made by the SRoll processing. Differences between 2015 and pre2016 power spectra are not particularly sensitive to the polarization mask. The crossspectra show the total signal level, dominated by polarized dust emission at ℓ ≲ 200 and, in EE, by the CMB at ℓ ≳ 200. This allows a direct comparison of the signal with noise plus systematic effects.
In Fig. 20, the 353GHz detset null test for the 2015 data (blue line) at 3 ≤ ℓ ≤ 55 is 30 times larger than the FFP8 noise, and is at a level larger than 10% of the dust foreground spectrum. In the 2015 data release, systematic effects in the 353GHz maps constitute the main uncertainty in the removal of dust emission from 100 and 143 GHz at low multipoles, dominating over statistical uncertainty in the dust removal coefficient (around 3%, see Sect. 5.2). In the pre2016 data detset differences (green line), the systematic effects are much lower, but not yet at the TOI noise level.
In summary, all known systematic residuals have been seen in at least one null test at the expected level. Conversely, there is no excess over noise seen in a null test that is not accounted for by a known systematic. This important conclusion fulfills the goal of this section. However, one systematic effect has not been corrected at all, namely the ADCinduced dipole distortion because it requires a better ADC model that can be applied at the TOI or ring level simultaneously with the correction of other systematic effects. This will be done in the next generation of data.
4. LFI lowℓ polarization characterization
The measurement of τ for the Planck 2015 release was based on the LFI 70GHz maps, cleaned of synchrotron and dust emission with 30GHz and 353GHz templates, respectively (Planck Collaboration XIII 2016). In this paper we use the 70, 100, and 143GHz polarization maps to calculate 70 × 100 and 70 × 143 crossspectra (see Sect. 6). Since the LFI and HFI instruments are based on very different technologies, the systematic effects they are subject to are largely independent. In particular, the dominant systematic effects in the two instruments (gain uncertainties for LFI and residual ADC effects for HFI) are not expected to be correlated. Therefore LFI × HFI crossspectra provide a crosscheck on the impact of certain systematic effects in the estimate of τ. There can be common mode systematics as well as chance correlations, however, so the crossspectra cannot be assumed to be perfectly free of systematic effects.
In this analysis we use LFI data from the 2015 release. A detailed discussion of the systematic effects in those data is given in Planck Collaboration III (2016). To assess the suitability of the LFI data for lowℓ polarization analysis, we analyse residual systematic effects in two ways, first using our model of all known instrumental effects (the “instrumentbased” approach; Planck Collaboration III 2016), and second using null maps of measured data as a representation of residual systematic effects (the “nulltestbased” approach). In each case, we evaluate the impact of systematic effects on the extraction of τ by propagating them through foreground removal, power spectrum estimation, and parameter extraction (Fig. 21).
Fig. 21 Schematic of the simulation plan to characterize the LFI polarization data at low multipoles. Both the instrumentbased and nulltestbased strategies are represented (see text). The lower left part of the diagram outlines the crossspectrum simulation analysis involving LFI and HFI data, including systematic effects (see Sect. 6). 
We then use our instrumentbased simulations to support a crossspectrum analysis between the LFI 70GHz channel and the HFI 100 and 143GHz channels. A summary of the effects that are most relevant for the present analysis is given in Appendix C.
4.1. Instrumentbased approach
Following Planck Collaboration III (2016), we produce a map of all systematic effects (the “systematics template”) at 70 GHz and add this to realizations of the fullmission noise and the CMB from FFP8. To quantify the impact of systematic effects in the foreground removal process, we use the corresponding systematic effects templates at 30 GHz and 353 GHz, and scale them to 70 GHz with spectral indexes α_{30/70} = 0.063 and β_{353/70} = 0.0077 (Planck Collaboration XI 2016). This is equivalent to assuming that our cleaning procedure leaves no synchrotron or dust contamination in the final 70GHz map, while we evaluate the impact of the rescaled noise and systematic effects. We then extract the power spectra for the temperature and polarization components at ℓ < 30. The spectra are calculated over the same sky region used to derive τ in Planck Collaboration XI (2016).
To calculate the bias introduced by systematic effects on τ, r, and log (A_{s}), we compare 1000 FFP8 realizations of the polarized CMB (for a fiducial value τ = 0.065), plus white and 1 /f noise including systematic effects, with 1000 similar simulations containing only CMB and noise. For each realization, we calculate the marginalized distributions for each of the three parameters X (= τ, r, or A_{s}), and calculate the differences ΔX = X_{syst.}−X_{nosyst.}, which represent the bias introduced in the estimates of parameter X by the combination of all systematic effects.
For the optical depth, we find a mean bias ⟨ Δτ ⟩ = 0.005 (Fig. 22), or approximately 0.25 times the standard deviation of the value of τ measured by LFI (Planck Collaboration XIII 2016). The positive measured bias may suggest a slight overestimate of the measured value of τ caused by systematics. We also tested the sensitivity of the measured bias with the fiducial value of τ. We repeated the analysis replacing the FFP8 CMB simulations with a set of realizations drawn from a cosmological model with τ = 0.05 (all other parameters were fixed to the FFP8 values) and find no measurable difference from the previous case: ⟨ Δτ ⟩ = 0.005.
In the above analysis, the values of the scaling coefficients α_{30/70} and β_{353/70} were fixed at the values measured in the real data, as appropriate in quantifying the impact of residual systematic effects in the released LFI maps and likelihood. For a more general assessment of the componentseparation and parameterestimation procedures in the presence of systematic effects, these quantities should be reestimated for each individual realization (e.g., Sect. 2.5 of Planck Collaboration XIII 2016).
Fig. 22 Bias on τ due to systematic effects, specifically showing the distribution of τ with (red/dashed) and without (blue/solid) the systematics template. The vertical dotted line shows the input value to the simulation, τ = 0.065. 
4.2. Nulltestbased approach
Fig. 23 Power spectra (D_{ℓ} in μK^{2}) of systematic effects at 30 GHz (left) and 70 GHz (right), with each effect coded as indicated in the legend. 
The above analysis only probes systematic effects that are known a priori, and is limited by the accuracy of our instrument model. Furthermore, Fig. 23 shows that the amplitude of systematic residuals at large scales is comparable to the noise at 30 GHz, and somewhat lower than the noise at 70 GHz. Therefore tests are also needed that do not rely on a model for the systematic effects. For this purpose, we use null maps of measurement data as “templates” of the LFI systematic effects, add to them simulated CMB realizations, and extract the value of τ. This nulltestbased approach has the advantage of using real data, instead of simulations, but it assumes that difference maps contain a level of contamination that is representative of that in the full maps. If systematic effects are correlated between Surveys or years, however, under or overestimates of some effects could result.
We test this assumption using our instrumentbased systematics model, comparing the power spectrum of the systematic effects predicted by our model for the full mission with the power spectra predicted for two combinations of halfmission null maps. Figure 24 shows results for the 70GHz channel. At ℓ< 15, the fullmission spectrum is higher on average by a factor of about 2.5 than the spectra of the difference maps, likely due to partial cancellation of common modes. However, the comparison indicates a reasonable consistency of the amplitudes of the three power spectra, within the observed scatter, at a level ≈ 10^{3} for ℓ> 2. This indicates that applying the nulltestbased approach to the real data provides a useful test of LFI data quality.
Fig. 24 Power spectra for 70 GHz of systematic effects from our instrumentbased systematics simulations (with no noise) for: the full 4year mission in green; the map difference (yr1 + yr2)−(yr3 + yr4) in red; and the map difference (yr2 + yr4)−(yr1 + yr3) in blue. 
Next, we compute the power spectra for the full mission and for the halfmission differences using the measurement data. Residuals in yearmap differences (e.g., year 1−year 2) are more consistent with noise than residuals in Survey differences, as expected from the increased efficiency of the destriping algorithm (Planck Collaboration VI 2016); however, we find only marginal improvement when using 2year map differences compared to combinations of singleyear maps. Thus we adopt the latter, which allows us to evaluate the two halfmission combinations where (3)is a weight based on the number of hits per pixel, N_{hits}, needed to equalize the noise in the two maps. We apply the same mask used for extraction of τ in the instrumentbased simulations. In order to correct for the mask modecoupling and the polarization leakage effect we used CrossSpect, a MASTERlike (Hivon et al. 2002) power spectrum estimator.
Figure 25 shows the 70GHz autospectra of the difference maps D_{12−34} and D_{24−13} using real data, after removal of the noise bias. For comparison, in the same figure we also show the same spectra of Fig. 24 from simulated systematic effects (we note the different scale). As expected, the scatter for the data is much larger than the scatter for the models of systematic effects. To verify whether the scatter observed in the data is consistent with instrument noise, we compute statistical error bars at each multipole from 1000 realizations of FFP8 noise. Over the multipole range 2 ≤ ℓ ≤ 29, we find reduced χ^{2} values of 0.94 for D_{12−34} (PTE = 0.55), and 1.06 for D_{24−13} (PTE = 0.38), accounting for ℓtoℓ correlations via the empirical noise C_{ℓ} covariance matrix. Considering only the lowest multipoles, 2 ≤ ℓ ≤ 9, we find reduced χ^{2} values of 0.61 (PTE = 0.77) and 1.25 (PTE = 0.26), respectively.
Fig. 25 Power spectrum of the measured 70GHz null maps for the differences D_{12−34} (orange) and D_{24−13} (light blue). The scatter is dominated by the measurement noise. The error bars are computed from 1000 Monte Carlo simulations from FFP8 noise. Negative values are due to the fact that the noise bias has been removed. For comparison we also plot the systematicsonly spectra (red and blue). 
The fact that the scatter is consistent with FFP8 noise confirms that the amplitude of systematic effects, while not negligible, is smaller than white and 1 /f noise, even at large scales. This level of systematic effects could still lead to small biases in the cosmological parameters, so we perform an analysis of the recovery of τ in the presence of systematic effects also for the nulltestbased templates. Since these are based on a single realization of the real data, we cannot use the approach of Sect. 4.1 and marginalize over the CMB and noise realizations to obtain the average systematic impact on parameters. Rather, we proceed similarly to what was done in the above analysis of nullmap power spectra. We start by calculating the D_{12−34} and D_{24−13} yearly differences for the 30 and 70 GHz instrumentbased systematic templates discussed above. We then generate two sets of 1000 CMB+noise and CMB+noise+simulated null maps, and analyse them as in the previous section. The resulting average bias on τ is ≈ 0.002 and ≈ 0.003, for the simulated D_{12−34} and the D_{24−13} combinations, respectively. In agreement with the power spectrum results of Fig. 24, we find that the null templates give rise to a lower bias than the full mission systematics templates.
Fig. 26 Bias on τ due for the nulltestbased systematics templates. The blue line shows the distribution of τ for the reference CMB+noise simulations, the red line shows the resulting distribution when including also the null systematic maps built from the instrumentbased simulations. The vertical magenta line shows the value of τ estimated for the single CMB realization to which the null maps from actual data, as opposed to the null maps from the instrumental simulations, were added. The dotted line shows the input value to the simulation, τ = 0.065. Top: D_{12−34} null combination. Bottom: D_{24−13} null combination. 
Fig. 27 EE and BB differences and crossspectra of the 2015 maps at 30, 44, and 70 GHz, for different data splits. These maps are masked so that 43% of the sky is used. This figure show similar plots to Fig. 20 for LFI frequencies. 
We then add a single CMB realization to the nulltestbased templates, and estimate parameters for these maps as well. Fig. 26 shows that the resulting estimates are well within the distributions. Note also that the bias in the case of D_{12−34} is somewhat lower than in the case of D_{24−13}, consistent with the behaviour expected from simulations (Fig. 24). These results suggest that the instrumentbased and nulltestbased null maps have a consistent impact on power spectra and parameters, and provide support for the use of the instrumentbased simulations in the crossinstrument analysis.
Figure 27 shows EE and BB pseudopowerspectra from the 2015 LFI data release used in this paper, showing the difference maps and crossspectra (signal) over 43% of the sky. Foregrounds are clearly visible in the crossspectra only at ℓ < 30 at 30 GHz and ℓ < 5 at 70 GHz. As will be shown in Sect. 5.2, although the removal of the modest synchrotron foreground component at 70 GHz may not be extremely accurate, its final uncertainty is well within the noise.
5. Interfrequency calibration, component separation, and power spectra
The lowmultipole polarization signals measured by Planck include Galactic foregrounds. In this section, we describe the combination of Planck frequency bands that we use to separate the CMB signal from these foregrounds. First we revisit dipole residuals in the singlefrequency maps to establish the level of precision of the intercalibration error between the frequency bands of Planck in both instruments. Then, we use the SMICA componentseparation code to show that the polarization signal at large angular scales consists mostly of CMB and two Galactic foreground components. A third component does not improve the fit, and we show that a simple internal linear combination (ILC) method is sufficient to clean the CMBdominated bands by regressing synchrotron with 30 GHz and dust with 353 GHz. Finally, we describe the angular power spectra at low multipoles using two estimators and subtracting the foreground components.
Solar dipole amplitude and direction by frequency, as well as for the average (AVG) of the 100 and 143GHz maps.
5.1. Interfrequency calibration
For the 44, 70, 100, 143, 217, and 353 GHz channels, the photometric calibration is based on the “orbital dipole”, i.e., the modulation of the solar dipole induced by the orbital motion of the satellite around the solar System barycentre (see details in Appendix B). The “solar dipole”, induced by the motion of the solar system barycentre with respect to the CMB, is stronger by an order of magnitude, and can be measured after masking the Galactic plane and removing foreground emission outside the mask. If the uncertainties associated with foreground removal are negligible, the amplitude of the solar dipole provides an excellent check on the relative calibration between different frequencies.
Most of the higherorder terms discussed in Notari & Quartin (2015) affect all CMB dipoles in the same way, and thus do not induce errors in the relative calibration of the solar dipole with respect to the orbital one but can contribute small factors to the frequency calibration differencies.
The frequencydependent kinematicdipoleinduced secondorder quadrupole and dipole from foreground monopoles (Kamionkowski & Knox 2003) have amplitude 2.0 × 10^{4}μK at 100 GHz, 3.9 × 10^{4}μK at 143 GHz, 8.2 × 10^{4}μK at 217 GHz, and 1.7 × 10^{3}μK at 353 GHz. These corrections have a negligible impact on the calibration of the CMB channels.
Table 2 summarizes the HFI measurements of the solar dipole from 100 to 353 GHz. At 100 GHz, the direction and amplitude of the solar dipole remain within two standard deviations as the sky fraction changes from 37% to 58%. At 143 GHz, the direction is the same as at 100 GHz within 2σ of the statistical noise. Higher frequencies show drifts in direction with changing sky fraction that are larger. The amplitude also drifts when the sky fraction is reduced.
The excellent stability of the solar dipole amplitude and direction with changes in the sky fraction implies that the orbital dipole calibration at each frequency is better than 0.1% for Planck CMB frequencies and WMAP Wband, as shown in Fig. 28 and Table 3, which give the properties of the solar dipoles measured at Planck frequencies and in the WMAP bands.
Table 3 gives the solar dipole and the relative amplitude as measured by Planck from 44 to 353 GHz, plus that of WMAP at 94 GHz, as well as the relative calibration determined from the first and second peak amplitude of the CMB spectrum. The relative amplitudes determined both ways are plotted in Fig. 28. The relative calibration determined from the first two peaks (specifically, over 110 ≤ ℓ ≤ 500) of the CMB power spectrum, is very close to that determined from the solar dipole. The difference in calibration between the solar dipole and the first two CMB peaks gives an upper limit on the residual transfer function error relative to 100 GHz, between ℓ = 1 and ℓ ≈ 300, also shown in Table 3.
Fig. 28 Relative calibration based on measurement of the first two acoustic peaks (blue) and on the solar dipole (red). The bottom panel is a zoom of the top panel. We use 100 GHz as the reference frequency for the dipole calibration method. In this plot, Planck 30GHz and WMAP Qband data may be affected by their relatively low angular resolution. 
Relative calibration between the solar dipole and the first and second peaks (ℓ = 110–500) in the CMB power spectrum with respect to 100 GHz.
The 545GHz channel is difficult to calibrate on the orbital dipole because the orbital dipole is so weak with respect to the dust emission. Instead, the 545GHz channel was calibrated on Uranus and Neptune (Planck Collaboration VIII 2016). Nevertheless, the (much stronger) solar dipole and acoustic peak amplitudes can be compared with those of the CMB channels, once again after cleaning dust using the 857GHz template. The result is that the planet calibrations agree with the CMB calibrations within 1.5% for 545 GHz. This implies that the 857 GHz is within 2.5% of the CMB calibration considering the uncertainty their respective planet calibration. The absolute calibration uncertainty using planets was estimated at 5%, and is now shown to be within 2% of the photometric CMB calibration based on the orbital dipole and the COBE/FIRAS absolute calibration of the CMB temperature (Fixsen et al. 1997; Fixsen 2009). Since the Herschel observatory was calibrated with the same planet models, intercalibration between Herschel and Planck is also at the 2% level (Bertincourt et al. 2016).
In summary, absolute calibration of the 100GHz channel on the orbital dipole is insensitive to sky fraction (and thus not affected by Galactic foreground removal). Based on that absolute calibration, we use the solar dipole to compare the calibration of the other Planck CMB channels and WMAP Wband. We also compare the calibration of the CMB channels using the CMB anisotropies in the multipole range 110–500. There is evidence of some potential orbital dipole calibration errors due to the dust foreground at frequencies above 100 GHz, but also of very stable behaviour of the transfer functions with Planck frequencies from 70 to 217 GHz. At 353 GHz, we still see a 0.3% discrepancy between various different calibrations, and for WMAP W band a 0.5% difference. Whether the origin of these discrepancies is in the calibration itself or is a result of transfer function errors cannot be determined by this analysis.
5.2. Diffuse component separation
In this section we discuss the effects of polarized diffuse component separation on the CMB spectrum (point sources can be ignored for lowmultipole work). To begin, we use the (blind) SMICA code (Planck Collaboration XII 2014), with no assumptions about the number and properties of the foreground components, to establish the number and properties of the foreground components that must be removed, and the fraction of sky f_{sky} to use. It turns out that two polarized foreground components are enough for polarization, i.e., including a third foreground component does not improve the fit. Not surprisingly, the two components are easily identified with dust and synchrotron emission. Changing f_{sky} from 0.4 to 0.5 makes essentially no difference in the resulting CMB power spectrum, but changing to 0.6 or 0.7 makes a substantial difference. This leads us to adopt f_{sky} = 0.5.
The two likelihood methods described in Sect. 6 use internal linear combination (ILC) or template fiting to remove synchrotron and dust emission. The 30GHz map is used as a template for synchrotron emission that is uncorrelated with dust, and the 353GHz map is used as a template for dust emission. Table 4 shows the “projection coefficients” used by the different likelihoods. Table 4 also gives the SMICA projection coefficients for two blind components, and shows that the diffusecomponentseparation procedure is very stable between the different approaches for 100 and 143 GHz. At 70 GHz, the dispersion is larger, but if no synchrotron component were explicitly removed, the synchrotron that is correlated with dust would still be taken care of; additionally, the uncorrelated part would have no effect to first order on the 100 × 143 GHz crossspectra because the synchrotron signal is very weak at 143 GHz, as shown in Fig. 29.
Fig. 29 Dust correction to C_{ℓ} in μK^{2}, using the 353GHz channel as a template (in blue), along with the synchrotron correction to C_{ℓ} in μK^{2} using the 30GHz channel as a template (in red). These are plotted for the 70–217 GHz channels as a function of frequency. We note that the points for synchrotron at 143 and 217 GHz are very low and can only be considered as upper limits. 
Table 5 gives the average value of the power spectrum removed for each foreground at ℓ = 4, at the peak of the EE reionization feature, together with the associated uncertainties computed with the ILC method. The uncertainties are always lower than the FEE signal, by more than an order of magnitude. Figure 29 shows this graphically.
At the two Planck frequencies with the lowest noise, 100 and 143 GHz, the level of synchrotron emission is lower than that of dust by factors of 10 and 300, respectively. For 100 GHz, not removing the synchrotron would introduce a bias of less than 6% of FEE. The 70GHz spectrum is limited by noise rather than the large uncertainty (16%) in the synchrotron foreground removal, even though this is the largest foreground removal uncertainty in the Planck CMB channels for the reionization peak.
For each CMB frequency, the amplitude of the power spectrum D_{ℓ} at ℓ = 4 that is removed for the dust and synchrotron foregrounds.
The spatial variation of the dust spectral energy distribution at high latitude has been analysed on 400 deg^{2} patches, and is smaller than the ILC uncertainty (Planck Collaboration Int. XXX 2016).
Figure 30 shows the level of residual foregrounds. It is clear that the componentseparation errors propagated to the power spectra are much lower than the sum of all systematic effects, and have a negligible impact on the determination of τ.
5.3. Power spectra
Likelihood analyses for CMB cosmological parameters were originally developed for CMB experiments that were dominated by CMB signal and detector noise, in contrast to more recent experiments with much lower noise in which residuals from instrumental systematic effects and foreground subtraction dominate over the detector noise, especially at low multipoles. Such analysis approaches rely on pixelpixel covariance matrices coming from a large number of simulations of CMB and noise, which are assumed to be two statistically independent Gaussian fields for the purposes of these likelihood codes. Pixelbased likelihoods can also account for noise correlations and masking of the sky. In the Planck likelihoods, the “noise” description combines the detector and readout chain white noise, some correlated fraction due to 1 /f noise, and residuals of systematic effects that have statistical properties that are not well understood. These systematic effect residuals thus need to be brought to levels much lower than the Gaussian noise.
Fig. 30 Residual errors in EE from component separation, estimated from the scatter of the componentseparation coefficients. The fiducial EE spectrum and the noise plus systematics residuals (green line) are also shown. 
5.3.1. Crossspectra at very low multipoles
⟨ D_{ℓ} ⟩ _{rms} over 2 ≤ ℓ ≤ 8 for autospectra simulations of the ADCinduced dipole distortion and noise at 100 and 143 GHz (Fig. 17), together with simulations of 100 × 143EE crossspectra (Fig. 31) and data BB crossspectra (bottom panel of Fig. 33) tracing final residuals and noise (signal is negligible).
Figure 9 shows that systematic effects in the 100GHz channel are below 10^{3}μK^{2} for 4 <ℓ < 100. At ℓ = 2 and 3, the ADCinduced dipole distortion dominates the systematic effects in simulated autospectra (Fig. 17). Figure 31 shows that the total level of simulated systematic effects (purple line) in the 100 × 143 crossspectrum is significantly reduced compared to the autospectra of both frequencies (Fig. 17) for most multipoles relevant for the reionization feature. In the first two columns of Table 6 we give ⟨ D_{ℓ} ⟩ _{rms} for 2 ≤ ℓ ≤ 8 for the autospectrum of the ADCinduced dipole distortion and noise from simulations (Fig. 17). The third column gives the same ⟨ D_{ℓ} ⟩ _{rms} for the simulated ADC induced dipole distortion EE crossspectra 100 × 143 (Fig. 31) tracing the final residuals. The fourth column gives this same quantity for the BB data crossspectrum (bottom panel of Fig. 33) for noise and residual systematics after removal of the average of the simulation of systematic effects (where the expected signal is very small) for QML2.
The steps in the reduction of the ADC residuals from the simulation autospectra (9 × 10^{3} and 13 × 10^{3}μK^{2}) to the crosspower spectrum (6 × 10^{3}μK^{2}) and finally the QML BB data crossspectrum of residual systematics and noise (3 × 10^{3}μK^{2}) lead us to only a small excess with respect to the expected noise (3.0 × 10^{3} to 2.0 × 10^{3}μK^{2}).
Fig. 31 Similarly to Fig. 17, residual EE crosspower spectra of systematic effects from the HFI pre2016 E2E simulations computed on 50% of the sky are shown for 100 × 143. Dotted lines are used to join two multipoles when one is negative. 
Crossspectra calculated from two different detsets within a single frequency, for both 100 and 143 GHz, have a higher level of residual systematics than those from crossfrequency detsets, as discussed below. As explored in Sect. 5.2, 100 × 143EE crossspectra can be cleaned of dust and of the correlated synchrotron fraction (60%) with the 353GHz template, and have negligible synchrotron power remaining because synchrotron emission is so weak at 143 GHz. We also described in Sect. 5.2 why we adopted f_{sky} = 0.5. For these reasons, we choose the full frequency 100 × 143 crossspectra and f_{sky} = 0.5 as our baseline.
5.3.2. Estimators
We use both pseudoC_{ℓ} (PCL) and quadratic maximum likelihood (QML) estimators to extract crossspectra from maps. Two versions of the PCL crossspectra were produced, using Xpol (Tristram et al. 2005) and Spice, a MASTERlike code (Hivon et al. 2002).
QML estimators have been widely discussed in the literature (e.g., Tegmark 1996; Bond et al. 1998; Efstathiou 2004, 2006). QML autospectra are close to optimal, but must be corrected for noise bias. The removal of this bias requires an accurate estimate of the pixel noise covariance matrix N_{ij}. In analogy with QML autospectra, it is straightforward to define a quadratic crossspectrum estimator that is unbiased. The resulting crossspectrum estimator will not have minimum variance, but nevertheless we retain the nomenclature “QML” to distinguish it from the PCL estimators described above.
For two maps “a” and “b,” the crossspectrum estimate is defined as (4)where i and j are pixel numbers, and ℓ and r form a paired index with ℓ denoting multipole number and r denoting spectrum type (e.g., TT, TE, EE or BB). The matrix E here is (5)where the covariance matrix C = S + N and is a “reshaped” covariance matrix of the form (6)that does not mix temperature estimates with polarization estimates. Provided that the noise between maps a and b is uncorrelated, the expectation value of Eq. (4) is (7)where (8)Equation (7) can be inverted to give a deconvolved estimator of C_{ℓr}. The variance of the crossspectrum estimator can be written in terms of the matrices S, N, and E as Inaccurate determinations of the noise covariance matrices N^{a} will therefore not bias the power spectrum estimates, but will lead to inaccurate estimates of the variance, via Eq. (10).
In Sects. 2 and 3 we discussed systematic effects in full resolution maps. Most residual systematic effects are negligible above ℓ = 50. For science analysis at ℓ< 20, we degrade the resolution to HEALPix N_{side} = 16 with the following smoothing: (11)This does not affect the analysis of systematic effects.
5.3.3. Bias estimate
To estimate the bias induced by residual systematic effects at the power spectrum level, we average simulations HFPS1, HFPS2, and HFPS3 (Sect. 2.2) and calculate crossspectra on the average data with both estimators (Fig. 32). These averages are then removed from the data. We note that the biases and their uncertainties are relatively small compared to the theoretical spectra shown in Fig. 17. This is due to the use of crossspectra, as mentioned above.
Fig. 32 100 × 143 crossspectra of noise and systematics from the average of HFPS1, HFPS2, and HFPS3 simulations, calculated using both PCL and QML estimators (with the covariance matrix from HFPS2). Theoretical spectra are plotted (in black) for τ = 0.05, 0.07, and 0.09. 
5.3.4. Building the crossspectra
For QML estimators, full mission pixelpixel noise covariance matrices at N_{side} = 16 were produced from the 2015 FFP8 simulations, which capture some aspects of correlations via pixel hit counts and the Planck scanning pattern, but which do not include errors from instrumental systematic effects. To account for instrumental systematic effects, a covariance matrix is derived from the HFPSs. The signal covariance matrix assumes the 2015 base ΛCDM model parameters with τ = 0.07. After adding CMB signal, this matrix can be inverted. We used the HFPS1 (83 realizations) and HFPS2 (100 realizations) simulations. One set is used to build the pixel covariance matrix. We compute the QML spectra using either the same set or the other one, to test the effect of overfitting and noise bias introduced by using the same set when the number of simulations is small. This bias is demonstrated in Appendix D, which shows that using the same set of simulations produces significant effects in the dispersion of the QML simulated spectra, as well as more discrepant PTEs. These biases decrease as the number of simulations increases. As an extra check, we swap the two independent sets and find fully compatible results.
Fig. 33 100 × 143 crossspectra used in this paper. Top: debiased XPol PCL (green), and the biased (red) and debiased (blue) SimBaL PCL spectra. Middle: QML power spectra (red biased, blue debiased) for HFPS1. Bottom: same for HFPS2. Model spectra for τ = 0.050, 0.070, and 0.090 are displayed in black (dashed) lines. 
The 100 × 143 crossspectra for science analysis are displayed in Fig. 33. Debiasing is shown for the SimBaL pseudoC_{ℓ} spectra, and is significant only at very low multipoles. The debiased Xpol and SimBaL PCL spectra are fully consistent.
The QML power spectra calculated from HFPS1 and HFPS2 differ only slightly. A quantitative estimate of the implications when propagated all the way to τ is given in Appendix D. The BB power spectra from both codes show negligible signal, as expected, and can be used as an estimate of the noise and an upper limit on residual bias.
The QML estimator produces more local estimates, i.e., with less covariance between multipoles, and better isolation of E and Bmodes. The latter is clearly seen in Fig. 33. We also expect the QML approach to give more optimal estimates than the PCL approach, and indeed the QML error bars are found to be significantly smaller than the PCL ones. In both cases, the impact of debiasing the spectra for the ADC nonlinearity distortion is small compared to the uncertainties.
Figure 34 compares statistical distributions of expected values of the EE crossspectra derived from the HFPSs to observed values (vertical lines), and gives the associated probabilitytoexceed (PTE) those values. For 2 ≤ ℓ ≤ 5, we see nonGaussian and asymetrical probability distributions coming from systematic residuals and cosmic variance. Figure 35 shows the PTE expressed as the equivalent χ^{2} of a Gaussian distribution, for a broader range of multipoles. The PCL and QML estimators are largely in agreement. There is a significant outliers: ℓ = 16 which is bad for both estimators. The PTEs are similar, demonstrating that the results from both procedures are in agreement with the expected statistics. Since the QML estimate gives significantly smaller error bars (Fig. 34), we use it for the final results.
Fig. 34 Expected EE crossspectrum statistical distributions for 2 ≤ ℓ ≤ 7, computed from the HFPS1 simulations, with a fiducial τ = 0.055. The blue (PCL) and red (QML) vertical lines show the observed values. The probability to exceed (PTE) is given in each panel. 
Fig. 35 PTE expressed as equivalent χ^{2} for both the QML and PCL estimators, derived from the statistic shown in Fig. 34. 
5.3.5. TE crossspectra
The TE crosscorrelation spectrum can in principle also contribute to the determination of τ, and has been used when the EE spectra were either too noisy (e.g., in early WMAP results) or limited by systematic effects (in earlier Planck data releases). Figure 36 shows that at present the TE spectrum is compatible with the range of τ values allowed by EE. Nevertheless the uncertainties are such that TE cannot bring any significant improvement to the determination of τ from EE, partly due to the large cosmic variance of the temperature signal at very low multipoles (see also Fig. 3 in Planck Collaboration Int. XLVII 2016). Furthermore, to include TE fully would require a more comprehensive analysis of component separation for the largescale temperature map, which has not yet been done.
Fig. 36 TT and TE100 × 143 crossspectra, plotted for the SimBaL results with and without the bias correction. The black lines shows the fiducial spectra for τ = 0.05, 0.07, and 0.09. 
6. Likelihoods for τ
6.1. Description of different likelihoods
In this section we estimate the reionization optical depth τ, using PCL and QML estimators of EE crossspectra between the two best Planck channels, 100 and 143 GHz, at very low multipoles. The τ parameter is strongly degenerate with the amplitude of the primordial spectrum A_{s}. The TT power spectrum constrains the combination A_{s}e^{− 2τ} at the subpercent level, 10^{9}A_{s}e^{− 2τ} = 1.875 ± 0.014 (Planck Collaboration XIII 2016). The τ–A_{s} degeneracy can be broken using PlanckTT data along with CMB lensing, or by combining with external data on largescale structure, both of which constrain A_{s}. However, the amplitude of the EE reionization feature depends quadratically on τ, with very little dependence on the other parameters of the ΛCDM model. A 10% constraint on τ from EE data constrains A_{s} at the 1% level (δA_{s}/A_{s} ≈ 2 δτ for τ ≈ 0.06), and can affect some of the tensions within the cosmological parameters.
We use the following two likelihoods with the crossspectra from Sect. 5.3.2 to estimate τ.

Lollipop (Mangilli et al. 2015) is based on a modification of the Hamimeche and Lewis approach (Hamimeche & Lewis 2008) for crossspectra at low multipoles. The offset is proportional to an effective noise term , derived from the HFPS set from which the covariance matrix is also computed. We use this likelihood only with PCL spectra in this paper (see also Planck Collaboration Int. XLVII 2016).

SimBaL is a likelihood code based on simulations, and targeted at estimating τ, as described in Appendix D. We apply it to both PCL and QML spectra. For QML spectra, we use either two independent subsets of HFPSs to determine the spectra and the covariance matrix, or (to limit the bias due to using the same small set of simulations) we use the full set, but with only a few eigenmodes to describe the systematic effects.
For QML, we define three versions of SimBal. SimBaL1 and SimBaL2 use covariance matrices determined from the simulation sets HFPS1 and HFPS2, respectively, while SimBaL3 uses the covariance matrix determined from all 283 simulations (HFPS1, HFPS2, and HFPS3) and only four eigenmodes. These three versions of the SimBal likelihood deal in two different ways with the difficulty of the limited number of simulations (as discussed in Appendix D).
The two likelihoods sample τ in the range 0.01–0.15 with a step Δτ = 0.001, and with all other parameters, except A_{s}, fixed to the Planck 2015 bestfit values. We keep A_{s} e^{− 2τ} at the fixed value from Planck 2015 (since it is tightly constrained by the higher multipoles). It was shown in Sect. 2.5 that the main systematic effect left in the maps is small, being significant only for ℓ = 2 and 3, and can be simulated and removed (see Sect. 6.2). Figure D.8 shows that taking ℓ_{min} to be 2, 3, or 4 shifts the τ posterior distributions by less than 4 × 10^{3}, confirming that the final removal of this last systematic effect (by the subtraction of the simulated effect, combined with the use of QML crossspectra) is very good. In Sect. 6.2, we test the robustness of our analysis with respect to the removal of residual systematics, and with respect to the two estimators of the power spectra described in Sect. 5.3. We also test the consistency of our results by using crossspectra between HFI and LFI, specifically 70 × 100 and 70 × 143.
Fig. 37 τ posteriors computed by SimBal from the 100 and 143GHz QML cross and autospectra, both with (solid curves) and without (dashed curves) debiasing by the mean of the simulated power spectra. While the 100 × 143 crosspower spectrum is not affected by debiasing, the 100 × 100 and 143 × 143 autopower spectra, which are still dominated by the systematic effects, are more affected. Nevertheless all debiased τ posteriors are consistent. 
6.2. τ determination from EE crossspectra only
Figure 37 shows τ posteriors computed by SimBal from the 100 and 143GHz QML cross and autospectra, both with and without debiasing. Debiasing makes essentially no difference for the 100 × 143 crossspectrum (blue full and dotted lines); however it makes a significant difference in the autospectra, as expected because the use of QML crossspectra removes part of the systematics. In the following, we always use debiased EE100 × 143 crossspectra to extract τ.
The accuracy of the SimBal likelihood in recovering τ is tested using the HFPS set and shown in Fig. 38 for both QML and PCL spectra. The input value of τ = 0.06 is recovered accurately as (QML, 83 realizations in HFPS1, with pixelpixel covariance matrix from 100 realizations in HFPS2; see Appendix D) and (PCL, 83 realizations in HFPS1). This demonstrates that our method does not remove signal and provides an essentially unbiased estimator for τ. The QML estimator produces smaller error bars, as expected, with a narrower and more symmetrical distribution.
Fig. 38 Left column: simulation of SimBaL posterior distributions using the HFPS sets for noise and systematic effects, with the fiducial value of τ = 0.060. Right column: distribution of the τ peak value from the left column. The top row is for the PCL estimator and HFPS1 simulations, and the bottom row for the QML estimator and HFPS2 simulations. 
Fig. 39 τ posteriors obtained with the two likelihood methods (Lollipop and SimBal) for the 100 × 143 QML and PCL crossspectra estimators, colour coded as indicated in the legend. 
Figure 39 shows τ posteriors for the data computed by Lollipop and SimBal from 100 × 143 PCL and QML crossspectra. For the PCL spectra, results from Lollipop and SimBal are fully consistent. The asymmetry of the posterior distribution at low τ is smaller for SimBaL than for Lollipop, due to SimBaL’s better handling of the statistics of the C_{ℓ}s (see Appendix D). For QML, the three SimBal estimates give consistent posterior distributions, with significantly narrower width than the PCL results. They therefore provide our tightest constraints on τ. The nearcoincidence of the highτ tails of the distributions would imply essentially identical upper limits on τ for the PCL and QML approaches, if the posteriors were used to provide upper limits on the reionization redshift. However, the posteriors with QML crossspectra show a clear detection of τ at the level of 3.5σ, as discussed quantitatively in Appendix D.
To test how efficiently crossspectra between the two Planck instruments suppress uncorrelated systematics, we used the HFI HFPS1 and the LFI endtoend simulations of lowmultipole systematic effects (83 for HFI, 10 for LFI). In the LFI 70GHz channel maps, residuals from calibration uncertainties dominate over systematic effects. They are at the same level as the dominant HFI systematic, but are lower than the noise (as shown in Sect. 4). To each of these simulations of systematic errors we add 100 CMB simulations, with an input value of τ = 0.05. We then calculate 70 × 100 and 70 × 143 PCL crossspectra and run the SimBaL likelihood.
Figure 40 (top) shows results from simulations. The peak values are not significantly biased with respect to the input value of τ = 0.05. This indicates that the crosscorrelation between LFI and HFI frequencies removes the simulated residual systematic effects rather well for both instruments. The distributions are asymetric and not very smooth because of the small number of LFI simulations. They show full width at half maximum 1.5 to 2 times larger than that obtained for the QML estimate using the HFI frequencies alone. Figure 40 (bottom) shows results from the data, specifically the HFI pre2016 maps at 100 and 143 GHz, and the LFI 70GHz 2015 released maps (Planck Collaboration VI 2016). As in the simulations, τ is extracted from the 70 × 100 and 70 × 143 PCL crossspectra using SimBal. The peak values and 68% upper and lower limits are Very low values are not excluded by these combinations of data, but the peak values are compatible with the baseline results from the HFI 100 × 143 crossspectrum discussed above.
Fig. 40 Posterior distributions of τ calculated by SimBaL from LFI × HFI crossspectra. Top: results from simulations using 10 LFI systematic effect realizations, 83 HFPS1 realizations, and 100 CMB realizations with input reionization parameter τ = 0.05. Bottom: results from data, specifically the HFI pre2016 maps at 100 and 143 GHz and the LFI 70GHz 2015 released maps (Planck Collaboration VI 2016). Results from 70 × 100 and 70 × 143 are consistent. 
Fig. 41 History of τ determination with WMAP and Planck. We have omitted the first WMAP determination (τ = 0.17 ± 0.04, Bennett et al. 2003), which was based on TE alone. 
6.3. Summary of results
Table 7 summarizes the results on the posterior distributions of τ based on the HFI 100 × 143EE crossspectra PCL and QML estimators shown in Fig. 39.
Peak values, 68% upper and lower limits, together with 95% upper limits, for the two likelihood methods (SimBaL and Lollipop) and the two crossspectra estimators (PCL and QML).
The QML results give a detection of τ at more than 3.5σ (see Appendix D), with the smallest uncertainties obtained so far from CMB data. For the same likelihood method and sky fraction, a cosmicvariancelimited measurement would have an uncertainty of 0.006. The peak values obtained with the PCL and QML methods agree to within 0.2σ. Crossspectra between the two Planck instruments (70 × 100 and 70 × 143) also give compatible results, but with larger uncertainties.
The most stringent results in Table 7 are obtained with the SimBal likelihood from 100 × 143 QML crossspectra. Taking a conservative uncertainty between the three QML results, we obtain (12)We will refer to this as the “lowE” data set and likelihood^{10}. There has been a significant decrease in the peak value of τ since its first determination from CMB TE measurements in 2003 and subsequent refinement using EE measurements from 2006. Figure 41 shows the history of τ estimates.
Reionization history models based on astrophysical observations of high redshift sources predict asymptotic values of τ at high redshift in the range 0.048−0.055 (Fig. 7 of Mashian et al. 2016) or 0.05−0.07 (Fig. 2 of Robertson et al. 2015). Our results are fully consistent with these expectations. For the first time, the upper limit on τ derived from CMB EE observations gives meaningful limits on how such models can be extrapolated to redshifts larger than 10.
7. Implications of a lower value of τ for cosmological parameters
The first Planck results on polarization at low multipoles (Planck Collaboration XI 2016; Planck Collaboration XIII 2016), based on the LFI maps, gave a lower value of τ than that given in the 9year WMAP analysis (Hinshaw et al. 2013). We also showed that cleaning the WMAP polarization maps for polarized dust emission using the HFI 353GHz maps led to a reduction in τ, consistent with the results from the LFIbased analysis. The results presented in this paper lead to even lower values of τ.
The amplitude of the reionization bump in the EE power spectrum at low multipoles scales approximately as τ^{2}. Low values of τ are therefore difficult to measure from largescale polarization measurements of the CMB. As this paper demonstrates, exquisite control of systematic errors, polarized foregrounds, and instrument noise are required in order to measure the small polarization signal induced by cosmic reionization. The main results presented in this paper are based on the 100 × 143EE crossspectrum, summarized in Fig. 39 and Table 7. The posteriors for τ in Fig. 39 show clear narrow peaks, with maxima at τ = 0.055, indicating detection of an EE reionization feature. Furthermore, in Table 7 we find a 95% upper limit of τ< 0.072. This limit is consistent but significantly lower than the Planck+LowP LFI results, the Planck dustcleaned WMAP results reported in Planck Collaboration XIII (2016), and the LFI × HFI results summarized in Fig. 40.
Current astrophysical observations exclude very low values of τ. If the Universe is abruptly reionized at redshift z_{re}, the optical depth caused by Thomson scattering is (e.g., Shull & Venkatesan 2008) (13)for the base ΛCDM model with a helium abundance by mass of Y_{P} (assuming that the helium remains neutral). The GunnPeterson test (Gunn & Peterson 1965; Fan et al. 2006) provides strong astrophysical evidence that the intergalactic medium was highly ionized by a redshift of z = 6.5. For the Planck 2015 base ΛCDM parameters, Eq. (13) gives τ = 0.039 for z_{re} = 6.5. This conservative astrophysical constraint on τ, which eliminates the low τ regions of the posteriors shown in Fig. 39, is in excellent agreement with our current constraint. The result also stands in good agreement with the PlanckTT+lensing+BAO constraints reported in Planck Collaboration XIII (2016), τ = 0.067 ± 0.016, which makes no use of CMB polarization at low multipoles.
Of course the true value of τ is important in understanding the formation of the first stars and the process of reionization, but it also has an impact on cosmological parameters that are of more fundamental significance. This is illustrated in Fig. 42. The principal degeneracies are between τ and A_{s} and between τ and σ_{8}. However, τ is also correlated with the spectral index of the fluctuation spectrum, n_{s}, which is set by early Universe physics. As we will see below, improving the constraint on τ helps break this particular degeneracy.
Parameter constraints for the base ΛCDM cosmology (as defined in Planck Collaboration XVI 2014), illustrating the impact of replacing the LFIbased lowP likelihood (used in the 2015 Planck papers) with the HFIbased SimLow likelihood discussed in the text.
Fig. 42 Parameter constraints for the base ΛCDM cosmology, illustrating the τ–n_{s} degeneracy and the impact of replacing the LFIbased lowP likelihood used in the 2015 Planck papers with the HFIbased SimLow likelihood discussed here. The values of τ and σ_{8} shift downwards. 
To assess the impact of our new polarization data on the cosmological parameters, we have constructed a simplified, lowmultipole likelihood based on simulations that include systematic effects. This is SimLow, based on the 100 × 143EE spectrum between ℓ = 2 and 20. SimLow is a simulationbased likelihood that gives the posterior distribution of fiducial C_{ℓ}s for each multipole; it is fully described in Appendix D.7. To explore the effects on cosmological parameters, we combine SimLow with the PlanckTT likelihood (i.e., Commander at ℓ< 30, and Plik at higher ℓ).
In the Planck 2015 analysis (Planck Collaboration XI 2016), we used an LFIbased low multipole polarization likelihood (referred to as “lowP”) in combination with low multipole temperature likelihood (referred as “lowTEB”). This gives (14)This result is mostly controlled by the lowP part of the likelihood combination. This constraint is compatible with the HFI results presented in Sect. 6, although it is statistically weaker. Combining the LFI polarization likelihood with the high multipole PlanckTT data, for base ΛCDM we obtain (15)Adding the PlanckTT likelihood drives the value of τ upwards by about 0.5σ. The A_{s}e^{− 2τ} degeneracy at intermediate multipoles (ℓ< 1500) is broken by the lensing effect seen in the higher part of the spectrum.
However, the ℓ ≳ 1000 part of the Planck spectrum is characterized by peaks that are slightly broader and smoother than what the ΛCDM model predicts. The highmultipole peak smoothing is compatible with a slightly stronger lensing amplitude, and translates into a roughly 2σhigh phenomenological parameter A_{L} value. The value derived from the lensing power spectrum (Planck Collaboration XIII 2016) supports that this would just be a statistical fluctuation, rather than a peculiar feature of the lensing power spectrum itself. Nevertheless, the preference for a larger lensing amplitude at high multipoles pushes the normalization and the optical depth values up. The lowP likelihood was not statistically powerful enough to counteract this trend, and so in the PlanckTT+lowP analysis τ is driven upwards compared to Eq. (14). This effect is discussed at length in Planck Collaboration XVI (2014) and Planck Collaboration XIII (2016).
Adding the Planck lensing measurements, which are compatible with lower values of A_{s}, drives τ down again, close to the original lowP value: (16)These shifts, and in fact the lowmultipole power deficit, are not of sufficiently high significance to suggest new physics. Moreover, the fact that adding Planck lensing causes τ to shift downwards suggests that τ is lower than the value in Eq. (14). Indeed SimLow alone gives the following constraint on τ: (17)We can anticipate what will happen if we replace the lowP likelihood with a statistically more powerful polarization likelihood favouring a low value of τ – the main effect will be to shift σ_{8} towards lower values, with a proportionately smaller shift of n_{s} also to smaller values. Furthermore this would not be consistent with solving the highmultipole peak smoothing through an underestimate of the effect of lensing. In fact, adding the Planck lensing measurement to SimLow and PlanckTT has a small impact on the value of τ, giving τ = 0.057 ± 0.0092.
Figure 42 compares the SimLow and lowP parameter constraints on τ, σ_{8}, and n_{s}, while Table 8 gives numerical results for parameters of the base ΛCDM model. The tighter constraint on τ brought by SimLow reduces the correlation between n_{s} and τ, and leads to slightly tighter bounds on n_{s}. However, larger parameter changes are seen for τ and A_{s}, each changing by about 1σ. We specifically find (18)in excellent agreement with the result from SimLow alone. The present day amplitude of the fluctuations, σ_{8}, decreases by about 1σ and its uncertainty shrinks by about 33%. This shift goes in the right direction to reduce the tensions with cluster abundance and weak galaxy lensing, as discussed in Planck Collaboration XIII (2016), although it is not yet sufficient to remove them entirely.
Changes in most of the other cosmological parameters are small, with deviations being less than 0.5σ. The largest deviation is in the spectral index n_{s}, due to its partial correlation with τ. This slight decrease in n_{s} means that we can now reject the scaleinvariant spectrum at the 6.7σ level (8.7σ when using the highℓ polarization data). The Hubble constant H_{0} decreases by 0.4σ; within the framework of the base ΛCDM model, this increases the tension with some recent direct local determinations of H_{0} (Riess et al. 2016) to around 3.2σ.
The use of SimLow in place of lowP has little effect on most of the usual extentions to the ΛCDM model, as can be seen in Table 9. The number of relativistic species, for example, remains compatible with 3. The phenomenological A_{L} parameter is essentially unchanged, and still discrepant at roughly the 2σ level from its expected value of 1. Additionally, the running of the spectral index is constrained to be even closer to zero.
Due to the lowering of the normalization, the CMB constraint on neutrino masses improves from ∑ m_{ν} < 0.72 eV (PlanckTT+lowP) to ∑ m_{ν} < 0.59 eV (PlanckTT+SimLow) and ∑ m_{ν} < 0.34 eV (PlanckTTTEEE+SimLow). When adding BAO information (see Planck Collaboration XIII 2016, for details), these constraints improve further to m_{ν} < 0.17 eV (PlanckTT+SimLow+lensing+BAO) and m_{ν} = 0.14 eV (PlanckTTTEEE+SimLow+lensing).
Constraints on 1parameter extensions of the base ΛCDM model obtained using the PlanckTT likelihood in combination with lowP and SimLow.
The HFI polarization measurements presented in this paper reduce some of the tensions for some of the cosmological parameters, but have no substantial impact on the qualitative conclusions of the Planck 2015 cosmology papers. In particular, the ΛCDM model remains an excellent description of the current data.
Further interpretation of these results in terms of reionization models can be found in the companion paper Planck Collaboration Int. XLVII (2016). A more complete quantitative description of all the cosmological consequences of the new HFIbased constraints, along with a complete reanalysis of the highℓ polarization using these new data, will be performed in a forthcoming Planck release.
8. Conclusions
This paper presents a value of the parameter τ derived solely from lowmultipole EE polarization, which is both the most accurate to date and the lowest central value. It depends on the highℓ multipoles of the TT power spectrum only through the constraint required to fix A_{s} e^{− 2τ}.
Measurements of polarized CMB anisotropies on these large (>10°) scales require allsky coverage and very low noise. Only space experiments with great stability, multiple redundancies, and enough frequencies to remove Galactic foregrounds, can achieve all of these simultaneously. The WMAP satellite was the first to reach most of these goals, with two telescopes to directly measure very large scales, passively cooled detectors, and nine years of observations; however, its lack of high frequencies did not allow a sufficiently accurate dust foreground subtraction to be carried out. The Planck mission achieved much lower noise, especially with the 100 mK bolometers of the HFI. Scanning the sky in nearly great circles at 1 rpm for 2.5 yr, Planck HFI required extreme stability and many levels of redundancy to measure the polarized CMB on the largest scales. The challenge was such that, for the first two cosmological data and science releases by the Planck team, there were still obvious signs of poorly understood systematic effects, which prevented the use of largescale HFI polarization data.
At 545 and 857 GHz, calibration is not as accurate as in the CMBcalibrated channels (70–353 GHz), but polarized dust emission can be removed well enough that τ is unaffected by dust residuals. Furthermore, most of the synchrotron foreground that is uncorrelated with dust is also mostly removed when using only 100 × 143 crossspectra. At CMB frequencies, polarized systematic effects have now been understood and modelled. One systematic effect is not yet corrected in the maps, but is removed as a small correction to the power spectra using simulations. Finally, the use of the 100 × 143 QML crossspectrum reduces even further the impact of the remaining systematic effects at CMB frequencies. These advances enable a measurement of τ to be made using only the low EE multipoles from HFI, along with the highmultipole TT constraint on A_{s} e^{− 2τ}. The path is now clear to the final Planck release in 2016, which will provide fullsky, polarized, CMB maps that are reliable on all scales.
This τ determination fully reconciles the CMB results with other astrophysical measurements of reionization from sources at high redshift. It also gives constraints on the level of reionization at redshifts beyond that of the most distant sources (z ≳ 10).
Largescale CMB polarization measurements ultimately contain more information on the reionization period than just the value of τ, and can also constrain turbulence in interstellar magnetic fields and Bmode polarization due to primordial gravitational waves, hence the physics of inflation and the early Universe. This last objective requires not only the knowledge of τ that is now achievable, but also data that are fully noiselimited at all frequencies.
Planck (http://www.esa.int/Planck) is a project of the European Space Agency (ESA), with instruments provided by two scientific consortia funded by ESA member states and led by Principal Investigators from France and Italy, telescope reflectors provided through a collaboration between ESA and a scientific consortium led and funded by Denmark, and additional contributions from NASA (USA).
As described in Planck Collaboration VII (2016), the downlink bandwidth allowed one and only one fastsampled detector signal to be sent to the ground at a time. One set of 80 samples was transferred to the ground for any given bolometer every 101.4 s. For all bolometers, the 40 fast samples from each half of the squarewavemodulated signal were summed onboard before being sent to the ground.
A small correlated noise component induced by the deglitching procedure and the TTF corrections leads to nonwhite noise spectra (as in Fig. 1) different from the noise levels predicted from the TOI analysis (Planck Collaboration VII 2016; Planck Collaboration VIII 2016). The differences are significant only at high multipoles.
Acknowledgments
The Planck Collaboration acknowledges the support of: ESA; CNES and CNRS/INSUIN2P3INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MINECO, J.A., and RES (Spain); Tekes, AoF, and CSC (Finland); DLR and MPG (Germany); CSA (Canada); DTU Space (Denmark); SER/SSO (Switzerland); RCN (Norway); SFI (Ireland); FCT/MCTES (Portugal); ERC and PRACE (EU). A description of the Planck Collaboration and a list of its members, indicating which technical or scientific activities they have been involved in, can be found at http://www.cosmos.esa.int/web/planck/planckcollaboration
References
 Bennett, C. L., Halpern, M., Hinshaw, G., et al. 2003, ApJS, 148, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Bertincourt, B., Lagache, G., Martin, P. G., et al. 2016, A&A, 588, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bond, J. R., Jaffe, A. H., & Knox, L. 1998, Phys. Rev. D, 57, 2117 [NASA ADS] [CrossRef] [Google Scholar]
 Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 811, 140 [NASA ADS] [CrossRef] [Google Scholar]
 Catalano, A., Ade, P., Atik, Y., et al. 2014, J. Low Temp. Phys. 176, 773 [Google Scholar]
 Delabrouille, J., Betoule, M., Melin, J.B., et al. 2013, A&A, 553, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dijkstra, M. 2014, PASA, 31, 40 [NASA ADS] [CrossRef] [Google Scholar]
 Efstathiou, G. 2004, MNRAS, 349, 603 [NASA ADS] [CrossRef] [Google Scholar]
 Efstathiou, G. 2006, MNRAS, 370, 343 [NASA ADS] [CrossRef] [Google Scholar]
 Fan, X., Strauss, M. A., Becker, R. H., et al. 2006, AJ, 132, 117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fixsen, D. J. 2009, ApJ, 707, 916 [Google Scholar]
 Fixsen, D. J., Weiland, J. L., Brodd, S., et al. 1997, ApJ, 490, 482 [NASA ADS] [CrossRef] [Google Scholar]
 Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Gunn, J. E., & Peterson, B. A. 1965, ApJ, 142, 1633 [NASA ADS] [CrossRef] [Google Scholar]
 Hamimeche, S., & Lewis, A. 2008, Phys. Rev. D, 77, 103013 [NASA ADS] [CrossRef] [Google Scholar]
 Hinshaw, G., Larson, D., Komatsu, E., et al. 2013, ApJS, 208, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Hivon, E., Górski, K. M., Netterfield, C. B., et al. 2002, ApJ, 567, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Kamionkowski, M., & Knox, L. 2003, Phys. Rev. D, 67, 063001 [NASA ADS] [CrossRef] [Google Scholar]
 Keihänen, E., KurkiSuonio, H., Poutanen, T., Maino, D., & Burigana, C. 2004, A&A, 428, 287 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kelsall, T., Weiland, J. L., Franz, B. A., et al. 1998, ApJ, 508, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Lamarre, J., Puget, J., Ade, P. A. R., et al. 2010, A&A, 520, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mangilli, A., Plaszczynski, S., & Tristram, M. 2015, MNRAS, 453, 3174 [NASA ADS] [CrossRef] [Google Scholar]
 Mashian, N., Oesch, P. A., & Loeb, A. 2016, MNRAS, 455, 2101 [NASA ADS] [CrossRef] [Google Scholar]
 Notari, A., & Quartin, M. 2015, JCAP, 1506, 047 [NASA ADS] [CrossRef] [Google Scholar]
 Pajot, F., Ade, P. A. R., Beney, J., et al. 2010, A&A, 520, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration 2015, The Explanatory Supplement to the Planck 2015 results, http://wiki.cosmos.esa.int/planckpla/index.php/Main_Page (ESA) [Google Scholar]
 Planck Collaboration VI. 2014, A&A, 571, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VII. 2014, A&A, 571, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VIII. 2014, A&A, 571, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IX. 2014, A&A, 571, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration X. 2014, A&A, 571, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XII. 2014, A&A, 571, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVI. 2014, A&A, 571, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration I. 2016, A&A, 594, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration III. 2016, A&A, 594, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IV. 2016, A&A, 594, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration V. 2016, A&A, 594, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VI. 2016, A&A, 594, A6 [Google Scholar]
 Planck Collaboration VII. 2016, A&A, 594, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VIII. 2016, A&A, 594, A8 [Google Scholar]
 Planck Collaboration IX. 2016, A&A, 594, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XI. 2016, A&A, 594, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XII. 2016, A&A, 594, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVI. 2016, A&A, 594, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXX. 2016, A&A, 586, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XLVII. 2016, A&A, 596, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck HFI Core Team. 2011a, A&A, 536, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck HFI Core Team. 2011b, A&A, 536, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Riess, A. G., Macri, L. M., Hoffmann, S. L., et al. 2016, ApJ, 826, 56 [NASA ADS] [CrossRef] [Google Scholar]
 Robertson, B. E., Ellis, R. S., Furlanetto, S. R., & Dunlop, J. S. 2015, ApJ, 802, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Rosset, C., Tristram, M., Ponthieu, N., et al. 2010, A&A, 520, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shull, J. M., & Venkatesan, A. 2008, ApJ, 685, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Szapudi, I., Prunet, S., Pogosyan, D., Szalay, A. S., & Bond, J. R. 2001, ApJ, 548, L115 [NASA ADS] [CrossRef] [Google Scholar]
 Tegmark, M. 1996, MNRAS, 280, 299 [NASA ADS] [CrossRef] [Google Scholar]
 Tristram, M., MacíasPérez, J. F., Renault, C., & Santos, D. 2005, MNRAS, 358, 833 [NASA ADS] [CrossRef] [Google Scholar]
 Tristram, M., Filliard, C., Perdereau, O., et al. 2011, A&A, 534, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zitrin, A., Labbe, I., Belli, S., et al. 2015, ApJ, 810, L12 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: HFI systematic errors
Here we describe tests of known systematic errors in the HFI instrument and their effects on lowmultipole polarization.
Appendix A.1: Glitches
Glitches induced by cosmic ray hits on the bolometers are a major source of systematic errors in HFI, and are described in detail in (Planck Collaboration VI 2014; Planck Collaboration X 2014; Planck Collaboration VII 2016). In the TOI processing software, glitches are detected via their transient nature, and are then masked and corrected.
Additionally, a longtail glitch template, built from stacking many events, is fitted to each detected event, and the tail is subtracted from the data. This procedure reduces the noise power by an order of magnitude at frequencies around 0.05 Hz. Uncertainties in this correction contribute to, but do not dominate, the low frequency noise of the detectors. The size of this contribution has been evaluated at the TOI level using simulations in Planck Collaboration VI (2014) and Planck Collaboration X (2014), where we have shown that it accounts for from a few percent up to 30% of the noise, depending on the detector, at frequencies between 0.016 and 0.1 Hz, and is negligible at higher frequencies.
Glitches remaining below the detection threshold contribute to the white noise at a level of 3 to 6%. Polarizationsensitive bolometers (PSBs) in the same module see many coincident events, so both detected and undetected glitches are a potential contaminant of the polarization measurements.
To check the effect of glitches on the polarization maps and power spectra, we used the endtoend simulations, which incorporate a physically motivated model of the glitches (Planck Collaboration X 2014; Catalano et al. 2014). We simulated the full mission data for four detectors at 143 GHz (the first detset), with and without glitches. The level of glitch residuals was evaluated by computing the power spectra of the map differences from the two sets, after running the processing pipeline (including the glitch removal) in simulations with glitches only, and using the same flags in both cases. The results, shown in Fig. A.1, indicate that, after correction, contamination of the E and Bmode polarization data by glitches accounts for up to 30% of the noise in the lowℓ regime (ℓ < 10), as discussed in Sec. 2.3. At low multipoles, residuals are dominated by glitchtail removal errors, while at high multipoles, residuals are dominated by undetected glitches below the threshold.
Figure A.1 gives predictions (from endtoend simulations) propagated to power spectra of all systematic effects and signal (blue curve), half ring difference of noise and glitches (green curve), and residuals after simulation of glitches and removal of the glitches in the processing pipeline (red curve).
Fig. A.1 Estimated contribution of glitch residuals to autopower spectra for simulations of the 143detset1, with the FEE model in black. The red curve shows the evaluation of the power spectra of the difference of maps of simulations with and without glitchdeglitch removal. Spectra here are averaged over four realizations. The blue curve shows the spectra after processing for the simulated signal, with all systematic effects included. The green curve is the autopower spectra of half ring differences, dominated by noise at high mutlipoles and by glitches at low multipoles. 
Appendix A.2: High energy cosmic ray events
We also detect rarer, higher energy events that affect many bolometers simultaneously. The correlated nature of these events is of particular concern for polarization measurements. Very high energy cosmic rays interact with the material in the HFI focal plane and in the Planck satellite, and create dense, energetic particle showers through a cascade of several hadronic interactions. The density of these particle showers generates glitches in many bolometers simultaneously, and heats up the bolometer environment. The densest showers also induce a small temperature decrease in the bolometer plate for 1 to 2 s, followed by a slow rise of the temperature for 1 to 2 min, (Planck Collaboration X 2014). This is due to release of trapped helium on the metal and charcoal surfaces, creating transient gas conduction between parts of the focal plane at different temperatures (bolometer plate, dilution cooler plate, and 1.4K box). The two detectors of a given PSB always show the same transient thermal behaviourin these events. The slow temperature behaviour mode, common to all bolometers, is removed well by subtracting the dark bolometer baseline, and does not affect the data after masking the initial parts of the events.
After deglitching, the global effect of the high energy events is estimated to be very weak. A small effect has been found in polarization power spectra when masking all events detected by 15 or more coincident glitches in the focal plane (120 000 events).
Appendix A.3: Baseline removal and nonlinearities
A low frequency baseline built from the dark bolometer TOI is subtracted from all other bolometers, removing the common effect of the temperature fluctuation of the bolometer plate. This has no significant residual effect on large angular scales, as shown using an endtoend simulation (see Fig. 20 of Planck Collaboration VII 2016). The effect at low multipoles is a small fraction of the signal (<0.3% for ℓ < 20 and 2.5% at ℓ = 4). This leads to levels an order of magnitude smaller than the EE fiducial model at all multipoles.
The bolometers have a small nonlinear response associated with the relatively large changes in the operating point associated with the slow drift of the temperature of the bolometer plate (with the instantaneous linear gain being given by the slope of the tangent at the operating point); these changes in the temperature of the bolometer plate come from Galactic cosmic ray rate modulation by the inhomogeneities of the solar wind, on timescales of hours to weeks (2 to 10 μK), as well as the longterm solar cycle modulation (80 μK) during the mission. The relative gain variation for a 100 μK temperature change is 2−3 × 10^{3}, and is taken into account as a first estimate of the nonlinear response term on the signal. The secondorder term is two orders of magnitude smaller for a μK temperature change (a few 10^{5}), a correction that is below other effects considered here. Planck HFI Core Team (2011a) predicts a 10^{5} gain change due to the optical loading variations coming from the CMB solar dipole at 143 GHz, the largest CMB signal.
A second small nonlinearity is due to the changes in the temperature of the warm electronics boxes on the satellite. The coefficients of this gain change have been measured in ground tests, and the effect has been propagated to the power spectrum using the endtoend simulations and the measured temperature history of the mission. The results are shown in Fig. A.2. The correction is less than 10^{4}μK^{2} in C_{ℓ} for all multipoles. The uncertainty on this correction is thus negligible.
A third nonlinearity, due to the ADC, will affect the thermal baseline differently for different detectors. This will induce a residual effect at low frequencies which contributes to the uncorrelated 1 /f noise.
Fig. A.2 Simulated contribution of warm readout electronics drifts on the autopower spectra of 143detset1. The red curves show the contribution of the electronics drift compared to the total map power (orange) and the FEE input to the simulations (black). The green curve is an estimate of the noise level based on the halfring differences. 
Appendix A.4: ^{4}HeJT cooler pickup
Electrical coupling between the ^{4}HeJT cooler drive and the bolometer readout (which are phaselocked) appears as very narrow lines (called “4K lines”) in the power spectrum of the TOI at harmonics of 10 Hz. These are removed from the TOI (see, e.g., Planck Collaboration VII 2016). Residuals from this removal process can affect the angular power spectra of the CMB, although due to the nearconstant spin rate of the Planck spacecraft, the 4K lines project to the sky at discrete angular scales. The lowest frequency line, seen at 10 Hz in the TOI, projects to ℓ ≃ 600. Since this feature is well away from the reionization structure in the angular power spectrum that we consider here, we do not attempt further cleaning of these features from the data.
Appendix A.5: Detector crosstalk
Electrical coupling, or crosstalk, between the pair of detectors in a PSB can create errors in the recovered polarization. The groundbased calibration of HFI and the inflight measurements of current crosstalk put an upper limit to this coupling at the level of 10^{3} (Pajot et al. 2010). Planck HFI Core Team (2011a) states this limit as −60 dB. Nevertheless, the square modulation of the bias current of the bolometers was kept at a constant amplitude throughout the mission.
The voltage crosstalk is more relevant than the current crosstalk for the coupling of the signal between detectors, and is estimated inflight using cosmic ray glitches. While the vast majority of the long glitches are coincident between the a and b detectors within a PSB pair, there is also a population of coincident short glitch detections, which have a shift in time on the order of tens of milliseconds that cannot be due to a single particle, and must come from crosstalk. These phaseshifted coincidences are stacked to solve for a crosstalk amplitude, which is typically in the range 1−3 × 10^{3} for pairs of PSB detectors. The crosstalk level is consistent from a to b bolometers and from b to a, in all cases.
Using the endtoend simulations, we propagate the measured voltage crosstalk between the PSBs into maps of the polarized sky using the glitchmeasured crosstalk amplitudes and phase shifts, and compare these to a simulated nocrosstalk case. The main effect of the crosstalk is to bias the relative calibration between detectors. The SRoll mapmaking algorithm recovers the calibration correction due to crosstalk to very good precision (as shown in Table A.1).
Figure A.3 shows the EE and BB angular power spectra of the simulations. Introducing detector crosstalk in the simulation changes the recovered polarization signal by a small amount; the absolute value of the difference is lower than the noise variance in the simulation. This might need to be corrected in future for the detection of very low level polarized signals.
Crosstalk simulation.
Fig. A.3 Autopower spectra EE (left) and BB (right) of the 217 GHz endtoend simulations carried out to test the effects of detector crosstalk. Red crosses indicate the FEE sky signal input to the simulation, to which a realization drawn from the noise variance is added (dashed black). The orange and blue crosses show the total power in the simulations, with and without simulated crosstalk, respectively. The absolute value of the difference (between the simulations including crosstalk and those without crosstalk) is shown in green. We note the difference in binning data between left and right plots. 
Appendix A.6: Instrumental polarization
The polarization efficiencies and angles of HFI detectors were measured on the ground, directly on the instrument (Rosset et al. 2010). The polarization efficiencies were determined to better than 0.2% for all PSBs. This error was shown to induce an uncertainty lower than 0.1% in the Emode power spectrum.
Detector polarization angles were propagated to the sky reference frame using mechanical modelling of the telescope and optical simulations. Errors were decomposed into two parts: a global error of , common to all detectors and due to the uncertainty in the angle of the calibrating polarizer; and an error for each detector, estimated to be lower than , dominated by systematic uncertainties in the ground apparatus. Such error levels, both for global focal plane angle and detector relative angles, were shown to induce an uncertainty below 10% of the cosmic variance on the Emode power spectrum in the range ℓ = 2−1000.
In order to confirm the ground calibration measurements on the sky data, we use the TB and EB spectra to check for a possible global rotation of the focal plane at 100, 143, and 217 GHz, independently. Rotation by an angle α induces E and B mixing, leading to spurious TB and EB signals given by where we have assumed that and equal zero for the CMB.
Furthermore, beam mismatch introduces intensitytopolarization leakage, which may be parametrized at the spectrum level (Planck Collaboration XIII 2016) as (A.3)where ϵ_{ℓ} and β_{ℓ} quantify the leakage into E and Bmodes, respectively. Assuming that differential ellipticity is the main source of beam mismatch, β_{ℓ} and ϵ_{ℓ} are dominated by m = 2 and m = 4 modes (Planck Collaboration XI 2016), and can be written as (A.4)Combining Eqs. (A.1)−(A.4), we can build the five parameter (α, β_{2}, β_{4}, ϵ_{2}, ϵ_{4}) model to be fitted to the TB and EB spectra at 100, 143, and 217 GHz.
To compute the angular power spectra, we subtract dust emission from the Q and U detset maps, using the 353 GHz maps as dust templates. Masking 50% of the sky as well as point sources, we use Xpol (Tristram et al. 2005) to compute the crossspectra between the two detsets in each channel, which efficiently suppresses the noise bias. We sample the model parameter space using a Markov chain Monte Carlo (MCMC) algorithm from which the probability densities of the parameters may be recovered. The angle α is constrained by the MCMC to the following values (68% CL, statistical only), depending on the frequency and spectrum examined.

100 GHz: ; .

143 GHz: ; .

217 GHz: ; .
The beam mismatch parameters are not constrained by the MCMC, but they have little impact (i.e., ≲ 1σ) on the values of α quoted above. Whatever the channel and spectra, the angle is never significantly different from zero, and always smaller than the average errors on PSB orientation (~ 1°) measured on the ground prior to launch (Rosset et al. 2010).
Comparison with the polarization orientation of the Crab nebula was discussed in Sect. 3.5.2 of (Planck Collaboration XXVI 2016) and shown to be consistent with the above numbers.
We conclude that the uncertainties of the instrumental polarization parameters as set up in the data processing pipeline have no significant effects on the lowℓEE and BB crossspectra, and are also in broad agreement with preflight determinations.
Appendix B: The SRoll global solution
This Appendix describes the SRoll method for constructing polarized sky maps from the HFI data. The method is divided into three steps, shown schematically in Fig. B.1. The first step bins the timeordered information into HEALPix pixels per period of stable satellite pointing to reduce the amount of data (Sect. B.1.3). The second step fits systematic effects, 1 /f noise, and calibration using differences between measurements in the same HEALPix pixel (Sects. B.1.6 and B.1.7). The third step cleans the data using the fitted parameters, and projects the signal to make polarized maps (Sect. B.1.8). The quality of the final maps depends on the efficiency with which the nonsky signals can be characterized during the fitting procedure in the second step. After describing the method, this appendix discusses masks in Sect. B.2, presents simulations in Sect. B.3, and tests and characterizes the efficiency with which SRoll can clean the signal in Sect. B.4.
Appendix B.1: Method
Appendix B.1.1: Systematic effects considered
The method corrects for baseline drifts due to 1 /f noise in a similar way as the standard unpolarized destriping method used for previous Planck products (Planck HFI Core Team 2011b; Planck Collaboration VIII 2014, 2016), but improves on that method by solving for the amplitudes of templates of additional systematic effects that are not fully corrected in the standard TOI processing, as well as leakage from temperature to polarization. The systematic effects considered are:
 1.
gain variation due to residual analoguetodigital converter nonlinearity (ADC nonlinearity);
 2.
pickup of Galactic emission and the orbital and solar dipoles in the far sidelobes of the telescope (i.e., the response at angles greater than 5° from the main beam axis);
 3.
temperaturetopolarization leakage due to bandpass mismatch between bolometers, from foreground emission with SEDs differing from that of the CMB, namely dust, freefree, and CO;
 4.
temperaturetopolarization leakage due to calibration mismatch on the CMB dipole;
 5.
largescale residual transfer functions that are not corrected in the deconvolution step of the standard TOI processing.
Fig. B.1 Overview and main functional tasks of SRoll. 
Appendix B.1.2: Data model
All PSB bolometers in a given frequency band are used to solve for the sky map at a given frequency. The data for each bolometer (indexed with b) are divided into pointing periods (indexed with i) and into periods of constant gain (indexed with k). The definitions of these data periods are described in detail in Sect. B.1.5. The data sample indexed with p in the equation below corresponds to a pixel within an HPR (HEALPix ring), not a sample in a time series.
We model the HPR data as (B.1)where:

g_{b,k} is the absolute gain of bolometer b during gain period k.

M_{b,i,p} is the bolometer signal sample data in pixel p of pointing period i.

I_{p}, Q_{p}, and U_{p} represent the common sky maps, excluding the dipole, seen by all bolometers in pixel p.

ρ_{b} is the polar efficiency, fixed at the ground measurement value.

φ_{b,i,p} is the detector angle of polarization compared to the northsouth axis.

V_{b,i,p,h} is the template of the empirical transfer function in bin h. The amplitude of this empirical transfer function is a very small correction (around 10^{3}), so the databuilt maps themselves are used to compute the templates using an FFT in the time domain. The source of the template thus suffers from the same deconvolution errors, leading to a 10^{6} level. On the dipole signal this induces a negligible 3.4 × 10^{3}μK_{CMB} shift. Similarly, using the data for foreground templates is more efficient than using a simulation model. The templates are computed for loworder harmonics of the spin frequency only in a small number of bins: h = 0 corresponds to the first harmonic; h = 1 to harmonics 2 and 3; h = 2 to harmonics 4–7, and h = 3 to harmonics 8–15.

γ_{b,h} is the amplitude of the transfer function template for each harmonic bin h.

C_{b,i,p,l} are the templates of foreground components l. The dust component is fitted at all frequencies, the CO component is fitted at all frequencies except 143 GHz, and the freefree component is fitted at 100 GHz only. The templates for the components, without bandpass mismatch leakage, are taken initially from the Planck 2013 results, and then iterated with the maps from previous iterations. The polarization of the foreground is taken into consideration only for the dust.

L_{b,l} is the amplitude of the bandpass mismatch leakage. For all l in the range 1−n_{comp}, we set .

is the total CMB dipole signal (the sum of solar and orbital dipole signals) from the Planck 2013 results: .

is the total dipole seen through the far sidelobes.

is the Galactic signal seen through the far sidelobes.

O_{b,i} is the 1 /f noise effect modelled as one offset per pointing period. We set since the Planck data provide no information on the monopole.

N_{b,i,p} is white noise.
Appendix B.1.3: HPR computation
The SRoll method is based on HPRs as for the Planck 2015 data release (see Sect. 6.1 of Planck Collaboration VIII 2016). This procedure reduces the amount of data, and in the limit that the spinaxis pointing, noise, and systematic effects are constant during a pointing period, the information loss is negligible.
Each detector’s TOI for each period of stable spinaxis pointing is binned into HEALPix pixels at N_{side} = 2048, compressing the data. Furthermore, all skymap templates fitted by SRoll are converted from TOI to HPR (dipoles, CO maps, dust map, and freefree map). The residual transfer function templates (see Sect. B.4.3) are computed at the timeline level and then projected into HPRs.
Appendix B.1.4: Empirical transfer functions
Differential time response on long timescales remaining in the data from the PSBs can create spurious largescale polarization residuals. Estimates of the time response of the bolometers and electronics, referred to as “very long time constants” (VLTCs), determined using planet crossings and stacked cosmic ray glitches, have been deconvolved from the TOI (Planck Collaboration VII 2014, 2016). This technique is not sensitive to time response on timescales longer than one second.
Residuals of far sidelobe (FSL) contributions after removal in the ringmaking phase, together with VLTCs, induce small changes in the phases and amplitudes of dipoles, as well as of all largescale structures in the maps. To remove these small residuals, the HPRmaker builds transfer function templates composed of harmonics of the spin frequency (with real and imaginary parts). These templates are four bands of low spin frequency harmonics of the bolometer signal, and provide a correction for all the residual transfer function effects, without attributing them to a specific physical mechanism through a physical model.
The imaginary part is easy to fit because it is not degenerate with the sky signal. At 100–217 GHz, a VLTC correction had already been included in the Planck 2015 data release to improve the knowledge of the lowest frequency portion of the time response (Planck Collaboration VII 2016). At 353 GHz, the VLTC correction had not been deconvolved inside the time response function. After correction in SRoll, the residual imaginary part of the empirical transfer function is at the 10^{4} level (see Fig. 11), which is subdominant to other systematic errors.
The real part of the transfer function is degenerate with the signal in a single pointing period. Nevertheless, the neargreatcircle scan of the telescope boresight, together with the spin axis precession, ensures that every point on the sky is crossed in at least two directions separated by more than 14°. The shape of the bolometer response along the scan is fitted in regions of the sky with enough structure. For low frequencies (100–217 GHz), the signaltonoise ratio of the Galactic signal is too low to support this fit, and the real part is then not removed at those frequencies. The relative consistency of the dipole and the first peak calibrations shown in Fig. 28 shows also that the real part of the transfer function has a negligible difference between 100, 143 GHz, and 217 GHz. At 353 GHz, the fit succeeds (see Fig. 10) and might be due to the contribution of the foregrounds to the dipole and very low mutlipoles that has not been accounted for.
Appendix B.1.5: Gain step optimization
SRoll assumes that the gain of each bolometer is constant over multiple rings. N_{gain} timesteps are chosen in such a way that the total dipole gain variance is constant over each period k. If there is significant variation of the effective gain within a defined period, then of course the method cannot account for it. This is a problem only for the very long gain periods required when the dipole or Galactic signal happens to be weak. These periods are, by definition, characterized by the relatively small dynamical range of the signal, and therefore a correspondingly small leakage from temperature to polarization. Furthermore, in these periods the signal samples only a small range in the ADC, so it is unlikely to exhibit much in the way of gain variations from the ADC nonlinearity.
Fig. B.2 Gains solved by SRoll for simulated ADC nonlinearities, using gain periods of uniform length (red) and optimized sampling (green). The blue curve shows the apparent gain variation fitted between a simulated total dipole and the same total dipole distorted by the simulated ADC nonlinearities. 
Figure B.2 shows the effects of the sampling on the variable gains solved by SRoll; gain periods of uniform length increase the instability of the gain determination, while optimized sampling improves stability.
Appendix B.1.6: Variable gain measurement
Solving for gain variability necessarily involves solving a nonlinear least squares equation. SRoll uses an iterative scheme to solve for the gains g_{i} of the bolometer at each ring i. At iteration n we have (B.2)the goal being to fit Δg_{i,n}. We can also remove the FSL effect as (B.3)Using Eqs. (B.2) and (B.3), Eq. (B.1) becomes (B.4)where:

i and p are the indices of the ring i and the HEALPix pixel p (recall that the gain is assumed to be constant over mutiple rings, as described in Sect. B.1.5, and the gain for ring i means the gain for the gain period k that contains the ring i);

g_{i} is the gain of the bolometer in ring i;

M_{i,p} is the measured bolometer signal for the HEALPix pixel p of ring i;

I_{p}, Q_{p}, U_{p} is the common sky map, excluding the dipole, seen by all bolometers in HEALPix pixel p;

, where is the solar dipole for the HEALPix pixel p and is the orbital dipole at ring i for the HEALPix pixel p.
Fig. B.3 Histogram of the η parameter measured on real data. For each gainstable period, the total dipole is fitted on the Galactic signal. 
Figure B.3 shows the histogram of this η parameter. Its dispersion is explained by the variation during the mission of the orientation of the ring with respect to the Galactic plane. Moreover, the amplitude of the dipolar component of the Galaxy increases with frequency. Thus the magnitude of the degeneracy increases, and then the dispersion of the η parameter increases with frequency as well. Then, Eq. (B.2) becomes (B.6)At all frequencies  η  < 1, and g_{i,n} → g_{i} when n → ∞. Thus the degeneracy between I_{p} and impacts only the convergence rate of the algorithm; SRoll always converges to the correct gain estimate when the data model is an accurate description of the sky signal.
Table B.1 shows how the gain convergence works in a simulation without noise for various η values. When η = 0, the gain is solved in one iteration. When η is nonzero, the gain convergence reaches 10^{5} after only four iterations. Thus, SRoll uses four iterations.
These simulations also test the degeneracy with bandpass leakage. They show that the degeneracy between gain determination and bandpass leakage is small. As we will see below, more detailed simulations that contain noise enable us to control more precisely the level of accuracy that SRoll reaches while fitting gain and other systematics.
Relative error on the gain given the number of iterations for different foreground levels.
We find that the η variance increases when smaller gainstable periods are used. Figure B.3 shows that the gainstable period cannot be decreased at 353 GHz, otherwise this would lead to  η  > 1. That is why we define longer periods at 353 GHz than at other frequencies. Masking pixels where the Galactic signal is degenerate with the dipole signal also helps to obtain  η  < 1.
An alternative solution could use a Galactic model to fit gain variations on shorter timescales. This procedure would work if the following requirements were met: first, the Galactic template must not have noise that correlates with the data, such as using a previous map version at the same frequency; second, the Galactic template has to have no residual dipole remaining from a previous calibration mismatch; and thus, SRoll currently solves for gain variations with dipole residual fitting.
Appendix B.1.7: Final minimization
During the last step of the Sroll algorithm, the minimization is the same as in the previous iteration, but the dipole is not removed. The minimization then adjusts the other parameters (offsets, VLTC, and bandpasses) with the gains as they were determined at the previous step.
Appendix B.1.8: Projection onto maps
SRoll projects timeordered information to pixel maps using gains and systematic effects fitted during the previous steps. From Eq. (B.1) we obtain (B.7)Thus, inside each pixel we can define (B.8)A linear system is used to find the values of I_{p}, Q_{p}, and U_{p} that minimize this χ^{2}.
Appendix B.2: Masks
Fig. B.4 Masks used during the mapmaking solution. The fraction of the sky used is f_{sky} = 0.862 at 100 GHz, f_{sky} = 0.856 at 143 GHz, f_{sky} = 0.846 at 217 GHz, and f_{sky} = 0.862 at 353 GHz. 
The solution for the map parameters is done with the sky partly masked. This masking has two goals, firstly, to avoid regions of the sky with timevariable emission, and secondly to avoid regions with a strong signal gradient. Thus, the brightest point sources are removed (e.g., the flux of 3C 273 is strong and it changes during the mission). The Galactic plane is also masked, due to the strong signal gradients there. Nevertheless, a relatively large sky coverage is needed to properly solve for the bandpass mismatch extracted from the Galactic signal. Figure B.4 shows the masks used and gives the sky fractions used in the map solution.
Appendix B.3: Endtoend simulations
Appendix B.3.1: General description
Simulations are used to characterize the performance of Sroll. We start from simulated TOI data sets after TOI processing, including the total CMB dipole (solar and orbital), CMB anisotropies, Galactic and extragalactic foregrounds, and noise. These are processed to simulate all systematics relevant for mapmaking. The SRoll mapmaking is then applied to these TOI to obtain the socalled E2E simulated frequency band I, Q, and U maps and associated TT, EE, and BB power spectra. For testing purposes, one can choose specific subsets of input signals, systematic effects, and processing modules.
Building a simulation data set with 100 noise and 100 CMB realizations was not feasible for this paper. We thus test if the CMB can be simulated separately from the noise and systematics and for this purpose two sets of TOIs are simulated, the first including all the inputs, and the second including all inputs except for CMB anisotropies. Maps are made with SRoll from both sets of TOIs. In the second case, the CMB anisotropies are added at the map level after the mapmaking. The difference between these two tests shows that we can avoid including the loop over many CMB cosmological signals at the TOI level.
Figure B.5 shows this difference at the power spectrum level. When the CMB is simulated at the TOI level, the calibration distortion is applied to the CMB and then solved by Sroll. All differences are much smaller than the signal, showing that the CMB can indeed be added to the residual map, i.e., (B.9)The difference is typically less than one percent of the fiducial EE spectrum for the CMB channels (it is higher at 353 GHz, but this channel is only used to remove the dust from CMB channels after multiplication by factors smaller than 10^{2}).
The TT, EE, and BB power spectra are then computed with Spice (Szapudi et al. 2001) from the output maps, using 50% of the sky if not otherwise specified. Then these power spectra are compared with the fiducial TT, EE, and BB spectra.
This also demonstrates that SRoll does not significantly affect the CMB anisotropies in any way and that the systematic effects can be treated as additive. This is an important result, but not entirely surprising because the dominant polarized systematics are the residual dipole leakage from calibration mismatch and the residual foreground leakage from bandpass mismatch, both of which produce spurious polarization signals that add to the CMB anisotropies. The distortions of the CMB anisotropies themselves are much lower than the noise and therefore negligible.
Fig. B.5 Power spectra of the differences between maps made from TOI including the CMB, and maps made from TOI including only noise and systematics with the same CMB added to the map. In this simulation run, the CMB power spectra are set to FEE. 
Appendix B.3.2: Foregrounds
The foreground model used in the simulations is based on the Planck Sky Model v1.9 (Delabrouille et al. 2013), in which each foreground template is convolved with the groundmeasured bandpass, leading to realistic bandpass leakage levels.
Appendix B.3.3: ADC nonlinearity
We build a simple fourparameter model of the ADC nonlinearity that remains after the correction at the TOI level using measurements performed during the warm phase of the mission. The most significant residual after this correction is the central discontinuity in the ADC. Our model of the residual nonlinearity is constructed as a function of the ADU value. The central discontinuity is a step, but it is smeared out by the noise in the data, so we model it as the derivative of a Gaussian. The left column of Fig. B.6 shows the gain as a function of ADU value for four representative bolometers. The values determined by SRoll from the data are shown as blue points. The red line shows the bestfit model for these gains, and the pink area surrounding the line shows the standard deviation computed from 100 realizations of the uncertainties in the fit. The right column of Fig. B.6 shows the measured and modelled gains as a function of pointing period, demonstrating that the model captures the broad features of the remaining gain variation observed during the mission. This model of the residual nonlinearity is used in the simulations. For each realization, we draw the parameters of the model from the fitted distribution, and the consequent nonlinearity is applied to the simulated data. The goal of the model is to simulate the distortion of the dipole that is present in the data (as discussed later in Sect. B.4.2), but not removed by SRoll. The simulations will be used to estimate the level of these systematic effects and to take them into account in the analysis of the data.
Appendix B.3.4: Low frequency transfer function
Lowfrequency transfer functions are not simulated, but are nevertheless solved by SRoll. The fitted parameters are thus expected to be zero. The residual lowfrequency transfer function patterns after SRoll are then interpreted as degeneracies with other simulated systematic errors. The main aim of the resolution of this transfer function is to model the very long time constants (VLTCs) that were not characterized at the TOI level.
Appendix B.3.5: Noise
The noise is assumed to be Gaussian and stationary throughout the mission (see Sect. 2.3). The noise is simulated by drawing a realization for each pointing period from the measured noise power spectrum of each bolometer and its readout chain. Correlations between bolometers are not included in the model.
Fig. B.6 Evolution of the gain solved by SRoll for four representative bolometers. The left column shows gain as a function of signal level in ADU, with values determined from the data in blue. The red line is the best fit, and the pink area surrounding it is the standard deviation of a fourparameter model of the residual ADC nonlinearity fitted to these data. The right column shows gain as a function of pointing period, and shows that the model is able to reproduce the broad features of the remaining gain variation observed during the mission. 
Appendix B.4: Validation and performance
Appendix B.4.1: Bandpasses
Fig. B.7 Error on the recovered leakage coefficients solved by SRoll, normalized to the average of all detectors in a frequency band, on a representative set of detectors. The error bars (in red) are the statistical distribution from 20 realizations of the simulations. 
Figure B.7 shows the error on the reconstructed leakage coefficients, normalized to the average of all detectors in a frequency band, for the three main foregrounds in the HFI channels (dust, CO, and freefree). Although the dispersion is in some cases larger than the statistical noise, the accuracy is better than one percent, and small enough to induce negligible effects on the gain determination and the sky maps produced by SRoll. The excess variance of the bandpass leakage distribution at 100 and 353 GHz is probably due to the degeneracy between the main leakage sources (CO and dust). This is further confirmed by the extreme stability of the intercalibration of the 100, 143, and 217 GHz detectors (see Fig. 13), which have rms dispersions in the range 0.5−2 × 10^{5}.
Synchrotron emission is weak in the HFI channels (see Sect. 5.2), and so the synchrotron leakage will be even weaker. Estimates using bandpasses from ground measurements show that it is negligible.
Fig. B.8 Residual autopower spectra of the CO bandpass leakage residuals. 
Fig. B.9 Residual autopower spectra of the dust bandpass leakage residuals. 
Fig. B.10 Residual autopower spectra of the freefree bandpass leakage residuals. 
Using the reconstructed leakage coefficients, the corresponding templates, and the pointing, we build maps of the residual bandpass leakage. Figures B.8−B.10 show the power spectra of those maps for CO, dust, and freefree foregrounds, respectively. These are negligible with respect to the fiducial EE spectrum at 100, 143, and 217 GHz. The 353 GHz channel is used only to clean the dust emission from the CMB channels, so the additional leakage induced in these channels will be scaled down by the coefficients used to do the cleaning (factors of 3.5 × 10^{4}, 1.7 × 10^{3}, and 1.6 × 10^{2} at 100, 143, and 217 GHz, respectively).
Appendix B.4.2: ADC nonlinearities
The full ADC nonlinearity induces distortions in the signal at low multipoles (as shown in Sect. 2.5), which are mostly corrected at the TOI level. We investigate here the residuals after the TOIlevel correction, as well as the secondary correction applied to them in the mapmaking, using the model described in Sect. B.3.3.
Fig. B.11 Gain variations solved by SRoll (red line; the light red band shows the variance of 20 noise simulation realizations) compared to the gain variation computed directly from the input ADC model (blue). 
Figure B.11 shows the results of ADC nonlinearity simulations for four representative bolometers. We simulate the nonlinearity using these simulations. The gains reconstructed by SRoll from 20 realizations are shown with a red line (mean) and pink band (standard deviation). These can be compared to the blue line, which shows the gain variation computed directly from the input ADC model. The difference between the blue and red lines averaged over the mission gives the average bias on the gain (shown in each panel as “AVG BIAS”) from 1 to 5 × 10^{5} in the CMB channels, and 4 × 10^{4} at 353 GHz.
Fig. B.12 Residual autopower spectrum of the ADC nonlinearity corrected by the gain variation model. 
The results of the simulations are the residual power spectra of the ADC nonlinearity corrected for the gain variation model, shown in Fig. B.12. The residuals in EE for 100 and 143 GHz are low for ℓ> 30, but close to the level of the fiducial EE spectrum for the reionization peak at ℓ = 2−3. At these frequencies, the apparent gain correction and interfrequency calibration from SRoll have been shown to be very accurate (see Fig. B.11 and Table B.2).
Fig. B.13 Residual polarization maps from simulations at 100 GHz. Maps in the left column are obtained by using a constant gain per ring, while maps in the right column are obtained when the simulations are run using the ADC nonlinearity model. 
Figure B.13 shows the residual polarization maps averaged from six simulations of the 100 GHz channel without noise. The left column shows the residuals obtained when the data are simulated using a linear model with a gain that is constant for each ring, but that varies over the mission. This is the model that SRoll fits to the data. It is able to recover the gains almost exactly, so the residuals are at a very low level. The right column shows a simulation including the complete effect of the residual ADC nonlinearity. The additional residuals seen in this case are produced by the distortion of the dipole signal that SRoll does not correct. The associated power spectra explain the lowℓ power spectrum residuals in Fig. B.12, as is demonstrated below for one specific detector.
Figure B.14 shows maps from the same simulations of the ADC nonlinearity (left column) compared to the data (right column) for the 100 GHz channel. The first row shows the maps made before correcting the ADC nonlinearity with the SRoll variable gain model, which demonstrates that the simulations can effectively reproduce the effect seen in the data. The second row shows maps made after correction with the SRoll variable gain model. The largescale patterns seen in the first row are mostly removed by this correction. The dipole distortion shown in Fig. B.13 is still present, but is partly hidden by the noise. The third row shows the map of the simulation when the ADC nonlinearity is not introduced, and contains only noise. The dipole distortion is not present in this map. This is a clear indication that the only residual remaining after the SRoll mapmaking is small and comparable to the noise only at large scales, and is the dipole distortion. A quantitative estimate of the accuracy of the simulation is made in Sect. 2.5 by comparing the power spectra of these maps. It shows indeed that the dipoledistortion residual is lower than the noise and other residual systematic effects at ℓ< 10, and exceeds those only at ℓ = 2−3.
Fig. B.14 Q maps of data (left column) and a simulation of the ADC nonlinearity effects (right column) at 100 GHz. The first row shows the total effect. The second row shows the improvement brought about by the correction of the linear gain variation. The third row shows the noise and residuals from other systematics. 
We have shown that maps made by SRoll contain a residual ADC nonlinearity effect from the distortion of the dipole, which cannot be corrected by the variable gain model. The distortion of the dipole can nevertheless be predicted using the gainvariationbased model of the nonlinearity shown in Fig. B.6 and described in Sect. B.3.3.
Appendix B.4.3: Residual complex transfer function
The correction applied to the VLTC only accounts for the shifts in the dipole that it produces, which had previously allowed us to use the orbital dipole for calibration. It is clear that there must be some time constants with low amplitudes affecting the data at low multipoles beside the dipole, which lead to additional small inscan shifts. The FSLs lead to a crossscan shift which is cancelled to first order in the sum of odd and evensurvey data, but which leaves an amplitude effect in the crossscan signal. Finally, the socalled baffle components of the FSLs is not correctly described by the firstorder GRASP model used to remove them. This produces an extra inscan shift and change in the amplitude. We therefore fit for empirical transfer functions to account for these small residuals. These transfer functions are fitted for four ranges of harmonics of the spin frequency (1, 2−3, 4−7, and 8−15). The accuracy of the reconstruction of the corresponding amplitudes is estimated by using the simulations.
Fig. B.15 Imaginary part of the complex transfer function residuals, simulated for all bolometers and for three ranges of spinfrequency harmonics. Errors bars are computed from 20 noise simulations. 
Figure B.15 shows the mean and standard deviation of the reconstructed amplitudes from 20 realizations of a simulation where no extra transfer functions were included, so the true value is zero in all cases. The results are consistent with zero within the statistical uncertainties. Therefore we conclude that the degeneracies with the other systematic effects are small and that this correction for the residual transfer function will be accurate at the level of the noise.
Fig. B.16 Simulated power spectra of the transfer function residuals compared to the fiducial CMB signal (black curve). For each bolometer, the amplitude of the empirical transfer function was drawn from the uncertainties shown in Fig. B.15. The residuals are negligible compared to the signal. 
Figure B.16 shows that the residual transfer functions in the maps lead to a negligible residual in the power spectrum for the CMB channels. The 353 GHz channel shows an effect close to the level of the fiducial signal, but this channel is used only to clean the CMB channels, so the residuals added to those channels by this process are scaled down by the corresponding cleaning coefficients (as discussed before for other simulations).
Appendix B.4.4: Noise
Figure B.17 shows power spectra of null tests on detector sets and half missions. The blue lines show the data, and the red lines and error bars show the mean and standard deviation of simulations containing noise alone. There is no evidence for unseen systematics in these null tests.
Appendix B.4.5: Calibration mismatch
Calibration mismatch from 20 simulations, in percent.
Fig. B.17 Power spectra of null tests on detsets (left column) and half missions (right column). The data are shown in blue and the mean and standard deviation from simulations containing noise alone are shown in red. 
Fig. B.18 Autopower spectra induced by the residual calibration mismatch in the maps. The residuals in the CMB channels are very low compared to the FEE signal. 
The difference between the mean of the input variable gain and the mean of the variable gain solved for by SRoll gives an estimate of the calibration mismatch in the final maps. Table B.2 shows the minimum, maximum, and rms of the detector calibration in each frequency band. The very low values of the rms of the 20 realizations (3−4 × 10^{5}) demontrate the excellent recovery of the gain by Sroll.
Figure B.18 shows the propagation to the power spectra of the calibration mismatch. For CMB channels, the levels are more than four orders of magnitude below the FEE spectrum.
Appendix B.5: Summary of systematic effects
Figure 17 shows the power spectra of all systematic effects, together with the EE spectrum of the fiducial CMB model. The noise dominates over all systematic effects at ℓ> 5. Bandpass and calibration mismatch are a factor of 100 below the signal. The effect of the dipole distortions induced by the ADC nonlinearity is comparable to the noise at ell ≤ 5 in the CMB channels, and at the level of the fiducial signal at 217 GHz. At 353 GHz, the ADCinduced apparent gain variations are only a factor of a few below the noise at very low multipoles. The bandpass leakage is only a small factor lower. The noise itself is low compared to the dust signal at high latitude (more than two orders of magnitude). It is therefore far below the fiducial signal when scaled with the appropriate coefficients to clean the 100 and 143 GHz maps. In the CMB channels, the dipole distortion (not removed from the data) is significant and particularly strong at very low multipoles. At 353 GHz, the systematic effects are negligible if this channel is only used to clean dust emission from the channels used for measuring τ.
Appendix C: LFI lowℓ polarization characterization
Here we provide a summary of the Low Frequency Instrument (LFI) systematic effects that are most relevant for the polarization analysis at low multipoles. A detailed discussion of all systematic effects for the LFI 2015 release is given in Planck Collaboration III (2016).
We have developed timedependent models for the known systematic effects in LFI. For each effect, we generate a timeline and project it into the map domain using the same pipeline used for the data, including the Planck scanning strategy and LFI mapmaking process. We call the resulting maps “systematics templates” and quantify the impact of these systematic effects by comparing their power spectra to the FFP8 simulations (which reproduce the measured LFI noise spectra, Planck Collaboration XII 2016), and to the recovered sky signal.
Figure 23 summarizes the results of our analysis for ℓ< 45 (resolution approximately 4°). For EE polarization, the most relevant effects come from uncertainties in the radiometer gains. ADC nonlinearity and sidelobe residuals (at 30 GHz) contribute at a lower level. We discuss each of these three sources of systematic effects in turn.
Appendix C.1: Gain reconstruction
The gain of a balanced differential coherent receiver may vary at the 1% level on a variety of timescales, depending primarily on changes in the thermal environment or the bias. For each LFI radiometer, we recover the gain changes as a function of time using the solar dipole (Planck Collaboration V 2016). Imperfect recovery of the true radiometer gain variation produces a systematic effect. The gain solutions are least accurate when the dipole signal is near its minimum, typically a period of several days, which produces effects at large angular scales. At multipoles ℓ< 30 for Stokes Q and U, the rms effect is 0.4 μK at 30 GHz, 0.3 μK at 44 GHz, and 0.15 μK at 70 GHz. These are the dominant systematic errors in the present analysis.
Appendix C.2: ADC correction
Deviations from linearity in the ADC response impact the LFI differently than the HFI. Linearity in the ADC requires that the voltage step size between successive binary outputs is constant over the entire signal dynamic range. Deviations from this ideal behaviour affect the calibration solution. Because the effect induces a variation of the detector white noise that is not matched by a corresponding variation in the absolute voltage level (which would be expected from LFI radiometers), the anomaly can be measured and corrected. We evaluate the residual error through simulations. Since the radiometer gain varies coherently over periods of days to weeks, the average output of the radiometer (and therefore the part of the ADC range being used) also varies on this timescale. Imperfect removal of ADC effects thus produces residuals at low multipoles. Additionally since the effect is independent for each ADC channel, it is not mitigated by differencing orthogonally polarized radiometers. As a consequence, the amplitudes of the residuals are comparable for temperature and polarization, implying of course a larger scientific impact on polarization. For the Q and U components at large scales, the rms residual effect is <0.1 μK at 30 and 44 GHz, and approximately 0.15 μK at 70 GHz.
Appendix C.3: Impact of sidelobes
Telescope sidelobe response can be relevant at low multipoles, both directly and through the calibration process; this is particularly the case at 30 GHz. To deal with this, we remove a model of the sidelobe signatures (“stray light”) in the LFI timelines by combining the GRASP models of the beams with the measured sky maps at each frequency (Planck Collaboration IV 2016). To test the impact of imperfectly modelled sidelobes on the calibration, we ran our calibration simulation pipeline in a mode where the stray light signal was added in the input, but not removed. This provides a worstcase limit (100% error) on the impact of the stray light on calibration. As expected, at 44 and 70 GHz the effects are small, while at 30 GHz we find a significant effect of approximately 1 μK. While this contamination is large compared to the expected CMB EE signal, it is clearly unimportant compared to the synchrotron component that the 30 GHz channel is used to measure, resulting in only a small effect in the 70 GHz map.
Appendix D: SimBaL
Appendix D.1: Description of the method
Within the ΛCDM paradigm, the statistics of low multipole polarization anisotropies depend principally on A_{s}, r, and τ. The highℓ analysis of the TT power spectrum yields a highly accurate constraint on A_{s}e^{− 2τ} (0.5% ), and one may assume negligible primordial tensors. The degeneracy between these two parameters can then be broken by an accurate measurement of the EElowℓ polarization feature since this depends roughly on A_{s}τ^{2}, which is linearly independent of the highℓ constraint. We thus develop a “simulationbased likelihood” code, named SimBaL, which is focused on τ estimation, and based on EE (pseudo) crosspower spectra , using either PCL or QML estimators.
Fig. D.1 Distribution of the SimBaL estimator from signalfree HPFS1 and HFPS1+HPFS3 simulations. The two distributions (blue and green dashed lines) are consistent with a Gaussian width σ = 0.015 for the relevant numbers of simulations. The value from the data (in red) is incompatible with these simulations at the 3.5σ level. 
Given our current understanding of systematic errors, in SimBaL we typically use a single QML crossspectrum constructed from a pair of maps. The Planck maps, at 100 and 143 GHz, are each ILCcleaned using 30 and 353 GHz maps as foreground tracers, as described in Sect. 5.2. The low dimensionality of the problem renders the likelihood analysis amenable to a scheme based on computing , where Ω stands for the other cosmological and noncosmological parameters. We formally have some freedom in the choice of Ω, but typically use bestfit values from the Planck cosmological parameter analysis (Planck Collaboration XIII 2016). More explicitly, we follow these steps.

1.
Compute the power spectra of simulations for various τ values, including noise and systematic effects.

2.
Compute , fitting ℓbyℓ a model to the simulated power spectra. This model is needed in order to reduce noise and to extrapolate beyond the measured range, given by the simulations (see Sect. D.4 for more details). Thus SimBaL (as is true for all likelihoods) needs an accurate model of the distribution at a given multipole. The SimBaL model uses:

, where A is a 3rdorder polynomial, for the central part of the fit;

, where B is a 1storder polynomial, for both tails of the fit.


3.
Compute the estimator, which it is closely related to τ, but easy to compute on pseudoC_{ℓ} spectra: (D.1)

4.
Use simulations to determine .

5.
Finally deduce from .
The simulation set is built from HFPS1 (or HFPS2) for the noise and systematic effects, to which we add 100 CMB realizations. We then use the simulations described above to construct the sampling distributions , for a grid of τ values from 0.010 to 0.180, in steps of 0.001. The relatively low number of noise simulations leads to some apparent features in the distributions, particularly at low τ where the CMB contribution is smaller; however, for a given number of simulations, this would affect the result only of the lowest τ values. Considering the number of simulations we have, this does not affect our results for the range of τ values considered in Sect. 6.
Figure D.1 shows the found with the data (red vertical line) compared to the distribution computed on simulations with noise and systematic effects, but without signal (green and blue curves). Despite the limited number of simulations available, the red vertical line is clearly well outside the histogram of signalfree simulations; this excludes τ = 0 with a significant probability. In order to quantify the significance of this result, we perform a Gaussian fit to the histograms (dashed lines), showing that the value from the data is incompatible with these simulations at approximately the 3.5σ level.
Fig. D.2 Distribution of the τ probability density as a function of the estimator (defined in the text) derived from the 100 × 143 crossspectra. The top panel uses simulations containing only white noise, while the bottom panel uses simulations with the systematic effects residuals added. As an example of a cut through this twodimensional plot, in the top panel an observed value of (horizontal black line) gives a probability distribution (red curve) that peaks at τ = 0.060. In the bottom panel, the same peak value of the probability distribution is obtained with an observed . 
Appendix D.2: Illustration of the method
Figure D.2 shows an illustration of the method for the 100 × 143 crossspectra. The top panel shows, for white noise only, the joint probability distribution in the plane, assuming a uniform prior on τ. For a given value of , the red curve shows the probability distribution of τ, which peaks at a value close to the value of by construction of the definition of . The bottom panel is built using the HFPS1 simulations (with systematic effects) and 100 CMB simulations. In this case, for low values of τ, the maximum of the probability distribution shifts slightly from the diagonal because of the inclusion of the systematic effects. In the top panel the value gives a probability distribution that peaks at τ = 0.06. When using the simulations that include systematic effects (bottom panel), the value gives a peak value of τ = 0.06. When a data power spectrum gives this value of , we determine with simulations that not including the systematic effects would produce a positive bias of about 0.008 in τ. When we include the systematic effects in the simulations, we show our method is able to recover unbiased estimates of τ with an uncertainty small compare to the noise.
Fig. D.3 Pixelpixel matrix χ^{2} values of the first 400 of the 143 GHz FFP8 simulations, evaluated against 18arbitraryeigenmode fits to either the first 50, 83, or 200 simulations. The χ^{2} of the simulation data set used to compute the pixelpixel covariance matrix has been normalized here. Due to the limited number of simulations, the real variance, estimated from the simulations that were not used to compute the pixelpixel matrix, is significantly higher (i.e., the red line jumps at about 50, the green line at about 83, and the blue line at about 200). 
Appendix D.3: Noise variance
Figure D.3 shows the χ^{2} values of low resolution maps from FFP8 noise simulations (see Planck Collaboration XII 2016) evaluated with covariance matrices built with an increasing fraction of the simulations. This illustrates the difficulty introduced by having only a limited number of simulations with which to estimate a lowresolution noise covariance matrix and then to evaluate properties of statistical distributions involving this covariance matrix. Each covariance matrix is computed from a subset of the available simulations, rescaling a fiducial covariance matrix and then adding extra modes to capture the additional systematic effects. It can be seen that the simulations that are not used to build the covariance matrix have a higher χ^{2} than those used to build the matrix. Hence using the same simulations twice gives a misleadingly low indication of the scatter (note that increasing the number of simulations used to fit the noise matrix decreases the χ^{2} discrepancy, as one would expect). Thus, to take into account this feature in the likelihood, SimBaL uses two different simulation sets, one to compute the pixelpixel matrix and the other to measure the noise variance.
Appendix D.4: Dependence on the model
SimBaL uses simulations of noise and systematic effects, to which are added 100 CMB realizations, in order to model the statistical distribution of for each C_{ℓ}. We fit this distribution for each multipole using an asymmetrical polynomial function detailed in Sect. D.1. Figure D.4 shows how the SimBaL polynomial approximation is efficient in capturing the shape of the probability distribution for low τ values when systematic noise effects dominate, while a SimBaL_HL approximation (based on a Hamimeche & Lewis model, Hamimeche & Lewis 2008) is a poor fit to the probability distribution tails for low values of τ. Figure D.5 shows the effect of these two models on the τ posterior computation. SimBaL_HL modelling overestimates the probability of small τ values, where the distribution is dominated by noise and systematic effects.
Fig. D.4 Probability slice at ℓ = 4 for two different τ values, τ = 0.025 (top panels) and τ = 0.06 (bottom panels). The left panels show the signalonly distribution, while the right ones show the signal plus noise and systematic effects. For small τ values, the distribution is dominated by noise and systematic effects and is fitted poorly by the HamimecheLewistype model (red line), but is wellcaptured by the SimBaL polynomial approximation (blue curve). 
Appendix D.5: Dependence on masks
Figure D.6 shows four sky masks with f_{sky} ranging from 0.30 to 0.60, built by thresholding the intensity of the Galactic polarized emission. Figure D.7 shows the effects of sky masking on PCL τ posteriors for the data. This shows that the effect of varying the retained sky fraction is small. In the main likelihood analysis, the mask used is the one with f_{sky} = 0.50, which is in fact the sky fraction that yields the best accuracy (smaller posterior width). Even for the more extreme sky fraction used there is a small bias (of around 0.4σ) towards lower peak values for the τ probability distribution. The generally very good consistency between these curves around the used f_{sky} = 0.50 shows that the componentseparation process is not a limitation in the determination of τ in this paper.
Fig. D.5 Posterior distributions of τ from SimBaL (blue curve) and SimBaL_HL (red curve), showing how they differ at low τ. The red line is an attempt to mimic in SimBaL the modified “HL” approach, and demonstrates that the divergence is due to residual nonGaussian behaviour of the statistics at low multipoles (see text). 
Fig. D.6 Masks for different values of f_{sky}. The mask shown in the bottomleft panel (50% sky fraction) is the one used in the Lollipop and the SimBaL analyses. The other three were used to evaluate the dependence of the results on changing sky fractions from 30 to 60%. 
Fig. D.7 Posterior distributions from SimBaL PCL 100 × 143 crosspower spectra obtained with different sky fractions, showing that variations caused by masking the Galaxy have little effect on the τ posterior likelihood for a range of mask sizes. 
Fig. D.8 Posteriors on τ using different multipole ranges in the 100 × 143 PCL crosspower spectrum. The ℓ ranges are indicated in the legend. 
Fig. D.9 Posteriors in τ when removing one multipole at a time in the 100 × 143 QML crosspower spectra. The discrepant curve is when removing ℓ = 5. 
Fig. D.10 Histograms of the peak values of the τ posterior distributions after removing one multipole each time, computed from 8300 simulations with τ = 0.055. For the red curve the maximum value of the 100 × 143 QML power spectrum in the range ℓ = 2−12 is removed, and for the blue curve ℓ = 5 is removed when it happens to be the maximum. The peak value of the τ posterior when removing ℓ = 5 from the data is shown in green, and is consistent with the blue curve, as expected. 
Fig. D.11 computed for the data using SimLow. Power spectrum values computed for different τ are shown in blue. For ℓ = 3 and ℓ = 5 we can see that very small values of τ are excluded, while for ℓ = 4 large values of τ have a low probability. 
Appendix D.6: Dependence on multipole range
Figure D.8 shows the effect of changing the range of multipoles used in the 100 × 143 PCL crosspower spectrum. As expected, using a minimum ℓ of 5 or 6 starts to bias the result because of the low level of the reionization feature above these multipoles. Although we know that the effects of ADC nonlinearity on the dipole are concentrated at ℓ = 2 and 3 before removal, including these multipoles does not affect the result substantially. This is an important result as it confirms that this systematic effect residuals in the cross power spectra are small as discussed in Sect. 5.3.1 and accurately simulated in the HFPS, and accounted for when used in the SimBaL likelihood.
Figure D.9 shows the effect of removing one multipole at a time from the QML results. Except for the case of removing ℓ = 5, which is the multipole where the EE crosspower spectrum is maximum and the probability to exceed near 2 σ (see Fig. 33), we obtain very stable results. Figure D.10 shows how removing one multipole from the τ posterior computation where the power spectrum is maximum does bias the result low as expected; however, the τ value obtained when removing ℓ = 5 is statistically consistent with what is expected from simulations.
As a consequence of this discussion, the likelihood analysis in this work uses the multipole range 2−20.
Appendix D.7: SimLow: a lowℓ likelihood based on SimBaL
SimBaL provides a τ posterior using the other parameters based on the best cosmology from the TT, TE, and EE power spectra at higher multipoles. Thus, SimBaL is not usable as a lowℓ likelihood for cosmological parameters that are not the ones used to build SimBaL. Because of this, we have developed SimLow to estimate a set of based on the SimBaL paradigm, by computing using simulations. The QML evaluation of the power spectra minimizes correlations between different ℓs, and SimLow computes the likelihood independently between multipoles.
In order to measure for the 100 × 143 power spectrum, we compute 100 realizations of each , with 3000 steps of 10^{4}μK^{2} added to the HFPS1 and HFPS3 simulations using the QML2 estimator for the power spectra. Figure D.11 shows the results of running SimLow for the first few multipoles of the data.
SimLow uses posteriors for multipoles in the range ℓ = 2–20. Using SimLow to compute the τ posterior with the same cosmological parameters as SimBaL (Fig. D.12) gives τ = 0.055 ± 0.009, which is fully consistent with the value from SimBaL and referred to as lowEH.
Fig. D.12 Posterior distribution for τ computed with the SimLow likelihood using the same cosmological parameters as for SimBaL. The posterior is consistent with the LowEH result. 
Appendix E: Glossary
This Appendix gathers definitions of acronyms and other terms widely used in this paper, in addition to those more general terms listed in the glossary and acronym list of Planck Collaboration (2015).
ADC nonlinearity:
analoguetodigital converter nonlinearities. The HFI bolometer readout electronics includes a 16bit ADC that has a very loose tolerance on the differential nonlinearity (the maximum deviation from one least significant bit, LSB, between two consecutive levels, over the whole range), specified to be not worse than one LSB. The implications of this feature for HFI performance proved to be a major systematic effect on the flight data. A wide dynamic range at the ADC input was needed to both measure the CMB sky and foregrounds, and properly characterize and remove the tails of glitches from cosmic rays. Operating HFI electronics with the necessary low gains increased the effects of the ADC scale errors on the CMB data (see Sect. 2 of Planck Collaboration VIII 2016, for further details).
ADU:
analoguetodigital unit.
Complex transfer function:
an empirical function that captures residuals of the bolometer/electronics time response deconvolution (including VLTC) as well as residuals from far sidelobe effects.
Detset or ds:
“detector set”, i.e., a combination of sets of polarizationsensitive bolometer pairs with both orientations. Specifically, 100 GHz ds1 combines 1001a/b and 1004a/b, 100 GHz ds2 combines 1002a/b and 1003a/b, 143 GHz ds1 combines 1431a/b and 1433a/b, and 143 GHz ds2 combines 1432a/b and 1434a/b.
Distorted dipole:
the difference between the actual dipole signal (which is affected by the ADC nonlinearity, like all signals) and the sine wave that would have been measured without the nonlinearity.
FTT, FEE, and FBB models:
CMB fiducial power spectra, based on bestfit Planck cosmological parameters (Planck Collaboration XIII 2016), with τ = 0.066 and r = 0.11.
FSL:
far sidelobe, i.e., the pickup of the sky signal, dominated by the spillover of radiation around the edges of the secondary and the primary telescope mirrors after being reflected by the secondary mirrors and main baffle (with the usual convention following the light from the detectors outwards).
HFPS:
HFI focal plane simulations, built with the pre2016 E2E software pipeline. They contain realizations of the noise described in Sect. B.3.5, with the systematics dominated by the additional ADC nonlinearity model described in Sect. B.4.2. HFPS1 contains 83 realizations, while HFPS2 and HFPS3 each contain 100 realizations.
HPR:
HEALPix ring, i.e., a partial map produced from projection onto a HEALPixpixelized sky data from a single pointing period.
PCL:
PseudoC_{ℓ} estimator, a specific method for deriving a crossspectrum.
Pre2016 E2E:
endtoend simulation pipeline, based on the 2015 E2E simulation pipeline, where the ADC nonlinearity correction in the mapmaking process has been added to the TOI. The 4K line and convolution effects are not simulated because that could affect the higher multipoles.
QML:
quadratic maximum likelihood estimator method, a specific method used to derive a crossspectra. The pixelpixel covariance matrices (QML1 and QML2) are built using three different QML estimators. For QML1 (or QML2), the pixelpixel matrix has been computed directly with the HFPS1 (or HFPS2) data set. For QML3, in order to minimize the effect of overfitting, the offdiagonal terms are reduced to the first four eigenmodes of a principal component analysis from HPFS1.
SimBaL1,2,3:
three posterior distributions for τ. They differ in the choice of the power spectrum estimator and the data set simulation used. SimBaL1 (or SimBaL2) uses QML1 (or QML2) as power spectrum estimator and HFPS2 (or QML1) as the simulation data set. SimBaL3 uses QML3 as power spectrum estimator and the full set of simulations (HFPS1+HFPS2+HFPS3) as the simulation data set.
SRoll:
a new polarized destriping algorithm, which removes residuals via template fitting, using HPR to compress the timeordered information.
VLTC:
very long time constant, i.e., bolometer time response longer than 1 s, which even at low amplitude may bias the calibration by distorting the dipole signal and causing leakage of the solar dipole into the orbital dipole.
All Tables
Estimates of the relative impact of the FSLs per bolometer within a frequency band.
Solar dipole amplitude and direction by frequency, as well as for the average (AVG) of the 100 and 143GHz maps.
Relative calibration between the solar dipole and the first and second peaks (ℓ = 110–500) in the CMB power spectrum with respect to 100 GHz.
For each CMB frequency, the amplitude of the power spectrum D_{ℓ} at ℓ = 4 that is removed for the dust and synchrotron foregrounds.
⟨ D_{ℓ} ⟩ _{rms} over 2 ≤ ℓ ≤ 8 for autospectra simulations of the ADCinduced dipole distortion and noise at 100 and 143 GHz (Fig. 17), together with simulations of 100 × 143EE crossspectra (Fig. 31) and data BB crossspectra (bottom panel of Fig. 33) tracing final residuals and noise (signal is negligible).
Peak values, 68% upper and lower limits, together with 95% upper limits, for the two likelihood methods (SimBaL and Lollipop) and the two crossspectra estimators (PCL and QML).
Parameter constraints for the base ΛCDM cosmology (as defined in Planck Collaboration XVI 2014), illustrating the impact of replacing the LFIbased lowP likelihood (used in the 2015 Planck papers) with the HFIbased SimLow likelihood discussed in the text.
Constraints on 1parameter extensions of the base ΛCDM model obtained using the PlanckTT likelihood in combination with lowP and SimLow.
Relative error on the gain given the number of iterations for different foreground levels.
All Figures
Fig. 1 Mean power spectra of the signalsubtracted, timeordered data from Survey 2 for each polarizationsensitive HFI frequency channel. The spectra are normalized at 0.25 Hz. Blue, green, red, and cyan represent 100, 143, 217, and 353 GHz, respectively. The vertical dashed line marks the spacecraft spin frequency. The sharp spikes at high frequencies are the socalled 4K cooler lines. These noise spectra are built before the time transfer function deconvolution. 

In the text 
Fig. 2 Noise crosspower spectra of the 143GHz bolometers, with the unpolarized spiderweb bolometers (SWBs) in red and the polarizationsensitive bolometers (PSBs) in blue. The lowlevel correlated white noise component of the PSB noise is associated with common glitches below the detection threshold. Autospectra are shown in black. The uncorrelated noise is in green. 

In the text 
Fig. 3 Histogram of the noise between 0.018 and 0.062 Hz (frequencies at which it is dominated by the uncorrelated 1 /f noise) for detector 1431a in blue, together with the bestfit Gaussian distribution in red. 

In the text 
Fig. 4 Autopower spectra, showing the level of the simulated FSL projected on the maps predicted using the GRASP model. At ℓ < 10, the FSL signal at all frequencies is at least two orders of magnitude smaller than the cosmological FEE signal. 

In the text 
Fig. 5 Relationship between input and output, for one spare ADC. The plot shows the difference between the measured digitized output signal level and the one with a perfectly linear ADC, as a function of the output level, over a signal range appropriate for the sky signals. Thus, on the one hand, a perfectly linear device would be a horizontal line at zero; in a real device such as shown here, on the other hand, the relationship between input and output is complicated and nonlinear everywhere, especially near the middle of the range around 0 ADU. 

In the text 
Fig. 6 Fastsampled signal from bolometer 1431a. Left: eighty samples from early in the mission, corresponding to one cycle of the 90Hz squarewave modulation of the bias current across the bolometer. Because of the squarewave modulation, the first and last sets of 40 samples are nearly mirror images of one another across the xaxis. The units are raw analoguetodigital units (ADU). Middle: normalized histogram of the fastsample signal values for day 91 of the mission, in ADU. The signal spread is dominated by the combination of the noise and the squarewave modulation of the bias current across the bolometer (left panel), with additional contributions from the CMB solar dipole. Right: normalized histograms for each day, starting with day 91, stacked left to right for the entire mission. Histogram values are given by colour, as indicated. The obvious symmetric trends during the mission are caused by drifts in temperature of the bolometer plate. Isolated days with large deviations from zero, seen as narrow black vertical lines, are due to solar flares. 

In the text 
Fig. 7 Simulation of the modulated noiseless sky signal per ring shows the variable amplitude of the dipole and thermal drift, which dominates the signal, as a function of time (expressed in stable pointing period number) during the entire mission for four representative bolometers. The different signs of the drifts derivative depends on the level of the compensation of the modulation, which can bring one state of the modulation on either side of the middle of the ADC. 

In the text 
Fig. 8 Dipole amplitude difference on the rings observed twice, one year apart, for each of the 143GHz polarizationsensitive bolometers. This detects the timedependent response associated with the excursions of the signal on the ADC, illustrated in Fig. 7. The blue curve shows the dipole differences in units of μK after ADC correction in the TOI processing, with no other processing. The red curve shows the measured dipole amplitude solved by Sroll, demonstrating the reliability of the model, which can then be applied to small signals. 

In the text 
Fig. 9 Power spectra of the 100GHz halfmission nulltest maps shown in Fig. B.14, which are dominated by ADC nonlinearity effects. We compare six simulations at 100 GHz drawn from the Sroll uncertainties (red lines), and the average for the full frequency 100GHz data (blue line). In the top panel only the TOI nonlinearity correction is applied. In the bottom panel, both the TOI correction and the timevarying gain correction are applied. This leaves only the dipole distortion. For both panels the green lines are for simulations without the ADC nonlinearity effect (only noise and other systematic residuals). 

In the text 
Fig. 10 Bestfit solutions for the real and imaginary parts of the empirical additional transfer function as a function of frequency, for the 353GHz bolometers. 

In the text 
Fig. 11 Ratio of the fitted data to simulated patterns detecting the residual imaginary part of the empirical transfer function, measured in odd minus even Survey difference maps averaged for sets of harmonics. The transfer function correction has been applied only over the four first sets of harmonic ranges (ℓ < 15); higher harmonics have not been corrected by the empirical transfer function. 

In the text 
Fig. 12 Residual solar dipole amplitude for each bolometer, by Survey. The average dipole at each frequency is subtracted. For 100 and 143 GHz (top panels), the variations are compatible with the relative calibration uncertainty of 10^{4}. At 353 GHz, the scale is expanded by a factor of five, and all detectors show an obvious odd/even pattern, which is marginally apparent at 217 GHz. 

In the text 
Fig. 13 Relative calibration measured by the dipole amplitude for each polarizationsensitive detector with respect to the mean dipole across a frequency band. 

In the text 
Fig. 14 Comparison of the response of each detector to dust as measured on the ground (blue) and as solved by SRoll (red). 

In the text 
Fig. 15 Simulation of the templatefitting tests for temperaturetopolarization leakage in the 100GHz maps, with (blue) and without (red) ADCnonlinearity dipole distortion. The green curve shows the result of the leakage fit in the HFI pre2016 data, which is comparable to the simulation with the ADC nonlinearities. The red curve is lower by a factor of 5 or more, showing that dipole distortions due to ADC nonlinearity are a significant contributor to the leakage. The black dashed line corresponds to the fiducial model FEE with τ = 0.066. 

In the text 
Fig. 16 EE auto and crossspectra of the global fit test of the temperaturetopolarization leakage, for 100, 143, and 217 GHz. Levels are similar to that for 100 GHz, for which Fig. 15 shows that this is dominated by the ADCnonlinearity dipole distortion. The dashed black line corresponds to the fiducial model FEE with τ = 0.066. 

In the text 
Fig. 17 Residual EE autopower spectra of systematic effects from the HFI pre2016 E2E simulations computed on 50% of the sky (colours specified in the top left panel apply to all panels). The purple line (ADC NL total residual) shows the sum of all effects associated with ADC nonlinearity. The dark blue line (ADC NL no distortion) shows the level without the dominant dipole distortion. The plots show also the FEE model (black curves). The 100GHz and 143GHz model scaled to 353 GHz with a dust SED is shown as dashed and dotted lines, respectively. 

In the text 
Fig. 18 EE autopowerspectra of detectorset, halfmission, and halfring difference maps for 100, 143, and 217 GHz. C_{ℓ} rather than D_{ℓ} is plotted here to emphasize low multipoles. Results from the Planck 2015 data release are on the left; results from the pre2016 maps described in this paper are on the right. Colourcoding is the same for all frequencies. We also show for reference an average of FFP8 simulations boosted by 20% to fit the halfring null tests. The black curves show the FEE model. 

In the text 
Fig. 19 EE and BB crosspower spectra of the residual effect computed from null tests between halfmission maps based on full mission minimization (red curve) or independent minimization for each half mission (blue curve). This second approach clearly shows a systematic effect. The sum of all systematic effects, dominated by the ADCnonlinearity dipole distortion shown in green in Fig. 17, is at the same level as the simulated null test with independent minimization. The FFP8 power spectrum is given again for reference. 

In the text 
Fig. 20 EE and BB spectra of the 2015 maps and the pre2016 maps used in this work at 100, 143, 217, and 353 GHz. The crossspectra of detset and halfmission maps and the autospectra of the detset and halfmission difference maps are shown. The maps are masked so that 43% of the sky is used. 

In the text 
Fig. 21 Schematic of the simulation plan to characterize the LFI polarization data at low multipoles. Both the instrumentbased and nulltestbased strategies are represented (see text). The lower left part of the diagram outlines the crossspectrum simulation analysis involving LFI and HFI data, including systematic effects (see Sect. 6). 

In the text 
Fig. 22 Bias on τ due to systematic effects, specifically showing the distribution of τ with (red/dashed) and without (blue/solid) the systematics template. The vertical dotted line shows the input value to the simulation, τ = 0.065. 

In the text 
Fig. 23 Power spectra (D_{ℓ} in μK^{2}) of systematic effects at 30 GHz (left) and 70 GHz (right), with each effect coded as indicated in the legend. 

In the text 
Fig. 24 Power spectra for 70 GHz of systematic effects from our instrumentbased systematics simulations (with no noise) for: the full 4year mission in green; the map difference (yr1 + yr2)−(yr3 + yr4) in red; and the map difference (yr2 + yr4)−(yr1 + yr3) in blue. 

In the text 
Fig. 25 Power spectrum of the measured 70GHz null maps for the differences D_{12−34} (orange) and D_{24−13} (light blue). The scatter is dominated by the measurement noise. The error bars are computed from 1000 Monte Carlo simulations from FFP8 noise. Negative values are due to the fact that the noise bias has been removed. For comparison we also plot the systematicsonly spectra (red and blue). 

In the text 
Fig. 26 Bias on τ due for the nulltestbased systematics templates. The blue line shows the distribution of τ for the reference CMB+noise simulations, the red line shows the resulting distribution when including also the null systematic maps built from the instrumentbased simulations. The vertical magenta line shows the value of τ estimated for the single CMB realization to which the null maps from actual data, as opposed to the null maps from the instrumental simulations, were added. The dotted line shows the input value to the simulation, τ = 0.065. Top: D_{12−34} null combination. Bottom: D_{24−13} null combination. 

In the text 
Fig. 27 EE and BB differences and crossspectra of the 2015 maps at 30, 44, and 70 GHz, for different data splits. These maps are masked so that 43% of the sky is used. This figure show similar plots to Fig. 20 for LFI frequencies. 

In the text 
Fig. 28 Relative calibration based on measurement of the first two acoustic peaks (blue) and on the solar dipole (red). The bottom panel is a zoom of the top panel. We use 100 GHz as the reference frequency for the dipole calibration method. In this plot, Planck 30GHz and WMAP Qband data may be affected by their relatively low angular resolution. 

In the text 
Fig. 29 Dust correction to C_{ℓ} in μK^{2}, using the 353GHz channel as a template (in blue), along with the synchrotron correction to C_{ℓ} in μK^{2} using the 30GHz channel as a template (in red). These are plotted for the 70–217 GHz channels as a function of frequency. We note that the points for synchrotron at 143 and 217 GHz are very low and can only be considered as upper limits. 

In the text 
Fig. 30 Residual errors in EE from component separation, estimated from the scatter of the componentseparation coefficients. The fiducial EE spectrum and the noise plus systematics residuals (green line) are also shown. 

In the text 
Fig. 31 Similarly to Fig. 17, residual EE crosspower spectra of systematic effects from the HFI pre2016 E2E simulations computed on 50% of the sky are shown for 100 × 143. Dotted lines are used to join two multipoles when one is negative. 

In the text 
Fig. 32 100 × 143 crossspectra of noise and systematics from the average of HFPS1, HFPS2, and HFPS3 simulations, calculated using both PCL and QML estimators (with the covariance matrix from HFPS2). Theoretical spectra are plotted (in black) for τ = 0.05, 0.07, and 0.09. 

In the text 
Fig. 33 100 × 143 crossspectra used in this paper. Top: debiased XPol PCL (green), and the biased (red) and debiased (blue) SimBaL PCL spectra. Middle: QML power spectra (red biased, blue debiased) for HFPS1. Bottom: same for HFPS2. Model spectra for τ = 0.050, 0.070, and 0.090 are displayed in black (dashed) lines. 

In the text 
Fig. 34 Expected EE crossspectrum statistical distributions for 2 ≤ ℓ ≤ 7, computed from the HFPS1 simulations, with a fiducial τ = 0.055. The blue (PCL) and red (QML) vertical lines show the observed values. The probability to exceed (PTE) is given in each panel. 

In the text 
Fig. 35 PTE expressed as equivalent χ^{2} for both the QML and PCL estimators, derived from the statistic shown in Fig. 34. 

In the text 
Fig. 36 TT and TE100 × 143 crossspectra, plotted for the SimBaL results with and without the bias correction. The black lines shows the fiducial spectra for τ = 0.05, 0.07, and 0.09. 

In the text 
Fig. 37 τ posteriors computed by SimBal from the 100 and 143GHz QML cross and autospectra, both with (solid curves) and without (dashed curves) debiasing by the mean of the simulated power spectra. While the 100 × 143 crosspower spectrum is not affected by debiasing, the 100 × 100 and 143 × 143 autopower spectra, which are still dominated by the systematic effects, are more affected. Nevertheless all debiased τ posteriors are consistent. 

In the text 
Fig. 38 Left column: simulation of SimBaL posterior distributions using the HFPS sets for noise and systematic effects, with the fiducial value of τ = 0.060. Right column: distribution of the τ peak value from the left column. The top row is for the PCL estimator and HFPS1 simulations, and the bottom row for the QML estimator and HFPS2 simulations. 

In the text 
Fig. 39 τ posteriors obtained with the two likelihood methods (Lollipop and SimBal) for the 100 × 143 QML and PCL crossspectra estimators, colour coded as indicated in the legend. 

In the text 
Fig. 40 Posterior distributions of τ calculated by SimBaL from LFI × HFI crossspectra. Top: results from simulations using 10 LFI systematic effect realizations, 83 HFPS1 realizations, and 100 CMB realizations with input reionization parameter τ = 0.05. Bottom: results from data, specifically the HFI pre2016 maps at 100 and 143 GHz and the LFI 70GHz 2015 released maps (Planck Collaboration VI 2016). Results from 70 × 100 and 70 × 143 are consistent. 

In the text 
Fig. 41 History of τ determination with WMAP and Planck. We have omitted the first WMAP determination (τ = 0.17 ± 0.04, Bennett et al. 2003), which was based on TE alone. 

In the text 
Fig. 42 Parameter constraints for the base ΛCDM cosmology, illustrating the τ–n_{s} degeneracy and the impact of replacing the LFIbased lowP likelihood used in the 2015 Planck papers with the HFIbased SimLow likelihood discussed here. The values of τ and σ_{8} shift downwards. 

In the text 
Fig. A.1 Estimated contribution of glitch residuals to autopower spectra for simulations of the 143detset1, with the FEE model in black. The red curve shows the evaluation of the power spectra of the difference of maps of simulations with and without glitchdeglitch removal. Spectra here are averaged over four realizations. The blue curve shows the spectra after processing for the simulated signal, with all systematic effects included. The green curve is the autopower spectra of half ring differences, dominated by noise at high mutlipoles and by glitches at low multipoles. 

In the text 
Fig. A.2 Simulated contribution of warm readout electronics drifts on the autopower spectra of 143detset1. The red curves show the contribution of the electronics drift compared to the total map power (orange) and the FEE input to the simulations (black). The green curve is an estimate of the noise level based on the halfring differences. 

In the text 
Fig. A.3 Autopower spectra EE (left) and BB (right) of the 217 GHz endtoend simulations carried out to test the effects of detector crosstalk. Red crosses indicate the FEE sky signal input to the simulation, to which a realization drawn from the noise variance is added (dashed black). The orange and blue crosses show the total power in the simulations, with and without simulated crosstalk, respectively. The absolute value of the difference (between the simulations including crosstalk and those without crosstalk) is shown in green. We note the difference in binning data between left and right plots. 

In the text 
Fig. B.1 Overview and main functional tasks of SRoll. 

In the text 
Fig. B.2 Gains solved by SRoll for simulated ADC nonlinearities, using gain periods of uniform length (red) and optimized sampling (green). The blue curve shows the apparent gain variation fitted between a simulated total dipole and the same total dipole distorted by the simulated ADC nonlinearities. 

In the text 
Fig. B.3 Histogram of the η parameter measured on real data. For each gainstable period, the total dipole is fitted on the Galactic signal. 

In the text 
Fig. B.4 Masks used during the mapmaking solution. The fraction of the sky used is f_{sky} = 0.862 at 100 GHz, f_{sky} = 0.856 at 143 GHz, f_{sky} = 0.846 at 217 GHz, and f_{sky} = 0.862 at 353 GHz. 

In the text 
Fig. B.5 Power spectra of the differences between maps made from TOI including the CMB, and maps made from TOI including only noise and systematics with the same CMB added to the map. In this simulation run, the CMB power spectra are set to FEE. 

In the text 
Fig. B.6 Evolution of the gain solved by SRoll for four representative bolometers. The left column shows gain as a function of signal level in ADU, with values determined from the data in blue. The red line is the best fit, and the pink area surrounding it is the standard deviation of a fourparameter model of the residual ADC nonlinearity fitted to these data. The right column shows gain as a function of pointing period, and shows that the model is able to reproduce the broad features of the remaining gain variation observed during the mission. 

In the text 
Fig. B.7 Error on the recovered leakage coefficients solved by SRoll, normalized to the average of all detectors in a frequency band, on a representative set of detectors. The error bars (in red) are the statistical distribution from 20 realizations of the simulations. 

In the text 
Fig. B.8 Residual autopower spectra of the CO bandpass leakage residuals. 

In the text 
Fig. B.9 Residual autopower spectra of the dust bandpass leakage residuals. 

In the text 
Fig. B.10 Residual autopower spectra of the freefree bandpass leakage residuals. 

In the text 
Fig. B.11 Gain variations solved by SRoll (red line; the light red band shows the variance of 20 noise simulation realizations) compared to the gain variation computed directly from the input ADC model (blue). 

In the text 
Fig. B.12 Residual autopower spectrum of the ADC nonlinearity corrected by the gain variation model. 

In the text 
Fig. B.13 Residual polarization maps from simulations at 100 GHz. Maps in the left column are obtained by using a constant gain per ring, while maps in the right column are obtained when the simulations are run using the ADC nonlinearity model. 

In the text 
Fig. B.14 Q maps of data (left column) and a simulation of the ADC nonlinearity effects (right column) at 100 GHz. The first row shows the total effect. The second row shows the improvement brought about by the correction of the linear gain variation. The third row shows the noise and residuals from other systematics. 

In the text 
Fig. B.15 Imaginary part of the complex transfer function residuals, simulated for all bolometers and for three ranges of spinfrequency harmonics. Errors bars are computed from 20 noise simulations. 

In the text 
Fig. B.16 Simulated power spectra of the transfer function residuals compared to the fiducial CMB signal (black curve). For each bolometer, the amplitude of the empirical transfer function was drawn from the uncertainties shown in Fig. B.15. The residuals are negligible compared to the signal. 

In the text 
Fig. B.17 Power spectra of null tests on detsets (left column) and half missions (right column). The data are shown in blue and the mean and standard deviation from simulations containing noise alone are shown in red. 

In the text 
Fig. B.18 Autopower spectra induced by the residual calibration mismatch in the maps. The residuals in the CMB channels are very low compared to the FEE signal. 

In the text 
Fig. D.1 Distribution of the SimBaL estimator from signalfree HPFS1 and HFPS1+HPFS3 simulations. The two distributions (blue and green dashed lines) are consistent with a Gaussian width σ = 0.015 for the relevant numbers of simulations. The value from the data (in red) is incompatible with these simulations at the 3.5σ level. 

In the text 
Fig. D.2 Distribution of the τ probability density as a function of the estimator (defined in the text) derived from the 100 × 143 crossspectra. The top panel uses simulations containing only white noise, while the bottom panel uses simulations with the systematic effects residuals added. As an example of a cut through this twodimensional plot, in the top panel an observed value of (horizontal black line) gives a probability distribution (red curve) that peaks at τ = 0.060. In the bottom panel, the same peak value of the probability distribution is obtained with an observed . 

In the text 
Fig. D.3 Pixelpixel matrix χ^{2} values of the first 400 of the 143 GHz FFP8 simulations, evaluated against 18arbitraryeigenmode fits to either the first 50, 83, or 200 simulations. The χ^{2} of the simulation data set used to compute the pixelpixel covariance matrix has been normalized here. Due to the limited number of simulations, the real variance, estimated from the simulations that were not used to compute the pixelpixel matrix, is significantly higher (i.e., the red line jumps at about 50, the green line at about 83, and the blue line at about 200). 

In the text 
Fig. D.4 Probability slice at ℓ = 4 for two different τ values, τ = 0.025 (top panels) and τ = 0.06 (bottom panels). The left panels show the signalonly distribution, while the right ones show the signal plus noise and systematic effects. For small τ values, the distribution is dominated by noise and systematic effects and is fitted poorly by the HamimecheLewistype model (red line), but is wellcaptured by the SimBaL polynomial approximation (blue curve). 

In the text 
Fig. D.5 Posterior distributions of τ from SimBaL (blue curve) and SimBaL_HL (red curve), showing how they differ at low τ. The red line is an attempt to mimic in SimBaL the modified “HL” approach, and demonstrates that the divergence is due to residual nonGaussian behaviour of the statistics at low multipoles (see text). 

In the text 
Fig. D.6 Masks for different values of f_{sky}. The mask shown in the bottomleft panel (50% sky fraction) is the one used in the Lollipop and the SimBaL analyses. The other three were used to evaluate the dependence of the results on changing sky fractions from 30 to 60%. 

In the text 
Fig. D.7 Posterior distributions from SimBaL PCL 100 × 143 crosspower spectra obtained with different sky fractions, showing that variations caused by masking the Galaxy have little effect on the τ posterior likelihood for a range of mask sizes. 

In the text 
Fig. D.8 Posteriors on τ using different multipole ranges in the 100 × 143 PCL crosspower spectrum. The ℓ ranges are indicated in the legend. 

In the text 
Fig. D.9 Posteriors in τ when removing one multipole at a time in the 100 × 143 QML crosspower spectra. The discrepant curve is when removing ℓ = 5. 

In the text 
Fig. D.10 Histograms of the peak values of the τ posterior distributions after removing one multipole each time, computed from 8300 simulations with τ = 0.055. For the red curve the maximum value of the 100 × 143 QML power spectrum in the range ℓ = 2−12 is removed, and for the blue curve ℓ = 5 is removed when it happens to be the maximum. The peak value of the τ posterior when removing ℓ = 5 from the data is shown in green, and is consistent with the blue curve, as expected. 

In the text 
Fig. D.11 computed for the data using SimLow. Power spectrum values computed for different τ are shown in blue. For ℓ = 3 and ℓ = 5 we can see that very small values of τ are excluded, while for ℓ = 4 large values of τ have a low probability. 

In the text 
Fig. D.12 Posterior distribution for τ computed with the SimLow likelihood using the same cosmological parameters as for SimBaL. The posterior is consistent with the LowEH result. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.