Radiation API

Radiation

Radiative transfer in curved spacetime: fluid orthonormal tetrads, synchrotron transfer coefficients (thermal, power-law and kappa; unpolarized and the 8 polarized), the analytic constant-coefficient Stokes coupling, the coherency-tensor polarized march, and the AART equatorial g-disk source. Follows blacklight's radiation/ module; deviations are noted in the individual files.

SimCoeffConfig{T<:Real}

Scalar configuration for the simulation synchrotron model (blacklight defaults). Mutable so precompute_power_law! / precompute_kappa! can fill the derived Gamma/2F1/Bessel constants in place; plasma-value cuts use the "-1 disables" convention. Parametric over the scalar type T (default Float64): a ForwardDiff Dual element type lets AD sensitivities w.r.t. any config scalar (e.g. d_unit, e_unit, b_unit) flow through the transfer pipeline.

electron_theta(c, rho_cgs, pgas_cgs, beta_inv) -> Float64

Electron temperature Thetae (dimensionless, k Te / me c^2) from the titebeta model. Returns NaN for the unsupported `plasmathermalfrac == 0` corner (blacklight parity — the NaN makes the thetae plasma cuts never trip).

plasma_cut_omits(c, rho_cgs, n_e_cgs, pgas_cgs, theta_e, bb_cgs, sigma, beta_inv) -> Bool

Per-sample plasma-value cut test (blacklight PlasmaCutOmits). Returns true when the sample should be OMITTED (treated as vacuum). Each bound is guarded by >= 0.0 so a negative threshold (canonically -1.0) disables it; 0.0 is active. All comparisons are false for NaN, so a NaN theta_e never trips the theta_e cut.

plasma_decode(c, rho, pgas, b_sq) -> NamedTuple

Shared code-unit -> cgs plasma decode for both the ray-traced imaging (eight_coefficients) and the GRMHD Monte Carlo (plasma_at): cgs conversion -> electron_theta (titebeta) -> plasma_cut_omits gate. Returns (rho_cgs, pgas_cgs, n_e_cgs, bb_cgs, sigma, beta_inv, theta_e, omit::Bool); omit == true means the sample is cut (treat as vacuum). Allocation-free and typed on the bare SimCoeffConfig UnionAll so AD Dual configs match. The GPU kernels re-implement this arithmetic against a flat isbits params struct (a mutable config cannot live on-device) and must stay bit-exact with it.

fluid_tetrad(ucon, ucov, kcon, kcov, up_con, gcov, gcon; renorm_e1=false) -> SMatrix{4,4}

Orthonormal tetrad aligned with the fluid 4-velocity. Returns the 4x4 matrix whose row a holds the contravariant components e_(a)^mu. up_con is the spacelike up-direction (e.g. the magnetic field b^mu). With renorm_e1 = true the spatial basis vector e_(1) is renormalized against gcov (needed off single CKS, where |det g| != 1).

thermal_synchrotron(thermal_frac, n_e, theta_e, B, nu, sin_th) -> (j_inv, alpha_inv)

Unpolarized thermal synchrotron coefficients in Lorentz-invariant form (M 28/30/31), with alpha from Kirchhoff's law. Returns (0.0, 0.0) on any degenerate input. Follows blacklight ThermalSynchrotronInvariant.

thermal_synchrotron_polarized(thermal_frac, n_e, theta_e, B, nu, sin_th, cos_th) -> NTuple{8}

Thermal synchrotron polarized coefficients (Marszewski et al. 2021), returned in the order (j_i, j_q, j_v, alpha_i, alpha_q, alpha_v, rho_q, rho_v). cos_th is SIGNED (sets the sign of the circular quantities). Includes the theta_e ~ 0 cold-plasma trap and the Kirchhoff / underflow guards. Follows blacklight ThermalSynchrotronPolarized.

power_law_synchrotron_polarized(c, n_e, B, nu, sin_th, cos_th) -> NTuple{8}

Power-law synchrotron polarized coefficients in the order (j_i, j_q, j_v, alpha_i, alpha_q, alpha_v, rho_q, rho_v) (accumulation terms carrying plasma_power_frac). Follows blacklight PowerLawSynchrotronPolarized.

kappa_synchrotron_polarized(c, n_e, B, nu, sin_th, cos_th) -> NTuple{8}

Kappa-distribution synchrotron polarized coefficients in the order (j_i, j_q, j_v, alpha_i, alpha_q, alpha_v, rho_q, rho_v) (accumulation terms carrying plasma_kappa_frac). Follows blacklight KappaSynchrotronPolarized.

precompute_power_law!(c, gamma_min, gamma_max) -> c

Fill the p-dependent power-law constants (unpolarized power_jj/power_aa plus the polarized set) in place. Follows blacklight PrecomputePowerLawCoefficients. Call once when plasma_power_frac != 0.

precompute_kappa!(c) -> c

Fill the kappa/w-dependent constants (unpolarized + polarized subsets, including the three-way branch on kappa and the Bessel rhoV factor) in place, using gamma(), besselk(), and the host 2F1. Follows blacklight PrecomputeKappaCoefficients. Call once when `plasmakappa_frac != 0`.

Modified Bessel function K_0(x), x > 0 (A&S 9.8.5/9.8.6).

Modified Bessel function K_1(x), x > 0 (A&S 9.8.7/9.8.8).

Modified Bessel function K2(x), x > 0, via K2 = K0 + (2/x) K1.

harmonic_bridge(a, b, x) = (a^-x + b^-x)^(-1/x)

Harmonic-mean bridge between low- and high-frequency synchrotron asymptotes with exponent x (blacklight HarmonicBridge). IEEE-safe: an underflowed branch gives pow(tiny, -x) -> inf so the bridge returns the other asymptote; both underflowing gives 0, the correct limit.

stokes_couple(S, j, α, ρ, ds; dtau_max=100.0) -> SVector{4}

Advance the orthonormal-frame Stokes 4-vector S (I,Q,U,V) over one segment of invariant path length ds through plasma with invariant coefficients j, α, ρ (each an SVector{4}; the U component is zero). The analytic constant-coefficient solution (StokesCoupleNoSplit): emission-only, total absorption, pure rotation, polarized absorption without rotation, or the full absorption+rotation branch with the degeneracy guard. The result is clamped to be physically admissible (I >= 0, |P| <= I).

couple_unpolarized(I, j, α, ds; dtau_max=100.0) -> Float64

One exact-segment update of the invariant intensity through plasma with emissivity j and absorptivity α over invariant path length ds (UnpolarizedTransferStep): I e^{-dτ} + (j/α)(1 - e^{-dτ}) with dτ = α ds, capped at dtau_max (saturate to j/α), and the α <= 0 emission-only limit I + j ds.

eight_coefficients(cfg, rho, pgas, fv, gcov, gcon, kcov, kcon, ν_obs, renorm_e1=false)
    -> (j::SVector{4}, α::SVector{4}, ρ::SVector{4}, tet, valid::Bool)

Invariant polarized coefficients at one sample: thermal + power-law + kappa additive mixing (blacklight CalculatePolarizedCoefficients). rho/pgas are code-unit scalars; fv is the NamedTuple from primitive_four_vectors. Pass renorm_e1 = true off single CKS so the returned fluid tetrad's e_(1) (hence the Q/U frame) is unit.

unpolarized_march(st, cfg, sampler, ys, dls, ν_obs, L_unit) -> Float64

Stokes-I fast mode: stable exponential propagator over the invariant intensity (blacklight unpolarizedtransfer). Returns Iν at the camera.

polarized_march(st, cfg, sampler, ys, dls, ν_obs, L_unit;
                cam_up, intensity_only=false) -> SVector{4}

March ONE ray's coherency tensor from source to camera. ys/dls are the tracerays outputs for this ray (source→camera order). sampler(x) returns the NamedTuple (rho, pgas, uu, bb, inside) at position x (code units). Returns the camera-frame Stokes (I, Q, U, V) at νobs (physical, scaled by ν_obs³).

EquatorialDisk(a; sub_kep=1, betar=1, betaphi=1,
               gamma_p=-1.5, mu_p=1-√(1-a²), sigma_p=0.5,
               redshift_power=3, rin=0.0)

isbits bundle of the AART g-disk knobs, threaded through render_analytic (and the equatorial_intensity seam) so the analytic tracer renders the g-disk photon-ring image. equatorial_disk_intensity(d, rs, νr, a, λ, η) evaluates the observed intensity for this disk. AART's canonical M87* case is EquatorialDisk(0.94).

equatorial_disk_intensity(rs, νr, a, λ, η;
    sub_kep=1, betar=1, betaphi=1,
    gamma_p=-1.5, mu_p=1-√(1-a²), sigma_p=0.5,
    redshift_power=3, rin=0.0) -> I_obs

Observed specific intensity of the AART optically-thin equatorial g-disk at source radius rs for a ray with emitted-p^r sign νr and conserved quantities (λ, η) (AART bright_radial, P1 Eq. 50):

I_obs = johnson_su_profile(rs) · gfactor_disk(rs)^redshift_power.

redshift_power is AART's gfactor exponent (σ = 3 for a frequency-integrated optically-thin disk: g³ combines the g² solid-angle boost with the g¹ intensity Lorentz-invariance I_ν/ν³). Zero inside the horizon or rs < rin.

gfactor_disk(rs, a, b, λ, η; sub_kep=1, betar=1, betaphi=1) -> g

The mixed-flow equatorial redshift factor: circular (gfactor_circular) for rs ≥ r_ISCO, plunging (gfactor_plunge) for rs < r_ISCO — AART's bright_radial gDisk/gGas selector (P1 Eq. B13). Continuous across the ISCO by construction; b is the emitted-p^r sign (νr). In the a→0, all-knobs=1 limit this is the Schwarzschild circular-orbit redshift g = 1/(uᵗ(1 − λΩ_K)); at large r, g → 1.

gfactor_circular(rs, a, b, λ, η; sub_kep=1, betar=1, betaphi=1) -> g

Redshift factor for a photon emitted OUTSIDE the ISCO by the (generally sub-Keplerian, mixed radial/azimuthal) circular-ish flow — AART gDisk (P1 Eq. B13). b is the sign of the emitted p^r (the redshift_sign / νr from analytic_emission). The velocity knobs interpolate between profiles: sub_kep scales the specific angular momentum (1 = Keplerian ℓ), betaphi blends Ω̂ (1) toward the frame-dragging Ω̄ (0), betar blends uʳ from 0 (1) toward pure radial infall (0). Defaults (all 1) = the equatorial Keplerian disk.

gfactor_plunge(rs, a, b, λ, η; sub_kep=1, betar=1, betaphi=1) -> g

Redshift factor for a photon emitted INSIDE the ISCO by the plunging flow — AART gGas (P1 Eq. B13). The plunging material conserves the ISCO angular momentum ℓms and energy Êms; its radial velocity ûʳ(r) is set by the geodesic radial potential (ν̂), continuously → 0 at the ISCO so gGas ≡ gDisk there.

johnson_su_profile(r; mu=0.0, sigma=0.5, gamma=-1.5) -> emissivity

The AART Johnson-SU radial emissivity envelope (P1 Eq. 58; cf. 2008.03879 Eq. 12): f(r) = exp(−½(γ + asinh((r−μ)/σ))²) / √((r−μ)² + σ²). mu sets the location, sigma the width, gamma the skewness; strictly positive and integrable. AART's canonical disk uses γ = −3/2, μ = 1 − √(1−a²), σ = 1/2.