Decomposition and embedding in the stochastic GW self-energy

First-principles treatment of electron excitation energies is crucial for guiding the development of new materials with tailored optoelectronic properties. Localized quantum states became the focal point for condensed systems^1–7 since the lifetime of localized excitations can be controllably modified by interfacial couplings.^8–11 Besides predicting experimental observables, theoretical investigations thus help elucidate the interplay of the electron–electron interactions and the role of the environment.

The excited electrons and holes near the Fermi level are conveniently described by the quasiparticle (QP) picture: the charge carriers are characterized by renormalized interactions and a finite lifetime limited by energy dissipation, which governs the deexcitation mechanism. Quantitative predictions of QPs necessitate the inclusion of non-local many-body interactions.

The prevalent route to describe QPs in condensed systems employs the Green’s function formalism.^12,13 The excitation energy and its lifetime are inferred from the QP dynamics. The many-body effects are represented by the non-local and dynamical self-energy, Σ. In practice, Σ is approximated by selected classes of interactions, forming a hierarchy of systematically improvable methods.¹² The formalism also allows constructing the self-energy for distinct states with a different form of Σ.

Here, we neglect the vibrational effects, and the induced density fluctuations dominate the response of the system to the excitation. The perturbation expansion is conveniently based on the screened Coulomb interaction, W, which explicitly incorporates the system’s polarizability. The neglect of the higher-order (vertex) terms,¹² canceling the “self-screening” error and capturing two-particles interactions,^14,15 leads to the popular GW method.^12,16–18 This approach predicts the QP gap, ionization potentials, and electron affinities, in good agreement with experiments.^16,18 Within the GW approximation, Σ is a product of W and the Green’s function, G, in the time domain.

Conventional implementations of the GW self-energy scale as O(N⁴) with the system size, and it is usually limited to few-electron systems. Recent developments have significantly reduced the prefactor of the conventional GW implementations.^19–22 Furthermore, a hybrid stochastic–deterministic approach attempted to reduce the cost by constructing the stochastic “pseudobands” when summing over quasi-degenerate states.²³ In contrast, our stochastic treatment recasts the full self-energy expression as a statistical estimator. The method employs sampling of electronic wavefunctions combined with the decomposition of operators.^24–26 In practice, this formulation decreases the computational time significantly and leads to a linear scaling algorithm.^24–26 

The stochastic approach employs real-space random vectors that sample the Hilbert space of the electronic Hamiltonian. The statistical error in Σ is governed by the number of vectors, which are used for the decomposition of G and the evaluation of W. Besides the calculation parameters, i.e., the number of stochastic vectors, the sampling critically depends on how uniformly the Hilbert space is represented. The statistical fluctuations are often large for localized orbitals in condensed systems (e.g., defects or moiré states) since it is problematic to sample both localized and delocalized states evenly in the real-space representation.²⁷ This translates to an increased computational cost for defects, impurities, and molecules.^26,27 Furthermore, too high fluctuations potentially hinder proper convergence and may result in sampling bias.

Here, we overcome this difficulty by constructing a hybrid deterministic–stochastic approach. We show how to efficiently decompose G and W in real-time and embed the strongly localized states. In this formalism, the problematic orbitals are treated explicitly without relying on their sampling by random vectors. A similar embedding scheme was so far employed only in static ground-state calculations.²⁸ The decomposition is general and can be used to sample arbitrary excitations. Naturally, it is especially well suited for spatially or energetically isolated electronic states. Furthermore, we employ the decomposition technique to compute Σ stochastically from an arbitrarily large subspace of interest. We show that the self-energy separation disentangles correlation contributions from different spatial regions of the system.

After deriving the formalism of the self-energy embedding and decomposition, we illustrate our methods numerically for nitrogen vacancy in large periodic cells of hBN. This system represents a realistic simulation of a prototypical single-photon emitter.^2,29–32 First, we demonstrate the reduction of the stochastic error for the defect states in the hBN monolayer. Next, we study the interaction of the defect in the hBN-graphene heterostructure.

The Green’s function (G), a time-ordered correlator of creation and annihilation field operators, describes the dynamics of an individual quasiparticle. The poles of G fully determine the single-particle excitations (as well as many other properties).

Solving directly for G is often technically challenging (or outright impossible). Alternatively, the Green’s function is often sought via a perturbative expansion of the electron–electron interactions on top of a propagator of non-interacting particles (i.e., the non-interacting Green’s function, G₀). The two quantities are related via the Dyson equation $G−1=G0−1−Σ$⁠, where Σ is a self-energy accounting, in principle, for all the many-body effects absent in G₀.

Calculations usually employ only a truncated expansion of the self-energy. Despite ongoing developments,^15,33,34 the most common approach is limited to the popular GW approximation to the self-energy, which is composed of the non-local exchange (Σ_X) and polarization (Σ_P) terms. In the time-domain, the latter operator is expressed as

where W_P is the time-ordered polarization potential due to the time-dependent induced charge density.²⁵ The potential is conveniently expressed using the reducible polarizability χ and the Coulomb kernel ν as

Evaluating the action of W_P on individual states is the practical bottleneck of the GW approach. Hence, despite that Eq. (1) together with the Dyson equation should be self-consistent, the QP energies are commonly computed by a “one-shot” correction, denoted as G₀W₀, where W₀ is computed from the non-interacting starting point.^12,17 In practice, Eq. (1) thus contains quantities constructed from the mean-field Hamiltonian, H₀, comprising one-body terms and local Hartree, ionic, and exchange-correlation potentials. For W_P, it is common to employ the random phase approximation (RPA). Beyond RPA, approaches are more expensive and, in general, do not improve the QP energies unless higher-order (vertex) terms are included in Σ.^15,35

Within this one-shot framework, the QP energies become¹² 

where ε⁰ are eigenvalues of H₀ and $vxc$ is the (approximate) exchange-correlation potential. In Eq. (3), Σ_P is in the frequency domain.

The stochastic G₀W₀ method seeks the QP energy via random sampling of wavefunctions and decomposition of operators in the real-time domain.^24–26 The expectation value of Σ is expressed as a statistical estimator. The result is subject to fluctuations that decrease with the number of samples as $1/N$⁠.

In this formalism, the polarization self-energy expression is separable.^24–26 Specifically, for a particular H₀ eigenstate ϕ, the perturbative correction becomes

where u_ζ(r, t) is an induced charge density potential and ζ is a random vector used for sampling of G₀ (discussed in detail in Sec. III and the  Appendix). The state ζ at time t is

where the U_0,t is time evolution operator, detailed in the  Appendix,

The projector P_μ(t) selects the states above or below the chemical potential, μ, depending on the sign of t. In practice, P_μ(t) is directly related to the Fermi–Dirac operator.^25,36–38 Since the Green’s function is a time-ordered quantity, the vectors in the occupied and unoccupied subspace are propagated backward or forward in time and contribute selectively to the hole and particle non-interacting Green’s functions.

The induced potential u_ζ(r, t) represents the time-ordered potential of the response to the charge addition or removal,

where $ζ¯$ spans the entire Hilbert space.^24–26 

In practice, we compute u_ζ from the retarded response potential, which is $ũζ=∫W̃P(r,r′,t)ζ¯(r′)ϕ(r′)d3r′$⁠; the time-ordering step affects only the imaginary components of the Fourier transforms of u_ζ and $ũζ$ (see the  Appendix).^13,25,26

The retarded response, $ũζ$⁠, is directly related to the time-evolved charge density $δn(r→,t)≡n(r,t)−n(r→,t=0)$ induced by a perturbing potential δ$v$⁠,

where we define a perturbing potential

Note that δ$v$ is explicitly dependent on the state ϕ and the $ζ¯$ vector; the latter is part of the stochastically decomposed G₀ operator.

The stochastic formalism further reduces the cost of evaluating $ũζ$⁠. Instead of computing $δn(r→,t)$ by a sum over single-particle states, we use another set of random vectors $η$ confined to the occupied subspace. Time-dependent density n(r, t) thus becomes^{24–26,37–39}

where η_l is propagated in time using U_0,t and H₀ in Eq. (6) is implicitly time-dependent. As common,^12,17 we resort to the density functional theory (DFT) starting point. Since H₀ is therefore a functional of $n(r→,t)$⁠, the same holds for the time evolution operator,

Furthermore, we employ RPA when computing $ũζ$⁠; this corresponds to evolution within the time-dependent Hartree approximation.^36,40,41

Practical calculations use only a limited number of random states. Consequently, the time evolved density exhibits random fluctuations at each space–time point. To resolve the response to δ$v$⁠, we use a two-step propagation whose difference is the δn that typically converges fast with N_η.^24–26 

The stochastic vectors $ζ$ and $η$⁠, introduced in Sec. II B, are constructed on a real-space grid and sample the occupied (or unoccupied) states. The number of these vectors (N_ζ and N_η) is increased so that the statistical errors are below a predefined threshold. Furthermore, the underlying assumption is that $ζ$ and $η$ sample the Hilbert space uniformly.

Here, we present a stochastic approach restricted only to a subset of states, while selected orbitals, $ϕ$⁠, are treated 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 approach for (i) the Green’s function, (ii) the induced potential, or (iii) both G₀ and u simultaneously. In the following, we present each case separately.

The non-interacting Green’s function is decomposed into two parts (omitting for brevity the space–time coordinates),

where $G0ϕ$ is the Green’s function of the constructed explicitly from $ϕ$ as

where θ is the Heaviside step function responsible for the time-ordering (corresponding to particle and hole contributions to G₀). The complementary part, $G̃0$⁠, is sampled with random states as in Eq. (4).

Unlike in the fully stochastic approach (where no $G0ϕ$ term is present), the sampling vectors are constructed as orthogonal to the {ϕ}-subspace, i.e.,

Here, $ζ¯0$ spans, in principle, the entire Hilbert space and P_ϕ is

The construction of $G̃0$ remains the same as in the fully stochastic case: $G̃0$ is decomposed by a pair of vectors $ζ¯$ and ζ(t), cf. Eqs. (5), (6), and (14).

Note that it is possible to generalize the projector P_ϕ, Eq. (15), to an arbitrary $ϕ$-subspace. The particular choice of ϕ does not affect the decomposition of the Green’s function [Eq. (12)] or the time evolution of ζ(t). However, the dynamics of $G0ϕ$ would require explicit action of U_0,t on ϕ that is, in principle, not an eigenstate of H₀. Furthermore, instead of an intrinsically localized $ϕ$-subspace, one can also employ localized orbitals constructed by a unitary transformation of orbitals.^42–44 We do not pursue this route here and select {ϕ}-subspace composed from the starting point eigenstates.

The retarded induced potential $ũζ$ is decomposed through partitioning of the time-dependent charge density, cf., Eq. (8) (omitting the space–time coordinates),

where the n_ϕ is the density constructed from occupied states ϕ (which we assume to be mutually orthogonal),

Here, f_j is the occupation of the jth state. Note that choosing ϕ within the unoccupied subspace is meaningless in this context since only occupied states contribute to the charge density n_ϕ. The complementary part, $ñ$⁠, is given by stochastic sampling [Eq. (10)] that employs random vectors $η̃$⁠,

where η spans the entire occupied subspace and

The time propagation of the charge density is similar to the fully stochastic case. The deterministic and stochastic vectors evolve as

where U_0,t is explicitly expressed as a functional of the two density contributions from Eq. (16). Note that these expressions are analogous to Eq. (11).

Both partitionings are used in conjunction. While G₀ contains contributions from both occupied and unoccupied states, only the former are included in $ũζ$⁠. The combined partitioning may use entirely different subspaces for the Green’s function and the induced potential. In Sec. V C, we employ the decomposition in both terms because it yields the best results and significantly reduces statistical fluctuations.

In Sec. III we partition retarded induced potential and the Green’s function, intending to decrease stochastic fluctuations in the self-energy. Here, we use the partitioning to achieve the second goal of this paper—to determine the contribution to Σ_P(ω).

Conceptually, we want to address quasiparticle scattering by correlations from a particular subspace. In the expression for Σ_P(ω) [Eq. (4)], this corresponds to accounting for selected charge density fluctuations in $ũζ$⁠. In practice, we construct the subspace polarization self-energy as

where we introduced the subspace induced potential $uζs$⁠, which is obtained from its retarded form [in analogy to Eq. (8)],

This potential stems from the induced charge density that includes contributions only from selected orbitals {ϕ}. Note that $nsr′,t$ is obtained either from individual single-particle states or from the stochastic sampling of the {ϕ}-subspace under consideration.

If the set of ϕ states is large, it is natural to employ the stochastic approach; the density is sampled according to Eq. (10) with vectors η^s prepared as

where the projector is in Eq. (19) and vectors η span the entire occupied subspace.

The time evolution of η^s follows Eq. (11), i.e., it is governed by the operator U_0,t that depends on the total time-dependent density,

Hence, despite $ΣPs(t)$ contains only fluctuations from a particular subspace, the calculation requires knowledge of the time evolution of both n^s and n.

In practice, we employ a set of two independent stochastic samplings: (i) vectors $η$ describing the entire occupied space and (ii) vectors $ηs$ confined only to the chosen {ϕ}-subspace. The first set characterizes the total change density fluctuation and enters U_0,t in Eq. (25). The stochastic–deterministic GW algorithm is in the  Appendix.

In this section, we will demonstrate the capabilities of the method introduced above. The starting-point calculations are performed with a real-space DFT implementation, employing regular grids, Troullier–Martins pseudopotentials,⁴⁵ and the Generalized gradient approximation (GGA)⁴⁶ functional for exchange and correlation. We investigate finite and 2D infinite systems using modified periodic boundary conditions with Coulomb interaction cutoffs.⁴⁷ 

The numerical verification for the SiH₄ molecule is in Sec. V B. To converge the occupied H₀ eigenvalues to <5 meV, we use a kinetic energy cutoff of 26 E_h and 64 × 64 × 64 real-space grid with the step of 0.3 a₀.

Large calculations for the V_N defect in hBN monolayer and in hBN heterostructure with graphene are in Secs. V C and V D. The heterostructure is built with an interlayer distance of 3.35 Å. In both cases, we consider relaxed rectangular 12 × 6 supercells containing 287 and 575 atoms. We performed the structural optimization in the QuantumESPRESSO code⁴⁸ together with Tkatchenko–Scheffler’s total energy corrections⁴⁹ and effective screening medium method.⁵⁰ Our band structure calculations in low dimensional systems have been validated against the QuantumESPRESSO calculations.²⁷ 

The GW calculations were performed using a development version of the StochasticGW code.^24–26 The calculations employ an additional set of 20 000 random vectors used in the sparse stochastic compression used for time-ordering of $ũζ$⁠.²⁵ The time propagation of the induced charge density is performed for a maximum propagation time of 50 a.u., with the time step of 0.05 a.u.

We verify the implementation of the methods presented in Secs. III and IV on the SiH₄ molecule. In particular, we test separately (i) the construction of the embedding schemes, (ii) decomposition of the time propagation, and (iii) the evaluation of the subspace self-energy.

The reference calculation employs only one level of stochastic sampling (for the Green’s function, while the rest of the calculation is deterministic). For small systems, this approach converges fast as the stochastic fluctuations are small.^25,26 Yet, we need N_ζ = 1500 vectors to decrease the QP energies errors below 0.08 eV. An illustration of the self-energy for the HOMO state of SiH₄ is in Fig. 1; everywhere in this figure, the stochastic error is below 0.13 eV for all frequencies.

We first inspect the results for an embedded deterministic subspace. Figure 1 shows that the different schemes for the explicit treatment of the HOMO state (Sec. III) produce the same self-energy curve. The inclusion of HOMO in W is combined with the stochastic sampling of the three remaining orbitals by N_η = 16 random vectors. Note that it is not economical to sample the action of W_P for small systems,^25,26 and this calculation serves only as a test case. This treatment yields a statistical error of 0.08 eV, i.e., the same as the reference calculation despite 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 the reference calculation (0.05 eV). Such a result is expected since the reference relies on the stochastic sampling of the Green’s function. The same happens when the frontier orbital is in both G and W (N_η = 16 random vectors sample the induced charge density). Tests for other states are not presented here but lead to identical conclusions.

Next, we verify that the density entering the time-evolution operator U_0,t[n(t)] can be constructed from different states than it is acting on. Specifically, the induced charge and the time-dependent densities may be sampled and built by different means. To demonstrate this, we propagate each of the H₀ eigenstates with U_0,t[n(t)], where n(t) is stochastically sampled by N_η = 16 random vectors. Only the induced charge density, entering Eq. (8), is computed from the {ϕ} eigenvectors. The agreement with the reference self-energy curve is excellent with differences smaller than the standard deviations at each frequency point [see Fig. 1(a)].

