Chebyshev polynomial method to Landauer–Büttiker formula of quantum transport in nanostructures

The Landauer–Büttiker formula plays an important role in the study of electronic quantum transport in nanostructures,^1–4 molecular systems,⁵ and even DNAs.⁶ It also plays an important role in calculating the thermal,^7–9 optical,^10,11 and phonon¹² transport in quantum structures. The Landauer–Büttiker formula relates the electronic conductance of a two-terminal or multi-terminal device to the quantum transmission.^1–3 The quantum transmission can be expressed in terms of Green’s functions, which is a standard numerical tool today.^4,13–15 Since this is a real space method, it is computationally demanding for a system related to a large number of orbital bases, e.g., large size systems,^16,17 biological molecules,⁶ and (quasi-)incommensurate systems.^18–20 

In the recent decade, the kernel polynomial method (KPM), a powerful numerical method working with Hamiltonians on large Hilbert spaces has attracted attention, such as the Chebyshev expansion.²¹ In most KPM calculations, the only matrix operations involved are the product between sparse matrices (Hamiltonians) and vectors, and matrix traces. For a sparse matrix with the 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.²¹ A direct application of the KPM is the calculation of the spectral function of an isolated system.^21–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 the KPM have been developed recently,^25–33 including the applications to superconductors,³⁴ topological materials,^35,36 quantum impurity problems,^37–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.⁴ 

In this paper, we will propose some KPMs 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.⁴ This is a typical context of an open system coupled to a bath.^42,43 We first fully describe this problem as generalization of the standard bath technique of KPM,⁴² where the needed self-energies and dressed Green’s function are 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 self-consistent process, which has a much less time and space consumption than that of the direct method of matrix evaluation without KPM.

This paper is organized as follows. After the general introduction, in Secs. II and III, we briefly introduce the basic knowledge of Chebyshev polynomials and Landauer–Büttiker formula, respectively. In Sec. IV, we describe the algorithm of calculating the 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 Sec. V. Some numerical examples are presented in Sec. VI. In Sec. VII, we provide a summary and some outlooks for future works.

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 given as²¹ 

which satisfy the recursion relations

The scalar product is defined as

with the weight function $(π1−x2)−1$⁠. It is thus easy to verify the orthogonality relation between Chebyshev polynomials,

In terms of this orthogonality relation (4), a piecewise smooth and continuous function f(x) with x ∈ [−1, 1] can be expanded as

with expansion coefficients

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).²¹ 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²¹ 

Among different kernels, here we adopt the Jackson kernel with the explicit expression²¹ 

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].²¹ For a general Hermitian operator, e.g., a Hamiltonian H with the 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),

with

so that the spectrum of $H̃$ 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. A larger N corresponds to a smaller η. Note that 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 the bold form of the corresponding operator.

In this section, we briefly review the Landauer–Büttiker formula represented as Green’s functions. Consider the two-terminal 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_β being the annihilation operator of the spin orbital β in the conductor (lead). Here, H is a matrix with elements H_αβ.

Due to the coupling to leads, now the (retarded) Green’s function of the conductor $GCr$ is, of course, not the original bare one $E+iη−HC−1$⁠. Thanks to the Dyson equation of Green’s functions, it can be expressed as the dressed one as^4,44

where Σ_L (Σ_R) is the self-energy of the left (right) lead. The technique of self-energy liberates 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 $Gpr$ is Green’s function of the lead p and H_pC is the coupling Hamiltonian between the lead p and the conductor. In the real space representation, $Gpr$ is an infinite dimensional matrix because the lead is semi-infinite. However, since only a few spin orbitals of the lead are connected to the conductor through H_pC, in the evaluation of Eq. (17), we only need to know the “surface” subset of the matrix $Gpr$⁠, i.e., those matrix elements $Gprαβ$ with α and β running over spin orbitals 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⁴⁴ or an iterative method.⁴⁶ Afterward, they are inserted into Eqs. (16) and (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 molecule.^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 the KPM to calculate Eq. (19), with slightly different steps.

Before evaluating the transmission function (19) from the Hamiltonian (14), two steps are essential: first, solving the self-energies [Eq. (17)] and, second, inclusion of them into the conductor’s Green’s function as Eq. (16). The numerical treatments of these steps by direct matrix calculations are very mature and well-known.^4,13 However, in the context of the 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 generalization of the bath technique of KPM,⁴² which will be described here in detail. In Sec. 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.⁴² The central task of obtaining the lead self-energy [Eq. (17)] is to calculate the surface Green’s function $Gprαβ$ of the lead p, with α and β running over the surface that will be connected to the conductor. In the context of the KPM, 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.

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α†$ being the creation operator in the lead at spin orbital 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. The series of Chebyshev vectors defined in Eq. (21) spans a Hilbert space $Hα$⁠. As can be seen from the definition, $Hα$ is a subspace of the Fock space for the lead operator $fα(†)$⁠. From the recursion relation, Eq. (2), it is easy to conclude the operation of H_p on $Hα$ as

In other words, in the subspace $Hα$⁠, the effect of H_p can be expressed in the matrix form as

with

Note that $(H^pα)mn≠⟨mα|Hp|nα⟩$ owing to the non-orthogonality of these Chebyshev vectors. For 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

Suppose the Chebyshev coefficient matrices μ_n of the lead p are 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 conductor–lead coupling to be in the following simple form:

where W is the effective “width” of the cross section, c_α (f_α) is the annihilation operator of the conductor (lead), and t_α denotes the hopping matrix elements. As in most practical cases, we have considered these W hopping bonds as 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 finite-dimensional representation with the basis ordered as

where |β⟩ (1 ≤ β ≤ M) are spin orbital basis states in the conductor and |n_α⟩ (0 ≤ n ≤ N − 1 and 1 ≤ α ≤ W) are Chebyshev vectors of the lead p defined in Eq. (21). It can be easily shown that the full Hamiltonian in this representation is an (M + W × N)-dimensional sparse matrix with nonzero blocks illustrated as follows:

Here, H_C is the Hamiltonian matrix of the isolated conductor with the size M × M and $H^pα$ are N × N matrices defined in Eq. (25) to be associated with the lead p. As for the conductor–lead coupling sub-matrices, each $Lpα$ is a W × N matrix in the form

where t_α are the coupling hopping integral defined in Eq. (28) and $μnβα$ of the lead are defined in Eq. (22), with {α, β} only running over the surface spin orbital states connected to the conductor. On the other hand, each $Mpα$ is an N × W matrix with only one nonzero element,

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.⁴² It has been shown that⁴² the dressed Green’s function $(E+iη)I−HC−Σp−1$ of the conductor coupled to the lead p can be obtained by using Eqs. (12) and (13), with H replaced by R.⁴² 

In most applications, the coefficients $μnβα$ associated with the lead surface are 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 spin orbitals that will be connected to the conductor. Now, we couple J unit cells to the left surface of this lead. Then, the Chebyshev coefficients $νnβα$ of the new left surface can be calculated by using the process introduced above. On the other hand, because these unit cells are just a natural extension of the lead, this new composite system (J unit cells coupled to the original lead) is also semi-infinite and is essentially equivalent to the original lead. As a result, the self-consistent condition

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) is satisfied within a given error. Although the choice $μnβα=0$ fails 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).

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, $GC=EI−HC−Σp−1$⁠. For a two-terminal device, the inclusion of left (L) and right (R) leads

can be similarly achieved by trivial generalization of the matrix in Eq. (30) as

once the Chebyshev coefficients μ_n associated with both leads are known.

The central task of the KPM is to obtain the corresponding Chebyshev 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. 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 Sec. IV was based on a mathematically rigorous realization of the standard bath approach of KPM,⁴² which is referred to as the “standard bath KPM.” In the practical simulations, however, calculating Chebyshev coefficient matrices from this standard method might be numerically demanding on central processing unit (CPU) based computers. For example, in the calculation of conductance in terms of the 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 notoriously expensive in many quantum transport simulations.^50–53 Now, we will present two practical improvements of the algorithm.

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 has been prohibited by protections from topology and/or symmetry.⁵⁶ See examples B and C in Sec. VI B for simulation results.

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 $NxL$ and $NxR$ are sufficiently large. Now, the conductor and leads are perfectly matched, so there will be no scattering at their boundaries.

If the lead lengths $NxL$ and $NxR$ needed to arrive within some precision are numerically acceptable, this algorithm will be numerically superior to 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 $Gijr(a)(E,H)$ (with i, j running over the conductor sites) is naturally the approximation of the dressed Green’s function (20).⁴ 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 such as R defined in Eq. (30). A non-Hermitian matrix has complex eigenvalues, 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 the KPM.

Similarly, the surface Green’s function of the lead can also be approximated by that of the finite one from Eq. (12), and 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 Sec. IV, especially the self-consistency calculation of the self-energy, which is very time-consuming.

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

The first example is the two-dimensional square lattice with the nearest hopping t,

where α and β are indices for sites (each with only one spin orbital) 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 as a function of the energy E is plotted as the solid lines, for different conductor sizes: (a) 25 × 25, (b) 60 × 60, and (c) 100 × 100. Different line colors corresponding to different Chebyshev terms N are also shown in each panel. For 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 $pe2h$⁠, with p being the number of active channels at energy E.⁴ A 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.⁵⁴ 

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.

The results are shown as the black solid line in Fig. 5, where the result of square lattice leads (red dashed line) is also presented for comparison. Similar to the popular method of wide band approximation,^50,51,55 this simplification largely reduces the time and space consumption of the self-consistent 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.

However, such scattering will be practically avoided if the conductor is a topological material with robust transport against backscattering, or a three-dimensional conductor with a 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,⁵⁶ defined on a two-orbital square lattice. The Hamiltonian in the k space can be written as⁵⁶ 

where h_αβ(k) is a 2 × 2 matrix defined as

with σ_x,y,z being 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 $(−M2,M2)$⁠. Due to the topological origin, the edge states will contribute a quantized conductance $1×e2h$ that is robust against elastic backscattering.

Figure 6 shows the numerical results (blue 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 the 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 need more Chebyshev terms to reach the perfect quantum transport value. Moreover, the transport near gap edges is more sensitive to scattering, which is a natural consequence of bulk–edge mixing.⁵⁷ 

In Fig. 7, we present the results from the finite lead approximation, as introduced in Sec. V B. Typically, the necessary length of the finite lead is less than 50 times the conductor length, $NxL,NxR≲50×L$⁠, to achieve a relative error less than 2% throughout most of the energy spectrum. This necessitates a sparse matrix with the 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³ to 10⁴). Moreover, the calculation of self-energies is also much simpler and faster than that 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 KPM. In Fig. 8, the computation time of these two methods is plotted as a function 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 contain information on the full energy spectrum, this improved KPM will be even superior when one needs data on more energy points. Furthermore, in the context of simulating irregular shaped conductors with an irregular structure of the Hamiltonian matrix, since the calculation cannot be reduced to a layer-to-layer recursive one,^13,47 the method we proposed will be applicable.

In summary, we introduced the Chebyshev polynomial method to the Landauer–Büttiker formula of the two-terminal transport device, by 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 self-consistent 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.

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 lead–conductor coupling,^51–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 Chebyshev forms to linear responses of other degrees of freedom of quantum transport structures.^7,10,12,59

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

We thank Professor S.-J. Yuan (Wuhan University) for beneficial discussions. This work was supported by the National Natural Science Foundation of China under Grant Nos. 11774336 and 61427901. Y.-Y.Z. was also supported by the Starting Research Fund from Guangzhou University under Grant No. RQ2020082.