# Cluster algorithms for general-𝐒 quantum spin systems

## Abstract

We present a general strategy to extend quantum cluster algorithms for $`S=\frac{1}{2}`$ spin systems, such as the loop algorithm, to systems with arbitrary size of spins. In general, the partition function of a high-$`S`$ spin system is represented in terms of the path integral of a $`S=\frac{1}{2}`$ model with special boundary conditions in the imaginary-time direction. We introduce additional graphs to be assigned to the boundary part and give the labeling probability explicitly, which completes the algorithm together with an existing $`S=\frac{1}{2}`$ cluster algorithm. As a demonstration of the algorithm, we simulate the integer-spin antiferromagnetic Heisenberg chains. The magnitude of the first excitation gap is estimated as to be 0.41048(6), 0.08917(4), and 0.01002(3) for $`S=1`$, 2, and 3, respectively.

The world-line quantum Monte Carlo (QMC) method is one of the most powerful tools in numerical investigations of quantum spin systems . One of the main advantages of the QMC method over other numerical ones, such as exact diagonalization and density matrix renormalization group (DMRG) method, is that it is applicable to rather large systems in any dimensions and can estimate their physical quantities statistically exactly. However, the conventional algorithm, based on local updates of spin configurations (world lines), suffers from strong correlations between successive configurations at low temperatures or in the vicinity of a second-order phase transition. The diverging auto-correlation time virtually makes simulations slower and slower, and finally it becomes practically impossible to simulate larger systems at lower temperatures. This drawback is called critical slowing down.

Recently, the inventions of the loop algorithm  and of its continuous-time variant  have led a great improvement of the QMC techniques for the $`S=\frac{1}{2}`$ XXZ model . The loop algorithm, which is a kind of cluster algorithms, realizes updates of the configuration by flipping non-local objects, referred to as loops. It has been shown that it is fully ergodic and drastically reduces the auto-correlation time, often by orders of magnitude, especially at low temperatures. Furthermore, by using the continuous-time version of the algorithm, one can completely eliminate discretization error originating from the Suzuki-Trotter decomposition; simulations can be performed directly in the so-called Trotter limit.

In general, it is a highly non-trivial task to construct an efficient cluster algorithm for a given system, because symmetry or special properties of the target system should be taken into account explicitly in its construction. As for the spin systems, the development of cluster algorithms for higher-$`S`$ models remains as an important and challenging problem. A cluster algorithm for the general-$`S`$ XXZ model in the discrete-time formulation has been proposed by Kawashima and Gubernatis. Unfortunately, their algorithm is rather complicated (105 different graphs appear even in the $`S=1`$ case), and moreover the Trotter limit is not well-defined in the algorithm. Another kind of algorithm in the discrete-time representation has been used in Ref. .

More recently, Harada et al. proposed a continuous-time loop algorithm for the $`S=1`$ antiferromagnetic Heisenberg model , in which the $`S=1`$ system is mapped into a path integral of an $`S=\frac{1}{2}`$ system with special boundary conditions in the imaginary-time direction. In the present letter, we generalize their method to construct cluster algorithms for systems with arbitrary size of spins. For simplicity, we consider a spin-$`S`$ antiferromagnetic Heisenberg model on a bipartite lattice as an example. Generalization to other models is straightforward as discussed later.

We consider the Hamiltonian of $`N_\text{s}`$ spins and $`N_\text{b}`$ bonds, defined as follows:

$$=\underset{i,j}{\overset{N_\text{b}}{}}\left(S_i^xS_j^x+S_i^yS_j^yS_i^zS_j^z\right),$$
(1)   
where $`S_i^\alpha `$ ($`\alpha =x`$, $`y`$, $`z`$ and $`i=1`$, $`\mathrm{}`$, $`N_\text{s}`$) is the $`\alpha `$-component of the spin-$`S`$ operator at site $`i`$.

In order to construct a cluster algorithm for $`S>\frac{1}{2}`$, it is crucial to represent each spin operator $`S_i^\alpha `$ as a sum of $`S=\frac{1}{2}`$ Pauli operators $`\{\sigma _{i,\mu }^\alpha \}`$. Following Ref. , we substitute $`S_i^\alpha =\frac{1}{2}_{\mu =1}^{2S}\sigma _{i,\mu }^\alpha `$ into the Hamiltonian (1), which yields the following Hamiltonian of $`2SN_\text{s}`$ spins and $`4S^2N_\text{b}`$ bonds:

$$\stackrel{~}{}=\frac{1}{4}\underset{i,j}{\overset{N_\text{b}}{}}\underset{\mu ,\nu =1}{\overset{2S}{}}\left(\sigma _{i,\mu }^x\sigma _{j,\nu }^x+\sigma _{i,\mu }^y\sigma _{j,\nu }^y\sigma _{i,\mu }^z\sigma _{j,\nu }^z\right).$$
(2)   
We refer to the $`S=\frac{1}{2}`$ spins in Eq. (2) as subspins hereafter. In terms of the Hamiltonian (2) defined on the extended phase space, the partition function of the original Hamiltonian (1) can be expressed as

$$Z=\text{Tr}(\mathrm{exp}[\beta \stackrel{~}{}]P).$$
(3)   
Here, we introduce a projection operator $`P`$, which is a direct product of local symmetrization operators $`\left\{P_i\right\}`$ ($`i=1,\mathrm{},N_\mathrm{s}`$). Each $`P_i`$ acts on $`2^{2S}`$-dimensional Hilbert space spanned by $`\{\sigma _{i,\mu }^z\}`$ ($`\mu =1,\mathrm{},2S`$), and projects out unphysical states with $`𝐒_i^2<S(S+1)`$. Note that $`P`$ commutes with $`\stackrel{~}{}`$, because $`\stackrel{~}{}`$ is invariant under the exchange of subspin indices at each site by definition.

Applying a Suzuki-Trotter decomposition for the exponential operator in Eq. (3) and inserting complete sets of eigenstates of $`\{\sigma _{i,\mu }^z\}`$ between the exponential factors, we obtain the following path-integral representation of the partition function:

$$Z=\underset{𝒞}{}W(𝒞)P(𝒞),$$
(4)   
where $`𝒞=\{𝒞_{i,\mu ,\tau }\}`$ ($`i=1,\mathrm{},N_\mathrm{s}`$, $`\mu =1,\mathrm{},2S`$, and $`0\tau \beta `$) denotes a world-line configuration of subspins and $`𝒞=\{𝒞_{i,\mu ,0},𝒞_{i,\mu ,\beta }\}`$ ($`i=1,\mathrm{},N_\mathrm{s}`$ and $`\mu =1,\mathrm{},2S`$) denotes a configuration of subspins at the boundaries in the imaginary-time direction. We take the continuous-imaginary-time representation, and $`𝒞_{i,\mu ,\tau }`$ denotes the $`\mu `$-th subspin direction (+1 or - 1) at $`i`$-th site and imaginary time $`\tau `$. The weight $`P(𝒞)=_iP_i(𝒞_i)`$ originates from the projection operator $`P`$, which can be interpreted as soft boundary conditions in the imaginary-time direction. Each $`(4S)`$-body local boundary weight $`P_i(𝒞_i)`$ takes a value of the inverse of the number of different configurations which are connected with $`𝒞_i`$ by permutation operations, i.e.,

$`P_i(𝒞_i)=\{\begin{array}{cc}{\displaystyle \frac{(2Sn_i)!n_i!}{(2S)!}}\hfill & \text{if }n_{i,0}=n_{i,\beta }\text{ (}n_i\text{)}\hfill \\ 0\hfill & \text{otherwise,}\hfill \end{array}`$ (7)   
where $`n_{i,\tau }=_\mu 𝒞_{i,\mu ,\tau }`$.

Note that apart from the boundary conditions, the weight $`W(𝒞)`$ is completely equivalent to that appears in the path integral representation of the system described by the Hamiltonian (2). Therefore, for that part, we adopt the same labeling rule as the original $`S=\frac{1}{2}`$ continuous-time loop algorithm , which assigns a ‘graph’ $`𝒢_W`$ to a configuration $`𝒞`$ with labeling probability

$$T_W(𝒢_W|𝒞)=\frac{V_W(𝒢_W)\mathrm{\Delta }_W(𝒞,𝒢_W)}{W(𝒞)}.$$
(8)   
Here, we follow the general framework of cluster algorithms presented in Ref. . The weight $`V_W(𝒢_W)`$ is $`𝒞`$-independent and non-negative, and $`\mathrm{\Delta }_W(𝒞,𝒢_W)`$ is a compatibility function, which takes 0 or 1. They satisfy

$$W(𝒞)=\underset{𝒢_W}{}V_W(𝒢_W)\mathrm{\Delta }_W(𝒞,𝒢_W)$$
(9)   
for any $`𝒞`$. By this procedure the $`(d+1)`$-dimensional space is decomposed into a set of loops. Note that at this stage, some of the loops remain opened, because we have not defined any graphs for the boundary part.

Next, for each local boundary weight $`P_i(𝒞_i)`$, we introduce $`(2S)!`$ types of graphs $`\{𝒢_{Pi}\}`$ ($`i=1,\mathrm{},N_\mathrm{s}`$), each of which consists of $`2S`$ edges connecting one of subspins at $`\tau =0`$ to a subspin at $`\tau =\beta `$ one by one (Fig. 1). We define a compatibility function $`\mathrm{\Delta }_P(𝒞_i,𝒢_{Pi})`$ as it takes 1 if every edge connects two subspins which have an identical direction, or it takes 0 otherwise. By using this compatibility function, we can decompose the local boundary weight as

$$P_i(𝒞_i)=\underset{𝒢_{Pi}}{}\mathrm{\Delta }_P(𝒞_i,𝒢_{Pi})/(2S)!,$$
(10)   
because the number of compatible graphs for a configuration $`𝒞_i`$ is given by $`(2Sn_i)!n_i!`$. We take the labeling probability for the local boundary weight as

$$T_P(𝒢_{Pi}|𝒞_i)=\frac{\mathrm{\Delta }_P(𝒞_i,𝒢_{Pi})/(2S)!}{P_i(𝒞_i)}=\frac{\mathrm{\Delta }_P(𝒞_i,𝒢_{Pi})}{(2Sn_i)!n_i!},$$
(11)   
that is, a graph is assigned out of the graphs compatible with $`𝒞_i`$ with equal probability. These boundary graphs make the remaining opened loops to be closed. By choosing the flipping probability to be free, that is, flipping spins on each loop simultaneously with probability $`\frac{1}{2}`$, one can show that the present stochastic process satisfies the detailed-balance condition .

The present algorithm includes the loop algorithm for the $`S=1`$ antiferromagnetic Heisenberg model by Harada et al  as a special case. However, as already seen, the present strategy does not depend so much on details of the model one considers; one can easily construct a cluster algorithm, if the mapped $`S=\frac{1}{2}`$ model has a cluster algorithm, which covers a model with general XYZ interaction or single-ion anisotropy , and the transverse-field Ising model . It can be applied also for a system with random size of spins and even for the classical Ising model . The resulting general-$`S`$ algorithm is ergodic, if the $`S=\frac{1}{2}`$ cluster algorithm lying at the base is ergodic. The details of the algorithm for these models and the proof of its ergodicity will be presented elsewhere . Note that the number of graphs introduced for the local boundary weight increases quite rapidly as $`S`$ increases. However, the computational time in selecting a graph to be assigned is merely proportional to $`S`$, because the procedure is nothing but the random-permutation generation.

As an application of the algorithm, we simulated the antiferromagnetic Heisenberg chains with $`S=1`$, 2, and 3. It is conjectured by Haldane  that the antiferromagnetic Heisenberg chain of integer spins has a finite excitation gap $`\mathrm{\Delta }(S)`$ above its unique ground state, and the antiferromagnetic spin correlation along the chain decays exponentially with a finite correlation length $`\xi _x(S)`$. For $`S=1`$ and 2, a number of numerical studies have been accomplished (e.g., see ) to confirm the validity of Haldane’s conjecture. However, estimation of the first excitation gap of higher-spin chains has not yet been successful, since the magnitude of $`\mathrm{\Delta }(S)`$ is considered to become exponentially small as $`S`$ increases .

Consider a spin-$`S`$ chain of $`L`$ sites at temperature $`T`$ ($`=1/\beta `$). We assume $`L`$ is even. The correlation function of the staggered magnetization in the imaginary-time direction:

$$C(\tau ;L,\beta )=\frac{1}{L^2\beta }\underset{i,j=1}{\overset{L}{}}_0^\beta 𝑑t(1)^{|ij|}S_i^z(t)S_j^z(t+\tau )$$
(12)   
is an even function of $`\tau `$, and satisfies $`C(\tau +\beta ;L,\beta )=C(\tau ;L,\beta )`$. At sufficiently low temperatures, the correlation function is expressed well as a sum of exponential functions, that is,

$$C(\tau ;L,\beta )=\underset{i=0}{}c_i\mathrm{cosh}\left[\frac{\tau \beta /2}{\xi _{\tau ,i}(L)}\right]\text{for }\beta \xi _{\tau ,0}\text{.}$$
(13)   
Here, we assume $`\xi _{\tau ,0}>\xi _{\tau ,1}>\xi _{\tau ,2}>\mathrm{}`$ without loss of generality. The coefficients $`\left\{c_i\right\}`$ are directly related to the dynamic structure factor at momentum $`k=\pi `$. In terms of $`\xi _{\tau ,0}(L)`$, the gap to be estimated is given by

$$\mathrm{\Delta }=\underset{L\mathrm{}}{lim}\frac{1}{\xi _{\tau ,0}(L)}.$$
(14)   
In general, to solve Eq. (13) directly is extremely ill-posed . However, as shown below, we can construct a systematic series of estimators at least for $`\xi _{\tau ,0}(\tau )`$, if the coefficient $`c_i`$ in Eq. (13) converges to zero rapidly enough for large $`i`$, and also if the difference between $`\xi _{\tau ,0}(L)`$ and $`\xi _{\tau ,1}(L)`$ remains finite even in the thermodynamic limit. For S=1, these two conditions are numerically shown to be satisfied . We expect similar situations for higher $`S`$, although there is no exact argument.

For a given $`L`$, the well-known second-moment estimator  of the correlation length,

$$\widehat{\xi }_\tau ^{(2)}=\frac{\beta }{2\pi }\sqrt{\frac{\stackrel{~}{C}(0)}{\stackrel{~}{C}(2\pi /\beta )}1},$$
(15)   
converges to $`\xi _{\tau ,0}(L)`$ in the low-temperature limit, besides systematic corrections of $`O(a_i\xi _{\tau ,i}/\xi _{\tau ,0})`$ ($`i=1,2,\mathrm{}`$). Here, $`\stackrel{~}{C}(\omega )`$ is the Fourier transform of the imaginary-time spin correlation function, i.e., $`\stackrel{~}{C}(\omega )=_0^\beta C(\tau )e^{i\omega \tau }𝑑\tau `$. In the present algorithm, $`\stackrel{~}{C}(\omega )`$ can be measured directly by means of the improved estimator, which reduces the variance of the data greatly. We also consider a fourth-moment estimator:

$$\widehat{\xi }_\tau ^{(4)}=\frac{\beta }{4\pi }\sqrt{3\frac{\stackrel{~}{C}(0)\stackrel{~}{C}(2\pi /\beta )}{\stackrel{~}{C}(2\pi /\beta )\stackrel{~}{C}(4\pi /\beta )}1},$$
(16)   
which has smaller corrections of $`O(a_i(\xi _{\tau ,i}/\xi _{\tau ,0})^3)`$. Construction of higher-order estimators is possible in a straightforward way.

