arXiv Vanity renders academic papers from arXiv as responsive web pages so you don’t have to squint at a PDF. Read this paper on arXiv.org.

# Large-scale Monte Carlo simulations of the isotropic three-dimensional Heisenberg spin glass

L. W. Lee Department of Physics, University of California, Santa Cruz, California 95064 Dept. of Mechanical Engineering, University of Washington, Seattle, WA 98195 Present address    A. P. Young Department of Physics, University of California, Santa Cruz, California 95064
May 10, 2020
###### Abstract

We study the Heisenberg spin glass by large-scale Monte Carlo simulations for sizes up to , down to temperatures below the transition temperature claimed in earlier work. The data for the larger sizes show more marginal behavior than that for the smaller sizes, indicating the lower critical dimension is close to, and possibly equal to three. We find that the spins and chiralities behave in a quite similar manner.

###### pacs:
75.50.Lk, 75.40.Mg, 05.50.+q

## I Introduction

Following the convincing numerical work of Ballesteros et al.Ballesteros et al. (2000) there has been little doubt that Ising spin glasses in three dimensions have a finite temperature transition. In this paper we shall study a related model for which the existence of a finite temperature transition is still controversial: the isotropic Heisenberg spin glass, which is composed of classical spins with three components. Early work McMillan (1985); Morris et al. (1986); Olive et al. (1986) indicated a zero temperature transition, or possibly a transition at a very low but non-zero temperature. However, following the pioneering work of VillainVillain (1977), which emphasized the role of “chiralities” (Ising-like variables which describe the handedness of the non-collinear spin structures), Kawamura and TanemuraKawamura and Tanemura (1987) proposed, for the XY case (which have two-component spins), that the spin glass transition only occurs at and that a chiral glass transition occurs at a finite temperature . This scenario requires that spins and chiralities decouple at long length scales. Kawamura and collaborators subsequently proposedKawamura (1998); Hukushima and Kawamura (2000, 2005) that this “spin-chirality decoupling” scenario also holds for Heisenberg spin glasses.

However, the absence of a spin glass transition in Heisenberg spin glasses has been challenged by Matsubara et al.Matsubara et al. (2001); Endoh et al. (2001), and Nakamura and EndohNakamura and Endoh (2002) who argued that the spins and chiralities order at the same low but finite temperature. Recently Picco and RitortPicco and Ritort (2005) also claimed a finite and inferred that probably , though they did not investigate the chiralities directly. In earlier workLee and Young (2003a), referred to as LY, we studied spin and chiral correlations on an equal footing, using the method of analysis that was the most successful for the Ising spin glassBallesteros et al. (2000), namely finite-size scaling of the correlation length. Considering a modest range of sizes, , we found that the behavior of spins and chiralities was quite similar and they both had a finite temperature transition, apparently at the same temperature.

However, quite recently Campos et al.Campos et al. (2006) were able to study larger sizes than LY, up to . They agreed with LY that there is a single transition at which both spins and chiralities order, but they also found evidence for crossover, at the largest sizes, to a “marginal” behavior, reminiscent of that at the Kosterlitz-Thouless-BerezenskiiKosterlitz and Thouless (1973); Berezinskii (1970) (KTB) transition in the two-dimensional XY ferromagnet where there is a finite transition temperature but no long-range order for . In fact the region below is a line of critical points in the KTB theory. If a line of critical points also exists in the three-dimensional Heisenberg spin glass, then is the lower critical dimension , below which there is no transition. However, given numerical uncertainties, Campos et al. cannot rule out the possibility that is slightly less than three, in which case there is spin glass order below .

Hukushima and KawamuraHukushima and Kawamura (2005) studied sizes up to and found more marginal behavior when comparing and 20, than for the smaller sizes. However, they argued that this effect is greater for the spin correlation length than for the chiral correlation length, and hence concluded that, while is finite, is zero or possibly non-zero but less than , i.e. there is spin-chirality decoupling.

In this paper we perform Monte Carlo simulations of the Heisenberg spin glass, along the lines of LY, but for larger sizes. Our main motivation is to investigate the claim of Campos et al.Campos et al. (2006) that there is a line of critical points for . In order to test whether there is a critical line, as proposed by Campos et al, or the usual situation of a single critical point at , it is necessary to investigate the behavior of the system below the estimated . Campos et al. were not able to do this and the evidence for the critical line was based on estimating corrections to scaling at . This is rather indirect, and perhaps not very reliable because the range of sizes and quality of the data, are not sufficient to disentangle the various corrections to scaling unambiguously. Very recently, the analysis of Campos et al. has been criticized by Campbell and KawamuraCampbell and Kawamura (2007). Here, we do not rely on corrections to scaling but obtain data below and so directly find evidence for more marginal behavior at larger sizes.

As we shall see, the transition temperature in the Heisenberg spin glass is very low and so it is surprising to us that Campos et al. Campos et al. (2006) did not use the technique of “parallel tempering”Hukushima and Nemoto (1996); Marinari (1998), which is the commonly used approach to speed up simulations of spin glasses, especially at very low temperatures. Equilibrating lattices as large as is very challenging, especially in the absence of parallel tempering. A second motivation of our study is therefore to use parallel tempering and to combine this with a very useful equilibration test described in the next section, to ensure that the data are equilibrated.

From our results, we see a strong crossover to much more marginal behavior for sizes , in agreement with Campos et al.Campos et al. (2006). Whether there is a KTB-like critical line as proposed by Campos et al., is, however, unclear. The lower critical dimension, seems to be close to, or possibly equal to, three. The behavior of the spins and chiralities is rather similar, and so, in contrast to Kawamura and collaboratorsHukushima and Kawamura (2005); Campbell and Kawamura (2007), we do not feel that our data supports spin-chirality decoupling.

## Ii Model and analysis

We use the standard Edwards-Anderson spin glass model

 H=−∑⟨i,j⟩JijSi⋅Sj, (1)

where the are 3-component classical vectors of unit length at the sites of a simple cubic lattice, and the are nearest neighbor interactions with a Gaussian distribution with zero mean and standard deviation unity. Periodic boundary conditions are applied on lattices with spins.

The spin glass order parameter, , at wave vector , is defined to be

 qμν(k)=1N∑iSμ(1)iSν(2)ieik⋅Ri, (2)

where and are spin components, and “” and “” denote two identical copies of the system with the same interactions. From this we determine the wave vector dependent spin glass susceptibility by

 χSG(k)=N∑μ,ν[⟨|qμν(k)|2⟩]av, (3)

where denotes a thermal average and denotes an average over disorder. The spin glass correlation length is then determinedPalassini and Caracciolo (1999); Ballesteros et al. (2000), from

 ξL=12sin(kmin/2)(χSG(0)χSG(kmin)−1)1/2, (4)

where .

For the Heisenberg spin glass, KawamuraKawamura (1998) defines the local chirality in terms of three spins on a line as follows:

 κμi=Si+^μ⋅Si×Si−^μ. (5)

The chiral glass susceptibility is then given by

 χμCG(k)=N[⟨∣∣qμc(k)∣∣2⟩]av, (6)

where the chiral overlap is given by

 qμc(k)=1N∑iκμ(1)iκμ(2)ieik⋅Ri. (7)

We define the chiral correlation lengths by

 ξμc,L=12sin(kmin/2)(χCG(0)χμCG(kmin)−1)1/2, (8)

in which is independent of . Note that will, in general, be different for along (the direction) and perpendicular to , though we expect that this difference will vanish for large sizes. We denote these two lengths by and respectively.

To equilibrate the system in as small a number of sweeps as possible, with the minimum amount of CPU time we perform three types of Monte Carlo move:

1. “Microcanonical” sweepsAlonso et al. (1986),
(also known as “over-relaxation” sweeps). We sweep sequentially through the lattice, and, at each site, compute the local field on the spin, . The new value for the spin on site is taken to be its old value reflected about , i.e.

 S′i=−Si+2Si⋅HiH2iHi, (9)

see Fig, 1. These sweeps are microcanonical because they preserve energy. They are very fast because the operations are simple and no random numbers are needed. For reasons that are not fully understood, it also seems that they “stir up” the spin configuration very efficientlyCampos et al. (2006) and the system equilibrates faster than if one only uses “heatbath” updates, described next. Figure 1: Hi is the local field on site i due to its neighbors. The spin at i is initially in direction Si. In a microcanonical (over-relaxation) move, the spin is reflected about Hi according to Eq. (9), and so ends up in direction S′i.
2. “Heatbath” sweepsOlive et al. (1986).
Since the microcanonical sweeps conserve energy they do not equilibrate the system. We therefore also include some heatbath sweeps since these do change the energy, typically doing one after every 10 microcanonical sweeps. As for the microcanonical case, we sweep sequentially through the lattice. Referring to Fig. 1, we take the direction of , the local field on site, to be the polar axis, denote the polar angle as , and define the azimuthal angle such that for the old spin direction. The new spin direction is characterized by angles and , relative to , as follows. The energy does not depend on the azimuthal angle, and so , where is a random number chosen uniformly between 0 and 1. The polar angle is chosen such that, after the move, the spin is in local equilibrium with respect to the local field, i.e. if , then

 P(x)=βHi2sinhβHieβHix, (10)

where . To determine with this probability, the procedureOlive et al. (1986) is to equate the cumulative distribution

 Q(x)=∫x−1P(x′)dx′=eβHix−e−βHieβHi−e−βHi, (11)

to a second random number (also with a uniform distribution between 0 and 1) so

 x=1βHiln[1+r2(e2βHi−1)]−1. (12)

It is then necessary to convert the new spin direction back to cartesian coordinates. Denoting the polar and azimuthal angles of by and relative the cartesian reference frame, and remembering that and are relative to , we havecom

 S′x = S′′xcosϕH−S′′ysinϕH, (13) S′y = S′′xsinϕH+S′′ycosϕH, (14) S′z = cosθcosθH−sinθsinθHcosϕ, (15)

where

 S′′x = cosθsinθH+sinθcosθHcosϕ, (16) S′′y = sinθsinϕ. (17)

We see that the calculations in the heatbath moves are quite involved, which is why we do mainly microcanonical moves, just including some heatbath moves to change the energy and thereby ensure the algorithm is ergodic. Note, though, that the acceptance probability for the heat bath moves, and also for the microcanonical moves, is unity, so no moves are wasted.

3. “Parallel tempering” sweeps.
At low temperatures spin glasses are easily trapped in minima (valleys) of the free energy. In order to ensure that the system visits different minima with the correct Boltzmann weight during the time of the simulation we use the method of parallel temperingHukushima and Nemoto (1996); Marinari (1998). One takes copies of the system with the same bonds but at a range of different temperatures. The minimum temperature, , is the low temperature where one wants to investigate the system (below in our case), and the maximum, , is high enough that the the system equilibrates very fast (well above in our case). A parallel tempering sweep consists of swapping the temperatures of the spin configurations at a pair of neighboring temperatures, and , for with a probability that satisfies the detailed balance condition. The Metropolis probability for this isHukushima and Nemoto (1996)

 P(T swap)={exp(−ΔβΔE),(if ΔβΔE>0),1,(otherwise), (18)

where and , in which is the energy of the copy at temperature . In this way, a given set of spins (i.e. a copy) performs a random walk in temperature space. Suppose that at some time in the simulation a copy is trapped in a valley at low-. Later on it will reach a high temperature where it randomizes quickly, so that when, still later, it is again at a low temperature, there is no reason for it to be in the same valley that it was in before. We do one sweep of temperature swaps after every ten microcanonical sweeps.

Table 1 gives the parameters of the simulations. It will be seen that the number of temperatures is very large, and it is instructive to discuss the reason for this. The difference between two neighboring temperatures, , must be sufficiently small that there is an overlap between the energy distributions at those temperatures, so there is a reasonable probability that the same spin configuration occurs for both temperatures. Otherwise, the probability for accepting the temperature swap, Eq. (18), will be very small. Relating the width of the temperature distribution to the specific heat in the normal way, one finds that

 δTT≲1√CN, (19)

where is the specific heat per spin. For the Heisenberg spin glass tends to a constant at low- because of the Gaussian (spinwave) fluctuations about the local equilibrium positions. Hence, for a given size, we choose temperatures which decrease in a geometric manner. The required number of temperatures to cover a fixed range of increases proportional to , and so, for large , many copies are needed to cover even a factor of 2 in . However, for the Ising spin glass with Gaussian interactions, the specific heat tends to zero as (presumably linearly in ). Hence much bigger steps in can be taken at low- than for the Heisenberg case, leading to fewer temperatures being needed.

To test for equilibrationKatzgraber et al. (2001) we require that data satisfy the relationKatzgraber and Young (2002)

 ql−qs=2zTU, (20)

which is valid for a Gaussian bond distribution. Here is the average energy per spin, is the “link overlap”, , is the number of nearest neighbor bonds, and is the lattice coordination number. Equation (20) is easily derived by integrating by parts the expression for the average energy with respect to , noting that the average is over a Gaussian function of the ’s. Figure 2: (Color online) Equilibration plot, testing Eq. (20) for L=24 at the T=0.144. It is seen that the data for 2TU/z come together at about 3×105 total sweeps (equilibration plus measurement) and then stay at their common value indicating that equilibration has been achieved. The lines are guides to the eye. It is seen that the energy comes close to its equilibrium value very quickly, whereas ql−qs takes much longer. Figure 3: (Color online) A plot of ξL/L as a function of the total number of sweeps for L=24 at T=0.144. It is seen the data flattens off at around 3×105 sweeps, the value where the two sets of data in Fig. 2 start to agree. This indicates that when the data in Fig. 2 agree within high precision, i.e. when Eq. (20) is satisfied, the correlation length has reached its equilibrium value. The line is a guide to the eye. Figure 4: (Color online) A plot of ξ∥c,L/L as a function of the total number of sweeps for L=24 at T=0.144. As for Fig. 3 it is seen the data flattens off at around 3×105 sweeps, the value where the two sets of data in Fig. 2 start to agree. The line is a guide to the eye.

The spins are initialized in random directions so the energy, the RHS of Eq. (20), is initially close to zero and decreases, presumably monotonically, to its equilibrium value as the length of the simulation increases. Hence the RHS of Eq. (20) will be too large if the simulation is too short to equilibrate the system. On the other hand, the LHS of Eq. (20), will be too small if the simulation is too short because it starts off close to zero and then increases with MC time as the two replicas start to find the same local minima. The quantity will be less dependent on Monte Carlo time than since it is a local variable for a single replica. (For the Ising case it is just a constant.) Hence if the simulation is too short the LHS of Eq. (20) will be too low. In other words, the two sides of Eq. (20) approach the common equilibrium value from opposite directions as the length of the simulation increases. Only if Eq. (20) is satisfied within small error bars do we accept the results of a simulation.

Figure 2 shows a test to verify that Eq. (20) is satisfied at long times. For the parameters used, , this occurs when the total number of sweeps () is about . Figures 3 and 4 show that the spin and chiral correlation lengths appear to become independent of , and hence are presumably equilibrated, when is larger than this same value. Hence, it appears that when Eq. (20) is satisfied to high precision, the data for the correlation lengths is equilibrated.

With the number of sweeps shown in Table 1, Eq. (20) was satisfied for all sizes and temperatures. The error bars are made sufficiently small by averaging over a large number of samples. We are simulating system sizes which are very large by spin glass standards (up to ), so it is crucial to have a stringent test like this for equilibration.

Since is dimensionless it has the finite size scaling formBallesteros et al. (2000); Palassini and Young (1999); Lee and Young (2003b)

 ξLL=˜X(L1/ν(T−TSG)), (21)

where is the correlation length exponent. Note that there is no power of multiplying the scaling function , as there would be for a quantity with dimensions. There are analogous expressions for the chiral correlation lengths. From Eq. (21) it follows that the data for for different sizes come together at . In addition, they are also expected to splay out again on the low- side if there is spin glass order below . In a marginal situation, with a line of critical points as in the KTB transition, the data for different sizes would come together at and then stick together at lower , see, for example, Fig. 3 of Ref. Ballesteros et al., 2000. Figure 5: (Color online) Data for ξL/L, the spin glass correlation length divided by system size, as a function of T for different system sizes. Note that there are very many data points for the larger sizes. This is because a large number of temperatures are needed for the parallel tempering algorithm, as discussed in Sec. II. Figure 6: (Color online) Enlarged view of a region of Fig. 5. Figure 7: (Color online) Data for ξ∥c,L/L, the “parallel” chiral glass correlation length divided by system size, as a function of T for different system sizes. Figure 8: (Color online) Enlarged view of a region of Fig. 7.

## Iii Results

We studied sizes from to , as shown in Table 1. The CPU time involved to get this data is about 15 Mac G5 CPU years.

We start with the data for shown in Fig. 5. As was found earlier by LY, the smaller sizes show a clear intersection, and also splay out at lower temperatures which would indicate spin glass order. However, the situation for the largest sizes is less clearcut, with the and data only coming together at a somewhat lower temperature than the temperature where the data intersect. Furthermore the data for the largest sizes at the lowest temperatures does not depend strongly on size, indicating close to marginal behavior. This is qualitatively in agreement with Campos et al.Campos et al. (2006). An enlarged view of the important region is shown in Fig. 6.

Data for the parallel chiral glass correlation length is shown in Fig. 7. The main features are the same as found for the spin glass correlation length in Fig. 5. The smaller sizes show clear intersections, but the data lies lower, though perhaps not quite to the same extent as in Fig. 5. As for the spin glass correlation length, the two larger sizes and 24 only come together at a lower temperature than the temperature (range) where the smaller sizes intersect, and do not splay out at still lower temperatures, at least in the range of studied and within the error bars. An enlarged view is shown in Fig. 8. In our view, the data in Figs. 5 and 7 are not very different; in particular, in both figures, the data for the largest two sizes merge at about the same temperature and stick together at lower temperatures. If we look at the largest three sizes there is a somewhat greater tendency for the chiral data to splay out at low-. Clearly larger sizes are needed to be sure of the trend in the thermodynamic limit.

With this data it is difficult to come to a firm conclusion about the nature of the transition. Clearly the smaller sizes have large corrections to scaling behavior, and so, at best, only the two largest sizes are in the asymptotic scaling region. The data is consistent with a line of critical points analogous to that in the KTB transition, as proposed by Campos et al., in which the data for large sizes would merge at a single temperature and remain independent of size at lower temperatures. This is a scenario in which the lower critical dimension is equal to three. Another scenario consistent with our data, in which is also equal to , is that but with an exponential divergence of the correlation lengths at . In that case, it is likely that the data for large sizes would join a common curve, but the temperature at which the common curve is joined would decrease as the system size increased. This scenario is found in the two-dimensional Heisenberg ferromagnet, whereas the critical line of the KTB theory is found in the two-dimensional XY-ferromagnet. Our data is also consistent with the possibility that, for the large sizes, the curves weakly intersect, implying a finite and a lower critical dimension slightly less than three.

Unfortunately, because of the crossover effects in the data, it does not appear possible to give a meaningful estimate of the critical exponents, and .

## Iv Conclusions

In agreement with Campos et al.Campos et al. (2006) and Hukushima and KawamuraHukushima and Kawamura (2005), we find a crossover in behavior for system sizes larger than about 16. Results from the larger sizes indicate a more marginal behavior than those from the smaller sizes. Unlike Kawamura and collaboratorsHukushima and Kawamura (2005); Campbell and Kawamura (2007), we do not find that this effect is very different for the spins and chiralities. Compare, for example, Figs. 5 and 7. However, larger sizes would be needed to confirm that the asymptotic behaviors of the spin and chiral glass correlation lengths are indeed similar. Our data extends up to , somewhat larger than the sizes () studied by Hukushima and Kawamura. Our range of sizes is the same as that of Campos et al., but we are able to go down to lower temperature, in particular below the putative spin glass transition temperature. Whereas Campos et al. argue in favor of a critical line below , in our view, other possibilities can not be ruled out, such as a transition at a lower value of or even an exponential divergence of the correlation length at .

To distinguish between these scenarios would require a study of still larger sizes. It will be difficult to go to very much larger sizes without a better algorithm, since the present study used a quite substantial amount of computer time. An unfortunate feature of the present algorithm is that parallel tempering for vector spin models requires a large number of temperatures. The large number arises from the temperature independent specific heat at low temperatures, which, in turn, comes from from a rather trivial feature: Gaussian spinwave fluctuations. However, it is difficult to see how to eliminate their effect, and thereby reduce the number of temperatures.

It is interesting to note that we (and also Campos et al.Campos et al. (2006)) have been able to study larger sizes for the Heisenberg spin glass () than has ever been done for the Ising spin glass. For example, Katzgraber et al.Katzgraber et al. (2006) studied the Ising spin glass using considerable CPU time, but were still not able to equilibrate sizes larger than near for the case of Gaussian interactions. With interactions, where “multispin coding” speeds the code up further, they could go up to .

That one can equilibrate larger sizes in the Heisenberg case is surprising bearing in mind that (i) the Heisenberg algorithm is more complicated and so one sweep takes more CPU time than for the Ising case, and (ii) the transition temperature is significantly lower for the Heisenberg case. Noting that the mean field transition temperature is , where is the number of neighbors and the number of spin components, if we take for the Heisenberg case, the ratio is about . For the Ising case, one finds  e.g. Ref. Katzgraber et al., 2006, so which is more than twice the corresponding ratio for the Heisenberg spin glass.

The fact that, despite all this, one can study larger sizes in Heisenberg spin glasses than in Ising spin glasses indicates that barriers are smaller in the Heisenberg case. Another way of putting this is that the extra degrees of freedom in the Heisenberg model allow the system to find a way round barriers, which would have to be gone over for the Ising spin glass.

Overall, our results indicate that spins and chiralities behave in a quite similar manner, and that the lower critical dimension of the Heisenberg spin glass is close to, and possibly equal to, three.

###### Acknowledgements.
We acknowledge support from the National Science Foundation under grant DMR 0337049. We are also very grateful to the Hierarchical Systems Research Foundation for a generous allocation of computer time on its Mac G5 cluster.
{"ft_lang_label":"__label__en","ft_lang_prob":0.8993534,"math_prob":0.97477293,"size":28773,"snap":"2020-34-2020-40","text_gpt3_token_len":7028,"char_repetition_ratio":0.15676597,"word_repetition_ratio":0.043669958,"special_character_ratio":0.24602926,"punctuation_ratio":0.15117304,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.97578424,"pos_list":[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16],"im_url_duplicate_count":[null,1,null,1,null,1,null,1,null,1,null,1,null,1,null,1,null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2020-08-10T10:46:15Z\",\"WARC-Record-ID\":\"<urn:uuid:ab407f0b-6034-4579-9ed3-19e4100d7eab>\",\"Content-Length\":\"441588\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:dd141732-9986-4e6a-a69f-5693d9c05b5e>\",\"WARC-Concurrent-To\":\"<urn:uuid:15ec58a5-76b9-4872-a76e-c27e14f44352>\",\"WARC-IP-Address\":\"172.67.158.169\",\"WARC-Target-URI\":\"https://www.arxiv-vanity.com/papers/cond-mat/0703770/\",\"WARC-Payload-Digest\":\"sha1:YQMQHUQQUKH6LGDF7VHJHOQTMKFYCSOP\",\"WARC-Block-Digest\":\"sha1:USA7MRQ6IAVT56PEOFROW2TABYIRKIJJ\",\"WARC-Identified-Payload-Type\":\"text/html\",\"warc_filename\":\"/cc_download/warc_2020/CC-MAIN-2020-34/CC-MAIN-2020-34_segments_1596439738674.42_warc_CC-MAIN-20200810102345-20200810132345-00484.warc.gz\"}"}