Steelman: dedusting B-modes with 3D dust, grain physics, and the projected field

1. Strongest physical case

Decorrelation is a superposition effect, not a grain property. The polarized dust signal in a CMB pixel is a sum of complex amplitudes from distinct emitters along the beam. Each emitter has its own real spectrum and its own polarization angle. Their relative weights change with frequency, so the net angle rotates and the net Stokes Q and U acquire different effective spectral indices. In the spin-moment language of Vacher et al. 2023, the net field carries a complex spectral index whose imaginary part is the frequency derivative of the angle:

$\psi(\nu)\simeq\psi(\nu_0)+\tfrac12\,\mathrm{Im}(\Delta\beta)\,\ln(\nu/\nu_0).$

A single modified blackbody with a fixed angle cannot represent this. Neither can a real-valued moment expansion. That is the precise sense in which pixel-wise and template methods are misspecified, and the draft is right to insist on it.

The three ingredients are the three factors of each term in the sum. Write the layered model as

$Q_\nu+iU_\nu=\sum_i \tau_i\;p_{0,i}\cos^2\gamma_i\;\Big(\tfrac{\nu}{\nu_0}\Big)^{\beta_i}\tfrac{B_\nu(T_i)}{B_{\nu_0}(T_i)}\;e^{2i\psi_i}.$

The column per layer comes from 3D dust maps. The emissivity index, temperature, and polarizing efficiency per layer come from the grain size distribution, alignment physics, and local conditions. The angle and the inclination factor come from the projected magnetic field per layer. None of these factors is measurable from the CMB channels at the required precision. All three are measured by data that contain no CMB. That independence, not the extra parameters, is the strongest argument for the approach.

The high-latitude sum has few terms, and they are now catalogued. The best evidence that the problem is finite rather than a continuum:

The Local Bubble shell dominates. Skalidis and Pelgrims 2019 used stars as distance candles and found the 353 GHz polarized sky at high southern and northern latitudes dominated by a magnetized structure at roughly the shell distance. O'Neill et al. 2025 built a 3D field model on the shell from Planck 353 GHz and the 3D-dust-derived shell geometry and found the Bubble likely dominates the polarized signal along most lines of sight.
The second layer is real and lies within the mapped volume. Halal, Clark and Tahani 2024 measured line-of-sight dust complexity from the Edenhofer et al. 2024 map and found lower polarization fractions on complex sightlines at fixed column. The effect requires dust beyond the shell and disappears when the map is truncated at the shell distance.
The clouds have distances and velocities. O'Neill et al. 2026 identified about seventeen hundred dust clouds within the Edenhofer volume and assigned radial velocities to several hundred by matching morphology to 21 cm emission, placing intermediate-velocity clouds in 3D for the first time.
Position-position-velocity now maps to position-position-position. McCallum et al. 2026 reconstructed a velocity-resolved 3D atomic hydrogen map within the same volume by jointly inverting HI4PI and the Edenhofer map with information field theory, separating local from distant emission.

The mechanism has been detected with exactly the proposed tracers. Pelgrims et al. 2021 selected target sightlines from H I cloud counts and H I-derived field orientations and control sightlines with one cloud, then compared 353 and 217 GHz angles. The target sample shows excess angle change and depolarization, robust to the choice of CMB map and pipeline, and the significance rises when residual systematics are corrected in SRoll2. A cloud-to-cloud spread of order fifteen percent in the ratio of 353 to 217 GHz polarized emission reproduces the signal. That is the size of spread a layered model with different temperatures per cloud produces without tuning.

The spectral variation is probably temperature-driven, which is the predictable part. Guillet et al. 2026 analysed complex residual maps in Planck polarization with spin moments and found the SRoll2 covariances compatible with temperature variations alone rather than index variations. Temperature per layer is the quantity that a 3D density plus radiation-field model predicts, and that Zelko et al. 2023 have already mapped in 3D in intensity. Index variations tied to composition are harder to predict; temperature variations tied to shielding and density are not.

The phases are predictable, not only the power. Clark and Hensley 2019 built H I-based Stokes maps with no free parameter linking H I to dust and found them highly correlated with Planck 353 GHz Q and U, reproducing the E/B asymmetry at large scales. BICEP/Keck XVI 2023 found the same template correlated with their deep-field B-modes down to 95 GHz. A moment expansion captures the statistics of deviations from a mean SED. A layered geometric model predicts where on the sky and at which multipoles the deviations sit. That is a qualitatively different kind of information, and it is what makes the residual falsifiable.

What the strongest version does not claim. It does not claim that Planck-level decorrelation is large; Planck 2018 XI walked back the Planck Int. L detection to an upper bound. It does not claim that the layered model replaces moments; unresolved sub-beam structure still needs them. It claims that the frequency dependence of the dust angle and SED at high latitude is generated by a small number of identifiable structures whose properties are measured by non-CMB data, and that this converts a nuisance of unknown spatial form into a template with predictable form.

2. Additional mechanisms and observational handles

3D dust, beyond extinction cubes

A 3D temperature map exists. Zelko, Finkbeiner, Lee and Green 2023 assigned a temperature and index to each Bayestar19 voxel, integrated to match Planck and IRAS emission, and resolved distinct temperatures for two Cepheus components on one sightline. This is the prototype of the "3D emission map" the question asks about. It also constrains the optical-depth-to-reddening conversion per region.
A 3D grain-size proxy exists. Zhang and Green 2025 measured the extinction curve slope for 130 million Gaia XP stars and mapped it in 3D. They find accretion-dominated growth in moderately dense gas and coagulation at higher density. Through the Serkowski relation, the slope also tracks the size of the aligned grains, which links directly to the alignment size threshold in radiative-torque theory.
All-sky reach. Wang et al. 2025 combined Gaia and LAMOST into an all-sky 3D reddening map reaching several kiloparsecs at high latitude at arcminute-to-degree resolution. Gaia DR4 will improve the local maps again.
Per-component SEDs from velocity decomposition. Planck Int. XVII 2014 correlated dust with H I components over the southern cap and found temperature anti-correlated with emissivity, which is direct evidence that layers differ in SED. Lenz, Hensley and Doré 2017 provide the column-to-reddening calibration in the diffuse regime. With the 2026 3D H I map these component SEDs can be assigned distances.

Grain physics and the polarized SED

Far-infrared polarimetry near the SED peak decides between astrodust and two-composition models. Hensley and Draine 2023 argue a single composition explains a polarization fraction flat within about ten percent from 250 micron to 3 mm. Forthcoming instruments: the BLAST Observatory super-pressure balloon for wide-field diffuse dust polarimetry, PRIMA with polarimetry from 80 to 300 micron now in NASA Phase A, and the FYST 850 GHz polarimeter from the ground.
High-frequency polarization at CMB depth from Taurus at 150 to 350 GHz over most of the sky at large scales, LiteBIRD to 402 GHz, Simons Observatory at 280 GHz, and CMB-S4 at 270 GHz. These replace Planck 353 GHz as the layer calibrator at higher signal-to-noise.
Alignment modelling is becoming sightline-testable. Tram et al. 2026 reproduce Serkowski and super-Serkowski starlight spectra with radiative torques plus paramagnetic relaxation under specified radiation fields. The same machinery predicts the aligned-grain fraction per layer from the 3D radiation field.
The low-frequency SED is not settled. Planck Int. XVII found an extra component at 100 GHz of order a quarter of the emission, attributable to magnetic dipole emission, carbon dust, or flattening silicate opacity. Hensley and Bull 2018 flagged iron grains as a leading parametric bias. This argues for marginalizing over grain scenarios rather than adopting one.

Projected field per layer

Starlight tomography has produced its first map. Pelgrims et al. 2024 inverted RoboPol polarimetry over four square degrees and found most sightlines cross more than one cloud, with the plane-of-sky field orientation changing with distance on most of them. Pelgrims et al. 2023 and Uppal et al. 2026 provide the inversion methods, the latter prior-free. PASIPHAE with the WALOP polarimeters is designed to deliver hundreds of polarized stars per square degree at high latitude.
The Local Bubble field model of O'Neill et al. 2025, validated against starlight polarization and against a simulated superbubble, supplies the orientation of the dominant layer over most of the high-latitude sky. Maconi et al. 2023 show from MHD that such a shell acts as a translucent polarization filter for emission behind it. Halal et al. 2024 find no imprint of the tangential-field assumption on the polarization fraction, so this layer's inclination remains a fitted quantity.
Line-of-sight field tracers constrain the inclination factor and the large-scale geometry. Erceg et al. 2024 found Faraday structures in LOFAR data at tens of parsecs, aligned with H I filaments, at the Bubble wall. Booth et al. 2025 place a local field reversal plane within a few hundred parsecs. Lei et al. 2026 find synchrotron depolarization enhanced where H I complexity is high, which ties the synchrotron model to the same neutral layers.
Filament population models from Hervías-Caimapo and Huffenberger 2022, extended to parity-odd spectra in 2025, populate sub-beam scales with the correct E/B, TE, and TB statistics given a layer's field orientation.

Statistics and inference

Spin moments as the interface. Vacher et al. 2025 show that second-order spin moments reproduce the six-layer PySM d12 model to a few parts in ten thousand across LiteBIRD bands, and that parametric and ILC methods respond to different moment properties. The layered forward model should be delivered as predicted moment maps with covariances, not only as band templates.
Complex residual covariances from Guillet et al. 2026 are the natural observable for validating the per-layer SED assignments against Planck.
Simulation-based inference with needlet ILC compression from Duivenvoorden et al. 2025 recovered unbiased r on the d12 model at a fifth of the sky. Its training set is the natural consumer of 3D-informed simulations.
Spectral imaging with bolometric interferometry, Regnier et al. 2024, detects the r systematic from unaccounted decorrelation in a CMB-S4-like setup. It is a detector, not a corrector, and pairs naturally with a predictive model.
Non-Gaussian small-scale extrapolation with scattering transforms or generative models fills the sub-beam structure the 3D maps cannot resolve.

3. Why this could matter for r

The bias scaling. After any linear or parametric clean, the residual dust B-mode power at a CMB band is the dust power there times the fraction of that power the model does not describe. Two features make this dangerous. The residual is positive-definite, so it biases r upward rather than adding noise. And a frequency-dependent rotation leaks the brighter dust E-modes into B, so the residual BB power is about twice the naive fraction.

The dust brightness in the BICEP/Keck field sets the scale. Assumptions: the BK18 amplitude at 353 GHz and multipole 80, a modified blackbody with index 1.55 at 19.6 K, and a tensor BB spectrum of roughly 0.07 microkelvin squared at multipole 80 for r of unity.

Band	Amplitude factor from 353 GHz	Dust D_ℓ at ℓ≈80	r-equivalent
150 GHz	≈0.047	≈0.010 µK²	≈0.15
95 GHz	≈0.019	≈0.0016 µK²	≈0.024

Each percent of unmodelled dust power at 150 GHz is therefore worth about 1.5×10⁻³ in r in the cleanest deep field. LiteBIRD's target of order 10⁻³ and the CMB-S4 statistical error of 3 to 5×10⁻⁴ found by Bianchini et al. 2025 require the unmodelled fraction to stay below roughly 0.2 percent. That is about one degree of coherent unmodelled rotation, or an unmodelled effective index dispersion near 0.05.

What the current data allow. Planck 2018 XI bounds the 217 by 353 GHz BB correlation ratio above 0.975 at high confidence over its largest region. Scaled by the squared log frequency ratio to 150 GHz, that permits a correlation loss of order 15 percent in power, which is why Planck Int. L quoted a bias near 0.02 in r for a ten percent loss in the joint analysis. BK18 finds its decorrelation parameter consistent with unity and reaches a statistical error of 0.009. SPIDER finds no line-of-sight decorrelation in its field with a Gaussian-mixture test, yet measures indices in two sub-regions that differ at almost four sigma, and a joint fit with Planck that rejects a simple MBB in one sub-region. Ritacco et al. 2023 find effective angle residuals of order a degree at 217 GHz and several degrees at 100 GHz at the largest scales. Osumi et al. 2021 show Planck pixel data allow index dispersions up to 0.45 at the highest latitudes. In short, the CMB data bound decorrelation at the few-percent level and cannot themselves reach the sub-percent level the next generation needs. Twenty-fold tighter control has to come from somewhere else.

The line-of-sight mechanism in numbers. Two clouds with equal polarized intensity at 353 GHz, temperatures 17 and 21 K, indices 1.50 and 1.65, and fields misaligned by 60 degrees:

Band	Net-angle rotation relative to 353 GHz	Naive residual fraction after 353 GHz template scaling	Δr if the whole 150 GHz field behaved this way
217 GHz	≈2.9°	≈1.0%
150 GHz	≈4.8°	≈2.8%	≈4×10⁻³ naive, ≈8×10⁻³ with E-to-B leakage
95 GHz	≈6.8°	≈5.5%	≈1.3×10⁻³ naive

The Pelgrims target selection covered about six percent of high-latitude pixels under strict cuts, and Panopoulou and Lenz 2020 find two to three H I clouds per pixel on average, so milder cases are common. If a fraction f of the polarized power in a deep field comes from such sightlines, the bias is f times the table value. For f between 0.05 and 0.3 that is a few 10⁻⁴ to about 10⁻³ at 150 GHz. This is exactly the next-generation error budget, and it is a bias.

Where a physical model beats a statistical one. Liu et al. 2025 found for Simons Observatory that at high complexity the bias on r is driven by spatial correlations between the dust index and the dust amplitude at higher orders, not by the amplitude structure itself. Those correlations are physical: colder, denser, higher-column layers have different temperatures and indices. A random-field moment model treats them as nuisance variance; a layered model with per-layer physics encodes them. Vacher et al. 2022 showed that moments in both index and temperature are needed at LiteBIRD depth to remove a bias of tens of sigma, at the cost of moment-CMB correlations that inflate the error. Priors on the moment maps from the forward model shrink that inflation. Sponseller and Kogut 2022 found that line-of-sight temperature dispersion above about 1.6 K biases r for PICO with a single MBB; a 3D temperature map predicts that dispersion per pixel.

The honest bottom line. The approach reduces r bias by the fraction of line-of-sight-mixed power it identifies, reduces error inflation by replacing free moments with constrained ones, and supplies null tests at 353 and 217 GHz that no CMB-only method offers. It does not remove the need for moments, and its gain is bounded by the accuracy of the layer assignment.

4. Pipeline sketch

Inputs. Planck 353 and 217 GHz Q, U in SRoll2 and PR4; GNILC dust intensity, temperature, and index; HI4PI; the Edenhofer 2024 and Wang 2025 extinction cubes; the Zhang and Green 2025 3D extinction-curve map; the Zelko 2023 3D temperature map; the McCallum 2026 velocity-resolved 3D H I map; the O'Neill 2026 cloud catalogue; starlight polarization compilations and the RoboPol and PASIPHAE fields; the O'Neill 2025 Bubble shell and field model; Faraday-depth moment maps as an optional inclination prior.
Layer segmentation. For each sightline, decompose the extinction cube into components with a dendrogram, as in Halal et al. 2024, and merge components that share velocity in the 3D H I map into magnetically coherent layers. Cap the number at three or four. Treat the Bubble shell as the first layer where it dominates. Treat intermediate-velocity gas as a distinct layer with a distance from the cloud catalogue where available, and as distance-free otherwise.
Per-layer intensity SED. Convert differential extinction to 353 GHz optical depth with a conversion that depends on the local extinction-curve slope and density, with the Zelko constraint as prior. Assign temperature from the 3D temperature map or from an attenuated interstellar radiation field, and an index from the adopted grain model conditioned on the slope. Fit the layered sum to Planck 353, 545, and 857 GHz intensity to obtain posteriors on per-layer column and temperature, with intra-layer moments from the size distribution.
Per-layer polarization. Assign the angle from the H I orientation in the distance-matched velocity window, overridden by starlight tomography where stars exist. Assign polarizing efficiency from the grain model. Assign the inclination factor from the shell geometry for the first layer and from H I velocity coherence or Faraday priors otherwise, or absorb it into a fitted amplitude with a prior. Predict 353 GHz Q and U and compare with Planck: correlation, E/B, TE, TB, and the polarization-fraction dependence on complexity. This is the first null test.
Frequency extrapolation and compression. Evaluate the layered sum in each experiment band with bandpass integration. Compress into spin-moment maps to second order with covariances, including the angle-rotation map. Deliver both band templates and moment maps.
Likelihood options. A pixel or needlet parametric fit where the dust model is the layered sum with per-layer amplitude nuisance and small index and temperature offsets under priors. A cross-spectrum likelihood where dust cross-spectra come from the forward model plus low-order residual moments whose power is bounded by the model's prediction. A constrained ILC that deprojects the predicted moment maps. A simulation-based inference pipeline trained on 3D-informed simulations. Synchrotron is modelled separately with its own angle; H I complexity informs its depolarization but not its orientation.
Validation. Predict the Pelgrims target-minus-control excess before looking at 217 GHz. Predict the 217 by 353 GHz correlation ratio and the frequency dependence of E/B and EB from Vacher et al. 2023. Test residual covariances against Guillet et al. 2026. Apply to the BICEP/Keck and SPIDER fields, where the current answer is "no decorrelation," and check that the model predicts that too. Run end-to-end on MHD volumes with a Bubble-like cavity, observing synthetic Gaia, H I, and starlight data with realistic incompleteness, to measure the residual bias of the whole chain.
Deliverables. A line-of-sight-complexity weight map for masking or down-weighting. Per-band dust templates with uncertainties. A moment prior covariance. A public data-driven layered sky model that replaces the phenomenological six-layer model in forecast suites.

Where it fails. Wrong layer matching creates spurious rotation as easily as it removes real rotation. The inclination remains the weakest factor. Extinction-to-emission conversion carries environment-dependent errors of order ten percent per layer. Sub-beam turbulence needs moments regardless. Starlight coverage is sparse until PASIPHAE delivers. The tangential-field assumption for the shell is not yet confirmed by polarization-fraction statistics.

5. What to keep from the draft

The separation of spatial and line-of-sight decorrelation, and the explicit two-cloud algebra showing that Q and U acquire different spectra.
The statement that pixel-wise MBB fits are misspecified rather than underconstrained, with the McBride, Hensley and Bull, and Vacher references.
The insight that alignment is a high-pass filter in grain size, so polarized dust is a large-grain phenomenon and the intensity-derived index need not apply to polarization.
The degeneracy table showing which ingredient breaks which degeneracy.
The hierarchy from a 3D-informed mask through layered fits to a harmonic-space hybrid, with the mask as the minimal, immediately usable product.
The validation section, especially predicting the target-sample signal before looking at 217 GHz.
The limitations list, in particular that velocity is not distance, extinction is not emission, the inclination is unmeasured, and grain composition must be marginalized.
The closing posture: the layered model as prior and diagnostic, not as an exact sky.

6. Upgrade list

Restate the observational status of decorrelation accurately. Planck Int. L's detection was superseded by the Planck 2018 XI bound above 0.975 for the 217 by 353 GHz ratio; BK18 finds no decorrelation; SPIDER finds none with a Gaussian-mixture test but finds index variation between sub-regions at almost four sigma. Present the Pelgrims detection with its sample sizes, its p-value range, and the statement that SRoll2 corrections raise the significance.
Correct the citation for spectral imaging. The bolometric-interferometry decorrelation test is Regnier et al. 2024 in A&A 686, A271, not Maffei et al.
Adopt the spin-moment formalism as the interface between the forward model and the likelihood, and quote the angle-rotation relation from Vacher et al. 2023. Add Vacher et al. 2025 on how Planck bounds the allowed 3D complexity and how parametric and ILC methods respond differently.
Add the temperature-dominance result of Guillet et al. 2026 and connect it to the predictability of per-layer temperature.
Replace the claim that 3D emission maps do not exist with the Zelko et al. 2023 3D temperature map, and name it as the prototype of the requested ingredient.
Add the 2026 data products that answer the draft's own limitations. The velocity-resolved 3D H I map for velocity-to-distance matching, the 3D cloud catalogue with intermediate-velocity clouds located in distance, and the all-sky Gaia plus LAMOST map for reach.
Add the Local Bubble as the dominant first layer, with Skalidis and Pelgrims 2019, Alves et al. 2018, Pelgrims et al. 2020, O'Neill et al. 2025, and Maconi et al. 2023, and Halal et al. 2024 as the caution on the tangential-field assumption.
Add starlight-tomography progress. The first degree-scale map, the BISP-1 and TRIShUL inversion methods, and the PASIPHAE survey specifications.
Add the 3D extinction-curve map as a grain-size proxy and link it through the Serkowski relation to the aligned-grain size and the alignment threshold.
Quantify the r budget as in Section 3: the dust r-equivalent per band, the rotation from a realistic cloud pair, the E-to-B leakage factor, and the fraction of mixed power that must be identified.
Add BICEP/Keck XVI as evidence that H I-based templates carry B-mode phase information in the deep field, and that filamentary and diffuse dust share an index there.
Add the physical-correlation argument from Liu et al. 2025 and the Sponseller and Kogut 2022 temperature-dispersion threshold, and the Osumi et al. 2021 point that Planck alone leaves high-latitude index dispersion nearly unconstrained.
Temper the index offset claim. Planck 2018 XI gives a polarized-minus-intensity index difference of 0.05 with an error of 0.03, so it is marginal, not established. Note that Vinzl et al. 2026 forecast LiteBIRD can measure E, B, and intensity indices separately.
Add the far-infrared and high-frequency roadmap and state what each instrument decides: composition from BLAST Observatory, PRIMA, and FYST; layer calibration from Taurus, LiteBIRD, Simons Observatory, and CMB-S4.
Add Faraday tomography as an inclination and geometry handle, and Lei et al. 2026 as the link between neutral-gas complexity and synchrotron depolarization.
Add the end-to-end requirement explicitly: a data-driven layered sky model, validated on MHD with synthetic ancillary data, consumed by both parametric pipelines and simulation-based inference.
Drop the unattributed claim that the fifteen percent ratio spread corresponds to a ten percent SED shape change per a named group; state only what Pelgrims et al. report.
Add the per-component SED evidence from Planck Int. XVII and the unexplained 100 GHz component, and the phase-resolved small-grain evidence from Hensley, Murray and Dodici 2022, as reasons to model grain populations per layer rather than per sky.