Atomistic Full-quantum and Full-band Transport Model for Zigzag Group-IVA Nanoribbon-based Structures with Noniterative Calculation Framework

1Department of Refrigeration, Air Conditioning and Energy Engineering, National Chin-Yi University of Technology, Taichung 411, Taiwan 2Department of Refrigeration & Air Conditioning and Energy, Far East University, Tainan 74448, Taiwan 3Department of Electrical Engineering, Far East University, Tainan 74448, Taiwan 4Quantum Engineering Laboratory, Department of Physics, Tamkang University, Tamsui, New Taipei 25137, Taiwan


Introduction
The graphene, which is regarded as a semimetal or zero-gap semiconductor, possesses a twodimensional (2D) honeycomb structure of carbon atoms, and it has attracted tremendous interest owing to its unique properties. (1)(2)(3) Therefore, graphene is considered as a promising material for future research. However, its zero energy gap has been widely regarded as an obstacle to its application in next-generation nanoelectronics. Therefore, material researchers have tried to solve the zero-gap problem by using new 2D honeycomb materials from other elements in group-IVA of the periodic table. Firstly, many investigations have shown that 2D silicene (4,5) and germanene, (4,5) which are the graphene-like analogues of silicon and germanium, respectively, have promising applications in nanoscale electronic devices. The last research in this field is related to the 2D honeycomb lattice of a Sn monolayer, commonly referred to as stanene. (4,6) Graphene possesses a planar 2D honeycomb lattice of carbon atoms with π-bonding. In contrast to graphene, other 2D layered group-IVA materials possess a buckled structure in the honeycomb lattice, which leads to partial sp 3 hybridization of an electronic configuration. (4) The buckling of a honeycomb lattice enhances the hybridization between π and σ orbitals. (6) Notably, owing to the enhancement of sp 3 hybridization with increasing group-IVA atomic number, out-of-plane buckling of the honeycomb lattice is enhanced. Buckled 2D materials possess a stronger spin-orbit coupling (SOC) interaction, which gives rise to a band gap between conduction and valence bands. (6,7) A small band gap is one of the most important differences of silicene and germanene from graphene. (8) For stanene, its relatively large buckling and thus its stronger SOC interaction open an obvious energy gap in the energy band diagram, (7) which solves the zero-gap problem for planar 2D nanostructures.
The nonequilibrium Green's function (NEGF) (9)(10)(11)(12)(13) method is the most powerful method for solving the quantum transport problems of nanoscale electronic devices. In the past, many studies have shown that the NEGF can be solved by a recursive (or iterative) technique. (14)(15)(16)(17)(18) A recursive technique may finish this job or not (divergence), while one of the disadvantages of this technique is the slow convergence in the iterative steps. Moreover, another technique for obtaining the NEGF is based on the calculation framework of the Dyson equation, which was proposed by Caroli et al. (19) and others. (20)(21)(22) One of the disadvantages of the Dyson equation is the difficulty to dispose the miscellaneous boundary conditions. However, the disadvantages of the recursive and Dyson techniques can be overcome by the method proposed in this paper, which a complex energy-band framework in association with the NEGF approach. (23)(24)(25)(26) In this paper, the properties of zigzag-edged group-IVA nanoribbon (z4ANR) structures such as the density of states (DOS), transmission coefficient, and conductance will be explored by the proposed technique.

Theoretical Methods
The z4ANR sample considered in this paper possesses N zigzag lines and has a central channel region composed of l atomic layers denoted by σ = 1, 2, ..., l, as shown in Fig. 1. It is supposed that flat-band cases occur in the left-side incoming (σ ≤ 0) and right-side outgoing (σ ≥ l + 1) regions, which are located outside of the central channel region.
In this study, the theoretical model 4 of a monolayer group-IVA element structure is a tightbinding calculation based on the τ β α (α = s, x, y, z) orbital, where β is an A or B sublattice, and the electron spin τ is along the up (↑) or down (↓) direction. Therefore, the Hamiltonian of these monolayer group-IVA elements can be expressed in 16 × 16 matrix form with the τ β α orbital, which is presented in Appendix A.

Hamiltonian matrix and state function of z4ANRs
In the flat-band case, the state function k ⊥ of a z4ANR, which has N zigzag lines, is a combination of the 16N tight-binding basis | , , k j τ β α ⊥ >, which can be written as , , , , where b j,α,β,τ specifies the expansion coefficient, τ β α denotes the α orbital of the group-IVA element with electron spin τ (= ↑ or ↓) located at sublattice β (= A or B), j is the index of the N adjacent zigzag lines in the transverse direction (||) of the z4ANR, and k ⊥ is the wave vector directed along the z4ANR direction (⊥). The electron energy E, which is shown as | , , is omitted for brevity here. Furthermore, the tight-binding basis of a z4ANR Hamiltonian can be expressed as where Ω denotes the normalization number, σ specifies a layer label increasing along the ⊥ direction, a' denotes the spacing of two neighboring layers, and | , , where , , , where H σ,σ and H σ,σ±1 are 16N × 16N matrices (see Appendix C), whose elements are written as , , , , ; , , , respectively.

Complex energy-band structure of z4ANRs
The Schrödinger equation ( ) | 0 H E k ⊥ − > = in a z4ANR structure can be expressed as (27)(28)(29) where I denotes a 16N × 16N identity matrix, a state function k ⊥ specifies the available planewave states in the left and right regions, and c σ can be expressed as a 16N-length column vector of coefficients whose components are denoted as c σ,j,α,β,τ , i.e., According to Bloch's theorem, the tight-binding coefficients must satisfy the relation To solve the eigenvalue ( ik a e ⊥ ′ ) of Eq. (6)

Hamiltonian matrix and wave function of z4ANR-based devices
In the z4ANR-based str ucture, the wave function of the Schrödinger equation where ( ) a k ⊥ denotes the amplitude coefficient of a z4ANR state function k ⊥ . Therefore, the wave function of the whole structure can be expressed in the form , , , , , , where f σ,j,α,β,τ is , , , , and thus we obtain the following combinative equation for the σth layer: (28)(29)(30) , where H σ,σ and H σ,σ±1 can be written in the form of 16N × 16N matrices, as shown in Appendix C, and f σ is a 16N-length column vector whose components are f σ,j,α,β,τ . Therefore, the Hamiltonian matrix of the whole structure can be expressed as

Boundary conditions solved by complex energy-band method
We rearrange the state functions , 1~32N k λ ⊥ = of the complex energy-band formalism, which are obtained from Eq. (6). Therefore, the index of λ = 1, 2, ..., 16N corresponds to the state functions, which either propagate (k ⊥ real) or decay ( complex) to the right-hand side. On the other hand, λ = 16N + 1, 16N + 2, ..., 32N correspond to those which either propagate or decay to the left-hand side. (29,30) The boundary conditions are such that we have a known incoming plane-wave state from the left contact, no incoming plane-wave state from the right contact, and unknown outgoing transmitted and reflected plane-wave states in the right and left contacts, respectively. For a given energy E and for a given amplitude (here unity) of an incoming planewave-like state (denoted by i) from the left, the wave functions in the left (L) and right (R) contacts must satisfy the boundary conditions of this example, which can be expressed in terms of the state functions of the complex energy-band structure as follows: (29,30) 16 , , are the unknown amplitude coefficients of the transmitted and reflected state functions, respectively. For convenience, we use ℜ and ℑ to denote the outgoing waves that propagate (or decay) to the left in the left contact and to the right in the right contact, respectively.

Applicative calculation using NEGF framework
The Green's function of the whole device is expressed as The Green's function G d is obtained, and then the transmission function ( ) T E is determined by the trace of (12,32) where ( ) T E is the product of the number of forward-propagating eigenstates M(E) and the transmission probability T(E), The conductance G(E) of the device is obtained from the Fermi energy E F , and it can be associated with the transmission function ( ) T E as (12,32-34) where h denotes the Planck constant, e denotes the electron charge, and 2e 2 /h represents the conductance quantum. Furthermore, the transmission function ( ) T E should be obtained at the Fermi energy E F using ( ) [ ] If the Green's function G d is obtained, the DOS can be determined via (12,15,32)

Results and Discussion
Figures 2(a)-2(d) show the energy-band structures of monolayer group-IVA elements, which are calculated by the tight-binding four-band technique with the SOC effect in the nearestneighbor approach. Furthermore, the tight-binding four-band results are in good agreement with the first-principles results only for the valence and conduction bands in the low-energy sections, (4) i.e., near the K points. Therefore, the four-band model used here can capture the essential physics of monolayer group-IVA elements in the low-energy region reasonably well.
Graphene possesses a planar (θ = 90°) honeycomb lattice of carbon atoms, where θ is the angle between the bond and the z direction. In contrast to graphene, other layered group-IVA materials possess a buckled (θ ≠ 90°) structure in the honeycomb lattice. θ is 101.7° for silicene, 106.5° for germanene, and 107.1° for stanene. (4) The buckling in the honeycomb lattice enhances the overlap between π and σ orbitals. Therefore, the buckling in the honeycomb lattice leads to partial hybridization between π and σ orbitals, and thus partial sp 3 hybridization exists in the orbital configuration. Furthermore, owing to the increasing sp 3 hybridization with increasing group-IVA atomic number, the out-of-plane buckling of the honeycomb lattice is enhanced. Owing to the low atomic mass of graphene, its intrinsic SOC interaction is negligible. As shown in Appendix A, the strength of the SOC interaction 0 ξ is 0.009 for graphene, 0.034 for silicene, 0.196 for germanene, and 0.8 for stanene in units of eV. (4) Compared with the planar graphene, the buckled 2D group-IVA elements possess a stronger intrinsic SOC interaction, which gives rise to a band gap in the K points. Possessing a small band gap is the most important difference of silicene and germanene from graphene, as shown in Figs. 2(a)-2(c). Compared with silicene and germanene, stanene possesses a larger buckled altitude and thus a stronger SOC interaction, which opens an obvious band gap in the energy-band diagram, as shown in Fig. 2(d).
A z4ANR, which has N zigzag lines (width), possesses 2N group-IVA atoms in the interior of its unit cell, as presented in Fig. 1. A 2N-  E k λ ⊥ = − complex energy-band diagram. (23)(24)(25) Moreover, some value of E and its corresponding real k ⊥ yield the traditional Ek ⊥ energy-band diagram of a z4ANR. For example, the traditional Ek ⊥ energy-band diagram of a zGNR with N (= 4) zigzag lines is presented in the left panel of Fig. 4. For a real value of k ⊥ , the propagating waves transmit to  Fig. 4), where g υ is the group velocity ]. Furthermore, for a complex value of k ⊥ , the evanescent waves decay exponentially in the right [Im( ) 0 k ⊥ > ] or left (Im( ) 0 k ⊥ < ) direction. Each of the eigenvectors of Eq. (6) corresponds to a pair of k ⊥ and k ⊥ − ; hence, half (16N) of the propagating or evanescent waves propagate rightward and the other half (16N) propagate leftward at a certain energy E, as shown in Fig. 4.
In the z4ANR-based devices considered in this section, the width of the channel is N (= 4) zigzag lines and the length of central channel is l (= 14) atomic layers. For these z4ANR-based samples, which assume a perfect flat-band potential profile, Figs. 5(a)-5(d) display the spectral diagrams of the transmission function (T ), conductance (G), and DOS as a function of energy E.
Furthermore, for these z4ANR-based samples, which now have a double-barrier structure (DBS) potential profile imposed on them (i.e., central channel length: l = 2 (barrier) + 10 (well) + 2 (barrier) = 14 atomic layers with barrier height 0.7 eV), Figs. 6(a)-6(d) display the spectral diagrams of T , G, and DOS as a function of E.
Assuming a perfect flat-band potential, we have T(E) = 1, and thus ( ) ( ) T E M E = , for the z4ANR-based samples. According to calculations using Eqs. (23) and (24), these flat-band samples have a staircase-like conductance, which is the number of propagating channels M(E) multiplied by the conductance quantum 2e 2 /h, as displayed in the right panel of Fig. 4 and Figs. 5(a)-5(d). (35)(36)(37) As presented in Figs. 5(a)-5(d) and Figs. 6(a)-6(d), it can be seen that in the DBS imposition, the quantization staircases of conductance spectra are destroyed and show more complicated oscillation behaviors due to the scattering from the DBSs. Namely, it is shown that the transition of conductance is from the quantized conductance in flat-band structures to resonant-tunneling conductance in DBSs.   A van Hove singularity can be seen as a singularity (nonsmooth point) in the DOS spectra of crystalline solids. The van Hove singularities occur where the derivative of the DOS with respect to energy E diverges, which yields a strikingly sharp peak in the DOS spectra. (38)(39)(40)(41)(42)(43) Therefore, in a perfect flat-band z4ANR, the DOS diverges at the onset of each subband, as shown in Figs. 5(a)-5(d) and Figs. 3(a)-3(d). Figures 5(a)-5(d) display the strikingly sharp and asymmetric peaks in the DOS spectra, which are also called van Hove singularities.
In quantum transport devices such as DBSs, a series of discrete and fine confined levels (or states) exist in the central channel, and a continuous distribution of states exists in the left and right contacts. When the central channel and the two contacts are joined, the discrete channel levels couple and then acquire a few states of the two contacts. (12,32) Therefore, the effect of coupling is to broaden the DOS in the central channel from its original discrete lines to a continuous spectrum, (12,32) as shown in Figs. 6(a)-6(d). With increasing DBS barrier height, the coupling strength of the two contacts to the central channel decreases, and thus the line-shape sharpness and depth of the DOS peaks and dips increase. Therefore, the DBS has been verified to be a suitable structure for exploring the coupling strength of contacts to a central channel, and an increase in DBS barrier height yields the DOS modification of the central channel. As shown in Figs. 6(a)-6(d), when energy E coincides with any resonant-tunneling confined levels of the DBS, irregularly spaced resonant-tunneling peaks occur in DOS spectra.

Conclusions
In this paper, we have developed an atomistic full-band and full-quantum model for calculating the electronic transport characteristics of z4ANR-based devices, which also involves the SOC effect. We have derived in detail the theoretical expressions for z4ANR-based devices such as their wave function, Hamiltonian matrix, and the Schrödinger-like equation in the NEGF form, which is based on the complex energy-band method. The proposed method is straightforward, nonrecursive, and thus computationally cost-efficient. Using the developed method, we have calculated and obtained important findings on z4ANR-based devices, such as their conductance quantization, van Hove singularities in the DOS, and the effect of contact interactions on the channel.     1  2  cos  ,  2  1  2  sin  ,  2  1  2  cos  ,  2  1  cos  ,  2  1 1 k is a wave vector (k ⊥ , || k ), θ is the angle between the bond and the z direction, a ( 3 b = ) is the lattice constant, b is the bond length projected on the layer plane, and a total of five interaction parameters exist, namely, Δ, V ppσ , V ppπ , V spσ , and V ssσ . Δ is related to the on-site energy difference between the s and p orbitals, while the remaining four parameters represent the nearest-neighbor interactions.
where σ x , σ y , and σ z are the Pauli spin matrices and ξ 0 is the strength of the SOC effect. The parameters in this appendix are shown in Tables I and II

Appendix B: Hamiltonian of Zigzag-edge Group-IVA Nanoribbon
The z4ANR Hamiltonian, which is calculated using the tight-binding model with the SOC effect and only nearest-neighbor interactions, in the | , , k j τ β α ⊥ > basis can be written in the N × N block matrix form as ( j = 1, 2, ..., N)   The g 2 matrix can be expressed in the 16 × 16 matrix form, which has the same form as the h 2 matrix, i.e., g 2 = h 2 .
The z ± matrix can be written in the 16 × 16 matrix form as