A real-time extension of density matrix embedding theory for non-equilibrium electron dynamics

Non-equilibrium electron dynamics is prevalent throughout chemistry and physics, for example, in electron transport through molecular junctions,^1–4 electron injection and transport following photoexcitation,^5,6 and driven electron dynamics in laser pulses.^7,8 The simulation of such processes is challenging due to the need to treat both large system sizes and electron correlation. A variety of methods have been developed for non-equilibrium electron dynamics, including non-equilibrium Green’s function approaches,^9–16 numerical path-integral techniques,^17–21 real-time Monte Carlo methods,^22–29 semiclassical approximations,^30–32 and wavefunction propagation methods.^33–53 In this work, we will present new developments in the latter class.

Most wavefunction-based methods fall into two categories: those which sacrifice accuracy for the ability to treat large system sizes, such as time-dependent density functional theory^{44,45,54–57} and time-dependent Hartree-Fock theory,^46–48 and methods which are highly accurate, but are limited to small system sizes, such as multi-configurational time-dependent Hartree theory^33,34 and the time-dependent density matrix renormalization group (DMRG).^{49–53,58,59} To improve the compromise between accuracy and efficiency, we here borrow an idea from electronic structure approximations, namely, that of quantum embedding. Embedding techniques work by dividing the total system into a small region of interest, termed the impurity, which is treated accurately, and the surrounding environment, which is treated in an approximate manner. This decomposition allows calculations on a large total system, while retaining a high level of accuracy in the region of interest.

One powerful embedding formulation that has been introduced for static electronic properties is the density matrix embedding theory (DMET).^60–62 In DMET, the surrounding environment is represented by a quantum bath, constructed to capture the entanglement between the environment and the impurity. The entanglement-based construction ensures that the size of the quantum bath is at most equal to the size of the impurity. The bath allows for strong coupling between the impurity and environment, while its small size ensures computational efficiency. In its single impurity formulation, the DMET wavefunction can be viewed as a kind of complete active space (CAS) wavefunction. In the more general multi-impurity DMET formulation, the total system is divided into multiple local impurities, each associated with its own embedded problem or related to each other by translation in a crystal environment. DMET in both its single impurity and multiple impurity versions has been successfully applied to fermion and spin lattice models^60,63–67 as well as ab initio molecular and condensed phase systems.^61,62,68,69

Here, we will develop a quantum embedding formulation for dynamics, by extending DMET to the real-time propagation of the electronic wavefunction. Specifically, we focus on a real-time extension of the single impurity formulation of DMET as a first step in the development of the multi-impurity time-dependent theory. The main conceptual difference with the static DMET is that the quantum bath becomes time-dependent. We use the time-dependent variational principle (TDVP)^70–72 to derive the dynamics of the quantum bath as well as that of the correlated wavefunction in the coupled impurity-bath problem. To maintain a constant Hilbert space for the impurity, we introduce a constraint in the TDVP that the impurity orbitals remain time independent. This provides the basis for a simple future generalization to a multi-impurity theory. The real-time DMET possesses analogous formal strengths to the original static formulation and provides an exact description of dynamics in the non-interacting, isolated cluster and large impurity size limits.

We demonstrate the strengths of the method by simulating several kinds of quantum quenches. These are defined by a sudden change in the Hamiltonian, which induces a non-equilibrium state in the system and thus subsequent dynamics. As a test system, we use the single impurity Anderson model (SIAM), and we compare our results, where possible, against time-dependent Hartree Fock (TDHF) theory, time-dependent complete active space self-consistent field (TD-CASSCF) theory, and numerically exact time-dependent density matrix renormalization group (TD-DMRG) benchmarks.^58,73 We find excellent numerical agreement between real-time DMET and the numerically exact TD-DMRG, including in the regime of the Kondo resonance. Overall, our results show that DMET offers an accurate treatment of the quantum dynamics in the impurity region, with a very affordable cost.

We begin by reviewing the static DMET algorithm to provide a foundation for the later presentation of real-time DMET. For simplicity, we will assume that the static problem of interest is the ground-state problem and will focus only on the “interacting bath” formulation of DMET.⁶² We also restrict our discussion to a single embedded impurity formulation, as described below.

Consider a full quantum system spanned by an orthonormal single-particle basis, indexed by p, q, r, and s. The starting, local, single-particle basis of the problem will be referred to as sites, while more general single-particle functions will be termed orbitals. The total size of the basis will be denoted N. The general second-quantized Hamiltonian for the full system can be written as

where h_pq = ⟨p|ĥ|q⟩ and $Vpqrs=⟨pq|V^|rs⟩$ are the one- and two-electron Hamiltonian matrix elements,

and

The operator $apσ†$ (a_pσ) creates (destroys) an electron of spin σ at site p.

The single particle basis of the full quantum system can be partitioned into a small subset of sites, i ∈ A, corresponding to the region of interest and termed the impurity; the number of these sites will be denoted N_imp. The remainder of the sites constitute the environment. Static DMET relies on the observation that, for any state of the full quantum system, the entanglement between the impurity and the surrounding environment can be exactly accounted for by a quantum bath that is the same size as the impurity. Specifically, given the exact ground state of the full quantum system, it can be written through its Schmidt decomposition as

where ψ_i is an expansion coefficient, $αi$ are (multi-electron) states in the Fock space spanned by the impurity A, and $βi$ are (multi-electron) states in the Fock space of the environment that constitute the quantum bath. Note that though the states $βi$ fully capture the entanglement with the environment, there are only M_A of them: the dimension of the Fock space of the impurity. Thus, in principle, the ground state can be determined by solving the Schrödinger equation with Ĥ projected into the small impurity plus bath Hilbert space. However, this is not a practical solution, as the definition of the environment states $βi$ requires knowledge of the exact solution. To circumvent this, in static DMET, the states $βi$ are calculated instead from the ground state of a simpler Hamiltonian, ĥ′. The static DMET approximations to the ground state and expectation values of the original interacting problem thus require self-consistently solving two coupled models: (i) for the ground state, $Φ$⁠, of the approximate Hamiltonian, ĥ′, in the full system Hilbert space and (ii) for the ground state, $|Ψimp$⁠, of the interacting problem, within the small embedding Hilbert space of the impurity coupled to the now approximate quantum bath. A self-consistency condition on the one-particle reduced density matrix links the two models.

In most applications of static DMET, the approximate Hamiltonian for the full quantum system is defined as a single-particle Hamiltonian of the form

where ĥ is most commonly chosen to be either the one-particle part of the total Hamiltonian, Ĥ, or the Fock operator derived from Ĥ. In the case of purely local interactions, as in the Anderson impurity model studied in this work, the two choices are equivalent. Here, û is the local correlation potential,

which approximates the effect of the local Coulomb interaction within the impurity. The sum in Eq. (6) is restricted to the N_imp sites within the impurity A. The elements u_pq are obtained through the self-consistency condition described below.

The quantum bath that defines the embedding problem is obtained from the ground state of ĥ′, $Φ$⁠, which takes the form of a simple Slater determinant. In this case, the multi-electron states, $βi$⁠, assume a particularly simple form: they constitute states in the Fock space spanned by a set of at most N_imp bath orbitals, multiplied by a core determinant. These embedding orbitals can be obtained in several mathematically equivalent ways. Here, we will assume that the bath orbitals are determined by diagonalizing part of the one-particle density matrix, ρ^Φ, computed from $Φ$⁠, with elements $ρpqΦ=⟨Φ|Eqp|Φ⟩$⁠. The one-particle density matrix can be partitioned into a N_imp × N_imp impurity block, a (N − N_imp) × (N − N_imp) environment block, and N_imp × (N − N_imp) off-diagonal coupling blocks,

Diagonalizing the environment block of the one-particle density matrix, $ρenv=RenvΛRenv†$⁠, yields three kinds of embedding orbitals: (i) a set of N_impbath embedding orbitals with eigenvalues between zero and two, that describe the entanglement between the environment and impurity, (ii) a set of N_occ − N_impcore embedding orbitals, where N_occ is the number of occupied orbitals in $Φ$⁠, with eigenvalues equal to two, that are thus not entangled with the impurity; these orbitals comprise the core determinant, and (iii) a set of N − N_imp − N_occvirtual embedding orbitals, with eigenvalues equal to zero, that are thus also not entangled with the impurity. The embedding problem thus consists of a complete active space (CAS) wavefunction calculation in which the impurity and bath orbitals comprise the active space, and the core determinant is comprised of the core embedding orbitals; the virtual embedding orbitals constitute the space external to the active and core spaces.

Before continuing, we introduce some notation for the embedding orbitals: impurity orbitals will be designated by indices i and j; embedding bath orbitals will be designated by y and z; the combination of all active space orbitals, which correspond to both the impurity and bath orbitals, will be designated by l, k, m, and n; embedding core orbitals will be designated by u, $v$⁠, w, and x; embedding virtual orbitals will be designated by indices a and b; the combination of all single-particle orbitals in the embedding basis will be designated by c, d, e, f, and g.

In the interacting-bath formulation of static DMET, the Hamiltonian of the embedding problem, Ĥ_imp, is obtained by projecting the original fully interacting Hamiltonian, Ĥ, into the active-space defined in the embedding basis, and including the contribution from the doubly occupied core determinant. This can be performed by a change of single-particle basis from the original site basis to the embedding basis while including a contribution from the core orbitals, such that

where

and

The rotation matrix from the site basis to the embedding basis is given by

where the identity matrix denotes that the impurity orbitals are the same in the original site basis and the embedding basis; R_env is defined above.

A wide range of solvers can be used to compute the correlated ground state, $|Ψimp$⁠, of the embedding Hamiltonian, Ĥ_imp, depending on the nature of the problem as well as the cost and accuracy requirements. In this work, we use exact diagonalization as the impurity solver, though previous work has also employed DMRG,^65–67 coupled cluster theory,^62,66 and auxiliary-field quantum Monte Carlo.⁶³ 

As described above, the elements of the correlation potential û are determined by a self-consistent procedure. Specifically, we minimize the difference between the impurity block of the one-body density matrices calculated from the uncorrelated wavefunction, $Φ$⁠, and correlated wavefunction, $|Ψimp$⁠,

and the elements $ρijΨimp=⟨Ψimp|Eji|Ψimp⟩$⁠. However, as in previous work,^61–63 the functional f(u) is not directly optimized, but instead a self-consistent iteration is used: f(u) is optimized with a fixed $|Ψimp$⁠; the optimal u is then used to update $|Φ$⁠, the embedding Hamiltonian, Ĥ_imp, and thus $|Ψimp$⁠.

In summary, the static DMET algorithm proceeds via the following steps:

1. we choose an initial guess for the correlation potential û;

2. we solve for the approximate Hamiltonian, ĥ′, to obtain the reference wavefunction $Φ$⁠;

3. we construct the embedding Hamiltonian using Eq. (8);

4. we use exact diagonalization to compute the ground state of the embedding problem, $|Ψimp$⁠, and construct the one-body density matrix $ρΨimp$⁠;

5. we minimize f(u) in Eq. (13), with $ρΨimp$ fixed, to obtain a new correlation potential u′;

6. if ||u − u′||_∞ > ε₀, the convergence threshold, we set u = u′ and go to step 1, otherwise the static DMET calculation is converged.

We now describe the central methodological contribution of the paper, namely, a real-time extension of DMET, which allows for the efficient time-propagation of the electronic wavefunction. We focus on the propagation of only a single impurity and its corresponding quantum bath as a first step in the development of a multi-impurity time-dependent framework.

Developing a real-time extension of DMET entails discerning how to appropriately propagate (i) the embedding orbitals, to give the approximate representation of the time-dependent environment, and (ii) the correlated CAS-like DMET wavefunction in the embedding problem, to describe the region of interest at a high level. Here, we utilize the TDVP to derive the equations of motion for both the embedding orbitals and the expansion coefficients for the determinants in the DMET CAS-like wavefunction. We introduce the constraint that the impurity orbitals are time independent in keeping with an embedding picture in which one can cleanly identify the contributions of different impurities to global observables. While such a constraint is not strictly necessary in a single impurity formulation (one could after all, just use a time-dependent CASSCF wavefunction), we keep this constraint in anticipation of the future development of the multiple impurity formalism and to see its effect on the single impurity dynamics. Our TDVP derivation intrinsically connects the low-level orbital dynamics and the high-level embedded dynamics and thus does not require a further self-consistency through a time-dependent correlation potential. We return to the question of the challenges of a self-consistent picture from a time-dependent correlation potential in the  Appendix. We now present the detailed derivation of the equations of motion, which constitute the working equations for the real-time DMET method.

To begin, we write the correlated wavefunction for the embedding problem as a time-dependent CAS wavefunction,

where $M(t)$ are time-dependent determinants in the active space defined by the impurity and bath orbitals coupled to a doubly occupied determinant comprised of the embedding core orbitals, and C_M(t) are the time-dependent expansion coefficients. The time-dependence of the determinants arises from the time-dependence of the embedding bath and core orbitals, as the impurity orbitals are kept time independent.

Following the TDVP,^{39,40,70–72} the equations of motion for both $M(t)$ and C_M(t) can be obtained by varying the Dirac-Frenkel action with fixed endpoints,

This procedure yields the variational equation

which must be satisfied for arbitrary variations of the wavefunction, δΨ_imp. We should note that the ground-state DMET wavefunction is not rigorously a stationary state of the TDVP since the ground-state DMET is not variationally optimized. However, numerical investigation seems to suggest that the difference is small.

The variation of the wavefunction with respect to the expansion coefficients can be written as

while the variation with respect to the embedding orbitals can be written as⁴⁰ 

where Δ_ab is an anti-Hermitian matrix. The complete time-dependence of the wavefunction can be expressed as

where we have introduced the single-particle Hermitian operator $X^$⁠, which governs the time-dependence of the embedding orbitals. The operator is defined as

where the elements $Xcd=iℏ⟨c|ḋ⟩$ are determined through the variational equation, Eq. (16), as shown below.

Inserting the variation with respect to the expansion coefficients, Eq. (17), and the time dependence of the wavefunction, Eq. (20), into the variational equation, Eq. (16), yields

which defines the equations of motion of the expansion coefficients. Inserting the variation with respect to the orbitals, Eq. (18), yields

where $Π=∑MMM$ is the projector into the CAS space defined by the impurity and embedding orbitals. Solving Eq. (23) defines the elements of the operator $X^$⁠.

The elements of $X^$ will now be derived for each type of orbital rotation in the embedding basis, utilizing the notation for the different kinds of embedding orbitals defined in Sec. II A. Equation (23) reduces to a trivial identity for an orbital pair {c, d} corresponding to the same orbital subspace (core, active, or virtual) due to the presence of the projector (1 − Π) out of the CAS space;^39,40 the determinants $ÊuvM=2δuvM$⁠, $ÊlkM$⁠, and $ÊabM=0$ are either zero or fall within the CAS space and are thus eliminated by the projector (1 − Π). These intraspace orbital rotations, Ê_cd = {$Êuv$⁠, Ê_lk, Ê_ab} are referred to as redundant, since the total wavefunction is invariant to such rotations if accompanied by the corresponding transformation of the expansion coefficients as seen by the presence of $X^$ in Eq. (22).^39,40,74,75 The elements of $X^$ for these redundant orbital pairs can then be freely chosen; in this work, we set these terms to zero such that $Xuv=Xvu*=Xlk=Xkl*=Xab=Xba*=0$⁠.

For the non-redundant orbital pairs, the projector, (1 −Π), in Eq. (23) can be dropped and the equation reduces to

where the 2-particle reduced density matrix has elements $ΓcdefΨimp=⟨Ψimp|Eefcd|Ψimp⟩$⁠. Equation (24) can now be solved for each non-redundant orbital rotation.

As mentioned above, the impurity orbitals in real-time DMET are restricted to be time independent. Therefore, all elements of $X^$ that include an impurity orbital are defined to be zero, such that $Xic=Xci*=0$⁠.

The other non-redundant orbital rotations can be obtained using the non-zero elements of the reduced density matrices for the embedding wavefunction, which are $ρuvΨimp=2δuv$⁠, $ρlkΨimp$⁠, $ΓuvwxΨimp=4δuwδvx−2δuxδvw$⁠, $ΓlukuΨimp=ΓulukΨimp=2ρlkΨimp$⁠, $ΓluukΨimp=ΓulkuΨimp=−ρlkΨimp$⁠, and $ΓlkmnΨimp$⁠. This then yields

and

where the matrix $ρbathΨimp$ corresponds to the bath block of the 1-electron reduced density matrix and $ρ¯bathΨimpyz=2δyz−ρyzΨimp$⁠.

