Anderson transition and multifractals in the spectrum of the Dirac operator of Quantum Chromodynamics at high temperature

# Anderson transition and multifractals in the spectrum of the Dirac operator of Quantum Chromodynamics at high temperature

László Ujfalusi Elméleti Fizika Tanszék, Fizikai Intézet, Budapesti Műszaki és Gazdaságtudományi Egyetem,
H-1521 Budapest, Hungary
Matteo Giordano Institute for Nuclear Research of the Hungarian Academy of Sciences,
Bem tér 18/c H-4026 Debrecen, Hungary
Ferenc Pittler MTA-ELTE Lattice Gauge Theory Research Group,
Pázmány P. sétány 1/A H-1117 Budapest, Hungary
Tamás G. Kovács Institute for Nuclear Research of the Hungarian Academy of Sciences,
Bem tér 18/c H-4026 Debrecen, Hungary
Imre Varga Elméleti Fizika Tanszék, Fizikai Intézet, Budapesti Műszaki és Gazdaságtudományi Egyetem,
H-1521 Budapest, Hungary
July 14, 2019
###### Abstract

We investigate the Anderson transition found in the spectrum of the Dirac operator of Quantum Chromodynamics (QCD) at high temperature, studying the properties of the critical quark eigenfunctions. Applying multifractal finite-size scaling we determine the critical point and the critical exponent of the transition, finding agreement with previous results, and with available results for the unitary Anderson model. We estimate several multifractal exponents, finding also in this case agreement with a recent determination for the unitary Anderson model. Our results confirm the presence of a true Anderson localization-delocalization transition in the spectrum of the quark Dirac operator at high-temperature, and further support that it belongs to the 3D unitary Anderson model class.

###### pacs:
71.23.An, 71.30.+h, 72.15.Rn, 12.38.Gc, 12.38.Mh, 11.15.Ha

## I Introduction

The Anderson metal-insulator transition is a genuine quantum phase transition, which has been widely investigated in condensed matter physics since the seminal paper of Anderson Anderson58 (). In the past years Anderson transitions were found in a wide range of physical systems, such as ultrasound in disordered elastic networks Hu08 (); Faez09 (), light in disordered photonic lattices in the transverse direction Segev13 (), or in an ultracold atomic system in a disordered laser trap AspectBouyer12 ().

A characteristic feature of Anderson transitions is the rich multifractal structure of critical eigenstates, which has been the subject of intense research activity in recent years (see Ref. EversMirlin08, for a review). Direct signs of multifractals at the metal-insulator transition point have been observed experimentally in dilute magnetic semiconductors Richardella10 (). Multifractality can moreover influence the behavior of various systems near criticality in different ways. For example, the large overlap of multifractal wave-functions can increase the superconducting critical temperature Feigelman10Burmistrov12 (). The multifractality of the local density of states may induce a new phase because of the presence of local Kondo effects induced by local pseudogaps at the Fermi energy Kettemann12 ().

The simplest model displaying a metal-insulator transition is the Anderson model, which describes non-interacting fermions in a disordered crystal. Disorder is usually introduced through a random on-site potential, while hopping elements are fixed (up to a random phase or rotation). In this case the system belongs to one of the Wigner-Dyson (WD) symmetry classes depending on the global symmetries of the system. On the other hand, systems with vanishing on-site terms and random hopping terms, if the lattice is bipartite, possess an additional chiral symmetry and belong to one of the chiral WD classes EversMirlin08 ().

Quite surprisingly, an Anderson transition has been shown to take place also in the spectrum of the Dirac operator in Quantum Chromodynamics (QCD) at high temperature GarciaGarcia:2006gr (); Pittler (); KP (); Kovacs12 (); Giordano14 () (see Ref. Giordano:2014qna, for a review). QCD is the quantum field theory governing strong interactions at the microscopic level, and operates on length and energy scales vastly different from the ones usually encountered in condensed matter physics. QCD is a non Abelian gauge theory, describing the interactions of quarks, which are fermions, and gluons, which are the vector bosons of the gauge symmetry. Although these are the fundamental degrees of freedom, they do not appear in the spectrum of the theory at low temperatures, which contains only hadrons, i.e., bound states of quarks and gluons, due to the phenomenon of confinement. However, at a (pseudo)critical temperature, , strongly interacting matter undergoes a crossover to the so-called quark-gluon-plasma phase, and at high temperatures quarks and gluons are deconfined. This transition is accompanied by the restoration of the approximate chiral symmetry that is spontaneously broken at low temperatures Aoki:2006we ().