The temperature dependence of $`\widehat{\xi }_\tau (L,T)`$ and also of the correlation length along the chain, $`\widehat{\xi }_x(L,T)`$, is shown in Fig. 2. In the present simulation, the system size is taken as $`L=2S/T`$. The integrated auto-correlation time of $`\stackrel{~}{C}(0)`$ is of order unity and no significant sign of its growth is observed for larger $`L`$ and $`1/T`$. Measurement of physical quantities is performed after discarding first $`10^3`$ Monte Carlo steps (MCS) for thermalization. One Monte Carlo step of the $`S=3`$ chain with $`L=5792`$ and $`T=0.0010359`$, which is mapped to the system of 34752 subspins and 208512 bonds prior to the simulation, takes about 16 seconds on 256 nodes of Hitachi SR-2201.

As seen clearly in Fig 2, $`\widehat{\xi }_\tau (L,T)`$ and $`\widehat{\xi }_x(L,T)`$ converge quite rapidly (probably exponentially) to a finite value for small $`T`$. We observe that for each $`S`$ the data which satisfy the conditions, $`1/T>6\widehat{\xi }_\tau `$ and $`L>6\widehat{\xi }_x`$, exhibit no temperature (and system-size) dependence besides statistical fluctuations. The difference between the second-moment estimate (Eq. (15)) and the fourth-moment one (Eq. (16)) is invisible in the vertical scale of Fig 2. The difference between these two estimates at the lowest temperature is about 0.2% and 0.1% for $`S=1`$ and $`S=2`$, respectively. For $`S=3`$, both coincide within the statistical errors. Furthermore, it is confirmed that there is no significant difference between the fourth-moment estimate and the sixth-moment one even for $`S=1`$.

Thus, by using the fourth-moment estimator, we conclude

$$\mathrm{\Delta }(S)=\{\begin{array}{cc}0.41048(6)\hfill & \text{ for }S=1\hfill \\ 0.08917(4)\hfill & \text{ for }S=2\hfill \\ 0.01002(3)\hfill & \text{ for }S=3\hfill \end{array}$$
(17)   
as the magnitude of the Haldane gap. The results for other physical quantities, such as the energy density and the staggered susceptibility, are presented in Table I. It should be emphasized that the present results are obtained without any extrapolation procedure; they are simply obtained by a single Monte Carlo run on the largest system at the lowest temperature for each $`S`$.

For $`S=1`$, the present estimate is completely consistent with $`0.41050(2)`$ and $`0.41049(2)`$ obtained by the DMRG calculation  and by the exact diagonalization , respectively. For $`S=2`$, on the other hand, the numerical uncertainty in the present estimate is much smaller than in the previous studies (see, e.g., TABLE I in Ref. ). Furthermore, our estimate is slightly larger than $`0.0876(13)`$, which is obtained by the most recent DMRG calculation . In the DMRG study, for some technical reasons, open boundary conditions have been used, which is known to give quite large systematic corrections compared to the periodic boundary conditions. The reason of the disagreement might be due to an inappropriate scaling assumption in the DMRG study. Finally, as for the $`S=3`$ case, it might be practically impossible to estimate the value of the gap by other numerical methods. The present result is completely new to our best knowledge.

The present authors would like to thank H. Takayama, N. Kawashima, H. G. Evertz, and K. Hukushima for stimulating discussions and comments. Most of numerical calculations for the present work have been performed on the CP-PACS at University of Tsukuba, Hitachi SR-2201 at Supercomputer Center, University of Tokyo, and the RANDOM at MDCL, Institute for Solid State Physics, University of Tokyo. The present work is supported by the “Large-scale Numerical Simulation Program” of Center for Computational Physics, University of Tsukuba, and also by the “Research for the Future Program” (JSPS-RFTF97P01103) of Japan Society for the Promotion of Science.