Dirac-type nodal spin liquid revealed by refined quantum many-body solver using neural-network wave function, correlation ratio, and level spectroscopy

Pursuing fractionalized particles that do not bear properties of conventional measurable objects, exemplified by bare particles in the vacuum such as electrons and elementary excitations such as magnons, is a challenge in physics. Here we show that a machine-learning method for quantum many-body systems that has achieved state-of-the-art accuracy reveals the existence of a quantum spin liquid (QSL) phase in the region $0.49\lesssim J_2/J_1\lesssim0.54$ convincingly in spin-1/2 frustrated Heisenberg model with the nearest and next-nearest neighbor exchanges, $J_1$ and $J_2$, respectively, on the square lattice. This is achieved by combining with the cutting-edge computational schemes known as the correlation ratio and level spectroscopy methods to mitigate the finite-size effects. The quantitative one-to-one correspondence between the correlations in the ground state and the excitation spectra enables the reliable identification and estimation of the QSL and its nature. The spin excitation spectra containing both singlet and triplet gapless Dirac-like dispersions signal the emergence of gapless fractionalized spin-1/2 Dirac-type spinons in the distinctive QSL phase. Unexplored critical behavior with coexisting and dual power-law decays of N\'{e}el antiferromagnetic and dimer correlations is revealed. The power-law decay exponents of the two correlations differently vary with $J_2/J_1$ in the QSL phase and thus have different values except for a single point satisfying the symmetry of the two correlations. The isomorph of excitations with the cuprate $d$-wave superconductors implies a tight connection between the present QSL and superconductivity. This achievement demonstrates that the quantum-state representation using machine learning techniques, which had mostly been limited to benchmarks, is a promising tool for investigating grand challenges in quantum many-body physics.


I. INTRODUCTION
Collective excitations such as magnons and phonons consist of many elementary particles and provide us with fundamental understanding beyond the non-interacting picture, where the spontaneous symmetry breaking and associated Nambu-Goldstone bosons are required in many cases. Fractionalization, on the other hand, offers another route to realize emergent particles manifesting even in the absence of the symmetry breaking and serves as one of the central concepts in modern physics. The conventional elementary particles themselves can often be viewed as a bound state of more elementary objects, namely, the fractionalized particles, and such exotic particles emerge through the deconfinement. A prominent example of the deconfinement occurs in quantum chromodynamics: The proton and neutron that had been considered to be elementary particles before have turned out each to be a composite particle of three quarks with fractionalized charges, though quarks are hardly detected in experiments directly because of the confinement. In condensed matter, though the electron is an elementary particle in the vacuum, such deconfinement of electrons can be seen at low energies in specific circumstances fol- * yusuke.nomura@riken.jp lowed by the ground-state structure of materials. Consequential emergent fractionalized particles were discovered in examples of polyacetylene soliton [1] and fractional quantum Hall states [2]. The expectation would be that the emergent particles arising from the fractionalization still have particle character as low-energy excitations distinct from the elementary particles in the vacuum and the collective excitations in the symmetry broken states, and then would have novel functions in their many-body states, which may be useful for future applications such as quantum computing.
The QSL is a potential platform of such a fractionalization, where suppressed magnetic order by geometrical frustration of the spin interaction is expected to drive the fractionalization. The QSL phase was theoretically proposed both through numerical supports and mean-field theories [3,4]. Experimental efforts also supported the existence [4].
However, theoretical and experimental efforts have not yet identified and established the nature of fractionalized particles in reality due to their hidden nature and various theoretical difficulties. So far, several different types of QSL have been proposed. One of the important properties to characterize the QSL is the excitation spectra: They are classified, first, by whether the excitations are gapped [as in the cases of gapped Z 2 spin liquids (shortranged resonating valence bond (RVB) states) [5,6] and chiral spin liquid [7]], or gapless [8][9][10][11][12][13][14]. In the gapless case, one candidate is the gapless continuum of both of the singlet and triplet in an extended region of the Brillouin zone [10,11], which may arise, for example, if spin-1/2 fermionic spinons emerging from the fractionalization constitute a large Fermi surface (or line) as in U (1) spin liquid [12]. Another proposal is the spinon nodal liquid, where a small number of spinon gapless points appear in the Brillouin zone, resulting in the discrete gapless points of the observable spin excitation as well [13,14] (see Fig. 6 shown later for illustration). At the gapless points, the dispersion may be either linear (Dirac dispersion) or quadratic.
To establish the real existence of the QSL and then narrow down the nature of the QSL, we need to identify excitation spectra connected to experimental indications for a proper Hamiltonian that really accommodates the QSL state. However, it remains a challenge because of highly competing energies of various quantum states. We need a highly accurate framework for both ground and excited states in a momentum-resolved fashion.
Such high accuracy is offered by a recently developed machine learning method for the ground state. Here, we extend this method to represent both the ground and excited states. To be more precise, we employ the restricted Boltzmann machine (RBM) combined with pair-product (PP) states [15]. The RBM+PP method is further supplemented by two independent state-of-the-art numerical procedures, namely the correlation ratio [16] and level spectroscopy [17] methods, to reach the thermodynamic limit quickly by reducing the finite-size effect.
We then apply the RBM+PP to a candidate Hamiltonian of the spin-1/2 antiferromagnetic (AF) Heisenberg model on the square lattice with the nearest-neighbor and next-nearest-neighbor exchange interactions, J 1 and J 2 , respectively, called the J 1 -J 2 Heisenberg model. We employ two independent analyses to settle down the controversy and obtain firm evidence for the QSL phase: A finite range of the QSL phase in the region 0.49 J 2 /J 1 0.54 is found. In the QSL phase, the singlet and triplet excitations are both gapless at four symmetric momenta in support of the nodal Dirac (or quadratic touching) dispersion of the fermionic spinon at (±π/2, ±π/2) in the Brillouin zone, which brings the coexisting powerlaw decay of spin-spin and dimer-dimer correlations. The isomorphic structure of the gapless excitations of spinons at (±π/2, ±π/2) with the d-wave superconducting state in the cuprate superconductors is suggestive of a mutual profound connection.

