# Nonlinear spatiotemporal instabilities in two-dimensional electroconvective
flows

Zhe Feng  Department of Mechanical Engineering, National University of
Singapore, 9 Engineering Drive 1, 117575 Singapore  Dongdong Wan  Department
of Mechanical Engineering, National University of Singapore, 9 Engineering
Drive 1, 117575 Singapore  Bo-Fu Wang  Shanghai Key Laboratory of Mechanics in
Energy Engineering, Shanghai Institute of Applied Mathematics and Mechanics,
School of Mechanics and Engineering Science, Shanghai University, Shanghai
200072, China  Mengqi Zhang  Department of Mechanical Engineering, National
University of Singapore, 9 Engineering Drive 1, 117575 Singapore

###### Abstract

This work studies the effects of a through-flow on two-dimensional
electrohydrodynamic (EHD) flows of a dielectric liquid confined between two
plane plates, as a model problem to further our understanding of the fluid
mechanics in the presence of an electric field. The liquid is subjected to a
strong unipolar charge injection from the bottom plate and a pressure gradient
along the streamwise direction (forming a Poiseuille flow). Highly-accurate
numerical simulations and weakly nonlinear stability analyses based on
multiple-scale expansion and amplitude expansion methods are used to unravel
the nonlinear spatiotemporal instability mechanisms in this combined flow. We
found that the through-flow makes the hysteresis loop in the EHD flow
narrower. In the numerical simulation of an impulse response, the leading and
trailing edges of the wavepacket within the nonlinear regime are consistent
with the linear ones, a result which we also verified against that in natural
convection. In addition, as the bifurcation in EHD-Poiseuille flows is of a
subcritical nature, nonlinear finite-amplitude solutions exist in the
subcritical regime, and our calculation indicates that they are convectively
unstable (at least for the parameters investigated). The validity of the
Ginzburg-Landau equation (GLE), derived from the weakly nonlinear expansion of
Navier-Stokes equations and the Maxwell’s equations in the quasi-electrostatic
limit, serving as a physical reduced-order model for probing the
spatiotemporal dynamics in this flow, has also been investigated. We found
that the coefficients in the GLE calculated using amplitude expansion method
can predict the absolute growth rates even when the parameters are away from
the linear critical conditions, compared favourably with the local dispersion
relation, whereas the validity range of the GLE derived from the multiple-
scale expansion method is confined to the vicinity of the linear critical
conditions.

##  I Introduction

In the current work, the electroconvection (EC) in a dielectric liquid with a
low electrical conductivity is studied. As the induced electric current is
weak in the dielectric liquid, the magnetic field can be neglected. Ions can
be generated via different methods (to be discussed below) and injected into
the liquid. Under a DC (i.e. direct current) electric field, the injected ions
can give rise to Coulomb force, which can influence and control the motion of
liquids in various industrial applications. The electrohydrodynamic (EHD)
conduction pumping technique can be used to electrically drive the dielectric
liquid flow in a single-phase loop (Feng and Seyed-Yagoobi, 2004). The charged
ions are generated in a dissociation process of the neutral electrolytic
species in the vicinity of electrodes. Aside from controlling the mass
transfer in dielectric liquids, the combined effects of an electrical field
and buoyancy can enhance the efficiency of heat transfer (McCluskey  _et al._
, 1991). For two miscible dielectric liquids under an electric field, the
Coulomb force can destabilise the interface by inducing fluid circulations
around the interface (Ward  _et al._ , 2019). The experimental results of Ref.
(Jalaal  _et al._ , 2013) indicate that such EHD mixers can significantly
increase the efficiency of mixing in liquid-liquid systems.

When the electroconvective flow is subjected to a through-flow, it can also
find numerous applications in various industrial processes. For example, Atten
& Honda Atten and Honda (1982) have investigated the EHD-Poiseuille flow to
understand the electroviscous (EV) phenomenon (referring to the increase of
apparent viscosity in a bulk flow as a result of the electrostatic effect).
The EV effect can be utilised to damp undesired vibration of devices in some
particular circumstances and such vibration dampers are useful in the case
where conventional mechanical dampers are inapplicable. Another relevant
application is the electrostatic precipitator (ESP), where a through-flow with
dusts and particles interacts with the electric field confined between two
plates with high-voltage wires in-between. ESPs serve as an effective and
reliable particulate control device for cleaning industrial and polluted gases
and can reduce the unwanted particles of various sizes Soldati  _et al._
(1993); Parker (1997); Mizuno (2000); Soldati (2000, 2003). Some ESP can
efficiently remove particles larger than
$2.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$, but the collection efficiency
decreases for smaller particles ($<$2.5\text{\,}\mathrm{\SIUnitSymbolMicro
m}$$) Jaworek  _et al._ (2018). Currently, we cannot find an effective
solution to this problem probably due to a lack of a systematic understanding
of how the electrified flow interacts with the through-flow.

In this work, we will study the problem of how a through-flow interacts with
the flow induced by an electric field from the perspective of flow
instability. Such knowledge can potentially in the end help us to, for
example, enhance the efficiency of ESP, design more efficient vibration damper
etc., as mentioned above. In order to analyse the physical problem with only
the most essential elements, in this work we simplify the problem by
considering the spatiotemporal instability of a Poiseuille flow subjected to
an electric field in a dielectric liquid confined between two infinite plane
walls. The ions are generated via the physicochemical process near the ion-
exchange membrane covering the walls. In subsequent work where a more similar
flow configuration to ESP is considered, the results may be applicable to the
flows in ESP more straightforwardly; we study the aforementioned academic
idealised flow in this work. In following sections, the linear and nonlinear
instabilities of this EHD flow with and without the through-flow, as well as
the linearly and nonlinearly convective and absolute instabilities of the EHD
flow and other flows, are reviewed.

###  I.1 Linear and weakly nonlinear instabilities of EHD flows

The EHD flow studied here is governed by the Navier-Stokes equations for the
flow field and the reduced Maxwell’s equation in the quasi-electrostatic limit
for the electric field Castellanos (1998). In the case of space-charge-limit
(SCL) injection (referring to an infinitely strong level of charge injection),
the linear stability criterion only depends on the electric Rayleigh number
$T$ Schneider and Watson (1970); Atten and Moreau (1972), which is a
dimensionless number quantifying the strength of external voltage added on the
top and bottom plates (see the problem formulation section II to follow). For
the EHD flow without a through-flow, Atten & Moreau Atten and Moreau (1972)
determined the linear stability criterion $T_{c}\approx 160.75$ (which has
also been subsequently confirmed in similar linear stability analyses (Wu  _et
al._ , 2015; Zhang  _et al._ , 2015; Wang and Sheu, 2016)). However, in the
EHD experiments using ion-exchange membranes to generate ions (Lacroix  _et
al._ , 1975; Atten and Lacroix, 1978), the onset of electroconvection took
place at about $T_{c}\approx 100\pm 10$. Continuous efforts have been devoted
to solving this now long-standing problem in the EC flows from several
perspectives, including the effect of neglected charge diffusion, the non-
modal energy growth and the validity of homogeneous charge injection
assumption. Pérez & Castellanos Pérez and Castellanos (1989) considered the
non-negligible effect of the charge diffusion on the linear stability
criterion and found that the value of $T_{c}$ decreases with a stronger charge
diffusion effect. Zhang et al. Zhang  _et al._ (2015) analysed the modal and
non-modal behaviours of this linearised EHD flow, but found that the transient
growth therein is small. Recently, Feng et al. Feng  _et al._ (2021) have
considered the possible deleterious effect of the inhomogeneity on ion
exchange membranes, which is inevitable in experiments but is conventionally
neglected in the stability analyses of EHD flows, by introducing some
stochasticity in the boundary conditions of their numerical simulations. Their
results showed that with stronger stochasticity, the electroconvection can
indeed take place at a lower electric Rayleigh number (rendering the numerical
prediction of $T_{c}$ closer to the experimental value).

In the case of EHD-Poiseuille flows, the combined effects of the through-flow
and the secondary flow induced by electric field have also been investigated.
This flow was experimentally studied by Atten & Honda Atten and Honda (1982).
At the onset of the EHD secondary motions, the apparent viscosity was found to
increase which can be indicated by the increased pressure gradient when
keeping a constant flow rate. The linear stability analysis of EHD-Poiseuille
flows was performed by Castellanos & Agrait Castellanos and Agrait (1992). At
low Reynolds numbers ($Re$), the most unstable modes are always transverse
rolls, and the transverse perturbations are stabilised. But they are
destabilised at high $Re$. As the decrease of ratio between hydrodynamic and
ionic mobilities, transverse rolls are also destabilised in the case of low
Reynolds numbers Lara  _et al._ (1997). The modal and non-modal stabilities of
EHD-Poiseuille flow were studied by Zhang et al. Zhang  _et al._ (2015). They
found that the linear instability criterion decreases compared with the no-
through-flow case, and the transient growth of energy increases with the
intensification of electric fields, which indicates that Poiseuille flows can
be destabilised by the electric field in a short time.

Beyond the linear instability $T_{c}$ in EHD flows, an infinitesimal
disturbance will continue to grow and the nonlinearity will start to manifest
itself when the disturbance amplitude is large enough. At this point, the
nonlinear effects should be analysed. A weakly nonlinear stability analysis is
an elegant mathematical tool that can unravel the effect of weak nonlinearity
in an analytical manner. This analysis approach was first conceptually
proposed by Landau Landau (1944), who, based on an argument of physical
symmetry, proposed the amplitude equation governing the evolution of
disturbance around the primary bifurcation at the linear criticality. When the
diffusion effect is considered in the derivation of the amplitude equation,
one can obtain the Ginzburg-Landau equation (GLE). In the context of
hydrodynamic stability, the GLE was first derived rigorously by Stuart using a
multiple-scale expansion method Stuart (1960) and by Watson using an amplitude
expansion method Watson (1960). These two methods have been proved to be
equivalent to each other at the critical point later Fujimura (1989). The
weakly nonlinear analysis has been applied successfully in canonical shear
flows and Rayleigh-Bénard convection Cross and Hohenberg (1993), etc. In EHD
flows, similar efforts have also been made to understand their nonlinear
effects using the weakly nonlinear stability theory, as reviewed in the
following.

In a study of two-dimensional EHD flow using a hydraulic model, Felici Felici
(1971) found two equilibrium solutions (one is of zero amplitude and the other
finite amplitude) at $T<T_{c}$, directly proving the existence of a
subcritical bifurcation in the EHD flow. The results of subcritical
bifurcation have been theoretically corroborated by Worraker & Richardson
Worraker and Richardson (1981) in their weakly nonlinear analysis of thermally
stabilized EHD flows (the flow becomes isothermal EHD flow when the thermal
difference between two plates is zero). A more systematic weakly nonlinear
analysis for the EHD flow has been conducted recently by Zhang Zhang (2016),
who investigated the primary bifurcation of EHD with and without through-flows
using the multiple-scale expansion method. His results show the universality
of the subcritical bifurcation in a large parameter space. This analysis is
confined in the vicinity of the linear critical conditions. The weakly
nonlinear behaviour of the EHD flow away from the linear criticality also
deserves to be investigated, which can be performed using the amplitude
expansion method (in which the linear growth rate is not necessarily small).
It is noted that the amplitude expansion method has recently been improved by
Suslov and his coworker Pham and Suslov (2018) (see also Ref. Cudby and
Lefauve (2021)) to be more rigorous and versatile and this method will be
utilised in the current work to probe the EHD flow when it is away from the
linear criticality.

###  I.2 Linear and nonlinear convective/absolute instabilities

When an infinitesimal impulse is added to a through-flow, the resultant
disturbance (taking the form of a wavepacket, that is, disturbance with a
broad spectrum structure) will grow in both temporal and spatial directions.
In which direction it grows faster is an important question in the
spatiotemporal stability analysis. In the linear phase, the concepts of linear
absolute (LA) and convective (LC) instabilities were first proposed by Briggs
Briggs (1964), among others, in the study of plasma instabilities. Later, this
methodology was adopted by the researchers in the field of hydrodynamic
instability, most notably Huerre, Monkewitz and Chomaz (Huerre and Monkewitz,
1985, 1990; Chomaz, 2005) among many others. If the initial small-magnitude
wavepacket grows with time but is convected away from the source location (so
that the flow will eventually be stable at all positions), the flow is
convectively unstable, and such disturbance cannot contaminate the upstream
region. If the disturbance is amplified and extends in both upstream and
downstream directions, the flow is called absolutely unstable.

By deriving the dispersion relation from the linearised governing equations,
Carrière & Monkewitz Carriere and Monkewitz (1999) studied the convective and
absolute instabilities in Rayleigh-Bénard-Poiseuille (RBP) convection of a
Newtonian fluid, which is constrained between two infinite parallel plates and
subjected to a thermal buoyancy force and external pressure gradient. In their
work, an impulse disturbance was imposed on the linearised flow (governed by
the linearised Navier-Stokes equations), which mathematically corresponds to
adding a force in the form of the Dirac delta function to the linearised
equations. The solution to the perturbed linearised equations is then their
Green’s function. Using the steepest decent method (Huerre and Monkewitz,
1985), Carrière & Monkewitz can approximate the long-time behaviour of the
Green’s function for each propagation ray and thus were able to determine the
absolute/convective instabilities of the disturbance. Besides, because of the
periodicity of the disturbance in the streamwise and spanwise directions in
their problem, Fourier transform can be applied in these directions and thus
each disturbance mode (with a specified combination of streamwise and spanwise
wavenumbers) can be studied separately for its absolute growth rate. They
found that the disturbance mode reaching zero group velocity at the
convective-absolute transition always corresponds to transverse rolls, while
the disturbance region occupied by the pure longitudinal rolls will remain
convectively unstable for all non-zero Reynolds numbers. This algorithm was
later applied to analyze the linear convective/absolute instabilities in EHD-
Poiseuille flows by Li et al. (Li  _et al._ , 2019). Their results confirmed
the transition from the convective instability to absolute instability in the
region occupied by the transverse rolls when the electric Rayleigh number
increases. The growth rate of oblique rolls is smaller than that of transverse
rolls. These conclusions regarding the linear absolute/convective
instabilities on EHD appear to be similar to those of RBP. In this work, we
will explore the similarities and dissimilarities of EHD-Poiseuille and RBP
flows in the nonlinear regime.

In addition to constructing above dispersion relations to study the
spatiotemporal instability of a flow, the large-time asymptotic response of an
impulse can also be retrieved from the results of numerical simulations by a
post-processing procedure (Delbende  _et al._ , 1998; Delbende and Chomaz,
1998). This method is proper for analysing nonlinear absolute and convective
instabilities of finite-amplitude disturbance. In this case, LA and LC are
replaced by their nonlinear counterparts, namely the nonlinear absolute (NA)
and convective (NC) instabilities which were first systematically studied by
Chomaz Chomaz (1992). If the system relaxes to the basic state everywhere in
the laboratory frame for all initial disturbances of finite extent and finite
amplitude, the instability is nonlinearly convective. While if the system is
not able to relax to the basic state everywhere in the laboratory frame for
some initial condition of finite extent and amplitude, the instability is said
to be nonlinearly absolute. To demonstrate the physical relevance and
significance of these nonlinear concepts, Chomaz studied the one-dimensional
subcritical GLE in a semi-infinite domain, and the hysteresis loop was shown
to be restricted to the nonlinear absolute instability range of the control
parameter Chomaz (1992). Later, the nonlinear convective and absolute
instabilities of parallel plane wakes were studied in the system of two-
dimensional Navier-Stokes equations using highly accurate numerical
simulations (Delbende and Chomaz, 1998). In the subcritical parameter range of
Görtler vortices with the asymptotic suction boundary layer over concave
walls, its instability is proved to be nonlinearly convective by both
theoretical study Park and Huerre (1988) and experiment Chomaz and Perrier
(1991). The front propagation of the linear wavepacket was calculated in this
work, and it was found that the nonlinearity can drive the amplitude of linear
impulse to saturate within the wavepacket as the system evolved from linear to
nonlinear regimes, but the trailing and leading edges of the wavepacket were
not affected, namely the trailing and leading velocities coincide in linear
and nonlinear regimes. From the viewpoint of wave propagation, the
linear/nonlinear convective and absolute instabilities can also be
discriminated from the signs of trailing and leading velocities. If they have
opposite signs, the system is absolutely unstable. Otherwise, it is
convectively unstable. In the end, the fully nonlinear dynamics of the
parallel wakes, including their speed selection of front, absolute and
convective instabilities, and global bifurcations were studied in the work of
Chomaz Chomaz (2003).

###  I.3 The position and structure of current work

In this work, we aim to study the nonlinear convective and absolute
instabilities in the EHD-Poiseuille flow following the methodology developed
by Chomaz and his coworkers in Refs. Delbende  _et al._ (1998); Delbende and
Chomaz (1998); Chomaz (2003). As explained above, knowledge on this idealised
flow can be potentially useful for the industrial applications, such as ESP
Soldati (2003) and EV damper devices Atten and Honda (1982). This idealised
flow retains the most salient flow dynamics in the EHD flow subjected to a
through-flow and at the same time reduces the modelling difficulty in
geometry, which is beneficial for nonlinear analyses. This manuscript is
closely related to and extends our previous works Li et al. Li  _et al._
(2019) on the LA/LC and Zhang Zhang (2016) on the derivation of the GLE for
describing the nonlinear dynamics in EHD-Poiseuille flows. Our extension
focuses on analysing the nonlinear effects in the absolute
instability/convective instability (or in short AI/CI) of EHD-Poiseuille flows
using the numerical simulation and testing the relevance of nonlinear GLE
(derived from a weakly nonlinear expansion of the Navier-Stokes equations and
the Maxwell’s equations in the quasi-electrostatic limit) in depicting the
nonlinear AI/CI, following and supplementing the similar theoretical research
conducted by Suslov & Paolucci Suslov and Paolucci (1997).

The paper is organised as follows. In section II, we introduce the nonlinear
governing equations with boundary conditions and their corresponding
linearisation around the base flow, the framework for the weakly nonlinear
analysis based on the multiple-scale expansion and amplitude expansion
methods, as well as the numerical method used to retrieve the large-time
asymptotics along a certain spatiotemporal ray from snapshots of numerical
simulations. In section III, the numerical code is validated against existing
results in the literature. Our results on the nature of first transition, the
weakly nonlinear analysis, as well as the linear and nonlinear impulse
responses are presented in detail in this section. Finally, in section IV, we
summarise our conclusions and discuss possible future works.

##  II Problem Formulation

###  II.1 Nonlinear equations

The 2D EHD-Poiseuille (or EHD-P) flow that we will study in this work is
confined between two parallel plates with a vertical half height $H^{*}$ and a
horizontal length $L^{*}$ (the superscript asterisk ∗ denotes dimensional
variables) in a $x$-$y$ plane as shown in Fig. 1 (with
$\hat{\mathbf{x}},\hat{\mathbf{y}}$ being unit vectors). A streamwise constant
pressure gradient is imposed to drive the flow in the lateral direction. The
unipolar charges are injected autonomously from the bottom plate with charge
density $Q_{0}^{*}$. An external DC voltage $\Delta\phi_{0}^{*}$ is applied to
the two plates. With the above dimensional variables, as well as the ionic
mobility $K^{*}$, the permittivity $\epsilon^{*}$ and density $\rho_{0}$ of
fluid, we can use the characteristic quantities $H^{*}$,
$H^{*2}/(K^{*}\Delta\phi_{0}^{*})$, $K^{*}\Delta\phi_{0}^{*}/H^{*}$,
$\rho_{0}^{*}K^{*2}\Delta\phi_{0}^{*2}/H^{*2}$, $\Delta\phi_{0}^{*}$,
$\Delta\phi_{0}^{*}/H^{*}$ and $Q_{0}^{*}$ to nondimensionalise the length,
time, velocity, pressure, electric potential, electric field and charge
density, respectively. Then the nondimensionalised governing equations of the
system can be given as

$\displaystyle\nabla\cdot\mathbf{U}=0,$ (1a)
$\displaystyle\frac{\partial\mathbf{U}}{\partial
t}+(\mathbf{U}\cdot\nabla)\mathbf{U}=-\mathbf{\nabla}P+\frac{M^{2}}{T}\nabla^{2}\mathbf{U}+CM^{2}Q\mathbf{E},$
(1b) $\displaystyle\frac{\partial Q}{\partial
t}+\nabla\cdot[(\mathbf{E}+\mathbf{U})Q]=\frac{1}{Fe}\nabla^{2}Q,$ (1c)
$\displaystyle\nabla^{2}\phi=-CQ,$ (1d) $\displaystyle\mathbf{E}=-\nabla\phi,$
(1e) $\displaystyle\text{with \ \ \ \
}M=\frac{(\epsilon^{*}/\rho_{0}^{*})^{1/2}}{K^{*}},\,\,T=\frac{\epsilon^{*}\Delta\phi_{0}^{*}}{K^{*}\mu^{*}},\,\,C=\frac{Q_{0}^{*}H^{*2}}{\Delta\phi_{0}^{*}\epsilon^{*}},\,\,Fe=\frac{K^{*}\Delta\phi_{0}^{*}}{D_{\nu}^{*}},$

where the dielectric liquid is assumed to be incompressible, Newtonian and
isothermal. The magnetic field is neglected in the reduced set of Maxwell’s
equations due to the low electric conductivity of the dielectric liquid
(Castellanos, 1998). Regarding the normalisation method, in the case of
turbulent EHD-Poiseuille flows, previous works Leonard  _et al._ (1983);
Soldati and Banerjee (1998) used the shear velocity as the characteristic
velocity. The current work intends to study the spatiotemporal instability of
EHD-Poiseuille flows, thus the control parameters are not strong enough to
trigger turbulence. The dimensionless parameter $M$ denotes the ratio between
the hydrodynamic mobility $(\epsilon^{*}/\rho_{0}^{*})^{1/2}$ and the ionic
mobility $K^{*}$ with $M<0.1$ for gases and $M$ at least of order one for
liquids (Castellanos and Agrait, 1992). The ratio of the Coulomb force to the
viscous force is described by the electric Rayleigh number $T$, which plays a
similar role as the Rayleigh number $Ra$ in RBC (short for Rayleigh-Bénard
convection). The third parameter $C$ quantifies the injection strength of
charges, and $C>1$ indicates strong injection, including the SCL injection,
which is the scenario considered in this work. The last dimensionless
parameter $Fe$ represents the reciprocal of the charge diffusivity
coefficient, and it typically ranges from $10^{3}$ to $10^{4}$ (Pérez and
Castellanos, 1989). The boundary conditions for field variables
$\mathbf{U}=(U,V)$, $Q$ and $\phi$ are

$\displaystyle y=-1$ $\displaystyle:\;\mathbf{U}=0,\,Q=1,\,\phi=1;$ (2a)
$\displaystyle y=1$ $\displaystyle:\;\mathbf{U}=0,\,\partial Q/\partial
y=0,\,\phi=0;$ (2b) $\displaystyle x=0\;\&\;x=L$
$\displaystyle:\;\text{periodic boundary conditions}.$ (2c)

impulseinjectorcollector$U(y)$

Figure 1: Schematic of the EHD-Poiseuille flow with an initially localized
impulse in the dielectric fluid which is confined between two infinite
parallel plates. Charged ions are injected into the fluid from bottom plate
(i.e. injector), and they are collected by top plate (i.e. collector).

A streamwise constant pressure gradient is provided to sustain the base
Poiseuille flow $U(y)=\mathcal{U}(1-y^{2})$ (with $y\in[-1,1]$), where
$\mathcal{U}$ denotes the maximum velocity magnitude of the Poiseuille flow,
and it can be used to control the strength of the through-flow. As $U(y)$ is
the steady solution to the momentum equation (1b) with no-slip boundary
conditions, the constant pressure gradient can be calculated as

$U(y)=-\frac{1}{2}\frac{dP}{dx}\frac{T}{M^{2}}\left(1-y^{2}\right),\,\
\frac{dP}{dx}=-2\mathcal{U}\frac{M^{2}}{T}.$ (3)

To evaluate the electric transfer in EHD-Poiseuille flows, the following
electric Nusseult number $Ne$ (Zhang, 2016) is defined

$Ne=\frac{I_{e}}{I_{0}},\;I_{e}=\iint_{\mathcal{V}}\bigg{[}Q(V+E_{y})-\frac{1}{Fe}\frac{\partial
Q}{\partial y}\bigg{]}dxdy,$ (4)

where $Ne$ denotes the ratio between total current $I_{e}$ (with convective
motions, i.e., $\mathbf{U}\neq 0$) and the conductive current $I_{0}$ (without
convective motions, i.e., $\mathbf{U}=0$).

###  II.2 Perturbative form of the governing equations

In order to study the linear and nonlinear instabilities in EHD-Poiseuille
flows, the field variables are decomposed into a summation of base states
(which are steady solutions to Eqs. 1) and perturbations, i.e.
$\mathbf{U}=\bar{\mathbf{U}}+\mathbf{u}$, $P=\bar{P}+p$,
$\mathbf{E}=\bar{\mathbf{E}}+\mathbf{e}$, $Q=\bar{Q}+q$ and
$\phi=\bar{\phi}+\varphi$, which can be substituted into the nonlinear
governing equations (1) to yield

$\displaystyle\nabla\cdot\mathbf{u}=0,$ (5a)
$\displaystyle\frac{\partial\mathbf{u}}{\partial
t}+(\mathbf{u}\cdot\nabla)\bar{\mathbf{U}}+(\bar{\mathbf{U}}\cdot\nabla)\mathbf{u}=-\nabla
p+\frac{M^{2}}{T}\nabla^{2}\mathbf{u}+CM^{2}(q\bar{\mathbf{E}}+\bar{Q}\mathbf{e})+\mathbf{N}_{\mathbf{u}},$
(5b) $\displaystyle\frac{\partial q}{\partial
t}+\nabla\cdot[(\bar{\mathbf{E}}+\bar{\mathbf{U}})q+(\mathbf{e}+\mathbf{u})\bar{Q}]=\frac{1}{Fe}\nabla^{2}q+N_{q},$
(5c) $\displaystyle\nabla^{2}\varphi=-Cq,$ (5d)
$\displaystyle\mathbf{e}=-\nabla\varphi,$ (5e)

where the perturbation velocity $\mathbf{u}=(u,v)$ and electric field
$\mathbf{e}=(e_{x},e_{y})$,
$\mathbf{N}_{\mathbf{u}}=-(\mathbf{u}\cdot\nabla)\mathbf{u}+CM^{2}(q\mathbf{e})$
and $N_{q}=-\nabla\cdot[(\mathbf{e}+\mathbf{u})q]$ are nonlinear terms of
perturbations. Discarding the nonlinear terms directly leads to the linearised
governing equations for the perturbations. $\bar{\mathbf{E}}$ and $\bar{Q}$
refer to hydrostatic solutions of the electroconvective flow (see, e.g., Ref.
Zhang  _et al._ (2015) for their profiles), and as mentioned above
$\bar{\mathbf{U}}=\mathcal{U}(1-y^{2})\hat{\mathbf{x}}$ is the Poiseuille
flow. The homogeneous boundary conditions are applied for the perturbed fields

$\displaystyle y=-1$ $\displaystyle:\;\mathbf{u}=0,\,q=0,\,\varphi=0;$ (6a)
$\displaystyle y=1$ $\displaystyle:\;\mathbf{u}=0,\,\partial q/\partial
y=0,\,\varphi=0;$ (6b) $\displaystyle x=0\;\&\;x=L$
$\displaystyle:\;\text{periodic boundary conditions}.$ (6c)

###  II.3 Numerical analysis of flow instability

####  II.3.1 Temporal instability

The temporal instability analysis of a flow concerns the development of the
disturbance on the flow in time. Time series of the disturbance can be
obtained either experimentally or numerically, based on which, useful
information, such as the most dangerous wavelength/wavenumbers and their
associated growth rates, can be calculated in the post-processing analysis
(Delbende and Chomaz, 1998). In order to get these pieces of spectral
information, Fast Fourier Transform (FFT) is conventionally applied; for
example, in our numerical simulations to be presented, FFT is used to
decompose the wavepacket with a broad spectrum structure into multiple waves
within the linear regime of the flow evolution. The growth rate of each wave
can be calculated using the temporal evolution of the total energy density; in
the EHD-Poiseuille flow to be studied, the norm of total energy in the
spectral space is defined as

$\tilde{\mathcal{E}}(\alpha,t)=\left[\frac{1}{2}\int_{-1}^{1}\left(\mathbf{u}(\alpha,y,t)^{2}+M^{2}\mathbf{e}(\alpha,y,t)^{2}\right)dy\right]^{1/2},$
(7)

at streamwise wavenumber $\alpha$. The total energy density of EHD-Poiseuille
flows consists of kinetic energy density $\mathcal{E}_{k}=\mathbf{u}^{2}/2$
and electric energy density $\mathcal{E}_{\varphi}=M^{2}\mathbf{e}^{2}/2$
(Castellanos, 1998; Zhang, 2016). When solving the nonlinear equations (1), we
use the perturbed velocity and electric field in equation (7) by subtracting
the base states from the nonlinear fields (i.e.,
$\mathbf{u}=\mathbf{U}-\bar{\mathbf{U}}$ and
$\mathbf{e}=\mathbf{E}-\bar{\mathbf{E}}$). In the linear regime,
$\tilde{\mathcal{E}}(\alpha,t)$ presents an exponential scaling at an
asymptotically large time according to the linear theory, that is,
$\tilde{\mathcal{E}}(\alpha,t)\propto e^{\omega(\alpha)t},$ with
$\omega(\alpha)\sim\frac{\partial}{\partial
t}\ln\tilde{\mathcal{E}}(\alpha,t)$ as $t\rightarrow\infty$. Thus, following
Ref. Delbende and Chomaz (1998), the leading growth rate at wavenumber
$\alpha$ can be computed from the norm of total energy density by
$\omega(\alpha)\approx\frac{\ln\left[\tilde{\mathcal{E}}(\alpha,t_{2})/\tilde{\mathcal{E}}(\alpha,t_{1})\right]}{t_{2}-t_{1}}$
where the wavenumber is constrained in the range $-\pi/\delta
x<\alpha<\pi/\delta x$ due to the spatial discretization ($\delta x$ refers to
the interval of adjacent grid points along the streamwise direction), and the
smallest increment of wavenumber is $\delta\alpha=2\pi/L$. The linear
stability theory also dictates that the flow field at an asymptotically large
time is proportional to the eigenvector. With this fact, we can verify our
calculation by checking the invariance of the quantity
$\tilde{\mathcal{E}}_{f}(\alpha,y,t)=\frac{\hat{\tilde{\mathcal{E}}}(\alpha,y,t)}{\mathcal{E}(\alpha,t)}$,
where the eigenvector $\hat{\tilde{\mathcal{E}}}(\alpha,y,t)$ at streamwise
wavenumber $\alpha$ is normalized by its amplitude
$\tilde{\mathcal{E}}(\alpha,t)$.

####  II.3.2 Spatiotemporal instability

The temporal stability analysis only concerns the development of the
disturbance in time. When a wave is propagating in a medium, its development
in space should also be studied. Considering both the spatial and temporal
development of the disturbance helps to differentiate two types of instability
mechanisms, i.e., convective instability and absolute instability. One can
characterise the spatiotemporal development of an unstable mode in the
physical space by the temporal growth rate $\omega$ and the spatial group
velocity $v_{g}$; the latter refers to the propagation speed of the overall
envelope of the wavepacket. As mentioned in the introduction section,
theoretical methods construct a dispersion equation to probe the asymptotic
spatiotemporal dynamics (see Ref. Huerre and Monkewitz (1985) for a general
procedure or our previous work Li  _et al._ (2019) for a more relevant
reference) and post-processing methods analyse the numerical data to extract
the spatiotemporal information of a flow. In this work, we adopt the latter.
To unambiguously extract the amplitude and the phase of a wavepacket, the
Hilbert transform is applied to the flow variables (Brancher and Chomaz,
1997), which reads

$\tilde{g}(x,y,t)=\hat{\tilde{g}}(x,y,t)e^{i\psi(x,y,t)},$ (8)

where $\hat{\tilde{g}}(x,y,t)$ represents the complex-valued amplitude
function of any flow variable, and $\psi(x,y,t)$ is the phase of the
wavepacket. The real amplitude function $\hat{\tilde{\mathcal{E}}}(x,t)$ of
energy norm can be obtained by integration along the $y$-axis

$\hat{\tilde{\mathcal{E}}}(x,t)=\left[\frac{1}{2}\int_{-1}^{1}\left(\hat{\tilde{\mathbf{u}}}(x,y,t)^{2}+M^{2}\hat{\tilde{\mathbf{e}}}(x,y,t)^{2}\right)dy\right]^{1/2}.$
(9)

It can be shown that the real amplitude function at an asymptotically large
time follows (see Ref. (Huerre and Monkewitz, 1985) for the derivation)

$\hat{\tilde{\mathcal{E}}}(x,t)\propto t^{-1/2}e^{\sigma(v_{g})t}\ \ \ \
\text{with}\ \ \ \,v_{g}=(x-x_{0})/t=\text{const},\,t\rightarrow\infty,$ (10)

where $x_{0}$ is the initial streamwise position of the impulse disturbance
and $\sigma(v_{g})$ refers to the dominant growth rate along the ray $v_{g}$.
Based on snapshots extracted within the linear regime, the growth rate can be
computed in the following discretized form

$\sigma(v_{g})\approx\frac{\ln\left[\hat{\tilde{\mathcal{E}}}(v_{g}t_{2},t_{2})/\hat{\tilde{\mathcal{E}}}(v_{g}t_{1},t_{1})\right]}{t_{2}-t_{1}}+\sigma_{0}(t_{1},t_{2}),\
\,\sigma_{0}(t_{1},t_{2})=\frac{\ln(t_{2}/t_{1})}{2(t_{2}-t_{1})},$ (11)

where $\sigma_{0}(t_{1},t_{2})$ is a finite-time correction for the
spatiotemporal growth rate, which arises from the $t^{-1/2}$ factor in
equation (10) (Delbende and Chomaz, 1998).

In some numerical instances to be followed, the nonlinear equations (1) are
solved to probe the spatiotemporal instability of the EHD-Poiseuille flows. It
has been noticed in Refs. Melville (1983); Delbende and Chomaz (1998) and we
also observed that the Hilbert transform (8) cannot accurately extract the
amplitude function and the phase of nonlinearly unstable wavepackets. In this
case, the enstrophy is used in the nonlinear analysis, following Refs.
Melville (1983); Delbende and Chomaz (1998),

$\eta(x,t)=\left(\int_{-1}^{1}\omega^{2}(x,y,t)dy\right)^{1/2},$ (12)

where $\omega=\partial v/\partial x-\partial u/\partial y$ denotes the
vorticity in 2D flows.

In the end, since the spatiotemporal instability properties of a flow can be
inferred from its response to an impulse disturbance, some physical
understanding of the impulse response is instructive. Edge velocities
$v_{\pm}$ refer to those at which neutral waves propagate (i.e.
$\sigma(v_{\pm})=0$), see Fig. 2(c) to follow, and can help discriminate
between linearly convective and absolute instabilities. The leading velocity
is always positive $v_{+}>0$ as the Poiseuille flow moves downwards along the
$x$-axis. Thus, if we have $v_{-}>0$ and the spatiotemporal ray $v_{g}=0$ lies
outside the range of exponentially growing modes (i.e., $\sigma(v_{g}=0)<0$),
the flow is convective instability; otherwise, if $v_{-}<0$ and the
spatiotemporal ray $v_{g}=0$ lies in the range of exponentially growing modes
(i.e., $\sigma(v_{g}=0)>0$), the flow is absolutely unstable whenever it is
unstable.

###  II.4 Weakly nonlinear stability analysis

In this section, we will present the theoretical derivation of the GLE in the
weakly nonlinear stability analysis. It can be used as a reduced-order model
to study the spatiotemporal development of a flow. In two-dimensional EHD
problems, by using the streamfunction formulation $\psi$ (with
$u=\frac{\partial\psi}{\partial y}$ and $v=-\frac{\partial\psi}{\partial x}$),
eliminating the pressure term, and expressing the variables related to the
electric field in terms of $\varphi$, we can recast the equation system (5)
into the following two-variable form Zhang (2016)

$\displaystyle\frac{\partial\nabla^{2}\psi}{\partial t}=$
$\displaystyle-\bar{U}\frac{\nabla^{2}\psi}{\partial
x}+\bar{U}^{\prime\prime}\frac{\partial\psi}{\partial
x}+\frac{M^{2}}{T}\nabla^{4}\psi-M^{2}\left(\bar{\phi}^{\prime}\frac{\partial\nabla^{2}\varphi}{\partial
x}-\bar{\phi}^{\prime\prime\prime}\frac{\partial\varphi}{\partial x}\right)$
$\displaystyle-\left(\frac{\partial\psi}{\partial
y}\frac{\partial\nabla^{2}\psi}{\partial x}-\frac{\partial\psi}{\partial
x}\frac{\partial\nabla^{2}\psi}{\partial
y}\right)-M^{2}\left(\frac{\partial\varphi}{\partial
y}\frac{\partial\nabla^{2}\varphi}{\partial x}-\frac{\partial\varphi}{\partial
x}\frac{\partial\nabla^{2}\varphi}{\partial y}\right),$ (13a)
$\displaystyle\frac{\partial\nabla^{2}\varphi}{\partial t}=$
$\displaystyle\bar{\phi}^{\prime\prime\prime}\left(\frac{\partial\varphi}{\partial
y}+\frac{\partial\psi}{\partial
x}\right)+\left(\bar{\phi}^{\prime}\frac{\partial}{\partial
y}-\bar{U}\frac{\partial}{\partial
x}\right)\nabla^{2}\varphi+2\bar{\phi}^{\prime\prime}\nabla^{2}\varphi+\frac{1}{Fe}\nabla^{4}\varphi$
$\displaystyle+\left(\frac{\partial\varphi}{\partial
x}-\frac{\partial\psi}{\partial
y}\right)\frac{\partial\nabla^{2}\varphi}{\partial
x}+\left(\frac{\partial\varphi}{\partial y}+\frac{\partial\psi}{\partial
x}\right)\frac{\partial\nabla^{2}\varphi}{\partial
y}+\nabla^{2}\varphi\nabla^{2}\varphi,$ (13b)

and the boundary conditions are periodic in $x$ direction and $\psi(\pm
1)=\psi^{\prime}(\pm 1)=0$ and $\varphi(\pm
1)=\varphi^{\prime\prime}(1)=\varphi^{\prime\prime\prime}(-1)=0$ in $y$
direction (and prime ′ means the derivative in $y$ direction). This equation
system can be further expressed in a compact form as

$\bm{M}\frac{\partial\bm{\gamma}}{\partial t}=\bm{L}\bm{\gamma}+\bm{N},\ \
\text{or}\ \ \left(\bm{M}\frac{\partial}{\partial
t}-\bm{L}\right)\bm{\gamma}=\bm{N},$ (14)

where $\bm{\gamma}=(\psi,\varphi)^{T}$, $\bm{M}$ is the weight matrix,
$\bm{L}$ is the linear operator, $\bm{N}$ is the nonlinear operator, whose
explicit expressions can be deduced straightforwardly by matching with
equation (13), see the appendix A and Ref. Zhang (2016).

There are various ways to perform a weakly nonlinear analysis and the two most
classical methods are the multiple-scale expansion Stuart (1960) around the
linear critical condition and the amplitude expansion Watson (1960) which can
be valid even away from the linear critical condition. They have been proved
to be equivalent to each other at critical conditions (Fujimura, 1989). In the
following section, we will briefly describe the key steps in the two
perturbation methods.

####  II.4.1 Multiple-scale expansion method

Our multiple-scale expansion method follows the schemes in Stewartson & Stuart
Stewartson and Stuart (1971) and Fujimura Fujimura (1989), where the time $t$,
streamwise spatial coordinate $x$, perturbation $\bm{\gamma}$ and the
governing parameter $T$ are expanded in series of a small quantity $\epsilon$
as

$\displaystyle\frac{\partial}{\partial t}$
$\displaystyle=\frac{\partial}{\partial
t_{0}}+\epsilon\frac{\partial}{\partial
t_{1}}+\epsilon^{2}\frac{\partial}{\partial t_{2}}+O(\epsilon^{3}),\ \ \ \
\frac{\partial}{\partial x}=\frac{\partial}{\partial
x_{0}}+\epsilon\frac{\partial}{\partial x_{1}}+O(\epsilon^{2}),$ (15a)
$\displaystyle\bm{\gamma}$
$\displaystyle=\epsilon\bm{\gamma}_{1}+\epsilon^{2}\bm{\gamma}_{2}+\epsilon^{3}\bm{\gamma}_{3}+O(\epsilon^{4}),\
\ \ \ \ \ \ \ \ T=T_{c}+\epsilon^{2}+O(\epsilon^{4}),$ (15b)

where $\epsilon^{2}=T-T_{c}$ is a measure of the value $T$ away from the
critical $T_{c}$. Accordingly, the operators introduced in equation (14)
should also be expanded

$\bm{M}=\bm{M}_{0}+\epsilon\bm{M}_{1}+\epsilon^{2}\bm{M}_{2}+O(\epsilon^{3}),\
\ \ \
\bm{L}=\bm{L}_{0}+\epsilon\bm{L}_{1}+\epsilon^{2}\bm{L}_{2}+O(\epsilon^{3}),\
\ \ \bm{N}=\epsilon^{2}\bm{N}_{2}+\epsilon^{3}\bm{N}_{3}+O(\epsilon^{4}).$
(16)

Substituting equations (15) and (16) into equation (14) and equating terms on
the same orders of $\epsilon$ give

$\displaystyle\left(\bm{M}_{0}\frac{\partial}{\partial
t_{0}}-\bm{L}_{0}\right)\bm{\gamma}_{1}=0,$ (17a)
$\displaystyle\left(\bm{M}_{0}\frac{\partial}{\partial
t_{0}}-\bm{L}_{0}\right)\bm{\gamma}_{2}=\left(\bm{L}_{1}-\bm{M}_{1}\frac{\partial}{\partial
t_{0}}-\bm{M}_{0}\frac{\partial}{\partial
t_{1}}\right)\bm{\gamma}_{1}+\bm{N}_{2},$ (17b)
$\displaystyle\left(\bm{M}_{0}\frac{\partial}{\partial
t_{0}}-\bm{L}_{0}\right)\bm{\gamma}_{3}=\left(\bm{L}_{2}-\bm{M}_{2}\frac{\partial}{\partial
t_{0}}-\bm{M}_{1}\frac{\partial}{\partial
t_{1}}-\bm{M}_{0}\frac{\partial}{\partial
t_{2}}\right)\bm{\gamma}_{1}+\left(\bm{L}_{1}-\bm{M}_{1}\frac{\partial}{\partial
t_{0}}-\bm{M}_{0}\frac{\partial}{\partial
t_{1}}\right)\bm{\gamma}_{2}+\bm{N}_{3}.$ (17c)

Once we assume the first-order solution $\bm{\gamma}_{1}$ to be normal-mode,
the solutions $\bm{\gamma}_{2}$ and $\bm{\gamma}_{3}$ can also be deduced and
they are

$\displaystyle\bm{\gamma}_{1}=A_{1}(x_{1},t_{1},t_{2})\tilde{\bm{\gamma}}_{1}(y)e^{i\alpha
x_{0}+\mu t_{0}}+c.c.,$ (18a)
$\displaystyle\bm{\gamma}_{2}=A_{1}^{*}A_{1}\tilde{\bm{\gamma}}_{20}(y)+\left[\frac{\partial
A_{1}}{\partial x_{1}}\tilde{\bm{\gamma}}_{21}(y)e^{i\alpha x_{0}+\mu
t_{0}}+c.c.\right]+\left[A_{1}^{2}\tilde{\bm{\gamma}}_{22}(y)e^{2i\alpha
x_{0}+2\mu t_{0}}+c.c.\right],$ (18b)
$\displaystyle\bm{\gamma}_{3}=\left[A_{1}\tilde{\bm{\gamma}}_{31,1}(y)e^{i\alpha
x_{0}+\mu t_{0}}+c.c.\right]+\left[\frac{\partial^{2}A_{1}}{\partial
x_{1}^{2}}\tilde{\bm{\gamma}}_{31,2}(y)e^{i\alpha x_{0}+\mu
t_{0}}+c.c.\right]+\left[A_{1}|A_{1}|^{2}\tilde{\bm{\gamma}}_{31,3}(y)e^{i\alpha
x_{0}+\mu t_{0}}+c.c.\right]+\cdots,$ (18c)

where $A_{1}$ is the complex amplitude of the perturbation, $\alpha$ the
streamwise wavenumber, $\mu$ the complex frequency with its real part being
the linear growth rate, and $c.c.$ the complex conjugate of the preceding one
term.

Evolution equations of $A_{1}$ can be derived by enforcing the solvability
condition Bender and Orszag (2013) on equation (17b and 17c), i.e., the
projection of the secular terms on the right hand side of these two equations
on the adjoint eigenvector $\bm{\gamma}_{1}^{\dagger}$ (as introduced below in
equation (20)) is zero. Following Ref. Zhang (2016), we introduce the inner
product

$\langle\bm{f},\bm{g}\rangle=\frac{1}{L_{p}}\int_{0}^{L_{p}}\int_{-1}^{1}\bm{f}\cdot\bm{M}_{\text{w}}\bm{g}\,dy\,dx_{0}={\langle\tilde{\bm{f}},\tilde{\bm{g}}\rangle}_{s}+c.c.,$
(19)

where $L_{p}=2\pi/\alpha$ is the wavelength,
$\bm{M}_{\text{w}}=\begin{bmatrix}1&0\\\ 0&\kappa^{2}\end{bmatrix}$ is the
weight matrix with $\kappa=M$. $\kappa$ also appears in the adjoint problem
resulting from the integration by part process below.
$\langle\cdot,\cdot\rangle_{s}$ is the corresponding inner product defined in
the spectral space. Then the adjoint problem corresponding to equation (17a)
can be defined through integration by part (Luchini and Bottaro, 2014) based
on the introduced inner product as

$\left\langle\left(\bm{M}_{0}\frac{\partial}{\partial
t_{0}}-\bm{L}_{0}\right)\bm{\gamma}_{1},\bm{\gamma}_{1}^{\dagger}\right\rangle=\left\langle\bm{\gamma}_{1},\left(\bm{M}_{0}^{\dagger}\frac{\partial}{\partial
t_{0}}-\bm{L}_{0}^{\dagger}\right)\bm{\gamma}_{1}^{\dagger}\right\rangle,\ \ \
\rightarrow\ \ \ \left(\bm{M}_{0}^{\dagger}\frac{\partial}{\partial
t_{0}}-\bm{L}_{0}^{\dagger}\right)\bm{\gamma}_{1}^{\dagger}=0.$ (20)

Here the explicit expressions of the adjoint operators $\bm{M}_{0}^{\dagger}$
and $\bm{L}_{0}^{\dagger}$ are given in the appendix A. This adjoint problem
is also an eigenvalue problem and it is solved similarly to the linear
eigenvalue problem in equation (17a) at the same parameter settings. Then the
projection (see Ref. Zhang (2016) for the detailed illustration of the
projection process) to eliminate secular terms leads to the following equation
involving the group velocity $c_{g}$

$\frac{\partial A_{1}}{\partial t_{1}}+c_{g}\frac{\partial A_{1}}{\partial
x_{1}}=0\ \ \ \text{with}\ \ \
c_{g}=-\frac{{\langle(\tilde{\bm{L}}_{1}^{\circ}-\mu_{c}\tilde{\bm{M}}_{1}^{\circ})\tilde{\bm{\gamma}}_{1},\tilde{\bm{\gamma}}_{1}^{\dagger}\rangle}_{s}}{{\langle\tilde{\bm{M}}_{0}\tilde{\bm{\gamma}}_{1},\tilde{\bm{\gamma}}_{1}^{\dagger}\rangle}_{s}},$
(21)

and the Ginzburg-Landau equation (GLE)

$\frac{\partial A_{1}}{\partial
t_{2}}=a_{1}A_{1}+a_{2}\frac{\partial^{2}A_{1}}{\partial
x_{1}^{2}}+a_{3}{|A_{1}|}^{2}A_{1},$ (22)

where the coefficients are

$a_{1}=\frac{{\langle\tilde{\bm{L}}_{2MT}\tilde{\bm{\gamma}}_{1},\tilde{\bm{\gamma}}_{1}^{\dagger}\rangle}_{s}}{{\langle\tilde{\bm{M}}_{0}\tilde{\bm{\gamma}}_{1},\tilde{\bm{\gamma}}_{1}^{\dagger}\rangle}_{s}},a_{2}=\frac{{\langle(\tilde{\bm{L}}_{2}^{\circ}-\mu_{c}\tilde{\bm{M}}_{2}^{\circ}+c_{g}\tilde{\bm{M}}_{1}^{\circ})\tilde{\bm{\gamma}}_{1}+(\tilde{\bm{L}}_{1}^{\circ}-\mu_{c}\tilde{\bm{M}}_{1}^{\circ}+c_{g}\tilde{\bm{M}}_{0})\tilde{\bm{\gamma}}_{21},\tilde{\bm{\gamma}}_{1}^{\dagger}\rangle}_{s}}{{\langle\tilde{\bm{M}}_{0}\tilde{\bm{\gamma}}_{1},\tilde{\bm{\gamma}}_{1}^{\dagger}\rangle}_{s}},a_{3}=\frac{{\langle\tilde{\bm{N}}_{3}^{\circ},\tilde{\bm{\gamma}}_{1}^{\dagger}\rangle}_{s}}{{\langle\tilde{\bm{M}}_{0}\tilde{\bm{\gamma}}_{1},\tilde{\bm{\gamma}}_{1}^{\dagger}\rangle}_{s}},$
(23)

where the symbol tilde $\tilde{}$ marks operators and eigenvectors in the
spectral space, $\mu_{c}$ is the frequency at the critical condition. All the
explicit expressions of these operators are provided in the appendix.

####  II.4.2 Amplitude expansion method

Amplitude expansion method was first used in the fluid mechanics community by
Watson Watson (1960) to study weakly nonlinear flows and has recently been
refined by Pham & Suslov Pham and Suslov (2018). The following version of the
amplitude expansion method is due to the latter. When the parameter of
interest is away from the linear critical condition, the above multiple-scale
expansion method is not guaranteed to be valid as the small quantify
$\epsilon$ introduced in equation (15) may be too large for the truncated
power series to be converged. In the amplitude expansion method, the magnitude
of the perturbation itself is taken as the small quantity, based on which the
expansion reads

$\bm{\gamma}=[A\tilde{\bm{\gamma}}_{11}E+c.c.]+|A|^{2}\tilde{\bm{\gamma}}_{20}+[\frac{\partial
A}{\partial
z}\tilde{\bm{\gamma}}_{21}E+c.c.]+[A^{2}\tilde{\bm{\gamma}}_{22}E^{2}+c.c.]+[A|A|^{2}\tilde{\bm{\gamma}}_{31}E+c.c.]+[\frac{\partial^{2}A}{\partial
z^{2}}\tilde{\bm{\gamma}}_{31z}E+c.c.]+\text{H. O. T.},$ (24)

where $A$ is the complex amplitude of the disturbance and $E=\exp(i\alpha x)$
(not to be confused with the symbol for the electric field $\mathbf{E}$).
Conventionally, in previous works (Watson, 1960; Fujimura, 1989; Pham and
Suslov, 2018; Cudby and Lefauve, 2021), $A$ is assumed to be only a function
of time. As suggested by one of the reviewers, we also consider its spatial
dependence for the analysis of spatiotemporal instability in the amplitude
expansion scheme. The only requirement for the above expansion to be valid is
that the disturbance amplitude is sufficiently small (even though the
governing parameters may be at a finite distance away from the linear critical
conditions). To facilitate the description of the solution process, we define
the following operators

$\tilde{\bm{M}}_{q\alpha}\tilde{\bm{\gamma}}_{pq}E^{q}\equiv\bm{M}(\tilde{\bm{\gamma}}_{pq}E^{q}),\tilde{\bm{L}}_{q\alpha}\tilde{\bm{\gamma}}_{pq}E^{q}\equiv\bm{L}(\tilde{\bm{\gamma}}_{pq}E^{q}),\mathcal{L}_{q\alpha,\sigma}\equiv\sigma\tilde{\bm{M}}_{q\alpha}-\tilde{\bm{L}}_{q\alpha},$
(25)

where subscripts $pq$ denote a mode with temporal amplitude of order $p$
($A^{p}$) and wavenumber $q\alpha$ (note that $p,q$ are integers). Then the
following Ginzburg-Landau equation describing the temporal evolution of the
amplitude arises

$\frac{\partial A}{\partial t}=-c_{gA}\frac{\partial A}{\partial x}+\sigma
A+a_{2A}\frac{\partial^{2}A}{\partial x^{2}}+KA|A|^{2}+\cdots.$ (26)

Substituting equations (26) and (24) into equation (14), and collecting the
like terms having the same order of amplitude and wavenumber, one can obtain

$\displaystyle\mathcal{L}_{\alpha,\sigma}\tilde{\bm{\gamma}}_{11}$
$\displaystyle=0,$ (27a)
$\displaystyle\mathcal{L}_{0,2\sigma_{r}}\tilde{\bm{\gamma}}_{20}$
$\displaystyle=\tilde{\bm{N}}_{20},$ (27b)
$\displaystyle\mathcal{L}_{\alpha,\sigma}\tilde{\bm{\gamma}}_{21}$
$\displaystyle=c_{gA}\tilde{\bm{M}}_{1\alpha}\tilde{\bm{\gamma}}_{11}+\tilde{\bm{n}}_{21},$
(27c) $\displaystyle\mathcal{L}_{2\alpha,2\sigma}\tilde{\bm{\gamma}}_{22}$
$\displaystyle=\tilde{\bm{N}}_{22},$ (27d)
$\displaystyle\mathcal{L}_{\alpha,\sigma}\tilde{\bm{\gamma}}_{31z}$
$\displaystyle=-a_{2A}\tilde{\bm{M}}_{1\alpha}\tilde{\bm{\gamma}}_{11}+\tilde{\bm{n}}_{31z},$
(27e)
$\displaystyle\mathcal{L}_{\alpha,\sigma+2\sigma_{r}}\tilde{\bm{\gamma}}_{31}$
$\displaystyle=-K\tilde{\bm{M}}_{1\alpha}\tilde{\bm{\gamma}}_{11}+\tilde{\bm{N}}_{31},$
(27f)

where the inhomogeneous terms
$\tilde{\bm{n}}_{21}=(\tilde{\bm{L}}_{1\alpha}^{\circ}-\sigma\tilde{\bm{M}}_{1\alpha}^{\circ})\tilde{\bm{\gamma}}_{11}$
and
$\tilde{\bm{n}}_{31z}=(\tilde{\bm{L}}_{2\alpha}^{\circ}-\sigma\tilde{\bm{M}}_{2\alpha}^{\circ}+c_{gA}\tilde{\bm{M}}_{1\alpha}^{\circ})\tilde{\bm{\gamma}}_{11}+(\tilde{\bm{L}}_{1\alpha}^{\circ}-\sigma\tilde{\bm{M}}_{1\alpha}^{\circ}+c_{gA}\tilde{\bm{M}}_{0\alpha})\tilde{\bm{\gamma}}_{21}$.
Explicit expressions of all the operators appearing here are given in the
appendix. Among these six equations,
$\mathcal{L}_{\alpha,\sigma}\tilde{\bm{\gamma}}_{11}=0$ is exactly the linear
eigenvalue problem corresponding to equation (17a). For $T>T_{c}$ near the
linear critical conditions, the linear growth rate $\sigma_{r}>0$, and thus,
the operators $\mathcal{L}_{0,2\sigma_{r}}$ and
$\mathcal{L}_{2\alpha,2\sigma}$ are non-singular (even when $T\to T_{c}$ where
$\sigma_{r}\to 0$), making the two corresponding equations readily solvable
Pham and Suslov (2018); Cudby and Lefauve (2021). It should be noted that when
$T<T_{c}$ (with $\sigma_{r}<0$), there will be resonances among decaying modes
and special treatments are required (Suslov, 2008; Pham and Suslov, 2018).
This is not a problem in the present analysis where we consider non-negative
$\sigma_{r}$. Near the linear critical condition (where the linear growth rate
$\sigma_{r}\to 0$), the operator
$\mathcal{L}_{\alpha,\sigma+2\sigma_{r}}\to\mathcal{L}_{\alpha,\sigma}$, which
is singular by construction. In this case, for the third-order solution, it
should be noted that if $\tilde{\bm{\gamma}}_{31}$ is a solution, following
the linear problem $\mathcal{L}_{\alpha,\sigma}\tilde{\bm{\gamma}}_{11}=0$, we
know that $\tilde{\bm{\gamma}}_{31}+C\tilde{\bm{\gamma}}_{11}$ will also be a
solution with $C$ being an arbitrary complex number. In this case, a
solvability condition resorting to the adjoint problem (as described in the
multiple-scale expansion) along with a normalization condition (Fujimura,
1989), is conventionally enforced to eliminate the ambiguity. However, away
from the linear critical condition, an orthogonality condition can be utilised
to make $\tilde{\bm{\gamma}}_{31}$ and $K$ uniquely determined Pham and Suslov
(2018), that is
$\langle\tilde{\bm{\gamma}}_{11},\tilde{\bm{\gamma}}_{31}\rangle_{\mathcal{M}}\equiv\tilde{\bm{\gamma}}_{11}^{H}\mathcal{M}\tilde{\bm{\gamma}}_{31}=0$,
where the superscript $H$ denotes conjugate (Hermitian) transpose and
$\mathcal{M}$ is a positive definite Hermitian matrix which in the present
work is the identity matrix $I$ (see Ref. Pham and Suslov (2018) for a general
discussion on the choice of the matrix $\mathcal{M}$). Such a simple choice
has the advantage of using all the information in the fundamental modes
$\tilde{\bm{\gamma}}_{11}$ and $\tilde{\bm{\gamma}}_{31}$.

With the above orthogonality condition defined, following the solution
procedure proposed in Crouch and Herbert (1993) and formally proved in the
recent study Pham and Suslov (2018), the extended equation system can be
arrived

$\begin{bmatrix}\mathcal{L}_{\alpha,\sigma+2\sigma_{r}}&\tilde{\bm{M}}_{1\alpha}\tilde{\bm{\gamma}}_{11}\\\
\tilde{\bm{\gamma}}_{11}^{H}\mathcal{M}&0\end{bmatrix}\begin{bmatrix}\tilde{\bm{\gamma}}_{31}\\\
K\end{bmatrix}=\begin{bmatrix}\tilde{\bm{N}}_{31}\\\ 0\end{bmatrix},$ (28)

which combines the solution of $\tilde{\bm{\gamma}}_{31}$ and $K$ into a
single problem and can be reliably solved away from the linear critical
condition. Similarly, to determine the coefficients $c_{gA}$ and $a_{2A}$ in
equation (27c,e), we also apply the orthogonality conditions
$\langle\tilde{\bm{\gamma}}_{11},\tilde{\bm{\gamma}}_{21}\rangle_{\mathcal{M}}=0$
and
$\langle\tilde{\bm{\gamma}}_{11},\tilde{\bm{\gamma}}_{31z}\rangle_{\mathcal{M}}=0$.
Thus, two more extended equation systems can be formulated

$\begin{bmatrix}\mathcal{L}_{\alpha,\sigma}&-\tilde{\bm{M}}_{1\alpha}\tilde{\bm{\gamma}}_{11}\\\
\tilde{\bm{\gamma}}_{11}^{H}\mathcal{M}&0\end{bmatrix}\begin{bmatrix}\tilde{\bm{\gamma}}_{21}\\\
c_{gA}\end{bmatrix}=\begin{bmatrix}\tilde{\bm{n}}_{21}\\\
0\end{bmatrix},\begin{bmatrix}\mathcal{L}_{\alpha,\sigma}&\tilde{\bm{M}}_{1\alpha}\tilde{\bm{\gamma}}_{11}\\\
\tilde{\bm{\gamma}}_{11}^{H}\mathcal{M}&0\end{bmatrix}\begin{bmatrix}\tilde{\bm{\gamma}}_{31z}\\\
a_{2A}\end{bmatrix}=\begin{bmatrix}\tilde{\bm{n}}_{31z}\\\ 0\end{bmatrix},$
(29)

As reported in Ref. Cudby and Lefauve (2021), the linear operators here may be
poorly conditioned near the critical point even if it is nonsingular but they
found that it still can be inverted to give converged results. This is also
the case in the present work.

####  II.4.3 Relation of the two expansion methods

The GLE (22) derived in the framework of multiple-scale expansion is in the
time scale $t_{2}$ and space scale $x_{1}$ (along with equation (21) in the
time scale $t_{1}$ and space scale $x_{1}$), while the GLE (26) obtained using
the amplitude expansion is in the time scale $t$ and space scale $x$.
Therefore, to make their coefficients comparable, we combine the GLE (22) with
equation (21) and transform them into

$\frac{\partial A}{\partial t}=-c_{g}\frac{\partial A}{\partial
x}+a_{1}\epsilon^{2}A+a_{2}\frac{\partial^{2}A}{\partial
x^{2}}+a_{3}{|A|}^{2}A,$ (30)

where the perturbation amplitude is also transformed from $A_{1}$ to $A$ using
the relationship $A=\epsilon A_{1}$. To make clear the physical meaning of the
perturbation amplitude $A$ for the Landau coefficient $a_{3}$ to be uniquely
determined, the linear eigenvectors ($\tilde{\bm{\gamma}}_{1}$ in the
multiple-scale expansion and $\tilde{\bm{\gamma}}_{11}$ in the amplitude
expansion) in the two expansion methods are normalized so that $\mathcal{A}$
denotes the global perturbation amplitude of total energy departing from the
base state (Thompson  _et al._ , 2001), which is defined by

$|\mathcal{A}|^{2}=\frac{1}{U_{\infty}^{2}\mathcal{V}}\iint_{\mathcal{V}}\frac{1}{2}\left[\mathbf{u}^{2}+M^{2}\mathbf{e}^{2}\right]dxdy,\
\ \ \text{with}\ \mathbf{u}=\mathbf{U}-\bar{\mathbf{U}},\
\mathbf{e}=\mathbf{E}-\bar{\mathbf{E}},$ (31)

where the velocity $U_{\infty}$ and system volume $\mathcal{V}$ are applied
for the purpose of nondimensionalization. In the current work, the maximum
velocity magnitude of the Poiseuille flow is used for defining $U_{\infty}$.

Then the group velocity $c_{g}$ should be compared with $c_{gA}$.
$a_{1}\epsilon^{2}$ (or the complex frequency $\mu$ in equation (18a) away
from the linear critical condition) in the multiple-scale expansion should be
compared with $\sigma$ in the amplitude expansion, both denoting the growth
rate of perturbed amplitude $\mathcal{A}$. The coefficients $a_{2}$ and
$a_{2A}$ share the same meaning. $a_{3}$ should be compared with $K$, both
referring to the Landau coefficient and are calculated within the weakly
nonlinear regime (i.e. not too far from the linear instability criterion
$T_{c}$). If they are positive, then the flow transition is subcritical;
otherwise, the flow transition is supercritical. We will compare the results
of these two methods in the result section.

In the end, we would like to mention in passing that Soldati & Banerjee in
Ref. Soldati and Banerjee (1998) have used a triple-decomposition method
Reynolds and Hussain (1972) to study the turbulent flows in EHD, with the flow
being decomposed into a sum of a time-average component, a coherent component
and the remaining turbulent component. The triple-decomposition method can be
applied to fully turbulent flows (with an emphasis on revealing the flow
dynamics of the coherent structures), whereas, the expansion methods in this
work are traditionally applied to transitional flows transitioning from the
laminar flows.

####  II.4.4 Spatiotemporal instability analysis based on GLE

The linear convective and absolute instabilities of the EHD-Poiseuille flow
can also be studied using the GLE. Due to its simplicity, this equation can be
considered as a reduced-order model to predict the absolute instability
criterion. Linearising the GLE (30) and the GLE (26) results in respectively

$\displaystyle\frac{\partial A}{\partial t}$
$\displaystyle=-c_{g}\frac{\partial A}{\partial
x}+a_{1}\epsilon^{2}A+a_{2}\frac{\partial^{2}A}{\partial
x^{2}}\,\,\,\,\text{(based on the multiple-scale expansion)},$ (32a)
$\displaystyle\frac{\partial A}{\partial t}$
$\displaystyle=-c_{gA}\frac{\partial A}{\partial x}+\sigma
A+a_{2A}\frac{\partial^{2}A}{\partial x^{2}}\,\,\,\,\text{(based on the
amplitude expansion)},$ (32b)

where the coefficients $c_{g}$, $a_{1}$ and $a_{2}$ should be evaluated near
the linear critical condition $T_{c}$, whereas the coefficients $c_{gA}$,
$\sigma$ and $a_{2A}$ can be evaluated at parameters away (but still not too
far) from the linear critical condition $T_{c}$. However, the critical
electric Rayleigh number $T_{ca}$ beyond which the linear instability changes
from convective to absolute may be not close to $T_{c}$, and thus, the
predicted value of these coefficients (and the resulting $T_{ca}$) from the
above linearised GLE may be inaccurate (see section III.3.1 for the discussion
of this point and the comparison in table 2). Instead, we use the following
form of the linearised GLE and the corresponding dispersion relation
$D(\omega,\beta)$ can be obtained by assuming the normal mode solution
$\tilde{A}e^{i\beta x-i\omega t}$ (Gao  _et al._ , 2013; Suslov and Paolucci,
2004):

$\frac{\partial A}{\partial t}=-c_{gT}\frac{\partial A}{\partial
x}+a_{1T}A+a_{2T}\frac{\partial^{2}A}{\partial
x^{2}},\,\,\,\,D(\omega,\beta)\equiv-\omega+c_{gT}\beta+ia_{1T}-ia_{2T}\beta^{2}=0.$
(33)

Here, the coefficients $c_{gT}$, $a_{1T}$ and $a_{2T}$ are directly evaluated
from the dispersion relation (Stewartson and Stuart, 1971; Suslov and
Paolucci, 2004) of the linear equation (17a) at each $T$ of interest, i.e.,

$c_{gT}=i\frac{\partial\mu}{\partial\alpha},\,\,\,a_{1T}=\mu,\,\,\,a_{2T}=-\frac{1}{2}\frac{\partial^{2}\mu}{\partial\alpha^{2}}.$
(34)

$\beta$ is the modulation wavenumber of the fundamental wave with wavenumber
$\alpha_{m}$ at which the maximum linear growth rate is attained. The
dispersion relation $D(\omega,\beta)$ has an equilibrium at $\beta_{0}$ given
by the condition $\partial\omega/\partial\beta|_{\beta_{0}}=0$, which results
in $\beta_{0}=c_{gT}/(2ia_{2T})$. Then the absolute complex frequency
$\omega_{0}$ and absolute growth rate $\omega_{0,i}$ (imaginary part of
$\omega_{0}$) can be expressed as

$\omega_{0}\equiv\omega(\beta_{0})=\frac{c_{gT}^{2}}{2ia_{2T}}+ia_{1T}+\frac{ic_{gT}^{2}}{4a_{2T}},\,\,\,\,\omega_{0,i}=a_{1T,r}-\frac{c_{gT}^{2}a_{2T,r}}{4|a_{2T}|^{2}}.$
(35)

In the framework of linear stability analysis, for $T>T_{c}$ a negative
$\omega_{0,i}$ corresponds to convective instability, while absolute
instability happens when $\omega_{0,i}$ is positive. It should be mentioned
that the linear AICI is conventionally determined by locating the saddle point
in a complex wavenumber plane (Li  _et al._ , 2019); in the present analysis
$(\beta_{0}+\alpha_{m})$ directly corresponds to such a saddle point.

###  II.5 Numerical method

The highly-accurate numerical simulations of EHD-Poiseuille flows are
performed using the open-source computational fluid dynamics solver Nek5000 by
Fischer et al. Fischer  _et al._ (2008), based on the spectral element method
(SEM) proposed by Patera Patera (1984). As the current work focuses on the
spatiotemporal instability of the EHD-Poiseuille flow, above governing
equations (1) and (5) are numerically solved without any turbulence model. The
SEM method combines the accuracy of spectral method with the geometrical
flexibility of finite element method. Within each spectral element, the
orthogonal basis of Legendre polynomials is used on Gauss-Lobatto-Legendre
(GLL) points, which further refine each element and approximate the solution
of governing equations as an expansion of Lagrange interpolants. The current
spatial discretisation follows the $\mathbb{P}_{N}-\mathbb{P}_{N-2}$
formulation (Maday and Patera, 1989) with polynomial order $N=7$. For temporal
integration, we use the second-order backward difference scheme implicit for
linear terms, coupled to a second-order extrapolation scheme explicit for
nonlinear terms. The varying step time is used with Courant number being
$0.2$. For the weakly nonlinear stability analysis, we use the spectral
collocation method to solve the resultant eigenvalue problems and impose the
solvability condition, following Zhang Zhang (2016).

##  III Result and Discussions

###  III.1 Grid independence and validation

In this section, the grid independence of current EHD-Poiseuille flows is
first tested by comparing the linear instability criteria $T_{c}$ and
spatiotemporal instability criteria $T_{ca}$ under different element
resolutions with parameters $C=10$, $M=10$, $Fe=10^{4}$, $\mathcal{U}=1$ and
domain range $x\in[0,100]$, $y\in[-1,1]$. At $T>T_{c}$, the electroconvection
can take place, and it will transition from convectively to absolutely
unstable state when $T>T_{ca}$. Their results are listed in the Table 1. Based
on the polynomial order $N=7$ in the current work, each spectral element
further consists of $8\times 8$ GLL nodes. For element size $640\times 24$,
the total number of grid points should be $4481\times 169=757289$ (adjacent
elements share the same edge point). From Table 1, it can be seen that the
differences of instability criteria between element sizes $640\times 24$ and
$800\times 28$ are both smaller than $1\%$. Thus, in the following numerical
simulations, the spectral element size $640\times 24$ is used for streamwise
length $L=100$ unless otherwise stated. If other values of $L$ are used, the
number of spectral elements along the streamwise direction will be
proportionally changed with $L$.

Element Size |  $480\times 20$ |  $640\times 24$ |  $800\times 28$
---|---|---|---
$T_{c}$ |  $150.16$ |  $150.77$ |  $150.81$
$T_{ca}$ |  $265.35$ |  $265.47$ |  $265.48$
Table 1: Grid independence test for EHD-Poiseuille flows by comparing the
linear instability criterion $T_{c}$ and spatiotemporal instability criterion
$T_{ca}$ at parameters $C=10$, $M=10$, $Fe=10^{4}$ and $\mathcal{U}=1$ within
the domain range $[0,100]$ in $x$ direction and $[-1,1]$ in $y$ direction.
Each spectral element is further refined by GLL points in Nek5000.

(a)$t$$\mathcal{A}$linear
regime(b)$\alpha$$\omega$(c)$v_{g}$$\sigma$$v_{-}$$v_{+}$$v_{g}=0$(d)$y$$\tilde{\mathcal{E}}_{f}$

Figure 2: Temporal and spatiotemporal instabilities of EHD-Poiseuille flows at
$C=10$, $M=10$, $Fe=10^{4}$, $\mathcal{U}=1$ and $T=190$. (a) temporal
evolution of global energy amplitude $\mathcal{A}$, the dotted line
corresponds to the linearly fitted curve, and two dashed lines indicate the
linear regime $t\in[6,20]$; (b) relation between wavenumber $\alpha$ and their
leading growth rates $\omega$; the solid line refers to current results, and
circles come from the linear stability analysis (Zhang  _et al._ , 2015), and
the dashed curve is the result obtained from the parabolic approximation of
the dispersion relation as in equation (33); (c) relation between the group
velocity $v_{g}=(x-x_{0})/t$ and $\omega$; (d) invariant eigenvectors
$\tilde{\mathcal{E}}_{f}$ at $\alpha=2.513$ (corresponds to the maximum growth
rate in panel b) at two different instants.

The accuracy of numerical simulations for electroconvective flows within the
framework of Nek5000 has been verified in our former work Feng  _et al._
(2021). Here, the capacity of the flow solver to study the temporal and
spatiotemporal instability is tested by analysing the snapshots of the EHD-
Poiseuille flows. To do so, we set the spatial domain range to $x\in[0,100]$
and $y\in[-1,1]$, and use the control parameters $C=10$, $M=10$, $Fe=10^{4}$,
$\mathcal{U}=1$ and $T=190$. Nonlinear governing equations (1) are evolved to
obtain snapshots of flow fields within the linear regime (about $t\in[6,20]$)
as shown in Fig. 2(a). The results of temporal instability analysis are shown
in panel b, where the leading growth rates (i.e. solid line) at different
wavenumbers are calculated from the data of numerical simulations and they
agree well with those from the local linear stability analysis (Zhang  _et
al._ , 2015) (circles). The maximum growth rate is $\omega=0.2409$ at
wavenumber $\alpha=2.513$. For the spatiotemporal instability analysis in
panel c, the maximum growth rate $\omega=0.2416$ is obtained at group velocity
$v_{g}=0.800$, which is consistent with above growth rate $\omega=0.2409$ with
a relative difference of $0.290\%$. At $v_{g}=0$, the growth rate is negative
$\sigma=-0.276$ which indicates a linearly convective type of instability for
the EHD-Poiseuille flow at $T=190$. The leading velocity $v_{+}=1.370$ and the
trailing velocity $v_{-}=0.290$ are also shown in panel c. They have the same
sign (i.e. positive) and the ray $v_{g}=0$ lies outside the range of
$[v_{-},v_{+}]$. Finally, the eigenvector $\mathcal{E}_{f}$ of wavenumber
$\alpha=2.513$ at $t=8$ and $t=10$ coincide with each other as shown in panel
d, which means that they are proportional to each other (as they are at the
asymptotically large time) and the current calculated growth rates are
faithful.

$T$$\sigma$

Figure 3: Variation of absolute growth rate with the electric Rayleigh number
for EHD-Poiseuille flows at $C=50$, $M=50$, $Fe=10^{4}$ and $\mathcal{U}=1$.
The circular symbols refer to current results, and solid line comes from the
work of Li  _et al._ (2019). $\sigma<0$: linearly convective instability;
$\sigma>0$: linearly absolute instability; $\sigma=0$: linearly neutral,
$T_{ca}=200.89$.

In addition, the absolute growth rates $\sigma(v_{g}=0)$ are computed for
discriminating the transition from linearly convective instability to absolute
instability, and they are compared with the linear results from Li et al. Li
_et al._ (2019). The parameters are $C=50$, $M=50$, $Fe=10^{4}$ and
$\mathcal{U}=1$. As the injection strength increases to $C=50$, the mesh
resolution in vertical direction is further refined to have $36$ spectral
elements. The current results are shown in Fig. 3. It can be found that the
absolute growth rates and criterion of absolute instability (i.e.
$T_{ca}=200.89$ in the current work) agree well with the linear results in
Ref. (Li  _et al._ , 2019) (whose $T_{ca}=200.37$), which verifies the
accuracy of the current algorithm for performing the spatiotemporal
instability analysis of the EHD-Poiseuille flows.

###  III.2 Weakly nonlinear stability analysis

In this section, the characteristics of transition from hydrostatic to
electroconvective state around the linear instability criterion $T_{c}$ are
first studied. In order to evaluate Landau coefficients from the data of
numerical simulations, we assume the following Stuart-Landau equation

$\frac{d\mathcal{A}}{dt}=c_{1}\mathcal{A}+c_{3}\mathcal{A}|\mathcal{A}|^{2},\,\,\,\frac{d\ln\mathcal{A}}{dt}=c_{1}+c_{3}|\mathcal{A}|^{2}.$
(36)

After obtaining the time series of global energy amplitude $\mathcal{A}(t)$ in
Eq. 31 from numerical simulations, $\frac{d\ln\mathcal{A}}{dt}$ can be plotted
as a function of $|\mathcal{A}|^{2}$. Then $c_{1}$ can be approximated by the
vertical intercept and $c_{3}$ by the slope. In terms of the coefficients
obtained from the weakly nonlinear analysis, $c_{1}$ should be compared with
$\mu$ and $\sigma$ and $c_{3}$ should be compared with $a_{3}$ and $K$.

(a)$\mathcal{A}$$t$$c_{1}=0.00717$(b)$\frac{d\ln\mathcal{A}}{dt}$$|\mathcal{A}|^{2}$$c_{3}=1.57320$

Figure 4: The first transition of EHD-Poiseuille flows characterized by the
global energy amplitude $\mathcal{A}$ at $C=10$, $M=10$, $Fe=10^{4}$,
$\mathcal{U}=1$ and $T=152$. (a) temporal evolution of $\mathcal{A}$; (b)
relation between $d(\ln\mathcal{A})/dt$ and $|\mathcal{A}|^{2}$.

With the parameters $C=10$, $M=10$, $Fe=10^{4}$, $\mathcal{U}=1$ and $T=152$
in a computational domain $x\in[0,2.456]$ and $y\in[-1,1]$, the results of the
growth rate $c_{1}$ and the Landau coefficient $c_{3}$ can be computed as
$c_{1}=0.00717$ and $c_{3}=1.57320$, as shown in Fig. 4, by performing the
linear-least-square fitting of the data in panels a and b (we will compare
these results with those in the weakly nonlinear stability theory to be
followed). Around $|\mathcal{A}|^{2}=0$ in panel b, the perturbation amplitude
is weak and the system stays within the linear regime. Thus, the $y$-intercept
corresponds the growth rate $c_{1}$ obtained in panel a. In the current case,
the Landau coefficient $c_{3}$ is found to be positve at $T=152$ close to the
linear instability criterion $T_{c}=150.77$, which indicates that the
transition from hydrostatic to electroconvective state is subcritical in the
EHD-Poiseuille flows Zhang (2016); Guan and Novosselov (2019) at the
aforementioned parameters. From the viewpoint of temporal evolution of
$\mathcal{A}$ in panel a, the deviation from the linear growth after $t=580$
is associated with an increased growth rate, which indicates that the higher-
order terms destabilises the flow (around the linear critical condition),
consistent with the subcritical bifurcation.

(a)$c_{1}$$T$(b)$c_{3}$$T$

Figure 5: Comparison of the linear growth rate and Landau coefficient between
the data of numerical simulations (circular symbols) and those from the
multiple-scale method (solid blue line) and amplitude expansion method (red
dashed line) at $C=10$, $M=10$, $Fe=10^{4}$, $\mathcal{U}=1$. (a) growth rate
$c_{1}$; (b) Landau coefficient $c_{3}$.

(a)$c_{1}$$M$(b)$c_{3}$$M$

Figure 6: Comparison of the linear growth rate and Landau coefficient between
current results (circular symbols) and those from the multiple-scale method
(solid blue line) and amplitude expansion method (red dashed line) at $C=10$,
$T=152$, $Fe=10^{4}$, $\mathcal{U}=1$. (a) growth rate $c_{1}$; (b) Landau
coefficient $c_{3}$.

In the above paragraph, the linear growth rate and the Landau coefficient are
approximately obtained from the time series of global energy amplitude in
numerical simulations. On the other hand, starting from the governing
equations (1) of EHD-Poiseuille flows, the GLE can be explicitly derived by
perturbation methods. We will consider both the multiple-scale method and
amplitude expansion method and compare their results with the linear growth
rate and the Landau coefficient obtained from the data of numerical
simulations. As shown in Fig. 5, the relation between $c_{1}$, $c_{3}$ and $T$
is investigated with other parameters being invariant. The linear growth rates
$c_{1}$ obtained from three methods are visually the same to each other. For
Landau coefficient $c_{3}$ in panel b, around the linear instability criterion
$T_{c}=150.77$, $c_{3}$ predicted from three methods agree well with each
other, which indicates the equivalence of the two perturbation methods around
the linear instability criterion. However, as the increase of electric
Rayleigh number from $T_{c}$, the value of $c_{3}$ obtained from the data of
numerical simulations gets closer to that of the amplitude expansion method.
This clearly confirms that the validity of amplitude expansion method is wider
than multiple-scale method when the governing parameter $T$ is away from the
linear critical value $T_{c}$. Except for the effect of electric Rayleigh
number on $c_{1}$ and $c_{3}$, the mobility ratio $M$ is also changed to
display the dependence of $c_{1}$ and $c_{3}$ on it, as shown in Fig. 6. The
electric Rayleigh number is fixed at $T=152$ which is very close to the linear
instability criterion $T_{c}=150.77$. Around this $T$, the multiple-scale
method and amplitude expansion method are approximately equivalent. Thus, the
results of $c_{1}$ and $c_{3}$ obtained from three methods are consistent with
each other. The reason may be that changing $M$ will not change the linear
critical condition Zhang  _et al._ (2015) and will not drive the flow away
from the linear critical condition (note that the $T$ used is close to
$T_{c}$).

$Ne$$T$$T_{c}$$T_{f}$

Figure 7: Hysteresis loops of EHD-Poiseuille flows at $C=10$, $M=10$,
$Fe=10^{4}$ and different $\mathcal{U}$. (1) Black solid line:
$\mathcal{U}=0$, $T_{c}=157.27$ and $T_{f}=107.4$; (2) Green solid line:
$\mathcal{U}=0.2$, $T_{c}=156.98$ and $T_{f}=107.8$; (3) Blue solid line:
$\mathcal{U}=0.5$, $T_{c}=155.48$ and $T_{f}=109.6$; (4) Red solid line:
$\mathcal{U}=1$, $T_{c}=150.77$ and $T_{f}=115.2$.

One salient feature of the subcritical transition is the existence of a
hysteresis loop in the bifurcation diagram, which has been studied and
discussed in the electroconvection without through-flows (e.g., in our
previous work Feng  _et al._ (2021)). When exposed to the Poiseuille flow, the
linear instability criterion $T_{c}$ and finite-amplitude criterion $T_{f}$
forming the hysteresis loop can present different results compared with the
case without through-flows. With parameters $C=10$, $M=10$ and $Fe=10^{4}$,
the hysteresis loops at $\mathcal{U}\in[0,0.2,0.5,1]$ are computed as shown in
Fig. 7. The electric Nusselt number $Ne$ within the range $T\in[T_{f},T_{c}]$
is obtained by gradually decreasing $T$. When the through-flow is intensified,
it can be seen that the electric transferring efficiency is reduced (i.e. $Ne$
becomes smaller) as the vertical motions are weakened by the through-flow Guan
and Novosselov (2019). The results of linear instability criteria $T_{c}$
agree well with former results (Zhang  _et al._ , 2015) and $T_{f}$ increases.
Furthermore, as the increase of $\mathcal{U}$ from $0$ to $1$, the range
between $T_{f}$ and $T_{c}$ decreases, which indicates that the subcritical
nature of the pure electroconvective flows can be suppressed by the Poiseuille
flow.

(a)$y$(b)$y$(c)$y$$x$

Figure 8: Distribution of the initially localised impulse defined by equations
(37) in the EHD-Poiseuille flow. Only part of the flow field is presented here
for a clear visualisation as the impulse is located around the central
position $x_{0}=30$ and $y_{0}=0$. (a) streamwise velocity; (b) normal
velocity; (c) charge density.

###  III.3 Spatiotemporal instability analysis

Next, the convective and absolute instabilities of EHD-Poiseuille flows are
analyzed by studying both linear and nonlinear responses of the flow field to
an initially localized impulse. For the profile of the impulse, we follow that
in Ref. (Delbende and Chomaz, 1998)

$\displaystyle u(x,y,t=0)$
$\displaystyle=-A_{m}(y-y_{0})\exp\left[-\left(\frac{(x-x_{0})^{2}}{2\sigma_{x}^{2}}+\frac{(y-y_{0})^{2}}{2\sigma_{y}^{2}}\right)\right],$
(37a) $\displaystyle v(x,y,t=0)$
$\displaystyle=A_{m}(x-x_{0})\frac{\sigma_{y}^{2}}{\sigma_{x}^{2}}\exp\left[-\left(\frac{(x-x_{0})^{2}}{2\sigma_{x}^{2}}+\frac{(y-y_{0})^{2}}{2\sigma_{y}^{2}}\right)\right],$
(37b) $\displaystyle q(x,y,t=0)$
$\displaystyle=A_{m}\exp\left[-\left(\frac{(x-x_{0})^{2}}{2\sigma_{x}^{2}}+\frac{(y-y_{0})^{2}}{2\sigma_{y}^{2}}\right)\right],$
(37c)

where $(x_{0},y_{0})$ refers to the central location of the impulse, and the
extent $\sigma_{x}=0.3$ along $x$-axis and $\sigma_{y}=0.2$ along $y$-axis.
The initial perturbation amplitude is set to be $A_{m}=10^{-3}$. The envelope
of this impulsed disturbance takes the Gaussian form, and two initial velocity
components $u$ and $v$ satisfy the requirement of continuity equations (5a).
This localized impulse initially added into the EHD-Poiseuille flow is shown
in Fig. 8 for reference.

(a)$v_{g}$$t^{\frac{1}{2}}\hat{\tilde{\mathcal{E}}}$$T=255$(b)$v_{g}$$t^{\frac{1}{2}}\hat{\tilde{\mathcal{E}}}$$T=275$(c)$v_{g}$$\sigma$$T=255$$v_{-}=0.028$$v_{+}=1.708$(d)$v_{g}$$\sigma$$T=275$$v_{-}=-0.025$$v_{+}=1.775$

Figure 9: Variation of $t^{\frac{1}{2}}\hat{\tilde{\mathcal{E}}}$ and $\sigma$
with the group velocity $v_{g}$ for EHD-Poiseuille flows at $C=10$, $M=10$,
$Fe=10^{4}$ and $\mathcal{U}=1$ for two values of $T$. (a-b) Temporal
evolution of $t^{\frac{1}{2}}\hat{\tilde{\mathcal{E}}}$; (c-d) spatiotemporal
growth rate $\sigma$ as a function of $v_{g}$.

After the onset of electroconvective flows (i.e. $T>T_{c}$), under the effects
of the electric field and the through-flow, the initial impulse in the flow
field can develop both temporally and spatially. To study the linear and
nonlinear responses of such impulse under periodic streamwise boundaries and
to avoid contamination, the spatial range along the streamwise direction is
extended to $x\in[0,100]$ to prevent the re-entrance of this impulse from the
outlet to inlet of this spatial domain. The parameters $C=10$, $M=10$ and
$Fe=10^{4}$ are kept constant in this section.

####  III.3.1 Linear impulse response

The transition from linearly convective instability to linearly absolute
instability can be discriminated by the absolute growth rate $\sigma$ of the
ray $v_{g}=0$, which varies almost linearly with the driving parameter $T$
(see Fig. 3). As defined in equation (10), the spatiotemporal growth rate
$\sigma$ can be determined by the exponentially temporal variation of
amplitude $t^{\frac{1}{2}}\hat{\tilde{\mathcal{E}}}$ at a given group velocity
$v_{g}$. Their relations are presented in Fig. 9(a-b) for two values of $T$.
It can be seen that the crossings of these amplitude curves define two group
velocities (i.e. $v_{-}$ and $v_{+}$) where the amplitude
$t^{\frac{1}{2}}\hat{\tilde{\mathcal{E}}}$ neither grows nor decays. Within
these two limiting edges $v_{g}\in[v_{-},v_{+}]$, the wavepacket always grows.
Based on the equation (11), the spatiotemporal growth rates can be calculated
at different group velocities as shown in panels c and d. At $T=255$, the
amplitude curve at $v_{g}=0$ in panel a decays, which indicates that
$\sigma(v_{g}=0)$ is negative and the flow is convectively unstable (or
$v_{-}$ and $v_{+}$ share the same direction). As the electric Rayleigh number
is increased to $T=275$, the amplitude curve in panel b grows at $v_{g}=0$ and
$\sigma(v_{g}=0)>0$, which indicates that the flow becomes absolutely unstable
(or the signs of $v_{-}$ and $v_{+}$ are different).

Results from | $T$ | $150.95\,\ (T_{c})$ | $190$ | $250$ | $300$
---|---|---|---|---|---
| $c_{g}$ | $0.7813$ | $0.8029$ | $0.8214$ | $0.8303$
multiple-scale | $a_{1,r}(T-T_{c})$ | $0$ | $0.2139$ | $0.4096$ | $0.5000$
expansion | $a_{2}$ | $0.2613-0.0177i$ | $0.2614-0.0252i$ | $0.2622-0.0275i$ | $0.2614-0.0274i$
| $\omega_{0M,i}$ | $-0.5814$ | $-0.3969$ | $-0.2267$ | $-0.1523$
| $c_{gA}$ | $0.7813$ | $0.8029$ | $0.8214$ | $0.8303$
amplitude | $\sigma_{r}$ | $0.0000$ | $0.2374$ | $0.5230$ | $0.7093$
expansion | $a_{2A}$ | $0.2613-0.0177i$ | $0.2614-0.0252i$ | $0.2622-0.0275i$ | $0.2613-0.0274i$
| $\omega_{0A,i}$ | $-0.5814$ | $-0.3734$ | $-0.1133$ | $0.0570$
| $c_{gT}$ | $0.7812$ | $0.8028$ | $0.8214$ | $0.8302$
dispersion | $a_{1T,r}$ | $-0.0000$ | $0.2374$ | $0.5230$ | $0.7092$
relation | $a_{2T}$ | $0.2615-0.0176i$ | $0.2616-0.0251i$ | $0.2624-0.0275i$ | $0.2614-0.0273i$
| $\omega_{0,i}$ | $-0.5809$ | $-0.3729$ | $-0.1129$ | $0.0573$
Relative | $|(\omega_{0M,i}-\omega_{0,i})/\omega_{0,i}|$ | $0.09\%$ | $6.44\%$ | $100.80\%$ | $365.79\%$
errors | $|(\omega_{0A,i}-\omega_{0,i})/\omega_{0,i}|$ | $0.09\%$ | $0.13\%$ | $0.35\%$ | $0.52\%$
Table 2: Coefficients in the linearised GLE for EHD-Poiseuille flows at
different $T$ and parameters $C=10$, $M=10$, $Fe=10^{4}$ and $\mathcal{U}=1$.
The coefficients $c_{g}$, $a_{1,r}$ and $a_{2}$ are calculated based on the
multiple-scale expansion using equations (21) and (23); the corresponding
absolute growth rate is
$\omega_{0M,i}=a_{1,r}(T-T_{c})-{c_{g}^{2}a_{2,r}}/({4|a_{2}|^{2}})$. The
coefficients $c_{gA}$, $\sigma_{r}$ and $a_{2A}$ are calculated based on the
amplitude expansion using equations (27a) and (29); the corresponding absolute
growth rate is
$\omega_{0A,i}=\sigma_{r}-{c_{gA}^{2}a_{2A,r}}/({4|a_{2A}|^{2}})$. The
coefficients $c_{gT}$, $a_{1T,r}$ and $a_{2T}$ are directly evaluated from the
dispersion relation using equation (34), and the corresponding absolute growth
rate is calculated using equation (35).The last two rows show the relative
errors between the absolute growth rate directly obtained from the dispersion
relation and those obtained using the two expansion methods.

As mentioned in section II.4.4, the absolute growth rate $\omega_{0,i}$ in
equation (35), which shares the same meaning as $\sigma(v_{g}=0)$ defined in
equation (11), can also be used to determine the AI/CI property (here we use
different notations for the absolute growth rate to differentiate it in
different methods, i.e., $\omega_{0,i}$ in the weakly nonlinear analysis and
$\sigma(v_{g}=0)$ in the post-processing analysis of snapshots from numerical
simulations). To calculate $\omega_{0,i}$, the three coefficients in the
linearised GLE should be evaluated first. As we have explained above, there
are three methods for doing this: the first is to calculate the coefficients
$c_{g}$, $a_{1}$ and $a_{2}$ in equation (32a) based on the formulae (21) and
(23); the second is to calculate $c_{gA}$, $\sigma$ and $a_{2A}$ in equation
(32b) according to the formulae (27a) and (29); and the third is to directly
use the dispersion relation given in equation (34) which corresponds to the
linear eigenvalue problem (17a). In order to quantify the differences between
these three methods, Table 2 compares the coefficients and the corresponding
absolute growth rates. At the linear critical condition $T_{c}=150.95$ (very
close to the value 150.77 obtained from numerical simulations in Fig. 7), a
favourable agreement can be observed and the relative errors between the
absolute growth rates are just 0.09%. However, in terms of the multiple-scale
expansion, with $T$ increased from $T_{c}$, the error increases rapidly, which
arises mainly from the discrepancy between the linear growth rate
$a_{1,r}(T-T_{c})$ and $a_{1T,r}$; as for the amplitude expansion, the error
increases much more slowly. Even at $T=300$, far away from the linear critical
condition, the calculated absolute growth rate in the amplitude expansion
method has only a relative error $0.52\%$ compared to the dispersion relation
method. This illustrates that the coefficients evaluated based on the
multiple-scale expansion indeed cannot be used for a reliable analysis when
$T$ is far from its critical value, while the results of the amplitude
expansion method remain quite accurate even far away from the linear critical
condition. To the best of our knowledge, this appears to be a new result
demonstrating the difference between the two expansion methods in the weakly
nonlinear phase in terms of the absolute instability.

In order to gain a global picture of how the absolute growth rate varies with
$T$, in Fig. 10, we plot the absolute growth rates $\sigma(v_{g}=0)$ as a
function of $T$ for $\mathcal{U}=[0.2,0.5,1]$ (red circles). Plotted in the
same figures are the absolute growth rate computed using the GLE (blue
squares). The critical electric Rayleigh number $T_{ca}$ indicating the onset
of linearly absolute instability is represented by solid symbols (when
$T<T_{ca}$, convectively unstable and when $T>T_{ca}$, absolutely unstable),
and they are computed by the linear least squares method. As the increase of
through-flow $\mathcal{U}$, it can be clearly seen that the spatiotemporal
instability criterion $T_{ca}$ also increases, which is consistent with the
linear result in our previous work (Li  _et al._ , 2019).

(a)$\mathcal{U}=0.2$$\sigma$$T$(b)$\mathcal{U}=0.5$$T$(c)$\mathcal{U}=1$$T$

Figure 10: Variation of the absolute growth rate $\sigma(v_{g}=0)$ with the
electric Rayleigh number $T$ for EHD-Poiseuille flows at $C=10$, $M=10$ and
$Fe=10^{4}$. The red circular symbols refer to current results computed by
post-processing the snapshots from numerical simulations, blue square symbols
are results predicted from the GLE, and solid lines denote the corresponding
linearly fitted results with solid dots indicating $T_{ca}$. At $T<T_{ca}$,
$\sigma<0$, the flow is convectively unstable; at $T>T_{ca}$, $\sigma>0$, the
flow is absolutely unstable; at $T=T_{ca}$, $\sigma=0$, the flow is linearly
neutral. (a) $\mathcal{U}=0.2$, $T_{ca}^{Simulation}=160.63$ and
$T_{ca}^{GLE}=160.74$; (b) $\mathcal{U}=0.5$, $T_{ca}^{Simulation}=178.96$ and
$T_{ca}^{GLE}=181.63$; (c) $\mathcal{U}=1$, $T_{ca}^{Simulation}=265.47$ and
$T_{ca}^{GLE}=282.06$.

At $\mathcal{U}=0.2$, the spatiotemporal instability criteria $T_{ca}$
predicted from current post-processing method of numerical simulations and the
GLE model agree well with each other as observed from Fig. 10, and their
difference is about $0.07\%$. While as the intensification of through-flow,
the difference between above two results gradually increases to $1.49\%$ at
$\mathcal{U}=0.5$ and $6.25\%$ at $\mathcal{U}=1$. This discrepancy could be
attributed to the fact that the dispersion relation of equation (17a) can not
be approximated well enough by a quadratic function over the range of
streamwise wavenumbers within which the linear growth rate is positive (Suslov
and Paolucci, 2004). For example, Fig. 2(b) illustrates the difference between
the parabolic approximation (dashed curve) and the true dispersion relation
(solid curve). Therefore, the validity of the linearised GLE is restricted,
leading to a reduced predictive ability of the convective and absolute
instabilities (Suslov and Paolucci, 2004). In addition, as the GLE is normally
derived around the linear instability criterion $T_{c}$ (Zhang, 2016), it is
only valid in a small region near $T_{c}$, and their results may be misleading
if $T$ is far from $T_{c}$. From Figs. 7 and 10, it can be seen that the
deviation between linear instability criterion $T_{c}$ and spatiotemporal
instability criterion $T_{ca}$ becomes larger as the increase of through-flow
strength $\mathcal{U}$.

(a)$\mathcal{A}$$t$

transient

exponential

growth

saturation

wavepacket

spreading

(b)$\eta$$x$$t=1$$t=10$$t=19$$t=28$(c)$\eta$$x$exponential
growth(d)$x$saturation(e)$x$wavepacket spreading

Figure 11: Temporal evolution of the global energy amplitude $\mathcal{A}$ and
the enstrophy $\eta$ describing the envelope of wavepackets for EHD-Poiseuille
flows at $C=10$, $M=10$, $Fe=10^{4}$, $\mathcal{U}=1$ and $T=220$. (a)
temporal evolutions of $\mathcal{A}$, it consists of four stages: transient;
exponential growth of amplitude; saturation; wavepacket spreading; (b)
waterfall of the wavepacket denoted by normalized enstrophy, the dotted line
denotes the initial position $x_{0}=20$ of the impulse; (c) wavepacket within
the regime of exponential growth; (d) wavepacket within the regime of
saturation; (e) wavepacket within the regime of streamwise spreading.

####  III.3.2 Nonlinear impulse response

When the initial local impulse becomes finite-amplitude and finite-extent,
nonlinearity becomes important. The complete temporal evolution of the global
energy amplitude $\mathcal{A}$ is shown in Fig. 11(a). The envelopes of
wavepackets within linear and nonlinear regimes are represented by the
enstrophy $\eta$ defined in equation (12), whose spatiotemporal evolution is
described by the waterfall plot in panel b. The wavepacket at each instant is
normalized by its maximum amplitude for a clear visualization. In panel a,
after an initial transient period, the local impulse starts to grow
exponentially as indicated by the linear dotted line. Within this linear
range, the amplitude and extent of the impulse both grow with time as shown in
panel b. After this linear regime, the wavepacket undergoes a saturation
process, during which the growth rate of global amplitude $\mathcal{A}$
decreases. Such a phenomenon can also be clearly demonstrated in panel c,
where the increment of impulse amplitude after every $\Delta t=4$ becomes
smaller. Finally, in the regime of wavepacket spreading, the growth rate of
$\mathcal{A}$ becomes even smaller because the impulse amplitude stops to grow
and it only expands along the streamwise direction as shown in panel d. At
$T=220<T_{ca}=265.47$, the EHD-Poiseuille flow is convectively unstable. Thus,
the wavepacket overall moves to the downstream region as shown in these
panels.

$x$$y$

Figure 12: Distribution of charge density and streamline pattern of
perturbated velocity $\mathbf{u}$ for EHD-Poiseuille flows at $C=10$, $M=10$,
$Fe=10^{4}$, $\mathcal{U}=1$ and $T=220$. Only part of the flow field is
presented here for a clear visualisation. This snapshot corresponds to the
wavepacket at $t=34$ in above Fig. 11. Red arrows represent the upward motion
of charges, and high-charge regions form around them. Blue arrows indicate
downward moving charges with low density.

(a)$\mathcal{A}$$t$(b)$\eta$$x$$t=8$(c)$\eta$$x$$t=19$(d)$\eta$$x$$t=30$

Figure 13: Comparison of the global energy amplitude $\mathcal{A}$ and the
enstrophy $\eta$ between linear (blue solid lines) and nonlinear (red solid
lines) cases for EHD-Poiseuille flows at $C=10$, $M=10$, $Fe=10^{4}$,
$\mathcal{U}=1$ and $T=220$. The EHD-Poiseuille flow is convectively unstable.
(a) temporal evolutions of $\mathcal{A}$; (b) $t=8$ within regime of
exponential growth; (c) $t=19$ within regime of saturation; (d) $t=30$ within
regime of wavepacket spreading.

Within the regime of wavepacket spreading in panels b and e of Fig. 11, it can
be seen that those local maxima peaks are not equally spaced. At $T>T_{c}$,
the electroconvective flow takes place around the upward moving charges (they
form the high-charge region) or downward moving charges (they form the low-
charge region) as shown in Fig. 12. As the high-charge region is much thinner
than the low-charge one, the streamwise distance between centers of cell 1 and
2 ($\Delta x_{12}=1$) is smaller than that between cells 2 and 3 ($\Delta
x_{23}=1.7$). Thus, the enstrophy $\eta$ after saturation regime does not
present the uniformly distributed maxima peaks. From a point of view of flow
physics, it is because of the charge-void regions (see the almost white
regions in Fig. 12) and the flow field that causes the non-uniform peaks in
EHD-Poiseuille flows. This is different from the case in RBC subjected to a
Poiseuille flow, see the results in the Supplemental Material.

Next, we will compare the impulse response in linear and nonlinear regimes.
When starting from the same initial condition consisting of one local impulse,
the time series of global energy amplitude $\mathcal{A}$ obtained from
nonlinear equations (1) and linear equations (5) at $T=220$ are compared in
Fig. 13. Within the linear stage, the global amplitudes agree well between the
linear and nonlinear simulations. With the increase of the perturbation
amplitude, the deviation of the nonlinear simulation from the linear one
becomes larger. The envelopes of the wavepackets in linear and nonlinear
regimes at different instants are compared using enstrophy $\eta$ in panels b,
c and d. It makes sense that the wavepackets coincide with each other in the
linear stage (at $t=8$), as shown in panel a. At the subsequent instants, the
mismatch between linear and nonlinear results becomes more visible, but
interestingly, the leading and trailing edges in linear (moves at velocity
$v_{\pm}$) and nonlinear (moves at velocity $v_{\pm}^{NL}$) simulations are
still roughly the same after the period of linear growth, which means that the
front separating the saturated wavepacket from the basic flow moves at the
same speed in the linear and nonlinear regimes (i.e. $v_{-}=v_{-}^{NL}$ and
$v_{+}=v_{+}^{NL}$). At $T=300$, the EHD-Poiseuille flow becomes absolutely
unstable based on the conclusion from Fig. 10, and it is also found that the
wavepackets obtained from linear and nonlinear equations at the same instant
moves at the same speed (results now shown due to space limit). This
conclusion is also valid for $\mathcal{U}=0.2$ and $\mathcal{U}=0.5$ (not
shown due to space limit). These results demonstrate that the linearly
convective/absolute instability criteria can be used to predict the transition
from nonlinearly convective to absolute instability, at least for the EHD-
Poiseuille flow that we have studied here. Due to the similar boundary
conditions and flow patterns with EHD-Poiseuille flows, the results of
nonlinear impulse response in RBC-Poiseuille flows are also studied and
provided in the Supplemental Material. In the latter flow, the nonlinearity
also does not influence the edge velocities of wavepackets.

It is known that the bifurcation of EHD-Poiseuille flows is subcritical, and a
hysteresis loop exists in the bifurcation diagram (see Fig. 7). Thus, within
the parameter range $T\in[T_{f},T_{c}]$, an initially nonlinear finite-
amplitude and finite-extent impulse should evolve without decaying to the
steady state. Different from the linear impulse explicitly given by equation
(37), the ultimately saturated solution of EHD-Poiseuille flows at $T=160$
(larger than $T_{c}$) can be directly adopted as the initially nonlinear
impulse (i.e. initial condition) of EHD-Poiseuille flows within the finite-
amplitude regime. In order to understand the instability nature of these
finite-amplitude impulses, the spatiotemporal evolutions of nonlinear
wavepackets at $T=125$ and $T=145$ are shown in Fig. 14. At the same instant,
the wavepacket envelope of $T=145$ is wider than that of $T=125$ due to the
stronger electroconvection in former one. As the increase of time, the
variation of impulse amplitude is very small, while the overall wavepacket
moves to the downstream region which clearly indicates the nonlinearly
convective instability of EHD-Poiseuille flows within the finite-amplitude
range $T\in[T_{f},T_{c}]$. In the end, the above linear and nonlinear results
are summarised in Fig. 15. The difference between linear and nonlinear impulse
responses stems from the characteristic of subcritical bifurcation in EHD-
Poiseuille flows.

(a)$\eta$$x$$T=125$(b)$\eta$$x$$T=145$

Figure 14: Temporal evolution of the enstrophy $\eta$ within the finite-
amplitude range $T\in[T_{f},T_{c}]$ for EHD-Poiseuille flows at $C=10$,
$M=10$, $Fe=10^{4}$ and $\mathcal{U}=1$. (a) $T=125$; (b) $T=145$. Figure 15:
Linear and nonlinear impulse responses of the EHD-Poiseuille flow at $C=10$,
$M=10$ and $Fe=10^{4}$. S: stable; LC: linearly convective; LA: linearly
absolute; NC: nonlinearly convective; NA: nonlinearly absolute. $T_{f}$:
finite-amplitude instability criterion; $T_{c}$: linear instability criterion;
$T_{ca}$: instability criterion between LC and LA. Coexistence of $T_{c}$ and
$T_{f}$ indicates the subcritical nature of bifurcations.

$\epsilon^{\prime}$$c_{f}$

Figure 16: Variation of the front-propagation velocity $c_{f}$ with the
parameter $\epsilon^{\prime}=(T-T_{c})/T_{c}$ for EHD-Poiseuille flows at
$C=10$, $M=10$, $Fe=10^{4}$ and $\mathcal{U}=1$. The circular symbols refer to
the results of numerical simulations, and solid line comes from the
theoretical prediction $c_{f}\approx\sqrt{\epsilon^{\prime}}$.

In Fig. 11, after the stage of linear growth, the initially localised impulse
finally evolves to the saturated state with finite amplitude and extent, and
the linear edges can transform into nonlinear fronts (van Saarloos and
Hohenberg, 1992), which indicates the formation of electroconvection in
downstream regions of the background through-flow. Such front separating the
stable and unstable regimes can propagate with a well-defined velocity, namely
the front-propagation velocity $c_{f}$ (Fineberg and Steinberg, 1987; Van
Saarloos, 2003). This velocity quantifies flow state of the wavepacket, and is
related to the linear edge velocities $v_{\pm}$ shown in Figs. 9(c-d). Their
difference is that the former is computed in the reference frame of the
background through-flow, and the latter are obtained in the laboratory frame.
In the current work, the background through-flow (Poiseuille flow) does not
change with time, and edges of wavepackets from linear to nonlinear regimes
always coincide with each other (see Fig. 13). Thus, the front-propagation
velocity can be computed using the leading and trailing velocities of linear
wavepackets, namely $c_{f}=(v_{+}-v_{-})/2$. The variation of $c_{f}$ with the
control parameter $\epsilon^{\prime}=(T-T_{c})/T_{c}$ (which is similar to
$\epsilon$ in the weakly nonlinear stability analysis in Eqs. 15) is shown in
Fig. 16, and $\epsilon^{\prime}$ also describes the deviation of $T$ from its
linear instability criterion ($T_{c}=150.77$ in Fig. 16). Such relation can
also be predicted from an Ginzburg-Landau equation, which indicates that the
$c_{f}\propto\xi\sqrt{\epsilon^{\prime}}$ (Dee and Langer, 1983) (with $\xi$
being a certain pre-factor). In different flows, the value of $\xi$ may vary
(Ahlers and Cannell, 1983). For the EHD-Poiseuille flow, near the linear
instability criterion in Fig. 16, it can be seen that the front-propagation
velocity approximately follow the relation of
$c_{f}\approx\sqrt{\epsilon^{\prime}}$, while the deviation increases at high
$\epsilon^{\prime}$.

##  IV Conclusion and Future Work

In this paper, we investigated the spatiotemporal instabilities and weakly
nonlinear properties of the EHD-Poiseuille flows by a post-processing method
of highly accurate numerical simulations. We aim to characterise the linear
and nonlinear AI/CI in EHD-Poiseuille flows, as well as to test the potential
of GLE in predicting the flow nature around its first bifurcation and the
transition from convective to absolute instability. As the EHD flow is
subjected to a through-flow, the effect of this through-flow strength is also
studied. Our new findings are summarised below.

In the weakly nonlinear analysis of EHD-Poiseuille flows, we use three methods
to compute the Landau coefficient $c_{3}$ which indicates the nature of the
first bifurcation. The first two methods are based on the amplitude equation,
derived from two different expansion methods, namely the multiple-scale method
and amplitude expansion method. The third method computes the Landau
coefficient from the results of numerical simulations. Near the linear
instability threshold, the Landau coefficients computed from these three
methods are all positive, confirming the subcritical bifurcation of EHD-
Poiseuille flows. As the increase of electric Rayleigh number $T$, the values
of the Landau coefficient obtained in the amplitude expansion method agree
well with those from the post-processing method, but the deviation of $c_{3}$
becomes larger between the multiple-scale method and the post-processing
method. This clearly indicates the limited convergence range of the multiple-
scale expansion method and the applicability of the amplitude expansion
method, even when $T$ is away from the linear instability threshold. Using the
highly accurate simulations, the subcritical bifurcation diagrams of EHD-
Poiseuille flows with different through-flow strengths $\mathcal{U}$ were also
computed with a diagnosis of the electric Nusselt number $Ne$. As the increase
of $\mathcal{U}$, the linear instability criterion $T_{c}$ decreases while the
finite-amplitude criterion $T_{f}$ increases, which indicates a narrower
hysteresis loop.

By studying the response of a localised impulse in the flow, we extracted
information on the spatiotemporal instabilities of EHD-Poiseuille flows from
the snapshots of numerical simulations. The GLE has also been derived to
predict the transition from CI to AI in order to test its predicability. When
the through-flow is weak, as the spatiotemporal transition takes place near
the linear critical condition $T_{c}$, the spatiotemporal instability
criterion $T_{ca}$ predicted by the GLE agrees well with that from the post-
processing method of numerical simulations. While as the increase of
$\mathcal{U}$, $T_{ca}$ becomes larger than $T_{c}$, and the applicability of
GLE becomes increasingly questionable. On the other hand, when $\mathcal{U}$
is fixed, by changing $T$, we found that the GLE coefficients calculated using
the amplitude expansion method can yield absolute growth rates very close to
those obtained from the dispersion relation, a conclusion which does not hold
for the multiple-scale expansion method, whose predicted absolute growth rates
are only close to those computed by the dispersion relation near the linear
critical conditions. This result demonstrates the difference between the two
expansion methods in the weakly nonlinear phase in terms of the absolute
growth rate.

Because of the nonlinearity in the flow, the localised impulse will, after the
linear growth phase, eventually evolve into a saturated wavepacket. It is
found that edges of the impulse and saturated wavepackets always coincide with
each other, which means that the nonlinearity does not influence the
propagation velocity of the wavepacket except for its amplitude (i.e., the
saturation). The same conclusion is also found for the RBC-Poiseuille flows in
the appendix. Within the finite-amplitude range $T\in[T_{f},T_{c}]$, the
spatiotemporal property of EHD-Poiseuille flows is proved to be nonlinearly
convective at least for the parameters we have investigated in this work. In
the nonlinear evolution process of the wavepacket, it is also found that the
front-propagation velocity increases with the control parameter
$\epsilon^{\prime}=(T-T_{c})/T_{c}$ according to a power law (i.e.
$c_{f}\propto\sqrt{\epsilon^{\prime}}$), which is consistent with the result
predicted by amplitude equations.

The current work studies the nonlinear spatiotemporal development of the
disturbance in EHD-Poiseuille flows on, e.g., its transition from nonlinear CI
to nonlinear AI and the power law for the front-propagation velocity. These
results can potentially provide some useful implications for studying the
fluid dynamics in industrial processes like ESP and EV damper devices. To this
end, the same flow configuration as that, for example, of ESP (wires between
two plates) should be used. Besides, here we mainly study the effect of
electric Rayleigh number on the nonlinear AI/CI, future works can consider a
more complete parameter study (including the injection strength $C$ and
mobility ratio $M$) and explore a larger parameter space in this flow.
Finally, since now we have demonstrated that the EHD flow can be influenced by
the Poiseuille flow, it may also be interesting to study how the EHD rolls can
change the Poiseuille flow, for example, delaying the turbulence transition in
the Poiseuille flow Soldati and Banerjee (1998) from the perspective of flow
stability.

###### Acknowledgements.

We acknowledge the financial support of a Start-up grant from the Ministry of
Education, Singapore to M.Z. (with the WBS No. R-265-000-689-114) and doctoral
research scholarships from the National University of Singapore to Z.F. and
D.W. D.W. is also partially supported by the China Scholarship Council. The
computational resources of the National Supercomputing Centre, Singapore
(https://www.nscc.sg) are acknowledged.

##  Appendix A Operators in the weakly nonlinear analysis

In this appendix, we provide the explicit expressions of the various operators
appearing in the weakly nonlinear analysis, including those in section II.4.1
for multiple-scale expansion and those in section II.4.2 for amplitude
expansion. All the below expressions are completely consistent with those in
our previous study (Zhang, 2016) (see the appendices therein); the only
difference is that in that paper the operators are given in physical space,
while in the present work we express them directly in the spectral space.

###  A.1 Weight matrices

The weight matrix $\bm{M}$ in physical space introduced in equation (14) is
$\bm{M}=\begin{bmatrix}\nabla^{2}&0\\\ 0&\nabla^{2}\end{bmatrix}$. To
facilitate the expression of various operators related to it, we introduce in
spectral space the operator
$\tilde{\bm{M}}^{(n)}=\begin{bmatrix}\nabla_{n}^{2}&0\\\
0&\nabla_{n}^{2}\end{bmatrix}$ where
$\nabla_{n}^{2}=-n^{2}\alpha^{2}+({\partial^{2}}/{\partial y^{2}})$ and
$\nabla_{n}^{4}=\nabla_{n}^{2}\nabla_{n}^{2}$ for later use. Then for the
multiple-scale expansion in section II.4.1, the $\bm{M}$–related operators
appearing in equation (21) and (23) can be expressed as
$\tilde{\bm{M}}_{0}=\tilde{\bm{M}}^{(1)}$,
$\tilde{\bm{M}}_{1}^{\circ}=\begin{bmatrix}2i\alpha&0\\\
0&2i\alpha\end{bmatrix}$ and $\tilde{\bm{M}}_{2}^{\circ}=\begin{bmatrix}1&0\\\
0&1\end{bmatrix}$. For the amplitude expansion in section II.4.2, the
$\bm{M}$–related operators appearing in equation (25) are
$\tilde{\bm{M}}_{0\alpha}=\tilde{\bm{M}}^{(0)}$,
$\tilde{\bm{M}}_{1\alpha}=\tilde{\bm{M}}^{(1)}$ and
$\tilde{\bm{M}}_{2\alpha}=\tilde{\bm{M}}^{(2)}$; the operators in the
expressions of $\tilde{\bm{n}}_{21}$ and $\tilde{\bm{n}}_{31z}$ in equation
(27) are $\tilde{\bm{M}}_{1\alpha}^{\circ}=\tilde{\bm{M}}_{1}^{\circ}$ and
$\tilde{\bm{M}}_{2\alpha}^{\circ}=\tilde{\bm{M}}_{2}^{\circ}$.

###  A.2 Linear operators

The linear operator $\bm{L}$ in equation (14) is in the physical space and it
is

$\bm{L}=\begin{bmatrix}-\bar{U}\frac{\partial\nabla^{2}}{\partial
x}+\bar{U}^{\prime\prime}\frac{\partial}{\partial
x}+\frac{M^{2}}{T}\nabla^{4}&-M^{2}\left(\bar{\phi}^{\prime}\frac{\partial\nabla^{2}}{\partial
x}-\bar{\phi}^{\prime\prime\prime}\frac{\partial}{\partial x}\right)\\\
\bar{\phi}^{\prime\prime\prime}\frac{\partial}{\partial
x}&\bar{\phi}^{\prime\prime\prime}\frac{\partial}{\partial
y}+\left(\bar{\phi}^{\prime}\frac{\partial}{\partial
y}-\bar{U}\frac{\partial}{\partial
x}\right)\nabla^{2}+2\bar{\phi}^{\prime\prime}\nabla^{2}+\frac{1}{Fe}\nabla^{4}\end{bmatrix}.$
(38)

Similar to $\tilde{\bm{M}}^{(n)}$, we introduce $\tilde{\bm{L}}^{(n)}$ for
convenience

$\tilde{\bm{L}}^{(n)}=\begin{bmatrix}-\bar{U}in\alpha{\nabla_{n}^{2}}+\bar{U}^{\prime\prime}in\alpha+\frac{M^{2}}{T}\nabla_{n}^{4}&-M^{2}\left(\bar{\phi}^{\prime}in\alpha{\nabla_{n}^{2}}-\bar{\phi}^{\prime\prime\prime}in\alpha\right)\\\
\bar{\phi}^{\prime\prime\prime}in\alpha&\bar{\phi}^{\prime\prime\prime}\frac{\partial}{\partial
y}+\left(\bar{\phi}^{\prime}\frac{\partial}{\partial
y}-\bar{U}in\alpha\right)\nabla_{n}^{2}+2\bar{\phi}^{\prime\prime}\nabla_{n}^{2}+\frac{1}{Fe}\nabla_{n}^{4}\end{bmatrix}.$
(39)

Then the $\bm{L}$–related operators in equations (21) and (23) for the
multiple-scale expansion read $\tilde{\bm{L}}_{0}=\tilde{\bm{L}}^{(1)}$, and

$\tilde{\bm{L}}_{1}^{\circ}=\begin{bmatrix}-\bar{U}(\nabla_{1}^{2}-2\alpha^{2})+\bar{U}^{\prime\prime}+\frac{4i\alpha
M^{2}}{T_{c}}\nabla_{1}^{2}&-M^{2}\left(\bar{\phi}^{\prime}\nabla_{1}^{2}-2\alpha^{2}\bar{\phi}^{\prime}-\bar{\phi}^{\prime\prime\prime}\right)\\\
\bar{\phi}^{\prime\prime\prime}&2i\alpha\left(\bar{\phi}^{\prime}\frac{\partial}{\partial
y}+2\bar{\phi}^{\prime\prime}\right)-\bar{U}(\nabla_{1}^{2}-2\alpha^{2})+\frac{4i\alpha}{Fe}\nabla_{1}^{2}\end{bmatrix},$
(40)
$\tilde{\bm{L}}_{2}^{\circ}=\begin{bmatrix}-3i\alpha\bar{U}+\frac{M^{2}}{T_{c}}(2\nabla_{1}^{2}-4\alpha^{2})&-3i\alpha
M^{2}\bar{\phi}^{\prime}\\\
0&\left(\bar{\phi}^{\prime}\frac{\partial}{\partial
y}+2\bar{\phi}^{\prime\prime}\right)-3i\alpha\bar{U}+\frac{1}{Fe}(2\nabla_{1}^{2}-4\alpha^{2})\end{bmatrix},\,\tilde{\bm{L}}_{2MT}=\begin{bmatrix}-\frac{M^{2}}{T_{c}^{2}}\nabla_{1}^{4}&0\\\
0&0\end{bmatrix}.$ (41)

Those appearing in equation (25) for the amplitude expansion are
$\tilde{\bm{L}}_{0\alpha}=\tilde{\bm{L}}^{(0)}$,
$\tilde{\bm{L}}_{1\alpha}=\tilde{\bm{L}}^{(1)}$ and
$\tilde{\bm{L}}_{2\alpha}=\tilde{\bm{L}}^{(2)}$. The operators in the
expressions of $\tilde{\bm{n}}_{21}$ and $\tilde{\bm{n}}_{31z}$ in equation
(27) are $\tilde{\bm{L}}_{1\alpha}^{\circ}=\tilde{\bm{L}}_{1}^{\circ}$ and
$\tilde{\bm{L}}_{2\alpha}^{\circ}=\tilde{\bm{L}}_{2}^{\circ}$.

###  A.3 Nonlinear operators

The nonlinear operator in physical space as introduced in equation (14) is
${\bm{N}}=(n_{\psi},n_{\varphi})^{T}$ with

$\displaystyle n_{\psi}=\left(\frac{\partial\psi}{\partial
x}\frac{\partial\nabla^{2}\psi}{\partial y}-\frac{\partial\psi}{\partial
y}\frac{\partial\nabla^{2}\psi}{\partial
x}\right)+M^{2}\left(\frac{\partial\varphi}{\partial
x}\frac{\partial\nabla^{2}\varphi}{\partial y}-\frac{\partial\varphi}{\partial
y}\frac{\partial\nabla^{2}\varphi}{\partial x}\right),$ (42a) $\displaystyle
n_{\varphi}=\left(\frac{\partial\varphi}{\partial
x}-\frac{\partial\psi}{\partial
y}\right)\frac{\partial\nabla^{2}\varphi}{\partial
x}+\left(\frac{\partial\varphi}{\partial y}+\frac{\partial\psi}{\partial
x}\right)\frac{\partial\nabla^{2}\varphi}{\partial
y}+\nabla^{2}\varphi\nabla^{2}\varphi.$ (42b)

The nonlinear operator $\tilde{\bm{N}}_{3}^{\circ}$ appearing in equation (23)
for the multiple-scale expansion is

$\tilde{\bm{N}}_{3}^{\circ}=\widetilde{\bm{N}}_{f}(\tilde{\bm{\gamma}}_{1},\tilde{\bm{\gamma}}_{20},1,0)+\widetilde{\bm{N}}_{f}(\tilde{\bm{\gamma}}_{20},\tilde{\bm{\gamma}}_{1},0,1)+\widetilde{\bm{N}}_{f}(\tilde{\bm{\gamma}}_{1}^{*},\tilde{\bm{\gamma}}_{22},-1,2)+\widetilde{\bm{N}}_{f}(\tilde{\bm{\gamma}}_{22},\tilde{\bm{\gamma}}_{1}^{*},2,-1),$
(43)

where the superscript ∗ denotes the complex conjugate of its argument and the
function $\widetilde{\bm{N}}_{f}$ is defined as
$\widetilde{\bm{N}}_{f}(\tilde{\bm{f}}_{1},\tilde{\bm{f}}_{2},q_{1},q_{2})=(\tilde{n}_{\psi},\tilde{n}_{\varphi})^{T}$
with

$\displaystyle\tilde{n}_{\psi}(\tilde{\bm{f}}_{1},\tilde{\bm{f}}_{2},q_{1},q_{2})=\left(iq_{1}\alpha\tilde{f}_{1\psi}\frac{\partial\nabla_{q_{2}}^{2}\tilde{f}_{2\psi}}{\partial
y}-\frac{\partial\tilde{f}_{1\psi}}{\partial
y}iq_{2}\alpha\nabla_{q_{2}}^{2}\tilde{f}_{2\psi}\right)+M^{2}\left(iq_{1}\alpha\tilde{f}_{1\varphi}\frac{\partial\nabla_{q_{2}}^{2}\tilde{f}_{2\varphi}}{\partial
y}-\frac{\partial\tilde{f}_{1\varphi}}{\partial
y}iq_{2}\alpha\nabla_{q_{2}}^{2}\tilde{f}_{2\varphi}\right),$ (44a)
$\displaystyle\tilde{n}_{\varphi}(\tilde{\bm{f}}_{1},\tilde{\bm{f}}_{2},q_{1},q_{2})=\left(iq_{1}\alpha\tilde{f}_{1\varphi}-\frac{\partial\tilde{f}_{1\psi}}{\partial
y}\right)iq_{2}\alpha\nabla_{q_{2}}^{2}\tilde{f}_{2\varphi}+\left(\frac{\partial\tilde{f}_{1\varphi}}{\partial
y}+iq_{1}\alpha\tilde{f}_{1\psi}\right)\frac{\partial\nabla_{q_{2}}^{2}\tilde{f}_{2\varphi}}{\partial
y}+\nabla_{q_{1}}^{2}\tilde{f}_{1\varphi}\nabla_{q_{2}}^{2}\tilde{f}_{2\varphi}.$
(44b)

Here, the second subscripts in $\tilde{f}_{1}$ and $\tilde{f}_{2}$ label the
corresponding components in the column vectors
$\tilde{\bm{f}}_{1}=(\tilde{f}_{1\psi},\tilde{f}_{1\varphi})^{T}$ and
$\tilde{\bm{f}}_{2}=(\tilde{f}_{2\psi},\tilde{f}_{2\varphi})^{T}$. That in
equation (28) for the amplitude expansion is in the same form as
$\tilde{\bm{N}}_{3}^{\circ}$ (but with different notations to differentiate
the two expansion methods)

$\tilde{\bm{N}}_{31}=\widetilde{\bm{N}}_{f}(\tilde{\bm{\gamma}}_{11},\tilde{\bm{\gamma}}_{20},1,0)+\widetilde{\bm{N}}_{f}(\tilde{\bm{\gamma}}_{20},\tilde{\bm{\gamma}}_{11},0,1)+\widetilde{\bm{N}}_{f}(\tilde{\bm{\gamma}}_{11}^{*},\tilde{\bm{\gamma}}_{22},-1,2)+\widetilde{\bm{N}}_{f}(\tilde{\bm{\gamma}}_{22},\tilde{\bm{\gamma}}_{11}^{*},2,-1).$
(45)

We also provide the corresponding expressions of the other two nonlinear terms
used in the solution process:
$\tilde{\bm{N}}_{20}=\widetilde{\bm{N}}_{f}(\tilde{\bm{\gamma}}_{11},\tilde{\bm{\gamma}}_{11}^{*},1,-1)+\widetilde{\bm{N}}_{f}(\tilde{\bm{\gamma}}_{11}^{*},\tilde{\bm{\gamma}}_{11},-1,1)$
and
$\tilde{\bm{N}}_{22}=\widetilde{\bm{N}}_{f}(\tilde{\bm{\gamma}}_{11},\tilde{\bm{\gamma}}_{11},1,1)$.

###  A.4 Adjoint operators

The weight matrix in the adjoint equation is self-adjoint
$\tilde{\bm{M}}_{0}^{\dagger}=\tilde{\bm{M}}_{0}$. The adjoint linear operator
is

$\tilde{\bm{L}}_{0}^{\dagger}=\begin{bmatrix}i\alpha\bar{U}\nabla_{1}^{2}+2i\alpha\bar{U}^{\prime}\frac{\partial}{\partial
y}+\frac{M^{2}}{T_{c}}\nabla_{1}^{4}&-\kappa^{2}\bar{\phi}^{\prime\prime\prime}i\alpha\\\
\frac{M^{2}}{\kappa^{2}}[i\alpha(\bar{\phi}^{\prime\prime\prime}+2\bar{\phi}^{\prime\prime}\frac{\partial}{\partial
y}+\bar{\phi}^{\prime}\frac{\partial^{2}}{\partial
y^{2}})-i\alpha^{3}\bar{\phi}^{\prime}-i\alpha\bar{\phi}^{\prime\prime\prime}]&S\end{bmatrix},$
(46)

where the coefficient $\kappa$ comes from the inner product in equation (19)
and the notation $S$ is short for

$S=i\alpha\bar{U}\nabla_{1}^{2}+i\alpha(\bar{U}^{\prime\prime}+2\bar{U}^{\prime}\frac{\partial}{\partial
y})-(\bar{\phi}^{\prime\prime\prime\prime}+\bar{\phi}^{\prime\prime\prime}\frac{\partial}{\partial
y})-(\bar{\phi}^{\prime}\nabla_{1}^{2}\frac{\partial}{\partial
y}+2\bar{\phi}^{\prime\prime}\frac{\partial^{2}}{\partial
y^{2}}+\bar{\phi}^{\prime\prime\prime}\frac{\partial}{\partial
y})+(\bar{\phi}^{\prime\prime}\nabla_{1}^{2}+2\bar{\phi}^{\prime\prime\prime}\frac{\partial}{\partial
y}+\bar{\phi}^{\prime\prime\prime\prime})+\frac{1}{Fe}\nabla_{1}^{4}.$ (47)

The adjoint eigenvector is
$\tilde{\bm{\gamma}}_{1}^{\dagger}=(\tilde{\psi}_{1}^{\dagger},\tilde{\varphi}_{1}^{\dagger})^{T}$.
The boundary conditions for $\tilde{\psi}_{1}^{\dagger}$ and
$\tilde{\varphi}_{1}^{\dagger}$ can be derived (from the integration by part
process) as

$\tilde{\psi}_{1}^{\dagger}(\pm
1)=\frac{\partial\tilde{\psi}_{1}^{\dagger}}{\partial y}(\pm 1)=0,$ (48a)
$\tilde{\varphi}_{1}^{\dagger}(\pm
1)=0,\,\,\,\frac{\partial\tilde{\varphi}_{1}^{\dagger}}{\partial
y}(1)=\frac{1}{Fe\bar{\phi}^{\prime}(1)}\frac{\partial^{2}\tilde{\varphi}_{1}^{\dagger}}{\partial
y^{2}}(1),\,\,\,\frac{\partial\tilde{\varphi}_{1}^{\dagger}}{\partial
y}(-1)=\frac{\tilde{\varphi}^{\prime}_{1}(-1)}{Fe\bar{\phi}^{\prime}(-1)\tilde{\varphi}^{\prime}_{1}(-1)+\tilde{\varphi}^{\prime\prime}_{1}(-1)}\frac{\partial^{2}\tilde{\varphi}_{1}^{\dagger}}{\partial
y^{2}}(-1).$ (48b)

More details regarding the derivation of boundary conditions here can be found
in our previous work (Zhang, 2016).

## References

  * Feng and Seyed-Yagoobi (2004) Y. Feng and J. Seyed-Yagoobi, “Understanding of electrohydrodynamic conduction pumping phenomenon,” Phys. Fluids 16, 2432–2441 (2004).
  * McCluskey  _et al._ (1991) F. M. J. McCluskey, P. Atten, and A. T. Perez, “Heat transfer enhancement by electroconvection resulting from an injected space charge between parallel plates,” Int. J. Heat Mass Transfer 34, 2237–2250 (1991).
  * Ward  _et al._ (2019) K. Ward, S. Matsumoto, and R. Narayanan, “The electrostatically forced Faraday instability: theory and experiments,” J. Fluid Mech. 862, 696–731 (2019).
  * Jalaal  _et al._ (2013) M. Jalaal, B. Khorshidi, and E. Esmaeilzadeh, “Electrohydrodynamic (EHD) mixing of two miscible dielectric liquids,” Chem. Eng. J. 219, 118–123 (2013).
  * Atten and Honda (1982) P. Atten and T. Honda, “The electroviscous effect and its explanation I - the electrohydrodynamic origin; study under unipolar DC injection,” J. Electrostat. 11, 225–245 (1982).
  * Soldati  _et al._ (1993) A. Soldati, P. Andreussi, and S. Banerjee, “Direct simulation of turbulent particle transport in electrostatic precipitators,” AIChE J. 39, 1910–1919 (1993).
  * Parker (1997) K. R. Parker, “Why an electrostatic precipitator?” in  _Applied Electrostatic Precipitation_ (Springer, 1997) pp. 1–10.
  * Mizuno (2000) A. Mizuno, “Electrostatic precipitation,” IEEE Trans. Dielectr. Electr. Insul. 7, 615–624 (2000).
  * Soldati (2000) A. Soldati, “On the effects of electrohydrodynamic flows and turbulence on aerosol transport and collection in wire-plate electrostatic precipitators,” J. Aerosol Sci. 31, 293–305 (2000).
  * Soldati (2003) A. Soldati, “Cost-efficiency analysis of a model wire-plate electrostatic precipitator via DNS based eulerian particle transport approach,” Aerosol Sci. Technol. 37, 171–182 (2003).
  * Jaworek  _et al._ (2018) A. Jaworek, A. Marchewicz, A. T. Sobczyk, A. Krupa, and T. Czech, “Two-stage electrostatic precipitators for the reduction of PM2.5 particle emission,” Prog. Energy Combust. Sci. 67, 206–233 (2018).
  * Castellanos (1998) A. Castellanos,  _Electrohydrodynamics_ , Vol. 380 (Springer Science & Business Media, 1998).
  * Schneider and Watson (1970) J. M. Schneider and P. K. Watson, “Electrohydrodynamic stability of space-charge-limited currents in dielectric liquids. I. theoretical study,” Phys. Fluids 13, 1948–1954 (1970).
  * Atten and Moreau (1972) P. Atten and R. Moreau, “Stabilité électrohydrodynamique des liquides isolants soumis à une injection unipolaire,” J. Méc. 11, 471–521 (1972).
  * Wu  _et al._ (2015) J. Wu, P. Traoré, A. T. Pérez, and P. A. Vázquez, “On two-dimensional finite amplitude electro-convection in a dielectric liquid induced by a strong unipolar injection,” J. Electrost. 74, 85–95 (2015).
  * Zhang  _et al._ (2015) M. Zhang, F. Martinelli, J. Wu, P. J. Schmid, and M. Quadrio, “Modal and non-modal stability analysis of electrohydrodynamic flow with and without cross-flow,” J. Fluid Mech. 770, 319–349 (2015).
  * Wang and Sheu (2016) B.-F. Wang and T. W.-H. Sheu, “Numerical investigation of electrohydrodynamic instability and bifurcation in a dielectric liquid subjected to unipolar injection,” Comput. Fluids 136, 1–10 (2016).
  * Lacroix  _et al._ (1975) J. C. Lacroix, P. Atten, and E. J. Hopfinger, “Electro-convection in a dielectric liquid layer subjected to unipolar injection,” J. Fluid Mech. 69, 539–563 (1975).
  * Atten and Lacroix (1978) P. Atten and J. C. Lacroix, “Electrohydrodynamic stability of liquids subjected to unipolar injection: Non linear phenomena,” J. Electrost. 5, 439–452 (1978).
  * Pérez and Castellanos (1989) A. T. Pérez and A. Castellanos, “Role of charge diffusion in finite-amplitude electroconvection,” Phys. Rev. A 40, 5844 (1989).
  * Feng  _et al._ (2021) Z. Feng, M. Zhang, P. A. Vazquez, and C. Shu, “Deterministic and stochastic bifurcations in two-dimensional electroconvective flows,” J. Fluid Mech. 922 (2021).
  * Castellanos and Agrait (1992) A. Castellanos and N. Agrait, “Unipolar injection induced instabilities in plane parallel flows,” IEEE Trans. Ind. Appl. 28, 513–519 (1992).
  * Lara  _et al._ (1997) J. L. Lara, A. Castellanos, and F. Pontiga, “Destabilization of plane Poiseuille flow of insulating liquids by unipolar charge injection,” Phys. Fluids 9, 399–406 (1997).
  * Landau (1944) L. D. Landau, “On the problem of turbulence.” C. R. Acad. Sci. URSS 44 (31), 1–314. (1944).
  * Stuart (1960) J. T. Stuart, “On the non-linear mechanics of wave disturbances in stable and unstable parallel flows Part 1. The basic behaviour in plane Poiseuille flow,” J. Fluid Mech. 9, 353–370 (1960).
  * Watson (1960) J. Watson, “On the nonlinear mechanics of wave disturbances in stable and unstable parallel flows,” J. Fluid Mech 9, 219 (1960).
  * Fujimura (1989) K. Fujimura, “The equivalence between two perturbation methods in weakly nonlinear stability theory for parallel shear flows,” Proc. R. Soc. London. A. Math. Phys. Sci. 424, 373–392 (1989).
  * Cross and Hohenberg (1993) M. C. Cross and P. C. Hohenberg, “Pattern formation outside of equilibrium,” Rev. Mod. Phys. 65, 851 (1993).
  * Felici (1971) N. Felici, “DC conduction in liquid dielectrics. Part II: electrohydrodynamic phenomena,” Direct Current Power Electron. 2, 147–165 (1971).
  * Worraker and Richardson (1981) W. J. Worraker and A. T. Richardson, “A nonlinear electrohydrodynamic stability analysis of a thermally stabilized plane layer of dielectric liquid,” J. Fluid Mech. 109, 217–237 (1981).
  * Zhang (2016) M. Zhang, “Weakly nonlinear stability analysis of subcritical electrohydrodynamic flow subject to strong unipolar injection,” J. Fluid Mech. 792, 328–363 (2016).
  * Pham and Suslov (2018) K. G. Pham and S. A. Suslov, “On the definition of Landau constants in amplitude equations away from a critical point,” R. Soc. Open Sci. 5, 180746 (2018).
  * Cudby and Lefauve (2021) J. Cudby and A. Lefauve, “Weakly nonlinear Holmboe waves,” Phys. Rev. Fluids 6, 024803 (2021).
  * Briggs (1964) R. J. Briggs,  _Electron-stream interaction with plasmas_ (MIT-Press, 1964).
  * Huerre and Monkewitz (1985) P. Huerre and P. A. Monkewitz, “Absolute and convective instabilities in free shear layers,” J. Fluid Mech. 159, 151–168 (1985).
  * Huerre and Monkewitz (1990) P. Huerre and P. A. Monkewitz, “Local and global instabilities in spatially developing flows,” Annu. Rev. Fluid Mech. 22, 473–537 (1990).
  * Chomaz (2005) J. M. Chomaz, “Global instabilities in spatially developing flows: Non-normality and nonlinearity,” Annu. Rev. Fluid Mech. 37, 357–392 (2005).
  * Carriere and Monkewitz (1999) P. Carriere and P. A. Monkewitz, “Convective versus absolute instability in mixed Rayleigh–Bénard–Poiseuille convection,” J. Fluid Mech. 384, 243–262 (1999).
  * Li  _et al._ (2019) F. Li, B.-F. Wang, Z.-H. Wan, J. Wu, and M. Zhang, “Absolute and convective instabilities in electrohydrodynamic flow subjected to a Poiseuille flow: a linear analysis,” J. Fluid Mech. 862, 816–844 (2019).
  * Delbende  _et al._ (1998) I. Delbende, J. M. Chomaz, and P. Huerre, “Absolute/convective instabilities in the Batchelor vortex: a numerical study of the linear impulse response,” J. Fluid Mech. 355, 229–254 (1998).
  * Delbende and Chomaz (1998) I. Delbende and J. M. Chomaz, “Nonlinear convective/absolute instabilities in parallel two-dimensional wakes,” Phys. Fluids 10, 2724–2736 (1998).
  * Chomaz (1992) J. M. Chomaz, “Absolute and convective instabilities in nonlinear systems,” Phys. Rev. Lett. 69, 1931 (1992).
  * Park and Huerre (1988) D. Park and P. Huerre, “On the convective nature of the Goertler instability,” Bull. Amer. Math. Soc. 33, 2552 (1988).
  * Chomaz and Perrier (1991) J. M. Chomaz and M. Perrier, “Nature of the Görtler instability: A forced experiment,” in  _The Global Geometry of Turbulence_ (Springer, 1991) pp. 23–32.
  * Chomaz (2003) J. M. Chomaz, “Fully nonlinear dynamics of parallel wakes,” J. Fluid Mech. 495, 57 (2003).
  * Suslov and Paolucci (1997) S. A. Suslov and S. Paolucci, “Nonlinear analysis of convection flow in a tall vertical enclosure under non-Boussinesq conditions,” J. Fluid Mech. 344, 1–41 (1997).
  * Leonard  _et al._ (1983) G. L. Leonard, M. Mitchner, and S. A. Self, “An experimental study of the electrohydrodynamic flow in electrostatic precipitators,” J. Fluid Mech. 127, 123–140 (1983).
  * Soldati and Banerjee (1998) A. Soldati and S. Banerjee, “Turbulence modification by large-scale organized electrohydrodynamic flows,” Phys. Fluids 10, 1742–1756 (1998).
  * Brancher and Chomaz (1997) P. Brancher and J. M. Chomaz, “Absolute and convective secondary instabilities in spatially periodic shear flows,” Phys. Rev. Lett. 78, 658 (1997).
  * Melville (1983) W. K. Melville, “Wave modulation and breakdown,” J. Fluid Mech. 128, 489–506 (1983).
  * Stewartson and Stuart (1971) K. Stewartson and J. T. Stuart, “A non-linear instability theory for a wave system in plane Poiseuille flow,” J. Fluid Mech. 48, 529–545 (1971).
  * Bender and Orszag (2013) C. M. Bender and S. A. Orszag,  _Advanced mathematical methods for scientists and engineers I: Asymptotic methods and perturbation theory_ (Springer Science & Business Media, 2013).
  * Luchini and Bottaro (2014) P. Luchini and A. Bottaro, “Adjoint equations in stability analysis,” Annu. Rev. Fluid Mech. 46, 493–517 (2014).
  * Suslov (2008) S. A. Suslov, “Two-equation model of mean flow resonances in subcritical flow systems,” Discrete Contin. Dyn. Syst. Series S, Selected Topics 1, 165–176 (2008).
  * Crouch and Herbert (1993) J.D. Crouch and Th. Herbert, “A note on the calculation of Landau constants,” Phys. Fluids A: Fluid Dynamics 5, 283–285 (1993).
  * Thompson  _et al._ (2001) M. C. Thompson, T. Leweke, and M. Provansal, “Kinematics and dynamics of sphere wake transition,” J. Fluids Struct. 15, 575–585 (2001).
  * Reynolds and Hussain (1972) W. C. Reynolds and A. K. M. F. Hussain, “The mechanics of an organized wave in turbulent shear flow. part 3. theoretical models and comparisons with experiments,” J. Fluid Mech. 54, 263–288 (1972).
  * Gao  _et al._ (2013) Z. Gao, A. Sergent, B. Podvin, S. Xin, P. Le Quéré, and L. S. Tuckerman, “Transition to chaos of natural convection between two infinite differentially heated vertical plates,” Phys. Rev. E 88, 023010 (2013).
  * Suslov and Paolucci (2004) S. A. Suslov and S. Paolucci, “Stability of non-Boussinesq convection via the complex Ginzburg–Landau model,” Fluid Dyn. Res. 35, 159 (2004).
  * Fischer  _et al._ (2008) P. F. Fischer, J. W. Lottes, and S. G. Kerkemeier, “Nek5000 web page,” (2008), available at: http://nek5000.mcs.anl.gov.
  * Patera (1984) A. T. Patera, “A spectral element method for fluid dynamics: laminar flow in a channel expansion,” J. Comput. Phys. 54, 468–488 (1984).
  * Maday and Patera (1989) Y. Maday and A. T. Patera, “Spectral element methods for the incompressible Navier-Stokes equations,” IN: State-of-the-art surveys on computational mechanics (A90-47176 21-64). New York , 71–143 (1989).
  * Guan and Novosselov (2019) Y. Guan and I. Novosselov, “Numerical analysis of electroconvection in cross-flow with unipolar charge injection,” Phys. Rev. Fluids 4, 103701 (2019).
  * van Saarloos and Hohenberg (1992) W. van Saarloos and P.C. Hohenberg, “Fronts, pulses, sources and sinks in generalized complex Ginzburg-Landau equations,” Phys. D 56, 303–367 (1992).
  * Fineberg and Steinberg (1987) J. Fineberg and V. Steinberg, “Vortex-front propagation in Rayleigh-Bénard convection,” Phys. Rev. Lett. 58, 1332 (1987).
  * Van Saarloos (2003) W. Van Saarloos, “Front propagation into unstable states,” Phys. Rep. 386, 29–222 (2003).
  * Dee and Langer (1983) G. Dee and J.S. Langer, “Propagating pattern selection,” Phys. Rev. Lett. 50, 383 (1983).
  * Ahlers and Cannell (1983) G. Ahlers and D. S. Cannell, “Vortex-front propagation in rotating Couette-Taylor flow,” Phys. Rev. Lett. 50, 1583 (1983).