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

Dedusting the CMB with three‑dimensional dust maps, grain physics and the plane‑of‑sky magnetic field

Summary of the argument

Polarized Galactic dust in every beam is a superposition of physically distinct emitting regions. CMB component separation assumes one modified blackbody and one polarization angle per pixel. When the regions differ both in their spectral energy distribution and in their magnetic‑field orientation, their relative complex weights change with frequency. The polarized SED is then not a modified blackbody, the polarization angle rotates with frequency, and maps at two frequencies decorrelate. The effect is usually called "spectral‑index decorrelation" because spatial variation of the dust index $\beta_d$ is its simplest cause. A more accurate name is frequency decorrelation of polarized dust, since temperature differences, polarization‑efficiency differences and, above all, angle differences between regions produce it too.

The three ingredients in the question supply exactly the structure that a per‑pixel fit lacks. Three‑dimensional dust maps give the number, distances, columns and environments of the regions. Grain physics gives each region's temperature distribution, which grains are aligned, the offset between polarized and total‑intensity SEDs, and the conversion between optical and submillimetre polarization. The sky‑projected field gives the complex phase of each region. Together they turn an arbitrary frequency‑dependent complex field into a low‑dimensional forward model whose parameters are constrained by data other than the CMB maps themselves. They act as priors and templates, not as a map to subtract. The chain is: layer structure, then polarized spectral weights and phases, then frequency‑dependent complex dust modes, then a constrained foreground likelihood, then a verified bound on the bias in $r$. The case is strongest at high Galactic latitude for deep‑field measurements of the recombination bump and weakest for full‑sky measurements of the reionization bump.

1. Polarized dust emission and the B‑mode budget

Aspherical grains spin with their short axis on average along the local magnetic field, through radiative‑torque alignment (Lazarian and Hoang 2007; Andersson, Lazarian and Vaillancourt 2015). Their far‑infrared to millimetre emission is linearly polarized perpendicular to the plane‑of‑sky field $\mathbf B_\perp$. Following Lee and Draine 1985 and Planck intermediate XX, the Stokes parameters of a region with source function $S_\nu=\sigma_\nu n_{\rm H}B_\nu(T_d)$, opacity $\sigma_\nu\propto\nu^{\beta_d}$, field inclination $\gamma$ to the plane of the sky and plane‑of‑sky angle $\chi$ are

For one region with uniform properties, $\mathcal P_\nu=p_0F\cos^2\gamma\,S_\nu e^{2i\psi}$, where $F\le1$ absorbs imperfect alignment and unresolved field dispersion within the region. Total intensity carries its own inclination factor, so $I$ and $\mathcal P$ do not share a common geometric prefactor. If the SED inside the region is uniform, $F$ is frequency independent and drops out of everything that follows. Only the frequency‑dependent part of intra‑region structure matters, and it is deferred to section 5.

Planck established the reference numbers for the diffuse sky.

Primordial tensors produce a recombination bump at $\ell\approx80$ and a reionization bump at $\ell<20$. The last row is the same magnitude as the excess BICEP2 attributed in 2014 to $r\approx0.2$, so the dust amplitude in that field expressed in $r$ units, call it $r_d$, is of order 0.1 to 0.2 at 150 GHz. The cleanest patches identified by Planck intermediate XXX are a few times lower. On the 70 per cent of sky LiteBIRD uses for the reionization bump, $r_d$ is one to two orders of magnitude higher.

A fraction $f_{\rm res}$ of dust $BB$ power left in a cleaned map near $\ell\approx80$ biases $r$ by roughly $f_{\rm res}\,r_d$ before marginalization. One per cent of residual dust is therefore $\delta r\approx10^{-3}$ in the BICEP/Keck field, several $\times10^{-4}$ in the cleanest patches, and far beyond budget on LiteBIRD's large sky fraction. The actual bias depends on how the residual projects onto the CMB SED in the joint fit, on the mask, on delensing and on the likelihood, so this scaling sets the required level of control, not the final number. That level is per‑cent in power for deep fields and per‑mille for the full sky.

2. The mechanism of frequency decorrelation

2.1 The single modified blackbody per pixel

Parametric methods (Commander, FGBuster, the BICEP/Keck multi‑component likelihood), template scaling, and moment‑based or constrained ILC methods all start from

with the same real scalar $\hat g$ multiplying $Q$ and $U$. This fixes the polarization angle to be frequency independent and the polarized SED to be a pure modified blackbody, often with $\beta^P_d=\beta^I_d$. Band‑pass integration and conversion to thermodynamic units are implied.

2.2 Superposition and its consequences

The sky in any beam is a sum over emitting regions $i$: distinct clouds along the line of sight, sub‑beam structures, and, for a power spectrum, all pixels in the mask:

The necessary condition for the dangerous case has two parts. If the regions have different SEDs but the same angle, the sum has spectral curvature but no rotation, and a scalar SED with one or two extra parameters absorbs most of it. If the angles differ but the SEDs are identical, the relative weights do not change with frequency and the summed pattern stays perfectly correlated across bands. Decorrelation of the kind no scalar SED can represent needs both frequency‑dependent relative weights and misaligned orientations (Tassis and Pavlidou 2015). Given both, four things break.

