Charge nonconservation of molecular devices in the presence of a nonlocal potential

TL;DR

This paper addresses charge nonconservation in molecular devices with nonlocal potentials, identifying its origin in approximations and pseudopotential definitions, and proposes a formulation to ensure charge conservation with practical application results.

Contribution

It introduces a nonlocal-potential formulation that guarantees charge conservation in molecular device simulations involving nonlocal potentials.

Findings

Charge nonconservation stems from approximation issues, not current definition.

The proposed formulation restores charge conservation in pseudopotential calculations.

Application to a double-barrier model demonstrates the effectiveness of the approach.

Abstract

In the presence of a nonlocal potential in molecular device systems, generally the charge conservation cannot be satisfied, and in literatures the modifications of the conventional definition of current were given to solve this problem. We demonstrate that, however, the nonconservation is not due to the invalidation of the conventional definition of current, but originates respectively from the improper approximations to electron-electron interactions and the inappropriate definition of current using pseudo wave functions in pseudopotential implementations. In this work, we propose a nonlocal-potential formulation of the interactions to fulfill the charge conservation and also give a discussion about the calculation of current when the pseudopotential is involved. As an example of application of our formulation, we further present the calculated results of a double-barrier model.

Equations

$${\mathbf{J}}_c=\frac{-i\hslash }{2m}\left[{\psi }^{\ast }\left(\mathbf{r}\right)\nabla \psi \left(\mathbf{r}\right)-\psi (\mathbf{r})\nabla {\psi }^{\ast }(\mathbf{r})\right].$$

$${\hat{H}}_0={\hat{T}}_s+{\hat{U}}_{ex},$$

$${\hat{T}}_s=\frac{-{\hslash }^2}{2m}\int d\mathbf{r}{\hat{\Psi }}^{†}\left(\mathbf{r},t\right){\nabla }^2\hat{\Psi }\left(\mathbf{r},t\right),$$

$${\hat{U}}_{ex}=\int d\mathbf{r}{\hat{\Psi }}^{†}\left(\mathbf{r},t\right)v\left(\mathbf{r},t\right)\hat{\Psi }\left(\mathbf{r},t\right).$$

$${\hat{\mathbf{J}}}_c(\mathbf{r},t)=\frac{-i\hslash }{2m}\left[{\hat{\Psi }}^{†}\left(\mathbf{r},t\right)\nabla \hat{\Psi }\left(\mathbf{r},t\right)-\left[\nabla {\hat{\Psi }}^{†}\left(\mathbf{r},t\right)\right]\hat{\Psi }\left(\mathbf{r},t\right)\right],$$

$$\frac{\partial \hat{n}\left(\mathbf{r},t\right)}{\partial t}=$$

$$=$$

$$-\left[{\nabla }^2{\hat{\Psi }}^{†}\left(\mathbf{r},t\right)\right]\hat{\Psi }\left(\mathbf{r},t\right)]$$

$$=$$

$${\hat{H}}_w={\hat{T}}_s+{\hat{U}}_{ex}+\hat{W},$$

$$\hat{W}=\int d{\mathbf{r}}^{\prime \prime }d{\mathbf{r}}^{\prime }w\left(\left|{\mathbf{r}}^{\prime }-{\mathbf{r}}^{\prime \prime }\right|\right){\hat{\Psi }}^{†}\left({\mathbf{r}}^{\prime }\right){\hat{\Psi }}^{†}\left({\mathbf{r}}^{\prime \prime }\right)\hat{\Psi }\left({\mathbf{r}}^{\prime \prime }\right)\hat{\Psi }\left({\mathbf{r}}^{\prime }\right).$$

$$\left[\hat{n}\left(\mathbf{r},t\right),\hat{W}\right]$$

$$=[{\hat{\Psi }}^{†}\left(\mathbf{r},t\right)\hat{\Psi }\left(\mathbf{r},t\right),\int d{\mathbf{r}}^{\prime \prime }d{\mathbf{r}}^{\prime }w\left(\left|{\mathbf{r}}^{\prime }-{\mathbf{r}}^{\prime \prime }\right|\right){\hat{\Psi }}^{†}\left({\mathbf{r}}^{\prime },t\right)$$

$$\times {\hat{\Psi }}^{†}\left({\mathbf{r}}^{\prime \prime },t\right)\hat{\Psi }\left({\mathbf{r}}^{\prime \prime },t\right)\hat{\Psi }\left({\mathbf{r}}^{\prime },t\right)]$$

$$=\int d{\mathbf{r}}^{\prime \prime }d{\mathbf{r}}^{\prime }w\left(\left|{\mathbf{r}}^{\prime }-{\mathbf{r}}^{\prime \prime }\right|\right)$$

$$\times [{\hat{\Psi }}^{†}\left(\mathbf{r},t\right)\delta \left(\mathbf{r}-{\mathbf{r}}^{\prime }\right){\hat{\Psi }}^{†}\left({\mathbf{r}}^{\prime \prime },t\right)\hat{\Psi }\left({\mathbf{r}}^{\prime \prime },t\right)\hat{\Psi }\left({\mathbf{r}}^{\prime },t\right)$$

$$-{\hat{\Psi }}^{†}\left(\mathbf{r},t\right){\hat{\Psi }}^{†}\left({\mathbf{r}}^{\prime },t\right)\delta \left(\mathbf{r}-{\mathbf{r}}^{\prime \prime }\right)\hat{\Psi }\left({\mathbf{r}}^{\prime \prime },t\right)\hat{\Psi }\left({\mathbf{r}}^{\prime },t\right)$$

$$+{\hat{\Psi }}^{†}\left({\mathbf{r}}^{\prime },t\right){\hat{\Psi }}^{†}\left({\mathbf{r}}^{\prime \prime },t\right)\delta \left(\mathbf{r}-{\mathbf{r}}^{\prime \prime }\right)\hat{\Psi }\left({\mathbf{r}}^{\prime },t\right)\hat{\Psi }\left(\mathbf{r},t\right)$$

$$-{\hat{\Psi }}^{†}\left({\mathbf{r}}^{\prime },t\right){\hat{\Psi }}^{†}\left({\mathbf{r}}^{\prime \prime },t\right)\hat{\Psi }\left({\mathbf{r}}^{\prime \prime },t\right)\delta \left(\mathbf{r}-{\mathbf{r}}^{\prime }\right)\hat{\Psi }\left(\mathbf{r},t\right)]$$

$$=0,$$

$${\hat{H}}_w={\hat{H}}_0+{\hat{V}}_{xc},$$

$${\hat{V}}_{xc}=\int d\mathbf{r}d{\mathbf{r}}^{\prime }{V}_{xc}({\mathbf{r}}^{\prime },\mathbf{r}){\hat{\Psi }}^{†}({\mathbf{r}}^{\prime },t)\hat{\Psi }(\mathbf{r},t).$$

$$\left[\hat{n}(\mathbf{r},t),{\hat{V}}_{xc}\right]$$

$$=\left[{\hat{\Psi }}^{†}\left(\mathbf{r}\right)\hat{\Psi }\left(\mathbf{r}\right),\int d{\mathbf{r}}^{\prime \prime }d{\mathbf{r}}^{\prime }{V}_{xc}({\mathbf{r}}^{\prime \prime },{\mathbf{r}}^{\prime }){\hat{\Psi }}^{†}({\mathbf{r}}^{\prime \prime },t)\hat{\Psi }({\mathbf{r}}^{\prime },t)\right]$$

$$=\int d{\mathbf{r}}^{\prime \prime }d{\mathbf{r}}^{\prime }{\hat{\Psi }}^{†}\left(\mathbf{r},t\right)V_{xc}({\mathbf{r}}^{\prime \prime },{\mathbf{r}}^{\prime })\delta \left(\mathbf{r}-{\mathbf{r}}^{\prime \prime }\right)\hat{\Psi }({\mathbf{r}}^{\prime },t)$$

$$-\int d{\mathbf{r}}^{\prime \prime }d{\mathbf{r}}^{\prime }{V}_{xc}({\mathbf{r}}^{\prime \prime },\mathbf{r}){\hat{\Psi }}^{†}({\mathbf{r}}^{\prime \prime },t)\delta \left(\mathbf{r}-{\mathbf{r}}^{\prime }\right)\hat{\Psi }\left(\mathbf{r},t\right)$$

$$=\int d{\mathbf{r}}^{\prime }{\hat{\Psi }}^{†}\left(\mathbf{r},t\right)V_{xc}(\mathbf{r},{\mathbf{r}}^{\prime })\hat{\Psi }({\mathbf{r}}^{\prime },t)$$

$$-\int d{\mathbf{r}}^{\prime \prime }{\hat{\Psi }}^{†}({\mathbf{r}}^{\prime \prime },t)V_{xc}({\mathbf{r}}^{\prime \prime },\mathbf{r})\hat{\Psi }\left(\mathbf{r},t\right)$$

$$\ne 0,$$

$${\hat{n}}_{nl}(\mathbf{r},t)=\frac{1}{i\hslash }$$

$$-{\hat{\Psi }}^{†}({\mathbf{r}}^{\prime },t)V_{xc}({\mathbf{r}}^{\prime },\mathbf{r})\hat{\Psi }\left(\mathbf{r},t\right)],$$

$$\frac{\partial \hat{n}\left(\mathbf{r},t\right)}{\partial t}+\nabla \cdot {\hat{\mathbf{J}}}_c\left(\mathbf{r},t\right)={\hat{n}}_{nl}(\mathbf{r},t).$$

$$\hat{W}$$

$$\approx \int d\mathbf{r}{\Psi }^{†}(\mathbf{r})\left[\int d{\mathbf{r}}^{\prime }\frac{⟨{\Psi }^{†}({\mathbf{r}}^{\prime })\Psi ({\mathbf{r}}^{\prime })⟩}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}\right]\Psi (\mathbf{r})$$

$$-\int d\mathbf{r}d{\mathbf{r}}^{\prime }{\Psi }^{†}(\mathbf{r})\frac{⟨{\Psi }^{†}({\mathbf{r}}^{\prime })\Psi (\mathbf{r})⟩}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}\Psi ({\mathbf{r}}^{\prime })$$

$$=\int d\mathbf{r}{\Psi }^{†}(\mathbf{r})U_H(\mathbf{r})\Psi (\mathbf{r})$$

$$+\int d\mathbf{r}d{\mathbf{r}}^{\prime }{\Psi }^{†}(\mathbf{r})U_{xc}({\mathbf{r}}^{\prime },\mathbf{r})\Psi ({\mathbf{r}}^{\prime }),$$

$$U_H(\mathbf{r})=\int d{\mathbf{r}}^{\prime }\frac{⟨{\Psi }^{†}({\mathbf{r}}^{\prime })\Psi ({\mathbf{r}}^{\prime })⟩}{|{\mathbf{r}-\mathbf{r}}^{\prime }|},$$

$$U_{xc}({\mathbf{r},\mathbf{r}}^{\prime })=-\frac{⟨{\Psi }^{†}({\mathbf{r}}^{\prime })\Psi (\mathbf{r})⟩}{|{\mathbf{r}-\mathbf{r}}^{\prime }|}.$$

$$\hat{F}(\mathbf{r})=\frac{1}{i\hslash }$$

$$-U_{xc}({\mathbf{r}}^{\prime },\mathbf{r}){\Psi }^{†}({\mathbf{r}}^{\prime })\Psi (\mathbf{r})],$$

$$\frac{\partial \hat{n}\left(\mathbf{r},t\right)}{\partial t}+\nabla \cdot {\hat{\mathbf{J}}}_c\left(\mathbf{r},t\right)=\hat{F}(\mathbf{r}),$$

Full text

Current density; Nonlocal potential; Nonconservation; Interaction

I Introduction

Over the past few decades, investigations on transport properties of mesoscopic systems and nanostructures have been extensively reported both on experimental advances1 ; 2 ; 3 ; 4 ; 5 ; 6 ; 7 ; 8 ; 9 ; 10 ; 11 ; 12 ; 13 and theoretical explorations.14 ; 15 ; 16 ; 17 ; 18 ; 19 ; 20 ; 21 ; 22 ; 23 ; 24 ; 25 ; 26 ; 27 It is widely acknowledged that these functional devices can be constituted by ultrasmall conjugated molecules, single-layer or multi-layer nanotubes, bulk organic molecules, etc., and plenty of interesting phenomena such as molecular field effects,1 Coulomb blockade,2 negative differential resistance3 and conductance switching effects4 have been revealed, which exhibit fundamental significance and potential microelectronic applications. In most of the works, considerable research efforts are focused on current-voltage ($I$-$V$) characteristics as the $I$-$V$ profiles provide opportunities for a deeper understanding of, e.g., the basic mechanism and structure properties, as well as promising guidance for future molecular nanoelectronics designs and manipulations.

On the theoretical side, calculations for the $I$-$V$ characteristics of molecular device systems are mostly performed by employing the self-consistent field (SCF) theory or nonequilibrium Green’s functions combined with density functional theory (NEGF-DFT),18 ; 19 ; 20 and the widely used DFT calculations at present can be vested in the SCF method. In comparison with conventional SCF,28 ; 29 ; 30 in addition to self-consistent Hartree potential DFT introduces the exchange-correlation potential, which is nonlocal if one wants to go beyond the local density approximation.31 It has been shown that if the Hamiltonian $H$ includes a general nonlocal potential $V(\mathbf{r},\mathbf{r}^{\prime})$, an extra term naturally appears in the continuity equation and the charge conservation will not be fulfilled.24 ; 25 ; 32 ; 33 Thus the calculations may give very incorrect results, even nonphysical results.

To resolve the problem, Li et al.32 proposed a scheme to modify the conventional definition of current density to include the additional current induced by the nonlocal potential $V(\mathbf{r},\mathbf{r}^{\prime})$, and therefore yield the charge conservation in a computationally efficacious way. However, either local or nonlocal exchange-correlation potential stems from the approximation to electron-electron interactions, and according to the conventional definition of current density we will demonstrate fundamentally that with the Hamiltonian including the exact electron-electron interactions the extra term does not appear and the conservation can be precisely satisfied. Hence the problem of charge nonconservation coming from the nonlocal exchange-correlation potential should not be settled by redefining the current density. Instead, it has to be resolved by finding a reasonable nonlocal-potential approximation to the electron-electron interactions to eliminate the extra term. On the other hand, norm-conserving pseudopotentials 34 are generally utilized to reduce the size of plane-wave basis sets in first-principles calculations, which is another origin of the nonlocal potential.35 Nevertheless, the pseudopotential implementations give the pseudo wave functions, while in the continuity equation $\partial_{t}\rho+\nabla\cdot\mathbf{J}=0$ the charge density $\rho$ and the current density $\mathbf{J}$ should be calculated by using the true wave functions instead of the pseudo ones. Therefore, when the pseudopotential method is employed, the extra term appearing in the continuity equation is due to the improper utilization of the pseudo wave functions for calculating the current, and the problem of charge nonconservation coming from the nonlocal pseudopotential should not be resolved by redefining the current density either. From a fundamental point of view, the continuity equation is a criterion regardless of any approximations brought in as long as the particles of the system are conserved, which can be easily proved with the original Hamiltonian. It is the purpose of this work to investigate the above problems.

The paper is organized as follows. In Sec. II, the origins of the nonlocal potential and the consequent issues of charge nonconservation are discussed. We lay special emphasis on the nonlocal exchange-correlation potential from the starting point of second quantization, and subsequently demonstrate that in DFT calculations the nonconservation is caused by the inappropriate definition of current using pseudo wave functions rather than the introduction of the pseudopotentials. As an example, the currents of a double-barrier model are numerically calculated in Sec. III to confirm our theoretical formulation. Section IV gives the conclusions and discussions.