Contributions of quarks to observables, as well as all quark correlation functions, are entirely encoded in the eigenvalues and the eigenvectors of the Dirac operator in the background of a non Abelian gauge field. Physical quantities are then obtained after averaging over the gauge field configurations with the appropriate path-integral measure. In this respect, the eigenmodes of the Dirac operator can be formally treated as the eigenstates of a random “Hamiltonian”, with disorder provided by the gauge field fluctuations. Among the eigenmodes, a prominent role is played by the low-lying ones: for example, they give the most important contributions to the quark correlation functions, and determine the fate of chiral symmetry through the Banks-Casher relation Banks:1979yr (). The low end of the spectrum looks completely different in the hadronic and in the quark-gluon-plasma phase. At low temperatures, the density of states is finite near the origin, and both low-lying and bulk eigenmodes are extended throughout the system. In contrast, at high temperatures, above , the density of states vanishes at the origin, and the low-lying eigenmodes are localized on the scale of the inverse temperature, while higher up in the spectrum, beyond a temperature-dependent critical “energy”, , the eigenmodes are again extended KP (); Kovacs12 (). The temperature dependence of the “mobility edge” was investigated in Ref. Kovacs12, , in which it was found that extrapolates to zero at a temperature compatible with . Typical Dirac eigenmodes in the localized, critical and delocalized regimes are shown in Fig. 1. The transition in the spectrum from localized to delocalized eigenmodes has been shown to be a second-order phase transition, with critical exponent compatible with the one found in the three-dimensional unitary Anderson model Giordano14 ().

It is rather surprising at first that the Anderson transition in the high-temperature QCD Dirac spectrum seems to belong to the same universality class as that of the three-dimensional unitary Anderson model. From the point of view of statistical systems, QCD at a finite temperature is in fact a four-dimensional Euclidean model, with the “temporal” dimension compactified on a circle of length . However, it has been argued that high-temperature QCD is an effectively three-dimensional disordered system with on-site disorder, the strength of which is set by the temperature Giordano:2015vla (); Bruckmann:2011cc (). While this makes it more plausible that the two models actually belong to the same universality class, it does not make it less important to look for further evidence. In this respect, finding the same multifractal structure in the critical eigenstates would give strong support to the claim of Ref. Giordano14, , and so to the broader universality of the critical properties of Anderson transitions. The study of this multifractal structure is precisely the aim of this work.

The plan of the paper is the following. In section II we give a brief discussion of multifractality, and of the method of multifractal finite-size scaling (MFSS). In section III we describe in some detail the Dirac operator and the numerical simulations of QCD employed in this paper. In section IV we study the correlations between eigenvectors of the Dirac operator around the critical energy, both for comparison to the 3D unitary Anderson model, and for their appropriate treatment in the statistical analysis. In section V we discuss the results of MFSS for the eigenvectors of the Dirac operator. Finally, in section VI we state our conclusions.

## Ii Finite-size scaling laws for generalized multifractal exponents

In this section we briefly review wave-function multifractality and the technique of multifractal finite-size scaling. The wave-function of a particle in is naturally associated to a local probability distribution, namely , giving the probability to find the particle in an infinitesimal neighborhood of . For smooth wave functions, the probability to find the particle in a small but finite neighborhood of of size scales as . For fractal wave functions, this probability scales as , where is called the fractal dimension. For strongly fluctuating wave-functions, however, this probability scales in general as , with an -dependent power called the local dimension. In turn, points with the same local dimension, , constitute a subset of characterized by its own fractal dimension, which generally depends on . The wave-function therefore defines not one, but many different fractals, and is therefore said to be multifractal. Multifractal wave-functions are strongly fluctuating on every length-scale, and their characterization requires an infinite number of fractal dimensions, called multifractal exponents (MFEs). An example of a multifractal wave-function is shown in Fig. 1(b).

Multifractality is a known feature of critical eigenfunctions at the Anderson metal-insulator transition EversMirlin08 (), that can be studied by means of multifractal finite-size scaling (MFSS) Rodriguez_prl (). In recent high-precision calculations Rodriguez11 (); Ujfalusi15 (); Ujfalusi14 (), MFSS has been successfully employed to determine the MFEs of critical eigenfunctions, as well as to obtain a more precise estimate of the critical disorder and of the critical exponents, for Anderson models in different symmetry classes. In this work we want to perform a similar MFSS analysis to study the Anderson localization-delocalization transition in the spectrum of the Dirac operator in QCD.

