Evidence for a Dirac Spin Liquid in the Kagome Heisenberg Antiferromagnet from High-Precision Neural-Network Quantum States

manuscript.pdf

Introduction

Introduction.—Quantum spin liquids (QSLs) represent a fundamental class of quantum phases beyond the Landau symmetry-breaking paradigm, exhibiting fractionalized excitations, long-ranged quantum entanglement, and emergent gauge fields [1, 2, 3, 4]. Understanding the nature of such phases of matter has become one of the central questions in condensed-matter physics. The long-standing experimental search for QSLs has placed kagome antiferromagnets at the forefront of two-dimensional quantum magnetism, with prominent candidate materials including herbertsmithite [5, 6, 7, 8], Zn-barlowite [9, 10, 11], and Y-based kagome antiferromagnets [12, 13, 14]. The spin-1/2 kagome Heisenberg model with nearest-neighbour interactions has become the minimal canonical model for kagome antiferromagnets, yet the nature of its ground state remains elusive. Early density-matrix renormalization-group (DMRG) calculations provided evidence for a $Z_2$ gapped QSL [15, 16, 17]. However, subsequent calculations based on infinite projected entangled simplex states found scaling behaviour consistent with a gapless state [18], and flux-threaded infinite-cylinder DMRG revealed Dirac-cone signatures [19]. Alongside tensor-network simulations, variational Monte Carlo (VMC) studies based on Gutzwiller-projected fermionic mean-field ansätze initially identified the $\pi$-flux $U(1)$ Dirac QSL as a leading variational candidate [20, 21, 22, 23, 24], although further optimizations have also revealed closely competing $Z_2$ states [25].

Despite the aforementioned numerical progress, both DMRG and VMC face intrinsic limitations. DMRG is largely restricted to quasi-one-dimensional cylinders, with a cost that grows rapidly with circumference, whereas VMC is constrained by the bias and limited expressivity of the underlying mean-field ansatz.

Neural-network quantum states (NQSs) provide a variational framework that can be directly applied to large two-dimensional lattices with periodic boundary conditions, thereby allowing studies on torus rather than quasi-one-dimensional cylinders. Moreover, NQSs can represent complex many-body correlations and sign structures without being constrained to a specific parton mean-field manifold [26, 27, 28]. Recent NQS-based studies have achieved ground-state energy estimates competitive with the best available results for $J_1$–$J_2$ square-lattice and triangular Heisenberg models [29, 30, 31, 32, 33], as well as two-dimensional Fermi-Hubbard models [34, 35, 36, 37].

The group convolutional neural-network (G-CNN) architecture of NQS has been applied to the kagome Heisenberg antiferromagnet, suggesting an evidence of spinon pair-density-wave (SPDW) ground state [38]. However, a subsequent analysis reveals that the unusually low variational energies reported in Ref. [38] are artifacts induced by nonergodic Monte Carlo sampling [39], and the signatures of SPDW require a further demonstration. These deficiencies leave the NQS study of the kagome Heisenberg ground state unresolved and demand a rigorous reassessment based on reliable sampling.

In this work, we study ground-state properties of spin-1/2 kagome Heisenberg antiferromagnet using deep neural quantum states whose architecture is tailored to the kagome geometry. Our variational ansatz extends the residual convolutional network architecture developed in Ref. [30, 31, 32], i.e., the ResNet, to the kagome model by replacing conventional grid-based kernels with message-passing convolutions that propagate features along the nearest-neighbour bonds while retaining an explicit representation of the three sublattices in each unit cell (see Fig. 1a-c). We benchmark that the ResNet with geometry-aware message-passing convolutions (KagomeResNet) outperforms its standard-convolution counterparts. To the best of our knowledge, the resulting variational energies for large system sizes on a torus achieve the lowest raw energies, apart from the disputed G-CNN results of Ref. [38]. More importantly, based on the analysis of static structure factors, correlation functions, and projections of the optimized NQS onto manifolds of Gutzwiller-projected mean-field states, we provide consistent evidence for a Dirac QSL ground state in the spin-1/2 kagome Heisenberg antiferromagnet.

Model and variational wave function. The Hamiltonian of kagome Heisenberg antiferromagnet studied in this work is

with $\mathbf{S}_{i} = (S_{i}^{x}, S_{i}^{y}, S_{i}^{z})$ being the spin-1/2 operators, and $\langle i,j \rangle$ indicating pairs of nearest-neighbor sites on a kagome lattice shown in Fig. 1a. Here, we consider the periodic boundary condition of the model, i.e., the torus, and the system sizes $L = 6$, 8 and 10, with number of sites $N = 108$, 192 and 300, respectively ($N = L \times L \times3$).

In the NQS approach of training the ground state of $\mathcal{H}$, a neural network provides a compact parametrization of the many-body wave function. For a spin-1/2 system with $N$ sites, the Hilbert space is spanned by the $z$-basis spin configurations $\lvert\mathbf{s} \rangle= \lvert s_{1}, \ldots, s_{N} \rangle$ with $s_{i} = \pm1$. The NQS parameterized by $\theta$ can map each spin configuration $\mathbf{s}$ to a complex wave-function amplitude $\psi_{\theta}(\mathbf{s})$. To provide a phase prior for the ground state of Hamiltonian (1), we incorporate a fixed reference phase corresponding to a coplanar $q = 0$ three-sublattice pattern $\phi_{q=0}(\mathbf{s})$ into the wave-function parametrization, thereby defining the variational state $\lvert\Psi_{\theta} \rangle= \sum_{\mathbf{s}} \psi_{\theta}(\mathbf{s})\phi_{q=0}(\mathbf{s})\lvert\mathbf{s} \rangle$. Within variational Monte Carlo, the parameters $\theta$ are optimized using MinSR [31], to minimize the variational energy $E_{\theta} = \langle\Psi_{\theta} \lvert\mathcal{H} \rvert\Psi_{\theta} \rangle/ \langle\Psi_{\theta} \vert\Psi_{\theta} \rangle$. We emphasize that the complex-valued network $\psi_{\theta}(\mathbf{s})$ can learn additional phase corrections beyond the reference factor $\phi_{q=0}(\mathbf{s})$. More details of the optimization are provided in the Supplemental Material.

The architecture of our NQS, the KagomeResNet, is plotted in Fig. 1a-c. First, the spin configuration $s_{i}$ is reshaped into three $L \times L$ matrices associated with the $A$, $B$, and $C$ sublattices of the kagome lattice. Then, as a key operation for our NQS, we perform message passing directly along the nearest-neighbor kagome connectivity. Specifically, in Fig. 1b, we show the operation for an $A$-sublattice site as an example, i.e.,

where $x_{\alpha,i}$ and $y_{\alpha,i}$ are the input and output feature vectors at site $i$ of sublattice $\alpha$, respectively. The matrices $S_A$ and $W_{1,\ldots,4}$ represent trainable transformations for the central site and its four directed nearest-neighbor bonds. Similar operations are performed for the $B$ and $C$ sublattices, and thus there are 12 nearest-neighbor message-passing channels in a complete kagome convolutional (KagomeConv) layer. The full structure of KagomeResNet is shown in Fig. 1c. The $l$-th residual block acts on $x_l \in\mathbb{R}^{3 \times L \times L \times c}$, where $c$ is the channel dimension, and maps it to $x_{l+1}$ via two KagomeConv layers with GELU activations and a residual connection. We provide more details for the KagomeResNet and conventional ResNet in the Supplementary Material. The network contains $n$ residual blocks, and the output is finally converted by normalization and a pair-complex activation into the wave-function amplitude $\psi_\theta(s)$.

