# Solution on the Bethe lattice of a hard core athermal gas with two kinds of particles

###### Abstract

Athermal lattice gases of particles with first neighbor exclusion have been studied for a long time as simple models exhibiting a fluid-solid transition. At low concentration the particles occupy randomly both sublattices, but as the concentration is increased one of the sublattices is occupied preferentially. Here we study a mixed lattice gas with excluded volume interactions only in the grand-canonical formalism with two kinds of particles: small ones, which occupy a single lattice site and large ones, which occupy one site and its first neighbors. We solve the model on a Bethe lattice of arbitrary coordination number . In the parameter space defined by the activities of both particles. At low values of the activity of small particles () we find a continuous transition from the fluid to the solid phase as the activity of large particles () is increased. At higher values of the transition becomes discontinuous, both regimes are separated by a tricritical point. The critical line has a negative slope at and displays a minimum before reaching the tricritical point, so that a reentrant behavior is observed for constant values of in the region of low density of small particles. The isobaric curves of the total density of particles as a function of (or ) show a minimum in the fluid phase.

###### pacs:

05.50.+q,64.60.Kw,67.70.D-## I Introduction

In the beginning of statistical mechanics, the interest was focused mainly on fluids, and the pioneering work on phase transitions concentrated on the liquid-gas transition, an example of this kind is the van der Waals equation of state vdw73 (). In his phenomenological theory, van der Waals considered the effect of hard core excluded volume interactions and attractive interactions in the thermodynamic behavior of the fluid. While attractive interactions between molecules are essential to produce the liquid-gas transition, it was later realized that even if only hard core interactions are considered, interesting effects arise in the models. Since in this case all allowed microscopic configurations of the system have the same energy, such models are athermal, and all thermodynamic effects are of entropic origin. Much is known about the continuous version of these models, usually called hard sphere systems. They were studied by a variety of techniques hd86 (), and a fluid-solid phase transition is found. It is worth recalling that in the pioneering work on the Monte Carlo simulational procedure mrrtt53 (), the physical system studied was a fluid of hard disks.

The fluid model with excluded volume interactions only may also be defined on a lattice, so that the positions occupied by the particles are restricted to sites of a lattice. In the simplest version of such models, the only constraint is that if a site is occupied by one particle, no other particles may be placed on it. In the grand-canonical ensemble, this model reduces to the Ising model without the interaction term, and is therefore trivially solved. No singularities are found in the thermodynamic functions, as expected. More interesting results are obtained if the range of excluded volume interactions is increased. This leads to a variety of models, and we refer to a recent paper where some of these models were studied using simulations for a comprehensive survey of the literature fal07 (). If a particle placed on a site excludes this site and its first neighbors, indeed a fluid-solid continuous transition is found. For bipartite lattices such as the hypercubic lattices, at lower densities, the particles occupy the lattice sites randomly, but above a critical density one of the two sublattices is occupied preferentially by the particles, so that the order parameter may be defined as:

(1) |

where is the number of particles in sublattice or , respectively, divided by the number of sites in the lattice, thus assuming a maximum value equal to . This model has a long history, it has been mentioned in the classical review by Domb d58 (), and has been studied using the virial expansion and the Bethe approximation by Burley shortly after b60 (). Since then, the thermodynamic behavior of the model has been investigated using a variety of analytical and numerical methods, which indicate a continuous phase transition in the Ising universality class. On the square lattice, precise estimates for the critical chemical potential and the density of particles at the transition were obtained using transfer matrix and finite size scaling extrapolation techniques, and the Ising critical exponents were verified with high precision gb02 (). A recent simulational study of this model on the square and cubic lattices may be found in cnd11 (). On the triangular lattice, this model is known as the hard hexagon model, and was solved exactly by Baxter b80 (). A continuous transition was found at the critical activity . The critical exponents are in the 3-state Potts universality class, as would be expected considering that in the high density phase the hexagons occupy preferentially one of the three sublattices.

An interesting generalization of the model is to consider a gas with both small and large particles, where the small particles occupy a single site and the large ones the site and its first neighbors. Let the activity of small (large) particles be (). On a square lattice, we may represent both particles as squares, as is shown in Fig. 1. Since a continuous transition occurs when only large particles are present and no transition is found for small particles only, one may ask what happens in intermediate situations, with both types particles on the lattice. This model was studied on the square lattice by Poland p84 () using high density series expansions in the grand-canonical formalism, and evidence was found that a tricritical point should be present in the phase diagram. In this paper we solve the model on a Bethe lattice of arbitrary coordination number in the grand-canonical ensemble. As the fugacity of the small particles is increased starting from zero, the fluid-solid transition remains continuous up to a certain value, above which it becomes discontinuous. Therefore, we find a tricritical point in the phase diagram of the model.

It is worth mentioning that a slight modification of the two-particle model on the square lattice makes it exactly solvable in a particular case, by allowing it to be mapped on the Ising model for which the exact solution is known for zero magnetic field. This particular case corresponds to fl92 (). In this lattice gas, the large particles are placed on the centers of the elementary squares of the lattice defined by the sites where the small particles are located, as is illustrated in figure 2. Both models are equivalent in the absence of small particles, but we notice that when we know the exact solution only for .

## Ii Definition of the mixed lattice gas model and its solution on the Bethe lattice

We study the grand-canonical version of the mixed lattice gas model defined on a lattice. In this model, two types of particles are present, and only excluded volume interactions between them are considered, so that the model is athermal. The particles of type 1 (small) are such that, when placed on a lattice site, they occupy this single site only, excluding other particles from it. Particles of type 2 (large), when placed on a site, occupy its first neighbor sites also. Figure 1 shows this model on a square lattice and the generalization to other hypercubic lattices is straightforward. While no phase transition is found for the case where only small particles are present, it is well established that for a pure system of large particles a continuous phase transition in the Ising universality class occurs fal07 ().

Here we will solve the model with both particles present in the grand-canonical ensemble on the Bethe lattice. The parameters of the model are the activities of small particles , where is the chemical potential of a small particle divided by , and of large particles . We proceed defining the model on a Bethe lattice, which is the core of a Cayley tree with general coordination number . We then consider partial partition functions (ppf’s) of the model on subtrees with fixed configurations of the root site, which may be empty (), occupied by a small particle () or by a large one (). Considering the operation of connecting subtrees with a certain number of generations of sites to a new root site, we may build a subtree with an additional generation and write down recursion relations for the ppf’s. If we call the partial partition function of a subtree with root site configuration , we have the following recursion relations for the pff’s:

(2a) | |||||

(2b) | |||||

(2c) |

The prime denotes subtrees with an additional generation. Let us define ratios of the ppf’s as , where now the configuration index assumes the values 1 and 2. From the recursion relations for the ppf’s, we may obtain the ones for the ratios, which are:

(3a) | |||||

(3b) |

The thermodynamic behavior of the model is defined by the values of the ratios after a large number of iterations of the recursion relations Eqs. (3). We find that, depending of the values of the activities, the recursion relations converge either to a fixed point or to a limit cycle of period 2. For convenience, we may define the variables and . The fixed point equations in these variables will be:

(4a) | |||||

(4b) |

When the recursion relations converge to a limit cycle, the values of the ratios in the core of the tree will display a layered structure, so that the ratios in sites in consecutive generations assume alternate values. It is convenient, in this case, to define two sublattices (A and B), whose sites are associated to the two values of the pair of variables . The equations defining this limit cycle values are:

(5a) | |||||

(5b) | |||||

(5c) | |||||

(5d) |

Although we were not able to find general solutions for these equations, it is not difficult to solve them numerically. This set of equations may be reduced to finding the roots of a polynomial in the variable , given by:

(6) |

where is given by:

(7) |

The fixed points or limit cycles will correspond to roots of this polynomial in the range with non-negative values for . Both stable and unstable fixed points will be found. Another numerical procedure to find the thermodynamic properties of the model is to iterate the recursion relations Eqs. (3) directly, this will lead only to the stable fixed points. To study the stability of the fixed points, it is useful to obtain the jacobian of the recursion relations. The elements of the jacobian matrix calculated at the fixed point, , are:

(8a) | |||||

(8b) | |||||

(8c) |

The stability limit of the fixed point may then be found requiring the dominant eigenvalue of this to have a unitary modulus. The jacobian for the limit cycle will be the product of two jacobian matrices defined above, calculated at the pair of variables at the limit cycle, so that .

The grand-canonical partition function of the model on the Cayley tree is obtained considering the operation of attaching subtrees to the central site of the tree. If the central site is in sublattice A, this leads to the following expression:

(9) |

and a similar expression with the sublattice indexes interchanged is obtained for sublattice B. It is easy then to write down expressions for the densities of sites with small and large particles in the center of the tree. They are:

(10a) | |||

(10b) |

and the densities on sublattice B are obtained permuting the sublattice indexes.

The free energy in the core of the Cayley tree, which corresponds to the Bethe lattice, may be obtained by a generalization of Gujrati’s argument g95 (), which may be found for a particular model in cg03 (). Here we present a simple derivation generalizing the one which is found in osb10 (). On the Cayley tree, we admit that the reduced free energy per site (which corresponds to the grand-canonical free energy divided by and the number of sites in the lattice, which is proportional to its volume) of sites in the ’th generation of the tree will be . For a tree with generations, starting to count at the surface (), we may then write the total free energy as:

(11) |

For a tree with one more generation, the free energy in terms of the free energies per site will be:

(12) |

Now we notice that:

(13) |

In the thermodynamic limit, we expect the free energies per site to reproduce the sublattice structure in the core of the tree, and since sites in consecutive generations belong to different sublattices, we notice that in both possible cases (A or B sublattice at the central site) we reach the conclusion:

(14) |

where the indexes stand for the sublattice. The free energy per site in the core of the tree is therefore given by:

(15) |

Using the Eq. (9) and the fixed point equations (4), after some algebra, we find the result:

(16) |

We notice that this expression is invariant under permutation of the sublattice indexes. For , where the model is solved trivially, the Bethe lattice calculation furnishes the exact solution. The fixed point value in this case is simply , the reduced free energy per site will be , and the density of particles becomes:

(17) |

## Iii Thermodynamic behavior of the model

For simplicity, we will start the study of the thermodynamic properties of the model in the limit , where the density of small particles is very small. In this region, the fluid-solid transition is continuous, and the stability limits of both phases are coincident. To obtain the critical line, we may consider the fixed point equations (4) and require the leading eigenvalue of the jacobian Eq. (8) to be equal to -1, so that:

(18) |

Now we solve these three equations up to first order in , supposing that , , and . We are thus lead to the following values of the expansion coefficients defined above: , , , , and . For the particular case , this solution has been obtained a long time ago b60 (). We also notice that the critical line has a negative initial slope. This may be understood physically realizing that the presence of few small particles does lead to an effective entropic attractive interaction between the large particles, thus favoring their ordering d10 (). At higher values of the slope becomes positive, and finally the transition becomes discontinuous, the critical line meets the coexistence line at a tricritical point.

We may also study the phase diagram in the limit . In the fluid phase, we have the asymptotic fixed point values and . In the solid phase we have , , , and . From these expressions, we may find the behavior of the coexistence line in this limit by requiring the free energies, Eq. (16), of both phases to be equal, and this leads us to the asymptotic behavior for the coexistence line for large values of .

The phase diagram in the plane is shown in figure 3 for a Bethe lattice with , similar diagrams are found for other values of . For larger than the tricritical value, there is an interval of values of where both fixed points are stable, thus characterizing a coexistence of both phases, at a value of for which both free energies are equal. The coexistence line is located between both stability limit lines, as expected. The numerical determination of the coexistence line has to be done carefully, particularly close to the tricritical point, where the range of values of for which both fixed points are stable becomes very narrow. The precise calculation of the localization of the tricritical point also demands some care. It is quite easy to calculate the stability limit of the symmetric fixed point, since the problem may be reduced to finding the solution of an equation in one variable. Once this line is found precisely, we solve the asymmetric fixed point equation on it, starting at a value of larger than the tricritical value. As the value of is lowered, the largest eigenvalue of the jacobian increases and the values of the variables and which solve the fixed point equations (5) become closer, as do the variables and . At the fixed point, the dominant eigenvalue is unitary and . This is illustrated in figures 4, thus leading to an estimated position of the tricritical point. Table (1) presents the locations of the tricritical point for several values of the coordination number , as well as the values of the densities at this point. An alternative procedure for determining the location of the tricritical point will be presented below. In all cases we studied, the critical value of for , given above, is smaller than at the tricritical point, therefore for below and above the minimum critical value of the solution displays a reentrant behavior as is increased, starting in the fluid phase, then getting into the solid phase and finally ending in the fluid phase again, with two continuous transitions.

The phase diagram in the plane defined by the densities of small and large particles is presented in figure 5. Again we may obtain the asymptotic behavior in the limits and using the results for these limits presented above. In the limit of low density of small particles, we find that the densities are given by and , the critical line for shows a linear behavior . In the high density limit, we find the following densities on the coexistence line: for the fluid phase and , so that for this phase we find . Therefore, the line reaches the point , with infinite slope. For the solid phase, we get and , so that . Qualitatively, these features of the phase diagram are similar to the ones in Fig. 3 of the paper by Poland p84 ().

