# Unconventional low temperature features in the one-dimensional frustrated
$q$-state Potts model

Yury Panov  Ural Federal University, 19 Mira street, 620002 Ekaterinburg,
Russia  Onofre Rojas  Department of Physics, Federal University of Lavras,
37200-900 Lavras-MG, Brazil

###### Abstract

Here we consider a one-dimensional $q$-state Potts model with an external
magnetic field and an anisotropic interaction that selects neighboring sites
that are in the spin state 1. The present model exhibits an unusual behavior
in the low-temperature region, where we observe an anomalous vigorous change
in the entropy for a given temperature. There is a steep behavior at a given
temperature in entropy as a function of temperature, quite similar to first-
order discontinuity, but there is no jump in the entropy. Similarly, second
derivative quantities like specific heat and magnetic susceptibility also
exhibit a strong acute peak rather similar to second-order phase transition
divergence, but once again there is no singularity at this point. Correlation
length also confirms this anomalous behavior at the same given temperature,
showing a strong and sharp peak which easily one may confuse with a
divergence. The temperature where occurs this anomalous feature we call
pseudo-critical temperature. We have analyzed physical quantities, like
correlation length, entropy, magnetization, specific heat, magnetic
susceptibility, and distant pair correlation functions. Furthermore, we
analyze the pseudo-critical exponent that satisfy a class of universality
previously identified in the literature for other one-dimensional models,
these pseudo-critical exponents are: for correlation length $\nu=1$, specific
heat $\alpha=3$ and magnetic susceptibility $\mu=3$.

##  I Introduction

The advantage of exactly solvable models is its easy handling to analyze
several properties, which can show interesting features despite it is
simplicity. In contrast, more detailed models are rarely exactly solvable, and
it would restrict us to performing only through numerical computations, which
prevents further analysis of these types of models. Some one-dimensional
models[1] can help us understand and predict leading behavior in more complex
models. From the experimental side, one-dimensional models accurately describe
several chemical compounds[2, 3]. That is why the one-dimensional models are
quite important to investigate, both from theoretical and experimental points
of view.

Earlier in the fifties, van Hove[4] proposed a theorem to verify the absence
of phase transition in uniform one-dimensional models with short-range
interaction. The validity of the proposed theorem follows the conditions: (i)
homogeneity, (ii) the Hamiltonian should not include particles positions terms
(like external fields). (iii) hard-core particles. Based on the Perron-
Frobenius theorem [5], condition, the van Hove theorem[4] is restricted to
limited one-dimensional systems. Later, Cuesta and Sanchez[6] tried to extend
the theorem of non-existence phase transition for a more general one-
dimensional system. Mainly, they included an external field and considered
point-like particles, which extends the theorem. Even with this extension, it
is still far from being a fully general theorem of non-existence phase
transition.

There are unusual one-dimensional models with a short-range coupling that
exhibit a phase transition at finite temperature. The zipper or Kittel model
[7], which is one of the simplest models with a finite size transfer matrix
that exhibits a first-order phase transition. Another model, considered by
Chui-Wicks[8], is typical of models called solid-on-solid for surface growth.
It has the infinite dimension transfer-matrix and is exactly solvable. Because
of the impenetrable condition of the subtract, the model shows the existence
of a finite temperature phase transition. One more model is that considered by
Dauxois-Peyrard[9], with an infinite dimension transfer matrix, which can be
explored numerically. Lately, Sarkanych et al.[10] proposed a one-dimensional
Potts model with invisible states and short-range coupling. The term invisible
means an additional energy degeneracy, which only contributes to the entropy,
but not the interaction energy. They named these states the invisible states,
which generate the first-order phase transition.

Motivated by low-dimensional systems, such as the simple zipper model[7] that
describes the long-chain nucleotides of deoxyribonucleic acid (DNA), Zimm and
Bragg[11] introduced an essentially phenomenological cooperative parameter,
which provides narrow helix-coil transitions. Since that several
investigations were driven in the literature[12, 13, 14, 15]. Cooperative
systems in one-dimension can be well represented by Potts-like models[12, 14],
where the helix-coil transition in polypeptides[13] can be studied, which is a
typical application of theoretical physics to macromolecular systems, the
results of which are quite appropriate to understand the physical properties
of the helix-coil transition. The polycyclic aromatic surface elements of the
carbon nanotube (CNT) and the aromatic DNA provide reversible adsorption.
Tonoyan et al.[15] adapted the Hamiltonian of the zipper model in order to
take into account the DNA-CNT interactions.

Potts model is a generalization of the Ising model to more than two
components, such as interacting spins in a crystalline lattice. Standard
$q$-state Potts model [16], with $q\geqslant 2$ has been assumed as an integer
denoting the number of states of each site. Potts model is quite relevant in
statistical physics. In some crystal-lattices would occur vacancies, which
leads to a site diluted Potts model. The site dilute $q$-state Potts model[17,
18] is equivalent to $q+1$-state standard Potts model[16]. Chaves and
Riera[19] investigated a particular case of dilute Potts chain. Recently, a
different dilute Ising spin-1 chain[20] was also studied in the framework of
the projection operator.

An unusual property called pseudo-transition was observed in some recent
works: Like in a double-tetrahedral chain of localized Ising spins with mobile
electrons showing a strong thermal excitation that easily suggests the
existence of a first-order phase transition[21, 22, *galisova20]. Similarly,
the frustrated spin-1/2, Ising-Heisenberg’s three-leg tube exhibited a pseudo-
transition[24]. In reference [25], this property was also observed when
studying the specific heat, where reported a sharp peak on the spin-1/2 Ising-
Heisenberg ladder with alternating Ising and Heisenberg inter-leg couplings.
This weird property was even observed in the spin-1/2 Ising diamond chain in
the neighboring of the pseudo-transition[26]. Besides, deeper investigations
were performed on this peculiar property in [27]. Additionally, the distant
correlation functions have been studied around pseudo-transition for a
spin-1/2 Ising-XYZ diamond chain[28]. A bit different proposal to identify the
pseudo-transition in one-dimensional models[29, *Rojas2020] was explored in
the framework of the phase boundary residual entropy relationship with the
finite temperature pseudo-transition. Further investigation around the pseudo-
transition was also focused on the universality and pseudo-critical exponents
of one-dimensional models[31].

We organize the article as follows. In Section II we present the model and
analyze the zero-temperature phase diagrams. In Section III, we study the
thermodynamics of the model and explore an anomalous phenomenon called pseudo-
transition for several physical quantities. In Section IV, we investigate the
manifest of the pseudo-transition in terms of the distant pair correlation
functions and correlation length. The pseudo critical exponents for the
correlation length, specific heat, and susceptibility we discuss in Section V.
Finally, in Section VI, we present our conclusions.

##  II The frustrated Potts chain

Let us consider a Potts model[16], here we assume the one-dimensional case,
whose Hamiltonian becomes

$H=-\sum_{i=1}^{N}\left\\{J\delta_{\sigma_{i},\sigma_{i+1}}+K\delta_{\sigma_{i},1}\delta_{1,\sigma_{i+1}}+h\delta_{\sigma_{i},1}\right\\},$
(1)

with $\sigma=\left\\{1,\ldots q\right\\}$. Whereas $J$ is bound coupling
parameter, $h$ denotes the external field aligned to state 1, and $K$ denotes
the parameter of an anisotropic interaction that selects neighboring sites
that are in the spin state 1.

The Hamiltonian (1) for the case $q=2$, like standard Potts model drops into
spin-1/2 Ising chain model. For the case of $q\geqslant 3$, the model becomes
a frustrated system for certain choices of parameters, as shown below.

It is worth notice that the Hamiltonian (1) can also be equivalent to a
diluted Potts model with $q-1$ states, as demonstrated by Wu [16]. The
features of the critical properties of the two-dimensional diluted Potts model
were studied earlier using Fortuin-Kasteleyn clusters [32] and the transfer
matrix method [33].

###  II.1 Zero temperature phase diagram

In order to analyze the phase diagram of $q$-state Potts model at zero
temperature, we identify four ground states assuming $q\geqslant 3$ for the
model (1), which read below

$\displaystyle|FM_{1}\rangle=$
$\displaystyle\prod_{i=1}^{N}\bigl{(}|1\rangle_{i}\bigr{)},$ (2)
$\displaystyle|FM_{2}\rangle=$
$\displaystyle\prod_{i=1}^{N}\bigl{(}|j\rangle_{i}\bigr{)},\quad
j=\\{2,\dots,q\\},$ (3) $\displaystyle|FR_{1}\rangle=$
$\displaystyle\prod_{i=1}^{N/2}\left(|1\rangle_{2i-1}|\sigma_{i}\rangle_{2i}\right),$
(4) $\displaystyle|FR_{2}\rangle=$
$\displaystyle\prod_{i=1}^{N}\left(|\sigma_{i}\rangle\right),\quad\sigma_{i}\neq\sigma_{i\pm
1},$ (5)

where $\sigma_{i}=\\{2,\dots,q\\}$. Additional information, concerning phases
and phase-boundaries, are listed in Table 1 for different physical quantities
at zero temperature, like magnetization, entropy, and pair distribution
functions (PDF), which can be obtained by taking the limit of $T\rightarrow 0$
of the quantities (23), (55) and (67).

In Fig. 1a is illustrated schematically the ground state phase diagram in the
plane $J-h$, for the Hamiltonian (1) assuming $K<0$. Here we observe four
phases illustrated above. Similarly, in panel (b) and (c) are illustrated
alternative phase diagrams for $K=0$ and $K>0$, respectively.

Figure 1: The ground state phase diagrams of the $q$-state Potts chain with
Hamiltonian (1) for $q\geqslant 3$. (a) $K<0$, (b) $K=0$, and (c) $K>0$. Here
the points are: $Q_{1}(0,0)$, $Q_{2}(-K/2,-K)$, $S(0,0)$, $P_{1}(-K,0)$,
$P_{2}(0,-K)$.

The first state is a type of ferromagnetic ($FM_{1}$) phase where all sites
are in state $\sigma_{i}=1$, with energy per spin $E_{FM_{1}}=-(J+K+h)$. The
second state is another type of "ferromagnetic" ($FM_{2})$ phase with energy
$E_{FM_{2}}=-J$. The average fraction of pairs of adjacent spins in the same
state $\mu$, where $\mu=2,\ldots q$, is PDF $g_{\mu,\mu}^{(1)}=1/(q-1)$ (see
Table 1). At the same time, $g_{\mu,\mu^{\prime}}^{(1)}=0$, where
$\mu\neq\mu^{\prime}$, which means that in the thermodynamic limit, the
fraction of pairs with different spin states and their contribution to the
energy of the system are equal to zero. This implies that, in the general
case, the $FM_{2}$ phase consists of $q-1$ sorts of equivalent macroscopic
ferromagnetic domains with spins in the $\mu$ state. Equation (3) corresponds
to the single-domain case. The entropy of both $FM_{1}$ and $FM_{2}$ phases is
zero. The third is a type frustrated ($FR_{1}$) phase, with alternating sites
are in states $\sigma_{2i-1}=1$, while the other sites $\sigma_{2i}$ can take
$\sigma_{2i}=2,\ldots q$, and the corresponding energy is $E_{FR_{1}}=-h/2$.
This phase state is frustrated at $q>2$. Due to every second site can be in
any of the $q-1$ states, the entropy per spin is equal to
$\mathcal{S}=\frac{1}{2}\ln(q-1)$. The fourth state is another type of
frustrated ($FR_{2}$) phase, in each site $\sigma_{i}$ can take independently
$\sigma_{i}=\\{2,\dots,q\\}$ but $\sigma_{i}\neq\sigma_{i\pm 1}$, whose
corresponding energy is $E_{FR_{2}}=0$. In this case, each subsequent site of
the chain can be in one of $q-2$ states, and the entropy is equal to
$\mathcal{S}=\ln(q-2)$. If $q>3$, the $FR_{2}$ phase is evidently frustrated
state. Although, for $q=3$, there is an alternation of sites in states
$\sigma_{i}=2$ and $\sigma_{i}=3$, and the entropy of the $FR_{2}$ phase is
zero.

It is important to note that two different situations can occur at phase
boundaries. The first case is when the states of two adjacent phases are mixed
at the microscopic level. For example, for the $FR_{1}$-$FM_{1}$ boundary, the
$FR_{1}$ state of any pair of sites can be changed to the $FM_{1}$ state, and
vice versa. Such a replacement does not lead to the appearance of microscopic
states from other phases, and the energy of the system does not change. As a
result, the entropy of such a mixed state at the phase boundary is greater
than the entropy of the adjacent phases. A similar situation is observed for
the boundaries $FR_{1}$-$FR_{2}$, $FR_{2}$-$FM_{2}$, and $FR_{1}$-$FM_{2}$ at
$K<0$. In the second case, it can be a pure state of one of the adjacent
phases, or the phase separation, when each phase is represented by the
macroscopic domains. So, on the $FR_{2}$-$FM_{1}$ boundary for $K>0$, the
replacement for a pair of neighboring nodes in state $FR_{2}$ by state
$FM_{1}$ leads to the appearance of $FR_{1}$ states, which are energetically
unfavorable. A similar situation occurs at the $FM_{1}$-$FM_{2}$ boundary.
Note that in the $FM_{1}$-$FM_{2}$ interface curve the residual entropy is
zero.

##  III Thermodynamics

The frustrated $q$-state Potts model Hamiltonian (1) can be solved through
transfer matrix technique, which results in a $q$-dimensional matrix, given by

$\mathbf{V}=\begin{pmatrix}d_{1}&t_{1}&t_{1}&\cdots&t_{1}&t_{1}\\\
t_{1}&d_{2}&t_{2}&\cdots&t_{2}&t_{2}\\\
t_{1}&t_{2}&d_{2}&\cdots&t_{2}&t_{2}\\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
t_{1}&t_{2}&t_{2}&\cdots&d_{2}&t_{2}\\\
t_{1}&t_{2}&t_{2}&\cdots&t_{2}&d_{2}\end{pmatrix},$ (6)

where $d_{1}=xkz$, $d_{2}=x$, $t_{1}=\sqrt{z}$, $t_{2}=1$, and
$x=\mathrm{e}^{\beta J}$, $k=\mathrm{e}^{\beta K}$, $z=\mathrm{e}^{\beta h}$.

Let us write the transfer matrix eigenvalues similarly to that defined in
reference [27], that being so, we have

$\displaystyle\lambda_{1}=$
$\displaystyle\,\frac{1}{2}\left(w_{1}+w_{-1}+\sqrt{(w_{1}-w_{-1})^{2}+4w_{0}^{2}}\right),$
(7) $\displaystyle\lambda_{2}=$
$\displaystyle\,\frac{1}{2}\left(w_{1}+w_{-1}-\sqrt{(w_{1}-w_{-1})^{2}+4w_{0}^{2}}\right),$
(8) $\displaystyle\lambda_{j}=$ $\displaystyle\,d_{2}-t_{2},\quad
j=\\{3,4,\dots,q\\},$ (9)

where

$\displaystyle w_{1}$ $\displaystyle=$ $\displaystyle d_{1},$ (10)
$\displaystyle w_{-1}$ $\displaystyle=$ $\displaystyle
d_{2}+\left(q-2\right)t_{2},$ (11) $\displaystyle w_{0}$ $\displaystyle=$
$\displaystyle\sqrt{q-1}\,t_{1}.$ (12)

The corresponding transfer matrix eigenvectors are

$\displaystyle|u_{1}\rangle=$
$\displaystyle\cos(\phi)|1\rangle+\tfrac{\sin(\phi)}{\sqrt{q-1}}\sum_{\mu=2}^{q}|\mu\rangle,$
(13) $\displaystyle|u_{2}\rangle=$
$\displaystyle-\sin(\phi)|1\rangle+\tfrac{\cos(\phi)}{\sqrt{q-1}}\sum_{\mu=2}^{q}|\mu\rangle,$
(14) $\displaystyle|u_{j}\rangle=$
$\displaystyle\sqrt{\tfrac{j-2}{j-1}}\Bigl{(}\tfrac{1}{j-2}\sum_{\mu=2}^{j-1}|\mu\rangle-|j\rangle\Bigr{)},\quad
j=\\{3,\cdots,q\\},$ (15)

where $\phi=\frac{1}{2}\cot^{-1}\left(\frac{w_{1}-w_{-1}}{2w_{0}}\right)$,
with $-\frac{\pi}{4}\leqslant\phi\leqslant\frac{\pi}{4}$.

By using the transfer matrix eigenvalues, we express the partition function as
follows

$\displaystyle Z_{N}=$
$\displaystyle\,\lambda_{1}^{N}+\lambda_{2}^{N}+(q-2)\lambda_{3}^{N}$
$\displaystyle=$
$\displaystyle\,\lambda_{1}^{N}\left\\{1+\Bigl{(}\tfrac{\lambda_{2}}{\lambda_{1}}\Bigr{)}^{N}+(q-2)\Bigl{(}\tfrac{\lambda_{3}}{\lambda_{1}}\Bigr{)}^{N}\right\\}.$
(16)

It is evident that the eigenvalues satisfy the following relation
$\lambda_{1}>\lambda_{2}\geqslant\lambda_{3}$. Hence, assuming $q$ finite, the
free energy per spin in thermodynamic limit reduces to

$f=-T\ln\left(\lambda_{1}\right).$ (17)

Note that the free energy is a continuous function with no singularity or
discontinuity, thus we do not expect any real phase transition at finite
temperature.

###  III.1 Pseudo-Critical temperature

Recently pseudo-critical temperature has been discussed in Ising and Ising-
Heisenberg spin models[21, 24, 25, 27], in several one-dimensional spin
models.

To find pseudo-critical temperature, we follow the same strategy to that used
in reference[27]. In our case the largest eigenvalues has the same structure
to that found in reference[27], so necessary conditions for a pseudo-
transition are met, if $w_{1}\sim w_{-1}\ggg w_{0}$, $|w_{1}-w_{-1}|\gg
w_{0}$. The pseudo-transition point can be obtained when the first term inside
the square root of $\lambda_{1}$ given by eq.(7) turns to zero, which gives

$e^{\frac{J+K+h}{T_{p}}}=q-2+e^{\frac{J}{T_{p}}}.$ (18)

In principle, using the above relation, one can find the critical temperature
as a function of some Hamiltonian parameters.

The $q$-state Potts chain does not exhibit a real spontaneous long range order
at finite temperature since its one-dimensional character. Therefore we define
a term “quasi” to refer low temperature regions mainly dominated by ground
state configuration. Hence $FR_{2}$ in low temperature region is called as
$qFR_{2}$, and so on. As shown in [29, *Rojas2020], pseudo-transitions occur
for states near those phase boundaries whose residual entropy is a continuous
function of the model parameters for at least one of the adjacent phases. As
discussed earlier, the state of the $FR_{2}-FM_{1}$ boundary coincides with
the $FR_{2}$ state, so for the $qFR_{2}-qFM_{1}$ boundary, we get from eq.(18)
the following relation

${\rm e}^{\frac{K+h}{T_{p}}}={\rm e}^{-\frac{J}{T_{p}}}(q-2)+1,$ (19)

which we can simplify and write approximately in the form

$T_{p}=\frac{J+K+h}{\ln\left(q-2\right)}.$ (20)

This is the known expression [21, 24, 27, 29, *Rojas2020] for the pseudo-
transition temperature:

$E_{FM_{1}}-E_{FR_{2}}=T_{p}\left(\mathcal{S}_{FR_{2}}-\mathcal{S}_{FM_{1}}\right),$
(21)

where the energy and entropy per unit cell are given at zero temperature.
Since $E_{FM_{1}}=E_{FR_{2}}$ at the $qFR_{2}$-$qFM_{1}$ boundary, $T_{p}$
tends to zero near to it.

Another phase boundary we focus is $qFM_{2}$-$qFM_{1}$. It is worth to mention
that, the entropy of the $FM_{1}$ and $FM_{2}$ phases is zero, so the entropy
is a continuous function for both adjacent phases. For the $qFM_{2}$-$qFM_{1}$
boundary the eq.(18) can be approximately written in the following form

$\frac{1}{T_{p}}\left(K+h\right)e^{J/T_{p}}=q-2.$ (22)

Figure 2: Density plot: (a) for entropy in the plane $T-h$, for $K=1$ and
$J=-0.5$. (b) Magnetization $m_{1}$ for the same set of parameters to the
panel (a). In (c) entropy in the plane $T-h$ for fixed $K=1$, $J=0.99$. In (d)
magnetization for the same set of parameter to the panel (c).

In Fig. 2a is reported the density plot of entropy in the plane $T-h$,
assuming fixed parameters $K=1$, $J=-0.5$. Dashed curve describes the boundary
$qFR_{2}$-$qFM_{1}$, which corresponds to the pseudo-critical temperature
$T_{p}$ as a function of $h$, according eq.(19). It can be seen that the curve
is an almost straight line well represented by (20). We can observe also how
the sharp boundary between quasi-phase melts smoothly for higher temperature.
Similar density plot is depicted in panel(b) for the magnetization $m_{1}$ in
the plane $T-h$ for the same set of parameters in panel (a). Analogously, we
analyze the phase boundary between $qFM_{2}-qFM_{1}$ in panel (c), assuming
fixed parameters $K=1$, $J=0.99$ and $T=0.01$. The dashed line is given by
eq.(19) and nicely approximated by (22), since there is no residual entropy in
the boundary the quasi-phase $qFM_{1}$ and $qFM_{2}$ leads to zero when
temperature vanishes, by looking entropy we cannot distinguishes the boundary
of quasi-phase. However, in panel (d) we illustrate the density plot of
magnetization $m_{1}$, for the same set of parameters to the panel (c), and we
observe clearly a sharp boundary between $qFM_{1}$ and $qFM_{2}$ regions, and
this boundary melts smoothly as soon as temperature increases.

###  III.2 Entropy and Specific heat

The entropy and the specific heat of the system can be obtained from the free
energy (17) by

$\mathcal{S}=-\frac{\partial f}{\partial T},\qquad
C=-T\frac{\partial^{2}f}{\partial T^{2}}.$ (23)

Of particular interest is the behavior of thermodynamic characteristics near
the pseudo-transition point. The general method for considering this issue was
developed in Ref. [31]. For the one-dimensional Potts model, it is possible to
find an explicit form of approximation of the free energy and other
thermodynamic quantities near $T_{p}$. Assuming that $\tau=(T-T_{p})/T_{p}\ll
1$ and taking into account the equation (18), we can write

$\displaystyle w_{1}=$
$\displaystyle\;e^{\beta(J+K+h)}=\tilde{w}_{1}^{\frac{1}{1+\tau}}\approx\tilde{w}_{1}\left(1-\ln\tilde{w}_{1}\,\tau\right),$
(24) $\displaystyle w_{-1}=$ $\displaystyle\;q-2+e^{\beta
J}=q-2+x_{p}^{\frac{1}{1+\tau}}\approx\tilde{w}_{1}-x_{p}\ln x_{p}\,\tau,$
(25)

where

$\tilde{w}_{1}=\left.w_{1}\right|_{\tau=0}=q-2+x_{p},$ (26)

and we define the parameters $x_{p}=e^{J/T_{p}}$ and $k_{p}=e^{K/T_{p}}$. The
value of $w_{1}-w_{-1}$ is zero at $\tau=0$ due to equation (18), so we get

$w_{1}-w_{-1}\approx-\tilde{w}_{1}\,a\,\tau,$ (27)

where we introduce the parameter $a$ as

$a=\ln\tilde{w}_{1}-\frac{x_{p}\ln x_{p}}{\tilde{w}_{1}}.$ (28)

Note that the parameter $a^{-1}$ define the slope of the pseudo-transition
curve in the plane $T-h$ (see Fig. 2a,c), $dT_{p}/dh=a^{-1}$.

The condition $|w_{1}-w_{-1}|\gg w_{0}$ which is met in the vicinity of
$T_{p}$ causes quasi-singular behaviour of the derivative of square root in
equations (7,8) at $\tau=0$. Assuming this is the case, we can write the
approximation

$4w_{0}^{2}=4(q-1)e^{\beta h}\approx 4(q-1)\frac{\tilde{w}_{1}}{x_{p}k_{p}}.$
(29)

This allows us to write eigenvalues (7,8) near $\tau=0$ as

$\lambda_{1,2}\approx\tilde{w}_{1}\left[1+\left(\tfrac{1}{2}a\tau-\ln\tilde{w}_{1}\right)\tau\pm\tfrac{1}{2}\sqrt{a^{2}\tau^{2}+b^{2}}\right],$
(30)

where $b^{2}=4w_{0}^{2}/\tilde{w}_{1}^{2}$, and yields the approximation for
the free energy (17) in the vicinity of $T_{p}$:

$f\approx-
T_{p}\left[\ln\tilde{w}_{1}-\tfrac{1}{2}a\tau+\tfrac{1}{2}\sqrt{a^{2}\tau^{2}+b^{2}}\right].$
(31)

The expressions (30,31) being exact at $\tau=0$ have a small deviation from
rigorous expansions of $\lambda_{1,2}$ and $f$ in the tiny vicinity
$0<|\tau|<\ln(x_{p}k_{p})\,b^{2}/a^{2}\lll 1$ due to the neglect of linear
terms in equation (29), but well approximate the functions $\lambda_{1,2}$ and
$f$ at $\ln(x_{p}k_{p})\,b^{2}/a^{2}<|\tau|\ll 1/|\ln x_{p}|$. This implies a
simple necessary condition for a pseudo-transition in the form $b/a\ll
T_{p}/|J|$.

For the entropy in the same region of $\tau$, using (23), we obtain the
following expression:

$\mathcal{S}\approx\frac{a}{2}\left(1+\frac{a\tau}{\sqrt{a^{2}\tau^{2}+b^{2}}}\right),$
(32)

The equation (32) describes the entropy jump in a small vicinity of $T_{p}$

$\Delta\mathcal{S}_{p}=\mathcal{S}\left(\tau>\tfrac{b}{a}\right)-\mathcal{S}\left(\tau<-\tfrac{b}{a}\right)=a,$
(33)

which may be related to the “latent heat” of pseudo-transition,
$Q=T_{p}\Delta\mathcal{S}_{p}=aT_{p}$.

Second derivation of equation (31) by temperature gives the approximation of
the specific heat near $T_{p}$,

$C\approx\frac{a^{2}b^{2}}{2\left(a^{2}\tau^{2}+b^{2}\right)^{\frac{3}{2}}}.$
(34)

This allows us to estimate the maximum value of the specific heat in $T_{p}$
as

$C_{p}=\tfrac{1}{2}a^{2}b^{-1}.$ (35)

We can qualify the peak of the specific heat near $T_{p}$ by its half-width at
half-maximum $\Psi_{\tau}$. From (34) we find that $\Psi_{\tau}=\gamma b/a$,
where $\gamma=\sqrt{2^{2/3}-1}\approx 0.7664$, and hence $\Psi_{\tau}\ll 1$
due to a necessary condition for a pseudo-transition.

Figure 3: The entropy and specific heat of the $qFM_{1}$ states near the
$FR_{2}$-$FM_{1}$ boundary (a, b) with $q=5$, $J=-0.5$, $K=1$, and (c, d) with
$q=3$, $J=-0.05$, $K=1$, in an external field $h$. The dotted lines in (b) and
(d) show the magnitude of $C_{p}$ at $T=T_{p}$ given by equations (37) and
(39) respectively.

The entropy and specific heat of the $qFM_{1}$ states having $q=5$, $J=-0.5$,
$K=1$ in a given external field $h$ are shown in Fig. 3a,b. At sufficiently
low temperatures, the $qFR_{2}$-$qFM_{1}$ pseudo-transition is observed, which
is accompanied by a jump in entropy and a narrow peak in the specific heat.

For the $qFR_{2}$-$qFM_{1}$ transition at $q>3$, the condition $x_{p}\ll q-2$
is met, since $J<0$ and $|J|/T_{p}\gg 1$, so approximately

$a=\ln(q-2),\qquad b=2\sqrt{\frac{q-1}{(q-2)x_{p}k_{p}}}.$ (36)

The entropy jump $\Delta\mathcal{S}_{p}=\ln(q-2)$ equals to the residual
entropy of the $qFR_{2}$ phase (see Table 1). An expression for the maximum
value of the specific heat,

$C_{p}=\frac{\ln^{2}(q-2)}{4}\sqrt{\frac{q-2}{q-1}}\,e^{\frac{K-|J|}{2T_{p}}},$
(37)

shows that $C_{p}$ decreases for the states close to the point $P_{1}$ in Fig.
1c. The magnitudes of the specific heat peaks (37) are shown in Fig. 3b with
dotted lines. It is interesting to note, that
$C_{p}\Psi_{\tau}=\tfrac{1}{2}\gamma\ln(q-2)={\rm const}(T_{p})$, so, having a
finite height, the specific heat peak tends to the delta function near the
$FR_{2}$-$FM_{1}$ boundary with the $T_{p}$ lowering.

In Fig. 3c,d the entropy and specific heat of the $qFM_{1}$ states having
$q=3$, $J=-0.05$, $K=1$ in a given external field $h$ are shown. The case
$q=3$ is special, since the residual entropy of the $FR_{2}$-$FM_{1}$ boundary
is a continuous function for both adjacent phases. For the parameter $a$ we
obtain

$a=\left(1-\ln x_{p}\right)x_{p},$ (38)

so the necessary condition for a pseudo-transition $b/a\propto
e^{(3|J|-K)/2T_{p}}\ll T_{p}/|J|$ is met if $K>3|J|$. Since in this case
$J<0$, the entropy jump (33) drops with decreasing of $T_{p}$ due to the
exponent in $x_{p}=e^{J/T_{p}}$, as can be seen from Fig. 3c. A maximal value
of the specific heat in $T_{p}$ approximately has the following form

$C_{p}=\frac{\left(1-J/T_{p}\right)^{2}}{4\sqrt{2}}\,e^{\frac{K-5|J|}{2T_{p}}}.$
(39)

From (39) we may conclude that if $K>5|J|$ the pseudo-transition is
accompanied with exponentially high peak of the specific heat. From a general
point of view, this issue was discussed in detail in [29, *Rojas2020].

Figure 4: The entropy and specific heat of the $qFM_{1}$ states near the
$FM_{1}$-$FM_{2}$ boundary (a, b) with $q=5$, $J=0.25$, $K=1$ in an external
field $h$, and (c, d) with $q=5$, $h=-0.995$, $K=1$ for the different values
of the coupling parameter $J$. The dotted lines in (b) and (d) show the
magnitude of $C_{p}$ at $T=T_{p}$ given by equation (41).

Figure 4 shows the temperature dependencies of the entropy and specific heat
for states having parameters close to the $FM_{1}$-$FM_{2}$ boundary for
$K=1$. One set of states, shown in Fig. 4a,b, has $J=0.25$ and different
values of the external field $h$, and the states in another set, shown in Fig.
4c,d, have the same $h$ and differ in the coupling constant $J$.

Near the $FM_{1}$-$FM_{2}$ boundary, for which the residual entropy is zero,
as for both adjacent phases, we get

$a=(q-2)\frac{1+\ln x_{p}}{x_{p}},\qquad
b=\frac{2}{x_{p}}\sqrt{\frac{q-1}{k_{p}}}.$ (40)

The condition for a pseudo-transition $b/a\propto e^{-K/2T_{p}}\ll T_{p}/|J|$
will met only if $K>0$. In this case $J>0$, so the entropy jump
$\Delta\mathcal{S}_{p}=a\propto e^{-J/T_{p}}$ decreases exponentially with
decreasing $T_{p}$. This effect is shown in Fig. 4a. The specific heat in
$T_{p}$ approximately is given by

$C_{p}=\frac{(q-2)^{2}\left(1+J/T_{p}\right)^{2}}{4\sqrt{q-1}}\,e^{\frac{K-2J}{2T_{p}}}.$
(41)

Equation (41) shows that for $K>2J$ the $qFM_{2}$-$qFM_{1}$ pseudo-transition
entails an exponentially high peak of the specific heat if the pseudo-
transition temperature is small enough. The decrease in the peak of the
specific heat near $T_{p}$ with increasing $J$ is shown in Fig. 4d.

###  III.3 Magnetization and Magnetic susceptibility

It is important to study the magnetization property of the present model. In
order to obtain a general average $\langle\delta_{\sigma,\mu}\rangle$, let us
define the following operator

$\bm{m}_{\mu}=\sum_{j=1}^{q}\delta_{j,\mu}|j\rangle\langle
j|=|\mu\rangle\langle\mu|.$ (42)

Therefore, the average of $\bm{m}_{\mu}$ becomes

$\langle\delta_{\sigma,\mu}\rangle=m_{\mu}=Z_{N}^{-1}\mathop{\mathrm{tr}}\left(\bm{m}_{\mu}\mathbf{V}^{N}\right).$
(43)

Expressing (42), in eigenstate basis given in (13-15), we have

$\tilde{\bm{m}}_{\mu}=\sum_{j_{1},j_{2}=1}^{q}\langle
u_{j_{1}}|\bm{m}_{\mu}|u_{j_{2}}\rangle|u_{j_{1}}\rangle\langle u_{j_{2}}|,$
(44)

with

$\langle u_{j_{1}}|\bm{m}_{\mu}|u_{j_{2}}\rangle=\langle
u_{j_{1}}|\mu\rangle\langle\mu|u_{j_{2}}\rangle,$ (45)

and denoting $\langle\mu|u_{j}\rangle=c_{\mu,j}$, we can write (44) as

$\tilde{\bm{m}}_{\mu}=\sum_{j_{1},j_{2}=1}^{q}c_{\mu,j_{1}}^{*}c_{\mu,j_{2}}|u_{j_{1}}\rangle\langle
u_{j_{2}}|.$ (46)

In a similar way, the transfer matrix in eigenstate basis (13-15), becomes

$\tilde{\mathbf{V}}=\sum_{j=1}^{q}\lambda_{j}|u_{j}\rangle\langle u_{j}|.$
(47)

Hence, substituting (46) and (47) in eq.(43), we obtain

$m_{\mu}=\frac{1}{Z_{N}}\\!\sum_{j_{1}=1}^{q}\langle
u_{j_{1}}|\negmedspace\sum_{j_{2},j_{3}=1}^{q}\negmedspace
c_{\mu,j_{2}}^{*}c_{\mu,j_{3}}|u_{j_{2}}\rangle\langle
u_{j_{3}}|\lambda_{j_{3}}^{N}|u_{j_{1}}\rangle.$ (48)

In thermodynamic limit ($N\rightarrow\infty$) we can simplify (48). Therefore,
writing the magnetization $m_{\mu}$ in terms of $c_{\mu,j}$, it reduces to

$m_{\mu}=c_{\mu,1}^{*}c_{\mu,1}=|c_{\mu,1}|^{2}.$ (49)

The explicit expression of coefficients are given by

$\displaystyle c_{1,1}=\cos\phi,\quad$ $\displaystyle
c_{\mu,1}=\frac{\sin\phi}{\sqrt{q-1}},$ (50) $\displaystyle
c_{1,2}=\sin\phi,\quad$ $\displaystyle c_{\mu,2}=\frac{\cos\phi}{\sqrt{q-1}},$
(51)

where $\mu=\\{2,\dots,q\\}$.

Analogously the remaining coefficients are written by

$\displaystyle c_{\mu,j}=$
$\displaystyle\begin{cases}0;&\mu=\\{1,j+1,\dots,q\\}\\\
\frac{1}{\sqrt{(j-1)(j-2)}};&\mu=\\{2,\dots,j-1\\}\\\
\sqrt{\frac{j-1}{j-2}};&\mu=j\end{cases},$ (52)

where we consider $j=\\{3,\cdots,q\\}$.

As a consequence, we can get the magnetization by using (50) in (49), which
becomes

$\displaystyle m_{1}=\cos^{2}\phi,\quad$ $\displaystyle
m_{\mu}=\frac{\sin^{2}\phi}{q-1}.$ (53)

Additionally, from above result, we get the following identity

$m_{1}+(q-1)m_{\mu}=1,\quad\mu=\\{2,\dots,q\\}.$ (54)

Alternatively one can obtain $m_{1}$ taking the derivative of free energy with
respect to external field $h$,

$\displaystyle m_{1}$ $\displaystyle=$
$\displaystyle\left\langle\delta_{\sigma,1}\right\rangle=-\frac{\partial
f}{\partial h},$ (55)

and $m_{\mu}$ we can obtain from (54).

On the other hand, the magnetic susceptibility $\chi_{1}$ can be obtained
deriving (55), which results in

$\chi_{1}=\frac{1}{4T}\sin(2\phi)^{3}\left(\frac{d_{1}+d_{2}+(q-2)t_{2}}{2t_{1}\sqrt{q-1}}\right),$
(56)

Similarly, one can obtain for $\mu=\\{2,\dots,q\\}$,

$\chi_{\mu}=\frac{\partial m_{\mu}}{\partial h},$ (57)

and we have the following relation

$\chi_{\mu}=-\frac{1}{q-2}\chi_{1},$ (58)

for $\mu=\\{2,\dots,q\\}$.

To find an approximate expression for the magnetization near $T_{p}$, we can
write it by using equations (17) and (55) in the form

$m_{1}=\frac{T}{\lambda_{1}}\frac{\partial\lambda_{1}}{\partial h}.$ (59)

When calculating the derivative with respect to $h$, we take into account that
$\tfrac{\partial}{\partial h}w_{1}=\beta w_{1}$, $\tfrac{\partial}{\partial
h}w_{-1}=0$, $\tfrac{\partial}{\partial h}w_{0}^{2}=\beta w_{0}^{2}$, and
obtain

$m_{1}=\frac{1}{2\lambda_{1}}\left(w_{1}+\frac{(w_{1}-w_{-1})w_{1}+2w_{0}^{2}}{(w_{1}-w_{-1})^{2}+4w_{0}^{2}}\right).$
(60)

Using equations (24-29) and leaving only the leading terms, we find an
approximation for $m_{1}$ in the following form:

$m_{1}\approx\frac{1}{2}\left(1-\frac{a\tau}{\sqrt{a^{2}\tau^{2}+b^{2}}}\right).$
(61)

Equation (61) describes the jump in magnetization from $m_{1}=0$ at $\tau>b/a$
to $m_{1}=1$ at $\tau<-b/a$ in the small vicinity of $T_{p}$ for both
$qFR_{2}$-$qFM_{1}$ and $qFM_{2}$-$qFM_{1}$ pseudo-transitions.

Similarly, for the susceptibility $\chi_{1}$ we may write, differentiating
(59) with respect to $h$:

$\chi_{1}=-\frac{m_{1}^{2}}{T}+\frac{T}{\lambda_{1}}\frac{\partial^{2}\lambda_{1}}{\partial
h^{2}}.$ (62)

The quasi-singular behavior in $T_{p}$ is caused by the second term in (62).
Using the same steps as in deriving (61) and leaving the main contributions,
we come to the following approximations for the susceptibility near $T_{p}$:

$\chi_{1}\approx\frac{b^{2}}{2T_{p}\left(a^{2}\tau^{2}+b^{2}\right)^{\frac{3}{2}}},$
(63)

and its maximum value:

$\chi_{1,p}=\frac{b^{-1}}{2T_{p}}.$ (64)

Figure 5: The specific heat (blue line) and susceptibility (orange line) peaks
in the vicinity of pseudo-transition points at (a) $q=3$, $J=-0.25$, $K=1$,
$h=-0.749995$, and (b) $q=5$, $J=1$, $K=1$, $h=-0.9999$.

Comparing (34) and (63), we see that dependencies of the specific heat and
susceptibility in the vicinity of $T_{p}$ are similar: they have the same
half-width at half-maximum $\Psi_{\tau}$ and differ by a scale factor
$\alpha$, which is the ratio of the specific heat and susceptibility in the
pseudo-transition point

$\alpha=\frac{C_{p}}{\chi_{1,p}}=a^{2}T_{p}.$ (65)

If for the $qFR_{2}$-$qFM_{1}$ pseudo-transition at $q>3$ the factor $\alpha$
decreases linearly with lowering $T_{p}$, for the $qFR_{2}$-$qFM_{1}$ pseudo-
transition at $q=3$ and the $qFM_{2}$-$qFM_{1}$ pseudo-transition, it
decreases exponentially due to dependence of $a$. In both the latter cases,
for some parameters, the giant magnitude of $\chi_{1,p}$ can be realized at
the low magnitude of $C_{p}$. Indeed, for the $qFR_{2}$-$qFM_{1}$ pseudo-
transition at $q=3$, we have $\chi_{1,p}\propto e^{(K-|J|)/2T_{p}}$, so the
giant peak of susceptibility exists at sufficiently low $T_{p}$ in the entire
pseudo-transition region defined, as it was found earlier (see equation (38)),
by the condition $K>3|J|$, but for the specific heat it becomes exponentially
high only if $K>5|J|$ due to equation (39). This case is shown in Fig. 5a,
where in the vicinity of the $qFR_{2}$-$qFM_{1}$ pseudo-transition at $q=3$,
$J=-0.25$, $K=1$, the peaks of the specific heat and susceptibility differ in
amplitude by seven orders of magnitude and coincide in shape. In turn, for the
$qFM_{2}$-$qFM_{1}$ pseudo-transition we have $\chi_{1,p}\propto
e^{(K+2J)/2T_{p}}$. Taking into account (40), we can conclude that the giant
peak of susceptibility exists for all $J>0$ at $K>0$, while for the specific
heat, the giant peak exists only at $0<2J<K$, as it follows from equation
(41). The similarity of the specific heat and susceptibility peaks in the
vicinity of the $qFM_{2}$-$qFM_{1}$ pseudo-transition at $q=5$, $J=1$, $K=1$,
is shown in Fig. 5b.

Temperature dependencies of the magnetic moment and susceptibility of the
$qFM_{1}$ states near the $FR_{2}$-$FM_{1}$ boundary are shown in Fig. 6.
These states both for $q=5$, $J=-0.5$ and $q=3$, $J=-0.05$ exhibit the
$qFR_{2}$-$qFM_{1}$ pseudo-transition with the continuous jump in the magnetic
moment and the exponentially high peak of susceptibility. Figure 7 shows the
magnetic moment and susceptibility for the same set of states as in Fig. 4.
The dotted lines in Fig. 6b,d and in Fig. 7b,d show an approximate value of
susceptibility in the pseudo-transition point defined by equation (64).

Figure 6: The magnetic moment and susceptibility of the $qFM_{1}$ states near
the $FR_{2}$-$FM_{1}$ boundary (a, b) with $q=5$, $J=-0.5$, $K=1$, and (c, d)
with $q=3$, $J=-0.05$, $K=1$, in an external field $h$. The dotted lines in
(b) and (d) show the magnitude of $\chi_{1,p}$ at $T=T_{p}$ given by equation
(64).

Figure 7: The magnetic moment and susceptibility of the $qFM_{1}$ states near
the $FM_{1}$-$FM_{2}$ boundary (a, b) with $q=5$, $J=0.25$, $K=1$ in an
external field $h$, and (c, d) with $q=5$, $h=-0.995$, $K=1$ for the different
values of the coupling parameter $J$. The dotted lines in (b) and (d) show the
magnitude of $\chi_{1,p}$ at $T=T_{p}$ given by equation (64).

Temperature dependencies of the magnetic moment, and susceptibility for states
near the $FM_{1}$-$FM_{2}$ boundary at $K=0$ and $K<0$ are shown in Fig. 8. It
can be seen the qualitative difference from the case $K>0$. At $K\leq 0$, the
magnetic moment changes with decreasing temperature from the value at the
$FM_{1}$-$FM_{2}$ boundary, $m_{1}=\tfrac{1}{q-1}$ at $K=0$ or
$m_{1}=\tfrac{1}{2}$ at $K<0$, to the values of the magnetic moment in the
$FM_{1}$ or $FM_{2}$ phase, which are equal to $1$ and $0$. Because of this,
the susceptibility peak is observed in both the $qFM_{1}$ and $qFM_{2}$
states.

Figure 8: The magnetic moment and susceptibility near the $FM_{1}$-$FM_{2}$
boundary for $qFM_{1}$ and $qFM_{2}$ states (green and blue lines) with
$q=10$, $J=1$, (a, b) for $K=0$, in an external field equal to $h=\Delta h$,
and (c, d) for $K=-1$, in an external field equal to $h=1+\Delta h$.

As shown if Fig. 8b,d, the peak of susceptibility at $K\leq 0$ becomes
arbitrarily large in height for the $qFM_{1}$ and $qFM_{2}$ states, which are
quite close to the $FM_{1}$-$FM_{2}$ boundary, but it is significantly wider
compared to the case of a pseudo-transition at $K>0$, or near the
$qFR_{2}$-$qFM_{1}$ pseudo-transition. The extremely narrow peak, which is
characteristic of a pseudo-transition, exists only for $K>0$. This follows
formally from equation (40) and the necessary inequality $b/a\ll T_{p}/|J|$.
The physical reason for this is the phase separation in the $FM_{1}$-$FM_{2}$
boundary at $K\leq 0$, which causes an intermediate value of the magnetic
moment at the $FM_{1}$-$FM_{2}$ boundary and its jump in both the $qFM_{1}$
and $qFM_{2}$ states. The one-dimensional ferromagnetic Ising model has the
same properties: zero magnetization of the ground state in the absence of an
external field is achieved by splitting into ferromagnetic domains with
opposite magnetization, and this state represents the boundary on the phase
diagram between phases with magnetizations equal $+1$ and $-1$ depending on
the sign of the external field. This situation can be considered as
intermediate between microscopic mixing of neighboring phases at the phase
boundary when there is no pseudo-transition, and a pure phase equal to one of
the neighboring phases when the pseudo-transition is realized.

##  IV Pair distribution correlation function

In order to accomplish our analysis, we consider the pair distribution
correlation function, which is defined as follows

$\Gamma_{\mu,\mu^{\prime}}^{(r)}=\langle\delta_{\sigma_{i},\mu}\delta_{\sigma_{i+r},\mu^{\prime}}\rangle-\langle\delta_{\sigma_{i},\mu}\rangle\langle\delta_{\sigma_{i+r},\mu^{\prime}}\rangle,$
(66)

where
$\Gamma_{\mu,\mu^{\prime}}^{(r)}\equiv\Gamma_{\mu,\mu^{\prime}}(\sigma_{i},\sigma_{i+r})$.

If the system is translationally invariant, we have
$\langle\delta_{\sigma_{i},\mu}\rangle=\langle\delta_{\sigma_{i+r},\mu}\rangle=\langle\delta_{\sigma,\mu}\rangle=m_{\mu}$,
and $\Gamma_{\mu,\mu^{\prime}}^{(r)}$ depends only on the distance $(r)$.

Now let us defined the average of PDF,

$g_{\mu,\mu^{\prime}}^{(r)}=\langle\delta_{\sigma_{i},\mu}\delta_{\sigma_{i+r},\mu^{\prime}}\rangle=\langle\bm{m}_{\mu}\bm{m}_{\mu^{\prime}}\rangle,$
(67)

so, in natural basis is expressed by

$g_{\mu,\mu^{\prime}}^{(r)}=Z_{N}^{-1}\mathop{\mathrm{tr}}\left(\bm{m}_{\mu}\mathbf{V}^{r}\bm{m}_{\mu^{\prime}}\mathbf{V}^{N-r}\right)$
(68)

whereas in eigenstate basis becomes

$g_{\mu,\mu^{\prime}}^{(r)}=Z_{N}^{-1}\mathop{\mathrm{tr}}\left(\tilde{\bm{m}}_{\nu}\mathbf{\tilde{V}}^{r}\tilde{\bm{m}}_{\tau}\tilde{\mathbf{V}}^{N-r}\right).$
(69)

Writing in terms of eigenstates

$\displaystyle g_{\mu,\mu^{\prime}}^{(r)}=$ $\displaystyle
Z_{N}^{-1}\sum_{j_{1}=1}^{q}\langle
u_{j_{1}}|\sum_{j_{2},j_{3}=1}^{q}c_{\mu,j_{2}}^{*}c_{\mu,j_{3}}|u_{j_{2}}\rangle\langle
u_{j_{3}}|\lambda_{j_{3}}^{r}$
$\displaystyle\times\sum_{j_{4},j_{5}=1}^{q}c_{\mu^{\prime},j_{4}}^{*}c_{\mu^{\prime},j_{5}}|u_{j_{4}}\rangle\langle
u_{j_{5}}|\lambda_{j_{5}}^{N-r}|u_{j_{1}}\rangle$ (70)

and simplifying (70), we have

$g_{\mu,\mu^{\prime}}^{(r)}=Z_{N}^{-1}\sum_{j_{1}=1}^{q}\sum_{j_{2}=1}^{q}c_{\mu,j_{1}}^{*}c_{\mu,j_{2}}\lambda_{j_{2}}^{r}c_{\mu^{\prime},j_{2}}^{*}c_{\mu^{\prime},j_{1}}\lambda_{j_{1}}^{N-r},$
(71)

in thermodynamic limit $(Z_{N}\rightarrow\lambda_{1}^{N})$. The above
expression reduces to

$g_{\mu,\mu^{\prime}}^{(r)}=\sum_{j_{1}=1}^{q}\sum_{j_{2}=1}^{q}c_{\mu,j_{1}}^{*}c_{\mu,j_{2}}\Bigl{(}\tfrac{\lambda_{j_{2}}}{\lambda_{1}}\Bigr{)}^{r}c_{\mu^{\prime},j_{2}}^{*}c_{\mu^{\prime},j_{1}}\left(\tfrac{\lambda_{j_{1}}}{\lambda_{1}}\right)^{N-r},$
(72)

manipulating conveniently, we have

$g_{\mu,\mu^{\prime}}^{(r)}=c_{\mu,1}^{*}c_{\mu,1}c_{\mu^{\prime},1}^{*}c_{\mu^{\prime},1}+\sum_{j=2}^{q}c_{\mu,1}^{*}c_{\mu,j}\Bigl{(}\tfrac{\lambda_{j}}{\lambda_{1}}\Bigr{)}^{r}c_{\mu^{\prime},j}^{*}c_{\mu^{\prime},1},$
(73)

or even in terms of eq.(43), we have

$g_{\mu,\mu^{\prime}}^{(r)}=m_{\mu}m_{\mu^{\prime}}+c_{\mu,1}^{*}c_{\mu^{\prime},1}\sum_{j=2}^{q}c_{\mu,j}c_{\mu^{\prime},j}^{*}\Bigl{(}\tfrac{\lambda_{j}}{\lambda_{1}}\Bigr{)}^{r}.$
(74)

Therefore, we can write the correlation function as follows

$\Gamma_{\mu,\mu^{\prime}}^{(r)}=c_{\mu,1}^{*}c_{\mu^{\prime},1}\sum_{j=2}^{q}c_{\mu,j}c_{\mu^{\prime},j}^{*}\Bigl{(}\tfrac{\lambda_{j}}{\lambda_{1}}\Bigr{)}^{r}=\\\
=c_{\mu,1}^{*}c_{\mu^{\prime},1}\Bigl{\\{}c_{\mu,2}c_{\mu^{\prime},2}^{*}\left(\tfrac{\lambda_{2}}{\lambda_{1}}\right)^{r}+\left(\tfrac{\lambda_{3}}{\lambda_{1}}\right)^{r}\sum_{j=3}^{q}c_{\mu,j}c_{\mu^{\prime},j}^{*}\Bigr{\\}}.$
(75)

Note that for $q=2$, the last term in (75) ceases to exist.

Taking into account the orthogonality relations for the coefficients in
(13,14,15),

$\sum_{j=1}^{q}c_{\mu,j}c_{\mu^{\prime},j}^{*}=\delta_{\mu,\mu^{\prime}},$
(76)

thus (75) becomes

$\Gamma_{\mu,\mu^{\prime}}^{(r)}=c_{\mu,1}^{*}c_{\mu^{\prime},1}\Bigl{\\{}c_{\mu,2}c_{\mu^{\prime},2}^{*}\left(\tfrac{\lambda_{2}}{\lambda_{1}}\right)^{r}+\\\
+\left(\delta_{\mu,\mu^{\prime}}-c_{\mu,1}c_{\mu^{\prime},1}^{*}-c_{\mu,2}c_{\mu^{\prime},2}^{*}\right)\left(\tfrac{\lambda_{3}}{\lambda_{1}}\right)^{r}\Bigr{\\}}.$
(77)

From (77), and after some algebraic manipulation, we can find the pair
correlation functions in terms of magnetization and transfer matrix
eigenvalues, which has the following form

$\displaystyle\Gamma_{1,1}^{(r)}$
$\displaystyle=m_{1}(1-m_{1})\left(\tfrac{\lambda_{2}}{\lambda_{1}}\right)^{r},$
(78) $\displaystyle\Gamma_{\mu,\mu}^{(r)}$
$\displaystyle=(1-m_{1})\left[m_{1}\left(\tfrac{\lambda_{2}}{\lambda_{1}}\right)^{r}+(q-2)\left(\tfrac{\lambda_{3}}{\lambda_{1}}\right)^{r}\right],$
(79) $\displaystyle\Gamma_{1,\mu}^{(r)}$
$\displaystyle=-\tfrac{m_{1}(1-m_{1})}{q-1}\left(\tfrac{\lambda_{2}}{\lambda_{1}}\right)^{r},$
(80) $\displaystyle\Gamma_{\mu,\mu^{\prime}}^{(r)}$
$\displaystyle=(1-m_{1})\left[m_{1}\left(\tfrac{\lambda_{2}}{\lambda_{1}}\right)^{r}-\left(\tfrac{\lambda_{3}}{\lambda_{1}}\right)^{r}\right],$
(81)

where $\mu,\mu^{\prime}=2,\ldots q$ and $\mu\neq\mu^{\prime}$. An alternative
expression of eqs.(78-81) are given in appendix A.

Figure 9: The pair distribution correlation function
$\Gamma_{\mu,\mu^{\prime}}^{(r)}$ as a function of temperature $T$ in
logarithmic scale, assuming fixed $q=5$, $J=0.9$, $K=1$, $h=-0.995$ and
several values of $r$. (a) For $\mu=\mu^{\prime}$, with $\mu=\\{2,3,\dots\\}$.
(b) For $\mu\neq\mu^{\prime}$ with $\mu,\mu^{\prime}=\\{2,3,\dots\\}$. (c) For
$\mu=\mu^{\prime}=1$.

In Fig. 9 we illustrate the pair correlation function
$\Gamma_{\mu,\mu^{\prime}}^{(r)}$ as a function of temperature in logarithmic
scale, here we consider the following fixed parameters $q=5$, $J=0.9$, $K=1$,
$h=-0.995$ and several distances $r=\\{1,5,10,20,50,100\\}$. In panel (a) we
report for the case of $\mu=\mu^{\prime}$ and $\mu=\\{2,3,\dots\\}$, given by
eq.(79). It is worth notice that for $T<T_{p}$ (enough below) this amount is
closely null, meaning that the pair distribution spins orientations have no
almost relationship, although for $T>T_{p}$ (enough above) the system exhibits
clearly the relationship between pair spins orientations and as expected
decreases with the distances as well as temperature. In panel (b) is depicted
the correlation function for the case $\mu\neq\mu^{\prime}$ and
$\mu,\mu^{\prime}=\\{2,3,4,\dots\\}$, given by eq.(81). For $T<T_{p}$ (enough
below) the $\Gamma_{\mu,\mu^{\prime}}^{(r)}$ become nearly null correlation
function, while for $T>T_{p}$ (enough above) the systems exhibits a
qualitatively different behavior, which in module decreases with with $r$ and
$T$. In panel (c) is displayed for the case of $\mu=\mu^{\prime}=1$ (according
eq.(78)), for this particular case we observe that the correlation function
becomes almost null far enough from the pseudo-critical temperature $T_{p}$,
while at $T_{p}$ the correlation function illustrate a peak, which decreases
as expected with $r$. The case $\mu=1$ and $\mu^{\prime}=\\{2,3,\dots\\}$
given by eq.(81) is simply the same to the case eq.(78) divided by $(q-1)$.

We can also notice that the four expression given by (78-81) satisfy the
following couple of identities

$\displaystyle\Gamma_{1,1}^{(r)}+(q-1)\Gamma_{1,\mu}^{(r)}=$ $\displaystyle
0,$ (82)
$\displaystyle\Gamma_{\mu,\mu}^{(r)}+\Gamma_{1,\mu}^{(r)}+(q-2)\Gamma_{\mu,\mu^{\prime}}^{(r)}=$
$\displaystyle 0.$ (83)

An equivalent expressions of the relations (82-83) are given in appendix A.

Similarly, we can obtain a couple of identities for $g$’s, by using (82) and
(83) which reduce to the following relations

$\displaystyle g_{1,1}^{(r)}+\left(q-1\right)g_{1,\mu}^{(r)}=$ $\displaystyle
m_{1},$ (84) $\displaystyle
g_{\mu,\mu}^{(r)}+g_{1,\mu}^{(r)}+(q-2)g_{\mu,\mu^{\prime}}^{(r)}=$
$\displaystyle m_{\mu},$ (85)

which is useful to studying the phase diagram, like illustrated in Table 1.

Table 1: The nearest neighbour pair distribution functions $g_{\mu,\mu^{\prime}}^{(1)}$, $\mu,\mu^{\prime}=2,\ldots q$, magnetization $m_{1}$, and entropy $\mathcal{S}$ of the frustrated Potts chain at zero temperature. Here $q_{1}=\sqrt{4q-3}$, $q_{2}=\sqrt{(q-1)(q+3)}$, $q_{3}=1+\sqrt{q-1}$, and $q_{4}=\sqrt{(q-1)^{2}+4}$. GS phase | $g_{1,1}^{(1)}$ | $g_{\mu,\mu}^{(1)}$ | $g_{1,\mu}^{(1)}$ | $g_{\mu,\mu^{\prime}}^{(1)}$ | $m_{1}$ | $\mathcal{S}$
---|---|---|---|---|---|---
$FM_{1}$ | $1$ | $0$ | $0$ | $0$ | $1$ | $0$
$FM_{2}$ | $0$ | $\dfrac{1}{q-1}$ | $0$ | $0$ | $0$ | $0$
$FR_{1}$ | $0$ | $0$ | $\dfrac{1}{2\left(q-1\right)}$ | $0$ | $\dfrac{1}{2}$ | $\frac{1}{2}\ln\left(q-1\right)$
$FR_{2}$ | $0$ | $0$ | $0$ | $\dfrac{1}{\left(q-1\right)\left(q-2\right)}$ | 0 | $\ln\left(q-2\right)$
$FR_{1}$-$FM_{1}$ | $\dfrac{1}{q_{1}}$ | $0$ | $\dfrac{q_{1}-1}{2q_{1}\left(q-1\right)}$ | $0$ | $\dfrac{2q-1+q_{1}}{q_{1}\left(1+q_{1}\right)}$ | $\ln\dfrac{1+q_{1}}{2}$
$FR_{2}$-$FM_{2}$ | $0$ | $\dfrac{1}{\left(q-1\right)^{2}}$ | $0$ | $\dfrac{1}{\left(q-1\right)^{2}}$ | $0$ | $\ln\left(q-1\right)$
$FR_{1}$-$FR_{2}$ | $0$ | $0$ | $\dfrac{1}{q\left(q-1\right)}$ | $\dfrac{1}{q\left(q-1\right)}$ | $\dfrac{1}{q}$ | $\ln\left(q-1\right)$
$FR_{1}$-$FM_{2}$ ($K<0$) | $0$ | $\dfrac{1}{q_{1}\left(q-1\right)}$ | $\dfrac{q_{1}-1}{2q_{1}\left(q-1\right)}$ | $0$ | $\dfrac{2q-2}{q_{1}\left(1+q_{1}\right)}$ | $\ln\dfrac{1+q_{1}}{2}$
$FR_{2}$-$FM_{1}$ ($K>0$) | $0$ | $0$ | $0$ | $\dfrac{1}{\left(q-1\right)\left(q-2\right)}$ | $0$ | $\ln\left(q-2\right)$
$FM_{1}$-$FM_{2}$ ($K>0$) | $0$ | $\dfrac{1}{q-1}$ | $0$ | $0$ | $0$ | $0$
$FM_{1}$-$FM_{2}$ ($K=0$) | $\dfrac{1}{q}$ | $\dfrac{1}{q}$ | $0$ | $0$ | $\dfrac{1}{q}$ | $0$
$FM_{1}$-$FM_{2}$ ($K<0$) | $\dfrac{1}{2}$ | $\dfrac{1}{2\left(q-1\right)}$ | $0$ | $0$ | $\dfrac{1}{2}$ | $0$
$Q_{1}$ | $0$ | $\dfrac{1}{q_{2}\left(q-1\right)}$ | $\dfrac{2}{q_{2}\left(q-1+q_{2}\right)}$ | $\dfrac{1}{q_{2}\left(q-1\right)}$ | $\dfrac{2}{q+3+q_{2}}$ | $\ln\dfrac{q-1+q_{2}}{2}$
$Q_{2}$ | $\dfrac{1}{2q_{3}}$ | $\dfrac{1}{2q_{3}\left(q-1\right)}$ | $\dfrac{1}{2q_{3}\sqrt{q-1}}$ | $0$ | $\dfrac{1}{2}$ | $\ln q_{3}$
$P_{1}$ | $\dfrac{3-q+q_{4}}{q_{4}\left(q-1+q_{4}\right)}$ | $0$ | $\dfrac{2}{q_{4}\left(q-1+q_{4}\right)}$ | $\dfrac{q-3+q_{4}}{q_{4}\left(q-1\right)\left(q-1+q_{4}\right)}$ | $\dfrac{q+1+q_{4}}{q_{4}\left(q-1+q_{4}\right)}$ | $\ln\dfrac{q-1+q_{4}}{2}$
$P_{2}$ | $0$ | $\dfrac{1}{\left(q-1\right)^{2}}$ | $0$ | $\dfrac{1}{\left(q-1\right)^{2}}$ | $0$ | $\ln\left(q-1\right)$
$S$ | $\dfrac{1}{q^{2}}$ | $\dfrac{1}{q^{2}}$ | $\dfrac{1}{q^{2}}$ | $\dfrac{1}{q^{2}}$ | $\dfrac{1}{q}$ | $\ln q$

###  IV.1 Correlation length

From transfer matrix eigenvalues, one can observe that
$1>\tfrac{\lambda_{2}}{\lambda_{1}}>\tfrac{\lambda_{3}}{\lambda_{1}}$, thus
for $r>1$, we can ignore
$\left(\tfrac{\lambda_{3}}{\lambda_{1}}\right)^{r}\rightarrow 0$, since
$\left(\tfrac{\lambda_{2}}{\lambda_{1}}\right)^{r}\gg\left(\tfrac{\lambda_{3}}{\lambda_{1}}\right)^{r}$.
Consequently, we can define the correlation length as follows

$\xi=\left[\ln\left(\frac{\lambda_{1}}{\lambda_{2}}\right)\right]^{-1}.$ (86)

An approximation for the correlation length in the vicinity of the pseudo-
transition point immediately follows from equations (30) for the eigenvalues:

$\xi\approx\frac{1}{\sqrt{a^{2}\tau^{2}+b^{2}}}.$ (87)

At the pseudo-transition point the correlation length reaches extremely high
values

$\xi_{p}=b^{-1},$ (88)

since the condition $b\ll 1$ is necessarily satisfied. Indeed, using the
expressions of $b$ in the given limit cases, we obtain

$\xi_{p}=\left\\{\begin{array}[]{ll}\tfrac{1}{2}\sqrt{\tfrac{(q-2)x_{p}k_{p}}{q-1}}\propto
e^{\frac{K-|J|}{2T_{p}}},&qFR_{2}{-}qFM_{1};\\\\[10.00002pt]
\tfrac{x_{p}}{2}\sqrt{\tfrac{k_{p}}{q-1}}\propto
e^{\frac{K+2J}{2T_{p}}},&qFM_{2}{-}qFM_{1}.\end{array}\right.$ (89)

The half-width at half-maximum for the peak (87) of the correlation length is
$\tilde{\Psi}_{\tau}=\sqrt{3}b/a$ and only by numeric factor differs from
$\Psi_{\tau}$.

Figure 10: The correlation length for $qFM_{1}$ states with (a) $q=5$,
$J=-0.5$, $K=1$, (b) $q=3$, $J=-0.05$, $K=1$, (c) $q=5$, $J=0.25$, $K=1$, and
for (d) $qFM_{1}$ and $qFM_{2}$ states with $q=10$, $J=1$, $K=0$, at a given
$h$.

Figure 10 shows the temperature dependences of the correlation length for the
same sets of states as in Fig. 6b, 6d, 7b, and 8b. As for the specific heat
and susceptibility, the correlation length near the pseudo-transition point in
Fig. 10a,b,c shows the giants peaks, which are qualitatively different from
other cases like shown in Fig. 10d.

##  V Pseudo-critical exponents

Now let us analyze the nature of the peaks of the specific heat,
susceptibility, and correlation length around $T_{p}$, whether follow some
critical exponent universality. A general technique to calculate the critical
exponents in the systems that exhibit was developed in Ref. [31]. For the one-
dimensional frustrated $q$-state Potts model, we can find these quantities
directly from approximations (34), (63), and (87). Considering the region of
$\tau$ where the curvature of the peaks becomes positive, $b/a<|\tau|\ll
T_{p}/|J|$, we obtain the following asymptotics

$\displaystyle\xi=$ $\displaystyle\;c_{\xi}|\tau|^{-1},$ $\displaystyle\qquad
c_{\xi}=a^{-1},$ (90) $\displaystyle C=$ $\displaystyle\;c_{f}|\tau|^{-3},$
$\displaystyle\qquad c_{f}=\frac{1}{2}a^{-1}b^{2},$ (91)
$\displaystyle\chi_{1}=$ $\displaystyle\;c_{\chi}|\tau|^{-3},$
$\displaystyle\qquad c_{\chi}=\frac{1}{2T_{p}}a^{-3}b^{2}.$ (92)

These found critical exponents are the same as for one-dimensional models of
the general class with pseudo-transitions [31]. Combining the critical
amplitudes, one can write the following relation

$\frac{c_{f}}{c_{\chi}}\,c_{\xi}^{2}=T_{p},$ (93)

which is fulfilled for all pseudo-transitions in the frustrated Potts model.

Figure 11: (a) The correlation as a function of $|\tau|$ in logarithmic scale,
for $q=5$, $J=0.25$, $K=1$, blue curves correspond to $\tau>0$, orange curves
denote the case $\tau<0$, and dotted lines correspond to $\xi(\tau)$ in
asymptotic limit. (b) The specific heat as a function of $|\tau|$ in
logarithmic scale for the same set of fixed parameters in (a). (c) The
magnetic susceptibility as a function of $|\tau|$ in logarithmic scale, for
the parameters assumed in panel (a).

In Fig. 11 we verify the power law behavior around the peak for some physical
quantities, these results are only valid for the ascending and descending part
of the peak, while this approach fails around the top of the peak. In panel
(a) we report the correlation length as a function of $|\tau|$ in logarithmic
scale assuming the parameters $q=5$, $J=0.25$, $K=1$. Blue curves correspond
for $\tau>0$, orange curves denote for $\tau<0$, and dotted lines describe the
asymptotic behavior given by eq. (90). We consider two values of the magnetic
field as indicated in Fig. 11a. For $h=-0.995$, we can observe clearly a
straight line with pseudo-critical exponent $\nu=1$ in a range of
$10^{-4}\lesssim\tau\lesssim 10^{-1}$ (blue curve). Suchlike behavior we also
observe for the case $-10^{-1}\lesssim\tau\lesssim-10^{-4}$ (orange curve).
Note that for $|\tau|\lesssim 10^{-4}$ the asymptotic behavior fails because
it corresponds to the peak of the curve. Similar behavior is illustrated for
$h=-0.9995$ the pseudo-critical exponent is accurately described by a straight
line with the same exponent $\nu=1$. Although the asymptotic approach is valid
roughly in the interval $10^{-1}\lesssim|\tau|\lesssim 10^{-6}$. In panel (b)
is depicted the specific heat as a function of $|\tau|$ in logarithmic scale,
assuming the same fixed parameters to that considered in panel (a), the
specific heat also exhibits in the ascending a descending part of the peak as
a straight line, which fits precisely to a straight dotted line with angular
coefficient $\alpha=3$ given by eq. (91), although evidenced in the shorter
interval and for smaller $|\tau|$, the straight lines fail for $|\tau|\lesssim
10^{-5}$ ($h=-0.9995$) and $|\tau|\lesssim 5\cdot 10^{-3}$ ($h=-0.995$),
because we are dealing with a peak and not a real singularity. Analogous
behavior is observed for the magnetic susceptibility in panel (c), assuming
the same set of fixed parameters given in panel (a). Once again, we observe
clearly a straight line with angular coefficient $\mu=3$ given by eq. (92),
which is valid for a wider interval compared to that of specific heat.

In summary, the pseudo-critical exponents are independent of the magnetic
field. We also conclude that these exponents satisfy the same universality
class found previously for other one-dimensional models.

##  VI Conclusion

Although one-dimensional systems could appear thoroughly studied, they still
surprise us by exhibiting unconventional new physics. There are some peculiar
one-dimensional models, that exhibit the presence of phase transition at
finite temperature, under the condition of nearest neighbor interaction[21,
22, *galisova20, 24, 25, 26, 27]. One-dimensional models of statistical
physics are attractive from both theoretical and experimental points of view,
often driven to develop new methods in order to solve more realistic models.

Here we investigated more carefully the one-dimensional Potts model with an
external magnetic field and anisotropic interaction that selects neighboring
sites that are in the spin state 1, by using the transfer matrix method. The
largest and the second largest eigenvalues are almost degenerate for a given
temperature, leading to the arising of pseudo-transition. The rise of
anomalous behavior in this model lies in the peculiar behavior of the transfer
matrix elements, all transfer matrix elements are positive, but some off-
diagonal matrix in the low-temperature region can be extremely small compared
to at least two diagonal elements. We have analyzed, the present model for
several physical quantities assuming $K=1$. The entropy and magnetization show
a steep function around pseudo-critical temperature, rather similar to a
first-order phase transition, while the correlation length, specific heat, and
magnetic susceptibility exhibit sharp peaks, around pseudo-critical
temperature, resembling a second-order phase transition, albeit there is no
true divergence. A further investigation of the pseudo-critical exponent
satisfy the same class of universality previously identified for other one-
dimensional models, these exponents are: for correlation length $\nu=1$,
specific heat $\alpha=3$ and magnetic susceptibility $\mu=3$. Whereas for
$K=0$ (standard Potts chain), we observe a qualitatively different behavior,
such as in entropy there is no step-like function, there is no sharp peak in
specific heat, while broad peak arises for magnetic susceptibility.

It is worthy to mention that the pseudo-transition is quite different from
that true phase transition because there is no jump in the first derivative of
free energy, nor divergence in the second derivative of free energy. In this
sense, it would be fairly relevant to observe this anomalous property
experimentally in chemical compounds.

###### Acknowledgements.

The work was partly supported by the Ministry of Education and Science of RF,
project No FEUZ-2020-0054, and Brazilian agency CNPq.

##  Appendix A Some addition relations

Here we give some additional alternative expressions, which would be useful
for further analysis, thus the eqs. (78-81) can be reduced to

$\displaystyle\Gamma_{1,1}^{(r)}$
$\displaystyle=\tfrac{\sin^{2}(2\phi)}{4}\left(\tfrac{\lambda_{2}}{\lambda_{1}}\right)^{r},$
(94) $\displaystyle\Gamma_{\mu,\mu}^{(r)}$
$\displaystyle=\tfrac{\Gamma_{1,1}^{(r)}}{(q-1)^{2}}+\tfrac{(q-2)\sin^{2}(\phi)}{(q-1)^{2}}\left(\tfrac{\lambda_{3}}{\lambda_{1}}\right)^{r},$
(95) $\displaystyle\Gamma_{1,\mu}^{(r)}$
$\displaystyle=-\tfrac{\Gamma_{1,1}^{(r)}}{q-1},$ (96)
$\displaystyle\Gamma_{\mu,\mu^{\prime}}^{(r)}$
$\displaystyle=\tfrac{\Gamma_{1,1}^{(r)}}{(q-1)^{2}}-\tfrac{\sin^{2}(\phi)}{(q-1)^{2}}\left(\tfrac{\lambda_{3}}{\lambda_{1}}\right)^{r}.$
(97)

Considering the identities given by (82) and (83), we obtain the following
identity which must satisfy:

$\Gamma_{1,1}^{(r)}+\left(q-1\right)\Gamma_{\mu,\mu}^{(r)}+2\left(q-1\right)\Gamma_{1,\mu}^{(r)}+{}\\\
+\left(q-1\right)\left(q-2\right)\Gamma_{\mu,\mu^{\prime}}^{(r)}=0.$ (98)

An equivalent relation we can verify, so the pair average distribution
functions obey the identity

$g_{1,1}^{(r)}+\left(q-1\right)g_{\mu,\mu}^{(r)}+2\left(q-1\right)g_{1,\mu}^{(r)}+{}\\\
+\left(q-1\right)\left(q-2\right)g_{\mu,\mu^{\prime}}^{(r)}=1.$ (99)

Below we simplify for the nearest neighbors PDF $g_{\mu,\mu^{\prime}}^{(1)}$,
which can be written as

$\displaystyle g_{1,1}^{(1)}$
$\displaystyle=\tfrac{\left(\lambda_{1}-\lambda_{3}-(q-1)t_{2}\right)\left(\lambda_{1}+\lambda_{2}-\lambda_{3}-(q-1)t_{2}\right)}{\lambda_{1}\left(\lambda_{1}-\lambda_{2}\right)},$
(100) $\displaystyle g_{\mu,\mu}^{(1)}$
$\displaystyle=\tfrac{\left(\lambda_{3}+t_{2}\right)\left(\lambda_{3}-\lambda_{2}+(q-1)t_{2}\right)}{\left(q-1\right)\lambda_{1}\left(\lambda_{1}-\lambda_{2}\right)},$
(101) $\displaystyle g_{1,\mu}^{(1)}$
$\displaystyle=\tfrac{\left(\lambda_{1}-\lambda_{3}-(q-1)t_{2}\right)\left(\lambda_{3}-\lambda_{2}+(q-1)t_{2}\right)}{\left(q-1\right)\lambda_{1}\left(\lambda_{1}-\lambda_{2}\right)},$
(102) $\displaystyle g_{\mu,\mu^{\prime}}^{(1)}$
$\displaystyle=\tfrac{\left(\lambda_{3}-\lambda_{2}+(q-1)t_{2}\right)t_{2}}{\left(q-1\right)\lambda_{1}\left(\lambda_{1}-\lambda_{2}\right)},$
(103)

where $\mu,\mu^{\prime}=2,\ldots q$, and $\mu\neq\mu^{\prime}$. This amounts
are useful to analysis the phase boundary properties illustrated in Table 1.

## References

  * Baxter [1982] R. J. Baxter,  _Exactly solved models in statistical mechanics_ (Academic Press, London ; New York, 1982).
  * Ghulghazaryan  _et al._ [2007] R. G. Ghulghazaryan, K. G. Sargsyan, and N. S. Ananikian, Partition function zeros of the one-dimensional Blume-Capel model in transfer matrix formalism, Physical Review E 76, 021104 (2007).
  * Hovhannisyan  _et al._ [2009] V. Hovhannisyan, R. Ghulghazaryan, and N. Ananikian, The partition function zeros of the anisotropic Ising model with multisite interactions on a zigzag ladder, Physica A: Statistical Mechanics and its Applications 388, 1479 (2009).
  * van Hove [1950] L. van Hove, Sur L’intégrale de Configuration Pour Les Systèmes De Particules À Une Dimension, Physica 16, 137 (1950).
  * Ninio [1976] F. Ninio, A simple proof of the Perron-Frobenius theorem for positive symmetric matrices, Journal of Physics A: Mathematical and General 9, 1281 (1976).
  * Cuesta and Sánchez [2004] J. A. Cuesta and A. Sánchez, General Non-Existence Theorem for Phase Transitions in One-Dimensional Systems with Short Range Interactions, and Physical Examples of Such Transitions, Journal of Statistical Physics 115, 869 (2004).
  * Kittel [1969] C. Kittel, Phase Transition of a Molecular Zipper, American Journal of Physics 37, 917 (1969).
  * Chui and Weeks [1981] S. T. Chui and J. D. Weeks, Pinning and roughening of one-dimensional models of interfaces and steps, Physical Review B 23, 2438 (1981).
  * Dauxois and Peyrard [1995] T. Dauxois and M. Peyrard, Entropy-driven transition in a one-dimensional system, Physical Review E 51, 4027 (1995).
  * Sarkanych  _et al._ [2017] P. Sarkanych, Y. Holovatch, and R. Kenna, Exact solution of a classical short-range spin model with a phase transition in one dimension: The Potts model with invisible states, Physics Letters A 381, 3589 (2017).
  * Zimm and Bragg [1959] B. H. Zimm and J. K. Bragg, Theory of the Phase Transition between Helix and Random Coil in Polypeptide Chains, The Journal of Chemical Physics 31, 526 (1959).
  * Badasyan  _et al._ [2010] A. V. Badasyan, A. Giacometti, Y. S. Mamasakhlisov, V. F. Morozov, and A. S. Benight, Microscopic formulation of the Zimm-Bragg model for the helix-coil transition, Physical Review E 81, 021921 (2010).
  * Ananikyan  _et al._ [1990] N. S. Ananikyan, S. A. Hajryan, E. S. Mamasakhlisov, and V. F. Morozov, Helix-Coil transition in polypeptides: A microscopical approach, Biopolymers 30, 357 (1990).
  * Badasyan  _et al._ [2013] A. Badasyan, A. Giacometti, R. Podgornik, Y. Mamasakhlisov, and V. Morozov, Helix-coil transition in terms of Potts-like spins, The European Physical Journal E 36, 46 (2013).
  * Tonoyan  _et al._ [2020] S. Tonoyan, D. Khechoyan, Y. Mamasakhlisov, and A. Badasyan, Statistical mechanics of DNA-nanotube adsorption, Physical Review E 101, 062422 (2020).
  * Wu [1982] F. Y. Wu, The Potts model, Reviews of Modern Physics 54, 235 (1982).
  * Aharony and Pfeuty [1979] A. Aharony and P. Pfeuty, Dilute spin glasses at zero temperature and the ${}^{\textrm{1}}$ / ${}_{\textrm{2}}$ -state Potts model, Journal of Physics C: Solid State Physics 12, L125 (1979).
  * Lubensky and Isaacson [1978] T. C. Lubensky and J. Isaacson, Field Theory for the Statistics of Branched Polymers, Gelation, and Vulcanization, Physical Review Letters 41, 829 (1978).
  * Chaves and Riera [1984] C. M. Chaves and R. Riera, Dilute Potts chain in a magnetic field, Canadian Journal of Physics 62, 69 (1984).
  * Panov [2020] Y. Panov, Local distributions of the 1D dilute Ising model, Journal of Magnetism and Magnetic Materials 514, 167224 (2020).
  * Gálisová and Strečka [2015] L. Gálisová and J. Strečka, Vigorous thermal excitations in a double-tetrahedral chain of localized Ising spins and mobile electrons mimic a temperature-driven first-order phase transition, Physical Review E 91, 022134 (2015).
  * Gálisová [2017] L. Gálisová, Magnetization plateau as a result of the uniform and gradual electron doping in a coupled spin-electron double-tetrahedral chain, Physical Review E 96, 052110 (2017).
  * Gálisová [2020] L. Gálisová, Pairwise Entanglement in Double-Tetrahedral Chain with Different Landé g-Factors of the Ising and Heisenberg Spins, Acta Physica Polonica A 137, 604 (2020).
  * Strečka  _et al._ [2016] J. Strečka, R. C. Alécio, M. L. Lyra, and O. Rojas, Spin frustration of a spin-1/2 Ising–Heisenberg three-leg tube as an indispensable ground for thermal entanglement, Journal of Magnetism and Magnetic Materials 409, 124 (2016).
  * Rojas  _et al._ [2016] O. Rojas, J. Strečka, and S. de Souza, Thermal entanglement and sharp specific-heat peak in an exactly solved spin-1/2 Ising-Heisenberg ladder with alternating Ising and Heisenberg inter–leg couplings, Solid State Communications 246, 68 (2016).
  * Strečka [2020] J. Strečka, Anomalous Thermodynamic Response in the Vicinity of a Pseudo-Transition of a Spin-1/2 Ising Diamond Chain, Acta Physica Polonica A 137, 610 (2020).
  * de Souza and Rojas [2018] S. de Souza and O. Rojas, Quasi-phases and pseudo-transitions in one-dimensional models with nearest neighbor interactions, Solid State Communications 269, 131 (2018).
  * Carvalho  _et al._ [2019] I. Carvalho, J. Torrico, S. de Souza, O. Rojas, and O. Derzhko, Correlation functions for a spin- 1 2 Ising-XYZ diamond chain: Further evidence for quasi-phases and pseudo-transitions, Annals of Physics 402, 45 (2019).
  * Rojas [2020a] O. Rojas, Residual Entropy and Low Temperature Pseudo-Transition for One-Dimensional Models, Acta Physica Polonica A 137, 933 (2020a).
  * Rojas [2020b] O. Rojas, A Conjecture on the Relationship Between Critical Residual Entropy and Finite Temperature Pseudo-transitions of One-dimensional Models, Brazilian Journal of Physics 50, 675 (2020b).
  * Rojas  _et al._ [2019] O. Rojas, J. Strečka, M. L. Lyra, and S. M. de Souza, Universality and quasicritical exponents of one-dimensional models displaying a quasitransition at finite temperatures, Physical Review E 99, 042117 (2019).
  * Janke and Schakel [2004] W. Janke and A. M. Schakel, Geometrical vs. Fortuin–Kasteleyn clusters in the two-dimensional q-state Potts model, Nuclear Physics B 700, 385 (2004).
  * Qian  _et al._ [2005] X. Qian, Y. Deng, and H. W. J. Blöte, Dilute Potts model in two dimensions, Physical Review E 72, 056132 (2005).


# Synchronization of Dust Acoustic Waves in a forced Korteweg-de Vries-Burgers
model

Ajaz Mir<EMAIL_ADDRESS>Indian Institute of Technology Jammu, Jammu,
J&K, 181221, India  Sanat Tiwari<EMAIL_ADDRESS>Indian Institute
of Technology Jammu, Jammu, J&K, 181221, India  Abhijit Sen  Institute for
Plasma Research, Gandhinagar, Gujarat, 382428, India  Chris Crabtree  Naval
Research Laboratory, Washington, DC 20375, USA  Gurudas Ganguli  Naval
Research Laboratory, Washington, DC 20375, USA  John Goree  Department of
Physics and Astronomy, University of Iowa, Iowa City, IA 52242, USA

###### Abstract

The synchronization of dust acoustic waves to an external periodic source is
studied in the framework of a driven Korteweg-de Vries-Burgers equation that
takes into account the appropriate nonlinear and dispersive nature of low
frequency waves in a dusty plasma medium. For a spatio-temporally varying
source term the system is shown to demonstrate harmonic (1:1) and super-
harmonic (1:2) synchronized states. The existence domains of these states are
delineated in the form of Arnold tongue diagrams in the parametric space of
the forcing amplitude and forcing frequency and their resemblance to some past
experimental results is discussed.

##  I Introduction

###

The nonlinear phenomenon of frequency synchronization is ubiquitous in many
physical, chemical, and biological systems and has been the subject of a large
number of studies over the past several years Pikovsky  _et al._ (2001);
Strogatz (2004); Balanov  _et al._ (2009). The simplest mathematical model
describing this phenomenon consists of an ensemble of globally coupled
nonlinear point oscillators that adjust their intrinsic frequencies to a
common collective frequency as the coupling strength is increased Acebrón  _et
al._ (2005); Kuramoto (1984); Aronson  _et al._ (1990); Kuramoto (2003). Such
a nonlinear phenomenon can also be observed in a continuum medium (a fluid)
where a self-excited oscillation or a wave can interact with a driving force
and adjust its oscillation or wave frequency Williams and Duff (2010); Block
_et al._ (2001); Klinger  _et al._ (1993, 1995); Gyergyek (1999); Nurujjaman
and Iyengar (2010). A plasma system with its wide variety of collective modes
and complex nonlinear dynamics provides a rich and challenging medium for the
exploration of synchronization phenomena. A number of past experimental
studies have examined the driven response of a plasma to an external frequency
source Ruhunusiri and Goree (2012); Pilch  _et al._ (2009); Klinger  _et al._
(1995); Williams (2014). There have also been a few studies devoted to an
investigation of mutual synchronization between two plasma devices Chaubey
_et al._ (2015); Fukuyama  _et al._ (2006); Chaubey  _et al._ (2016).

###

More recently, synchronization phenomena have been experimentally explored in
dusty plasma devices where it is easy to visualize the low-frequency wave
activity using fast video imaging. A dusty plasma is a four-component plasma
of electrons, ions, neutral gas atoms, and micron-size particles of solid
matter Shukla and Eliasson (2009); Morfill and Ivlev (2009); Shukla and Mamun
(2001). It can be produced in a laboratory device like a glow discharge
plasma, by introducing micron sized solid particles Thompson  _et al._ (1997);
Schwabe  _et al._ (2007); Barkan  _et al._ (1995); Heinrich  _et al._ (2009).
These small solid particles (dust) get negatively charged by absorbing more
electrons which have a higher mobility than ions. Such a charged medium
consisting of dust, ions and electrons, can sustain a variety of collective
modes Bandyopadhyay  _et al._ (2008); Deka  _et al._ (2017); Piel and Melzer
(2001); Shukla and Mamun (2001). The dust acoustic wave (DAW) or dust density
wave (DDW) is one such well known low frequency compressional mode Liu  _et
al._ (2018); Shukla and Mamun (2001). A DAW can be spontaneously excited due
to the onset of an ion-streaming instability. The DAW has a very low frequency
(typically 10–100 Hz) Ruhunusiri and Goree (2012); Thompson  _et al._ (1997)
due to the large mass of the dust particles and can consequently be visually
observed; through its images and video recording Thompson  _et al._ (1999);
Don  _et al._ (2014); Feng  _et al._ (2007); Schwabe  _et al._ (2007).

###

Synchronization of dust acoustic waves has been studied in an anodic plasma
(Pilch  _et al._ , 2009), radio-frequency (RF) and direct-current (DC) plasmas
(Ruhunusiri and Goree, 2012; Deka  _et al._ , 2020; Williams, 2018). Pilch et
al. (Pilch  _et al._ , 2009) reported the entrainment of DAWs through a
driving modulation to the anode. Ruhunusiri et al. (Ruhunusiri and Goree,
2012) reported observation of harmonic, super-harmonic, and sub-harmonic
synchrony of self-excited cnoidal DAWs. This was achieved through the driven
modulation of the streaming ions in the dust cloud. Their experiments showed
parametric regions for the occurrence of such synchrony in the form of Arnold
tongue diagrams in the state space of the driving frequency and driving
amplitude. They also observed features like the branching of the tongues and
the existence of an amplitude threshold for synchronization to occur. Williams
et al. (Williams, 2018) compared DAW synchronization in RF and DC generated
plasmas. Their results suggested that in a RF plasma, synchronization was
restricted to a part of the dust cloud volume unlike the complete dust cloud
synchrony in a DC discharge plasma. Deka et al. Deka  _et al._ (2020) observed
the synchronization of self-excited DDW, through the suppression mechanism, by
modulating ion streaming using an external sinusoidal driver. Recently, Liu et
al. (Liu  _et al._ , 2021) carried out experiments in the Plasma Kristall-4
(PK-4) device on board the International Space Station (ISS) under micro-
gravity conditions and reported phase locking for harmonic synchronization.

###

Theoretical efforts towards interpretation and physical understanding of these
experimental results have so far been limited to providing qualitative
comparisons with results obtained from very simple dynamical models. One of
the commonly employed mathematical model is the periodically driven van der
Pol oscillator Van Der Pol (1927); Pikovsky  _et al._ (2001); Balanov  _et
al._ (2009),

$\frac{d^{2}x}{dt^{2}}+(c_{1}-c_{2}x^{2})\frac{dx}{dt}+\omega_{0}^{2}x=A_{dr}\cos(2\pi
f_{dr}t)$ (1)

which describes the displacement $x$ of a harmonic oscillator with a natural
frequency $\omega_{0}$, with terms for a nonlinear damping $c_{2}x^{2}dx/dt$,
a source of energy for self-excitation $c_{1}dx/dt$, and a periodic driving
source at a frequency $f_{dr}$. The van der Pol oscillator can exhibit
synchronization not only at $f_{dr}/f_{0}\approx 1$, which is called
“harmonic” synchronization, but at ratios that are rational numbers. If
$f_{dr}/f_{0}>1$, the synchronization is said to be “super-harmonic”, whereas
if $f_{dr}/f_{0}<1$ it is “sub-harmonic”. Often, the van der Pol oscillator
model has been used as a reference for characterizing synchronization
phenomena in plasmas and other media that support the propagation of waves
Kovalev and Kovalev (2017); Chaubey  _et al._ (2015); Ruhunusiri and Goree
(2012); Deka  _et al._ (2020). However, as a point oscillator model its
dynamics is restricted to nonlinear oscillations and it cannot correctly
represent nonlinear waves. This is also evident from the fact that the van der
Pol model is an ordinary differential equation in time and therefore has no
spatial dynamics that characterizes a propagating wave. In addition, for
nonlinear ion acoustic or dust acoustic waves dispersion plays an important
role in defining their propagation characteristics. Hence there is a need to
develop appropriate theoretical models to describe nonlinear synchronization
of such dispersive waves.

###

In this paper, we present such a model and use it to demonstrate
synchronization of dust acoustic waves. The model is a generalization of the
forced Korteweg-de Vries (fKdV) model that was developed by Sen et al. Sen
_et al._ (2015) for driven nonlinear acoustic waves and subsequently modified
and extensively used to study nonlinear precursor solitons in dusty plasma
experiments Kumar Tiwari and Sen (2016); Jaiswal  _et al._ (2016a). In the
original fKdV model Sen  _et al._ (2015) the driving term was a constant (non-
oscillatory). For our synchronization application we convert this to an
oscillatory form with both a temporal and spatial periodicity. We also
introduce a viscosity term to take account of dissipative processes that are
important in most laboratory studies of dusty plasmas (Nakamura  _et al._ ,
1999; Nakamura and Sarma, 2001). The resulting equation has the form of a
forced Korteweg-de Vries-Burgers (fKdV-B) equation. This equation provides a
proper theoretical framework for the study of synchronization in a realistic
dispersive plasma system that includes natural growth and dissipation of
waves.

###

Numerical solution of this model equation, show clear signatures of harmonic
(1:1) and super-harmonic (1:2) synchronization. The characteristic features of
the synchronization are delineated using power spectral density (PSD) plots,
phase space plots and Lissajous plots obtained from the time-series data
collected at one spatial location. A parametric plot in the form of an Arnold
tongue diagram shows multiple tongues, each corresponding to the existence
region of a harmonic or a higher-order super-harmonic state. The tongues also
show a branching behaviour.

###

The rest of the paper is organized as follows. Section II briefly describes
the fKdV-B model and the numerical approach adopted to solve it. The section
also presents some numerical results for the undriven KdV and KdV-B equations
as background information on the characteristic nonlinear features of the
waves and to describe the diagnostic tools to be used for identifying
synchronization phenomena. Section III presents our main results on harmonic
and super-harmonic synchronization using the fKdV-B model. A brief summary and
some concluding discussion are provided in section IV.

##  II The fKdV-B equation and the numerical approach

###

The fKdV-B equation, a one-dimensional driven nonlinear partial differential
equation, is of the form:

$\displaystyle\frac{\partial n(x,t)}{\partial t}+\ \beta n(x,t)\frac{\partial
n(x,t)}{\partial x}$ $\displaystyle+\ \gamma\frac{\partial^{3}n(x,t)}{\partial
x^{3}}-\ \eta\frac{\partial^{2}n(x,t)}{\partial x^{2}}=\ F_{s}(x,t).$ (2)

Figure 1:  PSD for the times-series of KdV equation with initial perturbation
Eq. (4). Inset (II) shows the time-series obtained from the evolution of the
KdV equation. Inset (I) is the phase space plot for the same time-series.
Figure 2:  PSD of the evolution of KdV-B equation with initial perturbation
Eq. (4) at $\eta=0.0025$. The time-series of the KdV-B equation is shown in
the inset (II). Inset (I) is the phase space plot of the same KdV-B time-
series. Figure 3:  The harmonic (1:1) synchronization in the fKdV-B model with
$f_{s}<f_{0}$ and $f_{s}>f_{0}$. The time-series of the fKdV-B model (blue
bold line) and the forcing (brown dash-dotted line) at driver frequency (a)
$f_{s}=6.8$ Hz with threshold amplitude $A_{s}=0.15A_{0}$ and (b) $f_{s}=7.0$
Hz with threshold amplitude $A_{s}=0.15A_{0}$. (c) PSD of times-series (a).
(d) PSD of time-series (b). The insets (I) phase space plot and (II) Lissajous
figure reflect the frequency locking at the driver frequency.

Here $n(x,t)$ is the dependent variable (the perturbed density in this case)
and $F_{s}(x,t)$ is an external spatio-temporal forcing term. $\beta$,
$\gamma$, and $\eta$ are positive quantities representing the strength of
nonlinearity, dispersion, and viscous damping, respectively. The spatial
coordinate $x$ and time $t$ are normalized by the plasma Debye length
$\lambda_{D}$ and the dust plasma period $\omega_{pd}^{-1}$, respectively.

###

The fKdV-B equation can be derived from the full fluid-Poisson set of
equations, in the weakly nonlinear, dispersive and dissipative regime by using
a reductive perturbation method (Nakamura  _et al._ , 1999; Veeresha  _et al._
, 2010). Such a derivation in the absence of the viscosity term has been given
in detail by Sen et al. Sen  _et al._ (2015). In their derivation the source
term was taken to be a constant; here it is replaced by a more general term
that is a periodic function in both space and time. The KdV-B equation (i.e.,
Eq. (2) in the absence of the driving term) is well known in the literature
(Nakamura  _et al._ , 1999; Veeresha  _et al._ , 2010; Nakamura and Sarma,
2001) and has been employed in the past to model oscillatory shocks in dusty
plasmas Nakamura  _et al._ (1999); Jaiswal  _et al._ (2016b). The model has
also been used to study temporal chaos or spatial chaos by using a randomly
time varying (Orłowski, 1994) or randomly space varying (Rech, 2013) driving
term. In this work, we use a spatio-temporally varying periodic source and
carry out a numerical investigation of Eq. (2) to study the synchronization of
DAWs.

###

The driving source is taken to be of the form of a cnoidal-square travelling
wave,

$F_{s}(x,t)=A_{s}cn^{2}[2K(\kappa_{s})\\{x/\lambda_{s}-f_{s}t\\};\
\kappa_{s}]$ (3)

where $cn$ is the Jacobi elliptic function, $A_{s}$ is the driving amplitude,
$\lambda_{s}$ is the spatial wave length and $f_{s}$ is the driving frequency.
$K(\kappa)$ is the complete elliptical integral of first kind and the elliptic
parameter $\kappa$ is a measure of the nonlinearity of the wave. The cnoidal-
square travelling wave is an exact solution of the KdV equation (i.e., a
solution of the model equation in the absence of viscosity and the driving
term). It can therefore mimic the driving of the system by a DAW arising from
an external (coupled) plasma source. For the numerical solution of Eq. (2) the
initial waveform is also taken to be of the form,

$n(x,t=0)=A_{0}cn^{2}[2K(\kappa_{0})\\{x/\lambda_{0}\\};\ \kappa_{0}].$ (4)

with the values of $A_{0}$, $f_{0}$ and $\lambda_{0}$ different from those of
the driving source. The idea is to see whether the final driven modes of the
system synchronize to the frequency of the driver. Equation (2) is solved for
various values of $f_{s}$ and $A_{s}$ in order to find the regions of
synchronization in the parameter space of $(A_{s},f_{s})$.

Figure 4:  The super-harmonic (1:2) synchronization in the fKdV-B model with
$f_{s}<2f_{0}$ and $f_{s}>2f_{0}$. The time-series of the fKdV-B model (blue
bold line) and the forcing (brown dash-dotted line) at driver frequency (a)
$f_{s}=13.6$ Hz with threshold amplitude $A_{s}=0.35A_{0}$ and (b)
$f_{s}=14.0$ Hz with threshold amplitude $A_{s}=0.45A_{0}$. (c) PSD of times-
series (a). (d) PSD of time-series (b). The insets (I) phase space plot and
(II) Lissajous figure reflect the frequency locking at half of the driver
frequency.

###

Our numerical investigation of the fKdV-B equation is based on the pseudo-
spectral method (Boyd, 2003) and uses periodic boundary conditions. The code
is first bench marked by reproducing earlier results Sen  _et al._ (2015); Mir
_et al._ (2022) obtained for the fKdV equation. The various parameter values
associated with the model are taken to be as follows: The Jacobi elliptic
parameters $\kappa_{0}=\kappa_{s}=0.98$ for Eqs. (3) and (4). The wave vector
of the initial perturbation i.e., $k_{0}=16k_{m}$ where $k_{m}=(2\pi)/L_{x}$
being the minimum wave vector associated with a system of length $L_{x}=6\pi$.
The corresponding wavelength i.e., $\lambda_{0}=(2\pi)/k_{0}$ and amplitude
$A_{0}=10$ are kept fixed throughout the analysis. We have taken
$k_{s}=16k_{m}$ and $k_{s}=2\times 16k_{m}$ for studying harmonic (1:1), and
super-harmonic (1:2) synchronization states. The corresponding forcing
wavelength is $\lambda_{s}=(2\pi)/k_{s}$. Throughout the analysis, we have
only varied the forcing amplitude, $A_{s}$ and forcing frequency, $f_{s}$. We
have also chosen the coefficient values of the model equation as $\beta=2.5$,
$\gamma=0.0667$ and $\eta=0.0025$ for all the runs. Based on the chosen
parameters $\alpha$, $\beta$, $\kappa_{0}$ and $k_{0}$, the value of the
initial amplitude $(A_{0}=10)$ is derived using the relationship provided in
Mir et al. Mir  _et al._ (2020). The amplitude of the initial perturbation
chosen in this fashion will be governed by the exact solution of the KdV and
will be a stable solution of KdV for this particular amplitude.

###

We evolve the initial perturbation in Eq. (2) over long times for these
various different parameter values. During the spatio-temporal evolution, we
collect a time-series of the density field at a fixed spatial location and use
it to calculate the power spectral density. The PSD provides a useful tool for
distinguishing between synchronized and un-synchronized states.

###

As an illustrative example, we show in Fig. 1 the PSD, the time-series and the
phase space plot of the solution, obtained for a pure KdV equation (i.e., for
$\eta=A_{s}=0$). The time-series data has been collected up to $t_{max}=40$
$\omega_{pd}^{-1}$ with a time step $dt=10^{-5}$ $\omega_{pd}^{-1}$. The
maximum sampling frequency $f_{S}=1/dt$ and the Nyquist frequency is
$f_{N}=f_{S}/2$. This leads to a frequency resolution of $df=1/t_{max}$ for
the collected time-series. The time-series data corresponding to the first 20
periods is discarded to remove transient effects while constructing the PSD.
In Fig. 1 the nonlinear character of the mode is evident from the presence of
the higher harmonics in the PSD and from the shape of wave form in the time-
series. The single cycle phase space plot with its form resembling a
separatrix curve indicates an undamped nonlinear periodic wave, in this case
the exact cnoidal-square wave.

###

For comparison, we present in Fig. 2 the corresponding results for the
undriven KdV-B equation (with $\eta=0.0025$ and $F_{s}=0$). The effect of
viscous damping is seen in the frequency shift of the fundamental component in
the PSD towards a lower value of $f_{0}^{\eta}=4.6$ Hz, the reduced amplitude
in the time-series and the spiralling of the phase space plot towards the
origin. It is clear that in the presence of finite viscosity the cnoidal-
square wave can no longer be sustained as a nonlinear solution of Eq. (2) with
$F_{s}=0$ and the initial perturbation decays in time. The question is whether
by driving the system with a periodic source one can revive and sustain a
nonlinear solution that is also synchronized with the driver. The answer is in
the positive and we next present our results on such a phenomenon.

Figure 5:  The Arnold tongue diagrams for (a) harmonic (1:1) synchronization
state with $\eta=0.005$ and (b) harmonic (1:1) and super-harmonic (1:2)
synchronization states with $\eta=0.0025$ in the fKdV-B model. In subplot (a)
the amplitude is varied from $A_{s}=0.25A_{0}$ to $A_{s}=0.50A_{0}$ for 1:1
synchronization. In subplot (b) the amplitude is varied from $A_{s}=0.15A_{0}$
to $A_{s}=0.50A_{0}$ for 1:1, and $A_{s}=0.35A_{0}$ to $A_{s}=0.50A_{0}$ for
1:2 synchronization.

##  III Synchronization in fKdV-B model

In this section we present the main results of our work, namely, the
synchronization of the solutions of Eq. (2) to an external driver of the form
given by Eq. (3). We begin by discussing harmonic (1:1) synchronization for
which we choose the driving frequency to be slightly away from the fundamental
frequency of $f_{0}=6.9$ Hz that is characteristic fundamental frequency of
the undriven system. Two cases are considered, namely, $f_{s}=6.8$ and
$f_{s}=7.0$ Hz. The driving amplitude in both cases is taken to be $A_{s}=15$
with $A_{0}=10$ and $\eta=0.0025$. Figure 3 shows the attainment of harmonic
(1:1) synchronization for both these cases with subplots (a, b) devoted to
$f_{s}=6.8$ Hz and (c, d) to $f_{s}=7.0$ Hz, respectively. As can be seen from
the time-series plots in (a) and (c) the driven solutions are indeed locked to
the driver. This is also clearly seen in the PSDs where the fundamental
frequencies of the driven solutions are indeed at the frequency of the driver.
Furthermore, the phase space plots in (b) and (d) show that these solutions
constitute undamped nonlinear periodic waves that are maintained by a balance
between the nonlinear steepening, dispersive broadening, viscous damping and
amplification due to the external pumping by the driving term. The resultant
phase space curve, that has the characteristic shape of a separatrix,
represents a stationary cnoidal wave solution. The presence of dissipation
seems to be necessary for sustaining this synchronized driven solution. We
have found that in the partial differential equation Eq. (2), including not
just nonlinear and dispersive terms, but also a linear dissipative term,
allowed achieving synchronization of a wave. When we turned off dissipation,
by setting the viscosity coefficient to zero in Eq. (2), we did not observe
synchronization of the wave, for the conditions that we studied here. This is
different from the case of a point oscillator, as described by the van der Pol
oscillator Eq. (1), which requires a nonlinear dissipation term to obtain
synchronization. In the absence of viscosity one only gets nonlinear mixing
from the model as has been reported earlier in Mir et al. Mir  _et al._ (2020,
2022). The amount of viscosity also determines the threshold condition for the
driver amplitude.

###

To explore super-harmonic (1:2) synchronization we again consider two cases of
$f_{s}=13.6$ Hz and $f_{s}=14$ Hz which are slightly below and above the first
harmonic frequency $2f_{0}=13.8$ Hz of the undriven system. The results are
shown in Fig. 4 where the subplots (a,c) are devoted to $f_{s}=13.6$ Hz and
(b,d) to $f_{s}=14$ Hz respectively. As in the previous case of harmonic
synchronization, we see clear evidence of super-harmonic (2:1) synchronization
in the time-series plots, the PSDs and the phase space plots. The Lissajous
figures have a number eight-like trajectory which is indicative of a (1:2)
synchronized state. One significant difference from the harmonic
synchronization case is that the minimum threshold amplitude for the driver to
achieve a 1:2 state is different for the cases $f_{s}<2f_{0}$ and
$f_{s}>2f_{0}$. They are $A_{s}=0.35A_{0}$ and $A_{s}=0.45A_{0}$ respectively.

###

Finally, in Fig. 5 we present a consolidated picture of the existence domain
of these synchronized states in the parameter space of the driver frequency
$f_{s}$ and driver amplitude $A_{s}$ in form of an Arnold tongue diagram. To
obtain the Arnold tongue diagram, $A_{s}$ is varied in steps of 0.5 (which is
$0.05A_{0}$) from 0 to 5 while $f_{s}$ is varied in steps of 0.1 Hz for
harmonic synchronization and 0.2 Hz for the super-harmonic case. Figure 5(a)
shows the 1:1 entrained state tongue for $\eta=0.005$ and Fig. 5(b) shows the
1:1 and 1:2 entrained state tongues for $\eta=0.0025$.

###

We observe several interesting features in the Arnold tongue diagrams. To
start with, there is always a threshold amplitude $A_{s}$ below which no
synchronization occurs. For the harmonic (1:1) synchronization it is
$A_{s}=0.25A_{0}=2.5$ for $\eta=0.005$, while it is $A_{s}=0.15A_{0}=1.5$ for
$\eta=0.0025$. This is unlike the harmonic synchronization phenomenon observed
in a driven van der Pol model where no such threshold is found Jensen (2002).
Another important feature is a distinctive branching of the Arnold tongue that
is clearly seen for the (1:1) states at low forcing amplitudes. Arrows mark
this branching both in Fig. 5(a) and Fig. 5(b). The branching gives rise to a
non-synchronized region between the frequencies $f_{s}=6.4$ Hz to $f_{s}=5.7$
Hz at driver amplitude $A_{s}=0.25A_{0}=2.5$ in Fig. 5(a). This branching
narrows down with the increase in $A_{s}$. Another branch is seen in Fig. 5(a)
between $f_{s}=5.4$ Hz and $f_{s}=5.3$ Hz which also narrows down with
increase in $A_{s}$. A third feature is the asymmetric nature of the tongue
structures about $f_{0}$. The frequency width over which synchronization can
be obtained is much broader for $f_{s}<f_{0}$ compared to $f_{s}>f_{0}$.

##  IV Summary and Conclusions

###

To summarize, we have studied the phenomenon of synchronization of dust
acoustic waves to an external periodic driver in a model system described by
the forced Korteweg-de Vries-Burgers equation. This equation provides a proper
theoretical framework and a better physical model compared to the van der Pol
oscillator model for studying the dynamics of nonlinear dust acoustic waves by
properly accounting for nonlinear, dispersive and dissipative influences on
the waves. Using the model, we have successfully demonstrated harmonic (1:1)
and super-harmonic (1:2) synchronization states of DAWs. Our numerical
findings show a strong qualitative resemblance to the experimental results on
synchronization reported by Ruhunusiri et al. (Ruhunusiri and Goree, 2012). In
particular, comparison of our theoretical Arnold tongue diagram with their
experimental one shows the following common features. As in the experimental
Arnold tongue diagram we see the existence of amplitude thresholds as well as
clear evidence of the branching phenomena. However there are also important
differences. With our model we have not been able to obtain sub-harmonic
synchronization that have been observed in the experiment. Furthermore, our
model uses an external driver that closely resembles a nonlinear natural mode
of the system whereas in the experiment a purely time varying external
sinusoidal driver has been used. However, it is not clear what form this
driver takes inside the plasma system and whether it manifests itself as a
spatio-temporally varying perturbation. These and other questions, such as the
absence of sub-harmonic synchronization in the equation, the neglect of
dissipation arising from gas friction on the dust particles, etc., remain to
be explored in order to further improve the model.

###### Acknowledgements.

Work done by ST and AM was supported by the Indian Institute of Technology
Jammu Seed Grant No. SG0012. AS is thankful to the Indian National Science
Academy (INSA) for the INSA Honorary Scientist position. CC and GG acknowledge
NASA-JPL subcontract No. 1573108 and NRL Base Funds. JG was supported by
United States Department of Energy Grant No. DE-SC0014566, NASA/JPL RSA No.
1672641, and National Science Foundation Grant No. PHY-1740379.

## References

  * Pikovsky  _et al._ (2001) A. Pikovsky, M. Rosenblum, and J. Kurths,  _Synchronization: A Universal Concept in Nonlinear Sciences_ (Cambridge University Press, Cambridge, 2001).
  * Strogatz (2004) S. Strogatz,  _Sync: The Emerging Science of Spontaneous Order_ (Penguin, UK, 2004).
  * Balanov  _et al._ (2009) A. Balanov, N. Janson, D. Postnov, and O. Sosnovtseva,  _Synchronization: From Simple to Complex_ (Springer-Verlag, Berlin, 2009).
  * Acebrón  _et al._ (2005) J. A. Acebrón, L. L. Bonilla, C. J. Pérez Vicente, F. Ritort, and R. Spigler, Rev. Mod. Phys. 77, 137 (2005).
  * Kuramoto (1984) Y. Kuramoto, Progr. Theoret. Phys. Suppl. 79, 223 (1984).
  * Aronson  _et al._ (1990) D. Aronson, G. Ermentrout, and N. Kopell, Physica D: Nonlinear Phenomena 41, 403 (1990).
  * Kuramoto (2003) Y. Kuramoto,  _Chemical Oscillations, Waves, and Turbulence_ (Dover, New York, 2003).
  * Williams and Duff (2010) J. D. Williams and J. Duff, Phys. Plasmas 17, 033702 (2010).
  * Block  _et al._ (2001) D. Block, A. Piel, C. Schröder, and T. Klinger, Phys. Rev. E 63, 056401 (2001).
  * Klinger  _et al._ (1993) T. Klinger, A. Piel, F. Seddighi, and C. Wilke, Phys. Lett. A 182, 312 (1993).
  * Klinger  _et al._ (1995) T. Klinger, F. Greiner, A. Rohde, A. Piel, and M. E. Koepke, Phys. Rev. E 52, 4316 (1995).
  * Gyergyek (1999) T. Gyergyek, Plasma Phys. Controlled Fusion 41, 175 (1999).
  * Nurujjaman and Iyengar (2010) M. Nurujjaman and A. N. S. Iyengar, Phys. Rev. E 82, 056210 (2010).
  * Ruhunusiri and Goree (2012) W. D. S. Ruhunusiri and J. Goree, Phys. Rev. E 85, 046401 (2012).
  * Pilch  _et al._ (2009) I. Pilch, T. Reichstein, and A. Piel, Phys. Plasmas 16, 123709 (2009).
  * Williams (2014) J. D. Williams, Phys. Rev. E 90, 043103 (2014).
  * Chaubey  _et al._ (2015) N. Chaubey, S. Mukherjee, A. N. Sekar Iyengar, and A. Sen, Phys. Plasmas 22, 022312 (2015).
  * Fukuyama  _et al._ (2006) T. Fukuyama, R. Kozakov, H. Testrich, and C. Wilke, Phys. Rev. Lett. 96, 024101 (2006).
  * Chaubey  _et al._ (2016) N. Chaubey, S. Mukherjee, A. Sen, and A. N. S. Iyengar, Phys. Rev. E 94, 061201 (2016).
  * Shukla and Eliasson (2009) P. K. Shukla and B. Eliasson, Rev. Mod. Phys. 81, 25 (2009).
  * Morfill and Ivlev (2009) G. E. Morfill and A. V. Ivlev, Rev. Mod. Phys. 81, 1353 (2009).
  * Shukla and Mamun (2001) P. K. Shukla and A. A. Mamun,  _Introduction to Dusty Plasma Physics_ (Institute of Physics, Bristol, 2001).
  * Thompson  _et al._ (1997) C. Thompson, A. Barkan, N. D’Angelo, and R. L. Merlino, Phys. Plasmas 4, 2331 (1997).
  * Schwabe  _et al._ (2007) M. Schwabe, M. Rubin-Zuzic, S. Zhdanov, H. M. Thomas, and G. E. Morfill, Phys. Rev. Lett. 99, 095002 (2007).
  * Barkan  _et al._ (1995) A. Barkan, R. L. Merlino, and N. D’Angelo, Phys. Plasmas 2, 3563 (1995).
  * Heinrich  _et al._ (2009) J. Heinrich, S.-H. Kim, and R. L. Merlino, Phys. Rev. Lett. 103, 115002 (2009).
  * Bandyopadhyay  _et al._ (2008) P. Bandyopadhyay, G. Prasad, A. Sen, and P. K. Kaw, Phys. Rev. Lett. 101, 065006 (2008).
  * Deka  _et al._ (2017) T. Deka, A. Boruah, S. K. Sharma, and H. Bailung, Phys. Plasmas 24, 093706 (2017).
  * Piel and Melzer (2001) A. Piel and A. Melzer, Plasma Phys. Controlled Fusion 44, 1 (2001).
  * Liu  _et al._ (2018) B. Liu, J. Goree, T. M. Flanagan, A. Sen, S. K. Tiwari, G. Ganguli, and C. Crabtree, Phys. Plasmas 25, 113701 (2018).
  * Thompson  _et al._ (1999) C. Thompson, A. Barkan, R. Merlino, and N. D’Angelo, IEEE Trans. Plasma Sci. 27, 146 (1999).
  * Don  _et al._ (2014) W. Don, S. Ruhunusiri, and J. Goree, IEEE Trans. Plasma Sci. 42, 2688 (2014).
  * Feng  _et al._ (2007) Y. Feng, J. Goree, and B. Liu, Rev. Sci. Instrum. 78, 053704 (2007).
  * Deka  _et al._ (2020) T. Deka, B. Chutia, Y. Bailung, S. K. Sharma, and H. Bailung, Plasma Sci. Technol. 22, 045002 (2020).
  * Williams (2018) J. Williams, IEEE Trans. Plasma Sci. 46, 806 (2018).
  * Liu  _et al._ (2021) B. Liu, J. Goree, S. Schütt, A. Melzer, M. Y. Pustylnik, H. M. Thomas, V. E. Fortov, A. M. Lipaev, A. D. Usachev, O. F. Petrov, A. V. Zobnin, and M. H. Thoma, IEEE Trans. Plasma Sci. 49, 3958 (2021).
  * Van Der Pol (1927) B. Van Der Pol, Philos. Mag. 3, 65 (1927).
  * Kovalev and Kovalev (2017) D. P. Kovalev and P. D. Kovalev, Int. J. Bifurc. Chaos 27, 1750195 (2017).
  * Sen  _et al._ (2015) A. Sen, S. Tiwari, S. Mishra, and P. Kaw, Adv. Space Res. 56, 429 (2015).
  * Kumar Tiwari and Sen (2016) S. Kumar Tiwari and A. Sen, Phys. Plasmas 23, 022301 (2016).
  * Jaiswal  _et al._ (2016a) S. Jaiswal, P. Bandyopadhyay, and A. Sen, Phys. Rev. E 93, 041201 (2016a).
  * Nakamura  _et al._ (1999) Y. Nakamura, H. Bailung, and P. K. Shukla, Phys. Rev. Lett. 83, 1602 (1999).
  * Nakamura and Sarma (2001) Y. Nakamura and A. Sarma, Phys. Plasmas 8, 3921 (2001).
  * Veeresha  _et al._ (2010) B. M. Veeresha, S. K. Tiwari, A. Sen, P. K. Kaw, and A. Das, Phys. Rev. E 81, 036407 (2010).
  * Jaiswal  _et al._ (2016b) S. Jaiswal, P. Bandyopadhyay, and A. Sen, Phys. Plasmas 23, 083701 (2016b).
  * Orłowski (1994) A. Orłowski, Phys. Rev. E 49, 2465 (1994).
  * Rech (2013) P. C. Rech, Eur. Phys. J. B 86, 1 (2013).
  * Boyd (2003) J. P. Boyd,  _Chebyshev and Fourier Spectral Methods_ (Dover, New York, 2003).
  * Mir  _et al._ (2022) A. Mir, S. Tiwari, and A. Sen, Phys. Plasmas 29, 032303 (2022).
  * Mir  _et al._ (2020) A. A. Mir, S. K. Tiwari, J. Goree, A. Sen, C. Crabtree, and G. Ganguli, Phys. Plasmas 27, 113701 (2020).
  * Jensen (2002) R. V. Jensen, Am. J. Phys. 70, 607 (2002).