Electron flow in water, one trajectory at a time

The state

Pick two orbitals φ₁, φ₂ with orbital energies ε₁, ε₂, and a weight w. The page follows one electron in the superposition (atomic units, ħ = m = 1):

ψ(r, t) = √(1 − w) φ₁(r) e^−iε₁t + √w φ₂(r) e^−iε₂t

Each orbital alone is a stationary state: a real function times a phase. The superposition is not, because the two phases turn at different rates. The orbitals are restricted Hartree–Fock solutions in the cc-pVDZ basis at the experimental geometry, computed with PySCF: five occupied (1a₁ core, 2a₁, 1b₂, 3a₁, 1b₁ HOMO) and three empty O–H antibonding orbitals (4a₁ LUMO, 2b₂, 3b₂).^{[1, 2, 3, 4]}

The density beats

With a = √(1 − w), b = √w, and θ = Δε t where Δε = ε₂ − ε₁:

ρ(r, t) = a²φ₁² + b²φ₂² + 2ab φ₁φ₂ cos θ

Only the cross term, the transition density φ₁φ₂, oscillates. The cloud alternates between (aφ₁ + bφ₂)² and (aφ₁ − bφ₂)² with period T = 2π/Δε, shown in the readout as h/Δε. Symmetry decides which way the charge moves: the oscillating dipole points along the axis that carries the symmetry of the product of the two orbitals.

How each particle moves

Each point moves with Bohm's guiding equation, v = Im(∇ψ/ψ) = j/ρ.^{[5, 6, 7]} For this state it works out to

v(r, t) = −ab sin θ · K(r) / ρ(r, t),   K = φ₁∇φ₂ − φ₂∇φ₁

• Fixed directions. The velocity always points along ±K(r), which does not change in time, so each particle runs back and forth along one fixed streamline.
• Exact return. ρ depends on cos θ, which is symmetric about θ = π, while the speed factor sin θ is antisymmetric. The second half of each period retraces the first, and every particle is back where it started at t = T.
• Fast swirl, slow beat. Typical speeds are about 1 bohr per atomic time unit even when the beat takes femtoseconds. Most of the current circulates: its curl, 2∇φ₁ × ∇φ₂, is large, while its divergence, which actually moves density, is proportional to Δε.
• Vortex lines. Where the nodal surfaces of both orbitals cross, ρ vanishes at all times. Near such a line K grows like r and ρ like r², so particles whirl around it with speed ∝ 1/r.

What the points are

They are not many electrons. Each point is where the same single electron could be in a separate run of the experiment. The points start distributed as |ψ(0)|², and Bohmian dynamics is equivariant: an ensemble that starts as |ψ|² stays |ψ(t)|², so the cloud at any moment is the quantum density.^[8] The full ten-electron problem lives in 30 dimensions; this is the one-electron picture.

Colour comes from position, not from a label the point carries: it is the share a²φ₁² / (a²φ₁² + b²φ₂²) at the point's location. A point changes colour when it travels into the other orbital's region.

Keeping the cloud exact: the Poisson correction

Equivariance needs continuity, ∂ρ/∂t + ∇·j = 0, which for this state requires

∇·K = −2Δε φ₁φ₂

That identity holds only for exact eigenstates of a local Hamiltonian. Hartree–Fock orbitals in a Gaussian basis miss it point by point by a factor of two to seven, because of nonlocal exchange and the finite basis. Plain Bohm flow still gets the beat right qualitatively, but between full periods the cloud drifts away from |ψ(t)|².

Bohm + Poisson adds the smallest density-weighted correction that restores continuity:

K → K + σ∇χ,   σ = ½(φ₁² + φ₂²),   ∇·(σ∇χ) = −2Δε φ₁φ₂ − ∇·K

• Both sides of continuity share the factor sin θ, so one time-independent χ fixes every moment of the beat. It is solved once per orbital pair, ahead of time.
• Weighting by σ keeps the added current where the electron is. Without it the correction leaks into empty space, where v = j/ρ blows up.
• χ is a 229-term expansion (Gaussians plus polynomials) solved on a molecular integration grid.^[9] Continuity then holds exactly for every density moment up to fifth order, though not point by point.

Checked by running the particles for half a period and comparing with fresh samples of |ψ(T/2)|²:

Single orbitals: Bohm and Nelson

A single orbital is real, so Im(∇ψ/ψ) = 0: in Bohmian mechanics every particle stands still. The same holds for the whole ten-electron ground state, a real Slater determinant. The Nelson option shows a different interpretation, stochastic mechanics, where each particle drifts along ∇ψ/ψ and diffuses with ħ/2m.^[10] The cloud stays exactly |ψ|² (each step is a Metropolis-adjusted Langevin move), but individual paths are random.^[11]

Numerics

• Starting positions: importance sampling from Gaussians centred on the nuclei, then Metropolis sweeps with random step sizes.
• Time stepping: adaptive midpoint steps of at most 0.02 bohr, needed for the 1/r speeds near vortex lines.
• Particles are split across several Web Workers; the display interpolates between Bohmian steps.
• 1× speed is one atomic time unit (24.2 as) per second. When the integrator cannot keep up, the readout shows the actual rate.

What it means, and where the model stops

• Two occupied orbitals make a hole wavepacket. Sudden ionization can leave H₂O⁺ in a coherent superposition of, say, 1b₁⁻¹ and 3a₁⁻¹; the cloud shows where the missing electron is. This is attosecond charge migration.^[12] Orbital-energy gaps are fair estimates here: 2.00 eV for 1b₁–3a₁ against about 2.1 eV from measured ionization energies (12.62 and about 14.7 eV).^[13]
• An occupied plus an empty orbital make an excitation wavepacket. The shapes are reasonable, but Hartree–Fock virtual energies are far too high: the 18.4 eV gap for 1b₁ → 4a₁ is well above water's first absorption band near 7.4 eV, so the true beat is slower.^[14] Virtual orbitals also depend strongly on the basis set.
• Left out: electron correlation (one electron in the mean field of the others), nuclear motion and the decoherence it causes within a few femtoseconds, and Pauli exclusion, which makes a superposition of two doubly occupied orbitals meaningful only as a hole.