propagator

propagator.py — Semiclassical (Van Vleck–Pauli–Morette) wavefunction

Overview

This module assembles a semiclassical (WKB / Van Vleck–Pauli–Morette) wavefunction from a fan of classical rays. It is the top-level physics layer of the psipy package and depends exclusively on riemannian.py for geometry and Jacobi-field integration, and on symplectic.py for symplectic ray tracing and action accumulation.

The public entry point is compute_wavefunction(). The result is a WKBResult dataclass that bundles the gridded wavefunction, all per-ray data, and the scattered raw data ready for plotting or further analysis. Four visualisation helpers are provided:

plot_wavefunction() — master figure (density, phase, Re/Im, rays)
plot_ray_fan() — rays coloured by action, caustics marked
plot_interference_detail()— fringes, density, and action scatter

The module also supports three types of first‑order‑in‑time (or second‑order) partial differential equations, selected by the equation argument:

Schrödinger (default) — −i ∂ψ/∂t = ψOp(p, ψ) (Van Vleck propagator)
Parabolic (heat) — ∂u/∂t = ψOp(p, u) (real exponent)
Wave (hyperbolic) — ∂²u/∂t² = ψOp(p, u) (two branches, ±√H)

Physical Background

The Van Vleck–Pauli–Morette (VVP) propagator gives the semiclassical approximation to the quantum propagator K(x, x₀; t) in the limit ℏ → 0. For a system whose classical Hamiltonian is H = ½ gⁱʲ(x) pᵢ pⱼ (geodesic motion on a Riemannian manifold with metric tensor g), the wavefunction emanating from a point source at x₀ is

ψ(x, t) = Σ_k A_k(x) · exp(i S_k(x, t)/ℏ − i μ_k π/2)

where the sum runs over all classical paths (rays) k that connect x₀ to x in time t, and:

S_k(x, t) = ∫₀ᵗ p · ẋ dt’ (Hamilton’s principal function / action) = ∫₀ᵗ g_{ij}(x) vⁱ vʲ dt’ (on a pure-metric Hamiltonian)

A_k(x) = 1 / √|det J_k(x, t)| (Van Vleck amplitude)

det J_k = det(∂x / ∂p₀) (Jacobi determinant, a.k.a. Van Vleck determinant)

μ_k = number of caustic crossings (Maslov index). Each sign change of det J contributes +1 to μ_k, adding a phase factor exp(−iπ/2) = −i per crossing.

The amplitude A_k diverges when det J = 0, i.e. at caustics, which are the envelopes of the ray family. Near a caustic the WKB approximation breaks down and is replaced by a uniform asymptotic approximation:

In 1D, a fold caustic is corrected with the Airy function Ai(ξ) where ξ(x) = (α/2ℏ)^{1/3} (x − x_c). The fringe spacing ∝ ℏ^{1/3} is physically correct, and the uniform amplitude is ψ(x) ≈ 2π a_c ℏ^{1/6} |α|^{-1/3} · Ai(ξ(x)) · exp(i S_c/ℏ).
In 2D, fold caustics are patched with Airy functions in the transverse direction n̂ = ∇det J / |∇det J|, blended with a 2D Gaussian taper. Cusp caustics (where ∇det J ≈ 0) are treated with the Pearcey integral via the asymptotic module, giving an O(ℏ^{1/4}) scaling.

For the parabolic (heat) equation ∂u/∂t = ψOp(p, u), the WKB ansatz has a real exponent, and fold caustics are corrected with the parabolic cylinder function D_{-1/2}(ζ), ζ = (α/ℏ)^{1/4} (x − x_c). For the wave equation ∂²u/∂t² = ψOp(p, u), the dispersion relation splits into two branches, and the wavefunction is the coherent sum of both families.

Connection to the Metric

For a kinetic Hamiltonian H = ½ gⁱʲ pᵢ pⱼ the canonical momentum is pᵢ = g_{ij} vʲ (covariant, lowered by the metric), so the action becomes ∫ p · v dt = ∫ g_{ij} vⁱ vʲ dt. The inverse metric gⁱʲ governs the Hamiltonian equations of motion (Hamilton’s equations), while the metric g_{ij} maps velocities to momenta. This distinction is crucial for the action fallback: when no explicit momentum is stored in the trajectory, the code reconstructs pᵢ = g_{ij}(x) vʲ and integrates ∫ g_{ij} vⁱ vʲ dt, which is exact for any Riemannian metric. The old ∫ v² dt fallback (valid only for the flat unit‑mass case) is no longer used.

The Jacobi equation governing the evolution of the Jacobi field J = ∂x/∂p₀ along a geodesic is the geodesic deviation equation, which is curvature‑ dependent:

D²J/dt² + R(J, ẋ)ẋ = 0

where R is the Riemann curvature tensor. On a flat metric (R = 0) the Jacobi field grows linearly: J(t) = t · K₀. Curvature causes focussing (det J → 0) and defocussing.

Module Architecture

Dependency tree:

compute_wavefunction          ← public entry point
├── _build_hamiltonian_sym    — builds H = ½ gⁱʲ pᵢ pⱼ symbolically
├── hamiltonian_flow          — symplectic ray integration  (symplectic.py)
├── _det_J_from_jacobi        — Jacobi determinant along each ray
│   ├── _det_J_1d             — 1D: variational ODE via solve_ivp
│   └── jacobi_equation_solver— 2D: two Jacobi fields (riemannian.py)
├── _cumulative_action        — action integral ∫ p · v dt (exact for metric)
├── _maslov_index             — count sign changes of det J
└── van_vleck_sum             — coherent sum onto output grid
    ├── _asymptotic_correction_1d  — Airy patch at 1D fold caustics
    │   └── _airy_argument         — ξ(x) = (α/2ℏ)^{1/3} (x − x_c)
    └── _asymptotic_correction_2d  — Airy / Pearcey at 2D caustics

The entry point accepts either a Metric object (geodesic motion) or a general SymPy Hamiltonian. For metric mode, the fan of initial conditions is specified as velocities vⁱ (contravariant) via the v_fan argument; the module automatically converts them to canonical momenta pᵢ = g_{ij} vʲ before integrating. For a general Hamiltonian the fan is given directly as canonical momenta via p_fan.

Result dataclasses:

RayData    — per-ray: trajectory dict, det J array, S array, Maslov μ
WKBResult  — full output: gridded ψ, raw scattered data, all RayData

Package Dependencies

riemannian.py

Metric: Encodes the Riemannian metric tensor g_{ij}(x) as a SymPy expression. Provides symbolic and numerical evaluation of g, g⁻¹, and their derivatives. Used to convert velocities to momenta (p = g · v) and to build the kinetic Hamiltonian.
geodesic_solver(): Integrates the geodesic equations ẍ + Γ vv = 0 forward in time, returning position and velocity arrays. Used as a lightweight alternative to hamiltonian_flow when symplectic accuracy is not required.
jacobi_equation_solver(): Integrates the Jacobi (geodesic deviation) equation along a given geodesic for a specified initial variation (J₀, DJ₀). Returns the Jacobi field components J_x(t), J_y(t). Called twice per ray in 2D to form the 2×2 Jacobi matrix whose determinant gives the Van Vleck amplitude.

symplectic.py

hamiltonian_flow(): Integrates Hamilton’s equations (ẋ = ∂H/∂p, ṗ = −∂H/∂x) using a symplectic integrator (Störmer–Verlet by default, or RK45). Returns a trajectory dict containing both positions and canonical momenta, which are used directly for the action integral.

asymptotic.py

Analyzer / AsymptoticEvaluator: Evaluate oscillatory integrals I(λ) = ∫ a(t) exp(iλ φ(t)) dt via stationary-phase methods. Used here only for cusp (Pearcey) caustics in 2D, where the quartic normal-form phase φ(t) = t⁴/4 requires the specialised Pearcey evaluator.

Typical Usage

import sympy as sp
import numpy as np
from riemannian import Metric
from propagator import compute_wavefunction, plot_wavefunction

# 1. Define the geometry via a Metric object
x = sp.Symbol('x', real=True)
metric = Metric(1, (x,))          # flat 1D metric, g = 1

# 2. Define the source point and a fan of initial velocities
source = (0.0,)
v_fan  = np.linspace(-4.0, 4.0, 80)

# 3. Run the full pipeline
result = compute_wavefunction(
    metric    = metric,
    source    = source,
    v_fan     = v_fan,
    t_max     = 2.0,
    hbar      = 0.1,
    n_steps   = 500,
    N_grid    = 400,
    integrator= 'verlet',
)

# 4. Visualise
import matplotlib.pyplot as plt
plot_wavefunction(result, log_scale=True)
plt.show()

For a curved metric built from a Hamiltonian H = p²/(2 m(x)):

x, p = sp.symbols('x p', real=True, positive=True)
metric = Metric.from_hamiltonian(p**2 / (2 / x**2), (x,), (p,))
# metric.g_expr == x**2

For a general Hamiltonian (e.g. harmonic oscillator):

x, xi = sp.symbols('x xi', real=True)
H = xi**2/2 + x**2/2
result = compute_wavefunction(
    hamiltonian = H,
    coords      = (x,),
    momenta     = (xi,),
    source      = (0.0,),
    p_fan       = np.linspace(-2, 2, 80),
    t_max       = 5.0,
    hbar        = 0.2,
)

References

Van Vleck, J.H. (1928). “The correspondence principle in the statistical interpretation of quantum mechanics”. Proc. Natl. Acad. Sci. 14, 178.
Morette, C. (1951). “On the definition and approximation of Feynman’s path integrals”. Phys. Rev. 81, 848.
Gutzwiller, M.C. (1990). Chaos in Classical and Quantum Mechanics. Springer, New York. (Chapter 12: the semiclassical Green’s function.)
Maslov, V.P. & Fedoriuk, M.V. (1981). Semi-Classical Approximation in Quantum Mechanics. Reidel, Dordrecht. (Maslov index and caustics.)
Berry, M.V. & Mount, K.E. (1972). “Semiclassical approximations in wave mechanics”. Rep. Prog. Phys. 35, 315. (Uniform Airy approximation.)

class propagator.EquationType[source]

Bases: object

Selector for the type of PDE solved by compute_wavefunction().

Three first-order-in-time (or second-order) equations are supported. Pass one of these three string constants as the equation argument:

EquationType.SCHRODINGER (default, original behaviour)

−i ∂u/∂t = ψOp(p, u)

WKB ansatz: u = A exp( i S/ℏ ) Phase factor at caustics: exp(−i μ π/2) (Maslov). Caustic correction: Airy function Ai(ξ).

EquationType.PARABOLIC

∂u/∂t = ψOp(p, u)

WKB ansatz: u = A exp( S/ℏ ) (real exponent — no i) The action S is accumulated with the same classical rays but enters the exponent without the imaginary unit. Consequently:

The solution is real-valued (when the initial data are real).
There are no oscillatory fringes; instead, the solution concentrates exponentially around rays with the largest action.
Caustic corrections use the parabolic cylinder function D_{-1/2} (the real-axis analogue of the Airy function for fold caustics).
No Maslov phase — sign changes of det J contribute a real factor |det J|^{-1/2} that diverges, patched by D_{-1/2}.

EquationType.WAVE

∂²u/∂t² = ψOp(p, u)

The dispersion relation is ω² = H(x, p), giving two branches: ω₊ = +√H and ω₋ = −√H. For each initial ray direction the code integrates two Hamiltonians, H₊ = +√H and H₋ = −√H, effectively doubling the ray fan. The wavefunction is the coherent sum:

u = Σ_{k,±} A_k exp( i S_k^±/ℏ − i μ_k^± π/2 )

where S^± = ∫ p · ẋ dt along rays driven by H₊ / H₋.

This covers: * The scalar wave equation □u = 0 (H = −|p|²) * Acoustics in an inhomogeneous medium (H = −c²(x)|p|²) * Any second-order hyperbolic operator.

Usage:

result = compute_wavefunction(
    ...,
    equation = EquationType.WAVE,
)

PARABOLIC = 'parabolic'

SCHRODINGER = 'schrodinger'

WAVE = 'wave'

class propagator.RayData(traj: dict, det_J: ndarray, S_cum: ndarray, mu: int)[source]

Bases: object

All data associated with a single classical ray.

A ray is a solution of Hamilton’s equations starting from the source point with one particular initial velocity v₀. This dataclass bundles the raw trajectory together with the derived semiclassical quantities needed to evaluate the Van Vleck amplitude and phase.

traj

Trajectory dictionary returned by symplectic.hamiltonian_flow(). Keys depend on dimension:

1D: 't', 'x' (or the SymPy coord name), 'xi' (canonical momentum p = g v).
2D: 't', 'x', 'y' (or coord names), 'xi', 'eta'.

All values are 1D NumPy arrays of length n_steps.

Type:: dict

det_J

Jacobi determinant det(∂x/∂p₀) along the ray.

Positive away from caustics.
Changes sign at each caustic crossing (focal point).
det_J[0] = 0 by construction (the ray fan starts at a point source, so all rays share the same initial position — the Jacobi field starts from zero separation).

Type:: np.ndarray, shape (n_steps,)

S_cum

Cumulative action S(t) = ∫₀ᵗ pᵢ ẋⁱ dt′, evaluated at each time step. For a pure-metric Hamiltonian this equals ∫₀ᵗ g_{ij} vⁱ vʲ dt′ = 2E t (twice the kinetic energy times elapsed time) on flat metrics, but differs in general on curved ones.

Type:: np.ndarray, shape (n_steps,)

mu

Maslov index: the number of caustic crossings (sign changes of det_J) accumulated along the ray from t=0 to t=t_max. Each crossing contributes a phase factor exp(−iπ/2) = −i to the semiclassical amplitude.

Type:: int

S_cum: ndarray

det_J: ndarray

mu: int

traj: dict

class propagator.WKBResult(rays: List[RayData], X: ndarray, Y: ndarray | None, psi: ndarray, x_pts: ndarray, y_pts: ndarray | None, S_pts: ndarray, det_J_pts: ndarray, mu_pts: ndarray, hbar: float, t_max: float, dim: int, equation: str = 'schrodinger')[source]

Bases: object

Full output of compute_wavefunction().

Bundles the gridded semiclassical wavefunction together with all per-ray data and the raw scattered point cloud, making it self-contained for plotting, post-processing, or archiving.

rays

One RayData per successfully integrated ray. Failed rays (e.g. due to numerical blow-up) are silently dropped.

Type:: list of RayData

X

x-coordinates of the output grid.

1D: shape (N_grid,) — a 1D array of x positions.
2D: shape (N_grid, N_grid) — the x-component of a meshgrid, as returned by np.meshgrid.

Type:: np.ndarray

Y

y-coordinates of the output grid (2D only); None in 1D.

Type:: np.ndarray or None

psi

Semiclassical wavefunction on the output grid.

1D: shape (N_grid,).
2D: shape (N_grid, N_grid).

Assembled by van_vleck_sum() as the coherent sum over all rays, with Airy corrections applied near caustics.

Type:: np.ndarray (complex)

x_pts

x-coordinates of all trajectory points from all rays, concatenated. These are the scattered source points fed to scipy.interpolate.griddata.

Type:: np.ndarray, shape (n_rays × n_steps,)

y_pts

y-coordinates of all trajectory points (2D only); None in 1D.

Type:: np.ndarray or None

S_pts

Cumulative action at each scattered point (concatenated from all rays).

Type:: np.ndarray, shape (n_rays × n_steps,)

det_J_pts

Jacobi determinant at each scattered point.

Type:: np.ndarray, shape (n_rays × n_steps,)

mu_pts

Maslov index broadcast to every step of each ray (constant per ray, since μ is the total caustic count for that ray’s trajectory).

Type:: np.ndarray of int, shape (n_rays × n_steps,)

hbar

Reduced Planck constant used in the computation.

Type:: float

t_max

Integration time of the ray fan.

Type:: float

dim

Spatial dimension, 1 or 2.

Type:: int

S_pts: ndarray

X: ndarray

Y: ndarray | None

det_J_pts: ndarray

dim: int

equation: str = 'schrodinger'

hbar: float

mu_pts: ndarray

psi: ndarray

rays: List[RayData]

t_max: float

x_pts: ndarray

y_pts: ndarray | None

propagator.animate_wavefunction(result: WKBResult, frame_times=None, n_frames: int = 100, interval: int = 50, save_path: str | None = None, log_scale: bool = True) → matplotlib.animation.FuncAnimation[source]

Animate the time evolution of the semiclassical wavefunction.

Overview

This function creates an animation of ψ(x, t) (or u(x, t) for parabolic equation) by re‑assembling the Van Vleck–Pauli–Morette sum at each frame from the already integrated ray data. Unlike a static snapshot that is merely rotated in phase, this recomputes the full coherent superposition at each time slice, correctly capturing the spreading, focusing, caustic births, and true |ψ|² dynamics.

The animation calls the appropriate summation routine (van_vleck_sum(), parabolic_sum(), or wave_sum()) for each frame. No ray integration is repeated – only the gridding step is performed, so the overhead is modest.

Physical Background

For the Schrödinger equation (equation = EquationType.SCHRODINGER) the wavefunction is

ψ(x,t) = Σ_k A_k(t) · exp(i S_k(t)/ℏ − i μ_k(t) π/2)

with A_k = 1/√|det J_k|. For the heat‑type equation (EquationType.PARABOLIC) the sum is real‑exponential

u(x,t) = Σ_k A_k(t) · exp(S_k(t)/ℏ)

and for the wave equation (EquationType.WAVE) it is the coherent sum over two branches H₊ = +√H and H₋ = −√H.

param result:: Output of compute_wavefunction(). Must contain the full per‑ray time series (traj, det_J, S_cum). The animation uses the raw data stored in result.rays to evaluate the wavefunction at arbitrary times within [0, result.t_max].
type result:: WKBResult
param frame_times:: Physical times at which to render frames. Each value is matched to the nearest stored time step. If None, n_frames equally spaced times between 0 and result.t_max are used.
type frame_times:: array‑like or None
param n_frames:: Number of frames (ignored when frame_times is supplied).
type n_frames:: int, default 100
param interval:: Delay between frames in milliseconds.
type interval:: int, default 50
param save_path:: If given, save the animation to this path. Format is inferred from the extension: .gif uses Pillow, .mp4 uses ffmpeg.
type save_path:: str or None, default None
param log_scale:: If True, display the density as log(1 + |ψ|²) (or log(1 + u²) for parabolic) to reveal low‑amplitude features. Otherwise show the linear density.
type log_scale:: bool, default True
returns:: anim – The animation object. Call plt.show() to display it, or use anim.save(...) directly (though the function already supports saving via the save_path argument).
rtype:: matplotlib.animation.FuncAnimation

Notes

Axis limits for Re(ψ), Im(ψ), and the density are pre‑scanned over all frames to ensure a stable scale throughout the animation.
For 1D the figure layout mirrors plot_wavefunction(): density, phase, real part, imaginary part (the latter hidden for parabolic).
For 2D the layout consists of real part, imaginary part, density, and phase (or log|u| for parabolic).
For the wave equation the ray fan is assumed to be split into two halves representing the H₊ and H₋ branches (the order used by compute_wavefunction()). Both are combined in the sum at each frame.
The function uses the appropriate summation method based on the result.equation attribute.

Examples

import numpy as np
import sympy as sp
from riemannian import Metric
from propagator import compute_wavefunction, animate_wavefunction

# 1D free particle on a flat metric
x = sp.Symbol('x', real=True)
metric = Metric(1, (x,))
source = (0.0,)
v_fan  = np.linspace(-4.0, 4.0, 100)

result = compute_wavefunction(
    metric=metric, source=source, v_fan=v_fan,
    t_max=2.0, hbar=0.05, n_steps=500, N_grid=400,
)

# Create animation and display
anim = animate_wavefunction(result, n_frames=80, interval=40)
plt.show()

# Save as GIF
animate_wavefunction(result, save_path="wavepacket.gif")

propagator.compute_wavefunction(metric: Metric | None = None, source: Tuple | None = None, v_fan: ndarray | None = None, t_max: float | None = None, hbar: float = 1.0, n_steps: int = 400, N_grid: int = 300, xlim: Tuple | None = None, ylim: Tuple | None = None, integrator: str = 'verlet', hamiltonian: Expr | None = None, coords: Tuple | None = None, momenta: Tuple | None = None, p_fan: ndarray | None = None, parallel: bool = True, equation: str = 'schrodinger') → WKBResult[source]

Compute the semiclassical (Van Vleck–Pauli–Morette) wavefunction.

This is the main public entry point. It accepts two distinct input modes depending on whether you supply a Metric object (pure kinetic, geodesic motion) or an explicit SymPy Hamiltonian expression (general T + V systems).

The returned wavefunction psi is evaluated only at the final time t = t_max using the endpoint of each classical ray. (Full ray histories are stored in result.rays and can be used for animation.)

Input Modes

Mode A — Metric (original interface, v_fan required):

Pass a riemannian.Metric object. The Hamiltonian is built internally as H = ½ gⁱʲ pᵢ pⱼ. Initial momenta are obtained by converting the supplied velocity fan: p₀ = g(x₀) · v₀.

result = compute_wavefunction(
    metric = Metric(1, (x,)),
    source = (0.0,),
    v_fan  = np.linspace(-3, 3, 60),
    t_max  = 2.0,
)

Mode B — General Hamiltonian (new interface, p_fan required):

Pass a SymPy expression H(coords, momenta) together with the coordinate and momentum symbol tuples. Initial conditions are specified directly as a fan of canonical momenta p₀.

x, xi = sp.symbols('x xi', real=True)
H = xi**2 / 2 + sp.cos(x)        # pendulum-type Hamiltonian
result = compute_wavefunction(
    hamiltonian = H,
    coords      = (x,),
    momenta     = (xi,),
    source      = (0.0,),
    p_fan       = np.linspace(-2, 2, 60),
    t_max       = 3.0,
)

param metric:

Riemannian metric (Mode A). Mutually exclusive with hamiltonian.

type metric:

Metric or None

param source:

Initial position of the point source: (x₀,) or (x₀, y₀).

type source:

tuple of float

param v_fan:

Fan of initial velocities (Mode A only).

1D: shape (n_rays,)
2D: shape (n_rays, 2)

type v_fan:

np.ndarray or None

param t_max:

Total integration time.

type t_max:

float

param hbar:

Reduced Planck constant.

type hbar:

float, default 1.0

param n_steps:

Number of time steps per ray.

type n_steps:

int, default 400

param N_grid:

Output grid resolution.

type N_grid:

int, default 300

param xlim:

x-extent of the output grid (auto-detected if None).

type xlim:

tuple or None

param ylim:

y-extent of the output grid (auto-detected if None, 2D only).

type ylim:

tuple or None

param integrator:

Symplectic integrator: 'verlet' or 'rk45'.

type integrator:

str, default 'verlet'

param hamiltonian:

General SymPy Hamiltonian H(coords, momenta) (Mode B). Mutually exclusive with metric.

type hamiltonian:

sp.Expr or None

param coords:

Position symbols, e.g. (x,) or (x, y). Required in Mode B.

type coords:

tuple of sp.Symbol or None

param momenta:

Momentum symbols, e.g. (xi,) or (xi, eta). Required in Mode B.

type momenta:

tuple of sp.Symbol or None

param p_fan:

Fan of initial canonical momenta (Mode B only).

1D: shape (n_rays,)
2D: shape (n_rays, 2)

type p_fan:

np.ndarray or None

param parallel:

If True, use multiprocessing to integrate rays in parallel.

type parallel:

bool, default True

param equation:

Type of PDE: 'schrodinger', 'parabolic', or 'wave'.

type equation:

str, default EquationType.SCHRODINGER

returns:

Dataclass containing: - psi : wavefunction on the grid at t = t_max. - X, Y : grid coordinates. - rays : list of RayData with full time histories. - x_pts, y_pts : scattered positions of all ray points

(useful for debugging, but note that the static psi uses only the final points).

S_pts, det_J_pts, mu_pts : corresponding action, Jacobi determinant and Maslov index.
hbar, t_max, dim, equation.

rtype:

WKBResult

raises ValueError:

If the input mode cannot be resolved.

raises RuntimeError:

If every ray in the fan fails to integrate.

propagator.parabolic_sum(pts: ndarray, S: ndarray, det_J: ndarray, xlim: Tuple[float, float], ylim: Tuple[float, float] | None = None, N: int = 300, hbar: float = 1.0, reg: float = 0.0001, method: str = 'linear', caustic_threshold: float = 0.05) → Tuple[ndarray, ndarray, ndarray | None][source]

Assemble the semiclassical solution for the parabolic (heat-type) equation ∂u/∂t = ψOp u on a regular grid.

WKB formula

Unlike the Schrödinger case there is no imaginary unit in the exponent:

u_k(x) = exp( S_k(x)/ℏ ) / √max(|det J_k|, reg)

The real exponential means:

Rays with larger action dominate exponentially.
There are no oscillatory fringes between ray contributions.
At caustics the amplitude diverges as in the Schrödinger case, but is patched with the parabolic cylinder function D_{-1/2} instead of Ai.
The Maslov index is irrelevant (no phase); sign changes of det J are handled by the absolute value and the caustic patch.
mu is not required for the parabolic equation.

The interpolation strategy is identical to van_vleck_sum().

param pts:: Same meaning as in van_vleck_sum().
param S:: Same meaning as in van_vleck_sum().
param det_J:: Same meaning as in van_vleck_sum().
param xlim:: Same meaning as in van_vleck_sum().
param ylim:: Same meaning as in van_vleck_sum().
param N:: Same meaning as in van_vleck_sum().
param hbar:: Same meaning as in van_vleck_sum().
param reg:: Same meaning as in van_vleck_sum().
param method:: Same meaning as in van_vleck_sum().
param caustic_threshold:: Same meaning as in van_vleck_sum().
returns:: u, X, Y
rtype:: same layout as van_vleck_sum().

propagator.plot_interference_detail(result: WKBResult, save_path=None) → Figure[source]

Three-panel diagnostic figure focussing on interference and phase structure.

Panels (left to right)

Re(ψ) — interference fringes The real part of the wavefunction, which directly shows the fringe pattern. In 1D: line plot with filled area. In 2D: pcolormesh with RdBu_r diverging colourmap.
|ψ|² — probability density The squared modulus, showing where the quantum particle is likely to be found. In 1D: filled area with inferno colourmap. In 2D: pcolormesh with inferno colourmap.
S(x) coloured by Maslov index μ A scatter plot of the raw action values S at each ray trajectory point versus position x, coloured by the Maslov index μ of the corresponding ray (RdYlGn colourmap: green = μ=0, yellow = μ=1, red = μ≥2). This reveals how multiple sheets of the Lagrangian manifold (rays with the same x but different action) contribute to the interference pattern.

param result:: Output of compute_wavefunction().
type result:: WKBResult
param save_path:: If given, save the figure to this path at 150 dpi.
type save_path:: str or None, default None
returns:: fig – Three-panel figure, 16 × 5 inches.
rtype:: matplotlib.figure.Figure

propagator.plot_ray_fan(result: WKBResult, save_path=None) → Figure[source]

Plot the ray fan coloured by accumulated action, with caustics highlighted.

Each ray is drawn as a thin line whose colour is taken from the viridis colourmap, mapped linearly from the minimum to the maximum final action S(t_max) across all rays. Points where det J changes sign (caustic crossings) are marked with yellow dots.

In 1D the horizontal axis is time t and the vertical axis is position x(t). In 2D the axes are the spatial coordinates x and y, showing the geometric ray pattern in configuration space.

A colourbar on the right indicates the action scale.

Parameters:

result (WKBResult) – Output of compute_wavefunction().
save_path (str or None, default None) – If given, save the figure to this path at 150 dpi.

Returns:

fig – Single-panel figure, 10 × 6 inches.

Return type:

matplotlib.figure.Figure

propagator.plot_wavefunction(result: WKBResult, log_scale=True, save_path=None) → Figure[source]

Master visualisation figure for a WKBResult.

Dispatches to _plot_1d() or _plot_2d() based on result.dim. Both produce a dark-themed multi-panel figure.

1D layout (4 panels, 16 × 8 inches)

Top-left — Probability density |ψ|² (or log(1 + |ψ|²)).
Top-right — Phase arg(ψ) in [−π, π].
Bottom-left — Re(ψ) and Im(ψ) overlaid.
Bottom-right — Ray fan x(t) coloured by mean |det J|.

2D layout (5 panels, 20 × 8 inches)

Top-left — Probability density (pcolormesh, inferno colourmap).
Top-right — Phase map (pcolormesh, hsv colourmap, range [−π, π]).
Bottom-left — Ray fan in (x, y) space; caustic points in yellow.
Bottom-centre — Scatter of log(1 + |det J|) over all ray points.
Bottom-right — Maslov index μ scatter over all ray points.

param result:: Output of compute_wavefunction().
type result:: WKBResult
param log_scale:: If True, display log(1 + |ψ|²) instead of |ψ|² to reveal low-amplitude features (shadow regions, secondary fringes).
type log_scale:: bool, default True
param save_path:: If given, save the figure to this path at 150 dpi with tight bounding box. The figure is returned regardless.
type save_path:: str or None, default None
returns:: fig
rtype:: matplotlib.figure.Figure

propagator.van_vleck_sum(pts: ndarray, S: ndarray, det_J: ndarray, mu: ndarray, xlim: Tuple[float, float], ylim: Tuple[float, float] | None = None, N: int = 300, hbar: float = 1.0, reg: float = 0.0001, method: str = 'linear', caustic_threshold: float = 0.05) → Tuple[ndarray, ndarray, ndarray | None][source]

Assemble the Van Vleck–Pauli–Morette wavefunction on a regular grid.

This function takes the scattered output of the ray tracing (positions, actions, Jacobi determinants, and Maslov indices) and produces the gridded semiclassical wavefunction via a two-pass hybrid scheme:

Pass 1 — WKB everywhere

For each scattered point k compute the complex WKB contribution:

ψ_k = exp(i S_k/ℏ − i μ_k π/2) / √max(|det J_k|, reg)

Then interpolate Re(ψ_k) and Im(ψ_k) separately onto the output grid.

1D: np.interp on the sorted ray positions (fast, O(M log M)).
2D: scipy.interpolate.griddata with Delaunay triangulation (method='linear' by default; 'cubic' or 'nearest' also supported). Requires at least 3 non-collinear scattered points.

The regularisation floor reg prevents division by zero at exact caustics; it has negligible effect away from caustics where |det J| ≫ reg.

Pass 2 — Airy / Pearcey corrections at caustics

1D fold — _asymptotic_correction_1d() evaluates the pointwise Airy profile Ai(ξ(x)) with ξ = (α/2ℏ)^{1/3}(x − x_c), where α is the local slope of det J. The patch is blended into the WKB grid wherever |patch| > 0.
2D fold — _asymptotic_correction_2d() applies the same Airy profile in the transverse direction n̂ = ∇det J / |∇det J|, with a 2D Gaussian taper.
2D cusp (Pearcey) — detected when |∇det J| < 1e-10; handled via asymptotic.Analyzer with quartic normal-form phase.

Performance note

The Delaunay triangulation in 2D is computed once per call and is O(M log M). The Airy corrections are applied only on the small caustic subset of scattered points, so the overhead is negligible for caustic_threshold ≤ 0.1. For problems with many dense caustic clusters, reducing caustic_threshold (e.g. to 0.01) limits the patching to the immediate caustic zones and speeds up Pass 2.

param pts:

Scattered ray positions, shape (M, 1) in 1D or (M, 2) in 2D. These are the raw trajectory points from all rays, concatenated.

type pts:

np.ndarray

param S:

Cumulative action at each scattered point.

type S:

np.ndarray, shape (M,)

param det_J:

Jacobi determinant at each scattered point. May be positive or negative; the WKB amplitude uses |det_J|.

type det_J:

np.ndarray, shape (M,)

param mu:

Maslov index at each scattered point (constant within a ray).

type mu:

np.ndarray of int, shape (M,)

param xlim:

x-extent of the output grid.

type xlim:

tuple (x_min, x_max)

param ylim:

y-extent of the output grid. None selects 1D mode.

type ylim:

tuple (y_min, y_max) or None

param N:

Grid resolution. Output has shape (N,) in 1D or (N, N) in 2D.

type N:

int, default 300

param hbar:

Reduced Planck constant, used in the phase and Airy argument.

type hbar:

float, default 1.0

param reg:

Regularisation floor for the WKB amplitude: amp = 1/√max(|det J|, reg). Should be much smaller than the typical |det J| away from caustics.

type reg:

float, default 1e-4

param method:

Interpolation method passed to scipy.interpolate.griddata (2D only). Options: 'linear', 'nearest', 'cubic'.

type method:

str, default 'linear'

param caustic_threshold:

Relative threshold for caustic detection: |det J| / max|det J| < caustic_threshold → apply Airy patch. Increase to broaden the caustic zone; decrease to restrict patching to the immediate singularity.

type caustic_threshold:

float, default 0.05

returns:

psi (np.ndarray (complex)) – Semiclassical wavefunction. Shape (N,) in 1D or (N, N) in 2D.
X (np.ndarray) – x-grid coordinates. Shape (N,) in 1D or (N, N) in 2D (meshgrid).
Y (np.ndarray or None) – y-grid coordinates (2D only); None in 1D.

propagator.wave_sum(pts_plus: ndarray, S_plus: ndarray, det_J_plus: ndarray, mu_plus: ndarray, pts_minus: ndarray, S_minus: ndarray, det_J_minus: ndarray, mu_minus: ndarray, xlim: Tuple[float, float], ylim: Tuple[float, float] | None = None, N: int = 300, hbar: float = 1.0, reg: float = 0.0001, method: str = 'linear', caustic_threshold: float = 0.05) → Tuple[ndarray, ndarray, ndarray | None][source]

Assemble the semiclassical wavefunction for the wave (hyperbolic) equation ∂²u/∂t² = ψOp u on a regular grid.

Two-branch structure

The dispersion relation ω² = H factors into two branches H₊ = +√H and H₋ = −√H, each generating its own family of classical rays. The wavefunction is the coherent superposition:

u(x) = Σ_{k∈H₊} A_k exp(i S_k⁺/ℏ − i μ_k⁺ π/2) + Σ_{k∈H₋} A_k exp(i S_k⁻/ℏ − i μ_k⁻ π/2)

where S_k^± = ∫₀ᵗ p · ẋ dt′ along the respective branch rays. Each branch is summed via van_vleck_sum(), then the two grids are added.

Caustic corrections (Airy / Pearcey) are applied independently on each branch before superposition.

param pts_plus:: Scattered ray data for the H₊ = +√H branch.
param S_plus:: Scattered ray data for the H₊ = +√H branch.
param det_J_plus:: Scattered ray data for the H₊ = +√H branch.
param mu_plus:: Scattered ray data for the H₊ = +√H branch.
param pts_minus:: Scattered ray data for the H₋ = −√H branch.
param S_minus:: Scattered ray data for the H₋ = −√H branch.
param det_J_minus:: Scattered ray data for the H₋ = −√H branch.
param mu_minus:: Scattered ray data for the H₋ = −√H branch.
param xlim:: Grid and numerical parameters (same as van_vleck_sum()).
param ylim:: Grid and numerical parameters (same as van_vleck_sum()).
param N:: Grid and numerical parameters (same as van_vleck_sum()).
param hbar:: Grid and numerical parameters (same as van_vleck_sum()).
param reg:: Grid and numerical parameters (same as van_vleck_sum()).
param method:: Grid and numerical parameters (same as van_vleck_sum()).
param caustic_threshold:: Grid and numerical parameters (same as van_vleck_sum()).
returns:: u, X, Y
rtype:: same layout as van_vleck_sum().