# The nonlinear Benjamin-Feir instability – Hamiltonian dynamics, primitive
breathers, and steady solutions

David Andrade & Raphael Stuhlmeier

(Centre for Mathematical Sciences, University of Plymouth, Drake Circus, PL4
8AA, Plymouth, UK)

###### Abstract

We develop a general framework to describe the cubically nonlinear interaction
of a unidirectional degenerate quartet of deep-water gravity waves. Starting
from the discretised Zakharov equation, and thus without restriction on
spectral bandwidth, we derive a planar Hamiltonian system in terms of the
dynamic phase and a modal amplitude. This is characterised by two free
parameters: the wave action and the mode separation between the carrier and
the side-bands. The mode separation serves as a bifurcation parameter, which
allows us to fully classify the dynamics. Centres of our system correspond to
non-trivial, steady-state nearly-resonant degenerate quartets. The existence
of saddle-points is connected to the instability of uniform and bichromatic
wave trains, generalising the classical picture of the Benjamin-Feir
instability. Moreover, heteroclinic orbits are found to correspond to
primitive, three-mode breather solutions, including an analogue of the famed
Akhmediev breather solution of the nonlinear Schrödinger equation.

##  1 Introduction

The linear theory of water waves, developed and refined over the course of the
19${}^{\text{th}}$ and early 20${}^{\text{th}}$ centuries has endured by
virtue of its simplicity, elegance, and utility. When the governing equations
are linearised, waves may be superposed to create the myriad complex patterns
observed on the surface of the sea. The absence of interaction among such
waves leads, however, to an inability to explain certain properties of ocean
wave fields, particularly if these waves are steep or the time-scales
considered are long.

This interaction may occur between a wave and itself, as discovered by G. G.
Stokes in 1847, or as a mutual interaction between different waves. A natural
way to think about such weakly-nonlinear interactions is as occurring between
Fourier modes, each of which represents a periodic wave train (or plane wave)
with a distinct frequency, amplitude, and phase. When considering gravity
waves in deep water, the fundamental interaction is between four such Fourier
modes, which may or may not be distinct. This fact is not obvious a priori, as
it depends on the possibility of resonance between modes, and therefore on the
dispersion relation of the problem. The deep-water dispersion relation
$\omega^{2}=g|\mathbf{k}|$ between radian frequency $\omega$ and wavenumber
$\mathbf{k}$ provides a geometric constraint which means resonance is possible
for cubic (or higher) nonlinearity. This pioneering discovery by Phillips
(1960) soon became the object of intense study, with work to elucidate its
consequences for everything from wave spectra (Hasselmann, 1962) to the
evolution of uniform wave trains (so-called Stokes waves) (Benjamin and Feir,
1967; Zakharov, 1968).

One of the most remarkable early results was the finding that initially
uniform wave-trains are unstable, and tend to disintegrate into wave groups –
known as modulational or Benjamin-Feir instability after Benjamin and Feir
(1967). Indeed, uniform wave-trains are often difficult to produce
experimentally, a fact attributed to this instability. The perturbations
required to start the process of disintegration are small-amplitude side-
bands, or Fourier modes to either side of the plane wave. Interaction between
the modes transfers energy to the side-bands, which grow at the expense of the
plane wave (or carrier) and distort the pattern on the free-surface. Shortly
after the discovery of the instability Zakharov (1968) was able to reduce the
problem to Hamiltonian form, and to derive a simplified equation governing the
evolution of such wave groups for unidirectional waves in deep water, provided
the wave modes are not widely separated in Fourier space. This nonlinear
Schrödinger equation (NLS) was then used to recover Benjamin and Feir’s
criterion, specifying which ranges of side-band wavenumbers and carrier
steepness give rise to instability.

The nonlinear Schrödinger equation has also received considerable recent
attention in connection with the deterministic and stochastic modelling of
extreme waves, large wave crests which appear without warning on the surface
of the sea (a recent review of progress is given by Onorato and Suret (2016)).
The temporally or spatially localised exact solutions of the NLS, called
“breathers,” are of particular importance to our understanding of extreme
waves. Their recent experimental observation (Chabchoub et al., 2014) make
them an attractive model to generate extreme waves in the wave flume, and shed
light on the behaviour of such waves in the ocean. Indeed, each breather is a
particular manifestation of the underlying Benjamin-Feir instability for
special initial conditions.

The success of the NLS as a simplified model equation naturally engendered
interest in its agreement with more general water wave equations. The NLS
breather solutions, for example, were found to be robust when compared to
numerical solutions of the Euler equations by Slunyaev and Shrira (2013), who
also gave an overview of previous work in this vein. The Benjamin-Feir
stability threshold derived from the NLS has also been re-examined by Crawford
et al. (1981), who sought to remove the restriction to small mode separation.
However, the route to obtaining an instability criterion in their work
remained one of linearisation and the derivation of an eigenvalue criterion.
The smallness assumptions intrinsic to linear stability analysis mean that it
is limited to describing only the initial stages of evolution around
particular exact solutions. In this setting, the eigenvalues of the linearised
system govern the growth rate of initially small disturbances, whose
subsequent evolution can be captured by numerical simulation. However, a
numerical treatment has the potential to obscure fundamental, underlying
mechanisms, particularly in cases with many (ostensibly independent)
parameters. For example, computing the evolution of a carrier and two side-
bands requires the specification of eight parameters: two wavenumbers
$\mathbf{k}_{a},\mathbf{k}_{b},$ (the third is given from the resonance
condition $2\mathbf{k}_{a}=\mathbf{k}_{b}+\mathbf{k}_{c}$), three wave slopes
$\epsilon_{a},\epsilon_{b},\epsilon_{c},$ and three phases
$\theta_{a},\theta_{b},$ and $\theta_{c}.$

We aim to present a unified theory of the cubically nonlinear interaction of
three waves in deep water, which generalises the Benjamin-Feir instability.
Our starting point, like that of Crawford et al. (1981), is the reduced
Zakharov equation (ZE), which is free from any restrictions on mode
separation. Indeed, the nonlinear Schrödinger equation itself, as well as
later generalisations due to Davey and Stewartson (1974) or Dysthe (1979) can
all be derived as narrow-band limits of the ZE (see Stiassnie (1984); Gramstad
and Trulsen (2011)). After suitable transformations, the interaction of three
wave modes satisfying $2\mathbf{k}_{a}=\mathbf{k}_{b}+\mathbf{k}_{c}$ can be
recast as a planar Hamiltonian dynamical system. Phase plane analysis of this
system allows for a complete description of the dynamics for all times and
arbitrary initial conditions: these initial conditions are reduced to
specifying a single “dynamic phase” and a parameter specifying the
distribution of the wave action among the three modes. The mode-separation
distance plays the role of a bifurcation parameter in the problem.

With this simplification we are able to fully classify the dynamics of such
“degenerate quartets” of unidirectional waves. Fixed points in the phase-plane
correspond to steady-state near-resonant cases, and heteroclinic orbits are
seen to correspond to a variety of primitive breather-type solutions. We
obtain general instability results for uniform and bichromatic wave trains as
criteria for the existence of fixed-points, and show that the well-known
results of linear stability analysis can be recovered. In what follows we
first present the reformulation of the discrete Zakharov equation for a
degenerate quartet of waves in section 2. In section 3 we classify the fixed
points and phase-portraits, using mode separation as a bifurcation parameter.
The nonlinear stability of wave trains is considered in section 4, and
special, heteroclinic solutions are discussed in section 5. Finally, a
discussion of our results is presented in section 6. Appendix A contains some
simplified expressions for integral kernels used in computations.

##  2 The Zakharov equation for a degenerate quartet of waves

Investigating third-order wave-wave interaction on deep water without a
bandwidth restriction means that our starting point will be the reduced
Zakharov equation

$i\frac{dB_{n}}{dt}=\sum_{p,q,r=1}^{N}T_{npqr}\delta_{np}^{qr}e^{i\Delta_{npqr}t}B_{p}^{*}B_{q}B_{r},\,n=1,\,2,\ldots,N$
(1)

derived by Zakharov (1968), and in Hamiltonian form by Krasitskii (1994), and
here written in a convenient discrete formulation. While the equation itself
has the generic form of a cubic interaction equation (see Zakharov et al.
(1992)), the physics of the water-wave problem are contained in the kernels
$T_{npqr}=T(\mathbf{k}_{n},\mathbf{k}_{p},\mathbf{k}_{q},\mathbf{k}_{r}).$ We
use $\delta_{np}^{qr}$ to denote a Kronecker delta function

$\delta_{np}^{qr}=\begin{cases}1\text{ for
}\mathbf{k}_{n}+\mathbf{k}_{p}=\mathbf{k}_{q}+\mathbf{k}_{r}\\\ 0\text{
else.}\end{cases}$

and $\Delta_{npqr}=\omega_{n}+\omega_{p}-\omega_{q}-\omega_{r}$ to denote the
frequency detuning (a measure of departure from exact resonance). Throughout

$\omega_{n}^{2}=g|\mathbf{k}_{n}|$ (2)

is the linear dispersion relation for gravity waves in deep water.

The simplest non-trivial interaction is between three waves in deep water
$\mathbf{k}_{a},\mathbf{k}_{b}$ and $\mathbf{k}_{c},$ where one wave is
counted twice to satisfy the resonance condition
$2\mathbf{k}_{a}=\mathbf{k}_{b}+\mathbf{k}_{c}.$ Indeed this case – termed the
“degenerate quartet” – corresponds to the Benjamin-Feir instability (see Yuen
and Lake (1982) or (Mei et al., 2018, Ch. 14.9)), where mode $\mathbf{k}_{a}$
is interpreted as a “carrier wave” and modes $\mathbf{k}_{b}$ and
$\mathbf{k}_{c}$ as “side-bands.”

Rather than handle three complex equations, it is a convenient first step
towards simplification to write the discrete equation in terms of amplitude
and phase variables (see, e.g. Craik (1986) and references therein). The
equations then become

$\displaystyle\frac{d|B_{a}|}{dt}$
$\displaystyle=2T_{aabc}|B_{a}||B_{b}||B_{c}|\sin(\theta_{aabc}),$ (3a)
$\displaystyle\frac{d|B_{b}|}{dt}$
$\displaystyle=-T_{aabc}|B_{a}|^{2}|B_{c}|\sin(\theta_{aabc}),$ (3b)
$\displaystyle\frac{d|B_{c}|}{dt}$
$\displaystyle=-T_{aabc}|B_{a}|^{2}|B_{b}|\sin(\theta_{aabc}),$ (3c)
$\displaystyle\frac{d\theta_{a}}{dt}$
$\displaystyle=\Gamma_{a}+2\frac{T_{aabc}}{|B_{a}|^{2}}|B_{a}|^{2}|B_{b}||B_{c}|\cos(\theta_{aabc}),$
(3d) $\displaystyle\frac{d\theta_{b}}{dt}$
$\displaystyle=\Gamma_{b}+\frac{T_{aabc}}{|B_{b}|^{2}}|B_{a}|^{2}|B_{b}||B_{c}|\cos(\theta_{aabc}),$
(3e) $\displaystyle\frac{d\theta_{c}}{dt}$
$\displaystyle=\Gamma_{c}+\frac{T_{aabc}}{|B_{c}|^{2}}|B_{a}|^{2}|B_{b}||B_{c}|\cos(\theta_{aabc}),$
(3f)

where $\theta_{aabc}=\Delta_{aabc}t+2\theta_{a}-\theta_{b}-\theta_{c}$ is
called the dynamic phase. Note the additional factor of two appearing in the
equations for $|B_{a}|$ and $\theta_{a}$ because the resonance condition
$2\mathbf{k}_{a}=\mathbf{k}_{b}+\mathbf{k}_{c}$ can be satisfied for either
the tuple $(\mathbf{k}_{a},\mathbf{k}_{a},\mathbf{k}_{b},\mathbf{k}_{c})$ or
$(\mathbf{k}_{a},\mathbf{k}_{a},\mathbf{k}_{c},\mathbf{k}_{b})$.

In equations (3d) - (3f)

$\Gamma_{i}=|B_{i}|^{2}T_{i}+2\sum_{j\neq i}|B_{j}|^{2}T_{ij},$ (4)

where we use an abbreviated notation for the symmetric kernels:
$T_{i}=T_{iiii},\,T_{ij}=T_{ijij}=T_{jiji}$.

The next significant simplification relies on the observation that, although
there are ostensibly three distinct phases in the equations governing the
degenerate quartet (3a)–(3f), they occur only in the single combination
$2\theta_{a}-\theta_{b}-\theta_{c},$ which makes it possible to drop
subscripts and write $\theta$ and $\Delta$ in place of $\theta_{aabc}$ and
$\Delta_{aabc}$ without risk of confusion. We can further exploit this fact to
combine (3d) - (3f) into a single equation for the dynamic phase

$\displaystyle\frac{d\theta}{dt}$
$\displaystyle=\Delta+2\Gamma_{a}-\Gamma_{b}-\Gamma_{c}+T_{aabc}\left(4|B_{b}||B_{c}|-\frac{|B_{c}||B_{a}|^{2}}{|B_{b}|}-\frac{|B_{b}||B_{a}|^{2}}{|B_{c}|}\right)\cos(\theta).$
(5)

Equations (3a) - (3c) and (5) now form an autonomous system of ordinary
differential equations. It is known that this system of equations admits
periodic solutions which are given in terms of Jacobi elliptic functions, see
Shemer and Stiassnie (1985). We shall see that generic solutions are indeed
periodic, though we shall focus our attention principally on special, non-
periodic solutions.

The system of equations (3a)-(3c) and (5) admits the following conserved
quantities:

$\displaystyle|B_{a}|^{2}+|B_{b}|^{2}+|B_{c}|^{2}=A,$ (6)

the total wave action, and

$\displaystyle\frac{|B_{b}|^{2}}{A}-\frac{|B_{c}|^{2}}{A}=\alpha,$ (7)

the wave-action difference in the two side-bands. These conserved quantities
allow us to further reduce the number of parameters, by introducing a
normalized wave-action variable $\eta=|B_{a}|^{2}/A$, which by (6) must take
values between $0$ and $1$. The amplitudes may then be rewritten in terms of
$\eta$ and $\alpha$ as $|B_{a}|^{2}=A\eta$, $|B_{b}|^{2}=A(1-\eta+\alpha)/2$
and $|B_{c}|^{2}=A(1-\eta-\alpha)/2$.

In terms of the two variables $\eta$ and $\theta$ we find that the system can
be described by the Hamiltonian

$\displaystyle
H(\theta,\eta)=-AT_{aabc}2\eta\sqrt{(1-\eta)^{2}-\alpha^{2}}\cos(\theta)-(\Delta+A\Omega_{0})\eta-\frac{A\Omega_{1}}{2}\eta^{2},$
(8)

with

$\displaystyle\frac{\partial H}{\partial\theta}$
$\displaystyle=\frac{d\eta}{dt}=2AT_{aabc}\eta\sqrt{(1-\eta)^{2}-\alpha^{2}}\sin(\theta),$
(9) $\displaystyle-\frac{\partial H}{\partial\eta}$
$\displaystyle=\frac{d\theta}{dt}=\Delta_{aabc}+A\Omega_{0}+A\Omega_{1}\eta+2AT_{aabc}\frac{1+2\eta^{2}-3\eta-\alpha^{2}}{\sqrt{(1-\eta)^{2}-\alpha^{2}}}\cos(\theta).$
(10)

The two new coefficients are

$\displaystyle\Omega_{0}$
$\displaystyle=2\left[T_{ab}(1+\alpha)+T_{ac}(1-\alpha)\right]-\frac{1}{2}\left[T_{b}(1+\alpha)+T_{c}(1-\alpha)\right]-2T_{bc}$
(11) $\displaystyle\Omega_{1}$
$\displaystyle=2\left(T_{a}+T_{bc}\right)-4\left(T_{ab}+T_{ac}\right)+\frac{1}{2}\left(T_{b}+T_{c}\right).$
(12)

Such a Hamiltonian approach to the three-wave discretisation of the nonlinear
Schrödinger equation was first explored by Cappellini and Trillo (1991) in the
context of optics, and our Hamiltonian (8) can be shown to reduce in the
narrow-band limit to one analogous to that presented therein.

A final simplification of the system can be achieved by imposing an
equidistribution of wave action among the side-bands, i.e. $\alpha=0.$ Under
this assumption equations (9) and (10), as well as the Hamiltonian (8), become

$\displaystyle\frac{d\eta}{dt}$
$\displaystyle=2AT_{aabc}\eta(1-\eta)\sin(\theta),$ (13)
$\displaystyle\frac{d\theta}{dt}$
$\displaystyle=\Delta+A\Omega_{0}+A\Omega_{1}\eta+2AT_{aabc}(1-2\eta)\cos(\theta),$
(14) $\displaystyle H(\theta,\eta)$
$\displaystyle=-2AT_{aabc}\eta(1-\eta)\cos(\theta)-(\Delta+A\Omega_{0})\eta-\frac{A\Omega_{1}}{2}\eta^{2}.$
(15)

In subsequent computations we shall normalise all wavenumbers by the carrier
wavenumber $k_{a},$ and all frequencies by the carrier frequency $\omega_{a}$,
which is equivalent to setting the gravitational acceleration $g=1$. For the
unidirectional cases considered we shall also explicitly specify our
degenerate quartets in terms of mode separation $p$ as follows:
$\mathbf{k}_{a}=[1,0]$, $\mathbf{k}_{b}=[1-p,0]$ and $\mathbf{k}_{c}=[1+p,0]$,
where $0<p<1.$ We henceforth drop the bold-script and write $k_{i}$ for the
wavenumbers, to emphasise that these are scalars. Indeed, for unidirectional
waves in deep water it is possible to significantly simplify the kernels, as
shown by Dyachenko et al. (2017); Kachulin et al. (2019) and detailed in
Appendix A.

The wave action $A$ of the three-mode system remains a free parameter, and the
conservation law (6) makes it clear that $A$ may be “distributed” among the
three modes in various configurations. The relationship between the complex
amplitudes $B_{i}$ of the Zakharov formulation and the amplitude of the free
surface displacement (e.g. (14.5.5) of (Mei et al., 2018, Sec. 14.5)) makes it
possible to interpret $A$ in terms of the carrier steepness, which is
sometimes convenient. If only the carrier wave is present

$A=|B_{a}|^{2}=\frac{2ga^{2}\pi^{2}}{\omega_{a}}=\frac{2\pi^{2}g^{1/2}\epsilon^{2}}{k_{a}^{5/2}},$
(16)

where $a$ is the surface-wave amplitude and $\epsilon=ak_{a}$ is the wave
slope.

Aside from the wave action $A,$ the system (13) and (14) contains a further
free parameter: the separation between the Fourier modes $p.$ Together these
can be used to determine the values of all coefficients in the equation.

##  3 Hamiltonian dynamics in the phase-plane

The Hamiltonian dynamical system (13) – (14) describes fully the nonlinear
dynamics of three interacting deep-water waves. We can gain both qualitative
and quantitative understanding of this system via phase–plane analysis, where
the phase space is the surface of the truncated cylinder
$\\{(\theta,\eta)\mid-\pi\leq\theta\leq\pi,\,0\leq\eta\leq 1\\}.$

###  3.1 Fixed points – steady state degenerate quartets

The first step in unraveling the phase dynamics is to consider fixed points of
our system. The trajectories, corresponding to solutions of the degenerate
quartet with particular initial conditions, occupy the level curves of the
Hamiltonian (15). We find that fixed points of (13) – (14) occur in four
classes: $\eta=0$ or 1 and $\theta=0$ or $\pm\pi.$ The upper and lower
boundaries $\eta=1$ and $\eta=0$ correspond to a single wave train and a
bichromatic wave train, respectively (recall $\eta=|B_{a}|^{2}/A$). These two
well-known solutions of the Zakharov equation exhibit no energy exchange (see
Leblanc (2009) or Mei et al. (2018, Sec. 14.5 & 14.6)), since $d\eta/dt=0$.
Indeed, the effect of nonlinear interaction in such cases is solely to induce
a frequency correction (called Stokes’ correction after Stokes (1847), and
found for bichromatic waves in Longuet-Higgins & Phillips (1962); for a
discussion in the context of the Zakharov equation see Stuhlmeier & Stiassnie
(2019)).

To avoid overly bulky expressions in what follows we introduce the following
notation

$\displaystyle\Delta^{\prime}=\frac{\Delta}{AT_{aabc}},\quad\Omega_{0}^{\prime}=\frac{\Omega_{0}}{T_{aabc}},\quad\text{
and }\quad\Omega_{1}^{\prime}=\frac{\Omega_{1}}{T_{aabc}}.$ (17)

In these variables we obtain fixed points at $\eta=1$ of the form
$(\theta,\eta)=(\theta_{\pm 1},1)$ whenever $\theta_{1}$ satisfies

$\displaystyle\cos(\theta_{1})=\frac{\Delta^{\prime}+\Omega_{0}^{\prime}+\Omega_{1}^{\prime}}{2}.$
(18)

Such fixed points exist when the right-hand side is between -1 and 1, and
$\theta_{-1}=-\theta_{1}.$ These fixed points lie on the contour
$H=-\Delta-A\Omega_{0}-\frac{A\Omega_{1}}{2}.$

Similarly, we obtain fixed points at $\eta=0$ of the form
$(\theta,\eta)=(\theta_{\pm 0},0)$ whenever $\theta_{0}$ satisfies

$\displaystyle\cos(\theta_{0})=-\frac{\Delta^{\prime}+\Omega_{0}^{\prime}}{2},$
(19)

again provided the right hand side is between -1 and 1 and
$\theta_{-0}=-\theta_{0}$. These fixed points lie on the contour $H=0.$

Another class of fixed points is obtained for $\theta=0$ or $\pm\pi.$ In the
former case $(\theta,\eta)=(0,\eta_{0})$ is a fixed point whenever

$\displaystyle\eta_{0}=\frac{2+\Delta^{\prime}+\Omega_{0}^{\prime}}{4-\Omega_{1}^{\prime}},$
(20)

and $0\leq\eta_{0}\leq 1$. Note that $\eta_{0}=1$ if and only if $\theta_{\pm
1}=0,$ in which case this fixed point coincides with (18) above. For
$\theta=\pm\pi$ we find fixed points $(\pm\pi,\eta_{\pi})$ provided

$\displaystyle\eta_{\pi}=\frac{2-(\Delta^{\prime}+\Omega_{0}^{\prime})}{4+\Omega_{1}^{\prime}},$
(21)

and $0\leq\eta_{\pi}\leq 1$. Similarly, $\eta_{\pi}=0$ if and only if
$\theta_{\pm 0}=\pm\pi,$ such that this fixed point coincides with (19).

###  3.2 Phase portraits and bifurcation

The dynamical system (13) – (14) contains two parameters – wave action $A$ and
mode-separation $p$ – which together govern the existence of the fixed-points
given in equations (18)–(21), and the attendant dynamics. The natural choice
is to fix $A$ – akin to specifying the total energy of waves – and to use mode
separation as a bifurcation parameter.

In figure 1 we depict a series of characteristic phase portraits, where mode
separation $p$ increases from $p=0$ in Case (a) to $p=0.3$ in Case (h). In
this figure we have selected $A$ equivalent to a carrier wave with steepness
$\epsilon=0.1,$ so that $A=2\pi^{2}/100$ (see (16) and recall $g=k_{a}=1$).
The figure thus represents both physically realistic, as well as
characteristic behaviour, in a sense to be specified in detail below.

Figure 1: Phase portraits for $A=2\pi^{2}/100$ and various values of mode
separation $p$ between 0 and 0.3, plotted on the cylinder
$\\{(\eta,\theta)\mid\eta\in[0,1],\theta\in[-\pi,\pi]\\}.$ Coloured curves
represent contours of the Hamiltonian (15). Dashed lines show separatrices,
while dots denote fixed-points of the system.

Deferring the details for the moment we can give an overview of the dynamics
of our degenerate quartet as obtained by phase-plane analysis, and shown in
figure 1. In the limit of vanishing mode separation, Case (a), we find four
fixed points: a centre at $\theta=0,$ a pair of saddle points on $\eta=0,$ and
a semi-stable fixed point at $\theta=\pm\pi,\,\eta=1.$ The first bifurcation
occurs at $p=0$, and as $p$ increases the semi-stable point at
$(\theta,\eta)=(\pm\pi,1)$ splits into three fixed points: a centre at
$\theta=\pm\pi$ and two saddles at $\eta=1$ (Case (b)). Further increase in
the mode separation causes the centre at $\theta=\pm\pi$ to descend until the
saddle connections between the two fixed points at $\eta=0$ and $\eta=1$ merge
(Case (c)).

Further increasing $p$ leads the centre at $\theta=\pm\pi$ to descend towards
$\eta=0$ (Case (d)) where it subsequently vanishes together with the saddle
points along $\eta=0$ (Case (e)). As the modes are separated further the
remaining three fixed points draw closer together in phase space (Case (f))
before coalescing (Case (g)) and disappearing entirely (Case (h)) – this last
bifurcation leads to complete stabilisation of the system.

##  4 Nonlinear stability

While the phase-plane analysis presented above is remarkably simple it
contains a wealth of information about the fully nonlinear dynamics of
interacting degenerate quartets of deep water waves. We focus first on a
discussion of nonlinear instability results.

###  4.1 Nonlinear instability of a uniform wave train

Many classical studies of the stability of a uniform wave train begin by
establishing that such a monochromatic wave is a solution of the relevant
governing equation. This solution is then used as a starting point for
linearisation, and an instability criterion derived from the eigenvalues of
the linear system, as in Crawford et al. (1981). In our framework, uniform
wave trains are described by $\eta=1,$ for arbitrary values of $\theta.$

The initial small disturbance used classically to investigate the Benjamin-
Feir instability consists in imposing small side-bands, i.e. a shift to a
contour $\eta<1.$ In cases (a) through (f) of figure 1 we readily observe that
such a shift leads to growth of the side-bands, as $\eta$ decreases along the
contours, followed generically by periodic energy exchange. Maximal energy
exchange occurs when $\eta$ is changing fastest, and coincides with the
smallest changes in the dynamic phase $\theta.$ Such phase coherence has also
been observed in simulations, e.g. Houtani et al. (2022); Liu et al. (2021).

In fact we can show that the classical linear stability criterion is
equivalent to the existence of fixed points at $\eta=1$ in the nonlinear
system. The condition (18) for such fixed points to exist is

$\displaystyle-1\leq\frac{\Delta^{\prime}+\Omega_{0}^{\prime}+\Omega_{1}^{\prime}}{2}\leq
1.$ (22)

Squaring these inequalities, substituting $\Omega_{0}$ and $\Omega_{1}$ from
(11) and (12), and using (16) yields, after simplification, the discriminant
criterion

$\displaystyle\left(\frac{\Delta}{2}+(T_{a}-T_{ab}-T_{ac})|B_{a}|^{2}\right)^{2}-|B_{a}|^{4}T_{aabc}^{2}\leq
0,$ (23)

see (Mei et al., 2018, Eq. (14.9.16)). It is rather surprising that this
eigenvalue condition, which arises from an approach based on small amplitude
perturbations (and which is oblivious to the existence of fixed-points of the
original nonlinear problem) can be recovered exactly with our approach.

####  4.1.1 Restabilisation

The existence of fixed points at $\eta=1,$ and thus the instability of uniform
wave trains to some perturbation, is the generic situation for degenerate
quartets. We can appreciate this by considering equation (18), from which we
establish that fixed points exist at $\eta=1$ for some mode separation,
provided $A\in\left[0,{2\pi(2-\sqrt{2})\sqrt{g}}{k_{a}^{-5/2}}\right].$
Employing (16) this can be expressed in terms of carrier wave slope
$\epsilon,$ and implies the existence of fixed points (and thus instability)
for

$0<\epsilon<\sqrt{2-\sqrt{2}},$

far beyond the wave breaking threshold.

Figure 2: Existence of $\eta=1$ fixed points in $(\epsilon,p)$-parameter space
according to (18). The shaded region denotes the nonlinear instability domain
for the degenerate quartet

In the limit of small mode separation, we see that the instability domain
remains bounded (see figure 2), with fixed points for
$A\in[0,\frac{\sqrt{gk_{a}}\pi^{2}}{3k_{a}}].$ This yields the threshold

$0<\epsilon<\frac{1}{\sqrt{6}},$

which has previously been obtained numerically (see e.g. (Yuen and Lake, 1982,
Sec. VI.B.1)). In fact, this restabilisation for nearby side-bands (or long-
wavelength disturbances) is a characteristic of the broad-banded Zakharov
equation. Indeed, imposing a limit of narrow spectral bandwidth (in which the
Zakharov equation reduces to the NLS) implies

$\Omega^{\prime}_{0}\rightarrow 1,\,\Omega^{\prime}_{1}\rightarrow-3,\text{
and }\Delta^{\prime}\rightarrow\frac{\omega_{a}p^{2}\pi^{2}}{Ak_{a}^{5}}.$

This can be used to reformulate the fixed-point critrion and recover the well-
known instability criterion for NLS (Mei et al., 2018, Eq. (14.9.21))

$\frac{p}{k_{a}}\in\left[0,2\sqrt{2}\epsilon\right].$

As mode-separation increases, larger values of carrier steepness are required
to generate instability, up to the limit $p=1$ when the instability domain
contracts to the point $\epsilon=\sqrt{2-\sqrt{2}}.$ Conversely, for a given
carrier steepness and increasing mode separation the system will eventually
exit the instability domain – this bifurcation is precisely Case (g) in figure
1, where fixed points at $\eta=1$ coalesce and disappear.

###  4.2 Nonlinear stability of a bichromatic wave train

Instabilities of bichromatic wave trains, while less well-known than those for
uniform wave trains, have also been studied, for example by Leblanc (2009),
Badulin et al. (1995), or Ioualalen and Kharif (1994). While our setting is
somewhat different – in particular, we have only three unidirectionl modes, in
contrast to the four modes used by (Leblanc, 2009, Sec. 4) to discuss class Ib
instabilities – we can nevertheless obtain the nonlinear evolution of such
instabilities by our phase-plane analysis, recalling that a bichromatic wave
train occupies the line $\eta=0.$

In the present case, where the Fourier amplitudes of the side-bands $k_{b}$
and $k_{c}$ are taken to be identical, we naturally have a bichromatic wave
train with waves of different steepness. Mode $k_{b}$ corresponds to a longer
wave than mode $k_{c},$ and their slopes are related as

$\epsilon_{b}^{2}=\left(\frac{k_{b}}{k_{c}}\right)^{5/2}\epsilon_{c}^{2},$

see (Mei et al., 2018, Sec. 14.6). Without loss of generality we depict the
instability domain in terms of $\epsilon_{b}$ in figure 3.

Figure 3: Existence of $\eta=0$ fixed points in $(\epsilon_{b},p)$-parameter
space, according to (19). The shaded region denotes the instability domain for
a bichromatic wave train with equipartitioned mode amplitudes.

The existence of fixed points – and consequent instability – of the
bichromatic wave train is again seen to be generic. For any wave train, there
exists a value of mode separation $p$ such that the solution is unstable
within the framework of the degenerate quartet. This wave train stabilises for
sufficiently large mode separation, as seen in the bifurcation in Case (e) in
figure 1. Moreover, we find a critical threshold for mode separation $p\approx
0.42$ beyond which bichromatic wave trains are stable.

##  5 Heteroclinic solutions

The existence of special solutions to the nonlinear Schrödinger equation (NLS)
has attracted considerable attention in recent years, a discussion of which
may be found in Dysthe and Trulsen (1999) or, with an emphasis on
hydrodynamics and experimental verification, Chabchoub et al. (2014) and the
references therein. We will see that several remarkable solutions –
corresponding to heteroclinic orbits in the phase plane – exist for the three
wave system. These include a primitive breather analogous to that found by
Akhmediev et al. (1987), which approaches a plane wave as
$t\rightarrow\pm\infty.$

The form of the free-surface and its envelope are of particular interest for
such solutions; these can be recovered from the solution of the Zakharov
equation (1) by defining the following complex amplitude function

$\displaystyle
A(x,t)=\frac{1}{\pi}\left(\sqrt{\frac{\omega_{a}}{2g}}B_{a}e^{i(k_{a}x-\omega_{a}t)}+\sqrt{\frac{\omega_{b}}{2g}}B_{b}e^{i(k_{b}x-\omega_{b}t)}+\sqrt{\frac{\omega_{c}}{2g}}B_{c}e^{i(k_{c}x-\omega_{c}t)}\right).$
(24)

The free surface elevation is then obtained as

$\displaystyle\zeta(x,t)=\Re\left[A(x,t)\right],$ (25)

and the envelope is $|A(x,t)|$.

###  5.1 Primitive Breather solutions

The existence of special heteroclinic solutions, with properties similar to
the well-known breather solutions of the nonlinear Schrödinger equation,
clearly depends on the fixed points of our problem.

Figure 4: Time evolution of the free surface (right panels) and their
envelopes (left panels) along different heteroclinic solutions. In all cases
the steepness of the wave $k_{a}$ is set to $\epsilon_{a}=0.2$. Top panels:
Primitive Akhmediev breather for $p=0.3012$. Middle panels: (1-0) breather for
the critical $p=0.1670$. Bottom panels: (0-0) breather for $p=0.1570$.
$\tilde{x}=x/\lambda_{a}$, $\tilde{t}=t\omega_{a}$ and
$\tilde{z}=z/\epsilon_{a}.$

####  5.1.1 The primitive (1-1) Akhmediev breather

For values of $p$ and $A$ satisfying (18) (see section 4.1 and figure 2), the
two saddle points $(\theta,\eta)=(\pm\theta_{1},1)$ are connected by a
heteroclinic orbit. Along this orbit $\eta<1,$ and a trajectory approaches
$(\pm\theta_{1},1)$ as $t\to\pm\infty$, respectively.

In fact, we can compute these solutions explicitly, by considering the
following set of equations

$\displaystyle 2\eta\cos(\theta)$
$\displaystyle=\Delta^{\prime}+\Omega_{0}^{\prime}+\frac{\Omega_{1}^{\prime}}{2}\left(1+\eta\right).$
(26a) $\displaystyle\frac{d\eta}{d\theta}$
$\displaystyle=\frac{\eta\sin(\theta)}{\cos(\theta)-\frac{\Omega_{1}}{4}}.$
(26b) $\displaystyle\frac{d^{2}\theta}{dt^{2}}$
$\displaystyle=-2\sin(\theta)\frac{d\theta}{dt}.$ (26c)

Equation (26a) follows from equating the Hamiltonian (15) to its value at
$\eta=1$. Equation (26b) is obtained from the implicit function theorem by
regarding $\eta$ as a function of $\theta$ instead of $t$. The last equation
(26c) is obtained by taking the time derivative of equation (14). Equation
(26a) is further used to simplify the form of (26b) and (26c).

Equation (26b) is separable, and can be integrated directly to give

$\displaystyle\eta=\frac{4\cos(\theta_{1})-\Omega_{1}}{4\cos(\theta)-\Omega_{1}},$
(27)

where the integration constant is chosen so that the limits as $\theta$ tends
to $\pm\theta_{1}$ are 1.

Equation (26c) can also be integrated to obtain

$\displaystyle\frac{d\theta}{dt}=2\cos(\theta)-2\cos(\theta_{1}),$ (28)

where the integration constant is chosen so that the limits as $\theta$ tends
to $\pm\theta_{1}$ are 0. Integrating again yields an expression for the
dynamic phase

$\displaystyle\tan(\theta/2)=\tan(\theta_{1}/2)\tanh(\sin(\theta_{1})t),$ (29)

where the integration constant was chosen so that $\theta(0)=0.$ Together,
(29) and (27) describe the heteroclinic orbits shown in Case (d)–(f) of figure
1. (Note that the heteroclinic orbits around $\theta=\pm\pi$ in Case (b) can
be treated similarly.)

The primitive (1-1) Akhmediev breather describes an orbit in phase space which
connects one plane-wave solution with another. During the evolution energy is
transferred to the side-bands ($\eta$ decreases) while the dynamic phase
undergoes little change. Subsequent to this focusing the energy-transfer
process reverses, resulting in a phase-shifted plane wave. This situation is
depicted in dimensionless coordinates in figure 4 (top panels). The left panel
shows the envelope, which is periodic in space and “breathes” once at time
${t}=0,$ akin to the Akhmediev breather solution of the NLS. The right panel
shows the free surface, which begins as a monochromatic plane wave, grows into
a strongly modulated wave-train at ${t}=0,$ and reverts to a plane wave with
an evident phase shift.

####  5.1.2 (1-0) breather and (0-1) breathers

A special type of heteroclinic solution appears when the Hamiltonian at
$\eta=1$ vanishes, i.e. when $p$ and $A$ are such that
$\Delta^{\prime}+\Omega_{0}^{\prime}+\Omega_{1}^{\prime}/2=0$. There are two
such heteroclinic solutions, one linking fixed points $(\theta_{0},0)$ and
$(\theta_{1},1)$ and the other linking $(-\theta_{1},1)$ with
$(-\theta_{0},0)$ (see Case (c) in figure 1).

In this case equation (13) simplifies to

$\displaystyle\frac{d\eta}{dt}=\pm\eta(1-\eta)\sqrt{4-\frac{\Omega_{1}^{2}}{4}},$
(30)

where the sign depends on which fixed point we are considering. Integrating
this equation yields

$\displaystyle\eta=\frac{e^{\pm\sqrt{4-\frac{\Omega_{1}^{2}}{4}}t+C}}{e^{\pm\sqrt{4-\frac{\Omega_{1}^{2}}{4}}t+C}+1},$
(31)

where $C$ is an integration constant depending on the initial conditions. For
instance $C=0$ for the initial condition $\eta(0)=\frac{1}{2}$. Assuming that
a positive sign is chosen, the solution tends to $1$ as $t\to\infty$ and it
tends $0$ as $t\to-\infty$. The limits are reversed if the negative sign is
chosen.

The special case described by these orbits corresponds to a complete energy
transfer from a uniform wave-train to a bichromatic wave-train (or vice
versa). The envelope and free surface are depicted in the middle left and
right panels of figure 4, respectively. These show how spatially-periodic
modulation appears “out of nowhere” with time ${t},$ or disappears as time is
traversed in the negative sense.

####  5.1.3 (0-0) breather

For a given (sub-critical) value of mode separation $p,$ and $A$ such that
(19) is satisfied (see 4.2 and figure 3), two distinct fixed points
$(\theta,\eta)=(\pm\theta_{0},0)$ are linked by a heteroclinic solution that
approaches $(\pm\theta_{0},0)$ as $t\to\mp\infty,$ and along which $\eta>0.$

This solution satisfies the following set of equations:

$\displaystyle 2(\eta-1)\cos(\theta)$
$\displaystyle=\Delta^{\prime}+\Omega_{0}^{\prime}+\frac{\Omega_{1}^{\prime}}{2}\left(1+\eta\right),$
(32a) $\displaystyle\frac{d\eta}{d\theta}$
$\displaystyle=\frac{(1-\eta)\sin(\theta)}{\frac{\Omega_{1}}{4}-\cos(\theta)},$
(32b) $\displaystyle\frac{d^{2}\theta}{dt^{2}}$
$\displaystyle=2\sin(\theta)\frac{d\theta}{dt}.$ (32c)

As above, equation (32a) has been used to simplify equations (32b) and (32c).

Integrating equation (32b) yields

$\displaystyle\eta=2\frac{2\cos(\theta)+\Delta+\Omega_{0}}{4\cos(\theta)-\Omega_{1}}=\frac{4\cos(\theta)-4\cos(\theta_{0})}{4\cos(\theta)-\Omega_{1}},$
(33)

where the integration constant is chosen so that the limits as $\theta$ tends
to $\theta_{0}$ are 0.

Integrating equation (32c) yields

$\displaystyle\frac{d\theta}{dt}=-2\cos(\theta)+2\cos(\theta_{0}),$ (34)

where the integration constant is chosen so the limits as $\theta$ tends to
$\pm\theta_{0}$ are 0. Further integration yields

$\displaystyle\tan\left(\theta/2\right)=-\tan(\theta_{0}/2)\tanh\left(\sin(\theta_{0})t\right).$
(35)

This solution is the natural counterpart to the primitive Akhmediev breather
presented in section 5.1.1 – however, rather than tending to a uniform wave
train with $t\rightarrow\pm\infty$ it tends to a bichromatic wave train, again
with an attendant phase shift. The envelope and free-surface for such a case
are shown in the bottom left and right panels of figure 4, respectively. At
${t}=0,$ the modulation is at a minimum, and the free surface elevation is
close to a uniform wave-train.

###  5.2 Limiting solutions

It is clear from equation (15) that the horizontal lines $\eta=1$ and $\eta=0$
are level lines of the Hamiltonian. Hence, any solution of equations (13) and
(14) with suitable initial conditions so that $\eta=1$ or $\eta=0$ will remain
on these lines.

As with the fixed points themselves, there is no energy exchange among the
wave-modes for these solutions. However, they are not steady (or stationary)
solutions in the usual sense, because $\theta$ still satisfies the following
differential equation:

$\displaystyle\frac{d\theta}{dt}=\begin{cases}\Delta+A\Omega_{0}+A\Omega_{1}-2AT_{aabc}\cos(\theta),\quad\text{for
$\eta=1$,}\\\ \Delta+A\Omega_{0}+2AT_{aabc}\cos(\theta),\quad\text{for
$\eta=0$.}\end{cases}$ (36)

A further differentiation of these equations with respect to $t$ yields
equations (26c) and (32c) respectively, which means that equations (29) and
(35) are solutions for each case. Both limiting configurations of the system
exhibit what is sometimes referred to as phase-locking or phase coherence;
whereby regardless of the initial phase the dynamic phase converges to the
fixed value $-\theta_{1}$ at $\eta=1$ and $\theta_{0}$ at $\eta=0$. It is
expected that if initial conditions are taken close to $\eta=1$, i.e. if some
small amount of energy is initially put in the side bands, then a similar
behaviour will be observed in the evolution of the combined phase, as has been
recently reported by Houtani et al. (2022).

From the mathematical point of view both limiting configurations $\eta=1$ and
$\eta=0$ of the system can be understood as the limit of the generic periodic
solutions, as all the energy goes to either the wave $k_{a}$, or is
equipartitioned among the side bands $k_{b}$ and $k_{c}$, respectively. From
the physical point of view the dynamic phase $\theta$ disappears at both
$\eta=1$ or $\eta=0$ configurations and one is left with a classical
monochromatic Stokes waves or two co-propagating Stokes waves, respectively.

##  6 Discussion and conclusions

The reformulation of the discrete Zakharov equation for three modes as a two-
dimensional Hamiltonian dynamical system provides a new and powerful
perspective on a classical problem. In fact, the possibility of reformulating
the Zakharov equations for a degenerate quartet in terms of only two auxiliary
variables is intimately connected to the existence of an exact solution for
this configuration, which can be given in terms of elliptic functions (see
Shemer and Stiassnie (1985)). With subsequent normalisations, it is possible
to effectively describe the entire dynamics in terms of the total wave action
$A$ and a single bifurcation parameter $p,$ the mode separation.

The mode separation governs the appearance of fixed points in the phase plane.
Saddle points at $\eta=1$ and $\eta=0$ correspond to cases where the entire
energy is concentrated in one or two Fourier modes, respectively; the centres
correspond to a degenerate quartet of nearly resonant waves which undergoes no
time evolution. These fixed points are simply another case of the steady-state
waves recently discovered by Liao et al. (2016); Yang et al. (2022). Our
phase-plane analysis provides a simple way to obtain these solutions, and
highlights for the first time their critical role in organising the overall
dynamics.

The existence of fixed-points on the boundaries $\eta=1$ and $\eta=0$ of our
phase-space is intimately connected with the stability of uniform and
bichromatic wave trains, respectively. When such fixed points are absent, the
trajectories progress along nearly horizontal level lines on the cylindrical
phase-plane: the dynamic phase changes, but there is little redistribution of
energy. Loosely speaking, we may say that the trajectories do not need to
avoid the separatrices connecting the fixed points. The stability results we
obtain for uniform wave-trains demonstrate the robustness of the Benjamin-Feir
condition, which we recover from the fully nonlinear system without any small-
amplitude assumptions. For bichromatic wave trains our setting appears to be
different from those previously studied, and no comparable linear results
appear to be available.

We also obtain numerous special solutions – including primitive breathers –
which are identified with heteroclinic orbits in phase space. While the
corresponding solutions of the nonlinear Schrödinger equation involve a
continuum of modes, and therefore exhibit quantitatively different features,
the two-dimensional phase space allows for an incisive and transparent
analysis. It should be noted that, although breathers have been identified by
Ablowitz and Herbst (1990) with homoclinic orbits, this terminology obscures
the phase shift which is apparent in the heteroclinic connections between the
saddle points.

The primitive breather solutions we obtain include an analogue of the
Akhmediev breather, which arises from a plane-wave background, as well as
limiting cases in which a uniform wave train becomes (as $t\rightarrow\infty$)
a bichromatic wave train, and vice versa. While primitive Akhmediev breathers
exist in large ranges of parameter space (and thus for many configurations of
total energy and mode separation), the limiting (0,1) and (1,0) breathers are
unique for given $A$ or $p.$ We also find a seemingly new type of breather-
like solution, which arises from, and returns to, a bichromatic background.
The maximal energy exchange for this (0,0)-breather is therefore a
“demodulation”, in which the envelope flattens out.

##  Appendix A Interaction kernels

The interaction kernels appearing in the Zakharov equation (3) are lengthy,
and expressions for these can be found in Mei et al. (2018) or Krasitskii
(1994). For the case of unidirectional waves in deep water ($k_{i}>0$ for all
$i$) treated in the present study, considerable simplifications to the kernels
are possible.

$\displaystyle T(k,k,k,k)$ $\displaystyle=\frac{k^{3}}{4\pi^{2}}$ (37)
$\displaystyle T(k_{a},k_{b},k_{a},k_{b})$
$\displaystyle=\frac{1}{4\pi^{2}}k_{a}k_{b}\min(k_{a},k_{b})$ (38)
$\displaystyle T(k_{a},k_{b},k_{c},k_{d})$
$\displaystyle=\frac{(k_{a}k_{b}k_{c}k_{d})^{1/4}}{32\pi^{2}}\left[\left(k_{a}k_{b}\right)^{1/2}+\left(k_{c}k_{d}\right)^{1/2}\right]$
$\displaystyle\cdot\left(k_{a}+k_{b}+k_{c}+k_{d}-\left[|k_{a}-k_{c}|+|k_{a}-k_{d}|+|k_{b}-k_{c}|+|k_{b}-k_{d}|\right]\right)$
(39)

The equation for $T_{abcd}$ is taken from Kachulin et al. (2019) and adjusted
by a factor of $2\pi$. The kernels can be further simplified by making use of
the homogeneity property $T(\alpha k_{a},\alpha k_{b},\alpha k_{c},\alpha
k_{d})=\alpha^{3}T(k_{a},k_{b},k_{c},k_{d}).$ For a degenerate quartet of
waves $2k_{a}=k_{b}+k_{c}$ this enables us to write the kernel

$T(k_{a},k_{a},k_{b},k_{c})=\frac{k_{a}^{3}}{8\pi^{2}}\left[(1-q)(1-q^{2})^{1/4}(1+\sqrt{1-q^{2}})\right],$
(40)

where $q=p/|k_{a}|$ and $p$ is the mode separation distance, such that
$k_{b}=k_{a}+p,\,k_{c}=k_{a}-p.$

## References

  * Ablowitz and Herbst (1990) M. J. Ablowitz and B. M. Herbst.  On homoclinic structure and numerically induced chaos for the nonlinear Schrodinger equation.  _SIAM J. Appl. Math._ , 50(2):339–351, 1990.
  * Akhmediev et al. (1987) N. N. Akhmediev, V. M. Eleonskiĭ, and N. E. Kulagin.  First-order exact solutions of the nonlinear Schrödinger equation.  _Teoret. Mat. Fiz._ , 72(2):183–196, 1987.
  * Badulin et al. (1995) S. I. Badulin, V. I. Shrira, C. Kharif, and M. Ioualalen.  On two approaches to the problem of instability of short-crested water waves.  _J. Fluid Mech._ , 303(II):297–326, 1995.
  * Benjamin and Feir (1967) T. B. Benjamin and J. E. Feir.  The disintegration of wave trains on deep water Part 1. Theory.  _J. Fluid Mech._ , 27(03):417–430, 1967.
  * Cappellini and Trillo (1991) G. Cappellini and S. Trillo.  Third-order three-wave mixing in single-mode fibers: exact solutions and spatial instability effects.  _J. Opt. Soc. Am. B_ , 8(4):824, 1991.
  * Chabchoub et al. (2014) A. Chabchoub, B. Kibler, J. M. Dudley, and N. Akhmediev.  Hydrodynamics of periodic breathers.  _Philos. Trans. R. Soc. A Math. Phys. Eng. Sci._ , 372(2027), 2014.
  * Craik (1986) A. D. D. Craik.  Exact solutions of non-conservative equations for three-wave and second-harmonic resonance.  _Proc. R. Soc. London A_ , 406(July):1–12, 1986\.
  * Crawford et al. (1981) D. R. Crawford, B. M. Lake, P. G. Saffman, and H. C. Yuen.  Stability of weakly nonlinear deep-water waves in two and three dimensions.  _J. Fluid Mech._ , 105:177–191, 1981.
  * Davey and Stewartson (1974) A. Davey and K. Stewartson.  On three-dimensional packets of surface waves.  _Proc. R. Soc. A Math. Phys. Eng. Sci._ , 338:101–110, 1974\.
  * Dyachenko et al. (2017) A. I. Dyachenko, D. I. Kachulin, and V. E. Zakharov.  Super compact equation for water waves.  _J. Fluid Mech._ , 828:661–679, 2017.
  * Dysthe (1979) K. B. Dysthe.  Note on a modification to the nonlinear Schrodinger equation for application to deep water waves.  _Proc. R. Soc. A Math. Phys. Eng. Sci._ , 369:105–114, 1979\.
  * Dysthe and Trulsen (1999) K. B. Dysthe and K. Trulsen.  Note on Breather Type Solutions of the NLS as Models for Freak-Waves.  _Phys. Scr. T_ , 82:48–52, 1999.
  * Gramstad and Trulsen (2011) O. Gramstad and K. Trulsen.  Hamiltonian form of the modified nonlinear Schrödinger equation for gravity waves on arbitrary depth.  _J. Fluid Mech._ , 670:404–426, 2011.
  * Hasselmann (1962) K. Hasselmann.  On the non-linear energy transfer in a gravity-wave spectrum Part 1. General theory.  _J. Fluid Mech._ , 12(04):481, 1962.
  * Houtani et al. (2022) H. Houtani, H. Sawada, and T. Waseda.  Phase convergence and crest enhancement of modulated wave trains.  _Fluids_ , 7(8):275, 2022.
  * Ioualalen and Kharif (1994) M. Ioualalen and C. Kharif.  On the subharmonic instabilities of steady three-dimensional deep water waves.  _J. Fluid Mech._ , 262(II):265–291, 1994.
  * Kachulin et al. (2019) D. Kachulin, A. Dyachenko, and A. Gelash.  Interactions of coherent structures on the surface of deep water.  _Fluids_ , 4(2):1–21, 2019.
  * Krasitskii (1994) V. P. Krasitskii.  On reduced equations in the Hamiltonian theory of weakly nonlinear surface waves.  _J. Fluid Mech._ , 272:1–20, 1994.
  * Leblanc (2009) S. Leblanc.  Stability of bichromatic gravity waves on deep water.  _Eur. J. Mech. B/Fluids_ , 28(5):605–612, 2009\.
  * Liao et al. (2016) S. Liao, D. Xu, and M. Stiassnie.  On the steady-state nearly resonant waves.  _J. Fluid Mech._ , 794:175–199, 2016.
  * Liu et al. (2021) S. Liu, T. Waseda, and X. Zhang.  Phase locking phenomenon in the modulational instability of surface gravity waves.  In _Proceedings of the 36th International Workshop on Water Waves and Floating Bodies (IWWWFB), Seoul, Korea_ , pages 25–28, 2021.
  * Longuet-Higgins and Phillips (1962) M. S. Longuet-Higgins and O. M. Phillips.  Phase velocity effects in tertiary wave interactions.  _J. Fluid Mech._ , 12(03):333–336, 1962.
  * Mei et al. (2018) C. C. Mei, M. A. Stiassnie, and D. K.-P. Yue.  _Theory and applications of ocean surface waves_.  World Scientific Publishing Co., 3rd edition, 2018.
  * Onorato and Suret (2016) M. Onorato and P. Suret.  Twenty years of progresses in oceanic rogue waves: the role played by weakly nonlinear models.  _Nat. Hazards_ , 84(2):541–548, 2016.
  * Phillips (1960) O. M. Phillips.  On the dynamics of unsteady gravity waves of finite amplitude Part 1\. The elementary interactions.  _J. Fluid Mech._ , 9(2):193–217, 1960.
  * Shemer and Stiassnie (1985) L. Shemer and M. Stiassnie.  Initial instability and long-time evolution of Stokes waves.  In Y. Toba &. H. Mitsuyasu, editors, _Ocean Surf. wave Break. Turbul. Mix. radio probing_ , pages 51–57. 1985.
  * Slunyaev and Shrira (2013) A. V. Slunyaev and V. I. Shrira.  On the highest non-breaking wave in a group: Fully nonlinear water wave breathers versus weakly nonlinear theory.  _J. Fluid Mech._ , 735:203–248, 2013.
  * Stiassnie (1984) M. Stiassnie.  Note on the modified nonlinear Schrödinger equation for deep water waves.  _Wave Motion_ , 6(4):431–433, 1984.
  * Stokes (1847) G. G. Stokes.  On the theory of oscillatory waves.  _Trans. Camb. Phil. Soc_ , 8:441–455, 1847.
  * Stuhlmeier and Stiassnie (2019) R. Stuhlmeier and M. Stiassnie.  Nonlinear dispersion for ocean surface waves.  _J. Fluid Mech._ , 859:49–58, 2019.
  * Yang et al. (2022) X. Yang, J. Yang, and Z. Liu.  On the steady-state exactly resonant, nearly resonant, and non-resonant waves and their relationships.  _Physics of Fluids_ , 34(8):082107, 2022.
  * Yuen and Lake (1982) H. C. Yuen and B. M. Lake.  Nonlinear Dynamics of Deep-Water Gravity Waves.  In _Adv. Appl. Mech._ , pages 68–229. Academic Press, 1982.
  * Zakharov (1968) V. Zakharov.  Stability of periodic waves of finite amplitude on the surface of a deep fluid.  _J. Appl. Mech. Tech. Phys._ , 9(2):190–194, 1968\.
  * Zakharov et al. (1992) V. E. Zakharov, V. S. L’vov, and G. Falkovich.  _Kolmogorov Spectra of Turbulence I_.  Springer-Verlag Berlin Heidelberg, 1992.


# Analytical solution of the fractional linear time-delay systems and their
Ulam-Hyers stability

Nazim I. Mahmudov
Eastern Mediterranean University
Department of Mathematics
Famagusta, 99628 T. R. Northen Cyprus
Mersin 10 Turkey

(20 March 2021)

###### Abstract

We introduce the delayed Mittag-Leffler type matrix functions, delayed
fractional cosine, delayed fractional sine and use the Laplace transform to
obtain an analytical solution to the IVP for a Hilfer type fractional linear
time-delay system $D_{0,t}^{\mu,\nu}z\left(t\right)+Az\left(t\right)+\Omega
z\left(t-h\right)=f\left(t\right)$ of order $1<\mu<2$ and type $0\leq\nu\leq
1,$ with nonpermutable matrices $A$ and $\Omega$. Moreover, we study Ulam-
Hyers stability of the Hilfer type fractional linear time-delay system.
Obtained results extend those for Caputo and Riemann-Liouville type fractional
linear time-delay systems and new even for these fractional delay systems.

##  1 Introduction and auxiliary lemmas

Khusainov et al. [12] studied the following Cauchy problem for a second order
linear differential equation with pure delay:

$\left\\{\begin{array}[c]{c}x^{\prime\prime}\left(t\right)+\Omega^{2}x\left(t-\tau\right)=f\left(t\right),\
\ t\geq 0,\ \tau>0,\\\ x\left(t\right)=\varphi\left(t\right),\ \
x^{\prime}\left(t\right)=\varphi^{\prime}\left(t\right),\ \ -\tau\leq t\leq
0,\end{array}\right.$ (1)

where $f:\left[0,\infty\right)\rightarrow\mathbb{R}^{n}$, $\Omega$ is a
$n\times n$ nonsingular matrix, $\tau$ is the time delay and $\varphi$ is an
arbitrary twice continuously differentiable vector function. A solution of (1)
has an explicit representation of the form [12, Theorem 2]:

$\displaystyle x\left(t\right)$ $\displaystyle=\left(\cos_{\tau}\Omega
t\right)\varphi\left(-\tau\right)+\Omega^{-1}\left(\sin_{\tau}\Omega
t\right)\varphi^{\prime}\left(-\tau\right)$
$\displaystyle+\Omega^{-1}\int_{-\tau}^{0}\sin_{\tau}\Omega\left(t-\tau-s\right)\varphi^{\prime\prime}\left(s\right)ds$
$\displaystyle+\Omega^{-1}\int_{0}^{t}\sin_{\tau}\Omega\left(t-\tau-s\right)f\left(s\right)ds,$

where $\cos_{\tau}\Omega:\mathbb{R}\rightarrow\mathbb{R}^{n\times n}$ and
$\sin_{\tau}\Omega:\mathbb{R}\rightarrow\mathbb{R}^{n\times n}$ denote the
delayed matrix cosine of polynomial degree $2k$ on the intervals
$\left(k-1\right)\tau\leq t<k\tau$ and the delayed matrix sine of polynomial
degree $2k+1$ on the intervals $\left(k-1\right)\tau\leq t<k\tau$,
respectively.

It should be stressed out that pioneer works [12], [13] led to many new
results in integer and noninteger order time-delay differential equations and
discrete delayed system; see [1]-[11], [14], [15], [18]-[21].

Introducing the fractional analogue delayed matrices cosine/sine of a
polynomial degree, see formulas (3) and (4), Liang et al. [16] gave
representation of a solution to the initial value problem (2):

######  Theorem 1.

[16] Let $h>0$, $\varphi\in C^{2}\left(\left[-h,0\right],R^{n}\right)$,
$\Omega$ be a nonsingular $n\times n$ matrix, and let
$f:\left[0,\infty\right)\rightarrow R^{n}$ be a given function. The solution
$x:\left[-h,\infty\right)\rightarrow R^{n}$ of the intial value problem

$\left\\{\begin{array}[c]{c}{}^{C}D_{-h}^{\alpha}\left({}^{C}D_{-h}^{\alpha}\right)x\left(t\right)+\Omega^{2}x\left(t-h\right)=0,\
t\geq 0,\ h>0,\\\ x\left(t\right)=\varphi\left(t\right),\ \
x^{\prime}\left(t\right)=\varphi^{\prime}\left(t\right),\ \ -h\leq t\leq
0,\end{array}\right.$ (2)

has the form

$\displaystyle x\left(t\right)$ $\displaystyle=\left(\cos_{h,\alpha}\Omega
t^{\alpha}\right)\varphi\left(-h\right)+\Omega^{-1}\left(\sin_{h,\alpha}\Omega\left(t-h\right)^{\alpha}\right)\varphi^{\prime}\left(0\right)$
$\displaystyle+\Omega^{-1}\int_{-h}^{0}\cos_{h,\alpha}\Omega\left(t-h-s\right)^{\alpha}\varphi^{\prime}\left(s\right)ds$

where $\cos_{h,\alpha}\Omega t^{\alpha}$ is the fractional delayed matrix
cosine of a polynomial of degree $2k\alpha$ on the intervals
$\left(k-1\right)h\leq t<kh$, $\sin_{h,\alpha}\Omega t^{\alpha}$ is the
fractional delayed matrix sine of a polynomial of degree
$\left(2k+1\right)\alpha$ on the intervals $\left(k-1\right)h\leq t<kh$
defined as follows

$\cos_{h,\alpha}\Omega
t^{\alpha}:=\left\\{\begin{tabular}[c]{ll}$\Theta,$&$-\infty<t<-h,$\\\
$I,$&$-h\leq t<0,$\\\ &\\\
$I-\Omega^{2}\dfrac{t^{2\alpha}}{\Gamma\left(2\alpha+1\right)}+...+\left(-1\right)^{k}\Omega^{2k}\dfrac{\left(t-\left(k-1\right)h\right)^{2k\alpha}}{\Gamma\left(2k\alpha+1\right)},$&$\left(k-1\right)h\leq
t<kh,$\end{tabular}\ \ \right.$ (3) $\sin_{h,\alpha}\Omega
t^{\alpha}:=\left\\{\begin{tabular}[c]{ll}$\Theta,$&$-\infty<t<-h,$\\\
$\Omega\dfrac{\left(t+h\right)^{\alpha}}{\Gamma\left(\alpha+1\right)}$
$,$&$-h\leq t<0,$\\\ &\\\
$\Omega\dfrac{\left(t+h\right)^{\alpha}}{\Gamma\left(\alpha+1\right)}-\Omega^{3}\dfrac{t^{3\alpha}}{\Gamma\left(3\alpha+1\right)}$&\\\
$+...+\left(-1\right)^{k}\Omega^{2k+1}\dfrac{\left(t-\left(k-1\right)h\right)^{\left(2k+1\right)\alpha}}{\Gamma\left(\left(2k+1\right)\alpha+1\right)},$&$\left(k-1\right)h\leq
t<kh,$\end{tabular}\ \ \ \right.$ (4)

and $I$ is the identity matrix and $\Theta$ is the null matrix.

Mahmudov in [17] studied the following R-L linear fractional differential
delay equation of order $1<2\alpha\leq 2$ by introducing the concept of
fractional delayed matrix cosine
$\cos_{h,\alpha,\beta}\left\\{A,\Omega;t\right\\}$ and sine
$\sin_{h,\alpha,\beta}\left\\{A,\Omega;t\right\\}$ [17, Definitions 2 and 3].

$\left\\{\begin{array}[c]{c}{}^{RL}D_{-h+}^{\alpha}\left({}^{RL}D_{-h+}^{\alpha}\right)x\left(t\right)+A^{2}x\left(t\right)+\Omega^{2}x\left(t-h\right)=f\left(t\right),\
t\geq 0,\ h>0,\\\ x\left(t\right)=\varphi\left(t\right),\ \ \
\left(I_{-h^{+}}^{1-\alpha}x\right)\left(-h^{+}\right)=\varphi\left(-h\right),\
-h\leq t\leq 0,\\\ \ {}^{RL}D_{-h+}^{\alpha}x\left(t\right)=\
^{RL}D_{-h+}^{\alpha}\varphi\left(t\right),\\\
\left(I_{-h^{+}}^{1-\alpha}\left(D_{-h+}^{\alpha}x\right)\right)\left(-h^{+}\right)=\
\ ^{RL}D_{-h+}^{\alpha}\varphi\left(-h\right),\ -h\leq t\leq
0,\end{array}\right.$ (5)

where ${}^{RL}D_{-h+}^{\alpha}$ stands for the R-L fractional derivative of
order $0<\alpha\leq 1$ with lower limit $-h$, $A,\Omega\in\mathbb{R}^{n\times
n}$, $f\in C\left(\left[0,\infty\right),\mathbb{R}^{n}\right)$, $\varphi\in
C^{1}\left(\left[-h,0\right],\mathbb{R}^{n}\right)$. Obviously, the derivative
can be started at $-h$ instead of $0$, since the function $x\left(t\right)$
governed by (5) actually originates at $-h.$ However, as is known, changing
the starting point of the derivative modifies the derivative and leads to a
different problem. In this article, we study the case when the derivative
started at $0$.

We study the following Hilfer type linear fractional differential time-delay
equation of order $1<\mu<2$ and type $0\leq\nu\leq 1$:

$\left\\{\begin{array}[c]{c}D_{0,t}^{\mu,\nu}z\left(t\right)+Az\left(t\right)+\Omega
z\left(t-h\right)=f\left(t\right),\ \ \ t\in\left[0,T\right],\\\
z\left(t\right)=\varphi\left(t\right),\ \ -h\leq t\leq 0,\\\
\left.D_{0,t}^{-\left(2-\mu\right)\left(1-\nu\right)+1,\nu}z\left(t\right)\right|_{t=0}=a_{1},\\\
\left.D_{0,t}^{-\left(2-\mu\right)\left(1-\nu\right),\nu}z\left(t\right)\right|_{t=0}=a_{0},\end{array}\right.$
(6)

where $D_{0,t}^{\mu,\nu}$ stands for the Hilfer fractional derivative of order
$1<\mu<2$ and type $0\leq\nu\leq 1$ with lower limit $0$,
$A,\Omega\in\mathbb{R}^{d\times d}$, $f\in
C\left(\left[0,T\right],\mathbb{R}^{d}\right)$, $\varphi\in
C^{1}\left(\left[-h,0\right],\mathbb{R}^{d}\right).$

Delayed perturbation of Mittag-Leffler matrix functions serves as a suitable
tool for solving linear fractional continuous time-delay equations. First, the
delayed matrix exponential function (delayed matrix Mittag-Leffler function)
was defined to solve linear purely delayed (fractional) systems of order one.
Then, the second order differential systems with pure delay were considered
and suitable delayed sine and cosine matrix functions were introduced in [12].
Later Liang et al. [16] introduced the fractional analogue of delayed
cosine/sine matrices and obtained an explicit solution of the sequential
fractional Caputo type equations with pure delay, the case $A=\Theta$ (zero
matrix). Recently, Mahmudov [17] introduced the fractional analogue of delayed
matrices cosine/sine in the case when $A$ and $\Omega$ commutes to solve the
sequential Riemann-Liouville type linear time-delay system. It should be
noticed that the model investigated here is not sequential and differs from
that of discussed in [16], [17]. For the sake of completeness, we also refer
to studies of discrete/continuous variants of delayed matrices used to obtain
exact solution to linear difference equations with delays [1]-[21].

We introduce a concept of delayed Mittag-Leffler type matrix function of two
parameters:

######  Definition 2.

Delayed Mittag-Leffler type matrix function of two parameters
$Y_{\mu,\gamma}^{h}:\left[0,\infty\right)\rightarrow\mathbb{R}^{d}$ is defined
as follows:

$Y_{\mu,\gamma}^{h}\left(A,\Omega;t\right)=Y_{\mu,\gamma}^{h}\left(t\right):={\displaystyle\sum_{m=0}^{\infty}}{\displaystyle\sum\limits_{k=0}^{\infty}}\left(-1\right)^{k}Q_{k,m}^{A,\Omega}\dfrac{\left(t-mh\right)_{+}^{k\mu+\gamma-1}}{\Gamma\left(k\mu+\gamma\right)},$

where $\left(t\right)_{+}=\max\left\\{0,t\right\\}$ and

$\displaystyle Q_{k,m}$
$\displaystyle=Q_{k,m}^{A,\Omega}={\displaystyle\sum\limits_{j=m}^{k}}A^{k-j}\Omega
Q_{j-1,m-1}^{A,\Omega},\ \ .$ (7) $\displaystyle Q_{0,m}^{A,\Omega}$
$\displaystyle=Q_{k,-1}^{A,\Omega}=\Theta,\ \ Q_{k,0}^{A,\Omega}=A^{k}\ \
Q_{0,0}^{A,\Omega}=I,\ \ k=0,1,2,...,m=0,1,2,...$

In this definition $Q_{k,m}^{A,\Omega}$ plays the role of a kernel, see [17],
[20]. It is clear that

$\displaystyle Q_{k+1,m}^{A,\Omega}$ $\displaystyle=AQ_{k,m}^{A,\Omega}+\Omega
Q_{k,m-1}^{A,\Omega},$ $\displaystyle Q_{0,m}^{A,\Omega}$
$\displaystyle=Q_{k,-1}^{A,\Omega}=\Theta,\ \ Q_{0,0}^{A,\Omega}=I,$
$\displaystyle k$ $\displaystyle=0,1,2,...,m=0,1,2,...$

We state the main novelties of our article as below:

  * •

We introduce a novel delayed Mittag-Leffler type matrix function
$Y_{\mu,\gamma}^{h}\left(A,\Omega;t\right)$.

  * •

If $\Omega=\Theta,$ then

$\displaystyle Y_{2,1}^{h}\left(A^{2},\Theta;t\right)$
$\displaystyle={\displaystyle\sum\limits_{k=0}^{\infty}}\left(-1\right)^{k}A^{2k}\dfrac{t^{2k}}{\left(2k\right)!}=\cos\left(At\right),\
\ \ $ $\displaystyle AY_{2,2}^{h}\left(A^{2},\Theta;t\right)$
$\displaystyle=A{\displaystyle\sum\limits_{k=0}^{\infty}}\left(-1\right)^{k}A^{2k}\dfrac{t^{2k+1}}{\left(2k+1\right)!}=\sin\left(At\right),$
$\displaystyle Y_{\mu,1}^{h}\left(A^{2},\Theta;t\right)$
$\displaystyle={\displaystyle\sum\limits_{k=0}^{\infty}}\left(-1\right)^{k}A^{2k}\dfrac{t^{k\mu}}{\Gamma\left(k\mu+1\right)},\
\ \ $ $\displaystyle Y_{\mu,2}^{h}\left(A^{2},\Theta;t\right)$
$\displaystyle={\displaystyle\sum\limits_{k=0}^{\infty}}\left(-1\right)^{k}A^{2k}\dfrac{t^{k\mu+1}}{\Gamma\left(k\mu+2\right)},\
1<\mu<2.$

$Y_{\mu,1}^{h}\left(A^{2},\Theta;t\right)$ and
$Y_{\mu,2}^{h}\left(A^{2},\Theta;t\right)$ can be called fractional cosine and
sine for $1<\mu<2.$ Similar cosine/sine matrix functions were defined in [16],
[19] to solve $1<2\mu<2$ order sequential fractional differential equations.

  * •

If $A=\Theta$ then we have

$Q_{m,m}^{A,\Omega}=\Omega^{m},\ \ \ \
Y_{\mu,\gamma}^{h}\left(A,\Omega;t\right)={\displaystyle\sum_{m=0}^{\infty}}\left(-1\right)^{m}\Omega^{m}\dfrac{\left(t-mh\right)_{+}^{m\mu+\gamma-1}}{\Gamma\left(m\mu+\gamma\right)}.$

Moreover,

$\displaystyle Y_{\mu,1}^{h}\left(\Theta,\Omega^{2};t\right)$
$\displaystyle={\displaystyle\sum_{m=0}^{\infty}}\left(-1\right)^{m}\Omega^{2m}\dfrac{\left(t-mh\right)_{+}^{m\mu}}{\Gamma\left(m\mu+1\right)}=\cos_{\mu}^{h}\left(\Omega;t\right),\
$ $\displaystyle\Omega Y_{\mu,2}^{h}\left(\Theta,\Omega^{2};t\right)$
$\displaystyle=\Omega{\displaystyle\sum_{m=0}^{\infty}}\left(-1\right)^{m}\Omega^{2m}\dfrac{\left(t-mh\right)_{+}^{m\mu}}{\Gamma\left(m\mu+2\right)}=\Omega\sin_{\mu}^{h}\left(\Omega;t\right).$

Similar delayed cosine/sine matrix functions were defined in [2], [19] to
solve $1<2\mu<2$ order sequential fractional linear differential equations
with pure delay.

  * •

We give an exact analytical solution of the Hilfer type fractional problem (6)
using delayed Mittag-Leffler type matrix function
$Y_{\mu,\gamma}^{h}\left(A,\Omega;t\right)$ and study their Ulam-Hyers
stability. Obtained results are new even for Caputo and Riemann-Liouville type
fractional linear time-delay systems.

  * •

Although the problem considered by us is fractional of order $1<\mu<2$ and
type $0\leq\nu\leq 1$, our approach is also applicable to the classical
second-order equations. Thus our results are new even for the classical second
order oscillatory system.

Before introducing properties of $Y_{\mu,\gamma}^{h}\left(A,\Omega;t\right),$
we recall the definition of the Hilfer fractional derivative and Ulam-Hyers
stability:

######  Definition 3.

Let $m\in\mathbb{N},$ $m-1<\mu<m$, $0\leq\nu\leq 1$, $a\in\mathbb{R}$, and
$f\in C^{m}\left[a,b\right]$. Then the Hilfer fractional derivative of $f$ of
order $\mu$ and type $\nu$ is given by

$D_{a,t}^{\mu,\nu}f\left(t\right):=I_{a,t}^{\nu\left(m-\mu\right)}\frac{d^{m}}{dt^{m}}I_{a,t}^{\left(1-\nu\right)\left(m-\mu\right)}f\left(t\right),$

where

$I_{a,t}^{\gamma}f\left(t\right):=\frac{1}{\Gamma\left(\gamma\right)}\int_{a}^{t}\left(t-s\right)^{\gamma-1}f\left(s\right)ds$

is the R-L fractional integral of $f$ of order $\gamma>0$.

The main tool we use in this paper is the Laplace transform
$F\left(s\right):=L\left\\{f\left(t\right)\right\\}=\int_{0}^{\infty}e^{-st}f\left(t\right)dt,$
$\operatorname{Re}s>a$, which is defined for an exponentially bounded function
$f.$ Here are some of properties of the Laplace transform.

######  Lemma 4.

The following equalities hold true for sufficiently large _Re_
$\left(s\right)$ and appropriate functions $f,g$:

  1. (i)

$L\left\\{af\left(t\right)+bg\left(t\right)\right\\}=aL\left\\{f\left(t\right)\right\\}+bL\left\\{g\left(t\right)\right\\},\
\ \ a,b\in\mathbb{R};$

  2. (ii)

$L^{-1}\left\\{e^{-sh}s^{-1}\right\\}=1,\ \ t\geq h\geq 0;$

  3. (iii)

$L^{-1}\left\\{F\left(s\right)G\left(s\right)\right\\}=\left(f\ast
g\right)\left(t\right);$

  4. (iv)

$L\left\\{D_{0,t}^{\mu,\nu}f\left(t\right)\right\\}=s^{\mu}L\left\\{f\left(t\right)\right\\}-{\displaystyle\sum\limits_{k=0}^{m-1}}s^{m\left(1-\nu\right)+\mu\nu-k-1}I_{0,t}^{\left(1-\nu\right)\left(m-\mu\right)-k}f\left(0\right);$

  5. (v)

$L^{-1}\left\\{1\right\\}=\delta\left(t\right),\ $where $\delta\left(t\right)$
is Dirac delta distribution;

  6. (vi)

$L^{-1}\left\\{e^{-nsh}s^{-n}\right\\}=\frac{\left(t-nh\right)_{+}^{n-1}}{\left(n-1\right)!},\
\ h>0$, $n\in\mathbb{N}$;

  7. (vii)

$L^{-1}\left\\{e^{-sh}F\left(s\right)\right\\}=f\left(t-h\right),\ \ h\geq 0$;

  8. (viii)

$L^{-1}\left\\{e^{-sh}s^{\alpha\gamma-\beta}\left(s^{\alpha}I-A\right)^{-\gamma}\right\\}=\left(t-h\right)^{\beta-1}E_{\alpha,\beta}^{\gamma}\left(A\left(t-h\right)^{\alpha}\right)$,
$t\geq h$, where
$E_{\alpha,\beta}^{\gamma}\left(z\right)={\displaystyle\sum\limits_{k=0}^{\infty}}\dfrac{z^{\alpha
k}}{\Gamma\left(\alpha k+\beta\right)}\dfrac{\left(\gamma\right)_{k}}{k!}$ is
the three parameter Mittag-Leffler function, $\alpha,\beta,\gamma>0$,
$t\in\mathbb{R}$ and
$\left(\gamma\right)_{k}:=\gamma\left(\gamma+1\right)...\left(\gamma+k-1\right).$

######  Definition 5.

System (6) is Ulam-Hyers stable on $\left[0,T\right]$ if there exists $C>0$
such that for any $\varepsilon>0$ and for any function
$z^{\ast}\left(t\right)$ satisfying inequality

$\left\|D_{0,t}^{\mu,\nu}z^{\ast}\left(t\right)+Az^{\ast}\left(t\right)+\Omega
z^{\ast}\left(t-h\right)-f\left(t\right)\right\|\leq\varepsilon$ (8)

and the initial conditions in (6), there is a solution $z\left(t\right)$of (6)
such that

$\left\|z^{\ast}\left(t\right)-z\left(t\right)\right\|\leq C\varepsilon$

for every $t\in\left[0,T\right].$

We reduce the notations of $Y_{\mu,\gamma}^{h}\left(A,\Omega;t\right)$,
$Q_{k,m}^{A,\Omega}$ to a mere $Y_{\mu,\gamma}^{h}\left(t\right)$, $Q_{k,m}$
in the sequel.

######  Theorem 6.

The following formulae hold:

  1. 1.

The function $Y_{\mu,\gamma}^{h}\left(\cdot\right)$ is continuous on
$\left(0,+\infty\right).$

  2. 2.

$\frac{d}{dt}Y_{\mu,\gamma+1}^{h}\left(t\right)=Y_{\mu,\gamma}^{h}\left(t\right),\
\
\frac{d}{dt}Y_{\mu,\gamma+2}^{h}\left(t\right)=Y_{\mu,\gamma+1}^{h}\left(t\right)$
for all $t\in\mathbb{R}$.

  3. 3.

$D_{0,t}^{\mu,\nu}Y_{\mu,\gamma}^{h}=-AY_{\mu,\gamma}^{h}\left(t\right)-\Omega
Y_{\mu,\gamma}^{h}\left(t-h\right).$

###### Proof.

The proofs of the properties $1.$ and $2.$ are obvious. Proof of the property
$3.$ is based on the following formula:

$D_{0,t}^{\mu,\nu}t^{\alpha}=\frac{\Gamma\left(\alpha+1\right)}{\Gamma\left(\alpha-\mu+1\right)}t^{\alpha-\mu},\
\ t>0,\ n-1<\mu\leq n,\ 0\leq\nu\leq 1,\ \alpha>-1.$

∎

The main tool we use in this paper is the Laplace transform
$F\left(s\right):=L\left\\{f\left(t\right)\right\\}=\int_{0}^{\infty}e^{-st}f\left(t\right)dt,$
$\operatorname{Re}s>a$, which is defined for an exponentially bounded function
$f.$

######  Lemma 7.

We have

$\displaystyle
L^{-1}\left\\{\left(e^{-hs}\left(s^{\mu}I+A\right)^{-1}\Omega\right)^{m}s^{\mu-\gamma}\left(s^{\mu}I+A\right)^{-1}\right\\}$
$\displaystyle={\displaystyle\sum\limits_{k=0}^{\infty}}\left(-1\right)^{k-m}Q_{k,m}\dfrac{\left(t-mh\right)_{+}^{k\mu+\gamma-1}}{\Gamma\left(k\mu+\gamma\right)},$

where $Q_{k,m}$ is defined in (7).

###### Proof.

For $n=0$ by Lemma 4(viii) we have

$\displaystyle
L^{-1}\left\\{s^{\mu-\gamma}\left(s^{\mu}I+A\right)^{-1}\right\\}$
$\displaystyle=t^{\gamma-1}E_{\mu,\gamma}\left(-At^{\mu}\right),\ \ \ $
$\displaystyle
L^{-1}\left\\{e^{-sh}\left(s^{\mu}I+A\right)^{-\gamma}\right\\}$
$\displaystyle=\left(t-h\right)_{+}^{\mu-1}E_{\mu,\mu}\left(-A\left(t-h\right)^{\mu}\right),\
\ \ t\geq h.$

Let $Q_{k,0}=A^{k}$. For $n=1$, we use the convolution property (Lemma 4(iii))
of the Laplace transform to get

$\displaystyle L^{-1}\left\\{e^{-hs}\left(s^{\mu}I+A\right)^{-1}\Omega
s^{\mu-\gamma}\left(s^{\mu}I+A\right)^{-1}\right\\}$
$\displaystyle=L^{-1}\left\\{e^{-hs}\left(s^{\mu}I+A\right)^{-1}\Omega\right\\}\ast
L^{-1}\left\\{s^{\mu-\gamma}\left(s^{\mu}I+A\right)^{-1}\right\\}$
$\displaystyle=\int_{0}^{t}\left(s-h\right)_{+}^{\mu-1}E_{\mu,\mu}\left(-A\left(s-h\right)^{\mu}\right)\Omega\left(t-s\right)^{\gamma-1}E_{\mu,\gamma}\left(-A\left(t-s\right)^{\mu}\right)ds$
$\displaystyle={\displaystyle\sum\limits_{k=0}^{\infty}}{\displaystyle\sum\limits_{j=0}^{\infty}}\frac{\left(-1\right)^{k}A^{k}\Omega\left(-1\right)^{j}A^{j}}{\Gamma\left(\mu
k+\mu\right)\Gamma\left(\mu j+\gamma\right)}\int_{h}^{t}\left(s-h\right)^{\mu
k+\mu-1}\left(t-s\right)^{\mu j+\gamma-1}ds$
$\displaystyle={\displaystyle\sum\limits_{k=0}^{\infty}}{\displaystyle\sum\limits_{j=0}^{\infty}}\left(-1\right)^{k}A^{k}\Omega\left(-1\right)^{j}A^{j}\frac{\left(t-h\right)_{+}^{\mu
k+\mu j+\mu+\gamma-1}}{\Gamma\left(\mu k+\mu j+\mu+\gamma\right)}$
$\displaystyle={\displaystyle\sum\limits_{k=0}^{\infty}}\left(-1\right)^{k}{\displaystyle\sum\limits_{j=0}^{k}}A^{k-j}\Omega
A^{j}\frac{\left(t-h\right)_{+}^{\mu k+\mu+\gamma-1}}{\Gamma\left(\mu
k+\mu+\gamma\right)}$
$\displaystyle={\displaystyle\sum\limits_{k=1}^{\infty}}\left(-1\right)^{k-1}{\displaystyle\sum\limits_{j=0}^{k-1}}A^{k-1-j}\Omega
A^{j}\frac{\left(t-h\right)_{+}^{\mu k+\gamma-1}}{\Gamma\left(\mu
k+\gamma\right)}.$

Now, to use the mathematical induction, suppose that it holds for $n=m$. Then
convolution property yields

$\displaystyle
L^{-1}\left\\{\left(e^{-hs}\left(s^{\mu}I+A\right)^{-1}\Omega\right)^{m+1}s^{\mu-\gamma}\left(s^{\mu}I+A\right)^{-1}\right\\}$
$\displaystyle=L^{-1}\left\\{e^{-hs}\left(s^{\mu}I+A\right)^{-1}\Omega\right\\}\ast
L^{-1}\left\\{\left(e^{-hs}s^{-\beta}\left(s^{\mu}I+A\right)^{-1}\Omega\right)^{m}s^{\mu-\gamma}\left(s^{\mu}I+A\right)^{-1}\right\\}$
$\displaystyle=\int_{h}^{t}\left(s-h\right)_{+}^{\mu-1}E_{\mu,\mu}\left(-A\left(s-h\right)^{\mu}\right)\Omega{\displaystyle\sum\limits_{j=0}^{\infty}}\left(-1\right)^{j}Q_{j+m,m}\frac{\left(t-s-
mh\right)_{+}^{\mu j+\mu m+\gamma-1}}{\Gamma\left(\mu j+\mu m+\mu\right)}ds$
$\displaystyle={\displaystyle\sum\limits_{k=0}^{\infty}}{\displaystyle\sum\limits_{j=0}^{\infty}}\left(-1\right)^{k}A^{k}\Omega\left(-1\right)^{j}Q_{j+m,m}\int\limits_{h}^{t}\dfrac{\left(t-s-h\right)_{+}^{k\mu+\mu-1}}{\Gamma\left(k\mu+\mu\right)}\frac{\left(s-mh\right)_{+}^{\mu
j+\mu m+\gamma-1}}{\Gamma\left(\mu j+\mu m+\gamma\right)}ds$
$\displaystyle={\displaystyle\sum\limits_{k=0}^{\infty}}{\displaystyle\sum\limits_{j=0}^{\infty}}\left(-1\right)^{k}A^{k}\Omega\left(-1\right)^{j}Q_{j+m,m}\int\limits_{mh}^{t-h}\dfrac{\left(t-s-h\right)^{k\mu+\mu-1}}{\Gamma\left(k\mu+\mu\right)}\frac{\left(s-mh\right)^{\mu
j+\mu m+\gamma-1}}{\Gamma\left(\mu j+\mu m+\gamma\right)}ds$
$\displaystyle={\displaystyle\sum\limits_{k=0}^{\infty}}{\displaystyle\sum\limits_{j=0}^{\infty}}\left(-1\right)^{k}A^{k}\Omega\left(-1\right)^{j}Q_{j+m,m}\dfrac{\left(t-\left(m+1\right)h\right)_{+}^{k\mu+j\mu+\left(m+1\right)\mu+\gamma-1}}{\Gamma\left(k\mu+j\mu+\left(m+1\right)\mu+\gamma\right)}$
$\displaystyle={\displaystyle\sum\limits_{k=m+1}^{\infty}}\left(-1\right)^{k-m-1}{\displaystyle\sum\limits_{j=0}^{k-m-1}}A^{k-j}\Omega
Q_{j+m,m}\dfrac{\left(t-\left(m+1\right)h\right)_{+}^{k\mu+j\mu+\left(m+1\right)\mu+\gamma-1}}{\Gamma\left(k\mu+j\mu+\left(m+1\right)\mu+\gamma\right)}$

what was to be proved. ∎

######  Lemma 8.

We have

$\displaystyle Y_{\mu,\gamma}^{h}\left(t\right)$
$\displaystyle:=L^{-1}\left\\{s^{\mu-\gamma}\left(s^{\mu}I+A+\Omega
e^{-hs}\right)^{-1}\right\\}$
$\displaystyle={\displaystyle\sum_{m=0}^{\infty}}{\displaystyle\sum\limits_{k=0}^{\infty}}\left(-1\right)^{k}Q_{k,m}\dfrac{\left(t-mh\right)_{+}^{k\mu+\gamma-1}}{\Gamma\left(k\mu+\gamma\right)}.$

###### Proof.

It is easy to see that

$\displaystyle L^{-1}\left\\{s^{\mu-\gamma}\left(s^{\mu}I+A+\Omega
e^{-hs}\right)^{-1}\right\\}$
$\displaystyle=L^{-1}\left\\{s^{\mu-\gamma}\left(\left(s^{\mu}I+A\right)I+\left(s^{\mu}I+A\right)\left(s^{\mu}I+A\right)^{-1}\Omega
e^{-hs}\right)^{-1}\right\\}$
$\displaystyle=L^{-1}\left\\{\left(I+\left(s^{\mu}I+A\right)^{-1}\Omega
e^{-hs}\right)^{-1}s^{\mu-\gamma}\left(s^{\mu}I+A\right)^{-1}\right\\}$
$\displaystyle=L^{-1}\left\\{{\displaystyle\sum_{m=0}^{\infty}}e^{-mhs}\left(-1\right)^{n}\left(\left(s^{\mu}I+A\right)^{-1}\Omega\right)^{m}s^{\mu-\gamma}\left(s^{\mu}I+A\right)^{-1}\right\\}$
$\displaystyle={\displaystyle\sum_{m=0}^{\infty}}L^{-1}\left\\{e^{-mhs}\left(-1\right)^{m}\left(\left(s^{\mu}I+A\right)^{-1}\Omega\right)^{m}s^{\mu-\gamma}\left(s^{\mu}I+A\right)^{-1}\right\\}.$

Hence, by Lemma 7 we have

$L^{-1}\left\\{s^{\mu-\gamma}\left(s^{\mu}I+A+\Omega
e^{-hs}\right)^{-1}\right\\}={\displaystyle\sum_{m=0}^{\infty}}{\displaystyle\sum\limits_{k=0}^{\infty}}\left(-1\right)^{k}Q_{k,m}\dfrac{\left(t-mh\right)_{+}^{k\mu+\gamma-1}}{\Gamma\left(k\mu+\gamma\right)}.$

∎

##  2 Exact analytical solution and Ulam-Hyers stability

We obtain the exact analytical solution of the Hilfer type fractional second
order problem (6) using delayed Mittag-Leffler type matrix function
$Y_{\mu,\gamma}^{h}\left(A,\Omega;t\right)$ and study their Ulam-Hyers
stability.

######  Theorem 9.

The analytical solution of the IVP problem (6) has the form

$\displaystyle z\left(t\right)$
$\displaystyle=Y_{\mu,\left(\mu-2\right)\left(1-\nu\right)+1}^{h}\left(t\right)\left(I_{0,t}^{\left(1-\nu\right)\left(2-\mu\right)}\varphi\right)\left(0\right)$
$\displaystyle+Y_{\mu,\left(\mu-2\right)\left(1-\nu\right)+2}^{h}\left(t\right)\left(I_{0,t}^{\left(1-\nu\right)\left(2-\mu\right)-1}\varphi\right)\left(0\right)$
$\displaystyle-\int_{-h}^{0}Y_{\mu,\mu}^{h}\left(t-s-h\right)\Omega\varphi\left(s\right)ds+\int_{0}^{t}Y_{\mu,\mu}^{h}\left(t-s\right)f\left(s\right)ds.$

###### Proof.

Assume that the function $f$ and the solution of (6) is exponentially bounded.
By applying the Laplace transform to the both sides of (6), we obtain the
following relation

$L\left\\{D_{0,t}^{\mu,\nu}z\left(t\right)\right\\}+AL\left\\{z\left(t\right)\right\\}+\Omega
L\left\\{z\left(t-h\right)\right\\}=L\left\\{f\left(t\right)\right\\}.$

It follows that

$\displaystyle\left(s^{\mu}I+A+\Omega e^{-hs}\right)Z\left(s\right)$
$\displaystyle=s^{2\left(1-\nu\right)+\mu\nu-1}\left(I_{0,t}^{\left(1-\nu\right)\left(2-\mu\right)}\varphi\right)\left(0\right)$
$\displaystyle+s^{2\left(1-\nu\right)+\mu\nu-2}\left(I_{0,t}^{\left(1-\nu\right)\left(2-\mu\right)-1}\varphi\right)\left(0\right)$
$\displaystyle-\Omega\int_{0}^{\infty}e^{-st}z\left(t-h\right)dt+F\left(s\right),$

where $Z\left(s\right)=L\left\\{z\left(t\right)\right\\},\
F\left(s\right)=L\left\\{f\left(t\right)\right\\}$. For sufficiently large
$s$, such that

$\left\|A+\Omega e^{-hs}\right\|<s^{\mu},$

the matrix $s^{\mu}I+A+\Omega e^{-hs}$ is invertible and it holds that

$\displaystyle Z\left(s\right)$
$\displaystyle=s^{2\left(1-\nu\right)+\mu\nu-1}\left(s^{\mu}I+A+\Omega
e^{-hs}\right)^{-1}\left(I_{0,t}^{\left(1-\nu\right)\left(m-\mu\right)}\varphi\right)\left(0\right)$
$\displaystyle+s^{2\left(1-\nu\right)+\mu\nu-2}\left(s^{\mu}I+A+\Omega
e^{-hs}\right)^{-1}\left(I_{0,t}^{\left(1-\nu\right)\left(m-\mu\right)-1}\varphi\right)\left(0\right)$
$\displaystyle-\left(s^{\mu}I+A+\Omega
e^{-hs}\right)^{-1}\Omega\Psi\left(s\right)$
$\displaystyle+\left(s^{\mu}I+A+\Omega e^{-hs}\right)^{-1}F\left(s\right).$

By Lemma 8

$\displaystyle z\left(t\right)$
$\displaystyle=Y_{\mu,\left(\mu-2\right)\left(1-\nu\right)+1}^{h}\left(t\right)\left(I_{0,t}^{\left(1-\nu\right)\left(m-\mu\right)}\varphi\right)\left(0\right)$
$\displaystyle+Y_{\mu,\left(\mu-2\right)\left(1-\nu\right)+2}^{h}\left(t\right)\left(I_{0,t}^{\left(1-\nu\right)\left(m-\mu\right)-1}\varphi\right)\left(0\right)$
$\displaystyle-\int_{-h}^{0}Y_{\mu,\mu}^{h}\left(t-s-h\right)\Omega\varphi\left(s\right)ds+\int_{0}^{t}Y_{\mu,\mu}^{h}\left(t-s\right)f\left(s\right)ds,$
(9)

since

$\displaystyle L^{-1}\left\\{\left(s^{\mu}I+A+\Omega
e^{-hs}\right)^{-1}\Omega\Psi\left(s\right)\right\\}=L^{-1}\left\\{\left(s^{\mu}I+A+\Omega
e^{-hs}\right)^{-1}\right\\}\ast
L^{-1}\left\\{\Omega\Psi\left(s\right)\right\\}$
$\displaystyle=\int_{0}^{t}Y_{\mu,\mu}^{h}\left(t-s\right)\Omega\psi\left(s-h\right)ds=\int_{0}^{h}Y_{\mu,\mu}^{h}\left(t-s\right)\Omega\varphi\left(s-h\right)ds$
$\displaystyle=\int_{-h}^{0}Y_{\mu,\mu}^{h}\left(t-s-h\right)\Omega\varphi\left(s\right)ds.$

Now the assumption on the exponential boundedness can be omitted. We can
easily check that (9) is a solution of (6). ∎

######  Theorem 10.

Let $1<\mu<2$, $0\leq\nu\leq 1,$ $f\in
C\left(\left[0,\infty\right),\mathbb{R}^{d}\right).$ System (6) is stable in
Ulam-Hyers sence on $[0,T]$.

###### Proof.

Let $z^{\ast}\left(t\right)$ satisfy the inequality (8) and the initial
conditions in (6). Set

$X\left(t\right)=D_{0,t}^{\mu,\nu}z^{\ast}\left(t\right)+Az^{\ast}\left(t\right)+\Omega
z^{\ast}\left(t-h\right)-f\left(t\right),\ t\in\left[0,T\right].$

It follows from definition 5 that
$\left\|X\left(t\right)\right\|<\varepsilon$. By Theorem 9 we have

$\displaystyle z^{\ast}\left(t\right)$
$\displaystyle=Y_{\mu,\left(\mu-2\right)\left(1-\nu\right)+1}^{h}\left(t\right)\left(I_{0,t}^{\left(1-\nu\right)\left(m-\mu\right)}\varphi\right)\left(0\right)$
$\displaystyle+Y_{\mu,\left(\mu-2\right)\left(1-\nu\right)+2}^{h}\left(t\right)\left(I_{0,t}^{\left(1-\nu\right)\left(m-\mu\right)-1}\varphi\right)\left(0\right)$
$\displaystyle-\int_{-h}^{0}Y_{\mu,\mu}^{h}\left(t-s-h\right)\Omega\varphi\left(s\right)ds+\int_{0}^{t}Y_{\mu,\mu}^{h}\left(t-s\right)\left(f\left(s\right)-X\left(s\right)\right)ds.$

Thus we can estimate the difference $z^{\ast}\left(t\right)-z\left(t\right)$
as follows

$\left\|z^{\ast}\left(t\right)-z\left(t\right)\right\|=\left\|\int_{0}^{t}Y_{\mu,\mu}^{h}\left(t-s\right)X\left(s\right)ds\right\|\leq\varepsilon\int_{0}^{T}\left\|Y_{\mu,\mu}^{h}\left(T-s\right)\right\|ds=C\varepsilon.$

Then the problem (6) is Ulam-Hyers stable on $[0,T]$. ∎

##  3 Conclusion

The article solves a problem of finding exact analytical solution of
continuous linear time-delay systems using the delayed Mittag-Leffler type
matrix functions of two variables. In articles [12], [4] delayed exponential
is suggested to obtain an exact solution of delayed first order continuous
equations. Similar results for sequential Caputo type and Riemann-Liouville
type fractional linear time-delay systems of order $1<2\alpha<2$ were obtained
in [16], [17]. These results are obtained either for systems with pure delay
or under the condition of commutativity of $A$ and $\Omega$. In this article
we drop the commutativity condition. The result has been obtained by defining
the new delayed Mittag-Leffler matrix function and employing the Laplace
transform. The work contained in this article will be useful for future
research on fractional time-delay systems.

## References

  * [1] A. Boichuk, J. Diblík, D. Khusainov, and M. Ružicková, “Fredholm’s boundary-value problems for differential systems with a single delay,” Nonlinear Analysis: Theory, Methods and Applications, vol. 72, no. 5, pp. 2251–2258, 2010, doi: 10.1016/j.na.2009.10.025.
  * [2] B. Bonilla, M. Rivero, and J. J. Trujillo, “On systems of linear fractional differential equations with constant coefficients,” Applied Mathematics and Computation, vol. 187, no. 1, pp. 68–78, 2007, doi: 10.1016/j.amc.2006.08.104.
  * [3] J. Diblík, M. Fečkan, and M. Pospíšil, “Representation of a solution of the cauchy problem for an oscillating system with two delays and permutable matrices,” Ukrainian Mathematical Journal, vol. 65, no. 1, pp. 64–76, 2013, doi: 10.1007/s11253-013-0765-y.
  * [4] J. Diblík, D. Y. Khusainov, J. Lukáčová, and M. Ružicková, “Control of oscillating systems with a single delay,” Advances in Difference Equations, vol. 2010, pp. 1–15, 2010, doi: 10.1155/2010/108218.
  * [5] J. Diblík, D. Y. Khusainov, and M. Ružicková, “Controllability of linear discrete systems with constant coefficients and pure delay,” SIAM Journal on Control and Optimization, vol. 47, pp. 1140–1149, 2008, doi: 10.1137/070689085.
  * [6] J. Diblík and D. Y. Khusainov, “Representation of solutions of discrete delayed system $x(k+1)=ax(k)+bx(k-m)+f(k$) with commutative matrices,” Journal of Mathematical Analysis and Applications, vol. 318, no. 1, pp. 63–76, 2006, doi: 10.1016/j.jmaa.2005.05.021.
  * [7] J. Diblík and D. Y. Khusainov, “Representation of solutions of linear discrete systems with constant coefficients and pure delay,” Advances in Difference Equations, vol. 2006, pp. 1–13, 2006, doi: 10.1155/ade/2006/80825.
  * [8] J. Diblík, D. Y. Khusainov, J. Baštinec, and A. S. Sirenko, “Exponential stability of linear discrete systems with constant coefficients and single delay,” Applied Mathematics Letters, vol. 51, pp. 68–73, 2016, doi: 10.1016/j.aml.2015.07.008.
  * [9] J. Diblík, D. Y. Khusainov, J. Lukácová, and M. Ružicková, “Control of Oscillating Systems with a Single Delay.” Advances in Difference Equations, vol. 2010, pp. 1–15, 2010, doi: 10.1155/2010/108218.
  * [10] J. Diblík and B. Morávková, “Discrete matrix delayed exponential for two delays and its property,” Advances in Difference Equations, vol. 2013, no. 1, p. 139, 2013, doi: 10.1186/1687-1847-2013-139.
  * [11] J. Diblík and B. Morávková, “Representation of the solutions of linear discrete systems with constant coefficients and two delays,” Abstract and Applied Analysis, vol. 2014, pp. 1–19, 2014, doi: 10.1155/2014/320476.
  * [12] D. Y. Khusainov, J. Diblík, M. Ružicková, and J. Lukácová, “Representation of a solution of the cauchy problem for an oscillating system with pure delay,” Nonlinear Oscillations, vol. 11, no. 2, pp. 276–285, 2008, doi: 10.1007/s11072-008-0030-8.
  * [13] D. Y. Khusainov and G. V. Shuklin, “Linear autonomous time-delay system with permutation matrices solving,” Stud. Univ. Zilina. Math. Ser., vol. 17, no. 1, pp. 101–108, 2003.
  * [14] M. Li and J. Wang, “Finite time stability of fractional delay differential equations,” Applied Mathematics Letters, 2016, doi: 10.1016/j.aml.2016.09.004.
  * [15] C. Liang, J. Wang, and D. O’Regan, “Controllability of nonlinear delay oscillating systems,” J. Qual. Theory Differ. Equ., no. 47, pp. 1–18, 2017, doi: doi.org/10.14232/ejqtde.2020.1.53.
  * [16] C. Liang, J. R. Wang, and D. O’Regan, “Representation of a solution for a fractional linear system with pure delay,” Applied Mathematics Letters, vol. 77, pp. 72–78, 2017, doi: 10.1016/j.aml.2017.09.015.
  * [17] N. I. Mahmudov, “Delayed perturbation of mittag-leffler functions and their applications to fractional linear delay differential equations,” Mathematical Methods in the Applied Sciences, pp. 1–9, 2018, doi: 10.1002/mma.5446.
  * [18] N. I. Mahmudov, “Representation of solutions of discrete linear delay systems with non permutable matrices,” Applied Mathematics Letters, vol. 85, pp. 8–14, 2018, doi: 10.1016/j.aml.2018.05.015.
  * [19] N. I. Mahmudov, “A novel fractional delayed matrix cosine and sine,” Applied Mathematics Letters, 2019, doi: 10.1016/j.aml.2019.01.001.
  * [20] N. I. Mahmudov, “Delayed linear difference equations: the method of $z$-transform.” J. Qual. Theory Differ. Equ., no. 53, 2020, doi: doi.org/10.14232/ejqtde.2020.1.53.
  * [21] M. Pospíšil, “Representation of solutions of delayed difference equations with linear parts given by pairwise permutable matrices via z-transform,” Applied Mathematics and Computation, vol. 294, pp. 180–194, 2017, doi: 10.1016/j.amc.2016.09.019.