University of Southern California Dissertations and Theses

Quantum phase transitions in disordered antiferromagnets

(USC Thesis Other)

Quantum phase transitions in disordered antiferromagnets

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content QUANTUM PHASE TRANSITIONS IN DISORDERED ANTIFERROMAGNETS by Rong Yu A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (PHYSICS) August 2007 Copyright 2007 Rong Yu Dedication To my family ii Acknowledgments First I would like to express my gratitude to my PhD advisor, Professor Stephan Haas for his guidance on my research. He is a great mentor who gives me tremendous supports and encouragements during my PhD study. I would also like to thank Dr. Tommaso Roscilde who gives me great help and a lot of insightful advices on my dissertation. I would like to thank Professor Hubert Saleur for his support and many fruitful dis- cussions. I also thank Dr. Noah Bray-Ali, Dr. Weifei Li, Dr. Omid Nohadani, Dr. Bruce Normand, Dr. Pinaki Sengupta, Dr. Chushun Tian, Dr. Stefan Wessel, Amy Cassidy, Letian Ding, Wen Zhang, and all my colleagues and friends for their many invaluable discussions and great support. I express my love and gratitude to my grandma, Yulan Cao, my parents, Yuanjiang Yu and Shuhua L¨ u, and my little sister, Li Yu, for their love and encouragement throughout my life. Finally I would like to thank Professor Gene Bickers, Professor Werner D¨ appen, Professor Vitaly Kresin, and Professor Chongwu Zhou for serving my PhD advisory committee. iii Table of Contents Dedication ii Acknowledgments iii List of Figures vi Abstract xiii Chapter 1 Introduction 1 Chapter 2 Monte Carlo Simulation 6 2.1 Classical Monte Carlo ........................................................................................ 6 2.2 Quantum Monte Carlo ....................................................................................... 9 2.3 Stochastic Series Expansion Algorithm ........................................................... 11 2.3.1 SSE with Loop Operator Update ............................................................ 12 2.3.2 Directed Loops ........................................................................................ 22 2.3.3 Estimators of Thermodynamic Quantities in SSE .................................. 28 Chapter 3 Percolative Quantum Phase Transition in Two-Dimensional Bond-Diluted Antiferromagnets 34 3.1 Introduction ...................................................................................................... 34 3.1.1 Classical Percolation ............................................................................... 37 3.1.2 The Quantum Percolation Scenario ........................................................ 41 3.2 Inhomogeneous Bond Dilution of the QHAF on the Square Lattice ............... 46 3.3 Details on Numerical Methods ........................................................................ 49 3.4 Magnetism on the Percolating Cluster ............................................................. 51 3.5 Magnetic Phase Diagram and Critical Properties ............................................ 55 3.6 Criticality Near the Ladder Limit .................................................................... 62 3.7 The Quantum-Disordered Phase: Correlations and Griffiths-McCoy Singularities .....................................................................................................69 3.7.1 Short-Range Correlations ....................................................................... 69 iv 3.7.2 Griffiths-McCoy Singularities ................................................................ 75 3.8 Summary .......................................................................................................... 83 Chapter 4 Field-Induced Quantum Disordered Phases in a Two-Dimensional Spin-Gapped System with Site Dilution 86 4.1 Introduction ...................................................................................................... 86 4.2 Models and Numeric Method .......................................................................... 92 4.3 The Phase Diagram in Weak-Coupling Limit ................................................. 96 4.4 Field-Induced DFM Phase ............................................................................. 101 4.4.1 Classical vs Quantum Views on DFM Phase ....................................... 102 4.4.2 Local Clusters of Free Moments ........................................................... 104 4.4.3 Statistics on Local Clusters and Pseudo-Plateaus ................................. 113 4.5 Phase Transitions in 2D Site-Diluted Coupled Dimers ................................. 123 4.5.1 The Phase Diagram with Site Dilution ................................................. 123 4.5.2 Quantum Scaling and The Critical Exponents ...................................... 127 4.6 Summary ........................................................................................................ 133 Chapter 5 Conclusions 136 Bibliography 139 v List of Figures 2.1 The vertex representation of local Hamiltonian matrix elementW(b p ). The local bondb p is associated with sitei andj, and at propagation level p in the operator sequenceS L . Occupied (open) circles refer to spin up (down) states. Here we takeα i = S z i , i.e., in the eigenstate ofS z i . Part of the operator loop passing through this vertex is also shown. After the loop update, both the local operatorH bp and the local spin states are changed. The updated vertex is shown in the right. . . . . . . . . . . . 20 2.2 The only possible six vertices in S = 1/2 XXZ model. Their corre- sponding matrix elements are given in Eq. 2.28. . . . . . . . . . . . . . 20 2.3 Four out of the eight subsets of vertices. Each subset shown here cor- responds to an independent set of directed-loop equations. Upper-left: k = 1; upper-right: k = 2; Lower-left: k = 3; lower-right: k = 4. The other four subsets give no more new directed-loop equations, and can be obtained by interchanging the up and down spins while keeping the arrows. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.1 The evolution of order parameterM as concentration factor P . Three possible scenarios of percolative phase transition in a disordered quan- tum system are shown. A: the quantum phase transition takes place at the classical percolation thresholdP cl c , and belongs to the same univer- sality class as the classical percolation; B: the quantum phase transition takes place at the classical percolation threshold but with different criti- cal exponents; C: the quantum phase transition takes place at a different concentration than the classical percolation threshold. . . . . . . . . . . 35 3.2 A sketch of the percolating cluster in the nodes-links-blobs model. The backbone lattice is homogeneous at length scaleL>ξ. . . . . . . . . . 40 vi 3.3 A segment of ladder in the percolating cluster leads to local RVB states (yellow ellipsis). The spins involved in the RVB state are very weakly correlated to the rest of the cluster, causing an effective fragmentation of the percolating cluster into non-percolating subclusters. . . . . . . . 44 3.4 (a) Decomposition of the square lattice into dimer (D) and ladder (L) bonds; (b) randomly coupled dimer limit; (c) randomly coupled ladder limit. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.5 Inhomogeneous bond dilution favors the formation of local strongly fluctuating quantum states. The green ellipses indicate such local dimer- singlet, plaquette-singlet and ladder-like RVB states: (a)P L > P D ; (b) P D >P L . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.6 Upper panel: classical phase diagram for the inhomogeneous bond- percolation problem on the square lattice. The angleθ parametrizing the critical curve is indicated. Lower panel: scaling plots for the network strengthN c /N in the two strongly inhomogeneous casesP D = 0.9 and P L = 0.9. An excellent data collapse is realized, using the 2D per- colation exponentsβ = 5/36 and ν = 4/3, and with P (cl) L,c = 0.4007 (P D = 0.9) andP (cl) D,c = 0.0107 (P L = 0.9). . . . . . . . . . . . . . . . 53 3.7 Upper panel: scaling of the disorder-averaged squared staggered mag- netization on the percolating cluster with fixed size N c . The differ- ent curves correspond to various points along the classical percolation transition. The continuous lines represent quadratic fits. Lower panel: extrapolated thermodynamic values of the staggered magnetization from the upper panel as a function of the inhomogeneity parameter θ. The dashed line is a guide to the eye. . . . . . . . . . . . . . . . . . . . . . 55 3.8 Phase diagram of the spin-1/2 QHAF on the inhomogeneously bond- diluted square lattice. . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.9 Percolating clusters for the bond-diluted square lattice with size24×24 and periodic boundary conditions: (a)P D = 0.15,P L = 0.8; (b)P D = 0.5, P L = 0.5 (homogeneous case); (c) P D = 0.9, P L = 0.45; (d) P D = 1.0,P L = 0.44 (randomly coupled dimers). . . . . . . . . . . . . 58 vii 3.10 (A) Scaling plot of ξ x (open symbols), ξ y (full symbols) for P D = 1 (randomly coupled dimers). (B) Critical exponentsz andν for the mag- netic transition in the inhomogeneously bond-diluted S=1/2 Heisenberg antiferromagnet. The angleθ parametrizes the critical line in Fig. 3.8. Open (full) symbols refer to estimates obtained through the scaling ofξ x (ξ y ). The shaded areas correspond to the regions of parameters in which the transition is quantum. All lines are guides to the eye. . . . . . . . . 59 3.11 Cartoon of the ordering mechanism in randomly coupled ladders. When P D = 1/hri≈ 1/ξ ladder (see text) percolating strings of block spins (rec- tangles) appear in the system. The divergence of the correlation length along such strings drives the onset of long-range correlations between the strings and ultimately of 2D long-range order. . . . . . . . . . . . . 63 3.12 Correlation length of the bond-diluted 128× 2 ladder compared with the average ladder lengthhli. The data forξ ladder are taken atβ = 2048; comparison with data atβ = 4096 (not reported) shows no significant deviation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.13 Dominant effects of bond dilution on a two-leg ladder: (a) local antifer- romagnetic modulation due to a missing rung bond; (b) effective cou- plings between spins missing their rung bond; (c) enhancement of the singlet component (indicated by an ellipse) of the two rung dimers adja- cent to a missing leg bond. . . . . . . . . . . . . . . . . . . . . . . . . 65 3.14 Phase diagram of the QHAF on the inhomogeneously bond-diluted square lattice (see also Fig. 3.8) close to the ladder limit and compared with the single-ladder predictionP D,c = 1/ξ ladder (P L ) (see text). . . . . . . . 68 3.15 Temperature and finite-size scaling of the static structure factorS(π,π) for a representative point in the quantum-disordered regime (P D = 0.15, P L = 0.85). Dashed lines are a quartic polynomial fit inT = β −1 [see Eq. (3.20)] to extrapolate to theβ → ∞ limit (when required). Inset: β → ∞ extrapolated values of the static structure factor plotted vs. system size, clearly showing saturation in the thermodynamic limit. . . 70 viii 3.16 Left panel: temperature and finite-size scaling of the correlation length along the two lattice directions (ξ x andξ y ) for a representative point in the quantum-disordered regime (P D = 0.15,P L = 0.85). Dashed lines are a quartic polynomial fit inT =β −1 to extrapolate toβ→∞ (when required). Right panel: β → ∞ extrapolated values of the correlation length(s) vs. system size; the dashed lines are polynomial fits up to second order in L −1 on the points for L ≥ 48 to compensate for the finite-size effects. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 3.17 Finite-size scaling of the T = 0 static structure factor and correlation length(s) for a point in the quantum-disordered regime (P D = 0.15, P L = 0.9) closer to the ladder limit. The dashed lines associated with the correlation length data are quadratic fitsξ(L) =ξ(∞)+a/L+b/L 2 . 72 3.18 Real-space images of the bond energy (both panels), local static stag- gered structure factorS(i) (see text) (a) and of the static staggered cor- relationsC(0,i) between a reference point on a dangling bond and the remainder of the system (b), for the percolating cluster on a 32× 32 sample withP D = 0.15 andP L = 0.85 atβ = 8192. The thickness of each bond is proportional to its energy, whereas the radius of the dots in (a) and (b) is proportional to the local value ofS(i) andC(0,i). In (a) red boundaries highlight the regions where local correlations exceed the average value, namely whereS(i)&S(π,π). In (b) the red bound- aries mark the correlation volume beyond which C(0,i) ≤ 8× 10 −3 , corresponding to a distance (in units of the correlation length)|r i |/ξ≈ ln[C(0,0)/C(0,i)]≈ 3.4 . . . . . . . . . . . . . . . . . . . . . . . . . 73 3.19 Low-temperature uniform susceptibility of the percolating cluster in the quantum-disordered regime: ladder regime (upper panel) and dimer regime (lower panel). The solid lines are power-law fits of the form χ u,c ∼ T −1+α , and the resulting fit coefficients α are indicated. In the upper panel, the percolating cluster is picked as the largest cluster in a64×64 lattice, whereas in the lower panel it is grown freely from a seed site up to a fixed site ofN c = 1024 (see Sec. 3.3). . . . . . . . . . . . . . . . . 76 3.20 Sketch of a low-energy excitation in the bond-diluted ladder system. (a) Four spins with missing rung bonds are coupled across a clean region of lengthl with weak effective couplingsJ eff ∼ exp[−σl]; (b) rotating the spins at one end of the clean region by the same angle, spin excitations occur on the weaker bonds only, leaving the stronger ones unchanged. . 78 ix 3.21 Probability distribution of the local susceptibility for a representative point in the quantum-disordered phase (P D = 0.15,P D = 0.9), taken at β = 128L. The slope of the linear tail of the distribution is proportional to theα exponent in the thermodynamic limit (see text). . . . . . . . . . 80 4.1 A sketch of the phase diagram of the site-diluted weakly coupled dimers in the presence of a magnetic field. The labels on theh axis denotes the critical fields at each phase transition (see text in Sec. 4.2 and Sec. 4.3 for details). OBD: order-by-disorder phase. . . . . . . . . . . . . . . . 88 4.2 The field dependence of the uniform magnetization, uniform suscepti- bility, and staggered magnetization in 2D site-diluted weakly coupled dimers. The coupling ratioJ ′ /J = 1/4, and the dilution concentration p = 1/8. In the lower plot the orange diamonds show the staggered magnetization data in the thermodynamic limit, which are extracted from the finite-size scaling results. SF: the transverse ordered phase, which corresponds to a superfluid phase in the bosonic language. . . . . 97 4.3 The temperature dependence of the uniform susceptibility in 2D site- diluted weakly coupled dimers at the first magnetization PP h/J = 0.175. The coupling ratioJ ′ /J = 1/4, and the dilution concentration p = 1/8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 4.4 (a): The classical S → ∞ free moments has a canted antiferromag- netic order at finite field. (b): In the quantum case, the close-lying free moments prefer to form local clusters, leading to a novel DFM phase. The black and red arrows distinguish spins located on different sublattices.103 4.5 Field dependence ofm u and transverse spin componentS x(y) in the clas- sical effective spin model withξ 0 = 1.0,p = 0.0625 andL = 40. The correspondingm u in the quantum model with the same model parame- ters is also given as a comparison. . . . . . . . . . . . . . . . . . . . . 105 4.6 Distribution of the z component of local moments S z i in the effective spin model for classical spins withξ 0 = 1.0,p = 0.0625, andL = 40. . 106 4.7 Distribution of the z component of local moments S z i in the effective spin model for quantum spins withξ 0 = 1.0,p = 0.0625, andL = 40. . 108 4.8 Upper panel: distribution of the logarithm of the reduced bond energy h ¯ E b i in the effective spin model withξ 0 = 1.0,p = 0.0625, and system sizeL = 40. Lower panel: a closer look in the vicinity oflnh ¯ E b i = 0 and−0.4, indicating the formation of close-lying dimers and trimers. . . 109 x 4.9 Finite-size study oflnh ¯ E b i ath/J ′ = 0.008 in the effective spin model withξ 0 = 1.0 and p = 0.0625. The inset shows finite-size scaling of P(lnh ¯ E b i); the red line is a linear fit. . . . . . . . . . . . . . . . . . . . 111 4.10 Histogram of bond energies (a) and local magnetizations for unpaired spins (b) in the diluted weakly-coupled-dimer model. The model para- meters areJ ′ /J = 1/4,p = 1/8, and system sizeL = 40. . . . . . . . . 112 4.11 The distribution of local cluster gaps Δ(C) in the effective bosonic model, withξ 0 = 1.0 andp = 1/8. Each peak in the distribution can be related with an effective finite-size cluster sketched in the plot. The peak centered at Δ(C) ≈ 1.9J ′ comes from more complicated cluster geometries, such as quadrumers. . . . . . . . . . . . . . . . . . . . . . 119 4.12 Field dependence of the normalized magnetizationm u /m sat u calculated by appying Eq. 4.24 in the effective model with ξ 0 = 1.0 and p = 1/8 (shown in blue triangles). The result is compared with QMC result from the effective model (shown in red circles) with the same model parameters, and with QMC result from the diluted coupled-dimer model (shown in black squares) withJ ′ /J = 1/4 andp = 1/8. . . . . . . . . 121 4.13 Field dependence of the normalized magnetizationm u /m sat u in the effec- tive model withξ = 1.0, and various dilution concentrationp. Dashed blue lines show the values of 1/(1+4p) as guidelines for the magneti- zation of the first PP. . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 4.14 Phase diagram of site-diluted coupled dimers with dilution concentra- tion p = 0.125. It coincides with the phase diagram of 2D interact- ing bosons in the presence of a random potential. The three distinct phases: gapped disordered phase (the blue regime), gapless disordered phase (the red regime), and transverse ordered phase (the white regime) corresponds to the Mott insulating (MI), Bose glass (BG), and super- fluid (SF) phases in the bosonic system, respectively. In weak cou- pling limit J ′ ≪ J, the gapless disordered phase is separated by the gapped phase into two branches, corresponding to the DFM phase and the intermediate-field BG phase. The phase boundary in the clean limit p = 0 (black dashed line) is also shown in the same plot as a comparison. 125 4.15 Finite size scaling of the correlation length in the vicinities of critical h DFM (left) andh BG (right). The inter-dimer couplingJ ′ = J/4, and dilution concentrationp = 0.125. . . . . . . . . . . . . . . . . . . . . 128 xi 4.16 Finite size scaling of the staggered magnetization and the spin stiffness in the vicinities of critical fieldsh DFM (left) andh BG (right). The model parameters are the same as in the caption of Fig. 4.15. . . . . . . . . . . 129 4.17 From top to bottom: finite size scaling of correlation length ξ, stag- gered magnetizationm s , and spin stiffnessρ s at critical fieldsh DFM = 0.05(1) (left) and h BG = 0.56(1) (right) in the site-diluted coupled dimers withJ ′ /J = 0.35 andp = 1/8. . . . . . . . . . . . . . . . . . . 132 xii Abstract Recently quantum phase transitions have attracted the interest of both theorists and experimentalists in condensed matter physics. Quantum magnets provide a perfect play- ground for studying these phase transitions since they can be triggered by many control parameters such as frustration, lattice dimerization, and magnetic field. Most previ- ous studies have focused on the magnetic properties in pure systems. In these systems, responses to the triggering parameters are found to be uniform, leading to homogeneous phases. However little progress has been made so far on the phase transitions and prop- erties in disordered quantum magnets because they are more complicated systems, and few theoretical tools can be applied. In this thesis we use the stochastic series expansion quantum Monte Carlo method to study quantum phase transitions in disordered magnets. We find that disordered mag- nets can behave quite differently from pure systems. The system inhomogeneity can strongly affect phase transitions by changing their universality class. We also find order- disorder transitions are often accompanied by the appearance of novel quantum disor- dered phases, in which magnetic properties behave highly nontrivial, even singular. xiii In this thesis two examples are studied in great detail. The first one is the phase diagram of an inhomogeneous, bond-diluted two-dimensional antiferromagnet near the percolation threshold. We show that the magnetic transition can be tuned by the inho- mogeneity of the dilution from a classical percolation to a quantum phase transition. Interestingly the quantum transition still takes the nature of a renormalized percola- tive transition, with continuously varying critical exponents. A gapless quantum dis- ordered phase with no magnetic long-range order but geometric percolation is found. The low-temperature uniform susceptibility diverges as a non-universal power-law of the temperature in this phase, indicating that this is a quantum Griffiths phase. In the second example, we study field-induced quantum phase transitions in two-dimensional site-diluted coupled dimers. We find that an extremely small field drives the order-by- disordered phase into a quantum disordered disordered-free-moment phase. The mag- netization curve in this phase shows a series of pseudo-plateaus, which can be explained by the highly inhomogeneous field response of the free moments. A spin-boson map- ping reveals that the disordered-free-moment phase can be understood as a low-field Bose glass phase, and hence the associated phase transition is a realization of a two- dimensional Bose-glass-to-disordered-free-moment transition. This is confirmed by the entire phase diagram and the study of the scaling behavior of correlation length, order parameter, and spin stiffness in the critical regime. xiv Chapter 1 Introduction Quantum phase transitions (QPTs) have recently become one of the most exciting research topics in condensed matter physics. Contrary to thermal phase transitions, QPTs occur at zero temperature by tuning physical control parameters. This offers new insights into the low-temperature properties of condensed matter. Quantum magnets provide a large variety of materials which can realize various QPTs. Moreover, there are many ways to tune quantum fluctuations which can drive the system through a QPT. For example, in systems of weakly antiferromagnetically coupled dimers there is a singlet ground state which is short-range ordered. Applying a magnetic field may overcome the spin gap. The field excites triplets which can be treated asS = 1 bosons. These bosons (excited triplets) are delocalized and thus condense to the lowest eigenstate at zero tem- perature, giving an example of Bose-Einstein condensation in quantum magnets. Although major effort and progress has been made in studying phase transitions in quan- tum magnets, we notice that most studies to date have considered only pure systems. This is not very surprising since the disordered magnetic systems displays much more 1 complicated behavior and thus few theoretical techniques can be applied. However, from the experimental side the demand for studying disordered spin models is consid- erable since these models are more realistic, and tightly related to experiments. Prac- tically, no material is in a sense pure due to defects, impurities, and chemical substi- tutions. These factors can lead to disordered ground states by introducing geometric randomness into the system. Disordered ground states cannot be treated perturbatively in general because the impurities can be strongly correlated. Sometimes the correlations among impurities are so strong that counterintuitive behavior occurs. A famous example for this may be the order-by-disorder (OBD) phenomenon in the doped coupled ladders. In the pure limit, two-leg ladder has a disordered ground state protected by a finite spin gap proportional to the rung coupling. Hence an array of weakly coupled two-leg lad- ders in three dimensions (3D) does not present any long-range order. Upon doping with non-magnetic impurities, one would expect that this system still be disordered since the impurities reduces the connectivity of the system. However, at small doping concentra- tion, a long-range antiferromagnetic ordered state has been predicted theoretically and later confirmed by numerics and experiments [28, 118, 110]. This surprising result can be understood as follows: each non-magnetic impurity induces an effective free moment localized around the impurity site, which can be regarded as an unpaired spin neighbor- ing to a vacant site. Although the wave function of the free moment is spatially decaying as ψ(r)∼ e −r/ξ 0 r , (1.1) 2 whereξ 0 is the correlation length of the clean system, the overlap of multiple impurity wave functions leads to an effective coupling between the two moments. This effective couplingJ eff is long-ranged and preserves the unfrustrated nature of the system. In a doped 3D system, the free moments form a random network throughJ eff , and present long-range antiferromagnetic order. We have seen through the above example that disordered systems may have totally dif- ferent ground state properties from pure ones. When the system is driven to a phase transition novel phases may appear. In general, these novel phases cannot be treated within a mean-field level. In this thesis we apply the stochastic series expansion (SSE) quantum Monte Carlo (QMC) method to study phase transitions in disordered quantum magnets. The stochastic series expansion algorithm has been proved to be a very pow- erful quantum Monte Carlo technique to investigate the low-temperature properties of quantum spin systems at a large scale. The non-local updates in the algorithm helps avoiding the usual critical slowing down. This makes it suitable for studying quantum phase transitions. We will give an introduction of the SSE QMC method in chapter 2 Two examples of quantum phase transitions in disordered quantum magnets are studied through SSE QMC simulations in this thesis. The first one is the phase transition in a inhomogeneous bond-diluted Heisenberg antiferromagnet in two spatial dimensions (2D) near the percolation threshold. The second one is the field-induced phase transi- tions in a 2D site-diluted spin-gap system. 3 In chapter 3 we investigate the phase diagram of an inhomogeneous bond-diluted Heisenberg antiferromagnet in 2D near the percolation threshold. Percolation occurs in a variety of contexts, [26] ranging from blood vessel formation to clusters of atoms deposited on substrate surfaces. A fundamental question is whether the classical picture of permeating networks [95] applies as well to strongly fluctuating quantum systems. In this respect, low-dimensional quantum magnets offer an ideal playground to probe the interplay between quantum fluctuations and geometric randomness both theoreti- cally and experimentally. Recent investigations [106, 83] have focused on the case of homogeneous site and bond dilution in the quantum antiferromagnet on the square lat- tice, reporting a classical geometric percolation transition between magnetic order and disorder. In this study we show how inhomogeneous bond dilution leads to percola- tive quantum phase transitions, which we have been studied extensively by quantum Monte Carlo simulations. Quantum percolation introduces a novel 2D quantum dis- ordered phase, characterized by an infinite percolating network with vanishing antifer- romagnetic order parameter. In this phase, the low-temperature uniform susceptibility diverges algebraically with non-universal exponents. This is a signature that the novel quantum-disordered phase is a quantum Griffiths phase, as also directly confirmed by the statistical distribution of local gaps. This study thus presents evidence of a genuine quantum Griffiths phenomenon in a two-dimensional Heisenberg antiferromagnet. In chapter 4 we study the phase transition induced by the applied magnetic field in a 2D site-diluted coupled dimer system. As introduced above, in the clean system, the 4 field leads to an order-disorder transition which can be explained within the picture of Bose-Einstein condensation of excited triplets. Doping with non-magnetic impurities completely changes this picture, by introducing novel quantum phases. At zero field, the system orders magnetically due to the order-by-disorder (OBD) mechanism. This is a fundamental difference between the diluted and pure cases. This OBD phase is highly inhomogeneous, and it is easily destroyed by a weak field. The field destruction of the OBD phase leads to a novel quantum disordered disordered-free-moment (DFM) phase. The magnetization process is found to be highly nontrivial, characterized by a series of magnetization pseudo-plateaus (PPs) in the uniform magnetization curve. To explain the PPs we propose a physical scenario based on the formation of local free- moment clusters. This scenario combined with a spin-boson mapping reveals the nature of this DFM phase to be a low-field Bose glass (BG) phase. Hence the OBD-to-DFM transition realizes a Bose-glass-to-superfluid transition in 2D quantum magnets. We confirm this point by investigating the entire phase diagram of the diluted couple-dimer system and studying the scaling behaviors of several thermal quantities in the quantum critical regime. 5 Chapter 2 Monte Carlo Simulation In this chapter we describe a powerful Monte Carlo method, which we will apply to study the 2D disordered antiferromagnets in the following chapters. In section 2.1 we first introduce the basic idea of Monte Carlo simulation and the implementation of Monte Carlo algorithms to classical spin systems. We then describe algorithms that can be applied to quantum systems in section 2.2. In section 2.3 we focus on the stochastic series expansion algorithm, the method we will use throughout this thesis. 2.1 Classical Monte Carlo There are two different methods to study the properties of a statistical system: the deter- ministic method and the stochastic method. In the deterministic method, a certain trajec- tory in the phase space of the system is generated according to the equation of motion. An interested quantity is calculated as a time average along this trajectory. Ergodicity 6 guarantees that this time average equals to the ensemble average. Contrary to the deter- ministic method, Monte Carlo samples the system in a stochastic way, and ensemble average of a certain quantity is obtained based on the sample average. In a classical system, the statistical average of observableO can be expressed as hOi = 1 Z X C O(C)W(C), (2.1) where Z = P C W(C), and W(C) = e −βE C is the weight of configuration C. The Monte Carlo method works as a random walk in the configuration space such that the number of times each configuration C is visited is proportional to the weight W(C). ThenhOi becomes a simple arithmetic average hOi =hOi MC ≡ lim N MC →∞ 1 N MC N MC X i=1 O(C i ), (2.2) whereN MC is the number of Monte Carlo steps. To satisfy the above equation two con- ditions must be fulfilled by the Monte Carlo algorithm. The first one is called ergodicity, i.e., there should be a finite probability to connect any two configurations within a finite number of steps. The second one is the detailed balance condition: W(C)P(C→C ′ ) =W(C ′ )P(C ′ →C), (2.3) whereP(C→C ′ ) is the transition probability from configurationC toC ′ . 7 The detailed balance condition leaves freedom of implementing the algorithm since P(C→C ′ ) is not fully determined by Eq. 2.3. A simple choice with P(C→C ′ ) = min[1,e −β(E C ′−E C ) ] (2.4) leads to the Metropolis algorithm [53]. This algorithm satisfies both ergodicity and detailed balance condition, and is easy to be implemented. The only drawback of it comes from its local nature, and shows up when dealing with phase transitions. Let us consider a phase transition in a spin system. Spins are highly correlated within a cor- relation lengthξ. The configuration is not updated unless aboutξ d correlated spins are simultaneously flipped. But since Metropolis algorithm is a local update, the probabil- ity of a configuration being updated decays exponentially with increasingξ. When the system approaches to the critical point,ξ→∞. This leads toP(C→C ′ )→ 0, i.e., the algorithm becomes very inefficient. This difficulty is called critical slowing down. It can be overcome by various cluster algorithms, such as Swendsen-Wang algorithm [97] and Wolff algorithm [114, 115]. 8 2.2 Quantum Monte Carlo In a quantum system the thermal expectation value of an observableO is expressed as hOi = 1 Z Tr[Oe −βH ], (2.5) where the quantum partition functionZ = Tre −βH . If the system HamiltonianH can be easily diagonalized, it is trivial to reduce the expression in Eq. 2.5 into Eq. 2.1 so that Metropolis algorithm is applicable. But in most cases the diagonalization ofH is not straightforward, and a weight function is not generally well defined. So we need to construct a special configuration space. In this space the Hamiltonian, after some transformation, can be mapped to an effective classical one where a weight function is well defined. This procedure is usually denoted as quantum to classical mapping. The mapping is not unique. It determines the configuration space and the efficiency of an algorithm. Two mapping procedures are often used. One is the stochastic series expansion (SSE) method, which we will introduce in great detail in next section. The other method is through the path-integral formulation. We introduce the mapping by applying it to a 1DS = 1/2 Heisenberg antiferromagnetic system. We first rewrite the Hamiltonian into a sum of several terms H = P n b=1 H b . 9 Each term includes many commutable pieces, but any two terms are not commutable. This can be realized in antiferromagnets by rewriting H =H odd +H even , (2.6) where H odd (H even ) represents the sum of all the local Hamiltonians applying on odd (even) bonds. We then apply the Trotter-Suzuki decomposition [102, 96] of the imag- inary time path-integral formulation of the partition function. With the aid of Trotter formula e A+B = lim K→∞ e A/K e B/K K , (2.7) whereK is the Trotter number, we get Z = Tre −βH = lim K→∞ Tr e −Δτ P n b=1 H b K = Tr   Y b e −ΔτH b ! K   +O(Δτ 2 ) = X α 1 ,···,α nK hα 1 |e −ΔτH 1 |α 2 ihα 2 |e −ΔτH 2 |α 3 i···hα n |e −ΔτHn |α n+1 i··· ···hα (n−1)K+1 |e −ΔτH 1 |α (n−1)K+2 i···hα nK |e −ΔτHn |α 1 i. (2.8) whereΔτ =β/K, and the summation is performed over a complete orthonormal set of states. In this way we map the original quantum mechanical system ind-dimension to an effective(d+1)-dimensional classical model, with one extra imaginary time dimension. The original Hamiltonian is transformed to the transfer matrix in the effective model. 10 Loop cluster algorithm has been developed by constructing closed loops on the(d+1)- dimensional lattice [20, 113]. From Eq. 2.8 we see that for any finiteK, the correspond- ing systematic error of the algorithm is in the order of Δτ 2 . So it becomes exact only when the Trotter numberK →∞. This limit can be realized by the algorithm working in continuous imaginary time, which is called world line algorithm [5]. 2.3 Stochastic Series Expansion Algorithm In this section we introduce a different quantum Monte Carlo (QMC) algorithm, the Stochastic Series Expansion (SSE) algorithm. It is an extension of Handscomb’s method [31, 44, 43, 11], which is based on the power expansion of the partition function. By treating the partition fuction in this way, it totally avoids the systematic error of the Trotter decomposition [96]. In many cases, even when the system is close to a quan- tum critical point, the autocorrelation time is dramatically reduced due to the non-local nature of the SSE algorithm. This makes it a very efficient algorithm in dealing with quantum phase transitions (QPTs). Since it was first introduced 15 years before [80], SSE algorithm has been applied to a wide class of quantum lattice models including general quantum spin-S models with magnetic fields [98], bosonic models [98], and even 1D spinful fermion models [87]. 11 2.3.1 SSE with Loop Operator Update Here we describe the SSE algorithm by applying it to the S = 1/2 XXZ model with an external magnetic field. A detailed description of the SSE method can be found in Ref. [99]. The Hamiltonian of the XXZ model reads H = X hi,ji J ij [S x i S x j +S y i S y j +ΔS z i S z j ]−h X i S z i . (2.9) The single-axis anisotropy is denoted by Δ. Here we consider a rather general case. The spin coupling strengthJ ij can vary from bond to bond. A bond is ferromagnetic (FM) ifJ ij < 0, and antiferromagnetic (AFM) ifJ ij > 0. The model is defined on a d-dimensional lattice with arbitrary geometry. The only restriction is that if not all the bonds are FM, the lattice must be bipartite to avoid any sign problem due to frustration. It is convenient to rewrite the Hamiltonian in Eq. 2.9 in terms of bond operators H b , whereb refers to the bond connecting two sitesi(b) andj(b). H =− N b X b=1 H b , (2.10) 12 whereN b is the total number of bonds on the lattice. Each local bond operator consists a diagonal partH 1b and a off-diagonal partH 2b , H b =H 1b ∓H 2b , (2.11) where the upper (−) sign is for the AFM bond, and the lower (+) sign is for the FM bond. We will keep this convenience in following formulas in this section. The off- diagonal operator H 2b = J b 2 [S + i(b) S − j(b) +S − i(b) S + j(b) ], (2.12) whereJ b =|J ij |> 0. The diagonal part H 1b =C∓J b ΔS z i(b) S z j(b) +h[ S z i(b) z i(b) + S z j(b) z j(b) ], (2.13) wherez i is the coordination number of sitei.C is a constant to ensure that for any arbi- trary bondb,H 1b has only non-negative matrix elements. SinceC is a simple constant in the Hamiltonian, it only shifts the zero point of the energy spectrum. So that it has no effect on the thermodynamics of the system. This constant is conveniently chosen as C = Δ 4 J max + h 2 ( 1 z i + 1 z j )+ǫ, (2.14) whereJ max = max{J 1 ,J 2 ,...,J N b }, and the small constantǫ> 0. 13 We can take the following standard basis |αi =|S z 1 ,S z 2 ,...,S z N i, (2.15) where N denotes the totally number of spins (sites). The partition function is then expanded in this basis as Z = Tr{e −βH } = X {α} ∞ X n=0 (−β) n n! hα|H n |αi (2.16) = X {α} ∞ X n=0 β n n! hα|[ N b X b=1 H b ] n |αi, where β = 1/T is the inverse temperature. Use Eq. 2.11 to expand the Hamiltonian as the sum of a series of operator products consisting of both diagonal and off-diagonal bond operators, we obtain Z = X {α} ∞ X n=0 X Sn (−) n 2 β n n! hα| n Y k=1 H a k ,b k |αi. (2.17) HereS n refers to a sequence indicating the relative order of the bond operators appearing in each operator product. S n = [a 1 ,b 1 ],[a 2 ,b 2 ],...,[a n ,b n ], (2.18) 14 where a k = 1(2) if the k’s operator in the product is diagonal (off-diagonal). b k ∈ {1,2,...,N b } is a bond index. The number n 2 in Eq. 2.17 is the total number of off- diagonal operators on AFM bonds in the sequence S n . In the FM model, it is just trivially zero, and in the AFM model on a bipartite lattice, it is always an even number so that the expansion is positive definite. The partition function can then be sampled without any sign problems. It is easy to see from Eq. 2.16 that the average energyE can be related to the average expansion orderhni through E =hHi =−hni/β. (2.19) So the average expansion orderhni∼βJ max N b . The width of the distribution of expan- sion order is approximatelyhni 1/2 . Thus the expansion in Eq. 2.17 can be truncated at a maximum powern =L. The partition function is thus Z = X {α} X S L β n (L−n)! L! hα| L Y k=1 H a k ,b k |αi. (2.20) To fix the expansion order to beL one may randomly insert unit operatorsH 0,0 into the operator products. The sum is then over all possible ways of inserting unit operators. The error introduced by this truncation can be well controlled to be negligible if L is chosen to be reasonable large such that during each simulation n < L is always 15 guaranteed. In practice one takes L = νn max , where n max is the largest n reached during the simulation, and usually the factorν = 1.2−1.5. From Eq. 2.20 we see that SSE is to sample not only in the Hilbert space{α} but also in the ”operator space” S L . Starting the simulation from an empty operator sequence [0,0],[0,0],...,[0,0] and a random spin configuration, we have to find efficient update algorithms working for bothS L and{α}. Two different types of updating procedures are thus introduced. The first type of update is the diagonal update. This is a local update between diagonal and unit operators [0,0] p ↔ [1,b] p . Each successful attempt changes the expansion ordern by±1. If we define the intermediate state at propagation levelp by acting the firstp operators inS L onto the state|αi: |α(p)i = p Y k=1 H a k ,b k |αi, (2.21) then with this convenience the Metropolis acceptance probabilities, which are based on the detailed balance [81], are P([0,0] p → [1,b] p ) = min[1, N b βhα(p)|H 1b |α(p)i L−n ], (2.22) P([1,b] p → [0,0] p ) = min[1, L−n+1 N b βhα(p)|H 1b |α(p)i ]. (2.23) 16 The presence ofN b in above two probabilities is a direct result of the fact that there are N b possible ways in the substitution[0,0]→ [1,b] but only one way for [1,b]→ [0,0]. The diagonal update are performed consecutively forp = 1,2,...,L. To calculate the matrix elements in Eq. 2.22 and Eq. 2.23, one needs to propagate the spin state|α(p)i during the update procedure. This means spins are flipped once a off-diagonal operator in the sequence is encountered. The second type of update, called operator loop update, tries to make substitutions between [1,b] p and [2,b] p during each loop-like random walk in the configuration space (both spin sate and operator sequence). Although the loop update is global, at each step it only changes the local states of two spins that are connected by a bond. Thus we could focus on a reduced two-spin Hilbert space at each step. For each 16p6L we define the reduced Hilbert space of two spins associated by bondb p in propagating state|α(p)i as |α bp (p)i =|S z i(bp) ,S z j(bp) i. (2.24) Then the partition function in Eq. 2.17 can be rewritten as Z = X {α} ∞ X n=0 X Sn W(α,S n ), (2.25) 17 where the SSE weight function is W(α,S n ) = β n n! n Y p=1 W(b p ), (2.26) and W(b p ) =hα bp (p)|H bp |α bp (p−1)i. (2.27) Here we find that the SSE weight depends only on the matrix element of local Hamil- tonian. For a given model, there are only limited numbers of matrix elementsW(b p ). For example in the spin-1/2 XXZ model we considered here, the non-zero matrix ele- ments are as follows: W 1b ≡h↑↑|H b |↑↑i = Δ 4 (J max ∓J b )+h( 1 z i + 1 z j )+ǫ, W 2b ≡h↓↓|H b |↓↓i = Δ 4 (J max ∓J b )+ǫ, W 3b ≡h↑↓|H b |↑↓i = Δ 4 (J max ±J b )+ h z i +ǫ, (2.28) W 4b ≡h↓↑|H b |↓↑i = Δ 4 (J max ±J b )+ h z j +ǫ, W 5b ≡h↓↑|H b |↑↓i =J b /2, W 6b ≡h↑↓|H b |↓↑i =J b /2. There is a way to visualize these matrix elements. We find that these elements are fully determined by the four spin states|α bp (p)i and|α bp (p−1)i associated with the bondb p . 18 So each element is represented as a vertex consisting of a bar with four legs. As shown in Fig. 2.1, the bar refers to the bond operator H bp , and each leg (open or occupied circle) represents the local spin state before or after acting the operator H bp on it. In Fig. 2.2 vertices corresponding to the matrix elements in Eq. 2.28 are illustrated. The SSE weight is then a network of these vertices on the (d+1)-dimensional lattice. Here d refers to the spacial dimension of the system, and the extra dimension corresponds to the propagation through operator sequenceS L , with linear sizeL. To construct the network, any two vertex legs on the same site should be linked if there is no other vertex in between. We note that any two linked legs share the same local spin state. Since the partition function is a trace over the spin states, we should enforce periodical boundary condition along the propagation direction. In such a way we establish a mapping from the quantum system ind-dimension to a classical system in (d + 1)-dimension. Note that this mapping does not introduce any systematic error of Trotter decomposition. Now we describe the operator loop update. One starts the update by randomly picking up a vertex leg on the d + 1-dimensional lattice of vertices as the initial entrance leg. Then one of the four legs belonging to the same vertex as the entrance leg is selected as the exit leg from that vertex. Both the entrance and the exit spins are flipped. At the same time the vertex itself also changes its type due to the spin flipping. An example of how this works is shown in Fig. 2.1. As an exit has been found the leg linked to the exit becomes the new entrance and the above procedure iterates until the initial entrance leg is reached. There are two possible ways to close a loop. Either at the last step the 19 |S z j (p-1)> |S z i (p-1)> |S z j (p)> |S z i (p)> H b p entrance exit Figure 2.1: The vertex representation of local Hamiltonian matrix elementW(b p ). The local bondb p is associated with sitei andj, and at propagation levelp in the operator sequence S L . Occupied (open) circles refer to spin up (down) states. Here we take α i = S z i , i.e., in the eigenstate of S z i . Part of the operator loop passing through this vertex is also shown. After the loop update, both the local operator H bp and the local spin states are changed. The updated vertex is shown in the right. Figure 2.2: The only possible six vertices inS = 1/2 XXZ model. Their corresponding matrix elements are given in Eq. 2.28. exit leg is the initial entrance leg, or the leg linked to the last exit is the initial entrance leg. Note that when we start the loop we generate a link discontinuity (the two linked spins are in different states). We find that this link discontinuity is healed when the loop closes. We then reach to a new configuration, i.e., a new term in the partition function in Eq. 2.25. Here we see how the loop update works. It samples not only the state space{|αi}, but at the same time the operator sequences{S n }. But since the loop update only makes 20 changes between diagonal and off-diagonal operators, it does not modify the expansion ordern. The sampling overn is completed by applying the diagonal update. There is a special case that cannot be sampled by a loop update. If a spin in state|αi is not acted on by any of the operators in the sequenceS L , its state will not be flipped by any loop. This case only happens whenn is small, or in another word, at high temperatures. Then we may perform an extra spin flipping update at the end of each loop update to flip those spins with probability1/2. A MC step then includes a sweep of diagonal updates at all propagation levels in S L followed by several operator loop updates. The number of loop updates in one MC step N l should be chosen such that a significant fraction of the operators in S L are visited during each MC step. It can be determined self-consistently during the thermalization: tuneN l periodically such that the average number of vertices visited during one MC step is about 2L. However, when taking dataN l must be kept constant in order not to bias the measurements. Typically each loop has a length (number of vertices visited during the loop update) smaller than the expansion ordern. But in very few cases, there could be very long loops. To make the MC algorithm efficient, we must terminate these long loops once they exceed a maximum length∼ 100L. In order not to violate the detailed balance, in fact we should discard the configuration generated during the incomplete MC step, and restart a MC step based on the configuration at last complete step. Termination of loops does introduce bias to observable such as correlation function, but the effect is usually pretty small, and in some models this termination is not necessary [99]. 21 2.3.2 Directed Loops Now the only thing we did not mention on SSE algorithm is how to choose an exit for a given entrance. We will see that different choices give algorithms with different efficiency. But there are some general rules based on the symmetry of the model and the detailed balance. In XXZ model, the symmetry in spin space requires S z i(b) (p)+S z j(b) (p) =S z i(b) (p−1)+S z j(b) (p−1) (2.29) at each vertex. This symmetry limits the only possible vertices to be those shown in Fig. 2.2. Considering this restriction, one may easily check that in spin-1/2 XXZ model only three legs are allowed to be an exit for each type of vertices. Especially ath = 0 and Δ = 1, the exit leg is uniquely determined because of the higher symmetry. The algorithm is thus greatly simplified there. The detailed balance add restriction to the acceptance probability of exiting from a leg. But this restriction is not local since now we have a loop update. Let us denote s a configuration (includes|αi and S L ) with corresponding SSE weight W(s). Then the detailed balance requires P(s→s ′ )W(s) =P(s ′ →s)W(s ′ ), (2.30) 22 whereP(s → s ′ ) is the transfer probability from the configurations tos ′ . Generally, froms tos ′ the loop will experience updates at a few vertices. This leads to P(s→s ′ ) = X {loops} P(l 0 ) n Y k=1 P(s k−1 ,l k−1 →s k ,l k ), (2.31) and P(s ′ →s) = X {loops} P(l 0 ) n Y k=1 P(s k ,l k →s k−1 ,l k−1 ), (2.32) where s 0 = s, s n = s ′ , P(l 0 ) is the probability of having vertex leg l 0 as the initial entrance, andP(s k−1 ,l k−1 → s k ,l k ) is the transfer probability at one vertex from con- figurations k−1 tos k with entrance legl k−1 and exitl k . Since the loop finally closes, we setl n = l 0 . To satisfy the detailed balance the loops in Eq. 2.32 must be the same as those in Eq. 2.31. By comparing Eq. 2.31 and Eq. 2.32 we find that the detailed balance condition in Eq. 2.30 is satisfied if W(s k−1 )P(s k−1 ,l k−1 →s k ,l k ) =W(s k )P(s k ,l k →s k−1 ,l k−1 ) (2.33) at each 1 6 k 6 n, and for all possible SSE configurations. We should mention that Eq. 2.33 is just a sufficient but not necessary condition. The restriction of detailed balance that any SSE algorithm must satisfy is actually Eq. 2.30. But the choice of Eq. 2.33 as detailed balance is convenient. Now we rewrite the left-hand side of Eq. 2.33 23 asW(s)P(s,e→s ′ ,x) =W(s,e,x), wheree andx refer to the entrance and exit legs of the involved vertex respectively. Then the detailed balance condition becomes W(s,e,x) =W(s ′ ,x,e). (2.34) Another restriction comes from the normalization of transfer probabilities X x P(s,e→s x ,x) = 1, (2.35) wherex runs over all legs on the vertex. This leads to X x W(s,e,x) =W(b p (s)), (2.36) where W(b p (s)) is just one of the weights in Eq. 2.28. The equations Eq. 2.36 and Eq. 2.34 are called the directed-loop equations. The solution of these equations com- pletes the SSE QMC algorithm. Now we discuss the solutions of directed-loop equations in spin-1/2 XXZ model. We first observe that all the vertices can be divided into eight subsets that do not transform into each other. Each subset corresponds to a set of directed-loop equations. Consider- ing the symmetries, there are only four sets of independent equations. The subsets of 24 Figure 2.3: Four out of the eight subsets of vertices. Each subset shown here corresponds to an independent set of directed-loop equations. Upper-left:k = 1; upper-right:k = 2; Lower-left: k = 3; lower-right: k = 4. The other four subsets give no more new directed-loop equations, and can be obtained by interchanging the up and down spins while keeping the arrows. vertices correspond to these equations are shown in Fig. 2.3. All the sets of equations have similar structure. For setD k        b k1 a k3 a k2 a k3 b k2 a k1 a k2 a k1 b k3               1 1 1        =        W k1 W k2 W k3        , (2.37) 25 where the diagonal elementsb’s are bounce probabilities, i.e., the probabilities for exit- ing on the same leg as the entrance, and the weights on the right-hand side satisfy            W 11 W 12 W 13 W 21 W 22 W 23 W 31 W 32 W 33 W 41 W 42 W 43            =            W 1b W 4b W 5b W 1b W 3b W 5b W 2b W 3b W 5b W 2b W 4b W 5b            . (2.38) Even when we require the weightsa kl andb kl to be non-negative, there are still infinite solutions of Eq. 2.37 since there are six unknown variables but only three equations. A symmetric solution is given as a kl =W km W kn /(W kl +W km +W kn ) b kl =W 2 kl /(W kl +W km +W kn ), (2.39) where{l,m,n} is a permutation of{1,2,3}. It is called the heat-bath solution [82]. The advantage of the heat-bath solution is it always satisfies Eq. 2.37. But in this solution the bounce probabilities get relatively large values. This becomes severe whenh → 0 and Δ → 1 in theS = 1/2 XXZ model since the bounce probability approaches 1/2. The bounce process does not change the SSE configuration. It reduces the efficiency of the algorithm, and should be maximally suppressed, or even excluded. 26 The idea of minimizing bounce probabilities leads to a optimization of the solutions of directed-loop equations. It results in the directed-loop solution [99]. We may first rewrite the Eq. 2.37 as a kl = 1 2 (−W kl +W km +W kn )+ 1 2 (b kl −b km −b kn ). (2.40) To optimize the solution we want to minimize the bounce probabilities but still keep all the probabilities non-negative. A general way to find the optimal solution is given in Ref. [98] and [2]. Here we give an optimal solution for S = 1/2 XXZ model in the isotropic limit,Δ = 1. In this case we may always setǫ = 0. ForJ ij > 0, i.e., the AFM case, the only nonzero bounce probabilities are b 11 = h z i(b) −J b for h>z i(b) J b , b 21 = h z j(b) −J b for h>z j(b) J b , b 32 = h z i(b) for h> 0, b 42 = h z j(b) for h> 0. (2.41) For the FM case whereJ ij < 0, the non-zero bounce probabilities are b 11 = h z i(b) for h>,0 b 21 = h z j(b) for h>,0 b 32 = h z i(b) −J b for h>z i(b) J b , b 42 = h z j(b) −J b for h>z j(b) J b . (2.42) 27 In both cases the other probabilities can be solved from Eq. 2.40. We note that in the above optimal solution, the bounce processes can be completely excluded in absence of magnetic field. This is called the deterministic case where all the non-zero matrix ele- ments in Eq. 2.28 get the same value 1/2. More efficient algorithm has been developed in this case [82, 32]. 2.3.3 Estimators of Thermodynamic Quantities in SSE In SSE algorithm, measurements of diagonal quantities and quantities which can be expressed as product of bond operatorsH b can be easily obtained. For a local diagonal operatorO i at sitei is measured through hO i i = 1 n n−1 X p=0 hα(p)|O i |α(p)i. (2.43) Especially the magnetization and the static correlation functions are m z = 1 N N X i=1 hS z i i = 1 nN n−1 X p=0 N X i=1 hα(p)|S z i |α(p)i, (2.44) C(i,j) =hS z i S z j i = 1 n n−1 X p=0 hα(p)|S z i |α(p)ihα(p)|S z j |α(p)i. (2.45) 28 If the operator is a bond operatorH bp then hH bp i = 1 Z X {α} ∞ X n=0 X Sn (−β) n n! * α|H bp n Y k=1 H b k |α + . (2.46) We note that in above equationH bp can be absorbed into the product so that the expres- sion is similar as the one for diagonal operators. In this way we see that to find the expectation value of H bp is equivalent to sample N(b p ), i.e., the number of appearing indexb p in the operator sequences: hH bp i =− 1 β hN(b p )i. (2.47) This immediately leads to the estimator for the system energy in Eq. 2.19. In general for the product ofm bond operatorsH b 1 ,···,H bm we have * m Y p=1 H bp + = 1 (−β) m (n−1)! (n−m)! N(b 1 ,···,b m ) , (2.48) whereN(b 1 ,···,b m ) denotes the number of ordered subsequencesb 1 ,··· ,b m in oper- ator sequenceS L . The imaginary-time-dependent operator O 2 (τ)O 1 (0) =e τH O 2 e −τH O 1 (2.49) 29 can also be measured after some mathematics. First we apply a Taylor expansion to the expectation valuehO 2 (τ)O 1 (0)i: hO 2 (τ)O 1 (0)i = 1 Z X {α} ∞ X n=0 ∞ X m=0 (τ−β) n (−τ) m n!m! hα|H n O 2 H m O 1 |αi (2.50) After playing a trick on summing over indexm, we get hO 2 (τ)O 1 (0)i = 1 Z X {α} ∞ X n=0 n X m=0 (τ−β) n−m (−τ) m (n−m)!m! hα|H n O 2 H m O 1 |αi = 1 Z X {α} ∞ X n=0 n X m=0 X Sn (τ−β) n−m (−τ) m (n−m)!m! ×hα| n Y k=m+1 H b k O 2 m Y l=1 H b l O 1 |αi. (2.51) If bothO 1 andO 2 are diagonal, hO 2 (τ)O 1 (0)i = n−1 X p=0 n X m=0 τ m (β−τ) n−m (n−1)! β n (n−m)!m! O 1 [p]O 2 [p+m], (2.52) where O[p] ≡ hα(p)|O|α(p)i. Usually we are more interested in the integrated form R β 0 dτhO 2 (τ)O 1 (0)i. The integral can be easily evaluated using the beta function [3], and the result is Z β 0 dτhO 2 (τ)O 1 (0)i = β n(n+1) n−1 X p=0 O 1 [p] n−1 X p=0 O 2 [p]+ n−1 X p=0 O 1 [p]O 2 [p] ! . (2.53) 30 As an example, the static susceptibility χ(i,j) = Z β 0 dτhS z i (τ)S z j (0)i (2.54) can be estimated through χ(i,j) = β n(n+1) n−1 X p=0 S z i [p] n−1 X p=0 S z j [p]+ n−1 X p=0 S z i [p]S z j [p] ! . (2.55) Another complicated example is the spin stiffnessρ s . It is highly related to the superfluid density in bosonic models [72]. If we introduce a twist Hamiltonian with the twist applied along thex-direction of the square lattice H(φ) = J X hi,jix [cos(φ)(S x i S x j +S y i S y j )+sin(φ)(S x i S x j −S y i S y j )+ΔS z i S z j ] +J X hi,jiy [S x i S x j +S y i S y j +ΔS z i S z j ], (2.56) then ρ s = 1 JL 2 ∂ 2 F(φ) ∂φ 2 φ=0 , (2.57) where F(φ) is the free energy associated to the twist Hamiltonian. The spin stiffness can be expressed as combination of bond operators [81]. Applying the above tricks, it can be written as ρ s = 1 2Jβ hw 2 x +w 2 y i. (2.58) 31 In above equationw x andw y are winding numbers in thex andy directions: w α = (N + α −N − α )/L (α =x,y), (2.59) whereN ± x = P b N ± b,x , andN + b,x andN − b,x are the numbers that operatorsS + i(b) S − j(b) and S − i(b) S + j(b) associated with bondb inx direction appear in the operator sequenceS n . To estimate general off-diagonal quantities needs more effort. The basic idea to measure the off-diagonal quantities is to map the original system to a modified one in which these off-diagonal quantities are certainly well defined. Then a MC run simulating both the original and the modified system simultaneously will measure both the diagonal (from simulating the original system) and off-diagonal (from simulating the modified system) quantities. The above idea can be formulated as follow: for an off-diagonal operator O od , its averagehO od i =Tr[ρO od ]/Z. Now we defineρ ′ =ρO od andZ ′ =Trρ ′ . Then hO od i =Z ′ /Z, i.e., the ratio between the modified and original partition functions. The partition functionZ = P C W(C), whereW(C) is the SSE weight of configurationC. If the modified system can be simulated within the same MC run, then one must have Z ′ = P ′ C W(C), where P ′ reflects the constraint added by operatorO od . So that hO od i = P ′ C W(C) P C W(C) . (2.60) 32 As an example, we consider the following operatorO od = S + i S − j . It can be evaluated during constructing the operator loops by introducing the concept of worm. The worm refers to the path of an incomplete loop in the SSE configuration space. We denote the worm tail as the initial entrance vertex leg and the worm head as the current visiting vertex leg. The worm grows during the construction of the loop, and disappears when the loop closes, i.e., when the worm head meets its tail. To measure S + i S − j , we may define the modified system by attachingS + andS − to the worm head and tail respec- tively. According to the ergodicity and the detailed balance condition, the operator loop samples the partition function Z; and at the same time, the worm with two attached operators at its head and tail samples the modified partition functionZ ′ . Once a loop starts from sitei, the worm head passes sitej, and the head and the tail are at the same propagating level, we get a measurement ofS + i S − j . Based on Eq. 2.60 we obtain hS + i S − j i = N[S + i S − j ] N loop , (2.61) whereN loop is the number of operator loops during the MC simulation andN[S + i S − j ] is the number of measurements ofS + i S − j during the construction of loops. More complicated operators can also be measured within the same framework. Detailed introduction on measuring off-diagonal quantities can be found in Ref. [19] and [2]. 33 Chapter 3 Percolative Quantum Phase Transition in Two-Dimensional Bond-Diluted Antiferromagnets 3.1 Introduction Quantum phase transitions in low-dimensional quantum Heisenberg antiferromagnets (QHAFs) have been the subject of extensive investigations during the last two decades [10, 79]. The transition from a renormalized classical to a quantum-disordered state can be triggered by various parameters, including lattice dimerization, frustration and applied field. More recently, special interest has focused on phase transitions driven by geometric randomness of the lattice [106, 83], including site and bond disorder. 34 In a disordered classical system the phase transition driven by geometric randomness of the lattice belong to the class of percolation transition. When we consider the phase transition at zero temperature in a quantum system defined on the same lattice, the geo- metric randomness is always a relevant factor. However, in a quantum mechanical sys- tem the interplay between geometric randomness and the quantum fluctuations may lead to phase transitions of different universality class. As illustrated in Fig. 3.1, three differ- ent scenarios may appear depending on the strength of the quantum fluctuations in the system. P P c cl P c M 0 A B C Disordered Ordered Figure 3.1: The evolution of order parameterM as concentration factorP . Three possi- ble scenarios of percolative phase transition in a disordered quantum system are shown. A: the quantum phase transition takes place at the classical percolation thresholdP cl c , and belongs to the same universality class as the classical percolation; B: the quantum phase transition takes place at the classical percolation threshold but with different crit- ical exponents; C: the quantum phase transition takes place at a different concentration than the classical percolation threshold. 35 If quantum fluctuations are weak, the quantum transition takes place at the classical percolation threshold, and belongs to the same universality class as the percolation tran- sition. In this scenario, the quantum phase transition actually takes the classical nature, since weak quantum fluctuations are irrelevant factor. At intermediate quantum fluc- tuations the quantum critical point (QCP) still coincides with the classical percolation threshold, but the critical exponents are affected. For sufficiently strong quantum fluc- tuations, the quantum phase transition completely decouples from the classical one. It is entirely dominated by the quantum fluctuations. In this chapter we show that the above scenarios can be realized within an inhomogeneous bond-diluted quantum anti- ferromagnetic Heisenberg model by continuously tuning the randomness. Particularly we show that a novel quantum disordered phase with singular properties appears for strong enough quantum fluctuations. In the following of this section we first give a brief introduction of the percolation tran- sition for lattice models. Then we review recent studies on the percolative phase tran- sitions in 2D diluted quantum antiferromagnets, and introduce a quantum percolation scenario. 36 3.1.1 Classical Percolation Percolation theory has been applied to a variety of fields such as physics, chemistry, and material science. In its mathematical origination, this theory studies the behavior of connected clusters in a random network. Consider a L×L square lattice in 2D. In the percolation problem the lattice is ran- domly diluted. This can be realized by randomly removing either some of the sites (site percolation) or the bonds (bond percolation) on the lattice. We consider the simplest case, the probability of removing each site (bond) is independent, and equals to 1−p. For any finitep, the random diluted system consists lots of nearest-neighbor connected clusters. It is interesting to ask: at a givenp and in the limit of L → ∞, what is the probabilityP ∞ to find a cluster connecting the bottom of the lattice from its top? The answer is even more interesting: the behavior ofP ∞ indicates the system experiences a phase transition. There is a critical thresholdp c , below whichP ∞ is always zero, but P ∞ = 1 whenp > p c . So beyond the percolation thresholdp c there is always a perco- lating cluster connecting the top and bottom of the lattice. We note that the percolating cluster becomes infinitely large whenL→∞. In a 2D percolation system, it is the only infinite cluster in the system. 37 Just like in the thermal phase transitions, scaling behaviors can be observed in perco- lation. The percolation threshold p c gives the position of the critical point. One may define the distribution functionn s (p) as the number of clusters of size s per site. For larges,n s (p) behaves as n s (p) =s −τ F[(p−p c )s η ], (3.1) whereF(x) is an unknown scaling function, andτ andη are critical exponents. It then follows that [61] P ∞ (p)∝ (p−p c ) β , for|p−p c |≪ 1. (3.2) Here the criticalβ is defined as β = τ−2 η . (3.3) Here we see the similarity between percolation and thermal phase transition. The occu- pation probabilityp plays the role as the temperature, andP ∞ behaves just as an order parameter. We may further define the correlation between two sites as the probability they belong to the same cluster. Then the correlation length ξ, defined as the root- mean-square distance between two occupied sites within the same cluster, diverges atp c according to ξ∝|p−p c | −ν , (3.4) where ν = τ−1 dη . (3.5) 38 In 2D d = 2, and τ = 187/91, η = 36/91. These exponents, which depend only on dimensionality of the system, define a universality class of percolation transition. Especially we get the exponents for the order parameter and correlation length asβ = 5/36 andν = 4/3. The critical exponents also determines the structure of the percolating cluster. Let us consider aL×L lattice. The number of sitesM(L) belonging to the percolating cluster on this lattice is given by M(L) =L d P ∞ (p→p c )∝L d ξ −β/ν . (3.6) Close top c , the correlation length diverges, so thatL becomes the unique length scale on the finite lattice. This implies that we can replaceξ byL in Eq. 3.6, and leads to M(L)∝L D f , (3.7) where D f =d− β ν . (3.8) Recall the above results of critical exponents in 2D, we then getD f = 91/48, which shows that the percolating cluster is not a 2D object but a fractal with fractal dimension D f < 2. Since percolating cluster is a random network, it is a fractal in statistical sense. 39 node blob link dead end Figure 3.2: A sketch of the percolating cluster in the nodes-links-blobs model. The backbone lattice is homogeneous at length scaleL>ξ. A useful model to understand the fractal structure of the percolating cluster is called the nodes-links-blobs model [92, 18]. Consider the percolating cluster as an infinitely large random resistor. It includes two parts: one contributes to the total resistance, the other does not. In the nodes-links-blobs model, we call the first part the backbone of the percolating cluster. The latter part which does not change the total resistance but is just attached to the backbone is called dead end. The backbone forms a lattice with lattice spacingξ(p). The lattice is translational symmetric beyond the length scaleξ(p). Nodes 40 refer to the sites on the backbone lattice. Two adjacent nodes are connected through narrow quasi-1D links and quasi-2D blobs. A typical nodes-links-blobs structure is sketched in Fig. 3.2. Links are the weakest part of a percolating cluster since once a link is broken, so does the whole percolating cluster. 3.1.2 The Quantum Percolation Scenario Unlike the classical system, quantum systems may experience phase transitions even at zero temperature. This transition driven by quantum fluctuations instead of thermal ones is known as QPT [79]. Various factors such as lattice dimerization, frustration and magnetic field may trigger a QPT in a magnetic system. Besides those factors, a recent interest is on the QPTs triggered by geometric disorder. Strong geometric disorder not only breaks translational invariance and perturbs the ground state of the pure system, but it can also destabilize renormalized classical phases with long-range order (LRO) and drive the system to novel disordered phases. Various one-dimensional QHAFs with bond disorder have been found to display uncon- ventional quantum phases[45, 22]. For example, the undimerized QHAF chain is driven into a random-singlet phase [22] with algebraically decaying spin-spin correlations. The 41 low temperature susceptibility diverges as 1/[Tlog 2 T] in this phase, independent of the details of the bond disorder. In contrast, a dimerized chain shows a quantum Griffiths phase[35, 116] beyond a critical disorder strength, with exponentially decaying spin- spin correlations and a non-universal power-law divergent susceptibility. In two dimen- sions, the clean QHAF develops antiferromagnetic (N´ eel) LRO at zero temperature [46]. This introduces the intriguing possibility of a genuine order-disorder transition, driven by lattice randomness, i.e. from the N´ eel phase into one of the above unconventional disordered phases 1 . Site and bond dilution of the square lattice QHAF have been the focus of several recent studies [106, 83, 14, 58, 6], motivated by experiments on antiferromagnetic cuprates doped with nonmagnetic impurities [106, 16, 8]. From a geometric point of view, bond and site dilution reduce the connectivity of the lattice, ultimately leading to a percola- tive phase transition [95] beyond which the system is broken up into finite clusters. In a classical spin system, this percolation transition is coupled to a magnetic transition with the same critical exponents since spontaneous magnetic order cannot survive beyond the percolation threshold. In a quantum spin system, on the other hand, a progressive reduc- tion of the lattice connectivity enhances quantum fluctuations in a continuous fashion, raising the possibility of quantum destruction of magnetic order before the percolation threshold is reached. However, recent studies of homogeneously site- and bond-diluted 1 A direct transition from an ordered phase into an unconventional quantum disordered phase driven by randomness has been demonstrated so far only in the case of random quantum Ising models in a transverse field, either in one or two dimensions. See, e.g., Refs. [21, 23], [88], [71]. 42 QHAFs on the square lattice [106, 83, 39, 101] found that the magnetic transition takes place exactly at the percolation threshold, due to the fact that the percolating cluster at threshold shows LRO [83]. The critical exponents of the correlation length and of the order parameter,ν andβ, are found to take their classical percolation values[117, 83]. Therefore, in this case the magnetic transition is completely dominated by classical per- colation 2 . Alternative to the above classical percolation picture is the quantum percolation mech- anism, recently demonstrated in a number of model systems [120, 84, 105, 93]. This scenario is based on the fact that spins involved in locally strongly fluctuating quan- tum states, such as dimer singlets and resonating valence bonds (RVBs), are weakly correlated with the remainder of the system. In a random network of spins, the local strongly fluctuating states create weak links with small spin-spin correlations. If these weak links are part of the backbone of the percolating cluster, they can prevent the per- colating cluster from developing long-range order. We can thus roughly consider the spins belongs to the weak links missing from the geometric cluster they belong to, as sketched in Fig. 3.3. This implies that a cluster may be effectively fragmented by the 2 A recent study (T. V ojta and J. Schmalian, Phys. Rev. Lett. 95, 237206 (2005)) shows that only the critical exponents β and ν coincide with those (β cl = 5/36 and ν cl = 4/3) of the classical universality class, while the other exponents are shifted (α = α cl − νz, γ = γ cl + νz, δ = δ cl + νz/β, and η = η cl − z) due to the appearence of the dynamical critical exponent in the scaling relations (see, e.g., M. Continentino, Quantum scaling in Many-Body Systems, World Scientific, Singapore, 2001). Nonetheless it is interesting to notice that the experimentally accessible exponents, through e.g. neutron scattering, are the unshifted ones, β and ν, along with the critical exponent of the static structure factor γ S = γ− νz which concides with the percolation exponentγ cl . 43 Figure 3.3: A segment of ladder in the percolating cluster leads to local RVB states (yellow ellipsis). The spins involved in the RVB state are very weakly correlated to the rest of the cluster, causing an effective fragmentation of the percolating cluster into non-percolating subclusters. formation of RVB/singlet states. Therefore, if lattice dilution favors the local formation of such states such that the effective fragmentation holds for the percolating cluster, it is possible to drive the system towards a quantum disordered state before the percolation threshold is reached, thus decoupling percolation from magnetic ordering. Moreover, according to this argument the magnetic quantum transition appears as a renormalized percolative transition, with an effective geometry of the clusters dictated by quantum effects. This scenario is referred to as quantum percolation. 44 In this work we numerically scrutinize this quantum percolation scenario in theS = 1/2 QHAF on a square lattice with inhomogeneous bond dilution 3 . The inhomogeneous character of lattice randomness is the key ingredient to enhance quantum fluctuations, favoring the local formation of RVB/dimer-singlet states. Based on large-scale QMC simulations, we observe that the classical percolation scenario is preserved for moderate inhomogeneity, up to a multicritical point beyond which the magnetic transition deviates from the percolation threshold. The critical properties along the critical curve are inves- tigated through extracting the critical correlation length exponentν and the dynamical exponentz. The determined critical exponents defines a new universality class of per- colative quantum phase transition. At the multicritical point, the quantum percolation scenario sets in, leading to a novel quantum-disordered phase. Such a phase has very unconventional features, i.e. a finite correlation length but a gapless excitation spectrum, and a diverging uniform susceptibility at zero temperature with a non-universal expo- nent. These signatures allow us to identify this phase with a genuine two-dimensional quantum Griffiths phase [30, 52, 35, 116]. The rest of this chapter is organized as follows. In Sec. 3.2 the model is described, and in Sec. 3.3, some technical aspects of the simulations are reviewed. In Sec. 3.4 we discuss in detail the fate of antiferromagnetic order on the percolating cluster upon tuning the inhomogeneity of bond disorder. The complete phase diagram of the model and the critical properties at phase transitions are discussed in Sec. 3.5. Sec. 3.6 deals with the 3 A short, partial account of this work has appeared previously [120]. 45 emergence of the quantum-disordered phase in the limit of randomly coupled ladders, whereas Sec. 3.7 is dedicated to the novel ground state properties of the quantum- disordered regime. Conclusions are drawn in Sec. 3.8. 3.2 Inhomogeneous Bond Dilution of the QHAF on the Square Lattice We investigate the ground state and thermodynamics of theS = 1/2 QHAF on a two- dimensional square lattice with inhomogeneous bond dilution. The Hamiltonian of this system is given by H =J X hiji∈D ǫ (ij) D S i ·S j +J X hlmi∈L ǫ (lm) L S l ·S m . (3.9) The sums run over the dimer (D) and ladder (inter-dimer) (L) bonds, as indicated in Fig. 3.4(a). The exchange couplings are taken to be equal and antiferromagnetic (J > 0). ǫ α (α = D,L) is a random variable drawn from a bimodal distribution taking values 1 (‘on’) or 0 (‘off’). The inhomogeneity of this model is of purely statistical nature, stemming from the different probabilities of assigning the ‘on’-state (ǫ i = 1) to the dimer bonds versus the ladder bonds. ForD bonds this probability is denoted byp(ǫ D = 1) =P D , whereas forL bonds it is denoted byp(ǫ L = 1) =P L . The special caseP D = 46 D L (a) (b) (c) P =1, P <1 D P <1, P =1 L D L Figure 3.4: (a) Decomposition of the square lattice into dimer (D) and ladder (L) bonds; (b) randomly coupled dimer limit; (c) randomly coupled ladder limit. P L represents the previously studied homogeneous bond-dilution problem [95, 83]. In the general inhomogeneous caseP D 6= P L , two limits are noteworthy (see Fig. 3.4(b)- (c)): the caseP D = 1,P L ∈ [0,1] corresponds to a system of randomly coupled dimers, and the opposite case P L = 1, P D ∈ [0,1] corresponds to randomly coupled ladders. The degree of inhomogeneity can be parametrized by introducing a variable θ = arctan 1−P L 1−P D (3.10) such that the limiting cases of randomly coupled ladders, homogeneous bond dilution and randomly coupled dimers correspond toθ = 0,π/4, andπ/2 respectively. 47 (b) (a) Figure 3.5: Inhomogeneous bond dilution favors the formation of local strongly fluc- tuating quantum states. The green ellipses indicate such local dimer-singlet, plaquette- singlet and ladder-like RVB states: (a)P L >P D ; (b)P D >P L . It is evident that the inhomogeneous nature of bond dilution in the lattice enhances quantum fluctuations in the magnetic Hamiltonian. In the strongly inhomogeneous lim- itsθ→ 0,π/2 the structure of the percolating cluster is geometrically built from weakly coordinated segments of ladders (θ → 0), and weakly coordinated dimers (θ → π/2), e.g., finite necklace-like structures and finite randomly decorated chains (see Fig. 3.5). Both, antiferromagnetic Heisenberg ladders[17] and decorated chains [36, 57, 90], have quantum-disordered ground states with a finite correlation lengthξ 0 . In antiferromag- netic ladders the ground state has a RVB nature [112], and for the necklace it has a dimer-singlet nature [36, 57]. When segments of such structures become part of a larger cluster, as it is the case in the inhomogeneous percolation model, they locally retain ground state properties similar to the thermodynamic limit if their lengthl is large com- pared to the correlation length ξ 0 . This represents the core mechanism of non-trivial enhancement of quantum fluctuations through bond dilution, leading to novel quantum disordered phases. 48 3.3 Details on Numerical Methods The purely geometric problem of inhomogeneous percolation is studied using a gener- alized version of a highly efficient classical Monte Carlo algorithm [62, 63]. In turn, the quantum magnetic Hamiltonian Eq. (3.10) is investigated by Stochastic Series Expan- sion (SSE) QMC simulations based on the directed-loop algorithm [99]. The ground state properties are systematically probed by efficiently cooling the system down to its physicalT = 0 behavior via a successive doubling of the inverse temperatureβ (Ref. [83]). This approach is necessary since in two-dimensional (2D) diluted systems 2D correlations are often mediated by narrow quasi-one-dimensional links, such that the temperature scale at which these correlations set in is much lower than in the clean 2D case. The scaling analysis of the quantum simulation results is carried out in two comple- mentary ways [83]. In the first approach, the lattice size is fixed at L×L, and the entire system is simulated, including percolating and non-percolating finite clusters. In the second approach, a fundamental theoretical tool is used to ascertain the presence or absence of order in the system: we focus exclusively on the percolating cluster, i.e. discarding the finite clusters that cannot carry long-range order. On a finiteL×L lat- tice with periodic boundary conditions the percolating cluster is identified as the largest cluster, winding around at least one of the two spatial dimensions. The advantage of this 49 approach is that a single run of QMC simulation evaluates the observables on both the full lattice and the largest (percolating) cluster. At the same time the largest cluster has non-negligible size fluctuations (of orderL), leading to a significant error in the QMC data. An alternative approach, eliminating these size fluctuations, is to grow the perco- lating cluster freely from an initial seed up to a fixed sizeN c but without any restriction of lattice boundaries. Starting from a single occupied bond, one visits the six neigh- boring bonds, activating them with probability P D (P L ) if they are D (L) bonds, etc. The procedure is repeated until no new bonds are generated. The cluster geometry is accepted only if its final number of sites matches exactly the requiredN c . A drawback of this procedure is that the cluster build-up becomes time consuming for largeN c and away from the percolation threshold because of the large rejection rate. In the full-lattice studies,L values up to 64 were used. Within theβ-doubling scheme, inverse temperatures as high asβ = 256L have proven to be necessary to observe the physicalT = 0 behavior. For the studies on fixed-N c clusters we consideredN c values up to 2048 and inverse temperatures as high as β = 8N c . For each lattice size and for each pointP =(P D ,P L ),10 2 -10 3 realizations of the diluted lattice/percolating cluster were generated independently to obtain disorder-averaged observables. 50 3.4 Magnetism on the Percolating Cluster In this section, QMC results are discussed which address the evolution of antiferro- magnetic order on percolating clusters upon tuning the nature of bond dilution from homogeneous (θ = π/4) to inhomogeneous (θ → 0,π/2). A fundamental observation [83] is that the presence or absence of magnetic order in a diluted system and the nature of the transition from order to disorder are determined by the magnetic behavior of the percolating cluster. In two dimensions, only one percolating cluster can exist in the sys- tem. We denote withhM c i the disorder-averaged (h...i) staggered magnetization per site of a percolating cluster, which is estimated by hM c i = p hm 2 c i = v u u t * 3 N 2 c X ij (c) (−1) i+j S z i S z j + . (3.11) Here the summation P (c) is restricted to sites contained in the percolating cluster, and N c is the total number of these sites. This is to be contrasted with the overall magneti- zation of the diluted lattice withN =L 2 sites, hMi = p hm 2 i = v u u t * 3 N 2 X ij (−1) i+j S z i S z j + . (3.12) 51 Given that only the percolating cluster contributes to long-range order in the thermody- namic limit, one finds that forN →∞ hm 2 i = * 3 N 2 X ij (c) (−1) i+j S z i S z j + = N 2 c N 2 m 2 c . (3.13) Moreover, exploiting the self-averaging property of the size distribution of the percolat- ing cluster [83], close to the percolation threshold one can use hM(P)i≈A|P −P (cl) c | β M c (P), (3.14) whereP = (P D ,P L ) is the control parameter of the inhomogeneous percolation prob- lem, andP = P (cl) c defines the classical critical line of percolation thresholds in the (P D ,P L ) space. In Eq. (3.14) we have used the critical behavior of the so-called net- work strengthhN c /Ni≈A|P−P (cl) c | β , whereA is a constant amplitude, andβ = 5/36 is the critical exponent of classical percolation [95]. From Eq. (3.14) it is evident that, ifM c (P (cl) c )6= 0, the dilution-driven transition from magnetic order to disorder coincides with the classical percolation transition, both in terms of its location and of its critical exponents. Therefore only the absence of anti- ferromagnetic order on the percolating cluster can lead to a decoupling between the magnetic transition and the geometric percolation threshold. 52 -1.5 -1 -0.5 0 0.5 1 1.5 (P-P (cl) ) L 1/ν 0 0.5 1 1.5 2 (N c /N) L β/ν -1.5 -1 -0.5 0 0.5 1 1.5 0 0.5 1 1.5 2 L=48 L=64 L=96 L=128 L=64 L=96 L=128 L=160 0 0.5 1 P D 0 0.5 1 P L 0 0.5 1 0 0.5 1 finite clusters infinite cluster θ P D =0.9 P L =0.9 c Figure 3.6: Upper panel: classical phase diagram for the inhomogeneous bond- percolation problem on the square lattice. The angleθ parametrizing the critical curve is indicated. Lower panel: scaling plots for the network strengthN c /N in the two strongly inhomogeneous casesP D = 0.9 andP L = 0.9. An excellent data collapse is realized, using the 2D percolation exponentsβ = 5/36 andν = 4/3, and withP (cl) L,c = 0.4007 (P D = 0.9) andP (cl) D,c = 0.0107 (P L = 0.9). The critical curveP = P (cl) c , determined from classical Monte Carlo simulations, is shown in Fig. 3.6. Despite the inhomogeneous nature of percolation for P L 6= P D , excellent scaling of the classical MC results is observed using the 2D percolation expo- nents [95]. This demonstrates that the universality class remains unchanged along the critical curve. 53 Starting from the classical Monte Carlo result, the evolution of the order parameter can be tracked on the percolating clusterhM c (P (cl) c )i along the critical lineP = P (cl) c to monitor the effect of quantum fluctuations enhanced by inhomogeneity. To this end, we perform a scaling study of the staggered magnetization on percolating clusters for fixed sizeN c , withN c up to 2048 andβ = 8N c . The obtained data are extrapolated to the thermodynamic limit by a three-parameter polynomial fithm 2 c (N c )i =hm 2 c (∞)i+ aN −1/2 c +bN c . As presented in Fig. 3.7, at the homogeneous bond percolation thresholdP D = P L = 0.5 the percolating cluster is found to have antiferromagnetic long-range order, in agree- ment with Ref. [83]. However, the finite value of the order parameter is strongly reduced as the inhomogeneity is turned on, both in the ”dimer” (θ → π/2) and the ”ladder” (θ → 0) directions, until it vanishes at non-trivial θ values, θ L ≈ 0.34 in the ladder limit and θ D ≈ 1.23 in the dimer limit. This means that inhomogeneous percolation exhibits non-linear quantum fluctuations that are able to destroy the antiferromagnetic long-range order, in contrast to the homogeneous percolation case[83] with a renormal- ized classical ground state. According to Eq. (3.14) this must be reflected in a profound change of the critical properties, since the vanishing of the staggered magnetization on the percolating cluster decouples the magnetic transition from the percolation threshold. Beyond the two crit- ical valuesθ L andθ D , a higher concentration of bonds will be required for the system 54 0 0.05 0.1 N c -1/2 0 0.02 0.04 0.06 θ=0.92 θ=1.08 θ=1.24 0 0.05 0.1 N c -1/2 0 0.02 0.04 0.06 <m c 2 > θ=0.33 θ=0.41 θ=0.50 θ=0.79 0.5 1 1.5 θ 0 0.005 0.01 <m c 2 > Figure 3.7: Upper panel: scaling of the disorder-averaged squared staggered magneti- zation on the percolating cluster with fixed sizeN c . The different curves correspond to various points along the classical percolation transition. The continuous lines represent quadratic fits. Lower panel: extrapolated thermodynamic values of the staggered mag- netization from the upper panel as a function of the inhomogeneity parameter θ. The dashed line is a guide to the eye. to magnetically order than for the lattice to percolate, thus opening up an intermediate phase with a novel quantum-disordered ground state. 3.5 Magnetic Phase Diagram and Critical Properties In this section, QMC results are discussed for the magnetic phase diagram in the (P D ,P L ) plane. The magnetic transitions are located via scaling of the spin-spin correlation length ξ, extracted from the structure factor through the second moment estimator [15]. The 55 inhomogeneity of bond dilution breaks the discrete rotation symmetry of the lattice, such that two distinct correlation lengthsξ x andξ y need to be considered along thex andy lattice directions. At critical points separating 2D long-range order from disordered phases, both correlation lengths must scale linearly with the lattice size, ξ x ,ξ y ∼ L, regardless of the universality class of the transition. In locating the magnetic transitions we have verified that this condition is satisfied. This leads to the conclusion that, despite its inhomogeneity, the magnetic system has unique well-defined phase transitions. The resulting phase diagram is shown in Fig. 3.8. The magnetic transition line coincides with the percolation line for moderate inhomogeneity. As discussed in the previous section, in this region of the phase diagram the percolating cluster at the percolation threshold is antiferromagnetically ordered. Therefore not only the magnetic transition and the geometric transition coincide, but also the critical exponents of the magnetic transition are those of 2D percolation [120], according to Eq. (3.14). For strong inhomogeneity. the scenario of a classical percolation transition in the mag- netic Hamiltonian is precluded by non-linear quantum fluctuations. In fact, two multi- critical points occur, in the ladder and in the dimer direction, beyond which the magnetic transition decouples from the percolation threshold. The location of these multicritical points is quantitatively consistent with the vanishing of antiferromagnetic order on the percolating cluster at two critical values of the inhomogeneity parameterθ (Fig. 3.7). 56 0 0.2 0.4 0.6 0.8 1 P D 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 P L classical quantum (S=1/2) Finite clusters <M> = 0 Infinite cluster Infinite cluster <M> = 0 <M> = 0 Figure 3.8: Phase diagram of the spin-1/2 QHAF on the inhomogeneously bond-diluted square lattice. Beyond the two multicritical points, intermediate quantum phases appear, in which anti- ferromagnetic order is absent on the percolating cluster both at and away from the per- colation threshold. These regions represent novel quantum-disordered phases on an infinite family of 2D percolated random lattices. Representative structures of this fam- ily are shown in Fig. 3.9, together with the homogeneous percolating cluster. It is remarkable that the homogeneous percolating cluster (Fig. 3.9(b)) has long-range anti- ferromagnetic order in the thermodynamic limit, whereas the other percolating clusters in Fig. 3.9(a),(c),(d) do not. In fact, the latter structures have a significantly higher average number of bonds per site than the homogeneous percolating cluster, which has a fractal dimension[95] lower than 2 (D f = 91/48). A detailed discussion of the nature of the quantum-disordered phase will be given in Sec. 3.7. 57 1 11 21 E 1 11 21 1 11 21 C 1 11 21 (a) (b) (c) (d) Figure 3.9: Percolating clusters for the bond-diluted square lattice with size 24× 24 and periodic boundary conditions: (a)P D = 0.15,P L = 0.8; (b)P D = 0.5,P L = 0.5 (homogeneous case); (c)P D = 0.9, P L = 0.45; (d)P D = 1.0, P L = 0.44 (randomly coupled dimers). Now we investigate the critical exponents associated with this magnetic phase transition. Main results are reported in Ref. [120]. The finite-temperature quantum critical behavior of the correlation lengthξ∼T −1/z determines the dynamical critical exponentz [79, ?], whereas the critical exponentν governs the divergence of the correlation length at the T = 0 critical point,ξ∼|P −P c | −ν , whereP = (P D ,P L ).P c can be obtained as the value of disorder at whichξ exhibits a power-law behavior as a function ofT , and it is in agreement with the estimates coming from theT = 0 analysis described above. From the scaling relation ξ(P;T)∼ 1 T 1/z F ξ (T |P −P c | νz ) (3.15) 58 10 0 10 1 T |P’−P’ c | −νz 10 −1 10 0 ξ T 1/z P’=0.50 0.52 0.53 0.54 0.56 0.57 0.58 0.60 0.62 y x y x A 0 0.5 1 1.5 θ (rad) 1 1.5 2 z ν ν=1 ν=4/3 z=1.9 classical quantum quantum B Figure 3.10: (A) Scaling plot of ξ x (open symbols), ξ y (full symbols) for P D = 1 (randomly coupled dimers). (B) Critical exponentsz andν for the magnetic transition in the inhomogeneously bond-diluted S=1/2 Heisenberg antiferromagnet. The angleθ parametrizes the critical line in Fig. 3.8. Open (full) symbols refer to estimates obtained through the scaling ofξ x (ξ y ). The shaded areas correspond to the regions of parameters in which the transition is quantum. All lines are guides to the eye. the exponentν can be determined as the one realizing the best collapse ofT 1/z ξ(P;T) curves plotted versus T|P −P c | νz for differentP values, as shown in Fig. 3.10(A). The values extracted forν andz are shown in Fig. 3.10(B) as a function of the angleθ which parametrizes the critical curve of Fig. 3.8. Theν exponent is found to satisfy the Harris criterion for disordered systems [12]ν ≥ 2/d = 1 within error bars. Moreover, we observe that z and ν are almost constant over a large portion of the critical line, taking the valuesz≈ 1.9 andν≈ 4/3, and they both show a drop toz,ν≈ 1 only near the edges of the critical line in the strongly quantum regimes (P D → 1 andP L → 1). In particular, over a large portion of the critical line the correlation length exponentν is clearly consistent with the classical valueν = 4/3 for two-dimensional percolation [95]. Thez exponent can also be related to the critical exponents of the classical phase transition. At and around the homogeneous percolation pointP =P ′ = 0.5, where the 59 transition is classical, the exponentz is the scaling dimension of the finite-size gap of the percolating cluster diverging in the thermodynamic limit, Δ∼ L −z . Since the number of sitesN c of this cluster diverges asN c ∼ L D withD = 91/48 [95], if the finite-size gap vanishes as Δ∼N −ρ c it is easy to see thatz =ρD. For clusters with the geometry of the chain or of the square lattice,ρ = 1 [112, 7]. Extending this result to the fractal percolating cluster we obtainz=D=1.89.., in agreement with the QMC result. Surprisingly, the two multicritical points do not mark an evident discontinuity in the universality class of the model, and the classical percolation exponents appear to persist also for intermediate inhomogeneity beyond the multicritical points, changing to radi- cally different values only in the extreme inhomogeneous limits. Therefore the magnetic quantum phase transitions close to the multicritical points appear to retain a percolative character. This result can be understood within the following argument. To destroy the long-range nature of a random network it is sufficient to cut a few links on its backbone. This implies that the non-linear quantum fluctuations leading to quantum-disorder on the percolating cluster need only have short-wavelength components, in contrast to con- ventional quantum phase transitions in translationally invariant lattices. Indeed, strongly fluctuating states appearing on links made of segmented ladders or decorated chains (see Sec. 3.2) have a markedly local effect on magnetic correlations, i.e. they magnetically decouple the portions of the percolating cluster connected through that link. This effect is equivalent to geometrically removing that link, thus leaving two portions of the clus- ter disconnected [120]. If weak links occur on the backbone of the percolating cluster, 60 this mechanism is sufficient to lead to quantum disorder. Only quantum fluctuations with wavelengths comparable to the length of the weak link are required. To reestablish magnetic correlations between the two portions of the percolating cluster disconnected by quantum fluctuations, it is then necessary to add more bonds to the cluster in order to find an alternative path for magnetic correlations to spread over the cluster. Since this process is of geometrical character it endows the quantum phase transition with a percolative nature. The above picture of a percolative quantum phase transition breaks down for strong inhomogeneity. In this limit, the building blocks of the percolating cluster are the strongly quantum fluctuating substructures depicted in Fig. 3.5, and hence the quan- tum phase transition is driven by the competition between the energy scale of the gap above the ground state of such substructures and the energy scale of the spatially random couplings between the substructures. In a real-space picture, the competition occurs between the finite correlation length of the substructures and the characteristic length set by the spacing between the spatially random couplings. This competition will be discussed in detail in the next section. The quantum phase transition is completely dis- connected from classical percolation of the lattice, and it is reasonable to expect for it to be in a different universality class, as concluded in Ref. [120]. 61 3.6 Criticality Near the Ladder Limit So far, the onset of the quantum disordered phases was discussed, approaching them from the homogeneous limit of Eq. (3.9). A complementary understanding of the novel quantum phase and the associated quantum phase transitions can be gained by starting from the extreme inhomogeneous limit of randomly coupled ladders (θ = 0). This limit lends itself to a very simple analysis in terms of the properties of the single spin-1/2 antiferromagnetic Heisenberg ladder. In the limit ofP L = 1, the system is composed of an array of two-leg ladders randomly coupled byD bonds. The probability of activation of a dimer bond is reflected in the average distancehri of twoD bonds along a column separating two neighboring ladders, namelyP D = 1/hri. A singleD bond is sufficient to connect two ladders geometrically. In other words,hri =L is sufficient to connect all ladder subsystems on aL×L lattice. Taking the limitL → ∞, one obtains immediatelyP D → 0. Thus geometrically the system percolates at any infinitesimal concentrationP D > 0 and the percolating cluster is the entire lattice itself. However the system does not develop antiferromagnetic order until a finite value P D ≈ 0.32 is reached. This is not surprising since each two-leg ladder has a quantum- disordered ground state with a finite excitation gap. This gap has to be overcome before 62 <r>~ ξ ladder Figure 3.11: Cartoon of the ordering mechanism in randomly coupled ladders. When P D = 1/hri≈ 1/ξ ladder (see text) percolating strings of block spins (rectangles) appear in the system. The divergence of the correlation length along such strings drives the onset of long-range correlations between the strings and ultimately of 2D long-range order. the N´ eel ordered state can become the ground state. A more quantitative argument can be formulated based on the fact that each isolated ladder has a finite correlation length [112]ξ ladder ≈ 3.19. One can imagine a renormalization group transformation that cre- ates effective block spins from all the spins within a correlation volume in such a way that all block spins are perfectly uncorrelated with each other (see Fig. 3.11). When turning on the D bonds, it becomes evident that, to establish long-range correlations in the system, it is necessary to have a string of block spins percolating from one side to the other of the system. Given that all block spins are independent within the same ladder, this is only possible if, on average, each block spin is connected by aD bond to both its right and left neighboring block spin. This requires an average spacing of theD 63 bonds ofhri≈ ξ ladder on each column, which in turn leads to the following estimate of the critical concentration of dimer bonds P D,c = 1 ξ ladder . (3.16) Despite the simplicity of the argument, this estimate is surprisingly good:P D,c ≈ 0.313 to be compared with the QMC resultP D,c = 0.32(1). 0 0.1 0.2 0.3 0.4 0.5 1-P L 0 2 4 6 < l >/2 ξ ladder Figure 3.12: Correlation length of the bond-diluted 128× 2 ladder compared with the average ladder lengthhli. The data forξ ladder are taken atβ = 2048; comparison with data atβ = 4096 (not reported) shows no significant deviation. Interestingly, this argument can be extended to the case of randomly coupled weakly diluted ladders, namely for P L . 1. The correlation length of the ladder remains the 64 dominant length scale of the problem, as long as it is much smaller than the new length scale introduced by the dilution hli = 1 (1−P L ) 2 , (3.17) corresponding to the average length of the ladder segments after dilution. (a) (b) (c) J eff Figure 3.13: Dominant effects of bond dilution on a two-leg ladder: (a) local antiferro- magnetic modulation due to a missing rung bond; (b) effective couplings between spins missing their rung bond; (c) enhancement of the singlet component (indicated by an ellipse) of the two rung dimers adjacent to a missing leg bond. To test to which extent the above argument is applicable we performed SSE-QMC sim- ulations on a single bond-diluted two-leg ladder to determine the evolution of the cor- relation length upon dilution. The results are shown in Fig. 3.12. For low doping concentrations, we observe that the average correlation length of the ladder increases moderately. For higher concentrations of missing bonds, the average ladder lengthhli 65 becomes comparable to the correlation length, which then crosses over to a decreasing behavior simply reflecting that ofhli. We now focus on the increasing tendency for small dilution. This behavior reflects various competing mechanisms. The dominant effects having the highest probabilities (O(P L ) andO(P 2 L )) are: (a) The removal of a rung bond results locally in two uncoupledS = 1/2 free moments which introduce a local staggered modulation of the adjacent spins [25, 55] within a correlation volume∼ ξ ladder (Fig. 3.13(a)). This leads to a local increase in the antifer- romagnetic correlations along the legs; (b) If two rung bonds are removed sufficiently close to each other, the induced local antiferromagnetic modulations of the ladder ground state “lock in phase”, leading to an effective coupling between the free moments, exponentially decaying with the distance in the limit of low doping [55, 91], J eff ∼ (−1) i−j exp[−|i−j|/ξ ladder ]. (3.18) These long-range interactions also lead to an increase of the correlation length in the system through an order-by-disorder mechanism (Fig. 3.13(b)); (c) The removal of one leg bond, on the contrary, leads to a local enhancement of the rung-singlet component of the state of the two adjacent rungs (Fig. 3.13(c)). This 66 mechanism weakens the correlations along the two legs. Effect (c) has twice the probability compared to the effects (a) and (b) because there are twice as many leg bonds as rung bonds. The non-trivial competition between cor- relation enhancement (a)-(b) and correlation suppression (c) has the combined effect that in bond-diluted ladders, the correlation length always remains finite down to zero temperature. This is to be contrasted with the case of site-diluted ladders, in which the correlation length diverges algebraically with decreasing temperature, ξ ∼ T −0.4 , as reported in Ref. [29]. (Obviously, even in the site-diluted case the average correlation length is eventually upper-bounded by the average length of the ladder segments). It is interesting to note that for the case of site dilution, correlation-suppressing effects are not present beside the simple fragmentation of the ladder. The absence of diverging correlations in bond-diluted ladders is crucial for the occur- rence of an extended quantum-disordered phase. The finite correlation length of a single ladder for P L . 1 allows us to repeat the RG argument sketched in Fig. 3.11 for the case of bond-diluted ladders. One needs a finite and non-classical concentration of dimer bondsP D,c to drive the system into a N´ eel ordered phase, directly related to the 67 0 0.2 0.4 P D 0.5 0.6 0.7 0.8 0.9 1 P L classical quantum (S=1/2) P D,c =1/ξ ladder (P L ) Finite clusters <M> = 0 Infinite cluster Infinite cluster <M> = 0 <M> = 0 (b) Figure 3.14: Phase diagram of the QHAF on the inhomogeneously bond-diluted square lattice (see also Fig. 3.8) close to the ladder limit and compared with the single-ladder predictionP D,c = 1/ξ ladder (P L ) (see text). quantum correlation length of the systemξ ladder , as long as it is significantly lower than the classical lengthhli: P D,c (P L )≈ 1 ξ ladder (P L ) . (3.19) The above estimate is a theoretical prediction for the quantum-critical curve P c = (P D,c ,P L,c ), which can be directly compared with the QMC results, shown in Fig. 3.14. The range of quantitative validity of Eq. 3.19 is surprisingly large, extending down to P L ≈ 0.875, which is interestingly also the value where the quantum correlation length ξ ladder (P L ) starts to cross over towards the decreasing classical behavior (Fig. 3.12). At this point the ground state of the ladder segments is strongly altered by the presence of the dimer bonds, and hence the argument leading to Eq. (3.19) breaks down. 68 3.7 The Quantum-Disordered Phase: Correlations and Griffiths-McCoy Singularities In this section, we discuss the nature of the quantum-disordered phase. Ground-state properties and thermodynamic observables are considered which reveal the nature of the low-lying excitation spectrum. A fundamental aspect of this phase, as it is typical for disordered systems, are large variations in the local properties of the system, reflecting the local geometric structure defined by disorder. In the previous sections, we have discussed how the presence of segments of ladders or decorated chains leads to the local formation of RVB/dimer-singlet states. At the same time, bond dilution of ladders leads to the appearence of S = 1/2 local degrees of freedom, effectively interacting across the strongly quantum fluctuating regions. The coexistence of such different subsets of spins within the same system raises the question of what the global nature of the ground state really is. 3.7.1 Short-Range Correlations We first characterize the ground state in terms of global correlation properties. The absence of magnetic order in this phase has already been discussed in Sec. 3.4. Fig. 69 10 0 10 1 10 2 10 3 10 4 10 5 10 6 β 2.1 2.2 2.3 2.4 2.5 S(π,π) L=16 L=24 L=32 L=40 L=48 L=56 L=64 L=72 0 20 40 60 80 L 2.2 2.4 2.6 2.8 3 S(π,π) T=0 Figure 3.15: Temperature and finite-size scaling of the static structure factor S(π,π) for a representative point in the quantum-disordered regime (P D = 0.15,P L = 0.85). Dashed lines are a quartic polynomial fit in T = β −1 [see Eq. (3.20)] to extrapolate to theβ → ∞ limit (when required). Inset: β → ∞ extrapolated values of the static structure factor plotted vs. system size, clearly showing saturation in the thermodynamic limit. 3.15 shows the staggered structure factorS(π,π) for a representative point (P D = 0.15, P D = 0.85) in the quantum disordered phase on the ”ladder” side of the phase diagram in Fig. 3.8. The finite-temperature data for S(π,π), obtained using β-doubling, are extrapolated to the limitT → 0 by accounting for power-law temperature corrections up to fourth order: S(π,π;T) =S(π,π) T=0 + 4 X i=1 a i β i (3.20) where the a i ’s are fitting coefficients. The extrapolated T = 0 values of S(π,π) are shown in the inset of Fig. 3.15 as a function of system size. They clearly display saturation forL → ∞, proving not only that the system is disordered, but also that it has finite-range correlations. 70 1 2 3 4 5 6 ξ x 0 50 100 L 0 5 10 ξ T=0 10 2 10 3 10 4 10 5 β 2 4 6 8 ξ y L=16 L=24 L=32 L=40 L=48 L=56 L=64 L=72 ξ y ξ x Figure 3.16: Left panel: temperature and finite-size scaling of the correlation length along the two lattice directions (ξ x and ξ y ) for a representative point in the quantum- disordered regime (P D = 0.15,P L = 0.85). Dashed lines are a quartic polynomial fit in T =β −1 to extrapolate toβ →∞ (when required). Right panel: β →∞ extrapolated values of the correlation length(s) vs. system size; the dashed lines are polynomial fits up to second order in L −1 on the points for L ≥ 48 to compensate for the finite-size effects. This fact is further reflected in the correlation lengthξ x(y) along the two lattice dimen- sions, shown in Fig. 3.16. A fitting procedure similar to Eq. (3.20) is used to eliminate polynomial finite-temperature correction. The ξ T=0 values so obtained (right panel of Fig. 3.16) suggest convergence towards a finite value forL→ ∞, although the corre- lation length is large for the particular point (P D = 0.15, P D = 0.85), and, even for the largest considered size (L = 72),ξ does not display full saturation. A polynomial fit ξ(L) = ξ ∞ +a/L +b/L 2 for L ≥ 48 yields saturation values ξ x,∞ = 9(1) and ξ y,∞ = 12.5(5), resulting in a considerable correlation volume of ∼ 100 sites. It is important to keep in mind that the finite correlation length in the quantum disordered 71 0 10 20 30 40 50 60 70 80 90 L 1 2 3 4 5 6 S(π,π) ξ x ξ y Figure 3.17: Finite-size scaling of the T = 0 static structure factor and correlation length(s) for a point in the quantum-disordered regime (P D = 0.15, P L = 0.9) closer to the ladder limit. The dashed lines associated with the correlation length data are quadratic fitsξ(L) =ξ(∞)+a/L+b/L 2 . system varies strongly with the degree of inhomogeneity, and it can become as small as ξ∼ 3 in the limit of randomly coupled ladders. For instance, Fig. 3.17 shows the finite- size scaling of the static structure factor and correlation lengths for a point closer to the ladder limit,P D = 0.15 andP L = 0.9. For this point one observes excellent temperature saturation of the data up toL = 64 forβ ≤ 16384 without the need of any fitting pro- cedure. Moreover, the zero-temperature data reported in Fig. 3.17 show a much better saturating behavior within the system sizes considered, leading to the asymptotic values ξ x,∞ = 2.5(1) andξ y,∞ = 5.4(2). A clear picture of the short-range nature of correlations also emerges from real-space images of 72 (a) E=-J/2 E=-J/4 S(i)=2.5 (b) E=-J/2 E=-J/4 C(0,0)=0.25 C(0,i)=0.01 Figure 3.18: Real-space images of the bond energy (both panels), local static staggered structure factorS(i) (see text) (a) and of the static staggered correlationsC(0,i) between a reference point on a dangling bond and the remainder of the system (b), for the perco- lating cluster on a 32× 32 sample withP D = 0.15 andP L = 0.85 atβ = 8192. The thickness of each bond is proportional to its energy, whereas the radius of the dots in (a) and (b) is proportional to the local value ofS(i) andC(0,i). In (a) red boundaries highlight the regions where local correlations exceed the average value, namely where S(i)& S(π,π). In (b) the red boundaries mark the correlation volume beyond which C(0,i) ≤ 8× 10 −3 , corresponding to a distance (in units of the correlation length) |r i |/ξ≈ ln[C(0,0)/C(0,i)]≈ 3.4 . • the local staggered structure factor, S(π,π;i) = X j6=i (−1) i+j hS z i S z j i 0 ; (3.21) • the correlation function C(0,i) = hS z 0 S z i i 0 between a reference site and all the other spins; • the energy of each bondE hiji =JhS i ·S j i 0 , 73 all shown in Fig. 3.18 for a given sample corresponding to a representative point of the quantum-disordered phase on the “ladder side”. Here the symbolh...i 0 indicates the expectation value associated with the ground state of the specific sample considered, not to be confused with the disorder averageh...i. Starting with the bond energy, we observe that quantum fluctuations introduce a wide range of variability in its values, energetically promoting weakly correlated dimers or quadrumers and rung bonds on ladder segments, at the expenses of the energy of the adjacent bonds (due to the fact that aS = 1/2 particle can form a singlet with only one otherS = 1/2 particle). In particular the rareD bonds have in general the weakest ener- gies, due to the strong tendency of the energy to be minimized on the ladder segments. This leads to a clear quantum suppression of correlations in the direction transverse to the ladders. The local staggered structure factor of Eq. (3.21) measures the “effective correlation volume” (integral of the antiferromagnetic correlation function) around each spin, thus identifying the spins that are most strongly correlated with the remainder of the system. In Fig. 3.18(a) we observe that only a portion of the spins has a sizable local structure factor, when they belong to regions with a larger local coordination number. In Fig. 3.18(b) we have picked a site on a dangling bond attached to a segment of a ladder, and calculated its correlations to the rest of the system. Compared with Fig. 3.18(a), we observe that the correlation volume around the chosen spin is considerably smaller than 74 the cluster of “correlated” spins identified by the local structure factor. This implies that spins that are most correlated with the remainder of the system are not all correlated with each other. The finite-range nature of correlations is also clearly shown in the picture, along with the non-monotonic decay of correlations with distance. The latter is due to the fact that a dangling spin is mostly correlated with spins not involved in a local sin- glet state, even over long distances, due to the effective long-range couplings discussed in Sec. 3.6. Nonetheless, even in the presence of such couplings, correlations decay significantly with distance. This picture is to be contrasted with the case in which long- range couplings give rise to a genuine order-by-disorder phenomenon, as for instance in the site-diluted dimerized systems treated in Ref. [118]. There the real-space image of C(0,i) shows very little spatial decay of the correlations between the free moments of the system. 3.7.2 Griffiths-McCoy Singularities In this section, the nature of the low-lying excitation spectrum in the quantum- disordered phase is discussed. In a clean system, the presence of a finite real-space correlation length implies finite correlations also in imaginary time, i.e. the existence of a gap in the excitation spectrum. In disordered systems, however, this is not neces- sarily the case. The disorder considered here is completely uncorrelated in real space 75 but perfectly correlated in imaginary time, such that long-range correlations in the time dimension (and thus absence of a gap) may coexist with short-range real-space correla- tions. 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 T/J 10 -3 10 -2 10 -1 10 0 10 1 χ u,c P L =0.75 P L =0.8 P L =0.85 P L =0.9 P L =0.95 α=0.199(5) α=0.206(9) α=0.202(6) α=0.175(5) α=0.173(5) P D =0.15 L=64 percolating cluster only 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 T/J 10 -1 10 0 10 1 χ u,c P D =0.85 P D =0.9 P D =0.95 α=0.082(9) α=0.090(8) α=0.110(5) P L =0.42 N c =1024 percolating cluster only Figure 3.19: Low-temperature uniform susceptibility of the percolating cluster in the quantum-disordered regime: ladder regime (upper panel) and dimer regime (lower panel). The solid lines are power-law fits of the form χ u,c ∼ T −1+α , and the result- ing fit coefficientsα are indicated. In the upper panel, the percolating cluster is picked as the largest cluster in a 64×64 lattice, whereas in the lower panel it is grown freely from a seed site up to a fixed site ofN c = 1024 (see Sec. 3.3). The nature of the low-lying excitation spectrum is reflected in the low-T uniform sus- ceptibility. The global uniform susceptibility of the system is typically dominated by the Curie contribution from odd-numbered finite clusters, which masks the behavior of the quantum-disordered percolating cluster. Hence, in the following discussion the 76 contributions from the finite clusters are neglected. Instead, we focus on the uniform susceptibility of the percolating cluster with an even number of sitesN c : χ u,c = β N c X ij∈c hS z i S z j i. (3.22) The results for χ u,c as a function of temperature, both in the ladder regime and in the dimer regime, are shown in Fig. 3.19. A gapped spectrum would imply an exponentially vanishing susceptibility asT → 0. However, we find a non-vanishing susceptibility in this limit, and observe thatχ u,c instead displays a power-law divergence χ u,c ∼T α−1 (3.23) withα6= 0, such that this divergence is not Curie-like. In particular,α = α(P) varies continuously as one scans across the quantum-disordered region:α≈ 0.09÷0.2. This clearly shows that the system is not only gapless, but that its gapless excitations lead to singularities in the response to an external field. Quantum-disordered phases with divergent response properties atT → 0 are familiar in the literature of disordered systems. They are classified as quantum Griffiths phases, and have been theoretically shown to exist in bond-disordered dimerized chains [35, 116], 77 in bond-disordered S = 1 (Haldane) chains [34, 100], and in disordered transverse- field Ising models [21, 23, 88, 71] (see Ref. [37] for a recent review). As in our model, their central feature is the coexistence of short-range correlations together with so-called Griffiths-McCoy singularities in the response functions. Such singularities are due to the anomalous role played by local low-energy excitations living on rare regions of the system. (b) (a) J eff l Figure 3.20: Sketch of a low-energy excitation in the bond-diluted ladder system. (a) Four spins with missing rung bonds are coupled across a clean region of lengthl with weak effective couplingsJ eff ∼ exp[−σl]; (b) rotating the spins at one end of the clean region by the same angle, spin excitations occur on the weaker bonds only, leaving the stronger ones unchanged. What are the excitations leading to the singular response in our model? To answer this question, we focus on the ladder limit, referring to the discussion of bond dilution in ladder systems given in the previous section. There, it was observed that bond dilu- tion leaves freeS = 1/2 moments interacting via effective long-range couplings,J eff , 78 exponentially decaying with the distance. The presence of this weak energy scale in the system immediately populates the low-lying spectrum down to the ground state. For any arbitrarily small energy, one can always find twoS = 1/2 moments which are suf- ficiently far from each other to be coupled with this energy. Nonetheless, to be able to excite the system at this energy scale only, one needs the two moments to be separated by a large clean region, such that one can rotate all moments at one end of this region of a given angle, leaving the higher-energy couplings at shorter range unchanged (see Fig. 3.20). Therefore on the one hand there is an exponentially small energy of the localized excitation, leading to an equally weak local gap Δ: Δ∼ exp[−σl], (3.24) but, at the same time, the presence of a large clean region of lengthl between the free moments is a rare event with exponentially small probability w(l)∼ (1−P L ) l ∼ exp[−cl]. (3.25) By substituting Eq. (3.24) into Eq. (3.25) with the appropriate metric term to preserve the normalization condition of the probability distribution, the two exponentials com- pensate each other, leaving a power-law distribution[75] of the local gapΔ: w(Δ)∼ Δ c/σ−1 . (3.26) 79 -2 0 2 4 6 -10 -8 -6 -4 -2 lnP(ln loc ) ln loc =128L L=16 , =0.77 L=24 , =0.42 L=32 , =0.35 L=48 , =0.30 L=64 , =0.27 P D =0.10, P L =0.90 Figure 3.21: Probability distribution of the local susceptibility for a representative point in the quantum-disordered phase (P D = 0.15,P D = 0.9), taken atβ = 128L. The slope of the linear tail of the distribution is proportional to theα exponent in the thermody- namic limit (see text). A very similar argument can be applied to the dimer limit of the model. Here, e.g., an exponentially rare long clean decorated chain with weakly interacting spins at its ends leads to a similar probability distribution of local gaps. The above result has immediate consequences for the thermodynamics of the system [21, 23], giving rise to the observed power-law behavior of the uniform susceptibility χ ∼ T α−1 where α = c/σ. In particular the constantc is clearly disorder-dependent [c = |ln(1−P L )| via Eq. (3.25)], such thatα is expected to be non-universal, as it is observed in our data. 80 In order to verify the microscopic interpretation of the QMC results for χ u,c in terms of a power-law distribution of the local gaps, Eq. (3.26), we study an observable that directly probes the local gap, i.e. the local susceptibility [75], χ loc,i ≡ Z β 0 dτhS z i (τ)S z i (0)i. (3.27) Hence, if the i-th spin is involved in a local quantum state with a local gap Δ, one obtains in the large-τ (low-T ) limit: hS z i (τ)S z i (0)i∼ exp[−Δτ], (3.28) such that, upon imaginary-time integration χ loc ,i ∼ 1 Δ . (3.29) By a simple change of variables in Eq. (3.26) one then finds a probability distribution for the local susceptibility, w(χ loc )∼χ −1−α loc , (3.30) which is more conveniently recast in the form lnw(lnχ loc )∼−αlnχ loc , (3.31) 81 i.e. the function lnw(lnχ loc ) has a tail linear in (lnχ loc ) and a slope that directly yields the exponentα. Fig. 3.21 shows results for the logarithmic distribution of the local susceptibility, Eq. (3.31), on sites of the percolating cluster only, for a representative point in the quantum-disordered phase (P L = 0.9,P D = 0.1) and for various system sizes at very low temperatures (β = 128L). Consistent with the Griffiths-singularity scenario, one observes a linear tail inlnw(lnχ loc ), reflecting the power-law tail ofw(χ loc ) for all con- sidered system sizes. This tail shows a strong size dependence, with a decreasing slope for increasing size, and crossing over into a non-linear finite-size tail for largeχ loc val- ues. The tail of the distributionw(χ loc ) is the result of rare events associated with clean regions and exponentially small gaps, such that finite-size corrections significantly affect the tail statistics. Nonetheless, the extracted values of the tail slopeα converge towards the rangeα≈ 0.2÷0.3, fully consistent with the independent estimate obtained from the low-temperature behavior of the global susceptibilityχ u,c in the ladder limit. This result completes the picture of a quantum Griffiths phase for the quantum-disordered region of the system. 82 3.8 Summary In this chapter, we have investigated the highly non-trivial interplay between quantum fluctuations and geometric disorder realized in the inhomogeneously bond-diluted quan- tum Heisenberg antiferromagnet on the square lattice. We have observed that the inho- mogeneous nature of bond dilution allows for a continuous tuning of quantum fluctu- ations from linear to non-linear. Inhomogeneity enables us to tune the system from a renormalized classical phase, in which the magnetic transition is purely driven by geo- metric percolation, to a genuine quantum regime in which the magnetic transition is decoupled from the percolation transition. Hence, a novel quantum-disordered phase appears between the magnetic and the percolation transition. This phase has the nature of a quantum Griffiths phase, characterized by finite-range correlations coexisting with a divergent uniform susceptibility. This represents, to our knowledge, the first evidence of a quantum Griffiths phase in disordered 2D Heisenberg antiferromagnets 4 , the previous observations being limited so far to one-dimensional Heisenberg antiferromagnets [35, 116, 34, 100] and one- and 4 A real-space renormalization group study of the dimerized 2D Heisenberg antiferromagnet with bond disorder (Y .-C. Lin, R. M´ elin, H. Rieger, and F. Igl´ oi, Phys. Rev. B 68, 024424 (2003)) suggests the occurrence of a quantum Griffiths phase in the dimerized 2D Heisenberg antiferromagnet with bond disorder. Nonetheless, a detailed quantum Monte Carlo study of the above model in presence of bond dilution (C. Yasuda, S. Todo, M. Matsumoto, and H. Takayama, J. Phys. Chem. Solids 63, 1607 (2002); C. Yasuda, S. Todo, M. Matsumoto, and H. Takayama, Prog. Theor. Phys. Suppl. 145, 339 (2002)) shows that the gapped dimer-singlet phase of the clean system is stable to disorder, and it is driven to a gapless antiferromagnetic phase at a finite value of bond dilution through an order-by-disorder mechanism. 83 two-dimensional quantum Ising models [21, 23, 88, 71, 75]. Given the relevance of the two-dimensional Heisenberg model to the physics of cuprate oxides, we believe that this finding represents an important step towards the possibility of the experimental observation of a quantum Griffiths phase. Obviously, the most natural form of disorder in real Heisenberg antiferromagnets is site dilution, realized by doping of non-magnetic impurities in the magnetic lattice [106]. Bond dilution can, in principle, be realized by vacancies in the non-magnetic lattice, in particular by vacancies on the sites of the non-magnetic ions involved in the super- exchange paths. The presence of such vacancies naturally alters the local electronic structure of the system, such that its effect on the magnetic Hamiltonian may be more elaborate than simple bond dilution, and it might even introduce excess charge carriers in the system. The doping of the non-magnetic ions involved in the superexchange paths with other non-magnetic species introduces instead bond disorder [69, 48], whose effect can in principle be arbitrarily close to that of bond dilution. Beside the technical challenges in realizing bond dilution (and bond disorder in general) in real antiferromagnets, a fundamental aspect of this work is the central role of inhomo- geneity in the search for novel quantum phases induced by disorder in two-dimensional Heisenberg antiferromagnets. Recent studies [59] show that the N´ eel ordered state of the square-lattice Heisenberg antiferromagnet is extremely stable towards homogeneous bond disorder, and homogeneous bond dilution destroys N´ eel order only at the classical 84 percolation threshold [83]. Hence, inhomogeneity is an essential ingredient to realize strong enhancement of quantum effects through disorder beyond the previously investi- gated scenarios. 85 Chapter 4 Field-Induced Quantum Disordered Phases in a Two-Dimensional Spin-Gapped System with Site Dilution 4.1 Introduction Quantum phase transitions in spin gapped antiferromagnets (AF) have been exten- sively studied both theoretically and experimentally. The mechanism driving the sys- tem from a quantum disordered ground state with a finite spin gap to a gapless long- range antiferromagnetic ordered (LRAFO) state is now well understood. The usual way to disturb and close the spin gap is realized through either applying a mag- netic field [74, 121, 33, 122, 103, 9, 107, 38, 86, 64, 78] or doping the system with non-magnetic impurities [104, 47, 4, 70, 68]. However the physics behind the two methods is quite different. In one of the spin-gap systems, the weakly coupled 86 dimers, the finite singlet-triplet gap can be overcame by an applied magnetic field. This leads to a field-induced quantum phase transition, which can be well under- stood as a Bose-Einstein condensation (BEC) ofS = 1 triplets excited by the applied field [64, 1, 27, 49, 50, 40, 56, 111, 109, 65]. The field suppresses correlations parallel to the field but induce transverse LRAFO. Alternatively doping with non-magnetic impuri- ties to the spin-gap system has dramatic effects. In weakly coupled dimer systems each dopant releases a spin 1/2 soliton confined to the impurity [94, 60, 67], leading to the formation of localS = 1/2 magnetic moment around each impurity. The effective local moment can be approximated to the unpaired spin neighboring to the impurity. The wave function of the local moment decays exponentially in space. However the over- lap of wave functions produces a non-frustrated effective coupling between two local moments that decays exponentially with distance between the two impurities [91]. Thus in two dimension (2D), doping introduces a random network ofS = 1/2 unpaired spins coupled through non-frustrated long-range interactions beyond the coupled dimer back- ground. This 2D random network of effective spins destroys the spin gap immediately by fill the gap with low-lying excitations, and thus induces LRAFO atT = 0, [41, 42] which is known as the order-by-disorder (OBD) phenomenon in spin liquids [89]. It is very interesting to ask what is the ground state of the weakly coupled dimer sys- tem and how it evolves when both relevant factors, the magnetic field and non-magnetic impurities, are present. We will show, through a spin-boson mapping, generally a Bose 87 glass (BG) phase appears when both factors are active. However, the weakly interact- ing local moments and the strongly interacting dimers response the field so differently that we must distinguish the two quantum disordered phases associated with these two degrees of freedom. O B D Disordered Free Moment Plateau Bose Glass Transverse Ordered 0 h DFM h PL 0 h BG h Figure 4.1: A sketch of the phase diagram of the site-diluted weakly coupled dimers in the presence of a magnetic field. The labels on theh axis denotes the critical fields at each phase transition (see text in Sec. 4.2 and Sec. 4.3 for details). OBD: order-by- disorder phase. In this chapter we investigate the phase diagram of a 2D site-diluted coupled dimer system in the presence of a magnetic field. A detailed description of the model is given in next section. The phase diagram of this model in the weak coupling limit is found to be very rich, as sketched in Fig. 4.1. In this system the energy scales associated with the interacting local moments and the intact dimers are well separated. The system is then expected to experience a series of order-disorder phase transitions: the destruction of the OBD phase to a gapless disordered-free-moment (DFM) phase at low field, then the emergence of a transverse ordered phase from a gapless BG phase at a higher field. 88 In between the DFM and BG phase, there is a local transition to a gapped plateau (PL) phase, where all the local moments are fully polarized by the field. Note that a certain strong field will finally align all the spins in the system. Hence there should be an extra phase transition from the transverse ordered phase to a fully polarized phase at high field. But the physics there is pretty similar to what happens in lower fields [76]. So here we will limit our discussions within the OBD phase, the field-induced transverse ordered phase, and the disordered phases in between. We reserve the discussions on the phase transition from the transverse ordered phase to the fully polarized phase at high field and the properties of related quantum disordered phases to future work. The DFM phase is a field-induced disordered phase immediately after the destruction of the OBD phase. In this phase the local moments induced by the site dilution are not fully polarized by the uniform magnetic field. At a higher field the full polarization of these free moments give rise to a magnetization plateau, and a finite spin gap. At even higher field close to the critical field in the clean limit, the gap closes because of the appearance of triplets excited by the field. However the excited triplets (S = 1 bosons) are still localized due to the site dilution. The system then stays in a BG phase [24], until the bosons can develop real long-range order at higher fields. The above phase diagram is quite similar to the one obtained in a previous study on a bilayer system [76]. But in this chapter we will more concentrate on and systematically study the the magnetic properties of the DFM phase and the OBD-to-DFM phase transition. The motivation comes from two sides. On the numerical side, to avoid the expensive cost of accessing 89 ultra-low energy scales previous studies mainly focused on the problem of BEC on disordered lattices and the related Bose glass. So the magnetic response of the diluted system to low fields is still unclear. On the experimental side, it is much easier to study the system at low fields than at high fields. This makes our current study more relevant to experimental observations. To understand the nature of the DFM phase we map the original diluted system to a system of randomly distributed effective spins with exponentially decaying long-range interactions [41, 42], and propose a quantum magnetization scenario in this phase. In this scenario, the inhomogeneous distribution of the exchange couplings stabilizes for- mation of strongly quantum fluctuating close-lying clusters, such as nearest-neighboring dimers and trimers. Numerical investigations on the local magnetization and bond ener- gies support this picture. The close-lying clusters have dramatic effect on the magnetic properties in the DFM phase, leading to a series of pseudo-plateaus (PPs) in the uniform magnetization curve. By further mapping the effective spin model to a interacting hard- core boson system, we reveal the explicit relation between the cluster gap statistics and the uniform magnetization, from which we may quantitatively calculate the magnetiza- tion value and the field position of each PP. After the detailed scrutiny of the disorder state we go back to the discussion on phase diagram of the diluted system. We will not limit to the weak coupling limit, but propose a phase diagram for general inter-dimer couplings. Different from the mean-field phase 90 diagram given in Ref. [54], novel quantum disordered phases, such as BG and DFM phases, are naturally included in our phase diagram. To construct a connection between the OBD-to-DFM transition at low field and the intermediate-field BG-to-SF one, we study the scaling behavior of correlation length, antiferromagnetic order parameter, and the spin stiffness. Corresponding critical exponentsν, β, andz are extracted from the scalings. The following values are obtained at both transitions: ν = 1.0± 0.1, β = 0.9±0.1,z = 2.0±0.1. Interestingly we find that a universal set of critical exponents apply to both transitions. This confirms the argument that the DFM phase is another low-field BG phase, and the OBD-to-DFM transition also follows the 2D BG-to-SF universality class. This chapter is organized as follows: in the next section we describe in detail the models appeared in this paper together with the numerical details on these models. In Sec. 4.3, we present our QMC results on the phase diagram of 2D site-diluted weakly coupled dimers. In Sec. 4.4 propose the quantum scenario in DFM phase and test its validity by studying distribution of local magnetization and bond energies. We emphases that the formation of local clusters is crucial to understand the inhomogeneous uniform magne- tization curve. We introduce a spin-boson mapping to explain the complicated magneti- zation curve, especially the PPs. In Sec. 4.5 we give the phase diagram of 2D site-diluted coupled dimer system at general coupling strength in the presence of a magnetic field. We also study the scaling behavior at OBD-to-DFM and BG-to-SF transitions, and show that they belong to the same universality class. 91 4.2 Models and Numeric Method To quantitatively investigate the field response of a site-diluted spin-gapped antiferro- magnet we consider a 2D S = 1/2 Heisenberg antiferromagnetic model of weakly coupled dimers on a square lattice [118, 51]. The Hamiltonian reads H =J X i,j ǫ 2i−1,j ǫ 2i,j S 2i−1,j · S 2i,j +J ′ X i,j ǫ 2i,j ǫ 2i+1,j S 2i,j · S 2i+1,j +J ′ X i,j ǫ i,j ǫ i,j+1 S i,j · S i,j+1 −h X i,j S z i,j . (4.1) ǫ takes 0(1) with probabilityp(1−p). The subscriptsi andj arex andy coordinates of a site on the lattice. J and J ′ refer to the intra-dimer and inter-dimer couplings respectively. In Sec. 4.3 and 4.4 we consider the situation of weak coupling, i.e.,J ′ ≪ J. But in Sec. 4.5 we investigate the phase diagram for generalJ ′ 6J. In the clean limit (p = 0) and zero magnetic field there is a criticalJ ′ c ≃ 0.523 separating the disordered phase from the ordered one [51]. For J ′ ≪ J the ground state is in a quantum disordered phase with a finite spin gap Δ 0 ∼ J ≫ J ′ . Any finite but arbitrarily small dilutionp drive the disordered ground state to an ordered one through an OBD mechanism described in Sec. 4.1. If interest only the low-lying excitations 92 below Δ 0 , for any finite doping and h ≪ Δ 0 , the above diluted model can be well described by an effective spin1/2 Hensenberg model. The Hamiltonian reads H = X i<j J ij S i · S j −h X i S z i . (4.2) Here theS = 1/2 spins are the effective local moments, which approximately equal to those unpaired spins, in the diluted model. Hence, on a L×L lattice the number of effective spins isN s =pL 2 . Since the diluted sites are randomly located in the system, in the effective model, spins are also randomly distributed on a square lattice. In the first term of Eq. 4.2 the summation runs over all the pairs of spins in the system, which means long-range pairwise interactions are considered. The effective pairwise coupling J ij takes the following asymptotic form in the long-distance limit,r ij →∞: J ij = (−) |i−j|−1 J 1 r α ij e −r ij /ξ 0 . (4.3) whereJ 1 is a parameter determined byJ andJ ′ in the original dilute model. ξ 0 is the correlation length of the gapped ground state in the clean limit. Here we have already set the lattice spacing to 1. J ij has an alternate sign depending on the Manhattan distance |i−j| between sitei andj. The alternate sign ofJ ij reflects the non-frustrated nature of the interactions between the effective spins, and give rise to the AF ordering at zero field.α is a model dependent exponent, ranging from 0 to 1 in 2D. 93 Now we discuss the values of some parameters appears in our effective model. This is very important since they tell how good the effective model can reflect the original dilute model. For weakly coupled dimers, ξ 0 is very small, usually ξ 0 < 1. But the interaction will lose its long-range character ifξ 0 is too small. So in this study we take ξ 0 = 1. We take α = 1 based on a second-order perturbation theory [91]. Although different theory [91, 76, 54, 108] gives different estimate ofα, in the long-distance limit the dominating factor is actually the exponential terme −r ij /ξ 0 . The algebraic decayed term just contributes a weak correction to this exponential decay. Eq. 4.3 works well only for large r ij . Corrections are necessary for r ij ≤ ξ 0 . Since we take ξ 0 = 1, and there is no magnetic structures smaller than one lattice spacing, we only need to treat J(r ij = 1) specially. From Eq. 4.3 we get J(r ij = 1) = J 1 e −1/ξ 0 . J 1 can be easily estimated from the second-order perturbation to be in the order of J ′2 /J. As a consequence J(r ij = 1) ∼ J ′2 /J. However in the diluted model the two unpaired spins separated by one lattice spacing couples to each other with strengthJ ′ . To loyally reflect the diluted model we have to enforceJ(r ij = 1) = J ′ . This constraint leads to a singular distribution of J ij at r ij = 1. To relieve the singularity and for simplicity, we takeJ 1 = J ′ e 1/ξ 0 such thatJ ij = J ′ if sitei and sitej are separated by one lattice spacing. This will overestimate the effective couplings at larger ij . But our numerical results show that by taking current parameters the effective spin model captures main characters of the original diluted coupled dimer model. Now Eq. 4.3 changes to J ij = (−) |i−j|−1 J ′ r ij e −(r ij −1)/ξ 0 . (4.4) 94 We apply the Stochastic Series Expansion (SSE) QMC [99] to both the original and effective Hamiltonian, given by Eq. 4.1 and 4.2, to investigate the evolution of the ground state of the system in the presence of a magnetic field. For the effective model with N s randomly distributed spins and non-frustrated long-range interactions among them the computational effort can be reduced from∼ N 2 s to N s lnN s . [85] In dealing with both models the simulations are performed using a β-doubling approach [83] up to inverse temperatureβ/J = 2 15 to obtain theT = 0 feature of the system. Disorder averaging is typically performed over≈ 300 realizations. The long-range AF ordering is reflected in the appearance of transverse staggered mag- netization m s = p S ⊥ (π,π)/L 2 , (4.5) where L is the linear system size and S ⊥ (π,π) is the transverse staggered structure factor defined as S ⊥ (π,π) = 1 2L 2 X i,j ′ (−) |i−j| hS x i S x j +S y i S y j i. (4.6) The summation runs over all nearest neighbor pairs in the site-diluted model, but runs over all possible pairs of spins in the effective spin model. We also estimate the spin 95 stiffness, which turns to the superfluid density in the bosonic language, through counting the winding numbersW x,y in the two lattice directions of SSE worldlines [73]: ρ s = 1 2βJ hW 2 x +W 2 y i. (4.7) 4.3 The Phase Diagram in Weak-Coupling Limit In this section we will discuss the field evolution of the ground state of the site-diluted weakly coupled dimers. We start by giving a field scan of the uniform and staggered magnetization, identify each phase induced by the field, and then focus on the quantum disordered phases. In our study of the weakly coupled dimer system we take J ′ /J = 1/4. In the clean limit, it is far away from the critical J ′ c ≈ 0.523 at zero field, and the correlation length ξ 0 ≈ 0.99. When a magnetic field is applied the clean system experiences an order-disorder QPT at h (0) c1 = Δ 0 = 0.603(5). Doping with non-magnetic impurities drastically changes the ground state. The phase diagram of the doped system can be obtained by analyzing the staggered magnetization together with the uniform magneti- zation (m u = P i hS z i i/L 2 ) data. 96 0.00 0.01 0.6 0.7 0.00 0.05 0.10 0.0 0.2 0.4 0.6 0.00 0.05 0.10 BG PL m u h u /J DFM SF PL SF BG DFM OBD L = 16 L = 24 L = 32 L = 40 L m s h/J OBD h/J Figure 4.2: The field dependence of the uniform magnetization, uniform susceptibility, and staggered magnetization in 2D site-diluted weakly coupled dimers. The coupling ratioJ ′ /J = 1/4, and the dilution concentrationp = 1/8. In the lower plot the orange diamonds show the staggered magnetization data in the thermodynamic limit, which are extracted from the finite-size scaling results. SF: the transverse ordered phase, which corresponds to a superfluid phase in the bosonic language. Our QMC result of the field dependence of uniform and staggered magnetization for the site-diluted coupled dimers withJ ′ /J = 1/4 and at a moderate dilution concentration p = 1/8 is shown in Fig. 4.2. For 0 ≤ h/J ≤ 0.725 the numerically found phase diagram is similar to the one sketched in Fig. 4.1. At zero field a finite staggered magnetization m s = 0.032(3) is extrapolated in the thermodynamic limit, indicating the system is in the OBD phase. We find that this 97 OBD phase also survives at finite fields up toh = h DFM = 0.007(1)J. Remarkably the critical field that destroys the OBD phase h DFM ≪ J ′ , implying that the average coupling strength responsible for the long-range AF order turns out to be much smaller thanJ ′ . At h DFM the uniform magnetizationm u = 0.0208(6), is much less than the value m sat u = pS = 1/16 corresponding to fully polarized free moments. It is not saturated tom sat u until a much higher fieldh =h PL ≈ 0.5J is reached. Hence the DFM phase, which existing ath DFM < h < h PL , has a striking feature: the local moments are not fully polarized in this phase; whereas in a homogeneous antiferromagnet, the destruction of long-range order is always accompanied by the full polarization of the local free moments. Another dramatic feature of the magnetization curve in the DFM phase is: beside the large plateau appearing ath ≥ h PL , one observes the presence of apparent intermediate plateaus at around73% and 95% of the saturation magnetization. We will see through a detailed study in next paragraph that these features are actually pseudo-plateaus (PPs), which retain an extremely small slope (Fig. 4.2). The first PP extends up to h ≈ 0.7J ′ ; a second PP markedly appears around h ≈ 1.2J ′ ; the true saturation plateau is only attained ath≈ 2J ′ . Beside the uniform magnetization, we also studied the uniform susceptibility, which is calculated by taking χ u = dm u /dh, The QMC result is also presented in Fig. 4.2. It clearly shows that the DFM phase retains a finite uniform susceptibility. We also 98 1 10 100 0.01 0.1 u J h/J=0.175 L=24 L=32 Figure 4.3: The temperature dependence of the uniform susceptibility in 2D site-diluted weakly coupled dimers at the first magnetization PPh/J = 0.175. The coupling ratio J ′ /J = 1/4, and the dilution concentrationp = 1/8. calculated the finite temperature behavior of χ u at the PPs. If those are real plateaus with a finite gapΔ, the susceptibility vanishes as χ u ∼e −Δ/T . (4.8) This is just the case for the saturation plateau at h = h PL . But the QMC result in Fig. 4.3 suggests at PPs, the dependence of χ u on T is power-law. Moreover, when T → 0, χ u → χ u,0 ≈ 0.005 > 0, which implies the DFM phase is gapless. These unconventional properties have also been observed in the magnetic BG phase of triplet 99 quasi-particles living on intact dimers [76, 24, 77, 66]. Later we will argue that these two phases bear indeed strong analogies, although involving different degrees of freedom. Forh>h PL the field completely aligns all the unpaired spins, but it is not strong enough to break the strongly coupled dimers. The saturation of the magnetization of FMs then leads to the full restoration of a gapped disordered phase [54]. The magnetization keeps to the saturated valuep/2 until the field reaches the value corresponding to the gap of the clean system,h = Δ 0 = 0.60(1)J. Whenh> Δ 0 the gap to excite a triplet closes just as in the clean system. The reason is in the diluted system one can always find an arbitrarily large clean region devoid of any impurities with an exponentially small but finite probability in the thermodynamic limit. The local gap of this clean region is arbitrarily close toΔ 0 . Once the field overcomes this local gap it generates a triplet excitation in the clean region. What completely different from the clean limit is that the excited triplets are localized in the clean regions. No long-range order is induced. The localization effects come from two aspects. The fully polarized local moments have a finite gap to spin-flips, and the intact dimers close to the impurities have relatively higher local gaps than in the clean regions. Both these two factors serve as energy barriers to the triplets in the clean regions. This picture can be represented in bosonic language as S = 1 interacting bosons localized by a random potential. The ground state of the bosonic triplets is denoted as the Bose glass. It corresponds to the field range 0.6 6 h/J 6 0.67(1) in Fig. 4.2. According to the 100 above argument, in BG phase although there is still no magnetic order, the magnetization slowly increases fromp/2, leading to a finite susceptibility. The field drive the system back to an transverse-ordered state through a QPT at h BG /J ≈ 0.67. In bosonic language the transverse-ordered state is a superfluid state where the triplets can propagate coherently. The transition is thus a BG-to-SF transi- tion. It has different critical properties from the Mott insulator (MI) to SF transition directly happens in the clean system. 4.4 Field-Induced DFM Phase In previous section we have shown that in the DFM phase many magnetic properties behave unconventionally. In this section, we systematically investigate several quantities in this phase and explain those unconventional behaviors within a novel picture of the DFM phase. 101 4.4.1 Classical vs Quantum Views on DFM Phase We notice that the original site-diluted coupled dimer model described by Eq. 4.1 is too complicated. There are two basic energy scales: the intact dimers with energy scale E ∼ J, and unpaired spins with energy scale E 6 J ′ . In the DFM phase since h < h PL ∼ J ′ ≪ J the dimers are in fact irrelevant. To simplify the problem and to better understand the physics in this phase, we need a low-field effective model with only relevant degrees of freedom which have energy scaleE6J ′ . It has been shown [41, ?, 119] that the effective model introduced in Eq. 4.2 captures the main features of the field dependencem u andm s in the diluted coupled-dimer model. It is the minimum model that gives the correct physics of DFM phase. Let us think about the classical limit of the effective model. WhenS→∞, the system of randomly distributed spins in a weak enough magnetic field has a canted antiferromag- netic ground state (Fig. 4.4(a)). At higher fields, although some of the spins may locally be fully polarized by the field, the system of remaining spins will preserve long-range order since their transverse components are still staggered because of the long-range exchange interaction among them. This picture holds until eventually the field exceeds the largest effective coupling, i.e., h > J ′ , where all of the spins are fully subject to the magnetic field. This means the destruction of long-range antiferromagnetic order is always accompanied by the fully saturation of the magnetization in the classical limit. 102 Classical J eff >h Quantum J eff <h (a) (b) Figure 4.4: (a): The classical S → ∞ free moments has a canted antiferromagnetic order at finite field. (b): In the quantum case, the close-lying free moments prefer to form local clusters, leading to a novel DFM phase. The black and red arrows distinguish spins located on different sublattices. Quantum mechanically there is a substantially different scenario [76, 119]. At zero field, a majority of free moments are spaced from each other by an average distance hr ave i = 1/ √ p, large in the small dilution limit, and interact via average couplings hJ eff i≪J ′ . However, fluctuations in the spatial distribution of the impurities also lead to small clusters of free moments with linear size r ≪ r eff , and thus interacting with much stronger couplingsJ eff ∼ J ′ . Examples are the nearest neighboring dimers and trimers. A magnetic fieldh&hJ eff i destroys the long-range order of the free moments. But a finite-size gap Δ i ∼ J ′ prevents the clusters of close-lying moments from being fully polarized ath.J ′ . This quantum mechanical picture (in Fig. 4.4(b))results in the 103 fundamental consequence that the antiferromagnetic order disappears at a very low field but the free moments are far from saturation. 4.4.2 Local Clusters of Free Moments In last section we proposed two scenarios based on the classical and quantum effective models. They lead to different field evolution ofm u . The difference between the clas- sical and quantum pictures is reflected by the statistics of several local quantities. By directly compare the results from the effective models with the one from the original diluted model, we will see that the quantum scenario captures the dominant physics in the DFM phase. Let us first take a look at the field evolution of m u and the transverse spin compo- nent S x(y) in the classical effective model. The QMC result is presented in Fig. 4.5. At h = J ′ , the uniform magnetization hits its saturated value, and S x(y) → 0. This indicates the destruction of long-range order takes place where all the free moments lose transverse components, just as predicted by the classical picture. Interestingly the magnetization curve clearly shows the spin response to the field is still highly inhomoge- neous even in the classical case.m u increases much slower aroundh∼ 0.1J ′ , develops an apparent plateau in logarithmic scale. The magnetization value is about the same as 104 0.01 0.1 1 0.4 0.6 0.8 1.0 0.01 0.1 1 0.0 0.2 0.4 0.6 m u /m u sat classical quantum <S x(y) >/S h/J’ Figure 4.5: Field dependence ofm u and transverse spin componentS x(y) in the classical effective spin model withξ 0 = 1.0,p = 0.0625 andL = 40. The correspondingm u in the quantum model with the same model parameters is also given as a comparison. the one of the first PP in the quantum model. It suggests the magnetization value of the PP is determined by the inhomogeneous distribution of local moments. However during canting the transverse correlations among the spins antiparallel to the field keep increasing, so that S x(y) is stabilized to a finite value for 0.1J ′ < h < J ′ , showing the long-range order is still preserved. Therefore in the classical model the transition is always from a gapless ordered phase to a gapped phase with full saturation ofm u , and the DFM phase is a consequence of quantum many-body effect. 105 -0.4 -0.2 0.0 0.2 0.4 1E-5 1E-4 1E-3 0.01 0.1 1 P(S i z ) <S i z > h/J’=0.1 h/J’=0.2 h/J’=0.4 h/J’=0.6 h/J’=0.8 Figure 4.6: Distribution of thez component of local momentsS z i in the effective spin model for classical spins withξ 0 = 1.0,p = 0.0625, andL = 40. The histogram of local spin componenthS z i i provides a way to visualize the local spin state and to check the validity of the proposed scenarios. The distribution ofhS z i i in the classical model is shown in Fig. 4.6. At zero fieldhS z i i is uniformly distributed between −0.5 and 0.5. The presence of a very weak field immediately breaks the symmetry and induces two peaks located athS z i i =±0.5. The peak athS z i i = 0.5 corresponds to spins fully polarized by the magnetic field, whereas the peak athS z i i =−0.5 corresponds to spins antiferromagnetically coupled to the fully polarized ones. Forh/J ′ > 0.1, an extra peak in between appears. This peak moves towardhS z i i = 0.5 with increasing magnetic 106 field, corresponds to the partially polarized spins. At the same time, the peak athS z i i = −0.5 drops. This is consistent with the picture that classical spins are continuously polarized. Whenh =J ′ almost all the spins are fully polarized by the field, and thus the uniform magnetization saturates. Remarkably, up toh = 0.6J ′ there is still a significant portion of the spins antiparallel to the field, which accounts for the stabilization of finite transverse spin componentsS x(y) . We then show the result of the quantum model as a comparison. According to the quan- tum scenario described in previous section, we expect to see effects from the formation of local clusters, especially from the nearest-neighboring dimer and trimers. This is clearly seen in Fig. 4.7. At zero field there is a huge peak centered athS z i i = 0. There is a clearly double-peak structure at finite fields. Besides the peak athS z i i = S = 0.5 coming from the contribution of fully polarized spins, what surprisingly different from the classical case is, the central peak is stable at a finite field h 6 J ′ . It implies that there is a strong singlet component in the ground state wave function, coming from the local clusters with even number of free moments. As another difference from the clas- sical casem we do not observe any peak athS z i i = −0.5 in the quantum spin model. However we can resolve two extra peaks at hS z i i = 1/3 and hS z i i = −1/6 at finite fields, which correspond to partially polarized spins in free-moment trimers. Local free- moment clusters have widely different local gaps to full polarization, both due to their geometric structure (dimers, trimers, quadrumers, etc.), and to the local field they expe- rience from the other free moments. Yet the distribution of rare free-moment clusters 107 -0.2 0.0 0.2 0.4 0.01 0.1 1 10 P(<S i z >) <S i z > h/J’=0.63 h/J’=0.91 h/J’=1.09 h/J’=1.36 h/J’=1.81 -0.2 0.0 0.2 0.4 0.01 0.1 1 10 P(<S i z >) h/J’=0 h/J’=0.001 h/J’=0.01 h/J’=0.09 h/J’=0.27 h/J’=0.45 Figure 4.7: Distribution of thez component of local momentsS z i in the effective spin model for quantum spins withξ 0 = 1.0,p = 0.0625, andL = 40. clearly assigns dominant statistical weight to the dimers, and this simple geometric fact is the reason for the appearence of the first pseudo-plateau in them u curve: the magne- tization process nearly stops until the local gap of the dominant free-moment dimers is overcome ath . J ′ . Nonetheless, the slope remains finite because the dimers have a broad distribution of local gaps. 108 1E-3 0.01 0.1 1 -12 -8 -4 0 -0.6 -0.4 -0.2 0.0 0.00 0.01 0.02 0.03 0.04 P(ln< E b >) h/J’=0 h/J’=0.02 h/J’=0.20 h/J’=0.40 h/J’=0.80 h/J’=1.00 h/J’=1.20 h/J’=1.40 P(ln< E b > ) ln< E b > h/J’=0 h/J’=0.4 h/J’=0.8 h/J’=1.0 h/J’=1.2 trimer dimer Figure 4.8: Upper panel: distribution of the logarithm of the reduced bond energyh ¯ E b i in the effective spin model withξ 0 = 1.0,p = 0.0625, and system sizeL = 40. Lower panel: a closer look in the vicinity oflnh ¯ E b i = 0 and−0.4, indicating the formation of close-lying dimers and trimers. We may confirm the formation of nearest-neighboring free-moment dimers and trimers by further investigating the distribution of bond energies. A local singlet with coupling J ′ has the bond energy E b = − 3 4 J ′ . We define a reduced bond energy on each bond, ¯ E b =−4E b /(3J ′ ), and plot the field dependence of the distribution oflnh ¯ E b i in Fig. 4.8. There is a broad continuous low-energy (in absolute value) tail and several discrete 109 higher energy peaks. With increasing field the low-energy tail move towards high energy end and gets narrower, indicating that the low-energy degrees of freedom are suppressed by the magnetic field. However the peaks at higher energy scales are more stable. We resolve two peaks atlnh ¯ E b i = 0 and−0.40, corresponds to dimers with bond energy E b =− 3 4 J ′ , and trimers withE b =− 1 2 J ′ . Furthermore we study the finite-size scaling of the peak values of the dimer peak. As shown in Fig. 4.9 this peak decays when increasing system size. But as shown in the inset of Fig. 4.9 it is inversely proportional topL 2 −1, the coordination number of a free moment. We notice that the total number of bonds in a system withN s spins is N s (N s − 1)/2, and N s = pL 2 . So the scaling behavior just confirms that the number of dimers with couplingJ ′ is proportional toN s . The study on distribution ofE b andS z i confirms the existence of local free-moment clus- ters in DFM phase proposed in the quantum scenario. To see this scenario is applicable to the original diluted coupled-dimer model, we study the distribution of local magneti- zation of unpaired spins and bond energies in the original model. The result ofhE b i at J ′ =J/4 is shown in Fig. 4.10(a). The spectrum is really rich. Besides the complicated structure at low energy (absolute value) end, we identify three peaks athE b i≈−3J/4, corresponding to singlets on intact dimers; hE b i ≈ −3J ′ /4, corresponding to free- moment dimer singlets; andhE b i≈−J ′ /2, corresponding to free-moment trimers. 110 -0.6 -0.4 -0.2 0.0 0.2 0.00 0.05 0.10 0.00 0.02 0.04 0.06 0.08 0.00 0.05 0.10 P(ln< E b >) ln< E b > h/J’=0.008 p=0.0625 L=16 L=24 L=32 L=40 L=48 P(ln< E b >=0) 1/(pL 2 -1) Figure 4.9: Finite-size study oflnh ¯ E b i ath/J ′ = 0.008 in the effective spin model with ξ 0 = 1.0 andp = 0.0625. The inset shows finite-size scaling ofP(lnh ¯ E b i); the red line is a linear fit. The result of distribution of unpaired spinshS z i i is shown in Fig. 4.10(b). We see exactly the same signature as the one plotted in Fig. 4.7. At zero field there is a sharp peak at hS z i i = 0 This peak still keeps a significant finite value in the presence of a magnetic field up toJ ′ . This shows that the formation of strongly interacting singlets still plays an important role in the diluted model. Besides this singlet peak, there are peaks at abouthS z i i = 1/3 andhS z i i = −1/6 in finite fields, clearly showing the existence of free-moment trimers. 111 -0.8 -0.7 -0.6 -0.2 -0.1 0.0 1E-5 1E-4 1E-3 0.01 0.1 FM trimers FM dimers P(E b ) E b /J h/J=0 h/J=0.05 h/J=0.2 spin dimers (a) -0.4 -0.2 0.0 0.2 0.4 1E-3 0.01 0.1 P(<S i z >) <S i z > h/J=0.005 h/J=0.05 h/J=0.15 h/J=0.25 (b) FM dimers FM trimers Figure 4.10: Histogram of bond energies (a) and local magnetizations for unpaired spins (b) in the diluted weakly-coupled-dimer model. The model parameters areJ ′ /J = 1/4, p = 1/8, and system sizeL = 40. Let us end this section with a little summary. In this section we verified the proposed quantum scenario in the DFM phaseby investigating the distribution of local magneti- zationhS z i i and bond energyhE b i for both the effective model and the original model. In both models we confirm that there are free-moment clusters with nearest-neighboring 112 coupling J ′ formed. These local clusters do have significant effects on the magnetic properties in the DFM phase, just as discussed in last section. Such local clusters are not observed in the classical spin model, indicating that the DFM phase is completely a quantum disordered phase with no classical counterpart. 4.4.3 Statistics on Local Clusters and Pseudo-Plateaus Here we are going to discuss the highly inhomogeneous magnetization process in the DFM phase. We first map the effective model to a model of interacting bosons. Due to the quantum scenario this effective Hamiltonian can be treated based on a cluster decomposition. We then take the mean-field approximation and show the magnetization value and the field location of the PPs are quantitatively related to the statistics of free- moment clusters. Finally we compare our mean-field result with the one from QMC simulations. We first map the effective spin model given in Eq. 4.2 into a model of interacting hard- core bosons. It is found that the effective bosonic model is more convenient for under- standing the highly inhomogeneous magnetization process in DFM phase. We describe the mapping as follow. The effective spin model describes a system of randomly dis- tributed interacting spins on a bipartite lattice. After a rotation ofπ about theS z axis in 113 the spin space for all the spins on one sublattice, it can be mapped to aS = 1/2 sys- tem with only ferromagnetic couplings. The quantum ferromagnetic system is equiv- alent to a interacting hard-core boson system if we perform S + → b † , S − → b, and S z →n−1/2 =b † b−1/2: H =− X i<j t ij (b † i b j +b † j b i )+ X i<j V ij n i n j − X i μ i n i , (4.9) where the first term describes the hopping; the second term is neighboring interaction which can be either attractive or repulsive, depending on the sign of spin-spin coupling in the effective spin model; andμ i in the last term is an effective chemical potential.t ij , V ij , andμ i are related with model parameters in the effective spin model through: t ij = 1 2 |J ij |; (4.10) V ij =J ij ; (4.11) μ i =h+ 1 2 X j V ij . (4.12) Note that both the hopping and the interaction are between any two sitesi andj. The chemical potential includes both the applied magnetic field and an effective local field (see the second term in Eq. 4.12). According to the quantum scenario, a majority of spins participates the long-range antiferromagnetic order are weakly coupled among them. Once the long-range order is destroyed ath≈J N´ eel ≪J ′ , they are fully polarized by the field. Then the only active degrees of freedom are the free moments belonging 114 to the local clusters with strong coupling J eff ∼ J ′ . Spins in these clusters are not fully polarized in the DFM phase. In bosonic language only the sites in the clusters are not entirely occupied by the bosons, or equivalently, holes can only be excited in these clusters. This leads tohn i i ≈ 1. So it is more convenient to use the language of hole excitations. Apply a particle-hole transformation,b † →a,b→a † ,n→ 1− ˜ n, Eq. 4.9 is rewritten as H =− X i,j t ij 2 (a † i a j +a † j a i )+ X i,j V ij 2 ˜ n i ˜ n j − X i ˜ μ i ˜ n i , (4.13) where ˜ μ i = P j V ij 2 −h, is the chemical potential to generate a hole. Since the coupling decays exponentially in real space (see Eq. 4.4), the hopping between two sites in a same cluster is much stronger than the inter-cluster hopping. As an approx- imation, we may decompose the Hamiltonian in Eq. 4.13 into a sum of contributions from independent clusters, and drop the inter-cluster hopping terms. This leads to H = X C H N C (C), (4.14) 115 whereN C = 1,2,3,··· denotes the size of the clusterC. We treat the fully polarized spins, or the sites in absence of holes, as a single-site cluster. The hopping terms are all dropped for these clusters, and the Hamiltonian reads H 1 (1) = X j V 1 j 2 ˜ n j ˜ n 1 − ˜ μ 1 ˜ n 1 . (4.15) For clusters withN > 1, holes can be excited due to the strong intra-cluster hopping, the Hamiltonian is then H N>1 (C) = − X i,j∈C t ij (a † i a j +a † j a i )+ X i,j∈C V ij ˜ n i ˜ n j − X i∈C ( X j∈C V ij 2 −h)˜ n i + X i∈C,j∈ ¯ C V ij ˜ n i ˜ n j − X i∈C ( X j∈ ¯ C V ij 2 )˜ n i . (4.16) Note that we keep the neighboring interactions V ij ˜ n i ˜ n j between sites from different clusters. We then treat these terms in a mean-field way: X j∈ ¯ C V ij ˜ n i ˜ n j ≈ ( X j∈ ¯ C V ij h˜ n j i)˜ n i . (4.17) However we have to rigorously solve the part of intra-cluster interactions to obtain the correct value of the local gap to annihilate a hole in each cluster. It becomes exponen- tially difficult when the cluster size keeps growing. But thank the inhomogeneity of the system, the probability of growing a large cluster is also exponentially small. We can 116 then set a cutoff of the cluster size without introduce any significant systematic errors. In fact, we find it has been a really good approximation to diagonalize up toN C = 3, i.e., the trimers, and simply set the local gap Δ C = 2J ′ for clusters withN C > 4. The choice of 2J ′ is based on the following consideration: in the original diluted coupled- dimer model, we find that the geometry of the close-lying free-moment clusters are not arbitrary, but limited by the lattice frame in real space; for each such cluster, we find the upper bound of the local gap is2J ′ , corresponding to the gap of an infinitely long chain. The entire Hamiltonian can then be solved in a self-consistent way. Given a first estimate of theh˜ n i i, exact diagonalize each local cluster to find a set of newh˜ n i i. Then use the new set of values as the input and repeat the same process untilh˜ n i i converges. But we hereby give a much simpler way. Considering thath˜ n i i ≈ 0 for most of the sites, we may assume thath˜ n j i = 0 in Eq. 4.17. Moreover the penalty to generate more than one hole in a cluster is relatively high at large h, so that we assume at most one hole per cluster can be accommodated. In this way we can also drop the second term in Eq. 4.16. Finally we find the Hamiltonian in Eq. 4.16 is reduced to a tight-binding Hamiltonian H N>1 (C) =− X i,j∈C t ij (a † i a j +a † j a i )− X i∈C ( X j∈C V ij 2 −h)˜ n i − X i∈C ( X j∈ ¯ C V ij 2 )˜ n i , (4.18) and can be trivially solved. 117 Since for close-lying clusters, P j∈C V ij 2 ∼ J ′ ≫ P j∈ ¯ C V ij 2 , the third term in Eq. 4.18 just give a little perturbation on the local gap. This is easy to see in the expressions of local gapsΔ N C forN C 6 3: Δ 2 (1,2) ≈ J ′ + 1 2 (δ 1 +δ 2 )−h, (4.19) Δ 3 (1,2,3) ≈ 3J ′ 2 + 1 3 (δ 1 +δ 2 +δ 3 )−h, (4.20) whereδ i = P j∈ ¯ C V ij 2 . This treatment can be also applied to theN C = 1 clusters, and the local gap, within the same approximation, is simply: Δ 1 (1) =δ 1 −h. (4.21) These expressions on the local gaps can be written in a uniform form: Δ N C (C)≈ ¯ Δ N C + ¯ δ(C)−h, (4.22) where ¯ Δ N C is determined only by the cluster structure, and ¯ δ(C) = 1 N C P i∈C δ i , depends on the local environment of the cluster. ¯ δ(C) leads to a distribution of local gaps. The distribution ath = 0 is presented in Fig. 4.11. The curve shows a multi-peak signature. Each peak is centered about ¯ Δ N C with a widthh ¯ δ(C)i≪J ′ . The inhomogeneous magnetization process, especially the magnetization values and the field positions of the PPs can be quantitatively related with this distribution. We rephrase 118 -0.5 0.0 0.5 1.0 1.5 2.0 0.01 0.1 1 10 P( (C)) (C)/J’ L=80 L=128 single spin dimer trimer Figure 4.11: The distribution of local cluster gapsΔ(C) in the effective bosonic model, withξ 0 = 1.0 andp = 1/8. Each peak in the distribution can be related with an effective finite-size cluster sketched in the plot. The peak centered at Δ(C)≈ 1.9J ′ comes from more complicated cluster geometries, such as quadrumers. the magnetization process in the bosonic language. Starting from a very low field, where we assume there is on average one hole per cluster, we increase the field to annihilate holes. Once the field overcomes the local gap of a cluster, it becomes fully occupied with bosons (in absence of holes). The hole population in the system keeps decreasing as the field increases, which can be estimated by calculating the area under the distribution curve in Fig. 4.11 h˜ ni = Z P(Δ(C))Θ(Δ(C)−h)dΔ(C), (4.23) 119 where Θ(x) = 1 for x > 0 and Θ(x) = 0 for x 6 0, is a step function. Recall that h˜ ni = 1−hni = 1 2 −m u , the normalized uniform magnetization is then m u m sat u = 2 Z P(Δ(C))Θ(h−Δ(C))dΔ(C)−1. (4.24) Hence the relative magnetization values and the field positions of the PPs can be obtained by analyzing the distribution of the local gaps of clusters. Differentiate with respect toh on both sides of Eq. 4.24, we see that the uniform susceptibilityχ u is pro- portional to the distribution of local gaps. In the DFM phase, χ u > 0 since the local gaps take a continuous distribution. However, χ u can be exponentially small because each peak of the distribution of local gaps has a finite width h ¯ δ(C)i which is much smaller thanJ ′ . In the field range whereχ u is exponentially small, a magnetization PP is observed. Note that the assumption one hole per cluster (h˜ n(C)i = 1) at low fields is no longer true ath = 0, where the system is half-filled withh˜ n(C)i = 1/2. So our expression in Eq. 4.24 may not be valid at low fields. But it should work well at higher fieldsh&J ′ since the approximation is quite reasonable there. We compare them u result given by Eq. 4.24 with the QMC ones in Fig. 4.12. Remark- ably, them u /m sat u calculated from Eq. 4.24 reproduces all the PP structure. The field 120 0.0 0.5 1.0 1.5 2.0 2.5 0.2 0.4 0.6 0.8 1.0 m u /m u sat h/J’ L=40, QMC diluted dimer model L=40, QMC effective spin model L=128, Statistics on cluster gaps Figure 4.12: Field dependence of the normalized magnetizationm u /m sat u calculated by appying Eq. 4.24 in the effective model with ξ 0 = 1.0 and p = 1/8 (shown in blue triangles). The result is compared with QMC result from the effective model (shown in red circles) with the same model parameters, and with QMC result from the diluted coupled-dimer model (shown in black squares) withJ ′ /J = 1/4 andp = 1/8. positions of PPS agree well with the QMC results in both effective and the original mod- els. The agreement confirms the picture of cluster formation in the DFM phase and the relation between the PPs and the cluster statistics. We notice that the relative magnetiza- tion calculated from different ways does not perfectly fit. The deviation from the origi- nal model QMC result is due to the relative large dilution concentrationp. In Fig. 4.12, p = 1/8, corresponds to an average inter-momentr avg ≈ 2.8, close to the lattice spac- ing. In this case one cannot exclusively drop the inter-cluster interactions, and must 121 consider the effect of a finite lattice spacing. Thus we expect to see better agreement for smallerp. The magnetization value of each PP can also be determined approximately by the relative heights of the peaks in theP(Δ(C)). They are only weakly dependent of the ratioJ ′ /J, or in effective model, the correlation lengthξ 0 in the clean limit. On the contrary, they are sensitive to the dilution concentrationp. This is quite easy to see. Since the probability to find a free moment on a lattice with dilution concentration p is P(0) = p, and the probability to find a nearest-neighboring pair of gap Δ = J ′ is P(J ′ ) = 4p 2 . So thatP(J ′ )/P(0) = 4p. If we use such an approximate bimodal local gap distribution, we findm u /m sat u = P(0)/(P(0)+P(J ′ )) = 1/(1+ 4p) for the first PP. This approximation gives pretty good result for smallp, as seen in Fig. 4.13. 0.0 0.5 1.0 1.5 2.0 0.2 0.4 0.6 0.8 1.0 m u /m u sat h/J’ L=40, p=1/16, QMC L=40, p=1/8, QMC L=128, p=1/16, local gaps L=128, p=1/8, local gaps Figure 4.13: Field dependence of the normalized magnetizationm u /m sat u in the effective model withξ = 1.0, and various dilution concentrationp. Dashed blue lines show the values of 1/(1+4p) as guidelines for the magnetization of the first PP. 122 4.5 Phase Transitions in 2D Site-Diluted Coupled Dimers In this section we will discuss theT = 0 phase transition of site-diluted coupled dimer system with general coupling ratio 06 J ′ /J 6 1. We will limit the discussion to the regionh6J since the physics is well known and not quite interesting for higher fields. 4.5.1 The Phase Diagram with Site Dilution Before we start to discuss the phase diagram for the site-diluted case, we first take a glimpse on the phase diagram in the clean system. At zero field, the ground state is disordered for weak inter-dimer couplings, but ordered for strong couplings. These two ground states are separated by a critical coupling ratio(J ′ /J) c = 0.523(1). The univer- sality class at this critical point belongs to the classical three-dimensional Heisenberg model [51, 13], with ν = 0.71(1), β = 0.36(1), and z = 1. In the presence of a magnetic field, the critical point at zero field extends to a critical curve in the 2Dh-J ′ space. This curve separates the quantum disordered phase at low field and low inter- dimer coupling from the long-range ordered phase. The dynamical exponentz jumps to 2, [79] leading to mean-field values for other exponents, ν = 1/2 and β = 1/2. 123 The quantum disordered phase has a finite spin gap Δ 0 . In the weakly coupled limit J ′ ≪ J, Δ 0 ≈ J−J ′ . The phase diagram is analogous to the phase diagram of inter- acting bosons, with mapping the ordered phase to a superfluid phase and the disordered phase to a Mott insulating phase. In the weakly coupled limit such a mapping actually exists. Then the order-disorder transition in this 2D antiferromagnets realizes the 2D superfluid-insulator transition. Doping with non-magnetic impurities immediately changes the zero field ground state to a ordered one even in the weak couplings. For the weak coupling caseJ ′ =J/4 we have shown in Sec. 4.3 that the system is driven back to a disordered phase ath DFM ≪ J ′ . At a higher field h BG > Δ 0 , it experiences a second quantum phase transition back to the ordered phase. The physics of the BG phase is quite obvious: the strongly coupled dimers have both singlet and triplet components. At h > Δ 0 triplet quasi- particles appear in rare clean regions of the diluted systems. But they are localized in those regions and cannot coherently propagate ath < h BG . So no long-range order is formed in this phase. Since the gap to accommodate a triplet quasi-particle in a rare region can be arbitrarily small, the spectrum is gapless, and the system retains finite compressibility. However the DFM phase has a different physics. From the analysis in Sec. 4.4, a clear picture of the DFM phase emerges: in this phase, a majority of the free moments are polarized, but antiparallel spins exist, corresponding to rare free moment clusters. Upon a spin-to-hardcore-boson transformation, these antiparallel spins take the nature of bosonic hole quasi-particles localized on rare regions of the lattice. This 124 0.0 0.2 0.4 0.0 0.5 1.0 gapped disordered h/J J’/J p=0 p=1/8 XY ordered gapless disordered Figure 4.14: Phase diagram of site-diluted coupled dimers with dilution concentration p = 0.125. It coincides with the phase diagram of 2D interacting bosons in the presence of a random potential. The three distinct phases: gapped disordered phase (the blue regime), gapless disordered phase (the red regime), and transverse ordered phase (the white regime) corresponds to the Mott insulating (MI), Bose glass (BG), and superfluid (SF) phases in the bosonic system, respectively. In weak coupling limitJ ′ ≪J, the gap- less disordered phase is separated by the gapped phase into two branches, corresponding to the DFM phase and the intermediate-field BG phase. The phase boundary in the clean limitp = 0 (black dashed line) is also shown in the same plot as a comparison. aspect connects with the ordinary picture of a BG [24], and it further endows the OBD- to-DFM quantum phase transition with the nature of a localization transition: in the OBD phase the holes form a superfluid condensate, which is progressively depleted by the applied field (acting as a chemical potential), up to the point where the holes 125 undergo localization into a BG state, losing superfluidity but retaining compressibility, which corresponds to a finite uniform susceptibility. We have argued that the DFM phase is just another BG phase. This is clearly seen from the numerical phase diagram of the diluted coupled dimers in Fig. 4.14. The critical curve separating the quantum disordered phases from the ordered one is quite different from the clean case. The critical coupling ratio is not a monotonic function of the magnetic field. Starting from the origin, it first increases with increasing field, touching its maximum (J ′ /J) c ≈ 0.46 at abouth/J ≈ 0.29, and then decreases when the field further increases. As a result, for any coupling ratio J ′ /J < (J ′ /J) c , there are two critical fieldsh DFM andh BG . ForJ ′ /J. 0.28 there exists a plateau phase with a finite spin gap and magnetization saturating to m sat u = p/2. In bosonic language it is just the Mott insulating phase in the doping case. Between this Mott plateau phase and the ordered phase, we find a gapless quantum disordered phase with finite compressibility. In the weak coupling limit, it develops into a lower field DFM phase and a higher field BG phase. 126 4.5.2 Quantum Scaling and The Critical Exponents Next we investigate both the BG-to-SF and the OBD-to-DFM transitions. By studying the quantum scaling properties at these two QTPs we further confirm the conclusion that the DFM phase is a low-field BG phase, so that the two QPTs belong to the same universality class. We then determine this novel 2D BG-to-SF universality class by extracting some critical exponents. We will start from the weak coupling limit, then investigate the critical properties along the entire critical curve. Previous study on a site-diluted bilayer model [76] finds two similar order-disorder phase transitions. The intermediate-field QPT is identified as a BG-to-SF transition. However it was not able to determine the lower-field OBD-to-DFM transition because the associated energy scale is too small so that the singular behavior near the QCP is overwhelmed by thermal fluctuations. In 2D systems true long-range order appears only atT = 0 in the thermodynamic limit. So as an appropriate way to eliminate the effects of thermal fluctuations, we first extrapolate the T = 0 data for finite systems, then perform a finite-size scaling. 127 -2 0 2 0.0 0.2 0.4 -0.4 -0.2 0.0 0.2 0.4 0.1 0.2 0.3 h DFM =0.0067 =1.0 /L (h-h DFM )L 1/ L=24 L=32 L=40 h BG =0.67 =1.0 /L (h-h BG )L 1/ Figure 4.15: Finite size scaling of the correlation length in the vicinities of criticalh DFM (left) andh BG (right). The inter-dimer couplingJ ′ = J/4, and dilution concentration p = 0.125. We hereby estimate the correlation length through the momentum space structure factors by a second-moment method [15, 51], assumingS ⊥ ( ~ k) follows a Lorentz close to ~ k = (π,π): ξ = 1 Δk s S ⊥ (π,π) S ⊥ (π +Δk,π) −1 (4.25) where Δk = 2π L . The scaling form of the correlation length is ξ =Lf ξ [L 1/ν (h−h c )]. (4.26) 128 -2 0 2 0 2 4 6 -0.4 -0.2 0.0 0.2 0.4 0.5 1.0 1.5 -2 0 2 0 2 4 6 -0.4 -0.2 0.0 0.2 0.4 0.0 0.2 0.4 0.6 0.8 h BG =0.67 =1.0, =0.9 m s 2 L 2 / (h-h BG )L 1/ L=24 L=32 L=40 h DFM =0.0067 =1.0, =0.9 m s 2 L 2 / (h-h DFM )L 1/ h BG =0.67 =1.0, z =2.0 (h-h BG )L 1/ s L d+z-2 L=24 L=32 L=40 h DFM =0.0067 =1.0, z =2.0 (h-h DFM )L 1/ s L d+z-2 Figure 4.16: Finite size scaling of the staggered magnetization and the spin stiffness in the vicinities of critical fieldsh DFM (left) andh BG (right). The model parameters are the same as in the caption of Fig. 4.15. It defines the correlation exponentν. By studying finite-size scaling of correlation length in the vicinity of two QCPs, we estimate the two critical fields to beh DFM = 0.007(1), andh BG = 0.67(1) respectively. The QMC results of scaling of correlation length are shown in Fig. 4.15. For each transition the data of different sizes collapse well in the ordered phase. Remarkably we find that the scaling plots ofξ at two QPTs give us the same exponentν≈ 1.0(1). 129 We then discuss the scaling of staggered magnetization and spin stiffness. The finite-size scaling forms are m s =L −β/ν f ms [L 1/ν (h−h c )]; (4.27) ρ s =L −(d+z−2) f ρs [L 1/ν (h−h c )]. (4.28) whereβ is the order parameter exponent,z is the dynamical exponent, andd = 2 is the spacial dimension of the system. We present our results in Fig. 4.16. From the plot we independently estimate the critical fields to beh DFM = 0.007(1) andh BG = 0.67(1), which confirm our previously estimated values from the scaling of correlation length. For the correlation exponent we getν = 1.0(1) at bothh DFM andh BG from bothm s andρ s data in Fig. 4.16. These results are also consistent with the extracted value ofν from the scaling plots of the correlation length. The extracted value of order parameter β is:β = 0.9(1) at bothh DFM andh BG . Moreover the scaling ofρ s allows us to extract the dynamical exponent of the transition. We obtainz = 2.0(1) at bothh DFM andh BG from the lower plots in Fig. 4.16. The transition ath BG can be mapped to a BG-to-SF transition [76, 77]. z = d = 2 is theoretically predicted for this transition [22]. Some bounds of other exponents has also been given by the theory [22]. The correlation length exponentν must satisfy the Harris 130 criterion [12]ν> 2/d = 1. The exponentη, which is associated to the order parameter correlation function, is related withβ through 2β/ν =d+z−2+η. (4.29) It should satisfyη6 2−d = 0. We estimateη through Eq. 4.29 by using other expo- nents from our result, and getη =−0.2(2). Together with directly estimatedz = 2.0(1) andν = 1.0(1), we find our results of the exponents are consistent with these theoretical predictions. The consistence shows that our mapping from the field-induced transition in site-diluted antiferromagnets to the BG-to-SF transition in interacting bosons is appro- priate. Furthermore we notice that all the exponents from our result are in full agreement with the exponents estimated from the same transition in a bilayer model [76]. Since no further theoretical prediction is given on the exact value of exponents exceptz, our numerical finding of the exponentsβ andν enriches the theory on the BG-to-SF transi- tion. For the transition ath DFM no information is given from previous studies. But the above estimated values of critical exponents from our QMC results support the following argu- ment: the two phase transitions ath DFM andh BG belong to the same universality class of a 2D BG-to-SF QPT identified by ν = 1.0(1), β = 0.9(1), and z = 2.0(1). This confirms our conclusion of taking the DFM phase as a low-field BG phase. 131 -1 0 1 0.1 0.2 0.3 /L h DFM =0.045 =1.0 -2 0 2 0.0 0.1 0.2 0.3 0.4 0.5 h DFM =0.56 =1.0 L=24 L=32 L=40 -1 0 1 1 2 3 m s 2 L 2 / h DFM =0.045 =1.0, =0.9 -2 0 2 0 2 4 h DFM =0.56 =1.0, =0.9 L=24 L=32 L=40 -1 0 1 0 2 4 s L d+z-2 (h-h DFM )L 1/ h DFM =0.045 =1.0, z=2.0 -2 0 2 0 2 4 6 8 h DFM =0.56 =1.0, z=2.0 (h-h BG )L 1/ L=24 L=32 L=40 Figure 4.17: From top to bottom: finite size scaling of correlation lengthξ, staggered magnetization m s , and spin stiffness ρ s at critical fields h DFM = 0.05(1) (left) and h BG = 0.56(1) (right) in the site-diluted coupled dimers withJ ′ /J = 0.35 andp = 1/8. Moreover we have performed quantum scaling of correlation length, staggered magneti- zation, and spin stiffness on all the other points on the critical curve shown in Fig. 4.14. The extracted critical exponents all fall into the universality class of the proposed 2D BG-to-SF transition, with ν = 1.0(1), β = 0.9(1), and z = 2.0(1). As an example, the scaling plots of the correlation length, staggered magnetization, and spins stiffness ath DFM andh BG forJ ′ /J = 0.35 are shown in Fig. 4.17. One does see the agreement 132 of exponents with the extracted ones from Fig. 4.15 and Fig. 4.16. Based on these data we identify this gapless quantum disordered phase as the BG phase. It develops into a lower-field DFM phase and a higher-field BG phase in the weak coupling limit. 4.6 Summary We conclude in this section as follows. In this chapter we systematically investigated the phase diagram of 2D site-diluted coupled dimer system in the presence of a mag- netic field with QMC simulations. The phase diagram is much richer than the one of the clean system. In the phase diagram there is a long-range antiferromagnetically ordered phase, and a gapped disordered phase characterized by a magnetization plateau. These two phases are separated by a gapless quantum disordered phase. The phase transi- tion between the gapless quantum disordered phase and the long-range ordered phase belongs to the same universality class of BG-to-SF transition in 2D. The predicted value of dynamical exponentz = d = 2 is verified by our finite-size scaling of the spin stiff- ness. We also determined the values of other exponents through finite-size scaling of the staggered magnetization and the correlation length. The properties at phase tran- sitions and in the quantum disordered phases suggest that we can map the site-diluted spin system into a random interacting bosonic system, so that the three phases in the phase diagram corresponds to the superfluid, Bose glass, and Mott insulating phases, 133 respectively. The BG phase develops to two branches in the weakly coupled dimer systems: a lower field DFM phase, and a higher BG phase. The physics in these two phases is completely different. In the higher field BG the dominating physics is the ran- dom localization of triplets; while the lower field DFM is caused by the inhomogeneous polarization of doping-induced local moments. We focused on studying the magnetic properties in the DFM phase by mapping the site- diluted coupled dimer model to a low-energy effective spin model. The critical field h DFM is found to be extremely low. The magnetization process in the DFM phase is highly inhomogeneous. Numerically we observe a series of magnetization PPs, which means a gapless region with exponentially small uniform susceptibility. We propose a quantum scenario of the DFM phase, based on formation of local clusters through effective coupling close toJ ′ , to explain the inhomogeneous magnetization process and the formation of the PPs. We then study the distribution of the local magnetization and the bond energies to verify the proposed scenario. In both the effective spin model and the site-diluted model we confirm the formation of local singlets with effective coupling close toJ ′ , which is a fundamental argument of the scenario. By further performing a spin-boson transformation, and the statistics on the local clusters we successfully repro- duce the nontrivial uniform magnetization curve, which is in agreement with our QMC result. We give the relation of the distribution of the local cluster gaps and the uniform susceptibility, which allow us to quantitatively determine the magnetization values and the field positions of the PPs. All above finding are completely contrary to the picture of 134 continuous polarization in classical spin model, indicating that the DFM phase is a pure quantum many-body effect. 135 Chapter 5 Conclusions In this thesis we investigated quantum phase transitions in disordered Heisenberg anti- ferromagnets. We observed that the interplay between quantum fluctuations and geo- metric randomness can dramatically change the ground state properties of quantum spin systems. The long-range ordered phase in these systems is typically highly inhomoge- neous, and sensitive to perturbations. When it is perturbed, the universality class of the order-disorder transition can be different from the clean case, and the transition is often accompanied by the appearance of novel quantum disordered phases. Here we showed these general properties by means of two nontrivial examples in two-dimensional quan- tum Heisenberg models. In the first example, we studied quantum phase transitions in an inhomogeneous two- dimensional antiferromagnet driven by bond dilution. The system was found to experi- ence a percolation transition. The case of weak inhomogeneity was found to be irrele- vant, and the associate transition takes the nature of classical percolation, i.e., the mag- netic transition coincides with the classical percolation, and the critical exponents are 136 found to be exactly the same as the classical ones. However, we also found that moderate to strong inhomogeneity becomes relevant. In this case the magnetic transition deviates from the classical one beyond two multi-critical points, manifesting quantum nature. We found that for moderate inhomogeneity the magnetic transition can be understood within a proposed quantum percolation scenario as a renormalized percolative transition. The critical exponents continuously deviate from the classical ones. Interestingly we find a novel quantum disordered phase in between the quantum and classical critical curves. In this phase an infinite percolating cluster is present, but no long-range antiferromagnetic order exists. The finite temperature uniform susceptibility data on percolating clusters show a clear non-universal power-law divergent, implying that this phase is gapless and has a Griffiths singularity. We confirmed that this phase is a novel quantum Griffiths phase in two-dimensional Heisenberg antiferromagnets by performing statistics on the local energy gaps. In the second example we investigated the magnetic properties of a two-dimensional site-diluted spin-gap system. We carefully studied the field response of this system. The phase diagram was found to be very rich, including several magnetically ordered phases separated by a series of quantum disordered phases. In the weak coupling limit, we found that an extremely small but finite field may destroy the order-by-disorder phase, driving the system into a quantum disordered disordered-free-moment phase. The disordered-free-moment phase has short-range spin-spin correlations and a gapless 137 spectrum. The magnetization process is highly inhomogeneous. The uniform suscep- tibility is always finite but can be exponentially small, leading to a series of magnetic pseudo-plateaus in the magnetization curve. We proposed a quantum magnetization picture in this phase based on the formation of local clusters of the free moments. This picture was verified by studying the distribution of local magnetization and bond ener- gies. With the aid of a spin-boson mapping, it also gives a quantitative explanation of the pseudo-plateaus in the uniform magnetization curve. The order-by-disorder-to- disordered-free-moment phase transition was also investigated. It was found that this transition belongs to the same universality class of a Bose-glass-to-superfluid transition in two dimension. With the spin-boson mapping, it reveals that the disordered-free- moment phase is another low-field Bose glass phase. Our study of the general phase diagram of the coupled dimer system confirms this conclusion. 138 Bibliography [1] AFFLECK, I. Phys. Rev. B 43 (1991), 3215. [2] ALET, F., WESSEL, S., AND TROYER, M. Phys. Rev. E 71 (2005), 036706. [3] ARFKEN, G. Mathematical Methods for Physicists, 3rd ed. Academic Press, Orlando, 1985. [4] AZUMA, M., FUJISHIRO, Y., TAKANO, M., NOHARA, M., AND TAKAGI, H. Phys. Rev. B 55 (1997), R8658. [5] BEARD, B. B., AND WIESE, U.-J. Phys. Rev. Lett. 77 (1996), 5130. [6] BRAY-ALI, N., AND MOORE, J. E. Phys. Rev. B 69 (2004), 184505. [7] CARLSON, J. Phys. Rev. B 40 (1989), 846. [8] CARRETTA, P., RIGAMONTI, A., AND SALA, R. Phys. Rev. B 55 (1997), 3734. [9] CHABOUSSANT, G., JULIEN, M.-H., FAGOT-REVURAT, Y., HANSON, M., L ´ EVY, L. P., BERTHIER, C., HORVATI ´ C, M., AND PIOVESANA, O. Eur. J. Phys. B 6 (1998), 167. [10] CHAKRAVARTY, S., HALPERIN, B. I., AND NELSON, D. R. Phys. Rev. Lett. 60 (1988), 1057. [11] CHAKRAVARTY, S., AND STEIN, D. B. Phys. Rev. Lett. 49 (1982), 582. [12] CHAYES, J. T., CHAYES, L., FISHER, D. S., , AND SPENCER, T. Phys. Rev. Lett. 57 (1986), 2999. 139 [13] CHEN, K., FERRENBERG, A. M., AND LANDAU, D. P. Phys. Rev. B 48 (1993), 3249. [14] CHERNYSHEV, A. L., CHEN, Y. C., AND NETO, A. H. C. Phys. Rev. B 65 (2002), 104407. [15] COOPER, F., FREEDMAN, B., AND PRESTON, D. Nucl. Phys. B 210 (1982), 210. [16] CORTI, M., RIGAMONTI, A., TABAK, F., CARRETTA, P., LICCI, F., AND RAFFO, L. Phys. Rev. B 52 (1995), 4226. [17] DAGOTTO, E., AND RICE, T. M. Science 271 (1996), 618. [18] DE GENNES, P. G. J. Phys. (Paris) Lett. 37 (1976), L1. [19] DORNEICH, A., AND TROYER, M. Phys. Rev. E 64 (2001), 066701. [20] EVERTZ, H. G., LANA, G., AND MARCU, M. Phys. Rev. Lett. 70 (1993), 875. [21] FISHER, D. S. Phys. Rev. Lett. 69 (1992), 534. [22] FISHER, D. S. Phys. Rev. B 50 (1994), 3799. [23] FISHER, D. S. Phys. Rev. B 51 (1995), 6411. [24] FISHER, M., WEICHMAN, P. B., GRINSTEIN, G., AND FISHER, D. S. Phys. Rev. B 40 (1989), 546. [25] FUKUYAMA, H., NAGAOSA, N., SAITO, M., AND TANIMOTO, T. J. Phys. Soc. Jpn. 65 (1996), 2377. [26] GENNES, P. G. D. La Recherche (Paris) 7 (1976), 919. [27] GIAMARCHI, T., AND TSVELIK, A. M. Phys. Rev. B 59 (1999), 11398. [28] GREVEN, M., AND BIRGENEAU, R. J. Phys. Rev. Lett. 81 (1998), 1945. [29] GREVEN, M., AND BIRGENEAU, R. J. Phys. Rev. Lett. 81 (1998), 1945. [30] GRIFFITHS, R. B. Phys. Rev. Lett. 23 (1969), 17. [31] HANDSCOMB, D. C. Proc. Cambridge Philos. Soc. 58 (1962), 594. [32] HENELIUS, P., AND SANDVIK, A. W. Phys. Rev. B 62 (2000), 1102. [33] HONDA, Z., ASAKAWA, H., AND KATSUMATA, K. Phys. Rev. Lett. 81 (1998), 2566. 140 [34] HYMAN, R. A., AND YANG, K. Phys. Rev. Lett. 78 (1997), 1783. [35] HYMAN, R. A., YANG, K., BHATT, R. N., AND GIRVIN, S. M. Phys. Rev. Lett. 76 (1996), 839. [36] IGARASHI, J., TONEGAWA, T., KABURAGI, M., AND FULDE, P. Phys. Rev. B 51 (1995), 5814. [37] IGL ´ OI, F., AND MONTHUS, C. Phys. Rep. 412 (2005), 277. [38] JAIME, M., CORREA, V. F., HARRISON, N., BATISTA, C. D., KAWASHIMA, N., KAZUMA, Y., JORGE, G. A., STERN, R., HEINMAA, I., ZVYAGIN, S. A., SASAGO, Y., AND UCHINOKURA, K. Phys. Rev. Lett. 93 (2004), 087203. [39] KATO, K., TODO, S., HARADA, K., KAWASHIMA, N., MIYASHITA, S., AND TAKAYAMA, H. Phys. Rev. Lett. 84 (2000), 4204. [40] KAWASHIMA, N. J. Phys. Soc. Jpn. 73 (2004), 3219. [41] LAFLORENCIE, N., POILBLANC, D., AND SANDVIK, A. W. Phys. Rev. B 69 (2004), 212412. [42] LAFLORENCIE, N., POILBLANC, D., AND SIGRIST, M. Phys. Rev. B 71 (2005), 212403. [43] LEE, D. H., JOANNOPOULOS, J. D., AND NEGELE, J. W. Phys. Rev. B 30 (1984), 1599. [44] LYKLEMA, J. W. Phys. Rev. Lett. 49 (1982), 88. [45] MA, S. K., DASGUPTA, C., AND HU, C. K. Phys. Rev. Lett. 43 (1979), 1434. [46] MANOUSAKIS, E. Rev. Mod. Phys. 63 (1991), 1. [47] MASUDA, T., UCHINOKURA, K., HAYASHI, T., AND MIURA, N. Phys. Rev. B 66 (2002), 174416. [48] MASUDA, T., ZHELUDEV, A., UCHINOKURA, K., CHUNG, J.-H., AND PARK, S. Phys. Rev. Lett. 93 (2004), 077206. [49] MATSUMOTO, M., NORMAND, B., RICE, T. M., AND SIGRIST, M. Phys. Rev. Lett. 89 (2002), 077203. [50] MATSUMOTO, M., NORMAND, B., RICE, T. M., AND SIGRIST, M. Phys. Rev. B 69 (2004), 054423. 141 [51] MATSUMOTO, M., YASUDA, C., TODO, S., AND TAKAYAMA, H. Phys. Rev. B 65 (2002), 014407. [52] MCCOY, B. M. Phys. Rev. Lett. 23 (1969), 383. [53] METROPOLIS, N., ROSENBLUTH, A., ROSENBLUTH, M., TELLER, A. H., AND TELLER, E. J. Chem. Phys. 21 (1953), 1087. [54] MIKESKA, H.-J., GHOSH, A., AND KOLEZHUK, A. K. Phys. Rev. Lett. 93 (2004), 217204. [55] MIKESKA, H.-J., NEUGEBAUER, U., AND SCHOLLW ¨ OCK, U. Phys. Rev. B 55 (1997), 2955. [56] MISGUICH, G., AND OSHIKAWA, M. J. Phys. Soc. Jpn. 73 (2004), 3429. [57] MOUKOURI, S., CHEN, L., AND CARON, L. G. Phys. Rev. B 53 (1996), R488. [58] MUCCIOLO, E. R., NETO, A. H. C., AND CHAMON, C. Phys. Rev. B 69 (2004), 214424. [59] N. LAFLORENCIE, S. WESSEL, A. L. H. R. Phys. Rev. B 73 (2006), 060403(R). [60] NAKAMURA, T. Phys. Rev. B 59 (1999), R6589. [61] NAKAYAMA, T., AND YAKUBO, K. Fractal concepts in condensed matter physics. [62] NEWMAN, M. E. J., AND ZIFF, R. M. Phys. Rev. Lett. 85 (2000), 4104. [63] NEWMAN, M. E. J., AND ZIFF, R. M. Phys. Rev. E 64 (2001), 016706. [64] NIKUNI, T., OSHIKAWA, M., OOSAWA, A., AND TANAKA, H. Phys. Rev. Lett. 84 (2000), 5868. [65] NOHADANI, O., WESSEL, S., AND HAAS, S. Phys. Rev. B 72 (2005), 024440. [66] NOHADANI, O., WESSEL, S., AND HAAS, S. Phys. Rev. Lett. 95 (2006), 227201. [67] NORMAND, B., AND MILA, F. Phys. Rev. B 65 (2002), 104411. [68] OOSAWA, A., FUJISAWA, M., KAKURAI, K., AND TANAKA, H. Phys. Rev. B 67 (2003), 184424. [69] OOSAWA, A., AND TANAKA, H. Phys. Rev. B 65 (2002), 184437. 142 [70] OOSAWA, A., AND TANAKA, H. Phys. Rev. B 65 (2002), 184437. [71] PICH, C., YOUNG, A. P., RIEGER, H., AND KAWASHIMA, N. Phys. Rev. Lett. 81 (1998), 5916. [72] POLLOCK, E. L., AND CEPERLEY, D. M. Phys. Rev. B 36 (1987), 8343. [73] POLLOCK, E. L., AND CEPERLEY, D. M. Phys. Rev. B 36 (1987), 8343. [74] RICE, T. Science 298 (2002), 760. [75] RIEGER, H., AND YOUNG, A. P. Phys. Rev. B 54 (1996), 3328. [76] ROSCILDE, T. Phys. Rev. B 74 (2006), 144418. [77] ROSCILDE, T., AND HAAS, S. Phys. Rev. Lett. 95 (2005), 207206. [78] R ¨ UEGG, C., CAVADINI, N., FURRER, A., G ¨ UDEL, H.-U., KRMER, K., MUTKA, H., WILDES, A., HABICHT, K., AND VORDERWISCH, P. Nature 423 (2003), 62. [79] SACHDEV, S. Quantum Phase Transitions. Cambridge University Press, Cam- bridge, 1999. [80] SANDVIK, A. W. J. Phys. A 25 (1992), 3667. [81] SANDVIK, A. W. Phys. Rev. B 56 (1997), 11678. [82] SANDVIK, A. W. Phys. Rev. B 59 (1999), R14157. [83] SANDVIK, A. W. Phys. Rev. B 66 (2002), 024418. [84] SANDVIK, A. W. Phys. Rev. Lett. 89 (2002), 177201. [85] SANDVIK, A. W. Phys. Rev. E 68 (2003), 056701. [86] SEBASTIAN, S. E., SHARMA, P. A., JAIME, M., HARRISON, N., CORREA, V., BALICAS, L., KAWASHIMA, N., BATISTA, C. D., AND FISHER, I. R. Phys. Rev. B 72 (2005), 100404. [87] SENGUPTA, P., SANDVIK, A. W., AND CAMPBELL, D. K. Phys. Rev. B 65 (2002), 155113. [88] SENTHIL, T., AND SACHDEV, S. Phys. Rev. Lett. 77 (1996), 5292. [89] SHENDER, E. F., AND KIVELSON, S. Phys. Rev. Lett. 66 (1991), 2384. 143 [90] SIERRA, G., MARTIN-DELGADO, M. A., WHITE, S., SCALAPINO, D., AND DUKELSKY, J. Phys. Rev. B 59 (1999), 7973. [91] SIGRIST, M., AND FURUSAKI, A. J. Phys. Soc. Jpn. 65 (1996), 2385. [92] SKAL, A. S., AND SHKLOVSKII, B. I. Sov. Phys. Semicond. 8 (1975), 1029. [93] SKNEPNEK, R., VOJTA, T., AND VOJTA, M. Phys. Rev. Lett. 93 (2004), 097201. [94] SØRENSEN, E. S., AFFLECK, I., AUGIER, D., AND POILBLANC, D. Phys. Rev. B 58 (1998), R14701. [95] STAUFFER, D., AND AHARONY, A. Introduction to Percolation Theory. Taylor and Francis, London, 1994. [96] SUZUKI, M. Prog. Theor. Phys. 56 (1976), 1454. [97] SWENDSEN, R. H., AND WANG, J. Phys. Rev. Lett. 58 (1987), 86. [98] SYLJU ˚ ASEN, O. F. Phys. Rev. E 67 (2003), 046701. [99] SYLJU ˚ ASEN, O. F., AND SANDVIK, A. W. Phys. Rev. E 66 (2002), 046701. [100] TODO, S., KATO, K., AND TAKAYAMA, H. J. Phys. Soc. Jpn. Suppl. 69 A (2000), 355. [101] TODO, S., YASUDA, C., KATO, K., HARADA, K., KAWASHIMA, N., MIYASHITA, S., AND TAKAYAMA, H. Prog. Theor. Phys. Suppl. 138 (2000), 507. [102] TROTTER, H. F. Proc. Am. Math. Soc. 10 (1959), 545. [103] TSUJII, H., HONDA, Z., ANDRAKA, B., KATSUMATA, K., AND TAKANO, Y. Phys. Rev. B 71 (2005), 014426. [104] UCHIYAMA, Y., SASAGO, Y., TSUKADA, I., UCHINOKURA, K., ZHELUDEV, A., HAYASHI, Y., MIURA, N., AND B ¨ ONI, P. Phys. Rev. Lett. 83 (1999), 632. [105] VAJK, O. P., AND GREVEN, M. Phys. Rev. Lett. 89 (2002), 177202. [106] VAJK, O. P., MANG, P. K., GREVEN, M., GEHRING, P. M., AND LYNN, J. W. Science 295 (2002), 1691. [107] WATSON, B. C., KOTOV, V. N., MEISEL, M. W., HALL, D. W., GRANROTH, G. E., MONTFROOIJ, W. T., NAGLER, S. E., JENSEN, D. A., BACKOV, R., PETRUSKA, M. A., FANUCCI, G. E., AND TALHAM, D. R. Phys. Rev. Lett. 86 (2001), 5168. 144 [108] WEBER, H., AND VOJTA, M. Eur. Phys. J. B 53 (2006), 185. [109] WESSEL, S., NORMAND, B., AND HAAS, S. Phys. Rev. B 69 (2004), 220402. [110] WESSEL, S., NORMAND, B., SIGRIST, M., AND HAAS, S. Phys. Rev. Lett. 86 (2001), 1086. [111] WESSEL, S., OLSHANII, M., AND HAAS, S. Phys. Rev. Lett. 87 (2001), 206407. [112] WHITE, S. R., NOACK, R. M., AND SCALAPINO, D. J. Phys. Rev. Lett. 73 (1994), 886. [113] WIESE, U.-J., AND YING, H.-P. Z. Phys. B 93 (1994), 147. [114] WOLFF, U. Phys. Rev. Lett. 62 (1989), 361. [115] WOLFF, U. Phys. Lett. B 222 (1989), 473. [116] YANG, K., HYMAN, R. A., BHATT, R. N., AND GIRVIN, S. M. J. Appl. Phys. 79 (1996), 5096. [117] YASUDA, C., TODO, S., HARADA, K., KAWASHIMA, N., MIYASHITA, S., AND TAKAYAMA, H. Phys. Rev. B 63 (2001), R140415. [118] YASUDA, C., TODO, S., MATSUMOTO, M., AND TAKAYAMA, H. Phys. Rev. B 64 (2001), 092405. [119] YU, R., ROSCILDE, T., AND HAAS, S. cond-mat/0611615. [120] YU, R., ROSCILDE, T., AND HAAS, S. Phys. Rev. Lett. 94 (2005), 197204. [121] ZHELUDEV, A. cond-mat/0507534. [122] ZHELUDEV, A., HONDA, Z., KATSUMATA, K., FEYERHEIM, R., AND PROKES, K. Europhys. Lett. 55 (2001), 868. 145

Abstract (if available)

Abstract Recently quantum phase transitions have attracted the interest of both theorists and experimentalists in condensed matter physics. Quantum magnets provide a perfect playground for studying these phase transitions since they can be triggered by many control parameters such as frustration, lattice dimerization, and magnetic field. Most previous studies have focused on the magnetic properties in pure systems. In these systems, responses to the triggering parameters are found to be uniform, leading to homogeneous phases. However little progress has been made so far on the phase transitions and properties in disordered quantum magnets because they are more complicated systems, and few theoretical tools can be applied.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Phase diagram of disordered quantum antiferromagnets

PDF

Properties of hard-core bosons in potential traps

PDF

Disordered quantum spin chains with long-range antiferromagnetic interactions

PDF

Modeling nanodevices: from semiconductor heterostructures to Josephson junction arrays

PDF

Entanglement in strongly fluctuating quantum many-body states

PDF

Entanglement in quantum critical and topological phases

PDF

One-dimensional nanomaterials: synthesis and applications

PDF

Magnetization study of two dimensional helium three

PDF

Properties of magnetic nanostructures

PDF

Topics in quantum information and the theory of open quantum systems

PDF

Synthesis and application of one-dimensional nanomaterials

PDF

Electronic correlation effects in multi-band systems

PDF

Dynamical error suppression for quantum information

PDF

Out-of-equilibrium dynamics of inhomogeneous quantum systems

PDF

Modeling graphene: magnetic, transport and optical properties

PDF

Low field study of A-like phase for superfluid 3He in aerogel

PDF

Topological protection of quantum coherence in a dissipative, disordered environment

PDF

Entanglement parity effects in quantum spin chains

PDF

Optoelectronic properties and device physics of individual suspended carbon nanotubes

PDF

Healing of defects in random antiferromagnetic spin chains

Asset Metadata

Creator Yu, Rong (author)

Core Title Quantum phase transitions in disordered antiferromagnets

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Physics

Publication Date 07/24/2009

Defense Date 05/01/2007

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag disordered system,OAI-PMH Harvest,percolation,quantum magnets,quantum Monte Carlo,spin glass

Language English

Advisor Haas, Stephan W. (committee chair), Bickers, Nelson Eugene, Jr. (committee member), Dappen, Werner (committee member), Kresin, Vitaly V. (committee member), Zhou, Chongwu (committee member)

Creator Email rongyu@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-m652

Unique identifier UC163190

Identifier etd-Yu-20070724 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-523066 (legacy record id),usctheses-m652 (legacy record id)

Legacy Identifier etd-Yu-20070724.pdf

Dmrecord 523066

Document Type Dissertation

Rights Yu, Rong

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu