Main

Optimization problems are ubiquitous in everyday life (think, for instance, of deciding the best delivery route or scheduling jobs for the different tools in a factory). These problems can be mathematically formalized: N entities (for example, jobs queuing for their appropriate tools) compete as they try to satisfy their mutually contradictory goals. The overall frustration produced by a particular solution is quantified through a cost function, which one attempts to minimize. This task is best solved with the help of a computer, even for quite small N. Research into computational complexity studies how the computational resources (memory, computing time and so on) grow with N (ref. 14). If, for all known algorithms, the necessary resources grow faster with N than any polynomial, for example, like N!, the problem is considered hard. A small subset of these problems, named NP-complete, is of particular interest: if an efficient algorithm (with resources scaling polynomially in N) were discovered for any of the NP-complete problems, then a vast family of hard optimization problems in this subset would become easy. For physicists, the most familiar example of an NP-complete problem is finding the minimal energy state—the ground state—of an Ising spin-glass Hamiltonian on a non-planar graph15,16. This explains the surge of hardware specifically designed for minimizing a spin-glass Hamiltonian through a variety of algorithms and physical principles (see, for example, refs. 17,18,19,20,21,22).

Specifically, the strategy that concerns us here is quantum annealing. Both in the original formulation11, and also in its hardware implementation1,20, the aim is to solve the situation for spin glasses. In particular, D-wave chips solve Ising spin glass in D = 2 spatial dimensions (\({D}^{{\prime} } > 2\) can be coded over a D-wave’s D = 2 graph23; see pages 13 and 49 in ref. 24 for more on the definition of D).

Spin glass is the paradigmatic statistical model for quenched disorder25. In a transverse field, the Hamiltonian for S = 1/2 spins is

$$H=-\frac{1}{2}\sum _{x,y}\,\left[\,{J}_{x,y}{\sigma }_{x}^{Z}{\sigma }_{y}^{Z}\right]\,-\,\varGamma \sum _{x}\,{\sigma }_{x}^{X},$$
(1)

where Jx,y are the random couplings that define the problem instance under consideration, Γ is the transverse field and \({\sigma }_{x}^{X}\) and \({\sigma }_{x}^{Z}\) are, respectively, the first and third Pauli matrices acting on the spin at site x. The phase diagram for a two-dimensional interaction matrix Jx,y is sketched in Fig. 1a. For Γ = ∞, in the ground state, all spins are as much aligned with the transverse field as quantum mechanics allows them to be. (Paradoxically, from the point of view of the computational basis that diagonalizes the \({\sigma }_{{\boldsymbol{x}}}^{Z}\) matrices, this ground state seems to be a totally random statistical mixture). As Γ is diminished at zero temperature, the ground state varies. In particular, at Γ = 0, the ground state encodes the solution of the optimization problem we are interested in. At some point during the annealing, Γ goes through the critical value Γc that separates the disordered ground state from the spin-glass ground state, which has a glassy order in the computational basis. This is not just theoretical daydreaming. In a recent experiment conducted on a D-wave chip23, some 5,000 qubits displayed coherent quantum dynamics as Γ went through Γc, for annealings lasting several nanoseconds.

Fig. 1: Phase diagram and critical scaling for the two-dimensional quantum spin glass.
figure 1

a, Phase diagram for a two-dimensional Ising spin glass in terms of temperature T and transverse field Γ. For all T > 0, the system is disordered when studied at large length scales, so that it is in the paramagnetic phase (PM). At T = 0, the ground state seems disordered for Γ > Γc (from the point of view of the computational basis). For Γ < Γc, we encounter the spin-glass phase (SG), which is different for every disorder realization (equation (1)). b, Our finite-size scaling analysis (see, for example, refs. 48,60) of the critical point at T = 0 and Γ = Γc, in terms of the parameter k that represents Γ in the Trotter–Suzuki formulation (‘The Trotter–Suzuki formula’ in Methods; k grows as Γ decreases). Left, correlation length ξ(3) in units of the lattice size L versus k. The curves for the different L’s intersect at the critical point kc ≈ 0.29. Right, data in the left-hand panel of b, when represented as a function of the scaling variable L1/ν(k − kc) with 1/ν = 0.7, converge to a limiting curve as L grows. Points in b are statistical averages, and errors are one standard deviation. Our data set is fully described in Extended Data Table 1.

Source Data

Full size image

A strong theoretical command of the phase transition at Γc is clearly necessary. A very powerful tool in the analytical study of phase transitions is the renormalization group, which helps to clarify which properties of the critical point are not modified by the different microscopic details of different experiments. Only very broad features, such as symmetries, matter (making it possible to classify problems into universality classes). In fact, the study of disordered systems was one of the early applications of the renormalization group (see, for example, refs. 26,27,28), a strategy that is firmly established for D = 2 (ref. 29). Yet, it has taken considerable time and effort to show that the renormalization group—and the accompanying universality—also applies to disordered systems for D > 2 (refs. 30,31,32,33,34; even for D = 2, this was a hard endeavour for spin glasses35). Unfortunately, the study of the quantum spin-glass transition at finite D is considerably behind its thermal counterpart. Essentially, only the D = 1 case is well understood36,37,38,39.

The second-simplest problem to analyse, a spin glass with D = 2, poses quite a challenge. Indeed, different approaches have produced mutually contradictory predictions for the crucial physical quantity that ultimately determines whether the quantum computational complexity of the problem to be considered is smaller than its classical counterpart. We are referring to the energy gap Δ that separates the ground state from the first excited state of the Hamiltonian (1). Indeed, the required annealing time grows with 1/Δ2 (ref. 40). In a spin glass with N ∝ L2 quantum spins at Γ = Γc, where L is the linear size of the system, Δ ∝ L−z, where z is the so-called dynamic critical exponent. Early Monte Carlo simulations5,6,8 and a series-expansion study9 found finite values of z (for example, z ≈ 1.5 for D = 2 spin glasses5). A finite z is also a crucial assumption of the droplet model for the quantum spin-glass transition41. On the other hand, a real-space renormalization group analysis concluded that z = ∞ for D = 2 or 3 spatial dimensions10. (A Monte Carlo simulation too claimed that z = ∞ for D = 2 (ref. 42)).

The starting assumption of refs. 5,6,8,23 was a finite value of exponent z. Yet, to clarify the aforementioned controversy, our analysis is completely agnostic about z. Just as Rieger and Young pushed to the very limit the computational capabilities of the time using special hardware (transputers)5, we have performed unprecedented large-scale simulations on graphical processing units (GPUs) using highly tuned custom codes. A careful consideration of the global spin-flip symmetry, implemented by the parity operator \(P={\prod }_{x}{\sigma }_{x}^{X}\), turns out to be crucial. Although the gap for same-parity excitations scales algebraically in the number of spins, the gap for parity-changing excitations does close super-algebraically. The coexistence of these two qualitatively different scalings at the critical point is extremely unusual and is probably caused by Griffiths–McCoy singularities36,37,38,39,43. (In the spin-glass phase, instead, the spontaneously broken parity symmetry naturally generates an exponentially small gap but only for parity-changing excitations.) Although Griffiths–McCoy singularities are strongly dependent on D, interestingly, ref. 44 hinted at an algebraically scaling gap for same-parity excitations at the quantum critical point of a three-regular graph (a D = ∞ problem; parity-changing excitations were not studied in ref. 44).

The ground state

Our aim here is to study the phase transition as seen from the ground state (so that the spectra of excitations and, hence, exponent z do not play any role in the analysis in ‘The phase transition’). This entails taking the limit T → 0.

In the Trotter formulation that we use (‘The Trotter–Suzuki formula’ in Methods), the original quantum spins on an L × L lattice are replaced by classical spins on an L × L × Lτ lattice, Sx,τ = ±1. The extra dimension τ is named the Euclidean time. Γ is replaced by a new parameter k that grows as Γ decreases. The energy gap Δ translates into a correlation length η = 1/(kΔ) over Euclidean time. In this formulation, the limits T → 0 and Lτ → ∞ are equivalent.

Although our main results stem from Monte Carlo simulations, a complementary exact-diagonalization effort on small systems (‘Exact diagonalization’ in the next section and ‘More about exact diagonalization’ in Methods) has been extremely useful, both in shaping our analysis and in providing an understanding for how the limit T → 0 is approached (‘At the limit of zero temperature’).

Exact diagonalization

The main lessons that exact diagonalization of systems with size L ≤ 6 (‘More about exact diagonalization’ in Methods and Extended Data Figs. 1 and 3) have taught us are the following.

The parity operator P splits the spectrum of the Hamiltonian (equation (1)) into even energy levels (E0,e < E1,e < …) and odd levels (E0,o < E1,o < …). The ground state is even and its energy is EGS = E0,e.

The first excited state is E0,o. The minimal gap Δ = E0,o − E0,e displays dramatic fluctuations among samples, up to the point that a statistical analysis should be conducted in terms of log Δ. Furthermore, log Δ varies notably with k. By contrast, the sample-to-sample fluctuations of the same-parity gaps, Δe ≡ E1,e − E0,e and Δo ≡ E1,o − E0,o are very mild (also their k-dependence is mild). For all our samples, Δe and Δo are of similar magnitude and, unless Δ turns out to be inordinately large, Δe, Δo ≫ Δ.

Thermal expectation values of even operators (operators \({\mathcal{A}}\) such that \({\mathcal{A}}=P\,{\mathcal{A}}P\)) reach their T = 0 limit for surprisingly small values of Lτ. The reasons for this benign behaviour are understood (‘The limit of zero temperature’ in Methods).

The phase transition

We turn now to Monte Carlo simulations. The standard spin-glass correlation function, when computed for the ground state, is afflicted by a very large anomalous dimension that makes the spin-glass susceptibility χ(2) barely divergent at the critical point5. We have circumvented this problem by considering instead the correlation matrix M (equation (9) in Methods and refs. 45,46,47). From M, one can compute not only χ(2) but also a better behaved susceptibility χ(3). The corresponding correlation length ξ(3) is suitable for a standard finite-size scaling study of the phase transition (see, for example, ref. 48), which is illustrated in Fig. 1b.

The analysis in ‘The critical point and critical exponents’ in Methods found that for the critical point kc, the correlation-length exponent ν and exponents γ(n) (\({\chi }^{(n)}({k}_{{\rm{c}}})\propto {L}^{{\gamma }^{(n)}/\nu }\); for a discussion see ‘One-time observables’ in Methods):

$${k}_{{\rm{c}}}=0.2905(5),\quad \frac{1}{\nu }=0.71(24)(9),$$
(2)
$$\frac{{\gamma }^{(2)}}{\nu }=0.27(8)(8),\quad \frac{{\gamma }^{(3)}}{\nu }=1.39(23)(11).$$
(3)

The first error estimate is statistical whereas the second error accounts for systematic effects. Note that the bound ν ≥ 2/D (ref. 49) is verified and that \({\chi }^{(2)}\propto {L}^{\frac{{\gamma }^{(2)}}{\nu }\approx 0.3}\), indeed, barely diverges5.

At the limit of zero temperature

The naive approach to the limit T → 0 would be to study a fixed set of samples for a sequence of growing Euclidean lengths Lτ and to check when the results become independent of Lτ (indeed, T > 0 effects die out as \({{\rm{e}}}^{-{L}_{\tau }/\eta }\)). However, according to ‘Exact diagonalization’ above, this is just wishful thinking. Indeed, some instances have an inordinately small gap (and, hence, a huge Euclidean correlation length η), and so \({{\rm{e}}}^{-{L}_{\tau }/\eta }\approx 1\) for all values of Lτ that we can simulate (one would like \({{\rm{e}}}^{-{L}_{\tau }/\eta }\ll 1\), instead).

Fortunately, considering simultaneously periodic boundary conditions (PBCs) and antiperiodic boundary conditions (APBCs) over Euclidean time offers a way out. A detailed analysis (‘The limit of zero temperature’ in Methods) shows that the sequence of results for growing Lτ converges to T = 0 from opposite sides. As Lτ grows (Fig. 2c), the PBC sequence monotonically decreases, whereas the APBC sequence increases. Thus, the statistical compatibility of both types of boundary condition ensures that the zero-temperature limit has been reached (within our statistical errors).

Fig. 2: Ensuring that the zero-temperature limit has been reached by comparing PBCs and APBCs over Euclidean time.
figure 2

a, Correlation length ξ(3) (‘One-time observables’ in Methods) versus k, as computed for our largest systems with L = 24 and Lτ = 2,048 and with both PBCs and APBCs for the same set of 1,280 samples. The statistical agreement for PBCs and APBCs indicates that the T → 0 limit has been effectively reached for this quantity. b, As a for the Binder cumulant (‘One-time observables’ in Methods). The dashed line represents the critical point, kc ≈ 0.29. c, The even correlation functions Q2(τ) (‘Two-times observables’ in Methods), as computed for a single sample of L = 20 at k = 0.29, rather quickly reach their large-τ plateau. The functions depend on both Lτ and the boundary conditions. The PBC plateau decreases upon increasing Lτ, whereas the APBC plateau notably increases. The reason behind the stronger sensitivity of Lτ for APBCs is understood (‘The limit of zero temperature’ in Methods). Points in a, b, and c are statistical averages, and errors are one standard deviation. Our data set is fully described in Extended Data Table 1.

Source Data

Full size image

Spectra of excitations at the transition

The main tool for investigating excitations is the τ dependence of the Euclidean correlation function of several operators (‘Two-times observables’ in Methods). It is crucial to distinguish even operators (such as \(P\,{\mathcal{A}}P={\mathcal{A}}\)) from odd operators (\(P\,{\mathcal{A}}P=-{\mathcal{A}}\)). For even operators, the decay with τ is sensitive only to same-parity gaps (such as Δe and Δo defined in ‘Exact diagonalization’ above). Instead, odd operators feel the different-parity energy gap Δ.

For both symmetry sectors, the correlation functions computed for a sample decay exponentially (to zero for odd operators or, Fig. 2c, to a plateau for even operators). In both cases, the correlation lengths η and the energy gaps of appropriate symmetry are related as η = 1/(kΔ). Therefore, what the average over samples of an Euclidean correlation function really features is the probability distribution function (as computed over the different samples) of the correlation lengths η.

