Noncoplanar and chiral spin states on the way towards Néel ordering in fullerene Heisenberg models

,


I. INTRODUCTION
Antiferromagnets on infinite bipartite lattices generally show Néel ordering in dimensions greater than one.This can be detected through the staggered magnetisation of the ground state, which spontaneously breaks SU(2) spin-rotation symmetry.In the spectrum, the Goldstone mode corresponding to this symmetry breaking gives rise to a proliferation of gapless states with a range of spin quantum numbers, known as the Anderson tower of states [1], as well as a branch of gapless spin-wave excitations.Such Néel ordering of the ground state can be proven rigorously for the Heisenberg model on the three-dimensional cubic lattice for spin  ≥ 1/2 [2,3] as well as on the two-dimensional square lattice for  ≥ 1 [4].
For two-dimensional spin-1/2 systems, however, and especially to study ordering in frustrated magnets, one has to rely on numerical studies that are almost always performed on finite patches of the lattice.Spontaneous symmetry breaking never occurs in these finite systems.Nevertheless, ordering tendencies in the thermodynamic limit are already indicated by sharp Bragg peaks (with intensity proportional to system size) in the static correlation functions, as well as incipient Anderson towers of states in the spectrum, at energies well below those of quasiparticle excitations.Even beyond conventional ordering, the symmetry-resolved low-energy spectrum is an invaluable diagnostic of phases of matter, both computationally [5][6][7][8][9] and experimentally [10].
An interesting alternative to the paradigm above is considering strongly correlated systems on highly symmetric molecular geometries, which also exhibit a wide range of unusual quantum magnetic properties, such as magnetization jumps and plateaus, or the proliferation of lowest singlet (rather than magnetic) excitations [11][12][13].A case in point are fullerene structures [14], made up of pentagonal and hexagonal faces: while there is no limit on the number of hexagons, Euler's formula implies that they all have 12 pentagonal faces.This allows interpolating from the limit of extreme frustration (C 20 , a dodecahedron with pentagonal faces only) to large molecules that resemble the bipartite honeycomb lattice with a vanishing fraction of frustrated defects.
As a starting point to understanding strong-correlation ef-fects in carbon fullerenes, exact-diagonalisation (ED) and quantum-Monte-Carlo studies were performed on the C 20 Hubbard model [15].These found that the spin-triplet ground state of the weakly-interacting Hückel limit switches to a nondegenerate singlet as interactions are made stronger, consistent with the Heisenberg model on the same geometry [16].This shows that the Heisenberg limit provides useful information about the physically more relevant [17][18][19] intermediate- Hubbard model, which would pose considerably greater computational challenge.The lowest-lying excitations of the C 20 Heisenberg model are also singlets, including one belonging to a five-dimensional irreducible representation (irrep) of the icosahedral point group: the absence of low-energy triplets is incompatible with incipient magnetic ordering, as one would expect for such a highly frustrated molecule.While ED results have been obtained for the Heisenberg model on fullerene allotropes up to C 36 [20][21][22], their high degree of frustration and varying degrees of symmetry obstructs the emergence of any systematic trend.Much attention has also been devoted to the famous C 60 buckminsterfullerene geometry [Fig.1(b)], motivated by interest in chemical, nanotechnology, and quantum-computing [25] applications, as well as superconductivity observed in alkalimetal-doped fullerene lattices [26].Early studies of the Hubbard model [18,19] argued in favour of a nondegenerate singlet ground state at strong and intermediate interaction on this geometry as well.More recently, this has been corroborated by DMRG studies of the Heisenberg model by Rausch et al. [27], which found the lowest few eigenstates of the  = 0, 1, 2 spin sectors.In particular, these authors uncovered that the lowestlying triplet excitation is threefold degenerate and breaks cubic rotation symmetry.However, a detailed analysis of the spatial symmetries of the low-lying spectrum, desirable for understanding the low-energy physics and any symmetry-breaking tendencies of large fullerenes, requires complementary numerical approaches since, in general, tensor-network methods struggle to resolve spatial symmetries.
Here, we address this demand using a variational Monte-Carlo approach based on group-convolutional neural-network (GCNN) wave functions [28,29].These allow us to resolve the lowest-lying states in every spatial symmetry sector with Figures made using VESTA [23]; coordinates of C 32 taken from Ref. [24].
modest computational resources, and thus reconstruct much of the low-energy spectrum.In particular, we study the spin-1/2 nearest-neighbour Heisenberg model on the highest-symmetry allotropes of C 32 , C 60 , and C 80 , shown in Fig. 1.In what follows, we use  = 1 as the unit of energy and consider the zero-temperature limit.The smallest molecule allows us to benchmark the method against ED: despite the high degree of frustration, we obtain variational energies very close to the lowest exact ones in every symmetry sector considered.Likewise, our variational energies for C 60 match those obtained from DMRG; however, we obtain dozens of additional energies and wave functions across all icosahedral symmetry sectors.
Most importantly, we are able to account for much of this low-lying spectrum by adapting arguments on towers of states and the ground states of classical ( → ∞) Heisenberg models developed for lattice models.In particular, we find that the lowest-lying  = 0, 1, 2 states are captured by a triplet of lowenergy  = 1 modes, which play the role of Goldstone modes (at gapless points of the magnon dispersion relation) of an incipient noncoplanar order.This order can be understood as the result of the competition between incipient Néel ordering on the hexagons and frustration introduced by the pentagons: this is highlighted by the wave functions of the Goldstone modes, which show Néel-like alternating signs on maximal bipartite segments of the C 60 geometry.We provide a recipe, analogous to the Luttinger-Tisza method for lattice magnets, to predict the symmetry properties of these modes, which matches the numerical calculations perfectly.
We also perform the same analysis for the next smallest icosahedrally symmetric fullerene, C 80 .Remarkably, its lowlying spectrum comes in nearly degenerate pairs of states, which only differ in their parity under spatial inversion.We again account for this behaviour in terms of an incipient symmetry-breaking order.The large- Heisenberg model on this geometry has a chiral ground state (that is, it breaks spatial inversion and time-reversal, but not their product), for which tower-of-states analysis predicts such a degeneracy.We also construct an explicit chiral operator in terms of the Goldstonemode operators of the incipient order to relate the pairs of states to one another.Detecting such a chiral ordering in C 80 would be an interesting target of future computational (e.g., DMRG) and experimental studies.
The rest of the paper is organised as follows.In Sec.II, we generalise methods to detect incipient ordering in finite systems to molecules without translation symmetry.In Sec.III, we describe our GCNN ansatz and its optimisation protocol in detail, and benchmark it against ED on the C 32 molecule in Sec.III D. We detail our numerical studies of C 60 and C 80 in Secs.IV and V, respectively.Conjectures on the low-energy spectra of larger fullerenes, based on semiclassical arguments, are presented in Sec.VI.Perspectives and conclusions are given in Sec.VII.An appendix on subspace projection of high-dimensional irreps (Appendix A) and tables of exact and variational energies (Appendix B) complete the paper.

II. INCIPIENT ORDERING IN MOLECULAR MAGNETS
On an infinite two-dimensional lattice, the ground state of a magnetic Hamiltonian may break spin-rotation symmetry.This is indicated by the emergence of Bragg peaks, divergences of the reciprocal-space correlation function ⟨ ì (− ì  0 ) • ì ( ì  0 )⟩ at some wave vector ì  0 , as well as a gapless branch of Goldstone modes (magnons) corresponding to rotating the order parameter direction.The magnons become gapless at the Bragg peak position; in a primitive lattice, repeated application of the magnon creation operator to the ground state creates a sequence of zero-gap states with different spin quantum numbers, known as the (Anderson) tower of states [1,5].Symmetry quantum numbers of the states in this tower can be derived from the above construction, or by decomposing symmetry-broken classical ground states into irreps of the full symmetry group  spatial × SU(2) [5,12,30].While no true phase transition is possible on a finite system, symmetry-breaking tendencies can readily be established from numerical simulations of finite lattices, either from the finitesize scaling of reciprocal-space correlators, or from the energy spectrum, which contains a distinct set of low-lying excitations with symmetry quantum numbers consistent with the Anderson tower of states [5].
It is reasonable to expect similar precursors to ordering on large fullerene geometries: Even though these are always frustrated due to having 12 pentagonal faces, in the largemolecule limit, almost every face is hexagonal, so we can regard the structure as a large honeycomb lattice with a finite number of defects.Therefore, physical properties away from these defects ought to approach those of the honeycomb lattice, which sustains Néel order [31,32].
Since the fullerene geometry has no translational symmetry, we cannot directly probe such incipient ordering in reciprocal space.Bragg peaks, however, can be extracted in real space as well, as the dominant eigenvector of the correlation matrix , with a diverging eigenvalue corresponding to the order parameter [33].Likewise, the leading eigenvector of    on the fullerene geometry can be thought of as a real-space Bragg-peak "wave function"   .This wave function can be used to construct the Goldstone-mode operator repeated application of which creates an ansatz "tower of states" that can be compared to the low-lying eigenstates of the full Hamiltonian.Just as in the case of lattice systems, the (point-group) symmetry quantum numbers of this tower of states can be deduced either from the repeated application of the bosonic [34] operators Ŝ , or from decomposing a symmetry-broken classical ground state into irreps of  point × SU(2), the latter of which can be made controlled in the large- limit, where such symmetry-breaking ground states may form even for finite systems [12].A further analogy with lattice magnets allows us to predict these symmetry quantum numbers directly from the Hamiltonian, without computing the correlation matrix    of the quantum many-body ground state, similar to the Luttinger-Tisza method for lattice magnets [35][36][37].In the large- limit underlying the above arguments, the Hamiltonian (1) on a lattice can be Fourier transformed, without having to worry about complicated commutation relations between the ì ().The minimum of  () predicts the position of Bragg peaks, subject to compatibility with the unit-length constraint on spins in real space, which may also determine whether the order is collinear, coplanar, or noncoplanar [37].Likewise, the lowest-energy eigenvector of the Hamiltonian matrix (in our case, the adjacency matrix of the fullerene graph) is expected to recover the Bragg-peak wave function   .

III. GROUP-CONVOLUTIONAL NEURAL-NETWORK STATES
In the following section, we describe our numerical method to obtain the low-energy spectrum.Group-convolutional neural networks (GCNNs) [28,29,38], which play a central role in our approach, are discussed in Sec.III A; for a general discussion of neural-network-based variational Monte Carlo, we refer the reader to Refs.[39,40].In Sec.III B and Appendix A, we explain how the symmetry of a GCNN wave function can be constrained beyond the (multi-dimensional) irreps of the point group, which we found to substantially improve our results.Specific details of the GCNN architecture and other hyperparameters are given in Sec.III C. Finally, benchmarks against ED on a C 32 allotrope are presented in Sec.III D.

