Disordered systems near quantum critical points

Computer Physics Communications 121–122 (1999) 505–509 www.elsevier.nl/locate/cpc

Disordered systems near quantum critical points Heiko Rieger 1 Institut für Theoretische Physik, Universität zu Köln, 50923 Köln and NIC, Forschungszentrum Jülich, 52425 Jülich, Germany

Abstract Various aspects of disordered systems at and near quantum phase transitions (QPT), which are transitions at zero temperature that are driven by quantum instead of thermal fluctuations, are discussed. We describe a cluster algorithm in continuous (imaginary) time for general (pure or random) ferromagnetic Ising models in a transverse field. It works directly with an infinite number of time-slices in the imaginary time direction, avoiding the necessity to take this limit explicitly. We present some results for the so-called Griffiths–McCoy regions surrounding the QPT in magnetic systems, where various susceptibilities are found to be divergent at T → 0 without the system being critical.  1999 Elsevier Science B.V. All rights reserved.

The presence of quenched disorder has drastic effect on various characteristics of physical systems. The critical properties and with them the universality class of a phase transition might be changed, the dynamics gets usually much slower due to pinning effects and glassy low temperature phases might occur with features completely different from those of a pure system. Close to a so called quantum phase transition (QPT), a transition at zero temperature that is not driven by thermal fluctuations but by quantum fluctuations, the influence of quenched disorder is even more enhanced: it can happen that in addition to the effects described above singular, if not divergent, behavior of thermodynamic quantities can be observed not only at but also near a critical point. The physical origin of such regions of divergent behavior of for instance susceptibilities close to a QPT are similar to the Griffiths phase [1] in classical systems: strongly coupled clusters that can be arbitrarily large (so-called rare events) tend to order locally giving rise to very large contributions to the response to an appropriate external field and to metastable 1 E-mail: [email protected].

states is an extremely long life time, leading to slow non-exponential dynamics. At finite temperatures the emergent static singularities are hard to observe (experimentally and numerically), only dynamical effects can be noticed, sometimes as stretched exponential decay, sometimes as enhanced power laws of time dependent correlation functions [2]. At zero temperature however, when quantum effects become dominant and statics and dynamics become inextricably linked and a slow dynamics leads automatically to strongly enhanced susceptibilities. Moreover, activated dynamics via quantum tunneling is even slower than classical, thermally activated dynamics, which in turn leads even to divergences in the susceptibilities. To be specific we consider here model for a disordered quantum magnet with a strong Ising anisotropy and transverse fields: X X Jij σiz σjz + Γi σix (1) H=− hij i

i

spin- 21

operators located on any dwhere σi are the dimensional lattice, hij i indicates all nearest neighbor pairs on this particular lattice, Jij > 0 ferromagnetic interactions and Γi transverse fields that are randomly

0010-4655/99/$ – see front matter  1999 Elsevier Science B.V. All rights reserved. PII: S 0 0 1 0 - 4 6 5 5 ( 9 9 ) 0 0 3 9 3 - 8

506

H. Rieger / Computer Physics Communications 121–122 (1999) 505–509

distributed for instance according to a uniform distribution Jij ∈ [0, 1] and Γi ∈ [0, Γ ] modeling a quenched (i.e. time-independent) disorder. In one space dimension (d = 1) the scenario described above, QPT and Griffiths–McCoy region, has been worked out in much detail [3,4]. Various powerful analytical and numerical tools, which are specific for one space dimension, facilitate the thorough study of this paradigmatic system. For the randomness specified above, the system posses a second order phase transition at zero temperature T = 0 and Γ = 1 from a paramagnetic to a ferromagnetic (ground) state. The critical behavior turns out to be very unusual, the dynamics being ultra-slow [5], i.e. the dynamical critical exponent zcrit being infinite. Moreover, away from the critical point the time dependent correlations still decay algebraically although the spatial correlation length is finite. In this so called Griffiths–McCoy region [1,6] near the QPT various susceptibilities are still divergent depending on a particular parameter z(δ), the dynamical exponent in the Griffiths–McCoy phase, which characterizes the strengths of the singularities encountered in a distance δ away from the QPT [7]. In higher dimensions the analytical tools that were successful in 1d fail and one has to rely on numerical methods, in particular quantum Monte Carlo methods. In what follows we describe a Monte Carlo cluster algorithm presented the first time in [8] that is particularly suited to handle the inherent difficulties in the study of such a random quantum system: The origin of the Griffiths–McCoy singularities are strongly coupled clusters (with strong ferromagnetic bonds and weak transverse fields) which are extremely hard to equilibrate in a conventional Monte Carlo algorithm, therefore we had to use a cluster method. Moreover, the exponent z(δ) is a non-universal quantity for which reason we really have to perform the so-called Trotterlimit explained below. A continuous time algorithm that incorporates this limit right from the beginning (in the spirit of Refs. [9,10]) is the most efficient method we can think of. Note however, that this algorithm also performs very well for the pure case [8]. To derive our continuous time algorithm we use the Suzuki–Trotter decomposition to represent the free energy F of the system (1) at inverse temperature β = 1/T as the limit of a (d + 1)-dimensional classical Ising model [11]:

F = −β −1 lim ln Tr exp(−Sclass ), 1τ →0 X Kij Si (τ )Sj (τ ) Sclass = − τ,hij i

−

X

Ki0 Si (τ )Si (τ + 1).

(2)

τ,i

Here the additional index τ = 1, . . . , Lτ of the now classical Ising spin variables Si (τ ) = ±1 labels the Lτ d-dimensional (imaginary) time slices within which spins interact via Kij = 1τ Jij and among which they interact with strength Ki0 = − 12 ln tanh 1τ Γi . The number of time slices Lτ is related to 1τ by 1τ = β/Lτ , so that the limit 1τ → 0 in (2) implies Lτ → ∞. Taking the limit 1τ → 0 consecutive spins with the same value along the imaginary time direction, e.g., Si (τ ) = Si (τ + 1) = · · · = Si (τ + N), form continuous segments S i {[τ, τ + t]} of length t = N · 1τ rather than individual lattice points. Since we are going to take this limit implicitly we will have to consider these imaginary time segments as the dynamical objects in a Monte Carlo algorithm, and not the individual spin values at discrete imaginary times any more. These segments correspond to continuous, uninterrupted pieces of a spin’s world line during which it has the same value, say +1, and its two ends, in the following called cuts, are those times when this spin changes to another value, say −1. As the next step we apply the scheme of the Swendsen–Wang cluster update method [12] within the aforementioned implicit continuous time limit. Remember that in this method in order to construct the clusters to be flipped at random, neighboring spins pointing in the same direction are connected with a certain probability p: for neighbors in the space direction, for instance Si (τ ) and Sj (τ ) with hij i being nearest neighbors, it is pij = 1 − exp(−2Kij ) = 21τ Jij + O(1τ 2 ) and for neighbors in imaginary time, for instance Si (τ ) and Si (τ +1), it is pi0 = 1 −exp(−2Ki0 ) = 1 − 1τ Γi + O(1τ 2 ). These spin connection probabilities are now translated into probabilities for creating (cutting) and connecting segments. The probability for connecting spins along the imaginary time direction at a particular site i over a finite time interval of length t < β is given by the probability to set t/1τ bonds, i.e. in the limit 1τ → 0 0t /1τ

pi

= (1 − 1τ Γi )t /1τ → exp(−Γi t) .

(3)

H. Rieger / Computer Physics Communications 121–122 (1999) 505–509

507

Fig. 1. Sketch of the cluster construction procedure: (a) Configuration of segments before update – full (broken) lines correspond to worldline segments with spin up (down). (b) Insertion of new cuts (crosses) in addition to old ones (bars) according to a Poisson process along each worldline described by (3). (c) Connection of segments with probabilities given by (4). (d) Configuration of segments after assigning randomly a spin value to each cluster of connected segments. Before starting a new cycle (with A) the redundant cuts within segments are removed.

This means for each site one generates new cuts in addition to the old ones from the already existing segments via a Poisson process with decay time 1/Γi along the imaginary time direction. Next one connects segments on neighboring sites i and j that have the same state and a non-vanishing time-overlap t, e.g., S i {[t1 , t2 ]} and S j {[t3 , t4 ]} with t now the length of the interval [t1 , t2 ] ∩ [t3 , t4 ]. The probability for not connecting the two segments is given by (1 − pij )t /1τ = (1 − 21τ Jij )t /1τ → exp(−2Jij t) (4) in the limit 1τ → 0. Finally, one identifies clusters of connected segments and assigns each of them (and herewith all of the segments belonging to it) one value +1 or −1 with equal probability. In Fig. 1 the individual steps of this cluster update procedure are illustrated. The measurement of observables is straightforward: for instance the local magnetization mi at site i for one particular configuration of segments is simply given by the difference between the total length of all + segments and the total length of all − segments, divided by β. The expectation value for the local susceptibility is Zβ χi (ω = 0) =



 dτ σiz (τ )σiz (0) QM = β m2i MC ,

(5)

0

where h· · ·iQM is the quantum mechanical expectation value for model (1) and h· · ·iMC is an average over all configurations generated during the Monte Carlo run. For an actual implementation of the algorithm one has to provide sufficient memory in which the information about the segments S i {[t1 , t2 ]} (site index

i, starting point t1 , end point t2 and state +1 or −1) is stored in linked lists. The number of segments, which of course fluctuates, is typically of the order O(βΓmax Ld ), whereas in discrete time algorithm the number of spin values to be stored is βLd /1τ , which diverges in the limit 1τ → 0. As an example for the application of our algorithm we present results for the Griffiths–McCoy region of the two-dimensional random transverse Ising ferromagnet (TIFM) described by the Hamiltonian (1). The QPT of this system turns out to be located at Γc ≈ 4.4 [8], so that for Γ larger than this value spatial correlations should be short ranged. However, due to the presence of strongly coupled regions in the system the probability distribution of excitation energies (essentially inverse tunneling times for these ferromagnetically ordered clusters) becomes extremely broad. As a consequence we expect the probability distribution of local susceptibilities to have an algebraic tail at T = 0 [4,7] Ω(ln χlocal ) ≈ −

d ln χlocal , z(Γ )

(6)

where Ω(ln χlocal ) is the probability for the logarithm of the local susceptibility χi at site i to be larger than ln χlocal . The dynamical exponent z(Γ ) varies continuously with the distance from the critical point and parameterizes the strengths of the Griffiths–McCoy singularities also present in other observables. At finite temperatures Ω is chopped off at β, and close to the critical point one expects finite size corrections as long as L or β are smaller than the spatial correlation length or imaginary correlation time, respectively. We used β 6 1000 and averaged over at least 512 samples. In Fig. 2 we show data of Ω(ln χlocal )

508

H. Rieger / Computer Physics Communications 121–122 (1999) 505–509

Fig. 2. The integrated probability distribution of the local susceptibility χi (ln ω = 0) of the 2D random TIFM in the Griffiths–McCoy region at Γ = 6.0 (left) and Γ = 7.0 (right). The systems sizes L are larger than the spatial correlation lengths, which are small than far away from the critical point (Γ > Γc ≈ 4.2). There is, however, a dependence on the system size in the imaginary time direction, i.e. β, the inverse temperature, since here the system has long ranged correlations. Therefore β has to be chosen large enough to read off the correct slope of the tail of the probability distribution.

