{'fraction_non_alphanumeric': 0.04433815230827108, 'fraction_numerical': 0.03192346966195518, 'mean_word_length': 4.206491326245104, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spin-1/2 J1-J2 Heisenberg antiferromagnet on a square lattice: a plaquette renormalized tensor network study

We apply the plaquette renormalization scheme of tensor network states [Phys. Rev. E, 83, 056703 (2011)] to study the spin-1/2 frustrated Heisenberg J1-J2 model on a L*L square lattice with L = 8, 16 and 32. By treating tensor elements as variational parameters, we obtain the ground states for different J2/J1 values, and investigate the staggered magnetic order parameters, the nearest-neighbor spin-spin correlations and the plaquette order parameters. In addition to the well-known N\'eel-order and collinear-order at low and high J2/J1, we observe a plaquette order at intermediate J2/J1. A continuous transition between the N\'eel order and the plaquette order near J2 = 0.40J1 is observed. The collinear order emerges at J2 = 0.62J1 through a first-order phase transition.


I. INTRODUCTION
The search for exotic states in quantum magnets has been the topic of intensive research for the past decades. An extremely important question is when the conventional Néel order is destroyed, what kind of states can emerge. Frustrated antiferromagnetic spin systems, where the frustration from either the lattice geometry, or the presence of competing interactions, are candidate systems to study these states. It is proposed that when the Néel order is destroyed by quantum fluctuations, only short-range correlations will survive, and the system enters a quantum paramagnetic state which can be described as a resonant valence bond (RVB) state. 1 The RVB state can either be a valence bond solid (VBS) phase, where some of the lattice symmetries are broken, 2 or a featureless spin liquid with strong shortrange correlations without any broken spin symmetry. 3,4 One archetypical model to study the effect of frustration from competing interactions is the antiferromagnetic (AF) J 1 -J 2 Heisenberg model on a square lattice. 3,[5][6][7][8][9][10][11][12][13][14][15] The Hamiltonian is given by, where J 1 > 0 and J 2 > 0 are the nearest-neighbor (NN) and next-nearest-neighbor (NNN) couplings, and the sums ij and ij run over NN and NNN pairs, respectively. Recent interests of this model have been revived by the discovery of Fe-based superconducting materials 16 where a weakened AF order can be described by this model with S > 1/2. [17][18][19] Properties of this model for S = 1/2 in 2d have been studied extensively by a variety of methods, such as spin wave theory, 5 exact diagonalization(ED), 6,7,14 series expansion, 12,[20][21][22][23] large-N expansion, 2 functional renormalization group, 15 Green's function method, 11 projected entangled pair states, 24 etc. It is generally believed that in the region J 2 /J 1 0.4, the ground state (GS) of the model is the Néel phase with magnetic long-range order (LRO). In the region J 2 /J 1 0.65, spins in the GS are ordered at wave vector (π, 0) or (0, π), showing so-called collinear magnetic LRO. The GS in intermediate region is proposed to be a quantum paramagnet without magnetic LRO, but the properties of this phase are still under intensive debate. There are several proposals for the GS, such as a columnar dimer state, 21,24 a plaquette VBS order, 8,13,25 or a spin-liquid. 3,4 In the mean time, precise determination of the phase transition points is also not conclusive. Earlier series expansion studies 21 estimate the quantum paramagnetic region is between 0.38 J 2 /J 1 0.62. Recent ED study 14 using results of up to N = 40 to perform finite-size extrapolation estimates the transition points at J c1 2 ≃ 0.35J 1 and J c2 2 ≃ 0.66J 1 . Meanwhile, studies by combination of random phase approximation and functional renormalization group find this nonmagnetic phase begins near J 2 /J 1 ≈ 0.4 ∼ 0.45 and ends around 0.66 ∼ 0.68. 15 Numerical studies of frustrated quantum spin systems present great challenges in dimensions greater than one. The ED method is hampered by the limitation of system size one can simulate. At present, the largest system size on the square lattice that can be simulated is N = 40. 14,26 Due to the minus sign problem, 27 the powerful quantum Monte Carlo (QMC) method is not applicable to highly frustrated systems. In 1d, the density matrix renormalization group (DMRG) 28 algorithm, which generates matrix product states (MPS), can reach very high accuracy even for frustrated spin systems; however, direct extension of the algorithm to higher dimensions remains difficult. One promising proposal is to generalize the MPS to higher dimensions, the tensor network states (TNS), 29-31 which can serve as potential candidates for studying these systems. In the TNSs, the matrices are replaced by tensors of rank corresponding to the coordination number of the lattice. On a 2d square lattice, the tensor T s ijkl (σ s ) on site s has four indices, in addition to the physical index, which in the current case corresponds to the z-component σ s of a spin.
Here, we should mention, according to the TNS representation, the rank of tensors is chosen according to the coordination number instead of the interaction pattern. In this way, the area law of entanglement entropy can be satisfied well if bond dimension D is big enough, especially when J 2 not very large.
Contracting over all bond indices gives the wave function coefficient for a given spin state σ 1 , . . . , σ N . [32][33][34] In these tensor network based methods, one of the major obstacles is the computational complexity involved in the tensor contraction, then usually some type of approximation is required to make the computation manageable. Several schemes have been proposed to facilitate the contraction of the tensor networks. [32][33][34][35][36] In particular, a contraction scheme based on the plaquette renormalization with auxiliary tensors is proposed to retain the variational nature of the method, and it is shown that for the transverse Ising model, even with the smallest possible bond dimension (D = 2), non-mean-field results can be obtained. 34 In this paper, we use the TNS with the plaquette renormalization scheme to study the J 1 − J 2 Heisenberg model on a square lattice. We find that even with a small bond dimension D = 2, it already provides a useful way to study the nature of the transition and estimate the value of the transition points. The rest of this paper is organized as follows. In the following section, we review the plaquette renormalization scheme of TNS, and how to apply the scheme to the current model. Main results will be presented in Sec. III, as well as some discussions. Sec. IV will give a brief summary.

II. METHOD
We investigate the ground state of frustrated Heisenberg J 1 -J 2 model on a square lattice, using the plaquette renormalized tensor network 34 . The trial wave function is written as where tT r indicates the tensor trace that all the tensor indices are summed over. T s is rank-4 tensor on site s, with bond dimension D for each rank and σ s =↑ or ↓ is the physical spin state. Explicit contraction of the tensor network is computationally intensive. To keep the computational complexity from growing exponentially, auxiliary rank-3 tensors A n ijk are added to each level of the contraction process A sequence of plaquette renormalizations, n = 1, 2, . . ., is carried out and the bond dimension of each rank is thus kept constant after every plaquette contraction. 34 In order to compute physical expectation values based on a TNS, one has to contract the tensors of a bra and ket state over their physical (e.g., spin) indices in addition to the bond indices of the tensors. Normally, one would first construct the double tensors by performing the sum over the physical indices, where the labels a, b, c, d is a suitable combination of the indices of the bra (T s * ) and ket (T s ) tensors, i.e., a = i 1 + D(i 2 − 1), etc. In the calculation of the matrix element Ψ|Ô|Ψ of some operator involving one or several sites, similar tensors are constructed for the sites at which operators act weighted with a local expectation value σ ′ s |Ô s |σ s . In addition, the renormalization double tensors can be also formed The bond dimension of each rank in the resulting double tensor becomes D = D 2 . This renormalization scheme reduces the maximum computational complexity 34 to D 8 = D 16 for a double tensor network. The ground state wave function can be obtained by optimizing the elements of tensors T, A for the ground state energy. Since the plaquette renormalization is introduced at the wave function level, instead of the constructed double tensor network, the method remains variational and the final energy will give a upper bound for the true ground state energy. We optimize the wave function using the derivative-free Brent's method. 37 Compared to previous methods involving singular value decomposition (SVD), 32,33 the environment of a given tensor is fully taken into account in the current scheme. However, the introduction of the renormalization A tensors at the wave function level effectively reduces the maximum support of the entanglement entropy area law in this tensor network. To reduce the number of free parameters, we impose symmetries on the trial wave function. We use a single plaquette, i.e. 2 × 2 = 4 sites as a unit cell (Fig. 1), wherein tensors T on each site and auxilliary tensors A 0 are assumed to be different. This unit is translated to generate a 4 × 4 unit and another set of auxilliary tensors A 1 are added. This procedure is repeated until the full lattice is generated. Finally, the periodic boundary condition is applied. 34

III. RESULTS AND DISCUSSIONS
We obtain the ground state wave function by varying the elements in the tensors T and A with D = 2, which describes a slightly entangled state beyond the product (mean-field) state (D = 1). Figure 2(a) shows the ground state energy with system sizes L = 8, 16, and 32. A clear cusp near J 2 /J 1 = 0.62 is observed, signaling a first-order phase transition. A continuous change of the slope is found near J 2 /J 1 = 0.4, probably indicating a continuous phase transition there.
To study the details of the magnetic orders and the transition points, we compute the magnetic structure factor, or the square of staggered magnetization at wave vector q, defined as where r i = (x i , y i ), and q = (π, π) for the Néel order, and (0, π) or (π, 0) for the collinear order. M 2 (q) tends to the square of the order parameter in the thermodynamic limit if there is magnetic ordering at wave vector q, and scales like 1/N in a magnetically disordered phase. Figure 2(b) shows the results of the square of staggered magnetizations M 2 (π, π) and M 2 (π, 0). From the small J 2 /J 1 side, the Néel order is smoothly suppressed as J 2 increases, until J 2 /J 1 ≃ 0.40, where a discontinuous jump of the Néel order is observed for L = 8, and the jumps become less pronounced as the system size increases. This strong size dependence of the jump is another example that in a finite-size tensor network state with finite bond dimensions, there exists two energy minima near the transition, rendering the transition first-order at small N . For a putative continuous transition, these two minima move closer to each other with increasing N and the transition becomes continuous at N → ∞. 38 From the large J 2 /J 1 side, the collinear order also decreases smoothly, until J 2 /J 1 ≃ 0.6 where a clear firstorder transition occurs. Unlike the previous case, the jumps in M 2 (π, 0) remain robust upon increasing N , strongly suggesting against a continuous transition here. This transition to the collinear order is consistent with previous numerical calculations. 12,14,[20][21][22][23] We now use our data from different sizes to extract the order parameters in the thermodynamic limit. This allows us to estimate the transition points between the Néel/collinear state and the non-magnetic (disordered) phase. The finite-size extrapolation rules for the two-dimensional antiferromagnetic Heisenberg model are well-known [39][40][41] . Following Refs. 14 and 41, we define the Néel order parameter as m 0 = 2 lim N →∞ M (π, π). This normalization is chosen so that m 0 = 1 in a perfect Néel state. The finite-size behavior of M 2 (π, π) is given by, 14,41 where c is the spin-wave velocity and ρ is the spin stiffness. The order parameter for the collinear order is defined as m 1 = √ 8 lim N →∞ M (π, 0). The finite-size behavior of M (π, 0) is given by, 14,41 The extra 1/2 factor comes from the fact that the ground state has an extra two-fold degeneracy q = (π, 0), (0, π), and this symmetry is broken in the thermodynamic limit. Figure 3(a) shows the extrapolated results for m 0 and m 1 as a function of J 2 /J 1 . We find that the GS near J 2 /J 1 = 0.5 is magnetically disordered, i.e., both m 0 and m 1 vanish. Figure 3(b) shows the finite-size scaling of M 2 (π, π) and M 2 (π, 0) at J 2 /J 1 = 0.5, which both shows a 1/N scaling with the zero intercept as N → ∞. The transition points are estimated to be J c1 2 = 0.40J 1 and J c2 2 = 0.62J 1 , consistent with estimates from series expansion 12,[20][21][22][23]42 where J c1 2 ≈ 0.38J 1 and J c2 2 ≈ 0.62J 1 , and slightly different from ED results J c1 2 ≈ 0.35J 1 and J c2 2 ≈ 0.66J 1 . 14 Near J c1 2 , we fit the Néel order parameter m 0 to a power law m 0 ∼ (J 2 − J c1 2 ) β , and an asymptotic mean-field behavior consistent with β = 1/2 is also observed. 38 For J 2 = 0, we obtain m 0 = 0.592 which is slightly lower than the best estimate from the quantum Monte Carlo (m 0 = 0.6140). 43 Although it is also possible to extract c and ρ from our data based on Eq. (6), it is argued that determination of these quantities by fitting the prefactors of the leading finite-size corrections (O(1/L)) can not reach the same accuracy as the magnetic order parameters. 14 Analogous to how mean-field theory produces symmetry-broken states, this method can produce solutions which break spin-rotation symmetry on a finite lattice. 38,43 We examine the spin-rotation symmetry of the ground state, with the focus in the nonmagnetic phase. Figure 4 shows z and xy components of the square of staggered magnetization at q = (π, π) for L = 32, defined as For reference, the sum of the two is also included. In the Néel phase, the spin-rotational symmetry is clearly broken. 43 Increasing J 2 through a phase transition to the strongly frustrated regime (i.e., 0.45 J 2 /J 1 0.60), the spin-rotation symmetry is restored with M 2 z = 1 2 M 2 xy = 1 3 M 2 , as expected. In order to clarify the possible new phase in the highly frustrated region around J 2 /J 1 = 0.5, we calculate the nearest-neighbor spin-spin correlations for L = 32. Figure 5 shows the results for J 2 /J 1 = 0.10, which is deep inside the Néel phase, and J 2 /J 1 = 0.50, which is in the magnetically disordered phase. The numbers in black near the bond are the NN spin-spin correlation, and the thickness of the bond is proportional to its magnitude. For J 2 /J 1 = 0.50 [ Fig. 5(b)], the NN spin-spin correlations within a single plaquette are much stronger than those between plaquettes. On the other hand, deep inside the Néel phase J 2 /J 1 = 0.10 [ Fig. 5(a)], the NN spin-spin correlations shows a more uniform pattern, although weaker correlations are present in some bonds between plaquettes. Overall, it is clear that the correlations inside a 2 × 2 plaquette become stronger upon increasing J 2 /J 1 , which indicates a possible plaquette order in the magnetically disordered phase.
We also investigate the plaquette order parameter, which distinguishes clearly a Néel ordered phase from a plaquette order, defined as 24 The results of the plaquette order parameter are shown also in Figs. 5 (numbers in red italic) for J 2 /J 1 = 0.10 and 0.50. In the most frustrated region, we observe signature of the plaquette order. For J 2 /J 1 = 0.50, the plaquette order parameter is much stronger within a plaquette, consistent with observation from the spin-spin correlations. This order parameter is small in Néel phase (J 2 /J 1 = 0.10), although some traces of the plaquette order is still present. This might be due to the inherent structure of the renormalization scheme, which explicitly breaks the translational invariance, or possibly the plaquette correlations already start to build up in this regime. It remains to further explore whether this plaquette order is favored due to our renormalization scheme. The plaquette renormalization scheme reduces the amount of entanglement support between plaquettes by a factor of D compared with the exact contraction. This may bias toward those correlations compatible with the plaquette structure.

IV. CONCLUSION
We use the plaquette renormalization scheme to study spin-1/2 frustrated Heisenberg J 1 -J 2 model on a square lattice with different sizes of L = 8, 16, and 32. Using the smallest possible bond dimension D = 2 for the underlying tensors, we are already able to obtain results beyond the mean-field theory. Since our method is variational, and the calculations are done on finite lattices, we are able to perform finite-size scaling to extrapolate the order parameters in the thermodynamic limit. We observe signatures of a continuous transition at J c1 2 ≃ 0.40J 1 , and a first-order phase transition at J c2 2 ≃ 0.62J 1 , consistent with previous numerical calculations. 14,21 Our calculations on the NN spin-spin correlation and the plaquette order parameter indicates a possible plaquette VBS order for J c1 2 < J 2 < J c2 2 . The effects of the plaquette renormalization scheme and the bond dimension D dependence of the physical observables require further studies and will be presented in a future work. 44,45 V. ACKNOWLEDGEMENTS We thank A. Sandvik for useful conversation and collaboration on related work. We are grateful to National Center for High-Performance Computing Computer and Information Networking Center, NTU for the support of high-performance computing facilities. This work was partly supported by the National Science Council in Taiwan through Grants No. 100-2112-M-002 -013 -MY3, 100-2120-M-002-00 (Y.J.K.), and by NTU Grant numbers 99R0066-65 and 99R0066-68 (J.F.Y., Y.J.K.). Travel support from National Center for Theoretical Sciences is also acknowledged.
note added.-After submitting this manuscript, we recently learned of the DMRG work by Jiang et al. 46 and the tensor product state approach by Wang et al. 47 on the same model, which argue that the ground state in the nonmagnetic regime near J 2 /J 1 ∼ 0.5 could be a Z 2 spin liquid.