Preserved original model output. Check claims against the independent audit →

Dedusting CMB B-modes with 3D dust, grain-size physics, and the sky-projected magnetic field

Polarized thermal emission from Galactic dust is the leading astrophysical contaminant of primordial CMB B-modes at frequencies $\nu\gtrsim 70\,\mathrm{GHz}$. Component-separation pipelines that treat each pixel (or each power-spectrum bin) as a single modified blackbody (MBB) with a spatially constant spectral index $\beta_d$ fail when the effective dust spectral energy distribution (SED) and/or the polarization angle $\psi$ become frequency-dependent. That failure is what the CMB literature calls spectral-index (or frequency) decorrelation. The physical origin is line-of-sight and beam averaging of dust with spatially varying temperature $T_d$, opacity index $\beta_d$, grain-size mix, alignment efficiency, and magnetic-field orientation. Three independent observational ingredients—3D maps of the dust column, knowledge of the grain size distribution and alignment physics, and maps of the plane-of-sky magnetic field—do not merely “help” in a qualitative sense. They supply the missing third dimension and the missing microphysics that turn an underconstrained 2D extrapolation into a layered, physically motivated forward model. This essay states that argument operationally: what is being measured, why a one-MBB-per-pixel model is biased, how each ingredient contributes on its own, how they combine into a concrete dedusting strategy, and where the residual systematics remain.

Throughout, Stokes parameters $(I,Q,U)$ are used in the IAU convention, the linear polarization fraction is

and the polarization angle (electric-vector position angle) is

For optically thin thermal dust, the polarization direction is orthogonal to the plane-of-sky projection of $\mathbf{B}$. The dust intensity is written as a modified blackbody

with $B_\nu(T)$ the Planck function and a conventional reference $\nu_0=353\,\mathrm{GHz}$. Polarized emission from a single homogeneous volume element is

where $p_0$ is the intrinsic polarizing efficiency of the aligned grain population, $\gamma$ is the angle between $\mathbf{B}$ and the plane of the sky, and $\psi$ is set by the plane-of-sky field orientation. The factor $\cos^2\gamma$ is the geometric projection of aspherical grain emission (Lee & Draine 1985). As soon as several such volume elements with different $\{\beta_d,T_d,p_0,\gamma,\psi\}$ occupy the same beam, neither $p$ nor $\psi$ nor the effective $\beta_d$ is frequency-independent.

1. Polarized thermal dust and why B-mode / tensor-to-scalar-ratio studies care

Inflationary gravitational waves source a curl (B-mode) pattern in CMB polarization on degree scales, with amplitude parameterized by the tensor-to-scalar ratio $r$ at a pivot $k=0.05\,\mathrm{Mpc}^{-1}$. The present observational bound is $r_{0.05}<0.036$ at 95% confidence from BICEP/Keck through the 2018 season (BK18), tightening to $r<0.032$ when Planck PR4, BAO, and CMB lensing are combined (BICEP/Keck Collaboration 2021; Tristram et al. 2022). LiteBIRD targets a total uncertainty $\delta r\sim 10^{-3}$, including foreground residuals; CMB-S4 aims at a $5\sigma$ detection of $r=0.003$ in its deep patch. At these levels the Galactic polarized foreground is not a small correction. In the BICEP/Keck field, Galactic dust B-modes at $150\,\mathrm{GHz}$ already exceed the $r=0.01$ tensor spectrum by an order of magnitude (BICEP2/Keck Collaboration 2018). Over most of the high-latitude sky the dust B-mode power at $353\,\mathrm{GHz}$ is larger still.

Thermal dust emission is the long-wavelength tail of radiation from large interstellar grains in equilibrium with the interstellar radiation field. The mean SED in both intensity and polarization is well approximated by an MBB over the Planck HFI range. Planck intermediate results XXII found

for a mean colour temperature $T_d=19.6\,\mathrm{K}$ on intermediate-latitude patches (Planck Collaboration Int. XXII 2015). Later Planck 2018 analyses quote a slightly lower polarized index, $\beta_d\simeq 1.53\pm 0.02$ at the same $T_d$. The $\sim 0.08$ offset between polarized and total-intensity indices is statistically significant and is a first hint that the grains that dominate $I$ are not identical to those that dominate $Q$ and $U$. The same paper found that the dust polarization fraction decreases by $(21\pm 6)\%$ from $353$ to $70\,\mathrm{GHz}$, again inconsistent with a single grain population of frequency-independent $p$.

The polarization itself is produced because aspherical grains align with their short axes preferentially parallel to $\mathbf{B}$. Radiative alignment torques (RATs) are the leading mechanism in the diffuse ISM (Dolginov & Mitrofanov 1976; Draine & Weingartner 1996, 1997; Lazarian & Hoang 2007; Andersson, Lazarian & Vaillancourt 2015). Alignment is strongly size-dependent: only grains larger than a critical alignment radius $a_{\mathrm{align}}$ (typically $\sim 0.05\,\mu\mathrm{m}$ in the diffuse ISM) rotate suprathermally and remain aligned against gas randomization. Consequently the polarized submillimetre SED is weighted toward the large-grain tail of the size distribution, whereas total intensity also receives emission from smaller, hotter, poorly aligned grains.

Planck 353 GHz $Q,U$ maps remain the highest-signal-to-noise full-sky tracer of polarized dust. They reveal several statistical facts that any B-mode foreground model must reproduce:

• A maximum observed polarization fraction $p_{353}^{\max}\sim 0.20$–$0.22$, which is a lower bound on the intrinsic grain polarizing efficiency after line-of-sight depolarization (Planck Collaboration XII 2018).
• An $E/B$ power asymmetry $C_\ell^{EE}/C_\ell^{BB}\simeq 2$ on intermediate scales, together with a positive $TE$ correlation, both of which arise because filamentary density structures are preferentially aligned with $\mathbf{B}_{\mathrm{POS}}$ (Planck Collaboration Int. XXXVIII 2016; Clark et al. 2015).
• A tight anti-correlation between $p$ and the polarization-angle dispersion function $\mathcal{S}$, consistent with depolarization by a turbulent field rather than by spatially variable grain physics (Planck Collaboration Int. XIX 2015; Planck Collaboration XII 2018).
• Residual spatial variation of the polarized SED at the level relevant to $\ell\sim 80$ B-mode analyses (Planck Collaboration Int. L 2017; Ritacco et al. 2023).

