{'fraction_non_alphanumeric': 0.05570928623748678, 'fraction_numerical': 0.02791568324245137, 'mean_word_length': 3.8423645320197046, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Ground-State Phase Diagram of Quantum Heisenberg Antiferromagnets on the Anisotropic Dimerized Square Lattice

The S=1/2 and S=1 two-dimensional quantum Heisenberg antiferromagnets on the anisotropic dimerized square lattice are investigated by the quantum Monte Carlo method. By finite-size-scaling analyses on the correlation lengths, the ground-state phase diagram parametrized by strengths of the dimerization and of the spatial anisotropy is determined much more accurately than the previous works. It is confirmed that the quantum critical phenomena on the phase boundaries belong to the same universality class as that of the classical three-dimensional Heisenberg model. Furthermore, for S=1, we show that all the spin-gapped phases, such as the Haldane and dimer phases, are adiabatically connected in the extended parameter space, though they are classified into different classes in terms of the string order parameter in the one-dimensional, i.e., the zero-interchain-coupling, case.


I. INTRODUCTION
Recently many low-dimensional antiferromagnets with excitation modes separated from a ground state by a finite energy gap have been synthesized, and effects of impurities or magnetic fields on those materials have been investigated experimentally in relation to the impurity-induced long-range order (LRO) and magnetic-field-induced LRO. For example, there are the S = 1/2 quasi-one-dimensional (Q1D) Heisenberg antiferromagnet (HAF) with bond dimerization, CuGeO 3 , 1 S = 1 Q1D HAF's, NENP, 2 NDMAZ, 3 NDMAP, 4 and PbNi 2 V 2 O 8 , 5 and S = 1 Q1D HAF's with bond alternation, NTEAP 6 and NTENP. 7 Those materials have attracted our interest since they reveal various aspects of the quantum phase transition between the quantumdisordered spin-gapped phase and the classical (Néel) long-range-ordered phase.
Intrachain spin interaction with or without bond alternation and interchain interaction are considered to be the most basic ingredients to understand the quantum phase transitions mentioned above. More explicitly, they are expected to be modeled effectively by the quantum HAF on the anisotropic dimerized square lattice, which is described by the Hamiltonian: Here S i,j is the quantum spin operator at site (i, j) on the square lattice. The first two terms in r.h.s. represent the one-dimensional antiferromagnetic (AF) Heisenberg chains with alternating coupling constants, 1 and α (0 ≤ α ≤ 1), and the last term does the AF interchain exchange interaction (J ′ ≥ 0). We choose the x-axis as being along the chain direction and the y-axis as in the perpendicular one. The bond arrangement of this model is shown in Fig. 1.
The ground state of decoupled chains, i.e., J ′ = 0, has been well established. In particular, that of the S = 1/2 chain 8 is the dimer state with a finite spin gap except for the uniform case (α = 1), which has a critical ground state. On the other hand, for the S = 1 chain there exist two spin-gapped phases: the Haldane phase 9 at α > α c and the dimer phase at α < α c . 10 At the critical point α = α c between these two phases the gap vanishes. The value of α c has been estimated to be 0.5879 (6). 10,11 In both cases, the critical point is considered to belong to the Gaussian universality class. 10 For the AF LRO to appear the higher-dimensionality effect, i.e., the interchain interaction J ′ , is indispensable. In most of the numerical works reported so far, the effect of interchain coupling has been examined in certain approximated or perturbed ways. For example, Sakai and Takahashi 12 estimated the critical strength, J ′ c , for the uniform case (α = 1) by the exact diagonalization method for the intrachain interactions combined with the mean-field approximation for the interchain interaction, and obtained J ′(S=1/2) c = 0 and J ′(S=1) c ≥ 0.025. More recently, Koga and Kawakami 13 investigated the S = 1 model by the cluster-expansion method, and obtained J ′ c = 0.056(1) for α = 1. However, there have been only a very limited number of numerical works, in which both of the interchain and intrachain interactions are treated on an equal footing. 14,15 Such numerical analyses are certainly required, since the mean-field-like approximation is not necessarily appropriate even in the Q1D regime. 16 In the present paper, we report the results of quantum Monte Carlo (QMC) simulations by using the continuousimaginary-time loop algorithm 17-20 on the S = 1/2 and S = 1 HAF model described by Eq. (1). The present paper is organized as follows. In Sec. II, the method of our numerical analyses is explained. In Secs. III and IV, the ground-state phase diagram parameterized by the strength of the bond alternation, α, and that of the interchain coupling, J ′ , is determined precisely for S = 1/2 and S = 1, respectively. Especially, for the S = 1 system with α = 1, we obtain J ′ c = 0.043648 (8), which is consistent with 0.040(5) suggested by the recent QMC work, 15 but is much more accurate. Furthermore, both in the S = 1/2 and S = 1 systems, the quantum phase transitions between the spin-gapped phases and the AF-LRO phase are confirmed to belong to the same universality class with that of the 3D classical Heisenberg model: the exponent of the correlation lengths is ν = 0.71(3) for S = 1/2 and ν = 0.70(1) for S = 1, which coincides fairly well with that of the latter model, ν = 0.7048 (30). 21 We also show the results on the correlation length and the gap in the spin-gapped phase. In Sec. V, the topology of the phase diagram is discussed in detail based on the result of the present QMC calculation. We show that all the spin-gapped phases, such as the Haldane and dimer phases, are adiabatically connected with each other in the extended phase parameter space. This is in a sharp constrast to the strict 1D case, in which the spin-gapped phases are classified into different classes in terms of the so-called string-order parameter. 23 It is of interest that the 1D spin-gap phases, which have different hidden symmetry, are connected without encountering any singularity in the 2D phase diagram. The final section is devoted to the concluding remarks.

II. METHOD
We consider the system descrived by the Hamiltonian (1) with S = 1/2 and S = 1. The real-space size is L x × L y and the inverse temperature, i.e., the imaginarytime size, is β = 1/T . Periodic boundary conditions are imposed in the x-and y-directions. We use the continuous-imaginary-time loop algorithm with multicluster update. 17,18 Especially, for the S = 1 system we adopt the subspin-representation technique, 19,20 in which the S = 1 system is represented by an S = 1/2 system with special boundary conditions in the imaginary-time direction. By using these techniques, we can perform the simulatation up to L x × L y = 336 × 48 with β = 100 for the S = 1 case without encountering any difficulty.
The imaginary-time dynamical structure factor is defined by where q = (q x , q y ) is the wave-number vector, S z i (t) is the z-component of the spin on site i at imaginary time τ , and · · · is the thermal average. By using S d (q x , q y , ω), the staggered correlation length along the x-direction, ξ x , is then evaluated by the second-moment method, 20,24 The correlation length in the y-direction, ξ y , and that in the imaginary-time (τ ) direction, ξ τ , which is related to the gap ∆ by ∆ = 1/ξ τ , are calculated similarly. Finally the staggered susceptibility, χ s , is evaluated by All the structure factors are calculated by using the improved estimators. 25 The period of 10 2 -10 3 Monte Carlo steps (MCS) is used for thermalization and that of 10 3 -10 5 MCS for the evaluation of physical quantities. Near the critical point (α c , J ′ c ) of the ground-state transition, the correlation lengths diverge as where t is the distance from the critical point and ν is the critical exponent for the correlation length. Here we have put z = 1 assuming the Lorenz invariance. 26 Furthermore, the following finite-size-scaling (FSS) formula 27 holds near (α c , J ′ c ) and T = 0 for systems with the fixed ratio L x : L y : β, and similar ones for ξ y and ξ τ , and Here f and g are scaling functions and γ the exponent for χ s (∼ t −γ ). Note that L x T is put constant in the above equations. We assume a polynomial up to the second order for the scaling functions. By using leastsquares fitting, we obtain the critical point (α c , J ′ c ) and the associated critical exponents ν and γ.
In addition, at some points in the spin-gapped phase, we explicitly evaluate the correlation lengths, ξ x and ξ y , and the gap, ∆, at T = 0 in the thermodynamic limit L x , L y → ∞. For this purpose we extrapolate the simulated data first to the ground state T → 0 and then to the thermodynamic limit L x , L y → ∞.
In Fig. 2 we show the ground-state phase diagram of the S = 1/2 system obtained by the FSS analysis explained in the previous section. As an example of the FSS analysis, we show in Fig. 3 that of ξ τ against α for J ′ = 1 (dotted line in Fig. 2). The aspect ratio of the (2 + 1)dimensional system is taken as L x : L y : β = 1 : 1 : 1. By the least-squares fitting, the exponent ν and the critical coupling α c are estimated as 0.71(1) and 0.31407(5), respectively. Here, the figure in parentheses denotes the statistical error (1σ) in the last digit. We also perform the same analyses on other lines in the α-J ′ plane, whose results are presented by the solid circles in Fig. 2. For example, on the line α = J ′ (dashed line in Fig. 2), we obtain α c = J ′ c = 0.52337(3) and ν = 0.71(3). In the phase diagram we can see that the ground state of the chain (J ′ = 0) is the dimer state with a spin gap except for α = 1, 12,28-30 and that the region of the AF-LRO phase enlarges monotonically as J ′ increases.
Our phase diagram is qualitatively the same as that of Katoh and Imada (KI), 14 but not quantitatively. In particular, we obtain the critical point, α c = 0.31407 (5), on the line J ′ = 1 (dotted line in Fig. 2), which is significantly smaller than their estimate α c = 0.398. More importantly, in the present simulation, the critical exponent ν on the transition points is evaluated as ν = 0.71(1), which is consistent with ν = 0.7048 (30) for the 3D classical Heisenberg model. 21 The similar results have been obtained for the 2D 1/5-depleted HAF model. 22 On the other hand, KI concluded ν = 1. The reason of the these discrepancies might be due to the smallness of the system sizes and the inverse temperature used in the study by KI.

B. Correlation lengths and the gap
We also evaluate explicitly the ground-state correlation lengths and the gap on some points in the dimer phase by using the dynamic structure factors. Unless the points are very close to the critical line, these quantities in each systems with L (= L x = L y ) saturate to the ground-state values at temperatures we have simulated. For example, the T -dependences of these quantities are not to be discernible at T = 0.05 and 0.01 for α = J ′ = 0.4 and 0.5, respectively. On the other hand, the L-dependence still remains in sizes we have calculated. The L-dependences of the ground-state spatial correlation lengths and the gap are shown in Fig. 4 for α = J ′ = 0.5, which is close to the critical point α c = J ′ c = 0.52337 (3). Their values in the thermodynamic limit are estimated by fitting ξ k (L) to ξ k (L) = ξ − b exp(−cL), where k = x, y, or τ , ξ is the value in the thermodynamic limit, and b and c are fitting parameters. As a result, we obtain ξ x = 3.0089(9), ξ y = 2.2097(6) and ∆ = 0.32261(4) for α = J ′ = 0.4 and ξ x = 11.998(9), ξ y = 9.312(10), and ∆ = 0.0913(2) for α = J ′ = 0.5. As α (= J ′ ) becomes smaller, i.e., the Before going into the QMC analysis, let us here summarize the ground-state phase diagram of the S = 1 system argued so far, which is shown in Fig. 5. For some points in the phase diagram the ground state is well understood by the previous theoretical and numerical studies: (1) (α, J ′ ) = (0, 0): The system consists of a set of the isolated antiferromagnetically-coupled spin pairs. The ground state is a trivial tensor product of dimer singlets sitting on each bond. (2) (α, J ′ ) = (1, 0): The system consists of isolated x-parallel Haldane chains. (3) (α, J ′ ) = (1, 1): The system is a uniform and isotropic 2D HAF. There exists an AF LRO in the ground state. 31 (4) J ′ = ∞: The system consists of y-parallel Haldane chains. Note that in this limit the value of α becomes irrelevant (see also discussions in Sec. V).
In their analysis by the cluster expansion method, Koga and Kawakami 13 derived three phases in the Q1D region, and they called the regions which includes point (1), (2) and (3)  To demonstrate our FSS analyses, let us begin with critical behavior near (1, J ′ c ), which separates the H-I and AF phases (J ′ c ≃ 0.04 due to Ref. 15). We sweep J ′ near supposed J ′ c with α fixed to unity. In Fig. 6, ξ x /L x , ξ y /L y , and ξ τ /β of the systems with L x = L y = β ≡ L (= 16, 24,· · ·, 64) are plotted. As one sees immediately, the data suffer from quite large corrections to scaling, i.e., the crossing point of the scaled correlation lengths with two different L's clearly shifts to larger J ′ as the system size increases. We attribute these large corrections to the strong spatial anisotropy in the coupling constants (J ′ ≪ 1). Indeed, the value of ξ y at J ′ = 0.0435 is 9.29(7) (see Fig. 6 (b)) even for the largest system size L = 64, which is quite smaller than those in the other directions (ξ x = 60.4(3) and ξ τ = 24.7(2)). Among the three correlation lengths, ξ x is the largest: the growth of the correlation is dominated only by the system size in the x-direction. This indicates that we need larger lattices, especially in the x-direction, in order to perform a precise FSS analysis.
To simulate larger lattices with minimal costs, we therefore optimize the aspect ratio, L x : L y : β, as explained below. Expecting that the scaled correlation lengths, ξ x /L x , ξ y /L y , and ξ τ /β, become nearly equal with each other at the critical point, we set the aspect ratio L x : L y : β as 7 : 1 : 2 based on the data presented in Fig. 6. With this ratio we simulate systems with L x = 168, 224, 280, and 336 and perform the FSS analyses. The raw data of ξ x /L x , ξ y /L y , and ξ τ /β with this aspect ratio are shown in Fig. 7. Now ξ y = 12.02(5) even for the smallest system size (L x = 168 and L y = 24) at J ′ = 0.0435. The ratios ξ x /L x , ξ y /L y and ξ τ /β in a common range of J ′ become nearly equal, and corrections to scaling become much smaller than in the previous ones as we expected. The FSS plot for ξ x is shown in Fig. 8. The resultant J ′ c and ν are as follows: (J ′ c , ν) = (0.043648(9), 0.69(1)) from ξ x , (0.043649(8), 0.71(1)) from ξ y , and (0.043648(7), 0.69(1)) from ξ τ . Averaging these three values we conclude with and ν = 0.70(1).
Fixing the value of J ′ c and ν thus determined, we next perform the FSS analysis on the staggered susceptibility χ s . The raw data of χ s vs. J ′ are shown in Fig. 9 and the fitting result is shown in Fig. 10. The latter yields The exponents ν and γ we have obtained for the S = 1 system again agree with those of the 3D classical Heisenberg model. 21 The critical point obtained just above is consistent with the previous result by the method involving the mean-field approximation, J ′ c ≥ 0.025, 12 and also with that by the recent QMC method, J ′ c = 0.040 (5). 15 On the other hand, the present result significantly differs from that obtained by the cluster-expansion method, J ′ c = 0.56 (1). 13 The present FSS analyses on our extensive QMC results make us possible to obtain the critical value with much higher accuracy than the other methods.

C. Ground-state phase diagram
In a similar way as described in the previous subsection, we obtain other critical points on the α-J ′ phase diagram as shown in Fig. 11. First, the scaled correlation lengths, ξ x /L x , ξ y /L y , and ξ τ /β, are calculated up to L x = 64 with either α or J ′ fixed. Sweeping J ′ or α with sufficiently high resolution, we regard a crossing point of these ξ's as the critical point. Note that the optimal aspect ratio depends strongly on the value of α and J ′ . However, we adopt L x : L y : β = 1 : 1 : 1 for simplicity. Although the results thus obtained suffer from relatively larger systematic corrections than those presented in the last subsection, the absolute magnitude of the systematic error in the estimates should be still smaller enough than the symbol size in Fig. 11.
For some critical points the FSS analysis as in the previous subsection is also carried out. We obtain the exponents ν and γ which are consistent with Eqs. (10) and (11), respectively. This supports that the quantum critical phenomena in the S = 1 system also belong to the same universality class as that of the 3D classical Heisen- berg model. An exception is the 1D critical point located at (α c , J ′ c ) = (0.5879(6), 0), 11 which separates the dimer phase from the Haldane phase. The apparent value of ν starts to deviate from (10) when α becomes closer to α c . This is attributed to the crossover to the critical phenomena belonging to the Gaussian universality class. 10 We confirm that the AF-LRO phase exists between the two spin-gap phases at least down to J ′ = 0.01 at α = α c . Although in the present simulation it is quite difficult to prove the existence of the AF-LRO phase at smaller J ′ , we believe that the the point (α c , 0) is tricritical: the 1D critical point is unstable against an infinitesimal interchain coupling and the AF LRO immediately appears as the same as in the S = 1/2 uniform chain. 12,[28][29][30] Interestingly, the H-II and D phases are adiabatically connected with each other. The gapless AF-LRO phase does not touch the line of α = 0 as seen in Fig. 12, where the part of the whole phase diagram (Fig. 11) near α = 0 is magnified. Indeed, on the α = 0 line, which corresponds to the S = 1 two-leg ladder, it is shown that there exists no critical point by the recent QMC study. 34 Thus, the D phase can be identified with the H-II phase, and there are only two distinct spin-gap phases, H-I and H-II in Fig. 11. The closeness of the critical line to the α = 0 line is due to the strong AF fluctuations, which already exist in the two-leg ladder system. 34

D. Correlation lengths and the gap
We obtain the explicit values of ξ x , ξ y , and ∆ in the ground state at (α, J ′ ) = (1, 0.04). They are calculated for systems with sizes L x = 168, 224, 280, and 336 and with the aspect ratio L x : L y = 7 : 1 at T regarded as zero temperature. Their T -dependences are negligible at T = 0.01. We extrapolate the finite-size data to the thermodynamic limit in the same way as explained for  the S = 1/2 system. We obtain ξ x = 39.2(1), ξ y = 5.67(1), and ∆ = 0.0632 (2). As J ′ becomes smaller, ξ x becomes smaller and ∆ larger to reach at J ′ = 0 the single chain values ξ x = 6.0153(3) and ∆ = 0.41048 (6). 20

V. DISCUSSIONS
The analyses presented in the preceding section revealed that the ground-state phase diagram of the S = 1 system has a rather complicated topology, i.e., the H-II and D phases are adiabatically connected with each other, though the channel between them is quite narrow (Fig. 11). On the other hand, as for the H-I and D phases, in the 1D system (J ′ = 0) these two spingapped phases are distinctively separated by the critical point at α c = 0.5879(6), 10,11 and they are distinguished by the string-order parameter, 23 which is zero in the former phase and finite in the latter one. The transition can be viewed as a rearrangement of dimer-singlet pattern between the (1,1)-and (2,0)-valence-bond-solid (VBS) states. 32,33 We emphasize that once J ′ is introduced, however, the string-order parameter should vanish even in the H-I phase, being similar to the S = 1 ladder. 34 Still one may consider that the two phases essentially differ with each other since they are separated by the AF-LRO phase. If, however, we introduce the bond alternation also in the y-direction, the two phases can be connected without passing the gapless state as explained below.
Let us consider the 2D HAF model defined in the extended parameter space: The original Hamiltonian (1) corresponds to the case with J x = 1, α x = α, J y = J ′ , and α y = 1. The Hamiltonian (12) is invariant under the exchange between (J x , α x ) and (J y , α y ).
To draw the phase diagram in this extended parameter space, it is convenient to introduce a parameter, R ≡ J y /(J x + J y ). The limits J y → 0 and J y → ∞ correspond to R = 0 and 1, respectively. Since α x (α y ) becomes irrelevant in the limit J x → 0 (J x → ∞), the whole phase diagram in the three-dimensional parameter space is shaped as a tetrahedron. In Fig. 13, we present the ground-state phase diagram parametrized by R, α x (1 − R), and α y R for S = 1/2 (a) and S = 1 (b). It should be noted that the phase diagram should be invariant under the transformation, (R, α x (1 − R), α y R) ↔ (1 − R, α y R, α x (1 − R)), reflecting the symmetry in the Hamiltonian explained above. In the phase diagram, the edge AB (CD) corresponds to isolated x-parallel (yparallel) decoupled chains, the edge AD isolated four-spin plaquettes, and the edge AC (BD) the two-leg ladders in y-direction (x-direction).
The face ABC (and also CDB) in Fig. 13 corresponds to the original phase diagram shown in Figs. 2 and 11, though the y-parallel-chain limit J ′ → ∞ in the original diagram is represented by one vertex C (B) in the new ones. In the extended phase diagram the shaded (unshaded) area represents the AF-LRO (spin-gapped) phase on the ABC-and CDB-faces.
It should be emphasized that on the ACD-and DBAfaces there is no AF-LRO phase, since the system is onedimensional dimerized two-leg ladder. There exist only the 1D critical points discussed already. Especially, there is no critical point on the edge AD. Therefore, in the S = 1 case, the three spin-gapped phases, H-I, D, and H-II, are connected by the path C→A→D→B. Similarly, in the S = 1/2 case, the two dimer phases, which correspond to the vertex A and D, respectively, are connected directly by the path A→D. Thus, in both cases, there are only two phases, namely, the spin-gapped phase and the AF-LRO one.

VI. CONCLUDING REMARKS
In this paper, we have investigated the ground-state phase diagram of S = 1/2 and S = 1 HAF on the anisotropic dimerized square lattice by means of the extensive QMC simulation with the continuous-imaginarytime loop algorithm and the FSS analyses. It is confirmed that, for both S =1/2 and 1, the quantum critical phenomena in the model belong to the same universality class as that of the 3D classical Heisenberg model, except for the 1D critical points, which belong to the same universality class as that of the Gaussian model. We have also demonstrated that the spin-gapped phases of the 1D chain are connected when we introduce the interchain couplings with bond alternation. In the 2D system, only one spin-gapped phase exists in both of the S = 1/2 and S = 1 systems.
The results obtained in the present work are considered to be the proper basis for investigation of peculiar phenomena observed in the Q1D HAF's mentioned at the beginning of this paper. We have already reported the QMC analysis on the site-dilution-induced AF LRO in these materials based on the Hamiltonian (1). 35 In order to discuss various experimental results quantitatively, it is certainly necessary to take into account other ingredients than in Eq. (1), such as the next-nearest-neighbor intrachain interaction and the single-ion anisotropy. They are beyond the scope of the present work. The present results, however, demonstrate the role of the higher dimensionality which has been overlooked so far.
Most of the numerical calculations in the present work have been performed on the DEC Alpha, SGI ORI-GIN 2000, SGI 2800, and RANDOM at the Materials Design and Characterization Laboratory, Institute for Solid State Physics, University of Tokyo and on the Hitachi SR-2201 at the Supercomputer Center, University of Tokyo. The program used in the present simulation was based on the library 'Looper version 2' developed by S.T. and K. Kato and also on the 'PARAPACK version 2' by S.T. The present work was supported by the "Research for the Future Program" (JSPS-RFTF97P01103) of the Japan Society for the Promotion of Science. S.T's work was partly supported by the Swiss National Science Foundation.