Haake-Lewenstein-Wilkens approach to spin-glasses revisited

Spin glass problem. Spin glasses (SG) have entered solid state and statistical mechanics in the 1970s, and from the very beginning were considered to be one of the most outstanding and challenging problems of classical statistical physics and theory of disordered and complex systems Mezard ; Chuck , not to mention their quantum version (cf. Sachdev ; Veronica ; Toby and references therein). The most important and elaborated models of spin glasses are: the Edwards-Anderson (EA) model with short range interactions EA , and the Sherrington-Kirkpatrick (SK) model with infinite range interactions SK .

Sherrington-Kirkpatrick model. The SK model was solved approximately by its inventors using replica trick and replica symmetric solution of the equations that ”minimise” the free energy. This solution was clearly physically incorrect, leading to negative entropy at low temperature. G. Parisi found an ingenious way to break the replica symmetry in a hierarchical way Parisi . Parisi’s solution of the SK model turned out to be exact, first as a local extremum of the free energy Dedominicis , and proven rigorously to be unique Talagrand . To deepen the understanding of this amazing results it is also recommended to consult the Ref. Panchenko . Parisi’s solution and Parisi’s order parameter, interpreted in terms of probability of overlaps between different frozen configurations of the spin glass is nowadays accepted commonly. For this achievement, and many others, G. Parisi was awarded the Nobel Prize in physics in 2021 for ”the discovery of the interplay of disorder and fluctuations in physical systems from atomic to planetary scales.”

Edwards-Anderson model. In the case of EA model, we are very far from a rigorous solution. Most of our knowledge is based on numerical simulations on special purpose classical computers, going back to the 1980s BY ; Ogielski . It is widely believed that for Ising EA model there is no SG transition at non-zero temperature in 2D, but there is in 3D and higher dimensions. It is not clear that Parisi’s picture applies in these low dimensions; an alternative is provided by the ”droplet model” of Ref. Huse , which predicts that there exist only one (up to the spin flip) ground state in EA model in 2D, 3D and 4D, but the domain walls, separating the flipped region, from not flipped one, are complex and might even have fractal dimension. Of course, there are some rigorous results concerning EA model (cf. Arguin ), but they are rather weak and very scarce. Thus, the question of the nature of the SG, as well as many other questions concerning the EA model, is open (cf. Chuck and references therein). In recent years various aspects have been studied: ground states in $J=\pm 1$ model Laundry , information theory approach to 3D EA models people1 , absence of Almeida-Thouless line in 3D SG Katz1 , or universality in such systems Katz2 , and several other. The goal of this paper is to look at EA model more than 25 years after the publication of HLW , revising the approach developed there.

HLW approach In 1984 Universität Essen GHS initiated the extremely successful Sonderforschung Bereich ”Unordnung und große Fluktuazionen” with several neighbouring centres. Fritz was a speaker of this initiative for the next 12 years. He convinced Maciej Lewenstein (his summer-time postdoc) and his new PhD student Martin Wilkens to study short-range spin glasses. They formulated a new approach to this problem, based on idea of studying disorder-averaged probability distribution for configurations of two replicas/copies of the system HLW . The idea of the HLW approach is as follows. We consider two replicas/copies of the system and evaluate the joint probability distribution of configurations averaged over the disorder:

where $H(\sigma,\{K_{ij}\}_{\langle ij\rangle})=-\sum_{<ij>}K_{ij}\sigma_{i}\sigma_{j}$, $<ij>$ denotes nearest neighbors, $\langle\langle\cdot\rangle\rangle$ denotes average over disorder, and $Z(K)$ is the partition function calculated for a given configuration of the quenched disorder variables $K_{ij}$. We assume that $K_{ij}$ are iidrv’s (independent identically distributed random variables) with a Gaussian distribution, $P(K)=\exp(-K^{2}/2\Delta^{2})/\sqrt{2\pi\Delta^{2}}$ or a binary distribution, $P(K=\pm\Delta)=1/2$. Note that both distributions are even, that is, invariant under the change of sign of $K_{ij}$. The idea is to absorb the sign of $\sigma_{i}\sigma_{j}$ into $K_{ij}\to K_{ij}\sigma_{i}\sigma_{j}$, and introduce the spin overlap variables $\tau_{i}=\sigma_{i}\sigma_{i}^{\prime}$. We obtain the effective probability distribution for $\tau$’s

Here the number of relevant variables is reduced as we summed over dummy variables. Note that magnetic order for $\tau$’s implies the non-zero EA order parameter $q_{\rm EA}$ and vice versa,

We term $\langle\cdot\rangle$ or $\langle\cdot\rangle_{T}$ the thermal average over possible configurations. Denoting $\alpha=\beta\Delta$, with $\Delta$ the parameter characterizing the probability distributions for the disorder, HLW used a convenient high temperature expansion to calculate (1) up to 12 order in the expansion parameters $\alpha^{2}/(1+\alpha^{2})$ for the Gaussian, and $\tanh^{2}(\alpha)/(1+\tanh^{2}(\alpha))$ for binary case. In effect, they calculated

where effective Hamiltonian contained nearest neighbors couplings $-K_{1}$, next nearest neighbors couplings $K_{2}$, and elementary plaquette terms, $K_{3}$. The coefficients of these terms were explicit functions of temperature ($\alpha$ in the notation of the present paper). The critical surface separating ferromagnetic from paramagnetic region was estimated then using (optimized) real space renormalization group approach. It turned out that in 2D the $H_{\rm eff}$ never enters the ferromagnetic region, in 4D it enters the ferromagnetic region for sure, and in 3D the situation was not clear, suggesting that $H_{\rm eff}$ touches the critical region in a quadratic manner. That would imply that the critical exponents of the spin glass model are two times bigger than those of the standard Ising model, in agreement with the best numerical simulation available at that time.

HLW followers The paper by HLW did not found too many followers, but some very prominent are worth mentioning. Indeed, R. Swendsen with collaborators published two papers on HLW method in Phys. Rev. B in the end of 1980s. In the first one by J.-S. Wang and R.H. Swendsen Swendsen1 , the authors studied Monte Carlo renormalization-group of Ising spin glasses. Application of this approach to the $\pm J$ Ising spin glass showed clear differences between 2D, 3D, and 4D models. The data were consistent with a zero-temperature transition in two dimensions, and non-zero temperature transitions in three and four dimensions. In another paper Swendsen2 Monte Carlo and high-temperature-expansion calculations of a spin-glass effective Hamiltonian were performed. The authors studied the quenched random-coupling spin-glass problem from the point of view of a nonrandom effective Hamiltonian, by Monte Carlo and high-temperature-expansion methods. It was found that the high-temperature series of the spin-glass effective Hamiltonian diverges below the ferromagnetic transition temperature. The Monte Carlo approach does give reliable results at low temperatures. The results were compared with the HLW picture of spin-glass phase transitions.

Present work. In this paper we revise HLW approach. The idea is to estimate $P(\tau)$, performing saddle point/steepest descent approximation in calculating the Gaussian average of the disorder, which should be correct in the limit $\alpha\to\infty$. We argue that the resulting spin model has couplings that are positive in the region where $\tau_{i}\tau_{j}=1$’s, so it has tendency to order ferromagnetically on islands/domains, separated from other domains by negative couplings. In effect, boundary/edge effects start to play a role in estimates of various quantities that may characterize the order in our system.

Laplace’s method. The idea is to calculate the asymptotic behavior of the disorder average using the Laplace method, also known as the saddle point/steepest descent (SPSD) method, which we expect to be asymptotically accurate for $\alpha\to\infty$. The SPSD equations equating to zero the first derivatives of the logarithm of the integrand with respect to $\kappa_{ij}$’s read:

where $\langle\sigma_{i}\sigma_{j}\rangle$ is the thermal average of the neighboring spins correlator, calculated according to the canonical distribution $P(\sigma)=\exp{[\alpha^{2}\sum_{\langle ij\rangle}\kappa_{ij}\sigma_{i}\sigma_{j}]}/Z(\alpha^{2}k)$. There are two possibilities:

• •

  $\tau_{i}\tau_{j}=1$. In this case:

  $$\kappa_{ij}=2(1-\langle\sigma_{i}\sigma_{j}\rangle)>0,$$

  (5)

  so that the corresponding coupling is clearly ferromagnetic.

• •

  $\tau_{i}\tau_{j}=-1$. In this case

  $$\kappa_{ij}=-2\langle\sigma_{i}\sigma_{j}\rangle,$$

  (6)

  and the situation is more delicate. For $\alpha$ large, if $\langle\sigma_{i}\sigma_{j}\rangle)>0$, we expect the coupling $\kappa_{ij}$ to be ferromagnetic, but the above equation implies the opposite. Likewise, if the correlation function is negative, the $\kappa_{ij}<0$ should be ferromagnetic. The contradiction could be avoided if $\kappa_{ij}=0$, but the true situation is more complex, as we will see below, by solving systematically mean field equations. This contradiction is really an expression of frustration in our system!

Hessian matrix. In the zeroth order one can calculate now $P(\tau)$, substituting for $\kappa_{ij}$’s their SPSD values. One can go one step further calculating the Gaussian correction to the SPSD. To this aim we calculate the Hessian matrix of the second derivatives of the logarithm of the integrand. Let us introduce the shortened notation $(ij)=\mu$, $(i^{\prime}j^{\prime})=\nu$, $\sigma_{i}\sigma_{j}=c_{\mu}$, $\sigma_{i^{\prime}}\sigma_{j^{\prime}}=c_{\nu}$, etc. The Hessian matrix reads

Note that the correlations matrix

i.e. it is explicitly positively semi-definite. In effect the Hessian matrix:

so that the logarithm of the integrated function, which we consider is a strictly convex function of many variables, is expected to have one maximum, corresponding to our SPSD solutions. Note also that eigenvalues of the Hessian matrix are all negative and will typically be of order $\alpha^{4}$, and they are bounded in modulus from below by $\alpha^{2}$. One should thus expect that SPSD method should become for $\alpha\to\infty$ asymptotically very precise, if not exact.

Peierls approach. Peierls considers domain walls in a square lattice in 2D, defining them in an unambiguous way. In a ferromagnetic Ising model with the uniform coupling $K$ (with $\beta$ absorbed into $K$), and with periodic boundary conditions on a square of side $L$, and number of sites $N=L^{2}$, with all spins $\tau_{i}=1$ on the boundary, all domain walls are closed. Let $b$ denote the length of the domain’s boundary; We classify them according to length $b$, and within a class of given length we give each a number $i$. A wall of the length $b$ fits into a square of the side $b/4$ and area $b^{2}/16$. Let $m(b)$ be the number of domain walls of length $b$; it is obviously bounded by $m(b)\leq 4N3^{b-1}$. The next step is to consider the quantity $X(b,i)=1$, if the domain wall $(b,i)$ occurs in that configuration, and $X(b,i)=0$ otherwise. Clearly, the number of spins down fulfills:

Peierls estimates then the thermal average of $X(b,i)$ in the Gibbs-Boltzmann ensemble, bounding the partition function from below by the contribution of the configuration, in which all spins inside the considered domain were flipped, and obtaining the bound $\langle X(b,i)\rangle\leq\exp[-2\beta Kb]$, which leads to the desired result. Notably, it can be generalized to higher dimensions, with a little extra effort to estimate the entropy of contours, see Bonati for an elementary discussion and references therein for original work.

Thouless argument. In order to investigate the sensitivity to boundary conditions, we calculate $P_{+}=P(\{\tau_{i}=1\}_{\rm for\ all\ i^{\prime}s})$ on a cylinder in $d$ dimensions of cross-section $L^{d-1}$ and length $L$, and compare it to $P_{-}=P(\{\tau_{i}=1\}_{\rm on\ the\ left},\{\tau_{i}=-1\}_{\rm on\ the\ right})$ with $\tau_{i}=\pm 1$ on the left (right ) of a domain wall (DW), correspondingly. We determine the parameter $\delta=\ln(P_{+}/P_{-})$; Ferromagnetic order for $\tau$’s (SG order for $\sigma$’s) is indicated by $\delta>0$

Local mean field theory. We assume translation symmetry in $d-1$ transverse dimensions, so that magnetisation depend only on one index, $i$, enumerating the layer, and the couplings depend on two indices, enumerating involved single layer (two neighboring layers). Using standard mean field theory (MF), we find the magnetization that is the thermal average of $\langle\sigma_{i}\rangle$ at the $i$-th layer as

with

and boundary conditions:

Quantities to be determined. Our aim is to calculate logarithm of the probability $P_{+}$, $P_{-}$ and $\delta=\ln(P_{+}/P_{-})$. We denote $\ln(P_{\pm})=H_{\pm}$, and call it ”energy” in the following, so that

Being an extensive quantity, the energy of the system divided by the volume in all but one dimension is

with $F_{i}\left(m\right)=\alpha^{2}(2(d-1)\kappa_{i,i}m_{i}+\kappa_{i-1,i}m_{i-1}+\kappa_{i,i+1}m_{i+1})$. Note that we have included in this expression the term coming from the Gaussian fluctuations around the SPSD solution. The above quantity in the leading order should be a linear function of the cylinder’s length,

Gaussian fluctuation terms Generally speaking, Gaussian fluctuation terms play a sub-leading role, as expected. We approximate $\ln({\rm det}(\hat{H}_{\pm}))=\sum_{\mu}\ln(\lambda_{\mu})\approx\sum_{\mu}\ln(\hat{H}_{\mu\mu})$, that is the sum of logarithms of eigenvalues by the sum of logarithms of diagonal elements of the Hessian matrix. Noting that

we obtain

where the SPSD solutions for $m$’s are calculated for the case $\pm$ accordingly. The above expression undergoes, obviously, further simplifications under the translation symmetry.

High $\alpha$ regime Before going to numerical solutions, we first analyze the asymptotic regime $\alpha\to\infty$, where $A(\alpha)$ can be estimated analytically. We consider MF equations in the bulk of the $d$-dimensional hyper-cubic lattice. The corresponding self-consistent equations in the bulk are:

We transform the first two into an equation of $x=2\alpha^{2}d\kappa$.

For large $\alpha$ we get $x=\ln(4\alpha^{2}d)/2$, and $\kappa=\ln(4\alpha^{2}d)/(4\alpha^{2}d)$. As expected, $\kappa\to 0$ as $\alpha\to\infty$, and $m\to 1$, but $2\alpha^{2}d\kappa$ diverges as $\ln(4\alpha^{2}d)$. Elementary analysis leads to the result:

i.e. as expected $\ln(P_{\pm})=H_{\pm}$ becomes negative at large $L$ (when our analysis makes sense) and at large $\alpha$ (when SPSD should work well); $A(\alpha)$ diverges with $\alpha$, but very slowly, only logarithmically.

”Phase transition” at moderate $\alpha$ The solution of the MF equations change character as $\alpha$ grows from small values (when all $m_{i}=0$) to larger values (when all $m_{i}\neq 0$). We infer the existence of this ”phase transition” at a finite $\alpha$ by imposing that solutions get trivial at that point, $\alpha_{T}$. This way we can approximate Eq. (IV) for temperatures close to $\alpha_{T}$ as a series expansion for small $m_{i}$ to get:

To first order, $\kappa_{i,j}=2$ and $m_{i}=m$ $\forall i$, so we find the critical temperature:

Numerical calculations By numerically solving the system of equations $F(m)=m_{i}$ for $1\leq i\leq L$ and taking into account that in positions $i=\{0,L+1\}$ there are no spins and therefore conditions Eq. (12) and Eq. (13) apply, we find non-trivial solutions above a certain temperature threshold, see Fig. 1.