Correlation functions of Ising spins on thin graphs
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] A p r Correlation functions of Ising spins on thin graphs
Piotr Bialas
1, 2, ∗ and Andrzej K. Ole´s † Marian Smoluchowski Institute of Physics, Jagellonian University, Reymonta 4, 30–059 Krakow, Poland Mark Kac Complex Systems Research Centre, Faculty of Physics, Astronomy and Applied Computer Science,Jagellonian University, Reymonta 4, 30–059 Krakow, Poland
We investigate analytically and numerically an Ising spin model with ferromagnetic couplingdefined on random graphs corresponding to Feynman diagrams of a φ q field theory, which exhibits amean field phase transition. We explicitly calculate the correlation functions both in the symmetricand in the broken symmetry phase in the large volume limit. They agree with the results for finitesize systems obtained from Monte Carlo simulations. PACS numbers: 89.75.Hc, 05.10.–a, 05.90.+m
I. INTRODUCTION
The definition of a correlation function on random ge-ometries is not so obvious as in the usual fixed geometrycase. The problem is that we cannot just take two fixedpoints at a given distance because the distance betweenthem is changing. One possible solution is to sum overall pairs of points at a given distance G AB ( r ) ≡ n *X ij A i B j δ d ( i,j ) ,r + . (1)Here d ( i, j ) is the geodesic distance, i.e. , the shortestpath between points i and j . A and B are some quanti-ties defined in each vertex whose correlation we want tomeasure. The average is taken over all instances of thegeometry (configurations) and n is the size of the system.Please note that the distance d ( i, j ) is a very nonlocalquantity: in principle it does depend on the whole con-figuration not only on the endpoints i and j . Thereforethe correlation (1) is not a two-point function, which canlead to some interesting and non-intuitive behavior.Due to this non-locality the calculation of function (1)is in general very difficult. As in other fields of statisti-cal physics an insight can be gained by examining simplesolvable models, the best known being probably the Isingmodel. In this paper we study the Ising model on a ran-dom geometry ensemble made of all Feynman diagramsof a φ q theory.The paper is organized as follows: in Sec. II we presentand analytically solve the model. We proceed on inSec. III calculating the spin-spin correlation functionsusing the similarity to Cayley trees in the infinite vol-ume limit. Section IV addresses the influence of spinson geometry. Final discussion and summary are given inSec. V. ∗ [email protected] † [email protected] II. THE MODEL
We consider an Ising ferromagnet spin model on q -regular random graphs corresponding to φ q Feynman di-agrams of a zero-dimensional field theory [1, 2]. Thepartition function of this model is defined as a sum overall φ q Feynman diagrams G and all the values that thespins s , . . . , s n on a graph G can take Z = X G ∈G X s ,...,s n e − βH ( G ; s ,...,s n ) . (2)If no external magnetic field is present the energy of thesystem on a single diagram G reads H ( G ; s , . . . , s n ) = − X h i,j i∈ G s i s j , (3)where the sum is over all nearest neighbor pairs. Twovertices are considered as nearest neighbors if they areconnected by a link. This includes loops, in which case avertex is its own neighbor, and multiple links when twovertices are counted multiple times in the sum (3).Associating with the “up” and “down” spins the φ + and φ − fields respectively, we can generate the requisiteensemble from the Feynman diagram expansion of thepartition function Z = Z d φ + d φ − exp (cid:20) − ~φ T ∆ − ~φ + 1 q ! (cid:0) φ q + + φ q − (cid:1)(cid:21) , (4)where ~φ T ≡ ( φ + , φ − ) and ∆ is the transfer matrix∆ ≡ (cid:18) e β e − β e − β e β (cid:19) . (5)Following [1] we define the coupling constant g ≡ e β .We will use the saddle-point approximation method tocalculate the partition function Z in the large n limit. Westart by performing binomial expansion of the φ q termsin Eq. (4)exp (cid:20) q ! (cid:0) φ q + + φ q − (cid:1)(cid:21) = ∞ X n =0 n ! (cid:20) q ! (cid:0) φ q + + φ q − (cid:1)(cid:21) n = ∞ X n =0 q !) n n X k =0 φ qk + φ q ( n − k ) − k !( n − k )! . (6)Substituting this back into Eq. (4) and taking the termfor a specific n only we get the partition function of anensemble of graphs with exactly n vertices Z n = 1( q !) n n X k =0 C ( k, n ) Z dφ + dφ − e − ~φ T ∆ − ~φ φ qk + φ q ( n − k ) − (7)where C ( k, n ) ≡ k !( n − k )! . (8)At this point it is convenient to introduce the rescaledfields ψ ± and the variable zψ ± ≡ √ n φ ± , z ≡ kn , and rewrite Eq. (7) as Z n = C ( n ) Z dz C ( z, n ) Z dψ + dψ − e nf ( ψ + ,ψ − ) , (9)where C ( n ) is a constant depending on n only and f ( ψ + , ψ − ) = − ~ψ T ∆ − ~ψ + q [ z ln ψ + + (1 − z ) ln ψ − ] . (10)The field integral in Eq. (9) can be asymptotically ap-proximated using the saddle-point approximation Z dψ + dψ − e nf ( ψ + ,ψ − ) ≈ e nf ( ψ + ,ψ − ) πn √ det J , (11)where J is the Jacobian matrix J xy = f ψ x ψ y . We usethe shorthand notation f ψ x ψ y ≡ ∂ f ( ψ + , ψ − ) ∂ψ x ∂ψ y (cid:12)(cid:12)(cid:12)(cid:12) ψ x = ψ x ,ψ y = ψ y . (12) ψ ± are given by the saddle-point equations ∂f ( ψ + , ψ − ) ∂ψ ± = 0 , (13)which can be written in matrix form∆ − (cid:18) ψ + ψ − (cid:19) = (cid:18) q z/ψ + q (1 − z ) /ψ − (cid:19) . (14) The above system of two quadratic equations has in gen-eral four solutions, from which only the one positive inboth ψ + and ψ − has physical meaning ψ + ( z ) = g − r q ·· h z (cid:0) g − (cid:1) + p g −
1) (1 − z ) z i ,ψ − ( z ) = ψ + (1 − z ) . (15)Because of Eq. (14) f ( ¯ ψ + , ¯ ψ − ) simplifies to f ( ¯ ψ + , ¯ ψ − ) = q (cid:20) z ln ¯ ψ + + (1 − z ) ¯ ψ − − (cid:21) . (16)The logarithm of the C ( z, n ) term can be approximatedusing the Stirling formulaln C ( z, n ) ≈ − n [ z ln z + (1 − z ) ln (1 − z )] −
12 ln z (1 − z ) + ln (2 πn ) + n ln n − n. (17)So finally Z n ≈ C ( n ) Z d z e nF ( z ) (18)and combining Eqs. (9), (16), and (17) we obtain F ( z )to the leading order in nF ( z ) ≈ − z ln z − (1 − z ) ln (1 − z )+ q (cid:2) z ln ¯ ψ + + (1 − z ) ln ¯ ψ − (cid:3) . (19)The next non-leading term is given in Appendix A.Function F ( z ) is symmetric around z = 1 / F ′′ (1 /
2) = 0, from which wecan calculate the critical coupling g c = qq − . (20)For g > g c the F ( z ) function has a minimum at z = 1 / z = 1 / ± m ). In saddle-pointapproximation m is the average magnetization h M i ≡ k − n = n (2 z −
1) = nm. (21)Finding the zeros of the first derivative of F ( z ) we getfor q = 3 m = gg − r g − g + 1 (22)and m = gg − p g − q = 4 case. The susceptibility χ = n ( (cid:10) M (cid:11) −h M i )can be calculated in saddle-point approximation from theintegral (cid:10) M (cid:11) = n R d z (2 z − e nF ( z ) R dz e nF ( z ) , (24)leading to χ = − F ′′ (cid:0) (1 + m ) (cid:1) . (25)In the symmetric phase where m = 0 we obtain χ = (cid:18) − q g − g (cid:19) − . (26)In the broken symmetry phase we get for q = 3 χ = 4 g ( g − g − ( g + 1) . (27)Actually, for comparison with Monte Carlo (MC) simu-lations we will use˜ χ = 1 n ( (cid:10) M (cid:11) − h| M |i ) (28)instead of χ because average magnetization h M i is notwell defined in numerical simulations. On a finite latticeit is in principle zero for all values of β . However, in thebroken symmetry phase its measured value will in generaldepend on the algorithm used and the duration of thesimulation. The average absolute value of magnetizationis given by the integral h| M |i = n R d z | (2 z − | e nF ( z ) R dz e nF ( z ) (29)and in the saddle-point approximation in the symmetricphase it is equal to 1 n h| M |i = 2 π χ. (30)In the broken symmetry phase h| M |i = n m . We plotthe resulting expression for ˜ χ (dashed line) together withdata obtained from MC simulations in Fig. 1.Instead of doing saddle-point approximation of the in-tegrals (24) and (29) we can integrate them numerically.The results are plotted in Fig. 1 with solid lines. As onecan see the agreement is very good already for small sys-tems even if we keep the leading order terms only givenby Eq. (19). Including higher order terms does not im-prove the result significantly. III. CORRELATIONS
In this section we will calculate the connected spin-spincorrelation function G ssc ( r ) ≡ n *X i,j ( s i − m )( s j − m ) δ d ( i,j ) ,r + . (31) χ n = 1000n = 8000 FIG. 1. ˜ χ as a function of the coupling g for q = 3. Cir-cles (blue) render MC data for graphs of size n = 1000 andsquares (red) for those ones of size n = 8000 vertices. Theerror bars show statistical uncertainty which increases nearthe transition. The solid lines plot the corresponding resultsof numerical integration whereas the saddle-point approxima-tion is given by the dashed line. This is not the only way to define the connected correla-tion function on random geometry [3]. We will use thisone because it has the usual property X r G ssc ( r ) = χ. (32)Expanding the expression under the average we can ex-press (31) as G ssc ( r ) = 1 n *X i,j s i s j δ d ( i,j ) ,r + − m n *X i,j s i δ d ( i,j ) ,r + + m n *X i,j δ d ( i,j ) ,r + ≡ G ss ( r ) − m G s ( r ) + m G ( r ) . (33)To calculate the correlation functions we will use thefact that in the leading order of n the solution of themodel coincides with the Bethe solution of the Isingmodel on Cayley trees with coordination number equalto q [1]. Contrary to Cayley trees the Ising model onFeynman diagrams studied here exhibits a genuine phasetransition.The G ( r ) correlation function is a volume-factor: itis the average number of vertices at the distance r froma given vertex. On a Cayley tree with fixed geometry itis easy to calculate G ( r ) = q ( q − r − . (34) G (r) / ( · r- ) n = 16 000n = 32 000n = 64 000n = 128 000n = 256 000 FIG. 2. Ratio of the correlation function G ( r ) obtainedfrom MC simulations of ensembles of various sizes with q = 3at β = β c = ln 3 to its infinite volume limit q ( q − r − . In Fig. 2 we plot the ratio of the measured correlationfunction to the infinite volume limit (34). Please notethe scaling relation evident from the plot2 G ( r ; n ) = G ( r + 1; 2 n ) . (35)This can be also written as G ( r ) = n F (cid:18) r n (cid:19) = n e F ( r ln 2 − ln n ) . (36)We have checked this relation and found that it is verywell fulfilled already for graphs as small as 1000 vertices.The Ising model on Cayley trees is generated by theequations Φ + = 1( q − (cid:18) √ g Φ q − + 1 √ g Φ q − − (cid:19) Φ − = 1( q − (cid:18) √ g Φ q − + √ g Φ q − − (cid:19) (37)whose graphical interpretation is shown in Fig. 3. Theseequations are analogous to the ones obtained for Isingmodel on branched polymers [4]. The difference is thathere the probability of ending the branch is zero. Thisformally means that we have only infinite trees. Hencethere is no chemical potential associated with every ver-tex.The partition function can be then obtained as Z = 1 q ! (cid:0) Φ q + + Φ q − (cid:1) . (38)The system (37) will have in general many pairs of so-lutions. However, for g ≤ g c the system is in symmetric == ++ FIG. 3. Graphical representation of the saddle-point equa-tions given by Eqs. (37) for q = 3. The bright bubbles corre-spond to Φ + and the dark ones to Φ − . Smaller bright anddark points stand for vertices carrying the “up” and “down”spins respectively.FIG. 4. Graphical representation of the correlation function Z ( r ) xy given by Eq. (42) for q = 3. Bubbles depict Φ ± andpoints stand for individual vertices connected by links associ-ated with ∆. Only the first and the last vertex have specificspins and each of them contributes by a Φ ± factor. The spinof the r − M . phase and Φ + = Φ − . The one resulting equation can beeasily solved yieldingΦ q − ± = √ g g ( q − β ( q − . (39)When g > g c the dominant solution of Eqs. (37) will haveΦ + = Φ − giving a non-zero magnetization [2] m = Φ q + − Φ q − Φ q + + Φ q − . (40)Finding this solution requires solving the system (37).Although this cannot be done analytically for general q for q = 3 and q = 4 we get in the broken phaseΦ ± = √ gg − (cid:18) ± r g − g + 1 (cid:19) for q = 3 , √ g (cid:16) g ± p g − (cid:17) g − for q = 4 . (41)Substituting this into Eq. (40) we obtain the magnetiza-tion which is identical to the expressions (22) and (23)calculated using the previous method.We will derive the correlation functions by the methoddescribed in [5]. The correlation function can be repre-sented graphically as in Fig. 4. To this picture corre-sponds the expression Z ( r ) xy = 1[( q − Φ q − x r − z }| { ∆ M · · · ∆ M ∆ xy Φ q − y , (42)where M ≡ q − (cid:18) Φ q −
00 Φ q − − (cid:19) . (43)If we rewrite matrix M as M = (Φ + Φ − ) q − ( q − (cid:16) Φ + Φ − (cid:17) q − (cid:16) Φ − Φ + (cid:17) q − (44)we can see that the correlation functions are related tothe correlation functions of the Ising spin chain in theeffective magnetic field ( q − h eff with [6] h eff ≡
12 ln Φ + Φ − . (45)In the symmetric case in particular h eff = 0.The required correlation functions can be expressedusing Z ( r ) xy as G ( r ) = 1 Z X xy Z ( r ) xy , (46) G s ( r ) = 1 Z X xy yZ ( r ) xy , (47) G ss ( r ) = 1 Z X xy xyZ ( r ) xy . (48)After introducing matrices M = 1 p ( q − Φ q − +
00 Φ q − − ! (49)and e Q ≡ M ∆ M (50) Z ( r ) xy given by Eq. (42) can be rewritten in the form Z ( r ) xy = ( q − q − x M − xx ( q − (cid:16) M ∆ M (cid:17) rxy M − yy Φ q − y ( q − q − q x ( q − e Q rxy Φ q y ( q − . (51) Assuming that the eigenvalues of e Q are λ and λ , andthat the corresponding normalized eigenvectors are ( a, b )and ( − b, a ) respectively, we can calculate the matrixpower by diagonalizing e Q .After some algebra we obtain: G ( r ) = 1 Z λ r ( a Φ q + + b Φ q − ) + λ r ( b Φ q + − a Φ q − ) ( q − q − , (52) G s ( r ) = 1 Z λ r ( a Φ q + − b Φ q − ) + λ r ( b Φ q + − a Φ q − )( q − q − , (53) G ss ( r ) = 1 Z λ r ( a Φ q + − b Φ q − ) + λ r ( b Φ q + + a Φ q − ) ( q − q − . (54)In Appendix B we show that λ = q − , (55) λ = ( q − (cid:20) √ g ( q − q − + Φ q − − ) − (cid:21) = ( q − (cid:20) √ g Φ q − ( q − q − h eff ) − (cid:21) ≡ ( q − e λ , (56)where Φ ≡ p Φ + Φ − . (57)The corresponding eigenvector is equal to( a, b ) = (Φ q + , Φ q − ) q Φ q + + Φ q − . (58)Substituting this into Eq. (52) we get G ( r ) = 1 Z ( q − r − ( q − q + + Φ q − ) . (59)Using the definition (38) of Z we finally obtain the result(34) which is a check of the consistency of the methodused.For spin-spin correlation functions we obtain G s ( r ) = q ! ( q − r ( q − q −
1) Φ q + − Φ q − Φ q + + Φ q − = m G ( r ) , (60) G ss ( r ) = q ! ( q − r (Φ q + − Φ q − ) + 4Φ q + Φ q − λ r ( q − q − q + + Φ q − ) = G ( r ) (cid:20) q + Φ q − (Φ q + + Φ q − ) e λ r + m (cid:21) = G ( r ) e λ r cosh ( q h eff ) + m ! . (61)In the symmetric case when m = 0 and Φ + = Φ − aregiven by Eq. (39) we readily get g ss ( r ) ≡ G ss ( r ) G ( r ) = (cid:18) g − g + 1 (cid:19) r = tanh r β. (62) g ss (r) β = 0.45 β = β c β = 0.60 FIG. 5. Correlation function g ss ( r ) for q = 3. Symbols renderMC data for graphs of size n = 256000 in the symmetric phasefor β = 0 .
45 (circles), at the transition β = β c = ln 3 (squares),and in the broken phase for β = 0 .
60 (diamonds). The solidlines plot the analytical predictions given by Eqs. (62) and(63).
This is as predicted the result obtained for correlation ofIsing spins on the chain [6]. It is easy to check that therelation (32) is satisfied by the above function.In the broken phase inserting Eqs. (41) into Eq. (60)we obtain for q = 3 g ss ( r ) = 4( g − ( g + 1) 1( g − r + m (63)and for q = 4 g ss ( r ) = 4( g − g − r + m . (64)Although similar formulas can be derived for higher val-ues of q they are much more complicated.Figure 5 plots the correlation function g ss ( r ) for q = 3and different values of β . As one can see the agreementwith MC results is very good. In Fig. 6 correlation func-tions for different q ’s at the respective transition points β c ( q ) = ln qq − are compared. Again, MC results matchthe asymptotic predictions. The small discrepancies di-minish with the increase of the graphs’ size. IV. INFLUENCE OF SPINS ON GEOMETRY
On a finite tree spins can be integrated out exactly andthe resulting factor does not depend on the shape of thetree. This means that the spins do not have any influenceon the geometry. The situation changes when cycles areallowed, which is especially obvious for loops, i.e. , links g ss (r) q = 3q = 4q = 6 FIG. 6. Correlation function g ss ( r ) at the transition β = β c ( q )for different values of q . Points represent MC results forgraphs of size n = 64000 while the solid lines plot the ana-lytical predictions. Please note that for q = 6 the maximumencountered distance (the diameter of the graph) is r = 9. β n l FIG. 7. Average number of loops n l in graphs of size n = 32000vertices and q = 3. The black points represent MC data andthe dashed red line marks the transition point β = β c . attached at both ends to the same vertex. Without anyspins or β = 0 one can show combinatorially that the ex-pected number of such loops on a graph equals ( q − / β , however, we should observe an enhance-ment as each loop contributes the e β factor to the par-tition function contrary to links joining different verticeswhich can have different spins. For β → ∞ all the spinshave the same sign and again the geometry decouples.In Fig. 7 we plot the average number of loops n l as afunction of β . The results agree qualitatively with theabove scenario. Nevertheless, one should note that whilelooking pronounced this is still only a 1 /n effect. Thenumber of loops is independent of n and negligible in thelarge volume limit. V. SUMMARY AND DISCUSSION
We have analyzed in detail a simple model of spins ona random geometry. In the large n limit it is formallyequivalent to the Ising model on an infinite Cayley tree.However, it is well defined, has a genuine phase tran-sition, and can be easily simulated using Monte Carlomethods.We have derived expressions for correlation functionsin both the symmetric and the broken phase with meth-ods less formal then in [6]. From these calculations weobtain a picture of the transition. The correlation func-tion g ss ( r ) does not exhibit any critical behavior: thecorrelation length is finite for any finite β . Neverthe-less, the volume factor G ( r ) grows exponentially andthis growth offsets the decay of the correlation function g ss ( r ). As we increase r by one, the influence of the spinsat this distance drops by ( g − / ( g + 1), however, thenumber of these spins increases by q −
1. When g + 1 g − q − g . ACKNOWLEDGMENTS
We would like to thank Zdzis law Burda for many help-ful discussions. This research was supported in part bythe PL-Grid Infrastructure. MC simulations were per-formed on the Shiva computing cluster at the Facultyof Physics, Astronomy and Applied Computer Science,Jagellonian University, and at the Academic ComputerCentre CYFRONET AGH using the Zeus cluster.
Appendix A: Next to leading corrections to F ( z ) The Jacobian matrix J xy = f ψ x ψ y is equal to J = ∆ − − q (cid:18) z/ψ +
00 (1 − z ) /ψ − (cid:19) , (A1)from whichdet J = gg − " q √ g zψ + 1 − zψ − ! + q z (1 − z ) (cid:0) ψ − ψ + (cid:1) . (A2)and the approximation of the function F ( Z ) given byEq. (19) takes the form F ( z ) ≈ − [ z ln z + (1 − z ) ln (1 − z )]+ q (cid:2) z ln ¯ ψ + + (1 − z ) ln ¯ ψ − (cid:3) − n { ln (det J ) + ln [ z (1 − z )] } . (A3) Appendix B: Eigenvalues
Equations (37) can be rewritten as (cid:18) Φ + Φ − (cid:19) = 1( q − √ g Φ q −
2+ 1 √ g Φ q − − √ g Φ q − √ g Φ q − − ! (cid:18) Φ + Φ − (cid:19) ≡ Q · (cid:18) Φ + Φ − (cid:19) . (B1)The above equation means that matrix Q has an eigen-value equal to one and that the corresponding eigenvectoris proportional to (Φ + , Φ − ).Matrix e Q defined in Eq. (50) can be expressed as e Q xy = ( q − Q xy (cid:18) Φ x Φ y (cid:19) q − . (B2)One can check that this implies that e Q has an eigenvalue λ = q − a, b ) = (Φ q + , Φ q − ) q Φ q + + Φ q − . (B3)The second eigenvector is perpendicular to this one. Mul-tiplying it by e Q we obtain the second eigenvalue. [1] C. F. Baillie, D. A. Johnston, and J-P. Kownacki, Nucl.Phys. B , 551 (1994).[2] C. Bachas, C. de Calan, P.M.S. Petroupoulos, J. Phys. A , 6121 (1994).[3] B. V. de Bakker and J. Smit, Nucl. Phys. B , 343(1995). [4] J. Ambjorn, B. Durhuus, T. Jonsson, and G. Thorleifsson,Nucl. Phys. B , 568 (1993).[5] P. Bialas, Nucl. Phys. B , 645 (2000).[6] T. Morita and T. Horiguchi, Prog. Theor. Phys.54