• The effective SED is not a modified blackbody. $|\mathcal P_\nu|$ is a vector sum, so the effective index between two bands depends on the pair and can lie outside the range of the $\beta_i$ when layers partially cancel. Across the sky this appears as SED curvature, the object of the moment expansion of Chluba, Hill and Abitbol 2017.
• The angle rotates with frequency. $\psi_{\rm eff}(\nu)=\tfrac12\arg\sum_i\mathcal P_i\hat g_i(\nu)$ changes with $\nu$. This is a spin‑2 effect, which is why Vacher et al. 2023 generalized moments to complex spin‑moments.
• The polarized SED differs from the intensity SED even with identical grains, because intensities add as scalars and polarizations as vectors.
• Maps at two frequencies decorrelate. With $R_\ell(\nu_1,\nu_2)=C_\ell^{\nu_1\times\nu_2}/\sqrt{C_\ell^{\nu_1\times\nu_1}C_\ell^{\nu_2\times\nu_2}}$, Gaussian scatter $\sigma_\beta$ of a single‑layer index gives $R=\exp[-\tfrac12\sigma_\beta^2\ln^2(\nu_1/\nu_2)]$, and two uncorrelated layers whose relative amplitude changes by a factor $\rho$ between the bands, with power ratio $f$, give $R=(1+\rho f)/\sqrt{(1+f)(1+\rho^2f)}$.

A two‑layer example makes the numbers concrete. Layer 2 carries 30 per cent of the 353 GHz polarized amplitude, both layers have the same intrinsic $p$ and $\beta^P=\beta^I$, and the bands are 150 and 353 GHz.

The same configuration gives $R(217\times353)\approx0.999$, which is below Planck's sensitivity, yet $1-R^2\approx0.7$ per cent at $150\times353$. $R$ is a second‑order statistic, so a non‑detection at the per‑cent level in $R$ does not exclude a residual at the level that matters for $r\sim10^{-3}$.

A rotation $\delta\psi$ of the dust pattern between the template frequency and the target frequency converts dust $E$ into $B$. For a rotation coherent over a patch, $\delta C_\ell^{BB}\approx4\,\delta\psi^2C_\ell^{EE,\rm dust}\approx8\,\delta\psi^2C_\ell^{BB,\rm dust}$ for $E/B\approx2$. For the spatially varying rotation that superposition produces, the mixing kernel depends on the coherence of the misalignment, and the same expression with the rms $\langle\delta\psi^2\rangle$ gives the order of magnitude. Rotations of 1°, 2° and 3° leave 0.2, 1 and 2 per cent of the dust $BB$ power. The first row of the table is therefore a per‑cent effect.

2.3 What survives scalar SED freedom

The one‑template residual $(1-R_\ell^2)\,C_\ell^{\rm dust}$ is the irreducible leftover when a single dust template per region is scaled to the target band, which is the situation of template cleaning and of a cross‑spectral likelihood with one dust component and a decorrelation parameter. Pipelines that float $\beta_d$ and $T_d$ per pixel absorb much of the curvature part. They cannot absorb the differential scaling of $Q$ and $U$, because a real scalar multiplies both. The surviving residual is then of order $4\langle\delta\psi^2\rangle C_\ell^{EE,\rm dust}$, plus higher‑order curvature. What finally biases $r$ is the projection of that residual onto the CMB SED in the joint fit, which can be larger or smaller than the raw residual and must be measured with end‑to‑end simulations. The rotation term is the part that only external structural information can remove.

2.4 Observational status

Planck intermediate results L in 2017 reported a departure of the $217\times353$ $BB$ correlation ratio from unity, marginal in each region but above 99 per cent confidence in combination, increasing toward low column density. Sheehy and Slosar 2018 and Planck 2018 XI found the evidence not robust once noise bias and systematics were treated fully. Planck 2018 XI concluded that, if the large‑region limit applies to small patches, "decorrelation might not be a problem for CMB experiments aiming at a primordial B‑mode detection limit on the tensor‑to‑scalar ratio $r\simeq0.01$". It made no statement about $r\sim10^{-3}$, and by the second‑order argument above it could not. BK18 marginalizes over a decorrelation parameter $\Delta_d$, defined at $217\times353$ and $\ell=80$ and scaled to other band pairs with $\ln^2(\nu_1/\nu_2)$ following Vansyngel et al. 2017; the data are consistent with $\Delta_d=1$. BICEP/Keck XVI found no decorrelation in the BICEP/Keck field between filamentary and non‑filamentary dust, or from dust associated with intermediate‑velocity HI.

The line‑of‑sight mechanism has nonetheless been seen directly where it should be. Pelgrims et al. 2021 used HI4PI cloud counts from Panopoulou and Lenz 2020 and HI‑fibre orientations from Clark and Hensley 2019 to select sightlines with several clouds and misaligned fields, rejected the no‑decorrelation null at high significance in the 217 and 353 GHz angles, found the effect stronger for more misaligned fields, and showed that an intrinsic 15 per cent variation between clouds in the 353‑to‑217 GHz polarized ratio suffices. Ritacco et al. 2023, with SRoll2 maps at $\ell=4$–32 on 80 to 97 per cent of the sky, found that first‑order modified‑blackbody moments fit $EE$ but not $BB$, and that residual maps at different frequencies are not fully correlated, so spatial variation of the polarization angle with frequency must be included. Vacher et al. 2023 showed that the same spin‑moments make the dust $E/B$ ratio and $EB$ correlation frequency dependent. Guillet et al. 2026 built complex residual maps around the mean SED and found the SRoll2 residual covariances compatible with temperature‑only variation, while Planck PR4 is compatible with neither temperature‑only nor index‑only variation. The three data sets are consistent with one picture: the effect is localized to sightlines with comparable, misaligned layers, is absent or small on large‑sky averages and in the BICEP/Keck field, and current Planck systematics are already at the level of the astrophysical signal.

2.5 Why more free parameters is not the answer

McBride, Bull and Hensley 2023 showed that single‑sightline SED fitting measures only the parameters of broad distributions, and that two‑parameter fits are biased at more than $10\sigma$ precisely when both the SED and the angle vary along the sightline. Vacher et al. 2022 found that free first‑order $\beta$ and $T$ moments recover an unbiased $r$ for LiteBIRD's 100–402 GHz bands only at $\sigma(r)=8.8\times10^{-4}$, because the moments are degenerate with $r$. Hensley and Bull 2018 and Sponseller and Kogut 2022 reached similar conclusions. Extra freedom buys robustness with sensitivity, and it still cannot represent a rotation unless $Q$ and $U$ are given separate SEDs, which doubles the degeneracy. External information that fixes the structure of the superposition is the alternative. Its value lies as much in reducing the variance penalty of flexible models as in reducing bias.

3. How each ingredient helps on its own

3.1 Three‑dimensional maps of thermal dust emission

What exists is mostly three‑dimensional extinction, not emission. Gaia parallaxes with photometry or XP spectra give the extinction density along each sightline: Bayestar19 (Green et al. 2019), Leike, Glatzle and Enßlin 2020, Lallement et al. 2019 and Vergely et al. 2022, and Edenhofer et al. 2024, which uses 54 million stars to map dust to 1.25 kpc at 14′ resolution in 516 distance bins. A 3D emission map adds a per‑voxel emission law,

Zelko, Finkbeiner, Lee and Green 2022 took this step, assigning $T$ and $\beta$ to 500 million Bayestar19 voxels over three quarters of the sky by matching five Planck and IRAS intensity bands. They resolved two components of distinct temperature along the Cepheus sightline and found the ratio of emission optical depth to reddening varying by a factor of two. Two caveats set the fidelity of this product. Five intensity maps contain far fewer independent data than the number of voxels, so per‑voxel SEDs are regularized interpolations, not measurements. And at high latitude the same bands contain a large cosmic‑infrared‑background fraction, so voxel temperatures there are only meaningful after CIB separation, for example by fitting the HI‑correlated part of the intensity. A complementary decomposition comes from HI velocity. The HI4PI Gaussian decomposition of Panopoulou and Lenz 2020 finds about three clouds per 0.2 deg² pixel in the north Galactic cap and 2.5 in the south, and the Milky Way Atlas of McCallum et al. 2026 now reconstructs a velocity‑resolved 3D HI map within 1.25 kpc by combining HI4PI with the Edenhofer map, separating local from distant emission. For the geometry that matters most, Skalidis and Pelgrims 2019, Pelgrims et al. 2020 and O'Neill et al. 2024 find that the Local Bubble wall at 200–300 pc likely dominates the polarized signal on most high‑latitude sightlines, and Halal, Clark and Tahani 2024 find that the observed $p$ decreases where the column is split among components at different distances within 1.25 kpc.

On its own this ingredient gives five things.

• Decorrelation risk maps. The number of layers and their column fractions $w_i=A_{V,i}/A_V$ say where one modified blackbody is adequate and where it is not. This is a mask and field‑selection tool available today, and in HI form it is what Pelgrims et al. 2021 used to find the effect.
• Layer SED shapes from intensity. Intensity has far higher signal‑to‑noise than polarization. A CIB‑controlled 3D emission fit gives $T_i$ and $\beta^I_i$ per layer, which fixes the shape of each layer's SED up to the intensity‑to‑polarization transfer of section 3.2.
• Layer templates. Per‑layer intensity maps can be fitted against $Q$ and $U$ with regional coefficients, as done with HI velocity layers by Ghosh et al. 2017 and Adak et al. 2020 and with HI templates in the BICEP/Keck field by BICEP/Keck XVI.
• Realistic simulations. Multi‑layer models, from Martínez‑Solaeche, Karakci and Delabrouille 2018 (PySM d12) through DustFilaments (Hervías‑Caimapo and Huffenberger 2022) to the spin‑moment models of Vacher et al. 2024, are how priors on $\Delta_d(\ell;\nu_1,\nu_2)$ and pipeline validations are built. Measured 3D structure replaces assumed layer statistics.
• Temperature structure. The distance of each layer from heating sources predicts where $T$ differs between layers, which Guillet et al. 2026 suggest may dominate the SED variance in polarization.

The map alone does not give emission per unit extinction, which layers polarize, or their angles.

3.2 The dust grain size distribution

Grain size does not fix the millimetre index. In the Rayleigh limit the opacity per unit mass of a homogeneous grain is size independent, so $\beta$ is set by the frequency dependence of the dielectric function, which depends on composition and structure: amorphous silicate versus carbon, the temperature‑dependent two‑level‑system behaviour of Meny et al. 2007, porosity and mantles, and magnetic‑dipole emission from iron inclusions (Draine and Hensley 2013). Size enters the polarized millimetre SED in four other ways. The equilibrium temperature scales as $T\propto a^{-1/(4+\beta)}$, so a factor of two in size shifts $T$ by about 2 K. Grains below 0.01–0.02 μm are stochastically heated and contribute little below 353 GHz. Only grains above an alignment threshold $a_{\rm align}$ polarize; in radiative‑torque theory $a_{\rm align}$ is 0.05–0.1 μm and depends on the radiation field $U$, gas density and temperature, and grain shape (Hoang and Lazarian 2016), and the empirical alignment function of Hensley and Draine 2023 has negligible alignment below 0.01 μm and complete alignment above 0.1 μm. And composition correlates with size in every model. Polarized and total emission therefore weight the size distribution differently:

The offsets $\beta^P-\beta^I$ and $T^P-T^I$ are predictions of a size distribution plus alignment function plus composition, and the models disagree. Draine and Fraisse 2009 found $p$ rising from about 3 per cent at 100 μm to 15 per cent at 1 mm if only silicates align, and flat at 6–10 per cent if both silicate and carbon grains align. Planck intermediate XXII measured $\beta^P=1.59\pm0.02$ against $\beta^I=1.51\pm0.01$ and a 21 ± 6 per cent decrease of $p$ from 353 to 70 GHz. Guillet et al. 2018 built models compatible with the Planck constraints on translucent sightlines. Hensley and Draine 2023 argue that a single composite "astrodust" above about 0.02 μm most naturally explains a frequency‑independent far‑infrared polarization fraction. The mean offset is unsettled, so it must stay a global nuisance parameter. What the size distribution supplies is its layer dependence and sign structure. The distribution also varies in space: Zhang and Green 2025 mapped $R(V)$ in three dimensions from 130 million Gaia XP spectra and found it varying with extinction, consistent with accretion at low extinction and coagulation at higher extinction. $R(V)$ is tied to the aligned‑grain size through the wavelength of maximum starlight polarization, $R_V\approx5.6\,\lambda_{\max}/\mu{\rm m}$ (Whittet and van Breda 1978), a relation with considerable scatter.

On its own this ingredient helps in five ways.

• The intensity‑to‑polarization transfer. It converts the high‑signal‑to‑noise $(\beta^I_i,T^I_i)$ of section 3.1 into a physically bounded $(\beta^P_i,T^P_i)$ per layer, with the layer‑to‑layer differences tied to environment rather than fitted freely.
• Which layers polarize. The mass above $a_{\rm align}(U_i,n_i)$ sets the relative $p_i$ of layers, so a dense, shielded layer with poor alignment can dominate $I$ but not $\mathcal P$. This is the physical origin of many $\beta^P\neq\beta^I$ sightlines, and it is a statement about weights, not about $\beta$.
• The optical‑to‑submillimetre bridge. Starlight polarization is produced by grains near 0.1–0.2 μm at $\lambda_{\max}$, while 353 GHz polarized emission is dominated by the largest aligned grains. The ratio $R_{P/p}$ and its dependence on $R(V)$ and temperature is a size‑distribution quantity. Without it, starlight tomography in section 3.3 cannot be turned into per‑layer emission amplitudes.
• Bounding the moments. A plausible range of size distributions and radiation fields bounds the variance of $(\beta,T)$ within a beam, giving priors for moment and spin‑moment expansions and a way to discriminate temperature from index variation, as Guillet et al. 2026 did.
• Band placement. If temperature variation dominates, high‑frequency leverage matters most; Fuskeland et al. 2023 found that extending LiteBIRD's upper cutoff from 400 to 500–600 GHz reduces $\sigma(r)$ by 30–50 per cent when $\beta$ and $T$ vary on the sky, with most of the gain by 500 GHz.

These quantities should enter as regional hyperparameters with population priors, since there is no per‑voxel $n(a)$, $U$ or gas density, and the alignment physics is still being revised even for optical polarization spectra (Tram et al. 2026).

3.3 The sky‑projected magnetic field

Each layer contributes $\mathcal P_i=p_0F_i\cos^2\!\gamma_i\,\epsilon_iA_{V,i}\,e^{2i(\chi_i+\pi/2)}$ times its SED. The plane‑of‑sky angle $\chi_i$ sets the phase of each term, so the misalignment between layers is what turns SED differences into rotation and vector cancellation. Planck's 353 GHz map gives only the emission‑weighted mean. Per‑layer angles come from three sources. Starlight polarization of stars at known Gaia distances is the direct one: Panopoulou et al. 2019 demonstrated tomography on a diffuse sightline; Pelgrims et al. 2023 published the Bayesian inversion BISP‑1, which returns the number of dust layers, their distances and their Stokes $q_i,u_i$; Pelgrims et al. 2024 produced the first degree‑scale tomography map over 4 deg² with RoboPol, finding that most sightlines intersect more than one cloud; Uppal et al. 2025 added the frequentist method TRIShUL, reliable for cloud polarizations above 0.1 per cent with a few hundred stars; and PASIPHAE (Tassis et al. 2018) is surveying the high‑latitude sky with the WALOP polarimeters for this purpose. HI fibres are the second source: Clark et al. 2015 showed they trace the dust polarization angle; Clark and Hensley 2019 built Stokes maps from HI4PI that correlate at $r>0.75$ with Planck 353 GHz $Q$ and $U$ at the map level without free parameters and reproduce the $E/B$ asymmetry at $\ell\lesssim120$; Clark 2018 showed that velocity coherence traces line‑of‑sight tangling; and BICEP/Keck XVI detected dust polarization correlated with the HI morphology template at the 10–65 per cent level in $B$ modes over $20<\ell<200$, with $\beta=1.54\pm0.13$. Local Bubble field models (Alves et al. 2018; Pelgrims et al. 2020; O'Neill et al. 2024) are the third, supplying the angle of the dominant high‑latitude layer over the whole sky. Faraday rotation measures and Galactic‑scale field models are not substitutes: rotation measures weight free electrons and the signed line‑of‑sight field, and large‑scale models describe the coherent disk and halo field, not the turbulent wall of a 200 pc bubble.

On its own this ingredient contributes the following.

• Misalignment maps. Knowing $\chi_i$ per layer identifies where the vector‑sum problem exists at all. Pelgrims et al. 2021 needed only this, plus cloud counts, to detect the effect.
• Per‑layer polarized amplitudes without knowing $\gamma$ or alignment. The same aligned grains produce optical dichroic extinction and submillimetre polarized emission, so a layer's starlight $q_i,u_i$ times $R_{P/p}$ gives its 353 GHz $Q_i,U_i$; $\cos^2\gamma_i$ and $F_i$ cancel. This is the correction Tassis and Pavlidou 2015 proposed. It needs per‑cloud calibration: Mehandiratta et al. 2026 measured $R_{P/p}=3.67\pm0.05$ MJy sr$^{-1}$ in a RoboPol field against the Planck global $5.42\pm0.05$, and found the ratio differing between the two dominant overlapping clouds. The width of that variation is the prior width on the layer amplitudes.
• The pivot for spin‑moments and the $E/B$ statistics. The mean field direction defines the pivot around which spin‑moments are expanded, and filament–field alignment explains the dust $E/B\approx2$ and positive $TE$ and $TB$ in the models of Huffenberger, Rotti and Collins 2020, Clark et al. 2021, Hervías‑Caimapo and Huffenberger 2022 and Konstantinou et al. 2022.
• Coverage where stars are sparse. HI‑based angles provide $\chi$ per velocity component over the full high‑latitude sky, and with the McCallum et al. 2026 map they can be placed in distance bins.

4. Combining the three into a dedusting strategy

4.1 The forward model

Per pixel and frequency, in specific‑intensity units,

with the layer quantities supplied externally and their 2026 status stated honestly.

The free parameters are the CMB $Q,U$, synchrotron, any polarized anomalous microwave emission, per‑layer amplitude corrections $\alpha_i$ with priors of a few tens of per cent from the $R_{P/p}$ scatter, angle corrections $\delta\psi_i$ with priors from the tomography errors, a global or regional offset $\Delta\beta^P$, and residual moment terms for what the layer sum does not capture: intra‑layer variance of $(\beta,T)$, temperature‑dependent opacity, magnetic‑dipole emission, and dust beyond the mapped volume. The model generates the discrete‑cloud part of the frequency‑dependent angle and curvature; it does not generate everything, and the residual terms must remain in the likelihood with their own priors.

4.2 Parameter counting and conditioning

The counting decides whether the external information helps. If the layer angles are supplied, each layer has one real amplitude, and at a single frequency the two equations for $Q$ and $U$ separate two layers with a condition number of order $1/\sin2\Delta\psi$. The frequency dependence then follows from the known SED shapes, and the CMB is separated by its distinct blackbody spectrum as usual. If the angles are free, each layer has two real amplitudes and the layers can be told apart only through their SED difference, with a condition number set by $\Delta\beta\ln(\nu_1/\nu_2)$ and $\Delta T$, which is the ill‑posed problem of section 2.5. The third ingredient therefore changes the conditioning of the whole problem, and its errors are the ones that matter most. The corrections $\alpha_i$, $\delta\psi_i$ and $\Delta\beta^P$ help only if their priors are tighter than the direction along which they are degenerate with $r$; where they are not, they should be fixed and the residual carried as a systematic. A hierarchical model, with per‑pixel amplitudes, regional grain and alignment hyperparameters, and population priors on the layer statistics, is the identifiable form. For LiteBIRD's 15 bands the layered fit with $N\le3$ and regional spectral parameters is feasible; for Simons Observatory or CMB‑S4 with 6–9 bands, $N\le2$ is realistic.

4.3 Three products

Layered templates. Per‑layer $Q_i,U_i$ maps with fixed SED shapes, fitted with regional amplitude coefficients, as BICEP/Keck XVI did with HI templates.

Predicted spin‑moment templates. Writing $\mathcal P_\nu=\hat g(\nu;\bar\beta,\bar T)\,[\mathcal W_0+\mathcal W_1^\beta\ln(\nu/\nu_0)+\mathcal W_1^T\,\partial_T\ln B_\nu|_{\bar T}+\dots]$ with $\mathcal W_0=\sum_i\mathcal P_i$ and $\mathcal W_1^\beta=\sum_i\mathcal P_i(\beta_i-\bar\beta)$, the first‑order moments are complex and their phase relative to $\mathcal W_0$ encodes the rotation. The forward model predicts these maps. Fitting a few global amplitudes of predicted moment templates, or deprojecting their SEDs in a constrained moment ILC (Remazeilles, Rotti and Chluba 2021), keeps the correct functional form while breaking the moment–$r$ degeneracy that free moments suffer. The delta‑map method of Ichiki et al. 2019 and the power‑spectrum moments of Mangilli et al. 2021 and Azzoni et al. 2021 admit the same treatment. A mistemplated deprojection removes the wrong mode and can increase the bias, so the templates must carry their posterior uncertainty and be tested against a phase error and a missing cloud.

Cross‑spectral prior. For BK‑style likelihoods, the model gives the multi‑frequency dust covariance $C_\ell^{ab}=\sum_{ij}\hat g_i(a)\hat g_j(b)\,C_\ell^{ij}$ from the layer maps, a region‑specific prediction for $\Delta_d(\ell;\nu_1,\nu_2)$ that replaces the assumed $\ell$ and $\ln^2(\nu_1/\nu_2)$ scalings. It should enter as a prior with an envelope, not as a hard constraint, until the predicted $217\times353$ correlation matches Planck 2018 XI on large sky, Pelgrims et al. 2021 on selected sightlines, and BICEP/Keck XVI in the deep field at once.

4.4 Validation, circularity and staged use

Each prediction is testable with existing data: the $\psi(217)-\psi(353)$ map, which Pelgrims et al. 2021 measured statistically; the residual covariances of Guillet et al. 2026; the frequency dependence of $E/B$, $TE$ and $TB$ of Vacher et al. 2023 and Ritacco et al. 2023; and the correlation of deep‑field $B$‑mode maps with per‑layer templates. Beyond these, the method has to pass the following before it is used on $r$.

