© The Authors 2026

1 Introduction

Observations from the visible to the submillimeter wavelengths (submm) provide extensive evidence on the variations in dust composition and optical properties through the Galaxy. The importance of dust processing in the interstellar medium (ISM) is demonstrated first by the large variations seen in the extinction curve and the depletions of refractory elements (Draine 2003; Jenkins 2009; Schlafly et al. 2016). Additional evidence comes from variations in the spectral energy distribution (SED) and the far-infrared and submm dust opacity of molecular clouds and the diffuse ISM, observed by the Herschel (Ysard et al. 2013) and Planck space missions (Planck Collaboration 2013 XI 2014; Planck Collaboration Int. XVII 2014). These variations of the dust SED include, but are not fully explained by, the dependence of dust temperature on the local intensity of the interstellar radiation field (Fanciullo et al. 2015).

The observational signatures of dust evolution in the dust SED in polarization are harder to identify due to the reduced amplitude of the signal compared to the total intensity (Pelgrims et al. 2021; Ritacco et al. 2023). The interpretation of polarization data considers the combined variations in the structure of the Galactic magnetic field and the dust SED, both in the plane of the sky and along the line of sight in the three-dimensional Galaxy (Tassis & Pavlidou 2015; Planck Collaboration Int. L 2017; McBride et al. 2023; Vacher et al. 2025). This is the specific problem that our work addresses. We present a modeling framework that incorporates variations in dust properties within a turbulent magnetic field, designed to quantify the variance in dust spectral properties from submm polarization data.

In submm polarization maps, the diffuse polarized emission from dust grains that are locally aligned with magnetic-field orientation is summed over the light cone, a small region of the sky observed by the instrumental beam containing both plane-of-the-sky and line-of-sight variations. Variations in the magnetic field orientation within the light cone lead to depolarization (Planck Collaboration Int. XLIV 2016; Clark & Hensley 2019), and their combined variations with changes in dust properties lead to the frequency-dependence of the polarization angle (Tassis & Pavlidou 2015; Planck Collaboration Int. L 2017; Vacher et al. 2023b). For studies of the cosmic microwave background (CMB) polarization, the variation over the sky of dust spectral parameters, and the non-linear distortion they induce in the SEDs, together with their impact on the frequency dependence of polarization angles, makes the separation of the dust and CMB signals a difficult task, especially given the high sensitivity of next-generation experiments (Remazeilles et al. 2016; LiteBIRD Collaboration 2022; Wolz et al. 2024). The modeling of spatial variations in the dust polarization SED is therefore of prime importance for the search for primordial B-modes from inflation and for the putative EB correlation from cosmic birefringence (Diego-Palazuelos et al. 2022, 2023; Vacher et al. 2023a; Jost et al. 2023; Sullivan et al. 2025).

Planck observations show that the mean SED of dust polarized emission in the high-latitude diffuse ISM is remarkably well fitted by a modified black body (MBB) law from 353 GHz to microwave frequencies (Planck Collaboration 2018 XI 2020). This observational result leads us to model spatial variations of the dust polarization SED using the moment expansion method around an MBB, as introduced by Chluba et al. (2017). Previous studies carried out in preparation for future experiments, such as the Simons Observatory and LiteBIRD, show that when spatial variations are present in the foreground SED, the moment-expansion formalism provides a powerful tool to recover the CMB signal without bias, and with minimal parameter addition (Azzoni et al. 2021; Remazeilles et al. 2021; Désert 2022; Vacher et al. 2022; Sponseller & Kogut 2022; Azzoni et al. 2023; Wolz et al. 2024; Carones & Remazeilles 2024). The moment expansion formalism can be extended to polarized dust emission using complex spin-2 fields, referred to as"spin-moment" coefficients (Vacher et al. 2023b,a). In this work, we use this formalism to introduce a modeling framework that, within some simplifying assumptions, accounts for the combined variations of dust properties and magnetic field structure. Applying it to Planck data, we demonstrate the relevance of this framework for analyzing dust polarization microwave observations for foreground modeling and component separation.

Our paper is structured as follows. In Sect. 2, we introduce the dust model and data modeling framework based on complex residual maps. We enumerate the predictions of our model in Sect. 3. In Sect. 4, building on the work of Ritacco et al. (2023), we test the validity of our simplifying assumptions in relation to sky observations by analyzing the Planck SRoll2 maps (Delouis et al. 2019). In Sect. 5, we interpret our results in terms of fluctuations in dust spectral properties and compare them with those obtained using the PR4 version of the Planck data (Planck Collaboration Int. LVII 2020). In Section 6, we discuss our results, which we then summarize in Sect. 7.

2 Model definition

Our dust model and methodology, based on complex numbers, are designed to quantify the variance in the emission properties of grains aligned with magnetic field lines¹ which emit polarized radiation in the far-infrared and submm. In Sect. 2.1, we set out the assumptions of our dust signal model and in Sect. 2.2 we show how these hypotheses enable accurate modeling of the polarized signal using complex numbers and the moment expansion. In Sect. 2.5, we define the so-called polarization residual maps, from which we derive various covariances in Sect. 2.6. These covariances are then used to build the predictions of our model in Sect. 3.

2.1 Dust model hypothesis

Our dust model relies on simplifying hypotheses of two kinds: astrophysical (HA) and methodological (HM), listed below.

1. HA1: the local dust polarized SED emitted at any given position in the Galaxy is a single MBB with local temperature T and local spectral index β.

2. HA2: the variations of the dust polarization spectral parameters T or β, from one point of the Galaxy to another, are not correlated with variations in the orientation of the magnetic field.

3. HM3: only two extreme cases are explored: the fluctuations of dust properties are either pure β fluctuations over the whole sky or pure T fluctuations.

4. HA4: the grain alignment efficiency is uniform over the sky; only the structure of the magnetic field controls the dust polarization fraction. This hypothesis was shown by Planck Collaboration 2018 XII (2020) to be justified at the scale of a degree or more.

5. HM5: following Planck Collaboration Int. XLVIII (2016), we modeled the structure of magnetic fields as a turbulent field close to equipartition between its ordered and random components.

6. HM6: our data analysis methodology accounts for contamination of the dust polarization signal in Planck data by CMB, noise, and systematics, but it neglects possible contribution from polarized CO emission or leakage of CO emission in the 217 and 100 GHz channels, as well as the presence of residual synchrotron emission after template removal.

While these hypotheses are quite restrictive, our model is able to generate a frequency-dependence of polarization angles. The simplicity of the model and the clarity of its building hypotheses allow one to infer its predictions unambiguously and to identify its limits precisely. To assist the reader with the numerous notations, Table A.1 summarizes all the quantities used in our model.

2.2 Spin-moment expansion of the complex polarized SED

Following the work of Ritacco et al. (2023), we developed an analytical framework to model the observed variations with the frequency of polarization intensity and angle. For this purpose, we did not use the E and B decomposition (Zaldarriaga & Seljak 1997). Instead, we quantified the total variance of the Stokes parameters. For a given frequency ν, we define the dust complex polarized SED P_ν from the Stokes parameters Q_ν and U_ν as follows (complex quantities are written in bold font):

$${\mathbf{P}}_{\nu }\equiv Q_{\nu }+\dot{\mathbb{i}}{U}_{\nu }=P_{\nu }{\mathrm{e}}^{\dot{\mathbb{i}}2{\psi }_{\nu }},$$(1)

where $\dot{\mathbb{i}}$ is the imaginary unit $({\dot{\mathbb{i}}}^2=-1),P_{\nu }$ is the polarized intensity ($P_{\nu }=\sqrt{Q_{\nu }^2+U_{\nu }^2}=\sqrt{{\mathbf{P}}_{\nu }{\mathbf{P}}_{\nu }^{\star }}$, with ^⋆ denoting complex conjugation), and ψ_ν = 0.5 arctan(U_ν, Q_ν) is the frequency-dependent polarization angle.

Following Vacher et al. (2023b), we note that the use of complex spin-moments is well suited to studying fluctuations in the dust complex polarized SED. The SED P_ν observed over a region of the sky is given by the integral over all the local elementary SEDs dP_ν present within the observed light cone ℒ:

$${\mathbf{P}}_{\nu }={\int }_{\mathcal{L}}\mathrm{d}{\mathbf{P}}_{\nu }.$$(2)

Variations in phases in this integral can partially or fully cancel the polarization signal, a phenomenon known as depolarization (Planck Collaboration Int. XIX 2015). The so-called light cone, ℒ, is a region of the sky over which all local signals dP_ν – contained both in the plane of the sky and along the line of sight – are averaged. Depending on the context of the analysis, ℒ can refer to a map pixel, the instrumental beam, or even the entire sky, possibly masked. In the present work, we considered ℒ to be the instrumental beam surrounding each pixel of the sky.

We used an MBB model of dust emission (hypothesis HA1, Sect. 2.1) characterized by a local spectral index β and a local temperature T to express the local dust emissivity ε_ν(T, β) at frequency ν,

$${\epsilon }_{\nu }(\beta ,T)=B_{\nu }(T){\left(\frac{\nu }{{\nu }_0}\right)}^{\beta },$$(3)

where B_ν(T) is Planck’s law at frequency ν and temperature T, and ν₀ is a reference frequency. Assuming that the local complex polarized intensity, dP_ν0, at frequency ν₀ is well-defined in amplitude and phase everywhere in the light cone, the frequency-dependence of the local complex polarized SED, dP_ν, is

$$\mathrm{d}{\mathbf{P}}_{\nu }=\mathrm{d}{\mathbf{P}}_{{\nu }_0}\frac{{\epsilon }_{\nu }(\beta ,T)}{{\epsilon }_{{\nu }_0}(\beta ,T)}=\mathrm{d}{\mathbf{P}}_{{\nu }_0}\frac{B_{\nu }(T)}{B_{{\nu }_0}(T)}{\left(\frac{\nu }{{\nu }_0}\right)}^{\beta }.$$(4)

Because the emissivity ε_ν is a real quantity, dP_ν has the same polarization angle at all frequencies.

To estimate the impact of fluctuations in the spectral parameter s ∈ {T, β} on the observed complex polarized intensity P_ν (under our hypothesis HM3, which assumes that either T or β fluctuates over the whole sky but not both), we performed a Taylor expansion of the power-law term and the blackbody in dP_ν in Eq. (4) around an arbitrary pivot spectral index $\overline{\beta }$ and pivot temperature $\overline{T}$. We find that the frequency dependence of the complex polarized intensity can be expressed as (Vacher et al. 2023b)

$${\mathbf{P}}_{\nu }={\mathbf{P}}_{{\nu }_0}\frac{{\overline{\epsilon }}_{\nu }}{{\overline{\epsilon }}_{{\nu }_0}}(1+{\sum }_{n=1}^{\mathrm{\infty }}{a}_{\nu ,n}^s{\mathcal{W}}_n^s),$$(5)

where

$${\overline{\epsilon }}_{\nu }\equiv {\epsilon }_{\nu }(\overline{\beta },\overline{T})$$(6)

is the pivot SED, the function $a_{\nu ,n}^s$ of ν denotes the nth order derivative of the dust emissivity with respect to the spectral parameter s,

$$a_{\nu ,n}^s\equiv \frac{1}{n!}\frac{{\overline{\epsilon }}_{{\nu }_0}}{{\overline{\epsilon }}_{\nu }}\frac{{\mathrm{\partial }}^n\left({\epsilon }_{\nu }/{\epsilon }_{{\nu }_0}\right)}{\mathrm{\partial }{s}^n}(\overline{\beta },\overline{T}),$$(7)

and ${\mathcal{W}}_n^s$ is the nth order spin-moment of s (Vacher et al. 2023b), given² by

$${\mathcal{W}}_n^s\equiv \frac{1}{{\mathbf{P}}_{{\nu }_0}}{\int }_{\mathcal{L}}(s-\overline{s}{)}^n\mathrm{d}{\mathbf{P}}_{{\nu }_0}.$$(8)

The spin-moments ${\mathcal{W}}_n^s$ encode the departure of the spectral dependence of P_ν from the pivot SED ${\overline{\epsilon }}_{\nu }$. The pivot spectral parameters $\overline{\beta }$ and $\overline{T}$ are constant over ℒ. For simplicity, we further considered a single pair of pivots for all light cones and hence the whole sky. In this analysis, we used $\overline{T}\equiv 20\mathrm{K}$. The specific value for $\overline{\beta }$ is not critical, as it does not appear in the data analysis.

The frequency dependence of $a_{\nu ,n}^T$ and $a_{\nu ,n}^{\beta }$ is presented in Fig. 1 for n = 1, 2, and 3. Their analytical expression for n = 1 is as follows:

$$d_{\nu ,1}^{\beta }\equiv \ln \left(\nu /{\nu }_0\right),$$(9)

$$a_{\nu ,1}^T\equiv \frac{h}{k{\overline{T}}^2}[\frac{\nu {\mathrm{e}}^{h\nu /k\overline{T}}}{{\mathrm{e}}^{h\nu /k\overline{T}}-1}-\frac{{\nu }_0{\mathrm{e}}^{h{\nu }_0/k\overline{T}}}{{\mathrm{e}}^{h{\nu }_0/k\overline{T}}-1}].$$(10)

The values of the $a_{\nu ,n}^s$ functions – hereafter referred to as the spectral gradients – are listed explicitly up to n = 5 for both β and T in Chluba et al. (2017).

The polarization-weighted moment ${\mathcal{W}}_n^s$ is distinct from the intensity-weighted moments ${\omega }_n^s$ (Chluba et al. 2017):

$${\omega }_n^s\equiv \frac{1}{I}{\int }_{\mathcal{L}}(s-\overline{s}{)}^n\mathrm{d}I,$$(11)

with s ∈ {T, β} as before. Like ${\mathcal{W}}_n^s,{\omega }_n^s$ characterizes the statistical properties of the spectral parameters s of aligned grains. However, ${\omega }_n^s$ is independent of the structure of the magnetic field, whereas ${\mathcal{W}}_n^s$ is not. The difference between these two moments,

$${\mathrm{\Delta }}_n^s\equiv {\mathcal{W}}_n^s-{\omega }_n^s,$$(12)

is a complex quantity of interest. In Appendix B, we demonstrate that, under our hypothesis HA2, which assumes no correlation between the fluctuations of dust spectral properties and variations in the magnetic field direction, ${\mathrm{\Delta }}_n^s$ has a statistically zero mean: $⟨{\mathbf{\Delta }}_n^s⟩=0$.

Fig. 1 Values of $a_{\nu }^{\beta }$ and $a_{\nu }^T$ for orders n = 1 (blue), n = 2 (orange), and n = 3 (red), shown for ν0 = 353 GHz and pivot temperature $\overline{T}=20\mathrm{K}$. The points represent the four Planck HFI channels.

2.3 Spectral rotation of polarization angles

To gain further insight into how the spin-moments, ${\mathcal{W}}_n^s$, are related to the spectral dependence of P_ν, we first consider the specific case of small distortions of the polarized SED around the pivot SED ($\left|{\sum }_{n=1}^{\mathrm{\infty }}{a}_{\nu ,n}^s{\mathcal{W}}_n^s\right|\ll 1$). Taking the logarithm of Eq. (5) and using the approximation ln(1 + x) ~ x for x ≪ 1, we obtain

$$\ln \left(\frac{{\mathbf{P}}_{\nu }}{{\overline{\epsilon }}_{\nu }}\right)-\ln \left(\frac{{\mathbf{P}}_{{\nu }_0}}{{\overline{\epsilon }}_{{\nu }_0}}\right)\simeq {\sum }_{n=1}^{\mathrm{\infty }}{a}_{\nu ,n}^s{\mathcal{W}}_n^s,$$(13)

from which we take the real and imaginary parts and, using Eq. (1),

$$\ln \left(\frac{P_{\nu }}{{\overline{\epsilon }}_{\nu }}\right)-\ln \left(\frac{P_{{\nu }_0}}{{\overline{\epsilon }}_{{\nu }_0}}\right)\simeq {\sum }_{n=1}^{\mathrm{\infty }}{a}_{\nu ,n}^s\mathrm{Re}{\mathcal{W}}_n^s,$$(14)

$$2({\psi }_{\nu }-{\psi }_{{\nu }_0})\simeq \sum _{n=1}^{\mathrm{\infty }}{a}_{\nu ,n}^s\mathrm{Im}{\mathcal{W}}_n^s.$$(15)

Equation (15) describes the spectral rotation of the polarization angle with frequency, caused by fluctuations of the dust spectral parameters within the light cone. Its frequency dependence is governed by the $a_{\nu ,n}^s$ parameters and its amplitude depends on the imaginary parts ${\mathrm{I}\mathrm{m}\mathcal{W}}_n^s$ of the spin-moments. Figure 2 provides a numerical illustration of such a rotation in the log-complex plane for the sum of seven MBBs with Gaussian-distributed spectral parameters in a turbulent magnetic-field model inspired by Planck Collaboration Int. XLVIII (2016) (see also our description in Appendix C). We consider the two cases of our hypothesis HM3: pure fluctuations in T or pure fluctuations in β. Under the assumption of small variations in T or β, the complex SED plotted in the log-complex plane (ln (P_ν/ϵ_ν), 2ψ_ν) is close to a straight line for T fluctuations and can be non-linear for β fluctuations when the magnetic field direction is close to the line of sight (as shown in Fig. 2 for illustration). This shows that, in a given light cone, the rotation of the polarization angle with frequency is strongly correlated with the distortion of the polarized intensity relative to the pivot SED.

2.4 Mean complex normalized SED

From here on, we simplify our notation for channel i of frequency ν_i and use P_i to denote the map of complex polarized intensity, ${\overline{\epsilon }}_i$ for the pivot emissivity (constant over the sky), and $a_{i,n}^s$ for the spectral gradient.

Following Ritacco et al. (2023), we cross-correlated the polarization maps to determine the mean complex SED ${\overline{\mathbf{r}}}_i$ of dust polarization, averaged over the sky and normalized to our reference channel, given by

$${\overline{\mathbf{r}}}_i\equiv \frac{⟨{\mathbf{P}}_i{\mathbf{P}}_0^{\star }⟩}{⟨{\mathbf{P}}_0{\mathbf{P}}_0^{\star }⟩}=\frac{{\overline{\epsilon }}_i}{{\overline{\epsilon }}_0}(1+{\sum }_{n=1}^{\mathrm{\infty }}{a}_{i,n}^s{\overline{\mathcal{W}}}_n^s),$$(16)

where ${\overline{\mathcal{W}}}_n^s$ is the mean value of ${\mathcal{W}}_n^s$ over the sky, weighted by $P_0^2={\mathbf{P}}_0{\mathbf{P}}_0^{\star }$ as

$${\overline{\mathcal{W}}}_n^s\equiv ⟨P_0^2{\mathcal{W}}_n^s⟩/⟨P_0^2⟩.$$(17)

Fig. 2 Rotation of the SED in the complex plane for a sum of seven MBB SEDs with distinct spectral indices drawn from a Gaussian distribution with mean 1.5 and standard deviation 0.1 (red), or with temperatures drawn from a Gaussian distribution with mean 20 K and standard deviation 5 K (blue). The 3D orientation of the magnetic field is a random realization of the turbulent magnetic field inspired by Planck Collaboration Int. XLVIII (2016). The dashed lines indicate frequencies from 40 to 400 GHz, and the four colored points represent the Planck HFI bands.

2.5 Residual maps for complex polarization

In this section, we use the pivot SED to construct complex residual maps at all frequencies and relate them to the ${\mathcal{W}}_n^s$ moments. We define the residual map R_i for the complex polarization P_i:

$${\mathbf{R}}_i\equiv P_0(\frac{1}{{\overline{\mathbf{r}}}_i}\frac{{\mathbf{P}}_i}{{\mathbf{P}}_0}-1).$$(18)

This complex map characterizes the residuals at frequency ν_i after removal of the components correlated with the polarization map at the reference frequency ν₀, and subsequently scaled to the reference frequency ν₀. Here, R_i quantifies the complexity of the dust signal by measuring its deviations from the mean complex SED ${\overline{\mathbf{r}}}_i$.

Using Eqs. (5) and (16), we obtain

$${\mathbf{R}}_i=P_0(\frac{1+\sum _{n=1}^{\mathrm{\infty }}{a}_{i,n}^s{\mathcal{W}}_n^s}{1+\sum _{n=1}^{\mathrm{\infty }}{a}_{i,n}^s{\overline{\mathcal{W}}}_n^s}-1).$$(19)

By choosing the pivot SED close to the mean SED, we ensure that $\left|{\sum }_{n=1}^{\mathrm{\infty }}{a}_{i,n}^s{\overline{\mathcal{W}}}_n^s\right|=|{\overline{\mathbf{r}}}_i\frac{{\overline{\epsilon }}_0}{{\overline{\epsilon }}_i}-1|\ll 1$ (this hypothesis is verified in the data analysis in Sect. 4), so that, by expanding the previous equation to first order, we obtain

$${\mathbf{R}}_i\simeq P_0{\sum }_{n=1}^{\mathrm{\infty }}{a}_{i,n}^s({\mathcal{W}}_n^s-{\overline{\mathcal{W}}}_n^s).$$(20)

We decompose R_i into its real and imaginary components:

$$R_i^P\equiv \mathrm{Re}{\mathbf{R}}_i=P_0{\sum }_{n=1}^{\mathrm{\infty }}{a}_{i,n}\mathrm{Re}({\mathcal{W}}_n^s-{\overline{\mathcal{W}}}_n^s),$$(21)

$$R_i^{\psi }\equiv \mathrm{Im}{\mathbf{R}}_i=P_0{\sum }_{n=1}^{\mathrm{\infty }}{a}_{i,n}\mathrm{Im}({\mathcal{W}}_n^s-{\overline{\mathcal{W}}}_n^s).$$(22)

For $R_i^{\psi }$ to be non-zero in a given light cone ℒ, both the magnetic field and the spectral parameters should vary inside ℒ, whereas a variation in spectral parameters alone suffices for $R_i^P$ to be non-zero. We use the labels P and ψ for these map residuals because, as demonstrated in Sect. 2.3, in the limit of small fluctuations, the real and imaginary parts of the complex polarized SED correspond directly to the polarized intensity P and polarization angle ψ.

Similarly, we define the map residual $R_i^{\epsilon }$ of the dust polarized emissivity that corresponds to the moments ${\omega }_n^s$ (Eq. (11)) as

$$R_i^{\epsilon }\equiv P_0{\sum }_{n=1}^{\mathrm{\infty }}{a}_{i,n}^s({\omega }_n^s-{\overline{\omega }}_n^s).$$(23)

Here, ${\overline{\omega }}_n^s\equiv ⟨P_0^2{\omega }_n^s⟩/⟨P_0^2⟩$ denotes the sky-averaged value of ${\omega }_n^s$. The map residual $R_i^{\epsilon }$ represents what would be observed if the polarization angles were constant across each light cone. It maps the variations of the spectral properties of aligned grains, independently of the magnetic field structure, while $R_i^P$ mixes both.

2.6 Covariance matrices of residual maps

Consider a complex map X (resp. Y) with sky-averaged values $\overline{\mathbf{X}}$ and $\overline{\mathbf{Y}}$, defined as

$$\overline{\mathbf{X}}=⟨P_0^2\mathbf{X}⟩/⟨P_0^2⟩.$$(24)

We use the notation Cov(X, Y) for the covariance of X and Y weighted by $P_0^2$ and computed over the sky:

$$\mathrm{Cov}(\mathbf{X},\mathbf{Y})\equiv ⟨P_0^2(\mathbf{X}-\overline{\mathbf{X}})(\mathbf{Y}-\overline{\mathbf{Y}}{)}^{\star }⟩,$$(25)

where the quantity inside the brackets must be averaged over the partial sky after application of the mask.

For all frequency pairs (i, j), we define covariances for the polarization map residuals R^P and R^ψ:

$$c_{ij}^P\equiv ⟨R_i^P{R}_j^P⟩={\sum }_{n,m=1}^{\mathrm{\infty }}{a}_{i,n}^s{a}_{j,m}^s\mathrm{Cov}(\mathrm{Re}{\mathcal{W}}_n^s,\mathrm{Re}{\mathcal{W}}_m^s),$$(26)

$$c_{ij}^{\psi }\equiv ⟨R_i^{\psi }{R}_j^{\psi }⟩={\sum }_{n,m=1}^{\mathrm{\infty }}{a}_{i,n}^s{a}_{j,m}^s\mathrm{Cov}(\mathrm{Im}{\mathcal{W}}_n^s,\mathrm{Im}{\mathcal{W}}_m^s).$$(27)

While the residual maps R represent values for each light cone of the map, the covariances c are averaged quantities that sum up the statistics of all light cones that make up the masked sky. Because we are working at the map scale and not in Fourier space, this mask does not need to be continuous.

We also define the mixed-covariance of R^ψ with R^P:

$$c_{ij}^{\psi P}\equiv ⟨R_i^{\psi }{R}_j^P⟩={\sum }_{n,m=1}^{\mathrm{\infty }}{a}_{i,n}^s{a}_{j,m}^s\mathrm{Cov}(\mathrm{Im}{\mathcal{W}}_n^s,\mathrm{Re}{\mathcal{W}}_m^s),$$(28)

and its companion $c_{ij}^{P\psi }\equiv c_{ji}^{\psi P}$.

Finally, we define the covariance matrix c^ε, which characterizes the fluctuations of the emissivity of aligned grains as

$$c_{ij}^{\epsilon }\equiv ⟨R_i^{\epsilon }{R}_j^{\epsilon }⟩={\sum }_{n,m=1}^{\mathrm{\infty }}{a}_{i,n}^s{a}_{j,m}^s\mathrm{Cov}({\omega }_n^s,{\omega }_m^s).$$(29)

This covariance is the only one not contaminated by variations in the magnetic field structure and so can be directly compared safely to dust models or observations in total intensity.

All these quantities are N × N matrices of real numbers, with N the number of channels. In contrast, the covariance c^P of the complex map residual R,

$${\mathbf{c}}_{ij}^{\mathbf{P}}\equiv ⟨{\mathbf{R}}_i{\mathbf{R}}_j^{\star }⟩={\sum }_{n,m=1}^{\mathrm{\infty }}{a}_{i,n}^s{a}_{j,m}^s\mathrm{Cov}({\mathcal{W}}_n^s,{\mathcal{W}}_m^s),$$(30)

is an N × N matrix of complex numbers, related to other real matrices through

$${\mathbf{c}}_{ij}^{\mathbf{P}}=⟨R_i^P{R}_j^P+R_i^{\psi }{R}_j^{\psi }⟩+\dot{\mathbb{i}}⟨R_i^{\psi }{R}_j^P-R_i^P{R}_j^{\psi }⟩,$$(31)

$$=(c_{ij}^P+c_{ij}^{\psi })+\dot{\mathbb{i}}(c_{ij}^{\psi P}-c_{ij}^{P\psi }).$$(32)

We therefore have the trivial equalities, already expected from the definition of R^P and R^ψ (Eq. (21) and (22)):

$$\mathrm{Re}{\mathbf{c}}^{\mathbf{P}}=c^P+c^{\psi },$$(33)

$$\mathrm{Im}{\mathbf{c}}^{\mathbf{P}}=c^{\psi P}-c^{P\psi }.$$(34)

By symmetry, $\mathrm{Im}{\mathbf{c}}_{ij}^{\mathbf{P}}=0$ if i = j. For i ≠ j, we have

$$\mathrm{Im}{\mathbf{c}}_{ij}^{\mathbf{P}}={\sum }_{m>n\ge 1}(a_{i,n}^s{a}_{j,m}^s-a_{i,m}^s{a}_{j,n}^s)\mathrm{Im}\mathrm{Cov}({\mathcal{W}}_n^s,{\mathcal{W}}_m^s).$$(35)

3 Predictions of our model

Here and throughout this work, predictions of our model are highlighted by a boxed frame.

3.1 Expected level of correlation between residual maps

Our dust model assumes that only a single spectral parameter, β or T, fluctuates in the light cone, but not both simultaneously (hypothesis HM3, Sect. 2.1). Fluctuations in T or β propagate into the residual maps. We next examine how well the residual maps correlate across frequencies.

Limiting the moment expansion of Eq. (20) to first order gives ${\mathbf{R}}_i\simeq a_{i,1}^s{P}_0({\mathcal{W}}_{\mathbf{1}}^{\mathit{s}}-{\overline{\mathcal{W}}}_1^s)$ and ${\mathbf{R}}_j\simeq a_{j,1}^s{P}_0({\mathcal{W}}_{\mathbf{1}}^{\mathit{s}}-{\overline{\mathcal{W}}}_1^s)$. Residuals maps become simply proportional to each other. The correlation of R_i with R_j is perfect. This conclusion also holds for R^P and for R^ψ. Introducing higher-order n weakens this correlation by adding the common signals ${\mathcal{W}}_n^s$ in channels i and j, but are weighted with frequency-dependent coefficients a_i,n and a_j,n, which differ from a_i,1 and a_j,1. From the values of a_i,n (also displayed in Fig. 1), we computed the Pearson correlation coefficients ρ(a_i,n, a_i,m) across frequencies ν_i between the spectral parameters taken at two distinct orders n and m. We find $\rho (a_{i,1}^{\beta },a_{i,2}^{\beta })=-0.958$ and $\rho (a_{i,1}^T,a_{i,2}^T)=-0.999$ between orders n = 1 and m= 2, and $\rho (a_{i,1}^{\beta },a_{i,3}^{\beta })=0.905$ and $\rho (a_{i,1}^T,a_{i,3}^T)=0.998$ between orders n = 1 and m = 3. Despite changes in amplitude or sign from one order to another, the spectral gradients a_i,n remain well correlated across the first, second, and third orders. This correlation is particularly strong for temperature fluctuations. Considering that the first order, which produces a perfect frequency-correlation, dominates the signal, we expect the sum ${\sum }_{n=1}^{\mathrm{\infty }}{a}_{i,n}^s({\mathcal{W}}_n^s-{\overline{\mathcal{W}}}_n^s)$ to remain very well correlated through frequencies ν_i. Within the framework of our hypothesis HM3, our model, even including orders higher than one in the moment expansion, predicts a perfect correlation of residual maps across frequencies if only the dust temperature varies, and a strong but not perfect correlation if only the dust spectral index varies. We therefore obtain the first prediction P1 of our model:

$$\boxed{{\displaystyle {\rho }_{ij}^{\mathbf{P}}({\mathbf{R}}_i,{\mathbf{R}}_j)\simeq {\rho }_{ij}^P(R_i^P,R_j^P)\simeq {\rho }_{ij}^{\psi }(R_i^{\psi },R_j^{\psi })\simeq 1.}}$$(P1)

A word of caution is necessary. We are stating that the residuals maps at different frequencies are correlated with each other, not the maps of the dust signal themselves. As such, dust complexity creates frequency decorrelation in the sense of Planck Collaboration Int. XXX (2016), meaning that the maps at different frequencies are no longer proportional. However, the complex distortions remain strongly correlated. This correlation can be understood geometrically as resulting from the rigid rotation of the complex polarized intensity in the complex plane (see Sect. 2.3 and Fig. 2), which arises from beam and line-of-sight integration effects.

3.2 Properties of covariance matrices of residual maps

We now examine the properties of the covariances of residual maps across frequencies to provide clear predictions that can readily be compared with Planck polarization data.

Using our assumption of no correlation between the magnetic field structure and the fluctuations of dust properties (hypothesis HA2), Eq. (25) yields $\mathrm{Cov}({\omega }_n^s,{\omega }_m^s)\propto ⟨P_0^2⟩$. From Eq. (29) we therefore obtain the following scaling relation:

$$\boxed{{\displaystyle c^{\epsilon }\propto ⟨P_0^2⟩.}}$$(P2)

This prediction (P2) is natural: the variance of dust polarized emissivity scales with the square of the mean polarized intensity.

To explore the properties of the variance of polarization angles c^ψ and intensity c^P, we used the moment difference ${\mathrm{\Delta }}_n^s={\mathcal{W}}_n^s-{\omega }_n^s$ (Eq. (12)). We defined its sky-averaged value as $\overline{{\mathbf{\Delta }}_n^s}\equiv ⟨P_0^2{\mathbf{\Delta }}_n^s⟩/⟨P_0^2⟩$. As ${\mathbf{\Delta }}_n^s$ is statistically zero-mean (Sect. 2.2 and Appendix C), so is ${\mathbf{\Delta }}_n^s-\overline{{\mathbf{\Delta }}_n^s}$ so that Eqs. (26) and (27) become

$$c_{ij}^P\simeq {\sum }_{n,m=1}^{\mathrm{\infty }}{a}_{i,n}^s{a}_{j,m}^s(\mathrm{Cov}({\omega }_n^s,{\omega }_m^s)+\mathrm{Cov}({\mathrm{R}\mathrm{e}\mathrm{\Delta }}_n^s,{\mathrm{R}\mathrm{e}\mathrm{\Delta }}_m^s)).$$(36)

$$c_{ij}^{\psi }={\sum }_{n,m=1}^{\mathrm{\infty }}{a}_{i,n}^s{a}_{j,m}^s\mathrm{Cov}({\mathrm{I}\mathrm{m}\mathrm{\Delta }}_n^s,{\mathrm{I}\mathrm{m}\mathrm{\Delta }}_m^s).$$(37)

Covariances of ${\mathrm{R}\mathrm{e}\mathrm{\Delta }}_n^s$ and ${\mathrm{I}\mathrm{m}\mathrm{\Delta }}_n^s$ involve the statistical properties of ${\mathrm{\Delta }}_n^s$ that combine the fluctuations of the magnetic field and the optical properties of aligned grains within the light cone. In Appendix C, we present a toy model for a turbulent magnetic field close to equipartition, which successfully explains various statistical properties of dust polarization (Planck Collaboration 2018 XII 2020), and demonstrate that

$$⟨{\mathrm{R}\mathrm{e}\mathrm{\Delta }}_n^s{\mathrm{R}\mathrm{e}\mathrm{\Delta }}_m^s⟩\simeq \overline{\eta }⟨{\mathrm{I}\mathrm{m}\mathrm{\Delta }}_n^s{\mathrm{I}\mathrm{m}\mathrm{\Delta }}_m^s⟩\propto \frac{1}{p_0^2},$$(38)

where p₀ = P₀/I₀ is the polarization fraction at the reference frequency with I₀ the total intensity, and $\overline{\eta }\simeq 0.7\pm 0.2$ is a statistical parameter³ characterizing our turbulent magnetic field model. This property is the spectral analog of what is shown in Planck Collaboration 2018 XII (2020) for the spatial variations of polarization angles: dispersion in polarization angles anti-correlate with p₀.

In calculating the covariances $\mathrm{Cov}({\mathrm{I}\mathrm{m}\mathrm{\Delta }}_n^s,{\mathrm{I}\mathrm{m}\mathrm{\Delta }}_m^s)$ through Eq. (25), the scaling of $P_0^2$ with $p_0^2$ exactly cancels the dependence of $⟨{\mathrm{I}\mathrm{m}\mathbf{\Delta }}_n^s\times {\mathrm{I}\mathrm{m}\mathbf{\Delta }}_m^s⟩$ with $1/p_0^2$, yielding $\mathrm{Cov}({\mathrm{I}\mathrm{m}\mathbf{\Delta }}_n^s,{\mathrm{I}\mathrm{m}\mathbf{\Delta }}_m^s)\propto ⟨I_0^2⟩$ and therefore

$$\boxed{{\displaystyle c^{\psi }\propto ⟨I_0^2⟩.}}$$(P3)

Unlike c^ε, the variance c^ψ of the polarization angle residuals is predicted, unexpectedly, to scale with the mean square of the total intensity in the mask, with no dependence on the mean polarization fraction.

In the first term of the right-hand side of Eq. (36), we recognize, from Eq. (29), the variance of dust emissivity $c_{ij}^{\epsilon }$. The second term is, on average and according to Eqs. (37) and (38), proportional to $c_{ij}^{\psi }$, leading to the statistical relationship:

$$\boxed{{\displaystyle c^P\simeq c^{\epsilon }+\overline{\eta }{c}^{\psi }.}}$$(P4)

The variance in polarized intensity c^P not only contains the variance c^ε of the dust polarization emissivity from one region of the sky to another, but also includes an additional variance due to distortions of the SED within the light cone, which can be estimated from the variance of the polarization angle residuals. Prediction P4 implies that c^P, unlike c^ε, does not converge toward zero when the polarization fraction p₀ tends towards zero. It also suggests a recipe for estimating the unknown variance c^ε of the dust polarized emissivity:

$${\hat{c}}^{\epsilon }\equiv c^P-\overline{\eta }{c}^{\psi }.$$(39)

Since c^ε ≥ 0, this implies an upper limit for the ratio c^ψ/c^P (prediction P5):

$$\boxed{{\displaystyle \frac{c^{\psi }}{c^P}\le \frac{1}{\overline{\eta }}.}}$$(P5)

We use this last property to test whether our hypothesis HM6, which states that the signal is dominated by dust, holds, aside from contamination by CMB E-modes, noise, and systematics that we attempt to remove using simulations.

We conclude by analyzing the consequences for the covariances of the quasi-perfect correlation between the residual maps predicted in our model. Eqs. (26) and (27) show that the frequency dependence of any covariance c_ij is encoded in the parameters a_i,n and a_j,n. These parameters are common to all covariances and remain well-correlated across frequencies, as shown in Sect. 3.1. This implies that the ratio between two distinct covariances sharing the same frequency pair is independent of the particular pair (i, j) of channels considered (prediction P6):

