Schrödinger equation as real and imaginary fields
A Gaussian wave packet that follows the Schrödinger equation, as its real and imaginary parts.
- Name:
schrodinger - Field:
physics/waves - Formalism: ordinary differential equations
- References: Schrödinger 1926; Schwichtenberg 2020
Factory
@register_model(
field='physics/waves',
tags=('waves', 'field', 'pde', 'quantum', 'spectral'),
references=('Schrödinger 1926', 'Schwichtenberg 2020'),
defaults={'until': 4.0, 'dt': 0.1, 'solver': ETDRK4Solver(dt=0.01)},
)
def schrodinger(
n_sites: int = 256,
*,
length: float = 40.0,
D: float = 0.5,
momentum: float = 2.0,
width: float = 1.0,
omega: float = 0.0,
formulation: str = 'fields',
gamma: float = 0.0,
) -> System:
"""Return a Gaussian wave packet that follows the Schrödinger equation, as its real and imaginary parts.
The wave function `psi = u + i w` follows `psi_t = i D laplacian(psi) - i V psi`, written as `u_t = -D laplacian(w)
+ V w` and `w_t = D laplacian(u) - V u`, on a periodic lattice with spectral derivatives. In units with `ħ = 1`,
`D = 1 / (2 m)`. The potential is the harmonic well `V = omega² x² / (4 D)`, or zero for a free packet. A free
packet of momentum `k` moves at the group velocity `2 D k` and spreads, and its norm `sum(u² + w²)` is conserved.
The operator formulation builds the same particle as one quantum subsystem whose basis states are the sites, with
the Hamiltonian `-D laplacian + V` of the same spectral Laplacian, normalised so that `|psi|²` is the probability of
every site. It runs with the quantum solvers, and it accepts position dephasing at the rate `gamma`, which turns
the state into a density matrix.
Args:
n_sites:
The number of sites of the periodic lattice.
length:
The length of the periodic box, centered on 0.
D:
The coefficient of the Laplacian, `ħ / (2 m)`.
momentum:
The mean wave number of the packet.
width:
The width of the packet.
omega:
The angular frequency of the harmonic well.
formulation:
`'fields'` for two real fields, or `'operators'` for a particle on the sites of the lattice.
gamma:
The rate of position dephasing of the operator formulation.
Returns:
system:
The model, with the real part `packet.u` and the imaginary part `packet.w` of the wave function, or the
quantum state `schrodinger.state` of the operator formulation.
Raises:
ValidationError:
The formulation is unknown, or the field formulation has a dephasing rate.
Examples:
>>> from skcomplex.models.physics.waves import schrodinger
>>> schrodinger(64).elements[0].check().prepare().size_
64
"""
spacing = length / n_sites
positions = spacing * np.arange(n_sites) - length / 2
envelope = np.exp(-(positions**2) / (4 * width**2)) / (2 * np.pi * width**2) ** 0.25
psi = envelope * np.exp(1j * momentum * positions)
potential = omega**2 * positions**2 / (4 * D)
space = LatticeSpace(n_sites, spacing=spacing)
if formulation not in FORMULATIONS or (formulation == 'fields' and gamma):
error_msg = (
f'Unknown formulation {formulation!r}, or dephasing in the field formulation. Use one of '
f'{", ".join(FORMULATIONS)}, and gamma > 0 only with the operator formulation.'
)
raise ValidationError(error_msg)
if formulation == 'operators':
amplitudes = psi / np.linalg.norm(psi)
interactions: list[BaseInteraction] = [
LatticeHamiltonianInteraction('kinetic', target='packet', D=D),
]
if gamma > 0:
interactions.append(
PositionDephasingInteraction('dephasing', target='packet', gamma=gamma, n_sites=n_sites),
)
return System(
'schrodinger',
elements=[SiteElement('packet', amplitudes, hilbert=site_space(space, potential, projectors=gamma > 0))],
interactions=interactions,
)
packet = FieldElement(
'packet',
element=WaveFunctionElement,
space=space,
values={'u': psi.real, 'w': psi.imag, 'D': D, 'V': potential},
discretization='spectral',
)
return System(elements=[packet])
Components
class WaveFunctionElement(Element):
"""Site of the wave function `psi = u + i w` of a particle with the potential `V`, by its real and imaginary parts.
The parts follow the Schrödinger equation `i psi' = -D laplacian(psi) + V psi`, so
`u' = -D laplacian(w) + V w` and `w' = D laplacian(u) - V u`.
Args:
label:
The label of the site.
visual:
How the site is drawn, or `None` for none.
u:
The initial real part.
w:
The initial imaginary part.
D:
The coefficient of the kinetic energy, `hbar / (2 m)`.
V:
The potential at the site.
Examples:
>>> from skcomplex.models.physics.waves import WaveFunctionElement
>>> WaveFunctionElement('site', V=1.0).V
1.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
u: Annotated[float, StateSpec()] = 0.0,
w: Annotated[float, StateSpec()] = 0.0,
D: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
V: Annotated[float, ParameterSpec()] = 0.0,
) -> None:
super().__init__(label, visual=visual)
self.u = u
self.w = w
self.D = D
self.V = V
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the Schrödinger equation of the real and imaginary parts.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The rate laws of the two parts.
"""
return [
RateLaw(s.u, -s.D * laplacian(s.w) + s.V * s.w),
RateLaw(s.w, s.D * laplacian(s.u) - s.V * s.u),
]
Equations
packet.u[i]: diff(u_packet(t), t) == -D_packet*laplacian(w_packet(t)) + V_packet*w_packet(t)
packet.w[i]: diff(w_packet(t), t) == D_packet*laplacian(u_packet(t)) - V_packet*u_packet(t)