physicskit.quantum #

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:

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:

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:

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 | None) – 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:

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:

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:

Returns:

float

state: ndarray | None = None#

class physicskit.quantum.Box2DEigenstate(nx, ny, energy)[source]#

Bases: object

One eigenstate of RectangularBox2D.

Parameters:

nx (int)
ny (int)
energy (float)

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:

X (ndarray) – Coordinate meshgrid (e.g. from RectangularBox2D.grid()).
Y (ndarray) – Coordinate meshgrid (e.g. from RectangularBox2D.grid()).
Lx (float) – Box dimensions.
Ly (float) – Box dimensions.

Return type:

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:

m (int)
n (int)
k (float)
energy (float)

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:

r (ndarray) – Polar coordinate meshgrid (e.g. from CircularBox2D.polar_grid()).
phi (ndarray) – Polar coordinate meshgrid (e.g. from CircularBox2D.polar_grid()).
R (float) – Dot radius (used to normalize \(N_{mn}\)).

Return type:

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:

x (ndarray)
energies (ndarray)
wavefunctions (ndarray)
splitting (float)
tunneling_period (float)

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:

result (DoubleWellResult) – Output of solve().
t (ndarray) – Times to evaluate at.
side (str) – Which well the state starts localized in.

Return type:

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:

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:

result (DoubleWellResult) – Output of solve().
t (ndarray) – Times to evaluate at.
side (str) – Which well the state starts localized in.

Return type:

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:

result (DoubleWellResult) – Output of solve().
t (ndarray) – Times to evaluate at.
side (str) – Which well the state starts localized in.

Return type:

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:

x (ndarray)
energies (ndarray)
wavefunctions (ndarray)

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:

result (EigenResult) – Output of bound_states().
n (int) – Which bound state to return (0 = ground state).

Return type:

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:

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:

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:

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:

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.GaussianS(exponents, coefficients, center)[source]#

Bases: object

A contracted s-type Gaussian \(\phi(\mathbf r) = \sum_p d_p\,(2\alpha_p/\pi)^{3/4} e^{-\alpha_p|\mathbf r - \mathbf A|^2}\).

Parameters:

exponents (ndarray) – Primitive exponents \(\alpha_p\).
coefficients (ndarray) – Contraction coefficients \(d_p\) of the normalized primitives.
center (ndarray) – Center \(\mathbf A\) (bohr), a 3-vector.

center: ndarray#

coefficients: ndarray#

exponents: ndarray#

property weights: ndarray#: Coefficients times primitive normalization constants.

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:

alpha (complex) – Coherent-state amplitude.
x (ndarray) – Positions to evaluate at.
t (float) – Time.
n_max (int | None) – Fock-basis truncation order; defaults to coherent_n_max().

Return type:

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:

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:

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:

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:

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:

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:

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:

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:

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:

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.HartreeFockResult(energy, electronic_energy, nuclear_repulsion, orbital_energies, coefficients, density, n_iterations, converged)[source]#

Bases: object

Output of restricted_hartree_fock().

Variables:

energy (float) – Total energy (electronic plus nuclear repulsion), hartree.
electronic_energy (float) – \(E_0 = \tfrac12\sum_{\mu\nu}P_{\nu\mu}(H^\mathrm{core}_{\mu\nu} + F_{\mu\nu})\).
nuclear_repulsion (float) – \(\sum_{A<B} Z_AZ_B/R_{AB}\).
orbital_energies (numpy.ndarray) – Canonical orbital energies \(\varepsilon_i\), ascending.
coefficients (numpy.ndarray) – MO coefficients (columns), \(C^TSC = I\).
density (numpy.ndarray) – Closed-shell density matrix \(P = 2C_\mathrm{occ}C_\mathrm{occ}^T\).
n_iterations (int) – SCF iterations used.
converged (bool)

Parameters:

energy (float)
electronic_energy (float)
nuclear_repulsion (float)
orbital_energies (ndarray)
coefficients (ndarray)
density (ndarray)
n_iterations (int)
converged (bool)

coefficients: ndarray#

converged: bool#

density: ndarray#

electronic_energy: float#

energy: float#

n_iterations: int#

nuclear_repulsion: float#

orbital_energies: ndarray#

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:

theta (ndarray)
phi (ndarray)

Return type:

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:

Returns:

float

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

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

Parameters:

r (ndarray)
theta (ndarray)
phi (ndarray)

Return type:

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:

Returns:

float

n: int#

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

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

Parameters:

r (ndarray)
theta (ndarray)
phi (ndarray)

Return type:

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.MPS(B, S)[source]#

Bases: object

An open-boundary matrix product state in right-canonical form.

Variables:

B (list of ndarray) – Site tensors of shape (chi_left, 2, chi_right), each right-canonical.
S (list of ndarray) – Singular values on the bond to the left of each site (S[0] = [1]).

Parameters:

B (list[NDArray])
S (list[NDArray])

Examples

>>> psi = MPS.product_state([0, 1, 0, 1])  # up, down, up, down
>>> v = psi.to_dense()
>>> int(np.argmax(np.abs(v))), float(abs(v).max())
(10, 1.0)

apply_two_site(j, U, chi_max, eps=1e-12)[source]#

Apply a \(4\times4\) gate to sites j, j+1 and truncate; returns the discarded weight.

Return type:

Parameters:

j (int)
U (NDArray)
chi_max (int)
eps (float)

property bond_dimensions: list[int]#: Bond dimensions \(\chi\) of the N - 1 internal bonds.

bond_energies(H_bonds)[source]#

\(\langle h_{j,j+1}\rangle\) for every bond.

Return type:: NDArray[double]
Parameters:: H_bonds (list[NDArray])

canonicalize()[source]#

Restore the normalized right-canonical form with exact Schmidt values on every bond.

A right-to-left sweep makes every tensor right-orthonormal, then a left-to-right sweep of SVDs reads off the Schmidt values \(S_{j}\) and sets \(B_j = S_j^{-1} A_j S_{j+1}\) (Vidal’s relation). Non-unitary gates, as in imaginary-time evolution, break the canonical form that bond_energies() and the TEBD truncation rely on.

Return type:: None

copy()[source]#

Deep copy.

Return type:: MPS

entanglement_entropy()[source]#

Von Neumann entropy \(-\sum_\alpha s_\alpha^2\ln s_\alpha^2\) of each internal bond.

Return type:: NDArray[double]

expectation_local(op)[source]#

\(\langle O_j \rangle\) of a single-site operator on every site.

Return type:: NDArray[double]
Parameters:: op (NDArray)

classmethod product_state(spins)[source]#

Product state with physical index spins[j] (0 up, 1 down) on site j.

Return type:: MPS

to_dense()[source]#

The full \(2^N\) state vector, in the basis order of physicskit.condensed.spin_chains.

Return type:: NDArray[cdouble]

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.QuantumCircuit(n_qubits)[source]#

Bases: object

A list of gates on n_qubits qubits, simulated on a state vector.

Gates are appended with the methods below, which return the circuit so calls can be chained, and run() applies them in order.

Parameters:: n_qubits (int)

Examples

A Bell pair:

>>> qc = QuantumCircuit(2).h(0).cx(0, 1)
>>> np.round(qc.run(), 4).tolist()
[(0.7071+0j), 0j, 0j, (0.7071+0j)]

append(other, offset=0)[source]#

Append the gates of other with its qubit j mapped to j + offset.

Return type:

Parameters:

other (QuantumCircuit)
offset (int)

cphase(theta, control, target)[source]#

Controlled phase \(\mathrm{diag}(1, 1, 1, e^{i\theta})\).

Return type:

Parameters:

theta (float)
control (int)
target (int)

cx(control, target)[source]#

Controlled NOT.

Return type:

Parameters:

control (int)
target (int)

diagonal(phases, qubits=None)[source]#

Diagonal unitary \(|x\rangle \to e^{i\phi_x}|x\rangle\) on qubits (default all).

Return type:

Parameters:

phases (ArrayLike)
qubits (tuple[int, ...] | None)

gate(U, *qubits, name='U')[source]#

Append a \(2^k \times 2^k\) unitary acting on qubits (first is most significant).

Return type:

Parameters:

U (ArrayLike)
qubits (int)
name (str)

gate_count()[source]#

Number of gates of each kind.

Return type:: dict[str, int]

h(q)[source]#

Hadamard.

Return type:: QuantumCircuit
Parameters:: q (int)

inverse()[source]#

The adjoint circuit: gates in reverse order, each inverted.

Return type:: QuantumCircuit

permutation(perm, qubits=None)[source]#

Reversible classical gate \(|x\rangle \to |\pi(x)\rangle\).

Return type:

Parameters:

perm (ArrayLike)
qubits (tuple[int, ...] | None)

phase(theta, q)[source]#

Phase gate \(\mathrm{diag}(1, e^{i\theta})\).

Return type:

Parameters:

theta (float)
q (int)

static probabilities(psi, n_qubits, qubits=None)[source]#

Measurement probabilities of qubits (default all), marginalizing the rest.

Parameters:

psi (ArrayLike)
n_qubits (int)
qubits (tuple[int, ...] | None)

Return type:

NDArray[double]

Returns:

ndarray of float, shape (2**len(qubits),)

run(initial=0)[source]#

Apply the circuit to a basis state (an integer) or a state vector.

Parameters:: initial (ArrayLike | int) – Basis-state index, or a vector of \(2^n\) amplitudes.
Return type:: NDArray[cdouble]
Returns:: ndarray of complex, shape (2**n,)

static sample(probabilities, shots, seed=None)[source]#

Counts of each outcome in shots projective measurements.

Return type:

NDArray[int_]

Parameters:

probabilities (ArrayLike)
shots (int)
seed (int | None)

swap(a, b)[source]#

Swap two qubits.

Return type:

Parameters:

a (int)
b (int)

unitary()[source]#

The full \(2^n \times 2^n\) matrix of the circuit (small n only).

Return type:: NDArray[cdouble]

x(q)[source]#

Pauli X (NOT).

Return type:: QuantumCircuit
Parameters:: q (int)

z(q)[source]#

Pauli Z.

Return type:: QuantumCircuit
Parameters:: q (int)

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:

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:

x (ndarray) – Positions to evaluate at.
t (float) – Time.
coeffs (ndarray) – Expansion coefficients from eigenbasis_coefficients().
psi0 (ndarray) – The initial state sampled on x.

Return type:

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:

x (ndarray) – Positions to evaluate at.
t (float) – Time.
coeffs (ndarray) – Expansion coefficients from eigenbasis_coefficients().

Return type:

Returns:

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

class physicskit.quantum.QuantumTrajectories(H, c_ops)[source]#

Bases: object

Monte Carlo wavefunction (quantum-jump) unravelling of a Lindblad equation.

Each trajectory evolves an unnormalized state under the non-Hermitian effective Hamiltonian

\[H_\mathrm{eff} = H - \frac{i}{2} \sum_k L_k^\dagger L_k,\]

whose squared norm decays from 1. A uniform random number \(r\) is drawn, and when \(\lVert\psi\rVert^2\) falls below \(r\) a quantum jump \(\psi \to L_k\psi\) is applied, channel \(k\) chosen with probability \(\propto \lVert L_k \psi \rVert^2\); the state is renormalized and a fresh \(r\) drawn (the waiting-time formulation of Dalibard, Castin and Molmer 1992 and Dum, Zoller and Ritsch, Phys. Rev. A 45, 4879 (1992); review: Plenio and Knight, Rev. Mod. Phys. 70, 101 (1998), Sec. III). The ensemble average of \(|\psi\rangle\langle\psi|\) obeys the Lindblad equation.

The no-jump propagator \(e^{-iH_\mathrm{eff}\,dt}\) is applied exactly, so the only discretization error is the \(O(dt)\) resolution of each jump time. All trajectories are advanced together as one vectorized array.

Parameters:

H (ndarray) – Hamiltonian, shape (d, d).
c_ops (Sequence[ndarray]) – Collapse operators \(L_k\).

Examples

Starting in the excited state with pure decay, every trajectory eventually makes exactly one jump:

>>> import numpy as np
>>> traj = QuantumTrajectories(np.zeros((2, 2)), [sigma_minus()])
>>> res = traj.run(np.array([0, 1]), np.linspace(0, 20, 5), n_traj=50, dt=0.01,
...                rng=np.random.default_rng(1))
>>> {len(j) for j in res.jump_times}
{1}

H: ndarray#

c_ops: Sequence[ndarray]#

property effective_hamiltonian: ndarray#: The non-Hermitian \(H_\mathrm{eff} = H - \frac{i}{2}\sum_k L_k^\dagger L_k\).

run(psi0, times, n_traj=500, dt=None, rng=None)[source]#

Generate an ensemble of quantum trajectories.

Parameters:

psi0 (ndarray) – Initial pure state, shape (d,); normalized internally.
times (ndarray) – Increasing output times; times[0] is the start time.
n_traj (int) – Number of trajectories.
dt (float | None) – Maximum internal step. Defaults to one hundredth of the fastest jump timescale \(1/\sum_k \lVert L_k^\dagger L_k\rVert\), capped by the output spacing.
rng (Generator | None) – Random generator (for reproducibility).

Return type:

TrajectoryResult

Returns:

TrajectoryResult

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:

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:

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:

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)
k1 (float)
k2 (complex)
T (float)
R (float)

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:

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:

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:

X (ndarray) – Coordinate meshgrid (e.g. from grid()).
Y (ndarray) – Coordinate meshgrid (e.g. from grid()).

Return type:

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:

n1 (int)
n2 (int)
m (int)
shift (float)

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:

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:

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:

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:

See HarmonicOscillator.thermal_position_distribution().

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

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.TrajectoryResult(times, states, jump_times=<factory>, jump_channels=<factory>)[source]#

Bases: object

Output of QuantumTrajectories.run().

Variables:

times (numpy.ndarray) – Output times, shape (n_t,).
states (numpy.ndarray) – Normalized state of every trajectory at every output time, shape (n_traj, n_t, d).
jump_times (list of numpy.ndarray) – For each trajectory, the times at which a quantum jump occurred.
jump_channels (list of numpy.ndarray) – For each trajectory, the index \(k\) of the collapse operator \(L_k\) applied at each jump.

Parameters:

times (ndarray)
states (ndarray)
jump_times (list)
jump_channels (list)

density_matrices()[source]#

Ensemble-averaged density matrix \(\overline{|\psi\rangle\langle\psi|}\) at each time.

Return type:: ndarray
Returns:: numpy.ndarray – Shape (n_t, d, d).

expectation(op, average=True)[source]#

\(\langle\psi|A|\psi\rangle\) along each trajectory, or its ensemble average.

Parameters:

op (ndarray) – Operator \(A\), shape (d, d).
average (bool) – If True, average over trajectories.

Return type:

Returns:

numpy.ndarray – Real part of the expectation value, shape (n_t,) if average else (n_traj, n_t).

jump_channels: list#

jump_times: list#

property n_trajectories: int#: Number of trajectories.

states: ndarray#

times: ndarray#

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:

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:

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:

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:

Returns:

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

physicskit.quantum.amplitude_damping_kraus(gamma)[source]#

Kraus operators of the qubit amplitude-damping channel.

