Two kinds of correlation, and only one is small

A mean field is not one thing, and neither is what it misses.

The correlation energy is defined as a subtraction: the exact energy minus the best single-determinant energy. It is the quantity every method beyond Hartree–Fock exists to recover, and it is quoted in millihartrees, in kilocalories per mole, and as a percentage recovered.

The definition hides that the thing being subtracted is not one thing. Practitioners separate it into dynamic and static — sometimes non-dynamic, sometimes strong — and the separation is usually explained by example rather than by a number. In the two-site model both kinds are present and the model is small enough to separate them exactly.

The dimer, written in orbitals #

The smallest many-electron calculation is two electrons on two sites with an on-site repulsion $U$ and a hopping $t$. Written in the site basis its singlet sector has three configurations; written in the bonding and antibonding orbitals it has two that matter, and the two-by-two matrix is

$$\begin{pmatrix} -2t + U/2 & U/2 \\ U/2 & 2t + U/2 \end{pmatrix}.$$

The diagonal entries are both electrons in the bonding orbital and both in the antibonding one; the off-diagonal is the repulsion coupling them. Its lower eigenvalue is $(U - \sqrt{U^2 + 16t^2})/2$, which is the closed form for the dimer, and its eigenvector says how much of each configuration the ground state is.

The restricted single determinant is the first configuration alone, at $-2t + U/2$.

Occupation numbers, and what a single determinant is #

The useful diagnostic is not an energy at all. The natural orbitals are the eigenvectors of the one-particle density matrix and their eigenvalues are the occupations, and a single determinant is precisely a state whose occupations are integers — two and zero for a closed-shell pair.

So how far the occupations are from integers is how far the state is from being describable by one determinant, and it is a statement about the state rather than about an energy.

At $U = 0$: 2.0000 and 0.0000, entropy zero, correlation energy zero. Everything agrees that there is nothing to recover.

At $U = 4t$: 1.707 and 0.293. At $U = 12t$: 1.316 and 0.684. At $U = 120t$: 1.033 and 0.967, and the entropy is 0.6926 against $\ln 2 = 0.6931$.

The convergence is monotone at every step, so the failure of the single determinant is a smooth, one-way process with no threshold in it. That is worth knowing, because dynamic and static are usually presented as two regimes with a boundary.

The entropy, and what it counts #

The occupations have a natural summary that is worth naming because it is what a modern calculation reports: the entropy of the occupation numbers, $-\sum_i p_i \ln p_i$ with $p_i = n_i/2$.

It is zero for a single determinant, and it rises to $\ln 2$ — one bit — for the strongly correlated limit of a two-electron singlet. Its value is a count, in the information-theoretic sense, of how many determinants the state really needs, and unlike the correlation energy it is bounded.

That boundedness is the point. The energy error grows without limit and the entropy saturates, so the two quantities disagree not merely in scale but in whether they have a maximum. A method judged on the energy has an unbounded problem to solve; a method judged on the description has a bounded one, and knowing which is which decides whether adding more terms to a series can ever be enough.

The number the entropy tends to is exactly $\ln 2$ and not approximately: at $U = 120t$ it is 0.6926 against 0.6931, and the shortfall is the residual occupation of the bonding orbital. Two determinants of equal weight is the most entangled a two-electron singlet can be, and it is what a molecule at its dissociation limit is — which is why molecular orbital theory dissociates a hydrogen molecule wrongly, and why the failure there is not a small error but a wrong state.

The energy behaves differently #

Now the correlation energy over the same range, and it does not behave the same way at all.

$E_{\text{corr}} = 2t - \tfrac{1}{2}\sqrt{U^2 + 16t^2}$, which is $-0.0156$ at $U = 0.5t$, $-0.828$ at $U = 4t$, $-4.32$ at $U = 12t$, and grows without limit. As a fraction of the exact energy it is worse than useless, because the exact energy tends to zero as the repulsion grows: the fraction reaches eighteen thousand per cent at $U = 40t$.

A natural first guess is that the correlation energy goes through a maximum as a fraction of the total and comes back down. It does not, and the eighteen thousand per cent refutes it. What separates the regimes is not the size of the error but its power.

At weak repulsion the correlation energy is second order in $U$: a perturbative correction to a good reference, which is what Møller–Plesset theory is and why it works when it works. This is the regime the correlation hole describes: electrons keeping out of each other’s way by a small, smooth adjustment to an otherwise good description. At strong repulsion it is first order, which no perturbation series about a single determinant produces at any order.

