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

Steelman: 3D dust, grain physics, and projected $B$ as structured control of polarized-dust frequency decorrelation

The draft’s core claim is already the right one: frequency “decorrelation” of polarized dust is not a mysterious scalar $\beta_{\rm d}$ that loses correlation. It is the failure of a single real mixing coefficient to map $(Q,U)$ from one frequency to another, because the observed Stokes field is a vector sum of layers with different SEDs and different plane-of-sky (POS) field orientations. Three-dimensional thermal-dust information, grain-size (and more generally grain-property) knowledge, and POS magnetic-field geometry are not three independent “templates.” They are the three physical factors in that vector sum: how much dust is where, how each parcel emits and polarizes versus frequency, and in which Stokes direction it adds. Used that way, they convert an uncontrolled phenomenological residual into a constrained, testable, and marginalizable model of the dust mixing matrix.

What follows is the strongest scientifically honest version of that argument.

1. Strongest physical case

1.1 What must be true for a one-parameter dust SED to fail as a $B$-mode foreground model

For optically thin emission, the observed polarized dust field is

where $A$ is a differential dust column (or millimetre optical depth), $s_\nu(\boldsymbol\theta)$ is the intensity SED of the local grain mix, $p_\nu$ is the intrinsic polarization fraction of that mix, $\psi$ is the emission polarization angle (perpendicular to the POS field for the usual Davis–Greenstein / radiative-torque alignment geometry), $\gamma$ is the field inclination to the POS, and $g$ collects unresolved turbulent depolarization.

A single modified blackbody (MBB) per pixel, or any model of the form $P_\nu=f_\nu[\beta,T]\,P_{\nu_0}$ with real $f_\nu$, assumes four things at once:

1. $Q$ and $U$ share one frequency law;
2. $p$ and $\psi$ are frequency-independent;
3. one $(\beta,T)$ pair describes every frequency pair;
4. dust occupies one column of the component-separation mixing matrix.

Those assumptions fail if, and in practice mainly when, two conditions are satisfied simultaneously (Tassis & Pavlidou 2015; Poh & Dodelson 2017):

• the relative polarized intensities of distinct contributions change with frequency (different $T$, $\beta$, composition, size mix, alignment, or radiation field);
• those contributions do not share a common POS field orientation.

Either condition alone is milder. Different SEDs with parallel $\psi$ produce spectral curvature and an effective, frequency-pair-dependent $\beta_{\rm d}^{\rm eff}$, but $Q$ and $U$ still scale together. Misaligned fields with identical SEDs still factorize: the frequency dependence pulls out of the vector sum. The structurally dangerous case is the product of the two: the complex sum rotates,

and no real scalar maps $(Q,U)$ from $\nu_1$ to $\nu_2$. Pixel-level MBB fits of CMB Stokes amplitudes remain approximately unbiased under temperature or $\beta$ mixing without angle mixing; they become strongly biased when both the SED and $\psi$ vary along the line of sight (McBride, Bull & Hensley 2023). That is the operational theorem the three ingredients are meant to enforce.

The same physics appears at the power-spectrum level as

Template cleaning $P_{\rm clean}(\nu_{\rm CMB})=P_{\rm obs}(\nu_{\rm CMB})-\alpha P_{\rm d}(\nu_{\rm high})$ then leaves a residual of order $(1-\rho^2)$ times the dust $BB$ power at $\nu_{\rm CMB}$. Likelihoods that factorize $C_\ell^{\rm d}(\nu_i,\nu_j)=A_{\rm d}f_if_j$ underestimate the auto-spectra relative to the cross-spectra and can bias $r$ high. Both failures are mixing-matrix failures, not merely “dust is bright.”

1.2 Three distinct physical channels, not one “decorrelation parameter”

It is a mistake to treat “spectral-index decorrelation” as a single number $\Delta_{\rm d}$ or $1-\rho$. At least three mechanisms produce frequency-dependent polarized morphology, and they have different observational fingerprints:

A. Line-of-sight (LOS) vector mixing. Multiple magnetically misaligned clouds with different SEDs. Fingerprint: frequency-dependent $\psi$, different effective $\beta$ in $Q$ versus $U$, and a correlation of those effects with independent tracers of cloud multiplicity and misalignment. This is the Tassis–Pavlidou channel. Pelgrims et al. (2021) isolated it by selecting high-latitude sightlines with multiple, magnetically misaligned H I clouds and detecting enhanced 217–353 GHz angle differences relative to a control sample. A $\sim 15\%$ cloud-to-cloud variation in the polarized 353/217 ratio was enough to reproduce the signal. That is a detection of channel A in selected environments, not a sky-averaged $\rho$.

B. Plane-of-sky SED spatial variation. $\beta_{\rm d}$ and $T_{\rm d}$ vary across the map or within the beam. Fingerprint: frequency-dependent morphology without requiring $\psi(\nu)$ rotation in a single pixel; first-order moment expansions in $\beta$ and $T$ capture much of the $EE$ residual (Chluba, Hill & Abitbol 2017; Mangilli et al. 2021; Vacher et al. 2022). Ritacco et al. (2023), using SRoll2 maps, detect residuals at 100, 143 and 217 GHz after subtracting a mean polarized dust SED scaled from 353 GHz, find that intensity-based models underestimate the polarized complexity, and note that $BB$ residuals are not fully described by first-order MBB moments and are not fully correlated across frequencies. Channel B is already a LiteBIRD-level threat even without LOS angle mixing: Vacher et al. (2022) found $\sim 21\sigma_r$ biases from spatially varying $\beta$ and $T$ if those variations are omitted.

