Theoretical Background¶

This document provides a comprehensive overview of the mathematical and physical foundations underlying the schr package. It is intended for researchers seeking a deeper understanding of the numerical methods and quantum mechanics principles implemented in the framework.

Quantum Mechanics Foundations ¶

Hilbert Space Formalism ¶

The state of a quantum system is represented by a vector \(|\psi\rangle\) in a complex Hilbert space \(\mathcal{H}\). For a single particle in one dimension, this is typically the space \(L^2(\mathbb{R})\) of square-integrable functions.

Inner Product:

\[\langle\psi|\phi\rangle = \int_{-\infty}^{\infty} \psi^*(x) \phi(x) \, dx\]

Normalization Condition:

\[\langle\psi|\psi\rangle = \int_{-\infty}^{\infty} |\psi(x)|^2 \, dx = 1\]

The quantity \(|\psi(x)|^2\) represents the probability density of finding the particle at position \(x\).

Operators and Observables ¶

Physical observables correspond to Hermitian (self-adjoint) operators. Key operators include:

Position Operator:

\[\hat{x} \psi(x) = x \psi(x)\]

Momentum Operator:

\[\hat{p} \psi(x) = -i\hbar \frac{\partial\psi}{\partial x}\]

Hamiltonian (Energy) Operator:

\[\hat{H} = \frac{\hat{p}^2}{2m} + V(\hat{x}, t) = -\frac{\hbar^2}{2m} \frac{\partial^2}{\partial x^2} + V(x, t)\]

Expectation Values:

The expected value of an observable \(\hat{A}\) in state \(|\psi\rangle\):

\[\langle\hat{A}\rangle = \langle\psi|\hat{A}|\psi\rangle = \int_{-\infty}^{\infty} \psi^*(x) \hat{A} \psi(x) \, dx\]

Time Evolution ¶

Time-Dependent Schrödinger Equation (TDSE)¶

The evolution of quantum states is governed by the TDSE:

\[i\hbar \frac{\partial|\psi(t)\rangle}{\partial t} = \hat{H}|\psi(t)\rangle\]

In position representation:

\[i\hbar \frac{\partial\psi(x,t)}{\partial t} = \left[-\frac{\hbar^2}{2m} \frac{\partial^2}{\partial x^2} + V(x,t)\right] \psi(x,t)\]

Formal Solution:

For time-independent Hamiltonians, the formal solution is:

\[|\psi(t)\rangle = e^{-i\hat{H}t/\hbar}|\psi(0)\rangle = \hat{U}(t)|\psi(0)\rangle\]

where \(\hat{U}(t) = e^{-i\hat{H}t/\hbar}\) is the unitary time evolution operator.

Conservation Laws ¶

Probability Conservation:

\[\frac{d}{dt}\int_{-\infty}^{\infty}|\psi(x,t)|^2\,dx = 0\]

Energy Conservation (time-independent H):

\[\frac{d}{dt}\langle\hat{H}\rangle = 0\]

Numerical Methods ¶

Split-Step Fourier Method ¶

The split-step Fourier method (SSFM) is a pseudo-spectral technique for solving the TDSE. It exploits the fact that kinetic and potential operators are simple in different representations.

Operator Splitting ¶

The Hamiltonian is split into kinetic and potential parts:

\[\hat{H} = \hat{T} + \hat{V}\]

where

\[\hat{T} = -\frac{\hbar^2}{2m}\nabla^2, \quad \hat{V} = V(\mathbf{r}, t)\]

The evolution operator is approximated using the Strang splitting formula:

\[e^{-i\hat{H}\Delta t/\hbar} \approx e^{-i\hat{V}\Delta t/(2\hbar)} \, e^{-i\hat{T}\Delta t/\hbar} \, e^{-i\hat{V}\Delta t/(2\hbar)} + O(\Delta t^3)\]

This is a symmetric splitting, providing second-order accuracy in time.

Implementation ¶

Step 1: Half potential step (position space)

\[\psi_1(\mathbf{r}) = e^{-iV(\mathbf{r},t)\Delta t/(2\hbar)} \, \psi(\mathbf{r}, t)\]

Step 2: Kinetic step (momentum space)

Transform to momentum space via FFT:

\[\tilde{\psi}_1(\mathbf{k}) = \mathcal{F}[\psi_1(\mathbf{r})]\]

Apply kinetic evolution:

\[\tilde{\psi}_2(\mathbf{k}) = e^{-i\hbar k^2\Delta t/(2m)} \, \tilde{\psi}_1(\mathbf{k})\]

Transform back to position space:

\[\psi_2(\mathbf{r}) = \mathcal{F}^{-1}[\tilde{\psi}_2(\mathbf{k})]\]

Step 3: Half potential step (position space)

\[\psi(\mathbf{r}, t+\Delta t) = e^{-iV(\mathbf{r},t+\Delta t)\Delta t/(2\hbar)} \, \psi_2(\mathbf{r})\]

Advantages ¶

Spectral accuracy in space: FFT provides exponential convergence for smooth functions
Explicit method: No matrix inversion required
Unconditionally stable: Preserves unitarity
Efficient: \(O(N\log N)\) complexity per time step
Automatic boundary conditions: Periodic boundaries from FFT

Stability and Accuracy ¶

CFL Condition:

For stability, the time step should satisfy:

\[\Delta t < \frac{2\pi\hbar m}{k_{\max}^2} = \frac{2\pi\hbar m\Delta x^2}{\pi^2}\]

where \(k_{\max} = \pi/\Delta x\) is the Nyquist wavevector.

Spatial Resolution:

To accurately resolve features with characteristic momentum \(k_0\):

\[\Delta x < \frac{\pi}{k_0}\]

A typical choice is \(\Delta x = \lambda/10\) where \(\lambda = 2\pi/k_0\).

Fourier Transform Conventions ¶

The package uses the following Fourier transform conventions:

Forward Transform:

\[\tilde{\psi}(k) = \mathcal{F}[\psi(x)] = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} \psi(x) e^{-ikx} \, dx\]

Inverse Transform:

\[\psi(x) = \mathcal{F}^{-1}[\tilde{\psi}(k)] = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} \tilde{\psi}(k) e^{ikx} \, dk\]

Discrete Fourier Transform (DFT):

For a grid with \(N\) points and spacing \(\Delta x\):

\[\tilde{\psi}_k = \sum_{n=0}^{N-1} \psi_n e^{-2\pi ikn/N}\]

The momentum grid is:

\[\begin{split}k_j = \begin{cases} \dfrac{2\pi j}{N\Delta x} & 0 \leq j \leq N/2 \\[8pt] \dfrac{2\pi (j-N)}{N\Delta x} & N/2 < j < N \end{cases}\end{split}\]

Parseval’s Theorem (Norm Conservation):

\[\int_{-\infty}^{\infty} |\psi(x)|^2 \, dx = \int_{-\infty}^{\infty} |\tilde{\psi}(k)|^2 \, dk\]

This ensures probability conservation in both representations.

Momentum Operator in Fourier Space ¶

The momentum operator becomes multiplicative in Fourier space:

\[\mathcal{F}[\hat{p}\psi(x)] = \mathcal{F}\left[-i\hbar \frac{\partial\psi}{\partial x}\right] = \hbar k \, \tilde{\psi}(k)\]

Similarly, the kinetic energy operator:

\[\mathcal{F}\left[\frac{\hat{p}^2}{2m}\psi(x)\right] = \frac{\hbar^2 k^2}{2m} \, \tilde{\psi}(k)\]

This is the basis for the efficient SSFM implementation.

Runge-Kutta Methods ¶

The package also implements fourth-order Runge-Kutta (RK4) integration as an alternative solver.

RK4 Algorithm:

To solve \(\partial_t|\psi\rangle = -i\hat{H}|\psi\rangle/\hbar\):

\[\begin{split}\begin{aligned} k_1 &= -\frac{i}{\hbar}\hat{H}|\psi(t)\rangle \\ k_2 &= -\frac{i}{\hbar}\hat{H}\left(|\psi(t)\rangle + \Delta t k_1/2\right) \\ k_3 &= -\frac{i}{\hbar}\hat{H}\left(|\psi(t)\rangle + \Delta t k_2/2\right) \\ k_4 &= -\frac{i}{\hbar}\hat{H}\left(|\psi(t)\rangle + \Delta t k_3\right) \\ |\psi(t+\Delta t)\rangle &= |\psi(t)\rangle + \frac{\Delta t}{6}(k_1 + 2k_2 + 2k_3 + k_4) \end{aligned}\end{split}\]

Advantages: General-purpose, handles non-separable Hamiltonians.

Disadvantages: \(O(N^2)\) for finite-difference derivatives, doesn’t preserve unitarity exactly.

Boundary Conditions ¶

Absorbing Boundary Conditions ¶

For open quantum systems, absorbing boundaries prevent unphysical reflections at domain edges.

Polynomial Mask Method ¶

Multiply the wavefunction by a smooth mask function:

\[\psi(x, t+\Delta t) \rightarrow \mathcal{M}(x)\psi(x, t+\Delta t)\]

where

\[\begin{split}\mathcal{M}(x) = \begin{cases} \cos^n\left(\dfrac{\pi(x - x_0)}{2w}\right) & x_0 < x < x_0 + w \\[8pt] 1 & \text{interior} \\[8pt] \cos^n\left(\dfrac{\pi(x - x_1)}{2w}\right) & x_1 - w < x < x_1 \end{cases}\end{split}\]

Here \(n\) is the order (typically 4-8) and \(w\) is the absorption layer width.

Properties:

Smooth transition minimizes reflections
Higher \(n\) gives sharper absorption
Width \(w\) should be several wavelengths

Complex Absorbing Potential (CAP)¶

Add an imaginary potential in the absorption region:

\[V(x) \rightarrow V(x) - i\alpha f(x)\]

where \(f(x)\) is a smooth function (e.g., polynomial or exponential) near boundaries.

Effect: Exponentially damps the wavefunction amplitude:

\[\psi(x,t) \sim e^{-\alpha t/\hbar}\psi_0(x)\]

Periodic Boundary Conditions ¶

The FFT naturally imposes periodic boundary conditions:

\[\psi(x + L) = \psi(x)\]

where \(L = N\Delta x\) is the domain length. Useful for:

Plane wave states
Bloch states in periodic potentials
Systems with periodic symmetry

Quantum Electrodynamics ¶

Second Quantization ¶

Photon Field Quantization ¶

The electromagnetic field is quantized in terms of photon creation and annihilation operators:

\[\hat{a}_i, \quad \hat{a}_i^\dagger\]

satisfying the commutation relations:

\[[\hat{a}_i, \hat{a}_j^\dagger] = \delta_{ij}, \quad [\hat{a}_i, \hat{a}_j] = 0\]

Number States:

The Fock space is spanned by number states \(|n_1, n_2, \ldots, n_M\rangle\) where \(n_i\) is the number of photons in mode \(i\).

Action of Operators:

\[\begin{split}\begin{aligned} \hat{a}_i^\dagger|n_i\rangle &= \sqrt{n_i + 1}|n_i + 1\rangle \\ \hat{a}_i|n_i\rangle &= \sqrt{n_i}|n_i - 1\rangle \\ \hat{n}_i &= \hat{a}_i^\dagger\hat{a}_i, \quad \hat{n}_i|n_i\rangle = n_i|n_i\rangle \end{aligned}\end{split}\]

Free Field Hamiltonian ¶

\[\hat{H}_{\text{field}} = \sum_{i=1}^{M} \hbar\omega_i\left(\hat{a}_i^\dagger\hat{a}_i + \frac{1}{2}\right)\]

where \(\omega_i\) is the frequency of mode \(i\).

Light-Matter Interaction ¶

Jaynes-Cummings Model ¶

The JC model describes a two-level atom coupled to a single cavity mode:

\[\hat{H}_{\text{JC}} = \frac{\hbar\omega_a}{2}\hat{\sigma}_z + \hbar\omega_c \hat{a}^\dagger\hat{a} + \hbar g(\hat{\sigma}_+\hat{a} + \hat{\sigma}_-\hat{a}^\dagger)\]

where:

\(\omega_a\) is the atomic transition frequency
\(\omega_c\) is the cavity mode frequency
\(g\) is the coupling strength
\(\hat{\sigma}_z = |e\rangle\langle e| - |g\rangle\langle g|\) is the Pauli z operator
\(\hat{\sigma}_\pm = \dfrac{1}{2}(\hat{\sigma}_x \pm i\hat{\sigma}_y)\) are raising/lowering operators

Rotating Wave Approximation (RWA):

The terms \(\hat{\sigma}_+\hat{a}^\dagger\) and \(\hat{\sigma}_-\hat{a}\) are neglected (counter-rotating terms) when \(g \ll \omega_a, \omega_c\).

Physical Phenomena:

Rabi Oscillations: Periodic energy exchange between atom and field
Vacuum Rabi Splitting: Energy level splitting in strong coupling regime
Collapse and Revival: In multi-photon states, non-classical dynamics

Dipole Interaction ¶

For a particle with charge \(q\) interacting with the electromagnetic field:

\[\hat{H}_{\text{int}} = -q\hat{\mathbf{r}} \cdot \hat{\mathbf{E}}(\mathbf{r}, t)\]

where \(\hat{\mathbf{E}}\) is the quantized electric field operator.

In the dipole approximation (wavelength \(\gg\) atomic size):

\[\hat{\mathbf{E}}(\mathbf{r}, t) \approx \hat{\mathbf{E}}(\mathbf{r}_0, t) = i\sum_i\sqrt{\frac{\hbar\omega_i}{2\epsilon_0 V}} \left(\hat{a}_i e^{-i\omega_i t} - \hat{a}_i^\dagger e^{i\omega_i t}\right) \hat{\mathbf{e}}_i\]

where \(\hat{\mathbf{e}}_i\) is the polarization vector.

Fock Space Truncation ¶

In numerical simulations, the infinite-dimensional Fock space is truncated:

\[\mathcal{H} \approx \text{span}\{|n_1, \ldots, n_M\rangle : 0 \leq n_i \leq N_{\max}\}\]

Dimension:

For \(M\) modes with maximum \(N_{\max}\) photons each:

\[\dim(\mathcal{H}) = (N_{\max} + 1)^M\]

Convergence:

Results should be checked for convergence with respect to \(N_{\max}\). Typically \(N_{\max} = 10-50\) is sufficient for weak-to-moderate coupling.

Example Systems ¶

Harmonic Oscillator ¶

The 1D quantum harmonic oscillator has Hamiltonian:

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

Energy Eigenvalues:

\[E_n = \hbar\omega\left(n + \frac{1}{2}\right), \quad n = 0, 1, 2, \ldots\]

Eigenfunctions:

\[\psi_n(x) = \frac{1}{\sqrt{2^n n!}}\left(\frac{m\omega}{\pi\hbar}\right)^{1/4}e^{-m\omega x^2/(2\hbar)}H_n\left(\sqrt{\frac{m\omega}{\hbar}}x\right)\]

where \(H_n\) are Hermite polynomials.

Coherent States:

Displaced Gaussian wavepackets are approximate eigenstates of the annihilation operator:

\[\hat{a}|\alpha\rangle = \alpha|\alpha\rangle\]

with \(|\alpha\rangle = e^{-|\alpha|^2/2}\sum_{n=0}^\infty \frac{\alpha^n}{\sqrt{n!}}|n\rangle\).

Tunneling Through a Barrier ¶

Rectangular Barrier ¶

Potential:

\[\begin{split}V(x) = \begin{cases} V_0 & 0 < x < a \\ 0 & \text{otherwise} \end{cases}\end{split}\]

For \(E < V_0\), the transmission coefficient is:

\[T = \frac{1}{1 + \dfrac{V_0^2\sinh^2(\kappa a)}{4E(V_0 - E)}}\]

where \(\kappa = \sqrt{2m(V_0 - E)}/\hbar\).

Tunneling Time:

The time for a wavepacket to tunnel through is approximately:

\[\tau \sim \frac{ma}{\hbar\kappa}\]

Double-Slit Interference ¶

Physical Parameters ¶

Slit separation: \(d\)
Slit width: \(w\)
Wavelength: \(\lambda = 2\pi/k\)
Screen distance: \(L\)

Fringe Spacing:

For \(w \ll d \ll L\):

\[\Delta y = \frac{\lambda L}{d}\]

Visibility:

The fringe visibility depends on the coherence length \(\ell_c\):

\[V = \frac{I_{\max} - I_{\min}}{I_{\max} + I_{\min}} \approx e^{-d^2/(2\ell_c^2)}\]

For Gaussian wavepackets, \(\ell_c \sim \sigma\) (wavepacket width).

Atomic Units ¶

Fundamental Constants ¶

Constant	Symbol	Value (SI)
Bohr radius	\(a_0\)	\(5.291772 \times 10^{-11}\) m
Hartree energy	\(E_h\)	\(4.359744 \times 10^{-18}\) J = 27.211 eV
Atomic time unit	\(\hbar/E_h\)	\(2.418884 \times 10^{-17}\) s
Electron mass	\(m_e\)	\(9.109384 \times 10^{-31}\) kg
Elementary charge	\(e\)	\(1.602176 \times 10^{-19}\) C

Unit Conversions ¶

Length:

\[1\,a_0 = 0.05292\,\text{nm}\]

Energy:

\[1\,E_h = 27.211\,\text{eV} = 4.360\times 10^{-18}\,\text{J}\]

Time:

\[1\,\hbar/E_h = 24.19\,\text{as} = 2.419\times 10^{-17}\,\text{s}\]

Electric Field:

\[1\,E_h/(ea_0) = 5.142\times 10^{11}\,\text{V/m}\]

Example: Hydrogen Ground State ¶

In atomic units:

Energy: \(E_0 = -0.5\,E_h = -13.6\,\text{eV}\)
Bohr radius: \(a_0 = 1\,a_0\)
Wavefunction: \(\psi_{100}(r) = \frac{1}{\sqrt{\pi}}e^{-r}\)

Numerical Considerations ¶

Discretization Errors ¶

Spatial Discretization:

FFT-based methods have spectral accuracy for smooth functions:

\[\epsilon_{\text{space}} \sim e^{-\alpha N}\]

for some \(\alpha > 0\). However, discontinuities cause Gibbs oscillations.

Temporal Discretization:

For the split-step method:

\[\epsilon_{\text{time}} \sim C\Delta t^3\]

where \(C\) depends on potential derivatives.

Aliasing ¶

High-frequency components can alias to low frequencies due to finite grid spacing. To avoid:

Ensure wavepacket has negligible amplitude at \(k > k_{\text{Nyquist}} = \pi/\Delta x\)
Filter high-frequency components if necessary
Use sufficiently fine grid: \(\Delta x < \lambda_{\min}/10\)

Memory and Performance ¶

Memory Scaling:

1D: \(O(N)\) complex numbers
2D: \(O(N^2)\) complex numbers
3D: \(O(N^3)\) complex numbers

With complex64 (8 bytes each), a \(1024^3\) 3D grid requires ~8 GB.

Computational Complexity:

Per time step:

FFT: \(O(N^d \log N)\) where \(d\) is dimension
Potential application: \(O(N^d)\)

GPU acceleration provides 10-100× speedup over CPU for large grids.