In the remainder of this section we describe MFSS in some detail. Our methods and notations are essentially the same as in Ref. Ujfalusi15, , to which we refer the reader for a more detailed discussion. There is however one important difference, concerning the way in which the transition is approached. In Ref. Ujfalusi15, the transition was studied by looking at wave functions at the band center and varying the amount of disorder, . In QCD the amount of disorder is effectively set by the temperature, and it is more convenient to keep it fixed and study the transition as a function of energy, , by looking at wave-functions near the mobility edge, . Therefore, has been replaced by in the expressions of Ref. Ujfalusi15, .

Let us consider a -dimensional cubic lattice of linear size , and a critical eigenfunction of a random Hamiltonian, , defined on the lattice sites and normalized to 1. We can divide the lattice into smaller boxes of linear size , and compute the probability corresponding to the -th box as

 μk=∑→x∈boxk|ψ(→x)|2, (1)

where the sum runs over the lattice sites contained in the -th box. The generalized inverse participation ratios (GIPRs) are the moments of the box probability. The GIPRs and their derivatives read

 Rq=λ−d∑k=1μqkSq=dRqdq=λ−d∑k=1μqklnμk, (2)

where , and the sum runs over all the boxes of size . For small , the averages of and over disorder realizations follow a power-law behavior as a function of , which leads one to define the following exponents:

 Dq=limλ→01q−1ln⟨Rq⟩lnλαq=limλ→0⟨Sq⟩⟨Rq⟩lnλ. (3)

and are generalized fractal dimensions, usually referred to as multifractal exponents (MFEs). One can similarly define MFEs for localized and delocalized states by substituting critical eigenfunctions with localized or delocalized eigenfunctions in Eq. (1). In the delocalized/metallic part of the spectrum, states extend over the whole lattice, so their effective size grows proportionally to the volume, thus leading to . On the other hand, in the localized/insulating regime, states are exponentially localized, so that their effective size does not change with the system size, resulting in for , and for . At criticality, , the eigenstates are instead expected to be multifractal, with nontrivial, -dependent and .

This jump of the MFEs at the critical point happens only in an infinite system. The main idea of MFSS is to use the MFEs as order parameters for finite size-scaling. In order to do that we have to define the finite size version of the MFEs at a given energy,

 ~αensq(E,L,ℓ) = ⟨Sq⟩⟨Rq⟩lnλ, (4) ~Densq(E,L,ℓ) = 1q−1ln⟨Rq⟩lnλ, (5)

where it is understood that wave-functions of energy around are used on the right-hand side, and where the superscript ens is to remind the reader that one has to perform ensemble averaging over the different disorder realizations. and are called generalized multifractal exponents (GMFEs). Every GMFE approaches the value of the corresponding MFE at the critical point, , only in the limit . One can also define typical MFEs,

 ~αtypq(E,L,ℓ) = ⟨SqRq⟩1lnλ, (6) ~Dtypq(E,L,ℓ) = 1q−1⟨lnRq⟩lnλ, (7)

which can be used as well in a finite-size scaling analysis. However, as we said above, MFEs are defined through ensemble averaging in principle [see Eq. (3)], and when computing MFEs in Sec. V we use ensemble averaged quantities only.

In the renormalization group language, the Anderson transition is characterized by a single relevant operator AALR (), and so in the vicinity of the critical point one can derive scaling laws for the GMFEs, which can be summarized in a single equation, using a common letter, , for the GMFEs:

 ~Gq(E,L,ℓ)=Gq+1lnλGq(Lξ,ℓξ). (8)

At the critical point, the localization length diverges as , where for . The system sizes employed in this paper, however, are not big enough to justify the use of one-parameter scaling, and so we included the contribution of an irrelevant operator, , which leads us to write

 ~Gq(E,L,ℓ)=Gq + 1lnλ[Grq(ϱL1ν,ϱℓ1ν)+ (9) + ηℓ−yGirq(ϱL1ν,ϱℓ1ν)].

A second irrelevant term, proportional to , is expected to be less important and will be neglected in the analysis Rodriguez11 (); Ujfalusi15 ().

Fits to the numerical data are performed by expanding and in the variables and up to order and , respectively. The number of parameters therefore grows as . Moreover, and must also be expanded in powers of up to order and , which further increases the number of fitting parameters. The fit provides all the physically interesting quantities, namely the critical point, , the critical exponent, , the irrelevant exponent , and the MFE, .

A simpler fit can be performed by setting to a fixed value. In this case, dropping from the notation, we can write

 ~Gq(E,L)=Gq(Lξ)=Grq(ϱL1ν)+ηL−yGirq(ϱL1ν), (10)

having absorbed and the factor into new functions . The main advantage is that since are now single-variable functions, the number of expansion parameters grows only as . On the other hand, with this method one can determine only , , and , while the value of the MFE, , cannot be obtained.

## Iii Properties of the Dirac operator and details of the simulations

In this section we give the relevant details about the Dirac operator and QCD, and about how the QCD Dirac spectrum can be studied by means of numerical simulations. The continuum Euclidean Dirac operator is

 D(A)=4∑μ=1γμ(∂μ+igAμ), (11)

where are the Euclidean Dirac matrices, is the coupling constant, and is the non Abelian gauge field. More precisely, is a Hermitian matrix, where is real and the sum runs over the generators of . The Dirac operator is thus anti-Hermitian, so admitting a straightforward interpretation as ( times) the Hamiltonian of a quantum system. The Dirac operator is a chiral operator with the following structure in spinor space,

 D=i(0WW†0), (12)

with a complex matrix with no further symmetry Verbaarschot00 (). As a random matrix model, the Dirac operator in a random gauge field belongs therefore to the chiral unitary class. Chiral symmetry is expressed by the anticommutation relation , which implies that the nonzero eigenvalues come in pairs . It is thus sufficient to consider the positive part of the spectrum only.

The partition function of QCD at temperature can be expressed as a functional integral,

 ZQCD=∫[dA]e−Sg[A]∏fdet[D(A)+mf], (13)

with the constraint . The product is over the six different types of quarks (“flavors”), with the mass of quark . Here is a positive functional of the gauge field, which together with the determinants provides the probability distribution of the disorder, i.e., of the gauge field configurations. Numerical simulations of QCD require the discretization of Eq. (13) on a finite lattice. For a review of lattice QCD see, e.g., Ref. MontvayMunster, . While the discretization of the gauge fields poses no particular problem, and can be performed preserving exact gauge invariance Wilson:1974sk (), fermion fields are known to be more problematic, and the discretization of the Dirac operator spoils some of the properties of its continuum counterpart. Nevertheless, the discretization that we employed, namely staggered fermions Susskind:1976jm (), preserves the anti-Hermiticity and the symmetry of the spectrum with respect to the origin, and moreover preserves the chiral unitary symmetry class Verbaarschot00 ().

It must be noted at this point that in the case of the Anderson model, chiral and non-chiral symmetry classes differ only in their properties near the band center Evangelou03 (), i.e., , while the properties of the bulk of the spectrum are similar. For example, the authors of Ref. Evangelou03, found Wigner-Dyson statistics in the bulk spectrum of a three-dimensional chiral orthogonal disordered model. Moreover, even the critical exponent of the orthogonal and of the chiral orthogonal class turn out to be the same, up to very high numerical precision Biswas00 (). We expect the same to be true for the multifractal exponents.

Let us now describe the numerical setting in some detail. QCD is discretized on a periodic hypercubic lattice , of spatial extent in each direction and temporal extent . The gauge fields are replaced by corresponding gauge links, i.e., parallel transporters along each link of the lattice, which are elements of the gauge group, . The functional is discretized and expressed in terms of the gauge links, and the integration over gauge fields is replaced by the integration with the Haar measure over gauge links, i.e., over the gauge-group valued variables on the links. Finally, the continuum Dirac operator is replaced by the staggered Dirac operator, which reads

 Dstagxy=124∑μ=1ημ(x)[δx+^μ,yUμ(x)−δx−^μ,yU†μ(x−^μ)], (14)

with , and the gauge link connecting the lattice site to the neighboring site along direction . The staggered Dirac operator carries only spacetime and color indices, i.e., it has no spinorial structure. The eigenvalue equation must be supplemented with the antiperiodic boundary condition for the quark eigenfunction.

As we have already remarked, the Dirac operator can be viewed as a random Hamiltonian, with disorder provided by the fluctuations of the gauge fields, and distributed according to the Boltzmann weight appearing in the partition function. In its discretized version, the Dirac operator is a large sparse matrix, with nonzero random elements only in the off-diagonal, nearest-neighbor hopping terms, which depend on the parallel transporter on the corresponding link of the lattice. This resembles an Anderson model with off-diagonal disorder, although here the fluctuations of the gauge links are correlated, rather than independent. However, since the theory has a mass gap, correlations decrease exponentially with the distance. Moreover, the strong correlation between the different time-slices makes the model effectively three-dimensional, with the fluctuations of the temporal links acting effectively as a three-dimensional diagonal disorder Bruckmann:2011cc (); Giordano:2015vla (). The size of the gauge field fluctuations are determined by the temperature, which therefore is expected to play the same role as the amount of disorder in the Anderson model. This is confirmed by the fact that the temperature governs the position of the mobility edge.

In the present work we have studied the spectrum of the Dirac operator by generating gauge link configurations, i.e., realizations of disorder, by means of Monte-Carlo methods. Numerical calculations were done on a GPU cluster. In our simulations we have included only the three lightest flavors (up, down, and strange), with equal masses for the up and down quark. For many purposes, this is a good approximation of the real world. The lattice spacing in physical units was set to and the temporal size was fixed to , resulting in the temperature , well above the crossover temperature (see Refs. Pittler, ; KP, ; Kovacs12, ; Giordano14, for more details). Technical details about the numerical implementation and the scale-setting procedure can be found in Refs. Aoki:2005vt, ; Borsanyi:2010cj, . We have computed the eigenpairs of the Dirac operator from the smallest eigenvalue up to the upper end of the critical region, on lattices of spatial sizes in the range (in lattice units). A detailed list is reported in Tab. 1 along with the corresponding number of samples.

The three-dimensional box probability, Eq. (1), required for the multifractal analysis, was constructed as follows. To have a gauge-invariant description we summed over the color components, labelled by . Moreover, due to the strong correlation between the lattice time-slices, the eigenvectors of the Dirac operator look qualitatively the same on each of them, so we can also sum over the time-slices, . The squared amplitude is then defined as , and provides the basic three-dimensional spatial probability distribution, from which the box probability distribution is then obtained in the usual way.

## Iv Correlations between eigenvectors

In this section we investigate the correlations between different eigenvectors of the Dirac operator in a given gauge configuration. Our motivation is twofold. On the one hand, we want to compare the eigenvector correlations in QCD with the ones found in the unitary Anderson model. On the other hand, these correlations have to be properly taken into account when fitting the numerical data to determine the various critical quantities, as we do in Sec. V.

Cuevas and Kravtsov CuevasKravtsov07 () showed that in the Anderson model there are non-negligible correlations between eigenvectors. Similar correlations are therefore expected also in other disordered systems, like the one under consideration. The relevant quantities are the density-density correlations, which are defined in terms of the overlap integral, which for the -th and -th eigenfunctions reads

 Kij2=∫d3x|ψi|2∣∣ψj∣∣2. (15)

In the case of QCD, has the meaning discussed above at the end of Sec. III. One then defines the joint probability distribution of and of the energy difference between eigenstates,

 P(ω,k)=⟨∑i,jδ(Ei−Ej−ω)δ(Kij2−k)⟩. (16)

To characterize the average behavior of the overlap integral as a function of energy, its conditional expectation value is the natural choice,

 C(ω)=∫dkkP(ω,k)∫dkP(ω,k). (17)

The quantity is expected to be of order along the whole spectrum, where is the volume of the system. Indeed, for two delocalized states , while for two localized states is nonzero only if they happen to be in the same region, in which case it is of order 1, and the probability that this happens is of order .

Fig. 2(a) shows the eigenvector correlation in the unitary Anderson model at criticality. One can see a large enhancement of the correlation at small , and decreasing behavior with growing energy separations, which is similar to the results of Ref. CuevasKravtsov07, for the orthogonal Anderson model. Examining the same correlation for critical eigenfunctions in QCD, we also find an enhancement at small energy separations, see Fig. 2(b). In the critical regime the behavior of the two systems is very similar, and even the approximate exponent of the power-law decay is close to in both cases. This is a nice example of the similarity of the two models, and in Sec. V we present further similarities in more detail.

## V MFSS for the eigenvectors of the Dirac operator

In this section we would like to characterize the Anderson phase transition in the spectrum of the Dirac operator of QCD in the frame of the MFSS, described in Sec. II. As discussed at the end of Sec. III, a three-dimensional spatial probability distribution was calculated from the eigenvectors. From that, the GMFEs and were then computed according to Eqs. (4)–(7). More precisely, we chose values of energy, , in the range , and for the -th energy value and the -th gauge configuration we computed and according to Eq. (2). In order to decrease the numerical noise we averaged over all the eigenvectors in an energy range of width around . The GMFEs and are then obtained by averaging and over the index , i.e., over configurations, or in other words, over different realizations of disorder.

An example of the resulting GMFEs at fixed is depicted in Fig. 3. As the system size grows, the curves shift to opposite directions on the two sides of the transition. At low energy they shift down, indicating a localized phase, while at high energy they shift up, suggesting a metallic phase, as expected. In between, the curves should cross at a common point, corresponding to the critical energy, but due to finite size effects originating from the irrelevant terms this is true only approximately.

