Comment on "Gapless spin liquid ground state of the spin-$\frac{1}{2}$ $J_1$-$J_2$ Heisenberg model on square lattices"
Bowen Zhao Jun Takahashi Anders W. Sandvik



 
bwzhao@bu.edu 
Department of Physics, Boston University, 590 Commonwealth Avenue, Boston, Massachusetts 02215, USA

 
Department of Physics, Boston University, 590 Commonwealth Avenue, Boston, Massachusetts 02215, USA
Beijing National Laboratory for Condensed Matter Physics and Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China

 
sandvik@bu.edu 
Department of Physics, Boston University, 590 Commonwealth Avenue, Boston, Massachusetts 02215, USA
Beijing National Laboratory for Condensed Matter Physics and Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China



Liu et al. [Phys. Rev. B 98, 241109 (2018)] used Monte Carlo sampling of the physical degrees of freedom of a Projected Entangled Pair State (PEPS) 
type wave function for the S=1/2 frustrated J_1-J_2 Heisenberg model on the square lattice and found a non-magnetic state argued to be a gapless spin 
liquid when the coupling ratio g=J_2/J_1 is in the range g ∈ [0.42,0.6]. Here we show that their definition of the order parameter for another candidate 
ground state within this coupling window—a spontaneously dimerized state—is problematic. The order parameter as defined will not detect dimer order when 
lattice symmeties are broken due to open boundaries or asymmetries originating from the calculation itself. Thus, a dimerized phase for some range of g 
cannot be excluded (and is likely based on several other recent works).


Comment on 
 "Gapless spin liquid ground state of the spin-1/2 J_1-J_2 Heisenberg model on square lattices"
    Anders W. Sandvik
    December 30, 2023
===========================================================================================================



§ OVERVIEW
 
		
In a recent Rapid Communication <cit.>, Liu et. al. argued that there is a gapless spin liquid phase in the ground state of the S=1/2 frustrated 
square-lattice J_1-J_2 Heisenberg model for g=J_2/J_1  ∈ [0.42,0.6]. At variance with other recent works <cit.>, 
they found no  spontaneously dimerized valence-bond-solid (VBS) phase within this range of coupling ratios (where other works have roughly placed the VBS
at g ∈ [0.52-0.61]). They reached their conclusions based on the method of Monte Carlo sampling of gradient-optimized tensor network states 
<cit.>, which they have further refined for the specific case of a tensor network of the Projected Entangled Pair State (PEPS) 
type. Open-boundary lattices with up to 16× 16 spins were used, and, taken at face value, the results appear to be well converged and reliable.
			    
In this Comment we point out that the definition of the VBS order parameter used by Liu et al. has a potential flaw and may not capture long-range order 
correctly on the open-boundary lattices considered. The squared columnar VBS order parameters for x and y oriented dimers, m_dx^2 and m_dy^2, 
were defined in Eq. (2) of <cit.> as follows (up to typographical errors):

    m_d α^2 = 1/N_b^2∑_𝐫𝐫' e^i 𝐪_α· (𝐫 - 𝐫')
    (⟨ B^α_𝐫 B^α_𝐫'⟩ - ⟨ B^α_𝐫⟩⟨ B^α_𝐫'⟩),

where α = x, y, B^α_ r= S( r) · S( r + α̂) is the bond operator along the α direction,
and N_b is the number of bonds summed over. The wave-vector corresponding to columnar order is 𝐪_α =  (π, 0) and (0, π) for 
α = x and α = y, respectively. We can rewrite this squared order parameter in the equivalent form:

    m_d α^2 = ⟨ D_α^2 ⟩ - ⟨ D_α⟩^2,

where

    D_α = 1/N_b∑_𝐫 e^i 𝐪_α·𝐫𝐒(𝐫) ·𝐒(𝐫 + α̂).

The problem with the definitions is the subtraction of the non-uniform ⟨ B^α_𝐫⟩⟨ B^α_𝐫'⟩ in 
Eq. (<ref>) or ⟨ D_α⟩^2 in Eq. (<ref>) when long-range order is induced by some symmetry-breaking mechanism, e.g., 
with certain open lattice boundaries or some imperfection in the method used. In essence, the baby is then thrown out with the bath water.

We will demonstrate this problem by considering a columnar VBS state which is four-fold degenerate on periodic L × L lattices with even 
L. The ground state is uniform in the absence of some symmetry-breaking mechanism, and the subtracted term ⟨ D_α⟩^2 in 
Eq. (<ref>) vanishes. However, on rectangular lattices with L_x × L_y spins (even L_x and L_y) the ground state is unique and hosts
a specific dimer pattern. The two terms then cancel each other in the limit of large system sizes, thus rendering the definitions Eq. (<ref>)
and Eq. (<ref>) unsuitable for detecting the dimer order. On square L × L lattices there is a two-fold symmetry left, and the definitions 
can in principle detect the dimerization (albeit with a reduced value of the order parameter). However, in practice the calculation itself may break the 
90^∘ lattice rotation symmetry, and then again the definition is not suitable. Due to the likely symmetry breaking of the PEPS calculations in 
Ref. <cit.>, a VBS phase in the J_1-J_2 Heisenberg model cannot be ruled out based on the results presented.

In the following we will use a specific example of a quantum spin model with a well established columnar VBS phase to illustrate our arguments; the S=1/2
square-lattice J-Q_3 model <cit.> 
defined by the Hamiltonian

    H = - J ∑_⟨ jl ⟩ P_ij - Q_3  -2mm∑_⟨ ijklmn ⟩-2mm P_ij P_kl P_mn.

Here P_ij=1/4- S_i · S_j is a singlet projector and the first term in Eq. (<ref>) is the standard antiferromagnetic Heisenberg 
exchange between nearest neighbor spins. In the second term, the three index pairs ij, kl, and mn correspond to parallel links forming 
columns  on 3× 2 and 2× 3 lattice cells. This correlated singlet interaction leads to the formation of a four-fold degenerate columnar 
VBS with a spontaneous Z_4 symmetry-breaking transition at a critical value of Q_3/J. Here we will consider the case J=0 and focus on the detection 
of the columnar order on open lattices, using a ground-state valence-bond projector quantum Monte Carlo (QMC) method <cit.> with which the 
spin-rotation invariant bond correlations in Eqs. (<ref>) and (<ref>) can be evaluated easily. Though we use a different model and a 
different method for obtaining the ground state, the order parameter is the same as in Ref. <cit.>, and the problem of subtracting off a boundary-induced 
expectation value when the symmetry is broken is exactly the same. In a previous work by one of us <cit.>, related issues were discussed 
in the context of cylindrical lattices (often used in DMRG calculations <cit.>), i.e., ones with periodic boundary conditions in one direction 
and open boundaries in the other direction. Here we focus specifically on the problems with the order parameter definitions in Eqs. (<ref>) and 
(<ref>) when all boundaries are open, as is more practical in PEPS calculations.

In Sec. <ref> we first consider L_x× L_y lattices with L_x=2L_y, for which the VBS pattern is unique (for L_y an even number) 
and our arguments can be illustrated most clearly. In Sec. <ref> we consider the slightly more subtle case of L× L lattices (again with L 
even), on which the columnar VBS pattern is two-fold degenerate. We summarize our conclusions and discuss implications in Sec. <ref>.





§ RECTANGULAR LATTICES



As mentioned above, on an open 2L× L lattice with even L there is a unique columnar VBS pattern in the ground state. In the case
of the Q_3 model, J=0, Q_3=1 in Eq. (<ref>), the boundaries favor dimers (bonds with a higher singlet density) perpendicular 
to the edges, as illustrated for a 32× 16 system in Fig. <ref>. The dimer orientation favored by the longer 
edge survives in the center of the system, and it is  possible to use the dimer order parameter D_y defined in Eq. (<ref>) without 
squaring (as noted previously, e.g., in Ref. <cit.>). 

To demonstrate that correct results are obtained with D_y in the thermodynamic limit, in Fig. <ref> we compare results for ⟨ D_y⟩^2 
computed on the central L× L/2 bonds of 2L× L lattices (to eliminate some of the boundary enhancements of the order, though this is
not necessary) with results for ⟨ D^2⟩ = ⟨ D_x^2⟩ + ⟨ D_y^2⟩ calculated on periodic L× L lattices.
We have fitted both data sets using exponentially convergent forms, as expected for VBS order <cit.>, but details of the fits are not 
important here. The extrapolated, clearly non-zero values are compatible with each other. In contrast, using the definition m_dy^2 by Liu et al., 
Eq. (<ref>), gives results approaching zero with increasing system size (we only show rough fits, but the trend is
clear), as expected when the symmetry-broken order parameter has been subtracted off. The x-oriented order parameters should of course vanish, 
on account of the rectangular lattice shape inducing only y columnar order in the thermodynamic limit.





§ SQUARE LATTICES



As an example more closely corresponding to the calculations in Ref. <cit.>, we next consider the same Q_3  model as above
but on L × L lattices (even L), again with all open boundaries. Since now there is no anisotropy between the x and y directions, the true ground 
state does not have a unique locked-in dimer pattern, but is two-fold degenerate with fluctuations between x- and y-oriented order. On a small lattice,
these fluctuations are fully sampled in our QMC simulations running for reasonable times, and when averaged the bond patterns look more like a
plaquette VBS state. This is shown in Fig. <ref>(a) for a 32× 32 lattice. Here it should be noted that the dimers at
the boundaries do not fluctuate much, and the central part of the system can be regarded as a kind of domain-wall state with de-facto plaquette
order. For system sizes larger than the domain wall thickness, the bonds in the center of the system fluctuate collectively between actual long-range 
ordered horizontal and vertical bond patterns. For very large lattices, the time scale of these fluctuations between the two different bond order 
realizations becomes too long (diverging exponentially with the system size) to observe in simulations, and the system may in practice become completely 
trapped in one of the sectors. This is seen in Fig. <ref>(b) for an 80 × 80 system, where the central part of the system only 
exhibits strong y dimers.



To further illustrate this symmetry breaking occurring in the simulations, 
in Fig. <ref> we show the probability distribution P(D_x, D_y) of the dimer order parameter as 
collected in the QMC process. For the smallest system, L=32, the peak in the distribution corresponds to equal D_x and D_y, i.e., resonating plaquette 
order or equal amounts of static x and y dimers. For a slightly larger system, L=48, we observe the peak splitting into two, indicating a state that is
now fluctuating between x and y oriented bond order. The splitting between the peaks grows with increasing system size as the two patterns become
more dominated by the majority order, and the tunneling probability decreases (reflected in smaller weight close to the line D_x=D_y). For L=64, the 
two peaks have unequal density, due to the long time scale of fluctuation for this system size, and for the lagest system, L=80, the simulation was 
completely locked into the y sector. As is typical in systems with a discrete symmetry of the order parameter, the time scale of tunneling between
sectors should grow exponentially with increasing L, and in practice it is not possible to sample equally both sectors for large systems. 




The method-related symmetry breaking is not a problem in practice, as long as computed quantities are insensitive to the symmetry breaking, e.g., with 
the definition ⟨ D_x^2⟩ + ⟨ D_y^2⟩ of the VBS order parameter. As shown in Fig. <ref>, results based on
this definition for the open system agrees with those for periodic boundary conditions in the limit L →∞, though the extrapolation to
infinite size is easier for the periodic systems. On the open systems the behavior is non-monotonic.
We also show results for m_dx^2 + m_dy^2, based on the definitions by Liu et al. in 
Eq. (<ref>). In this case we see a sharp change in the behavior at a system size corresponding the the de-facto symmetry breaking of the QMC 
simulations for system sizes above L=64. For the smaller systems, the results appear to extrapolate to a non-zero value, but for the larger sizes 
the values drop rapidly toward zero. The latter behavior reflects the cancelation of the terms in the order-parameter definition, Eq. (<ref>), 
when the ground state is unique (in practice, due to the broken symmetry). For the smaller sizes the cancelation is not complete because of the two-fold 
degeneracy. A similar discontinuous behavior arising from symmetry breaking is seen in a symmetric definition based on the induced order parameter
(i.e., squaring the components after the mean value has been computed), ⟨ D_x⟩^2 + ⟨ D_y⟩^2, where
the results for the larger systems exhibit a jump up toward the results for periodic boundary conditions when the symmetry breaking takes place.
Before symmetry breaking we have ⟨ D_x⟩ = ⟨ D_y⟩, and, because the squares are taken, a value 1/2 of 
⟨ D_x^2⟩ + ⟨ D_y^2⟩ obtains (in the limit L →∞ without symmetry breaking).

A properly symmetrized version of the definition (<ref>), in its equivalent form  (<ref>), is 
⟨ D_x^2⟩ + ⟨ D_y^2⟩ - 1/2(⟨ D_x⟩ + ⟨ D_y⟩^2. In Fig. <ref> it can
be seen that the results for this quantity concide with m_dx^2 + m_dy^2 when there is no symmetry breaking, while after symmetry 
breaking the two definitions diverge sharply. The properly symmetrized definition should be 3/4 of the standard squared VBS order parameter 
for periodic boundaries in the thermodynamic limit.





§ CONCLUSION



We have discussed why the quantities m_dx^2 and m_dy^2 [Eq. (<ref>)] used in Ref. <cit.> 
may not capture long-range VBS order properly. It is clear that, in a system where the VBS order parameter symmetry is fully broken, the terms 
subtracted in Eq (<ref>) correspond to the actual order parameter of interest, and what is left vanishes for large distances (large systems).
In calculations with tesor network states, such as the PEPS used in Ref. <cit.>, symmetry breaking can take place due to unequal treatment of 
the x and y directions or imperfect optimization (even on periodic lattices). As in Monte Carlo simulations, which may be trapped in one out of 
two or more sectors of the order parameter, this kind of “artificial” symmetry breaking may not be a problem in practice, as long as the consequences 
are understood and taken into account properly.

A VBS phase in the J_1-J_2 Heisenberg model cannot be excluded by the results 
presented in Ref. <cit.>. Judging from other recent calculations with a variety of methods <cit.>, we expect VBS order 
in a narrow range of coupling ratios g=J_2/J_1 (roughly for g ∈ 0.52,0.61]). According to the same calculations, a gapless spin liquid may exist 
for slightly smaller values of g (roughly for g ∈ 0.45,0.52]). It would be very interesting to see the VBS order parameter from the calculations
in Ref. <cit.> without the subtraction of the crucial boundary induced contributions, as well as the boundary-induced order parameter itself.
We also point out that it may be advantageous to use rectangular lattices in PEPS calculations, as is evident from the behavior of ⟨ D_y⟩ 
in Fig. <ref>.

