3b · Mie, Rayleigh, and the Phase Matrices

For: users adding aerosols to their scene; method developers extending the scattering codebase; retrieval developers reasoning about polarization sensitivity.

Prev: 3a · Gas Absorption · Next: 3c · Mixing & δ-M Truncation

This page covers three of the four layer-optics arrays plus the aerosol contribution to the optical depth: the Fourier-moment phase matrices and , the aerosol single-scattering albedo , and the extinction cross-section (which maps microphysics to ).

What this page produces

For each aerosol type and each Fourier moment :

Symbol	Meaning	Used for
,	phase matrix Fourier moments	layer scattering
	aerosol single-scattering albedo	layer SSA
	extinction cross-section	aerosol optical depth from microphysics
	δ-M forward-truncation factor	consumed in Concepts/03c

For Rayleigh (and Cabannes), the same are computed from the Rayleigh / Cabannes Greek source.

The Mie pipeline

   Aerosol(size_dist, n_r, n_i)
        │
        ▼
    MieModel
        │
        ▼
   compute_aerosol_optical_properties
        │
        ▼
   AerosolOptics(GreekCoefs, ω̃, k, fᵗ)
        │
        ▼
   compute_Z_moments  (per Fourier moment m)
        │
        ▼
   (Z⁺⁺, Z⁻⁺) per (μ', μ, m)

The user-facing entry is compute_aerosol_optical_properties(model::MieModel) in src/Scattering/. The dispatcher branches on the decomposition algorithm (NAI2() or PCW()) carried in MieModel.

The Greek matrix B_l

Polarized Mie scattering is fully captured by the Greek expansion of the   phase matrix. Sanghavi & Stephens (2015, Eq. 22) print it in Siewert's form, which the vSmartMOM paper (Sanghavi et al. 2014) uses by reference:

Relative to that (Siewert/VLIDORT) convention the code's Stokes vector is    : construct_B_matrix stores and with the opposite sign and writes at , so its block is numerically the one above and only the elements differ. The conventions page gives the exact mapping, and the import/export checklist negates γ and ε for coefficients in Siewert's (2000) or VLIDORT's convention.

The six coefficients fully determine the polarized scattering pattern — is the phase function expansion (the only one a scalar code uses); the other five carry polarization coupling. They live in the GreekCoefs struct at src/Scattering/types.jl:462–475:

struct GreekCoefs{FT}
    α::Vector{FT}
    β::Vector{FT}
    γ::Vector{FT}
    δ::Vector{FT}
    ϵ::Vector{FT}
    ζ::Vector{FT}
end

The output of compute_aerosol_optical_properties packages this with the scalar bookkeeping (src/Scattering/types.jl:512–533):

struct AerosolOptics{FT}
    greek_coefs::GreekCoefs
    ω̃::Union{FT, Array{FT}}      # single-scattering albedo
    k::Union{FT, Array{FT}}      # extinction cross-section
    fᵗ::Union{FT, Array{FT}}     # δ-M truncation factor
end

Effective aerosol Fourier support

The maximum particle size parameter determines a conservative Mie-series length, but the extreme-radius tail of an integrated size distribution can carry negligible weight. The optional radiative_transfer.greek_beta_cutoff therefore applies a second test after size-distribution integration:

This is evaluated separately for each aerosol a; their maximum contributes to the band-level m_max alongside Rayleigh, surface, and source support. The absolute value is intentional because β can be signed. Only β is tested: the other five Greek families can be zero or change sign and are not used as the scalar phase-function support benchmark. The last passing degree is retained, so the implementation does not assume monotonic decay.

The maximum-size result remains the allocation ceiling, and stream/user caps are still applied after component aggregation. This criterion shortens the Fourier loop but does not itself reduce nstreams or angular-operator size.

Phase matrix Z from Greek + μ

Once the Greek coefficients are in hand, the per-Fourier-moment phase matrices come from , the generalized spherical-function matrices (Sanghavi 2014 Mie-Fourier Eq. 12–14):

The coded entry point is compute_Z_moments(pol_type, μ, greek, m, arr_type=arr_type) in src/Scattering/compute_Z_matrices.jl. It returns of shape (NquadN, NquadN) (or the pre-broadcast scalar-spectral version that gets expanded across nSpec in expandOpticalProperties).

NAI-2 vs PCW — two ways to compute the Greek expansion

Sanghavi 2014 (the Mie-Fourier paper) compares two decomposition methods:

• NAI-2 (Numerical Angular Integration, version 2; Siewert's formalism): computes the scattering matrix on a Gauss-Legendre quadrature in scattering cosine , then projects to Greek coefficients via integrals against associated Legendre functions. The   quadrature roots needed for accurate integration are determined upfront. Default path; entry at src/Scattering/compute_NAI2.jl:82.

• PCW (Pre-computed Wigner; Domke's formalism corrected): computes Greek coefficients directly from products of size-averaged Mie coefficients etc. weighted by precomputed Wigner-3j symbols. No reconstruction of the angular scattering matrix. Entry at src/Scattering/compute_PCW.jl:37.

Tradeoffs (Sanghavi 2014 Mie-Fourier §5):

Regime	Faster method
Single particle,    (size parameter > 228)	NAI-2 (no angular reconstruction needed)
Polydispersion over a finite size range up to ~80 μm	PCW (6–8× faster than NAI-2)
Single small particle	PCW (precomputed Wigner-3j makes individual particles cheap)

NAI-2 is the safer default. PCW pays off for big polydisperse populations that dominate aerosol retrievals over ocean. Selection is via the DECOMP_MAP in src/IO/Parameters.jl ("NAI2" or "PCW").

Rayleigh and the Cabannes choice

Pure-gas no-redistribution Rayleigh data

New science workflows select the calculation explicitly. Two methods are available, selected by dispatch on an AbstractMolecularRayleighMethod:

Method	Physics	Species
UsingRefInd()	refractive index plus (effective) King factor	He, Ar, N₂, O₂, CO₂
UsingPolTen()	molecular polarizability tensor (mean polarizability , anisotropy ); the tensor dispersion is a dispatch parameter	N₂, O₂, H₂

method = UsingRefInd()            # RefractiveIndexRayleigh() is an alias
p = molecular_rayleigh_properties(method, :CO2, wavelength_nm)
ratio = molecular_rayleigh_cross_section_ratio(method, :CO2, λ, λ_ref)
source = molecular_rayleigh_provenance(method, :CO2)

q = molecular_rayleigh_properties(UsingPolTen(), :N2, wavelength_nm;
                                  temperature_K = 288.15)
r = molecular_rayleigh_properties(UsingPolTen(ConstantAnisotropy()), :N2,
                                  wavelength_nm)   # Raman-module tensor

Both return the same named-tuple fields: species, wavelength_nm, refractive_index, king_factor, depolarization, and cross_section_cm2. The two-argument forms molecular_rayleigh_properties(species, λ) and molecular_rayleigh_cross_section_ratio(species, λ, λ_ref) remain the refractive-index method. RefractiveIndexRayleigh is kept as an alias of UsingRefInd, so existing calling code is unaffected.

Refractive index plus King factor (UsingRefInd)

The RI helper returns a refractive index, King factor, solver depolarization ratio, and scattering cross section in cm² molecule⁻¹ for He, Ar, N2, O2, and CO2. Its Ar and CO₂ coefficients were reviewed for the Habitable Worlds Observatory (HWO) exo-Earth studies and carry a machine-readable provenance record; He, N₂ and O₂ have no such record yet. The method argument is intentionally visible so an RI calculation cannot be confused with the polarizability-tensor implementation.

The cross-section-ratio helper is useful for a matched-reference-optical-depth experiment,

The cross section follows the conventional refractive-index/King-factor expression used by the cited coefficient compilations:

Here is exactly the depolarization convention consumed by get_greek_rayleigh. He and Ar are monatomic, so   and  ; they still have different cross-section dispersion. These helpers deliberately do not infer pressure, mean molecular mass, or layer column. They also do not turn on Raman redistribution. A caller comparing gases at fixed must set the optical-depth profile explicitly and build the Greek phase source from the returned depolarization at each wavelength (or place wavelengths with distinct depolarization in separate bands).

For Ar and CO₂, the coefficient provenance and distinct evidence intervals are machine-readable through molecular_rayleigh_provenance:

• The Peck–Fisher Ar measurement spans 17 vacuum wavelengths from 467.9465 to 2058.6909 nm. Sneep and Ubachs adopt its dispersion formula over 5000–33000 cm⁻¹ (303.0303–2000 nm). Ar has the monatomic identity  ; it has no empirical King-factor interval.

• The Bideau-Mehu CO₂ measurement spans 180.7–253.7 nm. The implementation is the corrected/scaled formula recorded by Sneep and Ubachs. Its broader 180.7–1694.5 nm envelope is a secondary database-declared formula range, not the primary measurement range. The source record is frozen at refractiveindex.info database commit 37773dce8128ccabcd9d5108d9958d6041244661, path database/data/main/CO2/nk/Bideau-Mehu.yml, SHA-256 d3b52d2a0b667a2e95aec110f35a7fced0c2ab3f02583d2e604aeb58c204f55c.

• The CO₂ King fit in Sneep and Ubachs is based on Alms–Burnham–Flygare measurements at 457.9, 488.0, 514.5, 568.2, and 647.1 nm. Therefore 457.9–647.1 nm is its empirical interpolation span; use outside that span is King-factor extrapolation, independently of whether the dispersion formula is in range.

• Wilmouth and Sayres (2019) directly check both gases continuously over 333–363 nm. Thalman et al. (2014), with its erratum, check Ar in selected regions between 350 and 660 nm, with corrected tabulated points at 370.0, 405.8, 532.2, and 660.0 nm. Those selected regions are not a continuous 350–660 nm validation interval. He et al. (2021) provide a continuous 307–725 nm direct check for CO₂.

The primary-measurement range, adopted dispersion envelope, King-factor support, and direct-scattering checks are not interchangeable; dispersion and King-factor extrapolation should be reported separately.

The N₂ and O₂ refractivities each join two dispersion fits at the wavenumber where they agree, so   and have no jump:

Gas	Short wavelengths	Long wavelengths	Join
N₂	Bates (1984), 288.15 K, fitted 254–468 nm (He et al. 2021, Table 1)	Peck and Khanna (1966), 0 °C, measured 467.9–2058.7 nm	436.309 nm
O₂	Bates (1984), 273.15 K, fitted 288–546 nm (He et al. 2021, Table 1)	Egan (2026), Eq. (8), 20 °C and 100 kPa, constrained at 0.254–1.542 µm	466.678 nm

Joining at the nominal limits instead would leave jumps of 2.8×10⁻⁴ (N₂, 468 nm) and 2.1×10⁻³ (O₂, 546 nm) in . Compared with the previous single fits extrapolated into the NIR, rises by 0.08 % (550 nm) to 0.21 % (2060 nm) and falls by 0.22 % to 0.89 %; for air the two nearly cancel (< 0.03 %). The King factors are those of He et al. (2021): N₂ from Bates (1984), and O₂ with the constant 1.09. Bates (1984) and Bodhaine et al. (1999) use 1.096 for O₂, which gives a 0.54 % larger O₂ cross section. Every fit is used with the ideal-gas number density of its reference state; gas compressibility (5–10×10⁻⁴ in   per molecule) is neglected, as in the cited compilations. Details and numbers are in docs/dev_notes/molecular_rayleigh_dispersion.md.

This API performs no explicit frequency redistribution. Its conventional Rayleigh cross section and effective King factor are not a separately computed Cabannes-only cross section. The RI method neither calls nor modifies the tensor/Raman machinery in InelasticScattering.

Polarizability tensor (UsingPolTen)

UsingPolTen(anisotropy) builds the Rayleigh properties from the mean polarizability and anisotropy of N₂, O₂ and H₂, through the Raman-module routines (getMolecularConstants, compute_effective_coefficents!, compute_σ_Rayl_coeff!). The anisotropy argument selects the wavelength dependence of the tensor:

• DispersiveAnisotropy() (default, UsingPolTen()): for N₂ and O₂, from the UsingRefInd refractivity through Lorentz–Lorenz,    , and   from its King factor ( at the fit's reference temperature), i.e. the anisotropy implied by the measured depolarization ratios behind the King factors (Bates 1984). The Buldakov temperature factor    scales relative to the reference temperature of the dispersion fit. At that temperature the cross section and depolarization equal the UsingRefInd values; between 200 and 300 K they differ by less than 4×10⁻⁴. H₂ has no UsingRefInd fit and keeps the Raman-module tensor.

• ConstantAnisotropy(): the Raman-module tensor, the Buldakov et al. single-pole dispersion and temperature dependence      and a wavelength- and temperature-independent anisotropy  . These constants produced the published Raman tables and remain the default of the Raman and Cabannes cross sections. Their equals the dispersive one near 488 nm (within 0.03 % for N₂ and 0.32 % for O₂); it is 2 % (N₂) and 4 % (O₂) low at 400 nm and 4 % and 8 % high at 2060 nm.

• Derived quantities:    (the same convention as the RI method),   , and  .

• refractive_index is the Lorentz–Lorenz equivalent at number_density_cm3 (default 2.546899×10¹⁹ cm⁻³). Pass number_density_cm3 = nothing to get refractive_index = nothing.

• temperature_K defaults to 288.15 K, the temperature of that default density and of the standard state of the He, Ar, N₂ (short-wavelength) and CO₂ RI fits and of the Bates (1984) and Bodhaine et al. (1999) cross sections, so the tensor is compared with those data at their own state. It enters only through the small temperature factor on .

Other species raise an ArgumentError that suggests UsingRefInd(). molecular_rayleigh_provenance is currently defined only for the RI method. The H₂ tensor is checked against the Peck and Huang (1977) refractivity with the King factor of the ab initio anisotropy of Raj, Hamaguchi and Witek (2018): the cross sections agree to 0.4 % over 400–1064 nm, the mean polarizabilities to 0.2 %, and the Raman-module H₂ anisotropy gives a King factor 0.35–0.73 % lower.

At 288.15 K over 350–800 nm the constant tensor and the RI method differ, as PolTen/RI − 1, by the amounts below. Refractivities are compared at the reference density of each gas: 2.546899×10¹⁹ cm⁻³ for N₂ and 2.68678×10¹⁹ cm⁻³ (273.15 K) for O₂. The O₂ RI refractive_index refers to that denser state, so comparing it with the PolTen default density would show a spurious ≈5 % difference in  . The cross sections, King factors and depolarizations are per molecule and do not depend on density.

Quantity	N₂	O₂
	−0.51 % (350 nm) to +0.30 % (800 nm)	−1.54 % (350 nm) to +0.63 %
	within ±0.12 %	within ±0.18 %
	−0.26 % to +0.18 %	−1.60 % to +0.98 %
	−7.3 % to +5.2 %	−15 % to +11 %

The depolarization differences arise mainly because the constant tensor anisotropy has no dispersion: rises with wavelength, while the RI King factors fall.    turns small King-factor differences into percent-level differences.

Model defaults

The automatic model_from_parameters elastic Rayleigh optical depth remains the terrestrial Bodhaine et al. (1999) air formula, scaled by surface pressure and rescaled by King factor to the configured depolarization. With the default depol: -1, the Rayleigh depolarization (which sets that King factor and the Rayleigh Greek coefficients) comes from the N₂/O₂ polarizability tensors at the band-centre wavenumber and the column-weighted mean temperature, with the tensor dispersion radiative_transfer.rayleigh_anisotropy (default DispersiveAnisotropy(); ConstantAnisotropy() reproduces the v2.2 pre-release depolarization bit for bit). Its variation across bands now follows the King factors and falls with wavelength; with the constant anisotropy only dispersed and ρ rose with wavelength. For air (VMR 0.8/0.2, 257 K) ρ is 0.0277 at 400 nm, 0.0264 at 550 nm, 0.0259 at 760 nm and 0.0254 at 1.6 µm (v2.2 pre-releases: 0.0263, 0.0271, 0.0276, 0.0280; v2.1.0, at 300 K and without 2π in the Buldakov pole: about 0.0281 from 400 nm to 2 µm). The Cabannes fraction always, and the Cabannes depolarization with depol: -1, use the Raman-module constants (getRamanAtmoConstants), as the Raman cross sections do. Under the dispersive default, a pure-elastic (noRS) run and a Raman run therefore no longer share the same elastic-plus-RRS decomposition (the implied ρ differs by up to ~10 % in the NIR); workflows that difference the two (Ring or filling-in spectra) should set rayleigh_anisotropy: ConstantAnisotropy(). Using the pure-gas helpers inside an RT run requires setting the optical-depth profile and the Greek source explicitly. Atmospheric composition is not yet connected to molecular columns and scattering-weighted phase-matrix mixing.

Elastic Rayleigh vs Cabannes in RT runs

For the molecular contribution (no aerosols, just air molecules):

• Pure-elastic runs (noRS) use the full Rayleigh Greek expansion. In this case rotational Raman is rolled into the effective depolarization ratio ; the phase matrix is the standard Rayleigh form with depolarization of ~0.028 for air at 532 nm.

• Raman-aware runs (RRS, VS_*) use the Cabannes Greek expansion. Rotational Raman is then handled explicitly by the inelastic kernel (Concepts/08), and the elastic Cabannes path uses a lower depolarization (~0.007).

The selector at src/CoreRT/LayerOpticalProperties/compEffectiveLayerProperties.jl:8–9:

_rayleigh_greek_source(::Union{noRS, noRS_plus}, greek_rayleigh, greek_cabannes) = greek_rayleigh
_rayleigh_greek_source(::AbstractRamanType, greek_rayleigh, greek_cabannes) = greek_cabannes

Mismatch — using Cabannes greek with noRS, or full Rayleigh greek with RRS — produces a ~1% bias on Stokes and a ~3% bias on , because the polarization-sensitive coefficients (, ) shift by the depolarization difference.

Aerosol microphysics → optics

The size distribution and refractive index are inputs to the Mie calculation. vSmartMOM ships with two microphysics schemes (in the Aerosols module):

• TOMAS15Scheme — TOMAS-15-bin sectional scheme used in CESM/GEOS-Chem.

• TwoMomentScheme — two-moment bulk scheme (AOD + effective radius, with fixed geometric width).

Plus a refractive-index database for common aerosol species (sulfate, dust, black carbon, organic carbon, sea salt). See Aerosols (microphysics) for the schemes (status: API still being stabilized — WIP).

Code anchors

Concept	Source
Mie entry (NAI-2)	src/Scattering/compute_NAI2.jl:82
Mie entry (PCW)	src/Scattering/compute_PCW.jl:37
Mie helper (a_n, b_n)	src/Scattering/mie_helper_functions.jl
GreekCoefs type	src/Scattering/types.jl:462–475
AerosolOptics type	src/Scattering/types.jl:512–533
Polarization types	src/Scattering/types.jl:142–215
Z from Greek + μ	src/Scattering/compute_Z_matrices.jl::compute_Z_moments
Truncation entry	src/Scattering/truncate_phase.jl (used by createAero)
Rayleigh / Cabannes selector	src/CoreRT/LayerOpticalProperties/compEffectiveLayerProperties.jl:8–9
GPU Mie kernel	src/Scattering/gpu_mie_kernels.jl:30–158
YAML decomposition map	src/IO/Parameters.jl::DECOMP_MAP

Hands-on tutorials

Runnable examples with Plotly figures:

• Scattering & Greek coefficients

• Mie deep dive

References

• Sanghavi (2014), Revisiting the Fourier expansion of Mie scattering matrices in generalized spherical functions, JQSRT 136:16–27, doi:10.1016/j.jqsrt.2013.12.015. NAI-2 vs PCW comparison; corrected Domke formalism.

• Sanghavi et al. (2014), JQSRT 133:412–433, doi:10.1016/j.jqsrt.2013.09.004, App. A.

• Hansen, J. E. & Travis, L. D. (1974), Light scattering in planetary atmospheres, Space Sci. Rev. 16:527. (Background on Mie + Greek matrix.)

• Sanghavi, S. & Stephens, G. (2015), Adaptation of the delta-m and δ-fit truncation methods to vector radiative transfer: Effect of truncation on radiative transfer accuracy, JQSRT 159:53–68, doi:10.1016/j.jqsrt.2015.03.007, Eqs. 21–22.

• Siewert, C. E. (1981, 1982), On the equation of transfer relevant to the scattering of polarized light, On the phase matrix basic to the scattering of polarized light. (NAI formalism.)

• Domke, H. (1974), The expansion of scattering matrices for an isotropic medium in generalized spherical functions, Astrophys. Space Sci. 29:379. (PCW basis.)

• de Rooij, W. A. & van der Stap, C. C. A. H. (1984), Expansion of Mie scattering matrices in generalized spherical functions, A&A 131:237. (Domke corrections that PCW further fixes.)

• Crib sheet: docs/dev_notes/theory_references.md §G.