The real-time DMET equations of motion can now be written as

for the embedding orbitals, where the elements X_dc are given in Eqs. (25)–(27), and

for the wavefunction coefficients, where Eq. (29) is obtained from Eq. (22) by noticing that the matrix elements $⟨M|X^|N⟩$ are non-zero only for intraspace rotations, and those components of $X^$ are all defined to be zero.

It is important to note that the equations of motion for real-time DMET are very similar to those derived for TD-CASSCF, which will be presented in Sec. II C.^35,39,40 The main difference is that in real-time DMET a subset of the active space orbitals, specifically the impurity orbitals, are restricted to be time independent. This allows us to cleanly identify correlated impurity states during the dynamics, which gives a natural route to extending the methodology to treat multiple impurities. Analogous to the static DMET, the real-time DMET is exact for the impurity properties in the non-interacting limit, in the limit when the size of the impurity becomes (half) the size of the full quantum system, and in the limit where there is no coupling between the impurity and the environment. The latter exact property is not ensured by TD-CASSCF.

We conclude Sec. II by introducing the working equations of motion for TD-CASSCF. The derivation of these equations has been worked out previously and bears close resemblance to the derivation of the real-time DMET equations of motion, so here we only provide the final results.^35,39,40

The time-dependent wavefunction in TD-CASSCF takes an analogous form to Eq. (14),

where $M(t)$ are time-dependent determinants comprised of the active space coupled to the doubly occupied core orbitals, and C_M(t) are the time-dependent expansion coefficients.

The equations of motion for the expansion coefficients and single-particle orbitals are

and

respectively, where, as in Sec. II B, the single-particle Hermitian operator $X^$ governs the time dependence of the orbitals. In comparison to real-time DMET, however, in TD-CASSCF, all of the active space orbitals are time dependent.

The non-redundant elements of the operator $X^$ are given by

and

where $ρijΨCAS=⟨ΨCAS|Eji|ΨCAS⟩$ and $ρ¯ΨCASyz=2δyz−ρyzΨCAS$⁠. All other elements of the operator $X^$ constitute redundant orbital rotations and are thus set to zero as in Sec. II B. Mirroring the notation from Secs. II A and II B, in this section, active space orbitals are designated by, j, k, l, m, and n; doubly occupied core orbitals are designated by u and $v$⁠; unoccupied virtual orbitals are designated by a; the combination of all single-particle orbitals is designated by c and d.

Equations (33)–(35) involve the inverses of $ρΨCAS$ and $ρ¯ΨCAS$⁠, which become numerically unstable when a subset of the single particle orbitals in the active space become fully unoccupied or occupied, respectively. To avoid this numerical instability, we regularize both matrices following a similar procedure to what is done in multiconfiguration time-dependent Hartree theory.^76,77 Specifically, any eigenvalue of the matrices below a small threshold ε is set to ε.

In this work, we will benchmark the real-time DMET on the non-equilibrium dynamics of the single impurity Anderson model (SIAM) following different quantum quenches. The SIAM is a model of an interacting impurity embedded in a non-interacting environment and can be realized in different physical systems, such as in quantum dots or molecules attached to metallic leads.^78–80 As the simplest example of a bulk interacting quantum problem, it provides a useful benchmark system for non-equilibrium electron dynamics in the presence of electron correlation. We emphasize, however, that real-time DMET is applicable to more general and complex systems.

The SIAM consists of a single quantum dot site, where Coulomb interactions are present, coupled to two non-interacting leads, Fig. 1. In this work, we use a real-space definition of the SIAM, in which the leads have nearest neighbor hopping terms, to allow for easy comparison with real-time DMRG calculations.⁵⁸ The leads are finite in size and the total system size including the leads and interacting site is N. The SIAM under a bias is then described by the Hamiltonian Ĥ = Ĥ_dot + Ĥ_leads + Ĥ_dot−leads + Ĥ_bias where

describes the quantum dot in isolation. The quantum dot is located at site d = N/2, V_g is the gate potential which controls the location of the energy level of the quantum dot, U is the local Coulombic interaction, n_d = n_d↑ + n_d↓, and $ndσ=adσ†adσ$⁠. The Hamiltonian of the leads in isolation is

where t_leads is the hopping amplitude of the lead and the subscript Lp (Rp) denotes site p in the left (right) lead. The quantum dot is coupled to the two leads through the term

where t describes the hopping amplitude between the surrounding leads and the quantum dot and the subscript L1 (R1) denotes the lead site that is closest to the quantum dot in the left (right) lead. Finally, a bias can be applied across the SIAM, of the form

In this work, we investigate the non-equilibrium dynamics following two types of quantum quenches, depicted in Fig. 1. In the first, the initial state for the subsequent dynamics is defined as the ground state of the SIAM under zero bias, ΔV = 0; the dynamics of the initial state are then propagated using a Hamiltonian including a finite bias, ΔV′ ≠ 0, which drives a current through the quantum dot. In the second, the initial state is defined as the ground state of the SIAM with a specific value of the local Coulomb interaction, U; the dynamics are then propagated using a Hamiltonian with a different interaction, U′ ≠ U.

The dynamics following the quantum quenches are characterized through several observables. Specifically, we investigate the time dependence of the occupancy on the quantum dot, n_d, and the time-dependent current through the dot, which is defined as the average of the current between the dot and the closest left lead-state and closest right lead-state, $J(t)=JL(t)+JR(t)/2$⁠, where^33,58

and

The use of the symmetrized current provides better numerical convergence to infinite system size, particularly when the left and right leads have a different number of sites.^33,58 The conductance, G, can be obtained by dividing the steady-state value of the current by the total bias applied across the leads during the dynamics, ΔV′.

For calculations utilizing real-time DMET, the initial state of the system is calculated using the static DMET algorithm described in Sec. II A using exact diagonalization as the impurity solver. The subsequent dynamics of the electronic wavefunction are propagated using the real-time DMET equations of motion, Eqs. (28) and (29), which are evaluated using the fourth-order Runge-Kutta method. All results are fully converged with a time step of 0.005.

We present real-time DMET results with varying impurity size, N_imp, which is kept the same between the initial static and subsequent real-time DMET calculations. The impurity always includes the quantum dot, followed by lead states in increasing distance from the quantum dot; a left-lead state is always included prior to a right-lead state, though results are relatively insensitive to this choice. Thus, an impurity of size N_imp = 3 includes the quantum dot and the closest lead-state from the right and left leads, while an impurity size of N_imp = 4 includes the quantum dot, the two closest lead-states from the left lead, and the closest lead state from the right lead.

For calculations utilizing TD-CASSCF, the initial state of the system is calculated using CASSCF as implemented in the PySCF package.⁸¹ The subsequent dynamics of the electronic wavefunction are propagated using the TD-CASSCF equations of motion, Eqs. (31) and (32), which are evaluated using the fourth-order Runge-Kutta method. All results are fully converged with a time step of 0.0001 and a value of the regularization parameter ε = 10⁻⁸.

We present TD-CASSCF results with varying active space size, N_CAS, which is kept the same between the initial static and subsequent time-dependent CASSCF calculations. The active space is always chosen to include the HOMO through HOMO-N_CAS/2 and LUMO through LUMO+N_CAS/2 orbitals from the Hartree-Fock calculation used to initialize the CASSCF calculation.

As a reference, we compare to results generated using TD-DMRG either computed using the ITensor library⁷³ with a bond dimension of 300 or from previous work.⁵⁸ Comparisons are also made to results generated using TDHF theory.^70,82,83

We now present our results utilizing real-time DMET to simulate the non-equilibrium electron dynamics in the SIAM following a variety of quantum quenches. In all cases, the parameter t_leads = 1.0 is taken as the energy scale and the parameter t = 0.4; the parameters V_g, U, U′, ΔV, and ΔV′ are varied to define a wide range of quantum quenches.

As mentioned in Sec. II, real-time DMET is exact in the non-interacting limit. This is numerically verified in Fig. 2 in which we present results for a non-interacting quantum quench in which a bias is suddenly switched on to drive current through the quantum dot; the parameters are U = U′ = 0, V_g = 0, ΔV = 0, and ΔV′ = −0.001. Figure 2 illustrates that the current through the quantum dot evaluated using real-time DMET (open circles) exactly matches results from exact dynamics (solid lines) for a range of total system sizes: N = 16 (green), N = 32 (blue), N = 64 (red), and N = 128 (cyan). The real-time DMET calculations use an impurity size of N_imp = 3; real-time DMET is exact in the non-interacting limit regardless of impurity size. The exact dynamics are obtained by integrating the equations of motion for the one-electron reduced density matrix of the total system, $iℏρ̇pq=∑rρprhrq−hprρrq$⁠.

Figure 2 also illustrates an important result regarding the non-equilibrium electron dynamics in a finite-size SIAM. Specifically, the current is seen to oscillate for small total system size. This can be attributed to a recurrence of the electron density following a reflection off of the end of the leads. The position and height of the recurrence provide a metric to benchmark the dynamics generated using real-time DMET when it is not exact, similar to using a Loschmidt echo. In addition, the steady-state behavior of the SIAM is defined as the plateau regime in between recurrences.

We now turn our attention to interacting quenches in which real-time DMET is only exact in the large impurity size limit. Figure 3 presents the time-dependent occupancy on the quantum dot, n_d(t), for a quantum quench in which the local Coulomb interaction on the quantum dot is suddenly switched on; the parameters are U = 0 and U′ ≠ 0 with V_g = ΔV = ΔV′ = 0. Such a quantum quench provides a useful benchmark for real-time DMET since the initial state is a non-interacting ground state; the initial state can thus be calculated exactly, such that the embedding orbitals are the exact embedding orbitals at time t = 0. The subsequent dynamics provide a test solely of the accuracy of the real-time DMET equations of motion. The figure compares results calculated using real-time DMET with an impurity size of N_imp = 3 (green), N_imp = 4 (blue), N_imp = 5 (red), and N_imp = 6 (cyan), and results calculated using TD-CASSCF with an active space size of N_CAS = 6 (dashed-green) and N_CAS = 8 (dashed-blue) to those obtained from time-dependent DMRG (black),⁷³ which can be taken as the exact answer. As stated in Sec. II A, the size of the active-space in real-time DMET is given by twice the size of the impurity, such that a real-time DMET calculation using an impurity size of N_imp is most fairly compared against a TD-CASSCF calculation using an active-space size of N_CAS = 2N_imp. Results calculated using TDHF (violet) are also presented as a reference mean-field calculation.

Figure 3(a) presents results for a small value of the Coulomb interaction, U′ = 1.0, and for a small total system size, N = 16. The numerically exact TD-DMRG dynamics are characterized by a rapid decrease in the quantum dot occupancy, n_d(t), followed by high-frequency, small amplitude, oscillations around a plateau value, with recurrence peaks occurring at t ≈ 20 and t ≈ 40 associated with the small total system size. Due to the severe approximations present in the method, the mean-field TDHF results overestimate the decrease in the time-dependent occupancy and do not correctly capture the high-frequency oscillations even for the relatively small value of the Coulomb interaction; TDHF, however, is able to correctly capture the position of the recurrence peaks. In comparison, real-time DMET quantitatively captures the decrease in the occupancy as well as the position and height of the recurrence peaks even for small impurity sizes, illustrating the ability of the real-time DMET method to capture the dominant time scales in the non-equilibrium electron dynamics. The high-frequency oscillations are captured at short times, t < 20, by real-time DMET, though the agreement becomes worse with increasing time. The worse agreement for the high-frequency oscillations can be attributed to finite-size effects associated with the impurity; since the dynamics of the embedding states are approximate, the embedding states are unable to fully damp the recurrence dynamics present within the impurity, which leads to artificial high-frequency oscillations. As seen in Fig. 3, and as will be illustrated further below, these artificial high-frequency oscillations disappear with increasing impurity size.

Figure 3(a) also demonstrates that TD-CASSCF is slightly more accurate than single-impurity real-time DMET, exhibiting near-exact agreement with TD-DMRG with a small active space. The TD-CASSCF equations of motion do not contain the constraint that the impurity orbitals must remain time independent as is the case for real-time DMET; this increased dynamic flexibility leads to a better approximation of the time-dependent wavefunction. This mirrors what is observed for static DMET; for the case of a single impurity, DMET does not optimize the partitioning between the impurity and environment, while this partitioning is fully optimized in a variational CASSCF. However, the presence of this fixed partitioning simplifies the construction of the multi-impurity form of DMET, where all sites can be correlated unlike in CASSCF. Importantly, Fig. 3(a) shows that the cost of the impurity constraint in real-time DMET is not large; thus calculations with larger impurities still rapidly converge to the exact answer.

Figure 3(b) presents results for the same small value of the Coulomb interaction, U′ = 1.0, but for a larger total system size, N = 64. For this system size, recurrence peaks are no longer present in the TD-DMRG dynamics for the time scale pictured. Instead, the time-dependent occupancy is characterized by a rapid decay to a plateau value. Once again, the TDHF results underestimate the plateau value. Both real-time DMET and TD-CASSCF show a marked improvement over TDHF, correctly capturing the time scale for the rapid decay and the magnitude of the plateau value. For smaller impurities, the real-time DMET dynamics exhibits artificial high-frequency oscillations in the plateau region, t > 5, for small impurity sizes. The magnitude of these oscillations can be clearly seen to diminish with increasing impurity size, and the N_imp = 6 results exhibit no oscillations whatsoever.

Figure 3(c) present results for a larger value of the Coulomb interaction, U′ = 3.0, and for the larger system size, N = 64. This provides a more stringent test for real-time DMET as a larger Coulomb interaction leads to stronger correlation between the impurity and the surrounding environment. The TD-DMRG results show similar behavior to the U′ = 1.0 results; the results are characterized by a rapid decay followed by a plateau. In comparison to the U′ = 1.0 results, though, the dynamics also exhibit high-frequency oscillations throughout the entire trajectory. Expectedly, the TDHF results are worse than those obtained in Figs. 3(a) and 3(b); the mean-field approximation becomes worse for higher values of the Coulomb interaction due to the presence of strong correlation. In comparison, both real-time DMET and TD-CASSCF are again able to correctly capture the rapid decay of the occupancy and the value of the plateau.

Taken together, the results in Fig. 3 clearly illustrate that real-time DMET with very small impurities provides a qualitatively correct picture of the non-equilibrium dynamics, and that accuracy increases rapidly with the impurity size. In fact, the N_imp = 6 results are in almost quantitative agreement with the TD-DMRG results for all of the presented cases. This illustrates the computational benefits of embedding since in comparison to the full exact calculation, which would involve N correlated sites, the real-time DMET method is able to provide the same description using only 2N_imp correlated orbitals (N_imp impurity and N_imp embedding bath orbitals).

In comparison to the calculations performed in Fig. 3, most realistic simulations will not involve an exactly computed initial state. As such, Fig. 4 presents the time-dependent occupancy on the quantum dot, n_d(t), for a quantum quench in which the local Coulomb interaction on the quantum dot is suddenly switched off; the parameters are U ≠ 0 and U′ = 0 with V_g = U/2 and ΔV = ΔV′ = 0. This quantum quench presents a situation in which the initial state obtained using static DMET (or CASSCF) is no longer exact; the subsequent dynamics are thus also not exact even though the dynamics are generated using a non-interacting Hamiltonian. The figure again compares results calculated using real-time DMET with an impurity size of N_imp = 3 (green), N_imp = 4 (blue), N_imp = 5 (red), and N_imp = 6 (cyan) to those we have calculated using TD-CASSCF with an active space size of N_CAS = 6 (dashed-green) and TD-DMRG (black),⁷³ which can be taken as the exact answer. Results using TDHF are no longer shown as even the initial static Hartree-Fock guess is too inaccurate.

The real-time DMET results in Fig. 4 show very similar behavior to that seen in Fig. 3. Figure 4(a) presents the results for a small total system size, N = 16, and a small value of the Coulomb interaction, U = 1.0. Again, the real-time DMET results are able to capture the position and amplitude of the recurrence peaks at t ≈ 20 and t ≈ 40 regardless of impurity size, while the small impurity size results show artificial high-frequency oscillations. Similarly, Fig. 4(b) illustrates that the real-time DMET results correctly capture the time scale of the rapid increase in occupancy on the dot and the correct height of the plateau region for a larger total system size, N = 64, regardless of impurity size. The small impurity size results exhibit oscillations, which disappear with increasing impurity size. Finally, Fig. 4(d) present results for U = 3.0 and a large, N = 64, total system size, respectively. As observed in Fig. 3, the agreement between real-time DMET and TD-DMRG is not as good as small impurity size for the larger value of U. However, the real-time DMET results are clearly seen to converge to the TD-DMRG reference results for N_imp = 6. As was seen in Fig. 3, Fig. 4 shows that TD-CASSCF converges faster than real-time DMET with respect to the size of the active space.

We conclude this section by investigating an interacting quantum quench in which a bias is suddenly switched on to drive current through the quantum dot. Figure 5 presents the time-dependent current through the quantum dot calculated using TD-DMRG from Ref. 58 (black) and real-time DMET with an impurity size of N_imp = 3 (green), N_imp = 4 (blue), N_imp = 5 (red), and N_imp = 6 (cyan); the parameters are U = U′ = 1.0, V_g = U/2, ΔV = 0, and ΔV′ = −0.005. Figure 5(a) presents results for a small total system size of N = 16. As seen in the non-interacting case, Fig. 2, the small total system size leads to oscillations of the current. The frequency of these oscillations are captured by real-time DMET regardless of impurity size; the amplitude of the initial oscillation is captured by all impurity sizes, while the amplitude of subsequent oscillations is accurately captured only by the larger impurity size calculations.

Figure 5(b) presents results for a larger total system size of N = 128. The TD-DMRG results exhibit a rapid increase of the current, followed by a plateau region characterized by low amplitude oscillations;⁵⁸ the oscillations in the plateau region have been extensively discussed in the literature^{58,59,84–86} and can be attributed to the level spacing within the leads associated with the finite size of the total system. The real-time DMET results are able to correctly capture both the time scale for the increase of the current and the plateau value of the current. It is important to emphasize that the capacity of real-time DMET to correctly capture the plateau value of the current under these conditions is a non-trivial result. In principle, the value of the gate-voltage, V_g = −U/2, should put the system in the conductance “valley,” such that no current should be observed. However, the Kondo effect, which arises from the coupling between the localized spin on the quantum dot and the conducting electrons in the leads, yields a non-zero value of the current.^87–90 The ability of real-time DMET to simulate this many-body effect illustrates the power of the method to treat dynamics in the presence of strong correlations between the impurity and its environment.

To further emphasize this point, Fig. 6(a) plots the conductance for the interacting quantum quench corresponding to Fig. 5(b) as a function of the gate potential V_g. As was done for the time-dependent DMRG calculations,⁵⁸ the conductance is calculated as the average over the oscillations in the plateau region of the time-dependent current; the error bars report on the error associated with this average. Figure 6(a) compares the conductance calculated using time-dependent DMRG⁵⁸ (black) and real-time DMET with an impurity size of N_imp = 3 (green), N_imp = 4 (blue), and N_imp = 5 (red). The figure illustrates that the conductance is correctly calculated using real-time DMET for all values of V_g even for the small impurity size of N_imp = 3. At V_g = −U/2, the exact conductance is 2e²/h; the real-time DMET reproduces this result just as accurately as the time-dependent DMRG.

To achieve the exact conductance, it is necessary to push to larger total system sizes. Figure 6(b) illustrates one of the main benefits of the real-time DMET method, the ability to treat significantly larger system sizes compared to non-embedding methods. Figure 6(b) presents the conductance as a function of total system size obtained from real-time DMET with an impurity size of N_imp = 3 (green), N_imp = 4 (blue), and N_imp = 5 (red). The conductance is obtained for the same interacting quantum quench as in Figs. 6(a) and 5(b) with V_g = −U/2. Figure 6(b) illustrates that by pushing to N = 256, the ideal limit of the conductance (black line) is recovered; results for smaller system sizes are also pictured to show the convergence with respect to the total system size.

Figure 6(c) illustrates that this behavior is not limited to small values of the Coulomb interaction. Figure 6(c) presents the conductance as a function of total system size for a large value of the Coulomb interaction U = U′ = 4.0; the color scheme is the same as above and the remaining parameters are V_g = −2.0, ΔV = 0, and ΔV′ = −0.005. The large value of the Coulomb interaction leads to finite size effects of the total system size due to the large size of the Kondo cloud. Figure 6(c) shows that the value of the conductance at N = 128, which was the size of the system used in Fig. 6(a) and corresponds to the largest system studied previously with time-dependent DMRG,⁵⁸ and is significantly below the ideal limit of 2e²/h. However, with DMET, we can easily increase the system size. The remaining points in Fig. 6(c) show that the conductance approaches the ideal limit as we increase the total system size from 128 to 512 sites.

Finally, Fig. 6(d) highlights why real-time DMET is able to capture the Kondo effect. Figure 6(d) presents the magnitude of the coefficients of an embedding bath orbital at a single point in time as a function of site index for the interacting quantum quench presented in Fig. 6(c); the total system size corresponds to N = 128. The figure shows that the bath orbital is delocalized across the leads, allowing for a proper treatment of the delocalized Kondo cloud.

In this work, we present an extension of the density matrix embedding theory (DMET) to simulate real-time non-equilibrium electron dynamics. Like in the static case, the real-time DMET method partitions the full system into an impurity and an environment. The environment is efficiently represented by a quantum-bath of the same size as the impurity. The dynamics of the embedding problem, which consists of the impurity coupled to the quantum bath, is obtained through the time-dependent variational principle, in which we introduce the constraint that the impurity single-particle orbitals are time independent; such a constraint maintains the definition of the impurity during the dynamics and ensures the highest fidelity representation of the Hilbert space is retained for the impurity region.

The accuracy of the real-time DMET method has been benchmarked through comparisons with time-dependent Hartree-Fock (TDHF) theory, time-dependent complete active space self consistent field (TD-CASSCF) theory, and time-dependent density matrix renormalization group (DMRG)^58,73 for a variety of quantum quenches in the single impurity Anderson model (SIAM). We have shown that real-time DMET shows a marked improvement over the mean-field TDHF, correctly capturing the rapid change in the occupancy of the quantum dot following the sudden switching on or off of the local Coulomb interaction on the dot regardless of the size of the impurity. Furthermore, real-time DMET rapidly converges with respect to the impurity size, to quantitatively capture recurrence peaks for small total system size, or the correct plateau behavior for large total system size, for a range of values of the Coulomb interaction. However, due to the constrained form of the DMET wavefunction in the impurity picture, the real-time DMET error is larger than that of TD-CASSCF with the same number of correlated orbitals.

Additionally, we illustrate that real-time DMET can describe the non-trivial Kondo behavior during an interacting quantum quench in which a bias is suddenly switched on across the leads in the SIAM. Real-time DMET, in agreement with time-dependent DMRG, exhibits a value of the conductance near the ideal limit, for small values of the Coulomb interaction and small total system sizes. However, in comparison to time-dependent DMRG, real-time DMET can also simulate the significantly larger system sizes necessary to recover the ideal limit of the conductance, even for large values of the Coulomb interaction.

Taken together, the results presented in this work illustrate the capability of the real-time DMET method for simulating non-equilibrium dynamics in which strong correlation plays an important role. In addition, we have established the starting point for the extension to a multi-impurity theory, similar to how static DMET has been most commonly employed. Other extensions can include the incorporation of finite temperature effects or the use of other wavefunction Ansätze, such as a coupled-cluster or matrix product state wavefunction. The latter approximations will allow for the simulation of much larger impurity sizes than those treated in this study.

This work was supported by the US National Science Foundation, through Nos. CHE-1265277 and CHE-1665333. We thank Miles Stoudenmire for help with ITensor. The PySCF modules used in this work were developed with support from the US National Science Foundation, through No. OAC-1657286.

As mentioned in Sec. II B, other real-time extensions of DMET are, in principle, conceivable. In this appendix, we introduce a real-time formulation that is more closely analogous to the static DMET, in which the equations of motion of the embedding orbitals are derived from a reference mean-field dynamics for the full quantum system, which is propagated in tandem with the correlated dynamics of the embedding system. A self-consistency condition in terms of the impurity density can then be introduced between the reference mean-field dynamics and the correlated dynamics of the embedding problem, enforced via a time-dependent correlation potential. One advantage of this picture is that it yields a strong global consistency condition when there are multiple impurities. However, we have found that $v$-representability problems almost always occur after sufficiently long-time propagations.

The derivation of this real-time extension of DMET entails deriving the equations of motion for (i) the mean-field reduced density matrix of the full quantum system, (ii) the correlated CAS-like DMET wavefunction in the embedding problem, and (iii) the embedding orbitals. The explicit time dependence of all terms is suppressed for clarity.

The time dependence of the mean-field reduced density matrix is given by

where ρ^Φ is the one-particle density matrix initially obtained from the static DMET calculation and analogous to static DMET, and ĥ″ is a single-particle Hamiltonian of the form

The time-dependent correlation potential

is distinguished from the correlation potential used in the static DMET calculation, and the elements $uppRT$ are obtained through the self-consistency condition described below.

The equations of motion for the correlated CAS-like DMET wavefunction, Eq. (14), are derived from the time-dependent Schrödinger equation such that

which yields

We have once again introduced a one-electron operator $X^$ which governs the time dependence of the embedding orbitals, such that the equation of motion of the embedding orbitals is

Here, as in Sec. II B, the impurity orbitals are restricted to be time independent, such that $Xic=Xci*=0$⁠. The remaining elements of $X^$ are derived using the condition that at all points in time, the embedding orbitals are eigenfunctions of the environment block of the mean-field one-particle density matrix, Eq. (7), such that

The elements of $X^$ are obtained by differentiating Eq. (A9) with respect to time and using Eqs. (A1) and (A8), yielding

for c ≠ d. The coupling block of the reduced density matrix, $ρcΦ$⁠, is given in Eq. (7), and the coupling block of the single-particle Hamiltonian, $hc″$⁠, is defined analogously. The diagonal elements, X_cc, are arbitrary and can be set to zero.

The elements of the time-dependent correlation potential, û^RT, are determined through a self-consistent procedure, such that, within the impurity, the mean-field density and the correlated density obtained from $|Ψimp$ are equivalent at all times. This condition is achieved by minimizing the difference between the second time derivatives of the two densities on the impurity; the first time derivative of the mean-field density is independent of the correlation potential. However, this condition leads to $v$-representability problems, in which the second time derivative of the mean-field density becomes independent of the correlation potential. The second time derivative of the mean-field density on the impurity is given by

Taking the derivative of Eq. (A11) with respect to an element of the correlation potential yields

which can clearly be seen to go to zero for specific values ρ^Φ. Numerically it is observed that these special values of ρ^Φ are almost always obtained for sufficiently long dynamics indicating that this formulation of real-time DMET has a $v$-representability problem.

Equation (A12) is analogous to the force equation used in the Runge-Gross time-dependent density functional theory derivation.⁴⁴ In that case, however, the right-hand side can be written as −ρ∇u under the assumption that ρ_r,r+δr = ρ_r,r for infinitesimal δr due to the continuity requirements on the single-particle density matrix, and thus, cannot vanish except when ρ = 0, ensuring $v$-representability. The more severe condition encountered in the lattice formulation is due to the lack of continuity requirement on the density matrix as has been observed previously.⁹¹