Monte Carlo API

MonteCarlo

BlackLightPlus Monte Carlo radiative-transfer module, following grmonty (Dolence et al. 2009) and its GPU descendants (gpumonty, igrmonty). Provides a counter-based Philox RNG, thermal-synchrotron emission physics, superphoton generation, tracking, Compton scattering, and spectrum binning. All kernel code is pure, allocation-free, isbits and StaticArrays-typed, with no global RNG state — one source for CPU, GPU (KernelAbstractions) and AD.

philox4x32(ctr, key) -> NTuple{4,UInt32}

Philox4x32-10 counter-based bijection: maps the 128-bit counter ctr (NTuple{4,UInt32}) and 64-bit key (NTuple{2,UInt32}) to four uniformly distributed 32-bit words. Pure, allocation-free, isbits-in/isbits-out.

philox_rand(key, ctr) -> Float64 in [0, 1)

One uniform double from the Philox stream key (NTuple{2,UInt32}) at counter ctr (NTuple{4,UInt32}). Two output words are combined into a ~53-bit mantissa so the result is strictly in [0, 1) (max value 1 - 2^-64). Pure/device-safe.

Convenience seeding: philox_rand(global_seed, stream_id, draw).

MCRng(key::NTuple{2,UInt32}, ctr::UInt32)

Immutable, isbits handle for a single superphoton's random stream. key identifies the stream (typically (global_seed, superphoton_id)); ctr is the next draw index. All samplers are pure functions of an MCRng and return the advanced handle, so no mutable/global state is ever touched — the same call runs on the CPU and inside a KA kernel.

mc_uniform(rng) -> (Float64 in [0,1), MCRng)

Draw one uniform and return it together with the counter-advanced stream.

mc_gauss(rng) -> (Float64, MCRng)

One standard-normal deviate via the Box–Muller transform (rejection-free, so no warp divergence). 1 - u1 keeps the logarithm finite for u1 -> 0.

sample_direction(rng) -> (SVector{3,Float64}, MCRng)

Isotropic unit vector on S^2. Uses the direct inversion z = 2u1 - 1, phi = 2 pi u2 (rejection-free, unlike gpumonty's polar method) so it never diverges across a warp. |v| = 1 exactly up to round-off.

standard_gamma(shape, rng) -> (Float64, MCRng)

A Gamma(shape, 1) deviate. shape >= 1 uses Marsaglia & Tsang (2000) — one normal + one uniform per trial, ~96%+ acceptance. 0 < shape < 1 uses the exact boosting identity Gamma(a) = Gamma(a+1) * U^{1/a} (rejection-free), which is friendlier to SIMT than gpumonty's nested Ahrens–Dieter branch. Clean-room, device-safe counterpart of curand.cu legacy_standard_gamma (algorithm per Marsaglia & Tsang 2000; no GPL source copied).

chi_square(df, rng) -> (Float64, MCRng)

A chi-square deviate with df degrees of freedom via 2 * Gamma(df/2, 1) (the same mathematical identity curand.cu's chi_square uses); used by the relativistic-Maxwellian electron sampler.

Superphoton

Monte Carlo superphoton state (isbits).

Fields

• x::SVector{4,Float64} : coordinate position x^mu.
• p::SVector{4,Float64} : COVARIANT wavevector pmu (me c^2 units).
• w::Float64 : statistical weight (photons per packet, >> 1).
• E::Float64 : conserved energy -p_t.
• tau_abs::Float64 : accumulated absorption optical depth.
• tau_scatt::Float64 : accumulated scattering optical depth.
• nscatt::Int32 : number of scatterings undergone.
• x1i,x2i::Float64 : emission-point coordinates (diagnostics).
• ne0,thetae0,b0::Float64 : plasma state at emission (ne, Thetae, B) (diagnostics).

Coordinate position x^mu.

Covariant wavevector p_mu.

Statistical weight (photons per packet).

Conserved photon energy E = -p_t.

Number of scatterings undergone.

j_nu(n_e, theta_e, B, nu, sin_theta) -> Float64

Non-invariant thermal-synchrotron emissivity jnu [erg s^-1 cm^-3 Hz^-1 sr^-1] in cgs, from `Radiation.thermalsynchrotron(jinv = jnu/nu^2) timesnu^2.sin_theta` is the pitch angle between the photon and the magnetic field.

alpha_nu(n_e, theta_e, B, nu, sin_theta) -> Float64

Non-invariant thermal-synchrotron absorptivity alphanu [cm^-1] from Kirchhoff's law, via `Radiation.thermalsynchrotron(alphainv = alphanu*nu), soalphanu = alphainv / nu`.

jmax(n_e, theta_e, B, nu) -> Float64

Rejection-sampling envelope for the pitch-angle distribution: j_nu at sin_theta = 1, which maximises the thermal-synchrotron emissivity over pitch angle for nu >~ nus. The direction is accepted with probability `jnu(theta)/j_nu(pi/2)` (as in grmonty/gpumonty).

int_jnu_zone(n_e, theta_e, B, nu) -> Float64

Pitch-angle-integrated emissivity ∫ j_nu dΩ [erg s^-1 cm^-3 Hz^-1] — the per-zone emission rate at frequency nu. Since j_nu depends on direction only through sin(theta) and is symmetric about the equator,

∫ j_nu dΩ = 4 pi ∫_0^{pi/2} j_nu(sin theta) sin theta d theta,

evaluated with a fixed composite Simpson rule (device-safe: constant loop bound, no allocation); validated against an independent fine quadrature to < 1e-6.

emission_tetrad(ucon, ucov, bcon, gcov, gcon) -> SMatrix{4,4}

(u, b-hat)-aligned orthonormal emission tetrad; row a holds e_(a)^mu (the layout of Radiation.fluid_tetrad): e(0) = uhat, e(1) = bhat (sets the pitch angle), e(2)/e(3) completed by Gram–Schmidt on x2, x3 with the metric inner product (after gpumonty tetrads.cu make_tetrad). ucov is accepted for signature parity with the imaging tetrad. A numerically null bcon falls back to the x1 direction.

EmissionZone

One emitting fluid cell, the host-side input to superphoton generation (isbits, StaticArrays-typed). Plasma state is cgs; the geometry is BlackLightPlus-native (signature (-,+,+,+), covariant/contravariant SVector{4}s and SMatrix{4,4} metrics).

Fields

• n_e : electron number density [cm^-3].
• theta_e : dimensionless electron temperature Thetae = k Te / m_e c^2.
• B : magnetic-field strength [Gauss].
• dV : proper emitting volume sqrt(-g) dx^3 L_unit^3 [cm^3] — the single factor multiplying int_jnu_zone to get an emission rate.
• x : coordinate position x^mu (emission point).
• ucon,ucov : fluid 4-velocity, contravariant / covariant.
• bcon : magnetic-field 4-vector b^mu (its spatial direction sets b-hat).
• gcov,gcon : metric and inverse metric at x.

WeightTable

Precomputed superphoton-weight function w(nu) on a log-frequency grid. lwgt[i] holds log w at node nu_i = exp(lnu_min + (i-1) dlnu), i = 1 … N_ESAMP+1. weight_at reads it back by log-nu interpolation. Host storage (Vector).

weight_at(nu, wgt::WeightTable) -> Float64

Superphoton weight w(nu) [photons/packet] by piecewise log-linear interpolation in ln nu (gpumonty linear_interp_weight), clamped to the grid endpoints. Host wrapper over the shared _mc_weight_at (mc_common.jl).

build_weight_table(zones, Ns; n_esamp, nu_min, nu_max) -> WeightTable

Build w(nu) = Sum_zones int_jnu_zone(nu) dV / (h Ns) on a log-nu grid of n_esamp+1 nodes over [nu_min, nu_max], storing log w. Ns is the target total superphoton count (sets the packet weight; the emergent physics is independent of it). After gpumonty init_weight_table.

build_weight_table(jt::ZoneEmissionTable, Ns; nthreads) -> WeightTable

Fast-path weight table from a prebuilt ZoneEmissionTable: w(nu_m) = (Sum_z jt.jnu[m,z] jt.dV[z]) / (h Ns), stored as log w. Threaded over zones into per-thread partial accumulators reduced in fixed thread-index order, so the result is deterministic but may differ from the serial method in the last bits (different zone summation order).

zone_nsuperphotons(zone, wgt) -> (N_s, dnmax)

Expected superphoton count and frequency-rejection envelope for zone:

N_s   = (dV / h) Int dln(nu) int_jnu_zone(nu) / w(nu)
dnmax = max_nu  int_jnu_zone(nu) / w(nu)

Both use the tabulated grid-node weights (parity with gpumonty init_zone).

sample_zone_photon(zone, wgt, dnmax, rng[, gt]) -> (Superphoton, MCRng)

Generate one superphoton from zone using the prebuilt weight table wgt and the zone's frequency-rejection envelope dnmax (from zone_nsuperphotons). Returns the packet and the counter-advanced RNG. Allocation-free; pure over (zone, wgt, dnmax, rng). Passing an IntJnuTable as gt switches the frequency-rejection numerator from the live pitch integral to the tabulated interpolant (opt-in speed path; dnmax must then be the matching tabulated_dnmax envelope so the acceptance stays <= 1).

sample_zone_photon(zone, wgt, rng) -> (Superphoton, MCRng)

Convenience method that derives the rejection envelope dnmax from the zone itself (via zone_nsuperphotons) before sampling. Prefer the 4-argument form in a hot loop and pass a dnmax computed once per zone.

ZoneEmissionTable

Cached per-(zone, frequency-node) pitch-angle-integrated emissivity for a GRMHD emitting-zone list: jnu[m+1, z] = int_jnu_zone(zone_z, nu_m) on the shared log-nu grid, dV[z] the zone's proper emitting volume. Built once and consumed by both build_weight_table and zone_counts. Host storage.

build_zone_emission_table(zones; n_esamp, nu_min, nu_max) -> ZoneEmissionTable

Compute the (n_esamp+1) x nzones emission-integral table jnu[m+1,z] = int_jnu_zone(zone_z, nu_m) once, threaded over zones. Each thread writes only its own zone columns, so the fill is race-free and bit-identical regardless of thread count.

zone_counts(jt::ZoneEmissionTable, wgt::WeightTable; nthreads)
    -> (Ns::Vector{Float64}, dnmax::Vector{Float64})

Per-zone expected superphoton count N_s[z] and frequency-rejection envelope dnmax[z] for every zone in jt, read from the cached jt.jnu. Threaded over zones (independent per-zone work, bit-identical regardless of thread count). jt and wgt must share the frequency grid (asserted).

kn_total(eps_e) -> Float64

Klein-Nishina TOTAL cross section normalised by the Thomson value, sigma_KN(eps_e) / sigma_T, for a photon of dimensionless energy eps_e in the electron rest frame. Uses the numerically stable low-energy expansion 1 - 2 eps_e for eps_e < 1e-3 (avoids catastrophic cancellation in the grouped form), else the closed KN form. -> 1 (Thomson) as eps_e -> 0 and decreases monotonically. Pure / allocation-free / device-safe.

boostcross(eps, mu_e, gamma_e) -> Float64

Doppler-boosted Klein-Nishina cross section (over sigmaT) for a lab-frame photon of energy eps colliding with an electron of Lorentz factor `gammaeat collision cosinemue. The photon energy in the electron rest frame isepse = eps gammae (1 - mue betae)(betae = sqrt(gammae^2-1)/gammae) and the result is weighted by the relative flux factor(1 - mue betae)`.

dNdgamma_e(theta_e, gamma_e) -> Float64

Normalised Maxwell-Juttner distribution dN/dgamma_e per unit electron number density (integrates to 1 over gamma_e in [1, inf)). The normalisation uses K2f = K2(1/Theta_e) exp(1/Theta_e) for Theta_e > 1e-2; below that the product overflows/underflows to NaN, so the large-argument asymptotic series K2f ~ sqrt(pi Theta_e/2)(1 + 15 Theta_e/8 + 105 Theta_e^2/128) (A&S 9.7.2) is used — carried through Theta_e^2 so the regimes agree to ~2e-7 at the switch.

total_compton_cross_num(eps, theta_e) -> Float64

Thermally-averaged hot Compton cross section sigma_h(eps, Theta_e) [cm^2], by 2D midpoint quadrature over mu_e in [-1,1] and gamma_e in [1, 1 + MAXGAMMA Theta_e], integrating (1 - mu_e beta_e) sigma_KN(eps_e) against dN/dgamma_e. Cold corners short-circuit: Theta_e < MINT && eps < MINW -> sigma_T; Theta_e < MINT -> kn_total(eps) sigma_T. Warp-uniform trip counts, so device-safe (after gpumonty total_compton_cross_num).

HotcrossTable

Precomputed hot-cross-section lookup table: the (NW+1) x (NT+1) matrix of log10(sigma_h [cm^2]) on the log-log (eps, Theta_e) grid, with logsigma[i+1, j+1] at eps = 10^(log10(MINW) + i dlw), Theta_e = 10^(log10(MINT) + j dlT). A plain Matrix: a shared read-only buffer the GPU uploads directly.

build_hotcross_table() -> HotcrossTable

Regenerate the full (NW+1) x (NT+1) table by evaluating total_compton_cross_num at every grid node (log10 of the cm^2 result). Host-side offline builder; production reads the cached file via load_hotcross. The energy loop is threaded over independent columns, so the fill is race-free and bit-identical regardless of thread count.

load_hotcross(path) -> HotcrossTable

Parse a gpumonty hotcross.dat file (i j log10(eps) log10(Theta_e) log10(sigma_h) per line, (NW+1)*(NT+1) rows) into a HotcrossTable. Only the integer grid indices and the log10(sigma_h) column are used, so the table is laid out independently of the file's row order.

hotcross_lookup(table, eps, theta_e) -> Float64

Hot Compton cross section sigma_h(eps, Theta_e) [cm^2] via the log-log table. Branch order (matches gpumonty total_compton_cross_lkup, ORDER MATTERS):

1. eps*Theta_e < 1e-6 -> Thomson sigma_T (cold, soft),
2. Theta_e < MINT -> cold Klein-Nishina kn_total(eps) sigma_T,
3. in-bounds -> bilinear interpolation in log10(eps), log10(Theta_e),
4. otherwise -> direct total_compton_cross_num fallback.

Device-safe (only reads the passed array + arithmetic).

kappa_es(table, nu, theta_e) -> Float64

Electron-scattering mass opacity sigma_h / m_p [cm^2 g^-1] at frequency nu (cgs) and electron temperature Theta_e. The dimensionless photon energy is eps = HPL nu / (m_e c^2). Matches grmonty kappa_es so that alpha_scatt = nu kappa_es n_e m_p.

alpha_inv_scatt(table, nu, theta_e, n_e) -> Float64

Lorentz-INVARIANT scattering absorption coefficient alpha_inv,s = nu n_e sigma_h at frequency nu (cgs), electron temperature Theta_e and number density n_e. Slots into the same invariant-alpha march as the synchrotron absorptivity.

lorentz_boost(v, u) -> SVector{4,Float64}

General special-relativistic boost of the 4-vector v into the frame whose 4-velocity (in the current frame) is u = (gamma, gamma*beta*n_hat). The inverse boost is lorentz_boost(v, (u[1], -u[2], -u[3], -u[4])) (reversed spatial velocity). Pure / allocation-free; _BOOST_SMALL guards the beta -> 0 direction normalisation (the boost then reduces to the identity). Signature (-,+,+,+).

klein_nishina(a, ap) -> Float64

The (unnormalised) Klein-Nishina differential cross section factor for incoming / outgoing rest-frame photon energies a / ap, (a/ap + ap/a - 1 + ch^2) / a^2 with ch = 1 + 1/a - 1/ap. Used as the rejection envelope target in sample_klein_nishina_energy.

sample_scattered_energy(k0, rng) -> (k0p, cos_theta, rng)

Scattered photon energy k0p and scattering cosine cos_theta in the ELECTRON REST frame for an incident rest-frame energy k0. Elastic Thomson for k0 <= 1e-4 (k0p = k0, cos_theta from the dipole), else Klein-Nishina (k0p from the KN spectrum, cos_theta = 1 - 1/k0p + 1/k0 — the Compton relation 1/k0p - 1/k0 = 1 - cos_theta).

sample_electron(theta_e, k_tet, rng) -> (p_tet::SVector{4}, rng)

Sample the 4-velocity p_tet = (gamma_e, gamma_e beta_e n_hat) of the electron a photon of tetrad-frame wavevector k_tet scatters off: (gamma_e, beta_e) from the Maxwell-Juttner distribution, the collision cosine mu from the (1 - beta_e mu) flux distribution, accepted with the Klein-Nishina probability kn_total(gamma_e (1 - beta_e mu) k_tet[1]). The direction is built in a random triad about the photon axis. Rejection-capped for device safety.

sample_electron_gamma(theta_e, rng) -> (gamma_e, beta_e, rng)

Draw an electron Lorentz factor / speed from the relativistic Maxwell-Juttner distribution at temperature Theta_e (Canfield-Howard-Liang y sampler: gamma_e = y^2 Theta_e + 1, beta_e = sqrt(1 - 1/gamma_e^2)). This is the BARE thermal distribution (no flux / Klein-Nishina weighting); its mean is <gamma_e> = 3 Theta_e + K1(1/Theta_e)/K2(1/Theta_e).

sample_mu(beta_e, r) -> Float64

Sample the cosine $\mu$ of the angle between the photon and electron momenta from the $(1 - \beta_e \mu)$ collision-flux distribution.

beta_e is the electron speed in units of $c$ and r is a uniform variate on $[0,1)$; the result is drawn by inverse CDF and clamped to $[-1, 1]$. The discriminant satisfies $(1-\beta_e)^2 \le \mathrm{det} \le (1+\beta_e)^2$, so the abs guarding the square root is defensive only.

See also sample_electron, compton_scatter.

compton_scatter(k_tet, theta_e, rng) -> (k_scattered_tet::SVector{4}, rng)

Scatter a photon of tetrad-frame wavevector k_tet off a thermal electron gas at temperature Theta_e, returning the scattered wavevector in the SAME tetrad frame. Samples the electron 4-velocity, boosts the photon into its rest frame, draws the scattered energy + angle (Thomson or Klein-Nishina), rebuilds the outgoing direction in a random triad about the incident rest-frame direction, and boosts back. The output is null (k'.k' = 0) to round-off.

scatter_superphoton(p_cov, ucon, ucov, bcon, theta_e, gcov, gcon, rng)
    -> (p_new_cov::SVector{4}, E_new, rng)

Coordinate-frame Compton scatter driver: builds the fluid (e_(0) = u) emission tetrad, transforms the incident covariant p_cov into it, scatters (compton_scatter), transforms back to the covariant p'_mu, and returns it with the new energy E' = -p'_t. Only e_(0) = u matters physically; the spatial tetrad axes cancel under the randomised scatter.

scatter_superphoton(sp::Superphoton, ucon, ucov, bcon, theta_e, gcov, gcon, rng)
    -> (Superphoton, rng)

Superphoton overload of the driver: scatters sp in place of its momentum, increments nscatt, and rebuilds the packet with the new covariant momentum (E = -p'_t recomputed by the constructor). Position and weight are unchanged (scattering is local).

track_superphoton(sph::Superphoton, model, rng) -> (Superphoton, recorded::Bool)

Transport sph through model until it escapes the sphere, decaying its weight by the accumulated absorption optical depth. Returns the escaped packet and recorded = true; a packet that hits the step cap (never, in practice) returns recorded = false. Absorption-only: rng is accepted but unused.

track_superphoton_scatter!(ws, sph, model, table, rng; max_scatt)
    -> (escaped::Superphoton, recorded::Bool, scatter_loss::Float64, rng::MCRng)

Transport one superphoton through the uniform flat sphere with the biased exponential-transform scatter decision (see the section header), spawning Compton-scattered children onto ws. The march is the shared _mc_march_scatter_core (mccommon.jl); this wrapper supplies the host opacity lookup and the ScatterWorkspace child sink. Returns the escaped incident photon (recorded=true if it left the sphere), the incident weight lost to scattering (`scatterloss`, checks unbiasedness) and the advanced RNG. Absorption is neglected (idealized Comptonization sphere).

bias_func(n_e, R, theta_e, w, nscatt) -> Float64

Importance bias for the exponential-transform scatter decision, sphere form (arXiv:2602.13198 §3.2; cf. gpumonty): b = max(1, n_e sigma_T R (n_sc+1)^2), growing with scattering generation n_sc so successively rarer higher-order scatterings stay well-sampled. theta_e/w are accepted for signature parity with gpumonty bias_func (unused in the sphere form).

ScatterWorkspace

Per-thread iterative work queue for the layered (non-recursive) multi-scatter loop: a LIFO stack of pending superphotons and their independent RNG streams, plus running diagnostics — total spawned child weight, first-order child weight (spawned by primaries) and child count. Children get a stream MCRng(child_seed, next_id) disjoint from the primaries'.

Push a photon (typically a freshly generated primary) and its stream onto the queue.

Pop the most recently pushed (photon, rng) (depth-first: bounded stack depth).

Spectrum

Binned escaped-superphoton accumulator. count/sumw/sumwE are n_thbins * n_ebins matrices holding, per cell, the raw superphoton count, Sum w (dNdlE) and Sum w E (dEdlE, E in me c^2 units). Frequency bin i (1-based) is centred at `exp(lnumin + (i - 0.5) dlnu); angle binjspans[(j-1) dth, j dth]of the equator-folded polar angle,dth = (pi/2)/n_thbins. Host storage (mutableMatrix`), reduced across threads after the run.

record!(spec, sph) -> spec

Bin the escaped packet sph by exit angle and energy, adding its count, w and w E to the cell. E = -p_t (me c^2 units) fixes `nu = E me c^2 / hand hence the frequency bin; the exit polar angle is taken about the x3 (z) axis and folded about the equator (mu = |k^3| / |kspatial|, valid in the flat escape region wherek^i = pi`). Packets outside the frequency grid are dropped.

record_energy!(spec, sph) -> spec

Bin the escaped packet sph by energy only into row 1 (nu = E m_e c^2 / h), ignoring direction. Used for the GRMHD per-inclination spectra, where the exit angle is computed externally from the asymptotic geodesic direction and each inclination band owns a single-row Spectrum. Off-grid packets are dropped.

merge_spectra!(dst, src) -> dst

Add src's accumulators into dst (grid-compatible). Used to reduce per-thread spectra after a threaded run.

Frequency bin centres [Hz].

nuLnu(spec) -> (nu, L)

Per-angle-bin nu L_nu [erg s^-1]: L[j, i] = (m_e c^2)(4 pi / dOmega_j)(1 / dln nu) Sum w E (gpumonty report_spectrum convention). nu are the bin centres.

nuLnu_total(spec) -> (nu, L, count)

Angle-summed (total) nu L_nu(nu) = (m_e c^2)(1 / dln nu) Sum_{all theta} w E [erg s^-1], with the per-bin raw superphoton count (summed over angle) for the population cut. This is Sum_theta nuLnu(theta) dOmega / (4 pi) — the convention matching analytic_nuLnu.

nuLnu_solid(spec, dOmega) -> (nu, L, count)

Isotropic-equivalent nu L_nu(nu) = (m_e c^2)(4 pi / dOmega)(1 / dln nu) Sum w E [erg s^-1] for a single-row (angle-selected) spectrum whose escaped photons subtend solid angle dOmega [sr] — the grmonty report_spectrum per-viewing-angle convention (the equivalent all-sky luminosity a distant observer at that inclination infers). count is the per-bin raw photon count.

total_power(spec) -> Float64

Bolometric recorded power Sum_all w h nu = (m_e c^2) Sum_all w E [erg s^-1] over every cell — the Monte Carlo estimate of the sphere's frequency-integrated luminosity (optically thin: all emitted energy escapes).

total_weight(spec) -> Float64

Total recorded superphoton weight Sum_all w over every cell (the packet-number content of the spectrum, used to normalise per-order populations).

mean_energy(spec) -> Float64

Weight-averaged escaped-photon energy <E> = Sum w E / Sum w [me c^2 units] over all cells (0 if the spectrum is empty). For the per-scattering-order spectra of `mcspectrum_scatter` this is the mean energy of that Compton bump, so the ratio between successive orders measures the up-scattering amplification.

SphereModel

Uniform optically-thin synchrotron sphere (isbits). Carries the plasma state, radius, the emitting volume dV = (4/3) pi R^3, and the (position-independent) flat-space kinematics: static fluid ucon/ucov, magnetic-field direction bcon (along x1), and the Minkowski metric gcov/gcon.

sphere_zone(model, x) -> EmissionZone

The emitting zone at coordinate x: the sphere's uniform plasma and flat-space kinematics with the whole-sphere volume dV. (For weight-table purposes one such zone represents the entire sphere; for per-photon generation only x varies.)

sample_sphere_photon(model, wgt, dnmax, rng[, gt]) -> (Superphoton, MCRng)

Generate one superphoton: draw a uniform-random position in the ball of radius R (r = R U^{1/3}, isotropic direction), build the local EmissionZone, and sample_zone_photon. Allocation-free; returns the counter-advanced RNG. gt (with its matching tabulated_dnmax as dnmax) selects the opt-in tabulated frequency rejection — see sample_zone_photon.

sample_sphere_photon(model, wgt, rng) -> (Superphoton, MCRng)

Convenience form that derives the frequency-rejection envelope dnmax from the sphere zone before sampling. Prefer the 4-argument form in a hot loop.

analytic_nuLnu(model, nu) -> Float64

Closed-form optically-thin thin-sphere spectrum [erg s^-1]:

nu L_nu = nu * V * int_jnu_zone(n_e, Theta_e, B, nu),   V = (4/3) pi R^3,

the frequency-times-volume-times-(angle-integrated)-emissivity, in the SAME angle-integrated (total) convention as nuLnu_total. This is the gpumonty python checker's Luminosity_analytic * nu * dv (its int_jnu equals int_jnu_zone).

mc_spectrum(model, Ns; nthreads, seed, nu_min, nu_max, n_esamp, n_ebins, n_thbins)
    -> Spectrum

Run the sphere Monte Carlo: build the one-zone weight table, generate ~Ns superphotons (uniform-random emission points in the ball), track each to escape and bin it into the returned Spectrum. Threaded with one accumulator per thread, reduced at the end; seed fixes the Philox streams. nu_min/nu_max set both the generation band and the spectrum grid. tabulated_gen = true opts into the tabulated frequency rejection (much cheaper generation, statistically identical but not guaranteed bitwise).

ScatteringSphere(; n_e=1e18, theta_e=3, B=100, R=3e5) -> SphereModel

The idealized Comptonization sphere (arXiv:2602.13198 §3.2 regime): a uniform thermal plasma that emits thermal-synchrotron seed photons and Compton-scatters them, sized so tau_es = n_e sigma_T R ~ 0.2 is moderate (successive up-scattering bumps visible). Returns a plain SphereModel; the scattering transport is selected by mc_spectrum_scatter. Synchrotron self-absorption is neglected in the scattering transport (it exceeds the Thomson opacity by ~1e9-1e12 at the emission peak, so a moderate-tau_es sphere has no closed-form seed); analytic_nuLnu is then the direct emission the unscattered nscatt=0 spectrum reduces to, up to exp(-tau_scatt).

scattering_tau_es(model) -> Float64

Radial Thomson scattering optical depth tau_es = n_e sigma_T R of the sphere (the soft-photon limit sigma_h -> sigma_T; the geometry-averaged escape depth over uniform isotropic emission is smaller, <tau> ~ 0.75 tau_es).

mc_spectrum_scatter(model, table, Ns; kwargs...) -> ScatterResult

Run the scattering-sphere Monte Carlo: build the synchrotron weight table over the generation band [gen_nu_min, gen_nu_max], generate ~Ns seed superphotons, and track+scatter each into per-order spectra on the wider spectrum grid [spec_nu_min, spec_nu_max] (which must reach the up-scattered bumps). table is the hot-cross lookup; max_scatt caps the scattering generations. Threaded over primaries with private per-order accumulators; children get RNG streams disjoint from the primaries'.

ScatterResult

Output of mc_spectrum_scatter. orders[n+1] is the escaped-Spectrum of photons with scattering order n (orders[end] lumps nscatt >= max_scatt); total sums all orders. Diagnostics: primary/child counts, the created primary and child weights, the first-order child weight (spawned by primaries) and scatter_loss, the total incident weight lost to scattering (the unbiasedness conservation law checks w_children ≈ scatter_loss).

GRMHDModel

The GRMHD Monte Carlo emission model: an AthenaK snapshot + a Cartesian Kerr-Schild spacetime + the cgs unit system and plasma configuration. NOT isbits (it wraps the snapshot arrays and the mutable SimCoeffConfig) — it is the HOST description consumed by emission_zones (to build the generation tables) and plasma_at (sampled along each superphoton geodesic during tracking).

Fields

• spacetime : KerrSchildCartesian(a).
• snap : AthenaKSnapshot (the fluid primitives + AMR geometry).
• cfg : SimCoeffConfig — units (d_unit/e_unit/b_unit), the titebeta electron model (rat_high/rat_low/use_p/gammas) and the plasma cuts (cut_sigma_max, cut_theta_e_min, ...).
• gamma : adiabatic index used for p_gas = (gamma-1) eint.
• L_unit : GM/c^2 [cm] — the code-length -> cm scale (fixes the optical-depth geometric factor and the proper cell volume).
• r_h : outer horizon KS radius (capture surface).
• R_em : emission-confinement radius [M]; a photon outbound past this streams freely to infinity (its escape direction is frozen — see the tracker).

plasma_at(model, x::SVector{4}) -> NamedTuple

Emitting-plasma state at CKS 4-position x, decoded EXACTLY as the imaging eight_coefficients does: trilinear sample_primitives -> primitive_four_vectors -> cgs conversion -> electron_theta (titebeta) -> plasma_cut_omits. Returns (n_e, theta_e, B, ucon, ucov, bcon, gcov, gcon, inside); inside = false (with zeroed plasma but valid gcov/gcon) for vacuum, off-grid, non-positive density, or any cut-omitted (e.g. sigma > sigmacut jet) sample. B is in Gauss, `nein cm^-3,theta_e` dimensionless.

emission_zones(model; max_zones=typemax(Int), stride=1, r_emit_max=Inf)
    -> Vector{EmissionZone}

Iterate the snapshot's cells and return one EmissionZone per emitting cell (plasma passes plasma_at's cuts, n_e > 0, theta_e > 0), with proper volume dV = sqrt(-g) dx1 dx2 dx3 L_unit^3. A cell is skipped when a finer AMR block owns its centre (de-duplication). stride > 1 subsamples the grid, inflating each retained cell's dV by stride^3 — a fast, unbiased estimator for tests. r_emit_max drops cells beyond that KS radius: the tenuous outer atmosphere passes the cuts but emits negligibly while forcing whole-domain tracking, so capping to the luminous torus is a large speedup at no luminosity cost (Inf disables). max_zones caps the list (diagnostic guard).

track_superphoton_grmhd(sph, model, C, R_record) -> (Superphoton, recorded::Bool)

Transport sph along the curved KS null geodesic (fixed-step RK4, affine step from _grmhd_dlam), decaying its weight by the deterministic absorption depth dtau_a = 1/2[(nu alpha_nu)_n + (nu alpha_nu)_{n+1}] (L_unit h / m_e c^2) dl (grmonty dtauK), sampling the fluid each step. Returns the escaped packet (recorded = true) on reaching R_record or streaming outbound past Rem; recorded = false on horizon capture or the step cap. Absorption-only.

track_superphoton_grmhd_scatter!(ws, sph, model, table, C, R_record, bias0,
                                 max_scatt, rng)
    -> (escaped::Superphoton, recorded::Bool, rng::MCRng)

Curved-geodesic tracker with the biased exponential-transform Compton scatter decision (grmonty track_super_photon): every step accumulates absorption and scattering optical depth, depletes the incident beam by the physical exp(-(dtau_a + dtau_s)), and on a trigger spawns a Compton-scattered child at a sub-step-interpolated event point (new energy/direction from scatter_superphoton, nscatt+1) onto ws. Escape/capture as in track_superphoton_grmhd.

grmhd_bias(bias0, theta_e, nscatt) -> Float64

Importance bias for the exponential-transform scatter decision on the GRMHD grid: b = max(1, bias0 (n_sc + 1)^2). bias0 > 1 boosts the rare scattering rate so the Compton tail is sampled at modest Ns; the (n_sc+1)^2 growth keeps higher orders populated. Unbiased for any b >= 1 (the child carries the reweighted w(1-e^{-dtau_s})/(1-e^{-b dtau_s}), see _mc_scatter_event). theta_e accepted for signature parity with the sphere bias_func.

mc_spectrum_grmhd(model::GRMHDModel, Ns; scattering=true, nthreads, seed,
                  incl_deg=(17.0,50.0,90.0), half_width_deg=10.0,
                  gen_nu_min=1e8, gen_nu_max=1e16, spec_nu_min=1e8, spec_nu_max=1e24,
                  n_esamp, n_ebins=160, C=80.0, R_record=1.0e3, max_scatt=3,
                  bias0=1.0, stride=1, r_emit_max=Inf, max_zones=typemax(Int))
    -> GRMHDSpectrumResult

Run the GRMHD Monte Carlo: build the emitting-zone list and grid weight table, distribute ~Ns seed superphotons over the zones (per-zone count ~ Ns, stochastic rounding), track each to escape/capture (curved geodesic, absorption always; Compton scattering when scattering=true), and bin the escaped photons by asymptotic inclination x log-nu. Threaded with one private accumulator set per chunk; each seed owns a Philox stream keyed by its global id, children a disjoint per-chunk stream. stride/`maxzonessubsample the grid for fast tests;remitmaxcaps the emitting region to the luminous torus.bias0boosts the rare scattering rate (unbiased for anybias0 >= 1);tabulatedgenopts into the tabulated frequency rejection (seemcspectrum`).

GRMHDSpectrumResult

Output of mc_spectrum_grmhd. spectra[k] is the single-row escaped Spectrum of photons whose asymptotic escape polar angle falls in inclination band incl_deg[k] +- half_width_deg; dOmega[k] is that band's solid angle 2 pi (cos th_lo - cos th_hi) (use nuLnu_solid(spectra[k], dOmega[k]) for the isotropic-equivalent nu Lnu). total is the angle-summed escaped spectrum. Diagnostics: `emittedpower(Sum w E m_e c^2 over generated seeds, erg/s),escaped_power` (over recorded photons), and the zone/photon counts.

incl_nuLnu(res::GRMHDSpectrumResult) -> (nu, L)

Per-inclination isotropic-equivalent nu L_nu [erg s^-1]: nu the shared bin centres, L[k, :] the spectrum for incl_deg[k] (via nuLnu_solid).