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

Three-dimensional, grain-aware magnetic modeling of polarized dust for CMB $B$-mode cleaning

Polarized thermal dust as a $B$-mode foreground

Large interstellar grains absorb optical and ultraviolet starlight and reradiate the energy in the far-infrared and millimetre bands. Because some grains are nonspherical and partially aligned with the Galactic magnetic field, their thermal emission is linearly polarized. For the usual alignment geometry, the electric-vector position angle of the emitted radiation is perpendicular to the plane-of-sky magnetic-field orientation; polarization produced by dichroic extinction of starlight is rotated by approximately $90^\circ$ relative to the emission polarization.

For optically thin dust, a common intensity model is a modified blackbody (MBB),

where $T_{\rm d}$ is an effective grain temperature and $\beta_{\rm d}$ describes the frequency dependence of the opacity. In the simplest polarization model,

up to the adopted Stokes-sign convention. Here $p_{\rm d}$ is the polarization fraction and $\psi_{\rm d}$ is the polarization angle. Planck found that a single MBB with $T_{\rm d}\simeq19.6$ K and mean polarized spectral index $\beta_{\rm d}^{P}\simeq1.53$ described its broad-region frequency spectra reasonably well, although such averages do not establish that every line of sight has a simple SED (Planck Collaboration XI 2020).

Dust matters for $B$-mode cosmology because its spin-2 $Q,U$ field contains substantial $B$-mode power at the angular scales of both the reionization and recombination signatures of primordial gravitational waves. The primordial contribution is proportional to the tensor-to-scalar ratio $r$, while the dust contribution has no reason to be negligible even in low-column-density fields. Changing the dust polarization amplitude or angle also changes its nonlocal $E/B$ decomposition. A small unmodeled dust residual can therefore mimic positive primordial $BB$ power or broaden—and sometimes bias—the posterior on $r$.

The CMB has a known blackbody spectrum and is frequency-independent when maps are expressed in thermodynamic-temperature units. Component separation exploits this fact. Dust cleaning fails when the assumed dust frequency mixing law is too restrictive, not simply because dust is bright.

What frequency or “spectral-index” decorrelation means

The phrase spectral-index decorrelation is somewhat imprecise: it is the polarized dust maps at different frequencies that decorrelate, rather than a scalar $\beta_{\rm d}$ somehow losing correlation. There are several related operational effects.

A single-MBB-per-pixel model assumes

with a real scalar $f_\nu$. Consequently:

1. $Q$ and $U$ have the same frequency scaling;
2. $p_{\rm d}$ and $\psi_{\rm d}$ are frequency-independent;
3. a unique pair $(\beta_{\rm d},T_{\rm d})$ describes all frequency pairs; and
4. dust occupies one spectral degree of freedom, or one column of the component-separation mixing matrix.

Actual emission is a vector sum over material along the line of sight:

Here $i$ denotes a cloud or distance layer, $A_i$ measures its dust column, $s_i(\nu)$ is its intensity SED, and $\gamma_i$ is the field inclination relative to the plane of the sky. If the layers have different $\beta_i$, temperatures, compositions, size distributions, or alignment efficiencies, their relative polarized weights change with frequency. If their $\psi_i$ also differ, the complex sum rotates as frequency changes:

In a two-cloud example, a frequency-dependent polarization angle requires both a change in the clouds’ relative polarized intensities and nonparallel projected fields. Different SEDs without field misalignment produce spectral curvature but no angle rotation; misaligned fields with identical SEDs still factorize into a single frequency law. This two-condition mechanism was emphasized by Tassis and Pavlidou and by Poh and Dodelson (Tassis & Pavlidou 2015; Poh & Dodelson 2017).

Even without angle rotation, sums of MBBs are generally not themselves exact MBBs. The inferred $\beta_{\rm d}^{\rm eff}$ can depend on the frequency pair, and $Q$ and $U$ can acquire different effective indices once differently oriented layers are mixed. Spatial variation of the SED within a beam or analysis region additionally changes the dust morphology under frequency scaling.

At the power-spectrum level the effect is commonly summarized by

A perfectly coherent dust template has $\rho_\ell^{BB}=1$, whereas genuine frequency-dependent morphology gives $\rho_\ell^{BB}<1$. Noise, CMB fluctuations, masks and instrumental systematics complicate estimation of this ratio.

This violates both template cleaning,

and likelihood models that factorize the dust cross-spectrum as $C_\ell^{\rm d}(\nu_i,\nu_j)=A_{\rm d}f_i f_j$. In the former case a dust residual remains; in the latter the cross-spectrum can underestimate the foreground power and bias $r$ upward.

The observational situation is nuanced. An intermediate Planck analysis reported broad-region decorrelation between 217 and 353 GHz, but subsequent analyses identified sensitivity to noise and systematics, and the final Planck polarized-foreground analysis found no significant broad-sky decorrelation at its sensitivity (Planck Intermediate Results L; Sheehy & Slosar 2018; Planck Collaboration XI 2020). In contrast, Pelgrims et al. selected sightlines independently using multiple, magnetically misaligned H I clouds and detected enhanced 217–353 GHz angle differences. They estimated that roughly $15\%$ cloud-to-cloud variation in the polarized 353/217 GHz ratio could reproduce the signal (Pelgrims et al. 2021). Thus decorrelation is a physically established possibility and has been detected in selected environments, but it is not accurately represented as a universal, already-measured scalar coefficient.

What each additional ingredient contributes

Three-dimensional thermal-dust information

A three-dimensional dust map replaces the projected column $A_{\rm d}(\hat{\boldsymbol n})$ by a differential distribution $A_{\rm d}(\hat{\boldsymbol n},s)$, or ideally a frequency-dependent emissivity $j_\nu^{\rm d}(\hat{\boldsymbol n},s)$. Its independent value is that it identifies how much emitting material lies in each distance interval.

This information can:

• distinguish genuinely simple sightlines from those containing several comparable clouds;
• assign separate amplitudes, temperatures or SED parameters to physically distinct structures;
• prevent one cloud’s fitted $\beta_{\rm d}$ from being interpreted as a property of the entire column;
• provide priors on which layers need independent polarization components;
• identify low-risk fields dominated by a single nearby structure; and
• connect dust temperature and grain evolution to environment, radiation field and gas phase.

Models explicitly allowing $T_{\rm d}$, $\beta_{\rm d}$, density and polarization properties to vary in three dimensions already demonstrate how such information generates more realistic component-separation tests (Martínez-Solaeche, Karakci & Delabrouille 2018).

A crucial qualification is that most present “3D dust maps” are maps of stellar reddening or extinction, not direct tomography of thermal emission. Gaia parallaxes combined with stellar photometry can localize reddening—for example, the Green et al. map extends over much of the sky to several kiloparsecs—but converting differential reddening into millimetre emissivity requires environment-dependent opacity and temperature assumptions (Green et al. 2019). Continuum thermal emission itself carries essentially no distance label. H I or molecular-line velocities, stellar distances and extinction must therefore supply the tomographic information.

On its own, a 3D intensity map cannot predict $Q$ and $U$: it gives the scalar weights of the layers but neither their polarization efficiencies nor their vector orientations.

Grain-size distribution and grain physics

For grain material $m$ and radius $a$, a physical emissivity is schematically

where $\mathcal U$ characterizes the local radiation field. The polarized emissivity replaces the absorption cross-section by a polarization cross-section and includes the alignment function $R_m(a)$:

Knowledge of $dn_m/da$ therefore helps determine:

• which grain sizes dominate total and polarized emission;
• the distribution of equilibrium or stochastic temperatures;
• which sizes are efficiently aligned;
• the frequency dependence of $p(\nu)=j_\nu^P/j_\nu^I$;
• departures from a pure opacity power law; and
• physically allowed correlations among $T_{\rm d}$, $\beta_{\rm d}$, composition and alignment.

The polarized SED need not equal the total-intensity SED because only a subset of the grain population may be aligned. Draine and Fraisse showed that plausible mixtures of silicate and carbonaceous grains produce significantly different polarization spectra even while satisfying extinction constraints (Draine & Fraisse 2009). Later DustEM-based models demonstrated strong degeneracies among composition, shape, porosity and alignment in reproducing Planck intensity and polarization (Guillet et al. 2018).

Grain information thus supplies a physical family of layer SEDs and informative priors, rather than treating every $\beta_i$ as an unrelated nuisance parameter. It is particularly valuable for avoiding the common assumption that $T_{\rm d}$ and $\beta_{\rm d}$ inferred from total intensity automatically apply to polarization.

Size distribution alone is nevertheless insufficient. At CMB wavelengths, large grains are usually in the Rayleigh regime, so the millimetre opacity slope may depend more strongly on dielectric properties, composition and internal structure than on radius itself. Grain size remains important through heating and size-dependent alignment, but it does not uniquely determine $\beta_{\rm d}$. Fanciullo et al. found that size changes alone could not explain observed emission-to-extinction polarization behavior; alignment and magnetic geometry also had to vary (Fanciullo et al. 2017).

The sky-projected magnetic field

The plane-of-sky field determines the phase $e^{2i\psi_i}$ of each layer’s polarized emission. Its inclination and unresolved fluctuations control geometrical depolarization through factors such as $\cos^2\gamma_i$ and through vector cancellation. Consequently, magnetic information predicts where a scalar dust-frequency rescaling is most likely to fail: sightlines containing comparable layers with strongly different $\psi_i$ are intrinsically vulnerable.

A projected-field map can come from:

• high-frequency dust polarization itself, after rotating the emission angle by $90^\circ$;
• optical or near-infrared starlight polarization;
• magnetically aligned H I structures in velocity-channel maps; or
• a regular-plus-turbulent Galactic magnetic-field model.

The first is high signal-to-noise where dust is bright but is already integrated along the line of sight. Stellar polarization becomes tomographic when stars at different Gaia distances bracket individual clouds. Bayesian methods developed for PASIPHAE explicitly infer the mean plane-of-sky field orientation and intrinsic scatter of distance-resolved layers (Pelgrims et al. 2023). Velocity-resolved H I morphology offers a complementary proxy and was central to the targeted Planck decorrelation detection.

Magnetic information alone does not determine a layer’s polarized amplitude or SED. Moreover, an ordinary two-dimensional field map gives only the net line-of-sight orientation—the quantity that changes under decorrelation—so distance resolution is highly desirable.

A concrete three-component dedusting strategy

The three ingredients naturally define a layer-based forward model. They provide, respectively, where and how much dust exists, how each layer emits and polarizes as a function of frequency, and in which Stokes direction that polarized radiation is added.

First, partition every line of sight into a modest number of physically motivated layers using differential extinction, thermal-emission morphology and H I or molecular velocity components. For layer $i$, infer a dust-column amplitude $A_i$, distance distribution and environmental variables such as radiation-field strength and gas phase. The number of statistically independent layers should be selected or regularized; it cannot exceed what the frequency data can support.

Second, use extinction curves, polarized extinction, total dust emission and polarized high-frequency measurements to constrain layer-level grain parameters

These parameters generate bandpass-integrated intensity and polarized emissivities $S^I_{bi}(\boldsymbol\theta_i)$ and $S^P_{bi}(\boldsymbol\theta_i)$. A practical analysis need not fit a full microphysical grain model in every pixel: one can construct a library of physical SEDs, reduce it to a few principal components, or translate it into priors on $(T_i,\beta_i)$, spectral curvature and $p_i(\nu)$.

Third, assign each layer a projected angle $\psi_i$ and priors on inclination, alignment and turbulent dispersion using starlight-polarization tomography, H I morphology and high-frequency polarization. The predicted dust signal in band $b$ is then

where $g_i$ includes inclination, alignment and sub-beam depolarization. The corresponding data model is

with beam convolution, colour corrections and the measured bandpass included explicitly.

