Phase Transition Analysis of the Potts-SOS Model with Spin Set {-1,0, 1} on the Cayley Tree

1. Introduction

The phenomenon of phase transitions is prevalent both in nature and in everyday life, making it a topic of great interest [1,2,3]. Particularly for lattice models defined on cubic lattices, this subject has captivated many scientists, showing an ever-growing trend in research [4,5]. Significant attention has been given to phase transitions in fields such as statistical mechanics and condensed matter physics. Due to the inherent structure of Cayley trees (CTs) and Bethe lattices (BLs), determining the phase transition behavior in lattice models (such as Ising, Potts, and SOS) on these lattices is more straightforward [6]. This research can provide valuable insights and potentially inspire further studies on cubic lattices. The Potts model is a generalization of the Ising model [7] used in statistical mechanics to study phase transitions, magnetism, and critical phenomena [8]. It is particularly useful in understanding systems with more than two possible states at each site, making it applicable to a wide range of physical, chemical, and even biological systems [9].

The solid-on-solid (SOS) model, a statistical mechanics model, is used to study the behavior of interfaces in materials [10]. The SOS model on the CT was introduced by Rozikov and Suhov [11]. In [10], the scenario of arbitrary conditions is addressed, and a vector-valued functional equation is derived to describe potential boundary laws of the model. Recently, in [12], Rahmatullaev and Abraev have provided the phase transition problem of the SOS model on the CT and the corresponding classification of Gibbs measures.

As discussed in our previous studies [13,14,15], we employed three primary methods to analyze the thermodynamic properties and phase transition phenomena of lattice models defined on the CTs:

(1)

Kolmogorov consistency condition,

(2)

Cavity method (utilizing the self-similarity of the Cayley tree),

(3)

Transition matrix method.

Recently, structures known as mixed-lattice models have become the focus of extensive research. Consequently, the Potts-SOS model can be viewed as a type of mixed-lattice model. Rahmetullaev et al. [16] studied the phase transition problem of the Potts-SOS model with the spin set $\{0,1,2\}$ on the CT using the Kolmogorov consistency condition. In [13], we demonstrated that there is no phase transition in the 1-D Potts-SOS model with the spin set $\{-1,0,+1\}$ and calculated various thermodynamic quantities using the transition matrix method. In this study, we aim to explore the uniqueness and construction of Gibbs measures associated with a mixed-type lattice model, characterized by Potts and SOS Hamiltonians with two distinct coupling constants, using the self-similarity of the second-order semi-infinite CT. To achieve this, we will first define the partial partition functions and then derive the Gibbs measures using the cavity method. To determine the number of these Gibbs measures, we will identify the fixed points of the two-dimensional dynamical system that we construct.

As is well known, the phase transition problem of lattice models on the CT has been a subject of study for a long time [2,17]. For a phase transition to occur, the number of Gibbs measures for the given model must exceed one. Generally, determining the exact number of Gibbs measures directly is a challenging problem [15]. Instead, researchers investigate the fixed points of the dynamical system formed by the iterative equations, as there is a one-to-one correspondence between these fixed points and the Gibbs measures [6]. Consequently, the existence and uniqueness of fixed points are analyzed to determine the presence of a phase transition for given lattice model. However, identifying all fixed points analytically, especially for high-dimensional dynamical systems, can be difficult [18,19].

The presence of repelling fixed points suggests the existence of additional fixed points [20,21]. There are several methods to determine whether a fixed point is repelling. One approach is to check if the derivative of the function at the fixed point is greater than one [22]. Another method involves identifying regions where the largest eigenvalue of the Jacobian matrix associated with the n-dimensional dynamical system exceeds one [22,23].

In our recent investigation, we have analyzed the stability of fixed points to address the phase transition problem for lattice models on CTs [24,25]. In this paper, we will employ both methods—examining the derivative and analyzing the eigenvalues of the Jacobian matrix—to determine the regions where phase transitions occur as a function of the coupling constants.

In our previous research, we employed both the boundary field condition [26,27] and the self-similarity property [14,15,28,29] of the CT to derive the partial partition functions for the Potts-SOS model on a uniform one-order CT. Our findings indicated that the partial partition functions obtained through both methods are equivalent, with the self-similarity method proving to be a more straightforward approach compared to the boundary field condition. In the current paper, we will build on this investigation by examining the thermodynamic properties of the Potts-SOS model using the self-similarity approach of the CT. Previous work [13] established that no phase transition occurs for the Potts-SOS model on a one-order CT, and therefore, on the $\mathbb{N}$ lattice. In this study, our goal is to demonstrate the existence of a phase transition by analyzing the partial partition functions associated with the spin set $\{-1,0,+1\}$.

2. Preliminaries and Construction of the Partial Partition Functions

In this section, we introduce key concepts and results that are fundamental to our findings. To avoid redundancy, we will not go into detailed information about CT in this article. Instead, we suggest readers consult relevant literature for more information [6,26,30]. Here, we will briefly outline a few characteristics of the CT. In a CT of order k, each node is connected to $k+1$ neighbors, except for the nodes at the boundary (outermost level), which have only one neighbor. The tree is organized into levels, starting with a single root node at the top (level 0), with each subsequent level branching out with k branches per node. Unlike square or cubic lattices, the CT contains no closed loops or cycles.

Since the vertices of the CT are loopless [26,27,30], the distance between any two vertices $x,y\in V$ is defined by the number of edges along the path connecting them. If $x^{\left(0\right)}$ is the root of the CT, then $d(x^{\left(0\right)},x)$ represents the distance from the root to the vertex x, measured by the number of levels separating them. Therefore, if we denote the set of vertices from the root up to the n-th level (or shell) as $V_n$, this set can be briefly expressed as $V_n=\{x\in V:d(x^{\left(0\right)},x)\le n\}$.

We analyze the Potts-SOS model on the second-order CT, defined by the Hamiltonian:

$$\begin{array}{c}\hfill H\left(\sigma \right)=-J_p\sum _{\langle x,y\rangle }{\delta }_{\sigma \left(x\right)\sigma \left(y\right)}-J\sum _{\rangle x,y\langle }|\sigma \left(x\right)-\sigma \left(y\right)|-h_x\sum _{x\in V}\sigma \left(x\right),\end{array}$$

(1)

where $\langle x,y\rangle$ represents the energy associated with nearest-neighbor bonds and the third term is the energy of each site.

Let V denote the set of vertices of the semi-infinite CT, and U be a finite subset of V. We denote the restriction of $\sigma$ to U as $\sigma \left(U\right)$. Given a fixed boundary configuration $\overline{\sigma }(V\setminus U)$, the total energy of $\sigma \left(U\right)$ is given by:

$$\begin{array}{cc}\hfill H_U(\sigma \left(U\right)\mid \overline{\sigma }(V\setminus U))& =-J_p\left(\sum _{\langle x,y\rangle :x,y\in U}{\delta }_{\sigma \left(x\right)\sigma \left(y\right)}+\sum _{\langle x,y\rangle :x\in U,y\notin U}{\delta }_{\sigma \left(x\right)\sigma \left(y\right)}\right)\hfill \\ \hfill & -J\left(\sum _{\langle x,y\rangle :x,y\in U}|\sigma \left(x\right)-\sigma \left(y\right)|+\sum _{\langle x,y\rangle :x\in U,y\notin U}|\sigma \left(x\right)-\sigma \left(y\right)|\right)-h_x\sum _{x\in U}\sigma \left(x\right).\hfill \end{array}$$

(2)

Given $\overline{\sigma }(V\setminus U)$, the partition function $Z_U\left(\overline{\sigma }(V\setminus U)\right)$ in volume U is defined as:

$$Z_U\left(\overline{\sigma }(V\setminus U)\right)=\sum _{\sigma \left(U\right)\in {\mathsf{\Phi }}^U}exp\left(-\beta H_U(\sigma \left(U\right)\mid \overline{\sigma }(V\setminus U))\right),$$

where $\beta ={\displaystyle \frac{1}{T}}$ is the inverse temperature.

For simplicity, $Z_{V_n}$ and the configuration $\sigma \left(V_n\right)$ in volume $V_n$ are referred to as ${\sigma }_n$ and $Z^{\left(n\right)}$, respectively. The total partition function $Z^{\left(n\right)}$ can be expressed as:

$$Z^{\left(n\right)}=\sum _{i\in \mathsf{\Phi }}{Z}_i^{\left(n\right)},$$

where:

$$Z_i^{\left(n\right)}=\sum _{{\sigma }_n\in {\mathsf{\Phi }}^{V_n}:\sigma \left(x^{\left(0\right)}\right)=i}exp\left(-\beta H_{V_n}({\sigma }_n\mid \overline{\sigma }(V\setminus V_n))\right).$$

2.1. Self-Similarity Approach to Determine the Partial Partition Functions

Consider a CT of order two with the root vertex $x^{\left(0\right)}$ (refer to Figure 1a). From the root vertex $x^{\left(0\right)}$, two edges extend to vertices $x_1$ and $x_2$. When viewing the lattice G as infinite, we obtain two infinite subgraphs $G_1$ and $G_2$ that are equivalent. Therefore, the tree exhibits self-similarity such that $G=G_1=G_2$ (see [30] for more details).

In this section, we will set the external field variable h to zero. For a finite configuration ${\sigma }_n\in {\{-1,0,+1\}}^{V_n}$, the Hamiltonian in Equation (1) can be rewritten as follows:

$$\begin{array}{c}\hfill H\left({\sigma }_n\right)=-J_p\sum _{\langle x,y\rangle }{\delta }_{\sigma \left(x\right),\sigma \left(y\right)}-J\sum _{\langle x,y\rangle }|\sigma \left(x\right)-\sigma \left(y\right)|,\end{array}$$

(3)

where ${\delta }_{\sigma \left(x\right),\sigma \left(y\right)}$ represents the Kronecker delta function.

For convenience, we use the following equation to compute partial partition functions, as illustrated in Figure 1(b):

$$\begin{array}{cc}\hfill Z_{i_0}^{\left(n\right)}& =\sum _{i_1=-1}^1\sum _{i_2=-1}^{1}exp\left(\beta J(|i_0-i_1|+|i_0-i_2|)+\beta J_p({\delta }_{i_0{i}_1}+{\delta }_{i_0{i}_2})\right)Z_{i_1}^{(n-1)}{Z}_{i_1}^{(n-1)}.\hfill \end{array}$$

(4)

After performing elementary algebraic operations, for $i_0=-1,0,1$, the partial partition functions are obtained from Equation (4) using the following formulas:

$$\begin{array}{cc}\hfill Z_{-1}^{\left(n\right)}& ={\left(e^{\beta J_p}{Z}_{-1}^{(n-1)}+e^{J\beta }{Z}_0^{(n-1)}+e^{2J\beta }{Z}_1^{(n-1)}\right)}^2,\hfill \\ \hfill Z_0^{\left(n\right)}& ={\left(e^{J\beta }{Z}_{-1}^{(n-1)}+e^{\beta J_p}{Z}_0^{(n-1)}+e^{J\beta }{Z}_1^{(n-1)}\right)}^2,\hfill \\ \hfill Z_1^{\left(n\right)}& ={\left(e^{2J\beta }{Z}_{-1}^{(n-1)}+e^{J\beta }{Z}_0^{(n-1)}+e^{\beta J_p}{Z}_1^{(n-1)}\right)}^2.\hfill \end{array}$$

(5)

Let $P\left({\sigma }^{\left(0\right)}\right)$ represent the probability that the spin at the root vertex $x^{\left(0\right)}$ of the semi-infinite CT has the value ${\sigma }^{\left(0\right)}$. For ${\sigma }^{\left(0\right)}\in \{-1,0,+1\}$, we obtain:

$$P\left({\sigma }^{\left(0\right)}\right)={\displaystyle \frac{Z_{{\sigma }^{\left(0\right)}}^{\left(n\right)}}{Z_{-1}^{\left(n\right)}+Z_0^{\left(n\right)}+Z_1^{\left(n\right)}}}.$$

(6)

Note that $\left(P(-1),P\left(0\right),P\left(1\right)\right)$ is a probability vector. After substitution, we have:

$$\left\{\begin{array}{c}{U}^{\left(n\right)}={\displaystyle \frac{P(-1)}{P\left(0\right)}}={\displaystyle \frac{Z_{-1}^{\left(n\right)}}{Z_0^{\left(n\right)}}},\hfill \\ V^{\left(n\right)}={\displaystyle \frac{P(+1)}{P\left(0\right)}}={\displaystyle \frac{Z_1^{\left(n\right)}}{Z_0^{\left(n\right)}}},\hfill \end{array}\right.$$

From the definition, we get the following system of recurrent equations:

$$\begin{array}{ccc}\hfill U^{\left(n\right)}& =& {\displaystyle \frac{{\left(e^{2J\beta }{V}^{(n-1)}+e^{J\beta }+e^{\beta J_p}{U}^{(n-1)}\right)}^2}{{\left(e^{\beta J_p}+e^{J\beta }(U^{(n-1)}+V^{(n-1)})\right)}^2}},\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill V^{\left(n\right)}& =& {\displaystyle \frac{{\left(e^{\beta J_p}{V}^{(n-1)}+e^{J\beta }+e^{2J\beta }{U}^{(n-1)}\right)}^2}{{\left(e^{\beta J_p}+e^{J\beta }(U^{(n-1)}+V^{(n-1)})\right)}^2}}.\hfill \end{array}$$

(8)

Thus, we have obtained a two-dimensional recursive dynamical system with

$$F:{\mathbb{R}}_{+}^2\to {\mathbb{R}}_{+}^2.$$

(9)

If we assume $u=\underset{n\to \infty }{lim}{U}^{\left(n\right)}$ and $v=\underset{n\to \infty }{lim}{V}^{\left(n\right)}$, then we get the following two-dimensional dynamical system;

$$\begin{array}{ccc}\hfill u& =& {\displaystyle \frac{{\left({\theta }^{2}v+\theta +\rho u\right)}^2}{{(\rho +\theta (v+u))}^2}},\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill v& =& {\displaystyle \frac{{\left(\rho v+\theta +{\theta }^{2}u\right)}^2}{{(\rho +\theta (v+u))}^2}},\hfill \end{array}$$

(11)

where $\theta =e^{\beta J}$ and $\rho =e^{\beta J_p}$.

Let us change the variables $u=t^2$ and $v=s^2$. From the equations (10) and (11), we have

$$\begin{array}{c}\hfill t^2={\displaystyle \frac{{\left({\theta }^2{s}^2+\theta +\rho t^2\right)}^2}{{(\rho +\theta (s^2+t^2))}^2}},\end{array}$$

(12)

$$\begin{array}{c}\hfill s^2={\displaystyle \frac{{\left(\rho s^2+\theta +{\theta }^2{t}^2\right)}^2}{{(\rho +\theta (s^2+t^2))}^2}}.\end{array}$$

(13)

1. For the first equation: Since $t>0$ and $s>0$, taking the square root of both sides of the equation (12), we have

$$t={\displaystyle \frac{{\theta }^2{s}^2+\theta +\rho t^2}{\rho +\theta (s^2+t^2)}}.$$

After some algebraic manipulation, we arrange to group like terms:

$$t\rho +t\theta s^2+t^3\theta ={\theta }^2{s}^2+\theta +\rho t^2.$$

2. Similarly, for the second equation: Since $t>0$ and $s>0$, taking the square root of both sides of the equation (13), we have

$$\begin{array}{c}\hfill s={\displaystyle \frac{\rho s^2+\theta +{\theta }^2{t}^2}{\rho +\theta (s^2+t^2)}}\end{array}$$

(14)

Similarly, after some algebraic manipulation, we arrange to group like terms:

$$s\rho +s^3\theta +s\theta t^2=\rho s^2+\theta +{\theta }^2{t}^2.$$

From the two simplified equations, we have:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}t\rho +t\theta s^2+t^3\theta ={\theta }^2{s}^2+\theta +\rho t^2,\\ s\rho +s^3\theta +s\theta t^2=\rho s^2+\theta +{\theta }^2{t}^2.\end{array}\right.\end{array}$$

(15)

Since $\rho \ne 0$ and $\theta \ne 0$, we need to find t and s such that both equations are satisfied.

We can set $s=t$ or $s+t=1$ as before to check for consistency:

Case 1: $s=t$

Substitute $s=t$ into the simplified equations:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}t\rho +t\theta t^2+t^3\theta ={\theta }^2{t}^2+\theta +\rho t^2,\\ t\rho +t^3\theta +t\theta t^2=\rho t^2+\theta +{\theta }^2{t}^2.\end{array}\right.\end{array}$$

(16)

Both equations simplify to:

$$t\rho +t^3\theta +t^3\theta ={\theta }^2{t}^2+\theta +\rho t^2.$$

Since the equations are symmetric, this case is consistent.

Case 2: $s+t=1$

Substitute $s=1-t$ into the simplified equations (15) and check for consistency:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}t\rho +t\theta {(1-t)}^2+t^3\theta ={\theta }^2{(1-t)}^2+\theta +\rho t^2,\\ (1-t)\rho +{(1-t)}^3\theta +(1-t)\theta t^2=\rho {(1-t)}^2+\theta +{\theta }^2{t}^2.\end{array}\right.\end{array}$$

(17)

Let us start by simplifying and solving the system of equations (17) with the assumption $s+t=1$. Substitute $s=1-t$ in both equations:

Equation 1:

$$t\rho +t\theta {(1-t)}^2+t^3\theta ={\theta }^2{(1-t)}^2+\theta +\rho t^2.$$

Expand and simplify each term:

$$\left\{\begin{array}{c}t\rho +t\theta (1-2t+t^2)+t^3\theta ={\theta }^2(1-2t+t^2)+\theta +\rho t^2,\\ t\rho +t\theta -2t^2\theta +t^3\theta +t^3\theta ={\theta }^2-2t{\theta }^2+t^2{\theta }^2+\theta +\rho t^2.\end{array}\right.$$

Combine like terms:

$$t\rho +t\theta -2t^2\theta +2t^3\theta ={\theta }^2-2t{\theta }^2+t^2{\theta }^2+\theta +\rho t^2.$$

Equation 2:

$$(1-t)\rho +{(1-t)}^3\theta +(1-t)\theta t^2=\rho {(1-t)}^2+\theta +{\theta }^2{t}^2.$$

Expand and simplify each term:

$$\left\{\begin{array}{c}(1-t)\rho +(1-3t+3t^2-t^3)\theta +(1-t)\theta t^2=\rho (1-2t+t^2)+\theta +{\theta }^2{t}^2,\\ \rho -t\rho +\theta -3t\theta +3t^2\theta -t^3\theta +\theta t^2-t\theta t^2=\rho -2t\rho +t^2\rho +\theta +{\theta }^2{t}^2.\end{array}\right.$$

Combine like terms:

$$\rho -t\rho +\theta -3t\theta +3t^2\theta -t^3\theta +\theta t^2-t^3\theta =\rho -2t\rho +t^2\rho +\theta +{\theta }^2{t}^2.$$

So, we have two simplified equations to solve for t:

$$\left\{\begin{array}{c}t\rho +t\theta -2t^2\theta +2t^3\theta ={\theta }^2-2t{\theta }^2+t^2{\theta }^2+\theta +\rho t^2,\\ \rho -t\rho +\theta -3t\theta +3t^2\theta -t^3\theta +\theta t^2-t^3\theta =\rho -2t\rho +t^2\rho +\theta +{\theta }^2{t}^2.\end{array}\right.$$

Let us isolate terms involving t:

• Rearrange Equation (1):

$$t\rho +t\theta -2t^2\theta +2t^3\theta ={\theta }^2-2t{\theta }^2+t^2{\theta }^2+\theta +\rho t^2,$$

$$t(\rho +\theta )-2t^2\theta +2t^3\theta ={\theta }^2-2t{\theta }^2+t^2{\theta }^2+\theta +\rho t^2.$$

2.

Rearrange Equation (2):

$$\rho -t\rho +\theta -3t\theta +3t^2\theta -t^3\theta +\theta t^2-t^3\theta =\rho -2t\rho +t^2\rho +\theta +{\theta }^2{t}^2,$$

$$\rho (1-t)+\theta (1-3t+3t^2+t^2-t^3-t^3)=\rho (1-2t+t^2)+\theta +{\theta }^2{t}^2.$$

We can solve these equations numerically to find the values of t. This substitution will lead to a more complex equation, which may have solutions depending on specific values of $\rho$ and $\theta$. In conclusion, the solutions to the system of equations can be either $s=t$ or $s+t=1$. These solutions ensure the consistency of the given equations after substituting $u=t^2$ and $v=s^2$.

Consequently, from the system of equations (12) and (13), we derive the following dynamical system:

$$\begin{array}{cc}\hfill f_1(t,s)& ={\displaystyle \frac{\left({\theta }^2{s}^2+\theta +\rho t^2\right)}{(\rho +\theta (s^2+t^2))}},\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill f_2(t,s)& ={\displaystyle \frac{\left(\rho s^2+\theta +{\theta }^2{t}^2\right)}{(\rho +\theta (s^2+t^2))}}.\hfill \end{array}$$

(19)

In Section 4, we will discuss the stability analysis of the fixed points of the system of Equations (18) and (19).

3. Phase Transition for the Potts-SOS Model

In this section we will find the fixed points satisfying the Equations (10) and (). Let us investigate the existence and multiplicity of fixed points of the given equation system by examining Case 1. Assuming $s=t$, from the Equation (14), we obtain:

$$\begin{array}{c}\hfill F\left(s\right):={\displaystyle \frac{(\rho +{\theta }^2)s^2+\theta }{\rho +2\theta s^2}}=s.\end{array}$$

(20)

Now, let us find the roots that satisfy this equation. First, multiply both sides by $\rho +2\theta s^2$ to clear the denominator:

$$s(\rho +2\theta s^2)=(\rho +{\theta }^2)s^2+\theta .$$

Next, we obtain a cubic equation in s as

$$\begin{array}{c}\hfill 2\theta s^3-(\rho +{\theta }^2)s^2+\rho s-\theta =0.\end{array}$$

(21)

To find the roots of Equation (21), we can apply various methods such as factoring (if possible), using the Rational Root Theorem, or numerical methods. The given equation can be solved numerically with Mathematica [31] or other computing tools.

For a cubic polynomial $a{s}^3+b{s}^2+cs+d=0$, the roots can be found using the cubic formula, but it is often simpler to use numerical methods or software tools for specific values of $\rho$ and $\theta$. From the Cardano’s formula, the discriminant of a cubic equation $a{x}^3+b{x}^2+cx+d=0$ is given by:

$$\mathsf{\Delta }=18abcd-4b^{3}d+b^2{c}^2-4a{c}^3-27a^2{d}^2.$$

Substituting $a=2\theta$, $b=-(\rho +{\theta }^2)$, $c=\rho$, and $d=-\theta$ into the discriminant formula, we can compute $\mathsf{\Delta }$ as

$$\mathsf{\Delta }=\theta ({\rho }^2-8\theta \rho +38{\theta }^2-12{\theta }^3+{\theta }^4)(2\rho -6\theta )-4{\theta }^7-108{\theta }^4.$$

Depending on the value of $\mathsf{\Delta }$:

• If $\mathsf{\Delta }>0$, the cubic equation has three distinct real roots.
• If $\mathsf{\Delta }=0$, the cubic equation has a multiple root and all its roots are real.
• If $\mathsf{\Delta }<0$, the cubic equation has one real root and two nonreal complex conjugate roots.

Analyzing the discriminant will help us determine the nature and multiplicity of the fixed points for the given equation. Alternatively, solving the cubic equation numerically for specific $\rho$ and $\theta$ will provide explicit values for the roots s.

Figure 2 illustrates the blue region where $\mathsf{\Delta }\ge 0$. This indicates that the cubic polynomial $2\theta s^3-(\rho +{\theta }^2)s^2+\rho s-\theta$ has real or multiple roots within this area. Consequently, this region signifies the occurrence of a phase transition. By choosing values of $\theta$ and $\rho$ from this region, we can numerically demonstrate that the function $F\left(s\right)$ given in (20) has three distinct fixed points.

Figure 3 shows various graphs of the $F\left(s\right)$ function to verify the presence of a phase transition. For $(\theta =0.49,\rho =22)$ and $(\theta =3.49,\rho =22)$, there are three distinct fixed points. In contrast, for $(\theta =3.49,\rho =2)$, there is only one fixed point. We know there is a one-to-one correspondence between Gibbs measures and the solutions of Equations (10) and (11) (see [6] for details). Thus, a phase transition is observed in the first two cases but not in the third. Therefore, we have examined the non-uniqueness of the solutions to these equations to determine if there is a phase transition in the Potts-SOS model on the second-order CT.

4. Stability Analysis of the Two-Dimensional Dynamical system

In this section, we will perform the stability analysis of the fixed points of the two-variable equation system (18) and (19). It is known that the presence of other fixed points in the region where a fixed point is repelling (unstable) indicates the existence of multiple Gibbs measures for the corresponding Potts-SOS model, signifying a phase transition.

In stability analysis, the eigenvectors of the Jacobian matrix are crucial for assessing the stability of a fixed point [32]. To compute the Jacobian matrix of the dynamical system given in (18) and (19) at the point $(t,s)=(2,2)$, we first need to find the partial derivatives of $f_1(t,s)$ and $f_2(t,s)$ with respect to t and s. So, we compute the Jacobian matrix J of the dynamical system given in (18) and () as follows:

$$J=\left[\begin{array}{cc}{\displaystyle \frac{\partial f_1}{\partial t}}& {\displaystyle \frac{\partial f_1}{\partial s}}\\ {\displaystyle \frac{\partial f_2}{\partial t}}& {\displaystyle \frac{\partial f_2}{\partial s}}\end{array}\right].$$

Let us calculate each partial derivative. For $f_1(t,s)$, we have

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial f_1}{\partial t}}& =& {\displaystyle \frac{2t\left(-{\theta }^2-s^2{\theta }^3+s^2\theta \rho +{\rho }^2\right)}{{\left(s^2\theta +t^2\theta +\rho \right)}^2}},\hfill \\ \hfill {\displaystyle \frac{\partial f_1}{\partial s}}& =& {\displaystyle \frac{2s\theta \left(-\theta +t^2{\theta }^2-t^2\rho +\theta \rho \right)}{{\left(s^2\theta +t^2\theta +\rho \right)}^2}}.\hfill \end{array}$$

Similarly, for $f_2(t,s)$, we have

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial f_2}{\partial t}}& =& {\displaystyle \frac{2t\theta \left(-\theta +s^2{\theta }^2-s^2\rho +\theta \rho \right)}{{\left(s^2\theta +t^2\theta +\rho \right)}^2}},\hfill \\ \hfill {\displaystyle \frac{\partial f_2}{\partial s}}& =& {\displaystyle \frac{2s\left(-{\theta }^2-t^2{\theta }^3+t^2\theta \rho +{\rho }^2\right)}{{\left(s^2\theta +t^2\theta +\rho \right)}^2}}.\hfill \end{array}$$

So, we obtain the Jacobian matrix as:

$$\begin{array}{cc}\hfill J(t,s):& =\left(\begin{array}{cc}{\displaystyle \frac{\partial f}{\partial t}}& {\displaystyle \frac{\partial f}{\partial s}}\\ {\displaystyle \frac{\partial g}{\partial t}}& {\displaystyle \frac{\partial g}{\partial s}}\end{array}\right)=\left(\begin{array}{cc}{\displaystyle \frac{2t\left(-{\theta }^2-s^2{\theta }^3+s^2\theta \rho +{\rho }^2\right)}{{\left(s^2\theta +t^2\theta +\rho \right)}^2}}& {\displaystyle \frac{2s\theta \left(-\theta +t^2{\theta }^2-t^2\rho +\theta \rho \right)}{{\left(s^2\theta +t^2\theta +\rho \right)}^2}}\\ {\displaystyle \frac{2t\theta \left(-\theta +s^2{\theta }^2-s^2\rho +\theta \rho \right)}{{\left(s^2\theta +t^2\theta +\rho \right)}^2}}& {\displaystyle \frac{2s\left(-{\theta }^2-t^2{\theta }^3+t^2\theta \rho +{\rho }^2\right)}{{\left(s^2\theta +t^2\theta +\rho \right)}^2}}\end{array}\right).\hfill \end{array}$$

(22)

Using vector calculus of a vector-valued function of several variables, we obtain the Jacobian matrix J given in (22) associated with the dynamical system given in (18) and (19) at the fixed point $(t,s)=(2,2)$ as

$$\begin{array}{cc}\hfill J(2,2)& =\left(\begin{array}{cc}{\displaystyle \frac{4\left(-{\theta }^2-4{\theta }^3+4\theta \rho +{\rho }^2\right)}{{(8\theta +\rho )}^2}}& {\displaystyle \frac{4\theta \left(4{\theta }^2+\theta (-1+\rho )-4\rho \right)}{{(8\theta +\rho )}^2}}\\ {\displaystyle \frac{4\theta \left(4{\theta }^2+\theta (-1+\rho )-4\rho \right)}{{(8\theta +\rho )}^2}}& {\displaystyle \frac{4\left(-{\theta }^2-4{\theta }^3+4\theta \rho +{\rho }^2\right)}{{(8\theta +\rho )}^2}}\end{array}\right).\hfill \end{array}$$

(23)

After elementary linear algebra operations, we compute the eigenvalues of the Jacobian matrix (23) as

$$\begin{array}{ccc}\hfill {\lambda }_1& =& {\displaystyle \frac{4\rho -4{\theta }^2}{\rho +8\theta }},\hfill \\ \hfill {\lambda }_2& =& {\displaystyle \frac{4{\rho }^2+4\rho {\theta }^2-8{\theta }^2}{{\rho }^2+16\rho \theta +64{\theta }^2}}.\hfill \end{array}$$

These eigenvalues help determine the stability and behavior of the system near the point $(t,s)=(2,2)$. As is known, in the region where the condition $|{\lambda }_{max}|-1>0$ is satisfied, the fixed point $(s,t)=(2,2)$ becomes repelling (unstable). Therefore, the presence of additional fixed points in this region indicates that a phase transition occurs. Let us determine the region to which the pair $(\theta ,\rho )$ that satisfies the inequality $|{\lambda }_2|-1>0$ belongs.

We consider the inequality

$$\begin{array}{c}\hfill 4\left|{\displaystyle \frac{-2{\theta }^2+{\theta }^2\rho +{\rho }^2}{{(8\theta +\rho )}^2}}\right|>1.\end{array}$$

(24)

Therefore, we analyze the following inequality:

$$\left|{\displaystyle \frac{-8{\theta }^2+4{\theta }^2\rho +4{\rho }^2}{{(8\theta +\rho )}^2}}\right|>1.$$

This inequality can be broken down into two separate inequalities:

Case I: ${\displaystyle \frac{-8{\theta }^2+4{\theta }^2\rho +4{\rho }^2}{{(8\theta +\rho )}^2}}>1.$

Multiplying both sides by ${(8\theta +\rho )}^2$, which is always positive, yields $-8{\theta }^2+4{\theta }^2\rho +4{\rho }^2>{(8\theta +\rho )}^2.$ Expanding and rearranging the terms, we obtain:

$$-8{\theta }^2+4{\theta }^2\rho +4{\rho }^2>64{\theta }^2+16\theta \rho +{\rho }^2,$$

which simplifies to $-72{\theta }^2+4{\theta }^2\rho -16\theta \rho +3{\rho }^2>0.$

Case II: ${\displaystyle \frac{-8{\theta }^2+4{\theta }^2\rho +4{\rho }^2}{{(8\theta +\rho )}^2}}<-1.$

Again, multiplying both sides by ${(8\theta +\rho )}^2$ results in $-8{\theta }^2+4{\theta }^2\rho +4{\rho }^2<-{(8\theta +\rho )}^2.$ Expanding and simplifying gives:

$$-8{\theta }^2+4{\theta }^2\rho +4{\rho }^2<-64{\theta }^2-16\theta \rho -{\rho }^2,$$

which simplifies to $56{\theta }^2+4{\theta }^2\rho +16\theta \rho +5{\rho }^2<0.$

Therefore, the regions in which the absolute value of the given expression is greater than 1 are characterized by the following inequalities:

$$\begin{array}{cc}\hfill & 3{\rho }^2+4(\theta -4)\theta \rho -72{\theta }^2>0,\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill & 56{\theta }^2+4{\theta }^2\rho +16\theta \rho +5{\rho }^2<0.\hfill \end{array}$$

(26)

These inequalities delineate the areas of interest in the $\theta$-$\rho$ plane where the absolute value of the expression exceeds 1. Since inequality (25) is a quadratic equation in $\rho$, we will solve the corresponding quadratic equation $3{\rho }^2+4(\theta -4)\theta \rho -72{\theta }^2=0$ to find the critical values of $\rho$. This will help us understand the regions where the inequality holds. So, we have

$$\begin{array}{ccc}\hfill {\rho }_1& =& {\displaystyle \frac{2}{3}}\left(4\theta -{\theta }^2-\theta \sqrt{70-8\theta +{\theta }^2}\right),\hfill \\ \hfill {\rho }_2& =& {\displaystyle \frac{2}{3}}\left(4\theta -{\theta }^2+\theta \sqrt{70-8\theta +{\theta }^2}\right).\hfill \end{array}$$

Regarding the solution set of the inequality in (), given that both $\theta$ and $\rho$ are positive, the solution set for this inequality is empty.

Now, let us proceed with the analysis of the inequality in (25) once more. We consider the inequality:

$$\rho >{\displaystyle \frac{2}{3}}\left(4\theta -{\theta }^2+\theta \sqrt{70-8\theta +{\theta }^2}\right).$$

To analyze this inequality, we perform the following steps:

Let us consider the function by

$$\begin{array}{cc}\hfill f\left(\theta \right)& ={\displaystyle \frac{2}{3}}\left(4\theta -{\theta }^2+\theta \sqrt{70-8\theta +{\theta }^2}\right).\hfill \end{array}$$

(27)

The term inside the square root is $70-8\theta +{\theta }^2$. To ensure that the square root is real and non-negative, we require $70-8\theta +{\theta }^2\ge 0.$ Rewriting this quadratic expression, we get:

$${\theta }^2-8\theta +70\ge 0.$$

Solving the quadratic equation ${\theta }^2-8\theta +70=0$, yields complex roots, indicating that the quadratic expression ${\theta }^2-8\theta +70$ is always positive for all real $\theta$. Hence, the square root term $\sqrt{70-8\theta +{\theta }^2}$ is always real and non-negative.

As $\theta \to 0$: For small values of $\theta$, the expression simplifies approximately to:

$${\displaystyle \frac{2}{3}}\left(4\theta -{\theta }^2+\theta \sqrt{70-8\theta +{\theta }^2}\right)\approx {\displaystyle \frac{2}{3}}·4\theta ={\displaystyle \frac{8\theta }{3}}.$$

As $\theta \to \infty$: For large values of $\theta$, the dominant term in the expression is:

$${\displaystyle \frac{2}{3}}\left(4\theta -{\theta }^2+\theta \sqrt{70-8\theta +{\theta }^2}\right)\approx {\displaystyle \frac{2}{3}}\left(\theta ·\theta \right)={\displaystyle \frac{2{\theta }^2}{3}}.$$

To gain a deeper understanding of the function $f\left(\theta \right)$, we plot:

$$f\left(\theta \right)={\displaystyle \frac{2}{3}}\left(4\theta -{\theta }^2+\theta \sqrt{70-8\theta +{\theta }^2}\right).$$

This plot illustrates how the function $f\left(\theta \right)$ varies with $\theta$ and helps determine the feasible range for $\rho$ (see Figure 4(Left)). The blue area in Figure 4(Right) illustrates the region where the inequality ${\displaystyle \frac{-8{\theta }^2+4{\theta }^2\rho +4{\rho }^2}{{(8\theta +\rho )}^2}}>1$ is satisfied. This blue region corresponds to the domain in which the model experiences a phase transition.

5. Conclusions

In this present study, we build upon the author’s previous research on the thermodynamic properties of the one-dimensional Potts-SOS model on the lattice of positive natural numbers $\mathbb{N}$ in [Akin H, Phys. Scr. 99 (2024), 055231] by exploring a mixed-type model on a second-order semi-infinite CT. We have determined the Gibbs measures through the construction of partition functions and utilized the self-similarity of the CT to analyze the uniqueness of these measures. Our stability analysis of the dynamical system at specific fixed points enabled us to identify the regions where phase transitions occur. Our results show that, while the one-dimensional Potts-SOS model does not exhibit a phase transition, the same model does display a phase transition when applied to the second-order CT.