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

Steelman: 3D dust emission, grain-size physics, and the plane-of-sky field as a structured prior for CMB B-mode dedusting

The useful claim is not that three extra datasets “solve dust.” It is that frequency decorrelation of polarized dust is a vector superposition problem with a small number of physically distinct layers, and that the three ingredients named in the question supply the three quantities that CMB data alone cannot pin down per pixel: layer weights and temperatures, the intensity-to-polarization transfer, and the relative Stokes phases. Once those are supplied externally, the leftover CMB-facing problem is almost linear, and the residual that a one-parameter modified blackbody (MBB) leaves at $r\sim10^{-3}$ becomes a predicted, spatially localized correction rather than an unconstrained systematic.

That is the argument worth keeping. What follows tightens it.

1. Strongest physical case

1.1 What “spectral-index decorrelation” actually is

Parametric component separation (Commander, FGBuster, BICEP/Keck-style multi-frequency likelihoods) and 353 GHz template extrapolation all assume that, in each pixel or each patch, the polarized dust sky factors as

with the same $(\beta_d,T_d)$ for $Q$ and $U$. That implies a frequency-independent polarization angle and a polarized SED that is a single MBB. “Decorrelation of the spectral index” is a slightly misleading name for the failure of this factorization.

The physical sky in a beam is a sum over emitting regions $i$ (distinct clouds along the line of sight, and sub-beam structure):

Each $\mathcal{P}_i$ is a spin-2 quantity whose phase is set by the plane-of-sky field angle $\chi_i$ of that region (emission polarized perpendicular to $\mathbf{B}_\perp$, so $\psi_i=\chi_i+\pi/2$). Four consequences follow, and they are not the same effect:

1. The effective SED is not an MBB. $|\mathcal{P}_\nu|$ is a vector sum. The pairwise index $\beta_{\rm eff}(\nu_1,\nu_2)$ depends on the band pair and can even lie outside the range of the $\beta_i$ when layers partially cancel. This is the object of the moment expansion (Chluba, Hill & Abitbol 2017) and of its polarized “spin-moment” extension (Vacher et al. 2023).

2. The angle rotates with frequency. $\psi_{\rm eff}(\nu)=\tfrac12\arg\sum_i\mathcal{P}_i\hat g_i(\nu)$ varies whenever both the $\psi_i$ and the $\hat g_i$ differ. Tassis & Pavlidou (2015) isolated the two necessary conditions: the ratio of the layers’ polarized intensities must change between the frequencies, and their fields must be misaligned. No scalar SED, however curved, represents a rotation. It is a spin-2 effect.

3. $\beta_{\rm eff}^P\neq\beta_{\rm eff}^I$ from geometry alone. Intensities add as scalars, polarizations as vectors. Even identical grain populations, viewed through misaligned fields, produce a polarized SED different from the intensity SED.

4. Maps at two frequencies are imperfectly correlated. With

a Gaussian scatter $\sigma_\beta$ in a spatially varying index gives $\mathcal{R}=\exp[-\tfrac12\sigma_\beta^2\ln^2(\nu_1/\nu_2)]$. Two uncorrelated layers with power ratio $f$ whose relative weight changes by a factor $\rho$ between the bands give $\mathcal{R}=(1+\rho f)/\sqrt{(1+f)(1+\rho^2 f)}$. After scaling a 353 GHz template to 150 GHz with the best single scalar, the residual power fraction is $1-\mathcal{R}^2\approx 2(1-\mathcal{R})$ for $\mathcal{R}$ near unity.

Effects (1) and (4) can be produced by either line-of-sight mixing or plane-of-sky $\beta$ variation. Effect (2) is diagnostic of line-of-sight mixing with misaligned fields. That distinction is the reason 3D information plus $\mathbf{B}_\perp$ is not interchangeable with a free $\beta(\hat n)$ map.

Pelgrims et al. (2021) detected (2) directly: sightlines with several HI clouds whose fibre orientations are misaligned by $\gtrsim 60^\circ$ show larger 217–353 GHz polarization-angle differences than control sightlines. An inter-cloud variation of $\sim 15\%$ in the 353/217 polarized-intensity ratio is enough to reproduce the measured effect. Planck 2018 XI, by contrast, found no robust power-spectrum detection of $\mathcal{R}^{BB}_\ell(217\times353)<1$, and judged its limit adequate for $r\simeq 0.01$. Those two results sit together: the effect is real, spatially sparse, and currently below the Planck $C_\ell$ noise except where one has an independent way to select the threatened pixels. That is exactly the selection that 3D maps and $\mathbf{B}_\perp$ provide.

1.2 Why CMB data alone cannot absorb the extra freedom

Adding layers or moments as free parameters per pixel does not automatically remove the bias. McBride, Bull & Hensley (2023) showed that a single line of sight constrains only the parameters of broad distributions of $T$ and $\beta$, and that two-parameter MBB fits are biased at $>10\sigma$ in illustrative cases precisely when both the SED and the angle vary along the sightline — i.e., in the Tassis–Pavlidou regime. Vacher et al. (2022) found that free $\beta$ and $T$ moments recover an unbiased $r$ for LiteBIRD only at $\sigma(r)\approx 8.8\times 10^{-4}$, because the moments are degenerate with $r$. Hensley & Bull (2018) and Sponseller & Kogut (2022) reached the same qualitative conclusion: extra spectral freedom costs sensitivity and does not break the degeneracy.

The missing information is therefore structural: which layers exist, how they are weighted, what SED each has, and how their Stokes vectors are oriented. CMB frequency maps measure only the sum.

1.3 How the three ingredients map onto that structure

A single emitting region with column corresponding to optical depth $\tau_{\nu_0}$, dust temperature $T_d$, opacity index $\beta_d$, plane-of-sky field angle $\chi$, and inclination $\gamma$ of $\mathbf{B}$ to the plane of the sky produces

where $p_0$ is the intrinsic polarization of the aligned grain population and $F\le 1$ accounts for imperfect alignment and unresolved field tangling. The three external ingredients attach to different factors in this expression, and to the sum over regions.

3D maps of thermal dust emission supply $N(\hat n)$, the layer columns $w_i$, and — once the 3D extinction is fitted to FIR intensities — the per-layer intensity SED $(T_i^I,\beta_i^I)$. “3D emission map” here means a distance-resolved emissivity, not just 3D extinction:

Edenhofer et al. (2024) map extinction to 1.25 kpc at up to 14′ in 516 distance bins from 54 million Gaia XP stars. Zelko et al. (2022) took the step to emission, proposing $T$ and $\beta$ for Bayestar19 voxels and matching Planck and IRAS bands; they resolved two distinct temperatures along the Cepheus sightline and found $\tau_{353}/E(B-V)$ varying by a factor of two. Complementary 3D structure comes from HI velocity: Panopoulou & Lenz (2020) find typically 2–3 Gaussian HI components per 0.2 deg² pixel at high latitude. At $|b|\gtrsim 60^\circ$, the Local Bubble wall at $\sim 70$–300 pc often carries most of the dust column (Pelgrims et al. 2020; O’Neill et al. 2024). The high-latitude B-mode problem is therefore frequently a two-layer problem — Bubble wall plus a subdominant distant screen — not an integral over a continuum of clouds.

On its own, the 3D emission map does five things that are already usable:

• Risk maps. $N$ and $w_i$ identify where a single MBB is adequate. Pelgrims et al. (2021) used the HI version of this as a target/control split and detected the effect. This is a mask and field-selection tool, available today.
• Layer SEDs from intensity. Intensity is measured at far higher S/N than polarization. If $(T_i^I,\beta_i^I)$ are taken from the 3D emission fit, the polarization problem reduces to $2N$ real amplitudes per pixel plus CMB, linear once the SED shapes are held fixed.
• Layer templates. Per-layer intensity (or HI) maps can be fitted against $Q$ and $U$ with regional coefficients, as Ghosh et al. (2017), Adak et al. (2020), and BICEP/Keck XVI already do with velocity layers.
• Physical $T$ structure. Distance from heating sources, and the Local Bubble geometry, predict where $T$ differs between layers. Guillet et al. (2026) find that SRoll2 residual covariances around the mean polarized SED are compatible with temperature-only variations (PR4 is compatible with neither temperature-only nor index-only, a warning about Planck systematics). If that SRoll2 result holds, the 3D $T$ field is the dominant SED lever.
• Simulations with real geometry. PySM d12, DustFilaments, and spin-moment skies replace ad-hoc layer statistics with the actual 3D field, which is how $\Delta_d(\ell;\nu_1,\nu_2)$ priors and pipeline validation should be built.

The 3D map does not say which layers polarize, what $\beta^P-\beta^I$ is, or what the $\psi_i$ are.

Grain-size distribution and alignment supply the intensity-to-polarization transfer. In the Rayleigh limit relevant at millimetre wavelengths, the opacity per unit mass is independent of grain radius, so $\beta$ is a property of the dielectric function (material, porosity, two-level systems, iron inclusions), not of size. Equilibrium temperature, by contrast, does depend on size: for grains large enough that optical absorption is geometric while FIR emission is still Rayleigh, $T\propto a^{-1/(4+\beta)}$, so a factor of two in radius shifts $T$ by $\sim 2$ K at 20 K. Grains smaller than $\sim 0.01$–$0.02\,\mu\mathrm{m}$ are stochastically heated and contribute little below 353 GHz. Radiative-torque alignment (Lazarian & Hoang 2007; Hoang & Lazarian 2016; review by Andersson, Lazarian & Vaillancourt 2015) sets a threshold $a_{\rm align}\sim 0.05$–$0.1\,\mu\mathrm{m}$ that depends on the radiation field $U$, gas density and temperature, and grain shape. Only the large end of $n(a)$ polarizes:

Thus $\beta^P-\beta^I$ and $T^P-T^I$ are predictions of a size distribution, an alignment function, and a composition, not free pixel parameters. Two further grain-physics facts are load-bearing:

• Which layers polarize. A dense, shielded layer with $a_{\rm align}$ pushed to large $a$ can dominate $I$ and not $\mathcal{P}$. That is a geometric-looking $\beta^P\neq\beta^I$ that is actually an alignment effect, and it is invisible to a 3D extinction map.
• The optical-to-submillimetre bridge. Starlight polarization is dichroic extinction by grains near $a\sim 0.1$–$0.2\,\mu\mathrm{m}$ at $\lambda_{\max}$. 353 GHz polarized emission is thermal emission from the same aligned large-grain population. The ratio $R_{P/p}\equiv P_{353}/p_V\approx 5.4\,\mathrm{MJy\,sr}^{-1}$ measured by Planck 2018 XII in the diffuse ISM is therefore a conversion from per-layer starlight $(q_i,u_i)$ to per-layer 353 GHz $(Q_i,U_i)$. In that transfer, $\cos^2\gamma_i$ and $F_i$ cancel. That is the correction Tassis & Pavlidou (2015) proposed, and it is the reason starlight tomography is not merely an angle prior.

Competing grain models are not a reason to discard this lever. Draine & Fraisse (2009) predicted a rising $p(\nu)$ if only silicates align; Hensley & Draine (2023) argue that a single “astrodust” composite for $a\gtrsim 0.02\,\mu\mathrm{m}$ most naturally explains the observed near-constancy of the FIR polarization fraction from 250 $\mu\mathrm{m}$ to 3 mm (Planck 2018 XI; BLASTPol). Planck intermediate XXII found $\beta^P=1.59\pm 0.02$ versus $\beta^I=1.51\pm 0.01$ and a $21\pm 6\%$ drop in $p$ from 353 to 70 GHz; Planck 2018 XI prefers a single MBB with $\beta_d^P=1.53\pm 0.02$ at $T_d=19.6$ K. The discrepancy is a prior of finite width, not a wall. A prior of width $\sim 0.05$–$0.1$ on $\beta^P-\beta^I$, and a physically bounded $R_{P/p}(R_V,U)$, is already far tighter than what polarization data alone provide per pixel.

The sky-projected magnetic field supplies the phases. Misalignment between layers is what turns SED differences into angle rotation and vector cancellation. Planck 353 GHz gives only the emission-weighted mean $\chi$. Per-layer $\chi_i$ come from four sources, in decreasing directness:

• Starlight polarization of stars at known Gaia distances: Panopoulou et al. (2019) on a diffuse sightline; BISP-1 Bayesian inversion (Pelgrims et al. 2023); the first degree-scale tomography map over 4 deg² with RoboPol (Pelgrims et al. 2024); PASIPHAE/WALOP as the high-latitude survey built for this purpose (Tassis et al. 2018).
• HI fibres: Clark et al. (2015); Clark & Hensley (2019) built PPV Stokes maps from HI4PI that correlate at $r>0.75$ with Planck 353 GHz $Q,U$ with no free parameters and reproduce $E/B$ at $\ell\lesssim 120$. BICEP/Keck XVI detected dust polarization correlated with that morphology template down to 95 GHz, with $\beta=1.54\pm 0.13$, and found no decorrelation between filamentary and non-filamentary dust in the BK field.
• Local Bubble field models (Alves et al. 2018; Pelgrims et al. 2020; O’Neill et al. 2025) for the dominant high-latitude layer.
• Faraday rotation measures (Hutschenreuter et al. 2022 and successors) for $B_\parallel$, hence a weak prior on $\gamma$. This last handle must be used as a prior, not a measurement of dust $\gamma$: RM weights $n_e B_\parallel$, dust weights $n_d$.

The phase information does three things that amplitudes cannot: it identifies where the vector-sum problem exists; it converts starlight $(q_i,u_i)$ into 353 GHz $(Q_i,U_i)$ via $R_{P/p}$; and it defines the pivot angle for a spin-moment expansion, so that the first-order complex moment’s phase relative to $\mathcal{W}_0=\sum_i\mathcal{P}_i$ is the frequency-dependent rotation.

1.4 The combined statement, tightened

The strongest version of the argument is therefore:

Polarized dust at CMB frequencies is a sum of a few discrete magnetized screens. 3D emission maps fix the screens’ columns and intensity SEDs. Grain-size and alignment physics convert those into polarized SEDs and into a starlight-to-submm transfer, and tell you which screens polarize. Plane-of-sky field measurements fix the relative Stokes phases. The CMB-facing model then has the correct functional form for both SED curvature and angle rotation. The remaining free parameters are small regional amplitude and angle corrections, not a new spectral index per pixel. That is what “dedusting” can mean at $r\sim 10^{-3}$.

Two important boundings keep this from overclaiming. First, a subdominant layer is enough: the first row of a two-layer calculation with layer 2 carrying 30% of the 353 GHz polarized intensity, $\Delta\beta=0.3$, and fields misaligned by 45° already rotates $\psi_{\rm eff}$ by $\sim 2^\circ$ between 353 and 150 GHz. Second, most high-latitude sightlines are not in this regime. The 3D data’s first job is to say which pixels are which.

2. Additional mechanisms and observational handles

These are real extra levers that belong in a stronger essay. None of them is speculative physics; several are already in use.

Line-of-sight versus plane-of-sky decorrelation. Both produce $\mathcal{R}<1$. Only LOS mixing with misaligned $\mathbf{B}$ produces frequency-dependent $\psi$, frequency-dependent $E/B$, and frequency-dependent dust $EB$ (Vacher et al. 2023 on spin-moments and $E/B$, $EB$). Ritacco et al. (2023) found that first-order MBB moments fit SRoll2 residual $EE$ but not $BB$, and identified angle rotation as the missing piece. A 3D+$\mathbf{B}_\perp$ model predicts all three; a spatially varying $\beta(\hat n)$ predicts only $\mathcal{R}<1$. That is a discriminant, not just a prior.

The Local Bubble as the default high-latitude two-layer geometry. O’Neill et al. (2024) reconstruct the dust-traced Bubble surface from the Edenhofer map: mean distance $\sim 170$ pc, thickness $\sim 35$ pc, a northern “chimney.” O’Neill et al. (2025) argue that the Bubble dominates the polarized signal on a majority of sightlines and reconstruct a 3D field assuming the field is tangent to the shell. Halal, Clark & Tahani (2024) did not find an imprint of wall inclination on $p_{353}$, but did find that splitting the column among multiple distance components within 1.25 kpc depolarizes. The honest synthesis is: the Bubble is often the dominant screen, the distant residual is the decorrelation threat, and 3D maps measure that residual’s weight. For recombination-bump experiments on clean southern patches, this is the right cartoon. For LiteBIRD’s reionization bump at $\ell<20$ over large $f_{\rm sky}$, intermediate-latitude sightlines have more distant disk layers, and the 1.25 kpc cutoff is a real limitation.

3D radiation field as a predictor of $T_i$, not just a fit. Dust temperature is set by the local interstellar radiation field $U$. 3D stellar catalogues plus 3D extinction predict $U(\mathbf{x})$, hence $T(a,\text{composition})$ and $a_{\rm align}(U,n_{\rm H})$. That closes a loop the draft leaves open: Zelko-style $T_i$ need not be fitted independently in every voxel if the heating is predicted. It also tells grain physics which layers should be poorly aligned.

3D $R(V)$ as a size-distribution map. Zhang & Green (2025) mapped $R(V)$ in three dimensions from 130 million Gaia XP spectra, finding variation with extinction consistent with accretion at low $A_V$ and coagulation at higher $A_V$. $R_V\approx 5.5\,\lambda_{\max}/\mu\mathrm{m}$ (Whittet & van Breda) ties $R(V)$ to the aligned-grain size that produces starlight polarization. This does not determine millimetre $\beta$ (a dielectric-function quantity), but it does modulate $R_{P/p}$ and $T^P-T^I$. At the low columns of B-mode fields one is in the accretion/diffuse regime, where a one-parameter $R(V)$ prior is most defensible.

HI PPV Stokes as the all-sky $\chi_i$ available now. Starlight tomography is sparse. Clark & Hensley (2019) already give a full-sky, velocity-resolved $\mathbf{B}_\perp$ tracer that matches Planck 353 GHz morphology. Combined with a velocity–distance assignment (3D extinction as the distance prior on HI components), this is the practical all-sky layer-angle map until PASIPHAE is complete. BICEP/Keck XVI is the existence proof that this template correlates with dust $B$ modes down to 95 GHz in a cosmological field.

Submillimetre polarimetry near the SED peak. BLASTPol, PILOT, JCMT POL-2, NIKA2-Pol, and any future balloon or CCAT-prime/FYST polarized continuum sit at 250–850 $\mu\mathrm{m}$, where temperature differences between layers are maximally expressed in the SED. Millimetre CMB bands are on the Rayleigh–Jeans side, where $\Delta T$ is partly degenerate with $\Delta\beta$. High-frequency polarized intensity is the empirical $T_i^P$.

CMB-experiment high-frequency channels as dust, not just foreground monitors. Fuskeland et al. (2023) found that extending LiteBIRD from 402 to 500–600 GHz reduces $\sigma(r)$ by 30–50% when $\beta$ and $T$ vary on the sky. SO/S4 280 GHz and LiteBIRD 337/402 GHz are $T^P$ measurements. If 3D maps predict that temperature variation dominates (Guillet et al. 2026 on SRoll2), this is where the extra bands pay off; if magnetic-dipole emission or index variation dominates, the 100–220 GHz slope is the lever instead. Grain physics plus 3D $T$ tell you which.

Molecular-gas tracers for layers HI misses. HI fibres fail in molecular and ionized gas. CO, OH, CH, and dust in emission at 850 $\mu\mathrm{m}$ pick up the denser screens. At high latitude these are rare, which is why the method works first there; at intermediate latitude they are required.

CIB-cleaned intensity for 3D emission fits. Zelko-style $T,\beta$ fits use Planck/IRAS intensities that contain CIB and zodiacal light. GNILC-style CIB deprojection, and a 3D model that is not asked to fit the CIB monopole, are part of making $\epsilon_\nu$ a dust quantity. Otherwise the 3D “emission” map is partly extragalactic.

MHD plus dust-polarization simulations as a prior on unresolved tangling. Within a named layer, turbulent tangling and sub-beam $(\beta,T)$ variation are not mapped. They produce residual high-$\ell$ decorrelation and the observed $p$–$S$ anticorrelation (Planck 2018 XII). Filament models (Huffenberger, Rotti & Collins 2020; Clark et al. 2021; Hervías-Caimapo & Huffenberger 2022; Konstantinou et al. 2022) already explain dust $E/B\approx 2$ and positive $TE$, $TB$ from the coupling of density filaments to $\mathbf{B}$. The 3D model should retain a statistical $F_i$ (or a polarization-angle dispersion $S_i$) per layer, informed by MHD, not pretend that $\psi_i$ is a single number down to the beam.

Faraday RM as a weak $\gamma$ prior, not a dust inclination map. POSSUM, GMIMS, and LOFAR RM grids constrain $B_\parallel$ weighted by thermal electrons. Combined with starlight $\mathbf{B}_\perp$ they give a 3D field direction in ionized gas. That is useful where the Local Bubble wall is the same structure seen in dust and in the local hot/warm ionized medium, and dangerous elsewhere. Keep it in the prior; do not put $\cos^2\gamma_i$ on the RM.

SPHEREx ice and PAH spectroscopy as a composition prior. SPHEREx’s 0.75–5 $\mu\mathrm{m}$ all-sky spectra map ices and PAHs along millions of sightlines. Ice mantles and PAH growth are grain-evolution diagnostics (accretion, coagulation) that correlate with $R(V)$ and with millimetre emissivity per $A_V$. This is most relevant at $A_V\gtrsim 1$, so it is a secondary handle for the cleanest B-mode patches and a primary one for intermediate-latitude LiteBIRD sky.

Predicted spin-moment maps as the interface to existing component-separation machinery. Writing

with $\mathcal{W}_0=\sum_i\mathcal{P}_i$ and $\mathcal{W}_1^\beta=\sum_i\mathcal{P}_i(\beta_i-\bar\beta)$, the first-order spin-moments are complex and their phase relative to $\mathcal{W}_0$ is the angle rotation. A 3D+$\mathbf{B}_\perp$+grain model predicts $\mathcal{W}_1^\beta$ and $\mathcal{W}_1^T$. Fitting a few global amplitudes of those predicted templates, or deprojecting their SEDs in a constrained moment ILC (Remazeilles, Rotti & Chluba 2021), keeps the correct functional form while breaking the moment–$r$ degeneracy that free moments suffer. The delta-map method (Ichiki et al. 2019) and power-spectrum moments (Mangilli et al. 2021; Azzoni et al. 2021) admit the same treatment. This is the practical contact point with FGBuster, NILC/cMILC, and BK-style likelihoods.

Dust $EB$ and cosmic birefringence. Frequency-dependent dust angle rotation produces a frequency-dependent dust $EB$. Isotropic cosmic birefringence produces an achromatic $EB$. A 3D model that predicts dust $\psi(\nu)$ is therefore a birefringence systematic control, not only an $r$ control (Clark et al. 2021; Ritacco et al. 2023; Diego-Palazuelos et al.). This is a second science payoff of the same ingredients.

