3D Network Model for Strong Topological Insulator Transitions
33D Network Model for Strong TopologicalInsulator Transitions
Jun Ho Son and S. Raghu , Stanford Institute for Theoretical Physics, Stanford University, Stanford, CA 94305, USA Stanford Institute for Materials and Energy Sciences, SLAC NationalAccelerator Laboratory, Menlo Park, CA 94025, USA
August 11, 2020
Abstract
We construct a three-dimensional (3D), time-reversal symmetric generalization of theChalker-Coddington network model for the integer quantum Hall transition. The novelfeature of our network model is that in addition to a weak topological insulator phasealready accommodated by the network model framework in the pre-existing literature, ithosts strong topological insulator phases as well. We unambiguously demonstrate thatstrong topological insulator phases emerge as intermediate phases between a trivial insula-tor phase and a weak topological phase. Additionally, we found a non-local transformationthat relates a trivial insulator phase and a weak topological phase in our network model.Remarkably, strong topological phases are mapped to themselves under this transfor-mation. We show that upon adding sufficiently strong disorder the strong topologicalinsulator phases undergo phase transitions into a metallic phase. We numerically deter-mine the critical exponent of the insulator-metal transition. Our network model explicitlyshows how a semi-classical percolation picture of topological phase transitions in 2D can begeneralized to 3D and opens up a new venue for studying 3D topological phase transitions. a r X i v : . [ c ond - m a t . d i s - nn ] A ug ontents A.1 Conversion formula between scattering matrices and transfer matrices . . . . . . 28A.2 Formula involving two scattering matrices . . . . . . . . . . . . . . . . . . . . . 29A.3 Constructing a global scattering matrix in practice . . . . . . . . . . . . . . . . 30
The Chalker-Coddington network model [1] has been a popular choice of models for studyingnon-interacting integer quantum Hall phase transitions with quenched disorder in two spatialdimensions. This model provides an intriguing semiclassical picture for the integer quantum Halltransitions: Systems with random but smooth chemical potentials are viewed as systems with“mesoscopic” droplets of local integer quantum hall phases in background of a topologicallytrivial insulator, mesoscopic droplets interpreted as local valleys of chemical potential. TheChalker-Coddington model strips away all degrees of freedom except those from edge modesbetween the trivial background and the droplets. Tunneling between edge modes in neighboringdroplets induces topological phase transitions. The critical point occurs when the edge modes“percolate” throughout the whole network. 2riginally devised for the integer quantum Hall transition (class A in the Altland-Zirnbauerclassification [2, 3]), the network model has been generalized to all other symmetry classes inwhich there are non-trivial topological phases [4, 5, 6, 7, 8]; the generalized network modelin each symmetry class describes critical phenomena between two topologically inequivalentgapped phases. Also, the network model was generalized into phase transitions involving atrivial insulator and a 3D weak topological insulator [9]. However, to our knowledge, a con-struction of a 3D network model whose phase diagram includes a strong time-reversal invarianttopological insulator does not exist.Our goal here is to construct such a 3D time-reversal invariant network model. In Sec. 2, wewill motivate and construct a 3D network model that can describe topological phase transitionsinvolving 3D strong topological insulators and point out some interesting dualities in the model.This section will make manifest how the semiclassical percolation picture of the 2D Chalker-Coddington model can be generalized to 3D. In Sec. 3, we present a detailed study of phasediagrams, both with and without quenched randomness. While the phase diagram for theclean system is largely a quantitative confirmation of what we anticipated in Sec. 2, the phasediagram for the dirty system has more intriguing features, most notably exotic phase transitionsfrom strong topological insulators to metallic phases induced by increasing disorder rather thandecreasing disorder as expected from the conventional localization-delocalization transition.
In this subsection, we briefly review the Chalker-Coddington model [1] and its time-reversalsymmetric generalization [7] to motivate our 3D construction and to set up some key notation.
In the Chalker-Coddington model, degrees of freedom are represented as complex number vari-ables ψ ( x,y ) associated with wavefunction amplitudes of fermionic states at site ( x, y ). For odd x , ψ ( x,y ) is associated with fermions moving in + y direction; ψ ( x,y ) with even x is associatedwith fermions moving in − y directions. Hence, fermion modes have chirality structure built inthe setup. Additionally, there are local unitary scattering matrices that endow local relationsamong neighboring ψ ( x,y ) ’s.A 2 × S A ( θ ) (cid:18) ψ ψ (cid:19) = (cid:18) ψ ψ (cid:19) , S A ( θ ) = (cid:18) i cos θ sin θ sin θ i cos θ (cid:19) (1)The key idea is that ψ and ψ ( ψ and ψ ) are associated with incoming modes (outgoingmodes) on the left side of Fig. 1(a); the role of incoming and outgoing modes are reversed on3 a)(b) Figure 1: The illustration of S A ( θ ) in Eq. (1) (a) in the bird’s eye view and (b) in the cross-sectional view.the right side of Fig. 1(a). In either cases, S A ( θ ) unitarily relates two incoming modes and twooutgoing modes.For the sake of later convenience when we introduce our 3D model, we also adopt cross-sectional graphical illustrations of S A ( θ ) in Fig. 1(b). The key observation behind this graphicalnotation is that the pairs ( ψ , ψ ) and ( ψ , ψ ) move in opposite vertical directions from eachother. Hence, we utilize ⊗ and (cid:12) symbols familiar from introductory electromagnetism coursesto represent these pairs in the cross-sectional view.The Chalker-Coddington model has one tuning parameter, which we call θ CC . There arethree types of local relations in the construction of the network model: scattering matrices S A ( θ CC ) that relate fermion amplitudes ψ ( x,y ) with x = 2 i + 1 , i + 2 and y = 4 j + 1 , j + 2(gold squares in Fig. 2), scattering matrices S A (cid:0) π − θ CC (cid:1) that relate fermion amplitudes ψ ( x,y ) with x = 2 i + 2 , i + 3 and y = 4 j + 3 , j + 4 (red squares in Fig. 2), and blue squares in Fig. 2between y = 4 j + 2 and y = 4 j + 3 and between y = 4 j + 4 and y = 4 j + 5. Blue squaresbetween ( x, y ) and ( x, y + 1) induce the following relation: ψ ( x,y ) = e iφ ( x,y ) ψ ( x,y +1) (2)Where e iφ ( x,y ) is a random U(1) phase. Hence, upon crossing the blue squares fermion modesacquire random U(1) phases, and these random phases encode the effect of disorder. Fig. 2(a)4 a)(b)(c) Figure 2: The illustration of the Chalker-Coddington model, (a) from the bird’s eye view (b,c)views from the cross sections. 5hows how these components are arranged in the network model. Alternatively, one can usethe cross-sectional view for the more compact illustration of how the scattering matrices areplaced, as in Fig. 2(b) and (c).From a numerical standpoint, thanks to how scattering matrices are placed in the model,one can naturally construct S y , a scattering matrix that relates fermion amplitudes at y and y + 1 from local scattering matrices illustrated in Fig. 2, for each y . Then, one may use S y ’sto build a global scattering matrix S that relates ψ ’s at the start of the network and ψ ’s atthe end of the network in an effective manner. From this global scattering matrix, one canextract information such as Landauer conductance and topological invariants [10], allowingone to study phase diagrams and critical properties [11] in detail. Alternatively, one canconvert scattering matrices S y ’s to transfer matrices T y ’s and use transfer matrices to computeLyapunov exponents to study localization properties [12, 13]. In Appendix A, we discuss howto construct transfer matrices from scattering matrices and how the global scattering matrixcan be efficiently evaluated numerically.While the Chalker-Coddington model is amenable to numerical methods, one can obtaina great deal of information, without explicit computation, by considering certain limits andinvoking a duality that is inherent to the model. Before discussing phase diagrams, it is usefulto have some intuition on S A ( θ ) in Eq. (1) from the cross-sectional view in Fig. 1(b). If θ = π ,in the setup represented by Fig. 1(b), the left modes and right modes represented as ⊗ and (cid:12) symbols are completely decoupled from each other, and S A ( θ ) is in the perfect transmissionlimit. Meanwhile, when θ = 0, the left modes coming in always reflect to the right modes, andvice versa. Hence, this limit maximizes backscattering between the left modes and the rightmodes. Essentially, in Fig. 1(b), tuning θ can be understood as tuning backscattering betweenfermion modes on the left and fermion modes on the right.Having this intuition at hand, it is useful to consider θ CC = 0 and θ CC = π limit. At θ CC = π , the gold squares in Fig. 2(b) are in the perfect transmission limit, but two modescoupled by a red square maximally backscatters with each other. This has an interestingconsequence: When open boundary condition is imposed along x -direction, the chiral modes at x = 1 and x = 2 L x ( L x represents width of the model, so x = 2 L x is the largest x -index) remaincompletely decoupled from bulk degrees of freedom. Existence of chiral edge modes indicatethat this phase should be identified as a ν = 1 quantum Hall phase. Meanwhile, at the limit θ CC = 0, there is maximal backscattering induced between modes within a unit cell ( x = 2 i + 1and x = 2 i + 2), and such backscattering tunes the system into a trivial insulator.The following duality puts additional constraints on the phase diagram: Now considerimposing periodic boundary condition along x -direction, effectively making the system into acylinder with circumference 2 L x . Then, the substitution ψ ( x,y ) → ψ ( x +1 mod 2 L x ,y ) maps thenetwork with θ CC = θ to the network with θ CC = π − θ . This duality relates the ν = 1quantum Hall phase near θ CC = π to the trivial insulator phase near θ CC = 0. Additionally,if the only two possible extended phases in the phase diagram are a ν = 1 quantum Hall stateand a trivial insulator (This is a natural expectation given that the Chalker-Coddington modelis essentially a 2D non-interacting fermion system with disorder), the critical point betweentwo phases is fixed to be at the self-dual point θ CC = π . Hence, the network model has theself-duality at the quantum Hall transition by construction.6 .1.2 The Time-Reversal Symmetric Generalization To construct the time-reversal symmetric generalization of the Chalker-Coddington model in[7], we will add time-reversed partners to the degrees of freedom in the original model forthe integer quantum Hall transition. Now, ψ ( x,y ) ,s = ↑ , ↓ carry spin labels in addition to spatialindices as those in the original Chalker-Coddington network model. ψ ( x,y ) , ↑ and ψ ( x,y ) , ↓ form aKramer’s pair. ψ ( x,y ) , ↑ may be understood as inherited from the original Chalker-Coddingtonmodel; especially, it has the same chirality structure as ψ ( x,y ) of the original model. In contrast, ψ ( x,y ) , ↓ ’s are associated with fermions flowing in the opposite direction from those representedby the s = ↑ counterpart. Under time-reversal symmetry, ψ ’s transform as: T : ψ ↑ , ( x,y ) → ψ ∗↓ , ( x,y ) , ψ ↓ , ( x,y ) → − ψ ∗↑ , ( x,y ) (3)Now, a basic building block of the time-reversal symmetric network, i.e., taking a role of S A ( θ ) in Eq. (1) in the Chalker-Coddington model, is a 4 × S T ( φ, θ ): S T ( φ, θ ) ψ , ↑ ψ , ↓ ψ , ↓ ψ , ↑ = ψ , ↓ ψ , ↑ ψ , ↑ ψ , ↓ , S T ( φ, θ ) = i cos φ cos θ cos φ sin θ sin φi cos φ cos θ φ cos φ sin θ cos φ sin θ − sin φ i cos φ cos θ − sin φ cos φ sin θ i cos φ cos θ (4)See Fig. 3(a) for the labels. We also adopt the cross-sectional illustration as in Fig. 3(b) aswell.The time-reversal symmetry requires that S T ( φ, θ ) should satisfy: − U T S ∗ T ( φ, θ ) U T = S † T ( φ, θ ) , U T = − − (5)One can explicitly check that the scattering matrix in Eq.(4) satisfies the above condition.In the time-reversal symmetric network model, there are two tunable parameters: φ QSH and θ QSH . The following modifications from the setup for the Chalker-Coddington modelillustrated in Fig. 2 will yield the time-reversal symmetric network model: First, to account fordoubled degrees of freedom, one should add the time-reversed partner modes moving in oppositedirections to the ones already in the setup. Gold squares now represent 4 × S T ( φ QSH , θ
QSH ), while red squares are now taken to be S T ( φ QSH , π − θ QSH ). Finally, the bluesquare at x and between y and y + 1 relates ψ ’s as the following: (cid:18) e iφ ( x,y ) e − iφ ( x,y ) (cid:19) (cid:18) ψ ( x,y ) , ↑ ψ ( x,y ) , ↓ (cid:19) = (cid:18) ψ ( x,y +1) , ↑ ψ ( x,y +1) , ↓ (cid:19) (6)Where e iφ ( x,y ) is a random U(1) phase once again that encode the effect of disorder. Note thatU(1) phase that s = ↑ and s = ↓ modes acquire upon crossing a blue square is chosen to beexactly complex conjugate to each other so that time-reversal symmetry is satisfied.7 a)(b) Figure 3: (a) The setup and the labels used in the definition of Eq. (4) (b) Cross-sectionalview of (a). In both figures and in future figures as well, we use red lines/symbols for spin-upcomponents and blue lines/symbols for spin-down components.8 a) (b)
Figure 4: (a) The naive generalization of the 2D approach to access a strong topological insulatorphase from the 3D network model. Here, one stacks the hypothetical 2D layers that encodephysics of a single Dirac cone and turns on interlayer backscattering between two layers inalternating pattern. (b) The approach explored in this paper. Each layer now has two Diraccones, one for σ z = +1 and one for σ z = −
1; half of the degrees of freedom, marked in red inthe figure, backscatter between layer 2 i + 1 and 2 i + 2, while the other half, colored in blue,backscatter between layer 2 i and 2 i + 1, roughly speaking.When φ QSH = 0, the spin-up sector and the spin-down sector are decoupled from oneanother. The model in this limit is simply two copies of the plain-vanilla network model stackedtogether and represents a physical system with S z conservation. Especially, ( φ QSH , θ
QSH ) =(0 ,
0) corresponds to a trivial insulator upon following the logic presented in reviewing thephases of the original Chalker-Coddington model; in the case ( φ QSH , θ
QSH ) = (0 , π ), the systemis equivalent to ν = 1 and its time-reversed partner ν = − quantum spin Hall insulator . φ QSH (cid:54) = 0 introduces coupling between modes with differentspins. This coupling explicitly breaks S z conservation symmetry and allows one to accessphysics unique to the symplectic symmetry class, such as symplectic critical metal that existsas an extended phase between a quantum spin Hall phase and a trivial phase while the systemremains to be a quantum spin Hall insualtor (a trivial insulator) near θ QSH = π ( θ QSH = 0).
We will start discussion of our 3D network model by giving intuition behind our construction.In the 2D network model construction reviewed in the previous paragraph, one could build thenetwork model by stacking spin-filtered helical fermions in which fermions move in oppositedirections if they have opposite spins and induce backscattering appropriately. The most naiveway to generalize this picture to describe the 3D strong topological insulator would be thefollowing: Assume that there is a 2D network that describes a single, time-reversal symmetricDirac cone . Now, one can imagine introducing strong backscattering between the layer 2 i and2 i + 1 to gap out this Dirac cone. Due to how we introduce backscattering, upon imposingopen boundary conditions along the direction toward which we stack the layers, a single Diraccone will be exposed and will serve as a surface state of the 3D strong topological insulator.See Fig. 4(a) for the illustration. 9ne may question existence of such 2D networks because of the fermion doubling problem.However, we note that the networks described by scattering matrices may avoid usual fermiondoubling problems. One way to see this is to take the scattering matrices as discrete-time evo-lution operators instead of local relations between static fermion amplitudes. This alternativeviewpoint allows one to interpret the network as a Floquet system which is known to evade somedoubling problems [14, 15]. We also observe that in the course of reviewing the 2D networkmodel, we clearly assumed that fermions of interest are chiral, which cannot be captured in alocal 1D lattice Hamiltonian.The prescription given in the previous paragraph is certainly tantalizing, but we are notaware of any time-reversal symmetric networks whose scattering matrices describe a singleDirac cone; at present, we have not succeeded in constructing such a network. Instead, toaccess 3D strong topological insulators from the network model viewpoint, we will take a moreroundabout route: We will imagine stacking 2D time-reversal symmetric network models thatdescribes the physics of two Dirac cones, i.e., the 2D network model for the quantum spin Halltransition, and figure out how to “gap out” half of the degrees of the freedom.To be more specific, let us start with stacks of the 2D networks described by scatteringmatrices given in Eq. (4), explicitly tuning θ = π and φ = 0 to keep things maximally simple.In this setup two fermions with different σ z eigenvalues are independent, and each componentindividually describe the physics of a single Dirac cone. Now, to expose only “half” of thelayer when the open boundary condition is imposed, we will imagine introducing interlayerscattering in the following way: Half of the fermions strongly backscatter between layer 2 i + 1and 2 i + 2, and the remaining half strongly backscatter between layer 2 i and 2 i + 1. Then, whenone imposes open boundary condition, half of the fermions on the top and the bottom layer donot strongly backscatter to neighboring layers, while all other modes are gapped out by strongbackscattering. Hence, roughly speaking only “half” of the layer will constitute gapless edgestates, yielding a strong topological insulator. See Fig. 4(b) for the illustration. We note thata similar idea has been employed in the context of the wire construction to study 3D fractionaltopological insulators [16].There is one caveat in the above intuition. While the 2D quantum spin Hall phase insome limit reduces to a ν = 1 topological phase for s = ↑ stacked with its time-reversed partnerforming a ν = − s = ↓ , such a limit does not exist for 3D strong topological insulators.Especially, to access physics of any 3D strong topological insulators, σ z = ↑ component and σ z = ↓ components should be mixed by inter-layer coupling. Hence, it is difficult to makerigorous claims about how surface states emerge as we spelled out somewhat hand-wavingly inthe previous paragraph. Nevertheless the intuition turns out to be largely correct, and we doobtain a strong topological insulator in the 3D network model built from this insight. Armed with the intuition, now we will explicitly construct our 3D network model whose phasediagram incorporates a 3D strong topological insulator phase. As before, we have complex-numbered variable ψ ( x,y,z ) ,s = ↑ , ↓ , at each site. Following the spirit of the 2D model, ψ ( x,y,z ) , ↑ withodd x and ψ ( x,y,z ) , ↓ with even x are associated with fermions moving in + z -direction, while10 a) (b) Figure 5: (a) The setup and labeling of fermion modes in definition of S L ( θ , θ ) (b) Therepresentation of S L ( θ , θ ) in the cross-sectional view (the view toward the black arrow in (a)) ψ ( x,y,z ) , ↓ with odd x and ψ ( x,y,z ) , ↑ with even x represent fermions flowing in − z -direction. Sincethe fully 3D illustration of the network can be confusing and difficult to visually understand, wewill mostly use the view from the cross-section normal to z -direction, following the graphicalnotations similar to the ones showcased previously in Fig. 1(b), Fig. 2(b), and Fig. 3(b).There are two scattering matrices that we will use as basic building blocks. The first is S T ( θ ) = S T ( φ = 0 , θ ) defined in Eq. (4). We use these scattering matrices to encode a pair oftwo-dimensional Dirac cones living on each layer, hence θ will be tuned to π . We note that φ can be tuned to be zero because “spin-orbit coupling” is present in our construction in adifferent type of scattering matrices.The second scattering matrix S L ( θ , θ ) controls interlayer scattering. They are defined as:(See Fig. 5 for labeling): (cid:126)ψ in = ( ψ , ↑ , ψ , ↓ , ψ , ↓ , ψ , ↑ , ψ , ↓ , ψ , ↑ , ψ , ↑ , ψ , ↓ ) (cid:126)ψ out = ( ψ , ↓ , ψ , ↑ , ψ , ↑ , ψ , ↓ , ψ , ↑ , ψ , ↓ , ψ , ↓ , ψ , ↑ ) S L ( θ , θ ) (cid:126)ψ in = (cid:126)ψ out , S L ( θ , θ ) = U † y (cid:18) S T ( θ ) 00 S T ( θ ) (cid:19) U y U y = ⊗ √ (cid:18) ii (cid:19) (7)The above expression for S L ( θ , θ ) looks complicated, yet there is an intuitive way tounderstand it. To see this, observe that U y mixes ψ k +1 , ↑ and ψ k +2 , ↓ or ψ k +1 , ↓ and ψ k +2 , ↑ – it11igure 6: Illustration of how the scattering matrix S L ( θ , θ ) acts in the limit ( θ , θ ) = (0 , π , π ), ( π , , π ).mixes two components with opposite σ z eigenvalues . Identifying ψ Lk,σ y = ↑ ∼ ψ k +1 , ↑ + iψ k +2 , ↓ √ , ψ Lk,σ y = ↓ ∼ ψ k +1 , ↑ − iψ k +2 , ↓ √ ,ψ Rk,σ y = ↑ ∼ ψ k +1 , ↓ − iψ k +2 , ↑ √ , ψ Rk,σ y = ↓ ∼ ψ k +1 , ↓ + iψ k +2 , ↑ √ , (8) L/R superscript account for the fact that the above equation takes linear combinationsof two fermion modes that move toward left/right in Fig. 5(a). In this basis that accountfor the basis transformation U y , S L ( θ , θ ) is a block diagonal matrix (cid:18) S T ( θ ) 00 S T ( θ ) (cid:19) . Theidea behind the expression S L ( θ , θ ) in Eq. (7) is that half of the degrees of freedom scatteraccording to S T ( θ ) and the remaining half according to S T ( θ ), but now the degrees of freedomare diagonal in σ y instead of σ z for other scattering matrices we have encountered so far.To give further intuition, we discuss how fermion modes scatter in special limits ( θ , θ ) =(0 , π , π ), ( π , , π ). When ( θ , θ ) = (0 , θ , θ ) = ( π , π ) (the upper rightcorner in Fig. 6), the modes on the upper layer and the lower layer are completely decoupled– in this limit, there is zero interlayer scattering. ( π , , π ) limit (the lower two panels12n Fig. 6) realize more exotic scenarios in which half of the fermion modes strongly backscatterbetween different layers while the remaining half do not scatter to different layers at all. Notethat this type of interlayer coupling is the one we would like to use to construct a strong 3Dtopological insulator phase, as we saw in the earlier subsection.
Whether certain fermion modesstrongly scatter to the neighboring layers or stay in the same layers are dependent on σ y -eigenvalues and chirality of the modes.Now we are ready to construct the 3D network model that we will study in this paper. Ournetwork model have two tuning parameters, θ and θ . We would be primarily interested in theparameter space 0 ≤ θ , θ ≤ π . There are three types of scattering matrices: S L ( θ , θ ) whichwe represent as green squares graphically, S L ( π − θ , π − θ ) represented by brown squares, and S T ( π ) illustrated as gold squares. Between z = 4 j + 1 and z = 4 j + 2, scattering matrices arearranged as in Fig. 7(a). Note that here, S T ( π )’s simply unitarily relate modes at ( x, y ) = (2 i, y )to ones at ( x, y ) = (2 i + 1 , y ). Scattering matrices between z = 4 j + 3 and z = 4 j + 4, depictedin Fig. 7(b), introduce more non-trivial relations involving interlayer scattering using S L ( θ , θ )and S L ( π − θ , π − θ ). Finally, fermion modes ψ ( x,y,z =4 j +2) ,s = ↑ , ↓ and ψ ( x,y,z =4 j +3) ,s = ↑ , ↓ are relatedby random U(1) phases that encode the effect of disorder, as in Eq. (6): (cid:18) e iφ ( x,y, j +2) e − iφ ( x,y, j +2) (cid:19) (cid:18) ψ ( x,y, j +2) , ↑ ψ ( x,y, j +2) , ↓ (cid:19) = (cid:18) ψ ( x,y, j +3) , ↑ ψ ( x,y, j +3) , ↓ (cid:19) (9)Similarly, random U(1) phases e ± iφ ( x,y, j +4) relates fermion modes ψ ( x,y, j +4) , ↑ / ↓ and ψ ( x,y, j +5) , ↑ / ↓ .To get some insight about the constructed network model, it is useful to consider how ournetwork behaves in certain limits. First, let us replace both green and brown rectangles inFig. 7(b) with S L ( π , π ). Note that in this limit, interlayer coupling is completely absent, andthe network reduces to the time-reversal symmetric network model at φ QSH = 0, θ QSH = π stacked along the y -direction. This makes it clear that our model is essentially consisted oflayers, each layer hosting two 2D Dirac cones.Next, we will consider ( θ , θ ) = (0 ,
0) and ( θ , θ ) = ( π , π ) limits. In these limits, interlayerscattering matrices (green and brown rectangles in Fig. 7(b)) reduce to either of two cases inthe upper panel of Fig. 6 . In Fig. 8(a) and (b), we pictorially represent what the cross-sectionalview between z = 4 j + 3 and z = 4 j + 4 reduces to in these limits. The key insight is that inthese limits, the network reduces to bundles of insulating 1D (in z -direction) networks . Hence,near these limits, insulating phases must arise.Since ( θ , θ ) = (0 ,
0) limit does not have any surface state upon imposing open boundaryconditions, we identify a phase near this limit as a trivial insulator. The situation is a bitdifferent for the limit ( θ , θ ) = ( π , π ). In this case, one can explicitly see that the fermionmodes at y = 1 and y = 2 L y are decoupled from the bulk, and the fermion modes alongthis layer once again reduce to the time-reversal symmetric network model tuned to criticality.This surface state indicates that near ( θ , θ ) = ( π , π ), our 3D network is in a weak topologicalinsulator phase . both species of fermions backscatter within their own unit cells, yielding trivialinsulators. Near ( θ , θ ) = ( π , π ), backscattering occurs across the neighboring unit cells. Whenthere is a boundary normal to the y -direction, the layer at the boundary will be decoupled fromthe bulk, giving a surface state with two Dirac cones. Hence, near ( θ , θ ) = ( π , π ), the systemis in a weak topological phase . 13 b)(a) = Figure 7: Illustration of the cross section between (a) z = 4 j + 1 and z = 4 j + 2 (b) z = 4 j + 3and z = 4 j + 4 showing how scattering matrices are placed in our 3D network model. S T ( π )are represented by gold rectangles, S L ( θ , θ ) by green rectangles, S L ( π − θ , π − θ ) by brownrectangles. In the below panel, what it means by two rectangles connected to each otherwith solid line: incoming/outgoing modes on the right of a gold square are identified withoutgoing/incoming modes on the left of a green square.14 b)(a) Figure 8: Illustration of the cross section between z = 4 j +1 and z = 4 j +2 at (a) ( θ , θ ) = (0 , θ , θ ) = ( π , π ). In both cases, the network reduces to bundles of 1D insulating networkswhose cross sections are marked as bold rectangle. In the case of ( θ , θ ) = ( π , π ), whenopen boundary condition is imposed, the top and bottom layer decouples from the bulk andconstitute surface state of a weak topological insulator.15inally, based on the discussion in the previous subsection, near ( θ , θ ) = (0 , π ) and ( π , x and y -directions, one can show that the following two transformations map ournetwork model to itself: A : ψ ↑ / ↓ , ( x,y,z ) → ψ ↓ / ↑ , ( − x, − y,z ) , ( θ , θ ) → ( θ , θ ) B : ψ ↑ / ↓ , ( x,y,z ) → ψ ↑ / ↓ , ( x − ,y − ,z ) , ( θ , θ ) → ( π − θ , π − θ ) (10) A is simply a spatial inversion transformation followed by a parameter change and does notchange underlying topology of the network. B is more akin to the duality we saw in the 2D modelin a sense that the non-local translation can change topological features. In particular, B mapsa trivial insulator near ( θ , θ ) = (0 ,
0) to a weak topological insulator near ( θ , θ ) = ( π , π ).As a final remark on our 3D network model, our 3D model has the structure in whichscattering matrices and transfer matrices along the z -direction can be efficiently constructed.Hence, one can use recursive methods to study properties of the 3D network model, exactly asin the 2D case. In the previous section, we gave intuition behind the 3D network model and explicitly con-structed it. While we already know some aspects of the phase diagram of our model byconsidering certain limits and using the transformations in Eq. (10), in this section, we willcarry out a more in-depth study of the phase diagram using numerical methods. In particular,computational results in this section confirm that strong topological insulators do appear inthe parameter space as anticipated from the crude physical picture in the previous section. Wewill also investigate how disorder affects the phase diagram.
Here, we are interested in phase diagram of the “clean network” in which we set all the randomU(1) phases that encapsulate the effect of disorder , i.e., e iφ x,y, j +2 / j +4 in Eq. (9), to 1. Studyingclean systems is much easier than studying systems with disorder for the following two reasons:First, there is no need to take any disorder averaging. Second, due to translation symmetrypresent in the clean network, one can study systems using a Fourier representation, whichdrastically reduces computational cost. 16ll the quantities we compute in this subsection are obtained from a global scattering matrix S associated with a system of size L x × L y × L z . Each unit cell has dimension 4 × × L x × L y × L z systems to have fermion modes ψ ( x,y,z ) ,s = ↑ / ↓ in which indiceshave ranges 1 ≤ x ≤ L x , 1 ≤ y ≤ L y , and 1 ≤ z ≤ L z . S relates the fermion modes at thetwo ends z = 1 and z = 4 L z + 1. : S (cid:18) Ψ in ,z =1 Ψ in ,z =4 L z +1 (cid:19) = (cid:18) Ψ out ,z =1 Ψ out ,z =4 L z +1 (cid:19) , S = (cid:18) r t ˜ t ˜ r (cid:19) (11) r and t refer to a reflection subblock and a transmission subblock respectively.To probe bulk physics, we will use periodic boundary condition on x and y direction. It isalso useful to twist boundary conditions by U(1) numbers e ik x or e ik y . The twisted boundaryconditions are formally defined as identifying: ψ (4 L x + x,y,z ) , ↑ / ↓ ≡ e ik x ψ ( x,y,z ) , ↑ / ↓ , ψ ( x,y +2 L y ,z ) , ↑ / ↓ ≡ e ik y ψ ( x,y +2 L y ,z ) , ↑ / ↓ (12)We adopt notations S ( k x , k y ), r ( k x , k y ), and t ( k x , k y ) for the global scattering matrix and itssubblocks obtained from twisting boundary conditions.To investigate physics of surface states, we employ systems with open boundary conditionsalong the y direction. This is formally implemented by replacing interlayer scattering matrices S L ( π − θ , π − θ ) controlling backscattering between fermion modes at y = 1 and y = 2 L y (Recallthat these are brown rectangles in Fig. 7(b) lying on y = 1 and y = 2 L y ) to S L ( π , π ). Since S L ( π , π ) does not introduce any backscattering between two layers, the suggested replacementremoves any direct coupling between y = 1 modes and y = 2 L y modes, realizing systems withthe desired boundary condition. The system is still periodic along the x direction, so one maytwist boundary conditions in the x -direction.We computed the following three quantities to map out the phase diagram: • Bulk conductance across z direction of the system with size L x × L y × L z and periodicboundary conditions in both the x and y directions, defined by G bulk = Tr t † t . Thanks totranslational symmetry, one can compute G bulk equivalently from S u ( k x , k y ) and t u ( k x , k y ),the global scattering matrix and its transmission block for the 1 × × L z network withtwisted boundary condition: G bulk = L x (cid:88) n x =1 L y (cid:88) n y =1 Tr (cid:20) t u (cid:18) πn x L x , πn y L y (cid:19)(cid:21) † t u (cid:18) πn x L x , πn y L y (cid:19) . (13)Computing G bulk from t u ( k x , k y ) sidesteps computationally costly operation of directlyconstructing the large matrix S and its subblock t . • G open , conductance of the system with size L x × L y × L z , periodic boundary conditionacross x direction and open boundary condition across y direction. G open probes surfacestates for choices of parameters in which the bulk is insulating. Similar to how we com-pute G bulk , we compute G open from S ux ( k x ) and t ux ( k x ), the scattering matrix and the17ransmission subblock of the system which only differs from the original system by choice L x = 1 and twisted boundary condition along the x -direction: G open = L x (cid:88) n x =1 Tr (cid:20) t ux (cid:18) πn x L x (cid:19)(cid:21) † t ux (cid:18) πn x L x (cid:19) . (14) • Strong topological invariant Q . Here, we follow the formulation from dimensional reduc-tion previously developed in [10]. Define U T to be an unitary part of the time-reversalsymmetry action on ψ out , z = in Eq. (11). U T is analogous to U T in Eq. (5) but generalizedto act on the whole 2D arrays of modes in the cross section of our 3D network model. If k x , k y = 0 , π , due to time-reversal symmetry, the reflection block satisfies, r ( k x , k y ) U T = − U T r T ( k x , k y ) , (15)implying that r ( k x , k y ) U T is skew-symmetric. The topological invariant associated with3D strong topological insulators is proposed to be: Q = Q k x =0 Q k x = π , Q k x =0 ,π = Pf r ( k x , U T (cid:112) det r ( k x , (cid:112) det r ( k x , π )Pf r ( k x , π ) U T (16)It is shown that Q and Q π is quantized to ± r ( k x , k y ) is a unitarymatrix; this limit is valid in insulating systems with large L z in which transmission blocksvanish. The branch of the square root of the determinant should be chosen so that (cid:112) det r ( k x , k y ) should be defined continuously along the line that starts at k y = 0 andends at k y = π . Because of this, the formula requires information along the lines in themomentum space although the expression only seems to take numbers computed fromtime-reversal invariant momenta.Due to translation symmetry, computing Q from systems with any L x and L y containsame topological information, as long as L z is large enough so that the reflection block isnearly unitary. Hence, one can compute Q from the system with size 1 × × L z in thecase of the clean system.We computed these three quantities for θ ∈ (cid:2) , π (cid:3) and θ ∈ (cid:2) , π (cid:3) , and plotted the valuesin Fig. 9 (a,b,c). For computing G bulk , we chose L x = L y = L z = 256; for computing G open , wechose L x = L z = 256, L y = 32. We fixed L z = 256 for evaluating topological invariant Q .In the plot for G bulk , G bulk is very close to to zero except for two lines θ = π and θ = π .This indicates that the θ = π and θ = π lines correspond to the phase transition lines. G bulk ≈ G bulk at the transitionlines is dependent on aspect ratios of systems in general.The two critical lines divide our parameter space θ , θ ∈ [0 , π ] into four regions, eachregion representing a distinct insulating phase. We will refer to insulating regions around thecorners ( θ , θ ) = (0 , π , , π ), and ( π , π ) as region 1,2,3, and 4 respectively. Fromearlier analysis, it is clear that region 1 and region 4 correspond to a trivial phase and a weaktopological insulator respectively. Meanwhile, the phase separation suggests that the insulating18 /4 /2 (a) (b) (c) (d) TrivialInsulatorTrivialInsulator WeakTI
StrongTI StrongTI
Figure 9: Color plots for numerically computed values of (a) G bulk (b) G open (c) Q in the two-dimensional parameter space. For the plot in (c), the yellow color stands for Q = 1, while thepurple color denotes Q = −
1. (d) Phase diagram of the clean 3D network model based oninformation from the plots in (a), (b), and (c)19hases in region 2 and 3 are topologically distinct from a trivial insulator in region 1 and aweak topological phase in region 4.The plots for G open and Q make nature of the insulating phases in region 2 and 3 clearer.The insulating phase in region 1 has a vanishing G open , signaling absence of any gapless surfacestate. Meanwhile, region 4 has finite conductance G open ≈
8. We already established that thisinsulating phase is a weak topological insulator, so finite G open demonstrates the anticipatedexistence of the gapless surface state. Region 2 and region 3 also have finite surface conductancewhich is roughly half of G open in the weak topological phase. This implies emergence of singleDirac cone surface states in these regions. Computing Q confirms that these two insulatingphases are strong topological insulators . We drew the phase diagram based on our findings fromthe three quantities in Fig. 9(d).We finish our discussion of the clean network model with implications of the transformationin Eq. (10) on the strong topological insulator phases in region 2 and region 3. A maps thenetwork with a parameter value in region 2 to another one with a value in region 3. Since A involves inversion-like transformation on fermion modes, this shows that the networks in region3 are “mirror images” of the networks in region 2. Any statement about the strong topologicalinsulator phase in region 2 applies identically to the phase in region 3 as well. B has a more interesting consequence: B transforms the network in region 2 to anothernetwork in the same parameter region! Even though B involves highly non-local transformation,it does not transmute topology of the strong topological insulator phases. As mentioned, thesame statement applies to region 3 as well. Additionally, there is a special parameter line θ + θ = π invariant under the transformation B in region 2 and region 3. This provides anintriguing picture of a topological phase transition in our 3D network model distinct from the2D network models: In the 2D network models, there are dualities that connect trivial andtopological insulators, and critical points/regions between two topologically distinct phases areself-dual. In our 3D network model, B , taking a similar role to the dualities in the 2D networkmodels, relates a trivial insulator and a weak topological phase; strong topological insulatorsare intermediate phases that are mapped to themselves under B . Now, we will introduce randomness in the U(1) phases e iφ x,y, j +2 / j +4 in Eq. (9) and see howquenched randomness modifies the phase diagram. Before embarking the numerical investi-gation, it is helpful to consider general consequence of adding disorder to a topological bandinsulator to blueprint our study. Adding small amount of disorder in general fills a band gapwith electronic states. However, if disorder is sufficiently weak, these electronic states pushedto the Fermi level are localized. Hence, the system will keep its identity as a strong topologicalinsulator although notion of the band gap disappears. In the strong disorder limit, all electronicstates undergo Anderson localization, and all information about band structure is lost. Hence,we arrive at a topologically-trivial Anderson insulator.In the above consideration, the two insulators that arise in the weak disorder limit and20he strong disorder limit are topologically distinct. Hence, at intermediate disorder, thereshould be a delocalized phase, in order to avoid a contradiction. The simplest scenario wouldbe the existence of a critical metal on the intermediate disorder regime; this metallic phaseis transformed into a strong topological insulator (a trivial insulator) upon making disorderweaker (stronger). The numerical survey on tight-binding models affirm the existence of suchmetallic phases. [18, 19]Since our interest is to investigate critical phenomena associated with strong topological in-sulator phases, we would like to access both strong topological insulators near the weak disorderlimit and intermediate-disorder critical metal with our 3D network model. This motivates usto choose the random U(1) phases φ ’s chosen from a normal distribution with zero mean andvariance σ . This choice introduces an extra parameter σ that tunes the disorder strength sothat we can interpolate between the weak disorder limit and the intermediate disorder regime. The above discussion also makes it evident that our primary interest in our exploration ofthe phase diagram is whether our network models in an insulating phase or a delocalized phasegiven a choice of the parameters. In a disordered system, it is not always clear whether a systemis insulating or conducting just by evaluating quantities at a single system size. To determinethis in a more clear and rigorous fashion, we study finite-size scaling of the localization lengthin a quasi-1D geometry [12, 13].To apply this method, consider a network with size L × L × L z under condition L (cid:28) L z , withperiodic boundary conditions imposed on both the x and y directions. There is a global transfermatrix T that relates fermion modes at the both ends of the network. Eigenvalues of L z ln T † T are associated with Lyapunov exponents along z -direction. There is an algorithm using QRdecomposition to compute these eigenvalues efficiently and without numerical overflow [13].Assuming the matrix L z ln T † T does not have zero eigenvalues, its eigenvalues always come inpairs with opposite signs, + λ i and − λ i . Additionally, due to the Kramers’ theorem, eigenvaluesalways have even multiplicity.Using the computed Lyapunov exponents, one may define a localization length of this quasi-1D system ξ as: ξ = 1 λ min , λ min : The smallest positive eigenvalue of 1 L z ln T † T . (17)Then, the idea is to utilize a dimensionless parameter Λ = ξL that takes a dimension of thecross section L into account and its finite-size scaling beta function β = d ln Λ d ln L . The sign of thebeta function determines whether the system is insulating or conducting in the thermodynamiclimit. In an insulating phase, β <
0, while in a metallic phase, β > β = 0 implies scaleinvariance and is associated with a critical point. In this subsection, as a quick diagnosis ofsigns of the beta functions, we compute and compare Λ’s for L = 6 , ,
10. All quantities in thissubsection are obtained from systems with length L z = L × , averaged over four disorderrealizations. The conventional probability distribution of random phases in the Chalker-Coddington network model isa uniform distribution in which φ ’s have equal probability to take any value in [0 , π ). In our 3D networkmodel context, the uniform distribution corresponds to σ → ∞ limit, and we checked that such strong disordercompletely suppresses strong topological insulator phases in our 3D network model. Hence, uniform distributionis not an optimal distribution choice in our study of the 3D network model. .1 0.2 0.3 0.4 0.510 L=6L=8L=10 (a) L=6L=8L=10 (b) L=6L=8L=10 (c)
Figure 10: Plots of numerically computed Λ’s versus different values of disorder strength σ ,upon fixing θ = 0 and setting (a) δ = 0 .
30 (b) δ = − .
30 (c) δ = 0 . In our numerical exploration of the network model with U(1) phases generated from a nor-mal distribution, we fix θ = 0 and tune θ to be in vicinity of π . For convenience, we willoccasionally use δ = θ − π instead of θ for notational simplicity.First, to observe the disorder-tuned metal-insulator transition, we choose δ from {− . , − . , − . , . , . , . , . } . For each choice of δ , we sweep through different values of σ to inspecthow the finite-size scaling properties of Λ change. For all δ choices except for δ = 0, weconfirmed that upon increasing σ , the network model undergoes insulator-metal transitions.Plots of Λ for δ = ± . β < β > δ ’s we investigated except for δ = 0, the insulating phases in the weak disorderregime are adiabatically connected to the insulators in the network with the same ( θ , θ )and zero disorder. Therefore, as promised earlier, topology of the insulating, dirty 3D networkmodel is naturally inherited from its clean counterpart. Especially, the observed insulator-metaltransition at δ > δ = 0 in Fig. 10 (c) is more puzzling. In the clean limit, this value of theparameter corresponds to a critical point between the strong topological insulator and thetrivial insulator. However, the plot suggests that β < δ = 0. Upon increasing σ , the metallic phase with β > δ = 0 is that while there are clear regions where β > <
0, it is difficult to pinpoint a critical point between the two regions where β is supposedto be zero. At the putative critical point, β = 0 implies L -independence of Λ, but the curvesfor L = 6, L = 8, and L = 10 do not meet nicely at one point. Such “mismatch” commonlyappears in numerical analysis of Λ in disordered systems [20, 21, 22] and is attributed toirrelevant variable effect. This effect is also expected to be present in data for Fig. 10(a) and(b) upon increasing precision of data but particularly afflicts data plotted in Fig. 10(c) severely.We comment on this large irrelevant variable effect at δ = 0 near the end of this section.To reveal the nature of the insulating phase at δ = 0 with finite disorder, we fix σ =0 . , . , .
40 and scan through different values of δ ∈ [ − . , . σ = 0 . L surroundedby β < β < β < β = 0 critical point is shifted away from δ = 0. Hence, we see insulating behaviors at δ = 0 with finite disorder simply because the transition point between the strong topologicalinsulator phase and the trivial insulator is renormalized away from δ = 0. Since the insulatorat δ = 0 is adiabatically connected to the insulator that appears at δ >
0, it is identified as astrong topological insulator. We see a similar behavior for σ = 0 .
35 as well. The plot for σ = 0 .
40 in Fig. 11(b) is perplexing – it seems that β < β ≥
0. However, from our numerics of scanning through different σ ’s atfixed δ , the two insulating phases at δ = ± . β ≥ β < σ = 0 .
4. We have one guess: This unexpected behavior may come from the fact that this pointis in the vicinity of a tricritical point [25, 26, 27, 28, 29] in which the three phase transition lines– the line between the strong topological insulator and the metallic phase, the line between thetrivial insulator and the metallic phase, and the line governing the direct transition betweenthe strong topological insulator and the trivial insulator – meet. The one-parameter finite-sizescaling ansatz of the renormalized length Λ is probably not valid near this point. We believethat the large irrelevant field effect in Fig. 10(c) can be also traced back to the same reason.We present the phase diagram obtained from our numerical data in Fig. 12. Blue dots areestimated from fixing δ = θ − π to certain values and sweeping through δ and finding the valueof σ with minimum of (cid:12)(cid:12)(cid:12) ln (cid:16) Λ L =10 Λ L =8 (cid:17)(cid:12)(cid:12)(cid:12) , taken as a proxy for | β | . Orange dots are obtained similarly Here, we implicitly assumed that the phase transition from a trivial insulator to a topological insulator isdirect when disorder is weak. This assumption has been believed to be true for the following reason: Typically,a 3D Dirac cone emerges at a critical point between a strong topological insulator and a trivial insulator inband theory context. Since 3D Dirac cone is perturbatively stable against disorder, one may argue that theaforementioned direct transition is generically possible. However, it was recently proposed that non-perturbativerare-region effect actually destabilizes Dirac semimetal [23, 24]. This implies that any strength of disorder turnsthe Dirac semimetal critical point into a diffusive metallic region, forbidding any direct transition between astrong topological insulator and a trivial insulator. This issue is not focus of our work, and this rare region effectis likely to be inaccessible from the small system numerics in this paper. Hence, we will presume that there is adirect transition between the strong topological phase insulator and the trivial insulator phase in our networkmodel at weak disorder throughout this paper. However, we do not rule out possibility that the “critical point”governing the phase transition between the two insulating phases actually belongs to a metallic region. a) (b) L=6L=8L=10 L=6L=8L=10
Figure 11: Plots of numerically computed Λ versus different values of δ , upon fixing θ = 0and setting (a) σ = 0 .
30 (b) σ = 0 . TrivialInsulator Strong TIMetal?
Figure 12: Estimated phase boundaries of our 3D network model upon fixing θ = 0 and varying( δ = θ − π , σ ). 24y fixing σ and scanning through δ . Near the area where the metal-insulator boundary andthe direct strong TI-trivial insulator boundary should meet, as we stated earlier, estimating thephase boundaries is difficult, and their location is left as blank in a question-marked region inFig. 12. Here, we estimate a critical exponent for the strong topological insulator-metal transition usingthe finite-size scaling properties of the renormalized length Λ. Here, we follow finite-size scalingframework developed in [20, 22] incorporating corrections from irrelevant variables. In thissection, we fix δ = 0 . θ = 0 and study strong topological insulator to metal transitiontuned by σ in this work.If we follow one-parameter scaling of Λ near a metal-insulator transition reviewed in theearlier section, at the critical point σ c , Λ is independent of the system size. One can alsoformulate a simple ansatz of how Λ behaves near critical point with a critical exponent ν associated with scaling of the correlation length ξ ∼ | σ − σ c | ν . Λ, or equivalently log Λ, is adimensionless number and hence we assume it takes a one-parameter scaling form,log Λ( σ, L ) = f (cid:0) ( σ − σ c ) L /ν (cid:1) (18)In addition, it is natural to expect that Λ( σ, L ) is an analytic function of σ for any finite L .Hence, one may Taylor-expand log Λ( σ, L ) near σ c as:log Λ( σ, L ) = log Λ c + A ( σ − σ c ) L /ν + A ( σ − σ c ) L /ν + · · · (19)where Λ c is defined as Λ at the criticality.In numerical practices, σ is a UV variable which includes potentially infinite number ofirrelevant variables in IR. While the effect of irrelevant variables should vanish in the ther-modynamic limit, in numerical calculation of Λ, due to the small system sizes, this effect willinadvertently show up in data, most commonly in a manner that points with ∂ Λ( σ,L ) ∂L = 0 seemto shift systematically to one direction as L increases; recall that this effect appeared in ourearlier studies of the phase diagram of our 3D network model. To accommodate this effect, weintroduce y , the smallest scaling dimension of irrelevant scaling variables, and some functionof σ , g ( σ ), which represent the dominant irrelevant scaling variable. Then, one can modify theansatz in Eq. (18) as: log Λ( σ, L ) = f (cid:18) ( σ − σ c ) L /ν , g ( σ ) L y (cid:19) (20)Starting from the above modified ansatz, once again one may assume analyticity and Taylor-expand, g first then f . In this section, we assume g ( σ ) is a constant and will fit our numericaldata to the following formula:log Λ( σ, L ) ≈ log Λ c + A ( σ − σ c ) L /ν + A ( σ − σ c ) L /ν + L − y (cid:2) A + A ( σ − σ c ) L /ν + A ( σ − σ c ) L /ν (cid:3) (21)25 .50 0.51 0.52 0.530.500.751.001.251.501.75 l n L=4L=6L=8L=10L=12
Figure 13: Plots of numerically computed Λ we used for extracting the critical exponent ν .For the non-linear fitting, we numerically computed Λ from the network with size L × L × L z for L = 4 , , , , L z = 10000 × L for L (cid:54) = 12 and L z = 8000 × L for L = 12, taking disorderaveraging for 20 independent realization. The estimated uncertainty in data points range from0 . − . L , total of 125 data points. The plot clearly shows theirrelevant variable effect in L = 4 and L = 6 data.We could obtain a fit with reasonable stability, checked by varying initial condition of thenon-linear fit and Monte-Carlo resampling. The fitting result yields: (95% confidence intervalon error estimate): ν = 1 . ± . y = 3 . ± . σ c = 0 . ± . c = 0 . ± . χ = 109 . ν is notcompatible with symplectic Anderson transition exponent computed from (non-topological)tight-binding models in [30], at least in 95% confidence interval. This is also at odds with theconjecture in [18] that the strong topological insulator-metal transition belongs to the sameuniversality class to the plain-vanilla symplectic Anderson transition. While it is tempting toclaim a new universality class based on our result, we believe that our data and fitting arerather crude. Hence, we advise readers to take results in this section rather as demonstrationthat it is feasible to extract critical exponents from the finite-size scaling analysis using ournetwork model than as a quantitative prediction. More detailed numerical studies involving26ata with smaller error bars and evaluations of critical exponents along different lines in theparameter space are needed to make any definite conclusion about the universality classes. In this paper, we constructed the 3D network model that exhibits both weak and strong time-reversal symmetric topological insulator phases. Our model, similar to the original 2D networkmodel by Chalker and Coddington, is amenable to numerical methods using quasi-1D geome-tries; this enables extensive numerical survey of our model. Through numerics, we observedthat some insulating phases, especially the strong topological phases, in the 3D network modeltransform into metallic phases upon adding sufficiently strong quenched randomness. Thisresult is consistent with previous results based on tight-binding Hamiltonians.Our model also has novel features that were not previously appreciated from the analysisbased on tight binding models: First, we found that in our network model, there is a non-local transformation that relates a weak topological insulator phase and a trivial phase. Strongtopological insulators are intermediate phases between the two aforementioned phases thatare mapped to itself under the transformation. Second, we showed that through finite-sizescaling analysis one can extract the critical exponent of the strong topological insulator-metaltransition. In the previous approach using tight-binding models technical difficulties associatedwith band gap scales complicated such analysis.There are several future directions. The most immediate extension of our work will beimproving our estimates on the critical exponent for the strong topological insulator-metaltransition by enhancing precision and computing exponents for different choices of parameters.This will allow one to make definite conclusions about whether the topological insulator-metaltransition is in the same universality class as the sympletic Anderson transition without anytopological feature.It would be also interesting to extend our construction to other topological phase transitionsin three spatial dimensions. For instance, the field-theory description of time-reversal invari-ant topological superconductor is largely similar to that of time-reversal invariant topologicalinsulators, except the complex fermions in the latter are replaced by real Majorana fermions.Hence, we expect that there is a fairly simple generalization of our model to the 3D topologicalsuperconductor transition. In contrast, a 3D network model for AIII topological insulator tran-sitions poses a more interesting challenge since there is no intrinsic 2D AIII topological insulatorand hence no weak 3D AIII topological insulator. Our approach inadvertently includes weaktopological insulator phases in the phase diagram. This leads us to believe that designing theAIII 3D network model involves more non-trivial modification of our methodology.Finally, we found the non-local transformation under which the strong topological insula-tor phases are invariant, but we did not comment on physical origin of this transformation.Connecting this non-local transformation to physical dualities/symmetries would be helpful inunderstanding topological phase transitions in greater depth.27 a) (b)
Figure 14: Graphical depiction of notations used in (a) formulas involving a single scatter-ing/transfer matrix and (b) formulas involving two scattering/transfer matrices .
Acknowledgement
We thank Jason Alicea, Jing-Yuan Chen, Andrew Essin, Prashant Kumar, Gil Refael, andXiao-Qi Sun for discussion. This work was supported in part by the US Department of Energy,Office of Basic Energy Sciences, Division of Materials Sciences and Engineering, under contractnumber DEAC02-76SF00515
Appendix A More on Scattering Matrices and TransferMatrices
In this appendix, we present formulas involving scattering matrices and transfer matrices usedfor the numerical implementation. We also discuss how these formulas can be used for con-structing a global scattering matrix in practice.
A.1 Conversion formula between scattering matrices and transfermatrices
In the main text, when reviewing the 2D network model and constructing our 3D networkmodel, we formulated the models using scattering matrices. Meanwhile, in numerical practices,one commonly uses transfer matrices. Here, we give formula on how scattering matrices and28ransfer matrices are mathematically related. Generally a single scattering matrix and a transfermatrix take the following form: S (cid:18) Ψ , + Ψ , − (cid:19) = (cid:18) Ψ , − Ψ , + (cid:19) , S = (cid:18) A BC D (cid:19) T (cid:18) Ψ , + Ψ , − (cid:19) = (cid:18) Ψ , + Ψ , − (cid:19) , T = (cid:18) a bc d (cid:19) (23) Ψ , + , Ψ , − , Ψ , + , and Ψ , − are length- n vectors, and each block of S and T matrix is a n × n matrix. This setup is illustrated in Fig. 14(a). Each block of a scattering matrix can beexpressed in terms of transfer matrix subblocks, and vice versa, as the following: S = (cid:18) − d − c d − a − bd − c bd − (cid:19) , T = (cid:18) C − DB − A DB − − B − A B − (cid:19) (24)Above two equations allow one to convert a scattering matrix to a transfer matrix and viceversa. A.2 Formula involving two scattering matrices
Now, we will consider a setup illustrated in Fig. 14(b) in which there are two scattering nodesdescribed by scattering matrices S and S , or equivalently two transfer matrices T and T .We label fermion modes and matrix subblocks analogously to Eq. (23), as the following: S (cid:18) Ψ , + Ψ , − (cid:19) = (cid:18) Ψ , − Ψ , + (cid:19) , S = (cid:18) A B C D (cid:19) S (cid:18) Ψ , + Ψ , − (cid:19) = (cid:18) Ψ , − Ψ , + (cid:19) , S = (cid:18) A B C D (cid:19) (25)Our goal is to construct a formula for a scattering matrix/transfer matrix of the combinedsystem S c / T c that gives relation between Ψ , ± and Ψ , ± : S c (cid:18) Ψ , + Ψ , − (cid:19) = (cid:18) Ψ , − Ψ , + (cid:19) , T c (cid:18) Ψ , + Ψ , − (cid:19) = (cid:18) Ψ , + Ψ , − (cid:19) (26)Building T c is easy – due to multiplicative structure of transfer matrices, T c is simply equalto T T . Meanwhile, S c can be built from subblocks of S and S according to the followingformula: S c = (cid:18) A + B ( I − A D ) − A C B ( I − A D ) − B C C + C D ( I − A D ) − A C D + C D ( I − A D ) − B (cid:19) (27)Note that above formula requires one to take matrix inverse, an operation more costly thanmatrix multiplication. Hence, constructing T c with matrix multiplication is generally fasterthan building S c from the above formula. As an alternative method, if one knows T − (note the29nverse) and S , S c may be computed from subblocks of T − and S according to the followingformula: T − = (cid:18) a b c d (cid:19) , S c = (cid:18) A c B c C c D c (cid:19) C c = ( a − D c ) − C D c = ( a − D c ) − ( D d − b ) A c = A + B c C c B c = B c D c + B d (28)The above expression also requires matrix inversion. However, it has a smaller number ofmatrix multiplications than Eq. (27) and may be evaluated more efficiently. A.3 Constructing a global scattering matrix in practice
Based on the result from the previous two subsections, we now consider a setting in which wehave fermion modes Ψ i, ± , i = 1 , , · · · n + 1 living on i th slice of the system and scatteringmatrices S j ( j = 1 , , · · · n ) that endow unitary relations between fermion modes Ψ j, ± and Ψ j +1 , ± . We will show how to construct a global scattering matrix S that relates Ψ , ± and Ψ n +1 , ± , the fermion modes that reside on the both ends of the system in numerical practice.We assume that cost of converting S j to T j or T − j is negligible; this is valid in the networkmodel context since scattering matrices across two neighboring slices are often composed ofsmaller blocks due to locality. Hence, one may convert scattering matrices and transfer matricesto each other by applying the formula Eq. (24) to the individual blocks. Numerical cost ofsuch operations is suppressed compared to matrix inversion or matrix multiplication operationsapplied to the full-sized transfer/scattering matrices when system size is large.If one considers computational cost only, the most efficient method is to compute a globaltransfer matrix T = T n T n − · · · T and convert T to S . However, this method does not workin practice, due to the following reason: T , especially when the system length n is large, oftenhas very large eigenvalues that lead to numerical overflow. Alternatively, one may start from S and iteratively apply Eq. (27) to S , S , · · · , S n or Eq. (28) to T − , T − , · · · , T − n . While thismethod is free from overflow, this method is clearly slower than multiplying transfer matrices.In numerical implementation, it is most practical to use the hybrid method – use matrixmultiplications to construct (inverse) transfer matrices across several slices, number of slicessmall enough to avoid overflow, and apply Eq. (28) iteratively to construct a global scatteringmatrix S from the transfer matrices for multi-slice chunks. References [1] J. T. Chalker and P. D. Coddington, “Percolation, quantum tunnelling and the integerhall effect,”
Journal of Physics C: Solid State Physics no. 14, (May, 1988) 2665–2679.302] A. Altland and M. R. Zirnbauer, “Nonstandard symmetry classes in mesoscopicnormal-superconducting hybrid structures,” Phys. Rev. B (Jan, 1997) 1142–1161.[3] S. Ryu, A. P. Schnyder, A. Furusaki, and A. W. Ludwig, “Topological insulators andsuperconductors: tenfold way and dimensional hierarchy,” New Journal of Physics no. 6, (2010) 065010.[4] S. Cho and M. P. A. Fisher, “Criticality in the two-dimensional random-bond isingmodel,” Phys. Rev. B (Jan, 1997) 1025–1031.[5] I. A. Gruzberg, A. W. W. Ludwig, and N. Read, “Exact exponents for the spin quantumhall transition,” Phys. Rev. Lett. (May, 1999) 4524–4527.[6] E. J. Beamond, J. Cardy, and J. T. Chalker, “Quantum and classical localization, thespin quantum hall effect, and generalizations,” Phys. Rev. B (May, 2002) 214301.[7] H. Obuse, A. Furusaki, S. Ryu, and C. Mudry, “Two-dimensional spin-filtered chiralnetwork model for the Z quantum spin-hall effect,” Phys. Rev. B (Aug, 2007) 075301.[8] I. C. Fulga, A. R. Akhmerov, J. Tworzyd(cid:32)lo, B. B´eri, and C. W. J. Beenakker, “Thermalmetal-insulator transition in a helical topological superconductor,” Phys. Rev. B (Aug, 2012) 054505.[9] J. T. Chalker and A. Dohmen, “Three-dimensional disordered conductors in a strongmagnetic field: Surface states and quantum hall plateaus,” Phys. Rev. Lett. (Dec,1995) 4496–4499.[10] I. C. Fulga, F. Hassler, and A. R. Akhmerov, “Scattering theory of topological insulatorsand superconductors,” Phys. Rev. B (Apr, 2012) 165409.[11] E. Abrahams, P. W. Anderson, D. C. Licciardello, and T. V. Ramakrishnan, “Scalingtheory of localization: Absence of quantum diffusion in two dimensions,” Phys. Rev. Lett. (Mar, 1979) 673–676.[12] A. MacKinnon and B. Kramer, “One-parameter scaling of localization length andconductance in disordered systems,” Phys. Rev. Lett. (Nov, 1981) 1546–1549.[13] A. MacKinnon and B. Kramer, “The scaling theory of electrons in disordered solids:additional numerical results,” Zeitschrift f¨ur Physik B Condensed Matter no. 1, (1983)1–13.[14] X.-Q. Sun, M. Xiao, T. c. v. Bzduˇsek, S.-C. Zhang, and S. Fan, “Three-dimensionalchiral lattice fermion in floquet systems,” Phys. Rev. Lett. (Nov, 2018) 196401.[15] S. Higashikawa, M. Nakagawa, and M. Ueda, “Floquet chiral magnetic effect,”
Phys. Rev.Lett. (Aug, 2019) 066403.[16] E. Sagi and Y. Oreg, “From an array of quantum wires to three-dimensional fractionaltopological insulators,”
Phys. Rev. B (Nov, 2015) 195137.3117] X.-L. Qi, T. L. Hughes, and S.-C. Zhang, “Topological invariants for the fermi surface ofa time-reversal-invariant superconductor,” Phys. Rev. B (Apr, 2010) 134508.[18] S. Ryu and K. Nomura, “Disorder-induced quantum phase transitions inthree-dimensional topological insulators and superconductors,” Phys. Rev. B (Apr,2012) 155138.[19] K. Kobayashi, T. Ohtsuki, and K.-I. Imura, “Disordered weak and strong topologicalinsulators,” Phys. Rev. Lett. (Jun, 2013) 236803.[20] A. MacKinnon, “Critical exponents for the metal-insulator transition,”
Journal ofPhysics: Condensed Matter no. 13, (Mar, 1994) 2511–2518.[21] Z. Wang, B. c. v. Jovanovi´c, and D.-H. Lee, “Critical conductance and its fluctuations atinteger hall plateau transitions,” Phys. Rev. Lett. (Nov, 1996) 4426–4429.[22] K. Slevin and T. Ohtsuki, “Corrections to scaling at the anderson transition,” Phys. Rev.Lett. (Jan, 1999) 382–385.[23] R. Nandkishore, D. A. Huse, and S. L. Sondhi, “Rare region effects dominate weaklydisordered three-dimensional dirac points,” Phys. Rev. B (Jun, 2014) 245110.[24] J. H. Pixley, D. A. Huse, and S. Das Sarma, “Rare-region-induced avoided quantumcriticality in disordered three-dimensional dirac and weyl semimetals,” Phys. Rev. X (Jun, 2016) 021042.[25] P. Goswami and S. Chakravarty, “Quantum criticality between topological and bandinsulators in 3 + 1 dimensions,” Phys. Rev. Lett. (Nov, 2011) 196803.[26] K. Kobayashi, T. Ohtsuki, K.-I. Imura, and I. F. Herbut, “Density of states scaling at thesemimetal to metal transition in three dimensional topological insulators,”
Phys. Rev.Lett. (Jan, 2014) 016402.[27] S. V. Syzranov, V. Gurarie, and L. Radzihovsky, “Unconventional localization transitionin high dimensions,”
Phys. Rev. B (Jan, 2015) 035133.[28] B. Roy and S. Das Sarma, “Diffusive quantum criticality in three-dimensional disordereddirac semimetals,” Phys. Rev. B (Dec, 2014) 241112.[29] J. H. Pixley, P. Goswami, and S. Das Sarma, “Disorder-driven itinerant quantumcriticality of three-dimensional massless dirac fermions,” Phys. Rev. B (Feb, 2016)085103.[30] Y. Asada, K. Slevin, and T. Ohtsuki, “Anderson transition in the three dimensionalsymplectic universality class,” Journal of the Physical Society of Japan74