at Γ = 6.0 and Γ = 7.0. From the asymptotic slopes 0.39 and 0.54, respectively, one derives via (6) the estimates z(Γ = 6) = 5.1 and z(Γ = 7) = 3.7, respectively. Collecting the estimates for z(Γ ) for other values of Γ one gets strong indications for a divergence limΓ →Γc z(Γ ) = ∞, as already observed in the onedimensional case. This concurs also with investigations at the critical point, where zcrit is also found to be infinite [14]. One is confronted with the intriguing question, why the one-dimensional and the twodimensional random bond ferromagnet in a random transverse field behave very similar, at the quantum critical point as well as in the Griffiths–McCoy phase. One might speculate, from the viewpoint of the real space renormalization group treatment carried through for the one-dimensional case [3], that the disorder in the two-dimensional case renormalizes via the decimation of bonds and sites to a kind of randomness that is reminiscent of bond or site dilution, for which a percolation transition occurs where z can be shown to diverge [13]. Thus, bearing such a RG picture in mind, also in the fully connected random bond ferromagnet the critical point is governed by a percolation fix-point as the diluted case and the physics emerging here (quantum activated dynamics at the critical point, z(δ → ∞) etc.) should be very similar. Work in this direction is in progress [16]. Concluding we also would like to mention that a similar phenomenon can be observed in another sys-

tem: the disordered Boson–Hubbard model for the superconductor-to-insulator transition in amorphous superconducting films [15]. A phase very similar to the Griffiths–McCoy phase in magnetic systems occurs here, it is the Bose glass phase that is not superconducting any more due to disorder induced localization of the bosons but still compressible (in contrast to the Mott-insulator). The one-particle Greens-function is actually also singular in this phase and the corresponding dynamical exponent z is identical to the critical exponent zcrit . Acknowledgements The author expresses his sincerest gratitude to his collaborators F. Iglói, N. Kawashima, J. Kisker, T. Ikegami, S. Miyashita and A.P. Young. This work was supported by the DFG. References [1] R.B. Griffiths, Phys. Rev. Lett. 23 (1969) 17. [2] M. Randeria, J.P. Sethna, R.G. Palmer, Phys. Rev. Lett. 54 (1985) 1321; A.J. Bray, Phys. Rev. Lett. 54 (1988) 720. [3] D.S. Fisher, Phys. Rev. Lett. 69 (1992) 534; Phys. Rev. B 51 (1995) 6411. [4] A.P. Young, H. Rieger, Phys. Rev. B 53 (1996) 8486; F. Iglói, H. Rieger, Phys. Rev. Lett. 78 (1997) 2473; F. Iglói, H. Rieger, Phys. Rev. B 57 (1998) 11404.

H. Rieger / Computer Physics Communications 121–122 (1999) 505–509 [5] H. Rieger, F. Iglói, Europhys. Lett. 39 (1997) 135; F. Iglói, H. Rieger, Phys. Rev. E 58 (1998) 4238. [6] B.M. McCoy, Phys. Rev. Lett. 23 (1969) 383. [7] M.J. Thill, D.A. Huse, Physica A 15 (1995) 321; H. Rieger, A.P. Young, Phys. Rev. B 54 (1996) 3328; M. Guo, R.N. Bhatt, D.A. Huse, Phys. Rev. B 54 (1996) 3336. [8] H. Rieger, N. Kawashima, Phys. Rev. Lett., submitted (1998); cond-mat/9802104. [9] B.B. Beard, U.J. Wiese, Phys. Rev. Lett. 77 (1996) 5130. [10] N.V. Prokof’ev, B.V. Svistunov, I.S. Tupistyn, JETP Lett. 64 (1996) 911; preprint cond-mat/9703200 (1997). [11] M. Suzuki, Progr. Theor. Phys. 56 (1976) 1454.

509

[12] R.H. Swendsen, J.-S. Wang, Phys. Rev. Lett. 58 (1987) 86. [13] T. Senthil, S. Sachdev, Phys. Rev. Lett.77 (1996) 5292; T. Ikegami, S. Miyashita, H. Rieger, J. Phys. Soc. Jpn. 67 (1998) 2761. [14] H. Rieger, N. Kawashima, in preparation; C. Pich, A.P. Young, cond-mat/9802108. [15] J. Kisker, H. Rieger, Phys. Rev. B 55 (1997) 11981R; Physica A 246 (1997) 348. [16] N. Kawashima, H. Rieger, in preparation; D. Huse, D. Fisher, private communication.