C. Intrinsic grain-population $p(\nu)$. Only a subset of grains is aligned. Silicate versus carbonaceous (or “astrodust” versus PAH) mixes, size-dependent radiative-torque alignment, porosity, and possible magnetic dipole emission all make the polarized SED differ from the total-intensity SED (Draine & Fraisse 2009; Guillet et al. 2018; Hensley & Draine 2023). Planck Intermediate Results XXII found $\beta_{\rm d}^P=1.59\pm0.02$ versus $\beta_{\rm d}^I=1.51\pm0.01$ (a $3.6\sigma$ difference at the time) and a $(21\pm6)\%$ decrease in dust polarization fraction from 353 to 70 GHz, interpretable as different polarization efficiencies among grain components. BLASTPol, by contrast, found $p$ roughly constant from $250\,\mu{\rm m}$ to $3\,{\rm mm}$ on selected molecular-cloud sightlines. Channel C is a frequency-dependent polarization fraction of a single layer, not LOS mixing. It does not rotate $\psi$ by itself, but it does mean that $T_{\rm d}$ and $\beta_{\rm d}$ fitted to $I_\nu$ must not be copied blindly onto $P_\nu$.

A phenomenological $\Delta_{\rm d}$ or a single extra mixing-matrix column conflates A–C, has the wrong frequency and $\ell$ dependence, and has the wrong sky pattern. The three-ingredient programme is the statement that A is diagnosed and weighted by 3D dust plus POS $B$, C (and the layer SEDs entering A and B) is regularized by grain physics, and B is the beam- and region-averaged remainder, which is then the natural target of a spin-2 moment expansion.

1.3 What each ingredient actually supplies — tightened

3D thermal-dust information supplies the scalar weights $A_i(\hat{\boldsymbol n})$, not $Q$ and $U$. Present “3D dust maps” (Green et al. 2019; Zhang, Green & Rix 2023; Edenhofer et al. 2024) are maps of stellar reddening. Continuum millimetre emission has no distance coordinate. The independent value of tomography is to name the layers: to distinguish a single-cloud sightline from several comparable columns, to stop one cloud’s fitted $\beta$ being attributed to the whole integral, and to provide a physically motivated truncation of the number of independent polarization components. At high Galactic latitude the relevant geometry is often simple and nearby: the Local Bubble wall at $\sim 100$–$300\,{\rm pc}$ plus a small number of intermediate-velocity clouds. That is exactly the two- or three-layer problem in which channel A is both most dangerous (comparable columns) and most solvable (few layers).

Conversion of differential reddening into millimetre emissivity is not optional bookkeeping. Optical extinction traces a UV/optical cross-section; millimetre emission traces the large-grain absorption cross-section times $B_\nu(T)$. The conversion factor varies with grain size, composition, $R_V$, radiation field $\mathcal U$, and gas phase. That conversion is the grain-physics ingredient. Without it, a 3D extinction cube is only a prior on where to put extra degrees of freedom, not a prediction of $s_\nu$.

Grain-size (more honestly: grain-property) information supplies a physical family of layer SEDs, not a unique $\beta_{\rm d}(a)$. For material $m$ and radius $a$,

and $j_\nu^P$ replaces $C_{\rm abs}$ by a polarization cross-section and multiplies by an alignment function $R_m(a)$. At CMB wavelengths the large grains that dominate the emission are typically in the Rayleigh regime, so the millimetre opacity slope depends more on dielectric function, composition, porosity and amorphous structure (including laboratory evidence that silicate $\beta$ is not a frequency-independent power law; Demyk et al. 2017) than on radius itself. Size remains indispensable through three other routes:

• equilibrium and stochastic temperatures (small grains run hotter; a size mix is a temperature mix, which is already an MBB integral);
• size-dependent alignment (radiative-torque alignment is inefficient for the smallest grains, so $p(\nu)=j_\nu^P/j_\nu^I$ is a population diagnostic);
• which material dominates $I$ versus $P$ (a two-component silicate/carbon mix, or an astrodust+PAH mix, generically makes the polarized SED differ from the intensity SED).

Fanciullo et al. (2017) already showed that size changes alone do not explain emission-to-extinction polarization behaviour; alignment and magnetic geometry must vary too. Guillet et al. (2018) showed that Planck $I$ and $P$ leave strong degeneracies among composition, shape, porosity and alignment. The honest use of grain information is therefore a prior family and a principal-component basis for layer SEDs, including the physically allowed $I$–$P$ mismatch, not a pixel-by-pixel microphysical inversion.

The POS magnetic field supplies the phase $e^{2i\psi_i}$ and the geometric depolarization, not the amplitude or the SED. A 2D field map from Planck 353 GHz is already the net LOS integral — precisely the quantity that changes under decorrelation — so it cannot by itself predict frequency-dependent $\psi$. Distance-resolved $\psi_i$ is the missing piece. It can come from:

• stellar polarization plus Gaia distances (PASIPHAE-style Bayesian layer inversion; Pelgrims et al. 2023, 2024);
• orientation of H I filaments in velocity channels (Clark, Peek & Putman 2014; Clark et al. 2015; Clark & Hensley 2019);
• a Galactic-field model only as a large-scale regular prior, never as a substitute for local turbulence.

The two-cloud theorem then becomes an observational selection: sightlines with comparable $A_i p_i\cos^2\gamma_i$ and large $\Delta\psi_i$ are the ones that must not be described by one real mixing coefficient. That is a stronger statement than “magnetic fields are important.”

1.4 The observational status, stated without overclaim

Planck Intermediate Results L reported a broad-region 217–353 GHz correlation ratio below unity. Subsequent work showed sensitivity to noise and residual systematics (Sheehy & Slosar 2018). Planck Collaboration XI (2020) found no significant broad-sky decorrelation at Planck sensitivity, with $\rho_\ell^{BB}(217,353)>0.991$ at 97.5% confidence over 71% of the sky at $50\le\ell\le160$, and a mean polarized index $\beta_{\rm d}^P=1.53\pm0.02$ for $T_{\rm d}=19.6\,{\rm K}$. BICEP/Keck BK18 likewise found $\Delta_{\rm d}$ peaking at unity in their patch, with $\Delta_{\rm d}>0.98$ (68% C.L.) in later H I-correlation work, and no significant evidence that allowing decorrelation improves the likelihood. Pelgrims et al. (2021) detected the targeted LOS effect. Ritacco et al. (2023) detected large-scale polarized SED residuals after a mean-SED subtraction.

The correct summary is: LOS frequency decorrelation is physically established and has been seen where it was predicted to be largest; it is not a measured, universal, already-calibrated coefficient at the level of current $r$ analyses; it is a next-generation residual class. That is exactly the regime in which an external physical model has leverage: Planck and BK18 do not already absorb it, and LiteBIRD / CMB-S4 / Simons Observatory will be limited by it if it is left phenomenological.

1.5 The compressed, still-physical statement

The strongest case is therefore not “3D maps plus grain sizes plus $B$ reconstruct the dust sky.” It is:

Polarized-dust frequency decorrelation is the sky- and frequency-dependent, complex (spin-2) residual of a one-column mixing matrix. Its leading physical generator is LOS superposition of a small number of magnetically misaligned layers with different grain SEDs. Three-dimensional extinction (matched to H I) identifies and weights those layers; grain-property information, including size, $R_V$, composition and alignment, supplies a low-dimensional family of polarized SEDs and the $I$–$P$ mismatch; distance-resolved POS $B$ supplies the phases. Together they determine which extra mixing-matrix columns or spin-moments are required, where they act, and what priors they carry, so that the residual dust $BB$ that leaks into $r$ can be marginalized rather than absorbed.

That is a steelman of a component-separation strategy, not of a deterministic ISM reconstruction.

2. Additional mechanisms and observational handles

These belong in a stronger essay because they are real, currently available, or imminent, and they map onto specific terms in the forward model.

2.1 H I as a 3D Stokes prior, not only a cloud counter

H I 21-cm spectroscopy (HI4PI, GALFA-H I, upcoming GASKAP / SKA pathfinders) does three separate jobs:

1. Cloud counting and velocity segmentation. Pelgrims et al. (2021) used this to select the target sample. It should be used to model every line of sight, not only to prove that an effect exists.
2. POS field orientation per velocity channel. Linear H I structures align with the field traced by starlight and by Planck 353 GHz (Clark et al. 2015). Clark & Hensley (2019) built parameter-free PPV Stokes maps from HI4PI whose velocity-integrated $Q,U$ correlate at $r>0.75$ with Planck 353 GHz and reproduce the $E/B$ asymmetry at $\ell\lesssim120$. BICEP/Keck XVI found the filamentary H I morphology template correlated with polarized dust from 95 to 353 GHz, with $\beta=1.54\pm0.13$ for that component, and no detected decorrelation inside the BICEP patch between filaments, the rest of the field, and intermediate-velocity gas. That is both a validation in a clean field and a ready-made layer-resolved $e^{2i\psi_i}$ prior.
3. LOS magnetic coherence. The velocity-space coherence of H I orientations traces LOS field tangling and therefore the geometric factor $g_i$ (Clark 2018). High $p$ at 353 GHz is a LOS-order tracer as much as an alignment tracer.

H I column also remains the best high-latitude total-intensity dust template in the diffuse ISM (Lenz, Hensley & Doré 2017), which is essential because CIB anisotropies contaminate faint 353–857 GHz intensity maps (Planck Intermediate Results XLVIII / GNILC). Polarization is largely immune to CIB; intensity-based $T_{\rm d}$ maps are not.

2.2 Stellar-polarization tomography is now a method, not a proposal

PASIPHAE’s Bayesian inversion (Pelgrims et al. 2023) infers per-layer mean POS orientation and intrinsic scatter from starlight polarization plus Gaia parallaxes, with a likelihood that carries polarization and parallax uncertainties. The first degree-scale map (Pelgrims et al. 2024) already shows that multiple dusty layers with distinct $\psi_i$ are recoverable. Scatter is not a nuisance: it is a prior on unresolved turbulent depolarization $g_i$ and, with a Davis–Chandrasekhar–Fermi or Skalidis–Tassis estimator, a handle on field strength. Optical polarization is rotated $\approx 90^\circ$ from millimetre emission polarization and weights the column by optical extinction, not millimetre emissivity. That mismatch is a feature if grain properties are modelled, and a systematic if they are not.

2.3 Inclination $\gamma$, not only POS angle

Decorrelation amplitudes depend on $p_i\cos^2\gamma_i$, not only on $\psi_i$. POS tracers do not measure $\gamma$. Extra handles:

• Faraday rotation measures (LOFAR, POSSUM, future SKA) constrain $\int n_e B_\parallel\,{\rm d}s$. Combined with H I or dust-derived path lengths they bound $\gamma_i$, especially in ionized envelopes of the same clouds.
• Total-intensity anisotropy. Hensley, Zhang & Bock (2019) detected that 353 GHz intensity per $N_{\rm H}$ depends on viewing geometry: $I/N_{\rm H}$ is smaller when $\boldsymbol B$ lies in the POS. That is the unpolarized counterpart of $\cos^2\gamma$ and a prior on inclination in the diffuse ISM.
• Polarization fraction versus column and angle dispersion ${\cal S}$. Planck’s $p$–$N_{\rm H}$ and $p$–${\cal S}$ relations already encode LOS and beam depolarization; they should be imposed as hierarchical constraints on $g_i$, not as a global $p_{\rm d}$ map.

Without $\gamma_i$, layer amplitudes $A_i p_i g_i$ remain degenerate. The steelman is not that Faraday or $I/N_{\rm H}$ measures every layer’s inclination; it is that they regularize the geometric factor that otherwise aliases into SED differences.

2.4 3D extinction curves, not only 3D extinction

Grain size and composition vary in space. Zhang & Green (2025) mapped $R(V)$ in three dimensions from Gaia XP spectra of $\sim 130$ million stars, finding systematic trends with column (accretion at low $A_V$, coagulation at higher $A_V$) and high $R(V)$ in star-forming regions. That is a direct, spatially resolved prior on ${\rm d}n/{\rm d}a$ and, more weakly, on composition. Combined with:

• Serkowski $\lambda_{\max}$ of polarized extinction (typical aligned-grain size);
• UV bump and far-UV rise (small-carbon / PAH fraction);
• FIR/submillimetre intensity SED from IRAS, Herschel, Planck 353–857 GHz, and GNILC $T_{\rm d}$ (radiation field $\mathcal U$ and large-grain temperature);
• laboratory opacity curves for amorphous silicates (Demyk et al.);
• physical dust codes (DustEM/THEMIS, astrodust+PAH: Hensley & Draine 2023),

one can build a library of layer SEDs $\{s_\nu^I,s_\nu^P,p(\nu)\}$ indexed by a few environmental parameters $(R_V,\mathcal U,$ gas phase, $A_V)$, then PCA-compress it. That is the practical meaning of “knowledge of the grain size distribution.” SPHEREx ice and dust features, and any future all-sky far-IR polarimetry, belong in the same library, not as CMB channels.

2.5 Spin-2 moment expansion as the identifiable compression

When the number of physical layers exceeds what $N_{\rm freq}$ can support, the correct compression is the spin-moment expansion of Vacher, Chluba, Aumont, Rotti & Montier (2022): moments are spin-2 objects, and the first-order moment in spectral parameters and polarization angle produces a frequency-dependent rotation in the complex plane. The draft already writes complex moments

The upgrade is to treat 3D+$B$+grain data as which moments to keep and what their complex priors are, not as an alternative to moments. Spatial $\beta,T$ variation (channel B) is $\mathrm{Re}\,{\cal M}_\beta,{\cal M}_T$; LOS angle mixing (channel A) is the imaginary, frequency-rotating part. Grain-induced $p(\nu)$ (channel C) is a mean-SED distortion of the reference $f_\nu$, not a moment of $\beta$. Conflating those three in one real $\Delta_{\rm d}(\ell)$ is how current parametric analyses both underfit and overfit.

Sponseller & Kogut (2022) showed that a moment expansion can mitigate $r$ bias from foreground SED complexity in simulations. Vacher et al. (2022) showed that, at LiteBIRD sensitivity, expansion in $\beta$ and $T$ is required to unbias $r$, at the cost of $\sigma_r$ rising because some moments correlate with the CMB. External 3D/$B$/grain priors are the way to keep the necessary moments and shrink the unnecessary ones, which is the only path that both unbiases $r$ and protects $\sigma_r$.

2.6 Other sky components that must be in the same model

A dust-only steelman that ignores the rest of the mixing matrix is incomplete:

• Synchrotron with spatially varying $\beta_{\rm s}$ and curvature, plus a dust–synchrotron spatial correlation $\epsilon$ (already in the BK18 likelihood). Low-frequency surveys (C-BASS, QUIJOTE, S-PASS, future SKA) are the analogue of 353 GHz for the other foreground.
• Polarized AME is not required by Planck Int. XXII, but a conservative analysis leaves a tightly prior-bounded residual rather than assuming zero.
• Lensing $B$ modes are a cosmological contaminant with a blackbody SED; they are not removed by dust modelling and set a floor unless delensed (CMB-S4).
• Polarized compact sources and bandpass/beam mismatch produce frequency-dependent $Q,U$ that mimic decorrelation and must be simulated with the same pipeline.
• Zodiacal / interplanetary dust is a subdominant but time-dependent large-scale term for a satellite (LiteBIRD); it is not Galactic tomography’s job, but it must not be absorbed into Galactic moments.

2.7 Diagnostics that are more powerful than $\rho_\ell^{BB}$

A stronger analysis uses:

• frequency-dependent polarization-angle maps, binned by H I-predicted misalignment (the Pelgrims test, run as a null test on the cleaned maps);
• $EE$, $BB$, $EB$, $TB$, and $TE$ of residual dust, including the Planck $E/B\simeq 2$ filamentary signature (Planck Intermediate Results XXXVIII);
• cross-spectra of CMB-cleaned maps with layer-resolved H I Stokes templates and with stellar-polarization layers;
• split tests: Local Bubble versus IVC velocity ranges; single-cloud versus multi-cloud masks;
• held-out frequencies and in-band spectral imaging if available (bolometric interferometry has been proposed as a residual detector; the physical model should make a quantitative prediction for the in-band $\psi(\nu)$ slope).

Statistical-anisotropy estimators that hunt residual filamentary dust at $r\sim 10^{-3}$ (e.g. Rotti & Huffenberger-type tools) are complementary null tests, not a replacement for the physical mixing model.

3. Why this could matter for $r$

3.1 The numerical target

