Chebyshev Polynomial Method to Landauer-B\"uttiker Formula of Quantum Transport in Nanostructures

Landauer-B\"uttiker formula describes the electronic quantum transports in nanostructures and molecules. It will be numerically demanding for simulations of complex or large size systems due to, for example, matrix inversion calculations. Recently, Chebyshev polynomial method has attracted intense interests in numerical simulations of quantum systems due to the high efficiency in parallelization, because the only matrix operation it involves is just the product of sparse matrices and vectors. Many progresses have been made on the Chebyshev polynomial representations of physical quantities for isolated or bulk quantum structures. Here we present the Chebyshev polynomial method to the typical electronic scattering problem, the Landauer-B\"uttiker formula for the conductance of quantum transports in nanostructures. We first describe the full algorithm based on the standard bath kernel polynomial method (KPM). Then, we present two simple butefficient improvements. One of them has a time consumption remarkably less than the direct matrix calculation without KPM. Some typical examples are also presented to illustrate the numerical effectiveness.


I. INTRODUCTION
Landauer-Büttiker formula plays an important role in the study of electronic quantum transports in nanostructures [1][2][3][4] , molecular systems 5 , and even DNAs 6 . It also plays an important role in calculating the thermal 7-9 , optical 10,11 and phonon 12 transports in quantum structures. Landauer-Büttiker formula relates the electronic conductance of a two-terminal or multi-terminal device to the quantum transmission [1][2][3] . The quantum transmission can be expressed in terms of Green's functions, which is a standard numerical tool today 4,[13][14][15] . Since this is a real space method, it is computationally demanding for a system related with large number of orbital basis, e.g., large size systems 16,17 , biological molecules 6 and (quasi-)incommensurate systems 18-20 . In the recent decade, a powerful numerical method treating with Hamiltonians on large Hilbert spaces has attracted attention, the kernel polynomial method (KPM), such as the Chebyshev expansion 21 . In most KPM calculations, the only matrix operations involved are product between sparse matrices (Hamiltonians) and vectors, and matrix traces. For a a) Electronic mail: yanyang@gzhu.edu.cn sparse matrix with dimension D, the matrix vector multiplication is only an order O(D) process. Thus the calculation of N moments of Chebyshev terms needs O(ND) operations and time 21 . A direct application of KPM is the calculation of the spectral function of an isolated system [21][22][23] . Taking advantage of appropriate analytical continuation, one can arrive at the evaluation of Green's functions 21,24 . Expressions of physical quantities in terms of KPM have been developed recently [25][26][27][28][29][30][31][32][33] , including the applications to superconductors 34 , topological materials 35,36 , quantum impurity problems [37][38][39] and ab initio calculations 40,41 . However, these methods are applicable to bulk or isolated systems 21,33 , not to scattering processes between leads in open systems, which corresponds to a realistic experimental setup 4 .
In this paper, we will propose some KPM methods to calculate the Landauer-Büttiker transmission in a two-terminal system: a conductor connected to left (L) and right (R) leads, as illustrated in Fig. 1. The transmission through the conductor can be written in terms of Green's functions, where the leads manifest themselves as self energies 4 . This is typical context of an open system coupled to a bath 42,43 . We first fully describe this problem as a generalization of the standard bath technique of KPM 42 , where the needed self energies and dressed Green's function as Chebyshev polynomials of some sparse matrices. To reduce the numerical consumption, we then propose two practical improvements. One of them can largely simplify the self-consistent calculation of the self energies, and the second of them can even avoid this selfconsistent process, which has a much less time and space consumption than those of the direct method of matrix evaluation without KPM. This paper is organized as follows. After the general introduction, in Sections II and III, we briefly introduce the basic knowledge of Chebyshev polynomials and Landauer-Büttiker formula, respectively. In Section IV, we describe the algorithm of calculating Landauer-Büttiker formula with Chebyshev polynomials, including the calculation of dressed Green's function and lead self energies, following the standard bath method of KPM. To reduce the numerical demanding, we propose two practical improvements in Section V. Some numerical examples are presented in Section VI. In Section VII, we provide a summary and some outlooks for future works.

II. CHEBYSHEV EXPANSION AND THE KERNEL POLYNOMIAL METHOD
In this section, we briefly summarize the definition and basic properties of Chebyshev polynomials that will be used. The Chebyshev polynomials T n (x) with x ∈ [−1, 1] are in the explicit form as 21 which satisfy the recursion relations, The scalar product is defined as with the weight function (π √ 1 − x 2 ) −1 . It is thus easy to verify the orthogonality relation between Chebyshev polynomials, In terms of these orthogonality relations (4), a piecewise smooth and continuous function f (x) with x ∈ [−1, 1] can be expanded as with expansion coefficients µ n Practically, the function f (x) should be numerically reconstructed from a truncated series with the first N terms in Eq. (5). However, experiences show that the numerical performance of this simple truncation is bad, with slow convergence and remarkable fluctuations (Gibbs oscillations) 21 . This can be improved by a modification of the expansion coefficients as µ n → g n µ n , where {g n } is the kernel. In other words, appropriate choices of the kernel {g n } will make the truncated series a numerically better approximation of the function 21 Among different kernels, here we adopt the Jackson kernel with the explicit expression as 21 which is suitable for the applications related to Green's functions.
Besides a numeric function f (x), the Chebyshev expansion can also be used to approximate the function of a Hermitian operator H (or equivalently its matrix H in an appropriate representation), if the eigenvalue spectrum of H is within the interval [−1, 1] 21 . For a general Hermitian operator, e.g. a Hamiltonian H with maximum (minimum) eigenvalue E max (E min ), this condition of spectrum can be satisfied by simply performing an appropriate rescaling on the matrix (and also on the energy scale), so that the the spectrum ofH is within [−1, 1]. Here the parameter ζ > 0 is a small cutoff to avoid numerical instabilities at the boundaries ±1. A proper rescaling, i.e., an appropriately small ζ will reduce the necessary N for a certain expansion to reach the same precision. Throughout this work, we fix ζ = 0.01. In practical uses, the lower and upper bounds of H can be estimated by using sparse matrix eigenvalue solvers, e.g., the FEAST algorithm of Intel MKL. After the calculation of physical properties with the help of Chebyshev polynomials, their correct dependence on the energy E can be restored by a simple inverse transformation of Eq. (9). Therefore in the following, we will always consider that the operator matrices have been rescaled according to Eq. (9) before they enter Chebyshev polynomials, and the tilde hats on the operators and eigenvalues will be omitted. It can be shown that 21,24,35 , the retarded (advanced) Green's function at energy E can be expanded in terms of Chebyshev kernel polynomials as with coefficient matrices Now the broadening η does not explicitly appear in the matrix elements. Rather, it is associated with N, the number of expansion moments. Larger N corresponds to a smaller η. Notice G and µ n are also operators. In a certain representation, these operators can be explicitly written as corresponding matrices H, G and µ n with the same size. Throughout this manuscript, all matrices will be written in a bold form of the corresponding operator. In this section, we briefly review the Landauer-Büttiker formula represented as Green's functions. Consider the twoterminal transport device illustrated in Fig. 1, with one conductor connected to two semi-infinite leads. Formally, the Hamiltonian of this combined system can be written as where H C is the Hamiltonian of the conductor, H L (H R ) is that of the left (right) lead, and H CL (H CR ) is the coupling from the conductor to the left (right) lead. It is convenient to write these Hamiltonians in the real space representation (tight binding model) as matrices. For example, the real space Hamiltonian of a conductor (lead) can be expressed in a generic second quantization form as with c β the annihilation operator of the spinorbital β in the conductor (lead where Σ L (Σ R ) is the self energy of the left (right) lead. The technique of self energy liberate one from inserting the full Hamiltonian (14) into Eq. (11) to obtain the dressed Green's function of the conductor. The self energy is the result of integrating out the degree of freedom of the lead 4,45 , i.e., where G r p is the Green's function of lead p, and H pC is the coupling Hamiltonian between lead p and the conductor. In the real space representation, G r p is an infinite dimensional matrix because the lead is semi-infinite. However, since only a few spinoribitals of the lead is connected to the conductor through H pC , in the evaluation of Eq. (17), we only need to know the "surface" subset of the matrix G r p , i.e., those matrix elements G r p αβ with α and β running over spinorbitals connected to the conductor. This subset will be called the surface Green's function.
At zero temperature, the two-terminal conductance G in Fig.1 is represented as the Landauer-Büttiker formula 1-4 , where e is the elementary charge, h is the Planck constant, and T is the transmission through the conductor. This transmission T at Fermi energy E can be expressed in terms of Green's functions as 4,13 , where Traditionally, the self energies (17) of the leads can be calculated explicitly by a direct diagonalization method 44 or an iterative method 46 . Afterwards, they are inserted into Eqs. (16), (20) and finally (19) for the evaluation of the transmission. In this process, the most time-consuming step will be the calculation of lead self energies, and the matrix inversion (which does not preserve the sparseness of the matrix) in Eq. (16). For a two-terminal device simulation where the conductor lattice can be well divided into layers of sites (layers should be defined in such a way that hoppings only exist between nearest layers), the simulation can be decomposed into a layer-to-layer recursive method, which is based on the Dyson equation for Green's functions 13,47 . This decomposition can remarkably reduce the time and space consumption in calculations. However, this recursive method will be technically tedious for a multi-terminal setup, and even impossible for, say, a twisted bilayer graphene 48,49 . In these examples, one still needs to calculate the full-size and dense matrices associated with Hamiltonians and Green's functions directly. In the following, we will investigate algorithms based on KPM to calculate Eq. (19), with slightly different steps.

IV. STANDARD BATH CHEBYSHEV POLYNOMIAL METHOD
Before evaluating the transmission function (19) from the Hamiltonian (14), two steps are essential: First, solving the self energies [Eq. (17)]; Second, inclusion of them into the conductor's Green's function as Eq. (16). The numerical treatments of these steps by direct matrix calculations have been very mature and well-known 4,13 . However, in the context of KPM, the realization of these steps is not easy nor straightforward, especially if one insists to avoid calculations related to large dense matrices. We achieve this goal by a generalization of the bath technique of KPM 42 , which will be described here in detail. In Section V B, another distinct algorithm will be introduced.
The lead connected to the conductor is semi-infinitely long and therefore it can be viewed as a bath 42 . The central task of obtaining the lead self energy [Eq. (17)] is to calculate the surface Green's function G r p αβ of lead p, with α and β running over the surface which will be connected to the conductor. In the context of KPM method, we need to calculate the Chebyshev coefficient matrix µ αβ n in Eq. (12) of the lead. This, of course, cannot be calculated by using Eq. (13) directly, as the lead Hamiltonian matrix H p is infinite dimensional. Instead, we will use a self consistent method as described below.

A. Basic Definitions
First, some useful mathematical structures related to an isolated lead p will be constructed. As suggested in Ref. 42 , we define the Chebyshev vectors as with |vac describing the lead vacuum, i.e., f † α |vac = 0, and f † α the creation operator in the lead at spinorbital state α. These Chebyshev vectors are not orthonormal and the scalar product By comparing with Eq. (13), one can see that, this matrix µ n is just the n-th Chebyshev coefficient matrix of the lead's Green's function.
In other words, in the subspace H α , the effect of H p can be expressed in a matrix form as with Notice that ( H α p ) mn = m α |H p |n α owing to the nonorthogonality of these Chebyshev vectors. For a truncation with N Chebyshev terms (7), the size of the matrix H α p is N × N.
Following Eq. (21), another useful relation can be obtained as

B. Dressed Green's Function
Suppose the Chebyshev coefficient matrices µ n of lead p have been known. Now we connect a conductor C to the lead. The Hamiltonian H C of the conductor is in the form of Eq. (15), and the size of the corresponding matrix H C is M × M. Without loss of generality, we consider the conductorlead coupling to be the following simple form where W is the effective "width" of the cross section, c α ( f α ) is the annihilation operator in the conductor (lead), and t α denotes the hopping matrix elements. As in most practical cases, we have considered these W hopping bonds are coupling sites between the conductor and the lead in a one-to-one way. Based on above definitions, now we can approximately express the full Hamiltonian H C + H p + H Cp + H.c. in the finitedimensional representation with basis ordered as |1 , · · · , |β , · · · |M , |0 1 , · · · , |n 1 , · · · , |N − 1 1 , · · · , |n α , · · · , |0 W , · · · , |n W , · · · , |N − 1 W , where |β (1 ≤ β ≤ M) are spinorbital basis states in the conductor, and |n α (0 ≤ n ≤ N − 1 and 1 ≤ α ≤ W ) are Chebyshev vectors of lead p defined in Eq. (21). It can be easily shown that, the full Hamiltonian in this representation is a (M + W × N)-dimensional sparse matrix with nonzero blocks illustrated as follows: Here, H C is the Hamiltonian matrix of the isolated conductor with size M × M, and H α p are N × N matrices defined in Eq. (25) associated with lead p. As for the conductor-lead coupling sub-matrices, each L α p is a W × N matrix in the form as where t α are the coupling hopping integral defined in Eq (28) This matrix R is determined if H C , H Cp and µ β α n are known, and it plays the central role in the Chebyshev approach to quantum systems coupled to a bath 42 . It has been shown that 42 , the dressed Green's function (E + iη)I − H C − Σ p −1 of the conductor coupled to lead p can be obtained by using Eq. (12) and Eq. (13), with H replaced by R 42 .

C. Self Energy
In most applications, the coefficients µ β α n associated with the lead surface is not known a priori, and practically they cannot be obtained by using Eq. (22). However, the above bath method also offers a practical algorithm of calculating the lead coefficients µ β α n in a self-consistent way, which will be described as follows. The lead is a semi-infinite crystal whose one-dimensional unit cell can be defined in a natural way. In Fig. 2, we illustrate a lead extending infinitely to the right direction, with the unit cell marked by the green dashed frame. The lead Chebyshev coefficients of the left surface are µ β α n , with α and β only running over the left surface spinorbitals which will be connected to the conductor. Now we couple J should hold. Practically, we can start from a guess of the lead Chebyshev coefficients, e.g., µ β α n = 0, and then repeatedly couple J unit cells to the left surface of the lead and calculate the Chebyshev coefficients ν β α n associated with the new surface, until the self-consistent condition (33) are satisfied within a given error. Although the choice µ β α n = 0 fail to meet the rule 0 α |0 α = vac| f α f † α |vac = µ αα 0 = 1, it does not affect the final convergence. A larger number J of unit cells will consume more time for each iteration step, but will reduce the number of iteration steps. Therefore an appropriate J should be carefully chosen for a concrete model. After the lead Chebyshev coefficients µ β α n are known, the surface Green's function can be obtained through Eq. (12), and the self energy through Eq. (17).

D. Counting Both Leads in
So far, we have shown that, replacing H in Eq. (13) with R defined in Eq. (30) will give rise to the dressed Green's function of the conductor when coupled to a single lead p, For a two-terminal device, the inclusion of left (L) and right (R) leads can be similarly achieved by a trivial generalization of the matrix in Eq. (30) as The central task of KPM is to obtain the corresponding Chebysheve coefficient matrices. One merit of the KPM is that these Chebyshev coefficients are independent of the energy E, i.e., the transport properties over the full energy spectrum are known if the corresponding Chebyshev coefficients have been calculated out. Particularly, when plotting Fig. 4 , one only needs to calculate the lead Chebyshev coefficients µ β α n for once, and then the energy dependence enters simply through Eq. (12) which is numerically cheap. In the traditional matrix inversion method, on the other hand, the full process of calculating the self energy (17) and the dressed Green's function (20) should be performed separately for different energies, which are numerically independent.
The algorithm described in Section IV was based on a mathematically rigorous realization of the standard bath ap-proach of the KPM 42 , which is referred as the "standard bath KPM". In the practical simulations, however, calculating Chebysheve coefficient matrices from this standard method might be very numerically demanding on central Q4 processing unit (CPU) based computers. For example, in the calculation of conductance in terms of KPM, the most time consuming process is the self-consistent calculation of Chebyshev coefficients µ β α n of the leads. As a matter of fact, the requirement of including all details of leads into the calculation like this is a notoriously expensive cost in many quantum transport simulations [50][51][52][53] . Now we will present two practical improvements of the algorithm.

A. Chain Shaped Leads
The first convenient simplification is to reduce the shape of leads into independent and semi-infinite 1D chains, as shown in Fig. 3 (a). Without transverse coupling in the lead, the coefficients are diagonal µ β α n = δ β α µ αα n and they are identical for different α. Now we only need to self-consistently calculate the Chebyshev coefficients of a 1D chain with width W = 1, and the dimension of the matrix R in Eq. (30) is reduced from W × J + W × N to M + N = J + N. However, the mismatch between the leads and the conductor will give rise to additional scattering at their boundaries. Therefore, this brute simplification is mostly suitable for topological materials where backscattering have been prohibited by protections from topology and/or symmetry 56 . See Examples B and C in the Section VI for simulation results.

B. Finite Lead Approximation
Here we propose another simple but efficient approximation to circumvent these difficulties. The original setup was that both leads should be semi-infinitely long, as illustrated in Fig. 1. Now we approximate both leads by two finite ones, as presented in Fig. 3 (b). It is reasonable to imagine that the result will approach the correct one when their lengths N L x and N R x are sufficiently large. Now the conductor and leads are perfectly matched, so there will be no scattering on their boundaries.
If the lead lengths N L x and N R x needed to arrive within some precision are numerically acceptable, this algorithm will be numerically more superior than the standard bath KPM described above. For instance, the dressed Green's function can be obtained from Eq. (12) directly, only if H is the coupled Hamiltonian matrix of the whole system, the conductor and two finite leads. Now, the sub-matrix G r(a) i j (E, H) (with i, j running over the conductor sites) is naturally the approximation of the dressed Green's function (20) 4 . In fact, due to the simple algebraic structure of Eq. (19), one only needs the matrix indices (i, j) running over the boundary sites connected to two leads. This process avoids the construction and calculation of complicated and non-Hermitian matrices like R defined in Eq. (30). A non-Hermitian matrix has complex eigenval- ues, leading to difficulties of scaling itself with Eq. (9) by its maximum and minimum eigenvalues, but an appropriate scaling of the matrix is key in the context of KPM.
Similarly, the surface Green's function of the lead can also be approximated by that of the finite one from Eq. (12), then the self energy is calculated with the help of Eq. (17). In brief, this method circumvents all complicated steps of the standard bath KPM described in Section IV, especially the self-consistency calculation of the self energy, which is very time-consuming.

VI. EXAMPLES
In this section, we present results from above KPM to calculate the two-terminal conductance of some example models.

A. Square Lattice Conductor with Square Lattice Leads
The first example is the two-dimensional square lattice with nearest hopping t, where α and β are indices for sites (each with only one spinorbital) in the conductor and leads, and α, β run over all nearest site pairs. The size of the conductor is L × W , and the widths of both leads are also W . The results from the standard bath KPM are shown in Fig.  4. Here the transmission T a function of energy E is plotted as the solid lines, for conductor sizes: (a) 25 × 25, (b) 60 × 60 (b) and (c) 100 × 100. Different line colors corresponding to different Chebyshev terms N are also shown in each panel. As a comparison, the red dashed line is the result from a direct matrix calculation of Eq. (19), without using KPM. Without any disorder in the conductor, the conductance is quantized as plateaus with values p e 2 h , with p the number of active channels at energy E 4 . Smaller N effectively corresponds to a stronger dephasing 21,36 , and therefore the conductance is not perfectly quantized. The largest deviation at smaller N (green lines in Fig. 4) happens around the band center E = 0, which is a van Hove singularity. This enhanced scattering is caused by the extremely large density of states around such a singularity 54 .
On the other hand, larger conductor sizes need more Chebyshev terms to reach the perfect conductance value of quantum transport. This is understandable since a larger conductor gives rise to a longer journey for the electron to experience the dephasing, which is induced by the finiteness of Chebyshev terms N. Due to the tedious process of the standard bath KPM, especially the self-consistent calculation of the self energies, it is even more time-consuming than the direct matrix calculation.

B. Square Lattice Conductor with Chain Shaped Leads
The results are shown as the black solid line of Fig. 5  (b), where the result of square lattice leads (red dashed line) are also presented as a comparison. Similar to the popular method of wide band approximation 50,51,55 , this simplification largely reduces the time and space consumptions of the selfconsistent calculation of the lead self energy. On the other hand, the mismatch between the lead and the conductor will give rise to remarkable additional scattering on the interface, leading to a distinct reduction of the conductance compared with the case of perfect leads.

C. Topological Material Conductor with Chain Shaped Leads
However, such scattering will be practically avoided if the conductor is a topological material with robust transport FIG. 5. The black line is the result for Example B, square lattice conductor with chain shaped leads as illustrated in Fig. 3 (a). The conductor size is 100 × 100, and the number of Chebyshev polynomial terms is N = 10000. The red line is identical to that in Fig. 4 (c), the result of square lattice leads (from direct matrix calculation), as a reference.
FIG. 6. Black solid lines are results for Example C, topological material conductor with chain shaped leads. The conductor is the quantum anomalous Hall model defined in Eqs. (36) and (37), and the chain shaped leads are illustrated in Fig. 3 (a). The red dashed lines mark the quantized conductance value in units of against backscattering, or three dimensional conductor with sufficiently large number of transport channels. Therefore this method is most applicable in these contexts. Here we adopt the typical model of the quantum anomalous Hall effect, the spin-up component of the Bernevig-Hughes-Zhang (BHZ) model 56 , defined on a two-orbital square lattice. The Hamiltonian in the k space can be written as 56 where h αβ (k) is a 2 × 2 matrix defined as with σ x,y,z the Pauli matrices acting on the space of two orbitals. The Chern number of this model is 1 when B/M > 0, so that there will be a pair of topological edge states in the bulk gap (− M 2 , M 2 ). Due to the topological origin, the edge states will contribute a quantized conductance 1 × e 2 h that is robust against elastic backscattering. Fig. 6 is the numerical results (black lines) of this topological model from our Chebyshev approach, with chain shaped leads as illustrated in Fig. 3 (a). Panels (a) and (b) are for conductor sizes 60 × 60 and 80 × 80, respectively, and the red lines mark the reference position of the quantized conductance. We can see that in most of bulk gap region, the simulated conductance is perfectly consistent with that predicted by the topological invariant theory. For example in Fig. 6 (a), the numerical match can be larger than 99.5% near the gap center E = 0. As in previous examples, larger conductor sizes needs more Chebyshev terms to reach the perfect quantum transport value. Moreover, the transport near gap edges are more sensitive to scattering, which is a natural consequence of bulk-edge mixing 57 .

D. Square Lattice Conductor with Finite Lead Approximation
In Fig. 7, we present the results from the finite lead approximation, as introduced in Section V B. Typically, the necessary length of the finite lead is less than 50 times of the conductor length, N L x , N R x 50 × L, to achieve a relative error less than 2% throughout most of the energy spectrum. This necessitates a sparse matrix with dimension 100 × L × W to store the Hamiltonian of the conductor and leads. This matrix is structurally simpler, and usually not larger than R [Eq. (34)] in the standard bath KPM, with an N dependent dimension (2N + L) ×W (N is the number of Chebyshev terms, typically ∼ 10 3 − 10 4 ) . Moreover, the calculation of self energies is also much simpler and faster than in the standard bath KPM, since it is now a one-time process and no self-consistent loops are needed here.
Moreover, here we will show that this method is also much faster than the direct matrix evaluation of Eq. (19) without KMP. In Fig. 8, the computation time of these two methods are plotted as functions of the length of the square shaped conductor. Here the "computation time" means the full time of obtaining one curve like those in Fig. 7, by calculating transmissions over 800 energy points. It can be seen that the KPM with finite leads is numerically much more advantageous than the traditional direct matrix calculation. As has been discussed above, since the Chebyshev coefficients have contained information on the full energy spectrum, this improved KPM method will be even more superior when one needs data on more energy points. Furthermore, in the context of simulating irregular shaped conductors with an irregular structure of Hamiltonian matrix, since the calculation cannot be reduced to a layer-to-layer recursive one 13,47 , this method we propose will be very applicable.

VII. SUMMARY AND OUTLOOK
In summary, we introduced the Chebyshev polynomial method to the Landauer-Büttiker formula of the two-terminal transport device, by a generalization of the standard bath technique of KPM. In this formula, the dressed Green's function can be expressed as Chebyshev polynomials of the matrix R defined in Eq. (30) or Eq. (34), and the self energies can be calculated through the dressed Green's function in a selfconsistent way. During this process, the most resource consuming step is the calculation of self energies of the leads. A simple solution is to reduce the topology of the leads to parallel and decoupled atomic chains, but the price is additional scattering on the interfaces. Another solution is to approximate the leads as finite ones with sufficient length. This algorithm avoids complicated matrix calculations in the standard bath KPM (especially the self-consistent process of obtaining self energies), and also avoids notable boundary scattering in the chain shaped lead method. The numerical experiments verified that this method has a much less numerical cost than that of the traditional method of direct matrix calculation without KPM.
Since the leads themselves are not the object of study and there is a wide freedom of choosing leads. One of the future efforts is to find an appropriate design of leads or leadconductor coupling [51][52][53]58 , with small resource demanding in the Chebyshev polynomial representation, while with small backscattering on the interfaces to the conductor. Furthermore, our method can also be generalized to Cheyshev forms to linear responses of other degree of freedoms of quantum transports structures 7,10,12,59 .
search Fund from Guangzhou University under Grant No. RQ2020082.

VIII. DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.