© ESO, 2016

1. Introduction

One of a set associated with the 2015 release of data from the Planck¹ mission, this paper describes the techniques we employed to calibrate the voltages measured by the Low Frequency Instrument (LFI) radiometers into a set of thermodynamic temperatures, which we refer to as photometric calibration. We expand on the work described in Planck Collaboration V (2014), henceforth Cal13; we try to follow the structure of the earlier paper as closely as possible to help the reader understand what has changed between the 2013 and the 2015 Planck data releases.

The calibration of both Planck instruments (for HFI, see Planck Collaboration VII 2016) is now based on the small (270 μK) dipole signal induced by the annual motion of the satellite around the Sun – the orbital dipole, which we derive from our knowledge of the orbital parameters of the spacecraft. The calibration is thus absolute and does not depend on external measurements of the larger solar (3.35 mK) dipole, as was the case for Cal13. Absolute calibration allows us both to improve the current measurement of the solar dipole (see Sect. 5) and to transfer Planck’s calibration to various ground-based instruments (see, e.g., Perley & Butler 2015) and other cosmic microwave background (CMB) experiments (e.g., Louis et al. 2014).

Accurate calibration of the LFI is crucial for ensuring reliable cosmological and astrophysical results from the Planck mission. Internally consistent photometric calibration of the nine Planck frequency channels is essential in component separation, where we disentangle the CMB from the various Galactic and extragalactic foreground emission processes (Planck Collaboration IX 2016; Planck Collaboration X 2016). In addition, the LFI calibration directly affects the Planck polarization likelihood at low multipoles, based on the LFI 70 GHz channel, which is extensively employed in the cosmological analysis of this 2015 release. Furthermore, a solid absolute calibration is needed to compare and combine Planck data with results from other experiments, most notably with WMAP. Detailed comparisons of the calibrated data from single LFI radiometers of the three LFI frequency channels, and between LFI and HFI, allow us to test the internal consistency and accuracy of our calibration.

In this paper, we quantify both the absolute and relative accuracy in the calibration of the LFI instrument and find an overall uncertainty of 0.35% (30 GHz map), 0.26% (44 GHz), and 0.20% (70 GHz). The level of the power spectrum near the first peak is now remarkably consistent with WMAP’s. Other papers in this Planck data release deal with the quality of the LFI calibration, in particular:

• Planck Collaboration X (2016)quantifies the consistency between the calibration of theLFI/HFI/WMAP channels in the context of foregroundcomponent separation, finding that the measured discrepanciesamong channels are a few tenths of a percent;

• Planck Collaboration XI (2016) analyses the consistency between the LFI 70 GHz low-ℓ polarization map and the WMAP map in pixel space, finding no hints of inconsistencies;

• Planck Collaboration XIII (2016) compares the estimate for the τ and z_re cosmological parameters (reionization optical depth and redshift) using either LFI 70 GHz polarization maps or WMAP maps and finds statistically consistent values.

To achieve calibration accuracy at the few-per-thousand level requires careful attention to instrumental systematic effects and foreground contamination of the orbital dipole. Much of this paper is devoted to discussing these effects and the means to mitigate them.

In this paper we do not explicitly discuss polarization-related issues. Although polarization analysis is one of the most important results of this data release, the calibration of the LFI radiometers is inherently based on temperature signals (Leahy et al. 2010). Estimates of the sensitivity in polarization, as well as the impact of calibration-related systematics on it, are provided by Planck Collaboration III (2016).

A schematic of the LFI pseudo-correlation receiver is shown in Fig. 1. We model the output voltage V(t) of each radiometer as $\begin{array}{ccc}\mathit{V}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathit{G}\mathrm{\left(}\mathit{t}\mathrm{\right)}\hspace{0.17em}[\mathit{B}\mathrm{\ast }\mathrm{(}{\mathit{T}}_{\mathrm{sky}}\mathrm{+}\mathit{D}\mathrm{)}]\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{+}\mathit{M}\mathrm{+}\mathit{N,}& & \end{array}$(1)where G is the gain (measured in VK^-1), B the beam response, D the thermodynamic temperature of the total CMB dipole signal (i.e., a combination of the solar and orbital components, including the quadrupolar relativistic corrections), which we use as a calibrator, and T_sky = T_CMB + T_Gal + T_other is the overall temperature of the sky (CMB anisotropies, diffuse Galactic emission and other² foregrounds, respectively) apart from D. Finally, M is a constant offset and N a noise term. In the following sections, we use Eq. (1) many times; whenever the presence of the N term is not important, it is silently dropped. The ∗ operator represents a convolution over the 4π sphere. We base our calibration on the knowledge of the spacecraft velocity around the Sun, which produces the orbital dipole, and use the orbital dipole to accurately measure the dominant solar dipole component. The purpose of this paper is to explain how we implemented and validated the pipeline that estimates the calibration constant K ≡ G^-1 (which is used to convert the voltage V into a thermodynamic temperature), to quantify the quality of our estimate for K, and to quantify the impact of possible systematic calibration errors on the Planck/LFI data products.

Several improvements were introduced in the LFI pipeline for calibration relative to Cal13. In Sect. 2 we recall some terminology and basic ideas presented in Cal13 to discuss the normalization of the calibration, i.e., what factors influence the average value of G in Eq. (1). Section 3 provides an overview of the new LFI calibration pipeline and underlines the differences with the pipeline described in Cal13. One of the most important improvements in the 2015 calibration pipeline is the implementation of a new iterative algorithm to calibrate the data, DaCapo. Its principles are presented separately in a dedicated section, Sect. 4. This code has also been used to characterize the orbital dipole. The details of this latter analysis are provided in Sect. 5, where we present a new characterization of the solar dipole. These two steps are crucial for calibrating LFI. Section 6 describes a number of validation tests we have run on the calibration, as well as the results of a quality assessment. This section is divided into several parts: in Sect. 6.1 we compare the overall level of the calibration in the 2015 LFI maps with those in the previous data release; in Sect. 6.2 we provide a brief account of the simulations described in Planck Collaboration III (2016), which assess the calibration error due to the white noise and approximations in the calibration algorithm itself; in Sect. 6.3 we describe how uncertainties in the shape of the beams might affect the calibration; Sects. 6.4 and 6.5 measure the agreement between radiometers and groups of radiometers in the estimation of the TT power spectrum; and Sect. 6.6 provides a reference to the discussion of null tests provided in Planck Collaboration III (2016). Finally, in Sect. 7, we derive an independent estimate of the LFI calibration from our measurements of Jupiter and discuss its consistency with our nominal dipole calibration.

2. Handling beam efficiency

In this section we develop a mathematical model that relates the absolute level of the calibration (i.e., the average level of the raw power spectrum $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{C}}_{\mathit{\ell }}\end{array}$ for an LFI map) to a number of instrumental parameters related to the beams and the scanning strategy.

The beam response B(θ,ϕ) is a dimensionless function defined over the 4π sphere. In Eq. (1), B appears in the convolution $\mathit{B}\mathrm{\ast }\mathrm{(}{\mathit{T}}_{\mathrm{sky}}\mathrm{+}\mathit{D}\mathrm{)}\mathrm{=}\frac{{\mathrm{\int }}_{\mathrm{4}\mathit{\pi }}\mathit{B}\mathrm{\left(}\mathit{\theta ,\varphi }\mathrm{\right)}\hspace{0.17em}\mathrm{(}{\mathit{T}}_{\mathrm{sky}}\mathrm{+}\mathit{D}\mathrm{)}\mathrm{\left(}\mathit{\theta ,\varphi }\mathrm{\right)}\hspace{0.17em}\mathrm{d\Omega }}{{\mathrm{\int }}_{\mathrm{4}\mathit{\pi }}\mathit{B}\mathrm{\left(}\mathit{\theta ,\varphi }\mathrm{\right)}\hspace{0.17em}\mathrm{d\Omega }}\mathit{,}$(2)whose value changes with time because of the change in orientation of the spacecraft. Since no time-dependent optical effects are evident from the data taken from October 2009 to February 2013 (Planck Collaboration IV 2016), we assume there is no intrinsic change in the shape of B during the surveys.

In the previous data release, we approximated B as a Dirac delta function (a pencil beam) when modelling the dipole signal seen by the LFI radiometers. The same assumption has been used for all the WMAP data releases (see, e.g., Hinshaw et al. 2009), as well as in the HFI pipeline (Planck Collaboration VIII 2016). However, the real shape of B deviates from the ideal case of a pencil beam because of two factors: (1) the main beam is more like a Gaussian with an elliptical section, whose FWHM (full width half maximum) ranges between 13′ and 33′ in the case of the LFI radiometers; and (2) farther than 5° from the beam axis, the presence of far sidelobes dilutes the signal measured through the main beam further and induces an axial asymmetry on B. Previous studies³ tackled the first point by applying a window function to the power spectrum computed from the maps to correct for the finite size of the main beam. However, the presence of far sidelobes might cause stripes in maps. For this 2015 data release, we used the full shape of B in computing the dipole signal adopted for the calibration. No significant variation in the level of the CMB power spectra with respect to the previous data release is expected, since we are basically subtracting power during the calibration process instead of reducing the level of the power spectrum by means of the window function. However, this new approach improves the internal consistency of the data, since the beam shape is taken into account from the very first stages of data processing (i.e., the signal measured by each radiometer is fitted with its own calibration signal B_rad ∗ D); see Sect. 6.6. The definition of the beam window function has been changed accordingly; see Planck Collaboration IV (2016).

In Cal13 we introduced the two quantities φ_D and φ_sky as a way to quantify the impact of a beam window function on the calibration⁴ and on the mapmaking process, respectively. Here we briefly summarize the theory, and we introduce new equations that are relevant for understanding the normalization of the new Planck-LFI results in this data release.

