text:
# Network models for localisation problems belonging to the chiral symmetry classes

## I Introduction

The classification by symmetry of Hamiltonians for disordered systems provides an important framework in the study of Anderson localisation. Three standard symmetry classes are long-established, and represented in the zero-dimensional limit by the three Wigner-Dyson random matrix ensembles.Mehta1991  The existence of additional symmetry classes has been implicit in work which also has an extensive history, but a complete classification has been set out only rather recently, by Altland and Zirnbauer.Zirnbauer1996 ; Altland1997  Hamiltonians belonging to one of these additional symmetry classes are characterised by a discrete symmetry which relates eigenvalues and eigenfunctions in pairs. This paper is concerned with three examples, known as the chiral symmetry classes. They can be realised as tight-binding Hamiltonians for systems in which sites can be divided into two sublattices and non-zero matrix elements connect only sites from opposite sublattices. In these terms, the discrete symmetry operation is multiplication of wavefunction components on one sublattice by a factor of minus one: applied to an eigenstate with energy $`E`$, this generates a second eigenstate, with energy $`E`$. Three distinct such symmetry classes are possible, termedAltland1997  AIII for systems without time-reversal invariance, BDI for spinless systems with time-reversal invariance, and CII for systems having time-reversal but not spin-rotation symmetry.

Early work on disordered Hamiltonians with chiral symmetry followed a variety of independent directions. A one-dimensional model of this type was solved exactly by Dyson Dyson1953 , who found a divergence in the density of states at the band centre, of a form shown by Ovchinnikov and Erikhman Ovchinnikov1977  to be characteristic for the symmetry class in one dimension. A closely related feature of this and similar models is that, while states away from the band centre are localised, both the mean free path and the localisation length diverge on approaching the band centre, as revealed in early treatments by Eggarter and Riedinger Eggarter1978 and by Ziman.Ziman1982  Two-dimensional systems with chiral symmetry were studied by Gade and Wegner, Gade1991 ; Gade1993  using a mapping to a non-linear sigma model: their calculations for the class AIII with weak disorder indicate a strongly divergent density of states at the band centre, also accompanied by a divergence of the localisation length. Chiral random matrix ensembles, investigated by Nagao and Slevin Nagao1993  and by Verbaarschot and Zahed,Verbaarschot1993  likewise show band centre anomalies in the density of states, albeit on energy scales comparable with the level spacing. And in a separate development, extensive analytical results have been obtained by Nersesyan et al Nersesyan1995  and by Ludwig et al Ludwig1994  for a particular two-dimensional model with chiral symmetry, the massless two-dimensional Dirac equation with random vector potential.

Much subsequent work has helped to establish the range of behaviour possible in these systems. In one dimensional models, a finite localisation length at zero energy and a modification of the Dyson singularity in the density of states can be induced Ovchinnikov1977 ; Comtet1995 ; Brouwer1998  by terms in the Hamiltonian that in the absence of disorder would generate a spectral gap. For two-dimensional systems, it is recognised that a random vector potential in the Dirac equation constitutes a special choice of disorder, and scaling flow from this towards a generic fixed point, with properties similar to those derived by Gade, has been studied in a field-theoretic framework Guruswamy2000  and by using results on random Gaussian surfaces. Motrunich2002  Moreover, while early numerical results Evangelou1986 ; Morita1997 ; Furusaki1999  were quite different from analytical predictions, careful tuning of model parameters and study of large systems has recently brought calculations closer to expected asymptotic behaviour.Motrunich2002 ; Ryu2002  Numerical work has also shown that a finite value for the zero-energy localisation length can be produced in two-dimensional chiral models, by the same mechanism that is effective in one dimension.Motrunich2002

Our aim in the following is to formulate and study new representatives of the chiral symmetry classes, in the form of network models.Chalker1988  In place of a Hamiltonian, these use the ideas of scattering theory and may be specified by a transfer matrix,Chalker1988  or by a unitary evolution operator for one step of discretised time.Klesse1995 ; Ho1996  The versions we set out here are connected in two distinct ways to particular network models without chiral symmetry which have been studied previously. For the class AIII, these connections are to the U(1) network model Chalker1988  investigated in the context of the integer quantum Hall plateau transition, while for CII and BDI, respectively, they are to the SU(2) and O(1) network models, which describe plateau transitions in dirty superconductors.Kagalovsky1997 ; Kagalovsky1999 ; Chalker2002a  In each case, the model with chiral symmetry is constructed by coupling two copies of the partner model. This coupling can be re-expressed after a transformation as equivalent to introduction of absorption and coherent amplification in the original models. This equivalence parallels the established link Simons1998  between chiral symmetry classes and non-Hermitian random operators.

The remainder of this paper is organised as follows. We show in Sec. II how to construct network models belonging to the chiral symmetry classes, and discuss in Sec. III some basic aspects, including the relation to non-Hermitian models, the continuum limit, and a lattice version of disorder analogous to a random vector potential in the Dirac equation. We present results from a numerical study of a two-dimensional model in the symmetry class AIII in Sec. IV, and summarise in Sec. V.

## II Construction of network models with chiral symmetry

In this section we formulate network models for each of the chiral symmetry classes, giving a detailed treatment of the chiral unitary class (AIII), and indicating in outline the equivalent steps for the chiral orthogonal (BDI) and chiral symplectic (CII) classes. Our strategy is simply to construct the two-dimensional internal space associated with chiral symmetry using two related copies of a network model without that symmetry. For each symmetry class we focus on two-dimensional models; models in quasi-one and three dimensions can be constructed in the same way, starting for example from Eq. (7).

The symmetry of a disordered system may be discussed in terms of a Hamiltonian $`H`$, a scattering matrix $`S`$, or a transfer matrix $`T`$. We are concerned with systems which conserve probability density, so that the Hamiltonian is Hermitian and the scattering matrix is unitary: $`H^{}=H`$ and $`S^{}=S^1`$. The equivalent condition for the transfer matrix involves the current operator $`J`$ and reads $`T^1=JT^{}J`$. Chiral symmetry is implemented on a two-dimensional internal space in which the Pauli operator $`\mathrm{\Sigma }_x`$ acts: for the Hamiltonian it is the requirement that $`\mathrm{\Sigma }_xH\mathrm{\Sigma }_x=H`$. Taking the scattering matrix to have the symmetry of $`e^{iH}`$, this implies that $`\mathrm{\Sigma }_xS\mathrm{\Sigma }_x=S^1`$. For the transfer matrix, it appears at first that one has a choice, taking flux in the pairs of scattering channels on which $`\mathrm{\Sigma }_x`$ acts to propagate either in opposite directions or in the same direction. In fact only the former is tenable. It implies that the anticommutator $`\{J,\mathrm{\Sigma }_x\}=0`$ and that the chiral symmetry condition is $`\mathrm{\Sigma }_xT\mathrm{\Sigma }_x=T`$, a property preserved under matrix multiplication. By contrast, the latter choice would lead to the commutator $`[\mathrm{\Sigma }_x,T]=0`$ and the chiral symmetry condition $`\mathrm{\Sigma }_xT\mathrm{\Sigma }_x=T^1`$. Since the last condition is not in general preserved under matrix multiplication, connection in series of two systems of this kind would generate a sample without chiral symmetry, and we therefore reject this alternative. Systems in class AIII have no other relevant discrete symmetry; those in classes BDI and CII are also invariant under time-reversal, in the absence and presence of Kramers degeneracy respectively.

