Ergodic and localized regions in quantum spin glasses on the Bethe lattice

TL;DR

This paper investigates the quantum dynamics of a transverse field Ising spin glass on the Bethe lattice, identifying a many-body localized phase at low transverse fields within the spin glass region, with implications for quantum algorithms.

Contribution

It reveals the existence of a many-body localized phase within the spin glass phase on the Bethe lattice and discusses its impact on quantum adiabatic algorithms.

Findings

Existence of a many-body localized region at low transverse fields.

Quantum dynamics split into localized and delocalized regions within the spin glass phase.

Implications for the performance of quantum adiabatic algorithms.

Abstract

By considering the quantum dynamics of a transverse field Ising spin glass on the Bethe lattice we find the existence of a many body localized region at small transverse field and low temperature. The region is located within the thermodynamic spin glass phase. Accordingly, we conjecture that quantum dynamics inside the glassy region is split in a small MBL and a large delocalized (but not necessarily ergodic) region. This has implications for the analysis of the performance of quantum adiabatic algorithms.

Figures

Equations

$$H=-\sum _{⟨i,j⟩}{J}_{ij}{\sigma }_i^z{\sigma }_j^z-\Gamma \sum _{i=1}^N{\sigma }_i^x,$$

$$q_{EA}\equiv \frac{1}{N}\sum _{i=1}^N⟨{\sigma }_i^z{⟩}^2,$$

$$C_{ij}\equiv ⟨{\sigma }_i^z{\sigma }_j^z⟩-⟨{\sigma }_i^z⟩⟨{\sigma }_j^z⟩,$$

$$\mathrm{Sp}(H_0)=\{0,2\Gamma ,4\Gamma ,\dots ,2N\Gamma \}.$$

$$|0⟩\equiv \underset{i=1}{\overset{N}{⨂}}|\to {⟩}_i,$$

$$|i\rangle\equiv\sigma_{i}^{z}|0\rangle=|\leftarrow\rangle_{i}\otimes\Big{(}\bigotimes_{j\neq i}|\rightarrow\rangle_{j}\Big{)},$$

$$H(J)=H_0+JV$$

$$E_0(J)=J{E}_0^{(1)}=0$$

$$E_n(J)=2\Gamma +J{E}_n^{(1)}$$

$$H_{\mathrm{hop}}=-\sum_{\langle i,j\rangle}J_{ij}\Big{(}|i\rangle\langle j|+|j\rangle\langle i|\Big{)}.$$

$$G_i={\left(E-\sum _{j\in \partial i}{J}_{ij}^2{G}_j^{(c)}\right)}^{-1},$$

$$G_j^{(c)}={\left(E-\sum _{k\in \partial j}{J}_{jk}^2{G}_k^{(c^{\prime })}\right)}^{-1},$$

$${\rho }_i(E)=\frac{1}{\pi }\underset{\Im E\to 0^{+}}{\lim }\Im G_i(E).$$

$$G^{(c)}={\left(E-K{J}^2{G}^{(c)}\right)}^{-1},$$

$$G^{(c)}=\frac{2}{E-\sqrt{E^2-4K{J}^2}}.$$

$$G_i^{-1}$$

$${\Sigma }_i=\frac{2(K+1)J^2}{E-\sqrt{E^2-4K{J}^2}}$$

$$H(\Gamma )=H_0+\Gamma V.$$

$$\psi (x)=⟨x|\psi (\Gamma )⟩$$

$$H=-\sum _a{E}_a|a⟩⟨a|-\Gamma \sum _{a,a^{\prime }}⟨a|V|a^{\prime }⟩|a⟩⟨a^{\prime }|.$$

$$\psi_{r}\equiv\max\Big{\{}\lvert\langle w|\psi\rangle\rvert:w\in\{\pm 1\}^{N},\operatorname{dist}(v,w)=r\Big{\}}\quad\quad r=1,2,\ldots$$

$$P(N,r,\xi)\equiv Pr\Big{(}\lvert\psi_{r}\rvert^{2}\leq Ce^{-r/\xi}\Big{)}.$$

$$\underset{r\to \infty }{\lim }\underset{N\to \infty }{\lim }P(N,r,\xi )=1.$$

$$P(N,r,\xi)=P\Big{(}\frac{\ln\lvert\psi_{r}\lvert^{2}}{2r}\leq-\frac{1}{2\xi}+O\Big{(}\frac{1}{r}\Big{)}\Big{)},$$

$$Z_N\equiv \frac{\ln |{\psi }_N|}{N}\to \ln (\Gamma /{\Gamma }_c)\text{as N→∞},$$

$${\psi }_{\alpha }(b)\approx \sum _{p\in \text{spaths}(a,b)}\prod _{i\in p}\frac{\Gamma }{E_a-E_i},$$

$$v_t=D\cdot A\cdot v_{t-1}$$

$$D_{bc}$$

$$A_{bc}$$

$$(v_0{)}_c$$

$${\psi }_{\alpha }(b)=(v_N{)}_b,$$

Full text

1 Introduction

\setstcolor

blue

In recent years the study of the different dynamical regimes of isolated quantum systems has received a lot of attention, due to improved experimental techniques [1, 2] and theoretical progress. In the latter, one can identify two distinguished but not independent lines of research: the first is the study of how ergodicity is realized in isolated quantum systems, a mechanism that goes under the name of eigenstate thermalization hypothesis [3]; The second is the study of the most typical mechanism for failure of ergodicity in presence of quenched disorder [4, 5, 8, 9, 10] (although some authors have suggested disorder in the initial state suffices [6, 7]) named many-body-localization. While dynamical phases satisfying ETH can be described by the usual tools of statistical mechanics and thermodynamics, MBL systems behave a lot like integrable systems [15] with local integrals of motions [9, 11, 12, 13, 14]: Transport is suppressed [16, 17], entanglement entropy grows slowly to its thermodynamic value [18], and some symmetry breaking phases can exist, protected from the Mermin-Wagner theorem, in low dimensions and at high temperature [19].

Soon after its inception, it has been pointed out that MBL phases can be detrimental [20] for the performance of Adiabatic Quantum Computation (AQC) protocol introduced in [21] (see also [22]). This has been contested in later works [23] and it remains a controversial claim. Since this protocol has proven to be the most promising for the realization of a quantum computer [24], sorting out this question is of paramount importance for both theoretical discussions and technological implications.

In a series of recent works, which involve one of the present authors [25, 26, 27], the question of the appearance of an MBL phase in some models of quantum spin glasses has been addressed with the result that, for realistic, mean-field glasses MBL can exist only in finite-connectivity models, while in fully-connected models only a weaker form, remnant of the clustering phase existing in the phase space of the classical model [28], exists. These earlier works point at the necessity to examine a finite-connectivity quantum spin glass, in search of MBL.

In this work we set to do exactly this: we analyze the quantum dynamics of an isolated quantum spin glass showing that there is an ETH phase at large transverse field, possibly extending all the way down to the spin-glass phase while at all system sizes we were able to study, an MBL phase exists for small transverse field. Therefore there should be an MBL-ETH dynamical transition in between.

2 Thermodynamics

The focus of the present paper is the transverse-field Ising spin glass model defined on a regular random graph (RRG) of degree $d=3$, whose Hamiltonian is defined by

[TABLE]

where $\Gamma>0$ is the strength of the transverse field, the disordered interaction couplings $J_{ij}$ take either of the two values in $\{\pm 1\}$ with equal probability and $\sigma_{i}^{a}$ (for $a=x,z$) is a Pauli matrix acting on the $i$-th spin of the system. The sum of the terms $\langle i,j\rangle$ is taken over the edges of a $3$-regular reandom graph $\mathcal{G}$. Consequently, the Hamiltonian of Eq. 1 has two sources of disorder: the disordered interactions $J_{ij}$ and the topology of the underlying regular random graph $\mathcal{G}$.

Here and in what follows, a $d$-regular111the degree of a regular graph is also customarily defined by referring to its “connectivity” $K$. The relations between degree and connectivity is given by $K=d-1$. random graph of size $N$ is defined as graph uniformly sampled from the set $\mathbb{G}_{N}^{d}$ of all connected graphs with $N$ vertices and fixed degree $d$. These graphs are related to the $d$-regular Bethe lattice, the unique (up to isomorphism) connected tree graph of fixed degree $d$ with denumerably-many vertices. The Bethe lattice was introduced in [29] to define models where the Bethe-Peierls approximation is exact. Regular graphs can be seen as finite-size approximation to the Bethe lattice in the sense that even though they contain loops (and the Bethe lattice does not) the fraction of loops of any fixed length $\ell$ vanishes as the size $N$ approaches infinity (*i.e. *almost all loops are of length greater than $\ell$). Thus RRGs are “tree-like” in the (finite-sized) neighbourhood of any of their vertices, and are locally indistinguishable from the Bethe lattice (see Fig. 1). For this reason, regular random graphs in the $N\rightarrow\infty$ limit have been used to model the statistical properties of the Bethe lattice and approximation schemes related to the Bethe-Peierls method (such as the cavity method and belief-propagation algorithms) are expected to perform well when applied to systems with local interactions that define a regular random graph structure.

The thermodynamical properties of the Hamiltonian of Eq. (1) were studied in a series of papers [30, 31] that collectively reconstructed the equilibrium phase diagram shown in Fig. 2. For large values of the transverse field, or for high temperatures, the system lies in a paramagnetic phase where the $z$-component of each spin fluctuates randomly, so that the average $z$-magnetization is zero. In the region of small values of $\Gamma$ and low temperatures $T$, each spin freezes independently in the $z$-direction and the system enters a phase where the Edwards-Anderson order parameter,

[TABLE]

becomes strictly positive with a continuous transition [30, 33]. The classical $\Gamma=0$ line of the phase diagram, where the transition is driven solely by the temperature, was studied in [34] where it was found that the critical temperature was given by the formula $T_{c}=1/\tanh^{-1}(1/\sqrt{2})\approx 1.13$. The interior of the $(T,\Gamma)$-plane was studied in Ref. [30]. The authors developed a quantum version of the cavity method to explore numerically the thermodynamic limit of the model and compute the critical line. The physics of the quantum $T=0$ line was studied by a series of papers [30, 31] that gave various estimates of the critical point using cavity-like approximation schemes affected by the systematic error that these methods exhibit when applied to loopy graphs. It was conclusively studied numerically in [33] where the critical value of the transverse field was computed by path-integral Monte Carlo methods using several different physical quantities. All were found to agree on an estimated value of $\Gamma_{c}=1.82\pm 0.02$ with a finite-size correction that disappears as $1/N$ for $N\rightarrow\infty$. Zero-temperature properties of the ground state were also studied in [33] where the ground state was found to have volumetric entanglement – as measured by the Rényi entropy of order two – in the entire part of the $T=0$ line which lies in the paramagnetic phase. The equal-time spin-spin connected correlation function in the $z$ direction,

[TABLE]

was also studied in [33] in order to analyze the average and the extremal behaviour of the correlations across the transition. They considered the spatial maximal correlation $C_{\max}(r)$, defined by picking the maximal value of $C_{ij}$ among all spins $j$ that are at fixed distance $r$ from a given spin $i$ (and averaging over $i$), and the mean correlation $C_{\mathrm{mean}}(r)$, defined by taking the average value of $C_{ij}$ among all spins $j$ that are at fixed distance $r$ from $i$. Both were shown to follow a stretched exponential decay, and the correlation lengths extracted from them were found to converge to a finite value at the critical point of the transition in the thermodynamic limit.

3 Ergodic region at large transverse field

In order to develop some intuition about the phase transition, the authors of Ref. [33] did a perturbative expansion by using the modulus of the interaction strength $J=|J_{ij}|$ as the perturbative parameter. For ease of reference we summarize here their results. The “free” Hamiltonian is the transverse-field term $H_{0}\equiv-\Gamma\sum_{i}\sigma_{i}^{x}$. By shifting the ground state energy $E_{0}=-\Gamma N$ to zero, its spectrum is given by

[TABLE]

Its unique ground state $|0\rangle$ (where we denote $\sigma^{x}|\rightarrow\rangle=|\rightarrow\rangle$)

[TABLE]

is taken as a pseudovacuum of quasiparticles. The operators $\sigma_{i}^{z}$ create an excitation on top of the ground state

[TABLE]

which is interpreted as a state containing a quasiparticle at site $i$. These form the $N$-fold degenerate eigenspace of $H_{0}$ with energy $2\Gamma$. Additional applications of the $\{\sigma_{j}^{z}\}$ operators either move the state upwards in the spectrum, creating states with two, three, or more quasiparticles $|i,j\rangle,|i,j,k\rangle,\ldots$, or downwards, by annihilating existing quasiparticles.

The perturbed Hamiltonian is

[TABLE]

where $V$ is the dimensionless spin-glass term $V=-(1/J)\sum_{\langle i,j\rangle}J_{ij}\sigma_{i}^{z}\sigma_{j}^{z}$

First-order perturbation theory gives a null correction to the ground-state energy

[TABLE]

since $JE_{0}^{(1)}=\langle 0\rvert JV\lvert 0\rangle=0$. The degenerate band of one-particle states is split by the perturbation into distinct levels

[TABLE]

for $n=1,\ldots,N$, where $JE_{n}$ are eigenvalues and $\{|\phi_{n}\rangle\}$ are the eigenvectors of the operator $\tilde{V}$, which is the perturbation operator $V$ restricted to the (unperturbed) one-particle subspace. The direct computation of the $E_{n}^{(1)}$ requires the diagonalization of an $N$ by $N$ matrix which, we will show, is a hopping matrix on the RRG.

Let us start by noticing that the term $\sigma_{i}^{z}\sigma_{j}^{z}$ applied to a quasiparticle state $|i^{\prime},j^{\prime},k^{\prime},\ldots\rangle$ can affect it in one of three ways: (i) it can create a pair of adjacent quasiparticles provided the sites $i$ and $j$ are devoid of them, (ii) it can move a quasiparticle from site $i$ to site $j$ provided site $j$ is empty and site $i$ is occupied (or vice versa, swapping the role of $i$ and $j$), or (iii) it can annihilate a pair of adjacent quasiparticles sitting on sites $i$ and $j$.

Note also that in the one-particle subspace annihilation processes cannot happen, while creation processes map a state into the three-particle subspace, which is orthogonal to the one-particle subspace. Therefore, only the hopping processes give a contribution to Eq. (2). For any state $|\psi\rangle$ in the one-particle subspace, the action of $J\tilde{V}$ is then equivalent to that of $H_{\mathrm{hop}}$, the Hamiltonian of a particle hopping on the same graph $\mathcal{G}$, with disordered hopping constants $J_{ij}=\pm J$

[TABLE]

For contrast, consider the Hamiltonian H_{\mathrm{hop}}^{(\mathrm{hom})}=-J\sum_{\langle i,j\rangle}\Big{(}|i\rangle\langle j|+|j\rangle\langle i|\Big{)} of a particle hopping on the graph $\mathcal{G}$ with homogeneous hopping coefficients $J>0$. Even though solving for the spectrum of the latter Hamiltonian is difficult for any finite graph $\mathcal{G}$, in the thermodynamic limit, our RRG $\mathcal{G}$ is the Bethe lattice and so spectral properties of this model can be computed exactly using an iterative method (see *e.g. *[35, 36] and Fig. 3). In particular, its spectral density is known to be supported on the set $[-2J\sqrt{K},2J\sqrt{K}]$, where $K$ is the constant connectivity of the graph. We give the proof here: one starts by writing down iteration equations for the diagonal Green’s function $G_{i}\equiv\left\langle i\rvert(E-H_{\mathrm{hop}})^{-1}\lvert i\right\rangle$ at the site $i$ for generic complex $E$

[TABLE]

where the cavity Green’s function $G^{(c)}_{j}=\left\langle j\rvert(E-H^{(i)}_{\mathrm{hop}})^{-1}\lvert j\right\rangle$ is the Green’s function of the operator $H^{(i)}_{\mathrm{hop}}$ obtained from the Hamiltonian $H_{\mathrm{hop}}$ by removing the hopping terms associated to the the edges incident to the site $i$. On the Bethe lattice the removal of these terms splits the system into $K+1$ isomorphic disconnected components that can be considered independently. Each component is an infinite rooted tree with a branching factor of $K$. These trees are isomorphic to each of their infinite descending subtrees rooted at any of their vertices. By writing the iteration equation (3) for $G_{j}^{(c)}$ one gets

[TABLE]

where $G^{(c^{\prime})}_{k}$ is a second-step cavity Green’s functions obtained from the $G^{(c)}_{j}$ Hamiltonian by further removing the hopping terms associated to the edges incident to $j$. By solving (4) and plugging the result in (3) one recovers $G_{i}$ and from $G_{i}$ the spectral density

[TABLE]

One can find $P(\Im\Sigma_{i})$ by $\Sigma_{i}(E)=E-G_{i}^{-1}$ and taking the limit $\Im{E}\to 0$. As we know from [4], the distribution of $\Im\Sigma$ will tend to a delta function for delocalized states and to a long-tailed distribution for localized states. However we can sidestep all this procedure recognizing that in all the equations only $J_{ij}^{2}$ appears and therefore the case $J_{ij}=\pm J$ is identical to the case $J_{ij}=J$ constant, which has only delocalized states. The latter statement follows from the observation that, given the constant connectivity of the graph and the constant value of the hopping $J$ the distribution of $G^{(c)}$ must be a delta function centered on the solution of the deterministic equation (3)

[TABLE]

which gives

[TABLE]

By inserting (6) into (3) one gets

[TABLE]

so

[TABLE]

irrespective of $i$ and therefore all the eigenstates are delocalized.

By extending this reasoning to any sector with $n$ particles, as soon as $n$ is finite when $N\to\infty$ we can conclude that all such few-particle states are delocalized. As $n/N$ becomes appreciable so does the interaction between the particle, we shall conjecture that the introduction of a small interaction between the particles (which is of $O(J^{2}/\Gamma)$), irrespective of its attractive or repulsive nature, does not localize excitations (this is true even in the presence of bound states) and the phase remains delocalized. The value of $\Gamma/J$ where this breaks down we cannot predict without considering quantitatively the interaction between the particles.

Having showed that there is an ergodic region for large $\Gamma/J$ we will show in the following section that starting from the opposite limit $\Gamma/J\ll 1$ we do have a localized region. Therefore there must be at least one dynamical phase transition in between. We unfortunately are not able to answer the question whether this transition is unique or a crossover through a sequence of ergodicity breaking transitions.

4 Many-body localization at small transverse field

4.1 Localization and the Forward Approximation

Let us consider the Hamiltonian (1) for small values of $\Gamma$. Then $H=H(\Gamma)$ can be treated as a small perturbation of the spin glass term $H_{0}\equiv-\sum_{\langle i,j\rangle}J_{ij}\sigma_{i}^{z}\sigma_{j}^{z}$ with the transverse-field operator $V=-\sum_{i}\sigma_{i}^{x}$ acting as the perturbation:

[TABLE]

Note that $H_{0}=-\sum_{\langle i,j\rangle}J_{ij}\sigma_{i}^{z}\sigma_{j}^{z}$ is diagonal in the $N$-fold $\sigma^{z}$ product basis, whose elements we label by strings $a\in\{\pm 1\}^{N}$ through the usual identification $|\uparrow\rangle\equiv|1\rangle$ and $|\downarrow\rangle\equiv|-1\rangle$. Now and throughout this section we use Latin letters $a,b,c,\ldots,x,y,z$ to label energies and eigenstates of the unperturbed Hamiltonian $H_{0}$ while Greek letters $\alpha,\beta,\gamma,\ldots$ are reserved as labels for the energies and the eigenstates of the perturbed Hamiltonian, so that in the $\Gamma\rightarrow 0$ limit the perturbed Greek labels “converge” to their corresponding Latin label222while this convergence does not hold in general for any choice of eigenbasis for the degenerate Hamiltonian $H_{0}$ we will later show that in our case this is the correct choice., *e.g. *for the energies $E_{\alpha}\rightarrow E_{a}$ and for the energy eigenstates $|\psi_{\alpha}\rangle\rightarrow|a\rangle$.

The eigenstates $|\psi(\Gamma)\rangle$ of the Hamiltonian $H(\Gamma)$ can be chosen to depend continuously on the parameter $\Gamma$, and will converge to the eigenstates of $H_{0}$ in the limit $\Gamma\rightarrow 0$. Consequently, if we write the wavefunctions of the states $|\psi(\Gamma)\rangle$ in the unpertubed energy eigenbasis

[TABLE]

we expect that these will converge to Kronecker delta functions as $|\psi(\Gamma)\rangle$ changes into the corresponding unperturbed eigenstate $|a\rangle$; this means that the state $|\psi(\Gamma)\rangle$ will be localized for small $\Gamma$. In the previous section we saw that in the large $\Gamma$ limit these states are delocalized so we expect that there will be a value $\Gamma_{c}$ where transition between these two behaviours occours. In fact, if we write the Hamiltonian is the eigenbasis of $H_{0}$ we get

[TABLE]

This is related to the Anderson model that describes a particle hopping on a lattice while under the effect of a disordered potential. This model is known to enter a localized phase when the disordered term dominates the hopping term. In Eq. (7) the unperturbed energies $\{E_{a}\}$ are the equivalent of a disordered potential while the values $\Gamma\left\langle a\rvert V\lvert a^{\prime}\right\rangle$ are the hopping coefficients for the particle. Note that here the underlying geometry of the hopping is defined by the matrix elements $\left\langle a\rvert V\lvert a^{\prime}\right\rangle$: the sites $a,a^{\prime}$ are adjacent (*i.e. *the particle can hop directly from one to the other) if and only if $\left\langle a\rvert V\lvert a^{\prime}\right\rangle\neq 0$. As $\Gamma\rightarrow 0$ the hopping is suppressed and the system is likely to enter a disorder-induced localized phase. In the transverse-field Ising spin glass Hamiltonian (1) the perturbation $V=-\sum_{i}\sigma_{i}^{x}$ defines a hopping over the $N$-dimensional Boolean hypercube $\mathcal{B}_{N}=\{\pm 1\}^{N}$.

A disorder-driven localization/delocalization transition is usually studied as a function of two quantities: the energy density $\epsilon$ of the eigenstates considered and the strength $W$ of the disorder in the Hamiltonian. One usually finds that states at a fixed energy density $\epsilon$ change from delocalized to localized as the disorder strength $W$ is increased, with an energy-dependent critial value $W_{c}(\epsilon)$ that marks the boundary between these two behaviours. The set of critical points $W_{c}(\epsilon)$ for different values of $\epsilon$ define the “mobility edge” of the system.

In our setup we keep the strength of the disordered interactions fixed $W\equiv\lvert J_{ij}\rvert=1$ so the relative strength of the disorder with respect to the ordering term $-\Gamma\sum_{i}\sigma_{i}^{x}$ is controlled by the parameter $W/\Gamma=1/\Gamma$. The mobility edge will consequently be defined by the critical values $\Gamma_{c}(\epsilon)\equiv 1/W_{c}(\epsilon)$.

Let us better define the kind of localization we will study. First, we fix a value $\epsilon$ for the energy density and $\Gamma$ for the strength of the transverse field. For a given realization of disorder of the Hamiltonian $H(\Gamma)$ (1) defined on a system of size $N$ and a state $|\psi\rangle$ with localization center $|v\rangle$ and energy density $\epsilon$ we define

[TABLE]

where the distance is taken over the Boolean hypercube. Note that in order to compute $\psi_{r}$ from the amplitudes $\lvert\langle w|\psi\rangle\rvert$ one needs to know (or find out) the localization center of the state $|\psi\rangle$. Next, for fixed $N,r\in\mathbb{N}$ and real numbers $C,\xi>0$ we define the quantity $P(N,r,\xi)$ as the probability (taken over all disorder realizations of $H(\Gamma)$ of size $N$ and states of energy density $\epsilon$) that the random variable $\lvert\psi_{r}\rvert^{2}$ satisfies $\lvert\psi_{r}\rvert^{2}\leq Ce^{-r/\xi}$

[TABLE]

We say that a disordered system is localized if there exist a real number $\xi>0$ such that

[TABLE]

This means that the probability distribution $p(x)=\lvert\psi(x)\lvert^{2}$ is under an exponentially-decaying envelope, except possibly for a region of finite radius.

Notice that we can write equivalently

[TABLE]

so we can study the distribution of the values of the random variable $Z_{r}\equiv\ln\lvert\psi_{r}\lvert/r$. It was observed [37] that if the system is localized then in the limit $N\rightarrow\infty$ the random variable $Z_{N}$ is peaked around the value

[TABLE]

which gives the value of $\Gamma_{c}(\epsilon)$. In practice one can calculate $\psi$ with $\Gamma=1$ and get $\Gamma_{c}$. By the way the relation $\xi=1/\ln(\Gamma/\Gamma_{c})$ tells us the critical exponent for the $\xi$ divergence is 1, since $\xi\sim 1/|\Gamma-\Gamma_{c}|$.

We have seen how one can study the localization properties of disordered systems by looking at the wavefunction values $\psi(x)$. Of course this is not always a simple task, as it usually requires the diagonalization of a matrix whose size grows exponentially with the size of the system. However, in a perturbative setup such as the one we described we can use a technique known as the “forward approximation” [4, 37], namely, we neglect the renormalization of the free energy of the unperturbed eigenstates $\{|a\rangle\}$. This is known to give an underestimate of the critical hopping strength both in the many-body and Anderson localization, it is then useful for our purpose of proving the stability of the phase (basically the same strategy is used in [5, 13]).

The forward approximation states that for a perturbed energy eigenstate $|\psi_{\alpha}\rangle$ that converges to $|a\rangle$ at $\Gamma=0$, the value of the wavefunction $\psi_{\alpha}(b)=\langle b|\psi_{\alpha}\rangle$ can be approximated by a sum of contributions associated to the shortest paths connecting the sites $a$ and $b$ in the Boolean hypercube:

[TABLE]

where the set $\text{spaths}(a,b)$ contains the shortest paths from $a$ to $b$.

Notice that the sum over paths of Eq. (10) can be computed numerically using a transfer matrix technique, in which case one computes iteratively the vector

[TABLE]

with

[TABLE]

and

[TABLE]

where $N$ is the distance between $a$ and $b$ in the Boolean hypercube. One can decrease memory requirements by noting that the vector $v$ is very sparse during most of the computation, so at each step most entries of the transfer matrix $\mathcal{T}=DA$ are irrelevant. This is because repeated applications of $\mathcal{T}$ to the initial state define a diffusion process on the Boolean hypercube where at each step $t$ one needs to propagate only the amplitudes of the vertices exactly at distance $t$ from the initial vertex. This means that in practice one does not need to store in memory the entire transfer matrix $\mathcal{T}$, but instead a new transfer matrix $\mathcal{T}_{t}$ is defined at each step that only propagates amplitudes from vertices actually relevant for that single step of propagation. This requires storing only $(N-t)\binom{N}{t}$ non-zero entries instead of $N2^{N}$ of the full transfer matrix $\mathcal{T}=DA$. The vector $v_{t}$ need to store only $\binom{N}{t}$ entries.

4.2 Numerical results

In this section we apply the previously-described methods to the transverse-field Ising spin glass Hamiltonian (1). Here we immediately face an issue: the computation of the forward approximation is obstructed by the fact that the spin glass term $H_{0}=-\sum_{i}J_{ij}\sigma_{i}^{z}\sigma_{j}^{z}$ has highly degenerate energy levels. This gives rise to diverging terms in Eq. (10) when $E_{a}=E_{i}$. In order to avoid this problem we add a weak and random longitudinal field term $H_{\mathrm{long}}$ to the Hamiltonian $H_{0}$:

[TABLE]

where each $h_{i}$ is distributed uniformly in $(-h,h)$ with $h=0.001$ (this has to be $\ll 1/N$ but $\gg e^{-aN}$). This has the effect of splitting the degeneracies while introducing only a negligible effect in the energies of the configurations (and therefore in the amplitudes $\psi_{\alpha}(b)$) and on the amplitudes of transition between non degenerate states.

We compute the many-body mobility edge for the system in the following way. For each system size $N=18,20,22,24,26$ we randomly generate a suitable number of realizations of disorder and for each of these we generate a set of initial states $a=a_{1},a_{2},\ldots,a_{k}$, making sure that their energy densities $\epsilon_{a}=E_{a}/N$ are (approximately) uniformly distributed in the range allowed by the model. Each of these states $a$ is then propagated to its $\mathbb{Z}_{2}$-symmetric state $b=-a$ (global spin flip) using the forward approximation algorithm with fixed $\Gamma=1$ in order to compute

[TABLE]

Notice that the configuration $-a$ is the only configuration $b\in\{\pm 1\}^{N}$ to satisfy $\operatorname{dist}(a,b)=N$, ergo for the state $|\psi_{\alpha}\rangle$ we read from Eq. (8) that (after setting $r=N$)

[TABLE]

The results are then binned according to the energy density of the initial state $a$ and the average of the random variable $Z_{N}=\ln\lvert\psi_{N}\lvert/N$ value was taken for each bin, obtaining $Z_{N}(\epsilon)=\langle Z_{N}\rangle$.

Using the formula $\Gamma_{c}(\epsilon)=\exp(-\ln\lvert\psi_{N}(\epsilon)\rvert/N)=\exp(-Z_{N}(\epsilon))$ from Eq. (9) we obtain a plot of the MBL critical point as a function of the energy density $\epsilon$ shown in Fig. 4. Note that this is the energy density of the unperturbed eigenstates $|a\rangle$, while usually one would write $\Gamma_{c}$ as a function of the energy density of the perturbed eigenstates $|\psi_{\alpha}\rangle$. However, the perturbed energies $E_{\alpha}=E_{a}+O(\Gamma^{2})$ coincide with the unperturbed ones up to second-order corrections in $\Gamma$, which we neglect.

In order to plot the MBL critical line in the $(T,\Gamma)$-phase diagram and compare it to the boundary of the glassy phase, we have to compute the relation $T=T(\epsilon)$ between temperature and (disorder-averaged) energy density. We used standard Monte Carlo methods to extract the (thermal) average energy density of different realizations of disorder at various temperatures and fixed $\Gamma=0$, then we took the average over the results333We note that as one approaches the ground state energy, small difference of energies translate to (relatively) large difference of temperatures due to small values of the heat capacity in the low-temperature regime. In order to effectively control this effect one would require better precision in the M.C. energy estimation.. In order to better understand the low energy regime we studied the ground state of the unperturbed (*i.e. *$\Gamma=0$) model. For each size $N=18,20,22,24,26$ we generated a large number ($\geq 1000$) of instances and extracted one of the ground states by performing a thermal annealing (whose results were checked against an exact solver for the smaller sizes). For each ground state we computed the $\Gamma_{c}$ value using the forward approximation. The disorder-averaged results are show in Fig. 5. Extrapolations give a value of $\Gamma_{c}=0.67$ in the thermodynamic limit, which seems consistent with Fig. 4.

Finally, we plotted a finite-size ($N=26$) estimate of the MBL critical line and the line of the glassy transition in the $(T,\Gamma)$-phase diagram (Fig. 6). The MBL phase seems to be strictly contained in the glassy phase, therefore there is a region of the phase diagram where the system is both glassy and delocalized.

5 Conclusion

We studied the localization properties of the transverse-field Ising spin glass model on the $3$-regular random graph in the limit where the trasverse-field is weak compared to the disordered interactions. This model is known to exhibit a transition from a paramagnetic to a glassy phase at low temperatures and weak transverse-field. The classical Ising spin glass model is widely believed to capture the complicated combinatorial structure of general $NP$-hard computational problems while the zero-temperature, weak transverse-field regime describes the final stage of a quantum annealing protocol designed to find the ground-state energy of the Ising spin glass. Many-body localization has been argued to be an obstacle to efficient quantum annealing due to the presence of exponentially-closing gaps in the localized phase.

We computed numerically the many-body mobility edge of the system in the forward approximation, finding that the energy eigenstates of the system indeed localize for small values of the transverse field at finite system sizes. When plotted against the equibilbrium phase diagram of the model, we discovered that the localized region does not coincide with the glassy phase. In particular, evidence points to the fact that the glassy phase is partitioned into a delocalized region and a localized one. We conjecture that the glassy, delocalized region will exhibit the same clustering of eigenstates observed in [27] for the $p$-spin model, where the eigenstates were found to form clusters inside of which the energies are distributed according to Wigner-Dyson while the global distribution of the energy levels of the model is Poissonian.

Moreover, we expect that classical methods that exploit the fine-tuning of thermal relaxation (such as simulated annealing) will perform poorly in the entire glassy phase while quantum annealers will perform poorly only once localization sets in. Therefore we conjecture that in the glassy, delocalized region of the phase space quantum annealing algorithms can outperform any classical thermal annealing protocol.

A natural future direction outlined by our work would be to check whether the same localization/delocalization transition is present when the disordered term of the Hamiltonian encodes a real-life computational problems such as $3$-SAT. In the affirmative case, a detailed comparison of the performance of *e.g. *simulated annealing and quantum annealing (either simulated numerically or by an actual experiment) inside of the region that is both glassy and delocalized would help shed light on the realistic capabilities of quantum annealers over classical thermal annealing and other algorithms based on stochastic local optimization.

6 Data accessibility

All numerical data are accessible upon request to the corresponding author.

7 Competing interests

The authors declare they have no competing interests.

8 Authors’ contributions

Both authors conceived and designed the study, G.M. carried out the numerical simulations. Both authors contributed in analyzing the data and drafting the manuscript. All authors give final approval for publication.

9 Acknowledgements

A.S. would like to thank the Google QAI group in Los Angeles (CA, USA), where part of this work was done.

10 Funding statements

A.S. acknowledges financial support in the form of a Google Faculty Research Award.