Transport through Quantum Anomalous Hall Bilayers with Lattice Mismatch

We theoretically investigate quantum transport properties of quantum anomalous Hall bilayers, with arbitrary ratio of lattice constants, i.e., with lattice mismatch. In the simplest case of ratio 1 (but with different model parameters in two layers), the inter-layer coupling results in resonant traversing between forward propagating waves in two layers. In the case of generic ratios, there is a quantized conductance plateau originated from two Chern numbers associated with two layers. However, the phase boundary of this quantization plateau consists of a fractal transitional region (instead of a clear transition line) of interpenetrating edge states (with quantized conductance) and bulk states (with unquantized conductance). We attribute these bulk states as mismatch induced in-gap bulk states. Different from in-gap localized states induced by random disorder, these in-gap bulk states are extended in the limit of vanishing random disorder. However, the detailed fine structure of this transitional region is sensitive to disorder, lattice structure, sample size, and even the configuration of leads connecting to it, due to the bulk and topologically trivial nature of these in-gap bulk states.

As an effect in 2D, the minimal model of QAH effect can be realized on a monolayer lattice with two orbitals at each site [1,13]. On the other hand, bilayer systems have been shown to be a platform for richer phenomena. For example, a graphene bilayer, even regularly stacked as a periodic lattice, has possessed electronic structures and transport properties significantly different from those of a monolayer one [14]. By twisting bilayers into generic non-commensurate or quasi-periodic lattices, more nontrivial physics emerge, including strong correlation, superconductivity and nontrivial topology [15][16][17][18][19].
Quasi-periodicity is a delicate pattern between periodicity and disorder, which may give rise to interesting behaviors of electronic wavefunctions, even in the single particle picture [20,21]. For example recently, plenty of interesting localization properties have been found in one-dimensional (1D) lattices with a quasi-periodic potential, such as welldefined mobility edges [22,23], and a critical region consisting of critical states [24]. Insightful predictions on other rich physics in 1D quasi-periodic lattices are proposed, for example, topology, non-Hermeticity, superconductivity and superfluidity [25][26][27][28]. Some interesting experimental realizations of quasicrystalline physics have also been discussed recently [29,30].
In this manuscript, for the QAH effect, we consider another method of creating quasi-periodicity, by coupling two QAH layers with arbitrary ratios of lattice constants. Experimentally, this can be realized by using state-of-theart fabrication technologies of topological hetero-structures [31][32][33], ultracold atomic [34][35][36] or photonic systems [37]. Firstly, by comparing the numerical and analytical results in the simple case of commensurate bilayer, we confirm that the fundamental physics underlying this system is governed by two coupled 1D Dirac equations, if no bulk states have been considered. Then, the quantum transports are investigated for different lattice constant ratios, by using a tight-binding model that naturally includes both bulk and edge states. The central plateau of the quantized conductance is not affected since the coupling has not been able to close and re-open the bulk gap. On the other hand, the phase boundary of this quantized conductance plateau consists of a fractal transitional region, with small and fractured islands of quantized and unquantized conductances penetrating into each other. This picture is distinct from that of the commensurate system, where the boundaries of the topological region are clear. Moreover, we find that these unquantized states in the transitional region are extending across the sample. This is also distinct from that of the disordered topological systems (topological Anderson insulators) where they are localized. As a result, the details of this transitional region is sensitive to disorder, shape of the sample, and even the configuration of leads connected to it. The effects from nonzero disorder and varying coupling strength are also discussed.
The rest of the paper is organized as follows. In Sections II, we introduce the model and the method we use to describe the electronic transport properties of bilayer QAH system. Section III is dedicated to compare the results of commensurate lattices obtained by two different approaches: a tight-binding model using the Landauer-Büttiker formalism and a mode-matching calculation in the continuum Dirac-like Hamiltonian approximation. Section IV is dedicated to present numerical results of the generic case with a full tight-binding simulation: mismatching lattice constants in two layers. Finally, we conclude in Sec. V summarizing our main results.