II. J1-J2 HEISENBERG MODEL ON SQUARE LATTICE
The two-dimensional (2D) J 1 -J 2 Heisenberg Hamiltonian reads where S i is the spin-1/2 operator at site i, whose α (α = x, y, z) component is S α i = 1 2 c † i σ α c i with the electron operator c † i = (c † i↑ , c † i↓ ) and the Pauli matrix σ α . We set J 1 = 1 as the energy unit and we restrict the parameter range as 0 ≤ J 2 ≤ 1. i, j and i, j denote nearestneighbor and next-nearest-neighbor bonds, respectively.
In this model, the J 1 and J 2 interactions compete with each other (the former favors the Néel-type AF configurations, whereas the latter favors the stripe-type AF configurations). Although it is clear that Néel-and stripe-type AF phases exist for small and large J 2 regions, respectively, around J 2 = 0.5, which is the classical boundary between the Néel and stripe phases, unconventional quantum phase(s) such as QSL may emerge. Despite many theoretical efforts [18][19][20][21][22][23][24][25][26][27][28][29][30][31], the intermediate phase(s) of this model is still controversial. Indeed, the studies performing a systematic investigation of J 2 dependence with modern numerical techniques have proposed different scenarios: QSL (either gapless [23,31] or gapped [21]), valence bond solid (VBS) (either columnar [30] or plaquette [25]), or both of them [27,32]. Deconfined quantum criticality was also proposed instead of the QSL phase [25,28], which is interpreted as the QSL phase shrunk to a point and the fractionalization occurs only at this continuous phase transition point between the two symmetry-broken states.
To settle the phase diagram after various highly controversial proposals, one needs to satisfy at least the following three requirements: 1. Systematic investigation on finely-resolved J 2 dependence must be performed to establish whether the QSL exists as a phase in a finite J 2 interval because the QSL region is not expected to be wide.
2. Calculation must be highly accurate because quite different states are highly competitive with near degeneracy of the energies.
3. Reliable estimate of the thermodynamic limit to ensure it in the realistic bulk systems. When finitesize systems are studied, methods that are sizeindependent or have small finite-size effects are required to allow reliable extrapolation to the thermodynamic limit.
Although recent rapid progress in variational numerical methods has contributed to better accuracy, previous studies satisfying all the points hardly exist. Most of the studies have argued whether the order parameter is finite or not in the thermodynamic limit; however, the order parameter is tiny around the continuous phase transitions, and it is hard to discuss whether the order parameter is really zero or not. Also, when finite-size systems are studied, the direct extrapolation of the order parameter has a large finite-size effect. Exceptionally, Ref. [32] employed the level spectroscopy analysis, which can mitigate the finite-size effect. However, it is an indirect method, which speculates phase transitions in the ground state indirectly from the excitation structure. Because there exists no rigorous proof for the one-to-one correspondence of the ground state and excitation structure, one needs to verify in ground-state quantities to settle the highly controversial issue. As is detailed below, we employ the RBM+PP method, which offers a unique way of calculating ground state and momentum-space excitation dispersion in a systematically improvable and tractable way, to satisfy the conditions 1 and 2. The high accuracy and tractable computational cost of the RBM+PP enable comprehensive correlation ratio (ground-state property) and level spectroscopy (excited-state property) analyses with small finite-size effects. To fulfill the condition 3, a crosscheck from the two independent analyses is essential.

A. Machine learning for quantum many-body systems
Physical properties of many-body systems are governed by the eigenstates of the many-body Hamiltonian. Therefore, once the eigenstates of the Hamiltonian in Eq. (1) are known, we can predict the nature of the J 1 -J 2 model precisely. However, there is difficulty in obtaining eigenstates because the dimension of the eigenstates grows exponentially as the system size increases. In the present case where we consider the J 1 -J 2 Hamiltonian on the L × L (= N site ) lattice with the periodic boundary condition, we cannot obtain the exact wave function when N site 50. However, by using machine learning techniques, we can compress the data of eigenstates and approximate the wave functions accurately with a finite number of parameters.
Here, we employ a newly developed machine learning method, RBM+PP [15], to obtain accurate representations for both the ground and excited states. The RBM is a type of artificial neural network having two (visible and hidden) layers [33]. Using the machine learning technique, one can construct accurate many-body wave functions, which are systematically improvable toward the exact solution [34]. Indeed, it has been shown both theoretically and numerically that the RBM variational state flexibly describes a variety of quantum states [34][35][36][37][38][39][40][41][42][43][44][45], including the states exhibiting the volume-law entanglement entropy [35,37], which is advantageous to represent not only the ground state but also the excited states. Indeed, the RBM is shown to accurately describe excited states of quantum spin Hamiltonians [44,46], for which existing numerical methods often encounter numerical difficulties. Meanwhile, the PP state (called "geminal" in quantum chemistry) is represented by fermion wave functions, which can also accommodate volume-law entanglement. The PP state mapped onto bosonic spin space can represent RVB states [47], serving as a powerful starting point of the ground state approximation for the quantum spin systems [48]. The combined wave function, RBM+PP, inherits advantages of both and acquires much better accuracy than those achieved by either of the RBM or PP state separately [15]. By the RBM+PP method with quantum number projections (see below), we can calculate momentum resolved excitations.
The RBM+PP wave function Ψ(σ) = σ|Ψ with |σ = i c † iσ |0 is given by (we neglect normalization factor) [15] for each spin configuration σ = (σ 1 , σ 2 , . . . , σ Nsite ) with σ i = 2S z i = ±1. The number of sites is given by N site = L × L and the periodic boundary condition is assumed. The RBM part is given by (we omit irrelevant bias term on the physical spins) with the spin state of hidden units h k = ±1, the interaction between physical and hidden variables W ik , and the bias on the hidden variables b k . The number of hidden units is taken to be 16. The sum over hidden variables can be evaluated analytically and Eq. (3) can be efficiently computed as φ RBM (σ) = k 2 cosh b k + i W ik σ i . To make it possible to express the sign change of the wave function, we take the b k and W ik variational parameters to be complex. The PP state mapped onto spin systems reads with real variational parameters f ↑↓ ij . ψ PP (σ) in Eq. (2) is related as ψ PP (σ) ≡ σ|ψ PP . Here, P G = i (1 − n i↑ n i↓ ) with n i↑ = c † i↑ c i↑ and n i↓ = c † i↓ c i↓ is the Gutzwiller projection prohibiting double occupancy.
We optimize the variational parameters {b k , W ik , f ↑↓ ij } to minimize the energy E = Ψ|H|Ψ Ψ|Ψ . The energy is a highly nonlinear function with respect to the parameters {b k , W ik , f ↑↓ ij }. Therefore, by interpreting the energy as a loss function, the task of obtaining the lowest-energy state can be recast as a machine-learning task, namely a high-dimensional optimization problem of the highly nonlinear function (RBM+PP) using the highly nonlinear loss function (energy) [49]. The details of the optimization method and the calculation conditions can be found in Appendix A.

B. Strategy to overcome numerical challenges
Various competing controversial scenarios have been proposed for the phase diagram and the nature of possible QSL as we mentioned above. The machine learning only is, despite its crucial importance, not enough to resolve these controversies. In fact, even when we obtain accurate representations of quantum states by the machine learning, (i) another challenge is how to reach quick convergence to the thermodynamic limit from available finite-size results (condition 3 listed in Sec. II). Furthermore, provided that the QSL phase exists, the next challenge is to elucidate its nature; (ii) it is essential to estimate the excitation gap structure and momentum resolved dispersion accurately. To overcome the challenge (i), the present paper employs an unprecedented combination of two methods and one supplementary analysis together and reaches quantitative agreements, which ensures the accuracy because the two methods are originally independent of each other. As a computational method to identify the quantum phases, this is the first attempt to use such combinations, and it successfully establishes a way to obtain the accurate phase diagram, which may serve as the standard method in the future. The first method is the correlation ratio method [16], which utilizes the ground state properties (see Sec. III B 1). The second is the level spectroscopy [17], which detects the signature of the phase transition in the excitation spectra (see Sec. III B 2). Both methods show small finitesize effects and quickly converge to the thermodynamic limit. These methods were developed independently and, in fact, measure the excited and ground-state properties, respectively, which are originally independent. However, the important point is that they have the one-to-one correspondence, conceptually similar to the fluctuationdissipation theorem and Kubo formula, between the equilibrium and non-equilibrium excited states. Such correspondence and match in the calculated results help to ensure the reliability of the phase diagram. Further, the obtained phase boundary is supported by the standard finite-size scaling method thanks to the universal scaling relations (see Sec. III C). For (ii), we use quantum number projection to reach the accuracy on the spectroscopy level [50] (see Sec. III B 3). Here, we address the advantages of employing these methods.

Correlation ratio
Correlation ratio R quantifies how sharp the structure factor peak is. R is given by R = 1 − S(Q + δq)/S(Q) [16,51], where S(q) is the structure factor, Q is the peak momentum, and Q+δq is the neighboring momentum. In the case of the square lattice, δq = (2π/L, 0) and (0, 2π/L). We see that the R value approaches 1 (0) when the peak becomes sharp (broad). Therefore, with increasing system size, R scales to 1 in the ordered phase with delta-function Bragg peak and 0 in the disordered phase. The crossing point of R curves for different system sizes does not depend sensitively on system size. Thus it is suitable for an accurate estimate of the phase boundary between ordered and disordered phases in the thermodynamic limit [16,51].
We examine R for both spin-spin and dimer-dimer cor-relations to detect Néel-AF and VBS transition points, respectively. The spin-spin correlation is given by C s (r i − r j ) = S i · S j . The dimer-dimer correlation is de- with the dimer operator D α i = S i · S i+α on the nearest-neighbor bonds for the α-direction (α = x, y). Hereafter, the subscripts "s", "d x ", and "d y " are used for spin-spin, dimer-dimer (α = x and α = y) correlations, respectively. Then, the structure factor is calculated from S γ (q) = 1 Nsite i,j C γ (r i − r j )e iq·(ri−rj ) with γ = s, d x , and d y . The two correlation ratios R Néel and R VBS are defined from S s (q) and S dx (q) [or equivalently S dy (q)] to determine the Néel-AF and VBS transition points, respectively. Close to the Néel-AF phase, the peak momentum is Q = (π, π) for S s (q). For VBS, Q = (π, 0) for S dx (q) and Q = (0, π) for S dy (q).

Level spectroscopy
Quantum phases are characterized by their unique structure in excitation spectra. At finite sizes, if the phases are different, low-lying excitations will be characterized by different quantum numbers. Therefore, the transition point can be estimated by the size extrapolation of the crossing point of the low-lying excitation energies [17]. This level spectroscopy method is known to have small system size dependence as well. Indeed, it has played an important role in precisely determining the Berezinskii-Kosterlitz-Thouless transition point for the sine-Gordon model [17]. This method offers an analysis completely different but complementary to the correlation ratio method.

Quantum number projection
The eigenstates of the Hamiltonian in finite size systems are labeled by quantum numbers. By optimizing the RBM+PP wave function for each quantum number sector, we can obtain both the ground state and lowlying excited states. We apply total-momentum and spinparity projections to the RBM+PP wave functions to specify the quantum number [50]: (double sign in the same order). S + (S − ) indicates even (odd) spin parity corresponding to even (odd) values of the total spin S. K is the total momentum. T R is a translation operator shifting all the spins by R. For each quantum number sector, we optimize the RBM+PP wave function to obtain the lowest-energy state. Although the spin-parity projection can only distinguish whether S is even or odd, we always obtain a singlet state for the even S sector and a triplet state for the odd S sector (we confirm it by calculating S expectation value for the obtained states). This is because the singlet (triplet) state is the lowest-energy state for each even (odd) S sector. We note that the quantum number projection is helpful not only to distinguish quantum numbers but also to lower the variational energy [52]. The ground state is given for K = (0, 0) and even S sector. The energies for other quantum number sectors measured from the ground state energy determine the excitation spectra. Then, we can obtain singlet and triplet excitations separately with momentum resolution. Exceptionally, we need special treatments to obtain S = 0 excited state at K = (0, 0) and S = 2 excited states, which are described in detail in Appendix A. As we mentioned above, the flexible representability of the RBM+PP gives accurate representations not only for ground states but also for excited states. The accurate estimate of momentum-resolved excitation gaps enables us to perform the above-mentioned level spectroscopy and also to elucidate the nature of the QSL phase.

C. Finite-size scaling method
Near the quantum critical point, the susceptibility χ at the ordering wave vector Q in finite-sized systems follow the following finite-size scaling form [53] where the universal scaling function f χ appears with the correlation length exponent ν and the susceptibility exponent γ. Here, t assumed to satisfy t ≪ 1 is the dimensionless distance to the critical point. In the present case Through the relation between χ and the structure factor S(t, Q, L) given by χ(t, Q, L) ∼ S(t, Q, L)L z with the dynamical exponent z, we find that the squared order parameter m 2 = S(t, Q, L)/L d for the d-dimensional system follows if the Fisher's scaling relation γ/ν = 2 − η holds for η associated with the anomalous dimension characterized by the power-law decay of the correlation, C(r) ∼ 1/r d+z−2+η for distance r = |r| at the critical point. Then the finite-size scaling plot should exhibit the universal scaling function f χ .

IV. RESULTS
First, we check the accuracy of the RBM+PP method in analyzing the J 1 -J 2 Heisenberg model (see Appendix B). We have confirmed that the RBM+PP achieves state-of-the-art accuracy not only among machine-learning-based methods but also among all available numerical methods. Indeed, the RBM+PP wave function marks the best precision for the ground state calculations among the compared methods for the  8 × 8 and 10 × 10 lattices (see Fig. 10 and Table II). We have also found that the RBM+PP represents excited states with unprecedented accuracy (Fig. 11).
Also, in our RBM+PP method, the computationally most demanding part is coming from the PP part, and the neural-network (RBM) part offers an efficient way of improving accuracy without increasing the scaling of computational cost, i.e., as compared to the PP only calculations with the computational cost of O(N 3 site ), the computational cost increases only by O(1). This is in contrast to the Lanczos method, which is also used to improve the variational energy (see Appendix B): if we apply the p-th order Lanczos step to improve the PP only calculations, it increases the computational cost by O(N p site ), and hence the total computational cost of the Lanczos-applied PP calculations scales as O(N p+3 site ). Although we calculated the ground state and various excited states independently, a tractable computationalcost scaling of the RBM+PP method allowed us to perform numerous independent calculations for large system sizes within given computational resources. Thus obtained high-quality data contribute to a reliable determination of the phase diagram consistently from both ground-state and excitation analyses (see below).

