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

Dedusting CMB B-modes with 3D dust columns, grain physics, and the sky-projected magnetic field

Polarized thermal emission from Galactic dust is the leading astrophysical contaminant of primordial CMB B-modes at frequencies $\nu\gtrsim 70\,\mathrm{GHz}$. Component-separation pipelines that treat each pixel, or each power-spectrum bin, as a single modified blackbody (MBB) with a frequency-independent polarization angle fail when the effective dust spectral energy distribution (SED) and/or the polarization angle $\psi$ become frequency-dependent. The CMB literature often calls this spectral-index decorrelation. That name is too narrow. The operational problem is polarized-dust frequency decorrelation: the dust sky at a CMB frequency is not a rescaling of the dust sky at $353\,\mathrm{GHz}$, either in amplitude or in angle. Spectral-index variation is only one cause.

Write the observed complex linear polarization as a sum over magnetically coherent layers (clouds) $c$,

Here $A_c$ contains dust column, alignment efficiency and field-inclination factors, $s_c(\nu)$ is the layer’s normalized polarized SED, and $\psi_c$ is its plane-of-sky polarization angle. This decomposition isolates the necessary physics:

• Multiple clouds do not by themselves cause frequency decorrelation. If every layer has the same $s_c(\nu)$, that SED factors out and the polarization pattern scales coherently with frequency despite field tangling.
• Different SEDs do not by themselves rotate the polarization angle if all layers share the same $\psi_c$.
• Line-of-sight (LOS) angle decorrelation becomes important when layers have both different polarized SEDs and non-parallel polarization vectors. Relative weights then change with frequency, rotating the net angle and changing the degree of cancellation in $Q$ and $U$.

That is the two-cloud mechanism of Tassis & Pavlidou (2015). Three observational ingredients—3D maps of the dust column (presently Gaia extinction cubes, not 3D thermal-emission maps), knowledge of the grain size distribution together with composition, shape and alignment, and maps of the plane-of-sky magnetic field—do not uniquely reconstruct the microwave sky. They constrain the latent variables in the sum above: the number, location and relative weights of layers ($A_c$); a family of priors on $s_c(\nu)$; and the layer angles $\psi_c$. They identify where extra spectral freedom is warranted. They do not supply per-layer temperature and opacity index as measured quantities. That gap is the central limitation of any “physical” forward model built from these data, and it must be carried through the pipeline rather than assumed away.

Throughout, Stokes parameters $(I,Q,U)$ follow the IAU convention, $p\equiv\sqrt{Q^2+U^2}/I$, and $\psi=\tfrac12\mathrm{atan2}(U,Q)$. For optically thin thermal dust the polarization direction is orthogonal to the plane-of-sky projection of $\mathbf{B}$. A single homogeneous volume element is an MBB,

with a conventional reference $\nu_0=353\,\mathrm{GHz}$. Polarized emission from that element is

where $p_0$ is the intrinsic polarizing efficiency of the aligned grain population, $\gamma$ is the angle between $\mathbf{B}$ and the plane of the sky, and the factor $\cos^2\gamma$ is the geometric projection of aspherical grain emission (Lee & Draine 1985). As soon as several such elements with different $\{\beta_d,T_d,p_0,\gamma,\psi\}$ occupy the same beam, neither $p$ nor $\psi$ nor the effective $\beta_d$ is frequency-independent. All sky models below are understood to be integrated through experimental bandpasses, with colour corrections, rather than evaluated at nominal centre frequencies.

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

Inflationary gravitational waves source a curl (B-mode) pattern in CMB polarization on degree scales, with amplitude parameterized by the tensor-to-scalar ratio $r$ at a pivot $k=0.05\,\mathrm{Mpc}^{-1}$. The present observational bound is $r_{0.05}<0.036$ at 95% confidence from BICEP/Keck through the 2018 season (BK18; BICEP/Keck Collaboration 2021). Combining BK18 with Planck PR4, BAO and CMB lensing tightens this to $r<0.032$ (Tristram et al. 2022). LiteBIRD targets a total uncertainty $\delta r\sim 10^{-3}$, including foreground residuals. The Simons Observatory (SO) and BICEP Array aim at similar or somewhat looser deep-field sensitivity. The US CMB-S4 project, often treated in older forecasts as the Stage-4 deep-field benchmark, was cancelled by DOE and NSF on 9 July 2025; the scientific requirement it encoded—control of Galactic residuals near $D_\ell^{BB}\sim 10^{-5}\,\mu\mathrm{K}^2$ on $\ell\lesssim 100$ in a clean patch (McBride, Bull & Hensley 2023)—remains the relevant target for whatever deep-field programme replaces it.

At these levels Galactic polarized dust is not a small correction. In the BICEP/Keck field, dust B-modes at $150\,\mathrm{GHz}$ already exceed the $r=0.01$ tensor spectrum by an order of magnitude. Experiments targeting $\sigma(r)\sim 10^{-3}$ cannot tolerate residuals merely because they are a small percentage of the original dust map. What matters is residual $BB$ power around the reionization peak ($\ell\lesssim 10$), which dominates LiteBIRD’s $r$ constraint over most of the sky, and around the recombination peak ($\ell\sim 80$), which dominates deep southern patches.

Thermal dust emission is the long-wavelength tail of radiation from large interstellar grains in equilibrium with the interstellar radiation field. The mean SED in both intensity and polarization is well approximated by an MBB over the Planck HFI range. Planck intermediate results XXII found

for a mean colour temperature $T_d=19.6\,\mathrm{K}$ on intermediate-latitude patches (Planck Collaboration Int. XXII 2015). Planck 2018 XI revises the polarized index to $\beta_d^{\mathrm{P}}=1.53\pm 0.02$ at the same $T_d$, and gives

a $1.7\sigma$ offset, not a statistically decisive detection that the grains dominating $I$ differ from those dominating $Q$ and $U$. The same revision implies a $353\to 70\,\mathrm{GHz}$ decline of the polarization fraction of only about $8\%\pm 5\%$ if the offset is attributed to $\beta_d^{\mathrm{P}}\neq\beta_d^{\mathrm{I}}$ at fixed $T_d$, weaker than the $(21\pm 6)\%$ decline reported in 2015. A $\beta_d^{\mathrm{P}}\neq\beta_d^{\mathrm{I}}$ offset can arise with a single grain composition if polarization preferentially weights LOS segments with ordered, sky-plane fields that have systematically different $T$ or $\beta$. Two-population models remain viable; they are not required by the 2018 numbers.

The polarization itself is produced because aspherical grains align with their short axes preferentially parallel to $\mathbf{B}$. Radiative alignment torques (RATs) are the leading mechanism in the diffuse ISM (Dolginov & Mitrofanov 1976; Draine & Weingartner 1996, 1997; Lazarian & Hoang 2007; Andersson, Lazarian & Vaillancourt 2015). Alignment is size-dependent: only grains larger than a critical alignment radius $a_{\mathrm{align}}$ rotate suprathermally and remain aligned against gas randomization. In the typical diffuse ISM this threshold is of order $0.05\,\mu\mathrm{m}$; the often-quoted $0.055\,\mu\mathrm{m}$ prefactor in Hoang & Lazarian (2014) is normalized to denser, cooler gas ($n_{\mathrm{H}}=10^3\,\mathrm{cm}^{-3}$, $T=10\,\mathrm{K}$, anisotropy factor $\gamma=0.1$) and is not directly the diffuse-ISM value. Polarized submillimetre emission is therefore a large, equilibrium-grain signal. PAHs and other nanoparticles are not a generic explanation for thermal-dust decorrelation at CMB frequencies. Anomalous microwave emission and possible magnetic dipole emission from ferromagnetic inclusions are separate, uncertain components.

Planck 353 GHz $Q,U$ maps remain the highest-signal-to-noise full-sky tracer of polarized dust. They reveal several statistical facts that any B-mode foreground model must reproduce:

• A maximum observed polarization fraction $p_{353}^{\max}\sim 0.20$–$0.22$, a lower bound on the intrinsic grain polarizing efficiency after LOS depolarization (Planck Collaboration XII 2018).
• An $E/B$ power asymmetry $C_\ell^{EE}/C_\ell^{BB}\simeq 2$ on intermediate scales, together with a positive $TE$ correlation, both of which arise because filamentary density structures are preferentially aligned with $\mathbf{B}_{\mathrm{POS}}$ (Planck Collaboration Int. XXXVIII 2016; Clark et al. 2015).
• A tight anti-correlation between $p$ and the polarization-angle dispersion function $\mathcal{S}$, consistent with depolarization by a turbulent field rather than by spatially variable grain physics (Planck Collaboration Int. XIX 2015; Planck Collaboration XII 2018).
• Residual spatial variation of the polarized SED at a level that is relevant to forecasts for $\ell\sim 80$ B-mode analyses, even where a mean MBB remains an adequate description of the sky-averaged SED (Planck Collaboration Int. L 2017; Ritacco et al. 2023).

Empirical status of frequency decorrelation. Current B-mode studies do not “suffer from” a measured bias on $r$ caused by decorrelation. They face a forecast problem: modest SED complexity that is allowed by Planck, and in some cases detected in selected pixels, would bias next-generation $r$ analyses if left unmodelled.