II. THE MODEL AND METHOD
The general form of a spinless bilayer QAH system can be expressed as where H ℓ (ℓ = 1, 2) is the Hamiltonian for the ℓ-th layer and H c is the coupling between them. Here, we choose each layer to be the spin up component of the Bernevig-Hughes-Zhang (BHZ) model [13] defined on a square lattice, with one s orbital and one p orbital on each site. In this manuscript, we fix the lattice constant of the first layer d 1 = 1 as the length unit, while vary that of the second layer d 2 continuously. The layer Hamiltonian H ℓ in the momentum space is where c † ℓ;kα (c ℓ;kβ ) creates (annihilates) an electron with wave number k and orbital α ∈ {s, p} in the ℓ-th layer. For the BHZ model, h ℓ;αβ is a 2 × 2 matrix defined as [13] h ℓ (k) = d 0 with σ x,y,z the Pauli matrices acting on the orbital space {s, p}. We assume that these model parameters of each layer can be tuned independently. In the absence of the inter-layer coupling H c , the ℓ-th layer is in the QAH state if B ℓ · M ℓ > 0 [38]. The real space version of the layer Hamiltonian H ℓ = ij,αβ h ℓ;αβ (i, j)c † ℓ;iα , c ℓ;jβ can be obtained from Eq.(2) and (3) by performing a straightforward inverse Fourier transformation c ℓ;kβ = 1 With the real space layer Hamiltonian at hand, we can define the inter-layer Hamiltonian H c . In most of this manuscript (i.e., except Section III), we adopt a simple form as [39,40]  where t is the strength of the inter-layer coupling and d ij is the distance between sites i and j in the first and second layer respectively. To simplify numerical calculations, we only count in inter-layer couplings t dij 0.01t. By varying the lattice constant of the second layer, the lattice mismatch can be felt by the electron through H c .
We concentrate on the transport properties of the QAH bilayer connected to two semi-infinite leads, as shown in figure 1. At zero temperature, the two-terminal conductance can be calculated within the Landauer-Büttiker formalism G = 2e 2 h T [41][42][43][44], where e is the elementary charge, h is the Planck constant and 2 is the spin degeneracy. This transmission T at Fermi energy E can be expressed in terms of Green's functions as [45,46], where the dressed retarded (advanced) Green's function G r(a) is the retarded (advanced) self-energy from lead p(= L, R), i.e., the left and right leads. The local current from site i to j along the bond is [47] where H ij the matrix element of the bare Hamiltonian H, and G n (E) = G r D (E)Γ L (E)G a D (E) is the correlation function. Since it is defined in the linear response regime, we simply take the voltage difference V S − V D between the source and drain leads to be unity.
To understand more physics underlying transport properties, we rely on the density of states (DOS) [48] and the normalized participation ratio (NPR) [49][50][51]. With the real space representation of the Green's function G r (E; i, j), the single particle local density of states (LDOS) ρ(i) at each site i, and the total density of states (DOS) can be calculated respectively as where N is the total number of sites, and the summation is over all sites of the sample. Here the Green's function G r can be the dressed one G whether one wants to include the effects from leads. With the LDOS normalized over the sampleρ(i), the NPR can be defined as follows [49][50][51]ρ The NPR is a significant diagnostic tool to characterize the localization transition. The localized (extended) phases are characterized by NPR = 0 (NPR = 0) in the large N limit [49][50][51].
To visualize and distinguish edge states from bulk ones, the edge-locality marker is introduced as [ where the summation runs over four edges of the bilayer QAH sample. In this work, we take the width of the edge as 10% of that of the bilayer sample.

III. THE COMMENSURATE BILAYER SYSTEM
In the topological phase of the bilayer, although the backscattering is forbidden, the inter-layer (forward) scattering is not. The nontrivial topology only guarantees the robustness of the total transmission through the bilayer, but the details of inter-layer scattering are not clear yet. As a starting point, let us gain some physical insights from the simplest bilayer with identical lattice constants but different model parameters. To further simplify the model so that an analytical continuum model can be easily obtained, we adopt a simpler form of inter-layer coupling in this section as where the inter-layer coupling t [green line in figure 1(a)] only exists between nearest sites from different layers.
In figure 2, we present the band structures of the QAH bilayer in the ribbon geometry, without (a) and with (b) the inter-layer coupling respectively. With remarkably different model parameters in two layers, the velocities of two groups of edge states are different too. In figure 2(b), the inter-layer coupling lifts the trivial degeneracy of edge states at the Γ point. For this commensurate system with the inter-layer coupling (12), the low-energy properties of the QAH edge states can be captured by a Dirac-like continuum model from k · p approximation [53,54]. This effective Hamiltonian has the form [55,56]: where k x is the x component of the momentum relative to the Dirac point, and v i > 0 is the velocity of edge states in the i-th layer. For each layer, only the forward channel should be considered here, due to the forbidden of backscattering. The inter-layer coupling coefficient F depends on the concrete realization of the bilayer. It can be verified that, when the model parameters of the bilayer QAH system satisfy A2 , F just equals the coupling strength t in Eq. (12). In the rest of this section, we adopt the model parameters satisfying this condition. In the low energy limit around k [53,54]. The expressions of eigenvalues and eigenfunctions of this coupled Dirac Hamiltonian are listed in the Appendix.
We use the transport configuration illustrated in figure 1(a), with one layer as the source lead while a decoupled bilayer as the drain lead, so that the intra-and inter-layer transports can be conveniently distinguished. Since the source lead is a monolayer carrying Chern number 1, the total transmission through the central sample will be quantized as 1, If it is in the topological state. Then due to the inter-layer coupling of the sample, an electron from the source lead of monolayer can be transmitted into either monolayer of the drain lead, corresponding to the intra-layer and inter-layer transmissions respectively. It is straightforward to solve this tunnelling problem. Due to one monolayer (with one pair of edge states) as the source, only one active channel (positive-velocity branch of this pair) is propagating in the central QAH bilayer, with the length L and width W . The wave functions in different spatial regions can be written as where k ℓ = E v ℓ , ℓ = 1, 2 at the incoming energy E. In the central bilayer sample the wave functions are linear combinations of the eigenfunctions given in Eq. (A2) in Appendix, with M and N the combination coefficients. Due to the perfect transmission of Dirac-like fermions, there is no reflection here [57][58][59]. We only need to match the wave functions at boundaries because the Dirac-like Hamiltonian is first-order differential [60,61]. From this continuity of the wave functions at the x = 0 and x = L, i.e., boundaries of the central bilayer, we obtain coefficients M , N , t 1 and t 2 . According to current conservation of the Dirac-like Hamiltonian, we have In figure 3, red curves are results calculated from this effective continuum model, for the intra-layer (left column) and inter-layer (right column) transmissions, with the central sample length L = 50 (upper row) and L = 100 (lower row) respectively. For comparison, corresponding results from the tight-binding simulations are also plotted as the black curves. The perfect quantization of the total transmission shown in figure 3 (c) and (f) confirm that the central sample is in the topological phase. Now let us scrutinize details of the layer resolved transmissions T 1→1 and T 1→2 . Curves in both colors agree very well, especially when E → 0. So in the low energy limit, the Dirac continuum approximation contains most of the physics. The intra-layer and inter-layer transmissions show a complementarity: the transmission minima in one configuration coincide with the maxima of the other one, since the sum of them is robustly quantized as 1.
Using the above mentioned analytical formalism, we can also obtain the dependence on the sample length L as where q 1 = F √ v1v2 by Eq. (A2). In figure 4, we can clearly see that for a fixed energy E, the transmission varies periodically with L from . At zero energy. the period π √ v1v2 F can be obtained from Eq. (15), which has an excellent agreement with the tight-binding result shown in figure 4(a).
In Fig. 5, we plot the spatial distribution of local currents [obtained from equation (6)] of the sample at a typical parameter setting (see the caption), which clearly show that the currents are carried by edge states.
In brief, in the topological phase, if the influence from bulk states are ignored, the transmission through the bilayer is quantized. Now the only interesting effects from the inter-layer coupling are just the resonant traversing between forward propagating waves in two layers.

IV. GENERIC RATIOS OF LATTICE CONSTANTS BETWEEN BILAYERS
In the rest of this manuscript, we will come back to the generic case of mismatching lattice constants between two layers, with the distance dependent inter-layer coupling Eq. (4). As mentioned in Section II, the lattice constant of the first layer d 1 is fixed to be the length unit, and that of the second layer d 2 will be varied. Now even in the topological phase, a low energy approximation based analytical treatment will be difficult, because, for example, the inter-layer coefficient F in the effective Hamiltonian (13) will have a very complicated dependence on model parameters and lattice configurations. Moreover, the effective Dirac Hamiltonian (13) only involves edge states and therefore it cannot describe any effects from bulk states, or the transition into topologically trivial phase, which we are focusing in the following. As a result, we will rely on numerical calculations based on the tight-binding model, as introduced in Section II.
The transport configuration will be taken to be a more "natural" setup as illustrated in figure 1(b). Each lead (source or drain) consists of two decoupled monolayers, otherwise the coupling between incommensurate bilayers will break the translation symmetry of the lead and make the self energies incomputable. Since we are not interested in what happens in the leads, this setup will not influence the main physics of the results. The model parameters of each monolayer lead are identical to those of the central sample layer it is connected to, to decrease unwanted scattering on the contact boundaries. The total transmission is a sum of all four possible transmissions between these monolayer leads, i.e., T total = T 1→1 + T 1→2 + T 2→1 + T 2→2 .
A. Phase Diagram on the E − d2 Plane Let us consider a system with where each monolayer is topologically nontrivial with Chern number 1. The topological phase of the bilayer can be identified as the quantized transmission T total = 2 with the contribution from edge states in both layers. When the lattice constant of the second layer d 2 is varied gradually, the bilayer system displays interesting transport properties in the presence of inter-layer coupling. figure 6(a)-(d) shows the phase diagram of the total transmission on the d 2 − E plane with the open boundary condition. In figure 6(a), the most prominent feature is the central plateau of the quantized conductance (red). The width of this central plateau grows larger with increasing d 2 . This can be understood as follows. With the stretching of the second layer, most of the inter-layer bonds becomes longer and therefore the associated coupling t becomes smaller. In other words, the effective "average" coupling between two layers becomes weaker. As a result, the bilayer somehow tends to approach the decoupling case, with a large bulk gap around E = 0 where the edge states live.
More interesting behaviors emerge on the boundary region of this central plateau of quantization. In the parameter space of a periodic lattice, the phase boundary of the topological state is a clear-cut line [13]. Here instead, as shown in figure 6(a) and its partial enlargement (b), the boundary consists of a fractal transitional region with interpenetrating quantized (red) and unquantized (blue) states, especially near the lower tip of the central plateau with d 2 ∼ 1. In this same enlarged region, figure 6(c) and (d) show results for different sample sizes. Although unquantized states (blue) flood in, a large number of small islands of quantization with T = 2 (red) still survive in a regular pattern of distribution. With increasing size [figure 6(c) and (d)], this picture of floods and islands will be more fragmented, with quantized and unquantized parts penetrating into each other with more fine details. To shed further light on the fractal nature of the transitional region, we calculate the dimension of the curve T total (E) by using a box-counting (BC) algorithm [66]. This algorithm counts the number N of squares of size δ which are necessary to continuously cover the graph of T (E) rescaled to a unit square. We choose an intermediate region (usually called the "scaling region") where scaling is linear in a log-log plot, i.e., where N ∼ δ −d . The slope d in the scaling region is the estimate of the BC dimension [67]. In Table 1, we show the results of the BC dimension for different sample sizes and different d 2 . With increasing size, the BC dimensions get larger, which is consistent with above knowledge obtained from the figure 6 that more finer details emerge. Now we will investigate the property and origin of this picture of transports. Let us first scrutinize this by relating transports with electronic wavefunction behaviors. In figure 6(e), the edge-locality marker B defined in Eq. (11) is plotted on the same E − d 2 plane, and also with figure 6(f) its partial enlargement. Here the LDOS was calculated from an isolated bilayer sample (i.e., without being connected to leads) in order to manifest the intrinsic property of the bilayer structure itself. It is interesting to notice the perfect correspondence between T and B [(a) and (e); (b) and (f)]: a quantized T = 2 corresponds to an edge state with B ∼ 1, while an unquantized T corresponds to a bulk states with B ≪ 1. In other words, the phase boundaries of the central plateau of conductance quantization consist of a transitional region, where a series of inter-crossing bulk states are imbedded in to separate the topological edge states into small islands.
To disclose more details of this picture, we plot two partly cross sections of figure 6(b) along E as the black curve in figure 7. The DOS and NPR are also plotted as the red and green curves respectively. Let us first compare the transmission (black) and the DOS (red). The transmission plateau with T = 2 always corresponds to a vanishing DOS, suggesting an edge state in the bulk gap. The peaks of the DOS correspond to trivial bulk states, which disrupt the transmission quantization. As for the NPR (green curve), that of the trivial bulk states is always larger than that of the edge states. This suggests that these trivial bulk states are quite extended. In figure 7, point A (B) represents a typical localized edge (extended bulk) state.
In order to have a more direct picture of these states, we present the distributions of LDOS in figure 8. The upper (lower) row corresponds to the energy point labeled as A (B) in figure 7. At point A, the wavefunctions are well localized at the edges. At point B, on the other hand, the wavefunctions are distributed in a well extended way throughout the sample bulk. This is consistent with above knowledge obtained from the NPR. As we have seen from figure 7 and figure 8, although these bulk states are narrow, they are quite extended. This is different from those extremely narrow and localized states in the bulk of the topological Anderson insulator [48,68]. The concrete configuration (positions, widths) of these bulk states depends on the details of the lattice sample, e.g., the sample size and lattice constant ratio d 2 /d 1 . Moreover, their transports are even sensitive to the device details, as will be presented in the following. So far in the transport calculations, the chemical potentials of leads µ are set to vary with that in the bilayer sample. In the energy region we have focused, they are in the bulk gap of leads, so that each lead just contributes two robust channel of topological edge states. Now in order to check the robustness of the transitional region, we try to supply more active channels into the bilayer sample. In the calculation, this can be simply realized by changing the chemical potentials µ of the leads into their bulk band. The transmission with this setup on the E − d 2 parameter plane is plotted in figure 9(b). For comparison, the same range of figure 6(a) (with lead chemical potentials in the bulk gap) is re-plotted here as figure 9(a). In figure 9(b), although the profile of those crossing bulk states can still be seen, the transmission is very asymmetric respect to E = 0. In the negative energy region, most of the quantization islands with T = 2 have been submerged by bulk states with T = 2. This can be attributed to the coexistence of edge states and in-gap bulk states. It is known that the coupling to bulk states can severely destroy the quantization and robustness of topological edge states [69]. Furthermore, the transports of bulk states depend sensitively on details of the device due to their spatial extensions over the bulk. Recently, similar quantum transport phase transition from varying lead chemical potential is found in one-dimensional long range systems [70]. To see impacts from leads more clearly, now we present the total transmission as a function of the lead Fermi energy µ (with respect to its half filling level), with all parameters of the central bilayer fixed. Three panels are for three different lead setups in figure 10. Three colors of curves correspond to three typical settings of the central bilayer sample shown in figure 6(a-d): on the central red plateau of quantization, on a red island of quantization, and in the blue sea of unquantization, respectively. In figure 10 (a), model parameters A ℓ , B ℓ , C ℓ , D ℓ and M ℓ are same for the sample and the leads. For a central bilayer sample on the central plateau of quantization (black dashed curve), the transmission is robustly quantized for all µ, suggesting an edge state deep inside the bulk gap. For a central bilayer sample on a red island of quantization (red curve), the transmission is only quantized until µ ∼ 1. After that, it rises up beyond 2, showing a mixing from bulk states. For the last case, a central bilayer sample in the blue sea of unquantization (blue curve), the transmission is not quantized at all, exhibiting typical resonance peaks and valleys of bulk states.
Similar phenomena (perfect quantization on the central island of quantization, partially quantization in the transitional region, and unquantization elsewhere) can be observed from figure 10 (b), with significantly different model parameter settings in the lead, as shown in the figure caption. Furthermore we use another completely different type of leads consisting of 1D chains with nearest hopping 1.52 and the results are presented in figure 10 (c). Now only the robust central quantization plateau (black curve) survive, and the transitional region is completely unquantized due to the strong interruptions from leads with a trivial topology.
The above results lead to a physical picture of incommensurate QAH bilayers briefly illustrated in figure 11(a). Here blue rectangles represent bulk bands and red curves are topological edge states between them. Since the incommensurateness breaks the standard Bloch theorem, here the horizontal axis does not necessarily stand for the conventional wavevector. Instead it can represent any other variable model parameters, e.g., lattice mismatch d 2 /d 1 .
Our above results conclude that incommensurateness between two layers gives rise to a cluster of bulk states (blue curves) penetrating into margins of the bulk gap. Here we intentionally plot them as tortuous curves to stress that their positions and energy widths may be sensitive to model parameters (including coupling of leads as observed just now in figures 9 and 10). These incommensurateness induced bulk states are narrow (but not flat), topologically trivial and spatially extended. We call them in-gap extended states (IGESs). Around the bulk gap center, there are only edge states (red lines), which contribute to a very robust transport, giving rise to the central red plateau of quantization in figure 6(a). In the marginal region, the existence of IGES makes things complicated. If the Fermi energy E R happens to be located in a sub-gap between IGESs, the quantization will survive, which results in a small red island of quantization in figure 6(a)-(d). Otherwise, on another Fermi energy E B with a coexistence of edge states and bulk IGESs, the quantization will be destroyed [69]. This corresponds to blue regions in 6 (a)-(d). Due to their bulk and topologically trivial nature, the concrete configurations (e.g., spatial distribution and width along energy axis) of IGESs will depend sensitively on details of the sample and device, as we have seen in figure 9. The parameter region of IGES corresponds to the transitional region outside the central quantization plateau in figure 6(a)-(d).
For comparison, the picture of purely disordered QAH effect is also shown as figure 11(b). Here, disorder induces in-gap localized states (blue horizontal lines), which are topologically trivial and flat [48]. In fact, they are so flat that their total measurement on the energy axis (sum of sub-band widths) even tends to zero in the thermodynamic limit [71]. Therefore, these in-gap localized states will not affect the robust transports of edge states.

B. Disorder Effect
So far the results were calculated without any disorder, with the electronic properties as a consequence of nontrivial topology and incommensurate lattice. In this subsection, we investigate the role of disorder on this system. The effect of non-magnetic impurities is included by introducing a random potential term to the Hamiltonian, where U i are random numbers uniformly distributed in (−0.5, 0.5) and V is a single parameter to control the disorder strength. Total transmissions in the presence of different disorder strengths are presented in figure 12, for d 2 = 2.5 (left column) and for d 2 = 1.2 (right column) respectively. At zero disorder, these two columns correspond to the central quantization plateau and transitional region respectively. In the former case, figure  12(a)-(d), the quantized transmission plateau is unaffected by weak disorder until V ∼ 6. In the latter case, figure  12(e)-(h), more unquantized dips appear even at weak disorder V ∼ 1. To characterize the stability and integrity of topological edge states in a more quantitative way, we introduce the quantization index defined on a certain energy interval [E 1 , E 2 ], where N Q is the number of quantized data points whose total transmission T total satisfies |T total − 2| < ξ with ξ > 0 a small tolerance error, and N T is the total number of data points within this energy interval. Here we fix ξ = 0.  The developments of the quantization index n Q with increasing disorder strength V at different d 2 are presented in figure 13. For the central quantization plateau, d = 2.5 (black curve), it is perfectly intact as a whole with n Q = 1 until strong disorder V ∼ 3. For the rest three cases, the energy region for calculating equation (13) contains considerable portions of transitional regions in the presence of IGESs. As a result, as shown in figure 13, all three curves of n Q show a rapid response and a non-monotonic dependence on V . They first decrease at weak disorder, which suggests some of those quantization islands in the transitional region are destroyed. This originates from bulk and mutable IGESs, consistent with previous conclusions. Disorder can drive more IGESs into the bulk gap, or widen the existing IGESs. In either case, the interference with bulk IGES will destroy the quantized transport of edge states [69]. Then, interestingly, they rise significantly at an intermediate disorder V ∈ (2, 5), indicating more quantized islands emerging. This can be directly confirmed from some snapshots of conductance at some V shown in figure 14, especially at V = 4.9. Now, many IGESs have been localized and flattened by sufficiently strong disorder and they lose their contribution to quantum transports [48,68,71]. Thus the topological edge states can manifest themselves on more Fermi energies. In other words, an intermediate disorder localizes IGESs and changes the picture from figure 11 (a) to (b). At the very end, everything is localized in the strong disorder limit V > 7, after the final touching of bulk bands and the annihilation of opposite Chern numbers [71,72].

C. Inter-layer Coupling
Before the end, now we study the effects from varying the inter-layer coupling. The developments of the total transmission with the increasing of inter-layer coupling strength at different d 2 are illustrated in figure 15. Let us start from a QAH phase at the inter-layer coupling t = 0. Both panels of figure 15(a) and (b) show that quantized transmission T total = 2 persists at weak inter-layer coupling but is destroyed for large inter-layer coupling. Similar to the above case of varying E, there is also a transitional region outside the central region of perfect quantization, with inter-crossing structures of quantized and unquantized islands. Similarly, they can also be attributed to narrow bulk states as a result of lattice incommensurateness.

V. SUMMARY
In this work, we studied the quantum transports of a bilayer QAH system. In the case of lattice match, we compared the transmissions from the tight-binding approach and the continuum Dirac model. Results from both methods have an excellent agreement especially in low-energy limit. We found that the intra-layer and inter-layer conductances show a complementarity since the sum of them is robustly quantized as 1 and the intra-layer transmission varies periodically with the central sample length. In this simple case, in the topological phase, the inter-layer coupling just leads to the resonant traversing between forward propagating waves in two layers.
In the case of lattice mismatch, the quantum transports are investigated by a full tight binding simulation. In the phase diagram on the E − d 2 plane, the central plateau of quantized conductance corresponds to the bulk gap which has not been influenced by the lattice mismatch. On the other hand, outside the central plateau of quantization, there is a transitional region with small and fractured islands of quantized and unquantized conductances penetrating each other. We identify this region as in-gap extended states penetrating with topological edge states. Details of this transitional region is sensitive to disorder, shape of the sample, and even the configuration of leads connected with it, due to the bulk and topologically trivial nature of these in-gap states. In phase diagram on the E − t plane, the quantized conductance persists at weak inter-layer coupling and there is also a similar transitional region. Our results offer a comprehensive view of transport through the QAH bilayer with lattice mismatch.