Distribution of partition function zeros of the model on the Bethe lattice
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 spinglass transition is related to the distribution of zeros directly in the thermodynamical limit. We clarify how the spinglass transition is characterized by the zeros of the partition function. It is also shown that in the spinglass phase a continuous distribution of singularities touches the axes of real field and temperature.
pacs:
05.50.+q, 75.50.Lk1 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 (LeeYang zeros) lie on the imaginary axis. This result is the famous circle theorem [1]. A ferromagnetic transition occurs when the zeros touch the realfield axis at , and zeros split the complexfield 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 [2].
The LeeYang circle theorem is valid for random ferromagnets (where all the interactions are greater than or equal to ) and all the LeeYang zeros lie on the imaginaryfield axis. Therefore the problem is reduced to finding the density of zeros on the imaginaryfield axis , which can be calculated from the analytic continuation of the magnetization as a function of the real positive field , as () [1]. For diluted ferromagnets, several works have explored the connection between the Griffiths singularity [3] and the density of zeros both in solvable models [4, 5] and experiments [6]. 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 [3], 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 imaginaryfield 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 finitesize systems on the complex external field [7, 8] and temperature [9, 10, 11] planes. The zeros were found to be distributed generally on twodimensional areas, not only along a line as in ferromagnets. Recently a detailed study of the density of zeros on the imaginaryfield axis suggested the existence of the Griffiths singularity also in spinglass models in the paramagnetic phase [8].
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 nontrivial distributions of zeros. We therefore investigate in this paper the relation between the distribution of zeros and the spinglass transition in the infinite system size for the spinglass 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 spinglass phase) using the cavity method [12], 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 spinglass phase, defined as the line where the spinglass susceptibility diverges, in accordance with previous studies. We also show that the system is singular everywhere as a function of the realvalued temperature or field in the spinglass 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 cyclefree graph where each site is connected to neighbours (figure 1). The Hamiltonian is
(1) 
where is a nearest neighbour pair, and the value of interaction distributes as
(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
(4)  
(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
(6) 
and
(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
(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 selfconsistent equation
(9) 
The existence of a unique solution of this equation is shown rigorously in [13] 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 [14].

The interior part of the Cayley tree: A lattice consisting of the interior part of an infinitesize 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 infinitesize Cayley tree.

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 spinglass 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 setup have a correspondence in the other. The model (i) and its phase diagram have been discussed in detail in [13], while for (ii) and its relations with the replica theory we refer the reader to reference [12].
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 finitesize 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 LeeYang zeros [1]. 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 twovariable polynomial,
(10) 
where is the number of sites and is the number of interactions. We can write the above equation with fixed temperature as
(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
(12) 
where is defined as the density of zeros on the complex plane. The complex magnetization is expressed as
(13) 
We can express the density of zeros in terms of the magnetization by an infinitesimal closed line integral,
(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
(16) 
Inserting the equation (15) into equation (14), we obtain
(17)  
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 to^{1}^{1}1Note 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.
(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 infinitesize 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 nonzero. 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) ^{2}^{2}2We can accept this equality by considering fact that the zeros are located in the fourdimensional complex temperaturefield space, and the twodimensional distributions with real and with real are just twodimensional 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 spinglass 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,
(19) 
where () is the imaginary part of . The onedimensional measure is restricted to the imaginary axis and the twodimensional 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 twodimensional part . In fact, the cavity field distribution always spreads in the twodimensional plane due to the initial conditions that we choose for the Bethe spin glass, as explained later. It is not possible to obtain the onedimensional part together with the twodimensional from equation (17) ^{3}^{3}3For nonfrustrated systems, we can obtain from the relation (17) (see Appendix A).. However we can calculate by a different method using the fact that the onedimensional 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, [1]. The analogy with electrostatics developed in reference [1] 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
(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;

Set the initial probability distribution of as .

Estimate the density of zeros by
(22) (23) for given and . The onedimensional density represents a density of zeros under a pureimaginary field while the twodimensional 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.
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 spinglass 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 nonzero along the real axis. Zeros touch the real axis all the interval between and . This suggests that the free energy is nonanalytic as a function of below . The onedimensional density vanishes to the numerical precision for both temperatures in figure 2. The twodimensional 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 spinglass 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 onedimensional 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 spinglass transitions seen from the distribution of LeeYang zeros. More details on the transition points will be discussed in later sections.
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 (twodimensional) zeros also touch the real temperature axis below the spinglass transition temperature. This confirms that the spinglass 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 spinglass 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 spinglass phase persists under a weak field (). The critical temperatures appearing in figure 4 are consistent with the phase diagram shown in the next subsection.
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 onedimensional curve near which most likely touches the real axis at . Since our method of equation (23) calculates the twodimensional density of zeros, a onedimensional 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 twodimensional 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 spinglass phase shown in figure 4. Indeed, this second critical temperature corresponds to the spinglass transition, as shown below.
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 twodimensional 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 nonzero density, and exists at the origin of the imaginary axis (). In contrast, it is difficult to determine the transition points caused by onedimensional 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 onedimensional 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).
Our result agrees well with the exact phase boundary [15, 16, 17, 18], (between para and ferro phases) for and (between para and spinglass phases) for . It is expected that the onedimensional distribution determines the ferromagnetic temperature and the twodimensional distribution determines the spinglass transition temperature . We discuss this hypothesis in detail below.
4.3.1 Behaviour of the twodimensional 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 nonzero density as the critical temperature (triangles in figure 6). The two results agree well with each other.
The boundary between the ferromagnetic and spinglass phases is harder to determine. As one sees in figure 6 (the left panel), the temperature where the onedimensional distribution ceases to touch the real axis is lower than the temperature where the twodimensional 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 [13]). It is shown in Appendix B that the spinglass 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 twodimensional distribution to the real axis correctly detects the onset of the spinglass 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 nonzero (twodimensional) density . The top panels of figure 7 are for the temperature dependence of the critical line for various .
For all biases , smoothly rises from at for , which is for the multicritical point on the Nishimori line [19] 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 [20] of the SherringtonKirkpatrick model [21]. The existence of the AT line on the Bethe lattice (regular random graph) was predicted in reference [15] 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 spinglass phase with replica symmetry breaking (RSB) [24, 25, 26]. Indeed, it is suggested that the spinglass 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 replicasymmetric solution in the Bethe lattice of our definition. The analysis of reference [13] shows that in the spinglass 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 [3]. The same is expected in spin glasses [8]. 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 onedimensional density of zeros on the imaginary axis as an indicator of possible Griffiths singularity for the Bethe spin glass.
For the Bethe spin glass our numerics does not show detectable signs of a Griffiths singularity above the spinglass phase, our estimated being zero within numerical precision in the paramagnetic phase above the spinglass 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 twodimensional 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
(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 leftbottom panel of figure 8, is seen to decrease when becomes finite around . This behaviour suggests that onedimensional zeros are buried under large cloud of twodimensional zeros far away from the real axis. Finally, for the onedimensional seems to vanish altogether above the spinglass 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 nonvanishing 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 finitedimensional 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 spinglass 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 twodimensional parts, where the onedimensional density is defined on the imaginary axis of the complex field plane and the twodimensional 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 twodimensional 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 onedimensional density determines the boundary between the paramagnetic and ferromagnetic phases since the onedimensional density is directly related to the spontaneous magnetization. The twodimensional density determines the boundary of the spinglass phase on the  and  planes. This phase boundary corresponds to the line where the spinglass susceptibility diverges as evaluated in Appendix B. Below the spinglass transition point, the twodimensional 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 spinglass phase. We have shown, by looking at the locations of the partition function zeros, that the system has a nonanalytic free energy for all below the spinglass 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 finitedimensional spin glasses. Further careful studies are required to understand the origin of differences between Bethe and Bravais lattices.
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.

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 ,
(25) Using the fixed point of equation (5) with and , we rewrite equation (25) as
(26) From equation (5) the fixed point for the pure ferromagnetic system in real field is found as
(27) where is the solutions of the following equation,
(28) 
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 purelyimaginary 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
(29) (30) where and and the initial is set as . We directly estimate the onedimensional density of the population on the imaginary field axis by using the relation,
(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 spinglass model (), this algorithm is used only in the estimation of part (equation (23)). Since the boundary condition of the Bethe spinglass 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 spinglass model.

The real part of complex magnetization (circles in figure 9): Complex magnetization in the vicinity of the imaginary field axis is estimated from
(32) where is the (numerically) fixed point of the complex cavity field. The onedimensional 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 spinglass (equation (22)).
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 deltapeaked at that fixed point. One finds that the density of zeros vanishes under this condition (see the right panel of figure 10).
The critical value of the imaginary field is given by the condition that the two curves in figure 10 touch:
(33) 
and
(34) 
We find that the critical imaginary field is determined as
(35) 
The critical temperature is given from this equation at as . This result is of course consistent with the wellknown critical temperature. Figure 11 shows the temperature dependence of for , , .
Appendix Appendix B The critical condition based on the spinglass susceptibility
It is quite possible that the divergence of the spinglass 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 spinglass susceptibility is defined as
(36) 
To derive the second identity, we assumed uniformity of the Bethe lattice and selected the central spin as . In a cyclefree 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
(37)  
(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 spinglass 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
(40) 
This yields the condition of the spinglass transition of the Bethe lattice and RRG.
In order to estimate the divergence points of the spinglass 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,
(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) cavitybias distribution , which is related to the convergent cavityfield distribution given in equation (9) as
(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 spinglass 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 spinglass 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 twodimensional zeros corresponds to the spinglass transitions.
References
References
 [1] Lee T D and Yang C N 1952 Phys. Rev. 87 410
 [2] Fisher M E 1964 The Nature of Critical Points (Lectures in Theoretical Physics vol 7C) (the University of Colorado Press, Boulder)
 [3] Griffiths R B 1969 Phys. Rev. Lett. 23 17
 [4] Bray A J and Huifang D 1989 Phys. Rev. B 40 6980
 [5] Laumann C, Scardicchio A and Sondhi S L 2008 Phys. Rev. E 77 61139
 [6] Chan P Y, Goldenfeld N and Salamon M 2006 Phys. Rev. Lett. 97 137201
 [7] Ozeki Y and Nishimori H 1988 J. Phys. Soc. Jpn. 57 1087
 [8] Matsuda Y, Nishimori H and Hukushima K 2008 J. Phys. A: Math. Theor. 41 324012
 [9] Bhanot G and Lacki J 1993 J. Stat. Phys. 71 259
 [10] Saul L and Kardar M 1993 Phys. Rev. E 48 R3221
 [11] Damgaard P H and Lacki J 1995 Int. J. Mod. Phys. 6 819
 [12] Mézard M and Parisi G 2001 Eur. Phys. J. B 20 217
 [13] Chayes J T, Chayes L, Sethna J and Thouless D J 1986 Commun. Math. Phys. 106 41
 [14] Obuchi T, Kabashima Y and Nishimori H 2009 J. Phys. A: Math. Theor. 42 5004
 [15] Thouless D J 1986 Phys. Rev. Lett. 56 1082
 [16] Carlson J M, Chayes J T, Chayes L, Sethna J and Thouless D J 1990 J. Stat. Phys. 61 987
 [17] Carlson J M, Chayes J T, Sethna J P and Thouless D J 1990 J. Stat. Phys. 61 1069
 [18] Kabashima Y 2003 J. Phys. Soc. Jpn. 72 1645
 [19] Nishimori H 1981 Prog. Theor. Phys. 66 1169
 [20] de Almeida J R L and Thouless D J 1978 J. Phys. A: Math. Gen. 11 983
 [21] Sherrington D and Kirkpatrick S 1975 Phys. Rev. Lett. 35 1792
 [22] Pagnani A, Parisi G and Ratiéville M 2003 Phys. Rev. E 68 046706
 [23] Jörg T, Katzgraber H G and Krza¸kała F 2008 Phys. Rev. Lett. 100 197202
 [24] Mézard M, Parisi G and Virasoro M A 1987 Spin Glass Theory and Beyond (Singapore: World Scientific)
 [25] Fisher K H and Hertz J A 1991 Spin Glasses (Cambridge: Cambridge University Press)
 [26] Nishimori H 2001 Statistical Physics of Spin Glasses and Information Processing: An Introduction (Oxford University Press)
 [27] Castellani T, Krzakala F and RicciTersenghi F 2005 Eur. Phys. J. B 47 99
 [28] Krza¸kała F 2005 Prog. Theor. Phys. Supp. No. 157 77
 [29] Martin O C, Mézard M and Rivoire O 2005 J. Stat. Mech. 2005 P09006
 [30] Rivoire O, Biroli G, Martin O C and Mezard M 2003 Eur. Phys. J. B 37 55
 [31] Zdeborová L 2009 Acta Physica Slovaca 59 169
 [32] Kwon C and Thouless D J 1988 Phys. Rev. B 37 7649