• Pelgrims et al. (2021) found excess 353-to-217 GHz polarization-angle differences and extra depolarization on H I-selected high-latitude sightlines with multiple, misaligned clouds. They explicitly caution that the paper “does not address whether the amplitude of the effect is large enough at the sky-map level to affect any particular experiment’s search for primordial B-modes.”
• Planck 2018 XI finds “no evidence for a loss of correlation” in multi-frequency cross-spectra over large high-latitude regions.
• BK18’s decorrelation parameter $\Delta_d$ (the 217$\times$353 correlation ratio, scaled in $\ell$ following Planck Int. L) has a posterior that peaks at unity. BK does not scale a 353 GHz template by an MBB; it fits a parametric dust-plus-synchrotron model to all auto- and cross-spectra of the BK, WMAP and Planck bands, measuring $\beta_d=1.49^{+0.13}_{-0.12}$ in-field (BICEP/Keck Collaboration 2021).
• BICEP/Keck XVI ran an IVC-versus-LVC two-cloud test in the BK field with H I morphology templates and found “no evidence for decorrelation … from the inclusion of dust associated with the intermediate velocity H I” (BICEP/Keck Collaboration 2023).
• SPIDER, in its southern patch, finds no significant LOS polarization-angle decorrelation, reports $\beta_d=1.45\pm 0.05$ ($1.47\pm 0.06$) for $E$ ($B$), and also reports that two sub-regions differ in $\beta_d$ at $3.9\sigma$ and that one sub-region is inconsistent with a simple MBB at more than $3\sigma$ (SPIDER Collaboration 2025).
• Guillet et al. (2026), using spin-moment residuals of Planck polarization, find that SRoll2 is compatible with a dust model that includes only temperature variations rather than spectral-index variations, while PR4 is incompatible with both—an important reminder that the inferred mix of $\Delta T$ versus $\Delta\beta$ is release-dependent.

Planck Int. L reported $\mathcal{R}_\ell^{BB}(217,353)$ departing from unity at $\ell\gtrsim 50$ in some large high-latitude regions, with an indicative $\Delta\beta_d^{\mathrm{P}}\sim 0.07$ “suggested by this simple analysis” for one region, and used $\mathcal{R}_\ell^{BB}(150,353)=0.85$ as an illustrative extrapolation, not as a measurement. Subsequent analyses (Sheehy & Slosar 2018; Planck Collaboration XI 2020) found $\mathcal{R}_\ell^{BB}(217,353)$ consistent with 1 within Planck noise over large sky fractions, while cautioning that the ratio can be spatially variable and that the effect is expected to grow toward smaller scales and larger frequency lever arms.

Component separation for $r$ is usually done in one of three families: (i) parametric pixel-domain fits (Commander, FGBuster); (ii) internal linear combination / needlet ILC methods; (iii) template multi-frequency likelihoods at the power-spectrum level, as in BICEP/Keck, which fit a parametric dust SED to all auto- and cross-spectra, sometimes with a decorrelation parameter $\Delta_d$. All three families become biased if the true dust SED is not a single MBB, or if $Q$ and $U$ do not share the same frequency scaling. That is the operational problem. It is presently a requirement on the model, not a demonstrated bias in published $r$ limits.

2. What frequency decorrelation means operationally

2.1 Two distinct, observationally mixed phenomena

Spatial (2D) decorrelation. The dust SED varies across the sky. If $\beta_d(\hat n)$ or $T_d(\hat n)$ is not constant, a map at frequency $\nu_1$ is not a scalar multiple of the map at $\nu_2$. At the power-spectrum level this is quantified by the correlation ratio

A single MBB with uniform $\beta_d$ implies $\mathcal{R}_\ell=1$. Gaussian spatial variations of $\beta_d$ with rms $\Delta\beta_d$ produce

on the scales where the $\beta_d$ map is coherent (Planck Collaboration Int. L 2017). BK18’s $\Delta_d$ is this correlation ratio evaluated at $217\times 353$ and scaled in $\ell$; it is not the same object as a free $\mathcal{R}_\ell$ in the Vansyngel et al. (2017) multi-layer LOS simulations.

Line-of-sight (LOS) decorrelation. Even a single pixel can fail to be a single MBB. If the LOS intercepts two or more clouds with different SEDs and misaligned $\mathbf{B}_{\mathrm{POS}}$, the Stokes parameters add. Because each cloud’s polarized-intensity ratio

is frequency-dependent whenever $T_1\neq T_2$ or $\beta_{d,1}\neq\beta_{d,2}$, the net polarization angle $\psi_{\mathrm{tot}}(\nu)=\tfrac12\mathrm{atan2}(U_\nu,Q_\nu)$ rotates with frequency, and the net polarization fraction is depolarized relative to either cloud.

The magnitude of that rotation is easy to overstate. Tassis & Pavlidou (2015) show that if the plane-of-sky fields “form an angle of less than $60^\circ$, then … the difference in polarization angle [is] less than $5^\circ$, independently of the other parameters,” while exceeding $20^\circ$ needs $\alpha\approx 85^\circ$. Direct computation for equal polarized intensities at $353\,\mathrm{GHz}$, $T=16$ and $22\,\mathrm{K}$, $\beta=1.5$, between 353 and $150\,\mathrm{GHz}$:

“Tens of degrees” requires near-orthogonal fields and a weight ratio that moves a substantial fraction of the way through unity. That is exactly where $p\to 0$ and the polarized power that biases $r$ is smallest. For a more realistic ISRF-driven $\Delta T\sim 1$–$2\,\mathrm{K}$ at fixed $\beta$, the 353-to-150 GHz weight ratio shifts by only $1$–$3\%$, and $\Delta\psi$ is a fraction of a degree at $\alpha=60^\circ$ and a few degrees even at $\alpha=80^\circ$. The 90° flip remains a measure-zero pathology of $\alpha\to 90^\circ$ at $\rho_\nu=1$.

LOS decorrelation is not equivalent to a spatially varying but frequency-independent $\beta_d$ per pixel. A pixel-wise MBB fit still assumes that $Q$ and $U$ share a common SED. Mixing of misaligned clouds violates that assumption: the effective spectral indices of $Q$ and of $U$ differ. McBride, Bull & Hensley (2023) make the corresponding point about parametric recovery: a two-parameter MBB fit to superposed MBBs “generally results in unbiased estimates of the CMB Stokes $Q$ and $U$ amplitudes in each pixel,” with large biases only “when there are significant changes in both the dust SED and polarization angle along the line of sight.” Intra-cloud temperature spread from a grain-size distribution, without angle variation, is not the dangerous case.

Pelgrims et al. (2021) provided the first direct evidence that the LOS mechanism is present in Planck data. Using H I to count clouds per sightline (Panopoulou & Lenz 2020) and H I filament orientation to estimate $\mathbf{B}_{\mathrm{POS}}$ per cloud (Clark & Hensley 2019), they selected a target sample of high-latitude pixels with column-weighted cloud count $\mathcal{N}_c\gtrsim 1.5$ and cloud–cloud misalignment $\gtrsim 60^\circ$, and a control sample of single-cloud sightlines. The target sample shows a statistically significant excess of 353-to-217 GHz polarization-angle differences and extra depolarization, with a significance that increases after residual Planck systematics are reduced (SRoll2). An intrinsic cloud-to-cloud variation of $\sigma_\rho\sim 0.15$ in the 353/217 polarized intensity ratio is sufficient to reproduce the measured effect in that selected sample. This is evidence for the mechanism, not a measurement of a universal all-sky decorrelation amplitude.

Translating 15% in $P_{353}/P_{217}$ into MBB parameters, for a conventional $T_d=19.6\,\mathrm{K}$, $\beta_d=1.53$:

These exceed typical Planck polarized-SED dispersions. Either part of the Pelgrims detection is selection or residual systematics, or composition (or porosity, or alignment efficiency) varies cloud-to-cloud at a level that a global grain-size distribution does not predict. A 3D extinction map plus a local ISRF estimate typically yields $\Delta T\sim 1$–$2\,\mathrm{K}$ and $\sim 1\%$ ratio changes, not 15%. Grain-size knowledge is then the wrong primary lever for the Pelgrims number; composition and environment are.

2.2 Why a single MBB per pixel fails

The MBB is a nonlinear function of $\beta_d$ and $T_d$. Averaging MBBs—along the LOS, inside the beam, or when computing spherical-harmonic coefficients over a large mask—does not return an MBB. Chluba, Hill & Abitbol (2017) formalized the resulting SED distortions as a moment expansion about a fiducial $(\bar\beta_d,\bar T_d)$. In polarization the moments are spin-2 complex quantities (Vacher et al. 2023, 2025). After the expansion pivot is chosen so that the real part of the first-order moment vanishes, that first-order moment is purely imaginary and is the leading term in the spectral rotation of $\psi$. Vacher et al. (2025) find that the imaginary parts of the first-order $\beta$ and $T$ moments are the main drivers of $BB$ decorrelation in a layered dust model. A residual stage that applies scalar moments identically to $Q$ and $U$ cannot represent the effect this essay is about.

A further degeneracy, not always stated, is intrinsic $\beta(\nu)$ curvature of a single grain composition, which is degenerate with LOS mixing of several MBBs (Kogut & Fixsen 2016). High-frequency polarimetry is needed to separate the two.

Operational consequences for component separation:

1. Template and harmonic-space methods. Scaling a 353 GHz $(Q,U)$ template by a single MBB, or fitting a single $\beta_d$ (and a single $\Delta_d$) to all auto- and cross-spectra, mis-predicts the dust at 90–150 GHz wherever $\psi(\nu)$ or the effective polarized $\beta_d$ differs from the model. The residual is a B-mode field that is spatially correlated with the Galaxy and spectrally intermediate between dust and CMB.

2. Pixel-domain parametric fits. A per-pixel MBB is misspecified whenever both the SED and the angle vary along the LOS (McBride et al. 2023). Hensley & Bull (2018) found that LOS mixing and iron grains are the two physical complexities that most stubbornly bias recovered CMB Stokes parameters, even when a generalized MBB with extra phenomenological parameters is fitted.

3. Harmonic-space likelihoods. Fitting $C_\ell^{BB}(\nu_i\times\nu_j)$ with a single $\beta_d$ and a single decorrelation parameter can absorb a mean loss of coherence but cannot represent frequency-dependent $\psi$ or the distinct $\beta$ and $T$ moment spectra. Mangilli et al. (2021) showed that ignoring SED distortions, or replacing them by one $\mathcal{R}_\ell$, biases $r$ above the sensitivity of next-generation experiments. Vacher et al. (2022) found that LiteBIRD-like simulations with spatially varying $(\beta_d,T_d)$ produce a $\sim 21\sigma_r$ bias if a single MBB is assumed ($\hat r\simeq 1.25\times 10^{-2}$ on $f_{\mathrm{sky}}=0.7$); expanding in $\beta_d$ alone is still insufficient; expanding in both $\beta_d$ and $T_d$ recovers an unbiased $r$ with $\sigma_r=8.8\times 10^{-4}$ after an $\ell$-cut that limits the $T$ moments to large scales. Those numbers depend on the simulated sky (PySM-like spatial variation, not a 3D cloud catalogue) and on the sky fraction; they establish that modest-looking SED complexity can matter at LiteBIRD sensitivity, not that any particular 3D reconstruction removes a stated fraction of that bias.

4. Blind ILC. Standard ILC nulls a single dust SED. Frequency-dependent mixing leaves a dust residual that ILC interprets as CMB. Moment-constrained ILC (Remazeilles, Rotti & Chluba 2021) deprojects the first few moment templates, at the cost of extra noise and of needing enough frequency channels. Different component-separation methods are sensitive to different properties of the spin-moments (Vacher et al. 2025).

Sponseller & Kogut (2022) give a complementary threshold: an LOS standard deviation $\sigma_T\gtrsim 1.6\,\mathrm{K}$ or $\sigma_\beta\gtrsim 0.045$ produces a nontrivial bias on $r$ in single-MBB fits, of order $\delta r\sim 1.4\times 10^{-4}$ for the temperature case in their PICO-like setup—negligible for current upper limits, relevant at $\delta r\sim 10^{-3}$. Ritacco et al. (2023) show that intensity-based models underestimate polarized complexity and that $BB$ spectra do not match first-order moments fitted to intensity.

In short: frequency decorrelation means that the dust sky at 150 GHz is not a rescaling of the dust sky at 353 GHz, either in amplitude or in polarization angle, and that the residual after a one-component subtraction has a B-mode spectrum that can be degenerate with $r$ at next-generation sensitivity. The rest of this essay is about constraining the missing layers rather than fitting their average—and about being precise on which layer properties the ancillary data actually measure.

3. How each ingredient helps on its own

3.1 3D maps of the dust column (extinction, not thermal emission)

A 2D dust map, however high its S/N, is a column integral. Decorrelation is generated by differences among the emitters that were integrated. A 3D column map converts that integral into a sum of discrete layers, each of which can be assigned its own SED prior and its own $\mathbf{B}$.

What actually exists today are primarily 3D extinction / reddening maps, not 3D maps of submillimetre emission. They are the right starting point only insofar as, in the diffuse ISM, dust mass, optical extinction, far-infrared emission, and H I column are tightly correlated (Boulanger et al. 1996; Lenz, Hensley & Doré 2017). Gaia parallaxes plus stellar photometry and low-resolution BP/RP spectra have produced a sequence of maps with steadily improving distance resolution:

