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

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.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.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().

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.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.

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.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).

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])

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: