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

Three-dimensional, grain-informed modelling of polarized dust for CMB $B$-mode cleaning

Here “spectral-index decorrelation” should be understood as frequency decorrelation of polarized Galactic dust: the Stokes fields $Q$ and $U$ at one frequency are not everywhere a single scalar multiple of those at another because dust temperature $T_{\rm d}$, opacity index $\beta_{\rm d}$, polarization efficiency, and polarization angle vary both across the sky and along each line of sight (LOS). Operationally, for dust $B$-modes one measures

often denoted $\Delta_{\rm d}$ in BICEP/Keck analyses. Perfect morphological correlation gives $\mathcal R_\ell^{BB}=1$. Importantly, broadband decorrelation at CMB frequencies is physically plausible but not yet established as a general property of clean high-latitude fields: Planck 2018, BK18, BICEP/Keck’s H I analysis, and SPIDER found no statistically significant decorrelation in their respective principal tests, although targeted Planck sightlines do show LOS-dependent polarization-angle rotation and later Planck analyses find polarized-SED residuals inconsistent with the simplest frequency scaling (Planck 2018 XI; BK18; BICEP/Keck XVI; SPIDER 2025; Pelgrims et al. 2021).

1. Polarized thermal dust and the $B$-mode problem

Non-spherical interstellar grains absorb and emit differently along their principal axes. When their short axes preferentially align with the Galactic magnetic field, their thermal emission is linearly polarized perpendicular to the plane-of-sky magnetic field, $\boldsymbol B_{\rm POS}$. Spatial variations in dust column and magnetic-field geometry generate both parity-even $E$-modes and parity-odd $B$-modes. Galactic dust has approximately twice as much $EE$ as $BB$ power and positive $TE$, largely because elongated neutral-gas structures tend to align with $\boldsymbol B_{\rm POS}$. Misalignment between structures and the field can also generate $TB$ and $EB$.

Dust is especially troublesome near the recombination peak, $\ell\simeq80$, where inflationary tensors would produce their most accessible $B$-mode signal. Current practice models each pixel or patch by a modified blackbody (MBB),

and extrapolates a high-frequency polarized template—usually Planck 353 GHz or an internally reconstructed dust map—to 90–150 GHz. This works if the same morphology is multiplied by one frequency law. It fails if different parts of a pixel or LOS have different SEDs and polarization orientations.

The required accuracy is severe. Using the usual

take $D_{80}^{BB,{\rm d}}(353)\simeq4.5\,\mu{\rm K}^2$ in the BICEP/Keck field. For $\beta_{\rm d}=1.6$, $T_{\rm d}=19.6$ K, and thermodynamic CMB units, the amplitude conversion is

where $G_\nu=\partial B_\nu(T_{\rm CMB})/\partial T_{\rm CMB}$. Hence

A representative tensor signal at $r=10^{-3}$ is only $D_{80}^{BB,{\rm tens}}\simeq7\times10^{-5}\,\mu{\rm K}^2$. For optimally rescaled template subtraction, decorrelation leaves approximately

An unmodelled $1-\mathcal R\simeq3.8\times10^{-3}$ therefore produces residual power comparable to an $r=10^{-3}$ tensor signal; reducing it to a 10% correction requires control at a few $10^{-4}$. Planck’s null result was adequate for $r\sim0.01$, but it does not establish the accuracy needed at $10^{-3}$.

2. Polarized emission is a vector LOS integral

For an optically thin LOS, define

Following the standard partially aligned-grain formulation of Lee & Draine, Wardle & Königl, and Planck Intermediate XX, the emitted Stokes parameters are

Here $\gamma$ is the angle between the magnetic field and the sky plane, while $\psi$ is the thermal-emission polarization angle: it is rotated by $90^\circ$ from $\boldsymbol B_{\rm POS}$. Starlight polarization is instead parallel to $\boldsymbol B_{\rm POS}$. The factor $p_0$ combines grain shape, polarized opacity, aligned mass fraction, and the Rayleigh reduction factor.

These equations use the IAU convention, $P\equiv Q+iU=|P|e^{2i\psi}$, with position angle measured north through east. Some HEALPix/CMB conventions reverse the sign of $U$; all external angle data must be placed in one convention before fitting.

The equations expose three distinct ways that geometry suppresses or changes polarization:

1. $\cos^2\gamma$ suppresses a cloud whose field points toward the observer.
2. Variations of $\psi$ cause vector cancellation between clouds or within the beam.
3. If the canceling components have different SEDs, the amount and direction of the cancellation change with frequency.

Cloud number alone is therefore not a decorrelation predictor. Multiple clouds whose fields and SEDs are identical remain perfectly correlated.

Superposition is not another MBB

Even for intensity, adding MBBs with different $T$ or $\beta$ does not give an exact MBB. A useful local representation is the moment expansion (Chluba, Hill & Abitbol 2017),

where the $\omega_n$ are amplitude-weighted moments of the spectral-parameter distribution; temperature mixing introduces the corresponding derivatives of $B_\nu(T)$.

For polarization, the moments are spin-2, and hence generally complex (Vacher et al. 2023):

Because $P_i$ carries the phase $e^{2i\psi_i}$, $W_n$ contains both spectral and geometrical information. At first order,

Thus a real moment changes polarized amplitude, whereas an imaginary moment rotates the polarization angle. This is the natural bridge between 3D cloud information and existing multifrequency component-separation pipelines.

Temperature mixing is relevant even below 353 GHz because the dust is not fully in the Rayleigh–Jeans limit. The logarithmic spectral slope is

At 150 GHz the slope is 1.741 for $T=15$ K and 1.849 for $T=25$ K: that temperature contrast changes the local slope by $0.108$, equivalent to $\Delta\beta\simeq0.1$. Around 20 K, a 1 K temperature change alters the 150/353-GHz colour by an amount equivalent to $\Delta\beta\simeq0.018$, conventionally rounded to 0.02. A $\beta$-only fit can therefore absorb temperature mixing and return a biased “spectral index” (Kogut & Fixsen 2016).

Worked two-cloud example

For two clouds write

Their combined angle obeys

and

Let $\beta_1=1.45$, $\beta_2=1.65$, $T_1=T_2=20$ K, $r(353)=0.5$, and $\Delta\psi=45^\circ$. The temperature term vanishes, so

Because $2\Delta\psi=90^\circ$,

The polarization rotates by $1.86^\circ$. The angular part of the normalized polarization-vector correlation loses

In addition, $|1+r e^{2i\Delta\psi}|$ changes with frequency, so the polarized amplitude has an SED different from the scalar intensity SED. Consequently, $\beta_P\ne\beta_I$ does not by itself prove that intensity and polarization arise from different grain species. Geometry alone can produce the difference. This is the LOS extrapolation effect emphasized by Tassis & Pavlidou (2015) and Poh & Dodelson (2017).

3. What decorrelation means statistically

Consider the simpler case of spatially varying $\beta=\bar\beta+\delta\beta$, with Gaussian variance $\sigma_\beta^2$, independent of dust amplitude. At frequency $\nu$, the fluctuating scale factor is $\exp[\delta\beta L]$. Gaussian averaging gives

This idealized result motivates the BICEP/Keck/Vansyngel ansatz

It is not a physical model for non-Gaussian, amplitude-correlated, or angle-rotating dust. It is a compact likelihood parameterization.

For frequency-dependent angle rotation with otherwise unchanged amplitude,

A power-spectrum correlation coefficient mixes spatial SED changes, LOS rotations, beam averaging, masking, noise, and calibration. A measured $\mathcal R<1$ is therefore not uniquely interpretable as variance of a scalar $\beta$.

The observational status reflects this difficulty. Planck Intermediate L initially reported 217×353-GHz decorrelation, but the result was sensitive to noise bias, correlated sky masks, and instrumental systematics (Sheehy & Slosar 2018); Planck 2018 found no significant evidence. BK18 likewise found no preference for a decorrelated dust model, and SPIDER found no significant LOS decorrelation, although one subregion disfavoured a common-morphology MBB over its full frequency range. In contrast, Pelgrims et al. selected sightlines using independent H I structure and detected frequency-dependent 217×353-GHz angle differences. A roughly 15% cloud-to-cloud spread in the polarized 353/217 ratio reproduced the observations. Ritacco et al. found low-$\ell$ polarized-SED residuals and evidence that angle variation must be represented, but Planck systematics remain limiting (Ritacco et al. 2023; Guillet et al. 2025).

4. What the three external ingredients contribute

4.1 Three-dimensional dust maps

Present 3D maps measure extinction or reddening toward stars, using stellar photometry/spectroscopy and Gaia parallaxes. They do not directly tomographically measure thermal emission. For example, the Edenhofer et al. map reaches 1.25 kpc, with angular resolution as fine as $14'$, parsec-scale radial sampling, and released posterior samples (Edenhofer et al. 2024). Such maps identify whether a LOS contains one dominant layer or several separated extinction peaks, estimate their distances and relative columns, and delineate clouds coherently across adjacent pixels.

Three-dimensional temperature products are an additional forward inference. Zelko et al. placed $T$, $\beta$, and an extinction-to-emission conversion in 3D reddening voxels, then required their LOS sums to fit Planck and IRAS intensity maps. Their product covers about three quarters of the sky at $27'$–$110'$ temperature resolution and finds up to factor-two variation in the conversion between reddening and emission optical depth (Zelko et al. 2022). It can distinguish temperature differences between strong distance components, but $T$, $\beta$, and opacity remain model-dependent and degenerate.

The gain for CMB cleaning is therefore not an optical prediction of the 150-GHz amplitude. It is a physically motivated basis: how many components should be fitted, where they lie, which pixels share a cloud, and where temperature mixing should produce non-zero spectral moments.

4.2 Grain size distribution and composition

Large equilibrium grains, approximately $a\gtrsim0.01\,\mu{\rm m}$, dominate far-infrared and millimetre emission. In a simplified heating balance,

up to weak size- and composition-dependent factors. Here $U$ is the interstellar radiation-field intensity. The size distribution therefore controls the temperature distribution and hence temperature-mixing moments; $U$ additionally depends on 3D environment.

Grain size also controls alignment. Radiative-torque models predict efficient alignment only above a threshold $a_{\rm align}$, commonly around $0.05$–$0.1\,\mu{\rm m}$ in diffuse gas and increasing approximately as $(n_{\rm H}T_{\rm gas}/U)^{2/7}$, with additional dependence on radiation anisotropy and magnetic properties. The Astrodust+PAH fit has negligible alignment for the smallest grains and nearly full alignment above about $0.1\,\mu{\rm m}$ (Hensley & Draine 2023). Polarization consequently weights larger, often cooler grains differently from total intensity.

By contrast, in the Rayleigh limit the millimetre opacity index is controlled primarily by dielectric material and structure, not grain radius. Silicate/carbon mixtures, composite “astrodust,” amorphous two-level-system absorption, ice mantles, coagulation, and accretion can alter $\beta$ and the opacity. Planck found $\beta_{\rm mm}\simeq1.53$ in atomic-dominated Galactic-plane regions and $\simeq1.65$ in molecular-dominated regions (Planck Intermediate XIV). Earlier diffuse-sky measurements gave $\beta_P\simeq1.59$ versus $\beta_I\simeq1.51$, while Planck 2018 obtained $\beta_P=1.53\pm0.02$ over its masks (Planck Intermediate XXII).

A single-composition astrodust model predicts a nearly flat far-infrared-to-millimetre polarization fraction, whereas a model with aligned silicates and relatively unaligned carbonaceous grains can produce a rising polarization fraction toward long wavelengths. BLASTPol observations of Vela C found a broadly flat submillimetre polarization spectrum after accounting for LOS temperature structure. Flatness supports template extrapolation; cloud-dependent departures provide another decorrelation channel.

Optical extinction constrains, but does not uniquely determine, this physics. $R(V)=A_V/E(B-V)$ responds to both size and composition. Gaia XP spectra now provide extinction-curve measurements toward about 130 million stars and a 3D $R(V)$ reconstruction; higher-order extinction-curve components contain information not captured by $R(V)$ alone (Zhang & Green 2025; Green et al. 2024). Grain models predict strong $R(V)$-dependent variation in $A_V/\tau_{353}$, but reproducing the observed $R(V)$–$\beta$ behaviour requires composition as well as size changes. Optical data should therefore provide priors on $\tau/A_V$, $T$, $p_0$, and $\beta$—never a deterministic $\beta(R_V)$ law.

4.3 The sky-projected magnetic field

Starlight polarimetry plus Gaia distances is the cleanest route to per-cloud $\psi_i$. Starlight Stokes parameters are cumulative; differencing their statistically inferred values across distance layers isolates the polarization induced by an individual cloud. Demonstrations now recover several clouds and changes in $\boldsymbol B_{\rm POS}$ out to approximately 2 kpc, but degree-scale tomography has so far been produced only in small fields (Pelgrims et al. 2023; Pelgrims et al. 2024). Simulations indicate that a cloud generally needs to induce $\gtrsim0.1\%$ starlight polarization for reliable recovery.

Where stars are sparse, H I velocity-channel morphology provides an orientation proxy because cold neutral filaments tend to align with $\boldsymbol B_{\rm POS}$. H I separates components kinematically, not unambiguously by distance, although the 2026 Milky Way Atlas combines H I velocities with a 3D dust prior to reconstruct local gas within 1.25 kpc (McCallum et al. 2026). Local-Bubble shell geometry is especially useful at $|b|\gtrsim60^\circ$: models place much of the relevant dust at roughly 100–300 pc and find that the shell likely dominates polarization over many high-latitude sightlines (O’Neill et al. 2024).

The field inclination $\gamma_i$ is not directly recovered by dust or starlight polarization. It can only be constrained statistically through polarization efficiency, alignment models, and geometrical assumptions such as a field tangent to a shell. Planck’s maximum $p\simeq22\%$ and the anticorrelation of $p$ with angle dispersion strongly constrain such models (Planck 2018 XII). Multiple distance components reduce $p$ at fixed column through cancellation, as directly found by Halal, Clark & Tahani (2024).

5. A concrete dedusting strategy

I advocate a hierarchical cloud-basis polarized-SED fit, with 3D data used to determine the basis and its priors while the millimetre data determine the amplitudes. Spin moments are derived from that fit and used both for compression and for robustness tests.

1. Segment and propagate map uncertainty. For every CMB pixel, decompose posterior samples of the 3D extinction profile into distance components using peaks or a spatially coherent dendrogram. Components must be associated across pixels so that a physical cloud, rather than every pixel, is the fundamental object. Retain alternative decompositions instead of selecting one cloud count irrevocably.

2. Build a grain-and-field prior layer. Assign each component a probabilistic conversion from extinction to $\tau_{353}$, temperature, $\beta$, polarization spectrum, and intrinsic polarization efficiency. Condition these priors on $R(V)$, extinction-curve components, density, radiation environment, and cloud class. Attach a per-cloud $\psi_i$ posterior from starlight tomography, supplemented by H I or a Local-Bubble model. Treat $\gamma_i$ and sub-beam tangling as nuisance parameters.

3. Fit the multifrequency polarization in a cloud basis. In thermodynamic units, use

The known or uncertain mixing weights $M_{pi}$ encode the 3D extinction morphology. A cloud SED $f_i(\nu)$ is shared across the pixels spanned by that cloud, with controlled spatial scatter. Crucially, $A_i$ remains a millimetre-wave fit parameter: optical extinction cannot supply absolute polarized amplitudes at percent precision. This construction has fewer effective degrees of freedom than fitting unrelated $(A,\beta,T,\psi)$ in every pixel, breaks part of the $\beta$–$T$–noise degeneracy, and explicitly represents LOS cancellation.

The posterior should provide a CMB map, cloud dust maps, predicted $\delta\psi(\nu)$, local frequency-correlation maps, and $W_1,W_2$. Cloud decompositions, extinction-to-emission conversion, and magnetic angles are marginalized, not inserted as exact templates.

4. Use the 3D model in four distinct ways.

• Primary separation: use the hierarchical cloud-basis likelihood above to infer CMB and foregrounds jointly.
• Moment control: translate the cloud posterior into priors on spin moments. Deproject only the moments predicted to be significant in each region using a constrained-moment ILC or fit them parametrically. This limits unnecessary noise amplification.
• Survey weighting: rank patches by predicted LOS complexity—such as comparable-strength misaligned clouds, large temperature contrast, or large ${\rm Im}(W_1)$—and optimize masks and observing depth independently of total dust intensity.
• Bias simulations: generate 3D-consistent skies in which dust column, SED, alignment, and magnetic geometry covary. These replace purely ad hoc emitting layers and complement maximally complex moment skies and MHD simulations such as Kritsuk et al. (2026).

5. Validate without reusing the same information circularly. Hold out frequency channels and clouds, then test predicted $Q/U$, 217×353 angle differences, $EE/BB$, $TE$, $TB$, $EB$, and the scale dependence of $\mathcal R_\ell$. Test first on H I- or 3D-selected complex sightlines, where LOS effects should be enhanced, and separately on simple sightlines. Extinction and stellar-angle templates contain no CMB, CIB, or millimetre-instrument noise, so their cross-correlations with each microwave band provide clean estimates of dust-correlated amplitude and colour, much as H I templates do in BICEP/Keck XVI. The fitted Planck/IRAS-derived temperature map is not independent in this sense and must not be used as its own validation.

Near-term leverage comes from PASIPHAE, whose design target is $\sim0.1\%$ polarization accuracy over more than $10^4$ square degrees of the polar caps, and from future Gaia releases. At millimetre wavelengths, Simons Observatory’s 225/280-GHz channels and LiteBIRD’s planned bands extending to approximately 400 GHz can distinguish $T$, $\beta$, and angle moments far better than 90/150/353 GHz alone (PASIPHAE; Simons Observatory; LiteBIRD).

6. Limitations and residual systematics

First, extinction maps are radially smoothed reconstructions from an inhomogeneous stellar sample. Angular and distance resolution degrade with distance and extinction; nearby clouds may be split artificially, while distant or diffuse components may be missed. A 1.25-kpc map omits high-latitude intermediate-velocity clouds known to extend toward 2 kpc. H I velocity components are not automatically distinct distances.

Second, extinction-to-emission conversion is uncertain by up to about a factor of two. $R(V)$ mixes size and composition, and the same extinction curve can admit different millimetre dielectric functions. Radiation-field, temperature, composition, and $\beta$ degeneracies remain. Optical data can regularize a fit but cannot predict the absolute polarized microwave amplitude.

Third, starlight tomography is limited by stellar density, intrinsic stellar polarization, parallax uncertainty, and a cloud-detection threshold. Starlight and emission weight the same material differently, and LOS structure can bias Serkowski-law grain constraints. Neither tracer measures the LOS magnetic-field component. Faraday rotation probes magnetized ionized gas rather than necessarily the dust-bearing phase; synchrotron angle is not a substitute for dust $\psi_i$. Conversely, using Planck 353-GHz polarization as the per-cloud angle prior is circular because it is already the emission-weighted LOS average.

Fourth, real $\beta$ statistics are non-Gaussian and correlated with dust amplitude and environment. Such correlations invalidate the simple Gaussian decorrelation formula. Simons Observatory forecasts find that high-order $\beta$-amplitude correlations can dominate component-separation bias, reaching an apparent $r\sim0.03$ in an extreme $\sigma(\beta_{\rm d})\simeq0.3$ simulation (Liu et al. 2025). A low-order moment truncation is therefore a hypothesis to test, not a guarantee.

Fifth, instrumental effects can masquerade as spectral physics: bandpass mismatch, calibration and polarization-angle errors, beam differences, temperature-to-polarization leakage, transfer-function uncertainty, and polarized CO or line emission all affect cross-frequency residuals. Current Planck conclusions depend on mapmaking: Guillet et al. find SRoll2 residual covariances compatible with predominantly temperature variation, while PR4 is compatible with neither a pure-$T$ nor a pure-$\beta$ model. CIB anisotropy contaminates total-intensity-derived dust templates, though not stellar extinction or polarization. Synchrotron complexity and residual lensing must be fitted simultaneously.

Finally, extra foreground freedom increases CMB variance and can absorb CMB modes if poorly constrained. Region selection can also induce a posteriori bias. The correct endpoint is therefore not a deterministic “3D dust subtraction,” but a joint posterior in which 3D information narrows physically plausible foreground behaviour and the remaining uncertainty is propagated into $r$.

Conclusion

Three-dimensional information helps because polarized dust is a distance-integrated vector, not because distance itself is a new cleaning frequency. Extinction tomography identifies the emitting components and their relative columns; grain physics supplies conditional priors on opacity, temperature, alignment, and polarization spectrum; starlight polarization and H I morphology constrain the per-cloud projected field angles. Their combination predicts which LOS have frequency-dependent cancellation and translates that prediction into the complex spin moments already usable by multifrequency pipelines.

