Maxym Dudka
email: maxdudka@icmp.lviv.ua

Olivier Bénichou
email: benichou@lptmc.jussieu.fr

Gleb Oshanin
email: oshanin@lptmc.jussieu.fr

Order-disorder transitions in lattice gases with annealed reactive constraints

Keywords: catalytically-activated reactions, lattice gases with athermal interactions, annealed constraints, order-disorder transitions, Bethe lattice

We study equilibrium properties of catalytically-activated A + A → ⊘ reactions taking place on a lattice of adsorption sites. The particles undergo continuous exchanges with a reservoir maintained at a constant chemical potential µ and react when they appear at the neighbouring sites, provided that some reactive conditions are fulfilled. We model the latter in two different ways: In the Model I some fraction p of the bonds connecting neighbouring sites possesses special catalytic properties such that any two As appearing on the sites connected by such a bond instantaneously react and desorb. In the Model II some fraction p of the adsorption sites possesses such properties and neighbouring particles react if at least one of them resides on a catalytic site. For the case of annealed disorder in the distribution of the catalyst, which is tantamount to the situation when the reaction may take place at any point on the lattice but happens with a finite probability p, we provide an exact solution for both models for the interior of an infinitely large Cayley tree -the so-called Bethe lattice. We show that both models exhibit a rich critical behaviour: For the annealed Model I it is characterised by a transition into an ordered state and a re-entrant transition into a disordered phase, which both are continuous. For the annealed Model II, which represents a rather exotic model of statistical mechanics in which interactions of any particle with its environment have a peculiar Boolean form, the transition to an ordered state is always continuous, while the re-entrant transition into the disordered phase may be either continuous or discontinuous, depending on the value of p.

Introduction

Catalytically activated reactions involve particles that react only in the presence of another agent -a catalyst -and remain chemically inactive otherwise. Such processes are widespread in nature and are also involved in a variety of technological and industrial applications [START_REF] Bond | Heterogeneous Catalysis: Principles and Applications[END_REF][START_REF] Avnir | Fractals in Science[END_REF]. Different kinetic, equilibrium and out-of-equilibrium properties of such reactions, as well as various theoretical concepts, available analytical and numerical approaches have been comprehensively reviewed in Ref. [START_REF] Liu | [END_REF].

The classical textbook approach to the kinetics of such processes focuses only on the average concentrations of the species involved. In consequence, one usually writes systems of differential rate equations of varying complexity with multiple parameters, as prescribed by the formal-kinetic "law of mass action" (see, e.g., Ref. [START_REF] Liu | [END_REF]). Ziff, Gulari and Barshad (ZGB) [4] and subsequently, Ziff and Fichthorn [5], Fichthorn, Gulari and Ziff [6], were apparently first to realise that fluctuations in adsorption/desorption events, fluctuations of coverages and spatial correlations can cause a severe departure from the deterministic descriptions based on the rate equations approach. Using computer simulations ZGB studied the so-called monomer-dimer model, introduced as an idealised description of the important process of CO oxidation on a catalytic surface, which revealed an emerging spectacular cooperative behaviour: it was found that the monolayers formed by the molecules adsorbed on the surface may undergo discontinuous or continuous phase transitions in different parameter ranges.

More specifically, ZGB considered a reaction occurring by the following three steps: a) an irreversible adsorption of CO molecules from their reservoir onto a single-crystal catalyst surface, modelled as a simple square lattice of adsorption sites, b) an irreversible adsorption of O 2 molecules onto the lattice from their reservoir with the subsequent dissociation of O 2 into two O atoms, each residing on a separate lattice site (while CO requires only a single site, which explains why the model is called the "monomerdimer" one) and c) an instantaneous (perfect) irreversible reaction between neighbouring absorbed CO and O, followed by an immediate desorption of the reaction product CO 2 , which does not interact further with the system. It was realised that on a twodimensional square lattice upon lowering the CO adsorption rate the system undergoes a first-order (discontinuous) phase transition from a CO saturated inactive phase into a reactive steady state, followed by a continuous transition into an O 2 -saturated inactive phase, which belongs to the same universality class as directed percolation and the Reggeon field theory [7]. While such a model evidently discards many important features of real physical systems, such as, e.g., molecular diffusion on the surface, reversibility of the adsorption of the species involved, imperfect reaction, when two neighbouring adsorbed reactants react only with a finite probability, and etc., it does reproduce two types of phase transitions which occur in realistic systems [START_REF] Liu | [END_REF].

A more simple, monomer-monomer model was introduced by Ziff and Fichthorn [5] and was subsequently studied in Refs. [8][9][10]. In this model one focuses on a general reaction scheme A + B → ⊘ where A and B are two chemically different species but both are monomers, in the sense that both require just a single lattice site for an irreversible adsorption. It was realised that the monomer-monomer model exhibits a first-order transition from a phase saturated with one of the species to the phase saturated with the other one. Allowing for an adsorption of one of the species, i.e., connecting the system to a reservoir, leads to a continuous transition which also belongs to the directed percolation universality class [10]. We note also that the fluctuationinduced behaviour of the catalytically-activated reactions between molecules with more complicated structures has been studied [11][12][13] and also the influence of the lateral diffusion of the adsorbed species on the reaction kinetics (see, e.g. Refs. [14][START_REF] Marro | Nonequilibrium Phase Transitions in Lattice Models[END_REF][START_REF] Argyrakis | [END_REF][17][18] and references therein) as well as of the heterogeneity in the spatial distribution of the catalyst (see, e.g., Refs. [19,20]) have been addressed.

As long as one is interested only in the thermodynamic equilibrium properties of the monolayers of the adsorbed molecules, formed in the course of such catalyticallyactivated reactions, one may realise that the latter should be very similar to the analogous properties of the adsorbates emerging in the models of hard-core lattice gases with athermal interactions, in which one studies a reversible deposition of "hard", nonoverlapping objects on two-dimensional lattices. Indeed, the constraint that any two particles appearing at the neighbouring sites have an infinite repulsion, (so that the input of such configurations into the partition function of the adsorbate is zero), and the condition, in case of catalytically-activated reactions, that once two particles appear at the adjacent adsorption sites they react and leave the system, play essentially the same role. A prominent example of such lattice gases, solved exactly by Baxter [21,[START_REF] Baxter | Exactly solved models in statistical mechanics[END_REF], is furnished by the so-called hard hexagon model, which is a two-dimensional lattice model of a gas, where identical particles are allowed to be on the vertices of a triangular lattice but no two particles may be adjacent. From the perspective of catalytically-activated reactions, this model can be interpreted as a simple reaction A + A → ⊘ between monomers A, being at contact with a reservoir maintained at a constant chemical potential µ, adsorbing onto empty lattice sites and desorbing back to the reservoir, and undergoing a chemical reaction as soon as any two A particles appear at the adjacent sites. Much progress has been made within the recent years in understanding equilibrium properties of such lattice gases, involving similar or different particles, which may also have a different shape (see, e.g., Refs. [START_REF] Bouttier | [END_REF]24]), revealing a rich critical behaviour, characterised by continuous and discontinuous phase transitions which are quite analogous to the ones observed in the irreversible ZGB model. We note also that a similarity between a simple A + A → ⊘ reaction and random sequential adsorption has been pointed out in Ref. [25], since the removal of two nearest-neighbour reactants is equivalent to a "deposition" of a "dimer" of two empty sites. Such a duality has been used in Ref. [25] to derive an exact solution of the model on the Bethe lattice.

In this paper we study equilibrium properties of an adsorbate formed in the course of the reaction

A + A → ⊘ (1)
between monomers A, undergoing continuous exchanges with their vapour phase -a reservoir maintained at a constant chemical potential µ, and taking place on the lattice with annealed reactive properties. Our goal here is to show that already this simplest possible reaction scheme involving only monomers of the same type exhibits quite a rich critical behaviour, depending moreover on the specific way how the reaction is modelled and also on the value of the particles' reactivity. The critical behaviour which we predict is characterised by a transition into a phase with a broken symmetry (ordered or an alternating state, in which layers with low and high densities are formed), which is continuous with a finite jump in compressibility. This transition is followed by, upon an increase of the chemical potential µ by a re-entrant transition into a disordered phase, which (depending on the reaction model) may be continuous with a finite jump in compressibility, or even discontinuous with a finite jump of density.

We will consider here two different ways of modelling the catalytic reaction: In the Model I, we stipulate that some fraction p of the bonds < ij > connecting nearest-neighbouring sites i and j possesses catalytic properties so that the A particles react when two of them appear on the neighbouring sites connected by such a bond. In the Model II, we suppose that not the bonds but rather some lattice sites themselves are the catalytic agents such that any two neighbouring As will react instantaneously and leave the system, when at least one of them resides on a special catalytic sites.

More specifically, in the Model I we associate with each bond < ij > a random variable ζ <ij> , such that ζ <ij> = 1, if < ij > is a catalytic bond (which event is chosen with probability p, independently of other bonds) and ζ <ij> = 0, otherwise, (with probability 1 -p). When any two A-s appear at the neighbouring sites i and j connected by a catalytic bond, they instantaneously react (and the product leaves the system). The A particles harmlessly coexist on neighbouring sites connected by noncatalytic bonds. For such a model, the corresponding grand canonical partition function of the adsorbate with a fixed distribution of the catalytic bonds can be written as

Z ({ζ <ij> }) = {n j } exp βµ j n j <ij> 1 -ζ <ij> n i n j , ( 2
)
where µ is the chemical potential, β is the inverse temperature, (measured in units of the Boltzmann constant k B ), n j is the Boolean variable describing the occupation of a given site j (n j = 1 for an occupied site and n j = 0, otherwise), the sum with the subscript {n j } denotes summation with respect to the states of all occupation variables of all sites, and lastly, the product sign with the subscript < ij > signifies that the product is taken over all bonds of the embedding lattice. For p ≡ 1, Z ({ζ <ij> }) in eq. 2 evidentlycoincides with the grand partition function of the gas of hard molecules with infinite nearest-neighbour repulsion [21,[START_REF] Baxter | Exactly solved models in statistical mechanics[END_REF].

In the Model II, we assign to each of the sites a random quenched variable χ j , which, in a similar fashion, is equal to 1 for a "catalytic" site (probability p, independently of other sites) and zero, otherwise, (with probability 1 -p). In this case, two A particles instantaneously react A + A → ⊘, (and the product ⊘ leaves the system), if at least one of them is on a catalytic site. If neither of them occupies a catalytic site, the particles do not react. For a given distribution of the catalytic sites, the grand canonical partition function Z ({χ i }) reads

Z ({χ i }) = {n j } exp βµ j n j i i,<ij> 1 -χ i n i n j , (3)
where the product with the subscript i runs over all sites of the lattice, while the product sign with the double subscript i, < ij > signifies that the product operation is taken over all the sites j neighbouring to the site i. Note that the expression in eq. 3 becomes identical to the one in eq. 2 for p ≡ 1, as it should be. We will focus in what follows on the annealed limit in the distribution of the catalytic bonds or sites. Such a limit is interesting in its own right since it corresponds to a physical situation in which the neighbouring adsorbed A particles may react at any place on the lattice, but only with a finite reaction probability p, so that the reaction between them is not perfect. The properties of the adsorbate in this limit will be described by the annealed grand canonical partition functions, i.e., averaged directly over the distributions of the random variables ζ <ij> and χ j . The annealed grand canonical partition function for the Model I reads

Z (bonds) (p) = {n j } exp βµ j n j (1 -p) <ij> n i n j , (4)
and describes a grand canonical partition function of a lattice gas with purely repulsive, soft interactions. Exact solution of the one-dimensional Model I for arbitrary p for both the cases of random quenched and of the annealed distribution of the catalytic bonds has been obtained in Refs. [26,27].

In turn, for the annealed Model II we have the following grand canonical partition function

Z (sites) (p) = {n j } exp βµ i n i (1 -p) i n i Ψ i , (5)
where Ψ i is a Boolean function of the form

Ψ i = 1 - i,<ij>
(1 -n j ) ≡ 0 when all n j = 0, 1 when at least one n j = 1.

Therefore, in the annealed Model II we again deal with a lattice gas with purely repulsive, soft interactions between the particles, but here the interactions are very peculiar: the total interaction energy of the particle occupying site i is not proportional to the number of neighbouring particles, as it usually happens, but has a Boolean form such that Ψ i = 0 if all the sites neighbouring to i are vacant, and Ψ i = 1 independently of the actual number of occupied neighbouring sites if at least one of them is occupied. Note that in this rather exotic model of statistical mechanics the interactions are therefore less restrictive than those appearing in the Model I, in which each pair of neighbouring particles contributes to the energy of the system. Exact solution of the Model II in one-dimensional systems was obtained in Refs. [28,29] for arbitrary p, for both the cases of random quenched and of the annealed distributions of the catalytic sites.

Analysing the critical behaviour of the annealed versions of the Model I and Model II, we will resort here to a mean-field-like approximation, solving both models exactly on the Bethe lattice -a deep interior of the so-called Cayley tree (see Fig. 1, left panel) well away of the boundary sites. This is a topologically simpler object than usual lattices, because such trees do not have closed loops (such that the approximation is equivalent to the standard Bethe-Peierls theory [START_REF] Baxter | Exactly solved models in statistical mechanics[END_REF]). We focus solely on the case of a Bethe lattice with the coordination number equal to three, which means that we provide here an approximate solution for the annealed Models I and II on a honeycomb lattice. We note that for many complicated lattice models of statistical mechanics such an approximation has been invoked first, since it is often amenable to a completely analytical analysis reducing the original many-body problem to the analysis of the limit solutions of some non-linear recursion schemes, in which a phase transition manifests itself by the spontaneous break up of the symmetry between recursion terms of an odd and an even order. In many cases, a detailed studies of the properties of the original lattice models can be carried out, showing no sign of a pathological behaviour [30]. To name but a few, we mention studies of the phase diagrams of athermal lattice gases [24,31,32], modulated phases of the Ising model with competing interactions [33][34][35], Potts models [36], lattice models of glassy systems [37,38], different aspects of the localisation transition [39,40] and phase diagram for ionic liquids in non-polarised nano-confinement [41].

Our main results are as follows: For the annealed Model I defined on the Bethe lattice with the coordination number three, we set out to show that there exists a) a critical value p = p (bonds) c = 8/9, such that for p < p (bonds) c no broken symmetry phase emerges, while for p (bonds) c < p < 1, upon increasing the chemical potential µ, one observes b) a transition from a disordered to an ordered (alternating) phase at µ = µ (bonds) c,1 (p). In this ordered phase the tree spontaneously partitions into alternating layers -the ones which are densely populated by the A particles and the ones which are almost completely devoid of them. For the original honeycomb lattice it will mean that (somewhere close to µ = µ (p) → ∞. d) We also calculate the mean density in the disordered and ordered phases, the staggered density in the ordered phase which serves as the order parameter for the system under study and the compressibility. For the annealed version of the Model II defined on the Bethe lattice with the coordination number three, we observe a critical behaviour, which turns out to be even somewhat richer than that of the Model I. We show that a) there exist two critical values of the reaction probability p: p , upon a gradual increase of µ, the system undergoes a continuous transition from a disordered into an ordered, alternating phase at µ = µ (p) the system undergoes a re-entrant continuous transition into a disordered phase, while for the latter this transition is discontinuous, with a finite jump of both the mean and the staggered densities. d) We also determine the mean density in the disordered and the ordered phases, the staggered density and the compressibility.

The paper is outlined as follows: In Sec. 2 we describe the system's geometry and present the details of the derivation of the annealed grand canonical partition functions (in what follows, we will call them for brevity just as "partition functions") for the Models I and II. In Sec. 3 we focus on the singular case p ≡ 1, in which both Models (as well as their annealed versions) are identical, and briefly recall the classical analysis due to Runnels of the critical behaviour of a gas of hard molecules on the Bethe lattice with the coordination number equal to three. This permits us to set up a unifying framework for what follows. In Sec. 4 we study the annealed version of the Model I for an arbitrary value of the reaction probability p, 0 < p ≤ 1. We first derive the recursion relations obeyed by the partition function in eq. 4 on a Bethe lattice and analyse their critical behaviour, as manifested by the breaking of the symmetry between the terms of odd and even order. Next, in Sec. 5, we focus on the annealed version of the Model II. We derive here the recursion relations obeyed by the annealed partition function defined in eqs. 5 and 6, and analyse their critical behaviour, and also determine the mean density in the disordered and ordered phases, the staggered density and the compressibility, as the functions of the chemical potential. Finally, we conclude in Sec. 6 with a brief recapitulation of our results.

Partition functions in the annealed limit

Consider a Cayley tree with the coordination number equal to three, (see Fig. 1, left panel), having N generations of sites (such that the "volume" M = 3(2 N -1) + 1 for N ≥ 0 ), in contact with the vapour phase of hard particles A. The A particles can adsorb onto the vacant sites and can desorb back to the reservoir. The reservoir is characterised by a chemical potential maintained at a constant value µ and measured relative to the binding energy of an occupied site, such that µ > 0 indicates a preference for adsorption. Correspondingly, the activity z is defined as z = exp(βµ). Hard-core interactions prohibit double occupancy of any node j and the occupation of this node is described by a Boolean variable n j such that n j = 1 for an occupied node and n j = 0, otherwise.

We will consider two ways of modelling the catalytic reaction between the neighbouring A particles, which differ in the definition of the catalytic agent. In Model I, we suppose that these are some bonds connecting neighbouring sites, which are deemed to have special catalytic properties promoting the reaction, while in Model II these will be some special catalytic sites. In the former case, only the A particles separated by a catalytic bond react, while in the latter, in order for the reaction to take place, at least one of two neighbouring As has to reside on a special catalytic site. In both cases, a given bond or a given site can be in a catalytic state with probability p, independently of the environment. The partition functions of the adsorbate corresponding to the Models I and II, Z ({ζ <ij> }) and Z ({χ i }), are defined in eqs. 2 and 3, respectively.

Here we will be concerned with the annealed versions of both Models whose partition functions are obtained by directly averaging Z ({ζ <ij> }) and Z ({χ i }) over the distributions of the catalytic bonds and sites.

The averaging of Z ({ζ <ij> }) in eq. 2 is very straightforward

Z (bonds) (p) = Z ({ζ <ij> }) ζ <ij> = {n j } z j n j <ij> 1 -ζ <ij> n i n j ζ <ij> = {n j } z j n j <ij> p (1 -n i n j ) + 1 -p , = {n j } z j n j (1 -p) <ij> n i n j , (7)
from which one can directly read off the result in eq. 4. Averaging of Z ({χ i }) in eq. 3

Order-disorder transitions in reactive lattice gases M. Dudka et al is only slightly more involved. Here we have

Z (sites) (p) = Z ({χ i }) χ i = {n j } z j n j i i,<ij> 1 -χ i n i n j χ i = {n j } z j n j i p i,<ij> (1 -n i n j ) + 1 -p . (8)
Further on, noticing that

p i,<ij> (1 -n i n j ) + 1 -p ≡ p i,<ij> (1 -n j ) + 1 -p n i , (9)
and next that, trivially,

p i,<ij> (1 -n j ) + 1 -p ≡ (1 -p) 1-i,<ij> (1-n j ) , (10)
we find that the disorder-average partition function Z (sites) (p) attains the form defined in our eqs. 5 and 6. Parenthetically, we note that the expressions in eqs. 4 and 5 are valid for any lattice, not necessarily for the Cayley tree. Next, focussing on the limit N → ∞, we will discard the surface effects considering the behaviour in the interior of the Cayley tree. In other words, we will turn to the thermodynamic limit of the Cayley tree, i.e., to the behaviour on the so-called Bethe lattice. Under such simplifying assumptions, the thermodynamics of the system is given by the annealed pressure (in units of the lattice cell area), defined as

P I (T, µ) = 1 β lim M →∞ ln Z (bonds) (p) M , (11)
for the Model I, and as

P II (T, µ) = 1 β lim M →∞ ln Z (sites) (p) M , (12)
for the Model II. Once P -s are known, all other thermodynamic quantities of interest can be obtained by differentiating P with respect to µ or T .

3. Singular case p ≡ 1. Runnels' analysis.

To set up the scene, we first consider the special singular case p ≡ 1, in which all the bonds and all the sites are catalytic, implying that the reactive constraint is imposed everywhere. Analysis in this case is much more simple than for an arbitrary p but will permit us to introduce all the basic concepts and to establish a unifying framework for what follows. In this case, eqs. 2 and 3 become identical, and also identical to the partition function of the lattice gas of hard particles, having an infinite repulsion between the neighbouring ones:

Z({1}) = {n j } z j n j <ij> 1 -n i n j . ( 13
)
Below we briefly recall some classical results due to Runnels [31] obtained for Z({1}) defined on the Bethe lattice with the coordination number three, using for further convenience a bit different settings and notations. The partition function Z({1}) of the entire Cayley tree writes:

Z({1}) = Z (0) N ({1}) + Z (1)
N ({1}) , (14)
where

Z (0) N ({1}) and Z (1)
N ({1} are the partition functions of the entire tree with a vacant and an occupied central site, respectively.

Define next a subtree with N generations (see Fig. 1, right panel), emanating from the central site O, and introduce two auxiliary partition functions -B N (1) and B N (0), where the former denotes the partition function of a subtree with an occupied central site, while for the latter the central site of the subtree is vacant. Clearly, we have

Z (0) N ({1}) = B 3 N (0) , (15)
because the entire Cayley tree with the vacant central site decomposes into three independent subtrees with a vacant central site, while Z

N ({1}) decomposes into three independent subtrees with an occupied central site :

Z (1) N ({1}) = z -2 B 3 N (1) , (16)
where the factor z -2 prevents the over-counting of the contribution of the central site.

We seek next the recursion relations obeyed by B N (1) and B N (0). For p ≡ 1, if the central site is occupied, the neighbouring one (generation 1) is always vacant, such that

B N (1) = z B 2 N -1 (0) , (17)
where the factor z stems from the contribution of the occupied central site. Next, if the central site is vacant, the neighbouring one (generation 1) can be either vacant or occupied, such that B N (0) can be represented via the partition functions of the subtrees with N -1 generations as

B N (0) = B 2 N -1 (0) + z -1 B 2 N -1 (1) , (18)
where the factor z -1 in the second term on the r-h-s of eq. 18 prevents the over-counting of the contribution of the occupied site, neighbouring the central one, from which two subtrees with N -1 generations emanate. Together with the evident "initial" conditions B 0 (1) = z and B 0 (0) = 1, eqs. 17 and 18 totally define B N (1) and B N (0) (and hence, Z({1})), for the Cayley tree with an arbitrary number of generations N. Equations 17 and 18 can be simplified by introducing an auxiliary parameter

x N (z) = B N (1) zB N (0) . ( 19
)
This parameter, multiplied by z, is simply the ratio of the number of the subtrees with an occupied central site and the number of the subtrees with a vacant central site. Dividing eq. 17 by eq. 18, we find that x N (z) obeys the following recursion :

x N (z) = 1 1 + z x 2 N -1 (z) , x 0 (z) = 1 . (20)
Once x N (z)-s are determined, one can readily find B N (0) via the recursion

B N (0) = 1 + zx 2 N -1 (z) B 2 N -1 (0) = B 2 N -1 (0) x N (z) , B 0 (0) = 1 , (21)
which follows from eq. 18 and the definitions of x N (z) in eqs. 19 and 20, and gives

B N (0) = N j=0 1 x 2 N-j j (z) . ( 22
)
In turn, the partition function of the entire Cayley tree in eq. 14 can be straightforwardly written using x N (z)-s as

Z N ({1}) = 1 + zx 3 N (z) N j=0 1 x 2 N-j j (z) 3 , (23)
so that the pressure P in the limit N → ∞ (and hence, M → ∞) will be given by

P (T, µ) = 1 β lim N →∞ 3 2 N M N j=0 2 -j ln 1 x j (z) + ln (1 + zx 3 N (z)) M = 1 β ∞ j=0 2 -j ln 1 x j (z) . (24)
Note that eq. 24 defines the pressure of the adsorbate on the entire Cayley tree, such that to get an analogous result for the Bethe lattice we have to subtract the contribution due to the boundary sites. For the case p ≡ 1, the corresponding procedure will be explained below. For the general case 0 < p ≤ 1 we will adapt the procedure elaborated in Ref. [42], which is explained in Appendix B.

Solutions of the recursion relations for p ≡ 1

The central question now is the behaviour of x N (z)-s, which are all functions of the activity z. Before we turn to analytical calculations, it might be expedient to numerically generate several consecutive terms x N (z), in order to get a clue to their actual behaviour.

In Fig. 2 we present first 100 terms x N (z) as functions of the activity z = exp(βµ).

One observes that for a relatively small βµ (or z), as N grows, all x N (z) converge to the same curve x(z). Therefore, for sufficiently small z one has a sequence convergence as N → ∞ to the unique limit x(z), which, in virtue of eq. 20, obeys

z x 3 (z) + x(z) -1 = 0 . ( 25
)
The discriminant △ = -27z 2 -4z of the cubic eq. 25 is strictly negative, meaning that eq. 25 has two non-real complex conjugate roots and one real root x(z) which is given explicitly by Two dashed blue curves starting at βµ = 2 ln(2) are the roots of the second term in eq. 29 defining two distinct limit curves x odd (z) and x even (z). The figure illustrates continuous phase transition from disordered state described by solution of eq. 25 (black line) to symmetry broken phase with alternating layers of different densities described by solution of eq. 29 (dashed blue lines). The phase transition occurs at z c = 4 at which all three curves meet.

x(z) = - 1 3z C + △ 0 C , (26
where

△ 0 = -3z , C = △ 1 + 3z √ -3△ 2 1/3 , △ 1 = -27z 2 . ( 27
)
This root is depicted in Fig. 2 by a thick black curve. Further on, Fig. 2 shows that there is an apparent "critical" value z c such that for z > z c the odd terms converge to one function, x 2N +1 (z) → x odd (z), while even terms converge to another one, x 2N (z) → x even (z). In this case of the so-called subsequence convergence, one may calculate x odd (z) and x even (z) by iterating eq. 20 once again, so that the the resulting recursion scheme will involve the terms of the same parity only; that is, we rewrite eq. 20 formally as

x N (z) = 1 + z 1 + z x 2 N -2 (z) 2 -1 . ( 28
)
Turning next to the limit N → ∞, we find that two distinct limit functions x odd (z) and x even (z) must satisfy the fifth-order equation, that conveniently factors into

z x 3 (z) + x(z) -1 zx 2 (z) -zx(z) + 1 = 0 . ( 29
)
One notices that the first term is precisely our previous eq. 25, which has a real root defined by eq. 26, while the second term has two roots:

x even (z) = z + √ z 2 -4z 2z , (30)
and

x odd (z) = z - √ z 2 -4z 2z , (31)
which we depict in Fig. 2 by blue dashed curves. These two roots are real only for z ≥ 4. Note now that the three curves defined by eqs. 26, 30 and 31, all meet at z = 4 each assuming the value 1/2. Since x odd (z) and x even (z) are complex-valued for z < 4, the unique limit x(z) in eq. 26 must obtain for z < 4. As shown in Ref. [31], the result in eq. 26 becomes unstable for z > 4. This signifies that the alternating limits x odd (z) and x even (z) obtain for z > 4.

One may notice that the special point z = 4, where the singularity occurs, can be obtained from the direct stability analysis of the fixed point solution of recursion 20. This recursion is stable if the slope of the function f (x N -1 (z)) in the right-handside of eq. 20 in the vicinity of the fixed point is smaller than unity; that being

|f ′ (x N -1 (z))| x N-1 (z)=x(z) < 1.
Therefore, solving equation:

f ′ (x N -1 (z))| x N-1 (z)=x(z) = -1 (32)
with x(z) defined by eq. 26 one gets z = 4 as the point where the solution 26 looses its stability and a bifurcation occurs. For z > 4 the limit cycle defined by the second term in eq. 29 is a stable solution. Therefore, z = 4 is the only candidate for a singularity or a phase transition. As shown by Runnels [31] (see also Ref. [32] for a more detailed discussion of the behaviour on the entire Cayley tree), in actual fact it is a rather subtle issue. The point is that on the Cayley tree the various derivatives of x j become indeed unbounded with increasing j, but not rapidly enough to spoil the convergence of eq. 24 and its derivatives, due to the presence of the factor 2 -j . The ensuing smooth overall behaviour attributable to this factor is easily understood : most of the sites of the Cayley tree are near its exterior surface, where they are essentially independent; the highly correlated sites are deep in the interior of the tree. In other words, zeros of the partition function do close in on z = 4 but with vanishing density, which allows the system to sneak through on the real axis with no transition. Conversely, if one discards the influence of the overwhelming majority of the surface sites and focuses on the behaviour of the interior of the Cayley tree far away the surface -the so-called Bethe lattice, one will indeed observe a transition at z = z c = 4. The transition occurring in this interior region far removed from the exterior surface is continuous, with a finite jump in the compressibility [31]. This is a transition of the disorder-order type so that at z = z c = 4 the interior region spontaneously partitions in alternating, highly occupied and almost devoid of particles layers. Note that these predictions on the existence and on the type of the phase transition obtained for the Bethe lattice, which is essentially a meanfield (albeit quite reliable [30]) approximation, appear to be qualitatively correct for models defined on corresponding (with the same coordination number) regular lattices, as evidenced later by the exact solution of the so-called hard-hexagons model [21].

Order parameter for p ≡ 1

We briefly recall next Runnel's calculations of the particle densities in the alternating layers. Differentiating eq. 24, one gets for the mean density

ρ(z) = βz ∂P ∂z = ∞ j=1 2 -j t j (z) , (33)
where t j (z) = -z∂ ln (x j (z)) /∂z satisfies the recursion of the form

t j (z) = b j (z) (1 -2t j-1 (z)) , t 0 (z) = 1 , (34)
in terms of b j (z) = 1 -x j (z). Straightforward calculations give [31]

t j (z) = j-1 i=0 (-2) j j n=j-i b n (z) , (35)
which, being inserted in eq. 33, allows to rewrite the latter equation as

ρ(z) = ∞ j=0 2 -j r j (z) , r j (z) = ∞ i=0 (-1) i j+i n=j b n (z) . (36)
Note now that b j (z) and r j (z) have simple interpretations [31]: From the preceding equations, it is clear that b j (z) is the fraction of the activity-weighted configurations of subtrees having an occupied root, or in other words, the average occupancy of the root at activity z. Noticing then that eq. 36 can be alternatively written as

r j (z) = b j (z) (1 -r j+1 (z)) , (37)
(which can be straightforwardly reiterated to give back eq. 36), one concludes that r j (z) is simply the average occupancy of a site at generation j. Indeed, eq. 37 can be interpreted as the definition of the average occupancy r j (z): it equals the factor b j (z), that being, the weighted fraction of configurations of a subtree with the selected site occupied, times the factor (1 -r j+1 (z)), which defines the fraction of configurations of a subtree with the inner site adjacent to the selected site vacant to permit occupancy of the selected site. Therefore, eq. 36 shows that the overall density is a weighted average of the densities of various generations, the weighting factor 2 -j being proportional to the number of sites in generation j.

One may now define the order parameter as the staggered density δρ = |r j (z)r j+1 (z)|, i.e., the difference of average occupations of adjacent layers deep in the interior of the tree in the limit j → ∞. Clearly, for z < z c all b j (z) converge to the single limit b(z) = 1 -x(z) regardless of the parity, where x(z) is the solution of eq. 25. In this case, in virtue of eq. 37, the average occupations of the layers converge as j → ∞ to the unique limit r j (z) → r(z) = b(z)/(1 + b(z)) such that δρ ≡ 0.

For z > z c the situation is different. As j → ∞, the parameter b j with odd j converges to b + (z) = 1 -x odd (z), eq. 31, while the parameter b j with even j converges to b -(z) = 1 -x even (z), eq. 30. Consequently, one finds from eq. 36 that the densities r j (z) in the layers with odd j converge to r + (z), while the average occupations of the layers with even j converge to r -(z), which are given explicitly by

r + (z) = b + (z)(1 -b -(z)) 1 -b + (z)b -(z) , r -(z) = b -(z)(1 -b + (z)) 1 -b + (z)b -(z) . (38)
In turn, for z > z c = 4 the order parameter obeys

δρ = b + (z) -b -(z) 1 -b + (z)b -(z) = x even (z) -x odd (z) x even (z) + x odd (z) -x even (z)x odd (z) . ( 39
)
and equals zero for z < z c . In the limit z → ∞, r + (z) → 1, r -(z) → 0, and δρ → 1.

Model I : Catalytic bonds

In this Section we analyse the critical behaviour of the annealed version of the Model I, whose partition function is defined by eq. 4. Here, similarly to the previously considered case p ≡ 1, Z (bonds) (p) can be formally decomposed as Z (bonds) (p) = Z

(p) = B 3 N (0, p) , Z (bonds,1) N
(p) = z -2 B 3 N (1, p) , (40)
where B N (0, p) and B N (1, p) are the partition functions of a subtree (see, Fig. 1, right panel) with a vacant and an occupied root, respectively. Further on, we find straightforwardly that B N (0, p) and B N (1, p) obey the following recursions :

B N (1, p) = zB 2 N -1 (0, p) + z(1 -p) B 2 N -1 (1, p) z , (41)
and

B N (0, p) = B 2 N -1 (0, p) + z -1 B 2 N -1 (1, p) . (42)
While the latter expression coincides with eq. 18 of the previous section, the former one, given by eq. 41, has a different form compared to that in eq. 17 since here it is possible to have two particles at the neighbouring sites.

The recursions in eqs. 41 and 42 can be simplified by introducing

x N (p, z) = B N (1, p) zB N (0, p) , (43)
which now obeys the recursion :

x N (p, z) = 1 + z(1 -p)x 2 N -1 (p, z) 1 + zx 2 N -1 (p, z) , x 0 (p, z) = 1 . (44)
Note that B N (0, p) still obeys B N (0, p) = (1+zx 2 N -1 (p, z))B 2 N -1 (0, p), like the analogous property in the case p ≡ 1, eq. 21, which implies that the pressure P I (T, µ) of the adsorbate is given by last line of eq. 24 with x j (z) replaced by (x j (p, z) -1 + p)/p.

Solutions of the recursion relations for p < 1

We proceed further by numerically iterating first 20 terms of the recursion defined in eq. 44 for two different values of the parameter p: p = 0.85 and p = 0.95, which we present in panels (a) and (b) of Fig. 3. For p = 0.85 (panel (a)) we observe a sequence convergence x N (p, z) → x(p, z), as N → ∞, to a single limiting curve x(p, z) for both even and odd N for any z. For this single limit case, x(p, z) obeys, in virtue of eq. 44, a cubic equation of the form,

zx 3 (p, z) -z(1 -p)x 2 (p, z) + x(p, z) -1 = 0 , (45)
which reduces to eq. 25 when p = 1. The discriminant △(p) of the cubic eq. 45,

△(p) = -4(1 -p) 3 z 3 -(8 + 20p -p 2 )z 2 -4z , (46)
is strictly negative, so that eq. 45 has a single real root given explicitly by

x(p, z) = 1 -p 3 - 1 3z C(p) + △ 0 (p) C(p) , (47)
where

△ 0 (p) = (1 -p) 2 z 2 -3z , C(p) = △ 1 (p) + 3z -3△(p) 2 1/3 , (48)
and

△ 1 (p) = -z 2 9(2 + p) + 2(1 -p) 3 z . ( 49
)
The real root in eq. 47 is depicted by a thick black line in Fig. 3, panel (a). Turning next to the case p = 0.95, we observe a more complicated behaviour. First of all, we notice that similarly to the case p ≡ 1, there is an apparent transition at a certain z = z (bonds) c,1 (p) from the sequence convergence to a single limit x(p, z), to a subsequence convergence to alternating limits, when the even and odd terms in the recursion in eq. 44 converge, respectively, to some well-defined curves x even (p, z) and x odd (p, z). Curiously enough, however, these two curves x even (p, z) and x odd (p, z) meet each other again at some z = z (p) the sequence convergence to a single limit x(p, z) is again restored. While, evidently, the transition at z = z (bonds) c,1 (p) is the same as observed by Runnels; that being, the transition from a disordered into an ordered phase, the second transition taking place at z = z (bonds) c,2 (p) can be identified as a re-entrant transition from the ordered into a disordered phase. To the best of our knowledge, we are unaware of any previous report on the presence of the latter re-entrant transition for a repulsive lattice gas on a Bethe lattice, although it is quite plausible from the physical point of view. Indeed, for p < 1 the possibility of having two particles at neighbouring sites is not strictly prohibited and the penalty for having such a pair can be paid by an increase of the chemical potential. Below we study the loci and the nature of these transitions in more details. (p), eq. 55, and z > z

(bonds) c,2 (p)
, eq. 56, are the single limit solution of eq.

45.

Two dashed blue curves in the intermediate region z

(bonds) c,1 (p) < z < z (bonds) c,2
(p) describe alternating state and are defined by two distinct limit curves x odd (p, z) and x even (p, z), eqs. 52 and 53. The figure illustrates appearing of broken symmetry phase for p > 8/9. Phase transitions between disordered and ordered state occur at z As in the previous section, we reiterate the recursion in eq. 44 to get a recursion scheme involving the terms of the same parity only:

x N (p, z) = 1 + z(1 -p) 1 + z(1 -p)x 2 N -2 (p, z) 1 + zx 2 N -2 (p, z) 2 × 1 + z 1 + z(1 -p)x 2 N -2 (p, z) 1 + zx 2 N -2 (p, z) 2 -1 , x 0 (p, z) = 1 . ( 50
)
Supposing next that x N (p, z) → x(p, z) as N → ∞, we find that x(p, z) obeys the fifth-order equation which, again, conveniently factors into

zx 3 (p, z) -z(1 -p)x 2 (p, z) + x(p, z) -1 × × z 1 + z(1 -p) 2 x 2 (p, z) -pzx(p, z) + 1 + z(1 -p) = 0 . ( 51
)
While the first factor is just our previous eq. 45, which has one real root defining the solution in the single limit case, the second factor has two solutions:

x even (p, z) = pz + D(p, z) 2z (1 + z(1 -p) 2 ) , (52)
and

x odd (p, z) = pz -D(p, z) 2z (1 + z(1 -p) 2 ) , (53)
where the discriminant D(p, z) is given by :

D(p, z) = -z 4(1 -p) 3 z 2 + (8 -12p + 3p 2 )z + 4 = 4(1 -p) 3 z z -z (bonds) c,1 (p) z (bonds) c,2 (p) -z (54) with z (bonds) c,1 (p) = - 8 -12p + 3p 2 + p 3/2 √ 9p -8 8(1 -p) 3 , (55)
and

z (bonds) c,2 (p) = -8 + 12p -3p 2 + p 3/2 √ 9p -8 8(1 -p) 3 . ( 56
)
We note now the following : both z (p), so that a critical behaviour can emerge only for such values of p and for such a range of z. For p < 8/9 the second factor in eq. 51 does not have real solutions and the only real root of eq. 51 is determined by the single limit solution, eq. 47. This is precisely the behaviour we observe in Fig. 3 in which the panel (a) corresponds to the single-limit case with p = 0.85 < p shows the subsequence convergence to alternating limiting curves x even (p, z) and x odd (p, z) for z

(bonds) c,1 (p) < z < z (bonds) c,2
(p), and the sequence convergence for z < z (p) ∼ 1/8(1-p) 3 → ∞, which explains why the re-entrant transition is absent in the Runnels' case p ≡ 1. We note lastly that z (p) can be obtained from the stability analysis of the derivative of the function f (x N -1 (p, z)), which defines the recursion in eq. 44. Substituting into equation f ′ (x N -1 (p, z)) = -1 the solution 52 or 53 instead of x N -1 (p, z), and solving the resulting quadratic equation with respect to z, one gets the values defined in eq.55 and eq. 56.

Order parameter for p < 1

The order parameter δρ is defined as the difference of average occupations of adjacent layers j and j + 1 in the limit j → ∞. For p < 1, the average density is still described by eq. 33 with t j (z) replaced by t j (p, z), which obeys the recursion

t j (p, z) = b j (p, z) 1 -2 p (1 -b j-1 (p, z)) 1 -p b j-1 (p, z) t j-1 (p, z) , t 0 (p, z) = 1 , (57)
where b j (p, z) = (1 -x j (p, z))/p. The recursion in eq. 57 can be solved exactly to give where δ n,j is the Kronecker-delta, such that δ n,n = 1 and is zero, otherwise. Further on, inserting eq. 58 into eq. 33 and re-arranging the series, we get the expansion in eq. 36 with r j (z) replaced by r j (p, z), which now obeys the recursion relation

t j (p, z) = j-1 i=0 (-2) j j n=j-i b n (p, z) p (1 -b n (p, z)) 1 -p b n (p, z) 1-δ n,j , (58)
r j (p, z) = b j (p, z) 1 - p (1 -b j (p, z)) 1 -p b j (p, z) r j+1 (p, z) , ( 59
)
whose solution is given by

r j (p, z) = ∞ i=0 (-1) i i+j n=j b n (p, z) p (1 -b n (p, z)) 1 -p b n (p, z) 1-δ n,i+j . ( 60
)
Similarly to the expression in eq. 36 for the p ≡ 1 case, eq. 60 defines the average occupation at generation j in the general case when 0 ≤ p ≤ 1.

We focus now on the single limit case obtained either for p < p (p). In this case all b j (p, z) converge as j → ∞ to b(p, z) = (1 -x(p, z))p -1 with x(p, z) defined by eq. 47. In virtue of eq. 59, we have that the average occupation at generation j converges to

r(p, z) = b(p, z) (1 -p b(p, z)) 1 -p b 2 (p, z) , (61)
such that the order parameter δρ ≡ 0. Further on, we consider the critical case when p > p

(bonds) c and z ∈ [z (bonds) c,1 (p), z (bonds) c,2 (p)].
Here, b j (p, z) with even j converges in the limit j → ∞ to b -(p, z) = (1-x even (p, z))p -1 , eq. 52, while b j (p, z) with odd j converges to b + (p, z) = (1 -x odd (p, z))p -1 , eq. 53. Correspondingly, the average occupations on the generation j deep in the interior of an infinitely large tree are given by . Green curve corresponds to p = 0.9, the red one -to p = 0.95 and the blue curve -to p = 0.99.

r -(p, z) = b -(p, z) (1 -p b + (p, z)) 1 -p b + (p, z)b -(p, z) (62
for even j, and

r + (p, z) = b + (p, z) (1 -p b -(p, z)) 1 -p b + (p, z)b -(p, z) (63)
for odd j, respectively. In consequence, the staggered density obeys

δρ = r + (p, z) -r -(p, z) = b + (p, z) -b -(p, z) 1 -p b + (p, z)b -(p, z) = x even (p, z) -x odd (p, z) 1 -p + p (x even (p, z) + x odd (p, z) -x even (p, z)x odd (p, z)) , (64)
where x even (p, z) and x odd (p, z) are defined by eqs. 52 and 53. In Fig. 5 we depict the behaviour of the order parameter in eq. 64 for several values of p > p (bonds) c . Using eqs. 52 to 54, we can formally rewrite eq. 64 in the form

δρ = 4(1 -p) 3 z z -z (bonds) c,1 (p) z (bonds) c,2 (p) -z (1 -p) 3 z 2 + (1 -2(1 -p)p)z -p . ( 65
)
For z sufficiently close to z

(p))/(1 -τ )) ≃ z (bonds) c,1 (p) + z (bonds) c,1 (p) ln(z (bonds) c,1 (p))τ, ( 66
)
where τ = (T -T c )/T c is the deviation of the reduced temperature from the critical point T c . Then, substituting eq. 66 into eq. 67 and expanding the resulting expression in powers of τ , we find

δρ ≃ 4(1 -p) 3 z (bonds) c,1 (p) z (bonds) c,2 (p) -z (bonds) c,1 (p) ln(z (bonds) c,1 (p)) (1 -p) 3 z (bonds) c,1 (p) 2 + (1 -2(1 -p)p)z (bonds) c,1 (p) -p √ τ + O(τ 3/2 ), ( 67
)
which implies that the order parameter has the form δρ ∼ τ b with the mean-field critical exponent b = 1/2. In a similar way, we analyse the behaviour of the order parameter in the vicinity of the inverted transition point; that being, z = z (bonds) c,2

(p), to get a scaling behaviour with the same critical exponent b = 1/2.

Mean densities and the compressibility

We address next the question of the order of the transitions which our system undergoes at z = z . Consider first the behaviour of the mean particle density at these points. In the sequence convergence limit, i.e., for z < z (p), the average occupation of the generations with odd and even j deep in the interior of the infinite tree is the same, such that in virtue of eqs. 36 and 61 we have

ρ out (p, z) = r(p, z) = b(p, z) (1 -p b(p, z)) 1 -p b 2 (p, z) , ( 68
)
where the subscript "out" signifies that we deal with the behaviour out of the critical region, b(p, z) = (1 -x(p, z))/p and x(p, z) is given by eq. 47.

Further on, within the critical region z ∈ [z

(bonds) c,1 (p), z (bonds) c,2 (p)
], the average occupations of all the sites at even and odd generations are defined by eqs. 62 and 63, respectively, such that the mean density within the critical region obeys

ρ in (p, z) = 1 2 (r + (p, z) + r -(p, z)) = 1 2 b -(p, z) + b + (p, z) -2pb + (p, z)b -(p, z) 1 -pb + (p, z)b -(p, z) = 2 + (2 -3p)z 2 + 2z (2 -3p + (1 -p) 3 z) . (69)
Next, one may readily notice that since b(p, z (p)), such that the density is a piece-wise continuous function of z.

(bonds) c,1 (p)) = b -(p, z (bonds) c,1 (p)) = b + (p, z (bonds) c,1 (p)
Consider next the compressibility κ T defined as

κ T = 1 ρ 2 ∂ρ ∂µ = β ρ 2 ∂ρ ∂ ln(z) . ( 72
)
Away of the critical region, i.e., for z < z (p), we use the expression in eq. 68 to get . Green curve corresponds to p = 0.9, the red one -to p = 0.95 and the blue curve -to p = 0.99.

κ T /β = x(p, z) -1 + p x(p, z) 1 -x(p, z) p (1 -2x(p, t)) -(1 -x(p, z)) 2 p (x(p, t) -2) + 2 (1 -x(p, z)) 2 . ( 73
)
On the other hand, within the critical region, i.e., for z ∈ [z

(bonds) c,1 (p), z (bonds) c,2 (p)
], we take advantage of the expression in eq. 69 to get

κ T /β = 2z 3p -2 -4(1 -p) 3 z -(2 -3p)(1 -p) 3 z 2 2 + (2 -3p)z 2 . ( 74
)
In Fig. 6 we plot the compressibility κ T defined by eqs. 73 and 74 versus βµ for different values of p > p (bonds) c

. We observe that the compressibility exhibits a finite upward jump when the system enters into the ordered phase at z = z (bonds) c,1 (p), and also a finite downward jump when it re-enters the disordered phase for z exceeding z

(bonds) c,2 (p)
. This implies that both transitions are continuous, or second order in the Ehrenfest nomenclature. In particular, for p = 0.9, κ/β → 5/16 when z approaches z (p) from above, respectively. We close this Section with the following remark: Exact calculation of the average occupation at generation j in the general case 0 ≤ p ≤ 1, culminating at our eq. 60, which we performed here by generalising very directly the Runnel's approach, appears to be quite cumbersome even for the simple model of a purely repulsive lattice gas (Model I) and thus we can hardly expect that an analogous analysis can be carried out for the Model II with its multi-site Boolean interactions. In this regard, it might be instructive to present a somewhat simpler derivation of the main results of this Section using a different approach, which may also shed some light on the physical meaning of the expressions for ρ in (z) in eq. 69 and for ρ out (z) in eq. 68. Let ρ 0 denote the mean density at the central site of a finite Cayley tree. This density is formally defined as

ρ 0 = Z (bonds,1) N (p) Z (bonds,0) N (p) + Z (bonds,1) N (p) , (75) where Z (bonds,1) N (p) and Z
(bonds,0) N (p) are given by eqs. 40. Taking into account the definition of x N (p, z) in eq. 43, we arrive at the following expression:

ρ 0 = zx 3 N (p, z) 1 + zx 3 N (z, p) . ( 76
)
Outside of the critical region, in the limit N → ∞, mean densities on each site of the tree are equal to the same quantity ρ out , and consequently, ρ out = ρ 0 . Plugging x N (z, p) → x(z, p) into eq. 76 and combining it with eq. 45, we recover the expression in eq. 68. Further on, within the critical region the situation is a bit more delicate, because here we have two alternating limits for N → ∞: x even (p, z) and x odd (p, z). Let ρ even (p, z) and ρ odd (p, z) obey

ρ even (p, z) = zx 3 even 1 + zx 3 even , ρ odd (p, z) = zx 3 odd 1 + zx 3 odd . (77)
Defining next the mean density in the critical region as ρ in = (ρ even + ρ odd )/2, and using eqs. 52 and 53, we obtain for ρ in exactly the same expression as the one in eq. 69.

Model II : Catalytic sites

We turn next to the annealed version of the Model II and examine the critical behaviour of the annealed partition function defined by our eqs. 5 and 6, as the function of the reaction probability p and of the chemical potential µ (or the activity z). Similarly to the previously considered case of catalytic bonds, we decompose formally Z (sites) (p) in eq. 5 as

Z (sites) (p) = Z (sites,0) N (p) + Z (sites ,1) N (p), where Z (sites,0) N (p) and Z
(sites,1) N (p) denote the disorder-averaged partition functions of the Cayley tree with a vacant and an occupied central sites. One immediately notices that the partition function of the entire tree with a vacant central site factorises into the product of the partition functions defined on the subtrees; that being, we still have Z

(sites,0) N (p) = G 3 N (0; p)
, where G N (0, p) is the partition function in eq. 5 defined on a subtree (see, Fig. 1, right panel) with a vacant root. Note that, evidently, G N (0; p = 1) = B N (0; p = 1) = B N (0).

The case of a tree with an occupied central site is more difficult, compared to the Model I, since here different subtrees emanating from the occupied central node are effectively coupled via the Boolean function Ψ 0 , eq. 6, and the way how they may get decoupled depends now on the occupations of the sites belonging to the first generation. In consequence, we have to represent Z

(p) = G 3 N (1, 1; p) z 2 (1 -p) 2 + 3 G 2 N (1, 1; p)G N (1, 0; p) z 2 (1 -p) + 3G N (1, 1; p)G 2 N (1, 0; p) z 2 + G 3 N (1, 0; p) z 2 , (78)
Further on, we find straightforwardly that G(n 0 , n 1 ; p) obey the following recursions:

G N (0, 0; p) = G 2 N -1 (0; p) , G N (1, 0; p) = zG 2 N -1 (0; p) , G N (0, 1; p) = 1 (1 -p)z (1 -p)G 2 N -1 (1, 0; p) + + 2(1 -p)G N -1 (1, 0; p)G N -1 (1, 1; p) + G 2 N -1 (1, 1; p) , G N (1, 1; p) = (1 -p)G N -1 (1, 0; p) + G N -1 (1, 1, p) 2 . ( 79
) Recalling next that G N (1; p) = G N (1, 1; p) + G N (1, 0; p) and G N (0; p) = G N (0, 0; p) + G N ( 0, 1; p)
, we find from eq. 79 the recursions obeyed by G N (1; p) and G N (0; p):

G N (0; p) = G 2 N -1 (0; p) + pzG 4 N -2 (0; p) + + 1 (1 -p)z G N -1 (1; p) -pzG 2 N -2 (0; p) 2 (80)
and

G N (1; p) = zG 2 N -1 (0; p) + G N -1 (1; p) -pzG 2 N -2 (0; p) 2 . (81)
Now, it is convenient to introduce new auxiliary functions

ψ N (p, z) = G N (1; p) zG 2 N -1 (0; p) ; ξ N (p, z) = G 2 N -1 (0; p) G N (0; p) , (82)
which obey, in virtue of eqs. 80 and 81, the following coupled recursion relations:

ψ N (p, z) = 1 + zξ 2 N -1 (p, z) ψ N -1 (p, z) -p 2 , ψ 0 (p, z) = 1 , ξ N (p, z) = 1 + pzξ 2 N -1 (p, z) + zξ 2 N -1 (p, z) (ψ N -1 (p, z) -p) 2 1 -p -1 , (83)
where ξ 0 (p, z) = 1. Note, that the parameter x N (p, z) used to simplify the recursion in the Model I is given now by x N (p, z) = ψ N (p, z)ξ N (p, z). Note, as well, that in the limit p → 1 all ψ N → 1, such that the second equation in eqs. 83 becomes ξ N = (1+zξ 2 N -1 ) -1 , which is just our previous eq. 20 written in term of the function ξ N (p, z).

We turn back now to our eq. 78. Using the representations

G N (1, 1; p) = G N (1; p) -G N (1, 0; p) = zG 2 N -1 (0; p) ψ N (p, z) -1 = zξ N (p, z)G N (0; p) ψ N (p, z) -1 , (84)
as well as our eqs. 79, we may rewrite eq. 78 formally as Finally, interpreting the definition of the function ξ N (p, z) in eq. 82 as a recursion relation with an evident "initial" condition G 0 (0; p) = 1, we find that G N (0; p) for an arbitrary N can be expressed via ξ N (p, z) as

Z (sites,1) N (p) = zξ 3 N (p, z)G 3 N (0; p) (1 -p) 2 ψ N (p, z) -p 3 + p(1 -p) 2 , (85)
G N (0; p) = N j=1 ξ -2 N-j j (p, z) , (86)
such that eq. 12 yields for the pressure of the adsorbate

P II (T, µ) = 1 β lim M →∞ ln Z (sites) (p) M = 3 β lim M →∞ ln G N (0; p) M + 1 β lim M →∞ 1 M ln 1 + zξ 3 N (p, z) (1 -p) 2 (ψ N (p, z) -p) 3 + p(1 -p) 2 = 1 β ∞ j=1 2 -j ln 1 ξ j (p, z) , (87)
where the sum in the last line defines the desired thermodynamic limit result for the entire Cayley tree. In Appendix B we explain how one may subtract the contribution due to the boundary sites and get an analogous expression for the deep interior of the Cayley tree -the Bethe lattice.

Solution of the recursion relations for p < 1

Before we proceed further, it seems again expedient to get first some general understanding of the behaviour of the recursion scheme defined by eqs. 83. To this end, we generate first 20 terms of ψ N (p, z) and of ξ N (p, z) for three different values of the parameter p: p = 0.6, p = 0.8 and p = 0.9. These results are plotted in Figs. 7 and8. We notice that the situation appears to be somewhat similar to the one encountered for the annealed Model I in the sense that also here two clearly distinct behaviours are observed : for p = 0.6 we have a convergence to the single limit curve for all values of z, while for larger p there is a range of z with an apparent single limit convergence, and a bounded region in which odd and even terms seemingly converge to alternating curves. This signifies that, first, there exists some critical value p (sites) c,1 of the parameter p, which lies somewhere in-between p = 0.6 and p = 0.8. Second, it shows that also for the annealed Model II for p > p (p), at which the systems re-enters into a disordered phase. There is, however, a notable distinction between the behaviour of ψ N (p, z) and ξ N (p, z) for p = 0.8 and p = 0.9: while for z close to z , where p

(sites) c,2
is calculated from eq. B.10.

similar smooth behaviour, in the vicinity of z (sites) c,2

(p) the functions ψ N (p, z) and ξ N (p, z) approach the single limit curve much more abruptly for p = 0.9 than for p = 0.8. This hints that for the former case we may encounter a phase transition of a different type.

We turn first to the single limit case. Supposing that ξ N (p, z) → ξ(p, z) and ψ N (p, z) → ψ(p, z) as N → ∞, we have that eqs. 83 become

ψ(p, z) = 1 + zξ(p, z) 2 ψ(p, z) -p 2 , ξ(p, z) =    1 + pzξ 2 (p, z) + zξ(p, z) 2 ψ(p, z) -p 2 1 -p    -1 . ( 88
)
Further on, it follows from eq. 88 that

ψ(p, z) -p = (1 -p) 1 ξ(p, z) -pzξ 2 (p, z) , (89)
which permits us to write down a closed-form equation determining ξ :

f 1 (ξ(p, z), p, z) = p 2 (1 -p)ξ 7 (p, z)z 3 -2p(1 -p)ξ 4 (p, z)z 2 + ξ N (p, z) =    1 + z(1 -p) p(1 -p) + K 2 N -2 (p, z) p(1 -p)zξ 2 N -2 (p, z) + K N -2 (p, z) 2    -1 , (91)
where we denote

K N -2 (p, z) = 1 -p + zξ 2 N -2 (p, z) ψ N -2 (p, z) -p 2 . ( 92
)
Further on, we assume that ξ N (p, z) → ξ(p, z), ψ N (p, z) → ψ(p, z) and K N (p, z) → K(p, z) as N → ∞, so that we may rewrite eqs. 91 and 92 as:

ψ(p, z) = 1 + z(1 -p) 2 K 2 (p, z) (p(1 -p)zξ 2 (p, z) + K(p, z)) 2 , (93)
ξ(p, z) = 1 + z(1 -p) (p(1 -p) + K 2 (p, z)) (p(1 -p)zξ 2 (p, z) + K(p, z)) 2 -1 (94)
where

K(p, z) = 1 -p + zξ 2 (p, z) ψ(p, z) -p 2 . ( 95
)
As a matter of fact, from our eq. 95 we can express ψ(p, z) in eq. 93 via K(p, z), which will result in a closed system of equations for K(p, z) and ξ(p, z). Next, expressing K(p, z) from the second equation via ξ(p, z) and plugging the result into the first equation, we find a general closed-form equation for ξ(p, z). Some straightforward (but rather tedious) calculations, which we omit here, show that the closed-form equation for ξ(p, z) can be cast into the form

f 1 (ξ(p, z), p, z)f 2 (ξ(p, z), p, z) = 0 , (96)
where function f 1 (ξ(p, z), p, z) is defined in eq. 90. The function f 2 (ξ(p, z), p, z) is a polynomial of the fourteenth order in ξ(p, z) and a polynomial of the eighth order in z, and is presented in an explicit form in Appendix A (see eq. A.1). Now, the root of the f 1 (ξ(p, z), p, z) = 0, eq. 90, defines the solution in the single limit case. On the other hand, similarly to the situation described in Sec. 4, the roots of the equation

f 2 (ξ(p, z), p, z) = 0 , (97)
define the solutions in the alternating limits case. Numerical analysis of a strongly nonlinear eq. 97 together with some complementary arguments presented in Appendix C permit us to draw the following conclusions:

• Two alternating limit solutions appear only for p > p (sites) c,1 ⋍ 0.794, which value is specific, of course, to the coordination number three of the Bethe lattice.

• the ordered phase can only exist for p > p

(p).

• There is another critical value of the reaction probability p = p , while the full diamond shows the location of the tricritical point, at which the lines of continuous and discontinuous transitions meet.

which are either less than z (p), mean densities at each site (including the central one) are equal. Consequently, supposing that for the Bethe lattice ψ N (p, z) = ψ(p, z) and ξ N (p, z) = ξ(p, z), and using eq. 89, we obtain

ρ (sites) out (p, z) = z [(1 -p)(1 -pzξ 3 (p, z)) 3 + pξ 3 (p, z)] 1 + z [(1 -p)(1 -pzξ 3 (p, z)) 3 + pξ 3 (p, z)] . (100)
Compressibility in this case can be found from the definition in eq. 72. Substituting eq. 100 into eq. 72, we find that the compressibility away of the critical region obeys

κ T /β = 1 z (1 -p)(1 -pzξ 3 (p, z)) 3 + pξ(p, z) 3 + + 3pξ 2 (p, z) (1 -p) (1 -pzξ 3 (p, z)) (3zξ ′ (p, z) + ξ(p, z)) + ξ ′ (p, z) ((1 -p)(1 -pzξ 3 (p, z)) 3 + pξ(p, z)) 3 2 , ( 101
)
where the derivative ξ ′ (p, z) = ∂ξ(p, z)/∂z is defined as:

∂ξ(p, z) ∂z = - ∂f 1 (ξ(p, z), p, z) ∂z ∂f 1 (ξ(p, z), p, z) ∂ξ(p, z) -1 . (102)
We turn next to the behaviour in the critical region. We note that here, in contrast to the Model I, the exact calculation of the densities in the alternating layers is hardly possible due to a very complicated form of the recursion schemes obeyed by ξ N (p, z) and ψ N (p, z). Therefore, we resort here to an approach described at the end of the previous Section. Within the critical region, i.e., for p > p even (p, z))/2. Using eq. 99, we then define the mean density of sites at the odd generation, ρ

(p, z) = zξ 3 odd (p, z) (ψ odd (p, z) -p) 3 + p(1 -p) 2 (1 -p) 2 + zξ 3 odd (p, z) (ψ odd (p, z) -p) 3 + p(1 -p) 2 , (103)
and

ρ (sites) even (p, z) = zξ 3 even (p, z) (ψ even (p, z) -p) 3 + p(1 -p) 2 (1 -p) 2 + zξ 3 even (p, z) (ψ even (p, z) -p) 3 + p(1 -p) 2 .( 104
)
Using eqs. 103 and 104, we can define the order parameter δρ = |ρ odd -ρ even |. The latter is depicted in Fig. 11 as the function of βµ for the case p = 0.8 (i.e., for p within the interval p

(sites) c,1 < p < p (sites) c,2
) and for the case p = 0.9 (i.e., for p such that p > p (sites) c,2

). One observes a markedly different behaviour at z = z (p) and smoothly approaches zero, while in the latter case the order parameter vanishes discontinuously at z = z (sites) c,2 (p). The compressibility in this case can be found if we replace ρ in eq. 72 by ρ in . Then, the only properties we need to know are the derivatives of ξ odd (p, z), ξ even (p, z), ψ odd (p, z) and ψ even (p, z) with respect to z. The latter can be calculated rather straightforwardly and we obtain the following result:

κ T /β = 2z A 1 A ′ 2 -A 2 A ′ 1 A 1 + 2A 2 2 , (105)
where the explicit expressions for the functions A A ′ 2 (s(p, z), σ(p, z)) are also presented in an explicit form in Appendix C (see eqs. C.17 and C.18).

Therefore, we obtain the compressibility for the whole range of βµ and arbitrary p. The results for p = 0.8 and p = 0.9 are presented in Figs. 12 and 13, respectively. One notices that for p = 0.8 the density is piece-wise continuous with cusps at z = z (p). This implies that both transitions are continuous transitions. On the other hand, for p = 0.9, (which value of p exceeds p (sites) c, 2), we observe that the density is piece-wise continuous at z = z (p), which signifies that in this range of values of the reaction probability p the re-entrant transition is discontinuous. Interestingly enough, behaviour of the compressibility after the re-entrant transition is markedly different depending whether the latter is continuous or discontinuous. In case of a continuous transition the compressibility jumps downwards at the transition point, passes through a maximum and then decreases monotonically with z. In case of a discontinuous transition the compressibility makes an upward jump and then monotonically decreases with z.

Conclusions

To conclude, we have studied here equilibrium properties of two lattice-gas models of catalytically-activated A + A → ⊘ reactions on a lattice of adsorption sites. In both models, the A particles are assumed to be in thermal contact with their vapour phase (a reservoir maintained at a chemical potential µ), adsorb onto empty adsorption sites and desorb from the lattice. We considered two different ways of modelling such reactions: in the Model I we assumed that some fraction p of the bonds connecting neighbouring adsorption sites possesses special catalytic properties so that any two As appearing on the sites connected by such a bond react instantaneously and desorb. In the Model II, we stipulated that some fraction p of the adsorption sites possesses such properties and the reaction takes place once at least one of the neighbouring A particles occupies a catalytic adsorption site.

We focused on the case of annealed disorder in the distribution of the catalytic bonds or sites, which is tantamount to the situation when the reaction between two As may take place at any point on the lattice but happens with a finite probability p, which means that the reaction is not perfect (instantaneous) but is characterised by some finite reaction constant. In this case the Model I describes a lattice gas with soft, purely repulsive nearest-neighbour interactions, while the Model II represents a lattice gas with particular multi-site interactions of particles: here, the interaction of an adsorbed particle with its nearest environment has a Boolean form -it is either zero, in case when the particle does not have any neighbour, or is a constant independent of the actual number of neighbouring adsorbed particles, if at least one of them is present. We provided exact analytical solutions for the annealed versions of both Models I and II on the interior of the Cayley tree -the Bethe lattice, and showed that they exhibit a rich "critical" behaviour with respect to µ and p, characterised by a transition into an ordered state, which for both models is continuous with a finite jump in compressibility, and a re-entrant transition into a disordered phase, which is continuous with a finite jump in compressibility for the Model I and, depending on the value of p, may be either continuous or discontinuous with a finite jump of density, for the Model II.

The wealth of critical phenomena which we observed for the Models I and II for a Bethe lattice geometry, certainly merits further, more deep investigation including a study of the kinetic behaviour in both models, as well as an analysis of both models on more realistic, regular or random adsorbent lattices with quenched (above and below the percolation threshold) and/or annealed distributions of the catalysts. Lastly, we would like to remark that it is interesting to extend our analysis over the case of a monomer-monomer model A + B → ⊘, which involves two types of particles [5,[8][9][10]. This model with a finite reactivity p has been already studied for one-dimensional finite systems in Ref. [43] and also for two-dimensional lattices with an annealed distribution of the catalytic bonds [44,45]. In Refs. [44,45] it was shown that, in particular, the grand canonical partition function for the A + B → ⊘ model can be mapped onto the partition function of the general spin S = 1 model [46,[START_REF] Lawrie | Phase Transitions and Critical Phenomena[END_REF], which permitted to exploit the large number of results available for the latter (see, e.g., Ref. [45]). For instance, for the symmetric case of equal chemical potentials for both species, a phase transition was predicted from the phase in which the particles of both sorts have the same densities, to the phase in which one sort of particles prevails. This transition can be of the first order or a continuous one, depending on the precise values of the system's parameters. It was also shown in Ref. [44] that in some parameter space the monomermonomer model defined on a honeycomb lattice with an annealed disorder in placement of the catalytic bonds, reduces to the original Blume-Emery-Griffiths model [START_REF] Blume | [END_REF], whose solution can be obtained in a closed form via a mapping to a zero-field Ising model on a regular honeycomb lattice, which can be solved exactly and exhibits a symmetrybreaking continuous transition [49,50]. For the A + B → ⊘ reaction on a lattice with catalytic sites no analytical results are available at present. and A 2 = z 2 (1-p) 2 (1 -pzσs + p 2 z 2 σ 3 ) 3 + p(1-p) s 3 -3σs -3pzσ(s 4 -4σs 2 + 2σ 2 )+ + 3(pzσ) 2 (s 5 -5σs 3 + 5σ 2 s)-(pzσ) 3 (s 6 -6σs 4 + 9σ 2 s 2 + 2σ 3 ) -p 2 σ 3 . (C. 14) In turn, the derivatives of A 1 and A 2 with respect to z obey zA ′ 1 = A 1 + z 2 3(1 -p) -p(1 -pzσs) 2 (σs + z (sσ ′ + σs ′ )) -2p 2 zσ 2 (2σ + 3zσ ′ ) + + p 3 z 2 σ 3 (3σs + z (4sσ ′ + σs ′ )) + 3p(s 2 s ′ -σs ′ -sσ ′ ) , (C.15) and zA ′ 2 = 2A 2 + z 3 3(1 -p) 2 (1 -pzσs + p 2 z 2 σ 3 ) 2 -p σs + z (sσ ′ + σs ′ ) + + 3p 2 z 2σ(σ+zσ ′ )(s 5 -5σs 3 +5σ 2 s)+5zσ 2 (s 4 s ′ -s 3 σ ′ -3s 2 σs ′ +2sσσ ′ +σ 2 s ′ ) --3p (σ + zσ ′ )(s 4 -4σs 2 + 2σ 2 ) + 4zσ(s 3 s ′ -s 2 σ ′ -2sσs ′ + σσ ′ ) + + p 2 zσ 2 (2σ + 3zσ ′ ) + 3p 2 σ 2 σ ′ + p(1 -p) 3s 2 s ′ -3 (sσ ′ + σs ′ ) + + 6zσ 3 (s 5 s ′ -s 4 σ ′ -4s 3 σs ′ + 3sσ(sσ ′ + σs ′ ) + σ 2 σ ′ ) --p 3 z 2 3σ 2 (σ + zσ ′ )(s 6 -6s 4 σ + 9σ 2 s 2 + 2σ 3 ) (C.16)

Lastly, we present explicit expressions for the derivatives s ′ and σ ′ , which enter the expressions for the compressibility :

s ′ = (σ+zσ ′ ) 1 pz 2 σ 2 - (1-p) 1-(1-p)pz 2 σ 2 1- 2(1-p)pz 2 σ 2 1-(1-p)pz 2 σ 2 (C.

  ) the lattice will partition into two sub-lattices -the one densely packed with A particles and the second -almost empty. c) Increasing µ further, a re-entrant (inverted) transition into a disordered phase at µ = µ should take place. Both transitions are continuous, with a finite jump in compressibility. When p → 1, the re-entrant transition disappears since µ (bonds) c,2

2 ≈

 2 0.813 and no critical behaviour emerges for p < p

  characterised by a finite jump of the compressibility. c) Our further analysis reveals that the behaviour for larger values of µ depends on whether p belongs to the interval p

Figure 1 .

 1 Figure 1. Left: Cayley tree with 3 generations emanating from the central site O. Right: A subtree with 3 generations emanating from the central site O.

Figure 2 .

 2 Figure 2. x N = x N (z) in eq. 20 for N = 1, 2, . . . , 100 plotted versus βµ = ln(z). Red curves are x N (z)-s with even N , while the green ones correspond to the terms with odd values of N . Arrows indicate the direction of growth of N . Thick black curve depicts the behaviour of the real root of eq. 25 (see eq. 26) on the interval [0, βµ = 2 ln(2)]. Two dashed blue curves starting at βµ = 2 ln(2) are the roots of the second term in eq. 29 defining two distinct limit curves x odd (z) and x even (z). The figure illustrates continuous phase transition from disordered state described by solution of eq. 25 (black line) to symmetry broken phase with alternating layers of different densities described by solution of eq. 29 (dashed blue lines). The phase transition occurs at z c = 4 at which all three curves meet.

  ) denotes the disorderaverage partition functions of the Cayley tree with a vacant (occupied) central site. In this case eqs. 15 and 16 read Z (bonds,0) N

Figure 3 .

 3 Figure 3. x N = x N (p, z) in eq. 44 for N = 1, 2, . . . , 20 plotted versus βµ = ln(z). Red curves are x N with even N , while the green ones correspond to odd values of N . Panel (a): p = 0.85. Thick black curve describe disordered state and is the real root of eq. 45 given explicitly by eq. 47. Panel (b): p = 0.95. Thick black curves in the regions z < z (bonds) c,1

  , where blue lines meet with black lines.

  in eqs. 52 and 53 are real and positive only for p > p (bonds) c = 8/9. The discriminant D(p, z) > 0 and hence, the solutions x even (p, z) and x odd (p, z), eqs. 52 and 53, are real and positive when p > p

  . In Fig.4we depict the parametric curve which solves the equation D(p, z) = 0. The region encircled by this curve corresponds to the critical region z

Figure 4 .

 4 Figure 4. Phase diagram for the annealed Model I. Blue line separating disordered and ordered (alternating) phases is the solution of equation D(p, z) = 0 with D(p, z) defined by eq. 54. The full circle indicates the location of the critical value of the parameter p = p (bonds) c = 8/9.

(bonds) c, 1 (

 1 p) or for z > z (bonds) c,2

Figure 6 .

 6 Figure 6. Compressibility κ T /β in eqs. 73 and 74 versus βµ for different p > p (bonds) c

  from below, and κ/β → 5/12 when z → z (bonds) c,1 (p) from above. For the reentrant transition for the same value of p, we have κ/β → 1/5 for z → z (bonds) c,2 (p) from below, and κ/β → 4/15 for z → z (bonds) c,2

(

  sites,1) N (p) as the sum of all configurations with different occupations of the sites of the first generation; namely, all configurations with an occupied central site and a) all three sites in the first generation occupied, b) two sites occupied and one -vacant, c) one site occupied and two -vacant, and c) all three sites in the first generation vacant. Introducing next auxiliary functions G(n 0 , n 1 ; p), Order-disorder transitions in reactive lattice gases M. Dudka et al which denote the partition function of a subtree with the root having an occupation number n 0 and the site at the first generation having the occupation number n 1 , we represent Z (sites,1) N (p) by directly counting all possible configurations as Z (sites,1) N

  which together with the relation Z (sites,0) N (p) = G 3 N (0; p) define the partition function Z (sites) (p) of the entire Cayley tree.

  (sites) c,1 there are some critical value z (sites) c,1 (p), at which the systems enters from a disordered phase into an ordered one, and some critical value z (sites) c,2 (p) > z (sites) c,1

Figure 7 .

 7 Figure 7. ψ N = ψ N (p, z) in eqs. 83 versus βµ for N = 1, 2, 3, . . . , 20. Green curves correspond to odd N , while the red ones -to even values of N . Panel (a) : p = 0.6, Panel (b): p = 0.8 and Panel (c): p = 0.9. Thick black lines are defined by eq. 89 with ψ(p, t) being the single limit solution of eq. 90. Thick blue dashed lines defining the solution in the alternating limits are the roots of eq. 97. Vertical thin dotted line in Panel (c) shows the discontinuity in ψ N (p, z) at z = z (sites) c,2 (p) emerging for

  (sites) c,1 and such z which obey the doublesided inequality z (sites) c,1 (p) < z < z (sites) c,2

Figure 10 .

 10 Figure 10. Phase diagram for the annealed Model II. Solid line demarcates a continuous phase transition between ordered (alternating) and disordered phases, while the dashed curve shows the discontinuous re-entrant transition between these two phases. The full circle corresponds to p (sites) c,1

∆ΡFigure 11 .

 11 Figure 11. Order parameter δρ versus βµ for different p > p c,1 . Blue curve corresponds to p = 0.8, the red one -to p = 0.9.

  z), and at the even generation, ρ (sites) even (p, z), as ρ (sites) odd

  -in the former case δρ attains a maximal value for z < z (sites) c,2

1 =Figure 12 .

 112 Figure 12. Mean density ρ, panel (a), and the compressibility κ T /β, panel (b), versus βµ for p = 0.8. The red curves in panels (a) and (b) correspond to eqs. 100 and 101, respectively. The blue curve in the middle of the panel (a) depicts ρ in = (ρ (sites) odd (p, z) + ρ (sites) even (p, z))/2 with ρ (sites) odd (p, z) and ρ (sites) odd (p, z) defined by eqs. 104 and 103. The blue curve in the panel (b) shows the compressibility in the critical region defined by eq. 105.

Figure 13 .

 13 Figure 13. Mean density ρ, panel (a), and the compressibility κ T /β, panel (b), versus βµ for p = 0.9. The red curves in panels (a) and (b) correspond to eqs. 100 and 101, respectively. The blue curve in the middle of the panel (a) depicts ρ in = (ρ (sites) odd (p, z) + ρ (sites) even (p, z))/2 with ρ (sites) odd (p, z) and ρ (sites) odd (p, z) defined by eqs. 104 and 103. The blue curve in the panel (b) defines the compressibility in the critical region defined by eq. 105.

  and at z = z (sites) c,2 (p)), while the compressibility exhibits finite jumps when the system enters into the ordered phase at z = z (sites) c,1 (p) and re-enters into the disordered phase at z = z (sites) c,2

  , and varies discontinuously at z = z (sites) c,2

1 .

 1 ) andσ ′ = -∂f σ (σ, p, z) ∂z ∂f σ (σ, p, z) ∂σ -(C.18)

Acknowledgments

The authors wish to thank Professor S. Dietrich for fruitful discussions. MD acknowledges the hospitality of the Max-Planck-Institute Stuttgart, where the most of work on this project has been done. OB is partially supported by the European Research Council Grant No. FPTOpt-277998. MD is supported in part by the FP7 EU IRSES project No. 612707 "Dynamics of and in Complex Systems".

is obtained from eq. B. 10.

which is a seventh-order equation in ξ(p, z). Note that for p ≡ 1, ψ(p, z) ≡ 1 such that eq. 90 reduces to a depressed cubic equation of the form zξ(p, z) 3 + ξ(p, z) -1 = 0, which is just our previous eq. 25. For general p < 1, we can only solve equation eq. 90 numerically. The real roots ξ(p, z), eq. 90, for fixed p = 0.6, 0.8, 0.9 are depicted by a thick black line in Fig. 8, (and correspondingly, the roots of ψ N (p, z) are depicted in Fig. 7), demonstrating a convergence of the recursion to this single limit curve for any z for p = 0.6, as well as for z < z (p) stay undefined, for the moment). Consider now the alternating limits case. As in the case of the Model I, (and also following the Runnels' analysis [31]), we reiterate eqs. 83 once more to get a recursion involving only the terms of the same parity with respect to N. This gives ⋍ 0.813. Green curves correspond to odd N , while the red ones -to even values of N . Thick black lines are the single limit solutions of eq. 90 while the thick dashed lines are the alternating limits solutions of eqs. 97. Note that ξ N (p, z) becomes multi-valued at z = z (sites) c,2 (p).

• For p > p (p) and the alternating limit solution of eqs. 97 become multi-valued and no longer represents the actual limits approached by ψ N (p, z) and ξ N (p, z) as N → ∞ (see Figs. (p) → ∞ when p → 1, such that the re-entrant transition disappears, as it should.

Mean density, compressibility and the order parameter for p < 1

We start with the calculations of the mean density and of the compressibility in the single limit case, i.e., outside of the critical region. The mean density at the central site of the Cayley tree is formally defined as

where Z

(sites,1) N (p) is defined by eq. 85 and Z (sites,0) N (p) = G 3 N (0; p). This implies that the mean density at the central site of a Cayley tree with N generations is given by

For a deep interior of an infinite Cayley tree, i.e., for the Bethe lattice, and outside of the critical region, that is, for any z and p < p

, or for an arbitrary p but z

We present here the polynomial f 2 (ξ(p, z), p, z) which enters eq. 96:

f 2 (ξ(p, z), p, z) = 1 + ξz ξ 13 q 3 p 4 z 6 -ξ 12 q 3 p 4 z 6 -2ξ 10 q 4 p 3 z 6 + 2ξ 11 q 3 p 3 z 5 (qz+1) + ξ 7 q 2 p 2 z 3 (3 -qz(5 -8qz))

where q = 1 -p.

Appendix B. Pressure for the Bethe lattice

To derive an explicit expression for the substrate pressure of a deep interior of the Cayley tree -the Bethe lattice, we use the procedure which is well described in Ref. [42].

Appendix B.1. Pressure for the annealed Model I

The procedure described in Ref. [42] consists of the following steps. First, using eq. 11 we can write the pressure of an adsorbate on an N-generation Cayley tree in the form : Rewriting the latter expression in terms of the pressure P I,N -1 for an N -1-generation Cayley tree we therefore arrive at

Reiterating this procedure n times, we find the following recursion relation for P I,N :

where P I,N n is now the pressure of an n-generation deep interior of the Cayley tree -an n-generation Bethe lattice, which reads

Lastly, we determine P I by dividing the above expression by βM, and taking next the limit N → ∞, (in which limit all x N -k (p, z) ≡ x(p, z)), which eventually yieldes the desired expression

We note next that the mean density can be obtained from eq. B.7 in the standard way by differentiating the pressure with respect to z, i.e., vρ(z) = βz∂P I /∂z, and excluding ∂x(p, z)/∂z and z from the resulting expression using eq. 45. In doing so, one arrives at the expression given in eq. 68. Pressure in the critical region, in which two alternating limits for x N (p, z) exist, is to be written in the form

The mean density in this case is derived using eqs. 52 and 53, to get eventually the expression in eq. 69.

Appendix B.2. Pressure for the annealed Model II Applying the analogous procedure for the derivation of the pressure for the Model II, we get following expression

Requiring the continuity of the pressure at the phase transition point, we determine the critical value of the activity z f (p) for the first order phase transition:

Together with eq. 88 and eqs. C.1 to C.4, the expression in eq. B.10 implicitly determines z f (p). The solution is depicted by a dashed line in Fig. 10.

Appendix C. Alternating limits solution for the annealed Model II

Here we present some complementary analysis of the alternating limits solution for the annealed Model II, which also turns out to be useful for the calculations of the staggered density and of the compressibility. We first revisit our eq. 91 and suppose that ξ N (p, z) and ψ N (p, z) with odd and even N converge as N → ∞ to ξ odd = ξ odd (p, z) and ψ odd = ψ odd (p, z), and to xi even = ξ even (p, z) and ψ even = ψ even (p, z), respectively. Then, in accord with eq. 83, we find that the latter obey:

Observe next that

such that ξ odd and ξ even obey, respectively,

Define next two auxiliary symmetric functions s = s(p, z) = ξ odd + ξ even and σ = σ(p, z) = ξ odd -ξ even . Multiplying both sides of eq. C.7 by ξ even and of eq. C.8 -by ξ odd , we arrive, after some slight rearrangements of both equations, to the following equations which define s and σ:

and

From eq. C.9 we readily find that Further on, ξ odd and ξ even are obtained from

Real values of ξ odd and ξ even obtained in this way reproduce the alternating limits solution, shown by the blue dashed curves in Figs. 7 and8. Further on, the density ρ in within the critical region can also be expressed through the symmetric auxiliary functions s and σ. Substituting next the resulting expression into eq. 72, we arrive at our eq. 105, which defines the compressibility within the critical region and is presented in the main text. In this equation the functions A 1 = A 1 (s, σ) and A 2 = A 2 (s, σ) are explicitly given by A 1 = z (1 -p) (1 -pzσs) 3 -6p 2 z 2 σ 3 + 3p 3 z 3 σ 4 s + p(s 3 -3σs) , (C. 13)