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

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 parameters 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 parameters 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.


Introduction
A significant international effort has been undertaken to deploy large-scale surveys of the comsic 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 τ (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 A5, page 1 of 16 Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0),which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.This article is published in open access under the Subscribe-to-Open model.Subscribe to A&A to support open access publication.A&A 669, A5 (2023) 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. 2016Remazeilles et al. , 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/Uintensity 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 selfconsistent 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.

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 n, which is described by an affine parameter s.Every voxel emits with an SED 1 : where Îν is referred to as the fundamental SED with N spectral parameters, p(s) = {p 1 (s), p 2 (s), . . ., p N (s)}.The amplitude or weight parameter, A(s), determines the relative contribution of each voxel to the total intensity 2 .The resulting total SED is given by an average along the line of sight, which we shall denote by ⟨. . .⟩.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 n (e.g., Chluba et al. 2017;Rotti & Chluba 2021): In the second definition, we introduced the distribution P( p, n) of the spectral parameters, p, along the fixed line of sight n, with the relative weights absorbed into the distribution itself.We note that the distribution P( p, 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 In order to provide a perturbative model of the average SED, the spectral dependence Îν in Eq. ( 2) can be expanded into a Taylor series with respect to p around the pivot p as: Here, we used the shorthand notation ∂ p j X( p) ≡ ∂X( p)/∂ p j .The pivot value p around which the series is carried out can be fixed by asking for the first term of the expansion to vanish upon averaging: ⟨A j (p j − p j )⟩ = 0.This minimizes the required terms in the Taylor series and leads to: 1 Hereafter, we use the shorthand notation for frequency dependent quantities X ν ≡ X(ν).
2 One simple example is the power law: I ν (A, β) = A (ν/ν 0 ) β where p = {β} has dimension one (N = 1) and Îν (β) = (ν/ν 0 ) β .In this, ν 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).
with α 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: Here, all the ω p i 1 are zero when using the value for 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 p and then correcting the solution by the values of ω p i 1 ≡ ∆p i using p j → p′ j = p j + ∆p j .This iterative process assumes that the moments are perturbative and convergence can be obtained with a finite number of terms.

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.

General introduction to polarized SED
A polarized signal is fully described by the four real Stokes parameters I ν , Q ν , U ν , V ν , all of them being frequency dependent quantities.As before, I ν describes the total (i.e., unpolarized + polarized) intensity, while the pair (Q ν , U ν ) and V ν respectively quantify the linearly and the circularly polarized part of the photon field.Both I ν and 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 ν with the moment expansion could introduce extra subtleties which go beyond the scope of this work.Henceforth, we assume V ν = 0 3 .
On the other hand, the Q ν and U ν 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 P ν : 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).
The spinor's modulus, P ν , is a real positive function called the linear polarization intensity and its argument defines the polarization angle γ ν : As a spin-2 quantity, when the frame in which Q ν and U ν are defined (e.g., by modifying the directions of the polarizers) is rotated by a right handed rotation around the n direction by an angle θ, P ν transforms as Zaldarriaga & Seljak (1997): Note that, unless stated otherwise, henceforth we use calligraphic variables for complex quantities (e.g., P ν , W, . . . ) and italic font for real quantities (e.g., P ν , A, Q ν , U ν , ω, . . .).

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 γ, 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 ′ ν P ν and a frequency dependent γ ν .
This consequence of polarized mixing is illustrated in Fig. 1 for a sum of two power laws A i (ν/ν 0 ) β i e 2iγ i with ν 0 = 300 GHz, A 1 = 2 Jy sr −1 , A 2 = 1 Jy sr −1 , β 1 = 1.8, β 2 = 1.2, 2γ 1 = 10 • and 2γ 2 = 80 • .One clearly sees that the resulting spinor P ν rotates in the complex plane with frequency.Modeling both the distorted SED P ′ and the frequency dependence of γ ν in a physically-motivated fashion is nontrivial, but can be achieved when generalizing the moment expansion to polarization.

Understanding the link to intensity moment expansion
To generalize the intensity moment expansion to polarization, we have to discuss how Q ν and U ν 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 ν,∥ and I ν,⊥ , which are orthogonal to each other, and I ν,× and I ν,⊗ , which are also orthogonal to each other but rotated by 45 • relative to the previous system.The total intensity I ν (polarized + unpolarized), and Stokes Q ν and U ν are then given by The color labels the frequency between 1 GHz (dark red) and 500 GHz (dark blue).Lower panel: polarization spinor P ν (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γ ν .The length of the colored bars and the lines of constant radius represent values of log 10 (P ν ).
This shows that both Q ν and U ν 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 ν,∥ , I ν,⊥ , I ν,× and I ν,⊗ as the cummulative 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 GB ν = A B ν (T ), where B ν (T ) is a blackbody spectrum.The linearly polarized radiation of a single voxel is then given by Starting with this voxel SED renders the problem quite complicated.For example, just considering the SED of Q ν , we can write 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 B ν (T ).In this case, the fundamental voxel SED is given by PGB ν ≈ B ν (T ), such that the single polarization state can be characterized by P GB ν ≈ A e 2iγ B ν (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.

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 γ ν 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, n, is then given by Eq. (2).For polarization, this is equivalent to performing a Taylor expansion of the spinor's modulus4 P ν = |P ν | with respect to the spectral parameters at each fixed polarization angle γ.However, we cannot use a perturbative approach to average over the polarization angles γ 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 Here, we again used P ν (A, p) = A Pν ( p), as for intensity.In analogy to the intensity moment expansion, one then finds which depends on the pivot p, as we specify in Sect.3.4.2.Since A e 2iγ can vanish, we cannot simply factor it out of the expressions.Instead, like for the intensity moments, we will again use Ā = ⟨|A e 2iγ |⟩ ≡ ⟨A⟩ as a way to normalize the distributions.This allows us to define the spin moments 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, Ω p j ...p l α and γ p j ...p l α .The latter defines average directions of polarization states associated with each of the moments of the SED.While considering the pair (Ω α , γ α ) or the real and imaginary parts of W α 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: which in this form can be interpreted as the sum of multiple SEDs with well-defined polarization states.It is this sum of welldefined single polarization states (i.e., defined by the complexvalued 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, p.However, the overall normalization Ā does not independently contribute, but was merely chosen to scale the moments.As such, it cannot be independently estimated, and only the values of Ā W α actually matter.

Average polarization angle
Since W 0 = Ω 0 e 2iγ 0 = A e 2iγ / Ā can generally vanish, there is no longer a trivially defined average polarization angle.In particular when W 0 ≈ 0, the average polarization angle can be fully determined by the higher order terms in Eq. ( 16) and also generally becomes frequency-dependent.
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 γ along the line of sight in terms of the distribution, P( p, γ, n).This then results in as the average polarization angle.This angle can also be used as a pivot when expanding the polarization state: Using this in Eq. ( 15), have The first term in the sum (i.e., k = 0), is the only non vanishing contribution if the distributions of γ and 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 γ and p.However, in terms of distinguishable parameters, only the total moments, W p j ...p l α , can really be constrained.

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 However, since ⟨A e 2iγ ⟩ can vanish, in general this cannot be a meaningful choice.Above, we already defined Ā = ⟨|A e 2iγ |⟩.In a similar manner, we can introduce the SED pivots as which is equivalent to the definition for the intensity moments.
Physically, this means that we disregard the geometrical properties of P ν and simply treat its modulus as an intensity.For our power-law example in Sect.3.4.1,this means β , 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, Ω 0 = |⟨Ae 2iγ ⟩| ≫ 0, one can safely choose the complex pivot While the spectral parameters 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 p can be interpreted as real correction of p, its complex part gives rise to the first order frequency dependence of the polarization angle γ ν and can add some spectral modulation to the polarized intensity.

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 5 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 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 γ and p is required to have a spectral dependence of the polarization angle.
If the angle distribution is Gaussian with average angle γ( n) and width σ γ ( n), then one finds 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 ∆γ = γ − γ, (see e.g.Eq. ( 18)) is unlikely to converge, but, as stressed already, does not add any new insight anyways.

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 χ 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 M ν and the simulated data signal is noted P S ν .We add Gaussian noise N ν to the simulation, with zero mean and standard deviation σ = σ Q + iσ U .The values of σ is chosen such that the signal to noise ratios 5 In doing so, we can use the normalizations P(p, n)d N p = Ā and

Q S
ν /Q M ν and U S ν /U M ν are constants over the whole frequency range (chosen arbitrarily to be 1 × 10 −5 )6 .
The χ 2 to minimize is given by χ The signal is considered over a frequency range going from 1 GHz to ν max in intervals of 1 GHz.The choice of ν max will depend on the example considered.We introduce the shorthand notation 'O(α)' to refer to the fit of the spin-moment expansion including all the terms up to order α.'O(0)' is the leading order/canonical SED.Two distinct routines are developed, fitting either the pair (Re(W p α ), Im(W p α )) or the pair (Ω p α , γ 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.

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 where δ(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 γ k .Inserting this into the definitions of the moments and pivots given in the previous section we trivially find the exact average Using the polarization moment expansion, we automatically have the normalization, pivot and complex-valued moments as where the ratios A k / Ā determine the probabilities to find p k and γ 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.

Power laws
As a first example of astrophysical relevance, we consider the simple case of power-law SEDs: The polarization state can then be characterized by P PL ν ≈ Ae 2iγ P PL ν .The only spectral parameter relevant for the moment expansion is the spectral index p = (β), normalized at a reference frequency ν 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 ν max = 150 GHz below which the synchrotron emission is dominant and ν 0 = 23 GHz as the WMAP frequency band (Bennett et al 2013) k , the spin-moment expansion in Eq. ( 16) can then be expressed as The choice of the reference frequency ν 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: ∂ ν P S ν = 0 or ∂ ν γ 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.

Hands-on example: Two power laws
Consider two power laws (M = 2) with different spectral indices (β 1 , β 2 ) and polarization angle (γ 1 , γ 2 ) along the same line of sight.The exact solution then reads Carrying out the intensity moment expansion of the two individual power laws with respect to their spectral indices, we obtain These expressions can be trivially extended to M power laws after extending the sums to M parameters A k , γ k and β 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 γ 1 = γ 2 = γ and β 1 β 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, ⟨P PL ν ⟩ = e 2iγ ⟨P PL ν (A, β)⟩, trivially describing the effect of spectral mixing only.If on the other  hand β 1 = β 2 = β and γ 1 γ 2 , we naturally have W β α α = 0 and everything is described by W 0 with a fixed SED.To obtain nontrivial consequences of polarized mixing, both γ k and β 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 , β 1 = −2.8,β 2 = −3.6,2γ 1 = 10 • and 2γ 2 = 80 • , implying Ā = 3 Jy sr −1 , β ≈ −3.06 and γ ≈ 16.6 • .In Fig. 2, the A5, page 7 of 16 A&A 669, A5 (2023) modulus P ν and the argument γ ν 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 β ≈ −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 Ā, β and W β α 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 γ ν .

Extreme cases and perturbative regime
To gain further insight, let us just consider the first order terms of the expansion: If W 0 ≈ 0, we indeed find the situation where we have a polarization angle fully determined by W β 1 , with a sign-flip at ν = ν 0 .The polarization SED is then determined by the first β 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 ν and −Q ν ).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 , β 1 = 1, β 2 = 2, 2γ 1 = 180 • and 2γ 2 = 0.1 • .One can see that both γ ν and P ν do not behave as smooth functions at the breaking point ν ≃ ν 0 , where one power law abruptly takes over the other one.Here, γ ν changes very rapidly from γ 1 to γ 2 .This is, however, not a problem for the spin-moment expansion, which allows us to recover P ν correctly.As mentioned above, one recovers a very small best-fit value for W 0 and W β 1 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 |W 0 | ≫ |W β 1 | > 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 With 7 W 0 = Ω 0 e 2iγ 0 , we can then write the expansion as 7 In general Ω 0 is expected to depart from unity.The corresponding polarized P ν (central panel) and polarization angle γ ν (lower panel).We choose logarithmic representations to emphasize the focus on the ν ∼ ν 0 point.
In the perturbative regime, we thus find a polarization SED that is again a power law with a spectral index and frequencydependent polarization angle given by 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, β, by a ∆ β.At higher order, the complexvalued 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 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 β fixed.As such, one breaks unwanted degeneracies and recovers a physically relevant value for the spectral parameters with a minimal dispersion.

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: x e x (e x − 1) 2 , with x = hν/k T and the usual natural constants.We also introduced the two temperature perturbations Θ Q = ∆T Q / T and Θ U = ∆T U / T , which are respectively defined for a coordinate system that is rotated by 45 • .At fixed Θ Q and Θ U , this means that (at lowest order in the temperature fluctuations) no spectral mixing happens, and hence γ remains frequency-independent.However, if we include terms at second order in Θ, one finds an additional frequency dependence that is characterized by a y-type distortion (e.g., see Appendix A of Chluba et al. 2015): Here, we introduced the total intensity temperature perturbation Θ I = ∆T I / T , which generally includes both polarized and unpolarized contributions.Depending on the ratio of Θ Q to Θ 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.

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 PGB ν ≡ B ν (T ), such that the single polarization state can be characterized by P GB ν ≈ A e 2iγ B ν (T ).We can then consider general GB superpositions.Following Chluba et al. (2017), we shall use β 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): x coth(x/2) (39b) x 2 cosh(x) + 2 cosh(x) − 1 (39c) with the frequency variable x = hν k βGB ≡ hν k T .The final GB moment expansion then takes the form: The functions Y GB k (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 To catch the domain on which β GB has a maximal impact, we choose ν max = 5000 GHz.We note the "loop" trajectory of P ν 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 ν , 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 | ≫ |W β GB 1 |, we can again use the split ∆ βGB = W β GB 1 /W 0 = a ∆ βGB + ib ∆ βGB into real and imaginary parts.With this, we can then write with T = 1/( βGB + ∆ βGB ).The leading order SED term, PGB ν ( T ), then depends on the function 1 with x R = hν ( βGB + a ∆ βGB )/k and x I = hν b ∆ βGB /k.One can see that polarized mixing leads to an imaginary photon chemical potential, µ = ix I .This causes a frequency-dependent rotation of the polarization plane and also modifications to the SED.At high frequencies, one finds ∆γ ν ≈ x I /2, while at low frequencies, one has the constant ∆γ ν ≈ 1 2 tan −1 (b ∆ βGB /a ∆ βGB ).

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 PmBB , such that the single polarization state can be characterized by Using the results for the power law and gray-body spectra of the previous sections, we then have 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 normalized9 modified blackbodies of parameters 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 graybody factor at high frequencies, we choose ν max = 800 GHz.Accordingly to the Planck high frequency bands, we choose ν 0 = 353 GHz (Planck Collaboration I 2016).One can see that the different power laws induce strong distortions at low frequencies ≤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 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 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 ν.

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.

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 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).

In intensity
While in the above, we described the signal along a given line of sight n, we now turn ourselves to averages between different lines-of-sights across sky patches.The average intensity, ⟨I ν ⟩, which generally depends on the line-of-sight moments and pivots of the SED expansion in the direction 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 ℓ m ( n).We shall denote the spherical harmonic coefficients of a quantity X as The spherical harmonic coefficients of the average intensity SED can then be expressed as: where P ℓm ( p) is the harmonic coefficient of the parameter distribution, P( p, n).Just as the single line-of-sight average, this superposition will introduce some additional mixing between intensities with different spectral parameters p( n) and introduces extra distortions.The spherical harmonic functions introduce a reweighting of the parameter distribution, which essentially A5, page 11 of 16 A&A 669, A5 (2023) implies that all moments can be directly computed using the harmonic coefficients of the parameter distribution, (P) ℓm ( p).
Expanding Îν (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): where δ i j is the Kronecker-δ and we introduced the (complexvalued) moment multipoles In Eq. ( 47), we used the average of all the parameters across the sky to define the average SED amplitude and pivot10 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.

In Polarization
The averaged polarized signal ⟨P ν (p( n))⟩ = ⟨P ν ( n)⟩ 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 ±2 Y ℓm (Newman & Penrose 1966;Goldberg et al. 1967).One can evaluate the spherical-harmonic coefficients of a general frequency dependent spinor field ±2 X ν for ℓ ≥ 2 as 11 For convenience, we then also define the harmonic expansion of the line-of-sight average of the polarization field as where we define the harmonic coefficient of the parameter distribution function, P( p, γ, n), as ±2 P ℓm ( p, γ).When applied to the spin-moment expansion of the polarization field, we then obtain Introducing the general instrumental transfer function W(ν, γ, 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.The total polarized signal in a frequency band of width ∆ν at average frequency ν c and within the spatial support of the beam, Ω, centered in an average line of sight nc then becomes: Here, ⟨P ν ⟩ W is now a function of ν c and the corresponding moments and pivots in the average direction nc .
Due to the ν integration in the above expression, it is not possible to simply factor the spectral shapes in the Taylor expansion of Pν 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: with the new averages defined as The pivot is chosen to be the same as in intensity: 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 B as W(ν, n) = B(γ, n)F(ν) 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 B. The result remains nontrivial since the spectral dependence of the moments need to be computed through the potentially complicated integral of F(ν)∂ p j ...p k Pν over ∆ν.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.

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, β.Considering only Q ν , from the discussion in Sect.3.3 we then have A similar expression follows for U ν , with the relevant parameters A × , β × , A ⊗ , β ⊗ .Let us again expand the fundamental SEDs around some average index β inside each voxel12 .This then yields 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.

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 where f ν = ν −2 and Γ(s) ∝ s 0 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: 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 with ĀFR = ⟨A⟩ 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 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 P( p, γ, Γ, n) ≈ P( p, γ, n) P(Γ, 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: (67) This expression shows that no new spectral mixing occurs in this case, as all moments are multiplied by the same frequencydependent factor.The average Faraday rotation coefficient is simply given by Using this as a pivot, we can then write series illustrating how a complicated frequency structure can be created from higher order moments of Γ − Γ.Given that these depend on powers of f ν , one can in principle determine these moments observationally.Finally, if the simple factorization of the parameter distribution function is not possible, one can write This generally introduces a complicated voxel SED reweighting by factors of ( f ν ) k and hence new SED shapes that in principle can all be separately accounted for.As is understood, the number of variables quickly becomes unmanagable unless highly perturbative situations are encountered.A more detailed discussion is, however, beyond the scope of this work.

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 powerspectrum level, especially for E-and B-mode applications.One could also think of links with cosmic birefringence, which A5, page 14 of 16 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.

Fig. 2 .
Fig.2.Illustration of the spinor P ν 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 ν (central panel) and polarization angle γ ν (lower panel).The exact result with associated (invisible) error bars (black) is compared to the best-fit moment representation at various orders.

Fig. 3 .
Fig.3.Illustration of the spinor P ν in the complex plane (Q, U) for a sum of two almost anti aligned power laws, with a close-up view around ν = ν 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 P ν (central panel) and polarization angle γ ν (lower panel).We choose logarithmic representations to emphasize the focus on the ν ∼ ν 0 point.
et al.: High precision modeling of polarized signals

Fig. 4 .
Fig. 4. Illustration of P ν 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 ν (central panel) and polarization angle γ ν (lower panel).

Fig. 5 .
Fig. 5. Illustration of the spinor P ν 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 ν (central panel) and polarization angle γ ν (lower panel).