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

Multi-layer polarized dust: 3D extinction, grain physics, and distance-resolved fields as a B-mode foreground model

A single 353 GHz dust template fails when different clouds along one sightline have different polarized spectral energy distributions (SEDs) and projected magnetic-field angles. The operational problem is not a mysterious “spectral-index decorrelation” of a single $\beta_d$, but frequency decorrelation: the polarized dust maps at two frequencies are not exact scalar multiples of one another. Three-dimensional dust information, grain-physics constraints, and distance-resolved field orientations can replace that single template with a physically informed, multi-layer model of Stokes $Q$ and $U$, whose uncertainties are propagated into the recovered CMB $B$-mode likelihood. The result is uncertainty reduction and principled down-weighting of unmodelable sightlines, not exact “dedusting.”

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

Inflation can source a stochastic background of primordial gravitational waves. Those tensor modes imprint a curl-like pattern in CMB polarization, the primordial $B$-mode, whose amplitude is customarily reported as the tensor-to-scalar ratio $r$. Current limits sit near $r \lesssim 0.03$ (95% CL, depending on data combination); LiteBIRD, CMB-S4, and the Simons Observatory target $\sigma(r)$ of order $10^{-3}$ or below. At that sensitivity the limiting systematic is not raw map noise. It is Galactic foregrounds, and at the frequencies used to hunt for the recombination peak ($\sim 90$–$220$ GHz) the dominant polarized foreground is thermal emission from aligned interstellar dust.

The relevant Galactic physics is distinct from the CMB $B$-mode. Aspherical dust grains in the diffuse ISM have their short axes statistically aligned with the local Galactic magnetic field (radiative-torque alignment being the leading mechanism in the diffuse gas). Thermal emission is then linearly polarized with the electric vector perpendicular to the plane-of-sky (POS) projection of $\mathbf{B}$. The resulting Stokes $Q$ and $U$ maps, after a spin-2 spherical-harmonic transform, contain both gradient-like $E$ modes and curl-like $B$ modes. Those Galactic $B$ modes are a property of the magnetized ISM’s projected geometry; they are not a CMB $B$-mode, and they are not a property of the magnetic field itself. At 353 GHz they already dominate the sky at high latitude on the angular scales of the recombination bump ($\ell \sim 80$). Planck Collaboration XI (2020) measure a mean polarized-dust spectral index $\beta_d^P = 1.53 \pm 0.02$ for a dust temperature of $19.6\,\mathrm{K}$, with $\beta_d^P - \beta_d^I = 0.05 \pm 0.03$, and an $E/B$ power ratio $\langle A^{BB}/A^{EE}\rangle \simeq 0.5$ together with a positive $TE$ correlation—signatures of the filament–field alignment of the dusty ISM, not of primordial tensors.

The conventional cleaning strategy is therefore: measure dust where it dominates (Planck 353 GHz, or a balloon/ground channel near $220$–$280$ GHz), assume a modified blackbody (MBB) SED with a single $\beta_d$ and $T_d$ per pixel or per patch, and scale that template into the CMB bands. In thermodynamic CMB units the CMB itself is achromatic, so a perfectly scaled dust template can in principle be separated from the CMB with a small number of frequencies. That strategy is exact only if the polarized dust sky at frequency $\nu$ is a scalar multiple of the sky at a reference frequency $\nu_0$. Whenever that fails, a residual with a dust-like angular spectrum leaks into the reconstructed CMB $B$ modes and biases $r$. BICEP2/Keck × Planck already showed that a $\sim 10\%$ loss of 150×353 correlation can shift $r$ by $\sim 0.02$. Poh & Dodelson (2017) estimated that line-of-sight (LOS) “extrapolation noise” of this kind, though subdominant to Planck’s uncertainties, becomes potentially important for experiments targeting $r \lesssim 0.0015$. Next-generation analyses cannot treat a one-template MBB as an adequate error model.

2. What “spectral-index decorrelation” means operationally

The phrase in the question is imprecise. $\beta_d$ is a fit parameter in an MBB,

and “spectral-index decorrelation” is often used as shorthand for the observation that a single $\beta_d$ does not describe all pixels, all frequencies, and both Stokes parameters at once. What component separation actually needs is a statement about maps. Write the complex polarized brightness

A single-template model asserts

with a real scalar $f$ (possibly pixel-dependent). Frequency decorrelation is the failure of that relation: $Q$ and $U$ at $\nu_1$ are not the same linear combination of $Q$ and $U$ at $\nu_2$, so the polarization angle

and the polarized fraction $p = |P|/I$ both become frequency-dependent. Equivalently, in harmonic space one defines the BB correlation ratio

Values $\rho_\ell < 1$ mean that the $B$-mode maps at the two frequencies are not perfectly correlated. That is a useful diagnostic, not a physical identification. The same statistic responds to (i) LOS mixing of misaligned clouds with different SEDs, (ii) spatial SED variations from pixel to pixel even if each pixel is a single MBB, (iii) CMB chance correlations, (iv) noise, calibration, and bandpass mismatch, (v) beam differences, and (vi) additional components (synchrotron, AME, magnetic dipole emission). Interpreting $\rho_\ell < 1$ as “the spectral index is decorrelated” collapses those mechanisms.

The residual that biases $r$ is a map-level object. If a parametric fit subtracts $\hat P_\nu = \hat f(\nu) P_{353}$ and the true $P_\nu$ is not of that form, the difference $\delta P_\nu = P_\nu - \hat P_\nu$ is a polarized field whose $B$-mode power adds to the CMB. Because dust $C_\ell^{BB}$ is red ($\alpha_{BB} \simeq -2.5$ near $\ell = 80$; Planck Collaboration XI 2020) a small fractional residual at 150 GHz can still exceed a primordial signal with $r \sim 10^{-3}$. The bias is worse if $Q$ and $U$ are forced to share a single SED parameter: angle rotation with frequency is then absorbed into an incorrect amplitude, and the residual is not even a rescaled copy of the 353 GHz template.

Empirical status is estimator- and region-dependent, and should not be quoted as a universal amplitude. Planck Intermediate Results L (2017) reported a departure of $\mathcal{R}_\ell^{BB}(353,217)$ from the no-variation expectation, stronger toward lower column density, and noted that either spatial SED variation or polarization-angle variation could produce it. Subsequent work (Sheehy & Slosar 2018; Planck Collaboration XI 2020) found that residual systematics and noise undermine that detection in a power-spectrum analysis of the PR3 maps: Planck 2018 reports no significant loss of correlation among the polarized HFI channels and sets lower limits on $\rho_\ell$ tighter than the 217×353 ratio alone. Those limits, if they applied to the small, clean patches used by ground-based experiments, would suggest that frequency decorrelation need not spoil an $r \simeq 0.01$ search at the recombination peak—while remaining inadequate for $r \sim 10^{-3}$. Independently, Pelgrims et al. (2021) selected high-latitude pixels using only H I: a target sample with multiple clouds and POS field angles misaligned by $\gtrsim 60^\circ$, and a control sample with a single dominant cloud. After CMB subtraction they detected a larger scatter in $\psi_{353}-\psi_{217}$ in the target sample, with a $\sim 15\%$ intrinsic variation in the 353/217 polarized-intensity ratio between clouds sufficient to match the signal. That is pixel-level LOS angle decorrelation, not a sky-averaged $\rho_\ell$. The two statements are compatible: a minority of sightlines can rotate while a large-area cross-spectrum remains consistent with $\rho_\ell \simeq 1$.

3. The line-of-sight mechanism

Discrete layers

For optically thin thermal dust, the observed polarized brightness is a vector sum along the LOS. Discretizing into clouds or voxels $k$,

Here $A_k$ is a column-density (or optical-depth) amplitude at a reference frequency $\nu_0$; $p_k$ is the intrinsic polarization fraction that the layer would have if $\mathbf{B}$ lay in the POS and the grains were as aligned as they are; $\gamma_k$ is the angle between $\mathbf{B}$ and the POS, so that the geometric factor $\cos^2\gamma_k$ vanishes for a field along the LOS (equivalently $\sin^2 i$ if $i$ is the inclination to the LOS); and $\psi_{E,k}$ is the emission polarization angle of that layer. Thermal-emission polarization is perpendicular to the POS field, so

(modulo $\pi$; orientations, not directions). The MBB factor in specific intensity is

with $B_\nu(T)$ the Planck function. Conversion into thermodynamic CMB units multiplies $f_\nu$ by the usual $\bigl(\partial B_\nu/\partial T|_{T_{\mathrm{CMB}}}\bigr)^{-1}$ and by the bandpass, but does not make the CMB chromatic: in those units the CMB Stokes parameters are independent of $\nu$.

