Contour integral method for obtaining the self-energy matrices of electrodes in electron transport calculations

We propose an efficient computational method for evaluating the self-energy matrices of electrodes to study ballistic electron transport properties in nanoscale systems. To reduce the high computational cost incurred in large systems, a contour integral eigensolver based on the Sakurai-Sugiura method combined with the shifted biconjugate gradient method is developed to solve an exponential-type eigenvalue problem for complex wave vectors. A remarkable feature of the proposed algorithm is that the numerical procedure is very similar to that of conventional band structure calculations. We implement the developed method in the framework of the real-space higher-order finite-difference scheme with nonlocal pseudopotentials. Numerical tests for a wide variety of materials validate the robustness, accuracy, and efficiency of the proposed method. As an illustration of the method, we present the electron transport property of the freestanding silicene with the line defect originating from the reversed buckled phases.

Figure 1

Schematic representation of a quasi-one-dimensional conductor sandwiched between two semi-infinite electrodes. Because of the localized feature of the real-space grid approach, semi-infinite crystalline electrodes can be divided into principal layers denoted by $L0,L1,...,Lm;R0,R1,...,Rm$. $L0$ and $R0$ are incorporated into the transition region as WFM planes.

Figure 2

Two equivalent contours in the complex $\lambda$ plane and $k$ plane. (a) A multiply connected domain between two concentric circles at the origin of radii ${\lambda }_{min}$ and ${\lambda }_{min}^{-1}$ in the $\lambda$ plane. (b) A simply connected domain obtained from (a) by variable conversion, $k=ln\lambda /ia$. The contour integral proceeds in a counterclockwise direction along a rectangular side path. For this type of contour, quadrature points are specified by a pair of polynomial orders of the Gauss-Legendre quadrature rule, that is, $N_q=(N_{q1},N_{q2})$. In the figure, there are $N_{q1}=6$ Gauss-Legendre quadrature points along the $\mathrm{Re}\left(k\right)$ axis and $N_{q2}=6$ points along the $\mathrm{Im}\left(k\right)$ axis. The total number of quadrature points is $N_{\mathrm{int}}=2N_{q1}+2N_{q2}$. The quadrature points are indicated by $\circ$ and $•$; however, the numerical integration can only be performed using the values indicated by $•$ owing to the symmetry of the Hamiltonian matrix $H\left(k\right)$.

Figure 3

Numerical results for fcc Au bulk at the Fermi energy of (a) distribution of eigenvalues within the domain enclosed by $\mathrm{\Gamma }$ and (b) residuals $\left|\right|[E-H\left(k_n\right)]{\varphi }_n{\left|\right|}_2$ when varying the order of the Gauss-Legendre quadrature rule, $N_q=(N_{q1},N_{q2})$. The number of target eigenvalues that do not include spurious eigenpairs is 54. The plots clearly show that the positions of the eigenpairs are almost unchanged and that the accuracy is straightforwardly improved by increasing $N_q$. Convergence is achieved at $N_q=(24,24)$: all target eigenvalues are below the convergence criterion, as indicated by the broken line.

Figure 4

Residuals $\left|\right|[E-H\left(k_n\right)]{\varphi }_n{\left|\right|}_2$ for (6,6)CNT at the Fermi energy calculated on the $k$ plane (EEP/SS) and $\lambda$ plane (QEP/SS). (a) For ${\lambda }_{min}=0.1$, the number of eigenvalues that do not include spurious eigenpairs is 52. (b) For ${\lambda }_{min}=0.01$, the number of eigenvalues that do not include spurious eigenpairs is 154. In both cases, the EEP/SS method uses $N_q=(24,24)$ as the order of the Gauss-Legendre quadrature rule; by contrast, the QEP/SS method uses the trapezoidal rule with the number of quadrature points being 256 per circle in Fig. 2. The other parameters are kept the same.

Figure 5

(a) Transition region of Au atomic chain with CO adsorption. Au, C, and O atoms are represented as gold, brown, and red balls, respectively. (b) Transmission spectra obtained using self-energy matrices calculated by the proposed method with four different ${\lambda }_{min}$ values: 0.999 (red line), 0.1 (blue line), 0.01 (green line), and 0.001 (black line). For clarity, transmission spectra are shifted by the amount of 2.5 with respect to the original values in descending order of the legend. (c) The transmission spectrum obtained with the proposed method (black line). The results obtained with the maximally localized Wannier function (WF) of Ref. [16] (red dashed line) and Ref. [58] (blue dashed line) are also shown in (c).

Figure 6

Differences between transmission spectra using approximated and exact self-energy matrices. The overall deviation is indicated in the parentheses following the legends.

Figure 7

Optimized interface structure of silicene with $\alpha -\beta$ interface.

Figure 8

Effects of the $\alpha -\beta$ interface on the transmission spectra. The solid line is the result of real-space grid calculation using the self-energy matrices obtained with ${\lambda }_{min}=0.01$. Empty dots denote the transmission spectrum without the defect.

Figure 9

(a) Band structure and (b) group velocity of silicene.

Figure 10

Charge densities of two bulk modes in III band for left electrode at $E=E_F+0.6$ and 0.91 eV.

Figure 11

An illustration of the scattering of silicene at $E=E_F+0.91$ eV. The symbols L, M, and U represent the lower buckled, nonbuckled, and upper buckled atoms, respectively.

Figure 12

An illustration of a one-dimensional system with square potential barriers. The parameters $a,b,V_0$, and $V_1$ represent the width of depths, width of barriers, barrier height in L and R regions, and barrier height in C region, respectively.

Figure 13

(a) Energy dispersion of the Kronig-Penny model with periodic square potentials and (b) transmission spectra. The parameters in atomic units are set as $a=2.0,b=0.2,V_0=10$, and $V_1=11$.