ebook img

Efficient evaluation of partition functions of frustrated and inhomogeneous spin systems PDF

0.22 MB·English
by  V. Murg
Save to my drive
Quick download
Download
Most books are stored in the elastic cloud where traffic is expensive. For this reason, we have a limit on daily download.

Preview Efficient evaluation of partition functions of frustrated and inhomogeneous spin systems

Efficient evaluation of partition functions of frustrated and inhomogeneous spin systems V. Murg1, F. Verstraete1,2, J. I. Cirac1 1Max-Planck-Institut fu¨r Quantenoptik, Hans-Kopfermann-Str. 1, Garching, D-85748, Germany 2Institute for Quantum Information, Caltech, Pasadena, US (Dated: February 2, 2008) 5 0 We present a numerical method to evaluate partition functions and associated correlation func- 0 tions of inhomogeneous 2–D classical spin systems and 1–D quantum spin systems. The method is 2 scalableandhasacontrollederror. Weillustratethealgorithmbycalculatingthefinite–temperature properties of bosonic particles in 1–D optical lattices, as realized in current experiments. n a PACSnumbers: 03.67.-a,03.75.Lm,02.70.-c,75.40.Mg J 0 2 The study of thermodynamic properties of 2–D classi- of strongly interacting bosonic atoms in optical lattices, cal spin systems has a very long and fruitful history. On a problem which has attracted a lot of interest in the ] r theonehand,ithasprovidedvaluableinsightsinthethe- atomic physics community in the last few years due to e h oryofmagnetismandphasetransitions. Ontheother,it therecentexperimentalachievements[12,13,14]. Inthis t hasallowedustodescribe1–Dquantumsystemsatfinite system, 20-100atoms are trapped by the combinationof o temperature by virtue of the Suzuki–Trotter decomposi- a periodic and a harmonic (i.e. inhomogeneous) poten- . at tion. tialcreatedby lasersin1–Dandat a finite temperature. m Asveryfewspinmodelsareexactlysolvable,manydif- In the so–called Tonks–Girardeau limit [15], the prob- - ferent approximate methods have been proposed to cal- lem can be exactly solved via fermionization [13], and d culate the associated partition functions. Monte Carlo thus this provides us with a reliable benchmark for our n seems to be the method of choice for non–frustrated method. Outside this limit, we are able to reproduce o systems, but fails in the description of frustrated and certain features experimentally observed [12]. c [ fermionic systems due to the notorious sign problem [1]. Our method relies in reexpressing the partition and Nosuchsignproblemoccursinthemethoddevelopedby correlation functions as a contraction of a collection of 1 Nishino [2], where the largest eigenvalue of the transfer 4-index tensors, which are disposed according to a 2–D v 3 matrix of the classical spin model can be approximated configuration. We will perform this task for both 2–D 9 byusingavariationofthedensitymatrixrenormalization classical and 1–D quantum systems. We will then show 4 group (DMRG) approach [3, 4]. By making use of the how the methods introduced in [16] can be used to ap- 1 Suzuki–Trotterdecomposition,thismethodhasalsobeen proximate these contractions in a controlled way, and 0 used to calculate the free energy of translational invari- thus lead to a scalable algorithm for the evaluation of 5 0 ant 1-D quantum systems [5, 6, 7]. The main restriction the quantities of interest. / of this method is that it cannot be applied in situations 1.- 2–D classical systems: Letusconsiderfirstthepar- t a in which the number of particles is finite and/orthe sys- titionfunctionofaninhomogeneousclassical2–Dn–level m tem is not homogeneous. The method may also become spin system on a L ×L lattice. For simplicity we will 1 2 - illconditionedwhenthetransfermatrixisnothermitian. concentrate on a square and nearest–neighbor interac- d Finally,inthecaseof1-Dquantumsystems,arecentde- tions, although our method can be easily extended to n o velopment[8,9]allowsustoovercometheseproblemsby other short–range situations. We have c extending the concept of matrix product states [10, 11] v: to mixed states, e.g. by using the idea of purification of Z = exp −βH(x11,...,xL1L2) , Xi states [8]. This method is, however, specially designed x11,.X..,xL1L2 (cid:2) (cid:3) for 1–D quantum systems, and cannot be extended to r a classical 2–D models. where Here we take a completely different approach which allows us to overcome the drawbacks of the above men- H x11,... = Hij xij,xi+1,j +Hij xij,xi,j+1 ↓ → tioned methods in both 2–D classical and 1–D quantum (cid:0) (cid:1) Xij h (cid:0) (cid:1) (cid:0) (cid:1)i systems. We achieve this by evaluating the associated partition and correlation functions directly. The main is the Hamiltonian, xij = 1,...,n and β is the inverse advantages of this method are that the approximations temperature. The singular value decomposition allows made are very well controlled, that it applies to frus- us to write trated, inhomogeneous and finite classical and quantum n systems, and that it can be generalizedto higher dimen- exp −βHij(x,y) = fij(x)gij (y), q qα qα sions. We will illustrate its performancewith the system α=1 (cid:2) (cid:3) X 2 with q ∈{↓,→}. Defining the tensors conditions in the vertical direction. A general expecta- tionvalueofanoperatoroftheformO=ZO1⊗···⊗ON n Xij = fij(x)gi−1,j(x)fij (x)gi,j−1(x), can also be reexpressed as a contraction of tensors with lrud ↓d ↓u →r →l the same structure: it is merely required to replace each x=1 X tensor X1j in the first row by thepartitionfunctioncannowbecalculatedbycontract- ing all 4-index tensors Xij arranged on a square lattice X1j Oj = OjTˆjSˆj Tˆj Sˆj . (ll′)(rr′)ud 1l 1r 2l′ 2r′ in such a way that, e.g., the indices l,r,u,d of Xij are [ud] (cid:0) (cid:1) h i contractedwiththeindicesr,l,d,uofthe respectiveten- 3.- Tensor contraction: In the following, we adapt the sorsXi,j−1,Xi,j+1,Xi−1,j,Xi+1,j. Inordertodetermine algorithmintroducedin[17]inordertocontracttheten- the expectation value of a general operator of the form sorsXij introducedaboveinacontrolledway. Themain O({xij}) = Z Oij(xij), one just has to replace each ij idea is to express the objects resulting from the contrac- tensor Xij by Q tionoftensorsalongthefirstandlastcolumninthe 2–D configurationas matrix product states (MPS) and those n Xij Oij = Oij(x)fij(x)gi−1,j(x)fij (x)gi,j−1. obtainedalongthecolumns2,3,...,L−1asmatrixprod- lrud ↓d ↓u →r →l uct operators (MPO) [8]. More precisely, we define x=1 (cid:0) (cid:1) X 2.- 1–D quantum systems: We consider the partition m functionofaninhomogeneous1–Dquantumsystemcom- hX1 | := Tr X1r11...XMrM1 hr1...rM| posed of L n-level systems, r1..X.rM=1 (cid:0) (cid:1) m Z =Trexp(−βH). | XLi := Tr X1L...XML |l ...l i l1 lM 1 M It is always possible to write the Hamiltonian H as a l1..X.lmM=1 (cid:0) (cid:1) sum H = kHk with each part consisting of a sum of Xj := Tr X1l1jr1...XMlMjrM |l1...ihr1... |, commuting terms. Let us, for simplicity, assume that P l1,rX1,...=1 (cid:16) (cid:17) H = H +H and that only local and 2-body nearest 1 2 neighbor interactions occur, i.e. H = Oi,i+1 and where m = n for 2–D classical systems and m = κ2 k i k Oi,i+1,Oj,j+1 = 0, with i,j = 1,...,L. The more for 1–D quantum systems. These MPS and MPOs are k k P associated to a chain of M m–dimensional systems and hgeneral case cain be treated in a similar way. Let us now their virtual dimension amounts to D = n. Note that consider a decomposition for 2–D classical systems the first and last matrices un- κ der the trace in the MPS and MPO reduce to vectors. β exp − Oi,i+1 = Sˆi ⊗Tˆi+1. (1) The partition function (and similarly other correlation M k kα kα (cid:18) (cid:19) α=1 functions) reads Z = hX1 |X2···XL−1| XLi. Evaluat- X ingthisexpressioniterativelybycalculatingstepbystep The singular value decomposition guarantees the exis- hXj |:=hXj−1 |Xj for j =2,...,L−1fails because the tence of such an expression with κ ≤ n2. As we will see virtual dimension of the MPS hXj | increases exponen- later, a smart choice of H = H can typically de- k k tially with j. A way to circumvent this problem is to crease κ drastically. Making use of the Suzuki–Trotter P replace in each iterative step the MPS hXj | by a MPS formula [25] hX˜j | with a reduced virtual dimension D˜ that approx- M imates the state hXj | best in the sense that the norm β 1 Z =Tr exp − Hk +O δK := khXj |− hX˜j |k is minimized. Due to the fact M M Yk (cid:18) (cid:19)! (cid:20) (cid:21) that this cost function is multiquadratic in the variables it can be readily seen that the partition function can of the MPS, this minimization can be carried out very againbecalculatedbycontractingacollectionof4-index efficiently [8, 16, 17]; the exponential increase of the vir- tensors Xij defined as tual dimension can hence be prevented and the iterative evaluation of Z becomes tractable, such that an approx- Xij ≡ TˆjSˆj Tˆj Sˆj , imation to the partition function can be obtained from (ll′)(rr′)ud 1l 1r 2l′ 2r′ h i[ud] Z ≃ hX˜L−1 | XLi. The accuracy of this approximation wheretheindices(l,l′)and(r,r′)arecombinedtoyielda depends only on the choice of the reduced dimension D˜ singleindexthatmayassumevaluesrangingfrom1toκ2. and the approximation becomes exact for D˜ ≥ DL. As Note that now the tensors Xij and Xi′j coincide, and thenormδK canbecalculatedateachstep,D˜ canbein- that the indices u of the first and d of the last row have creaseddynamicallyifthe obtainedaccuracyisnotlarge to be contracted with each other as well, which corre- enough. In the worst case scenario, such as in the NP- spondstoaclassicalspinsystemwithperiodicboundary complete Ising spin glasses [18], D˜ will probably have to 3 grow exponentially in L for a fixed precision of the par- tition function. But in less pathological cases it seems that D˜ only has to grow polynomially in L; indeed, the success of the methods developed by Nishino [2] in the translationalinvariantcaseindicatethatevenaconstant D˜ will produce very reliable results. 4.- Illustration: Bosons in optical lattices: A system of trapped bosonic particles in a 1–D optical lattice of L sites is described by the Bose-Hubbard Hamiltonian [19] L−1 L L U H =−J (a†a +h.c.)+ nˆ (nˆ −1)+ V nˆ , i i+1 2 i i i i FIG. 1: Density and (quasi)-momentum distribution in the i=1 i=1 i=1 Tonks-Girardeau gas limit, plotted for βJ = 1, L = 40, X X X where a†i and ai are the creation and annihilation oper- Nthe=nu2m1eraincadlVre0s/uJlts=fo0r.0D˜34=. 2Th(De˜d=ot8s)(carnodsstehs)e rseoplirdesleinnet ators on site i and nˆi = a†iai is the number operator. illustrates theexact results. From theinsets, theerror ofthe ThisHamiltoniandescribestheinterplaybetweentheki- numerical results can be gathered. netic energy due to the next-neighbor hopping with am- plitude J and the repulsive on-site interaction U of the particles. The last term in the Hamiltonian models the and Vˆ = e−2βMHRe−2βMHSe−2βM(HT−µNˆ). a˜†i, a˜i and n˜i harmonicconfinementofmagnitudeV =V (i−i )2. The thereby denote the projections of the creation, the an- i 0 0 variation of the ratio U/J drives a phase-transition be- nihilation and the number operators a†i, ai and ni on tweentheMott-insulatingandthesuperfluidphase,char- the q-particle subspace. The decomposition (1) of all acterized by localized and delocalized particles respec- two-particleoperatorsin expression(2) then straightfor- tively [20]. Experimentally, the variation of U/J can be wardly leads to a set of 4-index tensors Xlirjud, with in- realizedbytuningthedepthoftheopticallattice[19,21]. dices l and r ranging from 1 to (q +1)3 and indices u On the other hand, one typically measures directly the andd rangingfrom1 to q+1. Note that the typicalsec- momentumdistributionbylettingtheatomicgasexpand ondorderTrotterdecompositionwithH =Heven+Hodd and then measuring the density distribution. Thus, we would make the indices l and r range from 1 to (q+1)6. will be mainly interested here in the (quasi)–momentum Let us start out by considering the limit U/J → ∞ distribution in which double occupation of single lattice sites is pre- vented and the particles in the lattice form a Tonks– L 1 Girardeau gas [13]. In this limit, the Bose-Hubbard n = ha†a iei2πk(r−s)/L. k L r s HamiltonianmapstotheHamiltonianoftheexactlysolv- r,s=1 X able(inhomogeneous)XX-model,whichallowstobench- Our goal is now to study with our numerical method mark our algorithm. The comparison of our numerical the finite-temperature properties of this system for dif- results to the exact results can be gathered from fig. 1. ferent ratios U/J. We thereby assume that the system Here, the density and the (quasi)-momentum distribu- is in a thermal state corresponding to a grand canonical tion are considered for the special case βJ = 1, L = 40, ensemble with chemical potential µ, such that the par- N =21 and V /J =0.034. The numerical results shown 0 tition function is obtained as Z = Tre−β(H−µNˆ). Here, have been obtained for Trotter-number M =10 and two Nˆ = L nˆ represents the total number of particles. different reduced virtual dimensions D˜ = 2 and D˜ = 8. For the ni=u1meirical study, we assume a maximal particle– The norm δK was of order 10−4 for D˜ = 2 and 10−6 numbePr q per lattice site, such that we can project the for D˜ =8 [26]. From the insets, it can be gathered that Hamiltonian H on the subspace spanned by Fock-states the error of the numerical calculations is already very with particle-numbers per site ranging from 0 to q. The small for D˜ = 2 (of order 10−3) and decreases signifi- projected Hamiltonian H˜ then describes a chain of L cantly for D˜ =8. This errorcanbe decreasedfurther by spins,witheachspinactingona Hilbert-spaceofdimen- increasing the Trotter-number M. sionn=q+1. A Trotterdecompositionthatturnedout As the ratio U/J becomes finite, the system becomes to be advantageous for this case is physically more interesting, but lacks an exact mathe- matical solution. In order to judge the reliability of our e−β(H˜−µNˆ) = Vˆ†Vˆ M +O 1 , (2) numericalsolutionsinthiscase,wechecktheconvergence M2 with respect to the free parameters of our algorithm (q, (cid:16) (cid:17) (cid:20) (cid:21) D˜ and M). As an illustration, the convergence with re- with H˜ = H +H +H , H = −J L−1R(i)R(i+1), R S T R 2 i=1 spect to the parameter q is shown in figure 2. In this HS = −J2 iL=−11S(i)S(i+1), HT = PLi=1T(i), R(i) = figure, the density and the (quasi)-momentum distribu- a˜† +a˜ , S(i) = −i(a˜† −a˜ ), T(i) = 1n˜ (n˜ −1)+V n˜ tion are plotted for q = 2,3 and 4. We thereby assume i i P i i 2Pi i i i 4 FIG.2: Densityand(quasi)-momentumdistributionsforin- FIG. 3: FWHM of the (quasi)-momentum distribution as a teraction strengths U/J = 4 and 8. Here, βJ = 1, L = 40, function of U/J, calculated for temperatures βJ =0.5 (plus- N = 21 and M = 10. Numerical results were obtained sign),βJ =1(crosses) andβJ =2(dots). Thecorresponding for q =2 (plus-signs), q =3 (crosses) and q = 4 (solid line). (quasi)-momentum distributions for U/J = 2 and U/J = 8 For comparison, the distributions for U/J =0 (dotted lines) are illustrated in theplots at the right-handside. and U/J →∞ (dash-dotted lines) are also included. easily observed in present experiments. thatβJ =1,L=40andN =21andconsiderinteraction In summary, we have presented a numerical method strengths U/J = 4 and 8. The harmonic potential V is 0 for the investigation of thermal states of inhomogeneous choseninawaytodescribeRb-atomsinaharmonictrap 2–D classical and 1–D quantum systems. We have illus- offrequencyHz(alongthelinesof[13]). Wenotethatwe trated the usefulness of this method by applying it to a havetakenintoaccountthatchangesoftheratioU/J are system of trapped bosonic particles in a 1–D optical lat- obtained from changes in both the on-site interaction U tice – which is of current experimental interest. For this and the hopping amplitude J due to variations of the system, we have studied the error and the convergence depth of the optical lattice. The numerical calculations of the method. In addition, we have used the method have been performedwith M =10and D˜ =q+1. From to study some physical properties of this system and we the figure it can be gathered that convergence with re- have obtained results which can be verified experimen- spect to q is achieved for q ≥3. tally. We now useour numericalalgorithmtostudy aphysi- We thank M. A. Martin-Delgado and D. Porras for cal property of interacting bosons in an optical lattice, discussions. Work supported by the DFG, european namely the full width at half maximum (FWHM) of projects(ISTandRTN)andtheKompetenznetzwerkder the (quasi)-momentum distribution. It has been pre- BayrischenStaatsregierung Quanteninformation. dicted that the FWHM shows a kink at zero tempera- ture [22, 23, 24]. This kink is an indication for a Mott- superfluidtransition,sincetheFWHMisdirectlyrelated to the inverse correlation length. Experiments have also revealed this kink [12, 14]; they are, however, performed [1] M. Troyerand U.Wiese (2004), cond-mat/0408370. at finite temperature, something we can study with our [2] T. Nishino, J. Phys. Soc. Jpn. 64, 3598 (1995). [3] S. R.White, Phys. Rev.Lett 69, 2863 (1992). algorithm. In figure 3, we plot the numerical results for [4] U. Schollw¨ock (2004), cond-mat/0409292. the FWHM as a function of U/J for three different (in- [5] R. J. Bursill, T. Xiang, and G. A. Gehring, J. Phys. verse) temperatures βJ = 0.5,1 and 2. The physical Conden. Matter 8, L583 (1996). parameters L, N and V0 are thereby chosen as in the [6] N. Shibata, J. Phys. Soc. Jpn. 66, 2221 (1997). previouscase. The numericalresults have been obtained [7] X. Wang and T. Xiang, Phys.Rev.B 56, 5061 (1997). forM =10,q =4andD˜ =q+1. Foreachtemperature, [8] F. Verstraete, J. J. Garc´ia-Ripoll, and J. I. Cirac, Phys. three different regions can be distinguished: the super- Rev. Lett 93, 207204 (2004), cond-mat/0406426. fluid regionwith constantFWHM, the Mott-regionwith [9] M. Zwolak and G. Vidal, Phys. Rev. Lett 93, 207205 (2004), cond-mat/0406440. linearlyincreasingFWHMandanintermediateregionin [10] M. Fannes, B. Nachtergaele, and R. F. Werner, Comm. which both phases coexist. The value U/J at which the Math. Phys.144, 443 (1992). Mott-regionstartsincreaseswithincreasingtemperature, [11] F. Verstraete, M. A. Martin-Delgado, and J. I. Cirac, which is consistent with the criteria U ≫kBT,J for the Phys. Rev.Lett 92, 087201 (2004), quant-ph/0311087. appearance of the Mott-phase. This behaviour could be [12] T. St¨oferle, H. Moritz, C. Schori, M. K¨ohl, and 5 T. Esslinger, Phys. Rev. Lett. 92, 130403 (2004), cond- [20] M. P. Fisher, P. B. Weichmann, G. Grinstein, and D. S. mat/9312440. Fisher, Phys. Rev.B 40, 546 (1989). [13] B. Paredes, A. Widera, V.Murg, O. Mandel, S. F¨olling, [21] H. P. Bu¨chler, G. Blatter, and W. Zwerger, Phys. Rev. I.Cirac,G.V.Shlyapnikov,T.W.H¨ansch,andI.Bloch, Lett. 90, 130401 (2003). Nature429, 277 (2004). [22] C.Kollath,U.Schollw¨ock,J.vonDelft,andW.Zwerger, [14] M. K¨ohl, T. St¨oferle, H. Moritz, C. Schori, and Phys. Rev.A 69, 031601 (2004), cond-mat/0310388. T. Esslinger, Appl. Phys. B 79, 1009 (2004), cond- [23] S. Wessel, F. Alet, S. Trebst, D. Leumann, M. Troyer, mat/0411473. and G. G. Matrouni (2004), cond-mat/0411473. [15] M. Girardeau, J. Math. Phys. 1, 516 (1960). [24] L. Pollet, S. M. A. Rombouts, and P. J. H. Denteneer, [16] F. Verstraeteand J. I. Cirac (2004), cond-mat/0407066. Phys. Rev.Lett. 93, 210401 (2004), cond-mat/0408596. [17] F.Verstraete,D.Porras,andJ.I.Cirac,Phys.Rev.Lett [25] Notethatinpractice,itwillbedesirabletousethehigher 93, 227205 (2004), cond-mat/0404706. order versions of theTrotter decomposition. [18] F. Barahona, J. Phys. A 15, 3241 (1982). [26] Wenotethatwehavestoppedouriterativealgorithm at [19] D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, the point the variation of δK was less than 10−8. and P. Zoller, Phys. Rev. Lett 81, 3108 (1998), cond- mat/9805329.

See more

The list of books you might like

Most books are stored in the elastic cloud where traffic is expensive. For this reason, we have a limit on daily download.