GOW Equilibrium

The first two tests compare JAFF against closed-form solutions. This one is different: there is no formula to check against, so instead it compares the generated network against an independent implementation — the equilibrium state of the Gong, Ostriker & Wolfire (2017) network, computed with a separate code and reproduced here.

The network

Files: networks/GOW/GOW.jet, GOW.jfunc, GOW.hdf5 — the network in KIDA format, its custom rate and heating/cooling functions, and the interpolation tables those functions read. See Predefined Networks for what each file contains.

It has eighteen species and about fifty reactions, and — unlike the other two tests — the gas temperature is part of what gets solved rather than something you hold fixed.

What it checks

The parts of the pipeline the two small tests do not touch:

• auxiliary functions — rate coefficients defined in .jfunc rather than as closed-form expressions (Auxiliary Functions);
• table interpolation — cooling functions that read grids out of GOW.hdf5, and their derivatives in the Jacobian (Table Interpolation);
• thermal coupling — chemistry and temperature solved together, each feeding the other;
• environment parameters — the radiation field, extinction, cosmic-ray rate, metallicity and dust-to-gas ratio reaching the rates correctly.

A network this size has no analytic solution, but it is still testable because its equilibrium state traces out a curve. The idea is to hold everything fixed except the density, integrate to steady state at each density in turn, and check that the resulting sequence of states follows the reference.

Setup

Environment

Reference tables are provided for all three radiation fields on a shared density grid. Running chi = 10 on its own already exercises the network; running all three adds a check that the field strength reaches the photorates with the right scaling, which a single value cannot show on its own.

Initial state

Start at a temperature of 100 K, with the abundances below given relative to hydrogen. Multiply each by the total hydrogen density n_H of the run to recover number densities in cm⁻³:

The scan

Repeat the run at a series of densities, holding everything else fixed:

The reference grid has 95 points spaced by a factor of 1.05, which is finer than the test needs — feel free to pick a subset, or use your own spacing and interpolate.

The end time is far longer than anything physical, and that is deliberate: it guarantees the endpoint is the steady state and not a snapshot part way there. A good way to confirm you have reached it is to check that the last decade of the integration changes nothing.

What to plot

Plot the final abundances against n_H on log-log axes, and put the gas temperature on a separate plot rather than a second y-axis. The species span roughly twenty decades, so it helps to split them: the ones above 1e-7 are what the comparison really rests on, and the rest are useful mainly as a diagnostic. The figures below show the reference tables — a natural way to compare is to overlay your own final states as markers on the same axes.

Temperature falls with density in all three radiation fields, offset roughly by the field strength:

The carbon balance is the clearest way to tell the fields apart. At chi = 0.1, neutral carbon overtakes C⁺ near n_H ≈ 600 cm⁻³; at chi = 1 it climbs to about a tenth of C⁺ by the top of the range; and at chi = 10 it stays at least two decades below throughout:

What to expect

With Av = 0 the gas is unshielded at every field strength, so CO never builds up to anything significant — the highest CO abundance anywhere in the three tables is only 6e-12. What does change with density is the hydrogen and carbon ionization structure, along with the temperature.

At chi = 10, between the ends of the density grid:

By the high-density end, the electron abundance has fallen to roughly the carbon abundance: most of the protons and He⁺ have recombined, which leaves C⁺ as the main carrier of charge.

The same quantities across the three fields, at the two ends of the grid:

A weaker field means colder gas, more H₂, and more neutral carbon, and the trend is monotonic across all three. If your run reproduces chi = 10 but not the ordering between the fields, that usually means the field strength is reaching some photorates but not others.

Suggested acceptance, per density point:

The band on the trace species is wide for a reason. Down at x ~ 1e-20 the reference itself carries noise from its own solver, and small differences in how the cooling tables are interpolated push those numbers around without any physical meaning. The carbon sum is a different kind of statement: it follows from the structure of the network, and it holds to better than 1e-7 relative in all three reference tables, so it should hold just as tightly for any integrator.

Reference data

benchmarks/GOW/equilibrium_chi0.1.txt, equilibrium_chi1.txt and equilibrium_chi10.txt — 95 rows each, with n_H running from 10 to 981.28 cm⁻³ in steps of a factor 1.05. They come from Munan Gong's PDR code with the parameters of Gong, Ostriker & Wolfire (2017). Since the three share a single density grid, they line up row by row. Rather than matching the reference spacing, interpolate the reference onto whatever grid your own run uses.

There are two quirks of the files to be aware of before comparing. First, the 2H2 column is twice the H₂ abundance, not the abundance itself. Second, a few entries for very rare species carry slightly negative values (of order 1e-18, all in the chi = 10 table); these are numerical noise from the reference solver, so it is safe to treat anything below about 1e-17 as zero.

When it fails