Extreme statistics of the excitations in the random transverse Ising chain

In random quantum magnets, like the random transverse Ising chain, the low energy excitations are localized in rare regions and there are only weak correlations between them. It is a fascinating question whether these correlations are completely irrelevant in the sense of the renormalization group. To answer this question, we calculate the distribution of the excitation energy of the random transverse Ising chain in the disordered Griffiths phase with high numerical precision by the strong disorder renormalization group method and - for shorter chains - by free-fermion techniques. Asymptotically, the two methods give identical results, which are well fitted by the Fr\'echet limit law of the extremes of independent and identically distributed random numbers. Given the finite size corrections, the two numerical methods give very similar results, but they differ from the correction term for uncorrelated random variables. This fact shows that the weak correlations between low-energy excitations in random quantum magnets are not entirely irrelevant.


I. INTRODUCTION
Many-body systems in the presence of quenched disorder have unusual dynamical properties due to rare-region effects. In these systems -due to extreme fluctuations of strong couplings -domains are formed, which can remain locally ordered even in the paramagnetic phase. The relaxation time, τ , associated with turning the spins in such domains can be extremely large and it has no upper limit in the thermodynamic limit. This type of Griffiths singularities are responsible for non-analytical behaviour of several physical quantities (susceptibility, specific heat, auto-correlation function) in the so called Griffiths phase, which is an extended part of the paramagnetic phase 1 .
In random quantum systems Griffiths-effects are stronger than in classical ones, which is manifested in power-law decay of the auto-correlation function, as well as power-law singularities of the susceptibility and the specific heat at low temperatures 2,3 . In random quantum magnets with discrete symmetry, such as in the random transverse Ising model (RTIM) the low-energy excitations are localised and their properties can be successfully studied by the so called strong disorder renormalization group (SDRG) method 4 . As initiated be Ma, Dasgupta and Hu 5 and further developed by Fisher 6 the SDRG is a local renormalization method, in which quantum and disorder fluctuations are treated at the same time and degrees of freedom with a large excitation energy are successively eliminated. In these random quantum magnets the SDRG method is expected to provide asymptotically exact results not only at the critical point, the properties of which are governed by an infinite disorder fixed point 7 , but in the Griffiths phase as well, at least for the dynamical singularities.
In phenomenological descriptions it is often assumed that the localized excitations in random quantum magnets are independent 8 . For example distribution of lowenergy excitations in these systems are well approximated by the Fréchet distribution 9 , which represents the limit law of the extremes of independent and identically distributed (iid) random numbers. Recently extreme value statistics (EVS) has been applied to several problems, we can mention earthquakes, tsunamis, extreme flooding, big wildfire, extremes of climate, stock market risks in finance, sport records, etc [10][11][12][13][14][15] . In a mathematical point of view complete understanding of extreme value statistics is known for uncorrelated random numbers, in which case also the convergence to the limit laws is derived by mathematical 16 and renormalization group [17][18][19][20] (RG) methods. The iid limit distributions also apply to weakly correlated random numbers, while new types of limit distributions appear for the strongly correlated cases 21 .
It is an intriguing problem to what extent the localised excitations in random quantum magnets are independent? Whether the weak correlations between the rare regions are completely irrelevant or these are manifested in some effects, such as in the form of finite-size corrections? In this paper we make an effort to answer this question and consider the paradigmatic model, the RTIM in one dimension (1D) and calculate numerically the distribution of the first energy gap by the asymptotically exact SDRG method with high accuracy. For moderate L system sizes, we use also free-fermion techniques [22][23][24] to calculate the gaps of the random samples exactly. The distribution of the gaps for finite L are compared with an appropriate Fréchet distribution and their difference is analysed through finite-size scaling. To check the potential form of numerical errors the same type of numerical test is repeated for iid random numbers, too.
The rest of the paper is organised in the following way. In Sec.II the RTIM model is introduced and in Sec.III the methods to calculate its excitation energy are described. In Sec.IV the distributions of gaps of the RTIM are calculated and the finite-size corrections are compared with the analytical results of uncorrelated variables. To test the numerical accuracy we repeat this analysis for uncorrelated Kesten variables 25 . Our paper concludes with a discussion in the final section. Related results about the distribution of extremes of uncorrelated variables are recapitulated in the Appendix.

II. MODEL AND KNOWN RESULTS
Here, we consider the RTIM in 1D defined by the HamiltonianĤ in terms of the σ x,z i Pauli matrices at site i and the nearest neighbour couplings, J i > 0 and the transverse fields, h i > 0, are taken from the distributions, π 1 (J) and π 2 (h), respectively. Generally we use open boundary conditions (OBC-s) and work at zero temperature, T = 0.
In the thermodynamic limit, L → ∞, the controlparameter is defined as 6 : where [· · · ] av denotes averaging over quenched disorder and var(x) stands for the variance of x. For δ < 0 (δ > 0) the system is in the ordered ferromagnetic (disordered paramagnetic) phase and at δ = 0 there is a random quantum critical point. According to SDRG calculations 6,26,27 and numerical results 28,29 the critical behaviour of the system is controlled by an infinite disorder fixed point 7 . For example, the energy scale, defined by the lowest gap, ε, and the length-scale are related as In the paramagnetic phase this relation is in a power-law form which is due to Griffiths singularities. Here, the dynamical exponent, z = z(δ), depends on the distance from the critical point and is given by the positive root of the equation 26,27,30 : which in the vicinity of the critical point diverges as: The distribution of the first gap has been calculated analytically through the SDRG method 31 : in terms of u = u 0 ε(L/ξ) z , ξ being the correlation length and u 0 is a constant, defined by the standardisation. This relation is valid for L ≫ ξ and for δ ≪ 1, in which limit ξ ∼ δ −2 . We observe that Eq.(6) is just the Fréchet distribution, see in Eq. (26). An approximate form of the distribution of the lowenergy excitations can be obtained from the assumption that these excitations are localised and are due to extreme fluctuations of say n ≪ L consecutive strong bonds 8 . The probability to find such a strongly connected cluster in the system is given by: P L (n) ∼ L exp(−an). At the same time the excitation energy due to such a cluster is exponentially small: ε ∼ exp(−bn). Combining these expressions we obtain P L (ln ε) ∼ Lε 1/z , with z = b/a being the dynamical exponent in agreement with Eq.(4). Then, the cumulated distribution of the relaxation times, τ ∼ 1/ε is obtained in this approach: If we assume that the excitations are uncorrelated following the reasoning in the Appendix we arrive at the distribution in Eq.(6), which is calculated by the SDRG approach in the given limits.
In this paper, we aim to study this problem in more detail consider the following points. i) To calculate the distribution function through numerical iteration of the SDRG approach and to check if the result in Eq.(6) is valid in the entire Griffiths phase. ii) To confirm the results of the SDRG calculation with the exact gaps obtained through free-fermion techniques. iii) To check the form of the finite-size corrections of the two methods and to compare with the analytical result for the extremes of uncorrelated random variables.

III. METHODS TO CALCULATE THE GAP IN THE RTIM
The energy scale of the model is given by the lowest excitation energy of the Hamiltonian in Eq.(1). This is calculated for finite chains by the asymptotically exact SDRG method through iteration and for shorter chains by free-fermion techniques.

A. The SDRG method
In the SDRG procedure 4 we perform consecutive decimation steps, each time considering the local excitations, say at position i. These excitations correspond to either couplings or sites, having the value of the associated gaps: 2J i and 2h i , respectively. These gaps are sorted in descending order and the largest one, denoted by Ω, which sets the energy-scale in the problem, is eliminated. Then, between remaining degrees of freedom, new terms in the Hamiltonian are generated through perturbation calculation. This procedure is successively iterated, during which Ω monotonously decreases. At the fixed point, with Ω * = 0, one makes an analysis of the distribution of the different parameters and calculates the scaling properties. In the following, we describe the elementary decimation steps.

Strong-coupling decimation
In this case, the largest local term in the Hamiltonian is a coupling, say Ω = J i , connecting sites i and i + 1 and the two-site Hamiltonian is given bŷ The spectrum ofĤ cp contains four levels, the lower two being separated from higher two by a gap of 2J i . We omit the higher two levels, corresponding to merging the two strongly coupled sites into a spin cluster in the presence of a (renormalized) transverse fieldh, the value of which is given by second-order perturbation calculatioñ

Strong-transverse-field decimation
In this case, the largest local term is a transverse field, say h i , and due to its large value this site does not contribute to the longitudinal magnetisation and therefore it is eliminated. The renormalised coupling between the remaining sites is given bỹ which is calculated by second order perturbation method.
To calculate the smallest gap of a given sample, ε, we perform (L − 1) decimation steps up to the last spin cluster having an effective transverse fieldh = ε/2.

B. Free-fermion technique
In this method,Ĥ is expressed in terms of spinless free fermions [22][23][24] . In the first step the spin operators σ x,y,z i are mapped to fermion creation (annihilation) operators c † i (c i ) by using the Jordan-Wigner transformation 22 :

and the Ising Hamiltonian in Eq.(1) is written in a quadratic form
In the second step, the Hamiltonian in Eq.(11) is diagonalized through a canonical transformation 23 , in terms of the new fermion creation (annihilation) operators η † The energies of free fermionic modes, ǫ k , are given by the eigenvalues of a 2L × 2L tridiagonal matrix and we consider only the ǫ k ≥ 0 part of the spectrum 32 .

IV. NUMERICAL RESULTS
In the numerical calculation, the parameters of the Hamiltonian, J j and h j are taken from box-like distributions and set the energy-scale with J 0 = 1. In this case the critical point of the RTIM is located at h 0 = 1 and in the disordered phase, h 0 > 1, the dynamical exponent satisfies the equation: see in Eq. (5). In the vicinity of the critical point z diverges as z ≈ 1 ln h0 , h 0 → 1 + , whereas for large h 0 it approaches 1 as z ≈ 1 + 1 2h0 , h 0 → ∞. In the actual calculations we considered three points of the paramagnetic phase, h 0 = 2 (z = 1.747655), h 0 = 3 (z = 1.33542) and h 0 = 4 (z = 1.2112289). In the freefermion calculation of the gaps the lengths of the systems were L = 16, 24, 32, 48 and 64, whereas with the SDRG method larger systems up to L = 512 are treated. In all cases 10 10 independent samples are investigated.

A. SDRG gaps and the Fréchet distribution
We start by analysing the data collected by the SDRG method and presenting the distributions of the log-gaps in Fig.1 for different sizes. To test their relation with the Fréchet distribution, we calculated in each case the best fit of the analytical function in Eq.(32), having the dynamical exponent, z L ≡ γ L and the position of the maximum, x 0 = x 0 (L) as fit parameters. These are shown in Fig.1, too. As seen in Fig.1, the Fréchet distribution describes the numerical data well and the difference between the numerical points and the fitting curve is decreasing for increasing values of L. At a given value of h 0 the fitted value of the dynamical exponent, z L , has a small variation with the size, and approaches the exact asymptotic value, as given in Eq. (15). At the value h 0 = 4 this is illustrated in Table I, where the second column shows the fitted, finite size values of z L , while the third column shows the difference from the exact value. As a first step, we used a power-law form z − z L ∼ L −a to fit the finite-size corrections, and estimates for the exponent a are calculated through two-point fit. These are listed in the fourth column of the table and have a slow convergence, which generally indicates a strong correction to scaling term. Having this possibility in mind we have used another functional form: In this case the effective exponents for ω are also calculated through two-point fit and are presented in the fifth column of the table. This type of fitting turned out to be more stable, the effective exponents seem to converge to ω ≈ 1.33. We have repeated the same analysis at the other two points of the Griffiths phase. In both cases the logcorrection form is found to provide the better fit, with the correction exponents ω ≈ 1.83 at h 0 = 2 and ω ≈ 1.50 at h 0 = 3.
Most importantly, at a given value of h 0 in Fig.1 the distributions are shifted with increasing size, and according to Eq.(4) the position of the maximum is expected to follow the rule: x 0 (L) ≈ const. + z L ln L. Comparing the position of the maximum of the distribution at two sizes one can obtain estimates for the dynamical exponent as: z(L, 2L) = (x 0 (2L) − x 0 (L))/ ln 2. We have checked that generally z L < z(L, 2L) < z 2L and these estimates approach the asymptotic exact value with the same type of corrections as noticed for the case of z L in the previous paragraph.
We note that analysis of the data obtained by freefermion calculation of the gap gives similar results, but due to smaller values of L the asymptotic region of the effective exponents, z L is more remote.

B. Finite-size corrections to the gap distributions
We considered the gaps calculated by SDRG iteration for L = 24, 32, 48, 64, 96 and 128 and for comparison we calculated those by free-fermion techniques as well, for shorter chains with L = 16, 24, 32, 48 and 64, except for h 0 = 4, where the corrections are the smallest and we went up to L = 48. We analysed finite-size scaling of distributions in two different ways. First, we used the first standardisation convention in Eq.(23) and calculated the difference with the Fréchet extreme distribution, with an effective z L , calculated from the relations below Eq.(32). This difference is then rescaled by a factor ∆z = z − z L , and the results for h 0 = 2, 3 and 4 are drawn in Fig.2. In these figures the analytical results calculated for iid random numbers in Eq.(30) with γ = z and γ ′ = −1 are also presented. Here, the correction to scaling exponent, γ ′ = −1 corresponds to the expected scaling form in Eq. (16), which is obtained through numerical analysis of the data in Table I. According to Fig.2, we can draw the following conclusions. i) The scaled finite-size difference of the distribution function seems to approach a limiting curve for large L, which depends on z(h 0 ). ii) At a given z(h 0 ) the limiting curves are similar (if not identical) for the SDRG and the free-fermion data. iii) The convergence to this limit curve is slow, much slower than that of the iid random numbers, see later in Sec.IV C and the figures in the third column of Fig.2. This slow convergence is probably related to the logarithmic correction to the scaling of the dynamical exponent. iv) Finally, the expected limit curve of the numerical distribution differs from that of the analytical result calculated for iid random numbers (having identical parameter, z). Even though the overall shapes of the curves are similar, there are noticeable differences. In particular, the low-energy part of the numerical curves are stronger represented in the numerical curves, which can be interpreted as the reduction of the value of the gap due to small, but relevant correlations between the rare regions.
We have repeated the analysis with the distribution of the log-gaps and using the second standardisation condition in Eq. (24). The results about the finite-size corrections of the distributions are shown in Fig.3 together with the analytical curves for iid random numbers. In Fig.3 the distributions of the log-gaps have a faster finite-size convergence, than those of the gaps in Fig.2. The curves in Fig.3 are almost indistinguishable for larger sizes. In addition, we notice a difference between the shape of the numerical curves and the analytical results, the former being somewhat shifted to the right around zero. This indicates a reduction of the value of the gap due to small, but relevant correlations between the rare regions.
C. Numerical test for uncorrelated Kesten variables In order to test the possible convergence of the finitesize corrections for uncorrelated variables we have repeated the analysis in the previous subsection for a parent distribution generated by Kesten random numbers. Here, we remark that Kesten-type random variables are defined as 25 : where the s j -s are iid random numbers. It is known, that in the limit of m → ∞ there is a limit distribution, ρ ∞ (u), provided [ln s] av < 0. This limit distribution has a power-law tail where the exponent is the positive root α > 0 of the equation Note that the relation for α is analogous to the equation for the dynamical exponent, z of the RTIM, see in Eq.  Fig.2. Here, the analytical results are obtained with the exponents: γ = z(h 0 ) and γ ′ = −1, the latter follows from the analytical results in Ref. [33]. It is seen in Fig.2 that there is an overall good agreement between the numerical and analytical results. There is some size-dependence of the corrections, which changes sign between h 0 = 2 and h 0 = 3. In comparison to the distribution of the gaps in the RTIM, the agreement with the analytical results is much better and the finite-size corrections are smaller.
We have repeated the analyses of the data by considering log-variables, as given in Eqs. (31) and (32) and using the second standardisation condition in Eq. (24). Results of the numerical analysis in this case are given in the third column of Fig.3. As seen in this figure, there are some finite-size dependences of the scaled numerical data, but the expected asymptotic curves agree well with the analytical results. We note that practically no finite-size dependence of the correction term is observed for a pure power parent distribution in Fig.4 in the Appendix. , and by the free-fermion method (middle panels) with the first standardisation in Eq.(23) and compared with the analytical results with γ = z(h0) and γ ′ = −1 (dashed lines). In the right panels the maximum of uncorrelated Kesten random variables are shown, having the same scaling exponents, see text. The scaling exponents z(h0) correspond to h0 = 4 (first row), h0 = 3 (second row) and h0 = 2 (third row), see in Eq. (15).

V. DISCUSSION
In this paper we have considered a paradigmatic model of random quantum magnets, the random transverse Ising model in 1D and studied the distribution of lowenergy excitations in the paramagnetic Griffiths phase, with extensive numerical methods. We have considered a large set of random samples (10 10 ) and the calculation is performed by the approximate, but asymptotically correct SDRG method (up to L = 512) and for comparison we also used the free-fermion method for shorter chains (up to a size L = 64). Analysing the distribution of the gaps we have demonstrated with high precision that -in agreement with previous expectations -the limit distribution in the thermodynamic limit is in the Fréchet form. The Fréchet distribution depends on the value of the dynamical exponent, z and we have shown that a powerful method of calculation is through fitting a Fréchet curve to the numerical gap distributions. In this way, effective, finite-size estimates are obtained for the dynamical exponent z L , which then are extrapolated to L → ∞. According to the numerical data this convergence is in the form: z − z L ∼ ln ω L/L, where the exponent of the logarithm, ω, depends on the distance from the critical point.
More interestingly, we have systematically studied the finite-size corrections to the limit law and showed that the difference between the numerical distribution and the asymptotic Fréchet form scales with z − z L , and this scaled difference, P 1 (ε; z) is a unique function of the value of the gap, ε. We have performed this type of analysis for the gap, using the first standardisation condition in Eq. (23), as well as for the log-gap, when the second standardisation condition in Eq. (24) was used. In both cases the asymptotic form of the scaled finite-size correction function are found similar (if not identical) for the data with the free-fermion calculation and that from the   SDRG iteration. We can thus conclude that the SDRG method provides not only the correct asymptotic form of the gap distribution function, but the correct finite-size correction function as well.
The measured finite-size correction functions are also compared with the analytical results of iid random numbers, having the same decay exponent (γ = z for the gaps and γ = 0 for the log-gaps) and correction to scaling exponent, γ ′ = −1. The two curves are found to have similar shape, but there are also differences in the asymptotic forms. We have checked that the observed differences are larger than the statistical error of the calculation. For this purpose, we have analysed with identical methods the same set of uncorrelated Kesten random variables having the same characteristic exponent, z. For these iid random numbers the asymptotic finite-size corrections are found to be well described by the analytical results.
The observed difference in the finite-size corrections between the numerical curves and the analytical iid results indicates that the weak correlations between lowenergy excitations in the RTIM are not completely irrel-evant. This is probably related to the fact that the linear extension of rare regions, n, with extreme fluctuations of the strong couplings, scales as n ∼ ln L, see the reasoning above Eq.(7). Since n is not limited, the finite-size dynamical exponents contain a logarithmic multiplicative factor as given in Eq. (16).
Our investigations have concluded on the RTIM in 1D, but the observed results are possibly valid for other random quantum systems as well with localised excitations. For example, we can mention the RTIM in higher dimensions [34][35][36] , the random quantum Potts 37 and Ashkin-Teller models 38 and generally random quantum magnets with short range interactions and with discrete symmetry. Similar conclusions apply to some nonequilibrium processes, such as to the random asymmetric exclusion process 39 and the random contact process 40 as well.
and similarly for the extreme density P N (v) = dM N /dv Here, the a N and b N are free up to an additive constant, the value of which is fixed by different standardisation conditions. For the analytical calculation the condition: is convenient to use, whereas for analysing numerical data it is often better to require provided the second moment of the distribution exists.
In the paper, we refer to Eq.(23) and Eq.(24) as first and second standardisation condition, respectively. According to EVS theory, the limit distributions can be of three forms, depending on the large u tail of the parent distribution. In a renormalization group treatment 17,18 with the first standardisation in Eq.(23), the limit distributions are the fixed-point solutions and are given in the form P (v; γ) = (1 + γv) −1/γ−1 exp −(1 + γv) −1/γ , for 1 + γv ≥ 0 and the different universality classes depend on the value of γ. For γ > 0 it corresponds to a parent distribution µ(u) ≈ 1 − Au −1/γ , u ≫ 1 , which represents the Fréchet universality class. For γ = 0 the asymptotic approach of the cumulative parent distribution is faster than a power and represents the Gumbel universality class. Finally, for γ < 0 the asymptotic value is approached at a finite upper border with a power −1/γ and is given as the Weibull distribution. Finite-size corrections in the RG treatment are given in the form where and ψ(v; γ, γ ′ ) = (1 + γv) + γ ′ v − (1 + γv) γ ′ /γ+1 γ ′ (γ ′ + γ) .
Consequently, the correction term to the cumulative distribution M 1 (v; γ, γ ′ ) and that of the density, P 1 (v; γ, γ ′ ) = dM 1 (v; γ, γ ′ )/dv depend on the decay exponent of the parent distribution γ, as well as the correction exponent γ ′ . In practical analysis of the data with a parent distribution in Eq. (27), it is convenient to use logarithmic variables, x = ln u, so that e x−x0 = 1 + γv. Since the parent distribution of x is exponential, the extreme distribution is given by the Gumbel form and similarly for the probability densitỹ so thatM (x 0 , γ) = 1/e andP (x 0 , γ) = 1/(eγ). In this case, the finite-size correction term is given bỹ P 1 (x; 0, γ ′ ). In the numerical examples, we have γ = 0 and γ ′ = −1, and the correction to the density in the second standardisation condition is given by 41 withx = ax + b, a = π/ √ 6 and b = γ E = 0.5772156649 being the Euler-Mascheroni constant. As an illustration, we consider a pure power parent distribution µ(u) = 1 − u −1/γ , u ≥ 1 and in Fig.4 show the numerically calculated finite-size corrections to the distributions of the log-extremes, denoted by ln(τ ), for L = 16, 24, 32, 48 and 64 with γ = z(h 0 = 3) in the second standardisation condition. These are to be compared with the analytical result in Eq. (33), since in this case γ ′ = −1. The agreement is almost perfect.