• Green et al. (2019), Bayestar19: Gaia + Pan-STARRS 1 + 2MASS, northern sky.
• Leike, Glatzle & Enßlin (2020) and subsequent information-field-theory reconstructions: parsec-scale structure in the local ISM.
• Edenhofer et al. (2024): 54 million stars with Gaia XP-based distances and extinctions (Zhang, Green & Rix 2023; ZGR23), reconstructed with variational inference on a spherical grid from 69 to 1250 pc in 516 logarithmically spaced distance bins, angular resolution up to $N_{\mathrm{side}}=256$ ($\sim 14'$).
• All-sky Gaia+LAMOST reconstructions extend the distance reach to several kpc at coarser angular resolution.

Zeroth-order 3D fact: the Local Bubble. Skalidis & Pelgrims (2019) showed that at $|b|\ge 60^\circ$ the Planck 353 GHz polarized sky is dominated by a magnetized structure at 200–300 pc that coincides with the Local Bubble shell. Alves et al. (2018) and Pelgrims et al. (2020) fitted a one-shell magnetic model to the polar caps. Any layered model must be shown to beat this one-shell baseline at high latitude, not merely to add layers.

Halal, Clark & Tahani (2024) performed the test that the 3D-map programme actually needs: they used the Edenhofer map, defined a column-weighted component count with thresholds $\le 1.1$ and $\ge 1.5$, and compared to Planck $p_{353}$ at $80'$. Multi-component sightlines are depolarized by about 1.5% in absolute $p$ relative to single-component sightlines of similar total column. They found no imprint of Local Bubble geometry on $p_{353}$, and the multi-component signal disappears when the dust is restricted to 270 pc. The extended 3D dust distribution, not solely the Bubble wall, drives the observed variation of $p$ at $80'$. That is a real, quantitative use of Gaia 3D maps for polarization, and it is a geometric depolarization result, not a measurement of frequency decorrelation.

At high Galactic latitude much of the dust that dominates 353 GHz emission lies within a few hundred parsecs. That statement must be qualified. LiteBIRD’s $r$ sensitivity is dominated by $\ell<10$ over most of the sky, where dust beyond 1.25 kpc and near the plane contributes even after masking. The 3D per-cloud model is primarily a small-field, $\ell\approx 80$ tool for deep southern patches (BK, SPIDER, SO). Completeness is also reversed at the lowest columns: ZGR23 per-star extinction precision is about 0.03 mag, while the BK field has $E(B-V)\approx 0.01$–$0.03$ mag. Edenhofer et al. note that the Gaussian-extinction assumption is poor at extremely low extinction because of the catalogue’s positivity prior. Decomposing a $0.05\,\mathrm{mag}$ column into distance layers with $0.03\,\mathrm{mag}$ per-star errors is prior-dominated. The 3D maps are not most complete precisely where CMB deep fields need them.

Panopoulou & Lenz (2020) find on average three H I clouds per $0.2\,\mathrm{deg}^2$ pixel in the northern cap and 2.5 in the south. The column-weighted complexity $\mathcal{N}_c$ is lower: $\mathcal{N}_c<1.5$ on 60% of northern and 90% of southern high-latitude pixels. “Most high-latitude pixels are one-cloud” is true of $\mathcal{N}_c$ in the south, not of the raw cloud count, and not of the northern cap where IVCs matter.

What a 3D column map gives a component-separation pipeline, even before grain physics or $\mathbf{B}$ are added:

1. A physically motivated layering of $A_c$. Instead of one MBB per pixel one writes a sum over distance slices, with $\tau_i$ proportional to the differential extinction $\Delta A_V$ (or $\Delta E_{B-V}$) in bin $i$, up to an opacity conversion that is not known per layer. PySM model d12 already implements a six-layer phenomenological version of this idea (Martínez-Solaeche, Karakci & Delabrouille 2018); it is a high-complexity simulation tool, not a data-driven Gaia product. The 3D extinction map replaces those phenomenological layers with data-driven ones, with posterior uncertainties, on the local volume.

2. A prior on where LOS complexity exists. Pelgrims et al. (2021) needed H I cloud counts as a proxy for multi-cloud sightlines. A 3D extinction map provides the same information in true distance, including clouds that share a radial velocity but sit at different distances. Sightlines with a single dominant cloud can be treated with a one-MBB model; sightlines with two or more comparable columns are the ones that must be given extra SED freedom, or masked, or modelled as a sum. This is the Tassis & Pavlidou selection in distance rather than velocity. Pelgrims already made the H I version quantitative; extinction adds distance and loses velocity separation. The net gain is unquantified and is not automatic.

3. A conversion from extinction to emission only after $T_d$ is assigned. Optical/NIR extinction measures dust mass weighted by the optical extinction cross-section. Submillimetre emission measures dust mass times $\kappa_\nu B_\nu(T_d)$. Extinction is temperature-blind. A 3D radiation-field estimate—or a $T_d$ prior from FIR colour, with CIB and $\beta$–$T$ degeneracy—must be supplied from outside the extinction cube. Residuals between a predicted 353 GHz intensity template per layer and Planck 353 GHz $I$ then constrain remaining opacity and temperature variations.

4. A handle on which moments should be large. Pixels whose distance distribution of $\tau_i$ is unimodal should have small second-order moments; pixels that pierce two clouds should have larger moments, with a geometric prior from the column contrast. The 3D map does not predict the spectral moments without $T_i$ and $\beta_i$.

The essential point is that 3D extinction converts an underdetermined integral into a sum over a modest number of discrete, spatially coherent clouds—the same objects that Tassis & Pavlidou showed can source LOS decorrelation—provided those clouds can be matched to magnetic orientations and provided SED differences are obtained from somewhere else.

3.2 Grain size, composition, shape, and alignment

The grain size distribution $n(a)\,\mathrm{d}a$ does not set $\beta_d$ in the submillimetre. For $a\ll\lambda$ the mass opacity and its slope $\beta$ are size-independent for a given composition and shape. Size enters the CMB problem through three narrower channels: (i) $T_d(a)$ at fixed radiation field, because smaller grains are hotter; (ii) the alignment threshold $a_{\mathrm{align}}$, which high-pass-filters the polarized emissivity; (iii) the optical-to-submillimetre opacity conversion that turns $A_V$ into $\tau_{353}$. The frequency dependence of $p$ and any $\beta_d^{\mathrm{P}}-\beta_d^{\mathrm{I}}$ difference are dielectric-function, composition and shape questions.

The relevant polarized emissivity of a layer is

Grain physics should therefore supply a family of priors and spectral derivatives, not a supposedly known per-cloud $\beta_d$.

Size distribution. The classical MRN distribution $n(a)\propto a^{-3.5}$ from $a_{\min}\sim 5\,\mathrm{nm}$ to $a_{\max}\sim 0.25\,\mu\mathrm{m}$ (Mathis, Rumpl & Nordsieck 1977), and its descendants (Weingartner & Draine 2001; Zubko, Dwek & Arendt 2004; Compiègne et al. 2011; Jones et al. 2017), put most of the dust mass in grains with $a\sim 0.1\,\mu\mathrm{m}$ and most of the number in nanoparticles. Far-infrared and submillimetre emission is dominated by the big-grain population, which is in thermal equilibrium. Very small grains and PAHs are stochastically heated and dominate the mid-infrared; they are at most weakly aligned and contribute negligibly to polarized submillimetre emission. Total-intensity dust at $\nu\lesssim 100\,\mathrm{GHz}$ can include additional, less polarized, components (AME from spinning nanoparticles; possibly magnetic dipole emission from ferromagnetic inclusions).

Alignment is a high-pass filter in grain size. In RAT theory the torque efficiency scales as a broken power law in $\lambda/a$ (Lazarian & Hoang 2007). Grains align stably when they rotate suprathermally, which defines a minimum alignment size that increases in denser, less irradiated gas. Magnetically enhanced RAT alignment (MRAT) with iron inclusions can drive high-$J$ attractors for a broader range of shapes (Hoang & Lazarian 2016). The observational consequence is that $f_{\mathrm{align}}(a)$ is a function of the local radiation field $U$, gas density, and grain magnetic susceptibility. A 3D dust map that also traces density and (crudely) $U$ therefore predicts a 3D alignment prior, not a 3D alignment map.

Two grain scenarios, to be marginalized, not chosen.

• Silicate + carbon, carbons poorly aligned (Draine & Fraisse 2009; Guillet et al. 2018). Two populations contribute to $I$ but only one to $P$, so $p(\nu)$ is the ratio of a polarized MBB to a sum of MBBs. This picture was motivated by the 2015 $\beta_d^{\mathrm{P}}-\beta_d^{\mathrm{I}}$ offset and $p(\nu)$ decline, both of which are weaker in 2018 data.
• Astrodust + PAH (Hensley & Draine 2021, 2023). A common large-grain composition for $a\gtrsim 0.02\,\mu\mathrm{m}$ predicts a nearly frequency-independent FIR polarization fraction and attributes the optical-to-FIR polarization ratio to grain shape (oblate axial ratio $\sim 1.4$). Residual frequency dependence of $p$ then comes from LOS mixing of temperatures and field orientations, not from composition.

A third, more dangerous, variant includes microwave magnetic dipole emission from ferromagnetic inclusions, whose polarization is orthogonal to the electric-dipole thermal emission of the same grains (Draine & Hensley 2013). Hensley & Bull (2018) identified iron grains, together with LOS mixing, as the two most pernicious biases for parametric CMB recovery. Until FIR polarimetry between $\sim 100$ and $850\,\mu\mathrm{m}$ distinguishes these models, a B-mode pipeline must marginalize over at least the astrodust and silicate+carbon families, and treat magnetic dipole emission as a discrete systematic.

Why this still matters for decorrelation.

• A two-population model produces a polarized SED that is not a rescaling of the intensity SED, so using GNILC $\beta_d$ and $T_d$ maps (derived from intensity) to scale polarized templates is a model error. A one-population astrodust model predicts that intensity-derived $(\beta_d,T_d)$ should apply to polarization, up to LOS mixing.
• A broad $n(a)$ at fixed $U$ is already a superposition of MBBs inside a single cloud. McBride et al. (2023) show that this biases recovered $(\hat T_d,\hat\beta_d)$ away from the mean of the distribution, but does not by itself bias CMB $Q,U$ unless the angle also varies. Intra-cloud $T_d(a)$ is a moment prior, not a substitute for LOS angle structure.
• The conversion $\tau_{353}/A_V$ depends on the size distribution, porosity and composition. A constrained $n(a)$ bounds the systematic error on that conversion; it does not eliminate the $\sim 10\%$ scatter in $N_{\mathrm{HI}}/E_{B-V}$ at high latitude (Lenz et al. 2017), which is already comparable to the 15% Pelgrims signal.

In a pipeline, grain physics is a mapping from a small set of competing microphysical prescriptions plus local $(U,n_{\mathrm{H}})$ to priors on $\bigl(T_d,\,\beta_d^{\mathrm{P}},\,p_0,\,\kappa_\nu,\,\text{permitted curvature}\bigr)$, including the difference between intensity and polarized SEDs when the model has one.

3.3 Sky-projected magnetic field, layer by layer

The plane-of-sky field $\mathbf{B}_{\mathrm{POS}}$ sets $\psi$ and, through $\gamma$, the local polarization fraction. A single projected-field map integrated through the Galaxy is insufficient: it is the quantity one wants to explain. What matters is the relative angle of the individual emitting layers. Without per-layer $\psi_c$, even a perfect 3D mass map cannot predict $Q$ and $U$: one would still have to fit two amplitudes per layer, and extra layers would make the fit less, not more, constrained.

Observational tracers of $\mathbf{B}_{\mathrm{POS}}$.

1. Submillimetre polarized emission (Planck 353 GHz). The most direct full-sky tracer of the dust-weighted, LOS-integrated $\mathbf{B}_{\mathrm{POS}}$. It remains the amplitude calibrator of any 3D model. At high latitude it is usable at $\sim 80'$ (GNILC); native-resolution maps are noisy in the faintest B-mode fields.

2. Starlight polarization. Dichroic absorption polarizes background starlight parallel to $\mathbf{B}_{\mathrm{POS}}$. With Gaia distances, a stack of stars along a sightline becomes a tomographic probe of $\psi(d)$ (Panopoulou et al. 2019; Pelgrims et al. 2023, 2024). Tassis & Pavlidou (2015) already proposed this as the correction for LOS decorrelation: identify sightlines with $\Delta\psi\gtrsim 60^\circ$ between distance bins and mask or model them. This is the only currently available distance-resolved measurement of $\mathbf{B}_{\mathrm{POS}}$ in dust. The bottleneck is the surface density of high-accuracy, low-$p$ optical polarimetry at high latitude (PASIPHAE and related surveys). Mehandiratta et al. (2026) find that the ratio of polarized dust emission to starlight polarization, $R_{P/p}$, is $3.67\pm 0.05\,\mathrm{MJy}\,\mathrm{sr}^{-1}$ in a densely sampled intermediate-latitude field, versus the Planck sky-average $5.42\pm 0.05\,\mathrm{MJy}\,\mathrm{sr}^{-1}$, and that the two overlapping clouds in that field appear to have different $R_{P/p}$. Starlight polarization cannot fix a per-layer polarized-emission amplitude with a universal factor.

3. H I filament orientation. Linear structures in narrow H I velocity channels are aligned with $\mathbf{B}_{\mathrm{POS}}$ (Clark et al. 2014, 2015). The Rolling Hough Transform applied to HI4PI and GALFA-H I yields full-sky cubes of Stokes $Q_{\mathrm{HI}}(v),U_{\mathrm{HI}}(v)$ whose velocity integrals correlate with Planck 353 GHz $Q,U$ at Pearson $r>0.75$ and reproduce the dust $E/B$ asymmetry at $\ell\lesssim 120$ (Clark & Hensley 2019). H I coherence across velocity is a predictor of $p_{353}$. This is a 3D (position–position–velocity) magnetic map of the neutral medium. PPV is not physical distance; matching H I components to Gaia distance slices is a cross-identification problem, not a solved mapping.

4. Phase structure, not just PPV. RHT filaments trace the cold neutral medium (CNM). Lei & Clark (2024) find that at $|b|>30^\circ$ and $N_{\mathrm{HI}}<4\times 10^{20}\,\mathrm{cm}^{-2}$, $p_{353}$ correlates with the CNM fraction, and that the WNM-associated dust column has mean polarization fraction $\approx 0.22\,p_{\mathrm{CNM}}$ if the CNM is taken as maximally polarized ($p_{\mathrm{CNM}}\sim 0.2$). Assigning the CNM filament angle to a whole distance layer misassigns the WNM dust in it. An additional $\sim 18\%$ of the H I mass appears as magnetically ordered as the CNM and may be the unstable medium.

5. Faraday rotation and synchrotron polarization. These trace $\mathbf{B}_{\mathrm{LOS}}$ and the cosmic-ray-weighted field. They are complementary constraints on $\gamma$ and on the large-scale Galactic field, not substitutes for the dust-weighted $\mathbf{B}_{\mathrm{POS}}$. The inclination $\gamma$ remains poorly measured.

What $\mathbf{B}_{\mathrm{POS}}$ contributes to dedusting, even without 3D dust or grain sizes.

• A polarization-angle distribution per layer, not an exact map. Given $\psi_i(\hat n)$ for each cloud, the only free polarized parameter per layer is an amplitude (or $p_{0,i}\cos^2\gamma_i$). Frequency-dependent $\psi_{\mathrm{tot}}$ is then computed from the Tassis–Pavlidou sum if the layer SEDs differ. With ISRF-derived $\Delta T\sim 1$–$2\,\mathrm{K}$ that prediction is a fraction of a degree unless composition varies. The geometric prior is still useful: it prevents the fit from using CMB frequencies to discover the angle structure after the fact.
• A data-driven mask or weight for LOS-decorrelated pixels. Pelgrims et al. (2021) demonstrated that an H I-only selection already isolates the pixels in which Planck 353/217 $\Delta\psi$ is anomalous. A B-mode analysis can down-weight those pixels, or give them extra SED parameters, without using the CMB frequencies to discover the problem.
• A geometric model of $p$. H I-derived $p_{\mathrm{HI}}$ correlates with $p_{353}$ at Pearson $r\sim 0.6$ over the full sky (Clark & Hensley 2019)—36% of the variance, not “a large fraction.” Combined with the Lei & Clark CNM/WNM split, this still absorbs the leading geometric cause of spatially variable $p$, leaving grain physics to explain the residual, including $p(\nu)$.
• A prior on $E/B$ and parity-odd spectra, to be used as validation observables, not as rigid constraints on the CMB solution. Filament–field alignment produces $C_\ell^{EE}/C_\ell^{BB}\simeq 2$ and $TE>0$. Misalignment at the $\sim 2^\circ$ level correlates with the observed dust $TB$ (Clark et al. 2021; Cukierman et al. 2023).

The sky-projected field is the ingredient that turns a 3D column model into a 3D polarization model. A 2D Planck field orientation alone cannot untangle the LOS.

4. A combined dedusting strategy

The three ingredients are complementary because they constrain different factors in $\mathcal{P}_\nu=\sum_c A_c s_c(\nu)e^{2i\psi_c}$. They do not, even together, determine $s_c(\nu)$. A concrete pipeline therefore has a probabilistic layer catalogue, a family of emissivity priors, a forward model with bandpass integration, a restricted residual fit in the spin-moment formalism, and a validation stage that is not tuned on the tested frequency pair. Masking, spatial weighting and adaptive model complexity are immediately usable, even before a full layer reconstruction is mature.

Two observational regimes must be separated. The deep-field, $\ell\approx 80$ programme (BK, SPIDER, SO SATs) observes a few hundred square degrees at high southern latitude, where local dust dominates, Gaia maps overlap the relevant volume, and H I IVC/LVC tests already exist. The LiteBIRD, $\ell<10$ programme observes most of the sky, including dust beyond 1.25 kpc and near the plane. The layered local model is built for the first regime and is only a partial prior for the second.

4.1 Assemble compatible inputs

Bring the 3D extinction posterior, H I channel maps, starlight polarization, high-frequency $Q,U$, gas tracers (H I, and CO where needed) and CMB-frequency data to controlled beams and pixelization. Intermediate- and high-velocity H I components remain explicit distance-uncertain layers rather than being forced into a Gaia bin. Integrate all sky models through the actual experimental bandpasses, including colour corrections, beam mismatch, and leakage terms.

4.2 Construct a probabilistic layer catalogue

Cluster distance bins into a small number of dust-column peaks. Associate H I velocity components where the match is supported, while retaining ambiguous and distant H I components as separate latent layers. Propagate extinction-map distance and column uncertainties. Do not claim a unique PPV-to-parsec mapping; sample over association hypotheses.

Use column-weighted complexity $\mathcal{N}_c$ (Panopoulou & Lenz 2020; Halal et al. 2024) and magnetic misalignment to set the local model order:

• $\mathcal{N}_c\lesssim 1.1$ and aligned H I / starlight angles: one MBB.
• Intermediate: one MBB plus a low-order complex moment.
• $\mathcal{N}_c\gtrsim 1.5$ and $\alpha\gtrsim 60^\circ$: two (or three) layers, or a down-weight / mask.

This is spatial adaptation, not a claim that $N=2$–$8$ is universally sufficient.

At $|b|\ge 60^\circ$, start from the Local Bubble one-shell model (Alves 2018; Pelgrims 2020; Skalidis & Pelgrims 2019) and add layers only where Halal-style component counts and H I / starlight tomography demand them.

4.3 Infer per-layer field orientations, with a CNM/WNM split

Use differential starlight polarization where stellar sampling permits it; use H I filament orientation elsewhere. Treat the angle as a distribution. Require the vector sum of the layers to be consistent with high-frequency Planck $Q,U$.

Within each distance or velocity layer, do not assign the CNM RHT angle to the whole column. Split the layer column into a CNM-associated part, aligned with the filament angle and carrying $p\sim p_{\mathrm{CNM}}\sim 0.2\times\cos^2\gamma$, and a WNM-associated part with a disordered-field prior and mean $p\approx 0.22\,p_{\mathrm{CNM}}$ (Lei & Clark 2024). Absorb the poorly known inclination $\gamma$ and intrinsic alignment efficiency into calibrated amplitude nuisance parameters unless Faraday or starlight $p/A_V$ genuinely constrains them.

4.4 Assign emissivity priors, not unique SEDs

For each layer, infer priors on $T_d$, polarized opacity, $\beta_d^{\mathrm{P}}$, alignment efficiency and permitted spectral curvature from the local environment (FIR colour, $A_V$-attenuated ISRF, density) and a small set of competing grain models (astrodust; silicate+carbon; optional magnetic dipole). A 3D ISRF model at $\pm 20\%$ gives $\pm 0.7\,\mathrm{K}$, hence $\sim 1\%$ ratio changes between 353 and 150 GHz—too small to produce Pelgrims’ 15% without composition variation. That is a prior, not a prediction of the Pelgrims amplitude.

4.5 Forward model

with the analogous sum for intensity, using the intensity SED of the adopted grain family (which need not equal the polarized SED). The unresolved term is a low-order complex spin-moment expansion for sub-beam structure, imperfect cloud segmentation and residual grain variability (Chluba et al. 2017; Vacher et al. 2023, 2025). After pivot correction, the first-order moments $\mathcal{W}_1^\beta$ and $\mathcal{W}_1^T$ are purely imaginary and encode the leading spectral rotation of $\psi$. Second-order moments encode variances and $\beta$–$T$ covariances. The layered prior is expected to reduce, not zero, those moments: it should fix the large-scale, cloud-to-cloud contribution to $\mathrm{Im}\,\mathcal{W}_1^{\beta,T}$ on sightlines where both $\Delta s_c$ and $\Delta\psi_c$ are measured, and leave intra-cloud MHD turbulence to the residual expansion.

This forward model can generate frequency-dependent $\psi_{\mathrm{tot}}$ and $\mathcal{R}_\ell^{BB}<1$. Whether it generates them at a scientifically relevant amplitude depends on the SED differences assigned to layers. With extinction-plus-ISRF temperatures it generates $\sim 1\%$ weight shifts. With grain-family differences or FIR-colour $T_d$ it can generate more. With a free 15% cloud-to-cloud SED scatter it can reproduce Pelgrims, but that scatter is then a fitted parameter, not a prediction of the three ingredients.

4.6 Parametric residual fit

The forward model will not be perfect. The residual is fitted to the multi-frequency CMB data with a restricted parameter set, using the 3D construction as a strongly informative prior rather than as a rigid template.

H0. Adaptive masking and weighting. Use the 3D+H I+starlight selection to exclude or deweight pixels that are LOS-decorrelated (Pelgrims target sample; Halal multi-component sample; $\alpha\gtrsim 60^\circ$). Apply this to whatever dust model the cosmology pipeline already uses: a 353 GHz template with a global or piecewise $\beta_d$, or a BK-style parametric multi-frequency likelihood. This is the minimal use of the three ingredients and does not require predicting $T_i$ or $\beta_i$. Tassis & Pavlidou proposed the starlight-polarization version of this mask in 2015. It is available now.

H1. Layered MBB with amplitudes fitted at a reference frequency. Fix $\{\psi_i\}$ (as distributions) from H I/starlight. Fix layer columns from Gaia/H I, up to a calibrated opacity conversion. Put Gaussian priors on $\{\beta_{d,i},T_{d,i}\}$ from the grain family and FIR colour, not delta-function constraints. Fit $\{A_i^P\}$ (and optionally a small $\delta\beta_{d,i}$) to the ensemble of frequency maps, with the vector sum of layers tied to 353 GHz $Q,U$. The CMB is a blackbody added in $(Q,U)$. The number of dust amplitudes per pixel is $N_{\mathrm{lay}}$, but $N_{\mathrm{lay}}$ is 1 on most southern high-latitude pixels with $\mathcal{N}_c<1.5$, and 2–3 only where the 3D map and H I demand it. Degeneracies are controlled because $\psi_i$ is not free.

H2. Layered model plus complex spin-moments. Allow $\mathrm{Im}\,\mathcal{W}_1^{\beta}$ and $\mathrm{Im}\,\mathcal{W}_1^{T}$ (and, if the data support them, second-order moments) on the residual after the layered sum, or per layer on multi-cloud pixels. At LiteBIRD sensitivity the Vacher et al. (2022) result implies that moments in both $\beta_d$ and $T_d$ are required for spatially varying dust; Vacher et al. (2025) imply that the imaginary parts are the ones that source $BB$ decorrelation. The 3D model predicts the expected scale of the cloud-to-cloud contribution to those moments, which can be imposed as a Gaussian prior, preventing the moments from absorbing CMB. Scalar $\omega_1^\beta,\omega_1^T$ applied equally to $Q$ and $U$ are not sufficient.

H3. Harmonic-space hybrid. For $r$ likelihoods of the BICEP/Keck type, replace the single-SED, single-$\Delta_d$ dust model by the cross-spectra of the layered forward model, plus a low-order complex-moment expansion of the residual. The layered model supplies the frequency dependence of $\mathcal{R}_\ell^{BB}$ and of $C_\ell^{BB}/C_\ell^{EE}$ if the layer SEDs differ; the moments absorb isotropic SED distortions from averaging. Deliver the model in the form the cosmology analysis needs: frequency-dependent dust templates, a mixing matrix with prior covariance, constrained moment maps, or predicted multi-frequency dust cross-spectra plus a residual covariance.

Synchrotron is treated in parallel: a power law with spatially varying $\beta_s$, optionally a curvature term, and an independent $\psi_s$. H I-based $\psi_i$ should not be imposed on synchrotron. Joint dust–synchrotron modelling is needed below $\sim 100\,\mathrm{GHz}$.

4.7 Why the combination is more than the sum of parts

Each ingredient removes a specific degeneracy, and each leaves one:

A 3D mass map without $\mathbf{B}$ overfits $Q,U$. A $\mathbf{B}$ map without 3D mass cannot assign SED differences to layers. Grain physics without either map can predict a mean $p(\nu)$ but not its spatial pattern. Together they specify the Tassis–Pavlidou two-cloud (and $N$-cloud) geometry on the real sky. The spectral weights $s_c(\nu)$ remain a prior plus a residual fit.

4.8 Validation, not only subtraction

A physical forward model is more tightly constrained than a phenomenological MBB; both are falsifiable by $\chi^2$, $\mathcal{R}_\ell$, $\Delta\psi(\nu)$ and $P_{\nu_1}/P_{\nu_2}$ tests. That is how Planck Int. L, Pelgrims, and SPIDER’s $3\sigma$ sub-region were found. Required null tests:

• On the Pelgrims control sample (single H I cloud), the layered model must collapse to a single MBB and $\psi(\nu)$ must be constant; any remaining 353/217 $\Delta\psi$ is then a bound on residual systematics and on intra-cloud grain physics. Pelgrims’ statistic is a circular standard deviation of $\Delta\psi$ that is noise-sensitive and shifted between PR3 and SRoll2; pass/fail criteria must be stated per data release and with a noise model.
• On the target sample, the model must predict the sign of $\Delta\psi(353,217)$ from the H I-measured $\alpha$ before looking at 217 GHz. The magnitude is a test of the SED prior: an ISRF-only prior predicts $\ll 15\%$; if the data require 15%, the extra scatter is a fitted composition/environment parameter, not a success of grain-size knowledge.
• The predicted $p_{353}$ from H I coherence plus CNM/WNM plus RAT $p_0\cos^2\gamma$ should match Planck $p_{353}$ with residuals that correlate with tracers of grain alignment, not with the CMB. Halal et al. (2024) is the 3D-column version of this test.
• Cross-spectra $\mathcal{R}_\ell^{BB}(\nu_i,\nu_j)$ of the model, evaluated on the actual experiment bandpasses, should be compared to Planck 217$\times$353 and to BICEP/Keck 95$\times$150$\times$220$\times$Planck 353 without introducing a free $\Delta_d$. A free $\Delta_d$ is reserved for the residual.
• Morphological and parity statistics ($EE$, $BB$, $TE$, $TB$, $EB$, $EE/BB$, filament–field alignment) are validation observables.
• Split-frequency ILC or spectral-imaging tests (Régnier et al. 2024, not Maffei et al. 2024) should return a consistent $r$ on the cleaned maps. The point of those tests is to detect unmodelled LOS decorrelation; a successful 3D model should make the detection vanish.
• Repeat across data splits, bandpass perturbations, grain families and cloud-association hypotheses.

4.9 What this can and cannot do for $r$

No end-to-end demonstration currently shows that a 3D-informed likelihood recovers $r$ without a bias larger than the statistical error when the reconstruction uses only the ancillary data that will actually exist. Present “high-complexity” PySM suites (layered MBBs, d12) capture LOS decorrelation phenomenologically but do not ingest Gaia 3D maps or RAT-based alignment.

Published numbers that do exist, and that bound the problem rather than the proposed solution:

• Vacher et al. (2022): spatially varying $(\beta_d,T_d)$ in a LiteBIRD-like simulation, $f_{\mathrm{sky}}=0.7$, produce $\hat r\simeq 1.25\times 10^{-2}$ ($\sim 21\sigma_r$) under a single MBB; a joint $\beta$–$T$ moment treatment is unbiased with $\sigma_r=8.8\times 10^{-4}$. This is a moment-expansion result, not a 3D-map result.
• Sponseller & Kogut (2022): $\sigma_T\gtrsim 1.6\,\mathrm{K}$ or $\sigma_\beta\gtrsim 0.045$ along the LOS biases $r$ at a level that is small for current limits ($\delta r\sim 10^{-4}$) and relevant at PICO/LiteBIRD.
• McBride et al. (2023): residual $D_\ell^{BB}\sim 10^{-5}\,\mu\mathrm{K}^2$ on $\ell\lesssim 100$ is the scale at which $\sigma(r)\sim 10^{-3}$ experiments cannot be indifferent to dust residuals; CMB $Q,U$ in a pixel-domain MBB fit remain unbiased unless both SED and angle vary along the LOS.
• Two-cloud arithmetic (this essay, §2.1): realistic ISRF $\Delta T$ produces $1$–$3\%$ weight shifts and sub-degree to few-degree $\Delta\psi$ except on nearly orthogonal, nearly cancelling sightlines, which are depolarized. In a BK-like field where 150 GHz dust is $\sim 10\times$ the $r=0.01$ spectrum, a 1% coherent residual is $r\sim 10^{-3}$; a 15% residual on a small selected sky fraction, diluted by the control-sample majority, can sit at or below that level. This is an order-of-magnitude sketch, not a forecast of $\Delta r$.

The 3D approach could help through three distinct routes, none of which is yet quantified as a $\Delta r$:

1. Bias reduction on the dangerous pixels, by replacing a false rank-one mixing law with the few layer SEDs actually supported by ancillary data, specifically where $Q$ and $U$ have different effective scalings.
2. Variance control, by fixing or narrowing layer columns and angles so the CMB fit need not estimate an independent amplitude, angle, temperature and index for every layer. That is the principal advantage over an unconstrained multi-MBB model.
3. Spatial adaptation, so that most pixels do not carry maximum foreground complexity.

Success should be reported as residual $BB$, $\Delta r$, $\sigma(r)$, and robustness across grain and cloud-association models, from blinded injections, separately for a BK/SO-like field and for a LiteBIRD-like full-sky analysis. Until that campaign exists, the honest claim is architectural, not numerical.

5. Limitations, missing data, and residual systematics

A layered physical model reduces, but does not eliminate, foreground bias on $r$. The leading limitations are as follows.

The three ingredients do not measure per-layer $T_i$ and $\beta_i$. Extinction is temperature-blind. H I gives angle, not SED. Grain physics gives a global family of dielectric functions. FIR colour and an ISRF model give $\Delta T\sim 1$–$2\,\mathrm{K}$ with CIB and $\beta$–$T$ degeneracy. That is the quantity that generates Tassis–Pavlidou rotation, and it is the quantity the proposal does not supply as data. A 3D-informed pipeline that “predicts rather than fits” $\psi(\nu)$ is predicting a sub-degree effect unless it is allowed additional SED freedom.

3D maps trace extinction, not submillimetre emission. The conversion $\delta A_V\to\tau_{353}$ assumes a dust-to-extinction ratio and a size distribution that may vary with environment. High-latitude sightlines are mostly diffuse, where Lenz et al. (2017) find $\sim 10\%$ scatter in $N_{\mathrm{HI}}/E_{B-V}$; 10% layer-to-layer opacity errors, if correlated with $T_d$ or $\psi$, are enough to leave a residual at the $r\sim 10^{-3}$ level. CIB anisotropies further contaminate intensity-based $\tau_{353}$ and $T_d$ maps at high latitude. Polarized CIB is expected to be small but is not measured at the required level.

Distance completeness and low-column fidelity. Gaia-based maps reach $\sim 1.25\,\mathrm{kpc}$ at high resolution and several kpc at lower resolution. That covers most high-latitude dust in deep southern fields, but not LiteBIRD’s $\ell<10$ sky, and not IVCs/HVCs that can contribute polarized emission at high $|b|$ and are poorly sampled by stars. At BK/SPIDER columns the per-star extinction error is comparable to the total column; layer splitting is prior-dominated. A comparison of integrated Edenhofer $E(B-V)$ against H I-based $E(B-V)$ in the BK and SPIDER fields, with the prior contribution reported, is a required test that this essay does not perform.

H I velocity is not distance. At high latitude Galactic rotation does not provide a unique kinematic distance. Errors in PPV-to-distance matching assign the wrong $\psi$ to the wrong SED, which can create a spurious $\psi(\nu)$ as easily as it corrects a real one. Starlight polarization at known distances is the only currently available external check. No published failure rate of a joint Gaia–H I association algorithm exists; the matching must be sampled as a systematic.

CNM versus WNM. Assigning RHT angles to whole layers misassigns WNM dust, which is a large fraction of the high-latitude column and is substantially less polarized (Lei & Clark 2024).

The inclination $\gamma$ is poorly measured. Polarized emission constrains $\mathbf{B}_{\mathrm{POS}}$, not $\mathbf{B}_{\mathrm{LOS}}$. Wrong $\gamma_i$ mis-weight the layers in $Q,U$ and therefore bias the predicted $\psi(\nu)$.

Cloud-to-cloud $R_{P/p}$. Mehandiratta et al. (2026) show that starlight polarization does not convert to polarized emission with a universal factor, even in a four-square-degree field.

Grain-model uncertainty is a systematic on $p(\nu)$ and $\beta_d^{\mathrm{P}}-\beta_d^{\mathrm{I}}$. The 2018 offset is $0.05\pm 0.03$. Astrodust versus silicate+carbon versus magnetic dipole emission make qualitatively different predictions below $\sim 100\,\mathrm{GHz}$. Intrinsic $\beta(\nu)$ curvature of a single composition is degenerate with LOS mixing.

Alignment physics is not a solved 3D field. RAT predictions for $a_{\mathrm{align}}(U,n_{\mathrm{H}})$ are tested statistically, not cloud-by-cloud at high latitude.

Unresolved sub-beam structure still generates moments. LiteBIRD’s 140 GHz beam is $\sim 31'$; even a perfect cloud catalogue does not resolve MHD turbulence inside that beam. The 3D model reduces the moments; it does not set them to zero. Vacher et al. (2022) showed that at LiteBIRD noise, omitting $T_d$ moments biases $r$ even when $\beta_d$ moments are included. That conclusion survives a 3D prior.

Angular-scale mismatch. Planck 353 GHz polarization is usable at $80'$ at high latitude, HI4PI at $16'$, Edenhofer at $14'$, starlight is sparse. Planck Int. L’s decorrelation, where claimed, grows toward $\ell>100$. Every input is at or coarser than the scale where the effect is expected to matter most for deep-field recombination-bump science.

Planck 353 GHz polarization S/N is limited at high latitude. Extrapolating a smoothed 353 GHz template cannot recover small-scale $\psi(\nu)$. Higher-resolution polarized dust data are needed to calibrate layers at SO/BICEP Array resolution.

Instrumental and analysis systematics mimic LOS decorrelation. Bandpass mismatch, $I\to P$ leakage, gain mismatch, and residual CMB in the 353 GHz template all produce frequency-dependent $(Q,U)$ residuals. Pelgrims et al. found that correcting Planck residual systematics increased the significance of LOS decorrelation, which is encouraging, but Guillet et al. (2026) show that SRoll2 and PR4 disagree on whether the residuals look like $\Delta T$ or $\Delta\beta$. A pipeline that applies a 3D correction derived from Planck to a ground-based deep field must re-estimate the correction on that field’s own ancillary data.

Other foregrounds. Polarized synchrotron with spatially varying $\beta_s$ and possible curvature; a possible polarized AME at the percent level; CO and other line polarization; polarized compact sources. A perfect dust model does not yield an unbiased $r$ if synchrotron is wrong. The 3D dust/B-field construction does not automatically help synchrotron.

Statistical cost. Extra layers and moments consume degrees of freedom. Vacher et al. (2022) already saw $\sigma_r$ increase because dust moments correlate with the CMB SED. A 3D model that adds unconstrained amplitudes per layer will degrade $\sigma_r$ more than it reduces bias. If the ancillary-data constraints are too tight and wrong, one reintroduces bias; if they are too loose, one reintroduces degeneracy. That bias–variance trade is the central design problem.

What is missing as data, not as theory.

• Dense, high-latitude optical polarimetry of stars with Gaia distances (PASIPHAE-like), sufficient to measure $\Delta\psi$ between the two or three local dust layers in every deep-field B-mode pixel. Tassis & Pavlidou (2015) argued for tens of stars per CMB beam; the precise PASIPHAE density required per $N_{\mathrm{side}}=256$ pixel is a survey-design calculation that depends on the target $p$ and on the number of layers, and is not derived here.
• A joint 3D reconstruction that uses Gaia extinction, H I cubes, and starlight polarization in one likelihood, rather than matching products ex post, with a documented association-failure rate.
• FIR polarimetry from $\sim 100$ to $850\,\mu\mathrm{m}$ to decide between grain models and to measure $T_d$ in polarization, where CIB contamination is weaker (Hensley & Draine 2023).
• Direct 3D emission maps—tomography of the FIR intensity itself—which do not yet exist.
• A well-characterized $\gamma$ tracer in the neutral medium.
• End-to-end simulations that start from an MHD+dust cube with a realistic $n(a)$ and RAT alignment, observe it with LiteBIRD and SO/BICEP Array bandpasses, reconstruct with only the ancillary data that will actually exist, and report $\Delta r$, $\sigma(r)$ and residual $BB$ for a BK-like field and for a LiteBIRD-like sky, with and without the 3D prior, including the fraction of decorrelation power removed.

Until those exist, the honest status is: the three ingredients specify a well-posed physical geometry for the Tassis & Pavlidou (2015) / Pelgrims et al. (2021) mechanism that makes a single-MBB-per-pixel subtraction fail; they can already be used to flag and down-weight the worst sightlines; they can supply priors on layer columns and angles that reduce the number of free maps relative to a blind multi-MBB or blind moment expansion; they do not, by themselves, predict the SED differences that generate the effect at a scientifically relevant amplitude; they are not a turn-key replacement for moment expansions and multi-frequency residual fits. The right operational posture is to implement the layered construction as an informative prior and as a diagnostic, to keep a low-order complex $\beta_d$–$T_d$ spin-moment expansion for unresolved averaging, and to treat grain-composition, cloud-association and $\gamma$ uncertainties as discrete systematics that are marginalized, not assumed away.

6. Summary

Polarized-dust frequency decorrelation is the statement that the effective MBB parameters and the polarization angle of Galactic dust can be frequency-dependent, because Stokes $I,Q,U$ are integrals of emitters that do not share a common $(T_d,\beta_d,\psi)$. A single modified blackbody per pixel—or a 353 GHz template scaled by a global $\beta_d$, or a BK-style single-SED likelihood—is then a misspecified model. Current experiments do not show a detected $r$-bias from this effect (Planck 2018 XI, BK18 $\Delta_d$, BK XVI, SPIDER). Next-generation analyses targeting $\delta r\sim 10^{-3}$ cannot assume the effect is negligible: published moment-expansion forecasts (Vacher et al. 2022) and LOS-variance thresholds (Sponseller & Kogut 2022) put spatially varying $(\beta_d,T_d)$ and $\sigma_T\sim 1.6\,\mathrm{K}$ in the dangerous range, and Pelgrims et al. (2021) have detected the LOS angle mechanism in a selected high-latitude sample.

LOS angle rotation requires two conditions at once: different layer SEDs and different layer angles. Three-dimensional dust column maps (Gaia extinction cubes, matched probabilistically to H I clouds) replace the integral by a sum over a few distance layers and identify the sightlines where the mechanism can operate. They do not measure layer temperatures. The grain size distribution, together with composition, shape and RAT alignment, determines which grains emit in polarization and supplies a family of priors on $s_c(\nu)$ and on the $A_V\to\tau_{353}$ conversion; it does not set $\beta_d$ in the Rayleigh–Jeans tail. The sky-projected magnetic field, from H I filaments (CNM, with a WNM disorder prior), starlight tomography, and Planck 353 GHz, assigns $\psi_c$ to each layer so that frequency-dependent $\psi_{\mathrm{tot}}$ is computed rather than fitted—at an amplitude set by the SED prior, which is small unless composition varies.

The combined pipeline is a layered, physically parameterized forward model whose remaining freedom is a small set of amplitudes and complex intra-layer spin-moments, fitted to the CMB frequency maps under priors taken from the ancillary 3D data, with masking of the worst sightlines available immediately. Its limitations are structural, not merely practical: extinction is not emission, velocity is not distance, $\gamma$ and grain composition are underconstrained, the Local Bubble is the high-latitude baseline, low-column 3D maps are prior-dominated, and beam averaging still requires moments. Used as a prior and a diagnostic rather than as an exact sky, the same construction is currently the most direct way to attack the specific astrophysical failure mode that frequency decorrelation names, in the deep-field $\ell\approx 80$ regime where the ancillary maps overlap the emitting volume. It is not yet a demonstrated control of $\Delta r$.

What I took from the critic

The critic’s three structural points are right, and they reorganized the essay.

First, the two-cloud rotation was overstated. Tassis & Pavlidou (2015) and direct Planck-function arithmetic give $\Delta\psi\sim 2.5^\circ$ at $\alpha=60^\circ$ for $\Delta T=6\,\mathrm{K}$, not tens of degrees; tens of degrees need $\alpha\sim 85^\circ$ and near-cancelling polarized intensities, where $p\to 0$. I replaced the original claim with that table and with the corresponding ISRF-scale estimate ($\Delta T\sim 1$–$2\,\mathrm{K}$ $\Rightarrow$ $1$–$3\%$ weight shifts). I also translated Pelgrims’ 15% into $(\Delta\beta,\Delta T)$ and stated that an extinction-plus-ISRF model cannot produce it.

Second, the forward model as originally written could not generate the effect it claimed to predict, because none of the three ingredients measures per-layer $T_i$ or $\beta_i$. I recast the ingredients as constraints on $A_c$, on a family of $s_c(\nu)$ priors, and on $\psi_c$, and I stopped claiming that $\psi(\nu)$ is “computed rather than fitted” at a scientifically interesting amplitude without extra SED freedom.

Third, there was no quantitative link to $r$. I still do not have an end-to-end $\Delta r$ for this pipeline—that would be a simulation campaign, not a paragraph—but I now report the published numbers that bound the problem (Vacher et al. 2022; Sponseller & Kogut 2022; McBride et al. 2023), give an order-of-magnitude residual sketch from the two-cloud arithmetic, and state explicitly that the 3D method has not been shown to recover unbiased $r$ at $\delta r\sim 10^{-3}$.

I also adopted the critic’s factual corrections: $\beta_d^{\mathrm{P}}-\beta_d^{\mathrm{I}}=0.05\pm 0.03$ (Planck 2018 XI), not a significant two-population detection; grain size does not set submillimetre $\beta$; McBride et al. (2023) find unbiased CMB $Q,U$ unless both SED and angle vary; polarized moments are complex spin-2 quantities, with imaginary first-order parts driving $BB$ decorrelation (Vacher et al. 2025); BK18 fits a parametric multi-frequency model and uses $\Delta_d$, not a 353 GHz MBB template; d12 is Martínez-Solaeche, Karakci & Delabrouille (2018); Régnier et al. (2024) not Maffei; Planck Int. L’s $\Delta\beta\sim 0.07$ and $\mathcal{R}(150,353)=0.85$ are indicative/illustrative; CMB-S4 was cancelled on 9 July 2025; SPIDER’s non-detection coexists with a $3\sigma$ MBB failure in a sub-region; the Local Bubble (Skalidis & Pelgrims 2019), Halal et al. (2024), Lei & Clark (2024), Mehandiratta et al. (2026), Panopoulou & Lenz cloud counts versus $\mathcal{N}_c$, Guillet et al. (2026), and the low-column failure of Gaia maps in the BK field all belong in the argument. I stated the empirical status of decorrelation as a forecast, in line with Pelgrims’ own caveat and the four power-spectrum non-detections. I separated the LiteBIRD $\ell<10$ regime from the deep-field $\ell\approx 80$ regime.

What I took from the optimist

The optimist supplied the equation that should have opened the original draft: $\mathcal{P}_\nu=\sum_c A_c s_c(\nu)e^{2i\psi_c}$, and the two necessary conditions for LOS angle decorrelation. That is now the organising identity of the essay.

I adopted the broader term polarized-dust frequency decorrelation; the recasting of current products as 3D extinction/column maps; the recasting of grain size as one part of an emissivity-and-alignment prior that also requires composition, shape, radiation field and environment; the insistence on a layer-resolved field, not a 2D Planck angle; the distinction between physical-distance tomography and H I PPV tomography, to be associated probabilistically; the use of complex spin-moments as the safety net, including both $\beta_d$ and $T_d$; the elevation of masking, spatial weighting and adaptive model complexity as immediately usable outcomes; the requirement to integrate through experimental bandpasses; and the instruction to marginalize over grain families and cloud assignments rather than select one microphysical prescription.

I kept, as the optimist recommended, the original’s distinction between 2D spatial variation and genuine LOS mixing; the explanation that sums of MBBs are not another MBB; the complementarity of the three ingredients; the Pelgrims target/control design; the hybrid of a physical forward model and a residual moment expansion; the limitations list; and the concluding posture that the 3D construction is an informative prior and diagnostic, not a turnkey replacement for multifrequency fitting.

The optimist’s three routes to $r$—bias reduction, variance control, spatial adaptation—and the criterion of success (blinded $\Delta r$, $\sigma(r)$, residual $BB$) are now the evaluation section, in place of an unsourced claim that the method controls $r$.

What I rejected, and why

I rejected the critic’s demand for a new d12- or Vacher-2025-style end-to-end forecast of $\Delta r$ and $\sigma_r$ computed inside this essay. That demand is scientifically correct as a referee request on a paper, but producing a number here without a simulation suite would be a fabrication. I also rejected the demand to “demonstrate 3D-map fidelity” by a new Edenhofer-versus-H I comparison in the BK and SPIDER fields, to give a numerical H I-to-distance failure rate, and to quote a PASIPHAE star density per pixel as if it had been derived. Those are real gaps; they are labelled as missing work, not filled with invented figures.

I rejected the critic’s stronger conclusion that the three ingredients are therefore idle because they do not supply $T_i$ and $\beta_i$. They still constrain the other two factors in $\mathcal{P}_\nu$, they still identify the pixels where the dangerous (SED + angle) case can occur, and they still reduce the number of free maps relative to a blind multi-MBB. That is a real, if more modest, role, and it is exactly the role the optimist defended.

I rejected treating Pelgrims’ 15% as likely pure systematics. The conversion to $(\Delta\beta,\Delta T)$ shows that a simple MBB-plus-ISRF model cannot produce it, and that is an important negative result, but the target/control detection of the mechanism survived a change of Planck processing that increased its significance. The right statement is that the amplitude in the selected sample is not predicted by the three ingredients and must be treated as an extra SED-scatter parameter or as a composition variation.

I rejected rebuilding the whole proposal around a one-shell Local Bubble model as the thing to “beat” in a numerical showdown I have not run. I did insert the Local Bubble as the zeroth-order high-latitude baseline and Halal et al. (2024) as the existing 3D-map-versus-Planck test, which is the part of that demand that can be met from the literature.

From the optimist, I rejected only the residual temptation to imply that the 3D programme is already close to a cosmology deliverable. The optimist’s own best line is the one I kept: the defensible claim is that ancillary maps constrain the latent variables responsible for non-factorizable polarized emission, not that they uniquely reconstruct the microwave sky. I also did not drop the original question’s phrase “spectral-index decorrelation” entirely; I subordinated it to the broader term, because that is the language of the prompt and of Planck Int. L, while making clear that it is too narrow.

I did not reject the original draft’s core architecture—layered geometry from 3D columns, emissivity priors from grain physics, angles from $\mathbf{B}_{\mathrm{POS}}$, residual moments, target/control tests. I rejected its overclaims about what that architecture currently predicts and how large the two-cloud rotation is.