Time coordinates, transport, and energy

The numerical calculation discretises the physical model without adding fitted temperatures, band amplitudes or empirical time delays. Its checks distinguish conservation, iterative convergence and resolution sensitivity.

1. Computational sequence

1. Validate configuration and the canonical 22-epoch table; load the selected STARS history, optical properties, WISE responses and Vega spectrum.
2. Build the scaled and smoothed luminosity history, then conservatively average it into transport time bins.
3. Construct finite-volume thermal cells and grain numbers from the radial optical depths.
4. Build direct-illumination, secondary-absorption and escape kernels with the same opacity geometry.
5. Solve the coupled grain equilibrium and delayed secondary-heating problem by iteration.
6. Check grain temperatures, probability closure, thermal closure, energy accounting and the residual tail.
7. Calculate directional observer transfer on the WISE wavelength nodes, add transmitted source light and apply Vega calibration.
8. Profile the event epoch against the two-band photometry and save the configuration, provenance, diagnostics and response.

2. Cells, sizes and wavelength quadrature

Each occupied radial population is split logarithmically in radius. Polar cells are uniform in μ = cos θ across −cos θ_open to +cos θ_open. Empty radial gaps are retained in the geometry but do not become thermal cells. The model is axisymmetric in its thermal state, not a fully independent three-dimensional temperature grid.

For radial boundaries r_k,− and r_k,+, the power-law shape is volume-averaged:

Lengths are converted consistently to physical units when the column is normalised. Density is constant within a cell; using its value at the radial midpoint instead would be a different discretisation. Cell volumes and grain counts are:

The 41 grain-size nodes are logarithmic, but w_ℓ is not a uniform weight. It integrates a^−q between logarithmic cell boundaries, with the first and last boundaries fixed at a_min and a_max; the weights sum to one. Absorption cross-sections already contain πa², so that area must not be included a second time in w_ℓ.

The transport grid spans 0.01–1000 μm. Trapezoidal wavelength quadrature is included in source spectral fractions and emitted bin powers. Those binned quantities are summed directly; multiplying them by another wavelength interval would double-count the quadrature.

3. Time bins and conservative driver integration

The driver is sampled at 10-day rest-frame resolution independently of the radiative-transfer grid. A piecewise-linear antiderivative is integrated across each transport bin to obtain its mean luminosity:

The integration limits are clipped to the supplied driver support. The code checks the sum of bin energies against the retained driver integral to relative precision 10⁻¹⁰. This preserves the supplied energy; it does not rescale the model to a target energy. The full-driver energy and the finite retained-window energy are recorded separately.

The stored time is u = t_local−r/c. Radially propagating direct source light has zero extra lag in u, even though its travel time in local coordinates is r/c. This avoids an acausal direct front caused by mixing different radii within a finite cell. Secondary paths have Δu = (|r₂−r₁|+r₁−r₂)/c ≥ 0; observer paths have delay (r−n̂·r)/c before multiplying by 1+z.

4. Three-dimensional ray kernels

path_segments finds intersections with the spherical radial boundaries and conical polar boundaries. Midpoint lookup identifies which intervals contain dust. Vacuum intervals contribute flight time but no optical depth. For each occupied segment:

Emission positions are sampled uniformly in cell volume. Scrambled Sobol sequences provide deterministic, seeded quasi-Monte Carlo quadrature, and thermal directions are isotropic. Axial symmetry allows the emitting azimuth to be fixed for the all-direction transport calculation; the observer calculation samples full positions. These samples represent a continuous medium, not discrete fitted clouds.

Absorption within a uniform segment is placed at its conditional mean distance, using the wavelength-dependent optical depth. A lag between adjacent time nodes is deposited with non-negative linear weights that sum to the original probability. This preserves probability but approximates the full within-segment flight-time distribution, so time and spatial resolution still require testing.

For every source cell and wavelength, summed absorption plus escape probability must equal one to within 10⁻¹¹. The direct-source fractions are checked separately. The observer operator is directional and is not interchangeable with the all-direction escape operator.

5. Grain equilibrium and nonlinear iteration

The solver starts with directly absorbed power. At each iteration it partitions total absorbed spectral power among sizes, finds each equilibrium temperature and reconstructs the full thermal emission. Schematically, using convolution over non-negative lags:

H is absorbed cell power and P is emitted cell power in a wavelength bin. The emission depends nonlinearly on temperature, so this is a coupled causal Λ iteration, including zero-lag self-absorption. FFT convolution accelerates the delayed sums; padding prevents cyclic wraparound. Small numerical negatives are clipped and pre-source secondary roundoff is explicitly removed.

The active Thermal class tabulates the Planck power for every grain size on 2400 logarithmic temperature nodes from 0.5 K to max(5000 K, 2T_sub). Logarithmic interpolation supplies an initial temperature and two Newton corrections enforce equilibrium. It does not use the older helper’s temperature lookup as the production solution. Heating below the lower table threshold produces zero transient emission.

There is no post-hoc luminosity correction. Reconstructed single-grain emitted power must match its absorbed power with a maximum error, normalised by the largest power for that size, below 2×10⁻⁷. Every temperature is also checked against T_sub; a sublimation failure stops evaluation rather than modifying the opacity.

6. What “converged” means

At fixed discretisation, two changes in the absorbed field are monitored: the maximum absolute change divided by the maximum new absorbed power, and the L1 change divided by the sum of new absorbed power. Both must be below 10⁻⁶. The default limit is 60 iterations. Failure to converge does not produce a valid fit.

Archived Case A required eight iterations. Its final maximum relative change is 1.25×10⁻⁷, L1 change 7.88×10⁻⁷, and thermal closure residual 8.61×10⁻¹⁶. These are properties of the saved calculation, not results of a new execution for this documentation.

Three checks answer different questions:

• Probability and energy closure: does the discretised calculation conserve radiation?
• Iterative convergence: has the coupled solution settled on this fixed mesh?
• Resolution sensitivity: would finer radial, angular, time or wavelength sampling materially change the observables?

A tiny energy residual alone cannot establish a resolved light curve. Resolution comparisons should hold physical parameters fixed, state whether t₀ was re-profiled, and compare both fluxes and χ². Parameter reoptimisation is a further, distinct test.

7. Archived fine settings

The numerics block in the saved summary.json is authoritative. Legacy configuration keys such as n_radii, n_mu, n_phi and seed do not describe all active settings. The driver’s 10-day step and transport’s 30-day step are deliberately different.

8. Directional transfer and synthetic photometry

For each observer position sample, the code follows the entire escape ray, including inward-facing or opposite-side material, and records its grain column and geometrical delay. Its weight is exp(−column × C̄_λ). The emitted spectra are recomputed from the saved temperatures at the rest wavelengths corresponding to the native WISE response nodes, rather than approximating the bands by single wavelengths.

The caller applies distance, redshift, unit conversion and the Vega equal-energy band operator. A conversion factor of 10⁻⁷ changes F_λ from W m⁻³ to erg s⁻¹ cm⁻² Å⁻¹, matching the CALSPEC convention. The response-weighted target and Vega integrals share the same wavelength convention. The transmitted source spectrum is attenuated along the direct central ray and added without an independent amplitude parameter.

9. The finite-window energy ledger

All reported energies are in erg; internally the thermal transport uses binned powers in watts. The ledger distinguishes emission events from net energy escaping the system:

The working closure includes direct and thermal radiation pending absorption beyond the computed interval:

“Eventual escape” does not mean every photon has already arrived at the observer by the end of the stored time grid. Pending absorption terms describe radiation still in flight; further re-emission would require extending the thermal solution. The archived pending terms are zero after the implementation’s non-negative reporting convention, with a very small recorded convolution-roundoff term. Its relative global residual is −1.09×10⁻⁶.

The final ten steps contain approximately 5.02×10⁻¹⁰ of the intrinsic emitted energy. The solver rejects runs when that fraction exceeds 10⁻³, indicating an insufficient response window. Finite-window bookkeeping and full-driver truncation remain separate checks.

10. Saved state, visualisation and reproducibility

Runs save configuration, active numerics, file hashes and numerical-library versions. Matching fingerprints permit reuse of compatible transport kernels or completed output; a mismatched output fingerprint requires a new output directory. This distinguishes a physical parameter change from a changed input table, code version or numerical environment.

When thermal arrays are retained, they contain the time grid, individual grain temperatures, emitted bin powers, wavelength grid, grain counts and mesh geometry. A number-weighted displayed temperature is Σ_ℓN_kℓT_kℓ/Σ_ℓN_kℓ; it is not a fitted colour temperature or an energy-weighted temperature.

The website’s point cloud is a sampled export of the archived solution. At common local time, each point samples the thermal history at u = t_local−r/c. A source-retarded map, a simultaneous local snapshot and an observer arrival-time slice answer different questions and must retain their time labels. No new transport calculation is performed by the browser.

See the pinned implementation map and quick commands to reproduce the archive. Reproducing one case is distinct from rerunning the full search or validating every modelling assumption.