Component separation for $r$ is usually done in one of three families: (i) parametric pixel-domain fits (Commander, FGBuster), in which each pixel is an MBB plus a synchrotron power law plus CMB; (ii) internal linear combination / needlet ILC methods, which assume that the CMB is the unique component with a blackbody SED; (iii) template multi-frequency likelihoods at the power-spectrum level, as in BICEP/Keck, which scale a 353 GHz dust template with an MBB and a single $\beta_d$ (sometimes with a decorrelation parameter $\mathcal{R}_\ell$). All three families become biased if the true dust SED is not a single MBB, or if $Q$ and $U$ do not share the same frequency scaling. That is the operational problem of spectral-index decorrelation.

2. What spectral-index decorrelation means operationally

2.1 Two distinct, observationally mixed phenomena

The phrase “spectral-index decorrelation” is used in two related but not identical ways.

Spatial (2D) decorrelation. The dust SED varies across the sky. If $\beta_d(\hat n)$ or $T_d(\hat n)$ is not constant, a map at frequency $\nu_1$ is not a scalar multiple of the map at $\nu_2$. At the power-spectrum level this is quantified by the correlation ratio

A single MBB with uniform $\beta_d$ implies $\mathcal{R}_\ell=1$. Gaussian spatial variations of $\beta_d$ with rms $\Delta\beta_d$ produce

on the scales where the $\beta_d$ map is coherent (Planck Collaboration Int. L 2017). Planck Int. L reported $\mathcal{R}_\ell^{BB}(217,353)$ departing from unity at $\ell\gtrsim 50$ in some large high-latitude regions, with an indicative $\Delta\beta_d^{\mathrm{P}}\sim 0.07$. Subsequent analyses (Sheehy & Slosar 2018; Planck Collaboration XI 2020) found $\mathcal{R}_\ell^{BB}(217,353)$ consistent with 1 within Planck noise over large sky fractions, but they explicitly cautioned that the ratio can be spatially variable and that the effect is expected to grow toward smaller scales and toward larger frequency lever arms (e.g. 95 vs 353 GHz). A correlation ratio of $\mathcal{R}_\ell^{BB}(150,353)=0.85$, scaled as $[\ln(\nu_1/\nu_2)]^2$, is already enough to bias a BICEP/Keck-like $r$ analysis at a level comparable to current statistical errors (Planck Collaboration Int. L 2017).

Line-of-sight (LOS) decorrelation. Even a single pixel can fail to be a single MBB. If the LOS intercepts two (or more) clouds with different SEDs and misaligned $\mathbf{B}_{\mathrm{POS}}$, the Stokes parameters add:

Because each cloud’s polarized intensity ratio

is frequency-dependent whenever $T_1\neq T_2$ or $\beta_{d,1}\neq\beta_{d,2}$, the net polarization angle

rotates with frequency, and the net polarization fraction is depolarized relative to either cloud. This is the Tassis & Pavlidou (2015) mechanism. It is not a small higher-order correction: for two clouds with $\mathbf{B}_{\mathrm{POS}}$ misaligned by $\alpha\gtrsim 60^\circ$ and temperatures differing by several kelvin, $\psi$ can rotate by tens of degrees between 353 and 150 GHz, and $p$ can change by tens of percent. In the extreme $\alpha\to 90^\circ$ case, $\psi_{\mathrm{tot}}$ flips by $90^\circ$ when $\rho_\nu$ crosses unity, and $p_{\mathrm{tot}}$ passes through zero.

LOS decorrelation is not equivalent to a spatially varying but frequency-independent $\beta_d$ per pixel. A pixel-wise MBB fit still assumes that $Q$ and $U$ share a common SED, i.e. that $\psi$ is independent of $\nu$. LOS mixing of misaligned clouds violates that assumption: the effective spectral indices of $Q$ and of $U$ differ. Parametric component separation that fits a single $(\beta_d,T_d)$ to $(Q,U)$ jointly therefore leaves a frequency-dependent residual whose $E/B$ content is not that of the CMB.

Pelgrims et al. (2021) provided the first direct evidence that this LOS mechanism is present in Planck data. Using H I to count clouds per sightline (Panopoulou & Lenz 2020) and H I filament orientation to estimate $\mathbf{B}_{\mathrm{POS}}$ per cloud (Clark & Hensley 2019), they selected a target sample of high-latitude pixels with $N_c\gtrsim 1.5$ and cloud–cloud misalignment $\gtrsim 60^\circ$, and a control sample of single-cloud sightlines. The target sample shows a statistically significant excess of 353-to-217 GHz polarization-angle differences and extra depolarization, with a significance that increases after residual Planck systematics are reduced (SRoll2). An intrinsic cloud-to-cloud variation of $\sim 15\%$ in the 353/217 polarized intensity ratio is sufficient to reproduce the measured effect. That 15% is a $\sim 10\%$ variation in the polarized SED shape between 217 and 353 GHz (Hensley, McBride & Bull, related analyses). The effect is stronger for more misaligned fields, as the two-cloud model predicts.

2.2 Why a single MBB per pixel fails

The MBB is a nonlinear function of $\beta_d$ and $T_d$. Averaging MBBs—along the LOS, inside the beam, or when computing spherical-harmonic coefficients over a large mask—does not return an MBB. Chluba, Hill & Abitbol (2017) formalized the resulting SED distortions as a moment expansion about a fiducial $(\bar\beta_d,\bar T_d)$:

where $\Theta(\nu,T)=(x/T)\,e^x/(e^x-1)$ and $x=h\nu/kT$. The first-order moments $\omega_1^\beta$ and $\omega_1^T$ are proportional to the column-weighted means of $\Delta\beta_d$ and $\Delta T_d$; the second-order moments encode variances and $\beta$–$T$ covariances. In the Rayleigh–Jeans regime relevant to CMB channels, $\beta_d$ and $T_d$ are highly degenerate, so a fit that expands only in $\beta_d$ can absorb some $T_d$ variance—but not enough at LiteBIRD sensitivity (Vacher et al. 2022).

Operational consequences for component separation:

1. Template methods. Scaling a 353 GHz $(Q,U)$ template by a single MBB (or even by a spatially varying $\beta_d$ map derived from intensity) mis-predicts the dust at 90–150 GHz wherever $\psi(\nu)$ or the effective $\beta_d$ of polarization differs from that of intensity. The residual is a B-mode field that is spatially correlated with the Galaxy and spectrally intermediate between dust and CMB.

2. Pixel-domain parametric fits. A per-pixel MBB has three dust parameters $(A_Q,A_U,\beta_d)$ if $T_d$ is fixed, or four if $T_d$ is free, plus CMB $(Q,U)$. With $\sim 10$ frequency bands this looks overconstrained. It is not: the model is misspecified. McBride, Bull & Hensley (2023) showed that when the true LOS is a distribution of MBBs, the best-fit $(\hat\beta_d,\hat T_d)$ is biased away from both the mean and the median of the true distribution. Hensley & Bull (2018) found that LOS mixing and iron grains are the two physical complexities that most stubbornly bias recovered CMB Stokes parameters, even when a generalized MBB with extra phenomenological parameters is fitted.

3. Harmonic-space likelihoods. Fitting $C_\ell^{BB}(\nu_i\times\nu_j)$ with a single $\beta_d$ and a single decorrelation parameter $\mathcal{R}_\ell$ (Vansyngel et al. 2017) can absorb a mean loss of coherence but cannot represent frequency-dependent $\psi$ or the distinct $\beta$ and $T$ moment spectra. Mangilli et al. (2021) showed that ignoring SED distortions, or replacing them by one $\mathcal{R}_\ell$, biases $r$ above the sensitivity of next-generation experiments. Vacher et al. (2022) found that LiteBIRD-like simulations with spatially varying $(\beta_d,T_d)$ produce a $\sim 21\sigma_r$ bias if a single MBB is assumed; expanding in $\beta_d$ alone is still insufficient; expanding in both $\beta_d$ and $T_d$ recovers an unbiased $r$ with $\sigma_r=8.8\times 10^{-4}$.

4. Blind ILC. Standard ILC nulls a single dust SED. Frequency-dependent mixing leaves a dust residual that ILC interprets as CMB. Moment-constrained ILC (Remazeilles, Rotti & Chluba 2021) deprojects the first few moment templates, at the cost of extra noise and of needing enough frequency channels.

In short: spectral-index decorrelation means that the dust sky at 150 GHz is not a rescaling of the dust sky at 353 GHz, either in amplitude or in polarization angle, and that the residual after a one-component subtraction has a B-mode spectrum that is degenerate with $r$. The rest of this essay is about reconstructing the missing layers rather than fitting their average.

3. How each ingredient helps on its own

3.1 3D maps of thermal dust emission (and the 3D extinction maps that currently exist)

A 2D dust map, however high its S/N, is a column integral. Decorrelation is generated by differences among the emitters that were integrated. A 3D map converts that integral into a sum of discrete layers, each of which can be assigned its own SED and its own $\mathbf{B}$.

What actually exists today are primarily 3D extinction / reddening maps, not 3D maps of submillimetre emission. They are nonetheless the right starting point, because in the diffuse ISM dust mass, optical extinction, far-infrared emission, and H I column are tightly correlated (Boulanger et al. 1996; Lenz, Hensley & Doré 2017). Gaia parallaxes plus stellar photometry and low-resolution BP/RP spectra have produced a sequence of maps with steadily improving distance resolution:

