Decomposition and embedding in the stochastic GW self-energy
DDecomposition and embedding in the stochastic GW self-energy Mariya Romanova and Vojtˇech Vlˇcek a) Department of Chemistry and Biochemistry, University of California, Santa Barbara, CA 93106-9510,U.S.A.
We present two new developments for computing excited state energies within the GW approximation. First,calculations of the Green’s function and the screened Coulomb interaction are decomposed into two parts:one is deterministic while the other relies on stochastic sampling. Second, this separation allows constructinga subspace self-energy, which contains dynamic correlation from only a particular (spatial or energetic) regionof interest. The methodology is exemplified on large-scale simulations of nitrogen-vacancy states in a periodichBN monolayer and hBN-graphene heterostructure. We demonstrate that the deterministic embedding ofstrongly localized states significantly reduces statistical errors, and the computational cost decreases bymore than an order of magnitude. The computed subspace self-energy unveils how interfacial couplingsaffect electronic correlations and identifies contributions to excited-state lifetimes. While the embedding isnecessary for the proper treatment of impurity states, the decomposition yields new physical insight intoquantum phenomena in heterogeneous systems.PACS numbers: 71.15.-m, 31.15.Md, 71.45.Gm, 71.10.-w, 71.15.Dx, 71.20.-b, 71.23.An, 71.55.-i, 73.20.-r,73.20.At, 73.21.-b, 73.43.Cd, 61.72.Ji I. INTRODUCTION
First-principles treatment of electron excitation ener-gies is crucial for guiding the development of new mate-rials with tailored optoelectronic properties. Localizedquantum states became the focal point for condensedsystems: due to their long-range electronic correlation,the localized excitations exhibit tunability by interfacialphenomena.
Besides predicting experimental observ-ables, theoretical investigations thus help to elucidate theinterplay of the electron-electron interactions and the roleof the environment.The excited electrons and holes near the Fermi levelare conveniently described by the quasiparticle (QP) pic-ture: the charge carriers are characterized by renormal-ized interactions and a finite lifetime limited by energydissipation, which governs the deexcitation mechanism.Quantitative predictions of QPs necessitate the inclusionof non-local many-body interactions.The prevalent route to describe QPs in condensed sys-tems employs the Green’s function formalism.
Theexcitation energy and its lifetime are inferred from theQP dynamics. The many-body effects are representedby the non-local and dynamical self-energy, Σ. In prac-tice, Σ is approximated by selected classes of interac-tions, forming a hierarchy of systematically improvablemethods. The formalism also allows constructing theself-energy for distinct states with a different form of Σ.Here we neglect the vibrational effects, and the induceddensity fluctuations dominate the response of the systemto the excitation. The perturbation expansion is conve-niently based on the screened Coulomb interaction, W ,which explicitly incorporates the system’s polarizability.The neglect of the induced quantum fluctuations (while a) Electronic mail: [email protected] accounting for the induced density) leads to the popu-lar GW method. This approach predicts the QPgap, ionization potentials, and electron affinities in excel-lent agreement with experiments.
Within the GW approximation, Σ is a product of W and the Green’sfunction, G , in the time domain.Conventional implementations of the GW self-energyscale as O ( N ) with the system size, albeit with a smallprefactor, and it is usually limited to few-electron sys-tems. Stochastic treatment recasts the self-energy asa statistical estimator and employs sampling of elec-tronic wavefunctions combined with the decompositionof operators. In practice, this formulation decreasesthe computational time significantly and leads to a linearscaling algorithm.
The stochastic approach uses real-space random vec-tors that sample the Hilbert space of the electronicHamiltonian. The statistical error in Σ is governed bythe number of vectors, which are used for the decom-position of G and the evaluation of W . While extendedsystems with delocalized states are treated efficiently, thestatistical fluctuations are more significant for localizedorbitals. This translates to an increased computationalcost for defects, impurities, and molecules.
Further,too high fluctuations potentially hinder proper conver-gence and may result in sampling bias.Here we overcome this difficulty by constructing a hy-brid deterministic-stochastic approach. We show how toefficiently decompose G and W in real-time and embedthe strongly localized states. In this formalism, the prob-lematic orbitals are treated explicitly without relying ontheir sampling by random vectors. A similar embeddingscheme was so far employed only in static ground-statecalculations. The decomposition is general and can beused to sample arbitrary excitations. Naturally, it is es-pecially well suited for spatially or energetically isolatedelectronic states. Further, we employ the decompositiontechnique to compute Σ stochastically from an arbitrarily a r X i v : . [ phy s i c s . c o m p - ph ] J un large subspace of interest. We show the self-energy sepa-ration disentangles correlation contributions from differ-ent spatial regions of the system.After deriving the formalism of the self-energy em-bedding and decomposition, we illustrate our methodsnumerically for nitrogen vacancy in large periodic cellsof hBN. This system represents a realistic simulationof a prototypical single-photon emitter. First, wedemonstrate the reduction of the stochastic error for thedefect states in the hBN monolayer. Next, we studythe interaction of the defect in the hBN-graphene het-erostructure.
II. COMPUTING QUASIPARTICLE ENERGIESA. Self-energy and common approximations
The Green’s function ( G ), a time-ordered correlator ofcreation and annihilation field operators, describes thedynamics of an individual quasiparticle. The poles of G fully determine the single-particle excitations (as well asmany other properties).Solving directly for G is often technically challeng-ing (or outright impossible). Alternatively, the Green’sfunction is often sought via a perturbative expansion ofthe electron-electron interactions on top of a propaga-tor of non-interacting particles (i.e., the non-interactingGreen’s function, G ). The two quantities are related viathe Dyson equation G − = G − − Σ, where Σ is a self-energy accounting, in principle, for all the many-bodyeffects absent in G .Calculations usually employ only a truncatedexpansion of the self-energy. Despite ongoingdevelopments, the most common approach islimited to the popular GW approximation to theself-energy, which is composed of the non-local exchange(Σ X ) and polarization (Σ P ) terms. In the time-domain,the latter operator is expressed as:Σ P ( r , r (cid:48) , t ) = iG ( r , r (cid:48) , t ) W P ( r , r (cid:48) , t + ) , (1)where W P is the time-ordered polarization potential dueto the time-dependent induced charge density . Thepotential is conveniently expressed using the reduciblepolarizability χ as W P ( r , r (cid:48) , t ) = (cid:90) (cid:90) ν ( r , r (cid:48)(cid:48) ) χ ( r (cid:48)(cid:48) r (cid:48)(cid:48)(cid:48) , t ) ν ( r (cid:48)(cid:48)(cid:48) , r (cid:48) ) d r (cid:48)(cid:48) d r (cid:48)(cid:48)(cid:48) (2)Evaluating the action of W P on individual states is thepractical bottleneck of the GW approach. Hence, despiteEq. 1 requires a self-consistent solution, it is commonlycomputed only as a “one-shot” correction. In practice,Eq. 1 thus contains only G and W . Both quantitiesare constructed from the mean-field Hamiltonian, H ,comprising one-body terms and local Hartree, ionic, andexchange-correlation potentials. For W P , it is common to employ the random phase approximation (RPA). BeyondRPA approaches are more expensive and, in general, donot improve the QP energies unless higher-order (vertex)terms are included in Σ. Within this one-shot framework, the QP energiesbecome ε QP = ε + (cid:104) φ | Σ X + Σ P ( ω = ε QP ) − v xc | φ (cid:105) , (3)where ε are eigenvalues of H and v xc is the (approxi-mate) exchange-correlation potential. In Eq. 3, Σ P is inthe frequency domain. B. The stochastic approach to the self-energy
The stochastic G W method seeks the QP energy viarandom sampling of wavefunctions and decomposition ofoperators in the real-time domain. The expectationvalue of Σ is expressed as a statistical estimator. Theresult is subject to fluctuations that decrease with thenumber of samples as 1 / √ N .In this formalism, the polarization self-energy expres-sion is separable. Specifically, for a particular H eigenstate φ , the perturbative correction becomes: (cid:104) φ | Σ P ( t ) | φ (cid:105) = (cid:104) φ | iG ( t ) W P ( t ) | φ (cid:105)(cid:39) N ¯ ζ (cid:88) ¯ ζ (cid:90) φ ( r ) ζ ( r , t ) u ζ ( r , t ) d r , (4)where u ζ ( r, t ) is an induced charge density potential, and ζ is a random vector used for sampling of G (discussedin detail in the next section). The state ζ at time t is | ζ ( t ) (cid:105) ≡ U ,t P µ | ζ (cid:105) , (5)where the U ,t is time evolution operator U ,t ≡ e − iH t . (6)The projector P µ selects the states above or belowthe chemical potential, µ , depending on the sign of t .In practice, P µ is directly related to the Fermi-Diracoperator. Since the Green’s function is a time-ordered quantity, the vectors in the occupied and un-occupied subspace are propagated backward or forwardin time and contribute selectively to the hole and particlenon-interacting Green’s functions.The induced potential u ( r , t ) represents the time-ordered potential of the response to the charge additionor removal: u ( r , t ) = (cid:90) W P ( r , r (cid:48) , t )¯ ζ ( r (cid:48) ) φ ( r (cid:48) ) d r (cid:48) , (7)where ¯ ζ spans the entire Hilbert space. In practice, we compute u from the retarded responsepotential, which is ˜ u = (cid:82) ˜ W P ( r , r (cid:48) , t )¯ ζ ( r (cid:48) ) φ ( r (cid:48) ) d r (cid:48) ; thetime-ordering step affects only the imaginary componentsof the Fourier transforms of u and ˜ u . The retarded response, ˜ u , is directly related to thetime-evolved charge density δn ( (cid:126)r, t ) ≡ n ( r , t ) − n ( (cid:126)r, t = 0)induced by a perturbing potential δv :˜ u ( r , t ) = (cid:90) (cid:90) (cid:90) ν ( r , r (cid:48)(cid:48) ) χ ( r (cid:48)(cid:48) , r (cid:48)(cid:48)(cid:48) , t ) δv ( r (cid:48)(cid:48)(cid:48) , r (cid:48) ) d r (cid:48) d r (cid:48)(cid:48) d r (cid:48)(cid:48)(cid:48) ≡ (cid:90) ν ( r , r (cid:48) ) δn ( r (cid:48) , t )d r (cid:48) (8)where we define a perturbing potential δv = ν ( r , r (cid:48) )¯ ζ ( r (cid:48) ) φ ( r (cid:48) ) . (9)Note that δv is explicitly dependent on the state φ andthe ¯ ζ vector; the latter is part of the stochastically de-composed G operator.The stochastic formalism further reduces the cost ofevaluating ˜ u . Instead of computing δn ( (cid:126)r, t ) by a sumover single-particle states, we use another set of randomvectors { η } confined to the occupied subspace. Time-dependent density n ( r , t ) thus becomes n ( r , t ) = lim N η →∞ N η N η (cid:88) l | η l ( r , t ) | , (10)where η l is propagated in time using U ,t , and H in Eq. 6is implicitly time-dependent. As common, we resortto the density functional theory (DFT) starting point.Since H is therefore a functional of n ( (cid:126)r, t ), the sameholds for the time evolution operator: | η ( t ) (cid:105) = U ,t [ n ( t )] | η (cid:105) . (11)Further, we employ RPA when computing ˜ u ; this cor-responds evolution within the time-dependent Hartreeapproximation. Practical calculations use only a limited number of ran-dom states. Consequently, the time evolved density ex-hibits random fluctuations at each space-time point. Toresolve the response to δv , we use a two-step propaga-tion whose difference is the δn that typically convergesfast with N η . III. EMBEDDED DETERMINISTIC SUBSPACE
The stochastic vectors { ζ } and { η } , introduced in thepreceding section, are constructed on a real-space gridand sample the occupied (or unoccupied) states. Thenumber of these vectors ( N ζ and N η ) is increased so thatthe statistical errors are below a predefined threshold.Further, the underlying assumption is that { ζ } and { η } sample the Hilbert space uniformly.Here, we present a stochastic approach restricted onlyto a subset of states, while selected orbitals, { φ } , aretreated explicitly and constitute an embedded subspace.We denote this set as the { φ } -subspace. In the context of the G W approximation, we use the hybrid approachfor (i) the Green’s function, (ii) the induced potential, or(iii) both G and u simultaneously. In the following, wepresent each case separately. a. First, the non-interacting Green’s function is de-composed into two parts (omitting for brevity the space-time coordinates): G ≡ ˜ G + G φ , (12)where G φ is the Green’s function of the constructed ex-plicitly from { φ } as G φ ( r , r (cid:48) , t ) = (cid:88) j ∈{ φ } ( − θ ( t ) iφ j ( r ) φ ∗ j ( r (cid:48) ) e − iε j t (13)where θ is the Heaviside step function responsible for thetime-ordering (corresponding to particle and hole contri-butions to G ). The complementary part, ˜ G , is sampledwith random states as in Eq. 4.Unlike in the fully stochastic approach (where no G φ term is present), the sampling vectors are constructed asorthogonal to the { φ } -subspace, i.e.: (cid:12)(cid:12) ¯ ζ (cid:11) = (1 − P φ ) (cid:12)(cid:12) ¯ ζ (cid:11) . (14)Here, ¯ ζ spans in principle, the entire Hilbert space, and P φ is: P φ = (cid:88) j ∈{ φ } | φ j (cid:105)(cid:104) φ j | . (15)The construction of ˜ G remains the same as in the fullystochastic case: ˜ G is decomposed by a pair of vectors ¯ ζ and ζ ( t ), cf. Eqs. (14), (5) and (6).Note that it is possible to generalize the projector P φ , Eq.(15), to an arbitrary { φ } -subspace. The particu-lar choice of φ does not affect the decomposition of theGreen’s function, Eq. (12), or the time evolution of ζ ( t ).However, the dynamics of G φ would require explicit ac-tion of U ,t on φ that is, in principle, not an eigenstateof H . We do not pursue this route here and select { φ } -subspace composed from the starting point eigenstates. b. Second, the retarded induced potential ˜ u is de-composed through partitioning of the time-dependentcharge density, cf., Eq. (8), (omitting the space-time co-ordinates): n ≡ ˜ n + n φ , (16)where the n φ is the density constructed from occupiedstates φ (which we assume to be mutually orthogonal): n φ ( r , t ) = (cid:88) j ∈{ φ } f j | φ j ( r , t ) | . (17)Here, f j is the occupation of the j th state. Note thatchoosing φ within the unoccupied subspace is meaning-less in this context. The complementary part, ˜ n , is givenby stochastic sampling, Eq. (10), that employs randomvectors ˜ η : | (cid:101) η (cid:105) = (cid:16) − ˜ P φ (cid:17) | η (cid:105) , (18)where η spans the entire occupied subspace and˜ P φ = (cid:88) j ∈{ φ } f j | φ j (cid:105)(cid:104) φ j | . (19)The time propagation of the charge density is simi-lar to the fully stochastic case. The deterministic andstochastic vectors evolve as: | φ j ( t ) (cid:105) = U ,t [ n φ ( t ) , ˜ n ( t )] | φ j (cid:105) (20) | ˜ η ( t ) (cid:105) = U ,t [ n φ ( t ) , ˜ n ( t )] | ˜ η (cid:105) , (21)where U ,t is explicitly expressed as a functional of thetwo density contributions from Eq. (16). Note that theseexpressions are analogous to Eq. (11). c. Third, both partitionings are used in conjunc-tion. While G contains contributions from both occu-pied and unoccupied states, only the former are includedin ˜ u . The combined partitioning may use entirely differ-ent subspaces for the Green’s function and the inducedpotential. In Section V C, we employ the decompositionin both terms because it yields the best results and sig-nificantly reduces statistical fluctuations. IV. DECOMPOSITION OF THE STOCHASTICPOLARIZATION SELF-ENERGY
In the preceding section, we partition retarded inducedpotential and the Green’s function, intending to decreasestochastic fluctuations in the self-energy. Here, we usethe partitioning to achieve the second goal of this paper– to determine the contribution to Σ P ( ω ).Conceptually, we want to address quasiparticle scat-tering by correlations from a particular subspace. Inthe expression for Σ P ( ω ), Eq. (4), this corresponds toaccounting for selected charge density fluctuations in ˜ u .In practice, we construct the subspace polarization self-energy as: (cid:104) φ | Σ sP ( t ) | φ (cid:105) = 1 N ¯ ζ (cid:88) ¯ ζ (cid:90) φ ( r ) ζ ( r , t ) u sζ ( r , t ) d r , (22)where we introduced the subspace induced potential u sζ which is obtained from its retarded form (in analogy toEq. (8)): ˜ u sζ ( r , t ) = (cid:90) ν ( r , r (cid:48) ) δn s ( r (cid:48) , t )d r (cid:48) . (23)This potential stems from the induced charge densitythat includes contributions only from selected orbitals { φ } . Note that n s ( r (cid:48) , t ) is obtained either from individ-ual single-particle states or from the stochastic samplingof the { φ } -subspace under consideration. If the set of φ states is large, it is natural to employ thestochastic approach; the density is sampled according toEq. (10) with vectors η s prepared as: | η s (cid:105) = ˜ P φ | η (cid:105) (24)where the projector is in Eq. (19) and vectors η span theentire occupied subspace.The time evolution of η s follows Eq. (11), i.e., it isgoverned by the operator U ,t that depends on the total time-dependent density: | η s ( t ) (cid:105) = U ,t [ n ( t )] | η s (cid:105) . (25)Hence, despite Σ sP ( t ) contains only fluctuations from aparticular subspace, the calculation requires knowledgeof the time evolution of both n s and n .In practice, we employ a set of two independentstochastic samplings: (i) vectors | η (cid:105) describing the en-tire occupied space, and (ii) vectors | η s (cid:105) confined only tothe chosen { φ } -subspace. The first set characterizes thetotal change density fluctuation and enters U ,t in Eq. 25. V. NUMERICAL RESULTS AND DISCUSSIONA. Computational details
In this section, we will demonstrate the capabilitiesof the method introduced above. The starting-pointcalculations are performed with a real-space DFT im-plementation, employing regular grids, Troullier-Martinspseudopotentials, and the PBE functional for ex-change and correlation. We investigate finite and 2Dinfinite systems using modified periodic boundary condi-tions with Coulomb interaction cutoffs. The numerical verification for the SiH molecule is insection V B. To converge the occupied H eigenvalues to < E h and64 × ×
64 real-space grid with the step of 0.3 a .Large calculations for the V N defect in hBN monolayerand in hBN heterostructure with graphene are in sectionsV C and V D. In both cases, we consider relaxed rectan-gular 12 × together with Tkatchenko-Scheffler’s total energycorrections. The heterostructure is built with an inter-layer distance of 3.35 ˚A.The GW calculations were performed using a develop-ment version of the StochasticGW code. The calcu-lations employ an additional set of 20,000 random vectorsused in the sparse stochastic compression used for time-ordering of ˜ u . The time propagation of the inducedcharge density is performed for maximum propagationtime of 50 a.u., with the time-step of 0.05 a.u. FIG. 1. Verification for the SiH molecule. (a) The compar-ison of the total self-energies of the highest occupied state(HOMO) obtained with reference and embedding schemes.The figure demonstrates that all five approaches yield identi-cal Σ( ω ). In the plot, G denotes the decomposition schemewhere HOMO is embedded in the Green’s function; W de-notes the calculation with HOMO state embedded in screenedCoulomb potential; G+W corresponds to simultaneous de-composition of G and W. The red-dashed line represents thesum over subspace contributions shown in the panel below.(b) Subspace self-energy of the HOMO state that containsonly a contribution of a specific state i , denoted as Σ i in thefigure’s legend. B. Verification using molecular states
We verify the implementation of the methods pre-sented in Sections III and IV on SiH molecule. In partic-ular, we test separately : (i) the construction of the em-bedding schemes, (ii) decomposition of the time propaga-tion, and (iii) the evaluation of the subspace self-energy.The reference calculation employs only one level ofstochastic sampling (for the Green’s function, while therest of the calculation is deterministic). For small sys-tems, this approach converges fast as the stochastic fluc-tuations are small. Yet, we need N ζ = 1 ,
500 vectorsto decrease the QP energies errors below 0.08 eV. An il-lustration of the self-energy for the HOMO state of SiH is in Fig. 1; everywhere in the figure, the stochastic erroris below 0 .
13 eV for all frequencies.We first inspect the results for an embedded determin-istic subspace. Fig. 1 shows that the different schemes forthe explicit treatment of the HOMO state, Section III,produce the same self-energy curve. The inclusion ofHOMO in W is combined with the stochastic samplingof the three remaining orbitals by N η =16 random vec-tors. Note that it is not economical to sample the actionof W P for small systems, and this calculation servesonly as a test case. This treatment yields a statistical er-ror of 0.08 eV, i.e., the same as the reference calculationdespite the additional fluctuations due to the η vectors. When the HOMO orbital is explicitly included in G ,the resulting statistical error is below the error of thereference calculation (0.05 eV). Such a result is expectedsince the reference relies on the stochastic sampling of theGreen’s function. The same happens when the frontierorbital is in both G and W ( N η =16 random vectorssample the induced charge density). Tests for other statesare not presented here but lead to identical conclusions.Next, we verify that the density entering the time-evolution operator U ,t [ n ( t )] can be constructed fromdifferent states than it is acting on. Namely, the in-duced charge and the time-dependent densities may besampled and built by different means. To demonstratethis, we propagate each of the H eigenstates with U ,t [ n ( t )], where n ( t ) is stochastically sampled by N η =16random vectors. Only the induced charge density, enter-ing Eq. (8), is computed from the { φ } eigenvectors. Theagreement with the reference self-energy curve is excel-lent with differences smaller than the standard deviationsat each frequency point; see Fig. 1a.Finally, we inspect the subspace self-energy in which U ,t [ n ( t )] employs the total charge density sampled by N η =16 random vectors. In Fig. 1b shows four differentΣ s ( ω ) curves corresponding to the contributions of indi-vidual H eigenstates. Since the eigenstates are orthogo-nal, the total self-energy is simply the sum of individualΣ s ( ω ) components. The additivity of the subspace self-energies is demonstrated numerically in Fig. 1a. The sub-space results illustrate that HOMO and the bottom va-lence orbital exhibit the largest amplitudes of Σ sP ; hence,these two states dominate the correlation near the ion-ization edge. C. Deterministic treatment of localized states
The deterministic subspace embedding should numer-ically stabilize the stochastic sampling and decreases thecomputational cost. To test the methodology on a real-istic system, we consider the electronic structure of aninfinitely periodic hBN monolayer containing a single ni-trogen vacancy ( V N ) per a unit cell with dimensions of3 . × . The relaxed monolayer geometry shows only mild re-structuring. The vacancy introduces three localizedstates with C and C spatial symmetry. The former( C ) is singly occupied and forms an in-gap state. Thelatter state is doubly degenerate and pushed high in theconduction region. Due to the enforced spin-degeneracy,the electron-electron interactions are increased and C appears higher than in the previous calculations .The C and C single-particle wavefunctions are illus-trated in Fig. 2.We first focus on the delocalized top valence and bot-tom conduction states. The fully stochastic calculationsconverge fast for both of them. The charge density fluc-tuations are sampled by N η = 8 vectors, and the Green’sfunction requires N ζ =1500 to yield QP energies with sta-tistical errors below 0.03 eV for the valence band max-imum (VBM). The error for the conduction band mini-mum (CBM) is < .
01 eV. The computed resulting quasi-particle band-gap is 6 . ± .
04 eV, in excellent agreementwith previous calculations and experiments (providing arange of values between 6.1 and 6.6 eV).
To investigate the electronic structure in greater detail,we employ the projector-based energy-momentum analy-sis based on supercell band unfolding.
In practice,the individual wavefunctions within our simulation cellare projected onto the Brillouin zone of a single hBN unitcell. The resulting bandstructure is shown in Fig. 2a.Since our calculations employ a rectangular supercells,the critical point K of the hexagonal Brillouin zone ap-pears between the Γ and X points (marked on the hori-zontal axis by (cid:63) ). The figure shows that the fundamentalband gap is indirect; the direct transition is 7 . ± . The defect states appear as flat bands, labeled in Fig. 2by their symmetry. As expected from the outset, thestochastic calculations for C exhibit large fluctuations.While each sample is numerically stable, random vectorsproduce a self-energy curve with a significant statisticalerror at each frequency point, i.e., the time evolution isstrongly dependent on the initial choice of { ζ } and { η } .In this example, only the C state exhibits such behavior.In Fig. 2 illustrates the self-energies for the band edgeand the defect states. For all the cases, the plots showthe spread of the Σ P ( ω ) curves: these correspond to theouter envelope for 15 distinct calculations, each employ-ing 100 ζ sampling vectors (combined with N η = 8 each).For the C defect state, the stochastic sampling is pos-sibly biased. The variation is three times as big as thespread of the VBM, and almost seven times larger com-pared to CBM. Away from the QP energy, the fluctu-ations increase even further; the spread of the samplesbecomes two times larger near the maximum of Σ P at −
20 eV. The convergence of the QP energy is poor, andthe low sample standard deviation (roughly twice as bigas for VBM) suggests incorrect statistics. In practice,each sampling yields a self-energy curve that lies outsideof the standard deviation of the previous simulation.The deterministic embedding remedies insufficientsampling without increasing the computational cost.Naturally, we select the C defect state and treat it explicitly (while randomly sampling the rest of the or-bitals). Hence, according to the notation of Section III,the { φ } -subspace contains only a single orbital.The decomposition of the non-interacting Green’sfunction follows Eq. (12) and stabilizes the sampling.The spread of the self-energy curves decreases approxi-mately three times for a wide range of frequencies. Withthe embedded C state, the statistical error of the QP en-ergy decreases smoothly and uniformly with the numberof samples. Each new sampling falls within the error ofthe calculations. Yet, the final statistical error (0.03 eV)is less than 10% larger than for the delocalized states.The decomposition of the induced charge density aloneis less promising. Fundamentally, δn ( r , t ) contains contri-butions from the entire system, and a small { φ } -subspacewill unlikely lead to drastic improvement. Indeed, thestatistical errors and convergence behavior remain thesame as for the fully stochastic treatment.Embedding of the localized state in both G and δn ( r , t ) is, however, the best strategy. If both decomposi-tions use the same (or overlapping) { φ } -subspace, the in-duced charge density and the potential δv , Eq. (9), share(at least some) φ states. Consequently, the sampling ofΣ P ( ω ) becomes less dependent on the particular choiceof ζ vectors, which sample only states orthogonal to { φ } .Indeed, the embedding of a single localized state in G and W results in a nearly four-fold reduction of the statisti-cal fluctuation that translates to more than an order ofmagnitude savings in the computational time. The errorof the QP energy of the C state is approximately half ofthe error for VBM (16 meV for N ζ = 1 , D. Stochastic subspace self-energy
Having an improved description of the localized statesin hand, we now turn to the decomposition of the self-energy. The real-time approach described in Section IVallows inspecting the many-body interactions from a se-lected portion of the system. In contrast to the previoussubsection, the subspace of interest contains a large num-ber of states (irrespective of their degree of localization),and it is randomly sampled. The goal of this decompo-sition is to understand the role of correlation, especiallyat interfaces.To test our method, we investigate a periodic hBNmonolayer containing a single V N defect placed ongraphene. Such heterostructure has also been realizedexperimentally. The structure contains 2299 valenceelectrons, and it is illustrated in Fig. 3 together with se-lected orbital isosurfaces. The C state is energeticallylower than C (as in the pristine hBN monolayer). Both FIG. 2. Deterministic embedding for the nitrogen vacancy inhBN monolayer. (a) Black points represent individual statesthe unfolded band structure; blue lines are a guide for theeyes. Panels b, c, and d depict the self-energy curves (solidline) for C , C , and VBM state respectively; the states aremarked in panel a by red dashed rectangles. The statisticalerrors are smaller than the thickness of the line. Grey shadedareas (on each plot) correspond to the “spread” of self-energycurves for 15 distinct fully stochastic calculations; solid bluelines are the statistical averages. The light red area in panel crepresents the spread after deterministic embedding of the C state; the solid red line is the self-energy curve. The third col-umn contains the electron density of the corresponding states. defects only weakly hybridize with graphene and remainlocalized within the hBN sublayer. However, graphenepresence leads to a slightly increased delocalization ofthe defects within the monolayer. By unfolding the wave functions of the bilayer, we ob-tain the band structure shown in Fig. 3a. The grapheneand hBN bands are distinguished by a state projectionon the densities above and below the center of the in-terlayer region. Note that this simple approach captureshybridized states as having dual character (i.e., they ap-pear having both graphene and hBN contributions).The graphene portion of the band structure reproducesthe well-known semimetallic features with a Dirac pointlocated at the K boundary of the hexagonal Brillouinzone. As discussed above, the K appears in between Γand X of the rectangular cell, and it is labeled by (cid:63) inFig. 3. The typical Dirac cone dispersion of graphene isonly a little affected by the hBN presence. However,the ordering of the hBN and graphene states is non-trivial. In contrast to previous DFT calculations for3-times higher defect density, we notice that the Diracpoint remains close to the Fermi level despite the chargetransfer from the C defect state. This is not surpris-ing: sparse charge defects lead to only weak doping of graphene. Further, the previous DFT results employedsmall supercells and may suffer from significant electronic“overdelocalization” that spuriously enhances chargetransfer.The hBN part of the band structure, Fig. 3, is onlyweakly affected by the heterostructure formation. Thedelocalized states are qualitatively identical to those inthe monolayer, Fig. 2. The fundamental band-gap re-mains indirect and reduced to 6 . ± .
04 eV. The screen-ing introduced by graphene thus leads to a small changein E g ( < . sP ( ω ) using Eq. (22). Here, the { φ } -subspace isconstructed from stochastic samples of the 576 occupiedgraphene states (distinguished by green color in Fig. 3).For each sampling of the Green’s function, the subspaceis described by eight random vectors in the { φ } -subspace.Additional eight vectors sample the complementary sub-space.Fig. 3 shows the decomposition of the self-energy, indi-cating the contribution of graphene. For the delocalizedstates, the Σ P ( ω ) curves appear similar to those in thehBN monolayer. Specifically, we see that Σ P ( ω ) for VBMhas practically identical frequency dependence between-15 and 0 eV. While in-plane contributions from hBNdominate the entire curve, screening from the graphenesubstrate is significant. A na¨ıve comparison between themonolayer hBN and the heterostructure suggests that theenhanced maximum in ReΣ P at −
25 eV is caused extrin-sically by graphene. Surprisingly, this is not the case,and the effect is only indirect: the presence of grapheneleads to shifting of the hBN spectral features. At the QPenergy, graphene contributes to the correlation by 36%.The situation is different for CBM. Towards the staticlimit ( ω → The self-energy curves of the C and C defect statesare significantly different from those computed for themonolayer. In both cases, we observe a negative shiftrelated to the QP stabilization by non-local correla-tions. The Σ P contributions to the defects’ QP ener- FIG. 3. Excitations in an hBN-graphene heterostructure containing a nitrogen vacancy. (a) Calculated band structure – thestates localized on graphene are in green, while blue points represent hBN. Lines serve as a guide for the eyes. The individualisosurfaces are shown on the right for VBM, C defect state, and C defect states. Second row contains averaged self-energycurves of the: (b) VBM, (c) C defect state, (d) C defect state, (e) CBM. Solid red lines are the full heterostructure self-energies; green and blue dashed lines are respectively the Σ sP ( ω ) graphene monolayer and hBN monolayer subspace self-energies.The black line represents the self-energy of an isolated hBN monolayer with V N . The third row of panels contains the ImΣ sP ( ω )of graphene. Dashed vertical lines indicate the features in the experimental EELS spectrum from Ref. 60. The experimentallyobserved features, labeled as i, ii, and iii, are attributed to π → π ∗ interband transitions, π , and π + σ plasmon excitations,respectively. gies are roughly three times larger in the heterostruc-ture. Na¨ıvely, we expect that the localized states arestrongly coupled to graphene polarization modes, whichlead to negative Σ P . Surprisingly, the intrinsic hBN in-teractions are the major driving force. The defect statesare mostly affected by the in-plane induced density fluc-tuations (constituting 64% and 76% of the total polariza-tion self-energy for C and C ). Like the VBM, grapheneindirectly acts on the defect QP states by enhancing thein-plane fluctuations, rather than direct coupling to thedefects.For all the states illustrated in Fig. 3, the graphenesubspace contributions are governed by density oscilla-tions induced by electron removal from or addition tothe hBN layer. Collective charge excitations, i.e., plas-mons, naturally correspond to the most prominent fea-tures. They are the poles of the screened Coulomb inter-action and appear as strong peaks in the imaginary partof the subspace self-energy (bottom panels in Fig. 3).Conceptually, the substrate plasmons are directly re-lated to QP energy dissipation. Due to the weak couplingbetween the monolayers, ImΣ sP ( ω ) of the substrate showsthe same spectral features at the experimental electronenergy loss spectra measured for graphene alone. As dis- tinct states couple to individual plasmon modes differ-ently, the intensity of peaks in ImΣ sP ( ω ) varies. Dueto changes in geometry and the presence of hBN, thespectral features are shifted from the experimental data.Still, we identify π + σ and π plasmons. The former isuniversally the most prominent feature in the subspaceself-energy. The π (possibly together with π → π ∗ ) plas-mon contributes much less, and it is appreciable only forthe C defect state. We surmise that the additional peaknear zero frequency corresponds to the “static” polar-ization of graphene due to the charge transfer from thedefect.As shown in this example, the stochastic subspace self-energy efficiently probes dynamical electron-electron in-teractions at interfaces. The random sampling allows se-lecting an arbitrary subspace to identify its contributionto the correlation energy and the excited state lifetimes. VI. CONCLUSIONS
Stochastic Green’s function approaches represent aclass of low-scaling many-body methods based onstochastic decomposition and random sampling of theHilbert space. The overall computational cost is im-pacted by the presence of localized states that increasethe sampling errors. Here, we have introduced a practi-cal solution for the G W method. The localized statesare treated explicitly (deterministically), while the rest issubject to stochastic sampling. We have further shownthat the subspace-separation can be applied to decom-pose the dynamical self-energy into contributions fromdistinct states (or parts of the systems).Using nitrogen vacancy in hBN monolayer, we demon-strate that the deterministic embedding dramatically re-duces the statistical fluctuations. Consequently, the com-putational time decreases by more than an order of mag-nitude for a given target level of stochastic error. Theembedding is made within the Green’s function and thescreened Coulomb interactions independently or simul-taneously; the latter provides the best results. For de-localized states, embedding is not necessary and per-forms similar to the standard uniform sampling with real-space random vectors. The embedding scheme presentedhere is general, and it is the first step towards hybridtechniques that treat selected subspace at the distinct(higher) level of theory. The development of hybrid tech-niques is currently underway.We further demonstrate that the subspace self-energycontains (in principle additive) contributions from the in-duced charge density. The self-energy contribution froma particular portion of the system is computed by con-fining sampling vectors into different (orthogonal) sub-spaces. We exemplify the capabilities of such calculationson defects state in the hBN-graphene heterostructure.Here, the charged excitation in one layer directly couplesto density oscillations in the substrate monolayer. Theresponse is dominated by plasmon modes, which are iden-tified as features in the imaginary part of the subspaceself-energy. Surprisingly, the electronic correlations forthe defects are governed by the interactions in the hostlayer; the substrate only indirectly affects their strength.This example serves as a stimulus for additional studyof the defect states in heterostructures. Here, the cou-pling strength between individual subsystems is tunableby particular stacking order and induced strains. Thesubspace self-energy represents a direct route to exploresuch quantum many-body interfacial phenomena. ACKNOWLEDGMENTS
This work was supported by the NSF through the Ma-terials Research Science and Engineering Centers (MR-SEC) Program of the NSF through Grant No. DMR-1720256 (Seed Program). In part, this work wassupported by the NSF Quantum Foundry through Q-AMASE-i program Award No. DMR-1906325. The cal-culations were performed as part of the XSEDE com-putational Project No. TG-CHE180051. Use was madeof computational facilities purchased with funds from theNational Science Foundation (CNS-1725797) and admin- istered by the Center for Scientific Computing (CSC).The CSC is supported by the California NanoSystemsInstitute and the Materials Research Science and En-gineering Center (MRSEC; NSF DMR-1720256) at UCSanta Barbara. H. Yu, G.-B. Liu, J. Tang, X. Xu, and W. Yao, “Moir´e ex-citons: From programmable quantum emitter arrays to spin-orbit–coupled artificial lattices,” Science Advances (2017),10.1126/sciadv.1701696. G. Grosso, H. Moon, B. Lienhard, S. Ali, D. K. Efetov, M. M.Furchi, P. Jarillo-Herrero, M. J. Ford, I. Aharonovich, andD. Englund, “Tunable and high-purity room temperature single-photon emission from atomic defects in hexagonal boron nitride,”Nature Communications , 705 (2017). M. Yankowitz, J. Jung, E. Laksono, N. Leconte, B. L. Chittari,K. Watanabe, T. Taniguchi, S. Adam, D. Graf, and C. R. Dean,“Dynamic band-structure tuning of graphene moir´e superlatticeswith pressure,” Nature , 404–408 (2018). Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo,J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxi-ras, R. C. Ashoori, and P. Jarillo-Herrero, “Correlated insulatorbehaviour at half-filling in magic-angle graphene superlattices,”Nature , 80–84 (2018). U. Zondiner, A. Rozen, D. Rodan-Legrain, Y. Cao, R. Queiroz,T. Taniguchi, K. Watanabe, Y. Oreg, F. von Oppen, A. Stern,E. Berg, P. Jarillo-Herrero, and S. Ilani, “Cascade of phase tran-sitions and dirac revivals in magic-angle graphene,” Nature ,203–208 (2020). M. E. Turiansky, A. Alkauskas, and C. G. Van de Walle, “Spin-ning up quantum defects in 2d materials,” Nature Materials ,487–489 (2020). A. Gottscholl, M. Kianinia, V. Soltamov, S. Orlinskii, G. Mamin,C. Bradac, C. Kasper, K. Krambrock, A. Sperlich, M. Toth,I. Aharonovich, and V. Dyakonov, “Initialization and read-outof intrinsic spin defects in a van der waals crystal at room tem-perature,” Nature Materials , 540–545 (2020). C. Wang, C. Axline, Y. Y. Gao, T. Brecht, Y. Chu, L. Frunzio,M. H. Devoret, and R. J. Schoelkopf, “Surface participationand dielectric loss in superconducting qubits,” Applied PhysicsLetters , 162601 (2015). D. Y. Qiu, F. H. da Jornada, and S. G. Louie, “Environmen-tal screening effects in 2d materials: Renormalization of thebandgap, electronic structure, and optical spectra of few-layerblack phosphorus,” Nano Letters , 4706–4712 (2017). A. Tartakovskii, “Moir´e or not,” Nature Materials , 581–582(2020). L. Yuan, B. Zheng, J. Kunstmann, T. Brumme, A. B. Kuc,C. Ma, S. Deng, D. Blach, A. Pan, and L. Huang, “Twist-angle-dependent interlayer exciton diffusion in ws2–wse2 heterobilay-ers,” Nature Materials , 617–623 (2020). R. M. Martin, L. Reining, and D. M. Ceperley,
Interacting Elec-trons (Cambridge University Press, 2016). A. L. Fetter and J. D. Walecka,
Quantum Theory of Many-Particle Systems (Dover Publications, 2003). L. Hedin, “New Method for Calculating the One-Particle Green’sFunction with Application to the Electron-Gas Problem,” Phys.Rev. , A796–A823 (1965). M. S. Hybertsen and S. G. Louie, “Electron correlation in semi-conductors and insulators: Band gaps and quasiparticle ener-gies,” Phys. Rev. B , 5390–5413 (1986). F. Aryasetiawan and O. Gunnarsson, “The GW method”ReportsProg. Phys. , 237–312 (1998). M. Govoni and G. Galli, “Large scale gw calculations,” Journalof Chemical Theory and Computation , 2680–2696 (2015). D. Neuhauser, Y. Gao, C. Arntsen, C. Karshenas, E. Rabani,and R. Baer, “Breaking the Theoretical Scaling Limit for Pre-dicting Quasiparticle Energies: The Stochastic G W Approach,”Phys. Rev. Lett. , 076402 (2014). V. c. v. Vlˇcek, W. Li, R. Baer, E. Rabani, and D. Neuhauser,“Swift gw beyond 10,000 electrons using sparse stochastic com-pression,” Phys. Rev. B , 075107 (2018). V. Vlcek, E. Rabani, D. Neuhauser, and R. Baer, “Stochasticgw calculations for molecules,” J. Chem. Theory Comput. ,4997–5003 (2017). J. Brooks, G. Weng, S. Taylor, and V. Vlcek, “Stochastic many-body perturbation theory for moir´e states in twisted bilayer phos-phorene,” Journal of Physics: Condensed Matter , 234001(2020). C. Li, Z.-Q. Xu, N. Mendelson, M. Kianinia, M. Toth, andI. Aharonovich, “Purification of single-photon emission from hbnusing post-processing treatments,” Nanophotonics , 2049 – 2055(2019). T. T. Tran, K. Bray, M. J. Ford, M. Toth, and I. Aharonovich,“Quantum emission from hexagonal boron nitride monolayers,”Nature Nanotechnology , 37–41 (2016). A. L. Exarhos, D. A. Hopper, R. R. Grote, A. Alkauskas, andL. C. Bassett, “Optical signatures of quantum emitters in sus-pended hexagonal boron nitride,” ACS Nano , 3328–3336(2017). N. Mendelson, Z.-Q. Xu, T. T. Tran, M. Kianinia, J. Scott,C. Bradac, I. Aharonovich, and M. Toth, “Engineering and tun-ing of quantum emitters in few-layer hexagonal boron nitride,”ACS Nano , 3132–3140 (2019). E. Maggio and G. Kresse, “GW vertex corrected calculations formolecular systems,” J. Chem. Theory Comput. , 4765 (2017). M. Hellgren, N. Colonna, and S. de Gironcoli, “Beyond the ran-dom phase approximation with a local exchange vertex,” Phys.Rev. A , 045117 (2018). V. Vlcek, “Stochastic vertex corrections: Linear scaling methodsfor accurate quasiparticle energies,” Journal of Chemical Theoryand Computation , 6254–6266 (2019). A. M. Lewis and T. C. Berkelbach, “Vertex corrections to thepolarizability do not improve the gw approximation for the ion-ization potential of molecules,” Journal of Chemical Theory andComputation , 2925–2932 (2019). R. Baer and D. Neuhauser, “Real-time linear response for time-dependent density-functional theory,” The Journal of chemicalphysics , 9803–9807 (2004). Y. Gao, D. Neuhauser, R. Baer, and E. Rabani, “Sublinearscaling for time-dependent stochastic density functional theory,”The Journal of Chemical Physics , 034106 (2015). D. Neuhauser, E. Rabani, Y. Cytter, and R. Baer, “Stochas-tic optimally tuned range-separated hybrid density functionaltheory,” The Journal of Physical Chemistry A , 3071–3078(2016). E. Rabani, R. Baer, and D. Neuhauser, “Time-dependentstochastic Bethe-Salpeter approach,” Phys. Rev. B , 235302(2015). S. Baroni, S. de Gironcoli, and A. Dal Corso, “Phonons andrelated crystal properties from density-functional perturbationtheory,” Rev. Mod. Phys. , 515–562 (2001). D. Neuhauser and R. Baer, “Efficient linear-response method cir-cumventing the exchange-correlation kernel: theory for molecu-lar conductance under finite bias.” J. Chem. Phys. , 204105(2005). N. Troullier and J. L. Martins, “Efficient pseudopotentials forplane-wave calculations,” Phys. Rev. B , 1993 (1991). J. P. Perdew and Y. Wang, “Accurate and simple analytic rep-resentation of the electron-gas correlation energy,” Phys. Rev. B , 13244–13249 (1992). C. A. Rozzi, D. Varsano, A. Marini, E. K. U. Gross, and A. Ru-bio, “Exact coulomb cutoff technique for supercell calculations,”Phys. Rev. B , 205119 (2006). P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B.Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Co-coccioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli,P. Delugas, R. A. D. Jr, A. Ferretti, A. Floris, G. Fratesi, G. Fu-gallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Kkbenli, M. Lazzeri,M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen,A. O. de-la Roza, L. Paulatto, S. Ponc, D. Rocca, R. Sabatini,B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Tim-rov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Ba-roni, “Advanced capabilities for materials modelling with quan-tum espresso,” Journal of Physics: Condensed Matter , 465901(2017). A. Tkatchenko and M. Scheffler, “Accurate molecular van derwaals interactions from ground-state electron density and free-atom reference data,” Phys. Rev. Lett. , 073005 (2009). S. Park, C. Park, and G. Kim, “Interlayer coupling enhancementin graphene/hexagonal boron nitride heterostructures by inter-calated defects or vacancies,” The Journal of Chemical Physics , 134706 (2014). C. Attaccalite, M. Bockstedte, A. Marini, A. Rubio, andL. Wirtz, “Coupling of excitons and defect states in boron-nitridenanostructures,” Phys. Rev. B , 144115 (2011). B. Huang and H. Lee, “Defect and impurity properties of hexag-onal boron nitride: A first-principles calculation,” Phys. Rev. B , 245406 (2012). N. L. McDougall, J. G. Partridge, R. J. Nicholls, S. P. Russo,and D. G. McCulloch, “Influence of point defects on the nearedge structure of hexagonal boron nitride,” Phys. Rev. B ,144106 (2017). M. E. Levinshtein, S. L. Rumyantsev, and M. S. Shur, “Prop-erties of advanced semiconductor materials : Gan, aln, inn, bn,sic, sige,” (2001). F. Fuchs, J. Furthm¨uller, F. Bechstedt, M. Shishkin, andG. Kresse, “Quasiparticle band structure based on a generalizedkohn-sham scheme,” Phys. Rev. B , 115109 (2007). F. H¨user, T. Olsen, and K. S. Thygesen, “Quasiparticle gw cal-culations for solids, molecules, and two-dimensional materials,”Phys. Rev. B , 235132 (2013). V. Popescu and A. Zunger, “Extracting e versus p(cid:126)k effective bandstructure from supercell calculations on alloys and impurities,”Phys. Rev. B , 085201 (2012). H. Huang, F. Zheng, P. Zhang, J. Wu, B.-L. Gu, and W. Duan,“A general group theoretical method to unfold band structuresand its application,” New Journal of Physics , 033034 (2014). P. V. C. Medeiros, S. Stafstr¨om, and J. Bj¨ork, “Effects of ex-trinsic and intrinsic perturbations on the electronic structure ofgraphene: Retaining an effective primitive cell band structure byband unfolding,” Phys. Rev. B , 041407 (2014). P. Cudazzo, L. Sponza, C. Giorgetti, L. Reining, F. Sottile, andM. Gatti, “Exciton band structure in two-dimensional materi-als,” Phys. Rev. Lett. , 066803 (2016). F. Paleari, T. Galvani, H. Amara, F. Ducastelle, A. Molina-S´anchez, and L. Wirtz, “Excitons in few-layer hexagonal boronnitride: Davydov splitting and surface localization,” 2D Materi-als , 045017 (2018). Z.-Q. Xu, N. Mendelson, J. A. Scott, C. Li, I. H. Abidi, H. Liu,Z. Luo, I. Aharonovich, and M. Toth, “Charge and energy trans-fer of quantum emitters in 2d heterostructures,” 2D Materials ,031001 (2020). O. Salihoglu, N. Kakenov, O. Balci, S. Balci, and C. Kocabas,“Graphene as a reversible and spectrally selective fluorescencequencher,” Scientific Reports , 33911 (2016). J. Lee, W. Bao, L. Ju, P. J. Schuck, F. Wang, and A. Weber-Bargioni, “Switching individual quantum dot emission throughelectrically controlling resonant energy transfer to graphene,”Nano Letters , 7115–7119 (2014). C. Bjelkevig, Z. Mi, J. Xiao, P. A. Dowben, L. Wang, W.-N. Mei,and J. A. Kelber, “Electronic structure of a graphene/hexagonal-BN heterostructure grown on ru(0001) by chemical vapor deposi-tion and atomic layer deposition: extrinsically doped graphene,”Journal of Physics: Condensed Matter , 302002 (2010). The distributions of the wavefunctions are governed by the localexternal and Hartree potentials (since the orbitals correspond toeigenstates of the mean-field Hamiltonian H ). P. Mori-S´anchez, A. J. Cohen, and W. Yang, “Localization anddelocalization errors in density functional theory and implica-tions for band-gap prediction,” Phys. Rev. Lett. , 146401(2008). A. J. Cohen, P. Mori-S´anchez, and W. Yang, “Insights into cur-rent limitations of density functional theory,” Science , 792–794 (2008). M. K. Kinyanjui, C. Kramberger, T. Pichler, J. C. Meyer,P. Wachsmuth, G. Benner, and U. Kaiser, “Direct probe of linearly dispersing 2d interband plasmons in a free-standinggraphene monolayer,” EPL (Europhysics Letters) , 57005(2012). J. Towns, T. Cockerill, M. Dahan, I. Foster, K. Gaither,A. Grimshaw, V. Hazlewood, S. Lathrop, D. Lifka, G. D. Peter-son, R. Roskies, J. Scott, and N. Wilkins-Diehr, “Xsede: Accel-erating scientific discovery,” Computing in Science & Engineering16