In general, the fixed point or limit cycle may be obtained solving equation (6) for the variable . For a given value of , at small values of we find a single root for the equation in the interval . For values of below the tricritical value, above the critical value of , two roots are found, one of them being a double root. They correspond to the fixed point variables and at sublattices A and B. If , for above the limit of stability of the solid phase and below the limit of stability of the fluid phase, five roots are found, such that the smaller and the larger ones correspond to the fixed points associated to the solid phase, while the intermediate root corresponds to the fluid phase. Between the extremal and the central roots, there are two additional roots, which are unstable. Finally, in the region where only the solid phase is stable, three roots are found, the intermediate one being unstable. These findings lead to an alternative procedure to locate the tricritical point, since there the polynomial , its first and second derivatives vanish. This is discussed in some more detail in the appendix.

Another interesting feature of the Bethe lattice solution of this model is that the isobaric curves of the total density of particles as a function of one of the fugacities show a non-monotonical behavior in the fluid phase. Let us define the total density of particles as , so that it will be in the interval . If we recall that the grand-canonical free energy per site is related to the pressure by:

(19) |

where is the volume per site, so that we may define to be the reduced pressure. For a fixed value of the pressure, we may now obtain the density of particles as a function of the activity . Some resulting curves may be seen in figure 6. We notice that the curves are not monotonic, starting with a negative slope at . The minima in the isobars are located on the dotted lines in the phase diagrams of figures 3 and 6. We notice that this line starts at a particular point of the axis and ends at the tricritical point. Some aspects of these minima may be discussed analytically. To find the point of the axis where the minimum of the isobars is located, we may obtain a solution of the model for , since we can solve it exactly for . The fixed point values of the ratios are, up to linear terms in :

(20a) | |||||

(20b) |

Using these approximate solutions and the Eqs. (10) for the densities, we may then find an approximate expression for the total density of particles:

(21) |

The minimum of the isobars correspond to the condition:

(22) |

This derivative may be calculated at noticing that:

(23) |

and since , the last derivative in the expression above is:

(24) |

Finally, we get:

(25) |

so that in the limit the minimum of the isobars is located at . A similar analysis may be done close to , showing that the slope of the isobars is negative there. As can be seen in Eq. (24), for constant pressure is a decreasing function of , a feature which is valid in the solid phase also. Therefore, any isobar which starts at , will end at such that . This is apparent in figure 6, where we notice that each isobar for finite pressures ends at a finite values of and . Isobars which start in the solid phase, at with , will cross either the critical line or the coexistence line before they end in the fluid phase. In the first case we notice a discontinuity in the slope of the isobars. In the second case, as expected the density is discontinuous as the coexistence line is crossed, so that the minimum is located on the coexistence line.

It is worth mentioning that such non-monotonic behavior of the density of particles for constant pressure as a function of a fieldlike variable is found in nature, one of the most studied examples of this kind is the density anomaly of liquid water, where a maximum is found in the isobaric curves of density as a function of the temperature close to the freezing pointd03 (). In many studies in the literature, simple models were proposed which show such anomalies, and it is believed that an interparticle potential with two length scales may be the source of the thermodynamic anomalies. A recent work of this kind, where also many earlier studies are referenced, may be found in ssob10 (). For such models, the solution on tree-like lattices may also be useful osb10 (), and in this particular example both maxima and minima of the isobaric curves of density as a function of temperature were obtained. Although of course the present model is very different from the lattice gases related to water, it is interesting that here also two length scales are present in the interparticle interactions.

## Iv Conclusion

The Bethe lattice solution of the athermal lattice gas with two kinds of particle, a small one occupying a single site and the other a site and its first neighbors, shows a phase transition between a fluid phase, at low values of the activity of the large particles and a solid phase, which appears at higher values of and where one of the sublattices is preferentially occupied by the large particles. The results of the Bethe approximation for the particular case where no small particles are present () are recovered b60 (). The transition remains continuous for small values of , the critical line starts with a negative slope as is increased, passes through a minimum and ends at a tricritical point, so that for larger than the tricritical value the transition is discontinuous.

The behavior of the densities of particles, as shown in figure 5, may be compared with similar results obtained using series expansions for the model defined on the square lattice, shown in figure 3 in reference p84 () by Poland. Besides the expected quantitative differences, we notice in general a qualitative agreement of both diagrams. A significant difference is that in our calculation the lines corresponding to the fluid and the solid phases meet at the tricritical point forming an angle, while in the diagram by Poland a smooth junction is suggested. Since the dotted lines in Poland’s diagram are the result of an extrapolation, it seems that this aspect may be worth to be studied in more detail using other techniques. However, it may be possible that the angle we find here is a characteristic of the mean field approximation implicit in Bethe lattice calculations.

## Acknowledgments

We thank Prof. Ronald Dickman for calling our attention to the model studied here, for helpful discussions, and a critical reading of the manuscript. JFS is grateful for partial financial support by the brazilian agency CNPq.

## Appendix A Determination of the tricritical point

As mentioned in the text, one way to determine the tricritical point is to solve the set of three equations for the polynomial defined in Eq. (6):

(26a) | |||||

(26b) | |||||

(26c) |

Although it seems to be a rather simple task to solve this system of nonlinear algebraic equations for , , and , standard numerical methods, based on Newton-Raphson procedures, often do not converge to the expected (physical) solution. This may be due to the fact that the first and second derivatives of with respect to vanish at the solution. Therefore, we used another procedure, taking advantage of the fact that is a polynomial in the variable :

(27) |

where . Now , the value of at the tricritical point, is a triple root of the polynomial, so that we may write:

(28) |

Comparing terms with the same powers of in the equation above, we may solve for the coefficients in terms of the coefficients and . The result is:

(29) |

for , and

(30) |

The remaining equations, corresponding to the powers , , and , are:

(31a) | |||||

(31b) | |||||

(31c) |

The solution of these equations leads to the activities and the value of the variable at the tricritical point.

## References

- (1) J. D. van der Waals, On the continuity of the Gaseous and Liquid States, Dover (2004).
- (2) J. P. Hansen and I. R. McDonald Theory of Simple Liquids, Academic Press (1986).
- (3) N. Metropolis et al, J. Chem Phys. 21, 1087 (1953).
- (4) H. C. M. Fernandes, J. J. Arenzon, and Y. Levin, J. Chem. Phys. 126, 114508 (2007) and references therein.
- (5) C. Domb, Nuovo Cim. Suppl., Series X, 9, 26 (1958).
- (6) D. M. Burley, Proc. Phys. Soc. 75, 262 (1960) and 77, 451 (1961).
- (7) W. Guo and H. W. J. Blöte, Phys. Rev E 66, 046140 (2002).
- (8) A. G. Cunha-Neto and R. Dickman, Comp. Phys. Commun. 182, 719 (2011).
- (9) R. J. Baxter, J. Phys. A 13, L61 (1980); Exactly Solved Models in Statistical Mechanics, Academic Press, London (1982).
- (10) D. Frenkel and A. L. Louis, Phys. Rev. Lett. 68, 3363 (1992).
- (11) D. Poland, J. Chem. Phys. 80, 2767 (1984).
- (12) P. D. Gujrati, Phys. Rev. Lett. 74, 809 (1995).
- (13) A. Corsi and P. D. Gujrati, Phys. Rev. E 68, 031502 (2003).
- (14) T. J. Oliveira, J. F. Stilck, and M. A. A. Barbosa, Phys. Rev. E 82, 051131 (2010).
- (15) R. Dickman, private communication (2010). This effect is also discussed by Poland p84 ().
- (16) P. G. Debenedetti, J. Phys.: Condens. Matter 15, R1669 (2003).
- (17) J. N. da Silva et al, J. Chem. Phys. 133, 244506 (2010).

3 | 2.74598 | 11.9231 | 0.3937(8) | 0.1476(5) |
---|---|---|---|---|

4 | 1.16956 | 3.02938 | 0.2985(7) | 0.1117(5) |

5 | 0.734355 | 1.539662 | 0.2407(3) | 0.0894(3) |

6 | 0.533384 | 0.995756 | 0.2013(3) | 0.0746(3) |

7 | 0.418259 | 0.725526 | 0.1729(2) | 0.0641(3) |

8 | 0.343861 | 0.566972 | 0.1515(2) | 0.0561(2) |

9 | 0.291885 | 0.463718 | 0.1349(2) | 0.0499(2) |

10 | 0.253531 | 0.391509 | 0.1215(2) | 0.0450(2) |