caustics

caustics.py — Catastrophe classification, ray caustic detection, and uniform asymptotic corrections

Overview

The caustics module provides a comprehensive toolkit for working with caustics in semiclassical analysis and catastrophe theory. It consolidates previously scattered functionality into a single, well‑tested module with a clean API. The main components are:

1. Arnold classification (algebraic) – Symbolic classification of critical points of a function H(ξ) (1D) or H(ξ,η) (2D) into elementary catastrophes: Morse, fold (A2), cusp (A3), swallowtail (A4), butterfly (A5), and umbilics (D4±). Uses successive derivatives and Hessian analysis.

2. Catastrophe detection (adaptive numerical solver) – Robustly finds all critical points of H by combining symbolic solving, coarse grid scanning, Newton refinement, and clustering. Then classifies them using the algebraic routines.

3. Ray‑based caustic detection – Operates on a bundle of rays (output of wkb.py) that contain the integrated stability matrix J = ∂x(t)/∂q. Caustics are detected when det(J) crosses zero. This corrects the common mistake of using velocity zeros or momentum minima as caustic indicators. Maslov indices and phase shifts are computed automatically.

4. Uniform special functions – Implements the Airy function (fold), the Pearcey integral (cusp), and a swallowtail integral (A4), together with the correct prefactors for uniform WKB approximations near caustics.

5. Visualisation helpers – Functions to plot catastrophe points on 1D curves or 2D surfaces, and to overlay caustic events on ray bundles.

The module is designed to be independent of the other packages (it does not import wkb.py, psiop.py, etc.), but it is used by wkb.py for caustic detection and correction.

Mathematical background

Caustics are envelopes of families of rays (bicharacteristics) where the ray density becomes infinite. In the semiclassical approximation, the amplitude diverges at a caustic, and the standard WKB ansatz breaks down. Uniform approximations replace the oscillatory exponential by special functions.

For a Hamiltonian system with n degrees of freedom, a family of rays is parameterised by an initial‑condition parameter q ∈ ℝⁿ. The stability matrix

J(t) = ∂x(t)/∂q

satisfies the variational equation dJ/dt = H_{px}·J along each ray. A caustic occurs when det(J(t)) = 0. The Maslov index μ is the signed number of such zero crossings; each crossing contributes a phase shift of π/2 to the wave function.

Arnold’s classification of catastrophes (singularities of gradient maps) provides a taxonomy of caustics:

• A₂ (fold): normal form ξ³ – the simplest caustic, described by the Airy function.

• A₃ (cusp): normal form ξ⁴ – described by the Pearcey integral.

• A₄ (swallowtail): normal form ξ⁵ – a three‑parameter integral.

• D₄± (umbilics): normal forms ξ³ + η³ (hyperbolic) and ξ³ – 3ξη² (elliptic) – degenerate cases with Hessian rank 0.

The module classifies a critical point (where ∇H = 0) by examining its Hessian rank and higher derivatives, using algebraic invariants.

References

[1]

Arnold, V. I. Catastrophe Theory, Springer‑Verlag, 1986.

[2]

Duistermaat, J. J. “Oscillatory integrals, Lagrange immersions and unfolding of singularities”, Comm. Pure Appl. Math. 27, 207–281, 1974.

[3]

Maslov, V. P. & Fedoriuk, M. V. Semi‑Classical Approximation in Quantum Mechanics, Reidel, 1981.

[4]

Kravtsov, Yu. A. & Orlov, Yu. I. Caustics, Catastrophes and Wave Fields, Springer, 1999.

[5]

Berry, M. V. & Howls, C. J. “High orders of the Weyl expansion for quantum billiards”, Phys. Rev. E **50**(5), 3577–3595, 1994.

[6]

Connor, J. N. L. “Practical methods for the uniform asymptotic evaluation of oscillatory integrals”, Mol. Phys. **31**(1), 33–55, 1976.

class caustics.AdaptiveCriticalPointSolver(H_expr: Expr, xi_vars: Sequence[Symbol], coords: Sequence[Symbol] | None = None, bounds: Dict | None = None, coarse_n: int = 30, cluster_tol: float = 0.0001, grad_tol: float = 1e-08)

Bases: object

Find critical points of H(xi_vars) robustly.

Strategy

1. Symbolic solve (sp.solve) — exact when polynomial.

2. Coarse grid scan — evaluate |∇H| on a grid, keep near-zero cells as seeds for Newton refinement.

3. Newton / fsolve refinement — high-accuracy convergence from seeds.

4. DBSCAN-style deduplication — merge solutions closer than cluster_tol.

This replaces the fragile uniform-grid + sp.nsolve approach that missed solutions and was slow.

param H_expr:

type H_expr:

sp.Expr

param xi_vars:

Frequency variables (the unknowns for ∇H = 0).

type xi_vars:

sequence of sp.Symbol

param coords:

Additional parameter symbols (treated as unknowns too).

type coords:

sequence of sp.Symbol, optional

param bounds:

Search bounds per symbol. Defaults to (-5, 5) for all.

type bounds:

dict {symbol: (lo, hi)}, optional

param coarse_n:

Points per axis for the coarse scan (default 30).

type coarse_n:

int

param cluster_tol:

Distance threshold for deduplication (default 1e-4).

type cluster_tol:

float

param grad_tol:

|∇H| threshold to accept a point as critical (default 1e-8).

type grad_tol:

float

solve(method: str = 'auto') → List[ndarray]

Find all critical points.

Parameters:

method ("symbolic" | "adaptive" | "auto") –

• “symbolic” : sp.solve only (fast for polynomials)

• ”adaptive” : coarse grid + Newton (always numerical)

• ”auto” : try symbolic first, fall back to adaptive

Returns:

Each array gives the values of (xi_vars + coords) at the critical point.

Return type:

list of np.ndarray, each of shape (dim,)

to_dict_list(solutions: List[ndarray] | None = None) → List[Dict]

Convert solutions to list of dicts {symbol: value}. If solutions is None, calls self.solve() first.

class caustics.CausticEvent(ray_idx: int, time: float, time_idx: int, position: ndarray, momentum: ndarray, det_J: float, sign_change: int = 0, maslov_contribution: int = 1, arnold_type: str = 'unknown')

Bases: object

A caustic crossing event along a ray.

ray_idx

Type:

int Index of the ray in the bundle.

time

Type:

float Integration time t* where det(J) = 0.

time_idx

Type:

int Index in the time array closest to t*.

position

Type:

np.ndarray Spatial position x(t*) (shape (dim,)).

momentum

Type:

np.ndarray Momentum ξ(t*) (shape (dim,)).

det_J

Type:

float Value of det(J) at t* (close to 0).

sign_change

Type:

int +1 or -1 (sign of d(det J)/dt at crossing).

maslov_contribution

Type:

int +1 (= π/2 phase shift).

arnold_type

Type:

str Classification from classify_arnold_* (if available).

arnold_type: str = 'unknown'

det_J: float

maslov_contribution: int = 1

momentum: ndarray

position: ndarray

ray_idx: int

sign_change: int = 0

time: float

time_idx: int

class caustics.CausticFunctions

Bases: object

Special functions for uniform asymptotic corrections near caustics.

All functions include the correct normalisation prefactors for the uniform WKB approximation (Maslov-Fedoriuk / Duistermaat).

Fold (A2) : Airy function Ai Cusp (A3) : Pearcey integral P(x,y) Swallowtail (A4) : SW integral (3-parameter oscillatory)

static airy_Ai(z: float | ndarray) → float | ndarray

Airy function Ai(z).

Fold caustic uniform approximation:

u(x) ≈ 2√π · ε^{1/6} · |∂_s J|^{-1/2} · Ai(-ε^{-2/3} ζ(x)) · e^{iS_c/ε}

where ζ(x) is the local coordinate measuring distance to the caustic.

static airy_Ai_prime(z: float | ndarray) → float | ndarray

Derivative Ai’(z) of the Airy function.

static airy_Bi(z: float | ndarray) → float | ndarray

Second Airy function Bi(z) (exponentially growing branch).

static cusp_uniform(x: ndarray, y: ndarray, x_c: float, y_c: float, epsilon: float, a_c: float, S_c: float) → ndarray

Uniform approximation near a cusp caustic (A3).

Scaled Pearcey:

u(x,y) ≈ ε^{1/4} · a_c · P(ε^{-1/2}(x-x_c), ε^{-3/4}(y-y_c))

· exp(i S_c / ε)

Return type:

np.ndarray of complex, same shape as x and y.

static fold_uniform(x: ndarray, x_c: float, epsilon: float, a_c: float, S_c: float, dJ_ds: float) → ndarray

Uniform approximation near a fold caustic (A2).

Formula (Duistermaat 1974, eq. 2.4):

u(x) = 2√π · ε^{1/6} · a_c · |dJ/ds|^{-1/2}

· Ai(-ε^{-2/3} (x - x_c)) · exp(i S_c / ε)

Parameters:

• x (evaluation grid)

• x_c (caustic position)

• epsilon (small parameter)

• a_c (WKB amplitude at caustic)

• S_c (phase at caustic)

• dJ_ds (∂(det J)/∂s at the caustic (slope of det J vs parameter s))

Return type:

np.ndarray of complex

static maslov_phase_shift(mu: int) → complex

WKB phase correction factor from Maslov index μ.

Each caustic crossing (each zero of det J) contributes a phase shift of π/2. The total correction factor is exp(i μ π/2).

Parameters:

mu (int Maslov index (number of signed caustic crossings).)

Return type:

complex exp(i μ π/2)

static pearcey(x: float, y: float, t_range: float = 6.0, n_pts: int = 500) → complex

Pearcey integral for cusp caustic (A3):

P(x, y) = ∫_{-∞}^{∞} exp(i(t⁴ + x t² + y t)) dt

Evaluated by adaptive quadrature on [-t_range, t_range].

Parameters:

• x (control parameters)

• y (control parameters)

• t_range (integration half-width (default 6.0; increase for large |x|,|y|))

• n_pts (number of quadrature points)

Return type:

complex

static pearcey_grid(X: ndarray, Y: ndarray, t_range: float = 6.0, n_pts: int = 500) → ndarray

Vectorised Pearcey integral over 2D grids X, Y. Returns complex array of same shape as X.

static swallowtail(x: float, y: float, z: float, t_range: float = 5.0, n_pts: int = 400) → complex

Swallowtail integral for A4 caustic:

SW(x,y,z) = ∫_{-∞}^{∞} exp(i(t⁵ + x t³ + y t² + z t)) dt

Parameters:

• x (control parameters)

• y (control parameters)

• z (control parameters)

• t_range (quadrature parameters)

• n_pts (quadrature parameters)

Return type:

complex

class caustics.RayCausticDetector(ray_bundle: List[Dict], dimension: int, det_threshold: float = 0.05, H_expr: Expr | None = None, xi_syms: Tuple | None = None, x_syms: Tuple | None = None)

Bases: object

Detect caustics in a ray bundle by tracking the stability matrix J.

Mathematical background

Given a Hamiltonian H(x, ξ) and a family of rays parameterised by initial condition s ∈ ℝ^d, the stability matrix is

J(t) = ∂x(t, s) / ∂s₀

and satisfies the variational equation along the ray:

dJ/dt = H_{px}(x(t), ξ(t)) · J

where H_{px} = ∂²H/∂ξ∂x is the mixed Hessian (a d×d matrix in d-D).

A caustic occurs at t* when det(J(t*)) = 0. The Maslov index μ counts the number of such crossings (with sign).

Correct approach vs old CausticDetector

OLD (wrong):

1D: caustic ↔ dx/dt = 0 (velocity turning point, NOT a caustic) 2D: caustic ↔ |ξ| ≈ 0 (small momentum, NOT a caustic)

NEW (correct):

1D: caustic ↔ J(t) = ∂x(t)/∂x₀ = 0 (requires integrating J) 2D: caustic ↔ det([[∂x/∂x₀, ∂x/∂y₀],[∂y/∂x₀, ∂y/∂y₀]]) = 0

param ray_bundle:

Each dict is an integrated ray as returned by wkb.py: 1D keys: ‘t’, ‘x’, ‘xi’, ‘S’, ‘a0’, … 2D keys: ‘t’, ‘x’, ‘y’, ‘xi’, ‘eta’, ‘S’, ‘a0’, … Each ray must additionally contain the stability matrix entries: 1D: ‘J’ (array, shape (n_steps,)) 2D: ‘J11’,’J12’,’J21’,’J22’ (arrays, shape (n_steps,)) If J is absent, RayCausticDetector will attempt to compute it numerically from neighbouring rays if a full bundle is provided.

