Localization Trajectory and Chern-Simons axion coupling for Bilayer Quantum Anomalous Hall Systems

TL;DR

This study explores how bilayer quantum anomalous Hall systems with high Chern numbers behave under strong disorder, revealing complex localization paths and significant magnetoelectric effects linked to the Chern-Simons axion coupling.

Contribution

It uncovers the rich localization trajectories and stable intermediate states in bilayer QAH systems, extending understanding beyond single-layer models and standard scaling theories.

Findings

Presence of stable intermediate C=1 state tunable by parameters

Multiple localization scaling patterns beyond Pruisken picture

Remarkable peak in Chern-Simons axion coupling during localization

Abstract

Quantum anomalous Hall (QAH) multilayers provide a platform of topological materials with high Chern numbers. We investigate the localization routes of bilayer QAH systems with Chern number C = 2 under strong disorder, by numerical simulations on their quantum transport properties and the Chern-Simons axion coupling. Compared to the single layer counterpart with C = 2, the localization trajectories present much richer behaviors, for example, the existence of the stable intermediate state with C = 1 can be tuned by model parameters. This state was always unstable in the single layer case. Furthermore, the two parameter scaling trajectories also exhibit multiple patterns, some of which were not captured by the standard Pruisken picture. During the process towards localization, the Chern-Simons axion coupling shows a surprisingly remarkable peak which is even higher and sharper in the large size limit. Therefore the disordered bilayer QAH system can be a good candidate for this nontrivial magnetoelectric coupling mediated by orbital motions.

Equations

$${\mathcal{L}}_{\theta }=\frac{e^2}{h}\frac{\theta }{2\pi }\mathbf{E}\cdot \mathbf{B},$$

$${\mathcal{H}}_{\mathrm{bi}}={\mathcal{H}}_1+{\mathcal{H}}_2+{\mathcal{H}}_c,$$

$${\mathcal{H}}_L=\sum _{\mathbf{k},\mathbf{\alpha \beta }}{c}_{L;\mathbf{k}\alpha }^{†}{H}_{L;\alpha \beta }(\mathbf{k})c_{L;\mathbf{k}\beta },$$

$$H_L(\mathbf{k})$$

$${\epsilon }_L(\mathbf{k})$$

$$d_L^1(\mathbf{k})$$

$$d_L^3(\mathbf{k})$$

$$C_{L}=\left\{\begin{array}[]{lll}+1,&\quad 0<M_{L}/2B_{L}<2\\ 0,&\quad M_{L}/2B_{L}<0\\ -1,&\quad\mathrm{otherwise}.\end{array}\right.$$

$$\mathcal{H}_{c}=\sum_{i\alpha}\big{(}tc^{\dagger}_{1;i\alpha}c_{2;i\alpha}+\mathrm{H.c.}\big{)},$$

$$\sum _{L;i\alpha }{c}_{L;i\alpha }^{†}{U}_{L;i}{c}_{L;i\alpha }$$

$$G=\frac{e^2}{h}\mathrm{Tr}\left[{\Gamma }_S{G}^r{\Gamma }_D{G}^a\right],$$

$$\sigma_{xx}(E_{n})=\frac{e^{2}}{h}\pi\rho(E_{n})\frac{\partial^{2}E_{n}}{\partial k_{x}^{2}}\Big{|}_{k_{x}=0},$$

$${\sigma }_{xy}(E_n)$$

$${\mathbf{F}}_n$$

$${\mathbf{A}}_n$$

$${\bar{\alpha }}_{\mathrm{CS}}=\frac{e^2}{h}\frac{{\theta }_{\mathrm{CS}}}{2\pi },$$

$$\theta_{\mathrm{CS}}=-\frac{1}{4\pi}\int d\bm{k}\epsilon_{ijl}\mathrm{Tr}\big{(}A^{i}_{\bm{k}}\partial_{k_{j}}A^{l}_{\bm{k}}-i\frac{2}{3}A^{i}_{\bm{k}}A^{j}_{\bm{k}}A^{l}_{\bm{k}}\big{)},$$

$${\theta }_0=-8{\pi }^2\mathrm{Im}\mathrm{Tr}[PxPyPz],$$

$$\frac{{\theta }_0}{N_x{N}_y}\to {\theta }_2,\mathrm{for}{N}_x,N_y\to \infty$$

$${\theta }_0(\mathbf{r})=-8{\pi }^2\mathrm{Im}{\mathrm{Tr}}_{\mathbf{r}}[PxPyPz],$$

$$U_{L;i}=W\cdot {\epsilon }_{L;i},$$

Full text

I I. Introduction

The research of topological materials recovers the interests of searching the quantum anomalous Hall (QAH) effects, i.e., Chern insulators without an external magnetic fieldHaldane . QAH effects in different materials have been theoretically proposed QAH_CXLiu ; QAH_HJZhang ; QAH_QZWang ; QAH_RYu ; QAH_SCWu and experimentally observed QAH_Exp1 ; QAH_Exp2 ; QAH_Exp3 ; QAH_Exp4 ; QAH_Exp5 . The Chern insulator, as a topological state, is robust against weak disorder, manifested as a plateau of Hall conductance $\sigma_{xy}=\frac{e^{2}}{h}C$ with $C$ the Chern number, and a vanishing longitudinal conductance $\sigma_{xx}$Klitzing1980 ; Thouless1982 ; QNiuRMP . Nevertheless, it can still be localized into Anderson insulators with $C=0$ ultimately, by sufficiently strong disorder. This can be regarded as a disorder induced topological transition with Chern number from nonzero to zeroQHEDisorder ; YYZhang2012 ; ZG2015 . On the other hand, this trajectory towards localization can also be understood in the context of scaling and renormalization group (RG), i.e., the Pruisken formalismPruisken1985 ; Huckestein1995 , where all the quantized Hall conductances (associated with all nonzero Chern numbers) correspond to stable fixed points attracting all nearby RG flows in the conductance plane $(\sigma_{xy},\sigma_{xx})$. Besides the analytical RG treatment, these trajectories along scaling transformation can also be observed by experimentsRG_Exp2 ; Checkelsky2014 ; SGrauer2017 or simulated by numerical calculationsXue2013 ; Werner2015 ; ZG2015 ; YXXing2018 .

QAH states with different Chern numbers are topologically distinct from each other. The localization trajectory of the monolayer (two-band) QAH state with Chern number $C=2$ has been numerically investigatedSong2015 ; ZG2015 . The novelty is that the intermediate state with $C=1$ is not stable during the process towards localization, no matter how the model parameters are tuned. Microscopically this has been attributed to a $C=1$ state rooted in broad statistical tails, which cannot manifests itself after statistical averaging and size scalingZG2015 .

In additional to transport properties, the orbital mediated magnetoelectric coupling attracts many attentions in topological materials. The axion electrodynamicsWilczek1987 introduces a non-Abelian Lagrangian termSpaldin2005 ; Fiebig2005 ; Hehl2008 ; Essin2009 ; Olsen2017 ,

[TABLE]

where $\theta$ is a dimensionless parameter which is determined by the electronic structure of the material, and $\bm{E}$ and $\bm{B}$ are the electric and magnetic fields respectively. This quantity reflects a deep magnetoelectric connection mediated by the electronic orbital motions in solids, especially remarkable in topological materials. Due to the involvement of three dimensional motions [see Eq. (13) below], a bilayer system, which will be investigated in this paper, is the simplest lattice with a nonzero Chern-Simons axion couplingOlsen2017 .

An important knowledge from (quasi) two-dimensional (2D) systems is that a bilayer system may possess quite different properties from its monolayer counterpartGrapheneRMP ; QHE_Bilayer . Few layer QAH systems have been experimentally fabricated recently by molecular beam epitaxyQAH_Multilayer . The bilayer QAH system with $C=2$ we will discuss is essentially different from a monolayer one in at least two aspects. First, it is a four band model, or, a model with two pairs of bands, with each pair contributing one Chern number. Second, each pair corresponds to one layer which is distinguishable in the real space, and their coupling is tunable.

In this paper, we will investigate the trajectory towards localization of such a bilayer system with increasing disorder, by performing numerical simulations on the two terminal conductance, two-parameter scaling, and the axion magnetoelectric coupling. We find that, there is a stable Hall conductance plateau with $C=1$, in a region of intermediate disorder strength. For some model parameters, this state may be hardly visible (imperfect plateau of $\sigma_{xy}$ and nonzero $\sigma_{xx}$ ) at finite sample size, but scaling shows that it is still on the way towards $(\sigma_{xy}=1\times\frac{e^{2}}{h},\sigma_{xx}=0)$. The existence of such a stable intermediate Chern state has been ruled out for single layer QAH systemsSong2015 ; ZG2015 . Through numerical scaling, we also find that some trajectories cannot be perfectly fitted into the Pruisken flows, especially in the limit of strong inter-layer coupling. Finally, the disorder induced collapse of the quantum Hall plateau is shown to correspond to a peak of axion magnetoelectric coupling.

II II. The Model

The general form of a spinless bilayer QAH system can be expressed as

[TABLE]

where $\mathcal{H}_{L}$ ($L=1,2$) is the Hamiltonian for the $L$-th layer, and $\mathcal{H}_{c}$ is the coupling between them. Here, we choose each layer to be the spin up component of the Bernevig-Hughes-Zhang (BHZ) modelBernevig2006 defined on a square lattice, with one $s$ orbital and one $p$ orbital on each site. This is one of the minimum models of the QAH state. In the $k$ space, the layer Hamiltonian reads

[TABLE]

where $c^{\dagger}_{L;\bm{k}\alpha}$ ($c_{L;\bm{k}\alpha}$) creates (annihilates) an electron with wavenumber $\bm{k}$ and orbital $\alpha\in\{s,p\}$ in layer $L\in\{1,2\}$. For BHZ model, $H_{L;\alpha\beta}(\bm{k})$ is a $2\times 2$ matrix asBernevig2006

[TABLE]

where $\sigma_{i}$ are the Pauli matrices acting on the orbital space $\{s,p\}$. We have included a layer dependent momentum shift $\bm{q}_{L}=(q_{L}^{x},q_{L}^{y})$ of the band structure in the Brillouin zone, so that two groups of edge states (originated from two layers) can be better distinguished visually, as can be seen in Fig. 1 (a). We have checked that this shift does not affect the transport properties investigated in the following. The real space version of the layer Hamiltonian $\mathcal{H}_{L}=\sum_{ij,\alpha\beta}c^{\dagger}_{L;i\alpha}H_{L;\alpha\beta}(i,j)c_{L;j\beta}$ can be obtained from Eqs. (3) and (4) by performing a straightforward inverse Fourier transformation $c_{L;\bm{k}\beta}=\frac{1}{\sqrt{V}}\sum_{j}c_{L;j\beta}e^{-i\bm{k}\cdot{\bm{x}_{j}}}$, where $i$ is the site index. The size of such a finite sample is characterized by the length $N_{x}$ and the width $N_{y}$ measured in units of the lattice constant.

In the absence of the inter-layer coupling $\mathcal{H}_{c}$, the band structure and Chern number of each layer can be tuned independently by varying the parameters. For example, the band gap is $2|M_{L}|$, and the layer Chern number

[TABLE]

The inter-layer coupling is considered to be the simple form in real space as

[TABLE]

where $t$ is the strength of the coupling. Without inter-layer coupling, the Chern number of this “bilayer” system is trivially $C=C_{1}+C_{2}$. Consequently, if the adiabatic turning on of the inter-layer coupling term (6) does not close the bulk gap, the Chern number $C$ of the bilayer system will not change. Throughout this paper, we focus on this case with $C=2$, constructed from coupling two layers with $C_{L}=1$.

In Fig. 1, we present the band structures of the bilayer systems in the ribbon geometry, with two typical groups of model parameters. Two groups of edge states (red lines) in the bulk gap can be clearly seen, reflecting the Chern number $C=2$ at half filling. To have a direct visual comparison with the case of a single layerSong2015 ; ZG2015 , we have appropriately tuned the model parameter $\bm{q}_{L}$ [see Eq. (4)], so that two groups of edge states (originated from two layers respectively) are centered at different locations with two local gaps respectively. Furthermore, the magnitudes of these two local gaps can be similar [Fig. 1 (a), with $A_{1}=1.0$] or different [Fig. 1 (b) $A_{1}=0.3$]. In the following, we will see that these two patterns correspond to different quantum transport properties.

The effect of disorder is modeled in real space by adding a random onsite potential term

[TABLE]

to the Hamiltonian $\mathcal{H}_{\mathrm{bi}}$, where $U_{L;i}$ are independent random numbers uniformly distributed in $(-\frac{W}{2},\frac{W}{2})$, and $W$ is the disorder strength.

III III. Methods

At zero temperature, the two-terminal conductance of a finite sample can be expressed by Green’s functions asDatta ; Khomyakov2005 ; Caroli2005

[TABLE]

where $G^{r/a}$ is the retarded/advanced Green’s function, and $\Gamma_{S(D)}=i(\Sigma_{S(D)}^{r}-\Sigma_{S(D)}^{a})$ with $\Sigma_{S(D)}^{r/a}$ being retarded/advanced self energies due to the source (drain) lead, respectively. For a Chern insulator, this two-terminal conductance just equals the Hall conductance, $G=\sigma_{xy}$Datta . To plot the two-parameter scaling trajectories, longitudinal and transverse conductances should be calculated respectively, which will be described below.

For a finite 2D sample with disorder, the concept of the wavenumber $\bm{k}=(k_{x},k_{y})$ can be restored if twisted boundary conditions are adopted along both directionsQNiu1984 ; QNiuRMP ; QHEDisorder ; YYZhang2012 ; ZG2015 . This is simply equivalent to treating this finite sample as the supercell of an infinite superlattice. In this sense, hereafter, we call these finite samples as supercells and the energy bands of this superlattice as subbands. Now, the longitudinal conductance $\sigma_{xx}$ of the supercell is calculated as the Thouless conductanceThouless1972 ; Braun1997

[TABLE]

where $\rho$ is the density of states, and $E_{n}$ is the subband at the Fermi energy.

In this context of disordered supercell, the (intrinsic) Hall conductance $\sigma_{xy}$ is found to be proportional to the integration of the Berry curvature $F^{z}_{n}$ over the Brillouin zone (BZ) under the Fermi energy $E_{F}$Simons1983 ; Kohmoto1985 ; QNiuRMP ; Werner2015 ; AHERMP

[TABLE]

where $f_{0}(E)$ is the Fermi distribution at zero temperature, and \big{|}n(\bm{k})\big{\rangle} is the normalized wave function of the $n$-th subband which satisfies H(\bm{k})\big{|}n(\bm{k})\big{\rangle}=E_{n}(k)|n(\bm{k})\rangle. If $E_{F}$ is within the bulk gap, $\sigma_{xy}=\frac{e^{2}}{h}C$ with $C$ the total Chern number of occupied (sub)bands, where $\sigma_{xy}$ carried by the edge states and $C$ determined by the bulk bands are related through the bulk-boundary correspondenceThouless1982 ; QNiuRMP . The Hall conductance $\sigma_{xy}$ may not be quantized if some subbands are only partly filledQNiuRMP ; AHERMP . This relation between Berry curvature and Hall conductance has been used to design and control the anomalous Hall effect in magnetic Heusler compoundsKManna2018 . Our numerical evaluations of Eqs. (9) and (10) are based on the methods used in previous worksFukui2004 ; Werner2015 ; ZG2015 . Hereafter, all conductances will be expressed in units of $\frac{e^{2}}{h}$.

In the presence of disorder, an ensemble average of these transport quantities over disorder configurations should be performed to obtain physically meaningful results, for a certain group of model parameters and supercell size. Then, by enlarging the supercell size, the ensemble averaged conductance vector $\bm{\sigma}\equiv(\sigma_{xy},\sigma_{xx})$ forms a flow in this 2D parameter spaceWerner2015 ; Huckestein1995 ; ZG2015 . Throughout this paper, we will call this a “scaling flow” (or a “scaling trajectory”), while that calculated from RG theoriesPruisken1985 an “RG flow” (or an “RG trajectory”). These flows offer the information of transport properties in the thermodynamic limit.

Now we introduce the calculation of orbital magnetoelectric coupling. Due to the close relation between $\theta$ in Eq. (1) and the linear orbital magnetoelectric response $\alpha_{ij}=(\partial P_{i}/\partial B_{j})_{\bm{E}}=(\partial M_{j}/\partial E_{i})_{\bm{B}}$, with $\bm{P}$ and $\bm{M}$ the macroscopic polarization and magnetizationOlsen2017 , it is adequate to focus on the Chern-Simons axion (CSA) couplingEssin2009 ; Olsen2017

[TABLE]

where $\theta_{\mathrm{CS}}$ is a part of $\theta$ with a $2\pi$ ambiguity. For an infinite three-dimensional (3D) crystalline material, the CSA coupling can be expressed asXLQi2008

[TABLE]

where $A^{i}_{\bm{k}}$ is the matrix of Berry connection between subbands in direction $i$, and the trace is over occupied bands.

Moreover, the expression of the CSA coupling has been carefully generalized to an infinite 2D slab as $\theta_{2}$, and to a finite (zero dimensional) crystal as Sinisa2009 ; Malash2010 ; Marzari1997 ; Olsen2017

[TABLE]

where $P=\sum_{n}|\psi_{n}\rangle\langle\psi_{n}|$ is the projection operator onto the occupied subspace in the ground state. It has been verified that in the thermodynamic limitOlsen2017 ,

[TABLE]

This $\theta_{0}$ defined by Eq. (13) is convenient for treating disordered and finite samples which we are interested here. Moreover, $\theta_{0}$ is demonstrated to be a more fundamental definition of CSA than those in $k$ spaceOlsen2017 . The CSA coupling [Eqs. (12) and (13)] involves responses in all three directions and thus it is meaningless for single layer systems.

IV III. Results

IV.1 A. Appearance of the intermediate state with $C=1$

A Chern insulator (with $C\neq 0$) manifests itself as a robust Hall conductance plateau at weak disorder, but it can be localized at strong disorder. For a monolayer QAH system (with two bands) with Chern number $C=2$, it has been found that, the intermediate state with $C=1$ is not a stable plateau on the route towards localizationSong2015 ; ZG2015 . Here we try to test the stability of the intermediate state $C=1$ of the bilayer system (3). Due to the adjustability of each constituent layer, we can make these two layers (especially their gaps) quite asymmetric to each other, by tuning model parameters distinct on two layers.

To have a first glance at the quantum transport property in the presence of disorder, we calculate the two-terminal conductance by using Eq. (8), which is equal to the Hall conductance for a Chern insulatorDatta . In Fig. 2, the disorder averaged conductance $G$ is plotted as a function of the disorder strength $W$, with different $A_{1}$ (tuning layer gaps), or different $t$ (tuning the inter-layer coupling). Fig. 2 (a) is the results with fixed $t=0.3$ but different $A_{1}$. When $A_{1}=1.0$ (blue line), corresponding to the case of similar local gaps shown in Fig. 1 (a), the conductance decays from 2 to 0 persistently, without a visible intermediate plateau with $G=1$, similar to the monolayer caseZG2015 . However, with decreasing $A_{1}$, so that two local gaps are more and more different [see Fig. 1 (b)], an intermediate plateau with $G=1$ emerges gradually. This plateau with $G=\sigma_{xy}=1$ seems shorter if the inter-layer coupling $t$ is very large, as shown in In Fig. 2 (b). On the other hand, in the case of similar local gaps ($A_{1}=1$), as one can imagine, changing $t$ will never give rise to a stable plateau with $\sigma_{xy}=1$.

IV.2 B. Scaling flows

In Pruisken’s pioneering work based on field theoryPruisken1985 , the RG flows of the integer quantum Hall effect on the $\bm{\sigma}\equiv(\sigma_{xy},\sigma_{xx})$ plane is found to form a periodic pattern, which is illustrated as grey curves in Figs. 3 and 4. For these RG flows, $\bm{\sigma}=(n,0)$ are stable fixed points attracting all surrounding flows, corresponding to the robust $n$-th QHE plateau states. On the other hand, $\bm{\sigma}=(n+\frac{1}{2},\sigma^{*}_{xx})$ (with $\sigma^{*}_{xx}\sim 1$) are saddle points (unstable fixed points) practically repelling most surrounding flows. There are main semicircles connecting neighboring stable and saddle points. It is observed that this nested hierarchical structure of fixed points can be described by the modular group $\Gamma(1)$, which is infinite, discrete, and non-Abelian. Indeed, a very wide range of experimental quantum Hall data of scaling flows have been found to perfectly match just a few universality classes, and each class corresponds to a certain subgroup of the modular groupKSOlsen2018 . In this context, the above Pruisken picture corresponds to the class of $\Gamma_{\mathrm{T}}$KSOlsen2018 .

Another class of RG flow, $\Gamma_{\mathrm{R}}$, is just a size doubling of the above $\Gamma_{\mathrm{T}}$ flow patternKSOlsen2018 , with even Chern number states $\bm{\sigma}=(2n,0)$ as stable fixed points, and with odd Chern number states $\bm{\sigma}=(2n+1,2\sigma^{*}_{xx})$ as saddle points [See the blue curves in the inset of Fig. 3 (a)]. Numerical simulations have shown that the scaling flows of the monolayer QAHE with $C=2$ match this $\Gamma_{\mathrm{R}}$ pattern, with flows going away from $\bm{\sigma}=(1,2\sigma^{*}_{xx})$ZG2015 .

In order to test the modular symmetry of our bilayer model, we plot the supercell size driven scaling flows on the $\bm{\sigma}$ plane, by numerical simulations of conductances using Eqs. (9) and (10), as in previous worksWerner2015 ; ZG2015 ; YXXing2018 . We remind again that in this paper, the flows from the RG techniques and from the conductance simulations are termed as “RG flows” and “scaling flows” respectively. The results of scaling flows for different $A_{1}$ are displayed in Fig. 3 as red arrows, where each data point is an average over at least 10000 disorder configurations (30000 configurations around transitions). For comparison, RG flow curves of the Pruisken patternPruisken1985 are also illustrated in Fig. 3. These RG flows are just plotted by fitting the original RG equationsPruisken1985 with a few free parameters (including zooming factors in two directions), instead of by a rigorous RG treatment of our bilayer model. Due to the critical slowing down and large fluctuations near the transitions, it is rather difficult to have a precise and global reproduction (or prediction) of the RG flows from numerical scaling flowsHuckestein1995 . However, intuitive insights can still be drawn.

Fig. 3 (a) is the scaling flows for the case of $A_{1}=1$, with two local gaps similar. Here, we simultaneously fit the same group of scaling flows (red arrows) to the RG flows of $\Gamma_{\mathrm{T}}$ class (grey curves in the main panel) and to those of $\Gamma_{\mathrm{R}}$ class (blue curves in the inset), respectively.

The profile of this large semicircle of red arrows seems to follow the RG pattern of the modular symmetry $\Gamma_{\mathrm{R}}$ [blue curves in inset of Fig. 3 (a)], i.e., with $\bm{\sigma}=(1,2\sigma_{xx}^{*})$ a saddle point KSOlsen2018 , as in the monolayer case of $C=2$Song2015 . However, something is different if we scrutinize the orientations of the red arrows around $\sigma_{xy}=1$: Although far away from $\bm{\sigma}=(1,0)$, they show an obvious tendency towards $\sigma_{xy}=1$, which is opposite to the $\Gamma_{\mathrm{R}}$ RG flows (blue curves in the inset) going away from $\sigma_{xy}=1$ on the main semicircle. Instead, the orientations of these red scaling arrows around $\sigma_{xy}=1$ comply more with the $\Gamma_{\mathrm{T}}$ RG flows in the main panel, within the waterfall converging towards $\bm{\sigma}_{1}\equiv(1,0)$. In brief, the scaling flows around $\sigma_{xy}=1$ are more likely to match the RG flows of class $\Gamma_{\mathrm{T}}$ instead of class $\Gamma_{\mathrm{R}}$, although some of the scaling flows are not on the main semicircles.

Here in Fig. 3 (a), we meet an example of badly observable $C=1$ state in the modular symmetry class $\Gamma_{\mathrm{T}}$, manifested as non-quantized $\sigma_{xy}$ and nonzero $\sigma_{xx}$ at finite sample size. This can be attributed to some residual bulk states in the system [see following analysis on Fig. 5 (b)]. Since this state is still on the route towards the stable fixed point $\bm{\sigma}=(1,0)$, a quantized Hall plateau should be observed in an extremely large sample, which may significantly exceed the quantum coherence length in practice. We call such Chern states hardly visible in a mesoscopic sample (although it is stable in the RG sense) as a “weak quantum (anomalous) Hall effect”.

Another issue should be addressed here. In Pruisken’s formalismPruisken1985 , there is always a definite flow vector at any position of the upper half of the $(\sigma_{xy},\sigma_{xx})$ plane, except at isolated fixed points. For an experiment or a numerical simulation, one cannot predict a priori where a certain flow will startHuckestein1995 ; KSOlsen2018 . In other words, the key information of the scaling flows is the orientation configuration at different positions, instead of the positions themselves, because the Pruisken theory does not prohibit or prefer any position for the flows to appear. Indeed, in many experimental results, the scaling flows are not on the main semicirclesHuckestein1995 ; Checkelsky2014 ; KSOlsen2018 .

When $A_{1}$ is decreased, in Fig. 3 (b) to (d) , with two local gaps different, the state with $\sigma_{xy}=1$ becomes more and more attractive for scaling flows, suggesting that it is a stable state during the route towards localization.

The emergence of this stable intermediate plateau with $\sigma_{xy}=1$ is the first finding of this work. Since this state is much more stable and significant when the magnitudes of two local gaps (each accommodating a group of edge states) are extremely different [Fig. 1 (b)], this can be illustrated as disorder induced band inversion YYZhang2012 ; ZG2015 at the smaller local gap. This is completely different from the monolayer QAH effect with $C=2$Song2015 ; ZG2015 . In that case, changing model parameters, e.g., making two local gaps (each of which carries a group of edge states) different will never give rise to a stable plateau with $\sigma_{xy}=1$, and the scaling flows around $\sigma_{xy}=1$ are going away from it. ZG2015 .

Another notable feature can be seen in the most asymmetric case with small $A_{1}=0.3$ [Fig. 3 (d)]: besides straying away from the main semicircles between $\sigma_{xy}=1$ and $\sigma_{xy}=2$, the orientations of these scaling flows seem to deviate from the typical Pruisken RG configurationPruisken1985 ; KSOlsen2018 . This tendency becomes more remarkable when the inter-layer coupling $t$ is increased, as can be seen in Fig. 4 (c) and (d). For example, the orientations of the scaling flows remarkably deviate from those of the standard Pruisken RG flows, especially when $\sigma_{xy}\gtrsim 1.5$. When $\bm{\sigma}$ is high above the semicircle, RG trajectories (grey lines) should be flowing down straightly, but the scaling flows (red arrows) in Fig. 4 (d) are not the case. We conjecture that when the inter-layer coupling is much larger than the intra-layer hopping, the in-plane electronic motions might be suppressed, which leads to some competition induced instability of in-plane electronic transport. This may have an non-ignorable impact on the in-plane transport, which was, however, only understood in a purely 2D regime in traditional RG analysisPruisken1985 .

IV.3 C. Distributions of $\sigma_{xy}$

The above results were based on transport quantities averaged over disorder ensembles. To obtain more insights from another angle, now we investigate the statistical distributions of the Hall conductance at different disorder strength $W$. Fig. 5 is the statistical histograms of $\sigma_{xy}$ at different disorder strength, corresponding to the scaling flows of Fig. 3 (a), with just one semicircle profile of scaling flows and with a weak $\sigma_{xy}=1$ state. On the other hand, Fig. 6 corresponds to the case of Fig. 4 (a), with two semicircles and with a visible $\sigma_{xy}=1$ plateau. When the disorder is weak, $W=1.5$ [panels (a) in both figures], both cases show isolated narrow peaks with quantized Hall conductances, among which the $\sigma_{xy}=2$ peak is the most prominent. This suggests that the system is in a well defined insulating state with an extremely high probability of $\sigma_{xy}=2$. Such a distribution gives rise to a visible $\sigma_{xy}=2$ plateau after ensemble average, as we have seen above.

The most remarkable difference between Figs. 5 and 6 appears at the intermediate disorder strength $W=5$ [panels (b)]. In Fig. 6 (b), similarly, an isolated and prominent narrow peak at $\sigma_{xy}=1$ results in a visible transport plateau of $\sigma_{xy}=1$ which is the stable fixed point in Fig. 4 (a). On the other hand, the pattern is different in Fig. 5 (b). Here although the $\sigma_{xy}=1$ peak is still very high, it is rooted in continuous and broad tails, where $\sigma_{xy}$ is not quantized. This suggests that there exist considerable bulk states at this fermi level, making $\sigma_{xy}=1$ only a “weak topological state” which is hardly visible in a mesoscopic system, as stated above.

IV.4 D. Chern-Simons axion coupling

In this subsection, we explore how CSA coupling strength $\theta_{\mathrm{CS}}$ behaves during the process towards localization. Due to the size relation, Eq. (14), what makes sense for different sizes is the finite size CSA coupling (13) normalized by the sample area, $\theta_{0}/(N_{x}N_{y})$. In Fig. 7, solid curves are the $\theta_{0}/(N_{x}N_{y})$ as a function of disorder strength $W$ for different $A_{1}$. For a comparison with the transport properties, we also plot the corresponding two-terminal conductance [identical to Fig. 2 (a)] as dashed curves. It is interesting to notice that $\theta_{0}$ starts to increase sharply just before the collapse of the conductance plateau with $G=\sigma_{xy}=2$. Then it climbs to a peak after which the tendency of localization begins to dominate. However, the decay of $\theta_{0}$ is slower than that of $G$ after $W>10$. Another distinct feature is that there is also only one peak for all these cases, even for the case of $A_{1}=0.3$ where there are two visible transitions of the Hall conductance: $\sigma_{xy}=2\rightarrow 1$ and $1\rightarrow 0$. These facts suggest that during the process of localization, the variation of CSA coupling is closely related to the disorder induced topological phase transition, but not a simple “one to one” correspondence with transport properties.

To obtain the information in the thermodynamic limit, we present the CSA coupling $\theta_{0}/(N_{x}N_{y})$ for increasing sample sizes in Fig. 8, for $A_{1}=0.3$ and $A_{1}=1$ respectively. In both cases, the peak shapes of CSA become sharper with the increasing of the sample sizes. The peak value $\theta_{0\mathrm{(max)}}/(N_{x}N_{y})$ as a function of the sample size is plotted in the corresponding inset respectively. These peak values do not seem to converge until the largest size we can calculate, $N_{x}\times N_{y}=30\times 30$. With the size axis in the logarithmic scale, one can see that its dependence on sample size $N$ is as fast as (or a little bit slower than) $\theta_{0\mathrm{(max)}}/(N^{2})\sim\ln N$. This scaling growth makes this CSA coupling experimental observation possible in realistic materials. Therefore, besides intrinsic materials previously proposedCoh2011 ; XGWan2012 ; Mogi2017 , we find that multilayer QAH systems with disorder are also candidates for high $\theta$ materials.

The unexpected peak of $\theta_{0}$ prior to localization leads to a question about the underlying physics: how the electronic states are re-organized in the system during this process? To answer this, it is helpful to map $\theta_{0}$ to the real space. The zero-dimensional (finite sample) expression of $\theta_{0}$, Eq. (13), is believed to be a more fundamental definition of the CSA coupling than those $k$-space versions, and is free from the ambiguity of $2\pi$ and the subtleties of boundary conditionsOlsen2017 . Therefore in the spirit of Refs. Resta2011 ; Resta2013 ; Resta2017 , it is reasonable to define the “local marker” of the CSA coupling for a real space site $\bm{r}$ as

[TABLE]

The only difference from Eq. (13) is that the trace $\mathrm{Tr}_{\bm{r}}$ is only over the freedom degrees within the site $\bm{r}$, which is two orbitals in our model (4). The sum of $\theta_{0}(\bm{r})$ over all $\bm{r}$ of the whole finite sample is just $\theta_{0}$.

To have a first impression, let us see the development of the real space distribution of the local marker $\theta_{0}(\bm{r})$ with the increasing of disorder strength $W$, for a single disorder configuration. Here, “a single disorder configuration” means that the onsite potentials (7) of the sample have the form

[TABLE]

where $\{\varepsilon_{L;i}\}$ are a specific realization of random numbers on $2N_{x}N_{y}$ sites of the sample, which are uniformly distributed in $(-0.5,0.5)$, and $W$ is the single variable to control the disorder strength. The results are shown in Figs 9 and 10.

Then afterwards, to smooth out disorder fluctuations and also for a more quantitatively and statistically reliable observation, similar to Ref. Resta2011 ; Resta2013 , we plot the disorder averaged $\theta_{0}(\bm{r})$ along the central cross section of the sample in Fig. 11, with different colors for different disorder strengths.

In the absence of disorder, $W=0$, as shown in Figs. 9 and 10 [panels (a) and (b)], and in Fig. 11, the local marker of the CSA coupling $\theta_{0}(\bm{r})$ on the edges is negative to that in the bulk, leading to a small positive value of total $\theta_{0}$ due to the compensation. Such negative contribution from the edges reflects the peculiar property of the edge statesResta2011 ; Resta2013 . With the increasing of disorder, the typical magnitude of $\theta_{0}(\bm{r})$ on the edges is suppressed much faster than that in the bulk, as can be clearly seen in panels (c) and (d) of Figs. 9 and 10. The statistically averaged data in Fig. 11 also tells the same story. For example in Fig. 11 (b) and (d) associated with the second layers (which contribute $|\theta_{0}(\bm{r})|$ dominantly) in both cases, when $W$ is increased from [math] (black curves) to $2$ (red curves), the local values of $\theta_{0}(\bm{r})$ in the bulk almost remain intact, but those on the edges are reduced by half. This reduction of compensation results in a rise of the total $\theta_{0}$ and a subsequent peak observed in Figs. 7 and 8. Finally in the strong disorder limit, in panels (e) and (f), all types of orbital motions will be localized, and $\theta_{0}(\bm{r})$ tend to vanish uniformly.

V V. Summary

In summary, the localization trajectories of the bilayer QAH systems with $C=2$ are investigated numerically. Distinct to the monolayer counterpart in modular class $\Gamma_{\mathrm{R}}$ZG2015 , this bilayer system is in class $\Gamma_{\mathrm{T}}$, where there is a stable state with $C=1$ before the complete localization. When the bulk gaps associated with two layers are similar, this intermediate state with $C=1$ is weak at finite size but still stable in the extremely large size limit. These pictures are confirmed by numerical simulations of two-parameter scaling flows on the conductance plane. However, some of the scaling flows deviate from the standard Pruisken pattern, especially in the case of strong inter-layer coupling. Microscopically, different performances of this intermediate state with $C=1$ correspond to different types of distribution of the Hall conductance, and the weakness of the intermediate state $C=1$ can be attributed to residual bulk states in the system. These numerical findings can be tested by state of the art experimental techniquesRG_Exp2 ; Checkelsky2014 ; SGrauer2017 and raise new topics for renormalization group theoriesPruisken1985 ; Olsen2017 . Finally, we find the CSA coupling $\theta$ of such systems can be significantly high in the transition window with medium disorder strength. This phenomenon is understood by a real space analysis of the corresponding local marker, by which we find that responses of the edge and the bulk to the disorder are different, and that the weak disorder leads to a suppression of edge-bulk compensation and gives rise to an increase of $\theta_{0}$. Therefore we propose a new way of finding candidate axion insulators.

VI Acknowledgements

This work was supported by National Natural Science Foundation of China under Grant Nos. 11774336, 61427901 and 11374294. YYZ was also supported by the Starting Research Fund from Guangzhou University under Grant No. 69-18ZX10055.