LiteBIRD’s design goal is $\delta r\sim 10^{-3}$ including foregrounds and systematics; CMB-S4 aims at $r<0.001$ (95% C.L.) if $r=0$, or a $\gtrsim 5\sigma$ detection of $r=0.003$. In the BICEP2 field, Planck Intermediate Results XXX already found dust $D_\ell^{BB}\simeq 1.32\times 10^{-2}\,\mu{\rm K}^2$ at 150 GHz over $40<\ell<120$, comparable to a primordial signal with $r$ of order a few $\times 10^{-2}$. BK18 still has $A_{\rm d,353}=4.4^{+0.8}_{-0.7}\,\mu{\rm K}^2$ at $\ell=80$. Even in the cleanest southern hole, residual dust at the percent-to-few-percent level is an $r\sim 10^{-3}$ problem. On the larger LiteBIRD sky the dust is brighter, and the reionization bump ($\ell\lesssim 10$) is additionally exposed to the largest-scale SED residuals (Ritacco et al. 2023).

A template with correlation $\rho$ leaves a residual of order $(1-\rho^2)D_\ell^{\rm d}(\nu_{\rm CMB})$. Planck XI’s broad-sky 97.5% bound $\rho(217,353)>0.991$ still allows $1-\rho^2\sim 0.018$, i.e. a $\sim 2\%$ uncorrelated dust fraction between 217 and 353 GHz. Extrapolated to 100–150 GHz, and evaluated on the full LiteBIRD sky rather than a deep hole, that is not automatically negligible at $r=10^{-3}$. Pelgrims’ $\sim 15\%$ cloud-to-cloud variation in the polarized 353/217 ratio is larger than the sky-averaged bound because it is concentrated on the sightlines where channel A operates. Next-generation experiments average those sightlines into $C_\ell^{BB}$ unless they are modelled or masked.

3.2 How the bias actually enters an $r$ posterior

Three leakage paths matter, and they are not equivalent:

1. Template residual. A real scaling of a 353 GHz (or 220 GHz) map leaves a frequency-dependent dust $B$-mode that is not traced by the template. Because dust $BB$ is red in $\ell$ ($\alpha\simeq -2.4$; Planck XXX) it hits both the reionization and recombination windows. Sign: typically positive bias on $r$ if the auto-spectrum excess is interpreted as CMB.

2. Factorized-likelihood mismatch. Models $C_\ell^{\rm d}(\nu_i,\nu_j)=A_{\rm d}f_if_j$ assume $\rho=1$. If $\rho<1$, cross-spectra are smaller than the geometric mean of autos. Depending on the frequency set, this can bias $r$ up (residual auto power attributed to CMB) or, if a free $\Delta_{\rm d}$ is added, down (BK15 found $\Delta_{\rm d}$ to be a biasing parameter: in no-decorrelation simulations, 72% of $r$ posteriors peaked at zero once $\Delta_{\rm d}$ was freed). Phenomenological $\Delta_{\rm d}$ therefore both fails to describe the sky pattern of channel A and correlates with $r$.

3. Moment–CMB degeneracy. Expanding in $\beta$ and $T$ unbiases $r$ in Vacher et al. (2022) but inflates $\sigma_r$ because some moments are partly degenerate with a blackbody. Unconstrained extra parameters protect against bias by destroying statistical power. This is the quantitative reason to introduce 3D/$B$/grain information: to put a physically localized, complex prior on ${\cal M}_\beta,{\cal M}_T,{\cal M}_\psi$ so that the moments that are not generated by the ISM are not left free to absorb CMB $B$ modes.

McBride et al. (2023) make the map-level counterpart: a pixel MBB is unbiased for CMB $Q,U$ unless SED and angle vary along the LOS, in which case they found $>10\sigma$ biases in an illustrative model. That is channel A, and it is not cured by fitting a more flexible real SED.

3.3 What a successful 3D programme actually buys for $r$

It does not buy a noiseless dust map or a guaranteed $\sigma_r$ at the statistical floor. It can buy four concrete reductions of the foreground systematic:

Field and pixel triage. The two-condition theorem is a mask and a weighting scheme. Sightlines with one dominant, magnetically coherent layer (a large fraction of the high-latitude sky, including much of the BICEP patch, consistent with BK18’s $\Delta_{\rm d}\approx 1$ and BICEP/Keck XVI’s non-detection of intra-patch decorrelation) can be analysed with a one-column dust model. Sightlines with several comparable, misaligned H I clouds are either down-weighted, given extra complex moments, or used only as validation. This is more efficient than a global extra parameter that spends $\sigma_r$ everywhere.

A physically localized residual template. The Pelgrims effect has a predicted morphology: it correlates with H I-derived $\Delta\psi_i$ and with layer amplitude ratios, and it produces $\psi(\nu)$ of a definite sign and frequency slope. Subtracting or marginalizing a template with that morphology removes the part of $1-\rho$ that is actually present, rather than a Gaussian isotropic $\Delta_{\rm d}(\ell)$.

Informative moment priors. 3D dust plus grain physics predict the variance of $\beta_i$ and $T_i$ per beam (channel B) and the covariance of those fluctuations with $I_{353}$ and with $R(V)$. POS $B$ tomography predicts the complex phase of ${\cal M}_\beta$. That is the difference between Vacher et al.’s unbiased-but-degraded $\sigma_r=8.8\times 10^{-4}$ and a potentially tighter error bar that is still unbiased. In CMB-S4 pipeline-comparison exercises, allowing a free dust decorrelation parameter can remove a $\sim 10^{-3}$ bias at the cost of a tens-of-percent to factor-of-two hit in $\sigma(r)$. External priors are how one takes the bias down without taking the full $\sigma(r)$ hit.

A null-test ladder matched to $r\sim 10^{-3}$. Frequency-dependent $\psi$ in multi-cloud pixels, layer-template cross-spectra of the CMB map, and $EB/TB$ of residuals are tests whose failure modes are not degenerate with a primordial Gaussian isotropic $BB$. That is worth as much as a smaller central-value residual: it is how a claimed $r=0.003$ detection would be defended.

3.4 Honest dynamic range

On the BICEP patch, current data do not require this machinery at BK18 noise. The case for it is next-generation and wide-sky: LiteBIRD’s reionization bump, CMB-S4’s deep patch plus delensing, and any analysis that combines them. The grain and 3D data are least complete exactly where CMB experiments observe (high latitude, low $A_V$, sparse stars, weak H I). The steelman must therefore include a high-latitude prior-transfer step: calibrate layer SED families and H I–dust polarization ratios on intermediate-latitude, well-sampled sightlines, then apply them at high latitude with inflated variance, not as delta-function templates.

If that transfer is wrong, the model can bias $r$ by imposing an incorrect complex moment. The mitigation is hierarchical shrinkage (the microwave data can override the prior) plus the held-out-frequency and Pelgrims-style null tests. The approach is not safer than a phenomenological model by construction; it is safer if and only if the extra constraints are shown, in end-to-end sims that include prior misspecification, to reduce net $r$ bias at fixed coverage.

4. Pipeline sketch

A concrete sequence that a next-generation analysis could actually run:

0. Define the data vector. Multi-frequency $I,Q,U$ from the CMB experiment (LiteBIRD and/or SO/CMB-S4), plus Planck/WMAP/NPIPE for 30–353 GHz coverage; bandpasses, beams, leakage kernels, and polarization-angle calibration covariances as part of the forward model, not as afterthoughts.

1. Build a probabilistic 3D extinction cube. Use Gaia XP stellar parameters and extinctions (Zhang, Green & Rix 2023), inverted maps (Edenhofer et al. 2024; subsequent all-sky cubes), and $R(V)$ (Zhang & Green 2025). Keep posterior samples, not only the mean: high-latitude $A_V$ is comparable to the map uncertainty. Convert to a small number of distance layers per pixel by change-point or Bayesian complexity control (Local Bubble wall; one or two more distant sheets; optional IVC layer). Do not use more layers than the stellar sampling supports.

2. Match layers to H I PPV. For each pixel, associate extinction layers with H I velocity components by morphological correlation (histogram of oriented gradients / RHT) and by requiring consistency of column. Store, for each layer $i$: $N_{{\rm H},i}$, a velocity centroid, an H I-derived POS angle $\psi_i^{\rm HI}$ and coherence, and a flag for distance–velocity degeneracy. Use GNILC or H I-correlated intensity to estimate $T_i$ and millimetre optical depth, not raw 857 GHz (CIB).

3. Attach magnetic and geometric priors. Combine $\psi_i^{\rm HI}$ with stellar-polarization tomography where PASIPHAE / RoboPol / archival Heiles–Berdyugin data exist. Prior: $\psi_i\sim{\cal VN}(\psi_i^{\rm ext},\sigma_{\psi,i}^2)$ with $\sigma_{\psi,i}$ from the inversion scatter, not a global 10°. Bound $\gamma_i$ from Faraday RM where available and from $I/N_{\rm H}$ anisotropy; elsewhere use a hierarchical prior on $p_i\cos^2\gamma_i$ informed by the Planck $p$–${\cal S}$ relation. Optical versus millimetre weighting is a modelled grain-dependent offset, not a 90° rotation with identical $\psi$.

4. Attach grain SED priors. Map $(R_V,\mathcal U,$ phase, $A_V)_i$ onto a precomputed library (THEMIS / DustEM / astrodust+PAH, plus laboratory silicate opacities). Output for each layer a mean polarized SED $S_{bi}^P$ in every bandpass $b$, a mean intensity SED $S_{bi}^I$, and a small PCA basis of allowed deviations (curvature, $p(\nu)$, silicate/carbon mix). Do not fit ${\rm d}n/{\rm d}a$ per pixel. Do not force $S^P\propto S^I$.

5. Write the forward model.

convolved with the measured beam and colour-corrected with the real bandpass. $I_b^{\rm d}$ uses $S_{bi}^I$ and no $e^{2i\psi}$. $N_{\rm lay}$ is 1–3 on the high-latitude mask, higher only where 3D extinction and H I both demand it. Equivalent compressed form: a reference MBB plus complex moments $\{{\cal M}_\beta,{\cal M}_T,{\cal M}_\psi,\ldots\}$ whose amplitudes and phases have Gaussian (or shrinkage) priors taken from the layer sum above.

6. Joint hierarchical inference, not template subtraction.

Microwave data update $\{A_i,S_i,\psi_i\}$ inside their 3D/grain/$B$ priors. CMB $Q,U$ (or $E,B$) and $r$ are sampled jointly. Two analysis products, both needed:

• a posterior-mean dust map and covariance for map-level science and null tests;
• a foreground-marginalized $r$ likelihood (the BK-lite idea, but with a physically structured rather than 7-parameter phenomenological dust model).

Dedusting in the sense of subtracting a fixed template is a fallback for jackknives, not the baseline.

7. Identifiability control. If $N_{\rm freq}$ cannot support $N_{\rm lay}$ independent $(T,\beta,p,\psi)$, freeze SED shapes to the grain library, freeze $\psi_i$ to the magnetic tomography (with the quoted $\sigma_\psi$), and leave only a few amplitudes free. That is the point of the external data: they are what make the extra physical degrees of freedom not free. Run a nested sequence: 1-MBB; 1-MBB + $\mathrm{Re}\,{\cal M}_\beta,{\cal M}_T$; + complex ${\cal M}_\psi$ only in multi-cloud pixels; full 2–3 layer model on a subset. Quote $r$ for each, with a predefined criterion for which is baseline.

8. Validation and stress tests, all of which must pass before $r$ is reported.

• Held-out frequency (e.g. train without 217 or without 353, predict it).
• Pelgrims test on the CMB map: residual $\Delta\psi(217,353)$ must not correlate with H I misalignment.
• Cross-spectra with layer-resolved H I Stokes and stellar-pol layers consistent with noise.
• $EB$, $TB$ of residuals consistent with zero after angle calibration.
• End-to-end simulations with: (i) the fiducial 3D model; (ii) a wrong grain library; (iii) missing distant layers; (iv) realistic bandpass, leakage, and angle errors; (v) MHD or PySM skies (d12-class) that were not generated from the same layer model. The figure of merit is bias on $r$ and coverage of the $r$ interval, not $\chi^2$ of dust maps.
• Explicit comparison to a phenomenological $\Delta_{\rm d}$ and to a $\beta+T$ moment expansion without 3D priors, to measure the $\sigma_r$ cost and bias reduction.

9. Publish the intermediate products. Layer amplitudes, $\psi_i$ maps, moment maps, and the multi-cloud mask are scientific deliverables. They are how the community tests whether the $r$ result is foreground-limited in a physically specified way.

5. What to keep from the draft

The draft is already stronger than a typical “more data will help” essay. The following pieces should survive almost unchanged.

1. The operational definition of decorrelation. Maps decorrelate; a scalar $\beta_{\rm d}$ does not “lose correlation.” The four assumptions of a real per-pixel mixing law, and the statement that dust cleaning fails when the mixing law is too restrictive, are exactly right.

2. The two-condition theorem. Different SEDs without misalignment $\Rightarrow$ curvature, no $\psi(\nu)$; misaligned fields with identical SEDs $\Rightarrow$ still one frequency law; both $\Rightarrow$ rotation of $P_\nu$. This is the load-bearing lemma. Tassis & Pavlidou (2015) and Poh & Dodelson (2017) are the correct citations.

3. The observational nuance. Planck Int. L versus Sheehy & Slosar versus Planck XI versus Pelgrims et al. (2021), including the $\sim 15\%$ polarized 353/217 variation, is the right status report. Do not replace it with “decorrelation has been measured at $X\%$.”

4. The three-way split of ingredients. 3D maps $\to$ weights and layer count; grains $\to$ SED family and $I$–$P$ mismatch; POS $B$ $\to$ phase and geometric depolarization. The explicit statement that none of the three is sufficient alone is essential and should be kept.

5. The qualification that 3D maps are extinction, not millimetre tomography. This is the most important honesty in the draft. Keep it in the main argument, not only in a limitations list.

6. The Rayleigh-regime warning on grain size. Size does not uniquely determine millimetre $\beta$; it matters through heating and alignment. Fanciullo et al. (2017) and Guillet et al. (2018) are the right brakes on grain-model enthusiasm. Draine & Fraisse (2009) is the right citation for polarized versus total SEDs.

7. The layer forward model $P_b^{\rm d}=\sum_i A_i S_{bi}^P g_i e^{2i\psi_i}$ and the insistence on hierarchical priors rather than an immutable template. Joint marginalization with $r$ is the correct statistical endpoint.

8. The complex moment expansion as a compression, with ${\cal M}_\beta$ carrying $e^{2i\psi_i}$. Mangilli et al. (2021) and Vacher et al. (2022) belong there. This is how the idea remains identifiable.

9. The limitations list, especially incomplete high-latitude tomography, extinction-to-emission conversion, grain degeneracies, POS-only magnetic tracers, parameter identifiability versus $N_{\rm freq}$, instrumental mimics of $\psi(\nu)$, and other sky components. A steelman that dropped these would be weaker, not stronger.

10. The closing epigram, which is the right one: the ingredients do not eliminate foreground uncertainty; they replace an uncontrolled phenomenological decorrelation parameter with a physically structured vector model whose degrees of freedom can be constrained, validated, and marginalized.

Keep also the Planck mean SED numbers ($T_{\rm d}\simeq 19.6\,{\rm K}$, $\beta_{\rm d}^P\simeq 1.53$), the MBB formulae, and the $\rho_\ell^{BB}$ definition. They are accurate and they prevent the essay from floating.

6. Upgrade list

Numbered improvements the primary draft should adopt.

1. Retitle the problem as a mixing-matrix failure with three channels (LOS vector mixing, POS SED variation, intrinsic $p(\nu)$), not as “spectral-index decorrelation.” Reserve $\rho_\ell^{BB}$ and $\Delta_{\rm d}$ for operational summaries. State the two-condition theorem as the necessary condition for channel A, and McBride, Bull & Hensley (2023) as the pixel-level demonstration that CMB $Q,U$ stay unbiased unless both SED and angle vary.

2. Correct the Sponseller & Kogut citation. The physics sentence about structural failure when no common real scaling of $Q$ and $U$ exists points to McBride et al., MNRAS 519, 4370 (2023), not to Sponseller & Kogut. Sponseller & Kogut (2022, ApJ 936, 8) is the moment-expansion paper on mitigating $r$ bias and should be cited for that.

3. Quantify the $r$ target. LiteBIRD $\delta r\sim 10^{-3}$; CMB-S4 $r<0.001$ (95%) if $r=0$. Quote Planck XXX’s $D_\ell^{BB}(150\,{\rm GHz})\simeq 1.32\times 10^{-2}\,\mu{\rm K}^2$ in the BICEP2 field, BK18’s $A_{\rm d,353}$ and $\Delta_{\rm d}$ tests, Planck XI’s $\rho>0.991$ (97.5%), and Vacher et al.’s $\sim 21\sigma_r$ bias from unmodelled $\beta,T$ variation plus $\sigma_r=8.8\times 10^{-4}$ after $\beta+T$ moments. Include the $(1-\rho^2)D_\ell^{\rm d}$ residual scaling and the fact that a free $\Delta_{\rm d}$ is a biasing parameter (BK15). Without these numbers the essay cannot claim relevance for $r$.

4. Add Ritacco et al. (2023) as the current large-scale polarized-SED residual result: mean-SED subtraction leaves 100/143/217 GHz residuals; intensity-based models underestimate polarized complexity; first-order moments fit $EE$ better than $BB$; residuals are not fully correlated across frequencies. That is the observational bridge from Planck XI’s “no broad-sky $\rho<1$” to “next-generation analyses still need extra polarized degrees of freedom.”

5. Promote H I from a parenthetical tracer to a first-class Stokes prior. Cite Clark et al. (2015) and Clark & Hensley (2019) for PPV Stokes maps, and BICEP/Keck XVI for a 95–353 GHz correlation with $\beta=1.54\pm0.13$ and a non-detection of intra-patch decorrelation. The pipeline should fit using velocity-resolved H I $Q,U$, not only select multi-cloud pixels as Pelgrims et al. did.

6. Add 3D extinction-curve maps (Zhang & Green 2025) as the actual grain-size observable. Pair them with Serkowski $\lambda_{\max}$, FIR $T_{\rm d}$, and a DustEM/THEMIS/astrodust library. State explicitly that millimetre $\beta$ is not a function of $a$ in the Rayleigh regime, so the library is a prior on $(T,p(\nu),I/P$ mismatch, curvature), reduced to a few PCA components.

7. Add inclination handles: Faraday RM (LOFAR/POSSUM/SKA) for $B_\parallel$; Hensley, Zhang & Bock (2019) $I/N_{\rm H}$ viewing-geometry effect; Planck $p$–${\cal S}$ as a hierarchical prior on $g_i$. Without $\gamma$, $A_ip_i\cos^2\gamma_i$ aliases into SED differences.

8. Upgrade moments to spin-2 moments (Vacher et al. 2022, A&A, spin-moment paper). Identify $\mathrm{Re}\,{\cal M}_{\beta,T}$ with channel B and the frequency-rotating complex part with channel A. Use 3D/$B$/grain data to choose which moments exist and to set complex priors, with the explicit goal of unbiasing $r$ without paying the full unconstrained-moment penalty in $\sigma_r$.

9. Write field triage as a science product. A multi-cloud / misaligned-$\psi$ mask, and a complementary single-layer mask, should be delivered and used both for the baseline $r$ likelihood and for stress tests. This is the cheapest, most robust use of the three ingredients and is already supported by Pelgrims et al. (2021) plus BICEP/Keck XVI.

10. Specify high-latitude prior transfer. Calibrate SED families and H I–dust polarization ratios at intermediate latitude, apply at high latitude with extra hierarchical variance. Do not pretend Gaia tomography is uniformly good in the CMB fields.

11. Put synchrotron, $\epsilon_{\rm dust-sync}$, lensing $B$, compact sources, and bandpass/beam/angle calibration in the same forward model. A dust-only mixing matrix is not a CMB likelihood. Instrumental $\psi(\nu)$ and bandpass mismatch are the leading false positives for channel A.

12. Replace “subtract the posterior mean” as the headline with joint marginalization of dust and $r$. Keep mean-map subtraction only for null tests. Require end-to-end simulations that include wrong grain libraries and missing layers, with coverage of $r$ as the acceptance test.

13. Add a diagnostic section beyond $\rho_\ell^{BB}$: $\psi(\nu)$ versus H I misalignment on cleaned maps; cross-spectra with layer-resolved templates; $EB/TB$; Local Bubble versus IVC splits; held-out frequencies. A detection of $r$ at $10^{-3}$ that fails these tests is not a detection.

14. Name the current surveys in the pipeline, not only in citations: Gaia XP + Edenhofer-scale 3D dust; HI4PI / GALFA-H I; Planck 353 and NPIPE; PASIPHAE / RoboPol; C-BASS, QUIJOTE; LOFAR RM; IRAS/Herschel/GNILC $T_{\rm d}$; BK / SO / CMB-S4 / LiteBIRD frequency sets. An implementable argument is more convincing than an abstract one.

15. Keep the limitations as first-class caveats in §1, not only at the end. In particular: extinction $\neq$ millimetre emissivity; Rayleigh-regime $\beta(a)$; POS-only $B$; $N_{\rm freq}$ versus $N_{\rm lay}$; the risk that an incorrect physical prior biases $r$ more cleanly than a free $\Delta_{\rm d}$. The steelman is that a hierarchical, validated, moment-compressed layer model is the least bad way to spend extra degrees of freedom at $r\sim 10^{-3}$, not that it eliminates foreground uncertainty.

Adopt those fifteen items and the draft’s existing physics becomes a complete argument: not that 3D dust, grain sizes and projected $B$ “help dedusting” in the abstract, but that they are the physically indexed, observationally available priors on the complex, layer-sum mixing matrix whose mismodelling is the concrete systematic between current $r<0.036$ (BK18) and a credible $r\sim 10^{-3}$ measurement.