Distribution of partition function zeros of the \pm J model on the Bethe lattice

# Distribution of partition function zeros of the ±J model on the Bethe lattice

Yoshiki Matsuda, Markus Müller, Hidetoshi Nishimori, Tomoyuki Obuchi and Antonello Scardicchio Department of Physics, Tokyo Institute of Technology, Oh-okayama, Meguro-ku, Tokyo 152-8551, Japan International Center for Theoretical Physics, Strada Costiera 11, 34014 Trieste, Italy INFN, Sezione di Trieste, Strada Costiera 11, 34014, Trieste, Italy.
###### Abstract

The distribution of partition function zeros is studied for the model of spin glasses on the Bethe lattice. We find a relation between the distribution of complex cavity fields and the density of zeros, which enables us to obtain the density of zeros for the infinite system size by using the cavity method. The phase boundaries thus derived from the location of the zeros are consistent with the results of direct analytical calculations. This is the first example in which the spin-glass transition is related to the distribution of zeros directly in the thermodynamical limit. We clarify how the spin-glass transition is characterized by the zeros of the partition function. It is also shown that in the spin-glass phase a continuous distribution of singularities touches the axes of real field and temperature.

###### pacs:
05.50.+q, 75.50.Lk

## 1 Introduction

The physical quantities of a spin system are completely determined by the location of the zeros of partition function on the complex parameter planes, such as the temperature and external field. This means that all the singularities of the system should be identified with zeros, and a phase transition occurs when a system parameter is driven across an accumulation point of zeros.

For ferromagnetic systems, Lee and Yang proved that the zeros on the complex plane (Lee-Yang zeros) lie on the imaginary axis. This result is the famous circle theorem . A ferromagnetic transition occurs when the zeros touch the real-field axis at , and zeros split the complex-field plane in two regions, and , where denotes the real part. The density of zeros at the origin is directly proportional to the spontaneous magnetization in the ferromagnetic phase. As for the temperature plane (Fisher zeros), Fisher found that zeros of a pure ferromagnetic system on the square lattice in the absence of an external field lie on the unit circle in the complex plane of , where is the coupling constant and the inverse temperature .

The Lee-Yang circle theorem is valid for random ferromagnets (where all the interactions are greater than or equal to ) and all the Lee-Yang zeros lie on the imaginary-field axis. Therefore the problem is reduced to finding the density of zeros on the imaginary-field axis , which can be calculated from the analytic continuation of the magnetization as a function of the real positive field , as () . For diluted ferromagnets, several works have explored the connection between the Griffiths singularity  and the density of zeros both in solvable models [4, 5] and experiments . These studies suggest that the density of zeros has an essential singularity on the imaginary field axis at the origin as , where is a positive constant. The essential singularity at the origin, in Griffiths’ interpretation , is caused by the presence of large clusters due to random fluctuations of the interactions.

For spin glasses, the situation is more complicated. The locations of zeros are not restricted to the imaginary-field axis. Since it is difficult to treat analytically the partition functions of these systems, zeros have been studied mainly by numerical evaluation of partition functions of finite-size systems on the complex external field [7, 8] and temperature [9, 10, 11] planes. The zeros were found to be distributed generally on two-dimensional areas, not only along a line as in ferromagnets. Recently a detailed study of the density of zeros on the imaginary-field axis suggested the existence of the Griffiths singularity also in spin-glass models in the paramagnetic phase .

It is therefore highly desirable to study the distribution of zeros of spin glasses in the limit of infinite system size since the complex behaviour of spin glasses should manifest itself in non-trivial distributions of zeros. We therefore investigate in this paper the relation between the distribution of zeros and the spin-glass transition in the infinite system size for the spin-glass model on the Bethe lattice. We here define the Bethe lattice as the interior part of the infinitely large Cayley tree. This definition enables us to investigate the zeros of infinite size systems without approximations (even within the spin-glass phase) using the cavity method , which is perfectly suited for our purpose of studying the distribution of zeros of spin glasses. From the distribution of zeros we successfully detect the phase boundary of the spin-glass phase, defined as the line where the spin-glass susceptibility diverges, in accordance with previous studies. We also show that the system is singular everywhere as a function of the real-valued temperature or field in the spin-glass phase.

The outline of this paper is as follows. We first introduce the models and the relationship between the density of zeros and the distribution of cavity field in the next section. Then we explain the actual numerical procedures in section 3. The results are described in section 4, where distributions of zeros on the complex field and temperature planes are investigated and the phase diagrams are determined. The existence of the Griffiths singularity is also discussed. Finally we conclude this paper in section 5.

## 2 Formulation

In this section, the main ideas of this paper are presented. It is shown that the zeros of the partition function for the model on the Bethe lattice are efficiently evaluated using the cavity method.

### 2.1 Cayley tree and the cavity method

Let us consider the model of spin glasses on a Cayley tree in a uniform magnetic field . The Cayley tree is a cycle-free graph where each site is connected to neighbours (figure 1). The Hamiltonian is

 H=−∑⟨i,j⟩JijSiSj−H∑iSi, (1)

where is a nearest neighbour pair, and the value of interaction distributes as

 P(Jij)=pδ(Jij−J)+(1−p)δ(Jij+J) (2)

with and .

On the Cayley tree we can efficiently calculate the partition function by iteratively summing over spin variables layer by layer from the outer boundary. The partial partition function , which represents the partition function in the absence of deeper layers than the site , is updated as

 (3)

where labels spins on the layer previous to the current site (figure 1). The independent effects of spins on the previous layer have been passed to the present layer in terms of the cavity field and cavity biases , the definitions of which are given by

 hi = H+1βc−1∑j=1tanh−1(tanh(βJij)tanh(βhj)) (4) = H+c−1∑j=1uj(Jij,hj). (5)

We hereafter assume that the function takes the principal branch, which restricts the value of the imaginary part of the bias to a range . From now on we assume the imaginary part of all the fields defined modulo and the principal branch of the is considered. The reader should keep this in mind, as we will not burden the notation with this condition explicitly.

At the final step of this procedure, we consider the contribution from the central site of the Cayley tree, and the partition function of the whole system is obtained as

 Z(c)=2cosh(βh(c))c∏j=1[cosh(βJij)cosh(βuj(Jij,hj))Zj], (6)

and

 h(c)=H+c∑j=1uj(Jij,hj), (7)

where the superscript represents the central site.

To consider the typical behaviour of the Cayley tree (where we assume uncorrelated, typical boundary conditions, to be specified explicitly below), we perform the average over the quenched randomness, which introduces the distribution of the cavity field at the th layer . The update rule of is given by

 Pl+1H,β(h)=∫[δ(h−H−c−1∑j=1uj(hj))]Jc−1∏j=1PlH,β(hj)dhj, (8)

where denotes the average over the random interactions . In the thermodynamic limit, and if the distributions converge as to , the limiting distribution satisfies the following self-consistent equation

 PH,β(h)=∫[δ(h−H−c−1∑j=1uj(hj))]Jc−1∏j=1PH,β(hj)dhj. (9)

The existence of a unique solution of this equation is shown rigorously in  for real values of the cavity field. Although this proof is not directly applicable to our case of complex cavity field, we observed numerically that equation (9) has a solution for . The special case of pure ferromagnet () in a purely imaginary external field is discussed in detail in Appendix A.

### 2.2 Definition of the Bethe lattice

The Cayley tree is peculiar because the number of surface sites is comparable to the total number of sites and hence the contribution from the surface cannot be neglected. This property is inconvenient for studying the bulk of the system. Thus the Bethe lattice is often considered instead. There are two different definitions of the Bethe lattice, and we have to clearly distinguish them .

1. The interior part of the Cayley tree: A lattice consisting of the interior part of an infinite-size Cayley tree with uncorrelated boundary conditions. Alternatively, we can define this lattice as a finite Cayley tree, the boundary conditions of which are given by the convergent cavity field distribution of the infinite-size Cayley tree.

2. The regular random graph (RRG): A randomly generated graph under the constraint of a fixed connectivity . Since there exist many cycles, we cannot exactly treat the finite system size. In the limit under appropriate conditions (and outside the spin-glass phase), however, the contribution coming from the loops is expected to become negligible and the problem can be solved exactly by the cavity method.

Both definitions have advantages and disadvantages, and typically the phenomena that occur in one set-up have a correspondence in the other. The model (i) and its phase diagram have been discussed in detail in , while for (ii) and its relations with the replica theory we refer the reader to reference .

In this paper, we refer to model (i) (with typical boundary conditions which are not correlated among each other or with the couplings in the interior) and call it the Bethe lattice. Although definition (ii) is useful in that it has no ambiguity about the boundary conditions, it is incompatible with the formulation of the partition function zeros that uses equation (3). On the other hand, according to definition (i), we can easily define a finite-size system of the Bethe lattice and look at its zeros in the formalism of the cavity method.

### 2.3 The partition function zeros of the Bethe lattice

Let us first consider the zeros of the partition function with respect to the external field , the Lee-Yang zeros . It is known that these zeros are related to the magnetization of the system. In general, the partition function of an Ising system in an external field is expressed as a two-variable polynomial,

 Z(H,β)=eNβHeNbβJN∑M=0Nb∑E=0Ω(E,M)e−2EβJe−2MβH, (10)

where is the number of sites and is the number of interactions. We can write the above equation with fixed temperature as

 Z(H)=ξeNβHN∏i(e−2βH−e−2βHi), (11)

where is a constant (to be ignored hereafter). Equation (11) means that the partition function is a polynomial of degree in the fugacity , and there are roots on the complex plane. Taking the logarithm and dividing by , we find

 −βf(H)=βH+∫∫d2H′g(H′)log(e−2βH−e−2βH′). (12)

where is defined as the density of zeros on the complex plane. The complex magnetization is expressed as

 m(H)=1+∫∫d2H′g(H′)2e2β(H−H′)−1. (13)

We can express the density of zeros in terms of the magnetization by an infinitesimal closed line integral,

 g(H)=β2πilimr→01πr2∮|H−H′|=rm(H′)dH′ (14)

as is easily seen from the representation (13). Indeed the line integral picks up all poles within the circle, and all poles have the residue .

For the Bethe lattice, due to the uniformity of the system, the disorder averaged complex magnetization is given by

 (15)

where the distribution of the central complex field is calculated from the convergent distribution as

 P(c)H,β(h(c))=∫∫[δ(h(c)−H−c∑j=1uj(hj))]Jc∏j=1PH,β(hj)d2hj. (16)

Inserting the equation (15) into equation (14), we obtain

 g(H) = β2πilimr→01πr2∮|H−H′|=rdH′∫∫d2h(c)P(c)H′(h(c))tanh(βh(c)(H′)) (17) = limr→01πr2∫∫|H−H′|≤rd2H′∫∫d2h(c)P(c)H′(h(c))δ(βh(c)(H′)−πi2) = limr→01πr2∫∫|H−H∗|≤rd2H′P(c)H′(πi/2β) = P(c)H(πi/2β).

In the second line, we used the residue theorem under the condition that the radius is sufficiently small. This result connects the density of zeros at with the density of cavity fields at fixed external field . If we iterate the cavity field population numerically for a given pair of complex values , the distribution function of the cavity field yields the density of zeros in the -plane through equation (17).

The relation (17) can also be interpreted as follows. The whole partition function of the Bethe lattice is formally the same as that of the Cayley tree (6), the difference only being whether or not we use the limiting distribution . This implies that the equation of zeros becomes identical to111Note that equation (6) may seem to diverge when the factor becomes , but this is not the case. The reason is that the condition is always accompanied by and yields a finite value.

 2cosh(βh(c))=0⇒βh(c)=π2i. (18)

Equations (17) and (18) show that the complex support of the zeros of the partition function can be obtained from the knowledge of the distribution of the field acting on the central spin. This exact relation represents an important advantage of considering the Bethe lattice as the interior part of the infinite-size Cayley tree. It is one of the main results of the present paper.

We now turn our attention to the density of zeros as a function of the complex temperature. The condition (18) for the vanishing of the partition function is still valid for complex values of temperatures. In order to find the correct density one should follow steps analogous to those which lead from (11) to (17), by treating instead of . We have however not followed this route here. Instead we contented ourselves with detecting the support of the zeros, i.e., the location on the temperature plane where is non-zero. This can be performed by using the fact that the value of the density on the field plane at a complex temperature , , should be proportional to the density on the temperature plane in a field , . This means that the equality holds, where is a normalization factor (to be dropped hereafter) 222We can accept this equality by considering fact that the zeros are located in the four-dimensional complex temperature-field space, and the two-dimensional distributions with real and with real are just two-dimensional cross sections of the same distribution function in four dimensions.. This equality is sufficient to understand the phase diagram in terms of the locations of zeros, while a quantitative knowledge of the density of zeros on the temperature plane cannot be obtained due to the nontrivial normalization factor . Anyway we do not discuss quantitatively the value of density on the temperature plane in the present work.

## 3 Numerical procedures

Hereafter, we discuss the density of zeros on the plane since always appears in this combination. A hat on a variable, like , will mean the same quantity multiplied by . The location of the zeros on the complex temperature plane will also be analyzed.

For the spin-glass model, the zeros on the complex field plane are not restricted to the imaginary field axis but a finite fraction spreads into the complex plane. The remaining fraction still lies on the imaginary field axis. The total density of zeros on the field plane thus naturally splits into two parts,

 g(^H)=δ(R^H)g1(θ)+g2(^H), (19)

where () is the imaginary part of . The one-dimensional measure is restricted to the imaginary axis and the two-dimensional measure is continuously distributed on the complex plane.

We base our analysis mostly on equation (17). However, it is not sufficient for our purpose because it only yields the two-dimensional part . In fact, the cavity field distribution always spreads in the two-dimensional -plane due to the initial conditions that we choose for the Bethe spin glass, as explained later. It is not possible to obtain the one-dimensional part together with the two-dimensional from equation (17333For non-frustrated systems, we can obtain from the relation (17) (see Appendix A).. However we can calculate by a different method using the fact that the one-dimensional density on the imaginary field axis corresponds to the jump of the real part of the complex magnetization as crosses the imaginary axis. In other words, . The analogy with electrostatics developed in reference  allows us to apply this formula even when . Therefore we can calculate from equation (15).

We now obtain the distribution of for complex field and temperature . Let us explain the algorithms to obtain and estimate the density of zeros. We write the recursion relation (5) as

 (20)

where the interaction connects the current site with the next site , and labels the previous sites. The central cavity field is

 ^h(c)=^H+2βc∑j=1uj. (21)

Our numerical solution for the probability densities and is based on the method of population dynamics. We represent the distribution functions and in terms of a large number of variables and , whose distributions are supposed to follow the respective probability distributions. The elements of these sets are updated by randomly choosing and and following equations (20) and (21). The initial condition of the cavity bias, which correspond to the boundary condition of the Cayley tree, is taken to be . These initial fields are uncorrelated random variables, which introduce frustration into the system. Note that in general the value of has both real and imaginary parts because the value of the next generation is calculated from equations (20) and (21) with complex and . The distribution of is calculated from the Monte Carlo (dynamical) average of the population in the two dimensional complex plane, after the dynamics has converged to a limiting distribution.

The actual steps of the algorithm are as follows;

1. Set the initial probability distribution of as .

2. Update the sets of and with the recursion relations (20) and (21) by randomly choosing and out of the set until converges.

3. Estimate the density of zeros by

 g1(θ,T) = 12πR m(^H=iθ+0+) (22) = 12πRlim^HR→0+∫∫d2^h(c)P(c)^H(^h(c))tanh(^h(c)/2)∣∣∣^H=iθ+^HR g2(^H,T) = P(c)^H(^h(c)=πi) (23)

for given and . The one-dimensional density represents a density of zeros under a pure-imaginary field while the two-dimensional density is defined on the whole complex planes.

For simplicity only the Bethe lattice with connectivity has been studied. We have chosen representative points in the population dynamics and have performed at least cavity iterations until the population converges. The data were collected from the average of additional iterations after the initial (or more) iterations. These conditions have been used throughout this work.

## 4 Distribution of zeros for the Bethe spin glass

### 4.1 Zeros on the complex field plane

First, we show the density of zeros on the complex plane for real temperature. Figure 2 is the distributions of zeros on the complex plane with probability at temperature (left) and (right). The complex plane has been split into cells by dividing the real axis from to with an increment and the imaginary axis from to with an increment of . The density outside this range is omitted, since this is sufficient to see zeros near the real axis which are essential for critical phenomena. Both and are plotted in the same figure and coloured in a logarithmic scale; a black dot shows a very high density. Figure 2: Distribution of zeros on the complex ^H plane with p=0.5 at T=1.43 (paramagnet, left panel) and T=0.5 (spin-glass, right panel). Densities are coloured in a logarithmic scale. On the right panel, zeros with finite real part touch the real axis, whereas g1 on the imaginary axis vanishes to the numerical precision.

The left panel of figure 2 () is for the paramagnetic phase. Neither nor has a finite value on the real axis and there is thus no phase transition as a function of the real field. The right panel of figure 2 is in the spin-glass phase (), where zeros reach the real axis at and away from the origin. Thus, there is a phase transition at some on the Bethe lattice, where is the point where the density changes its value from zero to non-zero along the real axis. Zeros touch the real axis all the interval between and . This suggests that the free energy is non-analytic as a function of below . The one-dimensional density vanishes to the numerical precision for both temperatures in figure 2. The two-dimensional part has finite values on the imaginary axis. It is likely that the edge of (the point where on the imaginary axis and is closest to the origin) touches the origin as the temperature is decreased to . At lower the density of the zeros approaches the real axis also away from the imaginary axis, up to a point on the real axis. This defines a line of spin-glass transitions in the - plane, as discussed below.

In the ferromagnetic phase (the right part of figure 3), on the other hand, only touches the origin but does not. In marked contrast to the left panel of figure 2, the one-dimensional density has a finite value on the imaginary axis away from the origin on the left panel of figure 3. This is an important difference of ferromagnetic and spin-glass transitions seen from the distribution of Lee-Yang zeros. More details on the transition points will be discussed in later sections. Figure 3: Distribution of zeros on the complex ^H plane with p=0.9 at T=1.5 (left: para) and T=0.5 (right: ferro). Only the one-dimensional zeros on the imaginary axis, g1 (thin lines), touch the origin at low temperature.

### 4.2 Zeros on the complex temperature plane

Using complex values of the temperature in equation (20), we can also obtain information on the support of Fisher zeros, using equation (23). This information is sufficient for our purpose of detecting phase transitions from the point of view of zeros.

Figure 4 shows the distribution of zeros on the temperature plane at with (left) and (right). The (two-dimensional) zeros also touch the real temperature axis below the spin-glass transition temperature. This confirms that the spin-glass transition differs from the ordinary phase transition where zeros touch the real axis only at the transition temperature. This may be interpreted as the system staying critical at all temperatures below the transition temperature. It may also be taken as a signature of temperature chaos, i.e. the instability of the randomly frozen spin configurations in the spin-glass phase with respect to arbitrarily small temperature changes. In accordance with the right panel of figure 2, the right panel of figure 4 shows that the spin-glass phase persists under a weak field (). The critical temperatures appearing in figure 4 are consistent with the phase diagram shown in the next subsection. Figure 4: Distribution of zeros on the complex T plane with p=0.5 and H=0 (left) and H=0.5 (right). Zeros touch the real axis below the critical temperature TSG≃1.13 (left) and 0.5 (right). The apparent absence of zeros near the origin may be due to numerical rounding errors of tanhβ.

The distributions at close to 1 are also interesting. Figure 5 shows the density at with (left) and (right). The transition temperature between paramagnetic and ferromagnetic phases is known to be . It seems that the zeros distribute on an almost one-dimensional curve near which most likely touches the real axis at . Since our method of equation (23) calculates the two-dimensional density of zeros, a one-dimensional density is hard to identify with a high precision. On the right panel of figure 5, a weak field is added, to test for a spontaneous magnetization: zeros have finite densities along the real axis below , since is finite under a pure imaginary field in the ferromagnetic phase.

On the other hand, the two-dimensional distribution touches the real axis again in the low temperature region. Below , the zeros approach the real axis on the whole interval , similarly as in the spin-glass phase shown in figure 4. Indeed, this second critical temperature corresponds to the spin-glass transition, as shown below. Figure 5: Distribution of zeros on the complex T plane for p=0.9 and ^H=0 (left) and ^H=10−4i (right). The zeros approach the real axis at T=Tc, where Tc≃1.36 is the ferromagnetic critical temperature. In addition, zeros reach the real axis also in the low temperature region. The second critical temperature corresponds to the spin-glass transition TSG shown in figure 6. On the right panel, a weak pure-imaginary field immediately makes apparent the presence of ferromagnetic order below a critical temperature Tc: a one-dimensional distribution of zeros lies on the real axis in the right panel.

### 4.3 Phase diagram

Based on the above observations, we investigate the phase diagram on the - and - planes with real and . Graphical representations shown so far suggest that there are two types of transitions where one- and two-dimensional distributions of zeros reach the real axis separately. In order to determine the two transition points, we added a very small imaginary field . Since zeros on the imaginary axis in the complex plane reach the real axis before those away from the imaginary axis as the temperature is decreased, we can restrict to the imaginary axis and ask for the highest temperature for which a non-zero density, and exists at the origin of the imaginary axis (). In contrast, it is difficult to determine the transition points caused by one-dimensional zeros on the complex plane in the absence of a field as shown in the left panel of figure 5. The transition temperature derived from the first method of the complex is close to the point in the left panel of figure 5 where the almost one-dimensional distribution is likely to touch the real axis. Therefore we mainly used the method of complex field to identify the phase boundary (the left panel of figure 6). Figure 6: Phase diagrams on the p-T plane from the density of zeros. The solid lines are analytical results as given as Tc=1/tanh−1[1/(4p−2)] (para-ferro boundary) and TSG=1/tanh−1[1/√2] (para-spin glass boundary) [15, 16, 17, 18]. Left: Circles and triangles are determined by the two-dimensional density of zeros g2 on the complex field plane and temperature plane, respectively, which corresponds to spin-glass critical temperatures. Squares are phase boundaries calculated by the one-dimensional density of zeros g1, which indicates the onset of the ferromagnetic order. Immediately below triangles there is a phase which is both magnetized and glassy. The dotted line denotes the Nishimori line , on which the multicritical point (where three phases merge) is located. Right: Blow up of the field-induced critical temperatures on the left panel (marked in circles and triangles) between p=0.84 and p=0.93.

Our result agrees well with the exact phase boundary [15, 16, 17, 18], (between para and ferro phases) for and (between para and spin-glass phases) for . It is expected that the one-dimensional distribution determines the ferromagnetic temperature and the two-dimensional distribution determines the spin-glass transition temperature . We discuss this hypothesis in detail below.

#### 4.3.1 Behaviour of the two-dimensional distribution of zeros

Below the line drawn in circles in figure 6, the ordered phase is expected to be stable under an external field as shown in figure 2. This line is also determined from the zeros on the temperature plane. We set the temperature as with , where is real, and we assume the highest having a non-zero density as the critical temperature (triangles in figure 6). The two results agree well with each other.

The boundary between the ferromagnetic and spin-glass phases is harder to determine. As one sees in figure 6 (the left panel), the temperature where the one-dimensional distribution ceases to touch the real axis is lower than the temperature where the two-dimensional distribution starts touching the real axis (for biases to ). This indicates that there is a phase where the system is still spontaneously magnetized, but also glassy (see also ). It is shown in Appendix B that the spin-glass susceptibility diverges along a line which coincides with the touching points of , marked in circles and triangles in figure 6 . This implies that our method using the approach of the two-dimensional distribution to the real axis correctly detects the onset of the spin-glass phase. Indeed the proximity of zeros to the real axis at implies the divergence of some (higher order) susceptibility.

We next show the phase diagram on the - plane with fixed (figure 7). We added a complex external filed as and we assumed that the critical value of is the maximum value of having a non-zero (two-dimensional) density . The top panels of figure 7 are for the temperature dependence of the critical line for various . Figure 7: Spin-glass transition lines in the T-H plane (corresponding to the AT line in the SK model). The bottom panel is for p=0.50. The boundary at p=0.5 rises from T=Tc as Hc∝(Tc−T)1.5. The exponent 1.5 agrees with that of the AT line for the SK model. pc=(2+√2)/4≃0.854 is for the multicritical point.

For all biases , smoothly rises from at for , which is for the multicritical point on the Nishimori line  where the three phases merge. At low temperature, it seems that approaches a finite value at . The bottom panel is for , and the points well fit to . The exponent is the same as in the AT line  of the Sherrington-Kirkpatrick model . The existence of the AT line on the Bethe lattice (regular random graph) was predicted in reference  where the AT line should behave like in the vicinity of , and was checked by the population method in references [22, 23].

The agreement of the exponent with the behaviour of the AT line might be interpreted as an indication that the phase below the phase boundary is the spin-glass phase with replica symmetry breaking (RSB) [24, 25, 26]. Indeed, it is suggested that the spin-glass phase in the regular random graph has the full RSB [27, 28]. However, we should be careful because we have not directly seen the instability of a replica-symmetric solution in the Bethe lattice of our definition. The analysis of reference  shows that in the spin-glass phase a system of two replicas on the Cayley tree exhibits a diverging susceptibility with respect to an infinitesimal repulsion between the replicas. However, at the same time the cavity field distribution remains essentially identical to the replica symmetric approximation for the regular random graph; model (ii). Further studies of this phase from different viewpoints may be necessary to fully clarify its nature.

### 4.4 Griffiths singularity

The density of zeros on the imaginary field axis is closely related to the Griffiths singularity in the diluted ferromagnet . The same is expected in spin glasses . The Griffiths singularity is expected to manifest itself, if it is present, in the form of an essential singularity of the density of zeros upon approaching the origin along the imaginary field axis [4, 5, 6, 8]. If such a tail is present (which is very difficult to detect numerically), its touching of the real axis would indicate the onset of a Griffith phase, but not the phase transition to a ferromagnetically ordered or glassy state. In the presence of a Griffith phase, our criterion for the detection of phase transitions should thus strictly speaking be refined to the condition that a ”substantial density of zeros” (which grows at least like a power law with away from the real axis) touches the real axis.

The above discussion suggests to study in more detail the form of the one-dimensional density of zeros on the imaginary axis as an indicator of possible Griffiths singularity for the Bethe spin glass. Figure 8: One- and two-dimensional density of zeros on the imaginary field axis for high p at T=1.5. The top panels are for the ferromagnetic phase and the bottom for the paramagnetic phase. The contribution of the one-dimensional density g1 decreases with decreasing p whereas g2 increases.

For the Bethe spin glass our numerics does not show detectable signs of a Griffiths singularity above the spin-glass phase, our estimated being zero within numerical precision in the paramagnetic phase above the spin-glass phase, as shown in figure 2. Above the ferromagnetic phase, it seems that is finite above some , while there seems to be a jump between and .

Figure 8 shows the one- and two-dimensional densities of zeros on the imaginary axis at fixed with various . The top panels correspond to biases which are in the ferromagnetic phase while the bottom panels correspond to lower biases, in the paramagnetic phase. At high , the part is dominant because the line integrals

 W1=∫π−πg1(θ)dθ (24)

for and are nearly equal to , where with should hold due to the normalization. On the other hand, rapidly decreases with decreasing , and vanishes for large . In the left-bottom panel of figure 8, is seen to decrease when becomes finite around . This behaviour suggests that one-dimensional zeros are buried under large cloud of two-dimensional zeros far away from the real axis. Finally, for the one-dimensional seems to vanish altogether above the spin-glass phase as shown in figure 2.

It is difficult to draw a definite conclusion about Griffiths singularities due to the limit to detect very small but non-vanishing values of near the origin. However, our numerics points toward the possible absence of a Griffiths phase in the spin glass on the Bethe lattice, contrary to what happens for finite-dimensional systems and diluted ferromagnets. A more detailed study will be required to settle this question.

## 5 Conclusion

We have investigated the distribution of the partition function zeros for the spin-glass model on the Bethe lattice. An important relation to connect the density of zeros and the density of cavity fields for complex field and temperature is found, which enables us to treat the zeros in the limit of infinite system size. The densities are split into one- and two-dimensional parts, where the one-dimensional density is defined on the imaginary axis of the complex field plane and the two-dimensional density spreads over the complex plane. We investigated the phase diagram of the Bethe spin glass by estimating the transition points where the one- or two-dimensional densities begin to have finite values in the immediate vicinity of the real axis. Our results agree well with the analytically exact phase diagram. The one-dimensional density determines the boundary between the paramagnetic and ferromagnetic phases since the one-dimensional density is directly related to the spontaneous magnetization. The two-dimensional density determines the boundary of the spin-glass phase on the - and - planes. This phase boundary corresponds to the line where the spin-glass susceptibility diverges as evaluated in Appendix B. Below the spin-glass transition point, the two-dimensional density of zeros continuously touches the real field and temperature axes, which may be related to chaotic behaviour of the system as a function of the field and temperature in the spin-glass phase. We have shown, by looking at the locations of the partition function zeros, that the system has a non-analytic free energy for all below the spin-glass temperature .

We have not observed any evidence for the existence of a Griffiths phase in our numerics. If such a phase is indeed absent, it implies that the Bethe lattice behaves distinctly from finite-dimensional spin glasses. Further careful studies are required to understand the origin of differences between Bethe and Bravais lattices.

This work was partially supported by CREST, JST and by the Grant-in-Aid for Scientific Research on the Priority Area ‘Deepening and Expansion of Statistical Mechanical Informatics’ by the Ministry of Education, Culture, Sports, Science and Technology of Japan. YM and TO are grateful for the financial support provided through the Japan Society for the Promotion of Science (JSPS) Research Fellowship for Young Scientists Program.

## Appendix Appendix A Zeros of the pure ferromagnetic system

In order to check our method itself and the precision of the numerical analyses, we estimated the density of zeros for the pure ferromagnetic system. The zeros of a ferromagnetic system is located only on the imaginary axis, and thus . In this case the density of zeros on the imaginary axis, , can be found by both methods which we used to calculate and in the case (equations (17) and (22)). We here compare the exact density of zeros and two algorithms using equations (17) and (22). Figure 9 is for comparison of numerical calculations and the exact solution. Our algorithms agree perfectly with the exact solution. The details of the analyses are as follows.

1. The exact density of zeros (the solid line in figure 9): The density of zeros on the imaginary field axis is obtained as the analytic continuation of the real from real positive values of to purely imaginary ,

 2πg1(θ)=Rm(^H:=iθ). (25)

Using the fixed point of equation (5) with and , we rewrite equation (25) as

 2πg1(θ)=Rtanh[12(32^hf(^H,β)+^H)]∣∣∣^H:=iθ. (26)

From equation (5) the fixed point for the pure ferromagnetic system in real field is found as

 ^hf(^H,β)=2logx′, (27)

where is the solutions of the following equation,

 e^Hx3−e^H+2βx2+e2βx−1=0. (28)
2. The relation between density of zeros and cavity fields (squares in figure 9): The algorithm shown as equation (17) connects the density of cavity fields and density of zeros. Since the zeros for the ferromagnetic system lie only on the imaginary field axis, the density is estimated from the distribution of cavity field in a purely-imaginary field. We thus consider the free boundary condition in order to make the cavity field pure imaginary.

For the cavity field in a pure imaginary field , the recursion relations (20) and (21) are rewritten as

 I^ui = 2tan−1[tanhβtan(θ/2+(c−1)I^uj/2)]R^ui=0, (29) I^h(c) = θ+cI^uiR^h(c)=0, (30)

where and and the initial is set as . We directly estimate the one-dimensional density of the population on the imaginary field axis by using the relation,

 g1(θ)=P(c)H(I^hc=π). (31)

It can be shown in full generality that this procedure yields a density which correctly describes the jump of the real part of the magnetization across the imaginary axis. For the spin-glass model (), this algorithm is used only in the estimation of part (equation (23)). Since the boundary condition of the Bethe spin-glass is fixed in order to introduce frustration to the system and the cavity fields are complex, this method is not applicable to the estimation of for the spin-glass model.

3. The real part of complex magnetization (circles in figure 9): Complex magnetization in the vicinity of the imaginary field axis is estimated from

 m(θ)=lim^HR→0+tanh[12(32^h′f(^H,β)+^H)]^H=iθ+^HR, (32)

where is the (numerically) fixed point of the complex cavity field. The one-dimensional density on the imaginary axis is also calculated from the real part of the complex magnetization as . This method is used to calculate for the Bethe spin-glass (equation (22)). Figure 9: The density of zeros on the imaginary axis for the pure ferromagnetic system with c=3. The numerical estimations are consistent with the exact calculation of equation (26). The edge of zeros corresponds to θ0 in the text.

Note that the imaginary cavity field (equation (29)) does not converge for the pure ferromagnetic system when the density of zeros has a positive value for a given (the left panel of figure 10). However we take the average over a large number of the iterating steps to calculate the density of cavity field, which naturally performs the sampling of bulk properties of the Cayley tree. If the return map of the recursion relation (5) has a fixed point, on the other hand, the cavity field population is delta-peaked at that fixed point. One finds that the density of zeros vanishes under this condition (see the right panel of figure 10). Figure 10: The return map of equation (21) for the pure ferromagnetic system in pure imaginary field. Left: If the curve y=4tan−1[tanhβtan((θ+I^h)/2)] does not intersect y=I^h, the imaginary cavity field does not converge; the density of zeros has a finite value under this condition. Right: If the two equations intersect, the imaginary cavity field recursion has fixed point; thus the density has no value.

The critical value of the imaginary field is given by the condition that the two curves in figure 10 touch:

 2(c−1)tan−1[tanhβtan(θ0/2+^h/2)]=^h (33)

and

 dd^h[2(c−1)tan−1[tanhβtan(θ0/2+^h/2)]]=1. (34)

We find that the critical imaginary field is determined as

 θ0(β)=2cos−1[√(c−1−tanhβ)sinhβcoshβ]−2(c−1)tan−1[√(1−(c−1)tanhβ)tanhβc−1−tanhβ]. (35)

The critical temperature is given from this equation at as . This result is of course consistent with the well-known critical temperature. Figure 11 shows the temperature dependence of for , , . Figure 11: The temperature dependence of θ0 for connectivities c=2, 3, ⋯ 9 from left to right. At θ0=0, the temperature is given by T=1/tanh−1[1/(c−1)].

## Appendix Appendix B The critical condition based on the spin-glass susceptibility

It is quite possible that the divergence of the spin-glass susceptibility , which in the SK model is identified with the AT instability [14, 29], characterizes the instability signaled by the zeros of the partition function approaching the real axis. Here we show that this is indeed the case for the present Bethe lattice.

The spin-glass susceptibility is defined as

 χSG=1N∑i,j⎡⎣(∂⟨Si⟩∂hj)2⎤⎦J=∑j⎡⎣(∂⟨S0⟩∂hj)2⎤⎦J. (36)

To derive the second identity, we assumed uniformity of the Bethe lattice and selected the central spin as . In a cycle-free graph, an arbitrary pair of sites are connected by a single path. Let us assign site indexes from the origin 0 of the graph to a site of distance along the path as . For a fixed set of couplings and boundary fields, the chain rule of the derivative shows that

 ∂⟨S0⟩∂hG = ∂⟨S0⟩∂h0∂h0∂u0∂u0∂h1⋯∂hG∂uG=∂⟨S0⟩∂h0∂h0∂u0G∏g=1∂ug−1∂hg∂hg∂ug (37) = ∂⟨S0⟩∂u0G∏g=1∂ug−1∂ug, (38)

as is the cavity field on site and depends linearly on as , where represents the sum of the cavity biases from the other branches that flow into site . In the limit , the relevant factor in equation (38) is only . The spin-glass susceptibility is hence evaluated as

 (39)

where the factor denotes the number of sites of distance from the central site . The divergence condition of is given by

 log(c−1)+limG→∞1Glog⎛⎜⎝⎡⎣G∏g=1(∂ug−1∂ug)2⎤⎦J⎞⎟⎠=0. (40)

This yields the condition of the spin-glass transition of the Bethe lattice and RRG.

In order to estimate the divergence points of the spin-glass susceptibility, we numerically implement the calculation of the factor . This factor can also be written as and this latter form is more tractable at finite temperatures.

The numerical evaluation of the factor is straightforward,

 ∂u0∂uG ≈ u0(uG+ΔuG)−u0(uG)ΔuG. (41)

The procedure to evaluate this equation is as follows [28, 30, 31]. We arrange two replicas of an identical population expressing the convergent (real) cavity-bias distribution , which is related to the convergent cavity-field distribution given in equation (9) as

 πH,β(u)=∫dhPH,β(h)[δ(u−1βtanh−1{tanhβJtanhβh})]J. (42)

In addition, we introduce a uniform perturbation () into only one of the two replicas and then observe the square average of the variation, , after a certain number of the cavity updates.

In particular, we update two populations by iterations with the same set of . A critical line of the divergence of the spin-glass susceptibility is determined by whether the square average is numerically zero or much larger than the perturbation. The result is shown in figure 12 where the spin-glass susceptibility diverges below this line. At zero temperature, the critical probability is calculated as in references [27, 32], which is reproduced by our numerical calculation at zero temperature as .

This result agrees well with the phase boundary drawn in figure 6. Thus it is suggested that our phase boundary estimated by the two-dimensional zeros corresponds to the spin-glass transitions. Figure 12: The divergence points of the spin-glass susceptibility calculated by the instability of the real cavity field distribution. The critical probability at zero temperature is suggested as pc=0.91665(5). The circles denote the phase boundary estimated by the two-dimensional distribution of zeros (figure 6). They agree well with each other.

## References

•  Lee T D and Yang C N 1952 Phys. Rev. 87 410
•  Fisher M E 1964 The Nature of Critical Points (Lectures in Theoretical Physics vol 7C) (the University of Colorado Press, Boulder)
•  Griffiths R B 1969 Phys. Rev. Lett. 23 17
•  Bray A J and Huifang D 1989 Phys. Rev. B 40 6980
•  Laumann C, Scardicchio A and Sondhi S L 2008 Phys. Rev. E 77 61139
•  Chan P Y, Goldenfeld N and Salamon M 2006 Phys. Rev. Lett. 97 137201
•  Ozeki Y and Nishimori H 1988 J. Phys. Soc. Jpn. 57 1087
•  Matsuda Y, Nishimori H and Hukushima K 2008 J. Phys. A: Math. Theor. 41 324012
•  Bhanot G and Lacki J 1993 J. Stat. Phys. 71 259
•  Saul L and Kardar M 1993 Phys. Rev. E 48 R3221
•  Damgaard P H and Lacki J 1995 Int. J. Mod. Phys. 6 819
•  Mézard M and Parisi G 2001 Eur. Phys. J. B 20 217
•  Chayes J T, Chayes L, Sethna J and Thouless D J 1986 Commun. Math. Phys. 106 41
•  Obuchi T, Kabashima Y and Nishimori H 2009 J. Phys. A: Math. Theor. 42 5004
•  Thouless D J 1986 Phys. Rev. Lett. 56 1082
•  Carlson J M, Chayes J T, Chayes L, Sethna J and Thouless D J 1990 J. Stat. Phys. 61 987
•  Carlson J M, Chayes J T, Sethna J P and Thouless D J 1990 J. Stat. Phys. 61 1069
•  Kabashima Y 2003 J. Phys. Soc. Jpn. 72 1645
•  Nishimori H 1981 Prog. Theor. Phys. 66 1169
•  de Almeida J R L and Thouless D J 1978 J. Phys. A: Math. Gen. 11 983
•  Sherrington D and Kirkpatrick S 1975 Phys. Rev. Lett. 35 1792
•  Pagnani A, Parisi G and Ratiéville M 2003 Phys. Rev. E 68 046706
•  Jörg T, Katzgraber H G and Krza¸kała F 2008 Phys. Rev. Lett. 100 197202
•  Mézard M, Parisi G and Virasoro M A 1987 Spin Glass Theory and Beyond (Singapore: World Scientific)
•  Fisher K H and Hertz J A 1991 Spin Glasses (Cambridge: Cambridge University Press)
•  Nishimori H 2001 Statistical Physics of Spin Glasses and Information Processing: An Introduction (Oxford University Press)
•  Castellani T, Krzakala F and Ricci-Tersenghi F 2005 Eur. Phys. J. B 47 99
•  Krza¸kała F 2005 Prog. Theor. Phys. Supp. No. 157 77
•  Martin O C, Mézard M and Rivoire O 2005 J. Stat. Mech. 2005 P09006
•  Rivoire O, Biroli G, Martin O C and Mezard M 2003 Eur. Phys. J. B 37 55
•  Zdeborová L 2009 Acta Physica Slovaca 59 169
•  Kwon C and Thouless D J 1988 Phys. Rev. B 37 7649
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters   