# On the stability of waves in classically neutral flows

Colin Huber1, Meaghan Hoitt1, Nathaniel S. Barlow1, Nicole Hill2, Kimberlee
Keithley1,2, Steven J. Weinstein1,2 1 Department of Chemical Engineering,
Rochester Institute of Technology, Rochester, NY 14623  2 School of
Mathematical Sciences, Rochester Institute of Technology, Rochester, NY 14623
<EMAIL_ADDRESS>
###### Abstract

This paper reports a breakdown in linear stability theory under conditions of
neutral stability that is deduced by an examination of exponential modes of
the form $h\approx{{e}^{i(kx-\omega t)}}$, where $h$ is a response to a
disturbance, $k$ is a real wavenumber, and $\omega(k)$ is a wavelength-
dependent complex frequency. In a previous paper, King et al (Stability of
algebraically unstable dispersive flows, Phys. Rev. Fluids, 1(073604), 2016)
demonstrates that when Im$[\omega(k)]$=0 for all $k$, it is possible for a
system response to grow or damp algebraically as $h\approx{{t}^{s}}$ where $s$
is a fractional power. The growth is deduced through an asymptotic analysis of
the Fourier integral that inherently invokes the superposition of an infinite
number of modes. In this paper, the more typical case associated with the
transition from stability to instability is examined in which
Im$[\omega(k)]$=0 for a single mode (i.e., for one value of $k$) at neutral
stability. Two partial differential equation systems are examined, one that
has been constructed to elucidate key features of the stability threshold, and
a second that models the well-studied problem of rectilinear Newtonian flow
down an inclined plane. In both cases, algebraic growth/decay is deduced at
the neutral stability boundary, and the propagation features of the responses
are examined.

##  1 Introduction

The current work is motivated by the need to characterize fluid flows involved
in the manufacture of a variety of products. All flow processes are subjected
to time-varying disturbances induced from their surroundings, while product
quality typically necessitates the requirement of time invariance. In
particular, each product has a manufacturing tolerance to perturbations that
requires accurate characterization. Mathematical models that relate
disturbances to system responses are commonly used in industry to enable
guided experiments, and once validated, can be accurate enough to replace
experiments. An essential feature of disturbance modeling is whether the
underlying fluid system is stable or unstable, i.e. whether a response to
disturbances grows or decays. If the fluid system is stable, then product
specifications may be met by minimizing the magnitude of process disturbances;
if the system is not, control of the process is much more complex, and in some
cases, impossible. The characterization of fluid flow stability then, is
motivated by practical need. Within stable fluid systems, tight product
tolerances often dictate that perturbations be small and lie in a linear
regime – a product often becomes unsalable well within the dictates of
linearity. Thus, a corresponding linear operator may be used to model the
fluid system response. The characterization of the stability of such a linear
operator, then, is of paramount importance to develop sustainable operating
parameters that meet product specifications.

If a flowing liquid system in a manufacturing process is unstable, i.e.,
initiated disturbances are magnified, the flow can often easily be disrupted
away from a uniform state. This could be a desired outcome, as is the case for
liquid fuel atomizers [Ibrahim, 1995, Lin, 2003, El-Sayed et al., 2015] where
instability leads to the breakup of a liquid sheet into droplets, or in the
case of instability-driven turbulence to enhance mixing processes [Paul et
al., 2004]. Instability is unwelcome in other situations where layer
uniformity is essential, such as in the thin films used to coat ink jet and
copier papers, printed electronics, and liquid crystal display screens [Cohen
& Gutoff, 1992, Kistler & Schweizer, 1997]. Thus, it is important for
practitioners to control the parameters that influence fluid instability in
order to produce a salable product. The most widely used stability assessment,
referred to here as classical stability theory, provides the basis for much of
the hydrodynamic literature [Chandrasekhar, 1968, Huerre & Rossi, 1998,
Huerre, 2000, Schmid & Henningson, 2001, Drazin & Reid, 2004].

Classical stability theory is built on the following ideas. Typically, the
analysis starts with a full nonlinear operator and boundary conditions, which
are linearized about a base fluid flow (oftentimes exact) whose stability is
to be assessed. The resulting linear partial differential equation (PDE)
system may be expressed in terms of an operator, $L$, and flow response,
$h(x,t)$, as $Lh=0$, where $x$ and $t$ are respective space and time
variables. Note that forcing and initial conditions are neglected when
examining the PDE system, as the general homogeneous response characterizes
the classical stability of the medium. A solution of this equation on an
unbounded domain ($-\infty<x<\infty$) may be expressed as

$h=C{{e}^{i(kx-\omega t)}},$ (1)

where $k$ is a real wavenumber, $\omega$ is a complex frequency, a $C$ is an
amplitude. It is assumed that the fundamental responses (1) for each value of
$k$, henceforth called modes, may be superimposed to build the flow response
to any disturbance; furthermore, it is assumed that the stability of a complex
flow response may be characterized by the stability of its constituent modes.
Substitution of (1) into the linearized PDE system leads to a dispersion
relation of the form

$D(k,\omega)=0,~{}\textrm{ also written }\omega=\omega(k),$ (2)

that assures a nontrivial solution of the equation $Lh=0$. In classical
stability theory, the complex-valued $\omega(k)=\omega_{r}(k)+i\omega_{i}(k)$
is examined, where the function $\omega_{i}(k)$ determines the exponential
growth in time of a mode. The maximum growth rate over the range of
$k\in(-\infty,\infty)$, denoted as $\omega_{i,\textrm{max}}$, is used to
characterize the stability of the flow. At large times, the exponential nature
of the responses (1) dictate that growth rates at other wavenumbers are
subdominant, and the system thus grows as

$h\sim Ce^{\omega_{i,\textrm{max}}t},\textrm{ as }t\to\infty,$ (3a)

where $C$ is a constant. Thus, once $\omega=\omega(k)$ is established from
(2), the linear (exponential) stability is determined as follows
[Chandrasekhar, 1968, Huerre & Rossi, 1998, Huerre, 2000]:

$\omega_{i,\textrm{max}}<0\textrm{: the flow is linearly stable}$ (3b)
$\omega_{i,\textrm{max}}>0\textrm{: the flow is linearly unstable}$ (3c)
$\omega_{i,\textrm{max}}=0\textrm{: the flow is neutrally stable}.$ (3d)

The classification (3d) was used by [Rayeigh, 1880], and this classical
stability theory has been further developed over the last 100+ years.

The focus of this paper is the classification of neutral stability according
to (3d), and is motivated by [King et al., 2016], henceforth referred to as
KRK (the acronym for its first author). In that work, the classical stability
of a fluid flow system yields modes that are neutrally stable for all values
of $k$, i.e., $\omega_{i,\textrm{max}}$ = 0 for all wave numbers. While the
classical stability assessment (3d) indicates there will be no growth or decay
in a system response, KRK demonstrates both numerically and analytically that
disturbances can grow algebraically. Algebraic growth is defined as a system
response that obeys $h\sim Ct^{s}$, where $C$ is a constant and the exponent
$s$ is a positive rational number. It should be noted that algebraic decay of
disturbances ($s<0$) has also been identified in previous work when
exponential modes are neutrally stable for all real $k$, having both integer
[Case, 1960] and non-integer [Whitham, 2011, Lighthill, 2001, Barlow et al.,
2010] character. It is apparent that the classification of neutral stability
via classical means is deficient and warrants further study.

As is evident in the work of KRK and many earlier studies [De Luca & Costa,
1997, Barlow et al., 2011] (see KRK for a comprehensive literature review),
algebraic growth in linear PDE operators may be examined via a spatio-temporal
formulation involving both Fourier and Laplace transforms. When the inverse
Laplace transform is taken first, the resulting Fourier inversion integral has
an integrand with an exponential term identical to (1) with $\omega=\omega(k)$
according to the dispersion relation (2). As such, the connection between
classical stability analysis and a spatio-temporal analysis is made — the
assumption being that the growth characteristics of the superposition of modes
of the form (1) invoked via integration mimic exactly that of the individual
modes (1) when taken separately. In fact, KRK shows that this is not always
the case. When the exponent $s$ in $h\sim Ct^{s}$ is fractional, algebraic
growth can only be deduced via a superposition of modes via integration, and
is thus in fact “non-modal”. An additional important feature of systems
exhibiting algebraic growth is their sensitivity to the type of perturbations
to the system. KRK shows that if perturbations in initial velocity and forcing
are invoked, a system response may grow algebraically; however, if the same
system is perturbed in location, the system response can decay algebraically.
The sensitivity to initial conditions makes it imperative to consider a broad
set of perturbations to a system when establishing system stability,
especially so when examining the possibility of algebraic growth.

KRK focuses on a breakdown in classical stability theory for a dispersion
relation (2) in which all modes exhibit neutral stability, yet growth or decay
is actually predicted. The more typical situation in which neutral stability
arises in the literature is where a single mode is neutrally stable, and modes
corresponding to all other real wavenumbers exhibit damping according to the
form (1) with dispersion relation (2). It is widely held in prior literature
[Chandrasekhar, 1968, Schlichting, 1979, Drazin & Reid, 2004, Schmid &
Henningson, 2001] that flows exist in a state that is neither stable nor
unstable at a neutral stability boundary from classical theory. A question
arises as to whether the breakdown in classical stability theory (3d) reported
by KRK extends to the situations shown in Fig. 1 when only one wavenumber mode
is neutrally stable.

Fig. 1a is a schematic showing a scenario where this neutral stability
configuration occurs (see, for example [Manneville, 1998, Joulin & Vidal,
1998]); here, as a parameter (here, $B$) is varied through its critical value,
the flow changes from stable, to neutral, to unstable in accordance with the
characterization (3d). The features of such a transition is considered in
section 2 of this paper. Fig. 1b shows another type of transition that occurs
in the well-studied stability of a single layer Newtonian fluid flowing down
an inclined plane, where the critical parameter is the Reynolds number,
$R\\!e$, defined in section 3. As indicated in Fig. 1b, unstable flows for
$R\\!e>R\\!e_{c}$ have a maximum modal growth rate at finite wave number, and
as $R\\!e\to R\\!e_{c}$ this maximum growth rate diminishes to zero and its
associated wavenumber simultaneous approaches $k=0$. For $R\\!e\leq
R\\!e_{c}$, the $k=0$ mode is neutrally stable, and all other modes are damped
— in this case, there is not a true transition from instability to stability
based on a classical stability analysis, but rather a transition from
instability to neutral stability. This latter assessment has not been explored
in prior literature, as it is widely accepted that when $R\\!e<R\\!e_{c}$,
inclined plane flow is stable [Yih, 1963, L. Brevdo & Bridges, 1999]. It is
possible that the $k=0$ wavenumber is neglected in prior stability work
because it coincides with an interfacial perturbation that is flat, and thus
may be viewed as degenerate. The objective of this paper is to consider the
prospect of algebraic growth and decay for the classically neutral stability
configurations shown in Figs. 1a and 1b (at $B=B_{c}$ and $R\\!e=R\\!e_{c}$,
respectively).

(a)

(a)

(b)

(b)

Figure 1: Classical stability curves for two illustrative problems where
transition occurs due to the variation of a critical parameter: (a) The
transition from exponential growth to exponential decay for the model problem
in section 2 as the parameter $B$ is varied. The curves (top to bottom) in
correspond to $B$=$-1$, 0, and 1 (b) The effect of Reynolds number (defined in
section 3) on the stability transition for flow down an inclined plane. The
curves (top to bottom) correspond to $R\\!e$=1.1 $\cot\theta$, $\cot\theta$,
and 0.5 $\cot\theta$ where (for all curves) the incline angle is
$\theta=\pi/4$ and Weber number (defined in section 2) is 0.1.

The paper is organized as follows. Section 2 introduces a model PDE for the
transition shown in Fig. 1a, where it is shown that algebraic growth can occur
on the threshold of neutral stability. The algebraic growth is extracted via
long-time asymptotic analysis of the Fourier integral solution, which is then
compared with the Fourier series solution of the PDE. The spatio-temporal
classification of algebraic growth as being absolutely or convectively
unstable is examined. In section 3, the same approach is used to examine the
well-known incline plane flow studied by [Yih, 1963]. Here algebraic decay is
shown to occur on the neutral stability threshold. In section 3.1, the
governing PDE is derived in non-dimensional form. In section 3.2, the
classical stability analysis is reviewed and key elements are extracted
regarding the neutral stability threshold. In section 3.3, the integral
solution is examined via asymptotic analysis and the algebraic decay rate is
deduced. A summary of our results and concluding remarks are given in section
4. Key supporting analyses are provided in Appendices A-C.

##  2 Model Problem that exhibits algebraic growth

In this section, we focus on waves with vertical displacement $h(x,t)$,
described by the following PDE,

$\displaystyle\frac{\partial^{2}h}{\partial
t^{2}}+2c\frac{\partial^{2}h}{\partial x\partial
t}+c^{2}\frac{\partial^{2}h}{\partial x^{2}}+\frac{\partial^{4}h}{\partial
x^{4}}-\frac{\partial^{3}h}{\partial x^{2}\partial
t}-c\frac{\partial^{3}h}{\partial x^{3}}+Bc\frac{\partial h}{\partial
x}+B\frac{\partial h}{\partial t}+B^{2}h=f_{0}\delta(x)\delta(t)$
$\displaystyle h(x,0)=h_{0}\delta(x),~{}~{}\frac{\partial h}{\partial
t}(x,0)=u_{0}\delta(x),~{}~{}-\infty<x<\infty,~{}~{}t>0,$ (4)

where $x$ is the horizontal coordinate, $t$ is time, $B$ is an instability
parameter and $c$ is a convective parameter. In (4), the parameters $h_{0}$,
$u_{0}$, and $f_{0}$ are magnitudes of the initial conditions and forcing. The
use of the delta function as an efficient means to initiate disturbances is
well-established [King et al., 2016]. The PDE (4) was not physically
motivated; rather, it was reverse engineered from a dispersion relation that
leads to the “textbook” [Manneville, 1998] depiction of instability transition
as shown in Fig. 1a. This led to the expression (4) but with $c=0$. The
parameter $c$ was later added to incorporate an element of convection without
affecting the stability character of Fig. 1a.

###  2.1 Classical stability analysis & comparison with Fourier series

Substituting $h=C{{e}^{i(kx-\omega t)}}$ into the homogeneous version of (4)
leads to the dispersion relation

$\omega=\frac{-2Ck+(k^{2}+B)i\pm\sqrt{3k^{4}-2k^{2}B+3B^{2}}}{-2},$ (5)

whose imaginary part $\omega_{i}$ is plotted versus real $k$ in Fig. 1 for
$C=0$ and $B$=$-1$ (top curve), 0 (middle), and 1 (bottom), indicating
(classical) temporal growth rates of $\omega_{i,\textrm{max}}$=0.5, 0, and
$-0.5$ respectively. These (exponential) growth rates are validated by the
Fourier series solution of (4), as shown in Fig. 2a where the peak of the
impulse response is plotted versus time. It is clear from the figure that the
classical predictions of exponential growth ($B=-1$,
$\omega_{i,\textrm{max}}$=0.5) and decay ($B=1$,
$\omega_{i,\textrm{max}}=-$0.5) accurately describe the long-time behavior of
the response. However, for the classically neutral case of $B=0$, an
exponential growth rate of $\omega_{i,\textrm{max}}$=0 does not, in fact,
yield zero growth, as shown by the inset in Fig. 2a. The log-log plot of Fig.
2b indicates that the peak grows like $\approx t^{1/2}$ when $B=0$. We now
proceed to show that the algebraic growth rate is indeed exactly 1/2 via
asymptotic analysis of the integral solution of (4).

(a)

(a)

(b)

(b)

Figure 2: (a) The evolution of the magnitude of the wave peak vs. time in the
Fourier series solution of (4) for $c=0$ and various $B$ values. The slopes on
this semi-log plot confirm the exponential growth rates predicted by classical
stability analysis. (b) The $B=0$ case from (a) instead shown on a log-log
plot, indicating algebraic growth as $t^{1/2}$.

###  2.2 Exact and asymptotic solutions

Here we focus solely the case of classical neutral stability as shown in Fig.
2 where $B=B_{c}=0$. In this case, the system (4) becomes:

$\displaystyle\frac{\partial^{2}h}{\partial
t^{2}}+2c\frac{\partial^{2}h}{\partial x\partial
t}+c^{2}\frac{\partial^{2}h}{\partial x^{2}}+\frac{\partial^{4}h}{\partial
x^{4}}-\frac{\partial^{3}h}{\partial x^{2}\partial
t}-c\frac{\partial^{3}h}{\partial x^{3}}=f_{0}\delta(x)\delta(t),$
$\displaystyle h(x,0)=h_{0}\delta(x),~{}~{}\frac{\partial h}{\partial
t}(x,0)=u_{0}\delta(x),~{}~{}-\infty<x<\infty,~{}~{}t>0,$ (6)

To solve the system (6), the Fourier transform and its inverse are utilized;
these are given respectively by equations (7) and (8) as follows

$\hat{h}\left(k,t\right)=\int\limits_{-\infty}^{\infty}{h(x,t){{e}^{-ikx}}dx},$
(7)
$h\left(x,t\right)=\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}{\hat{h}(k,t){{e}^{ikx}}dx}.$
(8)

The Fourier transform of (6) yields

$\displaystyle\frac{d^{2}\hat{h}}{dt^{2}}+\left(k^{2}+2ick\right)\frac{d\hat{h}}{dt}+\left(k^{4}+ick^{3}-c^{2}k^{2}\right)\hat{h}=f_{0}\delta(t),$
$\displaystyle\hat{h}(0)=h_{0},~{}~{}\frac{d\hat{h}}{dt}(0)=u_{0},~{}~{}t>0,$
(9)

where $\hat{h}(k,t)$ denotes the Fourier transform of $h(x,t)$. Equation (9)
is a linear constant coefficient ordinary differential equation,a whose
solution is

$\displaystyle\hat{h}(k,t)=e^{-ickt}\left[F_{1}(k)e^{-r_{1}k^{2}t}+F_{2}(k)e^{-r_{2}k^{2}t}\right],$
$\displaystyle
F_{1,2}(k)=\frac{\mp(u_{0}+f_{0})i}{\sqrt{3}k^{2}}\pm\frac{h_{0}c}{\sqrt{3}k}\mp\frac{h_{0}}{6}\left(i\sqrt{3}\mp
3\right),~{}~{}r_{1,2}=\frac{1\mp i\sqrt{3}}{2}$ (10)

The inverse Fourier transform of (10) is

$\displaystyle\displaystyle h(x,t)$
$\displaystyle=\frac{1}{2\pi}\int_{-\infty}^{\infty}\hat{h}(k,t)e^{ikx}~{}dk,$
(11)
$\displaystyle=\frac{1}{2\pi}\int_{-\infty}^{\infty}\left[F_{1}(k)e^{-r_{1}k^{2}t}+F_{2}(k)e^{-r_{2}k^{2}t}\right]e^{ik\left(\frac{x}{t}-c\right)t}~{}dk,$

where the quantity $t$ is factored out of the exponential and the “velocity”
quantity $x/t$ is introduced to highlight that the structure of (11) and
associated integration technique to follow depends on whether $x/t=c$ or
$x/t\neq c$. Note that we obtain the same result given by (11) if we utilize
the Fourier-Laplace double inversion, as done in [King et al., 2016] (and
references therein). As is standard for methods such as steepest descent and
stationary phase, we construct a long-time asymptotic solution of (11) by
examining the integrals along a fixed velocities $x/t$ as $t\to\infty$. Before
doing this, we first split (11) into 3 integrals, sorted by order of the $k$
singularities in the integrand:

$\displaystyle\displaystyle h(x,t)$
$\displaystyle=\frac{i(u_{0}+f_{0})}{2\sqrt{3}\pi}\int_{-\infty}^{\infty}\frac{e^{-r_{2}k^{2}t}-e^{-r_{1}k^{2}t}}{k^{2}}e^{ik\left(\frac{x}{t}-c\right)t}~{}dk$
(12)
$\displaystyle+\frac{h_{0}c}{2\sqrt{3}\pi}\int_{-\infty}^{\infty}\frac{e^{-r_{1}k^{2}t}-e^{-r_{2}k^{2}t}}{k}e^{ik\left(\frac{x}{t}-c\right)t}~{}dk$
$\displaystyle+\frac{h_{0}}{12\pi}\int_{-\infty}^{\infty}\left[\left(3-i\sqrt{3}\right)e^{-r_{1}k^{2}t}+\left(3+i\sqrt{3}\right)e^{-r_{2}k^{2}t}\right]e^{ik\left(\frac{x}{t}-c\right)t}~{}dk.$

As evidenced by the coefficient of the first integral in (12), the delta
function forcing utilized yields a response equivalent to that for a
perturbation in velocity at $t=0$. It is useful to decompose the imaginary
exponentials of the first and third integrals of (12) into sines and cosines
via Euler’s relation while leaving the second integral as is, since identities
for these integrals are readily available [Gradshteyn & Ryzhik, 2007]. After
simplifying the first and third integrals, (12) is rewritten as

$\displaystyle\displaystyle h(x,t)=$
$\displaystyle\frac{2(u_{0}+f_{0})}{\sqrt{3}\pi}\int_{0}^{\infty}\frac{e^{-\frac{1}{2}k^{2}t}}{k^{2}}\sin\left(\frac{\sqrt{3}}{2}k^{2}t\right)\cos\left[k\left(\frac{x}{t}-c\right)t\right]~{}dk$
$\displaystyle+\frac{h_{0}c}{2\sqrt{3}\pi}\int_{-\infty}^{\infty}\frac{e^{-r_{1}k^{2}t}-e^{-r_{2}k^{2}t}}{k}e^{ik\left(\frac{x}{t}-c\right)t}~{}dk$
$\displaystyle+\frac{h_{0}}{3\pi}\int_{0}^{\infty}e^{-\frac{1}{2}k^{2}t}\left[3\cos\left(\frac{\sqrt{3}}{2}k^{2}t\right)+\sqrt{3}\sin\left(\frac{\sqrt{3}}{2}k^{2}t\right)\right]\cos\left[k\left(\frac{x}{t}-c\right)t\right]~{}dk.$

The three integrals in (LABEL:eq:trig) may be evaluated exactly using the
identities given by (LABEL:eq:Igeneral), (49), and (A.1) provided in Appendix
A to yield

$\displaystyle h(x,t)=$
$\displaystyle~{}\frac{\left(u_{0}+f_{0}\right)i\left(\frac{x}{t}-c\right)t}{4\sqrt{3\pi}}\left\\{\Gamma\left[-\frac{1}{2},\left(1+\sqrt{3}i\right)\frac{\left(\frac{x}{t}-c\right)^{2}t}{8}\right]-\Gamma\left[-\frac{1}{2},\left(1-\sqrt{3}i\right)\frac{\left(\frac{x}{t}-c\right)^{2}t}{8}\right]\right\\}$
$\displaystyle+~{}\frac{h_{0}c}{\sqrt{3}}~{}\textrm{Im}\left[\textrm{erf}\frac{\left(\frac{x}{t}-c\right)\sqrt{t}\left(\sqrt{3}-i\right)}{4}\right]+~{}\frac{h_{0}\cos\left[\frac{\left(\frac{x}{t}-c\right)^{2}t\sqrt{3}}{8}\right]}{\exp\left[\frac{\left(\frac{x}{t}-c\right)^{2}t}{8}\right]\sqrt{3\pi
t}},$ (14)

where $\Gamma$ is the upper incomplete gamma function.

The response (14) drastically simplifies along the ray $x/t=c$. Applying the
limit as $x/t\to c$ to the aggregate first term of (14) leads to
$(u_{0}+f_{0})\sqrt{t/(3\pi)}$. The second term in (14) is zero for $x/t=c$,
and one may directly substitute $x/t=c$ into the third term of (14) to obtain
$h_{0}/\sqrt{3\pi t}$; collecting these results, (14) reduces to

$\displaystyle
h(x,t)|_{\frac{x}{t}=c}=\frac{1}{\sqrt{3\pi}}\left[(u_{0}+f_{0})t^{\frac{1}{2}}+h_{0}t^{-\frac{1}{2}}\right].$
(15)

The large time behavior of $h$ for $x/t\neq c$ is found by expanding the first
two terms of (14) as $t\to\infty$; these expansions are found in [Abramowitz &
Stegun, 1972] as Eqns. 6.5.32 and 7.1.23, respectively. After replacing the
first and second terms of (14) with their leading-order behavior as
$t\to\infty$ and then factoring out the decaying exponential appearing in all
three terms, we obtain

$\displaystyle h(x,t)|_{\frac{x}{t}\neq
c}\sim\frac{e^{-\left(\frac{x}{t}-c\right)^{2}t/8}}{\sqrt{3\pi t}}$
$\displaystyle\left\\{\left[\frac{4(u_{0}+f_{0})}{\left(\frac{x}{t}-c\right)^{2}}-\frac{h_{0}c}{\left(\frac{x}{t}-c\right)}+h_{0}\right]\cos\left[\frac{\sqrt{3}\left(\frac{x}{t}-c\right)^{2}t}{8}\right]\right.$
(16)
$\displaystyle~{}\left.-~{}\frac{h_{0}c\sqrt{3}}{\left(\frac{x}{t}-c\right)}\sin\left[\frac{\sqrt{3}\left(\frac{x}{t}-c\right)^{2}t}{8}\right]\right\\},~{}~{}t\to\infty.$

To summarize equations (15) and (16), the response grows algebraically like
$t^{\frac{1}{2}}$ for $x/t=c$ (first term of (15), assuming $u_{0}+f_{0}\neq
0$) and decays exponentially for $x/t\neq c$ (all terms of (16)). The Fourier
series solution to (4) is shown in Fig. 3 for $u_{0}+f_{0}=1$. All Fourier
series solutions displayed in this paper are constructed on a periodic domain,
following the approach given in [Barlow et al., 2010]. In Fig. 3 the peak is
shown to be growing in accordance with the exact solution given by (15) along
$x/t=c$, as indicated by a black line in the figure.

(a)  Figure 3: (a) Fourier series solution (shaded surface) to (6) compared
with the exact solution (15) (black curve) along the line $x/t=c$.
$h_{0}=u_{0}=c=1$, $f_{0}$=0.

Fig. 3 gives the appearance of growth that is spreading as opposed to being
confined to a single peak moving at $x/t=c$. However, the growing response
that we see is, in fact, mostly made of transient (“short”-time) growth,
which, along any given ray $x/t\neq c$ will eventually damp as $t\to\infty$.
Points moving near the peak of the response at velocities closer to $c$ will
take a longer time to damp, as can be seen in Fig. 4, where the response is
tracked along specific $x/t$ rays and compared with the exact solution
($\bullet$’s) given by (15) for $x/t=c$ and the long-time asymptotic solution
(dashed lines) given by (16).

(a)  (b)

Figure 4: Fourier series solution (solid curves) to (6) along specific $x/t$
rays, compared with exact solution given by (15) (top, denoted by $\bullet$)
along $x/t=c$ and leading-order asymptotic solution given by (16) (bottom
dashed curves) for $x/t\neq c$. $h_{0}=u_{0}=c=1$, $f_{0}$=0.

Although the growth/decay behavior shown here in Figs. 3 and 4 is similar to
the algebraic growth behavior of KRK (where are all classical modes are
neutral), there is one key difference; Here, the response along non-growing
rays decays exponentially, allowing for the non-oscillatory dome-like
structure shown in Fig. 3. For the problem of KRK, such non-growing waves
decay algebraically, allowing for oscillatory effects to persist, leading to a
waveform that has many ripples that spread away from the growing peak (see
figures 1-5 in [King et al., 2016]). Despite these differences, it is worth
pointing out that the algebraic growth/decay character in both (6) and the
problem of KRK arises in the evaluation of Fourier integral solutions that
contain removable singularities, such as those shown in (LABEL:eq:trig). Such
singularities (both here and in KRK) require special attention, as it is often
not possible to interpret the integrals as principal values when the
integrands are even. Appendix A provides key steps to obtain either exact or
asymptotic solutions containing algebraic growth/decay.

Note that if $u_{0}+f_{0}=0$ and $h_{0}\neq 0$, then the response decays
algebraically for $x/t=c$ (second term of (15)) and decays exponentially for
$x/t\neq c$. Thus, an impulse can be introduced to a system, either through
$u_{0}$, $f_{0}$, or $h_{0}$, and this choice will affect the algebraic
stability or instability of the response. As shown previously in [Barlow et
al., 2011] and [King et al., 2016], algebraically growing systems are
particularly sensitive to the choice of initiating disturbance, as further
evidenced by the results above.

Besides being a model problem that illustrates algebraic growth, the PDE given
by (6) is also useful as a model for either convective instability (if $c\neq
0$, see Fig. 3) or absolute instability (if $c=0$), provided that growth is
enabled (i.e., $u_{0}+f_{0}\neq 0$). Convective instability is defined by a
response $h$ that grows and convects for all time (i.e., along a ray $x/t\neq
0$) but also decays at any fixed location for large time (i.e., along the ray
$x/t=0$) [Huerre & Rossi, 1998]; this is described, respectively, by (15) and
(16) for $c\neq 0$ and $u_{0}+f_{0}\neq 0$. Absolute instability is defined by
a response $h$ that grows at any fixed location for large time (i.e., along
the ray $x/t=0$) [Huerre & Rossi, 1998]; this is described, respectively, by
(15) and (16) for $c=0$ and $u_{0}+f_{0}\neq 0$. A comparison between
algebraic absolute/convective instability behavior and exponential
absolute/convective instability behavior is given in [King et al., 2016].

##  3 Algebraically decaying waves in a thin film along an inclined solid

###  3.1 Governing equations valid for small interfacial slope

We consider a Newtonian liquid of constant density, $\rho$, and constant
viscosity, $\mu$, flowing under the influence of gravity, $g$, along a solid
surface inclined to horizontal with angle $\theta$ as shown in Fig. 5. The
liquid layer is exposed to air having a constant atmospheric pressure, and the
air-liquid interface has a constant surface tension, $\sigma$. Under ideal
conditions, the liquid flows with a steady-state constant thickness $H_{0}$
(Fig. 5) and constant volumetric flow rate per unit width, $Q_{0}$. The film
is perturbed away from uniform due to external disturbances (to be specified
later) while remaining invariant in the $Z$ direction (out of Fig. 5), and the
total film thickness deviates from $H_{0}$ in accordance with the
parameterization $Y=H_{0}h(X,T)$, where $h$ is a dimensionless multiplier, $Y$
is the distance perpendicular to the wall into the fluid domain, $X$ is the
distance down the incline, and $T$ is time. The dynamics of the air are
neglected and the external pressure remains atmospheric for all time.

Figure 5: Schematic of flow of a thin film along an inclined solid. The flow
is invariant in the $Z$ direction, oriented into and out of the figure. The
steady state interface location is denoted as $H_{0}$, the perturbed interface
location denoted as $H$, and the perturbation magnitude given as $h$, a
fraction of $H_{0}$.

A set of approximate dynamical equations that govern this configuration may be
developed as follows [S. V. Alekseenko & Pokusaev, 1994, Weinstein & Ruschak,
2004]. It is assumed that the slope of the perturbed air-liquid interface and
underlying fluid trajectories are small. However, in contrast to lubrication
theory that utilizes such assumptions, inertial effects are retained since the
flow may be rapid. The simplified time-dependent Navier-Stokes equations and
continuity equation are integrated across the film thickness $H$, and boundary
conditions at the wall (no slip and kinematic conditions), and interface
(dynamic condition simplified for small slope, kinematic condition) are
applied. The result is an integral equation equivalent to the simplified
system of equations and boundary conditions. Following the approach pioneered
by Von Karmen/Polhausen in their treatment of boundary layer theories [von
Kármán, 1921, Pohlhausen, 1921, Schlichting, 1979], a velocity profile
parabolic in $Y$ (see coordinate system in Fig. 5) is assumed that satisfies
the wall boundary conditions, a no shear condition at the unknown interface
location $H(X,T)$, and integrates to yield the local volumetric flow per
width, $Q(X,T)$. The resulting system of equations may be written in
dimensionless form as [S. V. Alekseenko & Pokusaev, 1994]:

$\frac{\partial\bar{Q}}{\partial x}+\frac{\partial\bar{H}}{\partial t}=0$
(17a) $\frac{\partial\bar{Q}}{\partial t}+\frac{6}{5}\frac{\partial}{\partial
x}\left(\frac{{{{\bar{Q}}}^{2}}}{{\bar{H}}}\right)=\frac{3}{R\\!e}\left(\bar{H}-\frac{{\bar{Q}}}{{{{\bar{H}}}^{2}}}-\bar{H}\frac{\partial\bar{H}}{\partial
x}\cot\theta\right)+\frac{1}{W\\!e}\bar{H}\frac{{{\partial}^{3}}\bar{H}}{\partial{{x}^{3}}}$
(17b)

where:

$\bar{H}=\frac{H}{{{H}_{0}}},~{}\bar{Q}=\frac{Q}{{{Q}_{0}}},~{}x=\frac{X}{{{H}_{0}}},~{}t=\frac{T{{Q}_{0}}}{{{H}_{0}}^{2}},~{}R\\!e=\frac{\rho{{Q}_{0}}}{\mu},~{}W\\!e=\frac{\rho{{Q}_{0}}^{2}}{\sigma{{H}_{0}}}.$
(17c)

Note that $H_{0}$ and $Q_{0}$ are not independent and are related as follows

$H_{0}=\left(\frac{3Q_{0}\mu}{\rho g\sin\theta}\right)^{1/3},$ (17d)

where $H_{0}$ is the exact expression for the steady-state thickness for film
flow calculated from the Navier-Stokes Equations [Pritchard, 2011].

The nonlinear system (5) can further be simplified by restricting attention to
small interfacial perturbations, which is justified in practical applications
where highly uniform films are desired [Weinstein & Ruschak, 2004]. The
following forms are assumed:

$\bar{H}\sim 1+h(x,t),~{}\bar{Q}\sim 1+q(x,t),~{}h<<1\textrm{ and }q<<1.$ (18)

Equation (18) is substituted into the system (5) and terms quadratic or higher
in the perturbation quantities $h$ and $q$ are neglected. The two linearized
equations corresponding to (17a) and (17b) are combined into a single equation
to yield:

$\frac{{{\partial}^{2}}h}{\partial{{t}^{2}}}+\frac{3}{R\\!e}\frac{\partial
h}{\partial t}+\frac{12}{5}\frac{{{\partial}^{2}}h}{\partial x\partial
t}+\frac{1}{W\\!e}\frac{{{\partial}^{4}}h}{\partial{{x}^{4}}}+\left(\frac{6}{5}-\frac{3\cot\theta}{R\\!e}\right)\frac{{{\partial}^{2}}h}{\partial{{x}^{2}}}+\frac{9}{R\\!e}\frac{\partial
h}{\partial x}=0.$ (19a)

This is the desired governing equation that will be examined in what follows.
In accordance with the geometry in Fig. 5, the spatial domain that will be
considered is $-\infty<x<\infty$, and the following boundary conditions and
initial conditions are applied:

$h={{h}_{0}}\delta(x),~{}\frac{\partial h}{\partial
t}={{u}_{0}}\delta(x)\textrm{ at }t=0$ (19b) $h\to 0\textrm{ as
}x\to\pm\infty.$ (19c)

In (18), $h_{0}$ and $u_{0}$ are constants and $\delta(x)$ is the Dirac delta
function. The system (18) is well posed to solve for the response $h(x,t)$.
Note that, although (19a) is homogeneous, including an impulsive pressure
forcing of $f_{0}\delta(x)\delta(t)$ would have the same response as the
initial velocity condition imposed in (19b) as seen for the problem of section
2; thus, we have not incorporated it here.

###  3.2 Classical stability analysis

We now proceed to examine the classical stability of equation (19a). To do so,
boundary conditions (19b) and (19c) are neglected, and the following
disturbance form is assumed:

$h=A{{e}^{i(kx-\omega t)}}.$ (20)

In (20), $k$ is a real wavenumber, $\omega=\omega_{r}+i\omega_{i}$ is a
generally complex frequency, and $A$ is a constant. Substituting (20) into
(19a) and rearranging to assure a non-trivial solution leads to the following:

${{\omega}^{2}}+i\beta\omega-\gamma=0.$ (21a)
$\beta=\left(\frac{3}{R\\!e}+\frac{12i}{5}k\right),~{}\gamma=\frac{{{k}^{4}}}{W\\!e}-{{k}^{2}}\left(\frac{6}{5}-\frac{3\cot\theta}{R\\!e}\right)+\frac{9ik}{R\\!e}$
(21b)

As shown by [Yih, 1963], the neutral stability condition may be deduced by
examining the long wavelength limit as a perturbation series about $k\to 0$
holding $c=\omega/k$ fixed, which assures that the speed of any disturbance
given by the real part of $c$ is finite. Equation (21b) is first rewritten for
fixed $c$ by inserting $\omega=ck$, and the following expansion is inserted
into the result:

$c\sim{{c}_{0}}+{{c}_{1}}k+c_{2}k^{2}+c_{3}k^{3}+O(k^{4})\textrm{ as }k\to
0,~{}c=\omega/k\textrm{ fixed }$ (22)

where equating like powers leads to

$c_{0}=3$ (23a) ${{c}_{1}}=i\left(R\\!e-\cot\theta\right).$ (23b)
${{c}_{2}}=\frac{6}{5}R\\!e\left(\cot\theta-R\\!e\right).$ (23c)
${{c}_{3}}=i\left[\frac{36}{25}R\\!e^{2}\left(\cot\theta-R\\!e\right)-\frac{12}{5}R\\!e^{3}\left(\cot\theta-R\\!e\right)^{2}-\frac{Re}{3W\\!e}\right].$
(23d)

Finally, the results (23d) can be rewritten using the definition of $c$ in
(22) to obtain the final form of $\omega$ as:

$\displaystyle\omega\sim$ $\displaystyle
3k+i(R\\!e-\cot\theta)k^{2}+\frac{6}{5}R\\!e\left(\cot\theta-R\\!e\right)k^{3}$
$\displaystyle+i\left[\frac{36}{25}R\\!e^{2}\left(\cot\theta-R\\!e\right)-\frac{12}{5}R\\!e^{3}\left(\cot\theta-R\\!e\right)^{2}-\frac{Re}{3W\\!e}\right]k^{4}+O(k^{5})\textrm{
as }k\to 0.$

According to the form (20), the response grows exponentially if
$\omega_{i}>0$, and thus according to (LABEL:eq:17) for $k\neq 0$, there is
growth when $R\\!e>\cot\theta$. The approximate nature of equation (18) is
revealed here, for when the full linearized Navier-Stokes system is analyzed
in the long wavelength limit, the exact result from [Yih, 1963] is that
instabilities arise for $R\\!e>5/6\cot\theta$; apart from the coefficient
difference, the interpretation of instability is identical to that of the
approximate analysis. Note that, to make his stability assessment, Yih only
retained terms to order $k^{2}$ in (LABEL:eq:17) (i.e., $O$($k$) in (22)).
While the additional terms in (LABEL:eq:17) do not alter the above stability
conclusions, they are required to enable the analysis in section 3.3 and in
supporting analysis found in Appendix C.

The classical interpretation of exponential growth in the long wavelength
limit warrants additional discussion relevant to the current work. In
particular, equation (LABEL:eq:17) indicates that $\omega_{i}=0$ when $k=0$,
and thus equation (20) shows that this wavenumber exhibits a neutrally stable
condition. Fig. 6 provides a schematic of the growth rate, $\omega_{i}$, as a
function of $k$, for $k\geq 0$, based on solutions of the Orr-Sommerfeld
equations (see, for example [Weinstein, 1990, L. Brevdo & Bridges, 1999]). As
indicated in curve (a), for $R\\!e>\cot\theta$, waves exponentially grow over
a range of $k$, with the maximum wave growth occurring at finite $k$. For
smaller $R\\!e$ that satisfies $R\\!e>\cot\theta$ given by curve (b), the
wavenumber associated with maximum growth is reduced, until at the neutral
condition shown by curve (c), this wavenumber becomes coincident with $k=0$;
for all $k\neq 0$, modes damp exponentially in accordance with (20). This same
structure persists for $R\\!e<\cot\theta$. For all curves shown in Fig. 6,
note that $k=0$ is in fact a wavenumber that exibits neither growth nor decay,
from the classical characterization.

Figure 6: Typical plot of growth rate, $\omega_{i}$, vs real wavenumber, $k$.
In the $k\to 0$ limit, $\omega_{i}$ is given by equation (LABEL:eq:17). Curve
a): $R\\!e>\cot\theta$, an unstable condition. Note that the maximum growth
rate occurs at a finite wavenumber. Curve b): Smaller $R\\!e$ than curve a,
but still satisfying $R\\!e>\cot\theta$. Note that the maximum growth rate is
reduced and has shifted to smaller wavenumber. Curve c): The first condition
at which there is no growth as $R\\!e$ is lowered is $R\\!e=\cot\theta$. Curve
d): $R\\!e<\cot\theta$. For curves c) and d) Modes exponentially damp for all
$k\neq 0$, and $k=0$ is the wavenumber of maximum (but zero) growth.

The Orr-Sommerfeld solutions in Fig. 6 show why a stability assessment for
$R\\!e>\cot\theta$ may be made based on the $k\to 0$ limit. Although the
maximum growth in this parameter range generally occurs at a finite value of
$k$, the concavity of the asymptotic solution (LABEL:eq:17) is always positive
when there is growth; this enables positive concavity in the $k\to 0$ limit to
assure instability. For situations where $R\\!e=\cot\theta$ or
$R\\!e<\cot\theta$, the concavity of (LABEL:eq:17) is zero or negative,
respectively. For these cases, however, the $k=0$ wavenumber precludes a
stability conclusion (as pointed out by [Yih, 1963]), since $k=0$ corresponds
to the wavenumber of maximum (zero) growth. We set out to make this assessment
in the analysis to follow, where we explicitly examine the system response for
cases where $R\\!e\leq\cot\theta$ via asymptotic analysis.

###  3.3 Integral solution and long-time asymptotic behavior

We now proceed to solve the system (18) and examine its stability features,
following the process discussed in section 2.2. The Fourier transform (7) of
(18) yields

$\frac{{{d}^{2}}\hat{h}}{d{{t}^{2}}}+\beta\frac{d\hat{h}}{dt}+\gamma\hat{h}=0.$
(25a)

It is implicit in the use of the Fourier transform that equation (19c) is
satisfied. The associated initial conditions (19b) become:

$\hat{h}={{h}_{0}},~{}\frac{\partial\hat{h}}{\partial t}={{u}_{0}}\textrm{ at
}t=0.$ (25b)

The ordinary differential system (25b) is solved to yield

$\hat{h}=C_{1}(k)e^{-i\omega_{1}(k)t}+C_{2}e^{-i\omega_{2}(k)t},$ (26)

where $\omega_{1}(k)$ and $\omega_{2}(k)$ are the two roots of the quadratic
equation (21b) given as

${{\omega_{1}}}(k)=-\frac{i}{2}\left(\beta-\sqrt{\beta^{2}-4\gamma}\right),~{}{{\omega_{2}}}(k)=-\frac{i}{2}\left(\beta+\sqrt{\beta^{2}-4\gamma}\right)$
(27a)

and

${{C}_{1}}(k)=\frac{\omega_{2}h_{0}-iu_{0}}{\omega_{2}-\omega_{1}},~{}{{C}_{2}}(k)=\frac{iu_{0}-\omega_{1}h_{0}}{\omega_{2}-\omega_{1}}$
(27b)

The final solution for $h(x,t)$ may be obtained by utilizing the inverse
Fourier Transform (8) to obtain

$h\left(x,t\right)=\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}\left[{C}_{1}(k){{e}^{{{\phi}_{1}}(k)t}}+{C}_{2}(k){{e}^{{{\phi}_{2}}(k)t}}\right]dk$
(27c)

where

${{\phi}_{1}}(k)=i\left[k\frac{x}{t}-{{\omega}_{1}}(k)\right],~{}{{\phi}_{2}}(k)=i\left[k\frac{x}{t}-{{\omega}_{2}}(k)\right].$
(27d)

The result (27d) provides the integral solution to the system (18), where
$\beta$ and $\gamma$ are given in (21b).

The long time asymptotic behavior of (27d) for $R\\!e$$\leq\cot(\theta)$ may
be used to establish the stability of the flow system governed by equation
(18). In Appendix B, it is shown that the second term in (27c) damps faster
than the first as $t\to\infty$ for all $k$, and so we may recast (27c) as

$h\left(x,t\right)\sim\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}{C}_{1}(k){{e}^{{{\phi}_{1}}(k)t}}dk,~{}t\to\infty,~{}R\\!e\leq\cot(\theta)$
(28)

Although we cannot determine a closed-form solution for the Fourier integral
in (28), it can be evaluated via the method of steepest descent [Bender &
Orszag, 1999], where we allow $k=k_{r}+ik_{i}$ to be complex and look for a
closed integration path (in the complex $k$-plane) that includes the real line
and passes through a saddle point of (27d); this enables the use of Cauchy’s
theorem. An $n^{\textrm{th}}$ order saddle $k_{s}$ of $\phi_{1}$ is defined by

$\left.\frac{d\phi_{1}}{dk}\right|_{k_{s}}=\dots=\left.\frac{d^{n-1}\phi_{1}}{dk^{n-1}}\right|_{k_{s}}=0,~{}\left.\frac{d^{n}\phi_{1}}{dk^{n}}\right|_{k_{s}}\neq
0,~{}n>1.$ (29)

Note that, depending on the direction of approach, saddle points can describe
where Real[$\phi_{1}]$ reaches a maximum or minimum. The deformed integration
path of steepest descent moves through the saddle such that the imaginary part
of $\phi_{1}(k)$ is constant and real part of $\phi_{1}(k)$ attains its
maximum value at the saddle. This enables an asymptotic expansion about the
saddle point, which can be used to deduce the dominant behavior as $t$ becomes
large. The path of ascent, where a minimum of Real[$\phi_{1}]$ is reached at
the saddle, is not useful in the long-time evaluation of Fourier integrals
such as (28). Note that (28) contains branch points when $\beta^{2}=4\gamma$
and so care must be taken such that the aforementioned path does not enclose
such points, as this would violate Cauchy’s theorem.

Substituting (27d) into (29) leads to the relations for
$\omega(k)=\omega_{r}(k)+i\omega_{i}(k)$

$\left.\frac{\partial\omega_{r}}{dk_{r}}\right|_{k_{s}}=\frac{x}{t}$ (30a)
$\left.\frac{\partial\omega_{i}}{dk_{r}}\right|_{k_{s}}=0.$ (30b)

Equation (30b) highlights that each saddle point $k_{s}$ is paired with an
$x/t$ ray; it provides a simultaneous set of equations to solve for the real
and imaginary parts of $k_{s}$ for a given $x/t$. We may deduce the long-time
behavior of (28) along specific $x/t$ values by expanding $\phi_{1}(k)$ about
the corresponding saddles. Using this expansion within the method of steepest
descent outlined above (see Appendix C), we arrive at the following long-time
behavior for (28)

$h\left(x,t\right)|_{\frac{x}{t}=3}\sim\frac{h_{0}+u_{0}~{}R\\!e/3}{2\sqrt{\pi\left[\cot(\theta)-R\\!e\right]}}~{}t^{-\frac{1}{2}},~{}~{}R\\!e<\cot(\theta)$
(31a)
$h\left(x,t\right)|_{\frac{x}{t}=3}\sim\frac{1}{4\pi}\Gamma\left(\frac{1}{4}\right)\left[\frac{3~{}W\\!e}{R\\!e}\right]^{\frac{1}{4}}\left(h_{0}+u_{0}~{}R\\!e/3\right)~{}t^{-\frac{1}{4}},~{}~{}R\\!e=\cot(\theta)$
(31b) $h\left(x,t\right)|_{\frac{x}{t}\neq 3}\sim
c_{1}t^{-\frac{1}{2}}e^{c_{2}t}\cos\left(c_{3}t+c_{4}\right),~{}~{}c_{2}<0,~{}R\\!e\leq\cot(\theta)$
(31c)

where the real-valued parameters $c_{1}$, $c_{2}$, $c_{3}$, and $c_{4}$ are
functions of $x/t$; a more specific form of (31c) for direct use is given by
(C.2) in Appendix C. Note that there is a structural change for $x/t=3$ and
$x/t\neq 3$ (for all $R\\!e$$\leq\cot(\theta)$). For any $R\\!e$, $x/t=3$ is
the least damped ray, only exhibiting algebraic decay, while all other rays
damp exponentially at $t\to\infty$. There is also a structural change for
$R\\!e$ $<\cot(\theta)$ and $R\\!e$$=\cot(\theta)$ along the ray $x/t=3$. For
$R\\!e$$<\cot(\theta)$, the peak of the wave packet (at $x/t=3$) decays like
$t^{-1/2}$. For $R\\!e$$=\cot(\theta)$, the peak of the wave packet decays
like $t^{-1/4}$. This is a ramification of the saddle point being
2${}^{\textrm{nd}}$-order for the former and 4${}^{\textrm{th}}$-order for the
latter, as shown in Appendix C; this change of order can also be seen in
(LABEL:eq:17), setting $R\\!e$$=\cot(\theta)$. Here, we have deduced the
stability of the system (18) for $R\\!e$$\leq\cot(\theta)$, in a parameter
space where classical stability analysis is inconclusive.

The Fourier series solution to (18) is shown in Fig. 7, where the peak is seen
to decay in accordance with the asymptotic solution given by (31b) along
$x/t=3$, as indicated by a black line in the figure. Away from the peak (for
$x/t\neq 3$) the series solution approaches the asymptotic solution given by
(31c), as indicated in Fig. 8.

(a)  Figure 7: Fourier series solution (shaded surface) to (18) compared with
(31b) (black curve) along the line $x/t=3$. $R\\!e$=$\cot(\theta)$,
$h_{0}=u_{0}=1$, $\theta=\pi/4$, and $W\\!e$=0.1. (a)  Figure 8: Fourier
series solution (solid curves) to (18) along specific $x/t$ rays, compared
with the asymptotic solution given by (31c) (dashed curves) for
$R\\!e$=$\cot(\theta)$, $h_{0}=u_{0}=1$, $\theta=\pi/4$, and $W\\!e$=0.1.

##  4 Conclusions

Flow transition from stability to instability is examined in which a single
mode (i.e., for one value of wavenumber) indicates classical neutral
stability. It is found that fluid flows characterized as being neutrally
stable via exponential modes in fact may exhibit either algebraic growth or
decay. Two partial differential equation systems are examined, one that has
been constructed to elucidate key features of the stability threshold, and a
second that models the well-studied problem of rectilinear Newtonian flow down
an inclined plane. In the former case, algebraic growth occurs on the neutral
stability threshold. In the latter case, algebraic decay occurs both at and
below the critical Reynolds number.

A key difference between the algebraic growth problem presented here and the
one in [King et al., 2016] is that here only one wave-number touches the real-
axis on a classical stability plot, whereas all wavenumbers are neutral in the
prior work. A ramification of this appears to be that the non-growing rays
(away from the peak of the wavepacket) decay exponentially instead of
algebraically. Whether or not this is true for all such problems remains to be
explored. In both problems, algebraic instability arises from removable
singularities in the Fourier integral solution of the governing PDE.

For the second problem reported here – that of incline plane flow –
exponential growth occurs above a critical Reynolds number. In his own work,
[Yih, 1963] also reached the conclusion that the flow is unstable above the
critical value, but could not reach a conclusion for when the Reynolds number
was below the critical value. The discovery of algebraic stability below the
critical Reynolds number carries with it the danger that the flow will not be
as stable as it would be if all values of $\omega_{i}$ were negative and
disturbances damped exponentially. In the context of coating operations that
utilize an inclined plane-flow geometry to form liquid films, our results show
that disturbances to a process, even if damped, will damp out more slowly than
expected in the physical domain. It is possible that, given larger
disturbances, even a damped response may lead to unacceptable perturbations in
coated liquid films.

## References

  * [Abramowitz & Stegun, 1972] Abramowitz, M. & Stegun, I. (1972) Handbook of Mathematical Functions, page 298. Dover.
  * [Barlow et al., 2011] Barlow, N. S., Helenbrook, B. T. & Lin, S. P. (2011) Transience to instability in a liquid sheet. J. Fluid. Mech., 666, 358–390.
  * [Barlow et al., 2010] Barlow, N. S., Helenbrook, B. T., Lin, S. P. & Weinstein, S. J. (2010) An interpretation of absolutely and convectively unstable waves using series solutions. Wave Motion, 47(8), 564–582.
  * [Barlow et al., 2017] Barlow, N. S., Helenbrook, B. T. & Weinstein, S. J. (2017) Algorithm for spatio-temporal analysis of the signaling problem. IMA J. Appl. Math., 82, 1–32.
  * [Bender & Orszag, 1999] Bender, C. M. & Orszag, S. A. (1999) Advanced Mathematical Methods for Scientists and Engineers I. Springer Science & Business Media.
  * [Bleistein, 1984] Bleistein, N. (1984) Mathematical methods for wave phenomena. Academic Press.
  * [Case, 1960] Case, K. (1960) Stability of inviscid plane Couette flow. Physics of Fluids (1958-1988), 3(2), 143–148.
  * [Chandrasekhar, 1968] Chandrasekhar, S. (1968) Hydrodynamic and Hydromagnetic Stability. Clarendon Press.
  * [Cohen & Gutoff, 1992] Cohen, E. & Gutoff, E. (1992) Modern Coating and Drying Technology. Wiley.
  * [De Luca & Costa, 1997] De Luca, L. & Costa, M. (1997) Instability of a spatially developing liquid sheet. Journal of Fluid Mechanics, 331, 127–144.
  * [Drazin & Reid, 2004] Drazin, P. G. & Reid, W. H. (2004) Hydrodynamic stability. Cambridge university press.
  * [El-Sayed et al., 2015] El-Sayed, M. F., Moatimid, G., Elsabaa, F. & Amer, M. (2015) Hydrodynamic instabilities of two viscoelastic liquid sheet models in an inviscid gas medium. Atomization and Sprays, 25(2).
  * [Gradshteyn & Ryzhik, 2007] Gradshteyn, I. S. & Ryzhik, I. M. (2007) Table of integrals, series, and products: seventh edition. Academic Press.
  * [Huerre, 2000] Huerre, P. (2000) Perspectives in Fluid Dynamics, chapter Open shear flow instabilities, pages 159–229. Cambridge University Press.
  * [Huerre & Rossi, 1998] Huerre, P. & Rossi, M. (1998) Hydrodynamic Instabilities in open flows. In Godrèche, C. & Manneville, P., editors, Hydrodynamics and Nonlinear Instabilities, pages 81–288. Cambridge University Press.
  * [Ibrahim, 1995] Ibrahim, E. (1995) Spatial instability of a viscous liquid sheet. Journal of Propulsion and Power, 11(1), 146–152.
  * [Joulin & Vidal, 1998] Joulin, G. & Vidal, P. (1998) An introduction to the instability of flames, shocks, and detonations. In Godrèche, C. & Manneville, P., editors, Hydrodynamics and Nonlinear Instabilities, pages 493–673. Cambridge University Press.
  * [King et al., 2016] King, K. R., Weinstein, S. J., Zaretzky, P. M., Cromer, M. & Barlow, N. S. (2016) Stability of algebraically unstable dispersive flows. Phys. Rev. Fluids, 1(073604).
  * [Kistler & Schweizer, 1997] Kistler, S. F. & Schweizer, P. M. (1997) Liquid film coating. Springer.
  * [L. Brevdo & Bridges, 1999] L. Brevdo, P. Laure, F. D. & Bridges, T. J. (1999) Linear Pulse Structure and Signaling in a Film Flow on an Inclined Plane. J. Fluid. Mech., 396, 37–71.
  * [Lighthill, 2001] Lighthill, J. (2001) Waves in fluids. Cambridge university press.
  * [Lin, 2003] Lin, S. (2003) Breakup of liquid sheets and jets. Cambridge University Press.
  * [Manneville, 1998] Manneville, P. (1998) Overview. In Godrèche, C. & Manneville, P., editors, Hydrodynamics and Nonlinear Instabilities, pages 1–24. Cambridge University Press.
  * [Paul et al., 2004] Paul, E. L., Atiemo-Obeng, V. A. & Kresta, S. M. (2004) Handbook of industrial mixing: science and practice. John Wiley & Sons.
  * [Pohlhausen, 1921] Pohlhausen, K. (1921) Zur näherungsweisen Integration der Differentialgleichung der laminaren Reibungsshicht. ZAMM, 1(252-268).
  * [Pritchard, 2011] Pritchard, P. J. (2011) Fox and McDonald’s Introduction to Fluid Mechanics. Wiley, 8 edition.
  * [Rayeigh, 1880] Rayeigh, L. (1880) On the stability of certain fluid motions. Proc. Math. Soc. Lond., 11, 57–70.
  * [S. V. Alekseenko & Pokusaev, 1994] S. V. Alekseenko, V. E. N. & Pokusaev, B. G. (1994) Wave flow of liquid films. Begell House.
  * [Schlichting, 1979] Schlichting, H. (1979) Boundary Layer Theory. McGraw-Hill, seventh edition.
  * [Schmid & Henningson, 2001] Schmid, P. J. & Henningson, D. S. (2001) Stability and transition in shear flows. Springer Science & Business Media.
  * [von Kármán, 1921] von Kármán, T. (1921) Über laminare und turbulente Reibung. ZAMM, 1, 233–252.
  * [Weinstein, 1990] Weinstein, S. J. (1990) Wave Propagation In The Flow Of Shear-Thinning Fluids Down An Incline. AIChE J., 36(12), 1873–1889.
  * [Weinstein & Ruschak, 2004] Weinstein, S. J. & Ruschak, K. J. (2004) Coating flows. Annu. Rev. Fluid Mech., 36, 29–53.
  * [Whitham, 2011] Whitham, G. B. (2011) Linear and nonlinear waves, volume 42. John Wiley & Sons.
  * [Yih, 1963] Yih, C.-S. (1963) Stability of Liquid Flow down an Inclined Plane. Phys. Fluids, 6(3), 321–334.

##  Appendix A Evaluation of integrals in Eq. (LABEL:eq:trig)

###  A.1

In this section, general formulas for convergent integrals of the following
form are established:

$I=\int_{0}^{\infty}\frac{e^{-\beta k^{2}}}{k^{2}}\sin(\alpha
k^{2})\cos\left(bk\right)~{}dk,~{}~{}\beta>0,~{}~{}\alpha>0,~{}~{}b>0.$ (32)

First, we introduce the additional variable $\gamma$ in (32) as follows

$I=\int_{0}^{\infty}\frac{e^{-\gamma\beta k^{2}}}{k^{2}}\sin(\gamma\alpha
k^{2})\cos\left(bk\right)~{}dk,$ (33)

where $\gamma=1$ in (33) yields (32). Next, (33) is differentiated with
respect to $\gamma$ to obtain

$\frac{dI}{d\gamma}=-\beta\int_{0}^{\infty}e^{-\gamma\beta
k^{2}}\sin(\gamma\alpha
k^{2})\cos\left(bk\right)~{}dk+\alpha\int_{0}^{\infty}e^{-\gamma\beta
k^{2}}\cos(\gamma\alpha k^{2})\cos\left(bk\right)~{}dk.$ (34)

In order to recover the general solution for $I$, we note that (33) indicates
that the following constraint must be satisfied:

$I=0~{}\textrm{ at }~{}\gamma=0.$ (35)

The integrals in (34) are provided in closed-form on pages 493-494, Eqs.
(3.922-3) and (3.922-4) of [Gradshteyn & Ryzhik, 2007] (for
$\tilde{\beta}>0,~{}a>0,~{}b>0$):

$\displaystyle\int_{0}^{\infty}e^{-\tilde{\beta}k^{2}}\sin(ak^{2})\cos\left(bk\right)~{}dk$
$\displaystyle=-\sqrt{\frac{\pi}{8\left(\tilde{\beta}^{2}+a^{2}\right)}}\textrm{
exp
}\left[-\frac{b^{2}\tilde{\beta}}{4\left(\tilde{\beta}^{2}+a^{2}\right)}\right]$
(36a)
$\displaystyle\times\left\\{\sqrt{\sqrt{\tilde{\beta}^{2}+a^{2}}+\tilde{\beta}}~{}\sin\left[\frac{b^{2}a}{4\left(\tilde{\beta}^{2}+a^{2}\right)}\right]\right.$
$\displaystyle\left.-\sqrt{\sqrt{\tilde{\beta}^{2}+a^{2}}-\tilde{\beta}}~{}\cos\left[\frac{b^{2}a}{4\left(\tilde{\beta}^{2}+a^{2}\right)}\right]\right\\}$

$\displaystyle\int_{0}^{\infty}e^{-\tilde{\beta}k^{2}}\cos(ak^{2})\cos\left(bk\right)~{}dk$
$\displaystyle=\sqrt{\frac{\pi}{8\left(\tilde{\beta}^{2}+a^{2}\right)}}\textrm{
exp
}\left[-\frac{b^{2}\tilde{\beta}}{4\left(\tilde{\beta}^{2}+a^{2}\right)}\right]$
(36b)
$\displaystyle\times\left\\{\sqrt{\sqrt{\tilde{\beta}^{2}+a^{2}}+\tilde{\beta}}~{}\cos\left[\frac{b^{2}a}{4\left(\tilde{\beta}^{2}+a^{2}\right)}\right]\right.$
$\displaystyle\left.+\sqrt{\sqrt{\tilde{\beta}^{2}+a^{2}}-\tilde{\beta}}~{}\sin\left[\frac{b^{2}a}{4\left(\tilde{\beta}^{2}+a^{2}\right)}\right]\right\\}.$

Applying (36a) and (36b) to (34), and then integrating (34) from $\gamma=0$ to
1, we obtain

$\displaystyle
I=\sqrt{\frac{\pi}{8\left(\beta^{2}+\alpha^{2}\right)}}\int_{0}^{1}\frac{e^{-\beta\theta/\gamma}}{\sqrt{\gamma}}\left[C_{+}\cos\left(\alpha\theta/\gamma\right)+C_{-}\sin\left(\alpha\theta/\gamma\right)\right]~{}d\gamma$
$\displaystyle
C_{\pm}=\alpha\sqrt{\sqrt{\beta^{2}+\alpha^{2}}\pm\beta}\mp\beta\sqrt{\sqrt{\beta^{2}+\alpha^{2}}-\beta},~{}~{}\theta=\frac{b^{2}}{4\left(\beta^{2}+\alpha^{2}\right)}.$
(37)

Upon making the substitution $u=\beta\theta/\gamma$, (37) becomes

$I=\frac{b\sqrt{\beta\pi/2}}{4\left(\beta^{2}+\alpha^{2}\right)}\int_{\beta\theta}^{\infty}u^{-3/2}e^{-u}\left[C_{+}\cos\left(u\alpha/\beta\right)+C_{-}\sin\left(u\alpha/\beta\right)\right]~{}du,$
(38)

whose exact solution is given by Eqs. (3.944-2) and (3.944-4) on page 498 of
[Gradshteyn & Ryzhik, 2007], thus providing an exact solution for (32):

$\displaystyle I=\int_{0}^{\infty}$ $\displaystyle\frac{e^{-\beta
k^{2}}}{k^{2}}\sin(\alpha k^{2})\cos\left(bk\right)~{}dk=$
$\displaystyle\frac{b\sqrt{\pi/2}}{8\left(\beta^{2}+\alpha^{2}\right)}\left\\{\left(C_{+}+iC_{-}\right)\sqrt{\beta+i\alpha}~{}~{}\Gamma\left[-\frac{1}{2},\left(\beta+i\alpha\right)\frac{b^{2}}{4\left(\beta^{2}+\alpha^{2}\right)}\right]\right.$
$\displaystyle\hskip
65.04256pt+\left.\left(C_{+}-iC_{-}\right)\sqrt{\beta-i\alpha}~{}~{}\Gamma\left[-\frac{1}{2},\left(\beta-i\alpha\right)\frac{b^{2}}{4\left(\beta^{2}+\alpha^{2}\right)}\right]\right\\},$

where $\Gamma$ is the upper incomplete gamma function. Note that
(LABEL:eq:Igeneral) may be evaluated in the limit as $b\to 0$ to obtain

$\int_{0}^{\infty}\frac{e^{-\beta k^{2}}}{k^{2}}\sin(\alpha
k^{2})~{}dk=\sqrt{\frac{\pi}{2}}\frac{\alpha}{\sqrt{\sqrt{\alpha^{2}+\beta^{2}}+\beta}},~{}~{}\beta>0,~{}~{}\alpha>0.$
(40)

###  A.2

In this section, general formulas for convergent integrals of the following
form are established:

$I=\int_{-\infty}^{\infty}\frac{e^{-ak^{2}}-e^{-ck^{2}}}{k}e^{bk}~{}dk,~{}~{}a>0,~{}c>0.$
(41)

Although we follow a similar procedure as in Appendix A.1, the variable $b$
has no restriction on its sign, and thus may be used as a differentiation
variable in lieu of introducing a new “dummy” variable for this purpose.
First, (41) is differentiated with respect to $b$ to obtain

$\frac{dI}{db}=\int_{-\infty}^{\infty}\left(e^{-ak^{2}}-e^{-ck^{2}}\right)e^{bk}~{}dk.$
(42)

Equation (41) indicates that the following constraint must be satisfied:

$I=0~{}\textrm{ at }~{}b=0,$ (43)

since the integrand of (41) is odd and $k=0$ is a removable singularity. By
defining

$J=\int_{-\infty}^{\infty}e^{-ak^{2}}e^{bk}~{}dk.$ (44)

we can solve the first piece of (42) and then combine with the second piece
afterwards. The exponentials in (44) are rewritten by completing the square,
leading to

$J=e^{\frac{b^{2}}{4a}}\int_{-\infty}^{\infty}e^{-a\left(k-\frac{b}{2a}\right)^{2}}~{}dk,$
(45)

and the substitution $u=\sqrt{2a}\left(k-\frac{b}{2a}\right)$ is made to
obtain

$J=\frac{e^{\frac{b^{2}}{4a}}}{\sqrt{2a}}\int_{-\infty}^{\infty}e^{-\frac{1}{2}u^{2}}~{}du,$
(46)

where the integral in (46) is exactly equal to $\sqrt{2\pi}$, as defined by
the Gamma function. We now have the expression

$J=\sqrt{\frac{\pi}{a}}e^{\frac{b^{2}}{4a}},$

and, after repeating the above analysis for the second piece of (42), we
obtain

$\frac{dI}{db}=\sqrt{\frac{\pi}{a}}\left(e^{\frac{b^{2}}{4a}}-e^{\frac{b^{2}}{4c}}\right),$

which may be integrated with respect to $b$ (applying condition (43)) to
obtain

$I=\sqrt{\frac{\pi}{a}}\left[\int_{0}^{b}e^{\frac{b^{2}}{4a}}~{}db-\int_{0}^{b}e^{\frac{b^{2}}{4c}}~{}db\right].$
(47)

After making the variable substitutions $v=\frac{b}{2i\sqrt{a}}$ and
$v=\frac{b}{2i\sqrt{c}}$ in the first and second integrals respectively, (47)
becomes

$I=2\pi
i\left[\int_{0}^{\frac{b}{2i\sqrt{a}}}e^{-v^{2}}~{}dv-\int_{0}^{\frac{b}{2i\sqrt{c}}}e^{-v^{2}}~{}dv\right].$
(48)

The integrals in (48) are error functions, thus providing an exact solution
for (41):

$\displaystyle
I=\int_{-\infty}^{\infty}\frac{e^{-ak^{2}}-e^{-ck^{2}}}{k}e^{bk}~{}dk=\pi
i\left[~{}\textrm{erf}\left(\frac{b}{2i\sqrt{a}}\right)-\textrm{erf}\left(\frac{b}{2i\sqrt{c}}\right)\right],~{}~{}a>0,~{}c>0~{}$
(49)

##  Appendix B Justification for Asymptotic Equivalence in (28)

We now justify that the 2${}^{\textrm{nd}}$ term in (27c) is subdominant to
the first, which leads to the asymptotic equivalence shown in (28) as
$t\to\infty$. The typical approach to follow is to rewrite the integral (27c)
in two pieces as:

$h(x,t)={{h}_{1}}(x,t)+{{h}_{2}}(x,t)$ (50a)

where

${{h}_{1}}(x,t)=\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}{{{C}_{1}}}(k){{e}^{{{\phi}_{1}}(k)t}}dk,~{}~{}{{h}_{2}}(x,t)=\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}{{{C}_{2}}}(k){{e}^{{{\phi}_{2}}(k)t}}dk$
(50b)

In (50a), the path of integration is taken along the real axis as written. It
suffices, then to show that ${{h}_{2}}\ll{{h}_{1}}$ as $t\to\infty$ to justify
equation (28) . Interestingly, the integral ${{h}_{2}}$ cannot be evaluated
via standard techniques in the $t\to\infty$ limit. In particular the two
typical methods at our disposal – the method of steepest descent and
integration by parts – fail as follows.

When the method of steepest descent is used, the integral ${{h}_{2}}$ along
the real $k$-axis is evaluated as part of a complex contour integral, and the
topology of the complex phase function ${{\phi}_{2}}(k)$ is configured such
that saddle points, $k_{s}$, have lower
$\textrm{Real}[{{\phi}_{2}}({{k}_{s}})]$ values than surrounding regions of
the plane. Fig. 9 shows a typical case that represents the situation for all
$x/t$ rays evaluated in this paper. The indicated saddle points cannot be
accessed via a contour integral that includes the real axis (or even a portion
of it denoted as the contour C in Fig. 9) in such a way that the saddle points
have a maximum Real$[{{\phi}_{2}}({{k}_{s}})]$ along the deformed contour, and
the method fails.

Figure 9: Contours of constant $\textrm{Real}[{{\phi}_{2}}(k)]$ vs complex
$k={{k}_{r}}+i{{k}_{i}}$ at neutral stability ($R\\!e=\cot\theta$) for
$x/t=3$, $\theta=\pi/4$, and $W\\!e$=0.1. There are 3 saddle points ($\circ$)
denoted as $k_{s,j}$ for $j\in[1,3]$. Branch points ($\ast$) and branch cuts
are taken in accordance with the principle values of the square roots in
(27a). The directions in which $\textrm{Real}[{{\phi}_{2}}(k)]$ decrease and
increase are indicated with arrows using the notation D and U, respectively.
The orientation of the saddles does not enable a closed contour to be drawn
that includes the contour C (the real axis) and also has
$\textrm{Real}[{{\phi}_{2}}(k_{s,j})]$ as a maximum along the contour. Thus,
all the saddle points are spurious, i.e., none may be used via the method of
steepest decent to determine the long time asymptotic behavior of the
2$\mathrm{nd}$ term in the integral of (27c). For reference, the location of
the saddles and associated contour values are as follows:
$k_{s,1}=-0.2957-0.2492i$, $\textrm{Real}[{{\phi}_{2}}(k_{s,1})]=-2.4192$;
$k_{s,2}=0.2957-0.2492i$, $\textrm{Real}[{{\phi}_{2}}(k_{s,2})]=-2.4192$;
$k_{s,3}=0.6976i$, $\textrm{Real}[{{\phi}_{2}}(k_{s,3})]=-4.3440$.

The method of integration by parts fails precisely because the integral
${{h}_{2}}(x,t)$ is convergent in the infinite domain of $k$. This requires
its integrand go to zero as $k\to\pm\infty$ (and it does, based on the
contours in Fig. 9), and as a result repeated integration by parts merely
yields a zero result and no asymptotic behavior can be extracted.

A different approach is thus necessary to prove the assertion that
${{h}_{2}}(x,t)$ is subdominant. To do so, we return to the form of the
integral in (27c) and compare the magnitude of the terms in the integrand
directly at each value of $k$ (recall that $k$ is real as the path of
integration lies along the real axis). We denote the two pieces of the
integrand as:

${{I}_{1}}={{C}_{1}}(k){{e}^{{{\phi}_{1}}(k)t}},~{}~{}{{I}_{2}}={{C}_{2}}(k){{e}^{{{\phi}_{2}}(k)t}},$
(51a)

and thus using equations (27a), (27b), (27d), we obtain

$\displaystyle\frac{{{I}_{2}}}{{{I}_{1}}}=\frac{{{C}_{2}}(k)}{{{C}_{1}}(k)}{{e}^{({{\phi}_{2}}(k)-{{\phi}_{1}}(k))t}},~{}~{}\frac{{{C}_{2}}(k)}{{{C}_{1}}(k)}=\left(\frac{i{{u}_{0}}-{{\omega}_{1}}{{h}_{0}}}{{{\omega}_{2}}{{h}_{0}}-i{{u}_{0}}}\right),~{}~{}{{\phi}_{2}}(k)-{{\phi}_{1}}(k)=-\sqrt{{{\beta}^{2}}-4\gamma}$
(51b)

where all parameters are defined in equations (21b) and (27d). As $t\to\infty$
the real part of the exponential in equation (51b) governs the magnitude of
the ratio $I_{2}/I_{1}$; a typical plot of the growth rate in the exponential,
Real$[{{\phi}_{2}}(k)-{{\phi}_{1}}(k)]$ vs. $k$ is given in Fig. 10.

Figure 10: Plot of Real$[{{\phi}_{2}}(k)-{{\phi}_{1}}(k)]$ vs. $k$ for
$R\\!e=\cot(\theta)$, $\theta=\pi/4$, and $W\\!e$=0.1.

Fig. 10 indicates that as $t\to\infty$, the magnitude of
${{I}_{2}}/{{I}_{1}}\to 0$ for finite $k$. However, as the domain of $k$ is
increased toward $\pm\infty$ (not shown here), the plot asymptotes to zero,
indicating that the growth rates are comparable there. Furthermore, the
magnitude of the ratio ${{C}_{2}}(k)/{{C}_{1}}(k)$ in (51b) approaches 1 in
these limits. Thus, for all finite values of $k$, we can certainly establish
that ${{I}_{2}}$ is subdominant to ${{I}_{1}}$ as $t\to\infty$, but this
result is not proven in the limit of infinite $k$. This also indicates that in
(50b), ${{h}_{2}}\ll{{h}_{1}}$ as $t\to\infty$ if finite bounds are taken on
the indicated integrals instead of the infinite limits indicated.

To complete the proof that ${{h}_{2}}\ll{{h}_{1}}$ as $t\to\infty$, we rewrite
the equations for ${{h}_{1}}$ and ${{h}_{2}}$ in (50b) as

${{h}_{1}}(x,t)=\frac{1}{2\pi}\int\limits_{-\infty}^{-L}{{{C}_{1}}}(k){{e}^{{{\phi}_{1}}(k)t}}dk+\frac{1}{2\pi}\int\limits_{-L}^{L}{{{C}_{1}}}(k){{e}^{{{\phi}_{1}}(k)t}}dk+\frac{1}{2\pi}\int\limits_{L}^{\infty}{{{C}_{1}}}(k){{e}^{{{\phi}_{1}}(k)t}}dk~{}~{}$
(52a)
${{h}_{2}}(x,t)=\frac{1}{2\pi}\int\limits_{-\infty}^{-L}{{{C}_{2}}}(k){{e}^{{{\phi}_{2}}(k)t}}dk+\frac{1}{2\pi}\int\limits_{-L}^{L}{{{C}_{2}}}(k){{e}^{{{\phi}_{2}}(k)t}}dk+\frac{1}{2\pi}\int\limits_{L}^{\infty}{{{C}_{2}}}(k){{e}^{{{\phi}_{2}}(k)t}}dk~{}~{}$
(52b)

where $L>0$ is a finite number. We have already established that the finite
bound integral in (52b) is subdominant to the corresponding finite bound
integral in (52a) from Fig. 10 and the preceding arguments (in Fig. 10, $L$=10
as indicated). Integration by parts may be used on the remaining semi-infinite
integrals in equation (52a) and (52b) to yield

$\displaystyle\frac{1}{2\pi}\int\limits_{-\infty}^{-L}{{{C}_{1}}}(k){{e}^{{{\phi}_{1}}(k)t}}dk\sim\frac{{{C}_{1}}(-L)}{{2\pi
t~{}{\left.\frac{d{{\phi}_{1}}(k)}{dk}\right|}_{k=-L}}}{{e}^{{{\phi}_{1}}(-L)t}},$
$\displaystyle\frac{1}{2\pi}\int\limits_{L}^{\infty}{{{C}_{1}}}(k){{e}^{{{\phi}_{1}}(k)t}}dk\sim-\frac{{{C}_{1}}(L)}{2\pi
t~{}{{\left.\frac{d{{\phi}_{1}}(k)}{dk}\right|}_{k=L}}}{{e}^{{{\phi}_{1}}(L)t}},$
(53a)
$\displaystyle\frac{1}{2\pi}\int\limits_{-\infty}^{-L}{{{C}_{2}}}(k){{e}^{{{\phi}_{2}}(k)t}}dk\sim\frac{{{C}_{2}}(-L)}{{2\pi
t~{}{\left.\frac{d{{\phi}_{2}}(k)}{dk}\right|}_{k=-L}}}{{e}^{{{\phi}_{2}}(-L)t}},$
$\displaystyle\frac{1}{2\pi}\int\limits_{L}^{\infty}{{{C}_{2}}}(k){{e}^{{{\phi}_{2}}(k)t}}dk\sim-\frac{{{C}_{2}}(L)}{{2\pi
t~{}{\left.\frac{d{{\phi}_{2}}(k)}{dk}\right|}_{k=L}}}{{e}^{{{\phi}_{2}}(L)t}}\textrm{
as }t\to\infty.$ (53b)

We thus see that

$\displaystyle\frac{\frac{1}{2\pi}\int\limits_{-\infty}^{-L}{{{C}_{2}}}(k){{e}^{{{\phi}_{2}}(k)t}}dk}{\frac{1}{2\pi}\int\limits_{-\infty}^{-L}{{{C}_{1}}}(k){{e}^{{{\phi}_{1}}(k)t}}dk}=O({{e}^{({{\phi}_{2}}(-L)-{{\phi}_{1}}(-L))t}}),$
$\displaystyle\frac{\frac{1}{2\pi}\int\limits_{L}^{\infty}{{{C}_{2}}}(k){{e}^{{{\phi}_{2}}(k)t}}dk}{\frac{1}{2\pi}\int\limits_{L}^{\infty}{{{C}_{1}}}(k){{e}^{{{\phi}_{1}}(k)t}}dk}=O({{e}^{({{\phi}_{2}}(L)-{{\phi}_{1}}(L))t}})\textrm{
as }t\to\infty.$ (54)

Since according to Fig. 10, which is representative for all cases examined in
this paper, Real$[{{\phi}_{2}}(k)-{{\phi}_{1}}(k)]<0$ for finite $k$, this
indicates that the ratios in (54) go to zero as $t\to\infty$. Thus, we see
that all integrals in (52b) are subdominant to those in (52a) for all real
$k$, which establishes that ${{h}_{2}}\ll{{h}_{1}}$ as $t\to\infty$ in (50b).
This furthermore establishes the asymptotic equivalence indicated in equation
(28) of the main text. This conclusion is also demonstrated by the agreement
between our numerical solutions and the $t\to\infty$ asymptotic behavior of
equation (28).

##  Appendix C Long time asymptotic solution to (28)

The following analysis is separated into two subsections based on a structural
change in the dispersion relation for flow down an incline plane. Classical
analysis tells us that, for $R\\!e$$\leq\cot(\theta)$, the wavenumber of
maximum growth is $k=0$ (see Fig. 6) and that the exponential growth rate is
zero. Note that, by definition (30bb), a maximum in $\omega_{i}(k_{r})$ is a
saddle point. Any such maximum is also a contributing saddle point, in that a
steepest descent path can be closed back to the real-axis (see [Barlow et al.,
2017], Appendix A). Applying definition (30ba) to (LABEL:eq:17), we find that
the corresponding ray of maximum growth is $x/t=3$. Thus, for
$R\\!e$$\leq\cot(\theta)$, the peak of a wave packet travels at a velocity
$x/t$=3 as it flows down the incline plane. The explicit knowledge of the
($k_{s}$, $x/t$) pair at the peak allows us to make simplifications outlined
above. For $x/t\neq 3$, the saddle locations and their corresponding steepest
descent paths are less straightforward to deduce. For this reason, we separate
the following evaluation of integral (28) into two subsections - one for the
peak along $x/t=3$ and one for the “off-peak” rays $x/t\neq 3$.

###  C.1 Evaluation of (28) for $\frac{x}{t}=3$ and $R\\!e\leq\cot(\theta)$

Here, we apply the method of steepest descent to the integral (28) for
$x/t$=3, which corresponds to the contributing saddle $k_{s}=0$. Expanding
$\phi_{1}(k)$ (given by (27d)) about $k=0$, and making direct use of
(LABEL:eq:17), we obtain the following

$\displaystyle\phi_{1}(k)\sim$
$\displaystyle\left[R\\!e-\cot(\theta)\right]k^{2}-\frac{6}{5}iR\\!e\left[R\\!e-\cot(\theta)\right]k^{3}$
$\displaystyle+\left\\{\frac{36}{25}R\\!e^{2}\left[R\\!e-\cot(\theta)\right]-\frac{12}{5}\textrm{Re}^{3}c\left[R\\!e-\cot(\theta)\right]^{2}-\frac{R\\!e}{3W\\!e}\right\\}k^{4}+O(k^{5}),~{}~{}k\to
0,$ (55)

and thus

$\phi_{1}(k)\sim\left[R\\!e-\cot(\theta)\right]k^{2},~{}~{}\textrm{ for
}R\\!e<\cot(\theta)$ (56a)
$\phi_{1}(k)\sim-\frac{\textrm{Re}}{3W\\!e}k^{4},~{}~{}\textrm{ for
}R\\!e=\cot(\theta),$ (56b)

which indicates that $k_{s}=0$ is a 2nd order saddle for
$R\\!e$$<\cot(\theta)$ and a 4th order saddle for $R\\!e$$=\cot(\theta)$,
according to definition (29). Substituting $k=0$ into $C_{1}(k)$ and (56b) for
$\phi_{1}(k)$ in (28) leads to

$\displaystyle
h\left(x,t\right)|_{\frac{x}{t}=3}\sim\frac{1}{2\pi}\left(h_{0}+u_{0}~{}R\\!e/3\right)\int\limits_{-\infty}^{\infty}e^{-\left[\cot(\theta)-R\\!e\right]k^{2}t}~{}dk,~{}~{}t\to\infty,~{}~{}\textrm{
for }R\\!e<\cot(\theta)$ $\displaystyle
h\left(x,t\right)|_{\frac{x}{t}=3}\sim\frac{1}{2\pi}\left(h_{0}+u_{0}~{}R\\!e/3\right)\int\limits_{-\infty}^{\infty}e^{-\frac{R\\!e}{3W\\!e}k^{4}t}~{}dk,~{}~{}t\to\infty,~{}~{}\textrm{
for }R\\!e=\cot(\theta)$ (57)

where the integration path remains along the real line, since the argument of
the exponential is purely real and thus no rotation through the saddle is
required. Since the integrands of (57) are even, the integrals may be
rewritten

$\displaystyle
h\left(x,t\right)|_{\frac{x}{t}=3}\sim\frac{1}{\pi}\left(h_{0}+u_{0}~{}R\\!e/3\right)\int\limits_{0}^{\infty}e^{-\left[\cot(\theta)-R\\!e\right]k^{2}t}~{}~{}dk,~{}t\to\infty,~{}~{}\textrm{
for }R\\!e<\cot(\theta)$ (58a) $\displaystyle
h\left(x,t\right)|_{\frac{x}{t}=3}\sim\frac{1}{\pi}\left(h_{0}+u_{0}~{}R\\!e/3\right)\int\limits_{0}^{\infty}e^{-\frac{R\\!e}{3W\\!e}k^{4}t}~{}dk,~{}~{}t\to\infty,~{}~{}\textrm{
for }=R\\!e\cot(\theta).$ (58b)

Upon making the variable substitutions
$v=\left[\cot(\theta)-R\\!e\right]k^{2}t$ and $v=\frac{R\\!e}{3W\\!e}k^{4}t$
in (58a) and (58b) respectively, we obtain

$\displaystyle
h\left(x,t\right)|_{\frac{x}{t}=3}\sim\frac{\left(h_{0}+u_{0}~{}R\\!e/3\right)}{\pi
p\sqrt{\left[\cot(\theta)-R\\!e\right]t}}\int\limits_{0}^{\infty}e^{-v}v^{\frac{1}{2}-1}~{}~{}dv,~{}t\to\infty,~{}~{}\textrm{
for }R\\!e<\cot(\theta)$ (59a) $\displaystyle
h\left(x,t\right)|_{\frac{x}{t}=3}\sim\frac{\left(h_{0}+u_{0}~{}R\\!e/3\right)}{\pi
p\left\\{\frac{\textrm{Re}}{3W\\!e}t\right\\}^{1/4}}\int\limits_{0}^{\infty}e^{-v}v^{\frac{1}{4}-1}~{}~{}dv,~{}t\to\infty,~{}~{}\textrm{
for }R\\!e=\cot(\theta)$ (59b)

where the integrals above evaluate to $\Gamma(1/2)$ and $\Gamma(1/4)$
respectively, leading the the results given in (31a) and (31b).

###  C.2 Evaluation of (28) for $\frac{x}{t}\neq 3$ and
$R\\!e\leq\cot(\theta)$

Saddles associated with the off-peak rays $x/t\neq 3$ lie in the complex
$k$-plane, off of the real-line. Thus we must take care to evaluate integral
(28) along a closed path that joins the original path (the $k_{r}$-axis) with
a path through the saddle, such that Real[$\phi(k_{s})$] is a maximum along
this path. This allows for the integral to be replaced with an approximation
near the saddle as $t\to\infty$. A representative example is given in Fig. 11
for $x/t=2$ (parameter values given in the caption).

Figure 11: Contours of constant $\textrm{Real}[{{\phi}_{1}}(k)]$ vs complex
$k={{k}_{r}}+i{{k}_{i}}$ at neutral stability ($R\\!e=\cot\theta$) for
$x/t=2$, $\theta=\pi/4$, and $W\\!e$=0.1. Notation is the same as in Fig. 9.
The orientation of the saddles enables a closed contour (indicated by dashed
lines) to be drawn through $k_{s,1}$ and $k_{s,2}$ that includes a portion of
the real axis, C1, and also has $\textrm{Real}[{{\phi}_{1}}(k_{s,j})]$ as a
maximum along the contour – thus, the method of steepest decent may be used to
determine the long time asymptotic behavior of the integral involving
$\phi_{1}(k)$ in (27c) along the contour C1. Alternatively, the orientation of
saddle $k_{s,3}$ makes it spurious (does not contribute to the asymptotic
behavior) and is not accessible by connection to any real contour. The path of
integration in (27c) is the infinite $k_{r}$-axis, i.e. consists of the
contour C=C${}_{\mathrm{L}}\cup$C${}_{1}\cup$CR; however, it can be shown via
integration by parts that the contributions over CL and CR are subdominant to
that of C1, such that the long time asymptotic behavior is determined by the
integral C1 – and thus the saddle points dominate the asymptotic behavior. For
reference, the location of the contributing saddles and related quantities
necessary to determine the asymptotic behavior using (62), (63), and (73a) are
as follows: $a$=0.2957, $b$=0.2492, $c$=0.8218, $d$=0.6805, $r$=5.9476, and
$q$=5.4227.

Note that there are 3 saddle points associated with $x/t$=2, as indicated by
$\bullet$’s in Fig. 11, but only the two lying in the lower half plane are
accessible as maxima along a closed path with the real axis. The steepest
descent contour $\int_{\textrm{SD}}$ is shown by a dashed line. Note that
neither the branch points (indicated by $\ast$’s) nor any other singularities
are enclosed between the real axis and the steepest descent path, and so we
may apply Cauchy’s theorem to obtain

$h\left(x,t\right)|_{\frac{x}{t}\neq
3}\sim\frac{1}{2\pi}\int_{\mathrm{-SD}}{C}_{1}(k){{e}^{{{\phi}_{1}}(k)t}}dk,~{}t\to\infty,~{}R\\!e\leq\cot(\theta)$
(60)

where $\int_{\textrm{-SD}}$ denotes a path from left to right (opposite that
shown in Fig. 11) and the long-time asymptotic casting in (60) is solely due
the omission of the $\phi_{2}(k)$ term in (27d), and not from the path
deformation. As Reynolds number, Weber number, $\theta$, and $x/t$ vary, the
location of the saddles differ from that shown in Fig. 11. However, there is
always either one or two contributing saddles where Cauchy’s theorem may be
applied such that (60) holds. Additionally, since all $x/t\neq 3$ saddles are
2nd-order, a general asymptotic solution may be formulated according to the
method of steepest descent by looking in the vicinity of the saddles,
replacing $C_{1}(k)$ with $C_{1}(k_{s})$, $\phi(k)$ with its 2nd-order
expansion about $k_{s}$ in (60), and rewriting the limits of integration as
follows

$\displaystyle h\left(x,t\right)|_{\frac{x}{t}\neq 3}\sim$
$\displaystyle\frac{1}{2\pi}\Bigg{[}{C}_{1}(k_{s1})e^{\phi_{1}(k_{s1})}\overbrace{\int_{k_{s1}-\epsilon}^{k_{s1}+\epsilon}e^{\frac{1}{2}\left.\frac{d^{2}\phi_{1}}{dk^{2}}\right|_{k_{s1}}(k-k_{s1})^{2}t}dk}^{I_{1}}$
$\displaystyle+~{}{C}_{1}(k_{s2})e^{\phi_{1}(k_{s2})}\underbrace{\int_{k_{s2}-\epsilon}^{k_{s2}+\epsilon}e^{\frac{1}{2}\left.\frac{d^{2}\phi_{1}}{dk^{2}}\right|_{k_{s2}}(k-k_{s2})^{2}t}dk}_{I_{2}}\Bigg{]},~{}t\to\infty,~{}R\\!e\leq\cot(\theta)$

where $\epsilon$ is a small positive constant. Note that in the above we have
assumed two contributing saddles $k_{s1}$ and $k_{s2}$. For cases with two
saddles, such as that shown in Fig. 11, the locations of quantities in
(LABEL:eq:Cauchy) in the complex $k$-plane are as follows

$\displaystyle k_{s1,2}=\mp a-ib,~{}\omega_{1}(k_{s1,2})=\mp c-id,$
$\displaystyle\left.\frac{d^{2}\phi_{1}}{dk^{2}}\right|_{k_{s1,2}}=-r\mp iq,$
$\displaystyle a,~{}b,~{}c,~{}d,~{}q,~{}r>0.$ (62)

It is useful to write the following quantity from (62) in complex polar form

$r+iq=\rho
e^{i\Theta},~{}\rho=\sqrt{r^{2}+q^{2}},~{}\Theta=\tan^{-1}\left(\frac{q}{r}\right),~{}\Theta\in\left(0,\frac{\pi}{2}\right)$
(63)

from which it follows from (62) that

$\left.\frac{d^{2}\phi_{1}}{dk^{2}}\right|_{k_{s1,2}}=-\rho e^{\pm i\Theta}.$
(64)

We now proceed to solve for $I_{1}$ (denoted in (LABEL:eq:Cauchy)) using the
orientation of the path through the saddle $k_{s1}$ shown in Fig. 11. After
substituting (64) into (62) and splitting $I_{1}$ into two integrals entering
and leaving the saddle $k_{s1}$, we obtain

$I_{1}=\int_{k_{s1}-\epsilon}^{k_{s1}}e^{-\frac{1}{2}\rho
e^{i\Theta}(k-k_{s1})^{2}t}dk+\int_{k_{s1}}^{k_{s1}+\epsilon}e^{-\frac{1}{2}\rho
e^{i\Theta}(k-k_{s1})^{2}t}dk.$ (65)

The split made above is done such that the quantity $(k-k_{s1})$ can be
written in complex polar form in accordance with Fig. 11 with the appropriate
angle entering and leaving the saddle as follows

$\displaystyle(k-k_{s1})=Ue^{i\Psi}$
$\displaystyle\Psi\in\left[\frac{\pi}{2},\pi\right]\textrm{ for
}k\in\left[k_{s1}-\epsilon,k_{s1}\right]$
$\displaystyle\Psi\in\left[-\frac{\pi}{2},0\right]\textrm{ for
}k\in\left[k_{s1},k_{s1}+\epsilon\right].$ (66)

Upon substituting (66) into (65), we obtain

$I_{1}=e^{i\Psi}\int_{-\epsilon e^{-i\Psi}}^{0}e^{-\frac{1}{2}\rho
U^{2}te^{i(\Theta+2\Psi)}}dU+e^{i\Psi}\int_{0}^{\epsilon
e^{-i\Psi}}e^{-\frac{1}{2}\rho U^{2}te^{i(\Theta+2\Psi)}}dU.$ (67)

Note that in order for the integrals of (67) to represent paths of steepest
descent, the argument of the outer exponential must be real and negative, such
that following condition must hold between $\Psi$ and $\Theta$ [Bleistein,
1984].

$\Theta+2\Psi=2n\pi$ (68)

where $n=1$ for the first integral of (67) (entering the saddle) and $n=0$ for
the second integral of (67) (leaving the saddle). In accordance with (68), we
substitute $\Psi=\pi-\Theta/2$ into the first integral of (67) and
$\Psi=-\Theta/2$ into the second integral of (67) to obtain (after some
simplification)

$I_{1}=2e^{-i\frac{\Theta}{2}}\int_{0}^{\epsilon^{\prime}}e^{-\frac{1}{2}\rho
U^{2}t}dU.$ (69)

where $\epsilon^{\prime}=\epsilon e^{i\Theta/2}$ is a real quantity (recall
that $U$ is the real-valued magnitude defined in (66)). Thus (69) is a real
integral, translated and rotated from its skewed path through saddle $k_{s1}$
in Fig. 11. We now replace (69) with its semi-infinite extension

$I_{1}\sim 2e^{-i\frac{\Theta}{2}}\int_{0}^{\infty}e^{-\frac{1}{2}\rho
U^{2}t}dU,~{}~{}t\to\infty,$ (70)

noting that (70) asymptotically approaches (69) as $t\to\infty$, due the to
dominant contribution at the saddle. The improper integral in (70) evaluates
exactly to $\sqrt{2\pi/(\rho t)}/2$ (see [Gradshteyn & Ryzhik, 2007] equation
3.321.3), and thus (70) becomes

$I_{1}\sim e^{-i\frac{\Theta}{2}}\sqrt{\frac{2\pi}{\rho t}},~{}~{}t\to\infty.$
(71)

If we apply the same technique as done above to $I_{2}$ in (LABEL:eq:Cauchy)
involving saddle $k_{s2}$ in Fig. 11, whose real imaginary part is the same
but real part opposite sign of $k_{s1}$ and steepest descent path is oriented
perpendicular to that of $k_{s1}$, we obtain

$I_{2}\sim e^{i\frac{\Theta}{2}}\sqrt{\frac{2\pi}{\rho t}},~{}~{}t\to\infty.$
(72)

Substituting (71) and (72) into (LABEL:eq:Cauchy), using (62) to recognize
that $\phi_{1}(k_{s1,2})=\left[b\left(\frac{x}{t}\right)-d\right]\pm
i\left[c-a\left(\frac{x}{t}\right)\right]$ and
$C_{1}(k_{s1})=\overline{C_{1}(k_{s2})}$, we obtain

$\displaystyle h\left(x,t\right)|_{\frac{x}{t}\neq
3}\sim\sqrt{\frac{2}{\rho\pi
t}}~{}e^{\left(b\frac{x}{t}-d\right)t}\textrm{~{}Real}\left\\{C_{1}(k_{s1})e^{i\left[\left(c-a\frac{x}{t}\right)t-\Theta/2\right]}\right\\},~{}~{}t\to\infty,~{}~{}a\neq
0,~{}~{}R\\!e\leq\cot(\theta)$ (73a)

where the $a\neq 0$ condition reminds us that we have assumed two contributing
saddles symmetric about the imaginary axis, as shown in Fig. 11 for saddles in
the lower-half plane for $x/t=2$. The result (73a) may be applied to
contributing saddles in the lower or upper half plane (which is the case for
$x/t=4$, see Fig. 12) if the definitions in (62) are modified as follows

$\displaystyle a=|\textrm{Real }(k_{s})|$ (73b) $\displaystyle b=-\textrm{
Imag }(k_{s})$ (73c) $\displaystyle c=|\textrm{Real }[\omega_{1}(k_{s})]|$
(73d) $\displaystyle d=-\textrm{ sgn }[\textrm{Imag }(ks)]~{}|\textrm{Imag
}[\omega_{1}(k_{s})]|$ (73e) $\displaystyle r=\left|\textrm{Imag
}\left(\frac{d^{2}\omega_{1}}{dk^{2}}\right)_{k_{s}}\right|$ (73f)
$\displaystyle q=-\textrm{ sgn }[\textrm{Imag }(ks)]\left|\textrm{Real
}\left(\frac{d^{2}\omega_{1}}{dk^{2}}\right)_{k_{s}}\right|.$ (73g)

Figure 12: Contours of constant $\textrm{Real}[{{\phi}_{1}}(k)]$ vs complex
$k={{k}_{r}}+i{{k}_{i}}$ at neutral stability ($R\\!e=\cot\theta$) for
$x/t=4$, $\theta=\pi/4$, and $W\\!e$=0.1. Notation is the same as in Fig. 10
and discussion is the same as in Fig. 11. For reference, the location of the
contributing saddles and related quantities necessary to determine the
asymptotic behavior using (62), (63), and (73a) are as follows: $a$=0.4528,
$b$=0.1584, $c$=1.4871, $d$=0.5053, $r$=0.7466, and $q$=5.1297.

For $R\\!e$$<\cot(\theta)$, as either $R\\!e$ is decreased or $\frac{x}{t}\to
3$ (from above or below) the saddles move closer together until they congeal
on the imaginary axis and thus only one saddle contributes, as shown in Fig.
13. Carrying out the above analysis for this case leads to

$\displaystyle h\left(x,t\right)|_{\frac{x}{t}\neq
3}\sim\frac{1}{\sqrt{2\rho\pi
t}}~{}e^{\left(b\frac{x}{t}-d\right)t}\textrm{~{}Real}\left[C_{1}(k_{s1})e^{i\left(ct-\Theta/2\right)}\right],~{}~{}t\to\infty,~{}~{}a=0,~{}~{}R\\!e<\cot(\theta).$
(73h)

where $b$, $c$, and $d$ correspond to $k_{s,1}$ in Fig. 13.

Figure 13: Contours of constant $\textrm{Real}[{{\phi}_{1}}(k)]$ vs complex
$k={{k}_{r}}+i{{k}_{i}}$ at neutral stability ($R\\!e=0.05\cot\theta$) for
$x/t=2$, $\theta=\pi/4$, and $W\\!e$=0.1. Notation is the same as in Fig. 9.
The orientation of the saddles enables a closed contour (indicated by dashed
lines) to be drawn through $k_{s,1}$ that includes a portion of the real axis,
C1, and also has $\textrm{Real}[{{\phi}_{1}}(k_{s,1})]$ as a maximum along the
contour – thus, the method of steepest decent may be used to determine the
long time asymptotic behavior of the integral involving $\phi_{1}(k)$ in (27c)
along the contour C1. Alternatively, the orientation of saddles $k_{s,2}$ and
$k_{s,3}$ makes them spurious and are not accessible by connection to any real
contour. The path of integration in (27c) is the infinite $k_{r}$-axis, i.e.
consists of the contour C=C${}_{\mathrm{L}}\cup$C${}_{1}\cup$CR and evaluated
in the same manner as discussed in Fig. 11. For reference, the location of the
contributing saddle, and related quantities necessary to determine the
asymptotic behavior using (62), (63), and (73h) are as follows: $a$=0,
$b$=1.3728, $c$=0, $d$=2.8420, $r$=1.8563, and $q$=0.