Because of the motion of the solar system in the CMB rest frame, the solar dipole D is given by $\mathit{D}\mathrm{\left(}{x}\mathit{,t}\mathrm{\right)}\mathrm{=}{\mathit{T}}_{\mathrm{CMB}}\left(\frac{\mathrm{1}}{\mathit{\gamma }\mathrm{\left(}\mathit{t}\mathrm{\right)}\hspace{0.17em}(\mathrm{1}\mathrm{-}{\beta }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{·}{x})}\mathrm{-}\mathrm{1}\right)\mathit{,}$(3)where T_CMB is the CMB monopole, β = v/c is the velocity of the spacecraft, and γ = (1−β²)^{− 1 / 2}. Each radiometer measures the signal D convolved with the beam response B, according to Eq. (2); therefore, in principle, each radiometer has a different calibration signal. Under the assumption of a Dirac delta shape for B, Cal13 shows that the estimate of the gain constant $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{G}\end{array}$ is related to the true gain G by the formula $\begin{array}{c}\mathrm{˜}\\ {\mathit{G}}^{\mathrm{pen}}\end{array}\mathrm{=}\mathit{G}(\mathrm{1}\mathrm{-}{\mathit{f}}_{\mathrm{sl}})(\mathrm{1}\mathrm{+}{\mathit{\phi }}_{\mathit{D}})\mathit{,}$(4)where ${\mathit{f}}_{\mathrm{sl}}\mathrm{=}\frac{{\mathrm{\int }}_{\mathit{\theta }\hspace{0.17em}\mathit{>}\hspace{0.17em}{\mathrm{5}}^{\mathrm{◦}}}\mathit{B}\hspace{0.17em}\mathrm{d\Omega }}{{\mathrm{\int }}_{\mathrm{4}\mathit{\pi }}\mathit{B}\hspace{0.17em}\mathrm{d\Omega }}$(5)is the fraction of power entering the sidelobes (i.e., along directions farther than 5° from the beam axis), and ${\mathit{\phi }}_{\mathit{D}}\mathrm{=}\frac{{\mathit{\partial }}_{\mathit{t}}{\mathit{B}}_{\mathrm{sl}}\mathrm{\ast }\mathit{D}}{{\mathit{\partial }}_{\mathit{t}}{\mathit{B}}_{\mathrm{main}}\mathrm{\ast }\mathit{D}}$(6)is a time-dependent quantity that depends on the shape of B = B_main + B_sl and its decomposition into a main (θ< 5°) and a sidelobe part, on signal D, and on the scanning strategy because of the time dependence of the stray light; the notation ∂_t indicates a time derivative. Once the timelines are calibrated, traditional mapmaking algorithms approximate⁵B as a Dirac delta (e.g., Hinshaw et al. 2003; Jarosik et al. 2007; Keihänen et al. 2010), thus introducing a new systematic error. In this case, the mean temperature $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{T}}_{\mathrm{sky}}\end{array}$ of a pixel in the map would be related to the true temperature T_sky by the formula ${\mathit{T}}_{\mathrm{sky}}\mathrm{=}\begin{array}{c}\mathrm{˜}\\ {\mathit{T}}_{\mathrm{sky}}^{\mathrm{pen}}\end{array}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\phi }}_{\mathrm{sky}}\mathrm{+}{\mathit{\phi }}_{\mathit{D}}\mathrm{)}\mathit{,}$(7)which applies to timelines and should only be considered valid when considering details on angular scales larger than the width of the main beam. Cal13 defines the quantity φ_sky using the following equation: ${\mathit{\phi }}_{\mathrm{sky}}\mathrm{=}\frac{{\mathit{B}}_{\mathrm{sl}}\mathrm{\ast }{\mathit{T}}_{\mathrm{sky}}}{{\mathit{B}}_{\mathrm{main}}\mathrm{\ast }{\mathit{T}}_{\mathrm{sky}}}\left(\frac{{\mathit{T}}_{\mathrm{sky}}}{\begin{array}{c}\mathrm{˜}\\ {\mathit{T}}_{\mathrm{sky}}\end{array}}\right)\mathrm{=}\frac{{\mathit{B}}_{\mathrm{sl}}\mathrm{\ast }{\mathit{T}}_{\mathrm{sky}}}{\begin{array}{c}\mathrm{˜}\\ {\mathit{T}}_{\mathrm{sky}}\end{array}}\mathrm{·}$(8)See Fig. 2 for an example showing how φ_sky is computed.

In this 2015 Planck data release, we take advantage of our knowledge of the shape of B to compute the value of Eq. (2)and use this as our calibrator. Since the term B ∗ T_sky = (B_main + B_sl) ∗ T_sky is unknown, we apply the following simplifications:

• 1.

  we apply the point source and 80% Galactic masks (PlanckCollaboration 2015), in order not to consider theB_main ∗ T_sky term in the computation of the convolution;

• 2.

  we assume that B_sl ∗ T_sky ≈ B_sl ∗ T_Gal and subtract it from the calibrated timelines, using an estimate for T_Gal computed by means of models of the Galactic emission (Planck Collaboration IX 2016; Planck Collaboration X 2016).

The result of such transformations is a new timeline ${\mathit{V}}_{\mathrm{out}}^{\mathrm{\prime }}$. Under the hypothesis of perfect knowledge of the beam B and of the dipole signal D, these steps are enough to estimate the true calibration constant without bias⁶ (unlike Eq. (4)): $\stackrel{}{{{\mathit{G}}^{\mathrm{4}\mathit{\pi }}}_{\mathrm{˜}}}\mathrm{=}\mathit{G,}$(9)which should be expected, since no systematic effects caused by the shape of B affect the estimate of the gain G. To see how Eq. (7)changes in this case, we write the measured temperature $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{T}}_{\mathrm{sky}}\end{array}$ as $\begin{array}{c}\mathrm{˜}\\ {\mathit{T}}_{\mathrm{sky}}\end{array}\mathrm{=}\mathit{B}\mathrm{\ast }{\mathit{T}}_{\mathrm{sky}}\mathrm{+}\mathit{M}\mathrm{=}{\mathit{B}}_{\mathrm{main}}\mathrm{\ast }{\mathit{T}}_{\mathrm{sky}}\mathrm{+}{\mathit{B}}_{\mathrm{sl}}\mathrm{\ast }{\mathit{T}}_{\mathrm{sky}}\mathrm{+}\mathit{M}\mathit{.}$(10)Since in this 2015 data release we remove B_sl ∗ T_Gal, the contribution of the pickup of Galactic signal through the sidelobes (Planck Collaboration II 2016), the equation can be rewritten as $\stackrel{}{{{\mathit{T}}_{\mathrm{sky}}^{\mathrm{4}\mathit{\pi }}}_{\mathrm{˜}}}\mathrm{=}{\mathit{B}}_{\mathrm{main}}\mathrm{\ast }{\mathit{T}}_{\mathrm{sky}}\mathrm{+}{\mathit{B}}_{\mathrm{sl}}\mathrm{\ast }\mathrm{(}{\mathit{T}}_{\mathrm{CMB}}\mathrm{+}{\mathit{T}}_{\mathrm{other}}\mathrm{)}\mathrm{+}\mathit{M}\mathit{.}$(11)If we neglect details at angular scales smaller than the main beam size, then ${\mathit{B}}_{\mathrm{main}}\mathrm{\ast }{\mathit{T}}_{\mathrm{sky}}\mathrm{\approx }\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{f}}_{\mathrm{sl}}\mathrm{)}\hspace{0.17em}{\mathit{T}}_{\mathrm{sky}}\mathit{,}$(12)so that $\stackrel{}{{{\mathit{T}}_{\mathrm{sky}}^{\mathrm{4}\mathit{\pi }}}_{\mathrm{˜}}}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{f}}_{\mathrm{sl}}\mathrm{)}\hspace{0.17em}{\mathit{T}}_{\mathrm{sky}}\mathrm{+}{\mathit{B}}_{\mathrm{sl}}\mathrm{\ast }\mathrm{(}{\mathit{T}}_{\mathrm{CMB}}\mathrm{+}{\mathit{T}}_{\mathrm{other}}\mathrm{)}\mathrm{+}\mathit{M}\mathit{.}$(13)We modify Eq. (8)in order to introduce a new term ${\mathit{\phi }}_{\mathrm{sky}}^{\mathrm{\prime }}$: ${\mathit{\phi }}_{\mathrm{sky}}^{\mathrm{\prime }}\mathrm{=}\frac{{\mathit{B}}_{\mathrm{sl}}\mathrm{\ast }\mathrm{(}{\mathit{T}}_{\mathrm{CMB}}\mathrm{+}{\mathit{T}}_{\mathrm{other}}\mathrm{)}}{\begin{array}{c}\mathrm{˜}\\ {\mathit{T}}_{\mathrm{sky}}\end{array}}\mathrm{·}$(14)When solving for T_sky, Eq. (13)can be rewritten as ${\mathit{T}}_{\mathrm{sky}}\mathrm{=}\stackrel{}{{{\mathit{T}}_{\mathrm{sky}}^{\mathrm{4}\mathit{\pi }}}_{\mathrm{˜}}}\hspace{0.17em}\frac{\mathrm{1}\mathrm{-}{\mathit{\phi }}_{\mathrm{sky}}^{\mathrm{\prime }}}{\mathrm{1}\mathrm{-}{\mathit{f}}_{\mathrm{sl}}}\mathrm{+}{\mathit{T}}_{\mathrm{0}}\mathit{,}$(15)where T₀ = M/ (1−f_sl) is a constant offset that is not very relevant for pseudo-differential instruments like LFI. Equation (15)is the equivalent of Eq. (7)in the case of a calibration pipeline that takes the 4π shape of B into account, as is the case for the Planck-LFI pipeline used for the 2015 data release.

Since one of the purposes of this paper is to provide a quantitative comparison of the calibration of this Planck data release with the previous one, we provide now a few formulae that quantify the change in the average level of the temperature fluctuations and of the power spectrum between the 2013 and 2015 releases. The variation in temperature can be derived from Eqs. (7)and (15): $\frac{\stackrel{}{{{\mathit{T}}_{\mathrm{sky}}^{\mathrm{2015}}}_{\mathrm{˜}}}}{\begin{array}{c}\mathrm{˜}\\ {\mathit{T}}_{\mathrm{sky}}^{\mathrm{2013}}\end{array}}\mathrm{=}\frac{(\mathrm{1}\mathrm{-}{\mathit{\phi }}_{\mathrm{sky}}\mathrm{+}{\mathit{\phi }}_{\mathit{D}})(\mathrm{1}\mathrm{-}{\mathit{f}}_{\mathrm{sl}})}{\mathrm{1}\mathrm{-}{\mathit{\phi }}_{\mathrm{sky}}^{\mathrm{\prime }}}\mathrm{\approx }\frac{\mathrm{1}\mathrm{-}{\mathit{f}}_{\mathrm{sl}}\mathrm{-}{\mathit{\phi }}_{\mathrm{sky}}\mathrm{+}{\mathit{\phi }}_{\mathit{D}}}{\mathrm{1}\mathrm{-}{\mathit{\phi }}_{\mathrm{sky}}^{\mathrm{\prime }}}\mathrm{·}$(16)If we consider the ratio between the power spectra $\stackrel{}{{{\mathit{C}}_{\mathit{\ell }}^{\mathrm{2015}}}_{\mathrm{˜}}}$ and $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathrm{2013}}\end{array}$, the quantity becomes $\frac{\stackrel{}{{{\mathit{C}}_{\mathit{\ell }}^{\mathrm{2015}}}_{\mathrm{˜}}}}{\begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathrm{2013}}\end{array}}\mathrm{\approx }{\left(\frac{\mathrm{1}\mathrm{-}{\mathit{f}}_{\mathrm{sl}}\mathrm{-}{\mathit{\phi }}_{\mathrm{sky}}\mathrm{+}{\mathit{\phi }}_{\mathit{D}}}{\mathrm{1}\mathrm{-}{\mathit{\phi }}_{\mathrm{sky}}^{\mathrm{\prime }}}\right)}^{\mathrm{2}}\mathrm{·}$(17)In Sect. 6.1 we provide quantitative estimates of f_sl, φ_D, φ_sky, and ${\mathit{\phi }}_{\mathrm{sky}}^{\mathrm{\prime }}$, as well as the ratios in Eqs. (16)and (17).

3. The calibration pipeline

In this section we briefly describe the implementation of the calibration pipeline. Readers interested in more detail should refer to Planck Collaboration II (2016).

Evaluating the calibration constant K (see Eq. (1)) requires us to fit the timelines of each radiometer with the expected signal D induced by the dipole as Planck scans the sky. This process provides the conversion between the voltages and the measured thermodynamic temperature.

As discussed in Sect. 2, we have improved the model used for D, since we are now computing the convolution of D with each beam B over the full 4π sphere. Moreover, we are considering the B_sl ∗ T_Gal term in the fit in order to reduce the bias due to the pickup of Galactic signal by the beam far sidelobes. The model of the dipole D now includes the correct⁷ quadrupolar corrections for special relativity. The quality of the beam estimate B has been improved as well: we are now using all seven Jupiter transits observed in the full four-year mission, and we account for the optical effects of the variation in the beam shape across the band of the radiometers. It is important to underline that these new beams do not follow the same normalization convention as in the first data release (now ^∫_4πB(θ,ϕ) dΩ ≠ 1) because numerical inaccuracies in the simulation of the 4π beams cause a loss of roughly 1% of the signal entering the sidelobes⁸: see Planck Collaboration IV (2016) for a discussion of this point.

As for the 2013 data release, the calibration constant K is estimated once per each pointing period, i.e., the period during which the spinning axis of the spacecraft holds still and the spacecraft rotates at a constant spinning rate of 1 / 60 Hz. The code used to estimate K, named DaCapo, has been completely rewritten; it is able to run in two modes, one of which (the so-called unconstrained mode) is able to produce an estimate of the solar dipole signal, and the other one (the constrained mode) which requires the solar dipole parameters as input. We used the unconstrained mode to assess the characteristics of the solar dipole, which were then used as input into the constrained mode of DaCapo for producing the actual calibration constants.

We smooth the calibration constants produced by DaCapo by means of a running mean, where the window size has a variable length. That length is chosen so that every time there is a sudden change in the state of the instrument (e.g., because of a change in the thermal environment of the front-end amplifiers) that discontinuity is not averaged out. However, this kind of filter removes any variation in the calibration constants, whose timescale is smaller than a few weeks. One example of this latter kind of fluctuation is the daily variation measured in the radiometer backend gains during the first survey, which was caused by the continuous turning on-and-off of the transponder⁹ while sending the scientific data to Earth once per day. To keep track of these fluctuations, we estimated the calibration constants K using the signal of the 4.5 K reference load in a manner similar to that described in Cal13 under the name of 4 K calibration, and we have added this estimate to the DaCapo gains after having applied a high-pass filter to them, as shown in Fig. 5. Details about the implementation of the smoothing filter are provided in Appendix A.

Once the smoothing filter has been applied to the calibration constants K, we multiply the voltages by K in order to convert them into thermodynamic temperatures and remove the term B ∗ D + B_sl ∗ T_sky from the result, thus removing the dipole and the Galactic signal captured by the far sidelobes from the data. The value for T_sky has been taken from a sum of the foreground signals considered in the simulations described in Planck Collaboration IX (2016); refer to Planck Collaboration II (2016) for further details.

4. The calibration algorithm

DaCapo is an implementation of the calibration algorithm we used in this data release to produce an estimate of the calibration constant K in Eq. (1). In this section we describe the model on which DaCapo is based, as well as a few details of its implementation.

4.1. Unconstrained algorithm