A. Ground-state phase diagram
The RBM+PP method combined with the state-ofthe-art numerical techniques convincingly uncovers the phase diagram of the J 1 -J 2 Heisenberg model as shown in Fig. 1. In the small (large) J 2 region, the Néel-type (stripe) AF long-range order appears as in the classical phase diagram. In between these two phases, nonmagnetic ground states, QSL and VBS, are found in the region J Néel ≈ 0.61, respectively. Whereas VBS breaks lattice symmetry, QSL does not break any. Clearly and notably, QSL is stabilized in a finite region of J 2 around J 2 = 0.5. The phase transition between VBS and stripe-AF at J V-S 2 is of 1st order, which is characterized by the kink in the ground state energy, while the other two transitions are continuous (Fig. 13 in Ap- pendix C). Below, we describe the procedure to determine the continuous phase transition points.

Phase boundary determined by correlation ratio
Results for the correlation ratios, R Néel and R VBS , are shown in Figs. 2(a) and 2(b), respectively (see Figs. 14 and 15 in Appendix C for the raw data of correlation functions). We see clear crossings of curves for three sizes at nearly the same points at J 2 = J Néel 2 ≈ 0.49 for R Néel and at J 2 = J VBS 2 ≈ 0.54 for R VBS . This standard procedure strongly supports that the two transitions associated with the Néel-AF and VBS ordering take place at the different points close to these system-size independent crossings. It supports the existence of an intermediate phase without any long-range ordering, i.e., QSL phase in the range 0.49 J 2 0.54 (see Appendix C for the discussion of the system-size dependence of the crossing points).

Phase boundary determined by level spectroscopy
The level spectroscopy method was applied to the 2D J 1 -J 2 Heisenberg model before [32]. They interpreted the crossing between the lowest singlet and triplet excitations as the VBS-order boundary, following Ref. 54. In addition, they found the singlet-quintuplet crossing and interpreted it as a signal of the disappearance of the AF long-range order, because the transition from the AF long-range order to quasi-long-range order in onedimensional Heisenberg model with long-range interaction shows a similar behavior [32,55]. These two crossings extrapolated to L → ∞ limit gave different J 2 val-ues: J 2 = 0.463(2) and J 2 = 0.519(2) for the singletquintuplet and singlet-triplet crossings, respectively.
To critically crosscheck the consistency with the above correlation ratio result, we also reexamine the level spectroscopy analysis as a complementary check. We here enjoy the advantage of the momentum resolution in addition (contrary to Ref. 32). Figure 3 shows J 2 dependence of the excitation energies ∆ for sizes (a) 12 × 12 and (b) 16 × 16 at high-symmetry momenta. The singletquintuplet and singlet-triplet crossings signaling the AF-QSL and QSL-VBS transitions, respectively, are highlighted by arrows. The size extrapolation of the crossing points is shown in Fig. 4(a). We use L −2 scaling as in Refs. 32 and 54. The extrapolated thermodynamic values are J 2 = 0.493(2) and J 2 = 0.532(2) for the singletquintuplet and singlet-triplet crossings, respectively. The values are close to those of Ref. 32 above. Quantitative differences may well be ascribed to the smaller system sizes calculated in Ref. 32 than ours. As for the singlettriplet crossing, our result is also consistent with a more recent estimate by the variational Monte Carlo (VMC) method, which gives J 2 = 0.542(2) [56].
More importantly, our phase boundary estimated by the level spectroscopy has a striking quantitative agreement with the correlation ratio result described above. It is of great significance to see the one-to-one correspondence between the ground-state phases and the excitation structures. We then safely conclude that a finite QSL region around J 2 = 0.5 emerges (see Supplementary Note 1 in Appendix D for additional noteworthy features found in the level spectroscopy). Figure 4(b) further shows the size dependence of the excitation gap ∆ at the crossing points. ∆ × L seems to converge at a finite value as L → ∞ for both crossings. Therefore, the two critical points corresponding to AF-QSL and QSL-VBS transitions become gapless in the thermodynamic limit with the scaling ∆ ∝ 1/L.

B. Excitation spectrum in quantum spin liquid phase
As we see in Fig. 4(b), the singlet excitation with K = (π, 0) and (0, π) becomes gapless at both AF-QSL and QSL-VBS critical points, implying that it is gapless through the QSL region sandwiched by these two critical points. In the QSL phase, the triplet excitation at K = (π, π) is the lowest excited state in finite-size systems [lower than the gapless singlet at (π, 0)] lending support to the vanishing gap also for (π, π) triplet in the thermodynamic limit. By the excitation involving the triplet at (π, π) and the singlet at (π, 0), one can construct the triplet (0, π), which must be gapless if these two elementary excitations are excited far apart in the thermodynamic limit, even when they are repulsively interacting. In a similar way, one can construct a gapless singlet excitation at (π, π) and (0, 0). Therefore, the singlet and triplet excitations are both gapless at (0, 0), (a) System-size dependence of singlet-quintuplet (red dots) and singlet-triplet (black squares) level crossings indicated by red and black arrows in Fig. 3. The extrapolation to the thermodynamic limit is done by the polynomial fit a + b/L 2 + c/L 4 (solid curves). (b) System-size dependence of the excitation gap ∆ at the two level crossings. For the singlet-quintuplet level crossing in (a), the 18 × 18 data are added to corroborate the result.
(Supplementary Note 2 in Appendix D). On the other hand, the gap stays nearly constant at the intermediate K points. By combining the gap analysis at the critical points (see above) and the size extrapolation of the gap by the scaling a + b/L at the intermediate K points, we draw dispersion expected in the thermodynamic limit. The excitation spectra in the thermodynamic limit exhibit unconventional behavior in which the gap vanishes at the four high-symmetry momenta. We find only these four points as the gapless excitations suggesting Diractype linear dispersion around these four points. To cor- (π, 0) (π, π) (0, π) 0 q x q y q x q y (π/2, π/2) (π/2, −π/2) (−π/2, −π/2) (−π/2, π/2) 6. (a) Weight of lowest branch in the dynamic spin structure factor for q = (π, 0) and (π, π) for J2/J1 = 0.5. At each q point, the weight is normalized by the total spectral weight ∞ 0 dωSs(q, ω). (b) Schematic picture for plausible spinon dispersion around gapless points (±π/2, ±π/2), illustrated both for particle (pink) and hole (green) sides above and below the spinon Fermi energy. Two examples of two spinon excitations (two red and two black circles) are illustrated (see below). (c) The observable spin excitation is constructed from the two spinon excitations, which generates the gapless points at (π, 0), (0, π), (0, 0) and (π, π). For instance, the red circle with the momentum around (−π, 0) is constructed from the two spinon excitations shown as the small red circles with the momenta around (−π/2, π/2) and (−π/2, −π/2) in (b). The black circle is another example of spin excitation originated from the two spinon excitations shown as the small black circles in (b). Continuum incoherent spin excitations inside the cones are generated from the combinations of the two spinon excitations on the pink or green cone surfaces in (b).
roborate the conclusion about the four Dirac-type gapless points in the QSL phase, we have also calculated the excitation energies at (mπ/3, nπ/3) with m, n = 0, 1, 2, 3 for 12 × 12 lattice (Fig. 17 in Appendix C). From the limited momenta we studied, although other possibilities such as the higher-order dispersion (e.g., quadratic band touching) or tiny but extended gapless regions rather than points are not excluded, the results in Fig. 17 also support the Dirac-type nodal QSL.

C. Signature of fractionalization in quantum spin liquid
In the present QSL phase, one can expect an exotic fractionalization of particles, where a charge-neutral spin-1/2 excitation, called spinon, emerges. Although the spinon excitation cannot be detected experimentally, the evidence of the fractionalization can be detected as an incoherent continuum in the dynamic spin structure factor S s (q, ω) (spin-1 excitation) [57], which is interpreted by the two-particle (two-hole) or particle-hole excitation continuum of the spinons. We here compute the weight in S s (q, ω) at q = (π, 0) and (π, π) for the lowest triplet excitation shown in Fig. 5. If the excitation were the conventional magnon branch of a magnetic phase, the weight would be the order 1. If the weight vanishes, most of the weight lies in incoherent continuum at higher energies, supporting the emergence of fractionalized spinons [57]. Figure 6(a) shows the weight of the lowest branch in S s (q, ω) for q = (π, 0) and (π, π). We indeed see that the weight decreases as the system size increases. In particular, the weight at q = (π, 0) rapidly decreases to zero, which means that the spectral weight is dominated by the incoherent continuum. [We do not analyze the behavior at q = (π, π) in detail because of a numerical challenge due to the proximity to AF(Néel)-QSL phase boundary J 2 = J Néel 2 ≈ 0.49 (Supplementary Note 3 in Appendix D)]. This is a strong evidence that the fractionalization indeed occurs in the QSL phase of the J 1 -J 2 Heisenberg model. As we will discuss in Sec. V, the dispersion of the emergent fractionalized spinon is expected to be gapless at the points (±π/2, ±π/2) [ Fig. 6(b)].
D. Real-space correlation functions in the quantum spin liquid phase Figure. 7 shows the real-space decay of spin-spin and dimer-dimer correlation functions, |C s (r)| and |C dx (r)|, respectively, for the diagonal direction (r x = r y ) for 16 × 16 lattice in the QSL phase (J 2 = 0.5125) (for the definition of the correlation function, see Methods). If the correlation function shows power-law decay, it is expressed by the exponent z + η, namely the spin-spin and dimer-dimer correlations should show C(r) ∼ r −(z+η) (r = |r|) in the real space as critical behavior. Both correlation functions indeed show consistent behaviors with the power-law decay. It evidences the dual critical nature of the VBS and Néel-AF correlations in the QSL phase, and this ground-state property is consistent with the gapless singlet and triplet excitations clarified independently (Sec. IV B). On top of the one-to-one correspondence between the ground-state phases and excitation structures revealed by the correlation ratio and level spectroscopy analyses (Secs. IV A 1 and IV A 2), we again demonstrate (solid) and 1.62 (dashed), in which we consider the effect of the periodic boundary condition [27]. The values of z + η are taken from the analysis in Fig. 9(c). The upturn at large |r| is due to the periodicity of the lattice. a nice correspondence between the ground-state and excitation properties.

FIG. 8. Finite-size scaling analysis. (a,c) Data collapse for
Nèel-AF order parameter. We assume J Néel 2 /J1 = 0.49 and estimate the critical exponents z+η and ν. The Bayesian scaling analysis [58,59] gives z + η = 1.410(4) and ν = 1.21 (5). (b,d) Data collapse for the VBS order parameter. The same analysis with assuming J VBS 2 /J1 = 0.54 gives z +η = 1.436 (6) and ν = 0.67 (2). Solid curves are the inferred scaling functions. In (a) and (c), the 18×18 data are added to corroborate the result. The figure shows that the conventional finite-size scaling analysis consistently supports the results obtained by the correlation ratio and level spectroscopy. spectively.

V. DISCUSSION
As is discussed in Sec. II, to settle the highly controversial situation on the phase diagram, the calculations need to fulfill three conditions: 1. systematic J 2 dependence survey, 2. high accuracy, and 3. reliable estimate of the thermodynamic limit. The RBM+PP data achieves the state-of-the-art accuracy level both for ground-state and excited states, satisfying the condition 2. With the high accuracy, we have performed a systematic investigation on the J 2 dependence both for ground-state and excited-state properties (condition 1). For the condition 3, both the correlation ratio and level spectroscopy have given consistent results, supporting the conclusion of the QSL phase in the region 0.49 J 2 /J 1 0.54. We do not find such quantitative consistency before, and we became convinced of the existence of the QSL phase only after finding their consistency. Nevertheless, we note that, at a qualitative level, an overall consensus on the existence of the QSL is being formed among the best accurate methods (Refs. 23 and 25 and ours) clarified in the benchmark shown in Appendix B (Note that Ref. [25] obtained essentially vanishing order consistent with our finite QSL region, though they considered alternative possibilities as well, which was not settled within their analyses of the size dependence of the order parameter correlation). The spin excitation dispersion has been rarely studied in the literature except for the studies obtained by assuming a priori a variational form of Z 2 nodal spin-liquid wave function [60,61]. In Ref. 61, the spin cluster perturbation method is also employed to draw the dispersion. Our gapless structure lends support to these variational and the spin cluster perturbation studies in qualitative features, though our results have been obtained without such assumptions and approximations. Together with the consideration on the stability of the QSL phase [62] and the reason discussed below, our unbiased analysis evidences the QSL phase in the J 1 -J 2 Heisenberg model characterized by Z 2 nodal QSL (rather than U (1) QSL) with gapless and fractionalized spin-1/2 spinon excitations at (±π/2, ±π/2), proposed in an earlier study [23] (we did not exclude the possibility of U (1) QSL just from the spin excitation spectra because the Z 2 and U (1) QSL give very similar S s (q, ω) [62]).
The real spin excitations measurable in experiments must be made of two-spinon excitations, and thus the singlet and triplet gapless points are (0, 0), (π, 0), (0, π) and (π, π) [ Fig. 6(c)]. The gapless Dirac-type excitations in both singlet and triplet sectors show an excellent consistency with the dual critical nature of the VBS and Néel-AF correlations, decaying algebraically in the real space, in the QSL phase (Sec. IV D).
The finite-size scaling analysis shown above suggests that the value for critical exponent z + η is about 1.4 for both of the AF-QSL and QSL-VBS critical points (Fig. 8), which is suggestive of an emergent symmetry between the spin-spin and dimer-dimer correlations, associated with the Néel-AF and VBS orders, respectively. If the U (1) QSL is realized as a phase, we will see the emergent symmetry within the whole QSL region as the critical phase [63,64]. However, the power-law exponent z + η seems to change in the QSL region: While it increases as J 2 increases for the spin-spin correlation, the dimer-dimer correlation shows the opposite trend [ Fig. 9(c)]. It supports that the QSL with the emergent U (1) symmetry is absent for an extended J 2 region and implies the extended region of the Z 2 QSL instead. . It will be of great interest to investigate this issue further in the future, especially by considering more detailed system size dependence to further establish the thermodynamic behavior.
Since the excitation structure is isomorphic with the charge and spin excitations of the d-wave superconducting state in the cuprate superconductors, it is suggestive of the connection of the two; the superconducting state could be borne out from the present QSL immediately when carriers are doped. The present accurate estimate of the spinon excitation, especially, incoherent nature of the spin excitations with continuum, will provide us with insights into the unsolved puzzles of the cuprate superconductors including the incoherent transport and charge dynamics.