\[\begin{split}E_0 = \begin{pmatrix} 1 & 0 \\ 0 & \sqrt{1-\gamma} \end{pmatrix}, \qquad E_1 = \begin{pmatrix} 0 & \sqrt{\gamma} \\ 0 & 0 \end{pmatrix}\end{split}\]

(Nielsen and Chuang, Eq. 8.108): the excited state \(|1\rangle\) decays to \(|0\rangle\) with probability \(\gamma\). Lindblad decay at rate \(\Gamma\) for time \(t\) is exactly this channel with \(\gamma = 1 - e^{-\Gamma t}\).

Parameters:: gamma (float) – Decay probability, \(0 \le \gamma \le 1\).
Return type:: list[ndarray]
Returns:: list of numpy.ndarray

Examples

>>> import numpy as np
>>> rho = apply_kraus(np.diag([0.0, 1.0]), amplitude_damping_kraus(0.3))
>>> np.real(np.diag(rho)).round(3)
array([0.3, 0.7])

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:

Returns:

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

physicskit.quantum.apply_kraus(rho, kraus)[source]#

Apply a channel in Kraus form, \(\rho \mapsto \sum_i E_i \rho E_i^\dagger\).

Parameters:

rho (ndarray) – Density matrix, shape (d, d).
kraus (Sequence[ndarray]) – Kraus operators \(E_i\).

Return type:

Returns:

numpy.ndarray – The output density matrix.

Examples

>>> import numpy as np
>>> rho = np.diag([0.5, 0.5]).astype(complex)
>>> np.allclose(apply_kraus(rho, dephasing_kraus(0.2)), rho)
True

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.born_amplitude(potential, q, mass=1.0)[source]#

First Born approximation to the scattering amplitude of a central potential.

\[f_B(q) = -\frac{m}{2\pi}\int d^3r\, e^{-i\mathbf q\cdot\mathbf r} V(r) = -\frac{2m}{q}\int_0^\infty r\,V(r)\sin(qr)\,dr\]

(Born 1926; Sakurai, Modern Quantum Mechanics, Eq. 7.2.10), with \(q \to 0\) limit \(-2m\int_0^\infty r^2 V\,dr\). The oscillatory integral is done with QUADPACK’s Fourier-integral routine (scipy.integrate.quad() with weight="sin"), which needs \(rV(r)\) to decay at large \(r\).

Parameters:

potential (Callable[[ndarray], ndarray]) – \(V(r)\), vectorized over r.
q (ndarray) – Momentum transfer(s), \(q \ge 0\).
mass (float) – Particle (reduced) mass.

Return type:

Returns:

numpy.ndarray – Real Born amplitude \(f_B(q)\).

Examples

A Gaussian well \(V = V_0 e^{-r^2}\) has \(f_B = -\tfrac{\sqrt\pi}{2} m V_0\,e^{-q^2/4}\):

>>> import numpy as np
>>> f = born_amplitude(lambda r: -0.3 * np.exp(-r**2), 1.2)
>>> bool(np.isclose(f, 0.5 * np.sqrt(np.pi) * 0.3 * np.exp(-0.36)))
True

physicskit.quantum.born_phase_shifts(potential, k, l_max, mass=1.0, r_max=30.0)[source]#

Phase shifts in the Born approximation, \(\delta_\ell \approx -2mk\int_0^\infty r^2 V(r)\,j_\ell(kr)^2\,dr\).

The small-\(V\) limit of the variable-phase equation (Sakurai, Eq. 7.6.22); valid when every \(|\delta_\ell| \ll 1\).

Parameters:

potential (Callable[[ndarray], ndarray]) – \(V(r)\).
k (float) – Wavenumber.
l_max (int) – Highest partial wave.
mass (float) – Particle (reduced) mass.
r_max (float) – Outer integration radius.

Return type:

Returns:

numpy.ndarray – Shape (l_max + 1,).

Examples

>>> import numpy as np
>>> weak = lambda r: -0.01 * np.exp(-r**2)
>>> d_born = born_phase_shifts(weak, k=1.0, l_max=2)
>>> d_exact = partial_wave_phase_shifts(weak, k=1.0, l_max=2, r_max=8.0)
>>> bool(np.allclose(d_born, d_exact, rtol=0.02))
True

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:

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.dephasing_kraus(p)[source]#

Kraus operators of the qubit dephasing (phase-flip) channel.

\[E_0 = \sqrt{1-p}\,I, \qquad E_1 = \sqrt{p}\,\sigma_z\]

(Nielsen and Chuang, Eq. 8.3.3): populations are untouched and the coherences shrink by \(1 - 2p\). Pure Lindblad dephasing \(L = \sqrt{\gamma_\varphi/2}\,\sigma_z\) for time \(t\) is this channel with \(p = (1 - e^{-\gamma_\varphi t})/2\).

Parameters:: p (float) – Phase-flip probability, \(0 \le p \le 1\).
Return type:: list[ndarray]
Returns:: list of numpy.ndarray

Examples

>>> import numpy as np
>>> plus = 0.5 * np.ones((2, 2))
>>> round(float(apply_kraus(plus, dephasing_kraus(0.1))[0, 1].real), 12)
0.4

physicskit.quantum.deutsch_jozsa_circuit(f)[source]#

Deutsch-Jozsa circuit for a Boolean function given by its truth table.

\(H^{\otimes n}\), the phase oracle \(|x\rangle \to (-1)^{f(x)}|x\rangle\) (one query), and \(H^{\otimes n}\) again. The probability of reading \(|0\dots0\rangle\) is \(|2^{-n}\sum_x(-1)^{f(x)}|^2\): 1 if \(f\) is constant and 0 if it is balanced (Deutsch and Jozsa, Proc. R. Soc. A 439, 553 (1992); Cleve et al., Proc. R. Soc. A 454, 339 (1998)).

Parameters:: f (ArrayLike)
Return type:: QuantumCircuit
Returns:: QuantumCircuit

Examples

>>> balanced = [0, 1, 1, 0, 1, 0, 0, 1]
>>> psi = deutsch_jozsa_circuit(balanced).run()
>>> round(float(abs(psi[0]) ** 2), 12)
0.0

physicskit.quantum.dirac_hamiltonian(p, m=1.0)[source]#

Free Dirac Hamiltonian \(H = \boldsymbol\alpha\cdot\mathbf p + \beta m\) for momentum \(\mathbf p\).

With \(\beta = \gamma^0\) and \(\alpha_i = \gamma^0\gamma^i\) (Dirac, Proc. R. Soc. A 117, 610 (1928)). Its eigenvalues are \(\pm\sqrt{p^2 + m^2}\), each doubly degenerate (spin).

Parameters:

p (ndarray) – Momentum 3-vector.
m (float) – Mass.

Return type:

Returns:

numpy.ndarray – Shape (4, 4), Hermitian.

Examples

>>> import numpy as np
>>> np.linalg.eigvalsh(dirac_hamiltonian([0.0, 0.0, 0.75])).round(12)
array([-1.25, -1.25,  1.25,  1.25])

physicskit.quantum.dirac_hydrogen_energy(n, j, Z=1.0, alpha=0.0072973525643)[source]#

Exact Dirac-Coulomb energy level, in units of \(mc^2\) (rest energy included).

\[\frac{E_{nj}}{mc^2} = \left[1 + \left(\frac{Z\alpha} {n - (j + \tfrac12) + \sqrt{(j + \tfrac12)^2 - (Z\alpha)^2}} \right)^2\right]^{-1/2}\]

(Darwin, Proc. R. Soc. A 118, 654 (1928); Gordon, Z. Phys. 48, 11 (1928); Bjorken and Drell, Eq. 4.29). It depends on \(n\) and \(j\) only, so \(2S_{1/2}\) and \(2P_{1/2}\) stay degenerate (the Lamb shift lifts this). The ground state is \(E_{1,1/2} = mc^2\sqrt{1 - (Z\alpha)^2}\).

Parameters:

n (int) – Principal quantum number, \(n \ge 1\).
j (float) – Total angular momentum, \(\tfrac12 \le j \le n - \tfrac12\).
Z (float) – Nuclear charge.
alpha (float) – Fine-structure constant.

Return type:

Returns:

float – \(E/(mc^2)\).

Raises:

ValueError – If \(j\) is out of range or \(Z\alpha \ge j + 1/2\).

Examples

>>> import numpy as np
>>> bool(np.isclose(dirac_hydrogen_energy(1, 0.5, Z=50), np.sqrt(1 - (50 * 7.2973525693e-3) ** 2)))
True

physicskit.quantum.dirac_plane_wave_spinor(p, m=1.0, spin=1, energy_sign=1)[source]#

Free Dirac plane-wave spinor, eigenvector of dirac_hamiltonian() with energy \(\pm E\).

\[\begin{split}u_s(\mathbf p) = \sqrt{E + m}\begin{pmatrix} \chi_s \\ \dfrac{\boldsymbol\sigma\cdot\mathbf p}{E + m}\chi_s \end{pmatrix}, \qquad w_s(\mathbf p) = \sqrt{E + m}\begin{pmatrix} -\dfrac{\boldsymbol\sigma\cdot\mathbf p}{E + m}\chi_s \\ \chi_s \end{pmatrix},\end{split}\]

for energies \(+E\) and \(-E\) respectively, \(E = \sqrt{p^2+m^2}\), normalized to \(\psi^\dagger\psi = 2E\) (Bjorken and Drell, Eq. 3.7; for positive energy \(\bar u u = 2m\)).

Parameters:

p (ndarray) – Momentum 3-vector.
m (float) – Mass.
spin (int) – Two-component spinor \(\chi_{+} = (1, 0)^T\) or \(\chi_{-} = (0, 1)^T\).
energy_sign (int) – Positive- or negative-energy solution.

Return type:

Returns:

numpy.ndarray – Shape (4,), complex.

Examples

>>> import numpy as np
>>> p = np.array([0.3, -0.2, 0.5])
>>> u = dirac_plane_wave_spinor(p)
>>> E = np.sqrt(p @ p + 1)
>>> bool(np.allclose(dirac_hamiltonian(p) @ u, E * u))
True

physicskit.quantum.dirac_step_scattering(E, V0, m=1.0, group_velocity=True)[source]#

Reflection and transmission probabilities of a Dirac particle at a potential step.

A spin-up plane wave of energy \(E > m\) hits \(V(z) = V_0 \Theta(z)\). Matching the Dirac spinor at \(z = 0\) gives

\[\kappa = \frac{q}{p}\,\frac{E + m}{E - V_0 + m}, \qquad R = \left|\frac{1 - \kappa}{1 + \kappa}\right|^2, \qquad T = \frac{4\,\mathrm{Re}\,\kappa}{|1 + \kappa|^2},\]

with \(p = \sqrt{E^2 - m^2}\) and \(q = \sqrt{(E - V_0)^2 - m^2}\), so that \(R + T = 1\) (Bjorken and Drell, Sec. 3.3; Greiner, Relativistic Quantum Mechanics, Sec. 13.1). For \(|E - V_0| < m\), \(q\) is imaginary and \(R = 1\). In the Klein zone \(V_0 > E + m\) the transmitted wave is a negative-energy (antiparticle) state: choosing \(q\) so that its group velocity \(q/(E - V_0)\) points away from the step gives \(\kappa > 0\) and \(0 < R < 1\) – a finite transmission even as \(V_0 \to \infty\), Klein’s paradox (Klein, Z. Phys. 53, 157 (1929); Dombey and Calogeracos, Phys. Rep. 315, 41 (1999)).

Parameters:

E (ndarray) – Energy (including rest mass), \(E > m\).
V0 (ndarray) – Step height.
m (float) – Mass.
group_velocity (bool) – If False, use the naive \(q > 0\) (phase-velocity) choice, which in the Klein zone gives \(R > 1\), \(T < 0\).

Return type:

tuple

Returns:

R, T (numpy.ndarray) – Reflection and transmission probabilities (particle-current ratios).

Examples

Far inside the Klein zone \(T \to 4\kappa_\infty/(1+\kappa_\infty)^2\) with \(\kappa_\infty = \sqrt{(E+m)/(E-m)}\):

>>> import numpy as np
>>> R, T = dirac_step_scattering(2.0, 1e8)
>>> kinf = np.sqrt(3.0)
>>> bool(np.isclose(T, 4 * kinf / (1 + kinf) ** 2))
True

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:

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.electron_repulsion_tensor(basis)[source]#

Two-electron integrals \((\mu\nu|\lambda\sigma)\) in chemists’ notation.

Per primitive quartet

\[\frac{2\pi^{5/2}}{pq\sqrt{p+q}}\,K_{AB}K_{CD}\, F_0\!\left(\frac{pq}{p+q}|P - Q|^2\right)\]

(Szabo and Ostlund, Eq. A.41), filled using the 8-fold permutational symmetry.

Parameters:: basis (Sequence[GaussianS])
Return type:: ndarray
Returns:: numpy.ndarray – Shape (n, n, n, n).

physicskit.quantum.even_tempered_s_basis(center, n, alpha_min, ratio)[source]#

Uncontracted even-tempered s basis, \(\alpha_k = \alpha_{\min}\,\beta^{k}\), \(k = 0,\dots,n-1\).

Even-tempered sequences (Ruedenberg, Raffenetti and Bardo, 1973) approach the complete s basis as \(n\) grows, and with it the Hartree-Fock limit of a two-electron atom.

Parameters:

center (Sequence[float]) – Center (bohr).
n (int) – Number of functions.
alpha_min (float) – Smallest exponent.
ratio (float) – Geometric ratio \(\beta > 1\).

Return type:

list[GaussianS]

Returns:

list of GaussianS

Examples

>>> len(even_tempered_s_basis([0, 0, 0], 8, 0.1, 3.0))
8

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.fine_structure_expansion(n, j, Z=1.0, alpha=0.0072973525643)[source]#

Binding energy through order \((Z\alpha)^4\), in units of \(mc^2\).

\[\frac{E_{nj} - mc^2}{mc^2} \approx -\frac{(Z\alpha)^2}{2n^2} - \frac{(Z\alpha)^4}{2n^4}\left(\frac{n}{j + 1/2} - \frac34\right)\]

– the Bohr level plus the fine-structure correction (relativistic kinetic energy, spin-orbit, and Darwin terms combined; Sommerfeld 1916, Bethe and Salpeter, Quantum Mechanics of One- and Two-Electron Atoms, 1957, Eq. 17.13).

Parameters:

n (int) – Principal quantum number.
j (float) – Total angular momentum.
Z (float) – Nuclear charge.
alpha (float) – Fine-structure constant.

Return type:

Returns:

float

Examples

The \(2P_{3/2}\)-\(2P_{1/2}\) splitting is \(\alpha^4/32\):

>>> import numpy as np
>>> a = 7.2973525693e-3
>>> split = fine_structure_expansion(2, 1.5, alpha=a) - fine_structure_expansion(2, 0.5, alpha=a)
>>> bool(np.isclose(split, a**4 / 32))
True

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:

Returns:

numpy.ndarray

physicskit.quantum.gamma_matrices()[source]#

Dirac gamma matrices \(\gamma^\mu\) in the Dirac (standard) representation.

\[\begin{split}\gamma^0 = \begin{pmatrix} I & 0 \\ 0 & -I \end{pmatrix},\qquad \gamma^i = \begin{pmatrix} 0 & \sigma_i \\ -\sigma_i & 0 \end{pmatrix},\end{split}\]

satisfying \(\{\gamma^\mu, \gamma^\nu\} = 2 g^{\mu\nu}\) with \(g = \mathrm{diag}(1,-1,-1,-1)\) (Bjorken and Drell, Relativistic Quantum Mechanics, 1964, Eq. 2.7).

Return type:: ndarray
Returns:: numpy.ndarray – Shape (4, 4, 4): gamma[mu] is \(\gamma^\mu\).

Examples

>>> import numpy as np
>>> g = gamma_matrices()
>>> bool(np.allclose(g[1] @ g[1], -np.eye(4)))
True

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.grover_circuit(n, marked, iterations=None)[source]#

Grover search for the marked basis states of n qubits.

Starting from the uniform superposition, each iteration applies the phase oracle \(O|x\rangle = -|x\rangle\) for marked \(x\) and the diffusion operator \(2|s\rangle\langle s| - I = H^{\otimes n}(2|0\rangle\langle 0| - I)H^{\otimes n}\), rotating the state by \(2\theta\) toward the marked subspace (Grover, Phys. Rev. Lett. 79, 325 (1997)).

Parameters:

n (int)
marked (ArrayLike) – Indices of the marked items.
iterations (int | None) – Defaults to grover_optimal_iterations().

Return type:

Returns:

QuantumCircuit

Examples

>>> psi = grover_circuit(5, [19]).run()
>>> bool(abs(psi[19]) ** 2 > 0.99)
True

physicskit.quantum.grover_optimal_iterations(n_items, n_marked=1)[source]#

\(\lfloor \pi/(4\theta) \rfloor\) with \(\sin\theta = \sqrt{M/N}\).

Return type:

int

Parameters:

n_items (int)
n_marked (int)

physicskit.quantum.grover_success_probability(k, n_items, n_marked=1)[source]#

Probability of measuring a marked item after k iterations, \(\sin^2((2k+1)\theta)\).

Return type:

NDArray[double]

Parameters:

k (ArrayLike)
n_items (int)
n_marked (int)

Examples

>>> round(float(grover_success_probability(1, 4)), 12)  # N = 4: one query suffices
1.0

physicskit.quantum.hard_sphere_phase_shifts(k, a, l_max)[source]#

Exact phase shifts of an impenetrable sphere of radius \(a\): \(\tan\delta_\ell = j_\ell(ka)/y_\ell(ka)\).

The radial wavefunction \(\propto \cos\delta_\ell\,j_\ell(kr) - \sin\delta_\ell\,y_\ell(kr)\) must vanish at \(r = a\) (Sakurai, Eq. 7.6.31); in particular \(\delta_0 = -ka\) (modulo \(\pi\)). The total cross section tends to \(4\pi a^2\) as \(ka \to 0\) and to \(2\pi a^2\) as \(ka \to \infty\) (the extra \(\pi a^2\) is the diffraction “shadow”), approached as \(2\pi a^2[1 + 0.9962\,(ka)^{-2/3}]\) (Nussenzveig, J. Math. Phys. 10, 82 (1969)).

Parameters:

k (float) – Wavenumber.
a (float) – Sphere radius.
l_max (int) – Highest partial wave.

Return type:

Returns:

numpy.ndarray – Shape (l_max + 1,), in \((-\pi/2, \pi/2]\).

Examples

>>> import numpy as np
>>> bool(np.isclose(hard_sphere_phase_shifts(0.4, 1.0, 0)[0], -0.4))
True

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:

A (ndarray) – The matrix to test.
atol (float) – Absolute tolerance passed to numpy.allclose().

Return type:

bool

Returns:

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

physicskit.quantum.is_trace_preserving(kraus, atol=1e-12)[source]#

Check the completeness relation \(\sum_i E_i^\dagger E_i = I\).

Parameters:

kraus (Sequence[ndarray]) – Kraus operators.
atol (float) – Absolute tolerance.

Return type:

bool

Returns:

bool

Examples

>>> is_trace_preserving(amplitude_damping_kraus(0.4))
True

physicskit.quantum.kinetic_matrix(basis)[source]#

Kinetic-energy matrix \(T_{\mu\nu} = \langle\phi_\mu|-\tfrac12\nabla^2|\phi_\nu\rangle\).

Per primitive pair \(\frac{ab}{p}\left(3 - \frac{2ab}{p}|A-B|^2\right)(\pi/p)^{3/2}e^{-ab|A-B|^2/p}\) (Szabo and Ostlund, Eq. A.11).

Parameters:: basis (Sequence[GaussianS])
Return type:: ndarray
Returns:: numpy.ndarray

physicskit.quantum.klein_gordon_dispersion(p, m=1.0)[source]#

Positive-energy branch \(E = \sqrt{p^2 + m^2}\) of the Klein-Gordon equation.

Substituting \(\phi = e^{i(\mathbf p\cdot\mathbf x - Et)}\) into \((\partial_t^2 - \nabla^2 + m^2)\phi = 0\) (Klein, Z. Phys. 37, 895 (1926); Gordon, Z. Phys. 40, 117 (1926)) gives \(E^2 = p^2 + m^2\), with solutions of both signs \(E = \pm\sqrt{p^2+m^2}\).

Parameters:

p (ndarray | float) – Momentum magnitude \(|\mathbf p|\).
m (float) – Mass.

Return type:

Returns:

numpy.ndarray

Examples

>>> float(klein_gordon_dispersion(0.75, m=1.0))
1.25

physicskit.quantum.klein_gordon_plane_wave(x, t, p, m=1.0, energy_sign=1)[source]#

Free 1D Klein-Gordon plane wave \(\phi(x,t) = e^{i(px - Et)}\), \(E = \pm\sqrt{p^2+m^2}\).

Parameters:

x (ndarray) – Positions.
t (float) – Time.
p (float) – Momentum.
m (float) – Mass.
energy_sign (int) – Sign of the energy.

Return type:

Returns:

numpy.ndarray – Complex amplitudes.

Examples

>>> import numpy as np
>>> phi = klein_gordon_plane_wave(np.array([0.0]), t=0.0, p=1.0)
>>> phi
array([1.+0.j])

physicskit.quantum.klein_gordon_step_scattering(E, V0, m=1.0, group_velocity=True)[source]#

Reflection and transmission of a Klein-Gordon (spin-0) particle at a potential step.

Continuity of \(\phi\) and \(\phi'\) at \(z = 0\) gives

\[R = \left|\frac{p - q}{p + q}\right|^2, \qquad T = \frac{4p\,\mathrm{Re}\,q}{|p + q|^2},\]

with the transmitted wave chosen to move away from the step (group velocity \(q/(E - V_0) > 0\)). In the Klein zone \(V_0 > E + m\) this makes \(q < 0\): the transmitted charge current is negative and \(R > 1\) – the bosonic “superradiant” Klein paradox, the counterpart of the fermionic \(R < 1\) of dirac_step_scattering() (Greiner, Sec. 1.13; Dombey and Calogeracos 1999). At \(V_0 = 2E\) the transmitted and incident momenta cancel, \(q = -p\), and \(R\) and \(T\) diverge.

Parameters:

E (ndarray) – Energy (including rest mass), \(E > m\).
V0 (ndarray) – Step height.
m (float) – Mass.
group_velocity (bool) – If False, use the naive \(q > 0\) choice.

Return type:

tuple

Returns:

R, T (numpy.ndarray) – Reflection and transmission coefficients (charge-current ratios), \(R + T = 1\).

Examples

>>> R, T = klein_gordon_step_scattering(2.0, 5.0)
>>> bool(R > 1 and T < 0 and abs(R + T - 1) < 1e-12)
True

physicskit.quantum.lindblad_rhs(rho, H, c_ops=None)[source]#

Right-hand side of the Lindblad master equation, \(d\rho/dt\).

\[\mathcal{L}\rho = -i[H, \rho] + \sum_k \left( L_k \rho L_k^\dagger - \tfrac12 L_k^\dagger L_k \rho - \tfrac12 \rho L_k^\dagger L_k \right)\]

(Lindblad 1976, Eq. 5.1, with \(\hbar = 1\)).

Parameters:

rho (ndarray) – Density matrix, shape (d, d).
H (ndarray) – Hamiltonian, shape (d, d).
c_ops (Sequence[ndarray] | None) – Collapse (jump) operators \(L_k\), rates absorbed into them.

Return type:

Returns:

numpy.ndarray – \(d\rho/dt\), shape (d, d).

Examples

Spontaneous decay drains the excited state at rate \(\gamma\):

>>> import numpy as np
>>> rho_e = np.diag([0.0, 1.0]).astype(complex)
>>> drho = lindblad_rhs(rho_e, np.zeros((2, 2)), [np.sqrt(0.5) * sigma_minus()])
>>> np.real(np.diag(drho))
array([ 0.5, -0.5])

physicskit.quantum.lindblad_steady_state(H, c_ops)[source]#

The stationary state \(\mathcal{L}\rho_{ss} = 0\), normalized to unit trace.

Solved as the null space of lindblad_superoperator(). For a unique steady state the null space is one-dimensional.

Parameters:

H (ndarray) – Hamiltonian, shape (d, d).
c_ops (Sequence[ndarray] | None) – Collapse operators.

Return type:

Returns:

numpy.ndarray – \(\rho_{ss}\), shape (d, d).

Raises:

ValueError – If the steady state is not unique.

Examples

Resonance fluorescence: a qubit driven at Rabi frequency \(\Omega\) and decaying at rate \(\gamma\) settles to excited population \(\Omega^2/(\gamma^2 + 2\Omega^2)\) (Loudon, The Quantum Theory of Light, 3rd ed., Eq. 2.8.21 on resonance):

>>> import numpy as np
>>> Omega, gamma = 1.0, 0.5
>>> H = 0.5 * Omega * np.array([[0, 1], [1, 0]])
>>> rho = lindblad_steady_state(H, [np.sqrt(gamma) * sigma_minus()])
>>> round(float(rho[1, 1].real), 10) == round(Omega**2 / (gamma**2 + 2 * Omega**2), 10)
True

physicskit.quantum.lindblad_superoperator(H, c_ops=None)[source]#

The Liouvillian \(\mathcal{L}\) as a \(d^2 \times d^2\) matrix acting on row-major \(\mathrm{vec}(\rho)\).

Uses \(\mathrm{vec}(A\rho B) = (A \otimes B^{T})\,\mathrm{vec}(\rho)\) for the row-major (C-order) flattening rho.reshape(-1), so

\[\mathcal{L} = -i\,(H \otimes I - I \otimes H^{T}) + \sum_k \left( L_k \otimes L_k^{*} - \tfrac12 L_k^\dagger L_k \otimes I - \tfrac12 I \otimes (L_k^\dagger L_k)^{T} \right).\]

Parameters:

H (ndarray) – Hamiltonian, shape (d, d).
c_ops (Sequence[ndarray] | None) – Collapse operators.

Return type:

Returns:

numpy.ndarray – Shape (d**2, d**2), complex.

Examples

Its action on rho.reshape(-1) agrees with lindblad_rhs():

>>> import numpy as np
>>> H = np.array([[0.0, 0.3], [0.3, 1.0]])
>>> ops = [sigma_minus()]
>>> rho = np.array([[0.4, 0.1j], [-0.1j, 0.6]])
>>> S = lindblad_superoperator(H, ops)
>>> np.allclose(S @ rho.reshape(-1), lindblad_rhs(rho, H, ops).reshape(-1))
True

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:

Returns:

float

Raises:

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

physicskit.quantum.modular_exponentiation_circuit(a, N, n_count)[source]#

Prepare \(2^{-t/2}\sum_x |x\rangle|a^x \bmod N\rangle\) for Shor’s period finding.

The counting register (n_count qubits) is put in uniform superposition, and the reversible map \(|x\rangle|y\rangle \to |x\rangle|y \oplus (a^x \bmod N)\rangle\) is applied as one permutation gate, which stands in for the arithmetic circuit a real device would compile.

Parameters:

a (int) – Base, coprime to N.
N (int) – Number to factor.
n_count (int) – Counting qubits \(t\), typically with \(N^2 \le 2^t\).

Return type:

tuple[QuantumCircuit, int]

Returns:

circuit (QuantumCircuit) – On n_count + m qubits, counting register first.
m (int) – Work-register size, \(\lceil\log_2 N\rceil\).

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:

Returns:

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

physicskit.quantum.momentum_transfer(k, theta)[source]#

Momentum transfer \(q = |\mathbf k' - \mathbf k| = 2k\sin(\theta/2)\) for elastic scattering.

Parameters:

k (float) – Wavenumber.
theta (ndarray) – Scattering angle (radians).

Return type:

Returns:

numpy.ndarray

Examples

>>> import numpy as np
>>> float(momentum_transfer(1.5, np.pi))
3.0

physicskit.quantum.nuclear_attraction_matrix(basis, charges, positions)[source]#

Electron-nuclear attraction \(V_{\mu\nu} = -\sum_C Z_C\langle\phi_\mu|1/|\mathbf r - \mathbf R_C||\phi_\nu\rangle\).

Per primitive pair \(-\frac{2\pi}{p}Z_C\,e^{-ab|A-B|^2/p}F_0(p|P - C|^2)\) (Szabo and Ostlund, Eq. A.33).

Parameters:

basis (Sequence[GaussianS])
charges (Sequence[float] | ndarray) – Nuclear charges \(Z_C\).
positions (ndarray) – Nuclear positions, shape (n_nuclei, 3).

Return type:

Returns:

numpy.ndarray

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:

orb_a (HydrogenOrbital)
orb_b (HydrogenOrbital)
r (ndarray)
theta (ndarray)
phi (ndarray)
t (float)
ca (complex)
cb (complex)

Return type:

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:

Returns:

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

physicskit.quantum.overlap_matrix(basis)[source]#

Overlap matrix \(S_{\mu\nu} = \langle\phi_\mu|\phi_\nu\rangle\), with \((\pi/p)^{3/2}e^{-ab|A-B|^2/p}\) per primitive pair.

Parameters:: basis (Sequence[GaussianS])
Return type:: ndarray
Returns:: numpy.ndarray

physicskit.quantum.partial_wave_amplitude(theta, k, deltas)[source]#

Scattering amplitude from phase shifts, \(f(\theta) = \frac1k\sum_\ell (2\ell+1)e^{i\delta_\ell}\sin\delta_\ell\,P_\ell(\cos\theta)\).

(Faxen and Holtsmark, Z. Phys. 45, 307 (1927); Sakurai, Eq. 7.6.18.)

Parameters:

theta (ndarray) – Scattering angles (radians).
k (float) – Wavenumber.
deltas (ndarray) – Phase shifts \(\delta_0, \delta_1, \dots\).

Return type:

Returns:

numpy.ndarray – Complex amplitude.

Examples

Pure s-wave scattering is isotropic:

>>> import numpy as np
>>> f = partial_wave_amplitude(np.array([0.0, np.pi]), 1.0, [np.pi / 2])
>>> bool(np.allclose(f, [1j, 1j]))
True

physicskit.quantum.partial_wave_cross_section(k, deltas)[source]#

Total cross section \(\sigma = \frac{4\pi}{k^2}\sum_\ell (2\ell+1)\sin^2\delta_\ell\).

Equal, by the optical theorem, to \(\frac{4\pi}{k}\,\mathrm{Im}\,f(0)\) (Sakurai, Eq. 7.6.20).

Parameters:

k (float) – Wavenumber.
deltas (ndarray) – Phase shifts.

Return type:

Returns:

float

Examples

The unitarity limit of a resonant s wave is \(4\pi/k^2\):

>>> import numpy as np
>>> bool(np.isclose(partial_wave_cross_section(2.0, [np.pi / 2]), np.pi))
True

physicskit.quantum.partial_wave_phase_shifts(potential, k, l_max, mass=1.0, r_max=30.0, r_min=1e-08, breakpoints=(), rtol=1e-10)[source]#

Phase shifts \(\delta_\ell(k)\), \(\ell = 0,\dots,\ell_{\max}\), of a short-range central potential.

Integrates Calogero’s variable-phase equation

\[\frac{d\delta_\ell}{dr} = -\frac{U(r)}{k}\left[\hat\jmath_\ell(kr)\cos\delta_\ell(r) - \hat n_\ell(kr)\sin\delta_\ell(r)\right]^2, \qquad U = 2mV,\]

from \(\delta_\ell(0) = 0\) out to r_max, with the Riccati-Bessel functions \(\hat\jmath_\ell(x) = x j_\ell(x)\) and \(\hat n_\ell(x) = x y_\ell(x)\). \(\delta_\ell(r)\) is the phase shift produced by the potential truncated at \(r\), so its limit is the full phase shift (F. Calogero, Variable Phase Approach to Potential Scattering, Academic Press, 1967, Eq. 3.4). Unlike shooting and matching, it never overflows at large \(\ell\) and gives a continuous \(\delta_\ell\) (not reduced modulo \(\pi\)).

Parameters:

potential (Callable[[ndarray], ndarray]) – \(V(r)\), finite for \(r > 0\) and negligible beyond r_max.
k (float) – Wavenumber, \(k > 0\).
l_max (int) – Highest partial wave.
mass (float) – Particle (reduced) mass.
r_max (float) – Outer integration radius.
r_min (float) – Inner starting radius.
breakpoints (Sequence[float]) – Radii where \(V\) is discontinuous (e.g. a square well’s edge). The integration restarts there instead of shrinking its step size to resolve the jump.
rtol (float) – Relative tolerance of the ODE solver (absolute tolerance is rtol / 100).

Return type:

Returns:

numpy.ndarray – Shape (l_max + 1,), in radians.

Examples

An attractive square well of depth \(V_0\) and radius \(a\) has \(\tan(ka + \delta_0) = (k/K)\tan(Ka)\), \(K = \sqrt{k^2 + 2mV_0}\):

>>> import numpy as np
>>> well = lambda r: np.where(r < 1.0, -0.8, 0.0)
>>> d0 = partial_wave_phase_shifts(well, k=0.7, l_max=0, r_max=1.5, breakpoints=[1.0])[0]
>>> K = np.sqrt(0.49 + 1.6)
>>> bool(np.isclose(np.tan(0.7 + d0), 0.7 / K * np.tan(K)))
True

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:

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:

x (ndarray) – Uniform spatial grids.
y (ndarray) – Uniform spatial grids.
x0 (float) – Initial center of the packet along \(x\) (should sit well upstream of wall_x).
k0 (float) – Central wavenumber along \(x\).
sigma_x (float) – Initial packet widths (sigma_y should span both slits).
sigma_y (float) – Initial packet widths (sigma_y should span both slits).
wall_x (float) – Passed to double_slit_potential().
wall_thickness (float) – Passed to double_slit_potential().
slit_separation (float) – Passed to double_slit_potential().
slit_width (float) – Passed to double_slit_potential().
wall_height (float) – Passed to double_slit_potential().
hbar (float) – Propagation parameters, passed to SplitOperatorSolver2D.
m (float) – Propagation parameters, passed to SplitOperatorSolver2D.
dt (float) – Propagation parameters, passed to SplitOperatorSolver2D.
n_steps (int) – Total number of time steps.
save_every (int) – Save a snapshot every this many steps.

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.qft_circuit(n, swaps=True)[source]#

Quantum Fourier transform on n qubits from Hadamards and controlled phases.

Qubit \(j\) gets a Hadamard followed by controlled phases \(2\pi/2^{k-j+1}\) from each later qubit \(k\); the final swaps reverse the qubit order (Coppersmith, IBM Research Report RC19642 (1994); Nielsen and Chuang, sec. 5.1).

Parameters:

n (int)
swaps (bool)

Return type:

Returns:

QuantumCircuit

Examples

>>> N = 8
>>> F = np.exp(2j * np.pi * np.outer(np.arange(N), np.arange(N)) / N) / np.sqrt(N)
>>> bool(np.allclose(qft_circuit(3).unitary(), F))
True

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:

Returns:

numpy.ndarray

Raises:

ValueError – If l >= n.

physicskit.quantum.restricted_hartree_fock(basis, charges, positions, n_electrons, max_iter=200, tol=1e-10, damping=0.0)[source]#

Closed-shell Hartree-Fock by the Roothaan self-consistent-field procedure.

With core Hamiltonian \(H^\mathrm{core} = T + V\), the Fock matrix

\[F_{\mu\nu} = H^\mathrm{core}_{\mu\nu} + \sum_{\lambda\sigma}P_{\lambda\sigma} \left[(\mu\nu|\sigma\lambda) - \tfrac12(\mu\lambda|\sigma\nu)\right]\]

is diagonalized in the symmetrically orthogonalized basis \(X = S^{-1/2}\), the \(N/2\) lowest orbitals are doubly occupied, and the density \(P = 2C_\mathrm{occ}C_\mathrm{occ}^T\) is iterated to self-consistency, starting from the core-Hamiltonian guess (Szabo and Ostlund, Sec. 3.4.6, Eqs. 3.154 and 3.184). By the variational principle the energy is an upper bound, which approaches the Hartree-Fock limit from above as the basis grows.

Parameters:

basis (Sequence[GaussianS]) – Basis functions.
charges (Sequence[float] | ndarray) – Nuclear charges.
positions (ndarray) – Nuclear positions (bohr), shape (n_nuclei, 3).
n_electrons (int) – Even number of electrons.
max_iter (int) – Maximum SCF iterations.
tol (float) – Convergence threshold on the density (RMS change) and energy.
damping (float) – Fraction of the previous density mixed into each update.

Return type:

HartreeFockResult

Returns:

HartreeFockResult

Raises:

ValueError – For an odd electron count or too small a basis.

Examples

\(\mathrm H_2\) at \(R = 1.4\) bohr in STO-3G (Szabo and Ostlund, Sec. 3.5.2: \(E = -1.1167\) hartree):

>>> import numpy as np
>>> R = np.array([[0.0, 0.0, 0.0], [1.4, 0.0, 0.0]])
>>> basis = [sto3g_1s(1.24, r) for r in R]
>>> res = restricted_hartree_fock(basis, [1, 1], R, n_electrons=2)
>>> round(res.energy, 4)
-1.1167

physicskit.quantum.shor_period_from_measurement(y, n_count, N)[source]#

Period candidate from a counting-register outcome, by continued fractions.

\(y/2^t \approx s/r\); the denominator of the best approximation with \(r < N\) is returned.

Return type:

int

Parameters:

y (int)
n_count (int)
N (int)

Examples

>>> shor_period_from_measurement(192, 8, 15)  # 192/256 = 3/4
4

physicskit.quantum.sigma_minus()[source]#

Qubit lowering operator \(\sigma_- = |0\rangle\langle 1|\).

Return type:: ndarray
Returns:: numpy.ndarray – Shape (2, 2), complex.

Examples

>>> import numpy as np
>>> sigma_minus() @ np.array([0, 1])  # |1> -> |0>
array([1.+0.j, 0.+0.j])

physicskit.quantum.solve_lindblad(rho0, H, c_ops, times)[source]#

Propagate a density matrix under a time-independent Lindblad equation.

The solution is exact up to round-off: \(\mathrm{vec}\,\rho(t) = e^{\mathcal{L}(t - t_0)}\,\mathrm{vec}\,\rho(t_0)\), with the matrix exponential of lindblad_superoperator() taken once per distinct time step (a uniform grid costs a single scipy.linalg.expm()).

Parameters:

rho0 (ndarray) – Initial density matrix at times[0], shape (d, d). A state vector of shape (d,) is accepted and converted to \(|\psi\rangle\langle\psi|\).
H (ndarray) – Hamiltonian, shape (d, d).
c_ops (Sequence[ndarray] | None) – Collapse operators.
times (ndarray) – Increasing output times.

Return type:

Returns:

numpy.ndarray – \(\rho(t)\) at each time, shape (len(times), d, d).

Examples

\(T_1\) decay of the excited state, \(P_1(t) = e^{-t/T_1}\):

>>> import numpy as np
>>> t = np.array([0.0, 1.0, 2.0])
>>> rho = solve_lindblad(np.array([0, 1]), np.zeros((2, 2)), [sigma_minus()], t)
>>> np.allclose(rho[:, 1, 1].real, np.exp(-t))
True

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:

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:

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.sto3g_1s(zeta, center)[source]#

STO-3G contraction of a Slater 1s orbital with exponent \(\zeta\).

The \(\zeta = 1\) exponents are scaled as \(\alpha = \zeta^2\alpha_{(1)}\) (Szabo and Ostlund, Eq. 3.224).

Parameters:

zeta (float) – Slater exponent (see STO3G_ZETA).
center (Sequence[float]) – Position (bohr).

Return type:

GaussianS

Returns:

GaussianS

Examples

>>> import numpy as np
>>> phi = sto3g_1s(1.24, [0.0, 0.0, 0.0])
>>> round(float(overlap_matrix([phi])[0, 0]), 5)
1.0

physicskit.quantum.t1_t2_collapse_operators(T1, T2)[source]#

Collapse operators giving a qubit energy-relaxation time \(T_1\) and coherence time \(T_2\).

Relaxation \(L_1 = \sqrt{1/T_1}\,\sigma_-\) damps the excited population as \(e^{-t/T_1}\) and, on its own, the coherences as \(e^{-t/(2T_1)}\). The remaining decay comes from pure dephasing \(L_\varphi = \sqrt{\gamma_\varphi/2}\,\sigma_z\) with

\[\gamma_\varphi = \frac{1}{T_2} - \frac{1}{2T_1},\]

so that \(|\rho_{01}(t)| = |\rho_{01}(0)|\,e^{-t/T_2}\) (Breuer and Petruccione, The Theory of Open Quantum Systems, 2002, Sec. 3.4). Physical rates need \(T_2 \le 2T_1\).

Parameters:

T1 (float) – Energy-relaxation time, \(T_1 > 0\).
T2 (float) – Coherence time, \(0 < T_2 \le 2T_1\).

Return type:

list[ndarray]

Returns:

list of numpy.ndarray – [L_1] if \(T_2 = 2T_1\) (no pure dephasing), otherwise [L_1, L_phi].

Raises:

ValueError – If \(T_2 > 2T_1\) or either time is not positive.

Examples

>>> ops = t1_t2_collapse_operators(T1=1.0, T2=2.0)
>>> len(ops)
1

physicskit.quantum.tebd(psi, H_bonds, dt, n_steps, chi_max=64, imaginary=False, eps=1e-12, measure_every=1, observables=None)[source]#

Evolve an MPS with second-order TEBD.

Parameters:

psi (MPS) – Modified in place.
H_bonds (list[NDArray]) – Two-site Hamiltonian terms, e.g. from xxz_bond_hamiltonians().
dt (float) – Time step (imaginary-time step if imaginary).
n_steps (int)
chi_max (int) – Maximum bond dimension.
imaginary (bool) – Imaginary-time evolution \(e^{-H\tau}\) toward the ground state.
eps (float) – Relative singular-value cutoff.
measure_every (int)
observables (dict | None) – Single-site operators whose site-resolved expectation values are recorded.

Return type:

dict[str, NDArray]

Returns:

dict – "t", the total "energy", the half-chain "entropy", the accumulated "discarded_weight", and one array per observable.

Examples

Imaginary-time TEBD finds the ground state of an open 8-site Heisenberg chain:

>>> from physicskit.condensed.spin_chains import lowest_eigenstates, xxz_hamiltonian
>>> H = xxz_bond_hamiltonians(8)
>>> psi = MPS.product_state([0, 1] * 4)
>>> out = tebd(psi, H, dt=0.05, n_steps=400, chi_max=16, imaginary=True)
>>> E0 = lowest_eigenstates(xxz_hamiltonian(8, periodic=False, n_up=4)[0], k=1)[0][0]
>>> bool(abs(out["energy"][-1] - E0) < 1e-6)
True

physicskit.quantum.tfim_bond_hamiltonians(N, J=1.0, h=1.0)[source]#

Two-site terms of the open transverse-field Ising chain, \(H = -J\sum\sigma^x\sigma^x - h\sum\sigma^z\).

Matches physicskit.condensed.spin_chains.tfim_hamiltonian() with periodic=False.

Parameters:

N (int)
J (float)
h (float)

Return type:

list[NDArray[cdouble]]

Returns:

list of ndarray, length N - 1

physicskit.quantum.xxz_bond_hamiltonians(N, J=1.0, delta=1.0, h=0.0)[source]#

Two-site terms of the open XXZ chain, \(H = J\sum(S^xS^x + S^yS^y + \Delta S^zS^z) - h\sum S^z\).

Each bond operator is a \(4\times4\) matrix on (site \(j\), site \(j+1\)); the field of each site is split between its bonds so that \(\sum_j h_{j,j+1}\) is exactly the Hamiltonian of physicskit.condensed.spin_chains.xxz_hamiltonian() with periodic=False.

Parameters:

N (int)
J (float)
delta (float)
h (float)

Return type:

list[NDArray[cdouble]]

Returns:

list of ndarray, length N - 1

physicskit.quantum.yukawa_born_amplitude(q, g, mu, mass=1.0)[source]#

Closed-form Born amplitude of the Yukawa potential \(V(r) = g\,e^{-\mu r}/r\).

\[f_B(q) = -\frac{2mg}{\mu^2 + q^2}\]

(Sakurai, Eq. 7.2.15). As \(\mu \to 0\) (unscreened Coulomb, \(g = Z_1Z_2\alpha\)) \(|f_B|^2\) becomes Rutherford’s classical cross section,

\[\frac{d\sigma}{d\Omega} = \frac{4m^2g^2}{q^4} = \left(\frac{g}{4E}\right)^2 \frac{1}{\sin^4(\theta/2)}, \qquad E = \frac{k^2}{2m}.\]

Parameters:

q (ndarray) – Momentum transfer.
g (float) – Coupling strength (\(g > 0\) repulsive).
mu (float) – Inverse screening length.
mass (float) – Particle (reduced) mass.

Return type:

Returns:

numpy.ndarray

Examples

>>> float(yukawa_born_amplitude(1.0, g=0.5, mu=1.0))
-0.5

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:

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: