Order-disorder transition in a two-dimensional associating lattice gas

TL;DR

This study investigates the order-disorder phase transition in a two-dimensional associating lattice gas model using Monte Carlo simulations and Husimi lattice solutions, revealing a continuous transition in simulations and a discontinuous one in mean-field approximations.

Contribution

It provides a detailed comparison of phase transition behaviors in an associating lattice gas model using different methods, highlighting the universality class and the nature of the transition.

Findings

Transition is continuous in simulations with 3-state Potts universality.

Husimi lattice predicts a discontinuous transition with higher transition temperatures.

Three scenarios proposed for the transition behavior at infinite chemical potential.

Abstract

We study an associating lattice gas (ALG) using Monte Carlo simulation and solutions on Husimi lattices. In this model, the molecules have an orientational degree of freedom and the interactions depend on the relative orientations of nearest-neighbor molecules. We focus on the transition between the high-density liquid (HDL) phase and the isotropic gas phase in the limit of full occupancy ($\mu \to \infty$). Simulation results show a continuous phase transition at $\tau_c=k_BT_c/\gamma=0.4763(1)$ (where $-\gamma$ is the bond energy) between the low-temperature HDL phase, with a non-vanishing mean orientation of the molecules, and the high-temperature isotropic phase. Results for critical exponents and the Binder cumulant indicate that the transition belongs to the three-state Potts model universality class, even though the ALG Hamiltonian does not have the full permutation symmetry of the Potts model. In contrast with simulation, the Husimi lattice results furnish a discontinuous phase transition, characterized by a discontinuity of the nematic order parameter. The transition temperatures ($\tau_c=0.51403$ and $0.51207$ for trees built with triangles and hexagons, respectively) are slightly higher than the one found via simulation. Since the Husimi lattice studies show that the ALG phase diagram features a discontinuous gas-HDL line for finite $\mu$, three possible scenarios arise. The first is that in the limit $\mu \to \infty$ the first-order line ends in a critical point; the second is a change in the nature of the transition at some finite chemical potential; the third is that the entire line is one of continuous phase transitions. Results from other ALG models and the fact that mean-field approximations show a discontinuous phase transition for the three-state Potts model (known to possess a continuous transition) lends some weight to the third alternative.

Tables

Table 1. Table 1: Estimates for the critical temperature. The first six columns (second line) list the estimates for T c subscript 𝑇 𝑐 T_{c} obtained from analysis of the quantities listed on the first line, while the seventh column contains the global estimate and its uncertainty.

$c$	${\chi }_{\theta }$	${\chi }_Q$	$U_{4,\theta }$	$U_{4,Q}$	$U_{4,Q}^{\prime }$	${\tau }_c$
0.47632(1)	0.47636(1)	0.4761(1)	0.47621(2)	0.47637(8)	0.4764(5)	0.4763(1)

Equations

$${\hat{e}}_1=\mathbf{i},$$

$${\hat{e}}_2=+\frac{1}{2}\mathbf{i}+\frac{\sqrt{3}}{2}\mathbf{j},$$

$${\hat{e}}_3=-\frac{1}{2}\mathbf{i}+\frac{\sqrt{3}}{2}\mathbf{j}.$$

$$u_{ij}=-\gamma {\left(\frac{4}{3}\right)}^2\left[1-{\left({\mathbf{\eta }}_i\cdot \mathbf{r}\right)}^2\right]\left[1-{\left({\mathbf{\eta }}_j\cdot \mathbf{r}\right)}^2\right],$$

$$\mathcal{H}\left(\mathbf{\eta },\mathbf{r}\right)=\sum _{⟨i,j⟩}{\sigma }_i{\sigma }_j[\epsilon +u_{ij}],$$

$$\theta =\frac{3}{2}\left(\frac{\max \left(n_1,n_2,n_3\right)}{N}-\frac{1}{3}\right)$$

$$Q=\frac{1}{N}\sqrt{n_1^2+n_2^2+n_3^2-n_1{n}_2-n_1{n}_3-n_2{n}_3}.$$

$$p\left[\Gamma \to {\Gamma }^{\prime }\right]=\min \left[1,\frac{\Omega \left[E\left({\Gamma }^{\prime }\right)\right]}{\Omega \left[E\left(\Gamma \right)\right]}\right].$$

$$⟨\mathcal{O}⟩\equiv \frac{1}{\mathcal{Z}}\sum _E\overline{\mathcal{O}\left(E\right)}\Omega \left(E\right)\exp \left(-\beta E\right),$$

$${\chi }_{\theta }\left(\tau \right)=\frac{L^2}{\tau }\left[⟨{\theta }^2⟩-{⟨\theta ⟩}^2\right],$$

$$c\left(\tau \right)=\frac{{\beta }^2}{L^2}\left[⟨E^2⟩-{⟨E⟩}^2\right].$$

$$\theta \approx L^{-{\beta }_{\theta }/\nu }\mathcal{F}\left(t{L}^{1/\nu }\right),$$

$$\chi \approx L^{\gamma /\nu }\mathcal{X}\left(t{L}^{1/\nu }\right),$$

$$c\approx c_0+L^{\alpha /\nu }\mathcal{C}\left(t{L}^{1/\nu }\right),$$

$$U_{4,\theta }(T)=1-\frac{⟨{\theta }^4⟩}{3{⟨{\theta }^2⟩}^2}.$$

$$\nu =\frac{2}{{\alpha }^{\prime }+2{\beta }^{\prime }+{\gamma }^{\prime }},$$

$$g_{1,j_1}^{(M+1)}$$

$$g_{2,j_2}^{(M+1)}$$

$$g_{3,j_3}^{(M+1)}$$

$$n_b($$

$$R_{i,j}=\frac{g_{i,j}}{g_{i,1}+g_{i,2}+g_{i,3}}.$$

$${\mathbf{R}}_{iso}^{\ast }=\left(\begin{array}{ccc}x& (1-x)/2& (1-x)/2\\ (1-x)/2& x& (1-x)/2\\ (1-x)/2& (1-x)/2& x\end{array}\right),$$

$${\mathbf{R}}_{nem,1}^{\ast }$$

$${\mathbf{R}}_{nem,2}^{\ast }$$

$${\mathbf{R}}_{nem,3}^{\ast }$$

$$J_{i,j}={\left(\frac{\partial R_i^{\prime }}{\partial R_j}\right)}_{R^{\ast }},$$

$${\varphi }_b=-\frac{1}{3}{k}_{B}T\ln \frac{Y_{M+1}}{Y_M^4},$$

$$Y_M=\sum _{j_1,j_2,j_3=1}^3{W}_{\{j_1,j_2,j_3\}}{g}_{1,j_2}^{(M)}{g}_{1,j_3}^{(M)}{g}_{2,j_1}^{(M)}{g}_{2,j_3}^{(M)}{g}_{3,j_1}^{(M)}{g}_{3,j_2}^{(M)}.$$

$$⟨n_i{⟩}^{(M)}=\frac{1}{3Y^{(M)}}\sum _{j_1,j_2,j_3=1}^3$$

$$⟨n_i⟩=\frac{{\sum }_{j_1,j_2,j_3=1}^3{n}_{i}f(j_1,j_2,j_3,\{R^{\ast }\})}{{\sum }_{j_1,j_2,j_3=1}^{3}f(j_1,j_2,j_3,\{R^{\ast }\})},$$

$$f(j_1,j_2,j_3,\{R^{\ast }\})$$

Full text

Phase transitions, critical phenomena, Monte Carlo simulations, lattice gas, Mean-field approximations

I Introduction

It is no news that water exhibits quite unusual thermodynamic behavior, characterized by a set of anomalies, among which the anomaly in density is the best known Chaplin (2018). In recent decades, many models were developed with the purpose of investigating the fundamental mechanisms which lead to the water anomalies. Among them, lattice models have attracted much attention due to their easy implementation and low computational cost. These models usually include soft-core potentials to take care of excluded volume effects, and orientational interactions to represent hydrogen bonding between molecules. So far, these models are only able to reproduce some of the anomalies qualitatively. In spite of this, they exhibit rich phase diagrams that provide an ideal environment for the study of phase transitions as well as the validation of new computational techniques.

The first orientational lattice model for water was proposed by Bell and Lavis (BL) bell and Lavis (1970); Bell and Lavis (1970); Lavis (1973). It is defined on a triangular lattice, in which each site can be either vacant or occupied by a molecule. The molecules possess three bonding directions with 120∘ between them, resulting in two orientational states per molecule. The model exhibits three phases: gas, low density liquid (LDL) and high-density liquid (HDL) Fiore et al. (2009) at low, intermediate and high chemical potentials, respectively. While the gas-LDL transition is known to be discontinuous, characterized by a jump in density, the nature of the liquid-liquid transition is still controversial. For instance, mean-field approximations (Bethe lattice solutions Barbosa and Henriques (2008) and cluster-variation methods Bruscolini, Pelizzola, and Casetti (2002)) point to a discontinuous phase transition, whereas Monte Carlo simulations Fiore et al. (2009); Šimėnas, Ibenskas, and Tornau (2014) show a continuous transition. There is no consensus regarding the universality class of this transition. Fiore et al. Fiore et al. (2009) assert that it falls the Ising universality class, while Šimėnas et al. Šimėnas, Ibenskas, and Tornau (2014) report three-state Potts critical exponents.

Following the ideas of Bell and Lavis, Henriques and Barbosa introduced a two-dimensional (2D) associating lattice gas (ALG) Henriques and Barbosa (2005). In this model, also defined on a triangular lattice, each molecule has four bonding arms and two inert ones, the latter taking opposite directions on the lattice. Two of the bonding arms are proton donors in hydrogen bonds, while two are receptors, leading to eighteen states per molecule. The model Henriques and Barbosa (2005) also exhibits gas, LDL and HDL phases, but in contrast to simulation results for the BL model, the LDL and HDL phases are separated by a discontinuous transition line that ends at a bicritical point. There is also a gas-HDL line of continuous transitions, which starts at the bicritical point and seems to extend to large values of the chemical potential. The LDL and gas phases are separated by a continuous and a discontinuous transition line, which connect at a tricritical point. The former line meets the LDL-HDL and gas-HDL lines at the bicritical point Szortyka et al. (2009).

Variants of the ALG model have been investigated, among them, three-dimensional (3D) Girardi et al. (2007); Szortyka et al. (2010); Furlan, Almarza, and Barbosa (2016) and symmetric versions Balladares, Henriques, and Barbosa (2007); Furlan, Fiore, and Barbosa (2015). The 3D model Girardi et al. (2007); Szortyka et al. (2010) is defined on a body-centered cubic lattice in which each molecule possesses four bonding, and four inert arms. Despite the differences in geometry and number of orientational states in relation to the original ALG model Henriques and Barbosa (2005), the three-phase (gas, LDL, HDL) behavior is preserved. While a phase diagram featuring two tricritical points was suggested by Buzanno et al. Buzzano, De Stefanis, and Pretti (2008) based on a cluster-variation approach, simulation results by Szortyka et al. Szortyka et al. (2010) indicate that there is in fact a tricritical and a bicritical point, similar to the 2D case Szortyka et al. (2009). In the 3D model, however, a gas-LDL coexistence line meets the continuous LDL-HDL and gas-HDL transition lines at the bicritical point Szortyka et al. (2010).

The symmetric ALG model makes no distinction between donor and receptor bonding arms. This leads to a simplification, since the number of states per particle is substantially reduced. Balladares et al. Balladares, Henriques, and Barbosa (2007) investigated this model on the triangular lattice and found only two discontinuous (gas-LDL and a LDL-HDL) transition lines in the phase diagram, each ending at a critical point. As an aside, let us remark that this very same scenario was reported in the former studies of the original 2D ALG model Henriques and Barbosa (2005) and of the 3D version Girardi et al. (2007). It turns out that they were incorrect, as demonstrated in more recent analyses Szortyka et al. (2009); Buzzano, De Stefanis, and Pretti (2008); Szortyka et al. (2010), as discussed above. In fact, the semi-analytical solution of the symmetric model Balladares, Henriques, and Barbosa (2007) on a Husimi lattice build with hexagons (which is a mean-field approximation for the triangular lattice) unveils a phase diagram with three coexistence lines (gas-LDL, gas-HDL and LDL-HDL) meeting at a triple point Oliveira, Stilck, and Barbosa (2010). Moreover, more recent simulations of this model have provided evidence that the critical points reported in Balladares, Henriques, and Barbosa (2007) are actually tricritical points Furlan, Fiore, and Barbosa (2015), so that the thermodynamic behavior of this model is closer to the other versions. In contrast with these cases and with the mean-field results, however, a gas-HDL transition has not yet been observed in simulations of the 2D symmetric ALG, so that the existence of a continuous order-disorder transition and its universality class remains unclear for this model. In fact, studies of the critical exponents at the continuous phase transitions of all versions of the ALG model are still missing.

Motivated by this issue, we study the phase transition in the full lattice limit of the symmetric 2D ALG model Balladares, Henriques, and Barbosa (2007), using MC simulations (employing both Wang-Landau Wang and Landau (2001a, b) and Metropolis algorithms), as well as obtaining the thermodynamic properties of the model in the core of Husimi cacti Husimi (1950). A fully occupied lattice corresponds to the limit of infinite chemical potential, for which only the HDL and an isotropic, disordered phase (corresponding to what is identified as the gas phase in previous studies of the AGL) are expected, since the LDL phase becomes metastable already for small values of the chemical potential Balladares, Henriques, and Barbosa (2007); Oliveira, Stilck, and Barbosa (2010); Furlan, Fiore, and Barbosa (2015). Indeed, we find a single transition between the isotropic and the ordered HDL phase, whose loci, nature and universality class will be addressed in detail in the following.

The remainder of this paper is organized as follows. In section II we detail the model. In Sec. III we explain the simulation and analytic methods used. In Secs. IV and V we report our simulation results and Husimi-cacti findings, respectively. Finally in Sec. VI, we discuss our conclusions and perspectives for future work.

II Model

We consider the associating lattice gas (ALG) introduced in Ref. Henriques and Barbosa (2005) in its symmetric version Balladares, Henriques, and Barbosa (2007). The model was proposed in the context of water-like anomalies; despite its simplicity, it captures some features of liquid water, such as density Henriques and Barbosa (2005); Balladares, Henriques, and Barbosa (2007) and diffusion Szortyka et al. (2009) anomalies. The model is defined on a triangular lattice (coordination number $z=6$) in which each site can be either empty or occupied by a molecule. Each molecule has six arms, four of which are bonding arms, while the other two are non-bonding (“inert”). In the symmetric version Balladares, Henriques, and Barbosa (2007), all the bonding arms interact in the same manner, there being no distinction between proton donors and receptors. The inert arms do not interact and assume diametrically opposed positions, giving rise to three orientational states, ${\bm{\eta}}_{i}$ with $i=1,2$ or $3$. In this work we adopt a different notation from that of Refs. Henriques and Barbosa (2005) and Balladares, Henriques, and Barbosa (2007). We denote the generators of the triangular lattice by $\hat{e}_{1}$, $\hat{e}_{2}$ and $\hat{e}_{3}$, with,

[TABLE]

[TABLE]

and

[TABLE]

We use the same set of vectors to label the orientational states. Consider, for example, state 1, with bonding arms along $\pm\hat{e}_{2}$ and $\pm\hat{e}_{3}$. As illustrated in Fig. 1, we associate the vector ${\bm{\eta}}_{1}=\hat{e}_{1}$ with state 1, ${\bm{\eta}}_{2}=\hat{e}_{2}$ with state 2 and ${\bm{\eta}}_{3}=\hat{e}_{3}$ with state 3. (Thus ${\bm{\eta}}$ points along one of the nonbonding directions.)

In the ALG, interactions are restricted to nearest-neighbor (NN) pairs, so that the separations ${\bm{r}}_{i}-{\bm{r}}_{j}$ between interacting pairs again fall in the set $\{\pm\hat{e}_{j}\}$. A particle at site $k$ in state ${\bm{\eta}}_{i}$ has no interaction with its neighbors at sites ${\bm{r}_{k}}\pm\hat{e}_{i}$ since it has no bonding arms pointing toward these sites. On the other hand, for $i\neq j$ we have $|\hat{e}_{i}\cdot\hat{e}_{j}|=1/2$. Thus the interaction between a pair of particles at sites $i$ and $j$, separated by ${\bm{r}}$ (a unit vector in the set $\{\pm\hat{e}_{i},i=1,2,3\}$), can be written so:

[TABLE]

where $-\gamma$ denotes the energy associated with a NN particle pair with bonding arms pointed toward each other. Using this, the energy of the ALG model may be written:

[TABLE]

where the summation $\left\langle i,j\right\rangle$ runs over pairs of nearest neighbors, the $\sigma_{i}=0,1$ are occupation variables, and $\varepsilon>0$ represents an orientation-independent repulsive NN interaction. In the present work, we consider full occupancy ($\sigma_{i}=1,\forall i$), so that the parameter $\varepsilon$ is meaningless and $\gamma$ is the only remaining parameter in the Hamiltonian. Therefore, from here on we measure energy in units of $\gamma$, and define a dimensionless temperature $\tau=k_{B}T/\gamma$. We remark that when comparing our results with those in the literature for the full model Balladares, Henriques, and Barbosa (2007); Oliveira, Stilck, and Barbosa (2010), one should keep in mind that in these papers the particular case $\gamma/\varepsilon=2$ was considered, and the parameter $\varepsilon$ was used to define the reduced temperature and chemical potential. Thus, in these studies the reduced variables are twice those used here.

At full occupancy, each particle interacts with four of its six neighbors if all particles have the same orientational state (see Fig. 1). Any configuration such that one or more pairs of nearest neighbors have distinct states must have a higher energy. Thus the ground state is threefold degenerate, with energy per particle (in units of $\gamma$), $e=-2$. By contrast, the mean energy per particle in a configuration in which each orientation is chosen at random, independently, is $e_{random}=-4/3$. Since the gain in entropy per site is $\Delta s/k_{B}=\ln 3$, a crude estimate for the critical temperature is $\tau_{c}\approx 2/(3\ln 3)\simeq 0.6$.

The model studied here corresponds to the limit $\mu\to\infty$ of the full model, and at low temperatures it will be in the HDL phase. Varying the temperature, we expect to observe a transition from the low-temperature ordered phase, with a majority of particles in one of the orientational states, to a high-temperature disordered one, with equal populations among the three states. This isotropic phase corresponds to the gas phase in the general model. As noted in the Introduction, continuous gas-HDL phase transitions have already been observed in simulations of the nonsymmetric two-dimensional Szortyka et al. (2009) and three-dimensional Buzzano, De Stefanis, and Pretti (2008); Szortyka et al. (2010); Furlan, Almarza, and Barbosa (2016) ALG models, while a discontinuous transition was found in the mean-field approach for the symmetric model considered here Oliveira, Stilck, and Barbosa (2010). In the 2D case Szortyka et al. (2009); Oliveira, Stilck, and Barbosa (2010), the transition is observed for reduced chemical potentials $\bar{\mu}\equiv\mu/\gamma\gtrsim 1$. Szortika et al. Szortyka et al. (2009) proposed an order parameter

[TABLE]

where $n_{i}$ is the number of particles in state $i$ and $N$ is the total number of particles. Evidently, $\theta=1$ when all particles are in the same state, and $\theta=0$ for equally populated states.

In view of the orientation-dependent interaction and the orientationally ordered ground state, we can interpret the transition as one between nematic and isotropic phases. Analysis of such a transition Plischke and Bergersen (2006) suggests an alternative definition of the order parameter,

[TABLE]

Of central interest is the nature of the order-disorder transition in the model at full occupancy. (In principle, we would expect this also to describe the transition at large, but finite $\mu$.) Since the ground state is threefold degenerate, it is tempting to suppose that the transition falls in the three-state Potts class. While this may in fact be the case, we note that the energy, Eq. (5), does not have the symmetry of the Potts model, that is, it is not invariant under permutations among the states. (The ground state is permutation invariant, trivially, but due to the orientational dependence of the interactions, an arbitrary configuration is not.)

In light of the above considerations, we regard it as an open question whether the model belongs to the 3-state Potts universality class. The simulation results reported below provide some insight on this issue.

III Methods

III.1 Simulation details

We performed extensive Monte Carlo (MC) simulations of the symmetric ALG model Balladares, Henriques, and Barbosa (2007) using both Wang-Landau (WL) simulations Wang and Landau (2001a, b) for triangular lattices with lateral sizes $24\leq L\leq 72$, and Metropolis simulations Metropolis et al. (1953) for larger systems. The simulations are performed at full occupancy ($N=L^{2}$). The WL algorithm is an entropic sampling method, designed to estimate $\Omega\left(E\right)$, the number of configurations with energy $E$. Starting from $\Omega\left(E\right)=1,\forall E$, the estimates for $\Omega$ are gradually refined in a series of iterations. Each iteration generates a sequence of configurations as described below. An energy histogram $H\left[E\right]$, records the number of configurations having energy $E$. These simulations use a random initial configuration $\Gamma=\{{\bm{\eta}}_{i}\}$ with energy $E\left(\Gamma\right)$. To generate a candidate for the next configuration in the sequence ($\Gamma_{new}$), a site $i$ is chosen at random and its state is altered (${\bm{\eta}}_{i}\to{\bm{\eta}}^{\prime}_{i}$), generating configuration $\Gamma^{\prime}$. This is accepted as $\Gamma_{new}$ with probability

[TABLE]

With complementary probability, the current configuration $\Gamma$ is taken as the new configuration. The histogram and logarithm of the number of states updated so: $H\left[E\left(\Gamma_{\rm new}\right)\right]\rightarrow H\left[E\left(\Gamma_{\rm new}\right)\right]+1$ and $\ln\Omega\left[E\left(\Gamma_{new}\right)\right]\rightarrow\ln\Omega\left[E\left(\Gamma_{new}\right)\right]+\ln f$. The parameter $f$, known as the modification factor, is set to $e$, the base of natural logarithms, on the first iteration. A sequence of configurations is generated by repeating this procedure, until the histogram is flat, i.e., there is no $E$ such that $H(E)$ is less than (or greater than) $X\%$ of the mean value of the histogram over all energies. Following the usual practice, we use $X=20\%$, that is, the $80\%$ flatness criterion.

Once the flatness criterion is satisfied, the current iteration ends and a new one is started with the histogram set to zero, but with the $\Omega\left[E\right]$ carried forward from the preceding iteration. The new iteration proceeds as before except for a smaller modification factor, taken as the square root of the previous value (i.e., $\ln f\to\ln f/2$). The usual procedure is to end the simulation once $\ln f<10^{-8}$ or equivalently after 27 iterations.

As is evident from Eq. (8), WL sampling is a kind of Metropolis importance sampling with target probability distribution $P(\Gamma)\propto 1/\Omega\left[E(\Gamma)\right]$ in place of the usual Boltzmann distribution, $P(\Gamma)\propto\exp[-\beta E(\Gamma)]$. Thus if iterations extended arbitrarily long and there were no sampling noise, the estimates for the $\Omega\left(E\right)$ would converge to their true values. For further details about the convergence (or lack thereof) of WL sampling, see Belardinelli et al. (2014) and references therein.

Recent studies Caparica (2014) show that it is not necessary to use all 27 iterations of WL sampling. As the modification factor is barely greater than unity, estimates for the $\Omega\left(E\right)$ are not significantly modified in the final few iterations. We used 27 iterations for system sizes $L\leq 32$, 22 for $32<L<72$, and 20 for $L=72$.

Using the estimates for the $\Omega\left(E\right)$, the canonical average of a given thermodynamic observable ${\cal O}\left(E\right)$ at temperature $T$ can be computed via,

[TABLE]

where $\overline{{\cal O}\left(E\right)}$ represents the microcanonical average, and $\cal Z$ is the canonical partition function. Microcanonical averages are computed as simple averages of property ${\cal O}$ over all configurations with energy $E$ generated in the final iteration of the WL procedure. We compute the canonical average of the energy $E$ and its variance, and of order parameters $\theta$ and $Q$ and their second through fourth moments.

Although WL sampling yields useful results for systems with $L\leq 90$, it is not effective for larger systems. (The time to achieve a flat histogram becomes excessive.) We therefore use standard Metropolis sampling Metropolis et al. (1953) for system sizes $L=128$ and $L=256$. In these studies, we use $10^{6}$ MCS for equilibration followed by $2\times 10^{6}$ MCS to generate data. In both simulation methods, we perform an average over 60 independent realizations, starting from randomly generated initial configurations.

III.2 Mean-field approximations

Generally, mean-field approximations of a particular model correspond also to the solution of this model on a complete graph, with properly rescaled interactions Baxter (1985). Since solutions of a model on the core of a Cayley tree (a graph with no loops, in which each node has the same number of neighbors) usually correspond to the Bethe approximation of this model on a regular lattice with the same coordination number as the tree, Baxter suggested that these treelike lattices be called Bethe lattices Baxter (1985). Additional correlations are taken into account considering Husimi trees, built with clusters (polygons or polyhedrons) rather than single sites, so that short closed paths are present. Analysis of a model on the core of a Husimi tree yields its behavior on a Husimi lattice Husimi (1950).

We study the ALG model on such treelike lattices, in the limit in which each site is occupied by a molecule. Although we could begin with the six-coordinated Bethe lattice, one readily verifies that the ground state of the model on this lattice, without any closed paths, is highly degenerate. The thermodynamic behavior of the model on the Bethe lattice is qualitatively distinct from that found on the triangular lattice: there is no ordered phase. We therefore study the model on a Husimi tree built with triangles, in which three triangles meet at each site, as shown in Fig. 3. We also study the model on a more elaborate Husimi lattice built with hexagons composed of six elementary triangles (see Fig. 4). While the complete model, at finite chemical potential, was studied using similar techniques on this Husimi lattice Oliveira, Stilck, and Barbosa (2010), the limit ($1/\mu\rightarrow 0$) of interest here was not analyzed in that work. As explained in Ref. Oliveira, Stilck, and Barbosa (2010), the use of hexagons in contrast to simple triangles as elementary clusters is imperative for capturing the symmetries of the ground states of the three phases of the full model (gas, LDL, and HDL) on the Husimi lattice; in the case of full occupancy, however, a tree of triangles (Fig. 3) already captures the symmetry.

As usual, to solve a model on treelike lattices we define rooted sub-trees and their partial partition functions (ppf’s) for the possible configurations of the site at the root. For the present model, there are three possible orientations of the molecule on the root site of a sub-tree, which may be represented by the variable $\eta$ assuming the values $1,2$ or $3$, depending on the orientation of the pair of non-bonding arms, as well as the orientation of the root polygon. For example, in the tree built with triangles, a sub-tree is obtained if we remove the two hatched triangles in Fig. 3, and the root site is opposite to the edge with orientation 1. This defines the orientation of the root triangle in this case as being equal to 1. A sub-tree with an additional generation of polygons can then be built by connecting two pairs of sub-trees to a new root triangle (or five sub-trees to a new root hexagon on the second Husimi lattice). This operation leads to recursion relations for the ppf’s. The configuration of the molecule at the root site is fixed and one sums over the configurations of the molecules on the other vertices of the root polygon. Upon iteration, the ppf’s usually grow without bound, and diverge in the thermodynamic limit. One thus consider ratios of ppf’s, which often reach finite fixed-point values. Thus, a thermodynamic phase is associated to one of the stable fixed points of these ratios, and expected values of densities in the central region of the tree may be obtained at the fixed points considering the operation of attaching three pairs of sub-trees to the central triangle (or six sub-trees to the central hexagon in the second lattice). Phase transitions are therefore related to distinct fixed points competing in a region of the parameter space. Finally, a free energy is needed to determine the transition if coexistence of distinct phases occurs. To find the bulk free energy we resort to an ansatz due to Gujrati Gujrati (1995); Oliveira, Stilck, and Barbosa (2010). The details of the calculations are presented in the appendix.

IV Simulation results

In this section we report simulation results for the model at full occupation.

The critical temperature $\tau_{c}$ is obtained via finite size scaling (FSS) Fisher and Barber (1972) analysis of the susceptibilities $\chi_{\theta}$,

[TABLE]

and $\chi_{Q}$, defined analogously, and of the specific heat $c$,

[TABLE]

In Eqs. (10) and (11) the terms $\left\langle\ldots\right\rangle$ represent canonical averages, $L$ the linear size of the system and $\tau$ the temperature. (As mentioned in Sec. III, our estimates for thermal averages such as $\theta$, $\chi$ and $c$ are averages over sixty independent realizations). The expected finite-size scaling forms for the order parameter, susceptibility, and specific heat are,

[TABLE]

and similarly for $Q$,

[TABLE]

and

[TABLE]

where $\beta$, $\gamma$, and $\alpha$ are critical exponents as usually defined, and $t\equiv(\tau-\tau_{c})/\tau_{c}$ is the reduced temperature. Additionally we determined Binder’s fourth-order cumulant of the order parameters Binder (1981),

[TABLE]

$U_{4,Q}$ is defined analogously. At the critical point, $U_{4}$ tends to a universal value, characteristic of the universality class, system shape and boundary conditions employed Selke (2007).

IV.1 Determining $\tau_{c}$

The results for the order-parameters $\theta$ and $Q$, plotted in Fig. 5, suggest that a continuous phase transition occurs at some temperature in the interval $0.45<\tau<0.5$. The lower inset shows the fraction of particles $\rho_{i}$ in the majority, minority and intermediate states, for system size $L=32$, illustrating a continuous variation between ground state (all particles in the same state) and nearly equal populations. Relative uncertainties are plotted in the upper inset. As expected, they are largest in the critical region. In all cases, the relative uncertainty in $\theta$ is less than 3%; for system sizes other than $L=64$, it is $<2$%.

The susceptibilities $\chi_{\theta}$ and $\chi_{Q}$ are plotted in Fig. 6. Although these quantities exhibit similar behaviors, there are slight differences in the critical region; the differences are more evident for larger systems. The discrepancies between the susceptibilities are related to their structures. While $Q$ (see Eq. 7) takes into account the densities of all three states, $\theta$ only involves the majority density (see Eq. 6). This may be why the uncertainties (inset Fig. 6) in $\chi_{Q}$ are approximately one-third those in $\chi_{\theta}$. For the sizes shown in Fig. 6, the maximum relative uncertainty $\chi_{\theta}$ in is about $6\%$.

The specific heat $c$, shown in Fig. 7, exhibits behavior consistent with that of the order parameters and susceptibilities. Relative uncertainties in $c$ are smaller than $4\%$ for $48\leq L\leq 74$ and smaller than $2\%$ for $L<48$.

The Binder cumulants of the order parameters are shown in Fig. 8.

The crossings of the Binder cumulants for $\theta$ and $Q$ provide the estimates $\tau^{(Q)}_{c}=0.476(4)$ and $\tau^{(\theta)}_{c}=0.476(5)$, respectively.

Using our results for $c$, $\chi_{\theta}$, $\chi_{Q}$, and the cumulants, we estimate the critical temperature, $\tau_{c}$. For $c$ and the $\chi$’s, we define a size-dependent pseudocritical temperature as the temperature associated with the maximum value. Pseudocritical temperatues for the cumulants are identifying as the crossing temperatures of $U_{4}$ between: 1) the smallest system size studied, $L=24$, and the others ($L=32,\ldots,L=256$) and 2) the crossing temperatures between a given system size $L$ and the next system size, for example, $L=24$ with $L=32$, $L=32$ with $L=40$ and so on. The results from the adjacent sizes crossings are represented as $U^{\prime}_{4}$.

The pseudocritical temperatures are plotted versus $1/L$ in Fig. 9. All six sets of pseudocritical temperatures appear to converge to similar values as $L\to\infty$. Of note is the relative insensitivity to system size of the pseudocritical temperatures derived from cumulant crossings. The crossings between Binder cumulant of two adjacent system sizes for the order parameter $\theta$ suffer from large uncertainties for $L\geq 64$, affecting the estimate of the pseudocritical temperature. For this reason, we disregard this property in the calculation of the critical temperature.

The resulting estimates for $\tau_{c}$ are listed in Table 1. From the six estimates for the pseudocritical temperatures we derive the global estimate $\tau_{c}=0.4763(1)$. The global estimate of $\tau_{c}$ was obtained through a weighted average of the properties with weights $1/\sigma^{2}$, where $\sigma$ represents the uncertainty of each quantity.

IV.2 Critical exponents

Extrapolating the cumulant crossings to infinite size yields two estimates for the critical cumulant: $U^{\ast}_{4,\theta}=0.604(2)$ and $U^{\ast}_{4,Q}=0.606(1)$. Since there can only be one value, a more conservative estimate is $U^{*}_{4}=0.605(2)$. It is of interest to compare this to the critical cumulant of the two-dimensional, three-state Potts model. The best available estimate of the latter, to our knowledge, is Tomé and Petri (2002) $U^{*}_{\rm Potts}=0.61$. (In Ref. Tomé and Petri (2002), no uncertainty estimate is provided; we may assume it is on the order of 0.01.) In efforts to improve on this estimate we simulated the 3-state Potts model on $L\times L$ square lattices with periodic boundaries. The results for $U_{4,{\rm Potts}}$, shown in Fig. 8 (lower-right insert), yield $U^{\ast}_{4}=0.6124(8)$ for the critical cumulant. Although strictly speaking incompatible with our estimate for the associating lattice gas, the two estimates are rather similar, leaving open the possibility of a common universality class for the two models.

We estimate the critical exponent ratios $\beta/\nu$ and $\gamma/\nu$ using fits to the data for $\theta$, $Q$, and the associated susceptibilities $versus$ system size (see Fig. 10).

From the linear fits shown in Fig. 10 (left) and including the effect of the uncertainty in $\tau_{c}$ we obtain $\beta_{\theta}/\nu=0.130(3)$ and $\beta_{Q}/\nu=0.126(1)$ The inset shows that the residuals are much smaller than the uncertainties, indicating that adding additional terms to the fitting function, in the form of corrections to scaling, would not yield better estimates. The value of this ratio for the two-dimensional, three-state Potts model is Wu (1982); Nagai, Okamoto, and Janke (2013); Caparica (2014) $\beta_{\rm Potts}/\nu_{\rm Potts}=2/15=0.1333\ldots$. In the worst case ($\beta_{Q}$), the discrepancy is about $4\%$. Linear fits to $\ln\chi_{\theta}$ and $\ln\chi_{Q}$ as functions of $\ln L$ (see Fig. 10(right)) yield $\gamma_{\theta}/\nu=1.70(2)$ and $\gamma_{Q}/\nu=1.73(9)$. Our estimate for $\gamma_{\theta}$ differs from the value of the three-state Potts model, $\gamma_{\rm Potts}/\nu_{\rm Potts}=26/15=1.733\ldots$, by $1.7\%$. Analysis of the residuals of fits to the susceptibilities again indicates that adding further terms to the fitting function is not necessary.

Different from the order parameter and susceptibility, a linear fit to $\ln c$ as a function of $\ln L$ does not yield a good description of the data. We therefore fit the data with $c=aL^{\alpha/\nu}+c_{0}$ with $a$, $\alpha/\nu$ and $c_{0}$ as adjustable parameters. This form is capable of fitting the simulation data, as can be seen in Fig. 11. The inset shows that the residuals are smaller than the uncertainties, and do not exhibit a systematic tendency. A least-squares procedure yields $a=2.4(6)$, $\alpha/\nu=0.38(3)$ and $c_{0}=-3\pm 1$. The exponent ratio agrees to within uncertainty with that of the three-state Potts model, $\gamma/\nu_{\rm Potts}=0.4$.

Given the scaling relation $\alpha+2\beta+\gamma=2$, we have,

[TABLE]

where exponents with primes denote the corresponding exponents divided by $\nu$, determined in the finite-size scaling analysis discussed above.

Using the values for $\beta^{\prime}$ and $\gamma^{\prime}$ obtained using the order parameter $\theta$, equation 16 yields $\nu_{\theta}=0.85(5)$, while the values associated with order parameter $Q$ furnish $\nu_{Q}=0.84(1)$, for an average of $\nu=0.85(1)$, close to the Potts model value, $\nu_{\rm Potts}=5/6=0.833\ldots$.

As a further test regarding the universality class, we perform a data collapse using the Potts critical exponents, as shown in Fig. 12. The quantities $\theta,Q,\chi_{\theta}$ and $\chi_{Q}$ exhibit a good collapse using these exponents. We observe that the specific heat exhibits the poorest collapse.

The results for the exponents $\alpha,\beta$, $\gamma$ and $\nu$, as well as for the Binder cumulant strongly suggest that the phase transition studied here belongs to the three-state Potts model universality class. Deviations of the critical exponents from the Potts class values may be due to the flatness criterion, limited sample size, and restricted system sizes.

V Results on Husimi lattices

The thermodynamic behavior of the model on the core of both treelike lattices will be presented here, after a brief description of the calculations which lead to it. Some more details of these calculations may be found in the appendix. As mentioned in Sec. III, recursion relations for the partial partition functions on a rooted sub-tree are obtained, and ratios of these fpp’s usually remain finite as the recursion relations are iterated. The stable (real and positive) fixed points attained iterating the recursion relations for the ratios correspond to the thermodynamic limit. The behavior on the complete tree may then be found connecting sub-trees to a central polygon. Expected values of densities at this central polygon may then be calculated, and may be seen as approximations to the results on the triangular lattice.

Two stable fixed points are found. At higher temperatures, only a disordered isotropic fixed point is stable, for which the densities of molecules in the three possible orientations are equal. At low temperatures, three ordered nematic fixed points are stable, where the density of molecules in one of the three orientations is larger than the ones in the other two orientations. The stability of these fixed points may be studied using the Jacobian of the recursion relations for the ratios. It is found that there is an interval of temperatures where both kinds of fixed point are stable, signalling that the transition between the phases is discontinuous. Thus, to find the coexistence temperature the free energy of both phases should be equated. It should be remembered that this free energy corresponds to the model on the core of the tree; the free energy of the whole tree is dominated by the surface Baxter (1985); Gujrati (1995). We find a coexistence temperature $\tau_{c}=0.51403$ for the triangle tree and $\tau_{c}=0.51207$ for the one with hexagons.

The nematic order parameter $Q$ is shown in Fig. 13 as a function of the temperature for both treelike lattices. We note that the transition temperature is reduced, and the jump in the order parameter becomes smaller as we move from the triangle tree to that with hexagons. Also, the coexistence temperature decreases and the temperature interval in which both fixed points are stable becomes narrower. These results are consistent, as the transition temperature on the triangular lattices estimated from the simulations is still lower and the transition is continuous.

Finally, let us discuss how transition we study here fits into the general phase diagram of the model for finite chemical potential. The phase diagram of the complete model on the hexagon-Husimi tree is shown in Fig. 14, in the $\tau\times{\bar{\mu}}=\mu/\gamma$ plane. In this diagram we see that the gas-HDL coexistence line, as well as the gas and HDL spinodals, approach vertical asymptotes, and are very close to them already for $\bar{\mu}\gtrsim 10$, whose $\tau$ values are very close to the ones presented above for the solution in the full-occupancy case. Hence, our results for the infinite chemical potential limit correspond to the final point of the coexistence line between the HDL and gas phases. For finite $\mu$ the coexisting phases differ in the order parameter $Q$, and in the particle density $\rho$; the density is lower and $Q$ vanishes in the gas phase. In the full-occupancy limit we consider here, $Q$ still vanishes in the gas phase and is finite in the HDL phase, but $\rho=1$ in both phases.

Figure 15 shows the discontinuity of the nematic order parameter at the gas-HDL coexistence line as a function of the reduced chemical potential ${\bar{\mu}}$. We note that $\Delta Q$ becomes smaller as the chemical potential grows, and reaches its minimum value in the limit studied here. We can imagine three possible scenarios for this part of the phase diagram on the triangular lattice. One possibility is that in the limit $1/\mu\to 0$ the gas-HDL coexistence line ends at a critical point. Another would be that the transition becomes continuous already at some finite value of $\mu$, and thus the coexistence line would end at a tricritical point. Finally, the whole transition may be continuous. In this case, if the other transition lines (LDL-HDL and gas-LDL), which meet with the gas-HDL transition line at a triple point on the hexagon tree (see Fig. 14), remain discontinuous, this point would become a critical end point.

VI Conclusions

We investigate the gas-HDL phase transition of the ALG modelHenriques and Barbosa (2005) in its symmetric versionBalladares, Henriques, and Barbosa (2007) in the limit $1/\mu\to 0$ (full occupation) by Monte Carlo simulation and solution on Husimi lattices mean-field approximations for the solution on the triangular lattice). The simulations reveal a continuous transition at $\tau_{c}=0.4763(1)$. The critical temperature is estimated using finite-size scaling estimates derived from the order parameters, their susceptibilities, the specific heat, and the fourth-order Binder cumulants. In addition, we derive an estimate of $U_{4}^{*}$ for the three-state Potts model on a square lattice: $U_{4}^{*}=0.6124(8)$. This value is close to, but more accurate than, one reported previously [$U_{4}^{*}=0.61$] Tomé and Petri (2002). The critical exponents reported in Sec. IV also support a 3-state Potts model universality class.

The mean-field calculations are performed on Husimi trees constructed with triangles and hexagons. Both approaches revel a discontinuous phase transition with coexistence temperatures somewhat higher than the simulation estimate, viz., $\tau_{c}=0.51403$ (triangles) and $\tau_{c}=0.51207$ (hexagons). In the second case the difference in relation to the simulation is $\lesssim 7\%$. We note that mean-field calculations for the three-state Potts model in two-dimensions also yield discontinuous transitions Mittag and Stephen (1974), which may explain the difference found here for the ALG model.

The complete (mean-field) phase diagram shows a gas-HDL coexistence line which extends from $\bar{\mu}\approx 1$ to $\bar{\mu}\to\infty$, along which a reduction of the discontinuity of the nematic order parameter $Q$ is observed as $\mu$ increases. The existence of the discontinuous line gives origin to three possible scenarios that connect the results of simulations and mean-field. The first is that at the limit $1/\mu\to\infty$ the line ends in a critical point, the second is a change in the order of transition at a finite chemical potential, and the third raises the possibility of the whole transition be continuous for the triangular lattice. We remark that the ALG model with finite $\bar{\mu}$ is somewhat similar to a diluted three-state Potts system, for which a continuous transition is also observed Nienhuis et al. (1979); Qian, Deng, and Blöte (2005). This suggests that the entire gas-HDL line might be continuous. Another indication of this is the fact that it is continuous in the original 2D ALG model Szortyka et al. (2009), as well as in the 3D version Buzzano, De Stefanis, and Pretti (2008); Szortyka et al. (2010). In future work, we intend to perform Wang-Landau sampling, as well as transfer matrix calculations in order to obtain the complete phase diagram and identify which of the three scenarios mentioned above is correct.

VII Acknowledgments

The authors acknowledge financial support from the CNPq, CAPES and Fapemig, Brazilian agencies. We thank to M.C. Barbosa and W. Selke for helpful discussions and also Statistical physics Laboratory of Universidade Federal de Minas Gerais for the computational support.

Appendix A Solution of the model on Husimi trees

We show in some detail the calculations which lead to the thermodynamic properties on a Husimi tree, particularly the one built with triangles, as shown in Fig. 3. As described above, we start by defining partial partition functions (ppf’s) for a sub-tree with a fixed configuration of the molecule on the root site. We thus define nine ppf’s $g_{i,\eta}$ of rooted sub-trees, where the first index assumes the values $1,2$ or $3$ of the direction of the edge opposite to the root site of the triangle. On the sub-tree obtained removing the hatched triangles in Fig. 3, this index is 1. The second index specifies the orientation $\eta$ of the non-bonding arms of the molecule at the root.

The recursion relations are obtained by considering the operation of building a sub-tree with an additional generation $(M+1)$ by connecting two pairs of sub-trees (with $M$ generations each) to the vertices which are opposite to the root vertex of the new root triangle. We call $j_{i}$ the orientation of the molecule placed on the vertex opposite to the edge of orientation $i$ of the root triangle. The ppf’s of a sub-tree with $M+1$ generations are thus obtained through the recursion relations:

[TABLE]

The function $W_{\{j_{1},j_{2},j_{3}\}}$ is the statistical weight associated to the edges of the root triangle, defined as $W=\exp[2n_{b}(j_{1},j_{2},j_{3})/\tau]$, where $n_{b}$ is the number of edges of the triangle with no non-bonding arms on them and $\tau$ is the temperature. Thus:

[TABLE]

where $\delta_{i,j}$ is the Kronecker delta. In this expression $n_{b,i}$ stands for the number of edges in direction $i$ with a bond on them, assuming the values 0 or 1, so that, as expected, $n_{b}\in[0,3]$. Thus, the system of 9 recursion relations [Eqs. (17)] allows us to obtain the ppf’s of sub-trees with an arbitrary number of generations of triangles and, most importantly, the thermodynamic limit (when $M\to\infty$). Usually, the ppf’s diverge in this limit, so it is convenient to define ratios of ppf’s

[TABLE]

These ratios generally converge to finite values in the thermodynamic limit; six of them are independent, since $\sum_{j}R_{i,j}=1$. Therefore, the thermodynamic behavior of the model will be determined by the stable fixed points of the recursion relations for the ratios $R_{i,j}$, which are obtained from Eqs. (17).

We find fixed points with two different symmetries. In the isotropic fixed point the values of the 9 ratios, represented as a $3\times 3$ matrix, are:

[TABLE]

where the parameter $x$ is a function of the temperature $\tau$, whose exact expression can be easily found with the help of an algebra software, but is too long to be presented here. The three nematic fixed points are:

[TABLE]

Again the parameters $x_{i}$, with $i=1,2,3$, are functions of the temperature. We note that the isotropic fixed point is stable at high enough temperatures, while at lower temperatures one of the nematic fixed points is reached. We will see below that in the isotropic phase the orientation of the non-bonding arms of the molecules has no preferred direction, while in the nematic phase one of the three directions is more probable than the other two. As expected for a discontinuous transition, the isotropic and nematic fixed points are both stable in an interval of temperatures around the coexistence temperature. The region of stability of a given fixed point may be found from the largest eigenvalue of the Jacobian of the recursion relations for the ratios:

[TABLE]

where $R_{i}^{\prime}$ and $R_{i}$ denotes the ratios in generations $M+1$ and $M$, respectively, and the derivative is calculated at the fixed point whose stability is being considered. We find that the isotropic fixed point is stable for $\tau\geq 0.47299$, while the nematic ones for $\tau\leq 0.51970$. Hence, for $\tau\in[0.47299,0.51970]$ both fixed points are stable, so that the isotropic and nematic phases coexist in this temperature interval. This signals a discontinuous transition between the phases and, then, the coexistence temperature can be determined as the one at which the free energy in the bulk of the tree is equal for both phases.

Following the ansatz proposed by GujratiGujrati (1995), we assume that the contribution to the free energy per triangle is different for triangles located at the surface of the tree ($\phi_{s}$) and the ones in the bulk ($\phi_{b}$). In a tree with $M$ generations of triangles, the number of triangles in the bulk is $N_{b}=1+6+6\cdot 4+...+6\cdot 4^{M-2}$, while the number of triangles on the surface is $N_{s}=4\cdot 4^{M-1}$. Therefore, if $\varphi_{M}$ is the free energy of a tree with $M$ generations of triangles, we find that $\varphi_{M+1}-4\varphi_{M}=3\phi_{b}$, and since $\varphi_{M}=-k_{B}T\ln Y_{M}$, we have that:

[TABLE]

where $Y_{M}$ is the partition function on the whole tree, which can be obtained by connecting six sub-trees to the central triangle. This procedure leads to:

[TABLE]

In the thermodynamic limit $M\to\infty$, using the recursion relations [Eqs. (17)], the ratios [Eqs. (19)] at the fixed points, and this expression for the partition function [Eq. (26)], the bulk free energy per triangle may be obtained through Eq. (25). This yields the coexistence temperature $\tau_{c}=0.51403$. Therefore, we indeed have a discontinuous transition between the isotropic and nematic phases in the solution of the model on this Husimi lattice.

To further confirm this, we calculate the order parameter $Q$ as defined in Eq. 7. From the partition function (26) it is simple to obtain the mean value of the number of molecules with orientation $i$ of non-bonding arms in the central triangle:

[TABLE]

where this number has been normalized so that $0\leq\langle n_{i}\rangle\leq 1$. If we divide the numerator and the denominator by $\prod_{i=1}^{3}(g^{(M)}_{i,1}g^{(M)}_{i,2}g^{(M)}_{i,3})^{2}$ we may express this mean values in terms of the ratios, which in the thermodynamic limit assume their fixed point values. The result is:

[TABLE]

where:

[TABLE]

As expected, in the isotropic gas phase, one has $\langle n_{1}\rangle=\langle n_{2}\rangle=\langle n_{3}\rangle$, so that the order parameter is $Q_{gas}=0$. On the other hand, in the nematic HDL phase we find a non-vanishing $Q_{HDL}$, whose variation with the temperature near the coexistence is displayed in Fig. 13. At the coexistence point the discontinuity in $Q$ is $\Delta Q=0.67056$.

We performed similar calculations for the tree built with hexagons (shown in Fig. 4), which are actually a particular case of the ones presented in Oliveira, Stilck, and Barbosa (2010), when all lattice sites are occupied by molecules. For brevity, we will not detail the calculations here, and present the main results only. Since the tree in this case is built with larger clusters, hexagons with one central site and six sites at the border, it is expected that such a calculation should lead to results which are closer to those on the triangular lattice. Again, we find a discontinuous transition, at a somewhat lower temperature, $\tau_{c}=0.51207$. The discontinuity in the nematic order parameter at the coexistence is $\Delta Q=0.6216$, and for $\tau\in[0.48569,0.51575]$ both fixed points are stable. In this case, since there are effectively four sites per hexagon on the treelike lattice ($6/2$ at the hexagon’s border and the central one), there are more than one way to define the order parameter at the central hexagon. One could, for example, consider just the central site or a mean value of the order parameter over all sites of the central hexagon. We have done both calculations and found out that these values are quite close. The deviation is maximum at the coexistence temperature and the relative difference is about 0.25% there.