physicskit.quantum#

physicskit.quantum hub · History · Examples · API reference

physicskit.quantum: a visual and computational tour of quantum mechanics.

  • physicskit.quantum.core.eigensolvers – the Numerov shooting-method bound-state solver, plus ready-made potential wells (infinite, finite, asymmetric step, linear-gravitational, double, harmonic).

  • physicskit.quantum.core.operators – spin, ladder, position, and momentum operators as explicit matrices, the Pauli matrices, and generic commutator/expectation-value/Hermiticity utilities.

  • physicskit.quantum.core.solvers – split-operator time propagation for the time-dependent Schrodinger equation, in 1D and 2D.

  • physicskit.quantum.chapters.wave_packets – free and dispersing Gaussian wave packets, the double-slit experiment, and quantum revivals.

  • physicskit.quantum.chapters.potentials – asymmetric/step wells, the quantum bouncer, the double well, the finite square well, and 2D quantum boxes (rectangle, circular dot, stadium billiard).

  • physicskit.quantum.chapters.hydrogen_am – the hydrogen atom’s exact radial and angular (spherical harmonic) wavefunctions.

  • physicskit.quantum.chapters.harmonic_spin – the quantum harmonic oscillator (ladder operators, coherent states) and thermal (mixed) states.

  • physicskit.quantum.chapters.perturbation – the Zeeman and Stark effects, and a Floquet-driven infinite square well.

  • physicskit.quantum.chapters.entanglement – Bell states and CHSH correlations, and the Aharonov-Bohm ring.

See physicskit.semiclassical for WKB/EBK quantization, Van Vleck/Herman-Kluk semiclassical propagators, the Gutzwiller trace formula, and quantum scarring.

class physicskit.quantum.AharonovBohmRing(R=1.0, m=1.0, hbar=1.0, q=1.0)[source]#

Bases: object

A charged particle confined to a 1D ring threading a magnetic flux.

\(B=0\) everywhere on the ring itself – only the vector potential \(A = \Phi/(2\pi R)\) is nonzero there. The eigenspectrum

\[E_n(\Phi) = \frac{\hbar^2}{2mR^2}\left(n - \frac{\Phi}{\Phi_0}\right)^2\]

shows the Aharonov-Bohm phase shift \(\Delta\phi = (q/\hbar)\oint \mathbf A\cdot d\mathbf l = 2\pi\Phi/\Phi_0\) purely through boundary conditions, without the particle ever entering a region where \(B\neq0\).

Parameters:

  • R (float) – Ring radius.

  • m (float) – Particle mass.

  • hbar (float) – Value of \(\hbar\) to use.

  • q (float) – Particle charge.

R: float = 1.0#
aharonov_bohm_phase(Phi)[source]#

The Aharonov-Bohm phase \(\Delta\phi = 2\pi\Phi/\Phi_0\).

The phase accumulated by a particle encircling the flux tube once.

Parameters:

Phi (float) – Enclosed magnetic flux.

Return type:

float

Returns:

float

eigenstate(n, phi)[source]#

Real-space ring wavefunction \(\psi_n(\phi) = e^{in\phi}/\sqrt{2\pi R}\).

Parameters:

  • n (int) – Angular-momentum quantum number.

  • phi (ndarray) – Angular positions around the ring.

Return type:

ndarray

Returns:

numpy.ndarray – Complex-valued.

energy(n, Phi)[source]#

Eigenenergy \(E_n(\Phi)\).

Parameters:

  • n (int or array_like) – Angular-momentum quantum number(s).

  • Phi (float) – Enclosed magnetic flux.

Return type:

ndarray

Returns:

numpy.ndarray

property flux_quantum: float#

The flux quantum \(\Phi_0 = 2\pi\hbar/q\) (\(h/q\) in these units).

Type:

float

hbar: float = 1.0#
m: float = 1.0#
persistent_current(n, Phi, dPhi=1e-06)[source]#

Equilibrium persistent current carried by level \(n\).

\[I_n(\Phi) = -\frac{dE_n}{d\Phi},\]

periodic in \(\Phi\) with period \(\Phi_0\) – the direct experimental signature of the Aharonov-Bohm phase.

Parameters:

  • n (int) – Angular-momentum quantum number.

  • Phi (float) – Enclosed magnetic flux.

  • dPhi (float) – Step size for the central-difference derivative.

Return type:

float

Returns:

float

q: float = 1.0#
spectrum(Phi, n_range=5)[source]#

Sorted energies for \(n \in [-n_\text{range}, n_\text{range}]\).

Parameters:

  • Phi (float) – Enclosed magnetic flux.

  • n_range (int) – Range of angular-momentum quantum numbers to include.

Return type:

ndarray

Returns:

numpy.ndarray

class physicskit.quantum.BellCorrelations(state=None)[source]#

Bases: object

EPR correlations for the spin-1/2 singlet state.

\[\lvert\psi^-\rangle = \frac{\lvert01\rangle - \lvert10\rangle}{\sqrt2}.\]

Reproduces the textbook result \(E(a,b) = -\cos(\theta_a - \theta_b)\) and the CHSH violation of the classical (local hidden-variable) bound \(\lvert S\rvert \le 2\), up to the quantum (Tsirelson) bound \(2\sqrt2\).

Parameters:

state (ndarray) – The two-qubit state to use; defaults to the singlet bell_state()('psi-').

chsh_S(a, a_prime, b, b_prime)[source]#

The CHSH combination \(S = E(a,b) - E(a,b') + E(a',b) + E(a',b')\).

Parameters:

  • a (float) – The two measurement angles for qubit 1.

  • a_prime (float) – The two measurement angles for qubit 1.

  • b (float) – The two measurement angles for qubit 2.

  • b_prime (float) – The two measurement angles for qubit 2.

Return type:

float

Returns:

float

chsh_optimal()[source]#

\(\lvert S\rvert\) at the standard optimal angles \((0, \pi/2, \pi/4, 3\pi/4)\).

Gives the Tsirelson bound \(\lvert S\rvert = 2\sqrt2 \approx 2.828\), violating the classical bound of 2 predicted by any local hidden-variable theory.

Return type:

float

Returns:

float

correlation(theta_a, theta_b)[source]#

The spin correlation \(E(a,b) = \langle\psi\rvert \hat A(\theta_a) \otimes \hat B(\theta_b) \lvert\psi\rangle\).

Parameters:

  • theta_a (float) – Measurement angle for qubit 1 (from \(z\), in the x-z plane).

  • theta_b (float) – Measurement angle for qubit 2.

Return type:

float

Returns:

float

measurement_probabilities(theta_a, theta_b)[source]#

Joint outcome probabilities for projective spin measurements.

Parameters:

  • theta_a (float) – Measurement angle for qubit 1.

  • theta_b (float) – Measurement angle for qubit 2.

Return type:

dict

Returns:

dict – Maps '++', '+-', '-+', '--' to their probabilities.

monte_carlo_correlation(theta_a, theta_b, n_trials=20000, rng=None)[source]#

Simulate joint measurements and return the empirical correlation.

Samples simulated joint measurement outcomes from the exact quantum joint distribution and returns the empirical correlation \(\langle AB\rangle\), statistically converging to correlation()(theta_a, theta_b).

Parameters:

  • theta_a (float) – Measurement angle for qubit 1.

  • theta_b (float) – Measurement angle for qubit 2.

  • n_trials (int) – Number of simulated measurement pairs.

  • rng (Generator | None) – Random number generator; a fresh default one is used if omitted.

Return type:

float

Returns:

float

state: ndarray = None#
class physicskit.quantum.Box2DEigenstate(nx, ny, energy)[source]#

Bases: object

One eigenstate of RectangularBox2D.

Parameters:
energy: float#

The eigenenergy \(E_{n_x,n_y}\).

Type:

float

nx: int#

Quantum number along \(x\) (\(\ge 1\)).

Type:

int

ny: int#

Quantum number along \(y\) (\(\ge 1\)).

Type:

int

psi(X, Y, Lx, Ly)[source]#

Evaluate the eigenfunction on a grid.

\[\psi_{n_x,n_y}(x,y) = \sqrt{\frac{4}{L_x L_y}} \sin\!\left(\frac{n_x \pi x}{L_x}\right) \sin\!\left(\frac{n_y \pi y}{L_y}\right).\]
Parameters:
Return type:

ndarray

Returns:

numpy.ndarray – \(\psi_{n_x,n_y}(x,y)\) on the grid.

class physicskit.quantum.CircularBox2D(R=1.0, hbar=1.0, m_mass=1.0)[source]#

Bases: object

Circular infinite well (quantum dot) of radius \(R\).

Eigenstates are Bessel functions \(\psi_{mn}(r,\phi) \propto J_m(k_{mn} r)\, e^{im\phi}\), vanishing at \(r=R\), showing nodal rings (radial quantum number \(n\)) and angular momentum quantization (\(m\)).

Parameters:

  • R (float) – Dot radius.

  • hbar (float) – Value of \(\hbar\) to use.

  • m_mass (float) – Particle mass (named to avoid clashing with the angular momentum quantum number m).

eigenstate(m, n)[source]#

Build the \((m,n)\) eigenstate.

Parameters:

  • m (int) – Angular momentum quantum number.

  • n (int) – Radial index (\(n=1,2,3,\dots\)); the state uses the \(n\)-th positive zero of the order-\(\lvert m\rvert\) Bessel function.

Return type:

CircularEigenstate

Returns:

CircularEigenstate

polar_grid(n_r=150, n_phi=150)[source]#

A polar meshgrid spanning the dot.

Parameters:

  • n_r (int) – Number of radial points.

  • n_phi (int) – Number of angular points.

Returns:

r, phi (numpy.ndarray) – Meshgrid arrays (indexing='ij').

spectrum(m_max=4, n_max=4)[source]#

All eigenstates with \(\lvert m\rvert \le m_\text{max}\), \(1 \le n \le n_\text{max}\), sorted by energy.

Parameters:

  • m_max (int) – Maximum \(\lvert m \rvert\).

  • n_max (int) – Maximum radial index.

Return type:

list

Returns:

list of CircularEigenstate

class physicskit.quantum.CircularEigenstate(m, n, k, energy)[source]#

Bases: object

One eigenstate of CircularBox2D.

Parameters:
energy: float#

The eigenenergy \(E_{mn} = \hbar^2 k_{mn}^2 / 2m\).

Type:

float

k: float#

Wavenumber \(k_{mn} = z_{mn}/R\).

Type:

float

m: int#

Angular momentum quantum number.

Type:

int

n: int#

Radial quantum number (the state uses the \(n\)-th positive zero of \(J_{\lvert m\rvert}\)).

Type:

int

psi(r, phi, R)[source]#

Evaluate the eigenfunction on a polar grid.

\[\psi_{mn}(r,\phi) = N_{mn}\, J_m(k_{mn} r)\, e^{i m \phi}.\]
Parameters:
Return type:

ndarray

Returns:

numpy.ndarray – Complex-valued \(\psi_{mn}(r,\phi)\).

class physicskit.quantum.DoubleWellResult(x, energies, wavefunctions, splitting, tunneling_period)[source]#

Bases: object

Container for the output of DoubleWellSimulator.solve().

Parameters:
energies: ndarray#

Eigenenergies \(E_n\), ascending.

Type:

numpy.ndarray

splitting: float#

The symmetric/antisymmetric doublet gap \(\Delta E = E_1 - E_0\).

Type:

float

tunneling_period: float#

The tunneling oscillation period \(T = 2\pi\hbar/\Delta E\).

Type:

float

wavefunctions: ndarray#

Normalized eigenfunctions, shape (n_states, len(x)).

Type:

numpy.ndarray

x: ndarray#

The spatial grid.

Type:

numpy.ndarray

class physicskit.quantum.DoubleWellSimulator(lam=0.5, a=2.5, x_extent=6.0, n_points=1600, hbar=1.0, m=1.0)[source]#

Bases: object

Symmetric double well \(V(x) = \lambda (x^2 - a^2)^2\).

The two lowest eigenstates form a nearly degenerate symmetric/ antisymmetric doublet whose splitting \(\Delta E = E_1 - E_0\) sets the period \(T=2\pi\hbar/\Delta E\) of left-right tunneling oscillations of a wavepacket initially localized in one well.

Parameters:

  • lam (float) – Overall strength \(\lambda\) of the quartic potential.

  • a (float) – Half-separation of the two minima.

  • x_extent (float) – Half-width of the numerical domain, x in [-x_extent, x_extent].

  • n_points (int) – Number of grid points.

  • hbar (float) – Value of \(\hbar\) to use.

  • m (float) – Particle mass.

left_well_probability(result, t, side='left')[source]#

Probability of finding the particle in the left well (\(x<0\)) vs. time.

\[P(t) = \int_{-\infty}^0 \lvert\psi(x,t)\rvert^2\, dx,\]

oscillating with the tunneling period \(T=2\pi\hbar/\Delta E\).

Parameters:
Return type:

ndarray

Returns:

numpy.ndarray – \(P(t)\), shape (len(t),).

localized_state(result, side='left')[source]#

Build a state localized in one well from the doublet.

\[\psi(x) = \frac{\psi_0(x) \pm \psi_1(x)}{\sqrt 2},\]

with + for the left well and - for the right well.

Parameters:

  • result (DoubleWellResult) – Output of solve() (must include at least 2 states).

  • side (str) – Which well to localize the state in.

Return type:

ndarray

Returns:

numpy.ndarray – The localized wavefunction at \(t=0\).

solve(n_states=6)[source]#

Compute the eigenstates and the tunneling splitting/period.

Parameters:

n_states (int) – Number of bound states to compute (at least 2, for the doublet).

Return type:

DoubleWellResult

Returns:

DoubleWellResult

tunneling_oscillation(result, t, side='left')[source]#

Time-dependent probability density of a state localized on side.

See tunneling_wavefunction().

Parameters:
Return type:

ndarray

Returns:

numpy.ndarray – \(\lvert\psi(x,t)\rvert^2\), shape (len(t), len(x)).

tunneling_wavefunction(result, t, side='left')[source]#

Time-dependent wavefunction of a state localized on side.

Evolves \(\psi(x,0) = (\psi_0 \pm \psi_1)/\sqrt2\) under the two doublet energies; \(\lvert\psi(x,t)\rvert^2\) (tunneling_oscillation()) periodically tunnels to the other well with period \(T=2\pi\hbar/\Delta E\). Shaped for animate_density().

Parameters:
Return type:

ndarray

Returns:

numpy.ndarray – Complex-valued \(\psi(x,t)\), shape (len(t), len(x)).

class physicskit.quantum.EigenResult(x, energies, wavefunctions)[source]#

Bases: object

Container for the output of NumerovSolver.solve().

Parameters:
energies: ndarray#

Eigenenergies \(E_n\), ascending, shape (n_states,).

Type:

numpy.ndarray

wavefunctions: ndarray#

Normalized eigenfunctions \(\psi_n(x)\), shape (n_states, len(x)), with \(\int \lvert\psi_n\rvert^2\,dx = 1\) and \(\psi_n(x[0]) = \psi_n(x[-1]) = 0\).

Type:

numpy.ndarray

x: ndarray#

The spatial grid the wavefunctions are sampled on.

Type:

numpy.ndarray

class physicskit.quantum.FiniteSquareWell(V0=20.0, width=2.0, hbar=1.0, m=1.0)[source]#

Bases: object

A finite square well of depth \(V_0\) and width \(a\),

\[\begin{split}V(x) = \begin{cases} -V_0 & \lvert x\rvert < a/2 \\ 0 & \text{otherwise} \end{cases}.\end{split}\]

Provides bound states (\(E<0\)) via the Numerov solver and exact analytic transmission/reflection coefficients for scattering states (\(E>0\)).

Parameters:

  • V0 (float) – Well depth (\(V_0 > 0\); the floor sits at \(-V_0\)).

  • width (float) – Full width \(a\) of the well.

  • hbar (float) – Value of \(\hbar\) to use.

  • m (float) – Particle mass.

bound_states(x_extent=6.0, n_points=1600, n_states=6)[source]#

Compute the bound states (\(E<0\)).

Parameters:

  • x_extent (float) – Half-width of the numerical domain.

  • n_points (int) – Number of grid points.

  • n_states (int) – Number of Numerov eigenstates to request (some may be discarded if they lie above the well rim; see Notes).

Return type:

EigenResult

Returns:

physicskit.quantum.core.eigensolvers.EigenResult – Only the states with \(E<0\) are kept – the solver also returns box-quantized “continuum” states above the well rim, which are artifacts of the finite numerical domain and are filtered out here.

evanescent_wavefunction(result, n=0)[source]#

The \(n\)-th bound-state wavefunction.

Its exponential tails outside \(\lvert x\rvert < a/2\) are evanescent by construction of the Numerov solver.

Parameters:
Return type:

ndarray

Returns:

numpy.ndarray

potential(x)[source]#

Evaluate \(V(x)\) (a barrier if V0 < 0).

Parameters:

x (ndarray)

Return type:

ndarray

Returns:

numpy.ndarray

scattering(E)[source]#

Exact 1D transmission/reflection coefficients for \(E>0\).

\[T(E) = \left[1 + \frac{V_0^2 \sin^2(k_2 a)}{4 E (E+V_0)}\right]^{-1}, \qquad R = 1 - T,\]

with \(k_1=\sqrt{2mE}/\hbar\) outside the well and \(k_2=\sqrt{2m(E+V_0)}/\hbar\) inside it. \(T(E)=1\) exactly whenever \(k_2 a\) is a multiple of \(\pi\) (the Ramsauer-Townsend resonances).

The same formula analytically continues correctly into a classically forbidden interior (\(E+V_0<0\), e.g. a barrier built by passing a negative V0): \(k_2\) becomes purely imaginary, \(\sin\) becomes \(i\sinh\) of a real argument, and \(T\) reduces to the standard sub-barrier tunneling formula. This method evaluates \(k_2\) with complex arithmetic so that regime is handled automatically.

Parameters:

E (float) – Incident kinetic energy (\(E>0\)).

Return type:

ScatteringResult

Returns:

ScatteringResult

transmission_spectrum(E_values)[source]#

Transmission probability \(T(E)\) over an array of energies.

Parameters:

E_values (ndarray) – Incident energies (\(E>0\)).

Return type:

ndarray

Returns:

numpy.ndarray – \(T(E)\) for each value in E_values.

wavepacket_scattering(x_extent=25.0, n_points=2048, x0=-10.0, sigma0=1.5, k0=None, dt=0.0002, n_steps=6000, save_every=30)[source]#

Propagate a Gaussian wavepacket incident on this well/barrier.

A single FFT split-operator run (SplitOperatorSolver1D) shared by the tunnelling (V0 < 0, a barrier) and the resonant-scattering (V0 > 0, a well) demos: the packet splits into a reflected and a transmitted piece as it crosses \(\lvert x\rvert < \text{width}/2\).

Parameters:

  • x_extent (float) – Half-width of the numerical domain.

  • n_points (int) – Number of grid points.

  • x0 (float) – Initial center of the packet (upstream of the well/barrier).

  • sigma0 (float) – Initial packet width.

  • k0 (float | None) – Central wavenumber; defaults to a value giving mean kinetic energy comparable to \(\lvert V_0\rvert\), so both partial reflection and partial transmission occur.

  • dt (float) – Time step.

  • n_steps (int) – Total number of time steps.

  • save_every (int) – Save a snapshot every this many steps.

Returns:

  • x (numpy.ndarray) – The spatial grid.

  • frames (numpy.ndarray) – Complex wavefunction snapshots, shape (n_saved, len(x)).

  • times (numpy.ndarray) – The time of each snapshot.

class physicskit.quantum.FloquetDrivenBox(L=1.0, V0=5.0, omega=20.0, hbar=1.0, m=1.0, wall_height=5000.0, pad=0.4, n_grid=2048)[source]#

Bases: object

An infinite square well \([0,L]\) driven by an AC dipole field.

\[V(x,t) = V_0 (x - L/2) \cos(\omega t),\]

simulated with the FFT split-operator propagator. The box walls are emulated with a steep confining potential outside \([0,L]\) padded into a larger FFT domain.

Provides both time-domain propagation (to see multiphoton Rabi oscillations / photon-assisted tunneling between box eigenstates) and the Floquet quasi-energy spectrum from the one-period propagator \(\hat U(T)\).

Parameters:

  • L (float) – Well width.

  • V0 (float) – Driving field amplitude.

  • omega (float) – Driving angular frequency \(\omega\).

  • hbar (float) – Value of \(\hbar\) to use.

  • m (float) – Particle mass.

  • wall_height (float) – Height of the confining walls outside \([0,L]\).

  • pad (float) – Fractional padding of the FFT domain beyond \([0,L]\).

  • n_grid (int) – Number of FFT grid points.

box_eigenstate(n)[source]#

The \(n\)-th unperturbed infinite-well eigenstate on the grid.

Parameters:

n (int) – Quantum number (\(n=1,2,\dots\)).

Return type:

ndarray

Returns:

numpy.ndarray

box_energy(n)[source]#

Unperturbed box eigenenergy \(E_n = n^2\pi^2\hbar^2/(2mL^2)\).

Parameters:

n (int)

Return type:

float

Returns:

float

floquet_quasienergies(n_levels=4, dt=0.0001)[source]#

Floquet quasi-energy spectrum from the one-period propagator.

