Critical Exponent of the Anderson Transition using Massively Parallel Supercomputing
aa r X i v : . [ c ond - m a t . d i s - nn ] A ug Journal of the Physical Society of Japan
Full Paper
Critical Exponent of the Anderson Transition using Massively ParallelSupercomputing
Keith Slevin ∗ and Tomi Ohtsuki Department of Physics, Graduate School of Science, Osaka University, 1-1 Machikaneyama, Toyonaka, Osaka560-0043, Japan Division of Physics, Sophia University, Chiyoda-ku, Tokyo 102-8554, Japan
To date the most precise estimations of the critical exponent for the Anderson transition have been made using thetransfer matrix method. This method involves the simulation of extremely long quasi one-dimensional systems. Themethod is inherently serial and is not well suited to modern massively parallel supercomputers. The obvious alternativeis to simulate a large ensemble of hypercubic systems and average. While this permits taking full advantage of bothOpenMP and MPI on massively parallel supercomputers, a straight forward implementation results in data that doesnot scale. We show that this problem can be avoided by generating random sets of orthogonal initial vectors with anappropriate stationary probability distribution. We have applied this method to the Anderson transition in the three-dimensional orthogonal universality class and been able to increase the largest L × L cross section simulated from L =
24 (New J. Physics, , 015012 (2014)) to L =
64 here. This permits an estimation of the critical exponent withimproved precision and without the necessity of introducing an irrelevant scaling variable. In addition, this approach isbetter suited to simulations with correlated random potentials such as is needed in quantum Hall or cold atom systems.
1. Introduction
The transfer matrix method was first introduced into thefield of Anderson localisation by Pichard and Sarma, andMacKinnon and Kramer.
2, 3)
These papers reported the firstnumerical evidence in favour of the scaling theory of locali-sation. Since that pioneering work the method has been em-ployed extensively and very successfully to estimate withhigh precision the critical exponents for the standard Wigner-Dyson symmetry classes in various dimensions and forthe quantum Hall e ff ect. However, the method is inherently serial and is not wellsuited to modern massively parallel supercomputers. Here, wedescribe an adaptation of the method that is better suited tosuch computers. We show that the critical exponent can be es-timated correctly by simulating an ensemble of cubes ratherthan the single very long quasi-one dimensional systems inthe serial method. The key point is to make the initial ma-trix used in the method a random matrix that is sampled froman appropriate stationary probability distribution. We demon-strate the e ff ectiveness of the method by applying it to esti-mate the critical exponent for the Anderson transition in thethree dimensional orthogonal universality class.
2. Model and Simulation Method
The Hamiltonian for Anderson’s model of localisation may be written in the form H = X i E i | i i h i | − V X h i j i | i i h j | . (1)Here, | i i is a localised orbital on site i of a three dimensionalcubic lattice. The first sum is over all sites on the lattice andthe second sum is over pairs of nearest neighbours. We mea-sure all energies in units of the hopping energy V betweennearest neighbour orbitals (so that in what follows V does not ∗ [email protected] appear explicitly, i.e. V = E i of these orbitalsare assumed to be identically and independently distributedrandom variables with a uniform distribution p ( E i ) = ( / W | E i | ≤ W / , . (2)The parameter W determines the degree of disorder. Since thisHamiltonian commutes with the complex conjugation opera-tor, i.e. a time reversal operator that squares to plus unity, thismodel is in the orthogonal symmetry class. We consider a system with a uniform square cross section L × L which we divide into layers labelled by their x coor-dinate. We then re-write the time independent Schr¨odingerequation for a state vector | Ψ i and energy EH | Ψ i = E | Ψ i , (3)in the following form ψ x + − ψ x ! = M x ψ x − ψ x − ! . (4)Here, ψ x is the vector of wavefunction amplitudes on layer x (in some suitable order)( ψ x ) y , z = h x , y , z | Ψ i . (5) H x is the following sub-matrix of the Hamiltonian( H x ) y , z , y ′ , z ′ = h x , y , z | H | x , y ′ , z ′ (cid:11) , (6)and M x is the following 2 N × N transfer matrix (with N = L ), M x = H x − E N N − N N ! (7)with 1 N and 0 N the N × N unit and zero matrices, respec-tively. The boundary conditions in the transverse directionsmust also be specified. Throughout this paper, we set the
1. Phys. Soc. Jpn.
Full Paper energy at the band centre, i.e. E =
0, and impose periodicboundary conditions in the transverse directions.Since solutions of the time independent Schr¨odinger equa-tion are stationary, the probability flux must be conserved bythe transfer matrix multiplication. This implies that the trans-fer matrix must satisfy the following relation M Tx Σ c M x = Σ c , (8)where Σ c = N − i N i N N ! . (9)Here, i N = i N . We consider a very long quasi-one-dimensional bar oflength L x . The wave-function amplitudes on the first two lay-ers are related to the wave-function amplitudes on the last twolayers as follows ψ L x + − ψ L x ! = M L x · · · M ψ − ψ ! . (10)This involves the product of L x independently and identicallydistributed random matrices M = M L x · · · M . (11)According to the theorem of Oseledec, the following limit-ing matrix exists Ω = lim L x →∞ ln M T M L x . (12)Note the limit depends on the particular sequence of randommatrices not just on the distribution. However, for the eigen-values { γ i } of Ω we obtain the same values for almost allsequences, i.e. with probability one. These values are calledLyapunov exponents. From Eq. (8) it can be shown that theseeigenvalues occur in pairs of opposite sign. It is usual to num-ber them as follows γ > γ > · · · γ N > γ N + = − γ N > · · · > γ N = − γ . (13)To estimate the Lyapunov exponents we start with a 2 N × N orthogonal matrix, truncate the matrix product at a very largebut finite L x , and perform a QR factorization of the result QR = MQ . (14)Here, Q is a 2 N × N orthogonal matrix and R is a 2 N × N upper triangular matrix with positive diagonal elements. Wethen define ˜ γ i = L x ln R i , i . (15)In the limit of infinite length γ i = lim L x →∞ ˜ γ i . (16)For su ffi ciently large L x , the { ˜ γ i } may be used to estimate theLyapunov exponents. In practice, it is not possible to imple-ment the calculation Eq. (14) straightforwardly. Instead, toavoid round o ff error, it is necessary to perform intermediateQR factorizations at regular intervals (say after every q trans-fer matrix multiplications). The details are described in Ref.8. Note that, for the purposes of finite size scaling, it is usual, though not essential, to focus on the smallest positive Lya-punov exponent γ N . The calculation of the negative Lyapunovexponents is then not necessary. In this case, it is su ffi cient tomake Q a 2 N × N matrix with orthogonal columns, and R then becomes an N × N matrix. This saves considerable com-putational time.We next turn to the question of how to estimate the criticalparameters. The critical parameters of main interest are thecritical disorder W c separating the di ff usive (metal) phase for W < W c from the localised (insulator) phase for W > W c ,and the critical exponent ν that describes the divergence ofthe correlation (localisation) length ξ at the critical point ξ ∼ | W − W c | − ν . (17)These, and other critical parameters, are estimated by fittingthe system size and disorder dependence of the dimensionlessquantity Γ = γ N L , (18)to a model derived from the following one-parameter scalinghypothesis h ˜ γ N i L x = f L x ξ , L ξ , L ξ ! , (19)with f a scaling function. For su ffi ciently long systems suchthat the e ff ect of the initial matrix Q becomes negligible, weexpect the right hand side of Eq. (19) to be proportional to L x .Thus, after taking the limit of infinite sample length, Eq. (19)reduces to Γ = γ N L = f Q1D L ξ ! , (20)where an ensemble average or tilde are no longer needed anda related scaling function f Q1D has been introduced.This method requires the simulation of a single very longsample. While this method has been employed very success-fully in numerous simulations over the preceding decades,it is an inherently serial calculation that does not allow usto take advantage of modern parallel computers particularlythose with hundreds of CPU nodes.
The alternative that we consider here is to simulate an en-semble of much shorter samples and consider an ensembleaverage. For simplicity we consider cubes with L x = L . Thescaling hypothesis Eq. (19) then becomes Γ = h ˜ γ N i L = f L ξ , L ξ , L ξ ! = f L ξ ! . (21)It seems that we can now take full advantage of modern com-puters by simulating an ensemble of samples in parallel. How-ever, when attempting to analyse the numerical data we runinto a di ffi culty. In writing Eq. (21) we have neglected thedependence on the initial matrix Q . The importance of thismatrix is immediately made clear by reference to Fig. 1 wherewe show data for the sample mean of ˜ γ N obtained with Q = N N ! . (22)The data do not exhibit a common crossing point in the vicin-ity of the critical disorder W c ≈ .
5. While at first glance
2. Phys. Soc. Jpn.
Full Paper G W L=24L=16L=12
Fig. 1.
Estimates of the quantity
Γ = h ˜ γ N i L for ensembles of L × L × L cubes with L = ,
16 and 24 for disorder W in the range [15 . , . Q given in Eq. (22) has been used. it might appear that the inclusion of an irrelevant correctionmight su ffi ce to restore a common crossing point, this is notthe case; after further contemplation of Fig. 1 we see that therequired correction would have to be relevant not irrelevant.Fortunately, this problem is easily solved. Instead of a fixed Q , we make Q a random matrix that is sampled from a prob-ability distribution that is invariant under convolution with thetransfer matrix distribution, i.e. by generating 2 N × N randommatrices with orthogonal columns with a distribution that isinvariant or stationary under the operation Q ′ R = M x Q . (23)For such a distribution, we see from Eq. (33) of Ref. 8 that˜ γ N becomes a sum of i.i.d. random variables, from which itfollows that the dependence of h ˜ γ N i on the length L x vanishes.In this case h ˜ γ N i = γ N . (24)It immediately follows that f Q1D = f . (25)The next question is how to generate such matrices. Wehave found that the following procedure works well. We startwith Q given by Eq. (22) and calculate Q ′ R = M q · · · M Q . (26)The matrix R is then discarded and we set Q = Q ′ . Thiscalculation is then repeated a total of, say, r times. Note that,for a given sequence of transfer matrices, the result of thiscalculation depends only on the total number of transfer ma-trix multiplications, i.e. on the product of q and r , not on q and r separately. For a given L , we have found that, whena su ffi cient number of randomizing multiplications are per-formed, the distribution of ˜ γ N becomes independent of thenumber of such multiplications. We judge this by applyingthe Kolmogorov-Smirnov test to the resulting data for ˜ γ N with di ff erent numbers of randomizing multiplications. Forsu ffi ciently large number of randomizing multiplications we G W Fig. 2.
A verification of Eq. (24) for L =
12. Estimates of
Γ = h ˜ γ N i L ob-tained using the parallel method (points without error bars) are comparedwith estimates of Γ = γ N L obtained using the serial method (error bars with-out points). For both sets of data the precision is approximately 0 . -0.1 0.0 0.1 0.2 0.3 0.4 N Fig. 3.
Histograms for L =
24 and W = . γ N obtained with random Q (no shading) generated using 64 transfer matrixmultiplications and with fixed Q (shaded) given by Eq. (22). The ensemblesize is 589824. find that the Kolmogorov-Smirnov test is unable to distin-guish the distribution of ˜ γ N obtained. When performing theKolmogorov-Smirnov test we think it is important that the en-semble size used for the test match that of the ensemble sizeused to accumulate data for finite size scaling since deviationswhich are not apparent for small ensemble sizes may well berevealed for larger ensemble sizes.We have verified Eq. (24) for L =
12 by comparing dataobtained using the parallel method with data obtained usingthe serial method (see Fig. 2).To demonstrate further the importance of randomizing Q ,we compare in Fig. 3 the distribution of ˜ γ N obtained withfixed Q given by Eq. (22) with the distribution obtainedwith random Q generated using 64 transfer matrix multipli-cations. Note that, unlike eigenvalues whose order is arbitrary,the ˜ γ i are in the order they are obtained in the QR factor-
3. Phys. Soc. Jpn.
Full Paper ization. This is not in general in decreasing order. Nor is ˜ γ N always positive, as is seen in Fig. 3 where the distributionsextend to negative values. Another important point to graspfrom Fig. 3 is that if too large an ensemble with insu ffi cientrandomization of Q is generated, the errors in the estimationof h ˜ γ N i will be systematic not random. This would make reli-able finite size scaling impossible.Though not necessary for the estimation of the criticalexponent, we can calculate all ˜ γ , · · · ˜ γ N by making Q a2 N × N orthogonal matrix. In this case, we have found thatthe ˜ γ i do not usually occur in pairs of opposite sign. The onlycondition they satisfy in general is that their sum is zero. Nev-ertheless, we have noticed that, if a su ffi ciently large numberof transfer matrix multiplications is used to generate random Q , the ˜ γ i again occur in pairs of opposite sign. This is sim-ilar to Eq. (13) but without the decreasing ordering. However,this seems to require a much larger number of transfer matrixmultiplications to generate random Q than is needed whenfocussing as above only on the distribution of ˜ γ N .
3. Numerical Simulation of Anderson’s Model of Local-isation
We have simulated ensembles of cubes with dimensions L × L × L with L = , , , ,
48 and 64 and disorder W in therange [15 . , . L =
64 the range ofdisorder was restricted to [15 . , . L and W an ensemble of 589824 sampleswas simulated and h ˜ γ N i estimated using the sample mean witha precision given by the standard error in the mean calcu-lated using the standard formulae. In percentage terms, theprecisions of the ensemble averages obtained varied between0 .
07% and 0 .
13% depending on the pair of W and L consid-ered.To avoid round-o ff error, QR factorizations were performedafter every 4 or 8 transfer matrix multiplications.To obtain a stationary distribution of initial matrices Q , 64transfer matrix multiplications were used. To check that thiswas su ffi cient we compared with data obtained with 32 and 96multiplications. The results of the Kolmogorov-Smirnov testfor L =
48 are shown in Table I. It can be seen that, while32 transfer matrix multiplications are not su ffi cient, the distri-butions of ˜ γ N obtained with 64 and 96 multiplications cannotbe distinguished with this number of samples. For L = L , inpractice virtually all the computer time is spent on the largestsystem size.
4. Finite Size Scaling Analysis
The critical disorder, critical exponent, and other criticalquantities are estimated by fitting the size and disorder de-pendence of the dimensionless quantity Γ to a one parameter scaling model Γ = h ˜ γ N i L = F ( φ ) . (27)Here, F is a scaling function and φ is a relevant scaling vari-able. The scaling function is approximated by a Taylor seriestruncated at order n F ( φ ) = n X j = F j φ j . (28)The scaling variable has the form φ = u ( W − W c ) L α , (29)where α is the inverse of the critical exponent ν = α > , (30)and W c is the critical disorder. The function u is approximatedby a Taylor series truncated at order m , u ( W − W c ) = m X j = u j ( W − W c ) j . (31)To avoid any ambiguity in the model we impose the condition F = . (32)The critical exponent ν is expected to be universal, i.e. itshould be the same for all Anderson transitions in three-dimensional systems in the orthogonal symmetry class. Thescaling function, and in particular the quantity Γ c = F = F (0) , (33)are expected to be somewhat less universal, i.e. they shouldbe the same for all Anderson transitions in three-dimensionalsystems in the orthogonal symmetry class but depend on theboundary conditions imposed in the transverse directions. To determine the best fit we perform a non-linear leastsquares fit, i.e. we minimize the χ statistic. The quality ofthe fit is assessed using the goodness of fit probability p . Boththe goodness of fit probability and the precision of the fittedparameters are determined by generating and fitting an en-semble of 500 pseudo-data sets. The details of this procedurehave already been described in Ref. 8.In Fig. 4 we show a fit to 117 data points with m = n =
3. The orders of the truncations are determined by requir-ing that the goodness of fit is greater than 0 . m and n . The de-tails of the fit and values of the fitted parameters are shown inTable II. To demonstrate one parameter scaling the collapse ofall data points onto a single curve when the data are re-plottedversus φ is shown in Fig. 5.
5. Discussion
We have described an adaptation of the transfer matrixmethod often employed in the field of Anderson localisationfor massively parallel supercomputers. We have illustrated theuse of this method by applying it to Anderson’s model of lo-calisation in three dimensions and estimated the critical expo-nent ν = . ± . . (34)
4. Phys. Soc. Jpn.
Full PaperTable I.
Example of the Kolmogorov-Smirnov test for cross section 48 ×
48. The ensemble size is 589824. The table shows the p-value returned by theKolmogorov-Smirnov test. The rows and columns are labeled by the number of randomizing multiplications. ν Γ c W c N D N P p all data 1.572[1.566,1.577] 1.7372[1.7359,1.7384] 16.543[16.541,16.545] 117 7 0.5restricted W L Table II.
The details of the finite size scaling fits to all the data, to data with the range of disorder W restricted to [16.2,16.9], and to data with larger systemsizes L = , , ,
64 only. The precisions are expressed as 95% confidence intervals. The values of the χ statistic for the best fits are 112.1, 38.5, and 59.8respectively. G W Fig. 4. (Color online) The estimates (symbols) of
Γ = h ˜ γ N i L for L = W together with the finite sizescaling fit (solid lines) described in the text. The error bars of the numericaldata are smaller than the symbol size. -1.0 -0.5 0.0 0.51.01.52.02.53.0 G f Fig. 5. (Color online) The same data as in Fig. 4 but plotted versus thevariable φ of Eq. (29) to demonstrate the collapse of the data onto a commoncurve. The line is the scaling function Eq. (27). (Note that the error here is a standard error not a 95% confi-dence interval.) The largest systems size considered here was L =
64 and the precision of the critical exponent obtainedis approximately 0.2%. This is compared to the largest sys-tem size of L =
24 in Ref. 8 where a precision of 0.35%was obtained. An important di ff erence with Ref. 8 is that wedid not need to consider corrections to scaling due to an ir-relevant scaling variable. The smallest transverse size usedhere is L =
12. As can be seen from Fig. 5 of Ref. 8, irrele-vant corrections are already less than the precision of our datafor L ≥
12. Our estimate Eq. (34) is consistent with that ob-tained by multi-fractal analysis of eigenstates and withboth numerical and experimental work on the quantum kickedrotor.
In this work we simulated an ensemble of cubes, i.e. an en-semble of systems with aspect ratio fixed to unity. However,this choice is not optimal. In our calculation, about half thetime was spent randomising the initial matrix Q and abouthalf the time estimating the ensemble average of ˜ γ N . A moreeconomical approach would be to simulate a smaller ensem-ble of a longer systems. More of the computer time wouldthen be devoted to estimating the ensemble average of ˜ γ N rather than being “wasted” on randomising Q . Indeed, pro-vided the matrices Q are randomised properly, we see fromequation Eq. (24) that there is no need to keep the aspect ratiofixed and any convenient length can be simulated. The appro-priate choice will depend on the time limits set in the queuingsystem of the supercomputer being used.In the serial method the length of the system is increaseduntil the desired precision of the Lyapunov exponent is ob-tained. Thus, the precise number of transfer matrix multipli-cations is usually not known in advance. This is inconvenientwhen we wish to simulate systems with correlated randompotentials such as quantum Hall systems or cold atomsystems. The method we describe here may be better suitedto such problems.By allowing full exploitation of current supercomputers,the method described here may also be useful when studyinghigher dimensional systems
9, 10, 35) where the time constraintsof the transfer matrix method become more severe.This work was supported by JSPS KAKENHI GrantsNo. 15H03700, 17K18763 and No. 26400393. The authorsthank the Supercomputer Center, the Institute for Solid State
5. Phys. Soc. Jpn.
Full Paper
Physics, the University of Tokyo for the use of System B.
1) J.-L. Pichard and G. Sarma: J. Phys.
C14 (1981) L127.2) A. MacKinnon and B. Kramer: Phys. Rev. Lett. (1981) 1546.3) A. MacKinnon and B. Kramer: Z. Phys.B (1983) 1.4) E. Abrahams, P. W. Anderson, D. C. Licciardello, and T. V. Ramakrish-nan: Phys. Rev. Lett. (1979) 673.5) B. Kramer and A. MacKinnon: Reports on Progress in Physics (1993) 1469.6) K. Slevin and T. Ohtsuki: Phys. Rev. Lett. (1997) 4083.7) K. Slevin and T. Ohtsuki: Phys. Rev. Lett. (1999) 382.8) K. Slevin and T. Ohtsuki: New Journal of Physics (2014) 015012.9) Y. Ueoka and K. Slevin: Journal of the Physical Society of Japan (2014) 084711.10) K. Slevin and T. Ohtsuki: Journal of the Physical Society of Japan (2016) 104712.11) Y. Asada, K. Slevin, and T. Ohtsuki: Phys. Rev. Lett. (2002) 256601.12) E. P. Wigner: Mathematical Proceedings of the Cambridge Philosophi-cal Society (1951) 790.13) F. J. Dyson: Journal of Mathematical Physics (1962) 140.14) F. J. Dyson: Journal of Mathematical Physics (1962) 1199.15) B. Huckestein and B. Kramer: Phys. Rev. Lett. (1990) 1437.16) K. Slevin and T. Ohtsuki: Phys. Rev. B (2009) 041304.17) M. Amado, A. V. Malyshev, A. Sedrakyan, and F. Dom´ınguez-Adame:Phys. Rev. Lett. (2011) 066402. 18) H. Obuse, I. A. Gruzberg, and F. Evers: Phys. Rev. Lett. (2012)206804.19) P. W. Anderson: Phys. Rev. (1958) 1492.20) K. Slevin, T. Ohtsuki, and T. Kawarabayashi: Phys. Rev. Lett. (2000)3915.21) V. I. Oseledec: Trans. Moscow Math. Soc. (1968) 197.22) K. Slevin and T. Ohtsuki: Phys. Rev. B (2001) 045108.23) K. Slevin, Y. Asada, and L. I. Deych: Phys. Rev. B (2004) 054201.24) A. Rodriguez, L. J. Vasquez, K. Slevin, and R. A. R¨omer: Phys. Rev.Lett. (2010) 046403.25) A. Rodriguez, L. J. Vasquez, K. Slevin, and R. A. R¨omer: Phys. Rev. B (2011) 134209.26) L. Ujfalusi and I. Varga: Physical Review B (2015) 184206.27) G. Lemari´e, J. Chab´e, P. Szriftgise, J. C. Garreau, B. Gr´emaud, andD. Delande: Phys. Rev. A (2009) 043626.28) G. Lemari´e, B. Gr´emaud, and D. Delande: EPL (Europhysics Letters) (2009) 37007.29) M. Lopez, J.-F. Cl´ement, P. Szriftgiser, J. C. Garreau, and D. Delande:Phys. Rev. Lett. (2012) 095701.30) T. Ando: Journal of the Physical Society of Japan (1984) 3101.31) T. Ando and H. Aoki: Journal of the Physical Society of Japan (1985)2238.32) J. Chalker and P. Coddington: J. Phys. C21 (1988) 2665.33) B. Kramer, T. Ohtsuki, and S. Kettemann: Phys. Rep. (2005) 211.34) D. Delande and G. Orso: Phys. Rev. Lett. (2014) 060601.35) A. M. Garca-Garca and E. Cuevas: Physical Review B75