Nonequilibrium Green’s-Function Approach to the Suppression of Rectification at Metal–Mott-Insulator Interfaces\abst
The suppression of rectification at metal–Mott-insulator interfaces, which was previously shown by numerical solutions to the time-dependent Schrödinger equation and experiments on real devices, is reinvestigated theoretically using nonequilibrium Green’s functions. The one-dimensional Hubbard model is used for a Mott insulator. The effects of attached metallic electrodes are incorporated into the self-energy. A scalar potential originating from work-function differences and satisfying the Poisson equation is added to the model. For electron density, we decompose it into three parts. One is obtained by integrating the local density of states over energy to the midpoint of the electrodes’ chemical potentials. The others, obtained by integrating lesser Green’s functions, are due to couplings with the electrodes and correspond to an inflow and an outflow of electrons. In Mott insulators, incoming electrons and holes are extended over the whole system, avoiding further accumulation of charges relative to that in the case without bias. This induces collective charge transport and results in the suppression of rectification. \kword metal-insulator interface, rectification, Poisson equation, nonequilibrium Green’s function, Mott insulator
Correlated electron systems can be candidate materials for novel functions of electronic devices. Most theories for electronic devices are, however, based on a conventional one-electron picture for band semiconductors or band insulators. For correlated electron systems, such a picture may need to be replaced. Devices are always surrounded by interfaces, where the match or mismatch of work functions or bands  must be considered in an appropriate manner. Generally, where two materials with different work functions are attached, bands are bent. A Schottky barrier is formed at a metal-insulator interface, [2, 3] which is governed by the long-range Coulomb interaction. Because the barrier height is modulated by external bias, the current is usually not an odd function of external bias. The application of forward (reverse) voltage lowers (raises) the interfacial barrier, leading to a larger (smaller) current, i.e., rectifying action at the metal–band-insulator interface. 
The importance of electron correlation in charge transport through metal-insulator interfaces is recognized in metal-insulator-semiconductor field-effect transistor (MISFET) device structures based on organic single crystals of the quasi-one-dimensional Mott insulator (BEDT-TTF)(FTCNQ) [BEDT-TTF=bis(ethylenedithio)tetrathiafulvalene, FTCNQ=2,5-difluorotetracyanoquinodimethane].  In Mott insulators, the carrier injections are ambipolar even if the work function of the crystal is quite different from that of the electrodes. Such characteristics are reproduced in the one-dimensional Hubbard model attached to a tight-binding model, where the formation of Schottky barriers is taken into account by added potentials satisfying the Poisson equation.  A connection was suggested between these ambipolar field-effect characteristics and the suppression of rectification at metal–Mott-insulator interfaces. 
We have shown that rectification at metal–Mott-insulator interfaces is indeed suppressed, compared with rectification at metal–band-insulator interfaces, even for large work-function differences.  This fact is demonstrated by numerical solutions to the time-dependent Schrödinger equation and also by experiments on real devices made of organic crystals of (BEDT-TTF)(FTCNQ). With the current-voltage characteristics = of the metal-insulator interface, it is shown that the drain current of the field-effect transistor is approximately proportional to for small ; denotes the drain voltage, the gate voltage in the symmetric-gate operation,  and is the derivative of . Thus, the ambipolar field-effect characteristics [ an even function of ] are the consequence of the suppressed rectification [ an odd function of ] at the metal–Mott-insulator interfaces.
In the above-mentioned theoretical studies, after the ground state is obtained, a finite voltage is suddenly applied and maintained at a constant value. Because we employ the periodic boundary condition for a finite system (with appropriate choice of a gauge), the current under a constant bias finally oscillates owing to the finite-size effect.  Thus, the current density is estimated by averaging the time-dependent one over a given time period. The systems that can be treated by this method are limited to those systems in which the band structure of the left electrode coincides with that of the right electrode by shifting a constant energy. Furthermore, the time-averaged quantities are generally different from the corresponding quantities in the steady state.
In order to avoid such artifacts, we need to treat steady states, which are free of the limitation on the metallic electrodes and the current oscillation due to the finite-seize effect. Without a heat bath, a steady state is reached by couplings with infinitely large, metallic electrodes. In studying the current-voltage characteristics in such a steady state, nonequilibrium Green’s functions are common tools. Because we employ metallic electrodes where electrons are noninteracting, the effects of the electrodes can be incorporated into self-energies. [10, 11, 12] They allow us to discuss the characteristics in terms of the electron density distribution.
In this study, we employ them to reinvestigate the current-voltage characteristics caused by metal–Mott-insulator interfaces. The suppression of rectification at these interfaces can be elucidated by observing the spatial dependence of a “nonequilibrium” part of electron density. We will reformulate the Green’s-function method in such a way that we distinguish an inflow of electrons from an outflow of electrons in the “nonequilibrium” part of electron density, which is crucial to the interpretation of numerical results. We will show that, in Mott insulators, incoming electrons and holes are extended over the whole system, giving rise to a collective charge transport responsible for the suppressed rectification.
2 Model and Method
We consider an insulator, to which metallic electrodes are attached on the left and right sides. For this central part, we use the one-dimensional Hubbard model (with on-site repulsion, 0) for a Mott insulator and the one-dimensional tight-binding model with alternating transfer integrals (0) for a band insulator, both at half filling:
where () creates (annihilates) an electron with spin at site , , and . is the total number of sites in the central part. The scalar potential is introduced above to account for the redistribution of electrons at interfaces to form barriers, compensating for the work-function differences and in equilibrium.  The applied voltage is defined such that it is positive when the right electrode has a lower potential (for electrons) than the left and the current (without multiplication of charge) flows to the right.
Although the potential is defined on lattice points, we solve the Poisson equation in the continuum space,
where the potential and the expectation value of the electron density per site are functions of , and comes from the long-range Coulomb interaction. In order to match the Fermi levels, we set the boundary condition, i.e., the potentials in the metallic electrodes, as
In order to solve the Poisson equation analytically, we assume with a constant compressibility , as in ref. \citenyonemitsu_prb07b where analytic formulas are shown. There are alternative approaches tested, as explained in Appendix A. We have confirmed using the expectation value obtained self-consistently that the current-voltage characteristics and the charge distributions are qualitatively unchanged.
The effects of the left and right (=, ) electrodes consisting of noninteracting electrons on the central part [eq. (1)] are described generally by retarded self-energies. In the wide-band limit, they are independent of energy and their matrix elements with the site indices and are given by 
where =1 for = and 0 otherwise, =1 denotes the site connected to the left electrode, and = denotes the site connected to the right electrode. Within the Hartree-Fock approximation, the retarded Green’s function for spin , , is given by
where the diagonal elements of are given by
with =, and the off-diagonal elements are
and =0 for .
We numerically solve the eigenvalue equation 
where ( and are real) is an eigenvalue of the complex symmetric matrix corresponding to the right eigenvector . We can set up the eigenvalue equation for the adjoint matrix 
where is the Hermitian conjugate of . The ’s and ’s are not identical because and are non-Hermitian matrices. Because and are symmetric matrices, we can take =. In other words, the left eigenvector is given by the complex conjugate of the right one. They are functions of densities (=1, , ; =, ) and must be determined self-consistently. The eigenvectors are normalized according to
For the numerical solutions presented below, we can always find a complete set of eigenvectors satisfying
Then, the retarded Green’s function is written as
For the density , we modify the frequently used formula in ref. \citenbrandbyge_prb02 so as to maintain symmetry concerning simultaneous particle-hole transformation () and space-inversion operation () for = and +=0 ( is assumed to be an even number). The modification is also useful for distinguishing between an inflow and an outflow of electrons, as will be shown later. To be more precise, we decompose the density into the “equilibrium” part and the parts due to the couplings with the left and right electrodes:
where the “equilibrium” part is defined by integrating the local density of states over energy to the midpoint of the left and right chemical potentials, . When the left (right) chemical potential is higher, () corresponds to the inflow, and the other corresponds to the outflow. The “equilibrium” part is then written as
where is the Fermi distribution function in a virtual system with the chemical potential . Considering zero temperature, we substitute the step function for the Fermi distribution function: . The present decomposition is general and useful irrespective of whether or not the wide-band limit is applied. In general cases, the energy dependence of should be kept below.  The definition in ref. \citenbrandbyge_prb02 corresponds to the setting of at either or .
For the “nonequilibrium” part of the density, remember that the lesser self-energy is given in the wide-band limit by 
with . It can be decomposed again into the “equilibrium” part and the parts due to the couplings with the left and right electrodes:
Because the contribution from the “equilibrium” part is regarded as included in , the “nonequilibrium” part of the density, , can be defined as
where the “nonequilibrium” part of the lesser Green’s function, , is given by the Keldysh equation
with in eq. (19) and being the Hermitian conjugate of .
When we substitute eq. (22) into eq. (15), the second term of eq. (22) gives a logarithmic term. This is an artifact of the energy-independent imaginary parts of the eigenvalues in eq. (9) (i.e., the wide-band limit), so that it is ignored below. The first term of eq. (22) gives
The terms in the bracket above are reduced to the step function in the limit of 0 (i.e., 0).
Finally, by using the formula in ref. \citenhaug_book08 with ==1, the current from the left electrode is given by
Therefore, the current is expressed by the “nonequilibrium” parts of the density. We need to substitute =, =, and =0 into all the equations above.
Before showing the density distributions of electrons, we compare the current-voltage characteristics directly obtained by time evolution through a constant potential difference with those obtained by the nonequilibrium Green’s functions in the wide-band limit. As in our previous study,  we have numerically solved the time-dependent Schrödinger equation for the insulator sandwiched between two metallic electrodes with different work functions and a common bandwidth. Figure 1(a) shows the current-voltage characteristics thus obtained for the Mott insulator.
The parameters are introduced and the current density is defined in ref. \citenyonemitsu_prb07b and . The rectifying action is shown to be suppressed. Because the bandwidths of the metallic electrodes are finite, the current density tends to become saturated for a large .
Figure 1(b) shows the current-voltage characteristics obtained by the present approach with the nonequilibrium Green’s functions for the Mott insulator. The rectifying action is shown to be suppressed. Because the bandwidths of the metallic electrodes are assumed to be infinite, no saturation is observed. The values of the other parameters , , , , , and are the same as those in Fig. 1(a), and is close to that used in Fig. 1(a). However, the finite metallic electrodes in Fig. 1(a) are replaced by infinite ones in Fig. 1(b), making direct comparison difficult. As a general trend, the current-voltage characteristics directly obtained by time evolution show smoother curves possibly because the current density is averaged over a given time period. As the system size increases, the present approach gives smoother characteristics, as shown in Fig. 2. In any case, the suppression of rectification for Mott insulators is observed by both of the methods discussed above in a wide parameter space spanned by the coupling strength , the system size , the Coulomb parameter , and the work-function differences and .
Figure 2 shows the current , the “nonequilibrium” part of the density due to the coupling with the left electrode at the rightmost site , and that due to the coupling with the right electrode at the leftmost site , which are related by in eq. (25).
In the antisymmetric case of = and = (not shown), the present approach guarantees the symmetry concerning simultaneous particle-hole transformation and space-inversion operation, which leads to = for any . In a more general case of but =, Fig. 2 shows that the approximate relation still holds for only Mott insulators with and used here. In the band insulator, the absolute value of is larger for [Fig. 2(b)], where the barrier is lower than that for , as expected. Later in Fig. 8, we will see how the rectification is realized in the band insulator.
From now on, unless otherwise stated, we compare Mott and band insulators with a gap =0.1, which is much smaller than the values used in previous studies. [6, 7, 8] If we numerically solved the time-dependent Schrödinger equation as before for such a small gap, we would need to adopt a small to calculate the evolution for a long time, which is proportional to . Here, we use =200 (larger than before) and =0.03 (slightly smaller than before), which would also need long calculations in the previous approach. Because the present study focuses on the density distributions of incoming and outgoing electrons, we fix the parameter set to show them. Unless otherwise stated, the work function of the left (right) electrode is set to match the bottom (top) of the upper (lower) Hubbard or conduction (valence) band, ==. Thus, the barriers at the two interfaces almost or completely disappear for the right-going bias (0), and they become prominent for the left-going bias (0), as shown in Fig. 3.
The absolute values of are large for 0 and small for 0 for the band insulator.
Observing the charge distribution will help to understand the mechanism for making such a characteristic difference between Mott and band insulators. Figure 4(a) shows the charge density, =, and its “equilibrium” part, =, for a Mott insulator with a left-going bias and barriers.
Electrons accumulate near the left interface, while holes accumulate near the right interface. This is because the bands are bent near these interfaces (Fig. 3). It should be noted that electrons and holes already accumulate for =0, which is responsible for the respective band bending at the interfaces and the resultant matching of the chemical potentials to reach equilibrium from the isolated case with =0. By giving a finite bias, a certain number of electrons come in and the same number of electrons (holes) go out (come in) in the steady state, so that the incoming electrons or holes cause no instability. If one takes a closer look, the charge density alternates between even and odd sites, although its amplitude is very small. This 2 oscillation is induced by the boundaries.
Figure 4(b) shows the “nonequilibrium” parts of the charge density, one of which is due to the coupling with the left electrode, =, and the other is due to that with the right electrode, =. Here, the voltage is negative so that electrons are left-going on average. Because of the coupling with the left electrode, electrons go out and 0. Owing to the coupling with the right electrode, electrons come in and 0. Thus, the present definition of is useful for distinguishing between these flows. The property unique to Mott insulators is that is extended over the whole system in such a manner that is almost constant except the 2 oscillation. Mott insulators disfavor further accumulation of charges relative to that in the case of =0, so that incoming electrons and holes are extended over the system.
Figure 5 shows , , , and for a Mott insulator with a right-going bias.
For this particular set of (, , )=(0.05, 0.05, 0.1), the potential is constant everywhere [eq. (3)] and the density is almost unity (note the vertical scale) except for the 2 oscillation [eq. (2)]. Here, the voltage is positive so that electrons are right-going on average. Electrons come in from the left, 0, and go out to the right, 0. Both and are almost constant except for the small 2 oscillation. In Mott insulators, the averaged (i.e., 2-oscillation smoothed out) densities of the incoming electrons and holes show not only a small spatial modulation but also insensitivity to the sign of . This fact ensures the suppression of rectification because the current is given by the difference between the density of the incoming electrons at the exit and that of the outgoing electrons at the entrance [multiplied by the respective , eq. (25)].
The behavior of the metal–Mott-insulator interfaces is contrasted with that of the metal–band-insulator interfaces. Figure 6 shows , , , and for a band insulator with a left-going bias and barriers.
It is clearly shown in that electrons accumulate near the left interface, while holes accumulate near the right interface [Fig. 6(a)]. They are much larger than the corresponding quantities at the metal–Mott-insulator interfaces. Thus, the accumulation is more sensitive to the band bending.
The difference from the metal–Mott-insulator interfaces is conspicuous in and [Fig. 6(b)]. Here, electrons are left-going on average. The quantity is thus negative, but its magnitude becomes almost zero at some point. On its left side, is largest at the left interface and decays with increasing distance from the interface. This behavior originates from the outgoing electrons. On its right side, shows a maximum away from the right interface. This behavior can be regarded as due to , which is higher on the right side. In other words, holes further accumulate on the right side at a negative . The spatial dependence of is obtained by the particle-hole transformation and the space-inversion operation of . The quantity is positive and becomes almost zero at another point. On its right side, is largest at the right interface and decays with increasing distance from the interface. This behavior originates from the incoming electrons. On its left side, shows a maximum away from the left interface owing to further accumulation of electrons at lower ’s. This further accumulation is realized by the absence of on-site repulsion.
Figure 7 shows , , , and for a band insulator with a right-going bias.
For this particular set of (, , ), the potential is constant. The quantity is unity, but largely deviates from it. Because the density of electrons is much more sensitive to the bias-induced change in the potential distribution for band insulators, and are now quite different. Here, electrons are right-going on average, so that 0 and 0. Because of the absence of barriers, a substantial number of electrons come in and go out.
In order to see how the rectification is realized in a band insulator, we show its and in Fig. 8 with parameters used in Fig. 2(b) and different biases, whose magnitudes are larger than the gap =0.18 here.
For , the absolute value of , , is smaller than that for [Fig. 2(b)]. In this case, the density of the incoming electrons decays with increasing distance from the right interface, while the density of the outgoing electrons decays with increasing distance from the left interface. With increasing , both and slightly increase at each site , but they decay in quite similar manners [Fig. 8(a)].
For , on the other hand, increases more rapidly with than for . Now, the current substantially flows through the system. The density of the outgoing electrons decays with increasing distance from the right interface. However, with increasing , a substantial number of electrons are shown to penetrate from the left interface into the system [Fig. 8(b)]. Both and finally decay toward the right and left interfaces, respectively, giving small and values compared with and , respectively. However, for is much larger than and for because of the penetration of a larger number of electrons into the system. These behaviors of and in a band insulator are in contrast to those in a Mott insulator, where the averaged (i.e., 2-oscillation smoothed out) densities of the incoming electrons and holes show not only little spatial modulation but also insensitivity to the sign of . The dependence of and that of in the Mott insulator closely follow that of and that of in Fig. 2(a), respectively.
4 Summary and Conclusions
To elucidate the suppression of rectification at metal–Mott-insulator interfaces, which is obtained by numerical solutions to the time-dependent Schrödinger equation in our previous study,  we employ nonequilibrium Green’s functions in this paper. We consider one-dimensional half-filled electron systems and use the mean-field Hubbard model for a Mott insulator and the transfer-alternating tight-binding model for a band insulator. Metallic electrodes consisting of noninteracting electrons are attached to the system, and their effects are incorporated into self-energies within the wide-band limit. We take account of work-function differences, which are responsible for the formation of Schottky barriers, by adding to the model a scalar potential that satisfies the Poisson equation.
In the mean-field approximation, the retarded and advanced Green’s functions are obtained by solving the eigenvalue equation for a complex symmetric matrix, whose imaginary part comes from the self-energies due to the couplings with the electrodes. For the electron density to be determined self-consistently, we modify the frequently used formula in ref. \citenbrandbyge_prb02 so as to maintain some symmetry in the particular case used here. The modification is useful for distinguishing between an inflow and an outflow in the “nonequilibrium” part of the density. The current is given by the difference between these flows, which are measured at appropriate sites and multiplied by couplings with electrodes.
By plotting the spatial dependence of the density of incoming electrons and holes, we clarify the difference between metal–Mott-insulator and metal–band-insulator interfaces. In Mott insulators, the incoming electrons and holes are extended over the whole system in such a manner that and are almost constant over =1, , . Electrons/holes do not accumulate at any place, so that the bias-induced change in the density of electrons is extended. Thus, charge transport becomes collective.
This work was supported by Grants-in-Aid for Scientific Research (C) (No. 19540381), for Scientific Research (B) (No. 20340101), and “Grand Challenges in Next-Generation Integrated Nanoscience” from the Ministry of Education, Culture, Sports, Science and Technology.
Appendix A Alternative Approach to Band Bending
As is well known in electromagnetics, the Poisson equation is equivalent to the long-range Coulomb interaction, , where is the position vector for site . Then, we can define the Hartree potential as
and add a linear term to it as
where the constants and are determined so that the potential satisfies the boundary condition [eq. (3)]. In this way, we can treat the band bending at interfaces and the effect of the long-range Coulomb interaction on the insulator (e.g., charge ordering if it is strong) on the same footing. This method can be extended to include the screening effect. In addition, we found that its numerical convergence is generally better than that of solving the Poisson equation.
Appendix B Densities in the Atomic Limit
In the limit of 0, i.e., in the atomic limit, one can simplify the expressions for the density. The matrix is diagonal and its elements are given by , so that the eigenvalues are written as for = and as otherwise. Then, the “equilibrium” and “nonequilibrium” parts of the density at = are then written as
The total density is thus written as
which is a reasonable expression in the atomic limit.
-  Y. Takahashi, T. Hasegawa, Y. Abe, Y. Tokura, and G. Saito: Appl. Phys. Lett. 88 (2006) 073504.
-  W. Schottky: Naturwissenschaften 26 (1938) 843.
-  N. F. Mott: Proc. Cambridge Philos. Soc. 34 (1938) 568.
-  S. M. Sze and K. K. Ng: Physics of Semiconductor Devices (John Wiley & Sons, New Jersey, 2007) 3rd ed.
-  T. Hasegawa, K. Mattenberger, J. Takeya, and B. Batlogg: Phys. Rev. B 69 (2004) 245115.
-  K. Yonemitsu: J. Phys. Soc. Jpn. 74 (2005) 2544.
-  K. Yonemitsu: inMultifunctional Conducting Molecular Materials, ed. G. Saito, F. Wudl, R. C. Haddon, K. Tanigaki, T. Enoki, H. E. Katz, and M. Maesato (RSC Publishing, Cambridge, 2007), p. 276.
-  K. Yonemitsu, N. Maeshima, and T. Hasegawa: Phys. Rev. B 76 (2007) 235118.
-  T. Oka, R. Arita, and H. Aoki: Phys. Rev. Lett. 91 (2003) 066406.
-  Y. Meir and N. S. Wingreen: Phys. Rev. Lett. 68 (1992) 2512.
-  N. S. Wingreen, A.-P. Jauho, and Y. Meir: Phys. Rev. B 48 (1993) 8487.
-  A.-P. Jauho, N. S. Wingreen, and Y. Meir: Phys. Rev. B 50 (1994) 5528.
-  S. Datta: Electronic Transport in Mesoscopic Systems (Cambridge University Press, Cambridge, 1995).
-  M. Brandbyge, J.-L. Mozos, P. Ordejón, J. Taylor, and K. Stokbro: Phys. Rev. B 65 (2002) 165401.
-  H. Haug and A.-P. Jauho: Quantum Kinetics in Transport and Optics of Semiconductors (Springer, Berlin, 2008) 2nd ed.