qscat.core¶

The model-independent electron–diatomic scattering engine. ScatteringProblem is the recommended entry point: it bundles the grid, model, and vibrational basis once and exposes the supported observables as methods (the class docstring lists what is deliberately functional-only). The functional solvers below are the low-level layer those methods call.

This module never imports qscat.model at runtime — it depends only on the ResonanceModel protocol, so a new molecule needs no change here.

The problem object¶

class qscat.core.ScatteringProblem(grid, model, n_vib, v_init=0)[source]¶

A fully-specified electron-diatomic scattering problem.

Bundles the grid, model, and vibrational basis shared by every observable. grid is the electronic x nuclear TensorGrid; model is any qscat.model.ResonanceModel; n_vib is the number of vibrational states to solve; v_init is the initial vibrational level (default 0). The basis is solved once on construction and exposed as .eps / .chi / .basis.

Parameters:

grid (TensorGrid)
model (ResonanceModel)
n_vib (int)
v_init (int)

property eps: NDArray[float64]¶: Vibrational energies (Ha), ascending.

property chi: NDArray[complex128]¶: Vibrational eigenvectors, one row per level.

ve_cross_section(vprimes, E, *, ordering='COLAMD', lam_scale=1.0, return_wavefunction=False)[source]¶

Vibrational-excitation cross section; same parameters, defaults and return convention as qscat.core.ve_cross_section (which this delegates to with the bundled grid/model/basis). Returns the plain sigma array, or (sigma, psi) when return_wavefunction=True; the wavefunction is None per energy below threshold.

Parameters:

vprimes (list[int])
E (float | ArrayLike)
ordering (Literal['NATURAL', 'MMD_ATA', 'MMD_AT_PLUS_A', 'COLAMD'])
lam_scale (float)
return_wavefunction (bool)

Return type:

NDArray[float64] | tuple[NDArray[float64], NDArray[complex128] | None | list[NDArray[complex128] | None]]

da_cross_section(E, *, n_channels=1, ordering='COLAMD', return_wavefunction=False)[source]¶

Dissociative-attachment cross section; see qscat.core.da_cross_section. Returns the plain sigma array, or (sigma, psi) when return_wavefunction=True; the wavefunction is None per energy below threshold.

Parameters:

E (float | ArrayLike)
n_channels (int)
ordering (Literal['NATURAL', 'MMD_ATA', 'MMD_AT_PLUS_A', 'COLAMD'])
return_wavefunction (bool)

Return type:

NDArray[float64] | tuple[NDArray[float64], NDArray[complex128] | None | list[NDArray[complex128] | None]]

dr_cross_section(E, *, n_channels=3, ordering='COLAMD', return_wavefunction=False, return_amplitude=False)[source]¶

Dissociative-recombination cross section (ionic target); see qscat.core.dr_solve for the physics. Returns the plain sigma array, (sigma, psi) under return_wavefunction=True, (sigma, amplitude) under return_amplitude=True, or (sigma, psi, amplitude) under both; the wavefunction is None per energy below threshold. Callers that want the whole DrResult in one object call qscat.core.dr_solve directly.

Parameters:

E (float | ArrayLike)
n_channels (int)
ordering (Literal['NATURAL', 'MMD_ATA', 'MMD_AT_PLUS_A', 'COLAMD'])
return_wavefunction (bool)
return_amplitude (bool)

Return type:

td_ve_cross_section(vprimes, E, *, dt, n_steps, wp_in, wp_out=None, order=3, subtract_free_reference=True, method='tw', position=None, surface=None)[source]¶

Time-dependent VE cross section; see qscat.core.td_ve_cross_section (same method/wp_out/position/surface contract).

Parameters:

vprimes (list[int])
E (float | ArrayLike)
dt (float)
n_steps (int)
wp_in (dict[str, float])
wp_out (dict[str, float] | None)
order (int)
subtract_free_reference (bool)
method (Literal['tw', 'delta', 'flow'])
position (int | None)
surface (int | None)

Return type:

NDArray[float64]

td_ve_cross_sections_all(vprimes, E, *, dt, n_steps, wp_in, wp_out, position, surface, order=3, subtract_free_reference=True)[source]¶

All three TD-VE extractors from ONE propagation; see qscat.core.td_ve_cross_sections_all.

Parameters:

vprimes (list[int])
E (float | ArrayLike)
dt (float)
n_steps (int)
wp_in (dict[str, float])
wp_out (dict[str, float])
position (int)
surface (int)
order (int)
subtract_free_reference (bool)

Return type:

dict[str, NDArray[float64]]

td_da_cross_section(E, *, dt, n_steps, wp_in, method='flow', surface=None, position=None, wp_out=None, n_channels=1, order=3)[source]¶

Time-dependent DA cross section; see qscat.core.td_da_cross_section.

Parameters:

E (float | ArrayLike)
dt (float)
n_steps (int)
wp_in (dict[str, float])
method (Literal['tw', 'delta', 'flow'])
surface (int | None)
position (int | None)
wp_out (dict[str, float] | None)
n_channels (int)
order (int)

Return type:

NDArray[float64]

td_da_cross_sections_all(E, *, dt, n_steps, wp_in, surface, position, wp_out, n_channels=1, order=3)[source]¶

All three TD-DA extractors from ONE propagation; see qscat.core.td_da_cross_sections_all.

Parameters:

E (float | ArrayLike)
dt (float)
n_steps (int)
wp_in (dict[str, float])
surface (int)
position (int)
wp_out (dict[str, float])
n_channels (int)
order (int)

Return type:

dict[str, NDArray[float64]]

lcp_da_cross_section(E, *, Vd, Gamma, ordering='COLAMD', return_wavefunction=False)[source]¶

LCP dissociative-attachment cross section on this problem’s NUCLEAR grid; see qscat.core.lcp_da_cross_section. The curve (Vd, Gamma) is per-call (compute it with resonance_levels(…, return_curve=True) – see that docstring for why not local_complex_potential directly); mu/eps/chi/v_init come from the bundle, which is what pays down the functional signature’s documented argument-order exception (docs/adr/0007). The LCP magnitude needs the FINE per-molecule nuclear deck – construct the problem on it for physical numbers. Returns the plain sigma array, or (sigma, psi) when return_wavefunction=True; the wavefunction is None per energy below threshold.

Parameters:

E (float | ArrayLike)
Vd (NDArray[complex128])
Gamma (NDArray[float64])
ordering (Literal['NATURAL', 'MMD_ATA', 'MMD_AT_PLUS_A', 'COLAMD'])
return_wavefunction (bool)

Return type:

NDArray[float64] | tuple[NDArray[float64], NDArray[complex128] | None | list[NDArray[complex128] | None]]

lcp_ve_cross_section(vprimes, E, *, Vd, Gamma, ordering='COLAMD', return_wavefunction=False)[source]¶

LCP vibrational-excitation cross section on this problem’s NUCLEAR grid; see qscat.core.lcp_ve_cross_section. The curve (Vd, Gamma) is per-call (compute it with resonance_levels(…, return_curve=True) – see that docstring for why not local_complex_potential directly); mu/eps/chi/v_init come from the bundle, which is what pays down the functional signature’s documented argument-order exception (docs/adr/0007). The LCP magnitude needs the FINE per-molecule nuclear deck – construct the problem on it for physical numbers. Returns the plain sigma array, or (sigma, psi) when return_wavefunction=True; the wavefunction is None per energy below threshold.

Parameters:

vprimes (list[int])
E (float | ArrayLike)
Vd (NDArray[complex128])
Gamma (NDArray[float64])
ordering (Literal['NATURAL', 'MMD_ATA', 'MMD_AT_PLUS_A', 'COLAMD'])
return_wavefunction (bool)

Return type:

NDArray[float64] | tuple[NDArray[float64], NDArray[complex128] | None | list[NDArray[complex128] | None]]

resonance_levels(nuclear_grid_b, elec_grid_b, *, re_half_width=0.05, im_half_width=0.05, resid_tol=0.001, window=None, n_levels=None, rel_tol=0.0001, atol=1e-08, golden_rule=True, return_curve=False)[source]¶

Born-Oppenheimer quasi-bound levels of this problem’s anion; see qscat.core.resonance_levels. This problem’s own electronic/nuclear grids are the _a partners; nuclear_grid_b/elec_grid_b are the angle-moved partners (share every real node, differ only in the ECS tail angle – qscat.core.grids.ecs_angle_family builds a valid family). return_curve=True also returns the (Vd, Gamma) curve the levels were computed in – the input lcp_da_cross_section needs. Returns the plain ResonanceLevels, or (levels, Vd, Gamma) when return_curve=True.

Parameters:

nuclear_grid_b (FemDvrEcsGrid)
elec_grid_b (FemDvrEcsGrid)
re_half_width (float)
im_half_width (float)
resid_tol (float)
window (_Window | None)
n_levels (int | None)
rel_tol (float)
atol (float)
golden_rule (bool)
return_curve (bool)

Return type:

ResonanceLevels | tuple[ResonanceLevels, npt.NDArray[np.complex128], npt.NDArray[np.float64]]

exact_resonance_states(grid_electronic, grid_nuclear, *, shifts, window, k=8, rel_tol=0.0001, atol=1e-08)[source]¶

Exact 2-D resonance states by two-angle ECS stability; see qscat.core.exact_resonance_states. This problem’s grid is the base; grid_electronic/grid_nuclear are the one-angle-moved partner TensorGrids (ecs_angle_family builds all three consistently). Seeds (shifts) are passed in – typically resonance_levels’s output – so the exact solver never depends on the approximation it measures.

Parameters:

grid_electronic (TensorGrid)
grid_nuclear (TensorGrid)
shifts (npt.ArrayLike)
window (_Window)
k (int)
rel_tol (float)
atol (float)

Return type:

ExactResonanceStates

nrm_ve_cross_section(phi_d, vprimes, E, *, ingredients=None, n_states=None, include_background=True)[source]¶

NRM vibrational-excitation cross section; see qscat.core.nrm.nrm_ve_cross_section (this problem’s nuclear and electronic grids fill its leading NUCLEAR-grid-first pair).

Parameters:

phi_d (DiscreteState)
vprimes (list[int])
E (float | npt.ArrayLike)
ingredients (NrmIngredients | None)
n_states (int | None)
include_background (bool)

Return type:

npt.NDArray[np.float64]

nrm_da_cross_section(phi_d, E, *, ingredients=None, n_states=None)[source]¶

NRM dissociative-attachment cross section; see qscat.core.nrm.nrm_da_cross_section.

Parameters:

phi_d (DiscreteState)
E (float | npt.ArrayLike)
ingredients (NrmIngredients | None)
n_states (int | None)

Return type:

npt.NDArray[np.float64]

Cross sections (time-independent)¶

qscat.core.ve_cross_section(tgrid: TensorGrid, model: ResonanceModel, eps: NDArray[float64], chi: NDArray[complex128], v_init: int, vprimes: list[int], E: float | ArrayLike, *, ordering: Literal['NATURAL', 'MMD_ATA', 'MMD_AT_PLUS_A', 'COLAMD'] = 'COLAMD', lam_scale: float = 1.0, return_wavefunction: Literal[False] = False) → NDArray[float64][source]¶

qscat.core.ve_cross_section(tgrid: TensorGrid, model: ResonanceModel, eps: NDArray[float64], chi: NDArray[complex128], v_init: int, vprimes: list[int], E: float | ArrayLike, *, ordering: Literal['NATURAL', 'MMD_ATA', 'MMD_AT_PLUS_A', 'COLAMD'] = 'COLAMD', lam_scale: float = 1.0, return_wavefunction: bool = False) → NDArray[float64] | tuple[NDArray[float64], NDArray[complex128] | None | list[NDArray[complex128] | None]]

sigma_{v_init->v’}(E) in bohr^2, exact 2-D driven-equation solution.

Provisional API (docs/adr/0004-public-api-stability-policy.md): this wide functional signature is the layer the context-object refactor targets and may change in a minor release; ScatteringProblem.ve_cross_section is the stable route.

model supplies the Hamiltonian (model.hamiltonian(tgrid)), the interaction diagonal (model.interaction_diag(tgrid)) and the fixed partial wave (model.ell) – the entire molecule-specific input.

model.charge is forwarded to channel_vector, so an IONIC target’s entrance/exit channels are the energy-normalized Coulomb functions (coulomb_f_en) rather than the free riccati_bessel_en; charge == 0 (every neutral model) takes the identical pre-existing free-function branch. This is what lets dissociation.dr_cross_section reuse this solver’s driven sweep for H2+ instead of re-inlining it.

E may be scalar or an array; scalar returns shape (len(vprimes),), array returns (len(E), len(vprimes)). One sparse LU per energy is reused across all vprimes.

lam_scale scales V_int ONLY, for the free-particle and first-Born validation limits. It is a test lever, never a physics knob.

If return_wavefunction, also returns psi_plus (or None when E <= 0, since no driven-equation solve happens below threshold): one array for scalar E, one list entry per energy for array E.

Parameters:

tgrid (TensorGrid)
model (ResonanceModel)
eps (npt.NDArray[np.float64])
chi (npt.NDArray[np.complex128])
v_init (int)
vprimes (list[int])
E (float | npt.ArrayLike)
ordering (Ordering)
lam_scale (float)
return_wavefunction (bool)

Return type:

_Sigma | tuple[_Sigma, _PsiOut]

qscat.core.da_cross_section(tgrid: TensorGrid, model: ResonanceModel, eps: NDArray[float64], chi: NDArray[complex128], v_init: int, E: float | ArrayLike, *, n_channels: int = 1, ordering: Literal['NATURAL', 'MMD_ATA', 'MMD_AT_PLUS_A', 'COLAMD'] = 'COLAMD', return_wavefunction: Literal[False] = False) → NDArray[float64][source]¶

qscat.core.da_cross_section(tgrid: TensorGrid, model: ResonanceModel, eps: NDArray[float64], chi: NDArray[complex128], v_init: int, E: float | ArrayLike, *, n_channels: int = 1, ordering: Literal['NATURAL', 'MMD_ATA', 'MMD_AT_PLUS_A', 'COLAMD'] = 'COLAMD', return_wavefunction: bool = False) → NDArray[float64] | tuple[NDArray[float64], NDArray[complex128] | None | list[NDArray[complex128] | None]]

sigma_DA(E) in bohr^2, exact 2-D driven-equation DA cross section.

Reuses ve_cross_section(…, return_wavefunction=True) for Psi+ (one SparseLU.refactor sweep across E), then reads the outgoing dissociation flux out of it, one n_channels anion channel at a time. E may be scalar (returns (n_channels,)) or an array (returns (len(E), n_channels)). sigma = 0 for a closed channel (E <= 0 or E_DR = E_tot - eps_e <= 0).

The extraction projects Psi+ onto the anion electronic state to get the channel’s nuclear wave, then takes its VALUE at the outermost real nuclear node X (the DVR coefficient divided by sqrt(w_b), never the coefficient itself):

psi_n(R) = <phi_e^(n) | Psi+>_r      [c-product over the electronic axis]
S_n      = sqrt(K_n / 2 pi mu) psi_n(X)
sigma_n  = 4 pi^3 |S_n|^2 / 2E

– algebraically the same amplitude qscat.core.lcp.lcp_da_cross_section takes off its 1-D nuclear wave, so the exact solver and the LCP are compared through identical arithmetic.

Note

Why not the post-form volume T-matrix. The textbook alternative, T = <phi_e^(n) F_K | V_DR | Psi+> with the rearrangement interaction V_DR = V_int(r,R) + v0(R) - V_int(r,R_inf) (still available as v_dr_diag, and what dr_cross_section uses), is formally exact and agrees with the flux above to 5e-4 on F2. It is nonetheless the WRONG tool whenever sigma_DA is small: V_DR does not decay in r, so the integrand’s magnitude is set by the interaction region while the answer is set by how completely that region cancels. On NO – sigma_DA ~ 1e-9 bohr^2 against an integrand summing to ~2.6 – the required cancellation is ~1e6-fold, and whatever has not decayed at either edge of the integration region survives instead of the physics. Both edges did: measured on the shipped decks, the T-matrix answer moved by four orders when the electronic real region went 16 -> 48 bohr and by a further ~800x when the nuclear R_inf went 9.0 -> 15.0 bohr (where NO’s Morse v0 finally reaches 1e-9 Ha), converging only there onto the flux value. The flux extraction needs no cancellation and is invariant under both to 4 digits. See docs/physics/diatomic-ve-cross-sections.md and validation/diatomic/test_no_da_thesis.py.

If return_wavefunction, also returns the driven Psi+ (the SAME solution the flux is read from; None for E <= 0): one array for scalar E, one list entry per energy for an array E – same convention as ve_cross_section.

Parameters:

tgrid (TensorGrid)
model (ResonanceModel)
eps (npt.NDArray[np.float64])
chi (npt.NDArray[np.complex128])
v_init (int)
E (float | npt.ArrayLike)
n_channels (int)
ordering (Ordering)
return_wavefunction (bool)

Return type:

_Sigma | tuple[_Sigma, _PsiOut]

qscat.core.dr_cross_section(tgrid, model, eps, chi, v_init, E, *, n_channels=3, ordering='COLAMD')[source]¶

sigma_DR(E) in bohr^2, exact 2-D driven-equation dissociative- recombination cross section for a CHARGED target (e.g. H2+, charge=-1). See dr_solve for the physics; this is a thin sigma-only wrapper around it – a caller that also wants the wavefunction or T-matrix amplitude calls dr_solve directly and reads the DrResult fields.

Parameters:

tgrid (TensorGrid)
model (ResonanceModel)
eps (npt.NDArray[np.float64])
chi (npt.NDArray[np.complex128])
v_init (int)
E (float | npt.ArrayLike)
n_channels (int)
ordering (Ordering)

Return type:

_Sigma

qscat.core.dr_solve(tgrid, model, eps, chi, v_init, E, *, n_channels=3, ordering='COLAMD', store_wavefunction=False, store_amplitude=False)[source]¶

Exact 2-D driven-equation dissociative-recombination solve for a CHARGED target (e.g. H2+, charge=-1); returns a DrResult.

da_cross_section GENERALIZED two ways: (1) the incident channel is Coulomb rather than free – channel_vector(…, charge=model.charge); (2) the exit channel is a LOOP over n_channels Rydberg electronic states of the neutral, anion_electronic_states(…, n_states=n_channels) (the SAME bound-electronic-state solver DA uses for the anion state – the Rydberg series is likewise bound, below the -1/r Coulomb continuum). V_DR (the rearrangement interaction, NOT V_int) is unchanged from DA.

The driven Lippmann-Schwinger solve for Psi+ reuses ve_cross_section(…, return_wavefunction=True) – the same analyze-once / SparseLU.refactor-per-energy sweep, with model.charge forwarded to channel_vector so the incident channel is Coulomb – exactly as da_cross_section does. Only the exit-channel read differs: DR projects the post-form volume T-matrix against V_DR (below), DA reads a boundary flux.

E may be scalar (returns (n_channels,) fields) or an array (returns (len(E), n_channels) fields). sigma_n = 0 for a closed channel (E <= 0 or E_DR = E_tot - E_ryd(n) <= 0, E_ryd(n) = eps_e[n]).

DrResult.psi is populated (the driven Psi+, None for E <= 0) only when store_wavefunction is set: one array for scalar E, one list entry per energy for an array E – same convention as ve_cross_section/da_cross_section.

DrResult.amplitude is populated only when store_amplitude is set: the complex transition amplitude t the T-matrix sum is built from, shaped exactly like sigma ((n_channels,) for scalar E, (len(E), n_channels) for an array E, zero for a closed channel): sigma = 4*pi**3 * abs(t)**2 / (2*E). This is the amplitude the solver already forms, NOT a literal unitary S-matrix element – the thesis’s S_DR differs from it by the standard S = -2*pi*i*T factor and its own normalization. That factor is a fixed rotation and rescale: it changes neither the zeros nor the shape of t’s real/imaginary crossings, so returning t as computed (rather than guessing a normalization to synthesize an “S”) is what downstream resonance-pole fitting needs.

Parameters:

tgrid (TensorGrid)
model (ResonanceModel)
eps (npt.NDArray[np.float64])
chi (npt.NDArray[np.complex128])
v_init (int)
E (float | npt.ArrayLike)
n_channels (int)
ordering (Ordering)
store_wavefunction (bool)
store_amplitude (bool)

Return type:

DrResult

class qscat.core.DrResult(sigma, psi, amplitude)[source]¶

Result of the exact 2-D dissociative-recombination solve (dr_solve).

sigma: sigma_DR(E) in bohr^2 – (n_channels,) for scalar E, (len(E), n_channels) for array E.
psi: the driven Psi+ per energy (None when not stored, and None per energy below threshold) – one array for scalar E, one list entry per energy for array E.
amplitude: the complex T-matrix amplitude, shaped like sigma (None when not stored). See dr_solve for the S-vs-T normalization note.

Parameters:

sigma (NDArray[float64])
psi (NDArray[complex128] | None | list[NDArray[complex128] | None])
amplitude (NDArray[complex128] | None)

Cross sections (time-dependent)¶

qscat.core.td_ve_cross_section(tgrid, model, eps, chi, v_init, vprimes, E, *, dt, n_steps, wp_in, wp_out=None, order=3, subtract_free_reference=True, method='tw', position=None, surface=None)[source]¶

sigma_{v_init->v’}(E) (bohr^2), Pade propagation + an energy-extraction method’s transform: method=”tw” (default, Tannor-Weeks – a propagated Gaussian test packet), “delta” (a fixed-point line projection, eMoScat’s DiracTestFunction2d), or “flow” (a fixed-surface Wronskian flux, eMoScat’s FluxTestFunction2d) – see td_extractors.py for the three Extractor implementations and their formulas.

E (collision energy, Hartree) may be scalar or array-like; scalar E returns shape (len(vprimes),), array E returns (len(E), len(vprimes)) – matching qscat.core.driven.ve_cross_section, the exact differential oracle this converges to as dt -> 0 and n_steps -> infinity (see module docstring).

wp_in = {“r0”: …, “p0”: …, “sigma”: …} are the SAME incident Gaussian parameters used to build Psi(0) (via initial_state) and eta_incident; wp_out = {“r0_out”: …, “p0_out”: …, “sigma_out”: …} are the outgoing test function’s parameters – required by method=”tw” (omitting them raises ValueError, mirroring td_da_cross_section’s contract) and unused by “delta”/”flow”, which need no propagated test packet. The propagation (the expensive part) happens ONCE regardless of how many energies E are requested, since the recorded per-step series does not depend on E.

position (required, keyword-only, for method=”delta”) and surface (required, keyword-only, for method=”flow”) are the fixed electronic DVR index the Dirac/Flux extractor analyzes at – an index in the real (unscaled) region past the interaction, mirroring wp_out’s asymptotic standoff (see Dirac/Flux’s own docstrings for the validity requirement). Omitting the one the selected method needs raises ValueError.

subtract_free_reference (default True): when the diagonal/elastic channel is requested (v_init in vprimes), a SECOND V_int=0 propagation is run to supply the free-particle reference S_free(E) that the elastic channel subtracts (instead of a literal 1) – required for a correct elastic cross section, see sigma_one_energy. It doubles the propagation cost and is a no-op (skipped) when the elastic channel is not requested; set False to force the old literal-1 behavior. The inelastic channels are identical either way.

All three methods share one code path: the selected extractor is built (twice, identically, when the free-reference run applies), driven by propagate(…, out_channels=[], extractors=[…]), and asked for sigma(E, free=…). method=”tw” reproduces this function’s pre-refactor implementation to machine precision (see libs/qscat/tests/test_td_extractors.py’s golden regression tests); any other method value raises ValueError.

Parameters:

tgrid (TensorGrid)
model (ResonanceModel)
eps (npt.NDArray[np.float64])
chi (npt.NDArray[np.complex128])
v_init (int)
vprimes (list[int])
E (float | npt.ArrayLike)
dt (float)
n_steps (int)
wp_in (_WpIn)
wp_out (_WpOut | None)
order (int)
subtract_free_reference (bool)
method (Method)
position (int | None)
surface (int | None)

Return type:

npt.NDArray[np.float64]

qscat.core.td_ve_cross_sections_all(tgrid, model, eps, chi, v_init, vprimes, E, *, dt, n_steps, wp_in, wp_out, position, surface, order=3, subtract_free_reference=True)[source]¶

{“tw”: sigma, “delta”: sigma, “flow”: sigma} (each bohr^2, same shape convention as td_ve_cross_section) from ONE shared propagation.

This is the HONEST three-way comparison: TannorWeeks, Dirac, and Flux are all built up front and driven by a SINGLE propagate(…, extractors=[tw, dirac, flux]) call (and a single companion V_int=0 propagation, when subtract_free_reference applies) – identical dynamics psi(t_n) feed all three transforms, so any spread between the returned cross sections reflects a genuine difference between the energy-extraction methods (or, at an under-converged grid/propagation length, a shared discretization/truncation residual all three inherit together – see docs/physics/td-extractors.md), never a difference in what was propagated.

position/surface are the Dirac/Flux fixed electronic DVR indices (see td_ve_cross_section’s docstring); both are required here (unlike td_ve_cross_section, which only needs whichever one its method selects).

Parameters:

tgrid (TensorGrid)
model (ResonanceModel)
eps (npt.NDArray[np.float64])
chi (npt.NDArray[np.complex128])
v_init (int)
vprimes (list[int])
E (float | npt.ArrayLike)
dt (float)
n_steps (int)
wp_in (_WpIn)
wp_out (_WpOut)
position (int)
surface (int)
order (int)
subtract_free_reference (bool)

Return type:

dict[str, npt.NDArray[np.float64]]

qscat.core.td_da_cross_section(tgrid, model, eps, chi, v_init, E, *, dt, n_steps, wp_in, method='flow', surface=None, position=None, wp_out=None, n_channels=1, order=3)[source]¶

sigma_DA(E) (bohr^2), Pade propagation + a NUCLEAR-axis energy- extraction method’s transform – the dissociative-attachment (DA) sibling of td_ve_cross_section, built on the SAME propagation engine but the axis=”nuclear” Flux/Dirac/TannorWeeks extractors (td_extractors.py) instead of the electronic-axis ones.

method=”flow” (default, eMoScat’s FluxTestFunction2d axis_==’y’ branch – the natural DA extractor: a fixed-surface Wronskian flux needs no propagated outgoing test packet), “delta” (a fixed nuclear- node point projection), or “tw” (a propagated nuclear outgoing Gaussian test packet) – see td_extractors.py’s module docstring for the three nuclear-axis transforms and the _C_DA = pi reconciliation with the TI dissociation.da_cross_section oracle’s 4*pi^3 T-matrix convention.

E (collision energy, Hartree) may be scalar or array-like; scalar E returns shape (n_channels,), array E returns (len(E), n_channels) – matching dissociation.da_cross_section’s per-anion-channel contract (DA has no v’ vibrational index; n_channels selects how many anion electronic bound states, dissociation. anion_electronic_states at R_inf = tgrid.grids[1].R0, are tracked as exit channels).

wp_in = {“r0”: …, “p0”: …, “sigma”: …} is the SAME incident Gaussian parameters used to build Psi(0) (initial_state) and eta_incident – identical to td_ve_cross_section’s wp_in. surface (required for method=”flow”), position (required for method=”delta”), and wp_out (required for method=”tw”, now the NUCLEAR outgoing test packet’s r0_out/p0_out/sigma_out in R) are the method-specific extractor parameters; omitting the one method needs raises ValueError. Unlike td_ve_cross_section, there is NO wp_out/free-reference for “flow”/”delta” (no propagated test packet needed there) and NO elastic free-reference subtraction for any method – DA is a pure rearrangement channel with no v’==v_init diagonal to subtract a reference from (see td_extractors.py’s Flux/Dirac/TannorWeeks(axis=”nuclear”)’s sigma docstrings, which all raise on a non-None free).

The propagation (the expensive part) happens ONCE regardless of how many energies E are requested, exactly as td_ve_cross_section.

Parameters:

tgrid (TensorGrid)
model (ResonanceModel)
eps (npt.NDArray[np.float64])
chi (npt.NDArray[np.complex128])
v_init (int)
E (float | npt.ArrayLike)
dt (float)
n_steps (int)
wp_in (_WpIn)
method (Method)
surface (int | None)
position (int | None)
wp_out (_WpOut | None)
n_channels (int)
order (int)

Return type:

npt.NDArray[np.float64]

qscat.core.td_da_cross_sections_all(tgrid, model, eps, chi, v_init, E, *, dt, n_steps, wp_in, surface, position, wp_out, n_channels=1, order=3)[source]¶

{“flow”: sigma_DA, “delta”: sigma_DA, “tw”: sigma_DA} (each bohr^2, per anion channel, same shape convention as td_da_cross_section) from ONE shared propagation – the DA sibling of td_ve_cross_sections_all.

Flux, Dirac, and TannorWeeks (all axis=”nuclear”) are built up front and driven by a SINGLE propagate(…, extractors=[flux, dirac, tw]) call – identical dynamics psi(t_n) feed all three transforms, so any spread between the returned cross sections reflects a genuine difference between the energy-extraction methods (or, at an under- converged grid/propagation length, a shared discretization/truncation residual all three inherit together – see docs/physics/td-da.md), never a difference in what was propagated. No elastic free-reference run (DA has no elastic diagonal, see td_da_cross_section).

surface/position/wp_out are all REQUIRED here (unlike td_da_cross_section, which only needs whichever one its method selects).

Parameters:

tgrid (TensorGrid)
model (ResonanceModel)
eps (npt.NDArray[np.float64])
chi (npt.NDArray[np.complex128])
v_init (int)
E (float | npt.ArrayLike)
dt (float)
n_steps (int)
wp_in (_WpIn)
surface (int)
position (int)
wp_out (_WpOut)
n_channels (int)
order (int)

Return type:

dict[str, npt.NDArray[np.float64]]

Grids¶

qscat.core.electronic_grid(*, r_max=30.0, angle_deg=35.0, order=8, n_complex=8, tail_alpha=0.2, tail_skip=2)[source]¶

Electronic radial grid: real region [0, r_max] + an ECS tail at angle_deg.

The ECS pivot is R0 == r_max by construction.

Parameters:

r_max (float, optional) – Outer edge of the real region (bohr); the ECS pivot. Must exceed the fixed inner segments.
angle_deg (float, optional) – ECS rotation angle of the complex tail, in degrees.
order (int, optional) – DVR quadrature order (points per element).
n_complex (int, optional) – Number of complex (ECS-tail) elements.
tail_alpha (float, optional) – Geometric growth factor of the tail element lengths.
tail_skip (int, optional) – Number of leading tail elements before growth begins.

Returns:

The assembled electronic radial grid.

Return type:

FemDvrEcsGrid

Raises:

GridError – If r_max does not exceed the fixed inner segments.

qscat.core.nuclear_grid(*, angle_deg=35.0, r_max=40.0, n_complex=10, quadrature=14)[source]¶

Build the nuclear FEM-DVR-ECS grid: real region [0, 12] + ECS tail.

Real elements tile [0, 12] bohr in four segments (see _REAL_SEGMENTS), finely resolved (0.15 bohr) around the equilibrium bond length so the low vibrational states converge. n_complex uniform complex elements at angle_deg tile [12, r_max]. The ECS pivot R0 == 12.0 bohr by construction (GridSpec computes R0 as x_min plus the sum of real element lengths).

Parameters:

angle_deg (float, optional) – ECS rotation angle of the complex tail, in degrees.
r_max (float, optional) – Outer edge of the complex tail (bohr). Must be finite and strictly greater than the 12.0 bohr real-region endpoint, so the tail spans a positive length.
n_complex (int, optional) – Number of complex (ECS-tail) elements tiling [12, r_max]. Must be an integer of at least 1: without a complex element the grid carries no ECS tail, hence no outgoing boundary condition.
quadrature (int, optional) – DVR quadrature order (points per element).

Returns:

The assembled nuclear radial grid.

Return type:

FemDvrEcsGrid

Raises:

GridError – If n_complex is not an integer of at least 1, or r_max is not finite and beyond the real-region endpoint.

qscat.core.fem_grid_exp_tail(real_segments, *, angle_deg, quadrature, tail_n, tail_alpha=0.2, tail_skip=2, x_min=0.0)[source]¶

A FEM-DVR-ECS grid whose real region is uniform-per-segment (like segmented_grid) but whose ECS tail grows EXPONENTIALLY rather than tiling uniformly – the H2+ electronic/nuclear deck layout (eMoScat’s uniform_increment/exp element spec), which segmented_grid cannot express.

real_segments are (n_elements, endpoint) pairs from x_min, tiled exactly like segmented_grid’s real region. The tail appends tail_n complex ElementSpec`s at `angle_deg whose lengths are _ecs_tail(base, tail_n, skip=tail_skip, alpha=tail_alpha), with base the length of the LAST real element. The ECS pivot R0 is the last real endpoint.

Parameters:

real_segments (Sequence[tuple[int, float]])
angle_deg (float)
quadrature (int)
tail_n (int)
tail_alpha (float)
tail_skip (int)
x_min (float)

Return type:

FemDvrEcsGrid

qscat.core.segmented_grid(real_segments, complex_segments, *, angle_deg, quadrature, x_min=0.0)[source]¶

A FEM-DVR-ECS grid from eMoScat’s grids.txt segment format.

real_segments / complex_segments are (n_elements, endpoint) pairs: from x_min, each segment tiles n uniform elements up to endpoint. The complex part is an ECS tail at angle_deg; the ECS pivot R0 is the last real endpoint. complex_segments may be empty (a pure real grid). This is the per-molecule discretisation route – see docs/physics/diatomic-ve-cross-sections.md (DA nuclear grids).

Parameters:

real_segments (Sequence[tuple[int, float]])
complex_segments (Sequence[tuple[int, float]])
angle_deg (float)
quadrature (int)
x_min (float)

Return type:

FemDvrEcsGrid

qscat.core.ecs_angle_family(electronic, nuclear, *, electronic_angles, nuclear_angles)[source]¶

The three tensor grids qscat.core.exact_resonance_states needs.

That function takes a base grid plus two grids differing from it in exactly one ECS angle each, and the correctness of its whole answer rests on that “exactly one” – which is easy to get wrong by hand and, until this builder, was got right by copy-paste at five separate call sites. Build the family here instead:

base, moved_el, moved_nu = ecs_angle_family(
    lambda a: fem_grid_exp_tail(segments, angle_deg=a, quadrature=8, tail_n=25),
    lambda a: nuclear_grid(angle_deg=a, r_max=14.0, n_complex=3, quadrature=8),
    electronic_angles=(30.0, 40.0),
    nuclear_angles=(18.0, 12.0),
)
res = exact_resonance_states(model, base, moved_el, moved_nu, ...)

Parameters:

electronic (callable) – angle_deg -> FemDvrEcsGrid. Called twice each, with the two angles; everything else about the grid must be identical between the two calls, which is exactly what a one-argument closure guarantees.
nuclear (callable) – angle_deg -> FemDvrEcsGrid. Called twice each, with the two angles; everything else about the grid must be identical between the two calls, which is exactly what a one-argument closure guarantees.
electronic_angles (tuple of float) – (base, moved) ECS angles in degrees. The two must differ.
nuclear_angles (tuple of float) – (base, moved) ECS angles in degrees. The two must differ.

Returns:

(base, electronic_moved, nuclear_moved), each [electronic, nuclear].

Return type:

tuple of TensorGrid

Raises:

GridError – If either angle pair is degenerate, or if a builder returns grids that do not share their real nodes.

qscat.core.assert_shared_real_nodes(grid_a, grid_b, *, what='grid_a and grid_b')[source]¶

Reject two grids that do not share every real node.

A two-angle ECS stability test compares eigenvalues of two discretizations that must differ ONLY in their tail angle. A different real-region mesh makes the residuals meaningless – and, worse, meaningless in the flattering direction: real-region discretization error no longer cancels between the two spectra, so states that should match do not, and the test silently becomes a convergence check rather than a stability check.

Raises GridError naming what, so a mismatch surfaces here rather than as a downstream shape error.

Parameters:

grid_a (FemDvrEcsGrid)
grid_b (FemDvrEcsGrid)
what (str)

Return type:

None

Channels¶

qscat.core.channel_vector(tgrid, k, chi_v, l, *, charge=0)[source]¶

DVR coefficients of F_{E,l}(r) chi_v(R), masked to the unscaled region.

chi_v is already a coefficient vector; F is a function and picks up sqrt(w_r). charge=0 (default) is the neutral-target case: the electronic factor is the fast riccati_bessel_en free radial function, bit-for-bit unchanged from before charge existed (every VE/DA caller). charge != 0 is the ionic-target case (e.g. H2+): the electronic factor is the energy-normalized Coulomb function coulomb_f_en instead (mass-1 electron).

Parameters:

tgrid (TensorGrid)
k (float)
chi_v (NDArray[complex128])
l (int)
charge (int)

Return type:

NDArray[complex128]

qscat.core.anion_electronic_states(g_r, model, R_inf, n_states=1)[source]¶

Anion bound electronic state(s) at the dissociation limit R_inf.

Diagonalizes -1/2 d^2/dr^2 + surface(r, R_inf) (electron mass 1) on the electronic grid. surface is the FULL electronic potential at R_inf (v0(R_inf) + ell(ell+1)/2r^2 + v_int(r, R_inf)), so eps_e shares the H_2D energy scale (it includes v0(R_inf)) and the DA threshold eps_e - eps[v_init] is correct.

Returns (eps_e, phi_e): eps_e the n_states lowest-Re eigenvalues with |Im(E)| < _IM_TOL_HA AND Re(E) < v0(R_inf) (the genuinely bound states), real, ascending; phi_e shape (n_states, g_r.n), each c-product-normalized over the electronic real region. Raises ValueError if fewer than n_states bound states exist (e.g. n_states reached past the finite bound spectrum).

Parameters:

g_r (FemDvrEcsGrid)
model (ResonanceModel)
R_inf (float)
n_states (int)

Return type:

tuple[npt.NDArray[np.float64], npt.NDArray[np.complex128]]

qscat.core.v_dr_diag(tgrid, model)[source]¶

The rearrangement interaction V_DR = V_int(r,R) + v0(R) - V_int(r, R_inf), flat (C-order), length tgrid.size. R_inf = tgrid.grids[1].R0 (the nuclear ECS pivot / real-region endpoint, eMoScat’s nu_inf).

This – not V_int – is the operator in the DA/DR T-matrix: H - H_final, where H_final is the asymptotic channel Hamiltonian (electron bound in V_int(r, R_inf), free nuclei on v0). As R -> R_inf the V_int terms cancel and V_DR -> v0(R).

Parameters:

tgrid (TensorGrid)
model (ResonanceModel)

Return type:

npt.NDArray[np.complex128]

qscat.core.outgoing_channel(tgrid, chi_v, *, r0_out, p0_out, sigma_out)[source]¶

Phi_{v’} = g_out(r) chi_{v’}(R), masked, energy-independent.

chi_v is already a c-product-normalized DVR coefficient vector (see vibrational_states’s docstring); no rescaling is applied here.

Parameters:

tgrid (TensorGrid)
chi_v (NDArray[complex128])
r0_out (float)
p0_out (float)
sigma_out (float)

Return type:

NDArray[complex128]

qscat.core.outgoing_channel_nuclear(tgrid, phi_c, *, r0_out, p0_out, sigma_out)[source]¶

Phi_c = phi_c(r) g_out(R), masked, energy-independent.

The NUCLEAR-axis transpose of outgoing_channel: the outgoing test packet g_out sits in the nuclear coordinate R (tgrid.grids[1]) while phi_c – one of the anion electronic bound states (qscat.core.dissociation.anion_electronic_states) – sits in the electronic coordinate r (tgrid.grids[0]). phi_c is already c-product-normalized (that function’s docstring); no rescaling is applied here, mirroring outgoing_channel’s treatment of chi_v.

Parameters:

tgrid (TensorGrid)
phi_c (NDArray[complex128])
r0_out (float)
p0_out (float)
sigma_out (float)

Return type:

NDArray[complex128]

Vibrational structure¶

qscat.core.vibrational_states(grid, mu, n, v0)[source]¶

The n lowest bound eigenpairs of T_nuc(mu) + diag(v0(R)).

v0 is the neutral molecule’s potential-energy curve, evaluated pointwise on grid.points – passed in by the caller (e.g. qscat.model.N2.v0) rather than hardcoded, so this solver stays model-independent. The eigenvectors are Hermitian-normalized (v^dagger v = 1) by qscat.dvr.eigen(), which coincides with the c-product norm (v^T v = 1) for these real bound-state vectors.

Parameters:

grid (FemDvrEcsGrid) – The nuclear radial grid.
mu (float) – Nuclear reduced mass (atomic units).
n (int) – Number of lowest bound vibrational states to return.
v0 (callable) – The neutral potential-energy curve v0(R), evaluated pointwise on grid.points.

Returns:

A (eps, chi) named tuple: eps the real part of the n lowest-Re(E) eigenvalues (Hartree), ascending; chi the (n, grid.n) array of eigenvectors, one row per level.

Return type:

VibrationalBasis

Raises:

GridError – If any of the n lowest-Re(E) eigenvalues has |Im(E)| > _IM_TOL_HA – i.e. n reached past the true bound spectrum into quasi-continuum/ECS states, which are not valid vibrational levels.

class qscat.core.VibrationalBasis(eps, chi)[source]¶

The vibrational eigenbasis returned by vibrational_states.

A NamedTuple, so the historical eps, chi = vibrational_states(…) unpacking keeps working unchanged, while .eps / .chi give named access.

Parameters:

eps (NDArray[float64])
chi (NDArray[complex128])

eps: NDArray[float64]¶: Alias for field number 0

chi: NDArray[complex128]¶: Alias for field number 1

Wavepacket and correlation¶

qscat.core.gaussian_coeffs(grid, *, r0, p0, sigma)[source]¶

DVR coefficients of g(r) on grid, zero on the ECS tail.

Parameters:

grid (FemDvrEcsGrid)
r0 (float)
p0 (float)
sigma (float)

Return type:

NDArray[complex128]

qscat.core.initial_state(tgrid, chi_v, *, r0, p0, sigma)[source]¶

Psi(0) = g(r) chi_v(R), flat, masked, renormalized to unit Hermitian L2 norm.

Renormalization uses np.linalg.norm (the true sqrt(sum |psi_j|^2) probability norm), NOT the c-product self-pairing: for a wavepacket carrying a momentum phase exp(i p0 r), c_product(g, g) is a small oscillatory complex number, not a norm, and using it here would silently produce a state far from unit probability. The c-product is reserved for the chi_v self-pairing below (its ECS-basis normalization convention, per vibrational_states’ docstring) and, downstream, for correlation functions and the S-matrix – never for this state’s overall scale.

Parameters:

tgrid (TensorGrid)
chi_v (NDArray[complex128])
r0 (float)
p0 (float)
sigma (float)

Return type:

NDArray[complex128]

qscat.core.eta_incident(grid, k, l, *, r0, p0, sigma)[source]¶

eta_in(E) = c_product(g_in_coeffs, F_{E,l}_coeffs) on the electronic grid.

Parameters:

grid (FemDvrEcsGrid)
k (float)
l (int)
r0 (float)
p0 (float)
sigma (float)

Return type:

complex

qscat.core.eta_outgoing(grid, kp, l, *, r0_out, p0_out, sigma_out, mass=1.0)[source]¶

eta_out(E’) = c_product(g_out_coeffs, F^out_{E’,l}_coeffs) on grid.

F^out is the outgoing Hankel half, NOT the regular function – see module docstring. mass defaults to 1.0 (electronic, byte-identical to the pre-mass code – see _outgoing_coeffs’s docstring); a nuclear DA caller (grid = tgrid.grids[1]) passes mass=model.mu, l=0.

Parameters:

grid (FemDvrEcsGrid)
kp (float)
l (int)
r0_out (float)
p0_out (float)
sigma_out (float)
mass (float)

Return type:

complex

qscat.core.hankel_point_value(z_position, k, l=0, charge=0, *, mass=1.0)[source]¶

H^{(1)}_{E,l}(z_position)/2 – the outgoing-Hankel-half VALUE at a single physical (real, unscaled) coordinate, e.g. z_position = grid.real_points[position] for some fixed DVR index position (Dirac’s analysis point, td_extractors.py).

The scalar sibling of _outgoing_coeffs: same energy-normalized outgoing function – riccati_hankel_en_mass(z_position, k, l, mass)/2 (neutral, charge == 0) or coulomb_h1_en(z_position, k, charge, mass, l)/2 (charged target) – but evaluated at ONE point rather than converted to a sqrt(w_r)-scaled, masked DVR coefficient VECTOR: a delta-distribution test function needs the outgoing function’s VALUE, not an integral against it.

mass defaults to 1.0 (the electronic reduced mass, a.u.) – every existing (electronic) call site is untouched: riccati_hankel_en_mass( …, 1.0) reproduces riccati_hankel_en(…) bit-for-bit (2.0*1.0 == 2.0 exactly), and coulomb_h1_en(…, 1.0, l) is the same literal 1.0 the pre-mass code passed. A nuclear (dissociation) caller passes mass=model.mu.

Parameters:

z_position (float)
k (float)
l (float | int)
charge (float | int)
mass (float)

Return type:

complex

qscat.core.outgoing_surface_wave(z_surface, k, l, charge=0.0, *, mass=1.0)[source]¶

(phi_out, dphi_out) = H^{(1)}_{E,l}(z_surface)/2 and its SPATIAL derivative at z_surface – the Flux extractor’s per-channel outgoing wave + derivative (eMoScat FluxTestFunction2d’s phi_out_/dphi_out_, confirmed by port-scout reading FluxTestFunction2d.cpp’s constructor: it samples sphHankel1En(…)/2 – or the Coulomb sH1_en(…)/2 for a charged target – over an element’s nodes and applies the SAME DVR derivative GridVector::derivative uses; that machinery is reproduced here directly against the ANALYTIC function instead, see below).

mass (eMoScat’s reduced_mass(), mu_x_=1.0 electronic / mu_y_=mass nuclear) is the mass entering the energy normalization of the outgoing function – F^{(1)}_{E,l}(r) = sqrt(2 mass k/pi) r h_l^{(1)}(kr) (riccati_hankel_en_mass’s definition; mass does NOT enter the momentum argument kr, matching riccati_bessel_en_mass’s convention). Defaults to 1.0 – every existing (electronic) call site is untouched: at mass=1.0 the formulas below reproduce the pre-mass code bit-for- bit (2.0*1.0 == 2.0 exactly). A nuclear (dissociation) caller passes mass=model.mu.

Neutral (charge == 0): the VALUE is computed via hankel_point_value (same multiplication order as the old inline formula – bit-identical); only the DERIVATIVE is computed analytically here via the product rule,

dF/dr = sqrt(2 mass k/pi) * [h_l(kr) + kr * h_l’(kr)]

using scipy.special.spherical_jn/spherical_yn’s derivative=True option for h_l^{(1)}{}’(x) = j_l’(x) + i y_l’(x) – qscat.special. radial does not itself expose a derivative primitive, but the underlying scipy pieces it is built on already support one, so no finite difference is needed for this branch (checked against a finite difference of riccati_hankel_en/riccati_hankel_en_mass in test_correlation.py, confirming the analytic formula).

Charged (charge != 0, coulomb_h1_en): qscat.special.coulomb has no derivative primitive for the Coulomb functions (mpmath’s coulombf/ coulombg expose no derivative= option, and the F_l’/G_l’ recurrence needs extra pieces this module does not carry) – this branch falls back to a high-accuracy CENTRAL finite difference (4th-order, 5-point stencil) of coulomb_h1_en(…, mass, l) itself. Kept structurally for a charged target (e.g. H2+); N2/F2 are neutral, so only the analytic branch is exercised by the neutral-molecule gates.

Parameters:

z_surface (float)
k (float)
l (int)
charge (float)
mass (float)

Return type:

tuple[complex, complex]

qscat.core.propagate(tgrid, psi0, out_channels, *, dt, n_steps, sample_period=0, snapshot_times=None, keep_psi_at=None, hamiltonian, order=3, extractors=None)[source]¶

Propagate and sample. See module docstring for the two cadences.

hamiltonian is the propagation Hamiltonian (REQUIRED here – the model-agnostic engine has no default to fall back to; the caller supplies model.hamiltonian(tgrid), or model.hamiltonian(tgrid) - diag(model.interaction_diag(tgrid)) for the elastic free-reference path, see propagate_wavepacket’s free argument).

order is the diagonal-Pade order of the evolution operator (qscat.evolution.make_pade_stepper): O(dt^(2*order+1)) per step. The default 3 matches eMoScat (order-1 Crank-Nicolson under-converges badly over a multi-thousand-step run – ~100% accumulated error at dt=0.5/1.0 – and is the reason an order-1 TD cross section only reached ~10-15% of the TI oracle; order 3 brings it to convergence). See docs/physics/n2-2d-td-cross-section.md.

This is the propagate-ONCE engine: the trajectory is computed a single time and ex.record(psi) is called on every extractor in extractors (see the Extractor protocol) at every step, alongside the legacy out_channels correlation bookkeeping that fills PropagationResult.c (kept for the existing callers/tests that read .c directly – e.g. the N2 project’s td_propagation/td_cross_section shims and qscat_run.artifacts, the serialization layer). Pass out_channels=[] when only extractors are wanted (e.g. td_ve_cross_section’s method=”tw” route).

Parameters:

tgrid (TensorGrid)
psi0 (NDArray[complex128])
out_channels (list[NDArray[complex128]])
dt (float)
n_steps (int)
sample_period (int)
snapshot_times (list[float] | None)
keep_psi_at (list[float] | None)
hamiltonian (spmatrix)
order (int)
extractors (list[Extractor] | None)

Return type:

PropagationResult

qscat.core.sigma_from_correlations(tgrid, model, result, eps, v_init, vprimes, E, *, dt, wp_in, wp_out, free_result=None)[source]¶

sigma_{v_init->v’}(E) (bohr^2) from an ALREADY-COMPUTED PropagationResult.

The “cheap sigma(E) from a stored propagation” primitive: unlike td_ve_cross_section, this never runs (or re-runs) the Pade propagation – it only transforms the c_{v’}(t) already sitting in result (e.g. loaded back from a saved .npz, or a truncated/reused trajectory from a convergence study).

E (collision energy, Hartree) may be scalar or array-like; scalar E returns shape (len(vprimes),), array E returns (len(E), len(vprimes)) – the SAME convention as td_ve_cross_section and qscat.core.driven.ve_cross_section.

dt, wp_in, wp_out must match the values used to produce result (the quadrature step and the incident/outgoing wavepacket parameters that eta_incident/eta_outgoing are evaluated with) – this function does not validate that consistency, it trusts the caller.

free_result is the V_int=0 reference propagation (same wavepacket/grid, from propagate_wavepacket(…, free=True)); when supplied, the diagonal/elastic channel subtracts its S_free(E) instead of a literal 1 (see sigma_one_energy). Leave None to reproduce the old behavior (correct for the inelastic channels; the elastic channel then needs S_free -> 1, which this transform does not satisfy).

Parameters:

tgrid (TensorGrid)
model (ResonanceModel)
result (PropagationResult)
eps (npt.NDArray[np.float64])
v_init (int)
vprimes (list[int])
E (float | npt.ArrayLike)
dt (float)
wp_in (_WpIn)
wp_out (_WpOut)
free_result (PropagationResult | None)

Return type:

npt.NDArray[np.float64]

Time-dependent energy extractors¶

All three share one propagate-once protocol, so a single propagation can drive every extractor. See docs/physics/td-extractors.md.

class qscat.core.Extractor(*args, **kwargs)[source]¶

A recorder+transform pair driven by one shared propagate trajectory.

record is called with the current psi(t_n) at EVERY propagation step (accumulating whatever per-step datum the extractor needs – e.g. TannorWeeks in td_extractors.py appends c_product(Phi_v’, psi) per v’); sigma transforms the accumulated series into a cross section, shape (len(E), len(vprimes)). Letting propagate drive a LIST of these means one Pade trajectory can feed several alternative energy-extraction routes (Tannor-Weeks, and the delta/flow methods – TannorWeeks, Dirac, Flux in td_extractors.py) without re-propagating.

record(psi)[source]¶

Accumulate this step’s datum from the current psi(t_n).

Called at every propagation step.

Parameters:: psi (NDArray[complex128])
Return type:: None

sigma(E, *, free=None, n_steps=None)[source]¶

Transform the recorded series into a cross section.

Parameters:

E (float | ArrayLike) – Energy or energies at which to compute the cross section.
free (Extractor | None) – Free-reference extractor for the elastic subtraction.
n_steps (int | None) – Truncate the transform to the first n_steps recorded samples. A convergence probe; both runs must share the step schedule.

Return type:

Cross section with shape (len(E), n_channels).

class qscat.core.TannorWeeks(tgrid, model, eps, chi, v_init, vprimes, wp_out, *, wp_in, dt, axis='electronic', n_channels=1)[source]¶

The Tannor-Weeks Extractor: records c_{v’}(t), transforms via the existing eta-deconvolution + elastic free-reference logic.

Construct one instance per propagation run (the full run, and – when subtract_free_reference applies – a second instance for the V_int=0 free-reference run); pass it to time_dependent.propagate(…, extractors=[tw]). sigma(E, free=…) accepts the companion free- reference TannorWeeks (built from ITS OWN, separately propagated, run) so the elastic channel can subtract S_free(E) exactly as sigma_one_energy documents – this is the “second TannorWeeks” the time_dependent.propagate design notes describe.

Nuclear (axis=”nuclear”): the DISSOCIATIVE ATTACHMENT (DA) sibling of Flux(axis=”nuclear”)/Dirac(axis=”nuclear”) (module docstring’s TannorWeeks(axis=”nuclear”) section) – wp_out is now the NUCLEAR outgoing test packet’s parameters (r0_out/p0_out/sigma_out in R), n_channels selects how many anion electronic bound states (qscat.core.dissociation.anion_electronic_states, at R_inf = tgrid.grids[1].R0) are tracked as exit channels; vprimes is unused (pass []). record is UNCHANGED (the same per-channel c-product loop – only the channel test functions differ). sigma(E) returns the DA transform’s sigma_DA,c(E) per channel, shape (n_channels,) for scalar E (matching Flux/Dirac(axis=”nuclear”)’s contract). free is not supported (DA has no elastic diagonal to subtract a reference from; passing it raises ValueError).

Parameters:

tgrid (TensorGrid)
model (ResonanceModel)
eps (npt.NDArray[np.float64])
chi (npt.NDArray[np.complex128])
v_init (int)
vprimes (list[int])
wp_out (_WpOut)
wp_in (_WpIn)
dt (float)
axis (Axis)
n_channels (int)

record(psi)[source]¶

Append this step’s c_{v’}(t_n) = c_product(Phi_{v’}, psi) row.

Same per-channel c-product, same channel order, as the legacy out_channels loop in time_dependent.propagate – so a run driven purely through this extractor reproduces PropagationResult.c to machine precision.

Parameters:: psi (NDArray[complex128])
Return type:: None

property result: PropagationResult¶: The recorded series as a PropagationResult (norm/snapshots are not tracked by this extractor and are left empty/zero – only t/c feed sigma_from_correlations (electronic axis) or _tw_da_sigma_one_energy (nuclear axis)). The full series (self._result(None)); see sigma’s n_steps for a truncated read.

sigma(E, *, free=None, n_steps=None)[source]¶

sigma_{v_init->v’}(E) (bohr^2, electronic axis) or sigma_DA,c(E) (bohr^2 per anion channel c, nuclear axis) via the Tannor-Weeks transform (class docstring).

Electronic: free, when given, must be a companion TannorWeeks extractor recorded from a SEPARATE V_int=0 propagation (same wavepacket/grid); its series supplies sigma_from_correlations’s free_result – the elastic (v’==v_init) channel then subtracts the free-particle S_free(E) instead of a literal 1 (see time_dependent.sigma_one_energy). Leave None to reproduce the literal-1 fallback. free is typed as the general Extractor protocol so all three extractors share one sigma signature, but only a TannorWeeks is meaningful here; anything else raises TypeError.

Nuclear: DA has no elastic diagonal to subtract a reference from – free must be None (ValueError otherwise).

n_steps (keyword-only, default None): transform only the FIRST n_steps recorded samples instead of the full series – None reproduces today’s behavior byte-identically (self._rows[:None] == self._rows[:]). When free is also given, ITS series is truncated to the SAME n_steps (both runs share the same step schedule).

Parameters:

E (float | ArrayLike)
free (Extractor | None)
n_steps (int | None)

Return type:

NDArray[float64]

class qscat.core.Dirac(tgrid, model, eps, chi, v_init, vprimes, position, *, wp_in, dt, axis='electronic', n_channels=1)[source]¶

The delta (Dirac) Extractor: eMoScat’s DiracTestFunction2d – Tannor-Weeks with a delta-distribution test function instead of the Gaussian test packet.

Electronic (axis=”electronic”, the default): record projects psi onto chi_{v’} in the nuclear coordinate at a FIXED electronic DVR index position – no propagated outgoing test function, unlike TannorWeeks. position must land in the real (unscaled) electronic region (grid.real_points[position] <= grid.R0), typically an element end well past the interaction range (mirroring TannorWeeks’s wp_out standoff). sigma(E, free=…) accepts the companion free-reference Dirac (built from ITS OWN, separately propagated, run), matching TannorWeeks.sigma’s free contract.

Nuclear (axis=”nuclear”): the DISSOCIATIVE ATTACHMENT (DA) sibling of Flux(axis=”nuclear”) (module docstring’s Flux(axis=”nuclear”) section) – position is a NUCLEAR DVR index, n_channels selects how many anion electronic bound states (qscat.core.dissociation. anion_electronic_states, at R_inf = tgrid.grids[1].R0) are tracked as exit channels; vprimes is unused (pass []). record projects psi onto phi_c(r) (the electronic anion state) at the FIXED nuclear node position – the SAME point-VALUE projection as the electronic path, transposed to the other axis (no derivative, unlike Flux: a delta test function needs only the point value). sigma(E) returns the delta transform’s sigma_DA,c(E) per channel (module docstring), eta_out -> hankel_point_value(…, mass=model.mu) at position – shape (n_channels,) for scalar E, matching Flux(axis=”nuclear”)’s contract. free is not supported (DA has no elastic diagonal to subtract a reference from; passing it raises ValueError).

Construct one instance per propagation run (the full run, and – when an elastic free reference is wanted, electronic axis only – a second instance for a SEPARATE V_int=0 run); pass it to time_dependent. propagate(…, extractors=[dirac]).

position must be a real (unscaled) DVR index, on the axis-appropriate grid, in the asymptotic region (past the interaction) – same requirement as Flux’s surface.

Parameters:

tgrid (TensorGrid)
model (ResonanceModel)
eps (npt.NDArray[np.float64])
chi (npt.NDArray[np.complex128])
v_init (int)
vprimes (list[int])
position (int)
wp_in (_WpIn)
dt (float)
axis (Axis)
n_channels (int)

record(psi)[source]¶

Append this step’s point-value projection row (class docstring): b_{v’}(t_n) = <chi_{v’}|psi(position, .)> / sqrt(w_r[position]) (electronic axis) or b_c(t_n) = <phi_c|psi(., R=position)> / sqrt(w_R[position]) (nuclear axis) – the DVR-coefficient-to-VALUE conversion qscat.core.lcp’s boundary flux also relies on.

Parameters:: psi (NDArray[complex128])
Return type:: None

property result: PropagationResult¶: The recorded series as a PropagationResult (same shape/role as TannorWeeks.result; norm/snapshots are left empty/zero). The full series (self._result(None)); see sigma’s n_steps for a truncated read.

sigma(E, *, free=None, n_steps=None)[source]¶

sigma_{v_init->v’}(E) (bohr^2, electronic axis) or sigma_DA,c(E) (bohr^2 per anion channel c, nuclear axis) via the delta transform (module docstring / class docstring).

Electronic: free, when given, must be a companion Dirac extractor recorded from a SEPARATE V_int=0 propagation – same elastic free-reference contract as TannorWeeks.sigma. free is typed as the general Extractor protocol (one shared sigma signature) but only a Dirac is meaningful here; anything else raises TypeError.

Nuclear: DA has no elastic diagonal to subtract a reference from – free must be None (ValueError otherwise).

n_steps (keyword-only, default None): transform only the FIRST n_steps recorded samples – see TannorWeeks.sigma’s n_steps for the exact contract (byte-identical to today when None; free’s series is truncated to the SAME n_steps when both are given).

Parameters:

E (float | ArrayLike)
free (Extractor | None)
n_steps (int | None)

Return type:

NDArray[float64]

class qscat.core.Flux(tgrid, model, eps, chi, v_init, vprimes, surface, *, wp_in, dt, axis='electronic', n_channels=1)[source]¶

The flow (flux) Extractor: eMoScat’s FluxTestFunction2d – the time-energy Fourier transform of the probability flux projected onto the outgoing channel, at a FIXED surface.

Electronic (axis=”electronic”, the default): record appends, per v’, BOTH the value b_{v’}(t) = <chi_{v’}|psi(surface,.)> (Dirac’s line projection) AND its electronic-coordinate derivative d_{v’}(t) = <chi_{v’}| d/dr psi(surface,.)> (via qscat.dvr. dvr_first_derivative_at_node applied along the electronic axis, then a nuclear c-product onto chi_{v’} – no extra /sqrt(w) needed there, unlike b_{v’}: dvr_first_derivative_at_node already converts coefficient->value internally, see its docstring). sigma(E, free=…) is the Wronskian-like flux transform (module docstring); same elastic free-reference contract as TannorWeeks.sigma/Dirac.sigma. vprimes is the VE exit-vibrational-level list (electronic-axis param only).

Nuclear (axis=”nuclear”): the DISSOCIATIVE ATTACHMENT (DA) extractor (module docstring’s Flux(axis=”nuclear”) section) – surface is a NUCLEAR DVR index, n_channels selects how many anion electronic bound states (qscat.core.dissociation.anion_electronic_states, at R_inf = tgrid.grids[1].R0) are tracked as exit channels; vprimes is unused (pass []). record appends b_c(t)/d_c(t) per anion channel c. sigma(E) returns sigma_DA,c(E) per channel, shape (n_channels,) for scalar E (matching dissociation.da_cross_section’s per-channel return contract) – free is not supported (DA has no elastic diagonal to subtract a reference from; passing it raises ValueError).

surface must be a real (unscaled) DVR index, on the axis-appropriate grid, in the asymptotic region (past the interaction) – same requirement as Dirac’s position.

Parameters:

tgrid (TensorGrid)
model (ResonanceModel)
eps (npt.NDArray[np.float64])
chi (npt.NDArray[np.complex128])
v_init (int)
vprimes (list[int])
surface (int)
wp_in (_WpIn)
dt (float)
axis (Axis)
n_channels (int)

record(psi)[source]¶

Append this step’s b(t_n)/d(t_n) rows (docstring): per v’ (electronic axis) or per anion channel c (nuclear axis).

Parameters:: psi (NDArray[complex128])
Return type:: None

property series: tuple[NDArray[float64], NDArray[complex128], NDArray[complex128]]¶: The full recorded flux series (t, boundary_value, derivative) – the raw material of the flux transform, the Flux analogue of TannorWeeks/Dirac’s .result (which carry .t/.c instead). Public so a caller can persist the raw series (e.g. the qscat-run CLI’s optional correlations artifact); see sigma’s n_steps for a truncated read.

sigma(E, *, free=None, n_steps=None)[source]¶

sigma_{v_init->v’}(E) (bohr^2, electronic axis) or sigma_DA,c(E) (bohr^2 per anion channel c, nuclear axis) via the flux transform (module docstring / class docstring).

Electronic: free, when given, must be a companion Flux extractor recorded from a SEPARATE V_int=0 propagation – same elastic free-reference contract as TannorWeeks.sigma/Dirac.sigma. free is typed as the general Extractor protocol (one shared sigma signature) but only a Flux is meaningful here; anything else raises TypeError.

Nuclear: DA has no elastic diagonal to subtract a reference from – free must be None (ValueError otherwise).

Parameters:

E (float | ArrayLike)
free (Extractor | None)
n_steps (int | None)

Return type:

NDArray[float64]

The LCP approximation¶

The local-complex-potential reduction and the Born–Oppenheimer resonance levels built on it. These are the approximation under test against the exact solvers above; see docs/physics/diatomic-ve-cross-sections.md and docs/physics/lcp-resonance-levels.md.

qscat.core.local_complex_potential(model, nuclear_grid, elec_grid_a, elec_grid_b, *, re_half_width=0.05, im_half_width=0.05, resid_tol=0.001)[source]¶

Local complex potential (V_d(R), Gamma(R)) for the LCP DA/VE model.

V_d(R) = Re(E_pole(R)), Gamma(R) = max(0, -2 Im(E_pole(R))), E_pole the two-angle-matched resonance pole of -1/2 d^2/dr^2 + model.surface(r,R) (surface includes v0(R), so V_d = Re(E_pole) directly). Seeded from the bound anion at R_inf = nuclear_grid.R0 (anion_electronic_states) and continued INWARD; small-R breakdown freezes the electronic shift; the complex tail clamps Gamma=0. See module docstring.

Parameters:

model (ResonanceModel)
nuclear_grid (FemDvrEcsGrid)
elec_grid_a (FemDvrEcsGrid)
elec_grid_b (FemDvrEcsGrid)
re_half_width (float)
im_half_width (float)
resid_tol (float)

Return type:

tuple[npt.NDArray[np.complex128], npt.NDArray[np.float64]]

qscat.core.lcp_da_cross_section(nuclear_grid: FemDvrEcsGrid, mu: float, Vd: NDArray[complex128], Gamma: NDArray[float64], eps: NDArray[float64], chi: NDArray[complex128], v_init: int, E: float | ArrayLike, *, ordering: Literal['NATURAL', 'MMD_ATA', 'MMD_AT_PLUS_A', 'COLAMD'] = 'COLAMD', return_wavefunction: Literal[False] = False) → NDArray[float64][source]¶

qscat.core.lcp_da_cross_section(nuclear_grid: FemDvrEcsGrid, mu: float, Vd: NDArray[complex128], Gamma: NDArray[float64], eps: NDArray[float64], chi: NDArray[complex128], v_init: int, E: float | ArrayLike, *, ordering: Literal['NATURAL', 'MMD_ATA', 'MMD_AT_PLUS_A', 'COLAMD'] = 'COLAMD', return_wavefunction: bool = False) → NDArray[float64] | tuple[NDArray[float64], NDArray[complex128] | None | list[NDArray[complex128] | None]]

LCP dissociative-attachment sigma_DA(E) (bohr^2), TI resolvent form.

Solve (E_tot I - H_res) psi_sc = d, H_res = T_nuc + diag(V_d - i Gamma/2), doorway d = sqrt(Gamma/2pi) chi_{v_init}; the DA amplitude is the outgoing dissociation flux at the boundary X (outermost real point): S_DA = sqrt(K/2pi mu) psi_sc(X), psi_sc(X) = psi_sc[b]/sqrt(w_b) (the wavefunction VALUE, not the DVR coefficient), sigma = 4 pi^3 |S_DA|^2/2E. The DA threshold eps_e = V_d(R_inf) = Vd[b].real (open iff E_tot > eps_e).

Requires the FINE per-molecule nuclear grid (the K~58 outgoing wave is unresolved on a coarse grid). The T->infty limit of eMoScat’s TD ModelLCP/SMatrix.cpp. The approximation under test vs the exact-2D da_cross_section oracle – validated at sigma_DA(F2,0.03)=1.47 vs ~1.66.

Argument-order note (docs/adr/0007): this solver deliberately takes (nuclear_grid, mu, Vd, Gamma, …) rather than a model – the LCP equation contains no model; its physics input IS the curve, which may come from resonance_levels(return_curve=True), a fit, or a file. ScatteringProblem.lcp_da_cross_section supplies mu/eps/chi/v_init from its bundle.

If return_wavefunction, also returns the 1-D nuclear resolvent psi_sc(R) = (E_tot I - H_res)^-1 d per energy (None when the DA channel is closed – E <= 0 or E_tot <= eps_e): one array for scalar E, one list entry per energy for array E, same convention as driven/dissociation. psi_sc is the DVR-coefficient vector on the full nuclear grid (length nuclear_grid.n).

Parameters:

nuclear_grid (FemDvrEcsGrid)
mu (float)
Vd (NDArray[complex128])
Gamma (NDArray[float64])
eps (NDArray[float64])
chi (NDArray[complex128])
v_init (int)
E (float | ArrayLike)
ordering (Literal['NATURAL', 'MMD_ATA', 'MMD_AT_PLUS_A', 'COLAMD'])
return_wavefunction (bool)

Return type:

NDArray[float64] | tuple[NDArray[float64], NDArray[complex128] | None | list[NDArray[complex128] | None]]

qscat.core.lcp_ve_cross_section(nuclear_grid: FemDvrEcsGrid, mu: float, Vd: NDArray[complex128], Gamma: NDArray[float64], eps: NDArray[float64], chi: NDArray[complex128], v_init: int, vprimes: list[int], E: float | ArrayLike, *, ordering: Literal['NATURAL', 'MMD_ATA', 'MMD_AT_PLUS_A', 'COLAMD'] = 'COLAMD', return_wavefunction: Literal[False] = False) → NDArray[float64][source]¶

qscat.core.lcp_ve_cross_section(nuclear_grid: FemDvrEcsGrid, mu: float, Vd: NDArray[complex128], Gamma: NDArray[float64], eps: NDArray[float64], chi: NDArray[complex128], v_init: int, vprimes: list[int], E: float | ArrayLike, *, ordering: Literal['NATURAL', 'MMD_ATA', 'MMD_AT_PLUS_A', 'COLAMD'] = 'COLAMD', return_wavefunction: bool = False) → NDArray[float64] | tuple[NDArray[float64], NDArray[complex128] | None | list[NDArray[complex128] | None]]

LCP vibrational-excitation sigma_{v_init->v’}(E) (bohr^2), TI resolvent form.

Provisional API (docs/adr/0004-public-api-stability-policy.md): this wide functional signature is the layer the context-object refactor targets and may change in a minor release.

Solve (E_tot I - H_res) xi = d_{v_init}, H_res = T_nuc(mu) + diag(V_d - i Gamma/2), doorway d_v = sqrt(Gamma/2pi) chi_v; S-matrix element S_{v’<-v_init} = <d_{v’}|xi> by the DVR c-product (no conjugate); sigma = 4 pi^3 |S|^2 / 2E, exactly zero for E <= 0 and for a closed final channel (E_tot - eps[v’] <= 0).

Graduated from projects/n2_ti_cross_section/cross_section.py’s ve_cross_section (the deliberately dense 1-D toy model). This version is SPARSE and sweep-reusing: A(E) = E_tot I - H_res has an E-independent sparsity pattern, so the symbolic analysis is done once and SparseLU.refactor re-runs only the numeric factor per energy – the same structure as lcp_da_cross_section and driven.ve_cross_section. xi depends only on (E, v_init), so one solve per energy serves every channel in vprimes.

If return_wavefunction, also returns xi(R) per energy (None when E <= 0): one array for scalar E, one list entry per energy for array E – the driven solution nuclear_density.lcp_driven_solution consumes.

Parameters:

nuclear_grid (FemDvrEcsGrid)
mu (float)
Vd (NDArray[complex128])
Gamma (NDArray[float64])
eps (NDArray[float64])
chi (NDArray[complex128])
v_init (int)
vprimes (list[int])
E (float | ArrayLike)
ordering (Literal['NATURAL', 'MMD_ATA', 'MMD_AT_PLUS_A', 'COLAMD'])
return_wavefunction (bool)

Return type:

NDArray[float64] | tuple[NDArray[float64], NDArray[complex128] | None | list[NDArray[complex128] | None]]

qscat.core.resonance_levels(model: ResonanceModel, nuclear_grid_a: FemDvrEcsGrid, nuclear_grid_b: FemDvrEcsGrid, elec_grid_a: FemDvrEcsGrid, elec_grid_b: FemDvrEcsGrid, *, re_half_width: float = 0.05, im_half_width: float = 0.05, resid_tol: float = 0.001, window: tuple[float, float, float, float] | None = None, n_levels: int | None = None, rel_tol: float = 0.0001, atol: float = 1e-08, golden_rule: bool = True, return_curve: Literal[False] = False) → ResonanceLevels[source]¶

qscat.core.resonance_levels(model: ResonanceModel, nuclear_grid_a: FemDvrEcsGrid, nuclear_grid_b: FemDvrEcsGrid, elec_grid_a: FemDvrEcsGrid, elec_grid_b: FemDvrEcsGrid, *, re_half_width: float = 0.05, im_half_width: float = 0.05, resid_tol: float = 0.001, window: tuple[float, float, float, float] | None = None, n_levels: int | None = None, rel_tol: float = 0.0001, atol: float = 1e-08, golden_rule: bool = True, return_curve: bool = False) → ResonanceLevels | tuple[ResonanceLevels, NDArray[complex128], NDArray[float64]]

Quasi-bound levels of model’s anion, straight from the model.

Runs the electronic pole walk ONCE (resonance_pole_walk, seeded from the asymptotic anion bound state exactly as local_complex_potential does), lays the resulting curve onto BOTH nuclear grids with _assemble_lcp, and diagonalizes (lcp_resonance_levels). E_res(R) at real R does not depend on the nuclear tail angle, so the second grid costs one extra nuclear diagonalization and nothing else.

nuclear_grid_b must share nuclear_grid_a’s real segments and quadrature and differ only in its ECS tail angle – conventionally a SMALLER angle, which is always safe against the model’s divergence bound.

If return_curve, also returns (Vd_a, Gamma) – the very curve the levels were computed in, on nuclear_grid_a. A caller that wants both the levels and the LCP curve (to solve lcp_da_cross_section in it, or to plot it under the levels) MUST take this route rather than calling local_complex_potential separately: that would repeat the expensive electronic walk AND, if any setting differed, report a curve that is not the one the levels came from.

Parameters:

model (ResonanceModel)
nuclear_grid_a (FemDvrEcsGrid)
nuclear_grid_b (FemDvrEcsGrid)
elec_grid_a (FemDvrEcsGrid)
elec_grid_b (FemDvrEcsGrid)
re_half_width (float)
im_half_width (float)
resid_tol (float)
window (tuple[float, float, float, float] | None)
n_levels (int | None)
rel_tol (float)
atol (float)
golden_rule (bool)
return_curve (bool)

Return type:

ResonanceLevels | tuple[ResonanceLevels, npt.NDArray[np.complex128], npt.NDArray[np.float64]]

qscat.core.lcp_resonance_levels(nuclear_grid_a, nuclear_grid_b, mu, Vd_a, Vd_b, Gamma, *, window=None, n_levels=None, rel_tol=0.0001, atol=1e-08, golden_rule=True)[source]¶

Quasi-bound levels of H_N = T(mu) + V_d(R) - i Gamma(R)/2.

The Born-Oppenheimer approximation to the 2-D model’s resonance energies: step 1 (the fixed-R electronic pole, local_complex_potential) supplies the complex curve; this is step 2, the nuclear eigenvalue problem in it. The thesis’s H_LCP (Vana 2017 Eq. 1.65).

nuclear_grid_a/nuclear_grid_b must share every real node and differ only in their ECS tail angle – that is what makes the two spectra comparable. Vd_a/Vd_b are the curve laid onto each grid (identical on the real nodes, differing in the continued tail); Gamma is real and tail-zero, so the same array serves both.

Physical levels are selected by two-angle stability (match_angle_stable); the rotated dissociative continuum fails that test and drops out. Levels with Im E > atol are unphysical and are dropped with a warning.

See docs/physics/lcp-resonance-levels.md.

Parameters:

nuclear_grid_a (FemDvrEcsGrid)
nuclear_grid_b (FemDvrEcsGrid)
mu (float)
Vd_a (NDArray[complex128])
Vd_b (NDArray[complex128])
Gamma (NDArray[float64])
window (tuple[float, float, float, float] | None)
n_levels (int | None)
rel_tol (float)
atol (float)
golden_rule (bool)

Return type:

ResonanceLevels

class qscat.core.ResonanceLevels(energies, widths, states, residuals, real_weight, golden_rule)[source]¶

The quasi-bound vibrational levels of the anion in the LCP curve.

The thesis’s omega_j (Vana 2017, Sec. 1.5/3.4), promoted from real levels of Re V_res to genuine complex eigenvalues. These are complex-scaled (ECS) resonance eigenstates – NOT Siegert pseudostates, which carry an outgoing-wave condition at a finite radius and a surface-corrected orthogonality relation (Hvizdos et al., Phys. Rev. A 97, 022704 (2018), App. A). ECS rotates rather than truncates, so the plain bilinear c-product is the complete inner product here.

energies: E_v - i Gamma_v/2 (Hartree), ascending in Re E.
widths: Gamma_v = -2 Im E_v (Hartree), UNCLAMPED – after the Im E <= atol physicality filter the most negative representable value is -2*atol; a small negative width is a round-off diagnostic, and hiding it behind a clamp is how it goes unnoticed. Same convention as ExactResonanceStates.widths. A level below the anion dissociation limit carries only the ELECTRONIC autodetachment width; one above it also carries a NUCLEAR (dissociative) width. Both come out of the one diagonalization.
states: shape (n_levels, grid.n) DVR COEFFICIENTS c_i (psi(R_i) = c_i / sqrt(w_i)), c-product-normalized: sum_i c_i^2 = 1.
residuals: the two-angle ECS-TAIL stability residual per level, |E_a - E_b| between the matched eigenvalues of the two rotation angles (qscat.ecs.match_angle_stable). A large residual means the level is contaminated by the rotated continuum, not a genuine pole. It is NOT a real-region convergence diagnostic: nuclear_grid_a and nuclear_grid_b are required to share every real node and quadrature, so real-region discretization error is common to both spectra and cancels out of the difference – residuals stays near machine precision even on a badly under-resolved real grid. Judge real-region convergence separately, by refining the shared real nodes and checking that energies itself does not move.
real_weight: fraction of |c|^2 inside the real region – a diagnostic, not a normalization. Near 1 for a well-localized level.
golden_rule: E_v^(0) - i Re<chi_v|Gamma|chi_v>/2, the perturbative comparator (the Gamma = 0 levels plus the first-order width). The expectation is taken with the bilinear c-product, so on an ECS grid it is complex in general and only its REAL part is a width; the discarded imaginary part is a tail-amplitude residue, negligible for a level localized in the real region and a sign that the comparator is inapplicable for one that is not. This is what eMoScat and the thesis actually computed. Agreement with energies means the level is perturbative; divergence means it is genuinely broad and the non-perturbative treatment is load-bearing. nan where no comparator level could be paired, and all-nan when golden_rule=False. The distance guard also produces nan when the Gamma-induced real shift exceeds half the local level spacing (the strongly non-perturbative regime, where nan is the honest answer) and when output levels are near-degenerate.

Parameters:

energies (NDArray[complex128])
widths (NDArray[float64])
states (NDArray[complex128])
residuals (NDArray[float64])
real_weight (NDArray[float64])
golden_rule (NDArray[complex128])

save(path)[source]¶

Write to an .npz archive under the dataclass’s own field names.

Shares its mechanism with ExactResonanceStates.save – see qscat.core._archive for why the pair round-trips through the dataclass’s own field names rather than a hand-rolled cache.

Parameters:: path (str | PathLike[str])
Return type:: None

classmethod load(path)[source]¶

Read back a save() file, checking every field is present.

Parameters:: path (str | PathLike[str])
Return type:: ResonanceLevels

Exact resonance states¶

Poles of the full 2-D S-matrix, with no Born–Oppenheimer separation and no local approximation — the objects resonance_levels above approximates. Identification is by ECS angle stability in both coordinates; see docs/physics/exact-2d-resonances.md.

qscat.core.exact_resonance_states(model, grid_base, grid_electronic, grid_nuclear, *, shifts, window, k=8, rel_tol=0.0001, atol=1e-08)[source]¶

Exact 2-D resonance states near shifts, by two-angle ECS stability.

Parameters:

model (ResonanceModel) – Supplies the sparse 2-D Hamiltonian via model.hamiltonian(tgrid).
grid_base (TensorGrid) – The tensor-product grid the returned energies and states belong to.
grid_electronic (TensorGrid) – The same grid with the ELECTRONIC ECS angle changed, and only that.
grid_nuclear (TensorGrid) – The same grid with the NUCLEAR ECS angle changed, and only that.
shifts (array_like of complex) – Seed shifts to search near. Physically motivated seeds are the Born-Oppenheimer levels from qscat.core.lcp.resonance_levels; they are passed in rather than computed here so that the exact solver never depends on the approximation it exists to measure.
window (tuple of float) – (re_lo, re_hi, im_lo, im_hi), the complex-plane box the angle-stability match is restricted to.
k (int, optional) – Eigenpairs per shift per grid (default 8).
rel_tol (float, optional) – Angle-stability acceptance: a state is kept when its partner satisfies |E_a - E_b| < max(rel_tol*|E_a|, atol), for BOTH partners.
atol (float, optional) – Angle-stability acceptance: a state is kept when its partner satisfies |E_a - E_b| < max(rel_tol*|E_a|, atol), for BOTH partners.

Returns:

Possibly empty. An empty result is a normal outcome – no angle-stable state near these shifts – and not an error.

Return type:

ExactResonanceStates

Raises:

ValueError – If window catches no eigenvalue at all in one of the three spectra, which means the window or the shifts are misplaced rather than that nothing is there.
GridError – If the three grids do not isolate exactly one ECS angle per partner – see _check_family. qscat.core.grids.ecs_angle_family builds a family that passes by construction.

Notes

Cost is one sparse factorization per (grid, shift): three grids times the number of seeds. On a production 2-D deck that factorization dominates everything else, so keep the seed list short and prefer the MUMPS backend.

class qscat.core.ExactResonanceStates(energies, widths, states, residual_electronic, residual_nuclear)[source]¶

Exact 2-D resonance states and the evidence that they are resonances.

Parameters:

energies (NDArray[complex128])
widths (NDArray[float64])
states (NDArray[complex128])
residual_electronic (NDArray[float64])
residual_nuclear (NDArray[float64])

energies¶

(m,). E_r - i*Gamma/2, ascending in Re E, taken from the base grid.

Type:: ndarray of complex128

widths¶

(m,). Gamma = -2 Im E.

Type:: ndarray of float64

states¶

(m, n). states[i] is the 2-D eigenvector of energies[i] on the base grid, flattened in TensorGrid order and c-product normalized – row-per-state, the same orientation as chi, anion_electronic_states, and ResonanceLevels.states.

Type:: ndarray of complex128

residual_electronic¶

(m,). |E_base - E_theta_r|: how far the eigenvalue moved when the ELECTRONIC ECS angle changed. Small means the state does not live in the electronic continuum.

Type:: ndarray of float64

residual_nuclear¶

(m,). |E_base - E_theta_R|: the same for the NUCLEAR angle.

Type:: ndarray of float64

Notes

Both residuals are reported rather than a single midpoint energy. With two independent angle partners there is no one midpoint, and the pair carries strictly more information: a state can be solidly bound in one coordinate and marginal in the other, which one number would hide.

A small residual is necessary evidence, not sufficient proof, and this is not a theoretical caveat: on the H2+ DR windows 4 of 57 poles that passed both angle tests turned out not to be resonances at all, scoring overlaps of 6e-4 to 7e-3 against a Born-Oppenheimer basis where genuine states score 0.87-0.99. Angle stability says the eigenvalue did not move when the contour did; a rotated-continuum state that happens to sit in a stable corner produces that too, and so can a grid too coarse to resolve the state.

Two separate checks answer the two separate questions. For “is this a resonance at all, and which quasi-bound state is it”, pair the state against a BO basis with qscat.core.assignment.pair_by_overlap. For “is the grid adequate”, refine it and see whether the energy moves. Neither is implied by the residuals here.

save(path)[source]¶

Write to an .npz archive under the dataclass’s own field names.

Shares its mechanism with ResonanceLevels.save – see qscat.core._archive for why the pair round-trips through the dataclass’s own field names rather than a hand-rolled cache (one call site once stored res_el/res_nuc where this dataclass calls them residual_electronic/residual_nuclear, a rename away from silently loading garbage).

Parameters:: path (str | PathLike[str])
Return type:: None

classmethod load(path)[source]¶

Read back a save() file, checking every field is present.

Also rejects an archive written before the row-per-state layout flip (states stored column-per-state) – a guard specific to this class, not shared with ResonanceLevels.load, which never had it.

Parameters:: path (str | PathLike[str])
Return type:: ExactResonanceStates

Born–Oppenheimer reference states¶

The picture an exact pole departs from: electronic eigen-curves over a nuclear grid, and the product states φ_j(r; R)·χ_v(R) built on them. electronic_curves covers the bound (ion / Rydberg) case and resonance_curve the resonant (neutral / anion) one, so one comparator serves both. See docs/physics/h2plus-resonance-states.md.

qscat.core.electronic_curves(model, g_r, g_R, *, n_curves, with_states=False)[source]¶

The n_curves lowest electronic eigen-curves over g_R.

At every nuclear grid point R (including the complex ECS-tail points, not just the real region) this diagonalizes the frozen-nucleus electronic problem -1/2 d^2/dr^2 + model.surface(r, R) (electron mass 1) on g_r and keeps the n_curves lowest-Re E eigenvalues. For an ion those are the Rydberg series; for a neutral, the bound electronic curves.

This reproduces qscat.core.anion_electronic_states wherever that function’s “genuinely bound” gate succeeds, and deliberately does NOT reuse that gate, because the gate fails on the ECS tail: a complex R_inf makes every eigenvalue there pick up an O(Im v0(R)) imaginary shift, already past the library’s 1e-6 bound-state tolerance a hair’s width into the tail. Gating on it would raise on perfectly good curve points.

Parameters:

model (ResonanceModel) – Supplies surface(r, R).
g_r (FemDvrEcsGrid) – Electronic and nuclear radial grids.
g_R (FemDvrEcsGrid) – Electronic and nuclear radial grids.
n_curves (int) – How many lowest-Re E curves to keep.
with_states (bool, optional) – Also return the eigenvectors, phase-aligned across R. Costs n_curves * g_r.n * g_R.n complex128 – ~100 MB on a production H2+ deck – so it is opt-in. Required for bo_basis to build product states; the energies alone are enough for a level table.

Return type:

ElectronicCurves

Raises:

GridError – If n_curves exceeds the electronic grid’s dimension.

qscat.core.resonance_curve(model, g_r_a, g_r_b, g_R, seed_window, *, re_half_width=0.05, im_half_width=0.05, resid_tol=0.001, with_states=True)[source]¶

The anion RESONANCE electronic curve over g_R, as a one-curve family.

The neutral-target counterpart of electronic_curves: the state a dissociative-attachment resonance is built on is not bound, so it cannot be picked out by “lowest Re E”. It is found instead as the angle-stable pole of the fixed-R electronic problem at two ECS angles (qscat.ecs.find_resonance_pole), continued inward over descending real R exactly as qscat.core.lcp.resonance_pole_walk does – and this function keeps the eigenvector that walk discards.

The grid-sizing sibling qscat.tuning.resonance.resonance_curve_arrays (renamed from resonance_curve in the 2026-08-25 API surface pass) runs the same underlying pole walk but returns plain (R, V_d, Gamma) arrays with no states – use that one to size a grid, this one to build BO basis states.

energies[0] is V_d(R) - i Gamma(R)/2, so qscat.core.lcp.local_complex_potential’s curve and this one are the same object viewed with and without its eigenfunctions.

The ECS tail is frozen, not walked. find_resonance_pole needs a real R. Tail points therefore take the outermost real node’s eigenvector and that node’s electronic shift added to v0(z) – the same analytic continuation qscat.core.lcp.curve._assemble_lcp applies to the potential, for the same reason. A product state’s tail behaviour is carried by chi_v, not by phi_res, so this is where the approximation belongs.

On breakdown the last accepted state is frozen inward, matching resonance_pole_walk’s behaviour rather than raising: a pole finder that loses the track at small R has usually run into a region where the resonance has merged with the continuum, and the frozen value is the honest continuation. ConvergenceError is raised only when the finder fails at the seed edge, where there is nothing to freeze.

Parameters:

model (ResonanceModel)
g_r_a (FemDvrEcsGrid) – Electronic grids differing ONLY in ECS tail angle. g_r_a is the one the returned eigenvectors live on.
g_r_b (FemDvrEcsGrid) – Electronic grids differing ONLY in ECS tail angle. g_r_a is the one the returned eigenvectors live on.
g_R (FemDvrEcsGrid) – Nuclear grid the curve is tabulated on.
seed_window (tuple of float) – (re_lo, re_hi, im_lo, im_hi) for the pole at the OUTERMOST real R. qscat.core.anion_electronic_states at g_R.R0 is the natural source.
re_half_width (float, optional) – Half-widths of the window recentred on each accepted pole.
im_half_width (float, optional) – Half-widths of the window recentred on each accepted pole.
resid_tol (float, optional) – Angle-stability residual above which the walk is considered broken.
with_states (bool, optional) – Keep the eigenvectors (default True – unlike electronic_curves, since a one-curve family is cheap and the states are the point).

Returns:

n_curves == 1.

Return type:

ElectronicCurves

Raises:

ConvergenceError – If no pole is accepted at the seed edge.

qscat.core.bo_basis(curves, g_R, mu, *, n_vib, allow_partial=False)[source]¶

Vibrational ladders in curves, and the product states they define.

Each curve is fed to qscat.core.vibrational_states as the nuclear potential. That solver builds T_nuc(mu) + diag(v0(grid.points)), calling v0 with EXACTLY grid.points, so a closure returning the curve already tabulated on those same points is an exact lookup rather than an interpolation.

When curves carries eigenvectors, the product phi_j(r; R) chi_v(R) is formed and c-product-normalized (the bilinear, non-conjugated ECS pairing – a conjugated norm would weight the exponentially growing ECS tail instead of cancelling it). Without eigenvectors only the level table is produced, which is the cheap path when no overlap is going to be taken.

Parameters:

curves (ElectronicCurves) – From electronic_curves or resonance_curve. Its energies must be tabulated on g_R.
g_R (FemDvrEcsGrid) – The nuclear grid the curves live on.
mu (float) – Nuclear reduced mass (atomic units) – model.mu.
n_vib (int) – Vibrational levels REQUESTED per curve; see allow_partial.
allow_partial (bool, optional) – Ask each curve for as many of n_vib as it supports, padding the row with NaN, instead of raising when one cannot supply all of them.

Return type:

BoBasis

Raises:

GridError – If curves.energies is not tabulated on g_R, or n_vib < 1.

qscat.core.bo_basis_from_levels(curves, level_energies, level_states, *, curve=0)[source]¶

Product states from an ALREADY-SOLVED nuclear problem – the neutral path.

For a resonance curve the nuclear levels are quasi-bound, and qscat.core.lcp.resonance_levels already solves for them properly: it diagonalizes T(mu) + diag(V_d - i Gamma/2) on two nuclear grids and keeps what is angle-stable, which is a strictly harder problem than vibrational_states’ bound-state solve. Rather than reimplement it, this takes its output and forms the products:

levels = resonance_levels(model, nu_a, nu_b, el_a, el_b, ...)
cur = resonance_curve(model, el_a, el_b, nu_a, seed_window)
basis = bo_basis_from_levels(cur, levels.energies, levels.states)

Parameters:

curves (ElectronicCurves) – Must carry states. Only index curve is used.
level_energies (array_like of complex) – (n_levels,). E_v - i Gamma_v / 2.
level_states (ndarray of complex) – (n_levels, n_R). Nuclear eigenvectors, one row per level, on the SAME nuclear grid the curve is tabulated on. ResonanceLevels.states has exactly this shape and normalization.
curve (int, optional) – Which electronic curve to multiply in (default 0, the only one a resonance_curve family has).

Returns:

Keyed (curve, v), with energies a (1, n_levels) row.

Return type:

BoBasis

Raises:

GridError – On a curve/level grid mismatch, or a missing curve index or states.

qscat.core.n_eff(e_tot, thresholds)[source]¶

Hydrogenic effective quantum number of a state at e_tot.

n_eff = 1/sqrt(2 * binding), with binding measured to the nearest threshold ABOVE e_tot – NOT to the lowest one. The per-level threshold is the physically relevant one: a Rydberg state is bound against the ion level it sits below, not against the incident channel.

Raises ValueError if e_tot lies above every threshold, where the quantity is undefined rather than large.

Parameters:

e_tot (float)
thresholds (ArrayLike)

Return type:

float

qscat.core.admissible_levels(e_tot, thresholds, *, n_eff_max=None)[source]¶

The (curve, vib) levels that can exist AT this energy, from energy alone.

A Rydberg series is attached to a CLOSED channel: above eps[v] that vibrational channel is open and states attached to it are continuum rather than bound. Only thresholds above e_tot therefore contribute, and each contributes exactly one index:

binding = eps[v] - e_tot
n_eff   = 1/sqrt(2 * binding)
curve   ~ n_eff - 1

(the last from the measured series, where Ry_j has n_eff ~ j+1).

The consequence is a strong constraint and the reason this function exists: at fixed energy a HIGHER vibrational level needs a LARGER binding and so a LOWER Rydberg index. The admissible set is finite and small, so a basis can be checked for covering it – which turns “is this pole spurious, or is its partner merely missing?” from a judgement call into a computation. See basis_covers, and qscat.core.assignment.pair_by_overlap which uses it.

The set is finite only away from an accumulation region: as e_tot -> eps[v] the binding tends to zero and the admissible index diverges. That happens in the last ~1 mHa below each threshold. Pass n_eff_max to cut the series there explicitly rather than returning an index no basis will hold.

Parameters:

e_tot (float)
thresholds (ArrayLike)
n_eff_max (float | None)

Return type:

list[tuple[int, int]]

qscat.core.basis_covers(e_tot, thresholds, basis, *, curve_tol=1, n_eff_max=None)[source]¶

Does basis contain every level energetically admissible at e_tot?

When it does, a low overlap means the state has no BO partner at all – spurious. When it does not, a low overlap is uninformative. Without this distinction the two are indistinguishable, and conflating them once nearly cost eight genuine Ry_12..Ry_16 states on H2+.

curve_tol accepts a curve index within +/-1 of the predicted one because the n_eff ~ j+1 mapping is the asymptotic hydrogenic relation and the low curves depart from it.

Parameters:

e_tot (float)
thresholds (ArrayLike)
basis (BoBasis)
curve_tol (int)
n_eff_max (float | None)

Return type:

bool

class qscat.core.ElectronicCurves(energies, states)[source]¶

Electronic eigen-curves tabulated on a nuclear grid.

energies: (n_curves, n_R) complex. energies[j, k] is curve j at nuclear point g_R.points[k]. Complex on the nuclear ECS tail even for a genuinely bound curve – feeding a complex R into model.surface makes the whole electronic problem complex, which is a curve-parametrization artifact of the nuclear rotation and NOT a resonance width.
states: (n_curves, n_r, n_R) complex, phase-aligned across R, or an empty (0, 0, 0) array when built with with_states=False. Each column carries qscat.dvr.eigen’s Hermitian normalization; bo_basis c-normalizes the finished product.

Parameters:

energies (NDArray[complex128])
states (NDArray[complex128])

property n_curves: int¶: Number of electronic curves (rows of energies).

property has_states: bool¶: True when the curves were built with_states=True.

class qscat.core.BoState(psi, energy, curve, vib)[source]¶

One Born-Oppenheimer product state, with the identity that labels it.

Parameters:

psi (NDArray[complex128])
energy (complex)
curve (int)
vib (int)

property label: str¶: omega^j_v as a LaTeX string, the published labelling convention.

class qscat.core.BoBasis(energies, states)[source]¶

A (curve, vib)-keyed family of BO product states plus its level table.

energies is the rectangular level table (NaN where a curve supports fewer than n_vib levels); states is empty when the basis was built from curves carrying no eigenvectors, which is the level-table-only path.

Parameters:

energies (NDArray[complex128])
states (dict[tuple[int, int], BoState])

items()[source]¶

All ((curve, v), BoState) pairs, dict-style.

Return type:: list[tuple[tuple[int, int], BoState]]

property has_states: bool¶: True when the basis holds at least one state.

levels()[source]¶

Every (curve, vib) whose level energy is finite, curve-major.

Return type:: list[tuple[int, int]]

flat()[source]¶

(energy, key) for the finite levels, ascending in Re E.

Return type:: tuple[NDArray[complex128], list[tuple[int, int]]]

Assignment and verification¶

Is a pole a resonance at all, and which quasi-bound state is it? Angle stability is necessary and not sufficient; overlap against a BO basis is the check that catches a rotated-continuum eigenvalue sitting in a stable corner.

qscat.core.overlap(a, b)[source]¶

|<a|b>| under the c-product, with both sides c-normalized.

The c-product (bilinear, NOT conjugated) is the ECS-correct pairing – the same convention qscat.core’s cross sections use. A conjugated dot would weight the exponentially growing ECS tail instead of cancelling it, and on a rotated grid that is not a small difference.

This can exceed 1, and that is not a bug. The c-product is a bilinear form, so Cauchy-Schwarz does not bound it. With both states c-normalized the value is exactly |c(a, b)|, and the inflation over the Hermitian intuition is 1/sqrt(rho_a rho_b), where:

rho = |c(psi, psi)| / ||psi||^2

measures how close to REAL-VALUED (up to one global phase) a state is: 1 for a real vector, falling as the state acquires an internal phase profile. Measured on N2’s broad anion resonances the six clean identifications score 1.02 to 1.19, rising monotonically with Gamma as rho falls 0.66 to 0.42; H2+’s narrow Rydberg resonances stay at 0.87-0.99 because they are nearly real.

rho is NOT real_weight, and substituting one for the other is wrong. N2’s poles have rho of 0.42-0.66 while sitting 96 % inside the unscaled region – well localized, and genuinely complex where they live, because a broad resonance is. H2+’s high Rydberg poles have both collapse together, for the different reason that their orbitals leave the box. rho explains the inflation; real_weight catches the escape.

Do NOT “fix” this by dividing by the Euclidean norms. |c(a,b)| / (||a|| ||b||) is bounded by 1 and is WRONG here, which was measured rather than argued. Its denominator weights the exponentially growing ECS tail – reintroducing exactly the contamination the c-product’s numerator cancels, which is the same error as using vdot. On H2+ window 0 it collapses to 0.03, 0.008 and 0.006 for three states whose node counts identify them unambiguously as Ry_14..Ry_16 v=1, and it re-ranks all three onto the wrong (compact) partner. It penalizes diffuse states for being diffuse.

The right response to a state whose real_weight is small is not a different overlap: it is to report the localization and treat the state as box-limited. pair_by_overlap does that.

Returns 0.0 when either vector is (numerically) self-orthogonal, which is the honest answer: the normalization it would need does not exist.

Parameters:

a (NDArray[complex128])
b (NDArray[complex128])

Return type:

float

qscat.core.real_weight(psi, tgrid)[source]¶

Fraction of |psi|^2 inside the UNSCALED region of tgrid.

How much of a state the grid actually holds, as opposed to pushes into the rotated tail. qscat.core.lcp.ResonanceLevels reports exactly this for 1-D levels; the 2-D path had no equivalent, and the gap hid a real defect.

Distinct from the rho = |c(psi,psi)| / ||psi||^2 that inflates overlap above 1 (see there): rho measures how close to real-valued a state is, this measures where it lives. N2’s poles sit at real_weight 0.96 with rho 0.42-0.66 – localized, and genuinely complex where they live. The two move together only when the cause is a tail.

The overlap cannot see this and is not supposed to. The c-product cancels the rotated tail by construction, so a state that is 97 % tail still pairs at 0.99 with the BO product it genuinely is – a correct answer to the question asked, and a badly misleading summary of the state. Measured on H2+ window 0, real_weight falls 0.998 -> 0.68 -> 0.29 -> 0.12 -> 0.031 -> 0.008 -> 0.0003 as the Rydberg series climbs toward its threshold: those orbitals are simply larger than the 300-bohr box, and their overlaps stay at 0.99 throughout.

Returns a value in [0, 1]. Near 1 the grid holds the state; small values mean the answer is about the box, not the physics.

Parameters:

psi (NDArray[complex128])
tgrid (TensorGrid)

Return type:

float

qscat.core.Verdict¶: alias of Literal[‘ok’, ‘spurious’, ‘basis-limited’, ‘box-limited’, ‘weak’, ‘mixed’, ‘distant’]

qscat.core.pair_by_overlap(pole_energy, pole_state, basis, thresholds=None, *, no_partner=0.1, weak=0.5, mixed_ratio=0.7, max_shift_ev=0.02, n_eff_max=None, basis_complete=None, localization=None, min_localization=0.5)[source]¶

Pair one exact pole to the BO level it most resembles, with checks.

Parameters:

pole_energy (complex) – The pole, E_r - i Gamma/2. Only Re E enters the shift.
pole_state (ndarray of complex) – Its eigenvector, flat over the same (r, R) layout as basis.
basis (BoBasis) – Reference states, from qscat.core.bo.bo_basis. Must carry states.
thresholds (array_like of float, optional) – Channel thresholds for the closed-channel admissibility check. Omit at your peril: without it a poor overlap cannot be attributed to the state rather than to the basis, and every such pole is reported basis-limited rather than spurious.
no_partner (float, optional) – Verdict thresholds; see the module docstring on calibration.
weak (float, optional) – Verdict thresholds; see the module docstring on calibration.
mixed_ratio (float, optional) – Verdict thresholds; see the module docstring on calibration.
max_shift_ev (float, optional) – Verdict thresholds; see the module docstring on calibration.
n_eff_max (float, optional) – Passed to the admissibility check to cut the Rydberg series short of an accumulation region.
basis_complete (bool, optional) – Assert directly that the basis holds every state the BO picture admits here, overriding the thresholds computation. This exists for the NEUTRAL case, where there is no Rydberg series to count: an anion resonance curve has one electronic state and the basis is complete once all its vibrational levels are built. It is a claim, not a check – asserting it wrongly converts every basis-limited verdict into spurious, which is the one verdict that rejects a pole.
localization (float, optional) – This pole’s real_weight. Supply it: the overlap is blind to a state that has left the box, by design, so without this a box-limited pole is reported as a clean identification.
min_localization (float, optional) – Below this real_weight, the verdict is box-limited.

Return type:

OverlapPair

Raises:

ValueError – If basis carries no states (built with with_states=False).

qscat.core.pair_one_to_one(pole_energy, level_energy, *, max_distance=0.001)[source]¶

Match poles to levels ONE-TO-ONE by POSITION, minimising total displacement.

Only real parts are compared; widths play no part in the matching.

Returns {pole index: level index}, omitting any pair further apart than max_distance (default 1 mHa).

A cross-check, not the assignment. pair_by_overlap is the primary criterion; where the two disagree, overlap wins. This is here because an independent pairing that agrees is evidence, and one that disagrees marks a row worth looking at by hand.

Why not nearest-neighbour. Assigning each pole to whichever level is closest treats the poles independently, so two poles can claim one level while another goes unclaimed – and it fails worst exactly where the physics is most interesting. On H2+ two BO levels 20 uHa apart correspond to poles 154 uHa apart; nearest-neighbour hands both poles to one level and calls the second pairing ambiguous, when it is fully determined once the matching is required to be a bijection. A minimum-total-cost assignment uses the constraint the physics supplies – distinct states pair with distinct levels.

It is not infallible: where a whole neighbourhood is denser than the shifts, the global optimum can still permute members within it, so a pairing is only as trustworthy as the local spacing. Report the per-pair distance so that stays visible.

Parameters:

pole_energy (ArrayLike)
level_energy (ArrayLike)
max_distance (float)

Return type:

dict[int, int]

qscat.core.peak_positions(energy, sigma, *, prominence=10.0)[source]¶

Energies of prominent local maxima in a cross-section sweep.

A local maximum counts when it exceeds both neighbours AND is prominence times the smaller of the two samples two steps out – a deliberately crude test, because the sweeps this is used on are resonance spectra where a genuine peak towers over its surroundings and a sophisticated peak finder would mostly add tunable knobs.

Returns an empty array for a sweep shorter than five samples.

Parameters:

energy (NDArray[float64])
sigma (NDArray[float64])
prominence (float)

Return type:

NDArray[float64]

qscat.core.peak_alignment(marks, peaks, *, width)[source]¶

Distance from each computed position to the nearest observed peak.

Distances are quoted in units of a resonance width, and that is the point. The energy axis is not the scale on which “lands on the peak” means anything – a 3 uHa disagreement is negligible against a 1 meV-wide feature and fatal against a 0.02 meV-wide one. width is the FWHM of a typical resonance in the set being compared.

Marks with no peak on the same side are not special-cased: every mark takes its distance to the nearest peak anywhere in peaks, so a mark far outside the peak range reports a large distance rather than being dropped.

Raises ValueError if width <= 0, or if either array is empty – an alignment of nothing against something is not a number worth returning.

Parameters:

marks (ArrayLike)
peaks (ArrayLike)
width (float)

Return type:

PeakAlignment

class qscat.core.OverlapPair(pole_energy, level, overlap, second_level, second_overlap, shift_ev, verdict)[source]¶

One exact pole paired to a BO level BY OVERLAP, with the checks.

Parameters:

pole_energy (float)
level (tuple[int, int] | None)
overlap (float)
second_level (tuple[int, int] | None)
second_overlap (float)
shift_ev (float)
verdict (Literal['ok', 'spurious', 'basis-limited', 'box-limited', 'weak', 'mixed', 'distant'])

property shift_mev: float¶: shift_ev in milli-eV – the unit the reference tables quote.

property label: str¶: omega^j_v, or - when no level was assigned.

property is_quotable: bool¶: True only for ok. Every other verdict withholds something.

property is_resonance: bool¶

False only for spurious – the one verdict that rejects a pole.

basis-limited and weak withhold a label without rejecting the state, which is the honest reading: the evidence is about the basis, not about the pole.

class qscat.core.PeakAlignment(n_marks, n_peaks, median_widths, within_one_width, distances)[source]¶

How well a set of computed positions lands on a set of observed peaks.

Parameters:

n_marks (int)
n_peaks (int)
median_widths (float)
within_one_width (int)
distances (NDArray[float64])

Nonlocal resonance model¶

qscat.core.nrm – the Feshbach-projection nonlocal resonance model (NRM) for dissociative attachment (Houfek, Rescigno & McCurdy, Phys. Rev. A 77, 012710 (2008)). It sits between the LCP approximation above and the exact 2-D solvers, and is kept as a subpackage that qscat.core does not import at runtime (see the hard boundary above) – import it explicitly, e.g. from qscat.core.nrm import nrm_da_cross_section. See docs/physics/nonlocal-resonance-model.md.

class qscat.core.nrm.DiscreteState(*args, **kwargs)[source]¶

Anything that supplies phi_d(r;R) as c-normalized DVR coefficients.

phi_d(R)[source]¶

The discrete state at nuclear coordinate R.

Parameters:: R (float)
Return type:: NDArray[complex128]

class qscat.core.nrm.AsymptoticDiscreteState(grid, model, R_inf)[source]¶

Choice B (Sec. VI B): the R-independent bound state phi_b.

One electronic eigenproblem for the whole calculation. Eq. (67) holds trivially. The paper finds this choice essentially exact for DA and for VE provided the background T-matrix terms are included (p. 012710-9-10).

Parameters:

grid (FemDvrEcsGrid)
model (ResonanceModel)
R_inf (float)

phi_d(R)[source]¶

The same state at every R, by construction.

Parameters:: R (float)
Return type:: NDArray[complex128]

class qscat.core.nrm.PhysicalDiscreteState(grid, model, R_descending, elec_grid_b, *, r_d=10.0, re_half_width=0.08, im_half_width=0.08, resid_tol=0.001, seed_window=None)[source]¶

Choice A (Sec. VI A): the scattering function at Re E_res(R).

At each R the resonance energy is taken from the two-angle ECS pole walk (qscat.core.lcp.resonance_pole_walk, which returns the ELECTRONIC shift s(R) = Re E_pole(R) - v0(R)). Where s(R) > 0 the state is the fixed-nuclei scattering function at that energy (Eq. 17’s H_el); where s(R) <= 0 the electron is bound and the state is the lowest GENUINELY bound eigenvector of H_el(R) instead (p. 012710-8): Re(E) < 0 and |Im(E)| < _BOUND_IM_TOL. If no such eigenvalue exists, the walk’s sign decision does not match the fresh spectrum at this R (most likely a frozen/stale continuation point from a R_descending too coarse for resonance_pole_walk to track – see re_half_width/im_half_width below) and ConvergenceError is raised rather than silently returning whatever eigenvector happens to have the smallest real part – that eigenvector is not necessarily real or bound. Only the SCATTERING branch is then truncated by Eq. (69) (the paper’s cutoff applies to the scattering-derived state only, p. 012710-8); the bound branch is already square-integrable and is used as-is, so its exact eigenrelation H_el phi_b = E_b phi_b survives – truncating it would spuriously prevent V_dk -> 0 at large R (Eq. 67). Both branches are then c-normalized, which supplies the paper’s exp[-i delta(R)].

R_descending must be descending – the pole walk is seeded at large R and continued inward. States are precomputed on those nodes; phi_d(R) returns the nearest one, which is how the ingredient builder consumes it (one call per nuclear DVR point, on the same node set). The per-node electronic shift and which branch was taken are recorded as shift, gamma, and used_scattering (all aligned with R_descending), so a caller can inspect what happened rather than infer it from the output state alone.

The pole walk’s seed window is centered on the bound ANION electronic state at the outermost (largest) R in R_descending (qscat.core.dissociation.anion_electronic_states), mirroring qscat.core.lcp.curve._walk_from_anion_seed: at that R the discrete state is expected to already be essentially the asymptotic bound state (Eq. 67), so its energy is a reliable seed for the pole finder. seed_window overrides this and is passed straight to resonance_pole_walk if given.

re_half_width/im_half_width default wider (0.08 Ha) than qscat.core.lcp’s own default (0.05 Ha, tuned for a FINE nuclear grid’s closely-spaced real nodes): R_descending here is an arbitrary, caller-supplied array and may be sampled far more coarsely (e.g. eMoScat’s hand-picked nuclear decks step by 0.4 bohr near F2’s bound/resonant crossing around R ~ 2.6-3.0 bohr, where the electronic shift moves by ~0.09 Ha over that step) – too coarse for the 0.05 Ha window to recenter onto without losing the pole, which freezes the walk exactly where the physics is most interesting (the crossing itself).

Parameters:

grid (FemDvrEcsGrid)
model (ResonanceModel)
R_descending (npt.NDArray[np.float64])
elec_grid_b (FemDvrEcsGrid)
r_d (float)
re_half_width (float)
im_half_width (float)
resid_tol (float)
seed_window (tuple[float, float, float, float] | None)

phi_d(R)[source]¶

The discrete state at the precomputed node nearest R.

Parameters:: R (float)
Return type:: NDArray[complex128]

qscat.core.nrm.nrm_ingredients(elec_grid, model, phi_d, R_values)[source]¶

Build E_n(R), V_dn(R) and V_d(R) on the given nuclear nodes.

Parameters:

elec_grid (FemDvrEcsGrid) – The electronic radial grid.
model (ResonanceModel) – Supplies surface, v0 and ell.
phi_d (DiscreteState) – The discrete state (choice A or B).
R_values (ndarray) – Nuclear nodes, strictly DESCENDING. The tracking walk starts at the largest R, where phi_d -> phi_b and the P-space states are cleanly ordered, and continues inward.

Return type:

NrmIngredients

Raises:

ValueError – If R_values is not strictly descending.
ConvergenceError – If the P H_el P null mode cannot be identified unambiguously, or an eigenvector is c-product self-orthogonal.

Warns:

UserWarning – If the minimum _sign_align overlap across the R walk falls below _MIN_OVERLAP_WARN – a likely tracking failure (see NrmIngredients.min_overlap).

class qscat.core.nrm.NrmIngredients(R, v_d_discrete, E_n, V_dn, min_overlap=1.0)[source]¶

The energy-independent inputs to the nonlocal potential.

Parameters:

R (NDArray[float64])
v_d_discrete (NDArray[complex128])
E_n (NDArray[complex128])
V_dn (NDArray[complex128])
min_overlap (float)

R¶

The nuclear nodes, descending, as supplied to nrm_ingredients.

Type:: ndarray

v_d_discrete¶

V_d(R) of Eq. (20), complex, INCLUDING V_0(R). This is the paper’s discrete-state potential, NOT qscat.core.lcp’s Vd.

Type:: ndarray

E_n¶

(n_R, n_states) projected electronic energies of Eq. (56), EXCLUDING V_0(R) (which Eq. 61 adds separately). Adiabatically tracked across R; the null mode is dropped.

Type:: ndarray

V_dn¶

(n_R, n_states) discrete-continuum couplings of Eq. (59), aligned with E_n.

Type:: ndarray

min_overlap¶

Diagnostic: the minimum |c_product(prev, cur)| seen by _sign_align across the whole R walk, over every tracked state. Near 1.0 means _track paired the same physical P-space state at every adjacent (R_j, R_{j+1}); a value near 0 means it did not – a tracking failure, silently sign-“corrected” rather than raised (see nrm_ingredients, which warns below _MIN_OVERLAP_WARN). Defaults to 1.0 for an NrmIngredients built directly (e.g. in tests) rather than through nrm_ingredients, where no walk was performed.

Type:: float

qscat.core.nrm.nrm_da_cross_section(nuclear_grid, elec_grid, model, phi_d, eps, chi, v_init, E, *, ingredients=None, n_states=None)[source]¶

sigma_DA(E) in the nonlocal resonance model (bohr^2).

Parameters:

nuclear_grid (FemDvrEcsGrid) – The nuclear and electronic radial grids.
elec_grid (FemDvrEcsGrid) – The nuclear and electronic radial grids.
model (ResonanceModel) – The molecule.
phi_d (DiscreteState) – The discrete-state choice under test.
eps (ndarray) – Neutral vibrational energies and states (qscat.core.vibrational).
chi (ndarray) – Neutral vibrational energies and states (qscat.core.vibrational).
v_init (int) – Initial vibrational level.
E (float or array) – Incident electron kinetic energy or energies (hartree).
ingredients (NrmIngredients, optional) – Precomputed ingredients; built here if omitted. Pass them in when sweeping energies or comparing discrete-state choices – they are energy-independent and dominate the cost.
n_states (int, optional) – Truncate the sum over projected electronic states. None uses all.

Returns:

sigma_DA per energy; scalar-shaped for a scalar E.

Return type:

ndarray

qscat.core.nrm.nrm_ve_cross_section(nuclear_grid, elec_grid, model, phi_d, eps, chi, v_init, vprimes, E, *, ingredients=None, n_states=None, include_background=True)[source]¶

sigma_{v_init->v’}(E) in the nonlocal resonance model (bohr^2).

sigma = 4 pi^3 |T^res + T^bg|^2 / k_i^2, Eq. (28)/(31)/(37), on qscat.core.driven’s own normalization so the exact 2-D and nonlocal curves compare directly rather than through two conventions that happen to agree.

T^res (Eq. 31) is a contraction of the SAME Psi_d^+ solution dissociation.nrm_da_cross_section computes (_psi_d_for_energy, built once per energy and reused across every vprimes entry). T^bg (Eq. 37) is the non-resonant background PRA 77 shows a bare LCP curve is missing – largest for the broadest resonance.

include_background=False drops T^bg, giving the paper’s “nonlocal” curve as against its “nonlocal + bg” one (Figs. 4 and 8 plot both). That difference is the method’s own argument, not a debugging switch.

Parameters:

nuclear_grid (FemDvrEcsGrid) – The nuclear and electronic radial grids.
elec_grid (FemDvrEcsGrid) – The nuclear and electronic radial grids.
model (ResonanceModel) – The molecule.
phi_d (DiscreteState) – The discrete-state choice under test.
eps (ndarray) – Neutral vibrational energies and states (qscat.core.vibrational).
chi (ndarray) – Neutral vibrational energies and states (qscat.core.vibrational).
v_init (int) – Initial vibrational level.
vprimes (list of int) – Final vibrational levels.
E (float or array) – Incident electron kinetic energy or energies (hartree).
ingredients (NrmIngredients, optional) – Precomputed ingredients; built here if omitted. Pass them in when sweeping energies or comparing discrete-state choices – they are energy-independent and dominate the cost.
n_states (int, optional) – Truncate the sum over projected electronic states in F(E). None uses all.
include_background (bool, default True) – Add T^bg (Eq. 37) to T^res before squaring. True is the paper’s “nonlocal + background” curve; False is its bare “nonlocal” curve.

Returns:

sigma_{v_init->v’} per energy; scalar E returns shape (len(vprimes),), array E returns (len(E), len(vprimes)) – driven.ve_cross_section’s convention. A closed channel (E_tot - eps_vf <= 0, or E <= 0) contributes 0.0.

Return type:

ndarray

Plotting¶

qscat.core.plot_cross_sections(E_grid, sigma, *, channels=None, reference=None, thresholds=None, title=None, path)[source]¶

Plot one sigma(E) curve per channel, save to path as a PNG.

E_grid is (M,) Hartree; sigma is (M, C) bohr^2, column c is channel channels[c] (default range(C) when channels is None).

reference, if given, is (E_ref, sigma_ref) with sigma_ref shaped (N, C) – the SAME column order/count as sigma – drawn as dashed, marked curves in the same color as their matching computed channel, so the eye pairs “solid vs. dashed of the same color” rather than reading a legend. This is a plain overlay: this function does no interpolation, resampling, or alignment between E_grid and E_ref – the two curves are simply drawn on the same axes.

thresholds, if given, is {channel: threshold_energy_Ha}; channels present in both channels and this dict get a dotted vertical line at their threshold, colored to match that channel’s curve.

Non-positive cross sections (e.g. a closed channel returning exactly sigma=0 below its own threshold) are masked to NaN before plotting on the log-scaled y-axis: matplotlib’s log locator silently drops them from the visible curve rather than raising, but leaving raw zeros in would still emit a RuntimeWarning at log10(0) – this masks that cleanly rather than suppressing the warning.

No physics lives here: this function does not know what a channel, a threshold, or a molecule is – it plots whatever arrays it is given.

Requires the optional qscat[plot] extra (matplotlib); raises ModuleNotFoundError with an actionable hint if it is not installed.

Parameters:

E_grid (NDArray[float64])
sigma (NDArray[float64])
channels (list[int] | None)
reference (tuple[NDArray[float64], NDArray[float64]] | None)
thresholds (dict[int, float] | None)
title (str | None)
path (str | PathLike[str])

Return type:

None

qscat.core.plot_resonance_levels(levels, *, path, curves=None, band=None, xlim=None, ylim=None, baseline=None, pairing=None, title=None, xlabel='R (bohr)', ylabel='energy (Hartree)')[source]¶

Plot complex resonance levels omega_i = E_r - i*Gamma/2, save to path.

The left panel is the level diagram in the conventional form: any curves are drawn as solid lines against the internuclear coordinate, band shades a width envelope around one of them, and every level in levels becomes a horizontal dashed line at its Re E, one colour per series. The right panel plots each series’ difference from baseline, in meV, against level index.

Parameters:

levels (dict of str to ndarray of complex) – One entry per series, e.g. {“exact 2-D”: …, “BO / LCP”: …}.
curves (dict of str to (ndarray, ndarray), optional) – Background potential curves as (x, y) pairs, e.g. the neutral V0(R) and the resonance curve E_res(R).
band (tuple of (ndarray, ndarray, ndarray), optional) – (x, centre, half_width) shaded envelope — the Gamma(R)/2 band around a resonance curve.
baseline (str, optional) – Series the others are differenced against. Omit for a single panel.
pairing (dict of str to sequence of (int, int), optional) – Explicit (series index, baseline index) pairs per non-baseline series. Supply this whenever the difference panel is meant to be read as a physical shift. Without it the series are paired level-by-level in ascending Re E and truncated to the shortest, which is only correct when both sets are complete and ordered alike – and that assumption fails exactly where the physics is interesting. Measured on H2+: two BO levels 20 uHa apart correspond to exact poles 154 uHa apart, so index pairing crosses them and reports two shifts that belong to neither level. qscat.core.assignment.pair_by_overlap produces a pairing that is defensible; pair_one_to_one produces one that is at least a bijection. When pairing is given the series are used in the order supplied, NOT sorted.
path (str | PathLike[str])
xlim (tuple[float, float] | None)
ylim (tuple[float, float] | None)
title (str | None)
xlabel (str)
ylabel (str)

Return type:

None

Notes

No physics: it takes arrays and labels, so it serves any pair of level sets on any coordinate.

class qscat.core.ComparisonPanel(series, title=None, ylim=None, yticks=None)[source]¶

One panel of a plot_route_comparison grid.

series maps a legend label to a (M,) array on the shared E_grid; every entry is drawn on the same axes in the order given. ylim/yticks exist so a caller replicating a published panel can pin its axis rather than letting matplotlib choose one – the whole point of such a figure is that a reader can lay it next to the printed panel.

No physics: this is a label, some arrays, and axis limits.

Parameters:

series (dict[str, NDArray[float64]])
title (str | None)
ylim (tuple[float, float] | None)
yticks (NDArray[float64] | None)

qscat.core.plot_route_comparison(E_grid, panels, *, styles=None, path, xlim=None, xlabel='E (Hartree)', ylabel='$\\sigma$ (bohr$^2$)', suptitle=None, panel_size=(6.0, 4.2))[source]¶

Plot a grid of panels, each overlaying several named sigma(E) routes.

plot_cross_sections’ sibling for the other comparison shape: that one puts one CHANNEL per curve on a single log-scaled axis; this one puts one METHOD per curve on a LINEAR axis, and takes a grid of panels so several transitions (or several model variants) can be read at once.

Parameters:

E_grid (ndarray) – (M,) shared abscissa; every series in every panel is drawn on it.
panels (sequence of sequence of ComparisonPanel) – Row-major grid, panels[row][col]. Rows need not be equal length; the figure is sized to the longest, and missing cells are left blank.
styles (mapping, optional) – {series label: matplotlib kwargs}, applied to that label wherever it appears. Labels with no entry fall back to matplotlib’s own cycle. This is how a caller reproduces a published figure’s line styles.
xlim (tuple[float, float] | None) – Axis framing. panel_size is per panel, in inches.
xlabel (str) – Axis framing. panel_size is per panel, in inches.
ylabel (str) – Axis framing. panel_size is per panel, in inches.
suptitle (str | None) – Axis framing. panel_size is per panel, in inches.
panel_size (tuple[float, float]) – Axis framing. panel_size is per panel, in inches.
path (str | PathLike[str])

Return type:

None

Notes

Linear y by design: an absolute difference on a linear axis is the scale on which a published “the curves are practically the same” claim is made, and a log axis would misrepresent it. Non-positive values are plotted as they are (no masking) – on a linear axis a zero is meaningful.

Requires the optional qscat[plot] extra (matplotlib); raises ModuleNotFoundError with an actionable hint if it is not installed.