\documentclass[11pt,reqno]{amsart}
\usepackage[cp1251]{inputenc}
\usepackage[english,russian]{babel}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{graphicx}
\usepackage{xcolor}
\usepackage[colorlinks]{hyperref}



\usepackage{titlesec}
\titlespacing{\section}{0pt}{2.5ex}{1.5ex}
\titlespacing{\subsection}{0pt}{1.5ex}{1ex}
\titlespacing{\subsubsection}{0pt}{1.5ex}{1ex}
\titleformat{\section}{\large\bfseries\centering}{\thesection}{1em}{}
\titleformat{\subsection}[runin]{\bfseries}{\thesubsection.}{0.5em}{}[.\mbox{\ }]
\titleformat{\subsubsection}[runin]{\bfseries}{\thesubsubsection.}{0.4em}{}[.\mbox{\ }]





%%%%%---------- ORCiD
\newcommand\orcidicon[1]{\href{https://orcid.org/#1}{\includegraphics[scale=0.02]{orcid.pdf}}}


\def\udcs{532.135} %Here the author places classificators of the paper according to Russian classification system
\def\mscs{76А10} %Here the author places  classificators according to the AMS classification list.
\setcounter{page}{144}

\newtheorem{lemma}{Lemma}
\newtheorem{theorem}{Theorem}
\newtheorem{definition}{Definition}
\newtheorem{corollary}{Corollary}
\newtheorem{proposition}{Proposition}



\hypersetup{linkcolor=blue,filecolor=black,urlcolor=blue, citecolor=blue}

\def\logo{\textcolor{blue}{{\bf\Huge S\raisebox{0.2ex}{\hspace{0.55ex}\raisebox{0.05ex}e\hspace{-1.65ex}$\bigcirc$}MR}}}

\def\semrtop
     {
  \vbox{
     \noindent\logo \hspace{-3mm}  \parbox{12cm}{     \begin{center}
     {\LARGE \textsf{СИБИРСКИЕ ЭЛЕКТРОННЫЕ}} \\[2mm]
     {\LARGE \textsf{МАТЕМАТИЧЕСКИЕ ИЗВЕСТИЯ}} \\[2mm]
     {\Large \textsf{Siberian Electronic Mathematical Reports}} \\[1mm]
     {\Large\tt{\textcolor{blue}{http://semr.math.nsc.ru}}}\\[1mm]
     \raisebox{1ex}{ISSN 1813-3304 }
     \vspace*{-3mm}
%     {\small Институт математики им. С.Л.Соболева СО РАН}\\[-1mm]
%     {\small Sobolev Institute of Mathematics SB RAS}
     \end{center}
}
     \vspace{-1mm}
     \noindent
      \begin{tabular}{c}
     \hphantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa!} \\
     \hline\hline
     \end{tabular}

     \vspace{3mm}
     {\flushleft\it Vol. ??, No. ?, pp. ??? (2025) \hspace{\fill}{\rm\small УДК \udcs \hspace{-4mm}}}
     \newline  %Volume, pages in Russian, filled by Editorial board
     	{\rm\small \textcolor{blue}{https://doi.org/10.33048/semi.2023.20.???}}\hphantom{aaaaaaaaaaaaaaaaaaaaaaaa\;}{\rm\small MSC\ \ \mscs }
  }%\vbox
}




\begin{document}
\renewcommand{\refname}{References}
\renewcommand{\proofname}{Proof.}
\renewcommand{\figurename}{Fig.}




\thispagestyle{empty}

\title[Stability of the flow of viscoelastic liquid]{\large On the stability of the plane flow of viscoelastic polymer liquid}


\author[R.E. SEMENKO]{{\bf R.E. Semenko}\protect\orcidicon{0000-0001-8547-5627}}



\address{Roman Evgenievich Semenko  % Контактные данные всех авторов и место работы указываются только на английском языке
\newline\hphantom{iii} Sobolev Institute of Mathematics,
\newline\hphantom{iii} pr. Koptyuga, 4,
\newline\hphantom{iii} 630090, Novosibirsk, Russia}%
\email{\textcolor{blue}{r.e.semenko@math.nsc.ru}}%

\author[G.N.SHUKUROV]{{\bf G.N. Shukurov}\\ \\
}%
\address{Giyosiddin N. Shukurov  % Контактные данные всех авторов и место работы указываются только на английском языке
\newline\hphantom{iii} Novosibrsk State University,
\newline\hphantom{iii} Pirogova st., 1,
\newline\hphantom{iii} 630090, Novosibirsk, Russia}%
\email{\textcolor{blue}{g.shukurov@g.nsu.ru}}%

\thanks{\sc Semenko R.E., Shukurov G.N.,
On the stability of the plane flow of viscoelastic polymer liquid}
\thanks{\copyright \ 2025 Semenko R.E., Shukurov G.N.}
\thanks{\rm The study of Semenko R.E. was carried out within the framework of the state contract of the Sobolev Institute of Mathematics (project no. FWNF-2022-0008)}
\thanks{\it Received, Published}%


\semrtop \vspace{1cm}
\maketitle {\small

\vspace{-12pt}
\centerline{{\it Communicated by} {\sc ?.?. ?? }
%({\it filled by Editorial board})
}

\bigskip
\bigskip

\begin{quote}
\noindent{\bf Abstract:} We have studied the spectral problem for plane Poiseuille-type flow
of viscoelastic polymer liquid. The flow was modeled with the equations
of rheological Vinogradov--Pokrovskii model. The numerical spectral
procedure was used to calculate eigenvalues of the problem and to
determine the critical values of parameters where the instability
of the flow occurs. It was shown that the critical Reynolds number
goes to well-known value of approximately 5772 (that is the critical
value for the viscous fluid) while the relaxation time goes to zero.
The slight increase of elastic properties of the fluid (or Weissenberg
number) leads to rapid increase of critical Reynolds number, while
further increase of Weissenberg number destabilizes the flow by pushing
the critical Reynolds number to zero. Similar to a number of other
rheological models, the Vinogradov--Porkovskii model demonstrates
the elastic instability effect, that is the spectral problem has unstable
modes that are not a continuation of unstable modes of the viscous
newtonean flow.\medskip

\noindent{\bf Keywords:} Vinogradov--Pokrovskii rheological model, Poiseuille flow, spectral problem.
 \end{quote}
}

\bigskip

\section{Introduction}
Polymer solutions and melts are viscoelastic fluids, and the viscoelasticity
causes the number of features not typical for viscous newtonean fluids.
Among these features are shear-thinning effects, growth of extensional
viscosity, stress relaxation in non-stationary flows and many others
\cite{DoiEdw,Bird}. One field of the particular interest within the
theory of viscoelastic fluids is the stability properties of the flows
of certain geometries. The recent experimental studies of polymer
solutions \cite{Samanta,Srinivas} showed that the addition of small
amounts of polymer in the solution causes complex influence on the
stability of the flows which include both transition delay to turbulence
(the transition happens with higher Reynolds numbers compared to pure
viscous flows), but also the early turbulence effect which on the contrary reduces the critical Reynolds number.

Among the rheological models of polymer dynamics probably the most well-researched one is the famous Oldroyd-B model \cite{Oldroyd}.
Over the years the extensive research was performed for it 
and its simplified forms like upper-convected Maxwell model (UCM)
\cite{Shafeq,Sanchez,Sureshkumar,Larson}. But Oldroyd-B model by
itself is not without flaws due to relatively simple linear constitutive
equations, which leads to shear-rate independent viscosity, unbounded
growth of extensional viscosity and several other properties which
does not fit well with known experimental data \cite{Ramachandran}.
There are large number of more advanced models which addresses these flaws but the knowledge of their properties is much more limited.

One way to improve the adequacy of the modeling of viscoelastic flows
is mesoscopic approach which takes into account both general experimental
observations and dynamics of separate long macromolecules of polymer
\cite{Pokrovskii,AltGusPysh}. The relatively recent modified Vinogradov--Pokrovskii
model (mVP) already demonstrated good agreement with experiment in
several visometric flows \cite{Makarova,Pysh2020,Koshelev,Merzlikina},
but the understanding of the stability properties for this model is
quite limited so far, partially due to mathematical complexity of
the governing equations of the model. Saying that, some existing results
indicate that even basic flows within this model could be notoriously
unstable. There are evidences that the stationary flows in the plane
channels and tubes and even the state of the rests could be linearly
unstable, which put a question regarding the adequacy of the model
\cite{BlokhYegTk,TkachBib,Tkachev2023}. On the other hand, numerical
simulations performed in \cite{Semenko2023} for plane Poiseuille-type
flow in mVP model showed good stabilization of non-stationary solutions
on equilibrium states over time at least for relatively low Reynolds
and Weissenberg numbers. Granted, that simulations were taken for
simplified problem formulation with fixed pressure gradient and transverse
velocity, but nonetheless this evidence promote more thorough study
of linear stability of the problem.

One way to check the adequacy of stability research is to study how
the equations behave for low Weissenberg numbers. Since the zero relaxation
time (and thus zero Weissenberg number) should transform the spectral
problem to the one for pure viscous flow, the results should also
correlate with known information regarding the stability of Poiseuille
flow for Navier--Stokes equations, including the critical Reynolds
number of approximately $\mathrm{Re}\approx$5772. In the present
research we show the numerical results of the stability problem, including
the limit case of low Weissenberg number. We discuss the nature of
discovered instabilities and compare the results with the ones for viscous
Poiseuille flow.

\section{Problem statement}

Let us introduce the problem for linear stability of the steady plane
flow between two parallel plates. Assuming vector $\boldsymbol{u}=(u,v)^{\top}$
to be the velocity of the flow in two-dimensional space $(x,y)^{\top}$
and $2l$ to be a width of the flat channel, consider the flow shown
on Fig. \ref{fig:Scheme}, driven by the pressure gradient along the
$x$-axis. The plates $x=\pm l$ are fixed and no-slip conditions
for the flow are imposed at them. The continuity equation and momentum
equations for incompressible fluid are naturally 
\begin{align*}
 & \mathrm{div}\boldsymbol{u}=u_{x}+v_{y}=0,\\
 & \rho\frac{\mathrm{d}\boldsymbol{u}}{\mathrm{d}t}=\mathrm{div}\sigma,
\end{align*}
here $\rho$ is the constant density and $\sigma$ is the viscoelastic
stress tensor. Also
\[
\frac{\mathrm{d}}{\mathrm{d}t}=\frac{\partial}{\partial t}+(\boldsymbol{u},\nabla).
\]
\begin{figure}
\noindent \centering{}\includegraphics{Figures/Channel_nd}\caption{Flow in plane channel\label{fig:Scheme}}
\end{figure}
The modified Vinogradov--Pokrovskii model (mVP) utilizes the mesoscopic
approach for the formulation of governing equations. The stress tensor
is split into the isotropic and anisotropic parts 
\[
\sigma=-pI+\Pi,
\]
where $p$ is the pressure, $I$ is the unity tensor and $\Pi=\{a_{jk}\},\;j,k=1,2$
is the tensor of anisotropy (tensor of additional stresses) \cite{Pokrovskii}.
The tensor $\Pi$ represents the effects of viscous and elastic forces
of the flow. 

We will use non-dimensional form of the governing equations with the
scales $l$ for length, $u_{H}$ for velocity and $l/u_{H}$ for time,
where $u_{H}$ is the velocity of the steady flow at the center line
of the channel. The pressure is scaled with $\rho u_{H}^{2}$, and
the anisotropic stresses are scaled with $\eta_{0}u_{H}/l$, where
$\eta_{0}$ is the initial viscosity \cite{Pokrovskii}. The equations of
mVP then have the form 
\begin{align}
 & \mathrm{div}\boldsymbol{u}=u_{x}+v_{y}=0,\label{continuity_eq}\\
 & \frac{\mathrm{d}\boldsymbol{u}}{\mathrm{d}t}+\nabla p=\mathrm{div}\Pi,\label{momentum_eq}\\
 & \frac{\mathrm{d}a_{11}}{\mathrm{d}t}-2(a_{11}+1/(\mathrm{ReW}))u_{x}-2a_{12}u_{y}+\mathcal{L}_{11}=0,\label{a11_eq}\\
 & \frac{\mathrm{d}a_{22}}{\mathrm{d}t}-2(a_{22}+1/(\mathrm{ReW}))v_{y}-2a_{12}v_{x}+\mathcal{L}_{22}=0,\label{a22_eq}\\
 & \frac{\mathrm{d}a_{12}}{\mathrm{d}t}-(a_{11}+1/(\mathrm{ReW}))v_{x}-(a_{22}+1/(\mathrm{ReW}))u_{y}+\mathcal{L}_{12}=0.\label{a12_eq}
\end{align}
Here the equations (\ref{a11_eq})-(\ref{a12_eq}) are differential
constitution equations, which connect the stress tensor with tensor
of velocity gradient. $\mathrm{Re}=\rho u_{H}l/\eta_{0}$ is the Reynolds
number, $\mathrm{W}=\tau_{0}u_{H}/l$ is the Weissenberg number, $\tau_{0}$
is the relaxation time. The non-linear terms of (\ref{a11_eq})-(\ref{a12_eq})
are
\[
\mathcal{L}_{jj}=(1/(\mathrm{W})+(k-\beta)\mathrm{Re}(a_{11}+a_{22})/3)a_{jj}+\beta\mathrm{Re}(a_{jj}^{2}+a_{12}^{2}),\;j=1,2,
\]
\[
\mathcal{L}_{12}=(1/(\mathrm{W})+(k+\beta)\mathrm{Re}(a_{11}+a_{22})/3)a_{12}.
\]
The values $k$ and $\beta$ are phenomenological parameters of the
mVP. Thy represents the impact of the size and the orientation of
polymeric macromolecules on the anisotropy of the flow respectively
\cite{Pokrovskii,AltGusPysh}. According to \cite{Golovicheva}, the sensible
ratio is $k=1.2\beta$ which fits well with experimental measurements.
The equations (\ref{continuity_eq})-(\ref{a12_eq}) are supplemented
with boundary conditions 
\begin{equation}
\boldsymbol{u}(\pm1)=0.\label{noslip_cond}
\end{equation}

The stability analysis first require the steady-state solution for
linearization of the problem. Assume that 
\begin{equation}
u=u_{0}(y),\;v=0,\;a_{jk}=a_{jk0}(y),\;j,k=1,2,\label{steady_sol}
\end{equation}
then from (\ref{momentum_eq}) we have 
\[
p=p_{0}(x,y)=P(y)-Ax,
\]
where $P(y)$ is an unknown function, $A$ is the constant pressure
gradient.

Following \cite{Semenko2023,BlokhBamb}, we have 
\[
a_{120}(y)=-Ay,
\]
\begin{multline}
\left(\frac{2\overline{k}}{3}+\beta\right)a_{220}^{3}+\left(1+\beta+\frac{2\overline{k}}{3}\right)a_{220}^{2}/(\mathrm{ReW})+\\
+\left((1/(\mathrm{ReW})^{2}+\frac{2\overline{k}}{3}a_{120}^{2}+(\beta/(\mathrm{ReW})a_{120}^{2}\right)a_{220}+\\+(\beta/(\mathrm{ReW}))a_{120}^{2}=0.
\label{a22_cubic}
\end{multline}
The (\ref{a22_cubic}) provides the values of $a_{220}(y)$ as the
roots of cubic equation with variable coefficients. It was shown in
\cite{Semenko2023} that only one root of this equation could lead
to potentially stable solution, and that root satisfies the condition
$a_{220}(0)=0$. Sticking with this solution, the $a_{110}(y)$ is
found from
\[
a_{110}=a_{220}+\frac{2a_{120}^{2}}{1/(\mathrm{ReW)}+a_{220}}.
\]
 Finally, $u_{0}(y)$ is calculated from 
\begin{equation}
u_{0}'=\frac{\hat{K}_{I}a_{120}}{a_{220}+1/(\mathrm{ReW})}\label{u0_eq}
\end{equation}
with boundary condition $u_{0}(-1)=0$. The condition $u(1)=0$ would
be satisfied automatically due to the symmetry of $u(y)$ \cite{Semenko2023}.
Here and further on the derivatives by $y$ are marked as strokes,
\[
u_{0}'=\frac{\mathrm{d}u_{0}}{\mathrm{d}y}.
\]

It has to be noted that (\ref{steady_sol}) in general could not be
found analytically, unlike the case of simpler models like Oldroyd-B
\cite{Sanchez,PortDenn,Chaudhary}. That means the value of $u_{H}$
is not known a-priori for a given set of parameters and has to be
calculated numerically. Technically that means we have to determine
the value of pressure gradient $A$ such that $u_{0}(0)=1$.

Having the steady solution, now we can formulate the spectral problem
for linear stability. Let 
\[
U_{0}=\left(\begin{array}{c}
u_{0}(y)\\
v_{0}(y)\\
p_{0}(x,y)\\
a_{110}(y)\\
a_{120}(y)\\
a_{220}(y)
\end{array}\right).
\]
We are looking for solution of (\ref{continuity_eq})-(\ref{noslip_cond})
of the form
\begin{equation}
\hat{U}(t,x,y)=\left(\begin{array}{c}
\hat{u}\\
\hat{v}\\
\hat{p}\\
\hat{a}_{11}\\
\hat{a}_{12}\\
\hat{a}_{22}
\end{array}\right)=U_{0}+U(y)e^{\lambda t+i\omega x}=U_{0}+\left(\begin{array}{c}
u(y)\\
v(y)\\
p(y)\\
a_{11}(y)\\
a_{12}(y)\\
a_{22}(y)
\end{array}\right)e^{\lambda t+i\omega x}.\label{solution_harm}
\end{equation}
The presence of the solutions with $\mathrm{real}(\lambda)>0$ would
indicate the temporal instability of the flow.

Abandoning nonlinear terms, we can obtain the following system for
$v$ and $a_{jk},\;k=1,2$.
\begin{multline}
\lambda(v''-\omega^{2}v)=-i\omega u_{0}v''+(i\omega^{3}u_{0}+i\omega u_{0}'')v+\\
+\omega^{2}a_{11}'-i\omega a_{12}''-i\omega^{3}a_{12}-\omega^{2}a_{22}',\label{v_linear_eq}
\end{multline}
\begin{multline}
\lambda\omega\mathrm{W}a_{11}+2\omega\left(\frac{1}{\mathrm{Re}}+\mathrm{W}a_{110}\right)v'-2i\mathrm{W}a_{120}v''+\mathrm{W}\omega a_{110}'v+\\
+\omega\left(2\mathrm{ReW}\frac{k-\beta}{3}a_{110}+\mathrm{ReW}\frac{k-\beta}{3}a_{220}+1+2\beta\mathrm{ReW}a_{110}+iu_{0}\mathrm{W}\omega\right)a_{11}+\\
+\left(2\beta\mathrm{ReW}\omega a_{120}-2u_{0}'\mathrm{W}\omega\right)a_{12}+\mathrm{ReW}\omega\frac{k-\beta}{3}a_{110}a_{22}=0,\label{a11_linear_eq}
\end{multline}
\begin{multline}
\lambda\omega\mathrm{W}a_{12}-i\left(\frac{1}{\mathrm{Re}}+\mathrm{W}a_{220}\right)v''+\left(\mathrm{W}\omega a_{120}'-\left(\frac{1}{\mathrm{Re}}+\mathrm{W}a_{110}\right)i\omega^{2}\right)v+\\
+\mathrm{ReW}\omega\frac{k+2\beta}{3}a_{120}a_{11}+\omega\left(1+\mathrm{ReW}\frac{k+2\beta}{3}\left(a_{110}+a_{220}\right)+i\omega\mathrm{W}u_{0}\right)a_{12}+\\
+\omega\mathrm{W}\left(\mathrm{Re}\frac{k+2\beta}{3}a_{120}-u_{0}'\right)a_{22}=0.\label{a12_linear_eq}
\end{multline}
\begin{multline}
\lambda\mathrm{W}a_{22}+\mathrm{W}\left(a_{220}'-2i\omega a_{120}\right)v-2\left(\frac{1}{\mathrm{Re}}+\mathrm{W}a_{220}\right)v'+\\
\mathrm{ReW}\frac{k-\beta}{3}a_{220}a_{11}+2\beta\mathrm{ReW}a_{120}a_{12}+\\
+\left(i\omega\mathrm{W}u_{0}+\mathrm{ReW}\frac{k-\beta}{3}\left(a_{110}+2a_{220}\right)+1+2\beta\mathrm{ReW}a_{220}\right)a_{22}=0.\label{a22_linear_eq}
\end{multline}
Since the continuity equation should be satisfied at the boundaries
of the channel, the system (\ref{v_linear_eq})-(\ref{a22_linear_eq})
is supplemented with boundary conditions
\begin{equation}
v(\pm1)=0,\;v'(\pm1)=0.\label{v_boundary_cond}
\end{equation}

\textbf{Remark 1.} We should point out that the validity of no-slip
conditions (\ref{v_boundary_cond}) for polymer solutions could be disputed since the slipping phenomenon for them is observed
under the certain conditions \cite{AltPyshPysh,Sochi}. But we believe the
conditions (\ref{v_boundary_cond}) are acceptable for the analysis
of the stability properties of the model, even more so that the same
conditions are usually assumed for similar research of other rheological
models, such as UCM \cite{PortDenn,Gorodtsov}, Oldroyd-B \cite{Chaudhary,Hansen},
FENE-P \cite{Buza}, Giesecus \cite{Azaiez} and others, and identical
conditions here are essential if the one is going to compare the results
of these studies with the present one.

\textbf{Remark 2.} Note that we never specified the pressure gradient
$A$ while presented the steady solution $U_{0}$. Since we chose
the velocity scale $u_{H}$ to be a maximum value of the velocity,
i.e. $u_{0}(0)$, the constant $A$ should be taken such as $u_{0}(0)=1$.
For Navier--Stokes equations, and even for some viscoelastic models,
this condition simply lead to $A=2/\mathrm{Re}$. But here the
solution $U_{0}$ is not known precisely, which does not allow us to
get the exact formula for $A$. Thus we determine $A$ by simple bisection
numerical method while numerically integrating the equation (\ref{u0_eq})
from $y=-1$ to $y=0$. 

The problem (\ref{v_linear_eq})-(\ref{v_boundary_cond}) has three
major parameters: $\mathrm{Re}$, $\mathrm{W}$ and $\beta$ (remember,
that the ratio $k=1.2$$\beta$ is postulated here according to \cite{Golovicheva}).
It is important to notice that the number $\mathrm{W}$ represent
the ratio of elastic and viscous forces. If the relaxation time of
the media $\tau_{0}=0$ then $\mathrm{W}=0$ also and the elastic
forces are absent. The problem thus would be reduced to the well-known
stability problem for the Navier--Stokes equations. Indeed, if $\mathrm{W}=0$
then from (\ref{a11_linear_eq})-(\ref{a22_linear_eq}) we have
\[
\begin{aligned} & a_{11}=-\frac{2}{\mathrm{Re}}v',\\
 & a_{12}=\frac{i}{\omega\mathrm{Re}}v''+\frac{i\omega}{\mathrm{Re}}v,\\
 & a_{22}=\frac{2}{\mathrm{Re}}v'.
\end{aligned}
\]
From (\ref{v_linear_eq}) we then have 

\begin{equation}
\begin{aligned} & \lambda(v''-\omega^{2}v)+i\omega u_{0}v''-i\omega\left(u_{0}''+\omega^{2}u_{0}\right)=\frac{1}{\mathrm{Re}}\left(v^{(4)}-2\omega v''+\omega^{4}v\right),\\
 & v(\pm1)=0,\;v'(\pm1)=0.
\end{aligned}
\label{Orr-Sommer}
\end{equation}
which is the known boundary-value problem for Orr-Sommerfeld equation \cite{Drazin}.
The steady solution $U_{0}$ for the case $\mathrm{W}=0$ would also
be classical Poiseuille parabolic flow. That fact would allow us to
compare the calculations for mVP model for small $\mathrm{W}$ with
the known results for Navier-Stokes equations. 

\section{Numerical results}

\subsection{Calculation of the spectrum}

The eigenvalue problem (\ref{v_linear_eq})-(\ref{v_boundary_cond})
is solved with the help of spectral method by decomposition of the
solution into polynomials on Chebyshev grid. The Chebyshev points
on the interval $-1\le y\le1$ are
\[
y_{j}=\cos(\pi j/N),\;j=0,\dots,N.
\]
The unknown functions $v(y)$ and $a_{jk},(y)\;k=1,2$ are discretized
on that grid, thus transforming into column vectors $\vec{v}$ and
$\vec{a}_{jk}$ respectively, and their derivatives are approximated
with Chebyshev differential matrix $D$ \cite{Boyd,Thefeten}:
\[
v^{(n)}\approx D^{n}\vec{v}.
\]
Let then $\vec{U}$ be a discretized vector of the length $4(N+1)$,
\[
\vec{U}=\left(\begin{array}{c}
\vec{v}\\
\vec{a}_{11}\\
\vec{a}_{12}\\
\vec{a}_{22}
\end{array}\right).
\]
The system of differential equations (\ref{v_linear_eq})-(\ref{a22_linear_eq})
then would be transformed into the generalized algebraic eigenvalue
problem 
\begin{equation}
M_{1}\vec{U}=\lambda M_{2}\vec{U},\label{general_eig}
\end{equation}
where $M_{1}$ and $M_{2}$ are some matrices which depend on $\omega$
and $U_0$. The eigenvalues and eigenvectors of (\ref{general_eig})
are our approximation of the solution of (\ref{v_linear_eq})-(\ref{v_boundary_cond}).
The problem (\ref{general_eig}) is solved by the procedure ``eig''
of MATLAB.

\textbf{Remark 3.} Here we have omitted the details of the way the
boundary conditions are approximated. Please refer to the Appendix
for that, along with the discussion of the adequacy of the method
and possible ways of it's verification. 

\subsection{Analysis of the unstable modes}

We will start the presentation of the results with  $\mathrm{Re}=5000$, which is slightly below the critical value for viscous newtonean fluid. 
It would allow us to demonstrate both stabilizing and destabilizing effects of elastic forces on the plane flow. We would also assume $\omega=1$ unless stated otherwise.

First let us test the numerical algorithm for $\mathrm{W}=0$ to compare
the results with the one for Navier-Stokes equations. Even if mathematically
problems (\ref{v_linear_eq})-(\ref{v_boundary_cond}) and (\ref{Orr-Sommer})
are identical for $\mathrm{W}=0$, numerically they are not, since the
former is the system of second-order equations, while the latter is one
equation of forth order. The Fig. \ref{fig:Spectra-of-Orr-Sommerfeld}
shows the calculated spectrum of these two problems, and it can be
seen that the spectrum is identical, which is the evidence for the
algorithm to perform correctly.
\begin{figure}
\noindent \centering{}\includegraphics{Figures/spectreW0_wide}\caption{Spectra of Orr-Sommerfeld equation (red circles) and mVP model at
$\mathrm{W}=0$ (blue dots)\label{fig:Spectra-of-Orr-Sommerfeld}}
\end{figure}

The Fig. \ref{fig:Spectra-of-mVP} shows the spectrum evolution with
gradual increase of the Weissenberg number. It can be seen that the
flow could become unstable for sufficiently large $\mathrm{W}$ even
for the values of $\mathrm{Re}$ which are below the critical value
for Navier--Stokes equations. Two eigenvalues there are marked as $A$ and $B$ and it can be seen that if $\mathrm W$ is sufficiently large, 
these two eigenvalues obtain positive real parts, thus indicating the instability of the flow. Note that  the spectrum for viscous
fluids contain only one unstable mode \cite{Drazin}, while here we have at least two. The Fig. \ref{fig:Unstable-modes}
shows two unstable modes for these eigenvalues, marked as mode A and mode B respectively. One of that
modes (mode A) is symmetrical, and other (mode B) is asymmetrical.
\begin{figure}
\noindent \begin{centering}
\includegraphics{Figures/spectreW0d3_marked}\includegraphics{Figures/spectreW0d6_marked}\\
\includegraphics{Figures/spectreW1d3_marked}\includegraphics{Figures/spectreW3_marked}
\par\end{centering}
\caption{Spectra of mVP with various $\mathrm{W}$. Red circles are the position
of eigenvalues at $\mathrm{W}=0$\label{fig:Spectra-of-mVP}}
\end{figure}
\begin{figure}
\noindent \centering{}\includegraphics{Figures/eig1}\includegraphics{Figures/eig2}\caption{Unstable modes $v(y)$\label{fig:Unstable-modes}}
\end{figure}

It appears that mode A is the continuation of Tollmien--Schlichting
mode \cite{Drazin,Schmidt}. The Fig. \ref{fig:Trajectories-of-eigenvalues}
shown the trajectory of that mode with change of $\mathrm{W}$. The
mode B on the other hand is the continuation of the mode which never
become unstable for Navier--Stokes equations, but could become unstable
in viscoelastic flow. The same behavior of viscoelastic flow was reported
for other rheological models \cite{PortDenn,Chaudhary}, and the fact
that this mode become unstable only for that type of fluids allows
the one to connect these modes with the instability caused by elastic
forces. 
\begin{figure}
\noindent \centering{}\includegraphics{Figures/eig_trajectories}\caption{Trajectories of eigenvalues $A$ and $B$ from $\mathrm{W}=0$ to
$\mathrm{W}=3$\label{fig:Trajectories-of-eigenvalues}}
\end{figure}

Let us now calculate the critical values of parameters for the stability
of plane steady flow. Due to the practical purposes we would study
the stability of the flow not in the $\mathrm{Re}-\mathrm{W}$ plane
but rather in the $\mathrm{Re}-E$ plane, where $E=\mathrm{W}/\mathrm{Re}$
is the elastic number. Unlike $\mathrm{W}$, the value of $E$ is determined
only by rheological properties of the fluid and the geometry of the
channel and does not depend on the velocity of the flow, which make
it fixed for experimental setup in the slit channel. We will limit
the numerical experiments with the maximum value of $\mathrm{Re}=2\cdot10^{4}$.
The higher number of $\mathrm{Re}$ decrease the accuracy of the calculations,
and in our opinion the chosen maximum value is enough to describe
the general nature of the problem in sufficient details.

The stability
curve (Fig. \ref{fig:Zones-of-stability}) shows that the critical Reynolds number for $E\approx 0$ is indeed close to the one for the newtonean fluid (approximately 5772) but the increase of $E$ causes the rapid grow of critical $\mathrm{Re}$, while even larger $E$ lead to its eventual decrease. Our calculations did not reveal a minimum value
of $\mathrm{Re}$ which guarantees the stable flow regardless of the
$E$. Apparently, the sufficiently large $E$ can destabilize the
flow of any predetermined velocity, which matches with the findings
for number of other rheological models \cite{Sanchez,PortDenn,Chaudhary,Buza}.
\begin{figure}
\noindent \begin{centering}
\includegraphics{Figures/Re_E}
\par\end{centering}
\caption{Zones of stability and instability, $\mathrm{Re}-E$ plane\label{fig:Zones-of-stability}}
\end{figure}

The sharp raise of $\mathrm{Re}-E$ curve at low $E$ happens due
to the disappearance of the fragment of neutral stability curve of
the eigenmode A (Fig. \ref{fig:Neutral-stability-curves-shrink}).
On this figure it is evident that the neutral stability curve closes
and shrinks to nothing with increasing $E$. The critical value of
$E\approx 26\cdot 10^{-6}$ is when the curve disappears which causes the critical value
of $\mathrm{Re}$ to jump instantaneously to the next fragment of
the stability curve situated somewhere higher than the upper limit
of $\mathrm{Re}$ at which the calculations were performed. Further
increase of $E$ lead to downward movement of neutral stability curves
for modes $A$ and $B$, Fig. \ref{fig:Neutral-stability-curves}.
\begin{figure}
\noindent \centering{}\includegraphics{Figures/Re_om_shrink}\caption{Neutral stability curves of mode $A$\label{fig:Neutral-stability-curves-shrink}}
\end{figure}

\begin{figure}
\noindent \centering{}\includegraphics{Figures/m2_Re_om_E0d0001}\includegraphics{Figures/m2_Re_om_E0d0002}\\
\includegraphics{Figures/m2_Re_om_E0d001}\includegraphics{Figures/m2_Re_om_E0d003}\caption{Neutral stability curves of modes $A$ (blue) and $B$ (orange)\label{fig:Neutral-stability-curves}}
\end{figure}


\subsection{Continuous spectrum}

The area of the complex plane shown on the Fig. \ref{fig:Cont-spec-num}
contains the unstructured cloud of eigenvalues. Numerical experiments
with various grid size $N$ reveal that the position of this eigenvalues
depend on $N$ and thus the accuracy of calculations for that part
of the spectrum is questionable. For better understanding of that
part of the spectrum let us rewrite the equations (\ref{a11_linear_eq})-(\ref{a22_linear_eq})
as
\begin{equation}
\begin{aligned} & (\omega\mathrm{W}\lambda+A_{1})a_{11}+B_{1}a_{12}+C_{1}a_{22}=V_{10}v+V_{11}v'+V_{12}v'',\\
 & A_{2}a_{11}+(\omega\mathrm{W}\lambda+B_{2})a_{12}+C_{2}a_{22}=V_{20}v+V_{22}v'',\\
 & A_{3}a_{11}+B_{3}a_{12}+(\mathrm{W}\lambda+C_{3})a_{22}=V_{30}v+V_{31}v'.
\end{aligned}
\label{aij_system}
\end{equation}
\begin{figure}
\noindent \centering{}\includegraphics{Figures/cont_spec_num}\caption{Fragment of the spectrum at $\mathrm{Re}=5000,\mathrm{W}=0.2$\label{fig:Cont-spec-num}}
\end{figure}
Here, as it follows from (\ref{a11_linear_eq})-(\ref{a22_linear_eq}),
the coefficients $A_{j},\,B_{j},\,C_{j}$, $V_{j0},\,V_{j1},\,V_{j2}$,
$j=1,2,3$ are functions of $y$ according to their dependence on the steady
solution $U_{0}$. Let 
\[
M=\left(\begin{array}{ccc}
\omega\mathrm{W}\lambda+A_{1} & B_{1} & C_{1}\\
A_{2} & \omega\mathrm{W}\lambda+B_{2} & C_{2}\\
A_{3} & B_{3} & \mathrm{W}\lambda+C_{3}
\end{array}\right).
\]
By Cramer's rule, 
\begin{equation}
a_{jm}=\frac{\det(M_{jm})}{\det(M)},\;j,m=1,2,\label{aij_solution}
\end{equation}
where $M_{jm}$ is the matrix formed by replacing corresponding column
of the matrix $M$ with the vector of right-hand sides of (\ref{aij_system}).
It is clear that $\det M$ is the polynomial of power three of $\lambda$
with functions of $y$ as coefficients. By substituting (\ref{aij_solution})
to (\ref{v_linear_eq}), we can see that the resulting equation for
$v$ will be linear of the fourth-order since (\ref{v_linear_eq})
contains $a_{12}''$ and from (\ref{aij_solution}) it follows that
$a_{12}$ in turn linearly depends on $v''$. The coefficient at $v^{(4)}$
of the resulting equation would be proportional to
\begin{multline}
d_{4}=\frac{(\omega\mathrm{W}\lambda+A_{1})(\mathrm{W}\lambda+C_{3})V_{22}}{\det(M)}- \\ -\frac{A_{2}(\mathrm{W}\lambda+C_{3})A_{2}V_{12}-V_{12}C_{2}A_{3}+V_{22}A_{3}C_{1}}{\det M}.\label{d4_coeff}
\end{multline}
Note that the coefficient at $v^{(3)}$ of that equation would contain
the term proportional to $d_{4}'$. By multiplying the equation by
$\left(\det M\right)^{2}$, we get that the coefficient at $v^{(4)}$
would be zero if either numerator of $d_{4}$ is zero or $\det M=0$,
while the coefficient of $v^{(3)}$ would be non-zero in both cases.
That makes the roots of 5th-order equation
\begin{multline}
((\omega\mathrm{W}\lambda+A_{1})(\mathrm{W}\lambda+C_{3})V_{22}-A_{2}(\mathrm{W}\lambda+C_{3})A_{2}V_{12}+\\+V_{12}C_{2}A_{3}-V_{22}A_{3}C_{1})\det M=0\label{cont_spec_eq}
\end{multline}
the eigenvalues of (\ref{v_linear_eq})-(\ref{v_boundary_cond}).
Since the coefficients of (\ref{cont_spec_eq}) are functions of $y$,
we have the continuous spectrum of (\ref{v_linear_eq})-(\ref{v_boundary_cond})
with five eigenvalues for any $y$, $-1\leq y\leq1$. By mapping these
roots on the Fig. \ref{fig:Cont-spec-full} and comparing them with
the spectrum calculated in previous subsection, we can see that the
correlation between the two is present. The ends of branches of continuous
spectrum are captured with numerical algorithm, while the rest of
numerical eigenvalues create the unstable numerical cloud around continuous
spectrum. This behavior of spectral methods used on the problems with
continuous spectra is known and appear even in the trivial problems.
The important thing here is that even numerical algorithm does not
capture the continuous spectrum with good accuracy, this spectrum
is located in the left complex half-plane and does not provide instability
of the steady state.

\begin{figure}
\noindent \begin{centering}
\includegraphics{Figures/cont_spec}
\par\end{centering}
\caption{Numerical spectrum (blue) and solution of \eqref{cont_spec_eq} (red)\label{fig:Cont-spec-full}}
\end{figure}

The continuous spectrum of rheological models has been observed before,
see \cite{Renardy,Wilson,Graham} for example. The exact shape of
the spectrum strongly depend on the nonlinearity of constitutive equation
and thus is unique for the particular model.

\section*{Appendix}

\subsection*{The numerical scheme and the approximation of the boundary conditions}

The idea of spectral method is to approximate the function on the
grid of $N+1$ nodes by interpolating polynomial of power $N$. The
derivatives of the function thus are calculated as derivatives of
said polynomial. That is, if the values of function $v(y)$ are known
in the nodes of the grid $y_{j}$, $j=0,\dots,N$ and$\vec{v}$ is
the vector of $v(y_{j})$, then the vector
\[
\vec{v}^{(n)}=D^{n}\vec{v}
\]
is the approximation of derivative $v^{(n)}$ in the same nodes, where
$D$ is the constant differentiation matrix, independent on $v$.
Naturally, the matrix $D$ is singular, since polynomial does not
uniquely defined by its derivative. So the uniqueness of the solution
of boundary-value problem is established by implementing boundary
conditions. In the present research it was done by decomposing the
unknown functions into the polynomials which satisfies the boundary
conditions of the problem. Namely, by taking conditions (\ref{v_boundary_cond})
into account, we approximate $v(y)$ by polynomial $(1-y^{2})q(y)$,
where $q(y)$ is the polynomial of power $N$ and
\[
q(y_{j})=v(y_{j})/(1-y_{j})^{2},\;j=1,\dots,N-1,
\]
\[
q(\pm1)=0.
\]
It is clear that the (\ref{v_boundary_cond}) would be satisfied with
this approximation \cite{Boyd} The derivatives of $v(y)$ in the
grid nodes $y_{j}$, $j=1,\dots,N$ are approximates as
\[
\begin{aligned} & v^{'}(y_{j})\approx\frac{\mathrm{d}}{\mathrm{d}y}\left((1-y^{2})q(y)\right)|_{y=y_{j}}=-2y_{j}q(y_{j})+(1-y_{j}^{2})q'(y_{j}),\\
 & v^{''}(y_{j})\approx\frac{\partial^{2}}{\partial y^{2}}\left((1-y^{2})q(y)\right)|_{y=y_{j}}=-2q(y_{j})-4y_{j}q'(y_{j})+(1-y_{j}^{2})q''(y_{j}),\\
 & v^{'''}(y_{j})\approx\frac{\mathrm{d}^{3}}{\mathrm{d}y^{3}}\left((1-y^{2})q(y)\right)|_{y=y_{j}}=-6q'(y_{j})-6y_{j}q''(y_{j})+(1-y_{j}^{2})q'''(y_{j}),\\
 & v^{^{(4)}}(y_{j})\approx\frac{\mathrm{d}^{4}}{\mathrm{d}y^{4}}\left((1-y^{2})q(y)\right)|_{y=y_{j}}=-12q''(y_{j})-8y_{j}q'''(y_{j})+(1-y_{j}^{2})q^{(4)}(y_{j}).
\end{aligned}
\]
The formulas for approximation of the derivatives of fourth order
are needed to solve the equations of type (\ref{Orr-Sommer}). Otherwise,
the derivatives of the order up to second are sufficient since the
problem written in the form (\ref{v_linear_eq})-(\ref{v_boundary_cond})
contains only those.

The boundary conditions $q(\pm1)=0$ are implemented by removing first
and last row and column from differentiation matrices $D^{n}$, $n\le4$,
thus obtaining the regular matrices $\tilde{D^{n}}$ of the size $N-1\times N-1$
\cite{Boyd,Thefeten}. As the result, the vectors of derivatives for
$v$ in the inner nodes $y_{j},\;j=1,\dots,N-1$ are
\[
\begin{aligned} & \vec{v}'=\left(-2\mathrm{diag}(\vec{y}\vec{v})+\mathrm{diag}(1-(\vec{y})^{2})\tilde{D}\vec{v}\right)\mathrm{diag}\left(\frac{1}{1-(\vec{y})^{2}}\right),\\
 & \vec{v}''=\left(-2\mathrm{diag}(\vec{v})-4\mathrm{diag}(\vec{y})\tilde{D}\vec{v}+\mathrm{diag}(1-(\vec{y})^{2})\tilde{D^{2}}\vec{v}\right)\mathrm{diag}\left(\frac{1}{1-(\vec{y})^{2}}\right),\\
 & \mbox{etc.}
\end{aligned}
\]
Using that approximation, the problem (\ref{v_linear_eq})-(\ref{v_boundary_cond})
would be transformed into (\ref{general_eig}), where 
\[
\vec{U}=(v(y_{1}),\dots,v(y_{N-1}),a_{11}(y_{0}),\dots,a_{11}(y_{N}),a_{12}(y_{0}),\dots a_{22}(y_{N}))^{\top},
\]
and the length of $\vec{U}$ is $4N+2$. Note that the components
$v(y_{0})=v(y_{N})=0$ are not included in $\vec{U}$ since they are
already taken cared of by its interpolation. The matrix $M_{1}$ of
the size $4N+2\times4N+2$ consist of the blocks $M_{1}=\{m_{ln}\}$,
$l,n=1,\dots,4$,
\[
\begin{aligned} & m_{11}=-i\omega\mathrm{diag}(\vec{u}_{0})\tilde{D^{2}}+i\omega^{3}\mathrm{diag}(\vec{u}_{0}+i\omega\mathrm{diag}\left(\tilde{D^{2}}\vec{u}_{0}\right),\\
 & m_{12}=\left(\omega^{2}D\right)_{v},\\
 & m_{13}=\left(-i\omega D^{2}-i\omega^{3}I\right)_{v},\\
 & m_{14}=\left(\omega^{2}I\right)_{v},\\
 & m_{21}=\left(2i\mathrm{W}\mathrm{diag}(\vec{a}_{120})D^{2}-2\omega\left(I/\mathrm{Re}+\mathrm{W}\mathrm{diag}\left(\vec{a}_{110}\right)\right)D-\mathrm{W}\omega\mathrm{diag}\left(D\vec{a}_{110}\right)\right)_{h},\\
 & m_{22}=\omega\left(\mathrm{ReW}\frac{k-\beta}{3}\left(2\mathrm{diag}\left(\vec{a}_{110}\right)+\mathrm{diag}\left(\vec{a}_{220}\right)\right)+I\right.+ \\
 & +\left.2\beta\mathrm{ReW}\mathrm{diag}\left(\vec{a}_{110}\right)+i\mathrm{W}\omega\mathrm{diag}\left(\vec{u}_{0}\right)\right)_{h},\\
 & m_{23}=\left(2\beta\mathrm{ReW}\omega\mathrm{diag}\left(\vec{a}_{120}\right)-2\mathrm{W}\omega\mathrm{diag}\left(D\vec{u}_{0}\right)\right)_{h},\\
 & m_{24}=\mathrm{ReW}\omega\frac{k-\beta}{3}\left(\mathrm{diag}\left(\vec{a}_{110}\right)\right)_{h},\\
 & \mbox{etc.}
\end{aligned}
\]
Here $I$ is the unity matrix of $N+1\times N+1$, subscripts $"h"$
and $"v"$ denotes the removal of first and last rows and columns
of matrix respectively. Similarly, the $M_{2}$ is
\[
M_{2}=\left(\begin{array}{cccc}
\tilde{D^{2}}-\omega^{2}\tilde{I} & 0 & 0 & 0\\
0 & \omega\mathrm{W}I_{h} & 0 & 0\\
0 & 0 & \omega\mathrm{W}I_{h} & 0\\
0 & 0 & 0 & \mathrm{W}I_{h}
\end{array}\right),
\]
where $\tilde{I}$ is the unity matrix of $N-1\times N-1$. The eigenvalues
of (\ref{general_eig}) are found with ``eig'' procedure of MATLAB,
which utilizes the QZ-method. The method does not rely on inverse
matrices and work reliably even if $M_{1}$ or $M_{2}$ are singular
\cite{Kressner}.

\subsection*{Validation of the eigenvalues by shooting method}

The spectral methods are known for the high accuracy since they use
the information in all nodes of the grid to calculate the derivatives
of the function. Still, the numerical spectrum could be significantly
different from the real one. The method is prone to create phantom
eigenvalues which does not present in the actual problem, and the
position of the rest of the eigenvalues could also be off. In certain
situations up to the half of found spectrum could be discarded as
imprecise \cite{Schmidt}. One way to filter the ill-calculated eigenvalues
is to check the results on the grids of different sizes. In this research
we used the grids of the size $N=200$, $N=400$ and $N=800$. Only
the eigenvalues independent of the grid size were used for stability
analysis. But even the independence of numerical eigenvalues on the
$N$ does not guaranteed that these eigenvalues are not numerical
artifacts of the algorithm. Another way to check their validity is
to calculate them by alternative numerical method of different nature.
The shooting method could be used since it does not rely on the calculation
of the eigenvalues of large matrices. Let us briefly describe the
idea and the methods used for calculations in the present research.

As it was mentioned, the system (\ref{v_linear_eq})-(\ref{a22_linear_eq})
could be reduced to one linear fourth-order equation for $v$ with
coefficients as functions of $y$ and $\lambda$. By denoting the
linear operator as $\mathcal{L}_{4}(y,\lambda)$, we got the eigenvalue
problem
\begin{equation}
\begin{aligned} & \mathcal{L}_{4}(y,\lambda)v=0,\\
 & v(\pm1)=v'(\pm1)=0.
\end{aligned}
\label{v_bvp}
\end{equation}
The idea of shooting method is straightforward: we are trying to find
the non-zero solution of (\ref{v_bvp}) as linear combination of two
linearly independent solutions of the equation $\mathcal{L}_{4}v=0$
with two left boundary conditions. That is, we numerically solve two
initial-value problems
\[
\begin{aligned} & \mathcal{L}_{4}(y,\lambda)v_{1}=0,\\
 & v_{1}(-1)=v_{1}'(-1)=0,\\
 & v_{1}''(-1)=1,\;v_{1}'''(-1)=0
\end{aligned}
\;\mbox{and}\;\mbox{ \ensuremath{\begin{aligned} &  \mathcal{L}_{4}(y,\lambda)v_{2}=0,\\
  &  v_{2}(-1)=v_{2}'(-1)=0,\\
  &  v_{2}''(-1)=0,\;v_{2}'''(-1)=1. 
\end{aligned}
} }
\]
Any solution of (\ref{v_bvp}) could be expressed as $v(y)=c_{1}v_{1}(y)+c_{2}v_{2}(y)$.
The right boundary conditions of (\ref{v_bvp}) means
\[
\begin{aligned} & c_{1}v_{1}(1)+c_{2}v_{2}(1)=0,\\
 & c_{1}v_{1}'(1)+c_{2}v_{2}'(1)=0.
\end{aligned}
\]
Solution $v(y)$ is not trivial only if 
\begin{equation}
g(\lambda)=\left|\begin{array}{cc}
v_{1}(1) & v_{2}(1)\\
v_{1}'(1) & v_{2}'(1)
\end{array}\right|=0.\label{shooting_equation}
\end{equation}
The solution of (\ref{v_bvp}) is thus reduced to numerical solve
of (\ref{shooting_equation}) where the evaluation of the function
$g(\lambda)$ for any $\lambda$ requires the calculation of two linear
initial-value problems. Since the solutions of the equation $\mathcal{L}_{4}(y,\lambda)v=0$
could be stiff, we used the procedure ``ode15s'' of MATLAB (variable-step,
variable-order solver, see \cite{Shampine}) for that purpose. Still,
even this procedure returns high numerical errors of $g(\lambda)$
due to large magnitudes of solutions $v_{1,2}$ (up to $10^{20}$
at $y=1$) and their non-orthogonality. To combat this problem we
used the known method of intermediate orthogonalization of solutions
introduced by Godunov \cite{Godunov}. At each step of the numerical
integration of $v_{1,2}$ we checked the scalar product $(v_{1},v_{1}')^{\top}\cdot(v_{2},v_{2}')^{\top}$
and performed the Gram-Schmidt orthogonalization if this product exceeded
$10^{2}$. 

We used the Newton--Raphson method to calculate the solutions of
(\ref{shooting_equation}). This method requires good initial guess
of the solution for reliable convergence, so we have used the eigenvalues
obtained by spectral method as such. We assumed the eigenvalues are
genuine if the shooting method stabilizes in the short vicinity of
the initial position. This validation was used for all eigenvalues
of unstable modes and the results showed that eigenvalues calculated
by shooting method are not further away than $10^{-6}$ from the ones
obtained by the spectral method. We consider this match as sufficient
to believe that our calculations are accurate.

\section*{Conclusions}

The main result of present research we see in calculating the critical
values of Reynolds number for Poiseuille-type plane flow modeled by
mVP equations. It was discovered that the critical value is finite
and greater than zero for all values of Weissenberg number we have
checked. The latter indicates that plane Poiseuille-like flow for
mVP is not unconditionally unstable. Even more so, the results of
our calculations for limit $E \to0$ provide us with critical
value of $\mathrm{Re}$ which is identical to the one for Navier--Stokes
equations, which both support the accuracy of performed calculations
and adequacy of the rheological model itself, since it is expected
that the stability of polymer solution with almost absent elastic
forces should closely resemble the one for non-elastic viscous fluid.
The two unstable modes (symmetric and asymmetric) were discovered,
which correlates with the results for UCM and Oldroyd-B. But, unlike
the latter, we showed that for mVP the slight increase of Weissenberg
number from zero provides stabilizing effect on the plane flow, even
if the further increase of $\mathrm{W}$ (or $E$) leads to eventual drop of
critical value of $\mathrm{Re}$, apparently to arbitrarily low values.
This effect could reflect the known experimental fact that the addition of small amount of polymer into the viscous solvent delays the appear of turbulence \cite{Samanta}, but we assume that further research is required to verify
the correlation between these phenomena.



\bigskip
\begin{thebibliography}{1} % Ссылки на литературу приводятся только на английском языке (даже если у источника нет англоязычной версии).


\bibitem{DoiEdw}M. Doi, S.F. Edwards, The Theory of Polymer Dynamics,
Clarendon: Oxford, 1986.

\bibitem{Bird}R.B. Bird, C.F. Curtiss, R.C. Armstrong, O. Hassager,
\href{https://doi.org/10.1002/pol.1987.140251210}{\it Dynamics of Polymeric Liquids, Vol. 2},
Wiley, New York, 1987.

\bibitem{Samanta}D. Samanta, Y. Dubief, M. Holzner, C. Schдfer, A.N.
Morozov, C.Wagner, B. Hof, \href{https://doi.org/10.1073/pnas.1219666110}{\it Elasto-inertial turbulence},
Proceedings of the National Academy of Sciences, \textbf{110}:26 (2013),
10557--10562.

\bibitem{Srinivas}S.S. Srinivas, V. Kumaran, 2017 \href{https://doi.org/10.1017/jfm.2016.839}{\it Effect of viscoelasticity on the soft-wall transition and turbulence in a microchannel}.
J. Fluid Mech., \textbf{812} (2017), 1076--1118.

\bibitem{Oldroyd} J.G. Oldroyd, \href{https://www.jstor.org/stable/98422}{\it On the formulation of rheological equations of state},
Proc. R. Soc., \textbf{200}:1063, 523--541.

\bibitem{Shafeq}E.S.G. Shaqfeh, \href{https://doi.org/10.1146/annurev.fl.28.010196.001021}{\it Purely elastic instabilities in viscometric flows},
Annual Review of Fluid Mechanics, \textbf{28}:1 (1996), 129--185.

\bibitem{Sanchez}H.A. Castillo Sanchez, M.R. Jovanovic, S. Kumar,
A. Morozov, V. Shankar, G. Subramanian, H.J. Wilson, \href{https://doi.org/10.1016/j.jnnfm.2022.104742}{\it Understanding viscoelastic flow instabilities: Oldroyd-B and beyond},
J. Non-Newtonian Fluid Mech., \textbf{302} (2022), 104742.

\bibitem{Sureshkumar}R. Sureshkumar, A.N. Beris, \href{https://doi.org/10.1016/0377-0257(94)01279-Q}{\it Linear stability analysis of viscoelastic Poiseuille flow using an Arnoldi-based orthogonalization algorithm},
J. Non-Newtonian Fluid Mech., \textbf{56} (1995), 151--182.

\bibitem{Larson}R.G. Larson, \href{https://doi.org/10.1007/BF00366504}{\it Instabilities in viscoelastic flows},
Rheola Acta, \textbf{31} (1992), 213--263.

\bibitem{Ramachandran}A. Ramachandran, D. T. Leighton, \href{https://doi.org/10.1017/S0022112008000980}{\it The influence of secondary flows induced by normal stress differences on the shear-induced migration of particles in concentrated suspensions},
J. Fluid Mech., \textbf{603} (2008), 207243.

\bibitem{Pokrovskii}V.N. Pokrovskii, \href{https://doi.org/10.1007/978-90-481-2231-8}{\it The mesoscopic theory of polymer dynamics},
Springer, London, 2010.

\bibitem{AltGusPysh}Yu.A. Altukhov, A.S. Gusev, G.V. Pyshnograi,
{\it Vvedeniye v mezoskopicheskuyu teoriyu tekuchikh polimernykh sistem},
Barnaul: AltGPA, 2012 (in Russian).

\bibitem{Makarova}M.A. Makarova, A.S. Malygina, G.V. Pyshnograi,
G.O. Rudakov, \href{https://doi.org/10.1134/S0021894421070142}{\it Simulation of rheological properties of polyethylene melts under uniaxial tension},
J. Appl. Mech. Tech. Phys., \textbf{62} (2021), 1063--1071.

\bibitem{Pysh2020}G.V. Pyshnograi, N. A. Cherpakova, H.N.A.A. Joda,
\href{https://doi.org/10.1007/s10891-020-02159-8}{\it Special features of nonlinear behavior of a polymer solution on large periodic deformations},
J. Eng. Phys. Thermophys., \textbf{93} (2020), 617--625.

\bibitem{Koshelev}K. Koshelev, A. Kuznetcov, D. Merzlikina, G. Pyshnograi,
I. Pyshnograi, M. Y. Tolstykh, \href{http://10.1088/1742-6596/790/1/012014}{\it Mesoscopic rhelogical model for polymeric media flows},
J. Phys.: Conf. Ser., \textbf{790} (2017), 012014.

\bibitem{Merzlikina}D. A. Merzlikina, G. V. Pyshnograi, R. Pivokonskii,
P. Filip, \href{https://doi.org/10.1007/s10891-016-1423-7}{\it Rheological model for describing viscometric flows of melts of branched polymers},
J. Eng. Phys. Thermophys., \textbf{89} (2016), 652--659.

\bibitem{BlokhYegTk}A.M. Blokhin, A.V. Yegitov, D.L. Tkachev, \href{https://doi.org/10.1134/S0965542515050073}{\it Linear instability of solutions in a mathematical model describing polymer flows in an infinite channel},
Comput. Math. and Math. Phys., \textbf{55}:5 (2015), 848--873.

\bibitem{TkachBib}D.L. Tkachev, E.A. Biberdorf, \href{https://doi.org/10.33048/semi.2023.20.076}{\it Spectrum of a problem about the flow of a polymeric viscoelastic fluid in a cylindrical channel (Vinogradov-Pokrovski model)},
Siberian Electronic Mathematical Reports, \textbf{20}:2 (2023), 1269--1289.

\bibitem{Tkachev2023}D.L. Tkachev, \href{https://doi.org/10.1016/j.jmaa.2022.126914}{\it Spectrum and linear Lyapunov instability of a resting state for flows of an incompressible polymeric fluid},
J. of Mathematical Analysis and Applications, \textbf{522}:1 (2023),
126914.

\bibitem{Semenko2023}R.E. Semenko, \href{https://doi.org/10.1063/5.0140477}{\it Stabilization of the at Poiseuille-type flow for viscoelastic polymeric liquid},
Phys. Fluids, \textbf{35} (2023), 033112.

\bibitem{Golovicheva}I.E. Golovicheva, S.A. Zinovich, G.V. Pyshnograi,
\href{https://doi.org/10.1007/BF02465279}{\it Effect of the molecular mass on the shear and longitudinal viscosity of linear polymers},
J. Appl. Mech. Tech. Phys., \textbf{41} (2000), 347--352.

\bibitem{BlokhBamb}A.M. Blokhin, N.V. Bambaeva, \href{https://doi.org/10.1134/S0965542514050054}{\it Stationary solutions equations of an incompressible viscoelastic polymer liquid},
Comput. Math. and Mach. Phys, \textbf{54}:5 (2014), 874--899.

\bibitem{PortDenn}K.C. Porteous, M.M. Denn, \href{https://doi.org/10.1122/1.549279}{\it Linear stability of plane Poiseuille flow of viscoelastic liquids},
Transactions of the Society of Rheology, \textbf{16}:2 (1972), 295--308.

\bibitem{Chaudhary}I. Chaudhary, P. Garg, V. Shankar, G. Subramanian,
\href{https://doi.org/10.1017/jfm.2019.759}{\it Elasto-inertial wall mode instabilities in viscoelastic plane Poiseuille flow},
J. Fluid Mech, \textbf{881} (2019), 119--163.

\bibitem{AltPyshPysh}Yu.A. Altukhov, G.V. Pyshnograi, I.G. Pyshnograi,
\href{http://dx.doi.org/10.4236/wjm.2011.16037}{\it Slipping phenomenon in polymeric fluids ow between parallel planes}
, World Journ. of Mech., \textbf{1} (2011), 294--298.

\bibitem{Sochi}T. Sochi, \href{https://doi.org/10.1080/15583724.2011.615961}{\it Slip at fluid-solid interface},
Polymer Rev., \textbf{51}:4 (2011), 309--340.

\bibitem{Gorodtsov}V. A. Gorodtsov, A. I. Leonov, \href{https://doi.org/10.1016/0021-8928(67)90156-6}{\it On a linear instability of a plane parallel Couette flow of viscoelastic fluid},
J. Applied Mathematics and Mechanics, \textbf{31}:2 (1967), 289--299.

\bibitem{Hansen}R.J. Hansen, R. Little, P.G. Forame, \href{https://doi.org/10.1252/jcej.6.310}{\it Experimental and theoretical studies of early turbulence},
Journal of Chemical Engineering of Japan, \textbf{6}:4 (1973), 310--314.

\bibitem{Buza}G. Buza, J. Page, R.R. Kerswell, \href{https://doi.org/10.1017/jfm.2022.222}{\it Weakly nonlinear analysis of the viscoelastic instability in channel flow for finite and vanishing Reynolds numbers},
J. Fluid Mech., \textbf{940} (2022), A11.

\bibitem{Azaiez}J. Azaiez, G. M. Homsy, \href{https://doi.org/10.1017/S0022112094001254}{\it Linear stability of free shear flow of viscoelastic liquids},
J. Fluid Mech., \textbf{268} (1994), 37--69.

\bibitem{Boyd}J.P. Boyd, \href{https://link.springer.com/book/9783540514879}{\it Chebyshev and Fourier Spectral Methods.},
Springer Berlin, Heidelberg, 1989.

\bibitem{Thefeten}L.N. Trefethen, \href{https://doi.org/10.1137/1.9780898719598}{\it Spectral methods in MATLAB},
SIAM, Philadelphia, 2000.

\bibitem{Drazin}P. G. Drazin, \href{https://doi.org/10.1017/CBO9780511809064}{\it Introduction to Hydrodynamic Stability},
Cambridge University Press, Cambridge, 2002.

\bibitem{Schmidt}P.J. Schmid, D.S. Henningson, \href{https://doi.org/10.1007/978-1-4613-0185-1}{\it Stability and Transition in Shear Flows},
Springer Science+Business Media, New York, 2001.


\bibitem{Renardy}M. Renardy, Y. Renardy, \href{https://doi.org/10.1016/0377-0257(86)80002-7}{\it Linear stability of plane Couette flow of an upper convected Maxwell fluid}.
J. Non-Newtonian Fluid Mech., \textbf{22} (1986), 23--33.

\bibitem{Wilson}H.J. Wilson, M. Renardy, Y. Renardy, \href{https://doi.org/10.1016/S0377-0257(98)00087-1}{\it Structure of the spectrum in zero Reynolds number shear flow of the UCM and Oldroyd-B liquids},
J. Non-Newtonian Fluid Mech., \textbf{80} (1999), 251--268.

\bibitem{Graham}M.D. Graham, \href{https://doi.org/10.1017/S0022112098008611}{\it Effect of axial flow on viscoelastic Taylor\textendash Couette instability},
J. Fluid Mech., \textbf{360} (1998), 341--374.

\bibitem{Kressner}D. Kressner, \textit{\href{https://doi.org/10.1007/3-540-28502-4}{\it Numerical Methods for General and Structured Eigenvalue Problems}},
Springer-Verlag Berlin, Heidelberg, 2005.

\bibitem{Shampine}L.F. Shampine, M.W. Reichelt, J.A. Kierzenka, \href{https://doi.org/10.1137/S003614459933425X}{\it Solving Index-1 DAEs in MATLAB and Simulink},
SIAM Review, \textbf{41} (1999), 538--552.

\bibitem{Godunov}S.K. Godunov, \href{https://www.mathnet.ru/php/archive.phtml?wshow=paper&jrnid=rm&paperid=6626&option_lang=rus}{\it O chislennom reshenii krayevykh zadach dlya sistem obyknovennykh differentsialnykh uravneniy},
Usp. Mat. Nauk., \textbf{16} (1961), 171--174 (in Russian).

\end{thebibliography}
\end{document}