Diagonalizes \(\hat U(T)\) (\(T=2\pi/\omega\)), restricted to the span of the lowest n_levels box eigenstates, to obtain quasi-energies

\[\epsilon_k = -\frac{\hbar}{T}\arg(\lambda_k),\]

where \(\lambda_k\) are the eigenvalues of \(\hat U(T)\).

Parameters:

  • n_levels (int) – Size of the truncated box-eigenstate basis.

  • dt (float) – Time step used to build \(\hat U(T)\).

Return type:

ndarray

Returns:

numpy.ndarray – Sorted quasi-energies, shape (n_levels,).

make_solver(dt=0.0001)[source]#

Build a SplitOperatorSolver1D for this drive.

Parameters:

dt (float) – Time step.

Return type:

SplitOperatorSolver1D

Returns:

physicskit.quantum.core.solvers.SplitOperatorSolver1D

transition_probability(n_initial, n_final, t_max, dt=0.0001)[source]#

Multiphoton transition probability \(P_{i\to f}(t)\).

Propagates the n_initial box eigenstate under the AC drive and returns its overlap probability with the n_final unperturbed box eigenstate – shows resonant multiphoton transitions when \(\hbar\omega\) (or a multiple) bridges the level gap.

Parameters:

  • n_initial (int) – Initial box quantum number.

  • n_final (int) – Final box quantum number to project onto.

  • t_max (float) – Total propagation time.

  • dt (float) – Time step.

Return type:

tuple[ndarray, ndarray]

Returns:

  • times (numpy.ndarray)

  • probs (numpy.ndarray) – \(P_{n_\text{initial}\to n_\text{final}}(t)\).

class physicskit.quantum.GaussianDispersion(x0=0.0, sigma0=1.0, k0=5.0, hbar=1.0, m=1.0)[source]#

Bases: object

Free-particle spreading of a Gaussian wavepacket.

The exact analytic solution of the free-particle Schrodinger equation for a Gaussian initial condition; used to benchmark the FFT split-operator propagator.

Parameters:

  • x0 (float) – Initial center of the packet.

  • sigma0 (float) – Initial width.

  • k0 (float) – Central wavenumber.

  • hbar (float) – Value of \(\hbar\) to use.

  • m (float) – Particle mass.

center(t)[source]#

The ballistic center \(x_0 + (\hbar k_0/m)\,t\).

Parameters:

t (ndarray)

Return type:

ndarray

Returns:

numpy.ndarray

density(x, t)[source]#

The probability density \(\lvert\psi(x,t)\rvert^2\).

Parameters:

  • x (ndarray) – Positions to evaluate at.

  • t (float) – Time.

Return type:

ndarray

Returns:

numpy.ndarray

hbar: float = 1.0#
k0: float = 5.0#
m: float = 1.0#
psi(x, t)[source]#

The exact free-particle wavefunction \(\psi(x,t)\).

Parameters:

  • x (ndarray) – Positions to evaluate at.

  • t (float) – Time.

Return type:

ndarray

Returns:

numpy.ndarray – Complex-valued.

sigma0: float = 1.0#
trajectory(x, t_values)[source]#

Stack of \(\psi(x,t)\) snapshots over a time grid.

Shaped for animate_density().

Parameters:

  • x (ndarray) – Positions to evaluate at.

  • t_values (ndarray) – Times to evaluate at.

Return type:

ndarray

Returns:

numpy.ndarray – Complex-valued, shape (len(t_values), len(x)).

width(t)[source]#

The spreading width \(\sigma(t)\).

\[\sigma(t) = \sigma_0\sqrt{1 + \left(\frac{\hbar t}{2m\sigma_0^2}\right)^2}.\]
Parameters:

t (ndarray)

Return type:

ndarray

Returns:

numpy.ndarray

x0: float = 0.0#
class physicskit.quantum.HarmonicOscillator(m=1.0, omega=1.0, hbar=1.0)[source]#

Bases: object

The 1D quantum harmonic oscillator,

\[\hat H = \frac{\hat p^2}{2m} + \frac{1}{2} m \omega^2 \hat x^2.\]
Parameters:

  • m (float) – Particle mass.

  • omega (float) – Angular frequency \(\omega\).

  • hbar (float) – Value of \(\hbar\) to use.

coherent_n_max(alpha, margin=8.0)[source]#

Fock-basis truncation order for a coherent state \(\lvert\alpha\rangle\).

Captures the Poisson-distributed photon number \(n \sim \lvert\alpha\rvert^2\) out to several standard deviations \(\sqrt{\lvert\alpha\rvert^2}\).

Parameters:

  • alpha (complex) – Coherent-state amplitude.

  • margin (float) – Number of standard deviations of safety margin.

Return type:

int

Returns:

int

coherent_trajectory(alpha, t)[source]#

Classical mean position/momentum of a coherent state.

\[\langle x\rangle(t) = x_0\sqrt2\,\mathrm{Re}\!\left[\alpha e^{-i\omega t}\right], \qquad \langle p\rangle(t) = \sqrt{2\hbar m \omega}\,\mathrm{Im}\!\left[\alpha e^{-i\omega t}\right].\]
Parameters:

  • alpha (complex) – Coherent-state amplitude.

  • t (ndarray) – Times to evaluate at.

Returns:

x_t, p_t (numpy.ndarray)

coherent_wavefunction(alpha, x, t=0.0, n_max=None)[source]#

The Glauber coherent state \(\lvert\alpha\rangle\) in position space.

\[\psi_\alpha(x,t) = \sum_n c_n\, \phi_n(x)\, e^{-i E_n t/\hbar}, \qquad c_n = e^{-\lvert\alpha\rvert^2/2}\frac{\alpha^n}{\sqrt{n!}},\]

the Poissonian Fock coefficients of \(\lvert\alpha\rangle\). The packet is a Gaussian of fixed shape that oscillates rigidly, matching the classical trajectory (see coherent_trajectory()).

Parameters:
Return type:

ndarray

Returns:

numpy.ndarray – Complex-valued \(\psi_\alpha(x,t)\).

eigenfunction(n, x)[source]#

The Fock (energy) eigenstate \(\phi_n(x)\).

Normalized so that \(\int \lvert\phi_n\rvert^2\,dx = 1\).

Parameters:

  • n (int) – Quantum number.

  • x (ndarray) – Positions to evaluate at.

Return type:

ndarray

Returns:

numpy.ndarray

eigenfunction_t(n, x, t)[source]#

The time-evolved stationary state \(\psi_n(x,t) = \phi_n(x)\, e^{-i E_n t/\hbar}\).

A pure phase-rotating stationary state: \(\lvert\psi_n(x,t)\rvert^2\) is time-independent.

Parameters:

  • n (int) – Quantum number.

  • x (ndarray) – Positions to evaluate at.

  • t (float) – Time.

Return type:

ndarray

Returns:

numpy.ndarray – Complex-valued \(\psi_n(x,t)\).

eigenfunctions(n_max, x)[source]#

Stack of eigenstates \(\phi_0, \dots, \phi_{n_\text{max}}(x)\).

Parameters:

  • n_max (int) – Highest quantum number to include.

  • x (ndarray) – Positions to evaluate at.

Return type:

ndarray

Returns:

numpy.ndarray – Shape (n_max+1, len(x)).

energy(n)[source]#

Eigenenergy \(E_n = \hbar\omega(n+1/2)\).

Parameters:

n (int or array_like) – Quantum number(s).

Return type:

ndarray

Returns:

numpy.ndarray

hbar: float = 1.0#
m: float = 1.0#
omega: float = 1.0#
squeezed_uncertainties(r)[source]#

Position/momentum uncertainties of the \(r\)-squeezed vacuum (\(\phi=0\) axis).

Parameters:

r (float) – Squeezing parameter.

Returns:

dx, dp (float) – \(\Delta x\), \(\Delta p\); their product stays at the minimum \(\hbar/2\) for all \(r\).

squeezed_vacuum_wavefunction(x, r, phi=0.0)[source]#

Position-space squeezed vacuum \(\lvert r,\phi\rangle = \hat S(re^{2i\phi})\lvert0\rangle\).

With \(\hat S(\xi)=\exp[(\xi^*\hat a^2-\xi\hat a^{\dagger2})/2]\) and \(z=e^{2i\phi}\tanh r\), the exact wavefunction is the complex Gaussian

\[\psi(x) \propto \exp\!\left[-\frac{x^2}{2x_0^2}\,\frac{1+z}{1-z}\right], \qquad x_0=\sqrt{\hbar/m\omega},\]

squeezed along the rotated quadrature \(x\cos\phi+(p/m\omega)\sin\phi\) (variance \(e^{-2r}\) times the zero-point value) and anti-squeezed along the orthogonal one. In particular

\[\Delta x^2 = \frac{\hbar}{2 m \omega}\left(e^{-2r}\cos^2\phi + e^{2r}\sin^2\phi\right),\]

which is \(\frac{\hbar}{2m\omega}e^{-2r}\) at \(\phi=0\) (with \(\Delta p^2 = (\hbar m \omega/2)\, e^{2r}\)); at \(r=0\) every \(\phi\) gives the vacuum.

Parameters:

  • x (ndarray) – Positions to evaluate at.

  • r (float) – Squeezing parameter.

  • phi (float) – Squeezing angle: rotates the uncertainty ellipse in phase space by \(\phi\) (\(\phi=\pi/2\) squeezes momentum instead).

Return type:

ndarray

Returns:

numpy.ndarray – Complex-valued (real when phi=0) squeezed-vacuum wavefunction.

superposition_trajectory(n_values, coeffs, x, t_values)[source]#

Stack of superposition_wavefunction() snapshots over a time grid.

Shaped for animate_density().

Parameters:

  • n_values (array_like of int) – Quantum numbers included in the superposition.

  • coeffs (array_like of complex) – Amplitudes \(c_k\).

  • x (ndarray) – Positions to evaluate at.

  • t_values (ndarray) – Times to evaluate at.

Return type:

ndarray

Returns:

numpy.ndarray – Complex-valued, shape (len(t_values), len(x)).

superposition_wavefunction(n_values, coeffs, x, t)[source]#

A general time-evolved superposition of Fock eigenstates.

\[\psi(x,t) = \sum_k c_k\, \phi_{n_k}(x)\, e^{-i E_{n_k} t/\hbar},\]

with beats at every pairwise frequency \((E_{n_k}-E_{n_j})/\hbar\) – unlike a single stationary state (eigenfunction_t()) or a coherent state (coherent_wavefunction(), a rigidly oscillating Gaussian), a generic superposition’s density genuinely reshapes over time.

Parameters:

  • n_values (array_like of int) – Quantum numbers included in the superposition.

  • coeffs (array_like of complex) – Amplitudes \(c_k\) (renormalized internally so \(\sum_k \lvert c_k\rvert^2=1\)).

  • x (ndarray) – Positions to evaluate at.

  • t (float) – Time.

Return type:

ndarray

Returns:

numpy.ndarray – Complex-valued \(\psi(x,t)\).

thermal_populations(T, n_max, k_B=1.0)[source]#

Diagonal Boltzmann populations of the thermal density matrix.

\[\rho_{nn} = (1-e^{-x})\, e^{-nx}, \qquad x = \frac{\hbar\omega}{k_B T},\]

for \(\hat\rho = e^{-\beta\hat H}/Z\) (T=0 returns the pure ground state).

Parameters:

  • T (float) – Temperature (T<=0 returns the ground state).

  • n_max (int) – Highest Fock level to include.

  • k_B (float) – Boltzmann constant.

Return type:

ndarray

Returns:

numpy.ndarray – Populations \(\rho_{nn}\), shape (n_max+1,), summing to 1.

thermal_position_distribution(x, T, k_B=1.0, n_max=200)[source]#

The (temperature-broadened) position-space probability density.

\[P(x) = \sum_n \rho_{nn}\, \lvert\phi_n(x)\rvert^2.\]
Parameters:

  • x (ndarray) – Positions to evaluate at.

  • T (float) – Temperature.

  • k_B (float) – Boltzmann constant.

  • n_max (int) – Highest Fock level to include.

Return type:

ndarray

Returns:

numpy.ndarray

thermal_position_variance_analytic(T, k_B=1.0)[source]#

Closed-form position variance of the thermal state.

\[\langle x^2\rangle = \frac{\hbar}{2m\omega} \coth\!\left(\frac{\hbar\omega}{2 k_B T}\right),\]

reducing to the zero-point value \(\hbar/2m\omega\) as \(T\to 0\).

Parameters:

  • T (float) – Temperature.

  • k_B (float) – Boltzmann constant.

Return type:

float

Returns:

float

property x0: float#

Natural length scale \(\sqrt{\hbar/(m\omega)}\).

Type:

float

zero_point_uncertainties()[source]#

Position/momentum uncertainties of the ordinary ground state.

Equivalent to squeezed_uncertainties(r=0), provided for comparison.

Returns:

dx, dp (float)

class physicskit.quantum.HydrogenOrbital(n, l, m, Z=1, a0=1.0)[source]#

Bases: object

The hydrogenic wavefunction \(\psi_{nlm}(r,\theta,\phi) = R_{nl}(r)\, Y_l^m(\theta,\phi)\).

Parameters:

  • n (int) – Principal quantum number.

  • l (int) – Orbital angular momentum quantum number (\(0 \le l < n\)).

  • m (int) – Magnetic quantum number (\(\lvert m\rvert \le l\)).

  • Z (int) – Nuclear charge.

  • a0 (float) – Bohr radius.

Raises:

ValueError – If the quantum numbers violate \(0 \le l < n\) or \(\lvert m\rvert \le l\).

Z: int = 1#
a0: float = 1.0#
angular(theta, phi)[source]#

The angular part \(Y_l^m(\theta,\phi)\); see spherical_harmonic().

Parameters:
Return type:

ndarray

Returns:

numpy.ndarray

check_normalization(r_max=None, n_points=20000)[source]#

Integrate radial_density() over \([0, r_\text{max}]\) (should be 1).

Parameters:

  • r_max (float | None) – Upper integration limit; defaults to a value scaled to the orbital’s expected extent.

  • n_points (int) – Number of quadrature points.

Return type:

float

Returns:

float

density(r, theta, phi)[source]#

The probability density \(\lvert\psi_{nlm}\rvert^2\).

Parameters:
Return type:

ndarray

Returns:

numpy.ndarray

property energy: float#

Energy \(E_n = -Z^2/(2n^2)\) Hartree (atomic units, \(\hbar = m_e = e = 1\)).

Type:

float

l: int#
m: int#
most_probable_radius(r_max=None, n_points=20000)[source]#

The radius at which radial_density() is maximal.

Parameters:

  • r_max (float | None) – Upper bound of the search range; defaults to a value scaled to the orbital’s expected extent.

  • n_points (int) – Number of grid points used for the search.

Return type:

float

Returns:

float

n: int#
psi(r, theta, phi)[source]#

The full wavefunction \(\psi_{nlm}(r,\theta,\phi)\).

Parameters:
Return type:

ndarray

Returns:

numpy.ndarray – Complex-valued.

radial(r)[source]#

The radial part \(R_{nl}(r)\); see radial_wavefunction().

Parameters:

r (ndarray)

Return type:

ndarray

Returns:

numpy.ndarray

radial_density(r)[source]#

The radial probability density \(P(r) = r^2 \lvert R_{nl}(r)\rvert^2\).

The probability of finding the electron in a shell \([r, r+dr]\), already integrated over solid angle.

Parameters:

r (ndarray)

Return type:

ndarray

Returns:

numpy.ndarray

class physicskit.quantum.IsingEntangler(J=1.0, hbar=1.0)[source]#

Bases: object

Two initially unentangled qubits, entangled by an Ising coupling.

Both qubits start in \(\lvert+\rangle = (\lvert0\rangle+\lvert1\rangle)/\sqrt2\) (an eigenstate of \(\sigma_x\), i.e. unentangled and unbiased in the measurement basis), and evolve under

\[\hat H = \hbar J\, \sigma_z^{(1)}\otimes\sigma_z^{(2)}.\]

Since \(\lvert+\rangle\otimes\lvert+\rangle\) is a superposition of all four \(\sigma_z\otimes\sigma_z\) eigenstates, \(e^{-i\hat Ht/\hbar}\) accumulates a different phase on each computational-basis term, generating genuine entanglement – maximal (a Bell-equivalent state) at \(Jt=\pi/4\).

Parameters:

  • J (float) – Ising coupling strength.

  • hbar (float) – Value of \(\hbar\) to use.

J: float = 1.0#
concurrence(psi)[source]#

The Wootters concurrence, a two-qubit entanglement measure.

For a pure state, \(C = 2\lvert\det\Psi\rvert\) where \(\Psi\) is the state reshaped into a \(2\times2\) matrix (equivalently \(C=\sqrt{2(1-\mathrm{Tr}\,\rho_1^2)}\)); \(C=0\) for a product state and \(C=1\) for a maximally entangled (Bell) state.

Parameters:

psi (ndarray) – The 4-component two-qubit state vector.

Return type:

float

Returns:

float

concurrence_trajectory(t_values)[source]#

The concurrence \(C(t)\) over a time grid.

Parameters:

t_values (ndarray) – Times to evaluate at.

Return type:

ndarray

Returns:

numpy.ndarray

hbar: float = 1.0#
initial_state()[source]#

The unentangled product state \(\lvert+\rangle\otimes\lvert+\rangle\).

Return type:

ndarray

Returns:

numpy.ndarray – The 4-component state vector.

purity(psi)[source]#

The purity \(\mathrm{Tr}(\rho_1^2)\) of the reduced state.

Equal to 1 for an unentangled (product) state and 1/2 (minimal, for a single qubit) at maximal entanglement.

Parameters:

psi (ndarray) – The 4-component two-qubit state vector.

Return type:

float

Returns:

float

reduced_density_matrix(psi)[source]#

The reduced density matrix \(\rho_1=\mathrm{Tr}_2\lvert\psi\rangle\langle\psi\rvert\) of the first qubit.

Parameters:

psi (ndarray) – The 4-component two-qubit state vector.

Return type:

ndarray

Returns:

numpy.ndarray – The \(2\times2\) reduced density matrix.

state(t)[source]#

The evolved two-qubit state \(e^{-i\hat Ht/\hbar}\lvert\psi_0\rangle\).

\(\hat H\) is diagonal in the computational basis with eigenvalues \(\hbar J\) on \(\lvert00\rangle,\lvert11\rangle\) and \(-\hbar J\) on \(\lvert01\rangle,\lvert10\rangle\), so the propagator is applied as a plain phase per basis state.

Parameters:

t (float) – Time.

Return type:

ndarray

Returns:

numpy.ndarray – The 4-component state vector, in the \(\lvert00\rangle,\lvert01\rangle,\lvert10\rangle,\lvert11\rangle\) basis.

state_trajectory(t_values)[source]#

Stack of state() snapshots over a time grid.

Parameters:

t_values (ndarray) – Times to evaluate at.

Return type:

ndarray

Returns:

numpy.ndarray – Complex-valued, shape (len(t_values), 4).

class physicskit.quantum.NumerovSolver(x, V, hbar=1.0, m=1.0)[source]#

Bases: object

Bound-state solver for an arbitrary 1D potential \(V(x)\) on a fixed grid.

Parameters:

  • x (ndarray) – Uniform spatial grid. The wavefunction is forced to vanish at x[0] and x[-1] (hard walls), so the domain should extend well into the classically forbidden region for soft potentials.

  • V (callable or array_like) – Potential energy, either a function V(x) or a precomputed array matching x.

  • hbar (float) – Value of \(\hbar\) to use.

  • m (float) – Particle mass.

Notes

The Numerov finite-difference relation for \(\psi'' = f\psi\) is recast as the generalized eigenvalue problem

\[\hat A \boldsymbol\psi = E \hat B \boldsymbol\psi, \qquad \hat A = -\frac{\hbar^2}{2m h^2}\hat T + \hat B\,\mathrm{diag}(V),\]

where \(\hat T\) and \(\hat B\) are the tridiagonal matrices \(\mathrm{tridiag}(1,-2,1)\) and \(\mathrm{tridiag}(1,10,1)/12\) and \(h\) is the grid spacing. Validated against the analytic infinite-well and harmonic-oscillator spectra to four decimal places.

solve(n_states=6)[source]#

Compute the lowest bound states.

Parameters:

n_states (int) – Number of lowest-energy eigenstates to return. Must be smaller than the number of interior grid points.

Return type:

EigenResult

Returns:

EigenResult – The energies and normalized wavefunctions.

Raises:

ValueError – If n_states is not smaller than the number of interior grid points.

class physicskit.quantum.QuantumRevival(L=1.0, hbar=1.0, m=1.0, n_max=400)[source]#

Bases: object

A localized wavepacket in an infinite square well of width \(L\).

Because \(E_n = n^2\pi^2\hbar^2/(2mL^2)\) is quadratic in \(n\), the packet disperses into apparent chaos yet exactly reassembles into replicas of \(\psi(x,0)\) at fractional multiples of the revival time

\[t_\text{rev} = \frac{4mL^2}{\pi\hbar}\]

(since \(E_n = n^2\cdot(\hbar\cdot 2\pi/t_\text{rev})/2\)), with a mirror-image replica at \(t_\text{rev}/2\), and reduced-scale “fractional revivals” (clones/anti-clones) at rational fractions \(t_\text{rev}\cdot(p/q)\).

Parameters:

  • L (float) – Well width.

  • hbar (float) – Value of \(\hbar\) to use.

  • m (float) – Particle mass.

  • n_max (int) – Number of box eigenstates used in the expansion.

L: float = 1.0#
eigenbasis_coefficients(psi0_func, n_grid=4000)[source]#

Box-eigenbasis expansion coefficients of an initial state.

\[c_n = \langle\phi_n\rvert\psi_0\rangle, \qquad \phi_n(x) = \sqrt{2/L}\,\sin(n\pi x/L).\]
Parameters:

  • psi0_func (callable) – Function psi0_func(x) -> array giving the initial state.

  • n_grid (int) – Number of quadrature points used to evaluate the overlap integral.

Return type:

ndarray

Returns:

numpy.ndarray – Coefficients \(c_n\), shape (n_max,).

energy(n)[source]#

Box eigenenergy \(E_n = n^2\pi^2\hbar^2/(2mL^2)\).

Parameters:

n (ndarray)

Return type:

ndarray

Returns:

numpy.ndarray

fidelity_to_initial(x, t, coeffs, psi0)[source]#

Fidelity \(\lvert\langle\psi_0\rvert\psi(t)\rangle\rvert^2\) of the evolved state to the initial state.

Approaches 1 at (fractional) revival times.

Parameters:
Return type:

float

Returns:

float

gaussian_initial_state(x0, sigma)[source]#

Build a normalized Gaussian initial-state function for the box.

Parameters:

  • x0 (float) – Center of the packet (should lie well inside (0, L)).

  • sigma (float) – Width of the packet.

Returns:

callable – A function psi0(x) -> array (zero outside [0, L]).

hbar: float = 1.0#
m: float = 1.0#
n_max: int = 400#
property revival_time: float#

The revival time \(t_\text{rev} = 4mL^2/(\pi\hbar)\).

Type:

float

wavefunction(x, t, coeffs)[source]#

Evaluate the evolved state \(\psi(x,t) = \sum_n c_n \phi_n(x) e^{-iE_n t/\hbar}\).

Parameters:
Return type:

ndarray

Returns:

numpy.ndarray – Complex-valued \(\psi(x,t)\).

class physicskit.quantum.RabiProblem(omega0=1.0, omega_d=1.0, Omega=0.2, hbar=1.0)[source]#

Bases: object

A driven two-level system in the rotating-wave approximation (RWA).

In the frame rotating at the drive frequency \(\omega_d\), the RWA Hamiltonian is time-independent,

\[\hat H_\text{RWA} = \frac{\hbar}{2}\left(\Delta\,\sigma_z + \Omega\,\sigma_x\right), \qquad \Delta = \omega_d - \omega_0,\]

giving the exact closed-form propagator (the Rabi formula)

\[U(t) = \cos\!\left(\frac{\Omega_R t}{2}\right) I - i\sin\!\left(\frac{\Omega_R t}{2}\right) \frac{\Delta\,\sigma_z + \Omega\,\sigma_x}{\Omega_R}, \qquad \Omega_R = \sqrt{\Delta^2+\Omega^2},\]

with \(\Omega_R\) the generalized Rabi frequency. On resonance (\(\Delta=0\)) the population fully cycles between the two levels; detuning caps the cycling amplitude at \(\Omega^2/\Omega_R^2\).

Parameters:

  • omega0 (float) – Bare qubit splitting \(\omega_0\).

  • omega_d (float) – Drive frequency \(\omega_d\).

  • Omega (float) – Drive (Rabi) coupling strength \(\Omega\).

  • hbar (float) – Value of \(\hbar\) to use.

Omega: float = 0.2#
property detuning: float#

The detuning \(\Delta = \omega_d - \omega_0\).

Type:

float

excited_state_population(t)[source]#

The Rabi formula \(P_e(t) = (\Omega^2/\Omega_R^2)\sin^2(\Omega_R t/2)\).

Parameters:

t (ndarray) – Times to evaluate at.

Return type:

ndarray

Returns:

numpy.ndarray

property generalized_rabi_frequency: float#

\(\Omega_R = \sqrt{\Delta^2+\Omega^2}\).

Type:

float

hbar: float = 1.0#
omega0: float = 1.0#
omega_d: float = 1.0#
propagator(t)[source]#

The RWA time-evolution operator \(U(t)\).

Parameters:

t (float) – Time.

Return type:

ndarray

Returns:

numpy.ndarray – The \(2\times2\) unitary propagator.

state(t, psi0=None)[source]#

The evolved qubit state \(\psi(t) = U(t)\,\psi_0\).

Parameters:

  • t (float) – Time.

  • psi0 (ndarray | None) – Initial state; defaults to the ground state \(\lvert0\rangle\).

Return type:

ndarray

Returns:

numpy.ndarray – The 2-component qubit state.

state_trajectory(t_values, psi0=None)[source]#

Stack of state() snapshots over a time grid.

Shaped for state_to_bloch_trajectory().

Parameters:

  • t_values (ndarray) – Times to evaluate at.

  • psi0 (ndarray | None) – Initial state; defaults to the ground state.

Return type:

ndarray

Returns:

numpy.ndarray – Complex-valued, shape (len(t_values), 2).

class physicskit.quantum.RectangularBox2D(Lx=1.0, Ly=1.0, hbar=1.0, m=1.0)[source]#

Bases: object

Infinite 2D box of size \(L_x \times L_y\).

\[E_{n_x,n_y} = \frac{\pi^2\hbar^2}{2m} \left(\frac{n_x^2}{L_x^2} + \frac{n_y^2}{L_y^2}\right).\]

Rational ratios \(L_x/L_y\) (e.g. the square box \(1{:}1\)) produce accidental degeneracies between distinct \((n_x,n_y)\) pairs.

Parameters:

  • Lx (float) – Box dimensions.

  • Ly (float) – Box dimensions.

  • hbar (float) – Value of \(\hbar\) to use.

  • m (float) – Particle mass.

degeneracies(n_max=6, tol=1e-09)[source]#

Group states by (rounded) energy to reveal degenerate levels.

Parameters:

  • n_max (int) – Maximum quantum number along each axis.

  • tol (float) – Energy values within tol are treated as degenerate.

Return type:

dict

Returns:

dict – Maps each degenerate energy to its list of (nx, ny) pairs; non-degenerate levels are omitted.

energy(nx, ny)[source]#

Eigenenergy \(E_{n_x,n_y}\).

Parameters:

  • nx (int) – Quantum numbers (each \(\ge 1\)).

  • ny (int) – Quantum numbers (each \(\ge 1\)).

Return type:

float

Returns:

float

grid(n_points=200)[source]#

A Cartesian meshgrid spanning the box.

Parameters:

n_points (int) – Number of points along each axis.

Returns:

X, Y (numpy.ndarray) – Meshgrid arrays (indexing='ij').

spectrum(n_max=6)[source]#

All eigenstates with \(1 \le n_x, n_y \le n_\text{max}\), sorted by energy.

Parameters:

n_max (int) – Maximum quantum number along each axis.

Return type:

list

Returns:

list of Box2DEigenstate

class physicskit.quantum.ScatteringResult(E, k1, k2, T, R)[source]#

Bases: object

Container for the output of FiniteSquareWell.scattering().

Parameters:
E: float#

Incident energy.

Type:

float

R: float#

Reflection probability (\(R = 1 - T\)).

Type:

float

T: float#

Transmission probability.

Type:

float

k1: float#

Wavenumber outside the well.

Type:

float

k2: complex#

Wavenumber inside the well/barrier. Purely imaginary in a classically forbidden interior (sub-barrier tunneling), real otherwise.

Type:

complex

class physicskit.quantum.SplitOperatorSolver1D(x, V, hbar=1.0, m=1.0, dt=0.001)[source]#

Bases: object

Propagate a 1D wavepacket under potential \(V(x)\) or \(V(x,t)\).

Parameters:

  • x (ndarray) – Uniform spatial grid.

  • V (callable or array_like) – Potential energy: a function V(x) (static) or V(x, t) (time-dependent, detected by inspecting the function’s signature), or a precomputed array matching x.

  • hbar (float) – Value of \(\hbar\) to use.

  • m (float) – Particle mass.

  • dt (float) – Time step.

norm(psi)[source]#

Total probability \(\int \lvert\psi(x)\rvert^2\,dx\).

Parameters:

psi (ndarray) – Wavefunction sampled on self.x.

Return type:

float

Returns:

float – The normalization integral (should stay at 1.0 under unitary evolution).

propagate(psi0, n_steps, t0=0.0, save_every=1)[source]#

Propagate an initial state and record its history.

Parameters:

  • psi0 (ndarray) – Initial wavefunction (renormalized internally).

  • n_steps (int) – Total number of time steps to take.

  • t0 (float) – Initial time.

  • save_every (int) – Save a snapshot every save_every steps (plus the initial state).

Return type:

ndarray | tuple[ndarray, ndarray]

Returns:

  • frames (numpy.ndarray) – Stacked wavefunction snapshots, shape (n_saved, len(x)).

  • times (numpy.ndarray) – The time of each snapshot, shape (n_saved,).

step(psi, t=0.0)[source]#

Advance the wavefunction by one time step dt.

Parameters:

  • psi (ndarray) – Complex wavefunction at time t, sampled on self.x.

  • t (float) – Current time (only used if the potential is time-dependent).

Return type:

ndarray

Returns:

numpy.ndarray – The wavefunction at time t + dt.

class physicskit.quantum.SplitOperatorSolver2D(x, y, V, hbar=1.0, m=1.0, dt=0.001)[source]#

Bases: object

Propagate a 2D wavepacket under potential \(V(x,y)\) or \(V(x,y,t)\).

Parameters:

  • x (ndarray) – Uniform spatial grids along each axis.

  • y (ndarray) – Uniform spatial grids along each axis.

  • V (callable or array_like) – Potential energy: a function V(x, y) (static) or V(x, y, t) (time-dependent, detected by inspecting the function’s signature), or a precomputed array matching the (x, y) meshgrid.

  • hbar (float) – Value of \(\hbar\) to use.

  • m (float) – Particle mass.

  • dt (float) – Time step.

norm(psi)[source]#

Total probability \(\iint \lvert\psi(x,y)\rvert^2\,dx\,dy\).

Parameters:

psi (ndarray) – Wavefunction sampled on the (x, y) grid.

Return type:

float

Returns:

float – The normalization integral (should stay at 1.0 under unitary evolution).

propagate(psi0, n_steps, t0=0.0, save_every=1)[source]#

Propagate an initial state and record its history.

