Estimate of the critical exponent of the Anderson transition in the three and four dimensional unitary universality classes

Disordered non-interacting systems are classified into ten symmetry classes, with the unitary class being the most fundamental. The three and four dimensional unitary universality classes are attracting renewed interest because of their relation to three dimensional Weyl semi-metals and four dimensional topological insulators. Determining the critical exponent of the correlation/localistion length for the Anderson transition in these classes is important both theoretically and experimentally. Using the transfer matrix technique, we report numerical estimations of the critical exponent in a U(1) model in three and four dimensions.


Introduction
More than half a century after its discovery, 1) the Anderson transition continues to attract attention. The recent experimental realisation 2) in a cold atom system of the quantum kicked rotor (a system which has been mapped to the Anderson localisation problem 3,4) ) has opened a new avenue for research of Anderson localisation and the Anderson transition. The measurement 5) of the critical exponent of the correlation/localisation length for the Anderson transition in this experiment is in very good agreement with high accuracy numerical analysis for the three dimensional (3D) orthogonal universality class. 6,7) Disordered non-interacting systems are classified into ten symmetry classes. 8) The same symmetry classification 9,10) is also applicable to topological insulators. 11,12) Amongst these ten symmetry classes, the unitary class (class A) is the most general, where no symmetries such as time reversal, spin-rotation, chiral or particle-hole symmetries are present. Systems in magnetic fields and systems with magnetic impurities, for example, belong to the unitary class. The critical behavior of the Anderson transition in 3D systems subject to uniform mag- * slevin@phys.sci.osaka-u.ac.jp netic fields was recently studied with high precision. 13) Studies of disordered 3D systems subject to random magnetic fields 14) and with classical magnetic impurities 15) have also been reported.
In the quantum kicked rotor system, the effective dimensionality is determined by the number of incommensurate frequencies used to modulate the amplitude of the periodic kick. 16) The 3D orthogonal universality class has already been realised experimentally in this system. 2) The addition of modulation by an extra incommensurate frequency would realise the 4D orthogonal universality class. In addition, a version of the quantum kicked rotor in which the quantum Hall effect might be realised has also been proposed. 17) Thus, it is plausible that the 4D unitary universality class might be realised in a suitable quantum kicked rotor system. This universality class is also of interest from the view point of three dimensional Weyl semi-metals. 18) Just as two dimensional Dirac electrons appear on the surface of 3D topological insulators, three dimensional Weyl electrons appear on the "surface" of 4D topological insulators belonging to the unitary symmetry class.
On the other hand, the Anderson transition in the 2D symplectic class and the quantum spin Hall transition show the same critical exponent. 24) Whether the critical exponent of the 4D quantum Hall transition is the same as that of the Anderson transition in the 4D unitary class is an open problem.
In this paper, we study the 3D and 4D unitary universality classes and estimate the critical exponents with high precision by studying numerically the U(1) model. In the following section, we introduce the model, and the method is explained in Section 3. The results of the finite size scaling analyses are shown in Section 4, followed by the discussion.

Model
The U(1) model is a variant of Anderson's model of localisation 1) in which time reversal symmetry is broken by multiplying the unit hopping elements of that model by complex phases. For the purpose of estimating critical exponents it is helpful that the length scale associated with the breaking of time reversal symmetry is as short as possible. This is achieved by using completely random phases. 14) The Hamiltonian of the U(1) model is Here, |i is an electron orbital centered on site i of a d-dimensional simple cubic lattice and the first sum is over all the sites of this lattice. The lattice constant (taken as unity) sets the unit of length. The second sum is over pairs of nearest neighbour sites on this lattice.
The magnitude of the hopping elements (taken as unity) between nearest neighbour sites on the lattice sets the unit of energy. The orbital energies E i are independently and identically distributed according to the distribution where The scale of the fluctuations of the orbital energies is set by the parameter W.
As already mentioned the distribution of phases is taken as completely random, i.e. the φ i j with i < j are independently and identically distributed uniformly on [0, 2π]. To ensure the Hamiltonian is Hermitian we set In

Transfer matrix method
The brief description of the transfer matrix method that we give here follows closely Ref.
7. We refer the reader to that reference for further explanation and for details that are omitted.
In the transfer matrix method, the transmission of an electron with an arbitrary energy E through a very long disordered wire is considered. We denote the length of the wire by L x and consider cross sections L × L in the 3D simulation and L × L × L in the 4D simulation.
In practice, L x is many orders of magnitude larger than L. We impose periodic boundary conditions on the wavefunction in the transverse directions. Starting from the time independent Schrödinger equation with energy E we derive a transfer matrix product For sufficiently long disordered wires, the transmission probability decreases exponentially with the length of the wire. 26,27) The associated exponential decay length is equal to the reciprocal of the smallest positive Lyapunov exponent γ of the matrix product Eq. (5). The Lyapunov exponents are the eigenvalues of the matrix and are estimated using the procedure described in detail in Ref. 7.

3D
We set the energy E at the band centre, i.e. E = 0 and simulated the disorder range versus the disorder W.

4D
We set the energy E at the band centre and simulated the disorder range 32 ≤ W ≤ 42 (except for the largest two system sizes where we used a narrower range primarily because  The error in the data is less than the symbol size so we omit error bars. We also show the finite size scaling fit (sold lines).

Finite size scaling analysis
In our simulations we varied the disorder W while keeping the energy E fixed, so there is a critical disorder W c , which is a function of energy, that separates the localised and extended phases. At this critical disorder the localisation (correlation) length ξ has a power law that is described by a critical exponent ν. We performed finite size scaling analyses to estimate the critical exponent ν, the critical disorder W c as well as other quantities (described below).
When corrections due to irrelevant variables could be neglected, either because the numerical data had insufficiently high precision to resolve such corrections or because small system sizes had been excluded, we fitted the data to the equation When corrections due to irrelevant variables could not be neglected, we fitted the data to the equation Here, φ 1 and φ 2 are scaling variables and The first of these variables φ 1 is a relevant scaling variable whose size dependence is related to the critical exponent ν The second of these variables is an irrelevant scaling variable that permits corrections to scaling due to irrelevant scaling variables to be taken into account in an approximate manner.
The size dependence is described by an irrelevant exponent The scaling variables are expanded in powers of w which allows us to take account of possible nonlinearity of the scaling variables in the disorder W. For the relevant scaling variable we must have so we fix b 1,0 = 0. The scaling function is expanded in powers of the scaling variables To avoid ambiguity in the definition of the fit we fix The constant term, the value of which is expected to be universal, is denoted The orders of the expansions are defined by 4 integers m 1 , m 2 , n 1 and n 2 . The quality of the fit to the data is assessed using the χ 2 statistic and by calculating the goodness of fit probability.
We systematically performed fits for various orders of the expansions. After rejecting fits for which the goodness of fit was too small, typically Q ≪ 0.1, we chose the fit with smallest number of parameters for which the estimation of the parameters was stable against increase in the orders of the expansions. In each case, to determine the precision of the estimates of the fitted parameters we generated 400 synthetic data sets and determined 95% confidence intervals from the fits to these synthetic data sets. We refer the reader to Sec.
and then subtract this from the data point We then plot the data versus φ 1 . In addition, we plot the curve  Table I. The range of data, the orders of the expansions, the total number of data N D , the number of parameters N P , the value of χ 2 min obtained for the best fit, and the corresponding goodness of fit probability Q.

3D
We show the finite size scaling fit in Fig. 1 and give the details in Tables I and II. We found that it was not possible to fit the full data set without including a correction due to

4D
We show the finite size scaling fit in Fig. 2 and give the details in Tables I and II. We found that it was not necessary to include corrections to scaling due to irrelevant variables.
Again we also performed fits on a narrower disorder range, and also with smaller system sizes excluded. In both cases we found results that were entirely consistent with the fit of the full data set. We present the demonstration of single parameter scaling in Fig. 4.

Discussion
For 3D, in a previous study 28) we found ν = 1.43 ± .04 for a system in a uniform magnetic field. The estimate for the critical exponent that we obtained here is consistent with that and also more precise. In that work, which was also a transfer matrix study of Lyapunov exponents, we also estimated the quantity Λ c which is the inverse of Γ c . Translating the result of our previous work for easier comparison, we find Γ c = 1.760 ± .004. There is a clear discrepancy (≈ 3%) with the estimate obtained here. This is unexpected since this value is expected to be universal, i.e. to depend only on dimensionality and symmetry class and not on the details of the model considered. A possible explanation for this is the smaller system sizes used, and the neglect of corrections to scaling, in our previous work. Another possibility is that Γ c is affected slightly by the anisotropy 29) introduced by the magnetic field. Another explanation would be a violation of universality but we think it unlikely.
We have also applied the scaling method proposed by Harada 30,31) to our data. In Harada's method the scaling function is expressed as a Gaussian process rather than a polynomial.
Since the non-linearity of the scaling variables (higher order terms in Eq. (15)) is neglected in that approach, we restricted the range of disorder as described in Table I. We have con-