This model should be fitted jointly and hierarchically. The external 3D, grain and magnetic data supply probabilistic priors rather than an immutable dust template. Microwave data update the layer amplitudes and SEDs, while uncertainties and covariances propagate into the CMB map and ultimately the $r$ likelihood. Dedusting then means either subtracting the posterior mean dust signal while retaining its covariance or, preferably, marginalizing over dust jointly with $r$.

When the number of physical layers is too large, the same model can be compressed into a complex moment expansion:

The complex moments contain angle-weighted sums such as ${\cal M}_{\beta}\propto\sum_i A_i(\beta_i-\bar\beta)e^{2i\psi_i}$. They therefore represent both SED curvature and different frequency behavior in $Q$ and $U$. Three-dimensional and grain data determine which moments should be retained and constrain their priors; magnetic tomography constrains their complex phases. Moment-expansion methods have been developed by Mangilli et al. and Vacher et al., with simulations showing that omitted $\beta_{\rm d}$ and $T_{\rm d}$ moments can bias $r$ at next-generation sensitivities (Mangilli et al. 2021; Vacher et al. 2022).

The result must be tested against frequencies excluded from the fit, frequency-dependent polarization-angle maps, $EE$, $BB$, $EB$ and $TB$ spectra, and cross-correlations of cleaned maps with layer-resolved dust templates. Regions predicted to contain several misaligned clouds provide especially powerful stress tests. Simulations should include the same beams, masks, bandpasses, calibration errors and mapmaking systematics as the observations.

Limitations and residual systematics

Several obstacles prevent this approach from being a deterministic reconstruction.

• Incomplete tomography: Stellar sampling becomes sparse with distance and extinction, while high-latitude low-column regions have intrinsically weak reddening and polarization. Distance resolution is finite, and thermal continuum emission has no direct distance coordinate. H I velocity is not always a unique distance indicator.

• Extinction-to-emission conversion: Differential reddening traces an optical cross-section, not millimetre emissivity. The conversion varies with grain composition, size, $R_V$, radiation field and gas environment. Cosmic infrared-background fluctuations also contaminate faint total-intensity dust maps.

• Grain-model uncertainty: Composition, optical constants, porosity, shape, magnetic inclusions and alignment physics remain degenerate. Temperature mixing can masquerade as changes in $\beta_{\rm d}$, and noise produces a familiar $T_{\rm d}$–$\beta_{\rm d}$ degeneracy. The aligned population may differ substantially from the population dominating $I$.

• Magnetic-field incompleteness: Most tracers constrain only the plane-of-sky orientation, not field inclination or strength. Starlight extinction and thermal emission weight the column differently. H I filaments are useful but not universal field tracers, and turbulence below the angular or distance resolution produces irreducible beam and line-of-sight depolarization.

• Parameter identifiability: A CMB experiment has a finite number of frequency bands. Allowing every cloud independent $T_i,\beta_i,p_i,\psi_i$ rapidly creates an underdetermined fit. Physical priors, layer compression and moment truncation are therefore essential, but incorrect priors can themselves bias $r$.

• Instrumental effects: Polarization-angle calibration, temperature-to-polarization leakage, bandpass mismatch, beam mismatch and correlated noise can imitate SED differences or angle rotations. These effects also bias empirical estimates of $\rho_\ell^{BB}$.

• Other sky components: Synchrotron spectral curvature and spatial variation, dust–synchrotron correlation, lensing $B$-modes, polarized compact sources and possibly weak polarized anomalous emission must be modeled simultaneously.

A single MBB does not necessarily become unusable whenever several temperatures exist: finite-band data may be insensitive to the resulting curvature, and simulations find that it can remain approximately unbiased in the absence of polarization-angle mixing. The structural failure occurs when no common real scaling of $Q$ and $U$ exists; line-of-sight depolarization combined with different layer SEDs is the clearest example (Sponseller & Kogut 2023).

Thus 3D dust maps, grain-size information and projected magnetic-field measurements do not eliminate foreground uncertainty. Their value is that they replace an uncontrolled, phenomenological decorrelation parameter with a physically structured vector radiative-transfer model whose degrees of freedom can be constrained, validated and marginalized in the inference of $r$.