geometry

geometry.py — Geometric visualisation of pseudodifferential symbols in 1D and 2D

Overview

The geometry module provides a comprehensive framework for the geometric and semiclassical analysis of pseudodifferential operator symbols H(x,ξ) (1D) or H(x,y,ξ,η) (2D). It combines symbolic differentiation (SymPy) with numerical integration (SciPy) to compute:

• Hamiltonian flows (geodesics) including variational equations for Jacobian matrices.

• Caustics – points where rays focus – detected via zero crossings of the Jacobian determinant.

• Periodic orbits at fixed energies, with stability analysis (Lyapunov exponents).

• Maslov indices and phase shifts associated with caustic crossings.

• Gutzwiller trace formula and semiclassical spectra (1D).

• Weyl’s law, Berry‑Tabor level spacing, and KAM tori classification (2D).

The module is designed for both pedagogical exploration (via richly annotated multi‑panel figures) and research (extracting quantitative data such as caustic networks, Poincaré sections, and action‑angle variables).

Key features:

• Automatic dimension detection – a single uniform interface for 1D and 2D problems.

• Exact symbolic derivatives of the Hamiltonian, lambdified for fast numerical evaluation.

• Augmented ODE systems that simultaneously integrate the Hamiltonian flow and the full 4×4 Jacobian (in 2D), enabling precise caustic detection.

• Periodic orbit search with energy‑preserving initial guesses, deduplication, and stability computation.

• Caustic classification (fold, cusp) using curvature analysis and hierarchical clustering.

• Rich visualisation atlas – 15‑panel figures for 1D, dynamically arranged panels for 2D, covering:

  • Hamiltonian surface, level sets, vector field.

  • Group velocity, spatial projections, Jacobian evolution.

  • Phase space portraits, Poincaré sections, momentum space.

  • Periodic orbits, action‑energy diagrams, EBK quantisation.

  • Gutzwiller trace, semiclassical spectrum, level spacing distributions.

  • Caustic curves, Maslov phase shifts, phase space volume (Monte Carlo).

• Spectral analysis utilities – Weyl’s law, integrability classification, Berry‑Tabor smoothed density.

• Theoretical summaries printed on demand (print_theory(), print_theory_summary()).

Mathematical background

The module studies a Hamiltonian H on phase space T*ℝⁿ (n = 1,2). The associated Hamiltonian vector field is

X_H = ( ∂H/∂p , –∂H/∂x ) (1D) or ( ∂H/∂ξ , ∂H/∂η , –∂H/∂x , –∂H/∂y ) (2D).

Integral curves (x(t),p(t)) are called bicharacteristics or rays. Their projection onto configuration space x(t) gives the physical path.

Caustics occur when nearby rays converge, i.e. when the Jacobian matrix

J_ij = ∂x_i(t)/∂p_j(0)

becomes singular. In 1D J = ∂x/∂p₀; in 2D the 2×2 block ∂(x,y)/∂(ξ₀,η₀) is monitored. Zero crossings of det(J) mark caustics. At a caustic the semiclassical amplitude diverges; the correct uniform approximation involves special functions (Airy for fold, Pearcey for cusp) and a phase shift of –π/2 times the Maslov index (number of caustics traversed).

The Gutzwiller trace formula expresses the quantum density of states as a sum over periodic orbits:

ρ(E) ≈ (1/πℏ) Σ_γ (T_γ/√|det(M_γ – I)|) cos(S_γ/ℏ – πμ_γ/2) .

Here T_γ is the period, S_γ the action, M_γ the monodromy matrix, and μ_γ the Maslov index.

For integrable systems (KAM tori), the Einstein–Brillouin–Keller (EBK) quantisation gives

I_k = (1/2π) ∮ p dq = ℏ (n_k + α_k/4),

with Maslov indices α_k. Level spacings follow a Poisson distribution P(s) = e^{-s}. Chaotic systems obey the Wigner surmise P(s) = (πs/2) e^{-πs²/4} (Bohigas–Giannoni–Schmit conjecture).

References

[1]

Arnold, V. I. Mathematical Methods of Classical Mechanics, Springer‑Verlag, 1989.

[2]

Gutzwiller, M. C. Chaos in Classical and Quantum Mechanics, Springer‑Verlag, 1990.

[3]

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

[4]

Berry, M. V. & Tabor, M. “Level clustering in the regular spectrum”, Proc. R. Soc. Lond. A 356, 375–394, 1977.

[5]

Bohigas, O., Giannoni, M. J., & Schmit, C. “Characterization of chaotic quantum spectra and universality of level fluctuation laws”, Phys. Rev. Lett. 52, 1–4, 1984.

[6]

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

class geometry.CausticStructure(x: ndarray, y: ndarray, t: float, energy: float, type: str, maslov_index: int, strength: float)

Bases: object

Caustic structure with classification (2D only).

energy: float

maslov_index: int

strength: float

t: float

type: str

x: ndarray

y: ndarray

geometry.Geodesic

alias of Geodesic1D

class geometry.Geodesic1D(t: ndarray, H: ndarray, color: str = 'blue', x: ndarray = <factory>, xi: ndarray = <factory>, J: ndarray = <factory>, K: ndarray = <factory>)

Bases: GeodesicBase

1D geodesic in phase space (x, ξ).

J: ndarray

K: ndarray

property caustics: ndarray

Indices where J ≈ 0 (caustic points).

x: ndarray

xi: ndarray

class geometry.Geodesic2D(t: ndarray, H: ndarray, color: str = 'blue', x: ndarray = <factory>, y: ndarray = <factory>, xi: ndarray = <factory>, eta: ndarray = <factory>, J_full: ndarray = <factory>, det_caustic: ndarray = <factory>, caustic_indices: ndarray = <factory>)

Bases: GeodesicBase

2D geodesic in phase space (x, y, ξ, η) with full 4×4 Jacobian.

J_full: ndarray

caustic_indices: ndarray

property caustic_points: Tuple[ndarray, ndarray]

(x, y) coordinates of caustic points along the trajectory.

det_caustic: ndarray

eta: ndarray

property momentum_trajectory: ndarray

Trajectory in momentum space (ξ, η).

property spatial_trajectory: ndarray

Trajectory in configuration space (x, y).

x: ndarray

xi: ndarray

y: ndarray

class geometry.GeodesicBase(t: ndarray, H: ndarray, color: str = 'blue')

Bases: object

Common fields for all geodesics.

H: ndarray

color: str = 'blue'

property energy: float

Initial energy (constant along the geodesic).

t: ndarray

geometry.PeriodicOrbit

alias of PeriodicOrbit1D

class geometry.PeriodicOrbit1D(period: float, action: float, energy: float, x_cycle: ndarray, xi_cycle: ndarray, t_cycle: ndarray, x0: float = 0.0, xi0: float = 0.0, stability: float = 0.0)

Bases: PeriodicOrbitBase

1D periodic orbit in phase space.

stability: float = 0.0

x0: float = 0.0

xi0: float = 0.0

class geometry.PeriodicOrbit2D(period: float, action: float, energy: float, x_cycle: ndarray, xi_cycle: ndarray, t_cycle: ndarray, x0: float = 0.0, y0: float = 0.0, xi0: float = 0.0, eta0: float = 0.0, y_cycle: ndarray = <factory>, eta_cycle: ndarray = <factory>, stability_1: float = 0.0, stability_2: float = 0.0, maslov_index: int = 0)

Bases: PeriodicOrbitBase

2D periodic orbit with Maslov index and KAM stability.

property bohr_sommerfeld_condition: float

Bohr-Sommerfeld quantization condition with Maslov correction.

eta0: float = 0.0

eta_cycle: ndarray

property is_stable: bool

True if both Lyapunov exponents are negative (KAM stable).

maslov_index: int = 0

stability_1: float = 0.0

stability_2: float = 0.0

x0: float = 0.0

xi0: float = 0.0

y0: float = 0.0

y_cycle: ndarray

class geometry.PeriodicOrbitBase(period: float, action: float, energy: float, x_cycle: ndarray, xi_cycle: ndarray, t_cycle: ndarray)

Bases: object

Common fields for periodic orbits.

action: float

energy: float

period: float

t_cycle: ndarray

x_cycle: ndarray

xi_cycle: ndarray

class geometry.SpectralAnalysis

Bases: object

Additional spectral analysis tools (1D).

static analyze_integrability(spacings: ndarray) → Dict

Classify system via level-spacing statistics.

Returns dict with ratio <s²>/<s>², mean spacing, std spacing, and classification string.

static berry_tabor_formula(periodic_orbits: List[PeriodicOrbit1D], energy: float, window: float = 1.0) → float

Berry-Tabor smoothed density of states from periodic orbits.

ρ(E) = Σ_γ T_γ exp(-(E_γ-E)²/2σ²) / (σ√(2π) · 2π)

static weyl_law(energy: float, dimension: int, hbar: float = 1.0) → float

Weyl’s law: asymptotic number of states below energy E.

N(E) ~ (1/2πℏ)^d × Vol{H ≤ E}

Simplified form assumes phase-space volume ~ E^d.

class geometry.Spectrum(energies: ndarray, intensity: ndarray, trace_t: ndarray, trace: ndarray)

Bases: object

Semiclassical spectrum (1D only).

energies: ndarray

intensity: ndarray

trace: ndarray

trace_t: ndarray

class geometry.SymbolGeometry(symbol: Expr, x_sym: Symbol, xi_sym: Symbol)

Bases: SymbolGeometryBase

1D geometric and semiclassical analysis of a symbol H(x, ξ).

Computes: - Hamiltonian flow (geodesics with scalar Jacobi fields J, K) - Caustics (J → 0) - Periodic orbits at fixed energy - Gutzwiller trace formula - Semiclassical spectrum (FFT of trace)

compute_geodesic(x0: float, xi0: float, t_max: float, n_points: int = 500) → Geodesic1D

Integrate the Hamiltonian flow with variational (Jacobi) equations.

System:

dx/dt = ∂H/∂ξ dξ/dt = -∂H/∂x dJ/dt = ∂²H/∂ξ² J + ∂²H/∂x∂ξ K dK/dt = -∂²H/∂x∂ξ J - ∂²H/∂x² K

Initial conditions: J(0)=0, K(0)=1.

find_periodic_orbits(energy: float, x_range: Tuple[float, float], xi_range: Tuple[float, float], n_attempts: int = 50, tol_period: float = 0.001) → List[PeriodicOrbit1D]

gutzwiller_trace_formula(periodic_orbits: List[PeriodicOrbit1D], t_values: ndarray, hbar: float = 1.0) → ndarray

Gutzwiller trace formula: Tr[exp(-iHt/ℏ)] ≈ Σ_γ A_γ exp(iS_γ/ℏ).

semiclassical_spectrum(periodic_orbits: List[PeriodicOrbit1D], hbar: float = 1.0, resolution: int = 4000) → Spectrum

Extract semiclassical spectrum via FFT of the Gutzwiller trace.

class geometry.SymbolGeometry2D(symbol: Expr, x_sym: Symbol, y_sym: Symbol, xi_sym: Symbol, eta_sym: Symbol, hbar: float = 1.0)

Bases: SymbolGeometryBase

2D geometric and semiclassical analysis of H(x, y, ξ, η).

Computes: - Hamiltonian flow with full 4×4 Jacobian (20-dim augmented system) - Caustic detection via sign change of det(∂(x,y)/∂(ξ₀,η₀)) - Caustic classification (fold / cusp) and clustering - Periodic orbits with Maslov index - Monte Carlo phase space volume

compute_geodesic(x0: float, y0: float, xi0: float, eta0: float, t_max: float, n_points: int = 500) → Geodesic2D

Integrate the 2D Hamiltonian flow with full Jacobian evolution.

compute_phase_space_volume(E_max: float, x_range: tuple, y_range: tuple, xi_range: tuple, eta_range: tuple, n_samples: int = 200000) → float

Monte Carlo estimation of Vol{H(x,y,ξ,η) ≤ E_max}.

detect_caustic_structures(geodesics: List[Geodesic2D], t_fixed: float) → List[CausticStructure]

Detect, classify and cluster caustic structures at time t_fixed.

find_periodic_orbits_2d(energy: float, x_range: Tuple[float, float], y_range: Tuple[float, float], xi_range: Tuple[float, float], eta_range: Tuple[float, float], n_attempts: int = 30) → List[PeriodicOrbit2D]

class geometry.SymbolGeometryBase

Bases: ABC

Abstract base class for the geometry engines.

Provides shared helpers: - _remove_duplicate_orbits: generic orbit deduplication - _wrap_solve_ivp: stability computation wrapper

class geometry.SymbolVisualizer(geometry)

Bases: SymbolVisualizerBase

1D symbol visualizer — 15-panel geometric and spectral atlas.

Panels:

1 Hamiltonian surface (3D) 9 Periodic orbits (phase space) 2 Level sets (foliation) 10 Period-energy diagram 3 Hamiltonian vector field 11 EBK quantization 4 Group velocity ∂H/∂ξ 12 Gutzwiller trace formula 5 Spatial projection + caustics 13 Semiclassical spectrum 6 Jacobian J(t) 14 Orbit stability 7 Sectional curvature 15 Level spacing distribution 8 Energy conservation

visualize_complete(x_range: Tuple[float, float], xi_range: Tuple[float, float], geodesics_params: List[Tuple], E_range: Tuple[float, float] | None = None, hbar: float = 1.0, resolution: int = 100) → Tuple

class geometry.SymbolVisualizer2D(geometry)

Bases: SymbolVisualizerBase

2D symbol visualizer — adaptive multi-panel atlas (up to 18 panels).

Panels are assembled dynamically depending on available data (geodesics, periodic orbits, caustics, E_range).

visualize_complete(x_range: Tuple[float, float], y_range: Tuple[float, float], xi_range: Tuple[float, float], eta_range: Tuple[float, float], geodesics_params: List[Tuple], E_range: Tuple[float, float] | None = None, hbar: float = 1.0, resolution: int = 50) → Tuple

class geometry.SymbolVisualizerBase(geometry)

Bases: ABC

Abstract base for 1D and 2D visualizers.

abstractmethod visualize_complete(**kwargs) → Tuple

class geometry.Utilities2D

Bases: object

Additional analysis tools for 2D systems.

static compute_rotation_numbers(geo: Geodesic2D) → Tuple[float, float]

Return rotation numbers (ω_x, ω_y) for the geodesic.

static compute_winding_number(geo: Geodesic2D) → float

Winding number around the origin in configuration space.

static detect_kam_tori(periodic_orbits: List[PeriodicOrbit2D], tolerance: float = 0.1) → Dict

Cluster periodic orbits into KAM tori by action proximity.

geometry.print_theory(dim: int = 1) → None

Print theoretical background for 1D (dim=1, default) or 2D (dim=2).

geometry.print_theory_summary() → None

Alias for print_theory(dim=2) — kept for backward compatibility.

geometry.visualize_symbol(symbol: Expr, x_range: Tuple[float, float], xi_range: Tuple[float, float], geodesics_params: List[Tuple], E_range: Tuple[float, float] | None = None, hbar: float = 1.0, resolution: int = 100, x_sym: Symbol | None = None, xi_sym: Symbol | None = None) → Tuple

1D entry point: complete visualization of H(x, ξ).

Parameters:

• symbol (sympy expression)

• x_range ((min, max))

• xi_range ((min, max))

• geodesics_params (list of (x0, xi0, t_max) or (x0, xi0, t_max, color))

• E_range (optional (E_min, E_max) for spectral analysis)

• hbar (reduced Planck constant)

• resolution (grid resolution)

• x_sym (sympy symbols (default: 'x', 'xi'))

• xi_sym (sympy symbols (default: 'x', 'xi'))

Return type:

fig, geodesics, periodic_orbits, spectrum

geometry.visualize_symbol_2d(symbol: Expr, x_range: Tuple[float, float], y_range: Tuple[float, float], xi_range: Tuple[float, float], eta_range: Tuple[float, float], geodesics_params: List[Tuple], E_range: Tuple[float, float] | None = None, hbar: float = 1.0, resolution: int = 50, x_sym: Symbol | None = None, y_sym: Symbol | None = None, xi_sym: Symbol | None = None, eta_sym: Symbol | None = None) → Tuple

2D entry point: complete visualization of H(x, y, ξ, η).

Parameters:

• symbol (sympy expression)

• x_range ((min, max))

• y_range ((min, max))

• xi_range ((min, max))

• eta_range ((min, max))

• geodesics_params (list of (x0,y0,xi0,eta0,t_max) or (..., color))

• E_range (optional (E_min, E_max))

• hbar (as for 1D)

• resolution (as for 1D)

• x_sym (sympy symbols (defaults: x,y,xi,eta))

• y_sym (sympy symbols (defaults: x,y,xi,eta))

• xi_sym (sympy symbols (defaults: x,y,xi,eta))

• eta_sym (sympy symbols (defaults: x,y,xi,eta))

Return type:

fig, geodesics, periodic_orbits, caustics

Example

>>> x, y   = sp.symbols('x y',   real=True)
>>> xi, eta = sp.symbols('xi eta', real=True)
>>> H = xi**2 + eta**2 + x**2 + y**2
>>> fig, geos, orbits, caustics = visualize_symbol_2d(
...     H, (-2,2), (-2,2), (-2,2), (-2,2),
...     [(1,0,0,1,2*np.pi,'red')], E_range=(0.5,4), hbar=0.1)
>>> plt.show()