Let V_i be the ith sample of an uncalibrated data stream, and k(i) the pointing period to which the sample belongs. Following Eq. (1)and assuming the usual mapmaking convention of scanning the sky, T_sky, using a pencil beam, we model the uncalibrated time stream as ${\mathit{V}}_{\mathit{i}}\mathrm{=}{\mathit{G}}_{\mathit{k}\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathrm{(}{\mathit{T}}_{\mathit{i}}\mathrm{+}\mathit{B}\mathrm{\ast }{\mathit{D}}_{\mathit{i}}\mathrm{)}\mathrm{+}{\mathit{b}}_{\mathit{k}\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathrm{+}{\mathit{N}}_{\mathit{i}}\mathit{,}$(18)where we write B ∗ D_i ≡ (B ∗ D)_i and use the shorthand notation T_i = (T_sky)_i. The quantity G_k is the unknown gain factor for kth pointing period, n_i represents white noise, and b_k is an offset¹⁰ that captures the correlated noise component. We denote the sky signal by T_i , which includes foregrounds and the CMB sky apart from the dipole, and the dipole signal as seen by the beam B by B ∗ D_i. The dipole includes both the solar and orbital components, and it is convolved with the full 4π beam. The beam convolution is carried out by an external code, and the result is provided as input to DaCapo.

The signal term is written with help of a pointing matrix P as $\begin{array}{ccc}{\mathit{T}}_{\mathit{i}}\mathrm{=}\sum _{\mathit{p}}{\mathit{P}}_{\mathit{ip}}{\mathit{m}}_{\mathit{p}}\mathit{.}& & \end{array}$(19)Here P is a pointing matrix that picks the time-ordered signal from the unknown sky map m. The current implementation takes only the temperature component into consideration. In radiometer-based calibration, however, the polarization signal is partly accounted for, since the algorithm interprets whatever combination of the Stokes parameters (I,Q,U) the radiometer records as temperature signal. In regions that are scanned in only one polarization direction, this gives a consistent solution that does not induce any error on the gain. A small error can be expected to arise in those regions where the same sky pixel is scanned in vastly different directions of polarization sensitivity. The error is proportional to the ratio of the polarization signal and the total sky signal, including the dipole.

We determine the gains by minimizing the quantity ${\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\sum _{\mathit{i}}\frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}}({\mathit{V}}_{\mathit{i}}\mathrm{-}{\mathit{V}}_{\mathit{i}}^{\mathrm{mod}}{)}^{\mathrm{2}}\mathit{,}$(20)where ${\mathit{V}}_{\mathit{i}}^{\mathrm{mod}}\mathrm{=}{\mathit{G}}_{\mathit{k}\mathrm{\left(}\mathit{i}\mathrm{\right)}}\left(\sum _{\mathit{p}}{\mathit{P}}_{\mathit{ip}}{\mathit{m}}_{\mathit{p}}\mathrm{+}\mathit{B}\mathrm{\ast }{\mathit{D}}_{\mathit{i}}\right)\mathrm{+}{\mathit{b}}_{\mathit{k}\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathit{,}$(21)and ${\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}$ is the white noise variance. The unknowns of the model are m, G, b, and n (while we assume that the beam B is perfectly known). The dipole signal D and pointing matrix P are assumed to be known.

To reduce the uncertainty that arises from beam effects and subpixel variations in signal, we apply a galactic mask and include only those samples that fall outside the mask in the sum in Eq. (20).

Since Eq. (21)is quadratic in the unknowns, the minimization of χ² requires iteration. To linearize the model, we first rearrange it as $\begin{array}{ccc}{\mathit{V}}_{\mathit{i}}^{\mathrm{mod}}& \mathrm{=}& {\mathit{G}}_{\mathit{k}\mathrm{\left(}\mathit{i}\mathrm{\right)}}\left(\mathit{B}\mathrm{\ast }{\mathit{D}}_{\mathit{i}}\mathrm{+}\sum _{\mathit{p}}{\mathit{P}}_{\mathit{ip}}{\mathit{m}}_{\mathit{p}}^{\mathrm{0}}\right)\mathrm{+}{\mathit{G}}_{\mathit{k}\mathrm{\left(}\mathit{i}\mathrm{\right)}}^{\mathrm{0}}\sum _{\mathit{p}}{\mathit{P}}_{\mathit{ip}}\mathrm{(}{\mathit{m}}_{\mathit{p}}\mathrm{-}{\mathit{m}}_{\mathit{p}}^{\mathrm{0}}\mathrm{)}\\ & & \mathrm{+}\mathrm{[}\mathrm{(}{\mathit{G}}_{\mathit{k}\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathrm{-}{\mathit{G}}_{\mathit{k}\mathrm{\left(}\mathit{i}\mathrm{\right)}}^{\mathrm{0}}\mathrm{)}\mathrm{(}{\mathit{m}}_{\mathit{p}}\mathrm{-}{\mathit{m}}_{\mathit{p}}^{\mathrm{0}}{\mathrm{)}}^{\mathrm{]}}\mathrm{+}{\mathit{b}}_{\mathit{k}\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathit{.}\end{array}$(22)Here G⁰ and m⁰ are the gains and the sky map from the previous iteration step. We drop the quadratic term in brackets and obtain $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{V}}_{\mathit{i}}^{\mathrm{mod}}\mathrm{=}{\mathit{G}}_{\mathit{k}\mathrm{\left(}\mathit{i}\mathrm{\right)}}\left(\mathit{B}\mathrm{\ast }{\mathit{D}}_{\mathit{i}}\mathrm{+}\sum _{\mathit{p}}{\mathit{P}}_{\mathit{ip}}{\mathit{m}}_{\mathit{p}}^{\mathrm{0}}\right)\mathrm{+}{\mathit{G}}_{\mathit{k}\mathrm{\left(}\mathit{i}\mathrm{\right)}}^{\mathrm{0}}\sum _{\mathit{p}}{\mathit{P}}_{\mathit{ip}}\begin{array}{c}\mathrm{˜}\\ {\mathit{m}}_{\mathit{p}}\end{array}\mathrm{+}{\mathit{b}}_{\mathit{k}\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathit{.}\end{array}& & \end{array}$(23)Here $\begin{array}{c}\mathrm{˜}\\ {\mathit{m}}_{\mathit{p}}\end{array}\mathrm{=}\mathrm{(}{\mathit{m}}_{\mathit{p}}\mathrm{-}{\mathit{m}}_{\mathit{p}}^{\mathrm{0}}\mathrm{)}$(24)is a correction to the map estimate from the previous iteration step. Equation (23) is linear in the unknowns , G and b. We run an iterative procedure, where at each step we minimize χ² with the linearized model in Eq. (23), update the map and the gains as ${\mathit{m}}^{\mathrm{0}}\mathrm{\to }{\mathit{m}}^{\mathrm{0}}\mathrm{+}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{m}\end{array}$ and G⁰ → g, and repeat until convergence. The iteration is started from G⁰ = m⁰ = 0. Thus at the first step we are fitting just the dipole model and a baseline G_k(i)B ∗ D_i + b_k, and we obtain the first estimate for the gains. The first map estimate is obtained in the second iteration step.

DaCapo solves the gains for two radiometers of a horn at the same time. Two map options are available. Either the radiometers have their own sky maps, or both see the same sky. In the former case the calibrations become independent.

4.1.1. Solution of the linear system

Minimization of χ² yields a large linear system. The number of unknowns is dominated by the number of pixels in map m. It is possible, however, to reformulate the problem as a much smaller system as follows.

We first rewrite the model using matrix notation. We combine the first and last terms of Eq. (23) formally into ${\mathit{G}}_{\mathit{k}\mathrm{\left(}\mathit{i}\mathrm{\right)}}\left(\mathit{B}\mathrm{\ast }{\mathit{D}}_{\mathit{i}}\mathrm{+}\sum _{\mathit{p}}{\mathit{P}}_{\mathit{ip}}{\mathit{m}}_{\mathit{p}}^{\mathrm{0}}\right)\mathrm{+}{\mathit{b}}_{\mathit{k}\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathrm{=}\sum _{\mathit{j}}{\mathit{F}}_{\mathit{ij}}{\mathit{a}}_{\mathit{j}}\mathit{.}$(25)The vector a_j contains the unknowns b and G, and the matrix F spreads them into a time-ordered data stream. The dipole signal B ∗ D seen by the beam B, and a signal picked from map m⁰, are included in F.

Equation (21)can now be written in matrix notation as ${{V}}^{\mathrm{model}}\mathrm{=}\begin{array}{c}\mathrm{˜}\\ {P}\end{array}\begin{array}{c}\mathrm{˜}\\ {m}\end{array}\mathrm{+}{F}{a}\mathit{.}$(26)Gains G⁰ have been transferred inside matrix $\begin{array}{c}{\mathrm{˜}}_{}\\ {P}\end{array}$, ${\begin{array}{c}\mathrm{˜}\\ {P}\end{array}}_{\mathit{ip}}\mathrm{=}{\mathit{G}}_{\mathit{k}\mathrm{\left(}\mathit{i}\mathrm{\right)}}^{\mathrm{0}}{\mathit{P}}_{\mathit{ip}}\mathit{.}$(27)Using this notation, Eq. (20)becomes ${\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\mathrm{(}{V}\mathrm{-}\begin{array}{c}\mathrm{˜}\\ {P}\end{array}\begin{array}{c}\mathrm{˜}\\ {m}\end{array}\mathrm{-}{F}{a}{\mathrm{)}}^{\mathit{T}}{\mathit{C}}_{\mathit{n}}^{-1}\mathrm{(}{V}\mathrm{-}\begin{array}{c}\mathrm{˜}\\ {P}\end{array}\begin{array}{c}\mathrm{˜}\\ {m}\end{array}\mathrm{-}{F}{a}\mathrm{)}\mathit{,}$(28)where C_n is the white noise covariance.

Equation (28)is equivalent to the usual destriping problem of map-making (Planck Collaboration VI 2016), only the interpretation of the terms is slightly different. In place of the pointing matrix P we have $\begin{array}{c}{\mathrm{˜}}_{}\\ {P}\end{array}$, which contains the gains from the previous iteration step, and a contains the unknown gains beside the usual baseline offsets.

We minimize Eq. (28)with respect to $\begin{array}{c}{\mathrm{˜}}_{}\\ {m}\end{array}$, insert the result back into Eq. (28), and minimize with respect to a. The solution is identical to the destriping solution $\hat{{a}}\mathrm{=}\mathrm{(}{\mathrm{F}}^{\mathrm{T}}{\mathit{C}}_{\mathrm{n}}^{-1}\mathrm{ZF}{\mathrm{)}}^{-1}{\mathrm{F}}^{\mathrm{T}}{\mathit{C}}_{\mathrm{n}}^{-1}\mathrm{Z}{V}\mathit{,}$(29)where $\mathrm{Z}\mathrm{=}\mathrm{I}\mathrm{-}\begin{array}{c}\mathrm{˜}\\ {P}\end{array}\mathrm{(}{\mathrm{P}}^{\mathrm{T}}{\mathit{C}}_{\mathrm{n}}^{-1}\begin{array}{c}\mathrm{˜}\\ {P}\end{array}{\mathrm{)}}^{-1}{\begin{array}{c}\mathrm{˜}\\ {P}\end{array}}^{\mathrm{T}}{\mathit{C}}_{\mathrm{n}}^{-1}\mathit{.}$(30)We use a hat to indicate that $\hat{{a}}$ is an estimate of the true a. We are here making use of the sparse structure of the pointing matrix, which allows us to invert matrix ${\begin{array}{c}{}_{\mathrm{˜}}\\ {P}\end{array}}^{\mathrm{T}}{\mathit{C}}_{\mathrm{n}}^{-1}\begin{array}{c}{}_{\mathrm{˜}}\\ {P}\end{array}$ through non-iterative methods. For a detailed solution of an equivalent problem in mapmaking, see Keihänen et al. (2010) and references therein. The linear system in Eq. (29)is much smaller than the original one. The rank of the system is similar to the number of pointing periods, which is 44 070 for the full four-year mission.

Equation (29)can be solved by conjugate gradient iteration. The map correction is obtained as $\hat{{m}}\mathrm{=}\mathrm{\left(}{\begin{array}{c}\mathrm{˜}\\ {P}\end{array}}^{\mathit{T}}{\mathit{C}}_{\mathit{n}}^{-1}\begin{array}{c}\mathrm{˜}\\ {P}\end{array}{\mathrm{)}}^{-1}\stackrel{}{{{{P}}^{\mathit{T}}}_{\mathrm{˜}}}{\mathit{C}}_{\mathit{n}}^{-1}\mathrm{\right(}{V}\mathrm{-}{F}\hat{{a}}\mathrm{)}\mathit{.}$(31)Matrix ${\begin{array}{c}{}_{\mathrm{˜}}\\ {P}\end{array}}^{\mathit{T}}{\mathit{C}}_{\mathit{n}}^{-1}\begin{array}{c}{}_{\mathrm{˜}}\\ {P}\end{array}$ is diagonal, and inverting it is a trivial task.

A lower limit for the gain uncertainty, based on radiometer white noise alone, is given by the covariance matrix ${\mathit{C}}_{\mathit{â}}\mathrm{=}\mathrm{(}{{F}}^{\mathit{T}}{\mathit{C}}_{\mathit{n}}^{-1}{F}{\mathrm{)}}^{-1}\mathit{.}$(32)

4.2. Constrained algorithm

4.2.1. Role of the solar dipole

The dipole signal is a sum of the solar and orbital contributions. The solar dipole can be thought of as being picked from an approximately¹¹ constant dipole map, while the orbital component depends on beam orientation and satellite velocity. The latter can be used as an independent and absolute calibration. As we discuss in Sect. 5, this has allowed us to determine the amplitude and direction of the solar dipole and decouple the Planck absolute calibration from that of WMAP.

The solar dipole can be interpreted either as part of the dipole signal B ∗ D or part of the sky map m. This has important consequences. The advantage is that we can calibrate using only the orbital dipole, which is better known than the solar component and can be measured absolutely. (It only depends on the temperature of the CMB monopole and the velocity of the Planck spacecraft.) When the unconstrained DaCapo algorithm is run with erroneous dipole parameters, the difference between the input dipole and the true dipole simply leaks into the sky map m. The map can then be analysed to yield an estimate for the solar dipole parameters.

The drawback from the degeneracy is that the overall gain level is weakly constrained, since it is determined from the orbital dipole alone. In the absence of the orbital component, a constant scaling factor applied to the gains would be fully compensated for by an inverse scaling applied to the signal. It would then be impossible to determine the overall scaling of the gain. The orbital dipole breaks the degeneracy, but leaves the overall gain level weakly constrained compared with the relative gain fluctuations.

The degeneracy is not perfect, since the signal seen by a radiometer is modified by the beam response B. In particular, a beam sidelobe produces a strongly orientation-dependent signal. This is, however, a small correction to the full dipole signal.

4.2.2. Dipole constraint

Because of the degeneracy between the overall gain level and the map dipole, it makes sense to constrain the map dipole to zero. For this to work, two conditions must be fulfilled: 1) the solar dipole must be known; and 2) the contribution of foregrounds (outside the mask) to the dipole of the sky must either be negligible, or it must be known and included in the dipole model.

In the following we assume that both the orbital and the solar dipole are known. We aim at deriving a modified version of the DaCapo algorithm, where we impose the additional constraint ${{m}}_{\mathit{D}}^{\mathit{T}}{m}\mathrm{=}\mathrm{0}$. Here m_D is a a map representing the solar dipole component. We are thus requiring that the dipole in the direction of the solar dipole is completely included in the dipole model D, and nothing is left for the sky map. We note that m_D only includes the pixels outside the mask.

It turns out that condition ${{m}}_{\mathit{D}}^{\mathit{T}}{m}\mathrm{=}\mathrm{0}$ alone is not sufficient, since there is another degeneracy in the model that must be taken into account. The monopole of the sky map is not constrained by data, since it cannot be distinguished from a global noise offset b_k = const. It is therefore possible to satisfy the condition ${{m}}_{\mathit{D}}^{\mathit{T}}{m}\mathrm{=}\mathrm{0}$ by adjusting the baselines and the monopole of the map simultaneously, with no cost in χ². To avoid this pitfall, we simultaneously constrain the dipole and the monopole of the map. We require ${{m}}_{\mathit{D}}^{\mathit{T}}{m}\mathrm{=}\mathrm{0}$ and 1^Tm = 0, and combine them into one constraint ${{m}}_{\mathit{c}}^{\mathit{T}}{m}\mathrm{=}\mathrm{0}\mathit{,}$(33)where m_c now is a two-column object.

We add an additional prior term to Eq. (28)$\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\mathrm{(}{V}\mathrm{-}\begin{array}{c}\mathrm{˜}\\ {P}\end{array}\begin{array}{c}\mathrm{˜}\\ {m}\end{array}\mathrm{-}{F}{a}{\mathrm{)}}^{\mathit{T}}{\mathit{C}}_{\mathit{n}}^{-1}\mathrm{(}{V}\mathrm{-}\begin{array}{c}\mathrm{˜}\\ {P}\end{array}\begin{array}{c}\mathrm{˜}\\ {m}\end{array}\mathrm{-}{F}{a}\mathrm{)}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array} \mathrm{+}{\begin{array}{c}\mathrm{˜}\\ {m}\end{array}}^{\mathit{T}}{{m}}_{\mathit{c}}{\mathit{C}}_{\mathit{d}}^{-1}{{m}}_{\mathit{c}}^{\mathit{T}}\begin{array}{c}\mathrm{˜}\\ {m}\end{array}\end{array}& & \end{array}$(34)and aim at taking ${\mathit{C}}_{\mathit{D}}^{-1}$ to infinity. This will drive ${{m}}_{\mathit{c}}^{\mathit{T}}{m}$ to zero.

Minimization of Eq. (34) yields the solution $\hat{{a}}\mathrm{=}\mathrm{(}{\mathit{F}}^{\mathit{T}}{\mathit{C}}_{\mathit{n}}^{-1}\mathit{ZF}{\mathrm{)}}^{-1}{\mathit{F}}^{\mathit{T}}{\mathit{C}}_{\mathit{n}}^{-1}\mathit{Z}{V}\mathit{,}$(35)with $\mathit{Z}\mathrm{=}\mathit{I}\mathrm{-}\begin{array}{c}\mathrm{˜}\\ {P}\end{array}\mathrm{(}\mathit{M}\mathrm{+}{{m}}_{\mathit{c}}{\mathit{C}}_{\mathit{D}}^{-1}{{m}}_{\mathit{c}}^{\mathit{T}}{\mathrm{)}}^{-1}\stackrel{}{{{{P}}^{\mathit{T}}}_{\mathrm{˜}}}{\mathit{C}}_{\mathit{n}}^{-1}\mathit{,}$(36)where for brevity we have written $\mathit{M}\mathrm{=}{\begin{array}{c}\mathrm{˜}\\ {P}\end{array}}^{\mathit{T}}{\mathit{C}}_{\mathit{n}}^{-1}\begin{array}{c}\mathrm{˜}\\ {P}\end{array}\mathit{.}$(37)This differs from the original solution (Eqs. (29), (30)) by the term ${{m}}_{\mathit{c}}{\mathit{C}}_{\mathit{D}}^{-1}{{m}}_{\mathit{c}}^{\mathit{T}}$ in the definition of Z.

Equation (36)is impractical owing to the large size of the matrix to be inverted. To proceed, we apply the Sherman-Morrison formula and let C_D → 0, yielding $\mathrm{(}\mathit{M}\mathrm{+}{{m}}_{\mathit{c}}{\mathit{C}}_{\mathit{D}}^{-1}{{m}}_{\mathit{c}}^{\mathit{T}}{\mathrm{)}}^{-1}\mathrm{=}{\mathit{M}}^{-1}\mathrm{-}{\mathit{M}}^{-1}{{m}}_{\mathit{c}}\mathrm{(}{{m}}_{\mathit{c}}^{\mathit{T}}{\mathit{M}}^{-1}{{m}}_{\mathit{c}}{\mathrm{)}}^{-1}{{m}}_{\mathit{c}}^{\mathit{T}}{\mathit{M}}^{-1}\mathit{.}$(38)The middle matrix ${{m}}_{\mathit{c}}^{\mathit{T}}{\mathit{M}}^{-1}{{m}}_{\mathit{c}}$ is a 2x2 block diagonal matrix, and is easy to invert.

Equations (35)−(38)are the basis of the constrained DaCapo algorithm. The system is solved using a conjugate-gradient method, similar to the unconstrained algorithm.

The map correction becomes $\hat{{m}}\mathrm{=}\mathrm{(}\mathit{M}\mathrm{+}{{m}}_{\mathit{c}}{\mathit{C}}_{\mathit{D}}^{-1}{{m}}_{\mathit{c}}^{\mathit{T}}{\mathrm{)}}^{-1}{\begin{array}{c}\mathrm{˜}\\ {P}\end{array}}^{\mathit{T}}{\mathit{C}}_{\mathit{n}}^{-1}\mathrm{(}{V}\mathrm{-}\mathit{F}\hat{{a}}\mathrm{)}\mathit{.}$(39)One readily sees that $\hat{{m}}$ fulfils the condition expressed by Eq. (33), and thus so does the full map m.

The constraint breaks the degeneracy between the gain and the signal, but also makes the gains again dependent on the solar dipole, which must be known beforehand.

4.3. Use of unconstrained and constrained algorithms

We have used the unconstrained and constrained versions of the algorithm together to obtain a self-consistent calibration and to obtain an independent estimate for the Solar dipole.

We first ran the unconstrained algorithm, using the known orbital dipole and an initial guess for the solar dipole. The results depend on the solar dipole only through the beam correction. The difference between the input dipole and the true dipole are absorbed in the sky map.

We estimated the solar dipole from these maps (Sect. 5); since the solar dipole is the same for all radiometers, we combined data from all the 70 GHz radiometers to reduce the error bars. (Simply running the unconstrained algorithm, fixing the dipole, and running the constrained version with same combination of radiometers would have just yielded the same solution.)

Once we had produced an estimate of the solar dipole, we reran DaCapo in constrained mode to determine the calibration coefficients K more accurately.

5. Characterization of the orbital and solar dipoles

In this section, we explain in detail how the solar dipole was obtained for use in the final DaCapo run mentioned above. We also compare the LFI measurements with the Planck nominal dipole parameters, and with the WMAP values given by Hinshaw et al. (2009).

5.1. Analysis

When running DaCapo in constrained mode to compute the calibration constants K (Eq. (1)), the code needs an estimate of the solar dipole in order to calibrate the data measured by the radiometers (see Sect. 3 and especially Fig. 3), since the signal produced by the orbital dipole is ten times weaker. We used DaCapo to produce this estimate from the signal produced by the orbital dipole. We limited our analysis to the 70 GHz radiometric data, since this is the cleanest frequency in terms of foregrounds. The pipeline was provided by a self-contained version of the DaCapo program, run in unconstrained mode (see Sect. 4.2), in order to make the orbital dipole the only source of calibration, while the solar dipole is left in the residual sky map.

We bin the uncalibrated differenced time-ordered data into separated rings, with one ring per pointing period. These data are then binned according to the direction and orientation of the beam, using a Healpix¹² (Górski et al. 2005) map of resolution N_side = 1024 and 256 discrete bins for the orientation angle ψ. The far sidelobes should prevent a clean dipole from being reconstructed in the sky model map, since the signal in the timeline is convolved with the beam B over the full sphere. To avoid this, an estimate for the pure dipole is obtained by subtracting the contribution due to far sidelobes using an initial estimate of the solar dipole, which in this case was the WMAP dipole (Hinshaw et al. 2009). We also subtract the orbital dipole at this stage. The bias introduced by using a different dipole is of second order and is discussed in the section on error estimates (5.2).

DaCapo builds a model sky brightness distribution that is used to clean out the polarized component of the CMB and foreground signals to only leave noise, the orbital dipole, and far-sidelobe pickup. This sky map is assumed to be unpolarized, but since radiometers respond to a single linear polarization, the data will contain a polarized component, which is not compatible with the sky model and thus leads to a bias in the calibration. An estimate of the polarized signal, mostly CMB E modes and some synchrotron, needs to be subtracted from the timelines. To bootstrap the process we need an intial gain estimate, which is provided by DaCapo constrained to use the WMAP dipole. We then used the inverse of these gains to convert calibrated polarization maps from the previous LFI data release (which also used WMAP dipole calibration) into voltages and unwrap the map data into the timelines using the pointing information for position and boresight rotation. This polarized component due to E modes, which is ~ 3.5 μK rms on small angular scales plus an additional large-scale contribution of the North Galactic Spur of amplitude ~ 3 μK, was then subtracted from the time-ordered data. Further iterations using the cleaned timelines were found to make a negligible difference.

5.2. Results

To make maps of the dipole, a second DaCapo run was made in the unconstrained mode for each LFI 70 GHz detector using the polarization-cleaned timelines and the 30 GHz Madam mask, which allows 78% of the sky to be used. The extraction of the dipole parameters (Galactic latitude, longitude, and temperature amplitude) in the presence of foregrounds was achieved with a simple Markov chain Monte Carlo (MCMC) template-fitting scheme. Single detector hit maps, together with the white noise in the LFI-reduced instrument model (RIMO), were used to create the variance maps needed to construct the likelihood estimator for the MCMC samples. Commander maps (Planck Collaboration X 2016) were used for synchrotron, free-free, and thermal dust for the template maps with the MCMC fitting for the best amplitude scaling factor to clean the dipole maps. The marginalized distribution of the sample chains between the 16th and 84th percentiles were used to estimate the statistical errors, which were 0.̊004, 0.̊009, and 0.16μK for latitude, longitude, and amplitude respectively. The 50% point was taken as the best parameter value, as shown in Table 1. To estimate the systematic errors on the amplitude in calibration process due to white noise, 1 /f noise, gain fluctuations, and ADC corrections, simulated time-ordered data were generated with these systematic effects included. These simulated timelines were then calibrated by DaCapo in the same way as the data. The standard deviation of the input to output gains were taken as the error in absolute calibration with an average value of 0.11%.

Plots of the dipole amplitudes with these errors are shown in Fig. 6, together with the error ellipses for the dipole direction. As can be seen, the scatter is greater than the statistical error. Therefore, we take a conservative limit by marginalizing over all the MCMC samples for all the detectors, which results in an error ellipse (±0.̊02, ±0.̊05) centred on Galactic latitude and longitude (48.̊26, 264.̊01). The dipole amplitudes exhibit a trend in focal plane position, which is likely due to residual, unaccounted-for power in far sidelobes, which would be symmetrical between horn pairs. These residuals, interacting with the solar dipole, would behave like an orbital dipole, but in opposite ways on either side of the focal plane. This residual therefore cancels out to first order in each pair of symmetric horns in the focal plane, i.e., horns 18 with 23, 19 with 22, and 20 with 21 (see inset in Fig. 6). Since the overall dipole at 70 GHz is calculated by combining all the horns, the residual effect of far sidelobes is reduced.

6. Validation of the calibration and accuracy assessment

In this section we present the results of a set of checks we have run on the data that comprise this new Planck release. Table 2 quantifies the uncertainties that affect the calibration of the LFI radiometers.

6.1. Absolute calibration

In this section we provide an assessment of the change in the absolute level of the calibration since the first Planck data release, in terms of its impact on the maps and power spectra. Generally speaking, a change in the average value of G in Eq. (1)of the form $⟨\mathit{G}⟩\mathrm{\to }⟨\mathit{G}⟩\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\delta }}_{\mathit{G}}\mathrm{)}\mathit{,}\mathrm{with}{\mathit{\delta }}_{\mathit{G}}\mathrm{\ll }\mathrm{1}$(40)leads to a change of ⟨T⟩ → ⟨T⟩(1 + δ_G) in the average value of the pixel temperature T, and to a change C_ℓ → C_ℓ(1 + 2δ_G) in the average level of the measured power spectrum C_ℓ, before the application of any window function. Our aim is to quantify the value of the variation δ_G from the previous Planck-LFI data release to the current one. We did this by comparing the level of power spectra in the ℓ = 100–250 multipole range consistently with Cal13.

There have been several improvements in the calibration pipeline that have led to a change in the value of ⟨G⟩:

• 1.

  The peak-to-peak temperature difference of the referencedipole D used in Eq. (1)has changed by + 0.27% (see Sect. 5), because we now use the solar dipole parameters calculated from our own Planck measurements (Planck Collaboration I 2016).

• 2.

  In the same equation, we no longer convolve the dipole D with a beam B that is a delta function, but instead use the full profile of the beam over the sphere (see Sect. 2).

• 3.

  The beam normalization has changed, since in this data release B is such that (Planck Collaboration IV 2016) ${\mathrm{\int }}_{\mathrm{4}\mathit{\pi }}\mathit{B}\mathrm{\left(}\mathit{\theta ,\varphi }\mathrm{\right)}\hspace{0.17em}\mathrm{d\Omega }\lneqq \mathrm{1.}$(41)

• 4.

  The old dipole fitting code has been replaced with a more robust algorithm, DaCapo (see Sect. 4).

Table 4 lists the impact of these effects on the amplitude of fluctuations in the temperature ⟨T⟩ of the 22 LFI radiometer maps. The numbers in this table have been computed by rerunning the calibration pipeline on all the 22 LFI radiometer data with the following setup:

• 1.

  A pencil-beam approximation for B in Eq. (1)has been used, instead of the full 4π convolution (“Beam convolution” column), with the impact of this change quantified by Eq. (16), and the comparison between the values predicted by this equation with the measured change in the a_ℓm harmonic coefficients shown in Fig. 8 (the agreement is excellent, better than 0.03%).

• 2.

  The old calibration code has been used instead of the DaCapo algorithm described in Sect. 4 (“Pipeline upgrades” column).

• 3.

  The B ∗ T_Gal term has not been removed, as in the discussion surrounding Eq. (11)(“Galactic sidelobe removal” column).

• 4.

  The signal D used in Eq. (1)has been modelled using the dipole parameters published in Hinshaw et al. (2009), as was done in Cal13 (“Reference dipole” column).

We measured the actual change in the absolute calibration level by considering the radiometric maps (i.e., maps produced using data from one radiometer) of this data release (indicated with a prime) and of the previous data release and averaging the ratio: ${\mathrm{\Delta }}_{\mathit{\ell }}^{\mathit{x,}{\mathit{x}}^{\mathrm{\prime }}}\mathrm{=}{⟨\frac{{\mathit{C}}_{{\mathit{x}}^{\mathrm{\prime }}\mathrm{\times }\mathit{y}}}{{\mathit{C}}_{\mathit{x}\mathrm{\times }\mathit{y}}}⟩}_{\mathit{y}}\mathrm{-}\mathrm{1}\mathit{,}$(42)where y indicates a radiometer at the same frequency as x and x′, such that y ≠ x. The average is meant to be taken over all the possible choices for y (thus 11 choices for 70 GHz radiometers, 5 for 44 GHz, and 3 for 30 GHz) in the multipole range 100 ≤ ℓ ≤ 250. The way that cross-spectra are used in Eq. (42)ensures that the result does not depend on the white noise level. Of course, this result quantifies the ratio between the temperature fluctuations (more correctly, between the a_ℓm coefficients of the expansion of the temperature map in spherical harmonics) in the two data releases, and not between the power spectra.

The column labelled “Estimated change” in Table 4 contains a simple combination of allthe numbers in the table: $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\mathrm{Estimatedchange}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{ϵ}}_{\mathrm{beam}}\mathrm{)}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{ϵ}}_{\mathrm{pipeline}}\mathrm{)}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{\times }\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{ϵ}}_{\mathrm{Gal}}\mathrm{)}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{ϵ}}_{\mathit{D}}\mathrm{)}\mathrm{-}\mathrm{1}\mathit{,}\end{array}& & \end{array}$(43)where the ϵ factors are the numbers shown in the same table and discussed above. This formula assumes that all the effects are mutually independent. This is, of course, an approximation; however, the comparison between this estimate and the measured value (obtained by applying Eq. (42) to the 2013 and 2015 release maps) can give an idea of the amount of interplay of these effects in producing the observed shift in temperature. Figure 9 shows a plot of the contributions discussed above, as well as a visual comparison between the measured change in the temperature and the estimate from Eq. (43).

6.2. Noise in dipole fitting

We performed a number of simulations that quantify the impact of white noise in the data on the estimation of the calibration constant, as well as the ability of our calibration code to retrieve the true value of the calibration constants. Details of this analysis are described in Planck Collaboration III (2016). We did not include such sorts of errors as an additional element in Table 2, because the statistical error is already included in the row “Spread among independent radiometers”.

6.3. Beam uncertainties

As discussed in Planck Collaboration IV (2016), the beams B used in the LFI pipeline are very similar to those presented in Planck Collaboration IV (2014); they are computed with GRASP, properly smeared to take the satellite motion into account. Simulations were performed using the optical model described in Planck Collaboration IV (2014), which was derived from the Planck radio frequency Flight Model (Tauber et al. 2010) by varying some optical parameters within the nominal tolerances expected from the thermoelastic model, in order to reproduce the measurements of the LFI main beams from seven Jupiter transits. This is the same procedure adopted in the 2013 release (Planck Collaboration IV 2014); however, unlike the case presented in Planck Collaboration IV (2014), a different beam normalization is introduced here to properly take the actual power entering the main beam into account (typically about 99% of the total power). This is discussed in more detail in Planck Collaboration IV (2016).

Given the broad use of beam shapes B in the current LFI calibration pipeline, it is extremely important to assess their accuracy and the way errors in B propagate down to the estimate of the calibration constants K in Eq. (1).

In the previous data release we did not use our knowledge of the bandpasses of each radiometer to produce an in-band model of the beam shape, but instead estimated B by means of a monochromatic approximation (see Cal13). In that case, we estimated the error induced in the calibration as the variation of the dipole signal when using either a monochromatic or a band-integrated beam, since we believe the latter to be a more realistic model.

In this data release, we have switched to the full bandpass-integrated beams produced using GRASP, which represents our best knowledge of the beam (Planck Collaboration IV 2016). We tested the ability of DaCapo to retrieve the correct calibration constants K for LFI19M (a 70 GHz radiometer) when the large-scale component (ℓ = 1) of the beam’s sidelobes is: (1) rotated arbitrarily by an angle −160° ≤ θ ≤ 160°; or (2) scaled by ± 20%. We find that such variations alter the calibration constants by approximately 0.1%. However, we do not list such a small number as an additional source of uncertainty in Table 2, since we believe that this is already captured by the scatter in the points shown in Fig. 10, which were used to produce the numbers in the row Inconsistencies among radiometers.

6.4. Inter-channel calibration consistency

In this section we provide a quantitative estimate of the relative calibration error for the LFI frequency maps by measuring the consistency of the power spectra computed using data from one radiometer at time. By relative error we mean any error that is different among the radiometers, in contrast to an absolute error, which induces a common shift in the power spectrum. We computed the power spectrum of single radiometer half-ring maps and have estimated the variation in the region around the first peak (100 ≤ ℓ ≤ 250), since this is the multipole range with the best S/N.

The result of this analysis is shown in Fig. 10, which plots the values of the quantity ${\mathit{\delta }}_{\mathrm{rad}}\mathrm{=}\frac{{\mathrm{⟨}{{\mathit{C}}_{\mathit{\ell }}^{\mathrm{HR}}}^{\mathrm{⟩}}}_{\mathrm{rad}}}{{\mathrm{⟨}{{\mathit{C}}_{\mathit{\ell }}^{\mathrm{HR}}}^{\mathrm{⟩}}}_{\mathrm{freq}}}\mathrm{-}\mathrm{1}\mathit{,}$(44)where ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{HR}}$ is the cross-power spectrum computed using two half-ring maps, and ⟨·⟩ denotes an average over ℓ. The quantity ⟨C_ℓ⟩_freq is the same average computed using the full frequency half-ring maps. The δ_rad slope is symmetric around zero in the 70 GHz radiometers; this might be caused by residual unaccounted-for power in the far sidelobes of the beam. The same explanation was advanced in Sect. 5 to explain a similar effect. It is interesting to note that the amplitude of the two systematics is comparable; the trend in Fig. 6 has a peak-to-peak variation (in temperature) of about 0.5%, while the trend in Fig. 10 has a variation (in power) of roughly 1.0%. We combine the values of δ_rad for those pairs of radiometers whose beam position in the focal plane is symmetric (e.g., 18M versus 23M, 18S versus 23S, 19M versus 22M, etc.), since in these pairs the unaccounted-for power should be balanced. We have found that indeed all the six combinations of δ_rad are consistent with zero within 1 σ (see the inset of Fig. 10).

Since the cross-spectrum of two half-ring maps does not depend on the level of uncorrelated noise, the fluctuations of δ_i around the average value that can be seen in Fig. 10 can be interpreted as relative calibration errors. If we limit our analysis to the multipole range 100 ≤ ℓ ≤ 250, we can estimate the error of the 70 GHz map as the error on the average height of the peaks (i.e., the value $\mathit{\sigma }\mathit{/}\sqrt{\mathit{N}}$, with σ being the standard deviation and N the number of points) that is, 0.25, 0.16, and 0.10 percent and 30, 44, and 70 GHz, respectively.

6.5. Inter-frequency calibration consistency

In this section we carry out an analysis similar to the one presented in Sect. 6.4, where we compare the absolute level of the maps at the three LFI frequencies, i.e., 30, 44, and 70 GHz. We make use of the full frequency maps, as well as the pair of half-ring maps at 70 GHz. Each half-ring map has been produced using data from one of the two halves of each pointing period. We quantify the discrepancy between the 70 GHz map and another map by means of the quantity ${\mathrm{\Delta }}_{\mathit{\ell }}^{\mathrm{70}\mathrm{GHz,other}}\mathrm{=}\frac{{\mathit{C}}_{\mathit{\ell }}^{\mathrm{HR}\mathrm{1}\mathrm{\times }\mathrm{HR}\mathrm{2}}}{{\mathit{C}}_{\mathit{\ell }}^{\mathrm{HR}\mathrm{1}\mathrm{\times }\mathrm{other}}}\mathrm{-}\mathrm{1}\mathit{,}$(45)where ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{HR}\mathrm{1}\mathrm{\times }\mathrm{HR}\mathrm{2}}$ is the cross-spectrum between the two 70 GHz half-ring maps, and ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{HR}\mathrm{1}\mathrm{\times }\mathrm{other}}$ is the cross-spectrum between the first 70 GHz half-ring map and the map under analysis. In the ideal case (perfect correspondence between the spectrum of the 70 GHz map and the other map) we expect Δ_ℓ = 0. As was the case for Eq. (44), this formula has the advantage of discarding the white noise level of the spectrum ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{other}}$ by using the cross-spectrum with the 70 GHz map, whose noise should be uncorrelated.

Over the multipole range 100 ≤ ℓ ≤ 250, the average discrepancy¹³ is 0.15 ± 0.17% for the 44 GHz map, and 0.15 ± 0.26% for the 30 GHz map, as shown in Fig. 11. Such numbers are consistent with the calibration errors provided in Sect. 6.4.

6.6. Null tests

In Cal13 we provided a study of a number of null tests, with the purpose of testing the quality of the calibration. In this new data release, we have moved the bulk of the discussion to Planck Collaboration III (2016). We just show one example here, which is particularly relevant in the context of the LFI calibration validation. Figure 12 shows the variation in the quality of the maps due to the use of the full 4π convolution versus a pencil beam approximation, as discussed in Sect. 2. The analysis of many similar difference maps has provided us with sufficient evidence that using the full 4π beam convolution reduces the level of systematic effects in the LFI maps.

7. Measuring the brightness temperature of Jupiter

The analysis of the flux densities of planets for this Planck data release has been considerably extended. We only use Jupiter data for planet calibration, so we now focus the discussion on observations of this planet. The new analysis includes all seven transits of Jupiter through each main beam of the 22 LFI radiometers. The analysis pipeline has been improved considerably by considering several effects not included in Cal13.

Planets provide a useful calibration cross-check; in particular, the measurement of the brightness temperature of Jupiter can be a good way to assess the accuracy of the calibration, since Jupiter is a remarkably bright source with a S/N per scan as high as 50 and a relatively well known spectrum. Furthermore, at the resolution of LFI beams, it can be considered a point-like source.

7.1. Input data

Table 5 lists the epochs when the LFI main beams crossed Jupiter, and Figs. 13 and 14 give a visual timeline of these events. The first four transits occurred in nominal scan mode (spin shift 2′, 1° per day) with a phase angle of 340°, and the last three scans in deep mode (shift of the spin axis between rings of 0.5′, 15′ per day) with a phase angle of 250° (see Planck Collaboration I 2016). The analysis follows the procedure outlined in Cal13, but with a number of improvements:

• 1.

  The brightness of Jupiter was extracted from timelines to fullyexploit the time dependence in the data.

• 2.

  Seven transits have been considered instead of two, which allowed us to analyse the sources of scatter better among the measurements.

• 3.

  All the data were calibrated simultaneously.

• 4.

  Different extraction methods were exploited to find the most reliable among them.

In the following discussion, we refer to a timeline (one for each of the 22 LFI radiometers) as the list of values (t,x_p,t,x_t,ψ_t,ΔT_ant,t) with t the epoch of observation, x_p,t the instantaneous apparent planet positions as seen from Planck, x_t and ψ_t the corresponding beam pointing directions and orientations, and ΔT_ant,t the measured antenna temperature. We took the the values of x_t and ψ_t, as well as the calibrated values of ΔT_ant,t cleaned from the dipole and quadrupole signals, from the output of the LFI pipeline. We recovered x_p,t from the Horizons¹⁴ on-line service.

Samples from each radiometer timeline were used in this analysis only if the following conditions were met: (1) the samples were acquired in stable conditions during a pointing period (Planck Collaboration II 2016); (2) the pipeline did not flag them as “bad”; (3) their angular distance from the planet position at the time of the measurement was less than 5°; and (4) they were not affected by any anomaly or relevant background source. We checked the last condition by visually inspecting small coadded maps of the selected samples.

7.2. Description of the analysis pipeline

In the following paragraphs, we describe how we improved the pipeline used to extract the brightness temperature T_B of Jupiter from the raw LFI data. This extraction goes through an initial estimation of the antenna temperature T_A and a number of corrections to take various systematic effects into account. We present the methods used to estimate T_A in Sect. 7.2.1, and then in Sect. 7.2.2 we discuss the estimation of T_B. Since the computation of T_B requires an accurate estimate of the planet solid angle Ω_p, we discuss the computation of this factor in a dedicated part, Sect. 7.2.3.

7.2.1. Estimation of the antenna temperature

Following Cal13 and Cremonese et al. (2002), the recovery of the instantaneous planet signal from a timeline is equivalent to the deconvolution of the planet shape from the beam pattern B_t at time t. Since the planet can be considered a point source, the most practical way is to assume $\mathrm{\Delta }{\mathit{T}}_{\mathrm{ant}\mathit{,t}}\mathrm{=}{\mathit{T}}_{\mathrm{A}\mathit{,}\mathrm{p}}{\mathit{B}}_{\mathit{t}}\mathrm{\left(}\mathit{\delta }{{x}}_{\mathrm{p}\mathit{,t}}\mathrm{\right)}\mathrm{+}\mathit{b,}$(46)where T_A,p is the unknown planet antenna temperature, b the background, and B_t(δx_p,t) the beam response for the planet at the time of observation. Of course, B_t depends on the relative position of the planet with respect to the beam, δx_p,t. If a suitable beam model is available, B_t can be determined and T_A,p can be recovered from least squares minimization. We use an elliptic Gaussian centred on the instantaneous pointing direction as a model for the beam, because it shows a very good match with the main beam of the GRASP model (Planck Collaboration IV 2016), with peak-to-peak discrepancies of a few tenths of a percent (the importance of far sidelobes is negligible for a source as strong as Jupiter). To compute B_t, the pointings are rotated into the beam reference frame, since this allows for better control of the beam pattern reconstruction¹⁵.

7.2.2. Estimation of the brightness temperature

In Cal13, we computed the brightness temperature T_B from the antenna temperature T_A by means of the following formula (assuming monochromatic radiometers): ${\mathit{T}}_{\mathrm{B}}\mathrm{=}{\mathit{B}}_{\mathrm{Planck}}^{-1}\left({\mathit{T}}_{\mathrm{A}}{\mathit{f}}_{\mathrm{sl}}\frac{{\mathrm{\Omega }}_{\mathrm{b}}}{{\mathrm{\Omega }}_{\mathrm{p}}}{\frac{\mathit{\partial }{\mathit{B}}_{\mathrm{Planck}}}{\mathit{\partial T}}|}_{{\mathit{T}}_{\mathrm{CMB}}}\right)\mathit{,}$(47)where B_Planck is Planck’s blackbody function, Ω_b and Ω_p the beam and planet solid angles, and T_CMB = 2.7255 K is the temperature of the CMB monopole. In this 2015 Planck data release, we have also introduced corrections to account for the bandpass.

The accuracy of the T_B determination is affected by confusion noise (i.e., noise caused by other structures in the maps), which we estimate from the standard deviation of samples taken between a radius of 1° and 1.5° (depending on the beam) and 5° from the beam centre. These samples are masked for strong sources or other defects. Since background maps are not subtracted, the confusion noise is greater than the pure instrumental noise. However, since the histogram is described well by a normal distribution, we used error propagation¹⁶ to assess the accuracy of T_B against the confusion noise.

In the conversion of T_A,p into T_B through Eq. (46), the pipeline implements a number of small corrections:

• 1.

  detector-to-detector differences in the beam solid angleΩ_b, accounting for ± 6%, which is probably the most important effect;

• 2.

  changes in the solid angle of the planet, Ω_p, due to the change of the Jupiter–Planck distance, which introduces a correction factor of up to 6.9% percent;

• 3.

  changes in the projected planet ellipticity, due to the planeto-centric latitude of the observer and the oblateness of the planet, to reduce observations as if they were made at Jupiter’s pole;

• 4.

  blocking of background radiation by the planet, changing from about 0.7% to 1.5%, depending on the ratio Ω_p/ Ω_b;

• 5.

  a φ_sl correction, which accounts for the fraction of radiation not included in the main beam (about 0.2%).

7.2.3. Determination of the solid angle of Jupiter

The solid angle Ω_p of Jupiter for a given planet-spacecraft distance Δ_p and planeto-centric latitude δ_P is given by ${\mathrm{\Omega }}_{\mathrm{p}}\mathrm{\left(}{\mathit{\delta }}_{\mathrm{P}}\mathrm{\right)}\mathrm{=}{\mathrm{\Omega }}_{\mathrm{p}}^{\mathrm{polar}\mathit{,}\mathrm{ref}}{\left(\frac{{\mathrm{\Delta }}_{\mathrm{p}\mathit{,}\mathrm{ref}}}{{\mathrm{\Delta }}_{\mathrm{p}}}\right)}^{\mathrm{2}}{\mathit{e}}_{\mathrm{p}}\sqrt{\mathrm{1}\mathrm{-}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{p}}^{\mathrm{2}}\mathrm{)}{\sin }^{\mathrm{2}}{\mathit{\delta }}_{\mathrm{P}}}\mathit{,}$(48)where Δ_p,ref is a fiducial planet-spacecraft distance (for Jupiter, Δ_p,ref ≈ 5.2 AU), ${\mathrm{\Omega }}_{\mathrm{p}}^{\mathrm{polar}\mathit{,}\mathrm{ref}}$ is the solid angle of the planet as seen from its pole at the fiducial distance, and e_p(< 1) the ratio between the polar and equatorial radii of the planet.

Our new pipeline does not use the ellipticity, FWHM, and orientation of the elliptical beam parameters provided by Planck Collaboration IV (2016), as they were derived from a marginalization over T_A,p on the same Jupiter data used in this analysis. However, because of a degeneracy between T_A,p and the beam parameters, the results were sensitive to details of the fitting procedure (up to about 1%), such as the radius of the area being analysed and the minimization method. Therefore, we have first determined a new set of beam parameters for each transit and computed the weighted average (using a numerical minimization), and then we used these parameters to determine T_A,p.

7.3. Results

We present here the results per radiometer and transit, and we also discuss how to combine the measurements performed using the 30, 44, and 70 GHz radiometers to obtain three estimates of T_B at these three nominal frequencies.

7.3.1. Brightness temperatures per radiometer and transit

In Cal13, based on the first two transits, we noted that confusion noise alone was not able to account for variations in T_B found for different radiometers belonging to the same frequency channel. We therefore assumed the presence of some unidentified systematics error, dominating the ultimate accuracy in the measures.

That residual systematics are more important than confusion noise and background is still true in the 2015 data release. The values of T_B measured by radiometers in the same frequency channel have a spread of 0.6%, 1.0%, and 0.6% of the average signal at 30, 44 and 70 GHz, respectively, and are not normally distributed (see Fig 15). These observed dispersions are a factor of ~ 3 larger than the confusion noise, and cannot be ascribed to the background, whose effect is only 1% of the observed scatter. The excess dispersion must be due to a small residual systematic effect such as pointing, beam model, or mismatch in the centre frequency.

One possible cause of the observed dispersion in the brightness temperature is some systematic effect in the estimation of the beam parameters (see Sect. 7.2.3).

Non-Gaussianities in the beam, as well as beam smearing were investigated by replacing the elliptical beam with band-averaged beam maps derived from GRASP calculations. The results are consistent with the elliptical beam, with residuals of at most 4 × 10^-3 K in T_A,p.

We computed T_A,p and T_B again using an analytical approximation based on the assumption of negligible background (which is quite a good approximation for Jupiter), and compared the results. The two methods agree at the level of 6 × 10^-7 K for T_A,p, and at the level of 6 × 10^-4 K for T_B.

7.3.2. Combination of the results and comparison with WMAP

We now want to combine the measurements of the 22 LFI radiometers in order to have three estimates of T_B at the three LFI nominal frequencies, 30, 44, and 70 GHz.

Such a determination of T_B depends on proper knowledge of the central frequency, ν_cen, for each detector. This parameter is derived from the bandpasses, but these are not known exactly. It is possible to remove most of the differences amongst the 30 GHz and amongst the 44 GHz radiometers by changing the ν_cen values of the radiometers by as little as ± 0.2 GHz. However, it is still not clear how the bandshapes would have to be modified to explain these changes in ν_cen. For this reason, we did not include these corrections for 30 GHz and 44 GHz data.

At 70 GHz the situation is complicated even more by the fact that the ν_cen values of the radiometers are significantly spread over the channel bandwidth, so that each of the 70 GHz radiometers samples a slightly different portion of the Jupiter spectrum. Indeed, comparing T_B at 70 GHz from one transit to the next after having ordered the radiometers for increasing ν_cen, it is possible to see quite a significant correlation (see Fig. 16). Using a linear regression of the T_B at 70 GHz against ν_cen, a cleaned list of T_B values was obtained. Their average is identical to the simple weighted average of T_B for all of the transits and of the detectors T_B = (171.558 ± 0.008) K. The standard deviation of the whole set of samples reduces to just 2.7%. Interestingly, the inferred slope dT_B/ dν_cen = (0.2570 ± 0.0058) K / GHz matches the one from WMAP data very well, dT_B/ dν_cen = (0.243 ± 0.025) K / GHz, and the correlation between transits loses most of its statistical significance. We attempted the same test for the 30 GHz and 44 GHz channels, but the spread in ν_cen is too small to produce meaningful results. These results open the possibility of including the spectral slope directly as a free parameter of the fit in a future analysis.

The result of this analysis produced the brightness temperatures T_B listed in Table 6. The errors are 0.6%, 0.9%, and 0.6% at 30, 44, and 70 GHz, respectively. As in Cal13, in Fig. 17 we compared Jupiter’s T_B averaged in each LFI band with the spectrum provided by WMAP. The agreement is quite good, with a difference that does not exceed 0.5%. In this comparison, we must note that WMAP and Planck-LFI are calibrated on slightly different dipoles, with Planck assuming an amplitude of 3364.5 μK (Planck Collaboration I 2016), while WMAP assumed an amplitude of 3355 μK (Hinshaw et al. 2009). In addition, WMAP central frequencies are different from those of Planck. So to make the comparison, the measures provided in Weiland et al. (2011) are scaled by 1.00268 and linearly interpolated to the averaged ν_cen of each LFI frequency channel.

7.4. What does the analysis of Jupiter tell about the LFI calibration?

The errors in our estimates of Jupiter’s brightness temperature are 0.6%, 0.9%, and 0.6%. These are consistent with (albeit slightly larger than) the numbers presented in Sects. 6.4 and 6.5, where our analysis of the consistency of inter-channel and inter-frequency cross-spectra in the range 100 ≤ ℓ ≤ 250 produced estimates between 0.1% and 0.3%. We believe that the larger error bars are a direct consequence of a number of facts: (1) bandpass uncertainties play a larger role in the analysis of a non-thermodynamic sky signal like Jupiter’s emission in the microwave range; (2) the amount of data usable for this analysis is only a limited fraction of the overall amount used to create LFI maps; and (3) uncertainties in the main beam profile have a larger impact on the study of a point-like source such as Jupiter, with respect to the analysis presented in Sect. 6, which dealt with angular scales corresponding to 100 ≤ ℓ ≤ 250.

8. Conclusions

We have described the method used to calibrate the Planck/LFI data for the 2015 Planck data release and provided a quantitative analysis that shows the amount of change in detail at the calibration level thanks to all the improvements we implemented in the pipeline since the previous data release. Compared to our 2013 release, we have improved the LFI calibration in several ways, most notably in the use of an internally consistent dipole signal as a calibrator, as well as more accurate data analysis algorithms. As a result, we have also improved the accuracy of the LFI calibration by a factor of 3, which for the 2015 release ranges from 0.20% to 0.35%, depending on the frequency.

An important byproduct of our analysis is a novel estimate of the solar dipole signal; using LFI data alone, we have estimated the amplitude of the dipole to be 3 365.6 ± 3.0 μK and the direction of its axis to be (l,90°−b) = (264.01 ± 0.05°,48.26 ± 0.02°) (Galactic coordinates). This result matches the numbers provided by Hinshaw et al. (2009) within 1σ. This slight difference in amplitude has the effect of a shift by ~ 0.3% in the overall level of the calibrated timelines. Together with the improved LFI calibration, we now find very good agreement between the level of the spectra estimated by WMAP.

---

Acknowledgments

The Planck Collaboration acknowledges the support of: ESA; CNES, and CNRS/INSU-IN2P3-INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MINECO, JA and RES (Spain); Tekes, AoF, and CSC (Finland); DLR and MPG (Germany); CSA (Canada); DTU Space (Denmark); SER/SSO (Switzerland); RCN (Norway); SFI (Ireland); FCT/MCTES (Portugal); ERC and PRACE (EU). A description of the Planck Collaboration and a list of its members, indicating which technical or scientific activities they have been involved in, can be found at http://www.cosmos.esa.int/web/planck/planck_collaboration

References

Appendix A: Smoothing calibration curves

In this appendix we explain in some detail the algorithm used to smooth the LFI calibration constant timelines. This task is performed by the smoothing filter discussed in Sect. 3 (see also Fig. 3).

Appendix A.1: Purpose of the smoother

Variations in the orientation of the Planck spacecraft during the mission caused its instruments to observe the calibration signal (the CMB solar dipole) with varying amplitude. This in turn induced variations in the accuracy of the reconstruction of the calibration constant K = G^-1 (Eq. (1)), as shown in Fig. 4. Such variations are mainly due to statistical errors and are in general unrelated to the true stability of the LFI detectors. Therefore, before applying such calibration constants to the raw data measured by the LFI radiometers, we have applied a low-pass filter that smooths most of the high-frequency fluctuations.

In Cal13 we employed a simple smoothing filter, which used wavelets to clean the stream of calibration constants of high-frequency fluctuations (note that no smoothing was applied to the 30 GHz radiometers, since the calibration algorithm used for them made this step unnecessary).

In the 2015 Planck data release, we have improved our smoothing filter in order to consider sudden jumps in the calibration constant that are caused by a genuine change in the state of the radiometer; we call these real jumps. Such real jumps are not statistical effects and it is therefore incorrect to include them when smoothing the data. The actual smoothing algorithm used in the 2015 Planck-LFI data release works as follows:

• 1.

  determine if there are sudden variations in the calibrationconstants that might have a non-statistical origin and make a list ofthem;

• 2.

  split the stream of calibration constants into sub-streams, using the jumps found in the previous step;

• 3.

  apply a low-pass filter to each sub-stream defined in the previous step as boundaries.

In the following paragraphs, we provide more details about the implementation of these steps.

Appendix A.2: Detecting jumps in the calibration constants

We discuss here how our data analysis code is able to determine the presence of sudden jumps in the calibration constant that are due to some real change in the state of the radiometer. This task is not trivial, since we must look for sudden variations in a stream of numbers (the K values) that is dominated by statistical noise; moreover, the rms of such data changes with time, since it is correlated with the amplitude of the dipole signal ΔT_dip (see Fig. 4). Within a region of high rms, it is therefore possible to mistake a sudden change in the value of K due to statistical fluctuations with an intrinsic change in the radiometer’s calibration. We have therefore developed a figure of merit that allows us to disentangle these two families of jumps.

To implement our figure of merit, we have defined a procedure to quantify the level of statistical fluctuations in the data. Such a procedure is similar to a smoothing filter and we describe it with the aid of Fig. A.1, which shows its application to some real data (the values of K for radiometer LFI-27M, calculated during the first two sky surveys). As already described in Sect. 3, the variation in the statistical noise in K (panel A) is related to the strength of the dipole signal, i.e., the amplitude of ΔT_dip (panel B). Therefore, in order to properly weight the importance of variations in the value of K, we apply a moving average to the stream of K values, where the amplitude of the window (panel C) depends on the value of ΔT_dip. The result of the moving average (K_sm) is shown in panel D of Fig. A.1. For each window, the code computes the rms, σ_K, of the N values. This quantity depends both on the statistical noise and on the presence of real jumps in the value of K. Therefore, we use the expression $\mathrm{\Delta }{\mathit{T}}_{\mathrm{dip}}\mathrm{\times }\frac{{\mathit{\sigma }}_{\mathit{K}}}{\mathit{N}}$(A.1)as a figure of merit for determining the presence of real jumps. This quantity is proportional to the rms of the moving average but is weighted by the amplitude of the calibration signal: the stronger the latter, the more likely that a real jump is present in the N samples.

The threshold used with Eq. (A.1)is defined in terms of percentiles, specifically, we consider all the data that are greater than the nth percentile to mark the presence of a real jump. The value for n depends on the radiometer: for 30 GHz radiometers it is 99; for 44 GHz it is 99.9; and for 70 GHz it is 99.5. Such values have been determined by considering the quality of the null¹⁷ test (see Sect. A.4).

Appendix A.3: The smoothing algorithm

Once the positions of the jumps have been determined, the code applies a smoothing algorithm to each subset of values of the original K (produced by Da Capo) between two consecutive jumps. The algorithm applies a low-pass filter in the Fourier

domain that retains only 5% of the lowest frequencies. After this step, to further reduce the noise, we apply a moving average to the result, where each sample is weighted by the value of ΔT_dip.

Appendix A.4: Validation of the algorithm

Since an incorrect identification of a real jump is likely to produce stripes in maps, we ran a number of null tests in order to optimize the free parameters of the smoothing algorithm (e.g., the threshold used to detect jumps and the widths of the moving average windows). We calibrated the data using a number of combinations of parameters and produced single-survey sky maps (i.e., maps obtained using six months of data). We then differenced them, under the hypothesis that a perfect calibration would produce a map where the value of each pixel is consistent with zero. Figure A.2 shows an example of this analysis and shows that the jump found by the code near day 257 is likely to be a real jump, because not considering it as such leads to a stronger stripe in survey-difference maps.

All Tables

All Figures