Using NRHJSur3dq8

This is the usage guide for the released surrogate waveform models NRHJSur3dq8_AA_v3 (action-angle) and NRHJSur3dq8_AdA_v3 (adiabatic-angle), together with the merger-ringdown artifact NRHJSur3dq8_MR_v3 that both of them use.

Everything on this page has a runnable counterpart in examples/usage/. The numbers quoted are outputs of those scripts on one machine.

Read the release names from nrhjsurrogate.RELEASE_AA and nrhjsurrogate.RELEASE_ADA rather than hard-coding one, and your code follows the release.

Contents:

1. Finding the model files
2. Waveform evaluation
3. Physical units and polarizations
4. Starting earlier than the artifact reaches
5. The three waveform routes
6. Parameter derivatives
7. The Fisher matrix
8. The model's own error estimate
9. Backends
10. The adiabatic-angle model
11. Dynamics and remnant
12. Things that catch people out

1. Finding the model files

The package ships code only. The model artifacts are HDF5 files distributed separately, and the package fetches them for you the first time you name a released model and no copy is found. It downloads from the deposit into your per-user cache, checks the SHA-256 against the registry, and moves the file into place only if it matches, so a failed or altered download is never left where the resolver would pick it up. After that first call it is a local file like any other.

You can turn that off. NRHJ_AUTO_FETCH=0 disables it globally, which is what an air-gapped machine wants, and resolve_model(..., fetch=False) does it for one call; either way you get the old error listing everywhere that was searched. Two things are never downloaded: an explicit path that does not exist, because that request named a file rather than a model, and a model that is not part of the public release, because there is nothing deposited to fetch.

You need the merger-ringdown artifact too. NRHJSur3dq8_AA_v3.h5 names NRHJSur3dq8_MR_v3.h5 and resolves it by name on every load. If it is not there, the load raises rather than returning a model that stops at the merger:

FileNotFoundError: could not locate the model artifact 'NRHJSur3dq8_MR_v3.h5'.
  file   : NRHJSur3dq8_MR_v3.h5
  size   : 19606296 bytes
  sha256 : dd1bf00860d01e16cf740bd227f5f7fb163284cb42309febd7194bc8d235059a
  ...
Searched, in order:
  - ...

Resolution order, first hit wins:

The cache directory is $NRHJ_MODEL_CACHE_DIR if set, else platformdirs' user cache if that package happens to be installed, else $XDG_CACHE_HOME/nrhjsurrogate, else ~/.cache/nrhjsurrogate.

$LAL_DATA_PATH sits below $NRHJ_MODEL_DIR so a project-specific choice beats a shared site install, and above the cache so a file you put somewhere on purpose beats one the package downloaded on its own. The package never writes into $LAL_DATA_PATH.

Inspecting and checking, without loading a model:

from nrhjsurrogate import RELEASE_AA, models

models.search_path()                       # the directories, in order
models.model_cache_dir()
models.resolve_model(RELEASE_AA)           # name -> path
models.verify_model(path, RELEASE_AA)      # the file's sha256
models.resolve_model(RELEASE_AA, verify=True)   # both, raises on mismatch

Read release names from RELEASE_AA and RELEASE_ADA rather than typing a file name, so your code follows the release.

models.fetch_model(name) downloads an artifact into the cache directory and checks its SHA-256 before moving it into place, so a truncated or altered download can never be left where a loader would find it. It downloads from the Zenodo deposit by default; set $NRHJ_MODEL_BASE_URL to point it at a mirror instead, which also accepts a local directory or a file:// URL.

Every environment variable this section mentions, in one place:

Those are read at RUN time. A separate set is read at INSTALL time only -- NRHJ_BUILD_KOKKOS, NRHJ_BACKEND, NRHJ_CUDA_ARCH, NRHJ_NVCC_WRAPPER, CMAKE_ARGS -- and they decide whether the compiled backend exists at all and which Kokkos spaces it carries. They do nothing if you set them here. INSTALL.md sections 1b and 5 are their home. To see what the install actually produced:

python -m nrhjsurrogate.compiled_backend

Earlier spellings of the first four still work and warn once, naming the variable to use instead. The new name wins when both are set.

Runnable: examples/usage/10_finding_the_artifacts.py.

Checking the install and the files

Once the files are in place, one command checks that the install and the model files work together, against reference values that ship inside the package:

nrhjsur-verify                      # or: python -m nrhjsurrogate.verify

It prints one PASS, FAIL or SKIP line per check (file integrity by SHA-256, load, waveform, parameter derivatives, the model's own error estimate, hybridisation, dynamics, remnant, physical units, the adiabatic-angle model) and exits nonzero on any failure. With no model files it skips every model check and says which file and which variable it needed. pytest --pyargs nrhjsurrogate.verify runs the same checks under pytest. Details and the tolerances are in INSTALL.md section 6.

2. Waveform evaluation

from nrhjsurrogate import NRHJSurAA

model = NRHJSurAA.load_h5("NRHJSur3dq8_AA_v3")
times, modes = model.get_td_waveform_modes_geometric(4.0, 0.3, -0.2)

There are four waveform methods, across two axes: what you want, and which units you want it in. The unsuffixed name is physical units and the _geometric suffix is geometric ones, and each takes only the arguments that mean something in its own units.

This section covers the geometric mode face; section 3 covers the physical ones. What comes back is

• times, the time samples, in units of the total mass M, with zero at the amplitude peak
• modes[(l, m)], the complex mode r h_lm / M, one entry per mode

The model is also still callable, model(4.0, 0.3, -0.2), which returns the same content as one dictionary with the times under "t". That spelling is superseded but retained, and it forwards to the same implementation, so the two cannot disagree.

The artifact carries seven modes with m > 0: (2,2), (2,1), (3,3), (3,2), (4,4), (4,3), (5,5). The modes with m < 0 follow from the reflection symmetry of an aligned-spin binary, h_{l,-m} = (-1)^l conj(h_lm).

At mass_ratio = 4, spin1z = 0.3, spin2z = -0.2 the model returns 15775 samples covering -3819.1 M to +124.4 M, with |h_22| peaking at 2.4156e-01.

Useful arguments:

with_derivatives is the one argument on this surface whose value changes what comes back: False gives times, modes and True gives times, modes, derivatives. Nothing else varies that way, and behind the method each of the two is a separate call with one fixed shape.

times is evaluated directly. The model does not build its own grid and resample onto yours, so a sparse grid is both faster and more accurate than doing the resampling yourself. Samples outside the model's support come back as exactly zero.

Arrays broadcast. Passing arrays of mass_ratio, spin1z and spin2z gives modes of shape (batch, samples) on one grid shared by the whole batch.

Parameter range. The model covers roughly 1 <= mass_ratio <= 8 and spin components from about -0.95 to +1.0. model.in_hull(mass_ratio, spin1z, spin2z) answers for a specific point. Outside, the prediction clamps to the boundary, which is an extrapolation, and you get an OutOfHullWarning.

Runnable: examples/usage/01_load_and_modes.py.

3. Physical units and polarizations

times, modes = model.get_td_waveform_modes(
    mass1=50.0, mass2=12.5, spin1z=0.3, spin2z=-0.2,
    delta_t=1.0 / 4096.0, f_min=30.0, distance=400.0)

h_plus, h_cross = model.get_td_waveform(
    mass1=50.0, mass2=12.5, spin1z=0.3, spin2z=-0.2,
    delta_t=1.0 / 4096.0, f_min=30.0, distance=400.0,
    inclination=np.pi / 3.0)

Units follow the LAL convention: masses in solar masses, delta_t in seconds, f_min and f_ref in Hz of the (2, 2) gravitational-wave frequency, distance in Mpc, inclination, azimuth and phi_ref in radians. times is in seconds with zero at the amplitude peak.

So does the frame. (inclination, phi_ref) here means what (inclination, phiRef) means to LAL's SimInspiralChooseTDWaveform: phi_ref is the orbital phase at f_ref (the modes are rotated so that arg h_22 = -2 phi_ref there, with the heavier body on +x), and the line of sight sits at azimuth pi/2 in the source frame, which is where LAL puts it (the wave-frame Y axis is the ascending node). azimuth is that harmonic azimuth and defaults to pi/2; you do not normally touch it. Measured against NRHybSur3dq8 through LALSimulation, the two models' polarizations agree at the same (inclination, phi_ref) to the level their modes agree. The polarizations of release 1.0.0 were built in the model's native source frame at azimuth 0; source_frame="native", azimuth=0.0 reproduces them bitwise.

The served MODES are in the model's native source frame on every face, geometric and physical alike. LAL's source frame is a quarter turn and a strain sign away: h_lm(LAL) = i**(m - 2) * h_lm(here) (nrhjsurrogate.driver.aa_api.lal_source_frame_factor), unity on the m = 2 modes, so a mode-by-mode comparison against SimInspiralChooseTDModes needs that factor and nothing else once phi_ref is matched.

The example above gives a 0.4241 s waveform from -0.385883 s to +0.038190 s, with |h_22| peaking at 1.8062e-21, and at inclination = pi/3 the polarizations peak at max |h_plus| = 9.0797e-22 and max |h_cross| = 7.5521e-22.

Two products, two methods. get_td_waveform_modes gives the modes and get_td_waveform gives the polarizations,

h_plus - i h_cross = sum_lm h_lm  {}_{-2}Y_lm(inclination, azimuth)

summed over all m, the reflection modes included. On the polarization method inclination and distance are required rather than optional, because a polarization without a distance is not a strain and without an inclination there is nothing to project onto.

get_td_waveform returns two arrays and no time axis, which is the pycbc spelling this name comes from, where the returned series carry their own epoch and spacing. Here they do not, so if you need the axis take it from get_td_waveform_modes at the same arguments; it is delta_t-spaced and peak aligned.

distance decides the scale on the modes face. With it, modes are strain at that luminosity distance. Without it they are the geometric r h_lm / M. The geometric face is get_td_waveform_geometric, which takes an inclination but no distance, because in geometric units there is nothing for a distance to mean. At mass_ratio = 4, spins (0.3, -0.2) and inclination = pi/3 it gives max |h_plus| = 1.1107e-01 and max |h_cross| = 9.5352e-02 over the same 15775 samples as section 2.

f_ref sets the phase convention only. The modes are rotated so that arg h_22 = -2 phi_ref at the sample where the (2, 2) frequency first reaches f_ref. The default f_ref is f_min; if both are absent the reference is the first sample. Aligned spins do not precess, so nothing else needs a reference frequency.

The mass arguments are order agnostic. The pair is sorted internally and each spin follows its own body, so (50, 12.5, 0.3, -0.2) and (12.5, 50, -0.2, 0.3) give the same waveform.

times works here too, in seconds and peak aligned, the same convention as the returned axis.

No with_derivatives on the physical face. The derivatives are with respect to mass_ratio, spin1z and spin2z, and this face rotates the modes to a reference frequency using the waveform itself, so a physical-units derivative carries a chain-rule term the geometric one does not. Rather than return something that looks like a derivative and is missing a term, the derivative product is on get_td_waveform_modes_geometric only.

The superseded model.physical(...) is retained and returns the whole result as one dictionary with "t", the modes, and "hp"/"hc" when an inclination is given.

For many parameter points on one shared physical grid there is model.physical_batch(...). It takes the same keywords and returns the rows zero-padded onto a common grid, plus a "mask" marking where each row actually has support, which is the shape a batched likelihood consumes. Its own docstring says it covers the artifact's span only, but measured behaviour is that a below-reach f_min extends as in section 4 (60 solar masses from 10 Hz returned a 5.94 s waveform against a 1.42 s native span). Do not rely on it raising.

Runnable: examples/usage/02_physical_units.py.

4. Starting earlier than the artifact reaches

The artifact carries a finite stretch of inspiral. What happens when you ask for an earlier start depends on whether it is already covered.

In geometric units the lever is f_low, the dimensionless M f of the (2, 2) mode. On the physical methods it is f_min, the same quantity in Hz.

The trim really is only a trim. At mass_ratio = 4, spin1z = 0.3, spin2z = -0.2 the samples returned for an in-reach f_low are equal, bit for bit, to the corresponding samples of the untrimmed waveform.

The extension is controlled by hybridization, which takes three named states:

• "extend_if_below_reach" (the default): trim when in reach, extend when below.
• "never": never extend. A below-reach request raises a ValueError that names the point's own start frequency, for example f_low=0.001 is below this point's native start (M f = 0.006311) and hybridization='never' guarantees NR-only content. Use this when you need to know the content is numerical-relativity based throughout.
• "require": an intent declaration, "I require an extension". It raises when the extension cannot be delivered, for instance when f_low is absent or already in reach.

These replace the older hybridize="auto", hybridize=False and hybridize=True, which were a string mixed with two booleans on one keyword where True did not mean "yes please" but "raise unless you actually extend". The old spellings still work, through either keyword, and warn.

Measured at that same point, whose own content starts at M f = 0.006311:

The result dictionary of model(...) and model.physical(...) always carries "hybrid_report" and "dynamics", and both are None when the call did not produce one. Pass dynamics=True and "dynamics" is filled in: {Phi, J, E, omega, t} over the native orbital-phase grid, or over the extended range when the extension engaged. "hybrid_report" is filled in only when the extension actually engaged: it records what the extension did (where the artifact's content starts and where the waveform now starts, how many orbits and how much time were added, how well the join matches). Both are scalar-input, NumPy-backend only, and neither costs anything unless dynamics=True is passed.

They used to be absent rather than None when they did not apply, which meant every reader had to write the membership test before the lookup, and the test is the part that gets forgotten.

Runnable: examples/usage/03_hybridisation.py.

5. The three waveform routes

One artifact serves three ways of building a waveform from its content, selected by timing_mode.

"Derived" means the clock is obtained from the stored action and energy content through the first law of binary black hole mechanics rather than read from a fitted clock. AAintPN goes one step further and takes the mode amplitudes from the post-Newtonian expressions instead of the fitted content.

The two derived routes are reference models. They cost accuracy, which is the point of having them: the gap between them and the default route is a measurement, not a defect.

Practical restrictions on the derived routes:

• inspiral only. They end at the plunge, and include_merger_ringdown=True is refused.
• one parameter point at a time. Arrays are refused.
• they carry their own extension, which is a different mechanism from section 4. Below the stored content the action and energy continue along their post-Newtonian description, and the clock is a single integration over the whole extended range, so there is no join in time. They report what they did under the "gate_report" key.

Measured at mass_ratio = 4, spin1z = 0.3, spin2z = -0.2:

With f_low = 0.003 both derived routes return 122179 samples from -30544.5 M, having added 52.75 orbits.

Runnable: examples/usage/04_derived_routes.py.

6. Parameter derivatives

times, modes, derivatives = model.get_td_waveform_modes_geometric(
    4.0, 0.3, -0.2, with_derivatives=True)

modes[(2, 2)].shape          # (N,)
derivatives[(2, 2)].shape    # (3, N): d/dmass_ratio, d/dspin1z, d/dspin2z

model.modes_grad(...) is the same call under its older name and is retained.

By parameter name. Nothing in a (3, N) array says which axis holds the parameters or which row is which, so every derivative face also has a _named sibling that returns the same numbers as a namedtuple per mode:

times, modes, named = model.modes_grad_named(4.0, 0.3, -0.2)

named[2, 2].spin1z           # d h_22 / d spin1z, shape (N,), same row as
                             # derivatives[(2, 2)][1]
named[2, 2]._fields          # ('mass_ratio', 'spin1z', 'spin2z')

The array form is retained unchanged. nrhjsurrogate.derivatives_by_parameter converts any face's array output on demand (parameter_axis=1 with the mode list for the batched (B, 3, K, N) arrays), and nrhjsurrogate.DERIVATIVE_PARAMETERS states the row order once. The siblings are modes_grad_named, hybrid_modes_grad_named, modes_grad_batch_named and, on the AdA class, modes_grad_named.

The waveform comes with the derivative, and that is not a convenience. The gradient kernel produces both in one fused pass, so asking for the derivative alone saves nothing: measured at mass_ratio = 4, spins (0.3, -0.2), 8192 samples, the waveform kernel is 84.10 ms and the fused gradient kernel is 406.43 ms, and a route that called one and then the other would pay the 84.10 ms twice over for a waveform it already had.

The derivatives are analytic. They come from the same stored coefficients as the waveform, so no finite differencing of the waveform is involved. Times are geometric and peak aligned. With the merger-ringdown arm attached, which is the default, the derivatives cover the same span the waveform does: at mass_ratio = 4, spins (0.3, -0.2), 15775 samples from -3819.1 M to +124.4 M. With include_merger_ringdown=False it is the inspiral alone, 15276 samples from -3819.1 M to -0.3 M.

times= takes your own grid, in M, exactly as in section 2. A mode subset and a time_step are refused rather than ignored when you ask for derivatives: the fused kernel produces all seven modes on its own grid, and silently dropping either argument would answer a different question from the one asked.

As a check anyone can repeat, a central finite difference of the waveform in q, on 4096 samples between -3000 M and -50 M, matches the analytic dh/dq on the (2, 2) mode to a median of 1.78e-05 of the peak of dh/dq, worst 1.17e-04 near the plunge. That residual does not shrink when the difference step is made smaller, so it is not a finite difference truncation error.

Restrictions:

• On the default route, a f_low request has no derivative path. It raises and names hybrid_gradients.hybrid_modes_grad as the alternative.
• timing_mode="aa-int-pn" accepts f_low, and timing_mode="aa-int" accepts an in-reach f_low only.
• Using the extended waveform and its derivative together is inconsistent by a stated amount, and the model warns and quotes it (1.7e-3 to 5.2e-3 of the peak, worst on the (4,4) mode). Treat that warning as information you have to carry into anything that consumes both, such as a Fisher matrix.

Derivatives with respect to the eight physical parameters, rather than the three intrinsic ones, come from the adiabatic-angle class (section 10).

Runnable: examples/usage/05_parameter_derivatives.py.

7. The Fisher matrix

The Fisher matrix is F_ij = <dh/dp_i | dh/dp_j> under a noise-weighted inner product, and its inverse is the leading-order parameter covariance. It is built from the analytic derivatives, so again nothing here finite differences a waveform.

import numpy as np
from nrhjsurrogate import RELEASE_ADA
from nrhjsurrogate.driver.api import NRHJSurAdA
from nrhjsurrogate.driver.fisher import (PARAM_NAMES, InnerProduct,
                                         PhysicalModel, fisher_covariance,
                                         fisher_matrix)

model = NRHJSurAdA.load(RELEASE_ADA, backend="numpy")
strain_model = PhysicalModel(model.insp_ev, model.mr_ev,
                             f_plus=1.0, f_cross=0.0)

params = np.array([60.0, 2.0, 0.1, 0.0, 400.0, 0.0, 0.3, 1.1])
grid = np.arange(-1.0, 0.03, 1.0 / 4096.0)
inner_product = InnerProduct(grid, f_min=20.0)

fisher = fisher_matrix(strain_model, params, grid, inner_product)
covariance = fisher_covariance(fisher)

PARAM_NAMES is (mtot, q, chi1z, chi2z, dist, t_c, phi_ref, incl): total mass in solar masses, distance in Mpc, t_c and the grid in seconds, phi_ref and incl in radians. f_plus and f_cross are the antenna response, held fixed.

InnerProduct defaults to an Advanced LIGO zero-detuned high-power design noise curve; pass your own callable as psd.

The grid must stay inside the waveform. If it does not, you get geometric times leave [inspiral start, ringdown end] with both spans printed, rather than silent zero padding.

At the parameters above the 1-sigma widths are

free= restricts the matrix to a subset of parameters, holding the rest fixed. With only (mtot, q, chi1z, chi2z) free the same point gives 0.6227120, 0.08473675, 0.1017546 and 0.2528099.

The physical-parameter Fisher matrix is strongly degenerate, distance with inclination and coalescence time with reference phase. Its condition number at this point is 6.1e+12, and fisher_covariance warns about it. Read that warning as a reason to prefer a sampled posterior, not as an error. fisher_covariance deliberately does not discard small eigenvalues, because a small Fisher eigenvalue is a genuinely poorly-measured direction with a large variance, not numerical noise.

Runnable: examples/usage/07_fisher_matrix.py.

8. The model's own error estimate

The artifact stores enough to report how uncertain the model is at a parameter point, not just a waveform. Ask the MODEL, which is the only thing that knows how its content is laid out:

point = (4.0, 0.3, -0.2)          # mass ratio, spin1z, spin2z

times, error = model.get_waveform_error_estimate(*point)
delta_amplitude, delta_phase = error[(2, 2)]

It returns (samples, {(l, m): (delta_amplitude, delta_phase)}) for every mode the model serves: a 1-sigma amplitude uncertainty in the same units as the mode, and a phase uncertainty in radians, both arrays as long as samples. delta_phase is the envelope divided by the mode amplitude, so it means something only where the mode stands well above its own envelope.

domain chooses what the answer is parameterised by:

Orbital phase ends at the amplitude peak, so the phase domain cannot say anything about merger or ringdown; that is why time is the default. On the time domain the inspiral envelope is evaluated at Phi(t) through the model's own clock inversion, the merger-ringdown arm supplies its own envelope, and the two are joined by the same cosine ramp the waveform's own join uses. Pass your own grid as samples= in whichever unit the domain selects.

At mass_ratio = 4.0, spin1z = 0.3, spin2z = -0.2, over the whole waveform (t from -3819.079 to +124.421 M):

"Near the peak" is the median over t in (-10, +30) M. The same seven modes on the inspiral alone, domain="orbital_phase" over model.phi:

The dominant mode is the best determined of the seven, by a factor of 5 over the next one.

What this is, and what it is not

It is the surrogate's own predictive uncertainty: a 1-sigma statement from the fits about interpolating between the simulations they were trained on. Outside the trained parameter range it is an extrapolation, the call warns, and the answer should be read as a lower bound:

predict_var queried outside the training hull; the envelope is a GPR
extrapolation and should be treated as a lower bound (check sur.in_hull).

It is NOT a measurement of how far the model is from numerical relativity. The merger-ringdown arm is the model's weakest part and is recorded as a 100 to 200 percent problem at parameter-estimation masses, regression dominated; the envelope over that arm reports 0.4 to 3.6 percent. Both statements are true, and they are about different things: a posterior variance describes interpolation uncertainty and cannot see a bias in the fit itself. Do not read a small envelope as an accurate arm.

Two channels are also missing from the per-mode number, and both are reachable but not folded in. The first is the clock: tau has its own 1-sigma, a timing error is a phase error, and on (2,2) the timing-induced phase uncertainty is about 28 times the mode-content one reported above. See get_dynamics_error_estimate below and docs/wiki/error-aware-likelihood.md limit 2. The second is coefficient-to-coefficient correlation, which the diagonal propagation drops.

One more thing to know about the inspiral half. Every mode of the action-angle model is served as a weighted combination of two routes, the direct fit of the mode content and the fit of its ratio to the post-Newtonian expression, with one scalar weight per route per mode. The reported uncertainty is that of the route carrying the weight, not of the combination that is served: it is the direct route's envelope over the direct route's amplitude. Measured at mass_ratio = 4, spins (0.3, -0.2), the direct route carries 0.9898 or more of the weight on six modes and 0.849 on (2,2), so on six modes the band describes what you get to one percent and on the dominant mode it describes 85 percent of it. It is not propagated through the combination because the two routes are not independent, both being fitted from the same 260 simulations, so the textbook inverse-variance result understates, and it would narrow a band already measured a factor 2.6 optimistic in sample. It is therefore documented rather than changed, with improvements riding a later release; docs/wiki/model-error-estimate.md carries the weights for every mode and the fix that is planned. The adiabatic-angle model has no route split, so the statement is about the action-angle model only.

The dynamics uncertainties

The clock and the other fitted dynamics fields have their own stored variances, reached through a separate call so nothing is computed unless you ask for it (4.9 ms against the 109 ms a waveform costs):

orbital_phase, sigma = model.get_dynamics_error_estimate(*point)
sigma["tau"]      # 1-sigma of the clock, in M, over orbital_phase

At the same point: tau 2.3796e-01 M median, x 1.6956e-05, J 8.5713e-05, E 1.5347e-06, k 7.1095e-03 (one number, constant, returned broadcast over the grid). omega and eps are not reported: they are derivatives of the fitted clock rather than fitted fields.

The content-group functions

Three lower-level functions take a content group rather than the model, and are what the model-level calls are built from:

from nrhjsurrogate.waveform_error_estimate import (predict_var,
                                                    error_envelope,
                                                    mode_errors)

predict_var(model.sur_h, *point)               # name -> coefficient variances
error_envelope(model.sur_h, point, model.phi)  # name -> 1-sigma envelope
mode_errors(model.sur_h, point, model.phi)     # (l,m) -> (delta_A, delta_phi)

A group only knows its own content, and the AA arm spreads one mode over three groups: the odd-m direct coefficients live in model.sur_h keyed by (l, m), the even-m ones in model.sur_e keyed by ("D", l, m), and the PN-ratio route in model.sur_o and model.sur_e keyed by ("R", l, m). So mode_errors(model.sur_h, ...) returns the four odd-m modes and mode_errors(model.sur_e, ...) returns nothing at all, because that group has no bare (l, m) keys to return. That is why the per-mode question is asked of the model. The AdA arm keeps all seven modes in one group, so mode_errors on its inspiral surrogate does return seven.

A merger-ringdown arm needs merger_ringdown_error_envelope from the same module, never error_envelope directly: under the SWSpH representation the regressed coefficients are the detrended spheroidal residual, so the plain envelope returns the uncertainty of a different quantity.

As a Fisher noise contribution

The same estimate can be folded into a Fisher matrix as an effective noise contribution. waveform_error_psd carries the per-mode envelope through the polarization sum and reads the result as a power spectral density S_w, and fisher_matrix(..., error_psd=...) then weights by 1 / (S_n + S_w), so parameter widths grow wherever the model is less sure:

from nrhjsurrogate.driver.fisher import waveform_error_psd

error_psd = waveform_error_psd(strain_model, inner_product, grid)
fisher = fisher_matrix(strain_model, params, grid, inner_product,
                       error_psd=error_psd)

It joins the two arms with the same two functions the model-level call uses, so the two routes cannot drift apart. How much it matters depends on how loud the signal is. At 60 solar masses, q = 2 and 400 Mpc the 1-sigma widths are unchanged to four figures. Move the same source to 40 Mpc and every width grows by a factor 1.001. This is the leading noise-inflation term only, under a stationary and diagonal approximation to the model-error covariance; the derivation and its limits are in docs/wiki/error-aware-likelihood.md.

Runnable: examples/usage/08_model_error.py.

9. Backends

model = NRHJSurAA.load_h5("NRHJSur3dq8_AA_v3", backend="numpy")   # default
model = NRHJSurAA.load_h5("NRHJSur3dq8_AA_v3", backend="kokkos")  # compiled

"numpy" is the reference path and needs nothing beyond the declared dependencies. "kokkos" routes the same evaluation through a compiled extension module.

Pointing at the compiled module. The package looks for an nrhjkokkos*.so under cpp/build* beside the package, newest first. Setting $NRHJ_KOKKOS_DIR to a directory replaces that search entirely: that one directory is used and nothing else is consulted.

export NRHJ_KOKKOS_DIR=/path/to/the/build/directory

nrhjsurrogate.execution.compiled_module_loader.load_module().__file__ tells you which module was actually loaded.

Agreement. On the same parameter point the two backends return the same time grid and modes agreeing to between 4.9e-09 and 2.2e-08 relative, the worst being the (2, 2) mode.

Speed. Measured on one machine at OMP_NUM_THREADS=8, quoted in milliseconds per waveform:

These are one machine's numbers, not a specification. Measure on yours.

Threads. Importing the package sets OMP_PROC_BIND=close and OMP_PLACES=cores, but only if you have not already set either of them, and it never sets OMP_NUM_THREADS, which is yours to choose. The OpenMP runtime reads the binding variables once when it starts, so setting them after the import is silently ignored. Set them in the shell. nrhjsurrogate.thread_binding() asks the runtime what it actually did, rather than echoing what was requested. To opt out of the defaults entirely, set NRHJ_NO_THREAD_DEFAULTS to anything non-empty.

Runnable: examples/usage/09_backends.py.

10. The adiabatic-angle model

from nrhjsurrogate import RELEASE_ADA
from nrhjsurrogate.driver.api import NRHJSurAdA

model = NRHJSurAdA.load(RELEASE_ADA, backend="numpy")

This class carries the eight-parameter physical surface: detector strain, its derivatives with respect to all eight parameters, and its time derivative. It is the class the Fisher matrix of section 7 goes through.

Two differences from the action-angle class matter immediately.

It is not callable. model(q, chi1z, chi2z) raises TypeError: 'NRHJSurAdA' object is not callable. Use .modes(...).

.modes(...) returns a three-item tuple unless you ask otherwise.

names, time, modes = model.modes(q, chi1z, chi2z)   # h: (batch, modes, samples)
waveform = model.modes(q, chi1z, chi2z, as_dict=True)  # {"t": ..., (l,m): ...}

Forgetting as_dict=True and then using the result as if it were a dictionary or a number gives errors that do not mention the real cause, such as float() argument must be a string or a real number, not 'tuple', or tuple indices must be integers or slices, not str. If you see one of those, this is why.

The physical surface, with params rows in (mtot, q, chi1z, chi2z, dist, t_c, phi_ref, incl) order and one shared time grid in seconds:

h_plus, h_cross = model.strain(params, grid)      # (B, N) each
derivatives = model.strain_grad(params, grid)     # (B, 8, N)
dh_dt = model.time_derivative(params, grid)       # (B, N)
_, _, mode_derivatives = model.modes_grad(q, chi1z, chi2z)  # (B, 3, K, N)

model.hull_mask(q, chi1z, chi2z) is the batched in-range check.

The default backend for this class is "kokkos", unlike the action-angle class. strain_grad and time_derivative always go through the compiled path.

The two artifacts are not interchangeable between the two classes. Handing an action-angle file to NRHJSurAdA.load raises unsupported format_version.

Runnable: examples/usage/06_ada_model.py.

11. Dynamics and remnant

These are separate products, so they are separate methods rather than flags on a waveform call.

model.get_orbital_dynamics(mass_ratio, spin1z, spin2z) returns the orbital quantities the waveform is built on, at one parameter point, as arrays over the model's own phase grid of 4000 points: "phi", "tau" (the clock), "omega" (orbital frequency), "x" (the post-Newtonian parameter), "eps" (a measure of how far from adiabatic the evolution is) and the scalar "k".

model.get_remnant_properties_geometric(mass_ratio, spin1z, spin2z) returns the final black hole's mass in units of the total mass and its dimensionless spin. model.get_remnant_properties(mass1, mass2, spin1z, spin2z) returns the mass in solar masses instead. At mass_ratio = 4, spins (0.3, -0.2) these are M_f/M = 0.973603 and chi_f = 0.624352, and for mass1 = 50, mass2 = 12.5 that is M_f = 60.8502 solar masses.

The older names model.dynamics, model.remnant and model.remnant_physical are retained and forward to these.

12. Things that catch people out

• The merger-ringdown artifact is required. Deposit NRHJSur3dq8_MR_v3.h5 next to the model file. Without it the load fails outright.
• The polarizations need inclination. get_td_waveform requires it and a distance; the modes methods do not take an inclination at all.
• distance decides whether you get strain or r h/M on the modes methods; the geometric methods have no distance argument.
• with_derivatives is the only argument whose value changes the arity, and only on get_td_waveform_modes_geometric.
• "dynamics" and "hybrid_report" are always present in a result dictionary and are None when they do not apply. Code that iterated a result and skipped only "t" needs to test isinstance(key, tuple) instead.
• An adiabatic-angle model is not callable, and its .modes needs as_dict=True to give you a dictionary.
• in_hull is the authority on the parameter range, not a remembered rule of thumb. q = 8.0 exactly reports False, because the range stops a hair short of it.
• f_low=0 means the same as f_low=None, the full span, following the same convention as the LAL mode generators.
• The merger-ringdown argument is merger_ringdown_path on every loader. It used to be mr_path on load and mr on load_h5, two names for one argument on one class. Both still work and warn.
• A time grid outside the waveform is an error, not zero padding, in the Fisher path.
• predict_var outside the trained range is an extrapolation and its answer is a lower bound.
• safe=False does not mean "less accurate". It selects the fast evaluation path instead of the one that reproduces bit for bit. Measured on the released artifact over 9 parameter points and 7 modes, the largest departure is 4.53e-10 relative on the waveform and 1.43e-09 on its derivatives, which is 1.4e-6 of the implementation divergence already present: accuracy 1.00, unchanged. What you give up is bit identity, so a bitwise gate will go red. The keyword becomes bitwise_reproducible in a later release, with safe= still accepted.
• safe is process-wide, not per call or per thread. A thread that needs the canonical path, running while another thread holds a fast scope open, took the fast path on 200 of 200 evaluations, and two interleaved scopes can leave the flag stuck fast with no scope open at all. Use one flavour per process, or keep threads out of it. A threaded likelihood is exactly the shape that hits this.
• Set OMP_PROC_BIND before importing anything, or the setting is silently ignored.