We present a novel method for precise numerical solution of the irreducible
two-body problem and apply it to excitons in solids. The approach is based on
the Monte Carlo simulation of the two-body Green function specified by
Feynman's diagrammatic expansion. Our method does not rely on the specific form
of the electron and hole dispersion laws and is valid for any attractive
electron-hole potential. We establish limits of validity of the Wannier (large
radius) and Frenkel (small radius) approximations, present accurate data for
the intermediate radius excitons, and give evidence for the charge transfer
nature of the monopolar exciton in mixed valence materials.Comment: 4 pages, 5 figure

A. S. Mishchenko

B. V. Svistunov

E. A. Burovski

G. Travaglini

I. Egri

J. H. Wannier

J. I. Frenkel

K. A. Kikoin

L. X. Benedict

M. H. Cohen

M. Jarrell

M. Rohlfing

N. V. Prokof'ev

P. Lemmens

R. Knox

S. Albrecht

S. Curnoe

T. Kasuya

English

arXiv

This is the pre-published version harvested from ArXiv. The published version is located at http://prl.aps.org/abstract/PRL/v87/i18/e186402We present a novel method for the precise numeric solution of the irreducible two-body problem and apply it to excitons in solids. The approach is based on the Monte Carlo simulation of the two-body Green function specified by Feynman’s diagrammatic expansion. Our method does not rely on the specific form of the electron and hole dispersion laws and is valid for any attractive electron-hole potential. We establish limits of validity of the Wannier (large radius) and Frenkel (small radius) approximations, present accurate data for the intermediate radius excitons, and give evidence for the charge transfer nature of the monopolar exciton in mixed valence materials

Burovski, E

Mishchenko, A

Prokof'ev, Nikolai

Svistunov, B

ScholarWorks@UMass Amherst

Diagrammatic Quantum Monte Carlo for Two-Body Problems: Applied to Excitons

Crossref

We present a novel method for the precise numeric solution of the irreducible two-body problem and apply it to excitons in solids. The approach is based on the Monte Carlo simulation of the two-body Green function specified by Feynman’s diagrammatic expansion. Our method does not rely on the specific form of the electron and hole dispersion laws and is valid for any attractive electron-hole potential. We establish limits of validity of the Wannier (large radius) and Frenkel (small radius) approximations, present accurate data for the intermediate radius excitons, and give evidence for the charge transfer nature of the monopolar exciton in mixed valence materials

Prokof\u27ev, Nikolai

