This page mirrors docs/athenak_rt.md in the repository, which is the authoritative version.

athenak_rt (the top-level athenak_rt/ directory) computes line-of-sight luminosities, spectra, image-plane maps, and light curves from AthenaK .bin dumps, for example of a TDE run. It is the post-processing described in the section "Radiative post-processing and synthetic observables" of Jiang et al. (2026) and follows the photospheric method of Yang et al. (2026, ApJ, 998, 118). Publications that use athenak_rt should cite both papers. The package reads a dump, resamples the finest AMR level of the dense gas onto a uniform Cartesian ray grid, recovers the gas temperature from the same tabulated equation of state (EOS) that the simulation used, and integrates the transfer equation along parallel rays. A run is described by a parameter file in the AthenaK input format.

This is a post-processing tool for a hydrodynamic calculation without radiation. Read Physics and limits before quoting luminosities. A shorter README is athenak_rt/README.md in the repository.

The package is pure Python and is not installed into the environment. Run it from the repository root, or put the repository root on PYTHONPATH.

python3 -m pip install numpy numba h5py      # add pytest to run the tests
cd /path/to/athenak
python3 -m athenak_rt --help

Data needed at run time.

The EOS table is looked up in this order (athenak_rt/eos.py).

1. <tables>/eos_table, if it is set to a path. The file must exist.
2. The path recorded in the dump header (<hydro>/table). It is tried as written (relative to the current directory) and, if it is relative, relative to the run directory (the parent of <input>/dump_dir) and to dump_dir. The example decks record ../chabrier2021_t13_helm_union_prad_640.table relative to the directory in which AthenaK was started, so the table is found while the dumps are still in that directory's bin/.
3. The basename of the header path inside $ATHENAK_EOS_TABLE_DIR.

If none exists the run stops with a list of the places it tried. Set <tables>/eos_table when the dumps have been moved.

Select a mode with <transfer>/mode. All modes share the same cell selection, the same ray grid, and the same BH excision mask. Luminosities are projected, isotropic-equivalent values for the chosen line of sight, L = 4 pi sum(I dA). For a blackbody pixel, L_bol = 4 sigma T^4 dA. The effective temperature of a pixel is T_eff = (pi I / sigma)^(1/4).

scattering = false turns multifreq into pure absorption. In all modes except tau1 the integration along a ray stops once the optical depth exceeds tau_stop (default 30). In multifreq mode this is checked per photon energy. The Rosseland tables are loaded in every mode, and in multifreq mode they are not used for the result. With populations = eos and a table without He or H2 fractions the run falls back to Saha with a warning, and the choice used is stored as populations_used in the output.

The multifrequency mode is the one adopted in the paper. The grey tau1 photosphere is kept as a reference. In the paper it overestimates the luminosity by roughly a factor of seven, because near the Planck peak the bound-free opacity is much larger than the Rosseland mean and the emission emerges from cooler, outer layers.

• Density cut. Cells with rho > factor * dfloor enter the ray grid (strict inequality, <selection>/density_threshold_factor, default 10), or cells above an explicit density_threshold_code. Other cells are treated as empty.
• Ray-grid box. By default the box is the bounding box of the selected cells, padded by the larger of auto_box_padding_code (0.25) and auto_box_padding_fraction (0.10) times the box width on every side, with a minimum width of auto_box_min_width (1.0), and clipped to the mesh (auto_box_clip_to_mesh). For the bounding box the threshold is factor * max(hydro/dfloor, problem/rho_floor) from the header. With the example decks both floors are equal. <image>/box fixes the box explicitly.
• BH excision. The sphere of radius problem/bh_excise_radius from the header is emptied around the BH position stored in the header. <selection>/bh_mask_radius changes the radius and bh_mask = false switches the mask off.
• Directions. <transfer>/directions lists one or more of x, y, z, -x, -y, -z, and the observer sits on that face of the box. The image rows and columns are (y, x) for z rays, (z, x) for y rays, and (z, y) for x rays.
• Resampling. The finest AMR level available at each location is sampled by nearest-cell lookup onto image_size x image_size x los_steps points. The dump is streamed block by block and is never held in memory in full.

python3 -m athenak_rt PARAMETER_FILE

The program takes one argument, the parameter file. Apart from --help and --version there are no command-line options. Two annotated examples that list every key are in athenak_rt/examples/: tde_snapshot.rtin (one dump) and tde_lightcurve.rtin (a series of dumps). The smallest useful file is

<input>
dump_dir = /path/to/run/bin
dumps    = 327

The format is that of AthenaK input files.

• A block starts with its name in angle brackets on a line of its own, for example <transfer>, and holds the key = value lines up to the next block. A block may appear more than once.
• # starts a comment anywhere on a line, also after a value, so a path that contains # cannot be given. Blank lines are ignored.
• Every key that is left out takes the default listed below. Only <input>/dump_dir and <input>/dumps are required.
• An unknown block, an unknown key, a key given twice, and a value that cannot be read stop the program with the file name, the line, and the closest valid name, so a typo is never ignored silently.
• auto selects the derived value of the keys whose default is auto.
• Booleans are true or false (yes/no, on/off, 1/0 are accepted). Lists (directions, box, dumps) are separated by commas or blanks.
• Relative paths are taken relative to the directory of the parameter file, not the current directory. ~ and $VARIABLES are expanded.

AthenaK writes its dumps as <basename>.<variable>.<NNNNN>.bin into the bin/ directory of the run. <input>/dump_dir is that directory, basename and variable complete the file name (defaults TDEExternalLTEPrad and hydro_w), and dumps selects the numbers.

The selection is sorted and duplicates are removed. The program checks that every selected file exists before it computes anything. With missing_dumps = error (the default) a missing dump stops the run, and the message lists the missing numbers. With missing_dumps = skip the missing numbers are printed as a warning and the other dumps are processed. If none of the selected dumps exists the run always stops.

Every selected dump is processed for every direction. The run also writes a light curve when at least two dumps qualify (Sec. 4).

The program prints one line per dump and direction with the bolometric luminosity and the name of the output file. Products that are reused with skip_existing are listed in the same way, marked (reused).

Each run writes into <output>/dir.

For example, the file for the -y view of dump 327 is TDEExternalLTEPrad.hydro_w.00327.rt_multifreq_-y.h5.

The bands are nir_J, nir_H, nir_K (1250, 1650, 2200 nm), optical_u, _g, _r, _i (355, 475, 622, 763 nm), uv_UVW1, uv_UVM2, uv_UVW2, uv_FUV150 (260, 224.6, 192.8, 150 nm), and xray_0p3keV, xray_1keV, xray_2keV (athenak_rt/bands.py).

rt_lightcurve_<mode>.csv is a plain CSV file with a header line and one row per dump in dump order. It has no comment lines. It contains every qualifying dump of the output directory, see below.

In column names the direction -y is written minus_y, for example L_bol_minus_y and nuLnu_optical_g_minus_y. The time in units of P_mb is not written, because the package does not know P_mb. Divide time_code by the P_mb of the run (88.04 code units for the example decks).

rt_lightcurve_<mode>.h5 holds the same rows as arrays.

Both files are rewritten after each dump, so a light curve can be inspected while a long run is going. They are assembled from the per-dump products and hold every dump, not only those of the current run. A product enters the light curve if it is in the output directory, has the same file-name prefix (<basename>.<variable>) and mode, exists for all directions of the run, and passes the test of skip_existing (below): the same settings and tables and, where the dump is in dump_dir, the same time and cycle. Rows are ordered by dump number. The files are written once at least two dumps qualify, so a run of one dump next to earlier matching products updates the light curve, and a run of one dump in an empty directory writes none. A light curve that no longer matches the current settings is not deleted, but it is not updated until two dumps qualify.

Every product is first written under a temporary name ending in .part and renamed when it is complete, so an interrupted run leaves no truncated product. With <run>/skip_existing = true a (dump, direction) whose product exists is not computed again. A product is reused only if it matches the current run, and the run stops before computing anything if an existing product does not, so a product made with other settings or from another dump is never reused silently. The comparison covers the following.

• The settings that enter the product of the mode. The output directory, the thread count and the compression are never compared. Settings that the mode does not read are ignored.

  Mode	Not compared
  tau1	tau_stop, nfreq, emin_ev, emax_ev, scattering, populations
  grey	tau_photosphere, nfreq, emin_ev, emax_ev, scattering, populations, the bands
  grey-therm	tau_photosphere, nfreq, emin_ev, emax_ev, scattering, the bands
  multifreq	tau_photosphere, the MESA tables

• The input dump: snapshot_file, snapshot_index, time_code and cycle stored in the product against the dump of the current run. A new dump_dir whose dump has the same name but another time or cycle is refused.

• The table files, by SHA-256 (eos_table_sha256, mesa_*_sha256), so moving the repository or the table does not invalidate a product. Products from versions that did not store the hashes are compared by path as before.

Reused products are reported on the command line and enter the light curve as if they had been computed.

import h5py, numpy as np

with h5py.File("rt_out/TDEExternalLTEPrad.hydro_w.00327.rt_multifreq_z.h5") as f:
    print("t =", f.attrs["time_code"], " L_bol,iso =", f.attrs["luminosity_bolometric_erg_s"])
    e_ev = f["spectra/energy_ev"][:]
    nu_lnu = f["spectra/nu_lnu_erg_s"][:]
    print("peak of nu L_nu at", e_ev[np.argmax(nu_lnu)], "eV")
    teff = f["maps/effective_temp_K"][:]          # (rows, columns), NaN where I <= 0
    print(f["parameters/resolved"].asstr()[()])   # the settings of this product

lc = np.genfromtxt("rt_lightcurve/rt_lightcurve_multifreq.csv", delimiter=",", names=True)
t_pmb = lc["time_code"] / 88.04                   # P_mb of the example decks
print(lc["dump"], t_pmb, lc["L_bol_z"], lc["L_bol_minus_y"])

The paper's main synthetic-observable result is the multifrequency H/He continuum transfer of the last high-resolution (HR) dump, 327 at t / P_mb = 2.93, viewed along +z (face-on) and along -y (in-plane, with the nozzle on the near side), on 1024^2 pixels with 1024 samples per line of sight. The HR run is not part of the example decks, so point dump_dir at a run of your own. The example FID chain of inputs/TDE_examples/README.md produces dumps that can be processed the same way. This is athenak_rt/examples/tde_snapshot.rtin.

# athenak_rt parameter file: one dump, the settings of the paper's Fig. (last HR
# snapshot, t/P_mb = 2.93, multifrequency transfer seen from +z and from -y).
#
# Copy this file to the directory that should receive the output, set
# <input>/dump_dir, and run from the AthenaK repository root (or with the root on
# PYTHONPATH)
#     python3 -m athenak_rt /path/to/work/tde_snapshot.rtin
# Relative paths are relative to this file's directory; ~ and $VARIABLES are
# expanded.  A key that is left out takes the default given in its comment;
# "auto" asks for the derived value.

<input>
dump_dir      = /path/to/run/bin    # directory with the .bin dumps (required)
basename      = TDEExternalLTEPrad  # files are <basename>.<variable>.<NNNNN>.bin
variable      = hydro_w             # output variable of the dumps (needs dens and eint)
dumps         = 327                 # N | N1, N2, ... | A:B | A:B:STEP (inclusive) | all
missing_dumps = error               # error | skip (warn and process the others)

<output>
dir              = rt_out           # output directory (default athenak_rt_output)
hdf5_compression = lzf              # lzf | gzip | none
bands            = true             # NIR/optical/UV/X-ray band L_nu

<transfer>
mode            = multifreq         # tau1 | grey | grey-therm | multifreq
directions      = z, -y             # lines of sight among x y z -x -y -z (default z)
tau_photosphere = 1.0               # tau1: optical depth of the photosphere
tau_stop        = 30.0              # grey*, multifreq: stop a ray past this depth
nfreq           = 73                # multifreq: log-spaced photon energies
emin_ev         = 0.1               # multifreq: lowest photon energy [eV]
emax_ev         = 1000.0            # multifreq: highest photon energy [eV]
scattering      = true              # multifreq: thermalization depth (false: pure absorption)
populations     = eos               # eos (EOS-table fractions) | saha (ideal Saha, X=0.7, Y=0.3)

<image>
image_size                = 1024    # pixels across the image plane
los_steps                 = 1024    # samples along the line of sight
grid_dtype                = auto    # auto | float32 | float64 (auto: float32 from 1024 samples)
box                       = auto    # auto | xmin, xmax, ymin, ymax, zmin, zmax (code units)
auto_box_padding_code     = 0.25    # auto box: absolute padding per side
auto_box_padding_fraction = 0.1     # auto box: padding per side as a fraction of the width
auto_box_min_width        = 1.0     # auto box: minimum width
auto_box_clip_to_mesh     = true    # auto box: clip to the mesh

<selection>
density_threshold_factor = 10.0     # cells with rho_code > factor * dfloor enter the grid
density_threshold_code   = auto     # auto | explicit threshold in code units (overrides the factor)
bh_mask                  = true     # empty the BH excision sphere
bh_mask_radius           = auto     # auto (<problem>/bh_excise_radius of the dump) | radius in code units

<tables>
eos_table   = auto                  # auto: <hydro>/table of the dump header, also tried relative
                                    # to the run directory (parent of dump_dir) and in
                                    # $ATHENAK_EOS_TABLE_DIR; or a path, e.g.
                                    # /path/to/chabrier2021_t13_helm_union_prad_640.table
mesa_high_t = auto                  # auto: athenak_rt/data/gs98_z0.02_x0.7.data
mesa_low_t  = auto                  # auto: athenak_rt/data/lowT_fa05_gs98_z0.02_x0.7.data

<run>
threads       = auto                # numba threads (auto: all CPUs available to the process)
skip_existing = false               # true: keep per-dump outputs that already exist

Copy the file to the directory that should hold the results, set dump_dir (and eos_table if the dumps have been moved away from the run directory), and run

cd /path/to/athenak
python3 -m athenak_rt /path/to/work/tde_snapshot.rtin

The products are rt_out/TDEExternalLTEPrad.hydro_w.00327.rt_multifreq_z.h5 and ..._-y.h5 next to the parameter file. For the grey tau = 1 photosphere as the reference, copy the file and change only mode = tau1. The products then end in .rt_tau1_z.h5 and .rt_tau1_-y.h5 and can share the output directory. For the last HR dump of the paper the published values are

The multifrequency spectrum peaks near 60 eV. The numbers depend on the dump, on the density threshold that separates debris from the numerical atmosphere, and on the image sampling. In the paper, going from 512^2 to 1024^2 pixels changes the grey luminosities by 0.4% and 6.1%, and doubling the number of photon energies changes the multifrequency luminosity by less than 3%. A quick test at image_size = 192 and los_steps = 192 takes a small fraction of the cost.

athenak_rt/examples/tde_lightcurve.rtin has the same keys. Its differences from the snapshot example are

<input>
dumps         = 300:327             # dumps 300 to 327 inclusive; 300:327:3 takes every third
missing_dumps = skip                # error | skip (warn and process the others)

<output>
dir              = rt_lightcurve    # output directory

<image>
image_size                = 512     # the paper's figure used 1024; 512 is 8 times cheaper
los_steps                 = 512     # samples along the line of sight

<run>
skip_existing = true                # reuse per-dump outputs made with the same settings

The run processes every dump from both directions and writes rt_lightcurve/rt_lightcurve_multifreq.csv and .h5 besides the 56 per-dump products. Resubmitting the same parameter file after an interruption computes only the products that are missing. The paper does not present light curves, because the simulation does not reach the peak of the flare.

import athenak_rt

results = athenak_rt.run_parameter_file("tde_snapshot.rtin")    # as the command line
print(results[0].luminosity, results[0].hdf5_path)

config = athenak_rt.load_parameter_file("tde_snapshot.rtin")    # parse and check only
print(config.settings, config.directions, config.find_dumps())

athenak_rt.run(settings, dumps, log=print, *, directions=None, skip_existing=False, provenance=None) runs an athenak_rt.RTSettings object on a list of dump paths (directions defaults to the one in settings), and athenak_rt.process_snapshot handles one dump and one direction. The kernels in athenak_rt.transfer accept any (nz, ny, nx) cube of density in g/cm^3 and temperature in K.