The quantities $p_k$ and $\gamma_k$ are not interchangeable with $\beta_k$. Alignment efficiency and field inclination set the amplitude of $P$ at a given column; they produce frequency-dependent $Q/U$ only if they differ among layers whose SEDs also differ.

Two clouds, two necessary ingredients

Take two clouds with polarized amplitudes

and emission angles $\psi_1,\psi_2$. Then

The ratio of polarized weights $\rho_\nu \equiv w_1(\nu)/w_2(\nu)$ is frequency-dependent whenever the clouds differ in $T_d$, $\beta_d$, grain mix, or (frequency-dependent) alignment. The resultant angle satisfies

and $\psi_{\mathrm{tot}}$ swings by as much as the misalignment $\alpha \equiv |\psi_1-\psi_2|$ as $\rho_\nu$ goes from $\ll 1$ to $\gg 1$. Tassis & Pavlidou (2015) showed that for $\alpha \to 90^\circ$ the Stokes vectors in the $(Q,U)$ plane are antiparallel: $p_{\mathrm{tot}}$ passes through zero at $\rho_\nu = 1$, and $\psi_{\mathrm{tot}}$ jumps by $90^\circ$ if $\rho_\nu-1$ changes sign between frequencies. For $\alpha \lesssim 60^\circ$ the same temperature contrasts produce only a few degrees of angle shift and $\lesssim 10\%$ changes in $p$ between 150 and 353 GHz, which is why a misalignment cut is a useful mask.

Two ingredients are necessary for strong LOS angle decorrelation:

1. The ratio $f_{\nu,1}/f_{\nu,2}$ (hence $\rho_\nu$) changes with frequency, because the clouds differ in $T_d$, $\beta_d$, grain mix, or alignment.
2. Their POS field (equivalently emission) angles are misaligned.

The causal distinction is sharp.

• If the SEDs are identical, $w_1/w_2$ is independent of $\nu$. A geometrically complicated sightline still depolarizes ($p_{\mathrm{tot}} < p_k$ whenever $\alpha \neq 0$), and a tangled field can convert $E$ into $B$, but $P_\nu = f_\nu \times$ (a frequency-independent complex number). A 353 GHz template still scales. This is LOS depolarization without frequency decorrelation.
• If the angles are identical, $\psi_{\mathrm{tot}}$ is fixed. Mixing of different MBBs produces curvature: the effective SED is not itself an MBB, so a one-parameter $\beta_d$ fit misestimates the amplitude. That is a calibration error in the scalar $f(\nu)$, not a rotating resultant. $Q$ and $U$ still share a SED, so the maps remain scalar multiples up to that curvature.
• Only the product of SED contrast and angle misalignment makes $Q(\nu)$ and $U(\nu)$ scale differently, so that $P_{\nu_1}$ is not a real multiple of $P_{\nu_2}$.

Temperature contrast is enough to move $\rho_\nu$ even at fixed $\beta_d$. Typical high-latitude dust is not in the Rayleigh–Jeans limit across 150–353 GHz. For $T_d \simeq 20\,\mathrm{K}$, $h\nu/kT_d \simeq 0.36$ at 150 GHz and $\simeq 0.85$ at 353 GHz, so

still depends on $\nu$. Intermediate-velocity clouds (IVCs) are observed to be warmer than local low-velocity gas; that is a concrete astrophysical source of $\rho_\nu(\nu)$ along high-latitude sightlines (Planck Collaboration XI 2014; Panopoulou & Lenz 2020). Additional contrast comes from $\beta_d$ differences, from a different mix of aligned materials, and from any frequency dependence of $p_k$ itself.

A third, more mundane, source of $\rho_\ell < 1$ is inter-pixel SED variation: even a sky of perfect single-layer MBBs with a spatially varying $\beta_d(\hat n)$ will not be a scalar multiple of one 353 GHz map once one averages over a patch. LOS angle mixing is therefore not the only physical path to decorrelation; it is the path that a 3D field model is uniquely able to address.

4. Grain physics: $\beta_d$ is an emergent number

The MBB is a fitting function. The optically thin polarized emission is

The sum runs over compositional populations $c$ (astrodust/silicate, carbonaceous, PAHs, …). $n_c(a,s)\,\mathrm{d}a$ is the number density of grains with equivalent radius in $[a,a+\mathrm{d}a]$; $R_c(a,s)$ is the aligned fraction (including imperfect Davis–Greenstein or radiative-torque alignment); $\Delta C_{\mathrm{abs},c}$ is the difference of absorption cross-sections for $E$ parallel and perpendicular to the grain’s long axis, which by Kirchhoff’s law is also the polarized emission cross-section; and $T_c(a,s)$ is the grain temperature. For large grains $T_c$ is an equilibrium temperature set by the local interstellar radiation field (ISRF); for nanoparticles it is a stochastic distribution.

Several facts follow, and they are what “knowledge of the grain size distribution” can and cannot provide.

Large equilibrium grains dominate polarized thermal emission in CMB bands. Alignment functions inferred from the Serkowski law and from UV polarization (Kim & Martin 1995; Draine & Fraisse 2009; Guillet et al. 2018; Hensley & Draine 2023) have $R_c(a) \approx 0$ for $a \lesssim 0.05$–$0.1\,\mu\mathrm{m}$ and rise toward the big-grain end. Very small stochastically heated grains and PAHs dominate the mid-infrared unpolarized SED and, if spinning, the anomalous microwave emission near $20$–$40$ GHz; they are not interchangeable with the aligned population that produces 100–353 GHz polarization. Magnetic dipole emission from ferromagnetic inclusions is a further, poorly constrained, polarized continuum that can flatten the millimetre SED (Draine & Hensley 2013; Planck Collaboration Int. XVII 2014). A component-separation model that uses a spinning-dust or PAH template as a proxy for polarized thermal dust is solving a different problem.

The size distribution $n_c(a)$, the dielectric function, and the grain shape determine both $\Delta C_{\mathrm{abs}}(\nu)$ and the equilibrium $T_c(a)$. In the Rayleigh limit $C_{\mathrm{abs}} \propto a^3\nu\,\mathrm{Im}(\alpha)$, so the millimetre opacity spectral index is set by the imaginary part of the dielectric function, not by a power-law slope of $n(a)$. Smaller big grains, at fixed ISRF, run hotter because they radiate less efficiently in the FIR; a size mix is therefore a temperature mix, and a temperature mix is already an SED distortion relative to a single MBB (Chluba, Hill & Abitbol 2017). Porosity, ice mantles, and coagulation in denser gas change optical constants and the FIR emissivity per hydrogen atom. Alignment physics (RATs; Lazarian & Hoang 2007) couples $R_c(a)$ to the local radiation anisotropy, field strength, and gas density, so the polarized SED can differ from the total-intensity SED even in one voxel. That is a microphysical reason for $\beta_d^P \neq \beta_d^I$.

Observational constraints on $n_c(a)$ are degenerate proxies, not inversions of a unique size distribution:

• Optical/NIR extinction curves and $R_V = A_V/E(B-V)$ are sensitive to grains with $a \sim \lambda/2\pi$, i.e. the same large grains that polarize, but also to the unaligned population and to composition.
• Infrared SEDs (IRAS, Herschel, Planck intensity) constrain $T_d$ and $\beta_d$ of the emitting mix, including unaligned grains and, at shorter wavelengths, stochastically heated particles.
• Wavelength-dependent starlight polarization, including the Serkowski peak $\lambda_{\max}$, traces the size of the aligned population; a shift in $\lambda_{\max}$ is often read as a shift in characteristic aligned size, but shape, composition, and the alignment function $R_c(a)$ enter as well.
• The ratios $P_{353}/p_V$ and $p_{353}/(p_V/\tau_V)$ (Planck Collaboration Int. XXI 2015) constrain the FIR polarized emissivity of the aligned grains relative to optical dichroism, and have already forced revisions of silicate-only models (Guillet et al. 2018; Hensley & Draine 2023).

The Astrodust+PAH model (Hensley & Draine 2023) is a useful prior rather than a unique truth: a single composite material for $a \gtrsim 0.02\,\mu\mathrm{m}$, plus a PAH nanoparticle population, is motivated by the near frequency independence of the FIR polarization fraction from $\sim 250\,\mu\mathrm{m}$ to $\sim 3\,\mathrm{mm}$. DustEM models with aligned astrosilicates (and optional aligned carbons) can be tuned to Planck $P_{353}/p_V$ by changing elongation, porosity, or carbon inclusions (Guillet et al. 2018), but they still over-predict the steepness of the polarized millimetre SED unless the optical constants are adjusted. The lesson for component separation is that a physically informed prior on $(T_d,\beta_d,p)$ per layer is tighter than an unbounded pixel-by-pixel MBB, and is not a delta function at one Galactic $\beta_d$. Size distribution alone does not determine $\beta_d$.

5. How each ingredient helps—and what it does not constrain

5.1 Three-dimensional dust maps

A hypothetical true 3D emissivity cube would supply the scalar volume emissivity $j_\nu(\mathbf{x})$ and, with a field model, the polarized emissivity. That cube does not exist at microwave frequencies. What exists are reconstructions of extinction or reddening as a function of distance, inferred from stellar photometry and Gaia parallaxes: Bayestar (Green et al. 2019 and successors, Pan-STARRS/2MASS/Gaia), the Gaia XP map of Edenhofer et al. (2024) with parsec-scale distance resolution out to $1.25\,\mathrm{kpc}$ at up to $14'$ angular resolution, and deeper plane maps that combine optical and infrared photometry. Differential extinction $\partial A_V/\partial s$ locates clouds and weights them by dust column. That is the $A_k(\hat n)$ in the layer sum, up to a conversion factor.

The conversion is the first major uncertainty. Optical reddening measures the extinction cross-section of the grain mix at $\sim 0.5\,\mu\mathrm{m}$, not the millimetre opacity and not the aligned fraction. Going from $A_V(s)$ to $j_\nu(s)$ requires (i) a dust-to-extinction calibration (itself $R_V$- and composition-dependent), (ii) an ISRF to set $T_c(a)$, and (iii) a grain model for $\kappa_\nu$ and $\Delta C_{\mathrm{abs}}$. The ISRF is not the Mathis field everywhere: it is anisotropic near OB associations, attenuated inside clouds, and enhanced in the Gould Belt. High-latitude sightlines are dominated by nearby dust (scale height $\sim 100$–$150\,\mathrm{pc}$), so a $1.25\,\mathrm{kpc}$ map captures most of the local column, but not intermediate- and high-velocity halo clouds, and not dust beyond the stellar tracer volume. Smoothing, unmodelled small-scale structure, and distance errors in the 3D map all leak into $A_k$.

Used correctly, a 3D extinction map does three things a 353 GHz template cannot: it counts clouds along each sightline, it weights them, and it provides a geometric prior for where SED parameters are allowed to jump. Used incorrectly—as if it were a 3D microwave emission map—it simply relabels the one-template problem.

H I is a complementary, not identical, tracer. In the diffuse high-latitude ISM, $N_{\mathrm{HI}}$ correlates tightly with far-infrared dust emission (Boulanger et al. 1996; Lenz, Hensley & Doré 2017). Velocity-resolved H I (HI4PI, GALFA-H I) decomposes the LOS into kinematic components (Panopoulou & Lenz 2020). Kinematic distance is not geometric distance, and molecular and ionized phases are under-represented, but the combination of an H I cloud catalogue with 3D extinction is currently the most practical layerization of the high-latitude sky.

5.2 Grain-size distribution and SED model

Given layers, a grain model assigns priors on the per-layer parameters $\theta_k = \{T_k,\beta_k,\kappa_k,p_k,R_k,\ldots\}$, rather than imposing one SED on the full column. Extinction curves and $R_V$ in front of stars that sit behind a given layer, Herschel/Planck intensity SEDs of that layer (when it can be isolated), and $\lambda_{\max}$ from stars in the same distance slice all tighten those priors. The polarized millimetre SED is then predicted from $\Delta C_{\mathrm{abs}}(\nu) B_\nu[T(a)]$ integrated over the aligned size distribution, not from a free $\beta_d \in [1,2.5]$.

What this does not do is replace the microwave data. Optical constants in the millimetre are not measured in situ; they are extrapolated from FIR and from laboratory analogues. Alignment efficiency in a given cloud is not given by $n(a)$ alone. Temperature mixing inside a “layer” that is still a finite box in distance remains. The grain model’s role in the likelihood is to (i) correlate $T_d$ and $\beta_d$ in a physically allowed region, (ii) correlate polarized and unpolarized amplitudes, (iii) forbid using PAHs or spinning dust as a polarized thermal template, and (iv) predict, at least statistically, how much SED contrast is expected between, e.g., a warm IVC and a colder local cloud. Pelgrims et al. (2021) needed only $\sim 15\%$ contrast in the 353/217 polarized ratio; that is within the range of equilibrium-temperature differences at fixed $\beta_d$ for $T_d \sim 15$–$25\,\mathrm{K}$.

5.3 Sky-projected magnetic field, distance-resolved

An integrated Planck 353 GHz polarization angle is the vector sum. It cannot unmix layers. The quantity that enters $\psi_{E,k}$ is the POS orientation in that layer.

Starlight polarimetry supplies it, with two non-negotiable caveats. First, dichroic absorption polarizes starlight parallel to the POS field, whereas thermal emission is perpendicular; a $90^\circ$ rotation (and a consistent IAU versus COSMO convention) must be applied before comparing $\psi_\star$ to $\psi_{\mathrm{dust}}$. Second, stellar polarization is cumulative. A star at distance $d$ measures $\sum_{s<d} P_{\mathrm{abs}}(s)$. Layer angles come from differencing stars (or a tomographic inversion) across distance bins, not from a 2D polarization map. Pelgrims et al. (2023, 2024) demonstrate a Bayesian LOS inversion of stellar $q,u$ plus Gaia distances, producing the first degree-scale tomography of POS $\mathbf{B}$ in a $\sim 4\,\mathrm{deg}^2$ field; PASIPHAE (Tassis et al. 2018) is designed to extend that to high Galactic latitude at the stellar density required for CMB patches. Sparse stars, polarization noise at the $0.1\%$ level of the cleanest high-latitude fields, and the unknown inclination $\gamma_k$ remain. A projected field is an orientation modulo $180^\circ$; it does not give $|\mathbf{B}|$ or $\gamma_k$. The unknown inclination sits in the polarized amplitude through $\cos^2\gamma_k$ and is a prior, not a measurement.

H I morphology is a statistical field proxy, not a per-cloud magnetometer. Cold H I filaments align with the POS field (Clark et al. 2014, 2015); Clark & Hensley (2019) turn Rolling Hough Transform orientations into velocity-resolved Stokes $Q_{\mathrm{HI}}(v), U_{\mathrm{HI}}(v)$ that reproduce much of the Planck 353 GHz pattern after velocity integration. Pelgrims et al. (2021) used exactly that product, together with Panopoulou & Lenz cloud counting, to select misaligned sightlines. The proxy is validated on CNM-like structures and on the integrated sky; it need not hold in every WNM parcel, in molecular gas, or where velocity caustics masquerade as filaments. Faraday rotation and synchrotron polarization trace ionized-gas fields and cosmic-ray electrons; they are useful for Galactic-field modelling but are not a substitute for the dusty-phase POS angle.

A two-dimensional starlight or dust polarization map, however high its signal-to-noise, does not solve a LOS problem. The useful information is distance- or velocity-resolved.

6. Combining the three into a mixing matrix

6.1 From the emissivity integral to a data model

Collect the ingredients. Segment each sightline into $K$ layers (clouds or distance voxels) using 3D extinction and/or H I components, with amplitudes $A_k$ proportional to differential reddening, converted to millimetre optical depth through a grain/ISRF model. Assign each layer a polarization angle $\psi_{E,k}$ from stellar tomography and/or velocity-resolved H I, after the $90^\circ$ convention. Assign each layer SED parameters $\theta_k$ with grain-informed priors. The predicted dust Stokes vector at frequency $\nu$ and pixel $p$ is then

where $\mathbf{P}_{kp}^{\mathrm{dust}}$ is the layer’s polarized amplitude at a reference frequency (a two-vector set by $A_k p_k \cos^2\gamma_k$ and $\psi_{E,k}$), and $M_{\nu k}$ is the scalar mixing coefficient that carries that layer from $\nu_0$ to $\nu$ in thermodynamic CMB units, including the MBB (or a more microscopic $\Delta C_{\mathrm{abs}} B_\nu[T]$) and the bandpass. For a single layer with one angle, $Q$ and $U$ share $M_{\nu k}$; frequency-dependent angle appears only after summing layers.

The microwave data model is

with $s_p^{\mathrm{CMB}}$ achromatic in $K_{\mathrm{CMB}}$. Synchrotron has its own spectral index and, in general, its own LOS structure; it is a separate component, not a dust layer. Noise $n_{\nu p}$ includes correlated $Q$–$U$ covariance, beams, and calibration.

This is the same parametric mixing matrix used by Commander-like Gibbs samplers and by FGBuster-style maximum-likelihood separators (Eriksen et al.; Stompor, Errard and collaborators), except that the dust block is layered and the layer parameters are constrained by auxiliary data rather than being free per pixel. The implementation that matches the error budget is hierarchical Bayesian:

• auxiliary maps (3D extinction, stellar $q,u$, H I orientations, extinction curves) enter as probabilistic priors, with their own noise, selection functions, and model error;
• all microwave frequencies are fitted simultaneously in $Q$ and $U$, not by scaling a 353 GHz template;
• one marginalizes over $\theta_k$, over the conversion from reddening to millimetre opacity, over residual calibration and bandpass, and over the 3D-map distance errors;
• the CMB Stokes parameters (or their $B$-mode power, in a hybrid map/harmonic analysis) are sampled jointly, so that dust-model uncertainty inflates the $r$ posterior instead of hiding as a point-estimate residual.

When the number of physical layers exceeds what the frequency data can identify, compress. Two practical compressions, both informed by the 3D model rather than being blind, are:

1. Cloud templates. Collapse voxels that share a distance slice, a common $\psi_E$, and a common SED prior into one effective cloud. Tassis & Pavlidou’s two-cloud example is the minimal version; real high-latitude sightlines often have $\mathcal{N}_c \sim 1$–$3$ H I components (Panopoulou & Lenz 2020).
2. Moment expansion. Taylor-expand the MBB in $\beta_d$ and $T_d$ about a pivot (Chluba et al. 2017; Mangilli et al. 2021; Vacher et al. 2022). The first moments are weighted means of $\delta\beta$ and $\delta T$; the second moments include $\mathrm{Var}(\beta)$, $\mathrm{Var}(T)$, and $\mathrm{Cov}(\beta,T)$. Polarized moments must be defined on the spin-2 field, because angle mixing is a complex average of $P$, not of $I$ (Vacher et al. 2023). The 3D model predicts which moments are large: a two-cloud sightline with $\Delta\psi \sim 70^\circ$ and $\Delta T \sim 5\,\mathrm{K}$ has a specific first-order $T$ moment and a specific frequency-dependent $\psi$. LiteBIRD-oriented forecasts show that expanding in $\beta$ alone is not enough at $\sigma(r)\sim 10^{-3}$; temperature moments are required to keep $r$ unbiased (Vacher et al. 2022). A single scalar “decorrelation parameter” multiplying $C_\ell(\nu_1\times\nu_2)$ is a cruder compression, and can still bias $r$ (Mangilli et al. 2021).

6.2 Two products, not one

Where the posterior on $\{A_k,\theta_k,\psi_{E,k}\}$ is constrained—typically bright, simple, or well-starred sightlines—the multi-layer model subtracts dust better than a one-MBB template, because it can follow a frequency-dependent $\psi_E$ and a curved effective SED. Where it is not constrained—faint pixels, many misaligned clouds, no background stars, dust beyond the map volume—the same posterior flags the sightline. The second product is a principled mask or inverse-noise weight: down-weight pixels with $\mathcal{N}_c$ large and $\Delta\psi_{B}$ large, exactly the target sample of Pelgrims et al. (2021), and exactly the $\alpha \gtrsim 60^\circ$ cut anticipated by Tassis & Pavlidou (2015). Masking is not a failure of the method; it is the correct response to a likelihood that is genuinely broad. An analysis that always produces a cleaned map with no increase in the $r$ uncertainty is overconfident.

A practical pipeline, in outline:

1. Build layer weights $A_k(\hat n)$ from a 3D extinction cube (Edenhofer/Bayestar) plus H I cloud columns, with a conversion prior from a grain model.
2. Build layer angles $\psi_{E,k}$ from PASIPHAE-style tomography where stars exist, and from H I orientations in velocity slices elsewhere, with a hierarchical prior that the two tracers agree where they overlap.
3. Place grain-informed priors on $(T_k,\beta_k,p_k)$, allowing IVC/LVC systematic offsets.
4. Fit the microwave $Q,U$ data with the layered mixing matrix, or with a 3D-informed moment expansion if $K$ is large.
5. Propagate the full posterior into $C_\ell^{BB}$ and $r$, and publish a 3D-complexity map used as a weight.

That is the joint prediction of frequency-dependent $Q$ and $U$. Listing “3D maps + grains + field” without this mixing step does not yet specify a foreground model.

7. Limitations, residual systematics, and validation

This programme reduces a systematic to a larger statistical error. It does not remove dust.

Missing and misinterpreted data. Stellar polarimetry at high latitude is sparse; many CMB-deep fields have too few stars per beam at $p \sim 0.1\%$ to split two clouds. 3D extinction maps are smoothed relative to Planck’s $5'$ 353 GHz beam and relative to forthcoming CMB experiments, and they end at a finite distance. Dust–reddening conversion varies with $R_V$, grain evolution, and the ISRF; using a single factor is an unaccounted mixing-matrix error. H I misses H$_2$ and H$^+$ and assigns kinematic, not geometric, distances. Inclination $\gamma_k$ is not observed; it aliases into $p_k$. Grain optical constants at millimetre wavelengths remain a model choice; different successful extinction models predict different $\beta_d^P(\nu)$ and different $p(\nu)$.

Decorrelation is not one number. A pipeline that reports a single $\rho$ for the survey footprint can miss the Pelgrims-like tail of misaligned pixels and can overstate the problem if that tail is masked. Conversely, a detection of $\rho_\ell < 1$ on a large mask is not automatically LOS angle mixing.

Other foregrounds and instrument effects are not solved by 3D dust. Synchrotron, with a spatially varying $\beta_s$ and its own Faraday and LOS structure, must be modelled with low-frequency data (WMAP, Planck LFI, C-BASS, S-PASS, QUIJOTE). Lensing $B$ modes are a cosmological foreground to primordial tensors; they are delensed with internal CMB lenses or with large-scale structure tracers, not with dust maps. Beams, bandpasses, polarization-angle calibration, and cut-sky $E/B$ leakage produce spurious $B$ modes that can mimic a dust residual. CO and other line contamination in specific bands, and a possible magnetic-dipole dust component, add chromatic residuals that a two-cloud MBB does not capture.

Validation must be end-to-end. Minimum tests, all relative to a conventional one-MBB dust model:

• Simulations: inject a 3D dusty MHD (or layered) sky with known grain mix and field, observe it with the true bandpasses and noise, and recover $r$ with the hierarchical pipeline. Posterior coverage of $r$ is the figure of merit, not the residual map’s visual cleanliness.
• Held-out frequency: fit without 217 GHz (say) and predict $Q_{217},U_{217}$; the prediction error should match the posterior.
• Cross-spectra and split-map nulls: detector-set and half-mission crosses, as in Planck, to keep noise from masquerading as $\rho_\ell < 1$.
• Complexity-split tests: the improvement over one-MBB should concentrate on the H I-selected multi-cloud, misaligned sample; the control sample should be consistent with a single template.
• Transfer of 3D uncertainty: repeating the analysis on posterior samples of the extinction cube (Edenhofer et al. release 12 variational samples) must broaden the $r$ posterior. If it does not, the auxiliary data have been treated as truth.

The realistic scope is this. On a large fraction of the high-latitude sky a single MBB remains an adequate description at current sensitivity, which is why Planck Collaboration XI (2020) found no significant power-spectrum decorrelation. On a minority of sightlines, identifiable in 3D extinction, H I, and stellar polarization, a 353 GHz template is not predictive at 150 GHz, at a level that matters for $r \sim 10^{-3}$. A physically layered model of $Q$ and $U$, with grain-informed SED priors and distance-resolved angles, is the correct replacement for that template. Its output is a CMB $B$-mode likelihood that knows where the dust subtraction is trustworthy and where it is not.

Deviations from brief

None of substance. Two conventions are stated explicitly because they are easy to misread: $\gamma$ is the angle between $\mathbf{B}$ and the POS (polarized amplitude $\propto \cos^2\gamma$), not the inclination to the LOS; and the mixing coefficient $M_{\nu k}$ is a scalar per layer, with frequency-dependent $\psi_E$ arising only after the vector sum over layers. Planck Intermediate Results L is treated as a reported indication later questioned by the 2018 cross-spectrum analysis and by Sheehy & Slosar (2018), consistent with the brief’s empirical caution.