• Null test on the BICEP/Keck field. The pipeline must recover $\Delta_d\approx1$ there and not inject a correction, with Pelgrims‑style target and control splits repeated inside the mask.
• Held‑out frequencies. Predict a band not used in building the model and compare.
• Per‑ingredient error propagation. Monte Carlo over $\sigma(\psi_i)$, $\sigma(w_i)$, $\sigma(\beta^I_i)$, $R_{P/p}$ scatter and CIB in $T_i$, reporting residual $BB$ at the target band and $\delta r$ for each, at the angle accuracy PASIPHAE actually delivers at high latitude.
• End‑to‑end bias. Inject layered 3D skies with realistic angle errors, missing layers, band‑pass mismatch, instrumental angle errors, synchrotron, AME and lensing; report $\langle\hat r\rangle-r_{\rm in}$ and $\sigma(r)$ against a single modified blackbody, free moments, and the 3D prior at 2026 data quality.
• Circularity. Planck 353 GHz enters the 3D emission fit and the HI‑template calibration and is also the test data. Use independent splits, frequencies or surveys for validation.
• Birefringence. A predicted frequency‑dependent dust angle is degenerate with an instrumental angle error and with isotropic cosmic birefringence; the dust model must be shown not to bias birefringence at the 0.1° level.

Deployment is staged. Now: risk maps, masks and $\Delta_d$ priors from cloud counts, HI angles and 3D extinction, plus HI‑layer templates in deep fields. As PASIPHAE, CIB‑controlled 3D emission and 3D $R(V)$ mature: per‑layer templates and priors with posterior samples rather than a single map. With LiteBIRD and CMB‑S4 data: joint hierarchical fits in which the structural priors are refined only along directions where they are tighter than the degeneracy. Where the external data are uninformative, the model degrades to a decorrelation prior or a mask, which is the graceful fallback.

4.5 What the payoff can be

The payoff is a range that the external data decide, not a forecast. If a 3D map shows one layer carrying 90 per cent of the polarized column and tomography shows the rest misaligned by less than 20°, the model predicts $|\delta\psi_{\rm eff}|\lesssim0.5°$ between 353 and 150 GHz even for $\Delta\beta=0.3$, and a residual below $10^{-4}$ in $r$; such sightlines can be cleaned with a single modified blackbody with confidence and full sensitivity. If instead the map shows the first row of the table in section 2.2, the model supplies a per‑cent‑level correction with an error bar instead of an inflated $\Delta_d$ envelope, or the sightline is masked. The credible goal is to reduce the decorrelating residual by a factor of several in the fields that matter and to demonstrate, by simulation and held‑out tests, that what remains is below the experiment's foreground allocation.

5. Limitations, missing data and residual systematics

Limits of the 3D maps, and the nuisance terms they imply. They measure extinction, not emission: the factor‑of‑two variation in $\tau_{353}/E(B-V)$ found by Zelko et al. 2022 is a 100 per cent amplitude effect that must be modelled through grain properties or absorbed in $\alpha_i$. Voxel SEDs are regularized interpolations of five intensity maps that contain CIB and zodiacal light and that constrain the total‑intensity SED, not the large‑grain SED that polarization samples; the test is whether high‑latitude voxel temperatures survive CIB separation. Resolution of 14′ corresponds to 0.8 pc at 200 pc and 4 pc at 1 kpc, so a voxel is not a single physical cloud. The Edenhofer map stops at 1.25 kpc, adequate at $|b|\gtrsim30°$ but not at the intermediate latitudes that LiteBIRD's reionization‑bump analysis needs. Stellar shot noise and stellar‑model systematics propagate into $w_i$; methods such as FLOWER (O'Callaghan, Mandel and Gilmore 2025) target the high‑latitude regime. HI components are velocity peaks, not slabs with one $T$, $\beta$ and $\psi$; the velocity‑to‑distance mapping is weak at high latitude; molecular and ionized gas are missed. Within a layer, sub‑beam variation of $(\beta,T)$ correlated with density and field tangling set a floor on high‑$\ell$ decorrelation that should be measured with MHD or DustFilaments skies degraded to the 3D‑map and tomography resolution.

Limits of the grain physics. $R(V)$ is a one‑parameter proxy for grains of 0.01–0.3 μm, while millimetre polarized emission is dominated by grains above 0.1 μm. Astrodust, THEMIS and silicate‑plus‑carbon mixtures predict different $p(\nu)$ and $\beta^P-\beta^I$, and Planck intermediate XXII, Planck 2018 XI and Ritacco et al. 2023 do not agree on the offset, so the mean offset stays a nuisance parameter that LiteBIRD's bands above 300 GHz can fit and Planck's cannot. Millimetre dielectric functions and their temperature dependence are poorly measured in the laboratory. Alignment thresholds depend on unmeasured local $U$ and density, and the empirical alignment function is constrained only on average.

Limits of the field data. Starlight tomography needs enough polarized stars behind each layer; at high latitude polarizations are a few tenths of a per cent, the resolution is of order a degree, and nothing beyond the last usable star is constrained. Only 4 deg² of degree‑scale tomography is published. HI‑based angles fail where the gas is molecular or ionized, and the 10–65 per cent $B$‑mode correlation in the BICEP/Keck field, not the map‑level $r>0.75$, is the relevant number for cleaning. The inclination $\gamma_i$ cancels in the starlight transfer but matters for how a layer's intensity maps to its polarized amplitude, and it is available only indirectly through Local Bubble models and polarization‑fraction statistics. A single $\psi_i$ per layer ignores the tangling that Planck 2018 XII's anticorrelation between $p$ and the angle dispersion shows to be a strong depolarizer.

