The six-membered aromatic ring is ubiquitous in biology. It is a component of DNA and proteins, as well as woody biomass and petroleum1. The aromatic ring has been observed in interstellar space2,3, and aromatic structures, in general, are thought to pervade the interstellar medium4,5. The conjugated six-membered ring is also the building block of graphene6, a material with astonishing electronic properties.

Benzene is the archetypal aromatic molecule, displaying exceptional chemical stability as compared with other unsaturated hydrocarbons. It was first reported in 1825 by Faraday7, and its electronic structure has been discussed for over a century. Many structural formulae have been proposed, with the most satisfying single structure being that with alternating single and double bonds, proposed by Kekulé8,9. However, this structure does not capture the sixfold symmetry required to satisfy experimental evidence. To account for this, Kekulé proposed an oscillation between these structures, despite knowing nothing of electrons or quantum mechanics.

The discovery of the electron in 1897 led Lewis to propose the two-electron chemical bond in 1902, although this idea was not published until 191610. Following the introduction of quantum mechanics, Pauling showed that both Kekulé structures, which have the same energy, would resonate to bring about a lower energy form11. Hückel advanced a different approach, in which electrons would occupy molecular orbitals (MOs), delocalised over multiple atomic centres12.

As chemical electronic structure theory advanced in the ensuing decades, there emerged a competition between the MO approach, which offers computational simplicity, and the valence bond approach, which accords with, and engenders chemical intuition13. It has been suggested that banana (bent) double bonds are a superior description to the \(\sigma -\pi\) structure which is usually offered in chemical textbooks14, and this has been taken as evidence that the Kekulé description is more authentic than the delocalised, MO picture15,16,17,18. Empedocles and Linnett went further, and suggested a decoupling of the two spin sets of electrons such that the two sets of electron spins occupy alternate Kekulé structures19,20. Nevertheless, the delocalised canonical MOs of benzene are usually invoked to explain spectroscopic phenomena.

It has long been known that, after sufficient computational effort, the MO and VB approaches converge to the same description of electronic structure, and are thus equally valid as computational formalisms21. The MO approach has won out, on the whole, for reasons of computational simplicity, and it has been argued that this is at the expense of chemical insight13.

But, interpreting the electronic structure of a molecule in terms of the one-electron spin orbitals generated by MO theory has several major and irreparable flaws. Since the seminal work of Slater22, it has been known that a spatial wavefunction must exhibit anti-symmetry with respect to the exchange of like-spin electrons: The sign of the spatial part of the wavefunction must change when two electrons of the same spin swap positions. This fundamental property of half-integer spin particles is the origin of the Pauli exclusion principle, and led to the development of Slater determinants to describe simple wavefunctions. The MOs are non-unique with respect to the anti-symmetrised wavefunction since the determinant is invariant upon any unitary transform of the occupied orbitals. Furthermore, even if one were to insist upon employing the canonical MOs which result from correctly anti-symmetrised Hartree–Fock theory, interpretations that view these orbitals individually ignore anti-symmetry: The Hartree product is not a correct wavefunction22.

However, we may inspect the 3\(N\)-dimensional electronic wavefunction (\(3N=126\) in the case of benzene) that results from any theoretical framework (including MO theory) to regain chemical insight. This should be done without imposing (possibly erroneous) chemical intuition on the problem at hand. The electrons being described by the electronic wavefunction are fundamentally indistinguishable particles. That means that in the 3\(N\)-dimensional space of the wavefunction there are regions that are equivalent, related through the permutation of electrons. Because these regions are equivalent and together span the space of the wavefunction, they are analogous to tiles. All of the information of the wavefunction is contained in just one of the permutable tiles23,24.

Recently, we reported a method to identify and visualise wavefunction tiles. The procedure, known as dynamic Voronoi Metropolis sampling (DVMS), uses an iterative algorithm to identify regions in 3\(N\)-dimensional space delineated by nodes, in the case of a boundary between differently signed tiles, and a Voronoi diagram in the case of like-signed tiles25,26,27,28,29. A Voronoi diagram is a partitioning of space into regions closest to a particular site (of several). Here, the Voronoi decomposition of the 3\(N\)-dimensional space of the wavefunction is according to \({N}_{\alpha }!{N}_{\beta }!\) sites generated by permuting the labels on the electrons of the same spin, where \({N}_{\alpha }\) and \({N}_{\beta }\) are respectively the number of \(\alpha\) and \(\beta\) electrons. The DVMS algorithm is designed to find a self-consistent position \(\overline{{\bf{x}}}\) for one permutable site (and thus of all the permutations defining the Voronoi decomposition) that is the centroid of a tile \(R\), such that \(\overline{{\bf{x}}}={\int }_{R}{\bf{x}}{\Psi }^{2}d{\bf{x}}\). Self-consistency is required as \(\overline{{\bf{x}}}\) depends on \(R\), while simultaneously \(R\) is defined by \(\overline{{\bf{x}}}\).

The DVMS procedure reproduces motifs such as core electrons, single bonds, multiple bonds (of the banana type) and lone pairs25. The method takes as its sole input the correctly anti-symmetrised wavefunction, and is thus agnostic to the theoretical approach. As such, multiconfigurational wavefunctions can be investigated to inspect the effects of electron correlation: A 2-configuration wavefunction for C\({}_{2}\) was found to exhibit a triple bond and singlet-coupled biradical electrons occupying positions on the ends of the molecule, in accord with Shaik et al.’s quadruple bond25,30. The DVMS procedure has also been shown to reproduce the curly arrow electron movements along a chemical reaction path26, and connect electronic spectroscopy with the concept of vibrational normal modes27.

In this work, we revisit the electronic structure of benzene using the DVMS procedure. It is found that a single-determinant wavefunction reproduces Kekulé structures, but with the two spins sets equivocal to the other’s occupation of one Kekulé structure or the other. As the number of excited configurations is increased, the wavefunction evolves to distinctly prefer a relative disposition of the spins such that their Kekulé structures are anti-correlated.