A. Ansatz
Space-group symmetries (that is, ones that map computational basis states onto one another) can be imposed on any variational ansatz  0 using the projection formula [41] where the ĝ are the elements of the space group  and the   are their characters in a given   -dimensional irrep of .Here  stands for a spin configurationin the computational   basis:  = (  1 ,   2 . . .   ) for  spins-1/2 with    ∈ {↑, ↓} ≡ {±1}.Applying (5b) directly, however, requires evaluating the ansatz  0 many times, which may be prohibitively computationally expensive.
Instead, we use group-convolutional neural networks (GC-NNs) [28,29,38], a generalisation of the well-known convolutional neural networks (CNNs) to non-translational symmetries, which are able to efficiently generate all symmetryrelated evaluations of a neural-network ansatz as hidden layers indexed by the symmetry elements.Our feed-forward GCNNs start with an embedding layer which converts the input spin configuration  into such a hidden layer, ℎ (1) .Then, further group-valued hidden layers are generated by alternating nonlinearities and equivariant linear layers of the form The trainable variables of the ansatz are the kernel entries  (ì ) and  ().One can show [42] that acting with a symmetry  [27] using DMRG (cf.Fig. 6).
element k on the input spin configuration  permutes the labels of all subsequent layers as Therefore, we can regard the entries of the last layer as the amplitudes of a neural-quantum-state (NQS) ansatz ℎ () 0 () evaluated for all spin configurations related to  by spacegroup symmetry: so a symmetric ansatz is obtained by combining all entries of the last layer according to the projection formula (5).
In addition to using spatial symmetries of the molecules and the conservation of the magnetization   =  =1    , the parity symmetry P =  =1    , an element of the SU(2) spinrotation group, can be implemented in the   computational basis by flipping the sign of all    .We therefore imposed eigenvalues of  = ±1 on our ansätze in addition to the spacegroup irreps.Sampling in the   = 0 magnetisation sector, this allows us to distinguish between states with even ( = +1) and odd ( = −1) total-spin quantum number.We also performed simulations in the   = 2 sector, which isolate total-spin quantum numbers  ≥ 2.

B. Irrep subspace projection
In our numerical experiments, we found that training ansätze projected on higher-dimensional (  > 1) irreps directly using ( 5) is slower, less reliable, and more liable to instabilities than one-dimensional irreps, as shown in Fig. 2.This may be caused either by the training "wandering" between different wave functions in symmetry-protected multiplets, or by the zero characters   = 0 typical in these irreps reducing the number of wave-function terms in the sum (5), which is known to limit the expressivity of NQS ansätze [43].We remedied this problem by imposing further symmetry constraints that select a unique representative of each symmetry multiplet.Effectively, we project our wave functions first on the trivial irrep of a subgroup of , followed by projecting on the desired irrep of  itself.As explained in Appendix A, the combined effect of these projections can still be written in the form ( 5) with an effective character χ, which is no longer an irrep character of , but has overlap with precisely one of them.The benefits of this approach are illustrated for the five-dimensional H g ( = +1) irrep of C 60 in Fig. 2, which shows that subspace projection allows the variational optimiser to reach lower energies in fewer iterations.

C. Details of the numerical experiments
To obtain the results reported below, we used the same GCNN architecture as Ref. [29], illustrated in Fig. 3.We use real-valued kernels , , interspersed with selu activation functions, which allow us to reliably train deep GCNNs [29,44].In the output layer, we combine pairs of feature maps into complex-valued features, exponentiate them, and project the result on the desired irrep: where  is the number of hidden layers and  is the number of (real-valued) feature maps.Including exponentiation in ( 9) is important to represent the wide dynamical range of wavefunction amplitudes.We used GCNNs with eight hidden layers, each composed of 32 (for the C 32 geometry) or 12 (for the C 60 , C 80 geometries) feature maps, containing 174 336, 243 456, and 243 936 real variational parameters for the C 32 , C 60 , and C 80 geometries, respectively.We also compare the performance of these networks with shallower (three-layer) ones in Fig. 2 for the second-lowest-energy spin-singlet state of the C 60 Heisenberg model, which transforms under the H g irrep of the  ℎ point group (cf.Sec.IV).Without the irrep subspace projection described in Sec.III B, increasing network depth leads to a significant improvement in variational energy; however, after applying the projection, both GCNNs equally outperform the unprojected eight-layer one.
The ansätze were trained on a single A100 GPU using the stochastic reconfiguration algorithm implemented in Best GCNN variational energies for the C 32 geometry (symbols), compared with the low-energy spectrum obtained from exact diagonalisation (red, blue, and green dots for  = 0, 1, 2, respectively).The exact and variational energies are also given in Tables IV and V of Appendix B. Lower panel: difference between the variational and lowest exact energies in each sector.
NetKet [42] with learning rate  = 0.02.The variational energies and spin correlation functions reported below were obtained from averaging VMC local estimators of the Hamiltonian and the operators ì   • ì   for every pair of sites , , obtained for the same set of 2 18 = 262 144 samples for all operators.For wave functions projected on one-dimensional irreps, we expect that spin correlators across symmetry-related pairs of spins are equal: therefore, we explicitly averaged these correlators for the plots below.In addition, we summed the local estimators of ì   • ì   for all pairs of sites to obtain an estimate of the total-spin expectation value ⟨ 2 ⟩.

D. C 32 : comparison to exact diagonalisation
We first benchmark our method on the highest-symmetry ( 3ℎ ) isomer of C 32 , labelled II in Ref. [20], where exact diagonalisation is still possible.We extended the exact spectrum in Ref. [20] to all eigenstates below energy −15: The energies and point-group and spin quantum numbers of these states are listed in Table IV in Appendix B. The best variational energies achieved using the GCNN ansatz are listed in Table V and plotted against the exact spectrum in Fig. 4. In every point-group and parity symmetry sector, we achieve excellent Ground-state spin-spin correlation functions ⟨ ì   • ì   ⟩ as a function of graph distance on the C 32 geometry from the GCNN simulation (symbols) and ED (dots).The size of the symbols is proportional to the number of symmetry-related paths with equal correlators.
agreement with the exact results (see bottom panel with difference between exact and variational energies), with variational energies approaching the exact ground states much closer than the first excited state in the given symmetry sector.
Estimates of the total spin ⟨ 2 ⟩ for our optimised wave functions are listed in Table I.In every symmetry sector, we obtain a value extremely close to ( + 1) for an integer spin quantum number , as expected for a fully spin-rotation-symmetric state.Every odd-parity and   = 2 simulation returned states consistent with  = 1 and  = 2, respectively; in the evenparity sector, we find  = 2 ground states in two point-group symmetry sectors.In these sectors, we also find that the optimal variational energies in the  = +1 and   = 2 sectors coincide to very good accuracy, indicating that the two wave functions are the   = 0 and   = 2 states of the same quintet.
Finally, spin correlation functions ⟨ ì   • ì   ⟩ for the GCNN estimate of ground state are shown in Fig. 5 compared to the correlators of the exact ground state (cf.Ref. [20]); the two again match excellently.Due to the relatively low symmetry of the C 32 geometry, there are many inequivalent lattice sites and paths at all graph distances, so we limit ourselves to Best GCNN variational energies for the C 60 geometry (symbols have the same meaning as in Fig. 4) compared to optimal energies in SU(2)-symmetric DMRG [27] (horizontal lines).The error bars on the  = 2 DMRG result are represented by the green background shading.plotting the correlators as a function of graph distance only.We find no clear pattern in the correlators: their signs deviate from Néel order already for next-nearest neighbours, and their magnitudes are spread over a wide range of values at any given graph distance.This is consistent with the strong frustration expected in this molecule with 12 pentagonal and only 6 hexagonal faces.

IV. INCIPIENT NONCOPLANAR ORDER IN C 60
Next, we consider the highest ( ℎ ) symmetry isomer of C 60 , the famous buckminsterfullerene geometry.The optimised VMC energies for both parities, as well as in the   = 2 sector, Spin correlator are shown in Fig. 6 (see also Table VI in the appendix), while the expectation values of  2 are listed in Table II.Similar to the C 32 geometry, they are consistent with fully spin-rotationsymmetric states.In four symmetry sectors, the lowest-lying even-parity state is not a singlet but a quintet: This is also evidenced by the near coincidence of the optimised energies for  = +1 and   = 2.
Our variational energies match those of the lowest-energy  = 0, 1, 2 states, as well as the first excited  = 0 state, found in a recent SU(2)-symmetric DMRG study [27].The ground state is found to transform under the trivial irrep of the  ℎ point group.In agreement with DMRG, the first excited state is a spin-triplet; the T 2g irrep we identify is also qualitatively consistent with the cubic-symmetry-breaking pattern seen in the DMRG wave function.The lowest-energy  = 2 state is in the H g irrep, with A g and H u states at only slightly slightly higher energies [45]. in Fig. 7. Unlike C 32 , every site of the buckminsterfullerene geometry is equivalent, allowing us to also display the spatial distribution of the correlators.Our results are very close to the correlation functions measured in DMRG [27].For short graph distances, they follow an alternating sign pattern consistent with Néel ordering, and their amplitudes are close to those on the unfrustrated honeycomb lattice at the same graph distances, which we computed using stochastic series expansion [46] with the ALPS library [47,48].Further away, frustration reduces correlators and introduces a nontrivial sign structure, which, somewhat surprisingly, matches that of the classical ground state discussed in Ref. [18] (cf.their Fig. 1).
[For eight sites, marked with symbols in Fig. 7(a), the classical correlator is zero, while the  = 1/2 correlator is anomalously low.]Remarkably, at the largest graph distances, we again recover a Néel-like pattern with amplitudes around ±0.02; however, their signs are inverted compared to the honeycomb lattice.
To extract signatures of ordering from our data, we computed the eigenvalues of the correlation matrix    = ⟨ ì   • ì   ⟩, as well as the the adjacency matrix of the C 60 graph, which plays the role of the Hamiltonian matrix in the Luttinger-Tisza method.These spectra are plotted in Fig. 8.The ground state of the Luttinger-Tisza Hamiltonian is threefold degenerate, forming a T 2g irrep of  ℎ .This is compatible with the unit-length constraint of classical spins: the   ,   ,   components of the ground state form an orthonormal basis of the irrep, leading to a noncoplanar ground state, which can indeed be found by numerically minimising the classical Hamiltonian [18].Due to quantum fluctuations, the spectrum of the spin-1/2 correlation matrix    does not only contain this irrep; however, the weight (that is, the eigenvalue times the irrep dimension) of other irreps is suppressed roughly exponentially in the classical energy cost (note the logarithmic scale in the bottom panel of Fig. 8).In particular,    is dominated by the two lowest-energy irreps of the classical Hamiltonian, T 2g and G u .The two lowest-lying spin-triplet states also belong to these irreps, as expected from the tower-of-states construction of Sec.II.In particular, we find that the overlap between the state Ŝ |GS⟩, generated by the Goldstone-mode operator (3) applied to the ground-state |GS⟩, and the lowest-energy T 2g triplet (  = 0, odd parity) state is ≈ 0.917, very high for two 60-spin many-body states.
Applying the tower-of-states analysis introduced in Refs.[5,12] to the noncoplanar classical ground state correctly predicts that the lowest-lying  = 2 state transforms under the H g irrep.We can also reach this conclusion by applying the bosonic [34] Goldstone-mode operators (3) twice to the ground state.There is a total of nine such operators (threefold spatial and spin degeneracy), so the two-Goldstone Hilbert space consists of 9 × 10/2 = 45 states, transforming under the symmetric square of the T 2g ⊗ ( = 1) irrep of  ℎ × SU(2): The A g singlet and the T 2g triplet cannot be distinguished from the ground state and the one-Goldstone state based on symmetry quantum numbers; in fact, we expect them to have a high overlap.The H g singlet and the A g quintet, however, appear in the spectrum nearly degenerate with the H g quintet, even though they are not predicted by the tower-of-states analysis.
The similarly low-energy H u quintet cannot be explained based on T 2g operators alone, but it may come from a combination of low-lying T 2g and G u triplet excitations.
A representative (maximally symmetric around the centre of the Schlegel plot) eigenvector of the dominant T 2g irrep of    is plotted in Fig. 9(a).This eigenvector follows a perfect Néel pattern on ten hexagons around the "equator" of the C 60 structure, which is in fact its largest unfrustrated portion: this indicates a clear tendency towards the Néel ordering expected in the limit of large molecules.Away from these hexagons, frustration causes the eigenvector components to vanish.The concentration of the Goldstone wave function around the equator of the molecule matches the distribution of "local spin" in the lowest-energy triplet state obtained in DMRG [27].Eigenvectors in the next-highest-weight G u irrep [Fig.9(b)] live mostly on the same set of hexagons and display a Néel pattern modulated with a standing wave, analogous to a long-wave magnon excitation on a lattice model.The structure of the leading correlation eigenvector explains some surprising features of the spin correlators in Fig. 7, in particular the sign inversion of the "Néel order" seen at the largest graph distances.These correspond to pairs of sites belonging to opposite hexagons: Correlations between them are mostly mediated through the "equatorial belt" where the leading correlation eigenvectors are located.Along this belt, antipodal sites are 10 steps apart [see, e.g., the green path in Fig. 9(a)], therefore, they have the same sign in the Néel pattern of the eigenvector, consistent with their positive correlation shown in Fig. 7.However, the shortest path [e.g., the yellow path in Fig. 9(a)] between the same points only takes nine steps, so we observe positive correlations at odd graph distance, in an apparent inversion of the Néel pattern.However, since all of these odd-length paths pass through the fully-frustrated central, low-weight, region of the Schlegel plot, they do not contribute to the correlation function.Likewise, neighbours of the antipodal point are nine (eight) steps apart along the equatorial belt (frustrated region), which extends the inverted Néel pattern to this graph distance too.

V. CHIRAL INVERSION-SYMMETRY BREAKING IN C 80
After buckminsterfullerene, the smallest fullerene structure with full icosahedral symmetry is the 80-site molecule shown  in Fig. 1(c).The converged variational energies in all symmetry sectors are shown in Fig. 10 and Table VII; interestingly, the ground-state is found in a nontrivial point-group irrep (namely A u ), similarly to smaller fullerenes [20] and other frustrated magnetic molecules [12,49].The expectation values of  2 for the GCNN wave functions are listed in Table III, which deviate from the expected values ( + 1) substantially more than in the C 60 case.This is to be expected, given the much smaller spacing between energy levels (e.g., we find a quintet state in every space-group irrep within an energy window of 0.2).Nevertheless, all odd-parity states can be identified as predominantly triplet, and most even-parity states as  = 0 or  = 2.An exception is T 2g ,  = +1, whose ⟨ 2 ⟩ is consistent with a 2 : 1 mixture of a singlet and a quintet.Likewise, the T 1g ,   = 2 calculation reproducibly converges to ⟨ 2 ⟩ ≈ 12, consistent with  = 3 rather than  = 2, even though the  = +1 calculation in the same sector converges to a quintet at a lower energy.In both cases, however, the energy difference between the states in question is extremely small and may not be enough to guide the optimisation algorithm to a perfect  2 eigenstate, allowing expressivity limitations of the GCNN ansatz to dominate the optimisation trajectory.To the best of our knowledge, there are no variational-energy benchmarks for C 80 to which our results could be compared.
The most striking features of the spectrum in Fig. 10 are the near-degenerate singlet "ground states" in the two onedimensional irreps A u , A g , as well as the presence of triplet and quintet excitations in almost all symmetry sectors within a very narrow energy range.We will not attempt to account for every state in this dense spectrum, but only highlight the lowest-energy pair of states in each spin sector (circled in Fig. 10), all of which follow the pattern of incipient inversionsymmetry breaking seen for the ground state: A u , A g for  = 0; T 2g , T 2u for  = 1; and H u , H g for  = 2.
Unlike the C 60 geometry, there are two symmetryinequivalent kinds of lattice site: 60 sites (labelled P) belong to one of the 12 pentagonal faces; the remaining 20 (labelled H) are surrounded by three hexagonal faces each.Fig. 11 shows the spatial structure of spin correlation functions around both kinds of site in the (A u ) ground state.(The correlation structure of the A g singlet is visually indistinguishable.)The sign of correlations again alternates with graph distance, reminiscent of Néel ordering; unlike C 60 , however, this alternating pattern persists without any anomalies all the way to the antipodal points.Similar to C 60 , the magnitude of correlators dips at intermediate graph distances before increasing and levelling off for the largest distances at a typical value of about ±0.05, well above the equivalent figure for C 60 .This is consistent with the diminishing effect of frustration expected for large fullerenes: In the limit of an infinitely large fullerene "molecule" (which, however, still has only 12 pentagonal faces), we expect to recover the ground-state behaviour of the honeycomb Heisenberg antiferromagnet, which forms a Néel order with spin correlator [32] in the long-distance limit [cf.solid dots in Fig. 11(c)].
These features of the ground state and the low-lying spectrum can again be accounted for in terms of incipient Néel ordering.The ground state of the large- Heisenberg Hamiltonian, obtained either by direct simulation or the Luttinger-Tisza method, is again noncoplanar, transforming under the T 2u irrep of  ℎ .The same irrep is also dominant in the spectrum of the correlator matrix (Fig. 12); the gap to subleading eigenvalues is increased compared to C 60 , as expected for an incipient Bragg peak.
A representative eigenvector from this leading T 2u irrep is plotted in Fig. 13.Similar to C 60 , it forms a Néel pattern on the largest unfrustrated subgraph of the fullerene structure: The 10 hexagons (green) on the "equator" of the molecule are, however, laid out differently in the two molecules, thus antipodal points in C 80 acquire opposite signs in the Néel pattern, which explains the inversion-odd irrep.Furthermore, the unfrustrated support of the Néel pattern includes 10 further hexagons (yellow); however, these are separated from one another by pentagons, which frustrate and reduce the ordering amplitude.
12. Top: eigenvalues of the classical Hamiltonian matrix (i.e., the adjacency matrix) in the C 80 geometry.Note that the axis is reversed, so the lowest eigenvalue (the Luttinger-Tisza ground state) is to the right.Bottom: log-scale plot of the weight (eigenvalue times degeneracy) of the eigenspaces of the spin correlator matrix The eigenvalue of one A g eigenvector, corresponding to the net magnetisation, is zero within Monte Carlo error and is not plotted.
Similar to C 60 , the low-lying spectrum can be predicted either by applying the tower-of-states formalism to the classical ground state, or by constructing T 2u Goldstone operators (3) from the leading eigenvectors of the correlator matrix.The first generates the inversion-broken low-energy spectrum seen numerically: A g , A u for  = 0; T 2g , T 2u for  = 1; and H g , H u for  = 2. On the other hand, acting with the Goldstone operators on the A u , A g ground states once yields T 2g , T 2u triplets, while the two-Goldstone manifold includes H u , H g quintets [50].The latter construction also accounts for the energy ordering of the nearly degenerate pairs: those derived from the A u ground state (T 2g triplet, H u quintet) are all lower in energy than the tower of the A g singlet (T 2u triplet, H g quintet).
The Goldstone-mode operators also give microscopic ac- count of the apparent inversion-symmetry breaking in the spectrum.Since the T 2u irrep is threefold degenerate, we can construct three independent triplet operators using (3), which we label as ì S 1 , ì S 2 , ì S 3 .Using all three of these, we can uniquely construct the singlet operator X = ì S 1 • ( ì S 2 × ì S 3 ): one can verify that this operator transforms under the A u irrep [51], thus it might map the A g and A u singlet "ground states" on one another.Indeed, the overlap of X|A u(g) ⟩ and |A g(u) ⟩ is 0.743 (0.693), high values for 80-site many-body states.
The operator X is odd under both inversion and time-reversal symmetry, which strongly suggests that the incipient breaking of inversion symmetry is chiral in nature.The situation is somewhat similar to that of the tetrahedral order found in the triangular-lattice  1 −  2 −   model [5,6]: its tower of states can be captured in terms of three gapless triplet operators, the triple product of which breaks both mirror and time-reversal symmetry.This order, however, is only stabilised on the triangular lattice by a substantial   coupling, which breaks these symmetries explicitly, while in C 80 , the incipient chiral order emerges spontaneously.

VI. LARGER FULLERENES
It is very natural to ask what determines whether a given fullerene geometry supports such incipient chiral-symmetry breaking.The mechanism proposed above severely restricts the number of suitable symmetry groups and allotropes, as requires at least a noncoplanar classical ground state (usually associated with a threefold degenerate ground state in the Luttinger-Tisza spectrum, or a threefold degenerate leading "Bragg-peak" eigenvalue in the quantum correlation matrix), as well as inversion symmetry (so it can be broken).In the following, we explore families of fullerene structures that satisfy these requirements; for a numerical exploration of classical ground-state degeneracy in other geometries, see Ref. [52].
Let us first focus on fullerenes with full icosahedral ( ℎ ) symmetry.These can all be constructed by covering the 20 faces of an icosahedron with patches of the honeycomb lattice [53].For  ℎ symmetry, this covering has to be symmetric under every symmetry of a single triangle, which leads to two distinct series of solutions (Fig. 14): (i) there are  2 lattice points on every face, none of which is shared between more than one face, for a total of 20 2 sites; (ii) some honeycomb edges are aligned with the edges of the icosahedron, yielding a total of 60 2 sites.We conjecture that the incipient Bragg peaks of both kinds of structure are threefold degenerate: more specifically, we expect that the leading correlation eigenvectors show a pronounced Néel pattern on 10 of 20 faces of the icosahedron, which fades away once frustration occurs, as sketched in Fig. 14.We can associate a Néel order parameter with each face, which remains well-defined (i.e., does not change sign) under proper rotations; therefore, the pattern transform under the same irrep of the chiral icosahedral group  as the classical ground state of the C 20 Heisenberg model, namely T 2 .Under inversion, however, the two series behave differently: in the C 20 2 case, inversion-related sites have opposite spins (in fact, each triangle has a net magnetisation, which is odd under inversion), while the C 60 2 ansatz is inversion-even.That is, the lowest triplet excitation of all C 20 2 fullerenes is expected to transform as T 2u (leading to incipient chiral order), while C 60 2 molecules have low-lying T 2g triplet excitations.We verified this prediction by numerically obtaining the classical ground states of the two sequences up to  = 6 (C 720 ) and  = 4 (C 960 ), respectively.
Beyond, icosahedral symmetry, the ingredients proposed above are also available in the largest achiral subgroup of  ℎ , the pyritohedral group  ℎ , which too has three-dimensional irreps T g and T u (corresponding to both T 1 and T 2 of  ℎ ).This is particularly important for crystalline compounds of fullerenes, which cannot maintain full icosahedral symmetry.However, in systems such as the superconductor K 3 C 60 [26], the fullerenes retain full pyritohedral symmetry [54], thus our results remain relevant for their behaviour.More generally, noncoplanar ground states may arise even if they are composed of more than one irrep of the point group, especially if the classical ground state predicted by the Luttinger-Tisza method has support only on parts of the molecule [55].Such classical ground states would still support towers of states with three low-lying spin-triplet excitations, which may combine to generate an operator that breaks inversion and time-reversal symmetry.

VII. CONCLUSION
In summary, we demonstrated that competition between the Néel ordering tendency of the honeycomb lattice and frustrated pentagonal faces leads to a feature-rich incipient noncoplanar order in the quantum Heisenberg model on large fullerene lattices.We generalised a number of techniques commonly used to diagnose ordering tendencies on finite lattices to the molecular case, which lacks translation symmetry: 1. Bragg peaks can be defined in real space as diverging eigenvalues of the spin correlator matrix; in a molecule, they are labelled by point-group irreps rather than wave vectors.
2. Even in the deep quantum limit  = 1/2, these irreps can be predicted from the ground state of the large- version of the Hamiltonian, which can be constructed analogously to the Luttinger-Tisza method.
3. The low-energy spectrum is characterised by a tower of states, the quantum numbers of which can be predicted using ansatz Goldstone-mode operators constructed from the incipient Bragg peaks.
We used these approaches to analyse the low-energy spectrum and ground-state wave function of the nearest-neighbour spin-1/2 Heisenberg model on the icosahedrally-symmetric C 60 and C 80 fullerene geometries.Our numerical results were obtained from variational Monte Carlo using groupconvolutional neural-network (GCNN) wave-function ansätze, which allowed us to construct the symmetry-resolved lowenergy spectrum in detail.We benchmarked the method against ED on a C 32 allotrope and DMRG on C 60 , achieving excellent variational energies in both cases.
For buckminsterfullerene, we found an incipient Bragg peak transforming under the T 2g irrep of the icosahedral point group, which allows for the formation of a Néel pattern on the largest unfrustrated subset of the fullerene graph.This matches the noncoplanar ground state of the classical Heisenberg model [18]; furthermore, the low-energy spectrum of  = 0, 1, 2 excitations follows a tower-of-states structure derived from a triplet Goldstone-mode excitation transforming under the same irrep, pointing towards an incipient noncoplanar order with pronounced Néel-like features on this relatively small and highly frustrated system.
We find similar ordering tendencies on the C 80 geometry.The structure of the incipient Bragg peak is again determined by maximally covering the graph with a Néel pattern, which leads to a triplet of leading correlation eigenvectors transforming under the T 2u irrep.The tower of states corresponding to this inversion-odd incipient order consists of pairs of nearly degenerate (multiplets of) states, distinguished by their inversion eigenvalue.We can relate the wave functions of these states using an explicit operator constructed from Goldstone modes, which is odd under both inversion and time-reversal, indicating that the incipient ordering in C 80 is chiral.Such inversion-symmetry breaking may have interesting ramifications for optical probes: for instance, by breaking down the rule of mutual exclusion [56], it may make Raman-active modes of C 80 visible in infrared spectroscopy or vice versa.Furthermore, the noncoplanar magnetic textures of both molecules may induce anomalous Hall response [57] that may in turn affect the superconductivity of such materials as K 3 C 60 .Here too, the chiral magnetic ground state of C 80 may open the door to more exotic superconducting behaviour, which will be interesting to explore in future theoretical and experimental work.
Finally, we argue that our findings are not limited to the C 60 and C 80 geometries, but are relevant for much wider class of fullerene geometries.A case in point are the two sequences of fullerenes with full icosahedral symmetry, shown in Fig. 14, whose differing geometries lead to chiral incipient ordering in one sequence but not in the other, a surprisingly persistent frustration effect even in the limit of very large molecules.Numerically exploring these large molecules, fullerenes with lower symmetry, other magnetic molecules (e.g. with icosidodecahedral symmetry and larger spins [58,59]), as well as the intermediate- Hubbard model more relevant for real fullerenes, are all exciting directions for future work, calling for improvements to our current neural-quantum-state techniques, which will in turn also benefit studies of lattice models.From a technical point of view, tensor-network studies of C 80 would also be interesting, as they may show time-reversal symmetry breaking explicitly [60]. of this work.NQS simulations were performed using the NetKet [42] library.Reference stochastic series expansion data for the honeycomb lattice was obtained using the ALPS [48] library.All heat maps use perceptionally uniform colour maps developed in Ref. [61].Computing resources were provided by STFC Scientific Computing Department's STFC Cloud service.This work was granted access to the HPC resources of CALMIP center under the allocation 2022-P0677 as well as GENCI (grant A0130500225).This study has been (partially) supported through the EUR grant NanoX No. ANR-17-EURE-0009 in the framework of the "Programme des Investissements d'Avenir."This work benefited from the support of the project QMAHT ANR-22-CE30-0032-03 of the French National Research Agency (ANR).A. Sz.gratefully acknowledges the ISIS Neutron and Muon Source and the Oxford-ShanghaiTech collaboration for support of the Keeley-Rutherford fellowship at Wadham College, Oxford.For the purpose of open access, the authors have applied a Creative Commons Attribution (CC-BY) licence to any author accepted manuscript version arising.

Appendix A: Projecting on subspaces of higher-dimensional irreps
As outlined in Sec.III B, we restricted ansätze transforming under higher-dimensional irreps of the space group  onto a one-dimensional subspace of the symmetry-protected multiplets by imposing additional symmetry constraints, as follows.Consider an Abelian subgroup  of the space group .Restricting any irrep  of  onto  gives a valid representation thereof; for multi-dimensional irreps, however, this is no longer an irrep of , but can be decomposed into them.If an irrep   of  appears with multiplicity 1 in this decomposition, we can select a unique representative of  by first projecting onto   using (5), then onto  itself: This has the same form as the original projection (5), with the "effective character" For simplicity, we choose  for each irrep such that their trivial irrep appears with multiplicity 1, so Eq.(A2) reduces to averaging characters in (right) cosets of .
The  3ℎ point group of C 32 has two two-dimensional irreps, E ′ and E ′′ .A suitable choice of  for both is the  2 subgroup generated by one of the three 180 • rotations.
The  ℎ point group of C 60 and C 80 has three-, four-, and five-dimensional irreps.We decomposed each of these using the following : • T irreps: any  3 rotation subgroup • G irreps: the  2 subgroup that fixes an edge of a dodecahedron with the same symmetry group.In C 60 , the same group fixes an edge lying between two hexagons.
In C 80 , it fixes the hexagon diagonal lying between two nearest-neighbour pentagons.
Each of these groups decompose the given irrep into every one of their own irreps, each with multiplicity 1. Figures 9 and 13 are exceptions: there, we decomposed irreps of both the correlation eigenvectors and the low-lying wave functions into irreps of the  5 rotation group around the centre of the Schlegel plot.Such a decomposition of T irreps contains the trivial irrep, yielding rotationally symmetric plots.The G u irrep in Fig. 9(b), by contrast, decomposes into the four nontrivial irreps of  5 : for the plot, we used the real part of the  2 /5 rotation-eigenvalue component (that is, a linear combination of the  ±2 /5 components).

FIG. 7 .
FIG. 7. (a) Ground-state spin-spin correlation functions ⟨ ì  • ì   ⟩ in the C 60 geometry.The reference point  is marked with a black triangle.Two values below 0.005 in magnitude (highlighted with coloured symbols) were truncated for visibility.(b) Spin-spin correlators as a function of graph distance.Red and blue symbols stand for positive and negative correlators, respectively.Coloured dots show the spin correlation functions of a 512-site honeycomb lattice, measured using QMC; the dashed line is a spline connecting these dots and is included as a guide to the eye.

FIG. 9 .
FIG. 9. Eigenvectors of the ground-state correlator matrix ⟨ ì   • ì   ⟩ in the C 60 geometry corresponding to its largest [(a); T 2g irrep; 20.7% of all correlations] and second largest [(b); G u irrep; 20.4%] eigenvalues.The ten hexagons on the "equator" of the C 60 structure are highlighted in green.The yellow line in panel (a) indicates the shortest path (9 steps) connecting two antipodal points; the green line is the shortest path (10 steps) passing through sites with nonzero amplitude in the eigenvector.

2 A 1 𝑆 𝑧 = 2 FIG. 10 .
FIG.10.Best GCNN variational energies for the C 80 geometry.Symbols have the same meaning as in Figs.4 and 6.The two lowestenergy states in each spin-quantum-number sector are circled with dashed lines.

FIG. 11 .
FIG. 11. (a,b) Ground-state spin-spin correlation functions ⟨ ì   • ì   ⟩ in the C 80 geometry for P-type (a) and H-type (b) reference points  (marked with a black triangle).(c) Spin-spin correlations as a function of graph distance.Red and blue symbols stand for positive and negative correlators, respectively.Coloured dots show the spin correlation functions of a 512-sites honeycomb lattice, measured using QMC; the dashed line is a spline connecting these dots and is included as a guide to the eye.One value below 0.01 in magnitude [1.90(8) × 10 −3 , at graph distance 6, highlighted with coloured symbols in (a)] was truncated for visibility.

FIG. 13 .
FIG. 13.Eigenvector of the ground-state correlator matrix ⟨ ì   • ì   ⟩ in the C 80 geometry corresponding to its largest eigenvalue.Green and yellow hexagons indicate a maximal unfrustrated portion of the C 80 structure.

FIG. 14 .
FIG.14.Structure of the ansatz for the dominant eigenvector of    on fullerenes of the  ℎ -symmetric (i) C 20 2 and (ii) C 60 2 series, shown on the flattened net of an icosahedron.They form a Néel pattern of positive (red) and negative (blue) amplitudes on an unfrustrated belt covering 10 of the 20 icosahedron faces, and have only small support further away (gray).Under proper rotations, both transform as the T 2 irrep of .Under inversion, the two green triangles map onto each other; the arrangement of spins transforms as if the triangle was horizontally flipped on the diagram.This flips the sign of the C 20 2 structure and preserves that of C 60 2 .
+1) irrep of the C 60 fullerene structure.The energies are compared to the second-lowest spin-singlet energy found in Ref.
[29]cture of a five-layer GCNN of the type used in this work.Red and blue boxes stand for the embedding layer (6a) and the group convolutions (6b), respectively.Green boxes indicate selu activation functions.Yellow boxes represent the output layer(9).Figure based on Ref.[29].

TABLE I .
Total ⟨ 2 ⟩ for the optimised GCNN wave functions on the C 32 geometry.All are close to ( + 1) for an integer spin quantum number , indicating an accurately spin-rotation-symmetric wave function.

TABLE II
. Total ⟨ 2 ⟩ for the optimised GCNN wave functions on the C 60 geometry.All but one (in italics) are very close to ( + 1) for an integer spin quantum number , indicating an accurate spin-rotationsymmetric wave function.
Top: eigenvalues of the classical Hamiltonian matrix (i.e., the adjacency matrix) in the C 60 geometry.Note that the axis is reversed, so the lowest eigenvalue (the Luttinger-Tisza ground state) is to the right.Bottom: log-scale plot of the weight (eigenvalue times degeneracy) of the eigenspaces of the spin correlator matrix⟨ ì   • ì   ⟩.The eigenvalue of one A g eigenvector, corresponding to the net magnetisation, is zero within Monte Carlo error and is not plotted.

TABLE III .
Total ⟨ 2 ⟩ for the optimised GCNN wave functions on the C 80 geometry.

TABLE IV .
Energies, spin quantum numbers, and point-group irreps of every eigenstate of the C 32 Heisenberg model below energy −15 from exact diagonalisation.

TABLE V .
Best GCNN variational energies for the C 32 geometry.Bold numbers correspond to the lowest energy in each spin sector.We note that the lowest  = 2 state found by exact diagonalisation belongs to the  ′ 1 irrep; however, the variational energy difference between the  ′ 1 and  ′′ 1 sectors is smaller than the difference of either to the exact value.

TABLE VI .
[27] GCNN variational energies for the C 60 geometry, compared to the DMRG variational energies of Ref.[27].Bold numbers correspond to the lowest energy in each spin sector.+1 simulation that returned an  = 2 state.b  = 2 simulation that returned an  = 3 state.c State not clearly dominated by one -sector.
a  =

TABLE VII .
Best GCNN variational energies for the C 80 geometry.Bold numbers correspond to the lowest energy in each spin sector.