Data were then fitted with the scaling laws Eqs. (9) and (10), minimizing the quantity , using the MINUIT library James75 (). Here is the number of degrees of freedom and is the distance between the numerical data, , and the fitting function, , in the appropriate metric, i.e.,

 χ2=∑i,j(yi−fi)(C−1)ij(yj−fj), (18)

where is the covariance matrix of the data points. In the light of the results of Sec. IV, which show that there are strong correlations between eigenvectors in a given gauge configuration, strong correlations are also expected among GMFEs at different energies, and so the inclusion of correlations in the fitting procedure is necessary to obtain accurate results. The error bars of the best fit parameters were estimated by Monte-Carlo simulation, generating sets of synthetic data, distributed according to a Gaussian distribution with means equal to the raw data points and covariance matrix equal to the covariance matrix of the sample. We then determined the error bars from the distributions of the resulting fit parameters, choosing the confidence level.

In order to perform best fits, the scaling laws Eqs. (9) and (10) need to be expanded in powers of , and this requires to set the expansion orders of the relevant/irrelevant scaling term , as well as the expansion orders and of and . Since the relevant operator is more important than the irrelevant one we always used and . We then repeated the fit for several choices of the expansion orders.

The quality of the best fits was judged according to two criteria. The first criterion was how close the ratio approached unity, and only fits with were considered acceptable. The second criterion was stability against changing the expansion orders, in order to keep under control the systematic effects due to the truncation of the scaling function. We estimated the systematic error due to truncation as twice the standard deviation of the critical parameters, in the sample comprising the stable fits and the essentially equivalent ones obtained by increasing or lowering the expansion orders. The factor of 2 is required by consistency with the 95% confidence level chosen for the statistical error.

We first performed the MFSS at fixed , as described in Sec. II, both for ensemble and typical averaging. We used , as this value is compatible with several of the system sizes listed in Tab. 1. The fixed method is more stable, since the number of parameters to fit grows only linearly with the expansion orders. Stability was a serious issue, because the largest system size available, , was only about half of the one used in Refs. Rodriguez11, ; Ujfalusi15, ; Ujfalusi14, . Due to this limitation, fits were stable for adding or removing an expansion parameter only in the range . The reason for this is that, for large , the -th power in Eq. (2) strongly enhances the contribution of the few spatial points with very large (if ) or very small (if ) wave-function amplitude squared, which therefore dominate the sum. This results in an effectively reduced statistics and so in a noisy dataset, and leads to a regime , where GMFEs behave numerically the best. Notice that, by construction, and , and moreover due to a symmetry relation derived in Ref. MFME, .

The resulting critical parameters are listed in Tab. 2 and shown in Fig. 4. The results are essentially independent of and the type of averaging, as expected. We also checked that the critical parameters do not depend on the width of the energy window, , used in the computation of the GMFEs. As we show in Fig. 5, the results for and are independent of within errors. Moreover, the choice is optimal, as it leads to the best accuracy.

To quote a final result for the critical parameters, we have averaged the values of , and , and of the corresponding errors, obtained with the various GMFEs. (A weighted average, using the inverse of the error band as weight, yields similar numbers.) Our result for the critical point, , is compatible with the value reported in Ref. Giordano14, at the 2- level. On average, the systematic error on is , so negligible compared to the statistical error. Our result for the critical exponent, , agrees at the 1- level with the result of Ref. Giordano14, , and with previous results for the critical exponent of the unitary Anderson model Ujfalusi15 (); SO (). For this quantity, one has also to take into account that on average the systematic error due to truncation, , is of the same size as the statistical error. On the other hand, our value for the irrelevant exponent, , is significantly different from the value of Ref. Ujfalusi15, , . It is well known that it is very difficult to determine irrelevant exponents accurately, and to explain this discrepancy further work and higher-quality data are needed. It is possible that for the system sizes presently available, more than one irrelevant term gives important contributions, so that our result for would be a sort of “effective” irrelevant exponent. In any case this point requires further analysis.

As a final remark, we note that since results for different are strongly correlated, there is no significance in the fact that our values for the critical point are systematically lower, and the ones for the critical exponent are systematically higher than the reference values.

The convergence of the fixed- MFSS confirms the presence of a critical point in the QCD Dirac spectrum where the system undergoes a true localization-delocalization transition, employing completely different observables than the ones used in Ref. Giordano14, . The results of our analysis also provide further evidence that the transition in the QCD Dirac spectrum belongs to the universality class of the 3D unitary Anderson model. Moreover, despite the fact that it does not provide the values of the MFEs, the convergence of this method also strongly indicates the presence of multifractality at the critical point.

We next procedeed to apply the variable- method, in order to try and determine the multifractal exponents, and compare them to the ones obtained for the unitary Anderson model. However, this method requires small values of to work properly, and is also more demanding as it is a two-variable fit. In practice, the ratio reached a value close to unity only if we left out the smallest system sizes, below , and if we used data corresponding to and only. Although using and improved the convergence, the fits were still unstable against changing the expansion orders. This can be understood, as a similar amount of independent data is available as in the fixed- method, but there are many more parameters to fit, as discussed in Sec. II. In order to be able to estimate the MFEs, we then fixed the critical energy and the critical exponent to the values obtained with the fixed- method, and , in this way stabilizing the fits. The systematic uncertainty corresponding to this procedure was estimated by repeating the fits with and fixed to one of the four possible combinations of the values and , which are the lower and upper boundaries of the confidence interval of and , respectively. The largest and smallest values obtained in this way were then used as upper and lower error bar on the MFEs. We experienced that the main source of uncertainty comes from the choice of , while fits are much less sensitive to the choice of . Moreover, statistical errors (estimated by Monte-Carlo) and systematic errors due to truncation were comparatively negligible.

The results of this procedure are depicted in Fig. 6. A set of nontrivial MFEs was obtained, thus providing direct evidence of the multifractality of the critical eigenfunctions of the QCD Dirac operator. Moreover, our results for the MFEs in QCD are compatible with the ones obtained in the unitary Anderson model, which further confirms that the transition belongs to the chiral unitary Anderson class.

## Vi Summary

We investigated the Anderson transition in the spectrum of the Dirac operator of QCD at high temperature, found by the authors of Ref. Giordano14, . While that work made use of spectral statistics, our aim in this paper was to examine the transition by studying the eigenvectors, and their multifractal properties at the critical energy. The results of Ref. Giordano14, for the correlation length critical exponent suggested that the Anderson transition in QCD belongs to the same universality class as the three-dimensional unitary Anderson model. We therefore looked for more similarities between these models.

First we examined the correlations between eigenvectors of a given gauge configuration. We found strong correlations between eigenmodes of the QCD Dirac operator, decreasing with energy separation in a similar way as in the unitary Anderson model. We then performed two multifractal finite-size scaling (MFSS) analyses, one with fixed ratio of the coarse-graining box size to the system size, and one with variable . MFSS with the fixed- method allowed an alternative determination of the critical point and of the critical exponent, which is in agreement with the findings of Ref. Giordano14, , and, for the critical exponent, with those of Refs. Ujfalusi15, ; SO, for the unitary Anderson model. To perform MFSS with the variable- method and determine the multifractal exponents (MFEs), we performed fits fixing the critical energy and the critical exponent to the values obtained with the fixed- method. The resulting MFEs are compatible with the MFEs found in the unitary Anderson model.

In conclusion, our work confirms the presence of an Anderson metal-insulator phase transition in the spectrum of the Dirac operator in high-temperature QCD, and provides further evidence that this transition belongs to the three-dimensional unitary Anderson model class. Morever, we have shown that the critical wave-functions of the Dirac operator are multifractals. The physical consequences of the QCD Anderson transition and of multifractality still largely need to be explored, and may lead in particular to a better understanding of the QCD chiral transition. Further work along these lines might prove beneficial for condensed matter physics as well, as it approaches the subject of localization/delocalization transitions from a broader perspective.

###### Acknowledgements.
We thank the Budapest-Wuppertal group for allowing us to use their code to generate the gauge configurations. Financial support to LU and IV from OTKA under Grant No. K108676, and from the Alexander von Humboldt Foundation is gratefully acknowledged. MG and TGK are supported by the Hungarian Academy of Sciences under “Lendület” grant No. LP2011-011. FP is supported by OTKA under the grant OTKA-NF-104034.

## References

• (1) P. W. Anderson, Phys. Rev. 109, 5 (1958).
• (2) H. Hu, A. Strybulevych, J. H. Page, S. E. Skipetrov, and B. A. van Tiggelen, Nature Physics 4, 945 (2008).
• (3) S. Faez, A. Strybulevych, J. H. Page, A. Lagendijk, and B. A. van Tiggelen Phys. Rev. Lett. 103, 155703 (2009).
• (4) M. Segev, Y. Silberberg, and D. N. Christodoulides, Nature Photonics 7, 197 (2013).
• (5) F. Jendrzejewski, A. Bernard, K. Müller, P. Cheinet, V. Josse, M. Piraud, L. Pezzé, L. Sanchez-Palencia, A. Aspect, and P. Bouyer, Nature Physics 8, 398 (2012).
• (6) F. Evers and A. D. Mirlin, Rev. Mod. Phys. 80, 1355 (2008).
• (7) A. Richardella, P. Roushan, S. Mack, B. Zhou, D. A. Huse, D. D. Awschalom, and A. Yazdani, Science 327, 665 (2010).
• (8) M. V. Feigel’man, L. B. Ioffe, V. E. Kravtsov, and E. Cuevas, Annals of Physics 325, 1390 (2010); I. S. Burmistrov, I. V. Gornyi, and A. D. Mirlin, Phys. Rev. Lett. 108, 017002 (2012).
• (9) S. Kettemann, E. R. Mucciolo, I. Varga, and K. Slevin, Phys. Rev. B 85, 115112 (2012).
• (10) A. M. García-García and J. C. Osborn, Phys. Rev. D 75, 034503 (2007).
• (11) F. Pittler, Poisson - Random Matrix transition in the QCD Dirac spectrum, Ph.D. thesis (2013).
• (12) T. G. Kovács and F. Pittler, Phys. Rev. Lett.  105, 192001 (2010).
• (13) T. G. Kovács and F. Pittler, Phys. Rev. D 86, 114515 (2012).
• (14) M. Giordano, T. G. Kovács, and F. Pittler, Phys. Rev. Lett. 112, 102002 (2014).
• (15) M. Giordano, T. G. Kovács, and F. Pittler, Int. J. Mod. Phys. A 29, 1445005 (2014).
• (16) Y. Aoki, G. Endrődi, Z. Fodor, S. D. Katz, and K. K. Szabó, Nature 443, 675 (2006).
• (17) I. Montvay and G. Münster, Quantum fields on a lattice (University Press, Cambridge UK, 1994).
• (18) T. Banks and A. Casher, Nucl. Phys. B 169, 103 (1980).
• (19) M. Giordano, T. G. Kovács, and F. Pittler, JHEP 04, 112 (2015).
• (20) F. Bruckmann, T. G. Kovacs and S. Schierenberg, Phys. Rev. D 84, 034505 (2011).
• (21) A. Rodriguez, L. J. Vasquez, K. Slevin, and R. A. Römer, Phys. Rev. Lett. 105, 046403 (2010)
• (22) A. Rodriguez, L. J. Vasquez, K. Slevin, and R. A. Römer, Phys. Rev. B 84, 134209 (2011).
• (23) L. Ujfalusi and I. Varga, Phys. Rev. B 91, 184206 (2015).
• (24) L. Ujfalusi and I. Varga, Phys. Rev. B 90, 174203 (2014).
• (25) E. Abrahams, P. W. Anderson, D. C. Licciardello, and T. V. Ramakrishnan, Phys. Rev. Lett. 42, 673 (1979).
• (26) F. James and M. Roos, Comput. Phys. Commun. 10, 343 (1975).
• (27) J. J. M. Verbaarschot and T. Wettig, Annu. Rev. Nucl. Part. Sci. 50, 343 (2000).
• (28) K. G. Wilson, Phys. Rev. D 10, 2445 (1974)
• (29) L. Susskind, Phys. Rev. D 16, 3031 (1977); J. B. Kogut and L. Susskind, Phys. Rev. D 11, 395 (1975); T. Banks et al. [Cornell-Oxford-Tel Aviv-Yeshiva Collaboration], Phys. Rev. D 15, 1111 (1977).
• (30) S. N. Evangelou and D. E. Katsanos, J. Phys. A: Math. Gen. 36, 3237â3254 (2003).
• (31) P. Biswas, P. Cain, R. A. Römer, and M. Schreiber, phys. stat. sol. b 218, 205 (2000).
• (32) Y. Aoki, Z. Fodor, S. D. Katz, and K. K. Szabó, JHEP 01, 089 (2006).
• (33) S. Borsányi, G. Endrődi, Z. Fodor, A. Jakovác, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabó, JHEP 11, 077 (2010).
• (34) E. Cuevas and V. E. Kravtsov, Phys. Rev. B 76, 235119 (2007).
• (35) A. D. Mirlin, Y. V. Fyodorov, A. Mildenberger, and F. Evers, Phys. Rev. Lett. 97, 046803 (2006).
• (36) K. Slevin and T. Ohtsuki, Phys. Rev. Lett. 78, 4083 (1997).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters