^{1}, Feng Jiang

^{2}, Heng Tian

^{1}, Xiao Zheng

^{3}, Yanho Kwok

^{1}, Shuguang Chen

^{1}, ChiYung Yam

^{1}, YiJing Yan

^{2,3}and Guanhua Chen

^{1,a)}

### Abstract

Basing on our hierarchical equations of motion for time-dependent quantum transport[X. Zheng, G. H. Chen, Y. Mo, S. K. Koo, H. Tian, C. Y. Yam, and Y. J. Yan, J. Chem. Phys.133, 114101 (2010)10.1063/1.3475566], we develop an efficient and accurate numerical algorithm to solve the Liouville-von-Neumann equation. We solve the real-time evolution of the reduced single-electron density matrix at the tight-binding level. Calculations are carried out to simulate the transient current through a linear chain of atoms, with each represented by a single orbital. The self-energy matrix is expanded in terms of multiple Lorentzian functions, and the Fermi distribution function is evaluated via the Padè spectrum decomposition. This Lorentzian-Padè decomposition scheme is employed to simulate the transient current. With sufficient Lorentzian functions used to fit the self-energy matrices, we show that the lead spectral function and the dynamics response can be treated accurately. Compared to the conventional master equation approaches, our method is much more efficient as the computational time scales cubically with the system size and linearly with the simulation time. As a result, the simulations of the transient currents through systems containing up to one hundred of atoms have been carried out. As density functional theory is also an effective one-particle theory, the Lorentzian-Padè decomposition scheme developed here can be generalized for first-principles simulation of realistic systems.

Authors would like to thank Yu Zhang and Hui Cao for many helpful discussions on transport theory and spectral integral. The authors also give thanks to Dr. Jian Sun for some help with computer. Support from the Hong Kong Research Grant Council (HKU700808P, HKU700909P, HKU700711P, HKUST9/CRF/08), AOE (AOE/P-04/08), National Science Foundation of China (NSFC) (21103157, 21033008), and Fundamental Research Funds for the Central Universities of China (2340000034) is gratefully acknowledged.

I. INTRODUCTION

II. THEORY

A. Hierarchical equation of motions for reduced single-electron density matrix

B. Self-energy decomposition method

C. The central equations of the HEOM method

III. NUMERICAL IMPLEMENTATION

A. Matrix-based Lorentzian expansion

B. Initial values and numerical integration in time domain

IV. RESULTS

A. Numerical test on a single level system

B. Simulation of transient current through a chain of atoms

C. Simulation of one-dimensional atom chain with different contacts

D. Tight-binding model with the next nearest neighbor interactions and asymmetric left-and-right transient current

V. NUMERICAL EFFICIENCY

VI. CONCLUSION

### Key Topics

- Lead
- 34.0
- Linewidths
- 25.0
- Green's function methods
- 10.0
- Density functional theory
- 7.0
- Differential equations
- 6.0

## Figures

Contours for the integral of linewidth function. The upper contour (solid line, denoted as ) and lower contour (dotted line, denoted as ) are used for different signs of τ − *t*. The filled circles (on y axis) and hollow circles represent the Padè poles (μ_{α} = 0) and Lorentzian poles, respectively. The arrows indicate the contour integral directions.

Contours for the integral of linewidth function. The upper contour (solid line, denoted as ) and lower contour (dotted line, denoted as ) are used for different signs of τ − *t*. The filled circles (on y axis) and hollow circles represent the Padè poles (μ_{α} = 0) and Lorentzian poles, respectively. The arrows indicate the contour integral directions.

(Left) Transient currents for single level system with different W values; (Right) Transient currents for the same system reported in Ref. 34.

(Left) Transient currents for single level system with different W values; (Right) Transient currents for the same system reported in Ref. 34.

Left: Lorentzian fitting curves (dotted line for 4 Lorentzian fitting and dashed line for 9 Lorentzian fitting) and the accurate curve (solid line) for the dimensionless linewidth function (); Right: transmission spectrum of 4-site chain TB system calculated with the 4-Lorentzian fitted self-energy (dotted line) and 9-Lorentzian fitted self-energy (solid line).

Left: Lorentzian fitting curves (dotted line for 4 Lorentzian fitting and dashed line for 9 Lorentzian fitting) and the accurate curve (solid line) for the dimensionless linewidth function (); Right: transmission spectrum of 4-site chain TB system calculated with the 4-Lorentzian fitted self-energy (dotted line) and 9-Lorentzian fitted self-energy (solid line).

Time-dependent current for a 4-site 1D system with two linewidth fitting schemes: 4-Lorentzian (dashed line) and 9-Lorentzian (solid line) fitting. The inset is the magnified figure. The coupling constant: in the leads and device, *t* _{0} = 2 eV; between leads and device *h* = 0.9*t* _{0}. Temperature is 300 K, bias voltage is 0.01 V, 10 Padè points are used.

Time-dependent current for a 4-site 1D system with two linewidth fitting schemes: 4-Lorentzian (dashed line) and 9-Lorentzian (solid line) fitting. The inset is the magnified figure. The coupling constant: in the leads and device, *t* _{0} = 2 eV; between leads and device *h* = 0.9*t* _{0}. Temperature is 300 K, bias voltage is 0.01 V, 10 Padè points are used.

I-V curve (left), DOS spectrum (middle), and the transient current (right) for a 2-atom device with good contact (upper row) and poor contact (lower row). The bias voltage is 2 V, the coupling constants *t* in leads and device are 2 eV; and *t* from leads to device is 1.8 eV (good contact) or 0.4 eV (poor contact).

I-V curve (left), DOS spectrum (middle), and the transient current (right) for a 2-atom device with good contact (upper row) and poor contact (lower row). The bias voltage is 2 V, the coupling constants *t* in leads and device are 2 eV; and *t* from leads to device is 1.8 eV (good contact) or 0.4 eV (poor contact).

The linewidth functions with the lead Green's functions obtained from the accurate iteration method (a) or from the Lorentzian fitting (b) scheme. (Λ_{12} and Λ_{21}correspond the same curves.) The system is a 4-site TB model with the next nearest neighbor interactions. The parameters: t = 2 eV; h_{2} = 0.4 eV. (c) The transient current of the system under a step bias voltage (2 V) symmetrically applied on two leads at time = 0. The dotted line is the mirror (opposite sign) of the dashed line (right current).

The linewidth functions with the lead Green's functions obtained from the accurate iteration method (a) or from the Lorentzian fitting (b) scheme. (Λ_{12} and Λ_{21}correspond the same curves.) The system is a 4-site TB model with the next nearest neighbor interactions. The parameters: t = 2 eV; h_{2} = 0.4 eV. (c) The transient current of the system under a step bias voltage (2 V) symmetrically applied on two leads at time = 0. The dotted line is the mirror (opposite sign) of the dashed line (right current).

CPU time dependence on the simulation time (left) and site number N (right). In the left figure, the system is a 1D TB chain with 20 sites, 4 Lorentzians and 10 Padè points. And in the right figure, the simulation time is fixed at 0.1 fs, and the time step is 0.002 fs for these two cases.

CPU time dependence on the simulation time (left) and site number N (right). In the left figure, the system is a 1D TB chain with 20 sites, 4 Lorentzians and 10 Padè points. And in the right figure, the simulation time is fixed at 0.1 fs, and the time step is 0.002 fs for these two cases.

Article metrics loading...

Full text loading...

Commenting has been disabled for this content