# Determination of the order of phase transitions in Potts model by the graph-weight approach

## 1 Introduction

One of the less trivial questions in the study of phase-transitions is the determination of the order of transition. A large amount of work, especially on the Potts  model, ranging from the exact solutions  to different kinds of Monte Carlo (MC) simulations  (to cite only some of them), has been done in order to establish reliable and at the same time applicable criteria for distinction between the first- and second-order phase transition.

The coexistence of two different phases at the first-order phase transition point and its absence at the second-order phase transition point is the main physical fact on which all of the above methods (including the present one) are based. Performing the MC simulations on finite models in order to obtain the probability distribution of some quantity (such as the energy or the order parameter) which has different values in different phases, is the usual way to investigate these phase transitions numerically. In the present paper we show that the quantity which we define as the graph weight also has different values in the coexisting phases at the first-order transition point. In Section 3 it is shown how the graph weight is related to the energy and the free energy of the system. The coexistence of phases with different energy at the first-order transition point then implies the coexistence of graphs with different weights. It can be identified by the two peaks in the graph-weight probability distribution.

In order to test the graph approach, we analyse the Potts model in two cases: the mean-field (MF) case defined by the interactions of equal strength between all particles of the model and the one-dimensional $`(1d)`$ model with long-range (LR) interactions decaying with distance $`r`$ as $`r^{(1+\sigma )}`$, $`\sigma >0`$. In both cases two regimes (of the first- and second-order transition) are present. The regimes are separated by a point $`q_c`$ in the MF case and by a line $`(q_c,\sigma _c)`$ in the LR case.

The present approach has also the advantage of dealing directly with the non-integer values of the Potts states $`q`$. (Note that in reference  the non-integer $`q`$ values were not obtained by a direct calculation, but by a histogram-like extrapolation from the results for integer-$`q`$.) This may be of interest in studying the threshold $`q_c`$, separating the first- from the second-order transition regime, which is not always an integer. For example, in three-dimensional short-range (SR) Potts model, $`q_c`$ is not an integer , and in $`(1d)`$ Potts model with LR interaction, the onset of the first order transition at $`q_c`$ seems to depend continuously on the interaction-range parameter $`\sigma `$ .

The plan of the paper is as follows. In the next section we define the model and introduce the corresponding graph expansion. The basic steps of MC algorithm are also explained in that section. In the third section we discuss separately the results for the MF and the LR interaction case. In the last section we summarize and discuss the advantages and open problems connected with the graph approach.

## 2 Model and method

The graph representation  was already used in MC investigation of critical behaviour of Potts models . In the present approach we follow the line of reasoning given in reference , but with the basic difference that we study directly the probability distribution of graph weights and not the cluster probability distribution.

We begin by rewriting the model in the graph language. The reduced Hamiltonian of the model, with periodic boundary conditions, has the form

$$\frac{H}{k_BT}=\underset{i=1}{\overset{N1}{}}\underset{j=1}{\overset{Ni}{}}K_j\delta (s_i,s_{i+j}).$$
(1)   
where $`s_i=0,1,\mathrm{},q1`$ is the Potts particle placed on the $`i`$-th site of the chain, while $`\delta `$ denotes the Kronecker symbol. $`K_j`$’s denote interactions between two particles at distance $`j`$. Due to boundary conditions in the LR case, each $`K_j`$ has two contributions involving interactions at distances $`j`$ and $`(Nj)`$. In the MF case, the interactions are unique and the boundary conditions have no meaning. We use the system of units where $`K_j`$ is proportional to the inverse temperature.

The usual substitution, $`\mathrm{exp}[K_j\delta (s_i,s_{i+j})]=1+v_j\delta (s_i,s_{i+j})`$, with $`v_j=\mathrm{exp}(K_j)1`$, leads to the partition function of the form

$$Z_N=\underset{l=1}{\overset{N}{}}\underset{s_l=0}{\overset{q1}{}}\mathrm{exp}(H/k_BT)=\underset{l=1}{\overset{N}{}}\underset{s_l=0}{\overset{q1}{}}\underset{i=1}{\overset{N1}{}}\underset{j=1}{\overset{Ni}{}}[1+v_j\delta (s_i,s_{i+j})].$$
(2)   
It is straightforward to establish one-to-one correspondence between each member of the r.h.s. of the above equation and a graph on the chain consisting of $`N`$ particles. Each square bracket can contribute to the above product in two ways: it gives $`1`$ when there is no connection between the $`i`$-th and the $`i+j`$-th particle (inactive link), and $`v_j\delta (s_i,s_{i+j})`$ when there is a direct connection (active link) between these two particles. The analytical expression attached to each graph $`𝒢`$ will be called the graph weight, $`W_N(𝒢)`$. The explicit expression for $`W_N(𝒢)`$ follows from the r.h.s. of eq. (2). We can describe every graph consisting of $`c(𝒢)`$ clusters. By cluster we mean a set of particles interconnected by any type of active links and disconnected from other particles. Single particles are considered as one-particle clusters. The products of $`\delta `$-functions from the expression for each graph will delete all except $`c`$ $`q`$-summations on the r.h.s. of (2), so that the analytical expression for $`W_N(𝒢)`$ is

$$W_N(𝒢)=v_1^{b_1(𝒢)}v_2^{b_2(𝒢)}\mathrm{}v_{N1}^{b_{N1}(𝒢)}q^{c(𝒢)}.$$
(3)   
The symbols $`b_j(𝒢)`$ denote the total number of active links of type $`j`$.

In this way, the summation over all $`2^{N(N1)/2}`$ possible graphs between $`N`$ particles corresponds to the partition function of the model

$$Z_N=\underset{\mathrm{all}𝒢}{}W_N(𝒢)=\underset{W}{}𝒩_{N,W}W_N.$$
(4)   
The $`𝒩_{N,W}`$ denotes the number of different graphs with the same weight $`W_N`$. The above equation is thus the graph analogy of the more familiar expression

$$Z_N=\underset{E}{}𝒩_{N,E}e^{E_N/T},$$
(5)   
where the summation runs over all different energies $`E_N`$ of the system, while $`𝒩_{N,E}`$ is their degeneracy. Indeed, behind the analogy, there is a connection between the average number of clusters and the active links and average of energy and free energy of the system. The derivative of the partition function, given by eqs. (3), (4) and (5), over temperature  leads to

$$E_N=T^2\underset{j=1}{\overset{N1}{}}\frac{v_j+1}{v_j}b_j\frac{K_j}{T},$$
(6)   
while the derivative over $`q`$ leads to

$$\frac{\mathrm{ln}𝒩_{N,E}(q)}{q}\frac{1}{T}\frac{E_N}{q}=\frac{1}{q}c.$$
(7)   
The coexistence of two phases characterized by different values of energy corresponds, by the above relations, to the coexistence of graphs with different weights.

According to (4), one introduces the graph-weight probability distribution $`P_N`$

$$P_N=\frac{𝒩_{N,W}W_N}{Z_N}.$$
(8)   
Numerically, $`P_N`$ is obtained by a simple MC simulation of Metropolis type. The basic steps are:

(a) Pick at random one link in the graph $`𝒢`$ and change its status from active to inactive or vice versa. The resulting graph is called $`𝒢^{}`$.   
(b) Compare the random number $`0<r1`$ with the ratio $`W_N(𝒢^{})/W_N(𝒢)`$. If the ratio is smaller than $`r`$ keep $`𝒢`$, otherwise save $`𝒢^{}`$ as $`𝒢`$.   
(c) Return to (a).

At each step, the number of clusters has to be counted. Since we deal with a model with interactions of infinite range, the determination of clusters could not be performed by determination of their boundaries, but one has to check for possible connections with all the particles of the considered cluster. To do so, during the simulations, we have to keep record of the cluster structure of the system, i.e. which particle belongs to which cluster. We start with a simple and known cluster configuration. In each of the following MC steps, when one link is changed, we examine whether this has produced a change in cluster structure: if the link is deactivated - whether the corresponding cluster is split in two or not; if the link is activated - whether two clusters get connected or not.

Counting graphs by their weights through a large number of steps gives the unnormalized graph-weight probability distribution. Unfortunately, the shortcoming of the graph approach, is that it takes a much longer time to obtain the same precision of results compared to the simulation techniques that can be applied only for integer $`q`$ (like algorithm by Luijten and Blöte ). On the other hand, the comparison of the present approach with the simple Metropolis single spin flip algorithm done for the model of N particles with integer $`q`$ shows that $`10^4`$ flips per link of everyone of $`N(N1)/2`$ links in graph simulations give approximately the same precision as $`10^6`$ one-particle flips per particle in the spin simulations. For the chain of 400 sites considered here, the graph approach requires about $`25`$ times larger CPU time than the Metropolis algorithm on spins.

## 3 Results

The simulations were applied to the two mentioned cases of the Potts model with $`N=100,150,200,250,300,350`$ and $`400`$ particles. In final extrapolations only the data for $`N200`$ were used. The precision of simulations is determined by performing $`10^4`$ flips per link.

### 3.1 Mean-field case

The MF case of Hamiltonian (1) is given by taking the equal strength of interactions among all the particles, i.e. $`K_j=K/N`$, where $`K`$ denotes the inverse temperature. The exact work of Kihara et al.  puts in evidence that the MF Potts model has a second-order phase transition for $`qq_c=2`$ and a first-order phase transition for $`q>2`$. It thus qualifies as a good example for testing the present approach. In our simulation we start with the case $`q=3`$. In continuation, we consider lower values of $`q`$, closer to the threshold $`q_c`$ in order to examine how efficient our approach can be in detecting weak first-order transitions.

The shape of the finite-$`N`$ graph-weight probability distribution $`P_N`$ depends on temperature in a similar way as the energy-probability distribution does. Far from the transition temperature, $`P_N`$ has a gaussian form. At the transition, the shape of $`P_N`$ depends on the order of the transition. For second-order transitions, the shape has a non-gaussian form, but it still has only one maximum. For first-order transitions, $`P_N`$ has two maxima which transform into two $`\delta `$-functions in the thermodynamic limit. We define the temperatures $`T_N`$ where the two maxima in $`P_N`$ are of equal height. The corresponding positions on the graph-weight axis are $`W_N^\mathrm{o}`$ (for the ordered phase) and $`W_N^{\mathrm{do}}`$ (for the disordered phase). When the two-peak structure of $`P_N`$ becomes more pronounced with increasing $`N`$, one may conclude that the transition is of the first order in the thermodynamic limit with two coexisting phases: the one described by the graph weight $`W_N^\mathrm{o}W^\mathrm{o}`$, stable below the transition temperature $`T_NT_t`$, and the other one described by $`W_N^{\mathrm{do}}W^{\mathrm{do}}`$ stable above $`T_t`$.

In Figure 1 are presented the results of simulations for $`q=3`$. It shows dependence of $`P_N`$ on $`w_N\mathrm{ln}W_N/N\mathrm{ln}N`$ at temperatures $`T_N`$. The peaks emerge for $`200N400`$ and become more pronounced with increasing $`N`$. By the above arguments, the $`N`$-dependence of $`P_N`$ confirms the existence of the first-order phase transition, in agreement with the exact solution .

The behaviour for smaller $`N`$ considered ($`N<200`$) points out that only one maximum in $`P_N`$ may have different origins: either the transition is of the second order in the thermodynamic limit and we observe its finite-$`N`$ behaviour, or the transition is of the first order in the thermodynamic limit, but the correlation length, although finite, is comparable with the system size used and we are not able to see the coexistence of the two phases. Arbitrarily close to the threshold which separates the first- from the second-order transitions, the correlation length is expected to become arbitrarily large and finally whatever system size is used, it will always be too small to show the first-order character of the transition. Such a situation is expected in the MF case when $`q`$ approaches 2. The simulations on $`q=2.8`$ model performed at the corresponding temperatures defined as $`T_N`$ give the distributions similar to those shown in Fig. 1. On the contrary, simulations performed for $`q=2.1`$ and $`2.5`$ do not show the two-peak structure in the graph-weight probability distribution for the considered system sizes.

To some extent, the above presented simulations can be compared to the exact solution. Since $`𝒩_{N,E}(q)`$ and $`E_N(q)`$ can be exactly calculated for the transition temperature $`T_{\mathrm{MF}}^1=\mathrm{\hspace{0.17em}2}[(q1)/(q2)]\mathrm{ln}(q1)`$, , it is easy to obtain from eqs. (6) and (7) the analytic expressions for the positions $`w_N^{\mathrm{o},\mathrm{do}}`$ of the two peaks at temperature $`T_{\mathrm{MF}}`$ in the large $`N`$ limit:

$`w_N^\mathrm{o}`$ $`=`$ $`{\displaystyle \frac{(q1)^3+1}{q^2(q2)}}\mathrm{ln}(q1)+{\displaystyle \frac{a_q^\mathrm{o}}{\mathrm{ln}N}}+𝒪({\displaystyle \frac{1}{N}}),`$   
$`w_N^{\mathrm{do}}`$ $`=`$ $`{\displaystyle \frac{q1}{q(q2)}}\mathrm{ln}(q1)+{\displaystyle \frac{a_q^{\mathrm{do}}}{\mathrm{ln}N}}+𝒪({\displaystyle \frac{1}{N}}),`$ (9)   
where   
$`a_q^\mathrm{o}={\displaystyle \frac{(q1)^3+1}{q^2(q2)}}\mathrm{ln}(q1)\mathrm{ln}\left[2{\displaystyle \frac{q1}{q2}}\mathrm{ln}(q1)\right]+{\displaystyle \frac{\mathrm{ln}q}{q1}}\left[1{\displaystyle \frac{\mathrm{ln}(q1)}{q(q2)}}\right],`$   
$`a_q^{\mathrm{do}}={\displaystyle \frac{q1}{q(q2)}}\mathrm{ln}(q1)\mathrm{ln}\left[2{\displaystyle \frac{q1}{q2}}\mathrm{ln}(q1)\right]+\mathrm{ln}q\left[1{\displaystyle \frac{(q1)\mathrm{ln}(q1)}{q(q2)}}\right].`$

In Table 1 we compare the data for $`q=3`$ obtained by MC simulations at $`T_{N=400}=0.3566`$, with those calculated from (3.1) for finite $`N=400`$ and for $`N\mathrm{}`$.

The values of $`w_N^{\mathrm{o},\mathrm{do}}`$ cited in Table 1 show that the simulated results fit well the exact values for finite $`N=400`$. Relatively high difference ($`24\%`$ for $`w_N^\mathrm{o}`$ and $`32\%`$ for $`w_N^{\mathrm{do}}`$) compared to the exact $`N\mathrm{}`$ values comes from the slow $`1/\mathrm{ln}N`$ convergence.

In Table 2 we present the corresponding data for $`T_N`$ (obtained in simulations) in function of size. Since for large $`N`$ the dependence of $`T_N`$’s is approximately linear in $`N^1`$, we use the extrapolation form $`T_N=T_{\mathrm{MF},t}+b/N`$. The coefficients $`T_{\mathrm{MF},t}`$ and $`b`$ were obtained in least-squares approximation (LSA). By taking into account only the data for sizes $`200N400`$, one obtains $`T_{\mathrm{MF},t}(q=2.8)=0.376`$ and $`T_{\mathrm{MF},t}(q=3.0)=0.360`$. In both cases, the difference between the extrapolated values and the exact ones $`T_{\mathrm{MF}}(q=2.8)0.378`$ and $`T_{\mathrm{MF}}(q=3.0)0.361`$, is less than $`1\%`$.

### 3.2 One-dimensional long-range case

In the case of one-dimensional Potts model with interactions decaying as $`1/r^{1+\sigma }`$, the interaction strength $`K_j`$ is given by $`K[1/j^{1+\sigma }+1/(Nj)^{1+\sigma }]`$, where $`K`$ denotes the inverse temperature. This model has a phase transition at finite temperature  for all $`q`$ and $`0<\sigma 1`$. Our recent MC simulations of the energy probability distribution on $`q=3`$ and $`q=5`$ models  have shown that the order of the transition depends on $`q`$ and $`\sigma `$. For the same fixed value $`q=3`$, it was shown that the transition changes from the first- to the second-order one with increasing $`\sigma `$.

Again, according to the relations (6) and (7), one can interpret the meaning of the graphs through their connection with the energy and the free energy of the system.

The purpose of simulations presented in this subsection is twofold: first, we wish to test the above described graph approach comparing the results with those obtained earlier by MC simulation of energy probability distribution (strong first-order transition for model with $`q=3`$ and $`\sigma =0.1`$); second, by using the same graph approach we wish to analyze an example of a non-integer $`q`$ model with a second-order phase transition  (we choose $`q=0.5`$ and $`\sigma =0.1`$ and $`0.8`$).

The result of simulations for the case with $`q=3,\sigma =0.1`$ is shown in Fig. 2. It presents $`P_N`$ versus $`w_N`$ for system sizes ranging from $`200`$ to $`400`$, at temperatures $`T_N`$ where two peaks of the approximately same height appear. The figure shows that the depth of the minima increases with $`N`$. According to the discussion in the preceding subsection, a behaviour like that is characteristic for a first-order phase transition. This conclusion is confirmed by the results of the MC simulations for the energy probability distribution.

The temperatures $`T_N`$, where $`200N400`$ are presented in Table 2. Compared to the MF results and also to those for higher values of $`\sigma `$, the results for $`\sigma =0.1`$ converge very slowly, so that, for sizes considered here, we do not expect so good accuracy in extrapolation to $`N\mathrm{}`$. The extrapolation which takes one correction term, of the form $`T_N=T_{LR,t}+b/N^x`$ gives $`T_{LR,t}=7.73`$ with the convergence exponent $`x=0.1`$, while taking the form with two correction terms, e.g. $`T_N=T_{LR,t}+b/N^x+c/N^{2x}`$ gives by LSA fit $`x=0.12`$ and $`T_{LR,t}=7.45`$, value which differs by $`4\%`$. The first result for $`T_{LR,t}`$ shows the discrepancy of $`8\%`$ compared to the improved finite-range scaling (FRS)  result $`(T^{\mathrm{FRS}}=7.14)`$ cited in , and even larger discrepancies when compared to the renormalisation-group (RG) result $`(T^{\mathrm{RG}}=6.72)`$  and the earlier MC result $`(T^{\mathrm{MC}}=6.25)`$ . (Notice however that, due to the same slow convergence, the difference between these previous results is also quite significant when $`\sigma =0.1`$.) Apart from the reduced precision due to the small convergence exponent, the additionnal source of deviation can come from the crossover effect and be the consequence of the limitation to relatively small sizes for which the system does not exhibit full qualities of the first-order transition.

In order to investigate the case with non-integer $`q`$, we also considered the case $`q=0.5`$ with $`\sigma =0.1`$ and $`0.8`$, where by  and  the transitions are of classical and non-classical second-order type, respectively. The simulations were performed for a wide range of temperatures around critical temperatures obtained earlier by FRS . As one can expect on the grounds of the discussion at the beginning of the preceding subsection, no double-peak structure in graph-weight probability distribution has been found.

## 4 Conclusion

The graph-weight probability distribution is introduced and applied to the analysis of the order of transition on two special cases of the Potts model: its mean-field case and the case in $`1d`$ with power-law decaying interactions. The mean-field case was analyzed for $`q=3`$ and several non-integer values approaching $`q_c=2`$. The long-range case was examined for $`q=3`$ and $`q=0.5`$ with $`\sigma =0.1`$ and $`\sigma =0.1,0.8`$, respectively. The simulations were limited by time to the systems of sizes $`N400`$.

It is shown that the graph weight is an appropriate quantity to distinguish the coexisting phases at the first-order transition point. The physical interpretation of the graph weight becomes transparent through the equations (6) and (7) which relate the average number of active links and clusters to the average of energy and the free energy of the system.

By analyzing the graph-weight probability distributions of the MF and $`1d`$ LR Potts models, we have observed the behaviour of distributions characteristic for first-order transitions in the case of $`q=3`$ and $`2.8`$ MF models and $`q=3,\sigma =0.1`$ LR model, while the transition of the second order was obtained in the case of $`q=0.5,\sigma =0.1,0.8`$ LR model.

Transition temperatures have also been calculated. The estimated values of those temperatures in the thermodynamic limit agree for the MF case values within $`1\%`$, while for the case with power-law decaying interactions, where only the approximate results are available, it agrees with the discrepancy of $`5\%`$ to $`24\%`$ depending on the method.

All of the above facts qualify the graph-weight probability distribution as an alternative quantity in investigations of the order of transitions in Potts models, capable to deal directly with non-integer values of $`q`$. This makes the graph approach interesting in the framework of efforts of determination of the border between the first- and second-order regions by continuously varying $`q`$ in the $`(q,d)`$ plane of SR models or $`(q,\sigma )`$ plane in $`1d`$ LR models.

Figure captions

Fig. 1: The simulation of the graph-weight probability distribution $`P_N`$ of the $`q=3`$ MF Potts model versus $`w_N\mathrm{ln}W_N/N\mathrm{ln}N`$ performed at temperatures $`T_N`$ (see in text).

Fig. 2: The simulation of graph-weight probability distribution $`P_N`$ of the $`q=3,\sigma =0.1`$ $`1d`$ LR Potts model versus $`w_N\mathrm{ln}W_N/N\mathrm{ln}N`$ performed at temperatures $`T_N`$ (see in text).