$$\boxed{{\displaystyle \frac{c_{ij}^{\psi }}{c_{ij}^P}={\mathrm{C}}^{\mathrm{t}\mathrm{e}}\mathrm{\forall }(i,j).}}$$(P6)

We conclude with our final prediction, a null test resulting from the strong correlation of the spectral gradients across frequencies. The matrix Im c^P (Eq. (35)) involves the cross-product a_i,n a_j,m − a_i,m a_j,n, which is a small factor due to the strong correlation of the a_i,n parameters across frequencies. Hence, within the uncertainties, we expect (prediction P7)

$$\boxed{{\displaystyle \mathrm{Im}{\mathbf{c}}^{\mathbf{P}}\simeq 0.}}$$(P7)

In the following section, we compare our predictions with Planck polarization data.

4 Model validation using Planck polarization data

We applied our dust polarized emission model to analyze Planck data. We extend the power-spectrum analysis presented by Ritacco et al. (2023) to our model framework by measuring the total variance of the residual maps at the map level. The SRoll2 version of the Planck data is introduced in Sect. 4.1, and its associated simulations, used to debias our estimations, are described in Sect. 4.2. In the subsequent sections, we present the mean complex polarized SED in Sect. 4.3 and test our model predictions in Sect. 4.4.

4.1 Presentation of Planck maps

Ritacco et al. (2023) performed a power spectra analysis of the Planck data to characterize spatial variations of the SED of dust polarized emission in terms of the E and B–modes. We follow this analysis using Planck full- and half-mission sky maps⁴ at 100, 143, 217, and 353 GHz and end-to-end data simulations⁵ produced by the SRoll2 software, which corrects the main systematic effects that impact polarization on large angular scales (Delouis et al. 2019). We applied the scaling factors ${\tilde{\rho }}_i$ at frequency ν_i, listed in Table 1 of Ritacco et al. (2023), which represent corrections to the polarization efficiencies. We also applied dust color corrections listed in this table.

We used HEALPix pixelization (Górski et al. 2005) at N_side=32 (map pixel size of 1.8°). The resolution of the original maps was degraded applying a cosine filter in harmonic space (see formula (1) in Ritacco et al. 2023). We followed the procedure described in Sect. 2.2 of Ritacco et al. (2023) to subtract the synchrotron polarized emission from the maps at 100 and 143 GHz. We verified that uncertainties in the assumed spectral index of polarized synchrotron emission (β_S = −3.19 ± 0.07, Ritacco et al. 2023) do not affect our results. We refer to the obtained maps as raw maps.

We applied the complex polarization formalism of Sect. 2 by computing the variances of residual maps. Our methodology follows that of Ritacco et al. (2023). We used half-maps (two independent versions of the same map) to compute covariances, c_raw, for the Planck data (see Appendix D). Our results cannot be compared directly with those of Ritacco et al. (2023) because the transformation from Q and U to E and B is non-local. For most of our data analysis, we used a single mask representing f_sky = 97.3% of the sky, similar to the 97% mask of Ritacco et al. (2023). Our mask at N_side=32 (Fig. 3) was obtained by removing pixels at Galactic latitude |b| ≤ 3° in the inner Galaxy (−90° ≤ l ≤ 90°), as well as all pixels within 4° of the Crab pulsar. For certain plots, lower values of f_sky were used (see Fig. 3); these masks were obtained by removing pixels from the 97.3% mask, ordered by decreasing optical depth, τ₃₅₃, at 353 GHz.

In Section 5.3, we briefly comment on the comparison of our results obtained with SRoll2 maps and those obtained by applying the same methodology to the Planck PR4 maps. The PR4 Planck dataset is the latest and final Planck data release, available from the Planck Legacy Archive⁶. These Planck maps were produced by applying the NPIPE pipeline to the intensity and polarization data from both Planck-LFI and Planck-HFI (Planck Collaboration Int. LVII 2020), to create I, Q, and U calibrated maps for each frequency band. The PR4 data release comes with a comprehensive set of “end-to-end” Monte Carlo simulations processed with NPIPE, aimed at identifying and quantifying potential biases and the errors introduced by the pipeline.

Fig. 3 HEALPix masks with growing fsky (Nside = 32): 85% (dark blue), 90% (light blue), 92% (green), 95% (orange), and 97.3% (brown). The 97.3% mask additionally excludes the inner Galactic plane within 3 degrees of latitude and pixels within 4 degrees of the Crab pulsar (gray).

4.2 Data simulations, debiasing, and error-bars estimation

To assess uncertainties associated with Planck satellite noise, including systematics, and the CMB signal, we computed a set of simulated maps based on various models of dust polarization from PYSM 3 (Thorne et al. 2017; Zonca et al. 2021; Group 2025). The d1 model is built from the Planck 353 GHz Legacy maps, extrapolated at lower frequencies with an MBB spectrum, using the spatially varying temperature and spectral index obtained from the Planck intensity data with the Commander code (Planck Collaboration 2015 X 2016). The d10 model is a refined version of d1, in which the templates for amplitude and spectral parameters have been extracted from the GNILC analysis of the Planck data (Planck Collaboration Int. XLVIII 2016; Planck Collaboration 2018 XII 2020), with smaller scales included using the so-called logarithm of the polarization fraction tensor formalism (for more details, see Group (2025)). Finally, the d12 model consists of up to six superposed layers of modified black bodies, each with spectral parameters drawn from Gaussian distributions (Martínez-Solaeche et al. 2018). The dust model maps were combined with the SRoll2 simulations of data noise and systematics (Delouis et al. 2019) and with CMB maps computed from the CMB power spectrum for the ΛCDM fiducial Planck model (Planck Collaboration 2018 VI 2020) as described in Sect. 3.2 of Ritacco et al. (2023). Overall, at each of the four Planck frequencies, we obtained N_sims = 200 sets of simulated Q and U full and half-mission maps, with independent realizations of data uncertainties and CMB anisotropies. The simulations were treated identically to the Planck data.

Our model is designed to study the frequency dependence of the variances, c, of dust polarized emission. The Planck data must be corrected for the bias introduced by noise, systematics, and the CMB in the raw data squared quantities, c_raw. For simulation number k, the bias introduced by noise, systematics, and the CMB on any squared quantity, c, is estimated as

$$c_{\text{bias }}(k)\equiv c_{\text{sims }}(k)-c_{\mathrm{d}\mathrm{m}},$$(40)

where c_sims(k) and c_dm represent the covariances of simulation k and of the chosen PYSM 3 dust model, respectively. Taking the raw biased data value c_raw, we obtain, for each simulation k, a debiased estimate, $\hat{c}(k)$, for the covariance, c(k), in the Planck SRoll2 data:

$$\hat{c}(k)\equiv c_{\mathrm{r}\mathrm{a}\mathrm{w}}-c_{\mathrm{b}\mathrm{i}\mathrm{a}\mathrm{s}}(k)=c_{\mathrm{r}\mathrm{a}\mathrm{w}}-c_{\mathrm{s}\mathrm{i}\mathrm{m}\mathrm{s}}(k)+c_{\mathrm{d}\mathrm{m}}.$$(41)

The debiased value, c, and its uncertainty, σ_c, are taken as the mean and standard deviation of $\hat{c}(k)$ over its 200 values, respectively.

From now on, we present only the results obtained using d1, our reference model. The results obtained using d10 and d12 are presented in the Appendix F; they do not differ significantly from those presented below, ensuring the robustness of our conclusions against the dust model used in our debiasing procedure.

4.3 Mean complex polarized SED

We extended earlier studies (Planck Collaboration Int. XXX 2016; Planck Collaboration 2018 XI 2020; Ritacco et al. 2023) to the complex plane using the cross-correlation of Planck maps, and determined the mean, debiased, complex SED ${\overline{\mathbf{r}}}_i$ of dust polarization, normalized to our reference channel,

$${\overline{\mathbf{r}}}_i\equiv \frac{{\tilde{\rho }}_i{\tilde{\rho }}_0{⟨{\mathbf{P}}_i{\mathbf{P}}_0^{\star }⟩}_{\mathrm{r}\mathrm{a}\mathrm{w}}-\frac{1}{N_{\text{sims }}}\sum _{k=1}^{N_{\text{sims }}}{⟨{\mathbf{P}}_i{\mathbf{P}}_0^{\star }⟩}_{\mathrm{s}\mathrm{i}\mathrm{m}k}+{⟨{\mathbf{P}}_i{\mathbf{P}}_0^{\star }⟩}_{\mathrm{d}\mathrm{m}}}{{\tilde{\rho }}_0^2{⟨{\mathbf{P}}_0{\mathbf{P}}_0^{\star }⟩}_{\mathrm{r}\mathrm{a}\mathrm{w}}-\frac{1}{N_{\text{sims }}}\sum _{k=1}^{N_{\text{sims }}}{⟨{\mathbf{P}}_0{\mathbf{P}}_0^{\star }⟩}_{\mathrm{s}\mathrm{i}\mathrm{m}k}+{⟨{\mathbf{P}}_0{\mathbf{P}}_0^{\star }⟩}_{\mathrm{d}\mathrm{m}}},$$(42)

where ${\tilde{\rho }}_i$ is the correction to the polarization efficiency⁷ at frequency ν_i (Ritacco et al. 2023, Table 1). Similarly, we computed ${\mathbf{r}}_i^{\mathrm{d}\mathrm{m}}\equiv {⟨{\mathbf{P}}_i{\mathbf{P}}_0^{\star }⟩}_{\mathrm{d}\mathrm{m}}/{⟨{\mathbf{P}}_0{\mathbf{P}}_0^{\star }⟩}_{\mathrm{d}\mathrm{m}}$ for the dust model, ${\mathbf{r}}_i^{\text{sim}}(k)\equiv {⟨{\mathbf{P}}_i{\mathbf{P}}_0^{\star }⟩}_k/{⟨{\mathbf{P}}_0{\mathbf{P}}_0^{\star }⟩}_k$ for simulation k, and ${\overline{\mathbf{r}}}_i^{\text{sims}}\equiv {\sum }_{k=1}^{N_{\text{sims}}}{⟨{\mathbf{P}}_i{\mathbf{P}}_0^{\star }⟩}_k/{\sum }_{k=1}^{N_{\text{sims}}}{⟨{\mathbf{P}}_0{\mathbf{P}}_0^{\star }⟩}_k$ for the mean of all simulations.

The uncertainty on the real and imaginary components of the mean simulated SED ${\overline{\mathbf{r}}}_i^{\text{sims}}$ is taken as the standard deviation of the corresponding real and imaginary components of the 200 simulation values ${\mathbf{r}}_i^{\text{sim}}(k)$. The uncertainty in the observed SED ${\overline{\mathbf{r}}}_i$ must also include the uncertainty in the recalibration factor, ${\tilde{\rho }}_i$, of the Planck polarization data: $\sigma \left({\tilde{\rho }}_i\right)/{\tilde{\rho }}_i=0.5\mathrm{\%}$. We repeated the procedure described in Sect. 4.2 200 times, each time drawing a distinct recalibration factor for each channel, i, from a random Gaussian distribution of expectation ${\tilde{\rho }}_i$ and standard deviation $\sigma \left({\tilde{\rho }}_i\right)=0.005{\tilde{\rho }}_i$. We debiased this realization with simulation k alone to compute the k^th debiased estimate ${\hat{\mathbf{r}}}_i(k)$ of ${\overline{\mathbf{r}}}_i$. The uncertainty in the real and imaginary parts of ${\overline{\mathbf{r}}}_i$ is then taken as the standard deviation of the real and imaginary components of this set of 200 values ${\hat{\mathbf{r}}}_i(k)$, respectively. To test the assumption made in Sect. 2.5 to obtain Eq. (20), we computed the module of ${\overline{\mathbf{r}}}_i\frac{{\overline{\epsilon }}_0}{{\overline{\epsilon }}_i}-1$ for a pivot SED, ${\overline{\epsilon }}_i$, fitted to ${\overline{\mathbf{r}}}_i$, and found it smaller than 0.03 at all frequencies.

Figure 4 shows the spectral dependence of the mean complex polarized SED in the log complex plane, for the raw SRoll2 data, the mean of simulations, the d1 dust model used for debiasing and the resulting debiased SED ${\overline{\mathbf{r}}}_i$. We observe a variation of the polarization angle with frequency in the data, which may appear significant. However, our estimate for the uncertainty on ${\overline{\mathbf{r}}}_i$ does not include the uncertainty in the orientation of the Planck polarization bolometers. The most stringent estimates of these uncertainties come from the analysis of EB cross-spectra carried out to constrain cosmic birefringence (Planck Collaboration Int. XLIX 2016; Minami & Komatsu 2020; Diego-Palazuelos et al. 2022, 2023). Diego-Palazuelos et al. (2022) list their results in Table 1, per bolometer for several Galactic masks. The scatter among Galactic masks and the statistical uncertainties both lie in the range 0.2–0.3°. This agrees with the systematic uncertainty of 0.28° quoted by Planck Collaboration Int. XLIX (2016). The magnitude of this uncertainty reduces the significance of the imaginary component of our mean SED.

While the mean complex SED, ${\overline{\mathbf{r}}}_i$, depends on calibration errors in amplitude (correction factors ${\tilde{\rho }}_i$) and in phase (absolute orientation of bolometers in the focal plane), the definition of map residuals (Eq. (18)), involving the ratio of P_i to ${\overline{\mathbf{r}}}_i$, renders the map residuals, and therefore their covariances, independent of any calibration error, whether in amplitude or in phase. This facilitates the comparison of different frequency maps or datasets with distinct systematics.

Fig. 4 Mean complex polarized SED ${\overline{\mathbf{r}}}_i$ in the log complex plane for our 97% mask: debiased SRoll2 data (solid black line), raw SRoll2 data (dashed-dotted blue line), simulations (mean shown by the dashed red line with error bars representing the standard deviation), and the model d1 (dotted line). Errors bars are based on 200 simulations (see Sect. 4.3). Planck-HFI frequencies are indicated.

4.4 Comparison of model predictions with Planck polarization data

In this section, we compare our model predictions and hypotheses with the Planck SRoll2 data. Using simulations, we have verified that the 0.5% in the recalibration factors ${\tilde{\rho }}_i$ does not have a noticeable effect on our results.

First, we tested model prediction P1 by computing the Pearson correlation coefficient ${\rho }_{ij}=c_{ij}/\sqrt{c_{ii}{c}_{jj}}$, between the residual maps i and j. Our results for the R, R^P, and R^ψ residual maps are presented in Table 1 for all frequency pairs, using our 97.3% mask. The error bars on ρ_ij are calculated from the standard deviation of the simulated values ρ_ij(k), following the approach described in Sect. 4.2. We find that the correlation coefficient between the R^P residual maps is consistent with a perfect correlation. The correlation between the angle residual maps, R^ψ, is strong for the (217,143) pair, but weak for the (143,100) and (100,217) pairs, although with large uncertainties. Therefore, prediction P1 is validated for R^P but not for R^ψ. The correlation coefficient of Re c^P is high but only marginally compatible with 100% for the (217,100) and (143,100) pairs. In the context of Eq. (33), this loss of correlation between the residual maps R_i is a consequence of the loss of correlation between the residual maps $R_i^{\psi }$.

Second, we examined the validity of predictions P2 and P3, which relate c^ε and c^ψ to $⟨P_0^2⟩$ and $⟨I_0^2⟩$, respectively. Fig. 5 shows how the variances, computed at 217 GHz for a high S/N ratio, depend on the polarization fraction, p₀, at 353 GHz. This correlation plot is obtained by binning pixels of the 97.3% mask in bins of p₀ and computing, for each p₀ bin, the variances c at 217 GHzand the mean squared values, $⟨P_0^2⟩$ and $⟨I_0^2⟩$ (a method only applicable at the map level). A model, $c^{\epsilon }\propto ⟨P_0^2⟩$, shown as a dashed orange line, confirms our prediction P2. Another model, $c^{\psi }\propto ⟨I_0^2⟩$ (dashed red line), also confirms the surprising prediction P3 of our model: c^ψ is reasonably proportional to $⟨I_0^2⟩$, rather than $⟨P_0^2⟩$. Consequently, while c^ε follows $⟨P_0^2⟩$ as the polarization fraction tends toward zero, c^P and Re c^P saturate at positive values (prediction P4) due to depolarization effects (the variance of which is quantified here by c^ψ) that persist even at low p₀.

Predictions P5 and P6 are compared to Planck SRoll2 data in Fig. 6, which shows how the ratio $c_{ij}^{\psi }/c_{ij}^P$ varies with f_sky, for all pairs (i, j) of residual maps. The error bars for this ratio are computed from simulations, following the approach described in Sect. 4.2. Prediction P5 ($c^{\psi }/c^P\le 1/\overline{\eta }$, unhashed region) is verified for all frequency pairs except for the variances at 100 and 143 GHz at lower f_sky. This indicates that either the hypotheses of our model or the debiasing method break down at lower f_sky. For the 95 and 97.3% masks, the c^ψ/c^P ratio is, within error bars, almost independent of the frequency pair considered (prediction P6), with the notable exception of the variance at 100 GHz.

Overall, Table 1 and Fig. 6 highlight a limitation of our approach for the 100 GHz channel. In Sect. 6, we discuss to what extent this distinctive behavior can be attributed to dust physics beyond our hypotheses, to residual contamination by a component other than dust and CMB (e.g., synchrotron or CO emission), or to limitations of our debiasing method.

Finally, in Table 2, we test prediction P7, which states that $\mathrm{Im}{\mathbf{c}}_{ij}^{\mathbf{P}}$ should be, within the uncertainties, close to zero for all frequency pairs. This prediction is not verified, even for the frequency pair (217,143) for which the residual maps show the strongest correlation (see Table 1). This result is puzzling, as our model does not predict any correlation between R^ψ and R^P, except for chance correlation. This issue is examined in more detail in the following section.

Fig. 5 Test of predictions P2, P3, and P4: dependence on the observed polarization fraction p0 = P0/I0 at 353 GHz of the variances Re cP, cP, cψ, and cε computed at 217 GHz. Our models for cε (P2) and cψ (P3) are overplotted as dotted orange and red lines, respectively, using the mean value of the slope derived from these equations, multiplied by $⟨P_0^2⟩$ and $⟨I_0^2⟩$ in each bin of p0 = P0/I0.

Fig. 6 Test of predictions P5 and P6 as a function of fsky. The hashed region represents forbidden values according to P5.

5 Inferring the variance of the spectral properties of aligned grains

Polarization observations trace the spectral properties of a specific population of dust grains that are aligned with the magnetic field lines. They can be compared with the overall spectral properties inferred from dust far-infrared (FIR) and submm total emission, which trace all grains, both aligned and unaligned. Assuming a common temperature⁸ for aligned and unaligned grains, Planck Collaboration 2018 XI (2020) found close values for the mean spectral index in polarization and intensity: β_P − β_I = 0.05 ± 0.03 in the high-latitude diffuse ISM. In this section, we extend this work by quantifying the variance of the fluctuations in the spectral properties of aligned grains for the 97.3% mask, where the signal is dominated by the Galactic plane. We trace these variations back to their origin in the fluctuations of β or T, and correlate them with the spectral fluctuations observed in total intensity.

Fig. 7 Comparison of SRoll2 data with the predictions of our model for cP (blue) and cψ (red), for pure β fluctuations (top), and pure T fluctuations (bottom).

5.1 Disentangling fluctuations of dust temperature and spectral index

As demonstrated in Sect. 3.1, the spectral gradients $a_{i,n}^s$ for n = 2 and n = 3 are strongly correlated with $a_{i,1}^s$. By correlating the covariance matrices $c_{ij}^P$ and $c_{ij}^{\psi }$ with $a_{i,1}^{\beta }{a}_{j,1}^{\beta }$ and ${\overline{T}}^4{a}_{i,1}^T{a}_{j,1}^T$ separately⁹, we aimed to identify whether these fluctuations originated mainly from variations in β or T, while simultaneously estimating their amplitude. This analysis is illustrated in Fig. 7 for pure β fluctuations (labeled (top)) and pure T fluctuations (labeled (bottom)). The spectral trends observed for c^P and c^ψ are not compatible with a scenario in which they originate from fluctuations in β, as this scenario predicts a steeper slope. A better agreement is obtained if only the dust temperature T fluctuates. In the Galactic plane, starlight heating may indeed dominate the fluctuations of dust emission properties. For the variance in angles (c^ψ), fluctuations in T may also be favored, although the significant contamination observed in the 100 GHz residual angle map prevents a firm conclusion. To quantify these trends, we performed a χ² minimization using the MPFITCOVAR IDL routine (Markwardt 2009). Our data combine the two data sets, $c_{ij}^P$ and $c_{ij}^{\psi }$, into an array of 12 values. The corresponding 12 × 12 noise covariance matrix was computed from the covariances between the 200 simulated $c_{ij}^P(k)$ and $c_{ij}^{\psi }(k)$ covariance matrices (see Sect. 4.2). Using predictions P2, P3, and P4, we define a linear model for fluctuations in the spectral parameter s ∈ {T, β}, where $c_{ij}^{\psi }=a_{i,1}^s{a}_{j,1}^s{\left({\sigma }_s^{\psi }\right)}^2⟨P_0^2⟩$ and $c_{ij}^P=a_{i,1}^s{a}_{j,1}^s{\left({\sigma }_s^{\epsilon }\right)}^2⟨P_0^2⟩+\overline{\eta }{c}_{ij}^{\psi }$, with ${\sigma }_s^{\psi }\equiv \sqrt{\mathrm{Cov}({\mathrm{I}\mathrm{m}\mathcal{W}}_{\mathbf{1}}^s,{\mathrm{I}\mathrm{m}\mathcal{W}}_{\mathbf{1}}^s)}$ and ${\sigma }_s^{\epsilon }\equiv \sqrt{\mathrm{Cov}({\omega }_1^s,{\omega }_1^s)}$. When fitting with pure fluctuations in T, we obtain a χ² of 15.7 (or 1.57 per degree of freedom). For pure β fluctuations, χ² increases to 31. Removing $c_{100\times 100}^{\psi }$ from the dataset clearly improves the fit for fluctuations in T (χ²/dof = 0.96), but not for fluctuations in β. Changing the value of $\overline{\eta }$ does not affect the χ² of the fit’ it modifies only the value of ${\sigma }_s^{\epsilon }$, while keeping ${\left({\sigma }_s^{\epsilon }\right)}^2+\overline{\eta }{\left({\sigma }_s^{\psi }\right)}^2$ constant. Our best fit provides a reasonable value for the typical dispersion in dust temperature in the Galactic plane, corresponding to ${\sigma }_T^{\epsilon }=(2.7\pm 0.1)\mathrm{K}$ for a reference temperature $\overline{T}=20\mathrm{K}$.

Finally, Fig. 8 shows the frequency dependence of the mixed covariances $c_{ij}^{\psi P}$ (Eq. (28)). This graph confirms that the 143 GHz channel behaves differently from the 217 and 100 GHz channels, possibly indicating contamination by CO emission in these two channels (see Sect. 6). A first-order model (dashed line) based on fluctuations of the dust temperature, $c_{ij}^{\psi P}=a_{i,1}^T{a}_{j,1}^T{\mathrm{\Sigma }}_T^{\psi P}⟨P_0^2⟩$, with ${\mathrm{\Sigma }}_T^{\psi P}\equiv \mathrm{Cov}({\mathrm{I}\mathrm{m}\mathcal{W}}_{\mathbf{1}}^{\mathit{T}},{\mathrm{R}\mathrm{e}\mathcal{W}}_{\mathbf{1}}^{\mathit{T}})$, cannot reproduce the observations. A model based on fluctuations in β also fails. The frequency dependence of the mixed covariances $c_{ij}^{\psi P}$, together with the values of $\mathrm{Im}{\mathbf{c}}_{ij}^{\mathbf{P}}$ (Table 2), warrants further detailed analysis. We defer this to future work.

Fig. 8 Covariances $c_{ij}^{\psi P}$ compared to a first-order model (dashed line) based on fluctuations of the dust temperature alone.

5.2 Correlating fluctuations in intensity and polarization

We next correlated the fluctuations observed in the dust polarized emissivity with the fluctuations observed in total intensity. We used the d1 dust model, based on Commander total intensity maps, as a proxy for the dust emissivity fluctuations in total intensity. This model does not include any significant decorrelation in angles (see Sect. 4.1), which implies that R^ψ(d1) ≃ 0 and c^P(d1) ≃ c^ε(d1). Correlating R^P(SR) from the SRoll2 data with R^P(d1) from the d1 dust model therefore quantifies the variance of dust emissivity common to both residual maps, independently of the variance in angles present in the Planck SRoll2 maps R^P(SR). Figure 9 compares the dependence with the frequency ν_i of the variance $c_{ii}^{\epsilon }$(SR) of the SRoll2 data, of the variance $c_{ii}^{\epsilon }$(d1) = $c_{ii}^P$(d1) of the d1 dust model, and their cross-covariance,

$$c_{ii}^P(\mathrm{d}1\times \mathrm{S}\mathrm{R})\equiv ⟨R_i^P(\mathrm{d}1)\times R_i^P(\mathrm{S}\mathrm{R})⟩.$$(43)

Fig. 9 shows that the variance of emissivity of aligned grains (as traced by c^ε(SR)) and that of all grains (aligned or unaligned, traced by c^P(d1)), are of compatible amplitude, except at 217 GHz. The SRoll2 polarization data follow a T trend, while the d1 intensity data lie between the pure T and pure β trends. The data-model correlation coefficients,

$${\rho }_i^{\epsilon }(\mathrm{d}1\times \mathrm{S}\mathrm{R})=c_{ii}^P(\mathrm{d}1\times \mathrm{S}\mathrm{R})/\sqrt{c_{ii}^{\epsilon }(\mathrm{S}\mathrm{R})\times c_{ii}^P(\mathrm{d}1)},$$(44)

are relatively high (above 85%) at all frequencies. However, the value at 100 GHz must be interpreted with caution, as the value derived for $c_{100\times 100}^{\epsilon }$ through Eq. (39) is underestimated due to contamination by $c_{100\times 100}^{\psi }$, as observed in Fig. 7. Overall, this indicates that the fluctuations observed in the total SED of a given light cone correspond closely to the fluctuations present in the polarized SED of the same light cone.

Fig. 9 Comparison of the spectral dependence of covariances cε from SRoll2 data, cP for the dust model d1, and cP for the cross-product of SRoll2 and d1. The Pearson correlation coefficient ρ(d1, SR) is indicated (purple). Predictions for models with T or β fluctuations are overplotted.

5.3 Results obtained with the Planck PR4 data release

In Appendix E, we present figures and tables obtained using the same methodology but applied to the last official release of Planck data, PR4. Our goal was not to show a preference for one data release over another, but to demonstrate that our methodology can quantify differences between datasets and track their origin. Despite the clear differences in the mean complex SED, our assessment of the seven predictions remains the same. The main difference between these two data releases appears in the spectral dependence of covariances. The PR4 and SRoll2 covariances, c^P and c^ψ, are compatible within uncertainties, except for $c_{217\times 217}^P$ and $c_{100\times 100}^{\psi }$, possibly indicating a distinct handling of contamination by CO emission. The PR4 and SRoll2 c^ψP mixed covariances are compatible for all channel pairs. Contamination of the 100 GHz channel appears to be significant in c^P in the PR4 data, unlike in SRoll2. Finally, while the SRoll2 data shows a preference for pure fluctuations in T, PR4 data, with χ²/dof ≃ 14, are not compatible with a dust model assuming either pure fluctuations in T or in β, even when excluding $c_{100\times 100}^{\psi }$ and $c_{100\times 100}^P$. Our methodology reveals clear differences between the SRoll2 and PR4 data, particularly in the mean properties and variances of the 217 and 100 GHz channels contaminated by CO leakage.

5.4 Dependence on the dust model used for debiasing

The influence of our debiasing method (Sect. 4.2), applied with the d1 dust model, is tested in Appendix F using the d10 and d12 PySM 3 dust models. This effect is weak and does not alter our conclusions. However, the fact that no dust model includes a frequency-rotation of polarization angles at the pixel level (with the exception of d12, which does not appear to behave realistically, as indicated by the negative Pearson coefficients in Fig. F.4) is clearly problematic. One of the goals of this article is to propose new diagnostics for developing such models. Coupling these diagnostics with model-building tools based on moment expansion, as proposed by Vacher et al. (2025), could provide a powerful combination to provide new, realistic, and complex foreground models anchored in data.

6 Discussion

Our model, based on simple hypotheses, makes strong predictions and testable predictions. Its success enables the measure of the variance in the spectral properties of aligned grains. Its limitations provide new insights into the origins of the imperfect correlation observed in polarization angle residual maps.

6.1 Insights into the properties of aligned grains

Correlating fluctuations in the polarized SED with those observed in total intensity, we find a strong, though not perfect, correlation (ρ > 85%; Fig. 9). This behavior is expected. What we call dust grains are actually a set of grain populations with distinct composition, optical properties, and alignment properties (Shiu et al. 2024; Siebenmorgen 2023). Typically, small grains and possibly unaligned large carbonaceous grains do not emit in polarization (Kim & Martin 1995; Chiar et al. 2006). Maps of total intensity collect the thermal emission from all grains, whether aligned or not with the magnetic field. Because the impact of magnetic fields on fluctuations in the polarized SED is statistically mitigated within our methodology, the strong correlation between intensity and polarization residual maps leads us to conclude that the emission from aligned grains is highly dominant in Planck channels, compared with that from unaligned grains.

As a second result on dust properties, we find that the frequency-dependence of the variance of these fluctuations is consistent only with variations in dust temperature T, showing no indication of a variation in the dust spectral index β. This result may be specific to the low Galactic latitudes that dominate our signal, where the light cone samples the systematic increase in the radiation-field intensity toward the molecular, star-forming ring of the Galaxy. However, we do not draw firm conclusions from this preliminary result. Firstly, it is at odds with those drawn from analyses of FIR and submm total-intensity SEDs (Ysard et al. 2013; Planck Collaboration 2013 XI 2014; Planck Collaboration Int. XVII 2014) and from a recent evolution of dust models (Ysard et al. 2024), and secondly, because it is not confirmed using the PR4 version of the Planck data (see Sect. 5.3). This apparent contradiction requires further investigation.

6.2 Investigating the origin of the imperfect correlation of polarization angle residuals

Our model predicts a perfect correlation between the polarization SEDs and angles residuals (but not of the polarized SEDs themselves). Analyzing the Planck SRoll2 data with this model, we find that fluctuations of the submm polarized SED are indeed strongly correlated (ρ^P > 90%; Table 1) across frequencies for the polarized intensity, P_ν, but only weakly correlated for the polarization angles, ψ_ν. The level of correlation between residual maps is lower with PR4, but remains similar (see Appendix E). This raises the question of which physical processes or components could be responsible for the weaker-than-expected correlation observed between the polarization angle residual maps.

Dust itself could be the origin of part of this loss of correlation. First, we restricted our MBB modeling to the simple situation (hypothesis HM3), where only T or β vary in the sky, not both. Observations show that, in total intensity, both T and β vary within any given region of the sky. The possible correlation between T and β has also been the subject of a long series of studies (Dupac et al. 2003; Anderson et al. 2010; Planck Collaboration 2011 XXIV 2011; Planck Collaboration 2013 XI 2014). All of these studies are based on intensity rather than polarization maps, which motivated our simple hypothesis HM3. Nevertheless, the impact of simultaneous variations, possibly correlated or anti-correlated, of T and β on our analysis remains to be explored. Second, the predicted 100% correlation between frequencies for fluctuations in T derives from the law of blackbody thermal emission and is therefore as strict as a physical law. In contrast, fluctuations in β represent a simplification of what may truly occur. This simplification assumes that variations in the FIR and submm dust emissivity, Q_abs(ν), are perfectly correlated across frequencies. Laboratory data do not support this assumption (Demyk et al. 2022). The prediction of a 100% correlation for fluctuations in β is therefore not derived from solid physics, but from our simplifying hypothesis HA1, the MBB model. Finally, through the same hypothesis HA1, our study neglects refinements of dust models, such as variations in dust cooling efficiency (Ysard et al. 2024), the magnetic dipole emission of dust grains (Draine & Hensley 2013), which may produce polarized emission around 100 GHz with a polarization direction perpendicular to the magnetic field, or the polarized emission by a population of very large (up to 1 μm) grains that may have distinct dust properties (Siebenmorgen 2023). In summary, the predicted 100% correlation of map residuals between frequencies strongly relies on our hypothesis HA1, and is therefore questionable.

Within the framework of our hypothesis HM6, we also ignore any component other than dust, the CMB, and Planck noise and systematics. The removal of synchrotron polarized emission from the 143 and 100 GHz maps (see Sect. 4.1) may have left some residuals. There is also a significant component that is not included in the simulations, as it is difficult to model: the possible leakage of CO line emission into polarized Stokes maps Q and U, in the 100 and 217 GHz channels. In the Galactic plane, CO emission arises from various regions of the Galaxy, with distinct radial velocities and therefore distinct CO line emission frequencies. Because each Planck HFI bolometer that measures linear polarization has its own spectral transfer efficiency, a leakage from intensity to polarization is observed. This artificial polarization, absent from the 143 GHz channel, contributes to the observed decorrelation of polarization angles, as well as to the distortion in the spectral dependence of their covariances shown in Fig. 7. It can also affect the spectral dependence of the mean SED (Fig. 4). Finally, CO molecules may also have an unexpected effect on the measurement of dust polarization by Planck. Anisotropic resonant scattering of background dust polarized emission by foreground aligned CO molecules (Houde et al. 2013; Hezareh et al. 2013) can induce a rotation of the dust polarization signal, potentially affecting the correlation of angles in the 217 and 100 GHz channels, while leaving the 143 GHz channel uncontaminated.

7 Summary and conclusions

We have presented an analytical framework for interpreting variations in the submillimeter polarized SED in terms of fluctuations in spectral properties, T and β, of aligned grains across the sky and within the light cones of observation. Following Vacher et al. (2023b), our approach is based on complex quantities, allowing for a consistent treatment of the fluctuations in the spectral parameters of the polarized SED and of the spectral rotation of angles, isolated in distinct residual maps. When either T or β fluctuates, our model predicts an almost perfect correlation across frequencies between all dust polarization residual maps. The decorrelation observed in the Stokes parameters Q_ν and U_ν translates into a perfect frequency-correlation of P_ν in the complex plane.

Following Ritacco et al. (2023), we validated our model using Planck SRoll2 data (Delouis et al. 2019), with a f_sky = 97.3% mask excluding the inner Galactic plane. To characterize fluctuations in dust polarization properties, covariances of residual maps were computed at the map level, and debiased from contamination by the CMB, noise, and systematics using 200 end-to-end simulations. We find that polarized-intensity residual maps are strongly correlated across frequencies, while polarization angle residual maps are weakly correlated for the (217,143) and (217,100) GHz pairs and uncorrelated for the (143,100) GHz pair. This departure from the predicted perfect correlation between residual maps is discussed in the context of dust physics and data analysis. Polarized SED fluctuations are dominated by variations in the dust temperature (σ_T ~ 2.7 K, assuming a mean temperature $\overline{T}=20\mathrm{K}$ for aligned grains), excluding any significant contribution from fluctuations in the dust polarization spectral index, β. In our analysis, dominated by the Galactic plane, fluctuations in the dust emissivity observed in the polarization SED at all frequencies are strongly correlated with those traced by Commander PR2 maps characterizing the dust SED in total intensity (Planck Collaboration 2015 X 2016), implying that emission from aligned grains greatly exceeds emission by non-aligned grains in the submillimeter. These results, based on SRoll2 maps, are compatible with those obtained using the latest official Planck release, PR4, following the same methodology, with the notable exception that the spectral dependence of the covariances computed with PR4 data is not compatible with a dust model assuming pure fluctuations in either T or β.

Observations of dust emission in polarization provide access to a new dimension of dust evolution, inaccessible to unpolarized observations: the statistics of dust properties variations within the light cone (2D within the beam and 3D along the line of sight). This information appears in what we call the “spectral rotation of polarization angles.” Our model captures a significant fraction of the decorrelation in angles found in Planck data, as well as the variance observed in dust unpolarized-emission SEDs. The spatial variation of the SED induced by rotation of the polarization angle along the line of sight is a critical aspect that could strongly affect the detection of CMB B-modes and any other attempts to detect cosmic birefringence. Our methodology, which is independent of calibration errors and can be used with any mask, enables correlation analyses with tracers of the ISM to trace the origin of decorrelation. This work is useful and timely for the development of component-separation methods for current and future generations of CMB experiments. It will aid the development of new, more realistic sky simulations that include SED decorrelation with amplitudes consistent with Planck data.

Acknowledgements

V.G. thanks M. Houde for interesting discussions. L.V. acknowledges partial support by the Italian Space Agency LiteBIRD Project (ASI Grants No. 2020-9-HH.0 and 2016-24-H.1-2018), as well as the Radio-ForegroundsPlus Project HORIZON-CL4-2023-SPACE-01, GA 101135036 and through the Project SPACE-IT-UP by the Italian Space Agency and Ministry of University and Research, Contract Number 2024-5-E.0. A.R. acknowledges financial support from the Italian Ministry of University and Research – Project Proposal CIR01_00010. A.B. acknowledges financial support from the INAF initiative “IAF Astronomy Fellowships in Italy” (grant name MEGASKAT).

References