On the stability of waves in classically neutral flows

TL;DR

This paper investigates the behavior of wave stability in neutral flows, revealing algebraic growth at the stability boundary through analysis of exponential modes in two PDE systems, challenging traditional stability assumptions.

Contribution

It demonstrates algebraic growth at the neutral stability boundary for single-mode conditions in two PDE models, expanding understanding of wave behavior in neutral flows.

Findings

Algebraic growth occurs at the neutral stability boundary.

Propagation features of responses are characterized.

Analysis applies to both idealized and realistic flow models.

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, \textit{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.

Equations

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

$$D(k,\omega )=0,\text{ also written }\omega =\omega (k),$$

$$h\sim C{e}^{{\omega }_{i,\text{max}}t},\text{ as }t\to \infty ,$$

$${\omega }_{i,\text{max}}<0\text{: the flow is linearly stable}$$

$${\omega }_{i,\text{max}}>0\text{: the flow is linearly unstable}$$

$${\omega }_{i,\text{max}}=0\text{: the flow is neutrally stable}.$$

$$\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)$$

$$h(x,0)=h_0\delta (x),\frac{\partial h}{\partial t}(x,0)=u_0\delta (x),-\infty <x<\infty ,t>0,$$

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

$$\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),$$

$$h(x,0)=h_0\delta (x),\frac{\partial h}{\partial t}(x,0)=u_0\delta (x),-\infty <x<\infty ,t>0,$$

$$\hat{h}\left(k,t\right)=\underset{-\infty }{\overset{\infty }{\int }}h(x,t)e^{-ikx}dx,$$

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

$$\frac{d^2\hat{h}}{d{t}^2}+\left(k^2+2ick\right)\frac{d\hat{h}}{dt}+\left(k^4+ic{k}^3-c^2{k}^2\right)\hat{h}=f_0\delta (t),$$

$$\hat{h}(0)=h_0,\frac{d\hat{h}}{dt}(0)=u_0,t>0,$$

$$\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],$$

$$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}$$

$$h(x,t)$$

$$h(x,t)=$$

$$+\frac{h_{0}c}{\sqrt{3}}\text{Im}\left[\text{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}},$$

$$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].$$

$$h(x,t){|}_{\frac{x}{t}\ne c}\sim \frac{e^{-{\left(\frac{x}{t}-c\right)}^{2}t/8}}{\sqrt{3\pi t}}$$

$$\frac{\partial \bar{Q}}{\partial x}+\frac{\partial \bar{H}}{\partial t}=0$$

$$\frac{\partial \bar{Q}}{\partial t}+\frac{6}{5}\frac{\partial }{\partial x}\left(\frac{{\bar{Q}}^2}{\bar{H}}\right)=\frac{3}{Re}\left(\bar{H}-\frac{\bar{Q}}{{\bar{H}}^2}-\bar{H}\frac{\partial \bar{H}}{\partial x}\cot \theta \right)+\frac{1}{We}\bar{H}\frac{{\partial }^3\bar{H}}{\partial x^3}$$

$$\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},Re=\frac{\rho Q_0}{\mu },We=\frac{\rho {Q_0}^2}{\sigma H_0}.$$

$$H_0={\left(\frac{3Q_0\mu }{\rho g\sin \theta }\right)}^{1/3},$$

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

$$\frac{{\partial }^{2}h}{\partial t^2}+\frac{3}{Re}\frac{\partial h}{\partial t}+\frac{12}{5}\frac{{\partial }^{2}h}{\partial x\partial t}+\frac{1}{We}\frac{{\partial }^{4}h}{\partial x^4}+\left(\frac{6}{5}-\frac{3\cot \theta }{Re}\right)\frac{{\partial }^{2}h}{\partial x^2}+\frac{9}{Re}\frac{\partial h}{\partial x}=0.$$

$$h=h_0\delta (x),\frac{\partial h}{\partial t}=u_0\delta (x)\text{ at }t=0$$

$$h\to 0\text{ as }x\to \pm \infty .$$

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

Full text

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

[TABLE]

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

[TABLE]

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

[TABLE]

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]:

[TABLE]

[TABLE]

[TABLE]

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).

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,

[TABLE]

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

[TABLE]

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).

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:

[TABLE]

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

[TABLE]

[TABLE]

The Fourier transform of (6) yields

[TABLE]

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

[TABLE]

The inverse Fourier transform of (10) is

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

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).

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.

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]:

[TABLE]

[TABLE]

where:

[TABLE]

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

[TABLE]

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

[TABLE]

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:

[TABLE]

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:

[TABLE]

[TABLE]

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:

[TABLE]

where equating like powers leads to

[TABLE]

[TABLE]

[TABLE]

[TABLE]

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

[TABLE]

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.

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

and

[TABLE]

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

[TABLE]

where

[TABLE]

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

[TABLE]

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

[TABLE]

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)$

[TABLE]

[TABLE]

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)

[TABLE]

[TABLE]

[TABLE]

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.

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.

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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$):

[TABLE]

[TABLE]

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

[TABLE]

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

[TABLE]

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):

[TABLE]

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

[TABLE]

A.2

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

where

[TABLE]

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.

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:

[TABLE]

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

[TABLE]

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.

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

[TABLE]

[TABLE]

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

[TABLE]

[TABLE]

We thus see that

[TABLE]

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

[TABLE]

and thus

[TABLE]

[TABLE]

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

[TABLE]

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

[TABLE]

[TABLE]

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

[TABLE]

[TABLE]

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).

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

from which it follows from (62) that

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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].

[TABLE]

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)

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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