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

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

1. Polarized thermal dust and why $B$-mode studies care

The primordial $B$-mode signal sought by CMB experiments is a spin-2 polarization pattern generated by tensor perturbations. Its amplitude is conventionally described by the tensor-to-scalar ratio $r$, with a recombination feature near $\ell\simeq80$ and a reionization feature at $\ell\lesssim10$. Galactic thermal dust also produces $B$-mode polarization on these scales. At frequencies above roughly $100$ GHz it is usually the dominant polarized Galactic foreground, while gravitational lensing produces a separate, cosmological $B$-mode contaminant that must be treated by delensing rather than “dedusting.”

Large, nonspherical grains are heated by the interstellar radiation field and radiate in approximate thermal equilibrium. For an optically thin, homogeneous dust slab, the intensity is commonly represented by a modified blackbody (MBB),

Here $\tau_{\nu_0}$ is the optical depth at reference frequency $\nu_0$, $T_d$ is a colour temperature, and $\beta_d$ is an effective opacity index. None is fundamentally guaranteed to be constant through a telescope beam or along a line of sight.

Radiative-alignment torques, rather than classical Davis–Greenstein relaxation alone, provide the leading explanation for grain alignment. The short axes—or, statistically, angular momenta—of sufficiently large grains tend to align with the local magnetic field. Their thermal emission is consequently polarized with electric-vector position angle perpendicular to the plane-of-sky magnetic field $B_{\rm POS}$. For a uniform slab,

and schematically

where $\gamma$ is the angle between the magnetic field and the line of sight. Tangling within the beam or along the line of sight reduces the observed $p$ further. Planck observations directly connect low polarization fraction with polarization-angle disorder and hence with three-dimensional magnetic geometry; they also confirm the $90^\circ$ relation between starlight extinction polarization and dust-emission polarization Planck 2018 XII.

Dust is not merely an intensity foreground whose amplitude can be removed before forming $E$ and $B$. Its spatially varying $Q,U$ field has intrinsic $E$- and $B$-mode power. Planck found an effective polarized-dust index $\beta_d^P\simeq1.53\pm0.02$ for fixed $T_d=19.6$ K, approximately $C_\ell^{EE}/C_\ell^{BB}\simeq2$, and $C_\ell^{EE,BB}\propto\ell^{-2.4}$ over broad high-latitude masks Planck 2018 XI. These are ensemble summaries, not a universal SED or a material constant. In the BICEP/Keck field, the measured dust $BB$ power at 150 GHz exceeds the primordial power associated with $r\sim10^{-3}$ by orders of magnitude in power. LiteBIRD’s formal requirement, for example, is a total uncertainty $\delta r<10^{-3}$, including foreground and systematic errors LiteBIRD Collaboration. Thus even percent-level errors in extrapolating dust between frequency bands can be cosmologically important.

2. What spectral-index decorrelation means operationally

“Spectral-index decorrelation” is used for several related but non-identical effects. Component-separation requirements make the distinctions important.

Map-level non-proportionality

If one dust component with one angle describes a pixel, its polarization maps obey

Thus the two maps differ only by a real scalar, locally or globally depending on the model. In intensity, an effective two-frequency index relative to an assumed temperature $T_\star$ can be defined as

Spatial variation of $\beta_{\rm eff}$ makes a single all-sky scaling fail. That is not yet necessarily line-of-sight decorrelation: a sufficiently resolved one-MBB-per-pixel model may absorb spatial variation across the sky.

For polarization, $Q$ and $U$ may have genuinely different frequency dependences. Because either Stokes parameter can cross zero, a logarithmic “$\beta_Q$” or “$\beta_U$” is only locally meaningful. The robust operational statement is that no common real $f_d(\nu)$ makes both $Q_\nu=f_d(\nu)Q_{\nu_0}$ and $U_\nu=f_d(\nu)U_{\nu_0}$.

Spectral rotation

The observed polarization angle

can vary with frequency. This “spectral rotation” immediately proves that the polarization at two frequencies is not related by a positive scalar. It is especially dangerous because even a small frequency-dependent rotation mixes the dust $E$ and $B$ patterns.

Power-spectrum decorrelation

For $X=E$ or $B$, one measures the correlation ratio

After correct removal of noise and CMB biases, $R_\ell<1$ means that the two dust maps do not share exactly the same morphology. Spatially varying SEDs, line-of-sight mixing, beam averaging, or instrumental errors can all lower the measured ratio. Conversely, $R_\ell\simeq1$ over a large mask does not exclude localized rotation or decorrelation below the sensitivity relevant to $r=10^{-3}$.

Mixing-matrix rank

In thermodynamic CMB units, a rank-one dust model has

Across frequency, all dust pixels lie in the spectral subspace spanned by one vector $f_d(\nu)$. Commander-like parametric fits, and template cleaning that scales 353-GHz $Q,U$ with one MBB, rely on this structure. Line-of-sight components with distinct SEDs introduce additional spectral vectors:

Dust is now higher rank. If a component-separation model supplies too few dust degrees of freedom, the unmatched spectral mode can project onto the frequency-independent CMB column and bias $r$. An ILC is not automatically immune: unmodelled dust rank consumes frequency degrees of freedom and can remain in the minimum-variance combination.

The two-cloud mechanism

Let two clouds contribute

Their Stokes parameters are

Define the local logarithmic SED slopes

if $p_i$ is frequency independent. Where $Q,U\neq0$,

and therefore

The angle rotates at the rate

where primes denote $d/d\ln\nu$. Thus rotation requires both magnetic misalignment and unequal spectral slopes. It can become large near cancellation, where $|P|$ is small. Geometrically, $P(\nu)$ traces a curve rather than a fixed ray in the $Q,U$ plane. Reusing the 353-GHz polarization angle when extrapolating to 150 GHz is then biased.

This mechanism was proposed explicitly by Tassis and Pavlidou 2015. Pelgrims et al. used H I cloud counts and cloud-dependent field orientations to select adverse, magnetically misaligned sightlines before examining Planck polarization. They detected a 217-to-353-GHz angle/depolarization difference in the selected sample and found that approximately $15\%$ cloud-to-cloud variation in the polarized 353/217 ratio could reproduce it Pelgrims et al. 2021. That is compatible with Planck’s near-unity $R_\ell^{BB}$ over broad high-latitude masks: the former is a targeted pixel-level test, while the latter averages over mostly less adverse sightlines and is sensitive to noise and systematics. SPIDER likewise reported no significant full-field line-of-sight decorrelation, although it found significant subregion-dependent SED behaviour SPIDER Collaboration. Decorrelation is therefore a spatially heterogeneous risk, not an established uniform loss of correlation everywhere.

3. What the three physical ingredients contribute

A. Three-dimensional dust maps: locating the emitting material

Thermal emission is optically thin at CMB frequencies, so an ordinary intensity map supplies a column integral, not a distance. Useful “3D emission maps” are consequently reconstructions: stellar reddening and Gaia distances locate dust in voxels, while FIR/submillimetre intensities constrain the temperatures, opacities, and integrated emission of those voxels.

Bayestar19, for example, infers three-dimensional reddening from Gaia, Pan-STARRS, and 2MASS stellar data Green et al. 2019. Zelko et al. assigned $T_d$, $\beta_d$, and emission-to-reddening conversion factors to Bayestar voxels and required their line-of-sight sums to reproduce Planck and IRAS intensity maps Zelko et al.. Such a reconstruction can identify whether a pixel is dominated by one nearby cloud or contains several comparably weighted layers, and provides initial optical-depth weights $\tau_i$.

At high latitude, H I 21-cm velocity cubes provide a complementary decomposition of cold diffuse gas. Velocity is not a unique distance coordinate locally, but distinct velocity components can still identify separate clouds and their relative column densities. Filament orientations within velocity slices also carry magnetic information. These products should be treated as probabilistic layer bases, not exact Cartesian maps.

The central gain is that an integrated one-MBB fit is replaced by a sum over physically located emitters. A 3D map alone, however, models only where dust exists. It cannot determine the polarization direction of each layer, and extinction alone does not say how strongly a voxel emits.

B. Grain size distribution: assigning polarized SEDs

For a layer with size distribution $n_i(a)$, the physically relevant emissivities are schematically

The temperature average permits stochastic heating for nanoparticles. The polarized SED is differently weighted from total intensity because alignment efficiency generally rises toward the large-grain tail. Extinction-curve shape, $R_V$, FIR opacity per hydrogen atom, gas phase, and radiation-field strength can therefore supply priors on $T_i$, opacity shape, and $f_{{\rm align},i}$. A mean Milky Way distribution such as MRN or Weingartner–Draine is only a population prior; it is not a measurement of $n_i(a)$ in each cloud.

In the Rayleigh limit,

with $\beta$ controlled mainly by the material dielectric function, internal structure, and composition. Diffuse-ISM grains of radius $\sim0.05$–$0.2\,\mu{\rm m}$ are extremely deep in the Rayleigh regime at CMB wavelengths, so their MBB index is not determined by geometrical “greying.” Grain size matters strongly through equilibrium temperature, stochastic heating, alignment, aggregation, porosity, and correlated changes in composition. It changes $\beta_d$ directly through finite-size greying only if grains grow toward tens to hundreds of micrometres at FIR/submillimetre wavelengths.

Coagulation, mantling, and shattering consequently do not imply a universal sign for $\Delta\beta_d$. For example, one detailed dense-cloud evolution model predicts that accretion and aggregation lower $T_d$ while increasing the FIR/submillimetre index from about $1.5$ toward $2$ Köhler, Ysard & Jones. The safer inference is hierarchical:

not a deterministic relation $n(a)\mapsto\beta_d$.

Modern astrodust-plus-PAH models illustrate why intensity and polarized SEDs need not coincide: nanoparticle emission contributes to intensity, whereas aligned larger grains dominate FIR/submillimetre polarization Hensley & Draine. Anomalous microwave emission from rotating nanoparticles is a separate, mostly unpolarized component and should not be folded into the thermal-dust $\beta_d$.

C. $B_{\rm POS}$: assigning polarization geometry

For layer $i$,

where $D_i\le1$ represents unresolved field tangling. The plane-of-sky field determines $\psi_i$; the inclination and alignment physics determine the amplitude. Misaligned layers depolarize one another, and unequal SEDs make that depolarization frequency dependent.

Starlight-polarization tomography is especially valuable because foreground stars sample only nearby dust while background stars include additional layers. Combining stellar polarization with Gaia distances permits a differential inference of cloud polarization. PASIPHAE’s Bayesian inversion explicitly fits thin dust layers, their mean angles, and turbulence-induced scatter while propagating parallax and polarization uncertainties Pelgrims et al. 2023. The $180^\circ$ field ambiguity is harmless for $e^{2i\psi}$.

H I filament orientation supplies another statistical $B_{\rm POS}$ tracer as a function of velocity. Clark–Hensley-style H I pseudo-Stokes maps reproduce much of the integrated 353-GHz morphology and enabled the Pelgrims sightline selection. Faraday rotation and synchrotron provide complementary information on $B_{\rm LOS}$, electron-weighted field structure, and sometimes inclination, but they do not replace per-cloud $B_{\rm POS}$: their emitting and rotating media weight the line of sight differently from aligned dust.

4. A concrete three-dimensional dedusting strategy

A practical implementation should be hierarchical and uncertainty-aware rather than a deterministic subtraction.

1. Construct a layer basis. For each analysis pixel or superpixel, combine stellar-extinction tomography, reconstructed 3D thermal emission, and H I velocity components to define $N$ statistically supported layers. Infer posterior distributions for distances and optical-depth weights $\tau_i$. Merge layers that cannot be distinguished by the available frequency channels; the aim is the recoverable spectral rank, not maximal spatial granularity.

2. Assign physical SED priors. Use local extinction-curve information, gas phase, radiation field, opacity per $N_H$, and a grain model to obtain correlated priors on $(T_i,\beta_i,\kappa_i,f_{{\rm align},i})$. Allow the polarized SED to differ from the intensity SED. Where grain information is only population-level, encode its scatter rather than fixing all layers to a Galactic mean.

3. Assign per-layer geometry. Infer $\psi_i$ and its turbulent dispersion from starlight tomography and H I orientation. Use 353-GHz $Q,U$ as an integrated constraint that the layer sum must reproduce—not as a truth map that forces the same angle at every frequency. Marginalize over $\gamma_i$, alignment efficiency, and unresolved depolarization unless independent inclination information exists.

4. Forward-model the instrument channels. For each band $b$, integrate the physical spectrum through the measured bandpass:

The full data model adds CMB, synchrotron, noise, calibration, bandpass mismatch, beams, and polarization-angle uncertainties. Intensity $I$ can constrain $\tau_i,T_i$, but should enter with explicit models for the cosmic infrared background, zodiacal emission, and intensity-to-polarization leakage.

5. Perform component separation in one of three equivalent representations.

   • In a multi-layer parametric fit, replace the usual single dust amplitude by a small number of layer amplitudes or physical spectral modes. External 3D information supplies informative—but not infinitely narrow—priors. The CMB posterior is obtained jointly with all dust parameters, as in a physical extension of Commander.

   • In a localized constrained ILC, derive from the 3D posterior the leading dust spectral vectors or tangent vectors, such as $f_d$, $\partial f_d/\partial\beta$, $\partial f_d/\partial T$, and a rotation mode. Choose weights satisfying $w^Ta_{\rm CMB}=1$ and $w^TF_d=0$ for the supported dust subspace. Spatial localization is essential because the layer mixture varies over the sky. Each deprojected column increases map noise, so the 3D prior is useful precisely because it identifies which extra ranks are required locally.

   • In a subtract-then-expand analysis, subtract the posterior mean 3D forward model and describe only its residual with MBB moments. Schematically,

 For polarization, the coefficients are complex spin moments, so they encode both SED distortions and frequency-dependent rotation. Moment expansion is a flexible phenomenological description of unresolved mixing [Chluba, Hill & Abitbol](https://arxiv.org/abs/1701.00274); its spin-2 generalization makes spectral rotation explicit [Vacher et al.](https://arxiv.org/abs/2205.01049). The physical 3D model should remove the leading moments and regularize the rest. Forecasts show that moments can remove $r$ bias but enlarge $\sigma(r)$ through CMB–moment degeneracies [Vacher et al.](https://ora.ox.ac.uk/objects/uuid%3Ad5c42dfd-02d3-44f8-999d-0f2959f74ecb).

6. Use geometry-aware weighting. Sightlines containing one dominant cloud with coherent $B_{\rm POS}$ can tolerate a simpler model. Multi-cloud, strongly misaligned sightlines require more spectral freedom. Pixels with insufficient stellar sampling, ambiguous H I decomposition, or near-cancellation in $Q,U$ should be downweighted or masked rather than forced into one MBB.

7. Propagate the foreground posterior into $r$. Simulations must draw distance errors, layer columns, grain properties, angles, turbulent scatter, calibration, and noise jointly. The required outputs are recovered-$r$ coverage, $\Delta r$, residual $BB$ spectra, and frequency null tests—not merely a visually improved dust map.

The joint necessity is explicit. A 3D mass map without $B_{\rm POS}$ cannot predict spectral rotation. $B_{\rm POS}$ without a distance decomposition supplies only an integrated angle. Grain physics without 3D structure gives a mean SED but not the cloud-to-cloud contrast that drives decorrelation. Finally, 3D structure plus field angles but one identical SED for all layers reproduces depolarization while missing its frequency dependence.

5. Limitations and residual systematics

The most important limitation is that none of the required three-dimensional fields is directly and completely observed. Extinction traces total dust column; thermal emission weights dust by opacity and temperature; polarized emission weights only the aligned, anisotropically emitting subset. Their voxel weights cannot be assumed identical.

High-latitude sightlines—the most important ones for primordial $B$-modes—also contain few stars and little reddening, making extinction and polarization tomography noisy. Distance resolution is generally coarser than physical cloud thickness. H I velocity is only a proxy for depth, molecular or ionized dusty gas may be missing, and filament orientation is not a perfect field tracer. $B_{\rm POS}$ gives an angle but not $\gamma$, leaving alignment efficiency, inclination, and dust column partly degenerate.

Per-cloud grain-size distributions and compositions are rarely measured. Extinction curves may be unavailable or dominated by different grain populations from polarized emission. MBB parameters are correlated fit summaries: temperature mixing itself produces a non-MBB spectrum, and noise creates $T$-$\beta$ degeneracy. Laboratory uncertainties in amorphous-grain optical constants and magnetic inclusions propagate directly into long-wavelength extrapolation.

A useful first-order residual diagnostic is

The imaginary angle-error term is qualitatively different from an amplitude or index error and can create $B$-mode residuals even when the polarized intensity is well fitted. Near line-of-sight cancellation, small errors in layer amplitudes or angles can become large fractional errors in the integrated polarization.

Instrumental effects can imitate the same observables: bandpass mismatch, frequency-dependent beams, polarization-efficiency errors, absolute-angle calibration, intensity-to-polarization leakage, and residual CMB in 217-GHz dust estimates. CIB fluctuations complicate intensity-derived SED priors. Synchrotron has its own spatial and possibly line-of-sight spectral complexity, while dust–synchrotron correlation changes the likelihood. Dust $EB$ caused by filament–field misalignment is related magnetized-ISM physics but is not synonymous with $\beta_d$ decorrelation.

Validation should therefore combine:

• posterior predictive tests in every frequency band;
• cross-spectra between independent data splits;
• angle-difference tests targeted by 3D geometry;
• results with and without the most complex sightlines;
• simulations containing more layers than the fitted model;
• independent pipelines using parametric, constrained-ILC, and moment methods;
• explicit tests for $E$-to-$B$ leakage and calibration errors.

For a BICEP-class patch, dense targeted starlight polarimetry and H I decomposition could identify a manageable number of adverse clouds and support a small-layer likelihood. CMB-S4 would need the method across multiple ground-accessible fields, probably using spatially varying spectral bases rather than literal voxel-by-voxel fitting. LiteBIRD’s full-sky, low-$\ell$ analysis would use 3D information primarily as priors, complexity masks, and moment regularization, because uniform high-quality tomography will not exist over the whole sky.

Three-dimensional information therefore does not make dust removal exact. Its value is that it turns decorrelation from an unexplained extra spectral rank into a physically predicted and uncertainty-quantified one. Success is demonstrated only if realistic end-to-end realizations recover unbiased $r$ with correct coverage and a residual foreground budget safely below the experiment’s allocated $\Delta r$.

Deviations from brief

• The brief’s claim that ordinary diffuse-ISM grain growth directly lowers $\beta_d$ by approaching $\lambda/2\pi$ is not applicable at CMB wavelengths: submicron grains remain deep in the Rayleigh limit. Composition, structure, mantling, and aggregation dominate $\beta_d$; coagulation can increase rather than decrease it.
• “3D thermal-emission maps” are not direct distance-resolved emission measurements. Present products generally combine stellar extinction tomography with line-of-sight-integrated FIR data and model-dependent voxel assignments.
• Current observations do not establish strong, uniform decorrelation in all $B$-mode fields. Broad Planck masks and the full SPIDER field are consistent with little decorrelation, while 3D-selected, magnetically adverse Planck sightlines show a localized effect. The problem is a precision-limiting systematic risk rather than a universally detected failure.