Hybrid Monte-Carlo simulation of interacting tight-binding model of graphene
aa r X i v : . [ h e p - l a t ] N ov Hybrid Monte-Carlo simulation of interactingtight-binding model of graphene
Dominik Smith ∗ Theoriezentrum, Institut für Kernphysik, TU Darmstadt, 64289 Darmstadt, GermanyE-mail: [email protected]
Lorenz von Smekal
Theoriezentrum, Institut für Kernphysik, TU Darmstadt, 64289 Darmstadt, GermanyInstitut für Theoretische Physik, Justus-Liebig-Universität, 35392 Giessen, GermanyE-mail: [email protected]
In this work, results are presented of Hybrid-Monte-Carlo simulations of the tight-binding Hamil-tonian of graphene, coupled to an instantaneous long-range two-body potential which is modeledby a Hubbard-Stratonovich auxiliary field. We present an investigation of the spontaneous break-ing of the sublattice symmetry, which corresponds to a phase transition from a conducting toan insulating phase and which occurs when the effective fine-structure constant a of the systemcrosses above a certain threshold a C . Qualitative comparisons to earlier works on the subject(which used larger system sizes and higher statistics) are made and it is established that a C is of aplausible magnitude in our simulations. Also, we discuss differences between simulations usingcompact and non-compact variants of the Hubbard field and present a quantitative comparison ofdistinct discretization schemes of the Euclidean time-like dimension in the Fermion operator. ∗ Speaker. c (cid:13) Copyright owned by the author(s) under the terms of the Creative Commons Attribution-NonCommercial-ShareAlikeLicence. http://pos.sissa.it/
MC simulation of interacting tight-binding model of graphene
Dominik Smith
1. Introduction
In recent years, much interest has arisen in the properties of graphene, a one atom thick sheetof carbon atoms arranged on a hexagonal "honeycomb" lattice. It has become clear that such asimple structure generates a magnitude of unusual quantum effects, which not only make graphenea promising candidate for a wide range of technological applications but also a great model systemto study processes commonly associated with high-energy physics. The later stems from the factthat the tight-binding Hamiltonian which describes electrons in the p -orbitals of the carbon atoms(with additional terms describing electromagnetic two-body interactions) is well approximated by avariant of Quantum Electrodynamics in 2 + v F ≈ c /
300 of the electrons generates an effective fine structure constantwhich is a ≈ / ≈ . a , and whether the critical value lies in a range which can be experi-mentally realized (suspended graphene giving an upper bound). Since recent experiments provideevidence that graphene in vacuum is in fact a conductor [3], it is important to establish modelcalculations which match the observation.Early attempts at simulating graphene on the lattice studied the low-energy limit only, since theapplication of QCD methods is most direct here. Most prominently, Refs. [4] simulated the low-energy theory using staggered Fermions, whereas Ref. [5] investigated a variant of the Thirring model which has many similarities with QED + . Both find a phase transition to a gapped phasefor a > a C ≈
1, which is well within the physical region. More recently, a path integral formula-tion of the partition function was derived directly from the tight-binding theory, which preservesthe hexagonal structure of graphene and employs unphysical discretization in the time dimensiononly [6], with the spatial lattice spacing as a free parameter that can be taken from experiment.In this formulation, interactions are modeled by a non-local potential which is generated by aHubbard-Stratonovich field. It is immediately clear that this has many advantages. Not only does italleviate issues concerning parameter matching, it also opens up new opportunities to study physicsbeyond low energies, such as the effect of interactions on the neck-disrupting Lifshitz transition,which occurs in the tight-binding theory at the Van Hove singularity at finite density [7]. Recently,two distinct simulations of the hexagonal lattice have been conducted: One where an unscreened For an extensive but far from exhaustive review of graphene phenomenology see Refs. [1]. MC simulation of interacting tight-binding model of graphene
Dominik Smith
Coulomb potential is modeled by gauge links [8] and another where the instantaneous two-bodypotential was used [9]. The later work chose a form of the potential which accounts for additionalscreening by electrons in the s -orbitals, which was calculated within a constrained random phaseapproximation (cRPA) in Ref. [10]. While simulations with the unscreened potential yielded aprediction for a C consistent with the previous simulations of the low-energy theory, Ref. [9] foundevidence that screening is a mechanism which can move a C to larger values, possibly outside ofthe physically accessible region.In this work, first results of a very recently completed Hybrid-Monte-Carlo code are presented,which simulates the hexagonal theory with a non-local potential as in Refs. [6] entirely on GPUs.We present an investigation of the conductor-insulator phase transition on a rectangular graphenesheet of N x = N y = a C which is of a plausible magnitudecompared to the previous works. We also conduct an investigation of the discretization errors ofthe second order discretization scheme developed in Ref. [9] and discuss distinct treatments of theHubbard field. Our long-term goal is to obtain a precise prediction for a C and to investigate thephysics of the Van Hove singularity for the interacting case.
2. The path-integral and Hybrid-Monte-Carlo
A detailed derivation of the path integral formulation of the partition function of the interactingtight-binding model is presented in Refs. [6]. We only review the crucial steps here. The startingpoint is the Hamiltonian of the model in second quantized form, which is H = H tb + H c + H m = (cid:229) h x , y i , s ( − k )( a † x , s a y , s + a † y , s a x , s ) + (cid:229) x , y q x V xy q y + (cid:229) x m S ( a † x , + a x , + + a x , − a † x , − ) (2.1)Here a † x , s , a x , s are Fermionic creation- and annihilation operators which respect the usual anti-commutation relations. s denotes the electron spin and takes the values ± k ≈ . m S = ± m is a “staggered” mass, which has a different sign on each sub-latticeand which is added to serve as a seed for sub-lattice symmetry breaking (simulations are extrapo-lated to m → V xy need notbe further specified at this point, other than that it be positive definite. q x = a † x , a x , + a † x , − a x , − − b † x = a x , − , b x = a † x , − and flip the sign of these on one sub-lattice.The functional integral for Z = Tr e − b H (here b is 1 / k B T ) can now be derived using a co-herent state formalism. The essential step is to factor the exponential into N t terms e − b H = e − d H e − d H . . . e − d H , ( d = b / N t ) and to insert between them complete sets of Fermionic coherent The temperature is that of the gas of electronic quasi-particles only and is not equal to the physical temperature ofthe graphene sheet. MC simulation of interacting tight-binding model of graphene
Dominik Smith states h y k , h k | = h | e − (cid:229) x ( a x y ∗ x , k + b x h ∗ x , k ) , | y k , h k i = e − (cid:229) x ( y x , k a † x + h x , k b † x ) | i . This leads to the expressionTr e − b H = Z N t − (cid:213) t = (cid:20) (cid:213) x d y ∗ x , t d y x , t d h ∗ x , t d h x , t (cid:21) e − (cid:229) x ( y ∗ x , t + y x , t + + h ∗ x , t + h x , t + ) h y t + , h t + | e − d H | y t , h t i . (2.2)Here y x , t and h x , t are Grassman valued field variables which replace the ladder operators. Anti-periodic boundary conditions in the time-direction are implied ( y x , N t = − y x , , h x , N t = − h x , ) . Wenow use the identity h y , h | F ( a † l , a l , b † l , b l ) | y ′ , h ′ i = F ( y ∗ l , y ′ l , h ∗ l , h ′ l ) e (cid:229) l y ∗ l y ′ l + h ∗ l h ′ l , whichuses coherent states to map normal ordered functions of ladder operators to functions of Grassmanvariables and obtainTr e − b H = Z N t − (cid:213) t = (cid:20) (cid:213) x d y ∗ x , t d y x , t d h ∗ x , t d h x , t (cid:21) exp n − d h (cid:229) x , y Q x , t + , t V xy Q y , t + , t − (cid:229) h x , y i k ( y ∗ x , t + y y , t + y ∗ y , t + y x , t + h ∗ y , t + h x , t + h ∗ x , t + h y , t ) + (cid:229) x m S ( y ∗ x , t + y x , t + h ∗ x , t + h x , t )+ (cid:229) x V xx ( y ∗ x , t + y x , t + h ∗ x , t + h x , t ) i − (cid:229) x (cid:2) y ∗ x , t + ( y x , t + − y x , t ) + h ∗ x , t + ( h x , t + − h x , t ) (cid:3)o . (2.3)Here we have applied normal ordering to H , which leads to the additional term ∼ V xx and then con-sidered e − d H to be normal ordered, which implies a discretization error O ( b / N t ) (which vanishesfor N t → ¥ ). Also, we have defined a charge field Q x , t , t ′ = y ∗ x , t y x , t ′ − h ∗ x , t h x , t ′ .Eq. (2.3) contains fourth powers of y , h . We must get rid of these if we wish to carry out Gaus-sian integration and obtain a Fermion determinant, which is necessary for HMC. This is achievedby applying the Hubbard-Stratonovich transformation e − d (cid:229) Nt − t = (cid:229) x , y Q x , t + , t V xy Q y , t + , t ∼ Z " N t − (cid:213) t = (cid:213) x d f x , t e − d (cid:229) Nt − t = (cid:229) x , y f x , t V − xy f y , t − i d (cid:229) Nt − t = (cid:229) x f x , t Q x , t + , t , (2.4)which eliminates fourth powers, at the expense of introducing an auxiliary field f x , t and introducingthe inverse of the matrix V xy . The constant factor before the integral in Eq. (2.4) can be dropped.Combining equations (2.3) and (2.4) and carrying out Gaussian integration finally yieldsTr e − b H = Z " N t − (cid:213) t = (cid:213) x d f x , t e − d (cid:229) Nt − t = (cid:229) x , y f x , t V − xy f y , t | det M ( f ) | . (2.5)This form is suitable for simulation via HMC. The determinant can be sampled stochastically byusing pseudo-Fermion sources and thus f remains as the only "dynamical" field. The matrix M isdefined in terms of its components as N t b M ( x , t )( y , t ′ ) = N t b d xy ( d tt ′ − d t − , t ′ ) − k (cid:229) ~ n d y , x + ~ n d t − , t ′ + m S d xy d t − , t ′ + V xx d xy d t − , t ′ + i f x , t d xy d t − , t ′ . (2.6)Here ~ n denotes vectors connecting nearest-neighbor sites. Since M is hermitian, there is no sign4 MC simulation of interacting tight-binding model of graphene
Dominik Smith problem.Simulations based on Eq. (2.6) have a severe problem: The term ∼ f makes det M ( f ) apolynomial of order f N where N is the number of lattice sites. Thus, rounding errors are amplifiedin an uncontrollable way. As is discussed further in Section 3, this indeed make simulations using(2.6) impossible. The solution was worked out in Refs. [6], where it is shown that one may replace b N t V xx d xy d t − , t ′ + i b N t f x , t d xy d t − , t ′ −→ e i b Nt f x , t d xy d t − , t ′ , (2.7)which transforms f into a phase that is additive in det M ( f ) and thus numerically stable.It should be noted that that (2.6) is by far not the only possible form of M . As Refs. [6]discuss, there is a great deal of freedom in discretizing the time-direction and this can be exploitedto construct improved actions which converge faster to the continuum limit. A second order actionwas presented in Ref. [9], which is constructed by factoring exp ( − d H ) such that the interactingpart is split off, and inserting an additional set of coherent states between the terms:Tr e − b H = Z " N t − (cid:213) t = (cid:213) x d y ∗ x , t d y x , t d h ∗ x , t d h x , t N t − (cid:213) t = e − (cid:229) x ( y ∗ x , t y x , t + h ∗ x , t h x , t + y ∗ x , t + y x , t + + h ∗ x , t + h x , t + ) × h y t , h t | e − d ( H tb + H m ) | y t + , h t + ih y t + , h t + | e − d H C | y t + , h t + i . (2.8)Carrying out the calculation in analogy to what was discussed above (using the compact version ofthe Hubbard field), one obtains the following Fermion operator: M ( x , t )( y , t ′ ) = d xy ( d tt ′ − d t + , t ′ ) − b N t k (cid:229) ~ n d y , x + ~ n d t + , t ′ + b N t m S d xy d t + , t ′ : t even d xy d tt ′ − d xy d t + , t ′ exp ( − i b N t f x , ( t − ) / ) : t odd (2.9)Note that the Hubbard field now only appears on odd time-slices. It was hypothesized in Ref. [9]that this version of the Fermion operator has reduced discretization errors.We wish to investigate spontaneous breaking of sub-lattice symmetry. Thus we require aproper order parameter. The obvious choice is to use the difference of number density operators onthe two sub-lattices (denoted here as A and B ). In the functional integral form, this is expressed as h D N i = ZN t Z D y D y ∗ D h D h ∗ h (cid:229) X A , t (cid:0) y ∗ x , t + y x , t + h ∗ x , t + h x , t (cid:1) − (cid:229) X B , t (cid:0) y ∗ x , t + y x , t + h ∗ x , t + h x , t (cid:1) i e − b H = − b Z Z D f (cid:20) ¶¶ m det (cid:0) MM † (cid:1)(cid:21) e − S [ f ] = − N t N t − (cid:229) t = h (cid:229) x ∈ A M − ( x , t + )( x , t ) − (cid:229) x ∈ B M − ( x , t + )( x , t ) i . (2.10)The above is valid for the first order discretization scheme above. For the second order Fermionoperator, one may insert b D N on even time-slices only to obtain an analogous expression.
3. Results
We have simulated rectangular N x = N y = MC simulation of interacting tight-binding model of graphene
Dominik Smith operators. We chose a potential which is constructed piecewise: The on-site ( V ), nearest-neighbor( V ), next-nearest neighbor ( V ) and crossing term (across the hexagon, V ) are the cRPA valuesof Ref. [10], while the long-range part is an unscreened Coulomb potential. We account for pe-riodicity by adding one set of eight mirror images, arranged along the rectangular boundary. Thismay be problematic due to the zero mode in the potential. Moreover, we have chosen V = . V = . b = . V − (asin Ref. [9]). An effective a = ( / )( / e ) is introduced as a free parameter which impliesa rescaling of V xy → V xy / e . Our first observation is that the non-compact Hubbard field is indeedproblematic. Foremost, we were not able to observe distinctly non-zero expectation values of theorder parameter. We thus use the compact version only. The table below shows a comparison ofthe discretization errors of the first and second order discretization schemes (labeled "standard"and "improved" in the following). These were obtained by measuring h D N i on roughly one hun-dred independent configurations for each combination of ( a , m , N t ) and fitting results obtained fordifferent N t ∈ [ , . . . , ] to h D N i = c ∗ ( / N t ) + c . b [ e V − ] a m [ e V ] c , std c , imp c , std c , imp . .
87 0 . − . ( ) − . ( ) . ( ) . ( ) . .
87 0 . − . ( ) − . ( ) . ( ) . ( ) . .
55 0 . − . ( ) − . ( ) . ( ) . ( ) . .
18 0 . − . ( ) − . ( ) . ( ) . ( ) These results provide strong evidence that both versions not only converge to the same continuumlimit but are in fact equal for any N t . It thus appears that nothing is to be gained by using the secondorder action.Subsequently we performed an investigation of the phase transition. We simulated for manydifferent choices of a , where masses were chosen as m = . , . , . , . , . N t = [ , . . . , ] for each set of ( a , m ) . Againroughly one hundred independent measurements were done for each set ( a , m , N t ) The continuumresults are shown in Fig. 1. We have extrapolated to m = a using h D N i = c m + c m + c . We find that our choice of potential generates an a C ≈
4. Conclusions and Outlook
We have demonstrated here that our code produces qualitatively the expected behavior, witha phase transition occurring at an a C ≈ a C ≈ . MC simulation of interacting tight-binding model of graphene
Dominik Smith < D n > a m = 0.1m = 0.2m = 0.3m = 0.4m = 0.5 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 < D n > m a = 5.45 a = 5.12 a = 4.83 a = 3.61 a = 2.55 Figure 1:
Phase transition of N x = N y = m → a . which now allows us to simulate at larger volumes. When choosing the same potential as Ref. [9]and properly addressing boundary conditions, it appears we are now able to reproduce their resultsalso quantitatively. We are thus confident that our code is now ready for large-scale operation. Acknowledgements
We have benefited from discussions with Pavel Buividovich, Maxim Ulybyshev and the lateMikhail Polikarpov. This work was supported by the Deutsche Forschungsgemeinschaft withinSFB 634, by the Helmholtz International Center for FAIR within the LOEWE initiative of the Stateof Hesse, and the European Commission, FP-7-PEOPLE-2009-RG, No. 249203. All results wereobtained using Nvidia GTX and Tesla graphics cards.
References [1] A. H. Castro Neto et al., Rev. Mod. Phys. , 109 (2009); V. N. Kotov et al., Rev. Mod. Phys. ,1067 (2012) [arXiv:1012.3484].[2] V. P. Gusynin et al., Int. J. Mod. Phys. B (2007) 4611 [arXiv:0706.3016].[3] D. C. Elias et al., Nature Phys. , 701 (2011) [arXiv:1104.1396]; A. S. Mayorov et al., Nano Lett. ,4629 (2012), [arXiv:1206.3848].[4] J. E. Drut and T. A. Lähde, Phys. Rev. Lett. , 026802 (2009) [arXiv:0807.0834]; Phys. Rev. B ,165425 (2009) [arXiv:0901.0584].[5] W. Armour, S. Hands and C. Strouthos, Phys. Rev. B , 125105 (2010) [arXiv:0910.5646].[6] R. Brower, C. Rebbi and D. Schaich, PoS LATTICE 2011, 056 (2012) [arXiv:1204.5424];[arXiv:1101.5131].[7] B. Dietz, L. von Smekal et al., [arXiv:1304.4764].[8] P. V. Buividovich and M. I. Polikarpov, Phys. Rev. B 86, 245117 (2012) [arXiv:1206.0619].[9] M. V. Ulybyshev, P. V. Buividovich, M. I. Katsnelson and M. I. Polikarpov, Phys. Rev. Lett. 111,056801 (2013) [arXiv:1304.3660].[10] T. O. Wehling et. al., Phys. Rev. Lett. 106, 236805 (2011) [arXiv:1101.4007]., 125105 (2010) [arXiv:0910.5646].[6] R. Brower, C. Rebbi and D. Schaich, PoS LATTICE 2011, 056 (2012) [arXiv:1204.5424];[arXiv:1101.5131].[7] B. Dietz, L. von Smekal et al., [arXiv:1304.4764].[8] P. V. Buividovich and M. I. Polikarpov, Phys. Rev. B 86, 245117 (2012) [arXiv:1206.0619].[9] M. V. Ulybyshev, P. V. Buividovich, M. I. Katsnelson and M. I. Polikarpov, Phys. Rev. Lett. 111,056801 (2013) [arXiv:1304.3660].[10] T. O. Wehling et. al., Phys. Rev. Lett. 106, 236805 (2011) [arXiv:1101.4007].