• The dump is streamed and never loaded in full. The dense arrays are the two resampled cubes (density and temperature) plus a small level cube. At 1024^3 with float32 storage this is about 4.3 GB for each of the two cubes plus 1 GB for the level cube. With float64 (used by auto below 1024 samples) the cubes double. The excision mask is applied on a small sub-cube so that it does not add a full-size temporary array.
• Plan for about 10 GB of host memory at 1024^3 and about 2.5 GB at 512^3 in double precision, plus the dump buffers. grid_dtype = float32 halves the cube memory below 1024 samples, at the storage precision used for the paper's production size.
• The cost grows with the number of pixels times los_steps, and for multifreq with the number of photon energies. Each direction is a separate pass over the dump. The multifreq kernel distributes the image rows over the numba threads, and the tau1 and grey kernels distribute the pixels, so <run>/threads should match the cores available. The first call includes compilation.
• Reading a large dump is limited by the file system, because the density and internal energy of every MeshBlock have to be read. Copy the dumps to fast storage if necessary, and point dump_dir there.
• Output files are moderate in size. A 1024^2 map takes 8 MB before compression. For tau1 the per-pixel band maps add 14 bands times that. bands = false skips them.

The tests are in athenak_rt/tests and need pytest, numpy, numba, and h5py. They do not need an EOS table or a simulation dump.

cd /path/to/athenak
python3 -m pytest athenak_rt/tests

Running pytest athenak_rt/tests from the repository root also works, because conftest.py adds the repository root to the path. The tests cover analytic uniform and two-layer slabs (formal solution, tau-stop truncation, the tau = 1 photosphere, all six directions), the multifrequency kernel against per-frequency analytic slabs with and without scattering, the grey-limit identity, Saha charge neutrality and ionization limits, the continuum edges, the Planck quadrature, and end-to-end runs on a small synthetic two-level AMR dump with a synthetic EOS table. The parameter-file tests check the parser and its error messages, that the defaults equal those of RTSettings, the dump selection (lists, ranges, strides, all, missing dumps), the example files, a single-dump run, a three-dump light curve with two directions, and the restart with skip_existing.

The results are LTE post-processing estimates, and the following limits apply.

• No radiation in the simulation. The hydrodynamic calculation is adiabatic and omits radiative cooling. athenak_rt does not evolve radiation energy or momentum, does not enforce a radiative energy budget, and does not include photon trapping or advection or radiative feedback on the dynamics. When the emitted power is comparable to the dissipation rate, as at the epoch of the paper's example, cooling would not be negligible. A prediction of the emergent luminosity at these super-Eddington rates needs radiation hydrodynamics. The luminosities that the paper quotes for its HR snapshot (about 25 to 40 Eddington luminosities for M_BH = 1000 Msun) are illustrative.
• LTE, continuum only. The temperature comes from the LTE EOS table, and the source function is a (modified) blackbody. There are no lines, no metals, and no non-LTE source functions. Scattering enters only through the thermalization-depth approximation of a homogeneous scattering atmosphere applied locally along each ray, and not through a full scattering solution.
• Composition. The grey modes use MESA opacities for X = 0.7, Z = 0.02, which is inconsistent with the metal-free H/He EOS of the hydrodynamics. The multifrequency mode takes its populations from the EOS table, which is consistent with the hydrodynamics. With populations = saha it assumes X = 0.7, Y = 0.3 (fixed in athenak_rt/constants.py), so do not use that option with a table of another composition without changing these constants.
• Debris versus atmosphere. The snapshots contain no passive scalar that separates debris from the numerical atmosphere, so the density threshold does this job. The luminosity depends on it, most strongly in the in-plane view.
• Rosseland photosphere. The grey tau1 photosphere is scattering inclusive, so it is not the thermalization surface. It overestimates the luminosity compared with the frequency-dependent calculation.
• No relativistic effects. Doppler boosting, aberration, and light-travel time are neglected. In the paper's example snapshot the maximum gas speed is about 0.03 c.
• Projection. The luminosities are isotropic-equivalent values for one line of sight, and they are view dependent.
• Not included. Photoionization codes and line emission, and the native-AMR streaming grey solvers, are not part of the package.

The photospheric method follows Yang et al. (2026, ApJ, 998, 118) (arXiv:2510.25547), adapted to Cartesian AthenaK output, and the multifrequency transfer is that of Jiang et al. (2026). Publications that use athenak_rt should cite both papers. Opacities come from MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019), including OPAL tables (Iglesias & Rogers 1996) and the low-temperature tables of Ferguson et al. (2005), for the Grevesse & Sauval (1998) mixture. Credit these sources when you use results derived from the bundled tables (see athenak_rt/data/NOTICE).