From now on, we focus on the critical point at kc ≈ 0.29. The uncertainty δkc about the value of the critical point, approximately 5 × 10−4 in our estimate, introduces systematic errors. However, finite-size scaling theory (see, for example, ref. 48) tells us that we can quantify these systematic effects by considering the scaling combination \(\delta {k}_{{\rm{c}}}{L}_{\max }^{1/\nu }\) where \({L}_{\max }=24\) is the largest lattice size we simulated. This scaling combination corresponds to a value of approximately 0.044. The smallness of this number and the smooth dependence of the scaling functions on the scaling variable \((k-{k}_{{\rm{c}}}){L}^{1/\nu }\) (see, for example, Fig. 1b, right) make it apparent that this source of errors is under control in our simulations. Furthermore, the excellent agreement between our results for same-parity excitations (see below) with a recent experiment (‘Conclusions and outlook’) reinforces this conclusion.

Even operators

This case is of utmost relevance because only even excited states can cause the system to leave its ground state in (ideal) quantum annealing for the Hamiltonian (equation (1)). Our approach is not entirely satisfying in this respect because, for a given sample, we obtain the smallest of the two same-parity gaps Δe and Δo (one would like to study only Δe). Fortunately, both gaps have similar magnitudes (Extended Data Fig. 3).

The first optimistic indication comes from the (subtracted) correlation function in Fig. 3a, which goes to zero (within errors) for a moderate value of τ. Indeed, the empirical distribution function for the correlation length ηe in Fig. 3b indicates mild sample-to-sample fluctuations and a relatively weak dependence on L. In fact, as shown Fig. 3b, for all L > 12, the probability distribution function turns out to depend on the scaling variable

$$u=\frac{{\eta }_{{\rm{e}}}-{\eta }_{{\rm{e}}}^{0}}{{L}^{{z}_{{\rm{e}}}}},\quad {\eta }_{{\rm{e}}}^{0}=2.2(3),\quad {z}_{{\rm{e}}}=2.46(17).$$
(4)

(Setting \({\eta }_{{\rm{e}}}^{0}=0\), the whole curve cannot be made to scale and the resulting estimate ze ≈ 1.7 is lower, see also the concluding paragraph of ‘Two-times observables’ in Methods). Thus, as anticipated, we conclude that the even symmetry sector shows algebraic scaling for its gap.

Fig. 3: Studying the spectra of even excitations at the critical point.
figure 3

a, Sample-averaged subtracted correlation function Q2,s(τ) (‘Fitting process and estimating the Euclidean correlation length’ in Methods) becomes compatible with zero for moderate values of τ, for all our system sizes. b, Left, after computing the Euclidean correlation length \({\eta }_{{\rm{e}}}^{(s)}\) for each sample, we computed for each L the empirical distribution function F(ηe), namely the probability F of finding a sample with \({\eta }_{{\rm{e}}}^{(s)} < {\eta }_{{\rm{e}}}\) (note the horizontal error bars). Right, the data in the left-hand panel of b, when plotted as a function of the scaling variable u (equation (4)) do not show any L residual L dependence other than for our smallest sizes L = 8 and 12. Points in a and b are statistical averages, and errors are one standard deviation. Our data set is fully described in Extended Data Table 1.

Source Data

Full size image

Odd operators

As would be expected from the exact results for D = 1 (refs. 36,38,39) and the approximate renormalization group study for D = 2 (ref. 10), the odd correlation function \(\overline{C(\tau )}\) shown in Fig. 4a has, for L = 24, a power-law decay \(\overline{C(\tau )}\propto 1/{\tau }^{\widetilde{b}}\) with \(\widetilde{b}=0.6\). This implies that the magnetic susceptibility—the linear response to a magnetic field aligned with the z axis—diverges at the critical point. Indeed, the susceptibility diverges if \(\widetilde{b} < 1\) (because it is twice the integral of \(\overline{C(\tau )}\) for τ going from 0 to ∞).

Fig. 4: Studying the spectra of odd operators at the critical point.
figure 4

a, The decay of the sample-averaged correlation function C(τ) (equation (14)) approaches a power law as L increases. The dashed line is a guide for the eyes. Indeed, we needed to represent C(τ) in terms of \(\widetilde{\tau }=({L}_{\tau }/{\rm{\pi }})\sin ({\rm{\pi }}\tau /{L}_{\tau })\) to avoid distortions due to the PBCs (\(\widetilde{\tau }\) and τ are almost identical for small τ/Lτ). b, Empirical distribution function F(η) as a function of \(\log \eta \) for all our system sizes. Note that we can compute only up to some L-dependent F because our largest Lτ is not large enough to allow for a safe determination of η in some samples. c, For large η, the asymptotic behaviour \(F(\eta )=1-B/{\eta }^{b}\) is evinced by the linear behaviour (in logarithmic scale) of 1 ��� F as a function of η. We fond b ≈ 0.8. The dashed line is a guide for the eyes. Points in a, b, and c are statistical averages, and errors are one standard deviation. Our data set is fully described in Extended Data Table 1.

Source Data

Full size image

Furthermore, \(\widetilde{b} < 1\) also for k < kc (Extended Data Fig. 6). We, therefore, conclude that the susceptibility diverges in the paramagnetic phase. This is exactly the same behaviour encountered for D = 1. Accordingly, the probability distribution function of the Euclidean correlation length η (recall that η = 1/(kΔ)) has an extremely fat tail (Fig. 4b,c). This is the behaviour of a Lévy flight, which strongly suggests that the scenario of an infinite-randomness fixed point10 is, indeed, realized for D = 2 Ising spin glass.

Conclusions and outlook

We have solved a decades-long controversy through an extreme-scale simulation on GPUs and a careful consideration of the main symmetries of the problem. Our main conclusion is very optimistic: there is no reason in principle why a quantum annealer cannot remain in the ground state while entering the spin-glass phase (Fig. 1a). However, as discussed below, this is not quite the same as solving our optimization problem. To adiabatically enter the spin-glass phase, the annealing time would just need to grow as a power law with the number of quantum spins (equation (4)), provided that parity-changing excitations are avoided (something that, at least nominally, is within the capabilities of current hardware). Let us emphasize that this protective mechanism is crucial for making it possible to study the dynamics at the quantum critical point, as done in recent experiments23,50. Entering the spin-glass phase adiabatically is also essential for getting good approximate solutions to optimization problems. Universality and the renormalization group suggest that our optimistic conclusion extends to a vast family of problems (all problems that share the spatial dimensions and the basic symmetries as our spin glass on a square lattice, which includes the interaction graph of D-wave devices).

One may worry that unavoidable experimental effects, for example, tiny but random longitudinal magnetic fields, could break the parity symmetry. Fortunately, the mechanism discussed in the previous paragraph still protects the annealer. Indeed, the gaps internal to the parity sectors, Δe and Δo, are similar and much larger than the parity-mixing gap Δ. Hence, undesired parity errors will mostly mix the ground states in the even and odd sectors (which are essentially equivalent, as far as the solution to the classical optimization problem is concerned).

We should mention as well that the combination of critical exponents zν has been measured in a very recent experiment50, which produced results in excellent agreement with those reported in this work.

However, our findings pose, as well, many questions. Let us list a few.

We have seen that entering the spin-glass phase with a quantum annealer should be doable with an effort polynomial in the number of qubits. However, to solve an optimization problem, one still needs to go adiabatically all the way from the critical point to Γ = 0. This is a difficult journey, at least for problems with D → ∞ spatial dimensions (refs. 44,51,52). However, it has been recently announced that an algebraic speed-up, compared to classical algorithms, is within reach23. As we said above, having a finite exponent z is a basic prerequisite also for an algebraic speed-up.

We know that D = 2 Ising glasses may be both hard and easy to solve on a classical computer. For instance, problems formulated for a square lattice with nearest-neighbour interactions can be solved quite efficiently (see, for example, ref. 53). However, adding second-neighbour interactions results in an NP-complete problem. As far as we know, it is still unclear whether or not these two problems belong to the same (quantum) computational complexity class. As the second-neighbour interactions should play no role at Γc, differences between the two kind of problems (if any) should arise for transverse fields Γ < Γc.

In this work, we have chosen problem instances with uniform probability, but this is not a necessity. One could focus, instead, on samples that are particularly hard to solve with a classic digital computer54,55,56. It would be interesting to test whether the benign scaling in equation (4) remains unchanged under these challenging circumstances. We know that these classically hard problems are even harder to solve for a D-wave annealer57, but there are many possible explanations for this poor performance of quantum hardware (see, for example, refs. 58,59).

Another possible avenue of research is concerned with three-dimensional systems. A recent experiment conducted on D-wave hardware suggests z ≈ 1.3 (ref. 23). Whether this finite dynamic exponent is for only even excitations (as would be the case in two dimensions) or it is unrestricted6 is, probably, worthy of investigation.

Methods

Model and simulations

Our quantum spins occupy the nodes of a square lattice of side L, endowed with PBCs. The coupling matrix Jx,y in equation (1) is non-vanishing only for nearest lattice neighbours. A problem instance, or sample, is defined by the choice of the Jx,y matrix. The non-vanishing matrix elements, Jx,y = Jy,x are random, independent variables in our simulations. Specifically, we chose Jx,y = ± J0 with 50% probability. We chose energy units such that J0 = 1.

Given an observable \({\mathcal{A}}\), we shall refer to its thermal expectation value in a given sample as \(\langle {\mathcal{A}}\rangle \) (see equation (8) below; the temperature is as low as possible, ideally T = 0). Thermal expectation values are averaged over the choice of couplings (quenched disorder; see, for example, ref. 25). We shall denote the second average (over disorder) as \(\overline{\langle {\mathcal{A}}\rangle }\).

Crucial symmetries

The most prominent symmetries in this problem are the gauge and the parity symmetries. Both symmetries are exact for the Hamiltonian (equation (1)) and for its Trotter���Suzuki approximation (‘The Trotter–Suzuki formula’).

The parity symmetry \(P=\prod _{x}{\sigma }_{x}^{X}\) is a self-adjoint, unitary operator that commutes with the Hamiltonian (equation (1)), as well as with the exact (equation (7)) and approximate (equation (8)) transfer matrices. The Hilbert space is divided into two subspaces of the same dimension according to the parity eigenvalues, which are either +1 (even states) or −1 (odd states). We also classify operators as either even (\(P\,{\mathcal{A}}P={\mathcal{A}}\)) or odd (\(P\,{\mathcal{A}}P=-{\mathcal{A}}\)). Matrix elements of even operators are non-vanishing only if the two states have the same parity (in contrast, the parity of the states should differ for odd operators). An oversimplified but enlightening cartoon of the spectra in our problem is represented in Extended Data Fig. 1 (see below some exact-diagonalization results that support this view).

Parity symmetry is just a particular type of gauge transformation. Let us arbitrarily choose for each site nx = 0 or 1. The corresponding gauge operator \({G}_{\{{n}_{x}\}}=\prod _{x}{({\sigma }_{x}^{X})}^{{n}_{x}}\) is self-adjoint and unitary. It transforms the Hamiltonian in equation (1) into a Hamiltonian of the same type but with modified couplings61: \({J}_{x,y}\to {J}_{x,y}{(-1)}^{{n}_{x}+{n}_{y}}.\) The gauge symmetry is enforced by the process of taking the disorder average. Indeed, the gauge-transformed coupling matrix has the same probability as the original one. Hence, meaningful observables should be invariant under an arbitrary gauge transformation. The parity operator is obtained by setting nx = 1 for all sites, which does not modify Jx,y (hence, parity is a symmetry for a given problem instance not just a symmetry induced by the disorder average).

The Trotter–Suzuki formula

We follow the Trotter–Suzuki approximation62,63, which replaces the original quantum spins on an L × L lattice by classical spins on an L × L × Lτ lattice, Sx,τ = ±1. The extra dimension τ is named Euclidean time. We shall write S as a shorthand for the LDLτ spins in the system (D = 2, here). Instead, Sτ will refer to the LD spins at time τ. We take equal-strength couplings along the time and space directions (see, for example, ref. 5). The probability of S is given by

$$p({\bf{S}})=\frac{{{\rm{e}}}^{-k{\mathcal{E}}({\bf{S}})}}{Z},\quad Z=\sum _{\{{\bf{S}}\}}{{\rm{e}}}^{-k{\mathcal{E}}({\bf{S}})},$$
(5)

where Z is the partition function, and

$$\begin{array}{l}{\mathcal{E}}({\bf{S}})\,=\,-\mathop{\sum }\limits_{\tau =0}^{{L}_{\tau -1}}\left[\frac{1}{2}\sum _{x,y}\,{J}_{x,y}{S}_{x,\tau }{S}_{y,\tau }+\sum _{x}{S}_{x,\tau }{S}_{x,\tau +1}\right],\\ \,\,\,\,\varGamma \,=\,-\frac{1}{2k}\log \,\tanh \,k.\end{array}$$
(6)

Although we have not done so, note that one could use a different coupling over Euclidean time, kτ, which would relate k and Γ as \(\varGamma =-\frac{1}{2k}\log \,\tanh \,{k}_{\tau }\) or, equivalently, \({k}_{\tau }=-\frac{1}{2}\log \,\tanh (k\varGamma )\) (the continuum-time limit at T = 0 would be reached by taking Lτ → ∞ first and, afterwards, k → 0). PBCs are assumed over Euclidean time. Below, we shall find it useful to consider as well APBCs along only the τ direction. Besides, as the reader may check, k is a monotonically decreasing function of Γ.

Possibly the most direct connection between the D + 1 classical spin system and the original quantum problem is provided by the transfer matrix \(\widetilde{{\mathcal{T}}}\) (refs. 64,65). Let us define \({{\mathcal{H}}}_{0}=-\frac{1}{2}\sum _{x,y}{J}_{x,y}{\sigma }_{x}^{Z}{\sigma }_{y}^{Z}\) and \({{\mathcal{H}}}_{1}=-\varGamma \sum _{x}{\sigma }_{x}^{X}\). The quantum thermal expectation value at temperature T = 1/(kLτ) is

$$\langle \langle {\mathcal{A}}\rangle \rangle =\frac{{\rm{Tr}}{\mathcal{A}}\,{\widetilde{{\mathcal{T}}}}^{{L}_{\tau }}}{{\rm{Tr}}{\widetilde{{\mathcal{T}}}}^{{L}_{\tau }}},\quad \widetilde{{\mathcal{T}}}={e}^{-k({{\mathcal{H}}}_{0}+{{\mathcal{H}}}_{1})}.$$
(7)

Now, for \({\mathcal{A}}={A}_{{\rm{cl}}}(\{{\sigma }_{{\bf{x}}}^{Z}\})\), \({A}_{{\rm{cl}}}\) being an arbitrary function, the Trotter–Suzuki approximation amounts to substituting the true transfer matrix in equation (7) by its proxy \({\mathcal{T}}\) (\({\mathcal{T}}=\widetilde{{\mathcal{T}}}+{\mathcal{O}}({k}^{3})\)):

$$\langle {\mathcal{A}}\rangle =\frac{{\rm{Tr}}{\mathcal{A}}{{\mathcal{T}}}^{{L}_{\tau }}}{{\rm{Tr}}{{\mathcal{T}}}^{{L}_{\tau }}},\quad {\mathcal{T}}={e}^{-k{{\mathcal{H}}}_{0}/2}{e}^{-k{{\mathcal{H}}}_{1}}{e}^{-k{{\mathcal{H}}}_{0}/2}.$$
(8)

\(\langle {\mathcal{A}}\rangle \) can be computed as well by averaging Acl(Sτ), evaluated over configurations distributed according to equation (5). (The value of τ is arbitrary; hence, one may gain statistics by averaging over τ).

Finally, let us emphasize that both \({\mathcal{T}}\) and \(\widetilde{{\mathcal{T}}}\) are self-adjoint, positive-definite, transfer matrices that share the crucial symmetries discussed in ‘Crucial symmetries’.

Observables

The quantities defined in ‘One-time observables’ were aimed at probing the ground state as k (and, hence, Γ (equation (6))) varies. These quantities will always be averaged over disorder before we proceed with the analysis.

Instead, the time correlations in ‘Two-times observables’ will probe the excitations. These time correlations will be analysed individually for each sample (sample-to-sample fluctuations are considered in ‘Fitting process and estimating the Euclidean correlation length’).

One-time observables

We consider the LD × LD correlation matrices M and \(\widehat{M}\) (refs. 45,46) with p = (2π/L, 0) or (0, 2π/L):

$${M}_{x,y}=\langle {\sigma }_{x}^{Z}{\sigma }_{y}^{Z}\rangle ,\quad {[\widehat{M}]}_{x,y}={M}_{x,y}{e}^{\text{i}{\bf{p}}\cdot (x-y)}.$$
(9)

The n-body spin-glass susceptibilities at both zero and minimal momentum are

$${\chi }^{(n)}=\frac{\overline{{\rm{Tr}}\left[{M}^{n}\right]}}{{L}^{D}},\quad {F}^{(n)}=\frac{1}{{L}^{D}}\overline{{\rm{Tr}}\left[\widehat{M}{M}^{n-1}\right]}.$$
(10)

χ(n) and F(n) give us access to the second-moment correlation length (see, for example, ref. 48):

$${\xi }^{(n)}=\frac{1}{2\sin ({\rm{\pi }}/L)}\sqrt{\frac{{\chi }^{(n)}}{{F}^{(n)}}-1}.$$
(11)

As L grows, χ(n) and ξ(n) remain of order 1 in the paramagnetic phase, whereas, in the critical region, they diverge as \({\chi }^{(n)}\propto {L}^{{\gamma }^{(n)}/\nu }\) and ξ(n) ∝ L. In the spin-glass phase, χ(n) ∝ LD(n−1) (ξ(n) ∝ La with some unknown exponent a > 1).

Our χ(n=2) and ξ(n=2) are just the standard quantities in the spin-glass literature66,67. In fact, in the simplest approximation (see ref. 47 for a more paused exposition) at criticality and for large separations r between x and y, Mx,y ∝ vxvy/ra with vx, vy of order 1 (so, γ(n)/ν = (n − 1)D − na in this approximation). Hence, if D > a, γ(n) grows with n. Indeed, n = 3 turns out to be a good compromise between statistical errors, which grow with n, and a strong enough critical divergence of χ(n) (χ(n=2) barely diverges5).

Besides, we computed the Binder cumulant using Q2 = LDχ(n=2) as

$$B=\frac{{Q}_{4}}{{Q}_{2}^{2}},\quad {Q}_{4}=\sum _{x,y,z,u}\overline{{\langle {\sigma }_{x}^{Z}{\sigma }_{y}^{Z}{\sigma }_{z}^{Z}{\sigma }_{u}^{Z}\rangle }^{2}}.$$
(12)

The Gaussian nature of the fluctuations in the paramagnetic phase causes B to approach 3 as L grows for fixed k < kc. B reaches different large-L limits for fixed k ≥ kc (for k > kc, different behaviours are possible, depending on the degree of replica symmetry breaking68).

Two-times observables

Let us start by defining the time-correlation function of an observable \({\mathcal{A}}\) (for simplicity, consider a product of σZ operators at some sites):

$${C}_{{\mathcal{A}}}(\tau )=\frac{{\rm{Tr}}\,{\mathcal{A}}\,{{\mathcal{T}}}^{\tau }\,{\mathcal{A}}\,{{\mathcal{T}}}^{{L}_{\tau }-\tau }}{{\rm{Tr}}\,{{\mathcal{T}}}^{{L}_{\tau }}}.$$
(13)

\({C}_{{\mathcal{A}}}(\tau )\) can be computed from our spin configurations distributed according to the classical weight (equation (5)) by averaging \({\sum }_{{\tau }_{1}=0}^{{L}_{\tau }-1}{A}_{{\rm{cl}}}({{\bf{S}}}_{{\tau }_{1}}){A}_{{\rm{cl}}}({{\bf{S}}}_{{\tau }_{1}+\tau })/{L}_{\tau }^{2}\) (notation as in ‘The Trotter–Suzuki formula’).

Specifically, we have considered

$$C(\tau )=\frac{\sum _{x}{C}_{{\sigma }_{x}^{Z}}(\tau )}{{L}_{\tau }^{D}},\quad {Q}_{2}(\tau )=\frac{\sum _{x,y}{C}_{{\sigma }_{x}^{Z}{\sigma }_{y}^{Z}}(\tau )}{{L}_{\tau }^{2D}}.$$
(14)

Let us briefly recall some general results64,65 for \({C}_{{\mathcal{A}}}(\tau )\) that follow from the spectral decomposition of the transfer matrix (to simplify the notation, let us first disregard the parity symmetry and consider PBCs).

The τ-dependence is given by the additive contribution of every pair of states En < Em (n = 0 is the ground state). Each pair generates an exponentially decaying term \({B}_{n,m}[{{\rm{e}}}^{-\tau /{\eta }_{n,m}}+{{\rm{e}}}^{-({L}_{\tau }-\tau )/{\eta }_{n,m}}]\) with correlation length ηn,m = 1/(kΔn,m), where Δn,m = (Em − En). The amplitude is \({B}_{n,m}={{\rm{e}}}^{-{L}_{\tau }/{\eta }_{0,n}}\,| \langle n| A| m\rangle {| }^{2}/\widehat{Z}\), with \(\widehat{Z}=1+\sum _{n\ > \ 0}{{\rm{e}}}^{-{L}_{\tau }/{\eta }_{0,n}}\). Hence if Lτ ≫ η0,n the contribution of this pair of states can be neglected. Besides, in the presence of parity symmetry, for even \({\mathcal{A}}\) we find Bn,m = 0 if the parity of \(\left|n\right\rangle \) and \(\left|m\right\rangle \) differ (for odd operators, Bn,m = 0 if both parities are equal). This is why the largest correlation length for Q2(τ) is the maximum of ηe and ηo, whereas the relevant correlation length for C(τ) is η (Extended Data Fig. 1).

Moreover, for even operators, every state \(\left|n\right\rangle \) provides an additive contribution to a τ-independent term (namely, the plateau in Extended Data Fig. 2): \({{\rm{e}}}^{-{L}_{\tau }/{\eta }_{0,n}}\,| \langle n| A| n\rangle {| }^{2}/\widehat{Z}\). Instead, for odd operators, \(| \langle n| A| n\rangle | =0\) (hence, odd operators lack a plateau). To manage a situation with APBCs, one just needs to add an extra parity operator as a final factor in both the numerator and the denominator of both equations (8) and (13). If parity is a symmetry, as is the case for our problem, \(\widehat{Z}\) is modified as \(\widehat{Z}=1+\sum _{n > 0}{p}_{n}{{\rm{e}}}^{-{L}_{\tau }/{\eta }_{0,n}}\) (pn = ±1 is the parity of the state) and the contribution to the APBC plateau gets an extra factor pn, as well.

One may wish to average over samples \({C}_{{\mathcal{A}}}(\tau )\). The dominant time decay for a given sample will be approximately Be−τ/η. Hence, the time decay for the averaged \(\overline{{C}_{{\mathcal{A}}}(\tau )}\) is an integral ∫dBdηρ(B, η)Be−τ/η, where ρ(B, η) is the corresponding probability density (over the samples). For simplicity, let us assume that fluctuations of the amplitude B are mild. Then, the scaling in equation (4) implies that, for large Lτ, the asymptotic time decay of \(\overline{{C}_{{\mathcal{A}}}(\tau )}\) is a function of the scaled time τ/Lz, where z is a dynamic exponent that applies to the parity of A. One just needs to change the integration variable as u = η/Lz and recall the scaling form ρ(η) ≈ f(η/Lz)/Lz, where f is a suitable scaling function.

The limit of zero temperature

We shall assume that we can reach Lτ large enough so that \({{\rm{e}}}^{-{L}_{\tau }/{\eta }_{{\rm{e}}}},\;{{\rm{e}}}^{-{L}_{\tau }/{\eta }_{{\rm{o}}}}\ll 1\) (the notation is explained in Extended Data Fig. 1). Moreover, we shall not assume that \({\epsilon }\equiv {{\rm{e}}}^{-{L}_{\tau }/\eta }\) is small (in fact, for some samples, one could even have ϵ ≈ 1).

Now, consider an even operator \({\mathcal{A}}\), and let us define \({{\mathcal{A}}}_{{\rm{e}}}=\langle {0}_{{\rm{e}}}| {\mathcal{A}}| {0}_{{\rm{e}}}\rangle \) and \({{\mathcal{A}}}_{{\rm{o}}}=\langle {0}_{{\rm{o}}}| {\mathcal{A}}| {0}_{{\rm{o}}}\rangle \) (the thermal expectation value at exactly T = 0 is \({{\mathcal{A}}}_{{\rm{e}}}\)). The plateau at τ ≫ ηe, ηo (Extended Data Fig. 2) is

$${C}_{{\mathcal{A}}}(\tau \gg {\eta }_{{\rm{e}}},{\eta }_{{\rm{o}}})={{\mathcal{A}}}_{{\rm{e}}}^{2}+[{{\mathcal{A}}}_{{\rm{o}}}^{2}-{{\mathcal{A}}}_{{\rm{e}}}^{2}]\frac{\zeta {\epsilon }}{1+\zeta {\epsilon }},$$
(15)

where ζ = 1 for PBCs and ζ = −1 for APBCs. Thus, we get for the plateau of Q2(τ)

$${Q}_{2}(\tau \gg {\eta }_{{\rm{e}}},{\eta }_{{\rm{o}}})={Q}_{2,{\rm{e}}}+[{Q}_{2,{\rm{o}}}-{Q}_{2,{\rm{e}}}]\frac{\zeta {\epsilon }}{1+\zeta {\epsilon }},$$
(16)

where Q2,e and Q2,o are, respectively, the average over all pairs (x, y) of \({{\mathcal{A}}}_{{\rm{e}}}\) and \({{\mathcal{A}}}_{{\rm{o}}}\) (\({\mathcal{A}}={\sigma }_{x}^{Z}{\sigma }_{y}^{Z}\); recall equation (14)). Let us give a few hints about the derivation of equations (15) and (16). The contribution of state \(\left|n\right\rangle \) to the plateau is \({{\rm{e}}}^{-{L}_{\tau }/{\eta }_{0,n}}\,| \langle n| A| n\rangle {| }^{2}{p}_{n}/\widehat{Z}\) where pn = 1 for PBCs whereas, for APBCs, pn = 1 for even states and pn = −1 for odd states. As explained before, we just keep the states \(\left|{0}_{{\rm{e}}}\right\rangle \) and \(\left|{0}_{{\rm{o}}}\right\rangle \) when estimating the plateau.

To excellent numerical accuracy, the left-hand side of equation (16) is also the value one gets for \({\rm{Tr}}{M}^{2}/{L}^{2D}\) (Extended Data Fig. 2). In fact, the difference between \({\langle {\mathcal{A}}\rangle }^{2}\) and its plateau is \(\zeta {\epsilon }{({{\mathcal{A}}}_{{\rm{e}}}-{{\mathcal{A}}}_{{\rm{o}}})}^{2}/{(1+\zeta {\epsilon })}^{2}\) (hence, quadratic in \(({{\mathcal{A}}}_{{\rm{e}}}-{{\mathcal{A}}}_{{\rm{o}}})\) rather than linear, as in equation (15)).

Now, despite their simplicity, two important consequences follow from equations (15) and (16).

First, the limit T → 0 (or Lτ → ∞) is approached monotonically. Furthermore, the systems with PBCs and APBCs (Extended Data Fig. 2) approach the limit from opposite sides. We have explicitly checked all our samples, finding no instance where the APBC plateau lies above the PBC one (it is intuitively natural to expect that the PBC system will be more ordered than the APBC one). Hence, we conclude that the samples with PBCs converge to T → 0 from above, whereas the APBC ones converge from below.

Second, as 0 < Q2(τ) < 1 also for APBCs, one has \(0 < {Q}_{2}^{{\rm{APBC}}}(\tau \gg {\eta }_{{\rm{e}}},{\eta }_{{\rm{o}}})\)\( < {Q}_{2,{\rm{e}}} < 1\). Thus, \(| {Q}_{2}^{{\rm{APBC}}}(\tau \gg {\eta }_{{\rm{e}}},{\eta }_{{\rm{o}}})-{Q}_{2,{\rm{e}}}| < 1\), and we conclude that ∣Q2,o − Q2,e∣ < (1 − ϵ)/ϵ. Hence, quite paradoxically, the particularly difficult samples with ϵ ≈ 1 generate a very small finite-temperature bias in the PBC estimator (compare the Lτ dependence of the PBC and the APBC plateaus in Extended Data Fig. 2). This is why we are able to reach the T → 0 limit for the even operators, even if a fraction of our samples suffers from a large value of ϵ.

Simulation details

We followed two approaches: exact diagonalization of the transfer matrix (equation (8)) and Markov chain Monte Carlo simulations of the classical weight (equation (5)). GPUs were crucial for both. We provide here only the main details (the interested reader is referred to ref. 69).

Exact diagonalization is limited to small systems (up to L = 6 in our case). Indeed, the transfer matrix has a size \({2}^{{L}^{2}}\times {2}^{{L}^{2}}\). Parity symmetry has allowed us to represent \({\mathcal{T}}\) as a direct sum of two submatrices of half that size69. Specifically, we computed the eigenvalues \({{\rm{e}}}^{-k{E}_{0,{\rm{e}}}}\), \({{\rm{e}}}^{-k{E}_{1,{\rm{e}}}}\), \({{\rm{e}}}^{-k{E}_{0,{\rm{o}}}}\) and \({{\rm{e}}}^{-k{E}_{1,{\rm{o}}}}\), as well as the corresponding eigenvectors \(| {0}_{{\rm{e}}}\rangle \), \(| {0}_{{\rm{o}}}\rangle \), \(\left|{1}_{{\rm{e}}}\right\rangle \) and \(\left|{1}_{{\rm{o}}}\right\rangle \), for 1,280 samples of L = 6 at k = 0.31 and 0.305 (the same samples at both k values). We repeated the calculations for a subset of 350 samples at k = 0.3 and 0.295. We managed to keep the computing time within an acceptable time frame of 20 min per diagonalization using 256 GPUs, thanks to a highly tuned custom matrix-vector product69. These computations have proven to be invaluable in the process of taking the limit Lτ → ∞ (‘More about exact diagonalization’ in Methods).

Our Monte Carlo simulations used the parallel tempering algorithm70, implemented over the k parameter in equation (5), to ensure equilibration. We equilibrated 1,280 samples for each lattice size (Extended Data Table 1). As a rule, we estimated errors using the bootstrap method71, as applied to the disordered average.

We simulated six real replicas of every sample (six statistically independent simulations of the system), for several reasons. Replicas allowed us to implement the equilibration tests based on the tempering dynamics56. They also provided unbiased estimators of products of thermal averages (equation (10)). Finally, fluctuations between replicas allowed us to estimate errors for the time-correlation functions (equation (14)), as computed in a single sample (‘Fitting process and estimating the Euclidean correlation length’).

The Monte Carlo code exploits a three-level parallelization (multispin coding, domain decomposition and parallel tempering), which kept the spin-update time below 0.5 ps (ref. 69), competitive with dedicated hardware22.

More about exact diagonalization

The schematic representation of the spectrum in Extended Data Fig. 1 is based on the distribution functions in Extended Data Fig. 3 (we typically compute the inverse distribution function; see ‘Fitting process and estimating the Euclidean correlation length’ for details).

Indeed, the correlation length η displays very large sample-to-sample fluctuations (to the point that a logarithmic representation is advisable) and a very strong k-dependence (Extended Data Fig. 3, left). In contrast, ηe is always a number of order one in our L = 6 samples (Extended Data Fig. 3, middle). Furthermore, ηo/ηe ≈ 1 in all cases (Extended Data Fig. 3, right).

In fact, the distribution for η is a Lévy flight (Fig. 4c; for large η, \(F(\eta )=1-B/{\eta }^{b}\)). The mechanism allowing exponent b to vary with k (hence, with transverse field (equation (6))) is sketched in Extended Data Fig. 4. Let us compare the value of η for the same sample at k1 and k2 (k1 < k2). With great accuracy, \(\eta ({k}_{2})=\alpha {[\eta ({k}_{1})]}^{1+\beta }\), where α and β are constants (for fixed k1 and k2) and β > 0. Thus, ordering samples according to their η(k1) is the same as ordering by η(k2), because one is a monotonically increasing function of the other. Hence, the same sample occupies percentile F in the distribution for both k1 and k2. It follows that b(k2) = b(k1)/(1 + β) for the exponent characterizing the Lévy flight. In other words, because b(k2) < b(k1), the tail at large η becomes heavier as k increases (see ‘On the magnetic susceptibilities’ for an extended discussion).

The critical point and critical exponents

After taking care of the Lτ → ∞ limit (within errors) in our study of the phase transition, we still need to cope with the finite spatial dimension L. We shall do so using finite-size scaling72,73,74,75 (Fig. 1c). The main questions we shall address are the computation of the critical exponents and the estimation of the critical point. Our main tool will be the quotients method48,60,76, which, surprisingly, keeps our two questions somewhat separate.

The quotients method starts by comparing a dimensionless quantity at two sizes La < Lb (in our case, ξ(3)/L as a function of k). First, we located a coupling k*(La, Lb) such that the curves for La and Lb cross (Fig. 1b). Now, for dimensionful quantities A, scaling in the thermodynamic limit as \({\xi }^{{x}_{A}/\nu }\), we consider the quotient \({Q}_{A}={A}_{{L}_{a}}/{A}_{{L}_{b}}\) at k*(La, Lb). Barring scaling corrections, \({Q}_{A}={({L}_{a}/{L}_{b})}^{{x}_{A}/\nu }\), which yields an effective estimate of xA/ν. Indeed, considering only the leading correction to the scaling exponent ω, we have for the effective exponent:

$${\left.\frac{{x}_{A}}{\nu }\right|}_{{L}_{a},{L}_{b}}=\frac{{x}_{A}}{\nu }\,+\,\frac{1}{\log ({L}_{b}/{L}_{a})}\log \frac{1+{D}_{A}{L}_{b}^{-\omega }}{1+{D}_{A}{L}_{a}^{-\omega }},$$
(17)

where DA is an amplitude. Our estimates for the effective exponents can be found in Extended Data Table 2. Yet, effective exponents need to be extrapolated to the thermodynamic limit through equation (17). Unfortunately, we have not been able to estimate exponent ω, as there were two difficulties. First, the range of L values at our disposal was small. Second the analytic background48 for the r = 2 observables and for the Binder parameter (‘One-time observables’) compete with the L−ω corrections. Hence, we followed an alternative strategy. We fitted our effective exponents to equation (17) with fixed ω (the fitting parameters were the extrapolated xA/ν and the amplitude DA). To account for our ignorance about ω, we made it vary in a wide range 0.5 ≤ ω ≤ 2. The central values in equations (2) and (3) were obtained with ω = 1, whereas the second error estimate accounts for the ω-dependence of xA/ν. Indeed, the first error estimate is the statistical error as computed for ω = 1, whereas the second error estimate is the semi-difference between the extrapolations to infinite size obtained with ω = 2.0 and ω = 0.5. To take into account the data correlation, we employed a bootstrap method77. We considered only the diagonal part of the covariance matrix in the fits and performed a new fit for every bootstrap realization. Errors were computed from the fluctuations of the fitting parameters. Fortunately, the systematic errors turned out to be comparable (for 1/ν smaller) with the statistical ones.

Like the critical point, we expected scaling corrections of the form k*(La, Lb) = kc + DkF(La, Lb), where Dk is an amplitude78:

$$F({L}_{a},{L}_{b})={L}_{a}^{-\left(\omega +\frac{1}{\nu }\right)}\frac{1-{s}^{-\omega }}{{s}^{1/\nu }-1},\quad s=\frac{{L}_{b}}{{L}_{a}}\,.$$
(18)

Unfortunately, this result is not of much use without a ω estimate. Fortunately, see Extended Data Table 2 and Extended Data Fig. 5, the values of k*(La, Lb) obtained from ξ(3)/L seem not to depend on size. In fact, our estimate for kc in equation (2) is an interval that encompasses all our results (the shaded area in Extended Data Fig. 5). Furthermore, the crossing points for B and ξ(2)/L (Extended Data Fig. 5) seem also reasonably well represented by equation (18).

Fitting process and estimating the Euclidean correlation length

Our aim in this section is to determine the relevant correlation lengths for C(τ) and Q2(τ) at a fixed k, for our NS = 1,280 samples. The results are characterized through their empirical distribution function (Figs. 3 and 4). Given that NS is large, we needed an automated approach.

The first step was estimating, for a given sample, C(τ) and Q2(τ), as well as their standard errors, by using our six replicas. Now, the analysis of a noisy correlation function (such as C(τ) and Q2(τ); see, for example, Extended Data Fig. 2) needs a fitting window79,80. We chose the window upper limit as \({\tau }_{{\rm{w}},f}\equiv \mathop{\min }\limits_{\tau }\{\tau | \,f(\tau )=3.5{\sigma }_{f(\tau )}\}\), with f(τ) either C(τ) or Q2,s(τ) = Q2(τ) − Q2,pl, where Q2,pl is the plateau (Extended Data Fig. 2) and σf(τ) is the corresponding standard error. We faced two problems. First, for odd C(τ), some samples have τw,C ≥ Lτ/2. For these samples, η > Lτ, and hence, it was impossible to estimate them (Fig. 4b). Q2,s(τ) was not affected by this problem (Fig. 3). Second, we needed to estimate the plateau Q2,pl. To do so, we fitted Q2(τ) for τ ∈ [Lτ/4, Lτ/2] to a constant Q2,pl. In the few exceptions where this fit was not acceptable (as determined by its figure of merit χ2/degrees of freedom (dof) computed with the diagonal part of the covariance matrix), we proceeded as explained below (we used \({\tau }_{{\rm{w}},{Q}_{2}}={L}_{\tau }/2\) in those cases).

We determined the correlation lengths through fits to \(C(\tau )\,=\)\(B[{\text{e}}^{-\tau /\eta }+{\text{e}}^{(\tau -{L}_{\tau })/\eta }]\) and \({Q}_{2}(\tau )={Q}_{2,\text{pl}}+{B}_{\text{e}}{e}^{-\tau /{\eta }_{\text{e}}}+{B}_{\text{o}}{e}^{-\tau /{\eta }_{\text{o}}}\). The fitting parameters were the amplitudes and the correlation lengths (and, for the above-mentioned exceptional samples, also Q2,pl). Actually, for Q2(τ) we considered fits with one and with two exponential terms, keeping the fit with the smallest χ2/dof, as we could not determine which of the two correlation lengths obtained in the fit corresponded to the even gap (Extended Data Fig. 1). Hereafter, ηe is the larger of the two. As for the lowest limit of the fitting window, we started from \({\tau }_{\min ,{Q}_{2}}=1\) and \({\tau }_{\min ,C}={\tau }_{{\rm{w}},C}/10\), and we kept increasing the corresponding \({\tau }_{\min }\) until χ2/dof went below 0.5 for Q2 (below 1 for C(τ)).

Finally, we determined the empirical distribution function for the correlation lengths. Let X be either \(\log \eta \) or ηe (see below for some subtleties regarding η). We actually computed the inverse function X[F] by sorting in increasing order the NS values of X and setting \(X[F=i/{N}_{{\rm{S}}}]\) as the ith item in the ordered list. We obtained X[F] at the value of k of interest through linear interpolation of X[F] computed at the two nearest values of k in the parallel tempering grid. To estimate the errors in X[F], we employed a bootstrap method with 10,000 as the resampling value. In each resampling, we randomly picked NS values of X. For the chosen sample, we extracted X from a normal distribution centred in X as obtained from the fit. The standard deviation was the fitting error for X.

For \(X=\log \eta \), we needed to cope with the problem that we could determine X for only NOK of our NS samples. We decided to determine X[F] only up to \({F}_{{\rm{safe}}}\equiv ({N}_{{\rm{OK}}}-4\sqrt{({N}_{{\rm{S}}}-{N}_{{\rm{OK}}}){N}_{{\rm{OK}}}/{N}_{{\rm{S}}}})/{N}_{{\rm{S}}}\) (the maximum possible F minus four standard deviations). We imposed for every bootstrap resampling that X could be obtained in at least FsafeNS samples (this limitation was irrelevant in practice).

Let us conclude by mentioning that the estimates in equation (4) were obtained through a joint fit for ηe[F], with F = 0.5, 0.6, 0.7, 0.8 or 0.9. Errors were estimated as explained in ‘The critical point and critical exponents’.

On the magnetic susceptibilities

The sample-averaged linear susceptibility to an external magnetic field at T = 0, \({\chi }_{{\rm{lin}}}^{(h)}\), can diverge only if \(\overline{C(\tau )}\) decays slowly for large τ (because \({\chi }_{{\rm{lin}}}^{(h)}=1+2{\sum }_{\tau =1}^{\infty }\overline{C(\tau )}\); Extended Data Fig. 6). Yet, the periodicity induced by the PBCs (Extended Data Fig. 6a) made it difficult to study the behaviour at large τ. Fortunately, representing \(\overline{C(\tau )}\) as a function of \(\widetilde{\tau }=\frac{{L}_{\tau }}{{\rm{\pi }}}\sin ({\rm{\pi }}\tau /{L}_{\tau })=\tau [1+{\mathcal{O}}({\tau }^{2}/{L}_{\tau }^{2})]\) greatly alleviated this problem (Extended Data Fig. 6b). Thus armed, we could study the long-time decay of \(C(\tau )\propto 1/{\widetilde{\tau }}^{\widetilde{b}}\) as a function of k (Extended Data Fig. 6c). Indeed, \(\widetilde{b}\) decreased as k increased. As C(τ) ≈ Be−τ/η for any sample, the mechanism discussed in ‘More about exact diagonalization’ in Methods is clearly at play. The heavy tail of F(η) became heavier as k increased, which resulted in a decreasing exponent \(\widetilde{b}\). In fact, the critical exponent \(\widetilde{b}=1\) was encountered at k ≈ 0.285, well into the paramagnetic phase (\({\chi }_{{\rm{lin}}}^{(h)}=\infty \) if \(\widetilde{b}\le 1\)).

The Lévy-flight perspective provides a simple explanation for the results in refs. 6,41. In a single sample, the different susceptibilities to a magnetic field (linear, third-order and so on) are proportional to increasing powers of η. Hence, the existence of the disorder average of a given (generalized) susceptibility boils down to the existence of the corresponding moment of the distribution F(η). As soon as F(η) decays for large η as a power law, some disorder-averaged susceptibility will diverge (probably a higher-order one). Lower-order susceptibilities diverge at larger values of k. Hence, it is not advisable to use this approach to locate the critical point.