VI. SUMMARY
We have studied the 2D J 1 -J 2 Heisenberg model using a highly accurate machine-learning method, RBM+PP. Our achievements are summarized into the following points: the quantitative estimate of the phase diagram, useful insights into the QSL property to understand its nature, and the establishment of one-to-one correspondence between ground-state and excitation structure.
First, by combining the RBM+PP with the correlation ratio and level spectroscopy methods, we have been able to extrapolate to the thermodynamic limit reliably by two independent analyses. The agreement reached between the two at an unprecedented level has given the firm evidence for a finite QSL region 0.49 J 2 /J 1 0.54. The phase diagram is summarized in Fig. 1.
The QSL is characterized as the dual nature of the algebraic and coexisting correlations of the antiferromagnetic (associated with the Néel order) and dimer (associated with the VBS order) correlations, which had been thought incompatible before by the symmetry difference. The elucidated dual nature is also seen consistently in the excitation property: We have identified the Dirac-type dispersion with gapless points (0, 0), (π, 0), (0, π), and (π, π) in both the singlet and triplet excitation sectors (related each to the dimer-dimer and spin-spin correlations, respectively). The excitation structure is consistent with the emergence of the fractionalized spin-1/2 spinons with gapless Dirac dispersion. Interestingly, the power-law decay exponents of these two correlations change as a function of J 2 /J 1 and do not coincide except for a single point around J 2 = 0.52, which imposes a substantial constraint on the gauge structure of the QSL.
Finally, our comprehensive calculations have revealed a fundamental "law of correspondence" between the ground-state and excitation structure in the J 1 -J 2 Heisenberg model. By establishing the phase diagram, we have demonstrated that the evolution of the ground state indeed maps to the change in the excitation structure induced by the level crossing in a fingerprint fashion with one-to-one correspondence. We have also shown that the coexisting power-law decay of the dimer-dimer and spin-spin correlation functions in the real space in the QSL phase (ground-state property) consistently corroborates the gapless structure of singlet and triplet excitations, respectively. Such one-to-one correspondence has a fundamental significance in physics, as the oneto-one correspondence between the equilibrium and nonequilibrium excited states addressed in the fluctuationdissipation theorem and Kubo formula gives a foundation for the understanding of the linear response.
Such accurate, systematic, and comprehensive elucidation of the QSL with insights into the duality of the gapless correlations and the law of correspondence has been enabled by the RBM wave function combined with the PP state and the quantum number projection that offers state-of-the-art accuracy within a tractable computational cost: The high accuracy and the tractable computational-cost scaling of the RBM+PP method [O(N 3 site )] were necessary to prepare comprehensive highquality data for large system sizes to accomplish our achievement.
So far, the machine learning methods had been applied mostly to benchmark problems with known solutions. By combining the RBM+PP wave function with cutting-edge methods to reduce finite-size corrections, we have succeeded in uncovering QSL in the long-standing challenging problem. This achievement opens a new avenue of numerical methods applicable to the grand challenges of quantum many-body systems.

ACKNOWLEDGMENTS
We acknowledge useful discussions with Satoshi Morita, Anders W. Sandvik, and Zi Yang Meng. We also thank Satoshi Morita for providing us with the raw data in Ref. 27. Y.N. is grateful for fruitful discussions with Ribhu Kaul, Hidemaro Suwa, Yoshitomo Kamiya, Kenji Harada, Zheng-Cheng Gu, Giuseppe Carleo, and Ryui Kaneko. The implementation of the RBM+PP scheme is done based on the mVMC package [65]. The computation was mainly done at Supercomputer Center, Institute for Solid State Physics, University of Tokyo, and RIKEN supercomputers K and Fugaku. The authors are grateful for the financial support by a Grant-in-Aid for Scientific Research (Grant No. 16H06345) from Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan. Y.N. was supported by Grant-in-Aids for Scientific Research (JSPS KAKENHI) (Grants No. 17K14336 and No. 18H01158). This work is financially supported by the MEXT HPCI Strategic Programs, and the Creation of New Functional Devices and High-Performance Materials to Support Next Generation Industries (CDMSI) as well as by "Program for Promoting Researches on the Supercomputer Fugaku" (Basic Science for Emergence and Functionality in Quantum Matter -Innovative Strongly-Correlated Electron Science by Integration of "Fugaku" and Frontier Experiments -) (Grants No. hp200132 and hp210163). We also acknowledge the support provided by the RIKEN Advanced Institute for Computational Science under the HPCI System Research project (Grants No. hp170263, hp180170, and hp190145).
where the local energy E loc (σ) is given by E loc (σ) = σ ′ σ|H|σ ′ σ ′ |Ψ σ|Ψ . The E value depends on the variational parameters. To optimize the variational parameters to minimize E, we employ the stochastic reconfiguration (SR) method [66], which is equivalent to the imaginary-time Hamiltonian evolution e −τ H |Ψ within the Hilbert space spanned by the RBM+PP wave function. Because the imaginary-time Hamiltonian evolution e −τ H |Ψ always stably gives the lowest-energy state for each quantum number sector (as far as the initial RBM+PP state is not orthogonal to the lowest-energy state), the SR method enables stable optimizations. For further technical details of the SR optimization, we refer to Ref. 15.
The number of complex variational parameters in the RBM part is N hidden for b k and N hidden × N site for W ik , respectively. As for the real variational parameters f ↑↓ ij in the PP part, to reduce the number of parameters and the computational cost, we impose 4 × 4 sublattice structure for the 8 × 8, 12 × 12, and 16 × 16 lattices, and 6 × 6 sublattice structure for the 18 × 18 lattice, whereas we do not employ sublattice structure for the 6 × 6 lattice. In the case of the 4 × 4 sublattice structure, the number of independent f ↑↓ ij parameters is reduced from N 2 site to 4 × 4 × N site = 16N site , and the other f ↑↓ ij parameters are defined by spatial translation operations. In the presence of the Gutzwiller factor to map the PP state onto the spin system, the onsite f ↑↓ ii parameters become completely redundant, i.e., the wave function does not depend on f ↑↓ ii at all. Then, the number of relevant f ↑↓ ij parameters is 16(N site −1). For the initial values for {b k , W ik , f ↑↓ ij }, we put random numbers in order not to introduce bias in the initial variational state. More specifically, for each real and complex part of b k and W ik parameters, we put small random numbers from the interval [−0.05, 0.05]. In the case of the triplet state calculation, b k parameters are multiplied by 10 (note that if b k is zero, the RBM part is completely symmetric with respect to the global spin inversion). For the initial f ↑↓ ij parameters, we put random numbers from [−f max , f max ] with f max depending on the distance R ij between ith and jth sites. We typically take f max to be proportional to R −a ij with a ∼ 2. Random f ↑↓ ij parameters allow various spin ordering patterns whose period is within the sublattice size. A longer period structure than the system size is beyond the scope of this study, as are all the other earlier finite-systemsize studies. For each J 2 point, we perform at least three independent optimizations of the RBM+PP wave functions from different initial variational parameters. We discuss the initial-parameter dependence in more detail in Appendix B.
The computational cost of the RBM+PP wave function employing sublattice structure in the PP part scales with O(N 3 site ). In the RBM+PP method, a computationally demanding part is coming from the calculation of the PP wave function part and the neural network (RBM) part offers an efficient way of improving accuracy.
Special treatments to obtain some specific excited states As we describe in Sec. III B 3, we apply the spin-parity projection to distinguish whether the total spin S is even or odd. Because the singlet (triplet) state is the lowest state for each even (odd) S sector in the present study, we obtain a singlet (triplet) state for the even (odd) S sector. Therefore, we can obtain the singlet (S = 0) and triplet (S = 1) excited states with momentum resolution. However, we need special treatment to obtain S = 0 excited state at K = (0, 0) because the lowest-energy state in S = 0 and K = (0, 0) quantum number sector is the ground state. We use additional simplified point-group projection on top of those in Eq. (5) to obtain excited states belonging to a different irreducible representation of the C 4v point group of the square lattice than that of the ground state as follows: where the R π/2 is an operator to rotate the spin configuration by 90 degrees. With this projection, we can distinguish whether the state belongs to A (either A 1 or A 2 ) irreducible representation or B (either B 1 or B 2 ) irreducible representation under the C 4v point group (to distinguish between A 1 and A 2 or between B 1 and B 2 , we need full point group projection with 0, π/2, π, 3π/2 rotations). The ground state corresponds to the former, while the excited state corresponds to the latter.
We also need special treatment to obtain S = 2 excited states. To this end, we use the mVMC (many-variable variational Monte Carlo method) [65] based on the PP wave function. In the mVMC, the full spin projection to specify the total spin is available, and we apply it to get S = 2 states. The full spin projection is time-consuming (at least about five times) compared to the spin-parity projection. At the cost of longer computational time for the full spin projection, the mVMC (only PP) gives comparable accuracy to the RBM+PP method.

Calculation conditions
In the present study, we fix the number of hidden units N hidden to be 16. We always apply the spin-parity and momentum projections during the optimization of the RBM+PP wave function. The special treatments to obtain S = 0 excited state at K = (0, 0) and S = 2 excited Comparison of the ground state energy for the J1-J2 Heisenberg model. The comparison is made among the variational energies under the periodic boundary condition. The system sizes are (a) 6 × 6 and (b) 8 × 8. Our RBM+PP results are compared with those obtained by the variational Monte Carlo (VMC) method combined with the p-th order Lanczos steps [23], the density-matrix renormalization group (DMRG) (with 8182 SU (2) states) [25], the convolutional neural network (CNN) [67], and the exact diagonalization (ED) [68]. The CNN and ED results are available only for the 6 × 6 lattice.
states are described above. To improve the quality of the data for the correlation function in Figs. 2,6,7,8,9,14,15, and 16 quantitatively, we apply the simplified pointgroup projection in Eq. (A2) to the optimized ground state RBM+PP wave function for the sector with S = 0 and K = (0, 0). The ground state energy in Fig. 10 is also produced with the simplified point-group projection (see Appendix B for the details).

Appendix B: Benchmark
Accuracy of the RBM+PP wave function By applying the RBM+PP method to the 2D J 1 -J 2 Heisenberg model on the square lattice, we confirm that the RBM+PP achieves state-of-the-art accuracy not only among machine-learning-based methods [67,[69][70][71] but also among all available numerical methods. Figure 10 shows the comparison of the ground-state energy among various methods for the 6 × 6 and 8 × 8 lattices (see Table I for the raw data). Here, the RBM+PP energy is obtained by optimizing the RBM+PP wave function with the momentum, spin-parity, and simplified-point-group projections. We do not employ the sublattice structure in the f ↑↓ ij parameters (the sublattice structure used in the actual calculations is discussed in Appendix A), and the number of hidden units is 16 as commonly employed in the paper. Up to the 6 × 6 lattice, the exact diagonalization result is available. At J 2 = 0.5, where the frustration is strong, the relative error of the RBM+PP energy is less than 0.01 %, demonstrating the high accuracy of the RBM+PP wave function. For the 8 × 8 lattice, the RBM+PP wave function gives the best accu-   (2) rate energy among the compared variational methods for all J 2 values we studied.
We have also performed the benchmark calculations for the 10 × 10 lattice because the benchmarks of neuralnetwork wave functions in the literature have mainly been performed using the 10 × 10 lattice. We optimized the RBM+PP wave function with 16 hidden units (as used in the other system sizes) without introducing a sublattice structure in the PP part. We apply the momentum, spinparity, and point-group projections. Table II shows the comparison of the ground-state energy at J 2 = 0.5 among different wave functions. As in the 8 × 8 lattice result, the RBM+PP gives the best accuracy among the various methods. From the systematic benchmarks on the 6 × 6, 8×8, and 10×10 lattices, we conclude that the RBM+PP achieves state-of-the-art accuracy.
In the actual calculations, we employ the sublattice structure in the f ↑↓ ij parameters for the 8 × 8, 12 × 12, 16 × 16, and 18 × 18 lattices, to reduce the computational cost from O(N 4 site ) to O(N 3 site ) (see Appendix A). The reduction of the computational time enables us to perform systematic calculations for various J 2 values and for different quantum-number sectors. By employing the sublattice structure, the accuracy becomes slightly worse compared to that without sublattice structure. For example, in the case of the 8 × 8 lattice at J 2 = 0.5, the ground-state energy with the 4 × 4 sublattice structure is −0.498460 (6), which is compared to −0.498886(1) obtained without a sublattice structure. The difference is less than 0.1 %; therefore, high accuracy is retained even with the sublattice structure [73]. As for the spin-spin and dimer-dimer correlations, the obtained values of the order parameters are m 2 Néel = 0.06955 (8) and m 2 VBS = 0.01720 (3) in the case of the 4 × 4 sublattice structure, and m 2 Néel = 0.06724 (8) and m 2 VBS = 0.01703 (3) in the case of no sublattice structure. The actual calculations with sublattice structure tend to slightly overestimate the order parameters in the frustrated regime. We see a similar tendency in the case of the benchmark calculations of RBM only wave functions for the 6 × 6 lattice at J 2 = 0.5 with changing the number of hidden units [74], where the Néel-AF order parameter tends to be overestimated for a small number of hidden units. With increasing the number of hidden units, the accuracy improves, and the order parameter shows an excellent agreement with the exact results [74]. Considering the fact that improving accuracy tends to suppress the order parameter, our statement of the existence of the QSL phase with vanishing order parameters in the thermodynamic limit should be valid.
Remarkably, we also find that the RBM+PP accurately represents excited states as well as the ground state. Figure 11 shows the comparison of excitation energies for singlet, triplet, and quintuplet excitations between the exact and RBM+PP results for the 6×6 lattice. The agreement is excellent, where the difference in energy between the exact and the RBM+PP results is less than 0.01. Previously, there have been several attempts to obtain the excitation gap of the J 1 -J 2 model [21,23,25,32]. In Ref. 23 using the combination of the VMC and Lanczos methods, the excited states are obtained by changing boundary condition, which limits the number of excited states that can be calculated [only S = 2 with the momentum (0, 0) and S = 0 with (π, 0) or (0, π)]. Also, the accuracy does not reach the level shown in Fig. 11 even with the 2nd-order Lanczos being applied [VMC(p = 2)]. In Refs. 21, 25, and 32 using the density-matrix renormalization group (DMRG), the open boundary condition is employed, and hence the dispersion is not available be- cause the momentum is ill-defined. In the present study, we can obtain accurate excitation energies with momentum resolution. The accurate estimate of excitation gaps enables us to perform the level spectroscopy to estimate the phase boundary and elucidate the nature of the QSL phase.
Initial-parameter dependence of the RBM+PP optimization As we describe in Appendix A, we perform several independent optimizations of the RBM+PP wave functions for each J 2 point. Here, using the 8 × 8 lattice, we show how the difference in initial variational parameters affects the optimization. Figure 12 shows the initial-parameter dependence of the RBM+PP ground-state optimization curves for J 2 = 0.5 and J 2 = 0.6. For J 2 = 0.5, we see that four independent optimizations converge to the same energy stably. On the other hand, at J 2 = 0.6, the RBM+PP wave function whose optimization curve is shown in blue color seems to be trapped in a local minimum. The green curve is also trapped at similar energy (dotted line), but it eventually gets out of the local minimum. The behavior seen at J 2 = 0.6 can be understood from the proximity to the 1st-order transition point around J 2 = 0.61 between the VBS and stripe-AF phases. At large system sizes, there exists an energy-level crossing between the states belonging to the same quantum-number sector (zero total momentum and singlet), which gives a kink in the J 2 dependence of the ground state energy (Fig. 13). Therefore, different solutions are competing in small energy scale in the same quantum-number sector at J 2 = 0.6, which makes the optimization more unstable as compared to that at J 2 = 0.5.
One of the reasons for the stable optimization is that finite-size systems we have treated [N site is O(100)] have a finite energy level spacing except for level crossing points. The order of the energy level spacing is on the order of 0.1 (see, e.g., Fig. 3), and the RBM+PP method has a finer energy resolution (note that the energy axis scale of Fig. 12 is 0.1). The results in Fig. 12 suggest that, though the optimized variational parameters may depend on the initial parameters, the optimized wave functions themselves are essentially identical (we have confirmed this by calculating the overlap using the Monte Carlo method among the optimized wave functions).
From this benchmark, we notice that it is important to perform several independent optimizations to avoid being trapped in local minima. In the present study, although the optimizations of the RBM+PP wave functions are done independently for different J 2 points, thanks to the several independent optimizations at each J 2 value, all the physical quantities change smoothly and continuously. The phase transition between the VBS and stripe-AF phases at J V-S 2 in Fig. 1 is of 1st order. To see this, we show the ground state energy as a function of J 2 in Fig. 13. As the system size increases, we see a clear kink S s (q) q x / π q y / π q y / π q x / π q x / π q x / π q y / π q y / π FIG. 14. Structure factor for spin-spin correlation Ss(q).
in the energy curve at J V-S 2 ≈ 0.61, giving evidence for the 1st-order phase transition.

Structure factors
In Sec. IV A 1, we discuss the crossing of the correlation ratio. The correlation ratio quantifies how sharp the structure factor peak is. In Figs. 14 and 15, we show the raw data of the structure factors for spin-spin and dimer-dimer correlations, respectively, which are used in the correlation ratio analysis.
System-size dependence of the crossing J2 points of the AF and VBS correlation ratios As described in Sec. IV A 1, we determine the AF-QSL and QSL-VBS phase boundaries from the correlation ratio analysis. Figure 16 shows the system-size dependence of the crossing points of the correlation-ratio curves. We see that the system-size dependence is small. The fits of the system-size dependence with a + b/L dependence give the estimates of AF-QSL and QSL-VBS phase boundaries as J Néel Excitation gap at 12 × 12 lattice -sublattice-size dependence in the PP part As we mentioned in Appendix A, we impose the 4 × 4 sublattice structure in the f ↑↓ ij parameters in the PP part. With this setting, we have momentum resolution of 4 × 4 K points: K = (mπ/2, nπ/2) with m, n = −1, 0, 1, 2. To investigate the sublattice-size dependence, for 12 × 12 lattice, we also calculate the excitation energies using 6 × 6 sublattice structure. Then, we can calculate the excitation gaps at K = (mπ/3, nπ/3) with m, n = −2, −1, 0, 1, 2, 3. Figure 17 shows the f ↑↓ ij -sublattice-size dependence of the excitation energies. We see that the excitation gaps at high-symmetry K points [(0, 0), (π, 0), and (π, π)] show good agreement between the 4 × 4 and 6 × 6 sublattice structures. At the intermediate K points, the excitation energies stay larger than those at high-symmetry K points. This fact supports the scenario of Dirac-type nodal QSL. System-size dependence of the crossing points of the correlation ratio for spin-spin (red dots) and dimer-dimer (black squares) correlations, which are used to determine the phase boundary of Néel-AF and VBS, respectively. We focus on the crossing of the curves between the L1 × L1 and L2 × L2 lattices with (L1, L2) = (8, 12), (8,16), (12,16) and (16,18) for spin-spin correlations, and (L1, L2) = (8, 12), (8,16), and (12,16) for the dimer-dimer correlations. L mid is defined as L mid = (L1 + L2)/2.