Background-galaxy polarimetry. Starlight tomography stops at the last usable star. Optical/NIR polarimetry of background galaxies (and of globular clusters at known distance) constrains the total column’s integrated $(q,u)$, which is the sum of all layers including those beyond 1.25 kpc. Sparse, but the right constraint for LiteBIRD large-scale modes.

3. Why this could matter for $r$

The scaling is arithmetic, not rhetoric.

In the cleanest few per cent of sky used by ground-based recombination-bump experiments, polarized dust $BB$ at 150 GHz and $\ell\approx 80$ is of order the primordial signal for $r\sim 0.05$–$0.1$ (Planck XXX / Planck 2018 XI; the BICEP field is in this range). LiteBIRD’s requirement is a total $\delta r<10^{-3}$ (Lv1.01), with statistical, systematic, and margin each allocated $\sim 0.6\times 10^{-3}$; CMB-S4’s design goal is $\sigma(r)\sim 5\times 10^{-4}$ after delensing. A residual of one per cent of the 150 GHz dust $BB$ in a patch whose dust is equivalent to $r=0.1$ is $\delta r\sim 10^{-3}$. That is the relevant unit.

Power-spectrum decorrelation. A two-layer sky in which the second layer carries 30% of the 353 GHz polarized amplitude (so $f\approx(0.3/0.7)^2\approx 0.18$) and $\Delta\beta=0.3$ produces $\mathcal{R}(217\times 353)\approx 0.998$–$0.999$, invisible to Planck’s $C_\ell$ analysis, while $\mathcal{R}(150\times 353)\approx 0.996$ and $1-\mathcal{R}^2\sim 0.7\%$ of the 150 GHz dust power. That is $\delta r\sim (0.4\text{–}1)\times 10^{-3}$ in a clean patch if a 353 GHz template is scaled with a single MBB. Poh & Dodelson (2017) quantified the same line-of-sight “extrapolation noise” from Planck plus Pan-STARRS cloud statistics and found that it becomes important at $r\lesssim 0.0015$. Pelgrims et al. (2021) need only $\sim 15\%$ inter-cloud variation in the 353/217 polarized ratio — not $\Delta\beta=0.3$ everywhere — to match the detected angle effect. The Planck 2018 XI statement that the data suffice for $r\simeq 0.01$ is correct and is not a statement about $r\simeq 10^{-3}$.

Angle rotation as $E\to B$ leakage. A coherent rotation $\delta\psi$ of the dust field mixes $E$ into $B$: $\delta C_\ell^{BB}\approx 4\,\delta\psi^2\,C_\ell^{EE,\mathrm{dust}}$. With dust $E/B\approx 2$, that is $\approx 8\,\delta\psi^2\,C_\ell^{BB,\mathrm{dust}}$. Rotations of $1^\circ$, $2^\circ$, and $3^\circ$ between the template frequency and 150 GHz leave $\sim 0.2\%$, $1\%$, and $2\%$ of the dust $BB$ unmodelled. In the two-layer example with 45° field misalignment and $\Delta\beta=0.3$, $\psi_{\rm eff}$ moves from $11.6^\circ$ at 353 GHz to $\sim 9^\circ$ at 150 GHz, a $2^\circ$ rotation, already a $\delta r\sim 10^{-3}$ effect in a clean patch. Orthogonal fields ($90^\circ$ misalignment) produce no net rotation but a large $\beta^P-\beta^I$ from cancellation; that is a different residual, equally missed by a scalar SED.

Why a 3D model changes the residual, not just the error bar. Free moments and a $\Delta_d$ parameter (BICEP/Keck, following Vansyngel et al. 2017) marginalize this residual, at a cost in $\sigma(r)$. SO forecasts already include a moments-style decorrelation penalty; the expanded-SAT 10-year pessimistic case is $\sigma(r)\approx 0.0014$ with decorrelation versus $\approx 0.0007$ in an optimistic no-decorrelation case. CMB-S4 pipeline comparisons (Abazajian et al. / 2025 foreground-cleaning study) find that a complex dust model (their Model 2) biases recovered $r$ at $\sim(5\text{–}8)\times 10^{-4}$ unless extra foreground freedom is introduced, which then inflates $\sigma(r)$. The 3D+$\mathbf{B}_\perp$+grain model is a proposal to spend external information instead of sensitivity: replace $\Delta_d(\ell)$ and free $\mathcal{W}_1$ with a predicted spatial template whose amplitude is a single regional parameter. If the template is right, the bias drops without the full $\sigma(r)$ penalty of unconstrained moments. If the template is wrong, the residual has a definite morphological form (the leftover layer) and can be tested, which a scalar $\Delta_d$ cannot.

Where the gain is largest. The method does not have to work on the whole sky to matter. Planck’s non-detection of $\mathcal{R}<1$ in broad averages, plus Pelgrims’ detection in a pre-selected sample, imply that the threatened sky fraction is small. A 3D risk map that removes or specially models that fraction — misaligned multi-cloud pixels — is already a $\delta r$ reduction for patch experiments (SO, S4, BICEP Array). For LiteBIRD, the same map becomes a spatially varying prior on the first spin-moment, concentrated where $N\ge 2$ and $\Delta\chi$ is large, rather than a sky-uniform $\Delta_d$.

A sanity bound, kept honest. If the 3D map shows one layer carrying $\gtrsim 90\%$ of the polarized column and tomography shows the remainder misaligned by $\lesssim 20^\circ$, then even $\Delta\beta=0.3$ produces $|\delta\psi_{\rm eff}|\lesssim 0.5^\circ$ between 353 and 150 GHz and a residual $\ll 10^{-3}$ in $r$. Halal et al. (2024) and the Local Bubble reconstructions say this is a common high-latitude configuration. Where the maps show otherwise, the model supplies the correction. The gain is therefore heterogeneous: most of the sky is already safe, and the method’s job is to not leave the unsafe 10–20% of pixels as a coherent bias.

What this does not do. It does not replace multi-frequency coverage, bandpass knowledge, or polarization-angle calibration. An instrumental polarization-angle error is degenerate with dust angle rotation and with cosmic birefringence. It does not remove synchrotron–dust correlation or a possible polarized AME tail. It does not, by itself, solve the reionization-bump foreground at intermediate latitudes where the Edenhofer volume is incomplete. Those are separate terms in the $\delta r$ budget.

4. Pipeline sketch

A concrete sequence from the three ingredients to a CMB dust model, with the residual named at each step.

Stage 0 — products that exist.

1. Build a per-pixel layer catalogue: distance-binned $A_V$ (Edenhofer / Bayestar / Vergely) merged with HI Gaussian components (Panopoulou & Lenz), using 3D extinction as the distance prior on velocity components. Output: $N(\hat n)$, $w_i$, a flag for pixels with $N=1$ versus $N\ge 2$. Residual: layers inside one distance/velocity bin; molecular gas without HI.
2. Fit per-layer intensity SEDs to CIB-cleaned Planck+IRAS (and, where available, Herschel, BLAST, DIRBE) intensities, with a Gaussian-process or voxel smoothness prior on $T_i,\beta_i^I$ (Zelko-style). Output: $T_i^I(\hat n)$, $\beta_i^I(\hat n)$, $\epsilon_i=\tau_{353}/A_V$. Residual: CIB, zodiacal light, $\beta$–$T$ degeneracy, the factor-of-two $\tau/E(B-V)$ variation.
3. Assign per-layer $\chi_i$ from, in priority order: starlight tomography where stars exist; HI-fibre PPV Stokes elsewhere; Local Bubble field model for the nearest layer at $|b|\gtrsim 50^\circ$. Output: $\psi_i=\chi_i+\pi/2$ and an uncertainty $\sigma_{\psi_i}$ that includes the Clark & Hensley morphological scatter. Residual: molecular/ionized sightlines; tangling within a layer.

Stage 1 — grain physics as a transfer, not a fit.

4. For each layer, take $R(V)_i$ from Zhang & Green 3D maps (or a diffuse-ISM default at high latitude), $U_i$ from a 3D ISRF, and $n_{\rm H}$ from HI or a dust-to-gas conversion. Evaluate a bounded ensemble of grain models (astrodust+PAH; THEMIS; silicate+carbon) through $n(a)$, $f_{\rm align}(a;U,n_{\rm H})$, and the millimetre dielectric function. Output: a prior, not a delta function, on $(\beta_i^P-\beta_i^I,\; T_i^P-T_i^I,\; R_{P/p,i},\; p_{0,i}F_i\cos^2\gamma_i)$. Residual: millimetre laboratory opacities, the Planck XXII versus 2018 XI polarized-to-intensity offset, magnetic-dipole / TLS contributions at $\nu\lesssim 100$ GHz.

Stage 2 — polarized forward model.

5. Form the 353 GHz layer Stokes two ways and require consistency:
   • emission path: $\mathcal{P}_i(\nu_0)=p_{0}F_i\cos^2\gamma_i\;\epsilon_i A_{V,i}\,e^{2i\psi_i}$;
   • starlight path, where tomography exists: $\mathcal{P}_i(\nu_0)=R_{P/p,i}\,(q_i+iu_i)$. The starlight path cancels $\gamma$ and $F$. The emission path is the fallback. The sum $\sum_i\mathcal{P}_i$ is constrained to match Planck 353 GHz $Q,U$ within noise and a regional amplitude $\alpha_i$.
6. Propagate in frequency:

Free parameters: CMB $Q,U$; synchrotron (with its own prior); regional $\alpha_i$ (tens of per cent); small $\delta\psi_i$ with tomography-error priors; a global or regional $\Delta\beta^P$ offset. Everything that produces frequency-dependent angles and SED curvature is generated, not parametrized.

Stage 3 — three ways to insert this into an $r$ pipeline.

7. Layered parametric fit (map domain). LiteBIRD’s 15 bands give 30 polarized data per pixel versus $2+2N$ amplitudes and a handful of regional spectral parameters. For $N\le 3$ with fixed layer SEDs the problem is linear in the amplitudes; conditioning comes from the different $Q/U$ ratios of misaligned layers — the same geometric fact that causes the problem. For SO/S4 with 6–9 bands, restrict to $N\le 2$ and regional rather than per-pixel spectral parameters.
8. Predicted spin-moment templates (map or needlet domain). Compute $\mathcal{W}_0$, $\mathcal{W}_1^\beta$, $\mathcal{W}_1^T$ from the layer model. Deproject those SEDs in cMILC/NILC, or fit one amplitude per template in FGBuster. This is the right interface when $N$ is uncertain or the layer catalogue is incomplete: the moments capture the leading distortion with the correct spatial morphology and phase.
9. Cross-spectral prior (BK-style likelihood). Replace $\Delta_d$ and its assumed $\ln^2(\nu_1/\nu_2)$ and $\ell$ scalings with

from the layer maps, including the layer–layer cross-spectra that vanish only if the layers are spatially uncorrelated. This is the honest covariance: LOS mixing is not a scalar $\mathcal{R}(\nu_1,\nu_2)$.

Stage 4 — validation before cosmology, then staged deployment.

10. Test the forward model on quantities that do not require a primordial $B$-mode:
    • the $\psi(217)-\psi(353)$ map versus Pelgrims-style HI-selected samples;
    • residual-map covariances of Guillet et al. (2026);
    • frequency dependence of dust $E/B$, $TE$, $TB$, and $EB$ (Ritacco et al. 2023; Vacher et al. 2023);
    • cross-spectra of per-layer templates with BK/SO deep-field $B$ modes (as in BICEP/Keck XVI);
    • $\sum_i\mathcal{P}_i$ versus Planck 353 GHz, and starlight-predicted versus emission-predicted $\mathcal{P}_i$ in the 4 deg² tomography patch, then in PASIPHAE fields.
11. Deploy in order of decreasing robustness:
    • Now: risk maps and $\Delta_d$ priors from cloud counts, HI angles, and 3D extinction, used for field selection, masking, and a spatially varying $\Delta_d(\hat n)$ rather than a scalar.
    • As PASIPHAE, 3D emission, and 3D $R(V)$ mature: the full per-layer model as a prior, with $\alpha_i$ and $\delta\psi_i$ still free.
    • With LiteBIRD and CMB-S4 data: joint fit, in which the microwave data refine $\alpha_i$, $\delta\psi_i$, and $\Delta\beta^P$, while the 3D structure remains the regularizer that prevents the moment–$r$ degeneracy.

Named residuals after the full pipeline. An error $\delta\alpha_i$ leaves a fraction of layer $i$ at the wrong SED. An error $\delta\psi_i$ leaks $\sim 8\,\delta\psi_i^2$ of that layer’s $BB$. A missing layer reproduces the original problem in miniature. An error in $\Delta\beta^P$ is a coherent, partly band-fit-able bias. Unresolved tangling inside a layer remains as a statistical high-$\ell$ $\Delta_d$. Bandpass mismatch and polarization-efficiency errors mimic SED curvature and must be carried in the instrument model. SRoll2 versus PR4 disagreement (Guillet et al. 2026) shows that current Planck systematics are already comparable to the astrophysical signal being modelled; the pipeline cannot be validated on a Planck split that is itself the uncertainty.

5. What to keep from the draft

The draft is already a serious physical essay, not a list of surveys. The following pieces should survive essentially intact.

1. The four-way failure of a single MBB per pixel — non-MBB effective SED, frequency-dependent angle, $\beta^P\neq\beta^I$ from geometry, and $\mathcal{R}<1$ — with the explicit statement that angle rotation is spin-2 and cannot be absorbed by a scalar SED. This is the right diagnosis of the question’s “decorrelation of the spectral index.”

2. The two-layer numerical table ($\Delta\beta=+0.3$ at 45° and 90°; $\Delta T=-4$ K at 45°) and the accompanying $E\to B$ conversion $\delta C_\ell^{BB}\approx 8\,\delta\psi^2 C_\ell^{BB,\mathrm{dust}}$. The 11.6° figure is the correct Stokes arithmetic for a 30%/70% split and 45° field misalignment (polarization angles differ by 45°, Stokes vectors by 90°). That table is the paper’s core pedagogic result.

3. The $\mathcal{R}(217\times 353)\approx 0.999 \Rightarrow \sim 0.7\%$ residual at 150 GHz calculation. It is the correct way to connect a Planck-invisible correlation ratio to a $\delta r\sim 10^{-3}$ residual, and it should be kept next to Planck 2018 XI’s “fine for $r\approx 0.01$” conclusion so the two are not misread as contradictory.

4. Why extra free parameters fail, citing McBride, Bull & Hensley (2023) and Vacher et al. (2022). This is the logical warrant for introducing external data at all.

5. The Tassis–Pavlidou starlight transfer $R_{P/p}\times(q_i,u_i)\to(Q_i,U_i)$, with the cancellation of $\cos^2\gamma$ and $F$. This is the unique reason grain sizes and the projected field belong in the same sentence as 3D emission maps.

6. Pelgrims et al. (2021) as the existence proof that independent 3D + $\mathbf{B}_\perp$ information can select the pixels where the effect is detectable, and the $\sim 15\%$ inter-cloud SED ratio as the empirical target, not an assumed $\Delta\beta=0.3$ everywhere.

7. The linear-in-amplitudes structure once layer SEDs and angles are supplied, and the observation that misaligned layers are well-conditioned because they have different $Q/U$ — the cause of the problem is also the source of the leverage.

8. The three usage modes (layered parametric fit, spin-moment templates, cross-spectral prior) and the staged deployment (risk maps now; full model as PASIPHAE matures; joint fit with LiteBIRD/S4). This is an honest operations plan, not a claim that the catalogue is already complete.

9. The limitations section’s specific floors: 1.25 kpc angular/distance resolution; $\tau_{353}/E(B-V)$ variation; $R(V)$ as a one-parameter proxy for a size range that is not identical to the millimetre-emitting grains; 4 deg² of published degree-scale tomography; HI failure in molecular/ionized gas; SRoll2 versus PR4; degeneracy of dust angle rotation with instrumental angle error and birefringence. These are the right caveats. They should not be softened.

10. Planck 2018 XII phenomenology used as calibration, not as decoration: $p$ up to $\sim 22\%$, $R_{P/p}\approx 5.4\,\mathrm{MJy\,sr}^{-1}$, $p$–$S$ anticorrelation, $BB/EE\sim 0.5$, $C_\ell\propto\ell^{-2.4}$. Those numbers set $p_0 F\cos^2\gamma$, the starlight transfer, and the $E\to B$ conversion.

6. Upgrade list

Numbered improvements the primary should adopt.

1. Rename the problem. Lead with “line-of-sight vector superposition” (Tassis & Pavlidou), not “decorrelation of the spectral index.” Then split $\mathcal{R}<1$ into LOS-plus-misalignment (frequency-dependent $\psi$, $E/B$, $EB$) versus plane-of-sky $\beta$ variation (no angle rotation). The three ingredients are necessary for the first and only partially relevant to the second.

2. Put the Local Bubble two-layer cartoon in front. At the latitudes of recombination-bump experiments the geometry is often Bubble wall plus a subdominant distant screen. Quote O’Neill et al. (2024, 2025) and Halal, Clark & Tahani (2024) together: dominance of the Bubble is common; inclination-only models of $p$ fail; 3D complexity still depolarizes. That is the regime in which $N\le 2$ is not a simplification but the sky.

3. State the subdominant-layer theorem early. Equal columns are not required. A 10–30% polarized second layer with $\Delta\chi\sim 45^\circ$ and a modest SED difference is the $\delta r\sim 10^{-3}$ threat. 3D maps earn their keep by measuring that second-layer weight, not by reconstructing a continuous emissivity cube at CMB resolution.

4. Make the I$\to$P transfer and $R_{P/p}$ the central grain-physics product. Compress the grain-model survey. The microwave analysis needs a bounded prior on $\beta^P-\beta^I$, $T^P-T^I$, $R_{P/p}(R_V,U)$, and which layers polarize. It does not need to settle astrodust versus THEMIS in the same paper. Say explicitly: millimetre $\beta$ is a dielectric-function quantity (size-independent in the Rayleigh limit); size enters through $T(a)$ and $f_{\rm align}(a)$.

5. Add 3D $U$ as a predicted $T_i$, and 3D $R(V)$ as a size map. These are the grain-physics observables that actually exist (Gaia XP; Zhang & Green 2025) and that couple to alignment. Do not leave $T_i$ as a free Zelko parameter if the heating can be predicted.

6. Treat HI PPV Stokes as the all-sky $\chi_i$ of 2026, and PASIPHAE as the calibrator. The draft lists them as equals. Operationally, Clark & Hensley (2019) plus BICEP/Keck XVI are what can be run on the full high-latitude sky this year; starlight tomography is the higher-fidelity subset and the $R_{P/p}$ calibration.

7. Do not overclaim Faraday $\gamma$. Move RM to a weak prior. State the $n_e$ versus $n_d$ weighting in one sentence.

8. Promote predicted spin-moments and $C_\ell^{ab}=\sum_{ij}\hat g_i(a)\hat g_j(b)C_\ell^{ij}$ to primary CMB interfaces. A layered Stokes cube that is not reduced to a template, a cMILC constraint, or a BK covariance will not enter an $r$ likelihood. Write those three interfaces at the same depth as the ISM model.

9. Add dust $EB$ / birefringence as a second payoff. Frequency-dependent dust $\psi$ is a chromatic $EB$ source. The same 3D model that protects $r$ protects isotropic birefringence. Clark et al. (2021) and Ritacco et al. (2023) already flag the degeneracy.

10. Keep the residual-to-$r$ arithmetic, but quote the patch dust amplitude. Write $\delta r \sim f_{\rm res}\times r_{\rm dust}(150\,\mathrm{GHz},\ell=80)$ with $r_{\rm dust}\sim 0.05$–$0.1$ in current deep fields, $f_{\rm res}\sim 1-\mathcal{R}^2$ or $8\,\delta\psi^2$. Cite Poh & Dodelson (2017) for the $r\lesssim 0.0015$ threshold and Pelgrims’ 15% for the empirical SED contrast. Distinguish “1% of a 353 GHz template extrapolation” from “1% after a 3D layered fit.”

11. Use Planck 2018 XI and Pelgrims 2021 as complementary, not competing. No $C_\ell$ detection in broad averages and a detection in a 3D-selected sample is the evidence that spatially adaptive modelling wins. That is also why risk maps are the Stage-0 product with immediate cosmological use.

12. Add the high-frequency and submillimetre polarized handles (LiteBIRD HF, SO 280 GHz, BLASTPol/PILOT/NIKA2-Pol, CCAT) as empirical $T^P$, and Fuskeland et al. (2023) as the quantitative case that $T$ variation (if dominant) is a high-frequency problem.

13. Name CIB, zodiacal light, and $\tau/E(B-V)$ variation as part of the 3D emission step, not as a late caveat. A 3D extinction map is not a 3D emission map until those are modelled; Zelko et al. (2022) already showed the conversion varies by $\times 2$.

14. Keep a statistical within-layer residual. Do not imply that a single $\psi_i$ per cloud removes high-$\ell$ decorrelation. Retain $F_i$ or $S_i$ from MHD/filament models; that residual is what a small $\Delta_d(\ell)$ is still for after the layered model is applied.

15. Add an explicit validation protocol on non-$r$ observables (the $\psi(217)-\psi(353)$ map, residual covariances, $E/B(\nu)$, layer-template $\times$ deep-field $B$, starlight versus emission $\mathcal{P}_i$ in the 4 deg² patch) before claiming a $\delta r$ improvement. The SRoll2/PR4 split is a prerequisite, not an appendix.

16. Quantify PASIPHAE’s current incompleteness without undercutting the method. 4 deg² of degree-scale tomography is a pathfinder. The steelman is: the form of the correction is known, HI gives the all-sky approximation, and PASIPHAE replaces HI angles with starlight angles plus $R_{P/p}$ amplitudes over the high-latitude sky. That is a staged reduction of $\sigma_{\psi_i}$ and $\sigma_{\alpha_i}$, not a binary “on/off” for dedusting.

17. Include synchrotron–dust spatial correlation and a possible polarized AME tail as out-of-model terms that live in the same bands, so the dust model is not asked to fit them.

18. Drop or qualify unverifiable precision on “status in 2026” for datasets that are in progress. Keep PASIPHAE, Edenhofer 1.25 kpc / 14′, Zelko 27′–110′, and 4 deg² tomography as the documented state. Do not let a status table imply that per-layer $\psi_i$ is already all-sky at degree scale.

19. Add the conditioning argument in one sentence in the pipeline. Misaligned layers are linearly independent in $(Q,U)$ even when their SEDs are close; aligned layers with different SEDs are independent in frequency. The data can separate layers when at least one of those is true. When neither is true, the 3D model correctly predicts that a single MBB is enough.

20. Close with a falsifiable claim. If, in PASIPHAE-overlap fields, the layered model does not reduce the 217–353 GHz angle difference in Pelgrims-selected pixels relative to a single-MBB fit, the I$\to$P transfer or the layer catalogue is wrong and the method must not be used to shrink $\sigma(r)$. Conversely, a successful reduction there is the green light to promote predicted $\mathcal{W}_1$ templates into the cosmological likelihood.