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

synthesis · 86.1 KB
Grok 4.6critic GPT-5.6 Soloptimist Fable 5.1

Dedusting CMB B-modes with 3D dust, grain-size 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 ν≳70 GHz\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 become frequency-dependent. That failure is what the CMB literature calls spectral-index (or frequency) decorrelation.

The physical origin is a superposition effect, not a grain-by-grain mystery. The observed Stokes parameters are integrals of emitters that need not share a common temperature TdT_d, opacity index βd\beta_d, grain mix, alignment efficiency, or magnetic-field orientation. Frequency-dependent net angles are not automatic. They require a covariance between SED differences and polarization orientation: different angles with identical SEDs leave the net angle fixed; different SEDs with identical angles change the effective index but do not rotate the polarization. The controlling object is therefore the joint, polarized-emissivity-weighted distribution

f(β,T,ψ,p0,γ,κν∣d,n^),f(\beta,T,\psi,p_0,\gamma,\kappa_\nu\mid d,\hat n),

especially covariances such as Cov(β,e2iψ)\mathrm{Cov}(\beta,e^{2i\psi}) and Cov(T,e2iψ)\mathrm{Cov}(T,e^{2i\psi}).

Three independent observational ingredients—3D maps of the dust column (and, where they exist, of dust temperature), knowledge of the grain size distribution and alignment physics, and maps of the plane-of-sky magnetic field—do not magically identify that joint distribution. They supply independent, CMB-free constraints on its factors: distance-resolved columns, a mapping from local conditions to polarized SEDs, and orientation proxies per layer. Used together they convert an unconstrained 2D extrapolation into a layered forward model with informative priors, a mask of the worst sightlines, and predicted residual moments. They do not, by themselves, make the inverse problem fully determined, and they have not yet been shown end-to-end to deliver unbiased rr at 10−310^{-3}. This essay states that argument operationally: what is being measured, why a one-MBB-per-pixel model is misspecified, how each ingredient contributes on its own, how they combine into a concrete dedusting strategy, and where the residual systematics remain.

Conventions. Stokes parameters (I,Q,U)(I,Q,U) are written in the IAU convention: the electric-vector position angle (EVPA) ψ\psi is measured from north toward east, and

Q=Pcos⁡2ψ,U=Psin⁡2ψ,P≡Q2+U2.Q = P\cos 2\psi,\qquad U = P\sin 2\psi,\qquad P\equiv\sqrt{Q^2+U^2}.

Planck polarization products used for spherical-harmonic analysis are delivered in the COSMO convention, UCOSMO=−UIAUU_{\mathrm{COSMO}}=-U_{\mathrm{IAU}}; that sign must be applied before quoting E/BE/B, TBTB, or EBEB. The observed polarization fraction p=P/Ip=P/I is biased by noise on PP and by the intensity zero level; Planck analyses use debiased estimators. For optically thin thermal dust, the emission EVPA is orthogonal to the plane-of-sky projection of B\mathbf{B}. Starlight polarization is parallel to BPOS\mathbf{B}_{\mathrm{POS}}, and H I filaments are statistically parallel to BPOS\mathbf{B}_{\mathrm{POS}}. Assigning an H I or starlight magnetic angle directly to the emission Stokes equations without a 90∘90^\circ rotation flips the signs of both QQ and UU. Throughout, ψ\psi denotes the emission EVPA unless ψB\psi_B is written explicitly; they differ by 90∘90^\circ. The geometric factor cos⁡2γ\cos^2\gamma, with γ\gamma the angle between B\mathbf{B} and the plane of the sky, is the projection of aspherical grain emission (Lee & Draine 1985). A projected orientation BPOS\mathbf{B}_{\mathrm{POS}} measures ψB\psi_B, not γ\gamma.

The dust intensity of a homogeneous volume element is written as an MBB

Iν=τν0(νν0)βdBν(Td),I_\nu = \tau_{\nu_0}\left(\frac{\nu}{\nu_0}\right)^{\beta_d} B_\nu(T_d),

with Bν(T)B_\nu(T) the Planck function and a conventional reference ν0=353 GHz\nu_0=353\,\mathrm{GHz}. Polarized emission from the same volume, in IAU coordinates, is

Qν+iUν=p0cos⁡2γ  τν0(νν0)βdBν(Td)  e2iψ,Q_\nu+iU_\nu = p_0\cos^2\gamma\;\tau_{\nu_0}\left(\frac{\nu}{\nu_0}\right)^{\beta_d} B_\nu(T_d)\;e^{2i\psi},

where p0p_0 is the intrinsic polarizing efficiency of the aligned grain population. Equivalently, one may normalize at the reference frequency,

Qν+iUν=AP(νν0)βdBν(Td)Bν0(Td) e2iψ,Q_\nu+iU_\nu = A^P\left(\frac{\nu}{\nu_0}\right)^{\beta_d}\frac{B_\nu(T_d)}{B_{\nu_0}(T_d)}\,e^{2i\psi},

with APA^P the polarized intensity at ν0\nu_0. These two amplitude conventions produce different temperature-moment kernels, and they must not be mixed (see §2.3). As soon as several such volume elements occupy the same beam, the effective SED of II, QQ, and UU need not remain an MBB. The net EVPA becomes frequency-dependent only if the layers differ in both SED and ψ\psi.


1. Polarized thermal dust and why B-mode / tensor-to-scalar-ratio 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 rr at a pivot k=0.05 Mpc−1k=0.05\,\mathrm{Mpc}^{-1}. The present observational bound is r0.05<0.036r_{0.05}<0.036 at 95% confidence from BICEP/Keck through the 2018 season (BK18), tightening to r<0.032r<0.032 when Planck PR4, BAO, and CMB lensing are combined (BICEP/Keck Collaboration 2021; Tristram et al. 2022). LiteBIRD targets a total uncertainty δr∼10−3\delta r\sim 10^{-3}, including foreground residuals; CMB-S4 aims at a 5σ5\sigma detection of r=0.003r=0.003 in its deep patch. At these levels the Galactic polarized foreground is not a small correction. In the BICEP/Keck field, Galactic dust B-modes at 150 GHz150\,\mathrm{GHz} already exceed the r=0.01r=0.01 tensor spectrum (BICEP2/Keck Collaboration 2018). Over most of the high-latitude sky the dust B-mode power at 353 GHz353\,\mathrm{GHz} is larger still.

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

βdP=1.59±0.02,βdI=1.51±0.01\beta_d^{\mathrm{P}}=1.59\pm 0.02,\qquad \beta_d^{\mathrm{I}}=1.51\pm 0.01

for a mean colour temperature Td=19.6 KT_d=19.6\,\mathrm{K} on intermediate-latitude patches (Planck Collaboration Int. XXII 2015). The same paper reported that the dust polarization fraction decreases by (21±6)%(21\pm 6)\% from 353353 to 70 GHz70\,\mathrm{GHz}, inconsistent with a single grain population of frequency-independent pp. Planck 2018 XI, using PR3 maps and end-to-end simulations, quotes a slightly lower polarized index, βdP=1.53±0.02\beta_d^{\mathrm{P}}=1.53\pm 0.02 at the same TdT_d, and a polarized-minus-intensity difference

βdP−βdI=0.05±0.03.\beta_d^{\mathrm{P}}-\beta_d^{\mathrm{I}}=0.05\pm 0.03.

That offset is only ∼1.7σ\sim 1.7\sigma: it is a hint, not an established detection, that the grains dominating II are not identical to those dominating QQ and UU. The later, more conservative number is the one that should be carried into a B-mode analysis.

The polarization itself is produced because aspherical grains align with their short axes preferentially parallel to B\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 strongly size-dependent: only grains larger than a critical alignment radius aaligna_{\mathrm{align}} (typically ∼0.05 μm\sim 0.05\,\mu\mathrm{m} in the diffuse ISM) rotate suprathermally and remain aligned against gas randomization. Consequently the polarized submillimetre SED is weighted toward the large-grain tail of the size distribution, whereas total intensity also receives emission from smaller, hotter, poorly aligned grains. That is a population statement. It does not, by itself, fix βdP−βdI\beta_d^{\mathrm{P}}-\beta_d^{\mathrm{I}} or p(ν)p(\nu) on a given sightline: those also depend on composition, optical constants, and the radiation field (see §3.2).

Planck 353 GHz Q,UQ,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 p353max⁡∼0.20p_{353}^{\max}\sim 0.20–0.220.22 at 80′80' (Planck Collaboration XII 2018), a lower bound on the intrinsic grain polarizing efficiency after line-of-sight depolarization. The uncertainty on pmax⁡p_{\max} is dominated by the intensity zero level.
  • An E/BE/B power asymmetry CℓEE/CℓBB≃2C_\ell^{EE}/C_\ell^{BB}\simeq 2 on intermediate scales, together with a positive TETE correlation, both of which arise because filamentary density structures are preferentially aligned with BPOS\mathbf{B}_{\mathrm{POS}} (Planck Collaboration Int. XXXVIII 2016; Clark et al. 2015). Planck 2018 XI found αEE=−2.42±0.02\alpha_{EE}=-2.42\pm 0.02 and αBB=−2.54±0.02\alpha_{BB}=-2.54\pm 0.02 on its largest high-latitude region, and reported a positive dust TBTB.
  • A tight anti-correlation between pp and the polarization-angle dispersion function S\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 could matter for ℓ∼80\ell\sim 80 B-mode analyses, but whose detection is not uniform across data sets (Planck Collaboration Int. L 2017; Planck Collaboration XI 2020; Ritacco et al. 2023; see §2.2).

Component separation for rr is usually done in one of three families: (i) parametric pixel-domain fits (Commander, FGBuster), in which each pixel is an MBB plus a synchrotron power law plus CMB; (ii) internal linear combination / needlet ILC methods, which enforce unit response to a blackbody CMB and minimize variance—ordinary ILC does not impose a dust-null constraint, though constrained or moment-deprojected ILCs add explicit foreground constraints (Remazeilles, Rotti & Chluba 2021); (iii) template multi-frequency likelihoods at the power-spectrum level, as in BICEP/Keck, which scale a 353 GHz dust template with an MBB and a single βd\beta_d (sometimes with a decorrelation parameter Rℓ\mathcal{R}_\ell). All three families become biased if the true dust SED is not a single MBB, or if QQ and UU do not share the same frequency scaling. Spatially varying mixing also increases ordinary-ILC residuals and can induce ILC bias, but the mechanism is variance-weighted leakage of a non-blackbody foreground, not the nulling of one assumed dust SED. That is the operational problem of spectral-index decorrelation.


2. What spectral-index decorrelation means operationally

2.1 Two distinct, observationally mixed phenomena—and the covariance that couples them

The phrase “spectral-index decorrelation” is used in two related but not identical ways.

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

RℓBB(ν1,ν2)≡CℓBB(ν1×ν2)CℓBB(ν1) CℓBB(ν2).\mathcal{R}_\ell^{BB}(\nu_1,\nu_2) \equiv \frac{C_\ell^{BB}(\nu_1\times\nu_2)}{\sqrt{C_\ell^{BB}(\nu_1)\,C_\ell^{BB}(\nu_2)}}.

A single MBB with uniform βd\beta_d implies Rℓ=1\mathcal{R}_\ell=1. Under additional assumptions—Gaussian βd\beta_d fluctuations that are statistically independent of the amplitude field, and a βd\beta_d map that is coherent on the scales contributing to a given ℓ\ell—one obtains the familiar expansion

Rℓ(ν1,ν2)≈1−12(Δβd)2[ln⁡(ν1/ν2)]2\mathcal{R}_\ell(\nu_1,\nu_2)\approx 1 - \tfrac12(\Delta\beta_d)^2\left[\ln(\nu_1/\nu_2)\right]^2

(Planck Collaboration Int. L 2017). That formula is not generally applicable. Multiplication by a spatially varying βd\beta_d field couples harmonic modes; Δβrms\Delta\beta_{\mathrm{rms}} alone does not determine Rℓ\mathcal{R}_\ell. Spatial correlations between βd\beta_d and the dust amplitude, especially at higher orders, can bias rr even when a low-order moment expansion is fitted (Liu et al. 2025, in a Simons Observatory forecast). Those correlations are physical: colder, denser, higher-column layers have different temperatures and indices. A random-field moment model treats them as nuisance variance; a layered model with per-layer physics can encode them, if the layers are correctly assigned.

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 emission angles, the Stokes parameters add as complex amplitudes:

Qν+iUν=∑iPi(ν) e2iψi.Q_\nu+iU_\nu=\sum_i P_i(\nu)\,e^{2i\psi_i}.

Rotation of the net EVPA occurs only when both Pi(ν)/Pj(ν)P_i(\nu)/P_j(\nu) varies with frequency and ψi≠ψj\psi_i\neq\psi_j. A third, often omitted, condition controls the size of the rotation: the polarized amplitudes must be comparable. A dominant cloud produces little rotation regardless of the second cloud’s angle; large rotations, including the extreme ψtot\psi_{\mathrm{tot}} flip when ρν=P1/P2\rho_\nu=P_1/P_2 crosses unity at α→90∘\alpha\to 90^\circ, live near cancellation. This is the Tassis & Pavlidou (2015) mechanism. It is not equivalent to a spatially varying but frequency-independent βd\beta_d per pixel. A pixel-wise MBB fit still assumes that QQ and UU share a common SED, i.e. that ψ\psi is independent of ν\nu. LOS mixing of misaligned, spectrally distinct clouds violates that assumption: the effective spectral indices of QQ and of UU differ. Parametric component separation that fits a single (βd,Td)(\beta_d,T_d) to (Q,U)(Q,U) jointly therefore leaves a frequency-dependent residual whose E/BE/B content is not that of the CMB.

In the spin-moment language of Vacher et al. (2023a, A&A 669, A5), the net field carries a complex spectral index whose imaginary part is the frequency derivative of the angle:

ψ(ν)≃ψ(ν0)+12 Im(Δβ) ln⁡(ν/ν0).\psi(\nu)\simeq\psi(\nu_0)+\tfrac12\,\mathrm{Im}(\Delta\beta)\,\ln(\nu/\nu_0).

A single MBB with a fixed angle cannot represent this. Neither can a real-valued intensity moment expansion. Scalar moments in β\beta and TT capture SED distortions of II (and of PP if the angle is fixed); they do not represent residual frequency-dependent polarization angles. Polarized residual moments must be spin-2, generally complex quantities.

An illustrative two-cloud calculation is useful for scaling, not as a typical-sky claim. Two clouds with equal polarized intensity at 353 GHz, temperatures 1717 and 21 K21\,\mathrm{K}, indices 1.501.50 and 1.651.65, and emission angles misaligned by 60∘60^\circ produce a net-angle rotation of a few degrees from 353 to 150 GHz and a naive residual of a few percent of the dust power after a 353 GHz template is scaled by a single MBB. Frequency-dependent rotation also leaks the brighter dust EE modes into BB, so the residual BBBB can exceed that naive fraction. Poh & Dodelson (2017) found that temperature-only extrapolation noise, using Planck-compatible temperature distributions, is more than an order of magnitude below then-current errors: possible extreme configurations are not the expected sky average. The B-mode question is whether a subset of mixed, misaligned sightlines, with realistic occurrence rates, can still bias rr at 10−310^{-3}. That is a question about the tails of f(β,T,ψ∣n^)f(\beta,T,\psi\mid\hat n), not about the mean SED.

2.2 Observational status: detected in a selected sample, not detected as a sky-average

The empirical status must be stated sharply, because “decorrelation is a problem for rr” and “decorrelation has been measured on the B-mode sky” are different claims.

Planck large-region analyses do not detect a loss of correlation. Planck Collaboration Int. L (2017) reported RℓBB(217,353)\mathcal{R}_\ell^{BB}(217,353) departing from unity at ℓ≳50\ell\gtrsim 50 in some high-latitude regions, with an indicative ΔβdP∼0.07\Delta\beta_d^{\mathrm{P}}\sim 0.07. Planck 2018 XI, using PR3 maps and end-to-end simulations, found that systematics and noise had overstated that evidence, in agreement with Sheehy & Slosar (2018). For the multipole range 50≤ℓ≤16050\le\ell\le 160, the multi-frequency MCMC fit of Planck XI gives RℓBB(217,353)=0.989±0.005\mathcal{R}_\ell^{BB}(217,353)=0.989\pm 0.005 on the largest region (LR71, fskyeff=71%f_{\mathrm{sky}}^{\mathrm{eff}}=71\%), with an end-to-end lower limit of 0.9910.991; the two-frequency 217-versus-353 analysis is weaker (E2E lower limit 0.9800.980 on LR71, 0.7560.756 on LR24). Planck XI explicitly cautions that its sensitivity “precludes identifying how difficult the component-separation problem will be for more ambitious experiments targeting lower limits on rr”. If the LR71 limit applied to sub-orbital patches, decorrelation might not be a problem at r≃0.01r\simeq 0.01; it does not follow that it is negligible at r∼10−3r\sim 10^{-3}.

BICEP/Keck finds no decorrelation in its deep field. BK18’s decorrelation parameter is consistent with unity, with a statistical error σ(r)=0.009\sigma(r)=0.009. BICEP/Keck XVI (2023) correlated BICEP/Keck and Planck maps with H I morphology templates, detected filamentary polarized dust down to 95 GHz, measured β=1.54±0.13\beta=1.54\pm 0.13 for the filamentary component, and found no evidence for decorrelation between the filaments and the rest of the dust field, nor from including intermediate-velocity H I. The H I morphology template correlates in BB modes at a ∼10\sim 10–65%65\% level over 20<ℓ<20020<\ell<200 with maps that also contain CMB and noise—not a 75% per-layer reconstruction of the dust BB modes used for rr.

SPIDER finds no significant full-region LOS-decorrelation detection, using a Gaussian-mixture test, yet measures spectral indices in two sub-regions that differ at 3.9σ3.9\sigma, and a joint fit with Planck that rejects a simple MBB in one sub-region (SPIDER Collaboration 2025). Spatial SED variation is present; a full-patch LOS-angle effect is not detected.

Pelgrims et al. (2021) detected the LOS mechanism in a selected sample. Using H I to count clouds per sightline (Panopoulou & Lenz 2020) and H I filament orientation to estimate BPOS\mathbf{B}_{\mathrm{POS}} per cloud (Clark & Hensley 2019), they selected a target sample of high-latitude pixels with Nc≳1.5N_c\gtrsim 1.5 and cloud–cloud misalignment ≳60∘\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. The probability that the effect is due to noise alone ranges from 5×10−25\times 10^{-2} to 4×10−74\times 10^{-7}, depending on the CMB map and the treatment of residual systematics; correcting residual systematics in SRoll2 increases the significance. The effect is stronger for more misaligned fields, as the two-cloud model predicts. An intrinsic cloud-to-cloud variation of ∼15%\sim 15\% in the 353/217 polarized-intensity ratio is sufficient to reproduce the measured effect. That is all the paper claims. It does not recover individual-cloud TdT_d, βd\beta_d, grain sizes, or 150 GHz residuals, and it does not convert 15% in a band ratio into a 10% “SED shape” change.

Pixel-level and residual-map studies leave room at high latitude. Osumi et al. (2021) find that Planck polarization at ≳3.7∘\gtrsim 3.7^\circ scales is consistent with a uniform βd=1.55\beta_d=1.55 in ∼20%\sim 20\%-of-sky template fits, but that pixel-level multi-frequency fits still allow a Gaussian σβd\sigma_{\beta_d} up to 0.450.45 at the highest latitudes, 0.300.30 at mid-latitudes, and 0.150.15 at low latitudes. Ritacco et al. (2023) find effective angle residuals of order a degree at 217 GHz and several degrees at 100 GHz at the largest scales after mean-SED subtraction. Guillet, Vacher et al. (2026) analyse residual-map covariances in Planck polarization with spin moments: for SRoll2 the spectral dependence is compatible with temperature variations alone rather than index variations; for PR4 it is incompatible with both simple models. Temperature per layer is the quantity a 3D density plus radiation-field model most directly predicts; that does not make temperature the unique residual in every data release.

The CMB data therefore bound sky-average decorrelation at the few-percent level and cannot themselves reach the sub-percent control next-generation experiments need. Twenty-fold tighter control has to come from somewhere else—or from more frequencies, more moments, and more sky. The three ingredients discussed here are a candidate for that “somewhere else”, not a demonstration that it has been achieved.

2.3 Why a single MBB per pixel fails—and when it does not bias CMB Q,UQ,U

The MBB is a nonlinear function of βd\beta_d and TdT_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 (βˉd,Tˉd)(\bar\beta_d,\bar T_d). The kernels depend on the amplitude convention. If the fundamental SED is

Iν=A(νν0)βˉBν(Tˉ),I_\nu=A\left(\frac{\nu}{\nu_0}\right)^{\bar\beta}B_\nu(\bar T),

the first temperature derivative is proportional to Θ(ν,Tˉ)\Theta(\nu,\bar T), with

Θ(ν,T)=x exT(ex−1),x=hνkT.\Theta(\nu,T)=\frac{x\,e^{x}}{T(e^{x}-1)},\qquad x=\frac{h\nu}{kT}.

If instead AA is a reference-frequency intensity and the SED contains Bν(T)/Bν0(T)B_\nu(T)/B_{\nu_0}(T), the first temperature derivative involves Θ(ν,Tˉ)−Θ(ν0,Tˉ)\Theta(\nu,\bar T)-\Theta(\nu_0,\bar T). Mixing these conventions produces an internally inconsistent expansion. In the Rayleigh–Jeans regime relevant to CMB channels, βd\beta_d and TdT_d are highly degenerate, so a fit that expands only in βd\beta_d can absorb some TdT_d variance—but not enough at LiteBIRD sensitivity (Vacher et al. 2022).

For polarization the expansion must be spin-2 (Vacher et al. 2023a). Each spin moment has real and imaginary components; the first-order complex moment already describes frequency rotation of ψ\psi. Vacher et al. (2023b) showed that this induces a frequency dependence of the dust E/BE/B ratio and of EBEB. Vacher, Carones et al. (2025) showed that second-order spin moments reproduce the six-layer PySM d12 model to a few parts in 10410^{4} across LiteBIRD bands, and that parametric and ILC methods respond to different moment properties. A layered forward model should therefore be delivered to likelihoods as predicted moment maps with covariances, not only as band templates.

Operational consequences for component separation:

  1. Template methods. Scaling a 353 GHz (Q,U)(Q,U) template by a single MBB (or even by a spatially varying βd\beta_d map derived from intensity) mis-predicts the dust at 90–150 GHz wherever ψ(ν)\psi(\nu) or the effective βd\beta_d of polarization differs from that of intensity. The residual is a B-mode field that is spatially correlated with the Galaxy and spectrally intermediate between dust and CMB. Frequency-dependent angle residuals leak dust EE into BB.

  2. Pixel-domain parametric fits. A per-pixel MBB has three dust parameters (AQ,AU,βd)(A_Q,A_U,\beta_d) if TdT_d is fixed, or four if TdT_d is free, plus CMB (Q,U)(Q,U). With ∼10\sim 10 frequency bands this looks overconstrained. It is not: the model is misspecified. McBride, Bull & Hensley (2023) showed that when the true LOS is a distribution of MBBs, the best-fit (β^d,T^d)(\hat\beta_d,\hat T_d) is biased away from both the mean and the median of the true distribution—so the physical interpretation of fitted dust parameters is generally wrong. Recovered CMB Q,UQ,U, however, remain close to unbiased unless SED and polarization angle vary together along the LOS. Their intermediate example gave only ∼0.12σ\sim 0.12\sigma CMB-amplitude bias; the large (≳10σ\gtrsim 10\sigma) bias was an illustrative turbulent-angle model. That distinction is the whole point of bringing in BPOS\mathbf{B}_{\mathrm{POS}} per layer. 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ℓBB(νi×νj)C_\ell^{BB}(\nu_i\times\nu_j) with a single βd\beta_d and a single decorrelation parameter Rℓ\mathcal{R}_\ell (Vansyngel et al. 2017) can absorb a mean loss of coherence but cannot represent frequency-dependent ψ\psi or the distinct β\beta and TT moment spectra. Mangilli et al. (2021) showed that ignoring SED distortions, or replacing them by one Rℓ\mathcal{R}_\ell, biases rr above the sensitivity of next-generation experiments. Vacher et al. (2022) found that LiteBIRD-like simulations with spatially varying (βd,Td)(\beta_d,T_d) produce a ∼21σr\sim 21\sigma_r bias if a single MBB is assumed; expanding in βd\beta_d alone is still insufficient; expanding in both βd\beta_d and TdT_d recovers an unbiased rr with σr=8.8×10−4\sigma_r=8.8\times 10^{-4}, at the cost of moment–CMB correlations that inflate the error. Priors on the moment maps from a 3D forward model are a candidate way to shrink that inflation; they have not been demonstrated.

  4. Blind ILC. Ordinary ILC enforces unit CMB response and minimizes variance. Frequency-dependent mixing leaves a dust residual that can leak into the CMB solution. 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. Predicted spin-moment maps from a layered model are the natural templates to deproject.

In short: spectral-index decorrelation means that the dust sky at 150 GHz is not in general 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 rr. The rest of this essay is about reconstructing the missing layers rather than fitting their average—while remaining explicit about what the ancillary data actually measure.


3. How each ingredient helps on its own

3.1 3D maps of thermal dust emission (and the 3D extinction and temperature maps that currently exist)

A 2D dust map, however high its S/N, is a column integral. Decorrelation is generated by differences among the emitters that were integrated, coupled to differences in orientation. A 3D map converts that integral into a sum of discrete layers, each of which can be assigned its own SED and its own B\mathbf{B}. That is the information gain. It is not a measurement of layer-resolved polarized amplitudes or of βdP\beta_d^{\mathrm{P}}.

What actually exists today are primarily 3D extinction / reddening maps, plus one large-scale 3D dust-temperature reconstruction, not 3D maps of polarized submillimetre emission. They are nonetheless the right starting point, because 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, Gaussian-process prior on the dust density.
  • 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), reconstructed with variational inference (MGVI) on a spherical grid from 69 to 1250 pc in 516 logarithmically spaced distance bins, angular resolution up to Nside=256N_{\mathrm{side}}=256 (∼14′\sim 14'). This map resolves the internal structure of hundreds of local molecular clouds at roughly the angular scale of CMB B-mode analyses. The 516 radial bins are grid sampling of a regularized posterior, not 516 independent measurements along every LOS; the released posterior has strongly correlated spatial uncertainties. Low-extinction high-latitude directions are not automatically the best-resolved layer decompositions merely because much of their dust is local.
  • Zhang & Green (2025) mapped the extinction-curve parameter RVR_V in 3D from Gaia XP spectra of ∼130\sim 130 million stars. RVR_V varies with extinction in a way consistent with accretion-dominated growth in moderately dense gas and coagulation at higher density. That is a 3D grain-size proxy, not a microwave emissivity map. Through the Serkowski relation the optical extinction curve also tracks the size of the aligned grains, which links to the RAT alignment threshold—again as a prior, not as a per-layer βdP\beta_d^{\mathrm{P}}.
  • Zelko, Finkbeiner, Lee & Green (2023) assigned a temperature and index to each Bayestar19 voxel, integrated to match Planck and IRAS emission, and produced the first large-scale 3D dust-temperature map (at 110′110', 55′55', and 27′27'). They resolved distinct temperatures for two Cepheus components on one sightline. This is the prototype of the “3D emission map” the question asks about. It also constrains the optical-depth-to-reddening conversion per region, which they find varies by a factor of ∼2\sim 2, contrary to the constant conversion assumed in many emission-based reddening maps. The reconstruction is constrained by integrated FIR emission: layer temperatures are inferred, not independently measured, and they inherit the usual β\beta–TT degeneracy.

At high Galactic latitude the dust that dominates 353 GHz emission lies mostly within a few hundred parsecs—inside the volume of these maps. That is a crucial practical point: B-mode fields are high-latitude fields, so the 3D dust maps are most complete precisely where CMB experiments need them. Completeness is not the same as layer identifiability. Intermediate- and high-velocity clouds can still contribute polarized emission at high ∣b∣|b| and are poorly sampled by stars.

The Local Bubble is a candidate for the dominant first layer, not a solved geometry. Skalidis & Pelgrims (2019) used stars as distance candles and found the 353 GHz polarized sky at high southern and northern latitudes dominated by a magnetized structure at roughly the Local Bubble shell distance. O’Neill et al. (2025) built a 3D field model on the shell from Planck 353 GHz and a 3D-dust-derived shell geometry, under the assumptions that the shell dominates the polarized signal along most lines of sight and that the field is tangent to the dust-traced surface; they found the Bubble likely dominates the polarized signal along most sightlines, and they validated the reconstruction on a simulated superbubble. Maconi et al. (2023) show from MHD that such a shell can act as a translucent polarization filter for emission behind it. Halal, Clark & Tahani (2024), using Edenhofer et al. (2024) and Planck 353 GHz, found no imprint of the tangential-field assumption on the polarization fraction. They did find that, at fixed column, the polarization fraction decreases when the dust column is distributed among multiple components at different distances—but only when dust out to 1.25 kpc1.25\,\mathrm{kpc} is included; the effect disappears if the map is truncated at the shell distance (∼270 pc\sim 270\,\mathrm{pc}). The extended 3D distribution, not solely the Bubble, sets p353p_{353} at 80′80'. The Bubble is therefore a physically motivated first layer whose orientation can be given a geometric prior; its inclination and its dominance are still fitted quantities.

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

  1. A physically motivated layering. Instead of one MBB per pixel one writes
Iν(n^)=∑i=1Nlay(n^)τi(n^)(νν0)βd,iBν(Td,i),I_\nu(\hat n)=\sum_{i=1}^{N_{\mathrm{lay}}(\hat n)} \tau_i(\hat n)\left(\frac{\nu}{\nu_0}\right)^{\beta_{d,i}} B_\nu(T_{d,i}),

with τi\tau_i related to the differential extinction δAV,i\delta A_{V,i} in distance bin ii. The conversion is not a universal scale factor. If δAV,i\delta A_{V,i} is in magnitudes, the V-band optical depth is τV,i=δAV,i/1.086\tau_{V,i}=\delta A_{V,i}/1.086; going from τV\tau_V to τ353\tau_{353} requires the ratio of submillimetre opacity to optical extinction, which depends on size distribution, composition, porosity, and RVR_V, plus an explicit dust-to-gas convention. An overall scale can be calibrated on the high-latitude τ353\tau_{353}–AVA_V correlation and on the NHIN_{\mathrm{HI}}–EB−VE_{B-V} relation of Lenz et al. (2017), which holds for NHI<4×1020 cm−2N_{\mathrm{HI}}<4\times 10^{20}\,\mathrm{cm}^{-2} with ∼10%\sim 10\% scatter. Layer-to-layer variations of that conversion, if correlated with TdT_d or ψ\psi, are a residual at the r∼10−3r\sim 10^{-3} level. PySM model d12 already implements a six-layer phenomenological version of this idea and is used as the “high-complexity” dust model for CMB-S4 / LiteBIRD / SO forecasts. The 3D extinction map replaces those phenomenological layers with data-driven ones, subject to the conversion systematic.

  1. 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 related information in true distance, including clouds that share a radial velocity (and would be blended in H I) but sit at different distances. It does not provide “the same information” as an H I cloud count: H I also carries velocity-coherent magnetic orientation, and 3D extinction does not. 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 of MBBs. Unimodal extinction versus distance does not imply small ω2β\omega_2^\beta or temperature moments: a single cloud may contain illumination gradients, composition variation, and beam-scale mixing. Conversely, multiple distance peaks with identical polarized SEDs and aligned fields need not produce decorrelation.

  2. A conversion from extinction to emission, once TdT_d is assigned. Optical/NIR extinction measures dust mass weighted by the optical extinction cross-section. Submillimetre emission measures dust mass times κνBν(Td)\kappa_\nu B_\nu(T_d). Zelko et al. (2023) is the existing implementation of this conversion in 3D. Residuals between the predicted integrated emission and Planck 353 GHz II then constrain remaining opacity and temperature variations, rather than having to infer the entire 3D structure from the SED alone. Polarized CIB is expected to be small but is not measured at the required level; CIB anisotropies further contaminate intensity-based τ353\tau_{353} and TdT_d maps at high latitude.

  3. A handle on which moments should be large. Moment expansion treats all averaging as unknown. A 3D map tells the pipeline where moments are expected to be large: pixels that pierce two clouds with different local densities (hence different TdT_d) and different field orientations. It does not set the moments to zero.

Planck Collaboration Int. XVII (2014) already found, from correlating dust with H I velocity components over the southern cap, that temperature is anti-correlated with emissivity and that an extra component at ∼100 GHz\sim 100\,\mathrm{GHz} is of order a quarter of the emission, attributable to magnetic dipole emission, carbon dust, or flattening silicate opacity. That is direct evidence that layers differ in SED, and that the low-frequency polarized SED is not settled by an MBB.

3.2 Dust grain size distribution

The grain size distribution n(a) dan(a)\,\mathrm{d}a and the size-dependent alignment function falign(a)f_{\mathrm{align}}(a) control three quantities that a single MBB cannot represent: the frequency dependence of pp, a possible difference between βdI\beta_d^{\mathrm{I}} and βdP\beta_d^{\mathrm{P}}, and the coupling between local physical conditions and the SED of each 3D layer. They do not uniquely predict those quantities. Composition, shapes, porosity, aggregation, dielectric functions, iron-inclusion abundance, and their environmental variation are missing degrees of freedom. Hensley & Draine (2023) argue that a single “astrodust” composition for a≳0.02 μma\gtrsim 0.02\,\mu\mathrm{m} explains a polarization fraction flat within about 10% from 250 μm250\,\mu\mathrm{m} to 3 mm and attributes the optical-to-FIR polarization ratio to grain shape (oblate axial ratio ∼1.4\sim 1.4). The existence of that viable one-population model is the demonstration that n(a)n(a) does not uniquely predict βdP−βdI\beta_d^{\mathrm{P}}-\beta_d^{\mathrm{I}} or p(ν)p(\nu).

Size distribution. The classical MRN distribution n(a)∝a−3.5n(a)\propto a^{-3.5} from amin⁡∼5 nma_{\min}\sim 5\,\mathrm{nm} to amax⁡∼0.25 μma_{\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∼0.1 μma\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. They are not a superposition of equilibrium MBBs. For large grains in the long-wavelength Rayleigh regime, βd\beta_d is controlled strongly by material dielectric response, disorder, inclusions, and structure—not by n(a)n(a) alone. The practical CMB-relevant statement is that polarized dust is a large-grain phenomenon, while total-intensity dust at ν≲100 GHz\nu\lesssim 100\,\mathrm{GHz} can include additional, less polarized, components (AME from spinning nanoparticles; possibly magnetic dipole emission from ferromagnetic inclusions; a carbonaceous population with a different βd\beta_d).

Alignment is a high-pass filter in grain size. In RAT theory the torque efficiency scales as a broken power law in λ/a\lambda/a (Lazarian & Hoang 2007). Grains align stably when they rotate suprathermally, which defines a minimum alignment size aaligna_{\mathrm{align}}. The actual approximation (Hoang, Lazarian & Schlickeiser 2021; Hoang & Lazarian 2014, 2016) depends on normalized gas density nHn_{\mathrm{H}}, gas temperature TgasT_{\mathrm{gas}}, mean radiation wavelength λˉ\bar\lambda, radiation energy density UU, radiation anisotropy γrad\gamma_{\mathrm{rad}}, and an IR-damping factor 1+FIR1+F_{\mathrm{IR}}, as well as grain helicity and magnetic susceptibility. Magnetically enhanced RAT alignment (MRAT) with iron inclusions can drive high-JJ attractors for a broader range of shapes. Radiative torque disruption can set amax⁡a_{\max} in strong radiation fields. A 3D extinction map, even with a fiducial cloud thickness to estimate nHn_{\mathrm{H}} and a TdT_d-based estimate of UU, does not supply γrad\gamma_{\mathrm{rad}}, FIRF_{\mathrm{IR}}, or the magnetic properties. Using TdT_d to infer UU while claiming to predict TdT_d is circular. Alignment modelling is becoming statistically testable (Tram and collaborators have reproduced Serkowski and super-Serkowski starlight spectra with RATs plus paramagnetic relaxation under specified radiation fields), but it is not a cloud-by-cloud 3D alignment map at high latitude.

Why this matters for decorrelation.

  • Frequency-dependent pp. If two grain populations contribute to II but only one is aligned, p(ν)p(\nu) is the ratio of a polarized MBB to a sum of (possibly non-MBB) components and is not constant. This is a natural interpretation of the Planck Int. XXII decline of pp from 353 to 70 GHz, and of a possible βdP−βdI\beta_d^{\mathrm{P}}-\beta_d^{\mathrm{I}} offset: silicates (or astrodust) aligned, a carbonaceous or PAH-related component unaligned (Draine & Fraisse 2009; Guillet et al. 2018). The astrodust picture predicts instead that intensity-derived (βd,Td)(\beta_d,T_d) should apply to polarization, and that residual frequency dependence of pp comes from LOS mixing of temperatures and field orientations, not from composition. For component separation the distinction is operational: a two-population model produces a polarized SED that is not a rescaling of the intensity SED, so using GNILC βd\beta_d and TdT_d maps (derived from intensity) to scale polarized templates is a model error. Until FIR polarimetry between ∼100\sim 100 and 850 μm850\,\mu\mathrm{m} distinguishes these models (the diagnostic proposed by Hensley & Draine), a B-mode pipeline must marginalize over at least two grain scenarios, not adopt one. Forthcoming wide-field FIR polarimetry (BLAST Observatory, PRIMA in NASA Phase A, FYST 850 GHz) is what decides composition; Taurus, LiteBIRD (to 402 GHz), Simons Observatory (280 GHz), and CMB-S4 (270 GHz) replace Planck 353 GHz as the layer calibrator at higher S/N.

  • Size-dependent temperature of big grains. Smaller equilibrium grains are hotter at fixed radiation field because they absorb optical/UV light more efficiently relative to their FIR opacity. A broad n(a)n(a) at fixed UU is already a superposition of MBBs, i.e. a source of intra-cloud temperature moments. McBride et al. (2023) showed that fitting a single MBB to a broad TdT_d distribution returns biased (T^d,β^d)(\hat T_d,\hat\beta_d). Knowing n(a)n(a) lets one compute the intra-cloud moment structure as a prior, rather than fitting it blindly. Sponseller & Kogut (2022) found that line-of-sight temperature dispersion above about 1.6 K1.6\,\mathrm{K} biases rr for PICO with a single MBB; a 3D temperature map predicts that dispersion per pixel, with the caveats of Zelko et al. (2023).

  • Iron and magnetic dipole emission. Ferromagnetic inclusions produce microwave magnetic dipole emission 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. A size- and composition-resolved dust model is the only way to predict whether this component is present at a given frequency and whether it rotates ψ(ν)\psi(\nu); the prediction is scenario-dependent and must be marginalized.

  • Connecting 3D extinction to 3D emission. The conversion factor τ353/AV\tau_{353}/A_V depends on the size distribution (optical extinction is dominated by a∼λ/2π∼0.1 μma\sim\lambda/2\pi\sim 0.1\,\mu\mathrm{m} grains; submillimetre opacity has a different size weighting), on porosity, and on RVR_V. Zhang & Green (2025) give a 3D RVR_V prior on that conversion. A constrained n(a)n(a) turns a 3D AVA_V cube into a 3D τ353\tau_{353} cube with a physically bounded systematic error, rather than an unknown scale factor—provided one also specifies composition.

In a pipeline, the grain-size ingredient is therefore not a free SED template. It is a prior mapping

(n(a),falign(a),U,γrad,nH,composition)  ⟶  (Td(a), βd(a), p0(a), κν(a))\bigl(n(a),f_{\mathrm{align}}(a),U,\gamma_{\mathrm{rad}},n_{\mathrm{H}},\text{composition}\bigr) \;\longrightarrow\; \bigl(T_d(a),\,\beta_d(a),\,p_0(a),\,\kappa_\nu(a)\bigr)

that predicts, for each 3D cell, both the intensity SED and the polarized SED, including their difference, under an adopted grain scenario. The scenario is a discrete systematic, not a solved microphysical field.

3.3 Sky-projected magnetic field

The plane-of-sky field BPOS\mathbf{B}_{\mathrm{POS}} sets the emission EVPA ψ=ψB+90∘\psi=\psi_B+90^\circ. It does not measure the inclination γ\gamma, and it does not make the polarized sky a derived quantity. Polarized amplitude still requires emissivity, alignment efficiency, inclination, field tangling, and layer association. Without an orientation proxy, even a perfect 3D mass map cannot predict QQ and UU: one would still have to fit two amplitudes per layer. With an orientation proxy, the number of free polarized parameters per layer drops from two Stokes amplitudes to one (an amplitude, or p0,icos⁡2γip_{0,i}\cos^2\gamma_i), and frequency-dependent ψtot\psi_{\mathrm{tot}} becomes a predicted rather than freely fitted effect—if the proxy is correctly associated with the SED of that layer.

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

  1. Submillimetre polarized emission (Planck 353 GHz). This is the most direct full-sky tracer of the dust-weighted, LOS-integrated BPOS\mathbf{B}_{\mathrm{POS}}. Its limitation is exactly the limitation of 2D dust maps: it is an integral. It is the quantity one wants to explain, not the 3D field itself. It remains the amplitude calibrator of any 3D model. In the faintest B-mode fields those maps are noisy at native resolution; GNILC 80′80' products are used instead, which already average over the small-scale structure that generates small-scale decorrelation.

  2. Starlight polarization. Dichroic absorption polarizes background starlight parallel to BPOS\mathbf{B}_{\mathrm{POS}}. Starlight polarization is cumulative. A star “whose Gaia distance falls in the layer” does not directly give that layer’s angle. Differential tomography requires foreground/background stellar pairs or a joint inversion of stellar q,uq,u, with intrinsic stellar polarization, distance uncertainty, selection effects, and angular sparsity propagated (Pelgrims et al. 2023). Pelgrims et al. (2024) inverted RoboPol polarimetry over four square degrees and produced the first degree-scale starlight-polarization tomography map; most sightlines cross more than one cloud, with the plane-of-sky field orientation changing with distance on most of them. PASIPHAE with the WALOP polarimeters is designed to deliver hundreds of polarized stars per square degree at high latitude. Where stars exist, this is the only distance-resolved measurement of BPOS\mathbf{B}_{\mathrm{POS}} in dust. Until PASIPHAE-scale coverage exists, it is a sparse external check, not a full-sky layer angle map.

  3. H I filament orientation. Linear structures in narrow H I velocity channels are aligned with BPOS\mathbf{B}_{\mathrm{POS}} (Clark et al. 2014, 2015). The Rolling Hough Transform applied to H I4PI and GALFA-H I yields, with no free parameter relating H I intensity to dust emission, full-sky cubes of Stokes QHI(v),UHI(v)Q_{\mathrm{HI}}(v),U_{\mathrm{HI}}(v) whose velocity integrals correlate with Planck 353 GHz Q,UQ,U at r>0.75r>0.75 and reproduce the dust E/BE/B asymmetry at ℓ≲120\ell\lesssim 120 (Clark & Hensley 2019). That result is a position–position–velocity proxy, not a distance-resolved magnetic-field measurement, and it applies to large-scale, velocity-integrated maps. H I coherence across velocity is a predictor of p353p_{353}: disordered orientation in vv–θ\theta space corresponds to LOS field tangling and depolarization (Clark 2018). H I traces only part of the dusty ISM. Optical-depth corrections, CO-dark molecular gas, ionized-gas dust, IVC dust-to-gas variation, and non-filamentary neutral material can all break the assumed H I–dust layer correspondence. Matching H I velocity components to 3D extinction slices is a cross-identification problem, not a solved mapping. Errors in that matching assign the wrong ψ\psi to the wrong SED and can create a spurious ψ(ν)\psi(\nu) as easily as they correct a real one.

  4. Faraday rotation and synchrotron polarization. These trace BLOS\mathbf{B}_{\mathrm{LOS}} and the cosmic-ray-weighted field, respectively, i.e. different phases and a different projection. They are complementary constraints on γ\gamma and on the large-scale Galactic field, not substitutes for the dust-weighted BPOS\mathbf{B}_{\mathrm{POS}}. Faraday structures in LOFAR data at tens of parsecs, aligned with H I filaments, have been reported at the Bubble wall (Erceg et al. 2024); a local field reversal plane within a few hundred parsecs has been placed by Booth et al. (2025). These are optional inclination priors, not a γ\gamma map of the neutral dusty medium.

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

  • A polarization-angle template that can be applied layer by layer. Given ψB,i(n^)\psi_{B,i}(\hat n) for each cloud, the emission angle is ψi=ψB,i+90∘\psi_i=\psi_{B,i}+90^\circ, and the only free polarized parameter per layer is an amplitude (or p0,icos⁡2γip_{0,i}\cos^2\gamma_i). Frequency-dependent ψtot\psi_{\mathrm{tot}} is then computed from the Tassis–Pavlidou sum, not fitted as a nuisance—provided the velocity-to-distance association is correct.

  • 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 in a parametric fit, or give them extra SED parameters, without using the CMB frequencies to discover the problem after the fact. Masking is not automatically an improvement: it reduces foreground bias but also increases sample variance, changes mode coupling, and can create EE-to-BB leakage. Improvement must be demonstrated in mean-squared error or rr-coverage.

  • A geometric model of pp. Because pHIp_{\mathrm{HI}} derived from H I morphology correlates with p353p_{353} at r∼0.6r\sim 0.6 over the full sky (Clark & Hensley 2019), a large fraction of the variance in polarization fraction is geometric (LOS tangling and γ\gamma), not grain-alignment physics. Assigning ψi\psi_i and a geometric depolarization factor per layer therefore absorbs the leading cause of spatially variable pp, leaving grain physics to explain the residual, including p(ν)p(\nu). The factor cos⁡2γi\cos^2\gamma_i remains degenerate with p0,ip_{0,i}. Two practical options, both used in the literature, are: (i) absorb cos⁡2γi\cos^2\gamma_i into a single amplitude AiPA_i^P fitted to 353 GHz Q,UQ,U; (ii) use H I coherence / RHT peak width as a proxy for LOS tangling and γ\gamma, as in Clark (2018). Neither is a precision γ\gamma map.

  • A prior on E/BE/B and parity-odd spectra. Filament–field alignment produces CℓEE/CℓBB≃2C_\ell^{EE}/C_\ell^{BB}\simeq 2 and TE>0TE>0. A list of angles does not automatically generate that asymmetry; it requires the correct joint statistics of density morphology and field orientation. Cukierman, Clark & Halal (2023) found a global misalignment of H I filaments and Planck dust polarization of ∼2∘\sim 2^\circ, scale-independent above the H I4PI beam, that correlates with the observed dust TBTB. A 3D field model that includes a small, spatially coherent misalignment angle can predict, rather than overfit, parity-odd dust spectra that would otherwise leak into cosmic birefringence estimators. Filament population models (Hervías-Caimapo & Huffenberger 2022) can populate sub-beam scales with the correct E/BE/B, TETE, and TBTB statistics given a layer’s field orientation.

The sky-projected field is the ingredient that turns a 3D intensity model into a 3D polarization model. Without it, extra layers add extra (Q,U)(Q,U) amplitudes and the fit becomes less, not more, constrained. With it, the fit is still underdetermined at one frequency if Nlay>2N_{\mathrm{lay}}>2, or if two layers are nearly aligned.


4. A combined dedusting strategy

The three ingredients are complementary because they constrain different factors in the polarized radiative transfer. They do not supply the joint distribution by concatenation. A concrete pipeline has a forward-model stage, a parametric (or moment) residual stage, and a validation stage. The association of extinction distance, H I velocity, and starlight polarization must live in one probabilistic likelihood, not in an ex post match. Bandpass integration, thermodynamic-unit conversion, beams, transfer functions, calibration covariance, and reference-frequency colour corrections belong in the generative model; a monochromatic Qν,UνQ_\nu,U_\nu equation is not an r∼10−3r\sim 10^{-3} likelihood.

4.1 What is actually observed per layer

ObservableWhat it measuresWhat it does not measure
Gaia / 3D extinctionCumulative AV(d)A_V(d) to sampled stars; a regularized differential-extinction posteriorLayer TdT_d, βdP\beta_d^{\mathrm{P}}, κ353\kappa_{353}, p0p_0, γ\gamma, polarized amplitude
Zelko-like 3D TdT_dVoxel temperatures constrained by integrated FIR IνI_\nuPolarized SED; layer-resolved βdP\beta_d^{\mathrm{P}}; uniqueness under the β\beta–TT degeneracy
Zhang–Green 3D RVR_VOptical/NIR extinction-curve slopeMicrowave opacity and alignment efficiency
H I4PI / RHTBrightness and morphology in vLSRv_{\mathrm{LSR}}; a BPOS\mathbf{B}_{\mathrm{POS}} orientation proxy for some neutral structuresDistance; ionized and molecular dust; γ\gamma; layer amplitude
Starlight q,uq,uCumulative dichroic absorption to the star, parallel to BPOS\mathbf{B}_{\mathrm{POS}}Layer-resolved angle without a differential inversion; emission EVPA without +90∘+90^\circ
Planck 353 GHz Q,UQ,UDust-weighted, LOS-integrated polarized emission at one frequencyAny layer-resolved quantity; two Stokes observables cannot identify Nlay>1N_{\mathrm{lay}}>1 amplitudes

Layer-resolved TdT_d, βdP\beta_d^{\mathrm{P}}, κν\kappa_\nu, γ\gamma, and polarized amplitude are latent. The forward model may predict them under priors; it may not treat them as data.

4.2 Forward model

Discretize each LOS into NN layers. In practice N∼2N\sim 2–44 is the most that high-latitude data can support: Edenhofer-like distance bins should be clustered into magnetically coherent clouds using H I velocity–orientation structure (as in the Halal et al. 2024 dendrogram), so that one does not fit hundreds of correlated distance slices independently. Treat the Local Bubble shell as the first layer where it dominates, with a geometric prior on orientation from O’Neill et al. (2025), and with γ\gamma still free or absorbed into the amplitude. Treat intermediate-velocity gas as a distinct layer with a distance where stars or a cloud catalogue provide one, and as distance-free otherwise.

Column and distance (3D dust). For layer ii in pixel n^\hat n,

τ353,i(n^)=κ353κV δAV,i(n^)1.086,\tau_{353,i}(\hat n) = \frac{\kappa_{353}}{\kappa_V}\, \frac{\delta A_{V,i}(\hat n)}{1.086},

where κ353/κV\kappa_{353}/\kappa_V is set by n(a)n(a), composition, porosity, and the local RVR_V, and may be given a Zelko et al. (2023) prior. An overall scale is calibrated on the high-latitude τ353\tau_{353}–AVA_V and NHIN_{\mathrm{HI}}–EB−VE_{B-V} relations, not assumed universal per layer.

Local physical state (grain sizes + environment). Associate to each layer a radiation-field strength UiU_i and anisotropy γrad,i\gamma_{\mathrm{rad},i} from a 3D radiation-field model (not from TdT_d if TdT_d is also being predicted), and a density nH,in_{\mathrm{H},i} from δAV,i\delta A_{V,i} plus a thickness prior, or from H I brightness temperature in the matching velocity interval. Then, under an adopted grain scenario:

  • Compute Td(a)T_d(a) for the local radiation field by equating absorption of the ISRF to FIR emission, using the adopted dielectric function (astrodust or silicate+carbon; both are retained).
  • Compute aalign,ia_{\mathrm{align},i} from RAT/MRAT given (Ui,γrad,i,nH,i,Tgas,i,FIR)(U_i,\gamma_{\mathrm{rad},i},n_{\mathrm{H},i},T_{\mathrm{gas},i},F_{\mathrm{IR}}).
  • Integrate over n(a)falign(a)n(a)f_{\mathrm{align}}(a) to obtain layer-averaged Tˉd,i\bar T_{d,i}, βˉd,i\bar\beta_{d,i}, pˉ0,i\bar p_{0,i}, and intra-layer spin moments if the size distribution is broad.

This step is the difference between a phenomenological multi-MBB and a physical multi-MBB. It predicts, rather than fits, βdP−βdI\beta_d^{\mathrm{P}}-\beta_d^{\mathrm{I}} and p(ν)p(\nu) within a layer, up to the discrete grain-scenario systematic.

Orientation (sky-projected B\mathbf{B}). Assign ψB,i(n^)\psi_{B,i}(\hat n) from a joint likelihood of (i) the H I-derived angle in the velocity range matched to layer ii and (ii) differential starlight tomography where stars exist. The emission angle is ψi=ψB,i+90∘\psi_i=\psi_{B,i}+90^\circ. The association of H I velocity with extinction distance is a latent variable with a prior, not a delta function. False cloud splits and mergers, H I optical depth, missing gas phases, and starlight selection are part of that likelihood.

Emission integral, band-integrated. In IAU coordinates, at each experiment band bb with bandpass Wb(ν)W_b(\nu) and in thermodynamic CMB units,

(Q+iU)bd(n^)=∫dν Wb(ν) fCMB(ν)×∑iAiP(n^)(νν0)βd,iPBν(Td,i)Bν0(Td,i) e2iψi(n^),\begin{aligned} (Q+iU)_b^{\mathrm{d}}(\hat n) &=\int\mathrm{d}\nu\,W_b(\nu)\,f_{\mathrm{CMB}}(\nu)\\ &\quad\times\sum_i A_i^P(\hat n) \left(\frac{\nu}{\nu_0}\right)^{\beta_{d,i}^{\mathrm{P}}} \frac{B_\nu(T_{d,i})}{B_{\nu_0}(T_{d,i})} \,e^{2i\psi_i(\hat n)}, \end{aligned}

with an analogous sum for IbdI_b^{\mathrm{d}}, intra-layer spin moments inside each layer, and beams applied in harmonic space. If the aligned population is a subset of the emitting population, βd,iP≠βd,i\beta_{d,i}^{\mathrm{P}}\neq\beta_{d,i} and the polarized amplitude does not scale with Iν,iI_{\nu,i}. That is a prediction of the grain scenario, not an extra free map.

This forward model generates, when the ingredients are correctly associated:

  • frequency-dependent ψtot(n^,ν)\psi_{\mathrm{tot}}(\hat n,\nu) on multi-cloud, misaligned, comparable-amplitude sightlines (Tassis & Pavlidou);
  • RℓBB(ν1,ν2)<1\mathcal{R}_\ell^{BB}(\nu_1,\nu_2)<1 with a specific frequency dependence set by the layer ΔTd\Delta T_d, Δβd\Delta\beta_d, and angle–SED covariance, not by a free Rℓ\mathcal{R}_\ell—or Rℓ=1\mathcal{R}_\ell=1 if layer SEDs are the same or morphologies remain proportional;
  • a difference between intensity-derived and polarization-derived effective spectral indices, if the grain scenario produces one;
  • geometric E/BE/B asymmetry from filament–field alignment, if the joint density–orientation statistics are right.

It does not generate those features “automatically” from a list of columns and angles.

Identifiability. At 353 GHz there are two Stokes observables per pixel. For Nlay=1N_{\mathrm{lay}}=1, a single polarized MBB with fixed TdT_d is identified, up to the usual β\beta–TT degeneracy. For Nlay=2N_{\mathrm{lay}}=2 with known, distinct angles, two amplitudes remain; one frequency is then formally sufficient at the reference band, and additional frequencies constrain SED differences. The mixing matrix becomes ill-conditioned when Δψ→0\Delta\psi\to 0 or when the layer SEDs are nearly identical, and CMB Q,UQ,U remain degenerate with a dust residual that is spectrally close to a blackbody. For Nlay≥3N_{\mathrm{lay}}\ge 3, 353 GHz cannot identify the amplitudes even with known angles; the extra frequencies of a LiteBIRD- or CMB-S4-like experiment help only to the extent that layer SEDs differ. Nearly aligned layers and nearly identical SEDs are the generic high-latitude case, not the exception. The 3D ingredients therefore work as priors that reduce the volume of the degenerate subspace, not as a unique inversion.

4.3 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 model as a strongly informative prior rather than as a rigid template. A workable hierarchy, from fewest to most free parameters:

H0. Template cleaning with a 3D-informed mask. Use Planck 353 GHz (Q,U)(Q,U) as the dust template, scale with a global or piecewise βd\beta_d, but exclude or deweight pixels that the 3D+H I selection flags as LOS-decorrelated (Pelgrims target sample). This is the minimal use of the three ingredients. Tassis & Pavlidou proposed the starlight-polarization version of this mask in 2015. It is a hypothesis to be tested in rr-coverage, not already an improvement.

H1. Layered MBB with amplitudes fitted at a reference frequency. Fix {ψi}\{\psi_i\} from the orientation likelihood. Put Gaussian priors on {βd,i,Td,i}\{\beta_{d,i},T_{d,i}\} from the grain+3D model, with a width that reflects grain-scenario uncertainty rather than a delta function. Fit {AiP}\{A_i^P\} (and optionally a small δβd,i\delta\beta_{d,i}) to the ensemble of frequency maps, with Nlay=1N_{\mathrm{lay}}=1 on most high-latitude pixels and 22–33 only where the 3D map and H I demand it. The CMB is a blackbody added in (Q,U)(Q,U). Degeneracies are controlled because ψi\psi_i is not free, but they are not eliminated.

H2. Layered MBB plus intra-layer spin moments. Allow complex spin-2 moments ω1β,i\omega_1^{\beta,i} and ω1T,i\omega_1^{T,i} per layer, or, cheaper, a single set of spin moments on the sum of layers to capture beam averaging and unresolved substructure (Chluba et al. 2017; Mangilli et al. 2021; Vacher et al. 2022, 2023a). At LiteBIRD sensitivity the Vacher et al. (2022) result implies that moments in both βd\beta_d and TdT_d are required. The 3D model predicts the expected size and the real/imaginary structure of those moments, which can be imposed as a Gaussian prior, preventing the moments from absorbing CMB. If the layer angles are taken as fixed and correct, the remaining complex moments of the sum should be small; that is a testable claim, not an assumption.

H3. Harmonic-space hybrid. For rr likelihoods of the BICEP/Keck type, replace the single-SED, single-Rℓ\mathcal{R}_\ell dust model by the cross-spectra of the layered forward model, plus a low-order spin-moment expansion of the residual. The layered model supplies the frequency dependence of RℓBB\mathcal{R}_\ell^{BB} and of CℓBB(νi×νj)/CℓEE(νi×νj)C_\ell^{BB}(\nu_i\times\nu_j)/C_\ell^{EE}(\nu_i\times\nu_j); the moments absorb isotropic SED distortions from averaging. This is the natural interface to LiteBIRD and CMB-S4 power-spectrum pipelines, and to constrained ILC that deprojects the predicted moment maps. Simulation-based inference trained on 3D-informed simulations is a consumer of the same product, not a substitute for the physics.

Synchrotron is treated in parallel: a power law with spatially varying βs\beta_s, optionally a curvature term, and an independent ψs\psi_s (the synchrotron-emitting volume is not the dust-emitting volume). H I-based ψi\psi_i should not be imposed on synchrotron. Neutral-gas complexity can inform synchrotron depolarization, but not its orientation. Joint dust–synchrotron modelling is needed below ∼100 GHz\sim 100\,\mathrm{GHz} because a dust residual with a wrong βd\beta_d is degenerate with synchrotron, and conversely.

4.4 Why the combination is more than the sum of parts

Each ingredient removes a specific degeneracy, incompletely:

Degeneracy in a 2D MBB fitConstrained byResidual freedom
One βd\beta_d vs. a sum of βd,i\beta_{d,i}3D layersLayer SED still latent; conversion δAV→τ353\delta A_V\to\tau_{353}
Intensity βd\beta_d vs. polarized βd\beta_dgrain-size / alignment scenarioComposition, optical constants; astrodust vs two-population
Frequency-independent ψ\psi vs. rotating ψ(ν)\psi(\nu)per-layer BPOS\mathbf{B}_{\mathrm{POS}} if correctly matched to layersVelocity ≠\neq distance; starlight sparsity
Geometric pp vs. grain p0p_0B\mathbf{B} tangling (H I coherence) vs. RAT falignf_{\mathrm{align}}γ\gamma unmeasured; p0cos⁡2γp_0\cos^2\gamma remains one number
ω1β\omega_1^\beta, Im(ω1)\mathrm{Im}(\omega_1) vs. CMB3D-predicted spin-moment priorUnresolved sub-beam turbulence
Which pixels need extra SED freedom3D cloud count + misalignment angleMasking-induced E/BE/B leakage

A 3D mass map without B\mathbf{B} overfits Q,UQ,U. A B\mathbf{B} map without 3D mass cannot assign SED differences to layers. Grain physics without either map can predict a mean p(ν)p(\nu) but not its spatial pattern. Together they specify the Tassis–Pavlidou two-cloud (and NN-cloud) model on the real sky as a prior. Liu et al. (2025) found that at high complexity the bias on rr for Simons Observatory is driven by spatial correlations between the dust index and the dust amplitude at higher orders, not by the amplitude structure itself. Those correlations are exactly what a layered model with per-layer physics encodes, and what a random-field moment model treats as nuisance variance.

4.5 Validation, not only subtraction

A physical forward model is falsifiable, as is a phenomenological MBB: both can be tested with held-out frequencies, residual spectra, and posterior-predictive checks. A high-dimensional physical model can be harder to falsify if flexible amplitudes and moments absorb discrepancies. The tests that matter are held-out and end-to-end.

Required tests:

  • Fit only ancillary data plus 353 GHz, freeze the model, and predict 217/220 and then 95/150 GHz Q,UQ,U, Δψ\Delta\psi, EEEE, BBBB, and all cross-frequency Rℓ\mathcal{R}_\ell. Perform this separately on Pelgrims target/control samples and on the BICEP/Keck and SPIDER fields, where the current answer is “no sky-average decorrelation”. The model must predict that too, not only the Pelgrims excess.
  • 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.
  • On the target sample, the model must predict the sign and the approximate magnitude of Δψ(353,217)\Delta\psi(353,217) from the H I-measured α\alpha and the 3D-inferred ΔTd\Delta T_d, before looking at 217 GHz. Pelgrims et al. estimated a 15% polarized-intensity-ratio variation; a grain+3D model should predict a number of that order from TdT_d and βd\beta_d differences, not fit it.
  • The predicted p353p_{353} from H I coherence plus RAT p0cos⁡2γp_0\cos^2\gamma should match Planck p353p_{353} with residuals that correlate with tracers of grain alignment (e.g. AVA_V, radiation field), not with the CMB. Halal et al. (2024) is the existing version of this test for 3D complexity versus p353p_{353}.
  • Cross-spectra RℓBB(νi,νj)\mathcal{R}_\ell^{BB}(\nu_i,\nu_j) of the model, evaluated on the actual experiment bandpasses, should be compared to Planck 217×\times353 and to BICEP/Keck 95×\times150×\times220×\timesPlanck 353 without introducing a free Rℓ\mathcal{R}_\ell. A free Rℓ\mathcal{R}_\ell is reserved for the residual. Residual covariances should be compared to Guillet et al. (2026).
  • Spectral-imaging tests with bolometric interferometry (Regnier et al. 2024, A&A 686, A271—not a Maffei et al. paper) detect the rr systematic from unaccounted decorrelation in a CMB-S4-like setup. They are a detector, not a corrector, and pair naturally with a predictive model: a successful 3D model should make the detection vanish.
  • Ablation tests: 3D extinction, H I orientation, starlight tomography, and grain priors added separately. A strong 353 GHz calibration may dominate the result and conceal that the purported physical inputs add little predictive information.
  • Blinded end-to-end rr simulations in which the reconstruction receives mock Gaia stars, H I spectra, and starlight-polarization catalogues—not the true 3D cube. Include MHD substructure, a Bubble-like cavity, multiple grain compositions, out-of-model dust, synchrotron, lensing, realistic bandpasses, beams, calibration/leakage systematics, and survey masks. Compare against one-MBB fits, independent Q/UQ/U SEDs, scalar and spin moments, standard and constrained ILC, and phenomenological decorrelation likelihoods. Report Δr\Delta r, σr\sigma_r, residual DℓBBD_\ell^{BB}, coverage, sky fraction, and failure rate—not only goodness of fit.

Until those exist, the construction is a well-posed physical model for a demonstrated mechanism, not a turn-key replacement for moment expansions and multi-frequency residual fits.

4.6 An rr-budget, as a scaling, not as a result

After any linear or parametric clean, the residual dust B-mode power at a CMB band is of order the dust power there times the fraction of that power the model does not describe. The residual is positive-definite, so it biases rr upward rather than adding noise. Frequency-dependent rotation leaks dust EE into BB.

Order of magnitude in a BICEP/Keck-like deep field, using a modified blackbody with βd≃1.55\beta_d\simeq 1.55 at 19.6 K19.6\,\mathrm{K} and a tensor BBBB spectrum of order 0.07 μK20.07\,\mu\mathrm{K}^2 at ℓ≃80\ell\simeq 80 for r=1r=1: dust at 150 GHz is worth ∼0.1\sim 0.1 in rr-equivalent, and at 95 GHz ∼0.02\sim 0.02. Each percent of unmodelled dust power at 150 GHz is therefore worth about 10−310^{-3} in rr in the cleanest deep field. LiteBIRD’s target of order 10−310^{-3} and a CMB-S4 statistical error of a few ×10−4\times 10^{-4} require the unmodelled fraction to stay well below a percent. That is the scale of a degree of coherent unmodelled rotation, or of an unmodelled effective index dispersion of several hundredths—numbers that Planck XI does not rule out on small patches, that Osumi et al. (2021) leave open at high latitude, and that Pelgrims et al. (2021) detect statistically in a selected few-percent of high-latitude pixels.

If a fraction ff of the polarized power in a deep field comes from Pelgrims-like sightlines, the bias is ff times the two-cloud residual. For ff between 0.050.05 and 0.30.3 that is a few 10−410^{-4} to about 10−310^{-3} at 150 GHz—if those sightlines have comparable polarized amplitudes and if the cleaning does not already absorb the effect into a free βd\beta_d or Rℓ\mathcal{R}_\ell. This is exactly the next-generation error budget, and it is a bias. It is not a measurement of the residual of any existing 3D pipeline.

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


5. Limitations, missing data, and residual systematics

A layered physical model reduces, but does not eliminate, foreground bias on rr. The leading limitations, several of them blockers for a claim of controlled BBBB residuals, are as follows.

The core inverse problem remains underdetermined. The proposed inputs do not identify layer-resolved polarized amplitudes or SEDs. Discretizing the LOS is not the same as making the inversion identifiable. At one frequency, more than two layer amplitudes are underdetermined; similar MBB SEDs preserve the degeneracy across bands; nearly aligned angles make even two amplitudes ill-conditioned.

No end-to-end demonstration connects the ancillary products to controlled BBBB residuals or unbiased rr. Masking, spin-2 residuals, synchrotron covariance, bandpasses, beams, instrumental systematics, lensing, and ancillary-map uncertainty intervene before an rr likelihood is reached. Present “high-complexity” PySM suites (layered MBBs, d12) capture LOS decorrelation phenomenologically but do not ingest Gaia 3D maps or RAT-based alignment as observed products.

The H I-velocity-to-extinction-distance association is assumed more often than it is inferred. Velocity is not distance at high latitude. Galactic rotation does not provide a unique kinematic distance; two clouds at the same velocity can sit at different distances, and a single cloud can be split across velocities by turbulence. A wrong match assigns an angle to the wrong SED. Joint reconstructions of H I and 3D dust exist at Galactic-disk scales; they are not a solved high-latitude layer catalogue.

3D maps trace extinction, not submillimetre polarized emission. The conversion δAV→τ353\delta A_V\to\tau_{353} assumes a dust-to-extinction ratio and a size distribution that may vary with environment (RVR_V variations, grain growth in molecular gas, carbon depletion). High-latitude sightlines are mostly diffuse, where Lenz et al. (2017) find ∼10%\sim 10\% scatter in NHI/EB−VN_{\mathrm{HI}}/E_{B-V}, but 10% layer-to-layer opacity errors, if correlated with TdT_d or ψ\psi, are enough to leave a residual at the r∼10−3r\sim 10^{-3} level. Zelko et al. (2023) already find a factor-of-two variation in the emission-to-reddening conversion. CIB anisotropies further contaminate intensity-based τ353\tau_{353} and TdT_d maps at high latitude.

Distance completeness. Gaia-based maps reach ∼1.25 kpc\sim 1.25\,\mathrm{kpc} at high resolution (Edenhofer et al. 2024). That covers most high-latitude dust, but IVCs/HVCs can contribute polarized emission at high ∣b∣|b| and lie outside the local volume or are poorly sampled by stars. Pelgrims et al. used IVC/LVC H I components precisely because they are common at high latitude. A pipeline that trusts only Gaia distances will mis-model those layers unless H I velocity components are included as additional, distance-free clouds.

The inclination γ\gamma is poorly measured. Polarized emission constrains BPOS\mathbf{B}_{\mathrm{POS}}, not BLOS\mathbf{B}_{\mathrm{LOS}}. The factor cos⁡2γ\cos^2\gamma is degenerate with p0p_0 and with grain alignment efficiency. H I filament detectability encodes some γ\gamma information (filaments end-on are less linear on the sky), and starlight p/AVp/A_V encodes some, but neither is a precision γ\gamma map. Wrong γi\gamma_i mis-weight the layers in Q,UQ,U and therefore bias the predicted ψ(ν)\psi(\nu). The tangential-field assumption for the Local Bubble shell is not confirmed by polarization-fraction statistics (Halal et al. 2024).

Grain model uncertainty is a systematic on p(ν)p(\nu) and βdP−βdI\beta_d^{\mathrm{P}}-\beta_d^{\mathrm{I}}. Silicate+carbon with unaligned carbons versus astrodust versus models with magnetic dipole emission make qualitatively different predictions for the polarized SED below ∼100 GHz\sim 100\,\mathrm{GHz}. Until FIR polarimetry distinguishes these models, a B-mode pipeline must marginalize over at least two grain scenarios.

Alignment physics is not a solved 3D field. RAT predictions for aalign(U,nH,γrad)a_{\mathrm{align}}(U,n_{\mathrm{H}},\gamma_{\mathrm{rad}}) are tested statistically, not cloud-by-cloud at high latitude. Dense, shielded gas can have inefficient alignment; locally enhanced radiation fields can raise TdT_d and aaligna_{\mathrm{align}} together. Mis-modelling alignment produces a spatially variable p0p_0 that the pipeline may try to absorb into γ\gamma or into extra layers.

Unresolved sub-beam structure still generates spin moments. Even a perfect cloud catalogue does not resolve MHD turbulence inside a CMB beam (LiteBIRD’s 140 GHz beam is ∼31′\sim 31'; CMB-S4 will be finer but not arbitrarily fine). Intra-beam B\mathbf{B} tangling and TdT_d variations remain, and they are exactly what the Chluba/Vacher moment expansion is for. The 3D model reduces the moments; it does not set them to zero. Vacher et al. (2022) showed that at LiteBIRD noise, omitting TdT_d moments biases rr even when βd\beta_d moments are included. That conclusion survives a 3D prior.

Planck 353 GHz polarization S/N is limited at high latitude. All 3D polarized models are calibrated on Planck 353 GHz Q,UQ,U. Extrapolating a smoothed 353 GHz template cannot recover small-scale ψ(ν)\psi(\nu). Higher-resolution polarized dust data are needed to calibrate layers at CMB-S4 resolution. Circular calibration and validation on the same 353 GHz maps must be avoided: held-out 217 GHz is the minimum.

The nonlocal nature of E/BE/B. Pixel masking and spatially varying weights create EE-to-BB mixing; angle and amplitude errors contribute differently to BBBB; reproducing an H I–Planck map correlation at large scales does not establish control of the degree-scale dust BB modes used for rr.

Instrumental and analysis systematics mimic LOS decorrelation. Bandpass mismatch, I→PI\to P leakage, gain mismatch, and residual CMB in the 353 GHz template all produce frequency-dependent (Q,U)(Q,U) residuals. Pelgrims et al. found that correcting Planck residual systematics increased the significance of LOS decorrelation, which is encouraging (the astrophysical effect is not an artefact), but the converse is also possible in other data sets. 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 βs\beta_s and possible curvature; a possible polarized AME at the percent level; CO and other line polarization at a poorly known level; polarized compact sources. A perfect dust model does not yield an unbiased rr if synchrotron is wrong. The 3D dust/B-field construction does not automatically help synchrotron, which requires low-frequency surveys (C-BASS, S-PASS, QUIJOTE, Planck LFI, future 5–20 GHz polarimetry).

Statistical cost. Extra layers and moments consume degrees of freedom. Vacher et al. (2022) already saw σr\sigma_r increase because dust moments correlate with the CMB SED. A 3D model that adds unconstrained amplitudes per layer will degrade σr\sigma_r more than it reduces bias. The whole point of the three ingredients is to replace free maps by constrained ones: ψi\psi_i from H I/stars, τi\tau_i from Gaia, βd,i\beta_{d,i} and Td,iT_{d,i} from grain physics with tight priors. If those 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 B-mode pixel.
  • A joint 3D reconstruction that uses Gaia extinction, H I cubes, and starlight polarization in one likelihood, rather than matching products ex post.
  • FIR polarimetry from ∼100\sim 100 to 850 μm850\,\mu\mathrm{m} to decide between grain models and to measure TdT_d in polarization, where CIB contamination is weaker (Hensley & Draine).
  • Direct 3D polarized emission maps, which do not exist; extinction and Zelko-like 3D temperature are proxies for the intensity side only.
  • 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)n(a) and RAT alignment, observe it with LiteBIRD/CMB-S4 bandpasses and with realistic ancillary-data incompleteness, and show that the pipeline in §4 recovers rr without a bias larger than the statistical error.

Until those exist, the honest status is: the three ingredients specify a well-posed physical model for the exact mechanism (Tassis & Pavlidou 2015; Pelgrims et al. 2021) that makes a single-MBB-per-pixel subtraction fail; they can already be used to flag and down-weight the worst sightlines and to predict the frequency dependence of Rℓ\mathcal{R}_\ell and ψ(ν)\psi(\nu) at a level that phenomenological decorrelation parameters cannot, provided the layer association is validated on held-out frequencies; they are not yet a turn-key replacement for spin-moment expansions and multi-frequency residual fits. The right operational posture for LiteBIRD and CMB-S4 is to implement the layered forward model as an informative prior and as a diagnostic, to keep a low-order complex βd\beta_d–TdT_d spin-moment expansion for unresolved averaging, and to treat grain-composition and γ\gamma uncertainties as discrete systematics that are marginalized, not assumed away.


6. Summary

Spectral-index decorrelation of polarized dust is the statement that the effective MBB parameters and, when SED and orientation covary, the polarization angle of Galactic dust are frequency-dependent, because Stokes I,Q,UI,Q,U are integrals of emitters that do not share a common (Td,βd,ψ)(T_d,\beta_d,\psi). A single modified blackbody per pixel—or a 353 GHz template scaled by a global βd\beta_d—is then a misspecified model. Recovered CMB Q,UQ,U are strongly biased only when SED and angle vary together (McBride, Bull & Hensley 2023); that is precisely the regime the three ingredients target. Sky-average Planck and BICEP/Keck analyses do not detect a loss of correlation at the level relevant to r≃0.01r\simeq 0.01; a selected multi-cloud, misaligned sample does (Pelgrims et al. 2021); next-generation r∼10−3r\sim 10^{-3} experiments need control below the percent residual that current CMB data cannot guarantee.

Three-dimensional dust maps (presently Gaia extinction cubes, Zelko-like 3D temperatures, and 3D RVR_V, matched probabilistically to H I clouds, with the Local Bubble as a candidate first layer) replace the integral by a sum over a few distance layers and identify the sightlines where the Tassis–Pavlidou mechanism can operate. The grain size distribution and RAT alignment determine which grains emit in polarization, why βdP\beta_d^{\mathrm{P}} need not equal βdI\beta_d^{\mathrm{I}} and p=p(ν)p=p(\nu), and how to convert extinction into a polarized SED given the local radiation field and density—under an adopted grain scenario that must be marginalized. The sky-projected magnetic field, from H I filaments, starlight tomography, and Planck 353 GHz, assigns ψi\psi_i (and a geometric proxy for depolarization, not γ\gamma) to each layer, so that frequency-dependent ψtot\psi_{\mathrm{tot}} is predicted rather than fitted, if velocity is correctly matched to distance.

The combined pipeline is a layered, physically parameterized, band-integrated forward model whose remaining freedom is a small set of amplitudes and intra-layer spin moments, fitted to the CMB frequency maps under priors taken from the ancillary 3D data. Its limitations are real and some are blockers: extinction is not emission, velocity is not distance, γ\gamma and grain composition are underconstrained, the inverse problem remains underdetermined, beam averaging still requires complex moments, and no end-to-end rr demonstration yet exists. 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 spectral-index decorrelation names.


What I took from the critic

The controlling physics is SED–angle covariance, not “any multi-layer beam”. I rewrote the opening mechanism, the two-cloud discussion, and the McBride citation around that point, including the comparable-amplitude condition and Poh & Dodelson’s result that temperature-only extrapolation noise is small compared with then-current errors.

I adopted the critic’s convention demands: IAU versus COSMO UU sign, emission EVPA versus BPOS\mathbf{B}_{\mathrm{POS}} (the 90∘90^\circ rotation), and the statement that BPOS\mathbf{B}_{\mathrm{POS}} does not measure γ\gamma and does not make the polarized sky a derived quantity.

I corrected the moment expansion (amplitude convention; spin-2 complex moments rather than scalar H2 moments), the ILC description, the Rℓ\mathcal{R}_\ell–Δβ\Delta\beta formula’s assumptions, and the τi\tau_i conversion (the factor 1.0861.086 and the non-universal κ353/κV\kappa_{353}/\kappa_V). I removed the unusable RAT formula and the unattributed “15% is a 10% SED-shape change” sentence.

I restated the observational status as the critic demanded: Planck 2018 XI and Sheehy & Slosar walked back Planck Int. L; BK18 and BK XVI find no decorrelation in the deep field; SPIDER finds spatial SED variation without a full-region LOS detection; Pelgrims is a selected-sample detection with a stated pp-value range.

I accepted the identifiability critique as a blocker, not a quibble: the three products do not supply the joint distribution; layer association must be probabilistic; H0 masking is not already an improvement; a physical model is not uniquely falsifiable relative to an MBB. The pipeline section now states what is data, what is prior, and what is fitted, and it lists held-out prediction, ablation, and end-to-end rr simulations as requirements rather than as accomplished facts.

What I took from the optimist

I kept the draft’s useful skeleton (spatial versus LOS decorrelation, two-cloud algebra, alignment as a size high-pass filter, the degeneracy table, the H0–H3 hierarchy, held-out Pelgrims tests) and upgraded it with the optimist’s strongest verified observational handles.

The Local Bubble as a candidate dominant first layer (Skalidis & Pelgrims 2019; O’Neill et al. 2025), with Halal, Clark & Tahani (2024) as the necessary caution that the tangential-field assumption is not seen in p353p_{353} and that complexity beyond the shell matters. Zelko et al. (2023) as the existing prototype of a 3D emission/temperature map. Zhang & Green (2025) as a 3D grain-size proxy via RVR_V. Pelgrims et al. (2024) as the first degree-scale starlight tomography. BICEP/Keck XVI as evidence that H I templates carry some BB-mode phase information down to 95 GHz, with the 10–65% correlation quoted honestly. Planck Int. XVII as per-component SED evidence, including the extra ∼100 GHz\sim 100\,\mathrm{GHz} component.

I adopted spin moments as the interface between the forward model and the likelihood (Vacher et al. 2023a,b, 2025), including the angle-rotation relation, and Guillet et al. (2026) as a residual-covariance diagnostic, with the SRoll2-versus-PR4 split preserved. Liu et al. (2025) and Osumi et al. (2021) supply the argument that high-latitude βd\beta_d variance is poorly bounded and that β\beta–amplitude correlations, not just amplitude structure, drive rr bias. I corrected the spectral-imaging citation to Regnier et al. (2024) and tempered the βdP−βdI\beta_d^{\mathrm{P}}-\beta_d^{\mathrm{I}} offset to Planck XI’s 0.05±0.030.05\pm 0.03. I included an rr-budget as a scaling, and the FIR/high-frequency roadmap as the data that decide grain composition versus layer calibration.

What I rejected, and why

I rejected the original draft’s claim that multiple layers automatically make pp, ψ\psi, and β\beta frequency-dependent, and the claim that with BPOS\mathbf{B}_{\mathrm{POS}} “the polarized sky becomes a derived quantity”. Both are physically false: rotation needs SED–angle covariance and comparable amplitudes, and a projected orientation is not an amplitude or an inclination.

I rejected presenting extreme two-cloud rotations of “tens of degrees” as generic, and I rejected converting Pelgrims’s 15% band-ratio spread into a 10% SED-shape change. I rejected the printed RAT formula, the mixed moment kernels, the statement that standard ILC nulls a single dust SED, and the assertion that a physical model is falsifiable while an MBB is not.

I rejected several of the optimist’s upgrades that I could not verify as published results, rather than folding them in on trust: a 2026 velocity-resolved 3D H I map attributed to McCallum, a 2026 O’Neill catalogue of ∼1700\sim 1700 dust clouds with IVCs in distance, and a set of 2026 alignment, Faraday, and LiteBIRD-index papers that did not resolve as citable. The scientific need those citations were meant to fill (velocity–distance matching, IVCs, alignment tests) is already stated as a limitation.

I rejected the optimist’s implication that the high-latitude sum is already a finite, catalogued, solved object whose phases are predicted rather than fitted. The Local Bubble is a strong first-layer prior; Halal et al. (2024) show that dust beyond the shell still sets p353p_{353}. I rejected treating the rr-budget table as a measured residual of a 3D pipeline: it is a scaling that shows why sub-percent control matters, not a demonstration that the control exists. I rejected temperature-dominance as a settled sky property: Guillet et al. (2026) find it for SRoll2 residual covariances and not for PR4.

I rejected “H0 is already an improvement” and “3D extinction provides the same information in true distance as an H I cloud count”. Masking must be scored on rr-coverage; extinction and H I are complementary, not equivalent. I also rejected the idea that the layered model replaces moments. Unresolved sub-beam structure still needs spin moments; the 3D model is a prior on those moments, not a substitute.