Dynamical configuration interaction: Quantum embedding that combines wave functions and Green’s functions

Dynamical configuration interaction: Quantum embedding that combines wave functions and Green’s functions

Marc Dvorak marc.dvorak@aalto.fi Department of Applied Physics, Aalto University School of Science, 00076-Aalto, Finland    Patrick Rinke Department of Applied Physics, Aalto University School of Science, 00076-Aalto, Finland
July 19, 2019

We present the concept, derivation, and implementation of dynamical configuration interaction, a quantum embedding theory that combines Green’s function methodology with the many-body wave function. In a strongly-correlated active space, we use full configuration interaction (CI) to describe static correlation exactly. We add energy dependent corrections to the CI Hamiltonian which, in principle, include all remaining correlation derived from the bath space surrounding the active space. Next, we replace the exact Hamiltonian in the bath with one of excitations defined over a correlated ground state. This transformation is naturally suited to the methodology of many-body Green’s functions. In this space, we use a modified /Bethe-Salpeter equation procedure to calculate excitation energies. Combined with an estimate of the ground state energy in the bath, we can efficiently compute the energy dependent corrections, which correlate the full set of orbitals, for very low computational cost. We present dimer dissociation curves for H and N in good agreement with exact results. Additionally, excited states of N and C are in excellent agreement with benchmark theory and experiment. By combining the strengths of two disciplines, we achieve a balanced description of static and dynamic correlation in a fully ab-initio, systematically improvable framework.

Accurate theoretical predictions for systems with strongly-correlated electrons remains one of the major challenges in condensed matter physics and quantum chemistry. Even though this topic has been intensely studied for decades, research continues because of its great importance to both fundamental physics and technological applications. In principle, quantum many-body theoryFetter and Walecka (1971); Altland and Simons (2010) allows one to study electronic properties at the most fundamental level. Of course, the exact theory is hopelessly complex for realistic systems. Developing less expensive, approximate methods which predict spectra of correlated systems is therefore essential for the theoretical discovery of new phenomena or new materials.

Of particular interest in condensed matter physics are systems which couple a small number of strongly-interacting electrons to a much larger bath of weakly-interacting electrons. - or -electron atoms on surfaces or point defects in solids Bockstedte et al. (2018) are possible examples of such impurities. In methods designed for impurity problems, the important impurity states are described with a theory that is much more accurate and computationally expensive than the theory treating the bath. Strongly-interacting electrons on the impurity are referred to as statically correlated, which means they are energetically near each other, while the continuum of states in the bath is dynamically correlated. The multi-reference character of open-shell molecules in quantum chemistry also requires a balanced treatment of static and dynamic correlation. Developing a theory which can treat both regimes of correlation on equal footing within one unified framework is difficult.

Such impurity problems are so extensively studied and present such a major reduction in computational cost that it can be advantageous to construct an artificial impurity from an otherwise homogeneous system. This is the central idea behind quantum embedding or active space (AS) theories an effective impurity is selected from the complete system and treated with high accuracy. If the relevant physics are determined primarily by the impurity states, then such a partitioning gives a good overall description of the system. -electron levels in solids or a chemical active space in an otherwise uniform molecule can be considered an effective impurity cut out from a homogeneous system. With a sensible and physically motivated choice of active space, the impurity concept can be applied to any system, not only those with an obvious physical defect. Quantum embedding theoriesSun and Chan (2016) are therefore powerful methods for strongly-correlated electrons.

By construction, electrons in the impurity and bath exist in different regimes of correlation: static and dynamic. Static correlation on the impurity is best described by brute force exact diagonalization of the many-body Hamiltonian. While the method describes all levels of static correlation, it is far too expensive for large systems. Correlation among high energy degrees of freedom in the bath is described efficiently by the many-body Green’s function (GF). When the long-lived quasiparticle concept applies, approximate Green’s function methods can be applied to much larger systems than wave function methods and give excellent results for comparatively low cost.

In this work, we derive and test a new quantum embedding theory that treats static correlation with the electronic wave function (WF) and dynamic correlation with the many-body Green’s function. Unlike previous embedding theories which treat a single central object (either the WF or GF) described at different levels of theory in each space, we consider different quantities in either space. The conceptual key to our approach is to selectively rewrite the many-body Hamiltonian, depending on the chosen perspective of the physics. The direct construction of the Hamiltonian in the WF picture is based on bare electrons in the vacuum. An equivalent representation is to describe excitations above a correlated ground state. By choosing different representations of the Hamiltonian in either space, we can choose to work with either the WF or GF.

I Introduction

i.1 Brief survey of current methods

We very briefly introduce many-body WF and GF methods that are relevant to the current work. Many modern approaches to strong correlation in physics study the many-body Green’s function.Nolting (1971); Fetter and Walecka (1971) We refer to such applications of quantum field theory in condensed matter as many-body perturbation theory (MBPT), which is distinct from the chemists’ MBPT based on perturbation in the residual potential (as in Møller-Plesset perturbation theory or related methods). The GF is a direct link to spectral information about the system. Particle addition/removal spectra, as well as the neutral excitation spectrum, are accessible by solving for the single-particle or two-particle . In contrast, quantum chemistry focuses on efficient approximations to the many-body wave function. All of the same spectral information as with GFs is accessible if one knows the many-body wave functions for all eigenstates of the system.

The approximation Hedin (1965); Aryasetiawan and Gunnarsson (1998); Hybertsen and Louie (1986) and its extension to optical properties, the Bethe-Salpeter equation Salpeter and Bethe (1951); Rohlfing and Louie (2000); Blase et al. (2018); Albrecht et al. (1998) (BSE), are very successful GF theories for predicting spectra of weakly- to moderately-correlated systems. and the common implementation of the BSE have their limitations, however. BSE with a static kernel cannot describe multiple excitations.Sangalli et al. (2011) BSE also has a tendency to underestimate triplet excitation energies and suffers from self-screening error.Bruneval et al. (2015) There is not yet a widely adopted route to improve these shortcomings of the BSE, though developing extensions to both and BSE is a very active area of research.Romaniello et al. (2009a); Sangalli et al. (2011); Attaccalite et al. (2011); Zhang et al. (2013); Lischner et al. (2012); Maggio and Kresse (2017); Starke and Kresse (2012); Hung et al. (2016); Pavlyukh et al. (2016); Isobe et al. (2018); Kuwahara et al. (2016); Grüneis et al. (2014); Romaniello et al. (2009b); Peng and Kowalski (2018); Lange and Berkelbach (2018)

Dynamical mean-field theory (DMFT) is a quantum embedding theory based on the electronic GF. In combination with the local density approximation (LDA), LDA+DMFT Georges and Kotliar (1992); Held et al. (2008); Kotliar et al. (2006) has been very successful in predicting spectral properties of strongly-correlated solids. By including cluster extensions to DMFT,Kotliar et al. (2001); Maier et al. (2005) one can treat even non-local correlation with high accuracy. However, LDA+DMFT is based on a model Hamiltonian, rather than a true ab-initio Hamiltonian, and does not recover the exact Hamiltonian as the embedded region increases in size. The +DMFT methodBiermann et al. (2003); Tomczak et al. (2012); Biermann (2014); Boehnke et al. (2016) solves many of these issues. Other Green’s function embedding methods,Chibani et al. (2016); Aryasetiawan et al. (2009) including self-energy embedding theory Zgid and Gull (2017); Lan et al. (2015); Kananenka et al. (2015); Nguyen Lan et al. (2016) (SEET), also improve upon the shortcomings of LDA+DMFT. While they are a great success overall, all GF embedding theories have a complicated frequency structure. They require the memory consuming storage of frequency dependent quantities (the GF and self-energy) and the evaluation of difficult frequency integrals.

In contrast with such GF techniques, quantum chemistry methods based on calculating the electronic wave function are free of frequency dependence but suffer from a combinatorial explosion in the basis. Szabo and Ostlund (1996); Helgaker et al. (2014) High level coupled cluster Loos et al. (2018) or configuration interaction (CI), including their multi-reference variants,Evangelista (2011); Szalay et al. (2012) are very accurate but cannot be applied to systems larger than a few electrons. A family of multi-configuration self-consistent-field methods (MC-SCF), including the complete active space self-consistent-field method (CASSCF),Olsen (2011); Li Manni et al. (2016) simultaneously optimize single-particle orbitals to a multi-configurational wave function at the same time as diagonalizing the many-body Hamiltonian. CASSCF optimizations are state dependent and include dynamic correlation beyond the orbital active space, which itself is limited to orbitals.

Perturbation theory applied to either a single- or multi-reference wave function can also give excellent results for strongly-correlated molecules. Methods based on the Rayleigh-Schrödinger (RS) variant of perturbation theory avoid the explicit energy dependence of the perturbation expansion, a feature common to GF methods or Brillouin-Wigner techniques, but suffer from their own issues like intruder states.Choe et al. (2001) The complete active space perturbation theory to second order (CASPT2) method (and similar methods based on restricted active spacesMa et al. (2016); Manni et al. (2011); Li Manni et al. (2013)) has enjoyed widespread success as a reasonable compromise between computational cost and accuracy. CASPT2 begins with orbital optimization simultaneous with full configuration interaction in an AS, followed by perturbative corrections to second order. Confusingly, these methods are commonly referred to as MBPT in the quantum chemistry community but are different formalisms than quantum field theory, which is commonly called MBPT in condensed matter physics.

i.2 Motivation

Our goal is to construct a quantum embedding theory which combines the best features of quantum chemistry, GF embedding, and /BSE theory. To recover the exact Hamiltonian as the embedded region increases in size, we use exact diagonalization (ED) or CI for the electronic wave function in the embedded region, which we denote . Both ED and CI are based on the same, exact Hamiltonian, with their only difference being the choice of basis. With this choice, the accuracy of the theory can be systematically improved by increasing the size of . In practice, a large embedded region is numerically intractable, but this limit is a useful theoretical consideration to test the theory. We must also choose a theory in which to embed the wave function calculation. Motivated by the success of GF theories for weakly-correlated systems, we describe surrounding degrees of freedom with the many-body GF. Because the WF and GF are fundamentally different quantities, the embedding framework must serve as a bridge between two disciplines that are essentially disconnected.

At a high, conceptual level, such an embedding scheme is shown in Fig. 1. By changing the size of the embedded region, one can interpolate between a normal CI calculation () or a standard calculation with MBPT (). In these two limits, we enforce that the embedding must match these two theories. This requirement is the major guiding principle behind our approach. At intermediate sizes of , a CI problem is embedded inside of MBPT. Ideally, the final result should be invariant to the choice of embedded region. This is extremely difficult in practice but is another useful consideration for constructing the theory.

Figure 1: A schematic showing the proposed embedding theory as an interpolation between two exact theories, CI and MBPT. At intermediate sizes of , a configuration interaction problem is embedded inside of MBPT.

The choice of embedding CI inside of MBPT is motivated by their different strengths and weaknesses. High energy degrees of freedom in , the complementary space to , are dominated by dynamic correlation which is included even with simple approximations in MBPT. Such algorithms, including those based on the approximation, include some dynamic correlation with polynomial scaling and are size consistent. These methods include correlation among the full set of orbitals through sums over intermediate states in the diagrammatic expansion. CI, however, very efficiently captures static correlation among degenerate configurations or low energy degrees of freedom with a single matrix diagonalization. The CI Hamiltonian is also frequency independent, making it conceptually and computationally simpler than frequency dependent quantities in MBPT. Such correlation is also accessible with Green’s function methods, in principle, but requires a challenging self-consistent solution with difficult frequency integrals. A sensible partitioning of the Hamiltonian with an embedding theory could capture most static correlation in a system with a modest sized CI calculation and avoid the difficulties of describing static correlation with the GF. Motivated by these considerations, our goal is to treat static correlation with CI and dynamic correlation with MBPT.