### II.1 Network model for class AIII

We recall first the essential features of the U(1) network model for the integer quantum Hall plateau transition. A wavefunction in this model takes complex values $`z_l`$ on the links $`l`$ of the lattice illustrated by the full lines of Fig. 2. The forms of the transfer matrix and of the evolution operator follow from the properties of the elementary building units, shown in Fig. 1. A particle acquires a phase $`\varphi _l`$ on traversing link $`l`$, so that in a stationary state amplitudes at either end are related by

$$z^{}=e^{i\varphi }z.$$
(1)
In a similar way, stationary state amplitudes on the four links which meet at a node are related by a $`2\times 2`$ transfer matrix

$$\left(\begin{array}{c}z_1\\ z_2\end{array}\right)=\left(\begin{array}{cc}\mathrm{cosh}(a)& \mathrm{sinh}(a)\\ \mathrm{sinh}(a)& \mathrm{cosh}(a)\end{array}\right)\left(\begin{array}{c}z_3\\ z_4\end{array}\right),$$
(2)
where $`a`$ is real and all phase factors are associated with links. This equation may re-written in terms of a scattering matrix as

$$\left(\begin{array}{c}z_3\\ z_2\end{array}\right)=\left(\begin{array}{cc}\mathrm{cos}(\alpha )& \mathrm{sin}(\alpha )\\ \mathrm{sin}(\alpha )& \mathrm{cos}(\alpha )\end{array}\right)\left(\begin{array}{c}z_1\\ z_4\end{array}\right),$$
(3)
with $`\mathrm{sin}(\alpha )=\mathrm{tanh}(a)`$. The transfer matrix that results from assembling these units is described in detail in Ref. Chalker1988, , and the time evolution operator in Refs. Klesse1995,  and Ho1996, .

Introduction of a two-dimensional internal space associated with chiral symmetry results in a doubling of the number of wavefunction components. In this way, starting from two copies of the U(1) model the link amplitudes become two-component complex numbers, $`𝐳_l`$. In place of Eq. (1), the scattering properties of a link are characterised by a $`2\times 2`$ transfer matrix $`T`$, with

$$𝐳^{}=T𝐳.$$
(4)
Requiring $`\mathrm{\Sigma }_xT\mathrm{\Sigma }_x=T`$, $`T`$ has the form

$$T=e^{i\varphi }\left(\begin{array}{cc}\mathrm{cosh}(b)& \mathrm{sinh}(b)\\ \mathrm{sinh}(b)& \mathrm{cosh}(b)\end{array}\right),$$
(5)
where $`\varphi `$ is a real phase and $`b`$ is a real hyperbolic angle. It remains to discuss scattering at nodes of the doubled system. We replace Eq. (2) by

$$\left(\begin{array}{c}𝐳_1\\ 𝐳_2\end{array}\right)=𝟙\left(\begin{array}{cc}\mathrm{cosh}(a)& \mathrm{sinh}(a)\\ \mathrm{sinh}(a)& \mathrm{cosh}(a)\end{array}\right)\left(\begin{array}{c}𝐳_3\\ 𝐳_4\end{array}\right),$$
(6)
where $`𝟙`$ is here the unit matrix in the two-component space introduced on links. This choice amounts to the most general one compatible with chiral symmetry, since all scattering within the two-component space may be included in the link transfer matrices, $`T`$.

Combining these elements to make a two-dimensional system, we arrive at the model shown schematically in Fig. 2. The transfer matrix for the system as a whole acts in the $`[1,1]`$ (or $`[1,\overline{1}]`$) direction and may be written as a product of factors relating amplitudes on successive slices of the system. Alternate factors in the product represent links and nodes, and consist respectively of repeated versions of the $`2\times 2`$ and $`4\times 4`$ blocks appearing in Eqs. (5) and (6). We introduce disorder by taking the phase $`\varphi `$ in Eq. (5) to be an independent, uniformly distributed random variable on each link, and take the parameter $`a`$, characterising scattering at nodes, to be non-random. We have considered two ways to set the value of the coupling between the chiral subspaces, which we parameterise in terms of either the hyperbolic angle appearing in Eq. (5) or the compact angle $`\beta `$ related to $`b`$ by $`\mathrm{sin}(\beta )=\mathrm{tanh}(b)`$: we take $`\beta `$ either uniformly distributed, or $`b`$ to have a normal distribution of variance $`g`$.footnote1  We ensure that the system is statistically invariant under $`\pi /2`$ rotations of the lattice, which places requirements on the node parameter $`a`$, exactly as in the U(1) model: nodes lie on two distinct sub-lattices, and the node parameter $`a`$ on one sublattice is related to the parameter $`a^{}`$ on the other sublattice by the duality relation $`\mathrm{sinh}(a)\mathrm{sinh}(a^{})=1`$.

The same system can alternatively be described using a time evolution operator, $`S`$, in place of a transfer matrix. In order to specify this operator, which is unitary and has the symmetry of a scattering matrix, it is convenient first to consider the special case $`b=0`$, in which the two copies of the network model are uncoupled. Let $`U`$ denote the time evolution operator for one copy. Then, since from Eq. (5) link phases are the same in both copies but propagation directions are opposite, the evolution operator for the other copy is $`U^{}`$. The dimension of $`U`$ is equal to the number of links in the system. It is useful to define a diagonal matrix of the same dimension, with the angles $`\beta _l`$ for each link $`l`$ as diagonal entries: we use $`\beta `$ to denote this matrix. The time-evolution operator for the system with chiral symmetry can then be written

$$S=\left(\begin{array}{cc}U\mathrm{cos}(\beta )& U\mathrm{sin}(\beta )U^{}\\ \mathrm{sin}(\beta )& \mathrm{cos}(\beta )U^{}\end{array}\right).$$
(7)
It is straightforward to check that $`\mathrm{\Sigma }_xS\mathrm{\Sigma }_x=S^1`$ and that $`S^{}=S^1`$.

### II.2 Network model for classes BDI and CII

A model in class BDI can be obtained from one in class AIII simply by imposing time reversal invariance as an additional symmetry. This condition is conventionally written in the form $`H^{}=H`$, but for a discussion based on scattering matrices it is more convenient to make the transformation $`HQHQ^1`$ with $`Q=\mathrm{exp}\left(i\frac{\pi }{4}\mathrm{\Sigma }_x\right)`$. This transformation leaves the chiral symmetry relation $`\mathrm{\Sigma }_xH\mathrm{\Sigma }_x=H`$ unchanged. In the transformed basis one has $`H^{}=H`$, $`S=S^{}`$ and $`T^{}=T`$. To ensure a real time evolution operator, we restrict the link phases $`\varphi `$ to the values $`0`$ and $`\pi `$. Choosing these values randomly, the BDI model consists of two coupled copies of the class D models studied recently in the context of disordered superconductors.Chalker2002a  Alternatively, one could set $`\varphi =0`$ on all links, and introduce disorder only through the chiral couplings $`\beta `$.

We turn next to the symmetry class CII. Kramers degeneracy is a defining feature of the class and necessitates the introduction of an additional two-dimensional space arising from spin. The time reversal operation includes reversal of spin direction. Defining $`𝒞=i\tau _y`$, where $`\tau _y`$ is a Pauli matrix acting in the additional space, it is conventionally written in the form $`𝒞H^{}𝒞^1`$=H. As for the class BDI, it is again convenient to make the transformation $`HQHQ^1`$ with $`Q=\mathrm{exp}\left(i\frac{\pi }{4}\mathrm{\Sigma }_x\right)`$. In the transformed basis one has $`𝒞H^{}𝒞^1=H`$, $`𝒞S^{}𝒞^1=S`$ and $`𝒞T^{}𝒞^1=T`$ as equivalent expressions of time-reversal invariance. Applying these ideas to a network model, four channels propagate on a single link which, generalising Eq. (4), has a $`4\times 4`$ transfer matrix $`T`$ with the generic form

$$T=v\left(\begin{array}{cc}\mathrm{cosh}(b)& \mathrm{sinh}(b)\\ \mathrm{sinh}(b)& \mathrm{cosh}(b)\end{array}\right),$$
(8)
where $`v`$ is an SU(2) matrix and $`b`$ is a real hyperbolic angle. Adopting this form, the time evolution operator for class CII has the structure given in Eq. (7), but with $`U`$ representing a class C network model, studied previously in connection with the spin quantum Hall effect.Kagalovsky1997 ; Kagalovsky1999  The links of a such a class C model carry two co-propagating channels, coming from two spin components, and the evolution operator satisfies $`𝒞U^{}𝒞^1=U`$.

## III Discussion of the models

In this section we discuss some basic aspects which are common to models from all three chiral symmetry classes. These include: symmetry of the spectrum of the time-evolution operator; the relation to network models which have absorption and amplification; the Green function; the continuum limit; behaviour for some special parameter values; and a lattice analogue of random vector potential disorder in the Dirac equation.

### III.1 Spectral symmetries

The eigenphases $`E_\alpha `$ and eigenvectors $`\psi _\alpha `$ of the time evolution operator $`S`$ play the same role for a network model as do energy eigenvalues and eigenvectors for a model defined using a Hamiltonian. They satisfy

$$S\psi _\alpha =e^{iE_\alpha }\psi _\alpha .$$
(9)
Chiral symmetry, $`\mathrm{\Sigma }_xS\mathrm{\Sigma }_x=S^1`$, has the consequence that $`\mathrm{\Sigma }_x\psi _\alpha `$ is an eigenvector of $`S`$ with opposite eigenphase:

$$S\mathrm{\Sigma }_x\psi _\alpha =e^{iE_\alpha }\mathrm{\Sigma }_x\psi _\alpha .$$
(10)
Because of this, the spectrum of eigenphases is symmetric around the point $`E=0`$ and their density $`\rho (E)`$ may develop a singularity there. In addition, for models from the class CII, each eigenphase is Kramers degenerate. Specifically, in addition to $`\psi _\alpha `$, $`S`$ has a second eigenvector, $`𝒞\mathrm{\Sigma }_x\psi _\alpha ^{}`$, with the same eigenvalue $`e^{iE_\alpha }`$.

The symmetries described in the preceding paragraph are the ones of fundamental interest in this work. However, all network models based on the lattice of Fig. 2 have in addition a bipartite structure. This structure leads to symmetry of the eigenphase spectrum under the translation $`E\pi +E`$. When combined with chiral symmetry, $`EE`$, the consequence is that behaviour near $`E=0`$ is mirrored exactly at $`E=\pi `$, and that similar behaviour is repeated at $`E=\pm \pi /2`$.

In detail, the bipartite structure (which is inheritedHo1996  from the network model for the quantum Hall plateau transition, and discussed further in Sec. III.4) has the consequence that the evolution operator $`U`$ appearing in Eq. (7) (or its counterparts for classes BDI and CII) can be written with a suitable ordering of the basis states as

$$U=\left(\begin{array}{cc}0& A\\ B& 0\end{array}\right).$$
(11)
Introducing in the same basis a diagonal matrix $`\sigma =\mathrm{diag}(+1,1)`$, bipartite lattice structure implies that $`\sigma U\sigma =U`$. Applying this to our models (and suitably extending the meaning of $`\sigma `$), one finds that the pairs of eigenvectors of $`S`$ identified above, $`\psi _\alpha `$ and $`\mathrm{\Sigma }_x\psi _\alpha `$, are related to pairs $`\sigma \psi _\alpha `$ and $`\sigma \mathrm{\Sigma }_x\psi _\alpha `$, with eigenvalues $`\pi +E_\alpha `$ and $`\pi E_\alpha `$ respectively. These relationships are illustrated in Fig. 3. They result in an eigenphase spectrum that is reflection-symmetric not only about the points $`E=0`$ and $`E=\pi `$ but also about the points $`E=\pm \pi /2`$.

### III.2 Chiral symmetry and network models with absorption and amplification

Any Hamiltonian $`H`$ which satisifes the chiral symmetry condition $`\mathrm{\Sigma }_xH\mathrm{\Sigma }_x=H`$ has a block structure that is most clearly displayed after making the rotation $`R=(\mathrm{\Sigma }_x+\mathrm{\Sigma }_z)/\sqrt{2}`$. In the rotated basis, the chiral symmetry condition involves $`\mathrm{\Sigma }_z`$ rather than $`\mathrm{\Sigma }_x`$, since $`\mathrm{\Sigma }_z=R\mathrm{\Sigma }_xR^1`$. We have

$$RHR^1=\left(\begin{array}{cc}0& 𝐡\\ 𝐡^{}& 0\end{array}\right).$$
(12)
This establishes a decomposition of the Hermitian operator $`H`$ into the operators $`𝐡`$ and $`𝐡^{}`$, which are in general non-Hermitian. We show in this section that for a network model with chiral symmetry there exists a similar decomposition into a pair of systems, each with half the number of degrees of freedom. A single member of the pair, taken in isolation, can be interpreted as a system in which there is absorption and amplification of flux. We discuss this separation first in the context of the transfer matrix, and then in terms of the time evolution operator.

The logic can also be applied in the opposite direction. In this case, just as Eq. (12) provides a useful relationhermitisation  between a non-Hermitian operator $`𝐡`$, and the Hermitian operator $`H`$, so the arguments below can be used to relate a network model without flux conservation to a doubled system in which flux is conserved.

#### III.2.1 Transfer matrices with absorption and amplification

For our models all factors entering the transfer matrix are diagonal in the two-dimensional space on which $`\mathrm{\Sigma }_x`$ acts, except for the $`2\times 2`$ matrices appearing in Eqs. (5) and (8). These matrices are diagonalised by the rotation $`R`$, which is disorder independent and which mixes counter-propagating channels:

$$R\left(\begin{array}{cc}\mathrm{cosh}(b)& \mathrm{sinh}(b)\\ \mathrm{sinh}(b)& \mathrm{cosh}(b)\end{array}\right)R^1=\left(\begin{array}{cc}e^b& 0\\ 0& e^b\end{array}\right).$$
(13)
Because this rotation is independent of disorder, we can apply it to the transfer matrix $`T`$ for the system as a whole. The result reads symbolically

$$RTR^1=\left(\begin{array}{cc}𝒯(b)& 0\\ 0& 𝒯(b)\end{array}\right).$$
(14)
In the case of our model with AIII symmetry, $`𝒯(b)`$ is obtained from the transfer matrix for the corresponding U(1) model by including imaginary parts in the link phases:

$$\varphi _l\varphi _lib_l.$$
(15)
In a similar way, for our models with BDI or CII symmetry, $`𝒯(b)`$ is obtained from the transfer matrices for class D or class C models respectively, by associating an extra factor of $`\mathrm{exp}(b_l)`$ with propagation along each link $`l`$.

If $`𝒯(b)`$ or $`𝒯(b)`$ is regarded as the transfer matrix of a system in its own right, then that system is one with absorption and amplification. In fact, of course, the rotation $`R`$ gives the current operator $`J`$ off-diagonal components which couple $`𝒯(b)`$ and $`𝒯(b)`$, so that $`T`$, properly considered, is current-conserving as it should be.

From Eq. (14), one sees that the spectrum of Lyapunov exponents for a network model with chiral symmetry consists of separate contributions, from $`𝒯(b)`$ and from $`𝒯(b)`$. Models in which the $`b_l`$ are random and distributed symmetrically about $`b_l=0`$ clearly have doubly degenerate Lyapunov exponents. These can be calculated by studying just one member of the pair $`𝒯(\pm b)`$, which is important because it halves the size of matrices involved in numerical calculations.

#### III.2.2 Time evolution operators with absorption and amplification

It is interesting to re-examine this separation of the chiral model into subsystems, starting from the time-evolution operator in place of the transfer matrix. We begin by writing the action of the time-evolution operator $`S`$ in a way that emphasises the $`2\times 2`$ block structure evident in Eq. (7):

$$\left(\begin{array}{c}Z_3\\ Z_2\end{array}\right)=S\left(\begin{array}{c}Z_1\\ Z_4\end{array}\right).$$
(16)
The dimensions of the vectors $`Z`$ and $`\overline{Z}`$ appearing in Eq. (16) are the same as those of the matrix $`U`$ in Eq. (7). This expression can be reorganised into the form

$$\left(\begin{array}{c}Z_3\\ Z_4\end{array}\right)=K\left(\begin{array}{c}Z_1\\ Z_2\end{array}\right),$$
(17)
where $`K`$ is

$$K=\left(\begin{array}{cc}U& 0\\ 0& U\end{array}\right)\left(\begin{array}{cc}\mathrm{cosh}(b)& \mathrm{sinh}(b)\\ \mathrm{sinh}(b)& \mathrm{cosh}(b)\end{array}\right).$$
(18)
Properly speaking, since $`Z_1`$ and $`Z_2`$ are amplitudes of oppositely propagating fluxes, $`K`$ should be interpreted as a transfer matrix, in the way illustrated in Fig. 4. Despite this, our approach is to apply to $`K`$ the rotation $`R`$, and to discuss the rotated matrix as if it were an evolution operator. Doing this, we find

$$RKR^1=\left(\begin{array}{cc}𝒰(b)& 0\\ 0& 𝒰(b)\end{array}\right).$$
(19)
where

$$𝒰(\pm b)=U\mathrm{exp}(\pm b).$$
(20)
For the special case $`b=0`$, the subsystems represented by full and dashed lines in Fig. 2 are uncoupled. In this case one has $`𝒰(0)=U`$, which is the time evolution operator for one subsystem, and its inverse for the other subsystem. For general $`b`$, one can regard $`𝒰(\pm b)`$ as time evolution operators for systems with absorption and amplification. However, it is only the stationary states of these systems which have a simple physical significance. Such states exist either for an open system with scattering boundary conditions (the situation represented by the transfer matrix $`T`$), or for a closed system if a pair of eigenphases are zero (which can be arranged by fine-tuning system parameters). In this case, we can define eigenvectors

$$Z_+=𝒰(b)Z_+\text{and}Z_{}=𝒰(b)Z_{}.$$
(21)
Undoing the rotation $`R`$, the stationary states of $`S`$ satisfy $`\psi _\pm =S\psi _\pm `$, with

$$\psi _+=\left(\begin{array}{c}Z_+\\ Z_+\end{array}\right)\text{and}\psi _{}=\left(\begin{array}{c}Z_{}\\ Z_{}\end{array}\right).$$
(22)
### III.3 Green function

The spectral symmetry of the time evolution operator for network models with chiral symmetry is also shown clearly by considering the Green function, defined for complex energy $`z`$ by

$$G(z)=\frac{1}{1zS}.$$
(23)
The combined consequences of the unitarity and chiral symmetry of $`S`$ together fix the form of $`G(z)`$, which is particularly simple at $`z=1`$. For a discussion of this, it is again convenient to make the rotation $`R`$. We set $`\stackrel{~}{S}=RSR^1`$ and $`\stackrel{~}{G}(z)=RG(z)R^1`$, so that $`\mathrm{\Sigma }_z\stackrel{~}{S}\mathrm{\Sigma }_z=\stackrel{~}{S}^1`$. Unitarity of $`\stackrel{~}{S}`$ implies for $`\stackrel{~}{G}(z)`$

$$\stackrel{~}{G}^{}(z^{})+\stackrel{~}{G}(z^1)=1,$$
(24)
while chiral symmetry is the requirement that

$$\mathrm{\Sigma }_z\stackrel{~}{G}(z)\mathrm{\Sigma }_z+\stackrel{~}{G}(z^1)=1.$$
(25)
At $`z=1`$, these two equations reduce to $`\stackrel{~}{G}^{}(1)+\stackrel{~}{G}(1)=1`$ and $`\mathrm{\Sigma }_z\stackrel{~}{G}(1)\mathrm{\Sigma }_z+\stackrel{~}{G}(1)=1`$. Their solution fixes the form for $`\stackrel{~}{G}(1)`$ to be

$$\stackrel{~}{G}(1)=\frac{1}{2}\left(\begin{array}{cc}1& D\\ D^{}& 1\end{array}\right),$$
(26)
where $`D`$ is a complex matrix.

Our aim, then, is to determine $`\stackrel{~}{G}(1)`$ explicitly for the models we have introduced. We start from Eq. (23), which we write as

$$G(z)=zSG(z)+1.$$
(27)
Following the step from Eq. (16) to Eq. (17), using the matrix $`K`$ in place of $`S`$, and defining

$$K(z)=\left(\begin{array}{cc}1& 0\\ 0& z^1\end{array}\right)K\left(\begin{array}{cc}z& 0\\ 0& 1\end{array}\right),$$
(28)
we obtain

$$G(z)=\frac{1}{2}\left\{1+\frac{1+K(z)}{1K(z)}\mathrm{\Sigma }_z\right\}.$$
(29)
This result simplifies further at $`z=1`$ and after rotation by $`R`$

$$\stackrel{~}{G}(1)=\frac{1}{2}\left(\begin{array}{cc}1& \frac{1+Ue^b}{1Ue^b}\\ \frac{1+Ue^b}{1Ue^b}& 1\end{array}\right).$$
(30)
Note that, since

$$\left[\frac{1+Ue^b}{1Ue^b}\right]^{}=\left[\frac{1+Ue^b}{1Ue^b}\right],$$
(31)
it is clear the Green function has the form expected from Eq. (26).

Finally, we remark that similar manipulations may be used to construct a Hermitian operator which has the same eigenvectors as $`\stackrel{~}{S}`$ and has the block structure

$$\frac{1}{2i}\frac{1\stackrel{~}{S}}{1+\stackrel{~}{S}}=\frac{1}{2i}\left(\begin{array}{cc}0& \frac{1Ue^b}{1+Ue^b}\\ \frac{1Ue^b}{1+Ue^b}& 0\end{array}\right).$$
(32)
### III.4 Continuum limit

The network model representing the quantum Hall plateau transition is knownHo1996  to have a continuum limit in which the time evolution operator is that for a Dirac Hamiltonian in two space dimensions, with random vector and scalar potentials and a uniform mass which varies through zero on crossing the transition. In this subsection we investigate the continuum limit for network models with chiral symmetry. As one might expect, this enables us to make links with previous work on the Dirac equation with randomness in the classes AIIILudwig1994 ; Nersesyan1995  and BDI.Hatsugai1997  A similar approach should also be of interest for the CII network model, but we do not explore it here.

We consider first the model for the class AIII. We obtain the continuum limit by taking all link phases and all scattering amplitudes within the chiral space to be small,footnote2  so that in Eq. (5) $`|\varphi |1`$ and $`|b|1`$, and by taking the node parameter close to its self-dual value, so that in Eq. (3) $`\alpha =\pi /4+m/2`$, with $`|m|1`$.

Before going further, it is necessary to deal with a technical point which arises because of the bipartite structure of the models, discussed above in Sec. III.1. This structure has its origin in the form of $`U`$ displayed in Eq. (11), but is obscured in $`S`$ because of the nature of the off-diagonal blocks in Eq. (7). For that reason, it is useful to make the rotation represented by

$$P=\left(\begin{array}{cc}U^{}& 0\\ 0& 1\end{array}\right),$$
(33)
so that the evolution operator is

$$\widehat{S}PSP^1=\left(\begin{array}{cc}\mathrm{cos}(\beta )& \mathrm{sin}(\beta )\\ \mathrm{sin}(\beta )& \mathrm{cos}(\beta )\end{array}\right)\left(\begin{array}{cc}U& 0\\ 0& U^{}\end{array}\right).$$
(34)
In this basis, the bipartite structure of $`U`$ is clearly also a feature of $`\widehat{S}`$. The price paid is that the chiral symmetry condition, involving $`\widehat{\mathrm{\Sigma }}_xP\mathrm{\Sigma }_xP^1`$, becomes disorder-dependent (though not at lowest order in $`\beta `$), and for this reason we do not use the basis elsewhere in the paper. The effect of the rotation can be represented by moving some of the points on the network at which wavefunction amplitudes are defined, as illustrated in Fig. 5.

In this basis $`\widehat{S}^2`$, like $`U^2`$, is block-diagonal, with separate blocks acting within each of the two sublattices. We obtain a Dirac Hamiltonian $`\widehat{H}`$ by interpreting the part of $`\widehat{S}^2`$ that acts within one sublattice as $`\mathrm{exp}(i\widehat{H})`$, and expanding to lowest order in $`\varphi _l`$, $`b_i\beta _i`$ and $`m`$. Using the numbering of links around a plaquette specified in Fig. 5, we define

$`V=(\varphi _1+\varphi _2+\varphi _3+\varphi _4)/2,`$
$`A_x=(\varphi _4\varphi _2)/2\mathrm{and}A_y=(\varphi _1\varphi _3)/2`$ (35)
as in Ref. Ho1996, . We also define

$`W=(\beta _1+\beta _2+\beta _3+\beta _4)/2,`$
$`B_x=(\beta _1\beta _3)/2\mathrm{and}B_y=(\beta _4\beta _2)/2.`$ (36)
Finally, to put the result in a standard form, we obtain $`H`$ from $`\widehat{H}`$ using a rotation which takes $`(\sigma ^x,\sigma ^y,\sigma ^z)`$ in $`\widehat{H}`$ to $`(\sigma ^y,\sigma ^z,\sigma ^x)`$ in $`H`$. We have finally

$`H=`$ $`\left[(i_xA_x)\sigma ^x+(i_yA_y)\sigma ^y+m\sigma ^z+V\right]\mathrm{\Sigma }_z`$ (37)
$`+\left[B_x\sigma ^x+B_y\sigma ^y+W\right]\mathrm{\Sigma }_y.`$
As expected, the Hamiltonian $`H`$ fulfills the chiral symmetry condition, $`\mathrm{\Sigma }_xH\mathrm{\Sigma }_x=H`$. It can be cast in the form of Eq. (12), with

$$𝐡=(i_x𝐀_x)\sigma ^x+(i_y𝐀_y)\sigma ^y+m\sigma ^z+𝐕.$$
(38)
Here, $`𝐀_x=A_x+iB_x`$, $`𝐀_y=A_y+iB_y`$ and $`𝐕=V+iW`$ are, respectively, complex random vector and scalar potentials. If disorder is also introduced in the mass $`m`$ and as an additional term $`B_z\sigma ^z\mathrm{\Sigma }_y`$, $`H`$ is the most general form for a four-component Dirac Hamiltonian with randomness in symmetry class AIII Guruswamy2000 .

Our network model for class BDI can be treated in a similar way, but for there to be a continuum limit it is necessary to restrict link phases to the value zero, instead of allowing either $`\varphi _l=0`$ or $`\varphi _l=\pi `$. We again obtain $`H`$ in the form of Eq. (12), but with

$$𝐡=(i_xiB_x)\sigma ^x+(i_yiB_y)\sigma ^y+m\sigma ^z+iW,$$
(39)
which is a Dirac Hamiltonian with purely imaginary random vector and scalar potentials and a real mass $`m`$. The same form is obtained from the continuum limit of the Hatsugai-Wen-Kohmoto model,Hatsugai1997  which is known to belong to the symmetry class BDI.

### III.5 Behaviour for special parameter choices

The behaviour of the models we have introduced can be simplified by making one of several alternative special choices for the parameters $`\beta _l`$ that characterise scattering within the two-component space associated with chiral symmetry, or for the parameter $`\alpha `$, specifying scattering at nodes. We discuss these choices in this subsection, restricting ourselves to the model from class AIII, as illustrated in Fig. 2.

#### III.5.1 No chiral scattering

On setting $`\beta _l=0`$ for every link $`l`$, one has two decoupled U(1) network models. With the choices for link phases $`\varphi _l`$ and node parameters $`\alpha `$ that are set out in Sec. II, these models are in localised phases for $`\alpha \pi /4`$, with a localisation length which diverges as $`\alpha `$ approaches $`\pi /4`$.

#### III.5.2 Uniform chiral scattering

There are two alternative and natural ways to pick the parameters $`\beta _l`$ or $`b_l`$ without disorder.

One is to set $`\beta _l=\beta _0`$, a constant. This results in a model that is statistically invariant under $`\pi /2`$ rotations of the lattice. Its time evolution operator develops a gap in its eigenphase spectrum near $`E=0`$ for $`\beta _00`$. To show this, we note from Eq. (7) that the eigenphases $`E_\alpha `$ of $`S`$ can be written for constant $`\beta _l`$ in terms of $`\beta _0`$ and those of $`U`$, which we denote by $`e_\delta `$. We find

$$\mathrm{cos}(E_\alpha )=\mathrm{cos}(\beta _0)\mathrm{cos}(e_\delta ),$$
(40)
so that there are no eigenphases $`E_\alpha `$ in the interval $`[\beta _0,+\beta _0]`$. Because of the gap, this choice for the $`\beta _l`$ is of limited interest.

An alternative choice is to take the hyperbolic angles constant, with $`b_l=b`$ in the transfer matrix. This breaks the spatial symmetry of the model, because after a $`\pi /2`$ rotation of the lattice, the sign of $`b_l`$ is changed in the transfer matrix for one half of the links of the system. For the transfer matrix before such a rotation, it is clear from Eq. (14) that the Lyapunov exponents of the model with chiral symmetry, which we denote by $`\nu _n`$, are related to those of the U(1) model, $`\mu _m`$, by

$$\nu _n=\mu _m\pm b.$$
(41)
We take the fact that the Lyapunov exponents of the two models are related in such a simple way as an indication that this choice for $`b_l`$ does not lead to behaviour generic for the chiral symmetry class. The result is nevertheless interesting, because small fluctuations in $`b_l`$ about an average value $`b`$ are likely to restore generic properties of the symmetry class without large changes in the values of the Lyapunov exponents. In particular, for suitable values of $`b`$ and the node parameter $`\alpha `$, the spectrum for $`\nu _n`$ expected from Eq. (41) has a gap around $`\nu =0`$, indicating that a localised phase is possible for the network model with chiral symmetry.

#### III.5.3 Strong chiral scattering

The choice $`|\beta _l|=\pi /2`$ for all $`l`$ breaks the system up into a set of disconnected loops, as illustrated in Fig. 6. In this limit, possible values for the eigenphases of $`S`$ are $`E=0,\pi ,\pm \pi /2,,\pm \alpha ,\pm \alpha +\pi `$ and $`\pm \alpha \pm \pi /2`$, where $`\alpha `$ is the scattering angle at nodes. The probability for each eigenvalue to occur depends on the distribution chosen for the signs of $`\beta _l`$. To determine the nature of eigenfunctions as this limit is approached, degenerate perturbation theory in $`|\beta |\pi /2`$ is required. At leading order, this generates a tight-binding model with nearest neighbour hopping and site dilution on a square lattice.

#### III.5.4 Decoupled plaquettes

Setting the node parameter to $`\alpha =0`$ or $`\alpha =\pi /2`$, the system breaks up into decoupled plaquettes. Eigenphases for one such plaquette can be calculated analytically, and vary continuously with link phases $`\varphi _l`$. Eigenstates are therefore non-degenerate and localised in this limit. It is also interesting to consider behaviour at the edge of the system in a localised phase. For a single copy of the U(1) network model which is defined on the plane rather than a torus, two distinct localised phases arise in the two limits for $`\alpha `$: one with current-carrying edge states and one without. Following the construction described in Sec. II, for models with chiral symmetry two oppositely propagating edge states appear, with scattering between them for $`\beta _l0`$. Viewed as a strictly one-dimensional system, these edge states are expected to be critical if the average $`\beta _l`$ is zero, and localised otherwise. Weak coupling between plaquettes, however, is likely to localise the edge states for all values of $`\beta _l`$. We therefore do not expect an analogue of the quantum Hall plateau transition in the models we study, a conclusion that one may also anticipate from the fact that the models are statistically parity-symmetric.

### III.6 Models with only vector potential disorder

There has been very extensive interest in Dirac Hamiltonians on two space dimensions, with only vector potential disorder. In this subsection we show how lattice versions of such problems can be realised, using the network model. A similar treatment for tight-binding Hamiltonians has been described recently by Motrunich et al.Motrunich2002

Our starting point is a special choice for the chiral scattering angles $`\beta _l`$, derived from a potential function as follows. The potential is defined at the points which are marked with filled circles on the lattice shown in Fig. 2. Suppose that, on the network drawn with (say) full lines in this figure, the link $`l`$ runs from point $`i`$ to point $`j`$, and let the potential take the values $`\mathrm{\Phi }_i`$ and $`\mathrm{\Phi }_j`$ at these points. The hyperbolic angle $`b`$ is then the gradient of the potential, with

$$b_l=\mathrm{\Phi }_j\mathrm{\Phi }_i.$$
(42)
Introducing a diagonal matrix $`\mathrm{\Phi }`$, with $`\mathrm{\Phi }_l`$ as entries, one finds that the evolution operator $`𝒰(b)`$ appearing in Eq. (20) is a similarity transform of the operator $`U`$

$$𝒰(b)=e^\mathrm{\Phi }Ue^\mathrm{\Phi }.$$
(43)
Likewise, one has

$$𝒰(b)=e^\mathrm{\Phi }Ue^\mathrm{\Phi }.$$
(44)
These facts enables us to determine the stationary states for a system with non-zero $`\mathrm{\Phi }`$ in terms of those for a system with $`\mathrm{\Phi }=0`$, if boundary conditions or fine-tuning of the disorder realisation allow such states. Specifically, suppose that the system without chiral scattering has a stationary state $`\psi _0`$, so that

$$U\psi _0=\psi _0.$$
(45)
Then the system with chiral scattering has stationary states that, in the notation of Eq. (21), are given by

$$Z_+=e^\mathrm{\Phi }\psi _0\text{and}Z_{}=e^\mathrm{\Phi }\psi _0.$$
(46)
A similar simplification applies to the Green function at $`z=1`$. The relevant part is the off-diagonal block appearing in Eq. (30), for which one has

$$\frac{1+Ue^b}{1Ue^b}=e^\mathrm{\Phi }\frac{1+U}{1U}e^\mathrm{\Phi }.$$
(47)
Also the Hamiltonian defined by Eq.(32) shows a decomposition reminiscent of the continuum version

$`_\mathrm{\Phi }`$ $`=`$ $`\left(\begin{array}{cc}e^\mathrm{\Phi }& 0\\ 0& e^\mathrm{\Phi }\end{array}\right)\left(\begin{array}{cc}0& _0\\ _0& 0\end{array}\right)\left(\begin{array}{cc}e^\mathrm{\Phi }& 0\\ 0& e^\mathrm{\Phi }\end{array}\right),`$ (54)
where

$$_0=\frac{1}{2i}\frac{1U}{1+U}.$$
(55)
There are two separate applications of these ideas which are of interest. The first is to take the operator $`U`$ to be without disorder, setting all link phases $`\varphi _l=0`$. Choosing in addition the node parameter to have its critical value $`\alpha =\pi /4`$, $`U`$ has a stationary state $`\psi _0`$ which is constant in space. Then the wavefunctions $`Z_+`$ and $`Z_{}`$ are lattice versions of the eigenstates known for the Dirac equation with random vector potential, whose statistical properties have been very extensively studied. This is supported by the fact that in this case the continuum limit of the model Eq. (37) simplifies to a pair of Dirac fermions with real vector potential, equivalent to the ones described in Ref. Motrunich2002, .

A second application is to include disorder both via $`\mathrm{\Phi }`$ and in $`U`$ itself. Taking all link phases to be random, and again setting $`\alpha =\pi /4`$, $`\psi _0`$ is a critical wavefunction for the quantum Hall plateau transition. The combinations $`Z_\pm =e^{\pm \mathrm{\Phi }}\psi _0`$ will therefore have multifractal fluctuations that are determined from an interplay of the distributions of $`\psi _0`$ and of $`\mathrm{\Phi }`$. It seems very likely that this will result in a new type of critical behaviour, but we leave details for a future investigation.

## IV Numerical results for the two-dimensional AIII model

### IV.1 Introduction

In this section we present results from a numerical investigation of the two-dimensional AIII network model, as defined in Fig. 2 and Eq. (7). Most of our results are for systems in which the hyperbolic angles $`b_l`$ have a Gaussian distribution of zero mean and variance $`g`$. For these systems, we examine behaviour as a function of the node parameter $`\alpha `$ and $`g`$, parameterising the latter with $`\gamma `$, defined by $`\mathrm{sin}(\gamma )=\mathrm{tanh}(\sqrt{g})`$. Hence we consider the parameter space spanned by $`0\alpha \pi /4`$ and $`0\gamma \pi /2`$. We also present some data for systems in which the angles $`\beta _l`$ are uniformly distributed between $`0`$ and $`2\pi `$. To set the scene, we first summarise what can be anticipated for behaviour of the model, from a consideration of limiting values for the parameters $`\alpha `$ and $`\gamma `$ and from existing analytical results.

We start by recalling that at $`\gamma =0`$ the system consists of two decoupled copies of the U(1) network model. States in each copy are localised for all $`\alpha `$ in the range $`0\alpha <\pi /4`$, with a localisation length that diverges at $`\alpha =\pi /4`$. One expects the localised phase to be stable against small perturbations, so that it will extend over a region in parameter space with $`0\alpha <\pi /4`$ and $`\gamma `$ small. Conversely, critical states in the two U(1) models at $`\alpha =\pi /4`$ are easily mixed by non-zero chiral coupling $`b_l`$, so that it is reasonable to guess that the critical phase predicted by Gade and WegnerGade1991 ; Gade1993  will be found in the region with $`\alpha `$ close to $`\pi /4`$ and $`\gamma >0`$.

Properties along the (critical) line $`\alpha =\pi /4`$ follow from studies of the non-linear sigma model:Gade1991 ; Gade1993 ; Altland1999  the renormalisation group equations

$$\frac{\text{d}\lambda }{\text{d}l}=0\frac{\text{d}\lambda _A}{\text{d}l}=\frac{\lambda ^2}{16},$$
(56)
are arguedGuruswamy2000  to be exact to all order of perturbation theory. Note that when $`\alpha `$ is close to but not exactly equal to $`\pi /4`$, a non-zero mean mass for the underlying random Dirac fermions is generated. In that case, Eq. (56) should be extended to incorporate it, which is beyond the scope of the present work. This non-zero mass leads to localisation when large enough on a scale set by $`\gamma `$.

The coupling constant $`\lambda `$ parameterises the conductivity $`\sigma `$, with $`\sigma =1/(8\pi \lambda ^2)`$. Since $`\lambda `$ does not flow under renormalisation, the critical phase is described by a line of fixed points and may have a range of values for its conductivity. Flow of the second coupling constant, $`\lambda _A`$, determines properties at non-zero energy, giving in particular for the density of states $`\rho (E)`$ at small energy $`E`$ the scaling lawGade1991 ; Gade1993

$$\rho (E)\frac{E_0}{E}\mathrm{exp}\left(\kappa \sqrt{\mathrm{ln}\left(E_0/E\right)}\right),$$
(57)
where $`E_0`$ is a microscopic energy scale. Using the results of Ref. Guruswamy2000, , $`\kappa =4\sqrt{2}/\lambda ^2`$. As emphasised in Ref. Motrunich2002,  such behaviour for $`\rho (E)`$ sets in only below a crossover energy scale $`E_c`$: from Ref. Guruswamy2000, , we estimate this scale to be given by

$$E_c=E_0\mathrm{exp}\left(\alpha \sigma ^2\right),$$
(58)
where $`\alpha =(32\sqrt{2}\pi )^2`$. We note also that Motrunich et al have argued that the power law in the exponent of Eq. (57) is modified at very small energy scales,Motrunich2002  and that their results have been reproduced within a field-theoretic approach.Mudry2002

Since the derivations of these results depend on a number of technical assumptions, some checks are very desirable. Unfortunately, computational work is challenging because it is hard to reach scales smaller than $`E_c`$, and at intermediate energies the same analysis predicts an approximately power-law dependence of the density of states on energy, with an effective exponent which depends on disorder strength.

In our numerical work testing these expectations, we use calculations of two types, both applied to quasi-one dimensional systems of width $`M`$ plaquettes and length $`L`$ plaquettes (so that the time evolution operator $`S`$ appearing in Eq. (7) is a $`4ML\times 4ML`$ matrix, and the transfer matrix $`T`$ appearing in Eq. (14) is a $`4M\times 4M`$ matrix). For such samples, we find the localisation length $`\xi _M`$ as the inverse of the smallest positive Lyapunov exponent for the transfer matrix, calculated using the standard approach.Pichard1981  Applied directly to the transfer matrix as formulated in Sec. II, the resulting localisation length is that for eigenstates of the time evolution operator at eigenphase zero; by associating an extra phase $`e^{iE}`$ with every link (breaking chiral symmetry), we are also able to find the localisation length for states with non-zero eigenphase. Separately, we determine the density of states for eigenphases of the time evolution operator $`S`$, applying a version of the recursion method,MacKinnon1981  adapted to treat a unitary operator in place of the Hermitian one for which the approach was originally formulated. We use the terms eigenphase and energy interchangeably in describing our results.

### IV.2 Phase diagram

We determine a phase diagram for the model in the $`\alpha `$, $`\gamma `$ plane from a study of the zero-energy localisation length $`\xi _M`$ in quasi-one dimensional geometry, and its dependence on $`M`$. (We do not apply this approach to models with uniformly distributed $`\beta _l`$, since our transfer matrix calculations develop numerical instabilities if $`|\mathrm{cos}(\beta _l)|1`$.) The quantity of central interest is the ratio $`\xi _M/M`$. In a localised phase, $`\xi _M`$ approaches a finite limit $`\xi `$, the bulk localisation length, for large $`M`$, and so the ratio $`\xi _M/M`$ decreases with increasing $`M`$, varying as $`\xi /M`$ for sufficiently large $`M`$. By contrast, in a critical phase, the ratio approaches for large $`M`$ a finite constant whose value can be identified with the conductivity of the model.Chalker1993

Representative data for $`\xi _M/M`$ at $`M=8,16,32`$ and $`64`$ are shown as a function of $`\gamma `$, in Fig. 7 for $`\alpha =7\pi /32`$, and in Fig. 8 for $`\alpha =\pi /8`$.

In both cases, two regimes of behaviour are observed: for small $`\gamma `$, $`\xi _M/M`$ decreases with increasing $`M`$, indicating a localised phase, while for $`\gamma `$ close to $`\pi /2`$, $`\xi _M/M`$ is independent of $`M`$ for the system widths studied, suggesting a critical phase. By identifying the value of $`\gamma `$ which divides the two regimes and studying its dependence on $`\alpha `$, we arrive at the phase diagram displayed in Fig.9.

Properties in the critical phase are illustrated in more detail by the behaviour of the ratio $`\xi _M/M`$ as a function of $`\gamma `$ with $`\alpha =\pi /4`$, shown in Fig. 10. On this line in parameter space, behaviour at all points is critical rather than localised, and the data show a finite limiting value $`\sigma `$ for the ratio at large $`M`$, with $`\sigma `$ dependent on $`\gamma `$: such behaviour is expected from Eq. (56) and the identification of the coupling constant $`\lambda `$ as being a function only of conductivity, which is given in turn by $`\xi _M/M`$. This figure also illustrates crossover of behaviour with increasing $`M`$ at small $`\gamma `$, from that of the U(1) model to that of the AIII model. The crossover is responsible for discontinuous variation of $`\sigma `$ with $`\gamma `$ at $`\gamma =0`$.

We show in Fig. 11 the energy dependence of the localisation length at a point in the localised phase (left panel) and at a point in the critical phase (right panel), with data at $`L=2.5\times 10^5`$ for a range of system widths $`M`$, as well as an extrapolation to the two-dimensional limit. It is clear from the data that the localisation length in the critical phase decreases very rapidly with $`E`$ away from $`E=0`$.

### IV.3 Density of states

Our interest in the (average) eigenphase density of the time evolution operator is focussed mainly on behaviour near zero energy, where we expect in the critical phase singular behaviour approaching that of Eq. (57), and in the localised phase either a density vanishing quadratically with energy, as for the random matrix ensemble in this symmetry class,Zirnbauer1996  or singularities arising from Griffiths strings, as discussed in Ref. Motrunich2002, . Before examining the small energy region in detail, we give an overview of behaviour for the entire range of eigenphases, shown in Fig. 12 for $`\alpha =\pi /4`$ and uniformly distributed $`\beta _l`$. (Similar results are obtained with Gaussian $`b_l`$ at $`g1`$.) From the discussion of Sec. III.1, we need to display data only over the range $`0E\pi /2`$. These results were obtained for a system size of $`M=16`$ and $`L=10^5`$, which is large enough for $`\rho (E)`$ to be self-averaging at the scale used here.

A divergence of $`\rho (E)`$ as $`E0`$ is apparent in Fig. 12, which we now investigate in more detail. In Fig.13 we show $`\rho (E)`$ as a function of $`E`$, with logarithmic scales for both axes, for a sequence of models, with $`\alpha =\pi /4`$ and $`b_l`$ Gaussian distributed in each case. Members of the sequence are chosen to have successively smaller values of $`\sigma `$, and hence larger values for the energy scale $`E_c`$, since by increasing $`\gamma `$ we decrease $`\sigma `$, as is evident from the inset to Fig. 7. The data were obtained using systems of width $`M=16`$, and lengths in the range $`L=10^5`$ to $`L=10^9`$, in order to achieve satisfactory self-averaging. Over the energy range accessible, the variation of $`\rho (E)`$ with $`E`$ is approximately power-law, and can be characterised following Ref. Motrunich2002,  using an effective dynamical exponent $`z`$, with

$$\rho (E)|E|^{1+2/z}.$$
(59)
As one approaches the asymptotic low-energy behaviour, the effective exponent $`z`$ increases, ultimately to infinity: the large values of $`z`$ reached here indicate close approach to the limiting low-energy behaviour.

By repeating such calculations for a range of $`M`$ values, we have checked that these results are not influenced by finite system width system. This is shown in Fig. 14, in which the behaviour with uniformly distributed $`\beta _l`$ is also illustrated.

The rather rapid convergence to the large $`M`$ limit is a consequence of the small values of the localisation length for energies away from $`E=0`$, illustrated above.

Finally, we examine the evolution of the density of states as one moves from the critical phase into the localised phase. Results are presented in Fig. 15. At the broad scale to which these calculations are restricted, behaviour in the localised phase is consistent with the power law of Eq. (59), with a power $`1+2/z`$ which, deep in the localised phase, is positive and increases as the localisation length decreases. Such behaviour was proposed as generic for localised systems with chiral symmetry in Ref.Motrunich2002, .

## V Summary

In summary, we have introduced network models which are realisations of the chiral symmetry classes. We have argued that they are interesting from several points of view. They are useful as a starting point for numerical studies of these symmetry classes, offering access in the two-dimensional case we have examined to both critical and localised phases, and displaying band-centre singularities in the critical phase which approach quite closely the expected asymptotic form. The models also have a striking connection to network models without chiral symmetry, but with absorption and amplification. Moreover, by imposing constraints on the disorder, they serve as lattice versions of problems with randomness entering only through a vector potential. It seems likely that disorder of this kind can generate types of critical behaviour different from those known previously for localisation problems in two dimensions.

## Acknowledgements

This work was supported in part by the EPSRC under grant GR/J78327. We thank N. Regnault for use of computer facilities.