That is the honest content of the dynamic/static distinction. Dynamic correlation is a correction; static correlation is a failure of the thing being corrected.

An identity nobody was looking for #

The two curves in that figure are the same curve.

$$\frac{d \ln |E_{\text{corr}}|}{d \ln U} = n_{\text{bonding}}$$

at every repulsion, to nine decimal places. The measured slope at $U = 4t$ is 1.707106781 and the occupation is 1.707106781. At $U = 12t$: 1.316227766 against 1.316227766.

It falls out of the closed forms. With $R = \sqrt{U^2 + 16t^2}$ the slope is $U^2 / R(R - 4t)$ and the occupation is $2U^2/(U^2 + (R-4t)^2)$, and the denominator of the second simplifies to $2R(R - 4t)$. Two routes with nothing in common — a numerical derivative of an energy, and an eigenvector of a two-by-two matrix — landing on the same number.

The agreement is limited by the difference quotient rather than by the identity: halving the step halves the disagreement, which is what a finite difference does and what an approximate relation would not do.

So “second order in the interaction” and “one determinant” are not two facts about weak correlation. They are one fact written twice. The exponent is the occupation number, and the smooth fall from two to one is the smooth departure of the state from a single determinant.

The threshold that is not a threshold #

There is a second diagnostic already in this collection, and the identity explains why it agrees.

A mean field can be made to lower its energy by breaking spin symmetry — allowing the two spins to occupy different orbitals — and the point at which that first pays is a well-known instability. For the four-site chain it happens at $U = 2t$.

By $U = 2t$ the antibonding natural orbital already carries a tenth of an electron. So the instability of the mean field and the departure of the occupations from integers are pointing at the same thing, from two sides: the first is the mean field noticing it cannot describe the state, the second is the state saying so.

That is a case worth keeping in view. The broken-symmetry solution is closer in energy and further from the truth in description, which is the sharpest possible statement that an energy is not a measure of how good a wavefunction is.

The measure that cannot tell an attraction from a repulsion #

The refusal written for this calculation turned into the sharpest thing in it.

The first version required that an attraction — a negative $U$, driving both electrons onto one site — make the state more like a single determinant. It was refused: the antibonding occupation at $U = -4t$ is exactly what it is at $U = +4t$.

Both the correlation energy and the natural occupations are even functions of $U$. Reversing the sign of the interaction only reverses the sign of the off-diagonal element in that two-by-two matrix, and no scalar built from magnitudes can see it.

The two states could hardly be more different. At $U = +4t$ the ground state is covalent — one electron on each site, ionic weight 0.146. At $U = -4t$ it is ionic — both electrons on one site or the other, weight 0.854. The two weights sum to exactly one, reflections of each other about a half.

Every measure of correlation gives them the same number. The thing that tells them apart is the sign of the mixing coefficient, which is what every scalar diagnostic throws away.

Two pictures of the same state #

There is a third reading of the same numbers and it belongs to a different vocabulary.

Written in the site basis the ground state is a mixture of ionic configurations — both electrons on one atom — and a covalent one, and the mixture is what valence bond theory has always described. Written in the orbital basis it is a mixture of a bonding pair and an antibonding pair, which is what configuration interaction describes.

They are the same state, and the exact state lies in the plane spanned by the two descriptions at every repulsion. The ionic weight above — 0.500 at $U = 0$, falling to 0.146 at $U = 4t$ and 0.026 at $U = 12t$ — is the valence bond reading of exactly the numbers the natural occupations give.

So static correlation in the orbital language is ionic-covalent resonance in the valence bond language, and neither is more fundamental. What the orbital language has is the occupations, which make the failure of a single determinant a number; what the valence bond language has is the sign, which the occupations throw away.

What this cannot say #

Two sites. The dimer has one bonding and one antibonding orbital, so its correlation problem has exactly one parameter’s worth of structure. A real molecule has many partially occupied natural orbitals and the distinction between the two kinds is a distinction between sets of them, which is genuinely harder and genuinely fuzzier.

On-site repulsion only. The Hubbard interaction acts only when two electrons are on the same site, so the long-range part of the Coulomb interaction — which is most of what dynamic correlation is in a real molecule — has no analogue here.

No basis set. Every number above is exact within the model — there is no contracted set of primitives and no fitting of any kind, so none of the questions about what a finite basis does to a correlation energy arises. In a real calculation the correlation energy converges towards its basis-set limit far more slowly than the Hartree–Fock energy does, and separating that convergence from the physics is most of the practical difficulty.

The threshold practitioners already use, and what the identity makes of it #

The occupation number is not an exotic diagnostic invented for this essay. It is close to the quantity working quantum chemists check before trusting a single-reference calculation, and the identity above gives their rule of thumb a mechanism it has never had.

The practical tests all measure the same thing by different routes. The weight of the leading determinant in a correlated wavefunction: below about 0.9 and the calculation is in trouble. The largest single-excitation amplitude, a number produced automatically by every coupled-cluster program: above about 0.02 and the result is flagged. The natural orbital occupations themselves: a bonding orbital below roughly 1.98, or an antibonding one above 0.02, and the same warning applies.

Every one of those is a statement about how far the state is from a single determinant, and not one of them is an energy. That is the essay’s point arriving in the form of a habit — a field that has learned by experience to stop trusting its energies at a threshold set on a quantity that is not an energy.

What the identity adds is why the threshold is where it is. If the local power of the correlation energy in the interaction equals the bonding orbital’s occupation, then a threshold on the occupation is the same threshold on the exponent. An occupation of 1.98 is an exponent of 1.98 — a correlation energy still very nearly quadratic in the interaction, which is what a perturbative treatment assumes and what makes recovering a percentage of it meaningful. An occupation of 1.5 is an exponent of 1.5, and a perturbation series with a non-integer effective power is not a perturbation series that will converge in the way its user expects.

So the rule that a single-reference method may be trusted while the leading occupation stays near two is not a piece of accumulated superstition. In this model it is exactly the statement that the correlation energy is still second order in the interaction, which is the condition under which every step of the standard argument — a small correction to a good reference, expanded in powers, truncated — is a valid step.

Two qualifications keep that honest.

The identity is the dimer’s. Whether the exponent tracks a single occupation in a larger system is not established here, and the reason for expecting a weighted average rather than one number is obvious. The threshold’s shape — a quantity that measures distance from a determinant, used to decide whether an energy expansion is legitimate — is what carries over, rather than the equality.

And the thresholds are conventions with an empirical origin. Nobody derived 0.02; it was arrived at by finding out which calculations went wrong. What the identity supplies is not a better number but an account of what kind of number it is: a bound on how far the correlation energy has drifted from the power it needs to have for the method being used to be the right method.

What a percentage of correlation energy is worth #

Quantum-chemical methods are routinely compared by the fraction of the correlation energy they recover, and the identity above says what that fraction is measuring.

Near the weakly correlated end, where the exponent is close to two, the correlation energy is a well-behaved perturbative quantity and recovering ninety-five per cent of it means the wavefunction is close to the reference with a small correction. The percentage is doing what it appears to.

Away from that end it is not. The exponent has fallen, the energy has become first order in the interaction, and the reference is qualitatively wrong — but the correlation energy has also become large, so recovering ninety-five per cent of it is easier in relative terms and less meaningful. A method can score well on the percentage precisely where the state it is describing has stopped resembling the reference.

That is not a hypothetical failure mode; it is the reason multireference methods exist, and the reason a single number quoted for a method’s accuracy is always accompanied by a class of systems it applies to. The occupations are the honest diagnostic because they are a property of the answer rather than of the distance from a starting point.

Still open: more orbitals, and size consistency #

The identity above is a statement about a two-level system and it invites the obvious question: does something like it survive with more orbitals? The natural generalisation would relate the scaling of the correlation energy to a weighted average of the natural occupations, and there is no reason to expect it to be exact — but the dimer’s version is exact, and exact relations in small models are the ones worth trying to break.

The other direction is the one size consistency points at. A method exact for one dimer can fail for three — size consistency — and the failure is worst where the exact answer is easiest — the truncation is 2.54 short for three dimers at no repulsion at all, against 0.49 at a repulsion of four. Putting that beside the identity here suggests a sharper statement of what a truncated method is doing wrong: not recovering too little correlation energy, but recovering it with the wrong scaling, so that the error per unit grows with the number of units. Testing that needs the size-consistency calculation and the occupation analysis run together, on the same systems, which is one calculation and has not been done.