# Dynamics of kink train solutions in deformed Multiple sine-Gordon models

Marzieh Peyravi ID<EMAIL_ADDRESS>Department of Physics,
School of Sciences, Ferdowsi University of Mashhad, Mashhad 91775-1436, Iran
Nematollah Riazi<EMAIL_ADDRESS>Physics Department, Shahid Beheshti
University, Evin, Tehran 19839, Iran  Kurosh Javidan<EMAIL_ADDRESS>Department of Physics, School of Sciences, Ferdowsi University of Mashhad,
Mashhad 91775-1436, Iran

###### Abstract

This paper examines the effects of a thin layer of inhomogeneity on periodic
solutions of the Multiple-sine-Gordon (MsG) model. We investigate the dynamics
of the perturbed Double-sine-Gordon (DsG) system as a significant and more
practical case of such configurations. The thin barrier acts as a potential
well (potential barrier) and causes critical deformations in kink train
solutions and some basic properties of the periodic solutions, such as the
type of sub-kinks, their amplitude, energy and wavelength. Stability of the
initial kink chain during the interaction with medium defects is analyzed
using their phase diagram. Sudden changes in the profile of kink trains due to
the disruption of their amplitude and wavelength are considered. The time
evolution of moving kink chain solutions while interacting with medium
fractures is also studied.

keywords: Double-sine-Gordon, Multiple-sine-Gordon, Soliton collision,
Integrability.

##  I Introduction

The formation, stabilization, and propagation of localized waves in non-linear
media such as solitary waves, solitons, and breathers are fundamental problems
in many branches of science. Scattering of solitons from point-like int1 ;
int2 ; int3  or extended int4 ; int5  defects and impurities as a fundamental
problem of soliton theory has been widely investigated. Indeed, such effects
play a crucial role in determining the essential physical properties of non-
linear systems, and this issue has attracted many theoretical and practical
interests int6 ; int7 ; int8 . The processes of soliton trains in quantum
optics soltrain1 ; soltrain2 , a chain of localized solutions in Josephson-
Junction Arrays soltrain3 ; fx2 , and the dynamics of soliton ratchets in
perturbed background medium 12 ; 13 ; 14  are important topics, both in theory
and also applications. It is clear from an experimental point of view, we have
to consider spatial defects to find a realistic perspective of the system.
Inhomogeneity can arise in a medium due to various defects, such as
dislocations, impurities, imperfect grain boundaries, etc. Thus, one can see
defects as local deviations from an ordered structure in the system int4 ;
int9 ; int10 .

The family of Sine-Gordon (SG) models including its modifications, have
appeared in different branches of physics and engineering. There are a lot of
investigations that focus on multi-sine-Gordon (MsG) models and especially the
double-sine-Gordon (DsG) equation. The DsG model has been used to explain the
non-linear nature of many physical phenomena like the spin dynamics in the B
phase of super-fluid 3He 2 ; 3 , charge density wave condensate models of the
organic linear conductors, and dynamics of the liquid crystals org ; lq1 ; lq2
, quark confinement 10 , self-induced transparency phenomena 11  and many more
situations. The non-linear vibrations arising in the engineering applications
by considering the noise and uncertain properties en1 , or propagation of
ultra-short optical soliton trains tr  are modelled by different MsG
equations. Family of sG potentials appear in different ways to describe
natural non-linear phenomena. The general form of these potentials is in the
form $V(\phi)=a+b\cos(n\phi)+c\cos(m\phi)$ , where $a,b,c,n$ and $m$ are
system parameters. In practical situations (such as soliton/solitary wave-
based telecommunication, applications related to crystal structures and …) one
needs soliton train solutions. However non-linear properties of multi-
frequency sG potentials have been studied extensively, unfortunately, there
exist analytical periodic solutions only for very special cases of such
potentials. For this reason, real models are often approximated to such known
forms. Due to the non-linear nature of these equations, the applied
approximations cannot support all the results of the real phenomenon. In
particular, the stability of solutions and their sensitivity respect to
variation of model parameters (especially for periodic solutions) in reduced
models, have not been investigated. Thus, here we solve the general form of
the equations numerically. To ensure the correctness of outcomes, we compare
our results (in certain cases) where there are previously published analytical
periodic solutions. According to our best knowledge, the periodic solution has
been found only for the single and double sine-Gordon equations.
Unfortunately, such solutions are in series expansion forms. The breather
train solution 1j  and periodic solution 2j ; 3j  for the standard sG have
been presented analytically. The general behaviour of multi sine-Gordon model
has been investigated in few papers. One of the best studies is presented in
4j . The periodic solution for stationary multi sine-Gordon has been
investigated numerically in (paper2, ) by some authors of this manuscript.

In real environments, medium defects in the form of disorders, impurities, or
dislocations cause serious changes in the character of propagated waves. The
stability of travelling waves, especially in non-linear media is not a simple
and straightforward issue. In such problems, we need a phase diagram analysis
as well as the frequency–amplitude relationship under the existence of medium
inhomogeneities to understand the actual vibration properties of the medium.
Investigation of the dynamics of these types of phenomena needs periodic
solutions of MsG equations (mostly the DsG model) in the presence of medium
defects and impurities. This means that we have to deal with two problems at
the same time: finding periodic solutions for MsG equations and also adding
medium defects to the model by considering suitable potentials and/or
modifications. Motivated by this question, we have studied the interaction of
periodic solutions with medium defects.

In this contribution, we consider solutions of the DsG and MsG equation in
$1+1$ space-time dimensions. Our goal is to study the effect of some kinds of
barriers, somewhere among the chains of the kinks and anti-kinks of the DsG
and MsG systems. The amplitude and wavelength of solutions in the presence of
different perturbations will be studied. We thus propose to study the modified
periodic DsG and MsG solutions.

The structure of our presentation is as follows: Section II belongs to review
some general properties of the DsG and MsG system, including the action, the
field equation, and the single kink solution. In Section III we study the
effects of thin barriers on the periodic solutions of the DsG and MsG systems.
In Section IV, we focus on the collision of the DsG with the wall. We close in
Section V by summarizing our results and pointing out some directions for
future work.

##  II Basic Properties of the Double sine-Gordon and the Multiple-sine-Gordon
Systems

In this section, we briefly review the DsG and MsG formalism by introducing
the related action, the field equation, the energy-momentum tensor, and the
topological current.

Following our previous study on the periodic and step-like solutions of DsG
equation and static properties of the MsG systems paper2 , we focus on the
following potential (see Fig. 1):

$V(\phi)=1+\epsilon-\cos\phi-\epsilon\cos(N\phi)$ (1)

where $\epsilon$ is a non-negative constant and $N$ is an integer with $N=0$,
$N=2$, and $N>2$ corresponding to the sG, DsG, and MsG systems, respectively
paper2 ; ri1 . The harmonic term in this potential can result from the Fourier
expansion of an arbitrary, periodic potential
$\widetilde{V}(\phi)=\widetilde{V}(\phi+2n\pi)$.

Figure 1: (a) The DsG Potential. The dashed curve is for $\epsilon=10$, the
dash-dotted curve is for $\epsilon=1$, and the solid curve is for $\epsilon=0$
(sG). (b) The MsG Potential for $N=5$. The dashed curve is for $\epsilon=10$,
the dotted curve is for $\epsilon=1$ and the solid curve is for $\epsilon=0$ .

As the Fig. 1 illustrates, the above potential find its minima at:
$\phi_{min}=2n\pi$ with $n=0,1,2,...$. Based on the nature of the potential,
more vacua will appear in the profile of the potential, as $N$ increases, and
thus, extra sub-kinks (up to $N$) are created in the solution, which leads to
a more complicated system paper2 .

While integrable equations are an essential part of the non-linear wave
theory, most non-linear wave equations that emerge in physics and engineering
are not integrable majhool . Non-integrable equations involve much richer and
more complex dynamical solutions. Non-integrable equations not only lead to
unstable solitary waves with complicated and even chaotic interactions via
fractal scattering phenomena but also these solitary waves can collapse in
more than one dimension majhool .

Several powerful methods have been presented for deriving localized travelling
wave solutions for non-linear equations in recent years, for example, the
inverse scattering transformation, Hirota bilinear method, Backlünd and
Darboux transformations, sine-cosine and tanh-function methods, homogeneous
balance and Lie group analysis, F-expansion method and so on. Through these
methods, solitary, periodic, multiply periodic, quasi- or non-periodic wave
solutions as well as kink (anti-kink) trains have been constructed as
polynomials in the triangle, hyperbolic, or Jacobian elliptic functions for
the DsG model (See Fig. 2) j.amc ; DsGJaccobi ; DsG periodic . Several kink
train solutions for the DsG model have been proposed using Jacobi’s elliptic
functions, such as:

$\phi=\pi+2\arctan\left[(k^{\prime}_{2})^{-1}JN\bigg{(}\frac{2k^{\prime}_{2}\sqrt{\epsilon}}{k_{1}\sqrt{-k^{2}_{2}}}(x-x_{0}-vt),k\bigg{)}\right].$
(2)

where, $JN(u,k)$ is a suitable generalized Jacobi elliptic function, related
to the potential parameters,
$k^{2}_{1,2}=\frac{1}{3}\left[(1-4\epsilon)\pm\sqrt{(1+\epsilon)^{2}+8\epsilon}\right]$,
$k^{2}=\frac{k_{1}^{2}-k_{2}^{2}}{1-k_{2}^{2}}$, and
$k^{\prime}_{2}=\sqrt{1-k_{2}^{2}}$ DsGJa . The Jacobi elliptic functions take
the parameters $u$ and $k$ as inputs. In the following integration:

$u=\int^{\phi}_{0}\frac{d\theta}{\sqrt{1-ksin^{2}\theta}}$ (3)

the parameter $\phi$ is called the Jacobi amplitude $A(u,k)=\phi$. In this
framework, the elliptic sine is introduced as $SN(u,k)=\sin(A(u,k))$ and the
elliptic cosine is given by $CN(u,k)=\cos(A(u,k))$. In similar definition we
can define the elliptic delta amplitude as
$DN(u,k)=\frac{d}{du}\left(A(u,k)\right)$. Figure 2 demonstrate These three
Jacobi elliptic functions respect to the parameter $u$ with $k=0.99$.

Explicit solutions of non-linear equations are very helpful to analyse all the
dynamical features of the system. Unfortunately, only periodic solutions of sG
and DsG potentials have been obtained using different approaches, which are
often not explicit 1jj ; 2jj . The periodic solutions of the DsG potentials
are presented as expansions of Jacobi functions, which can be converted into
closed form only in certain values of model parameters, which are often not
physical 3jj . Explicit solutions for the MsG equation have not been found
yet. Thus, we can only rely on numerical calculations. To control the
accuracy, the results of numerical calculations on the MsG potential is
reduced to the DsG model and then compare them with analytical solutions of
the DsG equation.

According to our best knowledge, there are no such solutions for the general
form of the MsG system. Therefore we calculated periodic solutions for the MsG
potential numerically. Indeed, most of the presented solutions for more
complicated potentials, have been given in infinite series expansions. In such
situations, we have to use numerical calculations/simulations to study the
behaviour of the solutions. Nowadays, numerical computations play an important
and inevitable role in non-linear science and soliton theory majhool . Since
non-integrable systems are not solvable analytically, numerical methods
provide powerful tools for non-linear wave studies.

Figure 2: Jacobi Elliptic Functions as a function of u for $k=0.99$.

The Lagrangian density of the system consisting of the MsG system for real
scalar field $\phi(x,t)$ in $(1+1)$ dimensions is in the following form paper2
:

${\cal
L}=\frac{1}{2}\partial^{\mu}\phi\partial_{\mu}\phi-\left[1+\epsilon-\cos\phi-\epsilon\cos(N\phi)\right].$
(4)

From this Lagrangian density, the field equation of motion is as follows
paper2 :

$\Box\phi=-\sin\phi-N\epsilon\sin(N\phi).$ (5)

The energy-momentum tensor of the MsG equation can be obtained by using
Noether’s theorem and the invariance of the action under the space-time
translation $x^{\mu}\rightarrow^{\mu}+a^{\mu}$ paper2 ; 22 ; 23 :

$T^{\mu\nu}=\partial^{\mu}\phi\partial^{\nu}\phi-g^{\mu\nu}{\cal L},$ (6)

in which $g^{\mu\nu}=diag(1,-1)$ is the metric of the Minkowski space-time.
Besides, a topological current can also be defined for a generalized form of
sG through the following current 24 :

$J^{\mu}=\frac{1}{2\pi}\epsilon^{\mu\nu}\partial_{\nu}\phi,$ (7)

where $\epsilon^{\mu\nu}$ is the antisymmetric tensor with $\epsilon^{01}=1$.
This current is identically conserved ($\partial_{\mu}J^{\mu}=0$) and the
total topological charge of any localized, finite-energy solution is conserved
and also quantized. We have used the topological charge and energy of the
system for controlling the validity of our numerical calculations in every
selected step.

As one can check, the MsG equation possesses kink and anti-kink solutions
which correspond to transitions between the spatial boundary conditions
$\phi(\pm\infty)=2n\pi$ paper2 . The first integral of the equation (5) for
static, localized solutions reads:

$\frac{1}{2}(\frac{d\phi}{dx})^{2}=V(\phi),$ (8)

in which we have used the boundary conditions $\phi(\pm\infty)=2n\pi$. We
therefore have:

$x-x_{0}=\int\frac{d\phi}{\sqrt{2V(\phi)}}.$ (9)

Multiple sine-Gordon systems admit soliton-like kink solutions with
interesting properties, although this integral cannot be carried out
analytically for general $N$ and $\epsilon$ specially for $N>2$. However for
$N=2$ one can show that the following exact, static, single kink (anti-kink)
solution exists ri1 ; paper2 :

$\displaystyle\phi(x)=2\arccos\left[\pm\frac{\sinh\sqrt{4\epsilon+1}x}{\sqrt{4\epsilon+\cosh^{2}\sqrt{4\epsilon+1}x}}\right].$
(10)

The $(+)$ and $(-)$ sign are correspond to the anti-kink ($\tilde{k}$) and the
kink ($k$) solutions respectively ri1 . Figs. 3 illustrate DsG kink/anti-kink
Eq. (10) for two values of the parameter $\epsilon(=6,10)$ as well as
solutions of $N=3$ which is obtained numerically. We can interpret and adjust
our numerical solution for the kink train in DsG by comparing it with the
exact solution (2).

Figure 3:  Single soliton solutions (a) For DsG with $\epsilon=10$ (dash-
dotted curve) and $\epsilon=6$ (dashed curve) (b) The MsG system ($N=3$) with
$\epsilon=10$ (dash-dotted curve) and $\epsilon=6$ (dashed curve) .

It is easy to show that solutions of the MsG equation are classified and
characterized by (x-independent $\frac{dP}{dx}=0$) constant of integration
$P$:

$P=\frac{1}{2}\left(\frac{d\phi}{dx}\right)^{2}-V(\phi).$ (11)

in two types, namely: Step-like and periodic solutions (See Fig. 4). Note that
this classification is independent of $\epsilon$ and therefore holds true also
for sG equation. The step-like solutions are a sequence of kinks ($kk$) or
anti-kinks ($\tilde{k}\tilde{k}$) which are characterized by $P>0$, and the
periodic solutions are a sequence of kink, anti-kink ($k\tilde{k}k\tilde{k}$)
solutions that are characterized for $-2<P<0$ paper2 . The general behaviour
of our numerical solutions is in agreement with reported solutions in Ref. baz
. Recently, a set of new kink solutions for special types of MsG models have
been presented, using the mentioned Jacobi elliptic functions joseph . In such
cases, our static solutions are successfully reduced to kink solutions in the
Ref. joseph . Also, there are localized solutions for the triple sine-Gordon
equations ($N=3$) for specific values of the model parameters TsG .

In fact, one can interpret $P$ as the “pressure” (or tension in this 1d case)
when we consider such systems as a many-body interacting system
Moradi.EPJB.2018 ; Gani.EPJC.2018 ; Gani.EPJC.2019 . Note that $P$ is
different from energy density

$\rho(x)=\frac{1}{2}\left(\frac{d\phi}{dx}\right)^{2}+V(\phi);$ (12)

which changes with position $x$ paper2 ; Moradi.EPJB.2018 . However MsG
potential contains infinite topological sectors, but the multiplicity of
solutions is equal to $N$ which are different in terms of the number of their
sub-kinks. baz . In addition to soliton train solutions, the model contains
breather solutions for $N^{2}<4\pi$ bth .

##  III The effect of thin barriers

In this section, we report the behaviour of kink train solutions of the DsG
equation in the presence of a thin interval of inhomogeneity. For this purpose
$\epsilon$ is considered to be different within a finite range $x_{1}<x<x_{2}$
in which $x_{1}$ and $x_{2}$ are known locations while $x_{2}>x_{1}$.

As stated before, the behaviour of the solution depends on the value of
parameter $P$, which is a function of field tension (stress)
$\frac{d\phi}{dx}$ and also the potential ($U(\phi)$). This means that we can
control the field solution by correctly choosing the value of the field
($\phi(x)$) and its spatial derivative ($\frac{d\phi}{dx}$) as initial
conditions. It is clear that medium dislocations and impurities change these
conditions, and therefore such disorders critically affect the field solution
in the medium.
The area of perturbation, where the barrier-field interaction starts
($x_{1}$), and terminates ($x_{2}$), barrier width ($x_{2}-x_{1}$), and its
amplitude ($\epsilon$) are important parameters. These parameters determine
the general characteristics of the solution in the perturbed area as well as
the kink train configuration after passing through the barrier.

The minimum and maximum value of parameters depend on the possible values of
$\frac{d\phi}{dx}$ and $V(\phi)$ as one can find from Eq. (11). These values
generally depend on the condition of the field at the starting point of the
perturbation area and the characteristics of the barrier. The constraints
imposed at the border points of the barrier can limit the possible values
after the kink train leaves the perturbation region. The set of initial
conditions for Fig. 4(b) causes the response to appear as a kink train in the
perturbation region. The possible values for the $\frac{d\phi}{dx}$ are always
greater than the value of the potential at the barrier $V(\phi(x))$ (see Fig.
5). Indeed, by changing the length of the perturbation region, it is not
possible to decrease the value of P to find anti-kink train solution after the
interaction. It is clear that by changing the characteristics of the barrier
as well as initial kink train solution, other different train solutions may
occur.

Kink and anti-kink solutions Eq. (10) are opposite to each other in terms of
spatial expansion. Thus, the field of a kink train ($kk$) solution is an
increasing function with respect to the space. The field solution of anti-kink
train ($\tilde{k}\tilde{k}$) is a spatial decreasing function, but the
solutions in the form of kink anti-kink trains ($k\tilde{k}$) are periodic
functions of the space coordinate. One can find all these configurations in
Fig. 4. The initial configuration at $x<x_{1}$ of the Fig. 4(a) is a
$k\tilde{k}$ solution with a periodic variation. The dotted line curve in the
Fig. 4(a) demonstrates a ($\tilde{k}\tilde{k}$) configuration, while the solid
line shows a ($kk$) kink train solution.

Figure 4: (a) $\phi$ diagram as a function of $x$ for $N=2$ and
$P=-0.009987500000000$; $\epsilon=10$ for $0\leq x\leq x_{1}$, $x_{2}\leq
x\leq 100$ and $\epsilon=5$ for $x_{1}<x<x_{2}$. For dotted, dash-dotted and
solid curves $x_{1}=26.74,28.01,29.05$ and $x_{2}=40$ respectively. (b) $\phi$
diagram as a function of $x$ for $N=2$ and $P=-0.009987500000000$;
$\epsilon=10$ for $0\leq x\leq x_{1}$, $x_{2}\leq x\leq 100$ and $\epsilon=5$
for $x_{1}<x<x_{2}$. For dotted, dash-dotted, dashed and solid curves
$x_{1}=30$ and $x_{2}$ are corresponded to $35,40,45$ and $50$ respectively.

Figure 4(a) shows critical changes related to the starting point of
perturbation ($x_{1}$). By a slightly changes of $x_{1}$, we can provide step-
like kink ($kk$) (anti-kink $\tilde{k}\tilde{k}$) trains or periodic
$k\tilde{k}$ configurations, in the perturbation region. It is due to changing
the initial condition for creating the kink solution in the perturbed region.

The most important point is related to the slope of arriving kink (anti-kink)
as interacted with the perturbation. Fig. 4(a), clearly shows that the kink
(anti-kink) train is created in the perturbation region if the initial kink
arrives with a positive (negative) slope. If the initial kink arrives with a
very small slope, we will find $k\tilde{k}$ configuration in the perturbed
region.

Figure 4(b) demonstrates the interaction of kink solution with perturbations
starting from the same initial position (here $x_{1}=30$ with $\epsilon=5$)
but with different widths. Emerging configuration from the perturbed region
depends on the characters of the solution at the end point of the perturbation
region ($x_{2}$). In general, the final state of the solution after
interacting with the barrier can be any of the possible solutions. Some kink
configurations may not occur depending on the possible values for the slope of
the field in the barrier. In Figure 4(b), we do not expect the formation of
the state ($\tilde{k}\tilde{k}$), because the $kk$ configuration which has
been established in the perturbation region, does not have a negative slope.
Thus, the final state of the chain configuration after passing the barrier can
be periodic ($k\tilde{k}$ ) or step-like ($kk$). We may observe the emergence
of $k\tilde{k}$ configuration if the kink reaches the end point of the
perturbation area with its flat tail (Please see final $k\tilde{k}$
configurations in the Fig. 4(b) for $x_{2}=35,40$ and $45$). On the other
hand, a kink train configuration may occur after the interaction, if the field
solution reaches the barrier endpoint with its positive slope parts.
Similarly, if the established configuration in the perturbed area does not
have a positive slope, the $kk$ train will not occur after the field-impurity
interaction.

Our other simulations are in agreement with this procedure. If the anti-kink
configuration is created within the barrier region, we can reach $k\tilde{k}$
or $\tilde{k}\tilde{k}$ trains. For a $k\tilde{k}$ configuration in the
barrier, we expect to find $k\tilde{k}$, $kk$ or $\tilde{k}\tilde{k}$ trains
according to the final state of solution where the soliton-barrier interaction
is completed. The potential barrier may be an impurity/spatial dislocation or
a defect in the crystal lattice, where the kink train propagates in real
physical situation. According to the Figure 4, due to a sudden change in the
texture of the background space, the characteristic of the kink train
experiences critical changes. As Figure 4 shows, the interaction of the kink
train with a structural fault causes critical changes in the topological
characteristics of the solution. Of course, we have to pay attention on the
dynamical behaviour of the system, such as stability of solution, kink train
shape, its velocity after collision. To understand such these characters, we
study the description of the barrier-field interaction in the phase space as
shown in the Fig. 5. The quantity $\frac{d\phi}{dx}$ actually refers to the
field stress through the term $\frac{1}{2}(\frac{d\phi}{dx})^{2}$ which
expresses the amount of exerted force by the field (per unit length, in $1+1$
dimensions).

Figure 5:  The phase diagram for $N=2$ and $P=-0.009987500000000$;
$\epsilon=10$ for $0\leq x\leq 30$, $x_{2}\leq x\leq 100$ and $\epsilon=5$ for
$x_{1}<x<x_{2}$. Dotted, dash-dotted, dashed and solid curves $x_{2}$ are
corresponded to $35,40,45$ and $50$ respectively.

Figure 5 demonstrates the phase space diagram, for configurations of Fig.
4(b). As one can see, for the shortest barrier ($30<x<35$, the dotted curve)
phase diagram has only one ring due to the absence of sub-kinks. For different
impurity widths ($30<x<40$ and $30<x<45$, dash-dotted, dashed curves) we find
other diagrams containing twofold closed loops due to the presence of sub-
kinks. It should be noted that these diagrams relate to the same type of
solution but in different sectors. The noticeable point in the Fig. 5, relates
to the solid open curve which corresponds to the $kk$ train with sub-kinks. We
learn from Fig. 5 that, the footprint of each kink (anti-kink) is a single or
double concave (convex)-like curve in the phase space diagram. $k\tilde{k}$
configurations create closed loops in a limited area of the phase space plane
if both kink and anti-kink are located in the same sector. But $kk$ and
$\tilde{k}\tilde{k}$ configurations are not closed curves. This means that,
each kink (anti-kink) in created $kk$ or $\tilde{k}\tilde{k}$ trains are not
located in the same sector. However the multiplicity of solutions in this
model is $N=2$, but topological charge of all kinks (anti-kinks) are the same.
Thus, all created objects in a train solution, show similar behaviour baz .

Dynamic behaviours depending on open and closed paths in the phase space can
be explained based on equation (12). The energy of a kink chain configuration
can be obtained by calculating potential energy $V(\phi)$ and stress
$\frac{1}{2}(\frac{d\phi}{dx})^{2}$ . Figure 1 shows that function $V(\phi)$
is a periodic function, and always changes in a limited range. In other words,
this quantity never tends to infinity. But the $\frac{d\phi}{dx}$ (at least
theoretically) has no limit. Phase space ($\phi-\frac{d\phi}{dx}$) clearly
demonstrates the variation of $\frac{d\phi}{dx}$ in terms of the field
solution. Phase space paths in figures 5-8 demonstrate the changing range of
$\frac{d\phi}{dx}$ in terms of field value $(\phi)$. The phase space curves
represent the exchange of energy between its two different forms: potential
energy $V(\phi)$ and stress $\frac{1}{2}(\frac{d\phi}{dx})^{2}$. The closed
paths in the phase space correspond to the ($k\tilde{k}$) solutions, while the
open paths regard to ($kk$) or ($\tilde{k}\tilde{k}$) train solutions. The
phase space diagram of Fig. 6 clearly shows that the amplitude of
$\frac{d\phi}{dx}$ changes during the interaction with barrier. As a result,
the energy of the incoming solution changes due to collision with the
perturbation. But the energy of the kink train returns into its initial value
after interaction. This scenario does not occur in all interactions. For
example, Fig. 7 shows that the energy of the related kink train changes during
the interaction with the perturbation, but it is not restored to its initial
value after leaving the barrier.

For finding better knowledge about the similarity and differences between
created kinks (anti-kinks) during the field-impurity interaction, we have
plotted Figs. 6 and 7 with more details. In these figures, similar
$k\tilde{k}$ trains hit the barriers, which are located in the same spatial
position, but with different values of $\epsilon=5$ (Fig. 6) and $\epsilon=15$
(Fig. 7). Initial conditions for both situations are the same, but the related
barrier potential of defined situations is different. For barrier with
$\epsilon=5$ a $kk$ chain is created in the perturbed area. The solution for
barrier with $\epsilon=15$ is a $k\tilde{k}$ train with different energy but
in the same sector ($0\leq\phi\leq 2\pi$) as compared with initial
$k\tilde{k}$ (see Fig. 6). One can find from the phase space diagrams in Figs.
6 and 7 that, both solutions in the perturbation region have positive, zero,
and negative slope parts. Thus, all types of kink chain solutions can be
created after the interaction with defect ($x\geq x_{2}$), according to the
$\frac{d\phi}{dx}$ at the endpoint of perturbation ($x_{2}$).

These figures show that the established solutions for these two scenarios are
different. Indeed $kk$ train ($k\tilde{k}$ solution) is created in the
perturbation region for $\epsilon=5$ ($\epsilon=15$). According to the phase
space diagram of the Fig. 6 (Fig. 7), created $k\tilde{k}$ solution after
(during) the field-impurity interaction is located in the different (similar)
sector as compared with initial $k\tilde{k}$ train. An interesting issue is
the creation of kinks that are located in different spatial solutions but in
similar topological sectors. According to the Fig. 6 (Fig. 7), kinks which are
located in the regions $x_{1}\leq x\leq x_{2}$ and $x\geq x_{2}$ ($x\leq
x_{1}$) are related to the same sector $8\pi\leq\phi\leq 10\pi$
($0\leq\phi\leq 2\pi$). Figure 7 demonstrates the collision of the kink train
with a spatial defect (potential well) with lower density (smaller
$\epsilon$). It can be clearly found that the stress decreases in the thinner
medium. A very interesting point is that the amount of stress density returns
to its initial situation (with a good approximation) after the termination of
the field-impurity interaction as one can find from Fig. 6(b).

Moreover, it is vital to emphasize the role of $\epsilon$ value via the
appearance of potential well/barrier in the medium, in addition to the initial
position of perturbation. As the Figs. 4-7 demonstrate, the $\epsilon$ value
in the perturbed area has a decisive effect on the slope of $kk$ or
$\tilde{k}\tilde{k}$ train in the barrier region, which leads to the emergence
of different $k\tilde{k}$ chains after the interaction with the barrier. In
other words, the amplitude, and wavelength of the $k\tilde{k}$ chain will
change as a consequence of sub-kinks.

Figure 6: (a) The $\frac{d\phi}{dx}$ diagrams as a function of $x$ and (b) The
phase diagram $N=2$ and $P=-0.009987500000000$; $\epsilon=10$ for $0\leq x\leq
30$ and $40\leq x\leq 100$ and $\epsilon=5$ for $30<x<40$.

Figure 7: (a) The $\frac{d\phi}{dx}$ diagram as a function of $x$ and (b) The
phase diagram $N=2$ and $P=-0.009987500000000$; $\epsilon=10$ for $0\leq x\leq
30$ and $40\leq x\leq 100$ and $\epsilon=15$ for $30<x<40$.

Figure 8: (a) $\phi$ diagram as a function of $x$ and (b) The phase diagram
for $N=3$ and $P=-0.220000000000000$; $\epsilon=10$ for $0\leq x\leq x_{1}$,
$x_{2}\leq x\leq 100$ and $\epsilon=5$ for $x_{1}<x<x_{2}$. For dotted, dash-
dotted and solid curves $x_{1}=29.2,29.8,30$ and $x_{2}=40$ respectively.

In our simulations, $P$ is calculated numerically by solving the relevant
differential equation. Indeed, we need to know how the $P$ value changes
abruptly across the barrier. Here, we try to obtain an analytical expression
for the change in $P$ across the potential ($\bigtriangleup P$). We obtain our
results for time-independent static solutions, however, they can be easily
generalized to time-dependent solutions.

Consider two regions of different $\epsilon$ joining each other at $x=0$.
Integrating the static equation from $x=-\delta$ to $x=+\delta$

$\int^{+\delta}_{-\delta}\frac{d^{2}\phi}{dx^{2}}dx=\int^{+\delta}_{-\delta}\frac{\partial
V(\phi)}{\partial\phi}dx,$ (13)

and letting $\delta\longrightarrow 0$, we obtain

$\frac{d\phi}{dx}|_{0^{-}}=\frac{d\phi}{dx}|_{0^{+}}$ (14)

since there is no singularity in the potential around $x=0$. Together with the
continuity of $\phi$, we are now able to obtain the relation between the
$P$-values on either side of the junction:

$P_{1}=\frac{1}{2}(\frac{d\phi}{dx})^{2}-V_{1};\qquad\text{for $x<0$,}$ (15)

and

$P_{2}=\frac{1}{2}(\frac{d\phi}{dx})^{2}-V_{2};\qquad\text{for $x>0$,}$ (16)

using the continuity of $\phi$ and $\frac{d\phi}{dx}$ across the junction, we
obtain

$P_{2}-P_{1}=V_{1}-V_{2}=\left(\epsilon_{1}-\epsilon_{2}\right)\left(1-\cos(2\phi(0))\right);$
(17)

for the DsG equation, when $\phi(0)$ is the value of $\phi$ at $x=0$(junction
point). One observes that even for $\epsilon_{2}\neq\epsilon_{1}$, we can
still have $P_{2}=P_{1}$, i. e., if $\phi(x=0)=n\pi$. Otherwise, the value of
$P$ will change abruptly across the junction. The minimum and maximum change
in the value of $P$ is 0 and $2|\epsilon_{1}-\epsilon_{2}|$, respectively.
Similar relations can be presented for the MsG models too.

Figure 9: A single soliton across a boundary at $x=0$.

For a single soliton, we have $P=0$, if we consider a media with two different
$\epsilon$ on both sides of the boundary (on the left $\epsilon_{2}$ and the
right $\epsilon_{1}$) paper2 . Suppose that there is a single soliton on the
left side of the boundary ($P_{2}=0$), so the Eq.(17) reduces to:

$P_{1}=\left(\epsilon_{2}-\epsilon_{1}\right)\left(1-\cos(2\phi(0))\right);$
(18)

In this case, it is obvious that the vital condition for having the single
soliton on the right side of the boundary is $\epsilon_{2}=\epsilon_{1}$ or
$\phi_{0}=n\pi$ (See Fig. 9).

According to Eq. (17), the difference value of the parameter $P$ while
entering the perturbation area ($x=x_{1}^{+}$) becomes
$P_{1}=P_{2}-(\epsilon_{1}-\epsilon_{2})V(\phi(x_{1}))$ , while $\epsilon_{1}$
( $\epsilon_{2}$ ) relates to the $x>x_{1}$ ($x<x_{1}$). We have chosen
$\epsilon_{2}=10$ ($x<x_{1}$) and $\epsilon_{1}=5$ ( $x>x_{1}$ ) in the Fig.
4. Therefore, $\Delta P_{x_{1}}=+5V(\phi(x_{1}))$ (at starting point) and
$\Delta P_{x_{2}}=-5V(\phi(x_{2}))$ (while kink train leaves the barrier).
Since $V(\phi(x))$ is always positive, the value of parameter $P$ increases at
the starting point of the barrier, while it will decrease at the end of
perturbed region. To show the effect of $P$ in the configuration of kink train
while interacting with the barrier, we have chosen the starting point of the
boundary ($x_{1}$) to be close to the minimum value of $V(\phi(x))$ for
plotting the Fig. 4(a). In this case, for an initial value (in the Fig. 4(a)
$x_{1}=6.74$ ) $P$ remains negative after the kink train enters the barrier,
while for initial starting point ($x_{1}=29.05$ ) $\Delta P_{x_{1}}$ is enough
great to change the $P$ into a positive value after arriving at the barrier.
It is clear that for a location between these two values (where
$V(\phi(x_{1}))=0$) there is no change in $P$ , as one can found in the Fig.
4(a).

It should be noted that dislocations and impurities as spatial defects are not
the only sources of perturbations. Dissipative effects, dispersion, different
types of thermal noise and stochastic agents also strongly affect the
properties and dynamical behaviour of periodic solutions. It is shown that
thermally induced perturbations have definite impacts on the stability of some
localized solutions of sG models 1r . In such studies the spatially-uniform
Gaussian noisy force $\gamma(t)$ (with $<\gamma(t)>=0$) and delta-correlated
Gaussian thermal noise $\gamma_{T}(x,t)$ (with $<\gamma_{T}(x,t)>=0$ have been
considered 2r . Generation of travelling waves in non-linear models in the
presence of thermal noise is another important phenomenon related to the
environmental effects 3r ; 4r . The environmental perturbations as source of
space-time dependent noise and random fluctuations are critical components of
any non-linear and complex systems. Such defects play definitely important
roles in classical and quantum description of physical systems, from
microscopic to macroscopic scales 5r ; 6r ; 7r ; 8r ; 9r ; 10r ; 11r .

##  IV  soliton train-barrier interaction

The time evolution of soliton train solution during the interaction with
medium disorders and impurities is an important issue due to its practical
applications. The stability of the solution against medium defects can be
studied by finding the emerged object after the interaction of the incident
solution with a potential barrier.

We have performed several simulations to investigate the time evolution of a
moving kink train after interaction with perturbation. Figure. 10(a)
demonstrates the time evolution of the field solution while interacting with a
specified potential as a medium defect. The initial condition is the same as
what has been chosen for Fig. 7(a). From this figure, we expect to observe
moving $k\tilde{k}$ train interacting with the barrier. Therefore the initial
conditions at the beginning point of perturbation (and thus the emerged
solution from the barrier) are changing over time. As explained in section III
the emerging solution is $kk$, $k\tilde{k}$, and $\tilde{k}\tilde{k}$ trains
according to the initial condition. The positive slope (ascending part) of the
solution creates a $kk$ train of interaction, while the negative slope
(descending part) produces $\tilde{k}\tilde{k}$ configuration ( See Fig. 4).
It is clear that the flat parts of the emerging kink create $k\tilde{k}$
train. A complete set of initial conditions are created by positive, zero, and
negative slope parts. As one can find from the Fig. 10, a sequence of $kk$,
$k\tilde{k}$ and $\tilde{k}\tilde{k}$ trains are observed repeatedly, while
their amplitudes are changing in time.

The higher the kink slope as an initial condition, the faster the changes in
the soliton amplitude after the interaction. Faster changes (compact parts) in
Fig. 10(a) relate to a higher (positive or negative) slope of the initial
condition at the border of perturbation. The flat part of the initial
condition provides a complete set of $k\tilde{k}$ solutions, as we can see
within $1000<t<1500$, which depends on the velocity of the moving initial kink
solution.

Figure 10: (a) The $\phi$ diagram as a function of $t$ and (b) The phase
diagram, with initial values of Fig. 7.

Figure 10(b) demonstrates the time evolution of the kink solution in phase
space. As one can find, the phase space plot is formed by several times
repeating the closed loop curve in the phase space (Fig. 7(b)), of course with
some rapid changes that are seen as jumps due to field-barrier interaction.
Figure 10(a) shows that, the changes in the slope of $\phi$ are not uniform
with respect to time. Sometimes the $\frac{d\phi}{dx}$ changes are sharp and
sometimes slow. In figure 10(b), the numerical data sets of figure 10(a) have
been ploted. The areas with high density of numerical points are related to
the times when the slope of the field $\phi$ has more drastic changes. This
figure contains important information about the stability of moving kink
solutions while interacting with spatially limited perturbations. We have
examined several initial conditions. All our simulations have indicated that
limited amplitude barriers are not able to destroy the stability of soliton
trains; which is an important outcome.

Figure 11: (a) The $\phi$ diagram as a function of $t$ and (b) The phase
diagram, considering initial values of the solid curve in Fig. 8.

Figure 11 demonstrates field-potential interaction for $N=3$, the triple sG
model. The general behaviour of interaction is the same as what we have found
for the DsG model, except the generation of three sub-kinks in the profile of
the field solution. We have not found different behaviour or instability
condition for this model. Also, we have examined $N=4$ and $N=5$ solutions. No
new phenomena have been observed.

Time evolution of the field energy density $\rho(x,t)$ has been analysed
during the interaction with medium perturbations too. It has been shown
previously that, the total energy of single soliton changes during the
interaction with potential k1 ; k2 . Recovery of Soliton energy into its
initial quantity (before the interaction) depends on the nature of the soliton
model and indeed the impact of the barrier. There is not much information
about the time evolution of soliton train-barrier interactions. We have set up
several simulations to examine the time evolution of the energy density for
soliton chain solutions of $N=2$ and $N=3$ models while interacting with
simple rectangular potentials of different amplitudes and widths. As expected,
the total energy of the soliton train can be recovered with a good
approximation if the barrier amplitude is sufficiently low. The final profile
of the energy density after the interaction is a simultaneous function of the
barrier amplitude and the state of the soliton at the initial and final points
of the interaction. Thus, it is a complicated (and also challenging) study. In
general, the energy of solitons, their sub-kink contents, amplitude, and
widths are seriously affected by the interaction. Figure 7 presents the kink
train-barrier interaction, where the barrier has greater density (higher
$\epsilon$). The field stress density increases in the barrier, but it
decreases again after finishing the field-potential interaction, when the kink
train leaves the impurity area. The final value of the stress density after
the interaction strongly depends on the field initial conditions and the
characteristics of the barrier. The stability of the field energy and the
shape of the kink train during and after the interaction with the impurity is
understood by studying the time evolution of field potential energy and stress
density. Figure 10 (11) demonstrate interaction with the potential well
(barrier) which determines by lower (higher) value for $\epsilon$. An
important result that can be obtained from these figures is that the collision
of the kink train with environmental perturbations may distort the shape of
the periodic solution, but it does not have the ability to impose instability
on the system.

Figure 12 demonstrates three profiles of soliton-potential interaction at
different times for a moving train with initial conditions of Fig. 8. The
moving soliton collides with the potential, and as a result, the initial
conditions for creating new kinks train inside the inhomogeneity region change
with time. For this reason, the initial conditions also change at the endpoint
of the potential region.

Figure 12: Profiles of the energy density as a function of scaled time ($t$)
for $N=3$ and $P=-0.220000000000000$; $\epsilon=10$ out of the interaction
region and $\epsilon=5$ in the barrier; (a) $t=1$, (b) $t=1.2$ and (c) $t=1.8$
.

##  V concluding remarks

Impurities, dislocations, and other environmental inhomogeneities cause
decisive effects on the behaviour of localized and quasi-periodic solutions of
non-linear models. Such structural defects have been added to the model by
different methods. We have investigated the effect of such thin barriers on
the creation and evolution of multiple sine-Gordon kink trains by focusing on
the dynamics of the double sine-Gordon and the triple sine-Gordon models. The
analytical solution for soliton chain solutions (in series expansion) exists
only for the double Sine-Gordon model. For this reason, we calculated the kink
train solutions for multiple Sine-Gordon models using numerical calculations.
We have checked the validity of numerical calculations by comparing our
obtained results with analytical solutions of the double sine-Gordon model. We
investigated the final configuration of kink train solutions after the initial
kink chain-impurity interaction with details. The relationship between the
final solution (also the kink train created in the impurity region) and the
characteristics of the initial kink chain as well as the parameters of the
potential barrier has been analysed. Our simulations show that regardless of
the absolute location of the impurity and its width, the magnitude of the
field and its slope at the start and end points of the interaction (the
beginning and end of the perturbation region) are the most determining factors
in shaping the final kink train solution after passing through the potential.
The change in the value of model parameters ($\epsilon$ in our study) which is
determined by the nature of the impurity, determines the amplitude and
wavelength of final solutions after the interaction. The dynamics and general
behaviour of interaction and the conditions for establishing the kink train
solutions have been the same for all the investigated models ($N=2,...,N=5$).
The only difference is the appearance of sub-kinks, which originate from the
natural properties of different models.

The dynamics of soliton-potential interaction can not be fully understood
without considering the internal mode effects. Energy transfer between the
soliton (as a ground state) and internal mode (as excited states) clearly
changes the soliton’s initial condition while interacting with medium
impurities (as a potential barrier). Although the numerical solution
automatically simulates the effects of internal modes, it does not show the
effect of this phenomenon separately. For this reason, a suitable analytical
model is necessary to understand the effect of internal modes on the dynamics
of system evolution. Although we think that the effect is small, but this
issue should be investigated in further works.

###### Acknowledgements.

M. P. and K. J. acknowledge the support of Ferdowsi University of Mashhad. N.
R. acknowledges the support of Shahid Beheshti University Research Council.

### Data Availability Statement

This manuscript has no associated data or the data will not be deposited.
[Authors’ comment: This is a theoretical study and no experimental data has
been listed.]

## References

  * (1) D. I. Pushkarov and R. D. Atanasov, Phys. Lett. A 149, 287 (1990).
  * (2) X. D. Cao and B. A. Malomed, Phys. Lett. A 206, 177 (1995).
  * (3) A. D. Boardman, V. Bortolani, R. F. Wallis, K. Xie and H. M. Mehta, Phys. Rev. B 52, 12736 (1995).
  * (4) K. T. Stoychev, M. T. Primatarowa and R. S. Kamburova, Phys. Rev. E 70, 066622 (2004).
  * (5) M. T. Primatarowa and R. S. Kamburova, J. Phys. : Conf. Ser. 1762, 012014 (2021).
  * (6) A. Javed, T. Uthayakumar, M. O. D. Alotaibi, S. M. Al-Marzoug, H. Bahlouli and U. Al Khawaja, Communications in Nonlinear Science and Numerical Simulation, 103, 105968 (2021).
  * (7) T. Haug, R. Dumke, L. C. Kwek and L. Amico, Quantum Sci. Technol. 4, No. 4, 045001 (2019).
  * (8) W. C. Mak, B. A. Malomed, and P. L. Chu, Phys. Rev. E 67, 026608 (2003).
  * (9) Y. V. Kartashov, V. A. Vysloukh and L. Torner, Optics Express, 12, No 13, 2831 (2004).
  * (10) H. C. Gurgov and O. Cohen, Optics Express, 17, No. 9, 7052 (2009).
  * (11) J. J. Mazo and A. V. Ustinov, The sine-Gordon model and its applications, Springer, 155 (2014).
  * (12) N. Hatakenaka, H. Takayanagi, Y. Kasai and S. Tanda, Physica B: Condensed Matter, 284, 563 (2000).
  * (13) V. A. Gani and A. E. Kudryavtsev, Phys. Rev. E 60, No. 3, 3305 (1999).
  * (14) M. Croitoru, J. Phys. A: Math. Gen. 22, No. 7, 845 (1989).
  * (15) C. A. Popov, Wave Motion. 42, No. 4, 309 (2005).
  * (16) S. K. Varbeva, R. S. Kamburova, and M. T. Primatarowa, AIP Conference Proceedings, 2075, No. 1, 020009 (2019).
  * (17) B. A. Malomed, J. Opt. Soc. Am. B 31, No. 10, 2460 (2014).
  * (18) K. Maki and P. Kumar, Phys. Rev. B 14, 118 (1976).
  * (19) Y. Shiefman and P. Kumar, Phys. Scr. 20, 435 (1979).
  * (20) A. R. Bishop and T. Schneider,  Solitons and condensed matter physics. Springer series in solid state science, 8. Berlin, Heidelberg, New York: Springer; 246 (1978).
  * (21) Li. Y, Lu. X and C. Hou, Journal of Modern Optics 65.17, 2006 (2018).
  * (22) M. Saravanan, Chinese Journal of Physics 59, 426 (2019).
  * (23) T. Uchiyama, Phys. Rev. D 14, 3520 (1976).
  * (24) S. Duckworth, R. K. Bullough, P. J. Caudrey and J. D. Gibbon, Phys. Lett. 57 A, 19 (1976).
  * (25) Y. Shen, and Y. O. El-Dib, Journal of Low Frequency Noise, Vibration and Active Control 40.2, 683 (2021).
  * (26) H. Leblond, H. Triki, and D. Mihalache, Phys. Rev. A 86, 063825 (2012).
  * (27) V.Yu. Novokshenov, Asymptotic Analysis 12, No. 2, 77 (1996).
  * (28) Yue Shen and Yusry O El-Dib, Journal of Low Frequency Noise, , Vibration and Active Control, 40, No. 2, 683 (2021).
  * (29) W. A. Clarke, R. Marangell, Journal of Nonlinear Science 30, No. 6, 3421 (2020).
  * (30) D. Bazeia, L. Losano, J.M.C. Malbouisson, J.R.L. Santos, Eur. Phys. J. C 71, 1 (2011).
  * (31) M. Peyravi, N. Riazi and A. Montakhab, Eur. Phys. J. B 76, 547 (2010).
  * (32) M. Peyravi, A. Montakhab, N. Riazi and A. Gharaati, Eur. Phys. J. B 72, 269 (2009).
  * (33) N. Riazi and A. R. Gharaati, Int. J. Theor. Phys. 37, No. 3, 1081 (1998).
  * (34) J. Yang, Nonlinear Waves in Integrable and Nonintegrable Systems, SIAM, Philadelphia (2010).
  * (35) Y. Sun, Applied Mathematics and Computation, 258, 100 (2015).
  * (36) M. Wang and X. Li, Chaos, Solitons and Fractals, 27, No. 2, 477 (2006).
  * (37) H. C. Hu, S. Y. Lou and K. W. Chow, Chaos, Solitons and Fractals, 31, Issue 5, 1213 (2007).
  * (38) M. Pawellek, On a generalization of Jacobi’s elliptic functions and the double sine-Gordon kink chain, J. Math. Phys. 52, 113701 (2011).
  * (39) H. Liu and J. Li, Journal of Computational and Applied Mathematics 257, 144 (2014).
  * (40) Y. Sun, Appl. Math. and Comput 258, 100 (2015).
  * (41) S. P. Joseph, Int. J. Appl. Comput. Math 8, No. 1, 42 (2022).
  * (42) M. Guidry, Gauge Field Theories, an Introduction with Applications, Wiley, NewYork (1991).
  * (43) L. H. Ryder, Quantum Field Theory, Cambridge University Press (1985).
  * (44) N. Riazi, Int. J. Theor. Phys. GTNO, 8, No. 2, 115 (2001).
  * (45) D. Bazeia, L. Losano, J. M. C. Malbouisson and J. R. L. Santos, Eur. Phys. J. C 71, 1767 (2011).
  * (46) S. P. Joseph, Advances in Mathematics: Scientific Journal 9, No. 4, 2293 (2020).
  * (47) Y. Jian-Song and L. Sen-Yue, Chinese Physics Letters, 21, No. 4, 608 (2004).
  * (48) V .A. Gani, A. Moradi Marjaneh and D. Saadatmand, Eur. Phys. J. C 79, 620 (2019).
  * (49) A. Moradi Marjaneh, A. Askari, D. Saadatmand and S.V. Dmitriev, Eur. Phys. J. B 91, 1 (2018) .
  * (50) V. A. Gani, A. Moradi Marjaneh, A. Askari, E. Belendryasova and D. Saadatmand, Eur. Phys. J. C 78, 345 (2018).
  * (51) G. Delfino and G. Mussardo, Nuclear Physics B 516, Issue 3, 675 (1998).
  * (52) D. De Santis, C. Guarcello, B. Spagnolo, A. Carollo and D. Valenti, Chaos, Solitons and Fractals 170, 113382 (2023).
  * (53) D. De Santis, C. Guarcello, B. Spagnolo, A. Carollo and D. Valenti, Chaos, Solitons and Fractals 168, 113115 (2023).
  * (54) D. De Santis, C. Guarcello, B. Spagnolo, A. Carollo and D. Valenti, Chaos, Solitons and Fractals 158, 112039 (2022).
  * (55) C. Guarcello, D. Valenti, B. Spagnolo, V. Pierro and G. Filatrella, Phys. Rev. Appl. 11, No. 4, 044078 (2019).
  * (56) G. Parisi, Nobel Lecture: Multiple equilibria. Rev. Mod. Phys. 95, 030501 (2023).
  * (57) B. Lisowski, D. Valenti, B. Spagnolo, M. Bier and E. Gudowska-Nowak, Phys. Rev. E 91, No. 4, 042713 (2015).
  * (58) D. Valenti, G. Denaro, A. La Cognata, B. Spagnolo, A. Bonanno, G. Basilone, S. Mazzola, S. Zgozi and S. Aronica, Acta Phys. Pol. B 43, No. 5, 1227 (2012).
  * (59) R. Stassi, S. Savasta, L. Garziano, B. Spagnolo and F. Nori, New J. Phys. 18, No. 12, 123005 (2016).
  * (60) D. Valenti, G. Fazio, and B. Spagnolo, Phys. Rev. E 97, No. 6, 062307 (2018).
  * (61) C. Guarcello, D. Valenti, G. Augello and B. Spagnolo, Acta Phys. Pol. B 44, No. 5, 997 (2013).
  * (62) Y. V. Ushakov, A. A. Dubkov and B. Spagnolo, Phys. Rev. Lett. 107, No. 10, 108103 (2011).
  * (63) V. A. Gani, A. Moradi Marjaneh and K. Javidan, Eur. Phys. J C 81, 1124 (2021).
  * (64) K. Javidan, Journal of Physics A: Mathematical and General 39, No. 33, 10565 (2006).