© The Authors 2023

1. Introduction

Among the most important systematic effects that must be accounted for when observing the cosmic microwave background (CMB) are stray light and sidelobes (e.g., Barnes et al. 2003; Planck Collaboration III 2014). These effects come from the non-zero response of the detector to areas of the sky outside the main beam, however that may be defined for each particular case. Because microwave telescopes typically work near their diffraction limit, the presence of a certain extent of sidelobe effects is inevitable. Furthermore, their structure can be complicated by many different physical effects, such as spurious optical reflections or manufacturing irregularities in the detectors or optical elements. These signal contributions can have far-reaching consequences on the observed signal, particularly at large angular scales, as they do not behave in the same sky-stationary manner as the main beam signal.

Sidelobe signals can produce many types of errors in CMB analysis pipelines and they represent a potent source of systematic contamination (e.g., Planck Collaboration III 2016; Watts et al. 2023). In particular, as the sidelobe response functions often are broadly distributed, this contamination can confuse important signals such as the CMB solar and orbital dipoles that are used for calibration. Sidelobes uncertainties are coupled directly with foreground emission from diffuse galactic components, producing a significant source of contamination. In some experiments a further contaminating signal can originate from a source not on the sky, such as ground pickup or radio-frequency (RF) noise. In all cases, sidelobe signals are detrimental to the quality of the final sky maps and parameter estimates, and a dedicated effort is required to remove them. The process of characterizing and correcting these spurious signals is therefore an important part of optimal CMB mapmaking, and it requires optimized algorithms to do so efficiently. One of the most important of these steps is to convolve a beam or sidelobe response function with a sky map or model to generate a re-observed map.

Full-sky convolution on the sphere is a problem that has been important in the CMB field since the earliest satellite measurements. Early experiments like COBE (Toral et al. 1989) either did not model the sidelobes or they used simple pixel-based convolution approaches that (even at low resolution) still required radially symmetric beam approximations (Wu et al. 2001), or limited the applications to large scales (Burigana et al. 2001).

Wandelt & Górski (2001) presented the first harmonic space convolution algorithm, often referred to as “total convolution”, which achieved a large performance gain over pixel-based methods, by as much as a factor of $O(\sqrt{N_{\mathrm{pix}}})$. This breakthrough allowed for the calculation of these convolutions easily enough that they could be applied to each simulation, instead of requiring a dedicated study necessitating months of runtime.

Next, Prézeau & Reinecke (2010) developed the Conviqt approach, which was used both by several official Planck analysis pipelines (Planck Collaboration III 2016, Planck Collaboration Int. XLVI 2016, Planck Collaboration Int. LVII 2020) and to generate the Planck Full Focal Plane (FFP) simulations (Planck Collaboration XII 2016). This approach was an improvement over the state of the art, speeding up the computation of the Wigner recursion relationships used in the original harmonic space algorithm, as well as providing a standardized, user-friendly library, libconviqt, which could be incorporated into numerous pipelines.

In this paper, we introduce a new formulation of the Conviqt algorithm that is based on spherical harmonic transforms (SHTs), rather than directly computing the Wigner matrix elements. We are thus able to leverage the highly optimized libsharp SHT library to perform the bulk of the calculations (Reinecke & Seljebotn 2013). Although this new approach was not developed specifically for BEYONDPLANCK, this paper is the first to explicitly derive, discuss, and benchmark the method.

2. Sidelobes, libconviqt, and libsharp

2.1. Total convolution through spin harmonics

Based on a given sky map, $s(\widehat{n})$, and beam, $b(\widehat{n})$, our task is to compute a quantity c(ϑ, φ, ψ) ∈ ℝ that represents the convolution of these two fields, with the beam oriented in polar coordinates (ϑ, φ)¹ and rotated around its own central axis by ψ,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill c(\vartheta ,\phi ,\psi )& :=s_{l{m}_s}\ast b_{l{m}_b}\hfill \\ \hfill & \equiv {\displaystyle {\int }_{4\pi }s(\widehat{n})b(\widehat{n}\prime (\vartheta ,\phi )-\widehat{n},\psi )\mathrm{d}{\mathrm{\Omega }}_{\widehat{n}}.}\hfill \end{array}\end{array}$$(1)

Here, s_lms denotes the spherical harmonic coefficients of the sky signal and b_lmb is the beam in the same representation. In this expression, care must been taken to distinguish between the sky and beam bandlimits, m_s and m_b, as the two indices will be treated separately in the following derivation.

A computationally efficient solution for this problem was derived by Prézeau & Reinecke (2010), who exploited fast recurrence relations for Wigner d matrix elements to evaluate Eq. (1) in harmonic space. In the following, we show that this equation can alternatively be expressed in terms of spin-harmonics. The resulting algebra is in principle identical to the recursion relations used by Prézeau & Reinecke (2010), but the euqations are simply repackaged in a format that is significantly easier to implement in practical computer code, since it may use existing and highly optimized spherical harmonic libraries, such as libsharp (Reinecke & Seljebotn 2013), to perform the computationally expensive parts.

As shown by Wandelt & Górski (2001), Eq. (1) can be evaluated efficiently in harmonic space as

$$\begin{array}{c}\hfill c(\vartheta ,\phi ,\psi )={\displaystyle \sum _{l,m_s,m_b}{s}_{l{m}_s}{b}_{l{m}_b}^{\ast }{[D_{m_s{m}_b}^l(\phi ,\vartheta ,\psi )]}^{\ast },}\end{array}$$(2)

where s_lms and b_lmb are the spherical harmonic coefficients of the signal and beam, respectively, and $D_{m_s{m}_b}^l$ is the Wigner D-matrix.

This may be expressed as (Goldberg et al. 1967)

$$\begin{array}{c}\hfill D_{-ms}^l(\phi ,\vartheta ,-\psi )={(-1)}^m{\sqrt{\frac{4\pi }{2l+1}}}_s{Y}_{\mathit{lm}}(\vartheta ,\phi )e^{is\psi },\end{array}$$(3)

where _sY_lm(ϑ, φ) is the spin-weighted spherical harmonic and the placement of the negative signs are an arbitrary historical convention. Inserting this expression into Eq. (2) yields:

$$\begin{array}{c}\hfill c(\vartheta ,\phi ,\psi )={\displaystyle \sum _{l,m_s,m_b}\sqrt{\frac{4\pi }{2l+1}}{s}_{l{m}_s}{b}_{l-m_b}{·}_{-m_b}{Y}_{l{m}_s}(\vartheta ,\phi )e^{i{m}_b\psi },}\end{array}$$(4)

where we have assumed that the beam is real-valued in position space and we have used the symmetry relations:

$$\begin{array}{cc}\hfill D_{-m_s,-m_b}^l(\vartheta ,\phi )& ={(-1)}^{m_s+m_b}{[D_{m_s,m_b}^l(\vartheta ,\phi )]}^{\ast },\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill b_{l,m_b}^{\ast }{(-1)}^{m_b}& =b_{l,-m_b}.\hfill \end{array}$$(6)

Pulling the summation over m_b in front of the other sums yields

$$\begin{array}{c}\hfill c(\vartheta ,\phi ,\psi )={\displaystyle \sum _{m_b}{e}^{i{m}_b\psi }\sum _{l,m_s}\sqrt{\frac{4\pi }{2l+1}}{s}_{l{m}_s}{b}_{l-m_b}{·}_{-m_b}{Y}_{l{m}_s}(\vartheta ,\phi ).}\end{array}$$(7)

The terms in this outer sum can be arranged in the form m_b = 0, ±1, ±2, …. The contribution from m_b = 0 can be interpreted as a spin-0 spherical harmonic transform of the quantity $\sqrt{4\pi /(2l+1)}{s}_{l{m}_s}{b}_{l0}$, which can be easily computed by a library such as libsharp (Reinecke & Seljebotn 2013).

Since c(ϑ, φ, ψ) ∈ ℝ, we know that the contributions from the pairs m_b = ±1, ±2, … must be complex conjugate with respect to each other, and their combined contribution is therefore

$$\begin{array}{cc}\hfill e_{m_b}^{i{m}_b\psi }& S(\vartheta ,\phi )+e_{m_b}^{-i{m}_b\psi }{S}^{\ast }(\vartheta ,\phi )=\hfill \\ \hfill & 2[cos(m_b\psi )\text{Re}{(}_{m_b}S(\vartheta ,\phi ))-sin(m_b\psi )\text{Im}{(}_{m_b}S(\vartheta ,\phi ))],\hfill \end{array}$$(8)

where we have defined

$$\begin{array}{c}\hfill {}_{m_b}S(\vartheta ,\phi )\equiv {\displaystyle \sum _{l,m_s}\sqrt{\frac{4\pi }{2l+1}}{s}_{l{m}_s}{b}_{l-m_b}{·}_{-m_b}{Y}_{l{m}_s}(\vartheta ,\phi ).}\end{array}$$(9)

This is a spherical harmonic transform of a quantity with a spin m_b, which can also be computed efficiently by libsharp.

In practice, the transforms in Eq. (9) are implemented by separating S into its gradient and curl (or E and B) coefficients, a_lm (Lewis 2005):

$$\begin{array}{c}\hfill {}_{m_b}{S}_{l{m}_s}=-({}_{m_b}{E}_{l{m}_s}+i{}_{m_b}{B}_{l{m}_s}),\end{array}$$(10)

using the symmetry relations ${}_{m_b}{E}_{l-m_s}={(-1)}^{m_s}{}_{m_b}{E}_{l{m}_s}^{*}$ and ${}_{m_b}{B}_{l-m_s}={(-1)}^{m_s}{}_{m_b}{B}_{l{m}_s}^{*}$, and where the overall minus sign is a convention. Again making use of the symmetry relation in Eq. (6), this gives us:

$$\begin{array}{cc}\hfill {}_{m_b}{E}_{l,m_s}& =-s_{l{m}_s}\text{Re}(b_{l,m_b})\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill {}_{m_b}{B}_{l,m_s}& =-s_{l{m}_s}\text{Im}(b_{l,m_b}).\hfill \end{array}$$(12)

To summarize, efficient evaluation of the convolution integral in Eq. (1) may be done through the following steps: Firstly, for each m = 0, …, m_b, we pre-compute the spin spherical harmonic coefficients in Eqs. (11)–(12) and then compute the corresponding spin-m_b SHT with an external library such as libsharp. This results in a three-dimensional data cube of the form c(ϑ, φ, m_b). Then for each position on the sky, (ϑ, φ), we perform a Fourier transform to convert these coefficients to c(ϑ, φ, ψ), as given by Eq. (7).

In practice, the resulting c(ϑ, φ, ψ) data object is evaluated at a finite pixel resolution typically set to match the beam band limit. To obtain smooth estimates within this data object, a wide range of interpolation schemes may be employed, trading computational efficiency for accuracy. This issue is identical to what has been faced in previous approaches (Wandelt & Górski 2001; Prézeau & Reinecke 2010).

2.2. Comparison with libconviqt

To compare the results of this new total convolution approach with the older libconviqt approach of Prézeau & Reinecke (2010), we evaluated the convolution between the beam for one of the LFI 30 GHz receivers (28M) and a Commander 30 GHz sky model (Andersen et al. 2023; Svalheim et al. 2023) using both methods. The resulting convolution cubes were then observed using LFI’s scanning strategy for the first year of the Planck flight. The resulting map differences are shown in Fig. 1. The convolution cubes were also directly compared for accuracy, and found to agree with an integrated difference at the 10⁻⁸ level, indicative of differences at the level of numerical precision.

Fig. 1. Map level difference in temperature of the new SHT convolution algorithm compared to the old Conviqt approach, observed using the identical pointing of the first year of the 30 GHz Planck detectors. The differences are at the level of machine precision, indicating full agreement between the two algorithms.

Figure 2 compares the runtime between the two approaches for a test configuration with an elliptical beam and a fixed sky model, and with m_max = 0 (using only the radially symmetric part of the beam) and m_max = 10, respectively. In both cases, the new implementation outperforms the old approach at all but the lowest values of l_max, where the data read time dominates. Additionally, for compatibility with the old libconviqt approach, this test was performed with an older version of libsharp, so we expect that the new algorithm scales even more favourably than this with the latest implementation. We note that this is a significant real-life advantage of the new approach: any improvement in SHT libraries, which typically are subject to intensive algorithm development and code maintenance, translates directly into a computational improvement for the convolution algorithm.

Fig. 2. Runtime comparison between the libconviqt approach and the new spin-SHT approach for the convolution of an elliptical Gaussian with a set of random sky al, ms. This work ties or outperforms the previous approach for all values of lmax from 256 to 8192 for both mmax values shown. We note the log scale on the y-axis.

3. Sidelobe models

Figure 3 shows characteristic sidelobe response functions evaluated at a fixed frequency on the sky for a detector in each Planck LFI band. The sidelobe response for each detector within a single Planck band look visually quite similar, so only these representative ones are shown here. Each one is stored on disk as a set of a_lm’s with l_max = 512 and m_max = 100.

