# High precision modeling of polarized signals: moment expansion method generalized to spin-2 fields

Léo Vacher, Jens Chluba, Jonathan Aumont, Aditya Rotti, Ludovic Montier

► To cite this version:

Léo Vacher, Jens Chluba, Jonathan Aumont, Aditya Rotti, Ludovic Montier. High precision modeling of polarized signals: moment expansion method generalized to spin-2 fields. Astronomy and Astrophysics - A&A, 2023, 669, pp.A5. 10.1051/0004-6361/202243913. [hal-03669893](https://hal.science/hal-03669893)

HAL Id: hal-03669893

https://hal.science/hal-03669893

Submitted on 13 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# High precision modeling of polarized signals: Moment expansion method generalized to spin-2 fields

L. Vacher<fnref n="1" />, J. Chluba<fnref n="2" />, J. Aumont<fnref n="1" />, A. Rotti<fnref n="2" />, and L. Montier<fnref n="1" />

Received 30 April 2022 / Accepted 29 October 2022

## ABSTRACT

The modeling and removal of foregrounds poses a major challenge to searches for signals from inflation using the cosmic microwave background (CMB). In particular, the modeling of CMB foregrounds including various spatial averaging effects introduces multiple complications that will have to be accounted for in upcoming analyses. In this work, we introduce the generalization of the intensity moment expansion to the spin-2 field of linear polarization: the spin-moment expansion. Within this framework, moments become spin-2 objects that are directly related to the underlying spectral parameter and polarization angle distribution functions. In obtaining the required expressions for the polarization modeling, we highlight the similarities and differences with the intensity moment methods. A spinor rotation in the complex plane with frequency naturally arises from the first order moment when the signal contains both spectral parameter and polarization angle variations. Additional dependencies are introduced at higher order, and we demonstrate how these can be accounted with several illustrative examples. Our new modeling of the polarized signals reveals to be a powerful tool to model the frequency dependence of the polarization angle. As such, it can be immediately applied to numerous astrophysical situations.

**Key words.** polarization – cosmic background radiation – cosmology: observations – dust, extinction – ISM: general

## 1. Introduction

A significant international effort has been undertaken to deploy large-scale surveys of the cosmic microwave background (CMB) signal in the present, near, and far future. Multiple telescopes, such as ACT (Aiola et al. 2020), SPT (Sayre et al. 2020), The Simons Observatory (The Simons Observatory collaboration 2019), and CMB-S4 (CMB-S4 Collaboration 2019), are or will be observing large portions of the sky from the ground. Similarly, from space we are eagerly awaiting LiteBIRD (LiteBIRD Collaboration 2022), and in the future possibly even more ambitious CMB imagers (PICO Collaboration 2019; Delabrouille et al. 2021). In addition, we can hope for a CMB spectrometer such as PIXIE (Kogut et al. 2011) to target CMB spectral distortions (Chluba et al. 2021). The scientific targets are manyfold and of prime importance to cosmology, astrophysics, and high-energy physics.

Ever-increasing instrumental sensitivities imply that one also becomes sensitive to faint and complex effects that need to be properly modeled. The fine characterization of polarized astrophysical signals thus becomes an increasingly complicated and important challenge to cosmological analyses. The stakes are twofold: first, we wish to reach a better understanding of the sources themselves and the complex physics at play in the emission. This is of importance for Galactic physics, physics of recombination at the last scattering surface with the primordial anisotropies of the cosmic microwave background (CMB) and their secondary sources from the physics of galaxy clusters with the Sunyaev-Zeldovich (SZ) effect (Carlstrom et al. 2002; Mroczkowski et al. 2019), or lensing (Stompor & Efstathiou 1999). Secondly, high-precision component separation is required to remove the diffuse foregrounds and access the faint

perturbations of the CMB polarized signal. Finding how to deal with them is an unavoidable step in addressing questions that concern the cosmic history and high-energy physics through inferring key cosmological parameters as the reionization depth $\tau$ (Wise 2019) or the scalar-to-tensor ratio $r$, which probes inflation (Brout et al. 1978; Starobinsky 1980; Guth 1981). The philosophy of the present work is that one cannot properly deal with the second point without a detailed understanding of the polarized Galactic emission, deeply rooted in the physics.

It is known that the emission properties of the diffuse interstellar medium (ISM) change across the Galaxy on small and large scales. This assertion is well motivated from the Galactic magnetic field physics and the distribution of dust grain shapes and composition (Ferrière 2001; Jaffe et al. 2013), which is supported by numerous observations (Ysard et al. 2013; Hutton et al. 2015; Schlafly et al. 2016; Planck Collaboration Int. L 2017; Planck Collaboration XI 2020; Pelgrims et al. 2021). In observational conditions, averages over Galactic voxels with different spectral parameters are then unavoidable. They occur in several situations: along the line of sight inside the Galaxy, which cannot be reduced or avoided with instrumental considerations; between lines of sight inside the instrumental beam; and over patches of the sky when doing a spherical harmonic decomposition of the signal or averaging the data otherwise. Two related consequences follow immediately for the signal in intensity: if the fundamental spectral energy distributions (SEDs) are nonlinear, the average SED of the total signal differs from the canonical SED of the voxel. We refer to this phenomenon as SED distortions. The SED is distorted differently from one point of the sky to another, breaking the correlation between the maps at different frequency bands, leading to inaccurate extrapolations from one to another: we refer to this phenomenon as frequency decorrelation (see

<footnote n="1">¹ IRAP, Université de Toulouse, CNRS, CNES, UPS, Toulouse, France
e-mail: leo.vacher@irap.omp.eu</footnote>

<footnote n="2">² Jodrell Bank Centre for Astrophysics, Alan Turing Building, University of Manchester, Manchester M13 9PL, UK</footnote>
e.g., Tassis & Pavlidou 2015; Planck Collaboration Int. L 2017; Pelgrims et al. 2021).

To treat these averages in connection with CMB foregrounds, the (Taylor) moment expansion formalism was proposed (Chluba et al. 2017). A similar idea had been applied to the modeling of Sunyaev-Zeldovich signals, showing how spatial and frequency information can be nicely separated (Chluba et al. 2013). The moment formalism has proven to be very powerful when applied to component separation at the map level (Rotti & Chluba 2021; Remazeilles et al. 2016, 2021). A straightforward generalization to harmonic space and cross-frequency power-spectra domain has also proven to be useful (Mangilli et al. 2021; Azzoni et al. 2021; Vacher et al. 2022).

While the original formulation of the moment method was focused on the intensity, it was stressed that an extension to polarization can be readily obtained (Chluba et al. 2017), which is what we intend to do in the present work. Indeed, several applications already used the moment expansion method for polarized signals, treating the B-mode signal as an intensity (e.g., Remazeilles et al. 2021; Azzoni et al. 2021; Vacher et al. 2022). One can also find a similar approach in the Delta-map method (Ichiki et al. 2019), which used first order terms of the Q/U-intensity moments and already suggested a common treatment for the pair (Q, U).

In this work, we plan to rigorously derive and extend the moment expansion method to polarized signals. To do so, the Stokes parameters must be treated together, as the components of a single complex object. In the most general cases, extra subtleties come into play, which were not captured or discussed before. There are numerous advantages from thinking of linear polarization as a spin-2 quantity. As such they are not only described by a scalar quantity but also by an angle: the polarization angle. The averaging processes listed above will have one extra consequence for polarized signal: additionally to the spectral parameters, multiple angles will be mixed along and between lines of sight. We refer to this phenomenon as polarized mixing. In the presence of polarized mixing, the total signal will exhibit a frequency-dependent polarization angle. Being able to accurately model this frequency-dependent rotation from physically-motivated considerations represents a thorny challenge. In this work, we attempt to provide this extension to linearly polarized signal in a formal, natural, complete and self-consistent way. We pay particular attention to the formulation in terms of SED parameter distribution functions, which really is the origin of the name “moment expansion”. We subsequently see that this rewriting offers a powerful framework to grasp polarized mixing and its consequences.

After a review of the intensity moment formalism in Sect. 2, we discuss the nature of linear polarization and introduce the spin-moment formalism in a single line of sight in Sect. 3. In Sect. 4, we explore different example of sums of canonical SEDs along a line of sight, that are of astrophysical relevance. We study them both analytically and through a fitting procedure, demonstrating the ability of the spin-moment formalism to grasp distortions of the polarized SED. In Sect. 5, we generalize the formalism to deal with other kind of averaging effects: spherical harmonic transforms and instrumental effects. In Sect. 6, we discuss cases with extra complexity as Faraday rotation and more general voxel SEDs. Finally, we conclude in Sect. 7.

## 2. Intensity moment expansion

Before we discuss the generalization of the moment expansion for polarized light, we briefly recall the logical steps followed in

Chluba et al. (2017) to obtain the moment expansion in intensity. For now, we neglect beam averaging effects or expansions into spherical harmonic, but we cover these in Sect. 5.

We start by considering various voxels along a line of sight in the direction $\hat{n}$, which is described by an affine parameter s. Every voxel emits with an SED<fnref n="1" />:

$$I_v(A(s), \mathbf{p}(s)) = A(s) \hat{I}_v(\mathbf{p}(s)), \qquad (1)$$

where $\hat{I}_v$ is referred to as the fundamental SED with N spectral parameters, $\mathbf{p}(s) = \{p_1(s), p_2(s), \dots, p_N(s)\}$. The amplitude or weight parameter, $A(s)$, determines the relative contribution of each voxel to the total intensity<fnref n="2" />. The resulting total SED is given by an average along the line of sight, which we shall denote by $\langle \dots \rangle$. This average can be explicitly written in terms of an integral over the affine parameters, s, or, alternatively, as an integral of the intensity over the spectral parameter distribution function in the direction $\hat{n}$ (e.g., Chluba et al. 2017; Rotti & Chluba 2021):

$$\langle I_v(A, \mathbf{p}) \rangle = \int \frac{dA(s)}{ds} \hat{I}_v(\mathbf{p}(s)) ds \equiv \int \mathbb{P}(\mathbf{p}, \hat{\mathbf{n}}) \hat{I}_v(\mathbf{p}) d^N p. \qquad (2)$$

In the second definition, we introduced the distribution $\mathbb{P}(\mathbf{p}, \hat{\mathbf{n}})$ of the spectral parameters, $\mathbf{p}$, along the fixed line of sight $\hat{\mathbf{n}}$, with the relative weights absorbed into the distribution itself. We note that the distribution $\mathbb{P}(\mathbf{p}, \hat{\mathbf{n}})$ is not necessarily normalized to unity, as it determines the relative weight of each SED shape to the total. For convenience we shall define the average amplitude parameter as $\bar{A} = \langle A \rangle = \int [dA(s)/ds] ds = \int \mathbb{P}(\mathbf{p}, \hat{\mathbf{n}}) d^N p$.

In order to provide a perturbative model of the average SED, the spectral dependence $\hat{I}_v$ in Eq. (2) can be expanded into a Taylor series with respect to $\mathbf{p}$ around the pivot $\bar{\mathbf{p}}$ as:

$$\begin{aligned} \hat{I}_v(\mathbf{p}) = \hat{I}_v(\bar{\mathbf{p}}) &+ \sum_j (p_j - \bar{p}_j) \partial_{\bar{p}_j} \hat{I}_v(\bar{\mathbf{p}}) \\ &+ \frac{1}{2} \sum_{j,k} (p_j - \bar{p}_j)(p_k - \bar{p}_k) \partial_{\bar{p}_j} \partial_{\bar{p}_k} \hat{I}_v(\bar{\mathbf{p}}) \\ &+ \frac{1}{3!} \sum_{j,k,l} (p_j - \bar{p}_j)(p_k - \bar{p}_k)(p_l - \bar{p}_l) \partial_{\bar{p}_j} \partial_{\bar{p}_k} \partial_{\bar{p}_l} \hat{I}_v(\bar{\mathbf{p}}) \\ &+ \dots \end{aligned} \qquad (3)$$

Here, we used the shorthand notation $\partial_{\bar{p}_j} X(\bar{\mathbf{p}}) \equiv \partial X(\bar{\mathbf{p}})/\partial \bar{p}_j$. The pivot value $\bar{\mathbf{p}}$ around which the series is carried out can be fixed by asking for the first term of the expansion to vanish upon averaging: $\langle A \sum_j (p_j - \bar{p}_j) \rangle = 0$. This minimizes the required terms in the Taylor series and leads to:

$$\bar{\mathbf{p}} = \frac{\langle A \mathbf{p} \rangle}{\bar{A}} = \frac{\int \mathbb{P}(\mathbf{p}, \hat{\mathbf{n}}) \mathbf{p} d^N p}{\int \mathbb{P}(\mathbf{p}, \hat{\mathbf{n}}) d^N p}. \qquad (4)$$

<footnote n="1">¹ Hereafter, we use the shorthand notation for frequency dependent quantities $X_v \equiv X(v)$.</footnote>

<footnote n="2">² One simple example is the power law: $I_v(A, \beta) = A (v/v_0)^\beta$ where $\mathbf{p} = \{\beta\}$ has dimension one ($N = 1$) and $\hat{I}_v(\beta) = (v/v_0)^\beta$. In this, $v_0$ is arbitrarily defined, but in practical applications the choice is normally data-driven (e.g. motivated by the location of sensitive bands in experiments such as for Planck or WMAP) and A is the overall weight, with dimension that depends on the situation (see Appendix B).</footnote>
Next we introduce the moment coefficients of order $\alpha$

$$ \begin{aligned} \omega_{\alpha}^{p_j \dots p_l} &= \frac{\langle A(p_j - \bar{p}_j) \dots (p_l - \bar{p}_l) \rangle}{\bar{A}} \\ &= \frac{\int \mathbb{P}(\mathbf{p}, \hat{\mathbf{n}}) (p_j - \bar{p}_j) \dots (p_l - \bar{p}_l) \, d^N p}{\int \mathbb{P}(\mathbf{p}, \hat{\mathbf{n}}) \, d^N p}, \end{aligned} \quad (5) $$

with $\alpha$ being the number of parameters over which the average is done and the maximal order of the derivative associated with the moment coefficient. One can then write the total intensity as an expansion in terms of these moments:

$$ \begin{aligned} \langle I_v(A, \mathbf{p}) \rangle = I_v(\bar{A}, \bar{\mathbf{p}}) &+ \sum_j^N \omega_1^{p_j} \partial_{\bar{p}_j} I_v(\bar{A}, \bar{\mathbf{p}}) \\ &+ \frac{1}{2} \sum_{j,k}^N \omega_2^{p_j p_k} \partial_{\bar{p}_j} \partial_{\bar{p}_k} I_v(\bar{A}, \bar{\mathbf{p}}) \\ &+ \frac{1}{3!} \sum_{j,k,l}^N \omega_3^{p_j p_k p_l} \partial_{\bar{p}_j} \partial_{\bar{p}_k} \partial_{\bar{p}_l} I_v(\bar{A}, \bar{\mathbf{p}}) + \dots \end{aligned} \quad (6) $$

Here, all the $\omega_1^{p_i}$ are zero when using the value for $\bar{\mathbf{p}}$ as given by Eq. (4), while higher order moments capture the complexities added by line-of-sight averaging effects. In applications, the pivot value can be obtained using an iterative process by starting with a reasonable guess for $\bar{\mathbf{p}}$ and then correcting the solution by the values of $\omega_1^{p_i} \equiv \Delta p_i$ using $\bar{p}_j \rightarrow \bar{p}_j' = \bar{p}_j + \Delta p_j$. This iterative process assumes that the moments are perturbative and convergence can be obtained with a finite number of terms.

### 3. Moment expansion of polarized signals

In this section, we generalize the intensity moment expansion to polarization. We start by summarizing a few general aspects about how to describe polarized light and then highlight some of the important differences between intensity and polarization.

#### 3.1. General introduction to polarized SED

A polarized signal is fully described by the four real Stokes parameters $I_v, Q_v, U_v, V_v$, all of them being frequency dependent quantities. As before, $I_v$ describes the total (i.e., unpolarized + polarized) intensity, while the pair $(Q_v, U_v)$ and $V_v$ respectively quantify the linearly and the circularly polarized part of the photon field. Both $I_v$ and $V_v$ are scalar fields meaning that they are invariant quantities under transformations of the frame in which they are evaluated. As such, they can be described using the intensity moment expansion of Sect. 2. However a general treatment of polarized light including $V_v$ with the moment expansion could introduce extra subtleties which go beyond the scope of this work. Henceforth, we assume $V_v = 0^3$.

On the other hand, the $Q_v$ and $U_v$ are coordinate-dependent quantities transforming under frame rotations as the components of a spin-2 object. Therefore, they can be more naturally combined into a single spinor field $\mathcal{P}_v$:

$$ \mathcal{P}_v = Q_v + i U_v = P_v e^{2i\gamma_v}, \quad \text{with } i = \sqrt{-1}. \quad (7) $$

The spinor's modulus, $P_v$, is a real positive function called the linear polarization intensity and its argument defines the polarization angle $\gamma_v$:

$$ |\mathcal{P}_v| = P_v = \sqrt{Q_v^2 + U_v^2} \quad (8a) $$

$$ \arg(\mathcal{P}_v) = 2\gamma_v = \tan^{-1}(U_v/Q_v). \quad (8b) $$

As a spin-2 quantity, when the frame in which $Q_v$ and $U_v$ are defined (e.g., by modifying the directions of the polarizers) is rotated by a right handed rotation around the $\hat{\mathbf{n}}$ direction by an angle $\theta$, $\mathcal{P}_v$ transforms as Zaldarriaga & Seljak (1997):

$$ (\mathcal{P}_v)' = e^{-2i\theta} \mathcal{P}_v. \quad (9) $$

Note that, unless stated otherwise, henceforth we use calligraphic variables for complex quantities (e.g., $\mathcal{P}_v, W, \dots$) and italic font for real quantities (e.g., $P_v, A, Q_v, U_v, \omega, \dots$).

#### 3.2. Origin of frequency-dependent polarization angle from polarized mixing

The canonical SEDs usually considered in astrophysics (e.g. power laws, blackbodies, gray-bodies, modified blackbodies) generally assume a constant value for $\gamma$, independent of the frequency. This behavior is motivated by the existence of a preferred direction in the physical mechanisms at the origin of polarized emission such as magnetic fields and dust grain shape. A single emitting voxel in the Galaxy is thus expected to emit with a constant polarization angle as a function of frequency.

Adding voxels with the same polarization angle but varying spectral parameters simply leads to spectral complexity, very much like for intensity. Mixing varying polarization states with the same SED simply leads to a change in the direction of the total polarization, but no extra spectral complexity. However, if one mixes various SEDs with different polarization angles and various spectral parameters, the resulting polarized signal will inherit a distorted SED $P'_v \neq P_v$ and a frequency dependent $\gamma_v$.

$$ P_v(A_1, \mathbf{p}_1) e^{2i\gamma_1} + P_v(A_2, \mathbf{p}_2) e^{2i\gamma_2} + \dots = P'_v e^{2i\gamma_v}. \quad (10) $$

This consequence of polarized mixing is illustrated in Fig. 1 for a sum of two power laws $A_i(v/v_0)^{\beta_i} e^{2i\gamma_i}$ with $v_0 = 300$ GHz, $A_1 = 2$ Jy sr$^{-1}$, $A_2 = 1$ Jy sr$^{-1}$, $\beta_1 = 1.8$, $\beta_2 = 1.2$, $2\gamma_1 = 10^\circ$ and $2\gamma_2 = 80^\circ$. One clearly sees that the resulting spinor $\mathcal{P}_v$ rotates in the complex plane with frequency. Modeling both the distorted SED $P'$ and the frequency dependence of $\gamma_v$ in a physically-motivated fashion is nontrivial, but can be achieved when generalizing the moment expansion to polarization.

#### 3.3. Understanding the link to intensity moment expansion

To generalize the intensity moment expansion to polarization, we have to discuss how $Q_v$ and $U_v$ are obtained, and linked to intensity. To characterize the polarization state of the photon field, we measure the intensity with linear polarizers in four directions, $I_{v,\parallel}$ and $I_{v,\perp}$, which are orthogonal to each other, and $I_{v,\times}$ and $I_{v,\otimes}$, which are also orthogonal to each other but rotated by $45^\circ$ relative to the previous system. The total intensity $I_v$ (polarized + unpolarized), and Stokes $Q_v$ and $U_v$ are then given by

$$ (I_v^{\text{tot}})^2 = (I_v^{\text{unpol}})^2 + Q_v^2 + U_v^2 \quad (11a) $$

$$ Q_v = \frac{I_{v,\parallel} - I_{v,\perp}}{2}, \quad U_v = \frac{I_{v,\times} - I_{v,\otimes}}{2}. \quad (11b) $$

<footnote n="3">³ This is justified for CMB signals and component separation since classical physics in the primordial plasma is not expected to be source of any significant circular polarization (Montero-Camacho & Hirata 2018; Inomata & Kamionkowski 2019). Note however that a faint primordial V signal is expected in some models, see e.g. Hoseinpour et al. (2020).</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Illustration of the polarization angle spectral dependence from polarized mixing. Upper panel: polarization angle $\gamma_v$ as a function of frequency (in GHz) for the sum of two power laws with parameters $A_1 = 2$ Jy sr$^{-1}$, $A_2 = 1$ Jy sr$^{-1}$, $\beta_1 = 1.8$, $\beta_2 = 1.2$, $2\gamma_1 = 10^\circ$ and $2\gamma_2 = 80^\circ$. The color labels the frequency between 1 GHz (dark red) and 500 GHz (dark blue). Lower panel: polarization spinor $\mathcal{P}_v$ (in Jy sr$^{-1}$) in the complex $(Q, U)$-plane, in the same configuration and with the same color coding. The phase of the spinor is $2\gamma_v$. The length of the colored bars and the lines of constant radius represent values of $\log_{10}(P_v)$. The two black lines represent the two power laws $A_1(v/v_0)^{\beta_1}e^{2i\gamma_1}$ and $A_2(v/v_0)^{\beta_2}e^{2i\gamma_2}$ at 500 GHz.</figcaption></figure>

This shows that both $Q_v$ and $U_v$ describe differences between intensities and as such can have positive and negative contributions, depending on which polarizer response dominates. Thinking of each of the intensities $I_{v,\parallel}$, $I_{v,\perp}$, $I_{v,\times}$ and $I_{v,\otimes}$ as the cumulative signal from various emitters, means that one can create net polarization by (i) varying the number of emitters, that is to say the weight parameter $A$, and (ii) changing the spectra of the emitters in the different directions.

To give an example, let us assume that in all directions we have a simple gray-body SED, $I_v^{GB} = A B_v(T)$, where $B_v(T)$ is a blackbody spectrum. The linearly polarized radiation of a single voxel is then given by

$$ \mathcal{P}_v^{GB} = \frac{A_{\parallel} B_v(T_{\parallel}) - A_{\perp} B_v(T_{\perp})}{2} + i \frac{A_{\times} B_v(T_{\times}) - A_{\otimes} B_v(T_{\otimes})}{2} \quad (12) $$

Starting with this voxel SED renders the problem quite complicated. For example, just considering the SED of $Q_v$, we can

write

$$ \frac{A_{\parallel} B_v(T_{\parallel}) - A_{\perp} B_v(T_{\perp})}{2} = \frac{(A_{\parallel} - A_{\perp}) B_v(T_{\parallel}) + B_v(T_{\perp})}{2} + \frac{A_{\parallel} + A_{\perp}}{2} \frac{B_v(T_{\parallel}) - B_v(T_{\perp})}{2} \quad (13) $$

This means that two fundamental SED shapes are required: the sum and difference of two blackbody spectra. These are generally not blackbody spectra again (Chluba & Sunyaev 2004) and the moment expansion requires two series. If the number of spectral parameters is extended (here it was only $T$), then the number of fundamental voxel spectra increases rapidly, which can quickly make the situation quite complicated.

In astrophysical applications, it is commonly assumed that the only source of polarization is through variations of the number of emitters (i.e., the weight parameter $A$) at fixed SED parameters. In our example, this means $T_{\parallel} \approx T_{\perp} \approx T_{\times} \approx T_{\otimes} = T$, such that $Q_v \approx \frac{(A_{\parallel}-A_{\perp})}{2} B_v(T)$. In this case, the fundamental voxel SED is given by $\hat{P}_v^{GB} \approx B_v(T)$, such that the single polarization state can be characterized by $\mathcal{P}_v^{GB} \approx A e^{2i\gamma} B_v(T)$. For the foreground examples treated below, we similarly assume that inside a given voxel the spectral parameters remain constants. We further discuss how to go beyond this assumption in Sect. 6.1.

### 3.4. Spin moments: Moment expansion for spin-2 quantities

We now generalize the moment expansion for intensity presented in Sect. 2 to polarized signals. We discuss how this new framework arises naturally from the previous one and provides a powerful tool allowing us to model the frequency dependence of $\gamma_v$ in the presence of polarized mixing.

The generalization is indeed quite straightforward. For the intensity moment expansion, we performed a Taylor series in the spectral parameters for each emitting volume element, in Eq. (3). The line-of-sight average in one direction, $\hat{n}$, is then given by Eq. (2). For polarization, this is equivalent to performing a Taylor expansion of the spinor's modulus<fnref n="4" /> $P_v = |\mathcal{P}_v|$ with respect to the spectral parameters at each fixed polarization angle $\gamma$. However, we cannot use a perturbative approach to average over the polarization angles $\gamma$ since, in physical situations, one expects them to vary widely in a nontrivial way, such that the situation would quickly become mathematically inconsistent. This means that the line-of-sight average has to be generalized to include the polarization state in the parameter distribution function

$$ \langle \mathcal{P}_v \rangle = \langle P_v(A, \mathbf{p}) e^{2i\gamma} \rangle = \int \mathcal{P}(\mathbf{p}, \gamma, \hat{n}) \hat{P}_v(\mathbf{p}) e^{2i\gamma} d^N p \, dy. \quad (14) $$

Here, we again used $P_v(A, \mathbf{p}) = A \hat{P}_v(\mathbf{p})$, as for intensity. In analogy to the intensity moment expansion, one then finds

$$ \langle \mathcal{P}_v(A, \mathbf{p}, \gamma) \rangle = \hat{P}_v(\bar{\mathbf{p}}) \langle A e^{2i\gamma} \rangle + \sum_j^N \langle A e^{2i\gamma} (p_j - \bar{p}_j) \rangle \partial_{\bar{p}_j} \hat{P}_v(\bar{\mathbf{p}}) \\ + \frac{1}{2} \sum_{j,k}^N \langle A e^{2i\gamma} (p_j - \bar{p}_j)(p_k - \bar{p}_k) \rangle \partial_{\bar{p}_j} \partial_{\bar{p}_k} \hat{P}_v(\bar{\mathbf{p}}) + \dots, $$

which depends on the pivot $\bar{\mathbf{p}}$, as we specify in Sect. 3.4.2. Since $\langle A e^{2i\gamma} \rangle$ can vanish, we cannot simply factor it out of the expressions. Instead, like for the intensity moments, we will again use

<footnote n="4">⁴ Here we make the assumptions discussed in Sect. 3.3.</footnote>
$\bar{A} = \langle |A e^{2i\gamma}| \rangle \equiv \langle A \rangle$ as a way to normalize the distributions. This allows us to define the spin moments

$$W_{\alpha}^{p_j \cdots p_l} = \frac{\langle A e^{2i\gamma} (p_j - \bar{p}_j) \cdots (p_l - \bar{p}_l) \rangle}{\bar{A}} \equiv \Omega_{\alpha}^{p_j \cdots p_l} e^{2i\gamma_{\alpha}^{p_j \cdots p_l}} \quad (15)$$

very much like for intensity but with an extra spinor weight. The moments are now complex-valued and in the second step we expressed them in terms of the real numbers, $\Omega_{\alpha}^{p_j \cdots p_l}$ and $\gamma_{\alpha}^{p_j \cdots p_l}$. The latter defines average directions of polarization states associated with each of the moments of the SED. While considering the pair $(\Omega_{\alpha}, \gamma_{\alpha})$ or the real and imaginary parts of $W_{\alpha}$ give perfectly equivalent descriptions, one could be favored over the other for parameters estimation or physical interpretation. From a numerical perspective, considering the pair $(Q, U)$ as the components of a single object instead of two independent intensities will add correlations between their moments, which is expected to improve the accuracy of the parameter inference. The final polarization moment expansion then takes the form:

$$\langle \mathcal{P}_v(A, \mathbf{p}, \gamma) \rangle = \mathcal{W}_0 P_v(\bar{A}, \bar{\mathbf{p}}) + \sum_j \mathcal{W}_1^{p_j} \partial_{\bar{p}_j} P_v(\bar{A}, \bar{\mathbf{p}}) \\ + \frac{1}{2} \sum_{j,k} \mathcal{W}_2^{p_j p_k} \partial_{\bar{p}_j} \partial_{\bar{p}_k} P_v(\bar{A}, \bar{\mathbf{p}}) + \dots \quad (16) \\ + \frac{1}{\alpha!} \sum_{j, \dots, l} \mathcal{W}_{\alpha}^{p_j \cdots p_l} \partial_{\bar{p}_j} \cdots \partial_{\bar{p}_l} P_v(\bar{A}, \bar{\mathbf{p}}) + \dots,$$

which in this form can be interpreted as the sum of multiple SEDs with well-defined polarization states. It is this sum of well-defined single polarization states (i.e., defined by the complex-valued moments) with varying SEDs (i.e., the derivative spectra) that leads to rotation of polarization planes.

We comment that the number of parameters in Eq. (16) depends on the moment order that is used in the modeling. For each moment, two degrees of freedom are added (i.e., the real and imaginary parts). In addition, one has to determine the spectral parameter pivot, $\bar{\mathbf{p}}$. However, the overall normalization $\bar{A}$ does not independently contribute, but was merely chosen to scale the moments. As such, it cannot be independently estimated, and only the values of $\bar{A} W_{\alpha}$ actually matter.

### 3.4.1. Average polarization angle

Since $W_0 = \Omega_0 e^{2i\gamma_0} = \langle A e^{2i\gamma} \rangle / \bar{A}$ can generally vanish, there is no longer a trivially defined average polarization angle. In particular when $W_0 \approx 0$, the average polarization angle can be fully determined by the higher order terms in Eq. (16) and also generally becomes frequency-dependent.

To illustrate this aspect, let us consider the simple example of two power law spectra with equal weights $A$ along the $\pm Q$ direction ($\gamma_1 = 0$ and $\gamma_2 = \pi/2$). For these we have $e^{2i\gamma_1} = 1$ and $e^{2i\gamma_2} = -1$ implying $\langle \mathcal{P}_v \rangle = A(\nu/v_0)^{\beta_1} - A(\nu/v_0)^{\beta_2}$. For $\beta_1 \neq \beta_2$, we find $\mathcal{P}_v \neq 0$ unless $\nu \equiv \nu_0$, which is reflected by the fact that $W_0 = \langle A e^{2i\gamma} \rangle / \bar{A} \equiv (A - A)/[2A] = 0$, implying that the leading order term in Eq. (16) vanishes. Also, no matter what the frequency, in our example the polarization state will remain $Q_v$, with a change of sign at $\nu = \nu_0$ and hence flip of $0 \leftrightarrow \pi/2$. In this situation, all higher order moments remain real and $\gamma_v$ is highly non perturbative (i.e., not differentiable) at $\nu = \nu_0$.

There must be a way to define a meaningful average polarization angle for each of the moment terms. Indeed, if we simply

think of the average of $\gamma$ along the line of sight in terms of the distribution, $\mathbb{P}(\mathbf{p}, \gamma, \hat{\mathbf{n}})$. This then results in

$$\bar{\gamma} = \frac{\langle A \gamma \rangle}{\langle A \rangle} = \frac{\int \mathbb{P}(\mathbf{p}, \gamma, \hat{\mathbf{n}}) \gamma \, d^N \mathbf{p} \, d\gamma}{\int \mathbb{P}(\mathbf{p}, \gamma, \hat{\mathbf{n}}) \, d^N \mathbf{p} \, d\gamma} \quad (17)$$

as the average polarization angle. This angle can also be used as a pivot when expanding the polarization state:

$$e^{2i\gamma} = e^{2i\bar{\gamma}} \left[ 1 + \sum_{k=1}^{\infty} \frac{(2i)^k}{k!} (\gamma - \bar{\gamma})^k \right]. \quad (18)$$

Using this in Eq. (15), have

$$W_{\alpha}^{p_j \cdots p_l} = e^{2i\bar{\gamma}} \sum_{k=0}^{\infty} \frac{(2i)^k \langle A(\gamma - \bar{\gamma})^k (p_j - \bar{p}_j) \cdots (p_l - \bar{p}_l) \rangle}{k! \bar{A}} \quad (19)$$

The first term in the sum (i.e., $k = 0$), is the only non vanishing contribution if the distributions of $\gamma$ and $\mathbf{p}$ factorize (i.e., the two are uncorrelated variables), as we discuss in Sect. 3.4.3. Adding term by term in the series of Eq. (19) allows us to include information from higher order correlations of $\gamma$ and $\mathbf{p}$. However, in terms of distinguishable parameters, only the total moments, $W_{\alpha}^{p_j \cdots p_l}$, can really be constrained.

### 3.4.2. Definition of the pivot

How do we determine the spectral parameter pivot? In the intensity case, we simply demanded the first moments to vanish to fix the pivot. For polarization, this naively yields the condition

$$\langle A e^{2i\gamma} \rangle \bar{\mathbf{p}} = \langle A e^{2i\gamma} \rangle \mathbf{p}. \quad (20)$$

However, since $\langle A e^{2i\gamma} \rangle$ can vanish, in general this cannot be a meaningful choice.

Above, we already defined $\bar{A} = \langle |A e^{2i\gamma}| \rangle$. In a similar manner, we can introduce the SED pivots as

$$\bar{\mathbf{p}} = \frac{\langle |A e^{2i\gamma}| \mathbf{p} \rangle}{\langle |A e^{2i\gamma}| \rangle} \equiv \frac{\langle A \mathbf{p} \rangle}{\langle A \rangle} = \frac{\int \mathbb{P}(\mathbf{p}, \gamma, \hat{\mathbf{n}}) \, \mathbf{p} \, d^N \mathbf{p} \, d\gamma}{\int \mathbb{P}(\mathbf{p}, \gamma, \hat{\mathbf{n}}) \, d^N \mathbf{p} \, d\gamma}, \quad (21)$$

which is equivalent to the definition for the intensity moments. Physically, this means that we disregard the geometrical properties of $\mathcal{P}_v$ and simply treat its modulus as an intensity. For our power-law example in Sect. 3.4.1, this means $\bar{\beta} = (A_1\beta_1 + A_2\beta_2)/(A_1 + A_2)$, which is fully analogous to the result of a simple intensity moment expansion. As we shall see below, this choice is well motivated and leads to a well-behaved polarization moment formalism.

In the perturbative regime however, $\Omega_0 = |\langle A e^{2i\gamma} \rangle| \gg 0$, one can safely choose the complex pivot

$$\bar{\mathbf{p}} = \frac{\langle A e^{2i\gamma} \mathbf{p} \rangle}{\langle A e^{2i\gamma} \rangle} \Rightarrow \Delta \bar{p}_j = \frac{\langle A e^{2i\gamma} (p_j - \bar{p}_j) \rangle}{\langle A e^{2i\gamma} \rangle} = \frac{W_1^{p_j}}{W_0}. \quad (22)$$

While the spectral parameters $\mathbf{p}$ are real quantities, correcting by a complex number might seem incoherent. However, as we will discuss with examples, doing so is deeply relevant. While the real part of $\bar{\mathbf{p}}$ can be interpreted as real correction of $\mathbf{p}$, its complex part gives rise to the first order frequency dependence of the polarization angle $\gamma_v$ and can add some spectral modulation to the polarized intensity.
### 3.4.3. Independent angle distribution and de-polarization

In the definition of the line-of-sight average and spin moments, Eqs. (14) and (15), we kept the parameter distribution function general. The discussion is greatly simplified if the probability distributions for the spectral parameters and the polarization angles can be considered as independent. In this case, one has<fnref n="5" /> $P(p, \gamma, \hat{n}) \approx P(p, \hat{n})P(\gamma, \hat{n})$, such that

$$ \langle \mathcal{P}_v \rangle = \int \mathcal{P}(p, \hat{n}) \hat{P}_v(p) \, d^N p \int \mathcal{P}(\gamma, \hat{n}) e^{2i\gamma} \, d\gamma. \quad (23) $$

As this expression shows, spectral mixing and polarization angle averaging become completely independent, such that no frequency-dependent polarization angle can be expected. However, when summing over different physical conditions along the line of sight, the probability distribution becomes

$$ \begin{aligned} \mathcal{P}(p, \gamma, \hat{n}) = & \int A(s)\delta(\gamma - \gamma(s))\delta^N(p - p(s))ds \\ & \neq \int \delta(\gamma - \gamma(s))ds \int A(s)\delta^N(p - p(s))ds, \end{aligned} \quad (24) $$

introducing an unavoidable dependence between the angles and the spectral parameters. This dependence disappears if either the polarization angle or the spectral parameters are constant in the line of sight, highlighting that a variation of both $\gamma$ and $p$ is required to have a spectral dependence of the polarization angle.

If the angle distribution is Gaussian with average angle $\tilde{\gamma}(\hat{n})$ and width $\sigma_{\gamma}(\hat{n})$, then one finds

$$ \int \mathcal{P}(\gamma, \hat{n}) e^{2i\gamma} d\gamma = e^{2i\tilde{\gamma}(\hat{n})} e^{-2\sigma_{\gamma}^2(\hat{n})}. \quad (25) $$

This expression highlights that the dispersion of the angles leads to damping of the net polarization amplitude and ultimately complete depolarization if the distribution becomes too wide. In this case, a general perturbative expansion in $\Delta\gamma = \gamma - \tilde{\gamma}$, (see e.g. Eq. (18)) is unlikely to converge, but, as stressed already, does not add any new insight anyways.

## 4. Canonical SEDs

In this section, we illustrate the spin-moment framework on some detailed analytical and numerical examples relevant to astrophysical applications. We consider discrete sums of polarized SEDs along a given line of sight, often focusing on very few contributions. For the moment formalism, this can lead to non-perturbative cases, since in the limit of many emitters, the moments are expected to become more Gaussian due to the central limit theorem. Still, in most cases only a few moments are required to capture the dominant effects.

To highlight the performance of the moment formalism, we treat the sum of SEDs with noise as data and then use the moment representations to finite order as model. We perform a parameter estimation by means of curve fitting with $\chi^2$ minimization in complex-variables using the LMFIT python library (Newville et al. 2016). Hereafter, the model of linear polarization given by the spin-moment expansion is $P_v^M$ and the simulated data signal is noted $P_v^S$. We add Gaussian noise $N_v$ to the simulation, with zero mean and standard deviation $\sigma = \sigma_Q + i\sigma_U$. The values of $\sigma$ is chosen such that the signal to noise ratios

$Q_v^S/Q_v^M$ and $U_v^S/U_v^M$ are constants over the whole frequency range (chosen arbitrarily to be $1 \times 10^{-5}$)<fnref n="6" />.

The $\chi^2$ to minimize is given by $\chi^2 = \frac{1}{2}|\mathcal{P}_v^M - \mathcal{P}_v^S|^2/|\sigma|^2$. The signal is considered over a frequency range going from 1 GHz to $\nu_{max}$ in intervals of 1 GHz. The choice of $\nu_{max}$ will depend on the example considered. We introduce the shorthand notation 'O($\alpha$)' to refer to the fit of the spin-moment expansion including all the terms up to order $\alpha$. 'O(0)' is the leading order/canonical SED. Two distinct routines are developed, fitting either the pair (Re($W_{\alpha}^p$), Im($W_{\alpha}^p$)) or the pair ($\Omega_{\alpha}^p$, $\gamma_{\alpha}^p$). In all the examples considered, both lead to identical results and we leave a further comparison between the two implementations for future work. We are interested in the SED distortions and their behavior in the complex plane, which are only driven by the relative contributions of the different emission points. As such, we use natural units of Jy sr$^{-1}$, for all the SEDs. A more detailed discussion on the relevance of weights, normalization and change of units can be found in Appendices A and B.

### 4.1. General discrete sums of canonical SEDs

For a discrete sum of M SEDs along a line of sight one can trivially write the distribution function as

$$ \mathcal{P}(p, \gamma, \hat{n}) = \sum_{k}^{M} A_k \delta(\gamma - \gamma_k) \delta^{N}(p - p_k), \quad (26) $$

where $\delta(x - x_0)$ denotes Dirac's distribution and the sum extends over the discrete emission points along the line of sight with SED vectors $p_k$ and polarization angles $\gamma_k$. Inserting this into the definitions of the moments and pivots given in the previous section we trivially find the exact average

$$ \begin{aligned} \langle \mathcal{P}_v \rangle &= \sum_{k}^{M} A_k \int \delta(\gamma - \gamma_k) \delta^N(p - p_k) \hat{P}_v(p) e^{2i\gamma} \, d^N p \, d\gamma \\ &= \sum_{k}^{M} A_k e^{2i\gamma_k} \hat{P}_v(p_k). \end{aligned} \quad (27) $$

Using the polarization moment expansion, we automatically have the normalization, pivot and complex-valued moments as

$$ \bar{A} = \sum_{k}^{M} A_k \int \delta(\gamma - \gamma_j) \delta^N(p - p_j) \, d^N p \, d\gamma = \sum_{k}^{M} A_k \quad (28a) $$

$$ \bar{\mathbf{p}} = \sum_{k}^{M} (A_k / \bar{A}) \mathbf{p}_k, \quad \bar{\gamma} = \sum_{k}^{M} (A_k / \bar{A}) \gamma_k \quad (28b) $$

$$ W_0 = \sum_{k}^{M} (A_k / \bar{A}) e^{2i\gamma_k} \quad (28c) $$

$$ W_{\alpha}^{p_j \dots p_l} = \sum_{k}^{M} (A_k / \bar{A}) e^{2i\gamma_k} (p_{k,j} - \bar{p}_{k,j}) \dots (p_{k,l} - \bar{p}_{k,l}), \quad (28d) $$

where the ratios $A_k/\bar{A}$ determine the probabilities to find $p_k$ and $\gamma_k$. These expressions can then be inserted into Eq. (16) to obtain the polarization moment expansion. The derivatives of the spectra have to be computed individually, but generally the moment expansion is expected to converge with only a few terms.

<footnote n="5">⁵ In doing so, we can use the normalizations $\int \mathcal{P}(p, \hat{n})d^N p = \bar{A}$ and $\int \mathcal{P}(\gamma, \hat{n})d\gamma = 1$.</footnote>

<footnote n="6">⁶ The error bars used in all the figures are respectively given by $\sigma_Q^2$, $\sigma_U^2$, $\sigma_P^2 = (Q^2\sigma_Q^2 + U^2\sigma_U^2)/P^2$ and $\sigma_\gamma^2 = 0.5(U^2\sigma_Q^2 + Q^2\sigma_U^2)/P^4$.</footnote>
## 4.2. Power laws

As a first example of astrophysical relevance, we consider the simple case of power-law SEDs:

$$ \hat{P}_{\nu}^{\text{PL}}(\beta) = \left(\frac{\nu}{\nu_0}\right)^{\beta}. \qquad (29) $$

The polarization state can then be characterized by $P_{\nu}^{\text{PL}} \approx Ae^{2i\gamma} \hat{P}_{\nu}^{\text{PL}}$. The only spectral parameter relevant for the moment expansion is the spectral index $p = (\beta)$, normalized at a reference frequency $\nu_0$. This SED plays a crucial role in the foreground modeling of synchrotron on large scales (Planck Collaboration IV 2020). In the following numerical applications, we choose $\nu_{\text{max}} = 150$ GHz below which the synchrotron emission is dominant and $\nu_0 = 23$ GHz as the WMAP frequency band (Bennett et al 2013). Using $\partial_{\beta}^k \hat{P}_{\nu}^{\text{PL}} = \partial_{\beta}^k (\nu/\nu_0)^\beta = (\nu/\nu_0)^\beta \ln(\nu/\nu_0)^k$, the spin-moment expansion in Eq. (16) can then be expressed as

$$ \langle \mathcal{P}_{\nu}^{\text{PL}} \rangle = P_{\nu}^{\text{PL}}(\bar{A}, \bar{\beta}) \times \left\{ \mathcal{W}_0 + \mathcal{W}_1^{\beta} \ln\left(\frac{\nu}{\nu_0}\right) + \frac{\mathcal{W}_2^{\beta^2}}{2} \ln\left(\frac{\nu}{\nu_0}\right)^2 + \frac{\mathcal{W}_3^{\beta^3}}{6} \ln\left(\frac{\nu}{\nu_0}\right)^3 + \dots \right\}. \qquad (30) $$

The choice of the reference frequency $\nu_0$, around which to make the expansion, can have an impact on the convergence rate of model, but otherwise leaves the moment expansion unchanged. One choice is to pick a local extremum where the SED changes shape: $\partial_{\nu} P_{\nu}^{S} = 0$ or $\partial_{\nu} \gamma_{\nu}^{S} = 0$ depending on the distortion type. In front of real data the choice has to be made also from instrumental considerations. As an example, we now consider the superposition of two power laws in more detail.

### 4.2.1. Hands-on example: Two power laws

Consider two power laws ($M=2$) with different spectral indices ($\beta_1, \beta_2$) and polarization angle ($\gamma_1, \gamma_2$) along the same line of sight. The exact solution then reads

$$ \mathcal{P}_{\nu} = A_1 \left(\frac{\nu}{\nu_0}\right)^{\beta_1} e^{2i\gamma_1} + A_2 \left(\frac{\nu}{\nu_0}\right)^{\beta_2} e^{2i\gamma_2}. \qquad (31) $$

Carrying out the intensity moment expansion of the two individual power laws with respect to their spectral indices, we obtain

$$ \bar{A} = A_1 + A_2 \qquad (32a) $$

$$ \bar{\beta} = \frac{A_1}{\bar{A}} \beta_1 + \frac{A_2}{\bar{A}} \beta_2, \quad \bar{\gamma} = \frac{A_1}{\bar{A}} \gamma_1 + \frac{A_2}{\bar{A}} \gamma_2 \qquad (32b) $$

$$ \mathcal{W}_0 = \frac{A_1}{\bar{A}} e^{2i\gamma_1} + \frac{A_2}{\bar{A}} e^{2i\gamma_2} \qquad (32c) $$

$$ \mathcal{W}_{\alpha}^{\beta\alpha} = \frac{A_1}{\bar{A}} e^{2i\gamma_1}(\beta_1 - \bar{\beta})^{\alpha} + \frac{A_2}{\bar{A}} e^{2i\gamma_2}(\beta_2 - \bar{\beta})^{\alpha}. \qquad (32d) $$

These expressions can be trivially extended to M power laws after extending the sums to M parameters $A_k, \gamma_k$ and $\beta_k$ to find the values of the pivot and spin moments (see Sect. 4.1). However, for illustrations the two power-law case is more intuitive.

If $\gamma_1 = \gamma_2 = \bar{\gamma}$ and $\beta_1 \neq \beta_2$, we naturally find that all spin moments are aligned in the same directions of the complex-plane and hence no change in the polarization direction can occur as a function of frequency. In this case, $\langle \mathcal{P}_{\nu}^{\text{PL}} \rangle = e^{2i\gamma} \langle P_{\nu}^{\text{PL}}(A, \beta) \rangle$, trivially describing the effect of spectral mixing only. If on the other

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 2. Illustration of the spinor $P_{\nu}$ in the complex plane $(Q, U)$ for a sum of two power laws (upper panel). Black crosses mark the steps of 50 GHz on the signal. The values of frequencies are indicated above the crosses in GHz. The corresponding polarized intensity $P_{\nu}$ (central panel) and polarization angle $\gamma_{\nu}$ (lower panel). The exact result with associated (invisible) error bars (black) is compared to the best-fit moment representation at various orders.</figcaption></figure>

hand $\beta_1 = \beta_2 = \bar{\beta}$ and $\gamma_1 \neq \gamma_2$, we naturally have $\mathcal{W}_{\alpha}^{\beta\alpha} = 0$ and everything is described by $\mathcal{W}_0$ with a fixed SED. To obtain non-trivial consequences of polarized mixing, both $\gamma_k$ and $\beta_k$ need to vary.

As an illustration, let us consider the highly non-perturbative example, with synchrotron-like behavior: $A_1 = 2$ Jy sr$^{-1}$, $A_2 = 1$ Jy sr$^{-1}$, $\beta_1 = -2.8$, $\beta_2 = -3.6$, $2\gamma_1 = 10^{\circ}$ and $2\gamma_2 = 80^{\circ}$, implying $\bar{A} = 3$ Jy sr$^{-1}$, $\bar{\beta} \approx -3.06$ and $\bar{\gamma} \approx 16.6^{\circ}$. In Fig. 2, the
modulus $P_v$ and the argument $\gamma_v$ of the signal together with the recovered spinor representation are displayed for various orders of the expansion going from $O(0)$ to $O(3)$. For such strong deviations of spectral parameters, one cannot expect to find $\bar{\beta} \approx -3.06$ for the expansion at leading order, and it has to be treated as a free parameter of the model. This is particularly important when only a few moment terms are included. The expansion at higher orders then allows us to gradually recover the nontrivial polarization signal over the frequency range. In all cases but $O(0)$, the best fit values for $\bar{A}, \bar{\beta}$ and $W_\alpha^\beta$ are all compatible within one standard deviation with those given by Eq. (32). While, by definition, the leading order cannot encompass any rotation of the spinor with frequency, we can see that the moment expansion allows us to correctly model the frequency dependence of $\gamma_v$.

### 4.2.2. Extreme cases and perturbative regime

To gain further insight, let us just consider the first order terms of the expansion:

$$ \langle \mathcal{P}_v^{\text{PL}} \rangle \approx \hat{P}_v^{\text{PL}}(\bar{A}, \bar{\beta}) \times \left\{ \mathcal{W}_0 + \mathcal{W}_1^{\beta} \ln\left(\frac{\nu}{\nu_0}\right) \right\}. \quad (33) $$

If $\mathcal{W}_0 \approx 0$, we indeed find the situation where we have a polarization angle fully determined by $W_1^\beta$, with a sign-flip at $\nu = \nu_0$. The polarization SED is then determined by the first $\beta$ derivative of the power law, and polarization rotation would stem from higher order moments which are not included here. In this situation, we are dealing with two dominant (and near degenerate) contributions to the polarization state that are rotated by 90° to each other (e.g., $+Q_v$ and $-Q_v$). The moment expansion then describes how much these two power-law terms differ.

This situation is illustrated in Fig. 3 by a sum of two power laws with parameters $A_1 = A_2 = 2$ Jy sr$^{-1}$, $\beta_1 = 1$, $\beta_2 = 2$, $2\gamma_1 = 180^\circ$ and $2\gamma_2 = 0.1^\circ$. One can see that both $\gamma_v$ and $P_v$ do not behave as smooth functions at the breaking point $\nu \approx \nu_0$, where one power law abruptly takes over the other one. Here, $\gamma_v$ changes very rapidly from $\gamma_1$ to $\gamma_2$. This is, however, not a problem for the spin-moment expansion, which allows us to recover $\mathcal{P}_v$ correctly. As mentioned above, one recovers a very small best-fit value for $\mathcal{W}_0$ and $W_1^\beta \neq 0$ is the dominant term of the expansion, and polarization rotation mainly stems from the second and higher order moments.

If on the other hand we are in the perturbative regime where $|\mathcal{W}_0| \gg |W_1^\beta| > 0$, we can correct for the first order term of the expansion as in intensity, interpreting it as a correction to the spectral parameters with this correction now being complex-valued. let us split its real and imaginary parts as

$$ \Delta\bar{\beta} = \frac{W_1^\beta}{W_0} = a_{\Delta\bar{\beta}} + ib_{\Delta\bar{\beta}}. \quad (34) $$

With $^7\mathcal{W}_0 = \Omega_0 e^{2i\gamma_0}$, we can then write the expansion as

$$ \begin{aligned} \langle \mathcal{P}_v^{\text{PL}} \rangle &\approx \mathcal{W}_0 \hat{P}_v^{\text{PL}}(\bar{A}, \bar{\beta}) \times \left\{ 1 + \Delta\bar{\beta} \ln\left(\frac{\nu}{\nu_0}\right) \right\} \\ &\approx \bar{A} \mathcal{W}_0 \left(\frac{\nu}{\nu_0}\right)^{\bar{\beta}+a_{\Delta\bar{\beta}}+ib_{\Delta\bar{\beta}}} \\ &\equiv \bar{A} \Omega_0 \exp\left(2i\left[\gamma_0 + \frac{b_{\Delta\bar{\beta}}}{2} \ln\left(\frac{\nu}{\nu_0}\right)\right]\right) \left(\frac{\nu}{\nu_0}\right)^{\bar{\beta}+a_{\Delta\bar{\beta}}} . \end{aligned} \quad (35) $$

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 3. Illustration of the spinor $\mathcal{P}_v$ in the complex plane $(Q, U)$ for a sum of two almost anti aligned power laws, with a close-up view around $\nu = \nu_0$ (upper panel). Black crosses mark the steps of 10 GHz on the signal.The values of frequencies are indicated above the crosses in GHz. The corresponding polarized intensity $P_v$ (central panel) and polarization angle $\gamma_v$ (lower panel). We choose logarithmic representations to emphasize the focus on the $\nu \sim \nu_0$ point.</figcaption></figure>

In the perturbative regime, we thus find a polarization SED that is again a power law with a spectral index and frequency-dependent polarization angle given by

$$ \bar{\beta}^{\text{PL}} \approx \bar{\beta} + \text{Re} \left( \frac{W_1^{\beta}}{W_0} \right) \quad (36a) $$

$$ \gamma_v^{\text{PL}} \approx \gamma_0 + \frac{1}{2} \text{Im} \left( \frac{W_1^\beta}{W_0} \right) \ln \left( \frac{\nu}{\nu_0} \right). \quad (36b) $$

<footnote n="7">⁷ In general $\Omega_0$ is expected to depart from unity.</footnote>
This result clearly shows that even at the lowest order, the superposition of linearly polarized power-law SEDs generally leads to a frequency-dependent rotation of the polarization angle. The rotation is purely driven by the imaginary part of the pivot value in the power-law index. At lowest order, no spectral curvature is added, even if the spectral index departs from that of the simple intensity superposition, $\bar{\beta}$, by $a_{\Delta\bar{\beta}}$. At higher order, the complex-valued moments given in Eq. (32) lead to additional spectral complexity, which in most relevant situations can be captured in a perturbative manner.

Just as for intensity, in the perturbative regime one can cancel $W_1^\beta$ and correct the leading order according to Eq. (35). Proceeding iteratively on numerical examples, one can witness the quick convergence of the first order moment toward zero. Doing so allows us to find the right pivot for the expansion while still keeping the pivot $\bar{\beta}$ fixed. As such, one breaks unwanted degeneracies and recovers a physically relevant value for the spectral parameters with a minimal dispersion.

### 4.3. Blackbodies

As the second example, we briefly consider the superposition of blackbody spectra, a case that is directly relevant to primordial CMB polarization. The only free parameter is the blackbody temperature, and a temperature difference in two orthogonal directions is required to obtain a net polarization. This is an example where the origin of the polarization is due to the spectral parameters only. The SED for polarized light, at leading order in the temperature perturbation around the average, is then given by the first temperature derivative of a blackbody:

$$ B_v(\bar{T}) = \frac{2h}{c^2} \frac{v^3}{e^{hv/k\bar{T}} - 1}, \quad G_v(\bar{T}) = \frac{\partial B_v(\bar{T})}{\partial \ln \bar{T}} = \frac{2hv^3}{c^2} \frac{x e^x}{(e^x - 1)^2}, \\ \mathcal{P}_v^{BB} \approx G_v(\bar{T}) [\Theta_Q + i\Theta_U] \quad (37) $$

with $x = hv/k\bar{T}$ and the usual natural constants. We also introduced the two temperature perturbations $\Theta_Q = \Delta T_Q/\bar{T}$ and $\Theta_U = \Delta T_U/\bar{T}$, which are respectively defined for a coordinate system that is rotated by 45°. At fixed $\Theta_Q$ and $\Theta_U$, this means that (at lowest order in the temperature fluctuations) no spectral mixing happens, and hence $\gamma$ remains frequency-independent.

However, if we include terms at second order in $\Theta$, one finds an additional frequency dependence that is characterized by a y-type distortion (e.g., see Appendix A of Chluba et al. 2015):

$$ Y_v(\bar{T}) = \frac{2hv^3}{c^2} \frac{x e^x}{(e^x - 1)^2} \left[ x \frac{e^x + 1}{e^x - 1} - 4 \right], \\ \mathcal{P}_v^{BB} = G_v(\bar{T}) [\Theta_Q + 2\Theta_I\Theta_Q + i(\Theta_U + 2\Theta_I\Theta_U)] \\ + Y_v(\bar{T}) [\Theta_I\Theta_Q + i\Theta_I\Theta_U]. \quad (38) $$

Here, we introduced the total intensity temperature perturbation $\Theta_I = \Delta T_I/\bar{T}$, which generally includes both polarized and unpolarized contributions. Depending on the ratio of $\Theta_Q$ to $\Theta_U$, this will cause a small frequency-dependent rotation of the polarization planes. However, since this effect is at second order in the (small) CMB temperature differences, we leave a more detailed discussion to future work.

### 4.4. Gray-body spectra

In contrast to the blackbody, a gray-body (GB) spectrum also has a free normalization, caused by imperfect reflectivity of the

material. As discussed in Sect. 3.3, we only allow variations of A to create polarization. The fundamental SED is thus given by $\hat{P}_v^{GB} \equiv B_v(T)$, such that the single polarization state can be characterized by $\mathcal{P}_v^{GB} \approx A e^{2i\gamma} B_v(T)$.

We can then consider general GB superpositions. Following Chluba et al. (2017), we shall use $\beta_{GB} = 1/T$ as the spectral parameter. The SED derivatives then have a closed form using Eulerian numbers (Chluba et al. 2013), with the first few terms given by (see Eq. (38) of Chluba et al. 2017):

$$ \beta_{GB} \partial_{\beta_{GB}} \hat{P}_v^{GB} = - \hat{P}_v^{GB} \frac{xe^x}{(e^x - 1)} \quad (39a) $$

$$ \beta_{GB}^2 \partial_{\beta_{GB}}^2 \hat{P}_v^{GB} = + \hat{P}_v^{GB} \frac{xe^x}{(e^x - 1)} x \coth(x/2) \quad (39b) $$

$$ \beta_{GB}^3 \partial_{\beta_{GB}}^3 \hat{P}_v^{GB} = - \hat{P}_v^{GB} \frac{xe^x}{(e^x - 1)} x^2 \frac{\cosh(x) + 2}{\cosh(x) - 1} \quad (39c) $$

$$ \beta_{GB}^4 \partial_{\beta_{GB}}^4 \hat{P}_v^{GB} = + \hat{P}_v^{GB} \frac{xe^x}{(e^x - 1)} \frac{x^3}{2} \frac{\cosh(x) + 5}{\sinh^2(x/2)} \coth(x/2). \quad (39d) $$

with the frequency variable $x = \frac{hv}{k}\bar{\beta}_{GB} \equiv \frac{hv}{k\bar{T}}$.

The final GB moment expansion then takes the form:

$$ \langle \mathcal{P}_v^{GB} \rangle = P_v^{GB}(\bar{A}, \bar{T}) \times \left\{ W_0 + W_1^{\beta_{GB}} Y_1^{GB}(x) + \frac{1}{2} W_2^{\beta_{GB}^2} Y_2^{GB}(x) + \frac{1}{6} W_3^{\beta_{GB}^3} Y_3^{GB}(x) + \dots \right\} \quad (40a) $$

$$ \bar{A} = \langle A \rangle, \quad W_0 = \frac{\langle A e^{2i\gamma} \rangle}{\bar{A}}, \quad \bar{\beta}_{GB} \equiv \frac{1}{\bar{T}} = \frac{\langle A \beta_{GB} \rangle}{\bar{A}} \quad (40b) $$

$$ W_{\alpha}^{\beta_{GB}} = \frac{\langle A e^{2i\gamma} (\beta_{GB} - \bar{\beta}_{GB})^{\alpha} \rangle}{\bar{A}}, \quad Y_{k}^{GB}(x) = \frac{1}{\hat{P}_{v}^{GB}} \frac{\partial^{k} \hat{P}_{v}^{GB}}{\partial \bar{\beta}_{GB}^{k}}. \quad (40c) $$

The functions $Y_k^{GB}(x)$ will also be relevant to the discussion of modified blackbody spectra in Sect. 4.5.

In Fig. 4, we fit the above model on a sum of two gray-bodies with parameters $A_1 = A_2 = 10^6$, $\beta_1 = 1/20 \text{ K}^{-1}$, $\beta_2 = 1/13 \text{ K}^{-1}$, $2\gamma_1 = 72^\circ$, $2\gamma_2 = -90^\circ$. To catch the domain on which $\beta_{GB}$ has a maximal impact, we choose $\nu_{max} = 5000 \text{ GHz}$. We note the "loop" trajectory of $\mathcal{P}_v$ in the complex plane $(Q, U)$ inherited from the combination of the shape of the black-body SED and the frequency rotation. By definition again, $O(0)$ can only be a straight line and fail to grasp this complexity. Even if this case again is non-perturbative, we see that the moment expansion up to third order allows us to gradually account for the SED distortions and recover $P_v$, the polarized mixing inducing a highly nontrivial rotation of the polarization angle.

Like in the power-law example, in the perturbative regime we can obtain a more general expression for the leading order terms. Assuming that $|W_0| \gg |W_1^{\beta_{GB}}|$, we can again use the split $\Delta\bar{\beta}_{GB} = W_1^{\beta_{GB}}/W_0 = a_{\Delta\bar{\beta}_{GB}} + ib_{\Delta\bar{\beta}_{GB}}$ into real and imaginary parts. With this, we can then write

$$ \langle \mathcal{P}_v^{GB} \rangle \approx \bar{A} W_0 \hat{P}_v^{GB}(\bar{T}) \times \left\{ 1 + \Delta \bar{\beta}_{GB} Y_1^{GB}(x) \right\} \\ \approx \bar{A} W_0 \hat{P}_v^{GB}(\bar{T}) \quad (41) $$
<figure><img src="image_9.png" /><figcaption>Fig. 4. Illustration of $P_v$ in the complex plane $(Q, U)$ for a sum of two gray-bodies (upper panel). Black crosses mark the steps of 1000 GHz on the signal. The values of frequencies are indicated above the crosses in GHz. The corresponding polarized intensity $P_v$ (central panel) and polarization angle $\gamma_v$ (lower panel).</figcaption></figure>

with $\tilde{T} = 1/(\bar{\beta}_{GB} + \Delta\bar{\beta}_{GB})$. The leading order SED term, $\hat{P}_v^{GB}(\tilde{T})$, then depends on the function

$$ \frac{1}{e^{\frac{hv}{kT}} - 1} = \frac{1}{e^{x_R + i x_I} - 1} = \frac{e^{2i\Delta\gamma_v^{GB}}}{\sqrt{(e^{x_R} - 1)^2 + 2e^{x_R} [1 - \cos(x_I)]}} \quad (42a) $$

$$ \Delta\gamma_v^{GB} = \frac{1}{2} \tan^{-1} \left( \frac{e^{x_R} \sin(x_I)}{e^{x_R} \cos(x_I) - 1} \right) \quad (42b) $$

with $x_R = hv(\bar{\beta}_{GB} + a_{\Delta\bar{\beta}_{GB}})/k$ and $x_I = hv b_{\Delta\bar{\beta}_{GB}}/k$. One can see that polarized mixing leads to an imaginary photon chemical potential, $\mu = i x_I$. This causes a frequency-dependent rotation of the polarization plane and also modifications to the SED. At high frequencies, one finds $\Delta\gamma_v \approx x_I/2$, while at low frequencies, one has the constant $\Delta\gamma_v \approx \frac{1}{2} \tan^{-1}(b_{\Delta\bar{\beta}_{GB}}/a_{\Delta\bar{\beta}_{GB}})$.

### 4.5. Modified blackbodies

Another highly relevant SED is given by the modified blackbody spectrum. It is expected to provide a good model for the thermal dust intensity and polarized signal (Planck Collaboration XVII 2014; Planck Collaboration XXII 2015). In principle, one should allow for amplitude, temperature and spectral index variations inside each voxel. This provides multiple ways of creating polarization (see Sect. 3.3 for discussion). However, we assume that the voxel polarization is again only given by amplitude variations in the four emission directions. The fundamental voxel SED then reads $\hat{P}_v^{mBB} \equiv (\nu/\nu_0)^{\beta_d} B_v(T)$, such that the single polarization state can be characterized by $P_v^{mBB} \approx A e^{2i\gamma}(\nu/\nu_0)^{\beta_d} B_v(T)^8$. Using the results for the power law and gray-body spectra of the previous sections, we then have

$$ \begin{aligned} \langle P_v^{mBB} \rangle = P_v^{mBB}(\bar{A}, \bar{T}, \bar{\beta}_d) \times & \{ W_0 + W_1^{\beta_d} \ln(\nu/\nu_0) + W_1^{\beta_{GB}} Y_1^{GB}(x) \\ & + \frac{1}{2} W_2^{\beta_d^2} \ln^2(\nu/\nu_0) + W_2^{\beta_{GB}\beta_d} \ln(\nu/\nu_0) Y_1^{GB}(x) + \frac{1}{2} W_2^{\beta_d^2} Y_2^{GB}(x) \\ & + \frac{1}{6} W_3^{\beta_d^3} \ln^3(\nu/\nu_0) + \frac{1}{2} W_3^{\beta_{GB}\beta_d^2} \ln^2(\nu/\nu_0) Y_1^{GB}(x) \\ & + \frac{1}{2} W_3^{\beta_{GB}\beta_d} \ln(\nu/\nu_0) Y_2^{GB}(x) + \frac{1}{6} W_3^{\beta_d^3} Y_3^{GB}(x) + \dots \} \end{aligned} \quad (43a) $$

$$ \bar{A} = \langle A \rangle, \quad W_0 = \frac{\langle A e^{2i\gamma} \rangle}{\bar{A}}, \quad (43b) $$

$$ \bar{\beta}_{GB} \equiv \frac{1}{\bar{T}} = \frac{\langle A \beta_{GB} \rangle}{\bar{A}}, \quad \bar{\beta}_d = \frac{\langle A \beta_d \rangle}{\bar{A}} \quad (43c) $$

$$ W_{\alpha+\delta}^{\beta_{GB}\beta_d^{\beta_d}} = \frac{\langle A e^{2i\gamma} (\beta_{GB} - \bar{\beta}_{GB})^{\alpha} (\beta_d - \bar{\beta}_d)^{\delta} \rangle}{\bar{A}} \quad (43d) $$

up to third order. Due to the dimensionality of the problem, the moment representation quickly becomes cumbersome, but can be easily handled using modern computers.

In Fig. 5, we applied this expansion on the sum of two normalized<fnref n="9" /> modified blackbodies of parameters $A_1 = A_2 = 1$, $\beta_1 = 2$, $\beta_2 = 1$, $T_1 = 5$ K, $T_2 = 70$ K, $2\gamma_1 = 72^\circ$, $2\gamma_2 = -90^\circ$. To simulate the thermal dust signal over the CMB missions frequency ranges and in order to witness the transition between the effect of the power-law factor at low frequencies and the gray-body factor at high frequencies, we choose $\nu_{max} = 800$ GHz. Accordingly to the Planck high frequency bands, we choose $\nu_0 = 353$ GHz (Planck Collaboration I 2016). One can see that the different power laws induce strong distortions at low frequencies $\le 100$ GHz while the temperature induces an additional bending at high frequencies. In this very extreme case, all the moments up to order 3 are required to correctly model the signal. However, even in this nontrivial situation the moment expansion

<footnote n="8">⁸ This is the standard way to model polarized dust emission locally, with $\bar{A} = \pi_{\text{dust}} \tau_{\text{dust}} \cos^2(\Gamma_{\text{dust}})$, $\pi_{\text{dust}}$ being the polarization fraction, $\tau_{\text{dust}}$ the opacity and $\Gamma_{\text{dust}}$ the angle between the Galactic magnetic field and the plane of the sky (Draine & Fraisse 2009). In practice one would have to consider a correlation between $\Gamma_{\text{dust}}$ and $\gamma$.</footnote>

<footnote n="9">⁹ The two modified blackbodies are here normalized by a reference blackbody at $\nu = \nu_0$, as further discussed in Appendix A. As such, $Q, U$ and $P_v$ are unitless.</footnote>
<figure><img src="image_10.png" /><figcaption>Fig. 5. Illustration of the spinor $P_v$ in the complex plane $(Q, U)$ for a sum of two modified-blackbodies (upper panel). Black crosses mark the steps of 100 GHz on the signal. The values of frequencies are indicated above the crosses in GHz. The corresponding polarized intensity $P_v$ (central panel) and polarization angle $\gamma_v$ (lower panel).</figcaption></figure>

performs extremely well. The leading order fit $O(0)$ interprets a local maximum of the polarized intensity and as the peak of the gray-body spectra, leading to a wrongly small value for the recovered temperature. In astrophysical situations one expects the central limit theorem to render the examples more moderate, with fewer moments required at a given precision.

In the perturbative regime, using the results from the previous sections, the leading order moment description then is

$$ \langle \mathcal{P}_v^{\text{mBB}} \rangle \approx \frac{2h\nu^3}{c^2} \frac{\bar{A} \Omega_0 (\nu/\nu_0)^{\bar{\beta}_d + a\Delta\bar{\beta}_d} e^{2i\gamma_v}}{\sqrt{(e^{x_R} - 1)^2 + 2e^{x_R} [1 - \cos(x_I)]}} \quad (44a) $$

$$ \gamma_v = \gamma_0 + \frac{1}{2} b_{\Delta\bar{\beta}_d} \ln(\nu/\nu_0) + \Delta\gamma_v^{\text{GB}}, \quad (44b) $$

with definitions as in the previous section. This expression demonstrates that the rotation of the polarization plane now has two contributions, one from the power-law modulation, one from the temperature terms. The rotation caused by temperature terms is particularly important at high frequencies and can become rapid due to a near linear scaling with $\nu$.

## 5. Generalizations of the formalism

In this section, we discuss additional averaging processes, covering the spherical harmonic decomposition and beam averaging effects. These cases all naturally lead to a redefinition of the meaning and values of the moments and spectral pivot, but they do not actually change the structure of the moment representation. For each case, we briefly recap how the problem is treated in intensity before generalizing to spin moments.

### 5.1. Generalization to spherical harmonics

As it was done for intensity in Chluba et al. (2017), one can immediately generalize the spin-moment expansion in harmonic space. For that we have to leave the above restriction of considering a single line of sight $\hat{n}$ and consider the intensities as fields over the celestial sphere. Using the moment expansion formalism at the power spectra level is especially useful for component separation on large sky fractions as it has been already shown for Planck data (Mangilli et al. 2021), The Simons Observatory telescope (Azzoni et al. 2021) and LiteBIRD (Remazeilles et al. 2021; Vacher et al. 2022).

#### 5.1.1. In intensity

While in the above, we described the signal along a given line of sight $\hat{n}$, we now turn ourselves to averages between different lines-of-sights across sky patches. The average intensity, $\langle I_v \rangle$, which generally depends on the line-of-sight moments and pivots of the SED expansion in the direction $\hat{n}$, is a scalar field on the $S^2$-manifold. As such, it can be expanded on the orthogonal basis of the spherical harmonic functions $Y_{\ell m}(\hat{n})$. We shall denote the spherical harmonic coefficients of a quantity $X$ as

$$ (X)_{\ell m} = \int Y_{\ell m}^*(\hat{n}) X(\hat{n}) \, d^2\hat{n}. \quad (45) $$

The spherical harmonic coefficients of the average intensity SED can then be expressed as:

$$ \begin{aligned} \langle I_v \rangle_{\ell m} &\equiv \langle (I_v)_{\ell m} \rangle = \int Y_{\ell m}^*(\hat{n}) \mathbb{P}(\hat{p}, \hat{n}) \hat{I}_v(\hat{p}) \, d^N \hat{p} \, d^2 \hat{n} \\ &= \int \mathbb{P}_{\ell m}(\hat{p}) \hat{I}_v(\hat{p}) \, d^N \hat{p}, \end{aligned} \quad (46) $$

where $\mathbb{P}_{\ell m}(\hat{p})$ is the harmonic coefficient of the parameter distribution, $\mathbb{P}(\hat{p}, \hat{n})$. Just as the single line-of-sight average, this superposition will introduce some additional mixing between intensities with different spectral parameters $\hat{p}(\hat{n})$ and introduces extra distortions. The spherical harmonic functions introduce a reweighting of the parameter distribution, which essentially
implies that all moments can be directly computed using the harmonic coefficients of the parameter distribution, ($P_{\ell m}(\mathbf{p})$.

Expanding $\hat{I}_v(\mathbf{p})$ in Eq. (46) and following the same steps introduced in Sect. 2, one can derive the moment expansion in harmonic space having the exact same structure as Eq. (6):

$$ \begin{aligned} \langle I_v \rangle_{\ell m} ={}& \delta_{\ell 0} I_v(\bar{A}, \bar{\mathbf{p}}) + \sum_j (\omega_1)_{\ell m}^{p_j} \partial_{\bar{p}_j} I_v(\bar{A}, \bar{\mathbf{p}}) \\ & + \frac{1}{2} \sum_{j,k} (\omega_2)_{\ell m}^{p_j p_k} \partial_{\bar{p}_j} \partial_{\bar{p}_k} I_v(\bar{A}, \bar{\mathbf{p}}) \\ & + \frac{1}{3!} \sum_{j,k,l} (\omega_3)_{\ell m}^{p_j p_k p_l} \partial_{\bar{p}_j} \partial_{\bar{p}_k} \partial_{\bar{p}_l} I_v(\bar{A}, \bar{\mathbf{p}}) + \dots, \end{aligned} \quad (47) $$

where $\delta_{ij}$ is the Kronecker-$\delta$ and we introduced the (complex-valued) moment multipoles

$$ (\omega_{\alpha})_{\ell m}^{p_j \dots p_l} = \frac{\langle A(p_j - \bar{p}_j) \dots (p_l - \bar{p}_l) \rangle_{\ell m}}{\bar{A}}. \quad (48) $$

In Eq. (47), we used the average of all the parameters across the sky to define the average SED amplitude and pivot<fnref n="10" />

$$ \bar{A} = \int \mathbb{P}(\mathbf{p}, \hat{\mathbf{n}}) d^N p \frac{d^2 \hat{\mathbf{n}}}{4\pi} \quad \text{and} \quad \bar{\mathbf{p}} = \frac{\int \mathbb{P}(\mathbf{p}, \hat{\mathbf{n}}) \mathbf{p} d^N p \frac{d^2 \hat{\mathbf{n}}}{4\pi}}{\int \mathbb{P}(\mathbf{p}, \hat{\mathbf{n}}) d^N p \frac{d^2 \hat{\mathbf{n}}}{4\pi}}. \quad (49) $$

This choice might be best suited for convergence (canceling the first order moment of the intensity expansion over the whole sky) when done in real space. This expansion is straightforward to generalize to angular power-spectra and cross-frequency power spectra (Mangilli et al. 2021) and is what has so far been used to describe B-modes signal.

### 5.1.2. In Polarization

The averaged polarized signal $\langle \mathcal{P}_v(\mathbf{p}(\hat{\mathbf{n}})) \rangle = \langle \mathcal{P}_v(\hat{\mathbf{n}}) \rangle$ is now a section of the spin-2 bundle on the $S^2$-manifold. It can be expanded on the orthogonal basis of the spin-2 weighted spherical harmonics $\pm 2 Y_{\ell m}$ (Newman & Penrose 1966; Goldberg et al. 1967). One can evaluate the spherical-harmonic coefficients of a general frequency dependent spinor field $\pm 2 \chi_v$ for $\ell \ge 2$ as<fnref n="11" />

$$ \pm 2(\chi_v)_{\ell m} = \int \pm 2 Y_{\ell m}^*(\hat{\mathbf{n}}) \chi_v(\hat{\mathbf{n}}) d^2 \hat{\mathbf{n}} \quad (50) $$

For convenience, we then also define the harmonic expansion of the line-of-sight average of the polarization field as

$$ \begin{aligned} \pm 2 \langle \mathcal{P}_v \rangle_{\ell m} &= \int \pm 2 Y_{\ell m}^*(\hat{\mathbf{n}}) \mathbb{P}(\mathbf{p}, \gamma, \hat{\mathbf{n}}) \mathcal{P}_v(\mathbf{p}, \gamma) d^N p \, d\gamma \, d^2 \hat{\mathbf{n}} \\ &= \int \pm 2 \mathbb{P}_{\ell m}(\mathbf{p}, \gamma) \mathcal{P}_v(\mathbf{p}, \gamma) d^N p \, d\gamma, \end{aligned} \quad (51) $$

where we define the harmonic coefficient of the parameter distribution function, $\mathbb{P}(\mathbf{p}, \gamma, \hat{\mathbf{n}})$, as $\pm 2 \mathbb{P}_{\ell m}(\mathbf{p}, \gamma)$. When applied to the spin-moment expansion of the polarization field, we then obtain

$$ \pm 2 \langle \mathcal{P}_v \rangle_{\ell m} = \sum_j \pm 2 (W_1^{p_j})_{\ell m} \partial_{\bar{p}_j} P_v(\bar{A}, \bar{\mathbf{p}}) $$

$$ \begin{aligned} &+ \frac{1}{2} \sum_{j,k} \pm 2 (W_2^{p_j p_k})_{\ell m} \partial_{\bar{p}_j} \partial_{\bar{p}_k} P_v(\bar{A}, \bar{\mathbf{p}}) \\ &+ \frac{1}{3!} \sum_{j,k,l} \pm 2 (W_3^{p_j p_k p_l})_{\ell m} \partial_{\bar{p}_j} \partial_{\bar{p}_k} \partial_{\bar{p}_l} P_v(\bar{A}, \bar{\mathbf{p}}) + \dots, \end{aligned} \quad (52) $$

with the spin-moment multipoles

$$ \pm 2 (W_{\alpha}^{p_j \dots p_l})_{\ell m} = \frac{\pm 2 \langle A e^{2i\gamma} (p_j - \bar{p}_j) \dots (p_l - \bar{p}_l) \rangle_{\ell m}}{\bar{A}}. \quad (53) $$

Since this expansion is not defined for $\ell \le 2$, one does not have any leading order (monopole) term. Nevertheless, the spectral pivot is again defined by the average over the full sky,

$$ \bar{A} = \int \mathbb{P}(\mathbf{p}, \gamma, \hat{\mathbf{n}}) d^N p \, d\gamma \frac{d^2 \hat{\mathbf{n}}}{4\pi} \quad (54a) $$

$$ \bar{\mathbf{p}} = \frac{\int \mathbb{P}(\mathbf{p}, \gamma, \hat{\mathbf{n}}) \mathbf{p} d^N p \, d\gamma \frac{d^2 \hat{\mathbf{n}}}{4\pi}}{\int \mathbb{P}(\mathbf{p}, \gamma, \hat{\mathbf{n}}) d^N p \frac{d^2 \hat{\mathbf{n}}}{4\pi}} \quad (54b) $$

which is the only physically motivated choice. It now has to include the average over $\gamma$.

### 5.1.3. E- and B-modes

From the harmonic coefficients of the polarized spinor, one can then obtain the E- and B-mode coefficients of the polarized SED as (Zaldarriaga & Seljak 1997; Kamionkowski et al. 1997):

$$ \langle \mathcal{P}_v \rangle_{\ell m}^E = -\frac{1}{2} [ +2 \langle \mathcal{P}_v \rangle_{\ell m} + -2 \langle \mathcal{P}_v^* \rangle_{\ell m} ] \quad (55a) $$

$$ \langle \mathcal{P}_v \rangle_{\ell m}^B = -\frac{1}{2i} [ +2 \langle \mathcal{P}_v \rangle_{\ell m} - -2 \langle \mathcal{P}_v^* \rangle_{\ell m} ] . \quad (55b) $$

One can then expand both $\mathcal{P}_v$ and its complex conjugate using Eq. (52) and insert them into Eq. (55). Going from the expansion of the previous section to the E- and B-modes adds no extra averaging effect. For each mode, we get a moment expansion with the same structure as Eq. (52) with spin moments:

$$ (W_{\alpha}^{p_j \dots p_l})_{\ell m}^E = -\frac{1}{2} [ +2 (W_{\alpha}^{p_j \dots p_l})_{\ell m} + -2 (W_{\alpha}^{p_j \dots p_l^*})_{\ell m} ] , \quad (56a) $$

$$ (W_{\alpha}^{p_j \dots p_l})_{\ell m}^B = -\frac{1}{2i} [ +2 (W_{\alpha}^{p_j \dots p_l})_{\ell m} - -2 (W_{\alpha}^{p_j \dots p_l^*})_{\ell m} ] . \quad (56b) $$

Since this new expansion is simply derived from the one of $\mathcal{P}_v$, it inherits its amplitude and pivot defined with the $\gamma$-weighted full sky averages. This justifies the approximation used in previous studies where the B-mode signal was treated as an intensity at the map (e.g., Remazeilles et al. 2021) or power-spectra level (e.g., Azzoni et al. 2021; Vacher et al. 2022). We leave a detailed discussion on the subtleties of E- and B- modes and power-spectra generalizations for future work.

## 5.2. Averages inside the beam and bandpass effects

Introducing the general instrumental transfer function $W(v, \gamma, \hat{\mathbf{n}})$ will add additional mixing between the lines of sight but also new nontrivial spectral dependencies from the mixing in frequency. In practice, this very general $W$ function could also account for some effects due to intensity to polarization leakage.

<footnote n="10">¹⁰ Note that in Mangilli et al. (2021) and its follow-up papers (e.g. Vacher et al. 2022), the existence of a scale dependent pivot $\bar{p}(\ell)$ in the harmonic space level has been assumed. Proving or discussing formally this assumption is left for future work.</footnote>

<footnote n="11">¹¹ The spin-weighted harmonics function $s_{Y_{\ell m}}$ are defined for $\ell \ge |s|$.</footnote>
The total polarized signal in a frequency band of width $\Delta v$ at average frequency $v_c$ and within the spatial support of the beam, $\Omega$, centered in an average line of sight $\hat{n}_c$ then becomes:

$$ \langle \mathcal{P}_v \rangle_W = \int_{\Omega} \mathbb{P}(\mathbf{p}, \gamma, \hat{n}) W(v, \gamma, \hat{n}) \hat{P}_v(\mathbf{p}) e^{2i\gamma} d^N p \, dy \, d^2 \hat{n} \, dv. \quad (57) $$

Here, $\langle \mathcal{P}_v \rangle_W$ is now a function of $v_c$ and the corresponding moments and pivots in the average direction $\hat{n}_c$.

Due to the $v$ integration in the above expression, it is not possible to simply factor the spectral shapes in the Taylor expansion of $\hat{P}_v$ out of the integral as we did before (Chluba et al. 2017). The band-pass function of the instrument therefore couples different spatial regions with different spectral forms. One can still generalize the spin-moment expansion to account for this extra averaging. The new expansion becomes:

$$ \begin{aligned} \langle \mathcal{P}_v \rangle_W ={}& \mathcal{W}_0 \langle \hat{P}_v(\bar{\mathbf{A}}, \bar{\mathbf{p}}) \rangle_W + \sum_j \mathcal{W}_1^{p_j} \langle \partial_{\bar{\mathbf{p}}_j} \hat{P}_v(\bar{\mathbf{A}}, \bar{\mathbf{p}}) \rangle_W \\ & + \frac{1}{2} \sum_{j,k} \mathcal{W}_2^{p_j p_k} \langle \partial_{\bar{\mathbf{p}}_j} \partial_{\bar{\mathbf{p}}_k} \hat{P}_v(\bar{\mathbf{A}}, \bar{\mathbf{p}}) \rangle_W \\ & + \frac{1}{\alpha!} \sum_{j, \dots, l} \mathcal{W}_{\alpha}^{p_j \dots p_l} \langle \partial_{\bar{\mathbf{p}}_j} \dots \partial_{\bar{\mathbf{p}}_l} \hat{P}_v(\bar{\mathbf{A}}, \bar{\mathbf{p}}) \rangle_W + \dots, \end{aligned} \quad (58) $$

with the new averages defined as

$$ \langle X \rangle_W = \int \mathbb{P}(\mathbf{p}, \gamma, \hat{n}) W(v, \gamma, \hat{n}) X d^N p \, dy \, d^2 \hat{n} \, dv. \quad (59) $$

The pivot is chosen to be the same as in intensity:

$$ \bar{p}_j = \frac{\langle A p_j \partial_{\bar{p}_j} P_v(\bar{p}) \rangle_W}{\langle A \partial_{\bar{p}_j} P_v(\bar{p}) \rangle_W} \quad (60) $$

and the spin moments become:

$$ \mathcal{W}_{\alpha}^{p_j \dots p_l} = \frac{\langle A e^{2i\gamma} (p_j - \bar{p}_j) \dots (p_l - \bar{p}_l) \partial_{\bar{p}_{j \dots p_l}} \hat{P}_v(\mathbf{p}) \rangle_W}{\langle A \partial_{\bar{p}_{j \dots p_l}} \hat{P}_v(\mathbf{p}) \rangle_W}. \quad (61) $$

These expressions assume that the various SED derivative averages do not vanish, however, the conclusions are not affected.

Generally, one cannot disentangle the mixing due to the physics of the source and the one due to the instrumental response, and in practice this moment expansion can become awfully complicated. A frequent assumption is to introduce a factorization of $W$ as a band-pass term $F$ times a spatial polarized beam shape $\mathfrak{B}$ as $W(v, \hat{n}) = \mathfrak{B}(\gamma, \hat{n})F(v)$ see e.g. Planck Collaboration IX (2014). In this case, it is possible to split the integral of Eq. (57) to factorize the frequency dependent terms, allowing the vanishing of the spectral averages such that the expressions for the pivot and the spin moments are similar to the ones derived for the single line-of-sight case, extending only the averages to the multiple lines of sight included in $\mathfrak{B}$. The result remains nontrivial since the spectral dependence of the moments need to be computed through the potentially complicated integral of $F(v) \partial_{p_j \dots p_k} \hat{P}_v$ over $\Delta v$. The expansion can also be different in each band of the instrument, since they can have a different instrumental response $W$. Treating these cases in more detail is beyond the scope of this work.

# 6. Additional aspects

## 6.1. Effect of voxel-level SED variations

As stressed in Sect. 3.3, in most of this work we considered that inside each voxel, the net linear polarization have a single SED of which $Q$ and $U$ are the projections. The average signal then inherits a frequency-dependent polarization angle solely from polarized mixing across voxels. Among our examples in Sect. 4, the only exception was the mixing of blackbody spectra, where the polarization degree was caused by variations of the main spectral parameter, the blackbody temperature.

Can we extend the moment formalism to include voxel-level SED mixing caused by variations of the spectral parameters? Let us use the power-law SED as the example. In addition to the weight parameter variations in the different directions, we now also have variations of the spectral index, $\beta$. Considering only $Q_v$, from the discussion in Sect. 3.3 we then have

$$ Q_v = \frac{A_{\parallel} \hat{P}_v^{\text{PL}}(\beta_{\parallel}) - A_{\perp} \hat{P}_v^{\text{PL}}(\beta_{\perp})}{2} = \frac{(A_{\parallel} - A_{\perp}) \hat{P}_v^{\text{PL}}(\beta_{\parallel}) + \hat{P}_v^{\text{PL}}(\beta_{\perp})}{2} + \frac{A_{\parallel} + A_{\perp}}{2} \frac{\hat{P}_v^{\text{PL}}(\beta_{\parallel}) - \hat{P}_v^{\text{PL}}(\beta_{\perp})}{2}. \quad (62) $$

A similar expression follows for $U_v$, with the relevant parameters $A_\times, \beta_\times, A_\odot, \beta_\odot$. Let us again expand the fundamental SEDs around some average index $\bar{\beta}$ inside each voxel<fnref n="12" />. This then yields

$$ \begin{aligned} Q_v ={}& \frac{(A_{\parallel} - A_{\perp})}{2} \hat{P}_v^{\text{PL}}(\bar{\beta}) \left[ 1 + \sum_{k=1}^{\infty} \frac{\Delta\beta_{\parallel}^k + \Delta\beta_{\perp}^k}{2k!} \ln^k(v/v_0) \right] \\ & + \frac{A_{\parallel} + A_{\perp}}{2} \hat{P}_v^{\text{PL}}(\bar{\beta}) \left[ \sum_{k=1}^{\infty} \frac{\Delta\beta_{\parallel}^k - \Delta\beta_{\perp}^k}{2k!} \ln^k(v/v_0) \right] \\ ={}& \hat{P}_v^{\text{PL}}(\bar{\beta}) \left( \frac{(A_{\parallel} - A_{\perp})}{2} + \sum_{k=1}^{\infty} \frac{A_{\parallel}\Delta\beta_{\parallel}^k - A_{\perp}\Delta\beta_{\perp}^k}{2k!} \ln^k\left(\frac{v}{v_0}\right) \right) \end{aligned} \quad (63) $$

with $\Delta\beta_k = \beta_k - \bar{\beta}$. We already considered the first term for which spectral index variations between voxels lead to the moment expansion. In addition to these voxel-to-voxel variations one now obtains additional terms that are related to variations within the voxel. All these will cause a redefinition of the SED parameter distribution, however, no new spectral shapes are introduced and hence the moment method equally describes both averaging processes, after the moments are reinterpreted. However, to explicitly make the link to the underlying spectral parameter distributions along the line of sight and within each voxel involves a more complicated description, which we do not provide here.

We close by remarking that, if even the fundamental SEDs in each of the directions differ, then one should simply perform two independent moment expansions accounting for the independent types of SEDS individually. This will quickly increase the number of moments that are required to describe the complexity of the polarized field, but it should nevertheless work even if the physical properties are not independent. The main hope then is that line-of-sight averaging effects reduce the dimensionality of the problem to a manageable level. A more detailed discussion is, however, beyond the scope of this work.

<footnote n="12">¹² This in fact is $\bar{\beta}_v = (A_{\parallel}\beta_{\parallel} + A_{\perp}\beta_{\perp} + A_{\times}\beta_{\times} + A_{\odot}\beta_{\odot})/(A_{\parallel} + A_{\perp} + A_{\times} + A_{\odot})$.</footnote>
## 6.2. Frequency dependence of the polarization angle without mixing: The case of Faraday rotation

In the presence of magnetic fields, photons experience Faraday rotation while they propagate through the interstellar medium (e.g., Heald 2015). The polarization angle for one voxel at affine parameter s then changes as

$$ \gamma_v \approx \gamma + \Gamma f_v, \quad (64) $$

where $f_v = v^{-2}$ and $\Gamma(s) \propto \int_0^s n_e(s')B(s')ds'$ is an integral along the line of sight, with the electron density $n_e(s)$ and the component of the magnetic field in the direction of propagation, $B(s)$. Averaging over various emission points experiencing Faraday rotation adds some extra complications to the moment expansion in a way that is truly unique to polarization. For simplification let us consider the average over a single line of sight. The averaged spinor now becomes:

$$ \langle \mathcal{P}_v \rangle_{\text{FR}} = \int \mathbb{P}(\mathbf{p}, \gamma, \Gamma, \hat{\mathbf{n}}) \hat{P}_v(\mathbf{p}) e^{2i\gamma} e^{2i\Gamma f_v} d^N p \, d\gamma \, d\Gamma, \quad (65) $$

where we added the level of Faraday rotation as another parameter to the distribution. Taylor expanding the polarized intensity, one can define the spin-moment expansion using the moments

$$ \mathcal{W}_{\alpha, \text{FR}}^{p_j \cdots p_l} = \frac{\langle A e^{2i\gamma} e^{2i\Gamma f_v} (p_j - \bar{p}_j) \cdots (p_l - \bar{p}_l) \rangle_{\text{FR}}}{\bar{A}_{\text{FR}}}. \quad (66) $$

with $\bar{A}_{\text{FR}} = \langle A \rangle_{\text{FR}}$. The spin moments now become highly nontrivial and generally frequency dependent.

As discussed in Sect. 3.4.3, under simplifying assumptions one can assume that the emission processes themselves (quantified by the spectral parameters $\mathbf{p}$) are de-correlated from the Faraday rotation experienced by light on its way along the line of sight. This would allow us to write $\mathbb{P}(\mathbf{p}, \gamma, \Gamma, \hat{\mathbf{n}}) \approx \mathbb{P}(\mathbf{p}, \gamma, \hat{\mathbf{n}}) \mathbb{P}(\Gamma, \hat{\mathbf{n}})$. Physically, this is indeed well-motivated unless the emission process in one voxel knows about the structure of the magnetic field in another (more distant) voxel. The Faraday rotation and the SED averaging in the moments then becomes separable:

$$ \begin{aligned} \mathcal{W}_{\alpha, \text{FR}}^{p_j \cdots p_l} &\approx \frac{\int \mathbb{P}(\mathbf{p}, \gamma, \hat{\mathbf{n}}) \hat{P}_v(\mathbf{p}) e^{2i\gamma} d^N p \, d\gamma}{\int \mathbb{P}(\mathbf{p}, \gamma, \hat{\mathbf{n}}) d^N p \, d\gamma} \frac{\int \mathbb{P}(\Gamma, \hat{\mathbf{n}}) e^{2i\Gamma f_v} d\Gamma}{\int \mathbb{P}(\Gamma, \hat{\mathbf{n}}) d\Gamma} \\ &= \mathcal{W}_{\alpha}^{p_j \cdots p_l} \frac{\int \mathbb{P}(\Gamma, \hat{\mathbf{n}}) e^{2i\Gamma f_v} d\Gamma}{\int \mathbb{P}(\Gamma, \hat{\mathbf{n}}) d\Gamma}. \end{aligned} \quad (67) $$

This expression shows that no new spectral mixing occurs in this case, as all moments are multiplied by the same frequency-dependent factor. The average Faraday rotation coefficient is simply given by

$$ \bar{\Gamma}(\hat{\mathbf{n}}) = \frac{\langle A\Gamma \rangle_{\text{FR}}}{\langle A \rangle_{\text{FR}}} \approx \frac{\int \mathbb{P}(\Gamma, \hat{\mathbf{n}}) \Gamma \, d\Gamma}{\int \mathbb{P}(\Gamma, \hat{\mathbf{n}}) \, d\Gamma}. \quad (68) $$

Using this as a pivot, we can then write series

$$ \frac{\langle A e^{2i\Gamma f_v} \rangle_{\text{FR}}}{\bar{A}_{\text{FR}}} = e^{2i\bar{\Gamma} f_v} \left[ 1 + \sum_{k=1}^{\infty} \frac{(2i \hat{f}_v)^k}{k!} \frac{\langle A(\Gamma - \bar{\Gamma})^k \rangle_{\text{FR}}}{\bar{A}_{\text{FR}}} \right], \quad (69) $$

illustrating how a complicated frequency structure can be created from higher order moments of $\Gamma - \bar{\Gamma}$. Given that these

depend on powers of $f_v$, one can in principle determine these moments observationally. Finally, if the simple factorization of the parameter distribution function is not possible, one can write

$$ \mathcal{W}_{\alpha, \text{FR}}^{p_j \cdots p_l} = e^{2i\bar{\Gamma} f_v} \left[ \mathcal{W}_{\alpha}^{p_j \cdots p_l} + \sum_{k=1}^{\infty} \frac{(2i \hat{f}_v)^k}{k!} \frac{\langle A e^{2i\gamma} (p_j - \bar{p}_j) \cdots (p_l - \bar{p}_l) (\Gamma - \bar{\Gamma})^k \rangle_{\text{FR}}}{\bar{A}_{\text{FR}}} \right]. \quad (70) $$

This generally introduces a complicated voxel SED reweighting by factors of $(f_v)^k$ and hence new SED shapes that in principle can all be separately accounted for. As is understood, the number of variables quickly becomes unmanageable unless highly perturbative situations are encountered. A more detailed discussion is, however, beyond the scope of this work.

## 7. Conclusion

In the present work, we introduced the spin moments, the natural generalization of the intensity moment expansion introduced in Chluba et al. (2017) to polarized signals. We developed the formalism from basic principles, showing that the moments are promoted to spin-2 complex coefficients that can be expressed in term of the SEDs parameter distribution (see e.g. Eq. (15)).

Thinking about the spin moments in the form of spinors allows us to treat several subtleties due to the geometrical nature of polarization. A clear interpretation of the polarized mixing distortions arises, as we show that a rotation of the spinor with frequency is naturally induced from the distribution of spectral parameters and polarization angles. We demonstrate that, no general pivot can be defined ensuring the vanishing of the first order in the non-perturbative regime. In the perturbative regime, however, such a pivot can be found and hides interesting physics. Correcting for this pivot naturally gives rise, in addition to a shift of the spectral parameters, to a frequency-dependent rotation of the polarization angle of predictable spectral dependence. It can also bring extra modulations to the polarized intensity, for example in the case of gray-bodies.

We explored scenarios of increasing complexity, considering several canonical SEDs examples of first importance for astrophysics. Even when dealing with highly complex signals along the line of sight, we showed that the use of spin moments allows us to model the distorted polarized intensity and the frequency dependent rotation of the polarization angle, including only a few terms in the expansion.

We also discussed the effect of more complex averaging processes such as spherical harmonics mixing and instrumental mixing as well as non trivial polarization specific situations like SED variations at the voxel-level and Faraday rotation effect. The spin-moment formalism still applies in all these scenarios. In these increasingly complex cases, however, the interpretation of the moment coefficients becomes blurry and the various mixing processes cannot be expected to be properly disentangled. Doing so, we also rederived formally the expressions used in previous works, treating B-mode signal as an intensity at the map (Remazeilles et al. 2016) and power-spectra levels (Azzoni et al. 2021; Vacher et al. 2022).

This work opens the door to several follow-up applications. On the theoretical side, several discussions remain open on the details of the generalization of the formalism at the power-spectrum level, especially for E- and B-mode applications. One could also think of links with cosmic birefringence, which
recently received increased attention (see e.g. Diego-Palazuelos et al. 2022). More broadly, a lot of room is left for application of the spin moments as we defined them. One can think first to component separation where the spin moments could be competitive to model foreground distortions at the map level on large scales where a lot of averaging is done, calling for a comparative study with other pixel-based methods of component separation. Interesting questions related to Galactic physics have also to be addressed with the spin moments, such as the possibility to disentangle dust composition and Galactic magnetic field effects or tackle the Faraday rotation. One can also think of applications to SZ effect or spectral distortions, all topics that are left for future explorations.

**Acknowledgements.** L.V. would like to thanks François Boulanger for numerous and fruitful discussions. J.C. was supported by the Royal Society as a Royal Society University Research Fellow at the University of Manchester, UK (No. URF/R/191023). This work was also supported by the ERC Consolidator Grant CMBSPEC (No. 725456) as part of the European Union's Horizon 2020 research and innovation program.

## References

Aiola, S., Calabrese, E., Maurin, L., et al. 2020, J. Cosmol. Astropart. Phys., 2020, 047

Azzoni, S., Abitbol, M., Alonso, D., et al. 2021, J. Cosmol. Astropart. Phys., 2021, 047

Bennett, C. L., et al. 2013, ApJS, 208, 20

Brout, R., Englert, F., & Gunzig, E. 1978, Ann. Phys., 115, 78

Carlstrom, J. E., Holder, G. P., & Reese, E. D. 2002, ARA&A, 40, 643

Chluba, J., & Sunyaev, R. A. 2004, A&A, 424, 389

Chluba, J., Switzer, E., Nelson, K., & Nagai, D. 2013, MNRAS, 430, 3054

Chluba, J., Dai, L., Grin, D., Amin, M. A., & Kamionkowski, M. 2015, MNRAS, 446, 2871

Chluba, J., Hill, J. C., & Abitbol, M. H. 2017, MNRAS, 472, 1195

Chluba, J., Abitbol, M. H., Aghanim, N., et al. 2021, Exp. Astron., 51, 1515

CMB-S4 Collaboration. 2019, ArXiv e-prints [arXiv:1907.04473]

Delabrouille, J., Abitbol, M., Aghanim, N., et al. 2021, Exp. Astron., 51, 1471

Diego-Palazuelos, P., Eskilt, J. R., Minami, Y., et al. 2022, Phys. Rev. Lett., 128, 091302

Draine, B. T., & Fraisse, A. A. 2009, ApJ, 696, 1

Ferrière, K. M. 2001, Rev. Mod. Phys., 73, 1031

Goldberg, J. N., Macfarlane, A. J., Newman, E. T., Rohrlich, F., & Sudarshan, E. C. G. 1967, J. Math. Phys., 8, 2155

Guth, A. H. 1981, Phys. Rev. D, 23, 347

Heald, G. 2015, Astrophys. Space Sci. Lib., 407, 41

Hoseinpour, A., Zarei, M., Orlando, G., Bartolo, N., & Matarrese, S. 2020, Phys. Rev. D, 102, 063501

Hutton, S., Ferreras, I., & Yershov, V. 2015, MNRAS, 452, 1412

Ichiki, K., Kanai, H., Katayama, N., & Komatsu, E. 2019, Prog. Theor. Exp. Phys., 2019, 033E01

Inomata, K., & Kamionkowski, M. 2019, Phys. Rev. D, 99, 043501

Jaffe, T. R., Ferrière, K. M., Banday, A. J., et al. 2013, MNRAS, 431, 683

Kamionkowski, M., Kosowsky, A., & Stebbins, A. 1997, Phys. Rev. D, 55, 7368

Kogut, A., Fixsen, D. J., Chuss, D. T., et al. 2011, J. Cosmol. Astropart. Phys., 2011, 025

LiteBIRD Collaboration. 2022, PTEP, 11443, 114432F

Mangilli, A., Aumont, J., Rotti, A., et al. 2021, A&A, 647, A52

Montero-Camacho, P., & Hirata, C. M. 2018, J. Cosmol. Astropart. Phys., 2018, 040

Mroczkowski, T., Nagai, D., Basu, K., et al. 2019, Space Sci. Rev., 215, 17

Newman, E. T., & Penrose, R. 1966, J. Math. Phys., 7, 863

Newville, M., Stensitzki, T., Allen, D. B., et al. 2016, Astrophysics Source Code Library [record ascl:1606.014]

Pelgrims, V., Clark, S. E., Hensley, B. S., et al. 2021, A&A, 647, A16

PICO Collaboration. 2019, BAAS, 51, 194

Planck Collaboration IX. 2014, A&A, 571, A9

Planck Collaboration XVII. 2014, A&A, 566, A55

Planck Collaboration XXII. 2015, A&A, 576, A107

Planck Collaboration I. 2016, A&A, 594, A1

Planck Collaboration Int. L. 2017, A&A, 599, A51

Planck Collaboration IV. 2020, A&A, 641, A4

Planck Collaboration XI. 2020, A&A, 641, A11

Remazeilles, M., Dickinson, C., Eriksen, H. K. K., & Wehus, I. K. 2016, MNRAS, 458, 2032

Remazeilles, M., Rotti, A., & Chluba, J. 2021, MNRAS, 503, 2478

Rotti, A., & Chluba, J. 2021, MNRAS, 500, 976

Sayre, J. T., Reichardt, C. L., Henning, J. W., et al. 2020, Phys. Rev. D, 101, 122003

Schlafly, E. F., Meisner, A. M., Stutz, A. M., et al. 2016, ApJ, 821, 78

Starobinsky, A. A. 1980, Phys. Lett. B, 91, 99

Stompor, R., & Efstathiou, G. 1999, MNRAS, 302, 735

Tassis, K., & Pavlidou, V. 2015, MNRAS, 451, L90

The Simons Observatory collaboration 2019, BAAS, 51, 147

Thorne, B., Dunkley, J., Alonso, D., & Næss, S. 2017, MNRAS, 469, 2821

Vacher, L., Aumont, J., Montier, L., et al. 2022, A&A, 660, A111

Wise, J. H. 2019, Contemp. Phys., 60, 145

Ysard, N., Abergel, A., Ristorcelli, I., et al. 2013, A&A, 559, A133

Zaldarriaga, M., & Seljak, U. 1997, Phys. Rev. D, 55, 1830
# Appendix A: Alternative approaches

For astrophysical applications, it is common to normalize the modified blackbody SED in every pixel as in Thorne et al. (2017)

$$ \hat{P}_{\nu}^{\text{mBB}}(\hat{\mathbf{n}}) = \frac{\hat{P}_{\nu}^{\text{mBB}}(T(\hat{\mathbf{n}}), \beta(\hat{\mathbf{n}}))}{\hat{P}_{\nu_0}^{\text{mBB}}(T(\hat{\mathbf{n}}), \beta(\hat{\mathbf{n}}))}. \qquad (A.1) $$

Since here both $\beta$ and $T$ are treated as spatially varying parameters, this is no longer a constant SED renormalization once line-of-sight effects are included. This choice therefore complicated the expansion, leading to rescaling of the moments (e.g., see Eq.(4) of Vacher et al. 2022), without any physical meaning or insight being added. At the pixel level, this complication can be avoided by setting

$$ \hat{P}_{\nu}^{\text{mBB}*}(\hat{\mathbf{n}}) = \frac{\hat{P}_{\nu}^{\text{mBB}}(T(\hat{\mathbf{n}}), \beta(\hat{\mathbf{n}}))}{\hat{P}_{\nu_0}^{\text{mBB}}(\bar{T}(\hat{\mathbf{n}}), \bar{\beta}(\hat{\mathbf{n}}))}, \qquad (A.2) $$

which simply takes out the leading order term in the expansion in Eq. (43). However, when extending to applications on the full sky (or when averaging pixels), pixel to pixel variations of $\bar{\beta}$ and $\bar{T}$ will come into play. In this case, one should better normalize using one constant $\bar{\beta}$ and $\bar{T}$ across the sky to avoid additional complications from the variation of the normalization.

# Appendix B: Choice of weights

As discussed already in Chluba et al. (2017), when deriving the moment expansion, one can always reexpress the SEDs with another choice of units<fnref n="13" /> or with respect to an alternative choice of spectral variables. Such a choice should be motivated by the physics and the numerical behavior of the problem considered. Ultimately, doing so will simply be equivalent to a re-scaling of the weight coefficients appearing in the moment expansion. Such a change can have a significant impact on the convergence rate of the expansion and the interpretations of the moments coefficients but is mathematically equivalent.

For example, in the case of the gray-body, one could equivalently define the new weights $A'(T) = AT^3$ and expand around $T$ instead of $\beta_{GB}$. The fundamental SED then reads

$$ \hat{P}_{\nu}^{\text{GB}} = A' \left( \frac{\bar{T}}{T} \right)^3 \frac{x^3}{e^{\frac{x}{T}} - 1}. \qquad (B.1) $$

The moment expansion will change accordingly as displayed in Eq. (43) of Chluba et al. (2017). This allows us to interpret the spin moments and pivot directly in term of temperatures:

$$ \bar{A}' = \langle A' \rangle = \langle AT^3 \rangle \qquad (B.2a) $$

$$ \bar{T} = \frac{\langle A'T \rangle}{\bar{A}'} = \frac{\langle AT^4 \rangle}{\langle AT^3 \rangle} \qquad (B.2b) $$

$$ W_{\alpha}^{T} = \frac{\langle A'(T - \bar{T})^{\alpha} e^{2i\gamma} \rangle}{\bar{A}'} = \frac{\langle AT^{3} (T - \bar{T})^{\alpha} e^{2i\gamma} \rangle}{\langle AT^{3} \rangle}. \qquad (B.2c) $$

The exact same reasoning applies to a change of the power law's reference frequency $\nu_0 \to \nu'_0$ with the simple reweighting $A' = A(\nu'_0/\nu_0)^\beta$. However, in this case, no new SED derivatives are created, so the two expansions are identical.

<footnote n="13">¹³ Some changes of units as e.g. Jy.sr⁻¹ → µK$_{CMB}$ are slightly subtler and introduce new frequency-dependent terms.</footnote>