Calculation of Thermodynamic Quantities of 1D Ising Model with Mixed Spin-(s,(2t − 1)/2) by Means of Transfer Matrix

Abstract

In this paper, we consider the one-dimensional Ising model (shortly, 1D-MSIM) having mixed spin-$(s,(2t-1)/2)$ with the nearest neighbors and the external magnetic field. We establish the partition function of the model using the transfer matrix. We compute certain thermodynamic quantities for the 1D-MSIM. We find some precise formulas to determine the model’s free energy, entropy, magnetization, and susceptibility. By examining the iterative equations associated with the model, we use the cavity approach to investigate the phase transition problem. We numerically determine the model’s periodicity.

1. Introduction

Lenz introduced the one-dimensional (1D) classical spin-1/2 Ising model [1]. Compared to the spin-1/2 model, the spin-1 Ising model, also known as the Blume–Capel model, is more appropriate; hence, it was employed in the investigation of phase transitions in systems of three states [2,3]. Then, the spin-3/2 Ising model was used to extend the spin-1 Ising model’s conclusions [4,5].

The magnetic characteristics of mixed-spin systems have recently attracted both theoretical and experimental attention [6,7,8]. These systems are especially well suited for studying the magnetic properties of a particular class of magnetic materials, in addition to some technical applications [9]. Therefore, mixed-spin Ising models have been extensively studied in the literature lately [10]. In Ref. [11], we classified the disordered phases corresponding to the Ising model with spin-1 and spin 1/2 on the semi-finite Cayley tree (CT). Then, on the second-order CT, we examined the phase transition problem of the Ising model with spin-2 and spin 1/2. Furthermore, we discussed the chaotic behavior of the corresponding dynamical system [12]. The main methods for examining the properties of Ising models with the mixed spin on the square lattice [10], Bethe lattice [13], and CT [11,12] are Monte Carlo simulations [10], the cavity method [14,15,16], the Kolmogorov consistency method [11,12], the iteration method, and transfer matrix method. In [17], the magneto thermal parameters that characterize the magnetocaloric effect (MCE) behaviors, such as entropy, entropy change, and adiabatic cooling rate, are precisely calculated using the transfer matrix approach.

In the literature, there are numerous methods for determining the free energy of a given model [13,18]. The free energy of the lattice models on the Bethe lattice is calculated while boundary conditions are taken into consideration [13,18,19]. In-depth research has been performed on the thermodynamic properties of the 1D Ising model with single spin-s using the transfer matrix [20,21,22].

The author is not aware of any research on the thermodynamic properties of 1D Ising models with mixed spin. Due to its widespread applications, the matrix transfer approach is one of the most commonly used methods in statistical mechanics. In [21], the transfer matrix approach is considered to conduct an analytical study of the 1D Ising model with spin-s. To grasp physical issues as well as statistical mechanics, it is essential to compute thermodynamic quantities like free energy, entropy, magnetization, and magnetic susceptibility [22].

In this paper, to determine the thermodynamic quantities related to 1D-MSIM with the mixed spin-$(s,(2t-1)/2)$, we consider the transfer matrix approach. We enlarge a specific bond and establish the transfer matrix that corresponds to the bond from site j to site $j+1$. By considering the trace of several transfer matrices multiplied simultaneously, the partition function is reconstructed. Inspired by the results given for the 1D Ising model with the single spin-s [21], the partition function and the free energy of 1D-MSIM with mixed spin-$(s,\frac{(2t-1)}{2})$ having the nearest neighbor interactions and external field are calculated. To calculate the model’s free energy, entropy, magnetization, and susceptibility, some exact formulas are established via the corresponding transfer matrix. By using the cavity approach to the corresponding iterative equations of the model, the phase transition issue is studied. The periodicity of the model is also estimated, numerically.

2. Preliminaries

Here, we review definitions and key findings in the construction of the partition function for the 1D-MSIM with mixed spin-$(s,\frac{(2t-1)}{2})$ on the one-dimensional lattice $\mathbb{Z}$. We denote the set of integers by $\mathbb{Z}$ and the set of strictly positive integers by ${\mathbb{N}}^{+}$. Two vertices x and y, $x, y\in {\mathbb{N}}^{+}$ are called nearest-neighbor (NN) if there exists an edge $\ell \in L$ connecting them, which is denoted by $\ell =\langle x,y\rangle$. For $x,y\in {\mathbb{N}}^{+}$, the distance $d(x,y)$ is the length (the number of edges) of the shortest path connecting x with y. For any $x\in {\mathbb{N}}^{+}$, the direct successor of the vertex x is defined by $S\left(x\right)=\{y\in {\mathbb{N}}^{+}:d(x,y)=1\}.$

Denote the set of even natural numbers by $\mathbb{E}=\{2n:n\in \mathbb{N}\}$ and the set of odd natural numbers by $\mathbb{O}=\{2n-1:n\in {\mathbb{N}}^{+}\}$. We consider the mixed spin-state spaces denoted by $\Phi =\{\pm s:s\in {\mathbb{Z}}^{+}\}\cup \left\{0\right\}$ and $\Psi =\{\pm \frac{(2t-1)}{2}:t\in {\mathbb{Z}}^{+}\}$. For $A\subset {\mathbb{N}}^{+}$, a spin configuration ${\xi }_A$ on A is defined as a function.

