**Home****Retrieve****Search****Browse****Feedback****Help**

**Previous article** **Next article** **Up to abstract**

- Efstratios Manousakis and Román Salvador
*Supercomputer Computations Research Institute, Florida State University, Tallahassee, Florida 32306*

Using an extended version of Handscomb's Monte Carlo method we study the spin-*(1/2* quantum antiferromagnetic Heisenberg model on square lattices of various sizes. The calculated correlation length associated with antiferromagnetic order grows very rapidly as the temperature is lowered suggesting an essential singularity. Our results are consistent with those recently observed by neutron scattering done on Lo* _{2}*CuO

©1988 *The American Physical Society.*

**PACS:** 75.10.Jm 74.20.-z

The recently discovered high temperature superconductors^{1} exhibit interesting magnetic properties. The La_{2}CuO_{4-y}_material has a susceptibility anomaly at a Neaael temperature T_{roman}N _ which is sensitive to the value of y, increasing from T_{roman}N apeq 0 _ for y=0 to T_{roman}N apeq 295 _ K for y = 0 .^{2} Subsequent neutron scattering experiments^{3} show that the materials order antiferromagnetically and in a La_{2}CuO_{4}_single crystal there are instantaneous\p two dimensional (2D) antiferromagnetic correlations^{4} exceeding 200 Aang for T apeq 200 ndash 300 K, with no average staggered magnetization. Neutron scattering^{4} and Raman scattering from magnon pairs^{5} provide a large value for the antiferromagnetic (AF) coupling J apeq 10^{3} K.

There are theoretical studies which examine the possibility of superconductivity mechanisms originating purely from electronic degrees of freedom. In some of these studies the magnetic properties of the materials are one of the central points. A common point of departure is the Hubbard model in its strong coupling limit.^{6,7} In this limit and at half filling this model is equivalent\p to the spin cf12 antiferromagnetic Heisenberg model (AFHM)

H = J sum from { langle ij rangle } S vec_{i}cdot S vec_{j}, (1)

where langle ij rangle denotes nearest neighbor unit cells in the Cu-O plane and S vec_{i}_ is the spin operator of the conduction band electron located at the i th cell. This model is expected to describe the dynamics of spin fluctuations in the undoped La Cu-O material.

In this paper we simulate the quantum spin cf12 AFHM (1) in 2D using Handscomb's method as it has been extended for antiferromagnets.^{8} We perform the calculation on various size lattices and measure among other quantities the spin correlation function. The extracted correlation length xi (T) grows very rapidly as we lower the temperature. For the temperature varying from app 0.7 J to app 0.4 J , xi grows from app 1.5 to app 13 unit cells. The behavior of the correlation length suggests an essential singularity. The detailed behavior of the correlation function is different from that obtained by perturbative renormalization group analysis of the classical Heisenberg model. A dramatic growth of correlations is also revealed by the analysis of the neutron scattering experiments^{4} on La_{2}CuO_{4}_at room temperatures. If we take the value of J of the order of 10^{3} K, extrapolation of our results at T apeq 200 ndash 300 K gives correlation lengths of the same magnitude as those observed.^{4}

The Hamiltonian (1) can be written as^{8}

- beta H = ( beta J / 2 ) sum from { langle ij rangle } ^ (h_{ij^{2}-}h_{ij}) + const, (2)

where const= NJ/2 (N =number of units cells). The operator h_{ij}_ is equal to S_{i^{+}S}sub j^{-}+ S_{i^{-}S}sub j^{+}_, and flips antiparallel spins and gives zero in the case of parallel spins. This Hamiltonian, when derived from the strong coupling Hubbard model at half filling (see Ref. 7), describes processes in which the electron hops to a near\%est neighbor cell occupied by an electron of opposite spin, making the site doubly occupied momentarily, and in the final state both return either to the original configuration (h_{ij^{2}_}term) or to the one with spins exchanged (h_{ij_}term).

Following Ref. 8, any observable O hat can be calculated as

langle O hat rangle == {Tr ( O hat e^{{}- beta H } )} over {Tr ( e^{{}- beta H } )} = { size -2 sum_{r=0^{inf}^}size -2 sum_{{C}sub r } pi (C_{r}) OMEGA (C_{r})} over { size -2 sum_{r=0^{inf}^}size -2 sum_{{C}sub r } pi (C_{r})} , (3a)

pi (C_{r}) = (-1)^{{r}sub 1 } {( beta J / 2 )^{r}} over {r !} ^ Tr ( Q_{{}i_{1}} Q_{{}i_{2}} cdot cdot cdot Q_{{}i_{r}} ),

OMEGA (C_{r}) = {Tr ( O hat Q_{{}i_{1}} Q_{{}i_{2}} cdot cdot cdot Q_{{}i_{r}} )} over {Tr ( Q_{{i}sub 1 } Q_{{}i_{2}} cdot cdot cdot Q_{{}i_{r}} )} , N (3b)

where i_{n}_ denotes a link, langle ij rangle for instance, and Q_{{i}sub n } = h_{ij^{2}_}or h_{ij}_. C_{r}= lcurl i_{1}, i_{2},..., i_{r}rcurl _ is a sequence of r = r_{1}+ r_{2}_ operators and r_{1}_ (r_{2})_ is the number of h's (h^{2}'s) in the sequence. If the h operators in C_{r}_ do not form closed loops the trace in zero. Hence for a square lattice the number of h's in a string must be even and consequently pi (C_{r})_ is nonnegative. The trace of any string of operators is the product of the traces of all the clusters. A cluster is a set of lattice sites connected by operators. An isolated site is also a cluster (monomer). Given an operator sequence forming a cluster, we can construct either two possible states contributing to the trace or none. In the latter case the particular sequence is not allowed. The trace over the entire lattice is equal to either zero or 2^{{}n_{c}} _ (n_{c}= _total number of clusters).

A Markov chain generating a distribution pi (C_{r})_ of sequences C_{r}_ is the following. At each iteration of the random walk we can add or remove any number n_{a}_ or n_{d}_ of operators, respectively. We decide to add or delete an operator with equal probability. We select a given operator to be added with probability 1/4N and the specific location in a string with r operators with probability 1/(r+1). We remove a given operator with probability 1/r . The acceptance probability for a transition from the state C_{r}_ having r operators to the state C_{{r^{sprime}}^{sprime}_}having r^{sprime}= r + n_{a}- n_{d}_ operators and satisfying the detailed balance is given by

P ( C_{r}-> C_{{}r^{sprime}}^{sprime}) = min[ 1 , w ( 1 over 4N right )^{{r}- r^{sprime}} {pi ( C_{{}r^{sprime}}^{sprime})} over { pi ( C_{r}) } right ] , (4)

where the factor w = cf12 , 2, or 1, for r = 0 , r^{sprime}= 0 , or r != 0 and r^{sprime}!= 0 , respectively.

We have tested our program by comparison to the exact one dimensional case. We have also calculated all the obseverables calculated in Ref. 8 and we agree completely. Our main interest is the calculation of the staggered spin correlation function

G ( tau vec ) = (-1)^{{}tau_{x}+ tau_{y}} N^{-1}sum_{n}langle S_{z}( n vec ) S_{z}( n vec + tau vec ) rangle .

We checked that the calculated correlation function satisfies the susceptibility sum rule.

We have performed calculations on lattices with sizes 10^{2}, 20^{2}, and 30^{2} with periodic and open boundary conditions (BC). The number of iterations that we performed depends on the temperature and lattice size. For the higher temperatures and smaller lattices, we performed 0.5 times 10^{6} iterations for thermalization, and 10^{6} iterations for measurements. For the lower temperatures and bigger lattices longer runs were required for both thermalization and measurements. The quantities which required more iterations to reach equilibrium were r_{1}_, and the G( tau ). For example, for the 20^{2} lattice and at temperature 0 we performed 8 times 10^{6} iterations for thermalization and 12 times 10^{6} iterations for measurements. Simulations at lower temperatures or larger lattices are beyond realistic computational time scales.

We fitted the staggered G by the forms

lim from {tau -> inf } G ( tau ) = { lpile { { A e^{{}- tau / xi (T) } , for open BC, } above {A ^ cosh( {tau - L/2 } over { xi (T)} right ) , for periodic BC, } } (5)

where L is the size of the lattice. Figure 1(a) shows G calculated for the above three lattices for T=J with periodic and open BC. The lines are obtained by our fitting all except the first few points of G by the forms of Eqs. (5). As may be observed in Fig. 1(a), G is approximated very well by Eqs. (5) for several orders of magnitude and the slopes are quite independent of the BC and finite size effects. As T is lowered and the correlation length becomes comparable to the lattice size, the calculation of xi requires larger lattices. We checked that the projected correlation function, G_{p}(x) = (1/L ) size -2 sum_{y}G (x,y) _, gives the same correlations lengths within error bars. We found that open BC are more appropriate for the measurement of correlation lengths. In Fig. 1(b) we present G for our 20^{2} lattice for various T. The correlations grow very rapidly in the range 1.0 J ndash 0 , giving xi (J) = 0.9 and xi (0.4J) = 13 .

In Fig. 2 we give an equilibrium configuration of clusters in the 20^{2} lattice for T = 0.5 J (top figure) where xi app 3.5 . The filled (open) circles denote up (down) spins. The clusters are drawn by solid lines. There is a large cluster involving most of the lattice sites and some other smaller ones. The lower part of Fig. 2 shows the clusters for 10^{2} lattice and temperatures 0 (left hand side) and 1 (right hand side). In the 0 case we see almost a Neaael configuration and only two clusters, one being monomer and another containing 99 sites. At T = 1.5 J we have a rather disordered state with lots of small clusters.

The correlation length xi (T) is plotted in Fig. 3 on a logarithmic scale (see also Table I). The dashed line gives xi (T) obtained from the leading order contribution to G at high temperatures G(n) = ( beta J / 4 )^{n} and xi (T) = 1 / ln (4T/J). Note that the agreement at high temperature is remarkable. In the range of temperatures 0.5 J < T < 10J the numerical results are independent of lattice size. The behavior is clearly not linear; in fact, the slope -d ^ ln xi (T) / d ( ln T ) increases rapidly with decreasing T and it is app 7 at T = 0.4 J .

Now we compare this behavior with existing theories. An exponential growth of the correlation length as a function of 1/T at low T is predicted by perturbative renormalization group (PRG) theory^{9} of the classical Heisenberg model, giving xi (T) = A(T) e^{{}2 pi g J/T }, where the preexponential factor at low temperatures is A(T) = CJ/ T. This form does not fit the data except in a small range of temperatures. Choosing this range to be 0 < T < 0 we obtain C = 0 and g = 0 and the fit is shown in Fig. 3 by the dashed dotted line. This value of g is reduced from its classical value g_{roman}cl = S^{2}_, which within PRG theory could be accounted for by quantum fluctuations. The fit is not very good which indicates either that the onset of the PRG behavior is at lower T or that nonperturbative effects such as topological configurations change the behavior of xi (T). A well known example is the XY model^{10} where a phase with topological order exists because of vortices. Topological excitations are known to exist in the 2D classical Heisenberg model^{11} also. It is believed, however, that in the classical model they do not alter the behavior predicted by PRG.^{9} In the quantum AFHM the structure of the ground state and the changes in xi (T) which such excitations may cause are theoretically unknown. We also tried to fit the correlation length by a Kosterlitz Thouless^{10} form xi (T) = A ^ exp (B/|T-T_{c}|^{nu})_. The fit is good (solid line), giving nu = 0.5 , the same as in the Kosterlitz Thouless case, and T_{c}= 0 _ (and A = 0.178 , B=1.338). As a result of the Mermin Wagner theorem,^{12} the average staggered magnetization is zero at T != 0 . However, the theorem does not exclude a transition to a phase where G decays algebraically ( xi = inf ) below some finite T_{c}_ because of topological order. This fit could be fortuitous; the reader, therefore, should take this \fIonly\fR as an indication that topological instanton con\%figurations similar to those found in Ref. 11 might play an important role in our case.

Regardless of any theoretical interpretation, our results are consistent with the neutron scattering data.^{4} The neutron- and Raman scattering data give a value for the AF coupling of J apeq 10^{3} K. In the neighborhood of room temperature, i.e., app 0.3 J , a reasonable extrapolation gives correlation lengths of the order of those observed^{4} or higher. The average staggered magnetization at finite temperature in 2D is zero as expected.^{12}

The three dimensional (3D) AF ordering of La_{2}Cu_- O_{4-y}_happens at a lower temperature scale than the AF coupling J. This can be explained as a result of both the weak layer coupling inherent in these materials and the special crystalline arrangement which tends to frustrate a 3D order. The orthorombic distortion pre\%sumably relieves some frustration and produces three dimensional Neaael order at T_{roman}N app 200 _ K.

The authors are grateful to Dr. B. Berg, Dr. G. Bhanot, Dr. D. Duke, Dr. K. Huang, Dr. D. Petcher, Dr. Y.-L. Wang, and Dr. P. Weisz for useful discussions. This work was supported in part by the Florida State University (FSU) Supercomputer Computations Re\%search Institute (SCRI) which is partially funded by the U. S. Department of Energy through Contract No. DE FC05 85ER250000. A portion of the calculation was performed on the ETA 10 supercomputer at the FSU SCRI.

\fINote added.\fRmdashWhile this work was in press we obtained the detailed behavior^{13} of the correlation length as a function of temperature. A detailed comparison with the data as well as implementation of 3D effects due to weak interlayer coupling (producing a Neaael temperature T_{roman}N apeq 200 _ K) will be given elsewhere.^{14}

- [1] J. G. Bednorz and K. A. Muumlller, Z. Phys. B 64, 188 (1986); S. Uchida, H. Takagi, K. Katasawa, and S. Tanaka, Jpn. J. Appl. Phys. Pt. 2 26, L1 (1987); C. W. Chu, P. H. Hor, R. L. Meng, L. Gao, Z. J. Huang, and Y. Q. Wang, Phys. Rev. Lett. 58, 405 (1987); R. Cava, R. B. van Dover, B. Batlogg, and\p E. A. Rietman, Phys. Rev. Lett. 58, 408 (1987).
- [2] H. Thomann, P. Tindall, and D. C. Johnston, to be published.
- [3] D. Vaknin, S. K. Sinha, D. E. Moncton, D. C. Johnston,\p J. M. Newsam, C. R. Safinya, and H. E. King, Jr., Phys. Rev. Lett. 58, 2802 (1987).
- [4] G. Shirane, Y. Endoh, R. J. Birgeneau, M. A. Kastner,\p Y. Hidaka, M. Oda, M. Suzuki, and T. Murakami, Phys. Rev. Lett. 59, 1613 (1987).
- [5] K. B. Lyons, P. A. Fleury, J. P. Remeika, A. S. Cooper, and T. J. Negran, Phys. Rev. B 37, 2353 (1988).
- [6] J. E. Hirsch, Phys. Rev. Lett. 54, 1317 (1985); P. W. Anderson, G. Baskaran, Z. Zou, and T. Hsu, Phys. Rev. Lett. 58, 2790 (1987); S. Kivelson, D. S. Rokhsar, and J. P. Sethna, Phys. Rev. B 35, 8865 (1987); A. E. Ruckenstein, P. J. Hirschfeld, and J. Appel, Phys. Rev. B 36, 857 (1987).
- [7] K. Huang and E. Manousakis, Phys. Rev. B 36, 8302 (1987); E. Kaxiras and E. Manousakis, Phys. Rev. B 37, 656 (1988).
- [8] D. H. Lee, J. D. Joannopoulos, and J. W. Negele, Phys. Rev. B 30, 1599 (1984).
- [9] S. H. Shenker and J. Tobochnik, Phys. Rev. B 22, 4462 (1980).
- [10] J. Kosterlitz and D. Thouless, J. Phys. C 6, 1181 (1973);\p J. Kosterlitz, J. Phys. C 7, 1046 (1974).
- [11] A. A. Belavin and A. M. Polyakov, Pis'ma Zh. Eksp. Teor. Fiz. 22, 503 (1975) [JETP Lett. 22, 245 (1975)].
- [12] N. D. Mermin and H. Wagner, Phys. Rev. Lett. 22, 1133 (1966).
- [13] Y. Endoh \fIet al.\fR, to be published.
- [14] E. Manousakis and R. Salvador, to be published.

- TABLE I. The correlation length for various lattices and temperatures, and for open and periodic BC.
- FIG. 1. (a) Staggered correlation function for lattices of size 10^{2}, 20^{2}, and 30^{2} for T=J. (b) Correlation function for a 20^{2} lattice at various temperatures. Error bars smaller than the diameter of the open circles are emitted.
- FIG. 2. Top: Clusters in a 20^{2} lattice at equilibrium and at T = J . Bottom: Clusters in a 10^{2} lattice at T = 0 (left hand side) and at T=1 (right hand side). We only indicate one of the two possible spin states of the cluster.
- FIG. 3. Correlation length as a function of temperature for lattices of sizes 10^{2}, 20^{2}, and 30^{2} and open BC. Error bars smaller than the diameter of the open circles are omitted. Dashed line is xi (T) = 1 / ln (4T / J ), the leading contribution at high T. Dashed dotted and solid lines are as explained in the text.

View Page Images, PDF (326 kB), or ASCII text (18 kB)

**Previous article** **Next article** **Up to abstract**

**Home****Retrieve****Search****Browse****Feedback****Help**