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

"Dedusting" is the PASIPHAE collaboration's term for removing polarized Galactic dust from CMB maps well enough to reach the tensor‑to‑scalar ratio $r\sim10^{-3}$. The obstacle that motivates the question is that the dust polarization sky is not one modified blackbody per pixel. It is a superposition of physically distinct emitting regions along and across each line of sight, and this superposition makes the effective dust spectrum, and the effective polarization angle, depend on frequency. Below I set out what this "spectral‑index decorrelation" is, how three external ingredients each constrain it, and how they combine into a forward model that can be fitted or used as a prior. I close with what is still missing.

1. Polarized thermal dust and why B‑mode searches care

Interstellar grains are aspherical and spin with their short axis aligned, on average, with the local magnetic field. The accepted mechanism is radiative‑torque alignment, as developed by Lazarian and Hoang and reviewed by Andersson, Lazarian and Vaillancourt in 2015. Their far‑infrared to millimetre emission is therefore linearly polarized perpendicular to the plane‑of‑sky field $\mathbf B_\perp$. For a single emitting region with column $N$, 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 Stokes parameters are

where $p_0$ is the intrinsic polarization of the grain population, and $F\le 1$ accounts for imperfect alignment and field tangling within the region. The polarization fraction is $p=|\mathcal P|/I$ and the angle $\psi=\tfrac12\arg\mathcal P$. Planck 2018 XII found $p$ up to about 22 % at 353 GHz, a ratio of submillimetre polarized emission to optical starlight polarization $R_{P/p}\approx5.4$ MJy sr$^{-1}$ that is nearly constant in the diffuse ISM, and an anticorrelation between $p$ and the angle‑dispersion function that signals depolarization by field tangling. Planck 2018 XI measured a mean polarized spectral index $\beta_d^P=1.53\pm0.02$ for $T_d=19.6$ K, fitted well by a single modified blackbody from 353 GHz down to below 70 GHz, dust $BB$ power falling as roughly $C_\ell\propto\ell^{-2.4}$, a $BB/EE$ ratio near 0.5, and positive $TE$ and $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 at $\ell\approx80$ is of order the signal for $r\approx0.1$, which is what BICEP2's 2014 detection turned out to be. The current state and the targets are:

So the dust model must be accurate at the per‑cent level in power, and any un‑modelled effect that leaves one per cent of dust $BB$ in the cleaned map produces $\delta r\sim10^{-3}$. That is exactly the size of the effects discussed next.

2. What spectral‑index decorrelation means operationally

2.1 The single modified blackbody per pixel

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

with the same $(\beta_d,T_d)$ for $Q$ and $U$, hence 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. Band‑pass integration and conversion to K$_{\rm CMB}$ are implied throughout.

2.2 Superposition and its four consequences

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

Four things then break.

The effective SED is not a modified blackbody. The polarized intensity $|\mathcal P_\nu|$ is a vector sum, so the effective index between any two bands, $\beta_{\rm eff}(\nu_1,\nu_2)$, depends on the pair. It can even lie outside the range spanned by the $\beta_i$ when the 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 versions.

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 showed with a two‑cloud model that the effect is large when both conditions hold: the ratio of the clouds' polarized intensities changes between the two frequencies, and their fields are misaligned. No scalar SED, however curved, can represent a rotation. It 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. Intensities add as scalars, polarizations as vectors, so $\beta_{\rm eff}^P\neq\beta_{\rm eff}^I$ arises from geometry alone, before any grain physics.

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}}$, a spatially varying index with Gaussian scatter $\sigma_\beta$ gives $R=\exp[-\tfrac12\sigma_\beta^2\ln^2(\nu_1/\nu_2)]$, and two uncorrelated layers whose relative weight changes by a factor $\rho$ between the bands, with power ratio $f$, give $R=(1+\rho f)/\sqrt{(1+f)(1+\rho^2f)}$. The residual after scaling a 353 GHz template to 150 GHz with the best single scalar is $1-R^2\approx2(1-R)$ of the 150 GHz dust power. The important arithmetic is that $R$ is a second‑order statistic: a layer carrying 30 % of the polarized amplitude with $\Delta\beta=0.3$ gives $R(217\times353)\approx0.999$, invisible to Planck, yet leaves about 0.7 % of the dust power at 150 GHz, which is $\delta r\sim10^{-3}$ in the cleanest fields.