Fig. 3. Maps of the sidelobe response on the sky from a representative detector at (left to right) 30 GHz, 44 GHz, and 70 GHz. The beam orientation is such that the main beam is pointed directly at the north pole in these maps. The intensities are normalized such that the main beams have unit power at l = 0.

3.1. Main beam treatment

In the BEYONDPLANCK analysis, the sidelobe and main beam components of the sky response are separated, and the sidelobes are treated as a nuisance signal similar to the orbital dipole correction term, as can be seen in the BEYONDPLANCK global parametric model of the data:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill d_{j,t}& =g_{j,t}[{\mathsf{P}}_{tp,j}{\mathsf{B}}_{p{p}^{\prime },j}{\displaystyle \sum _c}{\mathsf{M}}_{\mathit{cj}}({\beta }_{p\prime },\mathrm{\Delta }b{p}^j)a_{p^{\prime }}^c+s_{j,t}^{\mathrm{orb}}+s_{j,t}^{\mathrm{fsl}}]\hfill \\ \hfill & +s_{j,t}^{1\mathrm{hz}}+n_{j,t}^{\mathrm{corr}}+n_{j,t}^{\mathrm{w}}.\hfill \end{array}\end{array}$$(13)

The other terms in this equation are discussed in detail in BeyondPlanck Collaboration (2023), but here the main beam signal is denoted as B_pp′, j and the sidelobe signal is extracted from the signal contribution and expressed as $s_{j,t}^{\mathrm{fsl}}$. This distinction allows the sidelobes to be treated separately from the main beam in all respects. Treating the main beam using the Conviqt formalism of this paper would be possible, but the additional precision needed to model it accurately would require much higher values of l_max, and, therefore, greatly increased computational time and memory requirements.

In the BEYONDPLANCK analysis, the main beam is used (in conjunction with the sidelobes) to compute the full 4π dipole response, as detailed in Sect. 3.3. Additionally, a Gaussian main beam approximation is used during component separation to smooth the sky model to the appropriate beam resolution for each channel. During mapmaking, beam effects are ignored and the beam is assumed to be pointed at the center of each pixel.

3.2. Sidelobe normalization

We adopted a normalization of the sidelobes that differs slightly from the normalization used within the Planck LFI collaboration. The Planck 2018 LFI beam products leave a small portion (around 1%) of known missing power within the system unassigned due to uncertainties about to which component it should be assigned (Planck Collaboration IV 2016). In the current analysis, we adopted the same approximation as for Planck DR4 (Planck Collaboration Int. LVII 2020) and renormalized the beam transfer function such that this power would be distributed proportionally at each l; that is, we rescaled the beam transfer function B_l so that its full sky integral would be B₀ = 1. This re-scaling is equivalent to assigning the unknown beam power uniformly over the entire beam. We note that this normalization is, in either case, always done before any higher-level analysis for both Planck 2018 and DR4; the only difference is whether the renormalization ought to be performed by external users through a deconvolution of a non-unity normalized main beam transfer function or not.

3.3. Orbital dipole and quadrupole sidelobe response

The treatment of the sidelobes is also important when generating orbital dipole and quadrupole estimates. Because Planck is calibrated primarily from the dipole measurements (Planck Collaboration I 2020; Planck Collaboration Int. LVII 2020; Gjerløw et al. 2023), the sidelobe’s contribution to the dipole can directly result in an absolute calibration error if not handled appropriately. While the CMB Solar dipole can easily be handled using the Conviqt approach described in Sect. 2.1, the orbital dipole is not sky-stationary and thus must be handled separately.

BEYONDPLANCK generates orbital dipole and quadrupole estimates directly from the Planck pointing information, using the satellite velocity data that has been stored at low resolution (one measurement per pointing period). With this information, it is possible to estimate the orbital dipole and quadrupole amplitude for each timestep, allowing the time-domain removal of the signal before it contaminates the final products with non-sky-stationary signal artifacts. Additionally, once this signal has been isolated from the raw data, it can be used as an aid in the calibration routines because of its highly predictable structure.

BEYONDPLANCK uses the same technique as Planck DR4 (see Appendix C of Planck Collaboration Int. LVII 2020) to generate the orbital dipole and quadrupole estimate. That is, we can express the signal $\stackrel{\sim }{D}$ seen by a detector observing a fixed direction, ${\widehat{\mathit{n}}}_0$, as the convolution of the dipole and quadrupole signal on the sky, $D(\widehat{\mathit{n}})$, with the full 4π beam response, $B(\widehat{\mathit{n}},{\widehat{\mathit{n}}}_0)$,

$$\begin{array}{c}\hfill \stackrel{\sim }{D}({\widehat{\mathit{n}}}_0)={\displaystyle \int }\mathrm{d}\mathrm{\Omega }B(\widehat{\mathit{n}},{\widehat{\mathit{n}}}_0)D(\widehat{\mathit{n}}).\end{array}$$(14)

Here it is useful to break the dipole signal up into three orthogonal components in the standard Cartesian coordinates (x, y, z) and we adopt the convention that the main beam points toward the north pole in our coordinate system.

The orbital CMB dipole and quadrupole can be expressed as a Doppler shift in each direction (Notari & Quartin 2015):

$$\begin{array}{c}\hfill D(\widehat{\mathit{n}})=T_0[\beta ·\widehat{\mathit{n}}(1+q\mathit{\beta }·\widehat{\mathit{n}})],\end{array}$$(15)

where β is the satellite velocity divided by the speed of light $\beta =\frac{\mathit{v}}{c}$; T₀ is the CMB temperature and q is the quadrupole factor dependent on the frequency ν, defined by

$$\begin{array}{c}\hfill q=\frac{a}{2}\frac{e^a+1}{e^a-1},\mathrm{where}a=\frac{h\nu }{k_{\mathrm{B}}{T}_0}.\end{array}$$(16)

Inserting these expressions into Eq. (15), we obtain:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \stackrel{\sim }{D}& =T_0{\displaystyle \int }\mathrm{d}{\mathrm{\Omega }}_{\widehat{\mathit{n}}}B(\widehat{\mathit{n}},{\widehat{\mathit{n}}}_0)[x{\beta }_x+y{\beta }_y+z{\beta }_z\hfill \\ \hfill & +q(x^2{\beta }_x^2+y^2{\beta }_y^2+z^2{\beta }_z^2\hfill \\ \hfill & +2xy{\beta }_x{\beta }_y+2xz{\beta }_x{\beta }_z+2yz{\beta }_y{\beta }_z)],\hfill \end{array}\end{array}$$(17)

where $\widehat{\mathit{n}}=(x,\mathit{y},z)$ is a unit direction vector that is also the integration variable and ${\widehat{\mathit{n}}}_0$ is the fixed direction of the satellite pointing for this timestep. We note that the geometric factors in this expression may be precomputed as

$$\begin{array}{c}\hfill S_x={\displaystyle \int }xB(\widehat{\mathit{n}},{\widehat{\mathit{n}}}_0)\mathrm{d}{\mathrm{\Omega }}_{\widehat{\mathit{n}}},\end{array}$$(18)

and, thus, we see that Eq. (17) may be written in the following form:

$$\begin{array}{cc}\hfill \stackrel{\sim }{D}& =T_0[S_x{\beta }_x+S_y{\beta }_y+S_z{\beta }_z\hfill \\ \hfill & +q(S_{\mathit{xx}}{\beta }_x^2+S_{\mathit{yy}}{\beta }_y^2+S_{\mathit{zz}}{\beta }_z^2\hfill \\ \hfill & +2S_{\mathit{xy}}{\beta }_x{\beta }_y+2S_{\mathit{xz}}{\beta }_x{\beta }_z+2S_{\mathit{yz}}{\beta }_y{\beta }_z)].\hfill \end{array}$$(19)

To compute $\stackrel{\sim }{D}$ for an arbitrary beam orientation, we simply need to rotate the satellite pointing and velocity vectors into the coordinate system used to define S, and then we can evaluate Eq. (19) very quickly.

BEYONDPLANCK further accelerates this operation by computing this rotation for only one point in twenty (chosen so as to still fully sample the dipole) and using a spline to interpolate between them. This saves the cost of calculating a new rotation matrix at each step, and instead relies on the smoothness of the signal to ensure continuity. The algorithm treats the final few points of each pointing period that do not divide evenly into the subsampling factor separately. This allows for the use of regular bin widths, which greatly speeds up the splining routines, while the final few points are calculated using the slower rotation matrix technique.

4. Sidelobe estimates

4.1. Posterior mean corrections

