GL(2,R) in the Finite One-Dimensional Ising Model with Nonuniform Couplings

1. Introduction

In the current article we demonstrate how some properties of the $\mathrm{GL}(2,\mathbb{R})$ matrices, together with those of the Chebyshev polynomials, can be used to study several Ising-type models with non-uniform couplings.

Variants of the one-dimensional Ising model with non-uniform couplings have been investigated extensively, both because of their physical relevance (see, e.g. [1] and refs. therein), and because of their rich mathematical structure ( see, e.g. [2], and refs. therein). The main difficulty stems from the noncommutativity of the transfer matrices subsequence associated with different lattice sites. Within the standard transfer-matrix approach these matrices generally do not commute and therefore cannot be diagonalized simultaneously, which prevents the derivation of a closed-form expression for the partition function. For constant or periodic couplings, the partition function can be expressed in closed form from the trace of a product of simultaneously diagonalizable transfer matrices. In the context of quasicrystals, quasi-periodic couplings present a particularly interesting and nontrivial case.

For ferromagnetic Ising model chains with aperiodic couplings (say, Fibonacci), the Lee–Yang theorem implies that the partition-function zeros lie on the unit circle. In contrast to periodic systems, however, these zeros are not continuously distributed, but form a Cantor-like structure with gaps. The associated trace-map dynamics has been studied in [3,4], while a rigorous thermodynamic-limit analysis for the Fibonacci substitution sequence is given in [5]. Related rigorous results have also been obtained using the Szegő recursion for orthogonal polynomials on the unit circle and the gap labeling theorem for CMV matrices [2]. For example, the trace-map technique typically provides only an efficient framework for the numerical evaluation of the free energy in the thermodynamic limit [3]. It should be emphasized, however, that exact closed-form expressions for standard thermodynamic quantities remain relatively rare.

In the present work, we do not address these important algebraic questions, which generally require rather sophisticated mathematical methods. Instead, by suitably relaxing the inhomogeneity conditions, one can still derive analytical results even for finite number of dynamical variables N in the chain. In particular, our approach is based on exploring the algebraic properties of Chebyshev polynomials. It should be emphasized, however, that exact analytical solutions leading to closed-form expressions for standard thermodynamic quantities still remain relatively rare in the literature. Within this framework, the algebraic structure of powers of $2\times 2$ matrices provides an efficient alternative to the direct diagonalization of transfer matrices. This simplification is a consequence of the fact that the algebra generated by a $2\times 2$ matrix $\mathbf{A}$ is fully governed by its two fundamental invariants, $\mathrm{Tr}\left(\mathbf{A}\right)$ and $\det \left(\mathbf{A}\right)$, via the Cayley–Hamilton theorem. Such an approach avoids the explicit spectral decomposition of the matrices, appearing under the trace, providing a closed algebraic framework for their evaluation. This observation naturally leads to a representation of matrix powers in terms of Chebyshev polynomials. In particular, the recursion relations satisfied by $2\times 2$ matrix powers can be mapped onto those of Chebyshev polynomials of the second kind, providing a closed algebraic framework for their evaluation.

The paper is organized as follows. In Section 2, the corresponding theorems make the above correspondence precise. Section 3 and Section 4 are devoted to the study of the thermodynamic quantities of the Ising model with periodically distributed (p)-site unit cells. These results provide the physical realization of the preceding mathematical framework and are developed for general (p)-site unit cells. Section 5 is devoted to a fluctuation induced well known Casimir force.Section 6 provides a brief overview of the Schottky anomaly, which is subsequently investigated for the models considered in this work. In Section 7, we consider the standard one-dimensional Ising model with periodic boundary conditions ((p=1))the sequence of single letter (AAA…), see Figure 2, uniform interaction coupling (K), and external magnetic field (h), allowing for ferromagnetic ((K>0)) or antiferromagnetic ((K<0)) interactions, as well as positive and negative magnetic fields. In Section 8, we investigate an Ising model with nearest-neighbor interactions modulated by a periodic word of period $p=2$, namely $ABAB\cdots$, over the alphabet $\{A,B\}$, see Figure 5. The model consists of alternating ferromagnetic and antiferromagnetic interactions, and we derive the corresponding free energy in the presence of a constant magnetic field.Finally, Section 9 summarizes the main results and offers some concluding remarks.

2. Algebraic Identities for Matrices in $\mathrm{GL}(2,\mathbb{R})$

In what follows, we will need the properties of the general linear group of $2\times 2$ invertible matrices $\mathbf{A}$ over the real numbers $\mathrm{GL}(2,\mathbb{R})$, i.e. $\mathrm{GL}(2,\mathbb{R})\equiv \left\{\mathbf{A}\in M(2,\mathbb{R})\left|\det \right(A)\ne 0\right\},$ where $M(2,\mathbb{R})$ is the set of all $2\times 2$ matrices with entries being real numbers. The current section delivers several algebraic identities for matrices $\mathbf{A}\in \mathrm{GL}(2,\mathbb{R})$. Furthermore, in applications, we are going to consider throughout the paper, we assume that all matrices are entrywise positive. Since the characteristic polynomial of every such matrix has positive discriminant, every matrix considered in this paper has two distinct real eigenvalues. Its determinant can be, however, both positive and negative.

Theorem 1

(Chebyshev representation of $2\times 2$ matrix powers). Let $\mathbf{A}\in \mathrm{GL}(2,\mathbb{R})$. For any choice of $\sqrt{\det \left(\mathbf{A}\right)}\in \mathbb{C}$, define

$$x={\displaystyle \frac{\mathrm{Tr}\left(\mathbf{A}\right)}{2\sqrt{\det \left(\mathbf{A}\right)}}}.$$

Let $U_n$ denote the Chebyshev polynomials of the second kind. Then, for all integers $n\ge 1$,

$${\mathbf{A}}^n={\left(\sqrt{\det \left(\mathbf{A}\right)}\right)}^n\left[U_{n-1}\left(x\right){\displaystyle \frac{\mathbf{A}}{\sqrt{\det \left(\mathbf{A}\right)}}}-U_{n-2}\left(x\right)\mathbf{I}\right],$$

(1)

where $U_{-1}\left(x\right)=0$. Moreover, the right-hand side is independent of the choice of the square root of $\det \left(\mathbf{A}\right)$.

Proof. 

By the Cayley–Hamilton theorem,

$${\mathbf{A}}^2=\mathrm{Tr}\left(\mathbf{A}\right)\mathbf{A}-\det \left(\mathbf{A}\right)\mathbf{I}.$$

(2)

Write

$${\mathbf{A}}^n=a_n\mathbf{A}+b_n\mathbf{I},$$

(3)

where $a_1=1$, $b_1=0$, and, by Equation (2),

$$a_2=\mathrm{Tr}\left(\mathbf{A}\right),and{b}_2=-\det \left(\mathbf{A}\right).$$

Assuming that Equation (3) holds for $n-1$, we obtain

$${\mathbf{A}}^n=a_{n-1}{\mathbf{A}}^2+b_{n-1}\mathbf{A}=(\mathrm{Tr}\left(\mathbf{A}\right)a_{n-1}+b_{n-1})\mathbf{A}-\det \left(\mathbf{A}\right)a_{n-1}\mathbf{I}.$$

Hence

$$a_n=\mathrm{Tr}\left(\mathbf{A}\right)a_{n-1}+b_{n-1},\mathrm{and}{b}_n=-\det \left(\mathbf{A}\right)a_{n-1},$$

(4)

and therefore

$$a_n=\mathrm{Tr}\left(\mathbf{A}\right)a_{n-1}-\det \left(\mathbf{A}\right)a_{n-2},n\ge 3.$$

(5)

Now define

$$c_n={\left(\sqrt{\det \left(\mathbf{A}\right)}\right)}^{n-1}{U}_{n-1}\left(x\right).$$

Since $U_0\left(x\right)=1$, $U_1\left(x\right)=2x$, and

$$U_n\left(x\right)=2x{U}_{n-1}\left(x\right)-U_{n-2}\left(x\right),$$

it follows immediately that

$$c_1=1,c_2=\mathrm{Tr}\left(\mathbf{A}\right),c_n=\mathrm{Tr}\left(\mathbf{A}\right)c_{n-1}-\det \left(\mathbf{A}\right)c_{n-2}.$$

Hence $c_n$ and $a_n$ satisfy the same recurrence with the same initial values, so

$$a_n={\left(\sqrt{\det \left(\mathbf{A}\right)}\right)}^{n-1}{U}_{n-1}\left(x\right).$$

Using Equation (4), we derive

$$b_n=-\det \left(\mathbf{A}\right)a_{n-1}=-{\left(\sqrt{\det \left(\mathbf{A}\right)}\right)}^n{U}_{n-2}\left(x\right).$$

Substituting these expressions into Equation (3) yields Equation (1). Finally, using the parity identity

$$U_n(-x)={(-1)}^n{U}_n\left(x\right),$$

we observe that replacing

$$\sqrt{\det \left(\mathbf{A}\right)}⟼-\sqrt{\det \left(\mathbf{A}\right)}$$

changes the factor ${\left(\sqrt{\det \left(\mathbf{A}\right)}\right)}^n$ by ${(-1)}^n$, while the argument of the Chebyshev polynomial changes sign, contributing the same factor. These factors cancel in both terms of Equation (1). Hence the right-hand side is independent of the choice of the square root of $\det \left(\mathbf{A}\right)$, completing the proof. □

Remark 1.

The Chebyshev representation avoids explicit diagonalisation and remains valid even when the eigenvalues are complex. As shown in Theorem 2, taking traces immediately yields a closed-form expression for $\mathrm{Tr}\left({\mathbf{A}}^n\right)$. In transfer-matrix methods of statistical mechanics, where $\mathbf{A}$ is a positive transfer matrix, the resulting trace formula gives exact finite-size expressions for the partition function.

Theorem 2

(Trace of powers of a $2\times 2$ matrix). Assume the hypotheses of Theorem 1. Then, for every integer $n\ge 0$,

$$\mathrm{Tr}\left({\mathbf{A}}^n\right)=2{\left(\sqrt{\det \left(\mathbf{A}\right)}\right)}^n{T}_n\left(x\right),$$

(6)

where $T_n$ denotes the Chebyshev polynomials of the first kind.

Proof. 

Taking the trace of Equation (1) yields

$$\mathrm{Tr}\left({\mathbf{A}}^n\right)=2{\left(\sqrt{\det \left(\mathbf{A}\right)}\right)}^n\left(x{U}_{n-1}\left(x\right)-U_{n-2}\left(x\right)\right).$$

The identity

$$T_n\left(x\right)=x{U}_{n-1}\left(x\right)-U_{n-2}\left(x\right),$$

together with the convention $U_{-1}\left(x\right)=0$, immediately gives Equation (6). □

We conclude this section with several remarks highlighting the broad applicability of the mathematical framework developed above to problems in physics. The same algebraic structure has been extensively employed for the special subgroup $\mathrm{SL}(2,\mathbb{R})\subset \mathrm{GL}(2,\mathbb{R})$, consisting of all $2\times 2$ real matrices with determinant equal to one. This group plays a fundamental role in a wide range of problems in mathematical physics, particularly in the analysis of superlattices, photonic crystals, and more general stratified media.

The use of Chebyshev polynomials in these contexts dates back to the pioneering works of Abelès [6] and Jones [7]. Their subsequent application to transfer-matrix methods and wave-propagation problems has been extensively documented in [8,9,10], as well as in the standard references [11] (§ 1.5.1, p. 69), [12] (§ 1.6.2, p. 55), and the comprehensive reviews [13,14]. In these works, the relevant matrix identities are typically derived using classical techniques, most notably the Cayley–Hamilton theorem, or direct mathematical induction.

In several physical applications, it is advantageous to work with representations of the algebra $\mathrm{SL}(2,\mathbb{C})$, which can be naturally extended to the full matrix algebra $\mathrm{Mat}(2,\mathbb{C})$. The most general formulation of this representation was stated, without proof, in [15]; see also [16], where the associated matrix polynomials are identified with Chebyshev polynomials. These polynomials play a central role in knot theory and arise in a wide variety of mathematical and physical contexts.

The proof presented here applies to an arbitrary invertible matrix $\mathbf{A}\in {\mathbb{R}}^{2\times 2}$ and is both elementary and direct. At the formal level, the representation is identical to that for $\mathrm{SL}(2,\mathbb{C})$; the only distinction is that, when $\det \left(\mathbf{A}\right)<0$, the quantity $\sqrt{\det \left(\mathbf{A}\right)}$ must be interpreted as a square root in $\mathbb{C}$. As established in Theorem 1, the resulting expression is independent of the choice of square root.

3. Ising Model with p Random Variables per Unit Cell

3.1. The Model

Let us consider a chain containing $N\in \mathbb{N}$ sites labelled by $i=1,\dots ,N$. At each site i we define a discrete random variable

$${\sigma }_i\in \{-1,+1\},$$

interpreted in physical sciences as a “spin” variable. The mutual dependence of the random variable is determined by interactions constants $J_i$ and bonds $(i,i+1)$, connecting the random variables ${\sigma }_i$ and ${\sigma }_{i+1}$.

Let us denote the configuration of the system as

$$\sigma =({\sigma }_1,\dots ,{\sigma }_N)\in {\{-1,+1\}}^N.$$

The mutual probability distribution of the N random variables is given by

$${\mathbb{P}}_N\left(\sigma \right|\left\{K_i\right\},h)={\displaystyle \frac{\exp \left[-\beta {\mathcal{H}}_N\left(\sigma \right|\left\{K_i\right\},h)\right]}{Z_N(\left\{K_i\right\},h)}}.$$

(7)

Explicitly, the Hamiltonian of the system is defined by

$${\mathcal{H}}_N\left(\sigma \right|\left\{K_i\right\},h)=-{\displaystyle \sum _{i=1}^N}\left(K_i{\sigma }_i{\sigma }_{i+1}+h{\sigma }_i\right).$$

(8)

depending on the dimensionless parameters

$$K_i=\beta J_i,\mathrm{and}h=\beta H.$$

The normalization factor in Equation (7)

$$Z_N(\left\{K_i\right\},h)=\sum _{\sigma }\exp \left[-\beta {\mathcal{H}}_N\left(\sigma \right|\left\{K_i\right\},h)\right]$$

(9)

is a quantity, used in physical sciences, to define some characteristics of the system of experimental interest, namely—the free energy, specific heat, average value of the random variables (i.e., the magnetisation), and also the Casimir force which is a quantity which study is currently very modern, due to its relation with nanotechnology; definitions of these quantities will be presented below. The parameter $\beta =1/\left(k_{B}T\right)$ has the physical meaning of inverse temperature, where $k_B$ is the Boltzmann constant. We note that the summation in Equation (9) is over all $2^N$ states of the random variables. Different expectation values, e.g., $\langle {\sigma }_i\rangle$ which is called the site magnetisation, are computed with respect to this distribution ${\mathbb{P}}_{\mathbb{N}}\left(\mathbf{\sigma }\right)$, namely

$$\langle {\sigma }_i(\left\{K_i\right\},h)\rangle =\sum _{\mathbf{\sigma }}{\sigma }_i{\mathbb{P}}_{\mathbb{N}}\left(\mathbf{\sigma }\right|\left\{K_i\right\},h).$$

(10)

Obviously, under periodic boundary conditions $\langle {\sigma }_i\rangle$ does not depend of the positions “i” of the random variable. We will call this, as usual in physical sciences, a magnetisation.

From $Z_N(\left\{K_i\right\},h)$ one derives the so-called free energy $f_N(\left\{K_i\right\},h)$ of the finite chain, which is of central importance in thermodynamics and statistical physics

$$\beta f_N(\left\{K_i\right\},h)\equiv -{\displaystyle \frac{1}{N}}\ln Z_N(\left\{K_i\right\},h),$$

(11)

as well as the so-called thermodynamics limit of it, defined as

$$\beta f_b(\left\{K_i\right\},h)\equiv \underset{N\to \infty }{\lim }\beta f_N(\left\{K_i\right\},h),$$

(12)

which is connected to the behaviour of the infinite chains (i.e., of the “bulk system”).

3.2. Transfer Matrix Formalism

The transfer-matrix formalism provides a convenient framework for the analysis of the partition function, allowing the free energy and other thermodynamic quantities to be expressed in terms of the eigenvalues of an appropriate transfer matrix.

Throughout this work, we consider Ising chains with periodically modulated couplings of period p, namely

$$J_{i+p}=J_i.$$

(13)

Such periodically modulated Ising chains (with ferromagnetic coupling constants) were introduced and analysed in [4]. Hence, the chain may be regarded as a periodic repetition of a unit cell of length p.

The corresponding local transfer matrices inherit the same periodicity,

$${\mathbf{T}}_{i+p}={\mathbf{T}}_i,$$

where

$${\mathbf{T}}_i=\left(\begin{array}{cc}{e}^{K_i+h}& e^{-K_i}\\ e^{-K_i}& e^{K_i-h}\end{array}\right),i=1,\dots ,p.$$

(14)

The transfer matrix associated with one period is therefore

$${\mathbf{M}}_p={\mathbf{T}}_1{\mathbf{T}}_2\cdots {\mathbf{T}}_p.$$

(15)

and is commonly referred to as the monodromy matrix, since it describes the evolution of the system over one complete period.

For a chain consisting of n identical unit cells, the total number of sites is $N=np$. The periodic structure then implies that the transfer matrix of the entire chain is simply the n-th power of the monodromy matrix:

$${\mathbf{T}}_1{\mathbf{T}}_2\cdots {\mathbf{T}}_N={\left({\mathbf{M}}_p\right)}^n.$$

(16)

Assuming periodic boundary conditions, the partition function can therefore be written as [4],

$$Z_N=\mathrm{Tr}\left[{\left({\mathbf{M}}_p\right)}^n\right].$$

(17)

Although the monodromy matrix is completely determined by the interaction pattern within a single unit cell, the frustration properties of the corresponding periodic chain depends on the number of repetitions n of that cell. The following proposition provides a precise characterization of the frustration property in the chain.

Proposition 1.

Let

$$\eta ={\displaystyle \prod _{i=1}^p}\mathrm{sgn}\left(J_i\right).$$

Consider a periodic Ising chain consisting of n repetitions of a unit cell of period p, so that $N=np$.Then the sign product around the entire ring is

$${\displaystyle \prod _{i=1}^N}\mathrm{sgn}\left(J_i\right)={\eta }^n.$$

Consequently, the chain is frustrated if and only if

$${\eta }^n=-1,$$

and unfrustrated if and only if

$${\eta }^n=+1.$$

In particular, if

$$\eta =-1,$$

then the chain is unfrustrated precisely when n is even, equivalently,

$$N=2kp,k\in \mathbb{N}.$$

(18)

An immediate consequence of the above proposition is the following corollary, which will be used in Section 8.

Corollary 1.

Consider a periodic Ising chain whose interaction pattern is determined by a repeating unit cell of period p. Assume that the unit cell contains an odd number of antiferromagnetic couplings, i.e.,

$${\displaystyle \prod _{i=1}^p}\mathrm{sgn}\left(J_i\right)=-1.$$

(19)

Let the total number of spins be $N=mp$, so that the chain consists of m repetitions of the unit cell. Then the chain is unfrustrated if and only if m is even. Equivalently,

$$N=2pk,k=1,2,3,\dots .$$

(20)

We now specialize this result to the interaction patterns considered in the present paper. In particular:

1.

For the interaction pattern

$$(J_a,-J_b),$$

which has period $p=2$, the chain is unfrustrated if and only if

$$N=4k.$$

2.

For the interaction pattern

$$(J_a,-J_b,-J_b,-J_b),$$

which has period $p=4$, the chain is unfrustrated if and only if

$$N=8k.$$

The above corollary provides a natural classification of the periodically modulated Ising chains studied in this work into frustrated and unfrustrated classes. Since frustration may substantially modify the transfer-matrix spectrum and, consequently, the thermodynamic behavior of the system, these two classes should generally be treated separately when analyzing finite-size effects and thermodynamic properties.

Although ${\mathbf{M}}_p$ is not unimodular in general, its powers admit an exact representation in terms of Chebyshev polynomials after normalization by $\sqrt{\det \left({\mathbf{M}}_p\right)}$. As a direct consequence of Theorem 2, we obtain for the partition function:

$$\mathrm{Tr}\left({\mathbf{M}}_p^n\right)=2{\left[\sqrt{\det \left({\mathbf{M}}_p\right)}\right]}^n{T}_n\left({\displaystyle \frac{\mathrm{Tr}\left({\mathbf{M}}_p\right)}{2\sqrt{\det \left({\mathbf{M}}_p\right)}}}\right).$$

(21)

In what follows, this representation will be systematically exploited for various values of the parameter p.

4. Free Energy Density for General p-Site Unit Cells

Let ${\lambda }_{p,\pm }$ denote the eigenvalues of the real, nonsingular matrix ${\mathbf{M}}_p$, ordered so that ${\lambda }_{p,+}\ge {\lambda }_{p,-}$. From the characteristic equation of ${\mathbf{M}}_p$, we obtain

$${\lambda }_{p,\pm }={\displaystyle \frac{1}{2}}\left[\mathrm{Tr}\left({\mathbf{M}}_p\right)\pm \sqrt{\mathrm{Tr}{\left({\mathbf{M}}_p\right)}^2-4\det \left({\mathbf{M}}_p\right)}\right].$$

(22)

Proposition 2.

Since every transfer matrix ${\mathbf{T}}_i$ has strictly positive entries, their product ${\mathbf{M}}_p$ also has strictly positive entries. Writing

$${\mathbf{M}}_p=\left(\begin{array}{cc}a& b\\ c& d\end{array}\right),$$

(23)

with $a,b,c,d>0$, we obtain

$$\mathrm{Tr}{\left({\mathbf{M}}_p\right)}^2-4\det \left({\mathbf{M}}_p\right)={(a-d)}^2+4bc>0.$$

(24)

Therefore, the eigenvalues ${\lambda }_{p,\pm }$ are real and distinct, and the corresponding bulk and finite-size free energies are real-valued, as expected.

Making use of Equation (21), one can calculate the bulk and finite-size free energies, employing the following well-known relations:

$$\begin{array}{ccc}\hfill \beta f_b(\left\{K_i\right\},h)& =& -{\displaystyle \frac{1}{p}}\ln \left({\lambda }_{p,+}\right)\hfill \end{array}$$

(25)

$$\begin{array}{ccc}\hfill \beta f_N(\left\{K_i\right\},h)& =& -{\displaystyle \frac{1}{N}}\ln \left[2{\left[\sqrt{\det \left({\mathbf{M}}_p\right)}\right]}^n{T}_n\left({\displaystyle \frac{\mathrm{Tr}\left({\mathbf{M}}_p\right)}{2\sqrt{\det \left({\mathbf{M}}_p\right)}}}\right)\right].\hfill \end{array}$$

(26)

Using the above free energies, one obtains:

• the variance of the Hamiltonian ${\mathcal{H}}_N\left(\sigma \right|\left\{K_i\right\},h)$ in the statistical mechanics is called the specific heat $c_N(\left\{K_i\right\},h)$ of the system, i.e.,

  $$c_N(\left\{K_i\right\},h)/k_B={\displaystyle \frac{1}{N}}{\beta }^2\left[〈{\mathcal{H}}_N{\left(\sigma \right|\left\{K_i\right\},h)}^2〉-{\langle {\mathcal{H}}_N\left(\sigma \right|\left\{K_i\right\},h)\rangle }^2\right]$$

  (27)

  where the averages are taken with respect to the probability distribution Equation (7). It is easy to show that

  $$c_b/k_B={\displaystyle \frac{1}{p}}{\beta }^2{\displaystyle \frac{{\partial }^2}{{\partial }^2\beta }}\ln \left({\lambda }_{p,+}\right)\mathrm{and}{c}_N(\left\{K_i\right\},h)/k_B={\displaystyle \frac{1}{p}}{\beta }^2{\displaystyle \frac{{\partial }^2}{{\partial }^2\beta }}\beta f_N(\left\{K_i\right\},h),$$

  (28)

  respectively for the bulk and finite-size systems. The factor $1/p$ appears because the dominant eigenvalue ${\lambda }_{+}$ corresponds to the transfer matrix of a unit cell containing p spins, rather than to a single spin. In the thermodynamic limit the partition function scales as

  $$Z_N(\left\{K_i\right\},h)\sim {\lambda }_{+}^{N/p},$$

  so that the free energy per spin is obtained by dividing $\ln {\lambda }_{+}$ by p. Consequently, all thermodynamic quantities defined per spin, including the specific heat, acquire the prefactor $1/p$.
• the average values of the random variables, i.e., the magnetisation

  $$m_N(\left\{K_i\right\},h)={\displaystyle \frac{1}{N}}〈{\displaystyle \sum _{i=1}^N}{\sigma }_i〉.$$

  (29)
  It directly follows that

  $$m_b={\displaystyle \frac{1}{p}}\beta {\displaystyle \frac{\partial }{\partial h}}\ln \left({\lambda }_{p,+}\right)and{m}_N(\left\{K_i\right\},h)=-{\displaystyle \frac{1}{p}}{\displaystyle \frac{\partial }{\partial h}}\left[\beta f_N(\left\{K_i\right\},h)\right],$$

  (30)

  respectively for the bulk and finite-size systems.

5. Fluctuation-Induced Casimir Force

The fluctuation-induced critical Casimir force arises from finite-size corrections to the free energy, caused by the restriction of thermal fluctuations. It is therefore closely related to the spectral properties of confined systems and serves as a useful probe of finite-size effects. For general introductions and overviews of the subject, see, e.g., Refs. [17,18].

The Casimir force in the one-dimensional Ising model has been extensively investigated in a number of works [19,20,21,22,23]. In the present work, we extend these studies to Ising chains with alternating interactions. Since alternating couplings lead to a nontrivial modification of the energy spectrum, they are expected to produce distinctive features in the corresponding Casimir force.

By definition, we have

$$\beta F_{\mathrm{Casimir}}(N,\left\{K_i\right\},h)=-\beta {\displaystyle \frac{\partial }{N}}N\left[f_N(\left\{K_i\right\},h)-f_b(\left\{K_i\right\},h)\right].$$

(31)

Let us first show that, in any of the models outlined above, the free energy density of the finite system approaches exponentially, in n, the corresponding one of the infinite system. From Equation (25) we obtain

$$\begin{array}{ccc}\hfill \beta f_N(\left\{K_i\right\},h)& =& -{\displaystyle \frac{1}{p}}\log {\lambda }_{p,+}-{\displaystyle \frac{1}{N}}\log \left[1+{\left({\displaystyle \frac{{\lambda }_{p,-}(\left\{K_i\right\},h)}{{\lambda }_{p,+}(\left\{K_i\right\},h)}}\right)}^n\right],\hfill \\ & =& \beta f_b(\left\{K_i\right\},h)-{\displaystyle \frac{1}{N}}\log \left[1+\exp \left(-{\displaystyle \frac{n}{\xi (\left\{K_i\right\},h)}}\right)\right],\hfill \end{array}$$

(32)

where

$$\xi (\left\{K_i\right\},h)\equiv \log \left[{\displaystyle \frac{{\lambda }_{p,+}(\left\{K_i\right\},h)}{{\lambda }_{p,-}(\left\{K_i\right\},h)}}\right].$$

(33)

Since the quantities defined above—the specific heat, the magnetisation and the Casimir force—are expressed via derivates of $\beta f_N(\left\{K_i\right\},h)$ it follows that they approach their bulk (infinite chain) value exponentially fast when $\xi =\mathcal{O}\left(1\right)$. For the Casimir force we obtain that it itself vanishes exponentially. In statistical physics the quantity $\xi$ is termed the correlation length and describes the rates at which the correlation between the random variables decays with the distance between them.

We close that general chapter by elucidating that the Casimir force is negative, i.e., attractive in physical sciences terminology, independent on the value of p. For establishing that it is enough to perform the explicit calculation of the force using Equation (32). Indeed, one has

$$\beta F_{\mathrm{Casimir}}(N,\left\{K_i\right\},h)=-{\displaystyle \frac{1}{N}}{\displaystyle \frac{n}{\xi (\left\{K_i\right\},h)}}{\displaystyle \frac{\exp \left(-\frac{n}{\xi (\left\{K_i\right\},h)}\right)}{1+\exp \left(-\frac{n}{\xi (\left\{K_i\right\},h)}\right)}}<0.$$

(34)

The above property holds independent on p and on the choice of values, either positive or negative, of coupling constants $K_1,\cdots ,K_p$. A visualisation of $X_{\mathrm{Casimir}}\equiv N\beta F_{\mathrm{Casimir}}(N,\left\{K_i\right\},h)$ as a function of the scaling variable $x=n/\xi (\left\{K_i\right\},h)$ is presented in Figure 1.

6. Heat Capacity and Schottky Anomaly

The Schottky anomaly is a physical phenomenon characterized by the occurrence of one or more maxima in the temperature dependence of the specific heat. Mathematically, it means that the variance of the sum of the random variables, defining the Hamiltonian of the system, possesses maximum, or maxima, as a function of the parameters characterising the joint probability distribution in the system. It originates from the structure of the energy spectrum and therefore provides useful information about the excitation energies of the system. Since the specific heat vanishes in both the low- and high-temperature limits, it must attain at least one maximum at a finite temperature. For general discussions of the Schottky anomaly, see, e.g., pp. 54–56 of Ref. [25], as well as Refs. [26,27,28].

From a physical perspective, a Schottky peak appears when the thermal energy becomes comparable to the energy gap separating the ground state from a set of excited states. In this temperature range, the occupation probabilities of the corresponding energy levels change most rapidly, leading to enhanced energy fluctuations and, consequently, to a maximum in the specific heat. If several well-separated energy scales are present, the system may exhibit multiple Schottky peaks, each associated with a different excitation gap.

In the present study, the occurrence of a two-peak Schottky anomaly in the one-dimensional Ising model with alternating ferromagnetic and antiferromagnetic interactions ($p=2$) can be understood directly from the spectral decomposition of the transfer matrix. Since the bulk free energy is determined by the logarithm of the dominant eigenvalue, the temperature dependence of thermodynamic quantities is governed by the structure of the transfer-matrix spectrum. When the leading eigenvalue factorizes or contains well-separated characteristic scales, the free energy decomposes into additive contributions associated with distinct excitation gaps. As a consequence, the specific heat exhibits Schottky-type maxima at temperatures determined by these gaps. This spectral interpretation will be made explicit in Sec. Section 8.4.

In what follows we derive exact specific expressions for several Ising type models with p random variables in their unit cell.

7. Ferromagnetic or Antiferromagnetic One-Dimensional Ising Models ($p=1$)

7.1. On the Case of Negative Determinant $\det \mathbf{T}<0$

We begin by considering the standard one-dimensional Ising model under periodic boundary conditions with interaction coupling K and external magnetic field h. We consider both the case ($K>0$), when the mutual influence of two adjacent random variables is positive, i.e., increasing the likelihood the two to have the same values, or the case ($K<0$), when the opposite is true. In physical sciences the chain with ($K>0$) is called a ferromagnetic one, while ($K<0$) is termed antiferromagnetic.

Figure 2. Topology of the standard Ising model with identical spins A under periodic boundary conditions. Red (blue) circles denote spins in the state $+1$ ($-1$) for the representative configuration shown.

Figure 2. Topology of the standard Ising model with identical spins A under periodic boundary conditions. Red (blue) circles denote spins in the state $+1$ ($-1$) for the representative configuration shown.

These models correspond to the special case $p=1$, for which $N=np$ reduces to $N=n$; see Equation (16). In this case, ${\mathbf{M}}_1={\mathbf{T}}_1\equiv \mathbf{T}\left(h\right)$, with $K_1=K$, is the standard transfer matrix of the one-dimensional Ising model. It is a real symmetric matrix with strictly positive entries. Consequently, its eigenvalues ${\lambda }_{1,\pm }\in \mathbb{R}$, with ${\lambda }_{1,+}>{\lambda }_{1,-}$. The eigenvalues are given by

$${\lambda }_{1,\pm }={\displaystyle \frac{1}{2}}\left[\mathrm{Tr}\left({\mathbf{M}}_1\right)\pm \sqrt{\mathrm{Tr}{\left({\mathbf{M}}_1\right)}^2-4\det \left({\mathbf{M}}_1\right)}\right].$$

(35)

Hence, Equation (21), together with

$$\mathrm{Tr}\left(\mathbf{T}\left(h\right)\right)=2e^K\cosh \left(h\right),\det \left(\mathbf{T}\left(h\right)\right)=2\sinh \left(2K\right),$$

(36)

gives

$$\mathrm{Tr}\left(\mathbf{T}{\left(h\right)}^N\right)=2{\left[\sqrt{2\sinh \left(2K\right)}\right]}^N{T}_N\left({\displaystyle \frac{e^K\cosh \left(h\right)}{\sqrt{2\sinh \left(2K\right)}}}\right),K\in \mathbb{R}\setminus \left\{0\right\}.$$

(37)

Since the determinant of the transfer matrix is negative whenever $K<0$, at this point the formula Equation (37) is an extension of our previous consideration in refs. [20,29].

Although $2\sinh \left(2K\right)<0$ for $K<0$, Remark 1 guarantees that the expression in Equation (37) is real and positive and hence

$$\ln Z_N\left(h\right)\in \mathbb{R}.$$

(38)

Thus, the partition function remains real and positive even for $K<0$, despite the fact that the Chebyshev representation formally involves the square root of the negative quantity $2\sinh \left(2K\right)$.

Equations (25) and (28) can be specialized to the present case by setting p = 1, which yields explicit expressions for the free-energy density and the specific heat. As these expressions are straightforward, we do not display them here. Instead, we focus on the corresponding numerical results, which are sufficient to highlight the main physical features of the system, in particular the Schottky anomaly.

7.1.1. Exact Symmetry with Respect to $K=0$ and $h=0$

Figure 3 shows that, for a periodic chain with an even number of sites N and in the absence of an external magnetic field, the specific heat and the associated Schottky anomaly are identical for ferromagnetic ($J>0$) and antiferromagnetic ($J<0$) interactions, depending only on $\left|J\right|$. In contrast, this symmetry is lost for a periodic chain with an odd number of sites N, as illustrated in Figure 3. This symmetry follows from the staggered spin transformation

$${\sigma }_i⟶{(-1)}^i{\sigma }_i,$$

(39)

which is well defined on a bipartite lattice, such as a one-dimensional periodic chain with an even number of sites, and maps the Hamiltonian with coupling J onto the equivalent Hamiltonian with coupling $-J$. Indeed,

$${\sigma }_i{\sigma }_{i+1}⟶{(-1)}^i{(-1)}^{i+1}{\sigma }_i{\sigma }_{i+1}=-{\sigma }_i{\sigma }_{i+1},$$

so that each nearest-neighbour interaction term changes sign.

Although this transformation acts locally on every bond, its global consistency under periodic boundary conditions, ${\sigma }_{N+1}\equiv {\sigma }_1$, requires the lattice to be bipartite, which is the case only for even N. For odd N, the lattice is not bipartite, and the staggered transformation Equation (39) cannot be defined consistently around the ring because the bond connecting the last and first sites joins two sites of the same parity. Consequently, the Hamiltonian with coupling J is no longer mapped onto that with coupling $-J$. The antiferromagnetic model is therefore geometrically frustrated, and its partition function and thermodynamic properties generally differ from those of the corresponding ferromagnetic model. Hence, for finite periodic chains with an odd number of sites, the free energy, specific heat, and other thermodynamic quantities are not functions of $\left|J\right|$ alone.

It is worth noting that the finite-size specific heat in the low-temperature regime, particularly under periodic boundary conditions, was investigated in [30]. The results presented here are fully consistent with those reported therein, thereby providing an independent confirmation of their findings.

7.1.2. Finite-Size Effects

The stronger finite-size effects observed for $J<0$ in Figure 4 originate from the competition between the antiferromagnetic exchange interaction and the external magnetic field. This competition generates low-lying excited states that remain close in energy to the ground state, thereby enhancing the contribution of the subleading eigenvalue of the transfer matrix. As a result, finite-size corrections decay slowly with the system size N, leading to a pronounced size dependence of the specific heat.

By contrast, for $J>0$ the exchange interaction and the magnetic field act cooperatively, resulting in a more robust energy gap and a faster suppression of subleading eigenvalue contributions. Consequently, finite-size corrections decay rapidly with increasing N, and the specific heat becomes only weakly dependent on the system size. This behaviour is reflected in Figure 4, where the finite-size curves quickly collapse onto the bulk result.

In summary, the finite-size behaviour is governed by whether the magnetic field competes with or reinforces the exchange interaction, which determines the effective spectral gap of the transfer matrix and controls the rate of convergence to the thermodynamic limit.

8. Ising Chain of N Sites with Period $p=2$ ($AB$) Interaction Structure

The alternating Ising model with with ferromagnetic coupling constants was first introduced in Ref. [31] and later analyzed in Refs. [32,33,34]. As a prototypical one-dimensional Ising system with spatially modulated exchange interactions, it has served as an important framework for investigating the effects of inhomogeneity and aperiodicity on magnetic properties. In this work we consider the simplest realization of this class, corresponding to spatial period $p=2$, in our notation. The couplings satisfy Equation (13), thereby defining a two-site unit cell. Each cell contains dynamical variables A and B, with nearest-neighbour interactions $K_a$ and $K_b$ assigned to the bonds $A-B$ and $B-A$, respectively (see Figure 5).

Figure 5. Alternating Ising model with two spin species A and B arranged periodically as $ABAB\dots$. The nearest-neighbour couplings $K_a$ and $K_b$ act on the bonds $A-B$ and $B-A$, respectively.

Figure 5. Alternating Ising model with two spin species A and B arranged periodically as $ABAB\dots$. The nearest-neighbour couplings $K_a$ and $K_b$ act on the bonds $A-B$ and $B-A$, respectively.

In the transfer-matrix formalism, we associate the matrices ${\mathbf{T}}_a$ and ${\mathbf{T}}_b$ with the couplings $K_a$ and $K_b$, respectively and corresponding monodromy matrix:

$${\mathbf{M}}_2:={\mathbf{T}}_a{\mathbf{T}}_b.$$

(40)

For a ring composed of n unit cells, the partition function (see Equation (17)) is

$$Z_N=\mathrm{Tr}\left[{\left({\mathbf{T}}_a{\mathbf{T}}_b\right)}^n\right]=\mathrm{Tr}\left({\mathbf{M}}_2^n\right).$$

(41)

Periodic boundary conditions require $N=2n$, ensuring an integer number of unit cells and compatibility between lattice periodicity and boundary closure, thus avoiding boundary-induced frustration.

8.1. Thermodynamic Properties

The thermodynamic properties of the system are determined by the eigenvalues ${\lambda }_{\pm }$ of ${\mathbf{M}}_2$. The matrix explicitly reads

$${\mathbf{M}}_2=\left(\begin{array}{cc}{e}^{K_a+K_b+2h}+e^{-K_a-K_b}& e^{K_a-K_b+h}+e^{-K_a+K_b-h}\\ e^{-K_a+K_b+h}+e^{K_a-K_b-h}& e^{-K_a-K_b}+e^{K_a+K_b-2h}\end{array}\right).$$

(42)

For the above $2\times 2$ matrix, the eigenvalues are

$${\lambda }_{2,\pm }={\displaystyle \frac{\mathrm{Tr}\left({\mathbf{M}}_2\right)}{2}}\pm {\displaystyle \frac{1}{2}}\sqrt{{\left[\mathrm{Tr}\left({\mathbf{M}}_2\right)\right]}^2-4\det \left({\mathbf{M}}_2\right)}.$$

(43)

In the present case,

$$\mathrm{Tr}\left({\mathbf{M}}_2\right)=2\left(e^{K_a+K_b}\cosh \left(2h\right)+e^{-(K_a+K_b)}\right),$$

(44)

and

$$\det \left({\mathbf{M}}_2\right)=\det \left({\mathbf{T}}_a\right)\det \left({\mathbf{T}}_b\right)=4\sinh \left(2K_a\right)\sinh \left(2K_b\right).$$

(45)

Inserting Equations (44) and (45) into Equation (43) yields

$$\begin{array}{ccc}\hfill {\lambda }_{2,\pm }& =& e^{K_a+K_b}\cosh \left(2h\right)+e^{-(K_a+K_b)}\hfill \\ & \pm & \sqrt{{\left(e^{K_a+K_b}\cosh \left(2h\right)+e^{-(K_a+K_b)}\right)}^2-4\sinh \left(2K_a\right)\sinh \left(2K_b\right)}.\hfill \end{array}$$

(46)

By a straightforward algebraic transformation, the discriminant in Equation (46) can be rewritten as

$$\begin{array}{cc}\hfill \Delta & ={\left(e^{K_a+K_b}\cosh \left(2h\right)+e^{-(K_a+K_b)}\right)}^2-4\sinh \left(2K_a\right)\sinh \left(2K_b\right)\hfill \\ & =e^{2(K_a+K_b)}{\sinh }^2\left(2h\right)+2\left(\cosh \left(2h\right)+\cosh \left(2(K_a-K_b)\right)\right),\hfill \end{array}$$

making it immediately clear that $\Delta >0$, so that ${\lambda }_{2,\pm }\in \mathbb{R}$. This also follows from Proposition 2, but the derivation is instructive.

Equations (44) and (45) generalize the results of [32,33], which were obtained in the special ferromagnetic case $K_a,K_b>0$. The explicit form of the determinant makes the dependence on the signs of $K_a$ and $K_b$ transparent, as it enters through the factors $\sinh \left(2K_a\right)$ and $\sinh \left(2K_b\right)$. In the following, we allow for the more general case $K_a,K_b\in \mathbb{R}$.

Next, it is easy to show that, as follows from Equation (25)

$$\beta f_b(\left\{K_i\right\},h)=-{\displaystyle \frac{1}{2}}\ln \left({\lambda }_{2,+}\right),$$

(47)

and

$$\beta f_N(\left\{K_i\right\},h)=-{\displaystyle \frac{1}{N}}\ln \left[2{\left[\sqrt{\det \left({\mathbf{M}}_2\right)}\right]}^n{T}_n\left({\displaystyle \frac{\mathrm{Tr}\left({\mathbf{M}}_2\right)}{2\sqrt{\det \left({\mathbf{M}}_2\right)}}}\right)\right].$$

(48)

The dependence on the signs of $K_a$ and $K_b$ is explicit through $\sinh \left(2K_a\right)\sinh \left(2K_b\right)$. Since $\sinh \left(x\right)$ is an odd function,

$$\sinh (-2K_{a,b})=-\sinh \left(2K_{a,b}\right),$$

the determinant changes sign whenever one coupling changes sign. Therefore:

$$K_a{K}_b>0⟹\det \left({\mathbf{T}}_{ab}\right)>0,$$

while

$$K_a{K}_b<0⟹\det \left({\mathbf{T}}_{ab}\right)<0.$$

Thus, for nonzero magnetic field $h\ne 0$, the free energy generally depends on the signs of the couplings, distinguishing ferromagnetic, antiferromagnetic, and mixed ferro–antiferromagnetic configurations.

8.2. Zero-Field Limit

For $h=0$, Equation (46) yields, after straightforward algebra, the largest eigenvalue

$${\lambda }_{2,+}=4\cosh \left(K_a\right)\cosh \left(K_b\right).$$

(49)

The corresponding zero-field bulk free energy is

$$\beta f_b(\left\{K_i\right\},0)=-{\displaystyle \frac{1}{2}}\ln \left({\lambda }_{2,+}\right).$$

Since $\cosh \left(x\right)$ is an even function, $\beta f_b(\left\{K_i\right\},0)$ depends only on the absolute values $|K_a|$ and $|K_b|$. In particular, all four combinations $(K_a,K_b)$, $(-K_a,K_b)$, $(K_a,-K_b)$, and $(-K_a,-K_b)$ share the same zero-field bulk free energy and the zero-field bulk thermodynamics is invariant under sign reversals of the couplings.

8.3. Specific Heath

The specific heat per spin as defined an Equation (28) the per-spin normalization yields

$$c_b={\displaystyle \frac{k_B{\beta }^2}{2}}{\displaystyle \frac{{\partial }^2}{\partial {\beta }^2}}\ln \left({\lambda }_{2,+}\right).$$

(50)

The temperature dependence of the bulk specific heat $c_b$, as obtained from Equation (50), is illustrated in Figure 6.

The specific heat thus appears as a superposition of two Schottky-type terms associated with the characteristic energy scales $J_a$ and $J_b$. When these scales are well separated, $|J_a|\ll |J_b|$, the corresponding activation temperatures $T_i\sim J_i/k_B$ are parametrically distinct, leading to two resolved maxima (see: brown, green and black dashed curves in Figure 6) As the interaction strengths become more similar (blue and red curves Figure 6), these maxima progressively coalesce into a single broadened peak.

The double-peak structure therefore directly reflects the coexistence of two distinct energy scales in the transfer-matrix spectrum.

The present model thus provides a minimal realization of a general mechanism: the appearance of multiple Schottky-type maxima as a thermodynamic manifestation of well-separated characteristic energy scales. Similar double-peak structures arise in a variety of systems with hierarchical excitation spectra (see Ref. [26] and references therein).

8.4. On the Emergence of the Double Schottky Structure

The bulk free energy per spin is determined by the dominant eigenvalue of the two-site transfer matrix. For $h=0$, the largest eigenvalue factorizes, see Equation (49).

Consequently, the bulk free energy becomes additive in the two coupling constants, and the specific heat assumes the exact form

$$c_b\left(T\right)={\displaystyle \frac{1}{2k_B{T}^2}}\left[J_a^2{\mathrm{sech}}^2\left(\beta J_a\right)+J_b^2{\mathrm{sech}}^2\left(\beta J_b\right)\right].$$

(51)

The specific heat is therefore a superposition of two Schottky-type contributions associated with the energy scales $|J_a|$ and $|J_b|$. Each contribution exhibits a maximum at $k_{B}T\simeq 0.83\left|J_i\right|$ ($i=a,b$). When the two energy scales are sufficiently different, the corresponding maxima become well separated, resulting in a characteristic double-peak (double-Schottky) structure. In the special case $|J_a|=|J_b|$, the two maxima coincide and merge into a single Schottky peak.

8.5. Frustration and Finite-Size Effects

Finite-size effects and frustration in a periodic chain must be analyzed simultaneously, as they are intrinsically intertwined. The finite system size together with the periodic boundary conditions restricts the admissible spin configurations and may modify the manifestation of frustration, while frustration itself can have a significant impact on the finite-size behavior of thermodynamic observables. Consequently, a consistent analysis must treat both effects on an equal footing. Using Equations (28) and (48), this interplay is illustrated in Figure 7, where the system sizes have been chosen in accordance with Corollary 1.

In the thermodynamic limit, the effects of frustration associated with the finite system size vanish, and both the frustrated and unfrustrated finite-size curves converge to the same bulk behavior, as required. The distinction between the two cases is manifested only in the direction of convergence: as shown in Figure 7, one sequence of finite-size curves approaches the bulk curve from above, whereas the other approaches it from below.

9. Summary

The analysis presented in this work relies on Theorems 1 and 2, which express powers and traces of ($\mathrm{GL}(2,\mathbb{R})$) matrices in terms of Chebyshev polynomials. These results furnish a convenient analytical framework for the study of one-dimensional Ising chains with periodically modulated nearest-neighbour interactions, which constitute the main subject of the present work. The class of models considered here is parametrized by an integer (p), representing the number of interactions within the monodromy matrix and thereby fixing the period of the interaction pattern along the chain. For a chain composed of (n) identical unit cells, the total number of sites is ($N=np$), and the transfer matrix of the entire system is given by the (n)-th power of the monodromy matrix.

Particular attention has been devoted to the representative cases $(p=1)$ and $(p=2)$. While $(p=1)$ corresponds to the standard ferromagnetic or antiferromagnetic Ising chain, the case $(p=2)$, associated with the interaction pattern $(J_a,-J_b$), involve competing ferromagnetic and antiferromagnetic couplings within a single period. Mathematically, if the the interaction is positive (i.e., for a ferromagnetic chain) the adjacent random variables tend to have equal signs (values), while the opposite is true if it is negative (i.e., for an antiferromagnetic chain).

We have performed a detailed study of the double-peak Schottky anomaly in the specific heat. The Schottky anomaly is called an anomaly because the specific heat curve shows a non-intuitive, unexpected peak, or peaks, that does not come from collective critical behaviour, but instead from a very simple microscopic mechanism: a two-level system with an energy gap, or gaps. It looks like a “critical peak”, but it is not reflecting a phase transition(s). That mismatch between appearance and physical origin is why it is historically called an anomaly. Mathematically, it means that the variance of the sum of the random variables, defining the Hamiltonian of the system, possesses maximum, or maxima, as a function of the parameters characterising the joint probability distribution in the system.

Our results confirm that the appearance of two specific-heat maxima is a universal consequence of the presence of two well-separated characteristic energy scales, differing by roughly an order of magnitude, in agreement with previous findings [26]. Unlike earlier investigations, which were largely based on phenomenological multi-level systems, the present work demonstrates in an exact manner the emergence of this phenomenon in a one-dimensional Ising model with competing ferromagnetic and antiferromagnetic interactions subject to an external magnetic field.

One of results of this work is the classification of the models into frustrated and unfrustrated subclasses, according to the parity of the number $\left(n\right)$ of unit-cell repetitions. The chain is unfrustrated if and only if $N=0\mathrm{mod}\left(2p\right)$ and frustrated, otherwise. Consequently, both the $p=1$ and $p=2$ models exhibit frustrated, as well as unfrustrated regimes. They are reflected by nontrivial thermodynamic behaviour and pronounced finite-size effects. We visualised that on the example of the variance of the Hamiltonian of the system, termed specific heat in statistical mechanics, in Figure 3, for $p=1$, and Figure 7, for $p=2$. We observe, that the specific heat, for the considered values of the parameters defining the joint probability distribution, possess two maxima as a function of T for $p=2$. For $p=1$ the specific heat shows one maximum when the sign of interaction is fixed - one for positive, and one for a negative one. The role of the magnetic field in the $p=1$ case is shown in Figure 4. There we see that the finite-size effects depend strongly on the cooperation ($J>0,h>0$), or, in contrast, in competition ($J<0,h>0$) between the exchange interaction and the magnetic field.

Finally, let us note that within our approach in terns out to be straightforward to determine the Casimir force within the considered models. Its scaling function is shown in Figure 1. It is remarkable, that this behaviour holds for any set of variables $\left\{K_i\right\}$ and the magnetisation field h. Note that the behaviour is the same for any p, which a valuable demonstration of the universality [17,24] hypothesis in the theory of critical Casimir force. We note that under periodic boundary conditions the force, as expected, is negative, having the physical meaning of an attraction.

To conclude, let us note that along the lines outlined in the current article, one can consider other, more complicated models, e.e., the case where a ferromagnetic bond of strength $J_a$ is followed by three antiferromagnetic bonds of strength $-J_b$. Repetition of this bond sequence generates a lattice with unit-cell period $p=4$. A related model was previously investigated in Ref. [35]. That study demonstrated that such a bond arrangement produces a characteristic multi-peak structure in the specific heat as a function of temperature and external magnetic field, as well as a terrace-like structure of the magnetisation. We hope to return to study of such, more complicated systems, in the near future.