Glassy phase of the frustrated Bethe lattice models

As a short range frustrated spin model, we study the ferromagnetic spin models on the Bethe lattice with next nearest neighbor (n.n.n.) interactions with negative sign. The iteration equations which are similar to belief propagation are derived and studied by population dynamics. At low temperatures, we find the ferromagnetic phase for zero or small n.n.n. interactions. As the strength of them increases, ferromagnetic phase transition temperature becomes zero and non-trivial distributions of the local fields arise, implying that there are glassy phases in these models.


Introduction
Frustrated spin models without randomness have been studied for a long time in the study of statistical physics [1][2][3][4]. The mixture of ferromagentic and antiferromagnetic interactions easily creates frustrations and bring about some complex low temperature phases [5]. One of the interesting possibility is the existence of the glassy states [6], similar to spin glass phase which has been studied for random spin models by using the mean field method [7,8]. About a decade ago, using the ideas from the spin glass theory, it was found that the long range antiferromagnetic(LRAF) spin models have glassy low temperature phase [9,10]. The modeling of LRAF is quite simple. By introducing quadratic constraints on the Fourier components of spin variables, we obtain the long range antiferromagnetic interactions if the number of constrained components are the same order of system size. Simulated annealing and replica theory imply that these models have really glassy low temperature phase as well as the ordered phase.
The properties of the LRAF can be very general since the key concept is the frustrations caused by the negative interactions. Various spin models may have this property in common. For example, in the ferromagnetic spin model on the square or simple cubic lattice, frustrations are easily created by introducing the next nearest neighbor(n.n.n.) interactions with negative sign. In contrast to the simplicity, the property of these models does not seem to be well understood, probably because of the possible complexity of the low temperature phase.
The studies of the LRAF are based on the replica mean field theory, which is applicable to the long range spin models. Since the idea of frustration is not limited to the long range models, it is highly desirable to clarify the situation of the spin models with short or middle range interactions. For this purpose, we should figure out the frustrated spin models with finite range interactions which can be treated by the mean field approach. In this respect, the spin models with finite connectivity or defined on the tree structure may the next subject to study [11,12].
We note that the frustration can be created in the lattice models in any dimension. It is known that the lattice spin models in high enough dimension is approximated to the spin models on the Bethe lattice [13]. Generally, on the Bethe lattice, due to the tree structure, the statistical summation can be performed by constructing the iteration equation, whose behavior reflects the property of the system. We note that, on the bipartite lattice, including the Bethe lattice, the frustrated spin models may be created by introducing the negative n.n.n. interactions. As we will see, introduction of n.n.n. interactions does not seem to destroy the basic property of the tree structure. Thus we can apply the iteration method, which is similar to the ones explored for statistical physics approach to various problems [12,14].
The purpose of this paper is to suggest and study the ferromagnetic Bethe lattice models with n.n.n. interactions with negative sign. In section 2, we introduce the iteration equations for this model. In section 3, we discuss the iteration equations to find the ferromagnetic phase transition. Sections 4 and 5 are devoted to the study of population dynamics for the iteration equations to study the glassy phase. Section 6 is devoted to some discussions.