Parameters:

  • psi0 (ndarray) – Initial wavefunction (renormalized internally).

  • n_steps (int) – Total number of time steps to take.

  • t0 (float) – Initial time.

  • save_every (int) – Save a snapshot every save_every steps (plus the initial state).

Returns:

  • frames (numpy.ndarray) – Stacked wavefunction snapshots, shape (n_saved, *X.shape).

  • times (numpy.ndarray) – The time of each snapshot, shape (n_saved,).

step(psi, t=0.0)[source]#

Advance the wavefunction by one time step dt.

Parameters:

  • psi (ndarray) – Complex wavefunction at time t, sampled on the (x, y) grid.

  • t (float) – Current time (only used if the potential is time-dependent).

Return type:

ndarray

Returns:

numpy.ndarray – The wavefunction at time t + dt.

class physicskit.quantum.StadiumBilliard2D(L=1.0, R=0.5, hbar=1.0, m=1.0)[source]#

Bases: object

Bunimovich stadium billiard: a rectangle of length L capped by two semicircles of radius R.

The classical dynamics is chaotic, and the quantum eigenstates of this box can exhibit scars – probability density concentrated along unstable classical periodic orbits – making this the natural bridge between physicskit.quantum and physicskit.chaos.

Eigenstates are obtained by direct diagonalization of the 2D finite-difference Laplacian on a grid, with the wavefunction forced to zero outside the stadium boundary (infinite-wall billiard).

Parameters:

  • L (float) – Length of the central rectangular section.

  • R (float) – Radius of the semicircular end-caps (also the half-height).

  • hbar (float) – Value of \(\hbar\) to use.

  • m (float) – Particle mass.

grid(n_points=120)[source]#

A Cartesian meshgrid enclosing the stadium.

Parameters:

n_points (int) – Number of points along the (longer) x-axis.

Returns:

  • (X, Y) (tuple of numpy.ndarray) – Meshgrid arrays (indexing='ij').

  • x, y (numpy.ndarray) – The 1D coordinate arrays used to build the meshgrid.

mask(X, Y)[source]#

Boolean mask, True inside the stadium.

The stadium is a central \(L \times 2R\) rectangle with two semicircular end-caps of radius \(R\).

Parameters:
Return type:

ndarray

Returns:

numpy.ndarray of bool

solve(n_points=160, n_states=6)[source]#

Diagonalize \(-\hbar^2/2m\,\nabla^2\) on the masked grid (Dirichlet boundary).

The curved end-caps are approximated by a “staircased” boundary (grid points are simply included or excluded by mask()), which converges only as \(O(h)\): low resolutions give eigenvalues with a few percent of resolution-dependent error and visibly blocky eigenstate edges. grid() keeps the grid isotropic (\(dx=dy\)) so the caps are at least staircased as circles rather than ellipses; beyond that, accuracy is improved simply by raising n_points – the sparse eigensolver used here stays fast (well under a second) even at n_points=260.

Parameters:

  • n_points (int) – Number of grid points along the x-axis (see grid()).

  • n_states (int) – Number of lowest-energy eigenstates to compute.

Returns:

  • energies (numpy.ndarray) – Eigenenergies, ascending, shape (n_states,).

  • wavefunctions (numpy.ndarray) – Normalized eigenfunctions, shape (n_states, *X.shape), zero outside the stadium.

  • X, Y (numpy.ndarray) – The coordinate meshgrid.

  • mask (numpy.ndarray of bool) – True inside the stadium (see mask()).

class physicskit.quantum.StarkLevel(n1, n2, m, shift)[source]#

Bases: object

One sublevel returned by stark_n2_quartet().

Parameters:
m: int#

Magnetic quantum number.

Type:

int

n1: int#

Parabolic quantum number.

Type:

int

n2: int#

Parabolic quantum number.

Type:

int

shift: float#

The first-order Stark energy shift.

Type:

float

class physicskit.quantum.SternGerlach(mu=1.0, grad_B=1.0, sigma0=1.0, sigma_x=2.0, k0=5.0, hbar=1.0, m=1.0, spin_weights=(0.5, 0.5))[source]#

Bases: object

Semiclassical Stern-Gerlach beam splitting.

A spin-1/2 particle beam enters an inhomogeneous field \(B_z(y) \approx B_0 + y\,\partial_y B_z\); the spin-dependent force \(F_y = \pm\mu\,\partial_y B_z\) (the sign set by the \(z\)-spin eigenvalue) pushes the two spin branches apart transversely while they drift along \(x\).

This implements the standard textbook simplified treatment: the two spin branches are modeled as two decoupled Gaussian wavepackets (rather than solving the fully entangled spin-position Hamiltonian), each in a uniform transverse force – exactly solvable, since a wavepacket under a constant force keeps the free-particle spreading law (GaussianDispersion) while its center follows the classical trajectory \(y(t) = y_0 + v_{y0}t + \tfrac12 a t^2\) (Kennard’s theorem / the Ehrenfest theorem being exact for a linear potential), picking up an extra “accelerated frame” phase \(\varphi(y,t) = m a t\, y/\hbar - m a^2 t^3/(6\hbar)\).

Parameters:

  • mu (float) – Magnetic-moment scale coupling the spin to the field gradient.

  • grad_B (float) – Field gradient \(\partial_y B_z\).

  • sigma0 (float) – Initial transverse (\(y\)) width of the beam.

  • sigma_x (float) – Initial longitudinal (\(x\), beam-propagation) width.

  • k0 (float) – Central longitudinal wavenumber (beam velocity \(\hbar k_0/m\)).

  • hbar (float) – Value of \(\hbar\) to use.

  • m (float) – Particle mass.

  • spin_weights (tuple) – Probabilities \((P_{+z}, P_{-z})\) of the incoming spin state in the measurement (\(z\)) basis – (0.5, 0.5) is the standard demo of an unpolarized or \(x\)-polarized incoming beam, split 50/50 by the apparatus.

acceleration(spin_sign)[source]#

Transverse acceleration \(a = \text{spin\_sign}\cdot\mu\,\partial_y B_z/m\).

Parameters:

spin_sign (float) – \(+1\) for spin-up, \(-1\) for spin-down.

Return type:

float

Returns:

float

grad_B: float = 1.0#
hbar: float = 1.0#
joint_density(x, y, t)[source]#

The joint spatial density on the \((x,y)\) plane at time \(t\).

\[P(x,y,t) = P_{+z}\,\lvert\psi(x,t)\rvert^2 \lvert\psi_+(y,t)\rvert^2 + P_{-z}\,\lvert\psi(x,t)\rvert^2 \lvert\psi_-(y,t)\rvert^2,\]

the incoherent (classically correlated, not entangled) mixture of the two spin branches – splitting into two lobes as \(t\) grows.

Parameters:

  • x (ndarray) – 1D position grids; combined into a meshgrid internally.

  • y (ndarray) – 1D position grids; combined into a meshgrid internally.

  • t (float) – Time.

Return type:

ndarray

Returns:

numpy.ndarray – Real-valued, shape (len(x), len(y)).

joint_density_stack(x, y, t_values)[source]#

Stack of joint_density() snapshots over a time grid.

Shaped for animate_density_2d().

Parameters:

  • x (ndarray) – 1D position grids.

  • y (ndarray) – 1D position grids.

  • t_values (ndarray) – Times to evaluate at.

Return type:

ndarray

Returns:

numpy.ndarray – Real-valued, shape (len(t_values), len(x), len(y)).

k0: float = 5.0#
longitudinal_wavefunction(x, t)[source]#

The (force-free) longitudinal packet \(\psi(x,t)\).

Parameters:

  • x (ndarray) – Positions along the beam-propagation axis.

  • t (float) – Time.

Return type:

ndarray

Returns:

numpy.ndarray – Complex-valued.

m: float = 1.0#
mu: float = 1.0#
sigma0: float = 1.0#
sigma_x: float = 2.0#
spin_weights: tuple = (0.5, 0.5)#
transverse_wavefunction(y, t, spin_sign)[source]#

The transverse packet \(\psi_{\pm}(y,t)\) for one spin branch.

Parameters:

  • y (ndarray) – Transverse positions.

  • t (float) – Time.

  • spin_sign (float) – \(+1\) for spin-up, \(-1\) for spin-down.

Return type:

ndarray

Returns:

numpy.ndarray – Complex-valued.

class physicskit.quantum.ThermalState(oscillator, T, k_B=1.0, n_max=200)[source]#

Bases: object

Convenience wrapper bundling an oscillator with a fixed temperature.

Parameters:

  • oscillator (HarmonicOscillator) – The underlying oscillator.

  • T (float) – Temperature.

  • k_B (float) – Boltzmann constant.

  • n_max (int) – Fock-basis truncation order used for the thermal sum.

T: float#
check_normalization(x)[source]#

Integrate position_distribution() over x (should be 1).

Parameters:

x (ndarray)

Return type:

float

Returns:

float

k_B: float = 1.0#
n_max: int = 200#
oscillator: HarmonicOscillator#
populations()[source]#

See HarmonicOscillator.thermal_populations().

Return type:

ndarray

Returns:

numpy.ndarray

position_distribution(x)[source]#

See HarmonicOscillator.thermal_position_distribution().

Parameters:

x (ndarray)

Return type:

ndarray

Returns:

numpy.ndarray

variance()[source]#

See HarmonicOscillator.thermal_position_variance_analytic().

Return type:

float

Returns:

float

class physicskit.quantum.TwinSlit(slit_separation=4.0, slit_width=0.4, k0=10.0, hbar=1.0, m=1.0)[source]#

Bases: object

Young’s double-slit for a massive particle.

Built from two coherent Gaussian point sources at the slit plane, propagated freely to a downstream screen – the quantum interference pattern is the same \(\lvert\psi_1+\psi_2\rvert^2\) mechanism as light, now for massive-particle probability amplitude.

Parameters:

  • slit_separation (float) – Distance between the two slits.

  • slit_width (float) – Width of each slit (sets the Gaussian source envelope).

  • k0 (float) – Central wavenumber of the incident beam.

  • hbar (float) – Value of \(\hbar\) to use.

  • m (float) – Particle mass.

amplitude(x, screen_distance, t=None)[source]#

Fraunhofer-style amplitude on a screen at distance \(L\).

Sums the two slit contributions, each a Gaussian aperture \(e^{-y^2/w^2}\) radiating \(e^{ikr}/\sqrt{r}\) (exact path length \(r\) to the screen point) under its single-slit diffraction envelope \(e^{-(kw\sin\theta)^2/4}\); \(x\) is measured along the screen.

Parameters:

  • x (ndarray) – Screen positions.

  • screen_distance (float) – Distance \(L\) from the slit plane to the screen.

  • t (float | None) – Unused; retained for API symmetry with time-dependent amplitudes.

Return type:

ndarray

Returns:

numpy.ndarray – Complex-valued amplitude.

hbar: float = 1.0#
intensity(x, screen_distance)[source]#

Normalized interference intensity \(\lvert\psi_1+\psi_2\rvert^2\).

Parameters:

  • x (ndarray) – Screen positions.

  • screen_distance (float) – Distance from the slit plane to the screen.

Return type:

ndarray

Returns:

numpy.ndarray – Normalized so that \(\int I\,dx = 1\).

k0: float = 10.0#
m: float = 1.0#
slit_separation: float = 4.0#
slit_width: float = 0.4#
physicskit.quantum.airy_bouncer_energies(alpha=1.0, n_states=8, hbar=1.0, m=1.0)[source]#

Analytic bound-state energies for the half-line gravitational bouncer.

For \(V(x) = \alpha x\) (\(x>0\), infinite floor at \(x=0\)),

\[E_n = -\alpha \left(\frac{\hbar^2}{2 m \alpha}\right)^{1/3} a_n = -\left(\frac{\hbar^2 \alpha^2}{2 m}\right)^{1/3} a_n,\]

where \(a_n < 0\) are the zeros of the Airy function, \(\mathrm{Ai}(a_n) = 0\), and \((\hbar^2/2m\alpha)^{1/3}\) is the bouncer’s natural length scale (see airy_wavefunction()).

Parameters:

  • alpha (float) – Slope of the linear potential.

  • n_states (int) – Number of energy levels to return.

  • hbar (float) – Value of \(\hbar\) to use.

  • m (float) – Particle mass.

Return type:

ndarray

Returns:

numpy.ndarray – The lowest n_states analytic energies, ascending.

physicskit.quantum.airy_wavefunction(x, alpha=1.0, n=0, hbar=1.0, m=1.0)[source]#

Normalized Airy-function eigenstate for the half-line bouncer.

Parameters:

  • x (ndarray) – Grid to evaluate the wavefunction on (values at \(x<0\) are forced to zero).

  • alpha (float) – Slope of the linear potential \(V(x) = \alpha x\).

  • n (int) – Quantum number (0 = ground state).

  • hbar (float) – Value of \(\hbar\) to use.

  • m (float) – Particle mass.

Return type:

ndarray

Returns:

numpy.ndarray – The normalized wavefunction \(\psi_n(x)\).

physicskit.quantum.annihilation_operator(n_max)[source]#

Lowering operator \(\hat a\) in a truncated Fock basis.

Acts as \(\hat a \lvert n\rangle = \sqrt{n}\,\lvert n-1\rangle\) on the truncated basis \(\{\lvert 0\rangle, \dots, \lvert n_\text{max}-1\rangle\}\).

Parameters:

n_max (int) – Dimension of the truncated Fock space.

Return type:

ndarray

Returns:

numpy.ndarray – The \(n_\text{max}\times n_\text{max}\) matrix representation of \(\hat a\).

physicskit.quantum.anticommutator(A, B)[source]#

The anticommutator \(\{\hat A, \hat B\} = \hat A \hat B + \hat B \hat A\).

Parameters:

  • A (ndarray) – Square matrices of matching shape.

  • B (ndarray) – Square matrices of matching shape.

Return type:

ndarray

Returns:

numpy.ndarray – \(\hat A \hat B + \hat B \hat A\).

physicskit.quantum.asymmetric_step_well(width, V_left, V_right)[source]#

A well of unequal left/right wall heights,

\[\begin{split}V(x) = \begin{cases} V_\text{left} & x < -\text{width}/2 \\ 0 & \lvert x\rvert \le \text{width}/2 \\ V_\text{right} & x > \text{width}/2 \end{cases}.\end{split}\]
Parameters:

  • width (float) – Full width of the zero-potential region.

  • V_left (float) – Wall height on the left side.

  • V_right (float) – Wall height on the right side.

Returns:

callable – A function V(x).

physicskit.quantum.asymmetric_well_states(width=2.0, V_left=40.0, V_right=15.0, x_extent=6.0, n_points=1200, n_states=6)[source]#

Bound states of a well with unequal wall heights.

The wavefunctions decay at different rates into the left vs. right forbidden regions (asymmetric evanescent tails), unlike the symmetric finite well.

Parameters:

  • width (float) – Full width of the zero-potential region.

  • V_left (float) – Wall height on the left side.

  • V_right (float) – Wall height on the right side.

  • x_extent (float) – Half-width of the numerical domain, x in [-x_extent, x_extent].

  • n_points (int) – Number of grid points.

  • n_states (int) – Number of bound states to compute.

Return type:

EigenResult

Returns:

physicskit.quantum.core.eigensolvers.EigenResult – The energies and wavefunctions.

physicskit.quantum.bell_state(kind='phi+')[source]#

One of the four maximally entangled two-qubit Bell states.

In the computational basis \(\lvert00\rangle, \lvert01\rangle, \lvert10\rangle, \lvert11\rangle\):

\[\lvert\Phi^\pm\rangle = \tfrac{1}{\sqrt2}(\lvert00\rangle \pm \lvert11\rangle), \qquad \lvert\Psi^\pm\rangle = \tfrac{1}{\sqrt2}(\lvert01\rangle \pm \lvert10\rangle).\]
Parameters:

kind (str) – Which Bell state to return. 'psi-' is the spin singlet.

Return type:

ndarray

Returns:

numpy.ndarray – The 4-component state vector.

Raises:

ValueError – If kind is not one of the four valid labels.

physicskit.quantum.commutator(A, B)[source]#

The commutator \([\hat A, \hat B] = \hat A \hat B - \hat B \hat A\).

Parameters:

  • A (ndarray) – Square matrices of matching shape.

  • B (ndarray) – Square matrices of matching shape.

Return type:

ndarray

Returns:

numpy.ndarray – \(\hat A \hat B - \hat B \hat A\).

physicskit.quantum.creation_operator(n_max)[source]#

Raising operator \(\hat a^\dagger\), the Hermitian conjugate of annihilation_operator().

Parameters:

n_max (int) – Dimension of the truncated Fock space.

Return type:

ndarray

Returns:

numpy.ndarray – The \(n_\text{max}\times n_\text{max}\) matrix representation of \(\hat a^\dagger\).

physicskit.quantum.double_slit_potential(X, Y, wall_x=0.0, wall_thickness=0.5, slit_separation=3.0, slit_width=0.8, wall_height=400.0)[source]#

An opaque wall at \(x=\text{wall\_x}\) pierced by two gaps.

\(V=\text{wall\_height}\) inside the wall slab except within two strips of width slit_width centered at \(y = \pm\text{slit\_separation}/2\), where \(V=0\).

Parameters:

  • X (ndarray) – Coordinate meshgrid (e.g. SplitOperatorSolver2D.X/.Y).

  • Y (ndarray) – Coordinate meshgrid (e.g. SplitOperatorSolver2D.X/.Y).

  • wall_x (float) – Position of the wall along \(x\).

  • wall_thickness (float) – Thickness of the wall along \(x\).

  • slit_separation (float) – Center-to-center distance between the two slits.

  • slit_width (float) – Width of each slit (gap in the wall).

  • wall_height (float) – Potential inside the opaque part of the wall.

Return type:

ndarray

Returns:

numpy.ndarray – \(V(x,y)\) on the grid.

physicskit.quantum.double_well(lam, a)[source]#

The symmetric quartic double well, \(V(x) = \lambda (x^2 - a^2)^2\).

Has minima at \(x = \pm a\), each of depth 0, separated by a barrier of height \(\lambda a^4\) at \(x=0\).

Parameters:

  • lam (float) – Overall strength \(\lambda\) of the quartic potential.

  • a (float) – Half-separation of the two minima.

Returns:

callable – A function V(x).

physicskit.quantum.expectation(operator, state)[source]#

Expectation value \(\langle\psi\rvert \hat A \lvert\psi\rangle\) of an operator in a discrete-basis state vector.

Parameters:

  • operator (ndarray) – The operator matrix \(\hat A\).

  • state (ndarray) – The (not necessarily normalized) state vector \(\lvert\psi\rangle\); flattened if not already 1D.

Return type:

complex

Returns:

complex – \(\langle\psi\rvert \hat A \lvert\psi\rangle\).

physicskit.quantum.finite_well(V0, width)[source]#

A finite square well,

\[\begin{split}V(x) = \begin{cases} 0 & \lvert x\rvert \le \text{width}/2 \\ V_0 & \text{otherwise} \end{cases}.\end{split}\]
Parameters:

  • V0 (float) – Height of the barrier outside the well.

  • width (float) – Full width of the well (zero-potential region).

Returns:

callable – A function V(x).

physicskit.quantum.free_gaussian_wavepacket(x, x0, sigma0, k0)[source]#

A minimum-uncertainty Gaussian wavepacket at \(t=0\).

\[\psi(x,0) = (2\pi\sigma_0^2)^{-1/4}\, \exp\!\left[-\frac{(x-x_0)^2}{4\sigma_0^2}\right] e^{ik_0 x}.\]
Parameters:

  • x (ndarray) – Positions to evaluate at.

  • x0 (float) – Center of the packet.

  • sigma0 (float) – Initial width (\(\Delta x = \sigma_0\)).

  • k0 (float) – Central wavenumber.

Return type:

ndarray

Returns:

numpy.ndarray

physicskit.quantum.gravitational_bouncer_states(alpha=1.0, x_extent=15.0, n_points=2000, n_states=8)[source]#

The ‘quantum bouncer’: a particle falling under uniform gravity.

Solves for \(V(x) = \alpha x\) (\(x>0\)) with a hard floor at \(x=0\). The grid’s lower boundary supplies the floor (\(\psi\) is forced to vanish there), so the domain spans only \(x \in [0, x_\text{extent}]\). Eigenstates are Airy functions \(\mathrm{Ai}(z)\); the energies returned here should closely match airy_bouncer_energies().

Parameters:

  • alpha (float) – Slope of the linear potential (e.g. \(mg\) for gravity).

  • x_extent (float) – Extent of the numerical domain, x in [0, x_extent].

  • n_points (int) – Number of grid points.

  • n_states (int) – Number of bound states to compute.

Return type:

EigenResult

Returns:

physicskit.quantum.core.eigensolvers.EigenResult – The energies and wavefunctions.

physicskit.quantum.harmonic_well(m=1.0, omega=1.0)[source]#

The quantum harmonic oscillator potential, \(V(x) = \tfrac{1}{2} m \omega^2 x^2\).

Parameters:

  • m (float) – Particle mass.

  • omega (float) – Angular frequency \(\omega\).

Returns:

callable – A function V(x).

physicskit.quantum.infinite_well()[source]#

The infinite square well: \(V(x) = 0\) inside the domain.

The grid boundaries passed to NumerovSolver act as the infinite walls, since the wavefunction is forced to vanish there.

Returns:

callable – A function V(x) returning zeros the same shape as x.

physicskit.quantum.is_hermitian(A, atol=1e-10)[source]#

Check whether a matrix is Hermitian, \(\hat A = \hat A^\dagger\).

Parameters:
Return type:

bool

Returns:

bool – True if A equals its conjugate transpose within atol.

physicskit.quantum.linear_gravitational_well(alpha)[source]#

The V-shaped ‘quantum bouncer’ well, \(V(x) = \alpha \lvert x\rvert\).

A particle in this potential is the quantum analogue of a ball bouncing under uniform gravity (\(\alpha = mg\) on the half-line); its eigenstates are Airy functions (see physicskit.quantum.chapters.potentials.airy_bouncer_energies()).

Parameters:

alpha (float) – Slope of the linear potential (e.g. \(mg\) for gravity).

Returns:

callable – A function V(x).

physicskit.quantum.linear_stark_shift(n, n1, n2, m, F, a0=1.0)[source]#

First-order Stark shift of a hydrogenic parabolic state.

For \(\lvert n, n_1, n_2, m\rangle\) under a uniform field \(F\) along \(z\),

\[\Delta E = \frac{3}{2}\, n\, (n_1 - n_2)\, a_0 F,\]

valid when \(n_1 + n_2 + \lvert m\rvert + 1 = n\). This linear-in- \(F\) splitting is a hallmark of hydrogen’s “accidental” \(l\)-degeneracy (removed by any core penetration).

Parameters:

  • n (int) – Principal quantum number.

  • n1 (int) – Parabolic quantum numbers.

  • n2 (int) – Parabolic quantum numbers.

  • m (int) – Magnetic quantum number.

  • F (float) – Electric field strength.

  • a0 (float) – Bohr radius.

Return type:

float

Returns:

float

Raises:

ValueError – If n1 + n2 + abs(m) + 1 != n.

physicskit.quantum.momentum_operator(n_max, m=1.0, omega=1.0, hbar=1.0)[source]#

Harmonic-oscillator momentum operator in the truncated Fock basis,

\[\hat p = i\sqrt{\frac{\hbar m \omega}{2}} \left(\hat a^\dagger - \hat a\right).\]
Parameters:

  • n_max (int) – Dimension of the truncated Fock space.

  • m (float) – Particle mass.

  • omega (float) – Angular frequency of the oscillator.

  • hbar (float) – Value of \(\hbar\) to use.

Return type:

ndarray

Returns:

numpy.ndarray – The \(n_\text{max}\times n_\text{max}\) matrix representation of \(\hat p\).

physicskit.quantum.number_operator(n_max)[source]#

Number operator \(\hat N = \hat a^\dagger \hat a\), diagonal with eigenvalues \(0, 1, \dots, n_\text{max}-1\).

Parameters:

n_max (int) – Dimension of the truncated Fock space.

Return type:

ndarray

Returns:

numpy.ndarray – The \(n_\text{max}\times n_\text{max}\) matrix representation of \(\hat N\).

physicskit.quantum.orbital_superposition_density(orb_a, orb_b, r, theta, phi, t, ca=0.7071067811865476, cb=0.7071067811865476)[source]#

The probability density of orbital_superposition_psi().

Parameters:
Return type:

ndarray

Returns:

numpy.ndarray

physicskit.quantum.orbital_superposition_psi(orb_a, orb_b, r, theta, phi, t, ca=0.7071067811865476, cb=0.7071067811865476)[source]#

A coherent superposition of two hydrogen eigenstates at time \(t\).

\[\psi(t) = c_a\,\psi_a\, e^{-iE_a t} + c_b\,\psi_b\, e^{-iE_b t},\]

beating at the Bohr frequency \(\omega_{ab}=E_a-E_b\) between the two stationary densities – genuine dynamics, unlike either eigenstate alone (whose density is static).

Parameters:

  • orb_a (HydrogenOrbital) – The two eigenstates in the superposition (should differ in energy for the beating to be visible).

  • orb_b (HydrogenOrbital) – The two eigenstates in the superposition (should differ in energy for the beating to be visible).

  • r (ndarray) – Coordinates to evaluate at (e.g. from a Cartesian grid, as in orbital_density_grid()).

  • theta (ndarray) – Coordinates to evaluate at (e.g. from a Cartesian grid, as in orbital_density_grid()).

  • phi (ndarray) – Coordinates to evaluate at (e.g. from a Cartesian grid, as in orbital_density_grid()).

  • t (float) – Time.

  • ca (complex) – Superposition amplitudes (renormalized internally).

  • cb (complex) – Superposition amplitudes (renormalized internally).

Return type:

ndarray

Returns:

numpy.ndarray – Complex-valued \(\psi(r,\theta,\phi,t)\).

physicskit.quantum.position_operator(n_max, m=1.0, omega=1.0, hbar=1.0)[source]#

Harmonic-oscillator position operator in the truncated Fock basis,

\[\hat x = \sqrt{\frac{\hbar}{2 m \omega}} \left(\hat a + \hat a^\dagger\right).\]
Parameters:

  • n_max (int) – Dimension of the truncated Fock space.

  • m (float) – Particle mass.

  • omega (float) – Angular frequency of the oscillator.

  • hbar (float) – Value of \(\hbar\) to use.

Return type:

ndarray

Returns:

numpy.ndarray – The \(n_\text{max}\times n_\text{max}\) matrix representation of \(\hat x\).

physicskit.quantum.propagate_double_slit(x, y, x0=-8.0, k0=8.0, sigma_x=0.7, sigma_y=4.0, wall_x=0.0, wall_thickness=0.5, slit_separation=3.0, slit_width=0.8, wall_height=400.0, hbar=1.0, m=1.0, dt=0.0001, n_steps=4000, save_every=40)[source]#

Propagate a wavepacket through double_slit_potential().

A genuinely time-propagated double-slit: a Gaussian packet, elongated along \(y\) (to illuminate both slits coherently) and moving along \(+x\) with wavenumber k0, is evolved with SplitOperatorSolver2D through the two-slit wall, producing the interference pattern building up downstream from genuine wave dynamics rather than the far-field Fraunhofer formula (TwinSlit).

Parameters:
Returns:

  • X, Y (numpy.ndarray) – The coordinate meshgrid.

  • frames (numpy.ndarray) – Complex wavefunction snapshots, shape (n_saved, len(x), len(y)).

  • times (numpy.ndarray) – The time of each snapshot.

physicskit.quantum.radial_wavefunction(n, l, r, Z=1, a0=1.0)[source]#

The hydrogenic radial wavefunction \(R_{nl}(r)\).

Normalized so that \(\int_0^\infty \lvert R_{nl}(r)\rvert^2 r^2\,dr = 1\).

Parameters:

  • n (int) – Principal quantum number.

  • l (int) – Orbital angular momentum quantum number (\(0 \le l < n\)).

  • r (ndarray) – Radii to evaluate at.

  • Z (int) – Nuclear charge.

  • a0 (float) – Bohr radius.

Return type:

ndarray

Returns:

numpy.ndarray

Raises:

ValueError – If l >= n.

physicskit.quantum.spherical_harmonic(l, m, theta, phi)[source]#

The spherical harmonic \(Y_l^m(\theta,\phi)\).

Parameters:

  • l (int) – Orbital angular momentum quantum number.

  • m (int) – Magnetic quantum number (\(\lvert m\rvert \le l\)).

  • theta (ndarray) – Polar angle, in \([0, \pi]\).

  • phi (ndarray) – Azimuthal angle, in \([0, 2\pi)\).

Return type:

ndarray

Returns:

numpy.ndarray – Complex-valued \(Y_l^m(\theta,\phi)\).

physicskit.quantum.spin_operator(axis, s=0.5, hbar=1.0)[source]#

Angular-momentum operator \(\hat S_\text{axis}\) for spin quantum number s.

For \(s=1/2\) this reduces to \((\hbar/2)\,\sigma_\text{axis}\). For general \(s\), the \((2s+1)\)-dimensional matrix representation is built from the ladder operators

\[\hat S_\pm \lvert s, m\rangle = \hbar\sqrt{s(s+1) - m(m\pm 1)}\, \lvert s, m\pm 1\rangle,\]

with \(\hat S_x = (\hat S_+ + \hat S_-)/2\) and \(\hat S_y = (\hat S_+ - \hat S_-)/2i\).

Parameters:

  • axis (str) – Which component (or ladder operator) to return.

  • s (float) – Spin quantum number (may be integer or half-integer).

  • hbar (float) – Value of \(\hbar\) to use.

Return type:

ndarray

Returns:

numpy.ndarray – The \((2s+1)\times(2s+1)\) operator matrix, in the \(\lvert s, m\rangle\) basis ordered \(m = s, s-1, \dots, -s\).

Raises:

ValueError – If axis is not one of 'x', 'y', 'z', '+', '-'.

physicskit.quantum.stark_n2_quartet(F, a0=1.0)[source]#

The linear Stark splitting of hydrogen’s \(n=2\) shell.

The 4-fold degenerate \(n=2\) level (2s, 2p0, \(2p_{\pm1}\)) splits under a field \(F\) into three levels with shifts \(\{-3a_0F, 0, 0, +3a_0F\}\) – the classic linear Stark quartet, since only the \(m=0\) pair (2s/2p0) is coupled by the perturbation.

Parameters:

  • F (float) – Electric field strength.

  • a0 (float) – Bohr radius.

Return type:

list

Returns:

list of StarkLevel – Sorted by energy shift.

physicskit.quantum.zeeman_spectrum(l, B, mu_B=1.0, g_s=2.0)[source]#

All \((m_l, m_s)\) sublevel shifts for orbital number \(l\), spin-1/2.

Parameters:

  • l (int) – Orbital angular momentum quantum number.

  • B (float) – Magnetic field strength.

  • mu_B (float) – Bohr magneton.

  • g_s (float) – Electron spin g-factor.

Return type:

ndarray

Returns:

numpy.ndarray – Sorted array of \(2(2l+1)\) energy shifts.

physicskit.quantum.zeeman_splitting(m_l, m_s, B, mu_B=1.0, g_s=2.0, hbar=1.0)[source]#

First-order (weak-field, no spin-orbit) Zeeman energy shift.

\[\Delta E = \mu_B B\,(m_l + g_s m_s),\]

in the \(\lvert l, m_l; s, m_s\rangle\) basis.

Parameters:

  • m_l (int) – Orbital magnetic quantum number.

  • m_s (float) – Spin magnetic quantum number.

  • B (float) – Magnetic field strength.

  • mu_B (float) – Bohr magneton.

  • g_s (float) – Electron spin g-factor.

  • hbar (float) – Value of \(\hbar\) (unused in this first-order formula; kept for API symmetry with other perturbation functions).

Return type:

float

Returns:

float

On this page
Show Source