A worked two‑layer example makes the angle and SED effects concrete. Layer 2 carries 30 % of the 353 GHz polarized intensity, both layers have the same intrinsic $p$ and the same $\beta^P=\beta^I$, and the pair of bands is 150 and 353 GHz.

A frequency‑dependent rotation $\delta\psi$ of the dust field converts $E$ into $B$: $\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$. Rotations of 1°, 2° and 3° between the template frequency and 150 GHz leave 0.2 %, 1 % and 2 % of the dust $BB$ power unmodelled. The first row of the table is therefore already a $\delta r\sim10^{-3}$ effect.

2.3 Observational status

Planck intermediate results L in 2017 reported $R^{BB}_\ell(217\times353)\approx0.95$ at $\ell\approx50$–160 in low‑column regions. Sheehy and Slosar 2018 and then Planck 2018 XI found the evidence not robust once noise bias and systematics were treated properly, and Planck 2018 XI concluded that its limit suffices for experiments targeting $r\approx0.01$, which pointedly excludes the next generation. 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 $\ell$ dependence; the data are consistent with $\Delta_d=1$ and allow at most a few per cent decorrelation. Meanwhile the line‑of‑sight mechanism has been seen directly. 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 whose fields are misaligned by at least about 60°, and found larger 217–353 GHz polarization‑angle differences there than on control sightlines. Ritacco et al. 2023, using SRoll2 maps at $\ell=4$–32 over 80–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 $EE$ but not $BB$, and identified the frequency dependence of polarization angles as the missing ingredient. 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 that SRoll2 residual covariances are compatible with temperature‑only variations, while Planck PR4 is compatible with neither temperature‑only nor index‑only variation, a warning that instrumental systematics are now at the level of the astrophysical effect.

2.4 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 a single line of sight only constrains the parameters of broad distributions, and that two‑parameter fits are biased at more than 10σ in illustrative cases precisely when both the SED and the 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 degenerate with $r$. Hensley and Bull 2018 and Sponseller and Kogut 2022 reached similar conclusions. Extra freedom costs sensitivity and does not remove degeneracy. What is needed is external information that fixes the structure of the superposition. That is what the three ingredients provide.

3. How each ingredient helps on its own

3.1 Three‑dimensional maps of thermal dust emission

By a 3D emission map I mean a distance‑resolved emissivity, built from a 3D extinction density and a per‑voxel emission law:

The extinction density now comes from Gaia parallaxes with photometry or XP spectra: 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, which uses 54 million stars to map dust to 1.25 kpc at 14′ resolution in 516 distance bins. Zelko, Finkbeiner, Lee and Green 2022 took the step to emission, proposing $T$ and $\beta$ for each of 500 million Bayestar19 voxels and integrating along the line of sight to match five Planck and IRAS bands. They resolved two clouds of distinct temperature along the Cepheus sightline and found the ratio of emission optical depth to reddening varying by a factor of two. A complementary "3D" comes from HI velocity: the HI4PI Gaussian decomposition of Panopoulou and Lenz 2020 gives about three clouds per 0.2 deg² pixel in the north Galactic cap and 2.5 in the south, and McCallum et al. 2026 now tie HI velocity components to distance by combining HI4PI with the Edenhofer map in an information‑field‑theory reconstruction. For the local geometry that matters most, Skalidis and Pelgrims 2019, Pelgrims et al. 2020 and O'Neill et al. 2024 show that the Local Bubble wall at 200–300 pc dominates the polarized sky at $|b|\gtrsim60°$, and Halal, Clark and Tahani 2024 show that the observed $p$ is lower where the column is split among components at different distances within 1.25 kpc.

On its own, this ingredient does five things for dedusting.

• 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 mask 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. The 3D emission fit gives $T_i$ and $\beta^I_i$ per layer, which fixes the shape of each layer's polarized SED up to the intensity‑to‑polarization transfer discussed in 3.2. The polarization problem then reduces to solving for $2N$ real amplitudes per pixel plus the CMB, a linear problem that current band counts can handle for $N\le3$ if the SED shapes are fixed.
• Layer templates. Per‑layer intensity maps can be fitted against $Q$ and $U$ with region‑level coefficients, the approach taken with HI velocity layers by Ghosh et al. 2017 and Adak et al. 2020, and by BICEP/Keck XVI with HI templates in the BK field.
• Realistic simulations. Multi‑layer 3D 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‑based models of Vacher et al. 2024, are how $\Delta_d(\ell;\nu_1,\nu_2)$ priors and pipeline validations are made. Real 3D structure replaces ad‑hoc layer statistics.
• Physical temperature structure. The distance of each layer from heating sources and the Local Bubble geometry predict where $T$ differs between layers, which the Guillet et al. 2026 result suggests may be the dominant source of SED variance.

The map alone does not give emission per unit extinction, nor which layers polarize, nor their angles. That is where the other two ingredients enter.

3.2 The dust grain size distribution

Modified‑blackbody parameters are not free numbers; they are set by grain physics. In the Rayleigh limit the far‑infrared opacity per unit mass is size independent, so $\beta$ is a property of the material's dielectric function. The equilibrium temperature of a grain of radius $a$ scales as $T\propto a^{-1/(4+\beta)}$, so a factor of two in size shifts $T$ by roughly 2 K. Grains smaller than 0.01–0.02 μm are stochastically heated and contribute little below 353 GHz. Crucially, 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 grain shape, as in Hoang and Lazarian 2016. The empirical alignment function of Hensley and Draine 2023 has negligible alignment below 0.01 μm and complete alignment above 0.1 μm. 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 plus alignment function plus composition. Draine and Fraisse 2009 showed that if only silicates align, $p$ rises from about 3 % at 100 μm to 15 % at 1 mm, while if both silicate and carbon grains align it is flat at 6–10 %. Planck intermediate XXII reported $\beta^P=1.59\pm0.02$ against $\beta^I=1.51\pm0.01$ and a 21±6 % decrease of $p$ from 353 to 70 GHz. Guillet et al. 2018 built models compatible with the Planck intensity and polarization constraints in translucent sightlines and showed that $R_{P/p}$ and the related emission‑to‑extinction ratios depend only weakly on column density. Hensley and Draine 2023 argue that a single composite "astrodust" for grains larger than about 0.02 μm most naturally explains the observed lack of frequency dependence of the far‑infrared polarization fraction. The discrepancy between these results is itself informative: the polarized SED offset is not settled, and further size‑ and composition‑dependent physics, 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 change the low‑frequency slope. The size distribution also varies across the sky and along the line of sight: Zhang and Green 2025 mapped $R(V)$ in three dimensions from 130 million Gaia XP spectra and find 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 Whittet and van Breda relation.

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)$ from 3.1 into $(\beta^P_i,T^P_i)$ with a physically bounded offset, instead of assuming equality or fitting freely. Even a prior of width 0.1 on $\beta^P-\beta^I$ is far more constraining than 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 the physical origin of many $\beta^P\neq\beta^I$ sightlines.
• The optical‑to‑submillimetre bridge. Starlight polarization is produced by grains of about 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}$ that converts the one into the other, and its dependence on $R(V)$, is a size‑distribution quantity. Without it, starlight tomography in 3.3 cannot be turned into per‑layer emission amplitudes.
• Bounding the moments. A physically 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 variations, as Guillet et al. 2026 did.
• Band placement. If $T$ variation dominates, high‑frequency leverage matters most; Fuskeland et al. 2023 found that extending LiteBIRD from 402 to 500–600 GHz reduces $\sigma(r)$ by 30–50 % when $\beta$ and $T$ vary on the sky. If magnetic‑dipole or index variation dominates, the 100–220 GHz bands carry the information.

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 field angle $\chi_i$ sets the phase of each term, so the misalignment between layers is what turns SED differences into angle rotation and vector cancellation. Planck's 353 GHz map gives only the emission‑weighted mean $\chi$. Per‑layer angles come from four 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 that 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, identifying several clouds; Uppal et al. 2025 added the frequentist TRIShUL method, reliable for cloud polarizations above 0.1 % with a few hundred stars; and PASIPHAE, described by Tassis et al. 2018, is surveying the high‑latitude sky with the WALOP polarimeters for exactly this purpose. HI fibres are the second source: Clark et al. 2015 showed they 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$ without free parameters and reproduce the $E/B$ ratio 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 down to 95 GHz with $\beta=1.54\pm0.13$ and found no decorrelation between the filamentary and non‑filamentary dust or from intermediate‑velocity gas in the BK field. Third, Local Bubble field 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 more distant layers. Fourth, Faraday rotation measures as reconstructed by Hutschenreuter et al. 2022 constrain the line‑of‑sight component and hence the inclination $\gamma_i$ that sets $\cos^2\gamma_i$.

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 in Planck.
• Per‑layer polarized amplitudes without knowing $\gamma$ or alignment. Because the same aligned grains produce optical dichroic extinction and submillimetre polarized emission, the starlight $q_i,u_i$ of a layer multiplied by $R_{P/p}$ gives that layer's 353 GHz $Q_i,U_i$ directly; $\cos^2\gamma_i$ and $F_i$ cancel in the transfer. This is the correction Tassis and Pavlidou 2015 proposed.
• The pivot for spin‑moments and the $E/B$ statistics. The mean field direction defines the pivot angle around which spin‑moments are expanded, and the alignment of filaments with the field explains the 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.
• 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 velocity‑distance 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:

The free parameters are the CMB $Q,U$, synchrotron, per‑layer amplitude corrections $\alpha_i$ with priors of a few tens of per cent, small angle corrections $\delta\psi_i$ with priors from the tomography errors, and a global or regional transfer offset $\Delta\beta^P$. Everything that produces frequency‑dependent angles and SED curvature is now generated by the model rather than parametrized.

4.2 Three ways to use it

Layered parametric fit. With LiteBIRD's 15 bands from 40 to 402 GHz there are 30 data per pixel, against $2+2N+2$ amplitudes and a handful of regional spectral parameters. With fixed layer SEDs the problem is linear in the amplitudes, and its conditioning comes less from the small SED differences than from the different $Q/U$ ratios of misaligned layers, which is the same fact that causes the problem in the first place. 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 feasible.

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 spin‑moments are complex, and their phase relative to $\mathcal W_0$ encodes the angle 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 as in Remazeilles, Rotti and Chluba 2021, retains the correct functional form while breaking the moment–$r$ degeneracy that free moments suffer from. 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.

Cross‑spectral prior. For BK‑style likelihoods, the model gives the full multi‑frequency dust covariance $C_\ell^{ab}=\sum_{ij}\hat g_i(a)\hat g_j(b)\,C_\ell^{ij}$ from the layer maps, which replaces the $\Delta_d$ ansatz and its assumed $\ell$ and $\ln^2(\nu_1/\nu_2)$ scalings with a region‑specific prediction.

4.3 Validation and staged use

