Driven Liouville-von Neumann Equation for Quantum Transport and Multiple-Probe Green's Functions

The so called Driven Liouville-von Neumann equation is a dynamical formulation to simulate a voltage bias across a molecular system and to model a time-dependent current in a grand-canonical framework. This approach introduces a damping term in the equation of motion that drives the charge to a reference, out of equilibrium density. Originally proposed by Horsfield and co-workers, further work on this scheme has led to different coexisting versions of this equation. On the other hand, the multiple-probe scheme devised by Todorov and collaborators, known as the hairy-probes method, is a formal treatment based on Green's functions that allows to fix the electrochemical potentials in two regions of an open quantum system. In this article, the equations of motion of the hairy probes formalism are rewritten to show that, under certain conditions, they can assume the same algebraic structure as the Driven Liouville-von Neumann equation in the form proposed by Morzan et al. [J. Chem. Phys. 2017, 146, 044110]. In this way, a new formal ground is provided for the latter, identifying the origin of every term. The performance of the different methods are explored using tight-binding time-dependent simulations in three trial structures, designated as ballistic, disordered, and resonant models. In the context of first-principles Hamiltonians the Driven Liouville-von Neumann approach is of special interest, because it does not require the calculation of Green's functions. Hence, the effects of replacing the reference density based on the Green's function by one obtained from an applied field are investigated, to gain a deeper understanding of the limitations and the range of applicability of the Driven Liouville-von Neumann equation.


Introduction
The interest in molecular conductance and electronic transport across nanostructures has inspired the development of a multiplicity of theoretical treatments to compute the current under an applied bias. The proposed schemes vary in complexity and computational cost, going from the original Landauer-Büttiker method 1,2 and different static models meant to describe steady state transport, [3][4][5][6][7] to dynamical methodologies that take into account the time evolution of the charge density. [6][7][8][9][10][11] The present article is concerned with the later class of approaches, and in particular with the subset of methods based on the so-called Driven Liouville-von Neumann (DLvN) equation. In this framework, the dynamics at the electrodes is modulated by an additional driving term that augments the standard Liouvillevon Neumann equation of motion. The role of these terms is to enforce part of the density matrix associated with the electrodes to remain close to a reference density matrix, thus introducing a charge imbalance between a "source" or left lead (L) and a "drain" or right lead (R). This driving term originally takes the form of a difference between the timedependent density matrixρ(t), and a reference density matrixρ 0 , whereĤ is the electron Hamiltonian, Γ is the driving rate parameter, and the matrixρ 0 can be defined as follows: This type of approach was first introduced by Horsfield and co-workers 12,13 as an intuitive way of including at the level of the density matrix the circulation of charge between the electrodes, and it was further enriched by the work of Nitzan 14 and Mazziotti. 15 These ideas were later taken on by Hod and co-workers [16][17][18][19] and reelaborated by deriving the driving terms from the theory of complex absorbing potentials. In doing so, they arrived at a modified form in which coherences were damped to zero, that improved the stability of the calculations and the steady-state convergence, and that was shown to satisfy Pauli's principle regardless of the initial conditions. 16 Franco and collaborators presented a formal derivation of this formula from the theory of non-equilibrium Green's functions, demonstrating that it can accurately capture time-dependent transport phenomena. 20 More recently, Zelovich et al. proposed a strategy to replace the single rate parameter Γ by diagonal matrices containing the broadening factors corresponding to the lead states. 21 These factors can be computed from the self-energies of the electrodes, thus providing a parameter-free version of the Driven Liouville-von Neumann approach. 21 Because of its conceptual simplicity and good compromise between computational cost and physical accuracy, the DLvN method has attracted significant attention, and several further refinements of its implementation and analysis of its theoretical foundation were made. Among these, of particular relevance is the adaptation to a first-principles real-time TDDFT framework carried out by Morzan and co-authors, 22 where an observed imbalance between injection and absorption of charge during the dynamics prompted a reformulation of the driving term. In an orthonormal basis, this equation of motion assumes the following structure:ρ where the subscripts indicate the corresponding blocks of the time-dependent density matrix, andρ 0 denotes a reference, time-independent density matrix. A major change with respect to previous expressions is that here all off-diagonal contributions are damped to their reference values (obtained e.g. from the polarized density matrix). This proved to be important for charge conservation and an appropriate balance between incoming and outgoing currents in the relatively small models tractable in TDDFT simulations. 22 Despite intuitively justifying this modification both from a numerical standpoint and interpreting it as a change in the boundary conditions, no formal theoretical derivation was provided at the time.
In parallel with these developments, a different driving term for the Liouville-von Neumann equation was put forward by Todorov and co-workers. 10 This new embedding method called multiple-probe or "hairy probes" (HP), formally derived a different structure of the driving term via the application of Green's functions and the Lippmann-Schwinger equation. 10,23 In spite of its good results and wide range of applicability, 24-26 its relationship with the reference density driving method remained undetermined.
The objective of the present study is to bridge the gap between the two methodologies: the heuristic formulation of the Driven Liouville-von Neumann equation, and the hairyprobes scheme. For that purpose, we have been able to rewrite the driving terms of the HP theory into two terms, one of which resembles the damping term of the reference density approach. In doing so we have found that the latest addition of Morzan and co-workers is the most compatible with HP, thus finally providing a formal framework for their correction of the non-diagonal blocks of the driving term. Through the application of these schemes to a series of model systems, we better characterize these methodologies for different situations and shed light on the reliability of the equation of motion including the driving term with the reference density matrixρ 0 .
In the next two sections we introduce the theoretical and methodological framework for this work: we first show how the HP theory can be rewritten to arrive to the Driven Liouville-von Neumann equation plus an additional term (section 2), and then provide a brief description of the models employed in the simulations (section 3). In section 4 we examine the impact that the progressive simplifications of the HP equation have on the accuracy of the physical description of these systems. After that we examine the effect that some of the relevant parameters of the models have on the baseline performance of these different approximations. In particular, we will focus on the Γ parameter (section 5) and the shape of the field used to generate the reference matrix in the DLvN scheme (section 6). Finally, in section 7, we discuss the differences between the two forms of the DLvN equation reported in the literature: the one emerging from the truncation of the HP method (eq. 4), and the one proposed by Hod and co-workers (eq. 3).

Recasting the multiple-probe equations of motions
We consider in this work a system formed by a central molecule or device, identified hereafter with the symbol C, coupled to a right and a left electrode (or lead), denoted as regions R and L respectively. In the multiple-probes framework this system is embedded in an implicitly represented environment via the coupling of each atom of the left and right leads to an external probe, which in isolation has a retarded and an advanced surface Green's function g ± j (E) and a surface local density of states d j (E) = −π −1 g + j (E). 10 These probes have fixed electrochemical potentials µ L and µ R , depending on whether they are connected to the left or right electrode, with Fermi-Dirac distributions f L (E) and f R (E). For this model, and assuming that the Hamiltonian of the system is time-independent, the HP theory provides the following equation of motion for the electronic density (see eq. 38 in reference 10 ): whereρ is the time-dependent density matrix,Ĥ is the electron Hamiltonian, subscript S refers to the full system (L, C and R sections), and theΣ andĜ matrices are the system's selfenergy and Green's function. We also assume that (i) all probes interact with the electrodes through the same coupling term γ, and (ii) the wideband limit applies to the probes so that their density of states becomes a constant, d j (E) = d. In these conditions, these matrices adopt the following form: In the above equations, Γ = 2πγ 2 d, and f L (E) and f R (E) correspond to the Fermi-Dirac distributions for the left and right probes. The matricesP L andP R are the corresponding projector operators, withP S being the projection over the whole explicit system (i.e. the identity).
We will now take equation 5 as the starting point for our derivation. Our goal in this section is to express it in terms of a difference between weighted density matrices-one corresponding to the current state of the system and the other to a reference system-that assumes the form of the driving term in the DLvN equations.
For this, we focus on the second (Σ +ρ terms on the right hand side of equation 5, and examine separately each of the sub-blocks corresponding to the orbitals of the electrodes (L, R) and of the device (C). The first thing to note is that, since allΣ matrices contain projections on the leads only, the second and third terms will vanish in the central block: In turn, the off-diagonal blocks of the second term can be expanded as: where for conciseness we have omitted the time-dependence ofρ. In matrix form, this can be written as: We now turn our attention to the third term on the right hand side of equation 5, considering separately each off-diagonal block. For the upper left one, we have: where we used the relation between the Green's function and the density of states matrix The integral in the last term of equation 12 defines a fictitious equilibrium density matrix with an electronic distribution f L corresponding to the left probes, that we will denoteρ 0L . After applying on the left and on the right the projection operator associated with the upper left block, we arrive at the final result, For the RR sub-matrix we get an analogous expression, but with the R projectors and the Fermi-Dirac distribution for the right probes, which defines a different reference density matrixρ 0R , To evaluate the off-diagonal terms involving the central region (blocks LC, CL, RC, and CR), we first observe that in all these cases part of the integrand vanishes when operating onΣ < due to the contiguous application of projectors belonging to different regions, Similarly, To recover the density of states from just one of the Green's matrices, we rewrite the latter in the following way: for which we have defined the matrixĜ S (E) = 1 2 (Ĝ − S (E) +Ĝ + S (E)). Inserting equations 17 and 18 in equations 15 and 16, we arrive at where we have introduced the matrixΩ L/R = f L/R (E)Ĝ S (E)dE. The terms involving the right lead can be treated on the same footing, to obtain an analogous result but incorporating the reference density andΩ corresponding to the right region, To complete the derivation, we consider the two remaining blocks that couple both leads together (LR, RL): The expression within the brackets can be rewritten introducing relations 17 and 18 to elicit the reference density and omega matrices, leading to the following result: which combines the left and right probes populated matrices. The RL block is treated analogously, to obtain: Thus, the third term in equation 5 can in matrix form be written as: where the first matrix on the right hand side will be referred to as the "region-weighted reference density matrix", or simply reference density matrix.
Finally, collecting all the pieces together, the HP master equation can be rewritten in the following way: We have thus arrived at a mathematical structure that is very similar to that of the (iii) same as (ii), but where the region-weighted reference density matrix is computed using a step-shaped field encoding the bias potential. In the latter, notice thatρ 0L andρ 0R are the same matrixρ 0 , corresponding to the polarized density in the presence of a step-field, and thus the equation of motion becomes identical to equation 4. These three schemes (i), (ii) and (iii), will be referred to as "full-HP", "partial-HP", and "step-potential" methods, or, for short, F-HP, P-HP, and ST-P, respectively. As mentioned above, the mathematical structure of the P-HP approach resembles closely the one of the DLvN equation 4, implemented in this study as the ST-P method. In particular, the driving term arising in the HP formula reinstates the damping of the coherences to the equilibrium values. In this sense, the ST-P scheme can also be thought of as an approximation to the P-HP method with an alternative reference density matrix. The motivation to explore it is to bypass the calculation of Green's functions, which becomes especially appealing in the context of ab-initio Hamiltonians.
For comparison purposes, we also consider in this work the alternative form of DLvN given by equation 3, which hereafter will be referred to as DLvN-z, whereas the version in equation 4 will be denoted as DLvN-e (since these two equations damp the lead-molecule coherences to zero and to their equilibrium values, respectively). It is important to note from this is that "ST-P" and "DLvN-e" refer to the exact same method: we will in general In all cases the systems are metallic, exhibiting a half-filled conduction band. These three systems are: • Ballistic: all atoms in the system (S) are equally spaced by 2.5Å, involving a hopping integral of -3.88 eV (which corresponds to the tight-binding reference parameters for gold 27 ). The on-site energy is the same for all atoms. There are 21 atoms in the central region C.
• Disordered: this system has the exact same geometry as the ballistic (same number of atoms all equally spaced), but with a different on-site energy for every atom belonging to the central region C. These different on-site energies were randomly generated using a homogeneous distribution with values between plus and minus the absolute value of the hopping term above.
• Resonant: all sites are equally spaced as in the previous cases, with the exception of two distances separating a group of 15 atoms in the center of the device, from two groups of 9 atoms on each end, connected to the electrodes. All of these belong to the C region, as can be seen in Figure 1. These two bonds were 0.5Å longer than the rest.
The on-site energies are the same for all atoms, but the hopping terms of the longer bonds are -1.88 eV.  The upper, middle and lower panels, correspond to the ballistic, the disordered, and the resonant models, respectively. Left and right panels show the currents for Γ equal to 0.1 eV and 0.6 eV. Figure 3 presents the steady state current-voltage curves obtained from the three methods deriving from HP. For the ballistic system all three approaches show only marginal differences between each other, which become negligible in the low bias region. The same is true for the disordered case. In both kinds of model system, the differences in the steady-state currents resulting from the three methods are below 5%, and typically much less than this at low biases. Somehow curiously, the scheme based on the step-potential reproduces more accurately the response of the full HP method, despite the difference in their reference densities. This suggests the existence of a compensation mechanism through which the stepgenerated reference density matrix balances the absence of theΩ contribution. A possible explanation can be found in its higher polarization, as discussed below.
The I − V curves in the resonant system display a more complex behavior. The P-HP and ST-P equations of motion still reproduce the currents of the full HP method, but with slightly larger deviations. Here, these appear to be comparable in the P-HP and in the ST-P method. The latter scheme produces in all three systems, for a given bias, higher currents than the first one. This trend was confirmed also in exploratory calculations for other resonant systems with variable separation between atoms: the step-potential yields systematically larger steady states currents than the P-HP method for a fixed potential difference. This might be related to the fact that in the case of a reference density generated from an abrupt step-field, the device is subject to a larger effective bias than in the case of a smoother reference density constructed from the Green's functions. In other words, the action of the field, through the direct modification of the on-site energies, affects more severely the density matrix of the electrodes when compared to the HP method, in which the electrochemical potential at the leads is controlled via the coupling with external probes according to γ. This is consistent with the eigenstates populations plotted in Figure 4. The full and the P-HP simulations produce a manifold of partially occupied states, both methods exhibiting essentially no differences from each other. When the reference density derives from the step-wise potential, instead, the eigenvalues of the density matrix are 0 or 1, with a

Effect of the Γ parameter
The question of the proper value of Γ is an interesting problem that has been investigated recently, 26 but has not been entirely settled. In the context of the DLvN approach, it has been argued in the literature that the magnitude of the driving rate should be somewhere between the energy spacing separating the lead levels, and the effective device-lead coupling. 20 If Γ is too high, the target density is enforced in the lead sections too tightly to allow for an appropriate response and relaxation at the interface. As a consequence, in this regime the current is a decreasing function of Γ (see e.g. refs. 20 , 16 , 22 ). Moreover, when Γ is too low, the relaxation process in the leads is not sufficiently fast to allow for the electronic density in the leads to reachρ 0 . The net effect in this case can be assimilated to a situation in which the leads are disconnected or weakly connected to the reservoirs and the system approaches a microcanonical regime.
On the other hand, in the context of hairy probes Γ depends on the system parameters d and γ, where d is the probe surface density of states and γ the matrix element coupling the sites at the electrodes with the external probes. While the value of d can be defined to fit the density of states of the device, there is no particular criterion to uniquely assign the value of γ. In principle, in the wideband limit it is expected to satisfy the second requirement listed above for Γ in DLvN: it must be larger than the mean energy level spacing in the leads, to ensure an effective broadening of the electronic states and hence a metallic behavior of the electrodes. In the present case this implies a lower bound for Γ of around 0.1 eV. Figure 5 illustrates the effect of Γ on the I − V plots in the case of the resonant system, also showing the derivatives of the current (conductance) in the insets. We chose this system because in its case the I − V curve displays a distinctive physical pattern, and the discrepancies between the three approaches are more significant than in any other model, both facts making improvements easier to identify. Indeed, these discrepancies tend to fade away as Γ decreases from 0.6 to 0.1 eV (within this range Γ remains larger than the mean energy level spacing in the leads, of 0.078 eV). Interestingly, the insets reveal that conductances derived from the ST-P methodology turn out to be more rugged than in the other methods. This can be attributed to the presence of the self-energy in the Green's functions, that, provided Γ is larger than the energy level spacing in the electrodes, screens the far ends of the system in the case of the F-HP and P-HP methods. In the absence of the self-energy, finite size effects manifest in the form of interference oscillations. As a matter of fact, for a small Γ the wiggles are visible in all three curves (see the inset on the right panel of Figure 5), because the screening effect dilutes as the Green's function carries a dependence on this parameter in the denominator. The insets also reveal that a large Γ value tends to damp the resonances in the ST-P case, which are otherwise intact for Γ = 0.1 eV.
The currents obtained for a bias of 1.5 V are depicted as a function of Γ in Figure 6, where it can be seen that the agreement between methods is progressive as the parameter gets smaller. Additionally, this Figure shows that all approaches are relatively robust against the variation of Γ, in particular the full HP scheme, for which the change in the steady state current is just 2% when this parameter is reduced by a factor of six.
In the F-HP and P-HP dynamics, the effect of Γ arises not only from the multiplication of the driving terms, but also because the reference matricesρ 0 andΩ are all functions of Γ through the self-energy and the Green's function (see equations 6 and 8). Therefore,  it becomes relevant to check how the structure of these matrices change when varying the value of Γ. Figures 7 and 8 show that the elements of these matrices exhibit a slight inverse dependence on Γ, which, as mentioned above, limits the spatial range of the Green's functions and reference density matrices in the electrode region, in a way that has no analogue in ST-P.
This spatial damping is visible on Figures 7 and 8. In particular, the magnitude of theΩ contribution is responsible for the difference between the F-HP and P-HP methods. All diagonal blocks of theΩ matrix are identically zero, and hence only the LC and LR blocks are depicted in Figure 7 (the remaining blocks, CL, RL, RC and CR are similar). The off-diagonal elements ofΩ are of the same order of magnitude as those of the reference matrixρ 0 . Yet, the impact of the former matrix on the dynamics is only secondary because it has zeros throughout its diagonal blocks, which is whereρ 0 gives the highest contribution. In spite of its seemingly small importance in numerical terms, though, the presence of theΩ matrix appears to balance the effect ofρ 0 in the current, furnishing the full HP method with a reduced sensitivity with respect to Γ. Figure 9 displays the difference between the reference density matrices in the P-HP and ST-P methods. Given that in the latter this matrix is independent of Γ, Figure 9 just reflects the effect of this parameter on the hairy probes reference density. Interestingly, since the decrease of Γ tends to relax the spatial damping of the P-HP reference density, which is in any case absent from the step-field generatedρ 0 , both matrices become with such a decrease more similar to each other. The agreement between these two approaches at low Γ is manifest in the behavior of the currents in Figures 5 and 6.
As mentioned in the introduction, Zelovich and co-workers proposed a method for the computation of a set of broadening factors that are applied to the lead states. In this way, the rate parameter Γ is replaced by diagonal matricesΓ L andΓ R , with dimensions given by the number of basis sets associated with the leads. The calculation of these matrices involves a self-consistent procedure to extract from the self-energy of the isolated lead, the broadening factors that afterwards must be assigned to the corresponding levels of the lead coupled to As a matter of fact, the authors reported that the adoption of the maximal broadening value calculated for such systems using their procedure as the single driving rate, produces current traces and steady-state occupations almost indistinguishable from those obtained using the parameter-free DLvN method.

Generation of the reference density
The encouraging results obtained with the ST-P implementation raise the question of whether it is possible to better reproduce the F-HP dynamics with a truncated equation of motion, through the optimization of the reference density. This pathway is considered in the present section, using the shape of the electric field as a tool to prepareρ 0 , although other strategies could be adopted as for example the state representation of Hod and co-workers. 16 We initially compare the results obtained with a step-like potential, with those corresponding to a ramp, a sigmoidal decay, and a double sigmoidal decay at the start and the end of the central region. These different patterns are illustrated in Figure 10. The analysis of the reference matrices generated in this way did not reveal any significant improvement in the description of the full hairy probes density, with the exception of a very slight increase in similarity for those matrix elements lying in the proximity of the central region. This improvement was most noticeable in the case of the sigmoid field (denoted as "SI-P"), as can be visualized in Figure 11.
As a matter of fact, the adoption of the sigmoidal field proved to be at least as good as the step potential, when not better, to reproduce the I − V plots. Figure 12 presents these curves for three resonant systems exhibiting different morphologies. Aside from the original structure displayed in Figure 1, two other models were examined in which the number of atoms in either the central or the lateral segments was extended. The value of Γ was fixed to 0.6 eV, for which the disagreement between methods was most noticeable. For the standard resonant system the performances obtained from the step or the sigmoidal potentials are comparable. The same is true for the alternative resonant model with a longer intermediate region, for which no significant differences are found between the results yielded by either method. However, in this case the description provided by these two approaches is manifestly worse. In this sense, it is noteworthy that the P-HP scheme is still able to capture the main features of the full hairy probes curve. On the other hand, in the third model where the device bears longer lateral segments, the reference calculated with the sigmoidal potential  to commute with the Green's functionsĜ ± S . In such a case, these contributions would vanish, suppressing ρ 0L and ρ 0R from the off-diagonal elements of the equation of motion, and thus giving rise to a driving term analogous to DLvN-z. Thus, the damping of the coherences either to zero or to equilibrium is related to this approximation in the Green's functions, which physically translates into an independency of the electron dynamics in the leads with respect to the device.
3. Finally, theΩ terms are dropped from the HP expression to arrive at DLvN-e. These terms arise from the Green's functions keeping track of the forward time-propagation of carriers coming from the probes. Thus, they introduce a propagation sense to the injected lead-device charge. Its suppression must imply a drop in the current, as is certainly observed in the simulations.
In ref. 22 it has been reasoned that while DLvN-z represents a device between electronic reservoirs at equilibrium with well-defined chemical potentials, DLvN-e resembles the state of a charged capacitor, where the target density represents the equilibrium state of the entire finite junction model and not just the leads. This seems to be consistent with the origins of each of these schemes, that have now come to light. Figure 13 depicts the current-voltage curves for the two DLvN approaches, together with the results of the HP method. At small couplings it can be seen that the performance of any of these methods is essentially indistinguishable from each other. As the Γ parameter is increased, discrepancies start to emerge.
DLvN-z has been shown to observe Pauli's principle regardless of the initial conditions. 16 It is interesting to note that, while the hairy probes scheme also fulfills Pauli's principle at long times (by construction), it does not necessarily obey it in the transient. This can be seen in Figure 14, which displays the occupations for F-HP, P-HP, DLvN-e, and DLvN-z in the resonant system. It must be recalled that in all these schemes, the dynamics is switched on smoothly, to avoid sudden jumps which could lead to numerical discontinuities. Specifically, the driving term and the Ω term are introduced in the first part of the simulation using a linear ramp. For F-HP, the positive and negative deviations from 1 and 0 respectivelywhich become smaller with a smoother ramp-disappear in the long term. When the Ω term is suppressed from the F-HP scheme, the dynamics and in particular the occupations are affected but the exclusion principle is still obeyed in the steady state. However, if the reference density matrix is replaced by the one obtained from an electric field to give the DLvN-e approach, the violation to the exclusion principle persists even in the steady state. This is not observed with the DLvN-z method, for which the exclusion principle is satisfied throughout the full dynamics, despite the adoption of the sameρ 0 as in DLvNe. Interestingly, this makes manifest that the observance of Pauli's principle is determined neither by Ω nor by the reference density, but by their combination. Finally, Figure 15 explores the effect of Γ and of the electrode size on the currents, using  of motion is that, at variance with the previous versions introduced in references 13 or 16 , the coherences are damped to the equilibrium density. It has been argued that in this approach the electrodes are not meant to represent infinite reservoirs with homogeneous and well defined chemical potentials. Instead, through the action of the driving operator, the leads are driven close to, rather than exactly to, the target density. 22  To neglect theΩ matrix in the driving term of the HP formula has minor effects on the dynamics and the steady state currents obtained for a variety of model systems. These effects are even less significant when the coupling between the probes and the leads is reduced by decreasing the Γ parameter. More specifically, our results show that the P-HP method can reproduce the behavior described by the full HP scheme for ballistic and disordered systems, and with some care it can also be tuned for more complicated resonant systems. In general, at least in the context of tight-binding Hamiltonians, it is possible to conclude that the DLvN method, incarnated here in the P-HP or ST-P equations of motion, converges to the hairy probes description in the limit of small couplings between the probes and the leads (providing Γ is still larger than the energy level spacing). The P-HP method can be thought of as a version of the DLvN approach in which the calculation of the reference density is based on Green's functions. The ST-P method, on the other hand, reproduces the formulation presented in reference 22 . It has proved to be a very good approximation to the full HP description, but the finite size effects which manifest in the absence of a self-energy imply a limitation in comparison with the P-HP approach that seems difficult to overcome.
The possibility to avoid the calculation of Green's functions acquires special interest for the applications in the context of TDDFT or other first-principles schemes. In that respect, we explored the substitution of the reference density computed from Green's functions in the P-HP method, by the one emerging from an electric field applied to the system. This strategy systematically produced higher currents than the P-HP method, presumably because the application of a field, in the present setting, amounts to an additive constant in the density matrix elements, whereas in the HP scheme the chemical potential is fixed at an external probe and not directly on the leads, whose polarization is mediated by a coupling parameter. The effective bias arising between the electrodes in operation conditions may thus not be as large as the one imposed by the field. Further sources of improvement for the region-weighted field-generated reference density method proved hard to find. Our results here seem to suggest that the step and sigmoideal shapes fields can better fit those parts of the reference matrices at the boundary with the device, providing the most accurate representations of the current in comparison with the other fields tested. Whereas the ballistic and disordered models could be described reasonably well with the ST-P and SI-P schemes, the resonant systems proved to be challenging for those methods using electric-field based reference densities. This limitation becomes particularly relevant for TDDFT applications, where the complexity of the chemical structures, far from the simple tight-binding models examined in the present work, may lead to stronger discrepancies between the results obtained using a reference density generated from Green's functions or from an electric field.
It would be interesting to establish how the different reference densities and the omission of theΩ matrix affect the dynamics in the case of first-principles simulations. That question will be the subject of future work.