"""Full QFT propagation: Klein‑Gordon and Dirac solvers with empirical parameters."""
from typing import Optional
import numpy as np
from scipy.fft import fft, fftfreq, ifft
from mercurial.params.empirical import QFTParams
[docs]
class KleinGordon2D:
"""
2D Klein‑Gordon field solver (real scalar field) using FDTD,
with parameters from empirical database.
"""
[docs]
def __init__(
self,
nx: int,
ny: int,
dx: Optional[float] = None,
m: Optional[float] = None,
params: Optional[QFTParams] = None,
initial_phi: Optional[np.ndarray] = None,
initial_pi: Optional[np.ndarray] = None,
absorbing_boundary: bool = True,
):
"""
Parameters
----------
nx, ny : int
Grid dimensions.
dx : float, optional
Spatial step. If None, uses params.dx.
m : float, optional
Mass (natural units). If None, uses params.mass_scalar.
params : QFTParams, optional
Empirical parameter dataclass. If None, uses default.
initial_phi, initial_pi : array, optional
Initial field and conjugate momentum.
absorbing_boundary : bool
Apply absorbing boundary conditions (PML‑like).
"""
if params is None:
params = QFTParams()
self.dx = dx if dx is not None else params.dx
self.m = m if m is not None else params.mass_scalar
self.absorb = absorbing_boundary
# Set time step from Courant condition (c=1)
self.dt = self.dx / np.sqrt(2) # stable for 2D Laplacian
self.nx = nx
self.ny = ny
# Field arrays (staggered in time)
self.phi = np.zeros((nx, ny))
self.pi = np.zeros((nx, ny))
if initial_phi is not None:
self.phi = initial_phi
if initial_pi is not None:
self.pi = initial_pi
# PML damping profile
self.pml_sigma = np.zeros((nx, ny))
if absorbing_boundary:
self._init_pml()
def _init_pml(self, width: int = 10, strength: float = 0.1):
for i in range(width):
sigma = strength * ((width - i) / width) ** 3
self.pml_sigma[i, :] += sigma
self.pml_sigma[self.nx - 1 - i, :] += sigma
self.pml_sigma[:, i] += sigma
self.pml_sigma[:, self.ny - 1 - i] += sigma
[docs]
def laplacian(self, field: np.ndarray) -> np.ndarray:
lap = np.zeros_like(field)
lap[1:-1, 1:-1] = (
field[2:, 1:-1]
+ field[:-2, 1:-1]
+ field[1:-1, 2:]
+ field[1:-1, :-2]
- 4 * field[1:-1, 1:-1]
) / self.dx**2
return lap
[docs]
def step(self, source: Optional[np.ndarray] = None) -> None:
self.phi += self.pi * self.dt
lap_phi = self.laplacian(self.phi)
dpi = (lap_phi - self.m**2 * self.phi) * self.dt
if source is not None:
dpi += source * self.dt
self.pi += dpi
if self.absorb:
self.pi *= 1 - self.pml_sigma * self.dt
[docs]
def energy_density(self) -> np.ndarray:
grad_sq = 0.5 * (
(np.gradient(self.phi, self.dx, axis=0)) ** 2
+ (np.gradient(self.phi, self.dx, axis=1)) ** 2
)
return 0.5 * (self.pi**2 + grad_sq + self.m**2 * self.phi**2)
[docs]
def correlation(self, r: int = 0) -> float:
phi_flat = self.phi.flatten()
if r > 0 and len(phi_flat) > r:
return np.corrcoef(phi_flat[:-r], phi_flat[r:])[0, 1]
return 1.0
[docs]
class Dirac1D:
"""
1+1D Dirac field solver using split‑operator method, with empirical mass.
"""
[docs]
def __init__(
self,
nx: int,
dx: Optional[float] = None,
m: Optional[float] = None,
params: Optional[QFTParams] = None,
initial_psi: Optional[np.ndarray] = None,
):
"""
Parameters
----------
nx : int
Number of spatial points.
dx : float, optional
Spatial step. If None, uses params.dx.
m : float, optional
Mass. If None, uses params.mass_dirac.
params : QFTParams, optional
Empirical parameter dataclass.
initial_psi : array, shape (2, nx), optional
Initial spinor (two components).
"""
if params is None:
params = QFTParams()
self.dx = dx if dx is not None else params.dx
self.m = m if m is not None else params.mass_dirac
self.dt = self.dx / 2.0 # stable for split‑operator
self.nx = nx
# Precompute momentum grid
self.k = 2 * np.pi * fftfreq(nx, self.dx)
if initial_psi is None:
# Gaussian wavepacket
x = np.linspace(-5, 5, nx)
sigma = 0.5
k0 = 5.0
psi0 = np.exp(-((x) ** 2) / (2 * sigma**2)) * np.exp(1j * k0 * x)
self.psi = np.array([psi0, psi0])
else:
self.psi = initial_psi.copy()
[docs]
def step(self, dt: Optional[float] = None) -> None:
if dt is None:
dt = self.dt
# Half mass
self.psi[0] *= np.exp(-1j * self.m * dt / 2)
self.psi[1] *= np.exp(1j * self.m * dt / 2)
# Kinetic step (Fourier)
for c in range(2):
psi_k = fft(self.psi[c])
psi_k *= np.exp(-1j * self.k * dt)
self.psi[c] = ifft(psi_k)
# Half mass again
self.psi[0] *= np.exp(-1j * self.m * dt / 2)
self.psi[1] *= np.exp(1j * self.m * dt / 2)
[docs]
def probability_density(self) -> np.ndarray:
return np.abs(self.psi[0]) ** 2 + np.abs(self.psi[1]) ** 2
[docs]
def current(self) -> np.ndarray:
return 2 * np.real(self.psi[0].conj() * self.psi[1])