Layered dust / B-modes

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

drive · 29.6 KB
Fable 5.1driver Grok 4.6

Three-dimensional dust, grain physics, and the plane-of-sky field as a cure for spectral-index decorrelation in CMB B-mode searches

1. The observational problem

Primordial tensor B-modes peak at the recombination bump, ℓ ≈ 80, where D_ℓ^BB ≡ ℓ(ℓ+1)C_ℓ^BB/2π ≈ 0.06 r μK²; gravitational lensing adds a nearly white ≈ 2×10⁻³ μK² at the same multipoles. The current state of the art is BK18 (BICEP/Keck XIII, Ade et al. 2021): σ(r) = 0.009 and r_0.05 < 0.036 at 95%, tightened to r < 0.032 with Planck PR4 and BAO (Tristram et al. 2022). The next generation targets σ(r) ≈ 3×10⁻³ (Simons Observatory) and δr ≲ 10⁻³ (CMB-S4, LiteBIRD).

Polarized thermal dust is the obstacle at every frequency above ~70 GHz. Planck Int. XXX established that there is no dust-free window: in the BICEP2 field the 353 GHz B-mode power extrapolated to 150 GHz is 1.32×10⁻² μK² over 40 < ℓ < 120, the same magnitude as the excess BICEP2 had attributed to r = 0.2. BK18 measures A_d,353 = 4.4 (+0.8/−0.7) μK² at ℓ = 80 and β_d = 1.49 (+0.13/−0.12) in that field. Scaled with a modified blackbody (MBB) at T_d = 19.6 K, that is about 1×10⁻² μK² at 150 GHz, equivalent to r ≈ 0.15, and about 1.5×10⁻³ μK² at 95 GHz, equivalent to r ≈ 0.02–0.03. Even in the cleanest patch, the dust must be modeled to a small fraction of itself.

The sky-averaged polarized dust SED is deceptively simple. Planck 2018 XI finds a single MBB with β_d^P = 1.53 ± 0.02 at T_d = 19.6 K from 353 GHz down to below 70 GHz, β_d^P − β_d^I = 0.05 ± 0.03, power-law angular spectra with α_EE = −2.42 ± 0.02 and α_BB = −2.54 ± 0.02 on the largest region (LR71), BB/EE = 0.524 ± 0.005, positive TE and TB, and a polarization fraction p whose frequency dependence is "not strong" (Planck Int. XXII had earlier reported a (21 ± 6)% drop of p from 353 to 70 GHz). SPIDER (Ade et al. 2025) finds β_d = 1.45 ± 0.05 (E) and 1.47 ± 0.06 (B) over its region.

Why does the shape matter at the few-percent level when the amplitude is fitted from the data? Because any cleaning that scales a dust template, or fits one MBB per pixel, removes only the part of the dust that is correlated with that template across frequency. If R is the correlation coefficient between the dust maps at the template and target frequencies, the residual power after optimal scaling is (1 − R²) C_ℓ^dust(ν_target), and the induced bias is

Δr ≈ (1 − R²) · D_ℓ^dust(ν_target, ℓ≈80) / D_ℓ^BB(r=1, ℓ≈80)
   ≈ 0.15 (1 − R²)            for BK-like dust at 150 GHz

so R = 0.99 costs Δr ≈ 3×10⁻³ and R = 0.999 costs Δr ≈ 3×10⁻⁴. A global mis-scaling of the SED is harmless; it is absorbed by the fitted amplitude. What is not harmless is any part of the dust whose frequency scaling differs from the rest within the same pixel or the same patch. That is what "decorrelation of the spectral index" means, and it is a bias on r, not merely an inflated error bar.

2. Decorrelation, operationally

A single MBB per pixel is a rank-1 statement about the frequency dimension:

Q(ν, n̂) = f(ν; β_d(n̂), T_d(n̂)) · Q_353(n̂)
U(ν, n̂) = f(ν; β_d(n̂), T_d(n̂)) · U_353(n̂)
f(ν; β_d, T_d) = (ν/353 GHz)^{β_d} · B_ν(T_d)/B_353(T_d)   (times unit conversion)

with the same scalar f for Q and U, hence a frequency-independent angle ψ = ½ arctan(U/Q) and a frequency-independent p. Spatial variation of β_d(n̂) already breaks the map-level proportionality between frequencies ("spatial decorrelation"), but it is still rank-1 per pixel and is absorbed by a spectral-index map (Planck Commander products) or by a moment expansion (Chluba, Hill & Abitbol 2017; Mangilli et al. 2021). Line-of-sight (LOS) decorrelation is qualitatively different: within one beam, Q and U acquire different SEDs and no pixel-wise scalar f(ν) exists. Operationally one sees it at four levels:

  • Map level: Q_ν2/Q_ν1 is not a scalar and Q_ν2/Q_ν1 ≠ U_ν2/U_ν1.
  • Angle level: ψ(ν) rotates with frequency ("spectral rotation"), the observable proposed by Tassis & Pavlidou (2015) and detected by Pelgrims et al. (2021).
  • Spectrum level: R_ℓ(ν1,ν2) = C_ℓ^{ν1×ν2} / √[C_ℓ^{ν1} C_ℓ^{ν2}] < 1. If a fraction ε² of the power at ν2 is uncorrelated with the map at ν1, R = √(1 − ε²) ≈ 1 − ε²/2. Planck Int. L parametrized R^BB_{50–160} = 1 − K (N_HI/10²⁰ cm⁻²)^γ with K = 0.40 ± 0.32, γ = −1.76 ± 0.69, predicted R ≈ 0.85 for the BICEP2 column (N_HI ≈ 1.6×10²⁰ cm⁻²) and a bias of +0.018 on r if ignored. BK15/BK18 encode the same idea as a parameter Δ_d multiplying dust cross-spectra.
  • Pipeline level: template cleaning leaks (1 − R²) of the dust power; a parametric per-pixel MBB fit pushes the misfit into the CMB amplitude (McBride, Bull & Hensley 2023 find >10σ biases on CMB Q,U when SED and angle both vary along the LOS); an internal linear combination (ILC) with a rank-1 dust column leaves the extra rank in the "CMB" channel. In mixing-matrix language, d_ν = Σ_c A_νc s_c + n, and LOS-complex dust occupies as many columns as there are distinct SEDs along the LOS. Moment expansions add columns ∂f/∂β, ∂f/∂T, ∂²f/∂β² with free spin-2 maps; they de-bias at the price of σ(r) (Vacher et al. 2022 find ≈ 21 σ_r bias from spatially varying β_d, T_d in LiteBIRD simulations, recovered by β and T moments at σ(r) = 8.8×10⁻⁴). Ritacco et al. (2023) show on Planck HFI at ℓ = 4–32 that a first-order moment fit reproduces EE but not BB, and that the residuals at 100, 143 and 217 GHz are not fully correlated, which they attribute to frequency-dependent polarization angles.

The observational record looks contradictory but is not. Planck Int. L claimed R^BB(217×353) < 1 at >99% confidence; Sheehy & Slosar (2018) traced this to noise bias in split cross-spectra and to correlations between patches and found R^BB ≈ 1.001 (ℓ = 55–90, 72% of sky); Planck 2018 XI, fitting four HFI frequencies jointly, found no evidence for a loss of correlation and concluded decorrelation "might not be a problem" for experiments targeting r ≈ 0.01. BK18 finds no evidence for Δ_d ≠ 1, yet keeps it fixed in the baseline because it is "a biasing parameter with respect to r". BK XVI finds no decorrelation between HI-filament dust and the rest of the dust field down to 95 GHz. SPIDER finds no significant LOS decorrelation, but its two sub-regions give β_d = 1.52 ± 0.06 versus 1.09 ± 0.09 (3.9σ apart), and a joint SPIDER+Planck fit in one sub-region is inconsistent with a common-morphology MBB at >3σ. Pelgrims et al. (2021), conversely, selected Planck pixels with two or more HI clouds (Panopoulou & Lenz 2020) whose HI-derived B_POS orientations are misaligned, and detected a larger ψ(217) − ψ(353) than in a single-cloud control sample; the probability of a noise-only origin ranges from 5×10⁻² to 4×10⁻⁷ depending on the CMB subtraction, and a ~15% cloud-to-cloud variation in the 353/217 polarized ratio reproduces the effect.

The resolution is that R_ℓ over tens of percent of the sky is a weak statistic for a localized, second-order effect. A 15% SED contrast between two layers of comparable polarized brightness gives, as shown in Section 4, an uncorrelated amplitude fraction of order 5–8% in those pixels, i.e. R ≈ 0.997–0.999, which averages toward unity over regions dominated by single-cloud sightlines. That level is invisible to Planck's R_ℓ and irrelevant at r ≈ 0.01, exactly as Planck 2018 XI states, and it is a bias of Δr ≈ 3×10⁻⁴–10⁻³ in a 150 GHz analysis of a BK-like field, exactly what a δr ≈ 10⁻³ experiment cannot tolerate. The question is therefore not whether decorrelation exists but how to predict where and how much, so it can be modeled instead of marginalized.

3. Microphysics: a polarized MBB is already an average over n(a)

The emissivity of a dust population is a sum over sizes and materials,

j_ν = ∫ da n(a) C_abs(a, ν) B_ν(T_d(a))

In the Rayleigh limit (2πa/λ ≪ 1, satisfied by a ≲ 1 μm at λ ≥ 0.85 mm by three orders of magnitude) a spheroidal grain has C_abs,j = k V Im[(ε − 1)/(1 + L_j(ε − 1))] for the field along axis j, with L_j the geometric depolarization factors. Two consequences follow. First, C_abs ∝ V ν · Im[…], so the submillimetre opacity per unit mass is independent of grain size; the spectral index β_d = 1 + d ln Im[…]/d ln ν is a property of the dielectric function ε(ν) of the material (silicate, carbon, or the composite "astrodust" of Draine & Hensley 2021), typically 1.5 ≲ β ≲ 2, with laboratory amorphous solids showing a temperature-dependent flattening toward millimetre wavelengths (two-level-system models, Meny et al. 2007). Second, the intrinsic polarization of a perfectly aligned grain, (C_⊥ − C_∥)/(C_⊥ + C_∥), is also size-independent and only weakly frequency-dependent for a single material. A flat p(ν) therefore points to a single composition among the aligned grains: Hensley & Draine (2023) argue that a one-composition model for grains larger than ≈ 0.02 μm most naturally explains Planck's weak frequency dependence of the far-infrared polarization fraction, whereas separately aligned silicate and unaligned carbon populations with different β and T (as in Draine & Fraisse 2009) predict a p(ν) that varies with frequency.

Temperature enters through radiative equilibrium with the interstellar radiation field of strength U:

∫ C_abs(a, ν) 4π J_ν(U) dν = 4π ∫ C_abs(a, ν) B_ν(T_d) dν

For a ≳ 0.1 μm the UV/optical absorption is geometric (∝ a²) while the FIR emission is ∝ a³ T^{4+β}, giving T_d ∝ (U/a)^{1/(4+β)}: a factor 2 in U or in a changes T_d by only 11–12%. Grains with a ≲ 0.01 μm are stochastically heated and have no T_d; they dominate λ ≲ 100 μm and the anomalous microwave emission below 60 GHz, are essentially unpolarized (PAHs), and are irrelevant to the polarized 100–353 GHz SED. A one-temperature MBB fitted to j_ν therefore returns an effective (β_eff, T_eff) that is a moment average over the T(a) and ε(ν) mixture, which is why fitted β and T are anticorrelated and why β_eff is generally below the material β.

So where does the grain size distribution actually enter the polarized SED? Through three channels. (i) Alignment is size-selective: radiative torques (Lazarian & Hoang 2007; Andersson et al. 2015) align grains with a ≳ 0.05–0.1 μm, so the polarized emission weights the mass in the largest grains, which are also the coolest; this is why β_d^P and β_d^I need not coincide, and the measured β_d^P − β_d^I = 0.05 ± 0.03 (Planck 2018 XI; 0.08 in Planck Int. XXII) is a direct measure of how different the aligned sub-population is. (ii) n(a) evolves with environment, and the evolution changes the material: coagulation and mantle accretion in the cold neutral and molecular medium produce aggregates with a lower β (by a few tenths) and a submillimetre emissivity per H higher by factors 2–3 (THEMIS: Jones et al. 2013, 2017; Köhler, Ysard & Jones 2015), consistent with the observed rise of τ_353/N_H and fall of T_d from diffuse to translucent gas (Planck 2013 XI; Planck Int. XVII); shattering in shocks replenishes small grains. (iii) The alignment efficiency itself depends on U/n_H, so p_i per cloud varies with 3D position and shielding.

For the CMB bands the relevant contrast is the ratio of two MBBs between 353 GHz and 150 GHz (or 217 GHz). Since hν/kT_d ≈ 0.4–0.9 here, this ratio is dominated by β:

d ln[S(353)/S(150)]/dβ = ln(353/150) = 0.86      Δβ = 0.2  gives 17%
d ln[S(353)/S(150)]/dT_d ≈ 0.015 K⁻¹ at 19.6 K   ΔT = 2 K  gives 3%

(d ln[S(353)/S(217)]/dβ = 0.49, d/dT_d ≈ 0.010 K⁻¹). Cloud-to-cloud β differences of 0.1–0.3, i.e. exactly the diffuse-versus-translucent evolutionary contrast, produce the ~15% polarized-ratio variation Pelgrims et al. (2021) need; temperature differences alone would require ΔT ≳ 10 K. Temperature, by contrast, controls the 545/857 GHz and IRAS bands that 3D emission maps fit. This division of labor between β (grain material and evolutionary state) and T_d (radiation field) is what allows the three ingredients to be combined rather than merely listed.

4. Two clouds: the Tassis–Pavlidou mechanism

Aligned oblate grains emit with the electric vector perpendicular to B_POS. With γ the angle between B and the LOS (Planck papers use the complementary angle to the plane of the sky), the standard integrals (Lee & Draine 1985; Fiege & Pudritz 2000; Planck Int. XX) are

I = ∫ S_ν [1 − p_0 (sin²γ − 2/3)] dτ
Q = ∫ p_0 S_ν sin²γ cos 2ψ dτ,    U = ∫ p_0 S_ν sin²γ sin 2ψ dτ,    ψ = ψ_B_POS + π/2

Collapse the LOS into two clouds with S_i(ν) = τ_i (ν/ν_0)^{β_i} B_ν(T_i) and polarized brightness P_i(ν) = p_i S_i(ν), p_i = p_0 f_align,i sin²γ_i:

Q(ν) + iU(ν) = P_1(ν) e^{2iψ_1} + P_2(ν) e^{2iψ_2}
             = P_1(ν) e^{2iψ_1} [1 + ρ(ν) e^{2iΔψ}],    ρ ≡ P_2/P_1,  Δψ ≡ ψ_2 − ψ_1

so that

P(ν)² = P_1² [1 + ρ² + 2ρ cos 2Δψ]
tan 2[ψ(ν) − ψ_1] = ρ sin 2Δψ / (1 + ρ cos 2Δψ)
dψ/d ln ν = ½ · sin 2Δψ · ρ/(1 + ρ² + 2ρ cos 2Δψ) · d ln ρ/d ln ν
d ln ρ/d ln ν = (β_2 − β_1) − [g(x_2) − g(x_1)],   g(x) = x eˣ/(eˣ − 1),  x_i = hν/kT_i

The derivative shows the two necessary conditions explicitly. If Δψ = 0 (aligned fields), sin 2Δψ = 0: different SEDs distort the total and polarized SEDs (the β and T moments of the LOS), but ψ is frequency-independent and the sky is still rank-1 in Q,U, curable by a spectral-index map or moments. If d ln ρ/d ln ν = 0 (identical SEDs), misaligned fields depolarize (P < P_1 + P_2) but the depolarization is the same at every frequency, so nothing decorrelates. Only the product of an SED contrast (set by n(a), material, U) and a field misalignment (set by B_POS geometry) rotates ψ with frequency. This is Tassis & Pavlidou's statement that the 350 GHz map is then not predictive of the 150 GHz map.

The consequence for cleaning follows from the same algebra. Scale the 353 GHz template with cloud 1's SED s(ν) = S_1(ν)/S_1(353); the residual at ν is

Δ(Q + iU) = s(ν) P_2(353) [ρ(ν)/ρ(353) − 1] e^{2iψ_2} ≈ s P_2(353) Δln ρ · e^{2iψ_2}

Its component along the template direction is a change of effective SED and is absorbed by a refitted amplitude; its component perpendicular to the template, s P_2 Δln ρ sin 2(ψ_2 − ψ_T), is the genuinely uncorrelated part. For two clouds of comparable polarized brightness with Δψ ≈ 45° and Δln ρ = 0.15, the uncorrelated amplitude is ≈ 5–8% of the total, 1 − R² ≈ 0.3–0.6%, and the angle rotates by ≈ 2° between 353 and 217 GHz. Both numbers match the Pelgrims et al. detection and the Planck 2018 XI non-detection simultaneously, and both are what a δr ≈ 10⁻³ analysis must remove.

5. Ingredient 1, alone: 3D maps of thermal dust emission

A 2D 353 GHz map cannot distinguish one cloud with p = 10% from two clouds with p = 15% partially cancelling. This is the where degeneracy, and 3D maps break it.

Three data classes exist. First, 3D reddening maps from stellar photometry and Gaia parallaxes: Bayestar19 (Green et al. 2019; 3π sr), Lallement et al. (2019, 2022) and Vergely et al. (2022) from Gaia+2MASS, Leike, Glatzle & Enßlin (2020; 400 pc at 1 pc resolution), and Edenhofer et al. (2024; 54 million Gaia BP/RP extinctions from Zhang, Green & Rix 2023, 14′ angular resolution, 516 distance bins from 69 pc to 1.25 kpc). These give the differential extinction density, i.e. dust column per distance bin, and they show that at |b| ≳ 30° most dust sits within ~300 pc in the Local Bubble wall plus a few nearby clouds. They are not emission maps: they know nothing about T_d or β. Second, 3D emission maps: Zelko & Finkbeiner (2022) fitted, for 500 million Bayestar19 voxels, a temperature, an emissivity slope, and a τ/E(B−V) conversion factor, constrained so that the LOS-integrated model matches five Planck and IRAS bands at 10′; the first 3π map of T_d in 3D, delivered at 110′, 55′ and 27′. They find the τ/E(B−V) factor varies by a factor of two, and they recover two distinct temperatures along the two-component Cepheus sightline at about 1σ. Third, HI as the third axis: HI4PI velocity structure decomposed into clouds (Panopoulou & Lenz 2020: on average 3 clouds per 0.2 deg² pixel in the northern cap and 2.5 in the south, including the BICEP/Keck field), and velocity-resolved HI filament orientation (Clark & Hensley 2019, whose HI-based Q,U correlate at r > 0.75 with Planck 353 GHz; used for the southern and northern caps by Ghosh et al. 2017 and Kim, Choi & Flauger 2019, and in BK XVI, where HI templates correlate with BK B-modes at the 10–65% level over 20 < ℓ < 200).

Alone, these maps (a) classify each pixel as single- or multi-layer, which is precisely the selection that made the Pelgrims detection possible; (b) provide the emission weights τ_i (or N_d,i) and, from Zelko-type fits, T_i and a Wien-side β_i per layer; (c) supply LOS statistics: Halal, Clark & Tahani (2024) show p_353 decreases when the column is split across several distance components, a 3D signature of LOS depolarization; Clark (2018) shows HI velocity coherence is linearly correlated with p_353. What they cannot do is predict Q and U: without ψ_i and f_align,i a 3D emission model is a model of I(ν) only.

6. Ingredient 2, alone: the grain size distribution

Section 3 showed that β_i, T_i, and p_i are not free parameters but functions of n(a), material, U, and phase. This is the what spectrum degeneracy: two Planck frequencies constrain only the emission-weighted β of a LOS, not the SED of each layer. Grain physics converts a per-layer dust column and environment into a per-layer SED prior:

  • diffuse atomic gas, U ≈ 1: astrodust-like β ≈ 1.5, T_d ≈ 19–20 K, τ_353/N_H ≈ 7×10⁻²⁷ cm²;
  • translucent/CNM: β lower by 0.1–0.3, T_d ≈ 16–18 K, emissivity per H up by 2–3, alignment reduced in shielded gas;
  • radiation field: T_i ∝ U(x_i)^{1/(4+β)} with U(x) from the 3D distribution of OB stars, so a cloud near Sco–Cen and one behind it differ in T_d predictably.

Observable proxies of the aligned-grain size distribution exist per cloud: the wavelength λ_max of the Serkowski curve (Kim & Martin 1995; Draine & Fraisse 2009) from multi-band stellar polarimetry, and R_V from Gaia XP extinction curves, which now exists in 3D (Zhang & Green 2025) and can be attached to the same distance bins as the dust density. Alone, a Galactic-mean n(a) (MRN, Weingartner & Draine 2001, Hensley & Draine 2023) fixes a mean SED and reproduces the sky-averaged β_d^P = 1.53; it says nothing about the cloud-to-cloud contrast that drives decorrelation. Using one n(a) for every layer would erase exactly the diversity the method needs.

7. Ingredient 3, alone: the sky-projected magnetic field

B_POS fixes ψ_i (emission perpendicular to B_POS) and, with γ_i and the turbulent dispersion, the depolarization within and between layers. This is the which Stokes mixing degeneracy: the integrated 353 GHz angle is a vector sum and cannot be decomposed into ψ_i from one frequency.

Data: Planck 353 GHz gives the integral Σ_i P_i e^{2iψ_i}. Starlight polarization traces the same aligned large grains in absorption, parallel to B_POS; Planck Int. XXI showed the emission-to-extinction ratio R_P/p = P_353/p_V ≈ 5.4 MJy sr⁻¹ and consistent angles. With Gaia distances, stars at increasing distance sample cumulative layers: Panopoulou et al. (2019) demonstrated two-cloud tomography toward Ursa Major, Pelgrims et al. (2023) formalized a Bayesian LOS inversion (BISP-1) that recovers per-cloud distance, mean (q,u), and intrinsic scatter, and Pelgrims et al. (2024) produced the first degree-scale (4 deg²) tomographic map with RoboPol. PASIPHAE (Tassis et al. 2018) plans ~10⁴ deg² over both polar caps at ~360 stars/deg² with two WALOP polarimeters (Skinakas, SAAO). HI filaments provide a statistical B_POS per velocity slab (Clark et al. 2014, 2015; Planck Int. XXXII, XXXV; Clark & Hensley 2019), which is how Pelgrims et al. (2021) obtained per-cloud orientations. Faraday rotation and Faraday tomography measure n_e B_LOS; they constrain γ and are valuable 3D magnetometry but they are not B_POS and must not be used as if they were. Alone, an integrated B_POS reproduces the 353 GHz angle by construction and tells nothing about how it will rotate at 150 GHz.

8. The joint method: a layer-by-layer polarized sky

Each ingredient removes one degeneracy: 3D emission maps break the where, grain physics breaks the what spectrum, B_POS tomography breaks the which Stokes mixing. Combined, they turn an under-constrained extra-rank problem into a constrained multi-layer forward model:

  1. Segment the LOS. For each CMB pixel, define layers i = 1…N from 3D dust density (Edenhofer/Bayestar distance bins) cross-matched with HI velocity components (Panopoulou & Lenz; Clark & Hensley), with column N_d,i, distance d_i, and phase indicators (HI linewidth, column, CO where present). At high latitude N ≤ 3 and everything is within ~500 pc.

  2. Assign an SED per layer. β_i, T_i, ε_i (emissivity per H) drawn from priors conditioned on phase, U(d_i) from a 3D radiation-field model, and size proxies (R_V, λ_max); Zelko-type 3D T_d fits and Planck 545/857 GHz enter as data on T_i. The astrodust+PAH model (Hensley & Draine 2023) or THEMIS supplies the mapping from n(a) and environment to (β, T, f_align), and the constraint β_d^P − β_d^I ≈ 0.05 fixes how the aligned sub-population differs from the mean.

  3. Assign ψ_i and p_i per layer. From starlight tomography, ψ_i = ψ_star,i + 90° and P_i(353) ≈ R_P/p · p_V,i, with R_P/p corrected for T_i and β_i; from HI filament orientation where stars are too sparse. Planck 353 (and 217, 143) Q,U enter as the integral constraint Σ_i P_i(ν) e^{2iψ_i} = (Q + iU)_ν, which pins the residual freedom in p_i sin²γ_i.

  4. Forward-model Q(ν), U(ν) at 90/150/220/280 GHz and across LiteBIRD's 40–402 GHz bands with bandpass integration. The Section 4 algebra shows that the model predicts, pixel by pixel, both the rotation ψ(ν) and the spin-2 "decorrelation template" s(ν) P_2 Δln ρ e^{2iψ_2}.

  5. Use it in cleaning in one of four ways, ordered by how much of the model one trusts:

    • as N polarized templates with tightly bounded (β_i, T_i) and fixed ψ_i, i.e. a layered Commander or a layered BK-style likelihood with amplitudes A_d,i;
    • as a spatially varying mixing matrix in a constrained ILC that deprojects the predicted decorrelation templates rather than generic ∂f/∂β maps (cf. the moment-ILC of Remazeilles et al. 2021);
    • as a subtracted sky whose residual complexity is absorbed by a low-order spin-moment expansion (Vacher et al. 2022, 2023) with moment maps predicted from the layers, W_β(n̂) = Σ_i (β_i − β̄) P_i e^{2iψ_i}, so only a few amplitudes are fitted and the σ(r) penalty of free moment maps is avoided;
    • as a weight or mask: pixels with large posterior variance of the predicted ψ(ν) (multi-cloud, misaligned, comparable P_i) are downweighted, the inverse of the Pelgrims selection.

  6. Validate. The model must predict where ψ(217) − ψ(353) is large in Planck (the Pelgrims target sample), must not overpredict decorrelation in SPIDER-like or LR71-like averages (R^BB ≈ 1 to within ~1%), must reproduce the BK XVI HI-dust correlation down to 95 GHz, and must reproduce the sign and amplitude of dust TB and EB from filament–field misalignment (Clark et al. 2021), including their predicted frequency dependence (Vacher et al. 2023).

  7. Propagate. Sample the 3D-map posteriors (Edenhofer et al. provide samples), the BISP-1 per-cloud posteriors, and the grain-model hyperparameters into an ensemble of dust skies; marginalize over it in the r likelihood. The residual bias is then set by the posterior width of Δln ρ_i and ψ_i, not by the unknown extra rank. Reducing the unmodeled SED contrast from 15% to 3% lowers the uncorrelated power by a factor 25, from Δr ≈ 10⁻³ to ≈ 4×10⁻⁵ in the Section 4 example.

9. Limitations and residual systematics

  • High-latitude stellar sparsity and polarimetric depth. At E(B−V) ≈ 0.02 and p_V/E(B−V) ≲ 13% (Panopoulou et al. 2019b), stellar polarizations are ≲ 0.2–0.3%; BISP-1 recovers clouds only when the signal exceeds ≈ 0.1%. PASIPHAE's degree-scale maps do not yet exist over CMB patches, and separating clouds 50 pc apart at 300 pc needs parallax precision at the edge of Gaia for the faint stars that dominate the counts.
  • Emission versus extinction weighting. Optical polarization weights layers by column and f_align; emission weights them by τ_i T_i^{4+β}-like factors; R_P/p scatters between sightlines (Planck Int. XXI, Planck 2018 XII). The relation between the two is a fixed 90° rotation, not an ambiguity, but the sign convention of U (IAU versus HEALPix/COSMO) is a genuine trap.
  • Inclination and turbulence. p_i degenerates between f_align,i and sin²γ_i; γ_i needs B_LOS (pulsar rotation measures, Zeeman) or statistical priors (Planck Int. XLIV: mean plus turbulent field, f_M ≈ 0.9, effective N ≈ 7 layers). Sub-layer dispersion of ψ within a distance bin is unresolved and depolarizes each layer differently.
  • Resolution and association. Zelko-type T_d is at 27′–110′; Edenhofer's 14′ is fine for ℓ ≲ 300 but distance bins are logarithmic; HI velocity is not distance, dark H₂ and ionized gas carry dust invisible in HI, and the N_HI-to-reddening ratio drifts above N_HI ≈ 4×10²⁰ cm⁻² (Lenz, Hensley & Doré 2017).
  • Grain-model dependence. Astrodust versus THEMIS give different β(ν) curvature toward millimetre wavelengths; a Wien-side β from 545/857 GHz need not equal the β at 100–353 GHz; stochastic heating contaminates IRAS 100 μm in 3D T_d fits; the T_d–β degeneracy (ΔT = 1 K trades against Δβ ≈ 0.02 at 217–353 GHz) persists per layer.
  • Contaminants that mimic rotation. CIB in intensity biases τ_i; residual CMB E-modes at 143/217 GHz rotate the measured dust angle (Pelgrims et al.'s significance range is set by the CMB subtraction); polarization-efficiency and bandpass-mismatch leakage, and absolute polarization-angle calibration errors between channels, are instrumental ψ(ν) rotations; synchrotron below ~100 GHz and weakly polarized AME below 60 GHz must be separated, not folded into β_d.
  • False friends. Dust EB from filament–field misalignment is a birefringence-like foreground with its own frequency dependence, not a β_d effect; moment expansions can de-bias without any 3D information but pay in σ(r); PySM d6 (Vansyngel-type frequency covariance with a single correlation length) and d12 (six-layer Martínez-Solaeche et al. 2018 model) are phenomenological stand-ins, not measurements. Planck's near-constant p from 353 to 70 GHz bounds the allowed SED diversity: the target is 10%-level layer contrasts, not order-unity redesigns of dust.

The three ingredients do not "solve" foregrounds. They convert an unknown-rank frequency mixing into a fit with a known number of layers, physically bounded SEDs, and externally measured angles, whose residual can be budgeted. That budget, with its stellar-density, γ, turbulence, calibration, and grain-model terms, is what a δr ≈ 10⁻³ analysis needs and what a 2D 353 GHz template cannot provide.

Deviations from brief

  • BICEP/Keck dust level. The brief states the BK field's dust BB is "comparable to r ~ 0.02 at 150 GHz". With A_d,353 = 4.4 μK² (BK18), the 150 GHz dust power is ≈ 1×10⁻² μK², equivalent to r ≈ 0.15; the r ≈ 0.02–0.03 equivalence holds at 95 GHz. Planck Int. XXX's 150 GHz value (1.32×10⁻² μK²) confirms the higher number. Corrected in Section 1.
  • "Large grains have greyer spectra". For interstellar sizes (a ≲ 1 μm) at λ ≥ 0.85 mm, 2πa/λ ≲ 0.007, so the Rayleigh limit holds and β is set by the dielectric function, not by size. The flattening of β with grain growth arises from the material and structure of aggregates and mantles, and from the temperature dependence of amorphous-solid absorption, as stated in Sections 3 and 6.
  • "180° ambiguity" between optical and submillimetre polarization. Q,U depend on 2ψ, so a 180° change of ψ is not an ambiguity. The relation is a deterministic 90° rotation; the practical pitfalls are the U sign convention and the different LOS weighting, which Section 9 lists instead.
  • Clarification, not correction. The brief leaves open whether T_d or β_d contrasts drive the SED ratio between clouds. At 100–353 GHz the ratio is dominated by Δβ (Δβ = 0.2 gives 17% between 353 and 150 GHz; ΔT = 2 K gives 3%), so the grain-evolution channel matters more than the radiation-field channel for decorrelation, while T_d is what 3D emission maps constrain best. Also, Zelko & Finkbeiner (2022) is a two-author paper, and PySM d12 (six-layer 3D model) is the community's layered stand-in alongside d6.