Finally, we inspect the subspace self-energy in which U_0,t[n(t)] employs the total charge density sampled by N_η = 16 random vectors. Figure 1(b) shows four different Σ^s(ω) curves corresponding to the contributions of individual H₀ eigenstates. Since the eigenstates are orthogonal, the total self-energy is simply the sum of individual Σ^s(ω) components. The additivity of the subspace self-energies is demonstrated numerically in Fig. 1(a). The subspace results illustrate that HOMO and the bottom valence orbital exhibit the largest amplitudes of $ΣPs$⁠; hence, these two states dominate the correlation near the ionization edge.

The deterministic subspace embedding should numerically stabilize the stochastic sampling and decreases the computational cost. To test the methodology on a realistic system, we consider the electronic structure of an infinitely periodic hBN monolayer containing a single nitrogen vacancy (V_N) per a unit cell with dimensions of 3.0 × 2.6 nm². The system comprises 1147 electrons with the defect state being singly occupied and hence positioned at the Fermi level.

In the current calculations, we enforce spin degeneracy. The reason is twofold: (i) half populated states are strongly polarizable and they exhibit stronger stochastic fluctuations in the time evolution. They are thus a more stringent test of the embedding. (ii) In Sec. V D, we compare the monolayer with a heterostructure to determine the role of substrate material on the self-energy. In the heterostructure, the magnetic splitting of spin-up/down components disappears.⁵¹ 

The relaxed monolayer geometry shows only mild restructuring. The vacancy introduces three localized states with the D_3h point group symmetry.^51,52 The singly occupied in-gap state is labeled by its irreducible representation $a2″$⁠.⁵² The second state (e′) is doubly degenerate and pushed high in the conduction region. Due to the enforced spin-degeneracy, the electron–electron interactions are increased and e′ appears higher than in the previous calculations.^29,53–55 The $a2″$ and e′ single-particle wavefunctions are illustrated in Fig. 2.

We first focus on the delocalized top valence and bottom conduction states. The fully stochastic calculations converge fast for both of them. The charge density fluctuations are sampled by N_η = 8 vectors, and the Green’s function requires N_ζ = 1500 to yield QP energies with statistical errors below 0.03 eV for the valence band maximum (VBM). The error for the conduction band minimum (CBM) is <0.01 eV. The computed resulting quasiparticle bandgap is 6.49 ± 0.04 eV, in excellent agreement with previous calculations and experiments (providing a range of values between 6.1 eV and 6.6 eV).^56–58 

To investigate the electronic structure in greater detail, we employ the projector-based energy–momentum analysis based on supercell band unfolding.^27,59–61 In practice, the individual wavefunctions within our simulation cell are projected onto the Brillouin zone of a single hBN unit cell. The resulting band structure is shown in Fig. 2(a). Since our calculations employ a rectangular supercells, the critical point K of the hexagonal Brillouin zone appears between the Γ and X points (marked on the horizontal axis by ⋆). The figure shows that the fundamental bandgap is indirect; the direct transition is 7.39 ± 0.04 eV, extremely close to the results for the pure hBN monolayer (reported to be in a range between 7.26 eV and 7.37 eV).^58,62,63

The defect states appear as flat bands, labeled in Fig. 2 by their irreducible representations. As expected from the outset, the stochastic calculations for $a2″$ exhibit large fluctuations. While each sample is numerically stable, random vectors produce a self-energy curve with a significant statistical error at each frequency point, i.e., the time evolution is strongly dependent on the initial choice of {ζ} and {η}. In this example, only the $a2″$ state exhibits such behavior.

Figure 2 illustrates the self-energies for the band edge and the defect states. For all the cases, the plots show the spread of the Σ_P(ω) curves: these correspond to the outer envelope for 15 distinct calculations, each employing 100 ζ sampling vectors (combined with N_η = 8 each). For the $a2″$ defect state, the stochastic sampling is possibly biased. The variation is three times as big as the spread of the VBM and almost seven times larger compared to CBM. The large spread is associated with increased statistical errors (50 meV), notably larger than those of delocalized VBM (30 meV) and CBM (9 meV) states. Away from the QP energy, the fluctuations increase even further; the spread of the samples becomes two times larger near the maximum of Σ_P at −20 eV. The convergence of the QP energy is poor, and the low sample standard deviation (roughly twice as big as for VBM) suggests incorrect statistics. In practice, each sampling yields a self-energy curve that lies outside of the standard deviation of the previous simulation.

The deterministic embedding remedies insufficient sampling without significantly impacting the computational cost. Naturally, we select the $a2″$ defect state and treat it explicitly (while randomly sampling the rest of the orbitals). Hence, according to the notation of Sec. III, the {ϕ}-subspace contains only a single orbital.

The decomposition of the non-interacting Green’s function follows Eq. (12) and stabilizes the sampling. The spread of the self-energy curves decreases approximately three times for a wide range of frequencies. With the embedded $a2″$ state, the statistical error of the QP energy decreases smoothly and uniformly with the number of samples. Each new sampling falls within the error of the 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 alone is less promising. Fundamentally, δn(r, t) contains contributions from the entire system, and a small {ϕ}-subspace will unlikely lead to drastic improvement. Indeed, the statistical errors and convergence behavior remain the same 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 decompositions use the same (or overlapping) {ϕ}-subspace, the induced charge density and the potential δ$v$ [Eq. (9)] share (at least some) ϕ states. Consequently, the sampling of Σ_P(ω) becomes less dependent on the particular choice of ζ vectors, which sample only states orthogonal to {ϕ}. Indeed, the embedding of a single localized state in G and W results in a nearly fourfold reduction of the statistical fluctuation that translates to more than an order of magnitude savings in the computational time. The error of the QP energy of the $a2″$ state is approximately half of the error for VBM (16 meV for N_ζ = 1500).

For completeness, we also applied the three types of embedding on VBM, which does not suffer from bias or large statistical errors. Unlike for the localized state, we observe only a negligible reduction in the fluctuations. Specifically, the stochastic error for VBM decreases from 30 meV to 29 meV regardless of whether G or G + W embedding is used. For delocalized states, the fully stochastic sampling is thus sufficient as expected from our previous work.^25–27 Note that the total cost of the calculation scales quadratically²⁵ with the number of states used to evaluate the retarded potential $ũζ$ [Eq. (8)]. An optimal computational cost is achieved when strongly localized states are treated via embedding, but the subspace of delocalized states is sampled stochastically.

Having an improved description of the localized states in hand, we now turn to the decomposition of the self-energy. The real-time approach described in Sec. IV allows inspecting the many-body interactions from a selected portion of the system. In contrast to Subsection V C, the subspace of interest contains a large number of states (irrespective of their degree of localization), and it is randomly sampled. The goal of this decomposition is to understand the role of correlation, especially at interfaces.

To test our method, we investigate a periodic hBN monolayer containing a single V_N defect placed on graphene. Such a heterostructure has also been realized experimentally.^64–67 The structure contains 2299 valence electrons, and it is illustrated in Fig. 3 together with selected orbital isosurfaces. The $a2″$ state is energetically lower than e′ (as in the pristine hBN monolayer). Both defects only weakly hybridize with graphene and remain localized within the hBN sublayer. However, the presence of graphene leads to an increased delocalization of the defect states within the monolayer⁶⁸ and decrease in the hBN fundamental bandgap, attributed to image-charge effects.⁶⁹ 

By unfolding the wave functions of the bilayer onto the hBN conventional rectangular unit cell, we obtain the band structure shown in Fig. 3(a). The graphene and hBN bands are distinguished by a state projection on the densities above and below the center of the interlayer region. Some states, however, extend appreciably across this boundary and lead to the appearance of spurious low-intensity peaks in the unfolded band structure.⁷⁰ Yet, individual bands are easy to distinguish, and one can clearly disentangle individual hBN and graphene states.

The graphene portion of the band structure reproduces the well-known semimetallic features with a Dirac point located at the K boundary of the hexagonal Brillouin zone. As discussed above, the K appears in between Γ and X of the rectangular cell, and it is labeled by ⋆ in Fig. 3. The typical Dirac cone dispersion of graphene is only a little affected by the hBN presence, and the Dirac point remains close to the Fermi level despite the charge transfer from the $a2″$ defect⁵¹ This is not surprising: sparse charge defects lead to only weak doping of graphene. Furthermore, the previous DFT results employed small supercells and may suffer from significant electronic “overdelocalization”^71,72 that spuriously enhances charge transfer.

The hBN part of the band structure (Fig. 3) is only weakly affected by the heterostructure formation. The delocalized states are qualitatively identical to those in the monolayer (Fig. 2). The fundamental bandgap remains indirect and reduced to 6.32 ± 0.04 eV. The screening introduced by graphene thus leads to a small change in E_g (<0.2 eV compared to the monolayer). The positions of the defect states, however, change notably. Here, both appear above the Fermi level due to the charge transfer from $a2″$ to the graphene layer.

Due to the charge transfer and Fermi level shift, it is clear that graphene is responsible for altering the charge fluctuations, i.e., the polarization part of the self-energy. In practice, graphene acts as a dielectric background inducing a significant screening of the Coulomb interaction in hBN. It stands to reason that the localized defect states would be strongly affected by such a polarizable layer and that the corresponding self-energy should be dominated by the spectral features originating in graphene.

To investigate the degree of coupling and the contribution to the self-energy from each monolayer, we compute $ΣPs(ω)$ using Eq. (22). Here, the {ϕ}-subspace is constructed from stochastic samples of the 576 occupied graphene states (distinguished by the green color in Fig. 3). For each sampling of the Green’s function, the subspace is described by eight random vectors in the {ϕ}-subspace. Additional eight vectors sample the complementary subspace.

Figure 3 shows the decomposition of the self-energy, indicating the contribution of graphene. For the delocalized states, the Σ_P(ω) curves of the heterostructure (black solid lines) appear similar to those of an hBN monolayer (blue solid lines). Specifically, we see that Σ_P(ω) for VBM [panel (b)] has practically identical frequency dependence between −15 eV and 0 eV. While in-plane contributions from hBN (blue dashed line) dominate the entire curve, screening from the graphene (green dashed line) substrate is substantial. A naïve comparison between the self-energy curves of the hBN monolayer and the heterostructure suggests that the enhanced maximum in ReΣ_P at −25 eV (marked by an asterisk) is caused extrinsically by graphene. However, this is not the case, and the effect is only indirect: the presence of graphene leads to shifting of the hBN spectral features. At the QP energy (marked by the dashed vertical line), graphene contributes to total Σ_P of the heterostructure by 36%.

The situation is different for CBM [panel (e)]. Toward the static limit (ω → 0), the self-energy curve becomes more negative. This shift is caused directly by the induced density fluctuations in the graphene substrate. The substrate screening renormalizes the quasiparticle bandgap, as has been shown in the literature of low-dimensional materials.^{9,27,69,73–75}

The self-energy curves of the $a2″$ panel (c) and e′ [panel (d)] defect states are significantly different from those computed for the monolayer. In both cases, we observe a negative shift related to the QP stabilization by non-local correlations. The Σ_P contributions to the defects’ QP energies are roughly three times larger in the heterostructure. Here, the intrinsic hBN interactions are the major driving force. The defect states are mostly affected by the in-plane induced density fluctuations (constituting 64% and 76% of the total polarization self-energy for $a2″$ and e′). Like the VBM, graphene indirectly acts on the defect QP states by enhancing the in-plane fluctuations, rather than direct coupling to the defects.

In the same vein, ImΣ(ω = ε^QP) of the graphene subspace is small, showing that there is only a small lifetime broadening due to the substrate. In detail, the graphene monolayer contributes to the broadening of the VBM state, $a2″$ state, and e′ around twice less than the hBN subspace. The only exception is the CBM state for which the substrate contribution dominates (by a factor of ∼3). As shown in this example, the stochastic subspace self-energy efficiently probes dynamical electron–electron interactions at interfaces. The random sampling allows selecting an arbitrary subspace to identify its contribution to the correlation energy and the excited state lifetimes, obtained from real and imaginary parts of Σ(ω = ε^QP), respectively.

Stochastic Green’s function approaches represent a class of low-scaling many-body methods based on stochastic decomposition and random sampling of the Hilbert space. The overall computational cost is impacted by the presence of localized states that increase the sampling errors. Here, we have introduced a practical solution for the G₀W₀ method. The localized states are treated explicitly (deterministically), while the rest is subject to stochastic sampling. We have further shown that the subspace-separation can be applied to decompose the dynamical self-energy into contributions from distinct states (or parts of the systems).

Using nitrogen vacancy in hBN monolayer, we demonstrate that the deterministic embedding dramatically reduces the statistical fluctuations. Consequently, the computational time decreases by more than an order of magnitude for a given target level of stochastic error. The embedding may be made within the Green’s function and the screened Coulomb interactions independently or simultaneously; the latter provides the best results. For delocalized states, embedding is not necessary and performs similar to the standard uniform sampling with real-space random vectors. The mechanism of embedding presented here is general, and it is the first step toward hybrid techniques that treat selected subspace at the distinct (higher) level of theory. The development of hybrid techniques employing various orbital localization procedures is currently underway.

We further demonstrate that the subspace self-energy contains (in principle, additive) contributions from the induced charge density. The self-energy contribution from a particular portion of the system is computed by confining sampling vectors into different (orthogonal) subspaces. We exemplify the capabilities of such calculations on defects state in the hBN-graphene heterostructure. Here, the charged excitation in one layer directly couples to density oscillations in the substrate monolayer. The electronic correlations for the defects are governed by the interactions in the host layer; the substrate only indirectly affects their strength.

This example serves as a stimulus for additional study of the defect states in heterostructures. Here, the coupling strength between individual subsystems is tunable by particular stacking order and induced strains. The subspace self-energy represents a direct route to explore such quantum many-body interfacial phenomena.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

This work was supported by the NSF through the Materials Research Science and Engineering Centers (MRSEC) Program of the NSF through Grant No. DMR-1720256 (Seed Program). We gratefully acknowledge support via the UC Santa Barbara NSF Quantum Foundry funded via the Q-AMASE-i program under Award No. DMR-1906325. The calculations were performed as part of the XSEDE⁷⁶ computational Project No. TG-CHE180051. Use was made of computational facilities purchased with funds from the National Science Foundation (Grant No. CNS-1725797) and administered by the Center for Scientific Computing (CSC). The CSC is supported by the California NanoSystems Institute and the Materials Research Science and Engineering Center (MRSEC; Grant No. NSF DMR-1720256) at UC Santa Barbara.

The new implementation is derived from the fully stochastic algorithm presented in Ref. 25. The steps (1), (2), and (a) are new, i.e., a result of the embedding introduced in this article. Steps (b)–(d) and (3) are similar to the fully stochastic scheme, and they employ hybrid quantities [i.e., decomposed time-dependent density Eq. (16) and/or the decomposed Green’s function Eq. (12)]. Step 4 is identical to the original algorithm.

The steps of the algorithm are as follows:

1. Select a set of deterministic states ϕ for embedding or subspace decomposition. Construct the projector P_ϕ [Eq. (15)] for the Green’s function and/or $P̃ϕ$ [Eq. (19)] for the retarded induced potential.

2. Create $Nζ¯0$ stochastic vectors (typically in the order of 10²–10³) for sampling the full Hilbert space. If including deterministic states in the Green’s function, construct $ζ¯$ via projection with P_ϕ [Eq. (14)]. For each sampling of the Green’s function, perform the next steps:

   • Create a set of $Nη̃$ (typically <10) random vectors in the occupied subspace. Use $P̃ϕ$ [Eq. (19)] and either

     • construct $η̃$ via Eq. (18) for the embedding scheme

     • or construct $η̃s$ via Eq. (24) for the subspace self-energy decomposition.

   • Calculate δ$v$(r) by Eq. (9) and perturb $η̃$ and/or ϕ as

with $λ≈10−3−10−5Eh−1$⁠.

Two propagations are needed: with λ = 0 and with a finite value of λ. We employ the random phase approximation, corresponding to the time-dependent Hartree approach. The time evolution thus depends only on the total charge density constructed using Eq. (16).

and then apply a time-ordering operation¹³ in frequency domain,

1. Compute ζ(r, t) for negative and positive times [for hole and particle components of the Green’s function—Eq. (5)]. If embedding in the Green’s function is applied, construct G^ϕ components using Eq. (13). Accumulate the matrix element of the self-energy via Eq. (4).

2. When steps (1)–(3) are performed, average the resulting ⟨ϕ|Σ(t)|ϕ⟩ from each sampling of the Hilbert space (which employed $ζ¯$ and ϕ orbitals). Perform Fourier transformation of the self-energy, and solve Eq. (3) in the frequency domain.

• (c)

  Perform time-propagation in discrete time steps dt (typically 0.05 a.u.) for states $η̃λ(r)$ and ϕ^λ(r),

• (d)

  Compute $ũζ(t)$ as the difference of perturbed and unperturbed Hartree potentials at each time step for λ = 0 and for a finite value of λ,