The Z_2 network model for the quantum spin Hall effect: two-dimensional Dirac fermions, topological quantum numbers, and corner multifractality

The quantum spin Hall effect shares many similarities (and some important differences) with the quantum Hall effect for the electric charge. As with the quantum (electric charge) Hall effect, there exists a correspondence between bulk and boundary physics that allows to characterize the quantum spin Hall effect in diverse and complementary ways. In this paper, we derive from the network model that encodes the quantum spin Hall effect, the so-called Z_2 network model, a Dirac Hamiltonian in two dimensions. In the clean limit of this Dirac Hamiltonian, we show that the bulk Kane-Mele Z_2 invariant is nothing but the SU(2) Wilson loop constructed from the SU(2) Berry connection of the occupied Dirac-Bloch single-particle states. In the presence of disorder, the non-linear sigma model (NLSM) that is derived from this Dirac Hamiltonian describes a metal-insulator transition in the standard two-dimensional symplectic universality class. In particular, we show that the fermion doubling prevents the presence of a topological term in the NLSM that would change the universality class of the ordinary two-dimensional symplectic metal-insulator transition. This analytical result is fully consistent with our previous numerical studies of the bulk critical exponents at the metal-insulator transition encoded by the Z_2 network model. Finally, we improve the quality and extend the numerical study of boundary multifractality in the Z_2 topological insulator. We show that the hypothesis of two-dimensional conformal invariance at the metal-insulator transition is verified within the accuracy of our numerical results.


Introduction
Spin-orbit coupling has long been known to be essential to account for the band structure of semiconductors, say, semiconductors with the zink-blende crystalline structure. Monographs have been dedicated to reviewing the effects of the spin-orbit coupling on the Bloch bands of conductors and semiconductors [1]. Electronic transport properties of metals and semiconductors in which impurities are coupled to the conduction electrons by the spin-orbit coupling, i.e., when the impurities preserve the time-reversal symmetry but break the spin-rotation symmetry, are also well understood since the prediction of weak antilocalization effects [2]. Hence, the prediction of the quantum spin Hall effect in two-dimensional semiconductors with time-reversal symmetry but a sufficiently strong breaking of spin-rotation symmetry is rather remarkable in view of the maturity of the field dedicated to the physics of semiconductors [3,4,5,6]. The quantum spin Hall effect was observed in HgTe/(Hg,Cd)Te quantum wells two years later [7]. Even more remarkably, this rapid progress was followed by the prediction of threedimensional topological insulators [8,9,10] and its experimental confirmation for Bibased compounds [11,12,13,14,15].
The quantum spin Hall effect, like its relative, the quantum (electric charge) Hall effect, can be understood either as a property of the two-dimensional bulk or as a property of the one-dimensional boundary. The bulk can be characterized by certain integrals over the Brillouin zone of Berry connections calculated from Bloch eigenstates. These integrals are only allowed to take discrete values and are examples of topological invariants from the mathematical literature. As is well known, the topological number ν takes integer values for the quantum (electric charge) Hall effect [16]. By contrast, it takes only two distinct values (ν = 0 or 1) for time-reversal invariant, Z 2 topological band insulators [4,8,9,10,17]. Because they are quantized, they cannot change under a small continuous deformation of the Hamiltonian, including a perturbation that breaks translation invariance, i.e., disorder.
The bulk topological quantum numbers are closely connected with the existence of stable gapless edge states along the boundary of a topological insulator, or more precisely along the interface between two insulators with different topological numbers. The number of gapless edge modes is determined by the difference of the topological numbers. On the edge of a two-dimensional Z 2 topological band insulator with ν = 1, there exists helical edge states, a Kramers' pair of counter propagating modes, which interpolates between the bulk valence band and the bulk conduction band. If one changes the Fermi energy from the center of the band gap to lower energies through the conduction band, one should observe a transition from a Z 2 topological insulator to a metal, and then from a metal to a trivial band insulator (ν = 0) without helical edge states. Since both helical edge states and a metallic phase are stable against (weak) disorder (due to the quantized topological number and to weak anti-localization, respectively), the same sequence of phases should appear as the Fermi energy is varied even in the presence of disorder, as confirmed recently by numerical simulations [18,19]. A question one can naturally ask is then whether there is any difference between the critical phenomena at the metal-to-Z 2 -topological-insulator transition and those at the metal-to-trivial-insulator transition. This is the question which we revisit in this paper, extending our previous studies [19,20]. It will become clear that one needs to distinguish between bulk and boundary properties in the universal critical phenomena.
For the quantum (electric charge) Hall effect, the Chalker-Coddington network model serves as a standard model for studying critical properties at Anderson transition between different quantum Hall states [21]. The elementary object in the Chalker-Coddington network model is chiral edge states. These edge states are plane waves propagating along the links of each plaquette which represents a puddle of a quantum Hall droplet formed in the presence of spatially slowly varying potential. They are chiral as they represent the mode propagating along equipotential lines in the direction determined by the external magnetic field. The Chalker-Coddington network model is a unitary scattering matrix that scales in size with the number of links defining the network, and with a deterministic parameter that quantifies the relative probability for an incoming mode to scatter into a link rotated by +π/2 or −π/2. By tuning this parameter through the value 1/2, one can go through a transition from one insulating phase to another insulating phase, with the topological number ν changed by one. This remains true even when the phase of an edge state along any link is taken to be an independent random number to mimic the effects of static local disorder. The Chalker-Coddington model is a powerful tool to characterize the effects of static disorder on the direct transition between two successive integer quantum Hall states. It has demonstrated that this transition is continuous and several critical exponents at this transition have been measured from the Chalker-Coddington model [21,22].
The present authors have constructed in [19] a generalization of the Chalker-Coddington model that describes the physics of the two-dimensional quantum spin Hall effect. We shall call this network model the Z 2 network model, which will be briefly reviewed in section 2. As with the Chalker-Coddington model, edge states propagate along the links of each plaquette of the square lattice. Unlike the Chalker-Coddington model there are two edge states per link that form a single Kramers' doublet, which corresponds to helical edge states moving along a puddle of a quantum spin Hall droplet. Kramers' doublets undergo the most general unitary scattering compatible with timereversal symmetry at the nodes of the square lattice. The Z 2 network model is thus a unitary scattering matrix that scales in size with the number of links defining the network and that preserves time-reversal symmetry. The Z 2 network model supports one metallic phase and two insulating phases, as we discussed earlier ‡. The metallic phase prevents any direct transition between the insulating phases and the continuous phase transition between the metallic and any of the insulating phases belongs to the two-dimensional symplectic universality class of Anderson localization [2].
Numerical simulations have shown that bulk properties at metal-insulator transition in the Z 2 network model are the same as those at conventional metal-insulator transitions in the two-dimensional symplectic symmetry class [19,20]. In fact, one can understand this result from the following general argument based on universality. The non-linear sigma model (NLSM) description is a very powerful, standard theoretical approach to Anderson metal-insulator transition [23]. A NLSM can have a topological term if the homotopy group of the target manifold, which is determined by the symmetry of the system at hand, is nontrivial. Interestingly, in the case of the symplectic symmetry class, as is called the statistical ensemble of systems (including quantum spin Hall systems) that are invariant under time reversal but are not invariant under SU(2) spin rotation, the NLSM admits a Z 2 topological term [24,25,26]. Moreover, the NLSM in the symplectic symmetry class with a Z 2 topological term cannot support an insulating phase. This can be seen from the fact that this NLSM describes surface Dirac fermions of a three-dimensional Z 2 topological insulator which are topologically protected from Anderson localization [27,28,29]. This in turn implies that any two-dimensional metalinsulator transition in time-reversal-invariant but spin-rotation-noninvariant systems should be in the same and unique universality class that is encoded by the NLSM without a topological term in the (ordinary) symplectic class.
Whereas bulk critical properties at the transition between a metal and a Z 2 topological insulator do not depend on the topological nature of the insulating phase, there are boundary properties that can distinguish between a topologically trivial and non-trivial insulating phases. Boundary multifractality is a very convenient tool to probe any discrepancy between universal bulk and boundary properties at Anderson transition [30,31]. To probe this difference, the present authors performed a multifractal analysis of the edge states that propagate from one end to the other in a network model at criticality with open boundary condition in the transverse direction [20]. It was found that boundary multifractal exponents are sensitive to the presence or absence of a helical ‡ The presence or absence of a single helical edge state in an insulating phase is solely dependent on the boundary conditions which one imposes on the network model.
Kramers' doublet propagating along the boundary.
The goal of this paper is 2-fold: (i) to establish a direct connection between the Z 2 network model and a Hamiltonian description of the Z 2 topological insulator perturbed by time-reversal symmetric local static disorder.
(ii) to improve the quality and extend the numerical study of boundary multifractality in the Z 2 topological insulator.
For item (i), in section 3, we are going to relate the Z 2 network model to a problem of Anderson localization in the two-dimensional symplectic universality class that is encoded by a stationary 4 × 4 Dirac Hamiltonian perturbed by static disorder that preserves time-reversal symmetry but breaks spin-rotation symmetry. This result is a natural generalization of the fact that the Chalker-Coddington network model can be related [32] to a 2 × 2 Dirac Hamiltonian with static disorder [33]. In the clean limit, we shall characterize the Z 2 insulating phases in the 4 × 4 Dirac Hamiltonian by a Z 2 topological invariant. In particular, we show that an SU(2) Wilson loop of Berry connection of Bloch wave functions is equivalent to the Z 2 index introduced by Kane and Mele [4]. The 4 × 4 Dirac Hamiltonian will allow us to make contact between the Z 2 network model and the NLSM description of two-dimensional Anderson localization in the symplectic universality class derived 30 years ago by Hikami et al. in [2]. In our opinion, this should remove any lingering doubts that the metal-insulator transition between a two-dimensional metallic state and a two-dimensional Z 2 insulator that is driven by static disorder is anything but conventional.
For item (ii), besides improving the accuracy of the critical exponents for onedimensional boundary multifractality in the Z 2 network model, we compute critical exponents for two zero-dimensional boundaries (corners) in section 4. We shall use these critical exponents to verify the hypothesis that conformal invariance holds at the metal-insulator transition and imposes relations between lower-dimensional boundary critical exponents.

Definition of the Z 2 network model for the quantum spin Hall effect
The Z 2 network model is defined as follows. First, one draws a set of corner sharing square plaquettes on the two-dimensional Cartesian plane. Each edge of a plaquette is assigned two opposite directed links. This is the network. There are two types S and S ′ of shared corners, which we shall call the nodes of the network. Second, we assign to each directed link an amplitude ψ, i.e., a complex number ψ ∈ C. Any amplitude ψ is either an incoming or outgoing plane wave that undergoes a unitary scattering process at a node. We also assign a 4 × 4 unitary matrix S to each node of the network. The set of all directed links obeying the condition that they are either the incoming or outgoing plane waves of the set of all nodal unitary scattering matrices defines a solution to the Z 2 network model. To construct an explicit representation of the Z 2 network model, the center of each plaquette is assigned the coordinate (x, y) with x and y taking integer values, as is done in figure 1. We then label the 8 directed links ψ nσ (x, y) of any given plaquette by the coordinate (x, y) of the plaquette, the side n = 1, 2, 3, 4 of the plaquette with the convention shown in figure 1, and the spin index σ =↑ or σ =↓ if the link is directed counterclockwise or clockwise, respectively, relative to the center of the plaquette. The 4 × 4 unitary S-matrix is then given by  at any node of type S or as  at any node of type S ′ , with Here, the 4 × 4 unitary matrix is presented with the help of the unit 2 × 2 matrix s 0 and of the 2 × 2 matrix (s 1 , s 2 , and s 3 are the 2 × 2 Pauli matrices) that are both acting on the spin indices σ =↑, ↓, together with the real-valued parameters r := tanh X, t := 1 cosh X , with For later use, we shall also introduce the real-valued parameter β ∈ [0, π] through r = cos β, t = sin β.
The Z 2 network model is uniquely defined from the scattering matrices S and S ′ . By construction, the S-matrix is time-reversal symmetric, i.e., and a similar relation holds for S ′ .
In [19], we obtained the phase diagram of the Z 2 network model shown schematically in figure 2(a). Thereto, (X, θ) are spatially uniform deterministic parameters that can be changed continuously. On the other hand, the phases χ n of all link plane waves in the Z 2 network model are taken to be independently and uniformly distributed random variables over the range [0, 2π). The line θ = 0 is special in that the Z 2 network model reduces to two decoupled Chalker-Coddington network models [19]. Along the line θ = 0, the point realizes a quantum critical point that separates two insulating phases differing by one gapless edge state or, equivalently, by one unit in the Hall conductivity, per spin. Alternatively, θ can also be chosen to be randomly and independently distributed at each node with the probability sin(2θ) over the range (0, π/2). This leaves X as the sole deterministic parameter that controls the phase diagram as shown in figure 2(b). When performing numerically a scaling analysis with the size of the Z 2 network model, one must account for the deviations away from one-parameter scaling induced by irrelevant operators. The Z 2 network model with a randomly distributed θ minimizes such finitesize effects (see [19]). The metallic phase is surrounded by the two insulating phases with the critical points X s and X l (> X s )) for 0 < θ < π/2. The fixed point denoted by a filled (green) square along the boundary θ = 0 is the unstable quantum critical point located at X CC = ln(1 + √ 2) separating two insulating phases in the Chalker-Coddington model. The fixed point denoted by the filled (blue) rhombus at the upper left corner is the unstable metallic phase. The shape of the metallic phase is controlled by the symmetry crossover between the unitary and symplectic symmetry classes. (b) The phase diagram for Z 2 network model with randomly distributed θ over the range (0, π/2).

Two-dimensional Dirac Hamiltonian from the Z 2 network model
The Chalker-Coddington model is related to the two-dimensional Dirac Hamiltonian as was shown by Ho and Chalker in [32]. We are going to establish the counterpart of this connection for the Z 2 network model. A unitary matrix is the exponential of a Hermitian matrix. Hence, our strategy to construct a Hamiltonian from the Z 2 network model is going to be to view the unitary scattering matrix of the Z 2 network model as a unitary time evolution whose infinitesimal generator is the seeked Hamiltonian. To this end, we proceed in two steps in order to present the Z 2 network model into a form in which it is readily interpreted as a unitary time evolution. First, we change the choice of the basis for the scattering states and select the proper unit of time. We then perform a continuum approximation, by which the Z 2 network model is linearized, so to say. This will yield an irreducible 4-dimensional representation of the Dirac Hamiltonian in (2 + 1)-dimensional space and time, a signature of the fermion doubling when deriving a continuum Dirac Hamiltonian from a time-reversal symmetric and local two-dimensional lattice model.

Change of the basis for the scattering states and one-step time evolution
Our goal is to reformulate the Z 2 network model defined in Sec. 2 in such a way that the scattering matrix maps incoming states into outgoing states sharing the same internal and space labels but a different "time" label. This involves a change of basis for the scattering states and an "enlargement" of the Hilbert space spanned by the scattering states. The parameter θ is assumed to be spatially uniform. We choose the plaquette (x, y) of the network.
At node S of the plaquette (x, y), we make the basis transformation and write the where we have defined and Here given n = 1, 2, 3, 4 and σ =↑, ↓, we have introduced the shift operators acting on ψ nσ (x, y), and similarly on the phases χ n (x, y) ∈ [0, 2π). We note that the scattering matrix N S is multiplied by the unitary matrix U from the left and the right in (12), because the Kramers' doublet acquires exactly the same phase χ n when traversing on the edge n of the plaquette (x, y) before and after experiencing the scattering N S at the node S. At node S ′ of the plaquette (x, y), we make the basis transformation and rewrite the scattering matrix S ′ (2) into the form  where we have defined As it should be Next, we introduce the discrete time variable l ∈ Z as follows. We define the elementary discrete unitary time evolution to be  Here, to treat on equal footing the nodes of type S and S ′ , we have enlarged the scattering basis with the introduction of the doublets Due to the off-diagonal block structure in the elementary time evolution, it is more convenient to consider the "one-step" time evolution operator defined by  The two Hamiltonians generating this unitary time evolution are then Evidently, the additivity of the logarithm of a product implies that From now on, we will consider H SS ′ exclusively since M S ′ S = exp(iH S ′ S ) merely duplicates the information contains in M SS ′ = exp(iH SS ′ ).

Dirac Hamiltonian close to θ = 0
In this section, we are going to extract from the unitary time-evolution (21) To this end and following [32], it is convenient to measure the link phases χ n (n = 1, 2, 3, 4) relative to their values when they carry a flux of π per plaquette. Hence, we redefine on all plaquettes. Our strategy consists in performing an expansion of defined in (22) to leading order in powers of with n = 1, 2, 3, 4 where ∂ x,y is the generator of infinitesimal translation on the network (the two-dimensional momentum operator). When θ = 0, the unitary time-evolution operator at the plaquette (x, y) is given by whereby with the 2 × 2 operator-valued matrices Observe that in the limit θ = 0, the Z 2 network model reduces to two decoupled U(1) network models where each time evolution is essentially the same as the one for the U(1) network model derived in [32].
In the vicinity of the Chalker-Coddington quantum critical point (24), we find the 4 × 4 block diagonal Hamiltonian where the 2 × 2 block are expressed in terms of linear combinations of the 2 × 2 unit matrix σ 0 and of the Pauli matrices σ x , σ y , and σ z according to and Thus, each 2×2 block Hamiltonian is of the Dirac form whereby the linear combinations enter as a scalar gauge potential and a vector gauge potential would do, respectively.
Any deviation of θ from θ = 0 lifts the reducibility of (36). To leading order in θ and close to the Chalker-Coddington quantum critical point (24), where we have set m = χ n = 0 and t x,y ± = 1. We obtain to this order. Next, we perform a sequence of unitary transformation generated by yielding with α = √ 2θ and The 2×2 matrices H + and H − describe a Dirac fermion with mass ±m in the presence of random vector potential ±(A x , A y ) and random scalar potential A 0 , each of which is an effective Hamiltonian for the plateau transition of integer quantum Hall effect [32,33].
The H ± sectors are coupled by the matrix element ασ 0 . The 4 × 4 continuum Dirac Hamiltonian H can be written in the form where τ 0 is a unit 2 × 2 matrix and τ x , τ y , and τ z are three Pauli matrices. The Hamiltonian (47) is invariant for each realization of disorder under the operation that implements time-reversal for a spin-1/2 particle. The Dirac Hamiltonian (47) is the main result of this subsection. It is an effective model for the Anderson localization of quantum spin Hall systems, which belongs to the symplectic class in view of the symmetry property (48). The Anderson transition in the Dirac Hamiltonian (47) should possess the same universal critical properties as those found in our numerical simulations of the Z 2 network model. In the presence of the "Rashba" coupling α, there should appear a metallic phase near m = 0 which is surrounded by two insulating phases. In the limit α → 0, the metallic phase should shrink into a critical point of the integer quantum Hall plateau transition.
The 4 × 4 continuum Dirac Hamiltonian H should be contrasted with a 2 × 2 Hamiltonian of a Dirac particle in random scalar potential, which has the minimal dimensionality of the Clifford algebra in (2 + 1)-dimensional space time and is invariant under time-reversal operation, σ y H * 2 σ y = H 2 . The 2×2 Dirac Hamiltonian (49) is an effective Hamiltonian for massless Dirac fermions on the surface of a three-dimensional Z 2 topological insulator. After averaging over the disorder potential V , the problem of Anderson localization of the surface Dirac fermions is reduced to a NLSM with a Z 2 topological term [25,26]. Interestingly, this Z 2 topological term prevents the surface Dirac fermions from localizing [27,28]. It is this absence of twodimensional localization that defines a three-dimensional Z 2 topological insulator [29]. In contrast, the doubling of the size of the Hamiltonian (47) implies that the NLSM describing the Anderson localization in the 4 × 4 Hamiltonian (45) does not come with a Z 2 topological term, because two Z 2 topological terms cancel each other. We can thus conclude that the critical properties of metal-insulator transitions in the Z 2 network model are the same as those in the standard symplectic class, in agreement with results of our numerical simulations of the Z 2 network model [19,20].
Before closing this subsection, we briefly discuss the Dirac Hamiltonian (47) in the clean limit where A 0 = A x = A y = 0. Since the system in the absence of disorder is translationally invariant, momentum is a good quantum number. We thus consider the Hamiltonian in momentum space where the wave number k = (k x , k y ). When α = 0, the Hamiltonian (50) becomes a direct sum of 2 × 2 Dirac Hamiltonian with mass of opposite signs. This is essentially the same low-energy Hamiltonian as the one appearing in the quantum spin Hall effect in HgTe/(Hg,Cd)Te quantum wells [6].

Z 2 topological number
We now discuss the topological property of the time-reversal invariant insulator which is obtained from the effective Hamiltonian (50) of the Z 2 network model in the absence of disorder. The topological attribute of the band insulator is intimately tied to the invarianceΘ under the operation of time-reversal represented bŷ where K implements complex conjugation. We are going to show that this topological attribute takes values in Z 2 , i.e., the Z 2 index introduced by Kane and Mele [4].
We begin with general considerations on a translation-invariant single-particle fermionic Hamiltonian which has single-particle eigenstates labeled by the wave vector k taking values in a compact manifold. This compact manifold can be the first Brillouin zone with the topology of a torus if the Hamiltonian is defined on a lattice and periodic boundary conditions are imposed, or it can be the stereographic projection between the momentum plane R 2 and the surface of a three-dimensional sphere if the Hamiltonian is defined in the continuum. We assume that (i) the antiunitary operation Θ = −Θ −1 = −Θ † that implements time-reversal leaves the Hamiltonian invariant, (ii) there exists a spectral gap at the Fermi energy, and (iii) there are two distinct occupied bands with the single-particle orthonormal eigenstates |uâ(k) and energies Eâ i.e., the overlaps between the occupied single-particle energy eigenstates with momentum −k and the time reversed images to the occupied single-particle energy eigenstates with momentum k, plays an important role [17]. The matrix elements (53) obey We used the fact thatΘ is antilinear to reach the second equality and that it is antiunitary withΘ 2 = −1 to reach the third equality. Hence, the 2 × 2 unitary sewing matrix w(k) with the matrix elements (53) can be parametrized as with the three complex-valued functions We observe that w(k) reduces to for some real-valued f (k) at any time-reversal invariant wave vector k ∼ −k (timereversal invariant wave vectors are half a reciprocal vector for a lattice model, and 0 or ∞ for a model in the continuum).
As we shall shortly see, the sewing matrix (53) imposes constraints on the U(2) Berry connection where the summation convention over the repeated index µ is understood (we do not make distinction between superscript and subscript). Here, at every point k in momentum space, we have introduced the U(2) antihermitian gauge field A µ (k) with the space index µ = 1, 2 and the matrix elements labeled with the U(2) internal indicesâ,b = 1, 2, by performing an infinitesimal parametric change in the Hamiltonian. We decompose the U(2) gauge field (58) into the U(1) and the SU (2) contributions where ρ 0 is a 2 × 2 unit matrix and ρ is a 3 vector made of the Pauli matrices ρ x , ρ y , and ρ z . Accordingly, Combining the identityΘ 2 = −1 with the (partial) resolution of the identity â=1,2 |uâ(k) uâ(k)| for the occupied energy eigenstates with momentum k yields where the proper restriction to the occupied energy eigenstates is understood for the unit operator on the right-hand side. Using this identity, we deduce the gauge transformation that relates the U(2) connections at ±k. For the U(1) and SU(2) parts of the connection, where we have decomposed w(k) into the U(1) (e iζ ) and SU(2) (w) parts according to (note that this decomposition has a global sign ambiguity, which, however, will not affect the following discussions).
Equipped with these gauge fields, we introduce the U(2) Wilson loop where the U(1) Wilson loop is given by while the SU(2) Wilson loop is given by The symbol P in the definition of the U(2) Wilson loop represents path ordering, while C is any closed loop in the compact momentum space.

By construction, the U(2) Wilson loop (67) is invariant under the transformation
the SU(2) Wilson loop W SU (2) [C] is quantized to the two values because of time-reversal symmetry. Furthermore, the identity which we will prove below, follows. Here, the symbol Pf denotes the Pfaffian of an antisymmetric matrix, and only the subset of momenta K ∈ C that are unchanged under K → −K contribute to the SU(2) Wilson loop. According to (54), the sewing matrix at a time-reversal symmetric wave vector is an antisymmetric 2 × 2 matrix. Consequently, the SU(2) part of the sewing matrix at a time-reversal symmetric wave vector is a real-valued antisymmetric 2 × 2 matrix (i.e., it is proportional to iρ y up to a sign). Hence, its Pfaffian is a well-defined and nonvanishing real-valued number. Before undertaking the proof of (74), more insights on this identity can be obtained if we specialize to the case when the Hamiltonian is invariant under any U(1) subgroup of SU(2), e.g., the z-component of spin σ z . In this case we can choose the basis states which diagonalize σ z ; σ z |u 1 (k) = +|u 1 (k) , σ z |u 2 (k) = −|u 2 (k) . Since the timereversal operation changes the sign of σ z , the sewing matrix takes the form which, in combination with (64) and (65) implies the transformation laws We conclude that when both the z component of the electron spin and the electron number are conserved, we can set and use the transformation law With conservation of the z component of the electron spin in addition to that of the electron charge, the SU(2) Wilson loop becomes We have used the fact that σ z is traceless to reach the last line. This line integral can be written as the surface integral by Stokes' theorem. Here, D is the region defined by ∂D = C, and covers a half of the total Brillouin zone (BZ) because of the condition (72). In turn, this surface integral is equal to the Chern number for up-spin fermions, where we have used the transformation law (79) to deduce that to reach the fourth equality.
To summarize, when the z component of the spin is conserved, the quantized SU(2) Wilson loop can then be written as the parity of the spin Chern number (the Chern number for up-spin fermions, which is equal to minus the Chern number for down-spin fermions) [3,4,5], Next, we apply the master formula (74) to the 4 × 4 Dirac Hamiltonian (50). To this end, we first replace the mass m by the k-dependent mass, and parametrize the wave number k as Without loss of generality, we may assume α > 0. The mass m k is introduced so that the SU(2) part of the sewing matrix is single-valued in the limit |k| → ∞. We then perform another series of unitary transformation with to rewrite the Hamiltonian (50) in the form The four eigenvalues of the Hamiltonian (90) are given by The occupied eigenstate with the energy E 1 (k) = −λ − k reads Table 1. θ ± k at the time-reversal invariant momenta k = 0 and k = ±∞, when m − Cα 2 < 0 (a) and m − Cα 2 > 0 (b). It is assumed that 0 ≤ arctan(|m|/α) ≤ π/2.
Notice that |u 2 (ϕ, k) = |u 1 (ϕ + π, −k) . The 2 × 2 sewing matrix w(k) is obtained from the eigenstates (92)-(93) as which is decomposed into the U(1) part, and the SU(2) part, of the sewing matrix. Here, we have defined θ ± k through the relation For the SU (2)  Pfaffian of the sewing matrix at the south and north poles of the stereographic sphere are respectively. Hence, for any time-reversal invariant path C passing through the south and north poles, where we have suppressed Cα 2 by taking the limit Cα 2 /|m| → 0. The value (100) taken by the SU (2) Wilson loop thus appears to be ambiguous since it depends on the sign of the mass m. This ambiguity is a mere reflection of the fact that, as noted in [20], the topological nature of the Z 2 network model is itself defined relative to that of some reference vacuum. Indeed, for any given choice of the parameters (X, θ) from figure 2 that defines uniquely the bulk properties of the insulating phase in the Z 2 network model, the choice of boundary conditions determines if a single helical Kramers' doublet edge state is or is not present at the boundary of the Z 2 network model. In view of this, it is useful to reinterpret the Z 2 network model with a boundary as realizing a quantum spin Hall droplet immersed in a reference vacuum as is depicted in figure 3(a). If so, choosing the boundary condition is equivalent to fixing the topological attribute of the reference vacuum relative to that of the Z 2 network model, for the reference vacuum in which the quantum spin Hall droplet is immersed also has either a trivial or nontrivial Z 2 quantum topology. A single helical Kramers' doublet propagating unhindered along the boundary between the quantum spin Hall droplet and the reference vacuum appears if and only if the Z 2 topological quantum numbers in the droplet and in the reference vacuum differ.
In the low-energy continuum limit (50), a boundary in real space can be introduced by breaking translation invariance along the vertical line x = 0 in the real space  (x, y) ∈ R 2 through the profile [see figure 3 for the mass.
We close Sec. 3 with a justification of the master formula (74). To this end, we regularize the continuum gauge theory by discretizing momentum space ( figure 4). We use the momentum coordinate i ∈ Z 2 on a rectangular grid with the two lattice spacings ∆k µ > 0. To each link from the site i to the nearest-neighbor site i + µ of the grid, we assign the SU(2) unitary matrix which is obtained by discarding U(1) part of the U(2) Berry connection. Consistency demands that We define the SU(2) Wilson loop to be where i n and i n+1 are nearest neighbors, i.e., their difference i n+1 − i n = η n is a unit vector η n . The Wilson loop is invariant under any local gauge transformation by which where the V i 's are U(2) matrices. Observe that the cyclicity of the trace allows us to write To make contact with the master formula (74), we assume that the closed path with vertices i ℓ parametrized by the index ℓ = 0, 1, · · · , N − 1 obeys the condition that . . .
i n is the wave vector + n m=1 η m ∆k µ m , . . . with η m = ±1 for m = 1, · · · , N/2 in order to mimic after discretization the condition that the closed path entering the Wilson loop is invariant as a set under the inversion (72).
On the discretized momentum lattice the sewing matrix (53) is defined by which obeys the condition i.e., the counterpart to the relation (54). This implies that w i 0 and w i N/2 are antisymmetric unitary 2 × 2 matrices. Furthermore, the sewing matrix w i must also obey the counterpart to (63), namely It now follows from (107) and (110) that . . .
In particular, we observe that since w i 0 is a 2 × 2 antisymmetric unitary matrix, i.e., w i 0 is the second Pauli matrix up to a phase factor, while holds for any three-vector n contracted with the three-vector ρ made of the three Pauli matrices. By repeating the same exercise a second time, one convinces oneself that the dependences on the gauge fields A i 0 ,i 1 and A i N−1 ,i 0 , A i 1 ,i 2 and A i N−2 ,i N−1 , and so on until A i n−1 ,i n and A i N−n ,i N−n+1 at the level n of this iteration cancel pairwise due to the conditions (107)-(110) implementing time-reversal invariance. This iteration stops when n = N/2, in which case the SU(2) Wilson loop is indeed solely controlled by the sewing matrix at the time-reversal invariant momenta corresponding to ℓ = 0 and ℓ = N/2, Since i 0 and i N/2 are invariant under momentum inversion or, equivalently, timereversal invariant, with α N/2 , α 0 = 0, π. Here, the Z 2 phases e iα N/2 and e iα 0 are none other than the Pfaffians respectively. Hence, is a special case of (74). (Recall that w i 0 and w i N/2 are real-valued.)

Numerical study of boundary multifractality in the Z 2 network model
In [20], we have shown that (i) multifractal scaling holds near the boundary of the Z 2 network model at the transition between the metallic phase and the Z 2 topological insulating phase shown in figure 2, (ii) it is different from that in the ordinary symplectic class, while (iii) bulk properties, such as the critical exponents for the divergence of the localization length and multifractal scaling in the bulk, are the same as those in the conventional two-dimensional symplectic universality class of Anderson localization. This implies that the boundary critical properties are affected by the presence of the helical edge states in the topological insulating phase adjacent to the critical point.
In this work, we improve the precision for the estimate of the boundary multifractal critical exponents. We also compute numerically additional critical exponents that encode corner (zero-dimensional) multifractality at the metal-to-Z 2 -topological-insulator transition. We thereby support the claim that conformal invariance is present at the metal-to-Z 2 -topological-insulator transition by verifying that conformal relations between critical exponents at these boundaries hold.

Boundary and corner multifractality
To characterize multifractal scaling at the metal-insulator transition in the Z 2 network model, we start from the time-evolution of the plane waves along the links of the network with the scattering matrices defined in (1)- (7) at the nodes S and S ′ . To minimize finite size effects, the parameter θ in (1)- (7) is chosen to be a random variable as explained in Sec. 2. We focus on the metal-insulator transition at X = X l = 0.971 as shown in figure 2(b). When we impose reflecting boundary conditions, a node on the boundary reduces to a unit 2 × 2 matrix. When the horizontal reflecting boundaries are located at nodes of type S ′ , as shown in figure 5(a), there exists a single helical edge states for X > X l . The insulating phase X > X l is thus topologically nontrivial.
For each realization of the disorder, we numerically diagonalize the one-step timeevolution operator of the Z 2 network model and retain the normalized wave function ψ σ (x, y), after coarse graining over the 4 edges of the plaquette located at (x, y), whose eigenvalue is the closest to 1. The wave function at criticality is observed to display the power-law dependence on the linear dimension L of the system, The anomalous dimension ∆ (ζ,ν) q , if it displays a nonlinear dependence on q, is the signature of multifractal scaling. The index ζ indicates whether the multifractal scaling applies to the bulk (ζ = 2), the one-dimensional boundary (ζ = 1), or to the zerodimensional boundary (corner) (ζ = 0), provided the plaquette (x, y) is restricted to the corresponding regions of the Z 2 network model. For ζ = 1 and 0, the index ν distinguishes the case ν = O when the ζ-dimensional boundary has no edge states in the insulating phase adjacent to the critical point, from the case ν = Z 2 when the ζ-dimensional boundary has helical edge states in the adjacent insulating phase. We ignore this distinction for multifractal scaling of the bulk wave functions, ∆ q , since bulk properties are insensitive to boundary effects. We will also consider the case of mixed boundary condition for which we reserve the notation ν = Z 2 |O.
It was shown in [31] that boundary multifractality is related to corner multifractality if it is assumed that conformal invariance holds at the metal-insulator transition in the two-dimensional symplectic universality class. Conversely, the numerical verification of Figure 5.
(a) Boundary multifractality is calculated from the wave function amplitudes near a one-dimensional boundary.
Periodic (reflecting) boundary conditions are imposed for the horizontal (vertical) boundaries.
(b) Corner multifractality is calculated from the wave function amplitudes near a corner with the wedge angle ϑ = π/2. Reflecting boundary conditions are imposed along both vertical and horizontal directions. The relationship between the scattering matrix at a node of type S ′ and the scattering matrix at a node of type S implies that it is a vertical boundary located at nodes of type S that induces an helical edge state when X > X l . this relationship between boundary and corner multifractality supports the claim that the critical scaling behavior at this metal-insulator transition is conformal. So we want to verify numerically if the consequence of the conformal map w = z ϑ/π , namely where ϑ is the wedge angle at the corner, holds. Equivalently, f (ζ,ν) (α), which is defined to be the Legendre transformation of ∆ (ζ,ν) q + dq, i.e., must obey if conformal invariance is a property of the metal-insulator transition.
To verify numerically the formulas (120), (123), and (124), we consider the Z 2 network model with the geometries shown in figure 5. We have calculated wave functions for systems with the linear sizes L = 50, 80, 120, 150, and 180 for the two geometries displayed in figure 5. Here, L counts the number of nodes of the same type along a boundary. The number of realizations of the static disorder is 10 5 for each system size.  derived analytically in [34]. Since the triangles and circles are consistent within error bars, our numerical results are reliable, especially between 0 < q < 1. If we use the numerical values of ∆ (1,Z 2 ) q as inputs in (120) with ϑ = π/2, there follows the corner multifractal scaling exponents that are plotted by the solid curve. Since the curve overlaps with the direct numerical computation of ∆ (0,Z 2 ) q within the error bars, we conclude that the relation (120) is valid at the metal-to-Z 2 -topological-insulator transition.
The value of α is consistent with that reported in [20], while its accuracy is improved. The solid curve obtained from the relations (123) and (124) by using f (1,Z 2 ) (α) as an input, coincides with f (0,Z 2 ) (α). We conclude that the hypothesis of conformal invariance at the quantum critical point of metal-to-Z 2 -topological-insulator transition is consistent with our numerical study of multifractal scaling. (a) The z dependence of ln |Ψ| 2 z,L at the metal-to-Z 2 -topologicalinsulator transition in the cylindrical geometry for L = 50, 80, 120, 150, 180 from the top to the bottom. (b) The z dependence ofα 0 (z) (•)and c(z) ( ) extrapolated from the system size dependence of ln |Ψ| 2 z,L averaged over a small interval of z's. α 0 (z) and c(z) at z = 0, 1 without averaging over a small interval of z's are shown by open circles and squares, respectively. The solid line represents the bulk value of α (2) 0 = 2.173 computed in [31]. The asymmetry with respect to z = 0.5 is due to statistical fluctuations.
At last, we would like to comment on the dependence on z of found in [20]. Here, z ≡ (x − 1)/2L, while x and y denote the positions on the network along its axis and along its circumference, respectively (our choice of periodic boundary conditions imposes a cylindrical geometry). The overline denotes averaging over disorder. Figure 7(a) shows the z dependence of ln |Ψ| 2 z,L for different values of L in this cylindrical geometry at the metal-to-Z 2 -topological-insulator transition. We observe that ln |Ψ| 2 z,L becomes a nonmonotonic function of z. We are going to argue that this nonmonotonic behavior is a finite size effect. We make the scaling ansatz where ζ = 1 if z = 0, 1 and ζ = 2 otherwise, while c(z) depends on z but not on L. To check the L dependence of ln |Ψ| 2 z,L in figure 7, we average ln |Ψ| 2 z,L over a narrow interval of z's for each L. Figure 7(b) shows the z dependence ofα 0 (z) (•) and c(z) ( ) obtained in this way. In addition,α 0 (z) and c(z) calculated for z = 0, 1 without averaging over the narrow interval of z's are shown by open circles and open squares, respectively.
We observe thatα 0 (z), if calculated by averaging over a finite range of z's, is almost constant and close to α (2) 0 = 2.173. In contrast,α 0 (z = 0, 1) ≈ 2.09, if calculated without averaging over a finite range of z's, is close to α (1,Z 2 ) 0 = 2.091. We also find that |c(z)| increases near the boundaries. We conclude that it is the nonmonotonic dependence of |c(z)| on z that gives rise to the nonmonotonic dependence of ln |Ψ| 2 z,L on z. This finite-size effect is of order 1/ ln L and vanishes in the limit L → ∞.

Boundary condition changing operator
Next, we impose mixed boundary conditions by either (i) coupling the Z 2 network model to an external reservoir through point contacts or (ii) by introducing a longrange lead between two nodes from the Z 2 network model, as shown in figure 8. In this way, when X > X l , a single Kramers' pair of helical edge states indicated by the wavy lines in figure 8 is present on segments of the boundary, while the complementary segments of the boundary are devoid of any helical edge state (the straight lines in figure 8). The helical edge states either escape the Z 2 network model at the nodes at which leads to a reservoir are attached [the green lines in figure 8(a) and figure 8(b)], or shortcut a segment of the boundary through a nonlocal connection between the two nodes located at the corners [figure 8(c)]. These are the only options that accommodate mixed boundary conditions and are permitted by time-reversal symmetry. As shown by Cardy in [35], mixed boundary conditions are implemented by boundary-conditionchanging operators in conformal field theory. Hence, the geometries of figure 8 offer yet another venue to test the hypothesis of two-dimensional conformal invariance at the metal-insulator quantum critical point.
When coupling the Z 2 network model to an external reservoir, we shall consider two cases shown in figure 8(a) and figure 8(b), respectively.
First, we consider the case of figure 8(a) in which only two nodes from the Z 2 network model couple to the reservoir. At each of these point contacts, the scattering matrix S that relates incoming to outgoing waves from and to the reservoir is a 2 × 2 matrix which is invariant under time reversal s 2 S * s 2 = S, and it must be proportional to the unit 2 × 2 matrix up to an overall (random) phase. Hence, the two-point-contact conductance in the geometry of figure 8(a) is unity, however far the two point contacts are from each other.
Second, we consider the case of figure 8(b) in which there are again two point contacts, however each lead between the Z 2 network model and the reservoir now supports two instead of one Kramers' doublets. The two point-contact scattering matrices connecting the Z 2 network model to the reservoirs are now 4 × 4 matrices, which leads to a non-vanishing probability of backscattering. Hence, this even channel two-point-contact conductance is expected to decay as a function of the separation between the two attachment points of the leads to the network.
To test whether the two-point-contact conductance in figure 8(a) and figure 8(b) do differ as dramatically as anticipated, we have computed numerically the two-pointcontact conductance at the quantum critical point X = X l for a single realization of the static disorder. The two-point-contact conductance is calculated by solving for the stationary solution of the time-evolution operator with input and output leads [36]. We choose the cylindrical geometry imposed by periodic boundary conditions along the  Figure 9(a) shows with the symbol • the dependence on r, the distance between the two contacts in figure 8(a), of the dimensionless two-point-contact conductance g. It is evidently r independent and unity, as expected. Figure 9(a) also shows with the symbol • the dependence on r of the dimensionless two-point conductance g for leads supporting two Kramers' doublet as depicted in figure 8(b). Although it is not possible to establish a monotonous decay of the two-point-contact conductance for a single realization of the static disorder, its strong fluctuations as r is varied are consistent with this claim.
We turn our attention to the closed geometry shown in figure 8(c). We recall that it is expected on general grounds that the moments of the two-point conductance in a network model at criticality, when the point contacts are far apart, decay as power laws with scaling exponents proportional to the scaling exponents ∆ (ζ,ν) q [36,37]. Consequently, after tuning the Z 2 network model to criticality, the anomalous dimensions at the node (corner) where the boundary condition is changed must vanish, since the two-point-contact conductance in figure 8(a) is r independent. (130) is another signature of the nontrivial topological nature of the insulating side at the Anderson transition that we want to test numerically. Thus, we consider the geometry figure 8(c) and compute numerically the corner anomalous dimensions. This is done using the amplitudes of the stationary wave function restricted to the links connecting the two corners where the boundary conditions are changed. Figure 9(b) shows the numerical value of the corner anomalous dimension ∆ (0,Z 2 |O) q . The linear sizes of the network are L = 50, 80, 120, 150, and 180 and the number of disorder realizations is 10 5 for each L. We observe that ∆ (0,Z 2 |O) q is zero within the error bars, thereby confirming the validity of the prediction (130) at the metal-to-Z 2 -topological-insulator transition.

Conclusions
In summary, we have mapped the Z 2 network model to a 4 × 4 Dirac Hamiltonian. In the clean limit of this Dirac Hamiltonian, we expressed the Kane-Mele Z 2 invariant as an SU(2) Wilson loop and computed it explicitly. In the presence of weak timereversal symmetric disorder, the NLSM that can be derived out of this Dirac Hamiltonian describes the metal-insulator transition in the Z 2 network model and yields bulk scaling exponents that belong to the standard two-dimensional symplectic universality class; an expectation confirmed by the numerics in [19] and [20]. A sensitivity to the Z 2 topological nature of the insulating state can only be found by probing the boundaries, which we did numerically in the Z 2 network model by improving the quality of the numerical study of the boundary multifractality in the Z 2 network model.