%\documentclass[twocolumn,showpacs,showkeys,preprintnumbers]{revtex4}
%\documentclass[twocolumn,showkeys,preprintnumbers]{revtex4}
\documentclass[14pt,showkeys,preprintnumbers]{revtex4}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{dcolumn}
\usepackage{bm}
\usepackage[colorlinks]{hyperref}
\usepackage{color}
\newcommand{\asinh}{\mbox{\textrm{\ensuremath{\,}asinh}}}
\newcommand{\sech}{\mbox{\textrm{\ensuremath{\,}sech}}}
\newcommand{\csch}{\mbox{\textrm{\ensuremath{\,}csch}}}
\setcounter{MaxMatrixCols}{10}


\begin{document}
%
%\maketitle

%\setlinespacing{1.75}
%\title{\textbf{\ Higher-order dispersion effect in the modulational instability spectrum
%of a relaxing nonlinear oppositely directed coupler with negative index material channel}}
\title{\textbf{\ Effects of saturable Function in three-core PIM-NIM-PIM coupler through Modulation instability}}


\author{P.H. Tatsing$^{1,2}$, A.C. Chamgoue$^{3}$, E. Kengne$^{4,5}$, A. Mohamadou$^{2,6,7}$ and T.C. Kofane$^{1,2,7,8}$\\
%\date{{
%\begin{flushleft}
%\begin{center}
${}^{1}${\footnotesize \textit{Faculty of Science, University of Yaounde I, P.O. Box 812, Yaounde, Cameroon.}}\\
${}^{2}${\footnotesize \textit{Centre d`Excellence Africain en Technologies de l`Information et de la Communication, Universite de Yaounde, Cameroon.}}\\
${}^{3}${\footnotesize \textit{School of Geology and Mining Engineering, University of Ngaoundere, P.O. Box , Ngaoundere, Cameroon.}}\\
${}^{4}${\footnotesize \textit{School of Physics and Electronic information Engineering, Zhejiang Normal University Jinhua, 321004, China.}}\\
${}^{5}${\footnotesize \textit{Department of Computer and Engineering, University of Quebec at Outaouais, 101 St-Jean-Bosco
 , Canada J8Y 3G5.}}\\
${}^{6}${\footnotesize \textit{Complex Systems, National Polytechnic School, University of Maroua, P.O. Box 46, Maroua, Cameroon.}}\\
${}^{7}${\footnotesize \textit{The Abdus Salam International Centre for Theoretical Physics, P.O. Box 586, Strada costiera 11, I-34014, Trieste, Italy.}}\\
${}^{8}${\footnotesize \textit{Botswana International University of Science and Technology, Private Bag 16, Palapye, Botswana.}}\\
%${}^{5}${\footnotesize \textit{Univ Lille 1, CNRS, UMR 8523, Lab Phys Lasers Atomes and Mol, F-59655 Villeneuve Dascq, France.}}\\
%${}^{6}${\footnotesize \textit{Department of Physics, Pondicherry University, Pondicherry 605014, India.}}\\
%${}^{5}${\footnotesize \textit{Centre for Plasma Physics, Department of Physics and Astronomy, Queen's University Belfast, BT7 1NN Northern Ireland, UK.}}\\
}
%\end{center}

%\end{flushleft}
%}
%}}

\begin{abstract}

We consider a model of three-core PIM-NIM-PIM coupler with Kerr-type saturable nonlinearity to study both analytically and numerically saturation effects on the modulational instability phenomenon. The analytical results show us that, in presence of saturable parameter in normal and anomalous dispersion regime it is formed some new instability bands compare to Shafeeque et $al$ (2015) in absence of saturable parameter. It is clearly also show that Power and nonlinears parameters can be used to control MI. Through numerical simulation, the generation of periodic soliton with a growing amplitude are obtained in absence of nonlinear saturation. However in presence of saturable parameter the train of soliton obtained turns into a turbulent state after certain propagations distance.
%We investigate modulational instability (MI) in nonlinear three-core coupler with negative-index metamaterial channel and modified kerr-type saturable nonlinearity (MSN). Using standard linear stability analysis we obtained the instability gain. We study in particular the combination of saturable nonlinearity with nonlinear parameters or forward to backward-propagating wave's power $\emph{f}$ effects on MI both normal and anomalous group velocity dispersion regimes. We observe that the instability gain exhibits significant changes due to the effects of the saturable nonlinearity. The control of the MI can be realized by adjusting the nonlinear parameters,backward-propagating wave's power $\emph{f}$, even in the presence of a saturable nonlinearity.
%At end, numerical simulations are carried out to explore nonlinear development of the MI, revealing the generation of periodic chains of localized peaks with growing amplitudes, which may transform into arrays of solitons.

\end{abstract}
%
\maketitle

{\textbf{Keywords:}}
Modulational instability, Nonlinear three-core PIM-NIM-PIM couplers (PIM: Positive Index Material; NIM: Negative Index Material); Saturable nonlinearity, Split-step Fourier method (SSFM).
%\pacs{000}
{\textbf{Electronic address:}}
 tatsingp@yahoo.fr


\section{Introduction}\label{sec:1}


Negative index material (NIM) is material which both the permittivity $(\varepsilon)$ and permeability $(\mu)$ parameters are set as negative at the same frequency \cite{veselago}.  The existence of such media was experimentally  demonstrated first in the microwave \cite{shelby} and then in the near-IR ranges \cite{linden,zhang1}.  NIM or metamaterial can be artificially designed and their properties  derive from their structures but not from the properties of the base materials. NIM research is interdisciplinary and include many fields such as optics, optoelectronics, nanoscience, antenna and electrical engineering \cite{said}.  In optics negative index material offer many potential applications such as optical filters, lenses for high-gain antennas, improving ultrasonic sensors, super lenses etc \cite{boratay}.

In nonlinear optics, the waveguiding structure composed of two adjacent waveguides preserving the direction of light propagation  is called a directed coupler. Jensen \cite{jensen} was the first researcher to introduce the notion of nonlinear directional coupler. Nonlinear directional coupler has applications in optical communication systems such as switching, wavelength-selective coupling, multi/demultiplexing and power splitting \cite{jensen,kengne}.  if one of the waveguides of the coupler is made from a material with a negative refractive index, the direction of input and output fields are exactly opposite in nature but the path of light propagation is preserves. This coupler is called the oppositely directed coupler in order to distinguish from the conventional one made from a material with a positive refractive index (PIM) \cite{litchinitser}. Such structure is first introduced by Halterman et $al$. in 2003 \cite{halterman}, then demonstrated  experimentally by Yuan et $al$. In 2006 \cite{yuan}. The interaction between the nonlinear effect and group-velocity dispersion (GVD)  is named MI. In nonlinear wave systems modulational instability is the most fundamental processes in nature, characterized by a continuous wave (CW) when it propagates together with a weak noise \cite{agrawal}. Recent years, some researchers  investigated the modulational instability (MI) in a nonlinear oppositely directed coupler with a negative-index material channel \cite{xiang, shafeeque1, zhang, shafeeque4, kengne2}. MI in nonlinear positive-negative index couplers with saturable nonlinearity \cite{tatsing, alves, tatsing2, tatsing3, tatsing4,houwe}. The influence of self-steepening and intrapulse Raman scattering on MI in oppositely directed coupler \cite{shafeeque1}.

 Another form of optical couplers are multi-core directional couplers that consist of more than two waveguides. Such couplers improved transmission characteristics that why they attracted much attention. Sharper power switching curves are well offer by three core couplers than conventional two-core coupler and for that raison three core couplers are really  significant \cite{langridge, shafeeque2, shafeeque3}. Many physical systems present nonlinear saturation that may influence them. Some important changes in instability band particulary in shape and/or amplitude can be observe due to the presence of nonlinear saturation in the systems. This has encouraged some researchers to study the MI in many saturable nonlinear systems  \cite{tatsing, alves, tatsing2, tatsing3, tatsing4,houwe}. The present work is to study the effects of a saturable nonlinearity on the MI in a three-core coupler with a negative material channel. The rest of the paper is organized as follows : Section2 present mathematical model and standard linear stability analysis. In Section3  we carried out in detail the influence of the parameters of the three channels of the coupler on the MI. In Section4, numerical simulation are presented and finally, in section 5 we concludes the paper.

\section{Theoretical background and linear stability analysis}
%\setlinespacing{1.9} \Large
%\bigskip
\subsection{Propagation equations and dispersion relation}

 Shafeeque et $al.$\cite{shafeeque2} introduced continuous wave propagation in a nonlinear three-core optical coupler containing a
 NIM channels by neglecting the cross-phase modulation (XPM) effect and high order time derivative terms.
 However the model present by Shafeeque et $al.$\cite{shafeeque2} does not include saturation effects, whereas saturation nonlinearity play a relevant role in the propagation of ultrashort pulse.  In presence of saturable nonlinearity function the above model is considered as follows.
  %which shall play a relevant role in the propagation of ultrashort
%pulses, in most studies of saturation behavior, the nonlinear term in the Shr$\ddot{o}$dinger equation is replaced by \cite{xianqiong}:



\begin{equation}\label{eq1}
    i\sigma_{1}\frac{\partial a_{1}}{\partial z} +
    i\frac{1}{v_{1g}}\frac{\partial a_{1}}{\partial t}
    +k_{12}a_{2}exp(-i\delta z) +
    \gamma_{1}g(\Gamma\left|a_{1}\right|^{2})a_{1} = 0 ,
\end{equation}
\begin{equation}\label{eq2}
    i\sigma_{2}\frac{\partial a_{2}}{\partial z} +
    i\frac{1}{v_{2g}}\frac{\partial a_{2}}{\partial t}
    +k_{21}a_{1}exp(i\delta z) +k_{23}a_{3}exp(i\delta z) +
    \gamma_{2}g(\Gamma\left|a_{2}\right|^{2})a_{2} = 0 ,
\end{equation}
\begin{equation}\label{eq3}
    i\sigma_{3}\frac{\partial a_{3}}{\partial z} +
    i\frac{1}{v_{3g}}\frac{\partial a_{3}}{\partial t}
    +k_{32}a_{2}exp(-i\delta z) +
    \gamma_{3}g(\Gamma\left|a_{3}\right|^{2})a_{3} = 0 .
\end{equation}

Here channel 1 and channel 3 are PIM and channel 2 is NIM. $\sigma_{1}$, $\sigma_{2}$ and $\sigma_{3}$ stand for sign of the refractive index. In this work, $\sigma_{1}=\sigma_{3}=1$ and $\sigma_{2}=-1$.
 $a_{1}$, $a_{2}$ and $a_{3}$  are the complex normalized amplitudes; The absolute values of
the group velocities for channel 1,2 and 3 are given by $v_{1g}$, $v_{2g}$ and $v_{3g}$. $k_{12}$, $k_{21}$, $k_{23}$ and $k_{32}$ are the coupling coefficients respectively. $\gamma_{1}$, $\gamma_{2}$ and $\gamma_{3}$ are nonlinear coefficient.
%$\gamma_{i}=\frac{\omega_{0}n_{2j}\mu_{j}(\omega_{0})P_{0}}{cA_{ff}}$
%(i=1,2 and 3) is the normalized nonlinearity coefficient,where
%$n_{2j}=\frac{12\Pi^{2}\chi_{j}^{(3)}}{\varepsilon_{j}(\omega_{_{0}})c}$,$\chi_{j}^{(3)}$
%is the refractive index.
%In this paper. $\delta=\beta_{1}-\beta_{2}$=$\beta_{3}-\beta_{2}$ represents the mismatch between the propagation
%constants in the individual channels.

The Modified Kerr type saturable nonlinearity function used here have the following expression:

\begin{equation}\label{eq4}
  g(\Gamma\left|a_{i}\right|^{2}) = \frac{\left|a_{i}\right|^{2}(2+\Gamma\left|a_{i}\right|^{2}/2)}{2(1+\Gamma\left|a_{i}\right|^{2}/2)^{2}},
\end{equation}

where $\Gamma = 1/P_{sat}$ is the saturation parameter with $P_{sat}$ the saturable power density and i= 1,2 and 3.
%For field intensities such that
%$\Gamma \left|a_{i}\right|^{2} \ll 1$, the system depicts the usual Kerr response. However, for a very intense field
%$\Gamma\left|a_{i}\right|^{2}\gg 1$, the dependence of the refractive index on the field intensity saturates.

The following equations are the solutions of Eqs.(1), (2) and (3).
\begin{equation}\label{eq8}
    a_{1}=u_{1}exp(-i\frac{\delta}{2}z)exp(iqz), \nonumber
\end{equation}
\begin{equation}\label{eq8}
    a_{2}=u_{2}exp(i\frac{\delta}{2}z)exp(iqz),
\end{equation}
\begin{equation}\label{eq8}
    a_{3}=u_{3}exp(-i\frac{\delta}{2}z)exp(iqz). \nonumber
\end{equation}

 %Formation of the bandgap in a uniform structure considered here is one of the
%unique properties of the PIM-NIM coupler arising from introducing
%of the NIM into one channel of the coupler.This feature is similar
%to the property of the Bragg gratings or distributed feedback
%structures. We derive the nonlinear dispersion relations of
%Eqs.(1), (2) and (3) as:
 Let’us now insert Eqs.(5) in Eqs. ((1) and (2)), we obtain the nonlinear dispersion relations.

\begin{equation}\label{eq9}
    \delta = -\frac{(k_{12}f^{2}+k_{21})l+k_{23}f}{fl}-Q_{1}(\gamma_{1}+\gamma_{2}f^{2})+
    \frac{3}{4}\Gamma Q_{1}^{2}(\gamma_{1}+\gamma_{2}f^{4})+ \frac{\Gamma^{2}}{4}Q_{1}^{3}(\gamma_{1}+\gamma_{2}f^{6}),
\end{equation}
%\begin{equation}\label{eq9}
%  \frac{\Gamma^{2}}{4}Q_{1}^{3}(\gamma_{1}+\gamma_{2}f^{6}),
%\end{equation}


\begin{equation}\label{eq10}
    q = -\frac{(k_{21}-f^{2}k_{12})l+k_{23}f}{2fl}-\frac{Q_{1}}{2}(\gamma_{1}-\gamma_{2}f^{2})+
    \frac{3}{8}\Gamma Q_{1}^{2}(\gamma_{2}f^{4}-\gamma_{1})+\frac{\Gamma^{2}}{8}Q_{1}^{3}(\gamma_{2}f^{6}-\gamma_{1}).
\end{equation}
%\begin{equation}\label{eq10}
%  \frac{\Gamma^{2}}{8}Q_{1}^{3}(\gamma_{2}f^{6}-\gamma_{1}).
%\end{equation}

Where $f=u_{2}/u_{1}$ and $l= \frac{u_{2}}{u_{3}}$, are the quantities that describe how the Power $P=u_{1}^{2}+ u_{2}^{2}+ u_{3}^{2}$ is divided between the forward-and backward propagating waves. Also $Q_{1}= \frac{P}{1+f^{2}+\frac{f^{2}}{l^{2}}}$.



\subsection{Linear stability analysis}
 %The MI gain spectrum was derive by the following standard procedure \cite{shafeeque2,litsinistser1}. We
%assume that the solutions to the governing Eqs. (1), (2) and (3) are
%perturbed slightly such that:

To observe verywell the effect of the saturable nonlinearity on MI, linear stability analysis are employ. The first step consists to add small
perturbations terms to a CW and then survey if the latter augments or disintegrate with propagation.
When we take into account the small perturbations, Eqs. (5) turn to
\begin{equation}\label{eq11}
a_{1}=(u_{1}+\xi_{1})exp(iqz)exp(-i\frac{\delta}{2}z),   \nonumber
\end{equation}
\begin{equation}\label{eq11}
 a_{2}=(u_{2}+\xi_{2})exp(iqz)exp(i\frac{\delta}{2}z),
\end{equation}
\begin{equation}\label{eq11}
    a_{3}=(u_{3}+\xi_{2})exp(iqz)exp(-i\frac{\delta}{2}z). \nonumber
\end{equation}

Where $\xi_{i}$ is a term of perturbation ($\left|\xi_{i}\right|\ll
u_{i}, i= 1, 2 $ and 3), $\xi_{i}$ is small.

To obtain the linear form, we substitute Eqs.(8) into Eqs. (\ref{eq1}), (\ref{eq2}) and (\ref{eq3}) and linearize in
$\xi_{i}$

\begin{equation}\label{eq12}
\begin{array}{l}
i\frac{\partial \xi_{1}}{\partial z} +
    i\frac{1}{v_{1g}}\frac{\partial \xi_{1}}{\partial t}
    +k_{12}\xi_{2}-k_{12}f\xi_{1} + N_{1}[(\xi_{1}+\xi_{1}^{\ast})]= 0,
\end{array}
\end{equation}


\begin{equation}\label{eq13}
\begin{array}{l} -i\frac{\partial \xi_{2}}{\partial z} +
    i\frac{1}{v_{2g}}\frac{\partial \xi_{2}}{\partial t}
    +k_{21}\xi_{1}-k_{21}f^{-1}\xi_{2} +k_{23}\xi_{3}-k_{23}l^{-1}\xi_{2}+ N_{2}[(\xi_{2}+\xi_{2}^{\ast})] = 0,
\end{array}
\end{equation}

\begin{equation}\label{eq14}
\begin{array}{l}
i\frac{\partial \xi_{3}}{\partial z} +
    i\frac{1}{v_{3g}}\frac{\partial \xi_{3}}{\partial t}
    +k_{32}\xi_{2}-k_{32}l\xi_{3} + N_{3}[(\xi_{3}+\xi_{3}^{\ast})]= 0.
\end{array}
\end{equation}



with: $N_{1}= R_{1}[1-\frac{3}{2}\Gamma Q_{1}-\frac{3}{4}\Gamma^{2}Q_{1}^{2}]$ ; $N_{2}= R_{2}[1-\frac{3}{2}\Gamma Q_{2}-\frac{3}{4}\Gamma^{2}Q_{2}^{2}]$ and $N_{3}= R_{3}[1-\frac{3}{2}\Gamma Q_{3}-\frac{3}{4}\Gamma^{2}Q_{3}^{2}]$

$R_{1} = \frac{P}{1+f^{2}+\frac{f^{2}}{l^{2}}}\gamma_{1}$ ;  $R_{2} = \frac{P}{1+\frac{1}{f^{2}}+\frac{1}{l^{2}}}\gamma_{2}$ and  $R_{3} = \frac{P}{1+l^{2}+\frac{l^{2}}{f^{2}}}\gamma_{3}$ ; $Q_{1} = R_{1}/\gamma_{1}$ ; $Q_{2} = R_{2}/\gamma_{2}$ and $Q_{3} = R_{2}/\gamma_{3}$

 The solution of the above linear equations (\ref{eq12}), (\ref{eq13}) and (\ref{eq14}) can be estimate as follows:
 %In order to solve the set of three linear Eqs. (\ref{eq12}), (\ref{eq13}) and (\ref{eq14}), We
%assume a plane wave ansatz constituting of both forward and
%backward propagation having the form,
\begin{equation}\label{eq15}
\xi_{j}(z,t)=m_{j}e^{i(Kz-\Omega t)}+ n_{j}e^{-i(Kz-\Omega t)}.
\end{equation}

Where $m_{j}$ and $n_{j}$ are real constants, K is the wave number and $\Omega$ is the perturbation frequency. By introducing
Eqs.(\ref{eq15}) into Eqs. (\ref{eq12}), (\ref{eq13}) and (\ref{eq13}), we obtain a matrix $6\times6$ having the following elements:
set of six linear coupled for $m_{i}$ and $n_{i}$.




\begin{equation}\label{eq16}
    \left(%
\begin{array}{cccccc}
  R_{11} & R_{12} & R_{13} & R_{14}& R_{15} & R_{16} \\
  R_{21} & R_{22} & R_{23} & R_{24}& R_{25} & R_{26} \\
  R_{31} & R_{32} & R_{33} & R_{34}& R_{35} & R_{36} \\
  R_{41} & R_{42} & R_{43} & R_{44}& R_{45} & R_{46} \\
  R_{51} & R_{52} & R_{53} & R_{54}& R_{55} & R_{56} \\
  R_{61} & R_{62} & R_{63} & R_{64}& R_{65} & R_{66} \\
\end{array}%
\right)\left(%
\begin{array}{c}
  m_{1} \\
  m_{2} \\
  m_{3} \\
  n_{1} \\
  n_{2} \\
  n_{3} \\
\end{array}%
\right) = 0.
\end{equation}
Where

\begin{equation}\label{eq17}
  R_{11}= N_{1};R_{12}= 0 ; R_{13}= 0 ; R_{14}= K-S-k_{12}f+N_{1} ; R_{15}= k_{12}; R_{16}= 0, \nonumber
\end{equation}
\begin{equation}\label{eq17}
  R_{21}= N_{1}-K-k_{12}f+S; R_{22}= k_{12}; R_{23}=0; R_{24}= N_{1}; R_{25}= 0; R_{26}= 0,   \nonumber
\end{equation}
\begin{equation}\label{eq17}
  R_{31}= 0; R_{32}= N_{2}; R_{33}= 0; R_{34}= k_{21}; R_{35}= N_{2}-K-k_{12}f^{-1}-k_{23}l^{-1}-S; R_{36}= k_{23},
\end{equation}
\begin{equation}\label{eq17}
  R_{41}= k_{21}; R_{42}= N_{2}+K-k_{12}f^{-1}-k_{23}l^{-1}+S; R_{43}= k_{23}; R_{44}= 0; R_{45}= N_{2}; R_{46}=0, \nonumber
\end{equation}
\begin{equation}\label{eq17}
  R_{51}= 0; R_{52}= 0; R_{53}= N_{3}; R_{54}= 0; R_{55}= k_{32}; R_{56}= N_{3}+K-k_{32}l-S,  \nonumber
\end{equation}
\begin{equation}\label{eq17}
  R_{61}=0; R_{62}= k_{32}; R_{63}= N_{3}-K-k_{32}l+S; R_{64}= 0; R_{65}= 0; R_{66}= N_{3}.\nonumber
\end{equation}
And


%This set have a nontrivial solution only if the determinant of the coefficient matrix is zero, which is given below

%However, in general $k_{12} \neq k_{21}\neq k_{23}\neq k_{32}$ and $v_{1g} \neq v_{2g}\neq v_{3g}$,
%but in this paper, in order to highlight the most important new
%physical effects associated with the PIM-NIM nonlinear coupler, in
%the following discussion
We assume in this paper that $ k_{12}=k_{21}=k_{23}=k_{32}=k; v_{1g}=v_{2g}=v_{3g}=v_{g}$  and  $S=\Omega/v_{g}.$

%We will use the matrix $6\times6$ to study the stability of the
%system that we considered. The determinant of the associated matrix in Eq.(\ref{eq16}) can be written as:
%
%\begin{equation}\label{eq20}
%    aS^{6} + bS^{5} + cS^{4} + dS^{3} + eS^{2} + fS+ g =0
%\end{equation}
%
%Where a, b,c, d, e,f and g are defined in Appendix.
%The determinant of the matrix R leads to a six order polynomial in S, where the roots should possess a
%nonzero and negative imaginary part that corresponds
%to a dispersion relation. MI occurs when the wave number possesses a nonzero imaginary part, which corresponds to an exponential growth of the perturbed amplitude. Assuming that roots $S_{i}$ with a negative part exist,
%the instability gain is given by the equation
The roots of the determinant of the matrix R should possess a nonzero and negative imaginary part. This allow to determine the stability of CW.
Then, the discussion of the MI depends on the power gain spectrum given by:

\begin{equation}\label{eq21}
  G = max_{i}[-Im(S_{i})].
\end{equation}

    %G= \left|Im(S_{Max})\right|
%where $Im(S_{Max})$ denotes the imaginary part of $S_{Max}$, and
%note that $S_{Max}$ is the root with the largest value.



%\section{Modulational instability gain spectrum}
\section{MODULATIONAL INSTABILITY IN THREE-CORE PIM-NIM-PIM OPTICAL COUPLER}


%We investigate now the study of MI in three-core coupler with NIM channel in presence of saturation parameter.
 We will discuss in detail now the dynamical behaviors of MI
in the oppositely directed three-core coupler. Six cases will be analyzed here:

\begin{equation}\label{eq21}
  (i)~ f>0 ~;~ l>0 ~and~ \gamma_{1}=\gamma_{2}=\gamma_{3}=1~;~ (iv)~ f<0 ~; ~l<0 ~and~ \gamma_{1}=\gamma_{2}=\gamma_{3} = 1.\nonumber
\end{equation}
\begin{equation}\label{eq17}
 (ii)~ f>0 ~;~ l>0 ~and~ \gamma_{1}=\gamma_{3} =1~;~ \gamma_{2}=0 ~;~(v)~ f<0~ ;~ l<0 ~and~ \gamma_{1}=\gamma_{3} =1~;~ \gamma_{2}=0.  \nonumber
\end{equation}
\begin{equation}\label{eq17}
 (iii)~ f>0 ~;~ l>0 ~and~ \gamma_{1}=\gamma_{3} =0~;~ \gamma_{2}=1~;~(vi)~ f<0 ~; ~l<0 ~and~ \gamma_{1}=\gamma_{3} =0~;~ \gamma_{2}=1.  \nonumber
\end{equation}

%\subsection{Influence of $\emph{f}$ and $\tau$ on MI spectrum}
\subsection{Influence of pump power and Saturable Nonlinearity on Modulational Instability}

We discuss the influence of the power and coupling coefficient on the modulation instability gain before discussing in detail the individual cases listed above. In normal dispersion regime, Fig.1a show us that in absence of saturation we observe the increases of initial gain with K and then is saturated \cite{shafeeque2}. In Fig.1b in presence of saturable parameter $\Gamma=0.2$ the same result are observe where the increase of P increase the MI and a new band appear for high value of K with a lower amplitude. %gain increases for relatively high values of K, reaches a maximum and then decreases.
However in anomalous dispersion regime, Fig.2a present in absence of saturation a single conventional MI band which increase with the increase of input power. Fig.2b show us the case where $\Gamma=0.2$, here we can see the presence of new stability region which appear for the high values of K due to the presence of saturation.
The two instability regions here are separated by a relatively wide stable region with lower pump power. The primary band at lower K values is due to the nonlinear PIM-NIM-PIM channel with $\Gamma=0$, while the second band is a result of Saturation. When power increases, the gain increases proportionally, and the two instability regions approach each other, thereby narrowing the stability region.

\subsection{Influence of coupling coefficient and Saturable Nonlinearity on Modulational Instability}

To understand the Influence of coupling coefficient and saturation on instability gain, Figs.3 and 4 were plotted.
Figs.3 and 4 shows the gain as a function of K for some values of coupling coefficient k. In normal dispersion region show in Fig.3a for $\Gamma=0$, MI gain increases with the increases of coupling coefficients and the instability band shifts towards higher value of K and saturate, but in presence of $\Gamma=0.2$ show by Fig.3b, MI gain increase in K and reaches a maximum and then decreases.
Fig.4 show the instability dependence on the coupling coefficient for anomalous dispersion regime. In Fig.4a where $\Gamma=0$ one single MI gain is observe and increases with K when the coupling coefficient increase. However the case is different when the saturation is taking into account, Fig.4b show the case where $\Gamma=0.2$, here a new band appear when the saturable parameter exist for the higher value of K. We also note here that the MI gain increase with the increase of coupling coefficient parameter. The two instability regions separated by a relatively wide stable region is observed with higher value of coupling coefficient contrary to the one of pump power. As coupling coefficient increases, the gain increases proportionally, but the wide stable region increase also. If the coupling coefficient decreases, the gain would decreases proportionally, and the two instability regions
approach each other, thereby narrowing the stability region.



\begin{figure}
  \begin{center}
  \includegraphics[width=8cm,height=9cm]{figa1.eps}\textbf{(a)}\includegraphics[width=8cm,height=9cm]{figa.eps}\textbf{(b)}\\
  \caption{The MI gain versus wave vector K in the normal dispersion regime for different values of input power and coupling coefficient in presence of saturable parameter $\Gamma$. (a) $\Gamma = 0$ and (b) $\Gamma = 0.2$.}\label{fig1}
  \end{center}
\end{figure}

\begin{figure}
  \begin{center}
  \includegraphics[width=8cm,height=9cm]{figa3.eps}\textbf{(a)}\includegraphics[width=8cm,height=9cm]{figb.eps}\textbf{(b)}\\
  \caption{The MI gain versus wave vector K in the anomalous dispersion regime for different values of input power and coupling coefficient in presence of saturable parameter $\Gamma$. (a) $\Gamma = 0$ and (b)  $\Gamma = 0.2$.}\label{fig1}
  \end{center}
\end{figure}


\begin{figure}
  \begin{center}
  \includegraphics[width=8cm,height=9cm]{figa2.eps}\textbf{(a)}\includegraphics[width=8cm,height=9cm]{figc.eps}\textbf{(b)}\\
  \caption{The MI gain versus wave vector K in the normal dispersion regime for different values of input power and coupling coefficient in presence of saturable parameter $\Gamma$. (a) $\Gamma = 0$ and (b) $\Gamma = 0.2$.}\label{fig1}
  \end{center}
\end{figure}
\begin{figure}
  \begin{center}
  \includegraphics[width=8cm,height=9cm]{figa4.eps}\textbf{(a)}\includegraphics[width=8cm,height=9cm]{figd.eps}\textbf{(b)}\\
  \caption{The MI gain versus wave vector K in the anomalous dispersion regime for different values of input power and coupling coefficient in presence of saturable parameter $\Gamma$. (a) $\Gamma = 0$ and (b)  $\Gamma = 0.2$.}\label{fig1}
  \end{center}
\end{figure}



\subsection{Influence of $\emph{f}$ on Modulational Instability}




We focus on the influence of $\emph{f}$ for different combination values of nonlinear parameters  $\gamma_{1}, \gamma_{2}$ and $\gamma_{3}$  and saturation parameter on MI  of. To this end, we know That $sgn(f)= 1$ stands for normal and $sgn(f)= - 1$ for anomalous dispersion, then both the sign and the value of f can influence MI.



\begin{figure}
  \begin{center}
  \includegraphics[width=7cm,height=7cm]{fig1a.eps}\textbf{(a)}\includegraphics[width=7cm,height=7cm]{fig1.eps}\textbf{(b)}\\
\includegraphics[width=7cm,height=7cm]{fig2a.eps}\textbf{(c)}\includegraphics[width=7cm,height=7cm]{fig2.eps}\textbf{(d)}\\
\includegraphics[width=7cm,height=7cm]{fig3a.eps}\textbf{(e)}\includegraphics[width=7cm,height=7cm]{fig3.eps}\textbf{(f)}\\
  \caption{(Color online) The MI gain spectra in normal dispersion regime vs $\emph{f}$ $(\emph{f}>0)$ and wave vector K under different nonlinear conditions for k=10 $cm^{-1}$ ; P=10$cm^{-1}$ and l = 2.~~(a ~~and~~ b)~~$\gamma_{1}=1,\gamma_{2}=1,\gamma_{3}=1$,  ~~(c~~ and~~ d)~~$\gamma_{1}=1,\gamma_{2}=0,\gamma_{3}=1 $ ~~(e~~and f)~~$\gamma_{1}=0,\Upsilon_{2}=1, \gamma_{3}=0.$ ~~For~~ (a);(c) and (e)~~$\Gamma
=0$ ~~and ~~for ~~ (b); (d) and (f)~~$\Gamma=0.1$.}\label{fig1}
  \end{center}
\end{figure}
%%
%%
\begin{figure}
  \begin{center}
  \includegraphics[width=7cm,height=7cm]{fig6a.eps}\textbf{(a)}\includegraphics[width=7cm,height=7cm]{fig6.eps}\textbf{(b)}\\
\includegraphics[width=7cm,height=7cm]{fig7a.eps}\textbf{(c)}\includegraphics[width=7cm,height=7cm]{fig7.eps}\textbf{(d)}\\
\includegraphics[width=7cm,height=7cm]{fig8a.eps}\textbf{(e)}\includegraphics[width=7cm,height=7cm]{fig8.eps}\textbf{(f)}\\
  \caption{(Color online) The MI gain spectra in normal dispersion regime vs $\emph{f}$ $(\emph{f}<0)$ and wave vector K under different nonlinear conditions for k = 10 $cm^{-1}$ ; P = 10$cm^{-1}$ and l = -2.~~(a and b)~~$\gamma_{1}=1,\gamma_{2}=1,\gamma_{3}=1$ ~~(c and d)~~$\gamma_{1}=1,\gamma_{2}=0,\gamma_{3}=1; $ ~~(e and f)~~$\gamma_{1}=0,\gamma_{2}=1, \gamma_{3}=0$, ~~For~~ (a);(c) and (e)~~$\Gamma=0$ ~~and ~~for ~~ (b); (d) and (f)~~$\Gamma=0.6$.}\label{fig1}
  \end{center}
\end{figure}



%\\newpage


 To analyse the impact of f on MI, we stand firstly in the normal dispersion regime where the parameter $\emph{f} >0$.
 Figure 5 show the dependence of the gain spectrum with respect to K and $\emph{f}$.
Figs(5a) and (5b) describes a coupler of case (i) where all the channels are nonlinear. In absence of saturable parameter we obtained four sidebands in Fig. 5(a)\cite{shafeeque2}. However, in presence of saturable parameter $\Gamma = 0.1$ Fig.5(b) show us more than four sidebands, the threshold condition also exist but not in all cases. It is clear here that the instability bands under saturation effects reduce the width  also its amplitude.
Figs.5(c) and 5(d) shows the case in which a NIM channels is linear ($\gamma_{2} = 0$) stand for case (ii). In Fig.5(c) two centered sidebands are observed around the zero propagation constant region and a nil value gain exists along the zero propagation constant K=0. Fig.5(d) show us the case where ($\Gamma = 0.1$) and there two sideband were observed. We also note here that by increasing the saturation parameter the gain amplitude is reduced also its width.
In Figs.5(e) and 5(f) we study the effect of NIM channel, here Channels 1 and 3 are linear and channel 2 nonlinear stand for case (iii) in presence and in absence of saturation. Note that without saturation ($\Gamma = 0$) there are two MI bands centered around the zero propagation constant region which are observed in Fig.5(e). In presence of saturation Fig.5(f) it is evident the influence of the value of the saturation parameter on the number of MI bands.
Then, by comparing the Fig.5(e) with Fig.5(f) one observes the appearance of a new instability bands for the highest value of $\emph{f}$ in Fig.5(f) due to the presence of saturation . However, the increase of saturation parameter reduced the maximum instability gain and its width.
From Fig.5, we can concluded that adjusting the value of $\emph{f}$, and $\Gamma$ can be used to control MI in nonlinear three core coupler. In normal dispersion regime threshold condition exists for $\emph{f}$, further the increase of saturation parameter can reduced the maximum gain and also reduce the width of sideband.


 In anomalous dispersion regime, MI can also be altered by varying $\emph{f}$ such as ($\emph{f} < 0$).
In Fig.6(a) and (b) standing for the case where the three media are nonlinear media (case (iv));  we note in Fig.6(a) that in absence of saturation parameter ($\Gamma = 0$) there exist two MI bands centered around the line K = 0 \cite{shafeeque2}. However, in presence of saturation parameter for ($\Gamma = 0.6$) we observe in Fig.6(b) four bands, in this figure the influence of saturable parameter is clearly observe on the number of MI bands.
Figs.6(c) and 6(d) When channels 1 and 3 standing for PIM channels are nonlinear and channel 2 standing for NIM channel is linear $\gamma_{2} = 0$ show the gain spectrum stand for case (v), two distinct sidebands appear on either side of zero propagation constant region. In Fig.6(c), in absence of saturation two sidebands are obtained and the maximum gain attained at lower value of $\left|f\right|$.
In Fig.6(d) the maximum gain and the band width are influenced by the presence of saturation parameter, when saturation is present MI only presents within a limited range of f and MI is obtained for lower value of $\left|f\right|$. Threshold condition for f can exist only in Fig.6(d) instead of the Fig.6(c).
In Figs.6(e) and 6(f) where channels 1 and 3 are linear ($\gamma_{1} = \gamma_{3} = 0$) and channel 2 nonlinear $\gamma_{2} = 1$ (case vi) two sidebands are obtained for the lower value of $\left|f\right|$. In presence of saturation, Fig.6(f) also show us that the maximum gain increase with the increase of saturable parameter and a large stability zone around the propagation constant K exist.
 Hence, from Fig.6 it can be concluded that in anomalous dispersion regime saturation enlarge the gain value and the generation region of the MI, and the MI generation is thresholdless for $\emph{f}$.


%Firstly we focus on the influence of s effect on MI in oppositely directed coupler for different


\section{Direct Numerical Simulations}

In this section, to verify the modulation stability/instability of weakly perturbed continuous waves studied in the analytical form, the nonlinear evolution of the MI was performed numerically. Pseudo-spectral method \cite{gottlieb,fornberg} was used to employ direct numerical simulations of Eqs.(1), (2) and (3). The initial conditions were taken in the form of the wave plane with imposed small periodic perturbation:

\begin{equation}\label{eqa1}
  a_{j}(0,\tau) = u_{j} + \varepsilon cos(\omega_{0}\tau),   (j= 1,2,3).
\end{equation}

Where $\varepsilon = 5*10^{-2}$ a relative value of the perturbation, and a frequency $\omega_{0}$.
%The parameters of numerical simulations are those obtained from the study of the linear stability analysis.
%The propagation of perturbed CW solutions for the case where the saturable parameter are taking into account are displayed in Figs.(7) to (10).

 Figs.(7) to (10) show the results for normal and anomalous dispersion regime, we set the parameters as  $k=10$ $m^{-1}$, $P=10$ kW, $\left|f\right| = \left|L\right| = 2$.

  Figs.(7) and (8) show us the case where the three channels are nonlinear  $\gamma_{1}= \gamma_{2}= \gamma_{3} = 1 ~~ (kW m)^{-1}$ (case i and iv). Fig.7 stand for the evolution MI in the normal dispersion regime and Fig.8 anomalous dispersion regime.
 We observe in Fig.(7a and 7c) in absence of saturation $\Gamma = 0$ a periodic chain of solitons like pulses are produced in both cores.
 Figs.(7b )and (7d) shows the influence of saturation parameter $\Gamma = 0.03$. In this case, main effects are oscillations of the background. In anomalous dispersion regime in presence of saturation show in Figs. (8b) and (8d) we clearly observe the stability at the initial stage of the evolution of initial perturbed CW, but the train of solitons turns into a chaotic pulse after certain propagations distance.

Fig.(9) show the case where the system is only influence by the PIM channels stand for case (ii) and (v). In absence of saturation Figs.(9a) and (9c) show us a generation of a periodic array of peaks with growing amplitude. However when saturation parameter is present $\Gamma = 0.2$; Figs (9b) and (9d) show there a chain of solitons with growing amplitudes are generated in normal and anomalous dispersion regime.

 Fig.(10) reveals the impact of NIM channel ($\gamma_{1}=0 (kW m)^{-1}$, $\gamma_{2}=1 ~~ (kW m)^{-1}$,~~ $\gamma_{3}=0$) and saturable parameter on MI stand for case (iii) et (vi). In absence of saturation Figs.10a and 10c show this case, the MI generates a chain of growing peaks narrow soliton oscillating with higher amplitude for high value  of propagation distance.
 In Figs.(10b) and (10d), we address the case when saturation $\Gamma = 0.2$ and NIM channel $\gamma_{1} = 1$ in normal and anomalous dispersion regime. In this case, we clearly observe the stability at the initial stage of the evolution of initial perturbed continuous wave , but after certain propagations distance a periodic array of peaks with growing amplitudes is generated. We clearly observe from Fig.(10) that the presence of NIM channel and saturation influence the generation of MI.


\begin{figure}
  \begin{center}
  \includegraphics[width=7cm,height=7cm]{fign1.eps}\textbf{(a)}\includegraphics[width=7cm,height=7cm]{fign3.eps}\textbf{(b)}\\
\includegraphics[width=7cm,height=7cm]{fign2.eps}\textbf{(c)}\includegraphics[width=7cm,height=7cm]{fign4.eps}\textbf{(d)}\\
  \caption{(Color online) Evolution of power wave showing the effects of saturation on the development of MI in normal dispersion regime under different nonlinear conditions ~~ $\gamma_{1}=1  ~~(kW m)^{-1}$, $\gamma_{2}=1 ~~ (kW m)^{-1}$,~~ $\gamma_{3}=1 ~~ (kW m)^{-1}$. For (a) and (c)~~$\Gamma =0~ kW^{-1}$, for (b) and (d)~~$\Gamma =0.03~ kW^{-1}$.}\label{fig1}
  \end{center}
\end{figure}


%f = 2, L = 2, $k=10$ $m^{-1}$, $P=10$ kW,~~


\begin{figure}
  \begin{center}
  \includegraphics[width=7cm,height=7cm]{fign1a.eps}\textbf{(a)}\includegraphics[width=7cm,height=7cm]{fign3a.eps}\textbf{(b)}\\
\includegraphics[width=7cm,height=7cm]{fign2a.eps}\textbf{(c)}\includegraphics[width=7cm,height=7cm]{fign4a.eps}\textbf{(d)}\\
  \caption{(Color online) Evolution of power wave showing the effects of saturation on the development of MI in anomalous dispersion regime under different nonlinear conditions ~~ $\gamma_{1}=1  ~~(kW m)^{-1}$, $\gamma_{2}=1 ~~ (kW m)^{-1}$,~~ $\gamma_{3}=1 ~~ (kW m)^{-1}$. For (a) and (c)~~$\Gamma =0~ kW^{-1}$, for (b) and (d)~~$\Gamma =0.03~ kW^{-1}$.}\label{fig1}
  \end{center}
\end{figure}
%. Also, the MI for the anomalous case of group
%velocity dispersion is illustrated in Fig. 2
%\includegraphics[width=7cm,height=7cm]{fign5.eps}\textbf{(a)}\includegraphics[width=7cm,height=7cm]{fign7.eps}\textbf{(b)}\\
\begin{figure}
  \begin{center}
  \includegraphics[width=7cm,height=7cm]{fign9.eps}\textbf{(a)}\includegraphics[width=7cm,height=7cm]{fign10.eps}\textbf{(b)}\\
\includegraphics[width=7cm,height=7cm]{fign9a.eps}\textbf{(c)}\includegraphics[width=7cm,height=7cm]{fign10a.eps}\textbf{(d)}\\
  \caption{(Color online) Evolution of power wave showing the effects of saturation on the development of MI under different nonlinear conditions ~~ $\gamma_{1}=1  ~~(kW m)^{-1}$, $\gamma_{2}=0 ~~ (kW m)^{-1}$,~~ $\gamma_{3}=1 ~~ (kW m)^{-1}$.(a and b) normal dispersion regime, (c and d) anomalous dispersion regime.  For (a) and (c)~~$\Gamma =0~ kW^{-1}$, for (b) and (d)~~$\Gamma =0.2~ kW^{-1}$.}\label{fig1}
  \end{center}
\end{figure}




\begin{figure}
  \begin{center}
\includegraphics[width=7cm,height=7cm]{fign6.eps}\textbf{(a)}\includegraphics[width=7cm,height=7cm]{fign8.eps}\textbf{(b)}\\
\includegraphics[width=7cm,height=7cm]{fign6a.eps}\textbf{(c)}\includegraphics[width=7cm,height=7cm]{fign8a.eps}\textbf{(d)}\\
  \caption{(Color online) Evolution of power wave showing the effects of saturation on the development of MI under different nonlinear conditions ~~ $\gamma_{1}=0  ~~(kW m)^{-1}$, $\gamma_{2}=1 ~~ (kW m)^{-1}$,~~ $\gamma_{3}=0 ~~ (kW m)^{-1}$. (a and b) normal dispersion regime, (c and d) anomalous dispersion regime.  For (a) and (c)~~$\Gamma =0~ kW^{-1}$, for (b) and (d)~~$\Gamma =0.2~ kW^{-1}$.}\label{fig1}
  \end{center}
\end{figure}

%\begin{figure}
%  \begin{center}
%\includegraphics[width=7cm,height=7cm]{fign6a.eps}\textbf{(a)}\includegraphics[width=7cm,height=7cm]{fign8a.eps}\textbf{(b)}\\
%  \caption{(Color online) Evolution of power wave showing the effects of saturation on the development of MI in anomalous dispersion regime under different nonlinear conditions ~~f =-2, L =-2, $k=10$ $m^{-1}$, $P=10$ kW,~~, $\gamma_{1}=0  ~~(kW m)^{-1}$, $\gamma_{2}=1 ~~ (kW m)^{-1}$,~~ $\gamma_{3}=0 ~~ (kW m)^{-1}$. For (a)~~$\Gamma =0~ kW^{-1}$ and for (b)~~$\Gamma =0.2~ kW^{-1}$.}\label{fig1}
%  \end{center}
%\end{figure}




\newpage

\section{Conclusion}
This paper investigated more the MI with modified Kerr-type saturable nonlinearity(MSN) in nonlinear three-core directional couplers.
Continuous wave solutions and linearizing firtly are used to obtain the dispersion relation. Secondly  the influence on MI of pump power, coupling coefficient in presence of saturation was analysed. The result show that, the presence of saturable nonlinearity can bring new sidebands or reduce the number of sidebands through the systems. However, the presence of quintic and septic nonlinearity bring by saturation (MSN) can also enlarge the generation of MI.
At end, a numerical simulation of the nonlinear development of the MI in different regimes which were studied analytically was performed. In absence of saturation, a generation of periodic soliton arrays with a growing amplitude are obtained . However in presence of nonlinear saturation, we clearly observe the stability at the initial stage of the evolution of initial perturbed CW, but the train of solitons turns into a turbulent state after certain propagations distance.
The present study provide a new way to generate solitons or ultrashort pulses in an oppositely directed three-core coupler with saturable nonlinearity.

%reinforces those results presented in \cite{shafeeque2},

%\section*{Acknowledgements}
%This work has been finalized during the stay of AM at the Abdus Salam International Center for Theoretical Physics (ICTP)
%Trieste-Italy through the Associate Program.

%\section*{Declaration of competing interest}
%
%The authors declare that they have no known competing financial interests or personal relationships that could have appeared to
%influence the work reported in this paper.

\newpage
\begin{thebibliography}{99}

\bibitem{veselago}  Veselago, V.G., "The electrodynamics of substances with simultaneously negative values of $\varepsilon$ and $\mu$." Sov. Phys. Uspekhi. 10(4), 509-514(1968). \label{veselago}
\bibitem{shelby}  Shelby, R.A., Smith, D.R., and  Schultz, S. "Experimental verification of a negative index of refraction," Science292,77-79 (2001). \label{shelby}
\bibitem{linden}  Linden, S., Enkrich, C. , Wegener, M., Zhou, J., Koschny, T., and  Soukoulis, C.M. "Magnetic response of metamaterials at 100 terahertz," Science306, 1351–1353 (2004). \label{linden}
\bibitem{zhang1}  Zhang, S., Fan, W., Panoiu, N.C., Malloy, K.J., Osgood, R.M., and Brueck, S.R. "Experimental demonstration of near-infrared negative-index metamaterials," Phys. Rev. Lett.95, 137404 (2005).\label{zhang1}
\bibitem{said}  Said Z., Sihvola, A., Vinogradov, A. "Metamaterials and plasmonics : Fundamentals, Modelling applications, New York : Springer-verlag. Pp. 3-10 Chap 3, 106. \label{said}
\bibitem{boratay} Boratay, A.K., Ekmel, O. "Radiation properties of a split ring resonator and monopole composite" Physica status solidi 224 (4) : 1192-96 (2007).\label{boratay}
\bibitem{jensen}  Jensen, S.M. "The nonlinear coherent coupler," IEEE. J. Quantum Electron. 18, 1580-1583 (1982). \label{jensen}
\bibitem{kengne} Kengne, E., Liu, WM.: Eur. Phys. J. B, 92, 235 (2019).\label{kengne}
 \bibitem{litchinitser}Litchinitser, N.M., Gabitov, I.R., Maimistov, A.I.: "Optical Bistability in a Nonlinear Optical Coupler with a Negative Index Channel" Phys. Rev. Lett. 99, 113902 (2007).\label{litchinitser}
 \bibitem{halterman}Halterman, K., Elson, J.M., Overfelt, P.L.: "Characteristics of boundmodes in coupled dielectric waveguides containing negative index media," Opt. Express 11, 521-529 (2003).\label{halterman}
  \bibitem{yuan} Yuan, Y., Ran, L., Chen, H., Huangfu, J., Grzegorczyk, T.M., Kong, J.A.: "Backward coupling waveguide coupler using left-handed material," Appl. Phys. Lett. 88, 211903 (2006).\label{yuan}
\bibitem{agrawal} Agrawal, G.P.: Nonlinear Fiber Optics, third ed., Academic Press Inc., San Diego, (2001). \label{agrawal}
\bibitem{xiang}  Xiang, Y.,  Wen, S.,  Dai, X., Fan, D.: "Modulation instability in nonlinear oppositely directed coupler with a negative-index metamaterial channel," Phys. Rev. E \textbf{82}, 056605 (2010).\label{xiang}

\bibitem{shafeeque1}  Shafeeque, A.K., Porsezian, K., Uthayakumar, T.: "Influence of self-steepening and intrapulse Raman scattering on modulation instability in oppositely directed coupler," Phys. Rev. E. \textbf{90}, 042910 (2014).\label{shafeeque1}
\bibitem{zhang} Zhang,  J.,  Dai, X., Zhang, L., Xiang,  Y., Li,  Y.: "Modulation instability in the oppositely directed coupler with a quadratic nonlinearity," \textbf{32}, 1-8 (2015).\label{zhang}
\bibitem{shafeeque4}  Shafeeque, A.K.,  Nithyanandan, K.,  Porsezian, K., Maimistov, A.I.: "Influence of birefringence in the instability spectra of oppositely directed coupler with negative index material channel," Phys. Rev. A. \textbf{93}, 023848 - 023858 (2016).\label{shafeeque4}
\bibitem{kengne2}  Kengne, E.: " Modulational instability and soliton propagation in an alternate right-handed and left-handed multi-coupled nonlinear dissipative transmission network, Chaos, Soliton and Fractals. \textbf{146}, 110866 (2021).\label{kengne2}
\bibitem{tatsing}  Tatsing, P.H., Mohamadou, A.,  Bouri, C., Tiofack, C.G.L.,  Kofane, T.C.: "Modulation instability in nonlinear positive-negative index couplers with saturable nonlinearity," J. Opt. Soc. Am. B \textbf{29}, 3218-3225 (2012).\label{tatsing}
\bibitem{alves}  Alves, E.O.,  Cardoso, W.B.,  Avelar, A.T.: "Modulation instability in a nonlinear oppositely directed coupler with saturable nonlinearities and higher-order effects," J. Opt. Soc. Am. B \textbf{33}, 1134-1142 (2016).\label{alves}
\bibitem{tatsing2}  Tatsing, P.H., Mohamadou, A., and   Tiofack, C.G.L.,: "Modified  Kerr-type  saturable  nonlinearity  effect  on  the  modulational instability  of  nonlinear  coupler  with  a  negative-index  metamaterial channel," Optik \textbf{127}, 4150-4155 (2016).\label{tatsing2}
\bibitem{tatsing3}  Tatsing,P.H., Mohamadou, A., Tiofack, C.G.L., T.C. Kofane,  " Modulational instability in nonlinear oppositely directed coupler with saturable delayed nonlinear response,"  Optical and  Quantum  Electronics \textbf{1493}, 1-20 (2018).\label{tatsing3}
\bibitem{tatsing4}  Tatsing, P.H.,  Mohamadou, A.,   Tiofack, C.G.L.,  Kofane, T.C.:  " Modulational instability in double-doped   directed couplers with  a non-kerr-like nonlinear refractive index change," Optical and  Quantum  Electronics \textbf{390}, 1-21 (2019).\label{tatsing4}
\bibitem{houwe}  Houwe, A.,  Abbagari, S.,  Saliou, Y., Mustafa Inc,  Doka, S.Y.,  Bouetou, T.B.,  Bayram, M.:  " Attitude of the Modulation Instability gain in Oppositely Directed Coupler with the effects of the Intrapulse Raman Scattering and Saturable Function," Results in Physics \textbf{31}, 104851 (2021).\label{houwe}
\bibitem{langridge}  Langridge, P.E., Firth, W.J.: Opt. Quantum Electron. 24315 (1992). \label{langridge}
\bibitem{shafeeque2}  Shafeeque, A.K.,  Nihyanandan, K., Porsezian, K.: "Theoretical investigation of modulation instability in a three-core coupler with negative index material channel," Phys. Lett. A. \textbf{11}, 223-229 (2015).\label{shafeeque2}
\bibitem{shafeeque3}  Shafeeque, A.K., Nihyanandan, K., Porsezian, K.,  Maimistov, A.I.:  "Modulation instability in a triangular threecore coupler with a negative-index material channel," J. Opt. \textbf{18}, 035102-035111 (2016).\label{shafeeque3}

\bibitem{xianqiong} X. Zhong, Ke Cheng, and K. S. Chiang, J. Opt. Soc. Am. B \textbf{31}, (2014) 7. \label{xianqiong}
\bibitem{litsinistser1} N. M. Litchinitser, C. J. Mc Kinstrie, C. M. de Sterke and G. P. Agrawal, "Instabilities and solitons in systems with spatiotemporal dispersion," J. Opt. Soc. Am B \textbf{18}, (2001) 45-54.\label{litsinistser1}
\bibitem{gottlieb}  Gottlieb, D.,  Orszag, S.A.,: Numerical Analysis of Spectral Methods (Society for Industrial and Applied Mathematics, (1977).\label{gottlieb}
\bibitem{fornberg}  Fornberg, B.: A Practical Guide to Pseudospectral Methods (Cambridge University, 1984).\label{fornberg}

%    \bibitem{kengne1}  Kengne, E.: " Envelope modulational instability in nonlinear dissipative transmission line. J. Nonlinear Oscil. \textbf{5(1)}, 23-31 (2002).\label{kengne1}
%\bibitem{kengne3}     Kengne, E., Liu, WM.: Eur. Phys. J. B, \textbf{92}, 235 (2019). \label{kengne3}
\end{thebibliography}
\end{document}