text:
# Low Energy Excitations in Spin Glasses from Exact Ground States

## I Introduction

There is still considerable debate about the nature of the spin glass state in finite dimensional spin glasses. Two principal theories have been investigated: the “droplet theory” proposed by Fisher and Husefh  (see also Refs. bm, ; mcmillan, ), and the replica symmetry breaking picture of Parisiparisi ; mpv ; bindery . In the droplet theory the lowest energy excitation of length scale $`ł`$ (a “droplet”) has energy of order $`ł^\theta `$ where $`\theta `$ is a positive exponent. Furthermore, the droplets have a surface with fractal dimension, $`d_s`$, less than the space dimension $`d`$.

Replica symmetry breaking (RSB) is well established in mean field theory, but it remains to be proven in finite dimensions. The precise nature of RSB in finite dimensions is not uniquely defined but it is generally agreed that a key feature of RSB is the existence of excitations whose energy, unlike that of droplets, remains of order unity even as their size tends to infinity. Furthermore, upon the creation of such a large scale, finite energy excitation, a finite fraction of the bonds change state (from satisfied to unsatisfied, or vice-versa) or, equivalently, the surface of these excitations is space filling, i.e. $`d_s=d`$.

Recently, Krzakala and Martinkm  (KM), and two of uspy-bulk  (PY), have argued, on the basis of numerical calculations at zero temperature, in favor of an intermediate scenario where there are large scale excitations whose energy does not increase with size, as in RSB, but which have a surface with $`d_s<d`$. Following KM we shall denote this the “TNT” scenario. In the TNT scenario it is necessary to introduce two exponents which describe the growth of the energy of an excitation of scale $`L`$: (i) $`\theta (>0)`$ such that $`L^\theta `$ is the typical change in energy when the boundary conditions are changed, for example from periodic to anti-periodic, in a system of size $`L`$, and (ii) $`\theta ^{}`$, which characterizes the energy of clusters excited within the system for a fixed set of boundary conditions ($`\theta ^{}`$ was called $`\theta _g`$ in Ref.km, ). The TNT picture has been challenged (although in opposite senses) by Marinari and Parisimapabulk  and by Middletonmiddleton . Subsequently, low temperature Monte Carlo simulationskpy  have found results consistent with the TNT scenario. The RSB, droplet, TNT and some other scenarios have been also studied by Newman and Steinns-old ; ns . For some recent related work, see Refs. mezard, ; moore, .

The work of KM and PY determined the ground state with and without a certain perturbation (which was different in the two cases), designed so that the ground state of the perturbed system is a large scale excitation of the original system. They used heuristic algorithms, i.e. algorithms which are not guaranteed to give the exact ground state, although both KM and PY argue that they do find the exact ground state in most cases.

In this paper, we reconsider the problem of determining $`\theta ^{}`$ and $`d_s`$, following the perturbation approach of PY, described in Section II, but we apply an exact algorithm, known as “branch and cut” junger-original , so we are guaranteed that the true ground state is reached every time. Exact optimization algorithms have been used before for spin glasses, see e.g. Refs. kobe, ; dombgreen, ; desim, , but, to our knowledge, their use in three-dimensions has been restricted to smaller sizes than studied here, and they were not used to investigate the real-space structure of the low-energy excitations.

Our implementation of the branch and cut technique can handle significantly larger sizes for free boundary conditions (bc) than for periodic bcliers , so we use free bc here. We consider a different (and enlarged) set of observables than PY, in the attempt to gain a fuller understanding of what picture fits better the whole set of observables. We also perform a similar analysis of the data of PY, who used periodic bc, in order to investigate the effects of different types of boundary conditions. The various pictures discussed refer to the large volume limit, while the sizes that can be currently reached are rather small. We will therefore pay particular attention to properly take into account corrections to scaling. In particular, we will try to determine what values of the parameters $`\theta ^{}`$ and $`dd_s`$ fit the data in the more “natural” way, namely with the smallest corrections to scaling for the range of sizes considered.

A summary of our results is as follows. We find that for periodic bc, a simple scaling ansatz fits the results in a natural way, i.e. with negligible corrections to scaling and no adjustable parameters besides $`dd_s`$ and $`\theta ^{}`$. This gives $`dd_s=0.42\pm 0.03`$; $`\theta ^{}=0.01\pm 0.03`$ (the meaning of the error bars will be explained later), which agrees with the results of PY, and is compatible with the TNT picture. We cannot rule out crossover to either the droplet or the RSB picture at length scales larger than our system sizes, but these scenarios, especially the latter, would require larger corrections to scaling than the TNT picture.

For free bc, all forms of fitting require some corrections to scaling. The most natural scenario, in the sense explained above, gives $`dd_s=0.45\pm 0.02`$; $`\theta ^{}=0.18\pm 0.03`$, with small corrections (of the order of 3$`\%`$), which is compatible with the droplet picture. Allowing somewhat larger corrections (of order 10$`\%`$), the data are also compatible with $`\theta ^{}=0`$, namely with the TNT picture. Finally, if we allow for very large corrections, the data are also consistent with the RSB picture.

In the second part of the paper, we analyze the performance of the branch and cut algorithm. We find that the number of elementary operations required to find the ground state increases exponentially with the size, as expected since computing a ground state of a three-dimensional spin glass system is an $`𝒩𝒫`$-hard problembarahona . We also find, interestingly enough, that the CPU time is larger for samples in which there is an excited state close in energy to the ground state energy, yet different from the ground state in the orientation of a large number of spins. We are not aware of any previous quantitative measures of this trend, which we expect to be common to other algorithms as well.

The rest of this paper is organized as follows. In Section II we describe the method of perturbing the ground states to get information about low energy excitations, introduced by PY. Our results for the nature of the large scale, low energy excitations are given in Section III. A short description of the branch and cut algorithm used is given in Section IV, and the performance of the algorithm is analyzed in Section V. Our conclusions are summarized in Section VI.

## II Ground state perturbation method

The Hamiltonian of the spin glass model is given by

$$=\underset{i,j}{}J_{ij}S_iS_j,$$
(1)
where the sites $`i`$ lie on a simple cubic lattice with $`N=L^3`$ spins in dimension $`d=3`$, $`S_i=\pm 1`$, and the $`J_{ij}`$ are nearest-neighbor interactions chosen from a Gaussian distribution with zero mean and standard deviation unity. Free boundary conditions are applied in all directions.

For a given set of bonds we determine the exact ground state using a branch and cut algorithm discussed in Section IV. Let $`S_i^{(0)}`$ be the ground state spin configuration. As in PY we then perturb the couplings $`J_{ij}`$ by an amount proportional to $`S_i^{(0)}S_j^{(0)}`$ in order to increase the energy of the ground state relative to the other states and therefore possibly induce a change in the ground state. This perturbation, which depends upon a positive parameter $`ϵ`$, is defined by

$$\mathrm{\Delta }_ϵ=\frac{ϵ}{N_b}\underset{i,j}{}S_i^{(0)}S_j^{(0)}S_iS_j,$$
(2)
where $`N_b=dL^{d1}(L1)`$ is the number of bonds in the Hamiltonian. We denote the unperturbed ground state energy by $`E^{(0)}`$ and the perturbed energy of the same state by $`E_ϵ^{(0)}`$. The energy of the unperturbed ground state will thus increase exactly by an amount $`\mathrm{\Delta }E^{(0)}E_ϵ^{(0)}E^{(0)}=ϵ.`$ The energy of any other state, $`\alpha `$ say, will increase by the lesser amount $`\mathrm{\Delta }E^{(\alpha )}E_ϵ^{(\alpha )}E^{(\alpha )}=ϵq_l^{(0,\alpha )},`$ where $`q_l^{(0,\alpha )}`$ is the “link overlap” between the states “0” and $`\alpha `$, defined by

$$q_l^{(0,\alpha )}=\frac{1}{N_b}\underset{i,j}{}S_i^{(0)}S_j^{(0)}S_i^{(\alpha )}S_j^{(\alpha )},$$
(3)
in which the sum is over all the $`N_b`$ nearest neighbor bonds. Note that the total energy of the states changes by an amount of order unity.

As we apply the perturbation, the energy difference between a low energy excited state and the ground state decreases by the amount

$$\mathrm{\Delta }E^{(0)}\mathrm{\Delta }E^{(\alpha )}=ϵ(1q_l^{(0,\alpha )}).$$
(4)
If there is at least one excited state such that $`E^{(\alpha )}E^{(0)}<\mathrm{\Delta }E^{(0)}\mathrm{\Delta }E^{(\alpha )}`$, then one of these excited states will become the ground state of the perturbed Hamiltonian. We denote the new ground state spin configuration by $`\stackrel{~}{S}_i^{(0)}`$, and indicate by $`q_l`$ and $`q`$, with no indices, the link- and spin-overlap between the new and old ground states $`S_i^{(0)}`$ and $`\stackrel{~}{S}_i^{(0)}`$, where $`q`$ is defined as usual by $`q=1/NS_i^{(0)}\stackrel{~}{S}_i^{(0)}`$.

Due to the spin flip symmetry of the Hamiltonian (1), the ground state is doubly degenerate, and therefore the distribution of $`q`$ is symmetric reality  around $`q=0`$. Hence, in the rest of the paper we will restrict ourselves to $`q0`$ without loss of information.

Consider the probability $`P(ϵ,L)`$ (with respect to the random couplings) that $`q`$ is less than unity, i.e. that $`S_i^{(0)}`$ and $`\stackrel{~}{S}_i^{(0)}`$ differ in a finite fraction of the spins. As discussed by PY, we assume that $`P(ϵ,L)`$ is dominated by those samples in which $`S_i^{(0)}`$ and $`\stackrel{~}{S}_i^{(0)}`$ differ by flipping a single connected cluster of spins, with linear size $`L`$. Deviations from this assumption give rise to corrections to scaling, as pointed out by Middletonmiddleton , and will be analyzed in Section III. There are two energy scales in the problem: the typical energy above the ground state of such an excitation, which scales as $`L^\theta ^{}`$ ($`\theta ^{}=\theta `$ in the droplet picture), and the threshold energy of Eq. (4), which scales as $`ϵL^{(dd_s)}`$ since $`1q_l`$ is proportional to the surface of the excitation, $`1q_lL^{(dd_s)}`$. Hence, the dimensionless probability $`P(ϵ,L)`$ is a function of the ratio of these two energy scales:

$$P(ϵ,L)=g(ϵL^\mu ),$$
(5)
where $`g(x)`$ is a scaling function and

$$\mu \theta ^{}+dd_s.$$
(6)
From this we obtain scaling relations for various observables. For example, since $`1qO(1)`$ and $`1q_lL^{(dd_s)}`$, we obtainpy-bulk :

$`1q`$ $`=`$ $`F_q(ϵL^\mu )`$ (7)
$`1q_l`$ $`=`$ $`L^{(dd_s)}F_{q_l}(ϵL^\mu ),`$ (8)
where $`\mathrm{}`$ is the average with respect to the random couplings. By measuring $`1q`$ and $`1q_l`$ we can then determine $`dd_s`$ and $`\theta ^{}`$, the two exponents discriminating the various pictures of the spin glass phase discussed in Section I.

For small $`ϵ`$, we expect the probability that the ground state changes to be proportional to $`ϵ`$ (for fixed $`L`$), which implies $`g(x)x`$ for $`x0`$. Hence $`F_q(x)`$ and $`F_{q_l}(x)`$ also vary linearly for small $`x`$, and the asymptotic scaling behavior for $`Lϵ^{1/\mu }`$ is

$`1q`$ $``$ $`ϵL^\mu ,`$ (9)
$`1q_l`$ $``$ $`ϵL^{\mu _l},`$ (10)
where

$$\mu _l\theta ^{}+2(dd_s).$$
(11)
In the RSB case, $`dd_s=\theta ^{}=0`$, and therefore $`\mu =\mu _l=0`$. The scaling relations in Eqs. (7, 8) reduce in this case to

$$1q=F_q(ϵ),1q_l=F_{q_l}(ϵ)\text{(RSB)},$$
(12)
and the asymptotic scaling behavior to

$$1qϵ,1q_lϵ\text{(RSB)}.$$
(13)
We see that both scaling and asymptotic scaling are in a sense trivial in RSB since the $`L`$ dependence disappears. Nevertheless, we will still use the term scaling.

It is also convenient to analyze just those samples in which the unperturbed and perturbed ground states are very different, i.e. where $`qq_{\mathrm{max}}`$, a threshold value. Denoting such restricted averages by $`\mathrm{}_c`$, we have

$$1q_l_c=L^{(dd_s)}F_{q_l}^c(ϵL^\mu ).$$
(14)
This is of the same form as in Eq. (8), but, for sufficiently small $`q_{\mathrm{max}}`$, the behavior of the scaling functions $`F_{q_l}(x)`$ and $`F_{q_l}^c(x)`$ at small argument will be different for the following reason. If we average over all samples we need to include the probability $`P(ϵ,L)`$ that the perturbation generates an excitation with $`q<1`$. This is proportional to $`ϵL^\mu `$ for $`ϵL^\mu 1`$, which is the reason why $`F_q(x)x`$ for small $`x`$. However, this factor is automatically taken into account in the selection of the samples in the restricted average in Eq. (14), and so should not be included again when performing the average. As a result, $`F_{q_l}^c(x)`$ tends to a constant for $`x0`$, therefore the asymptotic scaling is

$$1q_l_cL^{(dd_s)}.$$
(15)
In particular, in RSB this becomes

$$1q_l_c\text{const.}\text{(RSB)}.$$
(16)
Note that in both cases the asymptotic scaling is independent of $`ϵ`$.

When analyzing the numerical data, we must be aware that there are corrections to both (simple) scaling and asymptotic scaling that occur when $`L`$ is not large enough. Corrections to simple scaling take the form of additive corrections to relations such as Eqs. (5), (7), (8), and (14), whose amplitude is characterized by a correction-to-scaling exponent $`\omega `$. For example, including the leading correction, Eq. (14) becomes

$$1q_l_c=\frac{1}{L^{dd_s}}\left\{F_{q_l}^c(ϵL^\mu )+\frac{1}{L^\omega }G_{q_l}(ϵL^\mu )\right\}.$$
(17)
For $`ϵL^\mu 0`$, this gives the correction to asymptotic scaling corresponding to Eq. (15)

$$1q_l_c=\frac{1}{L^{dd_s}}\left(a+\frac{b}{L^\omega }\right).$$
(18)
For the RSB case, this goes over to

$$1q_l_c=a+\frac{b}{L^\omega },$$
(19)
rather than Eq. (16).

Even when these corrections to (simple) scaling are negligible and the scaling form, such as Eq. (14), is valid, the argument of the scaling function may not be sufficiently small for asymptotic scaling to hold. In this regime, when fitting the data to asymptotic scaling we have to consider further corrections to (asymptotic) scaling, whose form is obtained by expanding the scaling function in its argument. For example, the leading correction to Eq. (15), coming from expanding the $`F_{q_l}^c`$ in Eq. (14) to second order, will be

$$1q_l_c=\frac{1}{L^{dd_s}}\left(a+b\frac{ϵ}{L^\mu }\right)$$
(20)
which goes over to $`1q_l_c=a+bϵ`$ in RSB. In general, both types of corrections need to be borne in mind when fitting the data.

## III Results

We applied the perturbation method described in the previous Section to systems of size $`L=4,6,8,10,`$ and $`12`$. For each size, we considered five values of the perturbation strength $`ϵ`$, namely $`ϵ/\tau =\frac{1}{4},\frac{1}{2},1,2,`$ and 4, where $`\tau =\sqrt{6}`$ is the mean field transition temperature, except for $`L=12`$ for which only $`ϵ/\tau =\frac{1}{4}`$ and $`1`$ were considered. We choose this value of $`\tau `$ so we can compare our results with the results of PY for periodic bc. In order to discriminate between the different pictures, it is important to have high statistics. Table 1 reports the number of samples computed for each size. Note that the number of samples necessary to achieve a given statistical error increases as $`ϵ`$ decreases, since the fraction of samples in which $`\stackrel{~}{S}^{(0)}S^{(0)}`$ decreases.

### III.1 Spin and link overlap

#### III.1.1 Qualitative analysis

We start with a qualitative analysis of the results. In Fig. 1, we show scatter plots in the $`(q,q_l)`$ plane for $`L=4,10`$ and $`ϵ/\tau =\frac{1}{4},1`$, where each point represents one of 2000 randomly generated samples. Clearly, the link- and spin-overlap are strongly correlated. We note that, as $`ϵ`$ decreases, there are less points with small $`q`$, and that as $`L`$ increases the data shift towards larger values of $`q_l`$. Similar plotsphdthesis  for periodic bc show that $`q_l`$ is significantly lower than for free bc. While $`q`$ has a large variance (the points are distributed along the whole $`q`$ axis), the link-overlap has a much smaller variance, which decreases as $`L`$ increases, suggesting that in the thermodynamic limit $`q_l`$ either tends to one or becomes a well defined function of $`q`$.

To quantify this, in Fig. 2 we show the standard deviation of $`q_l`$

$$\sigma =\sqrt{q_l^2_cq_l_c^2},$$
(21)
restricted over samples with $`qq_{max}`$, as a function of the system size for $`ϵ/\tau =1`$. We take $`q_{\mathrm{max}}=0.2`$, since we are interested in the region of small $`q`$, which corresponds to large-scale excitations. A power law $`\sigma =aL^\delta `$ fits well the data with $`\delta =0.52\pm 0.03`$ ($`\chi ^2=1.80`$, the best fit is shown in Fig. 2). Here and in the following, unless stated otherwise, the error bars on the fit parameters are purely statistical in relation to the fitting form considerederrors . Restricting the average in Eq. (21) to different intervals of $`q`$ gives results also compatible with a power law. A vanishing $`\sigma `$ in the thermodynamic limit is consistent with RSB, which predicts that $`q_l`$ is a (nontrivial) function of $`q`$. It is also trivially consistent with the droplet model or the TNT picture, where $`q_l=1`$ for all $`q`$.

We also measured how the standard deviation of $`q`$ varies with $`L`$, finding that it varies between 0.28 and 0.32. It can be fitted both to a constant (as expected in RSB) or to a power law with a small exponent around 0.1. However the error bars are very large hence the fits are not very informative.

Under the RSB hypothesis, it is interesting to study the functional relationship between $`q`$ and $`q_l`$. In Fig. 3 we show the average value of $`q_l`$, restricted to intervals $`q[q_{\mathrm{min}},q_{\mathrm{max}}]`$, as a function of the mean value of $`q`$ in each intervalinterval  for $`ϵ/\tau =1`$. For fixed $`L`$, a quadratic form $`q_l=\alpha (L)+\beta (L)q^2`$, motivated by the infinite-range Sherrington-Kirkpatrick model where $`q_l=q^2`$, fits well the data for $`q`$ less than some cut-off value which increases with $`L`$ (see Fig. 3). The quadratic fit works well also for other values of $`ϵ`$, and $`\alpha (L)`$ and $`\beta (L)`$ show a modest variation with $`ϵ`$. We also tried global fits including data for all values of $`q`$ and $`L`$, obtaining similar results.

Extrapolating $`\alpha (L)`$ and $`\beta (L)`$ to $`L\mathrm{}`$ with fits of the form $`\alpha (L)=\alpha +b/L^c`$, $`\beta (L)=\beta +b^{}/L^c^{}`$, we obtain

$$q_l=(0.77\pm 0.02)+(0.27\pm 0.03)q^2,$$
(22)
where again the errors are purely statistical for the functional form considered. This nontrivial relation between $`q`$ and $`q_l`$ in the thermodynamic limit is consistent with RSB, while in the droplet or TNT pictures the data in Fig.3 would shift to $`q_l1`$ in this limit.

The power law form $`1\alpha (L)=b/L^c`$, $`\beta (L)=b^{}/L^c^{}`$, which implies $`q_l=1`$ in the large volume limit, fits poorly the data if we include all sizes. However, if we exclude the $`L=4`$ data, the quality of the fit becomes as good as that of the RSB fit just discussed. Hence, allowing for small corrections to scaling, the droplet or TNT scenario are consistent with the data.

This already shows that care must be taken to properly consider corrections to scaling when comparing the merits of the fits to various pictures. In the following, we will investigate in detail the validity of the various pictures by considering several observables and explicitly discussing the corrections to scaling for each picture.

#### III.1.2 Determination of $`dd_s`$

We start with the determination of $`dd_s`$ from various observables. We will show that for all observables, a wide range of values of $`dd_s`$ fits well the data when allowing corrections to scaling, but that for all observables the smallest corrections are attained for a value of $`dd_s`$ around 0.44 as in PY.

The main part of Fig. 4 plots $`1q_l_c`$ as a function of $`L`$, for various values of $`ϵ`$ and $`q_{\mathrm{max}}=0.4`$ ($`q_{\mathrm{max}}=0.2`$ gives essentially the same results). First, note that, independently of what picture holds in the $`L\mathrm{}`$ limit, the data deviate significantly from asymptotic scaling, see Eq. (15), in which the various $`ϵ`$ values should collapse on a single curve. Second, the data have a noticeable positive (upward) curvature for all values of $`ϵ`$. In Section III.3 we will show that in the case of periodic bc the data have a much smaller dependence on $`ϵ`$ and a much smaller curvature (see Fig. 12).

In order to determine how the various pictures fit the data of Fig. 4, we start by considering, following Ref.mapabulk, , the following three functional forms:

$`\text{Form I:}1q_l_c`$ $`=`$ $`a+b/L^c`$
$`\text{Form II:}1q_l_c`$ $`=`$ $`a+b/L+c/L^2`$ (23)
$`\text{Form III:}1q_l_c`$ $`=`$ $`b/L^c`$
Form I corresponds to the RSB prediction including the leading correction to scaling, see Eq. (19), with $`c\omega `$. Form II is a different parameterization of the corrections to scaling. Form III corresponds to the asymptotic behavior of both the TNT and droplet pictures without corrections to scaling, see Eq. (15), with $`cdd_s`$.

The results of these fits (performed by $`\chi ^2`$ minimization) are reported in Table 2. From the Table we see that Forms I and II, appropriate to RSB, fit well the data with a low $`\chi ^2`$ and $`a>0`$ outside the error bars. The best fits with Form I are shown by the dashed lines in Fig. 4. The variation of $`a`$ between Forms I and II provides a measure of the systematic error associated with the unknown corrections to scaling. Within this error, $`a`$ is independent of $`ϵ`$, as predicted by RSB. Therefore, the data for $`1q_l_c`$ are compatible with RSB, and our central estimate under the RSB hypothesis is

$$\underset{L\mathrm{}}{lim}1q_l_c=0.20\pm 0.02\text{(RSB)},$$
(24)
where the error takes into account also the uncertainty in the form of the corrections to scaling, assuming that the corrections considered in either Form I or II describe well the data in the whole range of sizes considered. Marinari and Parisi mapabulk  fitted Form I (resp. II) to their data for periodic boundary conditions, $`L14`$, and $`ϵ/\tau =4`$, and obtained $`a=0.24`$ (resp. $`a=0.30`$), from which we estimate a central value $`a=0.27\pm 0.05`$. This is just in agreement with our estimate above for free bc, suggesting that the infinite volume limit of $`1q_l_c`$, if nonzero, may be independent of the boundary conditions, although we do not have an argument why this should be the case.

The power law Form, III, appropriate to the droplet model or the TNT scenario, does not fit well the data if we include all the sizes, but if we exclude $`L=4`$, it does fit well for $`ϵ/\tau <2`$, and the fit parameters $`b`$ and $`c`$ are almost independent of $`ϵ`$. The quality of the fit of Form III is worse than that of Forms I and II, but still acceptable. The main point we want to stress, however, is that the worse fit of Form III alone does not necessarily favor the RSB picture, since Form III does not include corrections to scaling, while Forms I and II do. In other words, Forms I and II are rather “forgiving” with the RSB picture, allowing corrections of magnitude 100$`\%`$ \- 200$`\%`$ of the predicted $`L`$-independent asymptote, while Form III demands that the power-law scenario fits with corrections smaller than the (very small) statistical errors. By looking at Figure 4, it is clear that the data are closer to a power law than to an $`L`$-independent behaviour.

Therefore, in order to try a comparison that puts the various pictures on an equal footing, we performed fits with the following more general functional form:

$$\text{Form IV:}L^{dd_s}1q_l_c=a+b/L^c$$
(25)
where we fix $`dd_s`$ and minimize the $`\chi ^2`$ with respect to $`a,b,c`$, repeating the procedure for different values of $`dd_s`$. For $`dd_s=0`$, Form IV reduces to Form I, while for $`dd_s>0`$, it corresponds to Form III plus a correction-to-scaling term, with correction-to-scaling exponent $`\omega =c`$. We find that, as we might have expected from the previous discussion, Form IV fits well the data for a wide range of values of $`dd_s`$. For example, for $`ϵ/\tau =\frac{1}{4}`$, a value of $`dd_s`$ between 0 and 0.45 gives a goodness-of-fit parameter $`Q0.43`$, which is entirely acceptable.

This shows that, when allowing for corrections to scaling for all pictures, the droplet or TNT pictures are as good as RSB as far as the statistical quality of the fits is concerned. However, within the interval of acceptable values of $`dd_s`$, clearly the larger is $`dd_s`$ the smaller are the corrections to asymptotic scaling. For example, for $`ϵ/\tau =\frac{1}{4}`$ and $`dd_s=0.42`$, the correction term $`b/L^c`$ in Form IV amounts to only $`0.1\%`$ of the total for $`L=12`$, while for $`dd_s=0`$ it amounts to $`43\%`$. Hence a large value of $`dd_s`$ may be regarded as more “natural” in this range of sizes.

If we impose that the correction to scaling for $`L8`$ is less than an (arbitrary) limit of $`3\%`$, we obtain the estimate

$$dd_s=0.44\pm 0.03$$
(26)
where the error is purely statistical within this assumption. In Table II we show the fits obtained with Form IV imposing this value. This agrees with the estimate $`dd_s=0.42\pm 0.02`$ of PY for periodic bc (see also Section III.3 of this paper), which is reassuring since $`dd_s`$ should not depend on the boundary conditions.

Note that for $`dd_s=0`$, corresponding to Form I, corrections within $`3\%`$ from the asymptotic limit would only be attained for a size $`L200`$. We also note that, as we discussed in Section II, even in the regime where corrections to scaling are negligible, asymptotic scaling sets in only for $`Lϵ^{1/\mu }`$. This explains why, if $`dd_s0.44`$, the quality of the power-law fit in Table II becomes progressively worse as $`ϵ`$ increases. In particular, the deviation from asymptotic scaling is very significant for $`ϵ/\tau =4`$, and hence from the data of $`ϵ/\tau =4`$ alone one should not necessarily concludemapabulk  that an asymptotic power-law behavior is ruled out. This is seen also in the inset of Fig. 4, which shows that, if we exclude $`L=4`$, the data are compatible with the scaling relation Eq. (8), where the exponent $`\mu `$ is independently determined below.

PY determined $`dd_s`$ from the ratio $`R=1q_l/1q`$ which has the same scaling behavior as the quantity $`1q_l_c`$ used here, namely $`R=L^{(dd_s)}F_R(ϵ/L^\mu )`$, with $`F_R(x)`$ const. as $`x0`$. Middleton middleton  observed that, in two dimensions, small droplets introduce significant corrections to scaling, and suggested that this may be the case also in three dimensions, possibly invalidating the conclusions of PY. The quantity $`1q_l_c`$ is less sensitive to these corrections since, by restricting to small $`q`$, small droplets should give a smaller contribution, because we have eliminated the part at large $`q`$ where one can have only small droplets. Hence, to investigate these corrections, we fitted our data for $`R`$ with Forms I-IV above (with $`R`$ replacing $`1q_l_c`$). The results we find are very similar to those for $`1q_l_c`$: Forms I and II fit well the data with a low $`\chi ^2`$, giving $`a=0.27\pm 0.03`$ independent of $`ϵ`$ within the error bars. A power law fits well the data if we exclude $`L=4`$, with an exponent $`dd_s=0.43\pm 0.03`$ nearly independent of $`ϵ`$ and in agreement with Eq. (26). The residual dependence on $`ϵ`$ is well accounted for by a scaling plot similar to the inset in Fig. 4. Form IV also fits well the data for a wide range of values of $`dd_s`$. Again, a power law is more natural in the sense that corrections to scaling are smaller, and the smallest corrections are obtained for $`dd_s`$ around 0.43 as for $`1q_l_c`$. We interpret the fact that the two quantities give the same value of $`dd_s`$ as evidence that corrections due to small droplets are indeed not important in three dimensions in this range of sizes. In Section III.3 we will show that this is also the case for periodic bc.

To summarize this part, the data for both $`R`$ and $`1q_l_c`$ are compatible with a wide range of values of $`dd_s`$ between zero (corresponding to RSB) and $`0.44`$, but a value at the higher end of this range describes the data in a more natural way, in the sense that the corrections to scaling are smaller.

#### III.1.3 Determination of $`\theta ^{}`$

Next, we turn to the exponent $`\mu `$ defined in Eq. (6), from which we will extract the exponent $`\theta ^{}`$ which is the other exponent, with $`dd_s`$, characterizing the spin glass phase. To this end we consider the ratio

$$B=\frac{1q_l^2}{(1q_l)^2}$$
(27)
which follows the scaling law

$$B=F_B(ϵ/L^\mu ).$$
(28)
The factor $`L^{dd_s}`$ does not appear here since it cancels between numerator and denumerator of Eq. (27), thus allowing us to determine $`\mu `$ independently of $`dd_s`$. Following the arguments in Section II, we expect $`F_B(x)x`$ for small $`x`$ since both $`L^{dd_s}1q_l`$ and $`L^{2(dd_s)}(1q_l)^2`$ vary as $`ϵ/L^\mu `$ for $`ϵ/L^\mu 0`$, hence the asymptotic scaling of $`B`$ is $`Bϵ/L^\mu `$.

To determine $`\mu `$, we fit the scaling law Eq. (28) to our data assuming a polynomial form of order $`n`$ for $`F_B(x)`$, namely $`F_B(x)=_{i=0,n}c_ix^i`$, with $`c_0`$ set to zero in order to satisfy the asymptotic scaling $`F_B(x)x`$ as $`x0`$. We repeat the fit in an interval of values for $`\mu `$, and determine the value of $`\mu `$ which gives the best fit, varying $`n`$ until the $`\chi ^2`$ of the best fit becomes approximately constant. In this way we obtain

$$\mu =0.63\pm 0.03,$$
(29)
where the error is purely statistical, under the assumption that the corrections to scaling are smaller than the statistical errors of the data. As shown in Fig. 5, scaling is quite satisfactory, with all the data collapsing on one curve, although the data for different $`ϵ`$ overlap only slightly. The best fit for $`n=6`$ is displayed by the continuous line. We emphasize that this scaling plot is obtained with only one adjustable parameter, $`\mu `$. Note that in the asymptotic scaling regime the data should follow a straight line (power-law), while the data in the figure show a pronounced curvature. Significant corrections to asymptotic scaling must be expected for large $`ϵ`$, since $`B`$ must satisfy the inequality $`B1`$. The dashed line in Fig. 5 represents the linear term of $`F_B(x)`$, corresponding to asymptotic scaling, and the deviation from it gives a measure of the corrections to asymptotic scaling. The $`ϵ/\tau =\frac{1}{4}`$ data are quite close to asymptotic scaling, while the data for large $`ϵ`$ deviate significantly from it. Another manifestation of these corrections is that, if we fit the data with a power-law $`B=b/L^{\stackrel{~}{\mu }}`$, the effective exponent $`\stackrel{~}{\mu }`$ varies strongly with $`ϵ`$, converging towards 0.63 for $`ϵ0`$.

In RSB, $`Bϵ`$ as $`L\mathrm{}`$ since $`\mu =0`$. To test the RSB prediction, we performed fits of $`B`$ using Form I and Form II in Eq. (23), (where $`1q_l_c`$ is replaced by $`B`$). Form I gives unphysical (negative) values of $`a`$, while Form II gives an acceptable fit with a positive $`a`$ roughly proportional to $`ϵ`$. Therefore, the data for $`B`$ cannot rule out RSB. Note that, if RSB holds asymptotically, the data in Fig.5 would deviate from the scaling curve for larger $`L`$, saturating to a constant value for small values of $`ϵ/L^\mu `$. The good data collapse we observed, therefore, would be entirely accidental.

We believe that the observed data collapse is a good indication towards the validity of a scaling scenario with large $`\mu `$. Certainly this scenario is more natural since it fits the data with (almost) no corrections to (simple) scaling, while the corrections for RSB are very large as apparent from Fig.5. A similar conclusion was reached in the determination of $`dd_s`$.

As a further test, we can obtain a second estimate of $`\mu `$ from the quantity $`1q_l`$ ($`q`$ unrestricted), whose scaling and asymptotic scaling are given in Eqs. (8) and (10). We find results similar to those for $`B`$: A power law fit $`1q_l=b/L^{\stackrel{~}{\mu _l}}`$ (Form III) gives acceptable fits for $`ϵ/\tau \frac{1}{2}`$ and for all values of $`ϵ`$ if we exclude $`L=4`$. As for $`B`$, the effective exponent $`\stackrel{~}{\mu _l}`$ changes significantly with $`ϵ`$, due to corrections to asymptotic scaling, and by extrapolating it to $`ϵ=0`$ we obtain $`\mu _l=1.15\pm 0.08`$. This gives $`\mu =\mu _l(dd_s)=0.60\pm 0.09`$ which agrees with the estimate $`\mu =0.63\pm 0.03`$ obtained from $`B`$. We also verified that, as for $`B`$, the data collapse reasonably well on one curve for $`\mu =0.64\pm 0.05`$ according to Eq. (8), although the quality of the scaling is somewhat worse than that of Fig.5. To check the RSB prediction, we fitted the data to Forms I and II (where now $`1q_l_c`$ is replaced by $`1q_l`$), finding that they both fit well the data, with $`a`$ roughly proportional to $`ϵ`$ as expected in RSB, although, for small $`ϵ`$, $`a`$ is also compatible with zero. Therefore, as for $`B`$, the data are also consistent with RSB, but this scenario requires large corrections to scaling, while the hypothesis $`\mu =0.63`$ fits the data with almost no corrections to (simple) scaling.

In the analysis so far, we have determined the exponents $`\mu `$ and $`dd_s`$ using just the link-overlap $`q_l`$. By contrast, PY determined $`\mu `$ (for periodic bc) from the scaling of the spin-overlap $`q`$. An advantage of $`q_l`$ is that its variance is much lower, as shown in Fig. 1. In any event, we have verified that the scaling relation Eq. (7) fits well the data for $`q`$, giving $`\mu =0.65\pm 0.02`$, in agreement with the estimates from $`B`$ and $`1q_l`$.

Summarizing this part, we find that the data for all the quantities considered, namely $`B`$, $`1q_l`$, and $`1q`$, are consistent with the RSB prediction that $`\mu =0`$ asymptotically, but large corrections to scaling are required in the fit, similarly to what we observed in the determination of $`dd_s`$. The data are also fitted very well by a scaling scenario with $`\mu 0.63`$, with almost negligible corrections to scaling (but with sizeable corrections to asymptotic scaling, which instead were small for the observables considered for $`dd_s`$). Under the “natural” assumption of small corrections to scaling, from the estimates of $`\mu `$ and $`dd_s`$ in Eqs. (29) and (26), we obtain

$$\theta ^{}=\mu (dd_s)=0.19\pm 0.06,$$
(30)
where, again, the error is purely statistical subject to the condition of having small (less than 3%) corrections to scaling. This result agrees with the droplet theory which predicts that $`\theta ^{}=\theta >0`$, and is compatible with the value of $`\theta 0.2`$ characterizing the energy of domain walls induced by a change in boundary conditionstheta-3d . By contrast, for periodic bc and under the same assumption of small corrections to scaling, $`\theta ^{}`$ is compatible with zero (see PY and Section III.3). In Section III.4 we will analyze the origin of this discrepancy, and show that by allowing small (of order 10$`\%`$) corrections to scaling the free bc data can be reconciled with $`\theta ^{}0`$.

### III.2 Box overlaps

So far we have analyzed the link- and spin-overlap which are computed on the whole system (bulk). We now turn to a different observable, the box-overlap defined as

$$q_B=\frac{1}{L_B^d}\underset{i}{}S_i^{(0)}\stackrel{~}{S}_i^{(0)}$$
(31)
where the sum runs over the sites contained in a central cubic box of fixed size $`L_B=2`$. In the following we will only consider the absolute value $`|q_B|`$, which we still call $`q_B`$ for simplicity. When a large-scale cluster of spins is flipped, for large $`L`$ the probability that its surface goes across the central box is proportional to the ratio of its surface area, $`L^{d_s}`$, to the volume, $`L^d`$. Therefore $`1q_BL^{(dd_s)}`$ from which we obtain the scaling laws

$`1q_B`$ $`=`$ $`L^{(dd_s)}F_{q_B}(ϵ/L^\mu )`$ (32)
$`1q_B_c`$ $`=`$ $`L^{(dd_s)}F_{q_B}^c(ϵ/L^\mu )`$ (33)
where, as for the corresponding scaling functions for $`q_l`$, $`F_{q_B}(x)x`$ and $`F_{q_B}^c(x)`$ const. for small $`x`$. Hence the asymptotic scaling for $`L\mathrm{}`$ is

$`1q_B`$ $``$ $`ϵL^{\mu _l}`$ (34)
$`1q_B_c`$ $``$ $`L^{(dd_s)}.`$ (35)
In RSB, this reduces to $`1q_Bϵ`$ and $`1q_B_c`$ const. An advantage of $`q_B`$ over $`q_l`$ is that the former, being measured away from the boundaries, should have smaller corrections to scaling and be less sensitive to boundary conditions. Indeed, Monte Carlo simulationsricci ; py\_unpub  show that $`q_B`$ has rather small corrections to scaling. This is likely to be particularly important for the free boundary conditions used here.

Fig. 6 shows the restricted average $`1q_B_c`$, with $`q_{\mathrm{max}}=0.4`$, as a function of $`L`$ for two values of $`ϵ`$. The data are clearly decreasing with $`L`$, are essentially independent of $`ϵ`$, as expected from Eq. (35), and are close to a straight line on the logarithmic plot. This indicates that the power law fit, Form III, appropriate to the droplet and TNT scenarios, should work well and indeed it does, even for the largest value of $`ϵ`$ (we note however that the statistical errors are larger than for the link-overlap, hence the fits are less sensitive to corrections to scaling). The exponent is almost independent of $`ϵ`$, varying between 0.48 and 0.52, and from this we obtain the estimate

$$dd_s=0.48\pm 0.03$$
(36)
which is in agreement with the estimates $`dd_s=0.44\pm 0.03`$ from $`1q_l_c`$ and $`dd_s=0.43\pm 0.03`$ from $`R`$.

Forms I and II (with $`q_B`$ replacing $`q_l`$) also fit well the data, with $`a`$ between 0.14 and 0.36 (with no discernible trend with $`ϵ`$). Hence the data are also compatible with RSB, and under the RSB hypothesis, we estimate

$$\underset{L\mathrm{}}{lim}1q_B_c=0.25\pm 0.10\text{(RSB)}.$$
(37)
As usual, we note that the RSB scenario requires rather large corrections to scaling, while the power law fits the data with no corrections.

Fig. 7 shows the unrestricted average $`1q_B`$ multiplied by $`\tau /ϵ`$, which asymptotically should be independent of $`ϵ`$. The data show a small curvature and a significant $`ϵ`$ dependence, indicating that for this quantity we are not yet in the asymptotic scaling regime (similarly to what we observed for the quantity $`B`$). The data are fitted well by a power law, with an exponent that changes with $`ϵ`$ and tends towards $`\mu 0.63`$ for $`ϵ0`$. Fits using Forms I and II give $`a`$ compatible with zero. We also determined $`\mu `$ from the scaling relation Eq. (32), by fixing $`dd_s=0.44`$ and using the same fitting procedure as for $`B`$ (which assumes no corrections to scaling), finding

$$\mu =0.62\pm 0.04$$
(38)
which agrees with the various estimates of $`\mu `$ obtained from $`B`$, $`1q_l`$, and $`1q`$. Fig. 8 shows the corresponding scaling plot, in which the data collapse is fairly good.

To conclude this subsection, the data for box overlaps can be fitted with smaller corrections to scaling than the data for the bulk link- and spin-overlap. A fit to the generic scaling picture, with no corrections to scaling, gives results for the exponents $`dd_s`$ and $`\mu `$ in agreement with those from the bulk quantities analyzed in the previous subsections. However, as with the bulk observables, assuming large corrections to scaling, the data can also be fitted to the RSB picture.

### III.3 Comparison with periodic boundary conditions

In order to assess the effect of different boundary conditions, we have repeated part of the analysis above (with the exclusion of box-overlaps) for the data of PY (Ref.py-bulk, ) for periodic boundary conditions and $`L8`$. The ground states were obtained using a hybrid genetic algorithm as described in PY. This does not guarantee to find the true ground state, but the systematic errors due to occasionally missing it are smaller than the statistical errorspy-bulk .

Figs. 9 and 10 show the equivalent for periodic bc of Figs. 4 and 5 for free bc. The data in Fig. 9 shows much less curvature and also a smaller dependence on $`ϵ`$ than the corresponding data for free bc in Fig. 4, indicating that corrections to scaling are smaller for periodic bc. Table 3 reports the best fits using the three functional forms of Eq. (23). Form I fits well the data, but $`a`$ varies significantly with $`ϵ`$, and for small $`ϵ`$ it is compatible with zero. Form II also fits well, with $`a`$ independent of $`ϵ`$ within the statistical errors. From this fit we estimate

$$\underset{L\mathrm{}}{lim}1q_l_c=0.28\pm 0.03\text{(RSB)}$$
(39)
(see the comment after Eq.(24) as to the meaning of the error bar) which agrees with the estimate 0.24 of Marinari and Parisimapabulk , and is just consistent with our estimate $`0.20\pm 0.02`$ for free bc.

The power-law fit with no corrections, Form III, fits well the data for the two smallest values of $`ϵ`$ and, if we exclude $`L=3`$, for all but the largest value of $`ϵ`$. The exponent $`cdd_s`$ is nearly independent of $`ϵ`$ and gives

$$dd_s=0.43\pm 0.02\text{(periodic bc)}.$$
(40)
This result agrees with the estimate $`dd_s=0.42\pm 0.02`$ of PY obtained from the ratio $`R`$ defined above, confirming that corrections due to small droplets should not be important in three dimensions, and with our estimate $`dd_s=0.44\pm 0.03`$ for free bc, indicating that $`dd_s`$ does not depend on boundary conditions.

We also performed fits with Form IV which includes corrections to scaling. As for free bc, a wide range of values of $`dd_s`$ from zero to around 0.44 give a good fit, with the largest values giving the smallest corrections to scaling. The results of the fit for $`dd_s=0.43`$ are shown in Table 3. For the two smaller values of $`ϵ`$, the fits are difficult because corrections to scaling are very small, hence they are not shown.

We determined the exponent $`\mu `$ from the ratio $`B`$ using Eq. (28) and the fitting procedure described for free bc, obtaining

$$\mu =0.42\pm 0.03\text{(periodic bc)}$$
(41)
where, as for the estimate of $`dd_s`$ above, the errors are purely statistical with the assumption that corrections to scaling are smaller than the statistical errors of the data. Scaling is rather satisfactory as shown in Fig. 10. This value agrees with the estimate $`\mu =0.44\pm 0.02`$ of PY from the scaling of the spin-overlap but incompatible, within the statistical error bars, with the result $`\mu =0.63\pm 0.03`$ for free bc. We will return in Section III.4 on the origin of the discrepancy between free and periodic bc. The inset of Fig. 9 shows that, with these values of $`\mu `$ and $`dd_s`$, the scaling form for $`1q_l_c`$, Eq. (14), is also well satisfied. Finally, we verified that, if $`dd_s=0.43`$, the unrestricted average $`1q_l`$ satisfies scaling, giving $`\mu =0.45\pm 0.02`$ in agreement with the estimate from $`B`$.

Combining Eqs. (40) and (41), we obtain the estimate of $`\theta ^{}`$ for periodic bc:

$$\theta ^{}=\mu (dd_s)=0.01\pm 0.03\text{(periodic bc)}.$$
(42)
This is compatible with zero and, within the error bars, incompatible with the value $`\theta ^{}=\theta 0.2`$, where $`\theta `$ characterizes the energy of domain walls induced by boundary condition changes. A scenario in which $`\theta ^{}=0`$ and $`dd_s>0`$ is consistent with the TNT picture. Finally, we note that, although our analysis of the PY data uses different quantities to extract exponents, our results agree with those given by PY.

### III.4 Discussion and summary of the results

In the previous sections, we have seen that for both free and periodic bc, the analysis of all the different observables considered gives consistent results for the exponents $`dd_s`$ and $`\theta ^{}`$ under the assumption of minimal corrections to scaling. However, while the results for $`dd_s`$ for free and periodic bc agree with each other, the results for $`\theta ^{}`$ apparently do not, having found $`\theta ^{}=0.01\pm 0.03`$ for periodic bc and $`\theta ^{}=0.19\pm 0.06`$ for free bc. Since $`\theta ^{}`$, like $`dd_s`$, should not depend on the type of boundary conditions, the discrepancy must be due to different corrections to scaling for the two boundary conditions.

Therefore, it is important to analyze further the corrections to scaling. First, we recall that the scaling plots for the quantity $`B`$ in Figs. 5 (free bc) and 10 (periodic bc), from which we have determined $`\mu `$ (and hence $`\theta ^{}`$), are obtained by imposing that the whole data set (namely all values of $`ϵ`$ and $`L`$) satisfies scaling with corrections to scaling smaller than the statistical errors, which are less than $`1\%`$. Clearly, this is a very stringent requirement. If we relax this requirement, allowing some corrections to (simple) scaling, we can accommodate a larger range of values for $`\mu `$.

This is shown in Fig. 11, which gives a scaling plot for free bc, analogous to Fig. 5 but assuming the value $`\mu =0.42`$ determined from periodic bc. The polynomial fitting curve was obtained by excluding from the fit the data points for $`L=4`$ and 6. One can see that for larger $`L`$ the data collapse reasonably well on the curve. The deviation of the $`L=4,6`$ data from the curve, less than $`10\%`$, is a measure of the corrections to scaling. Therefore, we see that corrections to scaling of less than $`10\%`$ for the two smallest sizes are sufficient to remove the discrepancy in the value of $`\theta ^{}`$ between free and periodic bc. We verified that the also the other quantities considered, namely $`1q_l_c`$, $`R`$, $`1q_l`$, $`1q`$, can be fitted in a similar way.

We also tried the converse operation, namely a scaling plot of the data for periodic bc but using the value $`\mu =0.63`$. We find that one can get a relatively good data collapse excluding the sizes $`L=3,4`$, and 5, which deviate from the scaling curve by less than $`10\%`$. However, now the data for a given $`ϵ`$ approach the scaling curve from the right side instead of from the left side as in Fig. 11, but since they have an upward curvature, the correction to scaling should change curvature twice as $`L`$ increases, which is not very plausible. Hence, we believe that it is more natural to conclude that the correct value of $`\mu `$ is closer to 0.42 than to 0.63, namely that corrections to scaling are smaller for periodic than for free.

Indeed, in general, it is reasonable to expect that corrections are larger for free bc, because these bc have a free surface on which lie a fraction of sites which is quite substantial for moderate sizes. In Fig. 12 we plot together the $`1q_l_c`$ data of Fig. 4 and 9 for free and periodic bc. The data for free bc lie significantly below those for periodic bc, indicating that the surface of the excitations is smaller for free bc. For periodic bc, the domain wall has to “bend” to return to the same point on the “top surface” as it had on the “bottom surface”. This may be the source of the extra surface area.

Under the hypothesis $`dd_s0.44`$, Fig. 12 then shows that the corrections to asymptotic scaling are larger for free bc, since the free bc data show a more marked deviation from the asymptotic $`ϵ`$-independent behavior, and display a larger curvature. This is further indication that free bc have larger corrections.

Evidence that free bc have larger corrections was also found recently in Monte Carlo simulationshelmut , where some evidence was observed that the free bc data might have a crossover from droplet-like to either TNT- or RSB-like behavior at large sizes.

Incidentally, note that if RSB is the correct asymptotic picture and the $`L\mathrm{}`$ limit of $`1q_l_c`$ is the same for periodic and free bc, then Fig. 12 would indicate that the corrections are smaller for free bc (since the data are closer to their non-zero asymptotic value) which is not very plausible. Note, however, that we do not have an argument why in the thermodynamic limit $`1q_l_c`$ should be independent of boundary conditions.

To summarize the first part of the article, we have analyzed several quantities for periodic and free bc. For both types of boundary conditions, all the data are well described by a general scaling picture involving only two scaling exponents, $`dd_s`$ and $`\theta ^{}`$, with only small corrections to scaling. Some observables show significant corrections to asymptotic scaling, which are larger for free boundary conditions. Fitting this scaling picture to our data, we obtain comparable values of $`dd_s`$ for periodic ($`0.43\pm 0.02`$) and free boundary conditions ($`0.44\pm 0.03`$).

By imposing that corrections to scaling are less than the statistical errors of $`1\%`$, for periodic boundary conditions we obtain $`\theta ^{}0`$, which fits well the TNT scenario ($`dd_s>0,\theta ^{}=0`$), while for free boundary conditions we obtain $`\theta ^{}=0.19\pm 0.06`$, which fits well the droplet picture ($`dd_s>0,\theta ^{}>0`$). By relaxing this requirement and allowing larger corrections to scaling of order 10$`\%`$, the data for free bc can be also fitted by a scenario with $`\theta ^{}0`$. Therefore the data for free bc are also consistent with the TNT picture provided moderate corrections to scaling are allowed, larger than those for periodic bc. We have also provided direct evidence that indeed free bc have larger corrections to scaling.

Data for the box overlap for free boundary conditions indicates smaller corrections to asymptotic scaling, which is reasonable since the box is away from the surface, and are consistent with the scenario described above.

For both free and periodic bc, the data are also fitted well by the RSB picture ($`dd_s=0,\theta ^{}=0`$), but only if we allow very large corrections to scaling. In this case, the good scaling behavior we observed for all the observable considered would only be a finite size artifact, and would disappear at larger sizes. To test this possibility, large system sizes will be needed.

This concludes the first part of the paper, dedicated to the physical results. In the second part, we will describe the branch and cut algorithm employed, and analyze its performance in our computations.

## IV Branch and cut algorithm

Branch and cut is, to our knowledge, the fastest exact method for determining ground states of spin glasses in three dimensions. To apply this technique, we transform the problem of minimizing the Hamiltonian in Eq.(1) into a standard combinatorial optimization problem known as the maximum cut problem. (For a detailed description of optimization and related topics, see Ref. handbook, .) Consider the interaction graph $`G=(V,E)`$ associated with the spin glass Hamiltonian, where $`G`$ contains vertices $`1,\mathrm{},L^3V`$ associated with the spin sites and edges $`(ij)E`$ with weight $`c_{ij}=J_{ij}`$ associated with the couplings.

Given a partition of $`V`$ into two sets, $`WV`$ and its complement $`VW`$, the cut $`\delta (W)`$ associated with $`W`$ is defined as the set of edges with one endpoint, $`i`$ say, in $`W`$ and the other endpoint, $`j`$ say, in $`VW`$. In formulas, $`\delta (W)=\{(ij)EiW,jVW\}`$. The weight of a cut $`\delta (W)`$ is defined as the sum of the weights of the cut edges $`_{(ij)\delta (W)}c_{ij}`$. A maximum cut is a cut with maximum weight among all partitions $`W`$. It is easy to show that minimizing the Hamiltonian Eq.(1) is equivalent to finding a maximum cut in $`G`$, see Ref. desim, . If we know a maximum cut with node partition $`W`$ and $`VW`$, the corresponding ground state spin configuration can be read off by assigning the value up to the spins in $`W`$ and down to the spins in $`VW`$, or vice versa.

The branch-and-cut algorithm solves the maximum cut problem through simultaneous lower and upper bound computations. By definition, the weight of any cut gives a lower bound on the optimal cut value. Thus, we can start from any cut and iteratively improve the lower bound using deterministic heuristic rules (local search and other specialized heuristics, see Ref.juen-vlsi,  for details). How do we decide when a cut is optimal? This can be done by additionally maintaining upper bounds on the value of the maximum cut. Upon iteration of the algorithm, progressively tighter bounds are found, until optimality is reached.

Since the availability of upper bounds marks the difference between a heuristic and an exact solution, we now summarize how the upper bound is computed (for more details, see Ref. juen-vlsi, .) To each edge $`(ij)`$ we associate a real variable $`x_{ij}`$ and to each cut $`\delta (W)`$ an incidence vector $`\chi ^{\delta (W)}^E`$ with components $`\chi _{ij}^{\delta (W)}`$ associated to each edge $`(ij)`$, where $`\chi _{ij}^{\delta (W)}=1`$ if $`(ij)\delta (W)`$ and $`\chi _{ij}^{\delta (W)}=0`$ otherwise. Denoting by $`P_C(G)`$ the convex hull of the incidence vectors, it can be shown that a basic optimum solutionchvatal  of the linear program

$$\mathrm{max}\{\underset{(ij)E}{}J_{ij}x_{ij}xP_C(G)\}.$$
(43)
is a maximum cut. In order to solve (43) with linear programming techniques we would have to express $`P_C(G)`$ in the form

$$P_C(G)=\{x^EAxb,0x1\}$$
(44)
for some matrix $`A`$ and some vector $`b`$. Whereas the existence of $`A`$ and $`b`$ are theoretically guaranteed, even subsets of $`Axb`$ known in the literature contain a huge number of inequalities that render a direct solution of (43) impractical.

Instead, the branch-and-cut algorithm proceeds by optimizing over a superset $`P`$ containing $`P_C(G)`$, and by iteratively tightening $`P`$, generating in this way progressively better upper bounds. The supersets $`P`$ are generated by a cutting plane approach. Starting with some $`P`$, we solve the linear program $`\mathrm{max}\{_{(ij)E}J_{ij}x_{ij}xP\}`$ by Dantzig’s simplex algorithmchvatal . Optimality is proven if either of two conditions is satisfied: (i) the optimal value equals the lower bound; (ii) the solution vector $`\overline{x}`$ is the incidence vector of a cut.

If neither is satisfied, we have to tighten $`P`$ by solving the separation problem. This consists in identifying inequalities that are valid for all points in $`P_C(G)`$, yet are violated by $`\overline{x}`$, or reporting that no such inequality exists. The inequalities found in this way are added to the linear programming formulation, obtaining a new tighter partial system $`P^{}P`$ which does not contain $`\overline{x}`$. The procedure is then repeated on $`P^{}`$ and so on.

At some point, it may happen that (i) and (ii) are not satisfied, yet the separation routines do not find any new cutting plane. In this case, we branch on some fractional edge variable $`x_{ij}`$ (i.e. a variable $`x_{ij}\{0,1\}`$), creating two subproblems in which $`x_{ij}`$ is set to 0 and 1, respectively. We then we apply the cutting plane algorithm recursively for both subproblems.

## V Performance of the branch and cut algorithm

In this section we study the performance of our current implementation of the branch and cut algorithm, in particular the dependence of the number of computer operations on system size. The results for size $`L=12`$ were obtained with a more efficient version of the code, so performance for this size cannot be compared with that for the smaller sizes. Hence, in this section, we shall just consider sizes up to $`L=10`$.

Finding the ground state of the Hamiltonian Eq. (1) in three dimensions is an $`𝒩𝒫`$-hard problembarahona , and all known algorithms to solve this class of problems require a number of operations that grows exponentially on the size of the input, in the worst case.

However, depending on the problem, the number of operations for typical instances (for the spin glass problem, an instance is a realization of the random couplings, or sample), can grow considerably more slowly than the worst case exponential behavior. Furthermore, the number of operations can vary significantly from one instance to another. It is therefore useful to investigate experimentally the performance of the algorithm for typical instances, in order to try to extrapolate the computational resources necessary to go to larger sizes, and possibly to identify which parameters of the problem affect most the performance. De Simone et al.desim  measured the average CPU time used by the branch and cut algorithm to find the ground state of the two-dimensional $`\pm J`$ spin glass with periodic bc, up to $`L=70`$, showing that the average CPU time was approximated by a function proportional to $`L^6`$.

Here we analyze the performance of the branch and cut algorithm for the three-dimensional spin glass with free bc and Gaussian couplings. In order to do this, we first need a good measure of the performance. For a complex algorithm such as branch and cut, a simple and absolute measure of the number of operations is not available. Two possible measures are the CPU time and the number of linear programs solved during the run of the algorithm. In Table 4, we summarize the average running time needed for calculating an unperturbed ground state for the different system sizes.

The CPU time is not an accurate measure since it depends on the machine architecture and load. Furthermore, our computations were carried out on several different machines, so the CPU time is not useful here. We take instead the number of linear programs solved, $`n_p`$, because (i) it is a well-defined and machine independent quantity; (ii) we have observed that about $`95\%`$ of the time is spent in solving linear programs; (iii) for a fixed system size, $`n_p`$ correlates strongly, and almost linearly, with the CPU time. This is shown in Fig. 13, which plots the CPU time versus $`n_p`$ for 1000 randomly generated samples with $`L=10`$, computed on the same machine. Note that since the size of the linear programs is also growing with the system size, the CPU time per linear program increases strongly with $`L`$: the average (resp. median) CPU time goes from 0.00770 (resp. 0.833) seconds for $`L=4`$ to 0.833 (resp. 0.784) seconds for $`L=10`$.

Hence, $`n_p`$ severely underestimates the rate at which the number of operations increases with $`L`$.

From Fig. 13, we also note that the distribution of $`n_p`$ (and CPU times) is very broad, extending over three orders of magnitude. The histogram distribution of $`n_p`$ for different system sizes is shown in Fig. 14. In addition to shifting to larger $`n_p`$, the distribution broadens as $`L`$ increases. Also, there is some evidence of a double-peak structure. For $`L=10`$, we verified that the peak at smaller $`n_p`$ corresponds to samples that could be solved without branching, while the other peak corresponds to samples where branching was necessary. Since in each branching step the number of subproblems to be solved doubles, the number of linear programs increases rapidly and the second peak is at significantly larger $`n_p`$.

In order to identify which parameters of the problem, in addition to the size, affect the performance, we ask whether $`n_p`$ correlates with the physical observables we measure. No significant correlation was observed with the ground state energy. Fig. 15 plotsnote\_scatter  $`\mathrm{log}_{10}n_p`$ for the unperturbed ground state ($`ϵ=0`$) and $`L=10`$ versus the overlap between this state and the perturbed ground state with $`ϵ/\tau =4`$. We observe a distinct correlation between $`n_p`$ and $`q`$: for small $`q`$, more linear programs are needed than for large $`q`$. The figure shows that the typical number of linear programs is close to an order of magnitude larger if $`q0`$ than if $`q1`$. We observed a similar correlation for other values of $`ϵ`$ as well, and also between the CPU time and $`q`$. Again, the distribution of $`n_p`$ is quite broad as shown by the data for the standard deviation of $`\mathrm{log}_{10}n_p`$ in Fig. 15.

In order to quantify how the correlation between $`n_p`$ and $`q`$ changes with the system size, we show in Fig. 16 the average and median of $`n_p`$ as a function of $`N_b`$, as well as the conditional averages of $`n_p`$ restricted to samples with large ($`|q|0.9`$) and small ($`|q|0.1`$) overlap. We take the number of bonds, $`N_b`$, as a measure of the input size, since the maximum cut problem is specified in terms of the edge variables in the graph. From Fig. 16 we see that, first, all measures show an approximately exponential increase with $`N_b`$, with corrections for small $`N_b`$, and second, the difference between the conditional averages with small and large $`q`$ seems to increase with the system size, and is about one order of magnitude for $`L=10`$.

A qualitative difference between samples with small and large overlap is that samples with a small $`|q|`$ have a rougher “energy landscape”, namely states with an energy close to the ground state energy yet a spin configuration very different from the ground state. It is then intuitively clear why one would observe a correlation between $`q`$ and the running time for a stochastic algorithm employing local search heuristics, such as simulated annealing, since when the algorithm encounters one of these configurations with small overlap, it must retrace its steps by a large amount.

For the branch and cut algorithm, the reason for the correlation between $`n_p`$ and $`q`$ is less obvious, but some insight is provided by an analysis of “reduced cost fixing”. This is a feature of the branch and cut algorithm speeding up the computations. In every iteration of the algorithm, reduced cost fixing gives us a sufficient condition to decide which variables (corresponding to the edges in the graph) have already attained their optimal value. Therefore, we can fix the variables with “optimal” status to their current value for all the subsequent iterations of the algorithm, resulting in less overall computational effort. The more variables that can be fixed, the faster the algorithm is in practice.

Since the samples with small overlap have “almost optimal” solutions with spin configurations very different from the ground state, a smaller number of variables can be fixed. Here we do not have the “correct” edge values available until the end. As an example, we checked that for $`L=10`$ and $`ϵ=\tau `$, for 100 randomly chosen samples with small overlap ($`|q|0.1`$), in average $`409\pm 39`$ of the 2700 edge variables could be fixed in the first sub problem, i.e. before branching takes place. In contrast, for 100 randomly chosen samples with big overlap ($`|q|0.9`$), $`921\pm 34`$ of the edge variables could be fixed in the first sub problem, about twice as many. Of course, the less variables that can be fixed in the first sub problem, the more overall branching is necessary, resulting in more overall computational effort for samples with small overlap.

A consequence of the the broad distribution of the CPU time and of its correlation with the physical observables of interest, is that a cutoff in the CPU time produces a systematic error in these quantities. One has therefore to ensure that the cutoff is large enough so that the systematic error is smaller than the statistical error.

It is interesting to try to extrapolate the running time needed to deal with larger sizes. The average CPU time in Table 4 varies approximately as $`\mathrm{exp}(\alpha N_b)`$ with $`\alpha `$ somewhere between 0.0024 and 0.003. Extrapolating to $`L=14`$ ($`N_b=7644`$), this gives an average CPU time of around $`10^{8\pm 1}`$ seconds per sample, which is clearly very demanding. Furthermore, memory limitations will set in before we can reach this size. Again, note that $`n_p`$ increases much more slowly with $`N_b`$. The data for $`|q|0.1`$ in Fig.16, for example, vary approximately as $`\mathrm{exp}(\alpha N_b)`$ with a smaller $`\alpha `$ around 0.0017, showing that the dominant limiting factor is the solution of the linear programs. Note that the program used for $`L=12`$ is significantly faster than that used in this extrapolation. This long extrapolated running time gives us further motivation to continue our research on the improvement of this algorithm.

## VI Conclusions

Using an exact “branch and cut” optimization algorithm, we have studied the large-scale, low-energy excitations in the Ising spin glass in three dimensions with free boundary conditions, and compared the results with those obtained earlier by PY for periodic boundary conditions.

In the first part of the paper, we have discussed in detail how the whole set of observables analyzed is fitted by a general scaling picture characterized by two exponents, $`dd_s`$ and $`\theta ^{}`$, and how the values of these parameters predicted by the various physical pictures proposed for the spin glass phase fit our data. Our conclusions have been summarized at the end of Section III.4.

In the second part of the paper, we have analyzed the performance of the branch and cut algorithm, finding that the performance is worse when there is a low energy excited state close in energy to the ground state but far away in configuration space, and have given a quantitative analysis of this effect.

###### Acknowledgements.

We would like to thank A. J. Bray, G. Parisi, M. Mézard, D. S. Fisher, and M. A. Moore for helpful discussions and correspondence. APY acknowledges support from the NSF through grant DMR 0086287. MP would like to thank A. J. Bray for a useful suggestion on the data analysis. We thank the Regional Centre of Computing of the University of Cologne for the allocation of computer time. Over the years, Giovanni Rinaldi and Gerd Reinelt contributed much to the algorithm.


format: text

text:
# Certain locally nilpotent varieties of groups

## 1. Introduction and results

Let $`c0`$, $`d2`$ be integers and $`𝒩_c`$ be the variety of nilpotent groups of class at most $`c`$. We denote by $`𝒩_c^{(d)}`$ the variety of groups in which every $`d`$-generator subgroup is in $`𝒩_c`$. In , Gupta posed the following question:
For what values of $`c`$ and $`d`$ it is true that $`𝒩_c^{(d)}`$ is locally nilpotent?
Then he proved that for $`c(d^2+2d3)/4`$, the variety $`𝒩_c^{(d)}`$ is locally nilpotent. In , Endimioni improved the latter result where he proved that for $`c2^d2`$, the variety $`𝒩_c^{(d)}`$ is locally nilpotent. Here we improve the number $`2^d2`$ to $`2^d+2^{d1}3`$. In fact we prove:

###### Theorem 1.1.

(1) For $`c2^d+2^{d1}3`$, the variety $`𝒩_c^{(d)}`$ is locally nilpotent.
(2) For $`c2^d+2^{d1}+2^{d2}3`$, every $`p`$-group in the variety $`𝒩_c^{(d)}`$ is locally nilpotent, where $`p\{2,3,5\}`$.

Note that the variety $`𝒩_c^{(2)}`$ is contained in the variety of $`c`$-Engel groups and it is yet unknown whether every $`c`$-Engel group is locally nilpotent, even, so far there is no published example of a non-locally nilpotent group in the variety $`𝒩_c^{(2)}`$. In the last section of this paper, by considering the problem of locally nilpotency of the variety $`𝒩_c^{(2)}`$, we study periodic groups in this variety. Note that since every two generator subgroup of a group in $`𝒩_c^{(2)}`$ is nilpotent, every periodic group in $`𝒩_c^{(2)}`$ is a direct product of $`p`$-groups ($`p`$ prime). We reduce the question of Gupta for periodic groups in $`𝒩_c^{(d)}`$ to the locally nilpotency of $`p`$-groups of finite exponent in this variety where the exponent depends only on the numbers $`p`$ and $`c`$. In fact we prove that

###### Theorem 1.2.

Let $`p`$ be a prime, $`c>1`$ an integer and $`r=r(c,p)`$ be the integer such that $`p^{r1}<c1p^r`$. Then the following are equivalent:
(1) every $`p`$-group in $`𝒩_c^{(2)}`$ is locally nilpotent.
(2) if $`p`$ is odd, every $`p`$-group of exponent dividing $`p^r`$ in $`𝒩_c^{(2)}`$ is locally nilpotent, and if $`p=2`$, every $`2`$-group of exponent dividing $`2^{r+1}`$ in $`𝒩_c^{(2)}`$ is locally nilpotent.

## 2. Groups in the variety $`𝒩_c^{(d)}`$

Let $`F_{\mathrm{}}`$ be the free group of infinite countable rank on the set $`\{x_1,x_2,\mathrm{}\}`$, we define inductively the following words in $`F_{\mathrm{}}`$:

$$W_1=W_1(x_1,x_2)=[x_1,x_2,x_1,x_2],$$

$$W_n=W_n(x_1,x_2,\mathrm{},x_{n+1})=[W_{n1},x_{n+1},W_{n1},x_{n+1}]n>1;$$

$$V_1=V_1(x_1,x_2,x_3)=[[x_2,x_1,x_1,x_1,x_1],x_3,[x_2,x_1,x_1,x_1,x_1],x_3],$$

$$V_n=V_n(x_1,x_2,x_3,\mathrm{},x_{n+2})=[V_{n1},x_{n+2},V_{n1},x_{n+2}]n>1.$$

For a group $`G`$ and a subgroup $`H`$ of $`G`$, we denote by $`HP(G)`$ the Hirsch-Plotkin radical of $`G`$ and $`H^G`$ the normal closure of $`H`$ in $`G`$. We use the following result due to Heineken (see Lemma 8 of  and see Lemma 2 of  for the left-normed version).

###### Lemma 2.1.

Let $`G`$ be a group and $`g`$ an element in $`G`$ such that $`[g,x,g,x]=1`$ for all $`xG`$. Then the normal closure of $`g`$ in $`g^G`$ is abelian. In particular, $`gHP(G)`$.

###### Lemma 2.2.

Let $`G`$ be a group satisfying the law $`W_n=1`$ for some integer $`n1`$. Then, $`G`$ has a normal series $`1=G_nG_{n1}\mathrm{}G_1=G`$ in which each factor $`G_i/G_{i+1}`$ is locally nilpotent ($`i=1,2,\mathrm{},n1`$).

###### Proof.

We argue by induction on $`n`$. If $`n=1`$, then Lemma 2.1 yields that $`x_1HP(G)`$ for all $`x_1G`$ and so $`G`$ is locally nilpotent. Now suppose that the lemma is true for $`n`$ and $`G`$ satisfies the law $`W_{n+1}=1`$. By Lemma 2.1, we have $`W_n(x_1,\mathrm{},x_{n+1})HP(G)`$ for all $`x_1,\mathrm{},x_{n+1}G`$ and so $`G/HP(G)`$ satisfies the law $`W_n=1`$. Thus by induction hypothesis $`G/HP(G)`$ has a normal series of length $`n`$ with locally nilpotent factors, it completes the proof. ∎

###### Lemma 2.3.

Let $`G`$ be a $`p`$-group satisfying the law $`V_n=1`$ for some integer $`n1`$ where $`p\{2,3,5\}`$. Thus $`G`$ has a normal series $`1=G_{n+1}G_n\mathrm{}G_1=G`$ in which each factor $`G_i/G_{i+1}`$ is locally nilpotent ($`i=1,2,\mathrm{},n`$).

###### Proof.

We argue by induction on $`n`$. If $`n=1`$, then by Lemma 2.1
$`[x_2,x_1,x_1,x_1,x_1]HP(G)`$ for all $`x_1,x_2G`$. Thus $`G/HP(G)`$ is a $`4`$-Engel group and since every $`4`$-Engel $`p`$-group is locally nilpotent where $`p\{2,3,5\}`$ (see Traustason  and Vaughan-Lee ), so $`G/HP(G)`$ is locally nilpotent. Now suppose that $`G`$ satisfies the law $`V_{n+1}=1`$. By Lemma 2.1, $`V_n(x_1,\mathrm{},x_{n+2})HP(G)`$ for all $`x_1,\mathrm{},x_{n+2}G`$, so $`G/HP(G)`$ satisfies the law $`V_n=1`$. Thus by induction hypothesis, $`G/HP(G)`$ has a normal series of length $`n+1`$ with locally nilpotent factors, it completes the proof. ∎

We use in the sequel the following special case of this well-known fact due to Plotkin  that every Engel radical group is locally nilpotent. (see also Lemma 2.2 of )

###### Lemma 2.4.

Let $`H`$ be a normal subgroup of an Engel group $`G`$. If $`H`$ and $`G/H`$ are locally nilpotent, then $`G`$ is locally nilpotent.

Proof of Theorem 1.1. One can see that $`W_n`$ and $`V_n`$ are, respectively, in the $`(2^{n+1}+2^n2)`$th term and $`(2^{n+2}+2^{n+1}+2^n2)`$th term of the lower central series of $`F_{\mathrm{}}`$, for all integers $`n1`$.
(1) every group $`G`$ in the variety $`𝒩_c^{(d)}`$ satisfies the law $`W_{d1}=1`$ and so it follows from Lemmas 2.2 and 2.4, that $`G`$ is locally nilpotent.
(2) if $`d=2`$ then $`c4`$ and so every group in the variety $`𝒩_c^{(d)}`$ is $`4`$-Engel. But as it is mentioned in the proof of Lemma 2.3, every $`4`$-Engel $`p`$-group is locally nilpotent, where $`p\{2,3,5\}`$. Now assume that $`d3`$, then every group $`G`$ in the variety $`𝒩_c^{(d)}`$ satisfies the law $`V_{d2}=1`$ and so it follows from Lemmas 2.3 and 2.4, that $`G`$ is locally nilpotent. $`\mathrm{}`$

## 3. $`p`$-groups in the variety $`𝒩_c^{(2)}`$

###### Lemma 3.1.

Let $`c1`$ be an integer, $`p`$ a prime number and $`G`$ a finite $`c`$-Engel $`p`$-group. Suppose that $`x,yG`$ such that $`x^{p^n}=y^{p^n}=1`$ for some integer $`n>0`$. Let $`r`$ be the integer such that $`p^{r1}<cp^r`$. Then:
(a) if $`p`$ is odd and $`n>r`$ then $`[x^{p^{n1}},y^{p^{n1}}]=1`$.
(b) if $`p=2`$ and $`n>r+1`$ then $`[x^{2^{n1}},y^{2^{n1}}]=1`$.

###### Proof.

Suppose that $`KH`$ are two normal subgroups of $`G`$ such that $`H/K`$ is elementary abelian and $`a`$ is an arbitrary element of $`G`$. Put $`t=aK`$ and $`V=H/K`$. Since $`[V,_ct]=1`$, we have that $`[V,_{p^r}t]=1`$ and $`0=(t1)^{p^r}=t^{p^r}1`$ in $`\text{End}(V)`$. Thus $`[H,a^{p^r}]K`$ for all $`aG`$. Now let $`N=x^{p^r},y^{p^r}`$. Then $`[H,N]K`$ and since $`K,H`$ are normal in $`G`$; $`[H,M]K`$ where $`M=N^G`$ the normal closure of $`N`$ in $`G`$. Thus $`M`$ is a normal subgroup of $`G`$ centralizing every abelian normal section of $`G`$. By a result of Shalev ; if $`p`$ is odd then $`M`$ is powerful and if $`p=2`$ then $`M^2`$ as well as all subgroups of $`M^2`$ which are normal in $`G`$, are powerful. Suppose that $`p`$ is odd. Since $`M`$ is generated by $`\{(x^g)^{p^r},(y^g)^{p^r}|gG\}`$ and $`M`$ is powerful, by Corollary 1.9 of , $`M^{p^{nr}}`$ is generated by $`\{(x^g)^{p^n},(y^g)^{p^n}|gG\}`$ and so $`M^{p^{nr}}=1`$. On the other hand $`M^{p^{nr1}}`$ is powerful by Corollary 1.2 of . Thus $`[M^{p^{nr1}},M^{p^{nr1}}](M^{p^{nr1}})^p`$. Now by Theorem 1.3 of , we have $`(M^{p^{nr1}})^p=M^{p^{nr}}=1`$. Thus $`M^{p^{nr1}}`$ is abelian and the part (a) has been proved. Now assume that $`p=2`$. As it is mentioned in the above, since the subgroup $`x^{2^{r+1}},y^{2^{r+1}}^G`$ of $`M`$ is normal in $`G`$, it is also powerful and the rest of the proof is similar to the latter case. ∎

###### Lemma 3.2.

Let $`c>1`$ be an integer, $`p`$ a prime number and $`G`$ a $`p`$-group in the variety $`𝒩_c^{(2)}`$ and let $`r`$ be the integer satisfying $`p^{r1}<c1p^r`$.
(a) if $`p`$ is odd then $`G^{p^r}`$ is locally nilpotent.
(b) if $`p=2`$ then $`G^{2^{r+1}}`$ is locally nilpotent.

###### Proof.

Note that $`G`$ is a $`c`$-Engel group. Suppose $`p`$ is odd. First we prove that if $`x`$ is an element of $`G`$ such that $`x^{p^n}=1`$ for some $`n>r`$, then $`[a,x^{p^{n1}},x^{p^{n1}}]=1`$ for all $`aG`$. Let $`y=(x^1)^a`$, then it is enough to show that $`[y^{p^{n1}},x^{p^{n1}}]=1`$. Since $`x,a𝒩_c`$ then $`y,x𝒩_{c1}`$; therefore by Lemma 3.1, we have that $`[y^{p^{n1}},x^{p^{n1}}]=1`$. Now let $`z`$ be an arbitrary element of $`G`$ such that $`z^{p^n}=1`$ for some positive integer $`n`$. We prove by induction on $`n`$ that $`z^{p^r}HP(G)`$ and so it completes the proof for the case $`p`$ odd. If $`nr`$, then $`z^{p^r}=1HP(G)`$. Assume that $`n>r`$ then by induction hypothesis, $`z^{p^{r+1}}HP(G)`$. Now by the first part of the proof, $`[a,z^{p^r},z^{p^r}]=1`$ mod $`HP(G)`$, for all $`aG`$. So the normal closure of $`z^{p^r}HP(G)`$ in $`G/HP(G)`$ is abelian. Therefore $`z^{p^r}^G`$ is (locally nilpotent)-by-abelian, hence Lemma 2.4 implies that $`z^{p^r}HP(G)`$.
The case $`p=2`$ is similar. ∎

Proof of Theorem 1.2. Suppose that $`p`$ is odd and every $`p`$-group of exponent dividing $`p^r`$ in $`𝒩_c^{(2)}`$ is locally nilpotent. Let $`G`$ be a $`p`$-group in $`𝒩_c^{(2)}`$, then by Lemma 3.2(a), $`G^{p^r}`$ is locally nilpotent. By assumption, $`G/G^{p^r}`$ is locally nilpotent and since $`G`$ is a $`c`$-Engel group, it follows from Lemma 2.4 that $`G`$ is locally nilpotent. The case $`p=2`$ is similar and the converse is obvious. $`\mathrm{}`$

Now we use this result to some special cases. In fact we prove:

###### Proposition 3.3.

Every $`2`$-group or $`3`$-group in the variety $`𝒩_5^{(2)}`$ is locally nilpotent.

###### Proof.

By Theorem 1.2, we must prove that every $`2`$-group of exponent dividing $`8`$ and every $`3`$-group of exponent dividing $`9`$ in $`𝒩_5^{(2)}`$ is locally nilpotent. Suppose that $`G`$ is a $`2`$-group of exponent $`8`$ in $`𝒩_5^{(2)}`$. Let $`x,yG`$, then $`x,y𝒩_5`$ and of exponent $`8`$. It is easy to see that $`[x^4,y,x^4,y]=1`$. So by Lemma 2.1, $`x^4HP(G)`$ for all $`xG`$. Therefore $`G/HP(G)`$ is of exponent dividing $`4`$. By a famous result of Sanov (see 14.2.4 of ), $`G/HP(G)`$ is locally nilpotent and so by Lemma 2.4, $`G`$ is locally nilpotent.
Now suppose that $`G`$ is a $`3`$-group of exponent $`9`$ in the variety $`𝒩_5^{(2)}`$. Let $`x,yG`$, it is easy to see that $`[x^3,y,x^3,y]=1`$. The rest of the proof is similar to latter case, but we may use this well-known result that every group of exponent $`3`$ is nilpotent (see 12.3.5 and 12.3.6 of ) ∎