The floor was in the bookkeeping

An anomaly that is not the first of a series ran the self-consistent bond-order relaxation up the acene series and found that naphthalene’s anomaly is not the first of anything: the move on the most interior shared bond falls at every step and changes sign. The same calculation can be pointed in another direction. A heteroatom is a perturbation with a location, and the relaxation carries it along the molecule, so the question is how far.

That question has an answer, the answer is a clean exponential over six decades, and the first thing the measurement produced was not it.

The measurement, and what it first said #

Two molecules are relaxed to a fixed point — a chain of twelve fused hexagons with a gap held open by a staggered site energy, and the same chain with one carbon of the end ring shifted by a further $\beta$ — and every bond’s π bond order is differenced. Ring $k$’s response is the mean of that difference over its six bonds.

The profile falls by a factor of about four per ring: $3.8 \times 10^{-2}$ at the ring carrying the heteroatom, $9.2 \times 10^{-3}$ at the next, $1.9 \times 10^{-3}$, $5.7 \times 10^{-4}$, $1.6 \times 10^{-4}$ — and then it stops. From the sixth ring to the twelfth it sits between $4.5$ and $5.8 \times 10^{-5}$, and at the far end it is rising.

A flat tail six orders of magnitude below the head reads as a floor, and a floor in a self-consistent calculation reads as convergence noise. It is not: the fixed point’s residual is $7.3 \times 10^{-14}$, nine orders of magnitude below the floor, and tightening the tolerance changes neither the iteration count nor the answer.

The floor is $2.2 \times 10^{-5}$, $3.4 \times 10^{-5}$, $4.5 \times 10^{-5}$, $5.7 \times 10^{-5}$ and $7.0 \times 10^{-5}$ at feedback strengths of 0.2 to 0.6 — proportional to the coupling, to within a per cent or two. A quantity that scales with a model parameter is the model’s.

Where it comes from #

The standard relaxation sets each bond’s resonance integral from its own bond order,

$$k_{ij} = 1 + \lambda\,(p_{ij} - \bar{p}),$$

and iterates to a fixed point. $\bar p$ is the mean bond order of the uniform calculation, and it is there so that $\lambda = 0$ reproduces the unrelaxed calculation exactly and so that nothing moves in a molecule whose bond orders are all equal. For a molecule considered alone that is exactly right, and the frame that was allowed to relax says so.

It is wrong for a comparison of two molecules. A heteroatom changes the mean bond order — from 0.5263998859 to 0.5282217612 on this chain — so every coupling in the second molecule is shifted against every coupling in the first by $\lambda$ times that difference. That shift has no distance in it. It reaches the far end of the molecule as surely as it reaches the near one, and differencing the two relaxed solutions leaves it under every bond.

The repair is one argument. Measure both runs’ couplings from the same reference, which is the plain molecule’s mean, and the difference between them is a difference between two molecules rather than between two conventions.

The profile then keeps falling: $3.5 \times 10^{-5}$ at the sixth ring, $7.5 \times 10^{-6}$, $2.2 \times 10^{-6}$, $5.9 \times 10^{-7}$, $1.4 \times 10^{-7}$, $2.9 \times 10^{-8}$, $7.9 \times 10^{-9}$ at the twelfth. Six decades of clean exponential where there were two, a coefficient of determination of 0.99966 where it was 0.9669, and a fitted reach of 0.7187 rings where it was 0.9126.

The reported reach was twenty-seven per cent too long, and the error came from a line of bookkeeping rather than from anything about the molecule. A uniform shift in every coupling is very nearly a rescaling of the whole Hamiltonian, which bond orders are blind to — which is why the floor is a hundredth of the shift that causes it, and why nothing before now had ever noticed.

What survives and what does not #

Nothing computed on a single molecule is affected, and that is the whole reason such a defect goes unnoticed. Every earlier comparison was of bonds within one relaxed system, where the reference is common by construction and cancels exactly. The acene result’s headline sequence — the move on the most interior cross bond, falling and changing sign — is a difference between a relaxed and an unrelaxed calculation of the same molecule, so its reference is the same number twice.

The first comparison of two different relaxed molecules is this one, and it is where the convention had to fail.

There is a second reason the floor is easy to accept and hard to see through, and it is worth naming because it will recur. The floor is below everything anybody would look at. A response of $4.5 \times 10^{-5}$ in a bond order is six orders of magnitude under the response at the perturbation, far under any chemical significance, and the natural reaction to a number that small is to call it zero and move on. Doing that gives the right qualitative answer — the heteroatom does not reach the far end — and the wrong length, because the fit was taken through the region where the floor and the decay are comparable.

A number too small to matter can still be large enough to ruin a fit, and the fit is what the calculation was for. The same trap caught a healing length for a different reason, and in both cases the symptom is a fitted length that is too long and looks perfectly reasonable.

A reach that grows, again #

With the defect out of the way the measurement can be trusted, and the first thing it says is the thing already found for a different response.

A bare acene’s reach comes out at 1.1677, 1.3875, 1.4623, 1.5078 and 1.5581 rings at four to twelve rings, growing with no sign of stopping, while its own gap falls from 0.5899 to 0.1102. That is exactly the shape the magnetic reach showed and it has the same cause: the gap that would set the length is closing as the molecule grows, so there is no length to converge to.

Held open by a stagger, the reach settles: 0.7212, 0.7031, 0.7144, 0.7146, 0.7187. Moving by four thousandths across a factor of three in size, which is what a converged quantity looks like.

Two different responses of the same molecules failing in the same way is a stronger statement than either failure alone. It says the diverging length is a property of the acene rather than of what is being measured on it — which is what the magnetic measurement claimed and could not check with one response.

The reach against the gap, and the argument it disagrees with #

The stagger gives a gap that can be set to anything, which turns the reach into a function of one variable.

The slope is −0.688 across a factor of six in the gap, with the profile spanning between three and ten decades at each point and a coefficient of determination above 0.997 everywhere. The argument that a response with a gap decays like an evanescent tail says −1.

That is a discrepancy of the same shape and nearly the same size as the one measured for Peierls chains and then explained away: a fitted −0.60 against an argument that says −1, which turned out to be a profile that was not an exponential at all. This one is not that. Every profile fitted here is a straight line on a logarithmic axis to better than three parts in a thousand over as much as ten decades, so there is no concavity to blame and no window to blame it on.

Two responses, one length #

The last thing to do with a length is to put it beside the other one.

Relaxed, the structural reach is 6 to 16 per cent longer than the magnetic one, and the ratio is not constant, so the two are not one length in two units.

Held rigid, they very nearly are. The rigid structural reach and the magnetic reach agree to 0.7 per cent at the widest gap, to 1.4 per cent at the next, and to within six per cent everywhere; and their exponents against the gap are −0.6798 and −0.6803, which agree to a part in a thousand.

So the picture that comes out is a single object. The π system has one length, a magnetic perturbation and a structural one both see it, and the geometric feedback adds about a twelfth to it — 0.6621 rigid against 0.7187 relaxed on the twelve-ring molecule. The feedback lengthens the reach; it does not create it.

That the two exponents agree to a part in a thousand while both disagree with the argument by a third is the sharpest thing here. Two responses that share nothing but the molecule agree on a number the theory does not give, which makes the number a property of the model rather than of either measurement.

It is also a rare case of two independent measurements of one length to put beside each other. Everything else it measures a length for — a healing length in a distorted chain, a defect state’s decay, the magnetic reach — has exactly one instrument pointed at it, so a systematic error in the instrument is indistinguishable from the physics. Here two instruments that share no arithmetic beyond the molecule agree to a part in a thousand on an exponent, which is the only kind of evidence that can separate the two.

What was computed, and how #

Every bond order is a sum over the occupied orbitals of the product of two coefficients, from an eigensolver checked against closed forms. Each relaxation runs to a fixed point at a residual of $10^{-13}$ or better, reached in sixty-odd sweeps with a mixing of a half, and the residual is returned rather than assumed.

The chain is built by placing hexagon centres $\sqrt{3}$ apart and deduplicating the vertices, with the carbon and bond counts checked against $4m+2$ and $5m+1$. The heteroatom goes on the outermost carbon of the end ring — the site with the fewest neighbours competing for the perturbation, so the decay has the most decades to be measured over — and which site it is is stated rather than left to whichever index the builder assigned.

The fit stops where the profile stops falling. That rule is what the defect above made necessary and it is worth keeping now that the defect is gone: a floor is not a decay, and fitting one as a decay is how a length becomes a weighted average of a real number and a flat tail.

The refusal is a heteroatom of nothing. With the site energy set to zero the two systems are one system, every difference must be at the level of the fixed point’s own residual, and no decay may be fitted to it.

Where the model stops #

There is no chemistry in the heteroatom. It is a shift in one site’s energy, which is what a heteroatom is in this model and is not the same object as an electronegativity difference — nothing here says which atom would do it or by how much, and the reach measured is a property of the model’s response rather than of nitrogen.

The stagger is artificial. It is the right control, being the same molecule with one number changed, and a real molecule whose gap stays open as it lengthens is a different molecule.

And there is no σ framework. A real heteroatom perturbs the σ bonds it sits in as directly as the π system, and this model has no σ bonds at all — so the reach here is the π system’s alone, and it is a lower bound on how far a real substitution is felt rather than an estimate of it.

The generalisation #

A quantity computed as a difference of two self-consistent solutions is only as clean as the reference each of them chose.

The floor here was not noise, was not a numerical tolerance and was not a physical tail. It was a convention — the mean coupling each relaxation measures for itself — appearing in both solutions with different values and failing to cancel in the subtraction. It sat four orders of magnitude above the effect being measured, and it was invisible from inside either solution, because within one system a reference that shifts every bond by the same amount changes nothing at all.

That is the general shape, and it is not special to bond orders. Anything self-consistent chooses something: a chemical potential, a normalisation, a mean field, an origin. Choosing it per system is almost always right, because the choice is what makes the solution well conditioned. Differencing two such solutions then silently differences the two choices as well, and the residue is a constant that does not decay, does not depend on distance, and looks exactly like the floor of a real measurement.

The test is cheap and it is the one worth running on every profile: set the perturbation to zero and demand that the difference vanish. A convention that fails to cancel cannot pass that, because with no perturbation the two systems are one system and every difference must be at the fixed point’s own residual. A defect like the one above can survive indefinitely, because no ordinary check compares a number to what it should be when nothing was done to it.

Who found it, and when #

Self-consistent bond-order–bond-length relaxation is Coulson and Golebiewski’s, from 1961, and the linear relation used here is theirs. The observation that a perturbation in a gapped one-dimensional system decays exponentially over a length set by the gap is older still and is the same statement as the exponential falloff of a density matrix, which Kohn made precise in 1959 and which has since been sharpened considerably — including the result that the decay in a one-dimensional insulator is controlled by a complex band structure rather than by the gap alone, which is one place a two-thirds power could come from.

The defect above is a bookkeeping one, and it belongs to an implementation rather than to the literature. What is worth carrying from it is the shape rather than the line: a convention that cancels within a system does not cancel between two of them, and a quantity computed as a difference of two self-consistent solutions inherits every reference either of them chose.

Still open: the exponent against a closed form, and a central heteroatom #

The obvious open question is the exponent. Two responses agreeing on −0.68 while an argument says −1 is either a fact about this model or a fact about the fits, and the way to tell is a system where the answer is known: a one-dimensional chain with a staggered site energy has a closed-form band structure, its complex band at the gap can be written down, and the decay it predicts can be compared against a measurement made exactly the way these were. If the closed form gives −1 the fits are wrong; if it gives two thirds, the argument is.

The nearer question is where the heteroatom sits. Everything here puts it on the end ring, because that gives the longest run to measure over, and the response of a molecule to a perturbation at its end is not the response to one in its middle — the whole of the ring-current scatter is an end effect. Running the same measurement with the heteroatom in the central ring costs one line and gives a profile in two directions at once, which would say whether the reach measured here is the molecule’s or the end’s.