A parameter that never finds a value

Fitting a molecular orbital model to two measured ionisation energies gives a whole family of parameter sets that reproduce both exactly. That leaves a question two measurements cannot answer: given two models with different numbers of parameters, both fitted to the same data, is there a computation that says which of them is doing more work?

There is, and it is the obvious one. Show the model more measurements than it has parameters. A model with two parameters fitted to two numbers is a solved equation; the same model fitted to six numbers is a claim that can fail, and what it fails at says what the extra parameters are for.

Done here on six ordinary molecules, the answer comes out in three parts. The third parameter buys almost nothing. It never finds a best value. And the reason for both is visible before any fitting is done at all.

The set, and the property it was chosen for #

Six neutral alternant hydrocarbons: ethene, butadiene, hexatriene, benzene, naphthalene and anthracene. Their first vertical π ionisation energies are measured — 10.51, 9.08, 8.29, 9.24, 8.14 and 7.41 eV — and their highest occupied Hückel eigenvalues are computed here by diagonalising each molecule’s own adjacency matrix.

The set was chosen for a property that has nothing to do with its size. Ethene and benzene have the same eigenvalue. Both come out at exactly $x = 1$, ethene because a two-site system has levels at $\pm 1$ and benzene because a six-ring’s levels are $2, 1, 1, -1, -1, -2$. And butadiene and naphthalene have the same eigenvalue as well, both at $(\sqrt5 - 1)/2 = 0.618034$, for reasons that have nothing in common: one is a four-site chain and the other is a ten-site fused pair.

Neither coincidence is a numerical accident. Both are exact, and they are checked to nine decimal places.

That matters because of what the model is. Every version of it — with overlap or without, at any α, any β, any S — predicts the ionisation energy as a function of the eigenvalue and of nothing else:

$$\text{IE} = -\frac{\alpha + x\beta}{1 + xS}$$

A function gives one output per input. So the model is committed, before any parameter is chosen, to predicting the same ionisation energy for ethene and benzene, and the same one for butadiene and naphthalene. The measurements differ by 1.27 eV and 0.94 eV.

A floor computed without fitting #

That commitment sets a lower bound on the error of every member of the family, and the bound can be worked out without fitting anything.

The best a function can do on a group of inputs that are all the same is to return the group’s mean, and what it is left with is the spread within the group. Doing that arithmetic over the six gives a root-mean-square error of 0.4561 eV, and no choice of α, β and S — nor of any further parameters of the same kind — can go below it.

That number is the yardstick for everything below. A model’s fit is not interesting on its own; what is interesting is where it sits between the floor its form imposes and the scatter of the raw data.

The two-parameter fit #

Least squares on the six, with the overlap set to zero, gives $\alpha = -6.355$ eV and $\beta = -3.551$ eV, and a root-mean-square error of 0.5093 eV.

That is 0.053 eV above the floor. Put the other way round: ninety per cent of this model’s error on this data is irreducible — it belongs to the form of the model rather than to the values of its parameters — and the remaining tenth is what any amount of further fitting has to work with.

The fitted resonance integral, −3.55 eV, is at the upper end of the range the literature quotes for β, which runs from about −2.4 to −3.6 eV depending on which property it was fitted to. That range is the subject of the underdetermined fit, where even a single spectrum leaves it undetermined by a factor of two and a half. So a fitted β landing inside the range is not evidence that anything went right: the range is wide enough to contain most answers, and it is wide for the same reason the fit here is bad.

The other quantities the same fit produces are worth separating out, because they do not all share the problem. The bond orders come out of the eigenvectors and carry no parameter at all; the level spacings in units of β are eigenvalue differences and carry only the scale; and it is the quantities in electronvolts that carry the whole of the parameter’s uncertainty. The ionisation energies fitted here are the third kind, which is why they are the ones that can refuse the model.

What the third parameter does #

Adding the overlap gives a third parameter, and the model is linear in α and β once S is fixed — so the honest thing to do is scan S and solve exactly at each value, rather than run a search that might stop somewhere arbitrary.

The result is the shape of the curve rather than a number.

Within its own range the parameter buys 0.0046 eV, which is nine per cent of the distance to the floor, and it pays for that by moving β from −3.55 to −8.16 eV. A resonance integral of eight electronvolts is outside every published range by a factor of more than two, and it is what the fit wants at an overlap of 0.4 — which is itself at the top of what a neighbouring pair of p orbitals could plausibly have.

The obvious next question is whether the curve turns round somewhere beyond the physical range, since a minimum just outside it would at least be a minimum. It does not.

The residual falls monotonically as S rises, from 0.5093 at zero, and it does not turn round. At the far end of the scan, at $S = 1$ — which is not an overlap between two different orbitals but the overlap of an orbital with itself — the error is 0.5029 eV.

So the third parameter, pushed to a value it cannot have, recovers 12 per cent of the distance between the two-parameter fit and the floor. In absolute terms it improves the fit by 0.0064 eV, which is one part in eighty of the error it started with.

And it moves the other two parameters enormously to do it. β goes from −3.55 eV to −15.03 eV, a factor of 4.2, and α from −6.35 to −4.71. A parameter that has to move by a factor of four to buy a one per cent improvement is not being determined by the data; it is absorbing an error it does not understand.

That is the answer to the question, in a form that generalises: a parameter is doing work when the fit has a minimum in it. A monotone approach to an edge means the parameter has been given a job that belongs to some other term, and the fitting procedure is quietly using it as a fudge. Nothing about the count of parameters says this; only the shape of the response does.

Where the error actually lives #

The floor said that most of the error is in the ties. Looking at the individual residuals says something sharper.

The three open chains — ethene, butadiene, hexatriene — are all predicted to ionise more easily than they do, by 0.605, 0.531 and 0.355 eV. The three cyclic molecules — benzene, naphthalene, anthracene — are all predicted to ionise harder than they do, by 0.665, 0.409 and 0.416 eV. Every chain below zero, every ring above it, and 0.764 eV of clear air between the two groups.

That is not a residual. A residual is what is left when the systematic part has been taken out, and this one sorts the data perfectly by a property the model does not contain. Whether a conjugated system closes on itself appears nowhere in α, β or S, so no value of any of them can respond to it — which is why the fit cannot be improved by fitting harder, and why the third parameter had nothing to find.

What the missing term is has a name and it has appeared twice. A cyclic π system’s closed shell is stabilised relative to the same number of electrons in a chain, which is the shell closure aromaticity is, and the extra stabilisation is precisely what makes the highest occupied level harder to reach than its eigenvalue says. The model has the eigenvalues right — they are symmetry statements and a ring’s levels are its group’s characters — and it has no term for the shell.

Why the fit looks better than it is #

There is one more thing worth saying, because it is what makes this failure survivable in practice.

Fitted over a range of eigenvalues from 0.41 to 1.0, the model’s straight line through six points has a correlation coefficient that would be reported as good, and it produces β within the accepted range. Nothing in the summary statistics announces that half the points are on one side of the line and half on the other by design.

The two things that do announce it are both cheap. One is the tie: two molecules with the same predictor and different measurements can be spotted by sorting the predictor column, and it takes no arithmetic at all. The other is the sign of the residual, which needs a variable to sort by and no fitting.

Both are refusals rather than improvements, and both are available before a line is drawn. The same test applied elsewhere in this collection — a basis set fitted to an energy and tested on a property it was not fitted to — found the same shape of answer: what a model gets wrong is more informative than how well it fits, and the informative part is the pattern of the wrongness.

There is a second reason the measurements and the eigenvalues cannot line up exactly, and it is not aromaticity. The energy needed to remove an electron is not the energy of the level it came out of, because the electrons left behind rearrange — a defect measured in a model that does have repulsion in it, and one that grows with the system. It is another term the fit above has no room for, and it pushes in the same direction as the missing shell closure for the rings.

What a modeller should take from it #

Count the ties before fitting. Sorting the predictor column and looking for repeated values costs nothing and bounds the achievable error exactly. If two rows share a predictor and differ in the measurement, the difference between them is error the model will carry however it is fitted, and quoting a fit without saying so overstates what was achieved.

Watch the shape of the response to each parameter, not the improvement. A parameter with a physical meaning has an optimum inside its physical range; one that improves the fit monotonically towards an edge is being used to absorb something else. This is a stronger diagnostic than counting parameters, because it says which parameter is the fudge.

And look at the sign of the residual against variables the model does not contain. Which numbers carry a frame is the same question asked of a computed quantity rather than a fitted one, and the answer takes the same form: a number whose error responds to a variable is reporting on that variable.

None of the three is statistics, and all three are available on a page of arithmetic.

It is worth being clear about what the third parameter was supposed to mean, because that is what makes its behaviour a diagnosis rather than a nuisance. An overlap between two neighbouring p orbitals raises the antibonding level more than it lowers the bonding one, and a value of S is a statement about the size of that asymmetry. A fit that wants S = 1 is not making that statement: an orbital has unit overlap only with itself.

What is quoted, and what is computed #

The six ionisation energies are quoted, and they are the only quoted numbers in the essay. They are standard vertical values from ultraviolet photoelectron spectroscopy.

The six eigenvalues are computed, by building each molecule’s adjacency matrix and diagonalising it — including anthracene’s, which is built by the same routine that builds any linear acene rather than typed in. Both ties are therefore computed coincidences and are checked as such.

Every fit, every residual, the floor, the scan and the separation are computed. Nothing about the conclusion depends on a fitting library: the model is linear in two of its three parameters, so each point on the scan is a closed-form least-squares solution rather than the output of a search.

What this cannot say #

Koopmans’ theorem is assumed throughout, and it is not exact for anything — as has been measured. Identifying an ionisation energy with an orbital energy introduces an error that grows with the size of the system, and part of the chain-versus-ring separation could be that rather than aromaticity. Distinguishing them would need a model with repulsion in it.

Six molecules is six molecules. The perfect sorting into chains and rings would happen by chance about once in thirty-two times if the signs were random, which is suggestive rather than conclusive; what makes it more than suggestive is the 0.764 eV gap, which is one and a half times the whole root-mean-square error.

Only the first band of each spectrum is used. A photoelectron spectrum has as many π bands as there are occupied π levels, and using all of them would multiply the data by four — but it would also bring in the fact that a band is not one level, and the assignment of the deeper bands is a judgement rather than a measurement.

The vertical ionisation energies are not adiabatic ones, and the difference is a relaxation of the geometry rather than of the electrons. Using the other set would move every number and would not move the ties, since those are equalities between measurements of the same kind.

And the floor is a floor for this form of model only. A model that read anything else about a molecule — a ring count, a total π energy, the stabilisation measured from some reference — is not subject to it.

What was checked #

The set contains at least two pairs of molecules with the same eigenvalue, and the eigenvalues agree to nine decimal places — computed from two matrices that share no structure.

And the members of each pair have different ionisation energies, by more than half an electronvolt. Both halves are needed: the first alone would be a curiosity and the second alone would be a coincidence.

Both fits sit above the floor their own form imposes, which is arithmetic and cannot be otherwise — so it is a check on the fitting rather than a result.

The third parameter buys back less than a quarter of the distance to the floor. Measured: 12 per cent.

The residual is still falling at every step across the range an overlap can occupy, so there is no best overlap inside its own meaning, and the best value on the whole scan sits at an end of the range rather than inside it.

Every chain’s residual is below every ring’s, six of six, with more than a third of an electronvolt between the two groups. Measured: 0.764 eV.

And the refusal is an overlap that makes a denominator vanish. At $S = -1/x$ the model divides by nothing, and the fit stops rather than returning a large number that would look like a very bad fit instead of an invalid one.

The same disease shows up in the other tradition and it is worth naming, because it is the same parameter. A valence bond wavefunction’s ionic weight depends on how the overlap between the structures is divided up, and the number moves by more than the thing it is meant to measure. Both traditions have a parameter the data does not fix, and in both cases the honest response is to say which quantities survive it rather than to quote a fitted value.

A parameter that never settles has a name #

The behaviour measured here — a parameter that improves the fit by nothing and wanders without converging — is a recognised condition with a recognised diagnostic, and naming it is worth a paragraph because the diagnostic is cheaper than the search.

A parameter is unidentifiable when the objective being minimised is flat along its direction: changing it, with the others re-optimised, leaves the residual where it was. A search in such a direction does not fail to converge because it is badly implemented; it fails because there is nothing to converge to, and every point on the flat direction is as good as every other.

The test is one line and needs no optimiser. Fix the parameter at two widely separated values, refit the rest, and compare the residuals. If they agree to the precision of the data, the parameter is not determined and no amount of searching will determine it — which is the result obtained here by watching a search wander, at more cost and with less certainty.

Still open: putting the missing term in #

There is now a family of models the data cannot separate, and a model the data can refuse. What has not been done is the obvious repair: put the missing term in.

A term proportional to the count of rings, or to the π energy per electron, or to a computed shell closure would give the fit something to respond to, and the question is whether one such term absorbs the whole 0.764 eV separation or whether the three rings need three different corrections. That is a computation needing nothing new, and its interest is that a term which works is evidence about what the missing physics is, in a way that a term which merely improves the fit is not.

The nearer question is about the tie itself. Two molecules with the same eigenvalue and different measurements is a specific and reusable instrument — it turns a model into a claim with no fitting anywhere — and there are other predictors of the same shape. A bond order predicts a bond length; a computed splitting predicts a band position; a ring’s response predicts a susceptibility. Finding the ties in each of those, and asking what the within-tie spread is, would say which of those predictors are functions of too little.