# Spatiotemporal linear stability of viscoelastic free shear flows: non-affine
response regime

D. Bansal$\dagger$, D. Ghosh$\dagger$ and S. Sircar$\dagger,\ddagger$
$\dagger$ Department of Mathematics, IIIT Delhi, India 110020
$\ddagger$ Corresponding Author (email<EMAIL_ADDRESS>
###### Abstract

We provide a detailed comparison of the two-dimensional, temporal and the
spatiotemporal linearized analyses of the viscoelastic free shear flows in the
limit of low to moderate Reynolds number and Elasticity number obeying four
different types of stress-strain constitutive equations: Oldroyd-B, Upper
Convected Maxwell, Johnson-Segalman (JS) and linear Phan-Thien Tanner (PTT).
The resulting fourth-order Orr-Sommerfeld Equation is transformed into a set
of six auxiliary equations that are numerically integrated via the Compound
Matrix Method. The temporal stability analysis suggest (a) elastic
stabilization at higher values of elasticity number (shown previously in the
dilute regime [S. Sircar and D. Bansal, “Spatiotemporal linear stability of
viscoelastic free shear flows: Dilute regime”, Phys. Fluids 31, 084104
(2019)]), (b) a non-monotonic instability pattern at low to intermediate
values of elasticity number for the JS as well as the PTT model. To comprehend
the effect of elasticity, Reynolds number and viscosity on the temporal
stability curves of the PTT model, we consider a fourth parameter, the
centerline shear rate, $\zeta_{c}$. The ‘JS behaviour’ is recovered below a
critical value of $\zeta_{c}$ and above this critical value the PTT base
stresses (relative to the JS model) is attenuated thereby explaining the
stabilizing influence of elasticity. The Briggs idea of analytic continuation
is deployed to classify regions of temporal stability, absolute and convective
instabilities, as well as evanescent modes, and the results are compared with
previously conducted experiments for Newtonian as well as viscoelastic flows
past a cylinder. The phase diagrams reveal the two familiar regions of
inertial turbulence modified by elasticity and elastic turbulence as well as
(a recently substantiated) region of elastoinertial turbulence and the
unfamiliar temporally stable region for intermediate values of Reynolds and
Elasticity number.

Keywords: Phan-Thien and Tanner model, Johnson-Segalman model, Orr-Sommerfeld
equation, Spatiotemporal stability, Compound Matrix Method

##  1 Introduction

The understanding of the hydrodynamic stability and the flow transition in
free shear flows of viscoelastic liquids has continued to receive prolonged
interest due to its practical applications in microfluidic mixing [1],
viscoelastic stabilization via polymer addition [2] and in shearing flows of
viscoelastic biofluids, like mucus [3], cartilage [4] and adhesion-
fragmentation transition in cells [5, 6, 7, 8]. Most investigations of low to
moderate Reynolds number ($Re$) have either focused on the linearized studies
[9], experimental analysis [10] or the full Direct Numerical Simulations (DNS)
[11] of the affine response of dilute polymeric liquids. However, in semi-
dilute or moderately concentrated polymeric liquids, instabilities may arise
due to flow induced inhomogeneities, and an improved understanding of the
transition pathway of the non-affine or the non monotonic flow response in
(but not limited to) strong elongational flow, is expedient.

Using the local spatiotemporal stability analysis where the base flow is
assumed to be locally parallel [12], we search for absolute instabilities
(perturbations which grow exponentially in time at the point of excitation),
convective instabilities (disturbances which are swept downstream from the
source and decay at any fixed position in space) and the evanescent modes (or
false modes) [13]. The local approach is best-suited for moderate to high
$Re$, especially for slowly spreading shear flows. Flow-spreading effects can
be fully captured with global methods though the computational cost is
substantially higher. Juniper _et al._ [14] discusses the complementary nature
of local and global methods in their study of the stability of confined
Newtonian wakes.

The concepts of absolute and convective instabilities are well established to
study the evolution of impulse disturbances in Newtonian flows. Absolute
instability was earlier experimentally verified by Shoji _et. al_ [15] in
liquid jets. There is also an extensive literature on Newtonian wakes and
mixing layers, including the blunt body experiments listed by Oertel [16], the
linear analysis by Chomaz [17] and Healy [18] and the DNS studies by Pier
[19]. Overall, these studies found a destabilizing effect of finite boundaries
and an existence as well as a transition to absolute instability in the near
wake region.

The spatiotemporal analysis of viscoelastic free shear layers are more recent
and scarce. An early experimental study by Vihinen _et. al_ [20] reported
absolute instability in viscoelastic liquid jets. Pipe reported a stabilizing
effect of polymer addition in his experiments on viscoelastic cylindrical
wakes, which is counteracted by shear thinning and a transition from
convective to absolute instability at higher polymer concentrations [21]. In
contrast, the linear analysis of dilute mixing layers [9] and dilute jets [22,
23, 24] relay a significant range of parameters where viscoelasticity was
found to be destabilizing. A recent DNS study of jets [25] necessitates the
use of an extra (convective) timescale to characterize the memory fading
property of viscoelastic fluids.

We have limited our focus on the linear spatiotemporal analyses of
viscoelastic free shear layer, mainly within the regime of non-affine / non
monotonic flow response. We characterize the fluid flows as being absolutely
or convectively unstable and utilize the method of spatiotemporal analysis by
progressive moving of the isocontours in the complex frequency and wavenumber
plane, as proposed by Kupfer [26]. The next section describes the model of the
viscoelastic free shear layer flow coupled with four constitutive relations
for the extra elastic stress tensor for simultaneous comparison of results:
Oldroyd-B, Upper Convected Maxwell (UCM), Johnson-Segalman (JS) at the slip
parameter value, $a=0.5$, and the linear Phan-Thien and Tanner (PTT) model at
$a=0.5$ and the elongation parameter explored within the range
$\varepsilon\in[0,0.99]$, although a detailed study is predominantly conducted
at $\varepsilon=0.5$ (§2.1) as well as the details of the linear stability
analysis leading upto the fourth order Orr-Sommerfeld equation (OSE) (§2.2).
§3 introduces the Compound Matrix Method (CMM) to numerically solve the
resultant system of stiff ordinary differential equations (ODE) emerging from
the OSE, rewritten in terms of the auxiliary variables. §4 showcases the
simulation results including numerical method validation (§4.1), the temporal
(§4.2) and the spatiotemporal stability analyses (§4.3), followed with a brief
discussion on the implication of these results as well as the focus of our
future direction (§5). The appendix (§A) lists all the coefficients of the
OSE.

##  2 Mathematical model and linear stability analysis

Unlike Newtonian solvents, the transition to instability in polymeric liquids
depends on the details of the equations relating stress to the shear rate. The
linear PTT model properly describes the non-affine behavior of polymeric
liquids composed of a low to moderate concentration of high molecular weight
polymer in a very viscous Newtonian solvent at moderate shear rates, including
the effect of shear thinning, non-zero normal stress coefficient and stress
overshoot in transient flows [11]. Further, this model predicts finite stress
at finite strain rates for strong elongation flow, a feature which is achieved
by constraining the length of the polymer chain to a maximum allowable length.
The numerical solution obtained by linearizing the Navier Stokes along with
the extra stress constitutive equation is reliable in describing the
characteristics of the initial stages of the mixing layer transition [27].
Hence, the linear stability results can be utilized as the initial-boundary
conditions for Direct Numerical Simulation [11] or the Large Eddy Simulation
of polymeric liquids [28] in order to further study the flow evolving process
downstream.

###  2.1 Mathematical model

The continuity and the momentum equations for an incompressible flow in a free
shear flow configuration are as follows,

$\displaystyle\nabla\cdot{\bf v}$ $\displaystyle=0,$ (1a)
$\displaystyle\rho\frac{D{\bf v}}{Dt}$ $\displaystyle=-\nabla
p+\nabla\cdot\tau,$ (1b)

where ${\bf v}$ is the velocity vector, $\rho$ the density, $p$ the isotropic
pressure and $\tau$ the extra stress tensor, which is written as the sum of
the viscous, Newtonian stress, $\tau^{s}$ (= $\eta_{s}{\bf D}$, where ${\bf
D}=\nabla{\bf v}+\nabla{\bf v}^{\text{T}}$ is the rate of strain tensor) and
the elastic stress, $\tau^{p}$ (= $\eta_{p}{\bf A}$). $\eta_{s},\eta_{p}$ are
the solvent viscosity and the polymeric contribution to the shear viscosity,
respectively. Introducing the parameters, $\eta(=\eta_{s}+\eta_{p})$ and
$\nu(=\eta_{p}/\eta)$, representing the total viscosity and the elastic
contribution to the total viscosity, respectively, the extra stress tensor can
be rephrased as follows,

$\tau=\eta[(1-\nu){\bf D}+\nu{\bf A}].$ (2)

The tensor ${\bf A}$ satisfies the linear PTT equation [11],

$\left(1+\varepsilon\lambda\text{tr}{\bf A}\right){\bf
A}+\lambda\frac{DA}{Dt}={\bf D},$ (3)

where the Gordon-Schowalter (GS) convected derivative,
$\displaystyle\frac{DA}{Dt}$, is,

$\frac{DA}{Dt}=\frac{\partial{\bf A}}{\partial t}+{\bf v}\cdot\nabla{\bf
A}-\frac{1}{2}({\bf W}^{\text{T}}+a{\bf D})\cdot{\bf A}-\frac{1}{2}{\bf
A}\cdot({\bf W}+a{\bf D}).$ (4)

${\bf W}=(\nabla{\bf v}\\!\\!-\\!\\!\nabla{\bf v}^{\text{T}})$ is the
vorticity tensor, $\lambda$ is the polymer relaxation time, $a$ is a slip
parameter characterizing the non-affine motion of the chains and
$\varepsilon\in[0,1]$ is a dimensionless parameter describing the maximum
elongation of the polymer chains. For the parameters, $\varepsilon=0,a=1$, the
motion is affine and the PTT model reduces to the Oldroyd-B model. The
Oldroyd-B model describes well the behavior of dilute polymeric liquids
composed of a low concentration of high molecular weight polymer in a very
viscous Newtonian solvent at moderate shear rates. The model predicts no shear
thinning, a constant first normal stress coefficient and a zero second normal
stress coefficient. But this constitutive relation fails in many
circumstances, e. g., it fails to predict the physical value of the viscosity
in extensional flow when the relaxation time of the polymer times the
extension rate exceeds 0.5 (or when the elongational flow is strong enough to
drive the two sides of the polymer dumbbell infinitely far apart from one
another) [27].

An particular case of the Oldroyd-B model is the UCM model
($\varepsilon=0,a=1,\nu=1$), a case where the extra stress tensor is purely
elastic. A UCM liquid predicts a quadratic shear rate dependence of the first
normal stress difference and zero second normal stress difference (which is a
realistic behavior of polymer melts at moderated shear rates [27]), but a
constant shear viscosity (or no shear thinning behavior).

The JS model (illustrated in detail at a fixed value of slip parameter,
$a=0.5$, and $\varepsilon=0$) allows for a non-monotonic relationship between
the shear stress and rate of shear in a simple shear flow, consequently
explaining the ‘spurt’ phenomena or a dramatic increase of the volumetric flow
rate (equivalently a spike in the strain rate) at a critical stress which is
independent of molecular weight of the polymer [29].

Finally, the linear PTT model (illustrated via the model parameters
$a=0.5,\varepsilon\in[0,0.99]$) is based on finitely extensible springs; it
does reproduce shear-thinning and generally captures the instability
transition more accurately than the Oldroyd-B as well as the JS model [11]. In
particular, this model is useful in obtaining a closer fit to the real world
normal stress differences including the study of contraction and re-entrant
corner flows [30].

###  2.2 Linear Stability Analysis

Consider the free-stream velocity in a frame moving with the average flow
velocity, i. e., $U$ ($=\frac{1}{2}$($\mathcal{U}_{1}$ \- $\mathcal{U}_{2}$),
where $\mathcal{U}_{1}$($\mathcal{U}_{2}$) is the free-stream velocity of the
upper(lower) flow), and the momentum thickness, $\delta$ [31]. Utilizing
$U,\delta$ as the reference velocity and the length scale, respectively, we
nondimensionalize equation (1) and introduce the the dimensionless numbers,
$Re=\nicefrac{{\rho\delta U}}{{\eta}}$ and the Weissenberg number,
$We=\nicefrac{{\lambda U}}{{\delta}}$. Assuming that the mean flow is two-
dimensional (with $x$ and $y$ being the streamwise parallel and transverse
directions in space, respectively) and quasi-parallel with its variation
entirely in the transverse direction, i. e.,

$U(y)=\tanh(y),\,\,\Omega(y)=\tanh^{2}(y)-1,\,\,\Psi(y)=\log(\cosh(y)),$ (5)

where $U(y),\Omega(y),\Psi(y)$ are the dimensionless streamwise parallel mean
velocity (with zero transverse component), mean vorticity and the associated
streamfunction, respectively. Further, assuming that the mean flow supports a
two-dimensional disturbance field, the streamfunction and the extra stress
tensor are represented by the base state profile $(\Psi(y),\tau_{0}(y))$ plus
a small perturbation, which is Fourier transformed in $x$ and $t$ as follows,

$\displaystyle\psi(x,y,t)$ $\displaystyle=\Psi(y)+\phi(y)e^{{\it i}(\alpha
x-\omega t)},$ $\displaystyle\tau(x,y,t)$
$\displaystyle=\tau_{0}(y)+\varphi(y)e^{{\it i}(\alpha x-\omega t)},$ (6)

where $\phi(y),\varphi(y)$ are the transverse perturbations in the
streamfunction and the extra stress tensor and $\alpha,\omega$ are the complex
wavenumber and angular frequency, respectively. We note that equation (5) is a
solution of the momentum equations for incompressible flow provided there is a
dimensionless body force term on the right-hand-side of equation (1b) [32]. We
rephrase equations (1-4) in the streamfunction-vorticity formulation and avail
equation (6) to arrive at the equation governing the perturbation of the
streamfunction, given by the fourth order OSE [31].

$\left\\{\\!\\!i\left[\\!(\alpha
U\\!\\!-\\!\\!\omega)(D^{2}\\!\\!-\\!\\!\alpha^{2})\\!\\!-\\!\\!\alpha
D^{2}U\\!\right]\\!\\!-\\!\\!\frac{1-\nu}{Re}(D^{2}\\!\\!-\\!\\!\alpha^{2})^{2}\\!\\!\right\\}\phi\\!\\!=\frac{\nu}{\mathcal{F}Re}\\!\\!\sum_{n=0}^{4}\\!\\!c_{n}D^{n}\phi.$
(7)

where $D^{n}(\cdot)$ denote the n${}^{\text{t}h}$ derivative of any variable
function with respect to $y$ and the coefficients $c_{i}$’s (equations (20))
are listed in §A. $D=\nicefrac{{d}}{{dy}}$,
$\mathcal{F}=4d^{3}_{0}\left[aS_{1}S_{2}+S_{1}S_{3}+(1-a^{2})S^{2}_{1}\right]^{2}\left[S_{3}+(1+a)S_{1}\right]^{2}$
(refer §A). Following the usual Newtonian development [31], we assume that
$\phi=\phi_{r}+{\it i}\phi_{i}$ (the subscript $r/i$ from this point onwards
will denote the real and the imaginary components) and $\phi_{r}$ and
$\phi_{i}$ (both real valued functions) are even and odd functions of $y$
respectively. Thereby, we restrict the domain of integration to the upper half
of the flow such that the boundary conditions at $y=0$ (or the so called
centerline conditions), are altered as follows,

$\displaystyle\phi_{i}=\phi^{\prime\prime}_{i}=0$ (8a)
$\displaystyle\phi^{\prime}_{r}=\phi^{\prime\prime\prime}_{r}=0.$ (8b)

The far stream boundary conditions (or the conditions at $y\rightarrow\infty$)
is given by

$\phi_{r/i}=\phi^{\prime}_{r/i}=0.$ (9)

In § 3, we detail a solution procedure for solving the eigenvalue problem
(7-9).

##  3 Solution to the eigenvalue problem

The primary step in evaluating the eigenvalues ($\alpha,\omega$) and the
eigenfunction, $\phi$, describing the disturbance field (satisfying equation
(7)) is to probe the ramification of the far stream boundary conditions (9) by
using the mean flow information at far stream, i. e., $\lim
y\rightarrow\infty:U(y)=1,U^{\prime\prime}(y)=0$, which simplifies the OSE (7)
to the following constant coefficient ODE [33],

$\phi^{(4)}-2\alpha^{2}\phi^{\prime\prime}+\alpha^{4}\phi=\frac{{\it
i}ReS_{\infty}}{\nu+(1-\nu)S_{\infty}}(\alpha-\omega)(\phi^{\prime\prime}-\alpha^{2}\phi),$
(10)

where $S_{\infty}=1+We(\alpha-\omega)$. The solution to equation (10) can be
derived by setting $\phi=e^{\lambda y}$, such that one gets the characteristic
roots as $\lambda_{1,2}=\mp\alpha$ and $\lambda_{3,4}=\mp Q$, where
$Q\\!\\!=\\!\\!{\displaystyle\left[\alpha^{2}\\!\\!+\\!\\!\frac{{\it
i}ReS_{\infty}(\alpha\\!\\!-\\!\\!\omega)}{\nu\\!\\!+\\!\\!(1-\nu)S_{\infty}}\\!\right]^{\frac{1}{2}}}$.
The fourth order OSE (7) has four fundamental solutions, i. e.
$\\{\phi_{i}\\}^{4}_{i=1}$, whose asymptotic variation with
$y\rightarrow\infty$ is: $\phi_{1,2}\sim e^{\mp\alpha y}$; $\phi_{3,4}\sim
e^{\mp Qy}$. Then the general solution which satisfies the far stream
conditions (9) for real $(\alpha,Q)>0$, is of the form

$\phi=a_{1}\phi_{1}+a_{3}\phi_{3}.$ (11)

Equation (11) admits a non-trivial solution of the OSE, satisfying the far
stream condition (9) and the centerline conditions (8a, 8b) if and only if the
determinant of the associated matrix of the linear algebraic system given by
equation (11), vanishes at $y=0$, or

$\left(\phi_{1}\phi^{\prime\prime}_{3}-\phi^{\prime\prime}_{1}\phi_{3}\right)|_{y=0}=0,$
(12)

for the odd component (equation (8a)) and

$\left(\phi^{\prime}_{1}\phi^{\prime\prime\prime}_{3}-\phi^{\prime\prime\prime}_{1}\phi^{\prime}_{3}\right)|_{y=0}=0,$
(13)

for the even component (equation (8b)), respectively. Equations (12, 13) serve
as the dispersion relation of the problem and will be solved simultaneously.
The stiffness of the OSE (7) (e. g., in the case of far stream eigenmodes in
the limit $Re\rightarrow\infty$, we see that $|Q|\gg|\alpha|$, leading to an
immense contrast between the two sets of characteristic roots of equation
(10)) thereby causing the solution components corresponding to the different
fundamental solutions to lose linear independence. This source of parasitic
error growth necessitates the use of CMM [34], where one works with a set of
following auxiliary variables,

$\displaystyle y_{1}=\phi_{1}\phi^{\prime}_{3}-\phi_{3}\phi^{\prime}_{1},\quad
y_{2}=\phi_{1}\phi^{\prime\prime}_{3}-\phi_{3}\phi^{\prime\prime}_{1},\quad
y_{3}=\phi_{1}\phi^{\prime\prime\prime}_{3}-\phi_{3}\phi^{\prime\prime\prime}_{1},$
$\displaystyle
y_{4}=\phi^{\prime}_{1}\phi^{\prime\prime}_{3}-\phi^{\prime\prime}_{1}\phi^{\prime}_{3},\quad
y_{5}=\phi^{\prime}_{1}\phi^{\prime\prime\prime}_{3}-\phi^{\prime\prime\prime}_{1}\phi^{\prime}_{3},\quad
y_{6}=\phi^{\prime\prime}_{1}\phi^{\prime\prime\prime}_{3}-\phi^{\prime\prime\prime}_{1}\phi^{\prime\prime}_{3},$
(14)

satisfying the initial value problem (IVP) [32],

$\displaystyle y^{\prime}_{1}\\!\\!$ $\displaystyle=\\!\\!y_{2},$
$\displaystyle y^{\prime}_{2}\\!\\!$ $\displaystyle=\\!\\!y_{3}+y_{4},$
$\displaystyle y^{\prime}_{3}\\!\\!$
$\displaystyle=\\!\\!y_{5}\\!\\!-\\!\\!\left[\frac{(1\\!\\!-\\!\\!\nu)(c_{1}y_{1}\\!\\!+\\!\\!c_{2}y_{2}\\!\\!+\\!\\!c_{3}y_{3})\\!\\!-\\!\\!\mathcal{F}({\it
i}Re(\alpha
U\\!\\!-\\!\\!\omega\\!\\!+\\!\\!2\alpha^{2}\nu)y_{2}}{(1\\!-\\!\nu)c_{4}\\!+\\!\mathcal{F}\nu}\right],$
$\displaystyle y^{\prime}_{4}\\!\\!$ $\displaystyle=\\!\\!y_{5},$
$\displaystyle y^{\prime}_{5}\\!\\!$
$\displaystyle=\\!\\!y_{6}\\!\\!+\\!\\!\left[\frac{(1\\!\\!-\\!\\!\nu)(c_{0}y_{1}\\!\\!-\\!\\!c_{2}y_{4}\\!\\!-\\!\\!c_{3}y_{5})\\!\\!+\\!\\!\mathcal{F}({\it
i}Re(\alpha
U\\!\\!-\\!\\!\omega)\\!\\!+\\!\\!2\alpha^{2}\nu)y_{4}}{(1\\!-\\!\nu)c_{4}\\!+\\!\mathcal{F}\nu}\right.$
$\displaystyle\left.+\frac{\mathcal{F}({\it i}Re[(\alpha
U\\!\\!-\\!\\!\omega)\alpha^{2}\\!\\!+\\!\\!\alpha
U^{\prime\prime}]\\!\\!+\\!\\!\alpha^{4}\nu)y_{1}}{(1\\!-\\!\nu)c_{4}\\!+\\!\mathcal{F}\nu}\right],$
$\displaystyle y^{\prime}_{6}\\!\\!$
$\displaystyle=\\!\\!\frac{(1\\!\\!-\\!\\!\nu)(c_{0}y_{2}\\!\\!+\\!\\!c_{1}y_{4}\\!\\!-\\!\\!c_{3}y_{6})\\!\\!+\\!\\!\mathcal{F}\\!({\it
i}Re[(\alpha U\\!\\!-\\!\omega)\alpha^{2}\\!\\!+\\!\\!\alpha
U^{\prime\prime}])y_{2}}{(1\\!-\\!\nu)c_{4}\\!+\\!\mathcal{F}\nu}$
$\displaystyle+\frac{\mathcal{F}\alpha^{4}\nu
y_{2}}{(1\\!-\\!\nu)c_{4}\\!+\\!\mathcal{F}\nu},$ (15)

where the initial conditions are estimated by substituting the free stream
values of the unknown (i. e., in the limit $y\rightarrow\infty$ substitute
$\phi_{1}\sim e^{-\alpha y}$ and $\phi_{3}\sim e^{-Qy}$ in equation (14)) and
normalizing with respect to one of the variables (e. g., $y_{1}$) to remove
stiffness. The rescaled initial conditions for solving equations (15) are

$\displaystyle y_{1}=1.0,\quad y_{2}=-(\alpha+q),\quad
y_{3}=\alpha^{2}+q\alpha+q^{2},$ $\displaystyle y_{4}=q\alpha,\quad
y_{5}=-q\alpha(\alpha+q),\quad y_{6}=(q\alpha)^{2}$ (16)

The numerical solution for IVP (15, 16) is obtained by marching backward from
the free stream to the centerline. A suitable value of the eigenpair
($\alpha,\omega$) is obtained by enforcing the dispersion relation (equations
(12, 13)) in auxiliary variables and solving simultaneously, i. e.,

$\displaystyle\text{Re}(y_{2})=0\quad\text{at}\quad y=0\quad\text{for odd
component},$ (17a) $\displaystyle\text{Re}(y_{5})=0\quad\text{at}\quad
y=0\quad\text{for even component},$ (17b)

respectively. $\text{Re}(\cdot)$ denotes the real part of the complex valued
function.

##  4 Results

The zeros of the dispersion relation (equation (17)) were explored within the
complex $\alpha-\omega$ plane inside the region $-0.02\leq\omega_{r}\leq
0.13,-1.6\leq\omega_{i}\leq 0.25,\alpha_{r}\leq 1.6$ and $|\alpha_{i}|\leq
0.02$. Previous results indicate that the influence of viscoelasticity is
fully captured by the modified elasticity number, $E=\frac{\nu We}{Re}$, a
parameter representing the ratio of the fluid relaxation time to the
characteristic time for vorticity diffusion [9]. We highlight our instability
results versus this parameter. The continuation curves in §4.2 and 4.3 are
depicted within the range $E\in[10^{-3},\,\,5.0]$ using a discrete step-size
of $\triangle E=10^{-3}$. The numerical solution of the IVP (15, 16) was
determined via the fourth order Runge Kutta integration with a step-size of
$\triangle y=2.2\times 10^{-3}$. The integration domain was truncated at
$\eta=12.0$ (a point at which the free stream boundary conditions (9) were
imposed), which leads to a value of the momentum thickness (the reference
length scale introduced in §2.2), $\delta=0.30685$. The results in §4.2 and
§4.3 are compared at two different values of $Re$ (i. e., $Re=40$ and
$Re=400$) as well as, at $\nu=0.3$ (the viscous stress dominated case) and at
$\nu=0.7$ (the elastic stress dominated case).

###  4.1 Numerical method validation

First, our numerical method outlined in §3 is validated by reproducing the
absolute and the convective instability results for (a) inviscid mixing layers
investigated by Huerre and Monkewitz [12] or the Rayleigh instability equation
(figure 1a), and (b) spatially developing viscoelastic Oldroyd-B mixing layers
probed by Ray and Zaki [9] (figure 1b), using the non-dimensional base-state
velocity profile given by

$U(y)=1+S\tanh\left(\frac{y}{2}\right),$ (18)

where the parameter $S$ is the ratio of the difference and the sum of the
free-stream velocities of the upper and the lower half. In figure 1a, we
recover the familiar curve of the cusp / pinch points for different complex
pairs $(\alpha,\omega)$ determined by the numerical integration of the
inviscid Rayleigh equation together with the exponentially decaying far stream
boundary condition, equation (9). Note the crossover from the real-$\omega$
axis at the critical value, $S=1.315$, highlighting a transition from
convective instability (i. e., $\omega_{i}^{\text{cusp}}<0$) to absolute
instability at this value. Figure 1b presents the evolution of the critical
value of $S$ versus the elasticity number, $E$, at $\nu=0.5$ and
$Re=50,100,400,1000$ for an Oldroyd-B fluid. First, notice the reduction at
low elasticity number regime, $E<0.1$, followed by an enhancement of absolute
instability, with increasing $E$ (i. e., the critical value of $S$ eventually
drops with increasing $E$). Second, note that the elasticity number is the
dominant flow parameter measuring viscoelasticity while the Reynolds number
has negligible influence beyond $Re\geq 400$. A more detailed outlook of the
influence of viscoelasticity is acquired by examining the temporal growth
rates, described in the next section.

Figure 1: (a) Locus of the cusp points of the Rayleigh equation in
$\omega$-plane at critical values of the parameter, $S$, and (b) the critical
values of $S$ versus elasticity number, $E$, of the OSE (equation (7)) for
Oldroyd-B fluids, at $\nu=0.5$ and $Re=50,100,400,1000$. The mean velocity
profile is given by equation (18).

###  4.2 Temporal stability analysis

A positive sign of the temporal growth rate indicates whether absolute
instability is possible. The temporal stability analysis for viscoelastic free
shear flows in the dilute regime, for low to moderate $Re$ and $We$, was
earlier studied by us [32]. In the limit of large $Re$ and $We$ such that
$We/Re\sim\mathcal{O}(1)$, Azaiez conducted the temporal stability analysis
through an elastic Rayleigh equation and concluded elasticity as the
controlling flow parameter within the dilute flow regime [31]. In this study,
we extend these ideas within the non-affine / non-monotonic response regime.
Figure 2 presents the solution of the OSE (7) for purely real wavenumbers,
$\alpha_{r}$ (third column), while allowing angular frequency to be complex
number. The other two columns in figure 2 are the temporal growth rate or the
largest positive imaginary component of any root of the dispersion relation
(17), $\omega_{i}^{\text{Temp}}$ (second column) and the corresponding real
part of the frequency (first column). The four models are shown using solid
line (Oldroyd-B fluids), dotted line (UCM fluids), dash-dot line (JS fluids,
$a=0.5$) and dashed line (PTT fluids, $a=0.5,\varepsilon=0.5$).

First, notice the instability curves (in particular the curves for the UCM
fluid) do not start at the same point near $E=0$. Readers are reminded that
the simulation starts at $E_{0}=10^{-3}$, i. e., $E\equiv 0$ (the Newtonian
case) is not considered in the present discussion. Further note that the
elastic stress dominated case at $Re=40$, with the exception of the UCM fluid,
is unstable in the limit $E\rightarrow 0$ (figure 2h). This instability occurs
at short wavelength (i. e., large $\alpha_{r}$, figure 2i) and low frequency
(i. e., small $|\omega_{r}|$, figure 2g). In the dilute regime, the shear flow
instabilities were found to arise at zero elasticity number, through a
combination of instability via normal stress anisotropy and elasticity [32].
In the non-affine regime, we surmise a similar operative mechanism. The
overall trend, for large values of $E$, is that elasticity is stabilizing for
UCM fluids (figures 2b,e,h,k), stabilizing for PTT fluids at higher $Re$
(figures 2e,k) but has negligible small (although stabilizing) influence on
Oldroyd-B and JS fluids ($\omega_{i}^{\text{Temp}}$ curve is nearly flat in
the range $E\geq 3.5$). The mechanism of this elasticity-induced stabilization
is well documented by Hinch in an appendix to Azaiez [31], and is akin to the
action of a ‘surface tension’. The stretched polymers contribute to an
effective tension along the vibrating membrane that is the shear layer, and
this tension damps the perturbations. This important analogy with the surface
tension helps to provide a physical explanation for the influence of
viscoelasticity on mixing layer stability, in the large-$E$ limit.

The JS as well as the PTT fluids at $Re=40$ shows a non-monotonic behavior at
low to intermediate values of $E$ (figures 2b,h) with alternating regimes of
stability followed by instability, a feature typical of semi-dilute or
moderately concentrated polymeric liquids. These fluctuations are
characterized at short wavelengths (figures 2c,i) and low frequency (figures
2a,g). However, the PTT model ‘overestimates’ (‘underestimates’) the
instability predicted by the JS model for viscous stress dominated (elastic
stress dominated) case at low to intermediate values of $E$ (i. e., compare
dashed and the dash-dot curves in figures 2b,e versus figures 2h,k). In
addition, within the range $E\geq 3.5$, while the stability of the Oldroyd-B
as well as the JS fluid is almost exclusively a function of $E$ (figures
2b,e,h,k), the PTT model produces a more complicated behavior, dependent on
all the three parameters, $E,Re,\nu$. The reason for these observations is
explained as follows. The JS model is derived from a kinetic theory in which
the polymers are represented as beads connected by springs and this underlying
assumption of infinitely extensible springs is limiting in the JS model. The
finitely extensible (nonlinear) springs of the PTT model arrests the infinite
stresses at finite strain rate in extensional flows, leading to a complex
pattern described above.

Figure 2: The frequency, $\omega_{r}$, (first column); most unstable mode,
$\omega^{\text{Temp}}_{i}$, (second column); wavenumber, $\alpha_{r}$ (third
column) for parameters (a-c) $\nu=0.3$, $Re=40$, (d-f) $\nu=0.3$, $Re=400$,
(g-i) $\nu=0.7$, $Re=40$ and (j-l) $\nu=0.7$, $Re=400$ versus the elasticity
number, $E$, using using solid, dotted, dash-dot and dashed curves for
Oldroyd-B, UCM, JS (at $a=0.5$) and PTT fluids (at $a=0.5,\varepsilon=0.5$),
respectively.

In particular, the PTT results (compared with the JS results at
$\nu=0.7,Re=400$, figure 3a) indicate that a realistic value of the maximum
elongation, $\varepsilon$, can induce a significant stabilizing influence of
elasticity for large values of $E$. For small values of $\varepsilon$ (i. e.,
$\varepsilon\leq 0.1$), results are close to those obtained with the JS model.
As $\varepsilon$ is increased, the (nonlinear) stiffness of the polymer
molecules is increased, and one expects a reduction in the (stabilizing)
influence of elasticity. We have considered several values of $\varepsilon$,
but are particularly interested in results for $\varepsilon=0.5$. Previous
analysis of viscoelastic mixing layers in high Weissenberg number flows have
used values between $\varepsilon=0.3$ and $\varepsilon=0.7$ [35].

We gain a better understanding of these results and construct a criteria for
the recovery of the ‘JS behavior’ by examining the base-state polymer
stresses, $\tau_{0}$ (equation (24)). In steady homogeneous shear flow, the
behavior of the base-state polymer stress in the PTT mixing layer is closely
related to the standard viscometric results which, in succession, depends on
the non-dimensional shear rate [36]. At low shear rates, the elongation of the
PTT springs is modest, and the PTT and JS models produce similar results. At
higher shear rates, the springs become more stretched, nonlinear elasticity
becomes important, and the PTT polymer stress components become attenuated
relative to their JS counterparts. These trends are present in the tanh mixing
layer as well, though the velocity gradient is not constant, and the polymer
stress depends on the local shear rate. The role of elongation must also be
considered, and using Eqs. (24), the ratio of the JS and PTT stress components
can be shown to be functions of the shear rate, $\zeta\equiv\varepsilon
EU_{y}$. Figure 3b presents the base stresses at $\nu=0.7,Re=400$, with
respect to this shear rate at the centerline, $\zeta=\zeta_{c}=\varepsilon E$.
Unsurprisingly, when the shear rate, $\zeta_{c}$, is small, there is little
difference between the PTT and the JS models, but at larger values of
$\zeta_{c}$, the JS stress is larger. We can foresee that the relative
attenuation of the base stresses as $\zeta_{c}$ is increased, is connected to
the stabilizing influence of elasticity which was observed in our temporal
stability results. Conversely, we would expect to recover the JS behavior if
$\zeta_{c}$ is less than some critical value, $\zeta^{*}_{c}=0.126$, i. e.,
the ratio of all the base stresses are less than 90% beyond this critical
value of $\zeta_{c}$. A better characterization of these instabilities are
revealed through the spatiotemporal analysis in §4.3.

Figure 3: (a) The temporal growth rate for linear PTT fluids,
$\omega_{i}^{\text{Temp}}$ versus $E$, and (b) the log-log plot of the ratio
of the PTT versus JS base stresses as a function of the centerline shear rate.
Other parameters are fixed at $a=0.5,\nu=0.7,Re=400$ and $\varepsilon=0.5$ for
plot (b).

###  4.3 Spatiotemporal stability analysis

Spatiotemporal analysis is typically relevant when one introduces an impulse
excitation locally in a flow and observes how that disturbance evolves. In an
effort to determine the range of $E$ (for fixed $Re$ and $\nu$) for which the
flow regimes are absolutely unstable, convectively unstable or temporally
stable, we recover the absolute growth rate (or the growth rate at the cusp
point, $\omega^{\text{cusp}}_{i}$, figure 4) starting from the most unstable
temporal mode, $\omega^{\text{Temp}}_{i}$. The cusp point in the
$\omega-$plane is a saddle point satisfying the criteria,
$D(\alpha,\omega^{\text{c}usp})\\!=\\!\frac{\partial
D(\alpha,\omega^{\text{c}usp})}{\partial\alpha}\\!=\\!0$ but
$\frac{\partial^{2}D(\alpha^{\text{p}inch},\omega^{\text{c}usp})}{\partial\alpha^{2}}\\!\neq\\!0$
(where $D(\alpha,\omega)=0$ is the dispersion relation). However, not all cusp
points are unstable and, in particular, the evanescent modes are segregated
from the regular cusp points using the Briggs idea of analytic continuation
[26]. While the Oldroyd-B model is presumed to represent the dilute polymeric
liquids, the linear PTT model (due to the finite, nonlinear elongation of the
polymer chains at $a=0.5,\varepsilon=0.5$) portrays the instability transition
for moderately concentrated polymeric liquids. A discontinuity in the curves
in figure 4 indicates a region of temporal stability.

The instability pattern shown in figure 4 is a result of a complex interplay
between the inertial forces (proportional to $Re$) and the normal stress
anisotropy through elasticity (proportional to $E$). For example, within the
lower elastic number regime (i. .e., $E<0.5$) the dilute polymeric liquids
display constricted regions of temporal stability at higher value of $Re$
(comparing the solid curve at $Re=40$ (figure 4a) versus the curve at $Re=400$
(figure 4b)). Similarly, the elastic stress dominated case reveals
(convective) instability at lower value of $Re$ and $E$ (figure 4c). Clearly,
while the former observation is the result of inertia, the latter is the
outcome of instability generated via polymer elasticity. Analogous with the
temporal stability analysis, we find that the absolute growth rate of the PTT
fluids in the high elastic number regime is always lower than the JS fluids.
This is because the finite elongation attribute of the PTT fluids is able to
rein in the large elastic stress gradient buildup.

To explore the nature of these instabilities, we compute the boundaries of the
temporally stable regions (S), convective instabilities (C), evanescent modes
(E) and absolute instabilities (A) within a selected range of flow-elasticity
parameter space, i. e., $Re\in[0.1,400],E\in[10^{-3},5]$ and $\nu=0.3$ (figure
5a) and $\nu=0.7$ (figure 5b). The boundaries of the temporally stable,
convectively unstable and absolutely unstable region, for dilute liquids in
the viscous stress dominated case (figure 5a) are numerically estimated to
reside within the range $E<0.3$, approximately $0.3\leq E\leq 1.5$ and $E>1.5$
(for all values of $Re$), respectively. In contrast, the concentrated
polymeric liquids exhibit temporal stability and absolute instability in a
confined region, i. e., $E<0.1$ and $E>4$, respectively. For the elastic
stress dominated case (figure 5b), the dilute as well as the concentrated
liquid reveals convective instability within the range $Re<150$ and $E<0.1$
(including in the limit $Re\ll 1$) which is followed with alternating regions
of (temporal) stability and (convective) instability for larger values of $E$.
Absolute instability is unveiled by the Oldroyd-B fluids for $E>1.5$ and by
the PTT fluids for $E>1.0$, for practically all values of $Re$.

Figure 4: The cusp point, $\omega^{\text{cusp}}_{i}$ versus $E$ evaluated at
the flow parameters, (a) $\nu=0.3,Re=40$, (b) $\nu=0.3,Re=400$, (c)
$\nu=0.7,Re=40$, and (d) $\nu=0.7,Re=400$, shown using solid, dotted, dash-dot
and dashed curves for Oldroyd-B, UCM, JS (at $a=0.5$) and PTT fluids (at
$a=0.5,\varepsilon=0.5$), respectively.

We discuss the significance of the temporal stability analysis (figure 2) as
well as the spatiotemporal phase diagram (figure 5) in relation to the
experiments for Newtonian as well as viscoelastic flow past a cylinder, since
the shear layer instability is closely related to the instability of the
cylinder wakes [21]. In the case of Newtonian wakes for $Re<1$, the flow past
a cylinder is steady and without recirculation [37]. The emergence of a
locally, convectively unstable wake is detected at $Re\approx 5$, whereby the
selected perturbation are amplified and convected downstream but ultimately
leave the flow undisturbed [38]. The locally most unstable part of the wake
becomes absolutely unstable at $Re\approx 25$, although this is not sufficient
for self-sustained global oscillations of the wake. At $Re\approx 47$ the
region of absolute instability in the wake is large enough for the wake to
sustain time-amplified oscillations, which is followed by the onset of the
laminar two-dimensional von Kármán instability [39]. This transition to global
instability has been ascribed to a supercritical Hopf bifurcation towards a
limit cycle [40] and the linear stability analysis appears to faithfully
describe the Newtonian wake dynamics, ‘unreasonably’ far above the global
instability threshold of $Re$ [16].

In contrast, first notice that all four stress constitutive relations in the
phase diagram in figure 5 indicates a convectively unstable region for
moderate values of $E$. For the viscous stress dominated case, figure 5a, this
region lies in the range $0.6\leq E\leq 0.7$ (Oldroyd-B), $E_{0}\leq E\leq
0.1$ (UCM), $0.5\leq E\leq 0.6$ (JS and PTT); and for the elastic stress
dominated case, figure 5b, this region lies in the range $0.8\leq E\leq 1.1$
(Oldroyd-B), $E_{0}\leq E\leq 0.1$ (UCM), $0.5\leq E\leq 0.6$ (JS and PTT);
for significantly larger values of $Re$ (i. e., $Re>47$). This observation is
in congruence with the early experimental studies of viscoelastic vortex
street highlighting the reduction in the vortex shedding frequency (related to
the temporal growth rate in the present study) [41], as well as in the
intensity of the vorticity [42], namely the origin of inertial turbulence
modified by elasticity. Second, notice the appearance of convectively unstable
region for the elastic stress dominated case, for small values of $E$ and in
the limit $Re\rightarrow 0$. This result corroborates the experimental
findings of McKinley _et al._ [43] showing instability induced by elasticity
for viscoelastic flow past a cylinder at $Re\ll 1$ as well as the findings of
Coelho and Pinho [44] showing a significant destabilization of the wake for
shear-thinning, elastic fluids at $Re<40$, in other words the presence of
elastic turbulence.

Two other observations are noteworthy: first, the Oldroyd-B, JS and the PTT
models suggest absolute instability for sufficiently large values of $E$ (i.
e., $E>4$ ($E>1.5$) for the viscous (elastic) stress dominated case, figure 5a
versus 5b) and second, the presence of a temporally stable region for the JS
and PTT models for the viscous stress dominated case (figure 5a),
approximately within the range $0.6<E<4$ and $Re<100$. While the former
observation signifies a pathway to elastoinertial turbulence (appearing at
moderate $Re$ and large $E$) which characterizes the maximum drag reduction
state and even originates in linear instability studies of pipe flows [45],
the latter observation indicates an intricate tug-of-war between the inertial
destabilization and elastic stabilization. While the presence of the first
region was recently substantiated [46], a comprehensive experimental
ratification of the second region is eagerly awaited.

Figure 5: Viscoelastic free shear layer stability phase diagram at (a)
$\nu=0.3$, and (b) $\nu=0.7$, in the $E$–$Re$ parametric space. The regions S,
C and A are denoted by temporally stable, convectively unstable and absolutely
unstable regions, respectively. The domains outline by (A, E) and (C, E) are
those where both the stable and evanescent modes (denoted by E) are found.

##  5 Conclusions

This investigation addresses the linear, temporal and the spatiotemporal
analyses of free shear flows of dilute as well as moderately concentrated
polymeric liquids for low to moderate Reynolds number and Elasticity number,
by capturing the non-affine flow response of a ‘tanh’ base flow mixing
velocity profile. Section 2 presented the viscoelastic free shear flow model
as well as the elements of the linear stability analysis via the solution of
the Orr-Sommerfeld equation. Section 3 demonstrated the steps of the Compound
Matrix Method, utilized to numerically solve the resultant system of stiff
differential equations. The temporal stability analysis in section 4.2,
indicates (a) elastic stabilization at higher values of elasticity number and
(b) a non-monotonic instability pattern at low to intermediate values of
elasticity number for the JS as well as the PTT model. The spatiotemporal
phase diagram in section 4.3 divulge the familiar regions of inertial and
elastic turbulence, a recently verified region of elastoinertial turbulence
and the unfamiliar temporally stable region for intermediate values of
Reynolds and Elasticity number.

Although this study provides an improved understanding of the linear dynamics
of mixing layers for dilute to moderately concentrated polymeric liquids, a
number of simplifying assumptions were made, and the relaxation of these
assumptions paves a way for further progress. Understanding the importance of
base flow spreading, confinement and nonlinearity as well as the consideration
of the shear flows of polymer melts (or fluids with very large viscosity) are
of substantial importance. Finally, we note that we have considered only the
two-dimensional instabilities which is a common simplification in the absolute
/ convective instability studies of Newtonian flows where Squire’s
transformation can be applied [1]. A modified Squire’s transformation for
Oldroyd-B fluids also exist [47], but we are not aware of a similar result for
the PTT model. Hence, a consideration of the three-dimensional instability
modes in future studies with the PTT model may be worthwhile. Acknowledgments
D. B. and S.S. acknowledges the financial support of the Grant CSIR 09 / 1117
(0004) / 2017-EMR-I and DST ECR / 2017 /000632, respectively.

##  Appendix A coefficients of OSE

The OSE, outlined in §2.2, coupled with the linear PTT stress constitutive
equation, is given by

$\left\\{i\big{[}(\alpha U-\omega)(D^{2}-\alpha^{2})-\alpha
D^{2}U\big{]}-\frac{1-\nu}{Re}(D^{2}-\alpha^{2})^{2}\right\\}\phi=\frac{\nu}{\mathcal{F}Re}\sum_{n=0}^{4}c_{n}D^{n}\phi,$
(19)

where the coefficients $c_{i}$’s are specified as follows

$\displaystyle c_{0}$
$\displaystyle=\frac{-4\left(S_{3}+(1+a)S_{1}\right)^{2}}{(4S_{0}^{2}d_{0}^{3}(S_{3}+(1+a)S_{1})^{2})}\left[\left[S_{0}^{2}\left\\{(D^{2}\tilde{A})d_{0}^{2}-\tilde{A}d_{0}(D^{2}d_{0})-2(D\tilde{A})d_{0}(Dd_{0})+2\tilde{A}(Dd_{0})^{2}+\alpha^{2}\tilde{A}d_{0}^{2}\right\\}\right.\right.$
$\displaystyle+2(1-a^{2})i\alpha\left\\{S_{0}\left((D^{2}U)(S_{1}+S_{3})+(DU)(DS_{1}+DS_{3})\right)-(DU)(S_{1}+S_{3})(DS_{0})\right\\}(\tilde{A})d_{0}^{2}+2(1-$
$\displaystyle\left.a^{2})i\alpha(DU)(S_{1}+S_{3})S_{0}\left\\{(D\tilde{A})d_{0}^{2}-2\tilde{A}d_{0}(Dd_{0})\right\\}\right]+d_{0}^{3}\left[2S_{0}(S_{3}+(1+a)S_{1})\left\\{2(S_{1}+S_{3})(S_{2}-(1+\right.\right.$
$\displaystyle\left.a)S_{1})+1+2(S_{3}+(1+a)S_{1})(S_{1}+S_{3})\right\\}i\alpha(1+a)(D^{2}\tau_{0}^{22})+\left\\{(2(S_{1}+S_{3})(S_{2}-(1+a)S_{1})+1)-\right.$
$\displaystyle\left.2(S_{1}+S_{3})(S_{3}+(1+a)S_{1})\right\\}i\alpha(1-a)(D^{2}\tau_{0}^{11})+\left\\{(2(S_{1}+S_{3})(S_{2}-(1+a)S_{1})+1)a+(S_{1}+S_{3})\right.$
$\displaystyle\left.(S_{3}+(1+a)S_{1})(1+a^{2})\right\\}4\alpha^{2}(D\tau_{0}^{12})+2\left\\{2S_{0}(S_{3}+(1+a)S_{1})((DS_{1}+DS_{3})(S_{2}-(1+a)S_{1})+\right.$
$\displaystyle(S_{1}+S_{3})(DS_{2}-(1+a)DS_{1}))-(2(S_{1}+S_{3})(S_{2}-(1+a)S_{1})+1)((DS_{0})(S_{3}+(1+a)S_{1})+S_{0}$
$\displaystyle\left.(DS_{3}+(1+a)DS_{1}))\right\\}(i\alpha((1+a)D\tau_{0}^{22}+(1-a)D\tau_{0}^{11})+4a\alpha^{2}\tau_{0}^{12})+4\left\\{S_{0}(S_{3}+(1+a)S_{1})^{2}\right.$
$\displaystyle\left.(DS_{1}+DS_{3})-(S_{3}+(1+a)S_{1})(S_{1}+S_{3})((DS_{0})(S_{3}+(1+a)S_{1})+2S_{0}(DS_{3}+(1+a)DS_{1}))\right\\}$
$\displaystyle(i\alpha((1+a)(D\tau_{0}^{22})-(1-a)(D\tau_{0}^{11}))+2(1+a^{2})\alpha^{2}\tau_{0}^{12})+4S_{0}(S_{3}+(1+a)S_{1})(S_{1}+S_{3})(DS_{3}$
$\displaystyle\left.\left.+(1+a)DS_{1})(i\alpha((1+a)D\tau_{0}^{22}-(1-a)D\tau_{0}^{11})+2a\alpha^{2}\tau_{0}^{12})\right]\right],$
(20a) $\displaystyle c_{1}$
$\displaystyle=\frac{-4(S_{3}+(1+a)S_{1})^{2}}{(4S_{0}^{2}d_{0}^{3}(S_{3}+(1+a)S_{1})^{2})}\left[\left[(S_{0})^{2}\left\\{2(D\tilde{A})d_{0}^{2}+(D^{2}\tilde{B})d_{0}^{2}-\tilde{B}d_{0}(D^{2}d_{0})-2(\tilde{A})d_{0}(Dd_{0})-2(D\tilde{B})\right.\right.\right.$
$\displaystyle\left.d_{0}(Dd_{0})+2(\tilde{B})(Dd_{0})+\alpha^{2}\tilde{B}d_{0}^{2}\right\\}+2(1-a^{2})(i\alpha)\left\\{S_{0}((D^{2}U)(S_{1}+S_{3})+(DU)(DS_{1}+DS_{3}))-\right.$
$\displaystyle\left.\left.(DU)(S_{1}+S_{3})(D(S_{0})\right\\}\tilde{B}d_{0}^{2}+2(1-a^{2})(i\alpha)(DU)(S_{1}+S_{3})S_{0}\left\\{(\tilde{A}+D\tilde{B})d_{0}^{2}-(\tilde{B})d_{0}(Dd_{0})\right\\}\right]+d_{0}^{3}$
$\displaystyle\left[2S_{0}(S_{3}+(1+a)S_{1})\left\\{(2(S_{1}+S_{3})(S_{2}-(1+a)S_{1})+1+2(S_{3}+(1+a)S_{1})(S_{1}+S_{3}))(i\alpha)(1+a)(1\right.\right.$
$\displaystyle-2a)D\tau_{0}^{22}+(2(S_{1}+S_{3})(S_{2}-(1+a)S_{1})+1-2(S_{1}+S_{3})(S_{3}+(1+a)S_{1}))i\alpha(1-a)(1+2a)D\tau_{0}^{11}$
$\displaystyle\left.+(2a(S_{1}+S_{3})(S_{2}-(1+a)S_{1})+(S_{1}+S_{3})(S_{3}+(1+a)S_{1})(1+a^{2}))4\alpha^{2}\tau_{0}^{12}\right\\}+2\left\\{(DS_{1}+DS_{3})\right.$
$\displaystyle(S_{2}-(1+a)S_{1})+(S_{1}+S_{3})(DS_{2}-(1+a)DS_{1})-((DS_{0})(S_{3}+(1+a)S_{1})+S_{0}(DS_{3}+(1+a)DS_{1})$
$\displaystyle\left.)\right\\}\\!\\!(2(S_{1}\\!\\!+\\!\\!S_{3})(S_{2}\\!\\!-\\!\\!(1+a)S_{1})\\!\\!+\\!\\!1)(2ai\alpha((1-a)\tau_{0}^{11}-(1+a)\tau_{0}^{22})-\frac{4ai\alpha\nu}{ReE})-4\left\\{((DS_{1}+DS_{3})(S_{3}\right.$
$\displaystyle+(1+a)S_{1})+(S_{1}+S_{3})(DS_{3}+(1+a)DS_{1}))-((DS_{0}(S_{3}+(1+a)S_{1})+S_{0}(DS_{3}+(1+a)DS_{1}))(S_{1}$
$\displaystyle\left.\left.+S_{3})(S_{3}+(1+a)S_{1})\right\\}(2ai\alpha((1-a)\tau_{0}^{11}+(1+a)\tau_{0}^{22})+\frac{4i\alpha\nu}{ReE})\right],$
(20b) $\displaystyle c_{2}$
$\displaystyle=\frac{-4(S_{3}+(1+a)S_{1})^{2}}{(4S_{0}^{2}d_{0}^{3}(S_{3}+(1+a)S_{1})^{2})}\left[\left[S_{0}^{2}\left\\{\tilde{A}d_{0}^{2}+2(D\tilde{B})d_{0}^{2}+(D^{2}\tilde{C})d_{0}^{2}-\tilde{C}d_{0}(D^{2}d_{0})-2\tilde{B}d_{0}(Dd_{0})-\right.\right.\right.$
$\displaystyle\left.2(D\tilde{C})d_{0}(Dd_{0})+2\tilde{C}(Dd_{0})^{2}+\alpha^{2}\tilde{C}d_{0}^{2}\right\\}+2(1-a^{2})(i\alpha)\left\\{S_{0}((D^{2}U)(S_{1}+S_{3})+DU(DS_{1}+\right.$
$\displaystyle\left.\left.DS_{3}))-DU(S_{1}+S_{3})DS_{0}\right\\}\tilde{C}d_{0}^{2}+2i\alpha(1-a^{2})DU(S_{1}+S_{3})S_{0}\left\\{(\tilde{B}+D\tilde{C})d_{0}^{2}-\tilde{C}d_{0}(Dd_{0})\right\\}\right]$
$\displaystyle+d_{0}^{3}\left[2S_{0}(S_{3}+(1+a)S_{1})\left\\{(2(S_{1}+S_{3})(S_{2}-(1+a)S_{1})+1-2(S_{3}+(1+a)S_{1})(S_{1}+S_{3})\right.\right.$
$\displaystyle
2ai\alpha(1-a)(\tau_{0}^{11})-(2(S_{1}+S_{3})(S_{2}-(1+a)S_{1})+1+2(S_{1}+S_{3})(S_{3}+(1+a)S_{1}))2a(i\alpha)$
$\displaystyle(1+a)\tau_{0}^{22}\\!-\\!\left\\{(2(S_{1}\\!\\!+\\!\\!S_{3})(S_{2}\\!\\!-\\!\\!(1\\!\\!+\\!\\!a)S_{1})+1)a+2(S_{1}+S_{3})(S_{3}+(1+a)S_{1}))\right\\}\frac{4i\alpha\nu}{ReE}\\!-\\!4(1\\!-\\!a^{2})$
$\displaystyle\left.(S_{1}+S_{3})(S_{3}+(1+a)S_{1})(D\tau_{0}^{12})\right\\}-8(S_{3}+(1+a)S_{1})^{2}(1-a^{2})\tau_{0}^{12}(S_{0}(DS_{1}+DS_{3})-(DS_{0})(S_{1}+$
$\displaystyle\left.\left.S_{3}))\right]\right],$ (20c) $\displaystyle c_{3}$
$\displaystyle=-\frac{1}{(S_{0}d_{0}^{2})}\left[S_{0}\left\\{(\tilde{B})d_{0}+2(D\tilde{C})d_{0}-2(\tilde{C})(Dd_{0})\right\\}+2(1-a^{2})(i\alpha)(DU)(S_{1}+S_{3})(\tilde{C})d_{0}+2d_{0}^{2}\right.$
$\displaystyle\left.(S_{1}+S_{3})(1-a^{2})(\tau_{0}^{12})\right],$ (20d)
$\displaystyle c_{4}$ $\displaystyle=-\frac{(\tilde{C})}{(d_{0})}.$ (20e)

where the following coefficients are used in the definition of $c_{i}$’s
above,

$\displaystyle\tilde{A}$
$\displaystyle=\left(\frac{(1-a)}{2}(DU)+\varepsilon\tau_{0}^{12}\right)(-\frac{A}{H})+(\varepsilon\tau_{0}^{12}-\frac{(1+a)}{2}(DU))(F_{1}+\frac{(A)(G)}{H})-(i\alpha)(D\tau_{0}^{12})+\frac{(1-a)}{2}$
$\displaystyle\alpha^{2}\tau_{0}^{22}-\frac{(1+a)}{2}\alpha^{2}\tau_{0}^{11}-\frac{\alpha^{2}\nu}{ReE},$
$\displaystyle\tilde{B}$
$\displaystyle=\left(\frac{(1-a)}{2}(DU)+\varepsilon\tau_{0}^{12}\right)(-\frac{B}{H})+(\varepsilon\tau_{0}^{12}-\frac{(1+a)}{2}(DU))(F_{2}+\frac{(B)(G)}{H}),$
$\displaystyle\tilde{C}$
$\displaystyle=\left(\frac{(1-a)}{2}(DU)+\varepsilon\tau_{0}^{12}\right)(\frac{C}{H})-(\varepsilon\tau_{0}^{12}-\frac{(1+a)}{2}(DU))(\frac{(C)(G)}{H})+\frac{(1-a)}{2}\tau_{0}^{11}-\frac{(1+a)}{2}\tau_{0}^{22}$
$\displaystyle-\frac{\nu}{ReE},$ $\displaystyle A$
$\displaystyle=\frac{2(1-a^{2})}{(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))}\left[(i\alpha(D\tau_{0}^{22})+(1+a)\alpha^{2}\tau_{0}^{12})\varepsilon\tau_{0}^{11}-\right.$
$\displaystyle\left.(i\alpha(D\tau_{0}^{11})-(1-a)\alpha^{2}\tau_{0}^{12})(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))\right],$
$\displaystyle B$
$\displaystyle=\frac{-4i\alpha(1-a^{2})}{(\varepsilon\tau_{0}^{11}(1\\!\\!-\\!\\!a)\\!\\!+\\!\\!(i\alpha(U\\!\\!-\\!\\!c)\\!\\!+\\!\\!\frac{\nu}{ReE}\\!\\!+\\!\\!\varepsilon(\tau_{0}^{11}\\!\\!+\\!\\!2\tau_{0}^{22}))(1\\!\\!+\\!\\!a))}\big{[}(\varepsilon\tau_{0}^{11})(a\tau_{0}^{22}+\frac{\nu}{ReE})+(a\tau_{0}^{11}+\frac{\nu}{ReE})$
$\displaystyle\left.(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))\right],$
$\displaystyle C$ $\displaystyle=2(1-a^{2})\tau_{0}^{12},$ $\displaystyle H$
$\displaystyle=\frac{2(1-a^{2})}{(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))}\\!\\!\left[(i\alpha(U\\!\\!-\\!\\!c)\\!\\!+\\!\\!\frac{\nu}{ReE})^{2}\\!\\!+\\!\\!6\varepsilon((i\alpha)\right.$
$\displaystyle\left.(U-c)+\frac{\nu}{ReE})(\tau_{0}^{11}+2\tau_{0}^{22})+2(\varepsilon)^{2}(\tau_{0}^{11}+\tau_{0}^{22})^{2}\right],$
$\displaystyle F_{1}$
$\displaystyle=\frac{(i\alpha((1+a)(D\tau_{0}^{22})+(1-a)(D\tau_{0}^{11}))+4a\alpha^{2}\tau_{0}^{12})}{(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))},$
$\displaystyle F_{2}$
$\displaystyle=\frac{(2ai\alpha((1-a)\tau_{0}^{11}-(1+a)\tau_{0}^{22}-\frac{2\nu}{ReE}))}{(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))},$
$\displaystyle G$
$\displaystyle=\frac{((i\alpha(U-c)+2\varepsilon\tau_{0}^{11}+\frac{\nu}{ReE})(1-a)+2\varepsilon\tau_{0}^{22})}{(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))},$
$\displaystyle d_{0}$
$\displaystyle=\frac{2(1-a^{2})(DU)}{H}\left[(\frac{1-a}{2})(DU)+\varepsilon\tau_{0}^{12}\right]\\!\\!+\\!\\!S_{1}\\!\\!+\\!\\!\frac{2(1-a^{2})(DU)(G)}{H}\\!\\!\left[(\frac{1+a}{2})(DU)\\!\\!-\\!\\!\varepsilon\tau_{0}^{12}\right],$
$\displaystyle S_{1}$
$\displaystyle=i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+\tau_{0}^{22}),$
$\displaystyle S_{2}$
$\displaystyle=\varepsilon(1-a)\tau_{0}^{11}-\varepsilon(1+a)\tau_{0}^{22},$
$\displaystyle S_{3}$
$\displaystyle=\varepsilon(1-a)\tau_{0}^{11}+\varepsilon(1+a)\tau_{0}^{22},$
$\displaystyle S_{0}$
$\displaystyle=aS_{1}S_{2}+S_{1}S_{3}+(1-a^{2})S_{1}^{2}.$ (21)

along with the first derivatives,

$\displaystyle D\tilde{A}$
$\displaystyle=((\frac{1-a}{2})(D^{2}U)+\varepsilon(D\tau_{0}^{12}))(-\frac{A}{H})+((\frac{1-a}{2})(DU)+\varepsilon\tau_{0}^{12})(\frac{(DH)A-H(DA)}{H^{2}})+(\varepsilon(D\tau_{0}^{12})$
$\displaystyle-(\frac{1+a}{2})(D^{2}U))(F_{1}+\frac{AG}{H})+(\varepsilon\tau_{0}^{12}-(\frac{1+a}{2})(DU))(DF_{1}+\frac{HA(DG)+HG(DA)-AG(DH)}{H^{2}})$
$\displaystyle-i\alpha(D^{2}\tau_{0}^{12})+(\frac{1-a}{2})\alpha^{2}(D\tau_{0}^{22})-(\frac{1+a}{2})\alpha^{2}(D\tau_{0}^{11}),$
$\displaystyle D\tilde{B}$
$\displaystyle=((\frac{1-a}{2})(D^{2}U)+\varepsilon(D\tau_{0}^{12}))(-\frac{B}{H})+((\frac{1-a}{2})(DU)+\varepsilon\tau_{0}^{12})(\frac{(DH)B-H(DB)}{H^{2}})+(\varepsilon(D\tau_{0}^{12})$
$\displaystyle-(\frac{1+a}{2})(D^{2}U))(F_{2}+\frac{BG}{H})+(\varepsilon\tau_{0}^{12}-(\frac{1+a}{2})(DU))(DF_{2}+\frac{HB(DG)+HG(DB)-BG(DH)}{H^{2}}),$
$\displaystyle D\tilde{C}$
$\displaystyle=((\frac{1-a}{2})(D^{2}U)+\varepsilon(D\tau_{0}^{12}))(\frac{C}{H})+((\frac{1-a}{2})(DU)+\varepsilon\tau_{0}^{12})(\frac{(DC)H-C(DH)}{H^{2}})-(\varepsilon(D\tau_{0}^{12})$
$\displaystyle-(\frac{1+a}{2})(D^{2}U))(\frac{CG}{H})-(\varepsilon\tau_{0}^{12}-(\frac{1+a}{2})(DU))(\frac{HC(DG)+HG(DC)-CG(DH)}{H^{2}})+(\frac{1-a}{2})$
$\displaystyle(D\tau_{0}^{11})-(\frac{1+a}{2})(D\tau_{0}^{22}),$
$\displaystyle DA$
$\displaystyle=\frac{2(1-a^{2})}{(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))^{2}}\left[(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\right.$
$\displaystyle\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))\left\\{(i\alpha
D^{2}\tau_{0}^{22}+(1+a)\alpha^{2}D\tau_{0}^{12})(\varepsilon\tau_{0}^{11})+(i\alpha
D\tau_{0}^{22}+(1+a)\alpha^{2}\tau_{0}^{12})\right.$
$\displaystyle(\varepsilon D\tau_{0}^{11})-(i\alpha
D^{2}\tau_{0}^{11}\\!\\!-\\!\\!(1-a)\alpha^{2}D\tau_{0}^{12})(i\alpha(U\\!\\!-\\!\\!c)\\!\\!+\\!\\!\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))-(i\alpha
D\tau_{0}^{11}-(1-a)$
$\displaystyle\left.\alpha^{2}\tau_{0}^{12})(i\alpha(DU)+\varepsilon(D\tau_{0}^{11}+2D\tau_{0}^{22})\right\\}-\left\\{(i\alpha
D\tau_{0}^{22}+(1+a)\alpha^{2}\tau_{0}^{12})(\varepsilon\tau_{0}^{11})-(i\alpha
D\tau_{0}^{11}-(1-a)\right.$
$\displaystyle\left.\left.\alpha^{2}\tau_{0}^{12})(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))\right\\}(2\varepsilon
D\tau_{0}^{11}+(1+a)(i\alpha(DU)+2\varepsilon D\tau_{0}^{22}))\right],$
$\displaystyle DB$
$\displaystyle=\frac{-4i\alpha(1-a^{2})}{(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))^{2}}\left[(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\right.$
$\displaystyle\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))\left\\{\varepsilon
D\tau_{0}^{11}(a\tau_{0}^{22}+\frac{\nu}{ReE})+a\varepsilon\tau_{0}^{11}D\tau_{0}^{22}+(aD\tau_{0}^{11})(i\alpha(U-c)+\frac{\nu}{ReE}\right.$
$\displaystyle\left.+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))+(a\tau_{0}^{11}+\frac{\nu}{ReE})(i\alpha(DU)+\varepsilon(D\tau_{0}^{11}+2D\tau_{0}^{22}))\right\\}-(2\varepsilon
D\tau_{0}^{11}+(i\alpha(DU)+2\varepsilon$
$\displaystyle\left.D\tau_{0}^{22})(1+a))\left\\{\varepsilon\tau_{0}^{11}(a\tau_{0}^{22}+\frac{\nu}{ReE})+(a\tau_{0}^{11}+\frac{\nu}{ReE})(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))\right\\}\right],$
$\displaystyle DC$ $\displaystyle=2(1-a^{2})(D\tau_{0}^{12}),$ $\displaystyle
DH$
$\displaystyle=\frac{2(1-a^{2})}{(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))^{2}}\left[(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\right.$
$\displaystyle\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))\left\\{2i\alpha(DU)(i\alpha(U-c)+\frac{\nu}{ReE})+6\varepsilon
i\alpha(DU)(\tau_{0}^{11}+\tau_{0}^{22})+6\varepsilon(i\alpha\right.$
$\displaystyle\left.(U-c)+\frac{\nu}{ReE})(D\tau_{0}^{11}+D\tau_{0}^{22})+(4\varepsilon^{2})(\tau_{0}^{11}+\tau_{0}^{22})(D\tau_{0}^{11}+D\tau_{0}^{22})\right\\}-(2\varepsilon(D\tau_{0}^{11})+(i\alpha(DU)+$
$\displaystyle\left.2\varepsilon(D\tau_{0}^{22})(1\\!+\\!a))\\!\\!\left\\{(i\alpha(U\\!\\!-\\!\\!c)+\frac{\nu}{ReE})^{2}\\!\\!+\\!\\!6\varepsilon(i\alpha(U\\!\\!-\\!\\!c)+\frac{\nu}{ReE})(\tau_{0}^{11}\\!+\\!\tau_{0}^{22})+2\varepsilon^{2}(\tau_{0}^{11}+\tau_{0}^{22})^{2}\right\\}\right],$
$\displaystyle DF_{1}$
$\displaystyle=\frac{1}{(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))^{2}}\left[(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\right.$
$\displaystyle\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))(i\alpha((1+a)D^{2}\tau_{0}^{22}+(1-a)D^{2}\tau_{0}^{11})+4a\alpha^{2}(D\tau_{0}^{12}))-(2\varepsilon(D\tau_{0}^{11})$
$\displaystyle\left.+(i\alpha(DU)+2\varepsilon(D\tau_{0}^{22})(1+a))(i\alpha((1+a)D\tau_{0}^{22}+(1-a)D\tau_{0}^{11})+4a\alpha^{2}\tau_{0}^{12})\right],$
$\displaystyle DF_{2}$
$\displaystyle=\frac{2ai\alpha}{(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))^{2}}\left[(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\right.$
$\displaystyle\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))((1-a)(D\tau_{0}^{11})-(1+a)(D\tau_{0}^{22}))-(2\varepsilon(D\tau_{0}^{11})+(i\alpha(DU)+2\varepsilon$
$\displaystyle\left.(D\tau_{0}^{22}))(1+a))((1-a)\tau_{0}^{11}-(1+a)\tau_{0}^{22}-\frac{2\nu}{ReE})\right],$
$\displaystyle DG$
$\displaystyle=\frac{1}{(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))^{2}}\left[(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\right.$
$\displaystyle\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))((i\alpha(DU)+2\varepsilon(D\tau_{0}^{11}))(1-a)+2\varepsilon(D\tau_{0}^{22}))-((i\alpha(DU)+2\varepsilon$
$\displaystyle\left.(D\tau_{0}^{22}))(1+a)+2\varepsilon(D\tau_{0}^{11}))((i\alpha(U-c)+\frac{\nu}{ReE}+2\varepsilon\tau_{0}^{11})(1-a)+2\varepsilon\tau_{0}^{22})\right],$
$\displaystyle Dd_{0}$
$\displaystyle=DS_{1}+2(1-a^{2})\left[((\frac{1-a}{2})(D^{2}U)+\varepsilon(D\tau_{0}^{12}))(\frac{(DU)}{H})+((\frac{1-a}{2})(DU)+\varepsilon(\tau_{0}^{12}))\right.$
$\displaystyle(\frac{(H(D^{2}U)-(DU)(DH))}{H^{2}})-((\varepsilon
D\tau_{0}^{12}-(\frac{1+a}{2})(D^{2}U))(\frac{(DU)G}{H})+(\varepsilon\tau_{0}^{12}-(\frac{1+a}{2})(DU))$
$\displaystyle\left.(\frac{H(D^{2}U)G+H(DU)(DG)-(DH)(DU)G}{H^{2}})\right],$
$\displaystyle DS_{1}$
$\displaystyle=i\alpha(DU)+\varepsilon(D\tau_{0}^{11}+D\tau_{0}^{22}),$
$\displaystyle DS_{2}$
$\displaystyle=\varepsilon(1-a)D\tau_{0}^{11}-\varepsilon(1+a)D\tau_{0}^{22},$
$\displaystyle DS_{3}$
$\displaystyle=\varepsilon(1-a)D\tau_{0}^{11}+\varepsilon(1+a)D\tau_{0}^{22},$
$\displaystyle DS_{0}$
$\displaystyle=a(DS_{1})S_{2}+aS_{1}(DS_{2})+(DS_{1})S_{3}+S_{1}(DS_{3})+2(1-a^{2})S_{1}(DS_{1}).$
(22)

and the second derivatives,

$\displaystyle D^{2}\tilde{A}$
$\displaystyle=((\frac{1-a}{2})(D^{3}U)+\varepsilon(D^{2}\tau_{0}^{12}))(-\frac{A}{H})+2((\frac{1-a}{2})(D^{2}U)+\varepsilon(D\tau_{0}^{12}))(\frac{(DH)A-H(DA)}{H^{2}})+((\frac{1-a}{2})$
$\displaystyle(DU)+\varepsilon\tau_{0}^{12})\frac{H((D^{2}H)A-H(D^{2}A))-2(DH)((DH)A-H(DA))}{H^{3}}+(\varepsilon
D^{2}\tau_{0}^{12}-(\frac{1+a}{2})(D^{3}U))(F_{1}+$
$\displaystyle\frac{AG}{H})+2(\varepsilon
D\tau_{0}^{12}-(\frac{1+a}{2})(D^{2}U))(DF_{1}+\frac{HA(DG)+HG(DA)-AG(DH)}{H^{2}})+(\varepsilon\tau_{0}^{12}-(\frac{1+a}{2})(DU))$
$\displaystyle(D^{2}F_{1}+\frac{1}{H^{3}}(H^{2}(A(D^{2}G)+2(DG)(DH)+G(D^{2}A))-AGH(D^{2}H))-2(DH)(H(A(DG)+G(DA))-$
$\displaystyle AG(DH)))-i\alpha
D^{3}\tau_{0}^{12}+(\frac{1-a}{2})\alpha^{2}D^{2}\tau_{0}^{22}-(\frac{1+a}{2})\alpha^{2}D^{2}\tau_{0}^{11},$

$\displaystyle D^{2}\tilde{B}$
$\displaystyle=((\frac{1-a}{2})(D^{3}U)+\varepsilon(D^{2}\tau_{0}^{12}))(-\frac{B}{H})+2((\frac{1-a}{2})(D^{2}U)+\varepsilon
D\tau_{0}^{12})(\frac{(DH)B-H(DB)}{H^{2}})$
$\displaystyle+((\frac{1-a}{2})(DU)+\varepsilon\tau_{0}^{12})(\frac{H(D^{2}H)B-H^{2}(D^{2}B)-2(DH)(B(DH)-H(DB))}{H^{3}})+(\varepsilon(D^{2}\tau_{0}^{12})$
$\displaystyle-(\frac{1+a}{2})(D^{3}U))(F_{2}+\frac{BG}{H})+2(\varepsilon
D\tau_{0}^{12}-(\frac{1+a}{2})(D^{2}U))(DF_{2}+\frac{HB(DG)+HG(DB)-BG(DH)}{H^{2}})$
$\displaystyle+(\varepsilon\tau_{0}^{12}-(\frac{1+a}{2})(DU))(D^{2}F_{2}+\frac{1}{H^{3}}(H^{2}B(D^{2}G)+2H^{2}(DG)(DB)+H^{2}G(D^{2}B)-BGH(D^{2}H)$
$\displaystyle-2(DH)(HB(DG)+HG(DB)-BG(DH)))),$ $\displaystyle D^{2}\tilde{C}$
$\displaystyle=((\frac{1-a}{2})(D^{3}U)+\varepsilon(D^{2}\tau_{0}^{12}))(\frac{C}{H})+2((\frac{1-a}{2})(D^{2}U)+\varepsilon(D\tau_{0}^{12}))(\frac{(DC)H-C(DH)}{H^{2}})+((\frac{1-a}{2})$
$\displaystyle(DU)+\varepsilon\tau_{0}^{12})\frac{H^{2}(D^{2}C)-CH(D^{2}H)-2(DH)((DC)H-C(DH))}{H^{3}}-(\varepsilon
D^{2}\tau_{0}^{12}-(\frac{1+a}{2})(D^{3}U))(\frac{CG}{H})$
$\displaystyle-2(\varepsilon
D\tau_{0}^{12}-(\frac{1+a}{2})D^{2}U)(\frac{HC(DG)+HG(DC)-CG(DH)}{H^{2}})-(\varepsilon\tau_{0}^{12}-DU(\frac{1+a}{2}))$
$\displaystyle(\frac{H^{2}(C(D^{2}G)+2(DG)(DC)+G(D^{2}C))-CGH(D^{2}H)-2(DH)(HC(DG)+HG(DC)-CG(DH))}{H^{3}})$
$\displaystyle+(\frac{1-a}{2})D^{2}\tau_{0}^{11}-(\frac{1+a}{2})D^{2}\tau_{0}^{22},$
$\displaystyle D^{2}G$
$\displaystyle=\frac{1}{(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))^{3}}\left[(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\right.$
$\displaystyle\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))\left\\{(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))\right.$
$\displaystyle((i\alpha(D^{2}U)+2\varepsilon
D^{2}\tau_{0}^{11})(1-a)+2\varepsilon D^{2}\tau_{0}^{22})-(2\varepsilon
D^{2}\tau_{0}^{11}+(i\alpha(D^{2}U)+2\varepsilon D^{2}\tau_{0}^{22})(1+a))$
$\displaystyle\left.((i\alpha(U-c)+\frac{\nu}{ReE}+2\varepsilon\tau_{0}^{11})(1-a)+2\varepsilon\tau_{0}^{22})\right\\}-2(2\varepsilon
D\tau_{0}^{11}+(i\alpha(DU)+2\varepsilon D\tau_{0}^{22})(1+a))$
$\displaystyle\left\\{(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))(2\varepsilon
D\tau_{0}^{22}+(i\alpha(DU)+2\varepsilon D\tau_{0}^{11})\right.$
$\displaystyle\left.\left.(1-a))-(2\varepsilon
D\tau_{0}^{11}+(i\alpha(DU)+2\varepsilon
D\tau_{0}^{22})(1+a))((i\alpha(U-c)+\frac{\nu}{ReE}+2\varepsilon\tau_{0}^{11})(1-a)+2\varepsilon\tau_{0}^{22})\right\\}\right],$
$\displaystyle D^{2}F_{1}$
$\displaystyle=\frac{1}{(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))^{3}}\left[(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\right.$
$\displaystyle\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))\left\\{(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))\right.$
$\displaystyle(i\alpha((1+a)D^{3}\tau_{0}^{22}+(1-a)D^{3}\tau_{0}^{11})+4a\alpha^{2}D^{2}\tau_{0}^{12})\\!-\\!(2\varepsilon
D^{2}\tau_{0}^{11}\\!+\\!(i\alpha(D^{2}U)+2\varepsilon
D^{2}\tau_{0}^{22})(1+a))$
$\displaystyle\left.(i\alpha((1+a)D\tau_{0}^{22}+(1-a)D\tau_{0}^{11})+4a\alpha^{2}D\tau_{0}^{12})\right\\}-2(2\varepsilon
D\tau_{0}^{11}+(i\alpha(DU)+2\varepsilon
D\tau_{0}^{22})(1+a))\left\\{(\varepsilon\right.$
$\displaystyle\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))(i\alpha((1+a)D^{2}\tau_{0}^{22}+(1-a)D^{2}\tau_{0}^{11})+4a$
$\displaystyle\left.\left.\alpha^{2}D\tau_{0}^{12})-(2\varepsilon
D\tau_{0}^{11}+(i\alpha(DU)+2\varepsilon
D\tau_{0}^{22})(1+a))(i\alpha((1+a)D\tau_{0}^{22}+(1-a)D\tau_{0}^{11})+4a\alpha^{2}\tau_{0}^{12})\right\\}\right],$
$\displaystyle D^{2}C$ $\displaystyle=2(1-a^{2})D^{2}\tau_{0}^{12},$
$\displaystyle D^{2}A$
$\displaystyle=\frac{2(1-a^{2})}{(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))^{3}}\left[(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}\right.$
$\displaystyle+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))\left\\{(2\varepsilon
D\tau_{0}^{11}+(i\alpha(DU)+2\varepsilon(D\tau_{0}^{22}))(1+a))((i\alpha
D^{2}\tau_{0}^{22}+(1+a)\alpha^{2}D\tau_{0}^{12})\right.$
$\displaystyle(\varepsilon\tau_{0}^{11})+(i\alpha
D\tau_{0}^{22}+(1+a)\alpha^{2}\tau_{0}^{12})(\varepsilon
D\tau_{0}^{11})-(i\alpha
D^{2}\tau_{0}^{11}-(1-a)\alpha^{2}D\tau_{0}^{12})(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon$
$\displaystyle(\tau_{0}^{11}+2\tau_{0}^{22}))-(i\alpha
D\tau_{0}^{11}-(1-a)\alpha^{2}\tau_{0}^{12})(i\alpha(DU)+\varepsilon(D\tau_{0}^{11}+2D\tau_{0}^{22})))+(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U$
$\displaystyle-c)\\!+\\!\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}\\!+\\!2\tau_{0}^{22}))(1+a))((i\alpha
D^{3}\tau_{0}^{22}+(1+a)\alpha^{2}D^{2}\tau_{0}^{12})(\varepsilon\tau_{0}^{11})+2(i\alpha
D^{2}\tau_{0}^{22}\\!+\\!(1\\!+\\!a)\alpha^{2}$ $\displaystyle
D\tau_{0}^{12})(\varepsilon D\tau_{0}^{11})+(i\alpha
D\tau_{0}^{22}+(1+a)\alpha^{2}\tau_{0}^{12})(\varepsilon
D^{2}\tau_{0}^{11})-(i\alpha
D^{3}\tau_{0}^{11}-(1-a)\alpha^{2}D^{2}\tau_{0}^{12})(i\alpha(U-c)+$
$\displaystyle\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))-2(i\alpha
D^{2}\tau_{0}^{11}-(1-a)\alpha^{2}D\tau_{0}^{12})(i\alpha(DU)+\varepsilon(D\tau_{0}^{11}+2D\tau_{0}^{22}))-(i\alpha
D\tau_{0}^{11}-$
$\displaystyle(1-a)\alpha^{2}\tau_{0}^{12})(i\alpha(D^{2}U)+\varepsilon(D^{2}\tau_{0}^{11}+2D^{2}\tau_{0}^{22}))-((i\alpha
D^{2}\tau_{0}^{22}+(1+a)\alpha^{2}D\tau_{0}^{12})(\varepsilon\tau_{0}^{11})+(i\alpha
D\tau_{0}^{22}$ $\displaystyle+(1+a)\alpha^{2}\tau_{0}^{12})(\varepsilon
D\tau_{0}^{11})-(i\alpha
D^{2}\tau_{0}^{11}-(1-a)\alpha^{2}D\tau_{0}^{12})(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))-(i\alpha$
$\displaystyle
D\tau_{0}^{11}-(1-a)\alpha^{2}\tau_{0}^{12})(i\alpha(DU)+\varepsilon(D\tau_{0}^{11}+2D\tau_{0}^{22})))(2\varepsilon
D\tau_{0}^{11}+(1+a)(i\alpha(DU)+2\varepsilon D\tau_{0}^{22}))-$
$\displaystyle((i\alpha
D\tau_{0}^{22}+(1+a)\alpha^{2}\tau_{0}^{12})(\varepsilon\tau_{0}^{11})-(i\alpha
D\tau_{0}^{11}-(1-a)\alpha^{2}\tau_{0}^{12})(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22})))$
$\displaystyle\left.(2\varepsilon
D^{2}\tau_{0}^{11}+(1+a)(i\alpha(D^{2}U)+2\varepsilon
D^{2}\tau_{0}^{22}))\right\\}-2(2\varepsilon
D\tau_{0}^{11}+(i\alpha(DU)+2\varepsilon
D\tau_{0}^{22})(1+a))\left\\{((i\alpha\right.$ $\displaystyle
D^{2}\tau_{0}^{22}+(1+a)\alpha^{2}D\tau_{0}^{12})(\varepsilon\tau_{0}^{11})+(i\alpha
D\tau_{0}^{22}+(1+a)\alpha^{2}\tau_{0}^{12})(\varepsilon
D\tau_{0}^{11})-(i\alpha D^{2}\tau_{0}^{11}-(1-a)\alpha^{2}D\tau_{0}^{12})$
$\displaystyle(i\alpha(U\\!\\!-\\!\\!c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))\\!-\\!(i\alpha
D\tau_{0}^{11}-(1-a)\alpha^{2}\tau_{0}^{12})(i\alpha(DU)+\varepsilon(D\tau_{0}^{11}+2D\tau_{0}^{22})))(\varepsilon\tau_{0}^{11}$
$\displaystyle(1-a)+(i\alpha(U\\!\\!-\\!\\!c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))\\!\\!-\\!\\!((i\alpha
D\tau_{0}^{22}+(1+a)\alpha^{2}\tau_{0}^{12})(\varepsilon\tau_{0}^{11})-(i\alpha
D\tau_{0}^{11}$
$\displaystyle\left.\left.-(1-a)\alpha^{2}\tau_{0}^{12})(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22})))(2\varepsilon
D\tau_{0}^{11}+(1+a)(i\alpha(DU)+2\varepsilon
D\tau_{0}^{22}))\right\\}\right],$ $\displaystyle D^{2}H$
$\displaystyle=\frac{2(1-a^{2})}{(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))^{3}}\left[(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}\right.$
$\displaystyle+\varepsilon(\tau_{0}^{11}\\!+\\!2\tau_{0}^{22}))(1\\!+\\!a))\left\\{(\varepsilon\tau_{0}^{11}(1\\!-\\!a)+(i\alpha(U\\!-\\!c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))(2(i\alpha(U-c)\right.$
$\displaystyle+\frac{\nu}{ReE})i\alpha(D^{2}U)\\!-\\!2\alpha^{2}(DU)^{2}\\!+\\!6\varepsilon
i\alpha(D^{2}U)(\tau_{0}^{11}\\!+\\!\tau_{0}^{22})\\!+\\!12\varepsilon
i\alpha(DU)(D\tau_{0}^{11}\\!+\\!D\tau_{0}^{22})+\\!6\\!\varepsilon(i\alpha(U\\!-\\!c)$
$\displaystyle+\frac{\nu}{ReE})(D^{2}\tau_{0}^{11}\\!\\!+\\!\\!D^{2}\tau_{0}^{22})\\!+\\!4\varepsilon^{2}(D\tau_{0}^{11}\\!\\!+\\!\\!D\tau_{0}^{22})^{2}\\!+\\!4\varepsilon^{2}(\tau_{0}^{11}\\!\\!+\\!\\!\tau_{0}^{22})\\!\\!\\!\left.(D^{2}\tau_{0}^{11}\\!\\!+\\!\\!D^{2}\tau_{0}^{22})\right\\}\\!-\\!(2\varepsilon(D^{2}\tau_{0}^{11}\\!+\\!D^{2}\tau_{0}^{22})$
$\displaystyle+i\alpha(D^{2}U)(1\\!\\!+\\!\\!a))((i\alpha(U\\!\\!-\\!\\!c)\\!\\!+\\!\\!\frac{\nu}{ReE})^{2}\\!\\!+\\!\\!6\varepsilon(i\alpha(U-c)+\frac{\nu}{ReE})(\tau_{0}^{11}+\tau_{0}^{22})+2\varepsilon^{2}(\tau_{0}^{11}+\tau_{0}^{22})^{2})-2(2\varepsilon$
$\displaystyle(D\tau_{0}^{11}\\!\\!+\\!\\!D\tau_{0}^{22})+i\alpha
DU)(1\\!\\!+\\!\\!a))\left\\{(\varepsilon\tau_{0}^{11}(1\\!\\!-\\!\\!a)\\!\\!+\\!\\!(i\alpha(U\\!\\!-\\!\\!c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))(2(i\alpha
DU\right.$ $\displaystyle(i\alpha(U-c)+\frac{\nu}{ReE})+6\varepsilon i\alpha
DU(\tau_{0}^{11}+\tau_{0}^{22})+6\varepsilon(i\alpha(U-c)+\frac{\nu}{ReE})(D\tau_{0}^{11}+D\tau_{0}^{22})+4\varepsilon^{2}(\tau_{0}^{11}+\tau_{0}^{22})$
$\displaystyle\left.(D\tau_{0}^{11}\\!\\!+\\!\\!D\tau_{0}^{22})\right\\}\\!\\!-\\!\\!(2\varepsilon
D\tau_{0}^{11}\\!\\!+\\!\\!(i\alpha(DU)+2\varepsilon
D\tau_{0}^{22})(1+a))((i\alpha(U\\!\\!-\\!\\!c)+\frac{\nu}{ReE})^{2}+6\varepsilon(i\alpha(U\\!\\!-\\!\\!c)+\frac{\nu}{ReE})$
$\displaystyle\left.\left.(\tau_{0}^{11}+\tau_{0}^{22})+2\varepsilon^{2}(\tau_{0}^{11}+\tau_{0}^{22})^{2})\right\\}\right],$

$\displaystyle D^{2}F_{2}$
$\displaystyle=\frac{2ai\alpha}{(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))^{3}}\left[(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}\right.$
$\displaystyle+\varepsilon(\tau_{0}^{11}\\!\\!+\\!\\!2\tau_{0}^{22}))(1+a))((\varepsilon\tau_{0}^{11}(1-a)\\!\\!+\\!\\!(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}\\!\\!+\\!\\!2\tau_{0}^{22}))(1+a))((1-a)(D^{2}\tau_{0}^{11})$
$\displaystyle-(1+a)(D^{2}\tau_{0}^{22}))\\!\\!-\\!\\!(2\varepsilon(D^{2}\tau_{0}^{11})+(i\alpha(D^{2}U)+2\varepsilon(D^{2}\tau_{0}^{22}))(1+a))((1-a)\tau_{0}^{11}\\!\\!-\\!\\!(1+a)\tau_{0}^{22}\\!-\\!\frac{2\nu}{ReE}))$
$\displaystyle-2(2\varepsilon D\tau_{0}^{11}\\!\\!+\\!\\!(i\alpha
DU+2\varepsilon
D\tau_{0}^{22})(1\\!\\!+\\!\\!a))((\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))$
$\displaystyle\left.((1-a)D\tau_{0}^{11}\\!\\!-\\!\\!(1\\!\\!+\\!\\!a)D\tau_{0}^{22})\\!\\!-\\!\\!(2\varepsilon
D\tau_{0}^{11}+(i\alpha DU\\!\\!+\\!\\!2\varepsilon
D\tau_{0}^{22})(1+a))((1-a)\tau_{0}^{11}\\!\\!-\\!\\!(1+a)\tau_{0}^{22}-\frac{2\nu}{ReE}))\right],$
$\displaystyle D^{2}d_{0}$
$\displaystyle=D^{2}S_{1}+2(1-a^{2})\left[((\frac{1-a}{2})(D^{3}U)+\varepsilon(D^{2}\tau_{0}^{12}))(\frac{(DU)}{H})+2((\frac{1-a}{2})(D^{2}U)+\varepsilon(D\tau_{0}^{12}))\right.$
$\displaystyle(\frac{(H(D^{2}U)-(DU)(DH))}{H^{2}})+(\frac{H(H(D^{3}U)-(DU)(D^{2}H))-2(DH)(H(D^{2}U)-(DH)(DU))}{H^{3}})$
$\displaystyle((\frac{1-a}{2})(DU)+\varepsilon\tau_{0}^{12})-\left\\{((\varepsilon
D^{2}\tau_{0}^{12}-(\frac{1+a}{2})(D^{3}U))(\frac{(DU)G}{H})+2(\varepsilon
D\tau_{0}^{12}-(\frac{1+a}{2})(D^{2}U))\right.$
$\displaystyle(\frac{H(D^{2}U)G+H(DU)(DG)-(DH)(DU)G}{H^{2}})+(\varepsilon\tau_{0}^{12}-(\frac{1+a}{2})(DU))\frac{1}{H^{3}}(H(H(D^{3}U)G+2H(D^{2}U)$
$\displaystyle\left.\left.(DG)+H(DU(D^{2}G)-G(DU)(D^{2}H)-2(DH)(H(D^{2}U)G+H(DU)(DG)-(DU)(DH)G)))\right\\}\right],$
$\displaystyle D^{2}B$
$\displaystyle=\frac{-4i\alpha(1-a^{2})}{(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))^{3}}\left[(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\right.$
$\displaystyle\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))(1+a))((2\varepsilon
D\tau_{0}^{11}+(i\alpha(DU)+2D\tau_{0}^{22})(1+a))(\varepsilon
D\tau_{0}^{11}(a\tau_{0}^{22}+\frac{\nu}{ReE})+$ $\displaystyle
a\varepsilon\tau_{0}^{11}D\tau_{0}^{22}+aD\tau_{0}^{11}(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))+(a\tau_{0}^{11}+\frac{\nu}{ReE})(i\alpha(DU)+\varepsilon(D\tau_{0}^{11}$
$\displaystyle+2D\tau_{0}^{22}))+(\varepsilon\tau_{0}^{11}(1-a)+(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}\\!\\!+\\!\\!2\tau_{0}^{22}))(1\\!\\!+\\!\\!a))(\varepsilon(D^{2}\tau_{0}^{11})(a\tau_{0}^{22}\\!\\!+\\!\\!\frac{\nu}{ReE})$
$\displaystyle+2a\varepsilon(D\tau_{0}^{11})(D\tau_{0}^{22})+a\varepsilon\tau_{0}^{11}(D^{2}\tau_{0}^{22})+(aD^{2}\tau_{0}^{11})(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))+2aD\tau_{0}^{11}$
$\displaystyle(i\alpha(DU)+\varepsilon(D\tau_{0}^{11}+2D\tau_{0}^{22}))+(a\tau_{0}^{11}+\frac{\nu}{ReE})(i\alpha(D^{2}U)+\varepsilon(D^{2}\tau_{0}^{11}+2D^{2}\tau_{0}^{22}))-(2\varepsilon
D^{2}\tau_{0}^{11}+$ $\displaystyle(i\alpha(D^{2}U)+2\varepsilon
D^{2}\tau_{0}^{22})(1+a))(\varepsilon\tau_{0}^{11})(a\tau_{0}^{22}+\frac{\nu}{ReE})+(a\tau_{0}^{11}+\frac{\nu}{ReE})(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}$
$\displaystyle+2\tau_{0}^{22})))-(2\varepsilon
D\tau_{0}^{11}+(i\alpha(DU)+2\varepsilon D\tau_{0}^{22})(1+a))(\varepsilon
D\tau_{0}^{11})(a\tau_{0}^{22}+\frac{\nu}{ReE})+a\varepsilon\tau_{0}^{11}D\tau_{0}^{22}+$
$\displaystyle
aD\tau_{0}^{11}(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))+(a\tau_{0}^{11}+\frac{\nu}{ReE})(i\alpha(DU)+\varepsilon(D\tau_{0}^{11}+2D\tau_{0}^{22}))))-$
$\displaystyle 2(2\varepsilon D\tau_{0}^{11}+(i\alpha(DU)+2\varepsilon
D\tau_{0}^{22})(1+a))(\varepsilon
D\tau_{0}^{11}(a\tau_{0}^{22}+\frac{\nu}{ReE})+a\varepsilon\tau_{0}^{11}D\tau_{0}^{22}+aD\tau_{0}^{11}(i\alpha(U-c)$
$\displaystyle+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))+(a\tau_{0}^{11}+\frac{\nu}{ReE})(i\alpha(DU)+\varepsilon(D\tau_{0}^{11}+2D\tau_{0}^{22}))-(2\varepsilon
D\tau_{0}^{11}+(i\alpha DU$ $\displaystyle\left.+2\varepsilon
D\tau_{0}^{22})(1+a))(\varepsilon\tau_{0}^{11}(a\tau_{0}^{22}+\frac{\nu}{ReE})+(a\tau_{0}^{11}+\frac{\nu}{ReE})(i\alpha(U-c)+\frac{\nu}{ReE}+\varepsilon(\tau_{0}^{11}+2\tau_{0}^{22}))))\right].$
(23)

The components of the base state extra elastic stress tensor, $\tau_{0}$, is
given by

$\displaystyle\tau_{0}^{11}$
$\displaystyle=\frac{1}{24\varepsilon^{2}a^{2}(\frac{ReE}{\nu})^{4}M_{1}^{\frac{1}{3}}}\\!\\!\\!\left[\\!\\!2(2)^{\frac{1}{3}}\\!\\!\varepsilon^{2}(\frac{ReE}{\nu})^{6}(1\\!\\!+\\!\\!a)^{4}\\!\\!\\!-\\!\\!4(2)^{\frac{1}{3}}\varepsilon^{2}(\frac{ReE}{\nu})^{6}(1\\!\\!+\\!\\!a)^{2}(1\\!\\!-\\!\\!a^{2})\\!\\!+\\!\\!2(2)^{\frac{1}{3}}\varepsilon^{2}(\frac{ReE}{\nu})^{6}(1\\!\\!-\\!\\!a^{2})^{2}\right.$
$\displaystyle-\\!\\!6(2)^{\frac{1}{3}}\varepsilon^{2}(\frac{ReE}{\nu})^{8}(1\\!\\!+\\!\\!a)^{4}(1\\!\\!-\\!\\!a^{2})U^{2}\\!\\!+\\!\\!12(2)^{\frac{1}{3}}\varepsilon^{2}(\frac{ReE}{\nu})^{8}(1\\!\\!+\\!\\!a)^{2}(1\\!\\!-\\!\\!a^{2})^{2}U^{2}\\!\\!-\\!\\!6(2)^{\frac{1}{3}}\varepsilon^{2}(\frac{ReE}{\nu})^{8}(1\\!\\!-\\!\\!a^{2})^{3}U^{2}$
$\displaystyle\left.-8a\varepsilon(1+a)(\frac{ReE}{\nu})^{3}M_{1}^{\frac{1}{3}}+(2M_{1})^{\frac{2}{3}}\right],$
$\displaystyle\tau_{0}^{22}$
$\displaystyle=-\frac{(1-a)}{(1+a)}\tau_{0}^{11},$
$\displaystyle\tau_{0}^{12}$
$\displaystyle=\frac{DU(1-(1-a)\tau_{0}^{11}(\frac{ReE}{\nu}))}{(1+\varepsilon(\frac{ReE}{\nu})(\tau_{0}^{11}+\tau_{0}^{22}))},$
(24)

where the coefficient $M_{1}$ in the base stress, is given by

$\displaystyle M_{1}$
$\displaystyle=16\varepsilon^{3}(\frac{ReE}{\nu})^{9}a^{3}(1\\!\\!+\\!\\!a)^{3}\\!\\!\\!\left[\\!\\!(27a\varepsilon+9(1\\!\\!-\\!\\!a^{2}))(\frac{ReE}{\nu})^{2}U^{2}\\!\\!+\\!\\!1\\!\\!+\\!\\!3(3)^{\frac{1}{2}}(\frac{ReE}{\nu})(27\varepsilon^{2}(\frac{ReE}{\nu})^{2}a^{2}U^{4}\\!\\!+\\!\\!(1\\!\\!-\\!\\!a^{2})\right.$
$\displaystyle\left.(U+(\frac{ReE}{\nu})^{2}(1-a^{2})U^{3})^{2}+2a\varepsilon(U^{2}+9(\frac{ReE}{\nu})^{2}(1-a^{2})U^{4}))^{\frac{1}{2}}\right].$
(25)

## References

  * [1] T. M. Squires and S. R. Quake.  Microfluidics: fluid physics at nanoliter scale.  Rev. Mod. Phys., 77:977–1026, 2005.
  * [2] C. M. White and M. G. Mungal.  Mechanics and prediction of turbulent drag reduction with polymer additives.  Ann. Rev. Fluid Mech., 40:235–256, 2008.
  * [3] S. Sircar and A. Roberts.  Ion mediated crosslink driven mucous swelling kinetics.  DCDS-B, 21(6):1937–1951, 2016.
  * [4] S. Sircar, E. Aisenbrey, S. J. Bryant, and D. M. Bortz.  Determining equilibrium osmolarity in poly(ethylene glycol)/chondrotin sulfate gels mimicking articular cartilage.  J. Theor. Biol., 364:397–406, 2015.
  * [5] S. Sircar and D. M. Bortz.  Impact of flow on ligand-mediated bacterial flocculation.  Math. Biosci., 245(2):314–321, 2013.
  * [6] S. Sircar, J. G. Younger, and D. M. Bortz.  Sticky surface: sphere-sphere adhesion dynamics, 2014.
  * [7] S. Sircar and A. J. Roberts.  Surface deformation and shear flow in ligand mediated cell adhesion.  J. Math. Biol., 73(4):1035–1052, 2016.
  * [8] S. Sircar, G. Nguyen, A. Kotousov, and A. J. Roberts.  Ligand-mediated adhesive mechanics of two static, deformed spheres.  Euro. Phys. J. E, 39(10):1–9, 2016.
  * [9] P. K. Ray and T. A. Zaki.  Absolute instability in viscoelastic mixing layers.  Phys. Fluids, 26:014103, 2014.
  * [10] A. Groisman and V. Steinberg.  Mechaism of elastic instability in couette flow of polymer solutions: experiments.  Phys. Fluids, 10(2451), 1998.
  * [11] Z. Yu and N. Phan-Thien.  Three-dimensional roll-up of a viscoelastic mixing layer.  J. Fluid Mech., 500:29–53, 2004.
  * [12] P. Huerre and P. A. Monkewitz.  Absolute and convective instabilities in free shear layers.  J. FLuid Mech., 159:151–168, 1985.
  * [13] R. Patne and V. Shankar.  Absolute and convective instabilities in combined couette-poiseuille flow past a neo-hookean solid.  Phys. Fluids, 29:124104, 2017.
  * [14] O. Tammisola M. P. Juniper and F. Lundell.  The local and global stability of confined planar wakes at intermediate reynolds number.  J. Fluid Mech., 686:218–238, 2011.
  * [15] T. Shoji, E. W. Harris, A. Bernard, S. G. Schein, and A. R. Karagozian.  On the origins of transverse jet shear layer instability transition.  J. Fluid Mech., 890:A7–1, 2020.
  * [16] Jr. H. Oertel.  Wakes behind blunt bodies.  Ann. Rev. Fluid Mech., 22:539–564, 1990.
  * [17] I. Delbende and J. M. Chomaz.  Nonlinear convective/absolute instabilities in parallel two-dimensional wakes.  Phys. Fluids, 10:2724, 1998.
  * [18] J. J. Healy.  Destabilizing effects of confinement on homogeneous mixing layers.  J. Fluid Mech., 623:241–271, 2009.
  * [19] B. Pier.  Local and global instabilities in the wake of a sphere.  J. Fluid Mech., 603:39–61, 2008.
  * [20] I. Vihinen, A. M. Honohan, and S. P. Lin.  Image of absolute instability in a liquid jet.  Phys. Fluids, 9:3117, 1997.
  * [21] C. Pipe.  Experiments investigating the effects of fluid elasticity on laminar vortex shedding from a cylinder.  Master’s thesis, EPFL, 2005.
  * [22] P. K. Ray and T. A. Zaki.  Absolute/convective instability of planar viscoelastic jets.  Phys. Fluids, 27:014110, 2015.
  * [23] A. Alhushaybar and J. Uddin.  Convective and absolute instability of viscoelastic liquid jets in the presence of gravity.  Phys. Fluids, 31:044106, 2019.
  * [24] A. Alhushaybari and J. Uddin.  Absolute instability of free-falling viscoelastic liquid jets with surfactants.  Phys. Fluids, 32:013102, 2020.
  * [25] F. T. Pinho M. C. Guimaraes, N. Pimentel and C. B. da Silva.  Direct numerical simulations of turbulent viscoelastic jets.  J. Fluid Mech., 899:A11, 2020.
  * [26] K. Kupfer, A. Bers, and A. K. Ram.  The cusp map in the complex-frequency plane for absolute instability.  Phys. Fluids, 30(10):3075–3082, 1987.
  * [27] R. B. Bird, R. C. Armstrong, and O. Hassager.  Dynamics of polymeric liquids, Volume 1: Fluid Mechanics.  Wiley Intersciences, 2nd edition, 1987.
  * [28] V. Steinberg.  Role of thermal noise in dynamics of non-equilibrium systems: Macro-, meso- and microscopic.  J. Stat. Phys., 175((3-4)):664–680, 2018.
  * [29] D. S. Malkus, J. A. Nohel, and B. J. Plohr.  Analysis of new phenomena in shear flow of non-newtonian fluids.  SIAM J. Appl. Math., 51(4):899–929, 1991.
  * [30] J. D. Evans.  Re-entrant corner behavior of the ptt fluid with a solvent viscosity.  J. Non Newt. Fluid Mech., 165(9):527–537, 2010.
  * [31] J. Azaiez and G. M. Homsy.  Linear stability of free shear flow of viscoelastic liquids.  J. Fluid. Mech., 268:37–69, 1994.
  * [32] S. Sircar and D. Bansal.  Spatiotemporal linear stability of viscoelastic free shear flows: Dilute regime.  Phys. Fluids, 31:084104, 2019.
  * [33] T. K. Sengupta.  Instabilities of Flows and Transition to Turbulence.  CRC Press Boca Raton, 2012.
  * [34] B.S. Ng and W. H. Reid.  The compound matrix method for ordinary differential equations.  J. Comp. Phys., 58:209–228, 1985.
  * [35] T. Hagen and M. Renardy.  Boundary layer analysis of the phan-thien-tanner and giesekus model in high weissenberg number flow.  J. Non-Newt. Fluid Mech., 73:181–189, 1997.
  * [36] L. L. Ferras, M. L. Morgado, M. Rebelo, G. H. McKinley, and A. M. Afonso.  A generalised phan-thien-tanner model.  J. Non Newt. Fluid Mech., 269:88–89, 2019.
  * [37] S. Taneda.  Experimental investigation of the wakes behind cylinders and plates at low reynolds numbers.  J. Phys. Soc. Japan, 11(3):302–307, 1956.
  * [38] P. A. Monkewitz.  The absolute and convective nature of instability in two-dimensional wakes at low reynolds numbers.  Phys. Fluids, 31(5):999–106, 1988.
  * [39] P. Huerre J. M. Chomaz and L. G. Redekopp.  Bifurcations to local and global modes in spatially developing flows.  Phys. Rev. Lett., 60(1):25–28, 1988.
  * [40] M. Provansal, C. Mathis, and L. Boyer.  Bénard-von kármán instability: transient and forced regimes.  J. Fluid Mech., 182:1–22, 1987.
  * [41] V. N. Kalashnikov and A. M. Kudin.  Kármán ortices in the flow of drag reducing polymer solutions.  Nature, 255:445–446, 1970.
  * [42] O. Cadot and S. Kumar.  Experimental characterization of viscoelastic effects on two- and three-dimensional shear instabilities.  J. Fluid Mech., 416:151–172, 2000.
  * [43] G. H. McKinley, R. C. Armstrong, and R. A. Brown.  The wake instability in viscoelastic flow past confined circular cylinders.  Phil. Trans. R. Soc. A, 344:265–304, 1993.
  * [44] P. M. Coelho and F. T. Pinho.  Vortex shedding in cylinder flow of shear-thinning fluids i. identification and demarcation of flow regimes.  J. Non Newt. Fluid Mech., 110:143–176, 2003.
  * [45] I. Chaudhary, P. Garg, G. Subramanian, and V. Shankar.  Linear instability of viscoelastic pipe flow.  J. Fluid Mech., 908(A11), 2021.
  * [46] D. Samanta, Y. Dubief, M. Holzner, C. Schafer, C. Wagner A. Morozov, and B. Hof.  Elasto-inertial turbulence.  Proc. Natl. Acad. Sci. USA, 110:10557, 2013.
  * [47] A. Bistagnino, G. Boffetta, A. Celani, A. Mazzino, A. Puliafito, and M. Vergassola.  Nonlinear dynamics of the viscoelastic kolmogorov flow.  J. Fluid Mech., 590:61–80, 2007.