Frustrated spin model on the Bethe lattice
Let us first describe the spin model defined on the Bethe lattice and introduce some notations [13]. Image a central site c and add p sites connected to c. which are numbered by i=1, 2 ..., p. This is the first shell. Taking a site in the first shell and connecting p−1 new sites to it, we have a lattice with the second shell. Repeating this procedure L times, we have a Bethe lattice with depth L which has + ---(( ) ) ( ) p p p 1 1 1 2 L sites. By this construction, interior sites have p neighbors, while boundary sites have only one neighbor. We should study the property of the interior system to understand the spin model with no boundary effect. This may be achieved by adopting the fixed point values of the iteration equation to the boundary sites.
The spin models are defined by assigning the interactions between a site in shell k and the connected sites in shell k +1 for all k. These interactions are nearest neighbor(n.n.) interactions. The next nearest neighbor pair is naturally defined by the distance two, i.e., having one intermediate site between the two sites of the pair. Assuming the interactions between them, we can define the spin model with n.n.n interactions. This structure is illustrated in figure 1 for p=3 and L=3.
Denoting the n.n. interactions by J ij and n.n.n. interactions by K kl , the energy function of the system is given by where we assume Ising spins σ i =±1. The n.n. interactions are ferromagnetic, i.e., J ij >0, while n.n.n. interactions are antiferromagnetic with K kl >0. We assume that they take constant values J and K respectively, yet keep the indices in the following equations to show the connected sites. We can take negative J ij to have full antiferromagnetic models. However, it will be convenient to have a well known ferromagnetic theory for small K to check the consistency of the theory.
For p=2, our model reduces to the spin model on the one dimensional lattice with n.n.n interactions, which only has a paramagnetic phase at all fintie temperatures. If we add the longer range interactions further, the situation may become similar to the spin model discussed in the previous paper [9], where glassy phase appears at low temperatures. For p>2, the connectivity increases even with only n.n.n interactions, and we expect the similar condensed phase arises at low temperatures Now let us describe the procedure to obtain the iteration equation for the system, For the sake of simplicity, we assume that the inverse temperature β=T −1 is implicitly included in interactions. Just as in the usual Bethe lattice models, summation over the spins on shell L gives a Bethe lattice with L−1 depth with modified interactions. Due to n.n.n interactions, this creates local fields not only on the site in shell L−1 but also on the connected site in shell L−2. In addition, there arises the additional interaction between these sites by this procedure. Thus the result of summation on shell L is described by the amplitude which is a function of the spin pair on shell L−1 and shell L−2. This property is also valid for any step of summation. This procedure will be described by the equations similar to belief propagation. The iteration equation is obtained as follows. Take a site in shell k, which is denoted by 0. Then there are p−1 n.n. sites which belong to shell k+1, numbered as j=2, 3, L, p and one n.n. site which belongs to shell k−1, which is denoted by a number 1. We introduce amplitudes n s s = for bond ( j0), which reflects the statistical summation coming from shells l (>k+1). Then new amplitude n s s  ( ) , 0 1 0 1 for the bond (01) from shell k to shell k−1 should satisfy the equation given by where the summation in the exponential is over the sites i, j=2, L, p with i<j, and  z 0 1 is a normalization constant. Note that the notations such as  j 0 simply mean the directed bonds from a shell to lower shell. Let us briefly comment on the expression of the free energy. Although we shall not complete the calculation, the resulting expression will be helpful to see the consistency of the theory. As described above, the summation over spins on shell k+1 creates a function n s s   ( ) z , 0 1 0 1 0 1 for the bond (01) between shell k and shell k−1. Thus a constant  z 0 1 arises to the partition function by the statistical sum for the polygonal unit made of j=2, 3, L, p. Then the partition function is expressed by the product of  z j i , where i, j are nearest neighbor pairs. In this way we obtain where bonds are directed n.n. pairs. To avoid the effect of boundary, all parameters should be these of fixed point of iteration.
, e x p , 6 01 , is the energy function of the polygonal unit made of σ i (i=0, 1, 2, L, p). The first equation can be regarded as the effective statistical sum of bond (01). The second equation can be regarded as an effective statistical sum of polygonal unit made of σ i (i=0, 1, 2, L, p). In this way, we have The a-sum is over all sites and (i j)-sum is over all n.n. bonds. This expression shows that effective free energy z ln a works as a mean field free energy of the model. However, the summation of them over all sites counts the contribution from each n.n. bond twice, and should be subtracted once, which is performed by the second term. To obtain the thermodynamic limit of the free energy, the expression above should be averaged over the distribution of n s s  ( ) , j i j i , which will be studied in sections 4 and 5. We shall not discuss the free energy further in this paper, and concentrate on the phase diagram in the following sections.

Solution of the iteration equation
For the Ising model, we can express the iteration equation in a concise form. For two Ising variables, we need three parameters, denoted by a j , b j , and c j , to express the normalized n s s  ( ) , Note that b j reduces to the local fields of usual Bethe lattice model for K=0. New parameters a j and c j are created by the n.n.n. antiferromagnetic interactions. a j contributes to the n.n. interactions and c j contributes to b of the next iteration. d j should be adjusted to normalize the amplitude. Substituting n s s  ( ) , It is convenient to define the partition function: Then we obtain the iteration equations as This completes the iteration equations for n.n.n. Bethe lattice model. So far we have discussed the formulation for general p. The sums over σ i in z(σ 0 , σ 1 ) are straightforward for small p. As expected, the expressions become complicated as p increases. We shall present the results for p=3 in this and next section and discuss the situation for p=4 in section 5.
For p=3, we have s s = - 2 exp  cosh  exp  cosh  17   3 0 1  2  3  2  3 for z(σ 0 , σ 1 ). The temperature T is introduced by assuming the relation J=1/T and K=r/T, where we have introduced a parameter which measures the strength of frustration. We first discuss the ferromagnetic transition point. Since the parameters take definite values in the ferromagnetic phase, we assume a 2 =a 3 =a, b 2 =b 3 =b and c 2 =c 3 =c. As seen by the high temperature expansion, we note that a is different from zero even in the paramagnetic phase when ¹ K 0. Substituting b=c=0, we obtain the self-consistent equation for a = + - To find the transition point, we first study the equations to the first order of b and c. To this order, the selfconsistent equations reduce to: The equation for a does not change to this order. Putting these equations together, we obtain the equation for the transition temperature: . 23 For K=0, this equation reduces to that of the pure n.n.model, where a=0 and the ferromagnetic transition temperature T f is given by = ( ) T 1 2 tanh 1 f , which yields T f =1.82.
For K>0, (23) reveals that the ferromagnetic transition temperature decreases quite rapidly as K increases. Actually ferromagnetic phase ceases to exist even at zero temperature for r>r c , where the right hand side of (23) does not become larger than 1 for all temperatures. We found r c =0.245 where the right hand side of (23) becomes 1 at T=0. For r>r c , we expect the glassy phase appears at low temperatures. In the next section, we shall discuss the population dynamics to confirm the result in this section and to find the glassy phase.

Population dynamics
In the previous section, we have shown that the ferromagnetic phase ceases to exist at r=r c for p=3. For larger r, we expect that there is a glassy phase at low temperatures, where the local fields will be distributed. Accordingly, we should assume ¹ b b 2 3 , ¹ c c 2 3 , and ¹ a a 2 3 in the iteration equations, expecting that there are no fixed point values for them, but various values appear in the iteration equations (14)-(16). Let us assume that the distribution of º ( ) a a b c , , is described by the distribution function ( ) a P . Our goal is to find the fixed point disribution ( ) a P which satisfies where ( ) f a a , 1 2 stands for a set of iteration equations, and d a b c . The numerical study of ( ) a P will be performed by the population dynamics [12], assuming the population a n , , where  is a large enough integer. Then we assume The procedure of the population dynamics is as follows: After assuming an initial distribution a n , (1) choose two a n , which are denoted by a k and a l .
(2) calculate a new a for a k and a l by using (14)-(16).
(3) randomly chosen a p is replaced by the new a.
Repeat the steps (1)-(3) until some convergence condition is satisfied. We take  =10 10 4 5 and the number of iterations is typically 10 10 in the numerical simulations. In the description above, we have not commented on the nature of the initial distributions. For p=3, we simply assume that initial b n are randomly distributed with zero mean and varience of order 1, and a n and c n are values of order 1.
To find the phase transition, we study the average and variance of b with ( ) a P , which are denoted byb and b 2 . In the ferromagnetic and paramagnetic phases, ( ) a P is simply given by one delta function with one a, which implies At high temperatures, we found thatb andb 2 tend to exponentially small values, which corresponds to a paramagnetic phase. At low temperatures, the behavior ofb andb 2 depends on r=K/J.
For p=3, the boundary from paramagnetic to ferromagnetic phase obtained by population dynamics is in good agreement with the equation (23) down to very low temperatures. Namely, for small enough 0<r<r c ∼0.245 and at low temperatures, we obtain precise relation quite rapidly by population dynamics. These values do not fluctuate after reaching a fixed point value. This implies that ( ) a P is a single delta function, which corresponds to a ferromagnetic phase. The ferromagnetic transition temperature is denoted by T f . Just above r c =0.245, we found a paramagnetic phase down to zero temperature. For r>r g ∼0.255, there appears a phase withb 0, > b 0 2 , that is, the distribution of b is different from a single delta function. The transition temperature is denoted by T g , above whichb 2 tends to exponentially small values. Below T g , in contrast to the paramagnetic phase,b 2 continues to fluctuate around some positive value and does not show regular behavior in the studied number of steps. Non-trivial distribution of b implies that we have a random spin phase below T g . although our model does not have any randomness. Thus our iteration equations distinguish glassy phase from paramagnetic phase. The r-dependence of the critical temperatures T f and T g is depicted in figure 2.
A lot of studies may be required to clarify the reason why b and c are distributed. Here we only remark that the right hand side of (23) becomes negative for large r, implying that the fixed point of iteration, if it exists, is unstable. This may be the reason why b and c take various values.

p=4 case
The study for p=4 case may require another section to describe. Due to the non-trivial crossover region from the ferromagnetic to spin glass phase, the situation seems more complicated, richer and harder to study than the case p=3. We shall mainly discuss this point in this section.
For p=4, we have s s = -+ + + -+ + + -+ + + - for z(σ 0 , σ 1 ). The iteration equations are given by (14)-(16). The ferromagnetic phase transition is found by studying the self-consistent equation to the first order of b and c, as was done in section 3. As expected, the equations are rather complicated. We first remark that they reduce to = ( ) T 1 3 tanh 1 f for r=0, yielding T f =2.88, and T f becomes zero at r=r c ∼0.165. However, even the ferromagnetic phase does not seem to so simple as in the p=3 case. Actually, the population dynamics implies that the transition is discontinuous for r>∼0.15 and T f is finite up to r∼0.2 (see figure 3). In addition, the crossover from ferromagnetic to spin glass phase looks rather complicated. These two phases overlap in the range r=0.1∼0.2, where there seem to be at least two fixed points of the population dynamics at low enough temperatures, one is ferromagnetic and other is glassy phase. Let us explain the details of the simulations.
In p=3 case, the initial distributions of b were simply assumed to be random values with zero mean. In p=4 case, after several inspections, we realize that there is a region that the evolution of b seems to depend on the initial conditions. For this reason, we study the evolution by using two types of random distribution for the initial distribution b 0 : Type A is such that~ , which is expected to tend to the ferromagnetic distribution, and Type B is such thatb 0 . O is a number of order one and tends to be larger as the temperature decreases to have a quick convergence to the fixed point distribution. In the most parameter regions of r, T, the difference of initial conditions does not matter much as long as the parameters are not close to the critical line. However, for p=4, this is not the case in the region r=0.1∼0.2 at low temperatures.  We first note that, for r>0.15, the ferromagnetic transition temperature is different from the one obtained by the continuous assumption presented in section 3. Actually, the transition becomes discontinuous for 0.15<r<0.2 and T f is finite up to r∼0.2, as described above. The evolution of b is rather simple below T f as long as the initial condition A is concerned, that is, the relation = > is achieved quickly, implying that ( ) a P becomes the ferromagnetic distribution. However, in the crossover region, taking the initial condition B, we sometimes find that b does not tends to ferromagnetic distribution even below T f . That is, there is a temperature, denoted by T g , below which the initial condition A give a ferromagnetic distribution, while condition B gives rather a distribution which is regarded as glassy phase at least in the studied number of steps. For this reason, we are obliged to have a finite T g . T g smaller than T f are obtained as an upper temperature where B gives a glassy distribution. On the contrary, for T g <T<T f , both initial distributions tend to a ferromagnetic distribution. These aspects hold for 0.1<r<0.16.
For r>0.16, T g becomes greater than T f . In this case, T g is studied by assuming either A or B and found as the lowest T which give exponentially smallb andb 2 . We find a ferromagnetic phase below T g in the narrow region 0.16<r<0.2, giving T f . Below T f , initial conditions A and B give different fixed point distributions. In either case, the behavior ofb andb 2 looks complicated in glassy phase. (¯) b 2 continues to fluctuate around values smaller thanb 2 , implying that the distribution is asymmeric around b=0.
The r-dependences of T g and T f are depicted in figure 3. The line for T f for 0<r<0.165 is obtained by the analytic method described in section 3, and T f for 0.1<r<0.2 is obtained by the population dynamics. There are two T g lines: an isolated line below T f for 0.1<r<0.16 is a temperature below which two different initial conditions give different distributions. For r>0.16, there is another T g which is similar to the p=3 case.
We should remark that the existence of glassy phase in T<T g <T f is not conclusive since the convergence of b to glassy distribution can be an artifact of the finite number of simulation steps, that is, the distribution may tend to that of ferromagnetic one by larger simulation. Actually, at T just above T g , it takes a rather large number of steps to reach ferromagentic distribution for initial condition B. Anyway, the crossover from the ferromagnetic phase to glassy phase is much more complicated and richer in the p=4 case than in the p=3 case.

Discussion
In this paper, we have studied the frustrated spin model on the Beth lattice by introducing the n.n.n. interactions. The belief propagation-like equation is presented and studied by population dynamics. The effects of the n.n.n interactions are taken into account by introducing two local fields and interaction between the two sites. The behavior of them implies that there are two low temperature phases: ferromagnetic phase and glassy phase. We have presented the study for the Bethe lattice with p=3 and 4 nearest neighbors.
For p=3, the phase diagram is simple. There are two low temperature phases, ferromagnetic phase and glassy phase, which are separated by a paramagenetic phase. However, it remains to study the reason why nontrivial distribution of b arises. The situation becomes rather complicated for p=4. Two low temperature phases have a crossover region, where two phase coexist. In this region, the resulting distributions are either ferromagnetic or glassy, depending on the initial distributions. Further, the glassy phase in the crossover region looks rather complicated in the sense that (¯) b 2 is neither zero nor b 2 . However, this result is not conclusive. It can be an artifact of the finite number of simulation steps. In this sense, it is highly desirable to perform the analytic studies of the iteration equations, especially for the glassy phase.
The coexistence of random phase and ordered phase does not seem to be so special. For example, the Hopfield model has two kinds of phase at low temperatures; a phase which remembers the embedded patterns and spin glass phase, which has no such information [15]. The LRAF models described in introduction also has ordered phase and disordered phase. In these models, the transition to ordered phase are discontinuous just like the situation presented in this paper. We suppose that the configuration space in these models is divided into that of ordered phase and that of glassy phase. Comparing the results for p=3 and 4, we expect that the coexistence region will become larger as p increases. Thus the study for larger p will be helpful to understand the crossover.
In this paper, we have restrict ourselves to the study of the phase diagram. The thermodynamic functions in these phases remain to be studied. The relation to replica method will be also an interesting subject, although our model does not have any randomness. Apart from the spin models, the lattice gas models in high dimensions will be discussed in a similar manner as was done in this paper. It will be interesting and meaningful if glassy states are found in the framework of the lattice model without random interactions.