© ESO, 2014

1. Introduction

This paper, one of a set associated with the 2013 release of data from the Planck¹ mission (Planck Collaboration I 2014), describes the impact of the optical system, detector response, analogue and digital filtering and the scan strategy on the determination of the High Frequency Instrument (HFI) angular power spectra. An accurate understanding of this response, and the corresponding errors, is necessary in order to extract astrophysical and cosmological information from cosmic microwave background (CMB) data (Hill et al. 2009; Nolta et al. 2009; Huffenberger et al. 2010).

Bolometers, such as those used in the HFI on Planck, are phonon-mediated thermal detectors with finite response time to changes in the absorbed optical power. The observations are affected by attenuation and the phase shift of the signal resulting from a detector’s thermal response, as well as the analogue and digital filtering in the associated electronics.

In the small signal regime (appropriate for CMB fluctuations and Galactic emission), the receiver response can be well approximated by a complex Fourier domain transfer function, termed the time response. The time-ordered data (also referred to as time-ordered information, or TOI) are approximately deconvolved by the time response function prior to calibration and mapmaking (for recent examples of CMB observations with similar semiconductor bolometers see Tristram et al. 2005; Masi et al. 2006).

Ideally, the deconvolved TOIs represent the true sky signal convolved with the optical response of the telescope (or physical beam) and filtered by the TOI processing. The combination of time domain processing and physical beam convolution is, in practice, degenerate, due to the nature of the scan strategy; the Planck spacecraft scans at 1 rpm, with variations less than 0.1% (Planck Collaboration I 2011). The beams reconstructed from the deconvolved planetary observations are referred to as the scanning beam.

These deconvolved data are then projected into a pixelized map, as discussed in Sect. 6. of Planck Collaboration VI (2014). To a good approximation, the effect of the mapmaking algorithm is to average the beam over the observed locations in a given pixel. This average is referred to as the effective beam, which will vary from pixel to pixel across the sky. The mapmaking procedure implicitly ignores any smearing of the input TOI; no attempt is made to deconvolve the optical beam and any remaining electronic time response. Thus, any further use of the resulting maps must take the effective beam into account.

To obtain an unbiased estimate of the angular power spectrum of the CMB, one must determine the impact of this effective beam pattern on a measurement of an isotropic Gaussian random signal in multipole (ℓ) space. The filtering effect is well approximated by a multiplicative effective beam window function, derived by coupling the scan history with the scanning beam profile, which is used to relate the angular power spectrum of the map to that of the underlying sky (Hivon et al. 2002).

The solar and orbital motion of Planck with respect to the surface of last scattering provides a 60-s periodic signal in the time ordered data that is used as a primary calibrator (Planck Collaboration V 2014; Planck Collaboration VIII 2014). This normalizes the window function at a multipole ℓ = 1; the effective beam window function is required to transfer this calibration to smaller angular scales.

These successive products (the scanning beam, the effective beam and the effective beam window function) must be accompanied by an account of their errors, which are characterized through ensembles of dedicated simulations of the planetary observations. The error on the effective beam window function is found to be sub-dominant to other errors in the cosmological parameter analysis (Planck Collaboration XVI 2014).

The scanning beam is thus measured with on-orbit planetary data, coupling the response of the optical system to the deconvolved time response function and additional filters in the TOI processing. As shown in the following, the main effect of residual deconvolution errors is a bias in the beam window function of order 10^-4 at ℓ > 100, due to a residual tail in the beam along the scan direction.

The scanning beam can be further separated into the following components;

• 1.

  The main beam, which is defined as extending to30′ from the beam centroid.

• 2.

  The near sidelobes, which extend between 30′ and 5°. These are typically features below −30 dB, and include the optical effects of phase errors, consisting of both random and periodic surface errors and residuals due to the imperfect deconvolution of the time response.

• 3.

  The far sidelobes, which extend beyond 5°. These features are dominated by spillover: power coupling from the sky to the feed antennas directly, or via a single reflection around the mirrors and baffles. The minimum in the beam response between the main beam and direct spillover of the feed over the top of the secondary mirror falls at roughly 5°, making such a division natural (Tauber et al. 2010).

This paper describes the main beam and the near sidelobes, the resulting effective beam patterns on the sky and the effective beam window function used for the measurement of angular power spectra, along with their errors. The effects of the far sidelobes are mainly discussed in Planck Collaboration VI (2014) and Planck Collaboration XIV (2014), and their effects on calibration described in Planck Collaboration X (2014). A companion paper (Planck Collaboration IV 2014) computes the effective beams and window functions for Planck’s Low Frequency Instrument (LFI). Despite using very different methods which depend more strongly on optical modelling, the two instruments produce compatible power spectra, providing a cross-check of both approaches (Planck Collaboration XXXI 2014).

The signal-to-noise ratio of Jupiter observations with HFI allows the measurement of near sidelobes to a noise floor of − 45(−55) dB relative to the forward gain at 100 (857) GHz. The HFI analysis does not rely on a physical optics model to constrain the behaviour of the beam in this regime.

On-orbit measurements of the Planck HFI scanning beam and temporal transfer function have been previously reported in Planck HFI Core Team (2011a) for the Early Release Compact Source Catalog (ERCSC; Planck Collaboration VII 2011) and early science from Planck. Section 2 presents an improved model for the time response and explains how its parameters are derived using on-orbit data, and how it is deconvolved from the data. Figure 1 provides a flowchart of the determination of the time response parameters. Figure 2 is the flowchart of the steps that lead from planet measurements to effective beams, effective beam window functions and assessment of uncertainties. Section 3 describes how the scanning beams are reconstructed from planet observations. Section 3.4 specifically describes the effects of long time scale residuals in the data due to imperfect deconvolution of the time response. Section 4 describes the propagation of the scanning beams to effective beams and effective beam window functions using the scanning history of Planck. Section 5 describes the techniques used to propagate statistical and systematic errors and to check the consistency of the beams and window functions. Section 6 describes the final error budget and the eigenmode decomposition of the errors in the effective beam window function.

2. Time response

In Planck’s early data release, a 10-parameter model TF10 describes the time response of the bolometer/electronics system (Planck HFI Core Team 2011a). This section describes an improved model of the time response based on the HFI readout electronics schematics which more accurately reproduces phase shifts in the system close to the Nyquist frequency. The improved model also provides more degrees of freedom for the bolometer’s thermal response in order to describe more accurately the low frequency response of the bolometer.

2.1. Model

The new model is named LFER4 (Low Frequency Excess Response with four time constants) and consists of an analytic model of the HFI readout electronics (Lamarre et al. 2010) and four thermal time constants and associated amplitudes for the bolometer: $\mathrm{TF}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{=}\mathit{F}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}{\mathit{H}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\omega }\mathrm{;}{\mathit{S}}_{\mathrm{phase}}\mathit{,}{\mathit{\tau }}_{\mathrm{stray}}\mathrm{\right)}\mathit{,}$(1)where TF(ω) represents the full time response as a function of angular frequency ω. The time response of the bolometer alone is modelled by $\mathit{F}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}\mathit{,}\mathrm{4}}\frac{{\mathit{a}}_{\mathit{i}}}{\mathrm{1}\mathrm{+}\mathrm{i}\mathit{\omega }{\mathit{\tau }}_{\mathit{i}}}\mathit{,}$(2)and H′(ω;S_phase,τ_stray) is the analytic model of the electronics transfer function (whose detailed equations and parameters are given in Appendix A) with two parameters.

The overall normalization of the transfer function is forced to be 1.0 at the signal frequency of the dipole, leaving a total of 9 free parameters for each bolometer.

The sum of single-pole low-pass filters represents a lumped-element thermal model with four elements. The thermal model underlying the temporal transfer function is described elsewhere (Spencer et al., in prep.); this work adopts an empirical approach to correcting the data.

The two parameters of H′ mainly affect the high frequency portion of the time response. S_phase represents the phase difference between the bias and the lock-in summation, and is fixed in the model as a readout electronics setting. The second parameter τ_stray is the time constant of the bolometer resistance and the parasitic capacitance of the wiring and is measured independently during the checkout and performance verification (CPV) phase of the mission prior to the sky survey. All resistance and capacitance values of the readout electronics chain are fixed at values from the circuit diagram.

The in-flight data are used to determine the remaining seven free parameters. Low frequency parameters are constrained by minimizing the difference between the first and second survey maps (Sect. 2.2), while planetary observations are used to constrain those parameters governing the high frequency portion of the time response (Sect. 2.3).

The fastest thermal time constants in the LFER4 model roughly correspond to the time constants measured during pre-launch tests of the bolometers (Holmes et al. 2008; Pajot et al. 2010). The slower time constants contribute lower frequency response at the several percent level. The time constants in the LFER4 model are not exactly identical to those measured on glitches (Planck Collaboration X 2014) due to additional filtering applied by the deglitching. A future publication (Spencer et al., in prep.) will relate the time constants and amplitudes to the thermal properties of the bolometer and module.

2.2. Fitting slow time constants with galaxy residuals

As Planck scans across the Galactic plane, the low frequency time response of HFI spreads Galactic signal away from the plane. Surveys 1 and 2 consist of roughly six months of data each and cover nearly the same sky, scanned at almost opposite angles. The difference between maps made with data from the two individual surveys highlights the effect since the Galactic power is spread in different directions in the two surveys, creating symmetric positive and negative residuals in the difference map. The LFER4 parameters are varied to minimize the difference between these surveys. For most bolometers, the fit is limited to the slowest time constant and its associated amplitude, and in others the fit is extended to the two slowest time constants and associated amplitudes.

Other systematic effects can confuse the measurement of LFER4 parameters by creating a similar positive/negative residual pattern in the survey difference. The philosophy employed here is to minimize the survey residuals fitting only for LFER4 parameters, but to use simulations and data selections to test the dependency of the results with other systematics.

2.2.1. Survey difference method 1

Two techniques are used to perform the fit. The first method is based on map re-sampling in the time domain, using the pointing to generate synthetic TOI. The synthetic TOI of each survey is compared with the synthetic TOI of both surveys combined. Before the production of these synthetic TOI, the maps are smoothed to 30′. Given the fact that in consecutive surveys the scan direction is nearly opposite, the survey 2 TOI is very similar to the time-reversal of the survey 1 TOI. This symmetry is assumed to be exact, ensuring that time reversal is equivalent to taking the complex conjugate in the frequency domain. Since the sky signal has been convolved by the true transfer function, TF_true(ω), and deconvolved by the estimated – not fully correct – transfer function TF₀(ω), the single survey re-sampled time-stream is ${\mathit{d}}_{\mathrm{1}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathrm{ℱ}}^{-1}\mathrm{\left(}{\mathrm{TF}}_{\mathrm{true}}\mathrm{\right(}\mathit{\omega }\mathrm{)}\mathit{/}{\mathrm{TF}}_{\mathrm{0}}\mathrm{(}\mathit{\omega }\mathrm{\left)}\mathrm{ℱ}\mathrm{\right(}\mathit{s}\mathrm{\left(}{\mathit{t}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{)}\mathit{,}$(3)where ℱ is the Fourier transform operator, s is the true sky map observed at time t′, and d₁ is the synthetic TOI. This can be written as, ${\mathit{d}}_{\mathrm{1}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathrm{ℱ}}^{-1}\mathrm{\left(}\mathit{\delta }\mathrm{TF}\mathrm{\right(}\mathit{\omega }\mathrm{\left)}\mathrm{ℱ}\mathrm{\right(}\mathit{s}\mathrm{\left(}{\mathit{t}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{)}\mathit{,}$(4)having defined δTF(ω) = TF_true(ω)/TF₀(ω) as the corrective factor which should be applied to the estimated transfer function. The synthetic TOI of survey 2, d₂, is similarly defined. For the TOI of both surveys combined, the average of the two maps is used.

Using the time reversal property, the synthetic TOI d(t) of the full sky map (survey 1 and survey 2 combined) may be written as $\begin{array}{ccc}\mathit{d}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathrm{ℱ}}^{-1}\left(\frac{\mathit{\delta }\mathrm{TF}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{+}\mathit{\delta }{\mathrm{TF}}^{\mathrm{\ast }}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}}{\mathrm{2}}\mathrm{ℱ}\mathrm{\left(}\mathit{s}\mathrm{\right(}{\mathit{t}}^{\mathrm{\prime }}\mathrm{\left)}\mathrm{\right)}\right)\mathrm{·}& & \end{array}$In the Fourier domain, the ratio of the TOIs is $\begin{array}{ccc}\mathit{R}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{=}\frac{\mathrm{ℱ}\mathrm{\left(}{\mathit{d}}_{\mathrm{1}}\mathrm{\right(}\mathit{t}\mathrm{\left)}\mathrm{\right)}}{\mathrm{ℱ}\mathrm{\left(}\mathit{d}\mathrm{\right(}\mathit{t}\mathrm{\left)}\mathrm{\right)}}\mathrm{=}\frac{\mathit{\delta }\mathrm{TF}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}}{\mathrm{\Re }\mathrm{\left(}\mathit{\delta }\mathrm{TF}\mathrm{\right(}\mathit{\omega }\mathrm{\left)}\mathrm{\right)}}\mathrm{=}\mathrm{1}\mathrm{+}\mathrm{i}\frac{\mathrm{\Im }\mathrm{\left(}\mathit{\delta }\mathrm{TF}\mathrm{\right(}\mathit{\omega }\mathrm{\left)}\mathrm{\right)}}{\mathrm{\Re }\mathrm{\left(}\mathit{\delta }\mathrm{TF}\mathrm{\right(}\mathit{\omega }\mathrm{\left)}\mathrm{\right)}}\mathrm{·}& & \end{array}$Since the real part of the transfer function at low frequency is close to 1, the imaginary part of the ratio of the synthetic beams is equal to the imaginary part of the corrective factor, δTF(ω): $\begin{array}{ccc}\mathrm{\Im }\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\mathrm{\Im }\mathrm{\left(}\mathit{\delta }\mathrm{TF}\mathrm{\right(}\mathit{\omega }\mathrm{\left)}\mathrm{\right)}\mathit{.}& & \end{array}$The LFER4 model of the transfer function is fitted to this measure of the imaginary part of the correction to infer the parameters of the true transfer function, in the low frequency regime, typically below a few Hz.

2.2.2. Survey difference method 2

The second method looks at one-dimensional slices through the Galactic plane for each survey independently. A sky signal slice is obtained by resampling a sky map for a single survey and a single bolometer. The slices are taken along the scan direction and the sky signal is averaged over 5° in longitude. Only the sky region close to the Galactic plane is considered (10° above and below the Galactic plane and longitudes between −40° and 60°). The slice is convolved with the ratio of the transfer function used to create the map and a new LFER4 transfer function with trial parameter values. New sky survey maps are obtained by re-projecting the slices into pixelized maps. As for the previous method, parameters of the low-frequency part of the LFER4 transfer functions are chosen so that they minimize the residuals in the difference between surveys 1 and 2.

In practice, the first method is used for the 100–353 GHz bolometers, and at 545 GHz and 857 GHz the second method is used, being better adapted to the maps having the most structure in the Galactic difference residuals. Figure 3 shows an example of the residual remaining in a HEALPix² map (Górski et al. 2005) of the survey difference, after fitting the long time response, showing reduced asymmetry in the Galactic plane.

2.2.3. Survey difference systematics

In the survey difference solution for the time response, any systematic effect that creates a residual in the survey difference can be confused with a time response effect, in particular affecting the low frequency time response. This section identifies a number of residual-producing systematics and quantifies the resulting bias in the transfer function. These residual-producing systematics include far sidelobes, zodiacal light, pointing offset uncertainty, gain drifts, main beam asymmetry and polarized sky signal.

Due to the very high signal-to-noise ratio of Galactic signal at sub-millimetre wavelengths, far sidelobe pickup of the Galactic plane is detected in the 545 GHz and 857 GHz channels. A physical optics model of the far sidelobe pickup is used to estimate the signal from the Galactic centre. The optical depth of the zodiacal dust cloud along the line of sight varies as Planck orbits the Sun, leaving a residual in the survey difference (Planck Collaboration XIV 2014). A model of the zodiacal dust emission is subtracted in the reconstruction of the time response. The reconstruction of the time response is then repeated without subtracting the models; these do not significantly affect the result (Fig. 4).

As a probe of the effect of far sidelobes on the time constant determination, the pipeline is run on a survey difference map obtained from the sidelobe model alone for each of the 100, 143, 217, and 353 GHz channels. The method did not find a long time constant, as the sidelobe effect on the survey map is very different from the time constant effect.

A systematic offset in the pixel pointing creates a residual in the survey difference. The pointing solution reduces the pointing error to a few arcseconds RMS in both the co-scan and cross-scan direction. With the 6° s^-1 scanning speed, this error corresponds to frequencies greater than 1 kHz, far from the range affected by the long time constant.

Gain variability can also bias the estimation; due to non-linearity in the analogue to digital converters (ADC), the HFI responsivity to the sky signal varies at a few tenths of a percent throughout the mission. As a probe of this effect, gain-corrected data for the 100, 143 and 217 GHz bolometers are used to reconstruct the long time response. This has a negligible impact on the fitted parameters. Residual gain errors tend to leave monopolar residuals that are not coupled to the long time constants in the fitting procedure.

Some residual is expected in the survey difference maps because the asymmetric beam scans the sky at different angles in the two surveys. This is especially an issue at 545 GHz where the beam is substantially asymmetric. To quantify the amplitude, simulated survey difference maps are generated using a realistic asymmetric beam model convolved with the Planck Sky Model (Delabrouille et al. 2013); the residuals in the Galactic plane at 545 GHz are dominated by the main beam asymmetry (see Fig. 5).

Polarization sensitive bolometers (PSBs) show an additional residual in the survey difference maps because of polarized sky signal observed at slightly different crossing angles in the two surveys. For the PSBs, the low frequency time response is therefore determined using different levels of polarization masking. These studies do not suggest the presence of any significant level of bias from polarization, but only higher noise with wider masking. As an additional check on the contribution of residual polarization to the low frequency response, the survey difference of the sum of the two arms of each PSB pair is input to the pipeline. This sum is not sensitive to polarization, and no bias is found in the determination of the time response.

2.3. Fitting fast time constants with planet crossings

Filtering of the TOIs and errors in the deconvolution kernel results in ringing along the scan direction. The planets Mars, Jupiter and Saturn are high signal-to-noise sources that can be used to minimize this smearing by adjusting the parameters of the LFER4 model. See Sect. 3.1 for a description of the planet data.

In solving for the high frequency portion of the time response, the beam profile is forced to be compact. The optical beam is modelled as a spline function on a two-dimensional grid (see Appendix C for details) and the LFER4 parameters are fit by forcing conditions on the resulting beam shape (Fig. 6).

The planet data are first deconvolved with a time response model derived from pre-launch data, to recover an initial estimate of the beam profile. Jupiter is used for the 100, 143 and 217 GHz channels, while Saturn is used for higher frequencies (see Sect. 3.4 for a discussion of the non-linearity of Jupiter at higher frequencies).

Rather than deconvolve the planet data, the model parameters are determined in the forward sense. Since the beam is decomposed into B-spline functions, this basis is convolved with the temporal transfer function to retrieve the coefficients for each basis function using the planet data. These coefficients are applied to the original deconvolved B-spline functions to recover an estimate of the optical beam. The convolution is made in the Fourier domain by re-sampling the B-spline function onto a timeline with a sample separation corresponding to 4.̋5. The typical knot separation length of the basis function is between 1′ and 2′.

In the Fourier domain, the convolution of the temporal transfer function with the planet signal is $\mathrm{𝒫𝒞}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{=}{\mathrm{ℬ}}^{\mathrm{0}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}{\mathrm{TF}}^{\mathrm{0}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathit{,}$(5)where ℬ⁰(ω) is the Fourier transform of the slice through the peak planet signal in the scan direction b⁰(t), which is obtained by de-convolving planet data using a fiducial estimate, TF⁰(ω) of the transfer function. The slice b⁰(t) is then symmetrized about the origin (defined by the location of the maximum of the B-spline representation)to obtain b^sym(t), and its Fourier transform ℬ^sym(ω). This operation aims to recover a beam that, by construction, is symmetric, within the limits allowed by the model of the temporal transfer function, $\mathrm{𝒫𝒞}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{=}{\mathrm{ℬ}}^{\mathrm{sym}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}{\mathrm{TF}}^{\mathrm{\ast }}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathit{.}$(6)Here ℬ^sym(ω) is entirely defined by and TF⁰(ω) and the new estimate of the time response TF^∗(ω) is derived from Eq. (6). This function is parameterized in terms of the three shortest time constants (τ₁, τ₂, τ₃) and their associated amplitudes (a₁,a₂,a₃).

2.4. Stationarity of the time response

The time response of each HFI detector/readout channel is a function of the cryogenic temperature of the bolometers and the ambient-temperature components of the readout electronics. Both cryogenic and ambient temperatures change throughout the mission as the Galactic particle flux and the Planck spacecraft solar distance are modulated. The seasonal consistency of the scanning beam sets an upper limit on changes in the time response through the mission, shown below in Sect. 5.2.1.

2.5. Deconvolution of the data

The time response transfer function is deconvolved from the data and not included as part of the scanning beam, because the low frequency response of the bolometer would give an extended scanning beam, stretching many degrees from the main lobe along the scan direction.

Since the time response amplitude decreases as a function of frequency, the deconvolution operation increases the noise at high frequency to an unacceptably high fraction of the RMS. During the deconvolution stage of the TOI processing, a phaseless low-pass filter is applied in order to suppress the high frequency noise and keep pixel aliasing at a manageable level.

In the early data release, a finite impulse response low-pass filter was used for this purpose (Planck HFI Core Team 2011b). In the 2013 cosmology data release, the low-pass filter is further tuned for the 100, 143, 217 and 353 GHz channels to reduce filter ripple produced by the lowpass filters used in the early-release data. The filter is now implemented in the Fourier domain, with a kernel consisting of the product of a Gaussian and a squared cosine, $\mathit{K}\mathrm{\left(}\mathit{f}\mathrm{\right)}\mathrm{=}{\mathit{K}}_{\mathrm{1}}\mathrm{\left(}\mathit{f}\mathrm{\right)}{\mathit{K}}_{\mathrm{2}}\mathrm{\left(}\mathit{f}\mathrm{\right)}\mathit{,}$(7)where ${\mathit{K}}_{\mathrm{1}}\mathrm{\left(}\mathit{f}\mathrm{\right)}\mathrm{=}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}\mathit{f}\mathit{/}{\mathit{f}}_{\mathrm{Gauss}}{\mathrm{)}}^{\mathrm{2}}}$(8)and ${\mathit{K}}_{\mathrm{2}}\mathrm{\left(}\mathit{f}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{if}\mathit{f}\mathit{<}{\mathit{f}}_{\mathrm{c}}\mathit{,}\\ & \mathrm{if}{\mathit{f}}_{\mathrm{c}}\mathit{<}\mathit{f}\mathit{<}{\mathit{f}}_{\max }\mathit{,}\\ & \mathrm{if}\mathit{f}\mathit{>}{\mathit{f}}_{\max }\mathit{.}\end{array}$(9)Here f_max = f_c + k(f_samp/2 − f_c) and f_samp is the sampling frequency of the data. The parameters of the filter are the same for all bolometers in the bands 100–353 GHz: f_Gauss = 65 Hz; f_c = 80 Hz; and k = 0.9. To first order, this filter widens the scanning beam along the scan in an equivalent way to convolving the optical beam with a Gaussian with full-width at half-maximum (FWHM) of 2.́07. The filter introduces some rippling along the scan direction at the −40 dB level at 217 and 353 GHz, where the beams allow harmonic signal content close to the filter edge. The rippling is captured by the B-spline beam representation (see Fig. 11).

The 3-point finite impulse response filter is still used for the 545 and 857 GHz channels³. This extends the scanning beams along the scan direction more than the Gaussian-cosine Fourier filter (the 545/857 GHz filter time scale corresponds to 3.́0 on the sky), but has the advantage of reducing rippling and signal aliased from above the Nyquist frequency.

The deconvolution kernel multiplied by the data in the Fourier domain is the product of the lowpass filter with the inverse of the bolometer/electronics time response, $\mathit{D}\mathrm{\left(}\mathit{f}\mathrm{\right)}\mathrm{=}\mathit{K}\mathrm{\left(}\mathit{f}\mathrm{\right)}\mathit{/}\mathrm{TF}\mathrm{\left(}\mathit{f}\mathrm{\right)}\mathit{.}$(10)Figure 7 shows a comparison of the deconvolution functions resulting from the LFER4 model and from the TF10 model used in the ERCSC data. The differences between the two models appear mainly in the phase at high frequency, mostly above the signal frequency corresponding to the beam size. Although the phase of LFER4 is a more accurate description of the system, in practice replacing LFER4 with TF10 had little effect on the data, because of the lowpass filter applied at the time of deconvolution and the empirical determination of an overall sample offset in the pointing reconstruction.

In the HFI data processing (Planck Collaboration VI 2014), data chunks of length 2¹⁹(≈5 × 10⁵) samples are Fourier transformed at a time, overlapping by half with the subsequent chunk to avoid edge effects.

3. Scanning beams

The filtering of the TOI and the accuracy of the deconvolution kernel affect the angular response of the HFI detectors. An accurate estimate of the scanning beam, resulting from the filtering of the physical (optical) beam by these time domain convolutions, is required to relate the angular power spectra of the maps to that of the underlying sky. As described in Planck HFI Core Team (2011a) and Planck HFI Core Team (2011b), the HFI scanning beam profiles are measured using the planetary observations. HFI uses two flat sky representations of the two-dimensional scanning beams, one using Gauss-Hermite polynomials, and another using B-spline functions.

Three selections of planetary data are used to derive estimates of different aspects of the beam:

• 1.

  the first two observing seasons of Mars (main beam, andwindow functions);

• 2.

  all available seasons of Mars [3], Jupiter [5] and Saturn [4] (near sidelobe); and

• 3.

  all five Jupiter observations (residual time response).

The effective beam window functions used in the CMB analysis are ultimately derived from the first of these. In each case, the statistical properties of the beam representations and the choice of planetary data are studied using ensembles of simulated planet observations (Sect. 5).

While the signal-to-noise ratio of the Jupiter and Saturn data is significantly greater than the Mars data, at this stage of the analysis a reconstruction bias results in the main beams recovered from simulated Jupiter and Saturn observations that is not present in the simulations of Mars. Additionally, the non-linear response of some HFI detectors to the Jupiter signal (Planck HFI Core Team 2011a) makes the normalization of the planet peak response uncertain at the few percent level (see Sect. 3.4). Therefore the main beam model is established using Mars data, while observations of Jupiter and Saturn are used only to estimate the near sidelobes and residuals in the deconvolution of the transfer function.

The B-spline representation of the joint Mars observations is used as input for the calculation of the effective beam and the effective beam window function. Simulations have shown the B-spline representation to better capture the features outside of the main lobe, due to the necessarily finite order of the Gauss-Hermite decomposition. However, the Gauss-Hermite model is used for other systematics checks, including the consistency of the planets and observing seasons.

Because the Jupiter and Saturn data allow measurement of the beam response below −45 dB from the peak, there is no need to rely on a model of the telescope to determine the main beam or near sidelobe structure.

3.1. Planetary data handling

The JPL Horizons package⁴ (Giorgini et al. 1996) is programmed with Planck’s orbit to calculate the positions of the planets. Table 1 shows the dates when the planets were within 2° of the centre of the focal plane. By the end of HFI operations Mars was observed three times, Saturn four times, and Jupiter five.

The planets Jupiter, Saturn and Mars are among the brightest compact objects seen by Planck HFI; the signal amplitude affects the data handling in a number of ways. Moving solar system objects are flagged and removed from the TOI in the standard processing pipeline. Specialized processing for the planet data is required, with two major differences from the nominal processing (see Planck Collaboration VI 2014 for details).

While the pickup from the ⁴He-JT cooler is removed from the data as usual, pointing periods containing very bright sources such as the planets cannot be used to reliably estimate the line amplitude. Instead, during the planet observations, these are extrapolated from neighbouring pointing periods.

To better detect glitches near the extremely bright planet crossings, an estimate of the planet signal is subtracted from the data prior to glitch detection. Glitch template subtraction is performed on the signal-subtracted timeline in the same way as during nominal observations.

In order to remove the (quasi-stationary) astrophysical background from the planetary data, a bilinear interpolation of the frequency averaged map is subtracted. The full mission map is used for the five-season Jupiter data, while the nominal survey sky maps are used in the processing of the other planetary data.

The planets are oblate spheroids, and appear as ellipses slightly extended along the direction of the ecliptic. The Planck planet range and Planck sub-latitude calculated by JPL Horizons are used in combination with the polar and equatorial radii of the planets reported by Horizons to compute the angular size and ellipticity of each planet. During HFI observations, the mean angular radii of Jupiter, Saturn and Mars are 20.̋44, 8.̋542, and 4.̋111, respectively. The ratio of the equatorial to polar radii are 1.069, 1.109 and 1.006, respectively.

The finite planetary disc size increases the apparent size of the scanning beam and biases the inferred effective beam window function. The filtering in multipole space of a circular disc of angular radius R can be written as B_disc(ℓ) = 2J₁(ℓR) /(ℓR), where J₁(x) is the Bessel function of order 1. Figure 8 shows the ${\mathit{B}}_{\mathrm{disc}}^{\mathrm{2}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}$ for the three planetary discs. In practice, the effects of the disc size are mitigated by merging observations from the three brightest planets. The effects of the large Jupiter disc size are greatest where the spatial gradient of the beam is greatest, between the –3 to −10 dB contours of the beam. By excluding the Jupiter observations above − 10 dB, the disc size smearing is reduced, and setting a −20 dB threshold results in a bias in the window function below 10^-3 at all multipoles.

Planck observes Saturn at a range of ring inclination angles: 6.̊03, 2.̊45, 12.̊6 and 9.̊4 for seasons 1, 2, 3 and 4, respectively. While emission from the solid angle of the ring increases the effective planetary disc area, the brightness temperature of the ring tends to be much less than the planetary disc temperature for the Planck bands. For example Weiland et al. (2011) fit a single 90 GHz ring brightness temperature 14% that of the Saturn disc. In our beam reconstruction the average of the first three Planck observations of Saturn is used, and even in the extreme limit where the ring brightness temperature is the same and only Saturn data are used, the multipole space correction is 2 × 10^-3 at ℓ = 3000; this correction is ignored.

The elliptical shape of the planetary disc gives a further bias in the inferred window function, depending on whether the long axis of the beam is aligned with or perpendicular to the long axis of the planet. In this case the planetary disc is approximated as an elliptical Gaussian with σ = 0.5R. The worst case is a 10^-3 effect in the case of 100 GHz beams measured on Jupiter at ℓ = 3000 (where the 100 GHz window function is vanishingly small). At 143 and 217 GHz, the 10^-3 level is reached only at ℓ = 4000. The effect is negligible in the range of Planck’s sensitivity. With Saturn and Mars, the ellipticity effects are <10^-4 and <10^-6, respectively, at all multipoles; the ellipticity of the planetary discs has a negligible effect on the estimate of the scanning beam.

3.2. Main beam model

Two representations are used to describe the (two dimensional) main beam; a Gauss-Hermite basis (described in Appendix B following Huffenberger et al. 2010 and Planck HFI Core Team 2011a) and a B-spline basis (Appendix C).

The Gauss-Hermite (GH) and B-spline bases have very different characteristics. The GH representation uses a relatively small number of parameters and, in practice, amounts to a perturbative expansion about a Gaussian shape. The B-spline is quite general, using many more degrees of freedom to fit the data on a defined grid, with the spline controlling the behaviour in between. The bias and correlation structure of the noise of these two representations are characterized using Monte Carlo simulations of the planetary data which include all the details of the beam-processing pipeline used on the data (these simulations are described more completely in Sect. 5). In each case, the parameters of the representation are derived directly from the time-ordered data, without recourse to a pixelized map.

Figure 9 shows the B-spline scanning beams for the entire HFI focal plane, as reconstructed from the Jupiter and Saturn data. Figure 10 shows the radially binned, frequency averaged beam profile for the HFI channels, comparing the B-spline representation of the Mars data (solid black lines) with the combined Mars, Jupiter and Saturn data (filled and open points, and the blue dashed line). The B-spline maps are apodized with a Gaussian at a radius beyond the signal-to-noise floor of the Mars data: 13.́4, 13.́0, 11.́4, 12.́9, 17.́8, 17.́8 on average at 100, 143, 217, 353, 545,and 857 GHz respectively.

The azimuthally averaged beam window function, ${\mathit{B}}_{\mathit{\ell }}^{\mathrm{2}}$, from each of these models is compared to the known input beam model. At 100, 143 and 217 GHz the two methods perform comparably, with the B-spline having slightly smaller bias and variance; at the higher frequencies the B-spline performs demonstrably better, especially at 545 and 857 GHz, as expected due to the highly non-Gaussian shape of these multi-mode detectors.

3.3. Near sidelobes

While the HFI beams are Gaussian at the −25 dB level, non-trivial structure in the beam is captured in the data at lower power. There are two distinct components to the near lobe response: a discrete pattern of secondary lobes evident at frequencies of 217 GHz and above; and a diffuse shoulder consistent with phase errors in the aperture plane.

The Planck reflectors suffer from print-through of the honeycomb structure that supports the carbon-fiber face sheets. The size of the deformation has been measured during thermal testing to be less than 20 μm (Tauber et al. 2010). While small in amplitude, the strict periodicity of this pattern results in a correspondingly periodic structure in the near lobes, seen clearly in Fig. 12, which is slightly larger than predicted based on the pre-launch measurements. A simple grating equation describes the angular positions of the resulting contributions to the near sidelobes, $\sin {\mathit{\theta }}_{\mathit{n}}\mathrm{=}\frac{\mathit{n\lambda }}{\mathit{Yd}}\mathrm{·}$(11)where θ_n is the angular position of the order n lobe from the central beam peak, λ is the wavelength of the radiation, d is the grating periodicity and Y is a factor that describes the position of each reflector along the optical path, with Y = 1.00 for the primary reflector and Y = 1.80 for the secondary reflector. Three possible periodicities (19.6 mm, 30 mm, 52 mm) in the honeycomb array dominate the Planck dimpling pattern for the 857 GHz detectors, though only those for the 52 mm periodicity can be seen for the 545 GHz and 353 GHz detectors. For the highest frequency detectors, only the weaker lobes due to the 19.6 mm and 30 mm periodicities are seen outside the 40′ main beam model, but they contribute at most (0.050 ± 0.008)% to the integrated beam solid angle. A forthcoming publication (Oxborrow et al., in prep.) will present an in-depth study of the mirror surface.

The second component is a beam shoulder becoming significant near the −30 dB contour, and extending to 2–4 FWHM from the beam centre. This shoulder is consistent with scattering from random surface errors on the primary and secondary reflectors, and is reasonably well-described by a spectrum of surface errors with correlation lengths ranging from 2 to 12 cm, with an RMS of order 10 μm (Ruze 1966). The contribution of each of these components is included in the radially binned profiles shown in Fig. 10.

While the B-spline parameterization captures both the main beam and near sidelobe structure, the extended features must be separately included in the Gauss-Hermite beam representation. Figure 11 shows contour plots of a B-spline beam using Mars, Jupiter and Saturn data at each frequency.

The B-spline representation of the scanning beam used to compute the window function includes only that portion of the near sidelobe structure that falls within the signal-to-noise of the Mars data; the Jupiter and Saturn data provide an estimate of the beam solid angle that is neglected in the scanning beam product. The near sidelobe solid angle, and the resulting window function error, are sensitive to the details of the analysis, including the sky subtraction, offset removal and masking of in-scan ringing at the part-per-million level. Although a comprehensive study of these effects in the Jupiter and Saturn data are underway, a conservative, and model independent upper limit is obtained by taking the envelope of the noise floor to define the maximum solid angle allowed by the data (the red dashed line in Fig. 10). A reasonable estimate of the true solid angle in the near lobes can be obtained by extrapolating the data below the noise floor (the blue dashed line in Fig. 10). By either measure, the grating lobes and diffuse shoulder account for a small fraction of the total beam solid angle; for the 100, 143 and 217 GHz channels this contribution represents less than 0.15% of the total solid angle (see Table 2). The amplitude of the impact on the window function is estimated by comparing the Legendre transform of the maximal envelope of the Jupiter and Saturn data with that of the nominal Mars derived scanning beam. The result is shown as the family of green curves in Fig. 19. Because the Monte Carlo ensembles that are used to derive the error envelope neglect this near sidelobe structure in the beam that is input to the simulations, the window function error amplitudes have been scaled to accomodate the upper limit defined by the noise floor of the Jupiter and Saturn data.

3.4. Residual time response

The Planck spacecraft spin rate is constant to 0.1%, making the time response of the electronics and the bolometer degenerate with the angular response of the optical system.

The B-spline beam model extends ±20′ from the centre of the beam. An error in the time response on fast timescales will thus be exactly compensated by the scanning beam. However, errors in the time response beyond the limit of the scanning beam reconstruction will not be accounted for, and will bias the recovered beam window function.

To look for residual long tails due to incomplete deconvolution of the time response, all five seasons of Jupiter observations are binned into a 2D grid of 2′ pixel size extending 6° from the planet (Fig. 12). These are background-subtracted using the Planck maps and stacked by fitting a Gaussian to estimate the peak amplitude and centroid of each observation. The data are binned as a function of pointing relative to the planet centre.

While a static non-linearity correction is included in the TOI processing, partial ADC saturation and dynamical non-linearity bias the normalization of these tails by underestimating the Jupiter peak signal. An estimate of the non-linear correction is derived by fitting a gridded map of all three Mars observations to a map of all five Jupiter observations for each detector. A signal reduction at the peak of Jupiter is ruled out at the 1% level at 100 and 143 GHz, but detected at higher frequencies. Relative to Mars, the peak Jupiter signal is reduced by non-linearity on average by 3 ± 3%, 12 ± 3%, 12 ± 4% and 65 ± 20% at 217, 353, 545 and 857 GHz, respectively. The tail normalizations are scaled by these factors.

An example of the residual tail is shown in Fig. 13. As well as a tail following the planet due to imperfect deconvolution of the time response, there is also a tail preceding the planet crossing; this is due to the lowpass filter applied in the Fourier domain. The residual beam tails have amplitudes typically at the level of 10^-4 of the peak but extend several degrees from the centre of the beam. The response beyond 6° on the sky is consistent with noise for all detectors.

The signal to noise of the tail measurement greater than 6° from Jupiter sets a model-independent limit on the knowledge of the time response at signal frequencies from 0.016 Hz to 0.1 Hz at <10^-4.

These stacked Jupiter data are integrated to determine the expected bias in total beam solid angle from this remaining uncorrected time response (see Table 2). The mean integral values are typically a few times 10^-4 of the main beam solid angle, an order of magnitude lower than the error in the beam solid angle due to noise and other systematics.

The residual scanning beam tails can also bias the effective beam window function. The spherical harmonic transform of the residual that is not included in the model for the main scanning beam is computed, ${\mathit{B}}_{\mathit{\ell m}}^{\mathrm{full}}\mathrm{=}{\mathit{B}}_{\mathit{\ell m}}^{\mathrm{main}}\mathrm{+}{\mathit{B}}_{\mathit{\ell m}}^{\mathrm{tail}}$, where ${\mathit{B}}_{\mathit{\ell m}}^{\mathrm{main}}$ is the multipole space representation of the main scanning beam model and ${\mathit{B}}_{\mathit{\ell m}}^{\mathrm{tail}}$ is the multipole space representation of the long tail model. In all cases, the m = 0 (symmetric) mode of the ratio of ${\mathit{B}}_{\mathit{\ell m}}^{\mathrm{full}}$ to ${\mathit{B}}_{\mathit{\ell m}}^{\mathrm{main}}$ dominates higher order modes by at least a factor of 1000, meaning that the coupling to the scan strategy is negligible and the bias in the effective beam window function can be approximated by $\mathit{\delta }{\mathit{W}}_{\mathit{\ell }}\mathrm{=}\sum _{\mathit{m}}\frac{\mathrm{|}{\mathit{B}}_{\mathit{\ell m}}^{\mathrm{full}}{\mathrm{|}}^{\mathrm{2}}}{\mathrm{|}{\mathit{B}}_{\mathit{\ell m}}^{\mathrm{main}}{\mathrm{|}}^{\mathrm{2}}}\mathrm{-}\mathrm{1.}$(12)The main effect on the effective beam window function is that at low ℓ, the bias δW_ℓ approaches a value of twice the fractional contribution to the total solid angle. When the window function is normalized to unity at the dipole frequency, the effect is a roughly constant bias in the window function at a level of a few × 10^-4 for ℓ > 100. The contribution of the residual tail to the window function is neglected in the error budget.

3.5. Far sidelobes

The far sidelobes (FSL) are defined as the response of the instrument at angles more than 5° from the main beam centroid. Tauber et al. (2010) describes the pre-launch measurements and predictions of the far sidelobe response using physical optics models. Figure 14 compares the measured beam profile of detector 353-1 with the FSL physical optics model. The way the FSL are treated in the dipole calibration and in the scanning beam model affects the effective beam window function, and care is needed to check whether the off-axis response could bias the window function at ℓ > 40 (angular scales 5°) relative to ℓ< 40. To the extent that the physical optics simulations correctly predict the far sidelobe response (which is shown to be roughly the case in the survey difference maps), they produce effects negligible in the effective beam window function of HFI. Appendix D presents the details of this calculation.

As a check of the quality of the physical optics model of the far sidelobes, Planck Collaboration XIV (2014) attempt a template fit of the physical optics model to the survey difference maps in combination with a zodiacal light model. The template fits are presented in Table 4 of Planck Collaboration XIV (2014). The FSL signature is clearly detected at 857 GHz at a level much higher than predicted. As these channels are multi-moded (Maffei et al. 2010), the differences are not that surprising; it is very difficult to perform the calculations necessary for the prediction. In addition, the specifications for the horn fabrication were quite demanding, and small variations, though still within the mechanical tolerances, could give large variations in the amount of spillover.

For the lower frequency, single-moded channels, there is no clear detection of primary (PR) spillover (i.e., pickup from close to the spacecraft spin axis; see Appendix D). While the significant negative values of the best fit template amplitude may indicate some low-level, large-scale systematic, there seems to be nothing with the distinctive signature of PR spillover at frequencies between 100 and 353 GHz. These values indicate that the PR spillover values in Table 2 of Tauber et al. (2010) may be overestimated.

For the direct contribution of the secondary (SR) spillover, the situation is similar at 353 GHz, but at 217 and 143 GHz there is a 3σ detection at about the level expected, while at 100 GHz the value is about 2.5 times higher than expected, though the signal-to-noise level of the detection is less than 2σ. The baffle contribution to the SR spillover appears higher than predicted, though Planck Collaboration XIV (2014) note that the baffle spillover is difficult to model and the diffuse signal is easily contaminated by other residuals in the survey difference.

4. Effective beams

Unlike WMAP (Jarosik et al. 2011), for large portions of the sky the scan strategy of Planck does not azimuthally symmetrize the effect of the beams on the CMB map. Treating the beams as azimuthally symmetric leads to a flawed power spectrum reconstruction. To remedy this, the effective beam takes the coupling between the azimuthal asymmetry of the beam and the uneven distribution of scanning angles across the sky into account.

The effective beam is computed for each HFI frequency scanning beam and scan history in real space using the FEBeCoP (Mitra et al. 2011) method, as in Planck’s early release (Planck HFI Core Team 2011b). A companion paper describes the details of the application of FEBeCoP to Planck data (Rocha et al., in prep.).

FEBeCoP calculates the effective beam at a position in the sky by computing the real space average of the scanning beam over all crossings angles of that sky position. The observed temperature sky $\begin{array}{c}{\mathrm{􏽥}}_{}\\ {T}\end{array}$ is a convolution of the true sky T and the effective beam ℬ, $\begin{array}{c}\mathrm{􏽥}\\ {T}\end{array}\mathrm{=}{\mathrm{\Omega }}_{\mathrm{pix}}\hspace{0.17em}\mathrm{ℬ}\mathrm{\otimes }{T}\mathit{,}$(13)where Ω_pix is the solid angle of a pixel, and the effective beam can be written in terms of the pointing matrix A_ti and the scanning beam as ${\mathrm{ℬ}}_{\mathit{ij}}\mathrm{=}\frac{{\sum }_{\mathit{t}}{\mathit{A}}_{\mathit{ti}}\hspace{0.17em}\mathit{P}\mathrm{\left(}\hat{{r}}\mathit{j}\mathit{,}\hat{{p}}\mathit{t}\mathrm{\right)}}{{\sum }_{\mathit{t}}{\mathit{A}}_{\mathit{ti}}}\mathit{,}$(14)where t is the time-ordered data sample index and i is the pixel index. A_ti is 1 if the pointing direction falls in pixel number i, else it is 0, p_t represents the pointing direction of the time-ordered data sample, and $\hat{{r}}{}_{\mathit{j}}$ is the centre of pixel number j, where the scanning beam $\mathit{P}\mathrm{\left(}\hat{{r}}{}_{\mathit{j}}\mathit{,}\hat{{p}}{}_{\mathit{t}}\mathrm{\right)}$ is being evaluated (if the pointing direction falls within the cut-off radius).

The sky variation of the effective beam solid angle and the ellipticity of the best-fit Gaussian to the effective beam at HEALPix N_side = 16 pixel centres are shown for 100 GHz in Fig. 15. The effective beam is less elliptical near the ecliptic poles, where more scanning angles symmetrize the beam.

The mean and RMS variation of the effective beam solid angle across the sky for each HFI map are presented in Table 3.

4.1. Effective beam window functions

The multipole space window function of one (or two) observed map(s) is defined, in the absence of instrumental noise and other systematics, as the ratio of the ensemble averaged auto- or cross-power spectrum of the map(s) to the true theoretical sky power spectrum ${\mathit{W}}_{\mathit{\ell }}^{\mathrm{eff}}\mathrm{=}\mathrm{⟨}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{obs}}\mathrm{⟩}\mathit{/}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{sky}}\mathit{.}$(15)It must account for the azimuthal asymmetry of the scan history and the beam profile. This is done in the HFI data processing pipeline using the harmonic space method Quickbeam (described in Appendix E.2), which allows a quick determination of the nominal effective beam window functions and of their Monte Carlo based error eigenmodes (Sect. 6) for all auto- and cross-spectra pairs of HFI detectors, the error determination being computationally intractable with FEBeCop.

In the FEBeCoP approach, many (approximately 1000) random realizations of the CMB sky are generated starting from a given fiducial power spectrum ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{in}}$. For each beam model, the maps are convolved with the pre-computed effective beams, and the pseudo-power spectra ${\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{C}}^{\mathrm{conv}}\end{array}}_{\mathit{\ell }}$ of the resulting maps are computed and corrected by the mode coupling kernel matrix M (Hivon et al. 2002) for a given sky mask: ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{obs}}\mathrm{=}{\mathit{M}}_{\mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}}^{-1}{\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{C}\end{array}}_{{\mathit{\ell }}^{\mathrm{\prime }}}^{\mathrm{conv}}\mathit{.}$ The Monte Carlo average of kernel-corrected power spectra compared to the input power spectrum then gives the effective beam window function (Eq. (15), with ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{in}}$ replacing ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{sky}}$).

In addition, another harmonic space method (FICSBell; Appendix E.1) was also tested and all three methods give consistent results for the nominal window functions at 100–353 GHz (see Fig. E.1).

For two different maps obtained with different detectors or combination of detectors X and Y, because of the optical beam non-circularity and Planck’s scanning strategy, the cross-beam window function is not the geometrical mean of the respective auto-beam window functions, i.e., ${\mathit{W}}^{\mathit{XY}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\ne {\mathrm{\left[}{\mathit{W}}^{\mathit{XX}}\mathrm{\right(}\mathit{\ell }\mathrm{\left)}{\mathit{W}}^{\mathit{YY}}\mathrm{\right(}\mathit{\ell }{\mathrm{\left)}}^{\mathrm{\right]}}}^{\mathrm{1}\mathit{/}\mathrm{2}} \mathrm{if}\mathit{X}\ne \mathit{Y,}$(16)as illustrated in Fig. 16, while of course W^XY = W^YX for any X and Y.

The effective beam window functions for the 2013 maps are shown in Fig. 17. Variations in the effective beam window function from place to place on the sky are significant; the published window functions have been appropriately weighted for the analysis of the nominal mission on the full sky. Analyses requiring effective beam data for more restricted data ranges or for particular regions of the sky should refer to the specialized tools provided in Planck Collaboration (2013).

5. Uncertainties and robustness

Ensembles of simulated planetary observations are used to propagate noise and other systematic effects in the time-ordered data into errors in the beam reconstruction. These simulations are also used to estimate bias in the reconstruction of the beam by comparing the ensemble average beam with the input beam.

To assess the completeness of the statistical and systematic noise model, the consistency of the beam reconstruction derived from different planetary observations is measured against corresponding Monte Carlo ensembles. The MC-derived beam statistics, including both the bias and the correlation structure of the errors, have been found to be somewhat sensitive to the near lobe structure that is included in the beam used as an input to the simulations. An investigation of this effect is ongoing.

In addition, several sources of systematic error that can potentially impact the beam determination, but are not included in the current simulation pipeline, have been investigated. The most significant of these are the beam colour correction and distortion due to ADC non-linearity. These errors are estimated separately and found to be small in comparison to the error bars estimated in the eigenmode analysis described in Sect. 6.

5.1. Simulation of planet observations

The bias and uncertainty in the scanning beams are determined using ensembles of simulated planet observations. The simulation pipeline uses the LevelS Planck simulation code (Reinecke et al. 2006) to generate simulations of the first two observations of Mars and the first three observations of Jupiter and Saturn for each bolometer. This pipeline is used to generate 100 simulations for each bolometer at 353 GHz, 545 GHz and 857 GHz, 200 at 143 GHz and 217 GHz, and 400 at 100 GHz. Each simulation is reconstructed into a beam model using the identical procedure and software as for the real data.

The components of the simulations are as follows.

• 1.

  The scanning beams used as input to the simulations are thesame Mars beams used to calculate the effective beam windowfunction. As shown in Appendix B, thereconstruction bias depends more on the beam representationthan the exact input beam used in a simulation.

• 2.

  The Planck spacecraft pointing and the Horizons ephemerides are used to calculate the pointing relative to the planets for the simulation. An additional random 2.̋5 RMS per sample pointing error is added in both the in-scan and cross-scan directions with a power spectral density given by the pointing error estimate, consistent with the estimated error in the spacecraft pointing (see Sect. 4, and Fig. 7 of Planck Collaboration VI 2014). Error in the beam centroid determination is not simulated.

• 3.

  Detector noise is generated in the timeline by a random realization of Gaussian noise with a power spectral density as reported in Planck Collaboration VI (2014).

For a small number of the simulations, cosmic ray glitches are added to the simulation with the energy spectrum and rate measured in the data (Planck Collaboration X 2014), and the deglitching procedure from the standard pipeline is applied to detect samples contaminated by glitch transients and to subtract the long tails (Planck Collaboration VI 2014).

Each simulation in the ensemble provides an estimate of the scanning beam, which is then used as input to Quickbeam (Appendix E.2) to estimate an effective beam window function (EBWF). The resulting ensemble of EBWFs are then used to compute the bias and errors on the EBWF derived from the data; the mean of the ensemble provides a measure of the reconstruction bias, and the distribution of the ensemble gives the uncertainty. The bulk of the RMS of the ensemble can be captured with a small number of eigenmodes (Sect. 6).

This procedure allows for the direct determination of the EBWF and associated errors for each of the detector cross-correlations input to the angular power spectrum likelihood, described in Planck Collaboration XV (2014).

5.2. Absolute consistency: comparison of systematics against simulations

5.2.1. Seasonal consistency

One important test of the consistency of the scanning beam measurement is the stability over time, as measured in different seasons and on different objects over the course of the mission. The difference of the scanning beam window function between the first two observations of Mars is shown in Fig. 18. The residuals are well within the estimated uncertainty that is derived from the simulation ensemble for 100–353 GHz. Due to the sparse cross-scan sampling, the B-spline beam representation does not converge for a single planetary observation; for this test the Gauss-Hermite (GH) formalism is used.

5.2.2. Time variability of planet flux density

Mars is known to have diurnal variability at HFI frequencies due to the non-uniformity of the surface albedo⁵. A single detector scans Mars on time scales of a few hours, during which the brightness temperature of Mars can vary by several percent. Additionally, the Mars planetary disc area varies by several percent during an observation as it moves relative to Planck. To assess the impact of Mars variability on the scanning beam determination, GH representations of the beam are derived for each bolometer, with and without a model for the diurnal variability of Mars. Using χ² tests for goodness of fit, these results are indistinguishable.

5.2.3. Beam colour correction

The scanning beams are measured using planets, whose spectral energy distribution (SED) roughly follows the Rayleigh-Jeans spectrum. However, the beam window functions from these measurements are used to correct the angular power spectra of CMB anisotropies, which is characterized by a very different spectrum, one that is falling in power as a function of frequency relative to a Rayleigh-Jeans spectrum across the HFI bands.

Physical optics simulations, using the GRASP⁶ software, are used to investigate this effect. For each HFI channel from 100–353 GHz, monochromatic simulations are generated at five frequencies across each band using the pre-launch telescope model (Maffei et al. 2010; Tauber et al. 2010). The solid angle of these simulated beams varies across the band, due to diffraction and focussing effects. For 100–217 GHz, the beam size comes to a minimum near the band centre, while at 353 GHz the beam size rises with frequency.

No attempt is made to colour-correct the planet-derived window functions, because a telescope solution has not yet been determined that agrees with the measured solid angles. Spot checks at 143 and 353 GHz with a defocussed telescope model improve agreement between the data and the model, but does not change the trend of solid angle with frequency across the band. The investigation of these effects suggests that numerical models can bound the uncertainty, but cannot reliably predict a bias. Rayleigh-Jeans-weighted average and CMB anisotropy-weighted average beams are produced and used to compute an effective beam window function (blue curves in Fig. 19). The largest bias results at 353 GHz. For the frequencies 100, 143, and 217 GHz, however, the band average effect is less than 0.1% across the entire multipole range used in the CMB likelihood.

The beam solid angle also varies as a function of source SED; the GRASP simulations constrain the size of the beam solid angle colour correction from a ν² source SED (like the planets) to the IRAS convention ν^-1. At 545 and 857 GHz, the GRASP models of Murphy et al. (2010) provide a measure of the effect, which is found to be significant at 353 GHz, but not at the other frequencies, at the level of a few tenths of a percent (Table 2).

5.2.4. ADC non-linearity

Non-linearity in the ADC in the HFI readout electronics mainly manifests itself as an apparent gain drift (Planck Collaboration VI 2014; Planck Collaboration VIII 2014), however, the measured beam shapes are also biased by the effect.

Correcting the ADC non-linearity relies on a model that predicts which ADC codes contribute to each data sample. The presence of large signal gradients, such as a planet, or pickup from the ⁴He-JT cooler, make modelling difficult. A model for correcting every detector is still under development.

A model of the behaviour of the ADCs is used to apply non-linearity to simulated Mars observations. The B-spline scanning beams reconstructed from these simulations predict biases that are typically under 2% at ℓ ≤ 2000, less than the RMS error (see the magenta curves in Fig. 19 for simulated bias of the 100 through 353 GHz channels).

Using the brighter planets in future Planck main beam models will tend to reduce the effect, as the higher signal amplitudes sample a broader range of the ADC, tending to average down the resulting bias (Mather et al. 1993).

5.2.5. Pointing errors

While the simulated planet observations account for random pointing errors, pointing drifts remain at the level of a few arcseconds per pointing period in the co-scan and cross-scan directions common to all detectors (Planck Collaboration I 2014; Planck Collaboration VI 2014). These have the effect of widening the beam. Even in the worst case, if remaining errors are Gaussian-distributed with a 2.̋5 residual, the window function ${\mathit{B}}_{\mathit{\ell }}^{\mathrm{2}}$ is 0.2% biased at multipole ℓ = 3000. Since this is a very small effect relative to the other beam errors, it is considered negligible.

5.2.6. Glitches

The HFI data are affected by transient signals due to interactions with high energy particles. The planet simulation tools allow for the addition of a random population of glitches to simulated Mars observations with a realistic rate and energy spectrum (Planck Collaboration X 2014), testing the performance of the deglitching algorithm near planet signals and the effects of undetected low energy glitches in the channel with the highest glitch rate at each frequency band.

The RMS of the noise increases by at most 20% in the bolometer with the highest glitch rate, but more typically less than 10%, and the reconstruction bias changes by a negligible amount. The effects are small enough that they are not included in the error budget.

5.3. Relative consistency: window function ratios vs. spectral ratios

The CMB sky itself allows an additional, and statistically powerful, check on the relative consistency (but not the absolute accuracy) of the effective beam window functions within a frequency band. See Fig. 33, and the associated discussion in Sect. 7.3 of Planck Collaboration VI (2014), showing the self-consistency of the window functions at a level better than 0.4%, across the full range of multipoles employed in cosmological parameter estimation.

6. Total error budget

The distribution of beam solid angles reconstructed from simulated planet reconstructions provides the statistical error budget for the beam solid angles (Table 2). Other systematic effects are small in comparison.

The uncertainty in the effective beam window function is determined with the distribution of window functions computed from the simulated beams. These errors are highly correlated in multipole space, and we find that they can be fully described by a small number of eigenmodes and their covariance matrix (Appendices G.1 and G.2), which can be retrieved, together with the nominal beam window functions, from Planck Collaboration (2013) and the Planck Legacy Archive⁷. These simulations also allow us to determine the bias induced by the beam reconstruction pipeline, and correct the final beam window functions from it (Appendix G.3).

Studies have been performed to probe the impact of assumptions regarding the near sidelobe response of the beam. These studies suggest the potential that the Monte Carlo ensembles used to derive the error eigenmodes may not capture aspects of the bias and the correlation structure of the errors. A visual summary of the total error budget, including the impact of the near sidelobe response, is provided in Fig. 19. For the initial data release a conservative scaling of the MC-derived beam errors, a factor of 2.7 in power, is applied in order to accommodate the maximum extent of the bias that is allowable, given the signal-to-noise ratio of the Jupiter and Saturn data.

7. Conclusion

A combination of Jupiter and Saturn observations, survey difference maps, and checkout and performance verification (CPV) phase data is used to derive the bolometer/readout electronics time response transfer function. This transfer function is deconvolved from the HFI time-ordered data. Residuals can be detected as long tails, stretching up to 6° from compact sources, but are at a level less than 10^-4 of the beam peak, and so are insignificant for the scanning beams.

The effective beams and window functions are estimated for all HFI detectors using a scanning beam derived from a combination of the first two seasons of Mars observations. The effective beam products account for the partial symmetrization of the scanning beam which results from the scan strategy. The FEBeCoP algorithm produces effective beams in real space, while the Quickbeam method is used to produce the effective beam window functions and errors in harmonic space.

The final error budget for each effective beam window function (both auto- and cross-beams) is well-described by five eigenmodes for each beam and a correlation matrix of eigenvalues. The error parameters are estimated with an ensemble of simulated Mars observations. The simulations are based on random realizations of noise and pointing errors.

No significant time variation of the scanning beam or the time response is found. Cross-power window functions are also produced that describe the beam filtering effects in cross-correlations between independent HFI maps.

The high signal-to-noise ratio of the Jupiter and Saturn data limit the contribution of the near sidelobe response (between 30′ and 5° of the beam centroid) to the total beam solid angle to less than 0.2% at 100 to 217 GHz. This portion of the beam, which is beyond the signal-to-noise of the Mars data, is not fully represented in the beam simulations. The Monte Carlo derived error amplitudes are scaled by a factor of 2.7 to account for the contribution of the near sidelobe response to the beam window function.

Far sidelobes contribute negligibly to the window function, but may potentially impact the primary calibration and result in pickup from bright Galactic sources. Improved physical optics models for the far sidelobes will be forthcoming for future releases.

The impact of non-linearity in the ADC, the sensitivity of beam shape to the spectral intensity of the source, and residual cosmic ray transients are found to be insignificant sources of systematic error.

Knowledge of the beam window functions allows Planck HFI to extrapolate the dipole calibration over more than three orders of magnitude in angular scale. For the current data release, beam errors are subdominant to noise, foreground marginalization and sample variance in cosmological parameter estimation. Future Planck data releases will fully exploit all of the available planetary data and create wide-field beam maps, allowing an even more precise, and accurate, measurement of the main beams and near sidelobes.

---

Acknowledgments

The development of Planck has been supported by: ESA; CNES and CNRS/INSU-IN2P3-INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MICINN, 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); and PRACE (EU). A description of the Planck Collaboration and a list of its members, including the technical or scientific activities in which they have been involved, can be found at http://www.sciops.esa.int/index.php?project=planck&page=Planck_Collaboration.

References

Appendix A: Electronics model

This Appendix derives the effect of the HFI electronics filtering on the TOI. If the input signal (power) on a bolometer is ${\mathit{s}}_{\mathrm{0}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathrm{e}}^{\mathrm{i}\mathit{\omega t}}\mathit{,}$(A.1)the bolometer physical impedance can be written as: $\mathit{s}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathrm{e}}^{\mathrm{i}\mathit{\omega t}}\mathit{F}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathit{,}$(A.2)where ω is the angular frequency of the signal and F(ω) is the complex intrinsic bolometer transfer function. For HFI the bolometer transfer function is modelled as the sum of 4 single pole low pass filters, $\mathit{F}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}\mathit{,}\mathrm{4}}\frac{{\mathit{a}}_{\mathit{i}}}{\mathrm{1}\mathrm{+}\mathrm{i}\mathit{\omega }{\mathit{\tau }}_{\mathit{i}}}\mathrm{·}$(A.3)The modulation of the signal is done with a square wave, written here as a composition of sine waves of decreasing amplitude, ${\mathit{s}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathrm{e}}^{\mathrm{i}\mathit{\omega t}}\mathit{F}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\sum _{\mathit{k}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\frac{{\mathrm{e}}^{\mathrm{i}{\mathit{\omega }}_{\mathit{r}}\mathrm{(}\mathrm{2}\mathit{k}\mathrm{+}\mathrm{1}\mathrm{)}\mathit{t}}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{\mathit{\omega }}_{\mathit{r}}\mathrm{(}\mathrm{2}\mathit{k}\mathrm{+}\mathrm{1}\mathrm{)}\mathit{t}}}{\mathrm{2}\mathrm{i}\mathrm{(}\mathrm{2}\mathit{k}\mathrm{+}\mathrm{1}\mathrm{)}}\mathit{,}$(A.4)where the Euler relation sinx = (e^ix−e^−ix) /2i is used, and ω_r is the angular frequency of the square wave. The modulation frequency is f_mod = ω_r/2π and was set to f_mod = 90.18759 Hz in flight. This signal is then filtered by the complex electronic transfer function H(ω). Setting $\begin{array}{ccc}{\mathit{\omega }}_{\mathit{k}}^{\mathrm{+}}\mathrm{=}\mathit{\omega }\mathrm{+}\mathrm{(}\mathrm{2}\mathit{k}\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{\omega }}_{\mathit{r}}& & \end{array}$and $\begin{array}{ccc}{\mathit{\omega }}_{\mathit{k}}^{\mathrm{-}}\mathrm{=}\mathit{\omega }\mathrm{-}\mathrm{(}\mathrm{2}\mathit{k}\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{\omega }}_{\mathit{r}}& & \end{array}$gives $\mathrm{\Sigma }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\frac{\mathit{F}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}}{\mathrm{2}\mathrm{i}\mathrm{(}\mathrm{2}\mathit{k}\mathrm{+}\mathrm{1}\mathrm{)}}\mathrm{\left[}\mathit{H}\mathrm{\right(}{\mathit{\omega }}_{\mathit{k}}^{\mathrm{+}}\mathrm{)}{\mathrm{e}}^{\mathrm{i}{\mathit{\omega }}_{\mathit{k}}^{\mathrm{+}}\mathit{t}}\mathrm{-}\mathit{H}\mathrm{(}{\mathit{\omega }}_{\mathit{k}}^{\mathrm{-}}\mathrm{\left)}{\mathrm{e}}^{\mathrm{i}{\mathit{\omega }}_{\mathit{k}}^{\mathrm{-}}\mathit{t}}\mathrm{\right]}\mathit{.}$(A.5)This signal is then sampled at high frequency (2f_modN_samp). N_samp is a parameter of the HFI electronics corresponding to the number of high frequency samples in each modulation semi-period. In order to obtain an output signal sampled every π/ω_r seconds, the waveform is integrated on a semi-period, as done in the HFI readout. To also include a time shift Δt, the integral is calculated between nπ/ω_r + Δt and (n + 1)π/ω_r + Δt (with T = 2π/ω_r period of the modulation). The time shift Δt is encoded in the HFI electronics by the parameter S_phase, with the relation Δt = S_phase/N_samp/f_mod.

After integration, a sample of a bolometer at time t_n can be written as $\mathit{Y}\mathrm{\left(}{\mathit{t}}_{\mathit{n}}\mathrm{\right)}\mathrm{=}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{n}}\mathit{F}\mathrm{(}\mathit{\omega }\mathrm{\left)}{\mathit{H}}^{\mathrm{\prime }}\mathrm{\right(}\mathit{\omega }\mathrm{)}{\mathrm{e}}^{\mathrm{i}{\mathit{t}}_{\mathit{n}}\mathit{\omega }}\mathit{,}$(A.6)where $\begin{array}{ccc}{\mathit{H}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\sum _{\mathit{k}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathrm{(}\frac{\mathit{\pi \omega }}{\mathrm{2}{\mathit{\omega }}_{\mathit{r}}}\mathrm{+}\mathit{\omega }\mathrm{\Delta }\mathit{t}\mathrm{)}}[\frac{\mathit{H}\mathrm{\left(}{\mathit{\omega }}_{\mathit{k}}^{\mathrm{+}}\mathrm{\right)}{\mathrm{e}}^{\mathrm{i}{\mathit{\omega }}_{\mathit{k}}^{\mathrm{+}}\mathrm{\Delta }\mathit{t}}}{\mathrm{(}\mathrm{2}\mathit{k}\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{\omega }}_{\mathit{k}}^{\mathrm{+}}}\left(\mathrm{1}\mathrm{-}{\mathrm{e}}^{\frac{\mathrm{i}{\mathit{\omega }}_{\mathit{k}}^{\mathrm{+}}\mathit{\pi }}{{\mathit{\omega }}_{\mathit{r}}}}\right)& & \\ \mathrm{-}\frac{\mathit{H}\mathrm{\left(}{\mathit{\omega }}_{\mathit{k}}^{\mathrm{-}}\mathrm{\right)}{\mathrm{e}}^{\mathrm{i}{\mathit{\omega }}_{\mathit{k}}^{\mathrm{-}}\mathrm{\Delta }\mathit{t}}}{\mathrm{(}\mathrm{2}\mathit{k}\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{\omega }}_{\mathit{k}}^{\mathrm{-}}}\left(\mathrm{1}\mathrm{-}{\mathrm{e}}^{\frac{\mathrm{i}{\mathit{\omega }}_{\mathit{k}}^{\mathrm{-}}\mathit{\pi }}{{\mathit{\omega }}_{\mathit{r}}}}\right)& & ]\mathit{.}\end{array}$(A.7)The output signal in Eq. (A.6) can be demodulated (thus removing the (− 1)ⁿ) and compared to the input signal in Eq. (A.1). The overall transfer function is composed of the bolometer transfer function and the effective electronics transfer function, H′(ω), $\mathit{TF}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{=}\mathit{F}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}{\mathit{H}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathit{.}$(A.8)The shape of H(ω) is obtained by combining low and high-pass filters with Sallen Key topologies (Sallen & Key 1955), with their respective time constants, and accounting also for the stray capacitance low pass filter given by the bolometer impedance combined with the stray capacitance of the cables. The sequence of filters that define the electronic band-pass function H(ω) = h₀ ∗ h₁ ∗ h₂ ∗ h₃ ∗ h₄ ∗ h₅ are listed in Table A.1, with parameters listed in Table A.2.

Appendix B: Gauss-Hermite beams

HFI’s scanning beams are described by an elliptical Gaussian shape to an accuracy of 2–5% in solid angle. A Gauss-Hermite representation of the beam uses Hermite polynomials to provide higher-order corrections to a Gaussian model (Huffenberger et al. 2010). There are fewer degrees of freedom than in a gridded beam map, allowing higher signal-to-noise on each mode. However, because the order of the decomposition is truncated, in practice the description is ill suited to a description of features much beyond the extent of the main beam. Larger scale features of the beam, including a beam shoulder (Ruze 1966) or the effect of the print-through of the backing structure (grating lobes), must be included separately.

The initial two-dimensional Gaussian is parametrized as $\mathit{P}\mathrm{\left(}{\mathit{x}}_{\mathrm{1}}\mathit{,}{\mathit{x}}_{\mathrm{2}}\mathrm{\right)}\mathrm{=}\frac{\mathit{A}}{\mathrm{|}\mathrm{2}\mathit{\pi }{\Sigma }{\mathrm{|}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\exp \left[\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\sum _{\mathit{i,j}\mathrm{=}\mathrm{1}}^{\mathrm{2}}\mathrm{(}{\mathit{x}}_{\mathit{i}}\mathrm{-}{\mathit{x̅}}_{\mathit{i}}\mathrm{)}{\mathrm{\Sigma }}_{\mathit{ij}}^{-1}\mathrm{(}{\mathit{x}}_{\mathit{j}}\mathrm{-}{\mathit{x̅}}_{\mathit{j}}\mathrm{)}\right]\mathit{,}$(B.1)where A is an overall amplitude, (x₁,x₂) are two-dimensional Cartesian coordinates (the pointing is projected to a flat sky), are the coordinates of the beam centre, and the correlation matrix is given by ${\Sigma }\mathrm{=}\left(\begin{array}{c}\\ {\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}& \mathit{\rho }{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}\\ \mathit{\rho }{\mathit{\sigma }}_{\mathrm{1}}{\mathit{\sigma }}_{\mathrm{2}}& {\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}\end{array}\right)\mathit{.}$(B.2)Hence, the Gaussian model parameters are fitted. These can also be expressed in terms of the ellipticity, e (defined here as the ratio between the major and minor axes), and rotation angle α, of the Gaussian ellipse.

The coefficients to the Gauss-Hermite polynomials multiplying the elliptical Gaussian are then fitted. The basis functions are defined as ${\mathrm{\Phi }}_{{\mathit{n}}_{\mathrm{1}}{\mathit{n}}_{\mathrm{2}}}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{\propto }{\mathit{H}}_{{\mathit{n}}_{\mathrm{1}}}\mathrm{\left(}{\mathit{x}}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{\right)}{\mathit{H}}_{{\mathit{n}}_{\mathrm{2}}}\mathrm{\left(}{\mathit{x}}_{\mathrm{2}}^{\mathrm{\prime }}\mathrm{\right)}\exp \left(\mathrm{-}{{x}}^{{\prime }}\mathrm{·}{{x}}^{{\prime }}\mathit{/}\mathrm{2}\right)\mathit{,}$(B.3)where H_n(x) is the order-n Hermite polynomial and the primed coordinates x^′ rotates into a system aligned with the axes of the Gaussian and scaled to the major and minor axes σ_i (i.e., to the principle axes of the correlation matrix Σ). Having determined the elliptical Gaussian separately, the subsequent Gauss-Hermite polynomial fit is a generalized least-squares procedure, solvable by the usual matrix techniques, under the assumption of white noise. Because the full data model includes further effects such as pointing error, glitch effects, etc, a full error analysis requires a broad suite of simulations, described in Sect. 5.

Appendix C: B-spline beams

In this model of the scanning beam, a two-dimensional B-spline surface is fit to the time-ordered planet data. A smoothing criterion is applied during the fit to minimize the effects of high spatial frequency variations due to noise. This representation of the beam is superior to a simple two-dimensional binning of the data in its ability to capture large signal gradients.

B-splines are a linear combination of polynomials of degree k and order k + 1. They are defined by the location of their control points (called knots) of which there are 5 for 3rd degree polynomials.

The k-degree B-spline built using the knots {λ_i,...,λ_{i + k + 1}} (de Boor 1972; Cox 1972) is given by ${\mathit{P}}_{\mathit{i,}\mathrm{1}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ \mathrm{1}\mathit{,}\mathrm{if}\mathit{x}\mathrm{\in }\mathrm{[}{\mathit{\lambda }}_{\mathit{i}}\mathit{,}{\mathit{\lambda }}_{\mathit{i}\mathrm{+}\mathrm{1}}\mathrm{]}\\ \mathrm{0}\mathit{,}\mathrm{if}\mathit{x}\notin \mathrm{[}{\mathit{\lambda }}_{\mathit{i}}\mathit{,}{\mathit{\lambda }}_{\mathit{i}\mathrm{+}\mathrm{1}}\mathrm{]}\end{array}$(C.1)with recursion relations ${\mathit{P}}_{\mathit{i,l}\mathrm{+}\mathrm{1}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{\mathit{x}\mathrm{-}{\mathit{\lambda }}_{\mathit{i}}}{{\mathit{\lambda }}_{\mathit{i}\mathrm{+}\mathit{l}}\mathrm{-}{\mathit{\lambda }}_{\mathit{i}}}{\mathit{P}}_{\mathit{i,l}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{+}\frac{{\mathit{\lambda }}_{\mathit{i}\mathrm{+}\mathit{l}\mathrm{+}\mathrm{1}}\mathrm{-}\mathit{x}}{{\mathit{\lambda }}_{\mathit{i}\mathrm{+}\mathit{l}\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{\lambda }}_{\mathit{i}\mathrm{+}\mathrm{1}}}{\mathit{P}}_{\mathit{i}\mathrm{+}\mathrm{1}\mathit{,l}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathit{.}$(C.2)where the index l runs from 1 through the B-spline degree k. The B-spline knots {λ_i,...,zλ_{i + k + 1}} are located on a regularly spaced grid in the detector coordinate frame. At the edge of the reconstructed beam map area, 4 coincident knots are added to avoid vanishing basis functions, allowing a unique decomposition.

In the solution to the B-spline coefficients P_i,l(x), a smoothing criterion is introduced as a constraint on the sum of the derivatives of the beam at each knot, motivated by the assumption that the true beam does not contain very high spatial frequencies and prevents noise from biasing the reconstruction at spatial frequencies smaller than the smoothing scale.

A smoothing criterion η (Dierckx 1993) is related to the sum of the order k derivative of the beam model (P^k) evaluated at the knot locations λ_i: $\mathit{\eta }\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{g}}{\mathrm{\left[}{\mathit{P}}^{\mathit{k}}\mathrm{\right(}{\mathit{\lambda }}_{\mathit{i}\mathrm{+}}\mathrm{)}\mathrm{-}{\mathit{P}}^{\mathit{k}}\mathrm{(}{\mathit{\lambda }}_{\mathit{i}\mathrm{-}}{\mathrm{\left)}}^{\mathrm{\right]}}}^{\mathrm{2}}$(C.3)where g is the total number of distinct knots and λ_{i +} and λ_{i −} are the left and right derivative of the beam model evaluated at the knot location. The smoothing criterion is introduced as an extra term in the score function ζ in the least-squares minimization of the beam map with respect to the data, $\mathit{\zeta }\mathrm{=}\mathit{\eta }\mathrm{+}\mathit{p}\mathrm{\times }\mathit{\delta }\mathit{.}$(C.4)where δ is the usual squared difference between the data points y_r and the model P(x_r) at pointings x_r, $\mathit{\delta }\mathrm{=}\sum _{\mathit{r}\mathrm{=}\mathrm{1}}^{\mathit{m}}{\left[{\mathit{y}}_{\mathit{r}}\mathrm{-}\mathit{P}\mathrm{\left(}{\mathit{x}}_{\mathit{r}}\mathrm{\right)}\right]}^{\mathrm{2}}\mathit{.}$(C.5)and p is a relative weighting factor for the smoothing criterion. The knot spacing and the smoothing criterion weight p are determined separately for each frequency band based on physical optics simulations and the coverage of the two Mars transits given by the scanning strategy (see Table C.1). Simulated planet observations (see Sect. 5) show that the choice of B-spline knot spacing and smoothing parameters do not significantly bias the beam reconstruction.

Appendix D: Far sidelobe effects on the effective beam window function

The far sidelobes (FSL) are defined as the response of the instrument at angles >5° from the main beam centroid. The FSL affect both 1. the gain calibration of the instrument with the dipole, and 2. how the calibration is transfered to higher multipoles with the effective beam window function. Here we discuss the interplay of these effects.

The FSL can be separated into three main components (see Fig. 5 of Tauber et al. 2010):

• 1.

  Primary Reflector Spillover (PR Spillover) is the response ofthe instrument to radiation from just above the primary mirror thatreflects off the secondary mirror and arrives at the feed horns. Thisresponse is nearly aligned with the spin axis of the spacecraft, andtherefore scans very little of the sky on 1 min time-scales.

• 2.

  Secondary Reflector Direct Spillover(SR Direct Spillover, or SRD) is the response from directly above the secondary mirror. The sidelobe peaks roughly 10° from the main beam, and as such scans the sky in nearly the same way as the main beam as the satellite rotates. The azimuthal extent is roughly 30°.

• 3.

  Secondary Reflector Baffle Spillover (SR Baffle Spillover, or SRB) is response from radiation reflecting off the baffles. This is difficult to model, being diffuse radiation reflecting off the poorly known inner baffle surfaces. It is spread over a large fraction of the sky.

The HFI dipole calibration is performed assuming a delta-function (pencil) beam (Planck Collaboration VIII 2014). This leads to a bias in the calibration described by the ratio of the dipole convolved with the full-sky beam to the dipole convolved by a pencil beam, $\begin{array}{ccc}\begin{array}{c}\mathrm{˜}\\ \mathit{g}\end{array}\mathrm{=}\mathit{g}\frac{\mathit{P}\mathrm{\otimes }\mathit{D}}{{\mathit{P}}^{\mathrm{pencil}}\mathrm{\otimes }\mathit{D}}\mathit{,}& & \end{array}$where $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{g}\end{array}$ is the estimated gain, g is the true gain, P is the true full-sky beam, P^pencil is the assumed pencil beam, and D is the dipole. The true full-sky beam is taken to consist of three portions, $\begin{array}{ccc}\mathit{P}\mathrm{=}{\mathit{P}}_{\mathrm{main}}\mathrm{+}{\mathit{P}}_{\mathrm{NSL}}\mathrm{+}{\mathit{P}}_{\mathrm{FSL}}\mathit{,}& & \end{array}$the main beam P_main (within 20′ of the centroid), the near sidelobes P_NSL (between 20′ and 5°), and the far sidelobes P_FSL (response further than 5° from the centroid). The quantity $\begin{array}{ccc}{\mathit{P}}_{\mathrm{MNSL}}\mathrm{=}{\mathit{P}}_{\mathrm{main}}\mathrm{+}{\mathit{P}}_{\mathrm{NSL}}& & \end{array}$is defined as the beam within 5° of the centroid (the main lobe and near sidelobes). The bias in the calibration can be rewritten as $\begin{array}{c}\mathrm{˜}\\ \mathit{g}\end{array}\mathrm{=}\mathit{g}\left(\mathrm{1}\mathrm{-}{\mathit{f}}_{\mathrm{FSL}}\mathrm{+}\frac{{\mathit{P}}_{\mathrm{FSL}}\mathrm{\otimes }\mathit{D}}{{\mathit{P}}^{\mathrm{pencil}}\mathrm{\otimes }\mathit{D}}\right)\mathit{,}$(D.1)where f_FSL is the integral of the far sidelobe response relative to the full beam integral. The first term, 1 − f_FSL, is due to the loss in response of the main lobe and near sidelobes to far sidelobes, while the second term represents the coupling of the dipole into the sidelobes.

For each pointing period the convolution of the pencil beam with the dipole is approximated by ${\mathit{P}}^{\mathrm{pencil}}\mathrm{\otimes }\mathit{D}\mathrm{\simeq }\sin {\mathit{\theta }}_{\mathrm{main}}{\mathit{d}}_{\mathrm{pencil}}\mathit{,}$(D.2)where θ_main is the angle between the main beam centroid and the spin axis, and d_pencil is the dipole amplitude that would be observed with a pencil beam. The components of the far sidelobes enter into Eq. (D.1) differently, and assuming the majority of the response to the dipole in each component comes from one direction on the sky, can be approximated as $\begin{array}{ccc}{\mathit{P}}_{\mathrm{FSL}}\mathrm{\otimes }\mathit{D}& \mathrm{\simeq }& {\mathit{ϵ}}_{\mathrm{PR}}\sin {\mathit{\theta }}_{\mathrm{PR}}{\mathit{d}}_{\mathrm{pencil}}\mathrm{+}{\mathit{ϵ}}_{\mathrm{SRD}}\sin {\mathit{\theta }}_{\mathrm{SRD}}{\mathit{d}}_{\mathrm{pencil}}\\ & & \mathrm{+}{\mathit{ϵ}}_{\mathrm{SRB}}\sin {\mathit{\theta }}_{\mathrm{SRB}}{\mathit{d}}_{\mathrm{pencil}}\mathit{,}\end{array}$where ϵ_PR, ϵ_SRD, and ϵ_SRB are the fraction of the total solid angle in the PR spillover, the SR direct spillover, and the SR baffle spillover, respectively, and θ_PR, θ_SRD, θ_SRB are the angles between the spin axis and the peak of the PR spillover response, the SR direct spillover response, and the SR baffle spillover response respectively. Equation (D.1) then simplifies to $\begin{array}{c}\mathrm{˜}\\ \mathit{g}\end{array}\mathrm{\simeq }\mathit{g}\left(\mathrm{1}\mathrm{-}{\mathit{f}}_{\mathrm{FSL}}\mathrm{+}{\mathit{ϵ}}_{\mathrm{PR}}\frac{\sin {\mathit{\theta }}_{\mathrm{PR}}}{\sin {\mathit{\theta }}_{\mathrm{pencil}}}\mathrm{+}{\mathit{ϵ}}_{\mathrm{SRD}}\frac{\sin {\mathit{\theta }}_{\mathrm{SRD}}}{\sin {\mathit{\theta }}_{\mathrm{pencil}}}\mathrm{+}{\mathit{ϵ}}_{\mathrm{SRB}}\frac{\sin {\mathit{\theta }}_{\mathrm{SRB}}}{\sin {\mathit{\theta }}_{\mathrm{pencil}}}\right)\mathrm{·}$(D.3)From Tauber et al. (2010), θ_pencil ≃ 85°, θ_PR ≃ 10°, θ_SRD ≃ 75°, and θ_SRB ≃ 45° making the PR spillover and SR baffle spillover negligible. Equation (D.3) reduces further to $\begin{array}{c}\mathrm{˜}\\ \mathit{g}\end{array}\mathrm{\simeq }\mathit{g}\left(\mathrm{1}\mathrm{-}{\mathit{f}}_{\mathrm{FSL}}\mathrm{+}{\mathit{ϵ}}_{\mathrm{SRD}}\frac{\sin {\mathit{\theta }}_{\mathrm{SRD}}}{\sin {\mathit{\theta }}_{\mathrm{pencil}}}\right)\mathrm{·}$(D.4)SR direct spillover scans the dipole, but with a slightly different amplitude, since it is offset by 10° from the main lobe. The PR spillover does not modulate the dipole; aligned with the spin axis, the PR spillover contributes a nearly constant signal during each stable pointing period. The predicted values of the solid angle fractions are f_FSL = 3.3 × 10^-3 and ϵ_SRD = 1.9 × 10^-3 at 100 GHz (Tauber et al. 2010), making the measured gain $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{g}\end{array}\mathrm{\simeq }\mathrm{0.9985}\mathit{g}$.

However, the effect of the gain bias on the angular power spectrum is further reduced by corrections to the effective beam window function due to the FSL. Not including the FSL in the beam model tends to bias the window function at very low multipoles relative to high multipoles. Considering the measured sky signal at small scales $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{T}}_{\mathrm{sky}}\end{array}$ as compared to the true sky T_sky one has $\begin{array}{ccc}\begin{array}{c}\mathrm{˜}\\ {\mathit{T}}_{\mathrm{sky}}\end{array}\mathrm{=}\frac{\mathit{S}}{\begin{array}{c}\mathrm{˜}\\ \mathit{g}\end{array}}\mathrm{=}\frac{\mathit{P}\mathrm{\otimes }{\mathit{T}}_{\mathrm{sky}}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{f}}_{\mathrm{FSL}}\mathrm{+}{\frac{{\mathit{P}}_{\mathrm{FSL}}\mathrm{\otimes }\mathit{D}}{{\mathit{P}}^{\mathrm{pencil}}\mathrm{\otimes }\mathit{D}}}^{\mathrm{)}}}\mathit{,}& & \end{array}$where S is the signal measured by the bolometer. Solving for the true sky signal gives $\begin{array}{ccc}{\mathit{P}}_{\mathrm{MNSL}}\mathrm{\otimes }{\mathit{T}}_{\mathrm{sky}}\mathrm{=}\left(\mathrm{1}\mathrm{-}{\mathit{f}}_{\mathrm{FSL}}\mathrm{+}\frac{{\mathit{P}}_{\mathrm{FSL}}\mathrm{\otimes }\mathit{D}}{{\mathit{P}}^{\mathrm{pencil}}\mathrm{\otimes }\mathit{D}}\right)\begin{array}{c}\mathrm{˜}\\ {\mathit{T}}_{\mathrm{sky}}\end{array}\mathrm{-}{\mathit{P}}_{\mathrm{FSL}}\mathrm{\otimes }{\mathit{T}}_{\mathrm{sky}}\mathit{.}& & \end{array}$Considering that P_MNSL ⊗ T_sky ≃ (1 − f_FSL)T_sky, $\begin{array}{ccc}{\mathit{T}}_{\mathrm{sky}}\mathrm{\simeq }\left(\mathrm{1}\mathrm{+}\frac{{\mathit{P}}_{\mathrm{FSL}}\mathrm{\otimes }\mathit{D}}{{\mathit{P}}^{\mathrm{pencil}}\mathrm{\otimes }\mathit{D}}\frac{\mathrm{1}}{\mathrm{1}\mathrm{-}{\mathit{f}}_{\mathrm{FSL}}}\right)\begin{array}{c}\mathrm{˜}\\ {\mathit{T}}_{\mathrm{sky}}\end{array}\mathrm{-}\frac{{\mathit{P}}_{\mathrm{FSL}}\mathrm{\otimes }{\mathit{T}}_{\mathrm{sky}}}{\mathrm{1}\mathrm{-}{\mathit{f}}_{\mathrm{FSL}}}\mathrm{·}& & \end{array}$Now also considering that P_MNSL ⊗ D ≃ (1 − f_FSL)P^pencil ⊗ D, the true sky signal can be written as $\begin{array}{ccc}{\mathit{T}}_{\mathrm{sky}}\mathrm{\simeq }\begin{array}{c}\mathrm{˜}\\ {\mathit{T}}_{\mathrm{sky}}\end{array}\left(\mathrm{1}\mathrm{+}\frac{{\mathit{P}}_{\mathrm{FSL}}\mathrm{\otimes }\mathit{D}}{{\mathit{P}}_{\mathrm{MNSL}}\mathrm{\otimes }\mathit{D}}\mathrm{-}\frac{{\mathit{P}}_{\mathrm{FSL}}\mathrm{\otimes }{\mathit{T}}_{\mathrm{sky}}}{{\mathit{P}}_{\mathrm{MNSL}}\mathrm{\otimes }{\mathit{T}}_{\mathrm{sky}}}\frac{{\mathit{T}}_{\mathrm{sky}}}{\begin{array}{c}\mathrm{˜}\\ {\mathit{T}}_{\mathrm{sky}}\end{array}}\right)\mathit{,}& & \end{array}$or ${\mathit{T}}_{\mathrm{sky}}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\phi }}_{\mathit{D}}\mathrm{-}{\mathit{\phi }}_{\mathrm{sky}}\mathrm{)}\begin{array}{c}\mathrm{˜}\\ {\mathit{T}}_{\mathrm{sky}}\end{array}\mathit{,}$(D.5)where the second-order correction term ${\mathit{T}}_{\mathrm{sky}}\mathit{/}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{T}}_{\mathrm{sky}}\end{array}\mathrm{\simeq }\mathrm{1}$ is dropped, and the following factors are defined: $\begin{array}{ccc}{\mathit{\phi }}_{\mathit{D}}\mathrm{=}\frac{{\mathit{P}}_{\mathrm{FSL}}\mathrm{\otimes }\mathit{D}}{{\mathit{P}}^{\mathrm{pencil}}\mathrm{\otimes }\mathit{D}}\mathrm{;}& & \end{array}$and $\begin{array}{ccc}{\mathit{\phi }}_{\mathrm{sky}}\mathrm{=}\frac{{\mathit{P}}_{\mathrm{FSL}}\mathrm{\otimes }{\mathit{T}}_{\mathrm{sky}}}{{\mathit{P}}^{\mathrm{pencil}}\mathrm{\otimes }{\mathit{T}}_{\mathrm{sky}}}\mathrm{·}& & \end{array}$Here φ_sky represents the relative pickup of anisotropy in the sidelobe beam as compared to the main beam and near sidelobes. The response of the FSL to CMB anisotropy is negligible, but Galactic response may not be. The PR spillover contributes only a constant value per pointing period, because it is not modulated with the scan. Only the SR direct spillover enters into the formula. Since the SR direct spillover is nearly aligned with the main beam, Galactic signal is not picked up in the CMB anisotropy except close to the Galactic plane. So for foreground-clean regions of the sky, φ_sky ≪ f_FSL.

The quantity φ_D is the bias due to the ratio of dipole response in the far sidelobes to dipole response in the main beam. Again the PR spillover contributes only an offset per pointing period, which is removed by the polkapix destriping algorithm (Tristram et al. 2011; Planck Collaboration VIII 2014), so we are left with the SR direct spillover and SR baffle spillover. Applying the same approximation as in Eq. (D.2) gives $\begin{array}{ccc}{\mathit{\phi }}_{\mathit{D}}\mathrm{\simeq }{\mathit{ϵ}}_{\mathrm{SRD}}\frac{\sin {\mathit{\theta }}_{\mathrm{SRD}}}{\sin {\mathit{\theta }}_{\mathrm{pencil}}}\mathit{,}& & \end{array}$and Eq. (D.5) becomes ${\mathit{T}}_{\mathrm{sky}}\mathrm{\simeq }\left(\mathrm{1}\mathrm{+}{\mathit{ϵ}}_{\mathrm{SRD}}\frac{\sin {\mathit{\theta }}_{\mathrm{SRD}}}{\sin {\mathit{\theta }}_{\mathrm{pencil}}}\right)\begin{array}{c}\mathrm{˜}\\ {\mathit{T}}_{\mathrm{sky}}\end{array}\mathit{.}$(D.6)This result implies that FSL bias of the effective beam window function tends to cancel the gain bias (Eq. (D.4)), and the response to CMB anisotropy is unaffected to first order.

Simulations of the sky scanned with the far sidelobe physical optics model (Planck Collaboration XIV 2014; Planck Collaboration VI 2014) confirm this.

Appendix E: Harmonic space computation of the effective beam window function

Two harmonic space techniques (FICSBell and Quickbeam), developed independently but closely related in their formalism, have been used to compute the effective beam window functions. They provide a valuable cross-check of the pixel-based results obtained with FEBeCoP (Fig. 2) and their low computational requirements allow fast calculation of the window function errors through the processing of Monte Carlo simulations of planet observations.

Appendix E.1: FICSBell

FICSBell is a harmonic space method for computing the effective beam window function directly from the scanning beam and the scan history. The steps of this method are as follows.

• 1.

  The statistics of the orientation of each detec-tor d within each map pixel p is computed first, andonly once for a given observation campaign:${\mathit{w}}_{\mathit{s}}^{\mathit{d}}\mathrm{\left(}{{r}}_{\mathit{p}}\mathrm{\right)}\mathrm{=}\sum _{\mathit{j}}{\mathrm{e}}^{\mathrm{i}\mathit{s}{\mathit{\psi }}_{\mathit{j}}}\mathit{,}$(E.1)where ψ_j is the orientation of the detector with respect to the local meridian during the measurement j occurring in the direction r_p. Note that the s = 0 moment is simply the hit count map. The orientation hit moments are computed up to degree s = 4. At the same time, the first two moments of the distribution of samples within each pixel (i.e., the centre of mass and moments of inertia) are computed and stored.

• 2.

  The scanning beam map of each detector d is transformed into spherical harmonics: ${\mathit{b}}_{\mathit{\ell s}}^{\mathit{d}}\mathrm{=}\mathrm{\int }\mathrm{d}{r}{\mathit{B}}_{\mathit{d}}\mathrm{\left(}{r}\mathrm{\right)}{\mathit{Y}}_{\mathit{\ell s}}\mathrm{\left(}{r}\mathrm{\right)}\mathit{,}$(E.2)where B_d(r) is the beam map centred on the north pole, and the Y_ℓs(r) are the spherical harmonic basis functions. Higher s indices describe higher degrees of departure from azimuthal symmetry and, for HFI beams, the coefficients ${\mathit{b}}_{\mathit{\ell s}}^{\mathit{d}}$ are decreasing functions of s at most multipoles considered. It also appears that, for ℓ< 3000, the coefficients with |s | > 4 account for much less than 1% of the beam solid angle. Spot checks where window functions are computed with |s| ≤ 6 show a difference of less than 10^-4 for ℓ< 2000 at 100 GHz and for ℓ< 3000 at 143 and 217 GHz. For these reasons, only modes with |s| ≤ 4 are considered in the present analysis. Armitage-Caplan & Wandelt (2009) reached a similar conclusion in their deconvolution of LFI beams.

• 3.

  For a given CMB sky realization t, described by its spherical harmonics coefficients a_ℓm = ^∫drt(r)Y_ℓm(r), the ${\mathit{b}}_{\mathit{\ell s}}^{\mathit{d}}$ coefficients computed above are used to generate s-spin weighted maps, ${\mathit{m}}_{\mathit{s}}^{\mathit{d}}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\sum _{\mathit{\ell m}}{\mathit{b}}_{\mathit{\ell s}}^{\mathit{d}}{\mathit{a}}_{\mathit{\ell m}}{\mathit{s}}_{}{\mathit{Y}}_{\mathit{\ell m}}\mathrm{\left(}{r}\mathrm{\right)}\mathit{,}$(E.3)as well as the first and second derivatives, using standard HEALPix tools.

• 4.

  The spin-weighted maps and orientation hit moments of the same order s are combined for all detectors involved, to provide an “observed” map $\mathit{m}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\left(\sum _{\mathit{d}}\sum _{\mathit{s}}{\mathit{w}}_{\mathit{s}}^{\mathit{d}}\mathrm{\left(}{r}\mathrm{\right)}{\mathit{m}}_{\mathit{s}}^{\mathit{d}}\mathrm{\left(}{r}\mathrm{\right)}\right)/\sum _{\mathit{d}}{\mathit{w}}_{\mathrm{0}}^{\mathit{d}}\mathrm{\left(}{r}\mathrm{\right)}\mathit{.}$(E.4)Similarly the local spatial derivatives are combined with the location hit moments to describe the effect of the non-ideal sampling of each pixel (see Appendix F). In this combination, the respective number of hits of each detector in each pixel is considered, as well as the detector weighting (generally proportional to the inverse noise variance).

• 5.

  The power spectrum of this map can then be computed, and compared to the input CMB power spectrum to estimate the effective beam window function over the whole sky, or over a given region of the sky.

Monte Carlo (MC) simulations in which the sky realizations are changed can be performed by repeating steps 3, 4 and 5. The impact of beam model uncertainties can be studied by including step 2 in the MC simulations.

Appendix E.2: Quickbeam

By decomposing the scanning beam into harmonic coefficients B_ℓm, each TOI sample can be modelled as (neglecting the contribution from instrumental noise, which is independent of beam asymmetry) ${\mathit{T}}_{\mathit{i}}\mathrm{=}\sum _{\mathit{\ell ms}}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{s}{\mathit{\alpha }}_{\mathit{i}}}{\mathit{B}}_{\mathit{\ell s}}\begin{array}{c}\mathrm{˜}\\ {\mathit{T}}_{\mathit{\ell m}}\end{array}{\mathit{s}}_{}{\mathit{Y}}_{\mathit{\ell m}}\mathrm{\left(}{\mathit{\theta }}_{\mathit{i}}\mathit{,}{\mathit{\phi }}_{\mathit{i}}\mathrm{\right)}\mathit{,}$(E.5)where the TOI samples are indexed by i, and $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{T}}_{\mathit{\ell m}}\end{array}$ is the underlying sky signal. The spin spherical harmonic ${\mathit{s}}_{}{\mathit{Y}}_{\mathit{\ell m}}$ rotates the scanning beam to the pointing location (θ,φ), while the e^{− isα_i} factor gives it the correct orientation. Equation (E.5)may be evaluated using techniques developed for convolution in Wandelt & Gorski (2001) and Prezeau & Reinecke (2010), although manipulating a TOI-sized object is necessarily slow.

On the small angular scales comparable to the size of the beam, it is a good approximation to assume that the procedure of mapmaking from TOI samples is essentially a process of binning: $\mathit{T}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{=}\sum _{\mathit{i}\mathrm{\in }\mathit{p}}{\mathit{T}}_{\mathit{i}}\mathit{/}\mathit{H}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathit{,}$(E.6)where H(p) is the total number of hits in pixel .

Given a normalized, rescaled harmonic transform of the beam B_ℓm, sky multipoles $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{T}}_{\mathit{\ell m}}\end{array}$ and a scan history object given by $\mathit{w}\mathrm{\left(}\mathit{n̂}\mathit{,s}\mathrm{\right)}\mathrm{=}\sum _{\mathit{j}\mathrm{\in }\mathit{p}}{\mathrm{e}}^{\mathrm{i}\mathit{s}{\mathit{\alpha }}_{\mathit{j}}}\mathit{/}\mathit{H}\mathrm{\left(}\mathit{n̂}\mathrm{\right)}\mathit{,}$(E.7)where the sum is over all hits j of pixel p at location , and α_j is the scan angle for observation j. The harmonic transform of this scan-strategy object is given by ${\mathit{s}}_{}{\mathit{w}}_{\mathit{LM}}\mathrm{=}\mathrm{\int }{\mathrm{d}}^{\mathrm{2}}\mathit{n̂}{\mathit{s}}_{}{\mathit{Y}}_{\mathit{LM}}^{\mathrm{\ast }}\mathrm{\left(}\mathit{n̂}\mathrm{\right)}\mathit{w}\mathrm{\left(}\mathit{n̂}\mathit{,s}\mathrm{\right)}\mathit{.}$(E.8)The beam-convolved observation is then given by $\begin{array}{c}\mathrm{˜}\\ \mathit{T}\end{array}\mathrm{\left(}\mathit{n̂}\mathrm{\right)}\mathrm{=}\sum _{\mathit{s\ell m}}\mathit{w}\mathrm{(}\mathit{n̂}\mathit{,}\mathrm{-}\mathit{s}\mathrm{)}{\mathit{B}}_{\mathit{\ell s}}{\mathit{T}}_{\mathit{\ell m}}{\mathit{s}}_{}{\mathit{Y}}_{\mathit{\ell m}}\mathrm{\left(}\mathit{n̂}\mathrm{\right)}\mathit{.}$(E.9)Taking the ensemble average of the pseudo-C_ℓ power spectrum of these T_lm gives (Hanson et al. 2010; Hinshaw et al. 2007) $\begin{array}{ccc}\stackrel{}{{{\mathit{C}}_{\mathit{L}}^{\mathit{TT}}}_{\mathrm{˜}}}& \mathrm{=}& \sum _{\mathit{S}{\mathit{S}}^{\mathrm{\prime }}}\sum _{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}}\frac{\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{1}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathrm{4}\mathit{\pi }}{}_{\mathrm{(}\mathrm{-}\mathit{s}\mathrm{-}{\mathit{s}}^{\mathrm{\prime }}\mathrm{)}}{\mathrm{𝒲}}_{{\mathit{\ell }}_{\mathrm{1}}}{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{2}}\mathit{S}}{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{S}}^{\mathrm{\prime }}}^{\mathrm{\ast }}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{2}}}^{\mathit{TT}}\\ & & \mathrm{\times }\left(\begin{array}{c}\\ & & \\ & & \end{array}\right)\left(\begin{array}{c}\\ & & \\ & & \end{array}\right)\end{array}$(E.10)where ${\mathrm{\left(}\mathit{s}{\mathit{s}}^{\mathrm{\prime }}\mathrm{\right)}}_{}{\mathrm{𝒲}}_{\mathit{L}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}}\sum _{\mathit{M}}{\mathit{S}}_{}{\mathit{w}}_{\mathit{LM}}{{\mathit{S}}^{\mathrm{\prime }}}_{}{\mathit{w}}_{\mathit{LM}}^{\mathrm{\ast }}$(E.11)is a cross-power spectrum of scan history objects. Note that the w(n,s) used here can also incorporate a position dependent weighting to optimize the pseudo-C_ℓ estimate, such as inverse-noise or a mask, the equations are unchanged. Writing the pseudo-C_ℓ in position space (following Dvorkin & Smith 2009) with Wigner-d matrices gives $\begin{array}{ccc}\stackrel{}{{{\mathit{C}}_{\mathit{L}}^{\mathit{TT}}}_{\mathrm{˜}}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{8}\mathit{\pi }}\sum _{\mathit{S}{\mathit{S}}^{\mathrm{\prime }}}{\mathrm{\int }}_{-1}^{\mathrm{1}}\mathrm{d}\mathit{z}{\mathit{d}}_{\mathrm{00}}^{\mathit{L}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\\ & & \mathrm{\times }\left[\sum _{{\mathit{\ell }}_{\mathrm{1}}}{\mathit{d}}_{\mathrm{-}\mathit{s}\mathrm{-}{\mathit{s}}^{\mathrm{\prime }}}^{{\mathit{\ell }}_{\mathrm{1}}}\mathrm{\left(}\mathit{z}\mathrm{\right)}{}_{\mathrm{(}\mathrm{-}\mathit{s}\mathrm{-}{\mathit{s}}^{\mathrm{\prime }}\mathrm{)}}{\mathrm{𝒲}}_{{\mathit{\ell }}_{\mathrm{1}}}\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{1}}\mathrm{+}\mathrm{1}\mathrm{)}\right]\\ & & \mathrm{\times }\left[\sum _{{\mathit{\ell }}_{\mathrm{2}}}{\mathit{d}}_{\mathit{s}{\mathit{s}}^{\mathrm{\prime }}}^{{\mathit{\ell }}_{\mathrm{2}}}\mathrm{\left(}\mathit{z}\mathrm{\right)}{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{2}}\mathit{S}}{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{S}}^{\mathrm{\prime }}}^{\mathrm{\ast }}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{2}}}^{\mathit{TT}}\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{)}\right]\mathit{.}\end{array}$(E.12)This integral can be implemented exactly using Gauss-Legendre quadrature, at a cost of $\mathrm{𝒪}\mathrm{\left(}{\mathit{\ell }}_{\max }^{\mathrm{2}}{\mathit{s}}_{\max }^{\mathrm{2}}\mathrm{\right)}$. For simplicity, the equations here are written for the auto-spectrum of a single detector, but the generalization to a map made by adding several detectors with different weightings is straightforward. The cost to compute all of the necessary terms exactly in that case becomes $\mathrm{𝒪}\mathrm{\left(}{\mathit{\ell }}_{\max }^{\mathrm{2}}{\mathit{s}}_{\max }^{\mathrm{2}}{\mathit{N}}_{\det }^{\mathrm{2}}\mathrm{\right)}$.

On the flat sky, beam convolution is multiplication in Fourier space by a beam rotated onto the scan direction. Multiple hits with different scan directions are incorporated by averaging (as the scan history objects above encapsulate).

A scan strategy which is fairly smooth across the sky is nearly equivalent to observing many independent flat-sky patches at high L. There is a fairly good approximation to the beam convolved pseudo-power spectrum which is essentially a flat-sky approximation. In the limit that L ≫ ℓ₁, with C_ℓ2 and B_ℓ2 is a slowly-varying function in ℓ₂, and using the equality $\sum _{{\mathit{\ell }}_{\mathrm{2}}}\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{)}\left(\begin{array}{c}\\ {\mathit{\ell }}_{\mathrm{1}}& {\mathit{\ell }}_{\mathrm{2}}& \mathit{L}\\ \mathit{s}& \mathrm{-}\mathit{s}& \mathrm{0}\end{array}\right)\left(\begin{array}{c}\\ {\mathit{\ell }}_{\mathrm{1}}& {\mathit{\ell }}_{\mathrm{2}}& \mathit{L}\\ {\mathit{s}}^{\mathrm{\prime }}& \mathrm{-}{\mathit{s}}^{\mathrm{\prime }}& \mathrm{0}\end{array}\right)\mathrm{=}{\mathit{\delta }}_{\mathit{s}{\mathit{s}}^{\mathrm{\prime }}}\mathit{,}$(E.13)the pseudo-C_ℓ sum above can be approximated as ${\begin{array}{c}\mathrm{˜}\\ \mathit{C}\end{array}}_{\mathit{L}}^{\mathit{TT}}\mathrm{=}{\mathit{C}}_{\mathit{L}}^{\mathit{TT}}\sum _{\mathit{M}}{⟨{\left|\mathit{w}\mathrm{\left(}\mathit{n̂}\mathit{p}\mathit{,M}\mathrm{\right)}\right|}^{\mathrm{2}}⟩}_{\mathit{p}}\mathrm{|}{\mathit{B}}_{\mathit{LM}}{\mathrm{|}}^{\mathrm{2}}\mathit{,}$(E.14)where the average ⟨⟩_p is taken over the full sky. It is illustrative to consider two limits of this equation. Firstly, for a “raster” scan strategy in which each pixel is observed with the same direction:${\mathrm{⟨}{\left|\mathit{w}\mathrm{\left(}\mathit{n̂}\mathit{,M}\mathrm{\right)}\right|}^{\mathrm{2}}\mathrm{⟩}}_{\mathit{p}}\mathrm{=}\mathrm{1}\mathit{,}$(E.15)and the predicted transfer function is just the power spectrum of the beam. Secondly, for a “best-case” scan strategy, in which each pixel is observed many times with many different orientation angles, ${\mathrm{⟨}{\left|\mathit{w}\mathrm{\left(}\mathit{n̂}\mathit{,M}\mathrm{\right)}\right|}^{\mathrm{2}}\mathrm{⟩}}_{\mathit{p}}\mathrm{=}{\mathit{\delta }}_{\mathit{M}\mathrm{0}}\mathit{,}$(E.16)and the transfer function is the azimuthally symmetric part of the beam. Note that this is for a full-sky observation; in the presence of a mask, the average above produces an f_sky factor, as expected but neglects the coupling between L multipoles (which can be calculated with the more complete equations above).

Appendix F: Pixelization artefacts

PlanckHFI maps are produced at HEALPix resolution 11 (N_side = 2048), corresponding to pixels with a typical dimension of 1.́7. With the resolution comparable to the spacing between scanning rings (Planck Collaboration I 2011) there is an uneven distribution of hits within pixels, introducing a complication in the analysis and interpretation of the Planck maps. A sample of the Planck distribution of sample hits within pixels is illustrated in Fig. F.

The three effective beam codes may also be used to simulate the effect of pixelization on the observed sky: LevelS/TotalConvoler/Conviqt (Reinecke et al. 2006; Wandelt & Gorski 2001; Prezeau & Reinecke 2010); FEBeCoP (Mitra et al. 2011; and FICSBell (Appendix E).

For the measurement of CMB fluctuations, the effects of pixelization may be studied analytically. On the small scales relevant to pixelization, the observed CMB is smooth, both due to physical damping and the convolution of the instrumental beam. Taylor expanding the CMB temperature about a pixel centre to second order, the typical gradient amplitude is given by $\mathrm{⟨}\mathrm{|}\mathrm{\nabla }\mathit{T}{\mathrm{|}}^{\mathrm{2}}\mathrm{⟩}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}\sum _{\mathit{\ell }}\mathit{\ell }\mathrm{(}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{\left)}\mathrm{\right(}\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{C}}_{\mathit{\ell }}^{\mathit{T}}{\mathit{W}}_{\mathit{\ell }}\mathrm{\approx }\mathrm{1}\mathrm{\times }{\mathrm{10}}^{\mathrm{9}}\mathit{\mu }{\mathrm{K}}^{\mathrm{2}}\mathit{/}{\mathrm{rad}}^{\mathrm{2}}\mathit{.}$(F.1)Here the approximate value is calculated for a ΛCDM cosmology with a 7′ FWHM Gaussian beam. The typical curvature of the observed temperature, on the other hand is given by $\mathrm{⟨}\mathrm{|}{\mathrm{\nabla }}^{\mathrm{2}}\mathit{T}{\mathrm{|}}^{\mathrm{2}}\mathrm{⟩}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}\sum _{\mathit{\ell }}\mathrm{[}\mathit{\ell }\mathrm{(}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}\mathrm{(}\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{C}}_{\mathit{\ell }}^{\mathit{T}}{\mathit{W}}_{\mathit{\ell }}\mathrm{\approx }\mathrm{7}\mathrm{\times }{\mathrm{10}}^{\mathrm{14}}\mathit{\mu }{\mathit{K}}^{\mathrm{2}}\mathit{/}{\mathrm{rad}}^{\mathrm{4}}\mathit{.}$(F.2)On the scales relevant to the maximum displacement from the centre of a 1.7′ pixel, the maximum displacement is of order 1′(3 × 10^-4rad)), and so the gradient term tends to dominate, although the curvature term is still non-negligible. For each observation of a pixel, the displacement from the pixel centre can be denoted as d = d_θ + id_φ. The average over all hits within a pixel gives an overall deflection vector for a pixel centre located at denoted as . This represents the centre of mass of the hit distribution; Fig. F shows these average deflections using black arrows. The deflection field may be decomposed into spin-1 spherical harmonics as ${\mathit{d}}_{\mathit{lm}}\mathrm{=}{\mathrm{\int }}_{\mathrm{4}\mathit{\pi }}{\mathrm{1}}_{}{\mathit{Y}}_{\mathit{\ell m}}^{\mathrm{\ast }}\mathrm{d}\mathrm{\left(}\mathit{n̂}\mathrm{\right)}\mathit{.}$(F.3)With a second-order Taylor expansion of the CMB temperature about each pixel centre, it is then possible to calculate the average pseudo-C_ℓ power spectrum of the pixelized sky. This is given by $\begin{array}{ccc}{\mathit{C}}_{\mathit{\ell }}^{\mathit{T}}\mathrm{=}\mathrm{[}\mathrm{1}\mathrm{-}\mathit{\ell }\mathrm{(}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{R}}^{\mathit{d}}\mathrm{]}{\mathit{C}}_{\mathit{\ell }}^{\mathit{T}}{\mathit{W}}_{\mathit{\ell }}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\sum _{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}}\frac{{\mathit{\ell }}_{\mathrm{1}}\mathrm{(}{\mathit{\ell }}_{\mathrm{1}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{1}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathrm{4}\mathit{\pi }}& & \\ \mathrm{\times }{\left(\begin{array}{c}\\ & & \\ & & \end{array}\right)}^{\mathrm{2}}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{1}}}^{\mathit{T}}{\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}}\mathrm{[}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{2}}}^{\mathit{d}\mathrm{+}}\mathrm{+}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{\ell }\mathrm{+}{\mathit{\ell }}_{\mathrm{1}}\mathrm{+}{\mathit{\ell }}_{\mathrm{2}}}{{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{2}}}^{\mathit{d}\mathrm{-}}}^{\mathrm{]}}\mathit{,}& & \end{array}$(F.4)where R^d = ⟨|d|² ⟩ /2 is half the mean-squared deflection magnitude (averaged over hits within a pixel, as well as over pixels), ${\mathit{C}}_{\mathit{\ell }}^{\mathit{d}\mathrm{+}}$ is the sum of the gradient and curl power spectra of d_ℓm, and ${\mathit{C}}_{\mathit{\ell }}^{\mathit{d}\mathrm{-}}$ is the gradient spectrum minus the curl spectrum. The R^d term describes a smearing of the observed sky due to pixelization. For uniform pixel coverage of N_side = 2048 pixels ⟨|d|² ⟩ ^1/2 = (2R^d)^1/2 = 0.725′, while, for the hit distribution of Planck frequency maps, R^d is within 0.2% of this value for CMB channels, and 0.4% for all channels. This term is therefore accurately described by the HEALPix pixel window function, which is derived under the assumption of uniform pixel coverage, and the resulting relative error on the beam window function is at most 4 × 10^-4 for ℓ ≤ 3000.

The effect of pixelization is degenerate with that of gravitational lensing of the CMB, with the difference that it: (1) acts on the beam-convolved sky, rather than the actual sky; and (2) produces a curl-mode deflection field as well as a gradient mode. This is discussed further in Planck Collaboration XVII (2014), where the sub-pixel deflection field constitutes a potential source of bias for the measured Planck lensing potential. Indeed, Eq. (F.4) is just a slightly modified version of the usual first order CMB lensing power spectrum (Hu 2000), Lewis & Challinor (2006) to accommodate curl modes.

A useful approximation to Eq. (F.4) which is derived in the unrealistic limit that the deflection vectors are uncorrelated between pixels, but in practice gives a good description of the power induced by the pixelization, is that the couples the CMB gradient into a source of noise with an effective level given by ${\mathit{\sigma }}^{\mathit{N}}\mathrm{\approx }\sqrt{{\mathit{R}}^{\mathrm{T}}\frac{\mathrm{4}\mathit{\pi }}{{\mathit{N}}_{\mathrm{pix}}}\mathrm{⟨}\mathrm{\left|}\mathit{d}\mathrm{\right(}\mathit{n̂}\mathrm{)}{\mathrm{|}}^{\mathrm{2}}\mathrm{⟩}}\mathit{,}$(F.5)where the average is taken over all pixels and R^T is half the mean-squared power in the CMB gradient: ${\mathit{R}}^{\mathit{T}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{8}\mathit{\pi }}\sum _{\mathit{\ell }}\mathit{\ell }\mathrm{(}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}\stackrel{}{{{\mathit{C}}_{\mathit{\ell }}^{\mathit{T}}}_{\mathrm{˜}}}\mathit{.}$(F.6)For frequency-combined maps, $\sqrt{\mathrm{⟨}\mathrm{\left|}\mathit{d}\mathrm{\right(}\mathit{n̂}\mathrm{)}{\mathrm{|}}^{\mathrm{2}}\mathrm{⟩}}$ is typically on the order of 0.́1, and so the induced noise σ^N is approximately 2 μK′. This is small compared to the instrumental contribution, although it does not disappear when taking cross-spectra, depending on the coherence of the hit distributions of the two maps in the cross-spectra.

Appendix G: Beam window function error

Appendix G.1: Error eigenmodes

Consider two individual detectors or detection units (weighted sum of detectors) X and Y for which sky maps are available. Here X and Y can be the same or different. Putting aside the instrumental noise and other contingencies for the time being, the cross- (or auto-) angular power spectrum measured of the observed maps ${\mathit{C}}_{\mathrm{obs}}^{\mathit{XY}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}$ is on average related to the real input signal ${\mathit{C}}_{\mathrm{sky}}^{\mathit{XY}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}$ through $\mathrm{⟨}{\mathit{C}}_{\mathrm{obs}}^{\mathit{XY}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{⟩}\mathrm{=}{\mathit{C}}_{\mathrm{sky}}^{\mathit{XY}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}{\mathit{W}}_{\mathrm{eff}\mathit{,}\mathrm{true}}^{\mathit{XY}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}$(G.1)where ${\mathit{W}}_{\mathrm{eff}\mathit{,}\mathrm{true}}^{\mathit{XY}}$ is the effective beam window function. Note that because of the optical beam non-circularity and Planck scanning strategy, ${\mathit{W}}^{\mathit{XY}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\ne {\mathrm{\left[}{\mathit{W}}^{\mathit{XX}}\mathrm{\right(}\mathit{\ell }\mathrm{\left)}{\mathit{W}}^{\mathit{YY}}\mathrm{\right(}\mathit{\ell }{\mathrm{\left)}}^{\mathrm{\right]}}}^{\mathrm{1}\mathit{/}\mathrm{2}} \mathrm{if}\mathit{X}\ne \mathit{Y,}$(G.2)as illustrated in Fig. 16, while W^XY = W^YX for any X and Y. It also appears that in the ℓ range of interest, W^XY(ℓ) ≥ 0; therefore ${\mathit{W}}^{\mathit{XY}}\mathrm{=}{{}^{\mathrm{(}}\mathit{B}^{\mathit{XY}}{}^{\mathrm{)}}}^{\mathrm{2}}$, analogous to what is usually done for simple (circular) beam models. In what follows, the XY pair superscript is dropped except when they are required for clarity. In most cases, scientific analyses will be conducted on a best guess C_est(ℓ) of the sky power spectrum, in which the empirical C_obs(ℓ) is corrected from a nominal effective window W_eff,nom(ℓ)${\mathit{C}}_{\mathrm{est}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{=}{\mathit{C}}_{\mathrm{obs}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathit{/}{\mathit{W}}_{\mathrm{eff}\mathit{,}\mathrm{nom}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{;}$(G.3)therefore, on average, $\begin{array}{ccc}\mathrm{⟨}{\mathit{C}}_{\mathrm{est}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{⟩}& \mathrm{=}& {\mathit{C}}_{\mathrm{sky}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}{\mathit{W}}_{\mathrm{eff}\mathit{,}\mathrm{true}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathit{/}{\mathit{W}}_{\mathrm{eff}\mathit{,}\mathrm{nom}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathit{,}\\ & \mathrm{=}& {\mathit{C}}_{\mathrm{sky}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}{\left({\mathit{B}}_{\mathrm{eff}\mathit{,}\mathrm{true}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathit{/}{\mathit{B}}_{\mathrm{eff}\mathit{,}\mathrm{nom}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\right)}^{\mathrm{2}}\mathit{.}\end{array}$(G.4)The ratio B_eff,true(ℓ) /B_eff,nom(ℓ) which determines the uncertainties on the angular power spectrum originating from our beam knowledge is studied using the planet transit MC simulations described in Sect. 5. The scanning beam map determined with the B-Spline code described in Appendix C on each of these simulations is turned into an effective beam window function W_i(ℓ) for i = 1,2,...,n_MC (where n_MC = 100 is the number of MC simulations) using the Quickbeam formalism described in Appendix E.2.

Defining the means B_mean(ℓ) and W_mean(ℓ) as ${\mathit{B}}_{\mathrm{mean}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathrm{MC}}}\mathrm{\left(}{\mathit{W}}_{\mathit{i}}\mathrm{\right(}\mathit{\ell }\mathrm{)}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{/}{\mathit{n}}_{\mathrm{MC}}\mathit{,}$(G.5)and, ${\mathit{W}}_{\mathrm{mean}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathrm{MC}}}{\mathit{W}}_{\mathit{i}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathit{/}{\mathit{n}}_{\mathrm{MC}}\mathit{,}$(G.6)one can build the matrix of the deviations around the mean ${\mathrm{\Delta }}_{\mathit{i}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{=}{\mathit{f}}_{\mathit{s}}\ln \left({\mathit{B}}_{\mathit{i}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathit{/}{\mathit{B}}_{\mathrm{mean}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\right)\mathit{,}$(G.7)where the factor f_s has been applied. As discussed in Sect. 6, f_s = 2.7. This scaling factor is applied throughout the rest of discussion and is included in the delivered products and plotted error modes.

Sect. G.3 contains a discussion of how the MC average B_mean and nominal beam B_eff,nom are related and focus here on the dispersion around the mean.

Since the relative dispersion of the simulated W_i(ℓ) generally is very small (less than 1%), then W_mean(ℓ) ≃ B_mean(ℓ)² to a very good approximation (as illustrated in Fig. G.1) and the matrix Δ is well approximated by $\begin{array}{ccc}{\mathrm{\Delta }}_{\mathit{i}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}& \mathrm{\simeq }& \mathrm{1}\mathit{/}\mathrm{2}{\mathit{f}}_{\mathit{s}}\ln \left[{\mathit{W}}_{\mathit{i}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathit{/}{\mathit{W}}_{\mathrm{mean}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\right]\mathit{,}\\ & \mathrm{\simeq }& \mathrm{1}\mathit{/}\mathrm{2}{\mathit{f}}_{\mathit{s}}\left[{\mathit{W}}_{\mathit{i}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathit{/}{\mathit{W}}_{\mathrm{mean}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{-}\mathrm{1}\right]\mathit{.}\end{array}$(G.8)The quantity B(ℓ) was preferred over W(ℓ) in order to remain consistent with the usual description of the beam in linear map space, instead of quadratic space.

The statistical properties of the MC based beam window functions can be summarized in the covariance of the deviations Δ, defined as ${C}\mathrm{\equiv }{{\Delta }}^{\mathrm{T}}\mathrm{·}{\Delta }\mathit{/}\mathrm{(}{\mathit{n}}_{\mathrm{MC}}\mathrm{-}\mathrm{1}\mathrm{)}\mathit{,}$(G.9)where Δ is a matrix with n_MC rows and ℓ_max + 1 columns, and C is a square symmetric positive semi-definite matrix with ℓ_max + 1 rows and columns. It can be diagonalized into $\begin{array}{ccc}{C}& \mathrm{=}& {V}\mathrm{·}{{D}}^{\mathrm{2}}\mathrm{·}{{V}}^{\mathrm{T}}\mathit{/}\mathrm{(}{\mathit{n}}_{\mathrm{MC}}\mathrm{-}\mathrm{1}\mathrm{)}\mathit{,}\\ & \mathrm{=}& \end{array}$(G.10)where D is a diagonal matrix, with at most n_MC positive eigenvalues, and the eigenmodes matrix ${E}\mathrm{\equiv }{D}\mathrm{·}{{V}}^{\mathrm{T}}\mathit{/}\mathrm{(}{\mathit{n}}_{\mathrm{MC}}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$(G.11)of the beam uncertainty. Alternatively, one can perform a singular value decomposition (SVD) of Δ, which reads $\begin{array}{ccc}{\Delta }& \mathrm{=}& \\ & \mathrm{=}& {M}\mathrm{·}{E}\mathrm{(}{\mathit{n}}_{\mathrm{MC}}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}\end{array}$where M is an orthogonal n_MC × n_MC matrix (i.e., M^T·M = M·M^T = I_nMC), D is a diagonal matrix with n_MC non-negative eigenvalues of decreasing amplitude, and V is a matrix with ℓ_max + 1 rows whose n_MC columns are orthonormal vectors (i.e., V^TV = I_nMC), with I_nMC being the identity matrix of size n_MC × n_MC. The diagonalization of the covariance matrix C (Eq. (G.10)) has a numerical complexity scaling like ${\mathit{\ell }}_{\max }^{\mathrm{3}}$, while the SVD of Δ (Eq. (G.12)) scales like ${\mathit{\ell }}_{\max }{\mathit{n}}_{\mathrm{MC}}^{\mathrm{2}}$. Since n_MC ≪ ℓ_max the latter approach was prefered because it is much faster, and it naturally provides the mixing matrix M. Equation (G.12) has a degeneracy on the sign of M and V, which was lifted by constraining the vectors of V (and therefore E) to all be positive at ℓ = 200.

It appears that most of the statistical content of Δ or C is described by the first few modes n_modes with the largest eigenvalues, in which case the E matrix is truncated to its first n_modes rows with little loss of information. For HFI, n_modes = 5 is chosen, as illustrated in Fig. G.1.

The statistics of the elements of the mixing matrix M, and therefore of the MC measured beam window function fluctuations, is shown in Fig. G.3 to be very close to Gaussian, justifying the current analysis in terms of a covariance matrix.

The beam uncertainty model therefore is $\begin{array}{ccc}\mathit{B}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}& \mathrm{=}& {\mathit{B}}_{\mathrm{mean}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\exp {\mathrm{(}{{g}}^{\mathrm{T}}\mathit{.}{{E}}^{\mathrm{)}}}_{\mathit{\ell }}\\ & \mathrm{=}& {\mathit{B}}_{\mathrm{mean}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\exp \left(\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathrm{modes}}}{\mathit{g}}_{\mathit{k}}{\mathit{e}}_{\mathit{k}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\right)\mathit{,}\end{array}$where g is a n_modes element vector of independant Gaussian variates of zero mean and unit variance and e_k(ℓ) is the kth row of E.

The SVD decomposition of the beam uncertainty was performed for the ℓ range [ ℓ_min,ℓ_max ] with ℓ_min = 0 and ℓ_max depending on the frequencies of the two detectors involved in the beam considered. Currently ℓ_max = 2500 when the lowest frequency is 100 GHz, ℓ_max = 3000 when that frequency is 143 GHz, and ℓ_max = 4000 at higher frequency.

If orthogonal error eigenmodes are needed for the range $\mathrm{\left[}{\mathit{\ell }}_{\min }^{\mathrm{\prime }}\mathit{,}{\mathit{\ell }}_{\max }^{\mathrm{\prime }}\mathrm{\right]}$, with ${\mathit{\ell }}_{\min }^{\mathrm{\prime }}\mathrm{\ge }{\mathit{\ell }}_{\min }$ and ${\mathit{\ell }}_{\max }^{\mathrm{\prime }}\mathrm{\le }{\mathit{\ell }}_{\max }$, the provided E must first be truncated to the new range to give the E_t matrix with n_modes rows and ${\mathit{\ell }}_{\max }^{\mathrm{\prime }}\mathrm{-}{\mathit{\ell }}_{\min }^{\mathrm{\prime }}\mathrm{+}\mathrm{1}$ columns, and then singular value decomposed into ${{E}}_{\mathit{t}}\mathrm{=}{{R}}^{\mathrm{\prime }}\mathrm{·}{{D}}^{\mathrm{\prime }}\mathrm{·}{{V}}^{\mathrm{\prime }\mathrm{T}}\mathit{,}$(G.16)where the new set of orthogonal eigenmodes is ${{E}}^{\mathrm{\prime }}\mathrm{\equiv }{{D}}^{\mathrm{\prime }}\mathrm{·}{{V}}^{\mathrm{\prime }\mathrm{T}}\mathit{.}$(G.17)This is illustrated in Fig. G.2, where the eigenmodes are truncated to the frequency dependent ℓ ranges used in the PlanckC(ℓ) likelihood (Planck Collaboration XV 2014).

Appendix G.2: Eigenmode covariance

The previous approach, dedicated to the study of the uncertainty on the beam window function associated with a single pair of detectors (or maps), can be generalized to the simultaneous characterization of any number of pairs. For instance, for the three disjoint pairs, a = { UV }, b = { XY } and c = { ZT }, one writes $\mathrm{(}{{\Delta }}^{\mathit{a}}{{\Delta }}^{\mathit{b}}{{\Delta }}^{\mathit{c}}\mathrm{)}\mathrm{=}\mathrm{(}{{M}}^{\mathit{a}}{{M}}^{\mathit{b}}{{M}}^{\mathit{c}}\mathrm{)}\mathit{.}\left(\begin{array}{c}\\ {{E}}^{\mathit{a}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {{E}}^{\mathit{b}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {{E}}^{\mathit{c}}\end{array}\right)\mathit{.}$(G.18)and the covariance matrix reads $\begin{array}{ccc}{{C}}^{\mathrm{abc}}& \mathrm{=}& {\mathrm{(}{{\Delta }}^{\mathit{a}}{{\Delta }}^{\mathit{b}}{{\Delta }}^{\mathit{c}}\mathrm{)}}^{\mathrm{T}}\mathit{.}\mathrm{(}{{\Delta }}^{\mathit{a}}{{\Delta }}^{\mathit{b}}{{\Delta }}^{\mathit{c}}\mathrm{)}\\ & \mathrm{=}& {\left(\begin{array}{c}\\ & & \\ & & \\ & & \end{array}\right)}^{\mathrm{T}}\mathit{.}\left(\begin{array}{c}\\ & & \\ & & \\ & & \end{array}\right)\mathit{.}\left(\begin{array}{c}\\ & & \\ & & \\ & & \end{array}\right)\mathit{,}\\ & \mathrm{=}& {{E}}^{\mathit{a,b,cT}}\mathrm{·}{{M}}^{\mathit{a,b,c}}\mathrm{·}{{E}}^{\mathit{a,b,c}}\mathit{,}\end{array}$(G.19)where M^a,b,c is a square symmetric matrix with 3n_modes rows, and if one denotes L^a,b,c its Cholesky “root,” such that M^a,b,c = L^a,b,c·L^a,b,cT, then $\begin{array}{ccc}{\mathit{B}}^{\mathit{a}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}& \mathrm{=}& \\ {\mathit{B}}^{\mathit{b}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}& \mathrm{=}& \\ {\mathit{B}}^{\mathit{c}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}& \mathrm{=}& \end{array}$where g is the 3n_modes element vector of independent Gaussian deviates and is the same for Eqs. (G.20) to (G.22).

The Planck-HFI Reduced Instrument Model (RIMO) available at Planck Legacy Archive⁸ and described in Planck Collaboration (2013) contains the correlation matrix M^a,b,c,d,..., where a,b,c,d... each are a different element of the set of pairs that can be built out of the detection units available. So, for n_d detection units, the number of pairs will be n_d(n_d + 1)/2 and the correlation matrix will be square and symmetric, with the value “1” on its diagonal and have n_modesn_d(n_d + 1)/2 rows and as many columns. The n_modes relative to the same detector pair form adjacent rows and columns in the matrix, and the order of appearance of the pairs in the matrix is specified in the header of the FITS extension containing the matrix. The nominal B_ℓ and eigenmodes E are provided for each pair in a different extension.

The results above were obtained in the basis of eigenmodes provided, if one wants to obtain a beam correlated model on a different ℓ-range, with the eigenmodes E′ defined in Eq. (G.17), then the covariance matrix becomes ${{{C}}^{\mathrm{\prime }}}^{\mathrm{abc}}\mathrm{=}{\left(\begin{array}{c}\\ {{{E}}^{\mathrm{\prime }}}^{\mathit{a}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {{{E}}^{\mathrm{\prime }}}^{\mathit{b}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {{{E}}^{\mathrm{\prime }}}^{\mathit{c}}\end{array}\right)}^{\mathrm{T}}\mathrm{·}{{M}}^{\mathrm{\prime }\mathit{a,b,c}}\mathit{.}\left(\begin{array}{c}\\ {{{E}}^{\mathrm{\prime }}}^{\mathit{a}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {{{E}}^{\mathrm{\prime }}}^{\mathit{b}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {{{E}}^{\mathrm{\prime }}}^{\mathit{c}}\end{array}\right)\mathit{,}$(G.23)with ${{{M}}^{\mathrm{\prime }}}^{\mathit{a,b,c}}\mathrm{=}{\left(\begin{array}{c}\\ {{{R}}^{\mathrm{\prime }}}^{\mathit{a}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {{{R}}^{\mathrm{\prime }}}^{\mathit{b}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {{{R}}^{\mathrm{\prime }}}^{\mathit{c}}\end{array}\right)}^{\mathrm{T}}\mathrm{·}{{M}}^{\mathit{a,b,c}}\mathit{.}\left(\begin{array}{c}\\ {{{R}}^{\mathrm{\prime }}}^{\mathit{a}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {{{R}}^{\mathrm{\prime }}}^{\mathit{b}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {{{R}}^{\mathrm{\prime }}}^{\mathit{c}}\end{array}\right)\mathit{,}$(G.24)where the R′ matrices are obtained from the SVD in Eq. (G.16).

Appendix G.3: Monte Carlo bias

As discussed in Sect. 3.2, the Monte Carlo average of the B-Spline reconstructed effective beam window function B_mean(ℓ) can be different from the effective beam that would be obtained directly out of the simulation input beam map B_eff,in(ℓ), introduces a bias $\mathrm{1}\mathrm{+}{\mathit{\epsilon }}_{\mathrm{bias}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{\equiv }{\mathit{B}}_{\mathrm{mean}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\hspace{0.17em}\mathit{/}\hspace{0.17em}{\mathit{B}}_{\mathrm{eff}\mathit{,}\mathrm{in}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathit{,}$(G.25)which is interpreted as a consequence of the imperfect sampling of the beam map by the planet, the effect of the instrumental noise and pointing uncertainty and the smoothing feature of the B-spline functions. It is found that |ε_bias(ℓ = 500)| ≤ 2 × 10^-4 and |ε_bias(ℓ = 1000)| < 5 × 10^-4 for all effective beams, making it smaller than the relative dispersion on the beam window function described above. It is assumed that the beam reconstruction on the real data suffers from the same bias, and its beam window function B_eff,raw(ℓ) is corrected for this to provide the nominal beam $\begin{array}{ccc}{\mathit{B}}_{\mathrm{eff}\mathit{,}\mathrm{nom}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}& \mathrm{=}& {\mathit{B}}_{\mathrm{eff}\mathit{,}\mathrm{raw}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\hspace{0.17em}\mathit{/}\hspace{0.17em}\left[\mathrm{1}\mathrm{+}{\mathit{\epsilon }}_{\mathrm{bias}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\right]\mathit{,}\\ & \mathrm{=}& {\mathit{B}}_{\mathrm{eff}\mathit{,}\mathrm{raw}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\hspace{0.17em}{\mathit{B}}_{\mathrm{eff}\mathit{,}\mathrm{in}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\hspace{0.17em}\mathit{/}\hspace{0.17em}{\mathit{B}}_{\mathrm{mean}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathit{.}\end{array}$In doing so, ε_bias is assumed to be perfectly determined by the simulations, and no error is associated with this correction.

All Tables

All Figures