Haldane Gaps of Large-S Heisenberg Antiferromagnetic Chains and Asymptotic Behavior

The one-dimensional Heisenberg antiferromagnets of large-integer-$S$ spins are studied; their Haldane gaps are estimated by the numerical diagonalization method for $S=5$ and $6$. We successfully obtain a monotonically increasing sequence of finite-size energy difference data corresponding to the Haldane gaps from the huge-scale parallel calculations of diagonalization under the twisted boundary condition and create a monotonically decreasing sequence within the range of system sizes treated in this study from the monotonically increasing sequence. Consequently, the gaps for $S=5$ and $6$ are estimated to be $0.000050 \pm 0.000005$ and $0.0000030 \pm 0.0000005$, respectively. The asymptotic formula of the Haldane gap for $S\rightarrow\infty$ is examined from the new estimates to determine the coefficient in the formula more precisely.


Introduction
The Haldane gap -the energy gap between the unique ground state and the first excited state for the integer-S Heisenberg antiferromagnets in one dimension -is now well known; however, the presence of the gap surprised many condensed-matter physicists when it was originally conjectured in Refs. 1 and 2 by mapping the Heisenberg chain to the nonlinear σ model. Various investigations concerning whether the gap is present or absent have been carried out since Refs. 1 and 2 were published. Currently, the existence of a nonzero gap is widely believed. Extensive studies concerning this phenomenon have contributed considerably to our understanding of the various properties of quantum spin systems.
To date, various approaches have been attempted to estimate the magnitude of the gap. In particular, the Haldane gap for S = 1 was studied 3,4) soon after the conjecture because the S = 1 gap is larger than those for larger S . It is now known that the density matrix renormalization group (DMRG), 5) quantum Monte Carlo (QMC), 6) and numerical diagonalization (ND) 7) calculations give estimates that agree with each other to five decimal places: ∆/J ∼ 0.41048, where ∆ and J represent the gap value and the strength of the interaction defined later, respectively. For S = 2, after various studies, 6,[8][9][10] the gap estimates from the three approaches agree with each other within errors; ∆/J ∼ 0.089.
For cases of even larger S , however, it is more difficult to estimate the gap values. Therefore, the number of studies for S larger than 2 is much smaller. The first report of the S = 3 case is Ref. 6, which reported ∆/J = 0.01002 ± 0.00003 from a QMC calculation. This estimate was confirmed by an ND study 7) reporting ∆/J = 0.0092 ± 0.0010. The ND study also provided the gap estimates for even larger S : ∆/J = 0.00072 ± 0.00009 for S = 4 and ∆/J = 0.000047 ± 0.000010 for S = 5. This ND study was the first report for S = 4 and S = 5. Among the estimates, a very recent QMC calculation 11) gave ∆/J = 0.000799 ± 0.000005 for S = 4, which agrees with the estimate by the ND calculation. Un- * hnakano@sci.u-hyogo.ac.jp fortunately, no other approaches have successfully estimated the gap value for S larger than 4 to the best of our knowledge. The estimation of the Haldane gap for large S is still one of the most challenging issues in condensed-matter physics, particularly, from the viewpoint of computational statistical physics.
One reason why it is quite difficult to estimate the gap values numerically is that a high-cost calculation is necessary to treat considerably large systems. The requirement is strongly related to the fact that the gap values for large S are extremely small although they are nonzero. It is notable that the ND study 7) succeeded in the estimation by using the twisted boundary condition even though the ND method can only treat clusters that are much smaller than those treated by the QMC and DMRG calculations. Therefore, the validity of the method used in Ref. 7 should be examined from various viewpoints.
Concerning the Haldane gaps for large S , the original studies by Haldane 1,2) derived their asymptotic formula as follows: for S → ∞, where |S| represents the amplitude of each spin included in the target Heisenberg chain. However, the coefficient β cannot be determined only by Haldane's argument.
To determine β, numerical approaches are required for sufficiently large S . The ND calculations in Ref. 7 gave an estimate of β = 12.8 ± 1.5.
No other estimates are known to the best of our knowledge. If one combines Eq. (1) and the estimate of β, one can infer the gap value for even larger S . For S = 6, for example, one can easily predict from Eq. (1). If one can directly estimate the gap value for S = 6, it will be possible to compare the estimate with Eq.  Under these circumstances, the purpose of this study is to verify the method used in Ref. 7 from the following two aspects. The first one is to obtain a direct estimate of the gap value for S = 6 by the method used in Ref. 7. The additional estimate allows it to be compared with the prediction (3). The second aspect is to obtain finite-size energy gaps for S = 5 with clusters that were not treated in Ref. 7. Although in Ref. 7, the authors were able to treat only clusters up to 10 sites, in this study, we additionally report results for clusters with 12 and 14 sites; we will confirm that the additional results show the common behavior of the results up to 10 sites and give an estimate that is more precise than that in Ref. 7.
This paper is organized as follows. In the next section, the model Hamiltonian will be introduced. Our numerical method will also be explained. The third section is devoted to the presentation and discussion of our results. We will first treat the case of S = 5. Next, we will study the case of S = 6. From the additional gap values for S = 5 and S = 6, finally, we will estimate the coefficient β more precisely. In the final section, we will summarize our results and give some remarks.

Model Hamiltonian and Numerical Method
The Hamiltonian studied here is given by where S i represents the spin-S spin operator at site i. In this study, we particularly focus our attention on the cases of S = 5 and S = 6. We consider the case of an isotropic interaction in spin space in this study. The label of a spin site is represented by i, which should be an integer. The number of spin sites is denoted by N, which is assumed to be an even integer. Energies are measured in units of J; hereafter, we set J = 1, which indicates that the system is an antiferromagnet. We treat finite-size clusters with system size N under the twisted boundary condition. This condition is given by which should be noted by the difference from the periodic boundary condition given by S N+1 = S 1 . Owing to the twisted boundary condition, the system size N should satisfy N ≥ 4. Note also that the twisted boundary condition in studies of quantum spin systems is not so strange because the condition is effectively used in studies using level-spectroscopy analysis. [12][13][14][15] In particular, in Ref. 13, the authors used the twisted boundary condition to study the S = 1 one-dimensional antiferromagnet with bond alternation and successfully determined the gapless point between the two gapped phases, namely, the Haldane phase and the dimer phase. To find the We determine a normalized system sizeÑ defined as N + N 0 so that the three data for N = 4, 6, and 8 reveal a linear dependence: N 0 = 6.25840221. Blue closed diamonds represent C(N). Equation (7) gives C(N) from the finite-size gaps of system sizes N − 2, N, and N + 2.
gapless point in finite-size systems, in this reference, the authors searched for a level-crossing point of the ground state of the system under the twisted boundary condition. This means that a gapless case is realized in finite-size systems as a situation of the doubly degenerate ground state. The degeneracy indicates that the energy difference between the two states vanishes. One thus finds that the twisted boundary condition can appropriately capture such a gapless case. In Ref. 7, on the other hand, the authors applied the twisted boundary condition in gapped cases to show that the condition can contribute to the gap estimation. The standing position of this paper is to confirm the validity of the method used in Ref. 7. The reason why the twisted boundary condition is used here will also be mentioned when our practical results are presented in the next section.
We carry out our numerical diagonalizations based on the Lanczos algorithm to obtain the lowest energies of H in the subspace characterized by j S z j = M. Note here that the zaxis is taken as the quantized axis of each spin. ND calculations can treat quantum effects without any approximations and provide us with very precise results within numerical errors. Thus, one can obtain reliable information about the system. We focus our attention on the lowest energy excitation; therefore, our calculations are carried out for M = 0. The lowest energy and the first excited energy for a given N are denoted by E 0 and E 1 , respectively. We evaluate the energy difference given by for each N. Some of the Lanczos diagonalizations were carried out using an MPI-parallelized code that was originally developed in Ref. 7. The usefulness of our program was confirmed in large-scale parallelized calculations. 10,[16][17][18][19][20][21][22][23][24][25][26] Note here that the largest scale calculations in this study have been carried out using either the K computer or Oakforest-PACS. We carry out our calculations up to N = 14 for S = 5 and up to N = 12 for S = 6. Among our calculations, the case of N = 14 for S = 5 is the largest 27) with respect to the dimension of the

Results and Discussion
3.1 Case for S = 5 Now, we study our results for S = 5; our numerical results under the twisted boundary condition are presented in Table I. One can observe that the finite-size gap ∆ N is monotonically increasing with respect to N. This monotonic increase is a decisive merit of using the twisted boundary condition. As mentioned in Ref. 7, if we use the periodic boundary condition, each finite-size gap has a significantly large magnitude and is monotonically decreasing with increasing N. From such a data sequence, it is extremely difficult to extrapolate the sequence to the thermodynamic limit and to find whether the nonzero gap is present or absent. On the other hand, a data sequence with a monotonic increase can easily enable us to find that the gap opens. In Ref. 7, the dependence for S = 5 was observed up to 10 sites. In this study, we successfully clarified that our additional results for N = 12 and 14 maintain the monotonic increase.
Next, we analyze our results and produce another sequence that is monotonically decreasing within the range of N treated here and that approaches the gap in the thermodynamic limit N → ∞ according to the method used in Ref. 7, which employs a two-step procedure.
The first step is to draw a plot of ∆ N as a function of 1/Ñ instead of raw N, where we introduce a renormalized system sizeÑ defined as N+N 0 so that the three initial data for N = 4, 6, and 8 reveal a linear dependence in the plot of ∆ N . The result for S = 5 is depicted in Fig. 1. If we skip the first step and draw a usual plot of ∆ N versus 1/N, the 1/N dependence of ∆ N shows concave upward behavior for small sizes. In Ref. 7, the authors reported that the data up to N = 10, which was the maximum of N in this reference, do not show the convex upward behavior in the plot of ∆ N versus 1/N. Under such a situation, it was difficult to obtain an appropriate extrapolated value for ∆ N from the analysis using the plot of ∆ N versus 1/N, which was why ∆ N versus 1/Ñ was plotted in Ref. 7. In this paper, on the other hand, we additionally report data for N = 12 and 14, in which convex upward behavior is observed for large system sizes in the plot of ∆ N versus 1/N. Since the purpose of this paper is to examine the validity of the method used in Ref. 7, we followed the same procedure and will later discuss the cases of usingÑ or raw N.
The second step is to create a decreasing sequence from ∆ N , which is an increasing sequence in the plot of ∆ N as a function ofÑ. We focus our attention on three neighboring We determine a normalized system sizeÑ defined as N + N 0 so that the three data for N = 4, 6, and 8 reveal a linear dependence: N 0 = 13.0448842. Blue closed diamonds represent C(N). Equation (7) gives C(N) from the finite-size gaps of system sizes N − 2, N, and N + 2. data points of system sizes N−2, N, and N+2 in Fig. 1. When we apply the fitting curve of to the neighboring three data points, we can determine the parameters C, D, and E uniquely for a given N. Thus, we use C(N), D(N), and E(N) hereafter. Note here that E(N = 6) is necessarily unity owing to the above first step. The result for C(N) is also depicted at the correspondingÑ in Fig. 1. One can observe in Fig. 1 that C(N) is monotonically decreasing when N is increased within the range of N treated here. The sequence C(N) is approaching the gap value in the thermodynamic limit from the side that is opposite to that of ∆ N . Note that C(N) approaches the gap value in the thermodynamic limit regardless of whether one usesÑ or N. The behavior of approaching the gap value from the opposite side suggests that C(N) for the largest system size is the closest to the gap value that we want to know finally; C(N) for the largest system size in the present paper is ∼0.0000542. Therefore, we obtain as a new estimate whose error is smaller than that in Ref. 7. Consequently, we find that the method used in Ref. 7 can be applied for the additional data for N = 12 and N = 14, which are larger than those treated in Ref. 7. Note here that C(N) depends on whether one uses the plot of ∆ N versus 1/Ñ or that versus 1/N. If one uses the plot of ∆ N versus 1/N, one obtains C(N) for the largest system size in this paper to be ∼0.0000595. The comparison of the two results for C(N) from N or N suggests that the usage ofÑ causes the error of an estimate for the gap value to become smaller.

Case for S = 6
Next, we study the case of S = 6; our numerical results under the twisted boundary condition are presented in Table II. This study is the first report giving finite-size gaps for S = 6 to the best of our knowledge. One can also observe that the finite-size gap ∆ N is monotonically increasing with respect to N and find that the behavior of the monotonic increase is certainly maintained for S = 6. Let us analyze the monotonically increasing ∆ N and obtain a monotonically decreasing sequence that is approaching the gap value in the thermodynamic limit according to the same method as in the S = 5 case; the results are depicted in Fig. 2. We successfully obtain a monotonically decreasing C(N) within the range of N treated here. We finally obtain an estimate for the gap value for S = 6 to be ∆(S = 6)/J = 0.0000030 ± 0.0000005. (9) One finds that the present Eq. (9) agrees with the prediction (3). The agreement strongly indicates the validity of the method used in Ref. 7 and this study. Our estimate for S = 6 will be investigated in the future if other methods become available.

Asymptotic behavior
Now, we examine the asymptotic formula of Eq. (1) for the Haldane gap for S → ∞ from our estimate of the S = 5 and 6 gaps. To do this, we introduce new parameters x = S −1 and y = S −1 log(S 2 J/∆(S )) when the amplitude of each spin is taken to be |S| = S in this analysis. The asymptotic formula (1) is rewritten as Let us input the obtained estimates of the Haldane gaps to ∆(S ) in y and plot the x dependence of y. The result is depicted in Fig. 3. One can find linear behavior for finite but large S up to S = 6. We have fitted our data for S = 5 and 6 using the straight line (10); the best fit is produced by β = 13.0 ± 1.2. The good linear behavior suggests that the asymptotic formula (1) holds well for large S .

Summary and Remarks
We have studied the Haldane gaps of the integer-S Heisenberg antiferromagnetic chain model when S is large by the Lanczos diagonalization method. We have successfully obtained an estimate for S = 5 that is more precise than that in the previous study. We have also succeeded in obtaining an estimate for S = 6 for the first time. This study allows us to confirm the validity of numerical diagonalization calculations under the twisted boundary condition and the analysis giving an estimate for the gap value in the thermodynamic limit. Cases for even larger S will be studied in future works. Such studies are expected to contribute much to our fundamental understanding of quantum magnetism and further development of our numerical techniques.
Finally, let us give some comments concerning experiments. To date, experimental studies of the Haldane gaps have also been carried out. For S = 1, NEMP [28][29][30] is a famous achievement as a good candidate material. In recent years, some materials have been reported for S = 2. 31, 32) Unfortunately, there are no experimental reports for S = 3 to the best of our knowledge. Since the gap magnitudes for S larger than 2 are quite small, it is necessary to remove effects from interactions other than the Heisenberg chain interaction as much as possible or to separate the intrinsic gap in some skillful way. Even in the presence of some difficulties, experimental attempts for large-S Haldane gaps should be tackled in the future.