Method‑level residuals. Every structural error maps onto a residual with a definite form: $\delta\alpha_i$ leaves part of layer $i$ at the wrong SED, $\delta\psi_i$ leaks about $8\,\delta\psi_i^2$ of that layer's $BB$ power, a missing layer reproduces the original problem in miniature, and an error in $\Delta\beta^P$ is a coherent bias. A wrong high‑confidence prior can bias $r$ more than a flexible phenomenological model, which is why posterior samples, model averaging over grain prescriptions, and explicit missing‑layer uncertainty are part of the method rather than refinements. Band‑pass mismatch and polarization‑efficiency errors mimic SED curvature, and an instrumental angle error is degenerate with a dust rotation and with birefringence. The SRoll2 versus PR4 disagreement in Guillet et al. 2026 shows Planck systematics already at the level of the astrophysical signal. Polarized anomalous microwave emission and the dust–synchrotron correlation reported by Planck 2018 XI sit inside LiteBIRD's 40–90 GHz channels and inside the same likelihood; a 3D dust model does not protect $r$ from them.

The two science cases are different. For Simons Observatory, CMB‑S4 and BICEP/Keck on high‑latitude deep fields near $\ell\approx80$, a Local Bubble wall plus HI fibres plus targeted tomography may be enough, and per‑cent control is the requirement. For LiteBIRD's reionization bump on 70 per cent of the sky, the 1.25 kpc ceiling, degree‑scale tomography, and synchrotron and AME dominate, and per‑mille control is required. The same residual fraction does not convert to the same $\delta r$ in the two cases.

What is missing. An all‑sky, per‑layer plane‑of‑sky field from starlight tomography, which PASIPHAE is building; CIB‑controlled 3D emission maps with $\beta$ and $T$ per layer at CMB‑relevant resolution and beyond 1.25 kpc; a measured, per‑cloud link between 3D $R(V)$, submillimetre emissivity and $R_{P/p}$, of which Mehandiratta et al. 2026 is the first step; a settled polarized‑versus‑intensity SED offset, which needs LiteBIRD‑class sensitivity above 300 GHz; the line‑of‑sight field per layer; laboratory millimetre opacities of realistic grain analogues; and, above all, an end‑to‑end bias study at real data quality. None of these gaps prevents the first stage of section 4.4. Each sets the floor on the residual that the later stages can reach.

References

• Pelgrims et al. 2021, A&A 647, A16, arXiv:2101.09291
• Ritacco et al. 2023, A&A 670, A163, arXiv:2206.07671
• Vacher et al. 2022, A&A 660, A111, arXiv:2111.07742; Vacher et al. 2023, $E/B$ and $EB$ frequency dependence, arXiv:2210.14768; Vacher et al. 2024, arXiv:2411.11649
• Guillet et al. 2026, A&A 705, A177, arXiv:2510.18305
• Mehandiratta et al. 2026, cloud‑to‑cloud variation of $R_{P/p}$, arXiv:2510.15389
• BICEP/Keck XIII (BK18), arXiv:2110.00483; BICEP/Keck XVI, arXiv:2210.05684; Tristram et al. 2022, arXiv:2112.07961
• LiteBIRD Collaboration 2023, PTEP 2023, 042F01, journal; Fuskeland et al. 2023, arXiv:2302.05228
• Planck intermediate XX, arXiv:1405.0872; Planck intermediate XXII, arXiv:1405.0874; Planck intermediate XXX, arXiv:1409.5738; Planck intermediate L, arXiv:1606.07335; Planck 2018 XI, arXiv:1801.04945; Planck 2018 XII, A&A
• Tassis & Pavlidou 2015, arXiv:1410.8136; McBride, Bull & Hensley 2023, arXiv:2207.14213; Hensley & Bull 2018, arXiv:1709.07897; Sponseller & Kogut 2022, ApJ; Vansyngel et al. 2017, arXiv:1611.02577
• Edenhofer et al. 2024, arXiv:2308.01295; Zelko et al. 2022, arXiv:2211.07667; Zhang & Green 2025, Science 387, 1209, journal; McCallum et al. 2026, Milky Way Atlas, arXiv:2607.07451; O'Callaghan, Mandel & Gilmore 2025, MNRAS
• Skalidis & Pelgrims 2019, arXiv:1908.08706; Pelgrims et al. 2020, A&A; O'Neill et al. 2024, arXiv:2410.17341; Halal, Clark & Tahani 2024, arXiv:2404.11009
• Hensley & Draine 2023, arXiv:2208.12365; Guillet et al. 2018, arXiv:1710.04598; Draine & Fraisse 2009, arXiv:0809.2094; Tram et al. 2026, GRADE‑POL II, arXiv:2606.13123
• Pelgrims et al. 2023, BISP‑1, arXiv:2208.02278; Pelgrims et al. 2024, arXiv:2404.10821; Uppal et al. 2025, TRIShUL, arXiv:2510.25911
• Clark & Hensley 2019, arXiv:1909.11673; Clark 2018, arXiv:1802.00011; Clark et al. 2021, ApJ; Panopoulou & Lenz 2020, arXiv:2004.00647
• Hervías‑Caimapo & Huffenberger 2022, arXiv:2107.08317; Konstantinou et al. 2022, arXiv:2204.13127

Cited from memory of the literature without links: Chluba, Hill & Abitbol 2017; Mangilli et al. 2021; Azzoni et al. 2021; Remazeilles, Rotti & Chluba 2021; Ichiki et al. 2019; Sheehy & Slosar 2018; Ghosh et al. 2017; Adak et al. 2020; Martínez‑Solaeche et al. 2018; Huffenberger, Rotti & Collins 2020; Clark et al. 2015; Panopoulou et al. 2019; Tassis et al. 2018 (PASIPHAE); Alves et al. 2018; Lee & Draine 1985; Lazarian & Hoang 2007; Hoang & Lazarian 2016; Andersson, Lazarian & Vaillancourt 2015; Meny et al. 2007; Draine & Hensley 2013; Whittet & van Breda 1978; Green et al. 2019; Leike, Glatzle & Enßlin 2020; Lallement et al. 2019; Vergely et al. 2022.

What I took from the critic

• Factual corrections, all verified against the sources. Whittet and van Breda's coefficient is 5.6, not 5.5. LiteBIRD's bands are 34–448 GHz, not 40–402. Vacher et al. 2022 used only the 100–402 GHz bands. Planck 2018 XI's statement is about $r\simeq0.01$ and says nothing about the next generation; my "pointedly excludes" was editorial and is gone. Planck intermediate L is now described by its confidence level and column‑density trend rather than a single correlation number. The "PASIPHAE's term" attribution is removed.
• The Stokes model. Total intensity carries its own inclination factor, so I now use the Lee and Draine form and state that $I$ and $\mathcal P$ do not share a common geometric prefactor.
• Faraday rotation measures and Galactic‑field models do not give per‑layer inclination or the angle of a bubble wall. Both are removed as sources and named as non‑substitutes.
• $\beta$ is not a grain‑size number. Section 3.2 is rewritten: in the Rayleigh limit $\beta$ comes from composition and structure, and size acts through temperature, alignment threshold, polarized‑versus‑total weighting and the optical bridge. The "prior of width 0.1" claim is gone; the mean offset is a nuisance parameter and the grain model supplies its layer dependence.
• Parameter counting. With supplied angles the layers have $N$ real amplitudes, not $2N$. The overclaim that "everything is generated by the model" is replaced by an explicit list of what is not generated and stays as residual moment terms.
• The residual arithmetic. $1-R^2$ is now labelled as the one‑template bound, with a separate statement of what survives per‑pixel scalar SED freedom and a note that the bias on $r$ is the projection of the residual in the joint fit. The $E\to B$ formula is presented for a coherent rotation, with the rms version as an order of magnitude.
• The $\delta r$ conversion is now written as $f_{\rm res}\,r_d$ with $r_d$ sourced from Planck intermediate XXX, and the deep‑field and full‑sky science cases are separated throughout.
• Weak links made explicit. CIB and regularization in the 3D emission fit, extinction versus emission, HI components not being slabs, velocity not being distance, per‑cloud rather than global $R_{P/p}$, intra‑layer tangling, polarized AME and the dust–synchrotron correlation, mistemplated deprojection, circularity, and the birefringence degeneracy are all now in the pipeline or limitations with a corresponding test or nuisance parameter. The BICEP/Keck XVI null and the 10–65 per cent $B$‑mode correlation are used instead of the map‑level $r>0.75$.
• The referee's tests. The end‑to‑end bias study, the BICEP/Keck null test, per‑ingredient error propagation, condition numbers, held‑out frequencies and the requirement that a cross‑spectral prior match all three data sets before replacing $\Delta_d$ are now the validation section.

What I took from the optimist

• Naming and the necessary condition. The problem is frequency decorrelation of polarized dust, of which index variation is one cause, and line‑of‑sight decorrelation needs both frequency‑dependent relative weights and misaligned orientations. Both statements now open section 2.
• Ingredients as priors, not maps. The summary and section 4 present the three inputs as complementary priors on layer structure, spectral weights and phases, with the value split between bias reduction and the variance penalty of flexible models.
• Extinction versus emission is distinguished throughout, and the recast of grain size as a constraint on temperature, alignment, weighting and the optical bridge follows the optimist's formulation.
• Hierarchical modelling and posterior propagation. Regional hyperparameters with population priors, posterior samples instead of a single deterministic map, model averaging over grain prescriptions, explicit missing‑layer uncertainty, and the warning that a wrong high‑confidence prior can bias $r$ more than a flexible model are all adopted.
• Circularity safeguards, held‑out frequency tests and graceful fallback are in section 4.4.
• Near‑term deliverables first. Risk maps, masks, layer templates and decorrelation priors lead the staged plan; the all‑sky deterministic subtraction is not the goal.
• The causal chain summarizing the argument in one sentence is the last line of the summary.
• Each limitation tied to a test or nuisance term is the organizing principle of section 5.

What I rejected, and why

• The critic's LiteBIRD requirement of $\delta r<0.002$. The PTEP 2023 mission paper states a total uncertainty $\delta r<0.001$ including foregrounds, systematics and margin. I kept $10^{-3}$.
• The critic's claim that layers cannot be separated at a single frequency. That is true only with free angles. With supplied, distinct angles, the $Q$ and $U$ equations separate two layers at each frequency with a condition number of order $1/\sin2\Delta\psi$, independent of $\Delta\beta$. This is the whole point of the third ingredient, so I made the distinction explicit rather than conceding the claim.
• The critic's implicit standard that the essay is invalid without an end‑to‑end bias demonstration. The question asks how the ingredients can help, and the answer is a research programme with a staged, testable path. I adopted the demanded tests as the validation section and stated plainly that no such study exists yet, rather than dropping the pipeline.
• The critic's reading that grain physics adds nothing beyond the prior already used. A global width on $\beta^P-\beta^I$ is indeed not new. What the size distribution and alignment physics add is the layer dependence: which layers polarize, the temperature distribution of the emitters, and the optical‑to‑submillimetre transfer. I reframed rather than removed that ingredient.
• The critic's suggestion that the only near‑term deliverable is masking. Layer‑correlated templates and cross‑spectral priors with envelopes are also usable now, as BICEP/Keck XVI has shown for HI templates, so they stay in stage one.
• The optimist's advice to trim the numerical claims substantially. The numbers are what make the argument checkable, and several of the critic's corrections were only possible because the draft stated them. I trimmed the digressions and the reference list instead, and I verified every number that either reviewer questioned.
• The optimist's mention of [C II] as a handle on dark gas. It is outside the three named ingredients and has no CMB‑relevant map, so I left it out and noted only that HI misses molecular and ionized gas.