Electron Dynamics in Open Quantum Systems: The Driven Liouville-von Neumann Methodology within Time-Dependent Density Functional Theory

A first-principles approach to describe electron dynamics in open quantum systems driven far from equilibrium via external time-dependent stimuli is introduced. Within this approach, the driven Liouville-von Neumann methodology is used to impose open boundary conditions on finite model systems whose dynamics is described using time-dependent density functional theory. As a proof of concept, the developed methodology is applied to simple spin-compensated model systems, including a hydrogen chain and a graphitic molecular junction. Good agreement between steady-state total currents obtained via direct propagation and those obtained from the self-consistent solution of the corresponding Sylvester equation indicates the validity of the implementation. The capability of the new computational approach to analyze, from first principles, non-equilibrium dynamics of open quantum systems in terms of temporally and spatially resolved current densities is demonstrated. Future extensions of the approach toward the description of dynamical magnetization and decoherence effects are briefly discussed.

To model electron dynamics in molecular junctions, theoretical methods have been developed, 16,24, which can be broadly divided into methods that use model Hamiltonians (usually formulated in energy space) treating general transport phenomena while circumventing detailed descriptions of specific junctions, 63,64,72-78 and methods that explicitly take into account the chemical composition and structure of the system and thus allowing for a direct comparison with experimental findings. 31,32,56,60,61,71,[79][80][81][82][83][84][85] The latter are often highly computationally demanding and thus limited to treating relatively small systems.

The Driven Liouville-von Neumann Formalism
The DLvN formalism relies on a unitary transformation from a finite real-space representation of the junction model to its spectral state-representation, 65 The general DLvN scheme can be divided into four steps: (i) Spatial partitioning: the molecular junction is represented by a fully atomistic finite model system that is formally partitioned into three sections: the left lead (L), the extended molecule (EM), and the right lead (R, see illustration in Figure 1). The EM includes the active molecule and its adjacent lead subsections that serve to buffer the molecule from the lead regions, where boundary conditions are applied. In a non-orthogonal, atom-centered basis-set representation of the Kohn-Sham (KS) molecular orbitals, the partitioned KS Hamiltonian and overlap matrices ( and , respectively) can then be written in the following block form: For simplicity, in what follows we assume that the lead sections are spatially well separated such that the , and , blocks (and their conjugate counterparts) are negligible and can be safely replaced by zero matrix blocks of appropriate dimensions.
(ii) Block orthogonalization: when a non-orthogonal basis-set representation is used, one must ensure that the boundary conditions applied at the far edges of the systems do not directly affect the dynamics of the extended molecule region. 69 To this end, the block orthogonalization procedure of Kwok et al. 108 is adopted to transform the localized basis-functions of the EM section into new EM basis-functions that are mutually orthogonal to those of the finite L/R lead models. This block orthogonalization procedure where and are null and unit matrices of the relevant dimensions, respectively. Following this transformation, the overlap matrix becomes block diagonal: Since the transformation leaves the diagonal lead blocks, , , and unaffected, the procedure for applying boundary conditions does not change. We note that if an orthogonal basis-set is used, this step can be skipped.
where = , , are the unitary matrix blocks that transform the generalized eigenvalue equations: to their diagonal representation, where and are the generalized eigenvalues and eigenvectors matrices of each block, respectively. Within this representation, ̃= †̃ is a diagonal matrix containing the eigenvalues of ̃, and ̃= †̃= . With this, the full overlap matrix becomes the identity: and the KS Hamiltonian adopts the following form: ̃= †̃= (̃̃, ,̃̃, where the off-diagonal ̃, blocks represent couplings between the eigenstates of sections and , and the lead-lead couplings remain zero. Following the site-to-state transformation, the atomistic representation of the junction is replaced by a state-representation, where the single-particle states of the EM section are coupled to the corresponding lead states (see Figure 2).
Here, ̃ is the single-particle density matrix of the entire system, given in the state-representation. Furthermore, approaches based on field-induced polarized boundary conditions are limited in practice to linear two-lead setups, where a uniform field is applied along the main axis of the junction model. Therefore, we opt to pursue a full-fledged implementation of the DLvN scheme, where at each time-step the following workflow is followed: (a) obtain the KS Hamiltonian, , in a general non-orthogonal atomic orbital basis.
(b) transform to the state-representation to obtain ̃.
(c) build the target density matrices ̃0 and ̃0 and construct the driving term.
(d) transform the driving term to the site representation.
Note that to avoid the implicit time-dependence of the ̃ transformation (via ̃[̃] ), the propagation of the density matrix in item (e) is performed in real-space (see SI sec. 1 for further details.). To this end, we use an implicit Euler propagation scheme (see SI sec. 3), where: ( + ) = ( ) + ⋅ ( + , ( + )), and ( + , ( + )) represents the right-hand-side of Eq. (11) in the site representation. To solve Eq. (12), an iterative fixed-point scheme is implemented at each time step, where ( + , ( + )) at a given iteration is evaluated using ( + ) of the previous iteration, keeping ( ) on the right- describes the dynamics of ̇ and not of ̃̇ which, in turn, involves the unknown explicit dynamics of the transformation matrix (see SI sec. 1 for further details). These issues are absent in the block diagonal representation, which we use to perform the particle current calculation (see SI sec. 4): where [. ] is the trace operation, ℑ{. } represents the imaginary part.
For steady-state current evaluations, the right-hand-side of Eq.
Here, ̃ is the steady-state density matrix in the state representation, are the projection matrices onto the left and the right lead states, respectively, and Eq. (14) is solved iteratively using the following scheme. Given the density matrix at iteration , ̃, the KS Hamiltonian is constructed, and the Sylvester equation is solved to give ̃+ 1 . A damping scheme admixing these two density matrices, with weights and (1 − ), respectively, is then used to construct the next step KS Hamiltonian matrix, ̃. At each iteration step, ̃ is transformed to the block diagonal basis, ̃, and the particle current is calculated using Eq. (13) The process is repeated until the steady-state current is converged, such that its RMS value during = 20 consecutive iterations is smaller than a preset value (chosen as 10 −10 a.u. for the hydrogen chain calculations).
The entire simulation scheme is implemented in Python, 111 which makes recurrent calls at each timestep to the Gaussian suite of programs 112 to evaluate the KS Hamiltonian (and overlap) matrix elements in the atomic orbital basis. We note that the current proof-of-concept implementation employs spincompensated electron densities.

Results
The developed methodology is first benchmarked using a simple linear hydrogen chain molecular  (13) for several bias voltages. To improve numerical stability, the simulation starts from the ground state density of the system and the bias voltage is turned on gradually using a hyperbolic tangent switching function (see SI sec. 6). We confirm that during the propagation, all state occupations remain bound to The excellent correspondence between the steady-state currents obtained using the dynamical and the Sylvester calculations indicates that the DLvN EOM indeed reached a stable stationary state. Going beyond the total current flowing through the system, our approach allows also to analyze the spatially resolved current density within the chain: where ( ) is a matrix defined in the atomic orbital basis as ( ) = ℏ ( ) ( ) (see SI sec. 4). Figure 4 illustrates the spatially resolved steady-state current density for the hydrogen chain junction, whose total current is shown in Figure 3, under a bias voltage of 0.3 V obtained using Eq. (17). The axial (z) component of the steady-state current density is mostly uniform along the EM section with some expected variations near the nuclei positions (see sec. 7 of the SI for a plot of the current density integrated over the xy plane along the EM section). In the seamline between the weakly coupled molecule and the lead sections of the EM region (vertical black dashed lines) the current density reduces, indicating that it spreads over a larger cross section. The driven lead sections exhibit a less uniform current density map. This can be attributed to the fact that the uniform driving rate, Γ, applied in the state-representation of the junction translates to a spatially varying driving rate in real-space, with a larger value near the lead/EM interface region, where electrons are absorbed or injected. As a consequence, near the interface region of the sink lead, the direction of the current is reversed (see streamlines on the right side of Figure 4). This artificial behavior in the (unphysical) driven lead regions, however, has minor influence on the EM currents and negligible effect on the current flowing through the molecule itself. Far away from the interface, the current in the leads decays indicating that the lead approaches the equilibrium state of the implicit bath to which it is coupled. Following the benchmark hydrogen chain calculations, we now turn to discuss a more realistic model junction consisting of two graphene nanoislands bridged by a benzene molecule (see Figure 1b, junction model coordinates can be found in SI section 9). The geometry of the system was optimized via the LAMMPS 120 software using the reactive empirical bond-order potential (REBO). 121,122 The energy minimization was performed using the Fast Inertial Relaxation Engine (FIRE) algorithm 123 with a force tolerance of 10 −6 eV/Å. The vertical coordinates of all carbon atoms were fixed to keep the structure planar, so as to mimic a substrate supported junction model. The EM section is chosen to include the benzene unit and its two adjacent CH2 groups. The total current was calculated using Eq. (13) with the PBE 113,114 functional and the STO-3G and 6-31G** atomic centered Gaussian type basis-sets for the lead and EM sections, respectively. 115 A driving rate of ℏΓ = 1.09 eV was used to yield a smooth lead section DOS (see SI sec. 6). Figure 5 shows the current dynamics flowing through the EM section for several bias voltages calculated using the DLvN EOM, as well as the corresponding steady-state currents obtained from the solution of the Sylvester equation (14). The good agreement between the DLvN and Sylvester steady-state currents indicates that the DLvN EOM has reached a true stationary state of the system.

Conclusions
The

The driven Liouville von Neumann equation of motion in the realm of time-dependent density functional theory
In Here, ℋ is the KS Hamiltonian operator and = √−1. Next, we span the KS orbitals within a localized basis-set representation {| ⟩}: where is the th expansion coefficient of KS orbital | ⟩. Plugging Eq. (S2) into Eq. (S1) and assuming that the basis orbitals are constant in time we obtain: Multiplying Eq. (S3) by ⟨ | we obtain: Defining the overlap and KS Hamiltonian matrix elements as ≡ ⟨ | ⟩ and ℋ KS ≡ ⟨ |ℋ KS | ⟩, respectively, Eq. (S4) becomes: Since this equation is valid for all values of the indices and it can be written in matrix form as: Multiplying by the inverse of the overlap matrix, −1 , on the left we obtain: Accordingly, one can write the EOM for the complex transpose coefficient matrix as follows: where we used the relation ( −1 ) † = ( † ) −1 ( = † = ( −1 ) † = ( −1 ) † † ⟹ ( −1 ) † = ( † ) − ) and the fact that the overlap and Kohn-Sham matrices are Hermitian, such that † = and = † .
The latter relation stems from the fact that the density matrix, upon which depends, is Hermitian by construction (see Eq. (S9) below) and so are all the operators within (kinetic energy, Hartree, exchange, correlation, and external potential).
We can now define the single-particle density matrix in the localized basis-set representation as: where is a diagonal matrix holding the occupation numbers of the different single-particle states on its diagonal. The time evolution of the density matrix is obtained by its time derivative: ̇=̇ † +̇ † +̇ †. (S10) Here, the first two terms on the right-hand-side correspond to pure orbital dynamics, whereas the third term represents the dynamics of the orbital occupations. Inserting Eqs. (S7) and (S8) into Eq. (S10) we obtain: To rationalize this choice, we now plug Eq.
(S12) in Eq. (S11) to obtain: This equation assumes the form of a Liouville-von Neumann equation for a microcanonical (or canonical) system but with a general Hamiltonian matrix, + , that is neither Hermitian nor anti-Hermitian. The latter can be viewed as a dressed Hamiltonian, where we identify as a selfenergy like term representing the effects of the coupling of the system to an implicit bath. Note, however, that is energy independent and hence should be viewed as an approximation of the self-energy within the wide band limit.
To obtain the explicit expression of within the DLvN EOM we divide the system into three sections comprising of the left lead, the (extended-)molecule, and the right lead. then serves to mimic the effect of coupling of the lead sections to implicit Fermionic baths, characterized by equilibrium Fermi-Dirac distributions with given chemical potentials and electronic temperatures. To this end, we first neglect , which is equivalent to neglecting the real-part of the implicit baths' self-energies that induce lead level shifts due to the lead/implicit-bath couplings. This approximation becomes valid for sufficiently large lead models, with a relatively uniform and dense manifold of states, such that the level shifts become small with respect to the inter-level spacing. The remaining imaginary part, , marked for brevity as , introduces a finite lifetime (broadening) to the various lead levels due to their coupling to the implicit single-particle states of the reservoir. Hence, within the DLvN approach, the dressed KS Hamiltonian acquires the form: In this form, can be identified as an imaginary absorbing potential added to the lead sections of the original KS system serving to absorb outgoing electrons near the system boundaries (thus preventing their back-reflection into the system). Naturally, in order to avoid complete electronic depletion of the system, complementary emitting potentials should also be introduced in order to inject thermalized electrons into the system, as shown below.
Using the dressed Hamiltonian form of Eq. (S14) in Eq. (S13) we obtain: We can now multiply Eq. (S15) by from left and from right to obtain: Next, we introduce a block diagonalization transformation, , into Eq. (S16) to nullify the off-diagonal overlap matrix blocks as follows: 1 Eq. (S17) can be rewritten as follows: Next, we multiply Eq. (S20) by † on the left and by on the right, to obtain: Introducing −1 and ( † ) −1 † in the last two terms, respectively, yields: Next, we define: We note that since is a fixed transformation (time-independent within the fixed nuclei Born-Oppenheimer approximation) the relation ̃̇=̇ holds. With these definitions we obtain: where we have used the fact that ̃ † = ( † ) † = † † .
Next, we introduce the site-to-state transformation: such that ̃= †̃ is diagonal and †̃= are unit submatrices of the appropriate dimensions. With this, Eq. (S24) can be rewritten as: Using the relation †̃= we get: Next, we define: to obtain: Multiplying by † on the left and on the right we arrive at: Defining: We finally obtain: where ̃ † = ( †̃) † = †̃ † . In its simplest form ̃ is written as: which represents uniform broadening of all left and right lead levels. Hence, the last two terms in Eq.
(S32) can be written as: (S36) as follows: . Hence, we obtain † † =̃− 1 † = −1 such that: This is the DLvN EOM-our working equation-given in the site representation, which we use for the time propagation.
To further simplify this expression and avoid the full matrix transformations appearing in the driving term of Eq. (S39) we may separate it into the sink and source terms and treat them separately. Using Eqs. (S18), (S25), and (S34) the sink term can be written as: Using Eqs. (S23) and (S28) we now obtain: We may now use the relation ( ) = − ( ) to write: Using Eq. (S25) for the transformation matrix we can write: Therefore, we have where we have used Eq. (S18) for the expressions of and −1 . Similarly, and therefore: where ( , ) = [ ( − )/ + 1] −1 is the Fermi Dirac distribution, is Boltmann's constant, is the electronic temperature, / is the ℎ eigenvalue of the / lead, and / is the corresponding chemical potential. We note that in the present implementation we resort to Eq. (S39) for the propagation without considering the above simplifications. In practice, the propagation is performed as follows: 1. Construct a junction model with predefined lead and extended molecule sections.
2. Perform a ground state calculation to obtain the overlap matrix, , and the initial and matrices in the site representation.

Driving rate sensitivity test
While the driving rate, Γ, appearing in Eq. (11) of the main text can, in principle, be determined from the self-energy of the semi-infinite lead models, 2 in the current implementation we use it as a free parameter. To determine the value to be used in the dynamical simulations, we broaden the discrete energy levels of the finite lead models with Lorentzian functions of different widths and adopt the Lorentzian width parameter that provides a density of states that represents well that of a semi-infinite system (not too narrow to provide a discrete spectrum and not too wide to artificially wash out the electronic structure features of the lead) as our Γ value for the time-dependent calculations. [3][4][5] Figure S1 compares the density of states of the hydrogen chain studied in Fig. 3 of the main text for several Lorentzian broadening widths, from which we select the value of ℏΓ = 0.61 eV for the dynamical calculations performed for the hydrogen chain junction in the main text. To verify that our results are relatively insensitive to this choice, we present in Figure S2: the steady-state current as a function of Γ showing that above a value of ~0.3 eV the steady-state current weakly depends on Γ within the relevant range, set by the density of states analysis discussed above. Figure S1: The artificially broadened density of states of the hydrogen chain junction model considered in Fig. 3 of the main text for four Lorentzian widths. The inset shows the full broadened spectra. The broadened lead density of states for the graphitic junction model depicted in Fig. 1b of the main text appears in Fig. S5. The adopted value of ℏΓ = 1.09 eV provides adequate broadening of the energy levels to results in a density of states that satisfactorily represents that of the infinite nanoribbon lead electronic structure (Fig. S6).

Propagators
The driven Liouville von Neumann equation of motion for the single-particle density matrix (Eq. (S39) above) could, in principle, be propagated using one of many available propagation schemes. [6][7][8][9][10] However, in this particular case, care must be taken, since

Total current calculation
Within We are interested in calculating the instantaneous particle current flowing between the and driven lead sections through the region. This can be obtained from the expression for the time derivative of the particle number in this region, ̇. To this end, the relation = ( ) can be used, where is the total number of electrons, is the overlap matrix, and is the single particle density matrix. Note where we used the fact that is Hermitian and is real and symmetric. We therefore see that ( ) is not necessarily real and therefore cannot represent the particle number in the section. This can be remedied by using the Löwdin symmetric version of the trace formula: whose partial trace is real: [ ( Since the block diagonalization transformation only rotates the basis to make it diagonal to the and bases without modifying the latter (while readjusting the / and / Kohn-Sham Hamiltonian coupling blocks), the partial sums over the , , and indices retain their spatial interpretation as belonging to the corresponding system sections. We can, therefore, write the full trace as the sum of partial traces over the separate system sections: (S62) We shall first transform the driving term on the right-hand side from the state-to the block diagonal representation. To this end, we rewrite it in the following form: The first term on the right-hand-side in Eq. (S63) reads: .
Similarly, the second term on the right-hand-side reads: . (S71) Substituting this in Eq. (S60) for the time derivative of the particle number in the section yields: .
The overlap matrix in the block diagonal basis assumes the following form: , whose block is: The second term reads: , whose block is: The fourth term reads: which is the final expression that we use for the evaluation of the current in the block-diagonal representation.

Effects of the basis sets and size of leads on the calculated current
In all hydrogen chain benchmark calculations presented in the main text we used 180 atom lead models in combination with the STO-3G atomic centered basis-set to represent the Kohn-Sham orbitals and the PBE functional approximation. To evaluate the sensitivity of our results towards these choices we present in Figure S6a the steady-state current versus the applied bias voltage for lead models of 120,

Driving rate switching function
To increase numerical stability of our graphene nanoribbon junction transport calculations, the driving rate, Γ, was gradually increased from zero to its full value, Γ . 12 To that end, we chose a hyperbolic tangent switching function of the form: where is the time, 0 = 0.05 fs, and = 0.36 fs.

Integrated current density along the extended molecule section
As a further consistency check of our current calculations, we compare in Figure S8 the steady-state current of the hydrogen chain model (Fig. 3 of the main text) evaluated from the partial trace of the partitioned single-particle density matrix (Eq. 13 of the main text) to that obtained by cross-sectional spatial integration of the current density (Eq. 18 of the main text). In the figure we plot the latter evaluated at several axial positions along the chain. As expected, the steady-state current is spatially uniform, increasing/decreasing only near the edges of the EM section, where the source/sink terms are applied. Furthermore, the integrated current matches well the one obtained via Eq. 13 of the main text (dashed orange line), further validating our methodology.

Total number of electrons for the case of fractional occupations
In Eq. (S60) of SI section 4 above, we used the expression = ( ⋅ ) for the particle number. Here, we demonstrate that this expression is valid also for the case of non-idempotent density matrices, representing fractionally occupied states, of the following form: Here, is the total number of basis functions, and are the expansion coefficients of the orthonormal Kohn-Sham orbitals, ( ), within the non-orthogonal atom-centered orbital basis { }: and 0 ≤ ≤ 1 are the occupation numbers. For brevity of the presentation, spin indices have been omitted herein.
The density at a point in space can be written as: Spatially integrating the density ( ), we obtain the total number of electrons: where the overlap matrix elements in the atom-centered orbital basis are given by = ∫ ( ) * ( ) 3 . Using Eq. (S99) equation (S102), we finally obtain: which has the same form as the standard expression for the total number of electrons of a closed ground state system with integer state occupations.

Cartesian atomic coordinates of the GNR/Benzene/GNR junction
All atomic coordinates are given in Å.