Figure 4 shows the mean sidelobe signal estimates at each of the three LFI frequencies for the entire mission, co-added across each frequency and projected into sky coordinates, identically to the way the true sky signal is treated. Each map is averaged over 90 Gibbs samples produced in the main BEYONDPLANCK analysis (BeyondPlanck Collaboration 2023), after discarding burn-in and a thinning of the remaining chain by a factor of ten.

Fig. 4. Maps of the sidelobes convolved with the sky at each of the three LFI frequencies. From top to bottom: 30 GHz, 44 GHz and 70 GHz. Left column: unpolarized sky signal, central column: Q polarization and right column: U. Note the difference in the colour scales required to see the same level of detail in all three channels.

We note that these maps follow the traditional Planck LFI method of sidelobe correction by producing these signals in the time domain during TOD processing. These templates are therefore exactly correct for the maps produced by these pipeline runs, but will not precisely match up with analyses that use different data cuts, flagging or channel selection.

These results appear visually similar to the corresponding Planck DPC results presented in Fig. 7 of Planck Collaboration III (2016). The main difference is that the current results also include the CMB dipole, whereas the LFI 2015 DPC analysis showed the sidelobe pickup of dipole-subtracted maps. We see that the sidelobe signal is strongest at 30 GHz and that the dominant features in the co-added sky maps consist of a series of rings created by the interplay between the sidelobe pickup and bright Galactic plane features.

Figure 4 also clearly indicates that the overall level of sidelobe pickup at 44 GHz is significantly lower than for the 30 and 70 GHz channels. This is due to the particular location in the focal plane of two of the three 44 GHz feedhorns, which results in a significant under-illumination of both the primary and secondary reflectors of the Planck telescope for those two horns (see Fig. 4 of Sandri et al. 2010).

4.2. Error propagation

In addition to the posterior mean sidelobe maps, the BEYONDPLANCK pipeline outputs also provide an estimate of the sidelobe stability and statistical variation. Figure 5 shows the rms maps generated from the same sample of sidelobe signal estimates used in Fig. 4. Clear evidence of the scanning pattern can be seen, as expected. The sharp vertical lines visible in polarization (clearest in 30 GHz Q and U at the top, and 44 GHz U at the top and bottom) have been previously examined by the Planck team and they are caused by a chance alignment between the non-dense Planck scanning strategy and the shape of the HEALPix pixels. For an example of this effect, see Fig. 15 of Planck Collaboration VII (2014).

Fig. 5. Sidelobe rms maps at each of the three LFI frequencies. From top to bottom: 30 GHz, 44 GHz, and 70 GHz. Left column: unpolarized sky signal, central column: Q polarization and right column: U. Again note the different colour scales.

These posterior rms maps cannot be considered true sidelobe error estimates, however, since they do not account for uncertainties in the sidelobe response itself. Rather, they only show the change in the estimated sidelobe signal due to sky model variations from component separation. Full sidelobe error propagation would require sampling over the physical parameters that determine the detectors’ sidelobe response on the sky. Sampling the full set of optical model parameters is likely to be infeasible due to excessive computational time, however, identifying a minimal parameter set that may account for the main potential variations in the sidelobe response functions, as well as precomputing response functions over a grid of such parameters, would result in physically motivated uncertainties for the sidelobe models. This approach will be developed for future applications such as the LiteBIRD mission (Hazumi et al. 2019).

5. Impact on frequency and component maps

To assess the importance of sidelobe corrections on frequency and component maps, we performed two runs of the Commander code, starting from the same input data as the main pipeline run and with the identical random seed 12345678. As a comparison, we removed the far sidelobe correction from one of these secondary pipeline executions, and we differentiated the results between these two pipelines. The random seed is required by the sampling approach of the BEYONDPLANCK pipeline to explore the posterior distribution of many of the sampled parameters; identical seeds will result in identical samples barring other differences.

Figure 6 shows the differences in the frequency maps between the two cases in temperature, where the effects are the most obvious. The only large-scale features that can clearly be seen are the small dipole differences (most clearly visible at 70 GHz). These are directly caused by the dipolar component seen in Fig. 4, since this contribution to the total sky signal that was in the sidelobe term is now unaccounted for. In previous analyses, these dipole contributions were handled through specific modeling of exactly these effects, but this test makes it explicitly clear that correct dipole measurements require accurate knowledge of the sidelobe pickup.

Fig. 6. Frequency map difference plots smoothed to one degree at (top to bottom) 30, 44, and 70 GHz, comparing two pipeline executions with the same seed, one of which has no sidelobe correction.

Next, we see two more features in the difference maps that are more localized. The first of these are the ring structures that match the actual sidelobe map structures quite closely, which are, of course, the same rings from Fig. 4, and are not accounted for in the second pipeline run without sidelobe corrections. Additionally, there are some uniform residuals that are visible in the Galactic plane regions of the difference maps. These are caused by a calibration mismatch between the detectors at a single frequency. Since each of the detectors now sees a slightly different dipole signal on the sky, depending on its specific sidelobe response, their calibrations do not agree with one another, which causes signal residuals that are most visible in the plane where the signal amplitude is highest.

Figure 7 shows the differences in component maps from this same comparison, again in temperature. The CMB as well as the three low-frequency foreground components are estimated using the standard Commander3 technique described in Andersen et al. (2023). The AME component sees similar issues to the ones seen by the frequency maps above: the dipole is slightly incorrect, there are sidelobe-esque stripes, and the Galactic plane shows a strong residual. All of these effects have been seen directly in the frequency maps. The dipole difference seen here is precisely the one that contributes to the difference in calibration between the two different pipeline executions.

Fig. 7. Component map difference plots for (top to bottom) CMB, synchrotron, AME and freefree emission, comparing two pipeline executions with the same seed, one of which has no sidelobe correction.

The other three components (synchrotron, CMB, and free-free) show relatively fewer structural differences. They have absorbed some of the sidelobe-like ring structures, but the primary difference can be seen most clearly in the Galactic plane. Here, we notice a large residual caused by the inaccurate model of the Galactic emission being altered slightly by the gain and calibration differences between the two runs. As the Galactic emission is significantly brighter than the rest of the sky, small changes in calibration produce large errors such as the ones seen here.

Finally, in Fig. 8, we show the correlated noise map difference at 30 GHz between the two runs. We see that some of the missing sidelobe signal has been accommodated by the correlated noise component. These structures mirror the strongest sidelobe-like signals in the 30 GHz difference map in the top panel of Fig. 6. The fact that n_corr can accommodate some sidelobe residuals is, in fact, helpful, as it allows for some leeway in the final sidelobe model, in the sense that small uncertainties and artifacts that are inconsistent between different frequency channels will be mostly absorbed into the correlated noise component, rather than in the sky maps. The differences that we see in the maps, however, indicate that this process is not perfect, as some of the spurious signal still makes it to the final maps without a perfect sidelobe model. For a real-world example of these issues, we refer the interested reader to the ongoing BEYONDPLANCK re-analysis of the WMAP data, for which far sidelobe contamination appears to be a dominant problem (Watts et al. 2023).

Fig. 8. Difference in correlated noise, projected into the map domain at 30 GHz comparing two pipeline executions with the same seed, one of which has no sidelobe correction.

The residual errors seen in Figs. 6 and 7 are also present in the BEYONDPLANCK analysis, albeit at much lower levels. We know that our knowledge and modelling of the sidelobes are imperfect, as they are based on limited measurements of the physical LFI sidelobes and some of the power is unaccounted for. Future applications that aim for a robust measurement of primordial B-modes (constraining the tensor-to-scalar ratio r ≤ 0.01), will be required to marginalize over the sidelobe uncertainties in some manner, either directly via a Gibbs sampling of a subset of the instrument parameters or by parameterizing and fitting sidelobe error estimates. We do not believe that the sidelobe contribution causes significant errors in the LFI sample sets produced by BEYONDPLANCK, as it is unlikely to be more than a 10% error on the sidelobe estimates of Fig. 4. At 30 GHz, this corresponds to at most a 0.05% error in our temperature maps and a 1% error in polarization. We do expect, however, that as instrumental sensitivities improve, especially in polarization, this sidelobe term will need to be modeled very accurately, and the corresponding uncertainties must be propagated properly, for instance, using methods similar to those presented in this paper.

6. Summary and conclusions

This paper presents a formulation of the Conviqt algorithm in terms of spin-weighted spherical harmonics. This algorithm has already been already implemented in the latest versions of the libconviqt library and it has now also been re-implemented directly into Commander, where it is used for sidelobe corrections for the BEYONDPLANCK analysis framework. Based on the Monte Carlo samples produced in that analysis, we have presented novel posterior mean and standard deviation maps for each of the three Planck LFI frequency bands.

The full-sky sidelobe treatment techniques presented here are easily generalizable to other experiments, and can be tuned to match the required spatial characteristics of other instruments simply by adjusting the spherical harmonic bandpass parameters, l_max and m_max, of the sidelobe description. The only requirement for using the code with a new instrument is a HEALPix-compatible description of the sidelobe response function per detector. The more accurate this characterization of the instrument is, the closer the sidelobe estimation will get to approximating the true sidelobe contamination in the timestream.

We note that the approach presented here is less useful for ground- or balloon-based experiments where the dominant sidelobe pickup contains radiation from an environmental source. This pickup is not sky-synchronous, and, thus, it cannot be modeled purely as a beam-sky convolution, but instead must include additional contributions from, for example, telescope baffles, ground pickup, or clouds. For these types of experiments, other techniques such as aggressive baffling are likely to be better suited.

We also stress that the current implementation only supports sidelobe error propagation for sky model uncertainties – and not the uncertainties in the actual sidelobe response function itself, which are very likely to dominate the total sidelobe error budget. Future works should therefore aim to introduce parametric models for the sidelobe response itself, and sample (or at least marginalize) over the corresponding free parameters, as these are typically among the most important unknowns for many experiments.

Finally, we note that future CMB experiments such as LiteBIRD, targeting low B-mode limits, may need to consider more complex ways of handling sidelobes and beams. The ultimate solution in this respect is 4π beam convolution for every single timestep, which could be achieved using a similar framework to the approach discussed here. This would remove the sidelobes as a nuisance signal from the data model of Eq. (13) and instead incorporate them directly into the beam term, B_pp′, j. This approach should be feasible for a relatively low-resolution experiment such as LiteBIRD and will be further investigated in a subsequent study.

Acknowledgments

We thank Prof. Pedro Ferreira and Dr. Charles Lawrence for useful suggestions, comments and discussions. We also thank the entire Planck and WMAP teams for invaluable support and discussions, and for their dedicated efforts through several decades without which this work would not be possible. The current work has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement numbers 776282 (COMPET-4; BEYONDPLANCK), 772253 (ERC; BITS2COSMOLOGY), and 819478 (ERC; COSMOGLOBE). In addition, the collaboration acknowledges support from ESA; ASI and INAF (Italy); NASA and DoE (USA); Tekes, Academy of Finland (grant no. 295113), CSC, and Magnus Ehrnrooth foundation (Finland); RCN (Norway; grant nos. 263011, 274990); and PRACE (EU).

References

All Figures

	Fig. 1. Map level difference in temperature of the new SHT convolution algorithm compared to the old Conviqt approach, observed using the identical pointing of the first year of the 30 GHz Planck detectors. The differences are at the level of machine precision, indicating full agreement between the two algorithms.
In the text

	Fig. 2. Runtime comparison between the libconviqt approach and the new spin-SHT approach for the convolution of an elliptical Gaussian with a set of random sky al, ms. This work ties or outperforms the previous approach for all values of lmax from 256 to 8192 for both mmax values shown. We note the log scale on the y-axis.
In the text

	Fig. 3. Maps of the sidelobe response on the sky from a representative detector at (left to right) 30 GHz, 44 GHz, and 70 GHz. The beam orientation is such that the main beam is pointed directly at the north pole in these maps. The intensities are normalized such that the main beams have unit power at l = 0.
In the text

	Fig. 4. Maps of the sidelobes convolved with the sky at each of the three LFI frequencies. From top to bottom: 30 GHz, 44 GHz and 70 GHz. Left column: unpolarized sky signal, central column: Q polarization and right column: U. Note the difference in the colour scales required to see the same level of detail in all three channels.
In the text

	Fig. 5. Sidelobe rms maps at each of the three LFI frequencies. From top to bottom: 30 GHz, 44 GHz, and 70 GHz. Left column: unpolarized sky signal, central column: Q polarization and right column: U. Again note the different colour scales.
In the text

	Fig. 6. Frequency map difference plots smoothed to one degree at (top to bottom) 30, 44, and 70 GHz, comparing two pipeline executions with the same seed, one of which has no sidelobe correction.
In the text

	Fig. 7. Component map difference plots for (top to bottom) CMB, synchrotron, AME and freefree emission, comparing two pipeline executions with the same seed, one of which has no sidelobe correction.
In the text

	Fig. 8. Difference in correlated noise, projected into the map domain at 30 GHz comparing two pipeline executions with the same seed, one of which has no sidelobe correction.
In the text