physics-schrodinger

Answer

"""
Derivation of the Time-Dependent Schrödinger Equation from the de Broglie Hypothesis
=====================================================================================

Step-by-step:
  1. de Broglie hypothesis       E = ħω,  p = ħk
  2. Plane wave ansatz            Ψ(x,t) = A·exp(i(kx - ωt))
  3. Partial derivatives →        energy operator:   Ê = iħ ∂/∂t
     quantum operators            momentum operator: p̂ = -iħ ∂/∂x
  4. Classical Hamiltonian        H = p²/(2m) + V(x) → substitute operators
  5. Time-dependent SE            iħ ∂Ψ/∂t = [ -ħ²/(2m) ∂²/∂x² + V(x) ] Ψ
  6. Probability current          j = (ħ/m) Im(Ψ* ∂Ψ/∂x),  continuity: ∂ρ/∂t + ∂j/∂x = 0
  7. Free particle                V=0 → iħ ∂Ψ/∂t = -ħ²/(2m) ∂²Ψ/∂x²
"""

import numpy as np

# ── 1. de Broglie relations & plane wave ansatz ──────────────────────────
hbar = 1.0   # reduced Planck constant (atomic units)
m    = 1.0   # particle mass
k    = 2.0   # wavenumber
omega = hbar * k**2 / (2 * m)   # free-particle dispersion: E = ħω = ħ²k²/(2m)

def psi_plane(x: float, t: float) -> complex:
    """Plane wave: the fundamental building block of the derivation."""
    return np.exp(1j * (k * x - omega * t))

# ── 2. Verify operator relations ────────────────────────────────────────
x_test, t_test = 1.0, 0.5
dx, dt = 1e-5, 1e-5

# Time derivative → energy operator
dpsi_dt  = (psi_plane(x_test, t_test+dt) - psi_plane(x_test, t_test-dt)) / (2*dt)
E_op     = 1j * hbar * dpsi_dt
E_actual = hbar * omega * psi_plane(x_test, t_test)
assert np.isclose(abs(E_op), abs(E_actual), rtol=1e-4), "Energy operator mismatch"

# Spatial derivative → momentum operator
dpsi_dx  = (psi_plane(x_test+dx, t_test) - psi_plane(x_test-dx, t_test)) / (2*dx)
p_op     = -1j * hbar * dpsi_dx
p_actual = hbar * k * psi_plane(x_test, t_test)
assert np.isclose(abs(p_op), abs(p_actual), rtol=1e-4), "Momentum operator mismatch"

# Second derivative → kinetic energy
d2psi_dx2 = (psi_plane(x_test+dx, t_test) - 2*psi_plane(x_test, t_test)
             + psi_plane(x_test-dx, t_test)) / dx**2
T_op      = -hbar**2 / (2*m) * d2psi_dx2
T_actual  = (hbar**2 * k**2 / (2*m)) * psi_plane(x_test, t_test)
assert np.isclose(abs(T_op), abs(T_actual), rtol=1e-4), "Kinetic operator mismatch"

# ── 3. TDSE verification: iħ ∂Ψ/∂t = -ħ²/(2m) ∂²Ψ/∂x² ──────────────
lhs = E_op
rhs = T_op
residual = lhs - rhs
assert abs(residual) < 1e-4, f"TDSE violated: |residual| = {abs(residual)}"

# ── 4. Probability current & continuity ────────────────────────────────
# j(x,t) = (ħ/m) Im(Ψ* ∂Ψ/∂x)
def current(x: float, t: float) -> float:
    dpsi_dx_h = (psi_plane(x+dx/2, t) - psi_plane(x-dx/2, t)) / dx
    return float((hbar / m) * np.imag(np.conj(psi_plane(x, t)) * dpsi_dx_h))

def continuity_residual(x: float, t: float) -> float:
    drho_dt = (abs(psi_plane(x, t+dt))**2 - abs(psi_plane(x, t-dt))**2) / (2*dt)
    dj_dx   = (current(x+dx, t) - current(x-dx, t)) / (2*dx)
    return drho_dt + dj_dx

# For plane wave: ρ=1 (constant), j=ħk/m (constant), ∂ρ/∂t+∂j/∂x=0
assert abs(continuity_residual(0.0, 0.0)) < 1e-4, "Continuity violated"

# ── 5. Free particle Hamiltonian form ──────────────────────────────────
print("TDSE derivation verified!")
print(f"  ψ(x,t) = exp(i({k}x - {omega}t))")
print(f"  Energy operator:   iħ ∂/∂t → eigenvalue ħω = {hbar*omega}")
print(f"  Momentum operator: -iħ ∂/∂x → eigenvalue ħk = {hbar*k}")
print(f"  Hamiltonian: Ĥ = -ħ²/(2m) ∂²/∂x²  (free particle, V=0)")
print(f"  TDSE residual: |iħ∂ψ/∂t - Ĥψ| = {float(abs(residual)):.2e}")
print(f"  Continuity:     |∂ρ/∂t + ∂j/∂x| = {abs(continuity_residual(0.0,0.0)):.2e}")