The forward model makes predictions testable with existing data: the $\psi(217)-\psi(353)$ map, which is what Pelgrims et al. 2021 measured statistically; the residual covariances between frequencies 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 as in BICEP/Keck XVI. The deployment is naturally staged. Now: risk maps and $\Delta_d$ priors from cloud counts, HI angles and 3D extinction, used for field selection and masking. As PASIPHAE, 3D emission and 3D $R(V)$ maps mature: the full per‑layer model as priors. With LiteBIRD and CMB‑S4 data: joint fits with the structural priors, in which the data themselves refine $\alpha_i$, $\delta\psi_i$ and $\Delta\beta^P$. A useful sanity check of the payoff: if the 3D map shows one layer carrying 90 % of the polarized column and the tomography shows the remainder 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$. Where the maps show otherwise, the model supplies the correction instead of an inflated error bar.

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. Distance resolution of parsecs comes with angular resolution of 14′, and the Edenhofer map stops at 1.25 kpc, which is adequate at $|b|\gtrsim30°$ but not for 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 by O'Callaghan, Mandel and Gilmore 2025 are addressing the mid‑ to high‑latitude regime specifically. The 3D emission fits use Planck and IRAS intensities that carry cosmic‑infrared‑background and zodiacal contamination, and they constrain the total‑intensity SED, not the large‑grain SED that polarization samples. Within a layer, turbulent tangling of the field and sub‑beam variation of $(\beta,T)$ are not resolved, so a residual decorrelation at high $\ell$ and a residual depolarization remain to be described statistically.

Limits of the grain physics. $R(V)$ is a one‑parameter proxy for grains of 0.01–0.3 μm, while the millimetre polarized emission is dominated by grains above 0.1 μm. 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‑wave dielectric functions and their temperature dependence are poorly measured in the laboratory. Alignment theory gives a threshold that depends on unmeasured local conditions, and the empirical alignment function is only constrained on average. The recent GRADE‑POL work of Tram et al. 2026 shows the alignment physics is still being revised even for optical polarization spectra.

Limits of the magnetic field data. Starlight tomography needs enough polarized stars behind each layer; at high latitude the 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 has been published. HI‑based angles fail where the gas is molecular or ionized and rely on the HI‑to‑dust correspondence, and the mapping of HI velocity to distance is ambiguous at high latitude. The inclination $\gamma_i$ is only indirectly available, through rotation measures that weight $n_e$ rather than dust, and it matters for the layer weights in intensity even though it cancels in the starlight transfer. A single $\psi_i$ per layer ignores tangling within the layer, which Halal et al. 2024 show is a significant depolarizer.

Method‑level residuals. Every structural error maps onto a residual with a definite form: an error $\delta\alpha_i$ leaves a fraction of layer $i$ at the wrong SED, an error $\delta\psi_i$ leaks $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 that LiteBIRD's bands can partly fit but Planck's cannot. Adding layers costs noise unless their SEDs and angles are well separated. Band‑pass mismatch and polarization‑efficiency errors mimic SED curvature, and an instrumental polarization‑angle error is degenerate with a dust angle rotation and with cosmic birefringence, as Ritacco et al. 2023 and Clark et al. 2021 both emphasize. The disagreement between SRoll2 and PR4 in Guillet et al. 2026 shows that current Planck systematics are already comparable to the astrophysical signal being modelled. Synchrotron–dust correlation and any polarized anomalous microwave emission are outside the dust model but inside the same bands.

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 at CMB‑relevant resolution and beyond 1.25 kpc; a measured link between 3D $R(V)$ and submillimetre emissivity and polarization efficiency; a settled polarized‑versus‑intensity SED offset, which needs LiteBIRD‑class sensitivity above 300 GHz; the line‑of‑sight field component per layer; and laboratory millimetre opacities of realistic grain analogues. None of these gaps prevents the staged use described in section 4, but each one sets the floor on the residual that a 3D dedusting 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, spin‑moments, A&A; Vacher et al. 2023, $E/B$ and $EB$, 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 L, arXiv:1606.07335; Planck 2018 XI, arXiv:1801.04945; Planck 2018 XII, A&A; Planck intermediate XXII, arXiv:1405.0874
• 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
• Edenhofer et al. 2024, arXiv:2308.01295; Zelko et al. 2022, arXiv:2211.07667; Zhang & Green 2025, Science; McCallum et al. 2026, arXiv:2607.07451; O'Callaghan et al. 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 without links, from memory of the literature: 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; Hutschenreuter et al. 2022; 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.