A projection-based reduced-order method for electron transport problems with long-range interactions

The availability of high-performance computing (HPC) has prompted substantial progress in electron structure calculations using sophisticated models. Tremendous efforts have been made to implement the time-dependent density-functional theory (TDDFT)^1,2 in the context of electron transport.^3–9 Numerous software packages have been developed to facilitate such effort.^10,11 The main computational challenge in such an approach is still the computational cost. A typical situation is a molecular junction, connected to two leads that can be regarded as a quantum bath.¹² Ideally, one should consider a sufficiently large quantum bath to correctly model the effect of the semi-infinite leads. This necessarily introduces a large number of electronic degrees of freedom to the system. In terms of a TDDFT model, this induces many more wave functions and, more importantly, their representations at many more grid points into the computation. Specifically, the degrees of freedom have roughly an $O(N2)$ scaling, where N denotes the number of electrons. In a crude tight-binding approximation, a relatively larger number of electrons can be considered, but the computational challenge still remains.

Numerous attempts have been made to construct models so that only the electronic degrees of freedom near the junction are involved in the computation. These efforts range from the theoretical approach to derive open quantum system models¹³ to absorbing boundary conditions (ABCs),^14–16 the reduced density-matrix (RDM) approach by Subotnik et al.,¹⁷ the driven Liouville–von Neumann (DLvN) approach,^18–21 and the reduced-order approach.²² In principle, all these methodologies lead to a reduced quantum model with much fewer variables.

This paper focuses on such a reduction scheme for systems with Coulomb interactions. Coulomb interactions play a central role in electron transport and are responsible for the important observation of Coulomb blockade.^23,24 However, as suggested in Ref. 25 in the context of ABCs, it is still an open challenge to derive a boundary condition that takes into account such long-range effects. The ABCs can only be applied with a leap of faith. The development of the DLvN approach^18,19 and the reduced-order method²² also neglect this effect. However, in a more careful approximation of TDDFT, e.g., the real-space approximation^26,27 or the self-consistent tight-binding approximation,¹¹ one needs to take into account the Coulomb interactions, which have to be computed at each time step, and it leads to considerable computational overhead.

In order to compute the Coulomb potential at the junction, we replace a direct truncation of the bath regions with a sparse representation, where the grid spacings are gradually increased. The use of non-uniform grids that concentrate more at the regions of interest is not foreign in quantum mechanics: It has been a popular technique for computing the radial component of a Schrödinger equation, e.g., the logarithmic grids in the computation of pseudo-potentials.²⁸ Overall, our selected representations form a “multi-connected” sub-domain. Following the previously developed reduced-order technique,²² we project the Liouville–von Neumann equation into appropriate subspaces. As a result, the effective Hamiltonian consists of a self-energy that couples the sub-domains. The RDM, with a function interpolation, provides the global electron density profile, which can subsequently be used to estimate the Coulomb potential at the junction. The treatment using the domain decomposition method is reminiscent of the partition density-functional theory.²⁹ However, the current approach targets the dynamics at the center region, rather than the electron density in the entire domain.

The fact that the wave functions from density-functional theory are often extended creates a challenge for the model reduction effort based on the wave functions. Although such ideas have been discussed^2,4 to reduce the Kohn–Sham orbitals to a subdomain, it is not clear how to derive a reduced model that involves fewer orbitals. Therefore, we will base our reduction approach on the Liouville–von Neumann equation for the density matrix. A convenient and robust approach to derive a reduced model, as suggested by reduced-order techniques,^30–32 is by the Petrov–Galerkin projection, which reduces the method to finding appropriate subspaces. The reduced model for the RDM will be referred to as the reduced Liouville–von Neumann equation (RLvN). An important observation is that the Lie group structure is no longer present due to the fact that the effective Hamiltonian is non-Hermitian. Rather, the RLvN should be viewed as an open quantum system,³³ where the influence of the bath is implicitly incorporated. An important difference from those models, e.g., the Lindblad equation,³⁴ is that the present model does not have the trace-preserving property due to the electron transport nature.

In Sec. II, we will outline the domain decomposition and model reduction methods to derive the RLvN. Details regarding the implementations of the reduced model are presented in Secs. II E–II G. Results from several numerical experiments are demonstrated in Sec. III.

We formulate our method and algorithms based on the LvN equation

Here, ρ(x, x′, t) represents the density matrix of the entire system. Following TDDFT models,^2,35 the Hamiltonian operator H will depend on the electron density n(x, t), which can be determined from the density matrix ρ. These properties are expressed as follows:

Although our model works with the density matrix, our primary interest is in the electric current induced by a time-dependent external potential. The potential bias is denoted here by δU(t), and it is switched on at t = 0₊. It triggers the time evolution of electron density n(r, t), which has an implicit influence on the Hamiltonian H of the entire system. Motivated by the theory of linear response,^36–38 we first incorporate the external potential and write the total Hamiltonian as follows:

The response of the system due to the external potential in the linear response regime can be represented in terms of the perturbed density matrix,

Here, ρ₀ is the ground state density matrix of the system in the absence of the biased potential. A variety of methods are available for computing the ground state.^29,39 The corresponding electron density is denoted here by n₀(r), with its perturbation δn(r, t) given by

With the ground state density, we define an unperturbed Hamiltonian as

where V_XC is the exchange-correlation potential. V_H is the Hartree potential, and it is expressed as follows:

As a result, the total Hamiltonian H_tot(t) can be split up as follows:

where

Traditionally, the perturbation of the Hamiltonian operator is treated under the linear-response framework.⁴⁰ Also known as Sternheimer's formalism, it linearizes the Hamiltonian as follows:

Here, the integral operator $K$ is given in terms of the functional derivative with respect to the density,

Although Sternheimer’s approach is important in determining electro-magnetic properties, it leads to a dense matrix due to the Coulomb term. Therefore, we will keep the nonlinear form (9) instead. As we will demonstrate, it does not complicate the numerical implementation. We will break δH in (9) into three contributions,

Note that the Hartree term is linear with respect to the charge density,

Meanwhile, the exchange-correlation potential is nonlinear. This nonlinearity will be kept in our formulation.

After substituting (4) into (1) and dropping the high order term, [δH, δρ], the perturbed density matrix satisfies a response equation, given by

The perturbation approach also leads to a non-homogeneous term,

which incorporates the influence of the external potential as well as the density fluctuation δn.

To better elaborate the method, we first establish some notations. To begin with, we use Ω to denote the entire region. As a concrete example, we first discretize the Hamiltonian operator using a real-space method and the domain Ω consists of grid points. We denote the dimension by N = dim(Ω), and it indicates the dimension of the full system. In the case when a tight-binding method, e.g., Ref. 41, is used, the domain Ω consists of the atoms. The setup of the molecular junctions naturally divides the domain Ω into a center region, denoted here by $ΩC$⁠, and a bath region, denoted by $ΩB$⁠; $Ω=ΩB∪ΩC$ and $ΩB∩ΩC=∅$⁠.

As alluded to in the Introduction, the main obstacle in the model reduction is that the Coulomb potential is non-local, and completely neglecting the electron density in the bath region will lead to inaccurate computation of the Hartree potential. To address this issue, we extend the “center” region $ΩC$ by including points in the bath that are distributed according to

Here, x_R refers to the right end point of $ΩC$⁠. b > 1 indicates how the grid spacing is increased, and a > 0 is a prefactor. Similarly, we can select such grid points in the bath on the left, e.g., $xα=xL−ab−α,α∈N$⁠. Such grid distribution is known as the logarithmic grids in the solution of the Schrödinger equation in polar or spherical coordinates. Figure 1 demonstrates the distribution of such grid points, where the system consists of 300 lithium atoms. It also shows a scenario of how the grid points can be selected. This approach assumes that the system is one-dimensional, but it can also be extended to systems that are quasi-one-dimensional, e.g., a nanowire.

By including these additional grid points, we introduce the computational domain as the union of the selected grid points with the center region $ΩC$⁠,

The symbol $A$ indicates the labels of those grid points in the bath that have been selected according to (16). The dimension of the new domain Ω_I will be denoted by m. Of our particular interest is the case when m ≪ N. The new “bath” region corresponds to the complement of Ω_I,

Motivated by its remarkable success of reduced-order modeling techniques^30–32 in handling large-dimensional dynamical systems, we apply the Petrov–Galerkin projection approach directly to the full system equation (14). The method uses two subspaces $V$ and $W$ as the solution and test spaces, respectively. Although conventional Petrov–Galerkin projections arise from approximations of PDEs,⁴² the technique can also be applied to discrete dynamics, e.g., classical particle systems.^43,44 For reduced-order problems, the approach is more robust than the traditional moment matching methods.³² To facilitate implementations, the space $V$ will be represented by a set of basis vectors, which are arranged into the columns of a matrix, denoted by V. Therefore, $V=Range(V)$⁠, and the same notations will be used for the space $W$ as well, $W=Range(W)$⁠. The dimension of the two subspaces is assumed to be m.

Overall, the Petrov–Galerkin projection for Eq. (14) consists of the following three steps:

1. Seek the solution in the subspace $V$⁠. In our case, the approximate density matrix is spanned by the matrix V as follows:

   $δρ(t)≈δρ̃(t)=VD(t)V*.$

   (19)

   In the subspace approximation, the m × m matrix D will be referred to as the nodal values.
2. Project the dynamics (14) to the orthogonal complement of a subspace $W$⁠. We choose the subspace spanned by the column vectors of a matrix W, which yields

   $iddtW*δρ̃(t)W=W*[H0,δρ̃(t)]+Θ(t)W.$

   (20)

   The first two steps offer a low-rank approximation of the density matrix. Such projections have been widely used in solving matrix equations.⁴⁵ 

3. Use quadrature formulas to approximate functions outside the domain Ω_I using the nodal values in D(t).

The matrix products in the first two steps of the Petrov–Galerkin projection give rise to two coefficient matrices,

In the Galerkin projection, these matrices are often referred to as the mass and stiffness matrices.⁴² Direct computations from the first step (19) and second step (20) yield

where H_eff is the reduced Hamiltonian,

The bracket (22) has been generalized to non-Hermitian matrices,

Furthermore, the non-homogeneous term Θ_eff(t) is a reduction of the non-homogeneous term Θ(t) in (14),

This matrix product will be treated in step (3), i.e., approximated by a quadrature formula.

We will choose the subspace $V$ by specifying a set of orthonormal basis following our previous work.²² To be more specific, using the partition Ω = Ω_I ∪ Ω_II, we express the basis as columns of a matrix V,

Here, we follow the notations in domain decomposition methods⁴⁶ and order the grid points so that the points in Ω_I appear first. The elements in the Hamiltonian matrix H can be arranged in the same way to a block matrix,

Thus, in (26), we essentially pick the standard basis in the subdomain Ω_I. Then, they are extended to Ω_II by filling in zeros. This choice is motivated by the reduced-order modeling.³² See Ref. 22 for details about how the density-matrix equation (14) can be put into a standard form of reduced-order problems. In general, the choice of V corresponds to the identification of quantities of interest by defining the output vector.³² As we will demonstrate later, it will single out the components of the density matrix in the subdomain Ω_I, leading to a description using the reduced density matrix (RDM). In turn, one may consider how the rest of the dynamics is controlled by the RDM. The vectors in W represent such coupling. Specifically, we choose

Here, we require Im(ɛ) < 0 for stability.²² 

With the choices (26) and (28), one finds that

Here, G_I,I refers to the first diagonal block of Green’s function,^47,48G(ɛ). From the block inversion formula, we obtain

Meanwhile, for the stiffness matrix, we have

Combining with (23) and (30) and noticing the orthogonality V * V = I, we get

where Σ_I,I is the self-energy,⁴⁹ given by

Therefore, the choices of the subspaces according to (26) and (28) have a remarkable connection with the well-known Green’s function approach.^47,50 In particular, the self-energy Σ_I,I describes the influence from the domain Ω_II. It is also worthwhile to point out that since Ω_I is multi-connected, as depicted in Fig. 1, the self-energy Σ_I,I is an operator that couples a domain with many fragments.

The selection of the domain Ω_I is with the purpose to take into account the density fluctuation δn in the bath. With the electron density in the domain Ω_I, we can reconstruct the correction to the Hartree part of the Hamiltonian (7) using a quadrature formula,

Here, $Q1$ denotes the set of quadrature points and w(x, x′) denotes the quadrature weights. In the numerical tests, we will use the three-point Gaussian quadrature for each interval in Ω_I.

Since x_α from the bath is included in the computational domain Ω_I, the density matrix at x_α is explicitly computed, and as its diagonal elements, δn(x_α) is accessible. In the computation, we use a spline interpolation based on the electron density at these points, which is then evaluated at those quadrature points $x∈Q1$ so that the quadrature formula (34) can be applied to update the Hartree potential.

We now discuss how the non-homogeneous term (25) is evaluated. Due to the dependence of δH on the electron density (12), this term has to be computed at each time step. The main issue here is that the matrix product with W involves a summation over the entire domain, which requires $O(m2N2)$ operations.

As with many other applications of Galerkin projection,^42,51 this term can be treated by quadrature approximations. For convenience, let us first work with a summation form. Namely, we consider

In the present setting, x_k’s are uniformly distributed grid points in the entire domain. We assume that f(x) has certain smoothness (e.g., twice continuously differentiable). In the quadrature approximation, the sum is approximated by

Here, $Q2$ denotes a subset of the points, which will be regarded as quadrature points and A_α will be the corresponding weights. The quadrature weights in our computation are determined in advance by a spline interpolation. A simple idea is the following: For each x_α, one can set f(x_α) = 1 but zero for all other quadrature points. Using a cubic spline interpolation, f(x) ≈ s(x), and from the quadrature formula above, one can obtain the quadrature weight A_α = I[s], i.e., directly summing up the spline interpolant s(x).

Returning to (25), with the quadrature formula, we can reduce the matrix product,

The overline denotes the complex conjugate. If the number of quadrature points is k, this procedure reduces the computation from $O(m2N2)$ to $O(m2k2)$⁠.

An interesting case is when $Q2=ΩI$⁠. Combining (25) and (37), we have

where S is a skew-Hermitian matrix with elements given by

In the last step, we have used Eq. (21).

In Secs. II D–II F, we explained how the terms in the RLvN equation (22) are computed at each time step, and they are summarized in Algorithm 1. Here, we discuss how the dynamics is integrated.

We first discuss the approximation of the full model (1). We assume that a finite-difference approximation has been applied so that the Hamiltonian operator, as well as the charge density, is represented at grid points. We adopt the operator-splitting scheme of Watanabe and Tsukada⁵² for TDDFT models. More specifically, we split the Hamiltonian as follows:

where the potential term may include Coulomb, exchange-correlation, and external biased potential.

Despite the non-linearity in (1), one can still express the solution with an evolution operator,

The evolution operator can be effectively approximated by an operator-splitting scheme,

Specifically, the operator $E0$ corresponds to the dynamics determined by the kinetic energy,

for which the evolution operator corresponds to a unitary dynamics,

The matrix exponential in U₀ can be computed from direct spectral decomposition, e.g., by using fast Fourier transform.

On the other hand, the operator $E1$ dictates the dynamics

The diagonal of the right hand side is zero, which implies that the electron density remains constant. Therefore, this is effectively a linear equation. As a result, we can write

In addition, since H₂ is only diagonal, the evaluation is straightforward.

We now turn to the RLvN (22). Using the variation-of-constant formula, we first obtain an integral form of the solution (with the spatial variable x and x′ suppressed),

Here, the matrix U(Δt) is given by

which only needs to be computed once. Since the term Θ_eff in Eq. (25) also depends on the diagonals of D(t), this remains as an integral equation. To obtain an approximation, one can approximate the integral using the values at t (τ = 0). This yields

which is an analog of the exponential-Euler method.⁵³ 

This method can be extended to form a class of integrators that resemble the exponential-Runge–Kutta methods.⁵⁴ As an example, one can apply the above formula for half of a step, which is then combined with the mid-point rule for the integral (47),

When the numerical stability is largely determined by the spectrum of H_eff, such exponential integrators often allow much larger step sizes.

As examples, we consider a 1d lithium (Li) system, consisting of three Li chains with a lattice spacing of 2.8 Å. Two atom chains, each with 148 atoms, are used to model the bath. They are separated from the middle chain with four Li atoms by adding an extra spacing of 2.8 Å. This is the largest system that we have studied for which the full system can be simulated within a realistic time frame to generate a reference solution. Similar systems have been used by Refs. 55 and 56 to study transport properties.

The units follow those in Ref. 57 using Bohr radius and eV as the length and energy units, which subsequently determine the time units. We use a real-space method to approximate the solution of the full LvN model (1) as well as the reduced model that follows. Specifically, we use a finite-difference spatial discretization with N = 4000 grid points and grid spacing Δx = 0.2114. We use the three-point finite-difference formula for the kinetic energy term. The electron current can be extracted from the off-diagonal entries of the density matrix,

The initial ground state is computed self-consistently using a mixing method^58,59 with four previous iterations at each step. The step size for the time integration is Δt = 0.0025, and Eqs. (1) and (22) are integrated up to t = 50. Due to the finite size of the bath, the boundary effect will occur. However, we have verified that such a finite-size effect will not affect the dynamics in the center region within this time period.

We follow the model by Baker et al.⁵⁷ for the Coulomb potential and the exchange and correlation potentials. Figure 2 shows the exponential form from Ref. 57 for the Coulomb potential, plotted together with the grid points in the real-space approximation to show the scale. Also shown in the figure is the regularized Coulomb potential (Yukawa),

In terms of the atomic spacing, the exponential model has an effective range up to second nearest neighbors. However, the regularized Coulomb potential extends much farther. Both models will be implemented in the computation.

For the reduced model (22), we first pick the subdomains as follows. The entire domain is positioned in the interval [0, L] with L = 453.6. We keep the center region in the interval [x_L, x_R] = [205.8, 247.8]. Starting from this center region, we gradually increase the grid size. For example, in the bath on the right, we define the logarithmic grid x_α = x_R + ⌊2 × 1.2^α⌋Δx, with ⌊·⌋ denoting the nearest integer. Out of the 4000 grid points, m = 462 points have been selected with this procedure. Three-point Gaussian quadrature formulas are used for the calculation of the Coulomb potential in (34). For the non-homogeneous term (25), since we are approximating a summation, we need to pick quadrature points that coincide with the grid points. In the center region, we already keep all the grid points. Meanwhile, in the bath, for each interval with x_α+1 − x_α > 9Δx, we introduce two additional quadrature points in that interval.

Using the operator-splitting method (42), we first run the full model (1) up to t = 50. We exert a potential bias U_L = 0.1 on the left bath and U_R = 0 in the other bath. We first show the charge density at t = 8, along with the charge fluctuations δn(x, t) around the center region in Fig. 3. We observe that the charge density n(r, t) exhibits oscillations around the nuclei in the bath, which extends to the entire domain. The two “dips” in the density profile occur at the junctions. In contrast, the charge perturbation δn(x, t) spreads out into the bath, and it exhibits decay toward the far side of the bath. Therefore, it is important to take into account the density perturbation in the bath to accurately compute the Coulomb potential in the device. The reduced model constructed in this paper works with this charge perturbation, which is represented at sparse grid points with increasing spacing.

We first examine the dynamics of the current predicted by the full model (1). As a comparison, we consider three cases: the exponential potential,⁵⁷ the regularized Coulomb (52), and the case where the Coulomb potential is turned off. Figure 4 shows several snapshots of the current around the device region. In addition, Fig. 5 shows the current in the middle of the device, x = 424.6465. Triggered by the potential bias, a current peak first develops at the left junction, which later spreads to the left bath and the device region. When the Coulomb potential is turned off, the current in the device region quickly develops into a plateau, reaching a steady state. However, with the Coulomb potential, the situation is more subtle. Within the time period [0, 50], we have not observed any steady state. At t = 50, with the regularized Coulomb potential (52), the current starts to flatten. However, with the exponential model, the current in the device is still uneven at this point.

Next, we implement the algorithm (50) and run the reduced model (22) using the regularized Coulomb potential (52). For the subspace $W$ (28), we pick the shift to be ɛ = μ − 5i, with μ being the Fermi level. Figure 6 depicts the electron current around the junction for various time instances. One can observe that the potential bias creates a current peak initially, which then develops into an oscillatory pattern. Subsequently, the current profile starts to widen and spread toward the bath region. The time history of the current at x = L/2 is shown in Fig. 7. In both cases, the results agree well with the full simulation results within the transient period [0, 15]. Over longer time, the approximate current shows a tendency to approach to the steady state.

We also implemented the reduced model on a larger system with 480 atoms and N = 6000 grid points. In this case, the direct simulation of the full model is projected to take over a month, and therefore, we did not make such a comparison. Instead, we compare to the result from the previous test with 300 atoms. The computed current is shown in Fig. 8. We observe that these calculations yield consistent results, which also indicates that the effect of the bath is adequately represented by 148 atoms within this time interval.

The numerical tests are repeated using the exponential Coulomb potential.⁵⁷ Compared to the previous case, the current profile exhibits more oscillations, as shown in Fig. 9. In this case, the reduced models still accurately predicted the current in the center region up to t = 15. The comparison of the current at x = L/2 is shown in Fig. 10. We observe that although the reduced model is quite accurate within the transient period [0, 15] and is also reasonable in final period [30, 50], it has significant error in the period [15, 30].

It is also worthwhile to point out that the reduced-order approach^30,32 is effective in predicting pre-selected quantities of interest (the output), which in our case is the density matrix in the region Ω_I, and it embodies the electron current in the center region. However, it does not involve the current in the bath regions. Therefore, the reduced model (22) will not be able to capture the current outside the center region.

Next, we monitored the central processing unit (CPU) time for solving the full LvN equation (1) and the RLvN (22). We run the models for systems of various sizes. To understand the reduction in computational cost, we studied two distinct cases. First, we hold the system fixed while reducing the mesh size Δx, which causes the dimension N (number of grid points) to increase. The results are summarized in Table I. As the dimension N (number of grid points) increases, the CPU time corresponding to the simulation of the full LvN equation (1) increases considerably. This can be understood from algorithm (44): each step consists of multiplications by the unitary matrices with dimension N × N, which contributes to the majority of the computational cost. On the other hand, for the RLvN (22), algorithm (50) involves the multiplication by U(Δt) (48), which is of dimension m × m, and it offers considerable reduction in this aspect. A close inspection reveals that, aside from the matrix multiplication, the main computation comes from the calculation of the Coulomb potential (34) and the non-homogeneous term (37); for the RLvN (22), both have a $O(m2)$ scaling due to the quadrature approximation.

Another important scenario is to hold the mesh size Δx while enlarging the bath size. This will also increase the size of the full system and hence cause the computational cost to grow, as can be seen in Table II. However, for the RLvN (22), since most of the grid points are around the device region, the dimension of the subspace increases very little. The CPU time for the RLvN (22) also remains remarkably low. Note that when N = 4800, although the dimension of the subspace is less than the previous setup, the total CPU time is still comparable. This is due to the quadrature approximation.

Finally, we check the occupation numbers by examining the eigenvalues of the density matrix obtained from the reduced model (22). Since the reduced model no longer follows a unitary dynamics, the occupation numbers associated with density matrix are not guaranteed to be between 0 and 1. Therefore, we resort to numerical tests to examine this property. As shown in Fig. 11, most of the occupation numbers lie between 0 and 1. About five of these are slightly larger than 1 with the largest being 1.0210. Therefore, for the reduced model, this property holds approximately.

In this paper, we have constructed a class of reduced-order models for simulating electron transport problems. Based on a decomposition of the physical domain and by partitioning the density matrix accordingly into different blocks, we regard such a dimension reduction problem as a reduce-order problem,²² where a subspace projection method can accurately predict certain quantities of interest.^31,32 This approach naturally involves the self-energy in the effective Hamiltonian. In this work, our focus is placed on the treatment of long-range interactions. Specifically, we choose subspaces so that the electron density in the bath is captured on some sparse grids, which, in turn, determines the Hartree potential. This enables the accurate computation of the Coulomb and exchange-correlation potentials at the molecular junction. It is possible to enlarge the subspaces to enable more accurate approximations, as demonstrated in our prior work²² for a tight-binding model. However, based on the current results, we observed that the present approach using subspaces (26) and (28) is quite easy to implement, and it has already shown good agreement with the full model. Since our formulation is primarily through algebraic manipulations, the method is not specific to the underlying geometry of the electronic system. So we expect that this method can be extended to high-dimensional systems, e.g., that in Ref. 60, with similar selection of the contact region and grid points.

X.L.’s research was supported by the NSF under Grant Nos. DMS-1819011 and DMS-1953120.

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