A Simulation-Based Method for Correcting Mode Coupling in CMB Angular Power Spectra

Modern CMB analysis pipelines regularly employ complex time-domain filters, beam models, masking, and other techniques during the production of sky maps and their corresponding angular power spectra. However, these processes can generate couplings between multipoles from the same spectrum and from different spectra, in addition to the typical power attenuation. Within the context of pseudo-$C_\ell$ based, MASTER-style analyses, the net effect of the time-domain filtering is commonly approximated by a multiplicative transfer function, $F_{\ell}$, that can fail to capture mode mixing and is dependent on the spectrum of the signal. To address these shortcomings, we have developed a simulation-based spectral correction approach that constructs a two-dimensional transfer matrix, $J_{\ell\ell'}$, which contains information about mode mixing in addition to mode attenuation. We demonstrate the application of this approach on data from the first flight of the SPIDER balloon-borne CMB experiment.

Modern cosmic microwave background (CMB) analysis pipelines regularly employ complex time-domain filters, beam models, masking, and other techniques during the production of sky maps and their corresponding angular power spectra. However, these processes can generate couplings between multipoles from the same spectrum and from different spectra, in addition to the typical power attenuation. Within the context of pseudo-C based, MASTER-style analyses, the net effect of the time-domain filtering is commonly approximated by a multiplicative transfer function, F , that can fail to capture mode mixing and is dependent on the spectrum of the signal. To address these shortcomings, we have developed a simulation-based spectral correction approach that constructs a two-dimensional transfer matrix, J , which contains information about mode mixing in addition to mode attenuation. We demonstrate the application of this approach on data from the first flight of the SPIDER balloon-borne CMB experiment.

INTRODUCTION
Producing well-characterized maps of the cosmic microwave background (CMB) from time-ordered data requires accurately accounting for the impact of instrumental effects and any signal processing on the underlying astrophysical signal (e.g., Jarosik et al. 2007;Planck Collaboration et al. 2016). These processing techniques typically also have a nontrivial impact on the fidelity of the cosmological signal in ways that are spatially anisotropic and inhomogeneous. These effects must be precisely and accurately characterized in order to avoid biasing the estimation of the CMB angular power spectra, and therefore the inference of cosmological parameters.
The impact of timestream processing and instrument response is commonly approximated in the harmonic-domain with a filter window function unique to the instrument, derived from the analysis of large ensembles of signal simulations. In its simplest formulation, the processing pipeline is modeled as a power attenuation mechanism in each multipole, i.e., the filter window function is a transfer functiona simple one-dimensional vector of ratios of output to input power.
This paper addresses the shortcomings of the onedimensional model and proposes an alternative approach through the construction of a two-dimensional transfer matrix. This simulation-based spectral correction approach takes into account both the mode mixing and attenuation from instrumental and data processing effects including beams, filtering, and masking. Section 2 of this paper describes the theoretical motivation and compares the transfer matrix to the one-dimensional transfer function approach. A concrete example is presented in Section 3 using data from the first flight of SPIDER, a balloon-borne telescope designed to measure CMB polarization on roughly degree angular scales (SPIDER Collaboration 2022). This section explores different techniques for constructing the transfer matrix to reduce the computational demand including using binned power spectra and performing Fourier-space interpolation. Section 4 presents several comparisons of these techniques, including tests of signal recovery on spectra with different shapes. While all tested approaches were found to accurately recover a target spectrum identical to that used for the trans-fer matrix construction, the performance varied when applied to a different target spectrum. As discussed in Section 5, this has implications for increasingly sensitive CMB polarization measurements where the cosmological signals are heavily obscured by Galactic foregrounds. As the foreground power spectra are less well constrained and vary substantially between different sky regions, understanding the signal dependence of potential analysis techniques becomes even more important.

THEORETICAL DESCRIPTION
Maps of the CMB temperature (T ) and linear polarization (Q, U) anisotropies over the partial sky can be decomposed into linear combinations of spherical harmonics: where W (r) is the window representing the relative weights of the partial-sky mask. To avoid using spin-weighted spherical harmonic components, the ±2ã m coefficients are frequently expressed as a combination of scalar and fixed parity E and B components where the sign convention follows Zaldarriaga & Seljak (1997). For eachã m , the pseudo-power spectrumC is defined asC Also known as a pseudo-C (PCL), it is related to the angular power spectrum specified by the theory of primordial perturbations, C , via where K XX is a mode coupling kernel (or "mixing matrix") that accounts for the mixing within and between the observables X ∈ {T T, EE, BB, T E, EB, T B} due to the partial-sky mask, and the brackets · denote an ensemble average over infinite spectrum realizations. Because K XX is entirely determined by the chosen pixel weighting and the geometry of the cut sky (Hivon et al. 2002), in the absence of instrumental effects and noise, PCL estimators can use the (finite) set of measuredC to solve for the underlying power spectrum C . Examples of PCL estimators include MASTER (Hivon et al. 2002), NaMaster (Alonso Monge et al. 2019), and PolSpice (Chon et al. 2004). A major challenge in interpreting CMB data is that the instruments cannot directly probe the true sky signal; the incoming signals are inevitably altered by instrumental systematics and noise. Therefore, as described in Section 1, an experiment's raw datastreams must be processed to remove many different types of spurious signals. The act of observing and time-domain filtering both distort the signal estimate and present additional sources of mode coupling. Assuming that the coupling is homogeneous, isotropic, and linear, the impact of the experiment is captured by introducing another coupling matrix, F XX , such that In general, F XX is unique to the experiment and cannot be determined analytically.

The Filter Transfer Function
To reduce complexity, F XX is often approximated as a diagonal matrix whose entries represent the one-dimensional transfer function F XX , such that (here and hereafter, the superscripts X are suppressed for brevity). Using a set of pseudo-power spectraC derived from an ensemble of simulations of a known power spectrum C , the transfer function F can be estimated through an iterative process to avoid inverting K (Hivon et al. 2002;Dutcher et al. 2021). Note that F is frequently decomposed further into a filter component f and a beam component b 2 , such that F = f b 2 ; additional instrumental effects can be inserted in a similar fashion. This one-dimensional approximation implicitly assumes that each -mode remains independent throughout the entire filtering process. As long as this assumption holds, the mapping from input to output is one-to-one: F is simply the ratio of output to input power for each . However, in practice, we expect modes to become coupled with one another, where the mapping becomes many-to-one and the contributions from the coupled input modes become impossible to disentangle using F alone. The consequence is that changing the input power spectrum C also changes the outputC in some nontrivial way due to this many-to-one mapping. In other words, F is inextricably tied to the input used to compute it.
More concretely, the one-dimensional transfer function formulation conflates mode mixing with the in-mode filter gain. This introduces a sensitivity to the power spectrum of the simulated sky used to calibrate F ; the final spectrum is correct only if the simulated sky is statistically similar to the true sky (i.e., a Gaussian sky realization with the same power spectra).

The Multipole-Multipole Transfer Matrix
To address this shortcoming of the F formulation, we introduce a two-dimensional linear operator that encodes the (asymmetric) coupling between each pair: We refer to this coupling operator J as the "transfer matrix" because it directly relates the input of the true power spectrum to the output of the spectrum estimator. This relation holds as long as the filtering process is approximately linear. Treating the entire pipeline as a single operator avoids the diagonal approximation and ensures that all multipole-multipole couplings induced by time-domain and map-domain filtering are properly taken into account, i.e., we are not locked into a specific input spectrum. Note that these couplings extend to those between the six power spectra; both temperature-to-polarization (T -to-P) and E-to-B leakage are automatically included. Because standard PCL spectrum estimators provide the T T , EE, BB, T E, EB, and T B power spectra, any input mode can be readily related to an output mode from any of these six power spectra, resulting in 6 2 coupling matrices. We find it convenient to compile these individual matrices into a single 6 × 6 block matrix encapsulating every coupling between the six power spectra. Following the rules of matrix multiplication (Equation 8), we arrange these blocks horizontally according to input spectra and vertically according to output spectra (see Figure  1). The ordering of the six spectra does not matter as long as it is consistent between the two axes; likewise, the six C s must be concatenated in the same order. We choose the above ordering based on convenience.
While the approach presented here has similarities to that used by the BICEP/Keck Collaboration for their CMB polarization analyses, the implementation varies due to key differences in the observing strategies and analysis pipelines. As described in Ade et al. (2016), the simplicity of the BI-CEP/Keck observing strategy allows for the construction of an observing matrix that renders large simulation ensembles more computationally tractable. Having determined their filtering operations to be linear, the observing matrix is used to construct the transfer matrix J , which is then used similarly to reconstruct and interpret the on-sky power spectra. Because the observing matrix formulation is less generally applicable to experiments at other sites with different observing strategies, the remainder of this paper is dedicated to estimating J through other means.
To illustrate the application of the transfer matrix, we present results from the SPIDER balloon-borne telescope. During its first Antarctic long-duration balloon (LDB) flight in January 2015, SPIDER mapped 4.8 % of the sky with polarimeters operating at 95 and 150 GHz to constrain the Bmode power spectrum from primordial gravitational waves (SPIDER Collaboration 2022). An upcoming flight with additional 280 GHz receivers will provide improved characterization of the polarized Galactic dust foregrounds (Shaw et al. 2020).
SPIDER's processing pipeline is described more fully in SPIDER Collaboration (2022) and SPIDER Collaboration (2022, in prep.); here we briefly highlight the most relevant steps and note that they are sufficiently linear to allow usage of Equation 8. Once features such as cosmic-ray hits, payload transmitter signals, and thermal transients have been removed from the raw detector timestreams, they are filtered to reduce quasi-stationary noise. Null test performance was used to identify the weakest filter that sufficiently removed quasi-stationary noise, resulting in a fifth-order polynomial fit to each detector's data as a function of azimuth angle between scan turnarounds. The impact of scanning, filtering, and flagging are determined by applying the entire processing pipeline to an ensemble of time-domain signal simulations in a procedure known as re-observation. The reobserved timestreams are produced at the full data sample rate without any downsampling or binning. Unlike the BI-CEP/Keck experiment, SPIDER's observing strategy does not allow for the creation of an observing matrix, so obtaining the transfer matrix J requires producing an appropriate set of re-observed CMB maps.
SPIDER's measured and simulated timestreams are converted into two-dimensional maps of the sky with a binned mapmaker. This approach assembles the detector data into spatial pixels based on the telescope pointing and polarization sensitivity as described further in SPIDER Collaboration (2022). The computational simplicity of this method is crucial for enabling the generation of large simulation ensembles, including those used in this work.
The main SPIDER cosmological results presented in SPI-DER Collaboration (2022) use two complementary pipelines for power spectrum estimation. Here we describe results only from the Noise Simulation Independent (NSI) pipeline, while the other pipeline is presented in Gambrel et al. (2021). The NSI pipeline is similar to Xspect (Tristram et al. 2005) and Xpol (Tristram 2006), as well as to those used by SPT (Lueker et al. 2010), CLASS (Padilla et al. 2020), and SPI-DER's circular polarization analysis (Nagy et al. 2017). It decomposes the 2015 flight data into 14 maps by splitting the timestreams into interleaved 3-minute segments for each of the six receivers. This timescale was chosen to maximize the number of complete and well-conditioned maps of the sky region, taking into account SPIDER's scan rate and observing strategy. These maps are then co-added by frequency to produce 14 maps in each of the two SPIDER bands, 95 and 150 GHz. The cross-power spectra for all pairs of maps (neglecting the auto-spectra) are produced with PolSpice (Chon et al. 2004), and the distribution of the cross-spectra allows for an empirical measurement of the uncertainty without reliance on an instrumental noise model. In SPIDER's NSI pipeline, these cross-spectra are used to estimate cosmological parameters by comparing them to theoretical models that depend on the parameters of interest. It is therefore SPI-DER's power spectra, rather than the maps, that we correct for the beam and filtering effects as described in the rest of this work. Since the PolSpice estimator is linear, the assumption of linearity in the definition of the transfer matrix (Equation 8) is satisfied. Although the transfer matrix can be used to correct the power spectra for temperature-to-polarization leakage, we perform this correction on the maps based on simulations of Planck temperature-only maps (Planck Collaboration et al. 2020). This avoids the sample variance from the simulation ensemble, which is relatively large due to the amplitude of the temperature signals, by using existing measurements of the temperature anisotropies in SPIDER's observing region.
As J is merely the linear operator linking the input and output power spectra (Equation 8), its computation is conceptually straightforward: pass a (simulated) CMB map with a known spectrum through the re-observation pipeline, then compute the power spectra of the output map. But to capture the asymmetric mode-mode coupling between each pair, this process must be performed on each mode individually. Thus we begin with a set of unit δ-functions -one at each of interest -and simulate a set of CMB maps using the synfast HEALPix utility 1 (Gorski et al. 2005). Upon obtaining the re-observed maps, we organize their power spectra into columns, arranged by , to form J . The whole process is illustrated schematically in Figure 1.
To reduce the effects of sample variance from the CMB map realization and partial-sky spectrum estimation procedures, we repeat the above process 100 times and average the 100 resulting matrices to obtain the final J . However, the computation cost is intensive: with 100 map realizations of a δ-function in each of the three CMB modes (T T , EE, BB) that must be re-observed by each of SPIDER's six receivers, we require 1800 simulated maps per ∈ [2, 250 + ∆ b ] (a ∆ b "buffer" is necessary to account for couplings of higher multipoles with lower multipoles). To minimize the computational burden, we use full-mission SPIDER maps instead of the 14 NSI pipeline submaps because this was demonstrated to have a negligible effect on SPIDER's signal recovery in simulations. Nevertheless, re-observing the entire set of required maps would take about 1000 core-years on the Niagara supercomputing cluster (Loken et al. 2010;Ponce et al. 2019), which is computationally infeasible given SPI-DER Collaboration resources. Since the fidelity of the computation needs to be balanced against the time required, two practical options are investigated in the following sections: 1) use input spectra with multipole bins rather than evaluating each multipole individually, or 2) compute the transfer matrix for only a sample of input s and interpolate between them.