II Theoretical Formulation

In the context of first-principles calculations, if one uses true wave functions $\psi\left(\mathbf{r}\right)$ of the system from the very beginning throughout the processes, the continuity equation $\partial_{t}\rho_{c}+\nabla\cdot\mathbf{J}_{c}=0$ can be easily realized according to the Schördinger equation, where $\rho_{c}=\left|\psi\left(\mathbf{r}\right)\right|^{2}$ is the conventional electron density and $\mathbf{J}_{c}$ is the conventional current density in the absence of magnetic field with the definition as

[TABLE]

One of the origins of the charge nonconservation comes from the approximation to electron-electron interactions, i.e., some improper exchange-correlation potentials are introduced. We will show that according to the conventional definition of current density, the conservation is still satisfied in the presence of the interactions, but generally can be violated by introducing the nonlocal exchange-correlation potentials. To see this, the simplest case of Hamiltonian of a finite many-electron system is considered, in which the electron-electron interactions are not taken into account firstly so that the second quantized nonrelativistic Hamiltonian (a quantity with a caret symbol denotes an operator) is

[TABLE]

where in terms of the field operators $\hat{\Psi}^{{\dagger}}\left(\mathbf{r},t\right)$ and $\hat{\Psi}\left(\mathbf{r},t\right)$, the single-particle kinetic energy operator and the external potential operator can be written respectively as36 ; 37

[TABLE]

and

[TABLE]

The density operator can be defined as $\hat{n}(\mathbf{r},t)=$ $\hat{\Psi}^{{\dagger}}\left(\mathbf{r},t\right)\hat{\Psi}\left(\mathbf{r},t\right)$ and the current density operator as the conventional form is

[TABLE]

by means of the Heisenberg’s equation and the anticommutation relation that we obtain

[TABLE]

where the charge conservation is accomplished, as expected. The corresponding Hamiltonian including the interactions can be written as

[TABLE]

where

[TABLE]

Since the foregoing equations present the continuity with respect to $\hat{T}_{s}$ and $\hat{U}_{ex}$, only interaction operator is taken into account hereafter, i.e.,

[TABLE]

indicating that the density operator $\hat{n}$ is commutable with the interaction Hamiltonian $\hat{H}_{w}$, and hence the continuity equation (6) is still satisfied when the interaction is involved. In fact, the conventional definition of current and particle densities are widely used in time-dependent current-density functional theory (TDCDFT), and the continuity equation is frequently utilized as a constraint condition between the two densities. 45 ; 45-1 ; 45-2 ; 45-3

Next, we replace the electron-electron interaction $\hat{W}$ by a local potential and a nonlocal exchange-correlation potential. Since the local potential can be absorbed by the external-field potential $\hat{U}_{ex}$, and the interaction $\hat{W}$ can be replaced only by nonlocal exchange-correlation potential $\hat{V}_{xc}$, thus the Hamiltonian $\hat{H}_{w}$ is approximated (see Appendix A) as

[TABLE]

where

[TABLE]

After some calculations we reach

[TABLE]

where the continuity condition in Eq. (6) is no longer satisfied, and it contains an extra term

[TABLE]

i.e.,

[TABLE]

Thus far an inference can be drawn that the charge nonconservation is not due to the invalidation of the conventional definition of current density, but originates from the improper approximations to electron-electron interactions, and a reasonable nonlocal potential arising from the interactions should make the extra term be zero. This is the central conclusion of this work, and the task of resolving the problem is to find such a “no-current”nonlocal potential. Now, we take a mean-field approximation to the electron-electron interactions and rewrite the potential energy as

[TABLE]

where $\langle...\rangle$ represents the ensemble statistical average,

[TABLE]

and

[TABLE]

Therefore we have obtained our nonlocal exchange-correlation potential. Introducing an auxiliary variable $\hat{F}$ as

[TABLE]

and similarly we have

[TABLE]

where an extra term is still contained in the meaning of operator. However, after taking the ensemble statistical average 38 of $\hat{F}(\mathbf{r}),$ we can obtain

[TABLE]

i.e.,

[TABLE]

which indicates that in the meaning of statistical average the conservation can be again fulfilled. This result is quite satisfactory, because an observable quantity is a statistical-average one. One may find that the nonlocal exchange-correlation potential in Eq. (17) is similar to the Fock term of the Hartree-Fock approximation. In the case of zero temperature, this proposed nonlocal exchange-correlation potential is just the Fock term, and in other cases they are different. Normal Hartree-Fock approximation is only applicable to the ground state, while Eq. (17) can be used in nonequilibrium states.

Another origin of the nonconservation is the inappropriate definition of current density in pseudopotential implementations. To conserve the current, one should calculate the current density $\mathbf{J}$ by using the true wave function $\psi\left(\mathbf{r}\right)$ rather than the generally obtained pseudo wave function $\phi\left(\mathbf{r}\right)$ (see Appendix B). Note that in practical calculations, the strict pseudopotential would be somewhat difficult to be used due to the explicit energy dependence, and approximate model pseudopotentials are often introduced instead. As shown in Eq. (40), the nonlocal part of the pseudopotential comes from the core wave functions, which is confined in a small core region surrounding the atom, and the true and pseudo wave functions are identical outside this small region.34 ; 39 Thus the additional current induced by the nonlocal pseudopotential is localized in the core of the atom, and its contribution to the transport current can be neglected.

III Numerical Implementation

To illustrate the theoretical formulation proposed above, we consider a quasi-one-dimensional double-barrier model confined transversely to simulate the device system, as shown in Fig. 1. The nonlocal potential is placed only in the region between $-a$ and $a$, and the barriers can be regarded as part of two ideal leads without nonlocal potentials.

We imagine a central region C (the dashed-line box) that encloses both the core device and the barriers, and the total Hamiltonian for this model can be written as

[TABLE]

where $\hat{H}_{0}=\frac{-\hbar^{2}}{2m}\nabla^{2}+U_{0}$ with $U_{0}$ the local potential within the barriers, $\hat{V}_{H}$ is the Hartree potential and $\hat{V}_{xc}$ is the exchange-correlation potential: $\langle\mathbf{r|}\hat{V}_{xc}|\psi\rangle=\int V_{xc}\left(\mathbf{r,r}^{\prime}\right)\psi(\mathbf{r}^{\prime})d\mathbf{r}^{\prime}$. We simplify the potential in Eq. (16) and Eq. (17) to the form of one dimension as

[TABLE]

and

[TABLE]

where $\psi_{k}\left(x\right)$ is the longitudinal wave function, and

[TABLE]

which is averaged over the transverse wave functions.40 A single transverse wave function $\varphi\left(\rho\right)=1/\sqrt{2\pi l^{2}}\exp\left(-\rho^{2}/4l^{2}\right)$ is chosen universally, and a transverse radius $2l$ characterizes the size of the confinement. To proceed, we first solve the energy eigenequation of the double-barrier model with only local potential $U_{0}$ and numerically calculate the wave functions $\psi\left(x\right)$, thus the corresponding Hartree potential $V_{H}$ and exchange-correlation potential $V_{xc}$ of the central region can be constructed from the wave functions* *by using Eq. (23) and Eq. (24). Then we return to the calculation of the wave functions $\psi\left(x\right)$, and self-consistent calculations are performed using the iterative procedure in the numerical implementation. In general, once all the wave functions with eigenenergies below the Fermi energy $E_{F}$ are known, we can obtain the current density at all points. When the system is under applied bias, the Fermi levels in the left and right leads are taken respectively as $E_{F}+eV$ and $E_{F},$ where $E_{F}$ is the Fermi level in the equilibrium state and $V$ is the bias, i.e., the sums in Eq. (23) and Eq. (24) are up to $E_{F}+eV$ and $E_{F}$ respectively for left-incident and right-incident modes of the functions. Energy unit $\hbar^{2}/2ma^{2}$ is used throughout the calculations as a reduced coefficient, where $a$ is the width of the barrier, $m$ is the mass of the electron. The Fermi energy here is set as $E_{F}=1.0$, and the magnitude of the barriers is fixed to be $U_{0}=6.0$. In addition, all the calculations are preformed in low-temperature limit, i.e., the temperature $T=0$.

We first study the influences of a model nonlocal potential and our proposed potential in Eq. (17) on the particle currents, respectively. The nonlocal potential is nonzero when $x$ is from $-1$ to $1$, and the barriers are located in $\left[-2,-1\right]$ and $\left[1,2\right]$. The model nonlocal potential here is chosen as $V\left(x,x^{\prime}\right)=\lambda\exp\left[\eta\left(x-x^{\prime}\right)^{2}\right]$, where $\lambda$ and $\eta$ are two independent coefficients. In Fig. 2, we present the calculated currents of left-incoming electrons with a fixed energy $E=0.8$ and below the Fermi energy $E_{F}$, respectively. Here, both the conventional particle current and additional current are given, and the additional currents coming from nonlocal potential $V(x,x^{\prime})$ are calculated by

[TABLE]

From Fig. 2(a) to Fig. 2(d), one can see that for single-particle with the fixed energy the left and right lead currents hold to be constant, while in the central region the currents turn out to be varying with $x$, and the additional currents from the nonlocal potential are nonzero for both the model potential and our proposed potential. It means that in the case of single-particle, the lead currents and the central current calculated with the proposed nonlocal potential behave similarly to that obtained by using the model nonlocal potential, which also violates the conservation. However, Figures 2(e) and 2(f) show that when currents of all incoming electrons below the Fermi energy $E_{F}$ are considered, in the case of the model nonlocal potential the conventional current still violates the continuity condition due to the nonzero additional current, while the conventional currents calculated with the proposed nonlocal potential are seen to be conservative in the whole simulation region and the corresponding additional current vanishes. Moreover, particle currents in all the three regions serve as a minor correction to the ones with Hartree potential only (see Fig. 3).

Figure 4 plots the dependence of the electric current on the selection of each region. The width of the barriers are fixed as the above while the whole central region is changeable, and the nonlocal potential is absent therein. As we shall see, when the core region is chosen as small as $x=[-0.2,0.2]$, it exhibits a sublinear relation between the current $I$ and the voltage $V$, which is similar to that of a general single-barrier model, and the scattering effect in the central box vanishes. With the increases of the width as $x=\left[-0.5,0.5\right]$ and $x=\left[-1.0,1.0\right]$, the curves shift downwards gradually and finally become nonlinear ones, indicating that careful considerations on relevant area selection must be taken into actual systems to ensure the computational accuracy.

In Fig. 5, we present the $I$-$V$ characteristics of different sites along the $x$ axis when the above model nonlocal potential $V\left(x,x^{\prime}\right)$ is included. Currents of $x=0$, $x=0.5$ and $x=2.0$ differ mutually, as expected, which implies that currents in the central region and the leads calculated with the model nonlocal potential cannot meet the charge conservation. We further explore the contrastive characteristics with our approximation. From Fig. 6, one can see that when only Hartree potential $V_{H}$ is taken into account, currents of the three sites are equal to each other with the increase of the voltage, while the whole of the curve bears an upward shift when the proposed exchange-correlation potential $V_{xc}$ is also considered, and the currents obtained can be still conservative.

IV Conclusions and Discussions

In summary, it is of great importance to give correct dynamic charge and potential distributions of the transport systems, which is the key point for the valid currents and prospective applications of molecular devices. Once we are able to calculate the current strictly according to the original Hamiltonian including electron-electron interactions, the problems of charge nonconservation induced by nonlocal potential would not exist, and the results are undoubtedly reasonable. However, it is impossible to obtain rigorous solutions under such intricate interactions. We demonstrate these issues and attest that the nonconservation stems respectively from the improper approximations to electron-electron interactions and the inappropriate definition of current using pseudo wave functions in pseudopotential implementations, and propose a nonlocal-potential formulation of the interactions to fulfill the conservation in the meaning of statistical average, as well as give a verification about the calculation of current in the presence of pseudopotentials. With this method, we have also studied a double-barrier model to simulate the molecular device system, which further confirms our formulation.

ACKNOWLEDGEMENTS

This work is supported by National Natural Science Foundation of China under Grant No. 11675051.

Appendix A

Here, we will attempt to give the second quantized form of the nonlocal potential. According to common methods, the wave function $\psi$ can be expanded in terms of a complete basis set $\varphi_{m}$ as

[TABLE]

where $c_{m}$ is a coefficient independent of $\mathbf{r}$. Substituting the above wave function into the following Schrödinger equation, which includes the nonlocal potential $V_{xc}(\mathbf{r},\mathbf{r}^{\prime})$,

[TABLE]

we have

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

is the matrix elements of $V_{xc}(\mathbf{r},\mathbf{r}^{\prime})$. Similarly, the field operator is

[TABLE]

with

[TABLE]

where $\hat{a}_{n}^{{\dagger}}$ $\left(\hat{a}_{n}\right)$ is the creation (annihilation) operator. Therefore, we obtain the second quantized form of the nonlocal potential as

[TABLE]

Appendix B

We first discuss that in introducing the nonlocal pseudopotential, the current calculated with pseudo wave functions does not satisfy the conservation. Pseudopotentials were originally introduced to simplify electronic structure calculations by adding some core functions to the true wave function $\psi\left(\mathbf{r}\right)$ to obtain a smooth pseudo wave function41 ; 42 ; 43

[TABLE]

where $\phi\left(\mathbf{r}\right)$ is the pseudo wave function and $\psi_{n}\left(\mathbf{r}\right)$ is the core function, and the general form of the pseudopotential (only its nonlocal part) is

[TABLE]

from which we obtain

[TABLE]

where $G=\left(E-H_{0}\right)^{-1}$ is the Green’s function, and $H_{0}$ is the original Hamiltonian. The true wave function can then be denoted by

[TABLE]

Therefore, the corresponding Schrödinger equation for the pseudo wave function is

[TABLE]

where

[TABLE]

It is easy to verify that the conventionally defined current density $\mathbf{J}=\frac{-i\hbar}{2m}[\phi^{\ast}\left(\mathbf{r}\right)\nabla\phi\left(\mathbf{r}\right)-\phi\left(\mathbf{r}\right)\nabla\phi^{\ast}\left(\mathbf{r}\right)]$ along with the obtained pseudo electron density $\rho_{ps}=\left|\phi\left(\mathbf{r}\right)\right|^{2}$ do not meet the charge conservation. Nevertheless, this nonconservation is not caused by introducing the pseudopotential, but is that the above definition of current density cannot be used with the pseudo wave function $\phi\left(\mathbf{r}\right)$. The correct current density must be defined by using the true wave function $\psi\left(\mathbf{r}\right)$. According to the definition $\mathbf{J}_{r}=$ $\frac{-i\hbar}{2m}\left[\psi^{\ast}\left(\mathbf{r}\right)\nabla\psi\left(\mathbf{r}\right)-\psi(\mathbf{r})\nabla\psi^{\ast}(\mathbf{r})\right]$, we get

[TABLE]

where

[TABLE]

and

[TABLE]

Here, $\Pi(\mathbf{r)=}\int d\mathbf{r}^{\prime}d\mathbf{r}^{\prime\prime}G(\mathbf{r,r}^{\prime})V_{nl}(\mathbf{r}^{\prime},\mathbf{r}^{\prime\prime})\phi\left(\mathbf{r}^{\prime\prime}\right)$. On the basis of Eq. (39), in steady states we can easily give $\nabla\cdot\mathbf{J}_{ps}=-\nabla\cdot\mathbf{J}_{nl}$, which verifies the conservation with the definition of current density in Eq. (41).