The recommended implementation is a hierarchical cloud-basis fit whose absolute amplitudes remain determined by microwave data. It should output and marginalize over $W_1$, $W_2$, and angle-rotation maps, guide patch weighting, and generate physically covariant validation skies. Its benefit is conditional: it will improve $r$ inference only if cloud association, grain priors, and magnetic tomography predict held-out frequency data better than a conventional MBB model. Even then, it constrains—not eliminates—the dust foreground. That conversion of an unmodelled extrapolation error into a testable, marginalized uncertainty is the realistic contribution of 3D ISM knowledge to $r\sim10^{-3}$ CMB $B$-mode science.

References

• BICEP/Keck Collaboration, 2021, Improved Constraints on Primordial Gravitational Waves using Observations through the 2018 Season.
• BICEP/Keck Collaboration, 2022, Characterizing Dust Polarization through Correlations with Neutral Hydrogen.
• Chluba, Hill & Abitbol, 2017, Rethinking CMB foregrounds.
• Edenhofer et al., 2024, A parsec-scale Galactic 3D dust map out to 1.25 kpc.
• Green, Zhang & Zhang, 2024, The Dust Extinction Curve: Beyond $R(V)$.
• Guillet et al., 2025, Variance of dust temperature and spectral index in Planck polarization data.
• Halal, Clark & Tahani, 2024, Imprints of the Local Bubble and Dust Complexity on Polarized Dust Emission.
• Hensley & Draine, 2023, The Astrodust+PAH Model.
• Kogut & Fixsen, 2016, Foreground Bias from Parametric Models of Far-IR Dust Emission.
• Kritsuk et al., 2026, Ab initio modeling of Galactic dust polarized CMB foreground.
• Liu et al., 2025, The Simons Observatory: assessing the impact of dust complexity.
• McCallum et al., 2026, Milky Way Atlas: a velocity-resolved 3D map of local H I.
• O’Neill et al., 2024, A 3D Model of the Local Bubble’s Magnetic Field.
• Pelgrims et al., 2021, Evidence for line-of-sight frequency decorrelation of polarized dust emission.
• Pelgrims et al., 2023, Starlight-polarization-based tomography of the magnetized ISM.
• Pelgrims et al., 2024, The first degree-scale starlight-polarization tomography map.
• Planck Collaboration, 2015, Frequency dependence of thermal emission from Galactic dust.
• Planck Collaboration, 2020, Planck 2018 XI: Polarized dust foregrounds.
• Planck Collaboration, 2020, Planck 2018 XII: Galactic astrophysics using polarized dust emission.
• Poh & Dodelson, 2017, Line-of-Sight Extrapolation Noise in Dust Polarization.
• Ritacco et al., 2023, Dust polarization spectral dependence from Planck HFI data.
• Sheehy & Slosar, 2018, No evidence for dust $B$-mode decorrelation in Planck data.
• SPIDER Collaboration, 2025, Analysis of Polarized Dust Emission Using Data from the First Flight of SPIDER.
• Tassis & Pavlidou, 2015, Can dust emission properties be extrapolated from 350 to 150 GHz?.
• Vacher et al., 2023, Moment expansion generalized to spin-2 fields.
• Zelko et al., 2022, The First $3\pi$ 3D Map of ISM Dust Temperature.
• Zhang & Green, 2025, Three-dimensional maps of the interstellar dust extinction curve.

Deviations from brief

• “3D maps of thermal dust emission” has been corrected to 3D extinction-density maps plus model-dependent 3D temperature/emissivity inferences; direct distance-resolved thermal-emission measurements do not currently exist.
• The premise that present CMB $B$-mode studies demonstrably “suffer” from decorrelation has been qualified: LOS rotation is detected in selected sightlines, but broadband high-latitude decorrelation at CMB frequencies is not yet established.
• Grain size was not treated as directly determining $\beta$ in the Rayleigh regime; it primarily affects temperature and alignment, while material and structure determine the opacity slope.
• The quoted powers at $\ell=80$ were explicitly interpreted as $D_\ell=\ell(\ell+1)C_\ell/(2\pi)$, as is standard for the cited dust and tensor amplitudes.