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

Dedusting CMB B-modes with distance-resolved dust, grain physics, and the projected magnetic field

The argument of this essay is one causal chain. A distance-resolved dust column tells where the emitters are along each line of sight. Grain physics tells how each emitter scales with frequency and how strongly it polarizes. Distance-resolved knowledge of the plane-of-sky magnetic field tells the direction in which each layer contributes to the complex polarization $Q+iU$. The vector sum of those layer contributions predicts a chromatic dust foreground, one whose polarization pattern changes shape with frequency, and that is precisely what a single-template extrapolation from 353 GHz cannot represent.

1. Polarized thermal dust and why B-mode studies care

Primordial gravitational waves imprint a curl-like B-mode pattern on the CMB polarization with an amplitude set by the tensor-to-scalar ratio $r$. The signal peaks at the recombination bump near multipole $\ell\approx 80$, an angular scale of about two degrees, and at the reionization bump at $\ell<10$. Scalar perturbations produce no B-modes at linear order, so a detection measures the energy scale of inflation directly. Current bounds are $r<0.036$ at 95% confidence from BICEP/Keck data through 2018 combined with Planck and WMAP, and $r<0.032$ when Planck large-scale polarization is added by Tristram et al. 2022. The next generation, Simons Observatory, LiteBIRD and CMB-S4, targets $\sigma(r)$ of order $10^{-3}$.

Two nuisance signals compete with that target. The first is the lensing B-mode, produced when large-scale structure lenses E-modes into B-modes. It is a guaranteed cosmological signal with a nearly white spectrum equivalent to roughly 5 μK-arcmin of noise, and near $\ell\approx 80$ it exceeds a primordial signal with $r\approx 10^{-2}$. It is removed by delensing, not by foreground cleaning. The second is Galactic emission: synchrotron dominating below roughly 70 GHz and thermal dust dominating above roughly 100 GHz.

Thermal dust is polarized because non-spherical grains spin with their short axes preferentially aligned with the local magnetic field. They emit preferentially with the electric vector along their long axes, which is perpendicular to the plane-of-sky field. Planck measured polarization fractions up to about 20% at 353 GHz, an E/B power ratio near 2, a positive TE correlation, and a mean polarized spectral index of 1.5 to 1.6 with a temperature near 19.6 K in the diffuse high-latitude interstellar medium. Planck Intermediate Results XXX also established that there is no dust-free window: even in the cleanest few percent of the sky the dust B-mode power extrapolated to 150 GHz is comparable to a primordial signal with $r$ of order 0.01 to 0.1 at $\ell\approx 80$, and in the BICEP/Keck field it corresponds to $r$ of order 0.1. Reaching $r\approx 10^{-3}$ therefore requires modelling the dust polarization to better than about 1% in power, which is a few percent in amplitude, at the CMB channels.

The standard approach does this by extrapolation. One assumes that the dust $Q$ and $U$ maps have the same shape at every frequency and differ only by a scalar factor $f(\nu)$ given by a modified blackbody with a single $\beta$ and $T$ per pixel or region. The Planck 353 GHz map, or a dust component fitted jointly across frequencies, then serves as one template. The contrast in this essay is with a multi-layer physical model, in which each pixel contains several emitters with different SEDs and different polarization angles, so that the sky pattern itself becomes frequency dependent.

2. What "spectral-index decorrelation" means operationally

2.1 Correcting the phrase, and the operational definition

The spectral index does not decorrelate. What decorrelates is the dust polarization map across frequencies, and the cause is spatial and line-of-sight variation of the SED parameters, combined in the line-of-sight case with variation of the polarization angle. The operational statistic is the correlation ratio of B-mode cross-spectra,

with the auto-spectra estimated from cross-spectra of independent data splits so that noise does not bias them. If the two maps are proportional, $\mathcal R_\ell=1$. If $\mathcal R_\ell<1$, a fraction $1-\mathcal R_\ell^2$ of the dust B-mode power at $\nu_1$ cannot be predicted by any rescaling of the $\nu_2$ map. After template subtraction, or after a fit with a single SED, that fraction survives as a residual that looks like B-mode power. In a BICEP/Keck-like field at 150 GHz a 1% residual in dust power is a bias in $r$ of order $10^{-3}$, so sub-percent decorrelation already matters for next-generation targets.

Equally important operationally: a value below unity is not automatically astrophysical. Noise bias in the auto-spectra, bandpass mismatch between detectors, polarization-efficiency and gain errors, beam mismatch, and correlated noise between splits all produce or mimic decorrelation, and the statistical significance depends on the correlation between overlapping sky regions. These effects explain why the early Planck claim did not survive, as discussed in section 2.5.

2.2 The minimal model: polarization is a spin-2 vector sum

Write the normalized SED of emitter $k$ as a modified blackbody,

In a real pipeline $f_k$ is integrated through each bandpass $\tau_\nu(\nu')$, converted from intensity units to thermodynamic temperature, and color-corrected, so the coefficient that multiplies a layer's map in band $\nu$ is

The polarized emission from a sightline made of $K$ layers is

where $A_k$ is the layer's intensity at $\nu_0$, $p_k$ its intrinsic polarization fraction set by grain shape and alignment, $\gamma_k$ the inclination of its field to the plane of the sky, and $\psi_{E,k}$ the polarization angle of its emission, which equals the field angle plus 90 degrees. The factor $e^{2i\psi}$ is the essential point. Stokes $Q$ and $U$ form a spin-2 quantity, a headless vector in the $(Q,U)$ plane that rotates by $2\psi$ when the polarization angle rotates by $\psi$. Adding two polarized signals means adding two such vectors, not two scalar intensities. Two layers with equal polarized amplitude and angles differing by 90 degrees cancel exactly; angles differing by 45 degrees add in quadrature. Total intensity has no such cancellation, which is one reason the polarized SED of a sightline need not equal its intensity SED.

2.3 Two clouds: SED contrast times angle contrast

Take two layers with polarized amplitudes $P_1$, $P_2$ at $\nu_0$, angles $\psi_1$, $\psi_2$, and SEDs $f_1$, $f_2$. Define the frequency-dependent amplitude ratio $\rho(\nu)=P_2f_2(\nu)/[P_1f_1(\nu)]$ and the misalignment $\Delta=\psi_2-\psi_1$. Then

The observed angle is $\psi(\nu)=\psi_1+\tfrac12\arg\!\left[1+\rho(\nu)e^{2i\Delta}\right]$, and its frequency derivative is

with $x_k=h\nu/kT_k$. The angle rotates with frequency only if both factors are nonzero. The factor $\sin 2\Delta$ requires the layers' fields to be misaligned, neither parallel nor perpendicular. The factor $d\ln\rho/d\ln\nu$ requires their SEDs to differ, through $\beta$, through temperature, or both. Misaligned fields with identical SEDs give a vector sum that scales uniformly, so no rotation. Different SEDs with parallel fields give a scalar sum whose amplitude follows a non-modified-blackbody effective SED, but the angle is fixed. This is the two-cloud mechanism of Tassis & Pavlidou 2015, who showed that large decorrelation between 350 and 150 GHz needs both a frequency-dependent ratio of cloud polarized intensities and significantly misaligned fields. Poh & Dodelson 2017 named the statistical consequence line-of-sight extrapolation noise. With Planck-era uncertainties they found it more than an order of magnitude below the errors, but potentially important once experiments become sensitive to $r<0.0015$.

A worked example: equal polarized amplitudes at 353 GHz, $\Delta=45^\circ$, and an effective SED contrast $d\ln\rho/d\ln\nu=0.2$, which could be a $\beta$ difference of 0.2 or a smaller $\beta$ difference combined with a temperature difference.

The angle rotates by about 2.5 degrees between the two frequencies. A one-template extrapolation cannot reproduce this. The best scalar rescaling of the 353 GHz vector leaves a residual equal to $\sin(2\,\delta\psi)$, about 8.5% of the polarized amplitude and about 0.7% of the dust power, which is an $r$ bias of order $10^{-3}$ in a BICEP/Keck-like field.

With $\Delta=90^\circ$ there is no rotation, but the cancellation between layers is frequency dependent. For $\rho(353)=0.8$ the net amplitude is 0.2 at 353 GHz and 0.33 at 150 GHz, in units of layer 1. The net polarized amplitude therefore falls from 353 to 150 GHz by a factor 1.6 less than layer 1 alone, equivalent to an effective $\beta$ lower by about 0.6 than either layer, while the same pixel's intensity SED is an ordinary average. Cancellation geometry amplifies SED contrast, and it is invisible in intensity.

2.4 Spatial variation lowers the correlation without any rotation

Even if every sightline contained a single layer, the SED parameters vary from pixel to pixel, and maps at two frequencies are then not proportional. For a Gaussian $\beta$ field with dispersion $\sigma_\beta$ that is statistically independent of the polarized amplitude field, the map-level correlation is

Representative values for $\sigma_\beta=0.1$:

The quadratic dependence on the logarithmic frequency separation means that 217 versus 353 GHz is a weak probe of what happens at 150 or 95 GHz. Two features distinguish this spatial mechanism from the line-of-sight one. No polarization angle rotates within a pixel, and the effect is fully captured by letting $\beta$ vary spatially, which is what moment expansions do. Line-of-sight decorrelation additionally rotates angles, makes polarized SEDs differ from intensity SEDs, and makes the apparent single-layer $\beta$ vary across the sky with the layers' mixing ratio. In real data the two mechanisms are therefore coupled, and the Simons Observatory study by Liu et al. 2025 found that the bias on $r$ at high complexity depends on how $\beta$ and the dust amplitude are spatially related, not on the amplitude structure alone.

2.5 What has actually been measured

• Planck Intermediate Results L (2017), using 2015 maps, reported a 217×353 GHz B-mode correlation ratio below unity, marginal per region but above 99% confidence when regions were combined, with decorrelation increasing as the mean column density of the region decreased.
• Sheehy & Slosar 2018 reanalysed the same data and found that the no-decorrelation hypothesis could not be excluded, once noise bias in the split cross-spectra, evidence for instrumental systematics, and the correlations between overlapping regions were accounted for.
• Planck 2018 XI, with the 2018 maps, found no evidence for decorrelation and noted that if its limit holds in the smaller regions targeted by sub-orbital experiments, decorrelation should not prevent a detection at $r\approx 0.01$.
• BICEP/Keck model decorrelation with a single parameter $\Delta_d$, the 217×353 GHz correlation at $\ell=80$, with an assumed scaling in frequency proportional to $\ln^2(\nu_1/\nu_2)$ and a linear scaling in $\ell$. Their data constrain $\Delta_d$ to be above 0.98 at 68% confidence, and BICEP/Keck XVI found no decorrelation between H I filaments and the rest of the dust field in their region, nor from dust associated with intermediate-velocity H I.
• Pelgrims et al. 2021 performed the first targeted test of the line-of-sight mechanism. Sightlines were sorted with H I cloud counts and orientation templates into a target sample with multiple misaligned clouds and a control sample, and the polarization-angle difference between 217 and 353 GHz was compared. The probability that the effect is noise ranged from $5\times10^{-2}$ to $4\times10^{-7}$ depending on the CMB subtraction, and the effect grows with the misalignment.

The lesson is that Planck has not established a universal decorrelation coefficient. The two Planck analyses used different data releases, splits and estimators, and a single number for 217×353 GHz at $\ell=80$ would not transfer to 95 or 150 GHz or to a particular deep field. Decorrelation is a sightline-dependent quantity, and a distance-resolved model is the only kind that can predict it.

3. How the three ingredients help

3.1 3D dust maps: where the emitters are

The 3D dust maps available today are reconstructions of differential extinction or dust density, built from stellar photometry or spectra together with Gaia parallaxes: the Bayestar maps from Pan-STARRS, 2MASS and Gaia, the Lallement and Vergely maps, the Leike and Enßlin reconstruction, and Edenhofer et al. 2024, which reaches 1.25 kpc with an angular resolution of 14 arcminutes, parsec-scale radial resolution near the Sun, 54 million stars with Gaia BP/RP-based extinctions, and 12 posterior samples for uncertainty propagation. These are not tomographic measurements of microwave emission. Turning them into emission requires a column-to-emission conversion per layer, an emissivity per unit reddening and a temperature, and both depend on environment at the 10 to 30% level. A 3D emissivity field obtained by combining a density map with a radiative-transfer model of grain heating is a model product and should be labelled as such.

Even as extinction maps they answer the question the one-template method cannot ask: how many distinct dust structures lie along each sightline, at what distance, and with what column. At Galactic latitudes above 30 degrees most of the dust column sits within a few hundred parsecs, dominated by the wall of the Local Bubble and a handful of intermediate-latitude clouds, so a 1.25 kpc map contains most of the emitting dust in the regions where B-mode experiments observe. Intermediate-velocity clouds at kiloparsec heights contribute in places, and Pelgrims et al. 2024 found starlight-polarizing clouds out to 2 to 3 kpc that appear to correspond to them. The map posterior gives, for each pixel, the layer columns $E_k(\hat n)$ with covariance, and the cross-pixel coherence of a layer, the same cloud appearing in adjacent pixels with its own morphology, is what allows its emission to be regressed out of far-infrared intensity maps.

H I velocity structure supplements this. Panopoulou & Lenz 2020 decomposed HI4PI spectra into Gaussian components and counted on average about three velocity components per 0.2 square-degree pixel in the North Galactic cap and 2.5 in the South. Clark & Hensley 2019 built velocity-resolved orientation templates. But velocity is not distance. H I channel maps separate components kinematically, they miss molecular and ionized gas, and they should not be called a geometric 3D map.

3.2 Grain size distribution: how each emitter scales with frequency and polarizes

Schematically, the polarized emission along a sightline is

Here $n_c(a,s)$ is the size distribution of composition $c$ at position $s$ along the sightline. $\Delta C_{\rm abs}$ is the difference between absorption cross sections for radiation polarized perpendicular and parallel to the grain's alignment axis, the polarization analogue of the intensity opacity. $T_c(a,s)$ is the temperature set by the balance between absorption of the local interstellar radiation field and thermal emission. $R_c$ is the Rayleigh reduction factor describing imperfect alignment, which radiative-torque alignment theory makes size and environment dependent. $\psi_E$ and $\gamma$ are the field-controlled polarization direction and inclination.

Several facts follow. The grains that emit most of the polarized power in the CMB bands are large, with sizes above roughly 0.05 to 0.1 μm, non-spherical, aligned, and in approximate radiative equilibrium at 18 to 20 K in the diffuse medium. Nanoparticles and PAHs are stochastically heated and unaligned, consistent with the absence of polarized anomalous microwave emission. At millimetre wavelengths the emitting grains are deep in the Rayleigh limit, where the cross section per unit volume is independent of size and set by the material's dielectric function and the grain shape. So knowing the size distribution does not give the spectral index. The value of $\beta$ between 100 and 400 GHz is a material property: it depends on composition, on optical constants that are themselves temperature dependent in amorphous solids, on shape and porosity, on aggregation in dense gas, and on the mixture of temperatures inside a beam, which flattens the effective SED and adds curvature.

What size information does buy is threefold. First, alignment is size dependent, so the size distribution together with an alignment model gives the fraction of the mass that is aligned and hence $p_k$, the polarized-to-total ratio of a layer, and how it differs between layers. Second, the equilibrium temperature depends weakly on size and strongly on the radiation field, and a 3D density map plus the 3D stellar distribution gives that field at each layer's position, so layers embedded in denser gas are predicted to be colder. Third, the polarized SED can differ from the intensity SED because the aligned grains are a subset with their own composition. Planck Intermediate Results XXII measured a polarization spectral index of 1.59 ± 0.02 against 1.51 ± 0.01 in intensity between 100 and 353 GHz, and Planck 2018 XI found 1.53 ± 0.02 for polarization at 19.6 K. Ritacco et al. 2023 argued from the lowest multipoles that models based on total intensity underestimate the complexity of the polarized foreground by a significant factor.

Grain models translate these facts into layer-dependent priors. Draine & Fraisse 2009 computed polarized far-infrared emission for silicate and carbon mixtures, where different $\beta$ for the two materials predicts a frequency-dependent polarization fraction. The observed near-constancy of the polarization fraction between 100 and 353 GHz, together with the ratio of optical to far-infrared polarization, led Hensley & Draine 2023 to the Astrodust+PAH model, in which all grains larger than about 0.02 μm share one composite material shaped as 1.4:1 oblate spheroids. Guillet et al. 2018 built DustEM models with polarized extinction and emission for translucent sightlines and found that standard astrosilicates are not emissive enough to reproduce the observed ratio of polarized emission to starlight polarization. THEMIS 2.0 by Ysard et al. 2024 uses laboratory silicate optical constants and does not require perfect alignment. These models disagree in detail, which is exactly why they are useful as a family: they define plausible relations between environment and the polarized $(\beta,T,{\rm curvature})$ of a layer, to be used as priors with widths, never as fixed inputs.

The observational handles on a layer's grain population are the extinction curve and $R_V$ for the size distribution, the peak wavelength of the Serkowski curve for the size of the aligned grains, which grows in denser gas, the starlight polarization efficiency $p_V/E(B-V)$ for alignment times inclination, with a maximum near 13% per magnitude found by Panopoulou et al. 2019, and the far-infrared SED for temperature and intensity $\beta$. None of these uniquely recovers $n(a)$. Together they bound which grains are aligned and how a layer's polarized SED is likely to differ from its neighbours'.

The empirical bridge between starlight and emission is the Planck Intermediate Results XXI ratio between the 353 GHz polarized intensity and the visible polarization fraction, $R_{P/p}=5.4$ MJy sr$^{-1}$ with statistical and systematic uncertainties of 0.2 and 0.3, and the ratio of polarization fractions $R_{S/V}=4.2$. Because the same aligned grains produce both dichroic extinction and polarized emission, both observables carry the factor $R\cos^2\gamma$ and the ratio is nearly independent of inclination and alignment efficiency. This is what allows a per-layer starlight polarization degree to be converted into a per-layer polarized intensity at 353 GHz.

3.3 The sky-projected field: in which direction each layer contributes

Each layer's polarization angle is set by its plane-of-sky field. Emission is polarized perpendicular to the projected field and transmitted starlight is polarized approximately parallel to it. The Planck 353 GHz angle map is the integrated result of all layers, and using it as the per-layer field solution would be circular. Three routes give distance-resolved angles.

Starlight polarization with Gaia distances. The polarization of a star is the vector sum of the contributions of all dust in front of it, so the cumulative Stokes parameters $q(d)$ and $u(d)$ along a sightline form a step function. Each step is a cloud, and the increment is that cloud's Stokes vector, which gives both its angle and its polarization degree. Recovering the steps requires many stars per resolution element with reliable parallaxes and a statistical inversion. Panopoulou et al. 2019 demonstrated this on one diffuse sightline with RoboPol and Gaia DR2, separating two clouds whose mean projected-field orientations differ by 60 degrees. Pelgrims et al. 2023 formalized the Bayesian line-of-sight inversion, with thin layers, nested sampling and model comparison over the number of clouds, and showed it works when a cloud induces polarization above roughly 0.1%. Pelgrims et al. 2024 applied it over four square degrees near Galactic coordinates (103.3°, 22.3°) with 1448 stars, finding one to three clouds per sightline, two over most of the area, and a projected-field orientation that changes with distance along most sightlines. PASIPHAE aims to scale this to more than 3.5 million stars over more than 10,000 square degrees at high latitude, about 360 stars per square degree, using wide-field polarimeters in both hemispheres, and its stated purpose is the tomography needed here.

Per-layer polarized emission from starlight. Combining a layer's incremental starlight polarization with $R_{P/p}$ gives that layer's polarized intensity at 353 GHz, and the incremental angle rotated by 90 degrees gives its emission angle. The 3D extinction map independently gives the same layer's column, so the two datasets cross-check each other. The layer's polarization fraction is their ratio, and its shortfall from the maximum bounds $\cos^2\gamma_k$ times the alignment efficiency.

H I orientation per velocity component. Clark & Hensley 2019 showed that H I filaments align with the projected field and built $Q$ and $U$ templates per velocity channel from HI4PI that correlate with Planck 353 GHz at better than 0.75 with no free parameters, reproducing the E/B ratio. Where a velocity component can be associated with a distance layer, its orientation is a proxy angle for that layer. It is a proxy, not a measurement.

What projected-field data do not give: the inclination $\gamma$, which scales each layer's polarized amplitude by $\cos^2\gamma$; the field strength; and the alignment efficiency. The 180-degree ambiguity of linear polarization is harmless because only $e^{2i\psi}$ enters, but the sign convention for $U$ (IAU versus HEALPix) must be handled explicitly. The magnetic-field map does not remove dust. It supplies the phases in the vector sum, and without SED contrast between layers those phases produce no chromatic effect at all.

4. A concrete dedusting strategy

The workflow below can be reconstructed from the equations above.

1. Segment sightlines into layers. From the 3D dust posterior, compute the differential extinction profile along every sightline in the region and identify layers $k=1\ldots K(\hat n)$ as structures that are coherent across neighbouring sightlines. Carry the posterior samples, so each layer column $E_k(\hat n)$ comes with covariance across pixels and layers. Reconcile with H I components where their kinematics are distinct, without equating velocity with distance.

2. Calibrate column to emission. Regress the far-infrared and submillimetre intensity maps, Planck 353 to 857 GHz, IRAS 100 μm and the GNILC dust maps, on the layer templates $E_k(\hat n)$, fitting per layer an emissivity per unit reddening, a temperature and an intensity $\beta$. The layers' different morphologies make this identifiable. Priors on $T_k$ come from the radiation field at the layer's position.

3. Set layer polarization priors. From starlight tomography obtain each layer's incremental Stokes vector: the emission angle $\psi_k=\psi_{{\rm stars},k}+90^\circ$ and the polarized intensity $P_{353,k}=R_{P/p}\,p_{V,k}$, with $R_{P/p}$ a per-layer nuisance parameter whose prior width reflects the measured scatter. Where stars are too few, use the H I orientation of the associated velocity component with a broad prior.

4. Set layer SED priors. Assign each layer a prior on its polarized $(\beta_{P,k},T_k,{\rm curvature})$ from the grain-model family conditioned on its environment: density from the 3D map, radiation field, $R_V$ and Serkowski peak wavelength where measured, and the intensity SED from step 2. Allow the polarized $\beta$ to differ from the intensity $\beta$. Prior widths must reflect the spread between grain models, not the precision of a single model.

5. Forward-project every layer into every band. Compute the mixing coefficients $M_{\nu k}(\theta_k)$ through the actual bandpasses and unit conversions, convolve each layer with the instrument beam, and predict

where $\mathbf t_k$ is the layer's $(Q,U)$ map at $\nu_0$, the CMB is frequency independent in thermodynamic units, and synchrotron is a power law with its own curvature. Because the layers enter as complex maps with distinct SEDs, the predicted dust cross-spectra are

with $C_\ell^{jk}$ the cross-spectra between layer maps. The correlation ratio, and its dependence on frequency pair, multipole and sky position, are outputs of the model rather than a single fitted $\Delta_d$. The same machinery yields, from the auxiliary data alone, a predicted per-pixel angle rotation $\delta\psi(\nu_1,\nu_2)$ and a predicted unpredictable fraction of the polarized amplitude after the best scalar rescaling. Both can be tested against Planck 217×353 GHz before any $r$ analysis.

6. Infer jointly in $Q$ and $U$. Sample the posterior of the CMB maps, the layer maps $\mathbf t_k$, the SED parameters $\theta_k$, the synchrotron component, and the calibration and bandpass nuisances, given the microwave data and the auxiliary priors. The auxiliary maps enter as priors with covariance, never as exact subtraction templates: the microwave data are allowed to move layer amplitudes and angles within the tomographic uncertainty. Parametric methods, Commander-style Gibbs sampling or FGBuster-style maximum likelihood, take the layer SEDs directly. Blind methods such as SMICA, GNILC and constrained ILCs can use the layer templates as deprojection constraints or to define weights. Marginalize the nuisances into the cleaned BB spectrum and the $r$ likelihood. In a power-spectrum likelihood the layer cross-spectra $C_\ell^{jk}$ become the nuisance parameters, with priors from the auxiliary data.

7. Reduce where layers are not identifiable. Identifiability is a counting and signal-to-noise problem. Each pixel supplies $2N_\nu$ numbers, which must constrain two CMB amplitudes, two synchrotron amplitudes and $2K$ dust amplitudes, plus SED parameters that can be shared across patches. Two layers are feasible with six or more polarized bands when their SEDs differ, and the auxiliary priors regularize the near-degenerate direction along which two similar-SED layers trade amplitude. When more layers are needed, either group nearby clouds into a small number of effective SED and angle components, or replace the layer sum by physically informed spectral moments,

with separate moment maps for $Q$ and $U$ as in Vacher et al. 2022, or preferably the spin-moment formalism of Vacher et al. 2023, in which the moments are spin-2 objects and a rotation of the polarization vector with frequency emerges from the first-order moment whenever spectral parameters and angles vary together. Vacher et al. 2022 showed on LiteBIRD-like simulations that moments in both $\beta$ and $T$ prevent biases on $r$ with $\sigma(r)\approx 9\times10^{-4}$. The 3D and field data then supply priors on the moment amplitudes rather than leaving them free.

8. Mask or down-weight. The predicted rotation and unpredictable fraction define a line-of-sight complexity score. Sightlines with several layers, misaligned fields with large $|\sin 2\Delta|$, and strong SED contrast get down-weighted or masked, and deep fields can be chosen where a single layer dominates. This is the mitigation Tassis & Pavlidou proposed and it is a legitimate part of the strategy, not a failure of it.

9. Validate. Compare predicted and observed 217 to 353 GHz angle differences, which is the Pelgrims et al. 2021 test, predicted and observed correlation ratios, polarized versus intensity SEDs per layer, and E/B ratios and TE correlations. Run the pipeline on multi-layer simulations of the kind built by Martínez-Solaeche et al. 2018 and McBride et al. 2023. The latter showed that two-parameter modified-blackbody fits recover the CMB without bias unless both the SED properties and the polarization angles vary along the line of sight, and that the fitted $\beta$ and $T$ are then unrepresentative of the physical dust, which is exactly the regime the layered model is built for.

5. Limitations and residual systematics

What improves: extrapolation bias is reduced because the frequency dependence of the pattern is modelled rather than assumed away; mixing parameters that a purely microwave fit cannot separate, such as two layers with similar SEDs or angle against amplitude, are constrained by independent data; the model predicts frequency-dependent angle rotation and decorrelation as functions of position, multipole and frequency pair; and dangerous sightlines are identified before they bias $r$.

What remains:

• Extinction is not emission. The conversion per layer carries environment-dependent uncertainty. The maps have finite angular resolution, 14 arcminutes for Edenhofer et al., and radial smearing. They miss very diffuse dust and dust beyond the last well-measured stars, and they assume an extinction law.
• Starlight tomography is expensive and incomplete. Resolving clouds that polarize at the 0.1 to 1% level needs hundreds of stars per resolution element at 0.1 to 0.2% polarimetric precision. Intrinsic stellar polarization, the unknown inclination, the scatter of $R_{P/p}$, and the difference between pencil-beam stars and beam-averaged emission, where turbulence inside the beam depolarizes the emission but not the individual stars, all limit the per-layer amplitude prior.
• Grain physics gives priors, not answers. The relations between environment and polarized $\beta$ and $T$ differ between models. Temperature-dependent opacities, several temperatures within a layer, possible magnetic dipole emission from iron inclusions below 100 GHz, which Hensley & Bull 2018 flagged as one of the two most damaging complexities alongside line-of-sight effects, and SED curvature are not fixed by any auxiliary dataset. A prior narrower than the model spread risks biasing $r$.
• Auxiliary surveys do not replace broad microwave frequency coverage. The SED parameters are still measured by the channels between about 30 and 400 GHz. The highest-frequency channels, Planck 353 GHz, the LiteBIRD 337 and 402 GHz bands, and submillimetre bands from CCAT-prime, remain essential to anchor the dust.
• The method does nothing for the other terms. Synchrotron and its correlation with dust, lensing B-modes, instrumental polarization, bandpass and gain mismatch, polarization-angle calibration, beam leakage, and E-to-B leakage from the cut sky are untouched. Several of these produce spurious decorrelation, so a measured $\mathcal R_\ell<1$ can be attributed to astrophysics only after the systematics budget is closed.
• Removing thermal dust does not produce a primordial B-mode map. What remains after dedusting is lensing plus residual foregrounds plus noise, and the constraint on $r$ comes from the likelihood after marginalization over every nuisance introduced above.

The scientific value of the approach lies in the width and honesty of the resulting posterior on $r$, and in the model's validated ability to predict the frequency dependence of the dust polarization pattern. Visual agreement between a layered model and the Planck 353 GHz map proves little, because the integrated map is what the model is constructed to reproduce. The recent Simons Observatory study makes the stakes concrete: spatial variation of $\beta$ with a dispersion near 0.3 can bias $r$ by about 0.03 in low-order moment expansions, and more flexible foreground models are needed at the noise levels expected. A distance-resolved, physically parametrized model is the most constrained form of that flexibility, because every extra degree of freedom it introduces corresponds to a structure that can be seen in stars, gas or extinction.

Deviations from brief

No scientific error was found in the brief. Four refinements were made:

• The brief's cited Vacher et al. paper (arXiv:2111.07742) is the $Q$/$U$ moment-expansion study. The spin-moment formalism that makes the polarized moments rotation invariant is Vacher, Chluba, Aumont, Rotti & Montier 2023 (arXiv:2205.01049); the essay cites both and recommends the latter.
• The statement that grain size does not set $\beta$ was made explicit through the Rayleigh limit, where the cross section per unit volume is size independent.
• The brief says projected-field data do not determine inclination or alignment efficiency, which is true. The essay adds that the empirical ratio $R_{P/p}$ from Planck Intermediate Results XXI makes the layer's polarized intensity at 353 GHz nearly independent of those two unknowns, so they need not be determined separately for the amplitude prior.
• The brief attributes the difference between Planck Intermediate Results L and Planck 2018 XI to data treatment and statistics. The essay adds the specific reasons identified by Sheehy & Slosar 2018: noise bias in the split cross-spectra, evidence for instrumental systematics, and neglected correlations between overlapping regions.

Key references

• Tassis & Pavlidou 2015, MNRAS Letters, arXiv:1410.8136. Poh & Dodelson 2017, PRD 95, 103511, arXiv:1606.08922.
• Planck Intermediate Results L, A&A 599, A51, arXiv:1606.07335. Sheehy & Slosar 2018, PRD 97, 043522, arXiv:1709.09729. Planck 2018 XI, A&A 641, A11, arXiv:1801.04945.
• Pelgrims et al. 2021, A&A 647, A16, arXiv:2101.09291. Pelgrims et al. 2023, A&A 670, A164, arXiv:2208.02278. Pelgrims et al. 2024, A&A, arXiv:2404.10821. Panopoulou et al. 2019, ApJ 872, 56, arXiv:1809.09804. PASIPHAE: Tassis et al. 2018, arXiv:1810.05652.
• Edenhofer et al. 2024, A&A, arXiv:2308.01295. Panopoulou & Lenz 2020, ApJ 902, 120, arXiv:2004.00647. Clark & Hensley 2019, ApJ 887, 136, arXiv:1909.11673.
• Planck Intermediate Results XXI, arXiv:1405.0873, and XXII, arXiv:1405.0874. Panopoulou et al. 2019, A&A 624, L8, arXiv:1903.09684. Ritacco et al. 2023, A&A 670, A163, arXiv:2206.07671.
• Draine & Fraisse 2009, ApJ 696, 1. Guillet et al. 2018, A&A 610, A16, arXiv:1710.04598. Hensley & Draine 2023, ApJ 948, 55, arXiv:2208.12365. Ysard et al. 2024, THEMIS 2.0, arXiv:2401.07739.
• Vacher et al. 2022, A&A 660, A111, arXiv:2111.07742. Vacher et al. 2023, A&A 669, A5, arXiv:2205.01049. Hensley & Bull 2018, ApJ, arXiv:1709.07897. Martínez-Solaeche et al. 2018, MNRAS, arXiv:1706.04162. McBride, Bull & Hensley 2023, MNRAS 519, arXiv:2207.14213.
• BICEP/Keck XIII (BK18), PRL 127, 151301, arXiv:2110.00483. BICEP/Keck XVI, ApJ 945, 72, arXiv:2210.05684. Liu et al. 2025, Simons Observatory dust complexity, arXiv:2508.00073.