Using Radix

Radix solves the ionization, heating and radiative transfer of a photoionized slab of gas, as XSTAR does (xstar, 1-dimensional, one pass). It reads the atomic data (atdb.fits) and the Compton table (coheat.dat) of an XSTAR installation ($HEADAS/refdata/atdb.fits, or ftools/xstar/data/ of the source), and does not need XSTAR to run.

A slab

using Radix

model = Model("atdb.fits", "coheat.dat")                  # the elements and processes (abundances = :xdef, the default table of XSTAR)

run = slab(model;
    density = 1e4,                # cm⁻³ (or `pressure = 1e-7` dyn cm⁻² for constant pressure)
    column = 1e21,                # cm⁻²: the zones go on until this hydrogen column
    logξ = 1.0,                   # log of L/(n r²) (of L/(4π c P r²) for the pressure)
    luminosity = 1e46,            # erg s⁻¹ between 13.6 eV and 13.6 keV
    α = -1.0,                     # the source is E^α (or `spectrum = (E, L) -> blackbody(E, 0.5, L)`, `thermal`, `tabulated`)
    equilibrium = true,           # iterate the temperature of each zone to the thermal equilibrium (`false` keeps `T`)
    T = 10.0, xee = 1.0,          # the first guess, in 10⁴ K, and of the electron fraction
    emult = 1.0, taumax = 5.0, steps = 3,    # the thickness of the zones (XSTAR's `emult`, `taumax` and `nsteps`)
    vturb = 0.0)                  # km s⁻¹: the widths of the lines

Model(atdb, coheat; abundances, multiplier, exclude, cfrac) takes an abundance table (composition(:angr), ... or a vector for the elements 1 to 30), multiplier = Dict(26 => 2.0) for factors, exclude for the elements that are left out (Li, Be and B by default, as the atomic data of XSTAR has no data for them) and the covering fraction. Under the hood slab is powerlaw (the source), source_distance, Base.iterate of a mixture (the zones and the transfer); they can be called on their own, and it takes the zones ((r, Δr)) as well.

What comes back

A named tuple:

• zones: one named tuple per zone, with r and Δr (cm), T (10⁴ K), xee and nₑ, ntot (the hydrogen density), heating, cooling and imbalance (the relative difference), fractions[k][i] (the fraction of ion i of the element k of model.mixture.Z), populations (of the levels), radiation (the field of the zone), opacity (total, continuum, emissivity, bremsstrahlung) and escape. The first zone has no thickness: it only gives the opacity from which the others are chosen.
• E (eV) and incident (erg s⁻¹ erg⁻¹), the source on XSTAR's energy grid, transmitted (incident × exp(-dpthcont), the column of the output of XSTAR) and r;
• spectrum: the spectrum after the slab (L), the diffuse emission of the zones (inward, outward, and with the continuum opacity inward_continuum, outward_continuum, the columns emit_inward and emit_outward of XSTAR, in erg s⁻¹ erg⁻¹ like the incident spectrum);
• luminosities (inward, outward): the luminosity that each line and recombination edge of each element adds (10³⁸ erg s⁻¹; the records are model.mixture.elements[k].rates);
• depths (the optical depths of the lines and edges), dpthc and dpthcont (the depth of the continuum with and without lines and free-free).

Temperatures are in units of 10⁴ K, energies in eV, as in XSTAR; the spectra are in erg s⁻¹ erg⁻¹ (XSTAR's files are in units of 10³⁸ erg s⁻¹).

Results as arrays, and one file

XSTAR writes a dozen files (xout_abund1, xout_cont1, xout_spect1, xout_lines1, xout_rrc1, the xo_detal files, a log, ...). Radix returns everything in memory, as the named tuple above, so that it can be called repeatedly as a model (a fit, for instance) without touching the disk. tables turns the run into Results, four tables of plain arrays, one for each of the essential outputs:

t = tables(model, run)
t.zones       # radius, depth, thickness, ion_parameter, x_e, n_h, n_e, pressure, temperature, heating, cooling, heat_error and one column per ion (h_i, h_ii, he_i, ...)
t.spectrum    # energy, incident, transmitted, transmitted_lines, emit_inward, emit_outward, diffuse_inward, diffuse_outward, lines_inward, lines_outward (erg/s/erg, on XSTAR's energy grid)
t.lines       # ion, lower, upper, wavelength, emit_inward, emit_outward (erg/s), depth_inward, depth_outward
t.edges       # ion, level, energy, emit_inward, emit_outward (erg/s), depth_inward, depth_outward

write("model.fits", t; meta=(; DENSITY=1e4, COLUMN=1e21)) writes the four tables to one FITS file (an empty primary HDU with the meta keywords and the extensions ZONES, SPECTRUM, LINES and EDGES, with units); it passes HEASoft's ftverify. Nothing is written unless it is asked for. The spectrum has the pieces of both spectral files of XSTAR: xout_cont1 has transmitted, emit_inward and emit_outward (the continuum), and xout_spect1 has transmitted_lines and the emission diffuse + lines (the lines binned with their Voigt profiles, binned=false skips them: they take 10 s). Not in the tables: the detailed level populations (run.zones[k].populations keeps them in memory).

Agreement with XSTAR

Radix has been compared with XSTAR 2.59j for the rates of every record type, the ionization balance, the heating and cooling, the transfer (escape, optical depths, continuum opacity, lines, two-photon continua, diffuse emission), the thickness of the zones, the sources and constant pressure (the field notes have the details and the numbers, test/reference/ the reference output). For the thick slab at log ξ = 1 it reproduces the zones of XSTAR (depths to 10⁻⁴), the temperatures to 0.05%, the ion fractions to 10⁻³ or better and the diffuse emission to 5%. What differs:

• XSTAR ends its temperature iteration early in some zones (its dsec stops when the false position repeats a temperature): Radix finds the equilibrium, which is 0.4% to 2% from XSTAR's temperature in those zones.
• The atomic database of the XSTAR package ($HEADAS/refdata/atdb.fits) is not the one of its source tree: hydrogen has other records, and the temperatures of the tree are 3-6% lower than the package's. Use the one that you want to compare with.
• With vturbi > 0 XSTAR 2.59j smooths its continuum so that it is not absorbed below 20 keV (gsmooth2): Radix does not, and the comparisons with it use vturbi = 0.
• To reproduce the rounded constants of XSTAR (0.861707 eV for k × 10⁴ K, 12.56 for 4π, ...) use constants!(ucalc()) first; the default is CODATA 2022.
• The covering fraction cfrac cuts the photoexcitation of the lines, the escape probabilities and Thomson scattering, as in XSTAR (it agrees to 0.1% in temperature for cfrac = 0.5).
• The blackbody and thermal bremsstrahlung sources agree with XSTAR with the temperature iterated: the first zone of a blackbody (kT = 431 eV) to 0.03% in T and 0.3% in the ion fractions, and of bremsstrahlung (kT = 86 eV) to 0.4%. Through a slab of 10¹⁹ cm⁻² the first four zones of the blackbody (T from 11.8 to 7.6×10⁴ K) and all of the bremsstrahlung agree to 0.4% in T and 1% in the fractions. The blackbody slab then reaches a thermal front, where T falls to 1.8×10⁴ K in one zone: T agrees to 0.3% in the first zone of it and 1.7% at the end, and the trace ions of C to Fe differ by factors. XSTAR is not converged there (its T at the same radius changes by 1% with the size of the steps, it needs 20-30 iterations per zone, and its final row has a heating imbalance of 0.6%), and Radix changes by 0.8% in the same way.
• At a density of 10⁹ cm⁻³ the temperature of a zone at log ξ = 2 is 1.2% below XSTAR's (19.74 against 19.98, in the zone where the heating and cooling of XSTAR balance to 5×10⁻⁵); at 10⁴ cm⁻³ the same case agrees to 0.02%. It is the aluminium, and in it the helium-like Al XII: the heating minus the cooling of every other element agrees with XSTAR's table to a few 10⁻⁹. The rates of all the records of aluminium agree with XSTAR's ucalc at that cell (test/reference/ucalc/generate.jl with CONDITION= samples them), the photoexcitation rates agree to 0.02%, and the matrix of Radix has the same entries as the one that msolvelucy prints at lprint=5 (both count the identical pairs of CollisionAPED records of Al XII that the package database has, and the 36 stale rates that ucalc leaves between levels of equal energy change nothing). What differs are the populations: XSTAR's 2³P levels are 8% lower than Radix's (4.81×10⁻¹¹ against 5.20×10⁻¹¹ for J = 0), so the three lines of the 2³S–2³P triplet (900, 943 and 954 Å), which are 80% of the line cooling of aluminium, are 8% weaker. XSTAR's populations are not a steady state of its own matrix: one sweep of the relaxation that msolvelucy runs after the solution of its condensed (one level for each ion) matrix, from XSTAR's final populations, changes the 2³P levels by +10%. That relaxation stops when the populations above 10⁻⁶ change by less than 1% (crit2 = 10⁻²), so the small populations that depend on a slowly converging level are left a sweep behind. Radix solves the whole matrix, and its populations satisfy it to 10⁻¹³.
• passes=3 (XSTAR's npass) repeats the march with the optical depths of the lines and edges between each zone and the outer edge, as the escape probabilities of the outer side. It changes the ion fractions by up to 2.5% on the thick slab and follows XSTAR's change to 0.15%. passes must be odd.

Speed

Radix runs the elements of the gas in parallel, so start Julia with several threads (julia -t auto): the default is one. Wall times on an Apple M-series laptop with the database of the XSTAR source tree (the first call also compiles):

The zones are not the same in the second row (the thickness rule of Radix gives fewer), so compare the time of a zone: 41 s for XSTAR, 20 s and 11 s for Radix. passes=3 takes about twice as long as passes=1, XSTAR's npass=3 four times. Most of the time is the integration of the photoionization cross sections over the radiation (once for each balance of the elements, the heating, the opacities and the emissivities of each zone).

Running the tests

Pkg.test("Radix") runs the tests of the rates and the toy models. The comparisons with XSTAR (the level balance, heating, transfer, spectra and the whole database) need its atomic database, which is too large for the repository (875 MB): give the path of atdb.fits in RADIX_ATDB (the one of the XSTAR 2.59j package, $HEADAS/refdata/atdb.fits, or of its source tree: both pass) and use several threads:

RADIX_ATDB=$HEADAS/refdata/atdb.fits julia -t 4 --project -e 'using Pkg; Pkg.test()'

It takes about 7 minutes with 8 threads. RADIX_QUICK=1 skips the five tests that take minutes (the thermal equilibrium of zones, the thick slab with its emission lines, the three passes with the temperature iterated and the model loaded from the files): the rest take about 3 minutes. The continuous integration runs that quick suite: it downloads the database from HEASARC once (https://heasarc.gsfc.nasa.gov/FTP/software/lheasoft/xstar/atdb.fits) and keeps it in the cache of the workflow, and the whole suite is run locally.