type ray_bundle:

list of dict

param dimension:

1 or 2.

type dimension:

int

param det_threshold:

|det(J)| < det_threshold to flag a caustic (default 0.05). Relative to the initial value |det(J(0))|.

type det_threshold:

float

param H_expr:

Symbolic Hamiltonian. If provided, classify_arnold_* is called at each caustic to give the Arnold type.

type H_expr:

sp.Expr, optional

param xi_syms:

Momentum symbols matching H_expr.

type xi_syms:

tuple of sp.Symbol, optional

param x_syms:

Spatial symbols matching H_expr.

type x_syms:

tuple of sp.Symbol, optional

detect() → List[CausticEvent]

Detect all caustic events in the ray bundle.

Return type:

list of CausticEvent

maslov_index(ray_idx: int) → int

Return the Maslov index for ray ray_idx.

The Maslov index is the signed count of caustic crossings:

μ = Σ_k sign(d det(J)/dt at t_k*)

Each crossing contributes a phase shift of π/2 to the WKB amplitude.

Parameters:

ray_idx (int)

Return type:

int (total signed count; use abs() for the number of crossings)

maslov_phase(ray_idx: int) → float

Total Maslov phase = μ · π/2 for ray ray_idx.

static stability_matrix_ode_1d(H_px_func)

Return the RHS for the 1D stability ODE to be integrated with the ray.

Usage in your ray integrator

Add a state component J (scalar, initially 1) to the ray ODE. The extra equation is:

dJ/dt = H_px(x(t), ξ(t)) · J

where H_px = ∂²H/∂ξ∂x (scalar in 1D).

param H_px_func:

Lambdified ∂²H/∂ξ∂x.

type H_px_func:

callable (x, xi) → float

rtype:

callable (t, J_val, x_val, xi_val) → dJ/dt

static stability_matrix_ode_2d(H_px_func)

Return the RHS for the 2D stability ODE.

The 2×2 matrix J satisfies dJ/dt = H_{px} · J where H_{px} = [[∂²H/∂ξ∂x, ∂²H/∂ξ∂y],[∂²H/∂η∂x, ∂²H/∂η∂y]].

Pack J as a flat vector [J11, J12, J21, J22] in the state.

Parameters:

H_px_func (callable (x, y, xi, eta) → 2×2 array) – Lambdified H_{px} matrix.

Return type:

callable (t, J_flat, x, y, xi, eta) → dJ_flat/dt (4-vector)

summary() → Dict

Return a summary dict of detected caustic events.

caustics.classify_arnold_1d(f: Expr, xi: Symbol, point: float, max_order: int = 8, tol: float = 1e-08) → dict

Classify a 1D singularity of f(xi) at xi = point.

Algorithm: evaluate successive derivatives f^(k)(point) until the first non-vanishing one determines the A_k type.

Parameters:

• f (sympy.Expr) – Function of xi.

• xi (sympy.Symbol)

• point (float) – Critical point (should satisfy f’(point) ≈ 0).

• max_order (int) – Maximum derivative order to check (default 8 covers up to A7).

• tol (float) – Threshold for “non-zero” derivative.

Returns:

‘type’ : str e.g. “A2 (Fold)”, “A3 (Cusp)”, … ‘order’ : int k such that f^(k)(point) ≠ 0 ‘derivatives’ : dict {k: value} for k = 1..max_order ‘point’ : float

Return type:

dict with keys

caustics.classify_arnold_2d(H: Expr, xi: Symbol, eta: Symbol, point: dict, tol: float = 1e-08) → dict

Classify a 2D catastrophe at a critical point of H(xi, eta).

Follows Arnold’s classification:

Morse rank 2 → non-degenerate A_k rank 1 → fold (A2), cusp (A3), swallowtail (A4),

butterfly (A5), A6+

D4± rank 0 → hyperbolic umbilic (I>0) / elliptic umbilic (I<0) higher rank 0, I=0 → E6 or beyond

Parameters:

• H (sympy.Expr) – Function of (xi, eta).

• xi (sympy.Symbol)

• eta (sympy.Symbol)

• point (dict) – Numeric values {‘xi’: float, ‘eta’: float}. Extra keys (space coordinates) are ignored.

• tol (float)

Returns:

‘type’ : str ‘hessian’ : 2×2 list ‘third_order_tensor’ : dict of H_xxx, H_xxy, H_xyy, H_yyy ‘directional_derivatives’: dict of D3…D6 (rank-1 case) ‘cubic_invariant_I’ : float (rank-0 case)

Return type:

dict with keys

caustics.detect_catastrophes(H_expr: Expr, xi_vars: Sequence[Symbol], coords: Sequence[Symbol] | None = None, method: str = 'auto', bounds: Dict | None = None, coarse_n: int = 30, max_order: int = 6, tol: float = 1e-08) → List[Dict]

Detect and classify catastrophes of H_expr(xi_vars).

Combines AdaptiveCriticalPointSolver (robust search) with classify_arnold_1d / classify_arnold_2d (typing).

Parameters:

• H_expr (sp.Expr)

• xi_vars (tuple of sp.Symbol — (xi,) in 1D, (xi, eta) in 2D)

• coords (tuple of sp.Symbol, optional) – Additional parameter symbols (e.g. spatial coordinates x, y). When provided, catastrophes are tracked as families.

• method ("symbolic" | "adaptive" | "auto")

• bounds (dict {symbol: (lo, hi)}, optional)

• coarse_n (int — coarse grid resolution per axis)

• max_order (int — max derivative order for 1D classification)

• tol (float)

Returns:

‘point’ : dict {symbol: value} ‘type’ : str Arnold type ‘details’ : dict hessian / derivatives / invariants

Return type:

list of dict, each with

caustics.find_critical_points_numerical(grad_func, initial_guesses: List[ndarray], tolerance: float = 1e-06, domain: List[Tuple[float, float]] | None = None, cluster_tol: float = 0.0001) → List[ndarray]

Shared numerical kernel for locating critical points where ∇φ = 0.

Minimises |∇φ(x)|² from each initial guess using L-BFGS-B, then deduplicates nearby solutions with a greedy DBSCAN-style clustering.

This function is the single source of truth used by:

• asymptotic.Analyzer.find_critical_points

• fio_bridge._BoundAnalyzer.find_critical_points

It can also be called directly for lightweight use cases.

Parameters:

• grad_func (callable) – Function grad_func(*x) -> array-like returning the gradient ∇φ evaluated at the point x (unpacked coordinates).

• initial_guesses (list of np.ndarray) – Starting points for the optimisation.

• tolerance (float) – Convergence threshold: a point is accepted when |∇φ|² < tolerance.

• domain (list of (lo, hi) tuples, optional) – If provided, only points satisfying lo <= x_i <= hi are kept.

• cluster_tol (float) – Distance threshold for deduplication; solutions closer than this are merged into a single representative point (cluster mean).

Returns:

Unique critical points found within the given tolerance and domain.

Return type:

list of np.ndarray

caustics.plot_catastrophe(H: Expr, xi_vars: Sequence[Symbol], points: List[Dict], xi_bounds: Tuple[float, float] = (-3.0, 3.0), eta_bounds: Tuple[float, float] = (-3.0, 3.0), n: int = 300, title: str | None = None) → None

Plot H and mark classified catastrophe points.

Parameters:

• H (sympy.Expr)

• xi_vars ((xi,) or (xi, eta))

• points (list of dicts as returned by detect_catastrophes()) – each must have keys ‘point’ and ‘type’

• xi_bounds (axis ranges)

• eta_bounds (axis ranges)

• n (grid resolution)

• title (optional plot title)

caustics.plot_caustic_events(ray_bundle: List[Dict], events: List[CausticEvent], dimension: int, n_rays_plot: int = 20, title: str | None = None) → None

Overlay caustic events on the ray bundle plot.

Parameters:

• ray_bundle (list of ray dicts (with 'x', 't' etc.))

• events (list of CausticEvent from RayCausticDetector.detect())

• dimension (1 or 2)

• n_rays_plot (max number of rays to draw)