Preserved original model output. Check claims against the independent audit →
Dedusting CMB B-modes with 3D dust, grain-size physics, and the sky-projected magnetic field
Polarized thermal emission from Galactic dust is the leading astrophysical contaminant of primordial CMB B-modes at frequencies . Component-separation pipelines that treat each pixel, or each power-spectrum bin, as a single modified blackbody (MBB) with a frequency-independent polarization angle fail when the effective dust spectral energy distribution (SED) and/or the polarization angle become frequency-dependent. That failure is what the CMB literature calls spectral-index (or frequency) decorrelation.
The physical origin is a superposition effect, not a grain-by-grain mystery. The observed Stokes parameters are integrals of emitters that need not share a common temperature , opacity index , grain mix, alignment efficiency, or magnetic-field orientation. Frequency-dependent net angles are not automatic. They require a covariance between SED differences and polarization orientation: different angles with identical SEDs leave the net angle fixed; different SEDs with identical angles change the effective index but do not rotate the polarization. The controlling object is therefore the joint, polarized-emissivity-weighted distribution
especially covariances such as and .
Three independent observational ingredients—3D maps of the dust column (and, where they exist, of dust temperature), knowledge of the grain size distribution and alignment physics, and maps of the plane-of-sky magnetic field—do not magically identify that joint distribution. They supply independent, CMB-free constraints on its factors: distance-resolved columns, a mapping from local conditions to polarized SEDs, and orientation proxies per layer. Used together they convert an unconstrained 2D extrapolation into a layered forward model with informative priors, a mask of the worst sightlines, and predicted residual moments. They do not, by themselves, make the inverse problem fully determined, and they have not yet been shown end-to-end to deliver unbiased at . This essay states that argument operationally: what is being measured, why a one-MBB-per-pixel model is misspecified, how each ingredient contributes on its own, how they combine into a concrete dedusting strategy, and where the residual systematics remain.
Conventions. Stokes parameters are written in the IAU convention: the electric-vector position angle (EVPA) is measured from north toward east, and
Planck polarization products used for spherical-harmonic analysis are delivered in the COSMO convention, ; that sign must be applied before quoting , , or . The observed polarization fraction is biased by noise on and by the intensity zero level; Planck analyses use debiased estimators. For optically thin thermal dust, the emission EVPA is orthogonal to the plane-of-sky projection of . Starlight polarization is parallel to , and H I filaments are statistically parallel to . Assigning an H I or starlight magnetic angle directly to the emission Stokes equations without a rotation flips the signs of both and . Throughout, denotes the emission EVPA unless is written explicitly; they differ by . The geometric factor , with the angle between and the plane of the sky, is the projection of aspherical grain emission (Lee & Draine 1985). A projected orientation measures , not .
The dust intensity of a homogeneous volume element is written as an MBB
with the Planck function and a conventional reference . Polarized emission from the same volume, in IAU coordinates, is
where is the intrinsic polarizing efficiency of the aligned grain population. Equivalently, one may normalize at the reference frequency,
with the polarized intensity at . These two amplitude conventions produce different temperature-moment kernels, and they must not be mixed (see §2.3). As soon as several such volume elements occupy the same beam, the effective SED of , , and need not remain an MBB. The net EVPA becomes frequency-dependent only if the layers differ in both SED and .
1. Polarized thermal dust and why B-mode / tensor-to-scalar-ratio studies care
Inflationary gravitational waves source a curl (B-mode) pattern in CMB polarization on degree scales, with amplitude parameterized by the tensor-to-scalar ratio at a pivot . The present observational bound is at 95% confidence from BICEP/Keck through the 2018 season (BK18), tightening to when Planck PR4, BAO, and CMB lensing are combined (BICEP/Keck Collaboration 2021; Tristram et al. 2022). LiteBIRD targets a total uncertainty , including foreground residuals; CMB-S4 aims at a detection of in its deep patch. At these levels the Galactic polarized foreground is not a small correction. In the BICEP/Keck field, Galactic dust B-modes at already exceed the tensor spectrum (BICEP2/Keck Collaboration 2018). Over most of the high-latitude sky the dust B-mode power at is larger still.
Thermal dust emission is the long-wavelength tail of radiation from large interstellar grains in equilibrium with the interstellar radiation field. The mean SED in both intensity and polarization is well approximated by an MBB over the Planck HFI range. Planck intermediate results XXII found
for a mean colour temperature on intermediate-latitude patches (Planck Collaboration Int. XXII 2015). The same paper reported that the dust polarization fraction decreases by from to , inconsistent with a single grain population of frequency-independent . Planck 2018 XI, using PR3 maps and end-to-end simulations, quotes a slightly lower polarized index, at the same , and a polarized-minus-intensity difference
That offset is only : it is a hint, not an established detection, that the grains dominating are not identical to those dominating and . The later, more conservative number is the one that should be carried into a B-mode analysis.
The polarization itself is produced because aspherical grains align with their short axes preferentially parallel to . Radiative alignment torques (RATs) are the leading mechanism in the diffuse ISM (Dolginov & Mitrofanov 1976; Draine & Weingartner 1996, 1997; Lazarian & Hoang 2007; Andersson, Lazarian & Vaillancourt 2015). Alignment is strongly size-dependent: only grains larger than a critical alignment radius (typically in the diffuse ISM) rotate suprathermally and remain aligned against gas randomization. Consequently the polarized submillimetre SED is weighted toward the large-grain tail of the size distribution, whereas total intensity also receives emission from smaller, hotter, poorly aligned grains. That is a population statement. It does not, by itself, fix or on a given sightline: those also depend on composition, optical constants, and the radiation field (see §3.2).
Planck 353 GHz maps remain the highest-signal-to-noise full-sky tracer of polarized dust. They reveal several statistical facts that any B-mode foreground model must reproduce:
- A maximum observed polarization fraction – at (Planck Collaboration XII 2018), a lower bound on the intrinsic grain polarizing efficiency after line-of-sight depolarization. The uncertainty on is dominated by the intensity zero level.
- An power asymmetry on intermediate scales, together with a positive correlation, both of which arise because filamentary density structures are preferentially aligned with (Planck Collaboration Int. XXXVIII 2016; Clark et al. 2015). Planck 2018 XI found and on its largest high-latitude region, and reported a positive dust .
- A tight anti-correlation between and the polarization-angle dispersion function , consistent with depolarization by a turbulent field rather than by spatially variable grain physics (Planck Collaboration Int. XIX 2015; Planck Collaboration XII 2018).
- Residual spatial variation of the polarized SED at a level that could matter for B-mode analyses, but whose detection is not uniform across data sets (Planck Collaboration Int. L 2017; Planck Collaboration XI 2020; Ritacco et al. 2023; see §2.2).
Component separation for is usually done in one of three families: (i) parametric pixel-domain fits (Commander, FGBuster), in which each pixel is an MBB plus a synchrotron power law plus CMB; (ii) internal linear combination / needlet ILC methods, which enforce unit response to a blackbody CMB and minimize variance—ordinary ILC does not impose a dust-null constraint, though constrained or moment-deprojected ILCs add explicit foreground constraints (Remazeilles, Rotti & Chluba 2021); (iii) template multi-frequency likelihoods at the power-spectrum level, as in BICEP/Keck, which scale a 353 GHz dust template with an MBB and a single (sometimes with a decorrelation parameter ). All three families become biased if the true dust SED is not a single MBB, or if and do not share the same frequency scaling. Spatially varying mixing also increases ordinary-ILC residuals and can induce ILC bias, but the mechanism is variance-weighted leakage of a non-blackbody foreground, not the nulling of one assumed dust SED. That is the operational problem of spectral-index decorrelation.
2. What spectral-index decorrelation means operationally
2.1 Two distinct, observationally mixed phenomena—and the covariance that couples them
The phrase “spectral-index decorrelation” is used in two related but not identical ways.
Spatial (2D) decorrelation. The dust SED varies across the sky. If or is not constant, a map at frequency is not a scalar multiple of the map at . At the power-spectrum level this is quantified by the correlation ratio
A single MBB with uniform implies . Under additional assumptions—Gaussian fluctuations that are statistically independent of the amplitude field, and a map that is coherent on the scales contributing to a given —one obtains the familiar expansion
(Planck Collaboration Int. L 2017). That formula is not generally applicable. Multiplication by a spatially varying field couples harmonic modes; alone does not determine . Spatial correlations between and the dust amplitude, especially at higher orders, can bias even when a low-order moment expansion is fitted (Liu et al. 2025, in a Simons Observatory forecast). Those correlations are physical: colder, denser, higher-column layers have different temperatures and indices. A random-field moment model treats them as nuisance variance; a layered model with per-layer physics can encode them, if the layers are correctly assigned.
Line-of-sight (LOS) decorrelation. Even a single pixel can fail to be a single MBB. If the LOS intercepts two or more clouds with different SEDs and misaligned emission angles, the Stokes parameters add as complex amplitudes:
Rotation of the net EVPA occurs only when both varies with frequency and . A third, often omitted, condition controls the size of the rotation: the polarized amplitudes must be comparable. A dominant cloud produces little rotation regardless of the second cloud’s angle; large rotations, including the extreme flip when crosses unity at , live near cancellation. This is the Tassis & Pavlidou (2015) mechanism. It is not equivalent to a spatially varying but frequency-independent per pixel. A pixel-wise MBB fit still assumes that and share a common SED, i.e. that is independent of . LOS mixing of misaligned, spectrally distinct clouds violates that assumption: the effective spectral indices of and of differ. Parametric component separation that fits a single to jointly therefore leaves a frequency-dependent residual whose content is not that of the CMB.
In the spin-moment language of Vacher et al. (2023a, A&A 669, A5), the net field carries a complex spectral index whose imaginary part is the frequency derivative of the angle:
A single MBB with a fixed angle cannot represent this. Neither can a real-valued intensity moment expansion. Scalar moments in and capture SED distortions of (and of if the angle is fixed); they do not represent residual frequency-dependent polarization angles. Polarized residual moments must be spin-2, generally complex quantities.
An illustrative two-cloud calculation is useful for scaling, not as a typical-sky claim. Two clouds with equal polarized intensity at 353 GHz, temperatures and , indices and , and emission angles misaligned by produce a net-angle rotation of a few degrees from 353 to 150 GHz and a naive residual of a few percent of the dust power after a 353 GHz template is scaled by a single MBB. Frequency-dependent rotation also leaks the brighter dust modes into , so the residual can exceed that naive fraction. Poh & Dodelson (2017) found that temperature-only extrapolation noise, using Planck-compatible temperature distributions, is more than an order of magnitude below then-current errors: possible extreme configurations are not the expected sky average. The B-mode question is whether a subset of mixed, misaligned sightlines, with realistic occurrence rates, can still bias at . That is a question about the tails of , not about the mean SED.
2.2 Observational status: detected in a selected sample, not detected as a sky-average
The empirical status must be stated sharply, because “decorrelation is a problem for ” and “decorrelation has been measured on the B-mode sky” are different claims.
Planck large-region analyses do not detect a loss of correlation. Planck Collaboration Int. L (2017) reported departing from unity at in some high-latitude regions, with an indicative . Planck 2018 XI, using PR3 maps and end-to-end simulations, found that systematics and noise had overstated that evidence, in agreement with Sheehy & Slosar (2018). For the multipole range , the multi-frequency MCMC fit of Planck XI gives on the largest region (LR71, ), with an end-to-end lower limit of ; the two-frequency 217-versus-353 analysis is weaker (E2E lower limit on LR71, on LR24). Planck XI explicitly cautions that its sensitivity “precludes identifying how difficult the component-separation problem will be for more ambitious experiments targeting lower limits on ”. If the LR71 limit applied to sub-orbital patches, decorrelation might not be a problem at ; it does not follow that it is negligible at .
BICEP/Keck finds no decorrelation in its deep field. BK18’s decorrelation parameter is consistent with unity, with a statistical error . BICEP/Keck XVI (2023) correlated BICEP/Keck and Planck maps with H I morphology templates, detected filamentary polarized dust down to 95 GHz, measured for the filamentary component, and found no evidence for decorrelation between the filaments and the rest of the dust field, nor from including intermediate-velocity H I. The H I morphology template correlates in modes at a – level over with maps that also contain CMB and noise—not a 75% per-layer reconstruction of the dust modes used for .
SPIDER finds no significant full-region LOS-decorrelation detection, using a Gaussian-mixture test, yet measures spectral indices in two sub-regions that differ at , and a joint fit with Planck that rejects a simple MBB in one sub-region (SPIDER Collaboration 2025). Spatial SED variation is present; a full-patch LOS-angle effect is not detected.
Pelgrims et al. (2021) detected the LOS mechanism in a selected sample. Using H I to count clouds per sightline (Panopoulou & Lenz 2020) and H I filament orientation to estimate per cloud (Clark & Hensley 2019), they selected a target sample of high-latitude pixels with and cloud–cloud misalignment , and a control sample of single-cloud sightlines. The target sample shows a statistically significant excess of 353-to-217 GHz polarization-angle differences and extra depolarization. The probability that the effect is due to noise alone ranges from to , depending on the CMB map and the treatment of residual systematics; correcting residual systematics in SRoll2 increases the significance. The effect is stronger for more misaligned fields, as the two-cloud model predicts. An intrinsic cloud-to-cloud variation of in the 353/217 polarized-intensity ratio is sufficient to reproduce the measured effect. That is all the paper claims. It does not recover individual-cloud , , grain sizes, or 150 GHz residuals, and it does not convert 15% in a band ratio into a 10% “SED shape” change.
Pixel-level and residual-map studies leave room at high latitude. Osumi et al. (2021) find that Planck polarization at scales is consistent with a uniform in -of-sky template fits, but that pixel-level multi-frequency fits still allow a Gaussian up to at the highest latitudes, at mid-latitudes, and at low latitudes. Ritacco et al. (2023) find effective angle residuals of order a degree at 217 GHz and several degrees at 100 GHz at the largest scales after mean-SED subtraction. Guillet, Vacher et al. (2026) analyse residual-map covariances in Planck polarization with spin moments: for SRoll2 the spectral dependence is compatible with temperature variations alone rather than index variations; for PR4 it is incompatible with both simple models. Temperature per layer is the quantity a 3D density plus radiation-field model most directly predicts; that does not make temperature the unique residual in every data release.
The CMB data therefore bound sky-average decorrelation at the few-percent level and cannot themselves reach the sub-percent control next-generation experiments need. Twenty-fold tighter control has to come from somewhere else—or from more frequencies, more moments, and more sky. The three ingredients discussed here are a candidate for that “somewhere else”, not a demonstration that it has been achieved.
2.3 Why a single MBB per pixel fails—and when it does not bias CMB
The MBB is a nonlinear function of and . Averaging MBBs—along the LOS, inside the beam, or when computing spherical-harmonic coefficients over a large mask—does not return an MBB. Chluba, Hill & Abitbol (2017) formalized the resulting SED distortions as a moment expansion about a fiducial . The kernels depend on the amplitude convention. If the fundamental SED is
the first temperature derivative is proportional to , with
If instead is a reference-frequency intensity and the SED contains , the first temperature derivative involves . Mixing these conventions produces an internally inconsistent expansion. In the Rayleigh–Jeans regime relevant to CMB channels, and are highly degenerate, so a fit that expands only in can absorb some variance—but not enough at LiteBIRD sensitivity (Vacher et al. 2022).
For polarization the expansion must be spin-2 (Vacher et al. 2023a). Each spin moment has real and imaginary components; the first-order complex moment already describes frequency rotation of . Vacher et al. (2023b) showed that this induces a frequency dependence of the dust ratio and of . Vacher, Carones et al. (2025) showed that second-order spin moments reproduce the six-layer PySM d12 model to a few parts in across LiteBIRD bands, and that parametric and ILC methods respond to different moment properties. A layered forward model should therefore be delivered to likelihoods as predicted moment maps with covariances, not only as band templates.
Operational consequences for component separation:
-
Template methods. Scaling a 353 GHz template by a single MBB (or even by a spatially varying map derived from intensity) mis-predicts the dust at 90–150 GHz wherever or the effective of polarization differs from that of intensity. The residual is a B-mode field that is spatially correlated with the Galaxy and spectrally intermediate between dust and CMB. Frequency-dependent angle residuals leak dust into .
-
Pixel-domain parametric fits. A per-pixel MBB has three dust parameters if is fixed, or four if is free, plus CMB . With frequency bands this looks overconstrained. It is not: the model is misspecified. McBride, Bull & Hensley (2023) showed that when the true LOS is a distribution of MBBs, the best-fit is biased away from both the mean and the median of the true distribution—so the physical interpretation of fitted dust parameters is generally wrong. Recovered CMB , however, remain close to unbiased unless SED and polarization angle vary together along the LOS. Their intermediate example gave only CMB-amplitude bias; the large () bias was an illustrative turbulent-angle model. That distinction is the whole point of bringing in per layer. Hensley & Bull (2018) found that LOS mixing and iron grains are the two physical complexities that most stubbornly bias recovered CMB Stokes parameters, even when a generalized MBB with extra phenomenological parameters is fitted.
-
Harmonic-space likelihoods. Fitting with a single and a single decorrelation parameter (Vansyngel et al. 2017) can absorb a mean loss of coherence but cannot represent frequency-dependent or the distinct and moment spectra. Mangilli et al. (2021) showed that ignoring SED distortions, or replacing them by one , biases above the sensitivity of next-generation experiments. Vacher et al. (2022) found that LiteBIRD-like simulations with spatially varying produce a bias if a single MBB is assumed; expanding in alone is still insufficient; expanding in both and recovers an unbiased with , at the cost of moment–CMB correlations that inflate the error. Priors on the moment maps from a 3D forward model are a candidate way to shrink that inflation; they have not been demonstrated.
-
Blind ILC. Ordinary ILC enforces unit CMB response and minimizes variance. Frequency-dependent mixing leaves a dust residual that can leak into the CMB solution. Moment-constrained ILC (Remazeilles, Rotti & Chluba 2021) deprojects the first few moment templates, at the cost of extra noise and of needing enough frequency channels. Predicted spin-moment maps from a layered model are the natural templates to deproject.
In short: spectral-index decorrelation means that the dust sky at 150 GHz is not in general a rescaling of the dust sky at 353 GHz, either in amplitude or in polarization angle, and that the residual after a one-component subtraction has a B-mode spectrum that can be degenerate with . The rest of this essay is about reconstructing the missing layers rather than fitting their average—while remaining explicit about what the ancillary data actually measure.
3. How each ingredient helps on its own
3.1 3D maps of thermal dust emission (and the 3D extinction and temperature maps that currently exist)
A 2D dust map, however high its S/N, is a column integral. Decorrelation is generated by differences among the emitters that were integrated, coupled to differences in orientation. A 3D map converts that integral into a sum of discrete layers, each of which can be assigned its own SED and its own . That is the information gain. It is not a measurement of layer-resolved polarized amplitudes or of .
What actually exists today are primarily 3D extinction / reddening maps, plus one large-scale 3D dust-temperature reconstruction, not 3D maps of polarized submillimetre emission. They are nonetheless the right starting point, because in the diffuse ISM dust mass, optical extinction, far-infrared emission, and H I column are tightly correlated (Boulanger et al. 1996; Lenz, Hensley & Doré 2017). Gaia parallaxes plus stellar photometry and low-resolution BP/RP spectra have produced a sequence of maps with steadily improving distance resolution:
- Green et al. (2019), Bayestar19: Gaia + Pan-STARRS 1 + 2MASS, northern sky, Gaussian-process prior on the dust density.
- Leike, Glatzle & Enßlin (2020) and subsequent information-field-theory reconstructions: parsec-scale structure in the local ISM.
- Edenhofer et al. (2024): 54 million stars with Gaia XP-based distances and extinctions (Zhang, Green & Rix 2023), reconstructed with variational inference (MGVI) on a spherical grid from 69 to 1250 pc in 516 logarithmically spaced distance bins, angular resolution up to (). This map resolves the internal structure of hundreds of local molecular clouds at roughly the angular scale of CMB B-mode analyses. The 516 radial bins are grid sampling of a regularized posterior, not 516 independent measurements along every LOS; the released posterior has strongly correlated spatial uncertainties. Low-extinction high-latitude directions are not automatically the best-resolved layer decompositions merely because much of their dust is local.
- Zhang & Green (2025) mapped the extinction-curve parameter in 3D from Gaia XP spectra of million stars. varies with extinction in a way consistent with accretion-dominated growth in moderately dense gas and coagulation at higher density. That is a 3D grain-size proxy, not a microwave emissivity map. Through the Serkowski relation the optical extinction curve also tracks the size of the aligned grains, which links to the RAT alignment threshold—again as a prior, not as a per-layer .
- Zelko, Finkbeiner, Lee & Green (2023) assigned a temperature and index to each Bayestar19 voxel, integrated to match Planck and IRAS emission, and produced the first large-scale 3D dust-temperature map (at , , and ). They resolved distinct temperatures for two Cepheus components on one sightline. This is the prototype of the “3D emission map” the question asks about. It also constrains the optical-depth-to-reddening conversion per region, which they find varies by a factor of , contrary to the constant conversion assumed in many emission-based reddening maps. The reconstruction is constrained by integrated FIR emission: layer temperatures are inferred, not independently measured, and they inherit the usual – degeneracy.
At high Galactic latitude the dust that dominates 353 GHz emission lies mostly within a few hundred parsecs—inside the volume of these maps. That is a crucial practical point: B-mode fields are high-latitude fields, so the 3D dust maps are most complete precisely where CMB experiments need them. Completeness is not the same as layer identifiability. Intermediate- and high-velocity clouds can still contribute polarized emission at high and are poorly sampled by stars.
The Local Bubble is a candidate for the dominant first layer, not a solved geometry. Skalidis & Pelgrims (2019) used stars as distance candles and found the 353 GHz polarized sky at high southern and northern latitudes dominated by a magnetized structure at roughly the Local Bubble shell distance. O’Neill et al. (2025) built a 3D field model on the shell from Planck 353 GHz and a 3D-dust-derived shell geometry, under the assumptions that the shell dominates the polarized signal along most lines of sight and that the field is tangent to the dust-traced surface; they found the Bubble likely dominates the polarized signal along most sightlines, and they validated the reconstruction on a simulated superbubble. Maconi et al. (2023) show from MHD that such a shell can act as a translucent polarization filter for emission behind it. Halal, Clark & Tahani (2024), using Edenhofer et al. (2024) and Planck 353 GHz, found no imprint of the tangential-field assumption on the polarization fraction. They did find that, at fixed column, the polarization fraction decreases when the dust column is distributed among multiple components at different distances—but only when dust out to is included; the effect disappears if the map is truncated at the shell distance (). The extended 3D distribution, not solely the Bubble, sets at . The Bubble is therefore a physically motivated first layer whose orientation can be given a geometric prior; its inclination and its dominance are still fitted quantities.
What a 3D column map gives a component-separation pipeline, even before grain physics or are added:
- A physically motivated layering. Instead of one MBB per pixel one writes
with related to the differential extinction in distance bin . The conversion is not a universal scale factor. If is in magnitudes, the V-band optical depth is ; going from to requires the ratio of submillimetre opacity to optical extinction, which depends on size distribution, composition, porosity, and , plus an explicit dust-to-gas convention. An overall scale can be calibrated on the high-latitude – correlation and on the – relation of Lenz et al. (2017), which holds for with scatter. Layer-to-layer variations of that conversion, if correlated with or , are a residual at the level. PySM model d12 already implements a six-layer phenomenological version of this idea and is used as the “high-complexity” dust model for CMB-S4 / LiteBIRD / SO forecasts. The 3D extinction map replaces those phenomenological layers with data-driven ones, subject to the conversion systematic.
-
A prior on where LOS complexity exists. Pelgrims et al. (2021) needed H I cloud counts as a proxy for multi-cloud sightlines. A 3D extinction map provides related information in true distance, including clouds that share a radial velocity (and would be blended in H I) but sit at different distances. It does not provide “the same information” as an H I cloud count: H I also carries velocity-coherent magnetic orientation, and 3D extinction does not. Sightlines with a single dominant cloud can be treated with a one-MBB model; sightlines with two or more comparable columns are the ones that must be given extra SED freedom, or masked, or modelled as a sum of MBBs. Unimodal extinction versus distance does not imply small or temperature moments: a single cloud may contain illumination gradients, composition variation, and beam-scale mixing. Conversely, multiple distance peaks with identical polarized SEDs and aligned fields need not produce decorrelation.
-
A conversion from extinction to emission, once is assigned. Optical/NIR extinction measures dust mass weighted by the optical extinction cross-section. Submillimetre emission measures dust mass times . Zelko et al. (2023) is the existing implementation of this conversion in 3D. Residuals between the predicted integrated emission and Planck 353 GHz then constrain remaining opacity and temperature variations, rather than having to infer the entire 3D structure from the SED alone. Polarized CIB is expected to be small but is not measured at the required level; CIB anisotropies further contaminate intensity-based and maps at high latitude.
-
A handle on which moments should be large. Moment expansion treats all averaging as unknown. A 3D map tells the pipeline where moments are expected to be large: pixels that pierce two clouds with different local densities (hence different ) and different field orientations. It does not set the moments to zero.
Planck Collaboration Int. XVII (2014) already found, from correlating dust with H I velocity components over the southern cap, that temperature is anti-correlated with emissivity and that an extra component at is of order a quarter of the emission, attributable to magnetic dipole emission, carbon dust, or flattening silicate opacity. That is direct evidence that layers differ in SED, and that the low-frequency polarized SED is not settled by an MBB.
3.2 Dust grain size distribution
The grain size distribution and the size-dependent alignment function control three quantities that a single MBB cannot represent: the frequency dependence of , a possible difference between and , and the coupling between local physical conditions and the SED of each 3D layer. They do not uniquely predict those quantities. Composition, shapes, porosity, aggregation, dielectric functions, iron-inclusion abundance, and their environmental variation are missing degrees of freedom. Hensley & Draine (2023) argue that a single “astrodust” composition for explains a polarization fraction flat within about 10% from to 3 mm and attributes the optical-to-FIR polarization ratio to grain shape (oblate axial ratio ). The existence of that viable one-population model is the demonstration that does not uniquely predict or .
Size distribution. The classical MRN distribution from to (Mathis, Rumpl & Nordsieck 1977), and its descendants (Weingartner & Draine 2001; Zubko, Dwek & Arendt 2004; Compiègne et al. 2011; Jones et al. 2017), put most of the dust mass in grains with and most of the number in nanoparticles. Far-infrared and submillimetre emission is dominated by the big-grain population, which is in thermal equilibrium. Very small grains and PAHs are stochastically heated and dominate the mid-infrared; they are at most weakly aligned and contribute negligibly to polarized submillimetre emission. They are not a superposition of equilibrium MBBs. For large grains in the long-wavelength Rayleigh regime, is controlled strongly by material dielectric response, disorder, inclusions, and structure—not by alone. The practical CMB-relevant statement is that polarized dust is a large-grain phenomenon, while total-intensity dust at can include additional, less polarized, components (AME from spinning nanoparticles; possibly magnetic dipole emission from ferromagnetic inclusions; a carbonaceous population with a different ).
Alignment is a high-pass filter in grain size. In RAT theory the torque efficiency scales as a broken power law in (Lazarian & Hoang 2007). Grains align stably when they rotate suprathermally, which defines a minimum alignment size . The actual approximation (Hoang, Lazarian & Schlickeiser 2021; Hoang & Lazarian 2014, 2016) depends on normalized gas density , gas temperature , mean radiation wavelength , radiation energy density , radiation anisotropy , and an IR-damping factor , as well as grain helicity and magnetic susceptibility. Magnetically enhanced RAT alignment (MRAT) with iron inclusions can drive high- attractors for a broader range of shapes. Radiative torque disruption can set in strong radiation fields. A 3D extinction map, even with a fiducial cloud thickness to estimate and a -based estimate of , does not supply , , or the magnetic properties. Using to infer while claiming to predict is circular. Alignment modelling is becoming statistically testable (Tram and collaborators have reproduced Serkowski and super-Serkowski starlight spectra with RATs plus paramagnetic relaxation under specified radiation fields), but it is not a cloud-by-cloud 3D alignment map at high latitude.
Why this matters for decorrelation.
-
Frequency-dependent . If two grain populations contribute to but only one is aligned, is the ratio of a polarized MBB to a sum of (possibly non-MBB) components and is not constant. This is a natural interpretation of the Planck Int. XXII decline of from 353 to 70 GHz, and of a possible offset: silicates (or astrodust) aligned, a carbonaceous or PAH-related component unaligned (Draine & Fraisse 2009; Guillet et al. 2018). The astrodust picture predicts instead that intensity-derived should apply to polarization, and that residual frequency dependence of comes from LOS mixing of temperatures and field orientations, not from composition. For component separation the distinction is operational: a two-population model produces a polarized SED that is not a rescaling of the intensity SED, so using GNILC and maps (derived from intensity) to scale polarized templates is a model error. Until FIR polarimetry between and distinguishes these models (the diagnostic proposed by Hensley & Draine), a B-mode pipeline must marginalize over at least two grain scenarios, not adopt one. Forthcoming wide-field FIR polarimetry (BLAST Observatory, PRIMA in NASA Phase A, FYST 850 GHz) is what decides composition; Taurus, LiteBIRD (to 402 GHz), Simons Observatory (280 GHz), and CMB-S4 (270 GHz) replace Planck 353 GHz as the layer calibrator at higher S/N.
-
Size-dependent temperature of big grains. Smaller equilibrium grains are hotter at fixed radiation field because they absorb optical/UV light more efficiently relative to their FIR opacity. A broad at fixed is already a superposition of MBBs, i.e. a source of intra-cloud temperature moments. McBride et al. (2023) showed that fitting a single MBB to a broad distribution returns biased . Knowing lets one compute the intra-cloud moment structure as a prior, rather than fitting it blindly. Sponseller & Kogut (2022) found that line-of-sight temperature dispersion above about biases for PICO with a single MBB; a 3D temperature map predicts that dispersion per pixel, with the caveats of Zelko et al. (2023).
-
Iron and magnetic dipole emission. Ferromagnetic inclusions produce microwave magnetic dipole emission whose polarization is orthogonal to the electric-dipole thermal emission of the same grains (Draine & Hensley 2013). Hensley & Bull (2018) identified iron grains, together with LOS mixing, as the two most pernicious biases for parametric CMB recovery. A size- and composition-resolved dust model is the only way to predict whether this component is present at a given frequency and whether it rotates ; the prediction is scenario-dependent and must be marginalized.
-
Connecting 3D extinction to 3D emission. The conversion factor depends on the size distribution (optical extinction is dominated by grains; submillimetre opacity has a different size weighting), on porosity, and on . Zhang & Green (2025) give a 3D prior on that conversion. A constrained turns a 3D cube into a 3D cube with a physically bounded systematic error, rather than an unknown scale factor—provided one also specifies composition.
In a pipeline, the grain-size ingredient is therefore not a free SED template. It is a prior mapping
that predicts, for each 3D cell, both the intensity SED and the polarized SED, including their difference, under an adopted grain scenario. The scenario is a discrete systematic, not a solved microphysical field.
3.3 Sky-projected magnetic field
The plane-of-sky field sets the emission EVPA . It does not measure the inclination , and it does not make the polarized sky a derived quantity. Polarized amplitude still requires emissivity, alignment efficiency, inclination, field tangling, and layer association. Without an orientation proxy, even a perfect 3D mass map cannot predict and : one would still have to fit two amplitudes per layer. With an orientation proxy, the number of free polarized parameters per layer drops from two Stokes amplitudes to one (an amplitude, or ), and frequency-dependent becomes a predicted rather than freely fitted effect—if the proxy is correctly associated with the SED of that layer.
Observational tracers of .
-
Submillimetre polarized emission (Planck 353 GHz). This is the most direct full-sky tracer of the dust-weighted, LOS-integrated . Its limitation is exactly the limitation of 2D dust maps: it is an integral. It is the quantity one wants to explain, not the 3D field itself. It remains the amplitude calibrator of any 3D model. In the faintest B-mode fields those maps are noisy at native resolution; GNILC products are used instead, which already average over the small-scale structure that generates small-scale decorrelation.
-
Starlight polarization. Dichroic absorption polarizes background starlight parallel to . Starlight polarization is cumulative. A star “whose Gaia distance falls in the layer” does not directly give that layer’s angle. Differential tomography requires foreground/background stellar pairs or a joint inversion of stellar , with intrinsic stellar polarization, distance uncertainty, selection effects, and angular sparsity propagated (Pelgrims et al. 2023). Pelgrims et al. (2024) inverted RoboPol polarimetry over four square degrees and produced the first degree-scale starlight-polarization tomography map; most sightlines cross more than one cloud, with the plane-of-sky field orientation changing with distance on most of them. PASIPHAE with the WALOP polarimeters is designed to deliver hundreds of polarized stars per square degree at high latitude. Where stars exist, this is the only distance-resolved measurement of in dust. Until PASIPHAE-scale coverage exists, it is a sparse external check, not a full-sky layer angle map.
-
H I filament orientation. Linear structures in narrow H I velocity channels are aligned with (Clark et al. 2014, 2015). The Rolling Hough Transform applied to H I4PI and GALFA-H I yields, with no free parameter relating H I intensity to dust emission, full-sky cubes of Stokes whose velocity integrals correlate with Planck 353 GHz at and reproduce the dust asymmetry at (Clark & Hensley 2019). That result is a position–position–velocity proxy, not a distance-resolved magnetic-field measurement, and it applies to large-scale, velocity-integrated maps. H I coherence across velocity is a predictor of : disordered orientation in – space corresponds to LOS field tangling and depolarization (Clark 2018). H I traces only part of the dusty ISM. Optical-depth corrections, CO-dark molecular gas, ionized-gas dust, IVC dust-to-gas variation, and non-filamentary neutral material can all break the assumed H I–dust layer correspondence. Matching H I velocity components to 3D extinction slices is a cross-identification problem, not a solved mapping. Errors in that matching assign the wrong to the wrong SED and can create a spurious as easily as they correct a real one.
-
Faraday rotation and synchrotron polarization. These trace and the cosmic-ray-weighted field, respectively, i.e. different phases and a different projection. They are complementary constraints on and on the large-scale Galactic field, not substitutes for the dust-weighted . Faraday structures in LOFAR data at tens of parsecs, aligned with H I filaments, have been reported at the Bubble wall (Erceg et al. 2024); a local field reversal plane within a few hundred parsecs has been placed by Booth et al. (2025). These are optional inclination priors, not a map of the neutral dusty medium.
What contributes to dedusting, even without 3D dust or grain sizes.
-
A polarization-angle template that can be applied layer by layer. Given for each cloud, the emission angle is , and the only free polarized parameter per layer is an amplitude (or ). Frequency-dependent is then computed from the Tassis–Pavlidou sum, not fitted as a nuisance—provided the velocity-to-distance association is correct.
-
A data-driven mask or weight for LOS-decorrelated pixels. Pelgrims et al. (2021) demonstrated that an H I-only selection already isolates the pixels in which Planck 353/217 is anomalous. A B-mode analysis can down-weight those pixels in a parametric fit, or give them extra SED parameters, without using the CMB frequencies to discover the problem after the fact. Masking is not automatically an improvement: it reduces foreground bias but also increases sample variance, changes mode coupling, and can create -to- leakage. Improvement must be demonstrated in mean-squared error or -coverage.
-
A geometric model of . Because derived from H I morphology correlates with at over the full sky (Clark & Hensley 2019), a large fraction of the variance in polarization fraction is geometric (LOS tangling and ), not grain-alignment physics. Assigning and a geometric depolarization factor per layer therefore absorbs the leading cause of spatially variable , leaving grain physics to explain the residual, including . The factor remains degenerate with . Two practical options, both used in the literature, are: (i) absorb into a single amplitude fitted to 353 GHz ; (ii) use H I coherence / RHT peak width as a proxy for LOS tangling and , as in Clark (2018). Neither is a precision map.
-
A prior on and parity-odd spectra. Filament–field alignment produces and . A list of angles does not automatically generate that asymmetry; it requires the correct joint statistics of density morphology and field orientation. Cukierman, Clark & Halal (2023) found a global misalignment of H I filaments and Planck dust polarization of , scale-independent above the H I4PI beam, that correlates with the observed dust . A 3D field model that includes a small, spatially coherent misalignment angle can predict, rather than overfit, parity-odd dust spectra that would otherwise leak into cosmic birefringence estimators. Filament population models (Hervías-Caimapo & Huffenberger 2022) can populate sub-beam scales with the correct , , and statistics given a layer’s field orientation.
The sky-projected field is the ingredient that turns a 3D intensity model into a 3D polarization model. Without it, extra layers add extra amplitudes and the fit becomes less, not more, constrained. With it, the fit is still underdetermined at one frequency if , or if two layers are nearly aligned.
4. A combined dedusting strategy
The three ingredients are complementary because they constrain different factors in the polarized radiative transfer. They do not supply the joint distribution by concatenation. A concrete pipeline has a forward-model stage, a parametric (or moment) residual stage, and a validation stage. The association of extinction distance, H I velocity, and starlight polarization must live in one probabilistic likelihood, not in an ex post match. Bandpass integration, thermodynamic-unit conversion, beams, transfer functions, calibration covariance, and reference-frequency colour corrections belong in the generative model; a monochromatic equation is not an likelihood.
4.1 What is actually observed per layer
| Observable | What it measures | What it does not measure |
|---|---|---|
| Gaia / 3D extinction | Cumulative to sampled stars; a regularized differential-extinction posterior | Layer , , , , , polarized amplitude |
| Zelko-like 3D | Voxel temperatures constrained by integrated FIR | Polarized SED; layer-resolved ; uniqueness under the – degeneracy |
| Zhang–Green 3D | Optical/NIR extinction-curve slope | Microwave opacity and alignment efficiency |
| H I4PI / RHT | Brightness and morphology in ; a orientation proxy for some neutral structures | Distance; ionized and molecular dust; ; layer amplitude |
| Starlight | Cumulative dichroic absorption to the star, parallel to | Layer-resolved angle without a differential inversion; emission EVPA without |
| Planck 353 GHz | Dust-weighted, LOS-integrated polarized emission at one frequency | Any layer-resolved quantity; two Stokes observables cannot identify amplitudes |
Layer-resolved , , , , and polarized amplitude are latent. The forward model may predict them under priors; it may not treat them as data.
4.2 Forward model
Discretize each LOS into layers. In practice – is the most that high-latitude data can support: Edenhofer-like distance bins should be clustered into magnetically coherent clouds using H I velocity–orientation structure (as in the Halal et al. 2024 dendrogram), so that one does not fit hundreds of correlated distance slices independently. Treat the Local Bubble shell as the first layer where it dominates, with a geometric prior on orientation from O’Neill et al. (2025), and with still free or absorbed into the amplitude. Treat intermediate-velocity gas as a distinct layer with a distance where stars or a cloud catalogue provide one, and as distance-free otherwise.
Column and distance (3D dust). For layer in pixel ,
where is set by , composition, porosity, and the local , and may be given a Zelko et al. (2023) prior. An overall scale is calibrated on the high-latitude – and – relations, not assumed universal per layer.
Local physical state (grain sizes + environment). Associate to each layer a radiation-field strength and anisotropy from a 3D radiation-field model (not from if is also being predicted), and a density from plus a thickness prior, or from H I brightness temperature in the matching velocity interval. Then, under an adopted grain scenario:
- Compute for the local radiation field by equating absorption of the ISRF to FIR emission, using the adopted dielectric function (astrodust or silicate+carbon; both are retained).
- Compute from RAT/MRAT given .
- Integrate over to obtain layer-averaged , , , and intra-layer spin moments if the size distribution is broad.
This step is the difference between a phenomenological multi-MBB and a physical multi-MBB. It predicts, rather than fits, and within a layer, up to the discrete grain-scenario systematic.
Orientation (sky-projected ). Assign from a joint likelihood of (i) the H I-derived angle in the velocity range matched to layer and (ii) differential starlight tomography where stars exist. The emission angle is . The association of H I velocity with extinction distance is a latent variable with a prior, not a delta function. False cloud splits and mergers, H I optical depth, missing gas phases, and starlight selection are part of that likelihood.
Emission integral, band-integrated. In IAU coordinates, at each experiment band with bandpass and in thermodynamic CMB units,
with an analogous sum for , intra-layer spin moments inside each layer, and beams applied in harmonic space. If the aligned population is a subset of the emitting population, and the polarized amplitude does not scale with . That is a prediction of the grain scenario, not an extra free map.
This forward model generates, when the ingredients are correctly associated:
- frequency-dependent on multi-cloud, misaligned, comparable-amplitude sightlines (Tassis & Pavlidou);
- with a specific frequency dependence set by the layer , , and angle–SED covariance, not by a free —or if layer SEDs are the same or morphologies remain proportional;
- a difference between intensity-derived and polarization-derived effective spectral indices, if the grain scenario produces one;
- geometric asymmetry from filament–field alignment, if the joint density–orientation statistics are right.
It does not generate those features “automatically” from a list of columns and angles.
Identifiability. At 353 GHz there are two Stokes observables per pixel. For , a single polarized MBB with fixed is identified, up to the usual – degeneracy. For with known, distinct angles, two amplitudes remain; one frequency is then formally sufficient at the reference band, and additional frequencies constrain SED differences. The mixing matrix becomes ill-conditioned when or when the layer SEDs are nearly identical, and CMB remain degenerate with a dust residual that is spectrally close to a blackbody. For , 353 GHz cannot identify the amplitudes even with known angles; the extra frequencies of a LiteBIRD- or CMB-S4-like experiment help only to the extent that layer SEDs differ. Nearly aligned layers and nearly identical SEDs are the generic high-latitude case, not the exception. The 3D ingredients therefore work as priors that reduce the volume of the degenerate subspace, not as a unique inversion.
4.3 Parametric residual fit
The forward model will not be perfect. The residual is fitted to the multi-frequency CMB data with a restricted parameter set, using the 3D model as a strongly informative prior rather than as a rigid template. A workable hierarchy, from fewest to most free parameters:
H0. Template cleaning with a 3D-informed mask. Use Planck 353 GHz as the dust template, scale with a global or piecewise , but exclude or deweight pixels that the 3D+H I selection flags as LOS-decorrelated (Pelgrims target sample). This is the minimal use of the three ingredients. Tassis & Pavlidou proposed the starlight-polarization version of this mask in 2015. It is a hypothesis to be tested in -coverage, not already an improvement.
H1. Layered MBB with amplitudes fitted at a reference frequency. Fix from the orientation likelihood. Put Gaussian priors on from the grain+3D model, with a width that reflects grain-scenario uncertainty rather than a delta function. Fit (and optionally a small ) to the ensemble of frequency maps, with on most high-latitude pixels and – only where the 3D map and H I demand it. The CMB is a blackbody added in . Degeneracies are controlled because is not free, but they are not eliminated.
H2. Layered MBB plus intra-layer spin moments. Allow complex spin-2 moments and per layer, or, cheaper, a single set of spin moments on the sum of layers to capture beam averaging and unresolved substructure (Chluba et al. 2017; Mangilli et al. 2021; Vacher et al. 2022, 2023a). At LiteBIRD sensitivity the Vacher et al. (2022) result implies that moments in both and are required. The 3D model predicts the expected size and the real/imaginary structure of those moments, which can be imposed as a Gaussian prior, preventing the moments from absorbing CMB. If the layer angles are taken as fixed and correct, the remaining complex moments of the sum should be small; that is a testable claim, not an assumption.
H3. Harmonic-space hybrid. For likelihoods of the BICEP/Keck type, replace the single-SED, single- dust model by the cross-spectra of the layered forward model, plus a low-order spin-moment expansion of the residual. The layered model supplies the frequency dependence of and of ; the moments absorb isotropic SED distortions from averaging. This is the natural interface to LiteBIRD and CMB-S4 power-spectrum pipelines, and to constrained ILC that deprojects the predicted moment maps. Simulation-based inference trained on 3D-informed simulations is a consumer of the same product, not a substitute for the physics.
Synchrotron is treated in parallel: a power law with spatially varying , optionally a curvature term, and an independent (the synchrotron-emitting volume is not the dust-emitting volume). H I-based should not be imposed on synchrotron. Neutral-gas complexity can inform synchrotron depolarization, but not its orientation. Joint dust–synchrotron modelling is needed below because a dust residual with a wrong is degenerate with synchrotron, and conversely.
4.4 Why the combination is more than the sum of parts
Each ingredient removes a specific degeneracy, incompletely:
| Degeneracy in a 2D MBB fit | Constrained by | Residual freedom |
|---|---|---|
| One vs. a sum of | 3D layers | Layer SED still latent; conversion |
| Intensity vs. polarized | grain-size / alignment scenario | Composition, optical constants; astrodust vs two-population |
| Frequency-independent vs. rotating | per-layer if correctly matched to layers | Velocity distance; starlight sparsity |
| Geometric vs. grain | tangling (H I coherence) vs. RAT | unmeasured; remains one number |
| , vs. CMB | 3D-predicted spin-moment prior | Unresolved sub-beam turbulence |
| Which pixels need extra SED freedom | 3D cloud count + misalignment angle | Masking-induced leakage |
A 3D mass map without overfits . A map without 3D mass cannot assign SED differences to layers. Grain physics without either map can predict a mean but not its spatial pattern. Together they specify the Tassis–Pavlidou two-cloud (and -cloud) model on the real sky as a prior. Liu et al. (2025) found that at high complexity the bias on for Simons Observatory is driven by spatial correlations between the dust index and the dust amplitude at higher orders, not by the amplitude structure itself. Those correlations are exactly what a layered model with per-layer physics encodes, and what a random-field moment model treats as nuisance variance.
4.5 Validation, not only subtraction
A physical forward model is falsifiable, as is a phenomenological MBB: both can be tested with held-out frequencies, residual spectra, and posterior-predictive checks. A high-dimensional physical model can be harder to falsify if flexible amplitudes and moments absorb discrepancies. The tests that matter are held-out and end-to-end.
Required tests:
- Fit only ancillary data plus 353 GHz, freeze the model, and predict 217/220 and then 95/150 GHz , , , , and all cross-frequency . Perform this separately on Pelgrims target/control samples and on the BICEP/Keck and SPIDER fields, where the current answer is “no sky-average decorrelation”. The model must predict that too, not only the Pelgrims excess.
- On the Pelgrims control sample (single H I cloud), the layered model must collapse to a single MBB and must be constant; any remaining 353/217 is then a bound on residual systematics and on intra-cloud grain physics.
- On the target sample, the model must predict the sign and the approximate magnitude of from the H I-measured and the 3D-inferred , before looking at 217 GHz. Pelgrims et al. estimated a 15% polarized-intensity-ratio variation; a grain+3D model should predict a number of that order from and differences, not fit it.
- The predicted from H I coherence plus RAT should match Planck with residuals that correlate with tracers of grain alignment (e.g. , radiation field), not with the CMB. Halal et al. (2024) is the existing version of this test for 3D complexity versus .
- Cross-spectra of the model, evaluated on the actual experiment bandpasses, should be compared to Planck 217353 and to BICEP/Keck 95150220Planck 353 without introducing a free . A free is reserved for the residual. Residual covariances should be compared to Guillet et al. (2026).
- Spectral-imaging tests with bolometric interferometry (Regnier et al. 2024, A&A 686, A271—not a Maffei et al. paper) detect the systematic from unaccounted decorrelation in a CMB-S4-like setup. They are a detector, not a corrector, and pair naturally with a predictive model: a successful 3D model should make the detection vanish.
- Ablation tests: 3D extinction, H I orientation, starlight tomography, and grain priors added separately. A strong 353 GHz calibration may dominate the result and conceal that the purported physical inputs add little predictive information.
- Blinded end-to-end simulations in which the reconstruction receives mock Gaia stars, H I spectra, and starlight-polarization catalogues—not the true 3D cube. Include MHD substructure, a Bubble-like cavity, multiple grain compositions, out-of-model dust, synchrotron, lensing, realistic bandpasses, beams, calibration/leakage systematics, and survey masks. Compare against one-MBB fits, independent SEDs, scalar and spin moments, standard and constrained ILC, and phenomenological decorrelation likelihoods. Report , , residual , coverage, sky fraction, and failure rate—not only goodness of fit.
Until those exist, the construction is a well-posed physical model for a demonstrated mechanism, not a turn-key replacement for moment expansions and multi-frequency residual fits.
4.6 An -budget, as a scaling, not as a result
After any linear or parametric clean, the residual dust B-mode power at a CMB band is of order the dust power there times the fraction of that power the model does not describe. The residual is positive-definite, so it biases upward rather than adding noise. Frequency-dependent rotation leaks dust into .
Order of magnitude in a BICEP/Keck-like deep field, using a modified blackbody with at and a tensor spectrum of order at for : dust at 150 GHz is worth in -equivalent, and at 95 GHz . Each percent of unmodelled dust power at 150 GHz is therefore worth about in in the cleanest deep field. LiteBIRD’s target of order and a CMB-S4 statistical error of a few require the unmodelled fraction to stay well below a percent. That is the scale of a degree of coherent unmodelled rotation, or of an unmodelled effective index dispersion of several hundredths—numbers that Planck XI does not rule out on small patches, that Osumi et al. (2021) leave open at high latitude, and that Pelgrims et al. (2021) detect statistically in a selected few-percent of high-latitude pixels.
If a fraction of the polarized power in a deep field comes from Pelgrims-like sightlines, the bias is times the two-cloud residual. For between and that is a few to about at 150 GHz—if those sightlines have comparable polarized amplitudes and if the cleaning does not already absorb the effect into a free or . This is exactly the next-generation error budget, and it is a bias. It is not a measurement of the residual of any existing 3D pipeline.
The honest bottom line: the approach reduces bias by the fraction of line-of-sight-mixed power it correctly identifies, reduces error inflation by replacing free moments with constrained ones, and supplies null tests at 353 and 217 GHz that no CMB-only method offers. It does not remove the need for spin moments, and its gain is bounded by the accuracy of the layer assignment.
5. Limitations, missing data, and residual systematics
A layered physical model reduces, but does not eliminate, foreground bias on . The leading limitations, several of them blockers for a claim of controlled residuals, are as follows.
The core inverse problem remains underdetermined. The proposed inputs do not identify layer-resolved polarized amplitudes or SEDs. Discretizing the LOS is not the same as making the inversion identifiable. At one frequency, more than two layer amplitudes are underdetermined; similar MBB SEDs preserve the degeneracy across bands; nearly aligned angles make even two amplitudes ill-conditioned.
No end-to-end demonstration connects the ancillary products to controlled residuals or unbiased . Masking, spin-2 residuals, synchrotron covariance, bandpasses, beams, instrumental systematics, lensing, and ancillary-map uncertainty intervene before an likelihood is reached. Present “high-complexity” PySM suites (layered MBBs, d12) capture LOS decorrelation phenomenologically but do not ingest Gaia 3D maps or RAT-based alignment as observed products.
The H I-velocity-to-extinction-distance association is assumed more often than it is inferred. Velocity is not distance at high latitude. Galactic rotation does not provide a unique kinematic distance; two clouds at the same velocity can sit at different distances, and a single cloud can be split across velocities by turbulence. A wrong match assigns an angle to the wrong SED. Joint reconstructions of H I and 3D dust exist at Galactic-disk scales; they are not a solved high-latitude layer catalogue.
3D maps trace extinction, not submillimetre polarized emission. The conversion assumes a dust-to-extinction ratio and a size distribution that may vary with environment ( variations, grain growth in molecular gas, carbon depletion). High-latitude sightlines are mostly diffuse, where Lenz et al. (2017) find scatter in , but 10% layer-to-layer opacity errors, if correlated with or , are enough to leave a residual at the level. Zelko et al. (2023) already find a factor-of-two variation in the emission-to-reddening conversion. CIB anisotropies further contaminate intensity-based and maps at high latitude.
Distance completeness. Gaia-based maps reach at high resolution (Edenhofer et al. 2024). That covers most high-latitude dust, but IVCs/HVCs can contribute polarized emission at high and lie outside the local volume or are poorly sampled by stars. Pelgrims et al. used IVC/LVC H I components precisely because they are common at high latitude. A pipeline that trusts only Gaia distances will mis-model those layers unless H I velocity components are included as additional, distance-free clouds.
The inclination is poorly measured. Polarized emission constrains , not . The factor is degenerate with and with grain alignment efficiency. H I filament detectability encodes some information (filaments end-on are less linear on the sky), and starlight encodes some, but neither is a precision map. Wrong mis-weight the layers in and therefore bias the predicted . The tangential-field assumption for the Local Bubble shell is not confirmed by polarization-fraction statistics (Halal et al. 2024).
Grain model uncertainty is a systematic on and . Silicate+carbon with unaligned carbons versus astrodust versus models with magnetic dipole emission make qualitatively different predictions for the polarized SED below . Until FIR polarimetry distinguishes these models, a B-mode pipeline must marginalize over at least two grain scenarios.
Alignment physics is not a solved 3D field. RAT predictions for are tested statistically, not cloud-by-cloud at high latitude. Dense, shielded gas can have inefficient alignment; locally enhanced radiation fields can raise and together. Mis-modelling alignment produces a spatially variable that the pipeline may try to absorb into or into extra layers.
Unresolved sub-beam structure still generates spin moments. Even a perfect cloud catalogue does not resolve MHD turbulence inside a CMB beam (LiteBIRD’s 140 GHz beam is ; CMB-S4 will be finer but not arbitrarily fine). Intra-beam tangling and variations remain, and they are exactly what the Chluba/Vacher moment expansion is for. The 3D model reduces the moments; it does not set them to zero. Vacher et al. (2022) showed that at LiteBIRD noise, omitting moments biases even when moments are included. That conclusion survives a 3D prior.
Planck 353 GHz polarization S/N is limited at high latitude. All 3D polarized models are calibrated on Planck 353 GHz . Extrapolating a smoothed 353 GHz template cannot recover small-scale . Higher-resolution polarized dust data are needed to calibrate layers at CMB-S4 resolution. Circular calibration and validation on the same 353 GHz maps must be avoided: held-out 217 GHz is the minimum.
The nonlocal nature of . Pixel masking and spatially varying weights create -to- mixing; angle and amplitude errors contribute differently to ; reproducing an H I–Planck map correlation at large scales does not establish control of the degree-scale dust modes used for .
Instrumental and analysis systematics mimic LOS decorrelation. Bandpass mismatch, leakage, gain mismatch, and residual CMB in the 353 GHz template all produce frequency-dependent residuals. Pelgrims et al. found that correcting Planck residual systematics increased the significance of LOS decorrelation, which is encouraging (the astrophysical effect is not an artefact), but the converse is also possible in other data sets. A pipeline that applies a 3D correction derived from Planck to a ground-based deep field must re-estimate the correction on that field’s own ancillary data.
Other foregrounds. Polarized synchrotron with spatially varying and possible curvature; a possible polarized AME at the percent level; CO and other line polarization at a poorly known level; polarized compact sources. A perfect dust model does not yield an unbiased if synchrotron is wrong. The 3D dust/B-field construction does not automatically help synchrotron, which requires low-frequency surveys (C-BASS, S-PASS, QUIJOTE, Planck LFI, future 5–20 GHz polarimetry).
Statistical cost. Extra layers and moments consume degrees of freedom. Vacher et al. (2022) already saw increase because dust moments correlate with the CMB SED. A 3D model that adds unconstrained amplitudes per layer will degrade more than it reduces bias. The whole point of the three ingredients is to replace free maps by constrained ones: from H I/stars, from Gaia, and from grain physics with tight priors. If those constraints are too tight and wrong, one reintroduces bias; if they are too loose, one reintroduces degeneracy. That bias–variance trade is the central design problem.
What is missing as data, not as theory.
- Dense, high-latitude optical polarimetry of stars with Gaia distances (PASIPHAE-like), sufficient to measure between the two or three local dust layers in every B-mode pixel.
- A joint 3D reconstruction that uses Gaia extinction, H I cubes, and starlight polarization in one likelihood, rather than matching products ex post.
- FIR polarimetry from to to decide between grain models and to measure in polarization, where CIB contamination is weaker (Hensley & Draine).
- Direct 3D polarized emission maps, which do not exist; extinction and Zelko-like 3D temperature are proxies for the intensity side only.
- A well-characterized tracer in the neutral medium.
- End-to-end simulations that start from an MHD+dust cube with a realistic and RAT alignment, observe it with LiteBIRD/CMB-S4 bandpasses and with realistic ancillary-data incompleteness, and show that the pipeline in §4 recovers without a bias larger than the statistical error.
Until those exist, the honest status is: the three ingredients specify a well-posed physical model for the exact mechanism (Tassis & Pavlidou 2015; Pelgrims et al. 2021) that makes a single-MBB-per-pixel subtraction fail; they can already be used to flag and down-weight the worst sightlines and to predict the frequency dependence of and at a level that phenomenological decorrelation parameters cannot, provided the layer association is validated on held-out frequencies; they are not yet a turn-key replacement for spin-moment expansions and multi-frequency residual fits. The right operational posture for LiteBIRD and CMB-S4 is to implement the layered forward model as an informative prior and as a diagnostic, to keep a low-order complex – spin-moment expansion for unresolved averaging, and to treat grain-composition and uncertainties as discrete systematics that are marginalized, not assumed away.
6. Summary
Spectral-index decorrelation of polarized dust is the statement that the effective MBB parameters and, when SED and orientation covary, the polarization angle of Galactic dust are frequency-dependent, because Stokes are integrals of emitters that do not share a common . A single modified blackbody per pixel—or a 353 GHz template scaled by a global —is then a misspecified model. Recovered CMB are strongly biased only when SED and angle vary together (McBride, Bull & Hensley 2023); that is precisely the regime the three ingredients target. Sky-average Planck and BICEP/Keck analyses do not detect a loss of correlation at the level relevant to ; a selected multi-cloud, misaligned sample does (Pelgrims et al. 2021); next-generation experiments need control below the percent residual that current CMB data cannot guarantee.
Three-dimensional dust maps (presently Gaia extinction cubes, Zelko-like 3D temperatures, and 3D , matched probabilistically to H I clouds, with the Local Bubble as a candidate first layer) replace the integral by a sum over a few distance layers and identify the sightlines where the Tassis–Pavlidou mechanism can operate. The grain size distribution and RAT alignment determine which grains emit in polarization, why need not equal and , and how to convert extinction into a polarized SED given the local radiation field and density—under an adopted grain scenario that must be marginalized. The sky-projected magnetic field, from H I filaments, starlight tomography, and Planck 353 GHz, assigns (and a geometric proxy for depolarization, not ) to each layer, so that frequency-dependent is predicted rather than fitted, if velocity is correctly matched to distance.
The combined pipeline is a layered, physically parameterized, band-integrated forward model whose remaining freedom is a small set of amplitudes and intra-layer spin moments, fitted to the CMB frequency maps under priors taken from the ancillary 3D data. Its limitations are real and some are blockers: extinction is not emission, velocity is not distance, and grain composition are underconstrained, the inverse problem remains underdetermined, beam averaging still requires complex moments, and no end-to-end demonstration yet exists. Used as a prior and a diagnostic rather than as an exact sky, the same construction is currently the most direct way to attack the specific astrophysical failure mode that spectral-index decorrelation names.
What I took from the critic
The controlling physics is SED–angle covariance, not “any multi-layer beam”. I rewrote the opening mechanism, the two-cloud discussion, and the McBride citation around that point, including the comparable-amplitude condition and Poh & Dodelson’s result that temperature-only extrapolation noise is small compared with then-current errors.
I adopted the critic’s convention demands: IAU versus COSMO sign, emission EVPA versus (the rotation), and the statement that does not measure and does not make the polarized sky a derived quantity.
I corrected the moment expansion (amplitude convention; spin-2 complex moments rather than scalar H2 moments), the ILC description, the – formula’s assumptions, and the conversion (the factor and the non-universal ). I removed the unusable RAT formula and the unattributed “15% is a 10% SED-shape change” sentence.
I restated the observational status as the critic demanded: Planck 2018 XI and Sheehy & Slosar walked back Planck Int. L; BK18 and BK XVI find no decorrelation in the deep field; SPIDER finds spatial SED variation without a full-region LOS detection; Pelgrims is a selected-sample detection with a stated -value range.
I accepted the identifiability critique as a blocker, not a quibble: the three products do not supply the joint distribution; layer association must be probabilistic; H0 masking is not already an improvement; a physical model is not uniquely falsifiable relative to an MBB. The pipeline section now states what is data, what is prior, and what is fitted, and it lists held-out prediction, ablation, and end-to-end simulations as requirements rather than as accomplished facts.
What I took from the optimist
I kept the draft’s useful skeleton (spatial versus LOS decorrelation, two-cloud algebra, alignment as a size high-pass filter, the degeneracy table, the H0–H3 hierarchy, held-out Pelgrims tests) and upgraded it with the optimist’s strongest verified observational handles.
The Local Bubble as a candidate dominant first layer (Skalidis & Pelgrims 2019; O’Neill et al. 2025), with Halal, Clark & Tahani (2024) as the necessary caution that the tangential-field assumption is not seen in and that complexity beyond the shell matters. Zelko et al. (2023) as the existing prototype of a 3D emission/temperature map. Zhang & Green (2025) as a 3D grain-size proxy via . Pelgrims et al. (2024) as the first degree-scale starlight tomography. BICEP/Keck XVI as evidence that H I templates carry some -mode phase information down to 95 GHz, with the 10–65% correlation quoted honestly. Planck Int. XVII as per-component SED evidence, including the extra component.
I adopted spin moments as the interface between the forward model and the likelihood (Vacher et al. 2023a,b, 2025), including the angle-rotation relation, and Guillet et al. (2026) as a residual-covariance diagnostic, with the SRoll2-versus-PR4 split preserved. Liu et al. (2025) and Osumi et al. (2021) supply the argument that high-latitude variance is poorly bounded and that –amplitude correlations, not just amplitude structure, drive bias. I corrected the spectral-imaging citation to Regnier et al. (2024) and tempered the offset to Planck XI’s . I included an -budget as a scaling, and the FIR/high-frequency roadmap as the data that decide grain composition versus layer calibration.
What I rejected, and why
I rejected the original draft’s claim that multiple layers automatically make , , and frequency-dependent, and the claim that with “the polarized sky becomes a derived quantity”. Both are physically false: rotation needs SED–angle covariance and comparable amplitudes, and a projected orientation is not an amplitude or an inclination.
I rejected presenting extreme two-cloud rotations of “tens of degrees” as generic, and I rejected converting Pelgrims’s 15% band-ratio spread into a 10% SED-shape change. I rejected the printed RAT formula, the mixed moment kernels, the statement that standard ILC nulls a single dust SED, and the assertion that a physical model is falsifiable while an MBB is not.
I rejected several of the optimist’s upgrades that I could not verify as published results, rather than folding them in on trust: a 2026 velocity-resolved 3D H I map attributed to McCallum, a 2026 O’Neill catalogue of dust clouds with IVCs in distance, and a set of 2026 alignment, Faraday, and LiteBIRD-index papers that did not resolve as citable. The scientific need those citations were meant to fill (velocity–distance matching, IVCs, alignment tests) is already stated as a limitation.
I rejected the optimist’s implication that the high-latitude sum is already a finite, catalogued, solved object whose phases are predicted rather than fitted. The Local Bubble is a strong first-layer prior; Halal et al. (2024) show that dust beyond the shell still sets . I rejected treating the -budget table as a measured residual of a 3D pipeline: it is a scaling that shows why sub-percent control matters, not a demonstration that the control exists. I rejected temperature-dominance as a settled sky property: Guillet et al. (2026) find it for SRoll2 residual covariances and not for PR4.
I rejected “H0 is already an improvement” and “3D extinction provides the same information in true distance as an H I cloud count”. Masking must be scored on -coverage; extinction and H I are complementary, not equivalent. I also rejected the idea that the layered model replaces moments. Unresolved sub-beam structure still needs spin moments; the 3D model is a prior on those moments, not a substitute.