Aside from the use of a potentially flawed 
VBS order parameter, the calculations in Ref. <cit.> are impressive and suggest that the method of Monte Carlo 
sampling of the physical degrees of freedom and gradient-based optimization <cit.> may indeed be one of the most powerful 
ways to compute with tensor-netweok states. Very recently, further progress along these lines were reported in the context of the same J_1-J_2 
Heisenberg model up to system size 24× 24 <cit.>. The VBS order was not discussed, however.

0mm
ACKNOWLEDGMENTS
0mm

We would like to thank Wenyuan Liu for discussions.
This work was supported by the NSF under Grant No. DMR-1710170 and by a Simons Investigator Award.
The numerical calculations were carried out on Boston University's Shared Computing Cluster.


99

Liu18
W.-Y. Liu, S. Dong, C. Wang, Y. Han, H. An, G.-C. Guo, and L. He, 
Gapless spin liquid ground state of the spin-1/2 J_1-J_2 Heisenberg model on square lattices
Phys. Rev. B 98, 241109 (2018).

Gong14
S.-S. Gong, W. Zhu, D. N. Sheng, O. I. Motrunich, and M. P. A. Fisher, 
Plaquette ordered phase and quantum phase diagram in the spin-1/2 J_1-J_2 square Heisenberg model, 
Phys. Rev. Lett. 113, 027201 (2014).

Morit15
S. Morita, R. Kaneko, and M. Imada, 
Quantum spin liquid in spin-1/2 J_1-J_2 Heisenberg model on square lattice: Many-variable variational Monte Carlo study 
combined with quantum-number projections, J. Phys. Soc. Jpn. 84, 024720 (2015).

Wang18
L. Wang and A. W. Sandvik, 
Critical level crossings and gapless spin liquid in the square-lattice spin-1/2 J_1-J_2 Heisenberg antiferromagnet,
Phys. Rev. Lett. 121, 107202 (2018).

Haghshenas18
R. Haghshenas and D. N. Sheng, 
U(1)-symmetric infinite projected entangled-pair state study of the spin-1/2 square J_1-J_2 Heisenberg model, 
Phys. Rev. B 97, 174408 (2018).

Sandvik07
A. W. Sandvik and G. Vidal, 
Variational quantum Monte Carlo simulations with tensor-network states,
Phys. Rev. Lett. 99, 220602 (2007).

Schuch08
N. Schuch, M. M. Wolf, F. Verstraete, and J. I. Cirac, 
Simulation of quantum many-body systems with strings of operators and Monte Carlo tensor contractions,
Phys. Rev. Lett. 100, 040501 (2008).

Wang11
L. Wang, I. Pižorn and F. Verstraete, 
Monte Carlo simulation with tensor network states,
Phys. Rev. B 83, 134421 (2011).

Mambrini06
M. Mambrini, A. Läuchli, D. Poilblanc, and F Mila,
Plaquette valence-bond crystal in the frustrated Heisenberg quantum antiferromagnet on the square lattice,
Phys. Rev. B 74, 144422  (2006). 

Lou09
J. Lou, A. W. Sandvik, and N. Kawashima, 
Antiferromagnetic to valence-bond-solid transitions in two-dimensional SU(N) Heisenberg models with multispin interactions,
Phys. Rev. B 80, 180414 (2009).

Sandvik12
A. W. Sandvik, 
Finite-size scaling and boundary effects in two-dimensional valence-bond solids,
Phys. Rev. B 85, 134407 (2012).

Sandvik10
A. W. Sandvik and H. G. Evertz, Phys. Rev. B 82, 024407 (2010). 

Sandvik02
A. W. Sandvik, S. Daul, R. R. P. Singh, and D. J. Scalapino,
Striped phase in a quantum XY model with ring exchange,
Phys. Rev. Lett. 89, 247201 (2002).

Liu19
W.-Y. Liu, Y.-Z. Huang, S.-S. Gong, and Z.-C. Gu,
Accurate simulation for finite projected entangled pair states in two dimensions,
arXiv:1908.09359.