• Green et al. (2019), Bayestar19: Gaia + Pan-STARRS 1 + 2MASS, northern sky, Gaussian-process prior on the dust density.
• Leike, Glatzle & Enßlin (2020) and subsequent information-field-theory reconstructions: parsec-scale structure in the local ISM.
• Edenhofer et al. (2024): 54 million stars with Gaia XP-based distances and extinctions (Zhang, Green & Rix 2023), reconstructed with variational inference (MGVI) on a spherical grid from 69 to 1250 pc in 516 logarithmically spaced distance bins, angular resolution up to $N_{\mathrm{side}}=256$ ($\sim 14'$). This is the first map that resolves the internal structure of hundreds of local molecular clouds at roughly the angular scale of CMB B-mode analyses.
• All-sky Gaia+LAMOST reconstructions (2025) extend the distance reach to several kpc at coarser angular resolution.

At high Galactic latitude the dust that dominates 353 GHz emission lies mostly within a few hundred parsecs—inside the volume of these maps. That is a crucial practical point: B-mode fields are high-latitude fields, so the 3D dust maps are most complete precisely where CMB experiments need them.

What a 3D column map gives a component-separation pipeline, even before grain physics or $\mathbf{B}$ are added:

1. A physically motivated layering. Instead of one MBB per pixel one writes

with $\tau_i$ proportional to the differential extinction $\Delta A_V$ (or $\Delta E_{B-V}$) in distance bin $i$. PySM model d12 already implements a six-layer version of this idea, each layer a different MBB, and is used as the “high-complexity” dust model for CMB-S4 / LiteBIRD / SO forecasts (The Python Sky Model; BICEP/Keck and CMB-S4 pipeline papers). The 3D extinction map replaces those phenomenological layers with data-driven ones.

2. A prior on where LOS complexity exists. Pelgrims et al. (2021) needed H I cloud counts as a proxy for multi-cloud sightlines. A 3D extinction map provides the same information in true distance, including clouds that share a radial velocity (and would be blended in H I) but sit at different distances. Sightlines with a single dominant cloud can be treated with a one-MBB model; sightlines with two or more comparable columns are the ones that must be given extra SED freedom, or masked, or modelled as a sum of MBBs. This is exactly the Tassis & Pavlidou selection, made quantitative.

3. A conversion from extinction to emission, once $T_d$ is assigned. Optical/NIR extinction measures dust mass (weighted by the optical extinction cross-section). Submillimetre emission measures dust mass times $\kappa_\nu B_\nu(T_d)$. Given a 3D radiation-field estimate—or even a simple scaling of $T_d$ with height above the plane and with local $A_V$—one can predict the 353 GHz intensity template per layer. Residuals between that prediction and Planck 353 GHz $I$ then constrain the remaining opacity and temperature variations, rather than having to infer the entire 3D structure from the SED alone.

4. A handle on beam and pixel averaging. Moment expansion treats all averaging as unknown. A 3D map tells the pipeline which moments should be large: pixels whose distance distribution of $\tau_i$ is unimodal should have small $\omega_2^\beta$ and $\omega_1^T$; pixels that pierce two clouds with different local densities (hence different $T_d$) should have large moments, with a predicted frequency dependence.

Limitations of 3D extinction maps as emission templates are taken up in §5. The essential point here is that they convert an underdetermined integral equation into a sum over a modest number of discrete, spatially coherent clouds—the same objects that Tassis & Pavlidou showed are the origin of LOS decorrelation.

3.2 Dust grain size distribution

The grain size distribution $n(a)\,\mathrm{d}a$ and the size-dependent alignment function $f_{\mathrm{align}}(a)$ control three quantities that a single MBB cannot represent: the frequency dependence of $p$, the difference between $\beta_d^{\mathrm{I}}$ and $\beta_d^{\mathrm{P}}$, and the coupling between local physical conditions and the SED of each 3D layer.

Size distribution. The classical MRN distribution $n(a)\propto a^{-3.5}$ from $a_{\min}\sim 5\,\mathrm{nm}$ to $a_{\max}\sim 0.25\,\mu\mathrm{m}$ (Mathis, Rumpl & Nordsieck 1977), and its descendants (Weingartner & Draine 2001; Zubko, Dwek & Arendt 2004; Compiegne et al. 2011; Jones et al. 2017), put most of the dust mass in grains with $a\sim 0.1\,\mu\mathrm{m}$ and most of the number in nanoparticles. Far-infrared and submillimetre emission is dominated by the big-grain population, which is in thermal equilibrium. Very small grains and PAHs are stochastically heated and dominate the mid-infrared; they are at most weakly aligned and contribute negligibly to polarized submillimetre emission. The practical CMB-relevant statement is that polarized dust is a large-grain phenomenon, while total-intensity dust at $\nu\lesssim 100\,\mathrm{GHz}$ can include additional, less polarized, components (AME from spinning nanoparticles; possibly magnetic dipole emission from ferromagnetic inclusions; a carbonaceous population with a different $\beta_d$).

Alignment is a high-pass filter in grain size. In RAT theory the torque efficiency scales as a broken power law in $\lambda/a$ (Lazarian & Hoang 2007). Grains align stably when they rotate suprathermally, $J_{\max}/J_{\mathrm{th}}\gtrsim 3$, which defines a minimum alignment size

(up to order-unity factors and IR damping; Hoang & Lazarian 2014, 2016). In the typical diffuse ISM, $a_{\mathrm{align}}\sim 0.05\,\mu\mathrm{m}$; in denser, less irradiated gas $a_{\mathrm{align}}$ increases and large grains may become the only aligned population. Magnetically enhanced RAT alignment (MRAT) with iron inclusions can drive high-$J$ attractors for a broader range of shapes (Hoang & Lazarian 2016). Radiative torque disruption (Hoang et al. 2019) can set $a_{\max}$ in strong radiation fields. The observational consequence is that $f_{\mathrm{align}}(a)$ is a function of the local radiation field $U$, gas density, and grain magnetic susceptibility. A 3D dust map that also traces density and (crudely) $U$ therefore predicts a 3D alignment map.

Why this matters for decorrelation.

• Frequency-dependent $p$. If two grain populations contribute to $I$ but only one is aligned, $p(\nu)$ is the ratio of a polarized MBB to a sum of MBBs and is not constant. This is the natural interpretation of the Planck Int. XXII decline of $p$ from 353 to 70 GHz, and of the $\beta_d^{\mathrm{P}}-\beta_d^{\mathrm{I}}$ offset: silicates (or astrodust) aligned, a carbonaceous or PAH-related component unaligned (Draine & Fraisse 2009; Guillet et al. 2018). Hensley & Draine (2023) argue instead for a single “astrodust” composition for $a\gtrsim 0.02\,\mu\mathrm{m}$, which predicts a nearly frequency-independent FIR polarization fraction and attributes the optical-to-FIR polarization ratio to grain shape (oblate axial ratio $\sim 1.4$). The two pictures are not settled. For component separation the distinction is operational: a two-population model produces a polarized SED that is not a rescaling of the intensity SED, so using GNILC $\beta_d$ and $T_d$ maps (derived from intensity) to scale polarized templates is a model error. A one-population astrodust model predicts that intensity-derived $(\beta_d,T_d)$ should apply to polarization, and that residual frequency dependence of $p$ comes from LOS mixing of temperatures and field orientations, not from composition.

• Size-dependent temperature. Smaller grains are hotter at fixed radiation field because they absorb optical/UV light more efficiently relative to their FIR opacity. A broad $n(a)$ at fixed $U$ is already a superposition of MBBs, i.e. a source of $\omega_1^T$ and $\omega_2^T$ moments inside a single cloud. McBride et al. (2023) showed that fitting a single MBB to a broad $T_d$ distribution returns biased $(\hat T_d,\hat\beta_d)$. Knowing $n(a)$ lets one compute the intra-cloud moment structure rather than fitting it blindly.

• Iron and magnetic dipole emission. Ferromagnetic inclusions produce microwave magnetic dipole emission whose polarization is orthogonal to the electric-dipole thermal emission of the same grains (Draine & Hensley 2013). Hensley & Bull (2018) identified iron grains, together with LOS mixing, as the two most pernicious biases for parametric CMB recovery. A size- and composition-resolved dust model is the only way to predict whether this component is present at a given frequency and whether it rotates $\psi(\nu)$.

• Connecting 3D extinction to 3D emission. The conversion factor $\tau_{353}/A_V$ depends on the size distribution (optical extinction is dominated by $a\sim\lambda/2\pi\sim 0.1\,\mu\mathrm{m}$ grains; submillimetre opacity has a different size weighting) and on porosity. A constrained $n(a)$ turns a 3D $A_V$ cube into a 3D $\tau_{353}$ cube with a physically bounded systematic error, rather than an unknown scale factor.

In a pipeline, the grain-size ingredient is therefore not a free SED template. It is a mapping

that predicts, for each 3D cell, both the intensity SED and the polarized SED, including their difference.

3.3 Sky-projected magnetic field

The plane-of-sky field $\mathbf{B}_{\mathrm{POS}}$ sets $\psi$ and, through $\gamma$, the local polarization fraction. Without it, even a perfect 3D mass map cannot predict $Q$ and $U$: one would still have to fit two amplitudes per layer. With it, the polarized sky becomes a derived quantity, and frequency dependence of $\psi$ becomes a predicted, rather than fitted, effect.

Observational tracers of $\mathbf{B}_{\mathrm{POS}}$.

1. Submillimetre polarized emission (Planck 353 GHz). This is the most direct full-sky tracer of the dust-weighted, LOS-integrated $\mathbf{B}_{\mathrm{POS}}$. Its limitation is exactly the limitation of 2D dust maps: it is an integral. It is the quantity one wants to explain, not the 3D field itself. It remains the amplitude calibrator of any 3D model.

2. Starlight polarization. Dichroic absorption polarizes background starlight parallel to $\mathbf{B}_{\mathrm{POS}}$. With Gaia distances, a stack of stars along a sightline becomes a tomographic probe of $\psi(d)$ (Panopoulou et al. 2019; Pelgrims et al., subsequent 3D optopolarimetric reconstructions). Tassis & Pavlidou (2015) already proposed this as the correction for LOS decorrelation: identify sightlines with $\Delta\psi\gtrsim 60^\circ$ between distance bins and mask or model them. The bottleneck is the surface density of high-accuracy, low-$p$ optical polarimetry at high latitude (PASIPHAE and related surveys). Where stars exist, this is the only distance-resolved measurement of $\mathbf{B}_{\mathrm{POS}}$ in dust.

3. H I filament orientation. Linear structures in narrow H I velocity channels are aligned with $\mathbf{B}_{\mathrm{POS}}$ (Clark et al. 2014, 2015). The Rolling Hough Transform applied to H I4PI and GALFA-H I yields, with no free parameter relating H I intensity to dust emission, full-sky cubes of Stokes $Q_{\mathrm{HI}}(v),U_{\mathrm{HI}}(v)$ whose velocity integrals correlate with Planck 353 GHz $Q,U$ at $r>0.75$ and reproduce the dust $E/B$ asymmetry at $\ell\lesssim 120$ (Clark & Hensley 2019). H I coherence across velocity is a predictor of $p_{353}$: disordered orientation in $v$–$\theta$ space corresponds to LOS field tangling and depolarization (Clark 2018). This is a 3D (position–position–velocity) magnetic map of the neutral medium that dominates high-latitude dust.

4. Faraday rotation and synchrotron polarization. These trace $\mathbf{B}_{\mathrm{LOS}}$ and the cosmic-ray-weighted field, respectively, i.e. different phases and a different projection. They are complementary constraints on $\gamma$ and on the large-scale Galactic field, not substitutes for the dust-weighted $\mathbf{B}_{\mathrm{POS}}$.

What $\mathbf{B}_{\mathrm{POS}}$ contributes to dedusting, even without 3D dust or grain sizes.

• A polarization-angle template that can be applied layer by layer. Given $\psi_i(\hat n)$ for each cloud, the only free polarized parameter per layer is an amplitude (or $p_{0,i}\cos^2\gamma_i$). Frequency-dependent $\psi_{\mathrm{tot}}$ is then computed from the Tassis–Pavlidou sum, not fitted as a nuisance.

• A data-driven mask or weight for LOS-decorrelated pixels. Pelgrims et al. (2021) demonstrated that an H I-only selection already isolates the pixels in which Planck 353/217 $\Delta\psi$ is anomalous. A B-mode analysis can down-weight those pixels in a parametric fit, or give them extra SED parameters, without using the CMB frequencies to discover the problem after the fact.

• A geometric model of $p$. Because $p_{\mathrm{HI}}$ derived from H I morphology correlates with $p_{353}$ at $r\sim 0.6$ over the full sky (Clark & Hensley 2019), a large fraction of the variance in polarization fraction is geometric (LOS tangling and $\gamma$), not grain-alignment physics. Assigning $\psi_i$ and a geometric depolarization factor per layer therefore absorbs the leading cause of spatially variable $p$, leaving grain physics to explain the residual, including $p(\nu)$.

• A prior on $E/B$ and parity-odd spectra. Filament–field alignment produces $C_\ell^{EE}/C_\ell^{BB}\simeq 2$ and $TE>0$. Misalignment between filaments and $\mathbf{B}$ at the $\sim 2^\circ$ level correlates with the observed dust $TB$ (Huffenberger / Cukierman et al. 2023). A 3D field model that includes a small, spatially coherent misalignment angle can predict, rather than overfit, parity-odd dust spectra that would otherwise leak into cosmic birefringence estimators.

The sky-projected field is the ingredient that turns a 3D intensity model into a 3D polarization model. Without it, extra layers add extra $(Q,U)$ amplitudes and the fit becomes less, not more, constrained.

4. A combined dedusting strategy

The three ingredients are complementary because they constrain different factors in the polarized radiative transfer. A concrete pipeline has a forward-model stage, a parametric (or moment) residual stage, and a validation stage. It is not a slogan about “using more data.”

4.1 Forward model

Discretize each LOS into $N$ layers. In practice $N\sim 2$–$8$ is sufficient at high latitude: Edenhofer-like distance bins can be clustered into magnetically coherent clouds using H I velocity–orientation structure, so that one does not fit hundreds of distance slices independently.

Column and distance (3D dust). For layer $i$ in pixel $\hat n$,

where $\delta A_{V,i}$ is the differential extinction in that distance slice, $K_V$ is the extinction per H nucleon (set by $n(a)$ and composition), and $\kappa_{353}$ is the 353 GHz mass opacity (also set by $n(a)$, porosity, and composition). An overall scale $\kappa_{353}/K_V$ can be calibrated on the high-latitude $\tau_{353}$–$A_V$ correlation and on the $N_{\mathrm{HI}}$–$E_{B-V}$ relation of Lenz et al. (2017), which holds for $N_{\mathrm{HI}}<4\times 10^{20}\,\mathrm{cm}^{-2}$.

Local physical state (grain sizes + environment). Associate to each layer a radiation-field strength $U_i$ (from the local $A_V$-attenuated ISRF, or from a $T_d$ prior) and a density $n_{\mathrm{H},i}$ (from $\delta A_{V,i}$ and a fiducial cloud thickness, or from H I brightness temperature in the matching velocity interval). Then:

• Compute $T_d(a)$ for the local $U_i$ by equating absorption of the ISRF to FIR emission, using the adopted dielectric function (astrodust or silicate+carbon).
• Compute $a_{\mathrm{align},i}$ from RAT/MRAT given $(U_i,n_{\mathrm{H},i},T_{\mathrm{gas},i})$.
• Integrate over $n(a)f_{\mathrm{align}}(a)$ to obtain layer-averaged $\bar T_{d,i}$, $\bar\beta_{d,i}$, $\bar p_{0,i}$, and intra-layer moments $\omega^{\beta,T}_{k,i}$ if the size distribution is broad.

This step is the difference between a phenomenological multi-MBB and a physical multi-MBB. It predicts, rather than fits, $\beta_d^{\mathrm{P}}-\beta_d^{\mathrm{I}}$ and $p(\nu)$ within a layer.

Orientation (sky-projected $\mathbf{B}$). Assign $\psi_i(\hat n)$ from the H I-derived angle in the velocity range matched to layer $i$, or from starlight polarization of stars whose Gaia distances fall in the layer, or from a combination (starlight where available, H I elsewhere). The geometric factor $\cos^2\gamma_i$ is not measured directly. Two practical options, both used in the literature, are: (i) absorb $\cos^2\gamma_i$ into a single amplitude $A_i^P=p_{0,i}\cos^2\gamma_i\,\tau_i$, fitted to 353 GHz $Q,U$; (ii) use H I coherence / RHT peak width as a proxy for LOS tangling and $\gamma$, as in Clark (2018), which already predicts a large fraction of $p_{353}$. Faraday $\mathbf{B}_{\mathrm{LOS}}$ can provide an independent prior on $\gamma$ in regions where the ionized and neutral fields are correlated, but this is uncertain and should be treated as optional.

Emission integral. The model Stokes vector at each CMB frequency is

If the aligned population is a subset of the emitting population, $\beta_{d,i}^{\mathrm{P}}\neq\beta_{d,i}$ and the polarized amplitude does not scale with $I_{\nu,i}$. That is a prediction of the grain model, not an extra free map.

This forward model automatically generates:

• frequency-dependent $\psi_{\mathrm{tot}}(\hat n,\nu)$ on multi-cloud, misaligned sightlines (Tassis & Pavlidou);
• $\mathcal{R}_\ell^{BB}(\nu_1,\nu_2)<1$ with a specific frequency dependence set by the layer $\Delta T_d$ and $\Delta\beta_d$, not by a free $\mathcal{R}_\ell$;
• a difference between intensity-derived and polarization-derived effective spectral indices;
• geometric $E/B$ asymmetry from filament–field alignment.

4.2 Parametric residual fit

The forward model will not be perfect. The residual is fitted to the multi-frequency CMB data with a restricted parameter set, using the 3D model as a strongly informative prior rather than as a rigid template.

A workable hierarchy, from fewest to most free parameters:

H0. Template cleaning with a 3D-informed mask. Use Planck 353 GHz $(Q,U)$ as the dust template, scale with a global or piecewise $\beta_d$, but exclude or deweight pixels that the 3D+H I selection flags as LOS-decorrelated (Pelgrims target sample). This is the minimal use of the three ingredients and is already an improvement on a blind 353 GHz extrapolation. Tassis & Pavlidou proposed the starlight-polarization version of this mask in 2015.

H1. Layered MBB with amplitudes fitted at a reference frequency. Fix $\{\psi_i,\beta_{d,i},T_{d,i}\}$ from the forward model. Fit $\{A_i^P\}$ (and optionally a small $\delta\beta_{d,i}$) to the ensemble of frequency maps. The CMB is a blackbody added in $(Q,U)$. This is a physically layered version of Commander/FGBuster: the number of dust amplitudes per pixel is $N_{\mathrm{lay}}$, but $N_{\mathrm{lay}}$ is 1 on most high-latitude pixels and 2–3 only where the 3D map demands it. Degeneracies are controlled because $\psi_i$ is not free.

H2. Layered MBB plus intra-layer moments. Allow $\omega_1^{\beta,i}$ and $\omega_1^{T,i}$ per layer, or, cheaper, a single set of moments on the sum of layers to capture beam averaging and unresolved substructure (Chluba et al. 2017; Mangilli et al. 2021; Vacher et al. 2022). At LiteBIRD sensitivity the Vacher et al. result implies that moments in both $\beta_d$ and $T_d$ are required. The 3D model predicts the expected size of those moments, which can be imposed as a Gaussian prior, preventing the moments from absorbing CMB.

H3. Harmonic-space hybrid. For $r$ likelihoods of the BICEP/Keck type, replace the single-SED, single-$\mathcal{R}_\ell$ dust model by the cross-spectra of the layered forward model, plus a low-order moment expansion of the residual. The layered model supplies the frequency dependence of $\mathcal{R}_\ell^{BB}$ and of $C_\ell^{BB}(\nu_i\times\nu_j)/C_\ell^{EE}(\nu_i\times\nu_j)$; the moments absorb isotropic SED distortions from averaging. This is the natural interface to LiteBIRD and CMB-S4 power-spectrum pipelines.

Synchrotron is treated in parallel: a power law with spatially varying $\beta_s$, optionally a curvature term, and an independent $\psi_s$ (the synchrotron-emitting volume is not the dust-emitting volume). H I-based $\psi_i$ should not be imposed on synchrotron. Joint dust–synchrotron modelling is needed below $\sim 100\,\mathrm{GHz}$ because a dust residual with a wrong $\beta_d$ is degenerate with synchrotron, and conversely.

4.3 Why the combination is more than the sum of parts

Each ingredient removes a specific degeneracy:

A 3D mass map without $\mathbf{B}$ overfits $Q,U$. A $\mathbf{B}$ map without 3D mass cannot assign SED differences to layers. Grain physics without either map can predict a mean $p(\nu)$ but not its spatial pattern. Together they specify the Tassis–Pavlidou two-cloud (and $N$-cloud) model on the real sky.

4.4 Validation, not only subtraction

A physical forward model is falsifiable, which a phenomenological MBB is not. Required null tests, all of which are already close to existing analyses:

• On the Pelgrims control sample (single H I cloud), the layered model must collapse to a single MBB and $\psi(\nu)$ must be constant; any remaining 353/217 $\Delta\psi$ is then a bound on residual systematics and on intra-cloud grain physics.
• On the target sample, the model must predict the sign and the approximate magnitude of $\Delta\psi(353,217)$ from the H I-measured $\alpha$ and the 3D-inferred $\Delta T_d$, before looking at 217 GHz. Pelgrims et al. estimated a 15% polarized-SED variation; a grain+3D model should predict that number from $T_d$ and $\beta_d$ differences, not fit it.
• The predicted $p_{353}$ from H I coherence plus RAT $p_0\cos^2\gamma$ should match Planck $p_{353}$ with residuals that correlate with tracers of grain alignment (e.g. $A_V$, radiation field), not with the CMB.
• Cross-spectra $\mathcal{R}_\ell^{BB}(\nu_i,\nu_j)$ of the model, evaluated on the actual experiment bandpasses, should be compared to Planck 217$\times$353 and to BICEP/Keck 95$\times$150$\times$220$\times$Planck 353 without introducing a free $\mathcal{R}_\ell$. A free $\mathcal{R}_\ell$ is reserved for the residual.
• Split-frequency ILC or spectral-imaging tests (e.g. the bolometric-interferometry argument of Maffei et al. 2024) should return a consistent $r$ on the cleaned maps. The point of those tests is to detect unmodelled LOS decorrelation; a successful 3D model should make the detection vanish.

5. Limitations, missing data, and residual systematics

A layered physical model reduces, but does not eliminate, foreground bias on $r$. The leading limitations are as follows.

3D maps trace extinction, not submillimetre emission. The conversion $\delta A_V\to\tau_{353}$ assumes a dust-to-extinction ratio and a size distribution that may vary with environment ($R_V$ variations, grain growth in molecular gas, carbon depletion). High-latitude sightlines are mostly diffuse, where Lenz et al. (2017) find $\sim 10\%$ scatter in $N_{\mathrm{HI}}/E_{B-V}$, but 10% layer-to-layer opacity errors, if correlated with $T_d$ or $\psi$, are enough to leave a residual at the $r\sim 10^{-3}$ level. CIB anisotropies further contaminate intensity-based $\tau_{353}$ and $T_d$ maps at high latitude; GNILC separation helps but is not perfect. Polarized CIB is expected to be small but is not measured at the required level.

Distance completeness. Gaia-based maps reach $\sim 1.25\,\mathrm{kpc}$ at high resolution (Edenhofer et al. 2024) and several kpc at lower resolution. That covers most high-latitude dust, but Intermediate- and High-Velocity Clouds (IVCs/HVCs) can contribute polarized emission at high $|b|$ and lie outside the local volume or are poorly sampled by stars. Pelgrims et al. used IVC/LVC H I components precisely because they are common at high latitude. A pipeline that trusts only Gaia distances will mis-model those layers unless H I velocity components are included as additional, distance-free clouds.

H I velocity is not distance. The Clark & Hensley cubes are in $v_{\mathrm{LSR}}$, not parsecs. At high latitude Galactic rotation does not provide a unique kinematic distance; two clouds at the same velocity can sit at different distances, and a single cloud can be split across velocities by turbulence. Matching H I components to 3D extinction slices is a cross-identification problem, not a solved mapping. Errors in that matching assign the wrong $\psi$ to the wrong SED, which can create a spurious $\psi(\nu)$ as easily as it corrects a real one. Starlight polarization at known distances is the only currently available external check.

The inclination $\gamma$ is poorly measured. Polarized emission constrains $\mathbf{B}_{\mathrm{POS}}$, not $\mathbf{B}_{\mathrm{LOS}}$. The factor $\cos^2\gamma$ is degenerate with $p_0$ and with grain alignment efficiency. H I filament detectability encodes some $\gamma$ information (filaments end-on are less linear on the sky), and starlight $p/A_V$ encodes some, but neither is a precision $\gamma$ map. Wrong $\gamma_i$ mis-weight the layers in $Q,U$ and therefore bias the predicted $\psi(\nu)$.

Grain model uncertainty is a systematic on $p(\nu)$ and $\beta_d^{\mathrm{P}}-\beta_d^{\mathrm{I}}$. Silicate+carbon with unaligned carbons (Draine & Fraisse 2009; Guillet et al. 2018) versus astrodust (Hensley & Draine 2021, 2023) versus models with magnetic dipole emission (Draine & Hensley 2013) make qualitatively different predictions for the polarized SED below $\sim 100\,\mathrm{GHz}$. Hensley & Bull (2018) showed that iron grains are among the worst parametric biases. Until FIR polarimetry between 100 and $850\,\mu\mathrm{m}$ (the diagnostic proposed by Hensley & Draine) distinguishes these models, a B-mode pipeline must marginalize over at least two grain scenarios, not adopt one.

Alignment physics is not a solved 3D field. RAT predictions for $a_{\mathrm{align}}(U,n_{\mathrm{H}})$ are tested statistically, not cloud-by-cloud at high latitude. Dense, shielded gas can have inefficient alignment; locally enhanced radiation fields can raise $T_d$ and $a_{\mathrm{align}}$ together. Mis-modelling alignment produces a spatially variable $p_0$ that the pipeline may try to absorb into $\gamma$ or into extra layers.

Unresolved sub-beam structure still generates moments. Even a perfect cloud catalogue does not resolve MHD turbulence inside a CMB beam (LiteBIRD’s 140 GHz beam is $\sim 31'$; CMB-S4 will be finer but not arbitrarily fine). Intra-beam $\mathbf{B}$ tangling and $T_d$ variations remain, and they are exactly what the Chluba moment expansion is for. The 3D model reduces the moments; it does not set them to zero. Vacher et al. (2022) showed that at LiteBIRD noise, omitting $T_d$ moments biases $r$ even when $\beta_d$ moments are included. That conclusion survives a 3D prior.

Planck 353 GHz polarization S/N is limited at high latitude. All 3D polarized models are calibrated on Planck 353 GHz $Q,U$. In the faintest B-mode fields those maps are noisy at the native resolution; GNILC 80′ products are used instead, which already average over the small-scale structure that generates small-scale decorrelation (Planck Int. L emphasized that decorrelation grows toward $\ell>50$). Extrapolating a smoothed 353 GHz template cannot recover that small-scale $\psi(\nu)$. Higher-resolution polarized dust data (e.g. from a LiteBIRD high-frequency channel itself, BLAST-like balloon data, or future FIR polarimeters) are needed to calibrate layers at CMB-S4 resolution.

Instrumental and analysis systematics mimic LOS decorrelation. Bandpass mismatch, $I\to P$ leakage, gain mismatch, and residual CMB in the 353 GHz template all produce frequency-dependent $(Q,U)$ residuals. Pelgrims et al. found that correcting Planck residual systematics increased the significance of LOS decorrelation, which is encouraging (the astrophysical effect is not an artefact), but the converse is also possible in other data sets. Spider found no evidence for LOS decorrelation in its patch (SPIDER Collaboration 2025), consistent with a region that is magnetically simpler, but also a reminder that the Pelgrims detection is not a uniform sky property. A pipeline that applies a 3D correction derived from Planck to a ground-based deep field must re-estimate the correction on that field’s own ancillary data.

Other foregrounds. Polarized synchrotron with spatially varying $\beta_s$ and possible curvature; a possible polarized AME at the percent level; CO and other line polarization at a poorly known level; polarized compact sources. A perfect dust model does not yield an unbiased $r$ if synchrotron is wrong. The 3D dust/B-field construction does not automatically help synchrotron, which requires low-frequency surveys (C-BASS, S-PASS, QUIJOTE, Planck LFI, future 5–20 GHz polarimetry).

Statistical cost. Extra layers and moments consume degrees of freedom. Vacher et al. (2022) already saw $\sigma_r$ increase because dust moments correlate with the CMB SED. A 3D model that adds unconstrained amplitudes per layer will degrade $\sigma_r$ more than it reduces bias. The whole point of the three ingredients is to replace free maps by constrained ones: $\psi_i$ from H I/stars, $\tau_i$ from Gaia, $\beta_{d,i}$ and $T_{d,i}$ from grain physics with tight priors. If those constraints are too tight and wrong, one reintroduces bias; if they are too loose, one reintroduces degeneracy. That bias–variance trade is the central design problem, and it is not solved by any current end-to-end demonstration on the real sky at $\delta r\sim 10^{-3}$.

What is missing as data, not as theory.

• Dense, high-latitude optical polarimetry of stars with Gaia distances (PASIPHAE-like), sufficient to measure $\Delta\psi$ between the two or three local dust layers in every B-mode pixel.
• A joint 3D reconstruction that uses Gaia extinction, H I cubes, and starlight polarization in one likelihood, rather than matching products ex post.
• FIR polarimetry from $\sim 100$ to $850\,\mu\mathrm{m}$ to decide between grain models and to measure $T_d$ in polarization, where CIB contamination is weaker (Hensley & Draine).
• Direct 3D emission maps—i.e. tomography of the FIR intensity itself—which do not yet exist; extinction is a proxy.
• A well-characterized $\gamma$ tracer in the neutral medium.
• End-to-end simulations that start from an MHD+dust cube with a realistic $n(a)$ and RAT alignment, observe it with LiteBIRD/CMB-S4 bandpasses, and show that the pipeline in §4 recovers $r$ without a bias larger than the statistical error when the reconstruction uses only the ancillary data that will actually exist. Present “high-complexity” PySM suites (layered MBBs, d12) capture LOS decorrelation phenomenologically but do not yet ingest Gaia 3D maps or RAT-based alignment.

Until those exist, the honest status is: the three ingredients specify a well-posed physical model for the exact mechanism (Tassis & Pavlidou 2015; Pelgrims et al. 2021) that makes a single-MBB-per-pixel subtraction fail; they can already be used to flag and down-weight the worst sightlines and to predict the frequency dependence of $\mathcal{R}_\ell$ and $\psi(\nu)$ at a level that phenomenological decorrelation parameters cannot; they are not yet a turn-key replacement for moment expansions and multi-frequency residual fits. The right operational posture for LiteBIRD and CMB-S4 is to implement the layered forward model as an informative prior and as a diagnostic, to keep a low-order $\beta_d$–$T_d$ moment expansion for unresolved averaging, and to treat grain-composition and $\gamma$ uncertainties as discrete systematics that are marginalized, not assumed away.

6. Summary

Spectral-index decorrelation of polarized dust is the statement that the effective MBB parameters and the polarization angle of Galactic dust are frequency-dependent, because Stokes $I,Q,U$ are integrals of emitters that do not share a common $(T_d,\beta_d,\psi)$. A single modified blackbody per pixel—or a 353 GHz template scaled by a global $\beta_d$—is then a misspecified model, and the residual B-modes are degenerate with $r$ at the $10^{-3}$ level relevant to LiteBIRD and CMB-S4.

Three-dimensional dust maps (presently Gaia extinction cubes, matched to H I clouds) replace the integral by a sum over a few distance layers and identify the sightlines where the Tassis–Pavlidou mechanism operates. The grain size distribution and RAT alignment determine which grains emit in polarization, why $\beta_d^{\mathrm{P}}\neq\beta_d^{\mathrm{I}}$ and $p=p(\nu)$, and how to convert extinction into a polarized SED given the local radiation field and density. The sky-projected magnetic field, from H I filaments, starlight tomography, and Planck 353 GHz, assigns $\psi_i$ (and a geometric proxy for depolarization) to each layer, so that frequency-dependent $\psi_{\mathrm{tot}}$ is predicted rather than fitted.

The combined pipeline is a layered, physically parameterized forward model whose remaining freedom is a small set of amplitudes and intra-layer moments, fitted to the CMB frequency maps under priors taken from the ancillary 3D data. Its limitations are real: extinction is not emission, velocity is not distance, $\gamma$ and grain composition are underconstrained, and beam averaging still requires moments. Used as a prior and a diagnostic rather than as an exact sky, the same construction is currently the most direct way to attack the specific astrophysical failure mode that spectral-index decorrelation names.