Bin-Bin Transfer Matrix
CMB angular power spectra are typically binned into multipole ranges of width ∆ in order to increase the signalto-noise ratio and mitigate the correlative effects of adjacent multipole leakage. As a result, we consider the case where the input power spectra are also binned into the same multipole ranges, thereby reducing the computational demands of the calculation. By analogy with Equation 8, we define a bin-bin transfer matrix J bb that describes the average power in output bin b given the average power in each of the input bins b : where P b is a binning operator. At first glance, the δ-function inputs to J might be replaced by unit-boxcars for J bb , but, as it turns out, that may not be the best choice. Unlike the case for the full J , the binning introduces a dependence on the shape of the input spectrum used to estimate the transfer matrix (see the Appendix). Consequently, the choice of input spectrum should be as close as possible to the anticipated output spectrum, as is further discussed in Section 4. Because our application is toward degree-scale CMB B-modes, we use a Λ cold dark matter (CDM) + (r = 0.03) spectrum, provided by CAMB (Lewis et al. 2000), as our base model. Windowing this model spectrum at each bin of interest -i.e., multiplying it by unit-boxcars centred on each bin -provides the set of input spectra needed to construct the bin-bin transfer matrix J bb .
Because the analysis of SPIDER data uses only 12 bins in the range 2 ≤ ≤ 300, the computational demand is decreased by a factor of ∼30; months of cluster time is reduced to days. The price is a dependence on an input model, but unlike F , the ability of J bb to capture couplings between multipole bins provides a more general framework.

Fourier-space Interpolation
An alternative method for reducing computation time is to compute the transfer matrix J only at select multipoles (i.e., columns) and interpolate in between. This requires the matrix to be smooth in the -range of interest, but, as we shall see, this assumption is often well justified in practice.
Because J is a linear operator (Equation 8), the output from a δ-function input can be interpreted as the impulse response of a linear filter. As shown in the left panel of Figure 2, SPIDER's transfer components (T T → T T , EE → EE, etc.) have impulse responses resembling displaced sinc functions. This allows us to recast the problem from the domain to its dual frequency-like domain (hereafter simply the "frequency" domain), where the interpolation turns out to be considerably easier. For our application, we use the sign and normalization conventions established by NumPy's fft module; 2 for a given input , the discrete Fourier transform of its impulse response x -the frequency response -is defined by: We use N = 801 points (0 ≤ ≤ 800) in our computations.
The magnitudes |X f | and phases φ = ∠X f of the resulting transforms are shown in the right panel of Figure 2. In general, the properties of the responses are unique to each experiment and therefore often determined empirically. SPIDER's response corresponds to a linear-phase low-pass filter: the magnitude rolls off until it hits the noise floor of the filter at f = f c , while the phase φ is nearly linear up to that same frequency. The slope of the linear-phase component is proportional to the shift (in ) of the impulse, namely, m = φ/ f = −2π when f is expressed in normalized frequency units (cycles per multipole). The magnitude of the noise floor encodes the level of background noise at the input multipole . This simple structure in the frequency domain permits a straightforward interpolation. The magnitudes, in particular, can be interpolated directly at all f , as can the unwrapped phases at f < f c . For simplicity, we interpolate piecewiselinearly. In the noise region f > f c , the wrapped phases exhibit a curious behaviour: they settle onto one of two possible values separated by π. Thus, rather than interpolating, we implement a nearest-neighbour scheme for the phases at f > f c . Re-obs. theory spectrum Figure 4. Although anafast does not automatically correct for E-to-B leakage, a transfer matrix built using this utility can still capture such mode mixing induced by SPIDER's filter and mask. When the matrix is applied to a ΛCDM + (r = 0.03) theory spectrum, the result can predict the output from re-observing a map with the same theory spectrum (blue points; error bars indicate the error on the mean of the 500 simulations). However, replicating the full effects of the re-observation procedure in the higher bins requires computing a larger matrix than anticipated due to the wide tails of the EE → BB response. Namely, a matrix with max = 300 (dotted line) underestimates the power in the higher bins; this is alleviated somewhat by increasing to max = 370 (solid line), but that is still not quite enough. As a result, we proceed with a matrix constructed using the PolSpice estimator, as its built-in E-to-B leakage correction algorithm avoids the need to generate larger matrices, thereby reducing the computational demand.  Figure 6. Average ratios of re-observed 150 GHz spectra to those obtained from applying the transfer function and transfer matrices to fiducial ΛCDM + (r = 0.03) and foreground input spectra (500 and 300 realizations, respectively). A ratio of 1 represents the ideal performance, signifying a perfect replication of SPIDER's filtering process. As the transfer function F was derived using the re-observed ΛCDM result, applying it to that same spectrum again returns 1 by definition (left panels). Error bars indicate the error on the mean of the simulation ensembles.
One might be tempted to filter out the noise at f > f c completely. This requires care: any operation in the frequency domain carries consequences into the domain. For example, a simple low-pass filter with a rigid cutoff is equivalent to convolution with a sinc function in -space; this operation will pollute the impulse response with heavy ringing, and should be avoided. Because good filter design merits a study of its own, we do not carry out any filtering of this sort.
To complete the process, the interpolated values are inverse Fourier transformed back into -space to form our interpolated J interp (Figure 3). The choice of nodes for the interpolation, i.e., which columns of J to compute rather than interpolate, reduces to a balance between fidelity and computation time. Naturally we will want to concentrate our efforts to regions where J changes rapidly in order to obtain a good interpolation. For us, that is the low-regime. We thus compute the columns at = {5, 10, 15} in addition to SPIDER's bandcenters.
We note that this particular choice of interpolation method is motivated by SPIDER's sinc-function-like impulse responses. However, any sufficiently smooth J can be interpolated with similar techniques.

Choice of Spectrum Estimator
The matrix-making process encodes the spectrum estimator into the matrix -that is, J depends on the choice of estimator by design (Equation 8). Nonetheless, this choice can be trivially changed because the spectrum estimation and reobservation steps are completely independent of each other.
The anafast utility (provided in the healpy package) is a simple spectrum estimator that directly computes the power spectra from a set of input maps using Equation 4. In general, the interpolated matrix constructed using anafast provides very good agreement with the expected results. However, we find that the "buffer" needs to be quite large (∆ b > 120) in order to properly account for E-to-B leakage (Figure 4). To avoid the additional computational demand, in the following work, we proceed using the PolSpice estimator, which includes as part of its algorithm a correction for E-to-B leakage built in. This reduces the "buffer" to a much more manageable ∆ b = 50.
We use the same PolSpice parameters as the SPIDER B-mode analysis for this work, including a cosine apodization with σ = 10°, θ max = 50°, and symmetric_cl (SPIDER Collaboration 2022).  Figure 7. Solutions to the linear system J interp x = J interp C theory for a ΛCDM EE spectrum (left) and a dust foreground spectrum (right) computed by the truncated singular value decomposition (TSVD) method with 50 singular values and the biconjugate gradient stabilized (BiCGSTAB) method. Near-optimal results can be obtained when using the TSVD matrix as a preconditioner and providing an initial guess close to C theory (details provided in the text). Though not shown here, T T and BB spectra return similar results.

Cross-spectra: T E, EB, and T B
An outstanding problem in our discussion so far is the treatment of the cross-spectra T E, EB, and T B. While the inter-spectral couplings from auto-to cross-spectra (T T → T E, etc.) are implicitly accounted for by the spectrum estimator (which returns all six power spectra), the method described here cannot be used to compute the matrix components with a cross-spectrum as input. To wit, it is impossible to generate maps of Gaussian-distributed a m s with a nonzero correlation (the input mode) in T E, EB, or T B while the T T , EE, and BB spectra are imposed to be zero. Consequently, we approximate the diagonal cross-spectra transfer components (T E → T E, etc.) as the geometric mean of the related auto-spectra components, J T E,T E = (|J T T,T T ||J EE,EE |) 1/2 , etc., and treat the off-diagonal components (T E → EB, etc.) as negligible. Although this is not ideal, it was found not to bias simulation results for our application. An alternative method could be to construct the transfer matrix using the harmonic modes, rather than the power spectra, of the maps. As an example, see Choi et al. (2020), who postulate a transfer matrix built from (two-dimensional) Fourier transforms of T , E, and B maps.

COMPARISON OF METHODS
We quantify the effectiveness of the two methods described above (binning and interpolating) by two metrics: their ability to replicate the end-to-end SPIDER pipeline, and to recover the original signal. In the context of Equation 8, this corresponds to producing the correctC when given C , and vice versa.

Case 1: Solve forC , Given C
We demonstrate the first of these by applying our computed transfer function/matrices directly to a theory spectrum. This result is then compared to the output spectrum generated by passing a set of map realizations of the same input theory spectrum through the SPIDER re-observation procedure. To ensure broad applicability, we use two qualitatively different theory spectra as input: a ΛCDM + (r = 0.03) spectrum (from CAMB; 500 realizations) and a power-law dust foreground spectrum (300 realizations). These two spectra are chosen for their contrasting behaviour ( Figure 5).
As discussed in Section 2.1, the input-dependence of F makes it suitable only when the target spectrum does not deviate far from the input spectrum that was used to construct F . This is demonstrated in Figure 6. Having derived our F using a fiducial ΛCDM + (r = 0.03) spectrum as input, applying it to the same spectrum again outputs a trivial result (left panels). But when the target is the foreground spectrum (right panels), large deviations occur wherever mode-mode coupling dominates over in-mode gain -i.e., wherever the transfer matrix is least diagonal-like. In our case, this occurs at low s (Figure 3). The simple power attenuation model of F is insufficient in capturing SPIDER's filtering scheme in this regime.
The same kind of input-dependence is seen in the binned transfer matrix J bb (see the Appendix). Like F , the J bb constructed from a fiducial ΛCDM+(r = 0.03) spectrum is highly effective when the target spectrum is also a ΛCDM + (r = 0.03) spectrum, and diverges when the target is the fore-ground dust spectrum ( Figure 6). But this time, all foreground bins up to = 100 are noticeably afflicted. This is due to the loss of granularity from binning the spectrum and the matrix, as we had used a bin width of ∆ = 25 to match SPIDER's B-mode analysis. In the limit where the bin width is 1 -a single multipole -we recover the full, theoretically ideal transfer matrix J .
The interpolated transfer matrix J interp thus functions as a compromise allowing us to have a bin width of 1 without being computationally intractable. Despite the simple piecewise-linear interpolation approximation used, it reliably recovers both target spectra.

Case 2: Solve for C , GivenC
To recover C fromC , the linear system defined by Equation 8 must be solved numerically; directly inverting a linear system is not recommended except in special cases. One such case is the transfer function F -computing its inverse is trivial since it is a diagonal matrix (a one-dimensional function). Likewise, SPIDER's binned transfer matrix J bb is also easily invertible because it is diagonally dominant and measures only 72 × 72 (12 bins in each of the six correlation spectra); its condition number is very close to 1. The main challenge, therefore, is solving the system for the interpolated transfer matrix J interp , whose condition number is O(10 17 ).
The large condition number of J interp -and also of J in general -means that it is highly susceptible to numerical noise and cannot be solved using elementary linear algebraic techniques. Numerous algorithms of varying complexities exist to solve such ill-conditioned linear systems; the study of these techniques is beyond the scope of this paper. As a demonstration, below we focus on the truncated singular value decomposition (TSVD) method and the biconjugate gradient stabilized (BiCGSTAB) method (van der Vorst 1992). To test the efficacy of these two techniques, we computeC = J interp C theory , solve the system J interp x =C , and compare the solution x with the original C theory . The results are shown in Figure 7. We use the same test spectra as the previous section ( Figure 5) as C theory .
TSVD approximates the inverse matrix by rank-reducing the factorized matrices such that only the N largest singular values remain. Because singular values represent magnitudes along orthogonal axes, this truncation removes the low-amplitude components most susceptible to numerical instabilities. For our J interp , we truncate the singular values at the quasi-arbitrary cutoff of 10 −3 , which leaves the top 50 singular values; below this cutoff, the singular vectors are inundated with numerical noise. As shown in Figure 7, TSVD with N = 50 does well for a ΛCDM spectrum but remains noisy for a dust power law.
On the other hand, BiCGSTAB is an iterative algorithm. While we find that its standard configuration can recover the C theory to an acceptable tolerance, its convergence can be aided by a suitable preconditioner matrix and initial guess. For convenience, we use the truncated matrix from our TSVD tests above as the preconditioner, and, to replicate real conditions, we use a noisy version of C theory (±5 % noise) as the initial guess. This produces near-optimal results ( Figure  7) in fewer than 10 iterations.

CONCLUSION AND DISCUSSION
The transfer function F is an Ansatz introduced by Hivon et al. (2002) to represent the effects of datastream filtering in the process of CMB mapmaking. It was noted at the time that F cannot fully capture such filtering, on account of inconsistencies with the statistics governing the power spectrum C (in particular, time-domain filtering violates the assumption of isotropy), but the authors showed that such inconsistencies had a negligible impact on their results. With filtering and analysis techniques having grown far more complex as CMB experiments have become more sensitive in the two decades since, and considering the crucial role F fills in recovering cosmological signals, this Ansatz merits a more thorough re-examination. This is particularly important for measurements of CMB polarization in the presence of large foreground components, where there is significant uncertainty in the shape of the power spectrum for a given sky region, and where foreground residuals may be a concern in component-separated maps.
In this study, we present a different approach, the transfer matrix J , that accounts for a key element missing in F : the multipole-multipole coupling induced by the instrument response and timestream filtering. These couplings are not restricted to those within the same spectrum; also included are inter-spectral coupling components such as T -to-P and E-to-B leakage. Because the full matrix is computationally infeasible for most modern experiments, we additionally explore two approximations: a binned transfer matrix J bb and an interpolated transfer matrix J interp . Along with F , we compare their performance by applying them to the sky-to-spectrum pipeline from SPIDER's first flight.
When the target is a ΛCDM-like spectrum, we find that all three approaches are effective at replicating the sky-tospectrum pipeline. However, when the target is switched to a power-law dust foreground model, F and J bb -both being dependent on the fiducial model used to construct them -falter at large scales, and only J interp is able to reliably recover the underlying signal. However, like the full transfer matrix J itself, J interp is ill-conditioned and care needs to be taken in solving the linear system. ΛCDM-like T T , EE, and BB spectra are easily recovered using a TSVD approach, but a dust foreground-like spectrum requires more sophisticated algorithms like the BiCGSTAB method to remain numerically stable.
Nonetheless, because SPIDER employs a map-based method of foreground removal, we expect the remaining signal to deviate, at most, only slightly from ΛCDM. As long as this is the case, the error from binning remains small (Appendix). We therefore chose to use the binned transfer matrix J bb in the NSI pipeline for SPIDER's main result (SPI-DER Collaboration 2022), as it reliably recovers a ΛCDMlike spectrum ( Figure 6) and is easily invertible. However, we caution that J bb should only be pursued when the binning scheme used in analysis has been established, as its input bins must match the output bins; any changes to the binning scheme must be accompanied by a re-computation of J bb . On the other hand, J interp is only as accurate as the interpolation allows, but is attractive for the capacity to increase its resolution as desired. Additionally, J interp may be most advantageous for large angular scale foreground analyses, where it performs significantly better than J bb .
Of course, there is no need to adhere strictly to one technique. If computational resources are limited, one may employ a hybrid approach and compute the transfer matrix only in regimes dominated by mode-mode coupling, i.e., where the transfer matrix has significant off-diagonal components. For example, one may choose to commit to fully computing every column of the lowest bin, where the discrepancy between the input and output spectra is greatest (see Figure 6), and use a less expensive strategy elsewhere. We also leave open the exploration of different interpolation schemes, as the implementation described in this paper has not been found to significantly affect SPIDER's B-mode results. Alternatively, the simple structure of the transfer matrix in Fourier space also suggests, in place of interpolation, a model with a polynomial fit in magnitude and a linear fit in phase. We note that the biggest bottleneck is often computational time for the production of re-observed maps; as long as this portion of the pipeline is fixed, downstream tasks such as masking and choice of spectrum estimator are trivial to change.
As CMB experiments become capable of more precise measurements, improving the accuracy of the signal recovery process becomes increasingly important. The transfer matrix decouples the estimation of power spectra -and, by extension, the cosmological parameters -from assumptions about the shape of the underlying sky signal. This is particularly useful when the cosmological signals are heavily obscured by foreground Galactic dust, and is thus relevant for upcoming CMB instruments searching for primordial Bmodes (e.g., Abazajian et al. 2019;Ade et al. 2019;Hazumi et al. 2020;Shaw et al. 2020) or measuring the E-modes on the largest angular scales (Allison et al. 2015;Watts et al. 2018;Hazumi et al. 2020). As foreground component separation becomes an increasingly important part of CMB analyses, future experiments can benefit from the development of new signal-independent analysis techniques to improve the power spectrum reconstruction.

APPENDIX A. RELATION BETWEEN BIN-BIN AND MULTIPOLE-MULTIPOLE TRANSFER MATRICES
The power spectrum binning procedure can be described by an operator P b that determines the relative weighting of each multipole within a bin b: Likewise, we can bin the pseudo-power spectrum and writẽ as a consequence of Equation 8. Combining Equations 9 and A2, b J bb P b C = P b J C .
From this equation, we can identify the binned transfer matrix J bb as taking the form Note that the inner sum in the numerator is confined to the input bin b -recall each column of J bb is the output of passing a single input bin through the re-observation pipeline (in analogy with the unit δ-functions used to construct J ). Taking the sum over all input bins b , as on the left-hand side of Equation A3, recovers the sum over all input multipoles , as on the right-hand side. For SPIDER, we employ a uniform weighting, which reduces Equation A4 to The quantity in parentheses in Equation A6 can be taken to be J b , the proportion of the power in input bin b that contributes to the power in output multipole . It is equal to the weighted average of all of the J elements that lie within input bin b , with the weights given by some input model spectrum C . This highlights an important point: while J is a function of the data analysis pipeline only, J bb is additionally dependent on the input spectrum used to create it. We can understand this point by noting how J and J bb are each applied to a target power spectrum c : If we let δ represent the deviation of c from C (i.e., c = C + δ ; not to be confused with a δ-function centred at ), then Equations A7a and A7b can be understood intuitively: whereasc b represents the true amount of filtered power in bin b,ĉ b is an estimate that differs from the true value by In other words, the error from using J bb is the difference between the average of transformed excess power in each multipole and transform of average excess power in each bin b. Note that as δ → 0, the deviation δ b → 0 as expected, i.e., J bb returns the same result as J when applied to the input spectrum C used to create it. Consequently, the choice of input spectrum should be as close as possible to the target spectrum c to be applied, such that δ b is minimized.