solprop.nonlinear.ionisation#

Strong-field ionisation: the rate protocol, its conventions and diagnostics.

The rate W(|E|) is the probability per unit time that one atom loses its outermost electron in a field of amplitude |E|. It is the input to the plasma response: the ionised fraction is 1 - exp(-int W dt) and the energy the field pays for each electron is the energy per event.

solprop keeps no ionisation-rate formula of its own (SEP 3, D1). Every rate of a built-in gas, and of a custom species from explicit parameters, is one of optmat’s rate models, adapted by solprop.nonlinear.optmat_rates.OptmatRate; a caller may also supply any object satisfying IonisationRate (or the narrower EvaluableRate), or a tabulated solprop.translate.ionisation.RateTable. This module holds the contract they meet, their optical-cycle conventions, the energy binding of the plasma loss (bind_rate_energy()), the rate ceiling MAX_RATE, and the diagnostics keldysh(), barrier_suppression_field() and ionisation_fraction().

The manual page is docs/guide/media/ionisation.md. Rates are per atom or molecule, in s^-1, for a field amplitude in V/m, and are exactly zero below the rate’s threshold field and clipped at its ceiling so that a diverged step cannot put a NaN into the plasma column (docs/dev/stepper.md section 10, trap 9).

Module Attributes

`RateConvention`	Optical-cycle convention of a magnitude rate, distinct from orientation averaging.
`MAX_RATE`	Default ceiling on an ionisation rate, s^-1.
`IntegrationPolicy`	How the plasma column integrates along t (SEP 18).
`DT_RTOL`	Relative tolerance on the spacing a band `RunningIntegral` is used at.

Functions

`barrier_suppression_field`(...)	Return the field at which the Coulomb barrier is suppressed below the level.
`bind_rate_energy`(gas, rate)	Return the species a plasma response takes, its energy the rate's own.
`ionisation_fraction`(rate, field, dt, *[, ...])	Return the ionised fraction along a time axis for a given field.
`keldysh`(gas, wavelength, field_amplitude)	Return the Keldysh parameter of a gas at a wavelength and field amplitude.
`rate_convention`(rate)	Read an explicit optical-cycle convention without guessing for legacy rates.
`rate_energy`(rate)	Read the energy per ionisation event a rate declares for the plasma loss.
`rate_species`(rate)	Read a rate's declared species name and binding energy.

Classes

`EvaluableRate`(args, *kwargs)	Anything that turns a field amplitude into an ionisation rate.
`IonisationRate`(args, *kwargs)	A field-dependent single-ionisation rate per atom or molecule, in s^-1.
`RunningIntegral`(policy, taps[, dt, ...])	The resolved running-integral rule for field-grid plasma (SEP 18).

solprop.nonlinear.ionisation.DT_RTOL: Final[float] = 1e-12#: Relative tolerance on the spacing a band RunningIntegral is used at.

class solprop.nonlinear.ionisation.EvaluableRate(*args, **kwargs)[source]#

Bases: Protocol

Anything that turns a field amplitude into an ionisation rate.

The narrow contract a consumer of a rate actually needs: call it with an amplitude in V/m and get a rate per atom in s^-1. Both an IonisationRate (optmat’s OptmatRate, or a caller’s) and the tabulated solprop.translate.ionisation.RateTable satisfy it, structurally and without inheriting from it.

It exists because the plasma response (solprop.nonlinear.responses.PlasmaResponse) accepts either, and a response is in the physics layer while a table is in the translation layer: imports run P to T and never back (docs/dev/architecture.md section 3), so the response cannot name RateTable even in an annotation, because the test session’s beartype hook checks annotations at run time.

Examples

>>> from solprop.nonlinear.ionisation import EvaluableRate
>>> from solprop.nonlinear.optmat_rates import OptmatRate
>>> isinstance(OptmatRate("Ar", 800e-9), EvaluableRate)
True

solprop.nonlinear.ionisation.IntegrationPolicy#

How the plasma column integrates along t (SEP 18).

"band" filters each integrand with the grid’s designed symmetric taps before the trapezoid, so that the running integral is exact to a certified tolerance across the active band; "trapezoid" is the plain cumulative trapezoid, kept permanently for the frozen witness comparisons.

alias of Literal[‘band’, ‘trapezoid’]

class solprop.nonlinear.ionisation.IonisationRate(*args, **kwargs)[source]#

Bases: Protocol

A field-dependent single-ionisation rate per atom or molecule, in s^-1.

Implementations are vectorised in the field amplitude, even in its sign, exactly zero below threshold and never above maximum, so the translation layer can tabulate them and a kernel can evaluate the table without a divergence branch.

property citation: str[source]#: One-line reference for the model and the options in force.

property maximum: float[source]#: Ceiling on the rate, s^-1; the value is clipped there.

property parameters: tuple[tuple[str, float | int | bool | str], ...][source]#

Every setting that changes the rate, name-sorted and hashable.

The translation layer builds its table cache key from this, so two rates that differ in any way that matters give different keys and no ordering of keyword arguments can give two keys for one rate.

property threshold: float[source]#: Field amplitude, V/m, below which the rate is exactly zero.

solprop.nonlinear.ionisation.MAX_RATE: Final[float] = 4.134137333517015e+16#

Default ceiling on an ionisation rate, s^-1.

One ionisation per atomic unit of time (4.13e16 s^-1). No quasistatic rate means anything at or beyond this – the electron would leave in less than an orbit – so clipping here cannot touch a physical result, and it is what stops a diverged step from propagating an overflow into the plasma column.

solprop.nonlinear.ionisation.RateConvention#

Optical-cycle convention of a magnitude rate, distinct from orientation averaging.

alias of Literal[‘instantaneous-isotropic’, ‘cycle-averaged’, ‘unspecified’]

class solprop.nonlinear.ionisation.RunningIntegral(policy, taps, dt=None, band_edge=None, certified_error=0.0, boundary='point-reflection', exposure='floor-at-zero')[source]#

Bases: object

The resolved running-integral rule for field-grid plasma (SEP 18).

Each running integral int y dt is the cumulative trapezoid of a filtered integrand,

\[\tilde y_n = c_0 y_n + \sum_{k=1}^{K} c_k (y_{n-k} + y_{n+k}),\]

with symmetric taps c whose response, times the trapezoid’s own multiplier (theta/2) cot(theta/2), is one to within certified_error across the active band. The taps are a property of the grid: the translation layer designs them (solprop.translate.plasma_integration.design_running_integral()), and this value type only carries and applies them, so the physics layer never fits anything itself. K samples beyond each end of the window are supplied by point reflection, y_{-m} = 2 y_0 - y_m, which is exact for constant and linear integrands. The exposure is floored at zero rather than clamped, so the fraction stays in [f_0, 1) and is non-decreasing to within the filter’s undershoot.

trapezoid() is the legacy rule: one identity tap, no floor, and the arithmetic of scipy.integrate.cumulative_trapezoid() unchanged.

Attributes:

policy{‘band’, ‘trapezoid’}: The rule.
tapsnumpy.ndarray: (c_0, c_1, ..., c_K), centre first, read-only float64, shape (K + 1,); (1.0,) for the trapezoid. The filter has 2K + 1 points.
dtfloat or None: The sample spacing the taps were designed for, s; None for the trapezoid, which is the same rule at every spacing.
band_edgefloat or None: The upper edge of the active band the taps were fitted to, rad/s; None for the trapezoid.
certified_errorfloat: The certified max |H - 1| across the band, where H is the composite multiplier relative to the exact integral; zero for the trapezoid, whose error is not designed.
boundarystr: How the filter reads past the window’s ends: 'point-reflection', or 'none' for the trapezoid.
exposurestr: The positivity policy of the exposure: 'floor-at-zero', or 'none' for the trapezoid.

Raises:

GridError: If the taps are not a non-empty one-dimensional finite array, the trapezoid carries anything but the identity tap, or a band rule has no positive finite dt and band_edge.

Examples

>>> import numpy as np
>>> from solprop.nonlinear.ionisation import RunningIntegral
>>> legacy = RunningIntegral.trapezoid()
>>> legacy.integrate(np.ones(4), 0.5).tolist()
[0.0, 0.5, 1.0, 1.5]

band_edge: float | None = None#

boundary: str = 'point-reflection'#

certified_error: float = 0.0#

check(dt)[source]#

Refuse a band rule used at a sample spacing it was not designed for.

Parameters:

dtfloat: The spacing of the samples about to be integrated, s.

Raises:

GridError: If the rule is band and dt differs from its own by more than DT_RTOL relative.

dt: float | None = None#

exposure: str = 'floor-at-zero'#

filter(values)[source]#

Return the filtered integrand, with its ends closed by point reflection.

Parameters:

valuesnumpy.ndarray: The integrand along the last axis, at least K + 1 samples long.

Returns:

numpy.ndarray: c_0 y_n + sum_k c_k (y_{n-k} + y_{n+k}), the shape of values.

Raises:

GridError: If the axis is shorter than K + 1, so the reflection has nothing to reflect.

floor(exposure)[source]#

Apply the exposure policy: the floor at zero, or nothing for the trapezoid.

The floor is pointwise, not a running maximum, so it accumulates no bias: the final exposure is the unfloored one wherever that is positive.

Parameters:

exposurenumpy.ndarray: The running integral of the rate, dimensionless.

Returns:

numpy.ndarray: max(exposure, 0) under the band rule, exposure itself otherwise.

property half_width: int[source]#

Return K, the number of samples the filter reads on each side.

Returns:

int: len(taps) - 1; zero for the trapezoid.

integrate(values, dt)[source]#

Return the running integral of values along their last axis.

The filter’s sum has a fixed order, which the compiled backends reproduce: the centre term first, then the pairs k = 1 ... K, each pair added before it is weighted. The virtual samples are formed as 2 y_0 - y_m. The integral starts at zero, as the trapezoid’s does.

Parameters:

valuesnumpy.ndarray: The integrand, sampled uniformly along the last axis, at least K + 1 samples long.
dtfloat: The sample spacing, s.

Returns:

numpy.ndarray: The running integral, the shape of values.

Raises:

GridError: If dt is not the band rule’s own, or the axis is shorter than K + 1.

Examples

>>> import numpy as np
>>> from solprop.nonlinear.ionisation import RunningIntegral
>>> taps = np.array([1.1, -0.05])
>>> rule = RunningIntegral("band", taps, dt=0.5, band_edge=1.0)
>>> rule.integrate(np.arange(4.0), 0.5).tolist()
[0.0, 0.25, 1.0, 2.25]

policy: Literal['band', 'trapezoid']#

taps: Float[ndarray, 'n_taps']#

classmethod trapezoid()[source]#

Return the legacy rule: the plain cumulative trapezoid.

Returns:

RunningIntegral: One identity tap, no boundary rule and no exposure floor.

Examples

>>> from solprop.nonlinear.ionisation import RunningIntegral
>>> RunningIntegral.trapezoid().taps.tolist()
[1.0]

solprop.nonlinear.ionisation.barrier_suppression_field(ionisation_potential_j, charge)[source]#

Return the field at which the Coulomb barrier is suppressed below the level.

Classically the electron is free once the field has pulled the barrier down to the bound-state energy, which happens at E_b = Ip^2 / (4 Z) in atomic units. It is the natural scale for the barrier-suppression corrections and for the upper end of a rate table.

Parameters:

ionisation_potential_jfloat: Ionisation potential in joules.
chargefloat: Charge seen by the departing electron, in units of the elementary charge.

Returns:

float: Barrier-suppression field amplitude in V/m.

Examples

>>> from solprop.nonlinear.ionisation import barrier_suppression_field
>>> from solprop.physdata import Gas, ionisation_potential
>>> f"{barrier_suppression_field(ionisation_potential(Gas.He), 1.0):.3e}"
'1.050e+11'

solprop.nonlinear.ionisation.bind_rate_energy(gas, rate)[source]#

Return the species a plasma response takes, its energy the rate’s own.

For a rate that declares its energy per event (rate_energy()), the species with its ionisation potential replaced by that energy: the plasma loss then pays the energy the rate’s exponent uses, and a response’s potential, its compiled channel and the envelope’s consistency check all agree with the rate. A rate that declares none leaves the species exactly as given.

Parameters:

gasGas, str or Species: The response’s species.
rateEvaluableRate: The response’s rate.

Returns:

Gas, str or Species: The species as given, or a Species whose ionisation potential is the rate’s energy in J.

Raises:

MaterialError: If the rate declares an energy for another species.

Examples

>>> from solprop.nonlinear.ionisation import bind_rate_energy
>>> def rate(field):
...     return 0.0 * field
>>> bind_rate_energy("Ar", rate)
'Ar'

solprop.nonlinear.ionisation.ionisation_fraction(rate, field, dt, *, preionisation=0.0, integral=None)[source]#

Return the ionised fraction along a time axis for a given field.

Each atom or molecule survives with probability exp(-int W dt). The fraction is f = 1 - (1 - f0) exp(-int W dt) where f0 is the fraction already ionised before the pulse. Carrying f0 in the survival factor keeps f at or below one.

The integral is the plasma column’s own: pass the problem’s RunningIntegral to get the fraction a propagation forms under the band rule (SEP 18). Without one this is the plain cumulative trapezoid.

Parameters:

rateEvaluableRate: The rate to integrate, s^-1 as a function of field amplitude in V/m: an IonisationRate or a tabulated solprop.translate.ionisation.RateTable.
fieldnumpy.ndarray: Field amplitude in V/m, shape (n_to,), sampled uniformly in time.
dtfloat: Sample spacing in seconds.
preionisationfloat, optional: Fraction already ionised before the pulse, in [0, 1], at the one position this field is sampled at. A response whose preionisation is a path along z is resolved to its value at a plane before this is called.
integralRunningIntegral or None, optional: The running-integral rule; None is RunningIntegral.trapezoid().

Returns:

numpy.ndarray: Ionised fraction in [preionisation, 1], shape (n_to,), non-decreasing in time (under the band rule, to within the filter’s undershoot).

Raises:

MaterialError: If dt is not positive and finite or preionisation is outside [0, 1].
GridError: If integral was designed for another dt.

Examples

>>> import numpy as np
>>> from solprop.nonlinear.ionisation import ionisation_fraction
>>> from solprop.nonlinear.optmat_rates import OptmatRate
>>> t = np.linspace(-2e-14, 2e-14, 129)
>>> field = 4e10 * np.exp(-((t / 5e-15) ** 2))
>>> rate = OptmatRate("Ar", 800e-9)
>>> fraction = ionisation_fraction(rate, field, float(t[1] - t[0]))
>>> f"{fraction[-1]:.4f}"
'0.2148'

solprop.nonlinear.ionisation.keldysh(gas, wavelength, field_amplitude)[source]#

Return the Keldysh parameter of a gas at a wavelength and field amplitude.

gamma = omega sqrt(2 Ip) / E in atomic units: the ratio of the laser period to the time the electron needs to tunnel. Ionisation is multiphoton for gamma >> 1 and quasistatic tunnelling for gamma << 1.

Parameters:

gasGas or str: The species; its first ionisation energy is optmat’s (ionisation_potential()).
wavelengthfloat: Driving wavelength in metres.
field_amplitudearray_like or float: Field amplitude in V/m; the sign is ignored.

Returns:

numpy.ndarray: Keldysh parameter, dimensionless, with the shape of field_amplitude.

Raises:

MaterialError: If the gas is unknown, has no ionisation potential, or the wavelength or field is not positive and finite.

Examples

>>> from solprop.nonlinear.ionisation import keldysh
>>> f"{keldysh('Ar', 800e-9, 5e10):.3f}"
'0.630'

solprop.nonlinear.ionisation.rate_convention(rate)[source]#

Read an explicit optical-cycle convention without guessing for legacy rates.

Parameters:

rateEvaluableRate: Magnitude rate, optionally declaring a convention attribute.

Returns:

str: 'instantaneous-isotropic', 'cycle-averaged' or 'unspecified'. Missing or unrecognised declarations are unspecified. Field-grid consumers retain their existing acceptance rules; envelope plasma requires the first.

solprop.nonlinear.ionisation.rate_energy(rate)[source]#

Read the energy per ionisation event a rate declares for the plasma loss.

A rate that carries its own energy per event (SEP 3, D8: the loss uses the energy the rate’s exponent uses, solprop.nonlinear.optmat_rates.OptmatRate) declares it as plasma_energy beside its species_identity. Every other rate declares none, and its response takes the species’ ionisation potential, as before.

Parameters:

rateEvaluableRate: Any rate.

Returns:

tuple of (str, float) or None: The species name and the energy per event in J, or None for a rate that declares no energy.

Examples

>>> from solprop.nonlinear.ionisation import rate_energy
>>> def rate(field):
...     return 0.0 * field
>>> rate_energy(rate) is None
True

solprop.nonlinear.ionisation.rate_species(rate)[source]#

Read a rate’s declared species name and binding energy.

Parameters:

rateEvaluableRate: Analytic rate with gas, or a table retaining species_identity and ionisation_potential. An anonymous callable has neither.

Returns:

tuple of (str or None, float or None): Chemical species name and first ionisation potential in J. Anonymous rates return two None values; their response supplies the material binding.