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.

.. contents:: Table of Contents
   :local:
   :depth: 3

Quantum Mechanics Foundations
------------------------------

Hilbert Space Formalism
~~~~~~~~~~~~~~~~~~~~~~~~

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

**Inner Product:**

.. math::

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

**Normalization Condition:**

.. math::

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

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

Operators and Observables
~~~~~~~~~~~~~~~~~~~~~~~~~~

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

**Position Operator:**

.. math::

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

**Momentum Operator:**

.. math::

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

**Hamiltonian (Energy) Operator:**

.. math::

   \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 :math:`\hat{A}` in state :math:`|\psi\rangle`:

.. math::

   \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:

.. math::

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

In position representation:

.. math::

   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:

.. math::

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

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

Conservation Laws
^^^^^^^^^^^^^^^^^

**Probability Conservation:**

.. math::

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

**Energy Conservation (time-independent H):**

.. math::

   \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:

.. math::

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

where

.. math::

   \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:

.. math::

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

.. math::

   \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:

.. math::

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

Apply kinetic evolution:

.. math::

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

Transform back to position space:

.. math::

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

**Step 3: Half potential step (position space)**

.. math::

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

Advantages
^^^^^^^^^^

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

Stability and Accuracy
^^^^^^^^^^^^^^^^^^^^^^

**CFL Condition:**

For stability, the time step should satisfy:

.. math::

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

where :math:`k_{\max} = \pi/\Delta x` is the Nyquist wavevector.

**Spatial Resolution:**

To accurately resolve features with characteristic momentum :math:`k_0`:

.. math::

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

A typical choice is :math:`\Delta x = \lambda/10` where :math:`\lambda = 2\pi/k_0`.

Fourier Transform Conventions
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

The package uses the following Fourier transform conventions:

**Forward Transform:**

.. math::

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

**Inverse Transform:**

.. math::

   \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 :math:`N` points and spacing :math:`\Delta x`:

.. math::

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

The momentum grid is:

.. math::

   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}

**Parseval's Theorem (Norm Conservation):**

.. math::

   \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:

.. math::

   \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:

.. math::

   \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 :math:`\partial_t|\psi\rangle = -i\hat{H}|\psi\rangle/\hbar`:

.. math::

   \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}

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

**Disadvantages:** :math:`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:

.. math::

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

where

.. math::

   \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}

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

**Properties:**

* Smooth transition minimizes reflections
* Higher :math:`n` gives sharper absorption
* Width :math:`w` should be several wavelengths

Complex Absorbing Potential (CAP)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Add an imaginary potential in the absorption region:

.. math::

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

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

**Effect:** Exponentially damps the wavefunction amplitude:

.. math::

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

Periodic Boundary Conditions
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

The FFT naturally imposes periodic boundary conditions:

.. math::

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

where :math:`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:

.. math::

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

satisfying the commutation relations:

.. math::

   [\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 :math:`|n_1, n_2, \ldots, n_M\rangle` where :math:`n_i` is the number of photons in mode :math:`i`.

**Action of Operators:**

.. math::

   \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}

Free Field Hamiltonian
^^^^^^^^^^^^^^^^^^^^^^

.. math::

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

where :math:`\omega_i` is the frequency of mode :math:`i`.

Light-Matter Interaction
~~~~~~~~~~~~~~~~~~~~~~~~~

Jaynes-Cummings Model
^^^^^^^^^^^^^^^^^^^^^^

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

.. math::

   \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:

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

**Rotating Wave Approximation (RWA):**

The terms :math:`\hat{\sigma}_+\hat{a}^\dagger` and :math:`\hat{\sigma}_-\hat{a}` are neglected (counter-rotating terms) when :math:`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 :math:`q` interacting with the electromagnetic field:

.. math::

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

where :math:`\hat{\mathbf{E}}` is the quantized electric field operator.

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

.. math::

   \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 :math:`\hat{\mathbf{e}}_i` is the polarization vector.

Fock Space Truncation
~~~~~~~~~~~~~~~~~~~~~~

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

.. math::

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

**Dimension:**

For :math:`M` modes with maximum :math:`N_{\max}` photons each:

.. math::

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

**Convergence:**

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

Example Systems
---------------

Harmonic Oscillator
~~~~~~~~~~~~~~~~~~~

The 1D quantum harmonic oscillator has Hamiltonian:

.. math::

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

**Energy Eigenvalues:**

.. math::

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

**Eigenfunctions:**

.. math::

   \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 :math:`H_n` are Hermite polynomials.

**Coherent States:**

Displaced Gaussian wavepackets are approximate eigenstates of the annihilation operator:

.. math::

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

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

Tunneling Through a Barrier
~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Rectangular Barrier
^^^^^^^^^^^^^^^^^^^

Potential:

.. math::

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

For :math:`E < V_0`, the transmission coefficient is:

.. math::

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

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

**Tunneling Time:**

The time for a wavepacket to tunnel through is approximately:

.. math::

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

Double-Slit Interference
~~~~~~~~~~~~~~~~~~~~~~~~~

Physical Parameters
^^^^^^^^^^^^^^^^^^^

* Slit separation: :math:`d`
* Slit width: :math:`w`
* Wavelength: :math:`\lambda = 2\pi/k`
* Screen distance: :math:`L`

**Fringe Spacing:**

For :math:`w \ll d \ll L`:

.. math::

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

**Visibility:**

The fringe visibility depends on the coherence length :math:`\ell_c`:

.. math::

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

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

Atomic Units
------------

Fundamental Constants
~~~~~~~~~~~~~~~~~~~~~

.. list-table::
   :header-rows: 1
   :widths: 30 30 40

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

Unit Conversions
~~~~~~~~~~~~~~~~

**Length:**

.. math::

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

**Energy:**

.. math::

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

**Time:**

.. math::

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

**Electric Field:**

.. math::

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

Example: Hydrogen Ground State
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

In atomic units:

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

Numerical Considerations
------------------------

Discretization Errors
~~~~~~~~~~~~~~~~~~~~~

**Spatial Discretization:**

FFT-based methods have spectral accuracy for smooth functions:

.. math::

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

for some :math:`\alpha > 0`. However, discontinuities cause Gibbs oscillations.

**Temporal Discretization:**

For the split-step method:

.. math::

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

where :math:`C` depends on potential derivatives.

Aliasing
~~~~~~~~

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

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

Memory and Performance
~~~~~~~~~~~~~~~~~~~~~~

**Memory Scaling:**

* 1D: :math:`O(N)` complex numbers
* 2D: :math:`O(N^2)` complex numbers
* 3D: :math:`O(N^3)` complex numbers

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

**Computational Complexity:**

Per time step:

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

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