University of Massachusetts AmherstScholarWorks@UMass AmherstPhysics Department Faculty Publication Series Physics2001Diagrammatic Quantum Monte Carlo for Two-Body Problems: Applied to ExcitonsE BurovskiA MishchenkoNikolai Prokof 'evUniversity of Massachusetts - Amherst, prokofev@physics.umass.eduB SvistunovFollow this and additional works at: https://scholarworks.umass.edu/physics_faculty_pubsPart of the Physical Sciences and Mathematics CommonsThis Article is brought to you for free and open access by the Physics at ScholarWorks@UMass Amherst. It has been accepted for inclusion in PhysicsDepartment Faculty Publication Series by an authorized administrator of ScholarWorks@UMass Amherst. For more information, please contactscholarworks@library.umass.edu.Recommended CitationBurovski, E; Mishchenko, A; Prokof 'ev, Nikolai; and Svistunov, B, "Diagrammatic Quantum Monte Carlo for Two-Body Problems:Applied to Excitons" (2001). Physics Review Letters. 1190.Retrieved from https://scholarworks.umass.edu/physics_faculty_pubs/1190arXiv:cond-mat/0106207v1  [cond-mat.str-el]  12 Jun 2001Diagrammatic Quantum Monte Carlo for Two-Body Problem: ExcitonE. A. Burovski1, A. S. Mishchenko2,1, N. V. Prokof’ev3, and B. V. Svistunov11RRC ’Kurchatov Institute’, 123182, Moscow, Russia2Correlated Electron Research Center, Tsukuba Central 4, Tsukuba 305-8562, Japan3Department of Physics, University of Massachusets, Amherst, Masachusets 01003We present a novel method for precise numerical solution of the irreducible two-body problem andapply it to excitons in solids. The approach is based on the Monte Carlo simulation of the two-bodyGreen function specified by Feynman’s diagrammatic expansion. Our method does not rely on thespecific form of the electron and hole dispersion laws and is valid for any attractive electron-holepotential. We establish limits of validity of the Wannier (large radius) and Frenkel (small radius)approximations, present accurate data for the intermediate radius excitons, and give evidence forthe charge transfer nature of the monopolar exciton in mixed valence materials.PACS numbers: 71.53.-y, 02.70.Ss, 05.10.LnAfter it was realized that under certain conditions theelectron dynamics in conduction band is of two-particlenature due to Coulomb attaction to the hole in the va-lence band left behind [1], the problem of exciton becamea model example of an irreducible (center-of-mass motiondoes not separate from the rest of degrees of freedom)two-body problem. The simplest (still rather general)exciton Hamiltonian [2,3] consists of conduction and va-lence band contributions, H0, and coupling He-h:H0 =∑kεc(k)e†kek +∑kεv(k)hkh†k, (1)He-h = −N−1∑pkk′U(p,k,k′)e†p+kh†p−khp−k′ep+k′. (2)Here ek (hk) is the electron (hole) annihilation operator,εc(k) (εv(k)) is the conduction (valence) band dispersionlaw, N is the number of lattice sites, and U(p,k,k′) is anattractive interaction potential.Despite numerous efforts over the years there is no rig-orous technique to solve for exciton properties even forthe simplest model given above which treats electron-electron interactions as a static renormalized Coulombpotential with averaged dynamical screening. The onlysolvable cases are the Frenkel small-radius limit [1] andthe Wannier large-radius limit [4] which describe molecu-lar crystals and wide gap insulators with large dielectricconstant, respectively. Much more frequently encoun-tered cases of intermediate radius excitons (e.g. interme-diate gap semiconductors, LiF, or mixed valence systems)have to be dealt with using approximate numerical ap-proaches. There are powerful ab initio modern methods[5–7] for band structure and effective electron-hole poten-tial calculations, but the real bottleneck is in numericalsolution of the two-particle problem for a bulk material.One can either solve the Bethe-Salpeter equation on afinite mesh in reciprocal/direct space [5–7], or employthe random-phase approximation decoupling [2]. How-ever, both methods suffer from systematic errors, andthe Bethe-Salpeter equation on finite mesh may lead toincorrect eigenstates for the Wannier case [7]. Therefore,even the limits of validity of the Wannier and Frenkelapproximations can not be established by existing meth-ods.Besides, an efficient and rigorous method for the studyof exciton properties, given the band structure, is of highvirtue for phenomenological models. As an example, werefer to the protracted discussion of numerous (and of-ten contradictory) models concerning exciton propertiesin mixed valence semiconductors [8]. In Ref. [9] un-usual properties of SmS and SmB6 were explained byinvoking the excitonic instability mechanism assumingcharge-transfer nature of the optically forbidden exciton.Although this model explains quantitatively the phononspectra [10], optical properties [11,12], and magnetic neu-tron scattering data [13], its basic assumption has beencriticized as being groundless [14].In this Letter we describe how ground state proper-ties of excitons in the model (1)-(2) can be obtained nu-merically without systematic errors for arbitrary disper-sion relations εc(k) and εv(k)), and attractive potentialU(p,k,k′). First, we show that the problem fits intothe diagrammatic Monte Carlo (MC) method [15–17]which sums positively-definite perturbation series, in ourcase Feynman diagrams, for the two-particle MatsubaraGreen function, G. We then describe the procedure of ex-tracting various physical properties from the asymptoticlong-time behavior of G. Next, we discuss our resultsfor a particular tight-binding model and electron-hole in-teraction potential to see under what conditions Frenkeland Wannier approximations remain accurate. Finally,we present evidence that the band structure of mixed va-lence materials results in the charge-transfer character ofthe optically forbidden exciton.The two-particle Green function with total momentum2p in imaginary time representation is defined asGkk′p (τ) = 〈0 | ep+k′(τ)hp−k′(τ)h†p−ke†p+k | 0〉, (3)where the vacuum state | 0〉 corresponds to empty1conduction and filled valence bands, and hp−k(τ) =eHτhp−ke−Hτ , τ > 0. In the interaction representationG can be written as a sum of ladder-type Feynman di-agrams, see Fig. 1: pairs of horizontal solid lines repre-sent free electron-hole pair propagators,G(0)p (k, τ2−τ1) =exp (−ε(k)(τ2 − τ1)), where ε(k) = εc(p+ k)−εv(p− k)is the energy of the pair, and dashed lines representthe interaction potential. For purely numerical rea-0 τ1 τ2 τ3 τp+ kp− kp+ k1p− k1p+ k2p− k2p+ k′p− k′V1(2p) V2(k2 − k1) V2(k′ − k2)FIG. 1. A typical diagram contributing to Gkk′p (τ ).sons explained below we split the potential into twoterms U(p,k,k′) = V1(2p) + V2(k − k1), and expandin both V1 and V2. [This can be done because [2]U(p,k,k′) = V0 −W (2p) + U(k− k′) where V0 is theon-site coupling,W (2p) is the dipolar term, U(k− k′) =∑λ6=0 exp(iqRλ)/(Rλǫ(Rλ)) is the monopolar term, andǫ(Rλ) is a static dielectric screening function (we set theelectric charge to unity). Since U(q) is not positive def-inite [in fact∑q U(q) = 0] we add and subtract someconstant U¯ to ensure that V1(2p) = V0−W (2p)− U¯ andV2(q) = U¯ + U(q) are both positive - this imposes theonly limitation on value V0 in our method].The final answer for G is given by the sum of all pos-sible diagrams. Formally we can write this as a series ofmulti-dimensional integralsGkk′p (τ) =∞∑m=0∑ξm∫dx1 · · · dxm Fkk′p (τ ; ξm;x1, . . . , xm).where x1, . . . xm are internal variables [times and mo-menta, xi = (τi,ki)] of the m-th order diagram, thesummation over ξm accounts for different diagrams oforder m, and the “weight” F is given by the productof electron-hole propagators and interaction vertices ac-cording to standard rules. For positive V1 and V2 allterms in the series are positive definite and one may ap-ply the diagrammatic Monte Carlo technique developedin Refs. [15,16], which evaluates such series without sys-tematic errors (by Metropolis-type sampling of diagramsaccording to their weight directly in the momentum-timecontinuum). Since the method itself is well described inthe literature we will concentrate on the problem specificdetails only.The crucial for the whole scheme update is theone which changes the number of interaction ver-tices by one. To render algorithm efficient onehas to propose new internal parameters as closeas possible to the distribution function R(xm+1) =F (ξm+1;x1, . . . , xm, xm+1)/F (ξm;x1, . . . , xm) defined bythe ratio of the new and old diagram weights. This isdone in order to maximize the acceptance ratio Pacc — ifproposed xm+1 = (τm+1,km+1) are distributed accord-ing to some normalized function W (xm+1), then Pacc ∝R(xm+1)/W (xm+1). Otherwise the choice of W (xm+1)is a matter of computational convenience [15,16].First, we select (with equal probabilities) which inter-action vertex, V1 or V2, will be inserted, and then select atrandom the time interval where it will be placed, τm+1 ∈(τa, τb), where τa,b are the interval boundaries determinedeither by the existing interaction vertices or the diagramends. All the momenta at τ < τa and τ > τm+1 are keptuntouched. In case when V1(2p) is inserted, the newmomentum km+1 is proposed uniformly in the Brillouinzone (BZ). When V2(kb−km+1) is inserted (kb is the rel-ative motion momentum to the left of point (τb) the newmomentum km+1 is proposed using distribution functionW (km+1) = (β/2π arctanβ)3∏α(1 + βk(α)m+1/π)2)−1,where α = x, y, z. The parameter β is uniformly seededon interval [βmin, βmax] at each step, and βmin, βmax arefurther tuned to maximize the acceptance ratio. Wenote, that different distribution functions used to pro-pose new momentum km+1 when dealing with V1 and V2vertices was the only motivation behind an artificial sep-aration U = V1 + V2 [The actual gain in efficiency wasabout three orders of magnitude!]. Finally, the time po-sition for the new vertex was seeded according to thedistribution function dictated by the diagram weightsratio W (τm+1) = δǫ · e−δǫτm+1/(e−δǫτa − e−δǫτb) where,δǫ = ε(km+1)− ε(kb).We also employ standard Metropolis updates changingthe values of internal momenta and times, which substan-tially enhances the efficiency of the algorithm.We now turn to the discussion of how exciton proper-ties are obtained from the G(τ → ∞) limit. An eigen-state | ν;p〉 with energy Eν can be written as| ν;p〉 ≡∑kξpk(ν)e†p+kh†p−k | 0〉. (4)where amplitudes ξpk(ν) = 〈ν;p | e†p+kh†p−k | 0〉 de-scribe the wave function of internal motion of the exci-ton. In terms of exciton eigenstates we have, Gk=k′p (τ) =∑ν | ξpk(ν) |2 e−Eντ , and if τ is much larger than inverseenergy difference between the ground and first excitedstates, the Green function projects to the ground state,Gk=k′p (τ → ∞) =| ξpk(g.s.) |2 e−Eg.s.τ . Due to nor-malization condition∑k | ξpk(ν) |2≡ 1 the asymptoticbehavior of the sum G˜p =∑kGk=k′p is especially sim-ple, G˜(τ) → e−Eg.s.τ . This asymptotic behavior allowssimulations of energy and amplitudes at fixed τ [largeenough to make the corresponding systematic error neg-ligible], using the technique of Monte Carlo estimators.20 20 40 60 80One−particle bandwidth−4−3−2−10log ( Binding energy )0 0.5 1One−particle bandwidth0.20.40.60.81One−particle bandwidthFIG. 2. The dependence of the exciton binding energy onthe bandwidth Ec = Ev. Statistical errors are less than5 · 10−3 in relative units. The dashed line corresponds tothe Wannier model. The solid line is the cubic spline, thederivatives at the right and left ends being fixed by the Wan-nier limit and perturbation theory, respectively. Insert: theinitial part of the plot.To this end we differentiate each diagram for G˜(τ) [18]with respect to τ and arrive at the result (compare withRef. [16])Eg.s. = τ−1〈m+1∑j=1εj(k)∆τj − m〉MC, (5)where 〈...〉MC stands for the MC statistical average, mis the diagram order, εj(k) and ∆τj are the electron-hole pair energy and duration of the j-th propagator,respectively. By definition, in the limit τ → ∞ wehave Gk=k′p /G˜p =| ξpk(g.s.) |2, i.e. the distribution overquasimomentum k is related to the wave function of in-ternal motion. The wave function of the bound state canbe chosen real, and the Fourier transform may be usedto obtain | g.s.〉 in direct space [19].In this Letter we focus on the study of exciton prop-erties in a simple cubic 3D lattice with tight binding dis-persion laws for the electron and hole bandsεc,v(k) = E˜c,v ± (Ec,v/6)∑α(1− cos kα). (6)The choice of interaction parameters was motivated bythe possibility to cover all regimes (from Wannier toFrenkel limit) by varying the ratio between the band-width and the gap only. Our simulations were done forE˜v = 0, Eg ≡ E˜c = 1, W (2p = 0) = −0.168, V0 = 0.778,U¯ = 0.578, and ǫ(R) = 10 [20]. The binding energy inthe Frenkel limit EFL (Ec,v ≪ Eg) is then less than thegap, EFL = V1(2p = 0) +∑q V2(q) = 0.946, thus ren-dering the exciton stability for all values of Ec,v. In the0 0.005 0.01 0.0150 0.1 0.2 0.3MomentumCharge density(arb.units)(a)(b)FIG. 3. The momentum dependence of the charge density| ξpk(g.s.) |2 k2 for Ec = Ev = 60 (a) and Ec = Ev = 10 (b).Solid lines are the Wannier model result. Statistical errorsare typically of order 10−4.Wannier limit of large bandwidth Ec,v ≫ Eg the bindingenergy approaches 3/(2ǫ2Ec) (assuming Ev = Ec). Ofcourse, our parameters satisfy the requirement that V1and V2 are positive definite functions.Our results for the binding energy and wave functionare shown in Figs. 2, 3, and 4. First we notice that themethod works equally well in all regimes, and statisticalerrors are much smaller than symbols sizes in all plots.An unexpected result is that extremely large bandwidthEc/Eg > 20 is necessary for the Wannier approximationto be adequate: both the binding energy EB (Fig. 2) andthe wave function [21] (Fig. 3 and Fig. 4 (b)) demonstratelarge deviations for smaller Ec/Eg. Most surprisingly,for 1 < Ec/Eg < 10 the wave function has a large (anddominating) on-site component [Fig. 4(b)], but the bind-ing energy is not even close to the Frenkel limit! ForEc/Eg = 0.4 the wave function is almost entirely local-ized [Fig. 4(c)] but Eg.s. is still 50% away from the small-radius limit. Noticing that Eg.s. ≈ EFL − (Ec + Ev)/2(Ec = Ev), which holds for localised functions whenEc < 0.4, we deduce that the deviation from Frenkel re-sult is determined by the electron and hole delocalizationenergy. Our conclusion is then that the intermediate-range regime is very broad and relevant in most practicalcases.To study the structure of optically forbidden excitonsin mixed valence compounds we choose typical for thesesemiconducting materials band spectra [8], i.e. an almostflat valence band separated by an indirect gap from thewide conduction band with maximum at k = 0 and mini-mum at the BZ boundary. One can see in Fig. 5 that thisleads to the charge transfer character of the optically for-bidden monopolar exciton (W (2p) = 0) when the wavefunction of internal motion has almost zero on-site com-ponent, maximal charge density at near neighbours, andlarge long-ranged oscillations at neighboring sites. The30 200 400 600 800 1000Relative e−h distance0 2 4 6 8 10Coordination sphere0.050.10.150.20.25Envelope function (arb. units)0 1 2 3 4Coordination sphere(a)(b)(c)FIG. 4. The wave function of internal motion in realspace: (a) Wannier [Ec = Ev = 60]; (b) intermediate[Ec = Ev = 10]; (c) near-Frenkel [Ec = Ev = 0.4] regimes.The solid line in the panel (a) is the Wannier model resultwhile solid lines in other panels are to guide an eye only. Sta-tistical errorbars are of order 10−4.0 5 10 15Coordination sphere number−0.1−0.0500.05Envelope functionFIG. 5. The wave function of internal motion in real spacefor the optically forbidden monopolar (W (2p) = 0) excitondefined by the following model parameters: E˜c = 1.5, E˜v = 0,Ec = −0.5, Ev = 0.05, ǫ = 10, V0 = 0.578. Statistical error-bars are of order 10−4.difference with the previously discussed Ev,c/Eg = 0.4case, see Fig. 4(b), is remarkable.Finally, we would like to note that diagrammatic MCtechnique not only gives properties of the ground statebut is also suitable for the study of excited states andoptical absorption [22]. This can be done by simulatingthe τ dependence of G(p = 0, τ) =∑kk′ Gkk′p=0(τ), andsolving numerically equationG(p = 0, τ) =∫ ∞0g(ω) exp(−ωτ)dωto obtain the spectral function g(ω) [16].We thank A. Sakamoto and N. Nagaosa for fruitfuldiscussions. We acknowledge the support of the Na-tional Science Foundation under Grant DMR-0071767and RFBR grant 01-02-16508.[1] J. I. Frenkel, Phys. Rev. 17, 17 (1931).[2] I. Egri, Phys. Reports 119, 364 (1985).[3] R. Knox, Theory of excitons, (Academic press, New York1963).[4] J. H. Wannier, Phys. Rev. 52, 191 (1937).[5] S. Albrecht, L. Reining, R. D. Sole, and G. Onida, Phys.Rev. Lett. 80, 4510 (1998).[6] L. X. Benedict, E. L. Shirley, and R. B. Bohn, Phys. Rev.Lett. 80, 4514 (1998).[7] M. Rohlfing and S. G. Louie, Phys. Rev. Lett. 81, 2312(1998).[8] S. Curnoe and K. A. Kikoin, Phys. Rev. B 61 15714(2000);[9] K. A. Kikoin and A. S. Mishchenko, J. Phys.: Condens.Matter, 2, 6491 (1990).[10] A. S. Mishchenko and K. A. Kikoin, J. Phys.: Condens.Matter, 3, 5937 (1991).[11] G. Travaglini and P. Wachter, Phys. Rev. B 29, 893(1984).[12] P. Lemmens, A. Hoffmann, A. S. Mishchenko, M. Yu.Talantov, and G. Gu¨ntherodt, Physica B 206&207, 371(1995).[13] K. A. Kikoin and A. S. Mishchenko, J. Phys.: Condens.Matter, 7, 307 (1995);[14] T. Kasuya, Europhys. Lett., 26, 277 (1994); T. Kasuya,Europhys. Lett., 26, 283 (1994)[15] N. V. Prokof’ev and B. V. Svistunov, Phys. Rev. Lett.81, 2514 (1998).[16] A. S. Mishchenko, N. V. Prokof’ev, A. Sakamoto, and B.V. Svistunov, Phys. Rev. B 62, 6317 (2000).[17] A. S. Mishchenko and N. Nagaosa, Phys. Rev. Lett. 86,4624 (2001).[18] Working with the function G˜ introduces a certain for-mal problem: the zero- and first-order (with respect toV2(q = 0)) diagrams contain macroscopically large factorN . Being interested in the ground-state properties only,we can safely omit them, since explicit analysis showsthat they are irrelevant in the τ →∞ limit.[19] If the quantity | ξpk(g.s.) |2 has no nodes, the amplitudesξpk(g.s.) can be chosen positive. Since Gkk′p (τ ) is givenby the sum of positive-definite diagrams we know thatξpk(g.s.) has no nodes.[20] The direct coupling term was evaluated by generalizedEvald method, see e.g. M. H. Cohen and F. Keffer, Phys.Rev. 99, 1128 (1955).[21] While for Wannier regime the wave function is sphericallysymmetric, it’s symmetry for the Frenkel and intermedi-ate case regimes bears the lattice structure thus makingthe coordination spheres number a natural unit for therepresentation of the wave function in direct space.[22] In the case when the interband optical matrix elementis momentum independent, the absorption coefficient isproportional to g(ω).4