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

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

"Dedusting" is the PASIPHAE collaboration's term for removing polarized Galactic dust from CMB polarization maps well enough to measure or bound the tensor-to-scalar ratio at $r\sim10^{-3}$. Current B-mode analyses model dust in each pixel, or each patch, as one modified blackbody with one spectral index and one polarization angle. "Decorrelation of the spectral index" names what happens when that model fails: dust maps at two frequencies stop being scaled copies of each other. In the regime that matters for the next decade, the dominant cause is that the polarized dust signal in a beam is a vector sum over a few physically distinct magnetized dust layers whose spectra differ. Polarization is a spin-2 quantity, so the sum has a spectrum that is not a modified blackbody and an angle that depends on frequency. A scalar spectral index, however it is varied over the sky, cannot represent the second effect.

The three ingredients named in the question address the three pieces of that superposition that CMB frequency maps cannot fix on their own. Three-dimensional maps of thermal dust emission say how much dust sits at which distance along each sightline and how warm each piece is. The grain size distribution, together with composition and alignment physics, says how each piece's intensity turns into polarized intensity and how its polarized spectrum differs from its intensity spectrum. The sky-projected magnetic field says how the pieces' polarization vectors are oriented relative to one another, which is what turns spectral differences into angle rotation and cancellation. Section 1 sets the B-mode budget. Section 2 defines decorrelation and separates the three distinct effects that share the name. Section 3 treats each ingredient on its own. Section 4 combines them into a forward model, derives how accurate the external information must be, and gives three ways to insert the model into an $r$ likelihood. Section 5 lists what is missing.

1. Polarized dust emission and the B-mode budget

Interstellar grains are aspherical and spin with their short axes aligned, on average, with the local magnetic field. Radiative torques are the accepted alignment mechanism, developed by Lazarian and Hoang and reviewed by Andersson, Lazarian and Vaillancourt in 2015. The far-infrared to millimetre emission of aligned grains is linearly polarized perpendicular to the plane-of-sky field $\mathbf B_\perp$. For one optically thin emitting region with optical depth $\tau_{\nu_0}$ at a reference frequency, dust temperature $T_d$, opacity index $\beta_d$, plane-of-sky field angle $\chi$ and inclination $\gamma$ of $\mathbf B$ to the plane of the sky, the standard leading-order model is

Here $p_0$ is the intrinsic polarization of the aligned grain population and $F\le 1$ absorbs imperfect alignment and unresolved field tangling. Two things are approximate. Grain orientation also modifies $I_\nu$ at order $p_0$, and $F$ is frequency independent only if the same grains dominate at every frequency. Both corrections are small compared with the effects treated below, but the equation is a leading-order model, not an identity.

Planck fixed the phenomenology. Planck 2018 XII found a maximum polarization fraction of 22 % at 80′ resolution, an anticorrelation between $p$ and the polarization-angle dispersion $\mathcal S$ such that $\mathcal S\times p$ varies by only about 25 % while $p$ varies by a factor of 3 to 4, and polarization emission-to-extinction ratios that depend only weakly on column density. Planck 2018 XI measured a mean polarized spectral index $\beta_d^P=1.53\pm0.02$ at $T_d=19.6$ K, found the polarized SED well fitted by a single modified blackbody from 353 GHz to below 70 GHz, and confirmed dust $BB$ power falling roughly as $\ell^{-2.4}$, a $BB/EE$ ratio near 0.5, positive $TE$ and evidence for positive $TB$.

The B-mode connection is quantitative. Primordial tensor modes produce a recombination bump at $\ell\approx80$ and a reionization bump at $\ell<20$. Even in the cleanest few per cent of the sky, dust $B$-mode power at 150 GHz and $\ell\approx80$ is comparable to a primordial signal with $r\approx0.05$–0.1. The excess that BICEP2 reported in 2014 was of this size, and the joint analysis with Planck showed it to be consistent with dust plus lensing with no significant primordial component. The current limits and the targets are:

The unit of account for everything below is therefore

where $f_{\rm res}$ is the fraction of the dust $BB$ power at 150 GHz left in the cleaned map. With $r_{\rm dust}\approx0.1$, a one per cent residual is $\delta r\approx10^{-3}$, which is the entire LiteBIRD budget. Poh and Dodelson made this point in 2017 from Planck and 3D dust statistics: line-of-sight "extrapolation noise" from clouds at different temperatures becomes important once experiments are sensitive to $r<0.0015$.

2. What "spectral-index decorrelation" means

2.1 The one-modified-blackbody-per-pixel assumption

Parametric component separation, as in Commander, FGBuster and the BICEP/Keck multi-component likelihood, and the template extrapolation used in Planck and BICEP analyses, all rest on

with the same $(\beta_d,T_d)$ for $Q$ and $U$. This implies a frequency-independent angle $\psi$ and a polarized SED that is a pure modified blackbody. Often $\beta_d^P=\beta_d^I$ is assumed as well. Bandpass integration and unit conversion are implied throughout.

2.2 Three different effects share the name

Three distinct things make dust maps at two frequencies imperfectly correlated, and they need separate treatment because they have different statistics and different remedies.

Line-of-sight superposition. Within one beam, dust layers at different distances have different $(\beta,T)$ and different field angles. This produces a non-modified-blackbody SED, a frequency-dependent polarization angle, a polarized SED different from the intensity SED, and a correlation ratio below unity. It is the only one of the three that rotates the angle with frequency within a pixel. The three external ingredients target this effect.

Sub-beam mixing. Turbulent structure within a layer produces the same vector arithmetic at scales the beam does not resolve. It contributes to depolarization and to decorrelation at high multipoles. External maps at 14′ to degree resolution do not resolve it, so it must be treated statistically.

Spatial variation of the SED across the analysis region. From pixel to pixel the effective $\beta$ and $T$ vary, so no single scaling factor maps the 353 GHz map onto the 150 GHz map over a patch. This produces a correlation ratio below unity and SED curvature in power spectra, but no angle rotation within a pixel. The widely used formula for Gaussian scatter $\sigma_\beta$ in a spatially varying index with fixed morphology, uniform temperature, and index fluctuations independent of amplitude is

Amplitude–index correlations, masking, beam smoothing, noise and line-of-sight cancellation all alter it. Grain physics bounds $\sigma_\beta$ for this effect, but 3D maps and field angles are only indirectly relevant to it.

The rest of this section concerns the first effect, because it is the one that free spectral parameters cannot absorb and the one that the three ingredients are built to constrain.

2.3 Line-of-sight superposition and its four consequences

In any beam, the polarized sky is a sum over emitting layers $i$ along the line of sight:

where each $\mathcal P_i$ carries its own phase $2\psi_i$ set by that layer's field. Four things then break.

The effective SED is not a modified blackbody. $|\mathcal P_\nu|$ is a vector sum, so the effective index between any two bands depends on the pair, and it can lie outside the range spanned by the $\beta_i$ when layers partially cancel. Across the sky this appears as SED curvature and is the object of the moment expansion of Chluba, Hill and Abitbol 2017 and its polarized generalization to spin-2 fields by Vacher, Chluba and collaborators in 2023.

The angle rotates with frequency. $\psi_{\rm eff}(\nu)=\tfrac12\arg\sum_i\mathcal P_i\hat g_i(\nu)$ changes with $\nu$ whenever the $\psi_i$ differ and the $\hat g_i$ differ. Tassis and Pavlidou 2015 isolated the two necessary conditions with a two-cloud model: the ratio of the clouds' polarized intensities must change between the frequencies, and their fields must be misaligned. No scalar SED, however curved, represents a rotation. In the spin-moment language, writing $\mathcal P_\nu=\hat g(\nu;\bar\beta,\bar T)\,[\mathcal W_0+\mathcal W_1\ln(\nu/\nu_0)+\dots]$ with $\mathcal W_0=\sum_i\mathcal P_i$ and $\mathcal W_1=\sum_i\mathcal P_i(\beta_i-\bar\beta)$, the first-order rotation is

The rotation is the phase of the complex first moment relative to the zeroth moment. It vanishes when all layers share one angle, whatever their SEDs.

The polarized SED differs from the intensity SED even with identical grains. Intensities add as scalars, polarizations as vectors, so $\beta^P_{\rm eff}\neq\beta^I_{\rm eff}$ arises from geometry before any grain physics enters.

Maps at two frequencies decorrelate. For two layers that are spatially uncorrelated as random fields, with power ratio $f$ at the first frequency and a relative amplitude change $\rho$ between the two frequencies, $\mathcal R=(1+\rho f)/\sqrt{(1+f)(1+\rho^2 f)}$. After scaling a 353 GHz template to 150 GHz with the best single scalar, the residual is $1-\mathcal R^2\approx2(1-\mathcal R)$ of the 150 GHz dust power.

A worked example fixes the sizes. Layer 2 carries 30 % of the summed 353 GHz polarized amplitude, so its amplitude is 0.3 against 0.7 for layer 1 in common units. Both layers have the same intrinsic $p$ and the same $\beta^P=\beta^I$ within each layer. The arithmetic is pixel-level Stokes addition between 353 and 150 GHz.

The first-order spin-moment formula gives $-2.7°$ for the first row against the exact $-2.4°$. A rotation $\delta\psi$ coherent over scales larger than the multipoles of interest converts $E$ into $B$ with $\delta C_\ell^{BB}\approx4\,\delta\psi^2 C_\ell^{EE,\rm dust}\approx8\,\delta\psi^2 C_\ell^{BB,\rm dust}$ for $E/B\approx2$, so rotations of 1°, 2° and 3° leave about 0.2 %, 1 % and 2 % of the dust $BB$ power unmodelled. When $\delta\psi$ varies over the sky as the layer ratio varies, the leakage couples multipoles and splits between $E$ and $B$, and the estimate holds only for the rms rotation and to within a factor of about two. A direct calculation of the residual after the best global scalar for the first row gives 0.8 % of the total polarized power, consistent with the rotation estimate at that level. Either way, the first row is a $\delta r\sim10^{-3}$ effect in a field with $r_{\rm dust}\approx0.1$.

The correlation-ratio arithmetic is the reason Planck cannot see this. With layer 2's rms amplitude 3/7 of layer 1's, so $f\approx0.18$, and $\Delta\beta=0.3$, the two-layer formula gives $\mathcal R(217\times353)\approx0.999$, well inside Planck's uncertainty, while $\mathcal R(150\times353)\approx0.996$ and $1-\mathcal R^2\approx0.7$ % of the dust power at 150 GHz. $\mathcal R$ is a second-order statistic, and the residual it hides scales with the square of the frequency lever arm. Pelgrims et al. 2021 found that an intrinsic variation of only about 15 % in the 353-to-217 GHz polarized intensity ratio between clouds suffices to reproduce the angle effect they detected, so $\Delta\beta=0.3$ is an illustrative contrast, not a requirement.

2.4 Observational status

Planck intermediate results L in 2017 reported $\mathcal R^{BB}_\ell(217\times353)\approx0.95$ at $\ell\approx50$–160 in low-column regions. Sheehy and Slosar 2018 and Planck 2018 XI found this not robust once noise bias and systematics were handled, and Planck 2018 XI concluded that it found no evidence for decorrelation and that decorrelation might not be a problem for experiments aiming at $r\approx0.01$. That conclusion says nothing about $r\approx10^{-3}$. BK18 marginalizes a decorrelation parameter $\Delta_d$ defined as the 217×353 correlation ratio at $\ell=80$, scaled to other band pairs with $\ln^2(\nu_1/\nu_2)$ following Vansyngel et al. 2017 and with an assumed flat or linear $\ell$ dependence, and finds the data consistent with no decorrelation while allowing a few per cent.

The line-of-sight mechanism has been seen directly by selecting the sightlines where it should occur. Pelgrims et al. 2021 used HI4PI cloud counts from Panopoulou and Lenz 2020 and HI-fibre orientations from Clark and Hensley 2019 to pick high-latitude sightlines with several clouds and misaligned fields, and found 217–353 GHz polarization-angle differences that grow with the misalignment. Read together, the Planck null result in broad averages and the Pelgrims detection in a pre-selected sample say the effect is real, spatially sparse, and below Planck's power-spectrum noise unless one has independent three-dimensional information to select the affected pixels. That is the first job the 3D ingredients do.

Three further results constrain the form of the effect. Ritacco et al. 2023, using SRoll2 maps at $\ell=4$–32 over 80 to 97 % of the sky, found residuals at 100, 143 and 217 GHz after scaling 353 GHz with the mean SED, showed that first-order modified-blackbody moments fit the $EE$ residual spectra but not the $BB$ ones, and concluded that frequency dependence of polarization angles must be included. Vacher, Aumont and collaborators showed in 2023 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 polarized SED and found that SRoll2 residual covariances are compatible with a model containing temperature variations only, while Planck PR4 residuals are compatible with neither temperature-only nor index-only variation. The two data releases disagree, which is a warning that Planck systematics are already at the level of the astrophysical effect and that no conclusion about which parameter dominates is yet safe.

2.5 Why more free parameters is not the answer

One could add layers or moments as free parameters per pixel. McBride, Bull and Hensley 2023 showed with Fisher analysis that SED fitting on a single line of sight constrains only the parameters of broad distributions of $\beta$ and $T$, and that two-parameter fits are biased at more than 10σ in illustrative cases precisely when both the SED and the polarization angle vary along the sightline. Vacher et al. 2022 found that free $\beta$ and $T$ moments recover an unbiased $r$ for LiteBIRD only at $\sigma(r)\approx8.8\times10^{-4}$, because the moments are partly degenerate with the CMB. Hensley and Bull 2018 and Sponseller and Kogut 2022 reached similar conclusions: extra spectral freedom costs sensitivity and does not break the degeneracy.

The missing information is structural and joint. One needs to know which layers exist, how their polarized amplitudes are weighted, what SED each has, and how their Stokes vectors are oriented, with those properties attached to the same layer, because the residual is controlled by the correlation between a layer's SED and its angle, not by their marginal distributions. CMB frequency maps measure only the sum. The three ingredients constrain that structure from outside the CMB bands. Today they constrain it rather than fix it, and Section 4.2 quantifies how well they must do so.

3. How each ingredient helps on its own

3.1 Three-dimensional maps of thermal dust emission

A 3D emission map is a distance-resolved emissivity, built from a 3D extinction density and a per-voxel emission law:

The extinction density is now a Gaia-era data product: Bayestar19 by Green et al. 2019, Leike, Glatzle and Enßlin 2020, Lallement et al. 2019 and Vergely et al. 2022, and Edenhofer et al. 2024, who used 54 million stars to map dust from 69 pc to 1.25 kpc at up to 14′ with parsec-scale distance resolution, sampled at 516 logarithmically spaced bins, with 12 posterior samples released so that the correlated uncertainties can be propagated. The emission step is younger. Zelko, Finkbeiner, Lee and Green 2022 assigned $T$ and $\beta$ to 500 million Bayestar19 voxels and fitted the integrated emission to Planck and IRAS bands at 27′ to 110′ resolution. It is a heavily regularized global inverse problem. It resolved two dust components of different temperature along the Cepheus sightline at the 1σ level and found the ratio of emission optical depth to reddening varying by a factor of two. It is a demonstration of feasibility, not yet a validated per-layer SED product at the precision B-mode work needs. Two more things belong in the emission step rather than in a caveat: the Planck and IRAS intensities contain the cosmic infrared background and zodiacal light, which must be removed or modelled before a 3D emissivity is a dust quantity, and the intensity SED is dominated by all grains while polarization is dominated by the large aligned ones, so the intensity-to-polarization transfer of Section 3.2 is needed before intensity $(T_i,\beta_i)$ can be used for polarization.

A complementary "third dimension" comes from HI velocity. The HI4PI Gaussian decomposition of Panopoulou and Lenz 2020 finds on average three clouds per 0.2 deg² pixel in the north Galactic cap and 2.5 in the south, with intermediate-velocity clouds prominent in the north. McCallum et al. 2026 now tie HI velocity components to distance by combining HI4PI with the Edenhofer map in an information-field-theory reconstruction within 1.25 kpc.

The 3D maps also fix the default high-latitude geometry. The Local Bubble wall reconstructed by O'Neill, Zucker, Goodman and Edenhofer 2024 from the Edenhofer map has a mean distance of about 170 pc and a typical thickness of 35 pc, with a northern chimney. Skalidis and Pelgrims 2019, Pelgrims et al. 2020 and the 3D field model of O'Neill, Goodman, Soler, Zucker and Han 2024 find that the wall dominates the polarized signal on the majority of high-latitude sightlines. Halal, Clark and Tahani 2024 found no imprint of the wall's inclination on the 353 GHz polarization fraction, but did find that $p$ decreases when the column is substantially split among components at different distances within 1.25 kpc, an effect that vanishes if only dust within 270 pc is counted. The degree-scale starlight tomography of Pelgrims et al. 2024 at $b\approx22°$ finds several clouds out to about 2 kpc with field orientation changing with distance on most sightlines. The picture for recombination-bump fields at high latitude is therefore one dominant nearby layer plus one or two subdominant, more distant layers. The subdominant layer is the decorrelation threat, and 3D maps earn their keep by measuring its weight, not by reconstructing a continuous emissivity cube at CMB resolution.

On its own, the 3D emission map does five things.

• Decorrelation risk maps. The number of layers and their column fractions $w_i=A_{V,i}/A_V$ tell you where a single modified blackbody is adequate and where it is not. This is a masking and field-selection tool available today, and it is what Pelgrims et al. 2021 used in HI form to find the effect.
• Layer SEDs from intensity. Intensity is measured at far higher signal-to-noise than polarization. Where the 3D emission fit is trustworthy it gives $T_i$ and $\beta^I_i$ per layer, which fixes the shape of each layer's polarized SED up to the transfer in Section 3.2 and reduces the polarization problem to amplitudes and phases. The precision available today is that of a regularized reconstruction, and the temperature can be predicted rather than fitted only if the interstellar radiation field is modelled in 3D from the stellar distribution and the extinction itself, including shielding. That model is a task, not yet a product.
• Layer templates. Per-layer intensity or HI maps can be correlated with $Q$ and $U$. Ghosh et al. 2017 and Adak et al. 2020 built physical models of the 353 GHz polarized sky on HI velocity layers with a few field parameters per region and reproduced Planck's polarization statistics; BICEP/Keck XVI cross-correlated an HI morphology template with deep-field maps. Neither is yet a demonstration of distance-resolved polarized component separation.
• Simulations with the real layer configuration. Multi-layer sky models, from Martínez-Solaeche, Karakci and Delabrouille 2018, which became PySM's d12, through DustFilaments by Hervías-Caimapo and Huffenberger 2022, to the spin-moment skies of Vacher et al. 2024, use layer statistics calibrated to projected data. With 3D maps one can build simulations of a specific target field with its measured layer configuration rather than a random draw, which is how a region-specific $\Delta_d(\ell;\nu_1,\nu_2)$ prior should be made.
• Temperature structure. Distance from heating sources and the Local Bubble geometry predict where $T$ differs between layers. If the SRoll2 result of Guillet et al. 2026 holds up, temperature is the SED parameter whose variation matters most, and the 3D map is the only external handle on it.

The map alone does not give polarized emissivity, which layers polarize, their angles, or their inclinations.

3.2 The dust grain size distribution

Modified-blackbody parameters are set by grain physics, and the size distribution enters in a specific way. In the Rayleigh limit that holds at millimetre wavelengths, the far-infrared opacity per unit mass is independent of grain radius, so $\beta$ is a property of the material's dielectric function, of composition, porosity, two-level-system physics and iron inclusions, not of size. Size enters through two other factors. The equilibrium temperature of a grain of radius $a$ that absorbs geometrically in the optical and emits in the Rayleigh regime scales as $T\propto a^{-1/(4+\beta)}$, so a factor of two in size shifts $T$ by about 2 K at 20 K, and grains smaller than 0.01–0.02 μm are stochastically heated and contribute little below 353 GHz. And 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$, the gas density and temperature and the grain shape, as in Hoang and Lazarian 2016; the empirical alignment function of Hensley and Draine 2023 rises steeply with size and saturates for the largest grains. The polarized and total emission therefore weight the size distribution differently:

Thus $\beta^P-\beta^I$ and $T^P-T^I$ are predictions of a size distribution, an alignment function and a composition. Draine and Fraisse 2009 showed that if only silicates align, $p$ rises from about 3 % at 100 μm to about 15 % at 1 mm, while if silicate and carbonaceous grains both align it stays between 6 and 10 % at wavelengths longer than 100 μm. Planck intermediate XXII reported $\beta^P=1.59\pm0.02$ against $\beta^I=1.51\pm0.01$ between 100 and 353 GHz, a 3.6σ difference, and a $21\pm6$ % decrease of $p$ from 353 to 70 GHz. Planck 2018 XI, with better data, prefers a single modified blackbody for the polarized SED. Planck intermediate XXI measured the ratio of 353 GHz polarized emission to optical starlight polarization, $R_{P/p}=5.4\pm0.2\pm0.3$ MJy sr$^{-1}$, on selected diffuse sightlines with nearly orthogonal optical and submillimetre angles, and found it about 2.5 times higher than the dust models of the time predicted. Hensley and Draine 2023 argue that a single composite "astrodust" material for grains larger than about 0.02 μm most naturally explains both the lack of frequency dependence of the far-infrared polarization fraction and that optical-to-far-infrared polarization ratio. Guillet et al. 2018 built silicate-plus-carbon models compatible with the Planck constraints in translucent sightlines. Temperature-dependent opacities in the two-level-system model of Meny et al. 2007 and magnetic-dipole emission from iron nanoparticles in Draine and Hensley 2013 can further change the low-frequency slope.

The size distribution also varies in three dimensions. Zhang and Green 2025 mapped $R(V)$ 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.5\,\lambda_{\max}/\mu{\rm m}$ in the relation of Whittet and van Breda. It must be read for what it is: a one-parameter proxy for the total extinction curve over 0.01–0.3 μm, not a measurement of the aligned-grain size distribution, and it does not determine millimetre $\beta$. It modulates $T(a)$, $a_{\rm align}$ and $R_{P/p}$.

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 $(\beta^P_i,T^P_i)$ with a physically bounded offset instead of assuming equality or fitting freely. The width of that bound is not yet established by any single inference; the disagreement between Planck intermediate XXII, Planck 2018 XI and the SED models is itself the current width, of order 0.05–0.1 in $\beta^P-\beta^I$. Section 4.2 shows that even a width of 0.1 is far below what polarization data alone provide per pixel.
• Which layers polarize. The mass above $a_{\rm align}(U_i,n_i)$ predicts the relative $p_i$ of layers, so a dense, shielded layer with poor alignment can dominate $I$ but not $\mathcal P$. This is invisible to a 3D extinction map and is the physical origin of many $\beta^P\neq\beta^I$ sightlines.
• The optical-to-submillimetre bridge. Starlight polarization is dichroic extinction by grains of about 0.1–0.2 μm near $\lambda_{\max}$; 353 GHz polarized emission is thermal emission from the same aligned large-grain population. $R_{P/p}$ and its dependence on $R(V)$ and $U$ are size-distribution quantities. Without them, starlight tomography in Section 3.3 gives angles but not per-layer emission amplitudes.
• Bounding the moments. A plausible range of size distributions, compositions and radiation fields bounds the variance of $(\beta,T)$ within a beam and across a region, providing priors for moment and spin-moment expansions and a way to discriminate temperature from index variation, which is what Guillet et al. 2026 attempted.
• Band placement. If temperature variation dominates, high-frequency leverage matters most. Fuskeland et al. 2023 found that extending LiteBIRD's high-frequency telescope from 400 to 600 GHz reduces $\sigma(r)$ after foreground cleaning by 30–50 % when the dust model has spatially varying $\beta$ and $T$, with most of the gain reached by 500 GHz. If magnetic-dipole emission or index variation dominates, the 100–220 GHz slope carries the information instead. Grain physics plus the 3D temperature structure say which.

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, and the misalignment between layers is what turns SED differences into angle rotation and cancellation. Planck's 353 GHz map gives only the emission-weighted mean angle. Per-layer angles come from four sources, listed in order of present availability.

HI fibres. Clark et al. 2015 showed that linear HI structures trace the dust polarization angle. Clark and Hensley 2019 built position-position-velocity Stokes maps from HI4PI that correlate at $r>0.75$ with Planck 353 GHz $Q$ and $U$ with no free parameters and reproduce the dust $E/B$ ratio at $\ell\lesssim120$. Clark 2018 showed that velocity coherence traces line-of-sight tangling. BICEP/Keck XVI detected polarized dust correlated with the filamentary HI morphology template down to 95 GHz in the BK field, with $\beta=1.54\pm0.13$, found no decorrelation between filamentary and non-filamentary dust, and none from including intermediate-velocity HI. The same paper found that the template correlates with only about 10–65 % of the $B$-mode power over $20<\ell<200$. HI-based angles are therefore the all-sky, velocity-resolved angle map available now, and an incomplete cleaning template. With the velocity-to-distance assignment of McCallum et al. 2026 they can be placed in distance bins.

Starlight polarization of stars at known Gaia distances. This is the direct measurement. 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, 4 deg² observed with RoboPol, identifying several clouds out to about 2 kpc. Uppal et al. 2025 added the prior-free TRIShUL method, reliable when a cloud's polarization exceeds 0.1 % and the background-star fraction exceeds 10 % in samples of about 345 stars. PASIPHAE, described by Tassis et al. 2018, is surveying the high-latitude sky with the WALOP polarimeters for this purpose. Starlight tomography is the higher-fidelity subset and the calibrator for the HI angles.

Field models. Local Bubble models from Alves et al. 2018, Pelgrims et al. 2020 and O'Neill et al. 2024 supply the angle of the dominant high-latitude layer over the whole sky, and Galactic-field models such as Unger and Farrar 2024 cover distant layers coarsely.

Faraday rotation. Rotation measures, as reconstructed by Hutschenreuter et al. 2022, constrain $\int n_e B_\parallel\,ds$, an electron-weighted integral without dust-layer separation. The polarization angle gives the direction of $\mathbf B_\perp$ but not its magnitude, so rotation measures do not determine the inclination $\gamma_i$ of any dust layer. They can serve only as a weak prior on the line-of-sight field direction where the ionized and dusty structures coincide, as they may at the Local Bubble wall. Inclination otherwise stays latent, and Section 3.3's starlight transfer is designed so that it cancels.

On its own, the field 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 in Planck.
• Per-layer polarized amplitudes without knowing $\gamma$ or $F$. Starlight polarization is parallel to $\mathbf B_\perp$ and emission is perpendicular to it, so in a common angle convention a layer's 353 GHz Stokes vector is

with the minus sign implementing the 90° rotation. The factors $\cos^2\gamma_i$ and $F_i$ appear identically in emission and dichroic extinction and cancel, but only for a uniform, optically thin layer in which the same aligned grains produce both, sampled over the same solid angle. Stellar pencil beams and the emission beam differ, $R_{P/p}$ has been calibrated on selected diffuse sightlines rather than per layer, and it depends on environment through $R(V)$ and $U$. This is the correction Tassis and Pavlidou 2015 proposed, and it is why the field and the grain physics belong in the same sentence: the angle tomography supplies phases, and only with $R_{P/p}$ does it supply amplitudes.

• Conditioning of the layer separation. Consider two layers with externally fixed phases. The design columns for their amplitudes are $(\cos2\psi_i\,\hat g_i(\nu),\ \sin2\psi_i\,\hat g_i(\nu))$ over the bands. When the phases differ, the columns are linearly independent even for identical SEDs, so the amplitudes are well determined and a small SED difference is attributed to the right layer without the frequency lever having to do the separation. When the phases are left free, the columns for $Q_1$ and $Q_2$ are $(\hat g_1(\nu),0)$ and $(\hat g_2(\nu),0)$, nearly collinear when the SEDs are close, and the fit becomes noisy or degenerate. The Q–U plane has two dimensions, so phases alone cannot separate a third layer, but layers with identical SEDs do not need separating: their sum behaves as one layer. The geometric fact that causes the problem is also the source of the leverage, provided the phases are known.
• The pivot for spin-moments and the $E/B$ statistics. The mean field direction defines the pivot around which spin-moments are expanded, and the alignment of density filaments with the field explains dust $E/B\approx2$, 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.

4. Combining the three into a dedusting strategy

4.1 The forward model and what is joint about it

Per pixel and frequency, in specific-intensity units,

The point of the layered form is that every property of layer $i$, its weight, temperature, index, polarizing efficiency and angle, is attached to the same physical object. The correlations between a layer's SED and its angle, which are what control the residual, are built in rather than parametrized. The layer quantities and their documented status are:

The free parameters are the CMB $Q,U$, synchrotron, per-layer amplitude corrections $\alpha_i$ with priors, small angle corrections $\delta\psi_i$ with priors from the tomography errors, and a global or regional transfer offset $\Delta\beta^P$. The mean dust SED is fitted from the CMB-experiment bands themselves; the external information is needed only for the differences between layers, and that is what makes the requirements lenient.

4.2 How accurate must the external information be?

Consider template extrapolation from 353 GHz with a layered split. The true sky is $\mathcal P_\nu=\sum_i\mathcal P_i\hat g_i(\nu)$ and the measured total is $\mathcal P_{353}=\sum_i\mathcal P_i$. The model splits $\mathcal P_{353}$ into $\hat{\mathcal P}_i$ using external weights and angles, so $\sum_i\delta\mathcal P_i=0$ with $\delta\mathcal P_i=\hat{\mathcal P}_i-\mathcal P_i$, and extrapolates each layer with $\hat g_i$. The residual at 150 GHz is

for any reference $\bar g$. Every term is a product of two small quantities: an error in the split times an SED contrast, or a layer amplitude times an error in its SED. Writing $\epsilon$ for the rms split error in units of the total polarized amplitude, after the best global scalar over the patch, and $\Delta g/g$ for the fractional SED contrast between layers between 353 and 150 GHz, the residual power fraction is $f_{\rm res}\approx\epsilon^2(\Delta g/g)^2$. For $\Delta\beta=0.3$, $\Delta g/g\approx0.23$; for $\Delta T=4$ K at 16–20 K it is 0.08.

The table evaluates this for the first row of the two-layer example, where layer 2 is 30 % of the summed amplitude, misaligned by 45°, so that $|\mathcal P_2|/|\mathcal P|\approx0.4$, with $\Delta\beta=0.3$ and $r_{\rm dust}=0.1$.

Three conclusions follow. First, the requirements on external information are lenient because the residual is second order: angles to about 10°, layer amplitudes to tens of per cent, polarized weights to a factor of two and the transfer offset to 0.1 each leave $\delta r\sim10^{-4}$ where a single modified blackbody leaves $\sim10^{-3}$. This is why HI-based angles at degree resolution and 3D column fractions with factor-of-two uncertainty in polarizing efficiency already carry weight. Second, the binding requirement is the completeness of the layer inventory. A missing layer contributes with $\epsilon$ equal to its misaligned amplitude fraction, so a missing 10 % layer is tolerable and a missing 30 % layer reproduces the original problem. Third, all entries scale with $(\Delta g/g)^2$: a temperature contrast of 4 K makes them eight times smaller, while a steeper low-frequency change from magnetic-dipole emission would make them larger. These are scaling estimates for one configuration, not forecasts. A forecast must pass the actual residual map through the experiment's bands, beams, mask, filtering, delensing and likelihood, and must propagate the covariance of the external inputs; that computation has not been published for this model and is the outstanding step.

4.3 Three interfaces to an $r$ likelihood

Layered parametric fit. With externally fixed phases, each layer has one real amplitude plus a small $\delta\psi_i$ held by a prior, and the problem is linear in the amplitudes once the SED shapes are fixed. If the phases are left free the count doubles and the conditioning argument of Section 3.3 is lost. Freeing $\delta\psi_i$, $\Delta\beta^P$ or $T_i$ makes the problem nonlinear; in practice one linearizes around the external model or samples. LiteBIRD's 15 bands between 40 and 402 GHz give 30 polarized data per pixel against $2+N$ amplitudes and a few regional spectral parameters; for Simons Observatory or CMB-S4 with 6–9 bands the same model with $N\le2$ and regional rather than per-pixel spectral parameters is the realistic form. Rank is not sensitivity. Whether the fit is well conditioned in the presence of noise, bandpass uncertainty, calibration, synchrotron and angle-prior widths is a Fisher or posterior calculation with realistic inputs, and it should report the failure modes as layer SEDs and angles become similar.

Predicted spin-moment templates, and model-informed cMILC. These are two different uses of the same model. In the first, the forward model predicts the complex moment maps $\mathcal W_0$, $\mathcal W_1^\beta=\sum_i\mathcal P_i(\beta_i-\bar\beta)$ and $\mathcal W_1^T$, and a few global amplitudes of these templates are fitted alongside the CMB in FGBuster or a Gibbs sampler. The templates carry the correct spatial morphology and phase, which free per-pixel moments cannot, and their few amplitudes are far less degenerate with $r$ than free moments. The cost is that template errors bias $r$ rather than inflating $\sigma(r)$, so the residual after the template fit must be validated on non-cosmological observables. In the second use, the constrained moment ILC of Remazeilles, Rotti and Chluba 2021 deprojects prescribed moment SEDs at a variance penalty; the layered model informs which moments and which pivot to deproject and does not use the spatial maps. 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 two treatments.

Cross-spectral prior. For BICEP/Keck-style likelihoods the layer maps give the multi-frequency dust covariance

including the layer–layer cross-spectra that vanish only if the layers are spatially uncorrelated. This replaces the $\Delta_d$ ansatz and its assumed $\ell$ and $\ln^2(\nu_1/\nu_2)$ scalings with a region-specific prediction. The expression holds when each layer's SED is constant over the patch; spatial variation of a layer's SED within the patch cannot be pulled outside the harmonic transform and enters as an additional moment-type term, which is where the third effect of Section 2.2 re-enters.

4.4 A hierarchical likelihood, and what remains to be computed

The critic's central demand is correct: none of the above is a dedusting method until the external inputs are treated as noisy tracers of latent layers inside one likelihood,

with the second factor carrying the correlated radial uncertainties of the extinction posterior, the selection function and intrinsic-polarization contamination of the stellar sample, the HI-to-dust correspondence, and the ensemble of grain models. The forecast that compares a single modified blackbody, a spatially varying modified blackbody, free moments, cMILC, HI templates and the layered prior on the same simulated skies, reporting $r$ bias, variance and coverage, is the piece of work that would turn Section 4.2's scaling into a result. The scaling argument makes it worth doing; it does not replace it.

4.5 Validation on non-cosmological observables, staged use, and a falsifiable test

The forward model predicts several things that need no primordial signal and that were not used to build it: the $\psi(217)-\psi(353)$ map, which Pelgrims et al. 2021 measured statistically; the residual covariances between frequencies of Guillet et al. 2026; the frequency dependence of dust $E/B$, $TE$, $TB$ and $EB$ of Vacher et al. 2023 and Ritacco et al. 2023; cross-spectra of per-layer templates with deep-field $B$-mode maps, as in BICEP/Keck XVI; and, in the 4 deg² tomography patch and then in PASIPHAE fields, agreement between the starlight path $-R_{P/p}(q_i+iu_i)$ and the emission path $\epsilon_iA_{V,i}$ with an alignment model for each layer. Validation must use data not used for training: polarization angles rather than the Planck intensities that built the emission map, sky regions and data splits held out, and the SRoll2-versus-PR4 disagreement treated as a prerequisite to resolve rather than a footnote.

Deployment is staged. Now: risk maps and a spatially varying $\Delta_d(\hat n)$ prior from cloud counts, HI angles and 3D extinction, used for field selection and masking. As PASIPHAE, the 3D emission fits and the 3D $R(V)$ maps mature: the full per-layer model as a prior, with $\alpha_i$ and $\delta\psi_i$ free within it. With LiteBIRD and CMB-S4 data: a joint fit in which the microwave data refine $\alpha_i$, $\delta\psi_i$ and $\Delta\beta^P$ while the 3D structure remains the regularizer that prevents the moment–$r$ degeneracy.

The method makes a falsifiable prediction. In fields with starlight tomography, the layered model should reduce the 217–353 GHz polarization-angle difference on Pelgrims-selected multi-cloud sightlines relative to a single-modified-blackbody fit. If it does not, the intensity-to-polarization transfer or the layer inventory is wrong and the model must not be used to shrink $\sigma(r)$. If it does, that is the warrant for promoting predicted moment templates into the cosmological likelihood. The same test protects a second measurement: frequency-dependent dust angle rotation produces a chromatic dust $EB$ signal, while isotropic cosmic birefringence is achromatic, so a model that predicts dust $\psi(\nu)$ is also a birefringence systematic control, as Clark et al. 2021 and Diego-Palazuelos et al. 2022 emphasize.

5. Limitations, missing data and residual systematics

Limits of the 3D maps. They measure extinction, not emission. The factor-of-two variation in $\tau_{353}/E(B-V)$ found by Zelko et al. 2022 must be modelled through grain properties, and the intensity fits carry cosmic-infrared-background and zodiacal contamination and constrain the total-intensity SED rather than the large-grain SED that polarization samples. The Edenhofer map stops at 1.25 kpc, adequate at $|b|\gtrsim30°$ and not for the intermediate latitudes that LiteBIRD's reionization-bump analysis needs; the Pelgrims et al. 2024 tomography finds clouds beyond 1 kpc even at $b\approx22°$. Radial anticorrelations in the extinction posterior can create or erase apparent layers, so the posterior samples, not the mean map, must be propagated. Stellar shot noise and stellar-model systematics enter the layer weights; methods such as FLOWER by O'Callaghan, Mandel and Gilmore 2025 target the high-latitude regime specifically. Distance resolution of parsecs comes with angular resolution of 14′ to a degree, so the layered model is a tool for $\ell\lesssim200$; sub-beam structure within a layer remains a statistical residual described by a small $\Delta_d(\ell)$ at high $\ell$ and by a per-layer angle dispersion $\mathcal S_i$ informed by MHD and filament models. Predicting $T_i$ rather than fitting it needs a 3D interstellar radiation field with shielding, which does not yet exist as a product.

Limits of the grain physics. $R(V)$ is a one-parameter proxy for a size range that overlaps but is not identical with the grains that dominate millimetre polarization. Competing models, 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 polarized-to-intensity offset. Millimetre dielectric functions and their temperature dependence are poorly measured in the laboratory. The alignment threshold depends on unmeasured local conditions and the empirical alignment function is constrained only on average; the GRADE-POL work of Tram et al. 2026 shows that alignment physics is still being revised even for optical polarization spectra. $R_{P/p}$ is calibrated on selected diffuse sightlines and not per layer or per environment, and the 2.5-fold shortfall of earlier models against it shows how model dependent the transfer can be.

Limits of the magnetic-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, nothing beyond the last usable star is constrained, and the stellar selection function, intrinsic stellar polarization and the mismatch between pencil beams and the emission beam all enter the inferred layer Stokes vectors. Only 4 deg² of degree-scale tomography has been published. HI-based angles fail where the gas is molecular or ionized, rely on the HI-to-dust correspondence, and explain only 10–65 % of the $B$-mode power in the one deep field where this has been measured; velocity-to-distance assignment is ambiguous at high latitude. The inclination $\gamma_i$ is not measured by any of the tracers; it cancels in the starlight transfer and remains a source of error in the extinction-based fallback. A single $\psi_i$ per layer ignores tangling within the layer, which Halal et al. 2024 show is a significant depolarizer.

Method-level residuals. Each structural error has a definite form, and Section 4.2 gives its size: a split error leaves a layer's SED contrast in the map, an angle error leaks $E$ into $B$, a missing layer reproduces the original problem in miniature with its own amplitude fraction, and an error in $\Delta\beta^P$ is a coherent bias that LiteBIRD's bands can partly fit and Planck's cannot. Template errors bias $r$ rather than widening $\sigma(r)$, which is the price of spending external information instead of sensitivity. Bandpass mismatch and polarization-efficiency errors mimic SED curvature, and an instrumental polarization-angle error is degenerate with dust angle rotation and with cosmic birefringence. Synchrotron curvature, spatially varying dust–synchrotron correlation and any polarized anomalous microwave emission live in the same bands and must sit inside the same likelihood. The SRoll2-versus-PR4 disagreement in Guillet et al. 2026 shows that current Planck systematics are already comparable to the astrophysical signal being modelled.

What is missing. An all-sky, per-layer plane-of-sky field from starlight tomography, which PASIPHAE is building; 3D emission maps with $\beta$ and $T$ per layer validated against held-out polarization data and extended beyond 1.25 kpc; a measured link between 3D $R(V)$ and submillimetre emissivity and polarizing efficiency; $R_{P/p}$ calibrated per environment with its scatter; a settled polarized-versus-intensity SED offset, which needs LiteBIRD-class sensitivity above 300 GHz; laboratory millimetre opacities of realistic grain analogues; and the end-to-end forecast of Section 4.4. None of these gaps prevents the staged use described above. Each sets the floor on the residual that 3D dedusting can reach, and the second-order structure of that residual is the reason the floor can plausibly sit near $\delta r\sim10^{-4}$ rather than $10^{-3}$.

References

• Pelgrims et al. 2021, A&A 647, A16, arXiv:2101.09291
• Poh & Dodelson 2017, PRD 95, 103511, arXiv:1606.08922
• Ritacco et al. 2023, A&A 670, A163, arXiv:2206.07671
• Vacher et al. 2022, A&A 660, A111, arXiv:2111.07742; Vacher, Chluba, Aumont, Rotti & Montier 2023, A&A 669, A5, spin-moment formalism; Vacher, Aumont et al. 2023, A&A 672, A146, arXiv:2210.14768; Vacher et al. 2024, arXiv:2411.11649
• Guillet et al. 2026, A&A, arXiv:2510.18305
• BICEP/Keck XIII (BK18), arXiv:2110.00483; BICEP/Keck XVI, arXiv:2210.05684; Tristram et al. 2022, arXiv:2112.07961
• Planck intermediate XXI, arXiv:1405.0873; Planck intermediate XXII, arXiv:1405.0874; Planck intermediate L, arXiv:1606.07335; Planck 2018 XI, arXiv:1801.04945; Planck 2018 XII, arXiv:1807.06212
• 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; Fuskeland et al. 2023, arXiv:2302.05228; Remazeilles, Rotti & Chluba 2021, arXiv:2006.08628
• Edenhofer et al. 2024, arXiv:2308.01295; Zelko, Finkbeiner, Lee & Green 2022, arXiv:2211.07667; Zhang & Green 2025, Science; McCallum et al. 2026, arXiv:2607.07451; O'Callaghan et al. 2025, MNRAS
• O'Neill, Zucker, Goodman & Edenhofer 2024, arXiv:2403.04961; O'Neill, Goodman, Soler, Zucker & Han 2024, arXiv:2410.17341; Skalidis & Pelgrims 2019, arXiv:1908.08706; Pelgrims et al. 2020, A&A; 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; Hutschenreuter et al. 2022, arXiv:2102.01709; Diego-Palazuelos et al. 2022, arXiv:2201.07682

Cited without links, from the literature: Chluba, Hill & Abitbol 2017; Mangilli et al. 2021; Azzoni et al. 2021; Ichiki et al. 2019; Sheehy & Slosar 2018; Ghosh et al. 2017; Adak et al. 2020; Martínez-Solaeche, Karakci & Delabrouille 2018; Huffenberger, Rotti & Collins 2020; Clark et al. 2015; Panopoulou et al. 2019; Tassis et al. 2018 (PASIPHAE); Alves et al. 2018; Unger & Farrar 2024; 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; BICEP2/Keck and Planck joint analysis 2015.

What I took from the critic

The critic's most useful demand was for a quantitative accuracy requirement on the external inputs. Section 4.2 now derives one. It shows the residual is a product of a split error and an SED contrast, which is why the earlier assertion of "tens of per cent" priors can be defended rather than asserted. Section 4.4 states plainly that the end-to-end forecast has not been done and what it must contain.

Specific corrections I adopted: the three sources of decorrelation are now separated, with the Gaussian $\sigma_\beta$ formula assigned to spatial variation with its assumptions listed. The single-layer equations are labelled leading-order approximations. Faraday rotation no longer yields per-layer inclination; it is a weak line-of-sight prior at best. The starlight-to-emission transfer carries its minus sign and its conditions, and I checked and now quote the Planck XXI selection and the 2.5-fold shortfall of the models it was compared with. The parameter count and linearity claim are consistent: one real amplitude per layer with fixed phases. Template fitting of predicted moments is separated from constrained-moment ILC. The cross-spectral formula carries its regional-constant-SED condition. Zelko et al., Edenhofer et al., Ghosh et al., Adak et al. and Guillet et al. are characterized as their abstracts support, including the 1σ Cepheus result, and BICEP/Keck XVI's 10–65 % figure is included. The BICEP2 history is phrased correctly. The joint-distribution point, the radiation-field requirement, stellar selection functions, radial posterior correlations and the scale limit of the external maps are all now in the text.

What I took from the optimist

I adopted the reframing of the problem as line-of-sight vector superposition, and the observation that frequency-dependent angle, $E/B$ and $EB$ discriminate line-of-sight mixing from plane-of-sky index variation. I made the subdominant-layer point central and used Pelgrims et al.'s 15 % contrast as the empirical target rather than $\Delta\beta=0.3$. The grain-physics section now states that millimetre $\beta$ is a dielectric-function quantity while size enters through temperature, alignment and the optical-to-submillimetre bridge, with the intensity-to-polarization transfer as the central product. HI position-position-velocity Stokes maps are presented as the all-sky angle source available now with starlight tomography as the calibrator. Planck 2018 XI and Pelgrims et al. 2021 are read as complementary. The $\delta r\approx f_{\rm res}\times r_{\rm dust}$ unit and Poh & Dodelson's threshold anchor the budget. I added the first-order relation between the spin-moment phase and the angle rotation, checked against the worked example. I kept the within-layer statistical residual, moved CIB and zodiacal light into the emission step, added the chromatic dust $EB$ payoff for birefringence, and closed with the optimist's falsifiable test in tomography fields.

What I rejected, and why

From the critic, I rejected the implication that the method is blocked until identifiability is proven. The scaling argument shows the requirements are lenient enough that current external data are already informative, so the honest status is "unforecast", not "unjustified". I also kept the rotation-to-$BB$ and $1-\mathcal R^2$ estimates, with their assumptions stated, because a direct residual calculation agrees with them within a factor of two, which is adequate for a budget argument.

From the optimist, I did not put the Local Bubble two-layer cartoon in front of everything, because Halal et al. find no inclination imprint on $p$ and Pelgrims et al. 2024 find clouds beyond 1 kpc at $b\approx22°$; it is the default high-latitude geometry, not the frame. I dropped the "at least 60°" misalignment threshold, which the Pelgrims abstract does not state, and corrected the attribution of the Local Bubble field model to O'Neill et al. 2024. I omitted the Simons Observatory and CMB-S4 forecast figures I could not verify, and SPHEREx ice spectroscopy and background-galaxy polarimetry as too marginal for the cleanest fields. I did not present a 3D radiation field as closing the temperature loop, because no such product exists yet.