# Shapes of a filament on the surface of a bubble

S Ganga Prasath<EMAIL_ADDRESS>School of Engineering and
Applied Sciences, Harvard University, Cambridge MA 02138, USA..  International
Centre for Theoretical Sciences (ICTS-TIFR) Shivakote,
Hesaraghatta Hobli, Bengaluru 560089, India.  Joel Marthelot  Aix-Marseille
University, CNRS, IUSTI (Institut Universitaire des Systémes Thermiques
Industriels), 13013 Marseille, France.  Rama Govindarajan  International
Centre for Theoretical Sciences (ICTS-TIFR) Shivakote,
Hesaraghatta Hobli, Bengaluru 560089, India.  Narayanan Menon  Department of
Physics, University of Massachusetts Amherst, Amherst, MA 01003, USA.

###### Abstract

The shape assumed by a slender elastic structure is a function both of the
geometry of the space in which it exists and the forces it experiences. We
explore by experiments and theoretical analysis, the morphological phase-space
of a filament confined to the surface of a spherical bubble. The morphology is
controlled by varying bending stiffness and weight of the filament, and its
length relative to the bubble radius. When the dominant considerations are
geometry of confinement and elastic energy, the filament lies along a geodesic
and when gravitational energy becomes significant, a bifurcation occurs, with
a part of the filament occupying a longitude and the rest along a curve
approximated by a latitude. Far beyond the transition, when the filament is
much longer than the diameter, it coils around the selected latitudinal
region. A simple model with filament shape as a composite of two arcs captures
the transition well and for better quantitative agreement with the subcritical
nature of bifurcation, we study the morphology by numerical energy
minimization. Our analysis of filament’s morphological space spanned by a
geometric parameter, and one that compares elastic energy with body forces,
may provide guidance for packing slender structures on complex surfaces.

## Introduction

The spatial conformation of a slender elastic structure is a function both of
the geometry of the space in which it exists as well as the forces it
experiences. An unconfined elastic filament does not bend or stretch, but when
the elastic filament is confined to a two dimensional surface, the geometrical
properties of the surface plays a role in determining its morphology Huynen et
al. (2016); Guven et al. (2014). Such surface confinement can arise in a
variety of situations such as a DNA wrapping around histone Swigon et al.
(1998); Tobias et al. (2000); Coleman et al. (2000), plants climbing along
complex topographies Goriely and Neukirch (2006), packaging of optic fibre
cables around curved surfaces Seemann (1996) to long drill strings beneath the
earth’s surface Van der Heijden (2001); Tan and Digby (1993); Wu et al.
(1993); Wu and Juvkam-Wold (1993). The role of surface geometry in determining
the shape is crucial in all of these situations, however in each example the
filaments are also subject to forces (electrostatic, active internal stresses,
gravity, friction) that affect their conformation. For example a filament with
intrinsic curvature held at one of its ends, such as a hair Miller et al.
(2014), takes shapes determined by the competition between its own weight and
the bending stiffness of the structure Mahadevan and Keller (1996, 1999),
where for small mass density the filament retains its natural shape, and for
denser filaments, the structure undergoes a shape transition to a morphology
determined by the weight.

In this article we explore the general question of how substrate geometry and
external force combine to determine the shape of a filament by considering the
specific case of an elastic filament on the surface of a spherical bubble.
When a pre-wet filament is introduced on the surface of a droplet or a bubble,
if the droplet’s Laplace pressure exceeds a critical value, the filament
buckles Elettro et al. (2016); Bico et al. (2004); Prasath et al. (2021); Brau
et al. ; Schulman et al. (2017). In the buckled state the shape of the drop
can change in order to minimize the surface energy of the droplet and the
elastic energy of the filament Py et al. (2007); Prasath et al. (2021).
However, since we wish to explore the mechanics of the filament in a fixed
confining geometry, we work with highly-bendable filaments so that the
spherical geometry of the substrate remains unchanged by the filament. In this
high-bendability limit, we perform experiments where the forces due to gravity
and bending compete to determine the morphology of a filament confined to the
surface of a spherical pendant bubble. When the gravitational potential energy
of the filament is negligible compared to its bending energy, we find that the
filament lives along a geodesic – a longitude of the sphere – where it adopts
the radius of curvature of the bubble. However when the effects of gravity
become important, the filament undergoes a transition into a complex shape.
Beyond this transition, the filament lies partly along a longitude and partly
along a latitude at a particular polar angle. A calculation which accounts for
the energy of these two segments of the filament correctly captures the shape
below the transition length as well as the threshold control parameters for
bifurcation from the geodesic shape. However, the approximate theory fails to
capture the quantitative details of the shape after the transition. A
numerical minimization of the competing energies in the system identifies
deviations from the approximate theory and captures accurately the shapes in
experiments. On further increase in the length of the filament, the filament
enters a coiled phase with the filament aligned on top of itself in layers. In
what follows, we explore experimentally and by numerical simulation, a two-
parameter phase diagram of the morphology of the filament in a curved space:
one parameter is geometric, and describes the length of the filament versus
the radius of confinement; the other parameter encapsulates the competition
between elasticity and the external forces (gravity, in this case).

Figure 1: $(a)$ A soap bubble hanging from the end of a capillary tube of
diameter $1{\rm mm}$ with gravity pointing downwards. The elastic filament
sits on top of the bubble with one end hinged to the capillary tube. $(b)$
Schematic of the setup in which the bubble is approximated by a sphere with
$u^{(1)}$ being the azimuthal angle, and $u^{(2)}$ being the polar angle on
the sphere. Gravity points downwards, towards the negative $z$-axis. The
filament shown in red is represented using $(u^{(1)}(s),u^{(2)}(s))$ where $s$
is the arc-length along the filament.

## Experiments and length-scales

In our experiments the spherical surface is provided by a pendant soap bubble
hanging at the end of a capillary tube as shown in Fig. 1$(a)$. The bubble is
made using a liquid soap (commercial DAWNTM) in a 1:4 water-glycerol solution.
The high concentration of glycerol with viscosity 1000 cSt helps reduce the
drainage rate providing longer bubble survival time. We then place a thin
elastic filament made out of silicone on top of the bubble. These filaments
are made by first melting a rod of silicone at $300^{\circ}$C using a silicone
glue gun (Aptech Crown) and then pulling a small droplet of the melt with
tweezers which then sets in a few seconds. The diameters of the filaments made
by this procedure can be varied between
$\mathcal{O}(10\mu\text{m})-\mathcal{O}(100\mu\text{m})$ by varying the
pulling rate manually. The diameter is measured under an optical microscope,
where we also confirm that we have a filament of uniform thickness. We also
soak the filaments in Sudan Red G dye (SIGMA 17373) solution and dry them to
make them visible in white light. The bubble is illuminated by a diffuse
background light and we use a Nikon DLSR 5000 camera to observe the
conformation of the filament on the bubble. We wet the filaments in the soap
solution before placing them on the bubble surface in order to eliminate the
effects of the filament surface energy on the filament morphology. This
procedure also helps us maintain longer bubble survival times as a dry
filament can cause sudden changes in the surface energy of the interface
resulting in rupture. When the wet filament is placed on the bubble, the
change in the surface energy of the bubble is only due to distortion of the
fluid film around the thickness of the filament. In the experiments one end of
the filament is held at the end of the capillary tube and the other end is
left free. The hinged end of the filament experiences no torque and can rotate
freely about the end of the capillary tube. For a given filament of length
$L$, we quasistatically increase or decrease the bubble radius $R_{b}$ using a
microfluidic pump attached to the other end of the capillary tube.

Figure 2: $(a-f)$ Filament morphology on the surface of the bubble for
different bubble sizes for a fixed filament length of $L=30{\rm mm}$. For
large bubble size the filament stays along a longitude and as the bubble size
is decreased, assumes more complex shapes where it bends smoothly from the
longitude along a chosen latitude. The solid line is traced along the filament
for higher contrast, and is shown in orange color along the longitude, and in
black where it deviates from a longitude. In $(f)$ we show the filament in
self-contact. Coiling occurs on further decrease in bubble size. The scale bar
is $3{\rm mm}$.

In order to describe the phenomenon we observe in our experiments, we collect
in this section the relevant length-scales that play a role in the determining
the filament morphology. There are several length-scales associated both with
the filament and the bubble; we attempt to reduce the complexity by working in
regimes where some of these are irrelevant. Since we expect gravity to be
relevant in our experiments, we define the elasto-gravity length-scale
$l_{{\rm eg}}\sim(2EI/(\varrho g\pi t^{2}))^{1/3}$, where $E$ is the Young’s
modulus of the filament material, $I$ the area moment of inertia of the
filament, $\varrho$ the material density, $t$ the filament thickness and $g$
the gravitational constant. This is the length scale associated with filament
deformation when the weight of the filament is balanced by bending forces. In
the limit of strong capillarity, where the bubble does not deform, the radius
of the bubble $R_{b}$ is the only length-scale associated with the 2D
spherical surface. Using these two length-scales we can define a non-
dimensional number that quantifies the relative importance of gravity and
bending as the elasto-gravity bendability, $\Omega^{-1}_{g}=(l_{{\rm
eg}}/R_{b})^{3}=2EI/(\varrho g\pi t^{2}R_{b}^{3})$. The larger the value of
$\Omega_{g}$, the stronger the effects of gravity and the smaller the
$\Omega_{g}$, the stronger the effects of bending. Another non-dimensional
number arising purely out of geometry is $\Phi=(L/R_{b})$. We call this non-
dimensional number the coiling parameter, since as we will see later this
number helps describe the state of coiling. The two non-dimensional numbers of
interest thus are the coiling parameter $\Phi$ and the elasto-gravity
bendability $\Omega^{-1}_{g}$ and the different morphologies of the elastic
filament on the spherical bubble are described in this phase space spanned by
these two nondimensional groups.

In identifying the relevant length-scales in our experiments, we have tacitly
made assumptions of scale separation in the filament-bubble system that we
detail here. For a spherical bubble made of a fluid of surface tension
$\gamma$, fluid density $\varrho_{f}$, radius $R_{b}$ and a filament of length
$L$, Young’s modulus $E$, material density $\varrho$, and thickness $t$, there
are six length-scales that at play in the experiments. These are $L,R_{b},t$,
capillary length $l_{c}\sim\sqrt{\gamma/\varrho_{f}g}$, bendo-capillary length
$l_{{\rm bc}}\sim(EI/\gamma)^{1/3}$ and finally, the elasto-gravitational
length scale $l_{{\rm eg}}$ described above. Since we want the bubble to
remain spherical throughout the experiments, we ensure this by keeping the
capillary length to be the largest length scale in the system. Moreover we
want the surface energy of the bubble to not affect filament bending i.e., we
want the bendo-capillary length to be very small such that the filament is
able to completely lie on the bubble without affecting the bubble. This is
taken care of by operating in a scale-separated regime: $t\ll l_{{\rm
bc}}\ll(L,R_{b},l_{{\rm eg}})\ll l_{c}$. When the relevant length-scales are
$(L,R_{b},l_{{\rm eg}})$, the filament morphology is governed by the geometry,
the bending force and the self-weight of filament. This separation of scales
in our experiments is established by choosing a thin
($t\sim\mathcal{O}(100\mu{\rm m})$) and soft filament ($E=1\text{MPa}$) with a
bubble size much larger than filament diameter ($L,R_{b}\sim\mathcal{O}({\rm
cm})$) leading to an bendo-capillary length-scale $l_{{\rm
bc}}\sim\mathcal{O}({\rm mm})$. The ratio of capillary force and bending force
can be quantified using inverse capillary bendability, $\Omega^{-1}=(l_{\rm
bc}/R_{b})^{3}$. Lastly, in our experiments the tangential strain along the
filament is small and thus are in the inextensible elastica limit Audoly and
Pomeau (2010). This is ensured by choosing
$l_{m}=\gamma/E\sim\mathcal{O}(\textrm{nm})\ll t$. We explore the phase space
of filament morphologies by changing the filament length $L$, the bubble
radius $R_{b}$ and the filament thickness $t$ while keeping all the other
experimental parameters fixed.

## Results

When one end of the filament is hinged to the capillary tube at the north pole
of the bubble, with gravity pointing from north to south, the filament takes
different morphologies, as shown in Fig. 2, when the ratio $\Phi$ of filament
length to bubble radius is increased. Initially, the filament lives along a
longitude (see Fig. 2$(a)$). For small filament lengths we expect gravity to
not play a role in determining the filament configuration, so the shape is
determined purely by minimizing the bending energy. For a sphere the great
circle has the smallest curvature and thus the filament lies along a
longitude, since bending energy $\mathscr{E}\sim\varkappa^{2}L$, $\varkappa$
being the curvature. As we increase the ratio $\Phi$, the conformation of the
filament changes from that of a longitude, as shown in Fig. 2$(b-d)$.
Intuitively we recognise that this effect is due to gravity: the weight of the
filament pulls it towards the bottom of the sphere. At the largest values of
$\Phi$, the filament coils around a fixed latitude in the bottom half of the
sphere.

We explore the phase space of morphologies, spanned by the two variables
$\Phi$ and $\Omega^{-1}_{g}$, to ensure that the physical mechanism we
describe is indeed what we observe in experiments. For a fixed length of the
filament hinged at the north pole a change in the radius of the bubble $R_{b}$
changes both variables, causing us to traverse the morphological phase-space
along trajectories $\Phi\sim 1/R_{b}$ and $\Omega^{-1}_{g}\sim 1/R_{b}^{3}$.
We perform these experiments with eleven different filament lengths $L=0.5{\rm
cm}-5{\rm cm}$ and two different filament thicknesses $t=70\mu{\rm
m},100\mu{\rm m}$. We show these trajectories for different filament lengths
in Fig. 3$(a)$ by dashed lines. These trajectories and configurations can be
reversed by re-inflating the bubble, except in the regime where the filament
coils and self-contact of the filament leads to irreversibility in
configurations. Along a given trajectory the green dots correspond to a
filament configuration along the longitude, and red dots correspond to
deviations from the longitude to more complex shapes. Since changing the
thickness of the filament changes the bending stiffness ($EI\sim t^{4}$) and
changing bubble size changes the elasto-gravity bendability
($\Omega^{-1}_{g}\sim 1/R_{b}^{3}$), we are able to span two orders of
$\Omega^{-1}_{g}$ in the experiments.

The experimental observations are qualitatively explained by our arguments
regarding the role of the geometry ($\Phi$) and gravitational potential energy
versus elastic energy ($\Omega_{g}$). However, a quantitative description of
the nature of the transition from the longitude is desired. In such a
description, we would also like to predict the polar position at which the
filament chooses to situate itself beyond the transition point from a
longitude to a complex shape. In order to gain a quantitative understanding we
approximate the shape of the filament as an elastic curve on a sphere and
calculate the transitions between morphologies by minimizing the total energy
due to bending deformations as well as gravity. We validate this with the
experiments and gain further understanding by performing a numerical
minimization of the total energy by approximating the filament as a
combination of discrete elastic rods.

### Geometry and mechanics of the filament

We describe the morphology of the filament in our experiments using the
elastic energy and the gravitational potential energy of the filament. Since
stretching and twisting are extremely high-energy deformations, we consider
only bending energy in the inextensible limit. The total energy of the
filament in this approximation is

$\displaystyle\mathscr{E}=\frac{EI}{2}\int_{0}^{L}\varkappa^{2}(s)\
\textrm{d}s+\varrho g\frac{\pi
t^{2}}{4}\int_{0}^{L}(\bm{\mathcal{S}}(s)-\bm{\mathcal{S}}_{o})\cdot\hat{\bm{z}}\
\textrm{d}s,$ (1)

where $\bm{\mathcal{S}}(s)\equiv\bm{\mathcal{S}}(u^{(2)}(s),u^{(1)}(s))$ is
the location of the filament center-line on the surface of the sphere.
$\bm{\mathcal{S}}_{o}$ is a reference location about which the potential
energy is measured, $\varkappa(s)$ the magnitude of curvature along the arc-
length $s$ and $\hat{\bm{z}}$ the direction against gravity (see Fig. 1$(b)$).
For the specific case of a sphere we can explicitly write the parameterization
of the surface as:

$\bm{\mathcal{S}}(u^{(2)},u^{(1)})=w(\cos u^{(2)}\cos u^{(1)},\cos u^{(2)}\sin
u^{(1)},\sin u^{(2)}).$

The curvature vector $\bm{\varkappa}(s)$ can be decomposed in the orthogonal
Darboux frame as:

$\displaystyle\bm{\varkappa}(s)=$ $\displaystyle\
\varkappa_{n}\bm{\hat{N}}+\varkappa_{g}(\bm{\hat{N}}\times\hat{\bm{d}}_{3}),$
(2)

$\varkappa_{n}$ being the normal curvature along the surface normal given by

$\bm{\hat{N}}=\frac{\bm{\mathcal{S}}_{u^{(1)}}\times\bm{\mathcal{S}}_{u^{(2)}}}{||\bm{\mathcal{S}}_{u^{(1)}}\times\bm{\mathcal{S}}_{u^{(2)}}||}.$

where $\bm{\mathcal{S}}_{u^{(j)}}$ denotes derivative with respect to
$u^{(j)},j=1,2$ and $\varkappa_{g}$ is the geodesic curvature along the
binormal direction $(\bm{\hat{N}}\times\hat{\bm{d}}_{3})$, with
$\hat{\bm{d}}_{3}$ being the tangent vector along the curve
$\bm{\mathcal{S}}(s)$. The arc-length and the curvature are non-
dimensionalized using the bubble radius $R_{b}$ and the total energy is non-
dimensionalized using $(EI/2R_{b})$ as the energy scale:
$\tilde{s}=s/R_{b},\tilde{\varkappa}=\varkappa R_{b}$ and
$\tilde{\mathscr{E}}=\mathscr{E}/(EI/2R_{b})$. For a sphere we know that the
normal curvature is a constant everywhere, $\varkappa_{n}=(1/R_{b})$ and using
this we can rewrite the energy expression as:

$\displaystyle\tilde{\mathscr{E}}$
$\displaystyle=\int_{0}^{\Phi}[\tilde{\varkappa}_{g}^{2}+\Omega_{g}(\tilde{\bm{\mathcal{S}}}(\tilde{s})-\tilde{\bm{\mathcal{S}}}_{o})\cdot\hat{\bm{z}}]\
\textrm{d}\tilde{s}+\Phi.$ (3)

where tildes denote non-dimensional variables.

Figure 3: $(a)$ Phase diagram of coiling represented by the non-dimensional
elasto-gravity bendability, $\Omega^{-1}_{g}$ and the coiling parameter,
$\Phi$. Circles are experimental data points for a given combination of
$\Phi,\Omega^{-1}_{g}$. The green region is when the filament stays along the
longitude, separated from the red region by the solid line representing the
theoretical critical coiling parameter $\Phi_{c}(\Omega^{-1}_{g})$, which
represents the regime of deviations from longitude to the more complex shapes
seen in Fig. 2. $\Phi^{*}$ is the coiling parameter at which coiling starts
for $\Omega^{-1}_{g}\rightarrow+\infty$ and we see that
$\Phi_{c}\rightarrow\Phi^{*}$ for large $\Omega^{-1}_{g}$. The white region,
with orange experimental points is the regime of coiling. $(b)$ Curvature
along the filament, $\tilde{\varkappa}(\tilde{s})$ as a function of arc-
length, $\tilde{s}$ from numerical solutions for different values of $\Phi$
(shown in the legend) for a fixed $\Omega^{-1}_{g}=0.1$. $(c)$ Bifurcation
diagram of coiling: Fraction of filament length along longitude, $\alpha$ as a
function of the coiling parameter for fixed values of $\Omega^{-1}_{g}$.
Symbols are from experiments and solid lines (dashed lines) from theory by
minimizing energy expression in Eq. 9 with (without) the constraint of
$\alpha\Phi<\pi$ for different $\Omega^{-1}_{g}$ values shown in the inset.
The data points in the gray region indicate that the length along the
longitude never crosses the south pole. Numerical minimization shown in black
for $\Omega^{-1}_{g}=0.1$ captures the experimental observation accurately. We
also show the critical coiling parameter from the model, $\Phi_{c}$, in
turquoise. Inset also shows different shapes we obtain from numerics for
different filament lengths.

### Limit of zero gravity

In the asymptotic limit of negligible gravity, $\Omega_{g}\rightarrow 0$, the
filament morphology is determined only by the bending energy and the solution
amounts to minimizing the bending energy:

$\tilde{\mathscr{E}}\approx\int_{0}^{\Phi}\tilde{\varkappa}_{g}^{2}\
\textrm{d}\tilde{s}+\Phi.$ (4)

We know that $\Phi$ remains fixed for a given filament length, $L$ and bubble
radius, $R_{b}$ and thus $\tilde{\varkappa}_{g}=0$ or essentially any geodesic
is a solution. Since we have hinged one end of the filament, this allows only
for longitudes: geodesics that travel from the north pole to the south pole on
the bubble. In the weak gravity limit, the boundary condition at the free end,
where the filament experiences a torque $\sim EI/R_{b}$, does not modify the
trajectory from the geodesics (as also pointed out in Huynen et al. (2016);
Guven et al. (2014)) which for a sphere are the great circles. Thus the
filament takes this shape for $\Phi\leq 2\pi$. When $\Phi>2\pi$, we expect
coiling along a geodesic, with deviations from the geodesic due to self-
intersections at finite $t/R_{b}$. We define the coiling parameter $\Phi^{*}$
at which coiling starts as the filament length at which self-intersection
occurs. For $\Omega^{-1}_{g}\rightarrow\infty$ the non-dimensional length for
coiling is $\Phi^{*}\rightarrow 2\pi$.

### Finite gravito-bendability effects

In the experiments we have non-zero values for elasto-gravity bendability and
as a consequence the filament shape deviates from the geodesic and the
critical coiling parameter $\Phi_{c}$ at which this occurs is expected to be a
function of $\Omega^{-1}_{g}$. The complex shape we see in the experiments
beyond the transition from a geodesic can be simplified by approximating it to
be a combination of a longitude and a latitude (which is no more a geodesic)
at a given polar angle. This approximation ignores the smooth transition
between latitude and longitude, treating it geometrically as a sharp kink,
with no elastic energy cost. (We evaluate the effects of this approximation
later in this article via numerical minimization.) The total energy of such a
configuration is given by

$\displaystyle\tilde{\mathscr{E}}=$ $\displaystyle\
\tilde{\mathscr{E}}_{l}(\alpha\Phi)+\tilde{\mathscr{E}}_{g}((1-\alpha)\Phi),$
(5)

where $\alpha$ is the fraction of the filament along the longitude,
$\tilde{\mathscr{E}}_{l}$ is the energy due to the filament along the
longitude and $\tilde{\mathscr{E}}_{g}$ is the energy due to the section along
the latitude. The parameterization of the filament shape lying along any
latitude with polar angle $\xi$ is given by,
$u^{(2)}=\xi,u^{(1)}=\tilde{s}/\lambda,\tilde{\bm{\mathcal{S}}}=(\cos\xi\
\cos(\tilde{s}/\lambda),\cos\xi\ \sin(\tilde{s}/\lambda),\sin\xi).$ The
inextensibility constraint is enforced by setting $\lambda=\pm\cos\xi$. The
energy of this configuration can be calculated to be,

$\displaystyle\tilde{\mathscr{E}}_{g}((1-\alpha)\Phi;\Omega_{g})$
$\displaystyle=\ (1-\alpha)\Phi\Omega_{g}(\sin\xi+z_{o})$ (6) $\displaystyle\
+\int_{0}^{(1-\alpha)\Phi}\tilde{\varkappa}_{g}^{2}\
\textrm{d}\tilde{s}+(1-\alpha)\Phi,$ $\displaystyle=\
(1-\alpha)\Phi\big{[}1+\tan^{2}(\xi)+\Omega_{g}(\sin\xi+z_{o})\big{]}.$ (7)

We have used the fact that the geodesic curvature of any latitude for a
spherical metric as $\tilde{\varkappa}_{g}=\tan\xi$ Kreyszig (2019). Along a
longitude the potential energy of any filament configuration depends on where
the starting point of the filament $\tilde{s}=0$ lies. We can denote an
arbitrary longitude by the parameterisation:
$u^{(2)}=\tilde{s}-\vartheta,u^{(1)}=0$ where $\vartheta$ is the polar angle
at which the filament is hinged, which we will set to $-\pi/2$ denoting the
north-pole. We can now write the potential energy of the configuration as:
$\tilde{\mathscr{E}}_{l}^{p}=\Omega_{g}[\cos\vartheta-\cos(\vartheta-\Phi)+z_{o}\Phi]$.
Since all longitudes are geodesics, $\tilde{\varkappa}_{g}=0$ along these
curves, leading to the total energy,

$\displaystyle\tilde{\mathscr{E}}_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}l}}$
$\displaystyle=\int_{0}^{\alpha\Phi}\Omega_{g}(\mathbf{x}-\mathbf{x}_{o})\cdot\hat{\bm{z}}\
\textrm{d}\tilde{s}+\int_{0}^{\alpha\Phi}\tilde{\varkappa}_{g}^{2}\
\textrm{d}\tilde{s}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\alpha}\Phi,$
$\displaystyle=\alpha\Phi(1+\Omega_{g}z_{o})+\Omega_{g}\big{[}\cos\vartheta-\cos(\vartheta-\alpha\Phi)\big{]}.$
(8)

Lastly the geometric constraint from inextensibility relates $\alpha$ and
$\xi$ as: $\alpha\Phi=(\pi/2-\xi)$ and after setting $\vartheta=-\pi/2$, the
resultant energy can be written as:

$\displaystyle\tilde{\mathscr{E}}=$ $\displaystyle\
\alpha\Phi+\Omega_{g}\sin(\alpha\Phi)+(1-\alpha)\Phi\big{[}\csc^{2}(\alpha\Phi)+\Omega_{g}\cos(\alpha\Phi)\big{]}.$
(9)

The equilibrium configuration is given by extremising the one parameter energy
expression 9, achieved by solving $\delta\tilde{\mathscr{E}}/\delta\alpha=0$
for $\alpha\in[0,1]$. Before we minimize this expression, we must regularize
the unphysical divergence in energy $\tilde{\mathscr{E}}$ for
$\alpha\Phi\rightarrow\pi$, which occurs due to the absence of a cutoff length
scale associated with capillary bendability. This is done by introducing a
cut-off length associated with the capillary bendability (see SI sec. S2 for
details). The phase boundary for deviation from the geodesic shape is
represented by the critical coiling parameter $\Phi_{c}$ when $\alpha$ first
deviates from $1$. We plot this in Fig. 3$(a)$ as the solid curve and find
that our simple geometric theory captures the transition seen in experiments
accurately, both at small values of $\Omega^{-1}_{g}$ as well as
$\Omega^{-1}_{g}\sim\mathcal{O}(1)$. We see that the increase in the value of
critical coiling parameter $\Phi_{c}$ with increase in $\Omega^{-1}_{g}$ due
to reduction in effects of gravity is also captured by the model. Moreover the
trend is consistent with the prediction that $\Phi_{c}\rightarrow\Phi^{*}$
when $\Omega^{-1}_{g}\rightarrow\infty$.

Figure 4: $(a-d)$ Sequence of images in the coiling phase beyond the point of
self-intersection where we see the filament starting to pack on the surface of
the bubble. For a very small bubble size as in $(d)$ the curvature of the
transition zone starts interacting with the hinge at the tip of the capillary
tube. Scale bar is $3{\rm mm}$. $(e)$ Self-intersecting shape of a filament in
numerical simulation beyond the transition threshold for coiling when $\Phi=7$
for $\Omega^{-1}_{g}=0.1$.

The transition we see in our experiments in Fig. 2 from a longitude happens
when $\alpha$ changes its value as we tune $\Phi$ for a fixed value of
$\Omega^{-1}_{g}$, however when we change the bubble radius we modify both
$\Phi$ as well as $\Omega^{-1}_{g}$. To compare the value of $\alpha$ from our
experiments with the calculation at fixed $\Omega^{-1}_{g}$, we hold the
bubble size fixed while we increase the filament length on the surface by
feeding the filament from the top of the bubble. We perform experiments using
this protocol for five different values of $\Omega^{-1}_{g}$ between
$0.04-0.1$ and plot the results in Fig. 3$(c)$. As the bubble is spherical, we
are able to calculate the length of the part of the filament along the
longitude from the 3D coordinates of the point where it deviates from a
longitude. In Fig. 3$(c)$, symbols correspond to experimental measurements,
dot-dashed lines for the theoretical minimum energy shapes and solid coloured
lines for the minimal energy shape with $\alpha\Phi$ constrained to be less
than $\pi$, i.e., the filament cannot cross the south pole. For small $\Phi$,
the filament lies perfectly along a longitude, i.e., along $\alpha=1$. This
holds true till just beyond $\Phi=\pi$ (dashed vertical line), where the
filament is long enough to touch the south pole. At higher $\Phi$, both
experiments and theory show deviations from a longitudinal shape, i.e.
$\alpha<1$. We emphasize that the theory corresponds to a part of the filament
lying perfectly on a longitude and the remainder lying perfectly on a
latitude, with a sharp transition between the two segments with no elastic
cost. In experiments, however, the transition between these two segments
occurs smoothly, and the deformation contributes to elastic energy. Despite
this missing energy at the transition zone, the theory predicts the $\Phi_{c}$
well. However for $\Phi>\Phi_{c}$, we see deviations in the values of $\alpha$
between the experiments and the theory. This deviation in magnitude of
$\alpha$ is evident in Fig. 3$(c)$ as all $\alpha$ values from experiments lie
below the $\alpha\Phi=\pi$ line while the theoretical predictions of minimal
energy without additional constraint on $\alpha$ (dashed colored lines) lie
above this line. But providing an additional constraint of $\alpha\Phi<\pi$ to
Eq. 9 (solid colored lines) results in reasonable agreement with experiments,
despite the simplifying assumption made at the transition zone. Since the
theory neglects the deformation at the transition zone between the geodesic
and the latitude, accounting for this missing length scale would reduce the
length of filament on the longitude. This thus would reduce $\alpha$,
providing better agreement with experiment without explicitly adding an
additional constraint. In order to understand the transition zone, we perform
numerical simulations of the filament shapes using the energy in Eq. 3.

### Numerical simulation of filament shapes

We explore the different shapes that the filament takes on the bubble by
numerically minimizing the energy in the Eq. 3 with the filaments segmented
into discrete linked rods. The total energy of the system can be written as

$\displaystyle\tilde{\mathscr{E}}=$ $\displaystyle\
\sum_{j}[\tilde{\varkappa}_{j}^{2}+\Omega_{g}z(\tilde{s}_{j})]h,$ (10)

where
$\tilde{\varkappa}_{j}=(\hat{\bm{t}}(\tilde{s}_{j})-\hat{\bm{t}}(\tilde{s}_{j+1}))/h$
and $h$ is the length of the discrete rods. Along with this, we further have
the constraint that one end of the filament is fixed:
$\hat{\mathbf{x}}(0)=(0,0,1)$ and that all the coordinates of the filament lie
along the surface of the sphere: $||\hat{\mathbf{x}}(\tilde{s}_{j})||=1$. We
use the fmincon function in Matlab’s optimization toolbox to find the minimal
energy shape of the filament. We show in Fig. 3$(c)$ inset the shapes we
obtain from the minimization and in Fig. 3$(b)$ the curvature of the filament
for different values of $\Phi$ when $\Omega^{-1}_{g}=0.1$. The first part of
the filament in the numerical solution has constant curvature, corresponding
to it lying entirely along a longitude, and beyond some arc length, deviates
from this constant curvature. The region deviating from
$\hat{\varkappa}\approx 1$ corresponds to the transition zone where effects of
the length scale missing in the theory are now accounted for in the numerics.
We compute the value of $\alpha$ as a function of $\Phi$ where we define
$\alpha$ as the length from the location where $\tilde{\varkappa}(\tilde{s})$
reaches its first minimum to the end of the filament. This is shown as black
circles in Fig. 3$(c)$; we find that the bifurcation is indeed sub-critical
and for larger values of $\Phi$ we find that $\alpha\sim 1/\Phi$ away from the
bifurcation point. The sub-critical nature of the transition can be attributed
to a similar mechanism as that seen in the wrapping of a droplet by an elastic
sheet (see ref. Brau et al. ; Kusumaatmaja and Lipowsky (2011); Py et al.
(2007)) where the external body force that the droplet applies on the sheet
makes the flat state unstable and the sheet goes to a wrapped state. The
effect of gravity breaks the up-down symmetry of the filament configuration
and thus leads to a state where the preferred shape is no more along the
longitude. We also look at the bifurcation when the ends of the filament are
free in SI sec. S1 and find a similar instability mechanism.

Since the contribution to the total energy of the filament comes from bending
and potential energy along the filament latitude and longitude parts as well
as the energy of the transition zone, we find the scaling of each of these
components. The energy contribution due to bending scales as
$\mathscr{E}_{l}\sim EIl/R_{b}^{2}$, gravitation to be
$\mathscr{E}_{g}\sim\varrho gt^{2}l^{2}$ and bending in the transition goes as
$\mathscr{E}_{t}\sim EIl_{{\rm eg}}(1/R_{b}^{2}+1/l_{\rm eg}^{2})$. These can
be written in non-dimensional terms in with $EI/R_{b}$ as the energy scale to
get:

$\tilde{\mathscr{E}}_{l}\sim\Phi,\
\tilde{\mathscr{E}}_{g}\sim\Omega_{g}\Phi^{2},\
\tilde{\mathscr{E}}_{t}\sim(\Omega_{g}^{1/3}+\Omega_{g}^{-1/3}).$

Though it might seem at first glance that the transition zone energy
$\tilde{\mathscr{E}}_{t}$ has singular contribution since it diverges in both
limits of $\Omega^{-1}_{g}\rightarrow 0,\infty$, this is however not the case
as the transition zone exists only when
$\Omega_{g},\Omega^{-1}_{g}\sim\mathcal{O}(1)$. It is also evident from the
scaling above that the critical length for transition
$\Phi_{c}\sim\Omega^{-1}_{g}$ (which is also the case for filaments with free
ends as shown in SI sec. S1).

### Coiled phase

In the coiled phase the packing fraction of the surface of the bubble is
independent of $\Omega^{-1}_{g}$ and is a purely geometric quantity determined
by the length, thickness of the filament and the bubble radius. In Fig. 4
$(a-d)$ we show the shape of coils in the experiments where we see that the
filament after self-contact forms complex shapes but at small bubble radius,
aligns with itself, forming multiple coils. When the size of bubble becomes
very small, the boundary layer close to the fixed end of the filament
determines the overall orientation and this results in the tilt of the coil
seen in Fig. 4$(d)$. However in the numerical simulations, in Fig. 4$(e)$, as
there is no interaction between filament points, we see that they self-
intersect and are all now localized at the bottom of the bubble.

## Conclusion

Our results can be posed within the broader framework of confinement and
bendability of elastic structures. Confinement prescribes the geometric
constraints placed on the elastic object, which in this case is the constraint
of conforming to a spherical substrate of given radius. Bendability determines
the competition between forces due to elastic deformation and external body
forces, gravity in this case. Elastic instabilities in thin sheets King et al.
(2012); Schroll et al. (2013); Davidovitch et al. (2011) have been explored
within this framework. The experiments presented here show that this paradigm
is also applicable to filaments: the coiling parameter, which depends
explicitly on the curvature of the bubble, is a confinement parameter, while
the elasto-gravity bendability quantifies the relative magnitude of the forces
due to gravity and bending deformation. We have only explored spherical
substrates in this work, however, other substrate geometries such as
negatively curved surfaces, and surfaces with non-uniform curvature could lead
to different instabilities and morphologies.

## References

  * Huynen et al. [2016] Alexandre Huynen, Emmanuel Detournay, and Vincent Denoël.  Surface constrained elastic rods with application to the sphere.  _Journal of Elasticity_ , 123(2):203–223, 2016\.
  * Guven et al. [2014] Jemal Guven, Dulce María Valencia, and Pablo Vázquez-Montejo.  Environmental bias and elastic curves on surfaces.  _Journal of Physics A: Mathematical and Theoretical_ , 47(35):355201, 2014.
  * Swigon et al. [1998] David Swigon, Bernard D Coleman, and Irwin Tobias.  The elastic rod model for dna and its application to the tertiary structure of dna minicircles in mononucleosomes.  _Biophysical journal_ , 74(5):2515–2530, 1998\.
  * Tobias et al. [2000] Irwin Tobias, David Swigon, and Bernard D Coleman.  Elastic stability of dna configurations. i. general theory.  _Physical Review E_ , 61(1):747, 2000.
  * Coleman et al. [2000] Bernard D Coleman, David Swigon, and Irwin Tobias.  Elastic stability of dna configurations. ii. supercoiled plasmids with self-contact.  _Physical Review E_ , 61(1):759, 2000.
  * Goriely and Neukirch [2006] Alain Goriely and Sébastien Neukirch.  Mechanics of climbing and attachment in twining plants.  _Physical review letters_ , 97(18):184302, 2006\.
  * Seemann [1996] W Seemann.  Deformation of an elastic helix in contact with a rigid cylinder.  _Archive of Applied Mechanics_ , 67(1):117–139, 1996.
  * Van der Heijden [2001] GHM Van der Heijden.  The static deformation of a twisted elastic rod constrained to lie on a cylinder.  _Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences_ , 457(2007):695–715, 2001.
  * Tan and Digby [1993] XC Tan and PJ Digby.  Buckling of drill string under the action of gravity and axial thrust.  _International journal of solids and structures_ , 30(19):2675–2691, 1993.
  * Wu et al. [1993] Jiang Wu, HC Juvkam-Wold, and R Lu.  Helical buckling of pipes in extended reach and horizontal wells-part 1: preventing helical buckling.  1993\.
  * Wu and Juvkam-Wold [1993] J Wu and HC Juvkam-Wold.  Helical buckling of pipes in extended reach and horizontal wells-part 2: Frictional drag analysis.  1993\.
  * Miller et al. [2014] JT Miller, Arnaud Lazarus, Basile Audoly, and Pedro M Reis.  Shapes of a suspended curly hair.  _Physical review letters_ , 112(6):068103, 2014\.
  * Mahadevan and Keller [1996] L Mahadevan and Joseph B Keller.  Coiling of flexible ropes.  _Proceedings of the royal society of london. Series A: mathematical, Physical and Engineering Sciences_ , 452(1950):1679–1694, 1996.
  * Mahadevan and Keller [1999] L Mahadevan and Joseph B Keller.  Periodic folding of thin sheets.  _Siam Review_ , 41(1):115–131, 1999.
  * Elettro et al. [2016] Hervé Elettro, Sébastien Neukirch, Fritz Vollrath, and Arnaud Antkowiak.  In-drop capillary spooling of spider capture thread inspires hybrid fibers with mixed solid–liquid mechanical properties.  _Proceedings of the National Academy of Sciences_ , 113(22):6143–6147, 2016.
  * Bico et al. [2004] José Bico, Benoit Roman, Loic Moulin, and Arezki Boudaoud.  Elastocapillary coalescence in wet hair.  _Nature_ , 432(7018):690–690, 2004.
  * Prasath et al. [2021] S Ganga Prasath, Joel Marthelot, Narayanan Menon, and Rama Govindarajan.  Wetting and wrapping of a floating droplet by a thin elastic filament.  _Soft Matter_ , 17(6):1497–1504, 2021.
  * [18] Fabian Brau, S Ganga Prasath, and Benny Davidovitch.  Wettability of bendable solids: Insights from a two-dimensional, inextensible model.  _(to be submitted)_.
  * Schulman et al. [2017] Rafael D Schulman, Amir Porat, Kathleen Charlesworth, Adam Fortais, Thomas Salez, Elie Raphaël, and Kari Dalnoki-Veress.  Elastocapillary bending of microfibers around liquid droplets.  _Soft matter_ , 13(4):720–724, 2017.
  * Py et al. [2007] Charlotte Py, Paul Reverdy, Lionel Doppler, José Bico, Benoit Roman, and Charles N Baroud.  Capillary origami: spontaneous wrapping of a droplet with an elastic sheet.  _Physical review letters_ , 98(15):156103, 2007\.
  * Audoly and Pomeau [2010] Basile Audoly and Yves Pomeau.  _Elasticity and geometry: from hair curls to the non-linear response of shells_.  Oxford university press, 2010.
  * Kreyszig [2019] Erwin Kreyszig.  Differential geometry.  In _Differential Geometry_. University of Toronto Press, 2019.
  * Kusumaatmaja and Lipowsky [2011] Halim Kusumaatmaja and Reinhard Lipowsky.  Droplet-induced budding transitions of membranes.  _Soft Matter_ , 7(15):6914–6919, 2011.
  * King et al. [2012] Hunter King, Robert D Schroll, Benny Davidovitch, and Narayanan Menon.  Elastic sheet on a liquid drop reveals wrinkling and crumpling as distinct symmetry-breaking instabilities.  _Proceedings of the National Academy of Sciences_ , 109(25):9716–9720, 2012.
  * Schroll et al. [2013] RD Schroll, M Adda-Bedia, E Cerda, J Huang, N Menon, TP Russell, KB Toga, D Vella, and B Davidovitch.  Capillary deformations of bendable films.  _Physical review letters_ , 111(1):014301, 2013\.
  * Davidovitch et al. [2011] Benny Davidovitch, Robert D Schroll, Dominic Vella, Mokhtar Adda-Bedia, and Enrique A Cerda.  Prototypical model for tensional wrinkling in thin sheets.  _Proceedings of the National Academy of Sciences_ , 108(45):18227–18232, 2011.

Supplemental Materials: Shapes of a filament on the surface of a bubble

##  S1 Coiling with two free ends

In our experiments we held one end of the filament fixed at the bubble’s north
pole. However when both the ends are free the calculation of the critical
length $\Phi_{c}$ for transition between different morphologies gets
simplified. The filament shape is given by either the geodesic or the polar
angle of the latitude based on the energy, which is a function of the coiling
parameter and elasto-gravity bendability. In order to calculate the minimum
energy, we use the expressions in Eq. 7, 8 (see Fig. S1$(a)$ for the
functional form when $\Omega^{-1}_{g}=1$). The minimum energy stays along the
geodesic up to a critical value of $\Phi\leq 5.13$ when $\Omega^{-1}_{g}=1$.
Beyond this critical value of coiling parameter the solution branch shifts to
that of a latitude (shown by a solid gray line), as the energy of the geodesic
exceeds that of the latitude of the same length. The solid line in Fig.
S1$(b)$ is the critical coiling parameter $\Phi_{c}$ when the filament
transitions from a geodesic to a different latitude calculated as a function
of $\Phi$ for different values of $\Omega^{-1}_{g}$. The filament starts
coiling only when the filament wraps a given latitude completely and comes
into self-contact. This is indicated by the dashed line in Fig. S1$(b)$, which
is calculated by using the relation $\Phi^{*}=2\pi\cos(\xi_{c})$, where
$\xi_{c}$ is the polar angle of the latitude configuration.

Figure S1: $(a)$ Total energy of a filament as a function of its non-
dimensional length, $\Phi$ when it lies on a longitude and a latitude for
$\Omega^{-1}_{g}=1$. Filaments choose the minimum of these and we find the
mode switch happen at $\Phi_{c}=5.13$. $(b)$ Phase boundary defined based on
the intersection of energies $\tilde{\mathscr{E}}_{l},\tilde{\mathscr{E}}_{g}$
in Eq. 7, 8. The filament switches modes beyond $\Phi_{c}$ to the minimum of
the energies, indicated here as boundary of the green region with a solid
line.

### Asymptotic behaviour

When both ends of the filament are free, we expect $\Phi_{c}\rightarrow 2\pi$
in the limit of zero gravity just as seen in the main text. Thus, when
$\Omega_{g}\rightarrow 0$ the selected latitude is described by polar angle
$\xi_{c}\rightarrow 0$. We have two small parameters $\xi_{c}$, $\Omega_{g}$
whose relationship is evaluated by expanding the solution close to
$\Phi\rightarrow 2\pi$ as $\Omega_{g}\rightarrow 0$. From the minimised
solution we know that the latitude’s polar angle satisfies,

$\Omega_{g}\cos^{4}\xi_{c}=-2\sin\xi_{c}.$

Expanding the above expression, we get the leading order behaviour:

$\xi_{c}\sim-\frac{\Omega_{g}}{2}.$

We find that $\xi_{c}$ is independent of $\Phi$ after the transition from the
geodesic i.e. $\Phi>\Phi_{c}$. Further as gravity becomes stronger the polar
angle to which filament migrates moves towards south pole. This comes from a
compromise between bending energy and potential energy where filament chooses
to bend more to reduce the potential energy with increase in $\Omega_{g}$. In
order to find the relationship between the critical filament length $\Phi_{c}$
at which bifurcation happens in this asymptotic limit, we expand the full
energy expression for a latitude in Eq. 7 to get,

$\displaystyle\tilde{\mathscr{E}}_{g}$
$\displaystyle=\Phi(\Omega_{g}(\sin\xi_{c}+z_{o})+1+\tan^{2}(\xi_{c})),$ (S1)
$\displaystyle\approx\Phi\bigg{(}1-\frac{\Omega_{g}^{2}}{4}\bigg{)}.$ (S2)

A similar expansion for the energy of a geodesic can be performed to get

$\displaystyle\tilde{\mathscr{E}}_{l}$
$\displaystyle=\Phi-2\Omega_{g}\sin\frac{\Phi}{2}.$ (S3)

Since we are interested in the region close to $\Phi\rightarrow 2\pi$, we have
the small parameter $\delta=(2\pi-\Phi)$. The geodesic energy can be written
in terms of $\delta$ as:

$\displaystyle\tilde{\mathscr{E}}_{l}$ $\displaystyle\approx
2\pi-\delta-\Omega_{g}\delta+\mathcal{O}(\delta^{3})$ (S4)

The criteria for transition is found by solving
$\tilde{\mathscr{E}}_{g}=\tilde{\mathscr{E}}_{l}$ and this gives,

$\displaystyle\frac{\Omega_{g}^{2}\delta}{4}-\frac{\pi\Omega_{g}^{2}}{2}$
$\displaystyle=-\Omega_{g}\delta,$ (S5) $\displaystyle\delta$
$\displaystyle\approx\frac{\pi\Omega_{g}}{2}.$ (S6)

We plot in Fig. S1$(c,d)$ the asymptotic expressions derived above for
$\xi_{c}$ and $\delta$ as dashed lines and compare it with the full solution
shown as solid curves.

##  S2 Regularizing singularities

In Eq. 9 we have singularities at the south pole of the bubble because the
curvature of the fraction $(1-\alpha)$ of the filament along the latitude
diverges as the square of curvature of the latitude portion, whereas the
fraction of length along latitude $(1-\alpha)\rightarrow 0$ contributes only
linearly. In order to avoid this divergence, we need a regularization length
scale arising from the missing length-scale in the model, the bendo-capillary
length $l_{\rm bc}$. To resolve this divergence we multiply the singular
contribution with a function $\Lambda(\beta)$ that suppresses the singularity.
The energy then becomes,

$\displaystyle\tilde{\mathscr{E}}=$ $\displaystyle\
\alpha\Phi+\Omega_{g}\sin(\alpha\Phi)+\Omega_{g}\Phi(1-\alpha)\cos(\alpha\Phi)$
$\displaystyle\
+(1-\alpha)\Phi\big{\\{}[1-\Lambda(0)-\Lambda(\pi)-\Lambda(2\pi)]\csc^{2}(\alpha\Phi)\big{\\}},$
(S7) $\text{where}\ \Lambda(\beta)=\
\exp\bigg{(}-\frac{(\alpha\Phi-\beta)^{2}}{\eta(\Omega)}\bigg{)}.$

Here $\eta(\Omega)$ is the cut-off non-dimensional length-scale associated
with the capillary bendability, as this is the length scale that acts at the
boundary of scale separation. However owing to the limits of scale-separation
when the radius of curvature along a latitude approaches the length-scale
associated with capillary bendability, this energy expression is no more
valid. These regularizing terms are relevant at three locations along the
filament i.e. $\alpha\Phi=0,\pi,2\pi$, as all these locations result in
diverging values of the energy.


# New multi-hump exact solitons of a coupled Korteweg-de-Vries system with
conformable derivative describing shallow water waves via RCAM

Prakash Kumar Das  Department of Mathematics, Trivenidevi Bhalotia College,
Raniganj, Paschim Bardhaman, West Bengal, India-713347
<EMAIL_ADDRESS>
###### Abstract

In this article, a modification of the rapidly convergent approximation method
is proposed to solve a coupled Korteweg-de Vries equations with conformable
derivative that govern shallow-water waves. Based on the Leibniz and chain
rule of conformable derivative, these equations reduced into ODEs with
integer-order using traveling wave transformation. Adopting the modified
scheme a new novel exact solution of the reduced coupled ordinary differential
equations is obtained in terms of exponential functions. Finally, by putting
them back into traveling wave transformation the solutions of the considered
partial differential equations with conformable derivative are derived. To
ensure the boundedness of the derived solutions few theorems have been
proposed and proved. The derived results of the theorems are utilized to plot
the solutions. Graphics exhibit that solutions have variant multi-hump soliton
peculiarities and their tails decay to zero exponentially in a monotonic
manner. These results not only show the efficiency of the modified scheme but
also establish that the solution is enriched with new multi-hump features.

  * August 2020

Keywords: A coupled KdV equations, Conformable derivative, Rapidly convergent
approximation method, Exact solutions, Boundedness, Multi-hump solitons

##  1 Introduction

One-pulse solutions are common and significant features of dispersive partial
differential equations recounting physical systems. Beyond that, there may
exist two-pulse, three-pulse, and generally n-pulse (multi-hump) solutions for
such physical systems having a higher-order dispersion ( involving higher-
order differential operators). Moreover, when a spatial soliton is built of
multiple modes of the initiated waveguide, its soliton intensity profile may
be compounded and display several peaks. The multi-pulse solutions constitute
with multiple copies of the one-pulse solutions separated by finitely many
oscillations close to the zero equilibrium. These solutions have been derived
numerically for various nonlinear models as discussed below.

Numerically it is showed that there may exist infinitely many multi-pulse
solutions to the fifth order Korteweg -de Vries (KdV) equation [1, 2]. Also,
the existence of an infinite (but countable) set of localized stationary
solutions, experimentally and numerically established for nonlinear
generalised KdV equation in [3]. For solitary waves governed by incoherent
beam interaction in a saturable medium Ostrovskaya et al., [4] disclosed that
two-hump solitary waves are linearly stable in a wide region of their
existence, but all three-hump solitons are linearly unstable. For
nonintegrable multi-component nonlinear models formation of stable multi-hump
solitary waves, a novel physical mechanism has been presented in [5]. These
multi-hump optical solitons observed numerically in the models describing
laser radiation copropagating with a Bose-Einstein condensate and used to shed
light on the phenomenon of jet emission from a condensate interacting with a
laser [6].

Besides numerical methods, few direct methods also exist in literature to
investigate multi-hump solitons. In [7] the Hirota bilinear method has been
employed to find one- and two-hump exact bright and dark soliton solutions to
a coupled system between the linear Schr$\ddot{\text{o}}$dinger and KdV
equations. Using the Lax pair and Darboux transformation (DT), the authors in
[8] have derived multi-peak soliton. Most of the multi-pick solitons studied
in literature are using numerical schemes. So it is of great importance to
propose a new method which can derive exact multi-hump soliton solutions
efficiently and a simple way than existing direct methods. To attain the goal
here for the first time the rapidly convergent approximation method (RCAM) has
been employed to construct exact multi-hump solutions of a coupled
differential equations with conformable derivative.

Non-linear frational differential equations (NLFDE) model abundant branches of
applied mathematics describing evolution of physical processes in science and
engineering [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]. The
availability of solutions of such equations contribute a lot in visualize the
fractional system involved. In lots of cases, it is arduous task to derive
exact solution to these equations. Despite that reserchers have proposed a few
direct methods [24, 25, 26, 27, 28, 29, 30, 31, 32, 33] to look for exact
travelling wave solutions of NLFDE. The above mention direct methods often
make few assumptions of the form of the solution, as a result, they always
unable to produce a profound type solution. In addition to that, these methods
methodology always require to solve a system of nonlinear algebraic equations,
which is a hard task and becomes improbable when nonlinearity increases.
Whenever these systems of equations are solvable, they lead to some particular
solutions. As a consequence, these methods yield a set of special solutions
instead of a general one. Also, there are many approximation methods [34, 35,
36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46] to deal with equations for which
the direct methods do not work. These approximate methods are often found to
be slowly convergent and unable to provide the close form of the series
solution. These problems can be easily tackled by the RCAM [47, 48, 49, 50,
51, 52, 53, 54]. In this article, this scheme is used to obtain a new multi-
hump travelling wave solution of a coupled Korteweg-de Vries equations with
conformable derivative.

The famous couple KdV equation is a nonlinear frequency dispersion equation,
which describes shallow long-wave and small-amplitude phenomena, ion acoustic
waves in a plasma, acoustic waves on a crystal lattice and so on. For the time
being its generalisation containing time-fractional derivative getting great
deal attention to the researchers whose form is:

$\displaystyle D_{t}^{\alpha}u(x,t)-p\ u_{xxx}(x,t)-p_{1}\ u(x,t)\
u_{x}(x,t)-p_{2}\ v(x,t)\ v_{x}(x,t)=f(x,t),$ $\displaystyle
D_{t}^{\alpha}v(x,t)+q\ v_{xxx}(x,t)+q_{1}\ u(x,t)\ v_{x}(x,t)=g(x,t),\
0<\alpha\leq 1,$ (1)

where $D_{*}^{\alpha}$ denote the Caputo type differential operator and
$f(x,t)$ and $g(x,t)$ are two source terms respectively. Eq. (1) play an
important role in the propagation of waves and have significant contribution
in many applied fields such as fluid, mechanics, plasmas, crystal lattice
vibrations at low temperatures, etc. The Eq. (1) without source term
($f(x,t)=g(x,t)=0$) and $p=a,\ p_{1}=6a,\ p_{2}=2b,\ q=-1,\ q_{1}=3$ was first
proposed and obtained its approximate series solution by using Adomian
decomposition method in [39]. Latter this equation was considered by many
authors [55, 40, 56, 57] and they have applied the generalised differential
transform method, homotopy decomposition method, fractional reduced
differential transform method and Haar wavelets respectively to obtain the
approximate solutions. Matinfar et al. [58] applied the functional variable
method for deriving the exact solution of the said equations. Three numerical
technique based on the shifted Legendre polynomials, Meshless spectral method
and spectral collection method were applied in [59, 60, 61] for obtaining the
numerical solution of the Eq. (1) with source terms. In a recent work S.
Biswas et al., [41] proposed another variant of coupled fractional KdV
equation containing both space-time fractional derivative. They derived it
with the help of a semi-inversion method, variational principle, and
Lagrangian of the KdV equation. In addition, the authors employed the homotopy
analysis method to derive the approximate solution of the equation. The said
equation enjoy the form

$\displaystyle D_{t}^{\alpha}u(x,t)+6\ a\ u(x,t)\ D_{x}^{\alpha}u(x,t)-2\ b\
v(x,t)\ D_{x}^{\alpha}v(x,t)+a\ D_{x}^{3\alpha}u(x,t)=0,$ $\displaystyle
D_{t}^{\alpha}v(x,t)+3\ u(x,t)\ D_{x}^{\alpha}v(x,t)+D_{x}^{3\alpha}v(x,t)=0,\
\ 2<3\alpha<3.$ (2)

Equations (1) and (1) are defined using Caputo derivative, which do not
satisfy some main principles of classical integer order derivative [62, 63]
such as Leibniz rule, chain rule and etc. So it is not straightforward to
derive exact solutions to these equations involving this derivative. Hence in
this article, we consider the equation (1) with conformable derivative in the
form

$\displaystyle u^{(\beta)}_{t}(x,t)+6\ a\ u(x,t)\ u^{(\alpha)}_{x}(x,t)-2\ b\
v(x,t)\ v^{(\alpha)}_{x}(x,t)+a\ u^{(\alpha)}_{xxx}(x,t)=0,$ $\displaystyle
v^{(\beta)}_{t}(x,t)+3\ u(x,t)\
v^{(\alpha)}_{x}(x,t)+v^{(\alpha)}_{xxx}(x,t)=0,\ \ 0<\alpha,\ \beta\leq 1.$
(3)

Since the system (1) contains conformal derivatives is equivalent to the
classical system of integer-order KdV-type equations with variable
coefficients of a special kind [64]. This coupled KdV equations with variable
coefficients represents a simple generalization of Hirota-Satsuma coupled KdV
equations [65, 66], describe the interaction of two long waves [66, 67, 68]
and also has many applications in the above mentioned fields. Also in the case
$\alpha=\beta=1$, the above equations reduces to conventional Hirota-Satsuma
coupled KdV equations. In this assignment, the author’s intention is to
construct the exact multi-hump solutions of Eq.(1) by utilizing RCAM.

Section 2 presents the basic properties of conformable derivative. Basic
methodologies of RCAM are introduced in section 3. Using this scheme a class
of new travelling wave solutions for coupled KdV equation with conformable
derivative were presented in section 4. The boundedness conditions of
solutions of the KdV equation with conformable derivative and its reduced
ordinary differential equations have been presented in section 5 and verified
through plotting them. The outlook of the present work have been summarized
and some concluding remarks given in section 6.

##  2 Properties of conformable derivative

This section presents few definitions and features of conformable derivative
[27, 69]:

######  Definition 2.1.

The conformable derivative of a function $f:[0,\infty)\rightarrow\mathbb{R}$
of order $\alpha$ is defined by

$\displaystyle T_{\alpha}(h)(t)=\lim_{\epsilon\rightarrow
0}\frac{h\left(t+\epsilon\ t^{1-\alpha}\right)-h(t)}{\epsilon},$ (4)

for all $t>0,\ \alpha\in(0,1).$ If f is $\alpha$-differentiable in some
$(0,a),a>0,$ and $\lim_{t\rightarrow 0+}f^{(\alpha)}(t)$ exists, then define
$f^{(\alpha)}(0)=\lim_{t\rightarrow 0+}f^{(\alpha)}(t).$

Sometimes, we use the notation $f^{(\alpha)}(t)$ in place of
$T_{\alpha}(f)(t)$, to prevail the conformable derivatives of $f$ of order
$\alpha$. Few significant features of conformable derivative are presented
below:
If $\alpha\in(0,1],$ and $f,g$ be $\alpha$-differentiable at a point $t>0$
then we have

  * 1.

$T_{\alpha}\left(a\ f+b\ g\right)=a\ T_{\alpha}(f)+b\ T_{\alpha}(g),$ for all
$a,b\in\mathbb{R}.$

  * 2.

$T_{\alpha}\left(t^{p}\ \right)=p\ t^{p-\alpha}$ for all $p\in\mathbb{R}.$

  * 3.

$T_{\alpha}(\lambda)=0$, for all constant functions $f(t)=\lambda.$

  * 4.

$T_{\alpha}\left(fg\right)=f\ T_{\alpha}(g)+g\ T_{\alpha}(f).$

  * 5.

$T_{\alpha}\left(\frac{f}{g}\right)=\frac{g\ T_{\alpha}(f)+f\
T_{\alpha}(g)}{g^{2}}$.

  * 6.

If $f$ is differentiable, then
$T_{\alpha}(f)(t)=t^{1-\alpha}\frac{df}{dt}(t).$

######  Theorem 2.2 (Chain Rule [69, 70, 71]).

Let $f,g:(0,\infty)\rightarrow\mathbb{R}$ be two differentiable functions and
also $f$ is $\alpha$-differentiable, then, one has the following rule:
$T_{\alpha}(fog)(t)=t^{1-\alpha}\ g^{\prime}(t)f^{\prime}(g(t)).$

##  3 The modified RCAM

Consider a system of equations in the form

$\displaystyle{\cal X}^{{}^{\prime\prime\prime}}-{\cal A}^{2}.\;{\cal
X}^{{}^{\prime}}={\cal N},$ (5)

where ${\cal X},{\cal A}$ and ${\cal N}$ are matrix of dependent variables,
constant coefficients and nonlinear terms respectively, given by

${\cal X}=\left[\begin{array}[]{c}U_{1}(x)\\\ U_{2}(x)\\\ \vdots\\\
U_{k}(x)\\\ \end{array}\right],\ \ {\cal
A}=\left[\begin{array}[]{c}\lambda_{1}\ \ 0\ \ \cdots 0\\\ 0\ \ \lambda_{2}\ \
\cdots 0\\\ \vdots\\\ 0\ \ 0\ \ \cdots\ \ \lambda_{k}\\\ \end{array}\right],\
\text{and}\ {\cal
N}=\left[\begin{array}[]{c}N_{1}\left(U_{1}(x),\cdots,U_{k}(x)\right)\\\
N_{2}\left(U_{1}(x),\cdots,U_{k}(x)\right)\\\ \vdots\\\
N_{k}\left(U_{1}(x),\cdots,U_{k}(x)\right)\\\ \end{array}\right].$

Note that here considered all $\lambda_{i},\ i=1,2,\cdots,k$ are distinct.
Which indicate the modification over the existing scheme [53]. To find the
solution of (5), we recast it in exponential matrix operator form

$\hat{{\cal O}}[{\cal X}](x)={\cal N},$ (6)

where linear exponential matrix operator can be recast in the form

$\hat{{\cal O}}[\cdot](x)=e^{{\cal A}.x}\frac{d}{dx}\left(e^{-2{\cal
A}.x}\frac{d}{dx}\left(e^{{\cal A}.x}\frac{d}{dx}[\cdot]\right)\right).$ (7)

The inverse operator $\hat{{\cal O}}^{-1}$ of the operator $\hat{{\cal
O}}[](x)$ is given by

$\displaystyle\hat{{\cal O}}^{-1}[\cdot](x)=\int^{x}e^{-{\cal
A}.x^{{}^{\prime}}}\int^{x^{{}^{\prime}}}e^{2{\cal
A}.x^{{}^{\prime\prime}}}\int^{x^{{}^{\prime\prime}}}e^{-{\cal
A}.x^{{}^{\prime\prime\prime}}}[\cdot]dx^{{}^{\prime\prime\prime}}dx^{{}^{\prime\prime}}dx^{{}^{\prime}}.$
(8)

Operating $\hat{{\cal O}}^{-1}$ on $\hat{{\cal O}}[{\cal X}](x)$ and using
integration by parts produces

$\displaystyle\hat{{\cal O}}^{-1}\left({\cal
X}^{{}^{\prime\prime\prime}}-{\cal A}^{2}.\;{\cal
X}^{{}^{\prime}}\right)={\cal X}-{\cal C}.e^{{\cal A}.x}-{\cal D}.e^{-{\cal
A}.x}-{\cal E},$ (9)

where ${\cal C}=\left[\begin{array}[]{c}c_{1}\\\ c_{2}\\\ \vdots\\\ c_{k}\\\
\end{array}\right]$, ${\cal D}=\left[\begin{array}[]{c}d_{1}\\\ d_{2}\\\
\vdots\\\ d_{k}\\\ \end{array}\right]$ and ${\cal
E}=\left[\begin{array}[]{c}e_{1}\\\ e_{2}\\\ \vdots\\\ e_{k}\\\
\end{array}\right]$ are integration constant matrices. Operating $\hat{{\cal
O}}^{-1}$ both sides of (6) and using (9), gives

${\cal X}={\cal E}+{\cal C}.e^{{\cal A}.x}+{\cal D}.e^{-{\cal A}.x}+\hat{\cal
O}^{-1}[{\cal N}](x),$ (10)

involving three arbitrary constants matrices ${\cal C}$, ${\cal D}$ and ${\cal
E}$. To derive the terms involving the unknown ${\cal X}$ in R.H.S of (10), we
express them in the form

$\displaystyle{\cal X}$ $\displaystyle\cong$
$\displaystyle\left[\begin{array}[]{c}U_{1}(x)\\\ U_{2}(x)\\\ \vdots\\\
U_{k}(x)\\\
\end{array}\right]=\sum_{m=0}^{\infty}\left[\begin{array}[]{c}U_{1,m}(x)\\\
U_{2,m}(x)\\\ \vdots\\\ U_{k,m}(x)\\\ \end{array}\right],$ (19)

and ${\cal N}=\sum_{m=0}^{\infty}\Delta_{m}(x),$ with

$\displaystyle\Delta_{m}\cong\left[\begin{array}[]{c}\Delta_{1,m}(x)\\\
\Delta_{2,m}(x)\\\ \vdots\\\ \Delta_{k,m}(x)\\\
\end{array}\right]=\frac{1}{m!}\left[\frac{d^{m}}{d\epsilon^{m}}\left[\begin{array}[]{c}N_{1}(x)\\\
N_{2}(x)\\\ \vdots\\\ N_{k}(x)\\\ \end{array}\right]\right]_{\epsilon=0},$
(28)

and $N_{j}(x),\ j=1,2,\cdots,k$ are given by

$N_{j}(x)=\left(\sum_{k=0}^{\infty}U_{1,k}\epsilon^{k},\sum_{k=0}^{\infty}U_{2,k}\epsilon^{k},\cdots,\sum_{k=0}^{\infty}U_{k,k}\epsilon^{k}\right).$

The terms
$\Delta_{i,m}(x)=\Delta_{i,m}(U_{1,0}(x),U_{1,1}(x),…..,U_{1,m}(x),\cdots,U_{k,0}(x),U_{k,1}(x)$
$,…..,U_{k,m}(x)),\ i=1,2,\cdots,k$ are Adomian polynomials [72, 73, 74, 75,
76, 77] turned out from the formula (28). Use of (28) in (47) gives

$\displaystyle{\cal
X}=\left[\begin{array}[]{c}c_{1}\;e^{\lambda_{1}x}+d_{1}\;e^{-\lambda_{1}x}+e_{1}\\\
c_{2}\;e^{\lambda_{2}x}+d_{2}\;e^{-\lambda_{2}x}+e_{2}\\\ \vdots\\\
c_{k}\;e^{\lambda_{k}x}+d_{k}\;e^{-\lambda_{k}x}+e_{k}\\\
\end{array}\right]+\hat{\cal
O}^{-1}\left[\sum_{m=0}^{\infty}\left[\begin{array}[]{c}\Delta_{1,m}(x)\\\
\Delta_{2,m}(x)\\\ \vdots\\\ \Delta_{k,m}(x)\\\ \end{array}\right]\right].$
(37)

We follow the steps given in [53], to obtain the higher order correction term
as

$\displaystyle{\cal X}_{0}\cong\left[\begin{array}[]{c}U_{1,0}(x)\\\
U_{2,0}(x)\\\ \vdots\\\ U_{k,0}(x)\\\
\end{array}\right]=\left[\begin{array}[]{c}c_{1}\;e^{\lambda_{1}x}+d_{1}\;e^{-\lambda_{1}x}+e_{1}\\\
c_{2}\;e^{\lambda_{2}x}+d_{2}\;e^{-\lambda_{2}x}+e_{2}\\\ \vdots\\\
c_{k}\;e^{\lambda_{k}x}+d_{k}\;e^{-\lambda_{k}x}+e_{k}\\\ \end{array}\right],$
(47) $\displaystyle{\cal X}_{n+1}\cong\left[\begin{array}[]{c}U_{1,n+1}(x)\\\
U_{2,n+1}(x)\\\ \vdots\\\ U_{k,n+1}(x)\\\ \end{array}\right]=\hat{\cal
O}^{-1}\left[\begin{array}[]{c}\Delta_{1,n}(x)\\\ \Delta_{2,n}(x)\\\ \vdots\\\
\Delta_{k,n}(x)\\\ \end{array}\right],\ n\geq 0.$ (56)

In the following, we present a few special cases of iterative formulas
(47)-(56).
Case-I In case $\lambda_{i}>0,\ i=1,\cdots,m,$ and $\lambda_{i}<0,\
i=m+1,\cdots,k$ for $0<m<k$, use of the vanishing boundary condition
$U_{i}(-\infty)=0$ in (69) for the localized solution impart $d_{i}=e_{i}=0,\
\text{for}\ i=1,2,\cdots,m$ and $c_{i}=e_{i}=0,\ \text{for}\
i=m+1,m+2,\cdots,k$. Consequently the leading and higher order rectification
terms of the series solution are given by (56) with

$\displaystyle{\cal X}_{0}\cong\left[\begin{array}[]{c}U_{1,0}(x)\\\ \vdots\\\
U_{m,0}(x)\\\ U_{m+1,0}(x)\\\ \vdots\\\ U_{k,0}(x)\\\
\end{array}\right]=\left[\begin{array}[]{c}c_{1}\;e^{\lambda_{1}x}\\\
\vdots\\\ c_{m}\;e^{\lambda_{m}x}\\\ d_{m+1}\;e^{-\lambda_{m+1}x}\\\ \vdots\\\
d_{k}\;e^{-\lambda_{k}x}\\\ \end{array}\right].$ (69)

Case-II In case all $\lambda_{i}>0$, use of the condition $U_{i}(-\infty)=0$
in (69) for the localized solution confers $d_{i}=e_{i}=0,\ i=1,2,\cdots,k$.
Hence the leading and higher order amendment terms of the series solution are
given by (56) with

$\displaystyle{\cal X}_{0}\cong\left[\begin{array}[]{c}U_{1,0}(x)\\\
U_{2,0}(x)\\\ \vdots\\\ U_{k,0}(x)\\\
\end{array}\right]=\left[\begin{array}[]{c}c_{1}\;e^{\lambda_{1}x}\\\
c_{2}\;e^{\lambda_{2}x}\\\ \vdots\\\ c_{k}\;e^{\lambda_{k}x}\\\
\end{array}\right].$ (78)

Case-III Likewise in case when all $\lambda_{i}<0$, for boundary condition
$U_{i}(-\infty)=0$ to obtain the localized solution we go along with
restraining the term involving $e^{\lambda_{i}x}$. In this case, the
commanding and subsequent terms of the solution are given by (56) with

$\displaystyle{\cal X}_{0}\cong\left[\begin{array}[]{c}U_{1,0}(x)\\\
U_{2,0}(x)\\\ \vdots\\\ U_{k,0}(x)\\\
\end{array}\right]=\left[\begin{array}[]{c}d_{1}\;e^{-\lambda_{1}x}\\\
d_{2}\;e^{-\lambda_{2}x}\\\ \vdots\\\ d_{k}\;e^{-\lambda_{k}x}\\\
\end{array}\right].$ (87)

Similarly for boundary condition $U_{i}(\infty)=0$, one needs to proceed
contrary to the above cases. Using above presented iteration schemes and
symbolic software, one can obtain the general term of the series (or
generating function). Gives the exact solution of the considered system of
equations.

##  4 Solution of coupled KdV Eq. (1)

Let us consider the solution of the equation (1) by using travelling wave
transformation [25, 70, 71] in the form

$\displaystyle\begin{cases}u(x,t)=U(\xi),\\\
v(x,t)=V(\xi),\end{cases}\xi=\frac{k}{\alpha}\ x^{\alpha}+\frac{c}{\beta}\
t^{\beta}+\xi_{0},$ (88)

where $c,\ k$ and $\xi_{0}$ are constants. Using the chain rule of conformable
derivatives, the transformation (88) permit us to convert (1) to an ordinary
differential equation in the form :

$\displaystyle\begin{cases}U^{{}^{\prime\prime\prime}}(\xi)-\lambda_{1}^{2}\
U^{{}^{\prime}}(\xi)+\frac{6}{k^{2}}U(\xi)U^{{}^{\prime}}(\xi)-\frac{2b}{ak^{2}}V(\xi)V^{{}^{\prime}}(\xi)=0,\\\
V^{{}^{\prime\prime\prime}}(\xi)-\lambda_{2}^{2}\
V^{{}^{\prime}}(\xi)+\frac{3}{k^{2}}U(\xi)V^{{}^{\prime}}(\xi)=0,\end{cases}$
(89)

where

$\displaystyle\lambda_{1}=\sqrt{-\frac{c}{ak^{3}}}\ \text{ and}\
\lambda_{2}=\sqrt{-\frac{c}{k^{3}}}.$ (90)

It is important to note here that throughout our discussion, we assume that
$\lambda_{1}$ and $\lambda_{2}$ are positive real and that leads us to the
conditions $a>0$ and $ck^{3}<0$. The above presented system of equations can
be recast in the form (5) with

${\cal X}=\left[\begin{array}[]{c}U(\xi)\\\ V(\xi)\\\ \end{array}\right],\ \
{\cal A}=\left[\begin{array}[]{c}\lambda_{1}\ \ 0\\\ 0\ \ \lambda_{2}\\\
\end{array}\right],\ \text{and}\ \ {\cal
N}=\left[\begin{array}[]{c}-\frac{6}{k^{2}}U(\xi)U^{\prime}(\xi)+\frac{2b}{ak^{2}}V(\xi)V^{\prime}(\xi)\\\
-\frac{3}{k^{2}}U(\xi)V^{\prime}(\xi)\\\ \end{array}\right].$

Thereafter, we pursue the imitating steps of RCAM presented in the section 3
to construct the solution of (89). To constitute localized solutions
satisfying boundary condition $U(-\infty)=V(-\infty)=0$, we use (56) with
(78). That yield the following correction terms

$\displaystyle U_{0}(\xi)=$ $\displaystyle c_{1}\ e^{\lambda_{1}\xi},$
$\displaystyle V_{0}(\xi)=$ $\displaystyle c_{2}\ e^{\lambda_{2}\xi},$
$\displaystyle U_{1}(\xi)=$
$\displaystyle-\left[e^{2\lambda_{2}\xi}\left(ac_{1}^{2}\left(\lambda_{1}^{2}-4\lambda_{2}^{2}\right)e^{2\left(\lambda_{1}-\lambda_{2}\right)\xi}+bc_{2}^{2}\lambda_{1}^{2}\right)\right]\text{/}\left[ak^{2}\left(\lambda_{1}^{4}-4\lambda_{1}^{2}\lambda_{2}^{2}\right)\right],$
$\displaystyle V_{1}(\xi)=$
$\displaystyle-\left[3c_{1}c_{2}\lambda_{2}e^{\left(\lambda_{1}+\lambda_{2}\right)\xi}\right]\text{/}\left[k^{2}\lambda_{1}\left(\lambda_{1}+\lambda_{2}\right)\left(\lambda_{1}+2\lambda_{2}\right)\right],$
$\displaystyle U_{2}(\xi)=$
$\displaystyle\left[3c_{1}e^{\lambda_{1}\xi}\left(ac_{1}^{2}\lambda_{2}\left(\lambda_{1}+\lambda_{2}\right)^{2}\left(\lambda_{1}^{2}-4\lambda_{2}^{2}\right)e^{2\lambda_{1}\xi}+2bc_{2}^{2}\left(\lambda_{1}^{2}+2\lambda_{2}^{2}\right)\lambda_{1}^{3}e^{2\lambda_{2}\xi}\right)\right]$
$\displaystyle\text{/}\left[4ak^{4}\lambda_{1}^{4}\lambda_{2}\left(\lambda_{1}+\lambda_{2}\right)^{2}\left(\lambda_{1}^{2}-4\lambda_{2}^{2}\right)\right],$
$\displaystyle V_{2}(\xi)=$
$\displaystyle\left[c_{2}e^{\lambda_{2}\xi}\left(12ac_{1}^{2}\left(\lambda_{1}-2\lambda_{2}\right)\lambda_{2}^{3}e^{2\lambda_{1}\xi}+bc_{2}^{2}\left(\lambda_{1}+\lambda_{2}\right)\lambda_{1}^{3}e^{2\lambda_{2}\xi}\right)\right]\text{/}\left[8ak^{4}\lambda_{1}^{3}\lambda_{2}^{2}\right.$
$\displaystyle\left.\times\left(\lambda_{1}+\lambda_{2}\right)\left(\lambda_{1}^{2}-4\lambda_{2}^{2}\right)\right],$
$\displaystyle U_{3}(\xi)=$
$\displaystyle-\left[2a^{2}c_{1}^{4}\left(\lambda_{1}-2\lambda_{2}\right)^{2}\lambda_{2}^{2}\left(\lambda_{1}+2\lambda_{2}\right)^{3}e^{4\lambda_{1}\xi}+12abc_{1}^{2}c_{2}^{2}\lambda_{2}\left(\lambda_{1}^{3}-2\lambda_{2}\lambda_{1}^{2}\right.\right.$
$\displaystyle\left.+2\lambda_{2}^{2}\lambda_{1}-4\lambda_{2}^{3}\right)\lambda_{1}^{3}e^{2\left(\lambda_{1}+\lambda_{2}\right)\xi}+\left.b^{2}c_{2}^{4}\left(\lambda_{1}+2\lambda_{2}\right)\lambda_{1}^{6}e^{4\lambda_{2}\xi}\right]$
$\displaystyle\text{/}\left[4a^{2}k^{6}\lambda_{1}^{6}\left(\lambda_{1}-2\lambda_{2}\right)^{2}\lambda_{2}^{2}\left(\lambda_{1}+2\lambda_{2}\right)^{3}\right],$
$\displaystyle V_{3}(\xi)=$
$\displaystyle-\left[3c_{1}c_{2}e^{\left(\lambda_{1}+\lambda_{2}\right)\xi}\left(2ac_{1}^{2}\lambda_{2}^{2}\left(\lambda_{1}+2\lambda_{2}\right)^{2}\left(\lambda_{1}^{2}-\lambda_{1}\lambda_{2}-2\lambda_{2}^{2}\right)e^{2\lambda_{1}\xi}\right.\right.$
$\displaystyle\left.\left.+3bc_{2}^{2}\left(\lambda_{1}^{2}+\lambda_{2}\lambda_{1}+2\lambda_{2}^{2}\right)\lambda_{1}^{4}e^{2\lambda_{2}\xi}\right)\right]\text{/}\left[8ak^{6}\lambda_{1}^{5}\left(\lambda_{1}-2\lambda_{2}\right)\lambda_{2}\left(\lambda_{1}+\lambda_{2}\right)^{2}\right.$
$\displaystyle\times\left.\left(\lambda_{1}+2\lambda_{2}\right)^{3}\right],$
$\displaystyle\ \ \ \ \vdots$

where $c_{1}$ and $c_{2}$ are integration constants. Likewise, the higher
order terms can be calculated using symbolic computations available in
software packages Mathematica and Maple. Afterward calculating up to eight or
higher order iterations, the said software packages can easily provide the
generating functions of those correction terms. For this problem, we have
obtained the following generating functions

$\displaystyle U(\xi,\epsilon)=$ $\displaystyle
4k^{4}\left(\lambda_{1}-2\lambda_{2}\right)\left(\lambda_{1}+\lambda_{2}\right)^{2}\left(\lambda_{1}+2\lambda_{2}\right)^{2}\left[64c_{1}a^{2}k^{8}\lambda_{1}^{4}\lambda_{2}^{4}\left(\lambda_{1}-2\lambda_{2}\right)\right.$
$\displaystyle\left.\times\left(\lambda_{1}+\lambda_{2}\right)^{2}\left(\lambda_{1}+2\lambda_{2}\right)^{4}e^{\lambda_{1}\xi}-16\epsilon
abk^{2}c_{2}^{2}\lambda_{2}^{2}\left(\lambda_{1}+2\lambda_{2}\right)\left\\{\epsilon
c_{1}k^{2}\lambda_{1}^{2}\right.\right.$
$\displaystyle\times\left.\left.\left(\lambda_{1}^{2}-4\lambda_{2}^{2}\right)^{2}\left(\lambda_{1}^{2}+\lambda_{2}^{2}\right)e^{\left(\lambda_{1}+2\lambda_{2}\right)\xi}+\epsilon^{2}c_{1}^{2}\lambda_{2}^{2}\left(\lambda_{1}^{2}-3\lambda_{1}\lambda_{2}+2\lambda_{2}^{2}\right)^{2}\right.\right.$
$\displaystyle\times\left.\left.e^{2\left(\lambda_{1}+\lambda_{2}\right)\xi}+4k^{4}\lambda_{1}^{4}\lambda_{2}^{2}\left(\lambda_{1}^{2}+3\lambda_{2}\lambda_{1}+2\lambda_{2}^{2}\right)^{2}e^{2\lambda_{2}\xi}\right\\}+\epsilon^{4}b^{2}c_{1}c_{2}^{4}\lambda_{1}^{4}\right.$
$\displaystyle\times\left.\left(\lambda_{1}-2\lambda_{2}\right)\left(\lambda_{1}-\lambda_{2}\right)^{2}e^{\left(\lambda_{1}+4\lambda_{2}\right)\xi}\right]\text{/}Q(\xi,\epsilon)^{2},$
(91) $\displaystyle V(\xi,\epsilon)=$ $\displaystyle
8ac_{2}k^{4}\lambda_{2}^{2}\left(\lambda_{1}-2\lambda_{2}\right)\left(\lambda_{1}+\lambda_{2}\right)\left(\lambda_{1}+2\lambda_{2}\right)^{2}e^{\lambda_{2}\xi}\left[\epsilon
c_{1}\left(\lambda_{1}^{2}-3\lambda_{2}\lambda_{1}\right.\right.$
$\displaystyle\left.\left.+2\lambda_{2}^{2}\right)e^{\lambda_{1}\xi}+2k^{2}\lambda_{1}^{2}\left(\lambda_{1}^{2}+3\lambda_{2}\lambda_{1}+2\lambda_{2}^{2}\right)\right]\text{/}Q(\xi,\epsilon),$
(92)

where

$\displaystyle Q(\xi,\epsilon)=$ $\displaystyle
8ak^{4}\lambda_{2}^{2}\left(\lambda_{1}-2\lambda_{2}\right)\left(\lambda_{1}+\lambda_{2}\right)^{2}\left(\lambda_{1}+2\lambda_{2}\right)^{3}\left(\epsilon
c_{1}e^{\lambda_{1}\xi}+2k^{2}\lambda_{1}^{2}\right)-b\epsilon^{2}c_{2}^{2}e^{2\lambda_{2}\xi}$
$\displaystyle\times\left\\{\epsilon
c_{1}\left(\lambda_{1}^{2}-3\lambda_{2}\lambda_{1}+2\lambda_{2}^{2}\right)^{2}e^{\lambda_{1}\xi}+2k^{2}\lambda_{1}^{2}\left(\lambda_{1}^{2}+3\lambda_{2}\lambda_{1}+2\lambda_{2}^{2}\right)^{2}\right\\}.$
(93)

One can easily check that series expansion of generating functions
$U(\xi,\epsilon),V(\xi,\epsilon)$ in $\epsilon$ about $0$ provides the
correction terms as coefficients of several powers of $\epsilon$.
Consequently, the exact solution of (89) can be obtained by setting
$\epsilon=1$ in (4)-(4) in the form

$\displaystyle U(\xi)=$ $\displaystyle
4k^{4}\left(\lambda_{1}-2\lambda_{2}\right)\left(\lambda_{1}+\lambda_{2}\right)^{2}\left(\lambda_{1}+2\lambda_{2}\right)^{2}\left[64c_{1}a^{2}k^{8}\lambda_{1}^{4}\lambda_{2}^{4}\left(\lambda_{1}-2\lambda_{2}\right)\right.$
$\displaystyle\left.\times\left(\lambda_{1}+\lambda_{2}\right)^{2}\left(\lambda_{1}+2\lambda_{2}\right)^{4}e^{\lambda_{1}\xi}-16abk^{2}c_{2}^{2}\lambda_{2}^{2}\left(\lambda_{1}+2\lambda_{2}\right)\left\\{c_{1}k^{2}\lambda_{1}^{2}\right.\right.$
$\displaystyle\times\left.\left.\left(\lambda_{1}^{2}-4\lambda_{2}^{2}\right)^{2}\left(\lambda_{1}^{2}+\lambda_{2}^{2}\right)e^{\left(\lambda_{1}+2\lambda_{2}\right)\xi}+c_{1}^{2}\lambda_{2}^{2}\left(\lambda_{1}^{2}-3\lambda_{1}\lambda_{2}+2\lambda_{2}^{2}\right)^{2}\right.\right.$
$\displaystyle\times\left.\left.e^{2\left(\lambda_{1}+\lambda_{2}\right)\xi}+4k^{4}\lambda_{1}^{4}\lambda_{2}^{2}\left(\lambda_{1}^{2}+3\lambda_{2}\lambda_{1}+2\lambda_{2}^{2}\right)^{2}e^{2\lambda_{2}\xi}\right\\}+b^{2}c_{1}c_{2}^{4}\lambda_{1}^{4}\right.$
$\displaystyle\times\left.\left(\lambda_{1}-2\lambda_{2}\right)\left(\lambda_{1}-\lambda_{2}\right)^{2}e^{\left(\lambda_{1}+4\lambda_{2}\right)\xi}\right]\text{/}Q(\xi)^{2},$
(94) $\displaystyle V(\xi)=$ $\displaystyle
8ac_{2}k^{4}\lambda_{2}^{2}\left(\lambda_{1}-2\lambda_{2}\right)\left(\lambda_{1}+\lambda_{2}\right)\left(\lambda_{1}+2\lambda_{2}\right)^{2}e^{\lambda_{2}\xi}\left[c_{1}\left(\lambda_{1}^{2}-3\lambda_{2}\lambda_{1}\right.\right.$
$\displaystyle\left.\left.+2\lambda_{2}^{2}\right)e^{\lambda_{1}\xi}+2k^{2}\lambda_{1}^{2}\left(\lambda_{1}^{2}+3\lambda_{2}\lambda_{1}+2\lambda_{2}^{2}\right)\right]\text{/}Q(\xi),$
(95)

where

$\displaystyle Q(\xi)=$ $\displaystyle
8ak^{4}\lambda_{2}^{2}\left(\lambda_{1}-2\lambda_{2}\right)\left(\lambda_{1}+\lambda_{2}\right)^{2}\left(\lambda_{1}+2\lambda_{2}\right)^{3}\left(c_{1}e^{\lambda_{1}\xi}+2k^{2}\lambda_{1}^{2}\right)-bc_{2}^{2}e^{2\lambda_{2}\xi}$
$\displaystyle\times\left\\{c_{1}\left(\lambda_{1}^{2}-3\lambda_{2}\lambda_{1}+2\lambda_{2}^{2}\right)^{2}e^{\lambda_{1}\xi}+2k^{2}\lambda_{1}^{2}\left(\lambda_{1}^{2}+3\lambda_{2}\lambda_{1}+2\lambda_{2}^{2}\right)^{2}\right\\}.$
(96)

Equations (4)-(4) with (88) and (90) constitutes the solution of (1). These
solutions have been verified through symbolic computation.

##  5 Boundedness of solution

Bounded solutions have exigent bags of contribution in recounting perspective
of physical systems modelled by differential equations. This section is
devoted to deriving the boundedness conditions of solution (4)-(4) with (88)
and (90) of (1). To attain the goal we propose the following theorems:

######  Theorem 5.1.

The solution (4)-(4) of equation (89) is bounded in $\mathbb{R}$ if the
parameters $a,\ b,\ k,\ \lambda_{1},\ \lambda_{2}$ involved in the equation
and arbitrary constants $c_{1},\ c_{2}$ present in the solution satisfy any
one of the following conditions.
In case of  $0<\frac{1}{2}\lambda_{1}<\lambda_{2}<\lambda_{1}$:

I.(a) $c_{1}>0,\ c_{2}>0,\ a>0,\ b>0,\ k<0,$

I.(b) $c_{1}>0,\ c_{2}>0,\ a>0,\ b>0,\ k>0,$

I.(c) $c_{1}>0,\ c_{2}<0,\ a>0,\ b>0,\ k<0,$

I.(d) $c_{1}>0,\ c_{2}<0,\ a>0,\ b>0,\ k>0.$

In case of  $0<\lambda_{1}<\lambda_{2}$ :

II.(a) $c_{1}>0,\ c_{2}>0,\ a>0,\ b>0,\ k<0,$

II.(b) $c_{1}>0,\ c_{2}>0,\ a>0,\ b>0,\ k>0,$

II.(c) $c_{1}>0,\ c_{2}<0,\ a>0,\ b>0,\ k<0,$

II.(d) $c_{1}>0,\ c_{2}<0,\ a>0,\ b>0,\ k>0.$

In case of  $0<\lambda_{2}<\frac{1}{2}\lambda_{1}$ :

III.(a) $c_{1}>0,\ c_{2}>0,\ a>0,\ b<0,\ k<0,$

III.(b) $c_{1}>0,\ c_{2}>0,\ a>0,\ b<0,\ k>0,$

III.(c) $c_{1}>0,\ c_{2}<0,\ a>0,\ b<0,\ k<0,$

III.(d) $c_{1}>0,\ c_{2}<0,\ a>0,\ b<0,\ k>0.$

###### Proof.

Solution (4)-(4) have a common factor in denominators given by $Q(\xi).$
Boundedness of this solution demands that $Q(\xi)$ never vanishes (ie.,
$Q(\xi)>0$ or $Q(\xi)<0$). For $c_{1}>0$, we have
$\left(c_{1}e^{\lambda_{1}\xi}+2k^{2}\lambda_{1}^{2}\right)>0$ and

$\left\\{c_{1}\left(\lambda_{1}^{2}-3\lambda_{2}\lambda_{1}+2\lambda_{2}^{2}\right)^{2}e^{\lambda_{1}\xi}+2k^{2}\lambda_{1}^{2}\left(\lambda_{1}^{2}+3\lambda_{2}\lambda_{1}+2\lambda_{2}^{2}\right)^{2}\right\\}>0.$

Above presented conditions with $a>0$ lead us to the conditions for non-
vanishing $Q(\xi)$ as:

Case $Q(\xi)<0$: $\ \left\\{\left(\lambda_{1}-2\lambda_{2}\right)<0,\
b>0\right\\}\Rightarrow\
\left\\{\frac{1}{2}\lambda_{1}<\lambda_{2}<\lambda_{1},\ b>0\right\\}$ or
$\left\\{0<\lambda_{1}<\lambda_{2},\ b>0\right\\}.$

Case $Q(\xi)>0$: $\ \left\\{\left(\lambda_{1}-2\lambda_{2}\right)>0,\
b<0\right\\}\Rightarrow\ \left\\{0<\lambda_{2}<\frac{1}{2}\lambda_{1},\
b<0\right\\}.$

Case $Q(\xi)<0$ provide the conditions I.(a)$-$I.(d) and II.(a)$-$II.(d)
whereas Case $Q(\xi)>0$ give conditions III.(a)$-$III.(d) of the theorem. That
completes the proof of the theorem. ∎

To check the correctness of the conditions presented in the above theorem, few
particular values of parameters satisfying above conditions have been
presented in the table 1. These values have been utilized in figure 1 to plot
the solution (4)-(4). It is clear from the figures that the solution $U(\xi)$
(given in (4)) have three-pick soliton like shape for the conditions
I.(a)$-$I.(d) and II.(a)$-$II.(d) respectively and for remaining cases, it
bears two-pick soliton form. However $V(\xi)$ (given in (4)) always restrain
two-pick soliton form. Besides that, it is important to note here that these
solutions may bear a lower number of picks for different parameters values
satisfying above conditions.

Table 1:  Values of parameters present in $U(\xi),\ V(\xi)$ satisfying
conditions presented in Theorem 5.1, used in Figure 1.

Case | $\lambda_{1}$ | $\lambda_{2}$ | $a$ | $b$ | $k$ | $c_{1}$ | $c_{2}$ | Figure
---|---|---|---|---|---|---|---|---
| $0.4$ | $0.399$ | $1.6$ | $0.1$ | $-0.25$ | $1.9$ | $0.5$ | I.(a)
$\frac{1}{2}\lambda_{1}<\lambda_{2}<\lambda_{1}$ | $1.43$ | $1.42$ | $1$ | $3.1$ | $0.5$ | $4.9$ | $0.6$ | I.(b)
$a>0,\ b>0$ | $0.4$ | $0.399$ | $1.6$ | $0.1$ | $-0.25$ | $1.9$ | $-0.5$ | I.(c)
| $1.43$ | $1.42$ | $1$ | $3.1$ | $0.5$ | $4.9$ | $-0.6$ | I.(d)
| $0.38888$ | $0.3889$ | $2.1$ | $0.001$ | $-0.02$ | $8.9$ | $9.4$ | II.(a)
$0<\lambda_{1}<\lambda_{2}$ | $1.08$ | $1.0809$ | $0.005$ | $0.0055$ | $0.2$ | $0.039$ | $0.003$ | II.(b)
$a>0,\ b>0$ | $0.38888$ | $0.3889$ | $2.1$ | $0.001$ | $-0.02$ | $8.9$ | $-9.4$ | II.(c)
| $1.08$ | $1.0809$ | $0.005$ | $0.0055$ | $0.2$ | $0.039$ | $-0.003$ | II.(d)
| $1.11$ | $0.42$ | $1.2$ | $-3.1$ | $-1.1$ | $2.3$ | $0.3$ | III.(a)
$0<\lambda_{2}<\frac{1}{2}\lambda_{1}$ | $2$ | $0.99$ | $1.2$ | $-3.1$ | $1.1$ | $1.3$ | $0.3$ | III.(b)
$a>0,\ b<0$ | $1.11$ | $0.42$ | $1.2$ | $-3.1$ | $-1.1$ | $2.3$ | $-0.3$ | III.(c)
| $2$ | $0.99$ | $1.2$ | $-2.1$ | $1.1$ | $5.3$ | $-0.6$ | III.(d)

Figure 1: Plots of solution (4)-(4) for the conditions presented in Theorem
5.1, using values of parameters given in Table 1.

######  Theorem 5.2.

The solution (4)-(4) with (88) and (90) of equation (1) is bounded in
$\mathbb{R}\times\mathbb{R}$ if the parameters $a,\ b,\ c,\ k$ involved in the
equation and arbitrary constants $c_{1},\ c_{2}$ present in the solution
satisfy any one of the following conditions.
In case of
$\frac{1}{2}\sqrt{-\frac{c}{ak^{3}}}<\sqrt{-\frac{c}{k^{3}}}<\sqrt{-\frac{c}{ak^{3}}}$:

I.(a) $c_{1}>0,\ c_{2}>0,\ a>0,\ b>0,\ k<0,\ c>0,$

I.(b) $c_{1}>0,\ c_{2}>0,\ a>0,\ b>0,\ k>0,\ c<0,$

I.(c) $c_{1}>0,\ c_{2}<0,\ a>0,\ b>0,\ k<0,\ c>0,$

I.(d) $c_{1}>0,\ c_{2}<0,\ a>0,\ b>0,\ k>0,\ c<0.$

In case of  $0<\sqrt{-\frac{c}{ak^{3}}}<\sqrt{-\frac{c}{k^{3}}}$ :

II.(a) $c_{1}>0,\ c_{2}>0,\ a>0,\ b>0,\ k<0,\ c>0,$

II.(b) $c_{1}>0,\ c_{2}>0,\ a>0,\ b>0,\ k>0,\ c<0,$

II.(c) $c_{1}>0,\ c_{2}<0,\ a>0,\ b>0,\ k<0,\ c>0,$

II.(d) $c_{1}>0,\ c_{2}<0,\ a>0,\ b>0,\ k>0,\ c<0.$

In case of  $0<\sqrt{-\frac{c}{k^{3}}}<\frac{1}{2}\sqrt{-\frac{c}{ak^{3}}}$ :

III.(a) $c_{1}>0,\ c_{2}>0,\ a>0,\ b<0,\ k<0,\ c>0,$

III.(b) $c_{1}>0,\ c_{2}>0,\ a>0,\ b<0,\ k>0,\ c<0,$

III.(c) $c_{1}>0,\ c_{2}<0,\ a>0,\ b<0,\ k<0,\ c>0,$

III.(d) $c_{1}>0,\ c_{2}<0,\ a>0,\ b<0,\ k>0,\ c<0.$

###### Proof.

This theorem can be proved in a straightforward way by using equations
(88)-(90) and theorem 5.1. ∎

To verify the conditions presented in theorem 5.2, some particular values of
parameters satisfying the above conditions have been exhibited in table 2.
Utilizing those values, 3D plots of the solution (4)-(4) with (88) and (90)
for different cases have displayed in figure 2. Figures reveal that the
solution $u(x,t)$ (given by (4)-(4) with (88) and (90)) always have three-hump
soliton like features for the conditions I.(a)$-$I.(d) and II.(a)$-$II.(d)
respectively whereas it bears two-hump soliton shape for remaining cases of
theorem 5.2. However $v(x,t)$ (given by (4)-(4) with (88) and (90)) always
restrain two-pick form. Besides that, it is important to mention that these
solutions may bear a lower number of humps. From the 3D plots it is clear that
solution can describe various multi-pick soliton states for diverse values of
the corresponding free parameters.

The existing multi-pulse solutions in the literature [78, 2], are constituted
with multiple copies of the one-pulse solutions separated by finitely many
oscillations close to the zero equilibrium (has oscillatory decaying tails at
infinity). It is worth noting that here obtained multi-pulse solutions
((4)-(4)) of equation (1) constitute multiple picks but do not oscillate close
to the zero equilibrium. This is due to the fact that auxiliary equations of
the linear parts of the considered coupled nonlinear differential equations
(89) have distinct real roots $(\lambda_{1}\neq\lambda_{2})$. That ensures the
nonexistence of any periodic functions in the derived solution. So in view of
study [2] and above property lead us to conclude that derived solutions (Eq.
(4)-(4)) are homoclinic and always have finite numbers of humps and their
tails decay to zero exponentially in a monotonic fashion, as depicted in
figures 1 and 2. This fact, differentiate the derived solutions (exact) of
this paper with the existing multi-hump stationary wave solutions (numerical)
in literature [78, 2]. Also, it is important to mention here that usually
soliton trains arrive straight in shape but in figure 2 they are appeared in
bent shape due to the existence of conformable derivative parameter $\alpha$
and $\beta$ in the solution.

Table 2:  Values of parameters parameters in $u(x,t),\ v(x,t)$ satisfying
conditions presented in Theorem 5.2, used in Figure 2.

Case | $a$ | $b$ | $c$ | $k$ | $c_{1}$ | $c_{2}$ | $\alpha$ | $\beta$ | $\xi_{0}$ | Figure
---|---|---|---|---|---|---|---|---|---|---
| $0.990$ | $0.09$ | $1.0$ | $-1.8$ | $4.6$ | $.58$ | $.68$ | $.69$ | $28$ | I.(a)
$\frac{1}{2}\sqrt{-\frac{c}{ak^{3}}}<\sqrt{-\frac{c}{k^{3}}}<\sqrt{-\frac{c}{ak^{3}}}$ | $0.980$ | $0.09$ | $-1.6$ | $2.1$ | $2.8$ | $.58$ | $.70$ | $.75$ | $28$ | I.(b)
$a>0,\ b>0$ | $0.989$ | $0.02$ | $1.7$ | $-1.9$ | $3.8$ | $-.9$ | $.82$ | $.80$ | $28$ | I.(c)
| $0.999$ | $0.02$ | $-1.7$ | $1.9$ | $3.0$ | $-.8$ | $.75$ | $.85$ | $28$ | I.(d)
| $1.001$ | $0.001$ | $1.0$ | $-1.49$ | $6.9$ | $3.4$ | $.79$ | $.69$ | $12$ | II.(a)
$0<\sqrt{-\frac{c}{ak^{3}}}<\sqrt{-\frac{c}{k^{3}}}$ | $1.002$ | $0.002$ | $-1.5$ | $2.9$ | $.0069$ | $.0035$ | $.80$ | $.92$ | $12$ | II.(b)
$a>0,\ b>0$ | $1.001$ | $0.0007$ | $2.0$ | $-2.5$ | $.0069$ | $-.0035$ | $.88$ | $.82$ | $12$ | II.(c)
| $1.003$ | $0.001$ | $-1.5$ | $2.9$ | $.0068$ | $-.0034$ | $.80$ | $.90$ | $12$ | II.(d)
| $0.20$ | $-3.0$ | $2.1$ | $-2.0$ | $3.1$ | $.3$ | $.70$ | $.80$ | $4$ | III.(a)
$0<\sqrt{-\frac{c}{k^{3}}}<\frac{1}{2}\sqrt{-\frac{c}{ak^{3}}}$ | $0.22$ | $-2.4$ | $-2.0$ | $1.4$ | $2.9$ | $.5$ | $.80$ | $.90$ | $4$ | III.(b)
$a>0,\ b<0$ | $0.23$ | $-2.6$ | $2.2$ | $-1.3$ | $3.3$ | $-.4$ | $.90$ | $.98$ | $4$ | III.(c)
| $0.24$ | $-2.5$ | $-2.0$ | $1.3$ | $2.9$ | $-.5$ | $.80$ | $.85$ | $4$ | III.(d)

Figure 2: Plots of solution $u(x,t),\ v(x,t)$ for the conditions presented in
Theorem 5.2, using values of the parameters submitted in Table 2.

##  6  Conclusions

In this study, the RCAM was applied to a coupled Korteweg-de Vries equations
with conformable derivative, which is the generalisation of the mathematical
model of waves of shallow water surface equations. A novel exact solution in
terms of exponential function has been derived. Additionally, a few theorems
have been presented to predicts the boundedness of the obtained solution. All
families of boundedness conditions on the parameters present in the equation
and integration constants present in the solutions were tested by plot,
affirming their correctness. It was found that there exist triple-pulse,
double-pulse, and single-pulse solitons for the considered equation. The
author believes that it is a new finding for the coupled Korteweg-de Vries
equations with conformable derivative. In addition to that here for the first
time, a modification of RCAM is proposed to deal with couple nonlinear
differential equations whose auxiliary equations of linear parts have distinct
roots ($\lambda_{1}\neq\lambda_{2}$), and successfully derived exact multi-
hump solutions of coupled KdV Eq. (1). That extends its applicability to deal
with complex nonlinear equations and produce new featured solutions like
multi-hump solitons.

The reported results of the present work contribute a new stationary wave
solution to the class of coupled nonlinear systems with conformable derivative
that may have a composition of the different multi-hump soliton-type features
depending on the several restrictions in the parameters present on the
equations. Here we have successfully derived three-hump, two-hump, and one-
hump soliton solutions of Eq. (1), but unable to obtain the solitons having
four and higher-order hump due to the unavailability of more wave parameters.
It is straightforward to conclude that the rapidly convergent approximation
method is a powerful tool to solve various nonlinear models involving
conformable derivative.

The modification of RCAM to reveal solitary wave and multi-soliton solutions
of constant and variable coefficient differential equations is significant in
the studied field and might have an important impact on future research. So in
the near future, the author wants to modify the scheme for exhaling solitary
wave and multi-soliton solutions of constant and variable coefficient
differential equations.

## Acknowledgements

The author thanks the editor and reviewers for their comments and suggestions
to improve the paper in the revised form.

## Conflict of interest

The authors declare that they have no conflict of interest.

## References

## References

  * [1] Buffoni B and Séré E 1996 Communications on pure and applied mathematics 49 285–305
  * [2] Groves M 1998 Nonlinearity 11 341
  * [3] Gorshkov K, Ostrovsky L, Papko V and Pikovsky A 1979 Physics Letters A 74 177–179
  * [4] Ostrovskaya E A, Kivshar Y S, Skryabin D V and Firth W J 1999 Physical review letters 83 296
  * [5] Ostrovskaya E A, Mingaleev S F, Kivshar Y S, Gaididei Y B and Christiansen P L 2001 Physics Letters A 282 157–162
  * [6] Cattani F, Kim A, Hansson T, Anderson D and Lisak M 2011 EPL (Europhysics Letters) 94 53003
  * [7] Parra Prado H and Cisneros-Ake L A 2019 Chaos: An Interdisciplinary Journal of Nonlinear Science 29 053133
  * [8] Wang L, Li S and Qi F H 2016 Nonlinear Dynamics 85 389–398
  * [9] Magin R L 2006 Fractional calculus in bioengineering (Begell House Redding)
  * [10] Engheta N 1996 IEEE Transactions on Antennas and Propagation 44 554–566
  * [11] Schneider W R and Wyss W 1989 Journal of Mathematical Physics 30 134–144
  * [12] Chen Y, Sun R and Zhou A 2005 Fractional order calculus day at Utah State University
  * [13] Meral F, Royston T and Magin R 2010 Communications in Nonlinear Science and Numerical Simulation 15 939–945
  * [14] He J H 1997 International Journal of Turbo and Jet Engines 14 23–28
  * [15] He J H 2004 Chaos, Solitons & Fractals 19 847–851
  * [16] Banerjee J, Ghosh U, Sarkar S and Das S 2017 Pramana 88 70
  * [17] Das T, Ghosh U, Sarkar S and Das S 2018 Journal of Mathematical Physics 59 022111
  * [18] Kilbas A A A, Srivastava H M and Trujillo J J 2006 Theory and applications of fractional differential equations vol 204 (Elsevier Science Limited)
  * [19] Kumar D, Singh J, Tanwar K and Baleanu D 2019 International Journal of Heat and Mass Transfer 138 1222–1227
  * [20] Kumar D, Singh J and Baleanu D 2020 Mathematical Methods in the Applied Sciences 43 443–457
  * [21] Prakash A and Kaur H 2019 Chaos, Solitons & Fractals 124 134–142
  * [22] Prakash A, Kumar M and Baleanu D 2018 Applied Mathematics and Computation 334 30–40
  * [23] Goyal M, Baskonus H M and Prakash A 2020 Chaos, Solitons & Fractals 139 110096
  * [24] Aslan E C and Inc M 2019 Optik 196 162661
  * [25] Inc M, Yusuf A, Aliyu A I and Baleanu D 2018 Optik 162 65–75
  * [26] Yusuf A and Inc M 2020 Physica Scripta 95 035217
  * [27] Khalil R, Al Horani M, Yousef A and Sababheh M 2014 Journal of Computational and Applied Mathematics 264 65–70
  * [28] Pandir Y and Yildirim A 2018 Waves in Random and Complex Media 28 399–410
  * [29] Ellahi R, Mohyud-Din S T, Khan U et al. 2018 Results in physics 8 114–120
  * [30] Sabi’u J, Jibril A and Gadu A M 2019 Journal of Taibah University for Science 13 91–95
  * [31] Houwe A, Sabi’u J, Hammouch Z and Doka S Y 2020 Physica Scripta 95 045203
  * [32] Rezazadeh H, Osman M, Eslami M, Mirzazadeh M, Zhou Q, Badri S A and Korkmaz A 2019 Nonlinear Engineering 8 224–230
  * [33] Osman M, Rezazadeh H and Eslami M 2019 Nonlinear Engineering 8 559–567
  * [34] Hashemi M S, Inc M and Yusuf A 2020 Chaos, Solitons & Fractals 133 109628
  * [35] Korpinar Z, Inc M and Bayram M 2020 Applied Mathematics and Computation 367 124781
  * [36] Korpinar Z, Inc M, Hınçal E and Baleanu D 2020 Alexandria Engineering Journal 59 1405–1412
  * [37] Korpinar Z and Inc M 2018 Optik 166 77–85
  * [38] Yusuf A, Bayram M et al. 2019 Physica Scripta 94 125005
  * [39] Yong C and Hong-Li A 2008 Communications in Theoretical Physics 49 839
  * [40] Atangana A and Secer A 2013 The time-fractional coupled-Korteweg-de-Vries equations Abstract and Applied Analysis vol 2013 (Hindawi)
  * [41] Biswas S, Ghosh U, Sarkar S and Das S 2020 Journal of the Physical Society of Japan 89 014002
  * [42] Osman M 2019 Pramana 93 26
  * [43] Arqub O A, Osman M S, Abdel-Aty A H, Mohamed A B A and Momani S 2020 Mathematics 8 923
  * [44] Kumar D, Singh J, Purohit S D and Swroop R 2019 Mathematical Modelling of Natural Phenomena 14
  * [45] Goswami A, Singh J, Kumar D et al. 2019 Physica A: Statistical Mechanics and its Applications 524 563–575
  * [46] Prakash A, Goyal M and Gupta S 2019 Pramana 93 28
  * [47] Das P K and Panja M 2015 An improved Adomian decomposition method for nonlinear ODEs Applied Mathematics (Springer) pp 193–201
  * [48] Das P K and Panja M 2016 IJSEAS 2 334–348
  * [49] Das P K 2018 Sohag J. Math. 5 29–33
  * [50] Das P K, Singh D and Panja M 2018 Optik 174 433–446
  * [51] Das P K, Mandal S and Panja M M 2018 Mathematical Methods in the Applied Sciences 41 7869–7887
  * [52] Das P K, Singh D and Panja M M 2019 Journal of Advances in Mathematics 16 8213–8225
  * [53] Das P K 2019 Optik 195 163134
  * [54] Das P K 2020 Optik 223 165293
  * [55] Jin-Cun L and Guo-Lin H 2010 Chinese Physics B 19 110203
  * [56] Ray S S 2013 International Journal of Nonlinear Sciences and Numerical Simulation 14 501–511
  * [57] Bulut F, Oruç Ö and Esen A 2015 Computer Modeling in Engineering & Sciences 108 263–284
  * [58] Matinfar M, Eslami M and Kordy M 2015 Pramana 85 583–592
  * [59] Bhrawy A, Doha E, Ezz-Eldien S and Abdelkawy M 2016 Calcolo 53 1–17
  * [60] Hussain M, Haq S and Ghafoor A 2019 Applied Mathematics and Computation 341 321–334
  * [61] Albuohimad B, Adibi H and Kazem S 2018 Ain Shams Engineering Journal 9 1897–1905
  * [62] Tarasov V E 2016 Communications in Nonlinear Science and Numerical Simulation 30 1–4
  * [63] Tarasov V E 2013 Communications in Nonlinear Science and Numerical Simulation 18 2945–2948
  * [64] Tarasov V E 2018 Communications in Nonlinear Science and Numerical Simulation 62 157–163
  * [65] Hirota R and Satsuma J 1981 Physics Letters A 85 407–408
  * [66] Zhou Y, Wang M and Wang Y 2003 Physics Letters A 308 31–36
  * [67] Xie Y 2004 Physics Letters A 327 174–179
  * [68] Singh K and Gupta R 2006 International journal of engineering science 44 241–255
  * [69] Abdeljawad T 2015 Journal of computational and Applied Mathematics 279 57–66
  * [70] Eslami M 2016 Applied Mathematics and Computation 285 141–148
  * [71] Chen C and Jiang Y L 2018 Computers & Mathematics with Applications 75 2978–2988
  * [72] Adomian G 1994 Solving Frontier Problems Decomposition Method
  * [73] Duan J S and Rach R 2011 Applied Mathematics and Computation 218 4090–4118
  * [74] Adomian G and Rach R 1983 Journal of Mathematical Analysis and Applications 91 39–46
  * [75] Adomian G and Rach R 1993 Journal of Mathematical Analysis and Applications 174 118–137
  * [76] Adomian G and Rach R 1993 Applied mathematics and computation 57 61–68
  * [77] Adomian G and Rach R 1994 Nonlinear Analysis: Theory, Methods & Applications 23 615–619
  * [78] Chugunova M and Pelinovsky D 2007 Discrete & Continuous Dynamical Systems-B 8 773