# Bispectral analysis of nonlinear mixing in a periodically driven Korteweg-de
Vries system

Ajaz Mir<EMAIL_ADDRESS>Indian Institute of Technology Jammu,
Jammu, J&K, 181221, India  Sanat Tiwari<EMAIL_ADDRESS>Indian
Institute of Technology Jammu, Jammu, J&K, 181221, India  Abhijit Sen
Institute for Plasma Research, Gandhinagar, Gujarat, 382428, India

###### Abstract

The nonlinear response of a periodically driven Korteweg-de Vries model system
is studied using a variety of nonlinear drivers and compared to previous
results obtained for a purely time-dependent sinusoidal driver [Phys. Plasmas
27, 113701 (2020)]. It is found that a nonlinear driver in the form of a
cnoidal square wave or a travelling wave driver produces a spectral response
that is closer to experimental observations [Phys. Rev. Lett. 92, 085001
(2004)] than that predicted by the simple sinusoidal driver. Using a
bispectral analysis, we also firmly establish that the nature of the nonlinear
oscillations, due to the interaction between the periodic source and the
inherent collective mode of the system, is predominantly governed by a three-
wave mixing process. Furthermore, by studying the variation in the mixing
pattern, from a broad to a sparse frequency spectrum, as a function of the
driver frequency and its functional form, we propose a means of tailoring the
nature of such patterns. Our results could find useful applications in the
experimental interpretation and manipulation of nonlinear wave mixing patterns
in weakly nonlinear and dispersive plasma systems or similar phenomena in
neutral fluids.

##  I Introduction

####

The Korteweg-de Vries (KdV) model equation has served as a standard paradigm
for the description of many low-frequency nonlinear phenomena in various
nonlinear dispersive media such as neutral fluids, plasmas, optical fibres,
vibrating lattices etc. (Washimi and Taniuti, 1966; Gasch  _et al._ , 1986;
Hao and Maris, 2001; Leblond, 2008; Peradze and Tsintsadze, 2019; Redor  _et
al._ , 2019; Congy  _et al._ , 2016; Kamchatnov and Shchesnovich, 2004). As a
fully integrable partial differential equation, the KdV equation admits a
variety of exact solutions, including solitons, cnoidal waves, dispersive
shocks etc. that have been used to describe experimental observations of
nonlinear wave phenomena in the aforementioned systems. In plasmas, the KdV
equation has successfully modelled the existence and evolution of ion and dust
acoustic shocks as well as solitary waves (Mo  _et al._ , 2013; Bandyopadhyay
_et al._ , 2008; Nakamura and Tsukabayashi, 1984). More recently, a forced KdV
(fKdV) model was used to explain the emission of precursor magneto-sonic
solitons, created by a charge bunch moving in a magnetized plasma (Kumar and
Sen, 2020).

####

The KdV model has also been investigated in the presence of an external time
and space varying perturbation in order to investigate the response of a
driven nonlinear medium. Such a forced KdV equation that had a travelling wave
source was used to determine the conditions for the resonant excitation of
nonlinear oscillatory waves in a dispersive medium by Vainberg et al.
(Vainberg  _et al._ , 1983). Using a chirped frequency driving perturbation,
Friedland et al. (Friedland  _et al._ , 2015) studied the anomalous
autoresonance threshold for the excitation of large amplitude travelling waves
in a KdV model system. Aranson et al. (Aranson  _et al._ , 1992) took a
similar approach of chirping the frequency of the driver to excite solitons in
the system. In particular, for plasmas, the fKdV model predicts existence of
precursor and pinned solitons (Arora  _et al._ , 2021; Jaiswal  _et al._ ,
2016). These nonlinear emissions can prove useful in space awareness
applications, such as in the indirect detection of charged space debris
orbiting in the earth’s ionosphere (Sen  _et al._ , 2015). More recently,  Mir
et al. Mir  _et al._ (2020) used the fKdV model to study nonlinear wave mixing
in a dusty plasma system and related their findings to past experimental
observations of Nosenko et al. (Nosenko  _et al._ , 2004). The form of the
driver, in their case, was chosen to be a simple time-varying sinusoidal
function for which an exact analytic solution of the model could be obtained.
Using such an exact solution, time series data of the system response was
obtained for parameters close to the experimental conditions and the power
spectrum of this model data was compared with the experimental power spectrum.
Nonlinear mixing (NLM) was revealed by the presence of the combination of
frequencies of various harmonics of the inherent mode of the system and the
driving wave.

####

While such identification of combination frequencies provides quick first-hand
evidence of nonlinear wave mixing, it does not conclusively establish the
origin of the combination due to a three-wave mixing process. This is because
the power spectrum does not explicitly provide the phase coupling information
about the interactions. A more precise tool for establishing the physical
origin of the wave mixing is a bispectral analysis (Rao and Gabr, 1984; Nikias
and Raghuveer, 1987; Kim and Powers, 1979) that looks at the triple-
correlation of the time series of any dynamical quantity. A finite correlation
is obtained for a frequency triad $F(f_{1})$, $F(f_{2})$ and $F(f_{1}+f_{2})$
when they are formed by a coherent phase coupling mechanism (where $F$ is the
discrete Fourier component at a given frequency). There will be no correlation
among frequencies if they are spontaneously excited modes, i.e., no coherent
phase coupling is involved. Bispectral analysis has been extensively used in
the past to investigate coherent nonlinear interactions in plasmas (van
Milligen  _et al._ , 1995; Kim  _et al._ , 1997; Nosenko  _et al._ , 2006) as
well as in many biomedical (Siu  _et al._ , 2008; Tacchino  _et al._ , 2020)
and engineering applications (Hillis  _et al._ , 2006; Hall and Giannakis,
1995).

####

In this paper, we extend the earlier spectral analysis of Mir et al. Mir  _et
al._ (2020) and carry out a bispectral analysis to firmly establish the nature
of the nonlinear mixing phenomenon. We also expand the scope of their model
calculations by going beyond the simple time-varying sinusoidal driver to
examine the influence of different functional forms of the driver on the
mixing process. These different functional forms of the driver can arise in a
variety of situations. For example, in a typical experimental situation, if
the driving wave is one of the normal modes of the system, then it is likely
to be of a nonlinear form, e.g. ion acoustic waves or dust acoustic waves
(DAWs) have been observed to grow nonlinearly into cnoidal or cnoidal square
waveforms in many experiments (Flanagan and Goree, 2010; Ruhunusiri and Goree,
2012; Chaubey  _et al._ , 2015; Teng  _et al._ , 2009). The presence and
evolution of nonlinear modes or disturbances are, in fact, generic to
nonlinear media such as cnoidal waves in metamaterials (Mo  _et al._ , 2019),
cnoidal-square waveforms of dust acoustic waves (Liu  _et al._ , 2018), and
sawtooth-like shock waves (Heinrich  _et al._ , 2009). When the driving wave
is applied externally, it can also take a variety of waveforms, such as using
a cnoidal waveform electrical or optical signals to generate frequency combs
in Kerr microresonators (Kholmyansky and Gat, 2019), use of a sawtooth wave
for higher harmonic generation in a plasma (Sharma  _et al._ , 2020) and the
use of cnoidal waves to generate nonlinear frequency combs in microring
resonators (Qi  _et al._ , 2017). The present work is aimed at studying the
impact of such waveforms on the nonlinear mixing phenomena in plasmas.

####

It is seen that changing the profile and functional dependence of the driver
can significantly alter the response pattern of the system. In particular, a
travelling wave source is seen to excite additional frequencies that are not
seen in response to purely time-varying sources.This enhanced spectrum is
shown to bear a closer resemblance to the dusty plasma experimental data of
Nosenko et al. (Nosenko  _et al._ , 2004) than the model system studied by Mir
et al. (Mir  _et al._ , 2020). In addition to changing the functional form of
the driver, we also examine the effect of changing the driver frequency
relative to the natural frequency of the non-driven system. We find a
significant difference in the resultant spectral patterns depending upon
whether the driver frequency is larger or smaller than the natural frequency.
This asymmetry in the response pattern is explained based on a higher-order
mixing process. Our findings thereby offer a means of tailoring the response
patterns of driven systems that may find practical applications.

####

The paper is organized as follows. The model fKdV equations with different
forms of purely time-dependent forcing terms $F_{s}(t)$ or with a spatio-
temporal forcing term $F_{s}(x,t)$ are described in section II. The effect of
different forcing forms on nonlinear mixing is first discussed in section III
in terms of the differences in the power spectral densities (PSD) of the time
series data of the numerical solutions of the fKdV model. This is followed up
in section IV by a bispectral analysis of the same data and its physical
interpretation. Section V is devoted to a discussion on the effect of tunable
drivers on the NLM process and a discussion on the physical origin of the
asymmetry in the driven response. A brief discussion and some concluding
remarks are given in section VI.

##  II The forced Korteweg-de Vries model

####

A generalized form of the fKdV equation for the perturbed density, $n(x,t)$
can be written down as (Mir  _et al._ , 2020):

$\frac{\partial n(x,t)}{\partial t}+\ \alpha\ n(x,t)\frac{\partial
n(x,t)}{\partial x}+\ \beta\frac{\partial^{3}n(x,t)}{\partial x^{3}}=\
F_{s}(x,t).$ (1)

####

The Eq. (1) can be derived from the full set of cold fluid equations for the
dust component in the weakly nonlinear and weakly dispersive regime Sen  _et
al._ (2015). The ions and electrons, as lighter species compared to the dust,
are assumed to obey Boltzmann relations. $\alpha$ and $\beta$ are the
characteristic parameters quantifying nonlinearity and dispersion of the
medium, respectively. $F_{s}(x,t)$ is the external spatio-temporal forcing.
For $F_{s}(x,t)=0$ one recovers the standard KdV equation that admits exact
nonlinear solutions in the form of solitons (Farina and Bulanov, 2001;
Bandyopadhyay  _et al._ , 2008; Mo  _et al._ , 2013; Peradze and Tsintsadze,
2019) and periodic cnoidal waves (Mahmood and Haas, 2014; Liu  _et al._ ,
2018).  In the context of dusty plasmas, Liu et al. (Liu  _et al._ , 2018)
have shown that a cnoidal wave solution provides an excellent fit to their
experimental observations of self-excited dust acoustic waves sustained in an
RF plasma (Flanagan and Goree, 2011).

###  II.1 Korteweg-de Vries model with a time-dependent forcing, $F_{s}(t)$

####

An analytic solution of Eq. (1) with a purely time-dependent sinusoidal
forcing that is based on Hirota’s approach (Salas, 2011) has been presented
earlier in Ref. Mir  _et al._ (2020) and is given by:

$\displaystyle n(x,t)=\ \phi(t)+\ \mu\
cn^{2}\bigg{[}\frac{\sqrt{\mu\alpha}}{2\sqrt{\beta\kappa(\kappa+2)}}\xi(x,t);\
\kappa\bigg{]}$
$\displaystyle\xi(x,t)=\bigg{(}x-\frac{\kappa+\kappa^{2}-1}{\kappa(\kappa+2)}\alpha\mu
t-\alpha\psi(t)\bigg{)}.$ (2)

In Eq. (2), $cn$ is the Jacobi elliptic function, $\mu$ is the amplitude of
the cnoidal wave, and $\kappa$ is the elliptic modulus that quantifies the
nonlinear nature of the cnoidal wave.  In the context of dusty plasma, the
elliptic modulus $\kappa$ and the total harmonic distortion of DAW (Flanagan
and Goree, 2010) were used by Liu et al. (Liu  _et al._ , 2018), to quantify
the nonlinearity of their experimentally observed waveforms. With
$\kappa\rightarrow 0$, the wave attains a cosine waveform, and as
$\kappa\rightarrow 1$, the waveform becomes nonlinear with the cosine wave
getting converted to a cnoidal waveform. At $\kappa=1$ when
$\text{cn}^{2}(x,t;\ \kappa)=\text{sech}^{2}(x,t)$ (Abramowitz and Stegun,
1965), the cnoidal wave takes the form of a single soliton. Furthermore, the
quantities $\phi$ and $\psi$ are defined as,

$\phi(t)=\ \int F_{s}(t)\ dt\hskip 17.07182pt\text{and}\hskip
17.07182pt\psi(t)=\ \int\phi(t)\ dt.$ (3)

Table 1: Functional form of different forcing profiles Forcing profiles |  $F_{s}(t)$ or $F_{s}(x,t)$
---|---
Sinusoidal $F_{\sin}(t)$ | $A_{s}\sin(2\pi f_{2}t)$
Cnoidal $F_{cn}(t)$ | $A_{s}cn[4K(\kappa_{s})f_{2}t;\kappa_{s}]$
Cnoidal-square $F_{cn^{2}}(t)$ |  $A_{s}cn^{2}[2K(\kappa_{s})f_{2}t;\kappa_{s}]$ $F_{cn^{2}}(t)=F_{cn^{2}}(t)-\left\langle F_{cn^{2}}(t)\right\rangle$
Travelling wave $F_{s}(x,t)$ |  $A_{s}\sin(k_{s}x-2\pi f_{2}t)$ Solving fKdV numerically

The $\phi(t)$ and $\psi(t)$ can be calculated either analytically or
numerically depending on the choice of the forcing profile. As a
generalization of the work reported in Ref. Mir  _et al._ (2020), where only a
sinusoidal form of $F_{s}(t)$ was used, we have considered a variety of
functional forms for the drivers in our present work. These different forcing
terms are listed in Table 1. For the case of a sinusoidal forcing profile,
$\phi(t)=-\left(A_{s}/(2\pi f_{2})\right)\cos(2\pi f_{2}t)$ and
$\psi(t)=-\left(A_{s}/(2\pi f_{2})^{2}\right)\sin(2\pi f_{2}t)$. Here $f_{2}$,
$A_{s}$ and $\kappa_{s}$ are the forcing frequency, forcing amplitude and the
forcing elliptic modulus corresponding to each forcing form $F_{s}(t)$
respectively.

Figure 1:  (a) Time series of KdV for $\mu=18.5$, $\kappa=0.7$ (such that
$f_{1}=0.7$ Hz) with $\alpha=\beta=1$. (b) Power spectrum of (a). The
nonlinear wave (non-sinusoidal) nature is evident from both the time series
(a) and the presence of harmonics in its PSD (b).

####

A truncated Fourier series expansion can approximate the non-sinusoidal
forcing profiles.  In the context of dusty plasma, Merlino et al. (Merlino
_et al._ , 2012) obtained excellent fit to their experimental DAW profiles by
retaining terms up to the second harmonic in the Fourier series expansion of
the square of the cnoidal function, namely,

$\displaystyle F_{s}(t)=\
F_{cn^{2}}(t)=A_{s}cn^{2}[2K(\kappa_{s})f_{2}t;\kappa_{s}]$ $\displaystyle
F_{cn^{2}}(t)\approx(A_{s}/2)\\{\cos(2\pi f_{2}t)$ $\displaystyle+\
0.02A_{s}\cos(2\pi(2f_{2})t)$ $\displaystyle+\
0.004A_{s}\cos(2\pi(3f_{2})t)\\}.$ (4)

Such a Fourier series representation of a non-sinusoidal forcing term allows
one to solve the fKdV system analytically by using the general solution
described by (2) and (3). This facilitates the study of the nonlinear mixing
phenomenon as a function of the various forms of the driver.

####

A nonlinear cnoidal time series solution of KdV equation and its PSD are shown
in Fig. 1(a) and Fig. 1(b) respectively. The PSD shows a fundamental frequency
at $f_{1}=0.7$ Hz along with the even and odd harmonics of the wave at
$2f_{1}$, $3f_{1}$ and so on.

Figure 2:  Power spectrum of time series obtained from fKdV model for (a) sinusoidal forcing $F_{s}(t)=F_{\sin}(t)=A_{s}\sin(2\pi f_{2}t)$ , (b) cnoidal wave forcing $F_{s}(t)=F_{cn}(t)=A_{s}cn[4K(\kappa_{s})f_{2}t;\kappa_{s}]$, (c) cnoidal-square wave forcing $F_{s}(t)=F_{cn^{2}}(t)=A_{s}cn^{2}[2K(\kappa_{s})f_{2}t;\kappa_{s}]$. Here $A_{s}=\mu$, $\kappa_{s}=0.9$ and $f_{2}=1.7$ Hz. (d) Travelling wave forcing $F_{s}(x,t)=A_{s}\sin(k_{s}x-2\pi f_{2}t)$ with no initial condition, i.e., $n(x,0)=0$. Here $A_{s}=5\mu$, $k_{s}=9k_{0}$ with $k_{0}=(2\pi)/L$, which corresponds to frequency $f_{1}=k_{s}^{3}/(2\pi)=0.7$ Hz for a system of length $L=11\pi$ and $f_{2}=1.7$ Hz. Also $\alpha=\ \beta=\ 1$ in each case.  Table 2: Dominant frequencies observed in the fKdV model for various forcing forms as shown in Fig. 2. | Fig. 2(a) | Fig. 2(b) | Fig. 2(c) | Fig. 2(d)
---|---|---|---|---
Frequency (Hz) |  $F_{\sin}(t)$ $A_{s}=\mu$ |  $F_{cn}(t)$ $A_{s}=\mu$ $\kappa_{s}=0.9$ |  $F_{cn^{2}}(t)$ $A_{s}=\mu$ $\kappa_{s}=0.9$ |  $F_{s}(x,t)$ $A_{s}=5\mu$ $k_{s}=9k_{0}$
$f_{1}=0.7$ | ✓ | ✓ | ✓ | ✓
$f_{2}=1.7$ | ✓ | ✓ | ✓ | ✓
$P_{1}=f_{2}-f_{1}$ | ✓ | ✓ | ✓ | ✓
$P_{2}=2f_{1}$ | ✓ | ✓ | ✓ | ✓
$P_{3}=3f_{1}$ | ✓ | ✓ | ✓ | ✓
$P_{4}=f_{1}+f_{2}$ | ✓ | ✓ | ✓ | ✓
$P_{5}=2f_{2}-f_{1}$ | ✓ | ✓ | ✓ | ✓
$P_{6}=4f_{1}$ | ✓ | ✓ | ✓ | ✓
$P_{7}=2f_{1}+f_{2}$ | ✓ | ✓ | ✓ | ✓
$P_{8}=2f_{2}$ |  |  | ✓ | ✓
$P_{9}=3f_{1}+f_{2}$ | ✓ | ✓ | ✓ | ✓
$P_{10}=2f_{2}+f_{1}$ | ✓ | ✓ | ✓ | ✓
$N_{1}=f_{2}-2f_{1}$ | ✓ | ✓ | ✓ | ✓
$N_{2}=3f_{1}-f_{2}$ | ✓ | ✓ | ✓ | ✓
$N_{3}=4f_{1}-f_{2}$ | ✓ | ✓ | ✓ | ✓
$N_{4}=2f_{2}-3f_{1}$ | ✓ | ✓ | ✓ |
$N_{5}=2(f_{2}-f_{1})$ | ✓ | ✓ | ✓ |
$N_{6}=3(f_{2}-f_{1})$ | ✓ | ✓ | ✓ | ✓
$N_{7}=5f_{1}$ | ✓ | ✓ | ✓ | ✓
$N_{8}=3f_{2}-2f_{1}$ | ✓ | ✓ | ✓ | ✓
$N_{9}=3f_{2}-f_{1}$ | ✓ | ✓ | ✓ | ✓
$N_{10}=4f_{1}+f_{2}$ | ✓ | ✓ | ✓ | ✓

###  II.2 Korteweg-de Vries model with spatio-temporal forcing, $F_{s}(x,t)$

####

The KdV equation with a spatio-temporal travelling waveform of the forcing is
given by:

$\displaystyle\frac{\partial n(x,t)}{\partial t}+\ \alpha\
n(x,t)\frac{\partial n(x,t)}{\partial x}+\
\beta\frac{\partial^{3}n(x,t)}{\partial x^{3}}$ $\displaystyle=\
A_{s}\sin(k_{s}x-\ 2\pi f_{2}t).$ (5)

We have solved Eq. (5) numerically using a finite difference scheme. The code
has been validated by reproducing the results of Sen et al. (Sen  _et al._ ,
2015). The frequency corresponding to the forcing wave-vector $k_{s}$ is
chosen to satisfy the linear dispersion relation obtained by setting
$\alpha=0$ in the standard KdV equation, namely $f_{1}=-\beta
k_{s}^{3}/(2\pi)$. Here, $k_{s}=nk_{0}$ with $k_{0}=2\pi/L$ being the minimum
wave-vector associated with a system of length $L$. For the numerical solution
of Eq. (5) we have taken an initial condition of $n(x,0)=0$.

Figure 3:  PSD (a) and bicoherence (b) of fKdV for sinusoidal forcing. Inset
(a) is the PSD of fKdV for $F_{s}(t)=F_{\sin}(t)=A_{s}\sin(2\pi f_{2}t)$ with
$A_{s}=18.5$ and $f_{2}=1.7$ Hz. Initial parameters are $\mu=18.5$,
$\kappa=0.7$ (such that $f_{1}=0.7$ Hz) and $\alpha=\beta=1$. The small inset
within (a) shows the form of sinusoidal forcing profile. The bicoherence map
shows patches $(\gamma^{2}\approx 1)$ indicating that the waves at $F(f_{1})$,
$F(f_{2})$ and $F(f_{1})+F(f_{2})$ are not only frequency coupled but are
phase coupled as well. This confirms coherent nonlinear interaction between
the waves at $F(f_{1})$ and $F(f_{2})$.

##  III Nonlinear mixing in the fKdV model with different forcing forms

####

To investigate the NLM process under various non-sinusoidal time-dependent
forcing forms, we have used semi-analytic solutions by using (2) and (3).
Exact analytic forms are possible for cases when $\phi$ and $\psi$ in (3) are
exactly integrable. Sinusoidal forcing is one such example. While the
solutions apply to all nonlinearity regimes and dispersion, we have taken one
parameter set throughout the paper for uniformity. We have used $\mu=18.5$,
$\kappa=0.7$ (such that $f_{1}=0.7$ Hz) and $\alpha=\beta=1$ for all cases in
the manuscript unless specified otherwise.

####

The frequencies in the power spectrum as displayed in Fig. 2 are for (a)
sinusoidal, (b) cnoidal wave, and (c) cnoidal-square wave driver in the fKdV
equation. Their values are tabulated in Table 2. All listed frequencies match
the sum and difference frequencies of $f_{1}$ (natural KdV mode), $f_{2}$
(fundamental forcing mode) and their harmonics. These additional modes are
generated via the three-wave mixing mechanism, which is further confirmed via
a bispectral analysis described in section IV.

####

The spectrum shown in Fig. 2 for each non-sinusoidal forcing form is
qualitatively similar to the sinusoidal forcing case in the range of
frequencies that carry significant power. A noticeable difference is for the
cnoidal-square forcing case where we obtain an additional frequency
$P_{8}=2f_{2}$ that also exists in the power spectrum of the fluctuations
measured in the dusty plasma experiment of Nosenko et al. (Nosenko  _et al._ ,
2004). However, this could be because the cnoidal-square forcing form
inherently contains the frequency mode at $P_{8}=2f_{2}$ in contrast to other
forcing forms, as can be seen in Table 2.

####

We have also observed NLM in fKdV model with a travelling wave forcing form
$F_{s}(x,t)=A_{s}\sin(k_{s}x-2\pi f_{2}t)$. The spectrum due to travelling
wave forcing is shown in Fig. 2(d) and is also tabulated in Table 2. The
spectrum obtained for this case is unique compared to all time-dependent
forcing cases in two ways. First, certain frequencies are missing in the
spectrum and are named as $N_{4}$ and $N_{5}$ in Table 2. The same has been
shown as red bands in Fig. 2(d). Second, we have also observed a frequency
peak $P_{8}=2f_{2}$ marked in green patch in Fig. 2(d). Both the missing and
the extra frequency are consistent with the observations in the dusty plasma
experiment (Nosenko  _et al._ , 2004). In this case, the extra frequency is
not due to the inherent property of the driver but is a genuine outcome of an
NLM process.

##  IV Bispectral analysis of nonlinear mixing

####

The bispectral analysis is a statistical tool that provides a quantitative
measure of a coherent nonlinear interaction process. If a coherent nonlinear
interaction exists between three oscillations at frequencies $F_{1}$, $F_{2}$
and $F_{1}+F_{2}$, a peak will be generated in the principal domain of the
bispectrum at the intersection between $F_{1}$ and $F_{2}$. The bispectrum of
a dynamic process is always a complex quantity, and is defined by (Rao and
Gabr, 1984; Kim and Powers, 1979; Nikias and Raghuveer, 1987; Raju  _et al._ ,
2003; Siu  _et al._ , 2008)

$B(F_{1},F_{2})=\ \langle F(f_{1})F(f_{2})F^{\star}(f_{1}+f_{2})\rangle$ (6)

where $\langle...\rangle$ is the ensemble average over multiple samples, $F$
is the Fourier transform, $\star$ is the complex conjugate and $f_{1}$,
$f_{2}$ are the two frequencies of the triad ($F_{1}$, $F_{2}$,
$F_{1}+F_{2}$). The bispectrum $B(F_{1},F_{2})$ is a function of two
frequencies and is a non-zero quantity only if a phase coupling exists between
the frequency triplet $F_{1}$, $F_{2}$ and $F_{1}+F_{2}$. $B(F_{1},F_{2})$ is
identically zero for spontaneously generated modes, i.e., the modes generated
without phase coupling. The bispectrum’s ability to retain the phase
information lost by the power spectrum (the Fourier transform of a signal)
makes it a useful tool for analyzing coherent nonlinear interactions.

Figure 4:  PSD (a) and bicoherence (b) of fKdV for cnoidal-square forcing.
Inset (a) is the PSD of fKdV for
$F_{s}(t)=F_{cn^{2}}(t)=A_{s}cn^{2}[2K(\kappa_{s})f_{2}t;\kappa_{s}]$ with
$A_{s}=18.5$, $f_{2}=1.7$ Hz and $\kappa_{s}=0.9$. Initial parameters are
$\mu=18.5$, $\kappa=0.7$ (such that $f_{1}=0.7$ Hz) and $\alpha=\beta=1$. The
small inset within (a) shows the form of cnoidal-square forcing profile. Both
the PSD (a) and bicoherence (b) show a peak at $P_{8}=2f_{2}$ (encircled),
which is inherent in the cnoidal-square wave forcing.

####

The normalized bispectrum of a time series (Kim and Powers, 1979; Nikias and
Raghuveer, 1987; Raju  _et al._ , 2003; Siu  _et al._ , 2008; Nosenko  _et
al._ , 2006) gives the bicoherence and is given by:

$\gamma^{2}(F_{1},F_{2})=\
\frac{|B(F_{1},F_{2})|^{2}}{\langle|F(f_{1})F(f_{2})|^{2}\rangle\langle|F^{\star}(f_{1}+f_{2})|^{2}\rangle}.$
(7)

Bicoherence gives a measure of phase coherence between the coupled modes. It
is a measure of the fraction of power retained by modes due to phase coupling.
Theoretically, bicoherence is 1 for phase coupled modes, i.e., modes generated
due to coherent nonlinear interaction and 0 for uncoupled modes, i.e., modes
generated spontaneously.

####

Bicoherence is symmetric about the line $F(f_{1})=F(f_{2})$ because
$\gamma^{2}(F_{1},F_{2})=\ \gamma^{2}(F_{2},F_{1})$. For our analysis, we used
$M=8192$ sampling data points and $N=16$ time series segments so that the
total length of a time series was $K=M\times N$. A statistically significant
correlation between the coherent modes is determined by the condition
$\gamma^{2}>\sqrt{6/{2N}}=0.433$ as discussed in Siu et al. (Siu  _et al._ ,
2008).

Figure 5:  (a) PSD, and (b) bicoherence of fKdV for travelling wave forcing
$F_{s}(x,t)=A_{s}\sin(k_{s}x-2\pi f_{2}t)$ with no initial perturbation, i.e.,
$n(x,0)=0$. Here $A_{s}=5\mu$, $k_{s}=9k_{0}$ with $k_{0}=(2\pi)/L$, which
corresponds to frequency $f_{1}=k_{s}^{3}/(2\pi)=0.7$ Hz for a system of
length $L=11\pi$ and $f_{2}=1.7$ Hz. Also $\alpha=\ \beta=\ 1$. The PSD and
the bicoherence show a peak at $P_{8}=2f_{2}$ (encircled) generated via a
coherent nonlinear interaction.

###  IV.1 Bispectral analysis for time-dependent sinusoidal forcing

####

In Fig. 3 we show the power spectrum (inset) and the bicoherence of the time
series obtained for a sinusoidal forcing in the fKdV model. The frequencies in
the power spectrum as shown in Fig. 3(a) [inset] are listed in column (II) of
Table 2. All the listed frequencies match the sum and difference frequencies
of $f_{1}$(fundamental KdV mode), $f_{2}$ and their harmonics. While the power
spectrum confirms the presence of new modes, it does not confirm any frequency
or phase coupling.

####

To confirm the origin of excited modes, we calculated the bicoherence of the
same time series as in Fig. 3(b). Peaks in the power spectrum map to the dark
patches in the bicoherence space shown in Fig. 3. For example, $P_{4}=2.4$ Hz
in the power spectrum generates a patch ($P_{4}^{I}$) that is at the
intersection of $F_{1}=f_{1}=0.7$ Hz and $F_{2}=f_{2}=1.7$ Hz confirming a
phase coupled and frequency-coupled mode excitation. $P_{4}=2.4$ Hz also
generates another patch in bicoherence ($P_{4}^{II}$) at the intersection of
$F_{1}=2f_{1}=1.4$ Hz and $F_{2}=f_{2}-f_{1}=1.0$ Hz. Table LABEL:Tab_3,
column (IV) lists all possible combinations of frequencies leading to patches
in the bicoherence space with $P$ and $N$ nomenclature. Peaks marked with $P$
and $N$ are the ones that were observed and missing (not observed),
respectively, in the dusty plasma experiment (Nosenko  _et al._ , 2004) which
was used for validating nonlinear mixing in the fKdV model (Mir  _et al._ ,
2020). Many patches in bicoherence space suggest multiple possible
combinations of sum and difference of frequencies leading to a single peak in
the power spectrum. Thus the present bispectral analysis confirms the three-
wave mixing results for a sinusoidal forcing of the fKdV as presented earlier
by Mir et al. (Mir  _et al._ , 2020).

####

It should be mentioned that some of the specific patches, marked with a
$\star$ in column (IV) of Table LABEL:Tab_3, do not follow the standard sum
rules but are present in the power spectrum as well as in the bicoherence
plot. We believe they represent one of the frequencies $F_{1}$ or $F_{2}$ that
appears on the coordinate axes of Fig. 3. We have independently confirmed
their source and presence by constructing a frequency and phase coupled time
series generated using all frequencies observed in the power spectrum and
switching them on and off to see their impact on the bicoherence plot.

Figure 6:  Tailoring the nonlinear mixing profiles of fKdV model with forcing
frequency $f_{2}<f_{1}$ (left panel) and $f_{2}>f_{1}$ (right panel). Left
panel: (a, c, and e) with $f_{2}=0.5$ Hz , $f_{2}=0.3$ Hz and $f_{2}=0.1$ Hz,
respectively. Right panel: (b, d, and f) with $f_{2}=0.9$ Hz , $f_{2}=1.1$ Hz
and $f_{2}=1.3$ Hz, respectively. Insets in each case display the power
spectrum of KdV equation (black bold line) and forcing form (red dash-dotted
line) i.e., $F_{s}(t)=A_{s}\sin(2\pi f_{2}t)$. Here $\alpha=\beta=1$,
$\mu=18.5$, $\kappa=0.7$ (such that $f_{1}=0.7$ Hz) and $A_{s}=\mu$ in each
case.

###  IV.2 Bispectral analysis for time-dependent non-sinusoidal forcing

####

In practical scenarios, non-sinusoidal forcing patterns are more probable. The
plasma-based nonlinear dynamical study carried out by Chaubey et al. (Chaubey
_et al._ , 2015) is one such example, where two nonlinearly excited ion
acoustic modes interact with each other. But the studies can easily be
generalized for any dispersive media.

####

This section provides the bispectral analysis for the time series obtained
from the fKdV model with cnoidal and cnoidal-square forms (Fig. 4). Columns
III-IV in Table 2 list frequencies observed in the power spectra for the two
cases as mentioned earlier. Similarly, Columns VI-VII in the Table LABEL:Tab_3
contains all possible frequency combinations due to coherent phase coupling.

####

For all non-sinusoidal forcing cases, the general appearance of the power
spectra and bicoherence diagrams appear nearly similar to the earlier
discussed sinusoidal forcing case within the frequency range of significant
amplitudes. An exception is the cnoidal-square case (Fig. 4), where we notice
an additional frequency at $P_{8}=2f_{2}$. Such a frequency has also been
observed in the referred dusty plasma experiment (Nosenko  _et al._ , 2004).
However, it should be noted that this frequency of $P_{8}=2f_{2}$ is also
inherently present in the power spectrum of the driver on account of the
nature of its profile. So, in this case, it is not possible to unambiguously
assert that the appearance of this frequency in the response spectrum is a
result of three-wave mixing.

###  IV.3 Bispectral analysis for travelling wave forcing

####

The bicoherence spectrum of the time series obtained with a travelling wave
forcing is shown in Fig. 5 with its quantitative values listed in column VII
of Table LABEL:Tab_3. It shows an absence of the $N_{4}$ and $N_{5}$
frequencies and the presence of $P_{8}=2f_{2}$ frequency in agreement with the
power spectrum analysis. We reiterate that the observation is in contrast to
the purely time-dependent sinusoidal forcing and agrees with the results
reported in the dusty plasma experiment (Nosenko  _et al._ , 2004). Also,
unlike in the case of a cnoidal-square forcing, the presence of the
$P_{8}=2f_{2}$ frequency is not an artifact of such a frequency being present
in the driver itself. In this case, it is a genuine result of a three-wave
mixing process.

Figure 7:  Tailoring the nonlinear mixing profiles by tuning the forcing form.
The approximate form of cnoidal-square forcing as in (Merlino  _et al._ ,
2012) $F_{s}(t)=\ (A_{s}/2)\\{\cos(2\pi f_{2}t)+\ A_{1}\cos(2\pi(2f_{2})t)+\
A_{2}\cos(2\pi(3f_{2})t)+\ A_{3}\cos(2\pi(f_{T})t)\\}$. (a) $A_{3}=0$, (b)
$A_{3}=0.01A_{s}$, and $f_{T}=1.25f_{2}$ Hz, (c) $A_{3}=0.01A_{s}$, and
$f_{T}=1.50f_{2}$ Hz, and (d) $A_{3}=0.01A_{s}$, and $f_{T}=1.75f_{2}$ Hz.
Here, $A_{s}=\mu$, $f_{2}=1.7$ Hz, $A_{1}=0.02A_{s}$ and $A_{2}=0.004A_{s}$
for each case.

##  V Tailoring the mixing pattern through frequency tuning of the driver

####

Based on our above-discussed analysis of the influence of the functional form
of the driver on the nonlinear mixing process, as seen in the power spectra
and bicoherence spectra, we now discuss a possible means of tailoring the
nonlinear mixing pattern through the driver parameters. We adopted two
approaches for this objective.

####

First, we varied the driver frequency $f_{2}$ towards lower or higher values
than the system’s natural frequency $f_{1}$. Second, we added specific
additional frequency signals to enhance the power in the harmonics of the
forcing profiles considered above. In the first case, we found that when the
forcing frequency $f_{2}<f_{1}$, where $f_{1}$ is the natural frequency of the
KdV model, the resultant mixing pattern has a broad spectral form. This is
seen in subplots (a), (c) and (e) of Fig. 6 for cases $f_{2}=f_{1}-0.1$,
$f_{2}=f_{1}-0.3$ and $f_{2}=f_{1}-0.5$ respectively. On the other hand, when
the forcing frequency $f_{2}>f_{1}$, the spectrum of the mixing pattern is
sparse and contains fewer frequencies. This is shown in subplots (b), (d) and
(f) of the Fig. 6 for cases $f_{2}=f_{1}+0.1$, $f_{2}=f_{1}+0.3$ and
$f_{2}=f_{1}+0.5$ respectively.

####

Such an asymmetry in the nonlinear response of a system to the variation of
the driver frequency towards or away from its natural frequency has been
observed in other systems such as microwave mixers (Pedro and Carvalho, 2003;
Lau and Yariv, 1984). The underlying physical mechanism responsible for this
behaviour is a higher-order wave mixing phenomenon called inter-modulation
distortion (IMD). Basically the output frequencies of the primary three-wave
mixing, namely, $f_{1}+f_{2}$, $f_{1}-f_{2}$, $2f_{1}$ and $2f_{2}$, further
mix with the primary modes $f_{1}$ and $f_{2}$ to give rise to $3^{\text{rd}}$
order inter-modulation products such as $2f_{1}-f_{2}$ and $2f_{2}-f_{1}$. The
$3^{\text{rd}}$ order products may further lead to IMD if they are in
proximity of the frequency of the primary modes $f_{1}$ and $f_{2}$. This
higher-order IMD condition is better fulfilled for the condition $f_{2}<f_{1}$
as compared to $f_{2}>f_{1}$. This results in a broad mixing spectrum for the
former case compared to the latter. The additional frequency content created
by the $3^{\text{rd}}$ order distortion is also known as ”Spectral Regrowth”
that has not been observed before for the forced KdV model. Controlled
regulation of such spectral regrowth through frequency tuning of the driver
could prove useful in exploring the nonlinear characteristics of weakly
dispersive and small amplitude nonlinear excitations in various plasma and
fluid media.

Table 3:  Phase coherent modes observed in bicoherence for different type of forcing forms. Column (I) | Column (II) | Column (III) | Column (IV) | Column (V) | Column (VI) | Column (VII) | Column (VIII)
---|---|---|---|---|---|---|---
$F_{1}$ (Hz) |  $F_{2}$ (Hz) |  $F_{1}+F_{2}$ (Hz) | Interpretation | Sinusoidal | Cnoidal | Cnoidal-square | Travelling wave
$f_{2}-f_{1}$ | $f_{1}$ | $f_{2}-f_{1}$ | $P_{1}^{\star}$ | ✓ | ✓ | ✓ | ✓
$f_{1}$ | $f_{1}$ | $2f_{1}$ | $P_{2}$ | ✓ | ✓ | ✓ | ✓
$2f_{1}$ | $f_{1}$ | $3f_{1}$ | $P_{3}$ | ✓ | ✓ | ✓ | ✓
$f_{2}$ | $f_{1}$ | $f_{1}+f_{2}$ | $P_{4}$ | ✓ | ✓ | ✓ | ✓
$2f_{1}$ | $f_{2}-f_{1}$ | $f_{1}+f_{2}$ | $P_{4}$ | ✓ | ✓ | ✓ | ✓
$f_{2}$ | $f_{2}-f_{1}$ | $2f_{2}-f_{1}$ | $P_{5}$ | ✓ | ✓ | ✓ | ✓
$3f_{1}$ | $f_{1}$ | $4f_{1}$ | $P_{6}$ | ✓ | ✓ | ✓ | ✓
$2f_{1}$ | $2f_{1}$ | $4f_{1}$ | $P_{6}$ | ✓ | ✓ | ✓ | ✓
$f_{1}+f_{2}$ | $f_{1}$ | $2f_{1}+f_{2}$ | $P_{7}$ | ✓ | ✓ | ✓ | ✓
$3f_{1}$ | $f_{2}-f_{1}$ | $2f_{1}+f_{2}$ | $P_{7}$ | ✓ | ✓ | ✓ | ✓
$f_{2}$ | $2f_{1}$ | $2f_{1}+f_{2}$ | $P_{7}$ | ✓ | ✓ | ✓ | ✓
$f_{2}$ | $f_{2}$ | $2f_{2}$ | $P_{8}$ |  |  | ✓ | ✓
$2f_{2}-f_{1}$ | $f_{1}$ | $2f_{2}$ | $P_{8}$ |  |  | ✓ | ✓
$f_{2}+f_{1}$ | $f_{2}-f_{1}$ | $2f_{2}$ | $P_{8}$ |  |  | ✓ | ✓
$2f_{1}+f_{2}$ | $f_{1}$ | $3f_{1}+f_{2}$ | $P_{9}$ | ✓ | ✓ | ✓ | ✓
$4f_{1}$ | $f_{2}-f_{1}$ | $3f_{1}+f_{2}$ | $P_{9}$ | ✓ | ✓ | ✓ | ✓
$f_{1}+f_{2}$ | $2f_{1}$ | $3f_{1}+f_{2}$ | $P_{9}$ | ✓ | ✓ | ✓ | ✓
$3f_{1}$ | $f_{2}$ | $3f_{1}+f_{2}$ | $P_{9}$ | ✓ | ✓ | ✓ | ✓
$2f_{2}$ | $f_{1}$ | $f_{1}+2f_{2}$ | $P_{10}$ |  |  | ✓ | ✓
$2f_{1}+f_{2}$ | $f_{2}-f_{1}$ | $f_{1}+2f_{2}$ | $P_{10}$ | ✓ | ✓ | ✓ | ✓
$2f_{2}-f_{1}$ | $2f_{1}$ | $f_{1}+2f_{2}$ | $P_{10}$ | ✓ | ✓ | ✓ | ✓
$f_{1}+f_{2}$ | $f_{2}$ | $f_{1}+2f_{2}$ | $P_{10}$ | ✓ | ✓ | ✓ | ✓
$f_{1}$ | $f_{2}-2f_{1}$ | $f_{2}-2f_{1}$ | $N_{1}^{\star}$ | ✓ | ✓ | ✓ | ✓
$2f_{1}$ | $f_{2}-2f_{1}$ | $f_{2}-2f_{1}$ | $N_{1}^{\star}$ | ✓ | ✓ | ✓ | ✓
$3f_{1}$ | $f_{2}-2f_{1}$ | $f_{2}-2f_{1}$ | $N_{1}^{\star}$ | ✓ | ✓ | ✓ | ✓
$f_{2}+f_{1}$ | $f_{2}-2f_{1}$ | $f_{2}-2f_{1}$ | $N_{1}^{\star}$ | ✓ | ✓ | ✓ | ✓
$4f_{1}$ | $f_{2}-2f_{1}$ | $f_{2}-2f_{1}$ | $N_{1}^{\star}$ | ✓ | ✓ | ✓ | ✓
$2f_{1}+f_{2}$ | $f_{2}-2f_{1}$ | $f_{2}-2f_{1}$ | $N_{1}^{\star}$ |  |  | ✓ | ✓
$3f_{1}+f_{2}$ | $f_{2}-2f_{1}$ | $f_{2}-2f_{1}$ | $N_{1}^{\star}$ | ✓ | ✓ | ✓ | ✓
$3f_{1}-f_{2}$ | $f_{2}-2f_{1}$ | $3f_{1}-f_{2}$ | $N_{2}^{\star}$ | ✓ | ✓ | ✓ | ✓
$f_{2}-f_{1}$ | $3f_{1}-f_{2}$ | $3f_{1}-f_{2}$ | $N_{2}^{\star}$ | ✓ | ✓ | ✓ | ✓
$f_{2}$ | $3f_{1}-f_{2}$ | $3f_{1}-f_{2}$ | $N_{2}^{\star}$ | ✓ | ✓ | ✓ | ✓
$f_{1}+f_{2}$ | $3f_{1}-f_{2}$ | $3f_{1}-f_{2}$ | $N_{2}^{\star}$ | ✓ | ✓ | ✓ | ✓
$2f_{2}-f_{1}$ | $3f_{1}-f_{2}$ | $3f_{1}-f_{2}$ | $N_{2}^{\star}$ | ✓ | ✓ | ✓ | ✓
$2f_{2}$ | $3f_{1}-f_{2}$ | $3f_{1}-f_{2}$ | $N_{2}^{\star}$ |  |  | ✓ | ✓
$f_{1}$ | $3f_{1}-f_{2}$ | $4f_{1}-f_{2}$ | $N_{3}$ | ✓ | ✓ | ✓ | ✓
$4f_{1}-f_{2}$ | $f_{2}-2f_{1}$ | $4f_{1}-f_{2}$ | $N_{3}^{\star}$ | ✓ | ✓ | ✓ | ✓
$4f_{1}-f_{2}$ | $f_{2}-f_{1}$ | $4f_{1}-f_{2}$ | $N_{3}^{\star}$ | ✓ | ✓ | ✓ | ✓
$f_{2}$ | $4f_{1}-f_{2}$ | $4f_{1}-f_{2}$ | $N_{3}^{\star}$ | ✓ | ✓ | ✓ | ✓
$2f_{2}-f_{1}$ | $4f_{1}-f_{2}$ | $4f_{1}-f_{2}$ | $N_{3}^{\star}$ | ✓ | ✓ | ✓ | ✓
$f_{2}-f_{1}$ | $f_{2}-2f_{1}$ | $2f_{2}-3f_{1}$ | $N_{4}$ | ✓ | ✓ | ✓ |
$2f_{2}-3f_{1}$ | $3f_{1}-f_{2}$ | $2f_{2}-3f_{1}$ | $N_{4}^{\star}$ | ✓ | ✓ | ✓ |
$2f_{2}-3f_{1}$ | $4f_{1}-f_{2}$ | $2f_{2}-3f_{1}$ | $N_{4}^{\star}$ | ✓ | ✓ | ✓ |
$2f_{1}$ | $2f_{2}-3f_{1}$ | $2f_{2}-3f_{1}$ | $N_{4}^{\star}$ | ✓ | ✓ | ✓ |
$3f_{1}$ | $2f_{2}-3f_{1}$ | $2f_{2}-3f_{1}$ | $N_{4}^{\star}$ |  |  | ✓ |
$4f_{1}$ | $2f_{2}-3f_{1}$ | $2f_{2}-3f_{1}$ | $N_{4}^{\star}$ | ✓ | ✓ | ✓ |
$f_{2}$ | $f_{2}-2f_{1}$ | $2(f_{2}-f_{1})$ | $N_{5}$ | ✓ | ✓ | ✓ |
$2(f_{2}-f_{1})$ | $3f_{1}-f_{2}$ | $2(f_{2}-f_{1})$ | $N_{5}^{\star}$ | ✓ | ✓ | ✓ |
$2f_{2}-3f_{1}$ | $f_{1}$ | $2(f_{2}-f_{1})$ | $N_{5}$ | ✓ | ✓ | ✓ |
$2(f_{2}-f_{1})$ | $f_{1}$ | $2(f_{2}-f_{1})$ | $N_{5}^{\star}$ | ✓ | ✓ | ✓ |
$f_{2}-f_{1}$ | $f_{2}-f_{1}$ | $2(f_{2}-f_{1})$ | $N_{5}$ | ✓ | ✓ | ✓ |
$2(f_{2}-f_{1})$ | $4f_{1}-f_{2}$ | $2(f_{2}-f_{1})$ | $N_{5}^{\star}$ | ✓ | ✓ | ✓ |
$2(f_{2}-f_{1})$ | $2f_{1}$ | $2(f_{2}-f_{1})$ | $N_{5}^{\star}$ |  |  | ✓ |
$3f_{1}$ | $2(f_{2}-f_{1})$ | $2(f_{2}-f_{1})$ | $N_{5}^{\star}$ | ✓ | ✓ | ✓ |
$2f_{2}-f_{1}$ | $f_{2}-2f_{1}$ | $3(f_{2}-f_{1})$ | $N_{6}$ | ✓ | ✓ | ✓ | ✓
$3(f_{2}-f_{1})$ | $3f_{1}-f_{2}$ | $3(f_{2}-f_{1})$ | $N_{6}^{\star}$ |  |  | ✓ | ✓
$2(f_{2}-f_{1})$ | $f_{2}-f_{1}$ | $3(f_{2}-f_{1})$ | $N_{6}$ | ✓ | ✓ | ✓ | ✓
$3(f_{2}-f_{1})$ | $4f_{1}-f_{2}$ | $3(f_{2}-f_{1})$ | $N_{6}^{\star}$ | ✓ | ✓ | ✓ | ✓
$f_{2}$ | $2f_{2}-3f_{1}$ | $3(f_{2}-f_{1})$ | $N_{6}$ | ✓ | ✓ | ✓ | ✓
$2f_{1}+f_{2}$ | $3f_{1}-f_{2}$ | $5f_{1}$ | $N_{7}$ | ✓ | ✓ | ✓ | ✓
$5f_{1}$ | $f_{2}-2f_{1}$ | $5f_{1}$ | $N_{7}^{\star}$ | ✓ | ✓ | ✓ | ✓
$4f_{1}$ | $f_{1}$ | $5f_{1}$ | $N_{7}$ | ✓ | ✓ | ✓ | ✓
$f_{1}+f_{2}$ | $4f_{1}-f_{2}$ | $5f_{1}$ | $N_{7}$ | ✓ | ✓ | ✓ | ✓
$3f_{1}$ | $2f_{1}$ | $5f_{1}$ | $N_{7}$ | ✓ | ✓ | ✓ | ✓
$2f_{2}$ | $f_{2}-2f_{1}$ | $3f_{2}-2f_{1}$ | $N_{8}$ |  |  | ✓ | ✓
$3f_{2}-2f_{1}$ | $3f_{1}-f_{2}$ | $3f_{2}-2f_{1}$ | $N_{8}^{\star}$ | ✓ | ✓ | ✓ | ✓
$3(f_{2}-f_{1})$ | $f_{1}$ | $3f_{2}-2f_{1}$ | $N_{8}$ | ✓ | ✓ | ✓ | ✓
$f_{1}+f_{2}$ | $2f_{2}-3f_{1}$ | $3f_{2}-2f_{1}$ | $N_{8}$ | ✓ | ✓ | ✓ | ✓
$2(f_{2}-f_{1})$ | $2f_{2}-3f_{1}$ | $3f_{2}-2f_{1}$ | $N_{8}$ | ✓ | ✓ | ✓ | ✓
$2(f_{2}-f_{1})$ | $f_{2}$ | $3f_{2}-2f_{1}$ | $N_{8}$ | ✓ | ✓ | ✓ | ✓
$2f_{2}+f_{1}$ | $f_{2}-2f_{1}$ | $3f_{2}-f_{1}$ | $N_{9}$ | ✓ | ✓ | ✓ | ✓
$3f_{2}-2f_{1}$ | $f_{1}$ | $3f_{2}-f_{1}$ | $N_{9}$ | ✓ | ✓ | ✓ | ✓
$2f_{2}$ | $f_{2}-f_{1}$ | $3f_{2}-f_{1}$ | $N_{9}$ |  |  | ✓ | ✓
$2f_{1}+f_{2}$ | $2f_{2}-3f_{1}$ | $3f_{2}-f_{1}$ | $N_{9}$ | ✓ | ✓ | ✓ | ✓
$3(f_{2}-f_{1})$ | $2f_{1}$ | $3f_{2}-f_{1}$ | $N_{9}$ | ✓ | ✓ | ✓ | ✓
$2f_{2}-f_{1}$ | $f_{2}$ | $3f_{2}-f_{1}$ | $N_{9}$ | ✓ | ✓ | ✓ | ✓
$f_{1}+f_{2}$ | $2(f_{2}-2f_{1})$ | $3f_{2}-f_{1}$ | $N_{9}$ | ✓ | ✓ | ✓ | ✓
$2f_{2}+f_{1}$ | $3f_{1}-f_{2}$ | $4f_{1}+f_{2}$ | $N_{10}$ | ✓ | ✓ | ✓ | ✓
$3f_{1}+f_{2}$ | $f_{1}$ | $4f_{1}+f_{2}$ | $N_{10}$ | ✓ | ✓ | ✓ | ✓
$5f_{1}$ | $f_{2}-f_{1}$ | $4f_{1}+f_{2}$ | $N_{10}$ | ✓ | ✓ | ✓ | ✓
$2f_{2}$ | $4f_{1}-f_{2}$ | $4f_{1}+f_{2}$ | $N_{10}$ |  |  | ✓ | ✓
$2f_{1}+f_{2}$ | $2f_{1}$ | $4f_{1}+f_{2}$ | $N_{10}$ | ✓ | ✓ | ✓ | ✓
$f_{1}+f_{2}$ | $3f_{1}$ | $4f_{1}+f_{2}$ | $N_{10}$ | ✓ | ✓ | ✓ | ✓

####

We next discuss our second approach towards influencing the nature of the
system response spectra by tuning the driver characteristics. As mentioned
earlier, we had found that the effects of typical non-sinusoidal (e.g., a
cnoidal wave and a cnoidal-square wave) drivers on mixing were very similar to
that of a sinusoidal driver (refer to Fig. 2), within a dominant frequency
range. The standard driver limits both the amplitude and frequency of their
harmonics, restricting the richness of nonlinear mixing. The non-sinusoidal
drivers show dominant interactions because of their fundamental mode only, as
the amplitudes of their harmonics drop significantly. Hence, in our second
approach, we considered tuning the driver by externally pumping up the power
in a frequency close to one of the harmonics of the driver. This is done by
adding an extra sinusoidal component at a specific frequency and amplitude to
the Fourier series representation of a non-sinusoidal forcing term.

####

Fig. 7(a) shows the time series spectrum obtained from fKdV with a cnoidal-
square forcing. The spectrum is similar to the pure sinusoidal forcing case
with an additional frequency at $P_{8}=2f_{2}$ due to the driver’s inherent
profile. The power at the first harmonic ($2f_{2}$) is significantly reduced
and lies far away from the natural mode ($f_{1}$). Hence, the fundamental
forcing mode ($f_{2}$) shows more interaction with the natural mode $(f_{1})$
and inhibits the interaction of the first harmonic with the natural mode.
Figs. 7(b), 7(c) and 7(d) show the spectrum of mixing due to an extra
tailoring mode ($f_{T}$) at $f_{T}=1.25f_{2}$, $f_{T}=1.50f_{2}$ and
$f_{T}=1.75f_{2}$ respectively, each with amplitude of $50\%A_{f_{2}}$, where
$A_{f_{2}}$ is the amplitude of the fundamental mode. This results in a
frequency bunching due to the fractional frequency mode (i.e., the mode in
between the fundamental mode and the first harmonic) close to the fundamental
mode leading to an altered frequency spectrum. Thus both these techniques
offer one a simple but effective means of tailoring the response spectrum of
the driven KdV system.

##  VI Discussion and Conclusions

####

To summarize, we have carried out a detailed numerical study of the driven
response of a model KdV equation to a variety of driving sources that range
from a simple time-varying sine wave to a spatio-temporally varying plane wave
to nonlinear cnoidal waveforms. One of the first objectives of this study is
to firmly establish the nature of the process underlying the wave mixing
taking place in the system. An earlier study (Mir  _et al._ , 2020), using a
simple time-varying sinusoidal driver, had concluded that three-wave couplings
were responsible for the wave mixing in the system. The conclusion was based
on identifying some of the combination frequencies in the power spectrum of
the nonlinear fluctuations. However, the power spectrum information alone is
not sufficient to establish the existence of a three-wave coupling event. This
is because the power spectrum does not have the phase coupling information
about the interactions. A more precise tool for establishing the existence of
three-wave coupling is a bispectral analysis that looks at the triple-
correlation of the time series of any dynamical quantity. A finite correlation
is obtained for a frequency triad when they are formed by a coherent phase
coupling mechanism. In this paper, we have subjected the earlier data of Ref.
(Mir  _et al._ , 2020) to bispectral analysis and have confirmed that the
mixing process arising from a simple time-varying sinusoidal driver is indeed
due to three-wave interactions.

####

We have next gone on to generalize the findings of Ref. (Mir  _et al._ , 2020)
by changing the nature of the driving source and studying the impact of such
changes on the mixing process. In particular, we have chosen time-varying
nonlinear drivers in the form of cnoidal waves and the square of cnoidal waves
and compared the resultant power spectra and bicoherence spectra with those of
the purely sinusoidal driver. The absence of certain spectral lines or the
presence of new ones in the response spectra have been identified, and their
origin is traced to the changing nature of the natural spectra of the
nonlinear driving terms. We have also considered a linear driving term with
both a temporal and spatial variation and constitutes a travelling waveform.
In this case, we find an additional response frequency that was not present in
the case of the sinusoidal driver but had been experimentally observed in Ref.
(Nosenko  _et al._ , 2004). This is not surprising as the sine wave driver, in
the earlier model calculation, was adopted as an approximation to the long
wavelength ($k\rightarrow 0$) dust acoustic wave that had been excited in the
experimental system by an external laser. Such external perturbations can also
give rise to a nonlinear excitation e.g., as a cnoidal wave that can then act
as a driver of the KdV system. Thus both the nonlinear waveforms and the
travelling waveform considered as drivers in our model calculation can find
useful applications in experimental scenarios.

####

Looking at the sensitivity of the response spectra to the nature of the
driver, we have further extended our explorations to alter not just the form
of the driver but also the frequency of the driver with respect to the natural
frequency of the KdV system. Our results point to a novel means of altering
the spectral density of the response spectra by manipulating the driver
frequency to be smaller or larger than the system’s natural frequency or by
artificially injecting power in the driver at a frequency that is somewhat
removed from its fundamental frequency. The physical origin of these spectral
changes lies in higher-order wave interactions that can create additional
frequencies and alter the spectrum’s nature.

####

Finally, we would like to remark that our results have broad applicability and
relevance for understanding nonlinear phenomena in plasmas. Nonlinear mixing
is at the heart of wave-wave interactions in plasmas (Wilhelmsson, 1976) and
is responsible for such processes as harmonic generation parametric
instabilities (Shukla, 2000), the onset of weak turbulence, etc.. Such
interactions have been widely studied in the context of laser heating of
plasmas (Brueckner and Jorna, 1974), radiofrequency heating in magnetic
confinement devices (Tsytovich and Stenflo, 1973; Ritz  _et al._ , 1989), and
understanding a variety of electrostatic and electromagnetic fluctuations in
space plasmas (Chian, 1996). One of the principal objectives in such studies
is identifying the precise nature of the wave-wave interaction responsible for
a particular phenomenon. This becomes a challenging task because many of the
nonlinear phenomena often take place simultaneously. Bispectral analysis is a
convenient and accurate tool for carrying out such an analysis. The KdV based
model used in our analysis is a convenient semi-analytic framework for
demonstrating the power and utility of this diagnostic. The conclusions are by
no means limited to this model. They are also applicable to any physical
scenario where an external (or internally generated) wave drives the plasma
and induces nonlinear wave mixing. Our present results can be the basis for
understanding nonlinear mixing in plasma systems that are weakly nonlinear and
weakly dispersive. Such plasma systems are widely found in space and
laboratory setups and have served as convenient media for the study of
solitons and other coherent structures. We hope that our work will stimulate
further theoretical work in extending the present studies to plasma models of
higher dimensions and stronger nonlinearities.

###### Acknowledgements.

Work done by AM and ST was supported by IIT Jammu Seed Grant No. SG0012. AM
and ST acknowledge the use of AGASTYA HPC for present studies. ST also
acknowledges SERB Grant No. CRG/2020/003653 for partial support for the work.
AS is thankful to the Indian National Science Academy (INSA) for their support
under the INSA Senior Scientist Fellowship scheme.

## References

  * Washimi and Taniuti (1966) H. Washimi and T. Taniuti, Phys. Rev. Lett. 17, 996 (1966).
  * Gasch  _et al._ (1986) A. Gasch, T. Berning, and D. Jäger, Phys. Rev. A 34, 4528 (1986).
  * Hao and Maris (2001) H.-Y. Hao and H. J. Maris, Phys. Rev. B 64, 064302 (2001).
  * Leblond (2008) H. Leblond, Phys. Rev. A 78, 013807 (2008).
  * Peradze and Tsintsadze (2019) G. M. Peradze and N. L. Tsintsadze, Low Temp. Phys. 45, 103 (2019).
  * Redor  _et al._ (2019) I. Redor, E. Barthélemy, H. Michallet, M. Onorato, and N. Mordant, Phys. Rev. Lett. 122, 214502 (2019).
  * Congy  _et al._ (2016) T. Congy, A. M. Kamchatnov, and N. Pavloff, Phys. Rev. A 93, 043613 (2016).
  * Kamchatnov and Shchesnovich (2004) A. M. Kamchatnov and V. S. Shchesnovich, Phys. Rev. A 70, 023604 (2004).
  * Mo  _et al._ (2013) Y. C. Mo, R. A. Kishek, D. Feldman, I. Haber, B. Beaudoin, P. G. O’Shea, and J. C. T. Thangaraj, Phys. Rev. Lett. 110, 084802 (2013).
  * Bandyopadhyay  _et al._ (2008) P. Bandyopadhyay, G. Prasad, A. Sen, and P. K. Kaw, Phys. Rev. Lett. 101, 065006 (2008).
  * Nakamura and Tsukabayashi (1984) Y. Nakamura and I. Tsukabayashi, Phys. Rev. Lett. 52, 2356 (1984).
  * Kumar and Sen (2020) A. Kumar and A. Sen, New J. Phys. 22, 073057 (2020).
  * Vainberg  _et al._ (1983) Y. R. Vainberg, B. I. Meerson, and P. V. Sasorov, Radiophys. Quantum Electron. 26, 1114 (1983).
  * Friedland  _et al._ (2015) L. Friedland, A. G. Shagalov, and S. V. Batalov, Phys. Rev. E 92, 042924 (2015).
  * Aranson  _et al._ (1992) I. Aranson, B. Meerson, and T. Tajima, Phys. Rev. A 45, 7500 (1992).
  * Arora  _et al._ (2021) G. Arora, P. Bandyopadhyay, M. G. Hariprasad, and A. Sen, Phys. Rev. E 103, 013201 (2021).
  * Jaiswal  _et al._ (2016) S. Jaiswal, P. Bandyopadhyay, and A. Sen, Phys. Rev. E 93, 041201 (2016).
  * Sen  _et al._ (2015) A. Sen, S. Tiwari, S. Mishra, and P. Kaw, Adv. Space Res. 56, 429 (2015).
  * Mir  _et al._ (2020) A. A. Mir, S. K. Tiwari, J. Goree, A. Sen, C. Crabtree, and G. Ganguli, Phys. Plasmas 27, 113701 (2020).
  * Nosenko  _et al._ (2004) V. Nosenko, K. Avinash, J. Goree, and B. Liu, Phys. Rev. Lett. 92, 085001 (2004).
  * Rao and Gabr (1984) T. S. Rao and M. M. Gabr,  _Lecture Notes in Statistics: An Introduction to Bispectral Analysis and Bilinear Time Series Models_ (Springer-Verlag, New York, 1984).
  * Nikias and Raghuveer (1987) C. L. Nikias and M. R. Raghuveer, Proc. IEEE 75, 869 (1987).
  * Kim and Powers (1979) Y. C. Kim and E. J. Powers, IEEE Trans. Plasma Sci. 7, 120 (1979).
  * van Milligen  _et al._ (1995) B. P. van Milligen, C. Hidalgo, and E. Sánchez, Phys. Rev. Lett. 74, 395 (1995).
  * Kim  _et al._ (1997) J. S. Kim, R. J. Fonck, R. D. Durst, E. Fernandez, P. W. Terry, S. F. Paul, and M. C. Zarnstorff, Phys. Rev. Lett. 79, 841 (1997).
  * Nosenko  _et al._ (2006) V. Nosenko, J. Goree, and F. Skiff, Phys. Rev. E 73, 016401 (2006).
  * Siu  _et al._ (2008) K. L. Siu, J. M. Ahn, K. Ju, M. Lee, K. Shin, and K. H. Chon, IEEE Trans. Biomedical Eng. 55, 1512 (2008).
  * Tacchino  _et al._ (2020) G. Tacchino, S. Coelli, P. Reali, M. Galli, and A. M. Bianchi, IEEE Trans. Biomedical Eng. 67, 2696 (2020).
  * Hillis  _et al._ (2006) A. J. Hillis, S. A. Neild, B. W. Drinkwater, and P. D. Wilcox, Proc. Math. Phys. Eng. Sci. 462, 1515 (2006).
  * Hall and Giannakis (1995) T. E. Hall and G. B. Giannakis, IEEE Trans. Image Processing 4, 996 (1995).
  * Flanagan and Goree (2010) T. M. Flanagan and J. Goree, Phys. Plasmas 17, 123702 (2010).
  * Ruhunusiri and Goree (2012) W. D. S. Ruhunusiri and J. Goree, Phys. Rev. E 85, 046401 (2012).
  * Chaubey  _et al._ (2015) N. Chaubey, S. Mukherjee, A. N. Sekar Iyengar, and A. Sen, Phys. Plasmas 22, 022312 (2015).
  * Teng  _et al._ (2009) L.-W. Teng, M.-C. Chang, Y.-P. Tseng, and L. I, Phys. Rev. Lett. 103, 245005 (2009).
  * Mo  _et al._ (2019) C. Mo, J. Singh, J. R. Raney, and P. K. Purohit, Phys. Rev. E 100, 013001 (2019).
  * Liu  _et al._ (2018) B. Liu, J. Goree, T. M. Flanagan, A. Sen, S. K. Tiwari, G. Ganguli, and C. Crabtree, Phys. Plasmas 25, 113701 (2018).
  * Heinrich  _et al._ (2009) J. Heinrich, S.-H. Kim, and R. L. Merlino, Phys. Rev. Lett. 103, 115002 (2009).
  * Kholmyansky and Gat (2019) D. Kholmyansky and O. Gat, Phys. Rev. A 100, 063809 (2019).
  * Sharma  _et al._ (2020) S. Sharma, N. Sirse, and M. M. Turner, Plasma Sources Sci. Technol. 29, 114001 (2020).
  * Qi  _et al._ (2017) Z. Qi, G. D’Aguanno, and C. R. Menyuk, J. Opt. Soc. Am. B 34, 785 (2017).
  * Farina and Bulanov (2001) D. Farina and S. V. Bulanov, Phys. Rev. Lett. 86, 5289 (2001).
  * Mahmood and Haas (2014) S. Mahmood and F. Haas, Phys. Plasmas 21, 102308 (2014).
  * Flanagan and Goree (2011) T. M. Flanagan and J. Goree, Phys. Plasmas 18, 013705 (2011).
  * Salas (2011) A. H. Salas, Nonlinear Anal. Real World Appl. 12, 1314 (2011).
  * Abramowitz and Stegun (1965) M. Abramowitz and I. A. Stegun,  _Handbook of Mathematical Functions : With Formulas, Graphs, and Mathematical Tables_ (Dover Publications, New York, 1965).
  * Merlino  _et al._ (2012) R. L. Merlino, J. R. Heinrich, S.-H. Hyun, and J. K. Meyer, Phys. Plasmas 19, 057301 (2012).
  * Raju  _et al._ (2003) D. Raju, O. Sauter, and J. B. Lister, Plasma Phys. Control. Fusion 45, 369 (2003).
  * Pedro and Carvalho (2003) J. C. Pedro and N. B. Carvalho,  _Intermodulation Distortion in Microwave and Wireless Circuits_ (Artech House, Boston, 2003).
  * Lau and Yariv (1984) K. Y. Lau and A. Yariv, Appl. Phys. Lett. 45, 1034 (1984).
  * Wilhelmsson (1976) H. Wilhelmsson, Physica B+C 82, 52 (1976).
  * Shukla (2000) P. K. Shukla, Phys. Rev. Lett. 84, 5328 (2000).
  * Brueckner and Jorna (1974) K. A. Brueckner and S. Jorna, Rev. Mod. Phys. 46, 325 (1974).
  * Tsytovich and Stenflo (1973) V. Tsytovich and L. Stenflo, Phys. Lett. A 43, 7 (1973).
  * Ritz  _et al._ (1989) C. P. Ritz, E. J. Powers, and R. D. Bengtson, Phys. Fluids B: Plasma Phys. 1, 153 (1989).
  * Chian (1996) A.-L. Chian, Astrophys. Space Sci. 242, 249 (1996).


# SCHeMa: Scheduling Scientific Containers on a Cluster of Heterogeneous
Machines

Thanasis Vergoulis<EMAIL_ADDRESS>0000-0003-0555-4128 IMSI, “Athena”
RC ,  Konstantinos Zagganas<EMAIL_ADDRESS>Univ. of the Peloponnese & IMSI,
“Athena” RC ,  Loukas Kavouras<EMAIL_ADDRESS>IMSI, “Athena” RC ,
Martin Reczko<EMAIL_ADDRESS>B.S.R.C. “Al. Fleming” ,  Stelios Sartzetakis
<EMAIL_ADDRESS>IMSI, “Athena” RC & GRNET and  Theodore Dalamagas
<EMAIL_ADDRESS>IMSI, “Athena” RC

###### Abstract.

In the era of data-driven science, conducting computational experiments that
involve analysing large datasets using heterogeneous computational clusters,
is part of the everyday routine for many scientists. Moreover, to ensure the
credibility of their results, it is very important for these analyses to be
easily reproducible by other researchers. Although various technologies, that
could facilitate the work of scientists in this direction, have been
introduced in the recent years, there is still a lack of open source platforms
that combine them to this end. In this work, we describe and demonstrate
SCHeMa, an open-source platform that facilitates the execution and
reproducibility of computational analysis on heterogeneous clusters,
leveraging containerization, experiment packaging, workflow management, and
machine learning technologies.

containerization, data-driven science, reproducibility

††copyright: none††ccs: Information systems Computing platforms

##  1\. Introduction

The computational analysis of large datasets has been established as an
important part of the daily routine of scientists in many disciplines, shaping
the field of _data-driven science_. Due to the large size of the datasets,
such computations are assigned to the nodes of computational clusters owned by
the academic or research institution to which the scientists belong. It is
very common that such computational clusters are heterogeneous, consisting of
machines of very diverse specifications (CPUs, memory, disk, etc.) or
capabilities (e.g., support for FPGAs and other accelerators). This
_heterogeneity_ is due to the fact that these infrastructures (a) have to
serve a variety of analysis tasks, each having its own special needs (e.g., to
exploit accelerators), and (b) are usually built incrementally, with equipment
units being procured at different (and maybe significantly distant) time
periods, based on the availability of funds.

Similarly to any other scientific experiment, replicating and reproducing the
results of a computational analysis is an important guarantee for its
credibility. This is especially important nowadays, since there is an
increasing concern in the research and academic community about the existence
of a large number of scientific works that cannot be reproduced (Baker, 2016).
Although it may be an exaggeration that we experience an ongoing
_reproducibility crisis_ , it is unarguable that this is an important
phenomenon that needs to be addressed (Fanelli, 2018). This is why
facilitating reproducibility has become an important topic of many research
and academic disciplines.

In the context of data-driven science, facilitating reproducibility can be
translated into making the datasets, the code, and the configurations used for
the analysis openly available. Motivated by this need, approaches to pack up
scientific computational experiments (e.g., RO-crate (Carragáin et al., 2019))
have been introduced. Although such packages are really useful, their true
potential is not easily unleashed due to the fact that computer environments
(e.g., software libraries, packages, programming languages) are complex and
rapidly evolving, making the reproducibility and extension of computational
analyses challenging and tedious (Boettiger, 2014). For instance, although the
code of a computational experiment may be openly available (e.g., on GitHub or
other similar repositories), in many cases its installation may require
searching for old or, even, deprecated versions of third-party software,
resolving conflicts between different versions of particular dependencies
required by different software units, adapting the code to work on a different
operating system, and so on. Given the fact that, usually, scientific software
lacks comprehensive documentation, tasks like the previous may require
significant technical skills that most scientists do not possess.

The missing piece in this puzzle is the use of _containerization technologies_
(e.g., Docker, Singularity), along with experiment packaging, which has the
potential to alleviate issues like these (Boettiger, 2014; Chamberlain and
Schommer, 2014; Jimenez et al., 2015). Such technologies allow the code of a
complex software unit to be packed up with its dependencies so it can be
easily and reliably executed in a variety of computing environments (laptops,
PCs, Cloud or HPC nodes, etc). This packaged software units are known as
_container images_ , and can further facilitate the reproducibility of
computational analysis tasks since they are already configured and ready-to-
use without requiring advanced technical skills for their installation. In
addition, in the last years, a large number of scientific containers have been
released in public repositories (e.g., in the time of writing, Biocontainers
(da Veiga Leprevost et al., 2017) currently contain more than $2076$
containers).

It is evident, from the previous discussion, that various technologies, that
could facilitate the work of scientists in the direction of data-driven
science, have been introduced in the recent years. Consequently, various
platforms that attempt to combine these technologies to provide useful
services for the research community exist. For example, two relevant platforms
are EOSC Life’s WorkflowHub111WorkflowHub: https://workflowhub.eu/ and Galaxy
Europe222Galaxy Europe: https://usegalaxy.eu/ (Afgan et al., 2018). The former
is an under-development, federated repository of workflows that is based on
the SEEK platform (Wolstencroft et al., 2015), however it does not support
workflow execution by itself. The latter is a feature-rich platform that
supports job scheduling on heterogeneous clusters.

In this work, we introduce SCHeMa (Scheduler for scientific Containers on
clusters of Heterogeneous Machines) a new open source333SCHeMa’s code
repository: https://github.com/athenarc/schema (GNU/GPL license) platform that
aims to increase diversity in the aforementioned ecosystem, focusing on the
implementation of GA4GH specifications, the native support for CWL workflows,
the facilitation of the reproducibility of computational experiments, and the
adaptive scheduling of jobs on heterogeneous clusters. The platform exploits
containerization, experiment packaging, and workflow management technologies
to ease reproducibility, while it leverages machine learning technologies to
automatically identify the type of node that is more suitable to undertake
each submitted computational task.

It is worth mentioning that a deployment of SCHeMa powers the on-demand
computations performed on the Cloud infrastructure of the ELIXIR-GR
community444ELIXIR-GR Cloud Infrastructure (in beta): https://egci-
beta.imsi.athenarc.gr/, consisting of 45 physical nodes with 2600 CPU cores,
24 TBs RAM memory and 1 PB of storage, in total.

##  2\. System Overview

###  2.1. Design Objectives

SCHeMa, our open-source platform, has been designed and implemented with the
aim to assist the work of scientists in the era of data-driven and
reproducible science. In this context, our design had two main objectives: (a)
to make the reproduction of any computational experiment performed in the
platform as easy as possible, and (b) to allocate the resources of the
underlying heterogeneous cluster as wisely as possible.

Regarding the first objective, as was mentioned in Section 1, a set of
technologies should be combined together to achieve the desired goal. In
addition, we had to identify the most appropriate technologies to be used from
a multitude from available options. Our selection was made taking into
consideration the maturity of the technologies to be used, their compatibility
to each other, and the level of their dissemination in the scientific
community. Based on these criteria, we selected to adopt CWL555CWL website:
https://www.commonwl.org/ to describe software interfaces and workflows. We
selected RO-crate (Carragáin et al., 2019) to create packages that represent
computational experiments by storing the CWL description and the particular
configuration of the software used along with the input and output datasets
involved. Finally, we used containerization (Docker in particular) as a
technology to enable the easy software execution, without requiring technical
knowledge about building the involved software packages.

Regarding the second objective, we approached the problem of selecting the
most appropriate type of machine in the cluster as a classification problem
where the input features are metadata relevant to the job to be executed
(i.e., characteristics of the inputs used), while each class represents a
particular type of machine (e.g., regular-memory machine, large-memory
machine, slow-disk machine, etc). In particular, we implemented a profiler
that, after a user request, is able to analyse the execution behavior of a
software of interest by monitoring its execution on a wide range of different
inputs and configurations. The profiler trains and evaluates the accuracy of
various classification approaches in different hyperparameter configurations
and selects the best performing one as the prediction model to be used for any
execution of the particular software in the future.

###  2.2. Architecture

SCHeMa implements a wide range of functionalities to assist scientists in the
data-driven and reproducible science era. Most notable are (a) the option to
upload custom-made scientific containers or container-based workflows, (b) a
wizard and an API that facilitate the execution of individual containers or
workflows, (c) a machine-learning-based classifier that (after a required
training phase) identifies the type of cluster node which is more appropriate
to undertake a particular computational task, (d) a monitor that informs the
users about the consumption of computational resources, (e) a wizard to
transform executed analyses into RO-crate-based “experiment packages”, and (f)
a wizard to facilitate interconnection with open data repository services.
Figure 1 summarises SCHeMa’s architecture, which supports these (and some
extra) functionalities. In the following sections we discuss SCHeMa’s external
dependencies, as well as the implementation details of its internal software
components.

####  2.2.1. External dependencies

SCHeMa’s function depends on the existence of a couple of external
installations, the most important of which are the following:

  * •

A _Kubernetes 666Kubernetes website: https://kubernetes.io/ installation_
should be deployed on the computational cluster to be used. Kubernetes
undertakes the low-level orchestration and monitoring of the computational
jobs and interacts with the _Job Classifier_ component that provides the feed
of requested jobs along with recommendations about the most suitable types of
node for each job.

  * •

A _distributed file system_ should be installed on the hard disks of the
machines of the cluster. This file system is used to store input/output data
required/produced by the computational tasks. Currently the implementation
supports NFS volumes.

  * •

A private _Docker image registry_ supporting TLS security and user
authentication is required. SCHeMa uses this registry to upload container
images. This ensures that user-uploaded images remain isolated from the
outside world (especially those meant to be private).

Figure 1. The architecture of SCHeMa.

####  2.2.2. User Interface

A Web-based user interface has been developed using the Yii2 PHP
framework777Yii2 website: https://www.yiiframework.com/. It comprises various
wizards that offer execution, reproducibility, and monitoring functionalities
for computational experiments (see also Section 2.3). Of course, the function
of all these wizards heavily relies on the functionalities provided by the
rest components, on which we elaborate in the next sections.

####  2.2.3. Software Manager

This component implements functionalities to upload (or update) container
images and workflows. First of all, any involved container should be loaded in
the _Docker Image Registry_. In addition, in both cases (i.e., individual
container or workflow), a CWL description is required and the corresponding
metadata, which describe the required inputs and dependencies of the involved
software packages (the latter only for workflows), are loaded in the _Software
Metadata DB_. These data are used by various components, e.g., they are used
by the _User Interface_ wizard to automatically display a form containing one
input field for each input parameter of the software (see also the example of
Section 2.2.2).

####  2.2.4. Job Submission API

Apart from manually submitting computational jobs through the UI wizards, the
users are able to also submit batches of jobs programmatically using an
implemented API. This API is based on GA4GH’s Task Execution Schemas
(TES)888TES specification: https://github.com/ga4gh/task-execution-schemas and
Workflow Execution Service (WES)999WES specification:
https://github.com/ga4gh/workflow-execution-service-schemas API
specifications. The API supports batch execution and monitoring of
computational tasks and can be used by user-implemented scripts and programs.

####  2.2.5. Software Profiler

This component leverages machine learning to produce (after user request)
“execution profiles” for software that has been uploaded on SCHeMa. For each
software, this profiler builds a classification model that attempts to map
candidate jobs of this software to one class of nodes, which corresponds to
the type of node that is appropriate to undertake the computations of the job.
As an indicative example, a cluster could have two types of nodes, one with
regular-sized main memory and another one with large memory (ideally, to be
used by memory-intensive tasks). In this case, given a particular software of
interest, the objective of _Software Profiler_ would be to train a (binary)
classification model to assign each job of this software to a regular-memory
node or to a large-memory one.

The profiling process goes as follows: first, the uploader of the software
provides a set of alternative values/files for its input parameters (a
relevant _User Interface_ wizard exists to collect this information). Then,
the system runs the software for all input combinations, collects the
resources consumption of each run, and creates a dataset of samples, where
each sample is the combination of provided inputs and the recorded resources
consumption. Based on this dataset, the system trains and optimizes different
classifier models (implemented by scikit-learn101010https://scikit-
learn.org/stable/) using the grid search approach. The model presenting the
best accuracy is selected and stored in the _Execution Profile DB_. All the
stored models are exploited by the _Job Classifier_ component, when the
execution of the corresponding software is requested.

####  2.2.6. Job Classifier

This component receives as input job submissions from the corresponding wizard
of the _User Interface_ and the _Job Submission API_. As a first step, it
searches if there are any trained (by the _Software Profiler_) models for the
involved software packages. If yes, it exploits this model to create a
suggestion for the most appropriate node type to undertake the job and then
propagates it to the Kubernetes scheduler through the Kubernetes API. The
scheduler takes into consideration this suggestion and schedules the
corresponding job ensuring that the determined node type will be used. If
there is no trained model, then the job request is propagated to the
Kubernetes shceduler without any indication for the node type to be used.

####  2.2.7. Data Repository Connector.

This component implements an interconnection with various open data
repositories. Currently, two repository services are supported:
Zenodo111111Zenodo open repository: https://zenodo.org/ and
HELIX121212Hellenic Data Service (HELIX) repository:
https://hellenicdataservice.gr. _Data Repository Connector_ takes advantage of
the APIs provided by the repository services to download/upload datasets
from/to them. A relevant _User Interface_ wizard exploits this functionality
to enable users to download existing datasets from one of the supported
repositories and then use them for their analyses or to directly upload the
output of a particular analysis on a selected data repository.

####  2.2.8. Experiment Packager.

This component is responsible to create “experiment packages” (according to
the RO-crate (Carragáin et al., 2019) specification) from previously executed
computational jobs, after user request. A relevant wizard is implemented in
the _User Interface_ and, based on it, the users can easily create packages
that incorporate several metadata for the selected experiment, such as the
software used, its configuration, the input and output dataset, or even the
DOI of a relevant publication. To collect the required information, the
_Experiment Packager_ component communicates with the _Job Statistics_ and the
_Software Metadata_ databases. The resulting packages are stored inside an
_Experiments DB_ but it is also possible for the users to download the
packages to their local computer storage. The easy creation of RO-crate
packages is one of the main SCHeMa functionalities that facilitate
computational experiment reproducibility.

####  2.2.9. Monitoring Component

This component aggregates data coming from the low-level logging and
monitoring mechanisms of Kubernetes to create insightful reports about the
jobs being executed in the cluster. All jobs, both those submitted through the
UI and those submitted through the API are being considered. It also provides
statistics about the load of the cluster. The component constantly
communicates with the _Job Statistics DB_ to update its recorded information
or to use it for the production of the aggregated statistics. Finally, it
propagates data to the job execution wizard so the user can monitor each job’s
output and status.

###  2.3. Quick Tour of the Interface

For the interest of space, in this section we indicatively describe only the
main functionalities around the execution of scientific containers (however,
the same functionalities for workflows are very similar). It should be
highlighted that it is practically impossible for all the offered
functionalities to be described.

A list of all containers in the _Docker Image Registry_ , for which the
connected user has access, can be displayed after clicking on the “Software”
top menu item of SCHeMa. In the top of Figure 2, a screenshot of SCHeMa’s
interface after this action took place is illustrated. By hitting the arrow-
shaped button of any entry, the user is redirected to the job submission
wizard of the corresponding software131313Hitting the diagram-shaped button
the user can start the machine learning profiling of the same software.. This
wizard consists of a form that contains one input field for each input
parameter of the software; the form is automatically generated based on the
CWL description of the selected software (which is stored in the _Software
Metadata DB_). After providing input values for all required fields and
hitting “Run” the execution starts and the progress can be monitored through
the UI. Closing the browser tab is possible without interrupting the
execution; the user may revisit the execution page for this job by selecting
the corresponding entry in the “Job History” page (accessible again through a
top menu item). In the same page, the user can also select to rerun a
previously completed execution or to create an RO-crate object based on it.
Finally, any output files, which are stored in the distributed storage space,
can be found in the “Data” page.

##  3\. Demonstration

During the conference, we will explain the concept of reproducible and data-
driven science and its requirements to the audience and we will demonstrate
SCHeMa’s relevant functionalities elaborating on how they help in this
context. For this demonstration, we will exploit SCHeMa’s deployment for the
ELIXIR-GR Cloud Infrastructure, which is based on a relatively large
computational cluster (see Section 1). We will examine SCHeMa’s capabilities
in real-time by following any audience-defined scenario, however we will also
demonstrate some interesting scenarios we have identified.

Figure 2. A screenshot of SCHeMa’s interface.

The main scenario, is based on executing a pre-loaded scientific container
using the corresponding wizard. We will run the container twice: once without
exploiting its pre-trained execution profile and once leveraging it and we
will prompt the users to observe any differences (e.g., without using the
profile a node with unnecessary large resources may be selected, thus its
“precious” CPU time may be wasted instead of being used for a more demanding
job). After the execution, we will also guide the audience through the process
of uploading the output files on an open data repository and packaging the
whole experiment into an RO-crate object.

##  4\. Conclusion

We introduced and demonstrated SCHeMa, an open-source platform that aims to
assist the work of scientists in the data-driven science era through
facilitating the execution, reproducibility, and monitoring of computational
experiments on heterogeneous computational clusters. To this end, it leverages
various technologies like containerization, experiment packaging, workflow
management, and machine learning.

###### Acknowledgements.

This work was partially funded by the “ELIXIR-GR: Managing and Analysing Life
Sciences Data (MIS: 5002780)” project (co-financed by Greece and the European
Union - European Regional Development Fund).

## References

  * (1)
  * Afgan et al. (2018) Enis Afgan, Dannon Baker, Bérénice Batut, Marius van den Beek, Dave Bouvier, Martin Cech, John Chilton, Dave Clements, Nate Coraor, Björn A. Grüning, Aysam Guerler, Jennifer Hillman-Jackson, Saskia D. Hiltemann, Vahid Jalili, Helena Rasche, Nicola Soranzo, Jeremy Goecks, James Taylor, Anton Nekrutenko, and Daniel J. Blankenberg. 2018.  The Galaxy platform for accessible, reproducible and collaborative biomedical analyses: 2018 update.  _Nucleic Acids Res._ 46, Webserver-Issue (2018), W537–W544.  https://doi.org/10.1093/nar/gky379
  * Baker (2016) Monya Baker. 2016\.  1,500 scientists lift the lid on reproducibility.  _Nature_ 533 (2016), 452–454.  https://doi.org/10.1038/533452a
  * Boettiger (2014) Carl Boettiger. 2014\.  An introduction to Docker for reproducible research, with examples from the R environment.  _CoRR_ abs/1410.0846 (2014).  arXiv:1410.0846 http://arxiv.org/abs/1410.0846
  * Carragáin et al. (2019) Eoghan Ó Carragáin, Carole Goble, Peter Sefton, and Stian Soiland-Reyes. 2019.  A lightweight approach to research object data packaging. In _Bioinformatics Open Source Conference (BOSC) 2019_.
  * Chamberlain and Schommer (2014) Ryan Chamberlain and Jennifer Schommer. 2014.  Using Docker to support reproducible research.  _DOI: https://doi. org/10.6084/m9. figshare_ 1101910 (2014), 44.
  * da Veiga Leprevost et al. (2017) Felipe da Veiga Leprevost, Björn A Grüning, Saulo Alves Aflitos, Hannes L Röst, Julian Uszkoreit, Harald Barsnes, Marc Vaudel, Pablo Moreno, Laurent Gatto, Jonas Weber, et al. 2017\.  BioContainers: an open-source and community-driven framework for software standardization.  _Bioinformatics_ 33, 16 (2017), 2580–2582.
  * Fanelli (2018) Daniele Fanelli. 2018\.  Opinion: Is science really facing a reproducibility crisis, and do we need it to?  _Proceedings of the National Academy of Sciences (PNAS)_ 115, 11 (2018), 2628–2631.  https://doi.org/10.1073/pnas.1708272114
  * Jimenez et al. (2015) Ivo Jimenez, Carlos Maltzahn, Adam Moody, Kathryn Mohror, Jay Lofstead, Remzi Arpaci-Dusseau, and Andrea Arpaci-Dusseau. 2015.  The role of container technology in reproducible computer systems research. In _2015 IEEE International Conference on Cloud Engineering_. IEEE, 379–385.
  * Wolstencroft et al. (2015) Katherine Wolstencroft, Stuart Owen, Olga Krebs, Quyen Nguyen, Natalie J. Stanford, Martin Golebiewski, Andreas Weidemann, Meik Bittkowski, Lihua An, David Shockley, Jacky L. Snoep, Wolfgang Müller, and Carole Goble. 2015\.  SEEK: a systems biology data and model management platform.  _BMC Systems Biology_ 9, 1 (11 Jul 2015), 33\.  https://doi.org/10.1186/s12918-015-0174-y