Possible Ergodic-nonergodic regions in the quantum Sherrington-Kirkpatrick spin glass model and quantum annealing

TL;DR

This study investigates the ergodic and nonergodic regions in the quantum Sherrington-Kirkpatrick spin glass model, revealing how quantum fluctuations influence ergodicity and annealing efficiency at different temperatures.

Contribution

It identifies and characterizes the finite-temperature ergodic region in the quantum spin glass model using numerical methods, highlighting its impact on quantum annealing.

Findings

Existence of a low-temperature ergodic region dominated by quantum fluctuations.

System size independence of annealing time in the ergodic region.

Nonergodic behavior persists at high temperatures with non-vanishing order parameter distribution.

Abstract

We explore the behavior of order parameter distribution of quantum Sherrington-Kirkpatrick model in the spin glass phase using Monte Carlo technique for the effective Suzuki-Trotter Hamil- tonian at finite temperatures and that at zero temperature obtained using exact diagonalization method. Our numerical results indicate the existence of low but finite temperature quantum fluc- tuation dominated ergodic region along with the classical fluctuation dominated high temperature nonergodic region in the spin glass phase of the model. In the ergodic region, the order parameter distribution gets narrower around the most probable value of the order parameter as the system size increases. In the other region, the Parisi order distribution function has non-vanishing value every- where in thermodynamic limit, indicating nonergodicity. We also show, that the average annealing time for convergence (to a low energy level of the model; within a small error range) becomes system size independent for annealing down through the (quantum fluctuation dominated) ergodic region. It becomes strongly system size dependent for annealing through the nonergodic region. Possible finite size scaling type behavior for the extent of the ergodic region is also addressed.

Equations

$$H=H_0+H_I;H_0=-\sum _{i<j}{J}_{ij}{\sigma }_i^z{\sigma }_j^z;H_I=-\Gamma \sum _{i=1}^N{\sigma }_i^x.$$

$$H_{eff}$$

$$P(q)=\overline{\frac{1}{t_1}\sum _{t=t_0}^{t_0+t_1}\delta (q-q^{\alpha \beta }(t))},$$

$$P(|Q|)=\overline{\frac{1}{N}\sum _{i=1}^N\delta (|Q|-Q_i)}.$$

Full text

I Introduction

Considerable amount of investigations have been made in studying the nonergodic behavior sudip-binder of the spin glass phase of the classical Sherrington-Kirkpatrick (SK) spin glass model sudip-sk . The phenomenon of replica symmetry breaking, induced by nonergodicity, occurs due to the appearance of macroscopically high free-energy barriers separating the local minima. Such highly rugged nature of free-energy landscape in spin glass phase causes the system to get trapped into any one (locally) self-similar region of the configuration space. Consequently one gets a broad order parameter distribution (or the replica symmetry breaking) in the spin glass phase as suggested by Parisi sudip-parisi . In this case, along with the peak at any non-zero value of the order parameter, its distribution also contains a long tail extended to the zero value of the order parameter in the thermodynamic limit. This localization due to nonergodicity has been identified to be responsible for the NP hardness of equivalent optimization problems (see e.g., sudip-das ).

The situation seems to be quite different when the SK spin glass is placed under a transverse field. Due to the presence of the quantum fluctuations, the system is able to tunnel through the tall (but narrow) free-energy barriers sudip-ray ; sudip-bikas ; sudip-ttc-book ; sudip-lidar ; sudip-troyer ; sudip-Katzgraber , inducing ergodicity (or absence of replica symmetry breaking). Consequently one would expect a narrowly peaked order parameter distribution in quantum SK spin glass model in the thermodynamic limit sudip-ray . This ergodicity has been identified to be responsible (see e.g., sudip-bikas ; sudip-ttc-book ) for the success of quantum annealing.

We have studied the nature of order parameter distribution of transverse field SK spin glass at finite temperature using Monte Carlo simulation of the effective Suzuki-Trotter Hamiltonian and using the exact diagonalization technique at zero temperature. In this numerical study we tried to identify the possible ergodic spin glass phase (due to quantum tunneling) of the system. We find a low temperature region in the quantum SK system, where the tails of the order parameter distribution vanishes in thermodynamic limit, suggesting convergence of the order parameter distribution to be a peaked one around the most probable value. Although the system sizes we studied are not very large, we believe our study clearly indicates the existence of a low temperature ergodic region in the spin glass phase of this quantum SK model. On the other hand, in other (high temperature) part of the spin glass phase, the order parameter distribution appears to remain Parisi type sudip-young which indicates lack of ergodicity in this part of the spin glass phase. We have already identified sudip-sudip the quantum fluctuation dominated part of the spin glass phase boundary of this model, crossing over at finite temperature to the classical fluctuation dominated part (see also sudip-yao ). Here we find that the line separating the ergodic and the nonergodic regions pass through the zero temperature-zero transverse field point and the above mentioned quantum-classical crossover point on the phase boundary.

We also study the variation of the average annealing time in the finite temperature Suzuki-Trotter Hamiltonian dynamics for the model in both the ergodic and nonergodic regions. For annealing down to a fixed low temperature and low transverse field point through the (quantum fluctuation dominated) ergodic region, we find the average annealing time to be independent of system size. On the other hand the average annealing time is observed to grow strongly with the system size, when similar annealing is performed through the (classical fluctuation dominated) nonergodic region.

II Model

The Hamiltonian of the quantum SK spin glass model with $N$ Ising spins is given by (see e.g., sudip-ttc-book )

[TABLE]

Here $\sigma_{i}^{z}$, $\sigma_{i}^{x}$ are the $z$ and $x$ components of Pauli spin matrices respectively and $\Gamma$ denotes the transverse field. The Hamiltonian in Eq. (1) becomes the classical SK spin glass Hamiltonian ($H_{0}$) for zero value of the transverse field. The spin-spin couplings ($J_{ij}$) are distributed following Gaussian distribution \rho(J_{ij})=\Big{(}\frac{N}{2{\pi}J^{2}}\Big{)}^{\frac{1}{2}}\exp\Big{(}\frac{-NJ_{ij}^{2}}{2J}\Big{)}, where the mean and standard deviation of the distribution are zero and $J/\sqrt{N}$ respectively (see e.g., sudip-ttc-book ). In this work we take $J=1$. To perform Monte Carlo simulation at finite temperature we map the Hamiltonian (1) into an effective classical Hamiltonian $H_{eff}$ by using Suzuki-Trotter formalism (see e.g., sudip-bkc-book ):

[TABLE]

where $\sigma_{i}^{m}$ ($=\pm 1$) represents the $i$-th (classical) Ising spin in the $m$-th replica. We have an additional dimension (in Eq. 2), namely the Trotter dimension. Here $M$ denotes the total number of Trotter slices and $\beta$ is the inverse of temperature $T$. ${M\to\infty}$ as $T\to 0$.

III Monte Carlo Results

For finite temperature study, we perform Monte Carlo simulation on $H_{eff}$ to obtain the order parameter distribution in spin glass phase of our model. To obtain such distribution function we first allow the system to equilibrate with $t_{0}$ Monte Carlo steps and the thermal averaging is made over next $t_{1}$ time steps. In one Monte Carlo step we update all the spins of the system once. After the equilibration in each Monte Carlo step $t$ we calculate the replica overlap $q^{\alpha\beta}(t)$, which is defined as $q^{\alpha\beta}(t)=\frac{1}{NM}\sum_{i=1}^{N}\sum_{m=1}^{M}(\sigma_{i}^{m}(t))^{\alpha}(\sigma_{i}^{m}(t))^{\beta}$. Here $(\sigma_{i}^{m})^{\alpha}$ and $(\sigma_{i}^{m})^{\beta}$ denote the spins of two replicas (in the $m$-th Trotter slice) having identical set of $J_{ij}$’s. The order parameter distribution $P(q)$ can be obtained as

[TABLE]

where the overhead bar denotes the configuration average over several sets of $J_{ij}$’s. The order parameter $q$ is defined as $q=\frac{1}{MN}\sum_{m=1}^{M}\sum_{i=1}^{N}\overline{{\langle\sigma_{i}^{m}\rangle}^{2}}$, where $\langle..\rangle$ denotes the thermal average for a given configuration of disorder. From numerical data we compute the distribution function $P(q)$ for a given set of $T$ and $\Gamma$ by considering both area normalization and peak normalization (where peaks of the distributions are normalized).

In our simulation we work with the system sizes $N=60,120,180,240$ and the number of Trotter slices is $M=15$. We have found that the equilibrium time of the system is not identical throughout the entire region of $\Gamma-T$ plane. The equilibrium time of the system (for $60\leq N\leq 240$) is typically $\lesssim 10^{6}$ within the region $T<0.25$ and $\Gamma<0.40$, whereas it becomes $\lesssim 10^{5}$ for the rest of the spin glass phase region. We take $t_{1}=1.5\times 10^{5}$ for Monte Carlo averaging and the configuration average is made over $1000$ samples (configurations). Because of its symmetry we have determined the distribution of $|q|$ instead of $q$. We observe that the value of $P(|q|)$ for $q=0$ has a clear system size dependence. We extrapolate the values of $P(0)$ with $1/N$ to get the value of $P(0)$ for infinite system size. We also calculate the width $W$ at half maximum of the distribution function. The width $W$ is define as $W=|q_{2}-q_{1}|$ where the value of $P(|q|)$ becomes half of its maximum value at $q=q_{1}$, $q_{2}$. Again we extrapolate the values of $W$ with $1/N$. We observe two distinct behaviors of such extrapolated values of both $P(0)$ and $W$ in spin glass phase. In the low temperature (and high transverse field) case, we notice that the values of $P(0)$ and $W$ both go to zero in the large system size limit (see Fig. 1a). Such an observation indicates $P(|q|)$ would approaches to Gaussian form in thermodynamic limit, suggesting ergodic behavior of the system. In contrast, for the other case (high temperature case) we find that $P(0)$ has finite value even in thermodynamic limit (see Fig. 1b)). There seems to be no possibility of $P(|q|)$ to approach the Gaussian form of distribution for infinite system size limit. It indicates that the system remains nonergodic in this region of spin glass phase. To identify the ergodic and nonergodic regions in spin glass phase more accurately, we also study the behavior of the peak normalized order parameter distribution. From such study again we find in low temperature and high transverse field the values of $P(0)$ and $W$ (extrapolated with $1/N$) become zero in thermodynamic limit (see Fig. 2 (a, b)). Again from peak normalized order parameter distribution we find that for high temperature and low transverse field the extrapolated values of the tail and width of distribution remains non zero in infinite system size limit (see Fig. 3 (a, b)).

IV Zero-temperature diagonalization results

For zero temperature study of our model, we have investigated the distribution of the spin glass order parameter using an exact diagonalization technique. The diagonalization of the quantum spin glass has been performed using Lanczos algorithm lanczos_52 to obtain its ground state. In this case we have considered system sizes ($N$) upto $20$. The Hamiltonian in Eq. (1) can be written in the spin basis states which are indeed the eigenstates of the Hamiltonian $H_{0}$. After performing diagonalization, the $n$-th eigenstate of the Hamiltonian in Eq. (1) is found out as $|\psi_{n}\rangle~{}=\sum_{\alpha=0}^{2^{N-1}}a_{\alpha}^{n}|\varphi_{\alpha}\rangle$, where $a_{\alpha}^{n}=\langle\varphi_{\alpha}|\psi_{n}\rangle$ and $|\varphi_{\alpha}\rangle$ denote the eigenstates of the Hamiltonian $H_{0}$. As the consequence of our interest in zero temperature analysis, we shall here mainly focus on the ground state ($|\psi_{0}\rangle$) averaging of different quantities of interest. One can define the order parameter for this zero temperature system as $Q=(1/N)\sum_{i}\overline{\langle\psi_{0}|\sigma_{i}^{z}|\psi_{0}\rangle^{2}}=(1/N)\sum_{i}\overline{Q_{i}}$ (note that $Q$ here for $T=0$ differs from $q$ defined earlier for $T\neq 0$, using replica average) sudip-sudip . Here also, the overhead bar indicates the configuration averaging. $Q_{i}$ denotes the site-dependent local order parameter value. The distribution of the local order parameter is then represented by

[TABLE]

Similar to the case of finite temperature, we here also have investigated the behavior of $P(|Q|)$ in the spin glass phase at different values of $\Gamma$. The variation of $P(|Q|)$ as a function $|Q|$ at $\Gamma=0.3$ is shown in Fig. 4 for four different system sizes. It may be noted that in this case also we have plotted the distribution curves for different system sizes normalized to their maximum values as well as area normalization under the curves. From both the plots in Fig. 4, we observe that $P(|Q|)$ shows a peak at a finite value of $|Q|$ along with non-zero weight at $Q=0$. However, the value of $P(0)$ decreases with the increase of the system size (although one can still detect an upward rise of $P(|Q|)$ for lower values of $|Q|$). To get the behavior of $P(|Q|)$ in the thermodynamic limit, we have computed the value of $P(|Q|)$ for infinite size system for each $|Q|$ by plotting $P(|Q|)$ as a function of $1/N$. The extrapolation of $P(|Q|)$ for infinite size system is shown in the inset of both plots in Fig. 4 for $Q=0.0$ and $0.1$. In addition, we have also calculated the width ($W$) at half of maximum: $W=|Q_{2}-Q_{1}|$ where at $Q_{2}$ and $Q_{1}$ the value of $P(|Q|)$ is the half of its maximum value. We plot $W$ as a function of $1/N$ to get its extrapolated value for infinite size system (see Fig. 4). Finally we have also plotted $P(|Q|)$ as a function of $|Q|$ with the extrapolated value of $P(|Q|)$ for infinite system size (see Fig. 4). One can observe that $P(|Q|)$ curve for infinite system becomes narrower as compared to the cases of finite system size. On the other hand, due to the limitation of the maximum system size we could consider in our numerics, we are here not able to get $P(|Q|)$ curve for very large system sizes showing results consistent with delta function form. The effect of the limitation of the system size is also present in the plot of $W$ with $1/N$ since the extrapolated $W$ does not acquire strictly zero value here. However, we infer from our extrapolated numerical analysis (from results of the small system sizes) that eventually the $P(|Q|)$ curve would become a delta function at a finite values of $|Q|$ in thermodynamic limit. This would suggest the system to become ergodic in the spin glass phase at zero temperature with a definite spin glass order parameter value.

V Annealing through ergodic and nonergodic regions

Our observations described in the earlier sections indicate clearly the existence of both ergodic (low temperature and high transverse field) region as well as nonergodic (high temperature and low transverse field) regions, separated by a line originating from $T=0,\Gamma=0$ and the passing through the quantum-classical crossover point obtained earlier sudip-sudip ; sudip-yao on the phase boundary of the model. In order to check the dynamical features of these two regions, we have studied the annealing behavior of the system, again using the Suzuki-Trotter effective Hamiltonian with time ($t$) dependent $T$ and $\Gamma$: $T(t)=T_{0}(1-\frac{t}{\tau})$ and ${\Gamma}(t)={\Gamma}_{0}(1-\frac{t}{\tau})$. Here $T_{0}$ and ${\Gamma}_{0}$ correspond to points in the para phase such that they are practically equidistant from the phase boundary line in different parts of the phase diagram. We look for the variation of the annealing time $\tau$ to reach a very low free-energy (corresponding to very small values of $T\simeq 10^{-3}\simeq\Gamma$ to avoid singularities in the effective interaction $H_{eff}$ and dynamics; putting by hand these values of $T$ and $\Gamma$ in the last $t=\tau$ step of the annealing schedule), starting from para phase (high $T$ and $\Gamma$ values). We study annealing of the system when the annealing path (schedule) passes either through ergodic or through nonergodic region (see Fig. 5a). We find that the annealing time $\tau$ remains fairly constant for any annealing path (schedule) passing entirely through the ergodic SG(E) region and becomes strongly system size ($N$) dependent as the path passes through the nonergodic region (see Fig. 5b). It may be mentioned, deep inside the classical region of the spin glass phase (for $S$ value $\gtrsim 1$; see Fig. 5) the annealing time becomes strongly configuration dependent and hence the $N$-dependence of the average value of $\tau$ becomes somewhat irregular.

VI Summary and Discussions

We have studied the order parameter distribution in the spin glass phase of the quantum SK model, both at finite temperature (using Monte Carlo simulation of the effective Hamiltonian (2)) and at zero temperature (using exact diagonalization). For Monte Carlo simulation we have taken system sizes $N=60,120,180,240$ along with the Trotter size $M=15$ (Figs. 2,3). It may be mentioned that we cheeked that the Monte Carlo results remain practically unchanged if for such system sizes we vary the number of Trotter slices $M$ with the system size $N$ keeping the value of the scaled variable $M/N^{z/d}$ constant, where $z$ denotes the dynamical exponent and $d$ is the effective dimension of the system (see sudip-sudip for details). For zero temperature analysis we considered the system sizes $N=10,12,16$ and $20$ (Fig. 4). In the ergodic region SG(E) (see Fig. 5a), the (extrapolated) order parameter distribution is found to converge to a Gaussian form around a most probable value with increase of the system size (see Figs. 2$(a,b)$ for $T\neq 0$). Although at $T=0$ the system sizes we considered are very small, it can be anticipated that we will get a single and narrow peak in the order parameter distribution (see Fig. 4) around a most probable value for thermodynamically large system indicating the ergodicity of the spin glass phase at zero temperature. On the other hand, in the nonergodic region SG(NE) (see Fig. 5a), we get Parisi-type order parameter distribution where the long tail extends upto zero value of order parameter (see Figs. 3$(a,b)$). This (nonzero) weight of the distribution near the origin remains non-vanishing with increase in the system size $N$. This behavior of the order parameter distribution indicates the absence of ergodicity in the system in the SG(NE) region.

These results indicate the different regions of the spin glass phase of the quantum SK model as shown in Fig. 5a. It may be noted that the line separating the low temperature (quantum fluctuation driven) ergodic region of the quantum spin glass phase from the high temperature nonergodic region passes through the quantum-classical crossover point on the spin glass phase boundary obtained earlier sudip-sudip ; sudip-yao . Apart from this low temperature part of the spin glass phase the entire para phase of course remains ergodic.

In order to test the role of this quantum fluctuation induced ergodicity in the spin glass phase here, we have also studied the variation of the annealing time $\tau$ in the finite temperature Suzuki-Trotter Hamiltonian dynamics for $T(t)=T_{0}(1-\frac{t}{\tau})$ and ${\Gamma}(t)={\Gamma}_{0}(1-\frac{t}{\tau})$ to reach a desired low value of the free-energy (corresponding to a very low, but finite, values of $T$ and $\Gamma$ to avoid singularities in the effective interaction flipping dynamics). Here $T_{0}$ and ${\Gamma}_{0}$ values of course correspond to the para phase. For such annealing through the ergodic region we have found $\tau$ to be fairly independent of the system size $N$. However, it clearly starts growing with $N$ as one enters the nonergodic region (see Fig. 5b).

We believe, the numerical results reported here for the quantum SK model establishes the nature of the earlier conjectured sudip-ray ergodicity in the model and its role in quantum annealing sudip-bikas ; sudip-ttc-book ; sudip_atanu_14 of the SK model. It is also possible that the crossover region shrinks as $N\to\infty$. Indeed there are several publications sudip_Goldschmidt ; sudip_sachdev ; sudip_young_17 which contradict our conjecture and suggest these results to be due to the finite size effects in the numerical simulations (of course, the paper by Read et al. suggests ergodicity or absence of replica symmetry breaking as $T\to 0$). The same criticisms are also applied in sudip_young_17 to the experimental and numerical observations sudip-yao for “scrambling” or ergodicity in this system at low enough temperatures. Even if these effects are due to finite system size and the ergodic region becomes narrower with increasing system size, it is important to study such “finite size scaling“ like behavior of the annealing dynamics, so that one can perhaps extrapolate properly the finite size annealing results as the system size approaches the macroscopic limit. Such an extrapolation scheme, if formulated properly, will be extremely useful for the quantum annealing machines (like D-wave sudip_Denchev ) already developed and for developing the quantum machine learning algorithms sudip_Zecchina .

Acknowledgements.

We are grateful to Arnab Chatterjee, Arnab Das, Sabyasachi Nag, Purusattam Ray and Parongama Sen for their comments and suggestions. Bikas Chakrabarti gratefully acknowledges his J. C. Bose Fellowship (DST) Grant. Atanu Rajak acknowledges financial support from the Israeli Science Foundation Grant No. 1542/14.