physicskit.semiclassical<a class="headerlink" href="#physicskit-semiclassical" title="Link to this heading">#

Returns:

ndarray, shape matching k_grid – \(\sigma_{\min}[1-T(k)]\).

Examples

The lowest Dirichlet level of the unit disk is the first zero of \(J_0\), \(k=2.4048\):

>>> import numpy as np
>>> N = 120
>>> phi = 2 * np.pi * np.arange(N) / N
>>> pts = np.column_stack([np.cos(phi), np.sin(phi)])
>>> ks = np.linspace(2.2, 2.6, 81)
>>> f = bogomolny_quantization_function(ks, pts, -pts, np.full(N, 2 * np.pi / N))
>>> round(float(ks[np.argmin(f)]), 1)
2.4

physicskit.semiclassical.bogomolny_transfer_operator(points, normals, ds, k, boundary='dirichlet')[source]#

Discretized Bogomolny transfer operator \(T(k)\) on the boundary of a convex billiard.

\[T_{ij} = \sqrt{\frac{k\,\cos\theta_i\cos\theta_j}{2\pi i\,\ell_{ij}}}\; e^{ik\ell_{ij} - i\pi\nu/2}\,\sqrt{\Delta s_i\,\Delta s_j},\]

with \(\ell_{ij}\) the chord between boundary points \(i\) and \(j\). The symmetric weight \(\sqrt{\Delta s_i\Delta s_j}\) turns the integral kernel into a matrix whose eigenvalues approximate those of the operator. Use at least a few points per wavelength \(2\pi/k\); extra points only resolve the evanescent part of \(T\) better.

Parameters:

points (ndarray) – Boundary points, ordered along the boundary.
normals (ndarray) – Unit inward normals at points.
ds (ndarray) – Arc-length weight of each point.
k (float) – Wavenumber.
boundary (str) – Wall boundary condition (Maslov phase \(\nu=2\) or 0 per bounce).

Return type:

Returns:

ndarray of complex, shape (N, N) – The matrix \(T(k)\).

Raises:

ValueError – If boundary is not recognized.

Notes

Every chord must lie inside the billiard, so the boundary must be convex. Grazing chords (\(\cos\theta\to0\)) carry whispering-gallery states, which this leading-order kernel describes least accurately.

Examples

On a boundary of length \(L\) there are about \(kL/\pi\) propagating directions. \(T\) is close to unitary on them and close to zero on the rest; for the unit circle at \(k=20\) that is about 40 singular values near 1:

>>> import numpy as np
>>> N = 200
>>> phi = 2 * np.pi * np.arange(N) / N
>>> pts = np.column_stack([np.cos(phi), np.sin(phi)])
>>> T = bogomolny_transfer_operator(pts, -pts, np.full(N, 2 * np.pi / N), k=20.0)
>>> T.shape
(200, 200)
>>> sv = np.linalg.svd(T, compute_uv=False)
>>> int(np.sum(sv > 0.5))
41

physicskit.semiclassical.bohr_sommerfeld_energies(V, m, x_min, x_max, n_max, hbar=1.0, E_min=0.001, E_max=100.0)[source]#

Bohr-Sommerfeld (EBK) bound-state energies from the WKB quantization condition.

Solves \(\int_{x_1(E_n)}^{x_2(E_n)}p(x)\,dx=(n+\tfrac12)\pi\hbar\) for each \(n=0,\dots,n_{max}-1\) by root-finding on \(E\), re-locating the turning points at every trial energy. For the harmonic oscillator this reproduces the exact spectrum \(E_n=(n+\tfrac12)\hbar\omega\) to machine precision, since the WKB approximation is exact whenever the potential is exactly quadratic.

Parameters:

V (callable) – Potential energy function V(x), with a single classically allowed region containing its minimum for every energy searched.
m (float) – Particle mass.
x_min (float) – Search domain for turning_points(); must extend well into the classically forbidden region at E_max.
x_max (float) – Search domain for turning_points(); must extend well into the classically forbidden region at E_max.
n_max (int) – Number of levels to compute (\(n=0,\dots,n_{max}-1\)).
hbar (float) – Value of \(\hbar\) to use.
E_min (float) – Lower bracket for the ground-state root search; must be small enough that its turning points still resolve on the turning_points() search grid.
E_max (float) – Upper bracket for the highest level’s root search.

Return type:

Returns:

ndarray, shape (n_max,) – Quantized energies, ascending.

See also

wkb_action: The action integral this quantizes.
wkb_wavefunction: The corresponding approximate wavefunctions.

Examples

>>> import numpy as np
>>> V = lambda x: 0.5 * x ** 2
>>> energies = bohr_sommerfeld_energies(V, m=1.0, x_min=-20, x_max=20, n_max=5)
>>> np.round(energies, 6)
array([0.5, 1.5, 2.5, 3.5, 4.5])

physicskit.semiclassical.bouncing_ball_energies(R, hbar=1.0, m=1.0, n_max=6)[source]#

Predicted energies of the “bouncing ball” orbit family in a stadium billiard of cap radius R.

A trajectory launched perpendicular to the flat top and bottom walls of a Bunimovich stadium (see physicskit.quantum.chapters.potentials.StadiumBilliard2D) bounces straight up and down forever, blind to the length L of the central rectangle – exactly the motion of a particle in a 1D infinite square well of width \(2R\). Quantizing that 1D motion gives a leading-order prediction for where “bouncing ball” states (eigenstates concentrated on this orbit family) appear in the full 2D spectrum:

\[E_n = \frac{(n\pi\hbar)^2}{2m(2R)^2}, \qquad n=1,2,3,\dots\]

Parameters:

R (float) – Radius of the stadium’s semicircular end-caps (half-height of the billiard).
hbar (float) – Value of \(\hbar\) to use.
m (float) – Particle mass.
n_max (int) – Number of levels to compute.

Return type:

Returns:

ndarray, shape (n_max,) – Predicted bouncing-ball energies, ascending.

See also

bouncing_ball_orbit_points: The classical orbit these energies quantize.

Examples

This is exactly the ordinary infinite-square-well spectrum for a well of width \(2R\):

>>> import numpy as np
>>> R = 0.5
>>> energies = bouncing_ball_energies(R, n_max=3)
>>> exact = (np.arange(1, 4) * np.pi) ** 2 / (2 * (2 * R) ** 2)
>>> bool(np.allclose(energies, exact))
True

physicskit.semiclassical.bouncing_ball_orbit_points(x0, R, n_bounces=4)[source]#

Trace the classical “bouncing ball” orbit at horizontal position x0 in a stadium billiard.

A perpendicular launch from the bottom wall at \((x_0,-R)\) reflects straight back and forth off the flat top (\(y=R\)) and bottom (\(y=-R\)) walls, tracing the same vertical segment forever – returned here as a zig-zag list of points suitable for overlaying on a density plot.

Parameters:

x0 (float) – Horizontal position of the orbit, with \(|x_0| < L/2\) (strictly inside the flat central section of the stadium).
R (float) – Radius of the stadium’s semicircular end-caps.
n_bounces (int) – Number of full up-down traversals to trace.

Return type:

Returns:

x, y (ndarray) – Coordinates tracing the orbit, shape (2 * n_bounces + 1,).

See also

bouncing_ball_energies: Quantized energies of this orbit family.

Examples

>>> x, y = bouncing_ball_orbit_points(x0=0.2, R=0.5, n_bounces=2)
>>> x
array([0.2, 0.2, 0.2, 0.2, 0.2])
>>> y
array([-0.5,  0.5, -0.5,  0.5, -0.5])

physicskit.semiclassical.classical_momentum(E, V, x, m=1.0)[source]#

Classical momentum \(p(x) = \sqrt{2m(E-V(x))}\) in the allowed region.

Parameters:

E (float) – Total energy.
V (callable) – Potential energy function V(x).
x (ndarray) – Positions at which to evaluate \(p(x)\).
m (float) – Particle mass.

Return type:

Returns:

ndarray – \(p(x)\), same shape as x; nan wherever \(V(x)>E\) (the classically forbidden region, where WKB oscillates into a real exponential instead).

Examples

>>> import numpy as np
>>> V = lambda x: 0.5 * x ** 2
>>> p = classical_momentum(E=2.0, V=V, x=np.array([0.0, 1.0, 3.0]))
>>> np.round(p, 4)
array([2.    , 1.7321,    nan])

physicskit.semiclassical.classical_period(E, V, m, x_min, x_max, hbar=1.0, dE=0.0001)[source]#

Classical (round-trip) oscillation period \(T(E)=2\,dS/dE\) of a 1D bound orbit.

The action \(S(E)\) (wkb_action()) is the one-way traversal action, so its energy derivative \(dS/dE\) is the time to cross the classically allowed region once; the full period (there and back) is twice that – differentiated numerically here via a central difference.

Parameters:

E (float) – Energy at which to evaluate the period.
V (callable) – Potential energy function V(x).
m (float) – Particle mass.
x_min (float) – Search domain for turning_points().
x_max (float) – Search domain for turning_points().
hbar (float) – Value of \(\hbar\) to use (unused by the classical period itself; kept for a uniform signature with the rest of this module).
dE (float) – Step size for the central-difference derivative.

Return type:

Returns:

float – The classical period \(T(E)\).

See also

gutzwiller_density_of_states: Uses this as the trace formula’s smooth prefactor.

Examples

The harmonic oscillator’s period is famously independent of energy, \(T=2\pi/\omega\):

>>> V = lambda x: 0.5 * x ** 2
>>> round(classical_period(E=3.0, V=V, m=1.0, x_min=-20, x_max=20), 4)
6.2832

physicskit.semiclassical.coherent_state_overlap(q1, p1, q2, p2, gamma, hbar=1.0)[source]#

Closed-form overlap \(\langle g_{q_1,p_1}|g_{q_2,p_2}\rangle\) of two equal-width frozen Gaussians.

\[\langle g_1|g_2\rangle = \exp\!\left[-\frac{\gamma}{2}(q_1-q_2)^2 - \frac{(p_1-p_2)^2}{8\gamma\hbar^2} + \frac{i}{2\hbar}(q_1-q_2)(p_1+p_2)\right].\]

Parameters:

q1 (float) – Phase-space center of the bra state.
p1 (float) – Phase-space center of the bra state.
q2 (float) – Phase-space center of the ket state.
p2 (float) – Phase-space center of the ket state.
gamma (float) – Shared width parameter, as in frozen_gaussian_1d().
hbar (float) – Value of \(\hbar\) to use.

Return type:

Returns:

complex – The overlap \(\langle g_1|g_2\rangle\).

See also

frozen_gaussian_1d: The states being overlapped.

Examples

Matches direct numerical integration of the two wavepackets:

>>> import numpy as np
>>> x = np.linspace(-40, 40, 20000)
>>> gamma = 1.0
>>> psi1 = frozen_gaussian_1d(x, qc=0.3, pc=0.7, gamma=gamma)
>>> psi2 = frozen_gaussian_1d(x, qc=-0.2, pc=1.1, gamma=gamma)
>>> numeric = np.trapezoid(np.conj(psi1) * psi2, x)
>>> closed_form = coherent_state_overlap(0.3, 0.7, -0.2, 1.1, gamma)
>>> bool(abs(numeric - closed_form) < 1e-6)
True

A state’s overlap with itself is 1:

>>> round(float(abs(coherent_state_overlap(1.0, 2.0, 1.0, 2.0, gamma=0.8))), 8)
1.0

physicskit.semiclassical.count_caustics(Mqp_history)[source]#

Count sign changes of \(\partial q_t/\partial p_0\) along a trajectory – the Maslov index.

Each time \(\partial q_t/\partial p_0\) passes through zero, the trajectory crosses a focal point (conjugate point / caustic), and the semiclassical propagator picks up an extra phase of \(-\pi/2\) (Gutzwiller 1967; see also the \(-\pi/4\) phase at a single WKB turning point in physicskit.semiclassical.core.wkb.wkb_wavefunction(), which is this same phenomenon in the time-independent picture).

Parameters:: Mqp_history (ndarray) – \(\partial q_t/\partial p_0\) sampled along the trajectory, from propagate_trajectory_monodromy_action(); the first entry (always 0, at \(t=0\)) is ignored.
Return type:: int
Returns:: int – The Maslov index \(\mu\).

See also

van_vleck_propagator_1d: Uses this to fix the propagator’s overall phase.

Examples

>>> count_caustics(np.array([0.0, 1.0, 0.5, -0.3, -0.8, 0.2]))
2

physicskit.semiclassical.frozen_gaussian_1d(x, qc, pc, gamma, hbar=1.0)[source]#

A normalized, fixed-width (“frozen”) Gaussian wavepacket (coherent state).

\[g_{q_c,p_c}(x) = \left(\frac{2\gamma}{\pi}\right)^{1/4} \exp\!\left[-\gamma(x-q_c)^2 + \frac{i}{\hbar}p_c(x-q_c)\right].\]

The building block of the Herman-Kluk method: every trajectory carries one of these, always at the fixed width set by \(\gamma\) (hence “frozen”), riding on top of its classical phase-space point \((q_c,p_c)\).

Parameters:

x (ndarray) – Positions at which to evaluate the wavepacket.
qc (float) – Center of the Gaussian in position and momentum.
pc (float) – Center of the Gaussian in position and momentum.
gamma (float) – Width parameter (inverse squared length); larger \(\gamma\) means a narrower, more position-localized packet.
hbar (float) – Value of \(\hbar\) to use.

Return type:

Returns:

ndarray of complex – \(g_{q_c,p_c}(x)\), same shape as x.

See also

coherent_state_overlap: The closed-form overlap of two such Gaussians.

Examples

>>> import numpy as np
>>> x = np.linspace(-20, 20, 4000)
>>> psi = frozen_gaussian_1d(x, qc=1.0, pc=2.0, gamma=0.5)
>>> round(float(np.trapezoid(np.abs(psi) ** 2, x)), 8)
1.0

physicskit.semiclassical.gutzwiller_amplitude_from_monodromy(M)[source]#

Gutzwiller stability amplitude \(1/\sqrt{|2-\operatorname{tr}M|}\) for an isolated periodic orbit.

In the full (multi-dimensional, generically chaotic) Gutzwiller trace formula, each isolated periodic orbit contributes with an amplitude set by how strongly nearby trajectories diverge from it over one period: an unstable (hyperbolic) orbit with monodromy eigenvalues \(\lambda,1/\lambda\) (\(|\lambda|>1\)) has \(\operatorname{tr}M=\lambda+1/\lambda\), so \(|2-\operatorname{tr}M|\) grows with the instability and the orbit’s contribution to the trace formula shrinks – highly unstable orbits matter less for the exact spectrum, but (as in physicskit.semiclassical.systems.scarring) can still leave a visible imprint on individual eigenstates.

Parameters:: M (ndarray) – Monodromy matrix of one period of the orbit, e.g. from physicskit.semiclassical.core.propagators.propagate_trajectory_monodromy_action().
Return type:: float
Returns:: float – The stability amplitude.

Examples

>>> import numpy as np
>>> M = np.array([[2.0, 0.0], [0.0, 0.5]])
>>> round(gutzwiller_amplitude_from_monodromy(M), 6)
1.414214

physicskit.semiclassical.gutzwiller_density_of_states(E_grid, V, m, x_min, x_max, hbar=1.0, r_max=40, broadening=0.05)[source]#

Exact 1D trace-formula density of states, summed over repetitions of the single periodic orbit.

For a bound one-dimensional system, EBK quantization places levels exactly where \(S(E_n)/\hbar=(n+\tfrac12)\pi\). Poisson-summing the resulting delta comb \(g(E)=\sum_n\delta(E-E_n)\) over the integer \(n\) converts it into a sum over an integer \(r\) – physically, the \(r\)-fold repetition of the single primitive periodic orbit at each energy – giving the exact identity

\[g(E) = \frac{T(E)/2}{\pi\hbar}\left[1 + 2\sum_{r=1}^{\infty} \cos\!\left(\frac{2rS(E)}{\hbar} - r\pi\right)\right],\]

the one-dimensional Gutzwiller trace formula: the \(r=0\) term is the smooth Weyl density of states, and each \(r\ge1\) term is one repetition of the orbit, carrying the Maslov phase \(-r\pi\) (\(\sigma=2\) soft turning points per traversal, repeated \(r\) times). Truncating the sum at finite r_max with a convergence factor exp(-r*broadening) turns each delta function into a finite (Lorentzian-like) peak, suitable for numerical peak-finding.

Parameters:

E_grid (ndarray) – Energies at which to evaluate the density of states.
V (callable) – Potential energy function V(x).
m (float) – Particle mass.
x_min (float) – Search domain for turning_points().
x_max (float) – Search domain for turning_points().
hbar (float) – Value of \(\hbar\) to use.
r_max (int) – Number of orbit repetitions to sum.
broadening (float) – Per-repetition convergence factor; larger values broaden (and damp) each reconstructed peak.

Return type:

Returns:

ndarray – \(g(E)\), same shape as E_grid.

See also

classical_period: Supplies \(T(E)\).
physicskit.semiclassical.core.wkb.bohr_sommerfeld_energies: The exact peak locations this reconstructs.

Examples

The reconstructed peaks land exactly on the harmonic oscillator’s Bohr-Sommerfeld spectrum:

>>> import numpy as np
>>> from scipy.signal import find_peaks
>>> V = lambda x: 0.5 * x ** 2
>>> E_grid = np.linspace(0.2, 4.5, 600)
>>> dos = gutzwiller_density_of_states(E_grid, V, m=1.0, x_min=-20, x_max=20)
>>> peak_idx, _ = find_peaks(dos, height=0.3 * dos.max())
>>> np.round(E_grid[peak_idx], 1)
array([0.5, 1.5, 2.5, 3.5])

physicskit.semiclassical.herman_kluk_prefactor(M, gamma, hbar=1.0)[source]#

Herman-Kluk prefactor \(C_t(q_0,p_0)\) built from the monodromy matrix.

\[C_t = \sqrt{\frac{1}{2}\left(M_{qq} + M_{pp} - 2i\hbar\gamma M_{qp} + \frac{i}{2\hbar\gamma}M_{pq}\right)},\]

the standard Herman-Kluk form (Herman & Kluk 1984; Kay 1994) written for this module’s frozen Gaussians \(e^{-\gamma(x-q)^2}\), whose width parameter in the usual \(e^{-\gamma_s(x-q)^2/2}\) convention is \(\gamma_s=2\gamma\). It reduces to \(C_t=1\) at \(t=0\) (\(M=\mathbb{1}\)), consistent with the frozen Gaussian basis’s resolution of the identity, \(\int \tfrac{dq\,dp}{2\pi\hbar}\,|g_{q,p}\rangle\langle g_{q,p}| = \hat{1}\).

Parameters:

M (ndarray) – Monodromy matrix from propagate_trajectory_monodromy_action().
gamma (float) – Frozen-Gaussian width parameter, as in frozen_gaussian_1d().
hbar (float) – Value of \(\hbar\) to use.

Return type:

Returns:

complex – The prefactor \(C_t\).

Examples

>>> import numpy as np
>>> complex(herman_kluk_prefactor(np.eye(2), gamma=1.0))
(1+0j)

physicskit.semiclassical.herman_kluk_propagate_wavepacket(qc0, pc0, gamma, dVdx, d2Vdx2, V, m, dt, steps, x_eval, hbar=1.0, n_grid=41, n_sigma=6.0, params=None)[source]#

Propagate an initial frozen-Gaussian wavepacket with the multi-trajectory Herman-Kluk method.

Launches one classical trajectory from every point of a regular grid covering the initial coherent state’s phase-space support, and sums their frozen-Gaussian contributions

\[\psi(x,t) \approx \int\!\frac{dq_0\,dp_0}{2\pi\hbar}\, C_t(q_0,p_0)\,e^{iS(t)/\hbar}\, \langle g_{q_0,p_0}|\psi_0\rangle\, g_{q_t,p_t}(x),\]

approximating the integral over initial conditions by a Riemann sum on a grid spanning \(\pm n_\sigma\) standard deviations of the initial coherent state’s own phase-space Gaussian (\(\sigma_q=1/(2\sqrt\gamma)\), \(\sigma_p=\hbar\sqrt\gamma\)). The per-trajectory classical propagation (propagate_trajectory_monodromy_action()) is Numba-compiled, since it is called once per grid point – the dominant cost for any reasonably fine grid.

In the \(t\to0\) limit this reduces to the frozen-Gaussian resolution of the identity and reconstructs the initial wavepacket essentially exactly; away from that limit, this quadrature-grid evaluation of the Herman-Kluk integral does not exactly conserve \(\int|\psi|^2\,dx\) (the true continuous phase-space integral does, for at-most-quadratic potentials) – convergence in n_grid, n_sigma, and gamma should be checked for any serious use.

Parameters:

qc0 (float) – Center of the initial frozen-Gaussian wavepacket.
pc0 (float) – Center of the initial frozen-Gaussian wavepacket.
gamma (float) – Width parameter shared by the initial state and every frozen Gaussian in the propagation.
dVdx (callable) – Numba-jitted potential derivatives and the potential itself, as in propagate_trajectory_monodromy_action().
d2Vdx2 (callable) – Numba-jitted potential derivatives and the potential itself, as in propagate_trajectory_monodromy_action().
V (callable) – Numba-jitted potential derivatives and the potential itself, as in propagate_trajectory_monodromy_action().
m (float) – Particle mass.
dt (float) – Time step for each trajectory’s RK4 integration.
steps (int) – Number of RK4 steps (propagation time is dt * steps).
x_eval (ndarray) – Positions at which to evaluate the propagated wavepacket.
hbar (float) – Value of \(\hbar\) to use.
n_grid (int) – Number of grid points along each of the \(q_0\), \(p_0\) axes.
n_sigma (float) – Half-width of the sampling grid, in standard deviations of the initial state’s phase-space Gaussian.
params (ndarray | None) – Parameter vector passed through to dVdx/d2Vdx2/V.

Return type:

Returns:

ndarray of complex, shape matching x_eval – The propagated wavepacket \(\psi(x,t)\).

See also

van_vleck_propagator_1d: The single-trajectory propagator this sums many copies of.
frozen_gaussian_1d: The initial state and the basis each trajectory carries.

Examples

At very short times the sum reduces to the frozen-Gaussian resolution of the identity and reproduces the (unpropagated) initial coherent state to high accuracy:

>>> import numpy as np
>>> from numba import njit
>>> m, omega = 1.0, 1.0
>>> params = np.array([m * omega ** 2])
>>> dVdx = njit(lambda q, params: params[0] * q, cache=False)
>>> d2Vdx2 = njit(lambda q, params: params[0], cache=False)
>>> V = njit(lambda q, params: 0.5 * params[0] * q ** 2, cache=False)
>>> x = np.linspace(-6, 6, 400)
>>> psi = herman_kluk_propagate_wavepacket(
...     1.0, 0.0, gamma=1.0, dVdx=dVdx, d2Vdx2=d2Vdx2, V=V, m=m, dt=1e-6, steps=1, x_eval=x, n_grid=61, n_sigma=7.0, params=params
... )
>>> psi0 = frozen_gaussian_1d(x, qc=1.0, pc=0.0, gamma=1.0)
>>> fidelity = abs(np.trapezoid(np.conj(psi0) * psi, x)) ** 2
>>> bool(fidelity > 0.999)
True

physicskit.semiclassical.husimi_projection_1d(psi, s, hbar=1.0, sigma=None, resolution=60, s0_range=None, p0_range=None)[source]#

Husimi (coherent-state) phase-space projection of a 1D wavefunction slice on an open interval.

The non-periodic counterpart of physicskit.chaos.quantum.husimi.husimi_function(): overlaps psi with ordinary (non-periodized) coherent states \(g_{s_0,p_0}(s)=\exp[-(s-s_0)^2/2\sigma^2+ip_0(s-s_0)/\hbar]\) at every point of an \((s_0,p_0)\) grid, suited to a slice taken along an open boundary – e.g. the flat top wall of a stadium billiard, where bouncing_ball_orbit_points() lives – rather than a periodic domain.

Parameters:

psi (ndarray) – Wavefunction values sampled on s.
s (ndarray) – Uniform 1D coordinate grid psi is sampled on.
hbar (float) – Value of \(\hbar\) to use.
sigma (float | None) – Coherent-state position width. Defaults to \(\sqrt{\hbar}\) (the minimum-uncertainty, equal-spread-in-natural-units choice).
resolution (int) – Number of grid points along each of the \(s_0\), \(p_0\) axes.
s0_range (tuple | None) – Range of \(s_0\) to scan. Defaults to (s.min(), s.max()).
p0_range (tuple | None) – Range of \(p_0\) to scan. Defaults to \(\pm\pi\hbar/\Delta s\) (the grid’s Nyquist momentum).

Return type:

Returns:

S0, P0 (ndarray, shape (resolution, resolution)) – Phase-space grid (indexing="ij").
husimi (ndarray, shape (resolution, resolution)) – The Husimi distribution, normalized to a peak value of 1.

See also

scar_enhancement: A simpler, single-number scarring diagnostic in position space rather than phase space.

Examples

A Gaussian wave packet’s Husimi projection peaks at its own position and momentum:

>>> import numpy as np
>>> s = np.linspace(-10, 10, 2000)
>>> s0_true, p0_true, w = 2.0, 3.0, 1.0
>>> psi = np.exp(-(s - s0_true) ** 2 / (2 * w ** 2)) * np.exp(1j * p0_true * s / 1.0)
>>> S0, P0, H = husimi_projection_1d(psi, s, hbar=1.0, sigma=w, resolution=80, s0_range=(-2, 6), p0_range=(-2, 8))
>>> i, j = np.unravel_index(np.argmax(H), H.shape)
>>> bool(abs(S0[i, j] - s0_true) < 0.2 and abs(P0[i, j] - p0_true) < 0.2)
True

physicskit.semiclassical.langer_corrected_wkb(V, l, m, r_max, n_max, hbar=1.0, langer=True, E_max=None, r_min=1e-08, n_search=4000)[source]#

Radial WKB bound-state energies, with or without the Langer correction.

The radial Schrodinger equation for \(u(r)=rR(r)\) is a 1D problem on \(r>0\) with effective potential \(V(r)+\hbar^2 l(l+1)/(2mr^2)\). Applying the WKB rule \(\int_{r_1}^{r_2}p\,dr=(n_r+\tfrac12)\pi\hbar\) to it directly gives wrong energies, because the WKB phase near \(r=0\) is not the phase of the true solution \(u\sim r^{l+1}\). Langer (1937) showed that the substitution \(r=e^x\) turns the radial problem into one on the whole line, where WKB applies cleanly, and that the result is the same as using

\[V_{\rm eff}(r) = V(r) + \frac{\hbar^2(l+\tfrac12)^2}{2mr^2}\]

in the ordinary WKB rule. With this correction WKB gives the exact hydrogen levels \(E=-mk^2/(2\hbar^2(n_r+l+1)^2)\) for \(V=-k/r\), and the exact 3D harmonic-oscillator levels \(\hbar\omega(2n_r+l+\tfrac32)\).

Parameters:

V (callable) – Radial potential V(r), vectorized over a NumPy array.
l (int) – Orbital angular momentum quantum number.
m (float) – Particle mass.
r_max (float) – Outer edge of the radial domain; every outer turning point must lie inside it.
n_max (int) – Number of radial levels to compute (\(n_r=0,\dots,n_{max}-1\)).
hbar (float) – Value of \(\hbar\) to use.
langer (bool) – Use \((l+\tfrac12)^2\); if False, use the uncorrected \(l(l+1)\).
E_max (float | None) – Upper bracket for the root search. Defaults to \(V_{\rm eff}(r_{max})\), the highest energy whose outer turning point is still inside the domain.
r_min (float) – Inner edge of the radial domain. If \(V_{\rm eff}\) has no inner turning point (\(l=0\) without the Langer correction, for a Coulomb potential), the action integral starts here.
n_search (int) – Number of log-spaced grid points used to bracket turning points.

Return type:

Returns:

ndarray, shape (n_max,) – Energies for \(n_r=0,\dots,n_{max}-1\), ascending.

Raises:

ValueError – If a level lies above E_max (increase r_max).

See also

bohr_sommerfeld_energies: The same rule for a potential on the whole line.

References

R. E. Langer, “On the Connection Formulas and the Solutions of the Wave Equation,” Phys. Rev. 51, 669-676 (1937).

Examples

With the correction, WKB reproduces the hydrogen 2p, 3p, 4p levels \(-1/(2n^2)\) exactly; without it, it does not:

>>> import numpy as np
>>> V = lambda r: -1.0 / r
>>> np.round(langer_corrected_wkb(V, l=1, m=1.0, r_max=200.0, n_max=3), 6)
array([-0.125   , -0.055556, -0.03125 ])
>>> np.round(langer_corrected_wkb(V, l=1, m=1.0, r_max=200.0, n_max=3, langer=False), 4)
array([-0.1365, -0.0589, -0.0326])

physicskit.semiclassical.plot_classical_trajectory_on_wigner(x, psi, q_hist, p_hist, hbar=1.0, ax=None)[source]#

Overlay a classical phase-space trajectory on the exact quantum Wigner function.

Compares the classical and quantum pictures directly: the Wigner function of a quantum state, with the classical trajectory (physicskit.semiclassical.core.propagators.propagate_trajectory_monodromy_action()) that WKB/Van Vleck theory builds that state’s semiclassical approximation from, drawn on top – the two agree closely wherever the semiclassical approximation is accurate.

Parameters:

x (ndarray) – Position grid the wavefunction is sampled on.
psi (ndarray) – Quantum wavefunction sampled on x.
q_hist (ndarray) – Classical trajectory phase-space coordinates.
p_hist (ndarray) – Classical trajectory phase-space coordinates.
hbar (float) – Value of \(\hbar\) to use.
ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes)

Examples

>>> import numpy as np
>>> from physicskit.semiclassical.core.propagators import frozen_gaussian_1d
>>> x = np.linspace(-8, 8, 400)
>>> psi = frozen_gaussian_1d(x, qc=1.0, pc=0.5, gamma=1.0)
>>> theta = np.linspace(0, 2 * np.pi, 100)
>>> q_hist, p_hist = np.cos(theta), np.sin(theta)
>>> fig, ax = plot_classical_trajectory_on_wigner(x, psi, q_hist, p_hist)
>>> isinstance(fig, plt.Figure)
True

physicskit.semiclassical.plot_density_of_states(E_grid, dos, exact_energies=None, ax=None)[source]#

Plot the Gutzwiller trace-formula density of states, with exact levels marked.

Parameters:

E_grid (ndarray) – Energy grid.
dos (ndarray) – Density of states, e.g. from physicskit.semiclassical.core.gutzwiller.gutzwiller_density_of_states().
exact_energies (ndarray | None) – Exact (or Bohr-Sommerfeld) energy levels to mark with vertical dashed lines.
ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes)

Examples

>>> import numpy as np
>>> from physicskit.semiclassical.core.gutzwiller import gutzwiller_density_of_states
>>> from physicskit.semiclassical.core.wkb import bohr_sommerfeld_energies
>>> V = lambda x: 0.5 * x ** 2
>>> E_grid = np.linspace(0.2, 3.5, 200)
>>> dos = gutzwiller_density_of_states(E_grid, V, m=1.0, x_min=-20, x_max=20)
>>> energies = bohr_sommerfeld_energies(V, m=1.0, x_min=-20, x_max=20, n_max=3)
>>> fig, ax = plot_density_of_states(E_grid, dos, exact_energies=energies)
>>> isinstance(fig, plt.Figure)
True

physicskit.semiclassical.plot_husimi_1d(S0, P0, husimi, ax=None)[source]#

Contour-plot a 1D Husimi phase-space projection.

Parameters:

S0 (ndarray) – Phase-space grid, from physicskit.semiclassical.systems.scarring.husimi_projection_1d().
P0 (ndarray) – Phase-space grid, from physicskit.semiclassical.systems.scarring.husimi_projection_1d().
husimi (ndarray) – Husimi distribution, same shape as S0.
ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes)

Examples

>>> import numpy as np
>>> from physicskit.semiclassical.systems.scarring import husimi_projection_1d
>>> s = np.linspace(-10, 10, 1000)
>>> psi = np.exp(-(s - 2.0) ** 2 / 2.0) * np.exp(1j * 3.0 * s)
>>> S0, P0, H = husimi_projection_1d(psi, s, sigma=1.0, resolution=40, s0_range=(-2, 6), p0_range=(-2, 8))
>>> fig, ax = plot_husimi_1d(S0, P0, H)
>>> isinstance(fig, plt.Figure)
True

physicskit.semiclassical.plot_scar_map(X, Y, density, mask=None, orbit_x=None, orbit_y=None, ax=None)[source]#

Heatmap an eigenstate density with a classical periodic orbit overlaid.

Parameters:

X (ndarray) – Coordinate meshgrid.
Y (ndarray) – Coordinate meshgrid.
density (ndarray) – Probability density \(|\psi|^2\), same shape as X.
mask (ndarray | None) – True inside the billiard; density outside is not drawn.
orbit_x (ndarray | None) – Classical orbit coordinates, e.g. from physicskit.semiclassical.systems.scarring.bouncing_ball_orbit_points().
orbit_y (ndarray | None) – Classical orbit coordinates, e.g. from physicskit.semiclassical.systems.scarring.bouncing_ball_orbit_points().
ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes)

Examples

>>> import numpy as np
>>> from physicskit.semiclassical.systems.scarring import bouncing_ball_orbit_points
>>> x = np.linspace(-2, 2, 100)
>>> y = np.linspace(-1, 1, 60)
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> density = np.exp(-(X - 0.3) ** 2 / (2 * 0.2 ** 2))
>>> orbit_x, orbit_y = bouncing_ball_orbit_points(x0=0.3, R=1.0)
>>> fig, ax = plot_scar_map(X, Y, density, orbit_x=orbit_x, orbit_y=orbit_y)
>>> isinstance(fig, plt.Figure)
True

physicskit.semiclassical.plot_scar_map_interactive(X, Y, density)[source]#

Interactive Plotly heatmap of an eigenstate density, for zooming into scarred structure.

Parameters:

X (ndarray) – Coordinate meshgrid.
Y (ndarray) – Coordinate meshgrid.
density (ndarray) – Probability density \(|\psi|^2\), same shape as X.

Returns:

plotly.graph_objects.Figure

See also

plot_scar_map: The static Matplotlib heatmap equivalent, with orbit overlay support.

Examples

>>> import numpy as np
>>> x = np.linspace(-2, 2, 60)
>>> y = np.linspace(-1, 1, 40)
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> density = np.exp(-(X - 0.3) ** 2 / (2 * 0.2 ** 2))
>>> fig = plot_scar_map_interactive(X, Y, density)
>>> isinstance(fig, go.Figure)
True

physicskit.semiclassical.plot_wkb_wavefunction(x, psi, V=None, ax=None)[source]#

Plot a WKB (or exact) 1D wavefunction, optionally with the potential on a twin axis.

Parameters:

x (ndarray) – Position grid.
psi (ndarray) – Wavefunction values, e.g. from physicskit.semiclassical.core.wkb.wkb_wavefunction().
V (callable, optional) – Potential energy function V(x), drawn on a secondary y-axis if given.
ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes)

Examples

>>> import numpy as np
>>> from physicskit.semiclassical.core.wkb import bohr_sommerfeld_energies, wkb_wavefunction
>>> V = lambda x: 0.5 * x ** 2
>>> E2 = bohr_sommerfeld_energies(V, m=1.0, x_min=-20, x_max=20, n_max=3)[2]
>>> x = np.linspace(-6, 6, 1000)
>>> psi = wkb_wavefunction(x, E2, V)
>>> fig, ax = plot_wkb_wavefunction(x, psi, V=V)
>>> isinstance(fig, plt.Figure)
True

physicskit.semiclassical.propagate_trajectory_monodromy_action(q0, p0, dVdx, d2Vdx2, V, m, dt, steps, params=None)[source]#

Integrate a classical trajectory together with its monodromy matrix and Hamilton principal function.

Advances the Hamilton equations \(\dot q=p/m,\ \dot p=-V'(q)\) together with the variational (tangent) equations for small deviations, \(\dot{\delta q}=\delta p/m,\ \dot{\delta p} =-V''(q)\delta q\), whose solution operator is the monodromy matrix \(M(t)=\begin{pmatrix}\partial q_t/\partial q_0 & \partial q_t/\partial p_0\\ \partial p_t/\partial q_0 & \partial p_t/\partial p_0\end{pmatrix}\), and accumulates the classical action (Hamilton’s principal function) \(S(t)=\int_0^t\left[\tfrac{p^2}{2m}-V(q)\right]dt'\) along the way – everything van_vleck_propagator_1d() and herman_kluk_propagate_wavepacket() need from one trajectory, computed in a single RK4 pass.

Because the flow is Hamiltonian, \(\det M(t)=1\) exactly for all \(t\) (Liouville’s theorem); this is a good numerical sanity check on any trajectory this function returns.

Parameters:

q0 (float) – Initial position and momentum.
p0 (float) – Initial position and momentum.
dVdx (callable) – Numba-jitted force law dVdx(q, params) -> float, \(V'(q)\).
d2Vdx2 (callable) – Numba-jitted curvature d2Vdx2(q, params) -> float, \(V''(q)\).
V (callable) – Numba-jitted potential V(q, params) -> float.
m (float) – Particle mass.
dt (float) – Time step.
steps (int) – Number of RK4 steps (total propagation time is dt * steps).
params (ndarray | None) – Parameter vector passed through to dVdx/d2Vdx2/V. Defaults to an empty array.

Returns:

q_t, p_t (float) – Final position and momentum.
M (ndarray, shape (2, 2)) – Monodromy matrix at time t = dt * steps.
S (float) – Classical action accumulated along the trajectory.
Mqp_history (ndarray, shape (steps + 1,)) – \(\partial q_t/\partial p_0\) at every step, for count_caustics().

See also

van_vleck_prefactor: Turns M into a propagator amplitude.
count_caustics: Turns Mqp_history into a Maslov index.

Examples

For the harmonic oscillator, the monodromy matrix has the exact closed form \(\begin{pmatrix}\cos\omega t & \sin(\omega t)/(m\omega)\\ -m\omega\sin\omega t & \cos\omega t\end{pmatrix}\):

>>> import numpy as np
>>> from numba import njit
>>> m, omega = 1.0, 1.0
>>> params = np.array([m * omega ** 2])
>>> dVdx = njit(lambda q, params: params[0] * q, cache=False)
>>> d2Vdx2 = njit(lambda q, params: params[0], cache=False)
>>> V = njit(lambda q, params: 0.5 * params[0] * q ** 2, cache=False)
>>> t = 1.3
>>> q_t, p_t, M, S, _ = propagate_trajectory_monodromy_action(1.0, 0.3, dVdx, d2Vdx2, V, m, dt=t / 4000, steps=4000, params=params)
>>> M_exact = np.array([[np.cos(t), np.sin(t)], [-np.sin(t), np.cos(t)]])
>>> bool(np.max(np.abs(M - M_exact)) < 1e-6)
True
>>> round(float(np.linalg.det(M)), 8)
1.0

physicskit.semiclassical.scar_enhancement(density, X, Y, mask, x0, half_width)[source]#

Density enhancement in a tube around a vertical bouncing-ball orbit, relative to the billiard average.

\[\eta = \frac{\langle|\psi|^2\rangle_{\text{tube}}}{\langle|\psi|^2\rangle_{\text{billiard}}}, \qquad \text{tube} = \{(x,y)\in\text{billiard} : |x-x_0|\le w\},\]

a simple, direct measure of scarring: \(\eta\approx1\) for an ergodically-spread (unscarred) state, since the tube then just samples a representative fraction of the whole density, while \(\eta\gg1\) signals density piled up specifically along the orbit at \(x_0\).

Parameters:

density (ndarray) – Probability density \(|\psi(x,y)|^2\) on a grid.
X (ndarray) – Coordinate meshgrid matching density.
Y (ndarray) – Coordinate meshgrid matching density.
mask (ndarray) – True inside the billiard (e.g. from physicskit.quantum.chapters.potentials.StadiumBilliard2D.mask()).
x0 (float) – Horizontal position of the bouncing-ball orbit to test.
half_width (float) – Half-width of the tube around x0.

Return type:

Returns:

float – The enhancement factor \(\eta\).

Examples

A density sharply concentrated in a narrow ridge at x0 is strongly enhanced inside a tube covering that ridge:

>>> import numpy as np
>>> x = np.linspace(-2, 2, 400)
>>> y = np.linspace(-1, 1, 200)
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> mask = np.ones_like(X, dtype=bool)
>>> density = np.exp(-(X - 0.3) ** 2 / (2 * 0.05 ** 2))
>>> eta = scar_enhancement(density, X, Y, mask, x0=0.3, half_width=0.15)
>>> bool(12.0 < eta < 14.0)
True

An (exactly) uniform density is not enhanced anywhere:

>>> uniform = np.ones_like(X)
>>> round(scar_enhancement(uniform, X, Y, mask, x0=0.3, half_width=0.15), 8)
1.0

physicskit.semiclassical.thawed_gaussian_propagate(qc0, pc0, gamma, dVdx, d2Vdx2, V, m, dt, steps, hbar=1.0, params=None)[source]#

Propagate a Gaussian wavepacket with Heller’s thawed Gaussian approximation.

Heller (1975) wrote the packet as

\[\psi(x,t) = \exp\!\left\{\frac{i}{\hbar}\left[\alpha_t(x-q_t)^2 + p_t(x-q_t) + s_t\right]\right\}\]

and expanded the potential to second order about the moving centre \(q_t\). The Schrodinger equation is then solved exactly by

\[\dot q = \frac{p}{m},\quad \dot p = -V'(q),\quad \dot\alpha = -\frac{2\alpha^2}{m} - \frac{V''(q)}{2},\quad \dot s = \frac{i\hbar\alpha}{m} + \frac{p^2}{2m} - V(q).\]

The centre follows a classical trajectory, the complex width \(\alpha_t\) responds to the local curvature \(V''(q_t)\) (unlike the fixed width of frozen_gaussian_1d()), and \(s_t\) carries the classical action, the zero-point phase and the normalization. The result is exact for potentials at most quadratic, and fails once the packet grows wide enough to feel the third derivative of the potential.

The Riccati equation for \(\alpha\) becomes stiff where the packet focuses, so it is integrated in the equivalent linear form \(\alpha=P/2Z\) with \(\dot Z=P/m\), \(\dot P=-V''(q)Z\) (the tangent dynamics of the monodromy matrix), and then \(s_t=s_0+S_t+\tfrac{i\hbar}{2}\ln(Z_t/Z_0)\), with \(S_t\) the classical action and the phase of \(Z_t\) followed continuously.

Parameters:

qc0 (float) – Initial centre in position and momentum.
pc0 (float) – Initial centre in position and momentum.
gamma (float) – Initial width parameter, as in frozen_gaussian_1d() (\(\alpha_0=i\hbar\gamma\)).
dVdx (callable) – Numba-jitted potential derivatives and the potential itself, as in propagate_trajectory_monodromy_action().
d2Vdx2 (callable) – Numba-jitted potential derivatives and the potential itself, as in propagate_trajectory_monodromy_action().
V (callable) – Numba-jitted potential derivatives and the potential itself, as in propagate_trajectory_monodromy_action().
m (float) – Particle mass.
dt (float) – Time step.
steps (int) – Number of RK4 steps.
hbar (float) – Value of \(\hbar\) to use.
params (ndarray | None) – Parameter vector passed through to dVdx/d2Vdx2/V.

Returns:

t (ndarray, shape (steps + 1,)) – Times.
q, p (ndarray, shape (steps + 1,)) – Centre of the packet.
alpha, s (ndarray of complex, shape (steps + 1,)) – Complex width and phase-normalization parameters, for thawed_gaussian_wavefunction().

See also

thawed_gaussian_wavefunction: Evaluates the packet on a grid.
herman_kluk_propagate_wavepacket: Many frozen Gaussians instead of one thawed one.

References

E. J. Heller, “Time-dependent approach to semiclassical dynamics,” J. Chem. Phys. 62, 1544-1555 (1975).

Examples

The ground state of a harmonic oscillator (\(\gamma=m\omega/2\hbar\)) keeps its width, and its phase advances by the zero-point energy:

>>> import numpy as np
>>> from numba import njit
>>> dVdx = njit(lambda q, params: q, cache=False)
>>> d2Vdx2 = njit(lambda q, params: 1.0, cache=False)
>>> V = njit(lambda q, params: 0.5 * q ** 2, cache=False)
>>> t, q, p, alpha, s = thawed_gaussian_propagate(0.0, 0.0, 0.5, dVdx, d2Vdx2, V, m=1.0, dt=0.01, steps=300)
>>> complex(np.round(alpha[-1], 8) + 0.0)  # + 0.0 clears a signed zero
0.5j
>>> round(float(s[-1].real), 6)  # -E_0 t = -0.5 * 3
-1.5

physicskit.semiclassical.thawed_gaussian_wavefunction(x, q, p, alpha, s, hbar=1.0)[source]#

Evaluate a thawed Gaussian \(\exp\{i[\alpha(x-q)^2+p(x-q)+s]/\hbar\}\).

Parameters:

x (ndarray) – Positions at which to evaluate the packet.
q (float) – Centre in position and momentum.
p (float) – Centre in position and momentum.
alpha (complex) – Width and phase-normalization parameters, one time slice of the output of thawed_gaussian_propagate().
s (complex) – Width and phase-normalization parameters, one time slice of the output of thawed_gaussian_propagate().
hbar (float) – Value of \(\hbar\) to use.

Return type:

Returns:

ndarray of complex – \(\psi(x)\), same shape as x.

Examples

At \(t=0\) the packet is the frozen Gaussian it started as:

>>> import numpy as np
>>> x = np.linspace(-5, 5, 201)
>>> gamma, hbar = 0.7, 1.0
>>> s0 = -0.25j * hbar * np.log(2 * gamma / np.pi)
>>> psi = thawed_gaussian_wavefunction(x, 0.5, 1.2, 1j * hbar * gamma, s0, hbar)
>>> bool(np.allclose(psi, frozen_gaussian_1d(x, 0.5, 1.2, gamma, hbar)))
True

physicskit.semiclassical.turning_points(E, V, x_min, x_max, n_search=4000)[source]#

Locate the classical turning points (roots of \(E=V(x)\)) in \([x_{min},x_{max}]\).

Scans a fine grid for sign changes of \(E-V(x)\) and refines each with Brent’s method; returns every turning point found, not just the innermost pair, so multi-well potentials are handled correctly.

Parameters:

E (float) – Total energy.
V (callable) – Potential energy function V(x).
x_min (float) – Search domain (should extend into the classically forbidden region on both sides of the state of interest).
x_max (float) – Search domain (should extend into the classically forbidden region on both sides of the state of interest).
n_search (int) – Number of grid points used to bracket the roots; increase for potentials with closely spaced turning points.

Return type:

list

Returns:

list of float – Turning points, ascending.

See also

wkb_action: Integrates classical_momentum() between two turning points.

Examples

>>> V = lambda x: 0.5 * x ** 2
>>> [round(t, 4) for t in turning_points(E=2.0, V=V, x_min=-10, x_max=10)]
[-2.0, 2.0]

physicskit.semiclassical.van_vleck_prefactor(Mqp, hbar=1.0)[source]#

Van Vleck-Morette amplitude \(\sqrt{1/(2\pi\hbar\,|\partial q_t/\partial p_0|)}\).

The prefactor is the square root of (minus) the mixed second derivative of the classical action, \(\left|\partial^2S/\partial q_0\partial q_t\right| = 1/|\partial q_t/\partial p_0|\) – large where nearby trajectories stay close together (many classical paths reinforce the same endpoint) and formally divergent at a focal point (\(\partial q_t/\partial p_0=0\), where infinitesimally different initial momenta all reconverge on the same final position).

Parameters:

Mqp (float) – The \((0,1)\) element of the monodromy matrix, \(\partial q_t/\partial p_0\), from propagate_trajectory_monodromy_action().
hbar (float) – Value of \(\hbar\) to use.

Return type:

Returns:

float – The (real, non-negative) prefactor amplitude.

Examples

>>> round(van_vleck_prefactor(Mqp=1.0, hbar=1.0), 6)
0.398942

physicskit.semiclassical.van_vleck_propagator_1d(q0, p0, dVdx, d2Vdx2, V, m, dt, steps, hbar=1.0, params=None)[source]#

Semiclassical (Van Vleck-Morette) propagator amplitude along one classical trajectory.

\[K(q_t,t;q_0,0) \approx \sqrt{\frac{1}{2\pi\hbar\,|\partial q_t/\partial p_0|}}\, \exp\!\left[\frac{i}{\hbar}S(t) - i\frac{\pi}{4} - i\mu\frac{\pi}{2}\right],\]

the leading-order (\(\hbar\to0\)) approximation to the quantum propagator, exact whenever the potential is at most quadratic (the free particle and the harmonic oscillator), since then the WKB expansion this formula comes from truncates exactly.

Parameters:

q0 (float) – Initial position and momentum. p0 fixes which trajectory (and hence which final position q_t) is used; this function does not solve the two-point boundary-value problem of finding the p0 that reaches a prescribed q_t.
p0 (float) – Initial position and momentum. p0 fixes which trajectory (and hence which final position q_t) is used; this function does not solve the two-point boundary-value problem of finding the p0 that reaches a prescribed q_t.
dVdx (callable) – Numba-jitted potential derivatives and the potential itself, as in propagate_trajectory_monodromy_action().
d2Vdx2 (callable) – Numba-jitted potential derivatives and the potential itself, as in propagate_trajectory_monodromy_action().
V (callable) – Numba-jitted potential derivatives and the potential itself, as in propagate_trajectory_monodromy_action().
m (float) – Particle mass.
dt (float) – Time step.
steps (int) – Number of RK4 steps.
hbar (float) – Value of \(\hbar\) to use.
params (ndarray | None) – Parameter vector passed through to dVdx/d2Vdx2/V.

Returns:

q_t (float) – The trajectory’s final position.
K (complex) – The semiclassical propagator amplitude \(K(q_t,t;q_0,0)\).

See also

propagate_trajectory_monodromy_action: Supplies the trajectory, monodromy, and action.
herman_kluk_propagate_wavepacket: Sums many such trajectories into a full wavepacket propagator.

Examples

A free particle’s semiclassical propagator matches the exact quantum propagator \(\sqrt{m/(2\pi i\hbar t)}\exp[im(q_t-q_0)^2/(2\hbar t)]\) exactly:

>>> import numpy as np
>>> from numba import njit
>>> zero = njit(lambda q, params: 0.0, cache=False)
>>> q_t, K = van_vleck_propagator_1d(q0=0.0, p0=1.0, dVdx=zero, d2Vdx2=zero, V=zero, m=1.0, dt=1.0 / 2000, steps=2000)
>>> K_exact = np.sqrt(1.0 / (2j * np.pi * 1.0)) * np.exp(1j * (q_t - 0.0) ** 2 / 2.0)
>>> bool(abs(K - K_exact) < 1e-6)
True

A harmonic-oscillator trajectory that has not yet crossed a focal point (\(t<\pi/\omega\)) matches the exact Van Vleck (Mehler) propagator just as closely:

>>> m, omega = 1.0, 1.0
>>> params = np.array([m * omega ** 2])
>>> dVdx = njit(lambda q, params: params[0] * q, cache=False)
>>> d2Vdx2 = njit(lambda q, params: params[0], cache=False)
>>> V = njit(lambda q, params: 0.5 * params[0] * q ** 2, cache=False)
>>> q0, p0, t = 1.0, 0.3, 1.3
>>> q_t, K = van_vleck_propagator_1d(q0, p0, dVdx, d2Vdx2, V, m, dt=t / 4000, steps=4000, params=params)
>>> K_exact = (np.sqrt(m * omega / (2j * np.pi * np.sin(omega * t)))
...            * np.exp(1j * m * omega / (2 * np.sin(omega * t)) * ((q_t ** 2 + q0 ** 2) * np.cos(omega * t) - 2 * q_t * q0)))
>>> bool(abs(K - K_exact) < 1e-5)
True

physicskit.semiclassical.wkb_action(E, V, m, x1, x2, hbar=1.0)[source]#

WKB action integral \(S(E)=\int_{x_1}^{x_2} p(x)\,dx\) between two turning points.

This is the classical action accumulated on one traversal of the classically allowed region – half of the full round-trip action \(\oint p\,dx\), and (as a function of \(E\)) exactly the Hamilton principal function whose energy derivative gives the classical period, \(T(E)=dS/dE\), used by physicskit.semiclassical.core.gutzwiller.classical_period().

Parameters:

E (float) – Total energy.
V (callable) – Potential energy function V(x).
m (float) – Particle mass.
x1 (float) – The two turning points (e.g. from turning_points()).
x2 (float) – The two turning points (e.g. from turning_points()).
hbar (float) – Value of \(\hbar\) to use (only needed by callers converting this to a phase; the integral itself is \(\hbar\)-independent).

Return type:

Returns:

float – The action \(S(E)\).

Examples

For the harmonic oscillator, \(S(E)=\pi E/\omega\) exactly:

>>> V = lambda x: 0.5 * x ** 2
>>> x1, x2 = turning_points(E=2.0, V=V, x_min=-10, x_max=10)
>>> round(wkb_action(E=2.0, V=V, m=1.0, x1=x1, x2=x2), 6)
6.283185

physicskit.semiclassical.wkb_wavefunction(x, E, V, m=1.0, hbar=1.0)[source]#

The (real-valued, standing-wave) WKB wavefunction in the classically allowed region.

\[\psi(x) \approx \frac{C}{\sqrt{p(x)}} \cos\!\left(\frac{1}{\hbar}\int_{x_1}^{x}p(x')\,dx' - \frac{\pi}{4}\right),\]

with the \(\pi/4\) phase set by the Airy connection formula at the left (soft) turning point \(x_1\), and \(C\) fixed by numerical normalization over the allowed region. Diverges (is not valid) within a few de Broglie wavelengths of either turning point; zero is returned in the classically forbidden region rather than the (also invalid) naive imaginary-momentum continuation.

Parameters:

x (ndarray) – Positions at which to evaluate the wavefunction.
E (float) – Energy (e.g. from bohr_sommerfeld_energies(), for a properly quantized bound state).
V (callable) – Potential energy function V(x).
m (float) – Particle mass.
hbar (float) – Value of \(\hbar\) to use.

Return type:

Returns:

ndarray – \(\psi(x)\), same shape as x, normalized so \(\int|\psi|^2\,dx=1\) over x’s range.

See also

bohr_sommerfeld_energies: Quantized energies this is evaluated at.

Examples

The WKB wavefunction at the Bohr-Sommerfeld energy for quantum number n has exactly n nodes, the standard node-counting theorem:

>>> import numpy as np
>>> V = lambda x: 0.5 * x ** 2
>>> E5 = bohr_sommerfeld_energies(V, m=1.0, x_min=-20, x_max=20, n_max=6)[5]
>>> x = np.linspace(-10, 10, 4000)
>>> psi = wkb_wavefunction(x, E5, V, m=1.0)
>>> nonzero = psi[psi != 0]
>>> int(np.sum(np.diff(np.sign(nonzero)) != 0))
5

1D WKB wavefunctions, classical turning points, and Bohr-Sommerfeld (EBK) quantization.

The Wentzel-Kramers-Brillouin (WKB) approximation writes the wavefunction as \(\psi(x) \sim p(x)^{-1/2}\exp(\pm i\int p\,dx/\hbar)\) – an \(\hbar\to 0\) asymptotic expansion valid wherever the local de Broglie wavelength \(2\pi\hbar/p(x)\) varies slowly compared to itself. The prefactor \(p(x)^{-1/2}\) is exactly the classical probability density of a particle oscillating in the potential (it spends more time, and so is more likely to be found, where it moves slowest); this is the one-dimensional seed of every semiclassical idea in physicskit.semiclassical – the Van Vleck-Morette determinant in propagators is the same classical-probability prefactor for a propagator instead of a stationary state, and the periodic-orbit sum in gutzwiller is built from the same action integral used here for Bohr-Sommerfeld quantization.

At each classical turning point the WKB approximation itself breaks down (\(p(x)\to 0\)); matching the oscillatory interior solution through that breakdown region (via the Airy function connection formulas) costs each soft (linear) turning point a phase of \(\pi/4\). For a bound state with two such turning points, the round-trip quantization condition \(\oint p\,dx = 2\pi\hbar(n+\tfrac12)\) – equivalently \(\int_{x_1}^{x_2}p\,dx=(n+\tfrac12)\pi\hbar\) – is the one-dimensional Einstein-Brillouin-Keller (EBK) rule implemented by bohr_sommerfeld_energies().

For the radial motion of a central-force problem the same rule fails unless the centrifugal term is changed: langer_corrected_wkb() applies Langer’s (1937) replacement \(l(l+1)\to(l+\tfrac12)^2\), which makes radial WKB give the exact hydrogen spectrum.

physicskit.semiclassical.core.wkb.bohr_sommerfeld_energies(V, m, x_min, x_max, n_max, hbar=1.0, E_min=0.001, E_max=100.0)[source]#

Bohr-Sommerfeld (EBK) bound-state energies from the WKB quantization condition.

Solves \(\int_{x_1(E_n)}^{x_2(E_n)}p(x)\,dx=(n+\tfrac12)\pi\hbar\) for each \(n=0,\dots,n_{max}-1\) by root-finding on \(E\), re-locating the turning points at every trial energy. For the harmonic oscillator this reproduces the exact spectrum \(E_n=(n+\tfrac12)\hbar\omega\) to machine precision, since the WKB approximation is exact whenever the potential is exactly quadratic.

Parameters:

V (callable) – Potential energy function V(x), with a single classically allowed region containing its minimum for every energy searched.
m (float) – Particle mass.
x_min (float) – Search domain for turning_points(); must extend well into the classically forbidden region at E_max.
x_max (float) – Search domain for turning_points(); must extend well into the classically forbidden region at E_max.
n_max (int) – Number of levels to compute (\(n=0,\dots,n_{max}-1\)).
hbar (float) – Value of \(\hbar\) to use.
E_min (float) – Lower bracket for the ground-state root search; must be small enough that its turning points still resolve on the turning_points() search grid.
E_max (float) – Upper bracket for the highest level’s root search.

Return type:

Returns:

ndarray, shape (n_max,) – Quantized energies, ascending.

See also

wkb_action: The action integral this quantizes.
wkb_wavefunction: The corresponding approximate wavefunctions.

Examples

>>> import numpy as np
>>> V = lambda x: 0.5 * x ** 2
>>> energies = bohr_sommerfeld_energies(V, m=1.0, x_min=-20, x_max=20, n_max=5)
>>> np.round(energies, 6)
array([0.5, 1.5, 2.5, 3.5, 4.5])

physicskit.semiclassical.core.wkb.classical_momentum(E, V, x, m=1.0)[source]#

Classical momentum \(p(x) = \sqrt{2m(E-V(x))}\) in the allowed region.

Parameters:

E (float) – Total energy.
V (callable) – Potential energy function V(x).
x (ndarray) – Positions at which to evaluate \(p(x)\).
m (float) – Particle mass.

Return type:

Returns:

ndarray – \(p(x)\), same shape as x; nan wherever \(V(x)>E\) (the classically forbidden region, where WKB oscillates into a real exponential instead).

Examples

>>> import numpy as np
>>> V = lambda x: 0.5 * x ** 2
>>> p = classical_momentum(E=2.0, V=V, x=np.array([0.0, 1.0, 3.0]))
>>> np.round(p, 4)
array([2.    , 1.7321,    nan])

physicskit.semiclassical.core.wkb.langer_corrected_wkb(V, l, m, r_max, n_max, hbar=1.0, langer=True, E_max=None, r_min=1e-08, n_search=4000)[source]#

Radial WKB bound-state energies, with or without the Langer correction.

The radial Schrodinger equation for \(u(r)=rR(r)\) is a 1D problem on \(r>0\) with effective potential \(V(r)+\hbar^2 l(l+1)/(2mr^2)\). Applying the WKB rule \(\int_{r_1}^{r_2}p\,dr=(n_r+\tfrac12)\pi\hbar\) to it directly gives wrong energies, because the WKB phase near \(r=0\) is not the phase of the true solution \(u\sim r^{l+1}\). Langer (1937) showed that the substitution \(r=e^x\) turns the radial problem into one on the whole line, where WKB applies cleanly, and that the result is the same as using

\[V_{\rm eff}(r) = V(r) + \frac{\hbar^2(l+\tfrac12)^2}{2mr^2}\]

in the ordinary WKB rule. With this correction WKB gives the exact hydrogen levels \(E=-mk^2/(2\hbar^2(n_r+l+1)^2)\) for \(V=-k/r\), and the exact 3D harmonic-oscillator levels \(\hbar\omega(2n_r+l+\tfrac32)\).

Parameters:

V (callable) – Radial potential V(r), vectorized over a NumPy array.
l (int) – Orbital angular momentum quantum number.
m (float) – Particle mass.
r_max (float) – Outer edge of the radial domain; every outer turning point must lie inside it.
n_max (int) – Number of radial levels to compute (\(n_r=0,\dots,n_{max}-1\)).
hbar (float) – Value of \(\hbar\) to use.
langer (bool) – Use \((l+\tfrac12)^2\); if False, use the uncorrected \(l(l+1)\).
E_max (float | None) – Upper bracket for the root search. Defaults to \(V_{\rm eff}(r_{max})\), the highest energy whose outer turning point is still inside the domain.
r_min (float) – Inner edge of the radial domain. If \(V_{\rm eff}\) has no inner turning point (\(l=0\) without the Langer correction, for a Coulomb potential), the action integral starts here.
n_search (int) – Number of log-spaced grid points used to bracket turning points.

Return type:

Returns:

ndarray, shape (n_max,) – Energies for \(n_r=0,\dots,n_{max}-1\), ascending.

Raises:

ValueError – If a level lies above E_max (increase r_max).

See also

bohr_sommerfeld_energies: The same rule for a potential on the whole line.

References

R. E. Langer, “On the Connection Formulas and the Solutions of the Wave Equation,” Phys. Rev. 51, 669-676 (1937).

Examples

With the correction, WKB reproduces the hydrogen 2p, 3p, 4p levels \(-1/(2n^2)\) exactly; without it, it does not:

>>> import numpy as np
>>> V = lambda r: -1.0 / r
>>> np.round(langer_corrected_wkb(V, l=1, m=1.0, r_max=200.0, n_max=3), 6)
array([-0.125   , -0.055556, -0.03125 ])
>>> np.round(langer_corrected_wkb(V, l=1, m=1.0, r_max=200.0, n_max=3, langer=False), 4)
array([-0.1365, -0.0589, -0.0326])

physicskit.semiclassical.core.wkb.turning_points(E, V, x_min, x_max, n_search=4000)[source]#

Locate the classical turning points (roots of \(E=V(x)\)) in \([x_{min},x_{max}]\).

Scans a fine grid for sign changes of \(E-V(x)\) and refines each with Brent’s method; returns every turning point found, not just the innermost pair, so multi-well potentials are handled correctly.

Parameters:

E (float) – Total energy.
V (callable) – Potential energy function V(x).
x_min (float) – Search domain (should extend into the classically forbidden region on both sides of the state of interest).
x_max (float) – Search domain (should extend into the classically forbidden region on both sides of the state of interest).
n_search (int) – Number of grid points used to bracket the roots; increase for potentials with closely spaced turning points.

Return type:

list

Returns:

list of float – Turning points, ascending.

See also

wkb_action: Integrates classical_momentum() between two turning points.

Examples

>>> V = lambda x: 0.5 * x ** 2
>>> [round(t, 4) for t in turning_points(E=2.0, V=V, x_min=-10, x_max=10)]
[-2.0, 2.0]

physicskit.semiclassical.core.wkb.wkb_action(E, V, m, x1, x2, hbar=1.0)[source]#

WKB action integral \(S(E)=\int_{x_1}^{x_2} p(x)\,dx\) between two turning points.

This is the classical action accumulated on one traversal of the classically allowed region – half of the full round-trip action \(\oint p\,dx\), and (as a function of \(E\)) exactly the Hamilton principal function whose energy derivative gives the classical period, \(T(E)=dS/dE\), used by physicskit.semiclassical.core.gutzwiller.classical_period().

Parameters:

E (float) – Total energy.
V (callable) – Potential energy function V(x).
m (float) – Particle mass.
x1 (float) – The two turning points (e.g. from turning_points()).
x2 (float) – The two turning points (e.g. from turning_points()).
hbar (float) – Value of \(\hbar\) to use (only needed by callers converting this to a phase; the integral itself is \(\hbar\)-independent).

Return type:

Returns:

float – The action \(S(E)\).

Examples

For the harmonic oscillator, \(S(E)=\pi E/\omega\) exactly:

>>> V = lambda x: 0.5 * x ** 2
>>> x1, x2 = turning_points(E=2.0, V=V, x_min=-10, x_max=10)
>>> round(wkb_action(E=2.0, V=V, m=1.0, x1=x1, x2=x2), 6)
6.283185

physicskit.semiclassical.core.wkb.wkb_wavefunction(x, E, V, m=1.0, hbar=1.0)[source]#

The (real-valued, standing-wave) WKB wavefunction in the classically allowed region.

\[\psi(x) \approx \frac{C}{\sqrt{p(x)}} \cos\!\left(\frac{1}{\hbar}\int_{x_1}^{x}p(x')\,dx' - \frac{\pi}{4}\right),\]

with the \(\pi/4\) phase set by the Airy connection formula at the left (soft) turning point \(x_1\), and \(C\) fixed by numerical normalization over the allowed region. Diverges (is not valid) within a few de Broglie wavelengths of either turning point; zero is returned in the classically forbidden region rather than the (also invalid) naive imaginary-momentum continuation.

Parameters:

x (ndarray) – Positions at which to evaluate the wavefunction.
E (float) – Energy (e.g. from bohr_sommerfeld_energies(), for a properly quantized bound state).
V (callable) – Potential energy function V(x).
m (float) – Particle mass.
hbar (float) – Value of \(\hbar\) to use.

Return type:

Returns:

ndarray – \(\psi(x)\), same shape as x, normalized so \(\int|\psi|^2\,dx=1\) over x’s range.

See also

bohr_sommerfeld_energies: Quantized energies this is evaluated at.

Examples

The WKB wavefunction at the Bohr-Sommerfeld energy for quantum number n has exactly n nodes, the standard node-counting theorem:

>>> import numpy as np
>>> V = lambda x: 0.5 * x ** 2
>>> E5 = bohr_sommerfeld_energies(V, m=1.0, x_min=-20, x_max=20, n_max=6)[5]
>>> x = np.linspace(-10, 10, 4000)
>>> psi = wkb_wavefunction(x, E5, V, m=1.0)
>>> nonzero = psi[psi != 0]
>>> int(np.sum(np.diff(np.sign(nonzero)) != 0))
5

Semiclassical propagators: the Van Vleck-Morette determinant and Herman-Kluk frozen Gaussians.

Both propagators here reconstruct the quantum time-evolution operator from classical trajectories alone. Van Vleck-Morette theory (van_vleck_propagator_1d()) uses a single trajectory connecting a fixed start and end point, dressed with a stability prefactor built from its monodromy matrix – the linearized map \((\delta q_0,\delta p_0)\mapsto(\delta q_t,\delta p_t)\) obtained by propagating small deviations alongside the trajectory (propagate_trajectory_monodromy_action()). Herman and Kluk’s 1984 frozen-Gaussian method (herman_kluk_propagate_wavepacket()) instead sums the contributions of many such trajectories, one launched from each point of a phase-space grid under the initial wavepacket, each one carrying its own rigid (frozen-width) Gaussian, monodromy-built prefactor, and classical action phase. Heller’s 1975 thawed Gaussian (thawed_gaussian_propagate()) sits between the two: a single trajectory carries a single Gaussian whose complex width evolves with the local curvature of the potential.

Right-hand-side functions passed to the trajectory/monodromy integrator must be module-level @njit functions with signature dVdx(q, params) -> float / d2Vdx2(q, params) -> float / V(q, params) -> float, with params a float64 array – the same convention used throughout this package’s Numba-accelerated integrators, letting a single compiled kernel serve any potential. For the same reason as physicskit.integrators.fixed_step, the kernels that take these functions as arguments are not cache=True.

physicskit.semiclassical.core.propagators.coherent_state_overlap(q1, p1, q2, p2, gamma, hbar=1.0)[source]#

Closed-form overlap \(\langle g_{q_1,p_1}|g_{q_2,p_2}\rangle\) of two equal-width frozen Gaussians.

\[\langle g_1|g_2\rangle = \exp\!\left[-\frac{\gamma}{2}(q_1-q_2)^2 - \frac{(p_1-p_2)^2}{8\gamma\hbar^2} + \frac{i}{2\hbar}(q_1-q_2)(p_1+p_2)\right].\]

Parameters:

q1 (float) – Phase-space center of the bra state.
p1 (float) – Phase-space center of the bra state.
q2 (float) – Phase-space center of the ket state.
p2 (float) – Phase-space center of the ket state.
gamma (float) – Shared width parameter, as in frozen_gaussian_1d().
hbar (float) – Value of \(\hbar\) to use.

Return type:

Returns:

complex – The overlap \(\langle g_1|g_2\rangle\).

See also

frozen_gaussian_1d: The states being overlapped.

Examples

Matches direct numerical integration of the two wavepackets:

>>> import numpy as np
>>> x = np.linspace(-40, 40, 20000)
>>> gamma = 1.0
>>> psi1 = frozen_gaussian_1d(x, qc=0.3, pc=0.7, gamma=gamma)
>>> psi2 = frozen_gaussian_1d(x, qc=-0.2, pc=1.1, gamma=gamma)
>>> numeric = np.trapezoid(np.conj(psi1) * psi2, x)
>>> closed_form = coherent_state_overlap(0.3, 0.7, -0.2, 1.1, gamma)
>>> bool(abs(numeric - closed_form) < 1e-6)
True

A state’s overlap with itself is 1:

>>> round(float(abs(coherent_state_overlap(1.0, 2.0, 1.0, 2.0, gamma=0.8))), 8)
1.0

physicskit.semiclassical.core.propagators.count_caustics(Mqp_history)[source]#

Count sign changes of \(\partial q_t/\partial p_0\) along a trajectory – the Maslov index.

Each time \(\partial q_t/\partial p_0\) passes through zero, the trajectory crosses a focal point (conjugate point / caustic), and the semiclassical propagator picks up an extra phase of \(-\pi/2\) (Gutzwiller 1967; see also the \(-\pi/4\) phase at a single WKB turning point in physicskit.semiclassical.core.wkb.wkb_wavefunction(), which is this same phenomenon in the time-independent picture).

Parameters:: Mqp_history (ndarray) – \(\partial q_t/\partial p_0\) sampled along the trajectory, from propagate_trajectory_monodromy_action(); the first entry (always 0, at \(t=0\)) is ignored.
Return type:: int
Returns:: int – The Maslov index \(\mu\).

See also

van_vleck_propagator_1d: Uses this to fix the propagator’s overall phase.

Examples

>>> count_caustics(np.array([0.0, 1.0, 0.5, -0.3, -0.8, 0.2]))
2

physicskit.semiclassical.core.propagators.frozen_gaussian_1d(x, qc, pc, gamma, hbar=1.0)[source]#

A normalized, fixed-width (“frozen”) Gaussian wavepacket (coherent state).

\[g_{q_c,p_c}(x) = \left(\frac{2\gamma}{\pi}\right)^{1/4} \exp\!\left[-\gamma(x-q_c)^2 + \frac{i}{\hbar}p_c(x-q_c)\right].\]

The building block of the Herman-Kluk method: every trajectory carries one of these, always at the fixed width set by \(\gamma\) (hence “frozen”), riding on top of its classical phase-space point \((q_c,p_c)\).

Parameters:

x (ndarray) – Positions at which to evaluate the wavepacket.
qc (float) – Center of the Gaussian in position and momentum.
pc (float) – Center of the Gaussian in position and momentum.
gamma (float) – Width parameter (inverse squared length); larger \(\gamma\) means a narrower, more position-localized packet.
hbar (float) – Value of \(\hbar\) to use.

Return type:

Returns:

ndarray of complex – \(g_{q_c,p_c}(x)\), same shape as x.

See also

coherent_state_overlap: The closed-form overlap of two such Gaussians.

Examples

>>> import numpy as np
>>> x = np.linspace(-20, 20, 4000)
>>> psi = frozen_gaussian_1d(x, qc=1.0, pc=2.0, gamma=0.5)
>>> round(float(np.trapezoid(np.abs(psi) ** 2, x)), 8)
1.0

physicskit.semiclassical.core.propagators.herman_kluk_prefactor(M, gamma, hbar=1.0)[source]#

Herman-Kluk prefactor \(C_t(q_0,p_0)\) built from the monodromy matrix.

\[C_t = \sqrt{\frac{1}{2}\left(M_{qq} + M_{pp} - 2i\hbar\gamma M_{qp} + \frac{i}{2\hbar\gamma}M_{pq}\right)},\]

the standard Herman-Kluk form (Herman & Kluk 1984; Kay 1994) written for this module’s frozen Gaussians \(e^{-\gamma(x-q)^2}\), whose width parameter in the usual \(e^{-\gamma_s(x-q)^2/2}\) convention is \(\gamma_s=2\gamma\). It reduces to \(C_t=1\) at \(t=0\) (\(M=\mathbb{1}\)), consistent with the frozen Gaussian basis’s resolution of the identity, \(\int \tfrac{dq\,dp}{2\pi\hbar}\,|g_{q,p}\rangle\langle g_{q,p}| = \hat{1}\).

Parameters:

M (ndarray) – Monodromy matrix from propagate_trajectory_monodromy_action().
gamma (float) – Frozen-Gaussian width parameter, as in frozen_gaussian_1d().
hbar (float) – Value of \(\hbar\) to use.

Return type:

Returns:

complex – The prefactor \(C_t\).

Examples

>>> import numpy as np
>>> complex(herman_kluk_prefactor(np.eye(2), gamma=1.0))
(1+0j)

physicskit.semiclassical.core.propagators.herman_kluk_propagate_wavepacket(qc0, pc0, gamma, dVdx, d2Vdx2, V, m, dt, steps, x_eval, hbar=1.0, n_grid=41, n_sigma=6.0, params=None)[source]#

Propagate an initial frozen-Gaussian wavepacket with the multi-trajectory Herman-Kluk method.

Launches one classical trajectory from every point of a regular grid covering the initial coherent state’s phase-space support, and sums their frozen-Gaussian contributions

\[\psi(x,t) \approx \int\!\frac{dq_0\,dp_0}{2\pi\hbar}\, C_t(q_0,p_0)\,e^{iS(t)/\hbar}\, \langle g_{q_0,p_0}|\psi_0\rangle\, g_{q_t,p_t}(x),\]

approximating the integral over initial conditions by a Riemann sum on a grid spanning \(\pm n_\sigma\) standard deviations of the initial coherent state’s own phase-space Gaussian (\(\sigma_q=1/(2\sqrt\gamma)\), \(\sigma_p=\hbar\sqrt\gamma\)). The per-trajectory classical propagation (propagate_trajectory_monodromy_action()) is Numba-compiled, since it is called once per grid point – the dominant cost for any reasonably fine grid.

In the \(t\to0\) limit this reduces to the frozen-Gaussian resolution of the identity and reconstructs the initial wavepacket essentially exactly; away from that limit, this quadrature-grid evaluation of the Herman-Kluk integral does not exactly conserve \(\int|\psi|^2\,dx\) (the true continuous phase-space integral does, for at-most-quadratic potentials) – convergence in n_grid, n_sigma, and gamma should be checked for any serious use.

Parameters:

qc0 (float) – Center of the initial frozen-Gaussian wavepacket.
pc0 (float) – Center of the initial frozen-Gaussian wavepacket.
gamma (float) – Width parameter shared by the initial state and every frozen Gaussian in the propagation.
dVdx (callable) – Numba-jitted potential derivatives and the potential itself, as in propagate_trajectory_monodromy_action().
d2Vdx2 (callable) – Numba-jitted potential derivatives and the potential itself, as in propagate_trajectory_monodromy_action().
V (callable) – Numba-jitted potential derivatives and the potential itself, as in propagate_trajectory_monodromy_action().
m (float) – Particle mass.
dt (float) – Time step for each trajectory’s RK4 integration.
steps (int) – Number of RK4 steps (propagation time is dt * steps).
x_eval (ndarray) – Positions at which to evaluate the propagated wavepacket.
hbar (float) – Value of \(\hbar\) to use.
n_grid (int) – Number of grid points along each of the \(q_0\), \(p_0\) axes.
n_sigma (float) – Half-width of the sampling grid, in standard deviations of the initial state’s phase-space Gaussian.
params (ndarray | None) – Parameter vector passed through to dVdx/d2Vdx2/V.

Return type:

Returns:

ndarray of complex, shape matching x_eval – The propagated wavepacket \(\psi(x,t)\).

See also

van_vleck_propagator_1d: The single-trajectory propagator this sums many copies of.
frozen_gaussian_1d: The initial state and the basis each trajectory carries.

Examples

At very short times the sum reduces to the frozen-Gaussian resolution of the identity and reproduces the (unpropagated) initial coherent state to high accuracy:

>>> import numpy as np
>>> from numba import njit
>>> m, omega = 1.0, 1.0
>>> params = np.array([m * omega ** 2])
>>> dVdx = njit(lambda q, params: params[0] * q, cache=False)
>>> d2Vdx2 = njit(lambda q, params: params[0], cache=False)
>>> V = njit(lambda q, params: 0.5 * params[0] * q ** 2, cache=False)
>>> x = np.linspace(-6, 6, 400)
>>> psi = herman_kluk_propagate_wavepacket(
...     1.0, 0.0, gamma=1.0, dVdx=dVdx, d2Vdx2=d2Vdx2, V=V, m=m, dt=1e-6, steps=1, x_eval=x, n_grid=61, n_sigma=7.0, params=params
... )
>>> psi0 = frozen_gaussian_1d(x, qc=1.0, pc=0.0, gamma=1.0)
>>> fidelity = abs(np.trapezoid(np.conj(psi0) * psi, x)) ** 2
>>> bool(fidelity > 0.999)
True

physicskit.semiclassical.core.propagators.propagate_trajectory_monodromy_action(q0, p0, dVdx, d2Vdx2, V, m, dt, steps, params=None)[source]#

Integrate a classical trajectory together with its monodromy matrix and Hamilton principal function.

Advances the Hamilton equations \(\dot q=p/m,\ \dot p=-V'(q)\) together with the variational (tangent) equations for small deviations, \(\dot{\delta q}=\delta p/m,\ \dot{\delta p} =-V''(q)\delta q\), whose solution operator is the monodromy matrix \(M(t)=\begin{pmatrix}\partial q_t/\partial q_0 & \partial q_t/\partial p_0\\ \partial p_t/\partial q_0 & \partial p_t/\partial p_0\end{pmatrix}\), and accumulates the classical action (Hamilton’s principal function) \(S(t)=\int_0^t\left[\tfrac{p^2}{2m}-V(q)\right]dt'\) along the way – everything van_vleck_propagator_1d() and herman_kluk_propagate_wavepacket() need from one trajectory, computed in a single RK4 pass.

Because the flow is Hamiltonian, \(\det M(t)=1\) exactly for all \(t\) (Liouville’s theorem); this is a good numerical sanity check on any trajectory this function returns.

Parameters:

q0 (float) – Initial position and momentum.
p0 (float) – Initial position and momentum.
dVdx (callable) – Numba-jitted force law dVdx(q, params) -> float, \(V'(q)\).
d2Vdx2 (callable) – Numba-jitted curvature d2Vdx2(q, params) -> float, \(V''(q)\).
V (callable) – Numba-jitted potential V(q, params) -> float.
m (float) – Particle mass.
dt (float) – Time step.
steps (int) – Number of RK4 steps (total propagation time is dt * steps).
params (ndarray | None) – Parameter vector passed through to dVdx/d2Vdx2/V. Defaults to an empty array.

Returns:

q_t, p_t (float) – Final position and momentum.
M (ndarray, shape (2, 2)) – Monodromy matrix at time t = dt * steps.
S (float) – Classical action accumulated along the trajectory.
Mqp_history (ndarray, shape (steps + 1,)) – \(\partial q_t/\partial p_0\) at every step, for count_caustics().

See also

van_vleck_prefactor: Turns M into a propagator amplitude.
count_caustics: Turns Mqp_history into a Maslov index.

Examples

For the harmonic oscillator, the monodromy matrix has the exact closed form \(\begin{pmatrix}\cos\omega t & \sin(\omega t)/(m\omega)\\ -m\omega\sin\omega t & \cos\omega t\end{pmatrix}\):

>>> import numpy as np
>>> from numba import njit
>>> m, omega = 1.0, 1.0
>>> params = np.array([m * omega ** 2])
>>> dVdx = njit(lambda q, params: params[0] * q, cache=False)
>>> d2Vdx2 = njit(lambda q, params: params[0], cache=False)
>>> V = njit(lambda q, params: 0.5 * params[0] * q ** 2, cache=False)
>>> t = 1.3
>>> q_t, p_t, M, S, _ = propagate_trajectory_monodromy_action(1.0, 0.3, dVdx, d2Vdx2, V, m, dt=t / 4000, steps=4000, params=params)
>>> M_exact = np.array([[np.cos(t), np.sin(t)], [-np.sin(t), np.cos(t)]])
>>> bool(np.max(np.abs(M - M_exact)) < 1e-6)
True
>>> round(float(np.linalg.det(M)), 8)
1.0

physicskit.semiclassical.core.propagators.thawed_gaussian_propagate(qc0, pc0, gamma, dVdx, d2Vdx2, V, m, dt, steps, hbar=1.0, params=None)[source]#

Propagate a Gaussian wavepacket with Heller’s thawed Gaussian approximation.

Heller (1975) wrote the packet as

\[\psi(x,t) = \exp\!\left\{\frac{i}{\hbar}\left[\alpha_t(x-q_t)^2 + p_t(x-q_t) + s_t\right]\right\}\]

and expanded the potential to second order about the moving centre \(q_t\). The Schrodinger equation is then solved exactly by

\[\dot q = \frac{p}{m},\quad \dot p = -V'(q),\quad \dot\alpha = -\frac{2\alpha^2}{m} - \frac{V''(q)}{2},\quad \dot s = \frac{i\hbar\alpha}{m} + \frac{p^2}{2m} - V(q).\]

The centre follows a classical trajectory, the complex width \(\alpha_t\) responds to the local curvature \(V''(q_t)\) (unlike the fixed width of frozen_gaussian_1d()), and \(s_t\) carries the classical action, the zero-point phase and the normalization. The result is exact for potentials at most quadratic, and fails once the packet grows wide enough to feel the third derivative of the potential.

The Riccati equation for \(\alpha\) becomes stiff where the packet focuses, so it is integrated in the equivalent linear form \(\alpha=P/2Z\) with \(\dot Z=P/m\), \(\dot P=-V''(q)Z\) (the tangent dynamics of the monodromy matrix), and then \(s_t=s_0+S_t+\tfrac{i\hbar}{2}\ln(Z_t/Z_0)\), with \(S_t\) the classical action and the phase of \(Z_t\) followed continuously.

Parameters:

qc0 (float) – Initial centre in position and momentum.
pc0 (float) – Initial centre in position and momentum.
gamma (float) – Initial width parameter, as in frozen_gaussian_1d() (\(\alpha_0=i\hbar\gamma\)).
dVdx (callable) – Numba-jitted potential derivatives and the potential itself, as in propagate_trajectory_monodromy_action().
d2Vdx2 (callable) – Numba-jitted potential derivatives and the potential itself, as in propagate_trajectory_monodromy_action().
V (callable) – Numba-jitted potential derivatives and the potential itself, as in propagate_trajectory_monodromy_action().
m (float) – Particle mass.
dt (float) – Time step.
steps (int) – Number of RK4 steps.
hbar (float) – Value of \(\hbar\) to use.
params (ndarray | None) – Parameter vector passed through to dVdx/d2Vdx2/V.

Returns:

t (ndarray, shape (steps + 1,)) – Times.
q, p (ndarray, shape (steps + 1,)) – Centre of the packet.
alpha, s (ndarray of complex, shape (steps + 1,)) – Complex width and phase-normalization parameters, for thawed_gaussian_wavefunction().

See also

thawed_gaussian_wavefunction: Evaluates the packet on a grid.
herman_kluk_propagate_wavepacket: Many frozen Gaussians instead of one thawed one.

References

E. J. Heller, “Time-dependent approach to semiclassical dynamics,” J. Chem. Phys. 62, 1544-1555 (1975).

Examples

The ground state of a harmonic oscillator (\(\gamma=m\omega/2\hbar\)) keeps its width, and its phase advances by the zero-point energy:

>>> import numpy as np
>>> from numba import njit
>>> dVdx = njit(lambda q, params: q, cache=False)
>>> d2Vdx2 = njit(lambda q, params: 1.0, cache=False)
>>> V = njit(lambda q, params: 0.5 * q ** 2, cache=False)
>>> t, q, p, alpha, s = thawed_gaussian_propagate(0.0, 0.0, 0.5, dVdx, d2Vdx2, V, m=1.0, dt=0.01, steps=300)
>>> complex(np.round(alpha[-1], 8) + 0.0)  # + 0.0 clears a signed zero
0.5j
>>> round(float(s[-1].real), 6)  # -E_0 t = -0.5 * 3
-1.5

physicskit.semiclassical.core.propagators.thawed_gaussian_wavefunction(x, q, p, alpha, s, hbar=1.0)[source]#

Evaluate a thawed Gaussian \(\exp\{i[\alpha(x-q)^2+p(x-q)+s]/\hbar\}\).

Parameters:

x (ndarray) – Positions at which to evaluate the packet.
q (float) – Centre in position and momentum.
p (float) – Centre in position and momentum.
alpha (complex) – Width and phase-normalization parameters, one time slice of the output of thawed_gaussian_propagate().
s (complex) – Width and phase-normalization parameters, one time slice of the output of thawed_gaussian_propagate().
hbar (float) – Value of \(\hbar\) to use.

Return type:

Returns:

ndarray of complex – \(\psi(x)\), same shape as x.

Examples

At \(t=0\) the packet is the frozen Gaussian it started as:

>>> import numpy as np
>>> x = np.linspace(-5, 5, 201)
>>> gamma, hbar = 0.7, 1.0
>>> s0 = -0.25j * hbar * np.log(2 * gamma / np.pi)
>>> psi = thawed_gaussian_wavefunction(x, 0.5, 1.2, 1j * hbar * gamma, s0, hbar)
>>> bool(np.allclose(psi, frozen_gaussian_1d(x, 0.5, 1.2, gamma, hbar)))
True

physicskit.semiclassical.core.propagators.van_vleck_prefactor(Mqp, hbar=1.0)[source]#

Van Vleck-Morette amplitude \(\sqrt{1/(2\pi\hbar\,|\partial q_t/\partial p_0|)}\).

The prefactor is the square root of (minus) the mixed second derivative of the classical action, \(\left|\partial^2S/\partial q_0\partial q_t\right| = 1/|\partial q_t/\partial p_0|\) – large where nearby trajectories stay close together (many classical paths reinforce the same endpoint) and formally divergent at a focal point (\(\partial q_t/\partial p_0=0\), where infinitesimally different initial momenta all reconverge on the same final position).

Parameters:

Mqp (float) – The \((0,1)\) element of the monodromy matrix, \(\partial q_t/\partial p_0\), from propagate_trajectory_monodromy_action().
hbar (float) – Value of \(\hbar\) to use.

Return type:

Returns:

float – The (real, non-negative) prefactor amplitude.

Examples

>>> round(van_vleck_prefactor(Mqp=1.0, hbar=1.0), 6)
0.398942

physicskit.semiclassical.core.propagators.van_vleck_propagator_1d(q0, p0, dVdx, d2Vdx2, V, m, dt, steps, hbar=1.0, params=None)[source]#

Semiclassical (Van Vleck-Morette) propagator amplitude along one classical trajectory.

\[K(q_t,t;q_0,0) \approx \sqrt{\frac{1}{2\pi\hbar\,|\partial q_t/\partial p_0|}}\, \exp\!\left[\frac{i}{\hbar}S(t) - i\frac{\pi}{4} - i\mu\frac{\pi}{2}\right],\]

the leading-order (\(\hbar\to0\)) approximation to the quantum propagator, exact whenever the potential is at most quadratic (the free particle and the harmonic oscillator), since then the WKB expansion this formula comes from truncates exactly.

Parameters:

q0 (float) – Initial position and momentum. p0 fixes which trajectory (and hence which final position q_t) is used; this function does not solve the two-point boundary-value problem of finding the p0 that reaches a prescribed q_t.
p0 (float) – Initial position and momentum. p0 fixes which trajectory (and hence which final position q_t) is used; this function does not solve the two-point boundary-value problem of finding the p0 that reaches a prescribed q_t.
dVdx (callable) – Numba-jitted potential derivatives and the potential itself, as in propagate_trajectory_monodromy_action().
d2Vdx2 (callable) – Numba-jitted potential derivatives and the potential itself, as in propagate_trajectory_monodromy_action().
V (callable) – Numba-jitted potential derivatives and the potential itself, as in propagate_trajectory_monodromy_action().
m (float) – Particle mass.
dt (float) – Time step.
steps (int) – Number of RK4 steps.
hbar (float) – Value of \(\hbar\) to use.
params (ndarray | None) – Parameter vector passed through to dVdx/d2Vdx2/V.

Returns:

q_t (float) – The trajectory’s final position.
K (complex) – The semiclassical propagator amplitude \(K(q_t,t;q_0,0)\).

See also

propagate_trajectory_monodromy_action: Supplies the trajectory, monodromy, and action.
herman_kluk_propagate_wavepacket: Sums many such trajectories into a full wavepacket propagator.

Examples

A free particle’s semiclassical propagator matches the exact quantum propagator \(\sqrt{m/(2\pi i\hbar t)}\exp[im(q_t-q_0)^2/(2\hbar t)]\) exactly:

>>> import numpy as np
>>> from numba import njit
>>> zero = njit(lambda q, params: 0.0, cache=False)
>>> q_t, K = van_vleck_propagator_1d(q0=0.0, p0=1.0, dVdx=zero, d2Vdx2=zero, V=zero, m=1.0, dt=1.0 / 2000, steps=2000)
>>> K_exact = np.sqrt(1.0 / (2j * np.pi * 1.0)) * np.exp(1j * (q_t - 0.0) ** 2 / 2.0)
>>> bool(abs(K - K_exact) < 1e-6)
True

A harmonic-oscillator trajectory that has not yet crossed a focal point (\(t<\pi/\omega\)) matches the exact Van Vleck (Mehler) propagator just as closely:

>>> m, omega = 1.0, 1.0
>>> params = np.array([m * omega ** 2])
>>> dVdx = njit(lambda q, params: params[0] * q, cache=False)
>>> d2Vdx2 = njit(lambda q, params: params[0], cache=False)
>>> V = njit(lambda q, params: 0.5 * params[0] * q ** 2, cache=False)
>>> q0, p0, t = 1.0, 0.3, 1.3
>>> q_t, K = van_vleck_propagator_1d(q0, p0, dVdx, d2Vdx2, V, m, dt=t / 4000, steps=4000, params=params)
>>> K_exact = (np.sqrt(m * omega / (2j * np.pi * np.sin(omega * t)))
...            * np.exp(1j * m * omega / (2 * np.sin(omega * t)) * ((q_t ** 2 + q0 ** 2) * np.cos(omega * t) - 2 * q_t * q0)))
>>> bool(abs(K - K_exact) < 1e-5)
True

The Gutzwiller trace formula: reconstructing a spectrum from classical periodic orbits alone.

Gutzwiller’s (1971) trace formula expresses the density of states \(g(E)=\sum_n\delta(E-E_n)\) as a sum over every classical periodic orbit of the system, with no reference to quantum eigenstates at all – each orbit contributes an oscillatory term whose frequency (in \(E\)) is set by its classical action and whose amplitude is set by its linear stability. For a bound, one-dimensional system this is not merely an approximation: the single family of periodic orbits (one per energy, since \(f=1\) degree of freedom always has exactly one oscillation) has an exact trace formula, obtained by Poisson-summing the discrete Einstein-Brillouin-Keller (EBK) spectrum of physicskit.semiclassical.core.wkb – gutzwiller_density_of_states() implements exactly that sum, reconstructing delta-function peaks at the Bohr-Sommerfeld energies purely from the classical action \(S(E)\) and period \(T(E)=dS/dE\).

The general Gutzwiller formula (for a system with genuinely isolated, unstable periodic orbits, as in chaotic scattering or 2D/3D billiards) instead weights each orbit by \(1/\sqrt{|2-\operatorname{tr}M|}\), its monodromy matrix’s instability – gutzwiller_amplitude_from_monodromy() implements that general ingredient, for use with orbits found in higher-dimensional chaotic systems (the stadium billiard scars analyzed in physicskit.semiclassical.systems.scarring live on exactly this kind of isolated unstable orbit).

The smooth (\(r=0\)) part of the trace formula is the Weyl term. For a two-dimensional billiard, balian_bloch_counting_function() and balian_bloch_level_density() give it with Balian and Bloch’s boundary, corner and curvature corrections.

physicskit.semiclassical.core.gutzwiller.balian_bloch_counting_function(k, area, perimeter, corner_angles=(), total_curvature=0.0, boundary='dirichlet')[source]#

Smoothed number of billiard levels below wavenumber \(k\), with boundary corrections.

Balian and Bloch (1970) expanded the smoothed counting function of the Helmholtz equation \((\nabla^2+k^2)\psi=0\) in a 2D domain in powers of \(1/k\):

\[\bar N(k) = \frac{A k^2}{4\pi} \mp \frac{L k}{4\pi} + C, \qquad C = \sum_i \frac{\pi^2-\theta_i^2}{24\pi\theta_i} + \frac{1}{12\pi}\oint\kappa\,ds.\]

The first term is Weyl’s law, the second is the boundary correction (minus for Dirichlet, plus for Neumann), and the constant comes from corners of interior angle \(\theta_i\) and from the curvature \(\kappa\) of the smooth parts of the boundary. Examples: a rectangle has \(C=4\times\tfrac1{16}=\tfrac14\), a disk has \(C=\tfrac16\).

Parameters:

k (float or ndarray) – Wavenumber(s) (\(E=\hbar^2k^2/2m\)).
area (float) – Billiard area \(A\).
perimeter (float) – Boundary length \(L\).
corner_angles (sequence of float, default=()) – Interior angles \(\theta_i\) of the boundary’s corners, in radians.
total_curvature (float) – \(\oint\kappa\,ds\) over the smooth parts of the boundary (\(2\pi\) for a disk or a stadium, 0 for a polygon).
boundary (str) – Boundary condition; sets the sign of the perimeter term.

Returns:

float or ndarray – \(\bar N(k)\), same shape as k.

See also

balian_bloch_level_density: Its derivative \(d\bar N/dk\).

References

R. Balian and C. Bloch, “Distribution of eigenfrequencies for the wave equation in a finite domain. I,” Ann. Phys. 60, 401-447 (1970); “… III. Eigenfrequency density oscillations,” Ann. Phys. 69, 76-160 (1972). H. P. Baltes and E. R. Hilf, Spectra of Finite Systems (Bibliographisches Institut, Mannheim, 1976).

Examples

The unit square with Dirichlet walls has levels \(k^2=\pi^2(n_x^2+n_y^2)\), \(n_x,n_y\ge1\). The exact staircase fluctuates about the expansion with zero mean:

>>> import numpy as np
>>> n = np.arange(1, 40)
>>> k_exact = np.sort(np.pi * np.hypot(n[:, None], n[None, :]).ravel())
>>> ks = np.linspace(10.0, 80.0, 2000)
>>> corners = [np.pi / 2] * 4
>>> residual = np.searchsorted(k_exact, ks) - balian_bloch_counting_function(ks, area=1.0, perimeter=4.0, corner_angles=corners)
>>> bool(abs(residual.mean()) < 0.05)
True

physicskit.semiclassical.core.gutzwiller.balian_bloch_level_density(k, area, perimeter, boundary='dirichlet')[source]#

Smoothed density of billiard levels per unit wavenumber, \(d\bar N/dk\).

\[\bar\rho(k) = \frac{A k}{2\pi} \mp \frac{L}{4\pi},\]

the derivative of balian_bloch_counting_function() (the corner and curvature constant drops out). In energy, with \(E=k^2\) (\(\hbar^2/2m=1\)), \(\bar\rho(E)=A/4\pi\mp L/(8\pi\sqrt E)\).

Parameters:

k (float or ndarray) – Wavenumber(s).
area (float) – Billiard area \(A\).
perimeter (float) – Boundary length \(L\).
boundary (str) – Boundary condition; sets the sign of the perimeter term.

Returns:

float or ndarray – \(\bar\rho(k)\), same shape as k.

Examples

>>> round(balian_bloch_level_density(10.0, area=np.pi, perimeter=2 * np.pi), 6)
4.5

physicskit.semiclassical.core.gutzwiller.classical_period(E, V, m, x_min, x_max, hbar=1.0, dE=0.0001)[source]#

Classical (round-trip) oscillation period \(T(E)=2\,dS/dE\) of a 1D bound orbit.

The action \(S(E)\) (wkb_action()) is the one-way traversal action, so its energy derivative \(dS/dE\) is the time to cross the classically allowed region once; the full period (there and back) is twice that – differentiated numerically here via a central difference.

Parameters:

E (float) – Energy at which to evaluate the period.
V (callable) – Potential energy function V(x).
m (float) – Particle mass.
x_min (float) – Search domain for turning_points().
x_max (float) – Search domain for turning_points().
hbar (float) – Value of \(\hbar\) to use (unused by the classical period itself; kept for a uniform signature with the rest of this module).
dE (float) – Step size for the central-difference derivative.

Return type:

Returns:

float – The classical period \(T(E)\).

See also

gutzwiller_density_of_states: Uses this as the trace formula’s smooth prefactor.

Examples

The harmonic oscillator’s period is famously independent of energy, \(T=2\pi/\omega\):

>>> V = lambda x: 0.5 * x ** 2
>>> round(classical_period(E=3.0, V=V, m=1.0, x_min=-20, x_max=20), 4)
6.2832

physicskit.semiclassical.core.gutzwiller.gutzwiller_amplitude_from_monodromy(M)[source]#

Gutzwiller stability amplitude \(1/\sqrt{|2-\operatorname{tr}M|}\) for an isolated periodic orbit.

In the full (multi-dimensional, generically chaotic) Gutzwiller trace formula, each isolated periodic orbit contributes with an amplitude set by how strongly nearby trajectories diverge from it over one period: an unstable (hyperbolic) orbit with monodromy eigenvalues \(\lambda,1/\lambda\) (\(|\lambda|>1\)) has \(\operatorname{tr}M=\lambda+1/\lambda\), so \(|2-\operatorname{tr}M|\) grows with the instability and the orbit’s contribution to the trace formula shrinks – highly unstable orbits matter less for the exact spectrum, but (as in physicskit.semiclassical.systems.scarring) can still leave a visible imprint on individual eigenstates.

Parameters:: M (ndarray) – Monodromy matrix of one period of the orbit, e.g. from physicskit.semiclassical.core.propagators.propagate_trajectory_monodromy_action().
Return type:: float
Returns:: float – The stability amplitude.

Examples

>>> import numpy as np
>>> M = np.array([[2.0, 0.0], [0.0, 0.5]])
>>> round(gutzwiller_amplitude_from_monodromy(M), 6)
1.414214

physicskit.semiclassical.core.gutzwiller.gutzwiller_density_of_states(E_grid, V, m, x_min, x_max, hbar=1.0, r_max=40, broadening=0.05)[source]#

Exact 1D trace-formula density of states, summed over repetitions of the single periodic orbit.

For a bound one-dimensional system, EBK quantization places levels exactly where \(S(E_n)/\hbar=(n+\tfrac12)\pi\). Poisson-summing the resulting delta comb \(g(E)=\sum_n\delta(E-E_n)\) over the integer \(n\) converts it into a sum over an integer \(r\) – physically, the \(r\)-fold repetition of the single primitive periodic orbit at each energy – giving the exact identity

\[g(E) = \frac{T(E)/2}{\pi\hbar}\left[1 + 2\sum_{r=1}^{\infty} \cos\!\left(\frac{2rS(E)}{\hbar} - r\pi\right)\right],\]

the one-dimensional Gutzwiller trace formula: the \(r=0\) term is the smooth Weyl density of states, and each \(r\ge1\) term is one repetition of the orbit, carrying the Maslov phase \(-r\pi\) (\(\sigma=2\) soft turning points per traversal, repeated \(r\) times). Truncating the sum at finite r_max with a convergence factor exp(-r*broadening) turns each delta function into a finite (Lorentzian-like) peak, suitable for numerical peak-finding.

Parameters:

E_grid (ndarray) – Energies at which to evaluate the density of states.
V (callable) – Potential energy function V(x).
m (float) – Particle mass.
x_min (float) – Search domain for turning_points().
x_max (float) – Search domain for turning_points().
hbar (float) – Value of \(\hbar\) to use.
r_max (int) – Number of orbit repetitions to sum.
broadening (float) – Per-repetition convergence factor; larger values broaden (and damp) each reconstructed peak.

Return type:

Returns:

ndarray – \(g(E)\), same shape as E_grid.

See also

classical_period: Supplies \(T(E)\).
physicskit.semiclassical.core.wkb.bohr_sommerfeld_energies: The exact peak locations this reconstructs.

Examples

The reconstructed peaks land exactly on the harmonic oscillator’s Bohr-Sommerfeld spectrum:

>>> import numpy as np
>>> from scipy.signal import find_peaks
>>> V = lambda x: 0.5 * x ** 2
>>> E_grid = np.linspace(0.2, 4.5, 600)
>>> dos = gutzwiller_density_of_states(E_grid, V, m=1.0, x_min=-20, x_max=20)
>>> peak_idx, _ = find_peaks(dos, height=0.3 * dos.max())
>>> np.round(E_grid[peak_idx], 1)
array([0.5, 1.5, 2.5, 3.5])

Bogomolny’s semiclassical transfer operator on a billiard boundary.

Bogomolny (1992) replaced the Gutzwiller sum over periodic orbits with a finite operator on a Poincare surface of section. Its kernel is the semiclassical propagator from one point of the section to the next,

\[T(q,q';E) = \frac{1}{\sqrt{2\pi i\hbar}} \sqrt{\left|\frac{\partial^2 S}{\partial q\,\partial q'}\right|}\, \exp\!\left[\frac{i}{\hbar}S(q,q';E) - i\frac{\pi\nu}{2}\right],\]

summed over the classical paths that go from \(q'\) to \(q\) without crossing the section in between. Quantized energies are the zeros of \(\det[1-T(E)]\). Expanding \(\log\det[1-T]\) in traces of powers of \(T\) gives back the Gutzwiller trace formula, but the determinant needs only short orbits (one crossing of the section) and a matrix whose size is set by the number of wavelengths across the section.

For a convex billiard the boundary itself is the natural section. With arc length \(s\) along it, the path from \(s'\) to \(s\) is the straight chord of length \(\ell\), \(S=\hbar k\ell\), and \(|\partial^2\ell/\partial s\,\partial s'|=\cos\theta\cos\theta'/\ell\), where \(\theta,\theta'\) are the angles between the chord and the inward normals. A Dirichlet wall adds a phase \(\pi\) per bounce (\(\nu=2\)), a Neumann wall none. bogomolny_transfer_operator() builds the discretized kernel and bogomolny_quantization_function() the quantity whose minima give the levels.

References

E. B. Bogomolny, “Semiclassical quantization of multidimensional systems,” Nonlinearity 5, 805-866 (1992).

physicskit.semiclassical.core.bogomolny.bogomolny_quantization_function(k_grid, points, normals, ds, boundary='dirichlet')[source]#

Smallest singular value of \(1-T(k)\) on a grid of wavenumbers.

Bogomolny’s quantization condition is \(\det[1-T(k)]=0\). The determinant of the discretized operator is never exactly zero, and its magnitude varies by orders of magnitude with \(k\); the smallest singular value of \(1-T(k)\) is a better-scaled measure of how close \(T\) is to having an eigenvalue 1. Its sharp minima are the semiclassical levels, and a near-degenerate pair shows up as a single minimum.

Parameters:

k_grid (ndarray) – Wavenumbers at which to evaluate.
points (ndarray) – Boundary discretization, as in bogomolny_transfer_operator().
normals (ndarray) – Boundary discretization, as in bogomolny_transfer_operator().
ds (ndarray) – Boundary discretization, as in bogomolny_transfer_operator().
boundary (str) – Wall boundary condition.

Return type:

Returns:

ndarray, shape matching k_grid – \(\sigma_{\min}[1-T(k)]\).

Examples

The lowest Dirichlet level of the unit disk is the first zero of \(J_0\), \(k=2.4048\):

>>> import numpy as np
>>> N = 120
>>> phi = 2 * np.pi * np.arange(N) / N
>>> pts = np.column_stack([np.cos(phi), np.sin(phi)])
>>> ks = np.linspace(2.2, 2.6, 81)
>>> f = bogomolny_quantization_function(ks, pts, -pts, np.full(N, 2 * np.pi / N))
>>> round(float(ks[np.argmin(f)]), 1)
2.4

physicskit.semiclassical.core.bogomolny.bogomolny_transfer_operator(points, normals, ds, k, boundary='dirichlet')[source]#

Discretized Bogomolny transfer operator \(T(k)\) on the boundary of a convex billiard.

\[T_{ij} = \sqrt{\frac{k\,\cos\theta_i\cos\theta_j}{2\pi i\,\ell_{ij}}}\; e^{ik\ell_{ij} - i\pi\nu/2}\,\sqrt{\Delta s_i\,\Delta s_j},\]

with \(\ell_{ij}\) the chord between boundary points \(i\) and \(j\). The symmetric weight \(\sqrt{\Delta s_i\Delta s_j}\) turns the integral kernel into a matrix whose eigenvalues approximate those of the operator. Use at least a few points per wavelength \(2\pi/k\); extra points only resolve the evanescent part of \(T\) better.

Parameters:

points (ndarray) – Boundary points, ordered along the boundary.
normals (ndarray) – Unit inward normals at points.
ds (ndarray) – Arc-length weight of each point.
k (float) – Wavenumber.
boundary (str) – Wall boundary condition (Maslov phase \(\nu=2\) or 0 per bounce).

Return type:

Returns:

ndarray of complex, shape (N, N) – The matrix \(T(k)\).

Raises:

ValueError – If boundary is not recognized.

Notes

Every chord must lie inside the billiard, so the boundary must be convex. Grazing chords (\(\cos\theta\to0\)) carry whispering-gallery states, which this leading-order kernel describes least accurately.

Examples

On a boundary of length \(L\) there are about \(kL/\pi\) propagating directions. \(T\) is close to unitary on them and close to zero on the rest; for the unit circle at \(k=20\) that is about 40 singular values near 1:

>>> import numpy as np
>>> N = 200
>>> phi = 2 * np.pi * np.arange(N) / N
>>> pts = np.column_stack([np.cos(phi), np.sin(phi)])
>>> T = bogomolny_transfer_operator(pts, -pts, np.full(N, 2 * np.pi / N), k=20.0)
>>> T.shape
(200, 200)
>>> sv = np.linalg.svd(T, compute_uv=False)
>>> int(np.sum(sv > 0.5))
41

Feynman path integrals and their reduction to classical mechanics as \(\hbar\to 0\).

Dirac’s 1933 paper “The Lagrangian in Quantum Mechanics” observed that the quantum propagator ought to be built from \(\exp(iS/\hbar)\), one factor per infinitesimal time step; Feynman (1948), and later Feynman & Hibbs (1965) in their textbook Quantum Mechanics and Path Integrals, turned that observation into a full reformulation of quantum mechanics: the amplitude to go from \((x_0,0)\) to \((x_f,T)\) is a sum over every continuous path \(x(t)\) connecting the two endpoints, each weighted by a phase built from the classical action functional \(S[x(t)]=\int_0^T L(x,\dot x)\,dt\) of that path,

\[K(x_f,T;x_0,0) = \int \mathcal{D}[x(t)]\; \exp\!\left(\frac{i}{\hbar}S[x(t)]\right).\]

This module implements that sum directly, by Monte Carlo sampling a large ensemble of candidate paths and adding up their phasors \(\exp(iS/\hbar)\) one at a time, ordered from the path closest to the true classical (stationary-action) trajectory outward. Two systems are exactly solvable end to end here – the free particle (free_particle_classical_path(), free_particle_classical_action()) and the harmonic oscillator (harmonic_oscillator_classical_path(), harmonic_oscillator_classical_action()) – so every piece of the machinery can be checked against a closed form.

The physical content is the emergence of the principle of stationary action itself. Near the classical path, \(S\) is stationary (\(\delta S=0\) to first order in the deformation), so nearby paths all have nearly the same action and their phasors stay nearly aligned: they interfere constructively. Far from the classical path, \(S\) varies rapidly from one path to a neighboring one, so \(\exp(iS/\hbar)\) spins rapidly through many full turns as the path is varied, and the contributions of an ensemble of such paths cancel almost completely by destructive interference (this is Feynman’s own “arrow and stopwatch” picture, reproduced here as the phasor spiral built by feynman_phasor_partial_sums()). Shrinking \(\hbar\) makes the phase \(S/\hbar\) more sensitive to path deformations, so the “nearly aligned” neighborhood of the classical path shrinks – in the \(\hbar\to 0\) limit only the classical path itself survives, and quantum mechanics’ sum over histories reduces exactly to the classical principle of stationary (least) action, \(\delta S=0\).

Paths here are represented by their positions at \(N+1\) uniformly spaced times over \([0,T]\) (\(N\) time slices, \(\mathrm{d}t=T/N\)), and the action of a discretized path is computed with the standard midpoint-rule time slicing,

\[S[x] \approx \sum_{i=0}^{N-1}\left[\frac{1}{2}m \left(\frac{x_{i+1}-x_i}{\mathrm{d}t}\right)^2 - V\!\left(\frac{x_i+x_{i+1}}{2}\right)\right]\mathrm{d}t \;\xrightarrow[N\to\infty]{}\; \int_0^T\!\left[\frac{1}{2}m\dot x^2 - V(x)\right]dt.\]

References

P. A. M. Dirac, “The Lagrangian in Quantum Mechanics,” Physikalische Zeitschrift der Sowjetunion 3, 64-72 (1933).

R. P. Feynman and A. R. Hibbs, Quantum Mechanics and Path Integrals, McGraw-Hill (1965).

physicskit.semiclassical.core.path_integral.build_phasor_diagram(x0, xf, T, m, hbar, potential='free', omega=None, n_slices=60, n_paths=400, sigma=0.5, seed=None)[source]#

Build a complete Feynman phasor-spiral dataset for the free particle or harmonic oscillator.

Ties together the classical trajectory/action, a random ensemble of candidate paths (sample_random_paths()), their discretized actions (discretized_action()), and the resulting phasor partial sums (feynman_phasor_partial_sums()) – everything animate_feynman_phasor_spiral() needs to draw. The exact classical path is explicitly included as one member of the returned ensemble, since it is the single most important path in the whole construction.

Parameters:

x0 (float) – Fixed initial and final positions.
xf (float) – Fixed initial and final positions.
T (float) – Total elapsed time (must be positive).
m (float) – Particle mass.
hbar (float) – Reduced Planck constant (or an effective value thereof); must be nonzero. Smaller hbar concentrates the phasor sum more tightly around the classical path.
potential (str) – Which exactly-solvable system to use.
omega (float | None) – Angular frequency; required if potential="harmonic".
n_slices (int) – Number of time slices per path.
n_paths (int) – Number of random paths to sample (in addition to the classical path itself).
sigma (float) – Amplitude of the random path wiggles, passed to sample_random_paths().
seed (int | None) – Seed for reproducible path sampling.

Returns:

paths (ndarray, shape (n_paths + 1, n_slices + 1)) – The random ensemble with the classical path appended as the last row.
actions (ndarray, shape (n_paths + 1,)) – Discretized action of every path in paths (including the classical path’s discretized action, which need not exactly equal classical_action at finite n_slices).
classical_path (ndarray, shape (n_slices + 1,)) – The exact classical trajectory.
classical_action (float) – The exact (closed-form) classical action.
partial_sums (ndarray of complex, shape (n_paths + 2,)) – Phasor partial sums from feynman_phasor_partial_sums(), ordered closest-to-classical-action first.

Raises:

ValueError – If potential is not one of "free"/"harmonic", or if potential="harmonic" and omega is not given.

Examples

>>> paths, actions, x_cl, S_cl, partial = build_phasor_diagram(
...     x0=0.0, xf=1.0, T=1.0, m=1.0, hbar=0.05, potential="free",
...     n_slices=20, n_paths=50, sigma=0.3, seed=0)
>>> paths.shape[0] == actions.shape[0] == 51
True
>>> partial.shape
(52,)

physicskit.semiclassical.core.path_integral.discretized_action(x_array, dt, m, V)[source]#

Midpoint-rule discretized action of a single time-sliced path.

\[S[x] = \sum_{i=0}^{N-1}\left[\frac{1}{2}m \left(\frac{x_{i+1}-x_i}{\mathrm{d}t}\right)^2 - V\!\left(\frac{x_i+x_{i+1}}{2}\right)\right]\mathrm{d}t.\]

Parameters:

x_array (ndarray) – Path positions \(x_0,\dots,x_N\) at \(N+1\) uniformly spaced times.
dt (float) – Time-slice width.
m (float) – Particle mass.
V (Callable[[ndarray], ndarray]) – Potential energy function V(x) -> array_like, evaluated at the midpoint of each slice; must accept a NumPy array.

Return type:

Returns:

float – The discretized action \(S[x]\).

Examples

>>> import numpy as np
>>> x = np.array([0.0, 0.5, 1.0, 1.5, 2.0])
>>> V = lambda x: 0.0 * x
>>> round(discretized_action(x, dt=0.25, m=1.0, V=V), 6)
2.0

physicskit.semiclassical.core.path_integral.feynman_phasor_partial_sums(paths, actions, hbar, classical_action=None)[source]#

Feynman’s phasor (“arrow and stopwatch”) partial sums, ordered closest-to-classical first.

Sorts the given ensemble of paths by \(|S-S_{\rm cl}|\) (ascending) and accumulates the running sum of unit phasors \(\exp(iS/\hbar)\), one path at a time, in that order. Because \(S\) is stationary at the classical path, paths near it have nearly equal action and their phasors reinforce (the sum’s magnitude grows quickly at first); paths far from it have rapidly varying action and their phasors point in essentially random directions, so later additions contribute a slowly-converging random walk instead of steady growth. This is exactly the data behind Feynman’s arrow-and-stopwatch phasor-spiral diagram.

Parameters:

paths (ndarray) – Ensemble of paths, used only to check that actions has one entry per path; the ordering and summation depend only on actions.
actions (ndarray) – Discretized action of each path, e.g. from discretized_action().
hbar (float) – Reduced Planck constant (or an effective value thereof); must be nonzero.
classical_action (float | None) – Reference action \(S_{\rm cl}\) used to sort paths by proximity. If None, the minimum action found in actions is used instead.

Return type:

Returns:

ndarray of complex, shape (n_paths + 1,) – Partial sums \(\sum_{j\le n}\exp(iS_j/\hbar)\), starting from 0 (index 0) up to the full sum over all paths (index n_paths), accumulated in closest-to-classical-first order.

Raises:

ValueError – If hbar == 0, actions is empty, or len(paths) != len(actions).

Examples

>>> import numpy as np
>>> actions = np.array([0.0, 0.1, -3.0])
>>> paths = np.zeros((3, 2))
>>> partial = feynman_phasor_partial_sums(paths, actions, hbar=1.0, classical_action=0.0)
>>> partial.shape
(4,)
>>> bool(np.isclose(partial[0], 0.0))
True
>>> bool(np.isclose(partial[-1], np.sum(np.exp(1j * actions))))
True

physicskit.semiclassical.core.path_integral.free_particle_classical_action(x0, xf, T, m)[source]#

Exact classical (Hamilton’s principal function) action of a free particle.

\[S_{\rm cl} = \frac{m(x_f-x_0)^2}{2T}.\]

Parameters:

x0 (float) – Initial and final positions.
xf (float) – Initial and final positions.
T (float) – Total elapsed time (must be positive).
m (float) – Particle mass.

Return type:

Returns:

float – The classical action \(S_{\rm cl}\).

Raises:

ValueError – If T <= 0.

Examples

>>> free_particle_classical_action(x0=0.0, xf=2.0, T=1.0, m=1.0)
2.0

physicskit.semiclassical.core.path_integral.free_particle_classical_path(x0, xf, T, m, n_slices)[source]#

Exact classical trajectory of a free particle between fixed endpoints.

With \(V=0\) the Euler-Lagrange equation \(\ddot x=0\) has the unique solution connecting \((x_0,0)\) and \((x_f,T)\): the straight line \(x_{\rm cl}(t) = x_0 + (x_f-x_0)\,t/T\).

Parameters:

x0 (float) – Initial position, at \(t=0\).
xf (float) – Final position, at \(t=T\).
T (float) – Total elapsed time (must be positive).
m (float) – Particle mass (unused for the free-particle trajectory itself, kept for signature symmetry with discretized_action() and the harmonic-oscillator counterpart).
n_slices (int) – Number of time slices \(N\); the returned arrays have \(N+1\) points.

Returns:

t_array (ndarray, shape (n_slices + 1,)) – Uniformly spaced times over \([0,T]\).
x_array (ndarray, shape (n_slices + 1,)) – Classical trajectory \(x_{\rm cl}(t)\).

Raises:

ValueError – If T <= 0 or n_slices < 1.

Examples

>>> t, x = free_particle_classical_path(x0=0.0, xf=2.0, T=1.0, m=1.0, n_slices=4)
>>> x
array([0. , 0.5, 1. , 1.5, 2. ])

physicskit.semiclassical.core.path_integral.harmonic_oscillator_classical_action(x0, xf, T, m, omega)[source]#

Exact classical (Hamilton’s principal function) action of a harmonic oscillator.

\[S_{\rm cl} = \frac{m\omega}{2\sin(\omega T)} \left[(x_0^2+x_f^2)\cos(\omega T) - 2x_0 x_f\right].\]

Parameters:

x0 (float) – Initial and final positions.
xf (float) – Initial and final positions.
T (float) – Total elapsed time (must be positive).
m (float) – Particle mass.
omega (float) – Angular frequency of the oscillator.

Return type:

Returns:

float – The classical action \(S_{\rm cl}\).

Raises:

ValueError – If T <= 0 or omega*T is (numerically) a multiple of \(\pi\) (a conjugate point; see harmonic_oscillator_classical_path()).

Examples

>>> round(harmonic_oscillator_classical_action(x0=0.0, xf=0.0, T=1.0, m=1.0, omega=0.5), 6)
0.0

physicskit.semiclassical.core.path_integral.harmonic_oscillator_classical_path(x0, xf, T, m, omega, n_slices)[source]#

Exact classical trajectory of a harmonic oscillator between fixed endpoints.

The Euler-Lagrange equation \(\ddot x = -\omega^2 x\) has the unique solution matching both endpoints

\[x_{\rm cl}(t) = \frac{x_0\sin(\omega(T-t)) + x_f\sin(\omega t)} {\sin(\omega T)},\]

valid as long as \(\sin(\omega T)\neq 0\); at a conjugate point (\(\omega T\) an integer multiple of \(\pi\)) every path through both endpoints has the same action to first order and the classical path is not unique, so this raises instead of silently dividing by zero.

Parameters:

x0 (float) – Initial and final positions.
xf (float) – Initial and final positions.
T (float) – Total elapsed time (must be positive).
m (float) – Particle mass (unused for the trajectory itself, kept for signature symmetry with harmonic_oscillator_classical_action()).
omega (float) – Angular frequency of the oscillator.
n_slices (int) – Number of time slices \(N\); the returned arrays have \(N+1\) points.

Returns:

t_array (ndarray, shape (n_slices + 1,)) – Uniformly spaced times over \([0,T]\).
x_array (ndarray, shape (n_slices + 1,)) – Classical trajectory \(x_{\rm cl}(t)\).

Raises:

ValueError – If T <= 0, n_slices < 1, or omega*T is (numerically) a multiple of \(\pi\) (a conjugate point).

Examples

>>> t, x = harmonic_oscillator_classical_path(x0=0.0, xf=1.0, T=1.0, m=1.0, omega=1.0, n_slices=4)
>>> bool(np.isclose(x[0], 0.0)) and bool(np.isclose(x[-1], 1.0))
True

physicskit.semiclassical.core.path_integral.sample_random_paths(x0, xf, T, n_slices, n_paths, sigma, n_modes=8, seed=None)[source]#

Sample random paths between fixed endpoints as sine-series Brownian bridges.

Each path is the straight-line interpolation between the endpoints plus a random wiggle built from a discrete sine series,

\[x_i = x_0 + (x_f-x_0)\frac{i}{N} + \sigma\sum_{k=1}^{K} \frac{a_k}{k}\sin\!\left(\frac{k\pi i}{N}\right), \qquad a_k\sim\mathcal{N}(0,1)\ \text{i.i.d. per path},\]

where \(N\) is n_slices and \(K\) is n_modes. Every term of the sine series vanishes at \(i=0\) and \(i=N\) identically, so every sampled path exactly satisfies the fixed endpoint conditions \(x_0\) and \(x_f\) regardless of the random coefficients drawn – this is what makes it a valid “Brownian bridge” construction. Paths are built around the straight line (not around any particular classical trajectory), so the ensemble generically explores the space of continuous paths between the endpoints, and any special concentration of the phasor sum near the true classical path is a genuine result rather than built in by construction.

Parameters:

x0 (float) – Fixed initial and final positions shared by every sampled path.
xf (float) – Fixed initial and final positions shared by every sampled path.
T (float) – Total elapsed time (must be positive).
n_slices (int) – Number of time slices \(N\); each path has \(N+1\) points.
n_paths (int) – Number of random paths to draw.
sigma (float) – Overall amplitude of the random wiggle.
n_modes (int) – Number of sine modes \(K\) used to build each path’s wiggle; higher modes are automatically suppressed by the \(1/k\) weighting, so this mainly controls the finest wiggle scale resolved.
seed (int | None) – Seed for the random number generator, for reproducibility.

Return type:

Returns:

ndarray, shape (n_paths, n_slices + 1) – Sampled paths, each row a full path \(x_0,\dots,x_N\).

Examples

>>> paths = sample_random_paths(x0=0.0, xf=1.0, T=1.0, n_slices=10, n_paths=5, sigma=0.3, seed=0)
>>> paths.shape
(5, 11)
>>> bool(np.allclose(paths[:, 0], 0.0)) and bool(np.allclose(paths[:, -1], 1.0))
True

Quantum scars: eigenstate density concentrated on unstable classical periodic orbits.

Heller (1984) discovered that individual eigenstates of a classically chaotic system are not always the featureless, ergodically-spread blobs random-matrix intuition suggests: a small but robust fraction instead show enhanced probability density in a tube around one particular unstable classical periodic orbit – a scar – even though almost every orbit nearby diverges from it exponentially fast. The bridge to physicskit.semiclassical.core.gutzwiller is direct: the same unstable periodic orbits that weight the Gutzwiller trace formula’s oscillating sum are exactly the orbits that scar individual eigenstates, since eigenstates are (schematically) superpositions of trace-formula terms.

bouncing_ball_energies() and bouncing_ball_orbit_points() concern the single most famous scarred family in the Bunimovich stadium billiard (physicskit.quantum.chapters.potentials.StadiumBilliard2D): the “bouncing ball” orbits bouncing straight up and down between the flat top and bottom walls, which – unlike a generic stadium orbit – are only marginally unstable, so many low-lying eigenstates concentrate along them. scar_enhancement() quantifies that concentration for any eigenstate density, and husimi_projection_1d() builds a non-periodic analog of physicskit.chaos.quantum.husimi.husimi_function() for viewing a 1D wavefunction slice (e.g. along the billiard’s flat wall) in phase space.

physicskit.semiclassical.systems.scarring.bouncing_ball_energies(R, hbar=1.0, m=1.0, n_max=6)[source]#

Predicted energies of the “bouncing ball” orbit family in a stadium billiard of cap radius R.

A trajectory launched perpendicular to the flat top and bottom walls of a Bunimovich stadium (see physicskit.quantum.chapters.potentials.StadiumBilliard2D) bounces straight up and down forever, blind to the length L of the central rectangle – exactly the motion of a particle in a 1D infinite square well of width \(2R\). Quantizing that 1D motion gives a leading-order prediction for where “bouncing ball” states (eigenstates concentrated on this orbit family) appear in the full 2D spectrum:

\[E_n = \frac{(n\pi\hbar)^2}{2m(2R)^2}, \qquad n=1,2,3,\dots\]

Parameters:

R (float) – Radius of the stadium’s semicircular end-caps (half-height of the billiard).
hbar (float) – Value of \(\hbar\) to use.
m (float) – Particle mass.
n_max (int) – Number of levels to compute.

Return type:

Returns:

ndarray, shape (n_max,) – Predicted bouncing-ball energies, ascending.

See also

bouncing_ball_orbit_points: The classical orbit these energies quantize.

Examples

This is exactly the ordinary infinite-square-well spectrum for a well of width \(2R\):

>>> import numpy as np
>>> R = 0.5
>>> energies = bouncing_ball_energies(R, n_max=3)
>>> exact = (np.arange(1, 4) * np.pi) ** 2 / (2 * (2 * R) ** 2)
>>> bool(np.allclose(energies, exact))
True

physicskit.semiclassical.systems.scarring.bouncing_ball_orbit_points(x0, R, n_bounces=4)[source]#

Trace the classical “bouncing ball” orbit at horizontal position x0 in a stadium billiard.

A perpendicular launch from the bottom wall at \((x_0,-R)\) reflects straight back and forth off the flat top (\(y=R\)) and bottom (\(y=-R\)) walls, tracing the same vertical segment forever – returned here as a zig-zag list of points suitable for overlaying on a density plot.

Parameters:

x0 (float) – Horizontal position of the orbit, with \(|x_0| < L/2\) (strictly inside the flat central section of the stadium).
R (float) – Radius of the stadium’s semicircular end-caps.
n_bounces (int) – Number of full up-down traversals to trace.

Return type:

Returns:

x, y (ndarray) – Coordinates tracing the orbit, shape (2 * n_bounces + 1,).

See also

bouncing_ball_energies: Quantized energies of this orbit family.

Examples

>>> x, y = bouncing_ball_orbit_points(x0=0.2, R=0.5, n_bounces=2)
>>> x
array([0.2, 0.2, 0.2, 0.2, 0.2])
>>> y
array([-0.5,  0.5, -0.5,  0.5, -0.5])

physicskit.semiclassical.systems.scarring.husimi_projection_1d(psi, s, hbar=1.0, sigma=None, resolution=60, s0_range=None, p0_range=None)[source]#

Husimi (coherent-state) phase-space projection of a 1D wavefunction slice on an open interval.

The non-periodic counterpart of physicskit.chaos.quantum.husimi.husimi_function(): overlaps psi with ordinary (non-periodized) coherent states \(g_{s_0,p_0}(s)=\exp[-(s-s_0)^2/2\sigma^2+ip_0(s-s_0)/\hbar]\) at every point of an \((s_0,p_0)\) grid, suited to a slice taken along an open boundary – e.g. the flat top wall of a stadium billiard, where bouncing_ball_orbit_points() lives – rather than a periodic domain.

Parameters:

psi (ndarray) – Wavefunction values sampled on s.
s (ndarray) – Uniform 1D coordinate grid psi is sampled on.
hbar (float) – Value of \(\hbar\) to use.
sigma (float | None) – Coherent-state position width. Defaults to \(\sqrt{\hbar}\) (the minimum-uncertainty, equal-spread-in-natural-units choice).
resolution (int) – Number of grid points along each of the \(s_0\), \(p_0\) axes.
s0_range (tuple | None) – Range of \(s_0\) to scan. Defaults to (s.min(), s.max()).
p0_range (tuple | None) – Range of \(p_0\) to scan. Defaults to \(\pm\pi\hbar/\Delta s\) (the grid’s Nyquist momentum).

Return type:

Returns:

S0, P0 (ndarray, shape (resolution, resolution)) – Phase-space grid (indexing="ij").
husimi (ndarray, shape (resolution, resolution)) – The Husimi distribution, normalized to a peak value of 1.

See also

scar_enhancement: A simpler, single-number scarring diagnostic in position space rather than phase space.

Examples

A Gaussian wave packet’s Husimi projection peaks at its own position and momentum:

>>> import numpy as np
>>> s = np.linspace(-10, 10, 2000)
>>> s0_true, p0_true, w = 2.0, 3.0, 1.0
>>> psi = np.exp(-(s - s0_true) ** 2 / (2 * w ** 2)) * np.exp(1j * p0_true * s / 1.0)
>>> S0, P0, H = husimi_projection_1d(psi, s, hbar=1.0, sigma=w, resolution=80, s0_range=(-2, 6), p0_range=(-2, 8))
>>> i, j = np.unravel_index(np.argmax(H), H.shape)
>>> bool(abs(S0[i, j] - s0_true) < 0.2 and abs(P0[i, j] - p0_true) < 0.2)
True

physicskit.semiclassical.systems.scarring.scar_enhancement(density, X, Y, mask, x0, half_width)[source]#

Density enhancement in a tube around a vertical bouncing-ball orbit, relative to the billiard average.

\[\eta = \frac{\langle|\psi|^2\rangle_{\text{tube}}}{\langle|\psi|^2\rangle_{\text{billiard}}}, \qquad \text{tube} = \{(x,y)\in\text{billiard} : |x-x_0|\le w\},\]

a simple, direct measure of scarring: \(\eta\approx1\) for an ergodically-spread (unscarred) state, since the tube then just samples a representative fraction of the whole density, while \(\eta\gg1\) signals density piled up specifically along the orbit at \(x_0\).

Parameters:

density (ndarray) – Probability density \(|\psi(x,y)|^2\) on a grid.
X (ndarray) – Coordinate meshgrid matching density.
Y (ndarray) – Coordinate meshgrid matching density.
mask (ndarray) – True inside the billiard (e.g. from physicskit.quantum.chapters.potentials.StadiumBilliard2D.mask()).
x0 (float) – Horizontal position of the bouncing-ball orbit to test.
half_width (float) – Half-width of the tube around x0.

Return type: