<figure><img src="image_1.png" /></figure>

HAL
open science

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

Attila Szabo, Sylvain Capponi, Fabien Alet

► To cite this version:

Attila Szabo, Sylvain Capponi, Fabien Alet. Noncoplanar and chiral spin states on the way towards Néel ordering in fullerene Heisenberg models. Physical Review B, 2024, 109 (5), pp.054410. 10.1103/PhysRevB.109.054410. hal-04566181

HAL Id: hal-04566181

https://hal.science/hal-04566181

Submitted on 3 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Noncoplanar and chiral spin states on the way towards Néel ordering in fullerene Heisenberg models

Attila Szabó ⚪,¹,²,³ Sylvain Capponi ⚪,⁴ and Fabien Alet ⚪,⁴

¹Rudolf Peierls Centre for Theoretical Physics, University of Oxford, Oxford OX1 3PU, United Kingdom

²ISIS Facility, Rutherford Appleton Laboratory, Harwell Campus, Didcot OX11 0QX, United Kingdom

³Max Planck Institute for the Physics of Complex Systems, Nöthnitzer Straße 38, 01187 Dresden, Germany

⁴Laboratoire de Physique Théorique, Université de Toulouse, CNRS, UPS, 31062 Toulouse, France

<figure><img src="image_3.png" /></figure>

(Received 21 November 2023; accepted 3 January 2024; published 8 February 2024)

Using a high-accuracy variational Monte Carlo approach based on group-convolutional neural networks, we obtain the symmetry-resolved low-energy spectrum of the spin-1/2 Heisenberg model on several highly symmetric fullerene geometries, including the famous C₆₀ buckminsterfullerene. We argue that as the degree of frustration is lowered in large fullerenes, they display characteristic features of incipient magnetic ordering: Correlation functions show high-intensity Bragg peaks consistent with Néel-like ordering, while the low-energy spectrum is organized into a tower of states. Competition with frustration, however, turns the simple Néel order into a noncoplanar one. Remarkably, we find and predict chiral incipient ordering in a large number of fullerene structures.

DOI: 10.1103/PhysRevB.109.054410

## I. INTRODUCTION

Antiferromagnets on infinite bipartite lattices generally show Néel ordering in dimensions greater than one. This can be detected through the staggered magnetization 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 $S \ge 1/2$ [2,3] as well as on the two-dimensional square lattice for $S \ge 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–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–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₂₀, 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 effects in carbon fullerenes, exact-diagonalization (ED) and quantum Monte Carlo (QMC) studies were performed on the C₂₀ 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–19] intermediate-U Hubbard model, which would pose a considerably greater computational challenge. The lowest-lying excitations of the C₂₀ 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₃₆ [20–22], their high degree of frustration and varying degrees of symmetry obstruct the emergence of any systematic trend.

Much attention has also been devoted to the famous C₆₀ buckminsterfullerene geometry [Fig. 1(b)], motivated by interest in chemical, nanotechnology, and quantum-computing [25] applications, as well as superconductivity observed in alkali-metal-doped fullerene lattices [26]. Early

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Open access publication funded by the Max Planck Society.</footnote>
<figure><img src="image_4.png" /><figcaption>FIG. 1. Connectivity of the three fullerene graphs considered in this paper. (a) C$_{32}$ ($D_{3h}$ isomer), (b) C$_{60}$ ($I_h$ isomer), and (c) C$_{80}$ ($I_h$ isomer). These illustrations were made using Visualization for Electronic Structural Analysis (VESTA) [23]; coordinates of C$_{32}$ were taken from Ref. [24].</figcaption></figure>

studies of the Hubbard model [18,19] argued in favor of a nondegenerate singlet ground state at strong and intermediate interaction on this geometry as well. More recently, this has been corroborated by density matrix renormalization group (DMRG) studies of the Heisenberg model by Rausch et al. [27], which found the lowest few eigenstates of the $S = 0, 1, 2$ spin sectors. In particular, these authors uncovered that the lowest-lying 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 (VMC) 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 modest computational resources and thus reconstruct much of the low-energy spectrum. In particular, we study the spin-1/2 nearest-neighbor Heisenberg model

$$H = J \sum_{\langle ij \rangle} \vec{S}_i \cdot \vec{S}_j \quad (1)$$

on the highest-symmetry allotropes of C$_{32}$, C$_{60}$, and C$_{80}$, shown in Fig. 1. In what follows, we use $J = 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 ($S \to \infty$) Heisenberg models developed for lattice models. In particular, we find that the lowest-lying $S = 0, 1, 2$ states are captured by a triplet of low-energy $S = 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 low-lying spectrum comes in nearly degenerate pairs of states, which only differ in their parity under spatial inversion. We again account for this behavior in terms of an incipient symmetry-breaking order. The large-$S$ 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 Goldstone-mode 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 organized as follows. In Sec. II, we generalize methods to detect incipient ordering in finite systems to molecules without translation symmetry. In Sec. III, we describe our GCNN ansatz and its optimization 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. Appendices 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 $\langle \vec{S}(-\vec{k}_0) \cdot \vec{S}(\vec{k}_0) \rangle$ at some wave vector $\vec{k}_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

$$ \vec{S}(\vec{k}_0) = \sum_r e^{i\vec{k}_0 \cdot \vec{r}} \vec{S}_r \quad (2) $$

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 $G_{\text{spatial}} \times \text{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 finite-size 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 large-molecule limit, almost every face is hexagonal, and 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 $C_{ij} = \langle \vec{S}_i \cdot \vec{S}_j \rangle$, with a diverging eigenvalue corresponding to the order parameter [33]. Likewise, the leading eigenvector of $C_{ij}$ on the fullerene geometry can be thought of as a real-space Bragg-peak “wave function” $\psi_i$. This wave function can be used to construct the Goldstone-mode operator

$$ \hat{S}_{\psi}^{\pm} = \sum_{i} \psi_{i} \hat{S}_{i}^{\pm}, \quad \hat{S}_{\psi}^{z} = \sum_{i} \psi_{i} \hat{S}_{i}^{z}, \quad (3) $$

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 $\hat{S}_{\psi}$ or from decomposing a symmetry-broken classical ground state into irreps of $G_{\text{point}} \times \text{SU}(2)$, the latter of which can be made controlled in the large-$S$ 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 $C_{ij}$ of the quantum many-body ground state, similar to the Luttinger-Tisza method for lattice magnets [35–37]. In the large-$S$ limit underlying the above arguments, the Hamiltonian (1) on a

lattice can be Fourier transformed,

$$ H = \sum_k J(k) \vec{S}(-k) \cdot \vec{S}(k), \quad (4) $$

without having to worry about complicated commutation relations between the $\vec{S}(k)$. The minimum of $J(k)$ 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 non-coplanar [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 $\psi_i$.

### 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 the neural-network-based variational Monte Carlo approach, 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 (multidimensional) 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 $\psi_0$ using the projection formula [41]

$$ |\psi\rangle = \frac{d_{\chi}}{|G|} \sum_{g \in G} \chi_g^* \hat{s} |\psi_0\rangle, \quad (5a) $$

$$ \psi(\sigma) \equiv \langle \sigma | \psi \rangle = \frac{d_{\chi}}{|G|} \sum_{g \in G} \chi_g^* \psi_0 [\hat{g}^{-1}(\sigma)], \quad (5b) $$

where the $\hat{g}$ are the elements of the space group G and the $\chi_g$ are their characters in a given $d_{\chi}$-dimensional irrep of G. Here, $\sigma$ stands for a spin configuration in the computational $S^z$ basis: $\sigma = (S_1^z, S_2^z, \dots, S_N^z)$ for N spins-1/2 with $S_i^z \in \{\uparrow, \downarrow\} \equiv \{\pm 1\}$. Applying (5b) directly, however, requires evaluating the ansatz $\psi_0$ many times, which may be prohibitively computationally expensive.

Instead, we use group-convolutional neural networks (GCNNs) [28,29,38], a generalization of the well-known convolutional neural networks (CNNs) to nontranslational symmetries, which are able to efficiently generate all symmetry-related evaluations of a neural-network ansatz as hidden layers indexed by the symmetry elements. Our feed-forward GCNNs start with an embedding layer

$$ h_g^{(1)} = \sum_{\vec{r}} K(\hat{g}^{-1}\vec{r})\sigma(\vec{r}), \quad (6a) $$
which converts the input spin configuration $\sigma$ into such a hidden layer, $h^{(1)}$. Then, further group-valued hidden layers are generated by alternating nonlinearities and equivariant linear layers of the form

$$h_g^{(i+1)} = \sum_{k \in G} W^{(i)}(k^{-1}g)h_k^{(i)}. \quad (6b)$$

The trainable variables of the ansatz are the kernel entries $K(\vec{r})$ and $W(g)$.

One can show [42] that acting with a symmetry element $\hat{k}$ on the input spin configuration $\sigma$ permutes the labels of all subsequent layers as

$$h_g[\hat{k}^{-1}(\sigma)] = h_{kg}(\sigma). \quad (7)$$

Therefore, we can regard the entries of the last layer as the amplitudes of a neural-quantum-state (NQS) ansatz $h_0^{(L)}(\sigma)$ evaluated for all spin configurations related to $\sigma$ by space-group symmetry:

$$h_g^{(L)}(\sigma) = h_0^{(L)}[\hat{g}^{-1}(\sigma)]. \quad (8)$$

Thus, 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 $S^z = \sum_{i=1}^N S_i^z$, the parity symmetry $\hat{P} = \prod_{i=1}^N S_i^x$, an element of the SU(2) spin-rotation group, can be implemented in the $S^z$ computational basis by flipping the sign of all $S_i^z$. We therefore imposed eigenvalues of $P = \pm 1$ on our ansätze in addition to the space-group irreps. Sampling in the $S^z = 0$ magnetization sector, this allows us to distinguish between states with even ($P = +1$) and odd ($P = -1$) total-spin quantum number. We also performed simulations in the $S^z = 2$ sector, which isolate total-spin quantum numbers $S \ge 2$.

## B. Irrep subspace projection

In our numerical experiments, we found that training ansätze projected on higher-dimensional ($d_\chi > 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 $\chi_g = 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 $G$, followed by projecting on the desired irrep of $G$ itself. As explained in Appendix A, the combined effect of these projections can still be written in the form (5) with an effective character $\tilde{\chi}$, which is no longer an irrep character of $G$ but has overlap with precisely one of them. The benefits of this approach are illustrated for the five-dimensional $H_g$ ($P = +1$) irrep of $C_{60}$ in Fig. 2, which shows that subspace projection allows the variational optimizer to reach lower energies in fewer iterations.

<figure><img src="image_5.png" /><figcaption>FIG. 2. Evolution of the variational energy during the training protocol for three- and eight-layer GCNN ansätze, both with and without the irrep subspace projection described in Sec. III B, for the $H_g$ ($P = +1$) irrep of the $C_{60}$ fullerene structure. The energies are compared with the second-lowest spin-singlet energy found in Ref. [27] using DMRG (cf. Fig. 6).</figcaption></figure>

## 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 $K, W$, 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:

$$\tilde{h}_{n,g} = h_{n,g}^{(L)} + ih_{n+F/2,g}^{(L)}, \quad (9a)$$

$$\psi(\sigma) = \sum_{g \in G} \chi_g^* \sum_{n=1}^{F/2} \exp(\tilde{h}_{n,g}), \quad (9b)$$

where $L$ is the number of hidden layers and $F$ is the number of (real-valued) feature maps. Including exponentiation in (9) is important to represent the wide dynamical range of wave-function amplitudes.

We used GCNNs with eight hidden layers, each composed of 32 (for the $C_{32}$ geometry) or 12 (for the $C_{60}$ and $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 $I_h$ 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

<figure><img src="image_6.png" /><figcaption>FIG. 3. Structure of a five-layer GCNN of the type used in this paper. Red and blue boxes stand for the embedding layer (6a) and the group convolutions (6b), respectively. Green boxes indicate scaled exponential linear unit (SELU) activation functions. Yellow boxes represent the output layer (9). This figure is based on Ref. [29].</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIG. 4. Best GCNN variational energies for the C$_{32}$ geometry (symbols), compared with the low-energy spectrum obtained from exact diagonalization (red, blue, and green dots for S = 0, 1, 2, respectively). The exact and variational energies are also given in Tables IV and V in Appendix B. Lower panel: difference between the variational and lowest exact energies in each sector.</figcaption></figure>

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 NETKET [42] with learning rate $\eta = 0.02$. To maximally exploit GPU parallelism, we used 1024 parallel Markov chains to generate 3072 Monte Carlo samples per training step. In most simulations, we performed 2000 training steps, which took between 7 (for C$_{32}$) and 28 (for C$_{80}$) GPU hours.

The variational energies and spin correlation functions reported below were obtained from averaging VMC local estimators of the Hamiltonian and the operators $\vec{S}_i \cdot \vec{S}_j$ for every pair of sites i, j, 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 $\vec{S}_i \cdot \vec{S}_j$ for all pairs of sites to obtain an estimate of the total-spin expectation value $\langle S^2 \rangle$.

### D. C$_{32}$: Comparison to exact diagonalization

We first benchmark our method on the highest-symmetry (D$_{3h}$) isomer of C$_{32}$, labeled isomer II in Ref. [20], where exact diagonalization 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 in Appendix B and plotted against the exact

<table><caption>TABLE I. Total $\langle S^2 \rangle$ for the optimized GCNN wave functions on the C$_{32}$ geometry. All are close to $S(S+1)$ for an integer spin quantum number S, indicating an accurately spin-rotation-symmetric wave function.</caption><thead><tr><th>Irrep</th><th>P = +1</th><th>P = -1</th><th>S<sup>z</sup> = 2</th></tr></thead><tbody><tr><td>A'<sub>1</sub></td><td>0.0075(6)</td><td>2.0054(5)</td><td>6.0065(5)</td></tr><tr><td>A''<sub>2</sub></td><td>5.9338(15)<sup>a</sup></td><td>2.0072(7)</td><td>6.0095(7)</td></tr><tr><td>A''<sub>1</sub></td><td>0.0036(4)</td><td>2.0131(9)</td><td>6.0032(4)</td></tr><tr><td>A''<sub>2</sub></td><td>5.9975(7)<sup>a</sup></td><td>2.0187(10)</td><td>6.0060(5)</td></tr><tr><td>E'</td><td>0.0116(7)</td><td>2.0066(5)</td><td>6.0125(8)</td></tr><tr><td>E''</td><td>0.0085(6)</td><td>2.0054(6)</td><td>6.0184(9)</td></tr></tbody></table>

<sup>a</sup>P = +1 simulation that returned an S = 2 state.

spectrum in Fig. 4. In every point-group and parity symmetry sector, we achieve excellent agreement with the exact results (see bottom panel of Fig. 4 with the difference between exact and variational energies), with variational energies approaching the exact ground states much more closely than the first excited state in the given symmetry sector.

Estimates of the total spin $\langle S^2 \rangle$ for our optimized wave functions are listed in Table I. In every symmetry sector, we obtain a value extremely close to $S(S+1)$ for an integer spin quantum number S, as expected for a fully spin-rotation-symmetric state. Every odd-parity and $S^z = 2$ simulation returned states consistent with S = 1 and S = 2, respectively; in the even-parity sector, we find S = 2 ground states in two point-group symmetry sectors. In these sectors, we also find that the optimal variational energies in the P = +1 and $S^z = 2$ sectors coincide to very good accuracy, indicating that the two wave functions are the $S^z = 0$ and $S^z = 2$ states of the same quintet.

Finally, spin correlation functions $\langle \vec{S}_i \cdot \vec{S}_j \rangle$ for the GCNN estimate of the 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

<figure><img src="image_8.png" /><figcaption>FIG. 5. Ground-state spin-spin correlation functions $\langle \vec{S}_i \cdot \vec{S}_j \rangle$ 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.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIG. 6. 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 $S = 2$ DMRG result are represented by the green background shading.</figcaption></figure>

paths at all graph distances, so we limit ourselves to 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 neighbors, 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 ($I_h$) isomer of C$_{60}$, the famous buckminsterfullerene geometry. The optimized VMC energies for both parities, as well as in the $S^z = 2$ sector, are shown in Fig. 6 (see also Table VI in Appendix B), while the expectation values of $S^2$ are listed in Table II. Similar to the C$_{32}$ geometry, they are consistent with fully spin-rotation-symmetric states. In four symmetry sectors, the lowest-lying even-parity state is not a singlet, but a quintet: This is also

<table><caption>TABLE II. Total $\langle S^2 \rangle$ for the optimized GCNN wave functions on the C$_{60}$ geometry. All but one (in italics) are very close to $S(S+1)$ for an integer spin quantum number $S$, indicating an accurate spin-rotation-symmetric wave function.</caption><thead><tr><th>Irrep</th><th><i>P</i> = +1</th><th><i>P</i> = -1</th><th><i>S</i><sup>z</sup> = 2</th></tr></thead><tbody><tr><td>A<sub>g</sub></td><td>0.0022(3)</td><td>2.0121(8)</td><td>6.0038(4)</td></tr><tr><td>A<sub>u</sub></td><td>5.530(3)<sup>a</sup></td><td>2.0541(15)</td><td>6.0279(10)</td></tr><tr><td>T<sub>1g</sub></td><td>6.0013(9)<sup>a</sup></td><td>2.0078(6)</td><td>6.0542(14)</td></tr><tr><td>T<sub>1u</sub></td><td>0.0119(7)</td><td>2.0067(6)</td><td>6.0283(11)</td></tr><tr><td>T<sub>2g</sub></td><td>5.9898(16)<sup>a</sup></td><td>2.0065(7)</td><td>6.0129(8)</td></tr><tr><td>T<sub>2u</sub></td><td>0.0105(7)</td><td>2.0347(12)</td><td>6.0295(12)</td></tr><tr><td>G<sub>g</sub></td><td>0.0241(10)</td><td>2.0356(14)</td><td>6.0105(7)</td></tr><tr><td>G<sub>u</sub></td><td>0.0230(10)</td><td>2.0104(7)</td><td>6.0088(6)</td></tr><tr><td>H<sub>g</sub></td><td>0.0524(14)</td><td>2.0480(15)</td><td>6.0101(7)</td></tr><tr><td>H<sub>u</sub></td><td>5.9873(13)<sup>a</sup></td><td>2.0095(7)</td><td>6.0166(8)</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>FIG. 7. (a) Ground-state spin-spin correlation functions $\langle \vec{S}_i \cdot \vec{S}_j \rangle$ in the C$_{60}$ geometry. The reference point <i>i</i> is marked with a black triangle. Two values below 0.005 in magnitude (highlighted with colored 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. Colored dots show the spin correlation functions of a 512-site honeycomb lattice, measured using QMC approach; the dashed line is a spline connecting these dots and is included as a guide to the eye.</figcaption></figure>

evidenced by the near coincidence of the optimized energies for $P = +1$ and $S^z = 2$.

Our variational energies match those of the lowest-energy $S = 0, 1, 2$ states, as well as the first excited $S = 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 $I_h$ 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 $S = 2$ state is in the $H_g$ irrep, with $A_g$ and $H_u$ states at only slightly higher energies [45].

The ground-state spin correlation function $\langle \vec{S}_i \cdot \vec{S}_j \rangle$ is shown 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

<footnote><sup>a</sup><i>P</i> = +1 simulation that returned an <i>S</i> = 2 state.</footnote>
<figure><img src="image_11.png" /><figcaption>FIG. 8. Top: eigenvalues of the classical Hamiltonian matrix (i.e., the adjacency matrix) in the C60 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 $\langle \vec{S}_i \cdot \vec{S}_j \rangle$. The eigenvalue of one $A_g$ eigenvector, corresponding to the net magnetization, is zero within Monte Carlo error and is not plotted.</figcaption></figure>

series expansion [46] with the Applications and Libraries for Physics Simulations (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 $S = 1/2$ correlator is anomalously low.] Remarkably, at the largest graph distances, we again recover a Néel-like pattern with amplitudes around $\pm 0.02$; however, their signs are inverted compared with the honeycomb lattice.

To extract signatures of ordering from our data, we computed the eigenvalues of the correlation matrix $C_{ij} = \langle \vec{S}_i \cdot \vec{S}_j \rangle$, as well as the the adjacency matrix of the C60 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 $I_h$. This is compatible with the unit-length constraint of classical spins: the $S^x$, $S^y$, $S^z$ 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 minimizing the classical Hamiltonian [18]. Due to quantum fluctuations, the spectrum of the spin-1/2 correlation matrix $C_{ij}$ 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, $C_{ij}$ 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 $\hat{S}^z|GS\rangle$, generated by the Goldstone-mode operator (3) applied to the ground state $|GS\rangle$ and the lowest-energy $T_{2g}$ triplet ($S^z = 0$, odd parity) state is $\approx 0.917$, very high for two 60-spin many-body states.

<figure><img src="image_12.png" /><figcaption>FIG. 9. Eigenvectors of the ground-state correlator matrix $\langle \vec{S}_i \cdot \vec{S}_j \rangle$ in the C60 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 C60 structure are highlighted in green. The yellow line in (a) indicates the shortest path (nine steps) connecting two antipodal points; the green line is the shortest path (ten steps) passing through sites with nonzero amplitude in the eigenvector.</figcaption></figure>

Applying the tower-of-states analysis introduced in Refs. [5,12] to the noncoplanar classical ground state correctly predicts that the lowest-lying $S = 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 are a total of nine such operators (threefold spatial and spin degeneracy), so the two-Goldstone Hilbert space consists of $9 \times 10/2 = 45$ states, transforming under the symmetric square of the $T_{2g} \otimes (S = 1)$ irrep of $I_h \times \text{SU}(2)$:

$$ \text{Sym}^2[\text{T}_{2g} \otimes (S = 1)] = (A_g \oplus H_g) \otimes [(S = 0) \oplus (S = 2)] \oplus \text{T}_{2g} \otimes (S = 1). \quad (10) $$

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 $G_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 center of the Schlegel plot) eigenvector of the dominant $T_{2g}$ irrep of $C_{ij}$ is plotted in Fig. 9(a). This eigenvector follows a perfect Néel pattern on ten hexagons around the “equator” of the C60 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
<figure><img src="image_13.png" /><figcaption>FIG. 10. Best GCNN variational energies for the C$_{80}$ geometry. Symbols have the same meaning as in Figs. 4 and 6. The two lowest-energy states in each spin-quantum-number sector are circled with dashed lines.</figcaption></figure>

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 ten 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, neighbors 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 in Appendix B, 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 $S^2$ for the GCNN wave functions (listed in Table III) deviate from the expected values $S(S+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

<table><caption>TABLE III. Total $\langle S^2 \rangle$ for the optimized GCNN wave functions on the $C_{80}$ geometry.</caption><thead><tr><th>Irrep</th><th>$P = +1$</th><th>$P = -1$</th><th>$S^z = 2$</th></tr></thead><tbody><tr><td>$A_g$</td><td>0.0185(9)</td><td>2.242(3)</td><td>6.0209(10)</td></tr><tr><td>$A_u$</td><td>0.0119(7)</td><td>2.0248(11)</td><td>6.0264(11)</td></tr><tr><td>$T_{1g}$</td><td>6.019(2)<sup>a</sup></td><td>2.233(3)</td><td>11.766(3)<sup>b</sup></td></tr><tr><td>$T_{1u}$</td><td>5.958(2)<sup>a</sup></td><td>2.153(3)</td><td>6.064(2)</td></tr><tr><td>$T_{2g}$</td><td>1.998(6)<sup>c</sup></td><td>2.0301(13)</td><td>6.093(2)</td></tr><tr><td>$T_{2u}$</td><td>6.010(2)<sup>a</sup></td><td>2.075(2)</td><td>6.164(2)</td></tr><tr><td>$G_g$</td><td>0.466(3)</td><td>2.118(2)</td><td>6.085(2)</td></tr><tr><td>$G_u$</td><td>0.704(4)</td><td>2.179(3)</td><td>6.122(2)</td></tr><tr><td>$H_g$</td><td>0.343(3)</td><td>2.0398(14)</td><td>6.147(2)</td></tr><tr><td>$H_u$</td><td>0.350(3)</td><td>2.093(2)</td><td>6.091(2)</td></tr></tbody></table>

<sup>a</sup>$P = +1$ simulation that returned an $S = 2$ state.

<sup>b</sup>$S^z = 2$ simulation that returned an $S = 3$ state.

<sup>c</sup>State that is not clearly dominated by one spin-quantum-number sector.

window of 0.2J). Nevertheless, all odd-parity states can be identified as predominantly triplet, and most even-parity states can be identified as $S = 0$ or $S = 2$. An exception is $T_{2g}$, $P = +1$, whose $\langle S^2 \rangle$ is consistent with a 2:1 mixture of a singlet and a quintet. Likewise, the $T_{1g}$, $S^z = 2$ calculation reproducibly converges to $\langle S^2 \rangle \approx 12$, consistent with $S = 3$ rather than $S = 2$, even though the $P = +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 optimization algorithm to a perfect $S^2$ eigenstate, allowing expressivity limitations of the GCNN ansatz to dominate the optimization 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 one-dimensional 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 inversion-symmetry breaking seen for the ground state: $A_u$, $A_g$ for $S = 0$; $T_{2g}$, $T_{2u}$ for $S = 1$; and $H_u$, $H_g$ for $S = 2$.

Unlike the $C_{60}$ geometry, there are two symmetry-inequivalent kinds of lattice site: 60 sites (labeled “P”) belong to one of the 12 pentagonal faces; the remaining 20 (labeled “H”) are surrounded by three hexagonal faces each. Figure 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 leveling off for the largest distances at a typical value of about $\pm 0.05$, well above the equivalent figure for $C_{60}$. This is consistent with the diminishing effect of frustration expected
<figure><img src="image_14.png" /><figcaption>FIG. 11. (a) and (b) Ground-state spin-spin correlation functions $\langle \vec{S}_i \cdot \vec{S}_j \rangle$ in the $C_{80}$ geometry for P-type (a) and H-type (b) reference points $i$ (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. Colored dots show the spin correlation functions of a 512-site honeycomb lattice, measured using the QMC approach; 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) \times 10^{-3}$, at graph distance 6, highlighted with colored symbols in (a)] was truncated for visibility.</figcaption></figure>

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 behavior of the honeycomb Heisenberg antiferromagnet, which forms a Néel order with spin correlator $\langle \vec{S}_0 \cdot \vec{S}_r \rangle = \pm 0.0717(3)$ [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-$S$ Heisenberg Hamiltonian, obtained either by direct simulation or by the Luttinger-Tisza method, is again noncoplanar, transforming under the $T_{2u}$ irrep of $I_h$. The same irrep is also dominant in the spectrum of the spin-1/2 correlator matrix (Fig. 12); the gap to subleading eigenvalues is increased compared with $C_{60}$, as expected for an incipient Bragg peak.

<figure><img src="image_15.png" /><figcaption>FIG. 12. Top: eigenvalues of the classical Hamiltonian matrix (i.e., the adjacency matrix) in the $C_{80}$ geometry. Note that the axis is reversed and 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 $\langle \vec{S}_i \cdot \vec{S}_j \rangle$. The eigenvalue of one $A_g$ eigenvector, corresponding to the net magnetization, is zero within Monte Carlo error and is not plotted.</figcaption></figure>

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 ten hexagons (green) on the “equator” of the molecule are, however, laid out differently in the two molecules; therefore, 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 ten further hexagons (yellow); however, these are separated from one another by pentagons, which frustrate and reduce the ordering amplitude.

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

<figure><img src="image_16.png" /><figcaption>FIG. 13. Eigenvector of the ground-state correlator matrix $\langle \vec{S}_i \cdot \vec{S}_j \rangle$ in the $C_{80}$ geometry corresponding to its largest eigenvalue. Green and yellow hexagons indicate a maximal unfrustrated portion of the $C_{80}$ structure.</figcaption></figure>
numerically: $A_g$, $A_u$ for $S = 0$; $T_{2g}$, $T_{2u}$ for $S = 1$; and $H_g$, $H_u$ for $S = 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 a microscopic account 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 $\vec{\bar{S}}_1, \vec{\bar{S}}_2, \vec{\bar{S}}_3$. Using all three of these, we can uniquely construct the singlet operator $\hat{\mathcal{X}} = \vec{\bar{S}}_1 \cdot (\vec{\bar{S}}_2 \times \vec{\bar{S}}_3)$. One can verify that this operator transforms under the $A_u$ irrep [51], so it might map the $A_g$ and $A_u$ singlet “ground states” on one another. Indeed, the overlap of $\mathcal{X}|A_{u(g)}\rangle$ and $|A_{g(u)}\rangle$ is 0.743 (0.693), high values for 80-site many-body states.

The operator $\mathcal{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 $J_1-J_2-J_\chi$ 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 stabilized on the triangular lattice by a substantial $J_\chi$ 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 it 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 ($I_h$) symmetry. These can all be constructed by covering the 20 faces of an icosahedron with patches of the honeycomb lattice [53]. For $I_h$ symmetry, this covering has to be symmetric under every symmetry of a single triangle, which leads to two distinct series of solutions (Fig. 14): (a) There are $n^2$ lattice points on every face, none of which is shared between more than one face, for a total of $20n^2$ sites; (b) some honeycomb edges are aligned with the edges of the icosahedron, yielding a total of $60n^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

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>FIG. 14. Structure of the ansatz for the dominant eigenvector of $C_{ij}$ on fullerenes of the $I_h$-symmetric (a) $C_{20n^2}$ and (b) $C_{60n^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 $I$. 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_{20n^2}$ structure and preserves that of $C_{60n^2}$.</figcaption></figure>

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 transforms under the same irrep of the chiral icosahedral group $I$ 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_{20n^2}$ case, inversion-related sites have opposite spins (in fact, each triangle has a net magnetization, which is odd under inversion), while the $C_{60n^2}$ ansatz is inversion-even. That is, the lowest triplet excitation of all $C_{20n^2}$ fullerenes is expected to transform as $T_{2u}$ (leading to incipient chiral order), while $C_{60n^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 $n = 6$ ($C_{720}$) and $n = 4$ ($C_{960}$), respectively.

Beyond icosahedral symmetry, the ingredients proposed above are also available in the largest achiral subgroup of $I_h$, the pyritohedral group $T_h$, which too has three-dimensional irreps $T_g$ and $T_u$ (corresponding to both $T_1$ and $T_2$ of $I_h$). This is particularly important for crystalline compounds of fullerenes, which cannot maintain full icosahedral symmetry. However, in systems such as the superconductor $K_3C_{60}$ [26], the fullerenes retain full pyritohedral symmetry [54], so our results remain relevant for their behavior. More generally,
noncoplanar ground states may arise even if they are com-
posed 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 generalized a number of techniques
commonly used to diagnose ordering tendencies on finite lat-
tices 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 labeled by point-group irreps rather than wave vectors.

(2) Even in the deep quantum limit $S = 1/2$, these irreps
can be predicted from the ground state of the large-$S$ version
of the Hamiltonian, which can be constructed analogously to
the Luttinger-Tisza method.

(3) The low-energy spectrum is characterized 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 analyze the low-energy
spectrum and ground-state wave function of the nearest-
neighbor spin-1/2 Heisenberg model on the icosahedrally
symmetric C₆₀ and C₈₀ fullerene geometries. Our nu-
merical results were obtained from the variational Monte
Carlo approach using group-convolutional neural-network
wave-function ansätze, which allowed us to construct the
symmetry-resolved low-energy spectrum in detail. We bench-
marked the method against ED on a C₃₂ allotrope and DMRG
on C₆₀, achieving excellent variational energies in both
cases.

For buckminsterfullerene, we found an incipient Bragg
peak transforming under the T<sub>2g</sub> irrep of the icosahedral
point group, which allows for the formation of a Néel pat-
tern on the largest unfrustrated subset of the fullerene graph.
This matches the noncoplanar ground state of the clas-
sical Heisenberg model [18]; furthermore, the low-energy
spectrum of S = 0, 1, 2 excitations follows a tower-of-states
structure derived from a triplet Goldstone-mode excita-
tion 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₈₀ geometry.
The structure of the incipient Bragg peak is again deter-
mined by maximally covering the graph with a Néel pattern,
which leads to a triplet of leading correlation eigenvectors
transforming under the T<sub>2u</sub> irrep. The tower of states cor-
responding 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 func-
tions 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₈₀ 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₈₀ visible in infrared spectroscopy or
vice versa. Furthermore, the noncoplanar magnetic textures of
both molecules may induce an anomalous Hall response [57]
that may in turn affect the superconductivity of such materials
as K₃C₆₀. Here too, the chiral magnetic ground state of C₈₀
may open the door to more exotic superconducting behavior,
which will be interesting to explore in future theoretical and
experimental work.

Finally, we argue that our findings are not limited to the
C₆₀ and C₈₀ geometries, but are relevant for a much wider
class of fullerene geometries. A case in point consists of
the two sequences of fullerenes with full icosahedral sym-
metry, 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. Numerical explorations of
these large molecules, fullerenes with lower symmetry, other
magnetic molecules (e.g., with icosidodecahedral symmetry
and larger spins [58,59]), and the intermediate-*U* Hubbard
model more relevant for real fullerenes are all exciting di-
rections 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₈₀ would also be interesting,
as they may show time-reversal symmetry breaking explic-
itly [60].

ACKNOWLEDGMENTS

We thank Alexander Wietek for helpful discussions. We
are especially grateful to Roman Rausch for providing us
with additional DMRG data for C₆₀ and to Christopher
Roth and Maxime Thumin for their participation in the
early stages of this work. NQS simulations were performed
using the NETKET [42] library. Reference stochastic series
expansion data for the honeycomb lattice were obtained us-
ing the ALPS [48] library. All heat maps use perceptually
uniform color maps developed in Ref. [61]. Computing re-
sources were provided by the STFC Scientific Computing
Department's STFC Cloud service. This work was granted
access to the high-performance computing resources of
the CALMIP center under Allocation No. 2022-P0677 as
well as GENCI (Grant No. A0130500225). This study has
been (partially) supported through EUR NanoX, Grant No.
ANR-17-EURE-0009, in the framework "Programme des In-
vestissements d'Avenir." This work benefited from the support
of the project QMAHT, Grant No. ANR-22-CE30-0032-03,
of the French National Research Agency (ANR). A.S. grate-
fully acknowledges the ISIS Neutron and Muon Source and
the Oxford-ShanghaiTech collaboration for support of the
Keeley-Rutherford fellowship at Wadham College, Oxford.
## 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 G onto a one-dimensional subspace of the symmetry-protected multiplets by imposing additional symmetry constraints, as follows. Consider an Abelian subgroup H of the space group G. Restricting any irrep $\chi$ of G onto H gives a valid representation thereof; for multidimensional irreps, however, this is no longer an irrep of H but can be decomposed into them. If an irrep $\chi_H$ of H appears with multiplicity 1 in this decomposition, we can select a unique representative of $\chi$ by first projecting onto $\chi_H$

<table><caption>TABLE IV. Energies, spin quantum numbers, and point-group irreps of every eigenstate of the $C_{32}$ Heisenberg model below energy -15J from exact diagonalization.</caption><thead><tr><th>Energy</th><th>S</th><th>Energy</th><th>S</th></tr></thead><tbody><tr><td colspan="2">Irrep A'<sub>1</sub></td><td colspan="2">Irrep E'</td></tr><tr><td>-15.7336814</td><td>0</td><td>-15.8119171</td><td>0</td></tr><tr><td>-15.4531726</td><td>1</td><td>-15.6373036</td><td>1</td></tr><tr><td>-15.3587584</td><td>2</td><td>-15.6016680</td><td>1</td></tr><tr><td>-15.3407979</td><td>0</td><td>-15.2579452</td><td>1</td></tr><tr><td>-15.2635666</td><td>0</td><td>-15.2381470</td><td>0</td></tr><tr><td>-15.1351979</td><td>1</td><td>-15.2372144</td><td>2</td></tr><tr><td>-15.1042729</td><td>1</td><td>-15.1730374</td><td>0</td></tr><tr><td>-15.0863465</td><td>0</td><td>-15.1526199</td><td>0</td></tr><tr><td>-15.0376961</td><td>2</td><td>-15.1257156</td><td>1</td></tr><tr><td colspan="2">Irrep A'<sub>2</sub></td><td colspan="2">Irrep E''</td></tr><tr><td>-15.7736580</td><td>1</td><td>-15.0484146</td><td>1</td></tr><tr><td>-15.4543681</td><td>1</td><td>-15.0369711</td><td>1</td></tr><tr><td>-15.1821700</td><td>1</td><td>-15.0329827</td><td>2</td></tr><tr><td>-15.1214539</td><td>1</td><td>-15.0088721</td><td>1</td></tr><tr><td>-15.0740699</td><td>0</td><td colspan="2">Irrep E''</td></tr><tr><td>-15.0340154</td><td>2</td><td>-15.5748511</td><td>0</td></tr><tr><td>-15.0087825</td><td>1</td><td>-15.5658877</td><td>1</td></tr><tr><td colspan="2">Irrep A''<sub>1</sub></td><td colspan="2">Irrep E''</td></tr><tr><td>-15.9372271</td><td>0</td><td>-15.3801966</td><td>1</td></tr><tr><td>-15.4928797</td><td>0</td><td>-15.2871791</td><td>1</td></tr><tr><td>-15.4621861</td><td>0</td><td>-15.1992247</td><td>2</td></tr><tr><td>-15.3572440</td><td>2</td><td>-15.1629586</td><td>1</td></tr><tr><td>-15.3195996</td><td>1</td><td>-15.1445861</td><td>0</td></tr><tr><td>-15.2629763</td><td>1</td><td>-15.0995806</td><td>1</td></tr><tr><td>-15.0493065</td><td>1</td><td>-15.0982993</td><td>0</td></tr><tr><td>-15.0184389</td><td>2</td><td>-15.0617120</td><td>0</td></tr><tr><td>-15.0035477</td><td>1</td><td>-15.0514162</td><td>1</td></tr><tr><td colspan="2">Irrep A''<sub>2</sub></td><td colspan="2">Irrep E''</td></tr><tr><td>-15.5004530</td><td>1</td><td>-15.0397817</td><td>2</td></tr><tr><td>-15.3936259</td><td>1</td><td>-15.0141799</td><td>2</td></tr><tr><td>-15.3823066</td><td>1</td><td colspan="2"></td></tr><tr><td>-15.1635423</td><td>1</td><td colspan="2"></td></tr><tr><td>-15.0894038</td><td>1</td><td colspan="2"></td></tr><tr><td>-15.0325816</td><td>2</td><td colspan="2"></td></tr><tr><td>-15.0287921</td><td>0</td><td colspan="2"></td></tr><tr><td>-15.0202782</td><td>2</td><td colspan="2"></td></tr></tbody></table>

using (5) and then onto $\chi$ itself:

$$ \begin{aligned} |\psi\rangle &= \frac{d_\chi}{|G|} \sum_{g \in G} \chi_g^* \hat{g} \frac{1}{|H|} \sum_{h \in H} \chi_{H,h}^* \hat{h} |\psi_0\rangle \\ &= \frac{d_\chi}{|G||H|} \sum_{g \in G} \left( \sum_{h \in H} \chi_{gh^{-1}} \chi_{H,h} \right)^* \hat{g} |\psi_0\rangle. \end{aligned} \quad (A1) $$

This has the same form as the original projection (5), with the “effective character”

$$ \tilde{\chi}_g = \frac{1}{|H|} \sum_{h \in H} \chi_{gh^{-1}} \chi_{H,h}. \quad (A2) $$

For simplicity, we choose H for each irrep such that their trivial irrep appears with multiplicity 1, so Eq. (A2) reduces to averaging characters in (right) cosets of H.

The $D_{3h}$ point group of $C_{32}$ has two two-dimensional irreps, E' and E''. A suitable choice of H for both is the $C_2$ subgroup generated by one of the three 180° rotations.

The $I_h$ point group of $C_{60}$ and $C_{80}$ has three-, four-, and five-dimensional irreps. We decomposed each of these using the following H:

(i) For T irreps, any $C_3$ rotation subgroup.

(ii) For G irreps, the $C_{2v}$ 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 connecting two nearest-neighbor pentagons.

(iii) For H irreps, any $C_5$ rotation subgroup.

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 $C_5$ rotation group around the center 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 $C_5$: For the plot, we used the real part of the $e^{2\pi i/5}$ rotation-eigenvalue component (that is, a linear combination of the $e^{\pm 2\pi i/5}$ components).

<table><caption>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 $S = 2$ state found by exact diagonalization belongs to the $A'_1$ irrep; however, the variational energy difference between the $A'_1$ and $A''_1$ sectors is smaller than the difference between either and the exact value.</caption><thead><tr><th>Irrep</th><th>$P = +1$</th><th>$P = -1$</th><th>$S^z = 2$</th></tr></thead><tbody><tr><td>$A'_1$</td><td>-15.7279(3)</td><td>-15.4481(3)</td><td>-15.3491(3)</td></tr><tr><td>$A''_2$</td><td>-15.0200(5)<sup>a</sup></td><td><b>-15.7648(4)</b></td><td>-15.0226(4)</td></tr><tr><td>$A''_1$</td><td><b>-15.9342(3)</b></td><td>-15.3053(5)</td><td><b>-15.3537(3)</b></td></tr><tr><td>$A''_2$</td><td>-15.0250(4)<sup>a</sup></td><td>-15.4834(5)</td><td>-15.0264(3)</td></tr><tr><td>$E'$</td><td>-15.8020(4)</td><td>-15.6198(4)</td><td>-15.2220(5)</td></tr><tr><td>$E''$</td><td>-15.5691(3)</td><td>-15.5591(3)</td><td>-15.1786(6)</td></tr></tbody></table>

<footnote><sup>a</sup>$P = +1$ simulation that returned an $S = 2$ state.</footnote>
<table><caption>TABLE VI. Best GCNN variational energies for the C₆₀ geometry, compared with the DMRG variational energies of Ref. [27]. Bold numbers correspond to the lowest energy in each spin sector.</caption><thead><tr><th>Irrep</th><th>P = +1</th><th>P = -1</th><th>S<sup>z</sup> = 2</th></tr></thead><tbody><tr><td>A<sub>g</sub></td><td><b>-31.1302(2)</b></td><td>-29.7548(4)</td><td>-30.2517(2)</td></tr><tr><td>A<sub>u</sub></td><td>-29.6807(6)<sup>a</sup></td><td>-29.8208(7)</td><td>-29.6962(4)</td></tr><tr><td>T<sub>1g</sub></td><td>-29.9447(5)<sup>a</sup></td><td>-30.3186(3)</td><td>-29.9378(6)</td></tr><tr><td>T<sub>1u</sub></td><td>-30.1802(3)</td><td>-30.2944(3)</td><td>-30.0219(5)</td></tr><tr><td>T<sub>2g</sub></td><td>-29.9906(9)<sup>a</sup></td><td><b>-30.7685(4)</b></td><td>-30.0210(4)</td></tr><tr><td>T<sub>2u</sub></td><td>-30.2487(3)</td><td>-30.1286(5)</td><td>-30.0230(6)</td></tr><tr><td>G<sub>g</sub></td><td>-30.3033(5)</td><td>-30.0876(7)</td><td>-30.0850(4)</td></tr><tr><td>G<sub>u</sub></td><td>-30.3551(6)</td><td>-30.6118(5)</td><td>-30.1110(3)</td></tr><tr><td>H<sub>g</sub></td><td>-30.4189(7)</td><td>-30.0871(8)</td><td><b>-30.3251(4)</b></td></tr><tr><td>H<sub>u</sub></td><td>-30.2494(6)<sup>a</sup></td><td>-30.4232(4)</td><td>-30.2619(5)</td></tr><tr><td></td><td>S = 0</td><td>S = 1</td><td>S = 2</td></tr><tr><td>DMRG [27]</td><td>-31.131(7)</td><td>-30.775(6)</td><td>-30.3(2)</td></tr><tr><td></td><td>-30.440(9)</td><td></td><td></td></tr></tbody></table>

<table><caption>TABLE VII. Best GCNN variational energies for the C₈₀ geometry. Bold numbers correspond to the lowest energy in each spin sector.</caption><thead><tr><th>Irrep</th><th>P = +1</th><th>P = -1</th><th>S<sup>z</sup> = 2</th></tr></thead><tbody><tr><td>A<sub>g</sub></td><td>-41.0041(4)</td><td>-40.3017(9)</td><td>-40.7204(5)</td></tr><tr><td>A<sub>u</sub></td><td><b>-41.0387(3)</b></td><td>-40.6590(5)</td><td>-40.6273(5)</td></tr><tr><td>T<sub>1g</sub></td><td>-40.6007(8)<sup>a</sup></td><td>-40.7405(11)</td><td>-40.5850(7)<sup>b</sup></td></tr><tr><td>T<sub>1u</sub></td><td>-40.7063(8)<sup>a</sup></td><td>-40.8184(9)</td><td>-40.7155(7)</td></tr><tr><td>T<sub>2g</sub></td><td>-40.6382(7)<sup>c</sup></td><td><b>-40.9314(7)</b></td><td>-40.5858(8)</td></tr><tr><td>T<sub>2u</sub></td><td>-40.5837(9)<sup>a</sup></td><td>-40.8920(9)</td><td>-40.6155(9)</td></tr><tr><td>G<sub>g</sub></td><td>-40.8578(8)</td><td>-40.7622(9)</td><td>-40.7197(8)</td></tr><tr><td>G<sub>u</sub></td><td>-40.7444(8)</td><td>-40.8820(8)</td><td>-40.7117(8)</td></tr><tr><td>H<sub>g</sub></td><td>-40.8904(12)</td><td>-40.8434(6)</td><td>-40.7264(10)</td></tr><tr><td>H<sub>u</sub></td><td>-40.8350(8)</td><td>-40.8011(9)</td><td><b>-40.7634(8)</b></td></tr></tbody></table>

## APPENDIX B: TABLES OF OPTIMAL VARIATIONAL ENERGIES

We performed exact diagonalization on the C₃₂ Heisenberg model to extend the results in Ref. [20] to every eigenstate

below energy −15*J*; these energies are listed in Table IV. GCNN variational energies in every symmetry sector we considered are listed for C₃₂, C₆₀, and C₈₀ in Tables V, VI, and VII, respectively.

[1] P. W. Anderson, An approximate quantum theory of the antiferromagnetic ground state, Phys. Rev. **86**, 694 (1952).

[2] F. J. Dyson, E. H. Lieb, and B. Simon, Phase transitions in quantum spin systems with isotropic and nonisotropic interactions, J. Stat. Phys. **18**, 335 (1978).

[3] T. Kennedy, E. H. Lieb, and B. S. Shastry, Existence of Néel order in some spin-1/2 Heisenberg antiferromagnets, J. Stat. Phys. **53**, 1019 (1988).

[4] E. J. Neves and J. Perez, Long range order in the ground state of two-dimensional antiferromagnets, Phys. Lett. A **114**, 331 (1986).

[5] A. M. Läuchli, M. Schuler, and A. Wietek, Studying continuous symmetry breaking with exact diagonalization, in *Quantum Materials: Experiments and Theory*, edited by E. Pavarini, E. Koch, J. van den Brink, and G. Sawatzky (Forschungszentrum Jülich GmbH, 2016), Chap. 8.

[6] A. Wietek and A. M. Läuchli, Chiral spin liquid and quantum criticality in extended S = 1/2 Heisenberg models on the triangular lattice, Phys. Rev. B **95**, 035141 (2017).

[7] A. Wietek and A. M. Läuchli, Valence bond solid and possible deconfined quantum criticality in an extended kagome lattice Heisenberg antiferromagnet, Phys. Rev. B **102**, 020411(R) (2020).

[8] A. Wietek, S. Capponi, and A. M. Läuchli, Quantum electrodynamics in 2+1 dimensions as the organizing principle of a triangular lattice antiferromagnet, arXiv:2303.01585.

[9] Y. Nomura and M. Imada, Dirac-type nodal spin liquid revealed by refined quantum many-body solver using neural-network wave function, correlation ratio, and level spectroscopy, Phys. Rev. X **11**, 031034 (2021).

[10] D. Wulferding, Y. Choi, W. Lee, and K.-Y. Choi, Raman spectroscopic diagnostic of quantum spin liquids, J. Phys.: Condens. Matter **32**, 043001 (2020).

[11] J. Schnack, Effects of frustration on magnetic molecules: a survey from Olivier Kahn until today, Dalton Trans. **39**, 4677 (2010).

[12] I. Rousochatzakis, A. M. Läuchli, and F. Mila, Highly frustrated magnetic clusters: The kagomé on a sphere, Phys. Rev. B **77**, 094420 (2008).

[13] A. Furrer and O. Waldmann, Magnetic cluster excitations, Rev. Mod. Phys. **85**, 367 (2013).

[14] A. J. Stace and P. O'Brien, Fullerenes: past, present and future, celebrating the 30th anniversary of Buckminster Fullerene, Philos. Trans. R. Soc. London A **374**, 20160278 (2016).

[15] F. Lin, E. S. Sørensen, C. Kallin, and A. J. Berlinsky, Strong correlation effects in the fullerene C₂₀ studied using a one-band Hubbard model, Phys. Rev. B **76**, 033414 (2007).

[16] N. P. Konstantinidis, Antiferromagnetic Heisenberg model on clusters with icosahedral symmetry, Phys. Rev. B **72**, 064453 (2005).

[17] S. Chakravarty, M. P. Gelfand, and S. Kivelson, Electronic correlation effects and superconductivity in doped fullerenes, Science **254**, 970 (1991).

[18] D. Coffey and S. A. Trugman, Magnetic properties of undoped C₆₀, Phys. Rev. Lett. **69**, 176 (1992).

[19] R. T. Scalettar, A. Moreo, E. Dagotto, L. Bergomi, T. Jolicoeur, and H. Monien, Ground-state properties of the Hubbard model on a C₆₀ cluster, Phys. Rev. B **47**, 12316 (1993).

[20] N. P. Konstantinidis, s = 1/2 antiferromagnetic Heisenberg model on fullerene-type symmetry clusters, Phys. Rev. B **80**, 134427 (2009).

<footnote><sup>a</sup>P = +1 simulation that returned an S = 2 state.</footnote>

<footnote><sup>a</sup>P = +1 simulation that returned an S = 2 state.</footnote>

<footnote><sup>b</sup>S<sup>z</sup> = 2 simulation that returned an S = 3 state.</footnote>

<footnote><sup>c</sup>State that is not clearly dominated by one spin-quantum-number sector.</footnote>
[21] N. A. Modine and E. Kaxiras, Variational Hilbert-space-truncation approach to quantum Heisenberg antiferromagnets on frustrated clusters, Phys. Rev. B **53**, 2546 (1996).

[22] N. Konstantinidis, Zero-temperature magnetic response of small fullerene molecules at the classical and full quantum limit, J. Magn. Magn. Mater. **449**, 55 (2018).

[23] K. Momma and F. Izumi, VESTA 3 for three-dimensional visualization of crystal, volumetric and morphology data, J. Appl. Crystallogr. **44**, 1272 (2011).

[24] D. Tománek, Guide Through the Nanocarbon Jungle (Morgan & Claypool, San Rafael, CA, 2014).

[25] S. C. Benjamin, A. Ardavan, G. A. D. Briggs, D. A. Britz, D. Gunlycke, J. Jefferson, M. A. G. Jones, D. F. Leigh, B. W. Lovett, A. N. Khlobystov, S. A. Lyon, J. J. L. Morton, K. Porfyrakis, M. R. Sambrook, and A. M. Tyrishkin, Towards a fullerene-based quantum computer, J. Phys.: Condens. Matter **18**, S867 (2006).

[26] A. F. Hebard, M. J. Rosseinsky, R. C. Haddon, D. W. Murphy, S. H. Glarum, T. T. M. Palstra, A. P. Ramirez, and A. R. Kortan, Superconductivity at 18 K in potassium-doped $C_{60}$, Nature (London) **350**, 600 (1991).

[27] R. Rausch, C. Plorin, and M. Peschke, The antiferromagnetic $S = 1/2$ Heisenberg model on the $C_{60}$ fullerene geometry, SciPost Phys. **10**, 087 (2021).

[28] C. Roth and A. H. MacDonald, Group convolutional neural networks improve quantum state accuracy, arXiv:2104.05085.

[29] C. Roth, A. Szabó, and A. H. MacDonald, High-accuracy variational Monte Carlo for frustrated magnets with deep neural networks, Phys. Rev. B **108**, 054410 (2023).

[30] B. Bernu, P. Lecheminant, C. Lhuillier, and L. Pierre, Exact spectra, spin susceptibilities, and order parameter of the quantum Heisenberg antiferromagnet on the triangular lattice, Phys. Rev. B **50**, 10048 (1994).

[31] J. D. Reger, J. A. Riera, and A. P. Young, Monte Carlo simulations of the spin-1/2 Heisenberg antiferromagnet in two dimensions, J. Phys.: Condens. Matter **1**, 1855 (1989).

[32] E. V. Castro, N. M. R. Peres, K. S. D. Beach, and A. W. Sandvik, Site dilution of quantum spins in the honeycomb lattice, Phys. Rev. B **73**, 054422 (2006).

[33] In a lattice geometry, $C_{ij}$ is translation invariant, and so these eigenvectors are plane waves, recovering the usual momentum-space treatment.

[34] Unlike an infinite-lattice model, the Goldstone-mode operators (3) have a finite support on each site, and so they do not perfectly commute. However, the norm of the commutator is inversely proportional to the system size, and so we do not expect it to give rise to a low-lying “antisymmetric two-Goldstone” excitation. That is, we can treat the Goldstone modes as commuting, bosonic operators.

[35] J. M. Luttinger and L. Tisza, Theory of dipole interaction in crystals, Phys. Rev. **70**, 954 (1946).

[36] D. H. Lyons and T. A. Kaplan, Method for determining ground-state spin configurations, Phys. Rev. **120**, 1580 (1960).

[37] T. A. Kaplan and N. Menyuk, Spin ordering in three-dimensional crystals with strong competing exchange interactions, Philos. Mag. **87**, 3711 (2007).

[38] T. S. Cohen and M. Welling, Group equivariant convolutional networks, arXiv:1602.07576.

[39] G. Carleo and M. Troyer, Solving the quantum many-body problem with artificial neural networks, Science **355**, 602 (2017).

[40] J. Carrasquilla and G. Torlai, How to use neural networks to investigate quantum many-body physics, PRX Quantum **2**, 040201 (2021).

[41] V. Heine, Group Theory in Quantum Mechanics, International Series in Natural Philosophy Vol. 91 (Pergamon, Oxford, 1960).

[42] F. Vicentini, D. Hofmann, A. Szabó, D. Wu, C. Roth, C. Giuliani, G. Pescia, J. Nys, V. Vargas-Calderón, N. Astrakhantsev, and G. Carleo, NetKet 3: Machine learning toolbox for many-body quantum systems, SciPost Phys. Codebases **7** (2022).

[43] M. Reh, M. Schmitt, and M. Gärttner, Optimizing design choices for neural quantum states, Phys. Rev. B **107**, 195115 (2023).

[44] G. Klambauer, T. Unterthiner, A. Mayr, and S. Hochreiter, Self-normalizing neural networks, in Advances in Neural Information Processing Systems, Vol. 30, edited by I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (Curran Associates, Red Hook, NY, 2017), pp. 971–980.

[45] This coincidence of 11 states at nearby energies may explain the large error bars of the corresponding DMRG result, together with the smaller bond dimension (compared with the $S^z = 0$ states) used in the calculation of Ref. [27].

[46] A. W. Sandvik and J. Kurkijärvi, Quantum Monte Carlo simulation method for spin systems, Phys. Rev. B **43**, 5950 (1991).

[47] F. Alet, S. Wessel, and M. Troyer, Generalized directed loop method for quantum Monte Carlo simulations, Phys. Rev. E **71**, 036706 (2005).

[48] B. Bauer, L. D. Carr, H. G. Evertz, A. Feiguin, J. Freire, S. Fuchs, L. Gamper, J. Gukelberger, E. Gull, S. Guertler, A. Hehn, R. Igarashi, S. V. Isakov, D. Koop, P. N. Ma, P. Mates, H. Matsuo, O. Parcollet, G. Pawłowski, J. D. Picon et al., The ALPS project release 2.0: open source software for strongly correlated systems, J. Stat. Mech. (2011) P05001.

[49] R. Rausch, M. Peschke, C. Plorin, and C. Karrasch, Magnetic properties of a capped kagome molecule with 60 quantum spins, SciPost Phys. **12**, 143 (2022).

[50] The Goldstone-operator construction also predicts $A_u$, $A_g$ quin-tets. These are indeed at a similar energy to the $H_u$, $H_g$ ones but are not as clearly separate from the rest of the $S = 2$ spectrum as in the $C_{60}$ case.

[51] This can be seen without detailed calculation: The product of three inversion-odd operators is inversion-odd, and the only inversion-odd one-dimensional irrep of $I_h$ is $A_u$.

[52] N. P. Konstantinidis, Competition between frustration and spin dimensionality in the classical antiferromagnetic n-vector model with arbitrary n, SciPost Phys. Core **6**, 042 (2023).

[53] R. B. King and M. V. Diudea, The chirality of icosahedral fullerenes: a comparison of the tripling (leapfrog), quadrupling (chamfering), and septupling (capra) transformations, J. Math. Chem. **39**, 597 (2006).

[54] K. M. Allen, W. I. F. David, J. M. Fox, R. M. Ibberson, and M. J. Rosseinsky, Molecular structure of the fulleride anions in superconducting $K_3C_{60}$ and insulating $K_6C_{60}$ determined by powder neutron diffraction, Chem. Mater. **7**, 764 (1995).
[55] An example is the $C_{70}$ rugbyballene molecule, whose Luttinger-Tisza ground state (in the $A'_2$ irrep of the $D_{5h}$ point group) vanishes on 20 sites near the ends of the molecule, requiring an additional $E'_2$ component to form a valid classical ground state. However, this molecule is not inversion symmetric, and so no chiral state in the strict sense may form on it.

[56] P. Bernath, *Spectra of Atoms and Molecules* (Oxford University Press, Oxford, 2005).

[57] T. Yamaguchi and A. Yamakage, Theory of magnetic-texture-induced anomalous Hall effect on the surface of topological insulators, *J. Phys. Soc. Jpn.* **90**, 063703 (2021).

[58] E. Neuscamman and G. K.-L. Chan, Correlator product state study of molecular magnetism in the giant Keplerate $Mo_{72}Fe_{30}$, *Phys. Rev. B* **86**, 064402 (2012).

[59] J. Ummethum, J. Schnack, and A. M. Läuchli, Large-scale numerical investigations of the antiferromagnetic Heisenberg icosidodecahedron, *J. Magn. Magn. Mater.* **327**, 103 (2013).

[60] J. Cirac, D. Pérez-García, N. Schuch, and F. Verstraete, Matrix product density operators: Renormalization fixed points and boundary theories, *Ann. Phys. (Amsterdam)* **378**, 100 (2017).

[61] P. Kovesi, Good colour maps: How to design them, arXiv:1509.03700.