Difference between revisions of "Beams"

From Planck Legacy Archive Wiki
Jump to: navigation, search
m
 
(239 intermediate revisions by 9 users not shown)
Line 1: Line 1:
== Scanning Beams ==
 
  
The scanning beams describe the instrument’s instantaneous beam profile. Due to the near constant spin rate of the spacecraft, time domain effects (including residual time response and lowpass filtering) are degenerate with the spatial response due to the optical system. The scanning beam reconstruction recovers both of these effects, aside from residual time domain effects on a longer time scale than can be captured with the extent of the scanning beam model.
+
= Introduction =
  
In the paper <font color=red> cite P03c</font> we consider two models of the beam in order to better understand systematics in the reconstruction. Here we describe only the B-Spline beams which are used to compute the delivered effective beam (see next section).
+
The reconstruction of the spatial response for the Planck High Frequency Instrument (HFI) is informed by observation of planets. Furthermore, observations of Uranus and Neptune are used to obtain an absolute calibration for the 545 and 857 GHz channels.  
  
=== B-Spline Beam construction ===
+
The spatial response, also known as the beam response or the point spread function, relates the absolute calibration to different angular scales.
 +
We follow the nomenclature of Planck Collaboration VII {{PlanckPapers|planck2013-p03c}} (2014), where the "scanning beam" is defined as the coupled response of the optical system, the deconvolved time response function, and the software low-pass filter, applied to the data. The "effective beam" represents the averaging of the signal due to the scanning of the telescope and the mapmaking procedure, and varies from pixel to pixel across the sky.
  
We use seasons 1 and 2 of the Mars observation to reconstruct the beam. The data are processed with the bigPlanets TOI processing. We use JPL Horizons ephemerides to determine the pointing of each detector relative to the planet. We subtract the astrophysical background in the time domain using a bicubic interpolation of the Planck maps.
+
= The 2015 and 2018 Data Releases =
  
The time ordered data are used to fit a two dimensional B-Spline surface using a least square minimization and a smoothing criterion to minimize the effects of high spatial frequency variations. We therefore assume the scanning beam to be smooth. The smoothing criterion as well as the locations of the nodes used to compute the B-Spline basis functions are set using GRASP physical optics simulations as inputs which are the best assumptions on the spatial frequency content of the in-flight beams.
+
Most of the content found on this page, as well as further information regarding the beam products used for the 2015 (PR2) release of ''Planck'' data, can be found in the HFI DPC Paper 1 {{PlanckPapers|planck2014-a08}}.
 +
 
 +
For the PR3 2018 data release, the scanning beams are the same as for the previous release. The effective beams (Febecop) are very slightly affected by the rejection of the data taken during the last 1000 stable pointing periods of the cold mission (see {{PlanckPapers|planck2016-l03}}). This affects a small fraction of the sky in a negligible way.
 +
 
 +
The main improvement in PR3 over previous releases concerns the computation of effective beam window functions relevant for cosmological analysis. As explained in [[Beam_Window_Functions]], beam matrices describing exactly the leakage from temperature to polarization and cross-talk between
 +
spectra induced by the beam mismatch and beam ellipticity have been computed with <code>QuickPol</code> for all relevant combination of detectors and data splits. Accounting for such systematic effects in this way greatly improves the inter-frequency consistency of the CMB TE power spectra, compared to the 2015 analysis (see {{PlanckPapers|planck2016-l05}}).
 +
 
 +
 
 +
== Planet and main beam description ==
 +
 
 +
Here we redefine the "main beam" to be the scanning beam out to 100 arcmin from the beam axis. The sidelobe structure at this radius is dominated by diffraction at the mirror edges and falls as &theta;<sup>-3</sup>, where &theta; is the angle from the main beam axis. The main beam is used to compute the effective beam and the effective beam window function, which describes the filtering of sky signals. The smearing of the main beam cannot be significantly reduced without boosting the high-frequency noise. The regularization function (a low-pass filter) chosen has approximately the same width as the instrumental transfer function. The deconvolution significantly reduces the long tail of the scanning beam (see Fig. 1). The deconvolution also produces a more symmetric time response, so residual “streaking” appears both ahead and behind the main beam, although ahead of the beam it is at a level of less than 10<sup>−4</sup> of the peak response. The far sidelobes are defined as the response for &theta; > 5&deg;, roughly the minimum in the optical response. The response begins to rise as a function of angle &theta; beyond this due to spillover. The far sidelobes are handled separately from the beam effects (see Paper 2 {{PlanckPapers|planck2014-a09}} for details).
 +
 
 +
[[Image:Jupiter_Deconv_143-1b.png|thumb|500px|center|'''Figure 1:''' Scan across Jupiter with bolometer 143-1b, with pre-deconvolution in black and post-deconvolution in red. The plotted <i>y</i>-axis
 +
is linear within the range &plusmn; 10<sup>-6</sup> and logarithmic elsewhere. Similarly, the <i>x</i>-axis range is linear within the range &plusmn; 0.16s, and logarithmic outside.]]
 +
 
 +
The observations of planets are also used to estimate their brightness as a function of frequency, see {{PlanckPapers|planck2017-LII}}.
 +
 
 +
== Key Changes between 2013 and 2015 ==
 +
 
 +
There are several key changes in the reconstruction of the main beam since the 2013 data release, including those listed below.
 +
 
 +
* The time-ordered data from Saturn and Jupiter observations are merged prior to B-spline decomposition, accounting for residual pointing errors and variable seasonal brightness.  This is achieved by determining a scaling factor and a pointing offset by fitting the time-ordered data to a template from a previous estimate of the scanning beams. We iterate the planet data treatment updating the template with the reconstructed scanning beam. The process converges in five iterations within 0.1% accuracy in the effective beam window function.
 +
* Steep gradients in the signal close to the planet reduce the completeness of the standard glitch-detection and subtraction procedure; therefore the planet timelines are deglitched a second time.
 +
* The beam pipeline destripes the planet data, estimating a single baseline between 3&deg; and 5&deg; before the peak for each scanning circle. Baseline values are smoothed with a 40 circle sliding window. The entire scanning circle is removed from the beam reconstruction if the statistic in the timeline region used to estimate the baseline is far from Gaussian.
 +
* The main beam is now recontructed on a square grid that extends to a radius of 100 arcmin from the centroid, as opposed to 40 arcmin in the 2013 data release.  The cutoff of 100 arcmin was chosen so that a diffraction model of the beam at large angles from the centroid predicts that less than 5&times;10<sup>-5</sup> of the total solid angle is missing.
 +
* The scanning beam is constructed by combining data from Saturn observations, Jupiter observations, and physical optics models using GRASP software.
 +
* No apodization is applied to the scanning beam map.
 +
 
 +
The update of the time response deconvolved data has slightly changed the scanning beam, the effective beam solid angles, and the effective beam window functions.
 +
 
 +
{| class="wikitable" border="1" cellpadding="3" cellspacing="0" align="center" style="float:right; margin-left: 10px;
 +
|+ <small>'''Table 1. Band-average scanning beam solid angle (<math>\Omega _\mathrm{SB} </math>) and Monte Carlo-derived errors (<math>\Delta\Omega _\mathrm{MC}</math>) including noise, residual glitches, and pointing uncertainty.'''</small>
 +
|-
 +
|'''Band'''
 +
|<math>\Omega _\mathrm{SB}</math>
 +
|<math>\Delta \Omega _\mathrm{MC}</math>
 +
|-
 +
|[GHz]
 +
|[arcmin<math>^2</math>]
 +
|
 +
|-
 +
|100
 +
|104.62
 +
|0.13 %
 +
|-
 +
|143
 +
|58.50
 +
|0.07 %
 +
|-
 +
|217
 +
|26.92
 +
|0.13 %
 +
|-
 +
|353
 +
|25.93
 +
|0.09 %
 +
|-
 +
|545
 +
|25.23
 +
|0.08 %
 +
|-
 +
|857
 +
|23.04
 +
|0.08 %
 +
|}
 +
 
 +
A portion of the near sidelobes was not accounted for in the effective beam window function of the 2013 data (see {{PlanckPapers|planck2013-p03}}). To remedy this, the domain of the main beam reconstruction is extended to 100 arcmin, with no apodization. Saturn data are used where they are signal-dominated. Where the signal-to-noise ratio of the Saturn data fall below 9, azimuthally-binned Jupiter data are used.  At larger angles, below the noise floor of the Jupiter data, we use a power law (&prop; &theta;<sup>-3</sup>), whose exponent is derived from GRASP to extend the beam model to 100 arcmin.  Figure [[#HybridBeamDiagram|Fig. 2]] shows a scanning beam map with a diagram of the regions handled differently in the hybrid beam model. A summary of the solid angles of the hybrid beams is shown in Table 1. [[#FocalPlane_Map_DX11d|Figure 3]] shows a contour plot of all the scanning beams referenced to the centre of the focal plane.
 +
 
 +
[[Image:HybridBeamDiagram.png|thumb|500px|center|'''Figure 2:''' Scanning beam map for detector 143-6, indicating approximately the regions that are handled with a different data selection or binning.]]
 +
 
 +
[[Image:FocalPlane_Map_DX11d.png|thumb|800px|center|'''Figure 3:''' B-spline hybrid scanning beams reconstructed from Mars, Saturn, and Jupiter. The beams are plotted in logarithmic contours of −3, −10, −20, and −30 dB from the peak. PSB pairs are indicated with the "a" bolometer in black and the "b" bolometer in blue.]]
 +
 
 +
== Beam Error Budget ==
 +
 
 +
As in the 2013 release, the beam error budget is based on an eigenmode decomposition of the scatter in simulated planet observations. A reconstruction bias is estimated from the ensemble average of the simulations. We generate 100 simulations for each planet observation, including pointing uncertainty, cosmic ray glitches, and the measured noise spectrum. Simulated glitches are injected into the timeline with the correct energy spectrum and rate ({{PlanckPapers|planck2013-p03e}}) and are detected and removed using the deglitch algorithm. Noise realizations are derived as described in
 +
{{PlanckPapers|planck2013-p03c}}.
 +
Pointing uncertainty is simulated by randomizing the pointing by 1.5 arcsec rms in each direction. The improved signal-to-noise ratio compared to 2013 leads to smaller error bars: for instance, at &#8467; = 1000 the uncertainties on <i>B</i><sub>&#8467;</sub> are now (2.2, 0.84, 0.81)&times;10<sup>-4</sup> for 100, 143, and 217 GHz frequency-averaged maps respectively, reduced from the previous uncertainties of (61, 23, 20)&times;10<sup>-4</sup>.
 +
 
 +
A “Reduced Instrument Model” (RIMO; see [[The_RIMO]]) containing the effective <i>B</i><sub>&#8467;</sub> for temperature and polarization detector sets (for auto- and cross-spectra at 100 to 217 GHz) is included in the release, for a sky fraction of 100 %; another RIMO is provided for a sky fraction of 75 %. They both contain the first five beam error eigenmodes and their covariance matrix, for the multipole ranges [0, &#8467;<sub>max</sub>] with &#8467;<sub>max</sub> = 2000, 3000, and 3000 at 100, 143, and 217 GHz, respectively (instead of 2500, 3000, and 4000 previously). These new ranges bracket more closely the ones expected to be used in the likelihood analyses, and ensure a better determination of the leading modes on the customized ranges. As described in appendix A7 of Planck Collaboration XV (2014) {{PlanckPapers|planck2013-p08}}, these beam window function uncertainty eigenmodes are used to build the <i>C</i><sub>&#8467;</sub> covariance matrix used in the high-&#8467; angular power spectrum likelihood analysis. It was found that the beam errors are negligible compared to the other sources of uncertainty and have no noticeable impact on the values or associated errors of the cosmological parameters.
 +
 
 +
== Consistency of Beam Reconstruction ==
 +
 
 +
To evaluate the accuracy and consistency of the beam reconstruction method, we have compared beams reconstructed using Mars for the main beam part instead of Saturn, and using data from Year 1 or Year 2 only. We compared the window functions obtained with these new beams to the reference ones; new Monte Carlo simulations were created to evaluate the corresponding error bars, and the results are shown in Figs. [[#Consistency_Yearly_Beams_DX11d|4]] and [[#Consistency_Mars_Beams_DX11d|5]]. Note that the effective beam window functions shown in these figures are not exactly weighted for the scan strategy, rather we plot them in the raster scan limit, where <i>B</i><sub>&#8467;</sub> is defined as the sum over <i>m</i> of the <i>B</i><sub>&#8467;</sub><sup>m</sup>. The eigenmodes were evaluated for these different data sets, allowing us to compare them with the reference beam through a &chi;<sup>2</sup> analysis. The discrepancy between each data set and the reference beam is fitted using <i>N</i><sub>dof</sub>=5 eigenmodes using &#8467;<sub>max</sub> = 1200, 2000, and 2500 for 100, 143, and 217GHz channels, respectively. The &chi;<sup>2</sup> is then defined as
 +
 
 +
<math>
 +
\chi^2 = \sum_{i=1}^{N_\mathrm{dof}} (c_i / \lambda_i)^2,
 +
\label{eq:chi2}
 +
</math>
 +
 
 +
where <i>c<sub>i</sub></i> is the fit coefficient for eigenvector <i>i</i>, with eigenvalue &lambda;<sub>i</sub>. For each data set, the <i>p</i>-value to exceed this &chi;<sup>2</sup> value for a &chi;<sup>2</sup> distribution with <i>N</i><sub>dof</sub> degrees of freedom is indicated in Table 2. A high <i>p</i>-value indicates that the given data set is consistent with the reference beam within the simulation-determined error eigenvectors, with 100 % indicating perfect agreement. We find excellent agreement between the yearly and nominal beams for the SWB bolometer channels; however, in order to obtain reasonable beam agreement for the PSB detector sets, we find that the yearly beam errors must be scaled by a factor of 4.  We conclude that there is an unknown systematic error in the beam reconstruction that is not accounted in the simulations used to estimate the beam accounting in the Monte Carlo error bars. We assign a scaling factor of 4 to the error eigenmodes to account for this systematic error. This additional error scaling has a negligible effect on the cosmological parameters.
 +
 
 +
[[Image:Consistency_Yearly_Beams_DX11d.png|thumb|800px|center|'''Figure 4:''' Comparison of Year 1 and Year 2 based beams with the reference beam (Full Mission). Window functions are calculated using Quickbeam in raster scan configuration. Error bars are computed using MC simulations for Year 1 and Year 2, taking the maximal one for each.]]
 +
 
 +
[[Image:Consistency_Mars_Beams_DX11d.png|thumb|600px|center|'''Figure 5:''' Comparison of Mars-based beams with the reference beam (Full Mission). Mars-based beams are constructed from Mars data for the main beam (below about 10 arcmin) and Saturn and Jupiter data for the larger scales. Window functions are calculated using Quickbeam in raster scan configuration. Error bars are computed using MC simulations for Mars.]]
 +
 
 +
{| class="wikitable" border="1" cellpadding="3" cellspacing="0" align="center"
 +
|+ <small>'''Table 2. <i>P</i>-values in percent for the &chi;<sup>2</sup> comparison of the nominal beam with Year 1 and Year 2 data sets, following the definition in the above equation.'''</small>
 +
|-
 +
|'''Band'''
 +
|'''Year 1'''
 +
|'''Year 2'''
 +
|-
 +
|100-ds1
 +
|66.95
 +
|92.44
 +
|-
 +
|100-ds2
 +
|85.12
 +
|75.99
 +
|-
 +
|143-ds1
 +
|72.24
 +
|99.83
 +
|-
 +
|143-ds2
 +
|1.14
 +
|16.81
 +
|-
 +
|143-5
 +
|94.88
 +
|97.57
 +
|-
 +
|143-6
 +
|96.14
 +
|98.78
 +
|-
 +
|143-7
 +
|93.78
 +
|95.67
 +
|-
 +
|217-ds1
 +
|78.12
 +
|69.74
 +
|-
 +
|217-ds2
 +
|27.33
 +
|30.87
 +
|-
 +
|217-4
 +
|95.99
 +
|68.44
 +
|-
 +
|217-1
 +
|99.08
 +
|97.49
 +
|-
 +
|217-2
 +
|97.18
 +
|98.86
 +
|-
 +
|217-3
 +
|97.59
 +
|99.07
 +
|}
 +
 
 +
== Colour Correction of the Beam Shape ==
 +
 
 +
In general, the measured beam shape is a function of the spectral energy distribution (SED) of the measurement source. This is of particular concern because we measure the beam on a source with a roughly Rayleigh-Jeans SED, yet we use the effective beam window function to correct the CMB.
 +
 
 +
Planck Collaboration VII (2014) {{PlanckPapers|planck2013-p03c}} described possible levels of this bias using GRASP calculations with the pre-launch telescope model (Maffei et al. 2010 {{BibCite|maffei2010}}; Tauber et al. 2010 {{BibCite|tauber2010b}}). The prelaunch calculations did not agree well enough with the data to allow a direct application of the colour correction. A new telescope model based on flight data, and presented in a forthcoming paper, predicts beams that agree better with the data at 100–217 GHz, but show worse agreement at 353 GHz. The new model predicts a different beam shape colour correction, though at a similar order of magnitude to that shown in Planck Collaboration VII (2014) {{PlanckPapers|planck2013-p03c}}. This work is ongoing and will be completed after the 2015 release.
 +
 
 +
We search for a possible signature of the colour correction effect in the data. First, we check for consistency in the CMB angular power spectra derived from different detectors and frequency bands using the SMICA algorithm, in order to find discrepancies in beam shape and relative calibration of different data subsets of the CMB (Planck Collaboration A13 2014 {{PlanckPapers|planck2014-a13}}). There are hints of differences that are orthogonal to the beam error eigenmodes and appear as changes in relative calibration, but the preferred changes are at the 0.1 % level. Because this correction is of the same order of magnitude as the relative calibration uncertainty between detectors, we treat it as insignificant.
 +
 
 +
Second, we look at the relative calibration between detectors within a band for compact sources as a function of the source SED. In this method, any discrepancies due to solid angle variation with SED are degenerate with an error in the underlying bandpass. Additionally, because of the lack of bright, compact sources with red spectra, this method is limited in its signal-to-noise ratio. In the limit that any detected discrepancy is completely due to solid-angle variation with colour, at the level of 1 % in solid angle we do not detect variations consistent with the physical optics predictions. Given the lack of measurement of this effect at the current levels of uncertainty, as well as uncertainties in modelling of the telescope, we note that this effect may be present at a small level in the data, but do not attempt to correct for it or include it in the error budget.
 +
 
 +
== Effective beam window functions ==
 +
===Changes in QuickBeam ===
 +
The effective beam window functions <i>B(&#8467;)</i> for HFI, were computed using [[Beams#QuickBeam|QuickBeam]] as in 2013, with a few improvements.
 +
 
 +
==== Higher spin ====
 +
Since the scanning 2015 beams are slightly more elongated than the 2013 ones, because of different transfer function deconvolution, the QuickBeam approach was extended to include the beam coefficients <i>b<sub>&#8467;s</sub></i> up to <i>|s| = 6 </i>, instead of <i>|s| &le; 4</i>.
 +
We checked that, for <i>&#8467; &le; 3000</i> (or <i>&#8467; &le; 2000</i> at 100GHz), the coefficients with <i>|s| > 6</i> account for much less of 0.1% of the beam solid angle.
 +
 
 +
==== Impact of sky cut ====
 +
The effective beam window function will be affected by the finite area of the pixels on which the sky signal
 +
is averaged to produce the map. While the HEALPix pixels all have the same surface area
 +
for a given resolution parameter <i>N</i><sub>side</sub>, their shape differ and the trace of their moment of inertia:
 +
<p>
 +
<math>
 +
{\rm Tr} ({\bf I}_p) = \int d{\bf x} d{\bf y} \left[ ({\bf x}-{\bf x}_p)^2 + ({\bf y}-{\bf y}_p)^2 \right],
 +
</math>
 +
</p>
 +
where (<b>x</b>, <b>y</b>) are the tangential plane coordinates in the pixel <i>p</i> (valid for typical pixel size of 1.7' for <i>N</i><sub>side</sub>=2048) and (<b>x</b><sub><i>p</i></sub>,<b>y</b><sub><i>p</i></sub>) is the pixel centre,
 +
varies across the sky, especially in the polar regions, as illustrated in Fig. 6 below.
 +
[[Image:trace_inertia_R2048.png|thumb|500px|center|'''Figure 6:''' Map of the variation of the trace of the HEALPix pixel moment of inertia, for <i>N</i><sub>side</sub>=2048.]]
 +
The full-sky average of this moment of inertia determines the nominal HEALPix pixel window function, which is not included in the effective beam window function provided and must be corrected for separately.
 +
When the [[Frequency_Maps#Galactic_Plane_masks|Galactic masks]] used
 +
in Planck blank out the equatorial region, the window function of the remaining pixels will differ from the full-sky average, and will affect the resulting beam window functions.
 +
We followed a semi-analytical approach in QuickBeam, where the beam window functions are corrected from the departure between the full-sky average of the tensor of inertia and its average over un-masked pixels, and checked that the results agree
 +
very well with those of Monte Carlo simulations of the exact Planck scanning of the pixels performed with FEBeCoP, for several Galactic masks with <i>f</i><sub>sky</sub>=100%, 80%, and 40%.
 +
 
 +
The effective beam window function are therefore provided for two different sky fractions (100% and 75%) in different [[The_RIMO#Beam_Window_Functions|RIMOs]]. For these two sky fractions the square of the beam window functions (relevant for the power spectra analysis) are identical at &#8467;=1 and have a relative difference of about 0.4% at &#8467;=2000.
 +
 
 +
== Data Products ==
 +
 
 +
Beam products are available from the instrument model ([[The_RIMO]]) and further described in [[Scanning_Beams]] and [[Effective_Beams]].
 +
 
 +
Effective beam window functions are available in specific [[The_RIMO#Beam_Window_Functions|RIMOs]].
 +
They do <i>not</i>' contain the nominal HEALPix pixel window function.
 +
 
 +
= The 2013 Data Release =
 +
 
 +
== Scanning beams ==
 +
 
 +
The scanning beams describe the instrument’s instantaneous beam profile. Due to the near constant spin rate of the spacecraft, time domain effects (including residual time response and lowpass filtering) are degenerate with the spatial response due to the optical system. The scanning beam reconstruction recovers both of these effects, except for residual time-domain effects on a longer timescale than can be captured with the extent of the scanning beam model.
 +
 
 +
In {{PlanckPapers|planck2013-p03c}} we consider two models of the beam in order to better understand systematic effects in the reconstruction. Here we describe only the "B-spline beams," which are used to compute the delivered effective beam (see next section).
 +
 
 +
=== B-spline Beam construction ===
 +
 
 +
We use seasons 1 and 2 of the Mars observations to reconstruct the beam. The data are processed with the "bigPlanets" TOI processing (see [[Frequency_Maps#HFI_processing | here]] or section 3.11 of {{PlanckPapers|planck2013-p03}} for more information). We use JPL Horizons ephemerides to determine the pointing of each detector relative to the planet. We then subtract the astrophysical background in the time domain using a bicubic interpolation of the Planck maps.
 +
 
 +
The time-ordered data are used to fit a two dimensional B-spline surface using a least squares minimization and a smoothing criterion to minimize the effects of high spatial frequency variations. We therefore assume the scanning beam to be smooth. The smoothing criterion, as well as the locations of the nodes used to compute the B-spline basis functions, are set using GRASP physical optics simulations as inputs, which are the best assumptions on the spatial frequency content of the in-flight beams.
  
 
The smoothing criterion is defined as follows:
 
The smoothing criterion is defined as follows:
  
<math>\eta = \displaystyle{\sum_{i=1}^{g}\left(b^{k}(\lambda_{i+})-b^{k}(\lambda_{i-})\right)^2}
+
<math>\eta = \displaystyle{\sum_{i=1}^{g}\left(b^{k}(\lambda_{i+})-b^{k}(\lambda_{i-})\right)^2},
 
\label{smoothcrit}</math>
 
\label{smoothcrit}</math>
  
<math>\begin{aligned}
+
where &eta; is the Smoothing Criterion and
\eta &: \mbox{ Smoothing Criterion}\\
+
<i>b<sup>k</sup></i> represents the <i>k</i>th beam derivative evaluated at the node locations.
b^k &: \mbox{ $k^{th}$ beam derivative evaluated on the nodes locations}\end{aligned}</math>
 
  
And the global inversion criterion :
+
The global inversion criterion is:
  
<math>\zeta = \eta + p\times \delta</math>
+
<math>\zeta = \eta + p\times \delta</math>,
  
with <math>\delta</math> usual least square estimator and <math>p</math> coefficient giving the relative weight to <math>\delta</math> with respect to the smoothing criterion.
+
with &delta; being the usual least squares estimator and the <i>p</i> coefficient giving the relative weight of &delta; with respect to the smoothing criterion
  
<math>\delta = \displaystyle{\sum_{r=1}^{m}}\left(y_{r} - b(x_{r})\right)^2\label{estimator}</math>
+
<math>\delta = \displaystyle{\sum_{r=1}^{m}}\left(y_{r} - b(x_{r})\right)^2\label{estimator}</math>.
  
<math>\begin{aligned}
+
Here &delta; is the usual least squares criterion,
\delta &: \mbox{ usual least square criterion}\\
+
<i>r</i> is the index relative to the <i>m</i> data points (<i>r</i> in the range 1 to <i>m</i>),
r &: \mbox{ indice relative to the m data points, } r \in \{1, \ldots, m\}\\
+
<i>y<sub>r</sub></i> is the planet data of sample <i>r</i>,
y_r &: \mbox{ planet data of sample r}\\
+
<i>x<sub>r</sub></i> is the pointing of sample <i>r</i>,
x_r &: \mbox{ pointing of sample r}\\
+
and <i>b</i> is the reconstructed beam.
b &: \mbox{ reconstructed beam}\end{aligned}</math>
 
  
The B-Spline nodes are located on a regular spaced grid in the detector coordinate framset. At the edge of the reconstructed beam map area, 4 coincident nodes are added to avoid vanishing basis functions.
+
The B-spline nodes are located on a regularly spaced grid in the detector coordinate frameset. At the edge of the reconstructed beam map area, four coincident nodes are added to avoid vanishing basis functions.
  
Let <math>B_{i, k+1}</math>, <math>k</math> degree B-Spline build using nodes {<math>\lambda_{i}, ..., \lambda_{i+k+1}</math>} (''De Boor &amp; Cox'', 1972) :
+
Let <i>B<sub>i, k+1</sub></i> be a <i>k</i>-degree B-spline built using nodes {&lambda;<i><sub>i</sub></i>, &hellip;, &lambda;<i><sub>i+k+1</sub></i>} (the De Boor-Cox method), then
  
 
<math>B_{i,1}(x) = \left\{
 
<math>B_{i,1}(x) = \left\{
 
     \begin{array}{l}
 
     \begin{array}{l}
     1, \mbox{ si } x \in \mbox{[} \lambda_{i}, \lambda_{i+1} \mbox{[}\\
+
     1, \mbox{ if } x \in \mbox{[} \lambda_{i}, \lambda_{i+1} \mbox{]},\\
     0, \mbox{ si } x \notin \mbox{[} \lambda_{i}, \lambda_{i+1} \mbox{[}
+
     0, \mbox{ if } x \notin \mbox{[} \lambda_{i}, \lambda_{i+1} \mbox{]},
 
     \end{array} \right.</math>
 
     \end{array} \right.</math>
  
<math>B_{i, l+1}(x) = \displaystyle{\frac{x - \lambda_{i}}{\lambda_{i+l} - \lambda_{i}}} B_{i,l}(x) + \displaystyle{\frac{\lambda_{i+l+1}-x}{\lambda_{i+l+1}-\lambda_{i+1}}} B_{i+1, l}(x)</math>
+
<math>B_{i, l+1}(x) = \displaystyle{\frac{x - \lambda_{i}}{\lambda_{i+l} - \lambda_{i}}} B_{i,l}(x) + \displaystyle{\frac{\lambda_{i+l+1}-x}{\lambda_{i+l+1}-\lambda_{i+1}}} B_{i+1, l}(x)</math>,
  
<math>l=1, \ldots, k</math>
+
<math>l=1, \ldots, k</math>.
  
  
[[Image:FocalPlane_Map_BSScanningBeams_v53.png|800px|frame|center|Focal plane plot of B-Spline scanning beams using in-flight pointing reconstruction.  The contours are -3,-10,-20,-30 dB from the peak, and for PSB pairs the "a" bolometer is plotted in black and "b" in blue.]]
+
[[Image:FocalPlane_Map_BSScanningBeams_v53.png|thumb|500px|center|'''Figure 1:''' Focal plane plot of B-spline scanning beams using in-flight pointing reconstruction.  The contours are -3, -10, -20, and -30 dB from the peak, and for PSB pairs the "a" bolometer is plotted in black and the "b" in blue.]]
  
 
=== Simulations and errors ===
 
=== Simulations and errors ===
  
We estimate the reconstruction bias and noise in the measurements using an ensemble of simulated planet observations for each channel. Kept fixed in each simulation are:
+
We estimate the reconstruction bias and noise in the measurements using an ensemble of simulated planet observations for each channel. Further details are discussed in {{PlanckPapers|planck2013-p03c}}. Kept fixed in each simulation are:
  
* the input beam assumed: we use a supersampled version of the reconstructed B-Spline beam (or whatever comes out of the current ongoing tests!)
+
* the input beam assumed, specifically we use a over-sampled version of the reconstructed B-Spline beam (or whatever comes out of the current ongoing tests!);
* Astrophyical background is the same as that subtracted from the real data.
+
* the astrophyical background is the same as that subtracted from the real data;
* StarTracker pointing (using the ptcor6 pointing model).
+
* StarTracker pointing (using the "ptcor6" [[Detector pointing|pointing model]]).
  
 
The following are varied in each simulation:
 
The following are varied in each simulation:
  
* detector noise realizations obtained by filtering randomly generated white noise with the measured noise PSDs
+
* detector noise realizations, obtained by filtering randomly -generated white noise with the measured noise PSDs;
* random pointing errors with 2 arcsecond rms, and a spectrum that replicates the real errors.
+
* random pointing errors with 2 arcsec rms, and a spectrum that replicates the real errors;
* simulated glitches and the deglitching procedure
+
* simulated glitches and the deglitching procedure;
* Mars brightness temperature variability
+
* Mars brightness temperature variability.
  
400 simulated timelines are generated for each bolometer and for each of the two seasons of Mars observations used in the beam reconstruction. The simulated timelines are made into beam maps, projecting onto the B-Spline basis in the same way as the real data.
+
400 simulated timelines are generated for each bolometer and for each of the two seasons of Mars observations used in the beam reconstruction. The simulated timelines are made into beam maps, projecting onto the B-Spline basis in the same way as for the real data.
  
The beam maps are propagated to effective beam window functions using the quickbeam approach (see effective beams below) and used to evaluate the reconstruction bias and to construct error eigenmodes in the effective beam window function.
+
The beam maps are propagated to effective beam window functions using the "QuickBeam" approach (see effective beams below) and used to evaluate the reconstruction bias and to construct error eigenmodes in the effective beam window function.
  
<font color=red>Figure: random pointing error PSD Figures: error envelope plots (or should those go under effective beams?)</font>
+
<!-- <font color=red>Figure: random pointing error PSD Figures: error envelope plots (or should those go under effective beams?)</font> -->
  
 
=== Residuals ===
 
=== Residuals ===
  
There are two known beam effects that are not included in the main beam model and are estimated as a separate bias in flux and angular power spectrum measurement: 1. long tails due to errors in low frequency time response deconvolution, and 2. near sidelobes.
+
There are two known beam effects that are not included in the main beam model and are estimated as a separate bias in flux and angular power spectrum measurement: (1) long tails due to errors in low frequency time response deconvolution; and (2) near sidelobes.
  
We stack all five observations of Jupiter to estimate the long time scale residuals due to incomplete deconvolution of the long time scale response.
+
We stack all five observations of Jupiter to estimate the long timescale residuals due to incomplete deconvolution of the long timescale response.
  
<font color=red>Add some kind of mean tail plot</font>
+
<!-- <font color=red>Add some kind of mean tail plot</font>
  
 
Near sidelobes are also evaluated using stacked Jupiter (hopefully they will just be part of the v53bis B-Spline beams). The main features in the near sidelobes include a wide beam skirt, and dimpling lobes
 
Near sidelobes are also evaluated using stacked Jupiter (hopefully they will just be part of the v53bis B-Spline beams). The main features in the near sidelobes include a wide beam skirt, and dimpling lobes
<font color=red>Add sidelobe plots and tables</font>
+
<font color=red>Add sidelobe plots and tables</font> -->
 +
 
 +
== Effect of mirror support print-through  ==
 +
 
 +
The Planck reflectors suffer from print-through
 +
of the
 +
honeycomb structures that support the carbon-fibre face sheets.
 +
While the size of the deformation has been measured during tests below room temperature to be less than 20&mu;m, it is the strict periodicity of the deformation that contributes most to the additional beam-shape change.  A simple grating equation has been shown to describe the angular positions of the resulting near side lobes very well:
 +
 
 +
 
 +
<math>  \displaystyle{ \sin \theta _n = \frac{n \lambda}{Yd} } \label{DimplingLobes}</math>,
 +
 
 +
where:
 +
* &theta;<sub><i>n</i></sub> is the angular position of the <i>n</i>th-order lobe from the central beam peak;
 +
* &lambda; is the wavelength of the radiation;
 +
* <i>d</i> is the grating spacing of the periodicity;
 +
* <i>Y</i> is a factor that describes the position of the each reflector along the optical path;
 +
* <i>Y</i>=1.00 for the primary reflector;
 +
* <i>Y</i>=1.80 for the secondary reflector.
 +
 
 +
 
 +
Three possible periodicities (19.6mm, 30mm, 52mm) in the honeycomb array dominate the Planck dimpling pattern for the 857-GHz detectors, although 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 arcminute beam window, but they contribute at most (0.05 ± 0.01)% to the integrated beam.
 +
 
 +
A map of Jupiter has been created for each 857-GHz detector, using the first four surveys, with the background subtracted using the same sky area for the survey taken 12 months before, in which the planet is not present. The background-subtracted maps for each survey were stacked to make a single map for each of the four detectors and these were further stacked to create a single 857-GHz band map, using the  standard detector weightings proportional to 1/(NET)<sup>2</sup>. All the maps have 96 pixels per side,  covering ±1.4º, and are scaled to the RMS noise of the background. Central beam values range from 205,000 to 256,000, depending on the noise of the detector.
 +
 
 +
For each of these five maps of Jupiter an elliptical Gaussian main beam was fitted, and the amount of saturation estimated by finding the peak increase that gives the best Gaussian fit, typically 10-20%. A circular ruze envelope (see {{PlanckPapers|planck2013-p03c}}) was also fitted, together with a tilted Gaussian component slightly offset from the beam in the cross-scan (optical <i>x</i>-axis) direction, using data more than 0.13&deg; from the beam, with the known positions of the dimpling lobes masked out. Once these strongest components were removed, a tilted elliptical Gaussian was fitted for each visible dimpling lobe in up to five different sets of lobes; however, the innermost lobe set (52-mm periodicity) is generally obscured by the ruze component and for the outermost set (19.6-mm periodicity) most lobes are too weak to be fitted for individual detectors. The contributing areas of the lobes were determined by three methods: from the fitted Gaussians; from summation within boxed areas, using offsets determined from surrounding boxed areas; and from summation in the boxed areas using an offset determined from Gaussian fitting.
 +
 
 +
The uncertainties in the beam component areas given here include the spread in values from these different methods. Similarly, the lobe peaks in decibels are calculated using the raw, fitted and estimated desaturated beam peak values, and the uncertainties reflect the variations produced by the three methods. Averaged over the five map fits, the dimpling lobes contribute  (0.47 ± 0.13)% of the beam area, while the ruze envelope accounts for (9.2 ± 0.7)% of the beam, with the main beam making up the remaining (90.3 ± 0.7)%.
  
 +
For lower frequencies the dimpling lobes are less visible and will not be discussed here. For example, the noise floor for the 545-GHz band map is about 47 dB below the observed Jupiter peak, so that no dimpling lobes are visible further out than the 52mm lobes at 0.61&deg;. For 353 GHz the noise floor is at about 38 dB so that even the 52-mm lobes are not observed.
 +
 +
[[image:CompleteFitWithPNContoursV53bis.jpg|700px|thumb|center|'''Figure 2:''' Stacked Jupiter map for the 857-GHz detectors. Left: the data with  main beam and ruze envelope subtracted. Right: the fitted elliptical Gaussian dimpling lobes with GRASP contours plotted on top.]]
 +
 +
 +
The figure above shows how the dimpling lobes seen for the 857GHz band Jupiter map correspond to the contours calculated by the physical optics GRASP package (produced by TICRA, Denmark). The GRASP simulation was performed for the 857-1 detector and assumed a uniform dimpling distortion of 10&mu;m. The biggest departure from the simple grating model is that the lobe pattern is elongated in the cross-scan (vertical) direction&mdash;lobes that should lie on a line 30&deg; from the vertical are consistently found around 25&deg; and those on the 60&deg; line cluster around 53&deg;. This is due to the offset geometry of the mirror system, whereby the dimpling print-though is foreshortened along the optical <i>x</i>-axis due to the tilt of the mirrors, producing greater lobe spacing as seen by the incoming photons. The amplitudes of the fitted dimpling lobes vary significantly within each set, whereas the GRASP model shows a constant amplitude for lobes within a set (except for the 30-mm periodicity). Similar amplitudes are seen in lobes lying directly across the beam centre from each other, indicating that "facesheet dimpling" has not occurred uniformly in all directions, as the GRASP model assumes.
 +
 +
For the 30-mm periodicity the vertical lobes are significantly weaker than those at the sides, just as the GRASP model shows. Generally, the fitted lobes are somewhat stronger than those seen by GRASP, indicating that the dimpling is indeed larger than 10&mu;m. While the dimpling lobes are measureable at the highest two frequencies,  with a strong source (i.e., Jupiter), 90% of the dimpling lobe area at 857GHz, and 10% for all other frequencies, already appears in the beam functions  described above and needs no further correction. However, the presence of the uniform dimpling lobes provides a useful window by which to investigate the mirror geometry.
  
 
== Effective Beams ==
 
== Effective Beams ==
  
The '''effective beam''' is the average of all scanning beams pointing at a certain direction within a given pixel of the sky map for a given scan strategy. It takes into account the coupling between azimuthal asymmetry of the beam and the uneven distribution of scanning angles across the sky. It captures the complete information about the difference between the true and observed image of the sky. They are, by definition, the objects whose convolution with the true CMB sky produce the observed sky map.  
+
The ''effective beam'' is the average of all scanning beams pointing at a certain direction within a given pixel of the sky map for a given scan strategy. It takes into account the coupling between azimuthal asymmetry of the beam and the uneven distribution of scanning angles across the sky. It captures the complete information about the difference between the true and observed images of the sky. They are, by definition, the objects whose convolution with the true CMB sky produce the observed sky map.  
  
Several methods of effective beams determination have been developped and cross-validated.
+
Several methods of effective beam determination have been developed and cross-validated. The main products are produced using FEBeCoP and details of the processing are given in the [[Effective Beams]] products page. See also the equivalent [[Beams_LFI | page discussing the LFI beams]]
  
The main products are produced using FEBeCoP and details of the processing are given in the [[Effective Beams]] products page. See also the equivalent [[Beams_LFI | page discussing the LFI beams]]
 
 
<font color=red>Need satisfactory comparison plot</font>
 
  
 
=== FEBeCoP ===
 
=== FEBeCoP ===
 
-------------------
 
-------------------
  
The full algebra for this method for the calculation of effective beams was presented in [[http://arxiv.org/pdf/1005.1929| Mitra, Rocha, Gorski et al.]] <cite>#mitra2010</cite>. Here we summarise the main results. The observed temperature sky <math>\widetilde{\mathbf{T}} </math> is a convolution of the true sky <math>\mathbf{T} </math> and the effective beam <math>\mathbf{B}</math>:
+
The full algebra for this method for the calculation of effective beams was presented in Ref.{{BibCite|mitra2010}}. Here we summarize the main results. The observed temperature sky <math>\widetilde{\mathbf{T}} </math> is a convolution of the true sky <math>\mathbf{T} </math> and the effective beam <math>\mathbf{B}</math>:
  
 
<math>
 
<math>
Line 108: Line 350:
  
 
<math>
 
<math>
B_{ij} \ = \ \left( \sum_t A_{ti} \, b(\hat{\mathbf{r}}_j, \hat{\mathbf{p}}_t) \right) / \left({\sum_t A_{ti}} \right) \, ,
+
B_{ij} \ = \ \left[ \sum_t A_{ti} \, b(\hat{\mathbf{r}}_j, \hat{\mathbf{p}}_t) \right] \bigg/ \left[{\sum_t A_{ti}} \right] \, .
 
\label{eq:EBT2}
 
\label{eq:EBT2}
 
</math>
 
</math>
  
<math>t</math> is time samples, <math>A_{ti}</math> is <math>1</math> if the pointing direction falls in pixel number <math>i</math>, else it is <math>0</math>, <math>\mathbf{p}_t</math> represents the exact pointing direction (not approximated by the pixel centre location), and <math>\hat{\mathbf{r}}_j</math> is the centre of the pixel number <math>j</math>, where the scanbeam <math>b(\hat{\mathbf{r}}_j, \hat{\mathbf{p}}_t)</math> is being evaluated (if the pointing direction falls within the cut-off radius of <math>\sim 2.5 \times</math> FWHM.
+
Here <math>t</math> labels time samples, <math>A_{ti}</math> is 1 if the pointing direction falls in pixel number <i>i</i>, otherwise it is 0, <math>\mathbf{p}_t</math> represents the exact pointing direction (not approximated by the pixel centre location), and <math>\hat{\mathbf{r}}_j</math> is the centre of pixel number <i>j</i>, where the scanning beam <math>b(\hat{\mathbf{r}}_j, \hat{\mathbf{p}}_t)</math> is being evaluated if the pointing direction falls within the cut-off radius of about 2.5 &times; FWHM.
  
The algebra is a bit more involved for polarised detectors. The observed stokes parameters at a  pixel <math>i</math>, <math>(\widetilde{I}, \widetilde{Q}, \widetilde{U})_i</math>, are related to the true stokes parameters <math>(I, Q, U)_i</math>, by the following relation:
+
The algebra is a bit more involved for polarized detectors. The observed Stokes parameters at a  pixel <i>i</i>, <math>(\widetilde{I}, \widetilde{Q}, \widetilde{U})</math><sub><i>i</i></sub>, are related to the true Stokes parameters <math>(I, Q, U)</math><sub><i>i</i></sub>, by the following relation:
  
 
<math>
 
<math>
( \widetilde{I} \quad \widetilde{Q} \quad \widetilde{U})_i^T \ = \ \Delta\Omega \sum_j \mathbf{B}_{ij} \cdot (I \quad Q \quad U)_j^T,
+
( \widetilde{I} \quad \widetilde{Q} \quad \widetilde{U})_i^{\rm T} \ = \ \Delta\Omega \sum_j \mathbf{B}_{ij} \cdot (I \quad Q \quad U)_j^{\rm T},
 
\label{eq:a1}
 
\label{eq:a1}
 
</math>
 
</math>
  
where the polarised effective beam matrix
+
where the polarized effective beam matrix
  
 
<math>
 
<math>
\mathbf{B}_{ij} \ = \  \left[ \sum_t A_{tp} \mathbf{w}_t \mathbf{w}^T_t \right]^{-1} \sum_t A_{ti} \, b(\hat{\mathbf{r}}_j, \hat{\mathbf{p}}_t) \, \mathbf{w}_t  \mathbf{W}^T(\hat{\mathbf{n}}_j,\hat{\mathbf{p}}_t) \, ,
+
\mathbf{B}_{ij} \ = \  \left[ \sum_t A_{tp} \mathbf{w}_t \mathbf{w}^{\rm T}_t \right]^{-1} \sum_t A_{ti} \, b(\hat{\mathbf{r}}_j, \hat{\mathbf{p}}_t) \, \mathbf{w}_t  \mathbf{W}^{\rm T}(\hat{\mathbf{n}}_j,\hat{\mathbf{p}}_t) \, ,
 
\label{eq:a2}
 
\label{eq:a2}
 
</math>
 
</math>
  
and <math>\mathbf{w}_t </math>and <math>\mathbf{W}(\hat{\mathbf{r}}_j, \hat{\mathbf{p}}_t) </math> are the the polarisation weight vectors, as defined in \cite{mitra2010}.
+
and <math>\mathbf{w}_t</math>and <math>\mathbf{W}(\hat{\mathbf{r}}_j, \hat{\mathbf{p}}_t)</math> are the the polarization weight vectors, as defined in Ref.{{BibCite|mitra2010}}.
  
The task is to compute <math>B_{ij}</math> for temperature only beams and the <math>3 \times 3</math> matrices <math>\mathbf{B}_{ij}</math> for each pixel <math>i</math>, at every neighbouring pixel <math>j</math> that fall within the cut-off radius around the the center of the <math>i^\text{th}</math> pixel.
+
The task is to compute <math>B_{ij}</math> for temperature-only beams and the 3&times;3 matrices <math>\mathbf{B}_{ij}</math> for each pixel <i>i</i>, at every neighbouring pixel <i>j</i> that falls within the cut-off radius around the the centre of the <i>i</i>th pixel.
  
The effective beam is computed by stacking within a small field around each pixel of the HEALPix sky map.  Due to the particular features of Planck scanning strategy coupled to the beam asymmetries in the focal plane, and data processing of the bolometer and radiometer TOIs, the resulting Planck effective beams vary over the sky.  
+
The effective beam is computed by stacking within a small field around each pixel of the HEALPix sky map.  Due to the particular features of Planck's scanning strategy, coupled to the beam asymmetries in the focal plane, and data processing of the bolometer and radiometer TOIs, the resulting Planck effective beams vary over the sky.  
  
FEBeCoP, given information on Planck scanning beams and detector pointing during a mission period of interest, provides the pixelized stamps of both the Effective Beam, EB, and the Point Spread Function, PSF, at all positions of the HEALPix-formatted map pixel centres.
+
FEBeCoP, given information on the Planck scanning beams and detector pointing during a mission period of interest, provides pixelized stamps of both the Effective Beam and the Point Spread Function, at all positions of the HEALPix-formatted map pixel centres.
 +
 
 +
=== FICSBell ===
  
  
=== FICSBell ===
 
For more details, see <cite>#planck2013-p03c</cite>.
 
  
Since the HFI beams are not azimuthally symmetric, the scanning strategy has to be taken into account in the effective beam response modelling. This is done using the FICSBell method <font color=red>(Hivon et al, in preparation)</font>, which generalizes to polarization and to include other sources of systematics the approach used for TT <math>C(l)</math> estimation in WMAP-3yr <font color=red>Hinshaw et al (2007)</font> and by <font color=red>Smith et al (2007)</font> in the detection of CMB lensing in WMAP maps. The different steps of the method used for this study can be summarized as follows:
+
Since the HFI beams are not azimuthally symmetric, the scanning strategy has to be taken into account in the effective beam response modelling (for more details, see {{PlanckPapers|planck2013-p03c}}). This is done using the "FICSBell" method (Hivon et al., in preparation).  This generalizes to polarization and to include other sources of systematics the approach used for <i>TT</i> <i>C&#8467;</i> estimation in WMAP-3
 +
{{BibCite|hinshaw2007}} and by Smith et al (2007) {{BibCite|smith2007}}
 +
in the detection of CMB lensing in WMAP maps. The different steps of the method used for this study can be summarized as follows.
  
 
<ol>
 
<ol>
<li><p>The scanning related information (i.e., statistics of the orientation of each detector within each pixel) is computed first, and only once for a given observation campaign. Those orientation hit moments are only computed up to degree 4, for reasons described in point 2 below. At the same time, the first two moments of the distribution of samples within each pixel (ie, their center of mass and moments of inertia) are computed and stored on disc.</p></li>
+
<!--  1 -->
<li><p>The scanning beam map or beam model of each detector <math>d</math> is analyzed into its Spherical Harmonics coefficients</p>
+
<li><p>The scanning-related information (i.e., statistics of the orientation of each detector within each pixel) is computed first, and only once for a given observation campaign:</p>
<p><math>b^d_{ls} = \int d{\bf r} B_d({\bf r}) Y_{ls}({\bf r})\label{scanningBlm}</math></p>
+
<p><math>
<p>where <math>B_d(\bf{r})</math> is the beam map centered on the North pole, and <math>Y_{ls}(\bf{r})</math> is the Spherical Harmonics basis function. Higher <math>s</math> indexes describes higher degrees of departure from azimuthal symmetry and, for HFI beams, the coefficients <math>b^d_{ls}</math> are decreasing functions of <math>s</math> at most <math>l</math> considered. It also appears that, for <math>l<3000</math>, the coefficients with <math>|s| > 4</math> account for <math>1\%</math> or less of the beam throughput. For this reason, only modes with <math>|s| \le 4</math> are considered in the present analysis. <font color=red>Armitage-Caplan and Wandelt (2009)</font> reached a similar conclusion in their deconvolution of Planck-LFI beams.</p></li><li><p>The <math>b^d_{ls}</math> coefficients computed above are used to generate <math>s</math>-spin weighted maps, as well as the first and second order derivatives, for a given CMB sky realization.</p></li>
+
w^d_s({\bf r}_p) = \sum_j e^{i s \psi_j},  \label{eq:e1}
<li><p>The spin weighted maps and orientation hit moments of the same order <math>s</math> are combined for all detectors involved, to provide an “observed” map. 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 [sec:pixelization]). In this combination, the respective number of hits of each detector in each pixel is considered, as well as the weighting (generally proportional to the inverse noise variance) applied to each detector in order to minimize the final noise.</p></li><li><p>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.</p></li></ol>
+
</math></p>
Monte-Carlo (MC) simulations in which the sky realisations 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 into the MC simulations.
+
<p>where &psi;<sub><i>j</i></sub> is the orientation of the detector with respect to the local
 +
meridian during the measurement <i>j</i> occurring in
 +
the direction <b>r</b><sub>p</sub>. Note that the <i>s</i>=0 moment is simply the
 +
hit count map. The orientation hit moments are computed up to
 +
degree <i>s</i>=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 on disc.
 +
 
 +
<!--  2 -->
 +
<li>The scanning-beam map or beam model of each detector <i>d</i> is analysed into its spherical harmonic coefficients</p>
 +
<p><math>
 +
  b^d_{\ell s} = \int d{\bf r} B_d({\bf r}) Y_{\ell s}({\bf r})\label{scanningBlm}
 +
</math></p>
 +
<p>where <i>B<sub>d</sub></i>(<b>r</b>) is the beam map centred on the North pole, and <i>Y<sub>&#8467;s</sub></i>(<b>r</b>) is the spherical harmonics basis function. Higher <i>s</i> indices describe higher degrees of departure from azimuthal symmetry and, for HFI beams, the coefficients <i>b<sup>d</sup><sub>&#8467;s</sub></i> are decreasing functions of <i>s</i> at most multipoles considered.  
 +
Spot checks, where window functions are computed with
 +
|<i>s</i>|&le;6, show a difference of less than 10<sup>-4</sup> for
 +
&#8467;<2000 at 100 GHz and for &#8467;<3000 at 143 and 217 GHz. For these reasons, only modes with |<i>s</i>|&le;4 are considered in the present analysis.
 +
Armitage-Caplan & Wandelt (2009){{BibCite|armitage-caplan2009}} reached a similar conclusion in their deconvolution of LFI beams.
 +
For these reasons, only modes with <math>|s| \le 4</math> are considered in the present analysis.
 +
 
 +
<!--  3 -->
 +
<li>
 +
For a given CMB sky realization <i>t</i>, described by its spherical harmonic coefficients
 +
<i>a<sub>&#8467;m</sub></i> = &int; <i>d</i><b>r</b> <i>t</i>(<b>r</b>) <i>Y<sub>&#8467;m</sub></i>(<b>r</b>), the <i>b<sup>d</sup><sub>&#8467;s</sub></i> coefficients computed above are
 +
used to generate <i>s</i>-spin weighted maps, </p>
 +
<p><math>
 +
m^d_{s}({\bf r}) = \sum_{\ell m} b^d_{\ell s}\  a_{\ell m}\ {}_sY_{\ell m}({\bf r}), \label{eq:e3}
 +
</math></p>
 +
<p>as well as the first and second derivatives, using standard HEALPix tools.
 +
 
 +
<!--  4 -->
 +
<li>The spin-weighted maps and orientation hit moments of the same order <i>s</i> are combined for all detectors involved, to provide an “observed” map</p>
 +
<p><math>
 +
m({\bf r}) = \left(\sum_d \sum_s  w^d_s({\bf r}) m^d_s({\bf r})\right) \bigg/ \sum_d  w^d_0({\bf r}). \label{eq:e4}
 +
</math></p>
 +
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 [[#Pixelization_Artifacts|Pixelization Artifacts section]]). In this combination, the respective number of hits of each detector in each pixel is considered, as well as the weighting (generally proportional to the inverse noise variance) applied to each detector in order to minimize the final noise.
 +
 
 +
<!-- 5 -->
 +
<li>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.
 +
</li></ol>
 +
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 into the MC simulations.
  
 
=== QuickBeam ===
 
=== QuickBeam ===
For more details, see <cite>#planck2013-p03c</cite>
+
We now describe the "QuickBeam" approach (for more details, see {{PlanckPapers|planck2013-p03c}}).
 +
==== Formalism ====
  
Planck observes the sky after convolution with a “scanning beam”, which captures its effective response to the sky as a function of displacement from the nominal pointing direction. Decomposing the scanning beam into harmonic coefficients <math>B_{lm}</math>, each time-ordered data (TOD) sample can be modelled as (neglecting the contribution from instrumental noise, which is independent of beam asymmetry) <math>%T_i = \sum_{lms} D^{l}_{-m s} (\phi_i, \theta_i, \alpha_i) b_{ls} (-1)^{m) T_{lm} + n_i,
+
Planck observes the sky after convolution with a “scanning beam,which captures its effective response to the sky as a function of displacement from the nominal pointing direction. Decomposing the scanning beam into harmonic coefficients <i>B<sub>&#8467;m</sub></i>, each time-ordered data (TOD) sample can be modelled as (neglecting the contribution from instrumental noise, which is independent of beam asymmetry)  
T_i = \sum_{lms} e^{-i s \alpha_i} B_{ls} \tilde{T}_{lm} {}_s Y_{lm}(\theta_i, \phi_i),
+
<p><math>
\label{eqn:tod_beam}</math> where the TOD samples are indexed by <math>i</math>, and <math>\tilde{T}_{lm}</math> is the underlying sky signal. The spin spherical harmonic <math>{}_s Y_{lm}</math> rotates the scanning beam to the pointing location <math>(\theta, \phi)</math>, while the <math>e^{-i s \alpha_i}</math> factor gives it the correct orientation. Eq.  may be evaluated with the “TotalConvolver” algorithm of <font color=red>Wandelt and Gorski (2001)</font>, accelerated using the “conviqt” recursion relations <font color=red>Prezeau and Reinecke (2010)</font> This approach is implemented in LevelS.
+
%T_i = \sum_{\ell ms} D^{l}_{-m s} (\phi_i, \theta_i, \alpha_i) b_{\ell s} (-1)^{m) T_{\ell m} + n_i,
</ref>, although because it involves working with a TOD-sized objected it is necessarily slow.
+
T_i = \sum_{\ell ms} e^{-{\rm i} s \alpha_i} B_{\ell s} \tilde{T}_{\ell m} {}_s Y_{\ell m}(\theta_i, \phi_i),
 +
\label{eqn:tod_beam}
 +
</math></p>
 +
where the TOD samples are indexed by <i>i</i>, and <math>\tilde{T}_{\ell m}</math> is the underlying sky signal. The spin spherical harmonic <i><sub>s</sub>Y<sub>&#8467;m</sub></i> rotates the scanning beam to the pointing location (&theta;, &phi;), while the <i>e</i><sup>-i<i>s</i>&alpha;<sub><i>i</i></sub></sup> factor gives it the correct orientation.  The above equation may be evaluated with the “TotalConvolver” algorithm of Wandelt and Gorski (2001){{BibCite|wandelt2001}}, accelerated using the “conviqt” recursion relations of Prezeau and Reinecke (2010){{BibCite|prezeau2010}}. This approach is implemented in Level S, although because it involves working with a TOD-sized object, it 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 TOD samples is essentially a process of binning: <math>T(p) = \sum_{i \in p} T_i / H(p),
+
On small angular scales comparable to the size of the beam, it is a good approximation to assume that the procedure of mapmaking from TOD samples is essentially a process of binning:  
\label{eqn:map_beam_full}</math> where <math>H(p)</math> is the total number of hits in pixel <math>\hat{n}</math>.
+
<math>T(p) = \sum_{i \in p} T_i / H(p),\label{eqn:map_beam_full}</math>
 +
where <i>H(p)</i> is the total number of hits in pixel <math>\hat{n}</math>.
  
Start with a normalized, rescaled harmonic transform of the beam <math>B_{lm}</math>, sky multipoles <math>\tilde{T}_{lm}</math> and a scan history object <math>w(\hat{n}, s)</math> given by <math>w(\hat{n}, s) = \sum_{j \in p} e^{i s \alpha_j} / H(\hat{n})</math> where the sum is over all hits <math>j</math> of pixel <math>p</math> at location <math>\hat{n}_p</math>, and <math>\alpha_j</math> is the scan angle for observation <math>j</math>. The harmonic transform of this scan-strategy object is given by <math>{}_{s} w_{L M} = \int d^2 \hat{n} {}_s Y_{LM}^*(\hat{n}) w(\hat{n}, s).</math> The beam-convolved observation is then given by <math>\tilde{T}(\hat{n}) = \sum_{slm} w(\hat{n}, -s ) B_{ls} T_{lm} {}_s Y_{lm}(\hat{n}).</math> Taking the ensemble average of the pseudo-Cl power spectrum of these <math>T_{lm}</math> we find
+
Starting with a normalized, rescaled harmonic transform of the beam <i>B<sub>&#8467;m</sub></i>, sky multipoles <math>\tilde{T}_{\ell m}</math> and a scan history object <math>w(\hat{n}, s)</math> given by <math>w(\hat{n}, s) = \sum_{j \in p} e^{i s \alpha_j} / H(\hat{n})</math>, where the sum is over all hits <i>j</i> of pixel <i>p</i> at location <math>\hat{n}_p</math>, and &alpha;<sub>j</sub> is the scan angle for observation <i>j</i>. The harmonic transform of this scan-strategy object is given by
 +
<p><math>
 +
{}_{s} w_{L M} = \int d^2 \hat{n} {}_s Y_{LM}^*(\hat{n}) w(\hat{n}, s).
 +
\label{eq:spin_hits}
 +
</math></p><p>
 +
The beam-convolved observation is then given by  
 +
</p><p><math>
 +
\tilde{T}(\hat{n}) = \sum_{slm} w(\hat{n}, -s ) B_{\ell s} T_{\ell m} {}_s Y_{\ell m}(\hat{n}).
 +
\label{eq:tilde_T}
 +
</math></p><p>
  
<math>\begin{gathered}
+
Taking the ensemble average of the pseudo-<i>C</i><sub>&#8467;</sub> power spectrum of these <i>T<sub>&#8467;m</sub></i>, we find
\tilde{C}_{L}^{TT} = \sum_{S S'} \sum_{l_1 l_2} \frac{(2l_1+1)(2l_2+1)}{4\pi}
+
 
{}_{(-s -s')}{\cal W}_{l_1} B_{l_2 S} B_{l_2 S'}^* C^{TT}_{l_2}
+
</p><p><math>
\\ \times\left(
+
\tilde{C}_{L}^{TT} = \sum_{S S'} \sum_{\ell_1 ell_2} \frac{(2\ell_1+1)(2\ell_2+1)}{4\pi}
 +
{}_{(-s -s')}{\cal W}_{\ell_1} B_{\ell_2 S} B_{\ell_2 S'}^* C^{TT}_{\ell_2}
 +
  \left(
 
     \begin{array}{ccc}
 
     \begin{array}{ccc}
         \! l_1\! & l_2\!  & L\!  \\
+
         \! \ell_1\! & \ell_2\!  & L\!  \\
 
         \! s\! & -s\!  & 0\!
 
         \! s\! & -s\!  & 0\!
 
       \end{array}
 
       \end{array}
 
     \right) \left(
 
     \right) \left(
 
     \begin{array}{ccc}
 
     \begin{array}{ccc}
         \! l_1\! & l_2\!  & L\!  \\
+
         \! \ell_1\! & \ell_2\!  & L\!  \\
 
         \! s'\! & -s'\!  & 0\!
 
         \! s'\! & -s'\!  & 0\!
 
       \end{array}
 
       \end{array}
     \right)
+
     \right),
 +
\label{eq:pseudo_Cl_matrix}
 +
</math></p><p>
  
\end{gathered}</math>
+
where <math>{}_{(s s')}{\cal W}_{L} = \frac{1}{2L+1} \sum_{M} {}_{S} w_{LM} {}_{S'} w_{LM}^*</math> is a cross-power spectrum of scan history objects. Note that the <math>w(\hat{n}, s )</math> that we have used here can also incorporate a position-dependent weighting to optimize the pseudo-<i>C</i><sub>&#8467;</sub> estimate; this might be inverse-noise or a mask&mdash;the equations are unchanged. Writing the pseudo-<i>C</i><sub>&#8467;</sub> in position space (a la Dvorkin and Smith, 2009{{BibCite|dvorkin2008}}) with Wigner-d matrices we have
  
where <math>{}_{(s s')}{\cal W}_{L} = \frac{1}{2L+1} \sum_{M} {}_{S} w_{LM} {}_{S'} w_{LM}^*</math> is a cross-power spectrum of scan history objects. Note that the w(n,s) which we have used here can also incorporate a position dependent weighting to optimize the pseudo-Cl estimate, such as inverse-noise or a mask– the equations are unchanged. Writing the pseudo-Cl in position space (a la <font color=red> Dvorkin and Smith (2009)</font>) with Wigner-d matrices we have
+
</p><p><math>
 
 
<math>\begin{gathered}
 
 
\tilde{C}_{L}^{TT} = \frac{1}{8\pi} \sum_{S S'} \int_{-1}^{1} dz \ d^{L}_{00}(z)
 
\tilde{C}_{L}^{TT} = \frac{1}{8\pi} \sum_{S S'} \int_{-1}^{1} dz \ d^{L}_{00}(z)
\\ \times
+
\left[\sum_{\ell_1} d^{\ell_1}_{-s -s'}(z) {}_{(-s -s')}{\cal W}_{\ell_1} (2\ell_1+1) \right]  
\left[\sum_{l_1} d^{l_1}_{-s -s'}(z) {}_{(-s -s')}{\cal W}_{l_1} (2l_1+1) \right]  
+
\left[\sum_{\ell_2} d^{\ell_2}_{s s'}(z) B_{\ell_2 S} B_{\ell_2 S'}^* C^{TT}_{\ell_2}(2\ell_2+1) \right].
\\ \times
+
\label{eq:pseudo_Cl_corr}
\left[ \sum_{l_2} d^{l_2}_{s s'}(z) B_{l_2 S} B_{l_2 S'}^* C^{TT}_{l_2}(2l_2+1) \right].\end{gathered}</math>
+
</math></p><p>
  
This integral can be implemented exactly using Gauss-Legendre quadrature, with a cost of $\cal 0(l_{\rm max}^2 s_{\rm max}^2)$. For simplicity, we’ve written all the equations here for the auto-spectrum of a single detector, but the generalization to a map made by adding several detectors with different weighting is straightforward. The cost to compute all of the necessary terms exactly in that case becomes <math>\cal 0(l_{\rm max}^2 s_{\rm max}^2 N_{\rm det}^2)</math>.
+
This integral can be implemented exactly using Gauss-Legendre quadrature, with a cost of order (&#8467;<sub>max</sub><i>s</i><sub>max</sub>)<sup>2</sup>. 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 weighting is straightforward. The cost to compute all of the necessary terms exactly in that case becomes of order (&#8467;<sub>max</sub><i>s</i><sub>max</sub><i>N</i><sub>det</sub>)<sup>2</sup>.
  
Are beams really so difficult? On the flat-sky beam convolution is easy: just 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). Does the sphere really require everything to be so complicated? For a scan strategy which is fairly smooth across the sky, we can pretend that we are observing many independent flat-sky patches at high-L with fairly good accuracy. There is in fact a fairly good approximation to the beam convolved pseudo-Cl power spectrum which is essentially a flat-sky approximation. In the limit that <math>L \gg l_1</math>, with <math>C_{l_2}</math> and <math>B_{l_2}</math> being slowly-varying function in <math>l_2</math> the pseudo-Cl sum above can be approximated as <math>{\tilde{C}}_L^{TT} = C_L^{TT} \sum_{M} \left< \left| w(\hat{n}_p, M) \right|^2 \right>_p |B_{L M}|^2,</math> where the average <math><>_p</math> is taken over the full sky. It’s illustrative to consider three limits of this equation: for a “raster” scan strategy in which each pixel is observed with the same direction, we have <math>\left< \left| w(\hat{n}, M) \right|^2 \right>_p = 1,</math> and the predicted pseudo-Cl is just the power spectrum of the beam. For a &quot;best-case&quot; scan strategy, in which each pixel is observed many times with many different orientation angles, we have &lt; | w(, M) |<sup>2</sup> &gt;<sub>p</sub> = <sub>M0</sub>, and the transfer function is just 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 fsky factor, as expected. It just neglects the coupling between L multipoles (which can be calculated with the more complete equations above).
+
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 (encapsulated in the scan
 +
history objects described above).
  
==== Effective beam window functions ====
+
A scan strategy that is fairly smooth across the sky is nearly
The effective beam window functions $B(l)$ for HFI, computed using Quickbeam, are available in the [[The RIMO|RIMO]].
+
equivalent to  observing many independent flat-sky patches at high
They do not contain the pixel window function.
+
<i>L</i>. There is a fairly good approximation to the beam-convolved
 +
pseudo-power spectrum, which is essentially a flat-sky approximation. In the
 +
limit that <i>L</i>&#8811;&#8467;<sub>1</sub>, with <i>C</i><sub>&#8467;<sub>2</sub></sub> and <i>B</i><sub>&#8467;<sub>2</sub></sub> slowly-varying
 +
functions in &#8467;<sub>2</sub>, and using the equality
  
=== Pixelization Artifacts ===
+
</p><p><math>
For more details, see <cite>#planck2013-p03c</cite>
+
\sum_{\ell_2} (2 \ell_2+1)
 +
\left(
 +
    \begin{array}{ccc}
 +
        \! \ell_1\! & \ell_2\!  & L\!  \\
 +
        \! s\! & -s\!  & 0\!
 +
      \end{array}
 +
    \right)
 +
\left(
 +
    \begin{array}{ccc}
 +
        \! \ell_1\! & \ell_2\!  & L\!  \\
 +
        \! s'\! & -s'\!  & 0\!
 +
      \end{array}
 +
    \right)
 +
= \delta_{ss'},
 +
</math></p><p>
  
<font color=red>
+
the pseudo-<i>C<sub>&#8467;</sub></i> sum above can be approximated as
* Several codes available to simulate effects of pixelization.
+
* Mixes the CMB gradient into a pixelization ``noise'' with a level comparable to that of $2\mu Karcmin$ instrumental noise.
+
<math>
* Quantitative estimate of effect should be included with each released map, but expect not to matter significantly for CMB analysis, as small compared to instrumental noise.
+
{\tilde{C}}_L^{TT} = C_L^{TT} \sum_{M} \left< \left| w(\hat{n}_p, M) \right|^2 \right>_p |B_{L M}|^2,
</font>
+
\label{eq:pseudo_Cl_approx}
[sec:pixelization]
+
</math>
 +
where the average <math>\langle \rangle_p</math> 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,
 +
 +
</p><p><math>
 +
\left< \left| w(\hat{n}, M) \right|^2 \right>_p = 1,
 +
</math></p><p>
 +
 +
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,
 +
 +
\begin{equation}
 +
\left< \left| w(\hat{n}, M) \right|^2 \right>_p = \delta_{M0},
 +
\end{equation}
 +
 +
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 <i>f</i><sub>sky</sub> factor, as
 +
expected, but neglects the coupling between <i>L</i> multipoles (which
 +
can be calculated with the more complete equations above).
  
Planck maps are produced at resolution 11 <math>(N_{\rm side} = 2048)</math>, corresponding to pixels with a typical dimension of <math>1.7'</math>, comparable to the spacing between scanning rings . This results in an uneven distribution of hits within pixels, which introduces some complications in the analysis and interpretation of the maps. A sample of the hit distribution is illustrated in Fig. [fig:pixcoverage]. Below we discuss the simulation and modeling of this pixelization effect in more detail.
+
=== Pixelization artefacts ===
 +
We now discuss pixelization; for more details, see {{PlanckPapers|planck2013-p03c}}.
  
 +
Planck-HFI maps are produced at HEALPix resolution 11 (<i>N</i><sub>side</sub>=2048), corresponding to pixels with a typical dimension of 1.7arcmin. With the resolution comparable to the spacing between scanning rings, 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 the figure below.
  
[[Image:pixcoverage.png|frame|none|alt=image]]
+
[[Image:pixcoverage.png|thumb|500px|center|'''Figure 3:''' Illustration of TOI samples near the Galactic plane (grey
 +
  dots), over-plotted on a simulated CMB realization which has been
 +
  convolved with a 7' Gaussian FWHM beam and pixelized at <i>N</i><sub>side</sub>=2048. Associated scanning rings (grey lines) as well as centres of mass for the hit distribution (black arrows) are also plotted.]]
  
[fig:pixcoverage]
+
The collaboration has produced three codes that may be used to simulate the effect of pixelization on the observed sky: LevelS/TotalConvoler/Conviqt{{BibCite|wandelt2001}}{{BibCite|prezeau2010}}{{BibCite|reinecke2006}}; FeBeCoP{{BibCite|mitra2010}}; and [[#FICSBell|FICSBell]].
  
The collaboration has produced 3 codes which may be used to simulate the effect of pixelization on the observed sky, LevelS/TotalConvoler/Conviqt, FeBeCoP, and FICSBell <font color=red> references and further discussion of the three methods
+
For the measurement of CMB fluctuations, the effects of pixelization
and how they each simulate the pixelization effect.</font>.
+
may be studied analytically. On the small scales relevant to
 +
pixelization, the observed CMB is smooth, both due to physical damping
 +
as well as 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
 +
<math>\langle |\nabla T |^2 \rangle = \frac{1}{4\pi} \sum_{\ell} \ell(\ell+1)(2\ell+1) C_\ell^{TT} W_\ell \approx 1\times10^9 \mu K^2 / {\rm rad}^2,</math>
 +
where the approximate value is calculated for a
 +
&Lambda;CDM cosmology with a 7' FWHM Gaussian beam. The typical curvature of the observed temperature, on the other hand is given by <math>\langle |\nabla^2 T |^2 \rangle = \frac{1}{4\pi} \sum_{\ell} [\ell(\ell+1)]^2(2\ell+1) C_\ell^{TT} W_\ell \approx 7\times10^{14} \mu K^2 / {\rm rad}^4.</math> On the scales relevant to the maximum displacement from the centre of a 1.7' pixel, the maximum displacement is of order
 +
1' (= 3&times;10<sup>-4</sup>rad), and so the gradient term tends to dominate, although the curvature term is still non-negligible. For each observation of a pixel, we can denote the displacement from the pixel centre as <i>d</i>=<i>d</i><sub>&theta;</sub>+i<i>d</i><sub>&phi;</sub>. The average over all hits within a pixel gives an overall deflection vector, which we will denote for a pixel centre located at <math>\hat{n}</math> as <math>d(\hat{n})</math>. This represents the centre of mass of the hit distribution; [[#pixcoverage|Fig. 2]] shows these average deflections using black arrows. The deflection field <math>d(\hat{n})</math> may be decomposed into spin-1 spherical harmonics as <math>d_{\ell m} = \int_{4\pi} {}_1 Y_{\ell m}^* d(\hat{n}).</math> With a second-order Taylor expansion of the CMB temperature about each pixel centre, it is then possible to calculate the average pseudo-<i>C</i><sub>&#8467;</sub> power spectrum of the pixelized sky. This is given by
  
For the measurement of CMB fluctuations, it is also possible to gain intuition for the effects of pixelization analytically. On the small scales relevant to pixelization, the observed CMB is smooth, both due to physical damping as well as the convolution of the instrumental beam. Taylor expanding the CMB temperature about a pixel center to second order, the typical gradient amplitude is given by <math>\langle |\nabla T |^2 \rangle = \frac{1}{4\pi} \sum_{l} l(l+1)(2l+1) C_l^{T} W_l \approx 1\times10^9 \mu K^2 / {\rm rad}^2.</math> where the approximate value is calculated for a <math>\Lambda CDM</math> cosmology with a <math>7'</math>fwhm Gaussian beam. The typical curvature of the observed temperature, on the other hand is given by <math>\langle |\nabla^2 T |^2 \rangle = \frac{1}{4\pi} \sum_{l} [l(l+1)]^2(2l+1) C_l^{T} W_l \approx 7\times10^{14} \mu K^2 / {\rm rad}^4.</math> On the scales relevant to the maximum displacement from the center of a <math>1.7'</math> pixel, the maximum displacement is , and so the gradient term tends to dominate, although the curvature term is still non-negligible. For each observation of a pixel, we can denote the displacement from the pixel center as <math>d = d_{\theta} + i d_{\phi}</math>. The average over all hits within a pixel gives an overall deflection vector which we will denote for a pixel center located at <math>\hat{n}</math> as <math>d(\hat{n})</math>. This represents the center of mass of the hit distribution; in Fig. [fig:pixcoverage] we have plotted these average deflections using black arrows. The deflection field <math>d(\hat{n})</math> may be decomposed into spin-1 spherical harmonics as <math>d_{lm} = \int_{4\pi} {}_1 Y_{lm}^* d(\hat{n}).</math> With a second order Taylor expansion of the CMB temperature about each pixel center, it is then possible to calculate the average pseudo-Cl power spectrum of the pixelized sky. This is given by
+
<math>
 
+
C_\ell^{TT} = [1-\ell(\ell+1)R^d] {C}_\ell^{TT} W_\ell +
<math>\begin{gathered}
+
\frac{1}{2} \sum_{\ell_1 \ell_2} \frac{\ell_1(\ell_1+1)(2\ell_1+1)(2\ell_2+1)}{4\pi}  
C_l^{T} = [1-l(l+1)R^d] {C}_l^{T} W_l + \\
+
\left(
\frac{1}{2} \sum_{l_1 l_2} \frac{l_1(l_1+1)(2l_1+1)(2l_2+1)}{4\pi} \\
 
\times \left(
 
 
     \begin{array}{ccc}
 
     \begin{array}{ccc}
         \! l_1\! & l_2\!  & l\!  \\
+
         \! \ell_1\! & \ell_2\!  & \ell\!  \\
         \! l\! & -l\!  & 0\!
+
         \! \ell\! & -\ell\!  & 0\!
 
       \end{array}
 
       \end{array}
     \right)^2 C_{l_1}^{T} W_{l_1} \left[ C_{l_2}^{d+} + (-1)^{l + l_1 + l_2} C_{l_2}^{ d-} \right],
+
     \right)^2 C_{\ell_1}^{TT} W_{\ell_1} \left[ C_{\ell_2}^{d+} + (-1)^{\ell + \ell_1 + \ell_2} C_{\ell_2}^{ d-} \right],
\label{eqn:clt_pixelized}\end{gathered}</math>
+
\label{eqn:clt_pixelized}
 +
</math>
  
where <math>R^{d} = \langle |d|^2 \rangle/2</math> is half the mean-squared deflection magnitude (averaged over hits within a pixel, as well as over pixels). <math>C_l^{d+}</math> is the sum of the gradient and curl power spectra of <math>d_{lm}</math>, and <math>C_l^{d-}</math> is the gradient spectrum minus the curl spectrum. The <math>R^{d}</math> term describes a smearing of the observed sky due to pixelization. For uniform pixel coverage of <math>N_{\rm side}=2048</math> pixels <math>\sqrt{ \langle |d|^2 \rangle } = 0.725'</math>. For the hit distribution of Planck frequency maps, <math>R^{d}</math> is typically within <font color=red> xxx. calculate for final
+
where <i>R<sup>d</sup></i>=&lang;|d|<sup>2</sup>&rang;/2 is half the mean-squared deflection magnitude (averaged over hits within a pixel, as well as over pixels). <math>C_\ell^{d+}</math> is the sum of the gradient and curl power spectra of <i>d<sub>&#8467;m</sub></i>, and <math>C_\ell^{d-}</math> is the gradient spectrum minus the curl spectrum. The <i>R<sup>d</sup></i> term describes a smearing of the observed sky due to pixelization. For uniform pixel coverage of <i>N</i><sub>side</sub>=2048 pixels &radic;&lang;|<i>d</i>|<sup>2</sup>&rang;=0.725', while, for the hit distribution of Planck frequency maps, <i>R<sup>d</sup></i> is typically 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&times;10<sup>-4</sup> for &#8467;&le;3000.
maps, looks like will be better than 10%</font>percent of this value, and so this term is accurately described by the pixel window function, which is derived under the assumption of uniform pixel coverage.
 
  
The effect of pixelization is essentially 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 the [<cite>#planck2013-p12</cite>|Planck gravitational lensing] paper, where the subpixel deflection field constitutes a potential source of bias for the measured lensing potential. Indeed, Eq. [eqn:clt<sub>p</sub>ixelized] is just a slightly modified version of the usual first order CMB lensing power spectrum (<font color=red>Hu (2000)</font>, <font color=red>Lewis and Challinor (2006)</font>) to accommodate curl modes.
+
The effect of pixelization is essentially 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 the Planck 2013 gravitational lensing paper ({{PlanckPapers|planck2013-p12}}), where the sub-pixel deflection field constitutes a potential source of bias for the measured lensing potential. Indeed, the above equation is just a slightly modified version of the usual first-order CMB lensing power spectrum (see Hu 2000{{BibCite|Hu2000}}; Lewis &amp; Challinor 2006{{BibCite|Lewis2006}}) to accommodate curl modes.
  
A useful approximation to Eq.  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 <math>d(\hat{n})</math> couples the CMB gradient into a source of noise with an effective level given by <math>\sigma^{N} \approx \sqrt{ R^T \frac{4\pi}{N_{\rm pix}} \langle | d(\hat{n}) |^2
+
A useful approximation to this equation, 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 <math>d(\hat{n})</math> couples the CMB gradient into a source of noise with an effective level given by <math>\sigma^{N} \approx \sqrt{ R^{TT} \frac{4\pi}{N_{\rm pix}} \langle | d(\hat{n}) |^2
 
\rangle }, % (\muKarcmin ),</math>
 
\rangle }, % (\muKarcmin ),</math>
  
where the average is taken over all pixels and <math>R^T</math> is half the mean-squared power in the CMB gradient: <math>R^{T} = \frac{1}{8\pi} \sum_{l} l(l+1)(2l+1) \tilde{C}_l^{T}.</math> For frequency-combined maps, <math>\sqrt{ \langle | d(\hat{n}) |^2 \rangle  }</math> is typically on the order of <math>0.1'</math>, and so the induced noise is at the level of <math>\sigma^{N} \sim 2 \mu K arcmin</math>. This is small compared to the instrumental contribution, although it does not disappear when taking cross-spectra, depending on how coherent the hit distributions of the two maps in the cross-spectrum are.
+
where the average is taken over all pixels and <i>R<sup>TT</sup></i> is half the mean-squared power in the CMB gradient: <math>R^{TT} = \frac{1}{8\pi} \sum_{\ell} \ell(\ell+1)(2\ell+1) \tilde{C}_\ell^{TT}.</math> For frequency-combined maps, <math>\sqrt{ \langle | d(\hat{n}) |^2 \rangle  }</math> is typically on the order of 0.1', and so the induced noise is at the level of &sigma;<sup>N</sup>&asymp;2&mu;K.arcmin. This is small compared to the instrumental contribution, although it does not disappear when taking cross-spectra, depending on how coherent the hit distributions of the two maps in the cross-spectrum are.
  
 
= References =
 
= References =
<biblio force=false>
+
<References />
#[[References]]  
+
</biblio>
+
 +
 
 +
[[Category:HFI data processing|003]]

Latest revision as of 15:07, 2 July 2018

Introduction[edit]

The reconstruction of the spatial response for the Planck High Frequency Instrument (HFI) is informed by observation of planets. Furthermore, observations of Uranus and Neptune are used to obtain an absolute calibration for the 545 and 857 GHz channels.

The spatial response, also known as the beam response or the point spread function, relates the absolute calibration to different angular scales. We follow the nomenclature of Planck Collaboration VII Planck-2013-VII[1] (2014), where the "scanning beam" is defined as the coupled response of the optical system, the deconvolved time response function, and the software low-pass filter, applied to the data. The "effective beam" represents the averaging of the signal due to the scanning of the telescope and the mapmaking procedure, and varies from pixel to pixel across the sky.

The 2015 and 2018 Data Releases[edit]

Most of the content found on this page, as well as further information regarding the beam products used for the 2015 (PR2) release of Planck data, can be found in the HFI DPC Paper 1 Planck-2015-A07[2].

For the PR3 2018 data release, the scanning beams are the same as for the previous release. The effective beams (Febecop) are very slightly affected by the rejection of the data taken during the last 1000 stable pointing periods of the cold mission (see Planck-2020-A3[3]). This affects a small fraction of the sky in a negligible way.

The main improvement in PR3 over previous releases concerns the computation of effective beam window functions relevant for cosmological analysis. As explained in Beam_Window_Functions, beam matrices describing exactly the leakage from temperature to polarization and cross-talk between spectra induced by the beam mismatch and beam ellipticity have been computed with QuickPol for all relevant combination of detectors and data splits. Accounting for such systematic effects in this way greatly improves the inter-frequency consistency of the CMB TE power spectra, compared to the 2015 analysis (see Planck-2020-A5[4]).


Planet and main beam description[edit]

Here we redefine the "main beam" to be the scanning beam out to 100 arcmin from the beam axis. The sidelobe structure at this radius is dominated by diffraction at the mirror edges and falls as θ-3, where θ is the angle from the main beam axis. The main beam is used to compute the effective beam and the effective beam window function, which describes the filtering of sky signals. The smearing of the main beam cannot be significantly reduced without boosting the high-frequency noise. The regularization function (a low-pass filter) chosen has approximately the same width as the instrumental transfer function. The deconvolution significantly reduces the long tail of the scanning beam (see Fig. 1). The deconvolution also produces a more symmetric time response, so residual “streaking” appears both ahead and behind the main beam, although ahead of the beam it is at a level of less than 10−4 of the peak response. The far sidelobes are defined as the response for θ > 5°, roughly the minimum in the optical response. The response begins to rise as a function of angle θ beyond this due to spillover. The far sidelobes are handled separately from the beam effects (see Paper 2 Planck-2015-A08[5] for details).

Figure 1: Scan across Jupiter with bolometer 143-1b, with pre-deconvolution in black and post-deconvolution in red. The plotted y-axis is linear within the range ± 10-6 and logarithmic elsewhere. Similarly, the x-axis range is linear within the range ± 0.16s, and logarithmic outside.

The observations of planets are also used to estimate their brightness as a function of frequency, see Planck-2017-LII[6].

Key Changes between 2013 and 2015[edit]

There are several key changes in the reconstruction of the main beam since the 2013 data release, including those listed below.

  • The time-ordered data from Saturn and Jupiter observations are merged prior to B-spline decomposition, accounting for residual pointing errors and variable seasonal brightness. This is achieved by determining a scaling factor and a pointing offset by fitting the time-ordered data to a template from a previous estimate of the scanning beams. We iterate the planet data treatment updating the template with the reconstructed scanning beam. The process converges in five iterations within 0.1% accuracy in the effective beam window function.
  • Steep gradients in the signal close to the planet reduce the completeness of the standard glitch-detection and subtraction procedure; therefore the planet timelines are deglitched a second time.
  • The beam pipeline destripes the planet data, estimating a single baseline between 3° and 5° before the peak for each scanning circle. Baseline values are smoothed with a 40 circle sliding window. The entire scanning circle is removed from the beam reconstruction if the statistic in the timeline region used to estimate the baseline is far from Gaussian.
  • The main beam is now recontructed on a square grid that extends to a radius of 100 arcmin from the centroid, as opposed to 40 arcmin in the 2013 data release. The cutoff of 100 arcmin was chosen so that a diffraction model of the beam at large angles from the centroid predicts that less than 5×10-5 of the total solid angle is missing.
  • The scanning beam is constructed by combining data from Saturn observations, Jupiter observations, and physical optics models using GRASP software.
  • No apodization is applied to the scanning beam map.

The update of the time response deconvolved data has slightly changed the scanning beam, the effective beam solid angles, and the effective beam window functions.

Table 1. Band-average scanning beam solid angle ([math]\Omega _\mathrm{SB} [/math]) and Monte Carlo-derived errors ([math]\Delta\Omega _\mathrm{MC}[/math]) including noise, residual glitches, and pointing uncertainty.
Band [math]\Omega _\mathrm{SB}[/math] [math]\Delta \Omega _\mathrm{MC}[/math]
[GHz] [arcmin[math]^2[/math]]
100 104.62 0.13 %
143 58.50 0.07 %
217 26.92 0.13 %
353 25.93 0.09 %
545 25.23 0.08 %
857 23.04 0.08 %

A portion of the near sidelobes was not accounted for in the effective beam window function of the 2013 data (see Planck-2013-VI[7]). To remedy this, the domain of the main beam reconstruction is extended to 100 arcmin, with no apodization. Saturn data are used where they are signal-dominated. Where the signal-to-noise ratio of the Saturn data fall below 9, azimuthally-binned Jupiter data are used. At larger angles, below the noise floor of the Jupiter data, we use a power law (∝ θ-3), whose exponent is derived from GRASP to extend the beam model to 100 arcmin. Figure Fig. 2 shows a scanning beam map with a diagram of the regions handled differently in the hybrid beam model. A summary of the solid angles of the hybrid beams is shown in Table 1. Figure 3 shows a contour plot of all the scanning beams referenced to the centre of the focal plane.

Figure 2: Scanning beam map for detector 143-6, indicating approximately the regions that are handled with a different data selection or binning.
Figure 3: B-spline hybrid scanning beams reconstructed from Mars, Saturn, and Jupiter. The beams are plotted in logarithmic contours of −3, −10, −20, and −30 dB from the peak. PSB pairs are indicated with the "a" bolometer in black and the "b" bolometer in blue.

Beam Error Budget[edit]

As in the 2013 release, the beam error budget is based on an eigenmode decomposition of the scatter in simulated planet observations. A reconstruction bias is estimated from the ensemble average of the simulations. We generate 100 simulations for each planet observation, including pointing uncertainty, cosmic ray glitches, and the measured noise spectrum. Simulated glitches are injected into the timeline with the correct energy spectrum and rate (Planck-2013-X[8]) and are detected and removed using the deglitch algorithm. Noise realizations are derived as described in Planck-2013-VII[1]. Pointing uncertainty is simulated by randomizing the pointing by 1.5 arcsec rms in each direction. The improved signal-to-noise ratio compared to 2013 leads to smaller error bars: for instance, at ℓ = 1000 the uncertainties on B are now (2.2, 0.84, 0.81)×10-4 for 100, 143, and 217 GHz frequency-averaged maps respectively, reduced from the previous uncertainties of (61, 23, 20)×10-4.

A “Reduced Instrument Model” (RIMO; see The_RIMO) containing the effective B for temperature and polarization detector sets (for auto- and cross-spectra at 100 to 217 GHz) is included in the release, for a sky fraction of 100 %; another RIMO is provided for a sky fraction of 75 %. They both contain the first five beam error eigenmodes and their covariance matrix, for the multipole ranges [0, ℓmax] with ℓmax = 2000, 3000, and 3000 at 100, 143, and 217 GHz, respectively (instead of 2500, 3000, and 4000 previously). These new ranges bracket more closely the ones expected to be used in the likelihood analyses, and ensure a better determination of the leading modes on the customized ranges. As described in appendix A7 of Planck Collaboration XV (2014) Planck-2013-XV[9], these beam window function uncertainty eigenmodes are used to build the C covariance matrix used in the high-ℓ angular power spectrum likelihood analysis. It was found that the beam errors are negligible compared to the other sources of uncertainty and have no noticeable impact on the values or associated errors of the cosmological parameters.

Consistency of Beam Reconstruction[edit]

To evaluate the accuracy and consistency of the beam reconstruction method, we have compared beams reconstructed using Mars for the main beam part instead of Saturn, and using data from Year 1 or Year 2 only. We compared the window functions obtained with these new beams to the reference ones; new Monte Carlo simulations were created to evaluate the corresponding error bars, and the results are shown in Figs. 4 and 5. Note that the effective beam window functions shown in these figures are not exactly weighted for the scan strategy, rather we plot them in the raster scan limit, where B is defined as the sum over m of the Bm. The eigenmodes were evaluated for these different data sets, allowing us to compare them with the reference beam through a χ2 analysis. The discrepancy between each data set and the reference beam is fitted using Ndof=5 eigenmodes using ℓmax = 1200, 2000, and 2500 for 100, 143, and 217GHz channels, respectively. The χ2 is then defined as

[math] \chi^2 = \sum_{i=1}^{N_\mathrm{dof}} (c_i / \lambda_i)^2, \label{eq:chi2} [/math]

where ci is the fit coefficient for eigenvector i, with eigenvalue λi. For each data set, the p-value to exceed this χ2 value for a χ2 distribution with Ndof degrees of freedom is indicated in Table 2. A high p-value indicates that the given data set is consistent with the reference beam within the simulation-determined error eigenvectors, with 100 % indicating perfect agreement. We find excellent agreement between the yearly and nominal beams for the SWB bolometer channels; however, in order to obtain reasonable beam agreement for the PSB detector sets, we find that the yearly beam errors must be scaled by a factor of 4. We conclude that there is an unknown systematic error in the beam reconstruction that is not accounted in the simulations used to estimate the beam accounting in the Monte Carlo error bars. We assign a scaling factor of 4 to the error eigenmodes to account for this systematic error. This additional error scaling has a negligible effect on the cosmological parameters.

Figure 4: Comparison of Year 1 and Year 2 based beams with the reference beam (Full Mission). Window functions are calculated using Quickbeam in raster scan configuration. Error bars are computed using MC simulations for Year 1 and Year 2, taking the maximal one for each.
Figure 5: Comparison of Mars-based beams with the reference beam (Full Mission). Mars-based beams are constructed from Mars data for the main beam (below about 10 arcmin) and Saturn and Jupiter data for the larger scales. Window functions are calculated using Quickbeam in raster scan configuration. Error bars are computed using MC simulations for Mars.
Table 2. P-values in percent for the χ2 comparison of the nominal beam with Year 1 and Year 2 data sets, following the definition in the above equation.
Band Year 1 Year 2
100-ds1 66.95 92.44
100-ds2 85.12 75.99
143-ds1 72.24 99.83
143-ds2 1.14 16.81
143-5 94.88 97.57
143-6 96.14 98.78
143-7 93.78 95.67
217-ds1 78.12 69.74
217-ds2 27.33 30.87
217-4 95.99 68.44
217-1 99.08 97.49
217-2 97.18 98.86
217-3 97.59 99.07

Colour Correction of the Beam Shape[edit]

In general, the measured beam shape is a function of the spectral energy distribution (SED) of the measurement source. This is of particular concern because we measure the beam on a source with a roughly Rayleigh-Jeans SED, yet we use the effective beam window function to correct the CMB.

Planck Collaboration VII (2014) Planck-2013-VII[1] described possible levels of this bias using GRASP calculations with the pre-launch telescope model (Maffei et al. 2010 [10]; Tauber et al. 2010 [11]). The prelaunch calculations did not agree well enough with the data to allow a direct application of the colour correction. A new telescope model based on flight data, and presented in a forthcoming paper, predicts beams that agree better with the data at 100–217 GHz, but show worse agreement at 353 GHz. The new model predicts a different beam shape colour correction, though at a similar order of magnitude to that shown in Planck Collaboration VII (2014) Planck-2013-VII[1]. This work is ongoing and will be completed after the 2015 release.

We search for a possible signature of the colour correction effect in the data. First, we check for consistency in the CMB angular power spectra derived from different detectors and frequency bands using the SMICA algorithm, in order to find discrepancies in beam shape and relative calibration of different data subsets of the CMB (Planck Collaboration A13 2014 Planck-2015-A11[12]). There are hints of differences that are orthogonal to the beam error eigenmodes and appear as changes in relative calibration, but the preferred changes are at the 0.1 % level. Because this correction is of the same order of magnitude as the relative calibration uncertainty between detectors, we treat it as insignificant.

Second, we look at the relative calibration between detectors within a band for compact sources as a function of the source SED. In this method, any discrepancies due to solid angle variation with SED are degenerate with an error in the underlying bandpass. Additionally, because of the lack of bright, compact sources with red spectra, this method is limited in its signal-to-noise ratio. In the limit that any detected discrepancy is completely due to solid-angle variation with colour, at the level of 1 % in solid angle we do not detect variations consistent with the physical optics predictions. Given the lack of measurement of this effect at the current levels of uncertainty, as well as uncertainties in modelling of the telescope, we note that this effect may be present at a small level in the data, but do not attempt to correct for it or include it in the error budget.

Effective beam window functions[edit]

Changes in QuickBeam[edit]

The effective beam window functions B(ℓ) for HFI, were computed using QuickBeam as in 2013, with a few improvements.

Higher spin[edit]

Since the scanning 2015 beams are slightly more elongated than the 2013 ones, because of different transfer function deconvolution, the QuickBeam approach was extended to include the beam coefficients bℓs up to |s| = 6 , instead of |s| ≤ 4. We checked that, for ℓ ≤ 3000 (or ℓ ≤ 2000 at 100GHz), the coefficients with |s| > 6 account for much less of 0.1% of the beam solid angle.

Impact of sky cut[edit]

The effective beam window function will be affected by the finite area of the pixels on which the sky signal is averaged to produce the map. While the HEALPix pixels all have the same surface area for a given resolution parameter Nside, their shape differ and the trace of their moment of inertia:

[math] {\rm Tr} ({\bf I}_p) = \int d{\bf x} d{\bf y} \left[ ({\bf x}-{\bf x}_p)^2 + ({\bf y}-{\bf y}_p)^2 \right], [/math]

where (x, y) are the tangential plane coordinates in the pixel p (valid for typical pixel size of 1.7' for Nside=2048) and (xp,yp) is the pixel centre, varies across the sky, especially in the polar regions, as illustrated in Fig. 6 below.

Figure 6: Map of the variation of the trace of the HEALPix pixel moment of inertia, for Nside=2048.

The full-sky average of this moment of inertia determines the nominal HEALPix pixel window function, which is not included in the effective beam window function provided and must be corrected for separately. When the Galactic masks used in Planck blank out the equatorial region, the window function of the remaining pixels will differ from the full-sky average, and will affect the resulting beam window functions. We followed a semi-analytical approach in QuickBeam, where the beam window functions are corrected from the departure between the full-sky average of the tensor of inertia and its average over un-masked pixels, and checked that the results agree very well with those of Monte Carlo simulations of the exact Planck scanning of the pixels performed with FEBeCoP, for several Galactic masks with fsky=100%, 80%, and 40%.

The effective beam window function are therefore provided for two different sky fractions (100% and 75%) in different RIMOs. For these two sky fractions the square of the beam window functions (relevant for the power spectra analysis) are identical at ℓ=1 and have a relative difference of about 0.4% at ℓ=2000.

Data Products[edit]

Beam products are available from the instrument model (The_RIMO) and further described in Scanning_Beams and Effective_Beams.

Effective beam window functions are available in specific RIMOs. They do not' contain the nominal HEALPix pixel window function.

The 2013 Data Release[edit]

Scanning beams[edit]

The scanning beams describe the instrument’s instantaneous beam profile. Due to the near constant spin rate of the spacecraft, time domain effects (including residual time response and lowpass filtering) are degenerate with the spatial response due to the optical system. The scanning beam reconstruction recovers both of these effects, except for residual time-domain effects on a longer timescale than can be captured with the extent of the scanning beam model.

In Planck-2013-VII[1] we consider two models of the beam in order to better understand systematic effects in the reconstruction. Here we describe only the "B-spline beams," which are used to compute the delivered effective beam (see next section).

B-spline Beam construction[edit]

We use seasons 1 and 2 of the Mars observations to reconstruct the beam. The data are processed with the "bigPlanets" TOI processing (see here or section 3.11 of Planck-2013-VI[7] for more information). We use JPL Horizons ephemerides to determine the pointing of each detector relative to the planet. We then subtract the astrophysical background in the time domain using a bicubic interpolation of the Planck maps.

The time-ordered data are used to fit a two dimensional B-spline surface using a least squares minimization and a smoothing criterion to minimize the effects of high spatial frequency variations. We therefore assume the scanning beam to be smooth. The smoothing criterion, as well as the locations of the nodes used to compute the B-spline basis functions, are set using GRASP physical optics simulations as inputs, which are the best assumptions on the spatial frequency content of the in-flight beams.

The smoothing criterion is defined as follows:

[math]\eta = \displaystyle{\sum_{i=1}^{g}\left(b^{k}(\lambda_{i+})-b^{k}(\lambda_{i-})\right)^2}, \label{smoothcrit}[/math]

where η is the Smoothing Criterion and bk represents the kth beam derivative evaluated at the node locations.

The global inversion criterion is:

[math]\zeta = \eta + p\times \delta[/math],

with δ being the usual least squares estimator and the p coefficient giving the relative weight of δ with respect to the smoothing criterion

[math]\delta = \displaystyle{\sum_{r=1}^{m}}\left(y_{r} - b(x_{r})\right)^2\label{estimator}[/math].

Here δ is the usual least squares criterion, r is the index relative to the m data points (r in the range 1 to m), yr is the planet data of sample r, xr is the pointing of sample r, and b is the reconstructed beam.

The B-spline nodes are located on a regularly spaced grid in the detector coordinate frameset. At the edge of the reconstructed beam map area, four coincident nodes are added to avoid vanishing basis functions.

Let Bi, k+1 be a k-degree B-spline built using nodes {λi, …, λi+k+1} (the De Boor-Cox method), then

[math]B_{i,1}(x) = \left\{ \begin{array}{l} 1, \mbox{ if } x \in \mbox{[} \lambda_{i}, \lambda_{i+1} \mbox{]},\\ 0, \mbox{ if } x \notin \mbox{[} \lambda_{i}, \lambda_{i+1} \mbox{]}, \end{array} \right.[/math]

[math]B_{i, l+1}(x) = \displaystyle{\frac{x - \lambda_{i}}{\lambda_{i+l} - \lambda_{i}}} B_{i,l}(x) + \displaystyle{\frac{\lambda_{i+l+1}-x}{\lambda_{i+l+1}-\lambda_{i+1}}} B_{i+1, l}(x)[/math],

[math]l=1, \ldots, k[/math].


Figure 1: Focal plane plot of B-spline scanning beams using in-flight pointing reconstruction. The contours are -3, -10, -20, and -30 dB from the peak, and for PSB pairs the "a" bolometer is plotted in black and the "b" in blue.

Simulations and errors[edit]

We estimate the reconstruction bias and noise in the measurements using an ensemble of simulated planet observations for each channel. Further details are discussed in Planck-2013-VII[1]. Kept fixed in each simulation are:

  • the input beam assumed, specifically we use a over-sampled version of the reconstructed B-Spline beam (or whatever comes out of the current ongoing tests!);
  • the astrophyical background is the same as that subtracted from the real data;
  • StarTracker pointing (using the "ptcor6" pointing model).

The following are varied in each simulation:

  • detector noise realizations, obtained by filtering randomly -generated white noise with the measured noise PSDs;
  • random pointing errors with 2 arcsec rms, and a spectrum that replicates the real errors;
  • simulated glitches and the deglitching procedure;
  • Mars brightness temperature variability.

400 simulated timelines are generated for each bolometer and for each of the two seasons of Mars observations used in the beam reconstruction. The simulated timelines are made into beam maps, projecting onto the B-Spline basis in the same way as for the real data.

The beam maps are propagated to effective beam window functions using the "QuickBeam" approach (see effective beams below) and used to evaluate the reconstruction bias and to construct error eigenmodes in the effective beam window function.


Residuals[edit]

There are two known beam effects that are not included in the main beam model and are estimated as a separate bias in flux and angular power spectrum measurement: (1) long tails due to errors in low frequency time response deconvolution; and (2) near sidelobes.

We stack all five observations of Jupiter to estimate the long timescale residuals due to incomplete deconvolution of the long timescale response.


Effect of mirror support print-through[edit]

The Planck reflectors suffer from print-through of the honeycomb structures that support the carbon-fibre face sheets. While the size of the deformation has been measured during tests below room temperature to be less than 20μm, it is the strict periodicity of the deformation that contributes most to the additional beam-shape change. A simple grating equation has been shown to describe the angular positions of the resulting near side lobes very well:


[math] \displaystyle{ \sin \theta _n = \frac{n \lambda}{Yd} } \label{DimplingLobes}[/math],

where:

  • θn is the angular position of the nth-order lobe from the central beam peak;
  • λ is the wavelength of the radiation;
  • d is the grating spacing of the periodicity;
  • Y is a factor that describes the position of the each reflector along the optical path;
  • Y=1.00 for the primary reflector;
  • Y=1.80 for the secondary reflector.


Three possible periodicities (19.6mm, 30mm, 52mm) in the honeycomb array dominate the Planck dimpling pattern for the 857-GHz detectors, although 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 arcminute beam window, but they contribute at most (0.05 ± 0.01)% to the integrated beam.

A map of Jupiter has been created for each 857-GHz detector, using the first four surveys, with the background subtracted using the same sky area for the survey taken 12 months before, in which the planet is not present. The background-subtracted maps for each survey were stacked to make a single map for each of the four detectors and these were further stacked to create a single 857-GHz band map, using the standard detector weightings proportional to 1/(NET)2. All the maps have 96 pixels per side, covering ±1.4º, and are scaled to the RMS noise of the background. Central beam values range from 205,000 to 256,000, depending on the noise of the detector.

For each of these five maps of Jupiter an elliptical Gaussian main beam was fitted, and the amount of saturation estimated by finding the peak increase that gives the best Gaussian fit, typically 10-20%. A circular ruze envelope (see Planck-2013-VII[1]) was also fitted, together with a tilted Gaussian component slightly offset from the beam in the cross-scan (optical x-axis) direction, using data more than 0.13° from the beam, with the known positions of the dimpling lobes masked out. Once these strongest components were removed, a tilted elliptical Gaussian was fitted for each visible dimpling lobe in up to five different sets of lobes; however, the innermost lobe set (52-mm periodicity) is generally obscured by the ruze component and for the outermost set (19.6-mm periodicity) most lobes are too weak to be fitted for individual detectors. The contributing areas of the lobes were determined by three methods: from the fitted Gaussians; from summation within boxed areas, using offsets determined from surrounding boxed areas; and from summation in the boxed areas using an offset determined from Gaussian fitting.

The uncertainties in the beam component areas given here include the spread in values from these different methods. Similarly, the lobe peaks in decibels are calculated using the raw, fitted and estimated desaturated beam peak values, and the uncertainties reflect the variations produced by the three methods. Averaged over the five map fits, the dimpling lobes contribute (0.47 ± 0.13)% of the beam area, while the ruze envelope accounts for (9.2 ± 0.7)% of the beam, with the main beam making up the remaining (90.3 ± 0.7)%.

For lower frequencies the dimpling lobes are less visible and will not be discussed here. For example, the noise floor for the 545-GHz band map is about 47 dB below the observed Jupiter peak, so that no dimpling lobes are visible further out than the 52mm lobes at 0.61°. For 353 GHz the noise floor is at about 38 dB so that even the 52-mm lobes are not observed.

Figure 2: Stacked Jupiter map for the 857-GHz detectors. Left: the data with main beam and ruze envelope subtracted. Right: the fitted elliptical Gaussian dimpling lobes with GRASP contours plotted on top.


The figure above shows how the dimpling lobes seen for the 857GHz band Jupiter map correspond to the contours calculated by the physical optics GRASP package (produced by TICRA, Denmark). The GRASP simulation was performed for the 857-1 detector and assumed a uniform dimpling distortion of 10μm. The biggest departure from the simple grating model is that the lobe pattern is elongated in the cross-scan (vertical) direction—lobes that should lie on a line 30° from the vertical are consistently found around 25° and those on the 60° line cluster around 53°. This is due to the offset geometry of the mirror system, whereby the dimpling print-though is foreshortened along the optical x-axis due to the tilt of the mirrors, producing greater lobe spacing as seen by the incoming photons. The amplitudes of the fitted dimpling lobes vary significantly within each set, whereas the GRASP model shows a constant amplitude for lobes within a set (except for the 30-mm periodicity). Similar amplitudes are seen in lobes lying directly across the beam centre from each other, indicating that "facesheet dimpling" has not occurred uniformly in all directions, as the GRASP model assumes.

For the 30-mm periodicity the vertical lobes are significantly weaker than those at the sides, just as the GRASP model shows. Generally, the fitted lobes are somewhat stronger than those seen by GRASP, indicating that the dimpling is indeed larger than 10μm. While the dimpling lobes are measureable at the highest two frequencies, with a strong source (i.e., Jupiter), 90% of the dimpling lobe area at 857GHz, and 10% for all other frequencies, already appears in the beam functions described above and needs no further correction. However, the presence of the uniform dimpling lobes provides a useful window by which to investigate the mirror geometry.

Effective Beams[edit]

The effective beam is the average of all scanning beams pointing at a certain direction within a given pixel of the sky map for a given scan strategy. It takes into account the coupling between azimuthal asymmetry of the beam and the uneven distribution of scanning angles across the sky. It captures the complete information about the difference between the true and observed images of the sky. They are, by definition, the objects whose convolution with the true CMB sky produce the observed sky map.

Several methods of effective beam determination have been developed and cross-validated. The main products are produced using FEBeCoP and details of the processing are given in the Effective Beams products page. See also the equivalent page discussing the LFI beams


FEBeCoP[edit]


The full algebra for this method for the calculation of effective beams was presented in Ref.[13]. Here we summarize the main results. The observed temperature sky [math]\widetilde{\mathbf{T}} [/math] is a convolution of the true sky [math]\mathbf{T} [/math] and the effective beam [math]\mathbf{B}[/math]:

[math] \widetilde{\mathbf{T}} \ = \ \Delta\Omega \, \mathbf{B} \cdot \mathbf{T}, \label{eq:a0} [/math]

where

[math] B_{ij} \ = \ \left[ \sum_t A_{ti} \, b(\hat{\mathbf{r}}_j, \hat{\mathbf{p}}_t) \right] \bigg/ \left[{\sum_t A_{ti}} \right] \, . \label{eq:EBT2} [/math]

Here [math]t[/math] labels time samples, [math]A_{ti}[/math] is 1 if the pointing direction falls in pixel number i, otherwise it is 0, [math]\mathbf{p}_t[/math] represents the exact pointing direction (not approximated by the pixel centre location), and [math]\hat{\mathbf{r}}_j[/math] is the centre of pixel number j, where the scanning beam [math]b(\hat{\mathbf{r}}_j, \hat{\mathbf{p}}_t)[/math] is being evaluated if the pointing direction falls within the cut-off radius of about 2.5 × FWHM.

The algebra is a bit more involved for polarized detectors. The observed Stokes parameters at a pixel i, [math](\widetilde{I}, \widetilde{Q}, \widetilde{U})[/math]i, are related to the true Stokes parameters [math](I, Q, U)[/math]i, by the following relation:

[math] ( \widetilde{I} \quad \widetilde{Q} \quad \widetilde{U})_i^{\rm T} \ = \ \Delta\Omega \sum_j \mathbf{B}_{ij} \cdot (I \quad Q \quad U)_j^{\rm T}, \label{eq:a1} [/math]

where the polarized effective beam matrix

[math] \mathbf{B}_{ij} \ = \ \left[ \sum_t A_{tp} \mathbf{w}_t \mathbf{w}^{\rm T}_t \right]^{-1} \sum_t A_{ti} \, b(\hat{\mathbf{r}}_j, \hat{\mathbf{p}}_t) \, \mathbf{w}_t \mathbf{W}^{\rm T}(\hat{\mathbf{n}}_j,\hat{\mathbf{p}}_t) \, , \label{eq:a2} [/math]

and [math]\mathbf{w}_t[/math]and [math]\mathbf{W}(\hat{\mathbf{r}}_j, \hat{\mathbf{p}}_t)[/math] are the the polarization weight vectors, as defined in Ref.[13].

The task is to compute [math]B_{ij}[/math] for temperature-only beams and the 3×3 matrices [math]\mathbf{B}_{ij}[/math] for each pixel i, at every neighbouring pixel j that falls within the cut-off radius around the the centre of the ith pixel.

The effective beam is computed by stacking within a small field around each pixel of the HEALPix sky map. Due to the particular features of Planck's scanning strategy, coupled to the beam asymmetries in the focal plane, and data processing of the bolometer and radiometer TOIs, the resulting Planck effective beams vary over the sky.

FEBeCoP, given information on the Planck scanning beams and detector pointing during a mission period of interest, provides pixelized stamps of both the Effective Beam and the Point Spread Function, at all positions of the HEALPix-formatted map pixel centres.

FICSBell[edit]

Since the HFI beams are not azimuthally symmetric, the scanning strategy has to be taken into account in the effective beam response modelling (for more details, see Planck-2013-VII[1]). This is done using the "FICSBell" method (Hivon et al., in preparation). This generalizes to polarization and to include other sources of systematics the approach used for TT Cℓ estimation in WMAP-3 [14] and by Smith et al (2007) [15] in the detection of CMB lensing in WMAP maps. The different steps of the method used for this study can be summarized as follows.

  1. The scanning-related information (i.e., statistics of the orientation of each detector within each pixel) is computed first, and only once for a given observation campaign:

    [math] w^d_s({\bf r}_p) = \sum_j e^{i s \psi_j}, \label{eq:e1} [/math]

    where ψj is the orientation of the detector with respect to the local meridian during the measurement j occurring in the direction rp. 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 on disc.

  2. The scanning-beam map or beam model of each detector d is analysed into its spherical harmonic coefficients

    [math] b^d_{\ell s} = \int d{\bf r} B_d({\bf r}) Y_{\ell s}({\bf r})\label{scanningBlm} [/math]

    where Bd(r) is the beam map centred on the North pole, and Yℓs(r) is the spherical harmonics basis function. Higher s indices describe higher degrees of departure from azimuthal symmetry and, for HFI beams, the coefficients bdℓs are decreasing functions of s at most multipoles considered. 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)[16] reached a similar conclusion in their deconvolution of LFI beams. For these reasons, only modes with [math]|s| \le 4[/math] are considered in the present analysis.

  3. For a given CMB sky realization t, described by its spherical harmonic coefficients aℓm = ∫ dr t(r) Yℓm(r), the bdℓs coefficients computed above are used to generate s-spin weighted maps,

    [math] m^d_{s}({\bf r}) = \sum_{\ell m} b^d_{\ell s}\ a_{\ell m}\ {}_sY_{\ell m}({\bf r}), \label{eq:e3} [/math]

    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

    [math] m({\bf r}) = \left(\sum_d \sum_s w^d_s({\bf r}) m^d_s({\bf r})\right) \bigg/ \sum_d w^d_0({\bf r}). \label{eq:e4} [/math]

    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 Pixelization Artifacts section). In this combination, the respective number of hits of each detector in each pixel is considered, as well as the weighting (generally proportional to the inverse noise variance) applied to each detector in order to minimize the final noise.

  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 into the MC simulations.

QuickBeam[edit]

We now describe the "QuickBeam" approach (for more details, see Planck-2013-VII[1]).

Formalism[edit]

Planck observes the sky after convolution with a “scanning beam,” which captures its effective response to the sky as a function of displacement from the nominal pointing direction. Decomposing the scanning beam into harmonic coefficients Bℓm, each time-ordered data (TOD) sample can be modelled as (neglecting the contribution from instrumental noise, which is independent of beam asymmetry)

[math] %T_i = \sum_{\ell ms} D^{l}_{-m s} (\phi_i, \theta_i, \alpha_i) b_{\ell s} (-1)^{m) T_{\ell m} + n_i, T_i = \sum_{\ell ms} e^{-{\rm i} s \alpha_i} B_{\ell s} \tilde{T}_{\ell m} {}_s Y_{\ell m}(\theta_i, \phi_i), \label{eqn:tod_beam} [/math]

where the TOD samples are indexed by i, and [math]\tilde{T}_{\ell m}[/math] is the underlying sky signal. The spin spherical harmonic sYℓm rotates the scanning beam to the pointing location (θ, φ), while the e-isαi factor gives it the correct orientation. The above equation may be evaluated with the “TotalConvolver” algorithm of Wandelt and Gorski (2001)[17], accelerated using the “conviqt” recursion relations of Prezeau and Reinecke (2010)[18]. This approach is implemented in Level S, although because it involves working with a TOD-sized object, it is necessarily slow.

On small angular scales comparable to the size of the beam, it is a good approximation to assume that the procedure of mapmaking from TOD samples is essentially a process of binning: [math]T(p) = \sum_{i \in p} T_i / H(p),\label{eqn:map_beam_full}[/math] where H(p) is the total number of hits in pixel [math]\hat{n}[/math].

Starting with a normalized, rescaled harmonic transform of the beam Bℓm, sky multipoles [math]\tilde{T}_{\ell m}[/math] and a scan history object [math]w(\hat{n}, s)[/math] given by [math]w(\hat{n}, s) = \sum_{j \in p} e^{i s \alpha_j} / H(\hat{n})[/math], where the sum is over all hits j of pixel p at location [math]\hat{n}_p[/math], and αj is the scan angle for observation j. The harmonic transform of this scan-strategy object is given by

[math] {}_{s} w_{L M} = \int d^2 \hat{n} {}_s Y_{LM}^*(\hat{n}) w(\hat{n}, s). \label{eq:spin_hits} [/math]

The beam-convolved observation is then given by

[math] \tilde{T}(\hat{n}) = \sum_{slm} w(\hat{n}, -s ) B_{\ell s} T_{\ell m} {}_s Y_{\ell m}(\hat{n}). \label{eq:tilde_T} [/math]

Taking the ensemble average of the pseudo-C power spectrum of these Tℓm, we find

[math] \tilde{C}_{L}^{TT} = \sum_{S S'} \sum_{\ell_1 ell_2} \frac{(2\ell_1+1)(2\ell_2+1)}{4\pi} {}_{(-s -s')}{\cal W}_{\ell_1} B_{\ell_2 S} B_{\ell_2 S'}^* C^{TT}_{\ell_2} \left( \begin{array}{ccc} \! \ell_1\! & \ell_2\! & L\! \\ \! s\! & -s\! & 0\! \end{array} \right) \left( \begin{array}{ccc} \! \ell_1\! & \ell_2\! & L\! \\ \! s'\! & -s'\! & 0\! \end{array} \right), \label{eq:pseudo_Cl_matrix} [/math]

where [math]{}_{(s s')}{\cal W}_{L} = \frac{1}{2L+1} \sum_{M} {}_{S} w_{LM} {}_{S'} w_{LM}^*[/math] is a cross-power spectrum of scan history objects. Note that the [math]w(\hat{n}, s )[/math] that we have used here can also incorporate a position-dependent weighting to optimize the pseudo-C estimate; this might be inverse-noise or a mask—the equations are unchanged. Writing the pseudo-C in position space (a la Dvorkin and Smith, 2009[19]) with Wigner-d matrices we have

[math] \tilde{C}_{L}^{TT} = \frac{1}{8\pi} \sum_{S S'} \int_{-1}^{1} dz \ d^{L}_{00}(z) \left[\sum_{\ell_1} d^{\ell_1}_{-s -s'}(z) {}_{(-s -s')}{\cal W}_{\ell_1} (2\ell_1+1) \right] \left[\sum_{\ell_2} d^{\ell_2}_{s s'}(z) B_{\ell_2 S} B_{\ell_2 S'}^* C^{TT}_{\ell_2}(2\ell_2+1) \right]. \label{eq:pseudo_Cl_corr} [/math]

This integral can be implemented exactly using Gauss-Legendre quadrature, with a cost of order (ℓmaxsmax)2. 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 weighting is straightforward. The cost to compute all of the necessary terms exactly in that case becomes of order (ℓmaxsmaxNdet)2.

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 (encapsulated in the scan history objects described above).

A scan strategy that 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≫ℓ1, with C2 and B2 slowly-varying functions in ℓ2, and using the equality

[math] \sum_{\ell_2} (2 \ell_2+1) \left( \begin{array}{ccc} \! \ell_1\! & \ell_2\! & L\! \\ \! s\! & -s\! & 0\! \end{array} \right) \left( \begin{array}{ccc} \! \ell_1\! & \ell_2\! & L\! \\ \! s'\! & -s'\! & 0\! \end{array} \right) = \delta_{ss'}, [/math]

the pseudo-C sum above can be approximated as

[math] {\tilde{C}}_L^{TT} = C_L^{TT} \sum_{M} \left\lt \left| w(\hat{n}_p, M) \right|^2 \right\gt _p |B_{L M}|^2, \label{eq:pseudo_Cl_approx} [/math] where the average [math]\langle \rangle_p[/math] 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,

[math] \left\lt \left| w(\hat{n}, M) \right|^2 \right\gt _p = 1, [/math]

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,

\begin{equation} \left< \left| w(\hat{n}, M) \right|^2 \right>_p = \delta_{M0}, \end{equation}

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 fsky factor, as expected, but neglects the coupling between L multipoles (which can be calculated with the more complete equations above).

Pixelization artefacts[edit]

We now discuss pixelization; for more details, see Planck-2013-VII[1].

Planck-HFI maps are produced at HEALPix resolution 11 (Nside=2048), corresponding to pixels with a typical dimension of 1.7arcmin. With the resolution comparable to the spacing between scanning rings, 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 the figure below.

Figure 3: Illustration of TOI samples near the Galactic plane (grey dots), over-plotted on a simulated CMB realization which has been convolved with a 7' Gaussian FWHM beam and pixelized at Nside=2048. Associated scanning rings (grey lines) as well as centres of mass for the hit distribution (black arrows) are also plotted.

The collaboration has produced three codes that may be used to simulate the effect of pixelization on the observed sky: LevelS/TotalConvoler/Conviqt[17][18][20]; FeBeCoP[13]; and FICSBell.

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 as well as 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 [math]\langle |\nabla T |^2 \rangle = \frac{1}{4\pi} \sum_{\ell} \ell(\ell+1)(2\ell+1) C_\ell^{TT} W_\ell \approx 1\times10^9 \mu K^2 / {\rm rad}^2,[/math] where 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 [math]\langle |\nabla^2 T |^2 \rangle = \frac{1}{4\pi} \sum_{\ell} [\ell(\ell+1)]^2(2\ell+1) C_\ell^{TT} W_\ell \approx 7\times10^{14} \mu K^2 / {\rm rad}^4.[/math] 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, we can denote the displacement from the pixel centre as d=dθ+idφ. The average over all hits within a pixel gives an overall deflection vector, which we will denote for a pixel centre located at [math]\hat{n}[/math] as [math]d(\hat{n})[/math]. This represents the centre of mass of the hit distribution; Fig. 2 shows these average deflections using black arrows. The deflection field [math]d(\hat{n})[/math] may be decomposed into spin-1 spherical harmonics as [math]d_{\ell m} = \int_{4\pi} {}_1 Y_{\ell m}^* d(\hat{n}).[/math] 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

[math] C_\ell^{TT} = [1-\ell(\ell+1)R^d] {C}_\ell^{TT} W_\ell + \frac{1}{2} \sum_{\ell_1 \ell_2} \frac{\ell_1(\ell_1+1)(2\ell_1+1)(2\ell_2+1)}{4\pi} \left( \begin{array}{ccc} \! \ell_1\! & \ell_2\! & \ell\! \\ \! \ell\! & -\ell\! & 0\! \end{array} \right)^2 C_{\ell_1}^{TT} W_{\ell_1} \left[ C_{\ell_2}^{d+} + (-1)^{\ell + \ell_1 + \ell_2} C_{\ell_2}^{ d-} \right], \label{eqn:clt_pixelized} [/math]

where Rd=〈|d|2〉/2 is half the mean-squared deflection magnitude (averaged over hits within a pixel, as well as over pixels). [math]C_\ell^{d+}[/math] is the sum of the gradient and curl power spectra of dℓm, and [math]C_\ell^{d-}[/math] is the gradient spectrum minus the curl spectrum. The Rd term describes a smearing of the observed sky due to pixelization. For uniform pixel coverage of Nside=2048 pixels √〈|d|2〉=0.725', while, for the hit distribution of Planck frequency maps, Rd is typically 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 essentially 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 the Planck 2013 gravitational lensing paper (Planck-2013-XVII[21]), where the sub-pixel deflection field constitutes a potential source of bias for the measured lensing potential. Indeed, the above equation is just a slightly modified version of the usual first-order CMB lensing power spectrum (see Hu 2000[22]; Lewis & Challinor 2006[23]) to accommodate curl modes.

A useful approximation to this equation, 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 [math]d(\hat{n})[/math] couples the CMB gradient into a source of noise with an effective level given by [math]\sigma^{N} \approx \sqrt{ R^{TT} \frac{4\pi}{N_{\rm pix}} \langle | d(\hat{n}) |^2 \rangle }, % (\muKarcmin ),[/math]

where the average is taken over all pixels and RTT is half the mean-squared power in the CMB gradient: [math]R^{TT} = \frac{1}{8\pi} \sum_{\ell} \ell(\ell+1)(2\ell+1) \tilde{C}_\ell^{TT}.[/math] For frequency-combined maps, [math]\sqrt{ \langle | d(\hat{n}) |^2 \rangle }[/math] is typically on the order of 0.1', and so the induced noise is at the level of σN≈2μK.arcmin. This is small compared to the instrumental contribution, although it does not disappear when taking cross-spectra, depending on how coherent the hit distributions of the two maps in the cross-spectrum are.

References[edit]

  1. 1.01.11.21.31.41.51.61.71.81.9 Planck 2013 results. VII. HFI time response and beams, Planck Collaboration, 2014, A&A, 571, A7.
  2. Planck 2015 results. VII. High Frequency Instrument data processing: Time-ordered information and beam processing, Planck Collaboration, 2016, A&A, 594, A7.
  3. Planck 2018 results. III. High Frequency Instrument data processing and frequency maps, Planck Collaboration, 2020, A&A, 641, A3.
  4. Planck 2018 results. V. CMB Power Spectra and Likelihoods, Planck Collaboration, 2020, A&A, 641, A5.
  5. Planck 2015 results. VIII. High Frequency Instrument data processing: Calibration and maps, Planck Collaboration, 2016, A&A, 594, A8.
  6. Planck intermediate results. LII. Planets flux densities, Planck Collaboration Int. LII A&A, 607, A122, (2017)
  7. 7.07.1 Planck 2013 results. VI. High Frequency Instrument Data Processing, Planck Collaboration, 2014, A&A, 571, A6.
  8. Planck 2013 results. X. HFI energetic particle effects: characterization, removal, and simulation, Planck Collaboration, 2014, A&A, 571, A10.
  9. Planck 2013 results. XV. CMB power spectra and likelihood, Planck Collaboration, 2014, A&A, 571, A15.
  10. Planck pre-launch status: HFI beam expectations from the optical optimisation of the focal plane, B. Maffei, F. Noviello, J. A. Murphy, P. A. R. Ade, J.-M. Lamarre, F. R. Bouchet, J. Brossard, A. Catalano, R. Colgan, R. Gispert, E. Gleeson, C. V. Haynes, W. C. Jones, A. E. Lange, Y. Longval, I. McAuley, F. Pajot, T. Peacocke, G. Pisano, J.-L. Puget, I. Ristorcelli, G. Savini, R. Sudiwala, R. J. Wylde, V. Yurchenko, A&A, 520, A12+, (2010).
  11. Planck pre-launch status: The optical system, J. A. Tauber, H. U. Nørgaard-Nielsen, P. A. R. Ade, J. Amiri Parian, T. Banos, M. Bersanelli, C. Burigana, A. Chamballu, D. de Chambure, P. R. Christensen, O. Corre, A. Cozzani, B. Crill, G. Crone, O. D'Arcangelo, R. Daddato, D. Doyle, D. Dubruel, G. Forma, R. Hills, K. Huffenberger, A. H. Jaffe, N. Jessen, P. Kletzkine, J. M. Lamarre, J. P. Leahy, Y. Longval, P. de Maagt, B. Maffei, N. Mandolesi, J. Martí-Canales, A. Martín-Polegre, P. Martin, L. Mendes, J. A. Murphy, P. Nielsen, F. Noviello, M. Paquay, T. Peacocke, N. Ponthieu, K. Pontoppidan, I. Ristorcelli, J.-B. Riti, L. Rolo, C. Rosset, M. Sandri, G. Savini, R. Sudiwala, M. Tristram, L. Valenziano, M. van der Vorst, K. van't Klooster, F. Villa, V. Yurchenko, A&A, 520, A2+, (2010).
  12. Planck 2015 results. XI. CMB power spectra, likelihoods, and robustness of cosmological parameters, Planck Collaboration, 2016, A&A, 594, A11.
  13. 13.013.113.2 Fast Pixel Space Convolution for Cosmic Microwave Background Surveys with Asymmetric Beams and Complex Scan Strategies: FEBeCoP, S. Mitra, G. Rocha, K. M. Górski, K. M. Huffenberger, H. K. Eriksen, M. A. J. Ashdown, C. R. Lawrence, ApJS, 193, 5-+, (2011).
  14. Three-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Temperature Analysis, G. Hinshaw, M. R. Nolta, C. L. Bennett, R. Bean, O. Doré, M. R. Greason, M. Halpern, R. S. Hill, N. Jarosik, A. Kogut, E. Komatsu, M. Limon, N. Odegard, S. S. Meyer, L. Page, H. V. Peiris, D. N. Spergel, G. S. Tucker, L. Verde, J. L. Weiland, E. Wollack, E. L. Wright, ApJS, 170, 288-334, (2007).
  15. Detection of gravitational lensing in the cosmic microwave background, K. M. Smith, O. Zahn, O. Doré, Phys. Rev. D, 76, 043510-+, (2007).
  16. PReBeaM for Planck: A Polarized Regularized Beam Deconvolution Map-Making Method, C. Armitage-Caplan, B. D. Wandelt, ApJS, 181, 533-542, (2009).
  17. 17.017.1 Fast convolution on the sphere, Benjamin D. Wandelt, Krzysztof M. Gorski, Phys. Rev., D63, 123002, (2001).
  18. 18.018.1 Algorithm for the evaluation of reduced Wigner matrices, Gary Prezeau, Martin Reinecke, Astrophys. J. Suppl., 190, 267-274, (2010).
  19. Reconstructing Patchy Reionization from the Cosmic Microwave Background, Cora Dvorkin, Kendrick M. Smith, Phys. Rev., D79, 043003, (2009).
  20. A simulation pipeline for the Planck mission, M. Reinecke, K. Dolag, R. Hell, M. Bartelmann, T. A. Enßlin, A&A, 445, 373-373, (2006).
  21. Planck 2013 results. XVII. Gravitational lensing by large-scale structure, Planck Collaboration, 2014, A&A, 571, A17.
  22. Weak lensing of the CMB: A harmonic approach, W. Hu, Phys. Rev. D, 62, 043007, (2000).
  23. Weak gravitational lensing of the CMB, A. Lewis, A. Challinor, Phys. Rep., 429, 1-65, (2006).

(Planck) High Frequency Instrument

Data Processing Center

Cosmic Microwave background

reduced IMO

(Hierarchical Equal Area isoLatitude Pixelation of a sphere, <ref name="Template:Gorski2005">HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere, K. M. Górski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, M. Bartelmann, ApJ, 622, 759-771, (2005).

Noise Equivalent Temperature

Full-Width-at-Half-Maximum

(Planck) Low Frequency Instrument