To clearly distinguish our theory from previous work, we point out that we do not use CI as a high accuracy calculation of the self-energy or vertex in .Zgid et al. (2012); Pavlyukh and Hübner (2007) Our goal is to keep the electronic WF in and never compute a vertex in the strongly-correlated subspace. Treating an embedded vertex in would introduce the issues we are trying to avoid storage of the GF and vertex, frequency dependent impurity solvers, etc. Our theory is also not an application of perturbation theory based on the residual potential for mean-field potential .

Ii Theory

ii.1 Subspace partitioning

In this work, the relevant Hamiltonian is the non-relativistic, Born-Oppenheimer electronic Hamiltonian,


for one-body matrix elements and Coulomb matrix elements for two-body interaction . The Hilbert space for this Hamiltonian is all Slater determinants generated from orbitals of the non-interacting Hamiltonian.

Partitioning the Hilbert space is intrinsic to any embedding theory. Because Eq. 1 is defined in a many-body Hilbert space, we first partition the many-body space and then make a connection to the single-particle picture. We define projection operators and that project onto the strongly- and weakly-correlated portions of the many-body Hilbert space.


We consider only configurations and with fixed particle number . We define configurations in as the low-energy excitations of the system. These are most easily defined in the single-particle picture with an orbital active space. The AS can be constructed with an energy cutoff around the Fermi energy (above and below) or based on chemical intuition of the problem. An excitation that promotes any number of AS occupied orbitals to AS virtual orbitals is considered low-energy, and its corresponding -particle configuration is in . Note that this definition includes all excitation levels, , for configurations with single, double, triple, and all higher excitations in . In quantum chemistry terms, is generated as in complete active space (CAS) theories. and the orbital AS are related but not identical: is a space of many-body configurations, while the AS is a collection of single-particle states.

All other configurations are placed in . Our separation implies that transitions which mix AS orbitals with orbitals outside the AS are placed in . Different projection techniques to define the AS and are a topic for further study. Defining based solely on low-energy transitions is not mandatory. Here, we focus on the formalism instead of details of the projection method.

Figure 2: Orbitals that fall within an energy cutoff around are placed in an orbital AS, shown here in magenta. Configurations with AS excitations, shown in blue, belong to the many-body subspace . Configurations that contain any other transition belong to , shown in red. We place configurations of the mixed transition type from AS to outside AS in . Excitations include all excitation levels as in complete active space (CAS) theories of quantum chemistry.

There are two different classifications shown in Fig. 2: one defining the single-particle AS and one based on configurations. In this partitioning, one cannot say if any given orbital belongs to or as one might do when embedding the single-particle GF. In this projection scheme, can be considered a high energy bath, though the bath is a set of many-particle configurations or transitions, not individual orbitals.

We use the Feshbach-Schur decomposition,Dzuba et al. (1996, 2017); Pavlyukh et al. (2015) shown in Fig. 3, to transform the exact Hamiltonian to an effective, downfolded Hamiltonian in the space. This decomposition is also known as the Löwdin projection.Manni et al. (2011); Li Manni et al. (2013); Löwdin (1962, 1968) The Schrödinger equation can be written in block form based on the two projectors


Here, is the Hamiltonian projected only into the space. This block of the Hamiltonian corresponds to the frozen core approximation in configuration interaction, where excitations are completely ignored. The component of any given eigenstate of the full in the space, , depends on the solution in through the offdiagonal blocks of the Hamiltonian. The energy is the total electronic energy of the system. The second line of the equation reads


We can solve this equation for the solution, , and insert it into the first line of the Schrödinger equation.


We have introduced the resolvent defined by Eq. 5, . We insert Eq. 5 for into the first line of the Hamiltonian,


The downfolding procedure therefore gives an effective Hamiltonian of the size . Eq. 6 is the effective Schrödinger equation in the space and is an exact rewriting of the eigenvalue problem. While the size of the matrix to diagonalize is much smaller than the original , is now energy dependent, and the energy eigenvalue must be found self-consistently. The exact solution (equivalent to full CI or ED) requires inverting the block of the Hamiltonian to compute . This is an extremely large space for any realistic problem and the inversion is numerically intractable. The remainder of our theory is to find a suitable approximation for to simplify the inversion.

In matrix notation, is constructed with the matrix elements


The indices and refer to configurations in , while and denote configurations. Matrix elements of are evaluated with the Slater-Condon rules, Slater (1929); Condon (1930); Szabo and Ostlund (1996) which are briefly presented in Appendix A. Based on the Slater-Condon rules, there are certain selection rules for the matrix elements and which depend on the differences in occupation numbers between states and . The matrix elements are the CI matrix in the frozen core approximation. The elements , therefore, contain all the correlation beyond the frozen core approximation and can be considered a dynamical core correction or a configuration self-energy.

Figure 3: Schematic showing the block structure of the Hamiltonian and the matrix multiplication involved in the downfolding. The full Hamiltonian (a) is partitioned into spaces and (b). The effect of the Löwdin downfolding is to add an energy dependent correction to (c). The resulting effective Hamiltonian is of size . In practice, can be several orders of magnitude larger than .

ii.2 Combining Green’s function- and wave function-based theories

The most important element of our theory is the combination of different methods based on wave functions and Green’s functions. To demonstrate how we connect these two theories, we discuss the separation of any total energy into a ground state and excitation energy. While this is trivial for eigenstates of the full Hamiltonian, we discuss it in detail since subtle aspects of this separation for a subspace of the full Hamiltonian are important for the embedding procedure.

First, consider exactly diagonalizing the full Hamiltonian in Eq. 1. For total energy eigenvalues and excitation energies , the Hamiltonian can be written in the eigenbasis as


We denote the excitation matrix with eigenvalues as simply . The lowest element of is for two reasons: the reference energy to define the excitation energies is the correlated ground state, and the correlated ground state is the lowest eigenvalue of .

All of the quantities in Eq. 8 are accessible with either WF- or GF-based methods. In WF-based methods, the eigenvalues of the Schrödinger equation are total energies, and excitation energies are computed as total energy differences. In contrast with such quantum chemistry methods, MBPT computes excitation energies directly as the response of the system to perturbation, e.g., particle addition/removal or optical excitations. In these GF theories, the degrees of freedom are quasiparticles electrons and holes propagating over the correlated ground state. The quasiparticle Hamiltonians of many-body Green’s functions are effective, particle-like equations of motion with eigenvalues related to excitation energies of the system.

We use the excitation energies as the connection between these two theories. When introducing Eq. 8, we assumed an excitation matrix computed as total energy differences so that the excitation matrix is . Instead, consider computing the matrix with GF. If we compute the excitation matrix with GF, we effectively replace bare electronic degrees of freedom with quasiparticles. This transformation is the quasiparticle renormalization of the many-body Hamiltonian,


where is now a quasiparticle (QP) excitation matrix. Note that the scalar energy on the diagonal, , is the same in Eqs. 8 and 9. By definition of the correlation functions used in MBPT, excitation energies are also defined with respect to the correlated ground state. The renormalization in Eq. 9 is set up for the language and methodology of many-body GF.

Most importantly, the effective quasiparticle Hamiltonian of is different than the “bare” excitation matrix defined by . Quasiparticle Hamiltonians contain sums over occupied and virtual states at intermediate times, the so-called diagrammatic expansion. For weak coupling (), these expansions involve products of Coulomb matrix elements that effectively weaken the strength of the interaction. Certain classes of diagrams can even be summed to infinite order. The end result is that matrix elements of the quasiparticle Hamiltonian are often weaker than those based on the bare Hamiltonian. Computing with a quasiparticle Hamiltonian could, therefore, converge faster with respect to the many-body basis or perform better in a diagonal approximation than the bare Hamiltonian.

MBPT does not give exactly the same eigenstates as wave function methods based on diagonalizing Eq. 1. Even exact eigenstates of MBPT have a finite lifetime derived from the decay of the bare excitation into many different states. In practice, however, quasiparticle excitations in many systems are rather long-lived. As long as excitations from MBPT have sufficiently long lifetime, it is a safe approximation to replace with . This is a point we discuss in detail in later sections. Here, we point out that numerous benchmark studies comparing the two methods suggest this replacement is a reasonable approximation in weakly- to moderately-correlated systems.Jacquemin et al. (2015); Bruneval et al. (2015); Tiago and Chelikowsky (2006); Körbel et al. (2014)

Now consider diagonalizing only the block . This is the block of the Hamiltonian needed for the resolvent, . The diagonalized is a different matrix than the block of the diagonal matrix . We want to write the projected Hamiltonian in a way similar to Eq. 8. Working in the eigenbasis, we are free to subtract some scalar energy from the diagonal,


where . Here, is the eigenvalue of from exact diagonalization, and is an as-of-yet undefined energy.

Next, we impose physical constraints on the energy and eigenvalues . A fundamental requirement of the embedding is to recover a normal MBPT calculation in the limit that . For this reason, we assign meaning as a ground state energy in and relate to subspace excitation energies. While eigenstates of are not physical excitations since they exist in only a subspace of the full , they can be connected to a physical excitation by enlarging . This is the defining criterion for the matrix . In the limit , must reach the physical excitation matrix . Similarly, the ground state must reach the physical, correlated ground state .

We point out that the lowest eigenvalue of , the matrix in Eq. 10, is not zero. Unlike in Eq. 8, even the lowest eigenvalue of is itself an excited state or, in other words, can be connected to a physical excitation. Even though no eigenstate can be connected to the physical ground state, Eq. 10 is still exactly true for

In Eq. 10, we now assume that the subspace excitation matrix is computed with MBPT, . Assume an exact diagonalization of to find . Then we redefine the ground state energy as the difference between the eigenvalue and the QP excitation energy .


Equation 11 is useful to define the problem, but computing this way requires the exact diagonalization of , which is extremely expensive. To avoid this expense, we must develop a different strategy to compute . From now on, we assume is computed with MBPT and drop the label. If we can compute the ground state energy and excitation matrix separately, and for less expense than exact diagonalization, we can assemble them to rewrite the Hamiltonian.

The goal of our theory is to perform the same quasiparticle renormalization as in Eq. 9 for only the weakly-correlated subspace of the full Hilbert space.


The high energy space is dominated by dynamic correlation, which is described very well by MBPT. In , we describe static correlation with the bare Hamiltonian; in , we treat each configuration as a quasiparticle excitation propagating above a ground state. We construct so that it matches the exact as closely as possible. is a reference energy to define the subspace QP excitation matrix so that both quantities connect to their physical values as .

After computing the renormalized Hamiltonian , we insert it into the denominator of the resolvent in place of in Eq. 6. We treat as “exact” so that Eq. 6 is not solved perturbatively. We do not take the route of Brillouin-Wigner (BW) perturbation theory in the residual operator for mean-field potential . Such an order-by-order construction of the resolvent is possible but introduces a difficult energy dependence at every term in the expansion. Nor do we apply the Rayleigh-Schrödinger (RS) variant of the perturbation expansion for . Our procedure to compute excitation energies has its own energy dependence as in BW theory, to be discussed in the next section and in Appendix C, but it is a different energy than entering . We apply perturbation theory to the subspace problem of diagonalizing , which has a separate energy dependence than . If we can diagonalize even by MBPT the connection between and is automatically set up through the resolvent.

While our framework has similarities to other methods based on the Löwdin partitioning, this is because there is just one way to exactly partition the many-body Hilbert space. Our concept and final theory are different than past work.Manni et al. (2011); Li Manni et al. (2013); Dzuba et al. (2017, 1996) To transform the Hamiltonian, we need the excitation matrix and ground state energy .