$$x\in A\to {\xi }_A\left(x\right)=\left\{\begin{array}{cc}{\xi }_A\left(x\right)=\sigma \left(x\right)\in \Phi ,& x\in \mathbb{E}\\ {\xi }_A\left(x\right)=s\left(x\right)\in \Psi ,& x\in \mathbb{O}.\end{array}\right.$$

Let us consider finite set ${\mathbb{N}}_{2N+1}=\{1,2,\cdots ,2N+1\}$. Let ${\Omega }_E={\Phi }^{\mathbb{E}}$ and ${\Omega }_O={\Psi }^{\mathbb{O}}$ be the spaces of infinite configurations and ${\Omega }_{+,N}={\Phi }^{\mathbb{E}\cap {\mathbb{N}}_{2N+1}}$ and ${\Omega }_{-,N}={\Psi }^{\mathbb{O}\cap {\mathbb{N}}_{2N+1}}$ be the spaces of finite configurations. $\Xi ={\Omega }_E\times {\Omega }_O$ represents the configuration space. An element of ${\Omega }_E$ is denoted by $\sigma \left(x\right)$, for $x\in \mathbb{E}$. Similarly, an element of ${\Omega }_O$ is denoted by $s\left(x\right)$, for $x\in \mathbb{O}.$

In this paper, the spins in $\Phi$ and $\Psi$ can be placed on the odd-numbered vertices and the even-numbered vertices, respectively, when mixed spins are designed on the lattice ${\mathbb{N}}^{+}$ (see Figure 1).

Construction of the Partition Function Associated with the Model

Let us construct the partition function corresponding to 1D-MSIM with mixed spin $(s,\frac{2t-1}{2})$ on the lattice ${\mathbb{N}}^{+}$.

Let $\xi \in \Xi$, so that one has

$$\xi \left(x\right)=\left\{\begin{array}{cc}\sigma \left(x\right);\hfill & x\in \mathbb{E}\hfill \\ s\left(x\right);\hfill & x\in \mathbb{O},\hfill \end{array}\right.$$

where $\sigma \in \Phi$ and $s\in \Psi$.

Fix a finite volume $\Lambda \subset {\mathbb{N}}^{+}$. Let $\mathcal{E}=\{\ell =\langle x,y\rangle :d(x,y)=1,x,y\in {\mathbb{N}}^{+}\}$ be the set of nearest-neighbor edges with at least one endpoint in $\Lambda \subset {\mathbb{N}}^{+}$. We examine in detail the 1D-MSIM with mixed spin $(s,\frac{2t-1}{2})$ built on the lattice ${\mathbb{N}}^{+}$, having the Hamiltonian

$$\begin{array}{c}\hfill H\left(\xi \right)=-\sum _{\langle x,y\rangle \in \mathcal{E}}J\xi \left(x\right)\xi \left(y\right)-\sum _{y\in \Lambda }{h}_{\xi \left(y\right)}\xi \left(y\right),\end{array}$$

(1)

where the energy of each link tying up nearby sites is represented by the first sum, and the energy of each site is represented by the second sum.

The partial partition functions associated with the Hamiltonian (1) can be found as follows:

$$\begin{array}{c}\hfill Z_{\xi \left(x^{\left(0\right)}\right)}^{(2N+1)}(\beta ,h_{\xi \left(x^{\left(0\right)}\right)})=\sum _{{\eta }_N\in {\Xi }^{{\mathbb{N}}_{2N+1}}}exp\left\{-\beta H\left({\eta }_N\right)\right\},\end{array}$$

(2)

where $\xi \in \sigma \times s=\Xi$, $\beta =\frac{1}{k_{B}T}$ and J is the coupling constant.

Given the finite set ${\mathbb{N}}_{2N+1}$, by considering the boundary spins, or the spins on the $(2N+1)$-th level, we can construct the summation of Equation (2). Consequently, we state that the density-free energy function is

$$F(\beta ,h_{\xi \left(x\right)})=-k_{B}T\underset{N\to \infty }{lim}\frac{1}{(2N+1)}ln{Z}^{(2N+1)}(\beta ,h_{\xi \left(x\right)}),$$

(3)

where $Z^{(2N+1)}(\beta ,h_{\xi \left(x\right)})$ is the total partition function.

As a result of the derivatives of the free energy function with respect to certain parameters, many authors have studied other thermodynamic properties [21,22]. If the largest eigenvalue of the transfer matrix corresponding to the given model is ${\lambda }_{max}$, the thermodynamic quantities of this system are obtained as follows:

• Bulk free energy

$$\begin{array}{c}\hfill F(T,\beta h_{\xi \left(x\right)})=-k_{B}Tln{\lambda }_{max}.\end{array}$$

(4)

• Entropy

$$\begin{array}{c}\hfill S(T,\beta h_{\xi \left(x\right)})=-\frac{\partial F(T,\beta h_{\xi \left(x\right)})}{\partial T}.\end{array}$$

(5)

• Internal energy

$$\begin{array}{c}\hfill U(T,\beta h_{\xi \left(x\right)})=-\frac{\partial ln{\lambda }_{max}}{\partial \beta }.\end{array}$$

(6)

3. Construction of the Partition Function for the Model via Transfer Matrices

The Kramers–Wannier transfer technique may be used to construct the partition function. This method requires building a transfer matrix and determining its eigenvalues.

3.1. The Partition Function and the Boltzmann Weight

It is well known that the total partition function $Z^{(2N+1)}(\beta ,h_{\xi \left(x\right)})=Tr(-\beta H_{2N+1})$ depends on the particular realization $\{J_i,h_i\}$ (see [21,22] for details). Also, it is well known that the total partition function is equal to the Boltzmann weight’s sum over all possible states. Here, we consider blocks or configurations with $(2N+1)$ length, and under zero boundary conditions, that is, the spins designed for $m>2N+1$ are regarded as $\xi \left(m\right)=0$.

If we want to rewrite the Hamiltonian (1) for all configurations on the configuration space ${\Xi }^{2N+1}$, we can obtain the Boltzmann weight in the form:

$$\begin{array}{c}\hfill e^{-\beta H\left(\eta \right)+{\displaystyle \sum _{x\in {\mathbb{N}}_{2N+1}}}\eta \left(x\right)h_{\eta \left(x\right)}}=e^{-\beta H\left(\sigma \right)}\end{array}$$

(7)

From (1), for $i=1,2,\cdots ,N$, we denote the energy of the bond between sequential vertices $2i-1$ and $2i$ by

$$\begin{array}{cc}\hfill E({\sigma }_{2i-1},s_{2i})& =-\sum _{\langle 2i-1,2i\rangle \in \mathcal{E}}J{\sigma }_{2i-1}{s}_{2i}-\sum _{\langle 2i-1,2i\rangle \in \mathcal{E}}{h}_{s_{2i}}{s}_{2i}.\hfill \end{array}$$

(8)

Similarly, for $i=1,2,\cdots ,N$ we denote the energy of the bond between sequential vertices $2i$ and $2i+1$ by

$$\begin{array}{cc}\hfill \tilde{E}({\sigma }_{2i},s_{2i+1})& =-\sum _{\langle 2i,2i+1\rangle \in \mathcal{E}}J{s}_{2i}{\sigma }_{2i+1}-\sum _{\langle 2i,2i+1\rangle \in \mathcal{E}}{\tilde{h}}_{{\sigma }_{2i+1}}{\sigma }_{2i+1}.\hfill \end{array}$$

(9)

Note that, as can be seen in Equations (8) and (9), for $i>1$, while the oscillating magnetic fields $h_{{\sigma }_{2i-1}}$ and ${\tilde{h}}_{s_{2i}}$ fluctuate depending on where the vertex is located, the coupling constant J, which determines the energy between two successive vertices, remains constant.

It is well known that the partition function is equal to the sum of the Boltzmann weight $e^{-\beta H\left(\eta \right)}$ over all possible states in the lattice $\mathbb{N}$ [21,22]. Therefore, for the configuration $\eta \in {\Xi }^{2N+1}$, we can write the Boltzmann weight $e^{-\beta H\left(\sigma \right)}$ in the form

$$\begin{array}{ccc}\hfill e^{-\beta H\left(\eta \right)}& =& e^{-\beta H({\sigma }_1,s_2,{\sigma }_3,\cdots ,{\sigma }_{2N-1},s_{2N-1},{\sigma }_{2N+1})}\hfill \\ \hfill & =& e^{-\beta E({\sigma }_1,s_2)-\beta \tilde{E}(s_2,{\sigma }_3)-\cdots -\beta E({\sigma }_{2N-1},s_{2N})-\beta \tilde{E}(s_{2N},{\sigma }_{2N+1})}\hfill \\ \hfill & =& \left(e^{-\beta E({\sigma }_1,s_2)}\right)\left(e^{-\beta \tilde{E}(s_2,{\sigma }_3)}\right)\cdots \left(e^{-\beta E({\sigma }_{2N-1},s_{2N})}\right)\left(e^{-\beta \tilde{E}(s_{2N},{\sigma }_{2N+1})}\right).\hfill \end{array}$$

(10)

For this system, the canonical partition function can be written as

$$\begin{array}{c}\hfill Z_{2N+1}=\sum _{\eta \in {\Xi }^{2N+1}}{e}^{-\beta H\left(\eta \right)}=\sum _{\eta \in {\Xi }^{2N+1}}{e}^{-\beta \left({\sum }_{i=1}^{N}E({\sigma }_{2i-1},s_{2i})+{\sum }_{i=1}^N\tilde{E}({\sigma }_{2i},s_{2i+1})\right)}.\end{array}$$

(11)

The next step is to decompose the bonds using the bond representation and factor the Boltzmann weights into pairwise factors.

3.2. The Transfer Matrices

Here, we construct the transfer matrices corresponding to 1D-MSIM with the $(s,\frac{(2t-1)}{2})$ mixed spin on the lattice ${\mathbb{N}}^{+}$. We obtain two different transfer matrices. First, let us place the spins of the set $\Phi$ on the odd-numbered vertices of the lattice ${\mathbb{N}}^{+}$, and the spins in the set $\Psi$ on the even-numbered vertices. So, we define the entries of the transfer matrix $\mathbb{P}$ as follows:

$$\begin{array}{c}\hfill P_{{\sigma }_{2i-1}{s}_{2i}}=exp\left({\sigma }_{2i-1}{s}_{2i}\beta J+\beta s_{2i}{\tilde{h}}_{s_{2i}}\right)=e^{-\beta E({\sigma }_{2i-1},s_{2i})},\end{array}$$

(12)

where ${\sigma }_{2i-1}\in \Phi$ and $s_{2i}\in \Psi$ for $i=1,2,\cdots$.

Similarly, we can define the entries of the transfer matrix $\mathbb{Q}$ by

$$\begin{array}{c}\hfill Q_{s_{2i}{\sigma }_{2i+1}}=exp\left(s_{2i}{\sigma }_{2i+1}\beta J+\beta {\sigma }_{2i+1}{h}_{{\sigma }_{2i+1}}\right)=e^{-\beta \tilde{E}(s_{2i},{\sigma }_{2i+1})},\end{array}$$

(13)

where $s_{2i}\in \Psi$ and ${\sigma }_{2i+1}\in \Phi$ for $i=1,2,\cdots$.

From (10)–(13), we obtain

$$\begin{array}{cc}\hfill e^{-\beta H\left(\eta \right)}& =P_{{\sigma }_1{s}_2}{Q}_{s_2{\sigma }_3}{P}_{{\sigma }_3{s}_4}{Q}_{s_4{\sigma }_5}\cdots P_{{\sigma }_{2N-1}{s}_{2N}}{Q}_{s_{2N}{\sigma }_{2N+1}}.\hfill \end{array}$$

(14)

One can easily see that the matrix $\mathbb{P}$ has Card$(\Phi )\times$ Card$(\Psi )$ dimensions while matrix $\mathbb{Q}$ has Card$(\Psi )\times$ Card$(\Phi )$ dimensions. Assume that $\mathbb{P}\mathbb{Q}=\mathbb{M}$. The total partition function can be written as the trace of the product of $2N+1$ random transfer matrices.

3.3. 1D-MSIM with Mixed Spin-$(1,1/2)$

Due to the difficulty of computation of the transfer matrix and its eigenvalues for $s>1$ and $t>1$, we restrict ourselves to the mixed spin $(1,1/2)$. In this subsection, we compute the thermodynamic quantities of 1D-MSIM with mixed spin $(1,1/2)$. So, we consider the mixed spins as $\Phi =\{-1,0,1\}$ and $\Psi =\{-\frac{1}{2},\frac{1}{2}\}$.

Firstly, let us construct the transfer matrix as $\mathbb{P}=\left(P_{{\sigma }_{2i-1},s_{2i}}\right)$, for ${\sigma }_{2i-1}\in \{-1,0,1\}$ and $s_{2i}\in \{-\frac{1}{2},\frac{1}{2}\}$. So, we have

$$\begin{array}{c}\hfill \mathbb{P}=\left(\begin{array}{cc}{e}^{\frac{\beta }{2}\left(J-h_{-\frac{1}{2}}\right)}& e^{-\frac{\beta }{2}\left(J-h_{\frac{1}{2}}\right)}\\ e^{-\frac{1}{2}\beta h_{-\frac{1}{2}}}& e^{\frac{1}{2}\beta h_{\frac{1}{2}}}\\ e^{-\frac{\beta }{2}\left(J+h_{-\frac{1}{2}}\right)}& e^{\frac{\beta }{2}\left(J+h_{\frac{1}{2}}\right)}\end{array}\right).\end{array}$$

(15)

Secondly, let us place the spins of the set $\{-\frac{1}{2},\frac{1}{2}\}$ on the even-numbered vertices of the lattice ${\mathbb{N}}^{+}$, and the spins in the set $\{-1,0,1\}$ on the odd-numbered vertices (see Figure 1). We can define the second transfer matrix $\mathbb{Q}=\left(Q_{s_{2i},{\sigma }_{2i+1}}\right)$, for $s_{2i}\in \{-\frac{1}{2},\frac{1}{2}\}$ and ${\sigma }_{2i+1}\in \{-1,0,1\}$ as

$$\begin{array}{c}\hfill \mathbb{Q}=\left(\begin{array}{ccc}{e}^{\frac{\beta }{2}\left(J-2h_{-1}\right)}& 1& e^{-\frac{\beta }{2}\left(J-2h_1\right)}\\ e^{-\frac{\beta }{2}\left(J+2h_{-1}\right)}& 1& e^{\frac{\beta }{2}\left(J+2h_1\right)}\end{array}\right).\end{array}$$

(16)

3.4. Investigation of Thermodynamic Quantities in the Translation-Invariant Case

Assume that $h={\tilde{h}}_j=h_i$ for all $j\in \{-\frac{1}{2},\frac{1}{2}\}$ and $i\in \{-1,0,1\}$. This case is called the translation-invariant property. Thus, the transfer matrices $\mathbb{P}$ and $\mathbb{Q}$ given in (15) and (16) are obtained as follows

$$\begin{array}{c}\hfill \mathbb{P}=\left(\begin{array}{cc}{e}^{\frac{\beta }{2}(J-h)}& e^{-\frac{\beta }{2}(J-h)}\\ e^{-\frac{1}{2}\beta h}& e^{\frac{1}{2}\beta h}\\ e^{-\frac{\beta }{2}(J+h)}& e^{\frac{\beta }{2}(J+h)}\end{array}\right),\end{array}$$

(17)

$$\begin{array}{c}\hfill \mathbb{Q}=\left(\begin{array}{ccc}{e}^{\frac{\beta }{2}(J-2h)}& 1& e^{-\frac{\beta }{2}(J-2h)}\\ e^{-\frac{\beta }{2}(J+2h)}& 1& e^{\frac{\beta }{2}(J+2h)}\end{array}\right).\end{array}$$

(18)

Case I

When we multiply the transfer matrices $\mathbb{P}$ and $\mathbb{Q}$ given in Equations (17) and (18), we obtain the square matrix:

$$\begin{array}{c}\hfill \mathbb{M}=\mathbb{P}\mathbb{Q}=\left(\begin{array}{ccc}{e}^{-\frac{1}{2}(3h-2J)\beta }+e^{-\frac{1}{2}(h+2J)\beta }& e^{-\frac{1}{2}(J-h)\beta }+e^{\frac{1}{2}(J-h)\beta }& e^{\frac{h\beta }{2}}\left(1+e^{h\beta }\right)\\ e^{-\frac{1}{2}(3h-J)\beta }+e^{-\frac{1}{2}(h+J)\beta }& e^{-\frac{h\beta }{2}}\left(1+e^{h\beta }\right)& e^{\frac{1}{2}(h-J)\beta }+e^{\frac{1}{2}(3h+J)\beta }\\ e^{-\frac{3h\beta }{2}}\left(1+e^{h\beta }\right)& e^{-\frac{1}{2}(h+J)\beta }+e^{\frac{1}{2}(h+J)\beta }& e^{\frac{1}{2}(h-2J)\beta }+e^{\frac{1}{2}(3h+2J)\beta }\end{array}\right).\end{array}$$

(19)

After some tedious and complicated algebraic manipulations, the eigenvalues of the matrix $\mathbb{M}$ are obtained as:

$$\left\{\begin{array}{c}{\lambda }_1=0,\hfill \\ {\lambda }_2=\frac{1}{2e^{\frac{\beta }{2}(3h+2J)}}\left(A(\beta ,h,J)-\sqrt{{\left(A(\beta ,h,J)\right)}^2-4B(\beta ,h,J)}\right),\hfill \\ {\lambda }_3=\frac{1}{2e^{\frac{\beta }{2}(3h+2J)}}\left(A(\beta ,h,J)+\sqrt{{\left(A(\beta ,h,J)\right)}^2-4B(\beta ,h,J)}\right),\hfill \end{array}\right.$$

where

$$\begin{array}{cc}\hfill A(\beta ,h,J)& =\left(1+e^{h\beta }\right)\left(e^{h\beta }+e^{2J\beta }+e^{h\beta +J\beta }-e^{h\beta +2J\beta }+e^{2h\beta +2J\beta }\right),\hfill \\ \hfill B(\beta ,h,J)& =e^{2h\beta }{\left(e^{J\beta }-1\right)}^2(e^{J\beta }\left(1+e^{2h\beta }\right)+e^{h\beta }{\left(1+e^{J\beta }\right)}^2).\hfill \end{array}$$

For the sake of simplicity, substitute the variables $\theta =e^{\frac{\beta J}{2}}$ and $v=e^{\frac{\beta h}{2}}$ to obtain

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\lambda }_1=0,\hfill \\ {\lambda }_2=\frac{\left(1+v^2\right)\left({\theta }^4+v^2\left(1+{\theta }^2-{\theta }^4\right)+v^4{\theta }^4\right)-\sqrt{{\left(\left(1+v^2\right)\left({\theta }^4+v^2\left(1+{\theta }^2-{\theta }^4\right)+v^4{\theta }^4\right)\right)}^2-4v^4{\left({\theta }^2-1\right)}^2\left({\theta }^2(1+v^4)+v^2{(1+{\theta }^2)}^2\right)}}{2v^3{\theta }^2},\hfill \\ {\lambda }_3=\frac{\left(1+v^2\right)\left({\theta }^4+v^2\left(1+{\theta }^2-{\theta }^4\right)+v^4{\theta }^4\right)+\sqrt{{\left(\left(1+v^2\right)\left({\theta }^4+v^2\left(1+{\theta }^2-{\theta }^4\right)+v^4{\theta }^4\right)\right)}^2-4v^4{\left({\theta }^2-1\right)}^2\left({\theta }^2(1+v^4)+v^2{(1+{\theta }^2)}^2\right)}}{2v^3{\theta }^2}.\hfill \end{array}\right.\end{array}$$

(20)

By computing the trace of the matrix product $\mathbb{M}$, we can determine the partition function as

$$\begin{array}{c}\hfill Z^{(2N+1)}(\beta ,h_{\xi \left(x\right)})=Tr\left[M^N\right]=\sum _{i=1}^3{\lambda }_i^N,\end{array}$$

(21)

where ${\lambda }_i$ denotes the eigenvalues of the transfer matrix $\mathbb{M}$ for $i=1,2,3.$

From (11), we obtain

$$\begin{array}{c}\hfill Z^{(2N+1)}(\beta ,h_{\xi \left(x\right)})=Tr\left({\mathbb{M}}^N\right)=\sum _{j=1}^t{\lambda }_{t,j}^N={\lambda }_{max}^N\left(1+\sum _{j=1}^{t-1}{\left(\frac{{\lambda }_{t,j}}{{\lambda }_{max}}\right)}^N\right),\end{array}$$

(22)

where ${\lambda }_{t,j}$ denote the eigenvalues of transfer matrix $\mathbb{M}$.

So, one has $Z^{(2N+1)}(\beta ,h_{\xi \left(x\right)})={\lambda }_{max}^N$ for $N\to \infty$. As it is well known, the model’s critical behavior manifests itself in the thermodynamic limit as $N\to \infty$; therefore, the largest eigenvalue ${\lambda }_{max}={\lambda }_3$ of the transfer matrix $\mathbb{M}$ given in (19) stands as the only indicator of the bulk free energy:

$$\begin{array}{cc}\hfill F(\beta ,h)& =-\frac{1}{\beta }\underset{N\to \infty }{lim}\frac{1}{(2N+1)}ln{Z}^{(2N+1)}(\beta ,h)\hfill \\ \hfill & =-\frac{1}{\beta }\underset{N\to \infty }{lim}\frac{1}{(2N+1)}ln{\lambda }_3^N=-\frac{1}{2\beta }ln{\lambda }_3\hfill \\ \hfill & =-\frac{1}{2\beta }ln\left(\frac{A(\beta ,h,J)+\sqrt{{\left(A(\beta ,h,J)\right)}^2-4B(\beta ,h,J)}}{2e^{\frac{\beta }{2}(3h+2J)}}\right).\hfill \end{array}$$

(23)

3.5. Behavior of the Thermodynamic Quantities of 1D-MSIM with Mixed Spin-(s,(2t − 1)/2) in the Absence of a Magnetic Field

For $h=0$, we obtain the transfer matrix $\mathbb{M}$ given in (19) as

$$\begin{array}{c}\hfill \widehat{\mathbb{M}}=\left(\begin{array}{ccc}2cosh\left(J\beta \right)& 2cosh\left(\frac{J\beta }{2}\right)& 2\\ 2cosh\left(\frac{J\beta }{2}\right)& 2& 2cosh\left(\frac{J\beta }{2}\right)\\ 2& 2cosh\left(\frac{J\beta }{2}\right)& 2cosh\left(J\beta \right)\end{array}\right).\end{array}$$

(24)

The set of eigenvalues of the matrix $\widehat{\mathbb{M}}$ is obtained as follows:

$$\{{\widehat{\lambda }}_1=0,{\widehat{\lambda }}_2=-2+2cosh\left(J\beta \right),{\widehat{\lambda }}_3=4+2cosh\left(J\beta \right)\}.$$

It is clear that ${\widehat{\lambda }}_3={\widehat{\lambda }}_{max}.$ From the formula in (23), we obtain

$$\begin{array}{c}\hfill F(\beta ,0)=-\frac{1}{2\beta }ln\left(4+2cosh\left(J\beta \right)\right).\end{array}$$

(25)

Figure 2 shows the graph of free energy $F(\beta ,0)$ given in (25) as a function of $\beta$ for $J=0.2$. From Figure 2, it can be seen that the free energy is approaching zero for $T\to 0^{+}$. In Figure 2, the oscillating magnetic fields $h_{s_{2i}}$ and ${\tilde{h}}_{{\sigma }_{2i+1}}$ are assumed to be zero. If $h_{s_{2i}}$ and ${\tilde{h}}_{{\sigma }_{2i+1}}$ are taken to be nonzero, then, by using the parametric representation, the free energy function can be plotted in the form $h\to \left(\beta \right(h),F(h\left)\right)$.

In the absence of an external magnetic field, from (5), the entropy of the model is obtained as

$$\begin{array}{cc}\hfill S(\beta ,0)& =-\frac{\partial F(\beta ,0)}{\partial T}=-\frac{\frac{\partial F(\beta ,0)}{\partial \beta }}{\frac{\partial T}{\partial \beta }}\hfill \\ \hfill & =\frac{(-2ln(4+2cosh\left(J\beta \right))-cosh\left(J\beta \right)ln(4+2cosh\left(J\beta \right))+J\beta sinh\left(J\beta \right))k_B}{2(2+cosh(J\beta \left)\right)}.\hfill \end{array}$$

(26)

Similarly, from (6), we obtain the internal energy as

$$\begin{array}{c}\hfill U(\beta ,0)=\frac{-2ln(4+2cosh(J\beta \left)\right)-cosh\left(J\beta \right)ln(4+2cosh(J\beta \left)\right])+J\beta sinh(J\beta \left]\right)}{2{\beta }^2(2+cosh\left(J\beta \right))}.\end{array}$$

(27)

The changes in entropy of 1D-MSIM with mixed spin $(1,1/2)$ are given in Figure 3. Results show that the internal energy converges to zero at the higher temperature values, while becoming a constant 0 as $T\to 0$. When the temperature increases, the entropy increases until it reaches $ln\left(\sqrt{6}\right)$ and goes to zero at low temperature values.

Case II

Note that while successively placing spins on the vertices of a one-dimensional lattice, if the elements of $\Psi$ are placed in the first vertex of the lattice $\mathbb{N}$, then a transfer matrix having $2\times 2$ dimensions is obtained, and the eigenvalues of the transfer matrix are the same as those of the previous matrix $\mathbb{M}$.

$$\begin{array}{c}\hfill \tilde{\mathbb{M}}=\mathbb{Q}\mathbb{P}=\left(\begin{array}{cc}{e}^{-\frac{h\beta }{2}}+e^{-\frac{1}{2}(3h-2J)\beta }+e^{\frac{1}{2}(h-2J)\beta }& e^{\frac{h\beta }{2}}+e^{-\frac{h\beta }{2}}+e^{\frac{3h\beta }{2}}\\ e^{-\frac{h\beta }{2}}+e^{-\frac{3h\beta }{2}}+e^{\frac{h\beta }{2}}& e^{\frac{h\beta }{2}}+e^{-\frac{1}{2}(h+2J)\beta }+e^{\frac{1}{2}(3h+2J)\beta }\end{array}\right).\end{array}$$

(28)

If we substitute the variables $\theta =e^{\frac{\beta J}{2}}$ and $v=e^{\frac{\beta h}{2}}$ for simplicity, we obtain

$$\begin{array}{cc}\hfill \tilde{M}=\mathbb{Q}\mathbb{P}& =\left(\begin{array}{ccc}\frac{\theta }{v^2}& 1& \frac{v^2}{\theta }\\ \frac{1}{\theta v^2}& 1& \theta v^2\end{array}\right).\left(\begin{array}{cc}\frac{\theta }{v}& \frac{v}{\theta }\\ \frac{1}{v}& vs.\\ \frac{1}{\theta v}& \theta v\end{array}\right)\hfill \\ \hfill & =\left(\begin{array}{cc}\frac{\left(v^2-v\theta +{\theta }^2\right)\left(v^2+v\theta +{\theta }^2\right)}{v^3{\theta }^2}& \frac{\left(1-v+v^2\right)\left(1+v+v^2\right)}{v}\\ \frac{\left(1-v+v^2\right)\left(1+v+v^2\right)}{v^3}& \frac{\left(1-v\theta +v^2{\theta }^2\right)\left(1+v\theta +v^2{\theta }^2\right)}{v{\theta }^2}\end{array}\right).\hfill \end{array}$$

(29)

We obtain the set of eigenvalues of the matrix $\tilde{\mathbb{M}}$ given in $\left(29\right)$ as

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\tilde{\lambda }}_1=\frac{\left(1+v^2\right)\left({\theta }^4+v^2\left(1+{\theta }^2-{\theta }^4\right)+v^4{\theta }^4\right)-\sqrt{{\left(\left(1+v^2\right)\left({\theta }^4+v^2\left(1+{\theta }^2-{\theta }^4\right)+v^4{\theta }^4\right)\right)}^2-4v^4{\left({\theta }^2-1\right)}^2\left({\theta }^2(1+v^4)+v^2{(1+{\theta }^2)}^2\right)}}{2v^3{\theta }^2},\hfill \\ {\tilde{\lambda }}_2=\frac{\left(1+v^2\right)\left({\theta }^4+v^2\left(1+{\theta }^2-{\theta }^4\right)+v^4{\theta }^4\right)+\sqrt{{\left(\left(1+v^2\right)\left({\theta }^4+v^2\left(1+{\theta }^2-{\theta }^4\right)+v^4{\theta }^4\right)\right)}^2-4v^4{\left({\theta }^2-1\right)}^2\left({\theta }^2(1+v^4)+v^2{(1+{\theta }^2)}^2\right)}}{2v^3{\theta }^2}.\hfill \end{array}\right.\end{array}$$

(30)

It should be noted here that the eigenvalues of the matrices $\mathbb{M}$ and $\tilde{\mathbb{M}}$ are the same, except for 0. Therefore, it is obtained that ${\lambda }_3={\tilde{\lambda }}_2$ (see Equations (20) and (30)). We can choose any spin on the starting vertex, while mixed spins are placed on consecutive vertices of the lattice. From (20) and (30), it is clear that ${\tilde{\lambda }}_{max}={\tilde{\lambda }}_2$.

3.6. Magnetization and Magnetic Susceptibility

In this subsection, we construct magnetization and magnetic susceptibility by means of the eigenvalue ${\tilde{\lambda }}_2$ in the Formula (30).

Let us consider the reduced nearest neighbor spin–spin coupling interactions $K=\frac{\beta J}{2}=ln\left(\theta \right)$ and the reduced magnetic field $H=\frac{\beta h_{\xi \left(x\right)}}{2}=ln\left(v\right)$, we write the magnetization as

$$\begin{array}{c}\hfill m(T,H)=-\beta \frac{\partial F(T,H)}{\partial a{\tilde{\lambda }}_2}\partial H,\end{array}$$

(31)

and susceptibility

$$\begin{array}{c}\hfill \chi (T,H)=-{\beta }^2\frac{{\partial }^{2}F(T,H)}{\partial H^2}=\beta \frac{\partial }{\partial H}\left(\frac{1}{{\tilde{\lambda }}_2}\frac{\partial {\tilde{\lambda }}_2}{\partial H}\right).\end{array}$$

(32)

We do not give the exact expressions of the formulas for the magnetization and susceptibility here, because the operations are excessively long and complex. We numerically examine the behavior of these two quantities as functions of h and T.

The Mathematica software (Version 8.0, Wolfram Research, Inc., Champaign, IL, USA) [23] has been used to perform the calculations and to plot the figures. A three-dimensional plot of $m(T,H)$ is given in Figure 4. Taking into account the eigenvalue ${\tilde{\lambda }}_2$ in the Formula (30), for $J=-2$, one can easily see that, as the temperature T increases, the smoothness of the function decreases and exhibits a behavior similar to the step function. On the other hand, for $J=2$, the surface of the function becomes smooth. In contrast to the single-spin Ising chain’s typical single step [21], in Figure 4 (left), the magnetization graph shows three distinct steps at low temperatures. This phenomenon may be explained by spins adopting distinct states within even and odd sublattices, or by all spins assuming the same state.

A three-dimensional plot of $\chi (T,H)$ is given in Figure 5. Figure 5 (Left) and (Right) show the behavior of the susceptibility function given in Equation (32). For $J=-2$, three stacks resembling a boot are observed, and a stack appears for $J=2$ in the chosen region $T\in [0.001,7]$, $h\in [-6,6]$. As seen in Figure 5 (Left), while three different susceptibility peaks appear for $J=-2$ at low temperatures, the susceptibility peaks disappear as the temperature increases. For $J=2$, only a single susceptibility peak is observed at low temperatures (see Figure 5 (Right)).

4. Nonexistence Phase Transition in the Absence of the External Magnetic Field

In this section, we study 1D-MSIM with the mixed spin-(1,1/2) employing the ERR.

Suppose that $W^{\left(n\right)}=\frac{{\tilde{Z}}_{-\frac{1}{2}}^{\left(n\right)}}{{\tilde{Z}}_{\frac{1}{2}}^{\left(n\right)}}$, $X^{(n-1)}=\frac{Z_{-1}^{(n-1)}}{Z_0^{(n-1)}}$ and $Y^{(n-1)}=\frac{Z_1^{(n-1)}}{Z_0^{(n-1)}}$, from Equations (43) and (44), we obtain

$$\begin{array}{cc}\hfill W^{\left(n\right)}& =\frac{{\theta }^2{X}^{(n-1)}+\theta +Y^{(n-1)}}{X^{(n-1)}+\theta +{\theta }^2{Y}^{(n-1)}},\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill X^{(n-1)}& =\frac{{\theta }^2{W}^{(n-2)}+1}{\theta (W^{(n-2)}+1)},\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill Y^{(n-1)}& =\frac{W^{(n-2)}+{\theta }^2}{\theta (W^{(n-2)}+1)},\hfill \end{array}$$

(35)

where $\theta =e^{\frac{\beta J}{2}}$.

If we assume $x=\underset{n\to \infty }{lim}{X}^{(n-1)}$, $y=\underset{n\to \infty }{lim}{Y}^{(n-1)}$ and $w=\underset{n\to \infty }{lim}{W}^{\left(n\right)}$, then we obtain a new three-dimensional rational dynamical system (3D-RDS) as

$$\begin{array}{cc}\hfill w& =\frac{{\theta }^{2}x+\theta +y}{x+\theta +{\theta }^{2}y},\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill x& =\frac{{\theta }^{2}w+1}{\theta (w+1)},\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill y& =\frac{w+{\theta }^2}{\theta (w+1)}.\hfill \end{array}$$

(38)

If the values of x and y in Equations (37) and (38) are substituted into Equation (36), one obtains the following rational recursive equation $f:{\mathbb{R}}_{+}\to {\mathbb{R}}_{+}$:

$$\begin{array}{c}\hfill f\left(w\right):=\frac{3{\theta }^2+w\left(1+{\theta }^2+{\theta }^4\right)}{\left(1+{\theta }^2+{\theta }^4\right)+3w{\theta }^2}=w.\end{array}$$

(39)

From (39), we obtain a second-order equation

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

(40)

From Equation (40), it is clear that the function f given in (39) has only one positive fixed point, so there is no phase transition for the given model.

The graphs of the function f for antiferromagnetic ($\theta =0.421$) and ferromagnetic ($\theta =3.421$) values are plotted in Figure 6. The diagrams show that the function in both cases has just two fixed points. Additionally, there is only one positive fixed point. As a result, there is just one Gibbs measure and no phase transition in the model. As it is well known, the classical single spin 1D Ising model’s phase transition issue has attracted the interest of statistical mechanics researchers for over a century, and it has been established that there is no phase transition [1]. In this present paper, we demonstrated that the mixed spin 1D Ising model has no phase transition as an example.

Chaoticity of the Model

A dynamical system’s chaotic behavior is determined by how sensitively it depends on the initial conditions [24,25]. It has long been a challenge to see whether a model’s phase transition and the chaotic behavior of the corresponding dynamical system are related [26]. In this subsection, we investigate the chaoticity of the 3D-RDS given in (36)–(38). With the help of the Lyapunov exponent, we numerically study the model’s chaotic behavior.

From the definition of the Lyapunov exponent, we obtain

$$\begin{array}{cc}\hfill \lambda & =\frac{1}{n}ln\left|\prod _{i=0}^{n-1}{f}_{\theta }^{\prime }\left(w_i\right)\right|=\frac{1}{n}\sum _{i=0}^{n-1}ln\left|f_{\theta }^{\prime }\left(w_i\right)\right|\hfill \\ \hfill & =\frac{1}{n}\sum _{i=0}^{n-1}ln\left|\frac{{\left({\theta }^2-1\right)}^2\left(1+4{\theta }^2+{\theta }^4\right)}{{\left(1+(1+3w_i){\theta }^2+{\theta }^4\right)}^2}\right|.\hfill \end{array}$$

(41)

From (41), we have

$$\begin{array}{cc}\hfill \lambda \left(\theta \right):& =\underset{n\to \infty }{lim}\frac{1}{n}\sum _{i=0}^{n-1}ln\left|\frac{{\left({\theta }^2-1\right)}^2\left(1+4{\theta }^2+{\theta }^4\right)}{{\left(1+(1+3w_i){\theta }^2+{\theta }^4\right)}^2}\right|.\hfill \end{array}$$

(42)

Figure 7 shows the Lyapunov exponent of the dynamical system corresponding to the 1D-MSIM with mixed spin-(1,1/2). It is seen that the Lyapunov exponent is always negative in the ferromagnetic region. Therefore, the rational function $f\left(w\right)$ given in (39) is periodic. The behavior of the Lyapunov exponent $\lambda \left(\theta \right)$ in the antiferromagnetic region can also be examined.

5. The Average Magnetization for the Mixed Spin-(1,1/2) Ising Model

In this section, by using the exact recursion relations (ERRs) (see [9]), we obtain the partial partition functions of the 1D-MSIM with mixed spin-$(s,\frac{(2t-1)}{2})$. Contrary to the previous sections, here we place the spins of $\Psi$ in the first vertex of the lattice and the spins of $\Phi$ in the second vertex of the lattice, while placing spins at the vertices of the lattice ${\mathbb{Z}}_{+}$, so for $n=0,1,2,\cdots$, we obtain $s_{2n-1}\in \Psi$ and ${\sigma }_{2n}\in \Phi$. With the help of the cavity method (see [14,15,16] for details), we obtain the partial partition functions as follows:

$$\begin{array}{cc}\hfill {\tilde{Z}}_{s_1}^{\left(n\right)}& =\sum _{{\sigma }_2\in \Phi }exp\left(\beta (J{s}_1{\sigma }_2+h_{{\sigma }_2})\right)Z_{{\sigma }_2}^{(n-1)},\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill Z_{{\sigma }_2}^{(n-1)}& =\sum _{s_3\in \Psi }exp\left(\beta (J{\sigma }_2{s}_3+{\tilde{h}}_{s_3})\right){\tilde{Z}}_{s_3}^{(n-2)}.\hfill \end{array}$$

(44)

The Average Magnetization

In this subsection, assuming there is an external magnetic field, we obtain the magnetization equations for spins s and $\frac{(2t-1)}{2}$, respectively, as follows.

$$\begin{array}{c}\hfill M_{1/2}^{\left(n\right)}=\frac{{\displaystyle \sum _{s_1\in \Psi }}{s}_{1}exp\left(\beta (J{s}_1{\sigma }_2+h_{{\sigma }_2})\right){\tilde{Z}}_{s_1}^{\left(n\right)}}{{\displaystyle \sum _{s_1\in \Psi }}exp\left(\beta (J{s}_1{\sigma }_2+h_{{\sigma }_2})\right){\tilde{Z}}_{s_1}^{\left(n\right)}},\end{array}$$

(45)

$$\begin{array}{c}\hfill M_1^{(n-1)}=\frac{{\displaystyle \sum _{{\sigma }_2\in \Phi }}{\sigma }_{2}exp\left(\beta (J{\sigma }_2{s}_3+{\tilde{h}}_{s_3})\right)Z_{{\sigma }_1}^{(n-1)}}{{\displaystyle \sum _{{\sigma }_2\in \Phi }}exp\left(\beta (J{\sigma }_2{s}_3+{\tilde{h}}_{s_3})\right)Z_{{\sigma }_2}^{(n-1)}}.\end{array}$$

(46)

We are only concerned with the magnetization for mixed spin-(1,1/2) here, since the operations for obtaining partition functions are excessively long for $s>1$ and $t>1$. So, from (45) and (46), we obtain

$$\begin{array}{cc}\hfill M_{1/2}^{\left(n\right)}& =\frac{e^{\frac{h\beta }{2}}\left(e^{\left(-\frac{h}{2}-\frac{J}{2}\right)\beta }-e^{\left(-\frac{h}{2}+\frac{J}{2}\right)\beta }\right)\left(e^{-h\beta }-W^{\left(n\right)}\right)}{e^{-3h\beta }\left(e^{\frac{1}{2}(4h-J)\beta }+e^{\frac{1}{2}(4h+J)\beta }+e^{2h\beta }\right)+e^{-\frac{J\beta }{2}}\left(1+e^{\frac{J\beta }{2}}+e^{J\beta }\right)W^{\left(n\right)}},\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill M_1^{(n-1)}& =\frac{\left(e^{J\beta }-1\right)\left(e^{2h\beta }{Y}^{(n-1)}-X^{(n-1)}\right)}{2\left(2e^{h\beta +\frac{J\beta }{2}}+\left(1+e^{J\beta }\right)\left(X^{(n-1)}+e^{2h\beta }{Y}^{(n-1)}\right)\right)}.\hfill \end{array}$$

(48)

The following equation results from the simple substitutions of the variables $\theta =e^{\frac{J\beta }{2}}$ and $v=e^{\frac{h\beta }{2}}$

$$\begin{array}{cc}\hfill M_{1/2}^{\left(n\right)}& =\frac{\left(v^2{W}^{\left(n\right)}-1\right)({\theta }^2-1)}{\left(1+v^2{W}^{\left(n\right)}\right)\left(1+\theta +{\theta }^2\right)},\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill M_1^{(n-1)}& =-\frac{\left(X^{(n-1)}-v^4{Y}^{(n-1)}\right)({\theta }^2-1)}{2\left(X^{(n-1)}\left(1+{\theta }^2\right)+2v^2\theta +v^4{Y}^{(n-1)}\left({\theta }^2+1\right)\right)}.\hfill \end{array}$$

(50)

Now, let us assume that the external magnetic field h ($v=1$) is zero. In this case, after some calculations, as in Ref. [14], we can show that the fixed point for the 3D-RDS given in (36)–(38) is only $x=y=\left(\frac{1+{\theta }^2}{2\theta }\right)$ and $w=1$. So, we obtain zero average magnetization around this fixed point. That is, we have

$$M_{1/2}^{\left(n\right)}=M_1^{(n-1)}=0.$$

6. Conclusions

There are many methods for calculating the free energy of lattice models on the Bethe lattices or the d-dimensional integer lattice ${\mathbb{Z}}^d$ ($d\ge 1$). Some of these are the exact recursion relations [9], the cavity method [18], and vector-valued boundary conditions. One of the most widely used methods to examine thermodynamic quantities corresponding to a 1D Ising model is the transfer matrix technique [21,22]. In this present work, we constructed the partial partition functions associated with 1D-MSIM having mixed spin-$(s,\frac{(2t-1)}{2})$ via the transfer matrix. We have computed some thermodynamic quantities such as the free energy, entropy, magnetization, and susceptibility of the model. To the best of my knowledge, the thermodynamic properties of the 1D Ising model with the single spin have been investigated in several studies to date. Using the transfer matrix method, we evaluated the thermodynamic features of the 1D-MSIM for the first time. For the single spin-Isign model, Amin et al. [21] plotted the magnetization and susceptibility in two dimensions.

In our previous studies, we proved that there is a phase transition for the mixed-spin Ising model on the CT using many methods [11,12,14]. As it is known, there is no phase transition for a single spin 1D Ising model [1]. In this study, we investigated how mixed spins affect thermodynamic quantities and the existence of the phase transition. We demonstrated the uniqueness of limiting Gibbs measures associated with 1D-MSIM having mixed spin-(1,1/2). We have shown that if the external magnetic field is zero, for the aforementioned model, there is no phase transition. Our next research will deal with other issues in statistical mechanics. By considering the approach given by Albayrak [9], the isothermal entropy change of the average magnetization for 1D-MSIM will be analyzed.

The findings obtained in our present article exhibit different behaviors from the results given in previous studies [21,22]. Frankly, we cannot comment on the physical meaning of such behavior. These topics are covered in introductory courses in statistical mechanics at the undergraduate level. Therefore, we believe that the results of our present study will be of interest to a wide readership of statistical mechanics.