Benchmarks.—We now benchmark the performance of KagomeResNet for training the ground state of the Hamiltonian (1) with $L=6$ in the $\Gamma$-point sector. We first perform $10^4$ optimization steps imposing only translational symmetry, and subsequently for 1800 steps including the even spin-inversion sector $\mathcal{I}_s=+1$ and $D_6$ point-group symmetry. In Fig. 1d, for the KagomeResNet with $c=n=32$, we show the final 1000 iterations of the translationally symmetric optimization, followed by the full 1800-step optimization after the additional symmetry projections are imposed, and demonstrate that the converged value of the ground-state energy is $E=-47.1597(6)$, which is highly competitive with the previously reported raw variational energies [16, 22, 23], and also lower than the conventional ResNet comparable number of variational parameters (see Supplemental Material).

Recent work using a G-CNN ansatz for the kagome model (1) suggests that non-ergodic Monte Carlo sampling can induce low energy during the optimization, and reports $E \simeq-45.6$ for G-CNN with ergodic sampling [39]. In the inset of Fig. 1d, we show the acceptance rate $R_{\mathrm{acc}} \simeq0.555$ over the final 500 iterations with full symmetries, demonstrating that our Monte Carlo sampling is ergodic. The acceptance rate is monitored throughout all subsequent optimizations in this work to ensure that the Monte Carlo sampling remains ergodic.

We then benchmark the low-lying energies in the odd spin-inversion sector $\mathcal{I}_s=-1$. For $L=4$, the exact diagonalization study shows that the lowest-energy with $\mathcal{I}_s=-1$ lies in the $M$-point sector, while the lowest $\Gamma$-point energy is nearly degenerate, only $\sim6.8 \times10^{-3}$ higher in energy [40]. In Fig. 1e, we plot the final 1000 iterations of the translationally symmetric optimization with $\mathcal{I}_s=-1$ for both $M$ and $\Gamma$-point sectors, followed by 1800 additional iterations imposing their respective little-group symmetries $D_6$ and $C_{2v}$. Our results show that, different from the $L=4$ case, the lowest-energy state within the odd spin-inversion sector $\mathcal{I}_s=-1$ occurs at the $\Gamma$ point for $L=6$. More importantly, upon comparing both spin-inversion sectors, we find the overall ground state at the $\Gamma$ point with $\mathcal{I}_s=+1$. This restores the same ground-state quantum numbers as found by exact diagonalization for $L=4$,

in stark contrast to Ref. [38], which reported a change of the ground-state sector to the $M$ point with $\mathcal{I}_s=-1$ for $L=6$. We therefore find no evidence for the size-driven change of ground-state quantum numbers that underlies the SPDW interpretation of Ref. [38].

$L$	\multicolumn{2}{c}{$c=n=32$}	\multicolumn{2}{c}{$c=n=40$}
	$e$	$\sigma^2$	$e$	$\sigma^2$
6	$-0.436664(6)$	$0.00350(4)$	$-0.436760(5)$	$0.00312(3)$
8	$-0.435835(4)$	$0.00519(5)$	$-0.435949(5)$	$0.00444(4)$
10	$-0.435239(4)$	$0.00506(6)$	$-0.435472(5)$	$0.00434(4)$

In Table I (PDF p. 19), we systematically represent the energy density $e=E/N$ and the energy density variance $\sigma^2$ obtained from KagomeResNet with $n=c=32$ and $n=c=40$ for the kagome model (1) with $L=6$, 8 and 10 at the $\Gamma$ point with $\mathcal{I}_s=+1$. For the $L=8$ torus, previous references of the raw energy density include $e=-0.4286749(4)$ and $-0.428783(1)$ for the projected mean-field VMC based on $U(1)$ Dirac [22, 23] and $Z_2$ [25] QSLs, respectively, and $e \sim-0.4349$ for the symmetric projected entangled-pair states calculations [41]. Our variational energies are substantially lower than the aforementioned benchmarks, improving the best known estimate for the energy density $e$ to the third decimal place.

$L$	$\mu Z_2^*$	$\rho$	$\eta$	$\phi_1$	$\phi_2$	$\phi_3$
6	0.4884	0.3126	0.3428	2.0412	2.8204	2.8421
8	0.5534	0.3577	0.3920	2.5210	2.9475	2.9692
10	0.5739	0.3788	0.4137	2.6851	3.0009	3.0246

Ground-state correlation functions.—To study the ground-state properties of the Hamiltonian (1), for the trained NQS with $n=c=40$, we first calculate the spin structure factors

with $\mathbf{q}=(q_x,q_y)$ the wave vector in the Cartesian reciprocal-space coordinates, and the dimer structure factors

with the sum running over pairs of nearest-neighbor bonds, $\mathbf{r}_{m k}=(\mathbf{r}_m+\mathbf{r}_k)/2$ being the midpoint of bond $(m,k)$, and $D_{m k,n l}=\langle B_{m k}B_{n l}\rangle-\langle B_{m k}\rangle\langle B_{n l}\rangle$, where $B_{m k}=S_m\cdot S_k$ is the nearest-neighbor bond-energy operator. In Fig. 2a, we plot the spin structure factor $S_f^{zz}(\mathbf{q})$ in an extended Brillouin zone. We show that no sharp peak develops inside the first Brillouin zone, where the maxima occur at the reciprocal-lattice images $G$ of $\Gamma$. Since $G$ differs from $\Gamma$ only by a reciprocal-lattice vector, the corresponding real-space correlations characterize a zero-Bravais-momentum three-sublattice fluctuation pattern repeated across unit cells, rather than a distinct finite-momentum modulation. The nontrivial structure resides within the unit cell, where the $A$, $B$, and $C$ sublattice components form predominantly the nonuniform $\Gamma-E_2$ sector through the kagome sublattice form factor (see Supplemental Material). As shown in Fig. 2b, the normalized peak values of spin structure factors $S_{\max}^{zz}/N$ tend to 0 with $N\to\infty$, indicating the absence of long-range magnetic order. We also provide the data of $S^{D}_{f}(\mathbf{q})$ and its normalized maxima $S^{D}_{\max}$ in the Supplementary Material, showing that $S^{D}_{\max}/N\to0$ with increasing $N$, which rules out the valence-bond-solid (VBS) order. Additionally, the dramatically large value of $S^{D}_{\max}\sim87.37$ for $L=6$ reported by Ref. [38] is not reproduced in our calculations. Instead, $S^{D}_{\max}$ systematically decreases with increasing $L$. Thus, we argue that the signatures of an SPDW ground state reported for Hamiltonian (1) in Ref. [38] require further validation, in particular through explicit checks of the ergodicity of the Monte Carlo sampling.

The dimer structure factors (4) consider a particular bond form factor fixed by the bond-center positions. To resolve the full bond-correlation structure, we then focus on the connected bond-correlation matrix

where $B_{\alpha}(\mathbf{R})$ is the nearest-neighbor bond-energy operator, $\alpha=1,\ldots,6$ labels the complete set of six nearest-neighbor bonds associated with each Bravais unit cell, and $\delta B_{\alpha}(\mathbf{R})=B_{\alpha}(\mathbf{R})-\langle B_{\alpha}(\mathbf{R})\rangle$. In Fig. 2c, we plot the largest eigenvalue of the matrix $K^{\mathrm{conn}}_{\alpha,\beta}(\mathbf{q})$, i.e., $\lambda_1[K^{\mathrm{conn}}(\mathbf{q})]$, which quantifies the strongest connected dimer fluctuation at the momentum $\mathbf{q}$. It is seen that $\lambda_1[K^{\mathrm{conn}}(\mathbf{q})]$ is maximal at the $M$ points. In real space, the dominant $M$-point dimer fluctuations correspond to period-doubled bond modulations, with the bond pattern alternating under translations along one or both primitive Bravais-lattice directions. The three $M$ points represent the three symmetry-related orientations of this modulation. Condensation in the $M$-point dimer channel implies broken lattice-translation symmetry, while $\lambda_1[K^{\mathrm{conn}}(M)]/N\to0$ with increasing $N$ precludes long-range VBS order (see Fig. 2b).

Next, we resolve the $M$-point dimer fluctuations into symmetry-adapted sectors of the kagome space group, including $F_2$, $B_2$, and $F_1$ (see Supplemental Material for their explicit construction and symmetry projectors). As shown in the inset of Fig. 2d, the $F_2$ sector is dominant for all $L$. Together with Fig. 2a, this identifies the $M$-$F_2$ dimer and $\Gamma$-$E_2$ spin channels as the dominant symmetry-resolved fluctuations. To further examine the relation between these two channels, we compare the real-space correlations in the aforementioned two dominant symmetry channels. We define the sublattice-resolved spin and bond correlations as

where $\alpha,\beta=A,B,C$ for the spin correlations, while $\alpha,\beta=1,\ldots,6$ label the nearest-neighbor bonds for the dimer correlations defined above. The $\Gamma$-$E_2$ spin correlation is obtained by projecting the three-sublattice correlation matrix onto the $E_2$ sector, i.e., $C_{\Gamma-E_2}^{\mathrm{zz}}(\mathbf{r})=\frac{1}{2}\sum_{\alpha,\beta=A,B,C}(P_{E_2})_{\alpha\beta}C_{\alpha\beta}^{\mathrm{zz}}(\mathbf{r})$ with $P_{E_2}=I-\mathbf{u}\mathbf{u}^{\dagger}$ and $\mathbf{u}=\frac{1}{\sqrt{3}}(1,1,1)^{T}$. Similarly, we define $C_{F_2}^{D}(\mathbf{r})=\frac{1}{3}\sum_{i=1}^{3}e^{-i\mathbf{M}_{i}\cdot\mathbf{r}}[\mathbf{v}_{F_2}^{(i)}]^{\dagger}C^{D}(\mathbf{r})\mathbf{v}_{F_2}^{(i)}$, where the explicit $F_2$ bond vectors and their symmetry construction are given in the Supplemental Material.

In the $N_f=4$ QED$_3$ description of the $U(1)$ Dirac QSL of the kagome Heisenberg antiferromagnet, the $M$-$F_2$ dimer and $\Gamma$-$E_2$ spin channels carry the quantum numbers of spin-singlet and spin-triplet minimal monopoles, respectively [42]. These monopoles are components of the same emergent $SO(6)$ multiplet and thus share a common scaling dimension [42, 43]. In Fig. 2d, we directly compare $C_{F_2}^{D}(r)$ and $C_{\Gamma-E_2}^{\mathrm{zz}}(r)$ for different distances $r$ and system sizes $L$. Remarkably, all data collapse onto a proportional function $C_{F_2}^{D}(r)=A C_{\Gamma-E_2}^{\mathrm{zz}}(r)$ with $A\simeq2.52$. This is precisely the behavior expected for different components of the same emergent $SO(6)$ monopole multiplet, showing that the symmetry-resolved correlations obtained from our NQS calculations are consistent with the emergent monopole structure of the kagome $U(1)$ Dirac QSL.

Overlap of ground state with candidate QSLs.

Previous studies of projected mean-field VMC for the kagome model (1) cannot unambiguously distinguish between the between $Z_2$ and the U(1) Dirac spin-liquid states because their variational energies are nearly degenerate [21, 25]. Our NQS is not biased toward either of the competing parton mean-field ansätze, providing an independent and more accurate representation of the ground-state wave function of the Hamiltonian (1). This motivates us to assess whether the optimized NQS is closer to the U(1) Dirac or $Z_2$ spin-liquid manifold.

For this assessment, we use the following procedure. We determine the fidelity between our optimized NQS $\lvert\Psi_\theta^{\mathrm{opt}}\rangle$ with representative Gutzwiller-projected fermionic states from these competing variational families $\lvert\Phi_{\mu,\nu}\rangle$ with $\mu$ and $\nu$ labeling the variational family and spinon boundary sector, respectively. The fidelity is given by

and can be evaluated by Monte Carlo sampling configurations from the NQS probability distribution. We adopt approximately $3 \times10^5$ Monte Carlo samples to obtain fidelity estimates with small statistical uncertainties. As the Dirac spin-liquid reference ($\mu= \mathrm{Dirac}$), we employ the optimized extended U(1) Dirac state in Ref. [21], as one of the lowest projected-mean-field variational energies reported for the Dirac state. For $Z_2$ candidates, we consider two complementary references. The first ($\mu= Z_2$) is obtained by introducing a finite perturbation away from the Dirac state considered in Ref. [21], without producing any detectable change in the variational energy. The second ($\mu= Z_2^*$) is the fully optimized $Z_2$ resonating-valence-bond state in Ref. [25], whose variational energy lies slightly below that of their optimized U(1) state. Following the boundary conditions used in the original VMC studies, we restrict $\nu$ to the two mixed periodic-antiperiodic spinon sectors. For the Dirac state, only nondegenerate closed-shell sectors are considered, however for the $Z_2$ states, both mixed sectors are retained. Finally, we compare the maximum fidelity $F_\mu^{\max}$ over all admissible boundary sectors for different $\mu$. More details of the parametrizations and boundary conditions for the chosen mean-field states are given in the Supplemental Material.

The results of $F_\mu^{\max}$ for the aforementioned three Gutzwiller-projected mean-field states and the optimized KagomeResNet NQS with $n = c = 40$ are displayed in Fig. 3. Although the variational energies for the Dirac and $Z_2$ ($\mu= \mathrm{Dirac}$ and $Z_2$) reference states in Ref. [21] are indistinguishable, $F_{Z_2}^{\max}$ is significantly smaller than $F_{\mathrm{Dirac}}^{\max}$. Moreover, for a more competitive $Z_2$ reference ($\mu= Z_2^*$) in Ref. [25] with lower variational energy, $F_{Z_2^*}^{\max}$ is close to $F_{\mathrm{Dirac}}^{\max}$ at $L = 6$, while it becomes systematically smaller than $F_{\mathrm{Dirac}}^{\max}$ with increasing $L$. The relative separation is highlighted by the ratio $F_{\mathrm{Dirac}}^{\max}/F_{Z_2^*}^{\max}$, which increases linearly with $L$ (see the inset of Fig. 3a).

We emphasize that although both the $Z_2$ and $Z_2^*$ reference states originate from a gapped $Z_2$ mean-field ansätze, the latter, after Gutzwiller projection, exhibits properties resembling those of a Dirac QSL, such as a gapless spin excitation spectrum [25]. The much larger fidelity $F_{Z_2^*}^{\max} \gg F_{Z_2}^{\max}$ shows

that the optimized NQS is substantially closer to this reference with Dirac-like spin properties. Combined with the direct fidelity preference for the Dirac reference at $L = 8$ and $10$ (Fig. 3a), these results support a Dirac spin-liquid interpretation of the NQS ground state.

Discussion. — We have studied the ground-state properties of spin-1/2 kagome Heisenberg antiferromagnet with nearest-neighbor interactions using the KagomeResNet, an NQS architecture developed in this work, featuring message-passing convolutions tailored to the geometry of kagome lattice. In stark contrast to the results in Ref. [38], we show that the ground state lies in the $\Gamma$ instead of $M$ point and no evidence of SPDWs states can be identified. The analysis of symmetry-resolved spin and dimer correlations and the wave-function fidelity comparisons support the U(1) Dirac QSL nature for the ground state of the kagome model (1).

Future work can extend the application of KagomeResNet to explore the ground-state properties of other kagome systems, such as the long-range $XY$ kagome model, which can be realized by dipolar Rydberg atoms [44, 45, 46]. In addition, we have demonstrated that the KagomeResNet outperforms the conventional ResNet, suggesting that the geometry-aware message-passing convolutions can play a key role in the NQS architecture. This can motivate the design of NQS for solving the ground states of models with more complex geometry, including the spin-1/2 Heisenberg antiferromagnet on the shuriken lattice [47, 48] and maple-leaf lattice [49].

Acknowledgments

The NQS simulations are performed using Quantax [50] based on JAX [51] and Equinox [52]. We gratefully acknowledge the scientific support and high-performance computing resources provided by the Erlangen National High Performance Computing Center (NHR) of the Friedrich-Alexander-Universität Erlangen-Nürnberg (NHR Project No. nqsQuMat), and by the Jülich Booster (JUPITER Booster) at the Jülich Supercomputing Centre, Forschungszentrum Jülich (project neuquass). This work was supported by the German Research Foundation (DFG) via project 492547816 (TRR 360). C.R. acknowledges support from the Stewart Blusson Quantum Matter Institute at the University of British Columbia.

References

Supplementary Materials for “Evidence for a Dirac Spin Liquid in the Kagome Heisenberg Antiferromagnet from High-Precision Neural-Network Quantum States”

Zheng-Hang Sun; Christopher Roth; Ao Chen; Markus Heyl

supp.pdf

Supplementary Materials for “Evidence for a Dirac Spin Liquid in the Kagome Heisenberg Antiferromagnet from High-Precision Neural-Network Quantum States”

Details of network architecture and results for conventional residual convolutional network

In this section, we present the details of the conventional residual convolutional network ResNet, which is employed as a baseline for comparison with KagomeResNet considered in the main text.

As shown in Fig. 1a, the kagome lattice is regarded as a Bravais lattice with three sublattice sites $A$, $B$, and $C$ in each unit cell. For an $L \times L$ system, the total number of spins is $N = 3L^2$, and a spin configuration can be written as $\mathbf{s} = \{s_{\alpha}(i,j)\}$ with $\alpha= A, B, C$, $s_{\alpha}(i,j) = \pm1$, and $(i,j)$ labelling the unit cell along the two primitive vectors $\mathbf{a}_1$ and $\mathbf{a}_2$. The conventional ResNet treats the three spins within each unit cell as a single local degree of freedom. Specifically, we first introduce the binary variables $b_{\alpha}(i,j) = \frac{s_{\alpha}(i,j)+1}{2} \in\{0,1\}$ and encode the local spin configuration by the integer

with $k(i,j) \in\{0,\ldots,7\}$. Therefore, the $2^3 = 8$ possible spin configurations of a kagome unit cell are represented as eight local states.

A trainable embedding matrix $E \in\mathbb{R}^{c \times8}$ then maps each local state to a $c$-dimensional feature vector,

The resulting input to the convolutional network therefore has shape $X^{(0)} \in\mathbb{R}^{c \times L \times L}$. After the embedding layer, all spatial transformations are performed by standard two-dimensional convolutions on the Bravais-cell grid, as illustrated in Fig. 1b. For a hidden feature tensor $X \in\mathbb{R}^{c \times L \times L}$, a convolutional layer can be written as

where

is the weight matrix associated with the cell displacement $\boldsymbol{\delta} = (\delta_1,\delta_2)$. For periodic boundary conditions, the cell indices in Eq. (3) are understood modulo $L$, corresponding to the circular padding used in the implementation. The $3 \times3$ stencil contains nine Bravais-cell positions, including the central cell $\boldsymbol{\delta} = (0,0)$. On the Bravais lattice, six of the eight surrounding cells, $(\pm1,0)$, $(0,\pm1)$, $(1,-1)$, and $(-1,1)$ are nearest-neighbour Bravais cells, whereas $(1,1)$ and $(-1,-1)$ are next-nearest-neighbour Bravais cells. Hence, the standard $3 \times3$ convolution treats these cell-index displacements as independent learnable kernel components.

Figure Fig. 1c shows the entire ResNet architecture. Starting from the embedded tensor $X_0 = X^{(0)}$, the network contains $n$ residual blocks. Each block consists of two standard convolutions with GELU nonlinearities and an identity shortcut. For the $l$th residual block ($l = 0,\ldots,n-1$), the transformation implemented in our calculations is

Here, the “Norm” operation shown in Fig. 1c is therefore a deterministic rescaling by $1/\sqrt{l+1}$. All residual blocks preserve both the number of feature channels $c$ and the spatial dimensions $L \times L$. Both convolutions in each residual block contain trainable biases, except for the second convolution in the final residual block, for which the bias is omitted.

After the final residual block, the feature tensor is rescaled once more,

The variational parameters of the convolutional network are real. When a complex-valued wave function is required, the $c$ real output channels are divided into two equal sets and combined pairwise according to

which corresponds to the “Pair complex” operation in Fig. 1c. We therefore choose an even number of channels $c$ whenever a complex output is used. The complex features are subsequently passed through the default exponential output activation,

which is evaluated in the code using a numerically stabilized scaled representation to avoid overflow. Finally, the output is aggregated over the feature channels and symmetrized with respect to lattice translations. For the default translationally invariant sector, the final layer implements the zero-momentum projection, yielding the variational wave-function amplitude $\psi(\mathbf{s})$.

We note that the KagomeResNet considered in the main text follows the same overall residual architecture as the conventional ResNet shown in Fig. 1c. The two differences lie in the input representation and the convolution operator. In the following, we provide more details of KagomeResNet.

The first difference from the conventional ResNet concerns the representation of the spin configuration. As discussed above, the conventional ResNet groups the three spins within each kagome unit cell into one of eight local states and subsequently maps this state to a $c$-dimensional trainable embedding. In contrast, as shown in the Fig. 1a of the main text, KagomeResNet does not employ such an embedding layer. Instead, the microscopic kagome structure is retained explicitly throughout the network. For a kagome lattice with $N = 3L^2$ spins, the input spin configuration $\mathbf{s} = \{s_{\alpha}(i,j)\}$ with $\alpha= A,B,C$, is directly reshaped according to

with $X^{(0)}_{1,\alpha,i,j} = s_{\alpha}(i,j)$. Here, the first dimension denotes the feature-channel index, while the second dimension explicitly labels the three kagome sublattices $A$, $B$, and $C$. The remaining two indices $(i,j)$ specify the Bravais unit cell. The reshape operation contains no trainable parameters and preserves the original spin variables $s_{\alpha}(i,j) = \pm1$ without coarse graining the three sites of a unit cell into a single local state. Consequently, whereas the conventional ResNet represents the kagome lattice by a tensor of shape $(c,L,L)$ after the embedding layer, KagomeResNet starts from a tensor of shape $(1,3,L,L)$, and retains the explicit sublattice dimension throughout all convolutional layers. The feature-channel dimension is expanded from 1 to $c$ only by the first kagome convolution.

The second, and more important, difference is the convolution operator. As shown in the Fig. 1b of the main text, instead of applying a standard $3 \times3$ convolution on the Bravais-cell grid, KagomeResNet employs a geometry-aware convolution that performs message passing directly along the microscopic nearest-neighbour bonds of the kagome lattice. For an input tensor $X \in\mathbb{R}^{c_{\mathrm{in}}\times3\times L\times L}$, the kagome convolution produces $Y \in\mathbb{R}^{c_{\mathrm{out}}\times3\times L\times L}$, while preserving the explicit $A$, $B$, and $C$ sublattice axis. The nearest-neighbour connectivity is encoded by a fixed set of directed bond classes. Using the sublattice labels $(A,B,C) \equiv(0,1,2)$ and the Bravais-cell displacement $\boldsymbol{\delta} = (\delta_1,\delta_2)$, the six undirected nearest-neighbour bond classes used in the implementation are

Each tuple $(\alpha,\beta,\delta_1,\delta_2)$ denotes a bond connecting sublattice $\alpha$ in one unit cell to sublattice $\beta$ in a unit cell displaced by $\boldsymbol{\delta} = \delta_1\mathbf{a}_1 + \delta_2\mathbf{a}_2$. Both directions of every bond are included in the message-passing operation, resulting in twelve directed nearest-neighbour channels per unit cell. For every directed bond $e = (\alpha\to\beta,\boldsymbol{\delta})$, KagomeResNet associates a learnable matrix

which transforms the feature vector on the source site before it is accumulated on the target site. For a directed bond class $e = (\alpha,\beta,\boldsymbol{\delta})$, with $\boldsymbol{\delta} = \delta_1\mathbf{a}_1 + \delta_2\mathbf{a}_2$, the contribution from sublattice $\alpha$ in unit cell $\mathbf{r}$ to sublattice $\beta$ in unit cell $\mathbf{r}+\boldsymbol{\delta}$ is defined as

where $W_e \in\mathbb{R}^{c_{\mathrm{out}}\times c_{\mathrm{in}}}$ is the trainable weight matrix associated with the directed bond class $e$. The output feature on each target site is obtained by summing the contributions from all directed nearest-neighbour bond classes terminating at that site, together with the corresponding on-site contribution. The same matrix $W_e$ is shared over all lattice translations of a given directed bond class, thereby preserving translational equivariance while explicitly encoding the microscopic connectivity of the kagome lattice.

In Fig. 1d, we directly compare the energy trajectory of the final 1,800 optimization iterations after imposing $\mathcal{I}_s = +1$ and the $D_6$ point-group symmetry for both the conventional ResNet and KagomeResNet. Here, we consider the conventional ResNet with $n = 32$, $c = 42$ and $3 \times3$ kernel, whose number of trainable parameters is $N_{\mathrm{para}} = 1019046$, which is slightly larger than $N_{\mathrm{para}} = 974208$ for the KagomeResNet with $n = c = 32$ considered in our main text. It is seen that KagomeResNet outperforms conventional ResNet.

Details of variational optimization

For both the conventional ResNet and the KagomeResNet, we supplement the complex network output $\psi_{\theta}(\mathbf{s})$ with a fixed three-sublattice phase factor associated with a coplanar pattern,

where $s_j = \pm1$ and $\alpha(j) \in\{A, B, C\}$ labels the sublattice of site $j$. Each down spin contributes the phase associated with its sublattice. This factor has unit modulus and is held fixed throughout training. The complex-valued network $\psi_\theta(s)$ can learn additional configuration-dependent phase corrections. Thus, the reference phase provides a prior without fixing the full phase structure of the variational state.

We initially optimize the neural-network quantum state (NQS) with translational invariance imposed. After $10^4$ optimization steps, we additionally impose spin-flip symmetry and the spatial symmetries appropriate to the target symmetry sector, and continue training from the previously optimized network parameters. For the optimization method, we adopt AdamSR in Quantax, which combines MinSR [1] with Adam-type momentum and adaptive rescaling. We adopt the training rate is $10^{-3}$.

Additional data for the spin and dimer structure factors

In the main text, we present the spin structure factor for $L = 10$ as a function of the wave vector $\mathbf{q} = (q_x,q_y)$ expressed in Cartesian reciprocal-space coordinates (see Fig. 2a of the main text). For a direct comparison with Ref. [2], we additionally present the structure factors using the reciprocal-lattice coordinates $(q_1,q_2)$ employed therein. We emphasize that $q_1$ and $q_2$ are not Cartesian components of the momentum, but dimensionless coefficients with respect to the two primitive reciprocal-lattice vectors.

Using the real-space primitive vectors $\mathbf{a}_1 = (1,0)$, $\mathbf{a}_2 = \left(\frac{1}{2},\frac{\sqrt{3}}{2}\right)$, the corresponding reciprocal-lattice vectors are $\mathbf{b}_1 = 2\pi\left(1,-\frac{1}{\sqrt{3}}\right)$, $\mathbf{b}_2 = 2\pi\left(0,\frac{2}{\sqrt{3}}\right)$, which satisfy $\mathbf{a}_i \cdot\mathbf{b}_j = 2\pi\delta_{ij}$. The two momentum representations are therefore related by $\mathbf{q} = q_1\mathbf{b}_1 + q_2\mathbf{b}_2$. For an $L \times L$ torus, the momenta in the first reciprocal unit cell are quantized as

Momentum points in an extended-zone representation are obtained by reciprocal lattice translations,

which corresponds in Cartesian coordinates to $\mathbf{q}\to\mathbf{q}+m_1\mathbf{b}_1+m_2\mathbf{b}_2$.

Figure 2 shows the corresponding spin and dimer structure factors for $L=6$, $8$, and $10$ in the $(q_1,q_2)$ representation. For $L=6$, Ref. [2] reported a dramatic enhancement of $S_f^D(\mathbf{q})$ relative to the $L=4$ result. In particular, while the maximum value for $L=4$ was approximately 0.933, the reported maximum increased to $S_{\max}^D\simeq87.37$ for $N=108$, accompanied by numerous pronounced peaks throughout the extended reciprocal space. This strong enhancement was interpreted in Ref. [2] as a signature of the proposed spinon pair-density-wave state and its associated secondary orders. We do not observe such behavior. Already for the same $L=6$ ($N=108$) torus (Fig. 2d), the dimer structure factor remains of order unity and exhibits no anomalously enhanced peaks. Upon increasing the system size to $L=8$ and $L=10$ (Figs. 2e and f), no growth analogous to that reported in Ref. [2] develops. Instead, the maximum value $S_{\max}^D$ decreases with increasing $L$, while $S_{\max}^D/N$ extrapolates toward zero, as shown in Fig. 2b of the main text. The absence of any size-enhanced dimer peak is therefore inconsistent with the anomalously large dimer response reported in Ref. [2] and provides no support for the associated SPDW-induced dimer ordering tendency.

In the Fig. 2 of the main text, we present the data of spin structure factors $S_f^{zz}(\mathbf{q})$ and largest eigenvalue of the connected bond-correction matrix $\lambda_1[K^{\mathrm{conn}}(\mathbf{q})]$ in Cartesian reciprocal-space coordinates $\mathbf{q}=(q_x,q_y)$ with $L=10$. Here, in Fig. 3, we also display the results of $S_f^{zz}(\mathbf{q})$ and $\lambda_1[K^{\mathrm{conn}}(\mathbf{q})]$ ($\mathbf{q}=(q_x,q_y)$) with other two system sizes, i.e., $L=6$ and $8$ For all three system sizes, the maxima of $S_f^{zz}(\mathbf{q})$ occur at the reciprocal-lattice images $G$ of the $\Gamma$ point, whereas the maxima of $\lambda_1[K^{\mathrm{conn}}(\mathbf{q})]$ are located at the $M$ points. The locations of these dominant spin and dimer fluctuations are therefore robust against changes in system size.

Symmetry-resolved spin and dimer correlations

In this section, we provide the symmetry construction underlying the $\Gamma$-$E_2$ spin and $M$-$F_2$ dimer correlations discussed in Fig. 2 of the main text. The purpose of this analysis is to go beyond the momentum dependence of the scalar structure factors and resolve the internal sublattice or bond structure of the dominant fluctuations. As we show below, the spin maxima at reciprocal-lattice images of $\Gamma$ predominantly probe the nonuniform $E_2$ sublattice sector, whereas the dominant dimer fluctuations at the three $M$ points belong to the $F_2$ bond sector. These are precisely the symmetry channels that can be compared with the monopole operator content of the $N_f=4$ QED$_3$ description of the kagome Dirac spin liquid (DSL).

-E2 spin channel

We first consider the spin correlations. For sublattices $\alpha,\beta=A,B,C$, we define the $3\times3$ real-space correlation matrix

where $\mathbf{R}$ labels the Bravais cells and $N=3L^2$. Correspondingly, the sublattice-resolved momentum-space correlation matrix is obtained by Fourier transforming $C_{\alpha\beta}^{zz}(\mathbf{r})$.

The physical spin structure factor defined in the main text contains, in addition, the kagome sublattice form factor. Denoting the position of sublattice $\alpha$ inside the unit cell by $\boldsymbol{\rho}_\alpha$, it can be written as

where $\overline{\mathbf{q}}$ denotes the momentum folded back into the first Brillouin zone. This expression also explains why the physical structure factor of a non-Bravais lattice is not generally periodic under $\mathbf{q} \to\mathbf{q}+\mathbf{G}$: although $\overline{\mathbf{q}}$ is unchanged, the sublattice form factors $e^{i\mathbf{G}\cdot\boldsymbol{\rho}_{\alpha}}$ need not be uniform.

At the $\Gamma$ point, the three-dimensional sublattice space decomposes into the uniform one-dimensional component and its two-dimensional orthogonal complement. The normalized uniform vector is

with projector

The complementary two-dimensional subspace transforms as the $E_2$ representation and is selected by

For example, an orthonormal basis of this subspace can be chosen as

such that $P_{E_2}=\mathbf{e}_1\mathbf{e}_1^{\dagger}+\mathbf{e}_2\mathbf{e}_2^{\dagger}$.

We therefore define the symmetry-resolved $\Gamma$-$E_2$ correlation used in the main text as

The factor $1/2$ averages over the two components of the $E_2$ representation. In contrast, contraction with $P_{A_1}$ gives the uniform $\Gamma$-point spin channel. Equation (27) provides the connection to Fig. 2a of the main text. At a reciprocal-lattice image $\mathbf{q}=\mathbf{G}$ of $\Gamma$, the Bravais momentum is again $\Gamma$, but the vector of sublattice form factors

contains a nonuniform component in the $E_2$ subspace. Consequently, the pronounced $\mathbf{G}$-point maxima of $S_f^{zz}(\mathbf{q})$ probe predominantly the nonuniform $\Gamma$-$E_2$ fluctuations rather than a distinct finite-momentum magnetic mode. The persistence of the same momentum and sublattice structure for $L=6$, $8$, and $10$, shown in Figs. 2 and 3 of the Supplemental Material and Fig. 2 of the main text, identifies $\Gamma$-$E_2$ as the dominant symmetry-resolved spin channel.

This identification also has a direct field-theory interpretation. For the kagome $U(1)$ DSL, Ref. [3] shows that the spin-triplet minimal monopoles carry zero crystal momentum and angular momentum $l=2$ under the sixfold lattice rotation. On the kagome lattice these quantum numbers correspond precisely to the $\Gamma$-$E_2$ sector. Thus, $C_{\Gamma-E_2}^{zz}$ provides a microscopic lattice observable sensitive to the spin-triplet monopole channel of the DSL.

M-point bond space and symmetry projectors

We then resolve the dimer correlations. Since a kagome unit cell contains six nearest-neighbour bonds, retaining the complete internal bond structure requires a six-dimensional bond basis. We use the same bond convention in the measurements and throughout the analysis:

We denote $\delta B_{\alpha}(\mathbf{R}) = B_{\alpha}(\mathbf{R}) - \langle B_{\alpha}\rangle$ and collect these operators into the vector $\delta\mathbf{B}(\mathbf{R}) = (\delta B_1,\ldots,\delta B_6)^{\mathsf{T}}$. The connected real-space bond-correlation matrix is

whose Fourier transform gives the connected matrix $K_{\alpha\beta}^{\mathrm{conn}}(\mathbf{q})$ defined in the main text.

The strongest eigenvalue of $K^{\mathrm{conn}}(\mathbf{q})$ occurs at the three symmetry-related $M$ points. In the reduced reciprocal-lattice coordinates introduced in Sec. II, we choose $M_1 = (\frac{1}{2},0)$, $M_2 = (0,\frac{1}{2})$, and $M_3 = (\frac{1}{2},\frac{1}{2})$. The three momenta form the $M$ star under the $D_6$ point group.

At an individual $M_i$, the little co-group is $C_{2v}$. We take its generators to be a twofold rotation $C_2$ about a hexagon center and a reflection $\tau_i$ whose axis leaves $M_i$ invariant. A lattice symmetry acts not only by permuting the six bonds but, whenever a bond is mapped to a neighbouring reference unit cell, also generates the corresponding Bloch phase at $M_i$. We denote the resulting six-dimensional bond-space representation by $D_{M_i}(g)$.

For a one-dimensional irreducible representation $\lambda$ of the little group, the symmetry projector is

where $\chi_{\lambda}(g)$ is the corresponding character. Equivalently, labeling a one-dimensional sector by its eigenvalues $\eta_2=\pm1$ and $\eta_{\tau}=\pm1$ under $C_2$ and $\tau_i$, respectively, Eq. (36) becomes

This construction includes the momentum-dependent phases and therefore classifies the full microscopic bond patterns, rather than only their orientations inside one unit cell.

For our six-bond basis, the representation at a representative $M$ point decomposes as

In the $M$-star notation used in the main text, the $A_1$ and $A_2$ sectors generate the $F_1$ and $F_2$ bond modes, respectively. The leading eigenmode of $K^{\mathrm{conn}}(M_i)$ belongs to the one-dimensional, mirror-odd $A_2$ sector, i.e. the $F_2$ sector. The next modes are classified analogously by applying the projectors in Eq. (37); the dominant sectors displayed in the inset of Fig. 2d of the main text are $F_2$, $B_2$, and $F_1$.

Explicit $F_2$ bond vectors

With the bond ordering in Eq. (34), a normalized choice of the $F_2$ vectors at the three members of the $M$ star is

The three vectors are related by the $C_6$ lattice rotations that permute the members of the $M$ star, including the Bloch phases associated with bonds crossing the unit-cell boundary. Their overall signs are conventional and have no physical consequence, since all projected correlations involve the rank-one projectors

The $F_2$ contribution to the connected bond covariance at a given $M_i$ is therefore

Numerically, the leading eigenvector of $K^{\mathrm{conn}}(M_i)$ has essentially unit overlap with this symmetry sector for all three system sizes. The ordering of the leading symmetry-resolved eigenvalues, $F_2 > B_2 > F_1$, also remains unchanged from $L = 6$ to $L = 10$. Thus, the identification of the dominant $M$-point fluctuation with $F_2$ is not a consequence of a particular system size or of an unconstrained numerical eigenvector, but follows from a fixed symmetry projector.

To compare this finite-momentum dimer fluctuation directly with the $\Gamma$-point spin channel, we remove the oscillating $M$-point phase and average over the three members of the $M$ star. This gives the real-space correlation used in the main text,

The symmetry analysis above identifies the microscopic structure underlying Fig. 2c,d of the main text. The dominant bond fluctuations are centered at the $M$ star and exhibit an internal six-bond pattern that lies overwhelmingly in the mirror-odd $F_2$ sector.

Connection to the QED3 monopole multiplet

Next, we discuss why the two symmetry channels singled out above are particularly relevant for the Dirac-spin-liquid interpretation. The low-energy theory of the kagome $U(1)$ DSL is $N_f = 4$ QED$_3$. The six lowest charge-one monopole operators transform in the antisymmetric representation $6$ of SU(4), or equivalently as a vector of the emergent SO(6) symmetry. They decompose into three spin-singlet and three spin-triplet components. Their lattice quantum numbers on the kagome lattice were clarified in Ref. [3]. In particular, the spin-triplet monopoles have zero crystal momentum and angular momentum $l = 2$ under $C_6$, identifying them with the $\Gamma$-$E_2$ spin sector considered above.

The spin-singlet monopoles instead carry the translation quantum numbers of the three $M$ points and correspond to VBS fluctuations. Matching the translation and little-group transformation laws in the Table 3 of Ref. [3] to the bond-space representation constructed above identifies the mirror-odd component with the $M$-$F_2$ sector. This observation is useful because the symmetry resolution separates the physical dimer response into operators with different continuum quantum numbers, which are mixed together in the scalar dimer structure factor $S_f^D(q)$.

A central consequence of the emergent SO(6) structure is that all six components of the minimal-monopole multiplet have the same scaling dimension $\Delta_\Phi$. In the continuum theory their leading two-point function takes the form

as discussed in Refs. [3, 4]. The microscopic spin and dimer operators generally have different nonuniversal overlaps with their corresponding continuum monopole operators. Consequently, their long-distance correlations need not have the same amplitude, but they are expected to share the same scaling dependence,

up to subleading lattice and finite-size corrections. It follows that in the regime controlled by the common monopole operator,

where $A$ is nonuniversal. This is the motivation for the comparison in Fig. 2d of the main text. Across different separations and the three system sizes $L = 6$, $8$, and $10$, the two independently constructed symmetry-resolved correlations collapse approximately onto a common proportional relation.

Additional data for spin and dimer correlation functions

For comparison with the symmetry-resolved channels, here we define the raw real-space spin correlation by averaging over all physical site pairs belonging to the same periodic-boundary-condition distance shell,

where $N_r^{(s)}$ denotes the number of site pairs in the shell at distance $r$. We also define the raw connected dimer correlation by averaging over all physical nearest-neighbour bond pairs with the same periodic-boundary-condition distance between their bond midpoints,

where $r_b$ denotes the midpoint of bond $b$, and $N_r^{(b)}$ is the number of bond pairs in the corresponding distance shell.

The raw spin and dimer correlations are shown in Figure 4a and d, respectively. Both quantities exhibit pronounced oscillations between positive and negative values as a function of distance $r$, together with an overall decrease in magnitude toward longer distances. Such behavior reflects the superposition of different internal sublattice or bond-symmetry components in the symmetry-unresolved correlators, and the raw correlations therefore do not exhibit a clear scaling form.

To further examine the scaling behavior of the correlation functions, we now consider the symmetry-resolved correlation functions $C^{zz}_{\Gamma-E_2}(r)$ and $C^{D}_{F_2}(r)$ defined by Eq. (32) and (42), respectively. The results of $C^{zz}_{\Gamma-E_2}(r)$ for $L = 10$ are shown in Figure 4b and c, with the power-law, i.e., Eq. (44), and exponential fits, respectively. We emphasize that for the fitting analysis, we restrict ourselves to the axial $r = 2, 4, 6, 8$ and $10$, which correspond to displacements by one to five Bravais unit cells along the principal lattice directions. Since the primitive Bravais vectors have length 2 in units of the nearest-neighbor spacing, these are precisely the nonzero axial separations available on the $L = 10$ torus up to the half-system distance. Restricting the fits to these axial shells allows us to compare the distance dependence along a fixed crystallographic direction and avoids mixing inequivalent lattice directions, which can exhibit different finite-size and short-distance corrections.

As shown in Figure 4b and Figure 4c, the power-law form provides the better fit. Similarly, based on the results of $C_{F_2}^{D}(r)$ shown in Figure 4e and Figure 4f, the power-law form describes the dimer correlation function better than the exponential form within the same fitting window. For the power-law fits, the decay rates of $C_{\Gamma-E_2}^{Z_2}(r)$ and $C_{F_2}^{D}(r)$ are $2\Delta_{\Phi}=3.13$ and $2.98$, respectively, which are close to one another.

Gutzwiller-projected spin-liquid reference states and wave-function fidelity

In this section, we provide the explicit construction of the three Gutzwiller-projected fermionic reference states used in Fig. 3 of the main text and describe the evaluation of their fidelities with the optimized KagomeResNet wave functions. Following the notation of the main text, the three reference states are labeled by

Here, $\mu=\mathrm{Dirac}$ denotes the optimized extended $U(1)$ Dirac state of Ref. [5]; $\mu=Z_2$ denotes a finite-pairing $Z_2[0,\pi]\beta$ state obtained by perturbing the same Dirac state as considered in Ref. [5]; and $\mu=Z_2^*$ denotes the fully optimized $Z_2[0,\pi]\beta$ resonating-valence-bond (RVB) state of Ref. [6]. The auxiliary-spinon boundary condition is denoted by $\nu$, so that the corresponding projected reference state is written as $\lvert\Phi_{\mu,\nu}\rangle$.

Unified fermionic parent Hamiltonian

We use the fermionic parton representation

subject to the local single-occupancy constraint

All three reference states in Eq. (48) can be generated from the same general real singlet-BCS parent Hamiltonian,

where $t_{ij}=t_{ji}$ and $\Delta_{ij}=\Delta_{ji}$ are real in the gauges considered below. The case $i=j$ is included, allowing for on-site particle-hole and pairing terms.

Introducing the Nambu spinor

the parent Hamiltonian may equivalently be expressed in terms of matrices

in the $\tau_1$–$\tau_3$ plane of the SU(2) gauge space.

For a given reference family $\mu$ and auxiliary-spinon boundary condition $\nu$, we denote the mean-field ground state of Eq. (51) by $\lvert\Psi_{\mathrm{MF}}^{(\mu,\nu)}\rangle$. The physical reference state entering the fidelity calculation is

where

projects onto exactly one spinon per kagome site.

Equation (44) provides a common parent-Hamiltonian framework for all three reference states. The specific hopping and pairing amplitudes, however, must be implemented in the gauge convention of the corresponding literature Ansatz. In particular, the parametrization of Ref. [6] should not be obtained by directly identifying its parameters with those of Ref. [5]; the two Ansätze use different bond sign conventions and the former additionally contains third-neighbour RVB amplitudes.

Dirac reference state: $\mu= \mathrm{Dirac}$

For $\mu= \mathrm{Dirac}$, we use the optimized extended $U(1)$ Dirac state of Ref. [5]. In the gauge used therein, the NN and second-neighbour mean-field matrices are

where $s_{ij}=\pm1$ and $\nu_{ij}^{I}=\pm1$ follow the bond-sign convention of Fig. 1 of Ref. [5]. The Ansatz has zero gauge flux through the elementary triangles and $\pi$ flux through the hexagons.

The NN hopping defines the overall energy scale,

The optimized extended Dirac state reported in Ref. [5] has

where $\mu_{I}$ denotes the uniform on-site chemical potential of the mean-field Hamiltonian. We use the subscript “I” here to distinguish this parameter from the reference-state label $\mu$ employed throughout the main text.

Because the state is number conserving and the projected wave function is constructed at fixed half filling, the uniform chemical-potential term does not change the occupied single-particle orbitals. The projected reference state is therefore fully determined by the hopping Ansatz in Eq. (56) at half filling.

Ref. [5] performed the variational calculation on a 432-site cluster with mixed periodic-antiperiodic spinon boundary conditions and obtained

Finite-pairing $Z_{2}$ reference state: $\mu= Z_{2}$

The reference state denoted by $\mu= Z_{2}$ in the main text is a finite-pairing $Z_{2}[0,\pi]_{\beta}$ state in the immediate variational neighbourhood of the Dirac state of Ref. [5]. Its mean-field matrices are

A nonzero $\Delta_{2}$ generates noncommuting $\mathrm{SU}(2)$ fluxes and breaks the emergent $U(1)$ gauge structure of the Dirac Ansatz down to $Z_{2}$.

The specific finite-pairing reference used in Fig. 3 is defined by

For $\Delta_{2}=0.05$, Ref. [5] reports the corresponding optimized on-site pairing $\zeta=0.1780(2)$. Within the statistical resolution of the VMC calculation, introducing this finite pairing produces no detectable change in the variational energy relative to the Dirac state. Therefore, on the same $N=432$ cluster,

where Ref. [5] does not quote a separate independent energy value for this particular finite-pairing point.

We stress that the state defined by Eq. (61) is not the fully optimized minimum of the $Z_2[0,\pi]\beta$ variational family. Upon complete stochastic-reconfiguration optimization, Ref. [5] finds $\Delta_2,\zeta\to0$, such that the optimized state returns to the $U(1)$ Dirac state. We retain the finite-pairing state as an independent reference in order to test whether the optimized NQS shows an enhanced affinity to an explicit $U(1) \to Z_2$ pairing perturbation whose energy is nearly indistinguishable from that of the Dirac state.

Optimized $Z_2$ reference state: $\mu= Z_2^*$

For the second $Z_2$ reference, denoted by $\mu= Z_2^*$, we use the fully optimized $Z_2[0,\pi]\beta$ RVB state of Ref. [6]. This Ansatz is defined in a different gauge from that of Ref. [5] and contains RVB amplitudes up to third-neighbour bonds. Following Ref. [6], we introduce

such that

With the NN amplitude fixed to unity, the parent matrices are

Here $s_{ij}^{Z}=\pm1$ follows the six-site spinon-unit-cell gauge shown in Fig. 1 of Ref. [6]. The spinon unit cell is doubled along $\mathbf{a}_1$ and therefore contains six kagome sites.

In terms of the hopping and pairing amplitudes of the unified parent Hamiltonian in Eq. (51), the second- and third-neighbour terms are

Thus, $\eta$ is the magnitude of the third-neighbour RVB parameter and should not be identified with the on-site pairing $\zeta$ in Eq. (60). Similarly, $\rho$ cannot be identified directly with the Iqbal parameter $\chi_2$ without an explicit gauge transformation of the complete bond-sign structure.

For the three system sizes studied in the present work, we use the finite-size optimized parameters shown in Fig. 3 of Ref. [6]. The values used in our numerical implementation are The optimized parameters are presented graphically rather than tabulated numerically in Ref. [6]; the values in Table I (PDF p. 19) were therefore obtained by digitizing the published Fig. 3. The retained digits are used only to define the numerical reference states reproducibly and should not be interpreted as table-level precision of the original publication.

The corresponding variational energies reported directly in Table I (PDF p. 19) of Ref. [6] are. The $L=18$ cluster is the largest size for which Ref. [6] reports well-converged optimized RVB parameters. At this size, the corresponding optimized $U(1)$ state has energy density $-0.428719(1)J$, while the optimized $Z_2$ state has $-0.428737(1)J$. Thus, within this projected-RVB variational family, the $Z_2$ state lies slightly below the optimized $U(1)$ state in energy.

$L$	$N = 3L^2$	$E_{Z_2^*}/N$
6	108	$-0.428899(1)J$
8	192	$-0.428783(1)J$
10	300	$-0.428753(1)J$
18	972	$-0.428737(1)J$

$L$	PA	AP
6	open shell	closed shell
8	closed shell	closed shell
10	open shell	closed shell

The three reference states used in Fig. 3 are therefore complementary. The state $\mu= \mathrm{Dirac}$ provides an optimized extended $U(1)$ Dirac reference; $\mu= Z_2$ probes an explicit finite-pairing perturbation in the immediate neighbourhood of this Dirac state; and $\mu= Z_2^*$ provides an independently optimized and energetically competitive $Z_2[0,\pi]\beta$ reference with additional third-neighbour variational freedom.

Auxiliary-spinon boundary conditions

The physical kagome Heisenberg Hamiltonian and the KagomeResNet NQS are always defined with periodic boundary conditions. The boundary conditions discussed here apply only to the auxiliary fermions used to construct the Gutzwiller-projected reference states.

For an $L \times L$ torus, we define

with $\theta_{1,2}=0$ or $\pi$. We denote periodic and antiperiodic boundary conditions by P and A, respectively. Following the mixed periodic-antiperiodic boundary conditions used in the original VMC calculations Refs. [5, 6], the fidelity comparison in the main text is restricted to

Here, the first and second letters refer to the spinon boundary conditions along $\mathbf a_1$ and $\mathbf a_2$, respectively. For $\mu=\mathrm{Dirac}$, the number-conserving mean-field state is used only when the half-filled spinon spectrum forms a nondegenerate closed shell. For the system sizes considered here, the two mixed sectors satisfy An open-shell sector is excluded because the degenerate single-particle level at the Fermi energy does not define a unique half-filled Slater determinant without an additional prescription for the zero-mode occupations.

For the paired references $\mu=Z_2$ and $\mu=Z_2^*$, the BdG ground state is well defined in both mixed sectors, and both PA and AP are retained. Consequently, the quantity reported in Fig. 3 of the main text is

where $\mathcal B_\mu(L)$ denotes the admissible set of mixed boundary sectors for reference state $\mu$. For $\mu=\mathrm{Dirac}$, $\mathcal B_\mu(L)$ contains only nondegenerate closed-shell sectors; for $\mu=Z_2$ and $Z_2^*$, both mixed sectors are included.

Projected amplitudes and Monte Carlo evaluation of the fidelity

For the number-conserving reference $\mu=\mathrm{Dirac}$, the mean-field ground state is obtained by filling the lowest $N/2$ single-particle orbitals for each spin species. For a spin configuration $\mathbf s$ with sets of up- and down-spin sites $\mathcal U$ and $\mathcal D$, respectively, the Gutzwiller-projected amplitude can be expressed as the product of the corresponding Slater determinants,

where $(-1)^{P(s)}$ accounts for the fermionic reordering required to express the projected state in the spin-configuration basis used for the NQS.

For $\mu= Z_2$ and $Z_2^*$, diagonalization of the corresponding BdG Hamiltonian yields the BCS pair matrix $F^{(\mu,\nu)}$. After projection to $N_{\uparrow}=N_{\downarrow}=N/2$, the amplitude is

Thus, all three projected references and the optimized NQS are evaluated in exactly the same physical spin basis.

For each admissible pair $(\mu,\nu)$, the normalized squared many-body fidelity is

This is the same quantity defined in the main text.

We evaluate Eq. (72) using configurations sampled from the optimized NQS probability distribution

For each reference state, we define the amplitude ratio

It then follows that

and

Therefore,

Neither the NQS nor the projected reference state therefore needs to be normalized explicitly.

The fidelity provides information complementary to the symmetry-resolved correlation functions discussed in the preceding section. Whereas the latter probe selected long-distance operators, $F_{\mu,\nu}$ compares the complete many-body amplitude and phase structure of the optimized NQS with a given projected spin-liquid reference state.

As shown in Fig. 3 of the main text, the finite-pairing state $\mu=Z_2$ has a substantially smaller maximum fidelity than $\mu=\mathrm{Dirac}$, despite the two reference states having indistinguishable variational energies at the resolution reported in Ref. [5]. The more competitive reference $\mu=Z_2^*$ has a fidelity close to that of the Dirac state for $L=6$, but its $F_{Z_2^*}^{\max}$ becomes systematically smaller than $F_{\mathrm{Dirac}}^{\max}$ with increasing system size. Importantly, the $Z_2^*$ reference was found in Ref. [6] to have a slightly lower projected-VMC energy than the corresponding optimized $U(1)$ state. The fidelity comparison is therefore not simply tracking the variational-energy ordering within the projected-fermion family, but provides an independent wave-function-level diagnostic supporting the Dirac interpretation of the optimized KagomeResNet ground state.

References