ii.3 Excitation matrix

To compute the excitation matrix with GF, there are three conceptual hurdles to overcome. First, we must match the basis for the excitation matrix to the CI basis. The CI Hamiltonian is naturally written up to all excitation levels, while the typical correlation function for neutral excitations is written in a basis of only single excitations. The two basis sets must match for the matrix multiplication of the resolvent to be meaningful.

Second, we must eliminate the frequency dependence of the GF. In general, any many-body GF and its equation of motion are frequency dependent. Using an exact GF equation of motion for would therefore give an energy dependent matrix. However, we want to avoid complicated frequency dependencies, so we must eliminate the frequency dependence in .

Last, we must constrain the calculation of the excitation matrix to the subspace. If one considers exactly diagonalizing , it is clear that this matrix contains only intra- correlation. The corresponding calculation of the excitation matrix with GF must also include only intra- correlation.

ii.3.1 Definition based on Green’s function

Here, for an easier discussion, we briefly abandon our subspace partitioning to discuss the formalism for the full Hamiltonian. The continuation to a subspace of the Hamiltonian will be discussed later.

We first review the standard theory in MBPT so that our method can be built as an extension to it. We will restrict ourselves to optical excitations in this section for illustrative purposes. The standard approach for computing optical excitations with Green’s functions is based on the electron-hole correlation function,Onida et al. (2002); Fetter and Walecka (1971) , defined as

describes the propagating portion of the two-particle GF, defined by the second line of Eq. II.3.1, and ignores the motion of independent pairs, given by the first line of Eq. II.3.1. Here, is the time ordering operator and is the -particle ground state. Nolting (1971) Each number in Eq. II.3.1 is a set of spatial, spin, and time coordinates, . In order to approximate exact excitation energies from quantum chemistry, we restrict possible time orderings to simultaneous creation/annihilation of one e-h pair () and instantaneous annihilation/creation of a second e-h pair (). can be expanded in a basis of noninteracting electron-hole pairs, also called single excitations, and its equation of motion is determined by the Bethe-Salpeter equation (BSE).Salpeter and Bethe (1951) All excitations are accessible, in principle, by solving the frequency dependent BSE.

While single excitations are clearly described by the correlation function in Eq. II.3.1, excitation energies of exact eigenstates which involve the creation of double or higher noninteracting excitations are hidden in the BSE. These so-called multiple excitations are contained in the vertex function, , of MBPT. This is most easily demonstrated in the Lehmann representation, in which is written as


where creates the single excitation and the sum runs over exact eigenstates of the many-body system. The sum over is independent of the outer indices any single matrix element of , , contains poles at every, exact excitation energy of the system. In principle, one needs only the single, frequency dependent matrix element to access all excitations. The excitation energies can be read off as the pole positions of . In practice, however, this is extremely difficult. The amplitudes in the numerator of , which are the amplitudes for each pole, are extremely small for most eigenstates . The overlap between and is appreciable only for the eigenstate which has the same principal configuration as , and perhaps a small number of additional . This makes it impossible to numerically produce most poles in at the remaining eigenstates with weak amplitude and extract any meaning as they relate to a principal noninteracting transition.

The full e-h correlation function, , has several poles with high amplitudes, with each state having high overlap with roughly one . Compared to the single element , this makes it much easier to find excitation energies and relate poles in to a noninteracting transition. Any exact eigenstate of strong single excitation character should have a large Lehmann amplitude for a corresponding noninteracting single transition. However, the correlation function is still in a basis of only single excitations, and excitation energies for states with strong multiple excitation character can still be difficult to compute. For example, the overlap between any with of strong double excitation character will be low, making the pole corresponding to difficult to produce.

If we need excitation energies for multiple excitations, we are free to choose any correlation function which may have stronger amplitudes at the relevant poles. While it is not formally necessary to use a different correlation function to compute multiple excitations, it may make the calculation much easier. Furthermore, we are not interested in the spectrum of the correlation function, only the pole positions. Only the excitation energy is needed for the excitation matrix . This grants us freedom in choosing which correlation function to compute. Just as changing from a single matrix element to the full makes it easier to find single excitations, we introduce a new correlation function in the -particle space to make it easier to compute multiple excitations. For arbitrary excitation level , we write any excitation in the Hilbert space as a string of creation and destruction operators acting on the reference configuration. We write configurations as


for reference configuration and single-particle creation (destruction) operators (). Higher excitation levels follow accordingly. Define the string of excitation operators to create configuration , which can take any excitation level, as .


The set of all generates the entire -particle Hilbert space by acting on the reference configuration. The time dependence of is inherited from the time dependence of the Heisenberg operators . We take all creation and annihilation operators in the string for to act at the same time. In this way, instantaneously creates the many-body configuration , and instantaneously annihilates it. We use without the operator hat to denote the excitation energy for the creation process .

We consider a new correlation function , defined in the representation and related to an -particle GF,


for ground state . If we expand in a basis of all possible and , the matrix for each time ordering covers the entire -particle Hilbert space. With this definition, the outer lines of can be any possible excitation, not only singles. In principle, no excitation is also an allowed state, , in order to complete the Hilbert space. In Eq. 17, is our generic notation for removing the non-propagating portion of the GF, , from the full -particle GF, . This leaves only the time evolving portion or -particle propagator, which we label . We do not make any connection between and experimental spectra, but only use it for an easier calculation of multiple excitations.

Now we return to the embedding problem and consider , which tracks the propagation of only excitations. Allowing the outer lines of to be multiple excitations is more meaningful than just a difference of convention it allows us to directly compute multiple excitation energies with a frequency independent kernel. The connection between and is through matrix elements of the exact Hamiltonian, . These matrix elements are evaluated using the Slater-Condon rules (Appendix A) for all excitation levels. In order to multiply by the resolvent , they must share a common basis. If we estimate the denominator of the resolvent with MBPT, we must use MBPT to compute multiple excitations in an efficient way. Because is of dimension , every excitation is accessible with a static kernel. This way, the excitation matrix can be calculated in a static framework, and the matrix multiplication between and is meaningful. This connection is essential for the embedding. If we base our MBPT calculation on , which is in the basis of only single excitations, instead of , the MBPT basis does not match the CI basis. In this case, the energy dependence of and is necessary to couple to multiple excitations and greatly complicates the embedding.

Returning to the general , we define the excitation matrix as the frequency space equation for based on Eq. 17. can be rewritten in the Lehmann representation by inserting a complete set of eigenstates and Fourier transforming as


where the sum runs over all -particle eigenstates of the exact Hamiltonian. has poles at the exact excitation energies of the system. This representation is not very useful, however, since one needs the exact ground state and all -particle eigenstates of the system.

In order to actually compute , we want to set up a perturbation expansion for with MBPT. We define in terms of independently propagating, noninteracting electrons and holes. The time dependence of the individual fermionic operators depends on only the one-body portion of . The individual phases associated with each and entering combine so that the entire excitation propagates with an energy given by


for excitation level and eigenvalues of the one-body Hamiltonian . The corresponding representation is


The expansion amplitudes for are all 1 since the non-interacting problem is diagonal in this basis. The one-body eigenvalues are exactly the particle addition/removal energies of the noninteracting system.

Next, we assume some frequency dependent kernel which connects the bare to the full as


where the kernel is the irreducible version of the full kernel . Diagramatically, we can represent blocks of as in Fig. 4. For each block of Fig. 4, the number of external lines is determined by the excitation levels, and , of the basis states. Each block of is meant to be taken to any desired order in the time evolution operator .

Figure 4: Schematic representation of . In general, both and the kernel are functions of frequency. The reference configuration is also included in the basis for , though we omit it here for comparison to the subspace which contains only excitations. All excitations in couple to all others so that the matrix is not sparse. Each block contains both time orderings of . The electron-hole correlation function is the block outlined in red.

and are complicated objects. In this work, our interest in them is mostly utilitarian: to efficiently compute multiple excitations and match the basis for CI. Rather than exploring the exact structure of the kernel that connects initial and final states in Fig. 4, we focus on physical approximations to . We explore the relationships among , , the vertex function , and the BSE in future work.

ii.3.2 MBPT for

We can now discuss a strategy for calculating , the frequency space equation for the subspace correlation function . The subspace correlation function of interest is


Here, is a fictitious -particle ground state that contains only intra- correlation the same ground state discussed previously. The precise meaning of this state is difficult to define since we assume the reference configuration always belongs to . However, we can still construct an approximate with physically motivated approximations and by enforcing the correct limits on and as or .

The basic idea is to apply many-body perturbation theory for in the full many-body basis from the exactly projected Hamiltonian


For this subspace Hamiltonian, each interaction line carries projectors that check the overall -particle configuration. For each new interaction in the perturbation expansion, the configuration must belong to , otherwise the entire diagram is killed by the projector. For an interaction at time , we take the two projectors to act at times vanishingly close to , for . The projectors, therefore, check the configuration of the system just before and just after each interaction. Applying the projectors this way requires us to mix the two conceptual pictures: acts on the -particle configuration, but the diagrammatic expansion for only shows GF lines for the excited particles. We only need to check the excitations above the reference configuration to know if a certain diagram/configuration belongs to , so the two pictures agree.

To order in the interaction, the excitation propagates indefinitely. This is a terrible approximation since it neglects all electron-electron interactions. Diagrammatically, each particle participating in the excitation is described by a noninteracting single-particle GF line, . The sum of lines for the excitation gives the bare , as shown in Fig. 5a. is a diagonal matrix, as in Eq. 20.

To improve this estimate, we allow each propagating particle to gain a self-energy describing its interaction with other electrons. At this level, we assume that is still separable between different particles. Each excited particle is now described by a dressed . Diagrammatically, lines do not connect to each other. Off-diagonal elements of the self-energy allow any given electron or hole to decay into any other within a block for a fixed excitation level. However, blocks of describing different excitation levels are still uncoupled so that the kernel is block diagonal. At this level, is approximated by separate dressed lines ( electrons and holes) for excitation level .

Figure 5: Conceptual hierarchy of diagrams approaching the true for single (left) and double (right) excitations. Here, outer lines can only belong to excitations. Non-interacting particles (a) contributing to the excitation are given a self-energy (b), then allowed to interact with each other (c). The dashed line separate electrons from holes for the shown time ordering.

Finally, we allow propagating quasiparticles to interact with each other. For diagrams describing their interaction, the overall configuration must at all times belong to . In principle, all possible time orderings and interactions among quasiparticles are allowed this generates the exact kernel . We relate the full to the bare by some kernel as


in analogy with Eq. 21. Here, internal frequency integrals are implied. As grows to contain the entire Hilbert space, , the subspace kernel must become the full kernel . Here, contains both levels of improvement over discussed above the dressing of bare particles by a self-energy and the interactions among all dressed lines. At this level of theory, an initial excitation of level is allowed to decay through all possible decay channels to excitation level , as in Fig. 4. The full matrix is not sparse, though elements which couple different excitation levels are expected to be small. The essential point of this subsection is that MBPT for must begin from Eq. 23.

After Fourier transforming in Eq. 22, it gains a frequency dependence, . In analogy with other subspace quantities, the frequency is not related to a physical time or energy. is the frequency variable for the Fourier transform of Eq. 22 and defined by the auxiliary eigenvalue problem of diagonalizing : for the continuous variable .

The meaning of is more obvious by considering the eigenvalues of , as well as the meaning of the ground state defining the correlation function . In the Lehmann representation for , only energy differences of the type enter the denominator. Therefore, the energy dependence of the perturbation expansion for is not the physical or in the denominator of . More discussion along this line is in Appendix C. In principle, the time variables in Eq. 22 also need to be reinterpreted. In practice, however, can be treated as a generic frequency parameter.

ii.3.3 Approximations and practical implementation for

Now, we focus on the approximations and practical aspects of computing the excitation matrix .

We entirely skip the approximation and dress each with a self-energy. Here, we adopt the approximation. In order to constrain the correlation to the subspace, we restrict the polarization to the constrained random phase approximation (cRPA) instead of the full RPA. In the cRPA, only transitions enter the polarization and screen the Coulomb interaction. The cRPA is already well known in strongly-correlated physics for determining effective Hubbard parameters in GF embedding theories.Aryasetiawan et al. (2004); Şaşıoğlu et al. (2011); Vaugier et al. (2012); Hirayama et al. (2017); Shinaoka et al. (2015) The resulting interaction interpolates between the fully screened and bare interaction as changes size from to . Accordingly, the self-energy interpolates between the full self-energy and bare exchange in these two limits.

Note that the internal line in contains all poles. Because and exist at the same times, it is sufficient to constrain only to satisfy the projectors. Essentially, only one line in the self-energy diagram must be outside the AS for the overall configuration to belong to . We assign this constraint to the polarization entering the interaction. This construction also depends sensitively on starting from a true , which has poles unambiguously assigned to single-particle states, instead of a mean-field . The projectors act on non-interacting states, not mean-field particles, and cannot be applied exactly to a starting point.

Figure 6: Schematic representation of . Each particle participating in a single (a), double (b), or higher excitation gains a self-energy based on the series. For a single excitation, interactions at odd numbered times contribute to the hole propagator, while even numbered times contribute to the electron. The generalization to double and higher excitations follows by alternately assigning interactions at different times to different particles. The diagram demonstrates the complexity of applying the projectors at every intermediate time. To avoid a complex procedure, we simply apply the cRPA at the single quasiparticle level to constrain correlation to .

For , there are many particles propagating simultaneously. We assign a constrained self-energy to each propagating particle, with each insertion taken at different times as in Fig. 6. One can roughly argue that the self-energy insertion for any one particle should be invariant to all possible time orderings of the overall expansion. To satisfy this invariance, it follows that the self-energy of each particle should be constrained in the same way. This is not a rigorous proof, but a detailed exposition of all possible time orderings, excitation levels, and diagrammatic insertions for is tedious.

Even though the application of at each interaction appears to be formally well-defined, the mixing of the many-body projector with a perturbation expansion based on the single-particle is extremely complicated. To circumvent the complexity of applying exactly, we take advantage of a known and successful result (the cRPA) and apply it to our perturbation expansion at the single quasiparticle level. Crucially, applying the approximation to each particle correctly interpolates between our two required limits. As , the self-energy correctly reaches the full self-energy. As , the quasiparticle energies correctly approach the HF eigenvalues. The special limit of is discussed in detail in Appendix B.

We denote the GF described at this level of approximation by . We assume that lines are dominated by well-defined quasiparticles with a long lifetime. By either making a diagonal approximation or working in the diagonalized-QP basis, the excitation matrix at this level of theory, , has only diagonal elements . Neglecting the imaginary part of the self-energy, the excitation energy is the sum of these noninteracting quasiparticle energies,


The Fourier transform of is


The poles of are determined by the excitation energies .

Next, we must include interactions among quasiparticles. Quasiparticle interactions between the electron and hole are known to be very important, for example, in optical excitations. We also expect them to be important for the multiple excitations needed here. We assume a Dyson equation relating to the full as


Eq. 27 can be formally solved as


By Fourier transforming the exact in Eq. 22, the poles of are determined by the energy differences . Equating these energy differences with the pole positions of in Eq. 28 sets up an effective -particle Hamiltonian Martin et al. (2016)


for subspace excitation energies .

We can now approximate the kernel . For a single excitation in , we can use exactly the result for the electron-hole kernel based on our approximation. The electron-hole kernel is based on the derivative of the self-energy. In the common approximation to the BSE, the electron and hole therefore attract via and repel via their exchange interaction (ignoring variation of with respect to ). For multiple excitations, we include pairwise electron-hole interactions using the same interaction between all possible pairs.

The multiple excitations in also include electron-electron and hole-hole interactions. For this, we extrapolate the result for electron-electron terms known for trion matrix elements.Deilmann et al. (2016); Deilmann and Thygesen (2017) Their interaction is similar to the electron-hole case, except that both their direct and exchange interactions are screened. We use screened direct and screened exchange interactions among all electron-electron pairs and hole-hole pairs.

Figure 7: Schematic representation of the approximate . At this level of approximation, the exact kernel is approximated by separate, two-body interactions for single (a), double (b), and higher excitations in . The double-wiggly line represents the partially screened Coulomb interaction, , while corresponding exchange diagrams are implied but not shown.

The diagram in Fig. 7 represents this approximation diagramatically. We approximate with separable, effective two-body interactions. We take each quasiparticle interaction to occur at different times so that can be approximated in this way.

Last, we make a global static approximation to the frequency dependent problem shown in Fig. 7. We evaluate lines at their corresponding quasiparticle energies, and take each quasiparticle interaction to be instantaneous. These approximations transform the frequency dependent problem into a static renormalized Hamiltonian. The same set of approximations is very commonly applied to the BSE and quite successful for predicting optical excitation energies of weak- to moderately-correlated systems.Rohlfing and Louie (2000); Bruneval et al. (2015, 2015); Jacquemin et al. (2015); Tiago and Chelikowsky (2006); Körbel et al. (2014) By placing low energy transitions in , is effectively a large band gap system, and we expect the correlation to be weak. As long as has well-defined and long-lived quasiparticles, we can justify the static approximation. The full frequency dependence of is included in the self-energy but neglected for inter-quasiparticle interactions.

We can finally write our approximation to the excitation matrix , which is equivalent to the Hamiltonian in Eq. 29. We limit to a diagonal approximation. The diagonal approximation in is the major computational savings of our approach, as demonstrated in Fig. 8. In such an approximation, there is no coupling between forward and backward time orderings of . The Tamm-Dancoff approximation (TDA) makes connections to the resolvent of the auxiliary eigenvalue problem (Appendix C) and projected Hamiltonian easier, since neither quantity has any time ordering. Even in a diagonal approximation to in the TDA, however, correlation is described at the quasiparticle level. Diagonal matrix elements of are


Here, we reintroduce the configuration notation to make the connection to the wave function picture clear. The basis sets for the WF and GF calculations are now connected. The sums in Eq. 30 run up to the excitation level of the configuration. and denote the excited electrons and holes of the configuration they are exactly the creation and destruction operators in Eq. 16 that define the configuration . The final result to compute excitation energies is relatively straightforward: the one-body part of the effective Hamiltonian is quasiparticles, and the two-body part is their interaction via the screened interaction . Eq. 30 is not a diagonal approximation to the BSE, though it is BSE-like.

Figure 8: The diagonal approximation to the transformation dramatically reduces the computational cost while keeping a quasiparticle description of correlation. Coupling between the two spaces is computed in the WF picture through matrix elements of the exact , and . The WF description of hybridization avoids frequency dependent hybridization between and . The MBPT calculation, and therefore the frequency dependence of the perturbation expansion, is contained in the subspace. The physical in the denominator of determines the coupling strength between the two spaces, but dependence is not part of the perturbation expansion.

ii.4 Ground state energy

To complete the transformation , we must also estimate the ground state energy . We discuss formal limits and definitions, then our procedure for calculating .

ii.4.1 Definition and embedding limits

Consider an eigenstate of given by with eigenvalue . is on the order of a total energy of the system. Assume has a corresponding long-lived excitation computed with MBPT (as outlined previously) with excitation energy . The ground state energy is chosen to match the total energy :


We assume a single value exists to match Eq. 31 for all . If we exactly diagonalize , we can compute as the difference . Of course, this defeats the purpose of the transformation since we want to avoid the expense of diagonalizing . We must develop a different, more efficient strategy to estimate .

is difficult to define beyond Eq. 31. Because contains only excitations, there simply is no ground state in that can be related to the physical ground state. It is still possible, however, to choose a scalar value to satisfy Eq. 31. While it is unintuitive that this value does not correspond to an eigenstate, such a correspondence is not necessary to satisfy Eq. 31. We can still place some constraints on . In the limit that , we must recover a standard MBPT calculation with the full . In this case, we know that excitation energies are defined with respect to the fully correlated ground state, so .

The opposite limit also gives some insight. While the case trivially recovers FCI, consider one configuration in . This special case is considered in more detail in Appendix B. In this case, the exact diagonalization result is simply one diagonal matrix element. It can be decomposed into a clearly defined ground state, with energy equal to the energy of the reference configuration, , and excitation energy. The details are presented in Appendix B, but this case provides a critical second constraint: as . At intermediate sizes of , we enforce that must interpolate between these two limits: .

ii.4.2 Partitioning correlation energy

For our estimate of , we assume a nondegenerate ground state and partition the correlation energy of the true ground state into three distinct parts.


Here, is the intra- correlation, is the inter-space correlation, and is the intra- correlation. As and change size, the total correlation energy () must stay constant the total energy is independent of the partitioning between and . However, the individual contributions are allowed to change as the correlation is transferred between spaces. In the limit , ; in the limit , .

The change in the correlation energy mimics the desired behavior for the energy . We propose using


to estimate the ground state energy . If we can isolate the intra- correlation, Eq. 33 correctly interpolates between our set limits on . The limits are summarized as


To this end, we introduce another total energy, , given by


We see that, if it is possible to compute the energy , we can isolate and calculate as


Here, we assume that the various correlation energies are the same for all three total energies. Generally, this is not true. It is not possible to simply “turn off” a given part of the correlation as it appears in the full problem. However, it is an approximation that can be aided by intuition about choosing the AS.

It will prove useful to identify where these three contributions to the correlation exist in the partitioned problem of Eq. 6. exists in the projected block of , the frozen core Hamiltonian. The hybridization exists in matrix elements of (and ). Finally, the intra- correlation is in matrix elements of . This portion of the correlation is contained in the resolvent, which depends on the inversion of the block .

ii.4.3 Uncorrelated resolvent

With these considerations in mind, our strategy to estimate depends on an estimate of the total energy . To calculate the energy , we self-consistently solve Eq. 6 using an uncorrelated resolvent . By using an uncorrelated Hamiltonian in place of , we approximately remove the intra- correlation. It is an approximation that the correlation removed by this procedure is the same as the correlation energy of the real ground state, .

We choose the uncorrelated Hamiltonian based on particle addition or removal to the reference configuration, with no interactions among the added particles or holes. The potential must be fixed at the field created by the reference configuration. In the case of a Hartree-Fock (HF) mean-field starting point, which we use exclusively in this work, Koopman’s theorem dictates that these particle addition energies are the HF eigenvalues. Matrix elements of the uncorrelated Hamiltonian (which is diagonal) are


With the uncorrelated resolvent, matrix elements of in Eq. 6 are not zero. They still depend on the inter-space elements of and . The purpose of this procedure is to remove the intra- correlation in the resolvent, so that contains only the hybridization . is the self-consistent solution of Eq. 6 with the uncorrelated resolvent,


While this describes the procedure to estimate in this work, there are other possible routes to calculate . It could be possible to estimate with diagrammatics. For example, one could compute a correlation energy based on the cRPA. For now, we favor a more diagonalization-based approach. The RPA is known to introduce spurious bumps in dimer dissociation curves that could impact our test calculations.

ii.5 Final equations

After calculating the excitation matrix and energy , the effective Hamiltonian becomes


With this transformed Hamiltonian, the resolvent is


Combining the different total energies, we rewrite the resolvent as


where and . now has the interpretation of an optical frequency or excitation energy. Here, is the physical frequency related to exciting the system. The ground state (