T>T: The Reality of Electron Density in a World of Orbitals
Orbitals Are Mathematical
Pick up almost any introductory chemistry textbook and you will find an orbital depicted as a tangible object: a lobed shape around an atom, a probability cloud you could in principle photograph. I entered the insane world of chemistry through the side door and found it interesting that chemists speak of electrons being in 1s orbitals, of bonds formed by the overlap of p orbitals, of hybridisation as if it were a physical process the atom performs. This language is useful but fundamentally false.
An orbital is a single-particle wavefunction, a mathematical object that describes one electron moving in the averaged field of all the others. The operative word is averaged. The orbital is constructed by erasing the moment-to-moment correlations between electrons and replacing the other $N-1$ electrons with a smooth, static background potential. The orbital never existed as a physical object that could be measured; it only appears when we choose to describe a many-electron system using a single-product approximation.
What nature actually provides which we can measure is the electron density $\rho(\mathbf{r})$: the probability of finding an electron at position $\mathbf{r}$, regardless of where all the others happen to be. The Hohenberg-Kohn theorem (1964) places this on firm mathematical ground: the ground-state electron density uniquely determines the external potential and therefore all properties of the system. The density is not derived from orbitals; it is primary.
This distinction carries real physical consequences, and the best way to see them is to look at what happens to the density when one electron moves, comparing the prediction from Hartree-Fock (the archetypal orbital-based method) with a fully correlated wavefunction.
An Orchestra in a Dark Room
I have heard many analogies used describing Hartree Fock’s handling of electron correlation, one of which being that it resembles two ships passing each other in the night. I raise issue with this as each ship knows the other ship is there via their wakes. If this analogy were to hold, then the two ships would also be able to pass directly through one another.
Imagine an orchestra performing in a very large, pitch-black room. No musician can see anyone else or track individual players directly. Instead, each musician plays according to a continuously updated score generated from the total sound of the orchestra itself.
This score is not arbitrary: it encodes both the overall harmony (the average sound of everyone else) and strict rules that prevent certain musicians, those playing the same instrument, from producing identical notes in the same way. These rules subtly reshape how each part is played, even though no musician is explicitly tracking another.
All musicians adjust their playing in response to this shared score, and the score updates in response to their playing. This continues until the performance stabilises. No single musician can change their part without making the overall sound worse.
The final performance is therefore the best possible arrangement given that the entire orchestra is constrained to follow a single, globally consistent score derived from itself.
This is Hartree-Fock, in essence, one might say that it is not possible to use a classical analogue for a quantum mechanical situation, which is true, but this is the way that I visualise it. Each electron moves in the mean field generated by all the other electrons. Electron 1 sees a smooth charge cloud representing electrons 2, 3, …, $N$ smeared out across space. It does not know where electron 2 is right now. The interaction is real, the Coulomb repulsion is included, but it is averaged over all possible positions of the other electrons.
A real orchestra, of course, does react to individual players in real time. When the oboist enters late, the conductor adjusts; when the violins drift sharp, everyone else compensates. The musicians are correlated; their precise moment-to-moment positions and actions are coupled together in ways that the mean sound alone cannot capture.
Fully correlated quantum chemistry methods are the analogue: the wavefunction depends explicitly on the instantaneous positions of all electrons simultaneously, including the inter-electron distance $r_{12}$. When electron 1 is close to electron 2, the wavefunction reduces the probability amplitude for them to be nearby, not because the orbital changed, but because the full wavefunction encodes their correlated motion directly.
What Correlation Actually Means
The definition used throughout this post is the one due to Löwdin (1959)
\[E_\text{corr} = E_\text{exact} - E_\text{HF}\]Electron correlation is the difference between the exact non-relativistic energy and the Hartree-Fock energy in the limit of a complete basis set.
This is a negative number: Hartree-Fock always overestimates the energy, because averaging the electron-electron repulsion is never as good as accounting for the fact that electrons actively avoid each other. The correlation energy is typically small relative to the total energy (a fraction of a percent) but it is chemically decisive: it governs bond dissociation energies, reaction barriers, van der Waals interactions, and the accurate prediction of thermodynamic quantities.
The concept also sharpens what we mean by “uncorrelated”. A wavefunction is uncorrelated if and only if it can be written as a product of single-electron functions:
\[\Psi_\text{HF}(\mathbf{r}_1, \mathbf{r}_2) = \psi(\mathbf{r}_1)\,\psi(\mathbf{r}_2)\](ignoring spin antisymmetry for clarity). For such a state, the probability of finding electron 1 at $\mathbf{r}_1$ and simultaneously finding electron 2 at $\mathbf{r}_2$ factorises:
\[P(\mathbf{r}_1, \mathbf{r}_2) = |\psi(\mathbf{r}_1)|^2\,|\psi(\mathbf{r}_2)|^2 = \rho(\mathbf{r}_1)\,\rho(\mathbf{r}_2)\]The positions of the two electrons are statistically independent. Where electron 1 sits carries zero information about where electron 2 is likely to be. This is the defining failure of HF.
A correlated wavefunction, by contrast, cannot be factorised. The canonical example is any function that depends explicitly on $r_{12} = |\mathbf{r}_1 - \mathbf{r}_2|$:
\[\Psi_\text{FC}(\mathbf{r}_1, \mathbf{r}_2, r_{12}) \neq f(\mathbf{r}_1)\,g(\mathbf{r}_2)\]The joint density $|\Psi_\text{FC}|^2$ does not factorise, so knowing where electron 1 is does change the distribution of electron 2. This is the correlation hole: a region near electron 1 where the probability of finding electron 2 is reduced relative to the uncorrelated prediction.
Two Wavefunctions for Helium
Helium is the ideal system for this comparison. It has exactly two electrons, a known nuclear charge $Z = 2$, and an infinite nuclear mass (to excellent approximation). The non-relativistic Hamiltonian in atomic units is
\[\hat{H} = -\frac{1}{2}\nabla_1^2 - \frac{1}{2}\nabla_2^2 - \frac{Z}{r_1} - \frac{Z}{r_2} + \frac{1}{r_{12}}\]where $r_1$, $r_2$ are the electron-nucleus distances and $r_{12}$ is the interelectronic distance.
The Hartree-Fock wavefunction
The HF ground state is a symmetric (singlet) spatial wavefunction:
\[\Psi_\text{HF}(\mathbf{r}_1, \mathbf{r}_2) = \psi(r_1)\,\psi(r_2)\]where the orbital $\psi$ is expanded in a Laguerre basis with nonlinear parameter $A$:
\[\psi(r) = e^{-Ar/2}\sum_{k=0}^{N-1} c_k\, L_k(Ar)\]$L_k$ are standard Laguerre polynomials. The coefficients $c_k$ and parameter $A$ are optimised variationally. For helium with $N=25$ basis functions, the optimised value is $A = 5.236$, giving energy $E_\text{HF} = -2.86168$ hartree.
The HF conditional density, the probability of finding electron 2 at position $\mathbf{r}_2$ given that electron 1 is at $\mathbf{r}_1$ is
\[\rho_\text{HF}(\mathbf{r}_2\,|\,\mathbf{r}_1) \propto |\psi(r_2)|^2\]The conditional density is independent of $\mathbf{r}_1$. Moving electron 1 anywhere in space has zero effect on the predicted distribution of electron 2. The two electrons are statistically decoupled.
The fully correlated wavefunction
The fully correlated (FC) wavefunction is the same Pekeris-style expansion described in detail in A Modern Implementation of Pekeris’ Helium Calculation; the derivation and variational optimisation are covered there. The brief summary here is sufficient to understand the visualisation.
It uses perimetric coordinates, which are combinations of the three inter-particle distances that are always non-negative (satisfying the triangle inequality):
\[u = A(r_{12} + r_2 - r_1), \quad v = B(r_{12} + r_1 - r_2), \quad w = C(r_1 + r_2 - r_{12})\]For the symmetric helium ground state, $A = B$ and $C = 2A$, giving $C = A + B$. The wavefunction is:
\[\Psi_\text{FC}(r_1, r_2, r_{12}) = e^{-A(r_1 + r_2)}\sum_{i,j,k} c_{ijk}\, L_i(u)\,L_j(v)\,L_k(w)\]This triple Laguerre expansion explicitly depends on $r_{12}$, so the wavefunction, and hence the density, knows about the instantaneous separation between the two electrons. With $\Omega = 30$ (sum of indices $i + j + k \leq 30$) the basis contains 2856 symmetry-unique functions, and the optimised energy is $E_\text{FC} = -2.903724$ hartree. The exact non-relativistic energy for helium is $E_\text{exact} = -2.903724$ hartree to this precision, so the FC wavefunction is essentially exact.
The correlation energy is:
\[E_\text{corr} = E_\text{FC} - E_\text{HF} = -2.903724 - (-2.861680) = -0.042044 \text{ hartree} \approx -26 \text{ kcal/mol}\]Not a small number: this is the energy of a moderately strong chemical bond.
The FC conditional density does depend on $\mathbf{r}_1$:
\[\rho_\text{FC}(\mathbf{r}_2\,|\,\mathbf{r}_1) \propto |\Psi_\text{FC}(r_1, r_2, r_{12})|^2\]As electron 1 moves away from the nucleus, the density of electron 2 responds; it is no longer spherically centred on the nucleus alone, but shifts to partially screen electron 1’s departure and avoid the region near electron 1.
Building the Visualiser
I have been on a bit of a visuliser vibe recently and want to create a browser-based tool where a reader drags a slider to move electron 1 and immediately sees the conditional density of electron 2 respond (or not respond, in the HF case). There is no server; the HTML file must be entirely self-contained. The approach has two phases: compute everything in Python at script time, embed the results in the HTML, then let the browser handle only lightweight visibility toggling at interaction time.
Reading the wavefunction files
Both wavefunctions come from Maple .sv (save) files produced by a variational optimisation. Writing a full Maple parser would be excessive, the relevant data appears in a predictable, repeating format. A targeted regular expression is enough.
The HF orbital is stored in a block headed F1:= and terminated by +0:, with each term written as (coefficient)*L(order,u). The exponent $A$ appears on its own line as A:=5.236...:. Parsing it is straightforward:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
def parse_hf(path):
with open(path) as fh:
content = fh.read()
A = float(re.search(r"\nA:=([\d.eE+\-]+):", content).group(1))
f1_str = re.search(r"F1:=(.*?)\+0:", content, re.DOTALL).group(1)
terms = re.findall(r"\(([-+]?[\d.eE+\-]+)\)\*L\((\d+),u\)", f1_str)
max_n = max(int(n) for _, n in terms)
coeffs = np.zeros(max_n + 1)
for c, n in terms:
coeffs[int(n)] = float(c)
return A, coeffs
The FC file has a subtly different block terminator, a bare 0: on its own line rather than +0:. This is a quirk of how Maple serialises single-variable versus multi-variable expressions, and it is easy to miss: using the wrong terminator would cause the regex to consume the rest of the file. Each FC term has three Laguerre indices instead of one:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
def parse_fc(path):
with open(path) as fh:
content = fh.read()
A = float(re.search(r"\nA:=([\d.eE+\-]+):", content).group(1))
B = float(re.search(r"\nB:=([\d.eE+\-]+):", content).group(1))
C = float(re.search(r"\nC:=([\d.eE+\-]+):", content).group(1))
f_str = re.search(r"F:=\n(.*?)\n0:\n", content, re.DOTALL).group(1)
pattern = r"\(([-+]?[\d.eE+\-]+)\)\*L\((\d+),u\)\*L\((\d+),v\)\*L\((\d+),w\)"
raw = re.findall(pattern, f_str)
coeffs = np.array([float(c) for c, _, _, _ in raw], dtype=np.float64)
iidx = np.array([int(i) for _, i, _, _ in raw], dtype=np.int32)
jidx = np.array([int(j) for _, _, j, _ in raw], dtype=np.int32)
kidx = np.array([int(k) for _, _, _, k in raw], dtype=np.int32)
return A, B, C, coeffs, iidx, jidx, kidx
The three index arrays are stored separately rather than as a list of tuples. This is a deliberate choice that matters significantly in the inner loop, as will become clear shortly.
Evaluating Laguerre polynomials efficiently
The FC wavefunction has 5456 terms. Each term is a product $L_i(u) \cdot L_j(v) \cdot L_k(w)$, and the visualisation grid has $100 \times 100 = 10{,}000$ points, so the inner loop requires on the order of $10^8$ multiplications. The key observation is that many terms share the same index: $L_{17}(u)$ might appear in dozens of terms, and computing it from scratch each time would discard most of the work.
The solution is to build a table $T[n, p] = L_n(x_p)$ for every polynomial order $n$ and every grid point $p$. Built once per coordinate, every term lookup is then a single indexed array read with no further computation.
The Laguerre three-term recurrence:
\[L_0(x) = 1, \quad L_1(x) = 1 - x, \quad L_n(x) = \frac{(2n-1-x)\,L_{n-1}(x) - (n-1)\,L_{n-2}(x)}{n}\]maps directly onto numpy row-by-row operations, no Python loop over the 10,000 grid points, just vector arithmetic on each row:
1
2
3
4
5
6
7
8
def laguerre_table(max_n, x):
T = np.empty((max_n + 1, len(x)), dtype=np.float64)
T[0] = 1.0
if max_n >= 1:
T[1] = 1.0 - x
for n in range(2, max_n + 1):
T[n] = ((2*n - 1 - x) * T[n-1] - (n-1) * T[n-2]) / n
return T
Polynomial order runs along the rows so that T[n] is always a contiguous array, the exact layout needed for the access pattern T_u[iidx[t]] in the inner loop.
The HF density: the trivial case that proves the point
For the HF wavefunction $\Psi_\text{HF}(\mathbf{r}_1, \mathbf{r}_2) = \psi(r_1)\psi(r_2)$, the conditional density is $|\psi(r_2)|^2$, with $r_1$ cancelling exactly. It is computed once and never recomputed regardless of where the slider is placed:
1
2
3
4
5
6
7
8
def hf_density_2d(x_grid, z_grid, A_hf, coeffs_hf):
r2 = np.sqrt(x_grid**2 + z_grid**2).ravel()
u = A_hf * r2
T_u = laguerre_table(len(coeffs_hf) - 1, u)
psi = coeffs_hf @ T_u # dot product over basis functions
psi *= np.exp(-0.5 * u) # attach exponential prefactor
return (4 * np.pi * r2**2 * psi**2).reshape(x_grid.shape)
The factor $4\pi r_2^2$ converts from probability density to radial probability, the likelihood of finding the electron in a thin spherical shell at radius $r_2$, integrated over all angles. Without it, the heatmap is dominated by the very high density near the nucleus where $|\psi|^2$ peaks but the shell volume is tiny. With it, the visual peak sits near $0.5$ bohr, which is where the electron actually spends most of its time, and which matches what chemists call the radial distribution function.
The FC density: coordinates chosen for the physics
The FC case requires more care. The wavefunction is expressed in perimetric coordinates, which are linear combinations of the three inter-particle distances $r_1$, $r_2$, $r_{12}$:
\[u = A(r_{12} + r_2 - r_1), \quad v = B(r_{12} + r_1 - r_2), \quad w = C(r_1 + r_2 - r_{12})\]The defining property of perimetric coordinates is that they are guaranteed non-negative for any physical electron configuration; they satisfy the triangle inequality by construction. This makes them natural arguments for Laguerre polynomials, which are defined on $[0, \infty)$. A Cartesian or spherical expansion would require explicit enforcement of the triangle inequality as a domain constraint; perimetric coordinates encode it automatically.
With electron 1 fixed on the $z$-axis at height $r_1$, the interelectronic distance is $r_{12} = \sqrt{x_2^2 + (z_2 - r_1)^2}$, and the coordinate transformation is:
1
2
3
4
5
6
7
def fc_density_2d(r1_val, x_grid, z_grid, A, B, C, coeffs, iidx, jidx, kidx):
r2 = np.sqrt(x_grid**2 + z_grid**2).ravel()
r12 = np.sqrt(x_grid**2 + (z_grid - r1_val)**2).ravel()
u_arr = A * (r12 + r2 - r1_val)
v_arr = B * (r12 + r1_val - r2)
w_arr = C * (r1_val + r2 - r12)
Floating-point arithmetic can produce small negative values near the triangular boundary, points where $r_{12} \approx r_1 + r_2$ to within rounding error. These are unphysical but we keep them rather than discarding them: zeroing them out preserves the heatmap grid structure, whereas dropping them would introduce visible gaps near the nucleus:
1
2
3
4
valid = (u_arr >= -1e-9) & (v_arr >= -1e-9) & (w_arr >= -1e-9)
u_safe = np.where(valid, np.clip(u_arr, 0.0, None), 0.0)
v_safe = np.where(valid, np.clip(v_arr, 0.0, None), 0.0)
w_safe = np.where(valid, np.clip(w_arr, 0.0, None), 0.0)
Now the Laguerre tables are built once per coordinate, and the wavefunction is evaluated as the sum over all 5456 terms:
1
2
3
4
5
6
7
8
9
10
11
12
T_u = laguerre_table(int(iidx.max()), u_safe)
T_v = laguerre_table(int(jidx.max()), v_safe)
T_w = laguerre_table(int(kidx.max()), w_safe)
psi = np.zeros(len(r2), dtype=np.float64)
for t in range(len(coeffs)):
psi += coeffs[t] * T_u[iidx[t]] * T_v[jidx[t]] * T_w[kidx[t]]
psi *= np.exp(-A * (r1_val + r2))
psi = np.where(valid, psi, 0.0)
return (4 * np.pi * r2**2 * psi**2).reshape(x_grid.shape)
This is where storing iidx, jidx, kidx as separate integer arrays pays off. The expression T_u[iidx[t]] retrieves a contiguous array of 10,000 values in one numpy operation. Had the indices been stored as a list of tuples [(i, j, k), ...], each access would require unpacking at the Python level, roughly an order of magnitude slower over 5456 iterations.
The for t in range(len(coeffs)) loop is still the bottleneck and could be eliminated with Cython or numba. At 8 $r_1$ positions with $N = 10{,}000$ grid points, the total evaluation time is a few seconds on a laptop, acceptable for a once-per-run preprocessing step, though not suitable for on-the-fly browser computation.
Making the browser tool work correctly
Once all densities are computed, the challenge is delivering an interactive figure with no server. The core strategy is to pre-bake every density as a separate Plotly trace, embed all of them in the HTML as JSON, and then have the browser switch between them by toggling which traces are visible. No data travels anywhere at interaction time; the transitions are instantaneous regardless of grid size.
Plotly provides built-in updatemenus buttons intended for exactly this purpose. In practice they are unreliable when the figure contains many pre-baked traces and Plotly’s animation scheduler is active: both buttons end up showing the same density regardless of which is clicked. The root cause is a conflict between the animation state machine and the restyle calls that updatemenus dispatches internally.
The fix is to bypass Plotly’s control widgets entirely. Plain HTML <button> elements call Plotly.restyle() directly from JavaScript onclick handlers, which runs the same underlying API but outside the scheduler:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
function visForHF(ri) {
var v = new Array(1 + 2*NR).fill(false);
v[0] = true; // HF heatmap
v[1 + NR + ri] = true; // e₁ marker
return v;
}
function visForFC(ri) {
var v = new Array(1 + 2*NR).fill(false);
v[1 + ri] = true; // FC heatmap for this r₁
v[1 + NR + ri] = true; // e₁ marker
return v;
}
function applyVis() {
var v = (mode === 'hf') ? visForHF(ri) : visForFC(ri);
Plotly.restyle('plot', {'visible': v});
}
All state lives in two variables: mode (either 'hf' or 'fc') and ri (the current slider index). Both the mode buttons and the slider call applyVis() after updating these variables, so the displayed traces are always consistent. The trace layout is:
| Index | Content |
|---|---|
| 0 | HF heatmap (visible on load) |
| 1 … NR | FC heatmap for each $r_1$ position (all hidden) |
| NR+1 … 2NR | Electron 1 marker for each $r_1$ position (all hidden) |
One deliberate choice: the electron 1 marker is shown in HF mode even though the HF density does not respond to it. Moving the marker while the density stays frozen makes the independence explicit and immediate; which is precisely the point the visualisation is trying to demonstrate. A reader who does not already know that HF is uncorrelated can see it directly.
2D heatmap - drag the slider to move electron 1 along the $z$-axis. The HF density never changes; the FC correlation hole tracks electron 1.
3D surface — density height encodes probability. The correlation hole is a literal dip that follows electron 1 (white ✕). Drag to rotate.
What the Visualisation Shows
When you move the slider on the left, the HF density sits motionless, a spherically symmetric blob centred on the nucleus. No matter where electron 1 is placed, the HF prediction for electron 2 does not shift. This is the mean-field approximation in action: the orbital does not respond to the other electron’s actual position.
On the right, the FC density, for lack of a less poetic term, ‘breathes’. When electron 1 is near the nucleus ($r_1 \approx 0$), both electrons occupy roughly the same region and the density is compressed inward. As electron 1 moves to $r_1 = 1$ bohr (near the density maximum), a depletion, the Coulomb hole, opens up around the position of electron 1: the probability of finding electron 2 there is suppressed because the full wavefunction encodes their mutual repulsion directly. Further out at $r_1 = 2$–$3$ bohr, the density of electron 2 partially follows electron 1 outward, consistent with the nucleus now being less effectively screened and electron 2 drawing closer to it.
The HF density cannot show any of this because its wavefunction, by construction, cannot encode it. The orbital $\psi(r_2)$ does not contain $r_{12}$; it cannot know where electron 1 is.
Why This Matters
The correlation hole is not a curiosity. It is the physical origin of several phenomena that HF gets wrong:
Dispersion forces. Van der Waals attraction between noble gas atoms is entirely a correlation effect. Two helium atoms, at HF level, have zero interaction energy at long range (they are closed-shell and repel). The actual London dispersion force which keeps helium a liquid at 4 K under pressure arises because instantaneous fluctuations in one atom’s density induce a correlated response in the other’s. HF’s mean-field averages these away.
Bond dissociation. As a bond stretches, the electrons in the bonding orbital increasingly localise on opposite atoms. HF forces both electrons to occupy the same spatial orbital, which means it assigns equal probability to the ionic configurations (both electrons on one atom) as to the covalent one (one on each). At dissociation, HF incorrectly predicts a 50/50 ionic-covalent mixture for a homolytic bond. The correlation energy, which enforces that the electrons avoid each other, corrects this.
Reaction barriers. Transition states involve partially formed and broken bonds, precisely the regime where electron positions are most strongly correlated. Reaction barrier heights computed at HF are typically 5-20 kcal/mol too high; DFT and post-HF methods recover most of the missing correlation.
The Density as Primary
There is something which I find clarifying about the density perspective. The electron density $\rho(\mathbf{r})$ can, in principle, be measured directly by X-ray crystallography. It integrates to the number of electrons. It is a function of three variables regardless of how many electrons the system has. Orbitals, by contrast, are gauge-dependent (a unitary rotation of the occupied orbital space leaves all observables unchanged), non-unique, and, in a relativistic treatment, strictly not even well-defined.
The more modern view cemented by the development of density functional theory is that the density is the fundamental quantity. Orbitals are a computational device: a convenient intermediate representation for evaluating the kinetic energy (Kohn-Sham orbitals) or for building systematic approximations to the correlation hole (Hartree-Fock as a starting point for perturbation theory or coupled cluster).
I hope the visualisation above makes this more concrete. The HF orbital is genuinely useful, it gives a compact, chemically interpretable picture, and its squared modulus gives a reasonable first approximation to where the electron actually is. But the shape of that density does not change when the other electron moves, and this is not a technical limitation to be patched later. It is a direct consequence of having chosen a description based on a product of orbitals, which by construction cannot represent correlation.
The fully correlated density does change, because the wavefunction it comes from depends on all three inter-particle distances simultaneously. This is the sense in which it is closer to what nature does: electrons are not independent players reading from individual scripts, but participants in a coupled system where each particle’s behaviour is inseparable from all the others’.
References
- P.-O. Löwdin, “Correlation problem in many-electron quantum mechanics I. Review of different approaches and discussion of some current ideas,” Adv. Chem. Phys. 2, 207 (1959).
- P. Hohenberg, W. Kohn, “Inhomogeneous electron gas,” Phys. Rev. 136, B864 (1964).
- C. L. Pekeris, “Ground state of two-electron atoms,” Phys. Rev. 112, 1649 (1958).
- T. Kato, “On the eigenfunctions of many-particle systems in quantum mechanics,” Commun. Pure Appl. Math. 10, 151 (1957).
- A. L. Baskerville, “Variational and perturbative methods for atomic systems,” PhD thesis, University of Sussex (2019).