1 Introduction
The onset of instability of a horizontal interface separating two miscible solutions in a porous medium has been the subject of investigation in a number of theoretical and experimental studies. The problem finds application in several fields of natural science due to its direct relevance in innumerable industrial applications, such as surface and colloid physical chemistry, soil physics, petroleum engineering and powder metallurgy, in addition to its academic interest from a fluid mechanics point of view regarding the evolution of instabilities and the subsequent nonlinear dynamics (Philip Reference Philip1970; Huppert & Sparks Reference Huppert and Sparks1984; Homsy Reference Homsy1987; Schmitt Reference Schmitt1994; Schmitt et al. Reference Schmitt, Ledwell, Montgomery, Polzin and Toole2005; Trevelyan, Almarcha & De Wit Reference Trevelyan, Almarcha and De Wit2011; Huppert & Neufeld Reference Huppert and Neufeld2014).
The conditions for the interface between two immiscible fluids with different densities and viscosities in a porous medium to destabilise and deform into finger-like structures when subject to a gravitational field or an imposed pressure gradient were laid out by Saffman & Taylor (Reference Saffman and Taylor1958). This followed from the well-known Rayleigh–Taylor (RT) instability that occurs when a denser solution overlies a lighter one, which deforms the interface between them (Manickam & Homsy Reference Manickam and Homsy1995; Fernandez et al. Reference Fernandez, Kurowski, Limat and Petitjeans2001, Reference Fernandez, Kurowski, Petitjeans and Meiburg2002; Martin, Rakotomala & Salin Reference Martin, Rakotomala and Salin2002; Gandhi & Trevelyan Reference Gandhi and Trevelyan2014; Gopalakrishnan et al. Reference Gopalakrishnan, Carballido-Landeira, De Wit and Knaepen2017; De Paoli, Zonta & Soldati Reference De Paoli, Zonta and Soldati2019b).
In an initially statically stable system with a less dense solution overlying a denser solution, the miscible interface separating the two solutions can destabilise due to differential diffusion of solutes, where the driving energy comes from the component having the higher or lower diffusivity. When the solute in the denser lower solution diffuses faster than the solute in the lighter upper solution, a double-diffusive (DD) instability develops, which triggers a convective motion deforming the interface into fingers that develop symmetrically across it (Turner Reference Turner1979; Huppert & Turner Reference Huppert and Turner1981; Green Reference Green1984; Huppert & Sparks Reference Huppert and Sparks1984; Schmitt Reference Schmitt1994; Cooper, Glass & Tyler Reference Cooper, Glass and Tyler1997; Pringle & Glass Reference Pringle and Glass2002; Schmitt et al. Reference Schmitt, Ledwell, Montgomery, Polzin and Toole2005; Trevelyan et al. Reference Trevelyan, Almarcha and De Wit2011; Radko Reference Radko2013). Apart from the differential diffusion of solutes such as salt and sugar (for instance), which result in double diffusion (Pringle & Glass Reference Pringle and Glass2002), such phenomena are responsible for the thermohaline convection observed in oceans, which is triggered due to the differential diffusion of salt and heat (Schmitt Reference Schmitt1994; Schmitt et al. Reference Schmitt, Ledwell, Montgomery, Polzin and Toole2005).
Instability also arises when the solute in the lighter upper solution diffuses faster than the one in the lower denser solution. Regions with locally unstable density gradients develop over time induced by a differential diffusion effect, and is termed as diffusive-layer convection (DLC) (Turner & Stommel Reference Turner and Stommel1964; Griffiths Reference Griffiths1981; Turner Reference Turner1985; Trevelyan et al. Reference Trevelyan, Almarcha and De Wit2011). The faster diffusion of the upper solute towards the denser lower solution results in a depletion zone above the interface, and an accumulation below, triggering a convective motion on either side of the initial interface.
In all the scenarios presented so far, the differential diffusion of the solutes can lead to non-monotonic density profiles, which significantly influence the flow dynamics (Trevelyan et al. Reference Trevelyan, Almarcha and De Wit2011; Carballido-Landeira et al. Reference Carballido-Landeira, Trevelyan, Almarcha and De Wit2013; Gopalakrishnan et al. Reference Gopalakrishnan, Carballido-Landeira, Knaepen and De Wit2018). The RT instability involving a single species seldom occurs in the various porous media flows observed in nature and in industrial processes, as double-diffusive effects usually coexist, due to the difference in the differential diffusion of solutes. This coupling between diffusive and convective processes can result in a double-diffusive flux that could be much larger than in a single solutal system (Turner Reference Turner1985).
The stability, or the instability, of these buoyancy-driven flows has been mainly addressed using linear base states for the time-evolving density profiles, or by means of numerical linear stability analyses. In the context of DD instabilities, the limiting conditions for salt fingering to occur at an interface were provided by Huppert & Sparks (Reference Huppert and Sparks1984) by considering linear base states. They showed that an infinitesimal perturbation will develop into fingers when 
$R<\unicode[STIX]{x1D70F}^{3/2}$, where 
$R$ is the buoyancy ratio and 
$\unicode[STIX]{x1D70F}$ is the diffusion coefficient ratio (that of heat to that of salt, in this instance). The buoyancy ratio 
$R$ expresses the relative solutal buoyancy contribution of each species to the density profile. Trevelyan et al. (Reference Trevelyan, Almarcha and De Wit2011) carried out a stability analysis of various buoyancy-driven flows in a two-dimensional porous medium and Hele-Shaw cells. Using a quasi-steady-state approximation, they carried out a linear stability analysis of time-evolving base-state profiles numerically and mapped out the various regimes of instability in the 
$(R,\unicode[STIX]{x1D70F})$ parameter space. The neutral stability curve obtained numerically was in agreement with that of Huppert & Sparks (Reference Huppert and Sparks1984). More recently, using analytical methods, the onset conditions for an RT instability to develop were investigated by Gandhi & Trevelyan (Reference Gandhi and Trevelyan2014), where they used a linear combination of step functions for the base-state concentration profiles.
In the present work we provide analytically the conditions for instability to develop in the various aforementioned buoyancy-driven flows in porous media. Specifically, we address ourselves to this problem and ask: Can we predict whether the interface separating two miscible solutions inside a vertically oriented two-dimensional porous medium or a Hele-Shaw cell is stable, or unstable, from the base-state density profile of the concentrations? Or, alternatively, can we provide the conditions for instability in these flows?
The analytical approach adopted in the present work has been previously used in the context of Rayleigh–Bénard convection, instability of shear layers and Taylor–Couette flows, among other hydrodynamic instabilities (Roberts Reference Roberts1960; Chandrasekhar Reference Chandrasekhar1961). Unlike in the earlier analytical works, where linear base states were considered, we have used time-evolving base-state profiles to extract information regarding the stability of the system. For instance, in the parameter space corresponding to 
$\unicode[STIX]{x1D70F}>1$ and 
$R>1$, Trevelyan et al. (Reference Trevelyan, Almarcha and De Wit2011) numerically found that the maximum growth rates were always real, which is confirmed in this study. Specifically, when the base-state density profiles are strictly monotonically increasing downwards, such that the density gradient of the base state is positive, we have shown that instability develops if and only if the second derivative of the base-state density profile is zero at a point on either side of the initial interface such that it changes sign. In these regimes, instability arises as stationary modes, with real eigenvalues. Whereas, when there is an interval in the spatial domain where the first derivative of the base-state density profile is negative, the instability may develop as a stationary or as an oscillatory mode (complex eigenvalue). The neutral stability curve obtained analytically is in agreement with the earlier findings. The fact that the eigenvalues are strictly real in the context of single-species RT instabilities also follows from the work presented here.
To this end, the geometry and the governing equations are introduced in § 2, along with the base-state density profiles in the non-dimensional 
$(R,\unicode[STIX]{x1D70F})$ plane. In § 3 we reduce this to a characteristic value problem from which the conditions for instability can be obtained, which is discussed in § 4. After considering some special cases, namely single-species RT instability and when 
$\unicode[STIX]{x1D70F}=0$, the paper finishes in § 5 with the main conclusions.
2 Geometry and governing equations
We consider a two-dimensional porous medium in which gravity points downwards and two miscible solutions are in contact along a horizontal interface such as that sketched in figure 1. The upper solution contains a solute A while the lower solution contains a solute B, with initial concentrations 
$A_{0}$ and 
$B_{0}$, respectively. The 
$y$-axis is defined as the horizontal axis where the solutions are initially in contact, and the 
$x$-axis as the vertical one increasing in the downwards direction, such that 
$x<0$ is the upper region and 
$x>0$ is the lower region. By considering dilute solutions, the diffusion coefficients 
$D_{A}$ and 
$D_{B}$ of species A and B, respectively, can be assumed constant. Moreover, the density varies linearly with the concentrations as 
$\unicode[STIX]{x1D70C}(A,B)=\unicode[STIX]{x1D70C}_{0}[1+\unicode[STIX]{x1D6FC}_{A}A+\unicode[STIX]{x1D6FC}_{B}B]$, where 
$\unicode[STIX]{x1D70C}_{0}$ is the density of the pure solvent, 
$A$ and 
$B$ are the concentrations of the respective species, and 
$\unicode[STIX]{x1D6FC}_{A}$ and 
$\unicode[STIX]{x1D6FC}_{B}$ are their corresponding solutal expansion coefficients defined as 
$\unicode[STIX]{x1D6FC}_{A}=(1/\unicode[STIX]{x1D70C}_{0})\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}/\unicode[STIX]{x2202}A$ and 
$\unicode[STIX]{x1D6FC}_{B}=(1/\unicode[STIX]{x1D70C}_{0})\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}/\unicode[STIX]{x2202}B$. The dynamics of the buoyancy-driven instabilities of this miscible interface are described by Darcy’s equation coupled to advection–diffusion equations for the concentrations 
$A$ and 
$B$. The set of equations governing the system can be non-dimensionalised by considering the characteristic velocity 
$\mathbb{U}=gK\unicode[STIX]{x1D6FC}_{A}A_{0}/\unicode[STIX]{x1D707}$, length 
$\mathbb{L}=D_{A}\unicode[STIX]{x1D719}/\mathbb{U}$ and time 
$\mathbb{T}=\mathbb{L}/\mathbb{U}$ scales, where 
$g$ is the magnitude of the acceleration due to gravity, 
$K$ is the permeability, 
$\unicode[STIX]{x1D707}$ is the dynamic viscosity and 
$\unicode[STIX]{x1D719}$ is the porosity (taken as 
$1$). The resulting non-dimensional equations read (Trevelyan et al. Reference Trevelyan, Almarcha and De Wit2011; Gopalakrishnan et al. Reference Gopalakrishnan, Carballido-Landeira, Knaepen and De Wit2018):
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}p=-\boldsymbol{u}+(A+RB)\hat{\boldsymbol{i}}_{x}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & A_{t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}A=\unicode[STIX]{x1D6FB}^{2}A, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & B_{t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B=\unicode[STIX]{x1D70F}\unicode[STIX]{x1D6FB}^{2}B, & \displaystyle\end{eqnarray}$$ where 
$p$ is the pressure, 
$\boldsymbol{u}$ is the flow velocity, 
$\hat{\boldsymbol{i}}_{x}$ is the unit vector in the direction of gravity, 
$R=\unicode[STIX]{x1D6FC}_{B}B_{0}/\unicode[STIX]{x1D6FC}_{A}A_{0}$ is the dimensionless buoyancy ratio and 
$\unicode[STIX]{x1D70F}=D_{B}/D_{A}$ is the ratio of the diffusion coefficients (Trevelyan et al. Reference Trevelyan, Almarcha and De Wit2011; Gopalakrishnan et al. Reference Gopalakrishnan, Carballido-Landeira, Knaepen and De Wit2018). The buoyancy ratio 
$R$ expresses the relative contributions of species B and A to the dimensional density profile and is also equal to 
$R=(\unicode[STIX]{x1D70C}_{B}-\unicode[STIX]{x1D70C}_{0})/(\unicode[STIX]{x1D70C}_{A}-\unicode[STIX]{x1D70C}_{0})$, where 
$\unicode[STIX]{x1D70C}_{A}=\unicode[STIX]{x1D70C}_{0}(1+\unicode[STIX]{x1D6FC}_{A}A_{0})$ and 
$\unicode[STIX]{x1D70C}_{B}=\unicode[STIX]{x1D70C}_{0}(1+\unicode[STIX]{x1D6FC}_{B}B_{0})$ are the densities of the upper and lower solutions, respectively. The concentrations are non-dimensionalised using 
$A_{0}$ while the density is scaled as 
$(\unicode[STIX]{x1D70C}(A,B)/\unicode[STIX]{x1D70C}_{0}-1)/\unicode[STIX]{x1D6FC}_{A}A_{0}$. As initial conditions, we impose in the top layer (
$x<0$) 
$A=1$, 
$B=0$, while in the bottom layer (
$x>0$) 
$A=0$, 
$B=1$, with no velocity in the entire domain (
$\boldsymbol{u}=\mathbf{0}$). The boundary conditions are 
$A=1$, 
$B=0$ for 
$x\rightarrow -\infty$, and 
$A=0$, 
$B=1$ for 
$x\rightarrow \infty$.
Figure 1. Sketch of the initial physical problem.
In the absence of convection, the base-state concentration profiles, assuming a domain of infinite height, are given by
$$\begin{eqnarray}\bar{A}(x,t)=\frac{1}{2}\,\text{erfc}\left(\frac{x}{2\sqrt{t}}\right),\quad \bar{B}(x,t)=\frac{1}{2}\,\text{erfc}\left(-\frac{x}{2\sqrt{\unicode[STIX]{x1D70F}t}}\right),\end{eqnarray}$$while the non-dimensional base-state density profile can be constructed as
$$\begin{eqnarray}\bar{\unicode[STIX]{x1D70C}}(x,t)=\bar{A}(x,t)+R\bar{B}(x,t).\end{eqnarray}$$
Figure 2. Typical base-state density profiles 
$\bar{\unicode[STIX]{x1D70C}}(x,t)$ in the 
$(R,\unicode[STIX]{x1D70F})$ parameter space. Base-state density profiles are monotonically increasing in the shaded region.
The base state of our problem follows a diffusive dynamics satisfying (2.3) and (2.4) with 
$\boldsymbol{u}=\mathbf{0}$ along with the initial conditions specified earlier. This is in contrast with the earlier analytical studies on the stability of buoyancy-driven flows, where linear base states were considered (Huppert & Sparks Reference Huppert and Sparks1984; Gandhi & Trevelyan Reference Gandhi and Trevelyan2014). As noted in Trevelyan et al. (Reference Trevelyan, Almarcha and De Wit2011), the density is antisymmetric about 
$x=0$, as 
$\bar{\unicode[STIX]{x1D70C}}(x,t)+\bar{\unicode[STIX]{x1D70C}}(-x,t)=1+R$. This also implies that the convective patterns in porous media flows evolve the same way on both sides of the initial interface (Trevelyan et al. Reference Trevelyan, Almarcha and De Wit2011).
The different buoyancy-driven instabilities discussed in § 1, namely the Rayleigh–Taylor (RT), double-diffusive (DD) and diffusive-layer convection (DLC), can be obtained depending on the values of parameters 
$R$ and 
$\unicode[STIX]{x1D70F}$. The RT instabilities develop when 
$R<1$, as we have a ‘heavy on top of light’ configuration; whereas for 
$R>1$, the system is initially in a stable stratification (‘light on top of heavy’). The single-species RT instability occurs when 
$\unicode[STIX]{x1D70F}=1$ and 
$R<1$, as the diffusion coefficients are the same. When the solutes diffuse at different rates 
$\unicode[STIX]{x1D70F}\neq 1$, DD and DLC mechanisms can destabilise the miscible interface. The DD instabilities are at play when 
$\unicode[STIX]{x1D70F}>1$, whereas DLC mechanisms develop over time when 
$\unicode[STIX]{x1D70F}<1$. Indeed, when 
$R<1$, DD and DLC instabilities are coupled with RT effects. Typical base-state density profiles corresponding to RT, DD and DLC are shown in figure 2 in the 
$(R,\unicode[STIX]{x1D70F})$ plane. A more complete description of the various features of the base-state density profiles pertaining to different instability mechanisms can be found in Trevelyan et al. (Reference Trevelyan, Almarcha and De Wit2011).
We now turn our attention to the problem at hand, i.e. to obtain the conditions for this system to be linearly unstable along with the various features of the instability, from the base-state density profiles.
3 Reduction to a characteristic value problem
As the flow is incompressible, we consider the streamfunction 
$\unicode[STIX]{x1D713}$ formulation, using 
$u=\unicode[STIX]{x1D713}_{y}$ and 
$v=-\unicode[STIX]{x1D713}_{x}$, where 
$\boldsymbol{u}=(u,v)$, thereby satisfying 
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0$. Taking the curl of (2.1) and substituting 
$\unicode[STIX]{x1D713}$ into (2.3) and (2.4), we have the following system of equations (Trevelyan et al. Reference Trevelyan, Almarcha and De Wit2011):
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}=A_{y}+RB_{y}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & A_{t}+\unicode[STIX]{x1D713}_{y}A_{x}-\unicode[STIX]{x1D713}_{x}A_{y}=\unicode[STIX]{x1D6FB}^{2}A, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & B_{t}+\unicode[STIX]{x1D713}_{y}B_{x}-\unicode[STIX]{x1D713}_{x}B_{y}=\unicode[STIX]{x1D70F}\unicode[STIX]{x1D6FB}^{2}A. & \displaystyle\end{eqnarray}$$Introducing perturbations in the form of normal modes to the base-state solutions (which are assumed to vary more slowly than the perturbations) we have
$$\begin{eqnarray}[\unicode[STIX]{x1D713},A,B]=[0,\bar{A},\bar{B}]+\unicode[STIX]{x1D716}\exp (\unicode[STIX]{x1D70E}t+\text{i}ky)[\text{i}k^{-1}\hat{\unicode[STIX]{x1D713}},\hat{A},\hat{B}],\end{eqnarray}$$ where 
$\unicode[STIX]{x1D716}$ is a small parameter, 
$k$ denotes the wavenumber and the complex frequency 
$\unicode[STIX]{x1D70E}~(=\unicode[STIX]{x1D70E}_{r}+\text{i}\unicode[STIX]{x1D70E}_{i})$ is the eigenvalue. The overbars indicate base-state variables. With a positive growth rate, 
$\unicode[STIX]{x1D70E}_{r}>0$, the perturbation grows exponentially in time and the system is linearly unstable. While adopting the normal mode approach to describe any arbitrary disturbance, we have made a quasi-steady-state approximation (QSSA) by assuming that the base state is varying more slowly than the perturbations (Farrell & Ioannou Reference Farrell and Ioannou1996). Although there is only one time scale in this problem, 
$\mathbb{T}$, which would imply that the base state and disturbances both evolve with this characteristic time scale, the base state under consideration here smears out diffusively in time, and hence the QSSA is justified as noted in the earlier works by Tan & Homsy (Reference Tan and Homsy1986) and Trevelyan et al. (Reference Trevelyan, Almarcha and De Wit2011). Here it is pertinent to point out that another approach to the QSSA has recently been adopted in the works of Kim & Choi (Reference Kim and Choi2011), Pramanik & Mishra (Reference Pramanik and Mishra2013) and Hota & Mishra (Reference Hota and Mishra2018) by using a similarity transformation 
$(\unicode[STIX]{x1D709},t)$, where 
$\unicode[STIX]{x1D709}=x/\sqrt{t}$, such that the base-state concentration profile does not explicitly depend on time, and hence can be translated along the flow direction as time passes. At a time instant near zero, such a self-similar QSSA method is shown to give more accurate results in comparison to the classical QSSA, though excellent agreement is observed for later times.
We now consider the base-state solution at a given time 
$t$ and treat this frozen profile as though it were steady. Linearising system (3.1)–(3.3) in 
$\unicode[STIX]{x1D716}$, we obtain
$$\begin{eqnarray}\displaystyle & \hat{\unicode[STIX]{x1D713}}_{xx}=k^{2}(\hat{\unicode[STIX]{x1D713}}+\hat{A}+R\hat{B}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D70E}\hat{A}=\hat{A}_{xx}-k^{2}\hat{A}+\bar{A}_{x}\hat{\unicode[STIX]{x1D713}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D70E}\hat{B}=\unicode[STIX]{x1D70F}(\hat{B}_{xx}-k^{2}\hat{B})+\bar{B}_{x}\hat{\unicode[STIX]{x1D713}}. & \displaystyle\end{eqnarray}$$ If we define 
${\mathcal{L}}=(D^{2}-k^{2})$, where 
$D=\text{d}/\text{d}x$, we have
$$\begin{eqnarray}\displaystyle & {\mathcal{L}}\hat{\unicode[STIX]{x1D713}}=k^{2}(\hat{A}+R\hat{B}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D70E}\hat{A}={\mathcal{L}}\hat{A}+\bar{A}_{x}\hat{\unicode[STIX]{x1D713}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D70E}\hat{B}=\unicode[STIX]{x1D70F}{\mathcal{L}}\hat{B}+\bar{B}_{x}\hat{\unicode[STIX]{x1D713}}. & \displaystyle\end{eqnarray}$$The boundary conditions with respect to which the above equations must be solved are given by
$$\begin{eqnarray}|\hat{A}|=|\hat{B}|=|\hat{\unicode[STIX]{x1D713}}|=0,\quad \text{for}\;x=\pm \infty ,\end{eqnarray}$$as we have considered an infinite domain with the perturbations decaying at the corresponding vertical limits.
We have 
$\bar{\unicode[STIX]{x1D70C}}=\bar{A}+R\bar{B}$ (2.6), and consequently 
$\bar{\unicode[STIX]{x1D70C}}_{x}=\bar{A}_{x}+R\bar{B}_{x}$. Multiplying (3.10) by 
$R$ gives
$$\begin{eqnarray}\unicode[STIX]{x1D70E}R\hat{B}=R\unicode[STIX]{x1D70F}{\mathcal{L}}\hat{B}+R\bar{B}_{x}\hat{\unicode[STIX]{x1D713}},\end{eqnarray}$$and summing with (3.9), we get
$$\begin{eqnarray}\unicode[STIX]{x1D70E}\hat{A}+\unicode[STIX]{x1D70E}R\hat{B}={\mathcal{L}}\hat{A}+R\unicode[STIX]{x1D70F}{\mathcal{L}}\hat{B}+\bar{\unicode[STIX]{x1D70C}}_{x}\hat{\unicode[STIX]{x1D713}}.\end{eqnarray}$$ With 
$R\hat{B}=(1/k^{2}){\mathcal{L}}\hat{\unicode[STIX]{x1D713}}-\hat{A}$ (from (3.8)), equation (3.13) reduces to
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D70E}\hat{A}+\unicode[STIX]{x1D70E}[(1/k^{2}){\mathcal{L}}\hat{\unicode[STIX]{x1D713}}-\hat{A}]={\mathcal{L}}\hat{A}+\unicode[STIX]{x1D70F}{\mathcal{L}}[(1/k^{2}){\mathcal{L}}\hat{\unicode[STIX]{x1D713}}-\hat{A}]+\bar{\unicode[STIX]{x1D70C}}_{x}\hat{\unicode[STIX]{x1D713}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x1D70E}}{k^{2}}{\mathcal{L}}\hat{\unicode[STIX]{x1D713}}={\mathcal{L}}\hat{A}+\frac{\unicode[STIX]{x1D70F}}{k^{2}}{\mathcal{L}}^{2}\hat{\unicode[STIX]{x1D713}}-\unicode[STIX]{x1D70F}{\mathcal{L}}\hat{A}+\bar{\unicode[STIX]{x1D70C}}_{x}\hat{\unicode[STIX]{x1D713}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left[\frac{\unicode[STIX]{x1D70E}}{k^{2}}{\mathcal{L}}-\frac{\unicode[STIX]{x1D70F}}{k^{2}}{\mathcal{L}}^{2}-\bar{\unicode[STIX]{x1D70C}}_{x}\right]\hat{\unicode[STIX]{x1D713}}=(1-\unicode[STIX]{x1D70F}){\mathcal{L}}\hat{A}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle [\unicode[STIX]{x1D70E}{\mathcal{L}}-\unicode[STIX]{x1D70F}{\mathcal{L}}^{2}-\bar{\unicode[STIX]{x1D70C}}_{x}k^{2}]\hat{\unicode[STIX]{x1D713}}=k^{2}(1-\unicode[STIX]{x1D70F}){\mathcal{L}}\hat{A}. & \displaystyle\end{eqnarray}$$
${\mathcal{L}}\hat{A}=\unicode[STIX]{x1D70E}\hat{A}-\bar{A}_{x}\hat{\unicode[STIX]{x1D713}}$ (from (3.9)), equation (3.14d) can be rewritten as 
$$\begin{eqnarray}\displaystyle & [\unicode[STIX]{x1D70E}{\mathcal{L}}-\unicode[STIX]{x1D70F}{\mathcal{L}}^{2}-\bar{\unicode[STIX]{x1D70C}}_{x}k^{2}]\hat{\unicode[STIX]{x1D713}}=k^{2}(1-\unicode[STIX]{x1D70F})(\unicode[STIX]{x1D70E}\hat{A}-\bar{A}_{x}\hat{\unicode[STIX]{x1D713}}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & [\unicode[STIX]{x1D70E}{\mathcal{L}}-\unicode[STIX]{x1D70F}{\mathcal{L}}^{2}-\bar{\unicode[STIX]{x1D70C}}_{x}k^{2}+k^{2}(1-\unicode[STIX]{x1D70F})\bar{A}_{x}]\hat{\unicode[STIX]{x1D713}}=k^{2}(1-\unicode[STIX]{x1D70F})\unicode[STIX]{x1D70E}\hat{A}. & \displaystyle\end{eqnarray}$$
$\hat{\unicode[STIX]{x1D713}}$ from (3.9) (
$\hat{\unicode[STIX]{x1D713}}=(\unicode[STIX]{x1D70E}\hat{A}-{\mathcal{L}}\hat{A})/\bar{A}_{x}$) in (3.15b), we get 
$$\begin{eqnarray}\unicode[STIX]{x1D70E}^{2}{\mathcal{L}}\hat{A}-\unicode[STIX]{x1D70E}[(1+\unicode[STIX]{x1D70F}){\mathcal{L}}^{2}+\bar{\unicode[STIX]{x1D70C}}_{x}k^{2}]\hat{A}+\unicode[STIX]{x1D70F}{\mathcal{L}}^{3}\hat{A}-k^{2}\unicode[STIX]{x1D701}(x){\mathcal{L}}\hat{A}=0,\end{eqnarray}$$ where 
$\unicode[STIX]{x1D701}(x)=\bar{A}_{x}(1-\unicode[STIX]{x1D70F})-\bar{\unicode[STIX]{x1D70C}}_{x}$.
The above equation along with the boundary conditions (3.11) constitute a characteristic value problem, which is indeed a formidable one if it has to be solved in any direct manner. Fortunately, it will appear that this is not necessary to obtain the conditions for instability from the characteristics of the base-state density profile 
$\bar{\unicode[STIX]{x1D70C}}(x,t)$. For this, a method based on a variational principle (Roberts Reference Roberts1960; Chandrasekhar Reference Chandrasekhar1961) can be used, which is described in the following section.
4 Conditions for instability
We multiply (3.16) by 
$\hat{A}^{\ast }$, the complex conjugate of 
$\hat{A}$, and integrate over the range of 
$x$, the vertical extent of the domain, which gives
$$\begin{eqnarray}\int _{-\infty }^{\infty }\unicode[STIX]{x1D70E}^{2}\hat{A}^{\ast }{\mathcal{L}}\hat{A}\,\text{d}x-\int _{-\infty }^{\infty }\hat{A}^{\ast }\unicode[STIX]{x1D70E}[(1+\unicode[STIX]{x1D70F}){\mathcal{L}}^{2}+\bar{\unicode[STIX]{x1D70C}}_{x}k^{2}]\hat{A}\,\text{d}x+\int _{-\infty }^{\infty }\hat{A}^{\ast }(\unicode[STIX]{x1D70F}{\mathcal{L}}^{3}\hat{A}-k^{2}\unicode[STIX]{x1D701}(x){\mathcal{L}}\hat{A})\,\text{d}x=0.\end{eqnarray}$$ It makes sense to write the above integrals, as the operator 
${\mathcal{L}}=(D^{2}-k^{2})$ is self-adjoint, implying that it is closed and bounded on 
$\mathbb{R}$, where 
$\mathbb{R}$ is the extended real line.
To simplify the terms corresponding to each integral, it is convenient to define the following:
$$\begin{eqnarray}\displaystyle & \displaystyle \mathbb{I}_{1}=\int _{-\infty }^{\infty }\unicode[STIX]{x1D70E}^{2}\hat{A}^{\ast }{\mathcal{L}}\hat{A}\,\text{d}x, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \mathbb{I}_{2}=\int _{-\infty }^{\infty }-\unicode[STIX]{x1D70E}\hat{A}^{\ast }[(\unicode[STIX]{x1D70F}+1){\mathcal{L}}^{2}+\bar{\unicode[STIX]{x1D70C}}_{x}k^{2}]\hat{A}\,\text{d}x, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \mathbb{I}_{3}=\int _{-\infty }^{\infty }\hat{A}^{\ast }(\unicode[STIX]{x1D70F}{\mathcal{L}}^{3}\hat{A}-k^{2}\unicode[STIX]{x1D701}(x){\mathcal{L}}\hat{A})\,\text{d}x. & \displaystyle\end{eqnarray}$$ We have 
$\mathbb{I}_{1}$, which can be rewritten as
$$\begin{eqnarray}\displaystyle \mathbb{I}_{1} & = & \displaystyle \int _{-\infty }^{\infty }\unicode[STIX]{x1D70E}^{2}\hat{A}^{\ast }(D^{2}-k^{2})\hat{A}\,\text{d}x\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \int _{-\infty }^{\infty }\unicode[STIX]{x1D70E}^{2}\hat{A}^{\ast }D^{2}\hat{A}\,\text{d}x-\int _{-\infty }^{\infty }\unicode[STIX]{x1D70E}^{2}\hat{A}^{\ast }k^{2}\hat{A}\,\text{d}x\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \unicode[STIX]{x1D70E}^{2}[\hat{A}^{\ast }D\hat{A}]_{-\infty }^{\infty }-\int _{-\infty }^{\infty }D\hat{A}^{\ast }D\hat{A}\,\text{d}x-\int _{-\infty }^{\infty }\unicode[STIX]{x1D70E}^{2}|k\hat{A}|^{2}\,\text{d}x\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \int _{-\infty }^{\infty }-\unicode[STIX]{x1D70E}^{2}|D\hat{A}|^{2}\,\text{d}x-\int _{-\infty }^{\infty }\unicode[STIX]{x1D70E}^{2}|k\hat{A}|^{2}\,\text{d}x\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \int _{-\infty }^{\infty }-\unicode[STIX]{x1D70E}^{2}[|D\hat{A}|^{2}+|k\hat{A}|^{2}]\,\text{d}x,\end{eqnarray}$$
$\mathbb{I}_{2}$ and 
$\mathbb{I}_{3}$ can be simplified as (see appendix A) 
$$\begin{eqnarray}\mathbb{I}_{2}=\int _{-\infty }^{\infty }-\unicode[STIX]{x1D70E}[(\unicode[STIX]{x1D70F}+1)[(D^{2}-k^{2})\hat{A}]^{2}+\bar{\unicode[STIX]{x1D70C}}_{x}|k\hat{A}|^{2}]\,\text{d}x,\end{eqnarray}$$
$$\begin{eqnarray}\mathbb{I}_{3}=\int _{-\infty }^{\infty }-\unicode[STIX]{x1D70F}[|D(D^{2}-k^{2})\hat{A}|^{2}+|k(D^{2}-k^{2})\hat{A}|^{2}]\,\text{d}x+\int _{-\infty }^{\infty }k^{2}\unicode[STIX]{x1D701}(x)[|D\hat{A}|^{2}+|k\hat{A}|^{2}]\,\text{d}x.\end{eqnarray}$$ Thus (4.1), which was broken down as 
$\mathbb{I}_{1}$, 
$\mathbb{I}_{2}$ and 
$\mathbb{I}_{3}$, reads as follows:
$$\begin{eqnarray}\displaystyle & & \displaystyle \int _{-\infty }^{\infty }-\unicode[STIX]{x1D70E}^{2}[|D\hat{A}|^{2}+|k\hat{A}|^{2}]\,\text{d}x+\int _{-\infty }^{\infty }-\unicode[STIX]{x1D70E}[(\unicode[STIX]{x1D70F}+1)[(D^{2}-k^{2})\hat{A}]^{2}+\bar{\unicode[STIX]{x1D70C}}_{x}|k\hat{A}|^{2}]\,\text{d}x+\cdots \nonumber\\ \displaystyle & & \displaystyle \quad +\int _{-\infty }^{\infty }\!-\unicode[STIX]{x1D70F}[|D(D^{2}\!-\!k^{2})\hat{A}|^{2}+|k(D^{2}-k^{2})\hat{A}|^{2}]\,\text{d}x+\!\int _{-\infty }^{\infty }k^{2}\unicode[STIX]{x1D701}(x)[|D\hat{A}|^{2}+|k\hat{A}|^{2}]\,\text{d}x=0.\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$Thus we have
$$\begin{eqnarray}\int _{-\infty }^{\infty }(\unicode[STIX]{x1D70E}^{2}{\mathcal{X}}+\unicode[STIX]{x1D70E}{\mathcal{Y}}+{\mathcal{Z}})\,\text{d}x=0,\end{eqnarray}$$where
$$\begin{eqnarray}\displaystyle & {\mathcal{X}}=[|D\hat{A}|^{2}+|k\hat{A}|^{2}], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\mathcal{Y}}=[(1+\unicode[STIX]{x1D70F})[(D^{2}-k^{2})\hat{A}]^{2}+\bar{\unicode[STIX]{x1D70C}}_{x}|k\hat{A}|^{2}], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\mathcal{Z}}=\unicode[STIX]{x1D70F}[|D(D^{2}-k^{2})\hat{A}|^{2}+|k(D^{2}-k^{2})\hat{A}|^{2}]-k^{2}\unicode[STIX]{x1D701}(x)[|D\hat{A}|^{2}+|k\hat{A}|^{2}]. & \displaystyle\end{eqnarray}$$ Note that the terms 
${\mathcal{X}},{\mathcal{Y}},{\mathcal{Z}}$ are real, whereas 
$\unicode[STIX]{x1D70E}$ may have real (
$\unicode[STIX]{x1D70E}_{r}$) and imaginary (
$\unicode[STIX]{x1D70E}_{i}$) parts. If 
$\unicode[STIX]{x1D70E}_{r}$ is positive, then the system would have perturbations that would grow exponentially in time, and is linearly unstable; whereas a negative 
$\unicode[STIX]{x1D70E}_{r}$ implies that the system is linearly stable. The perturbations could be oscillatory or stationary, with a non-zero value of 
$\unicode[STIX]{x1D70E}_{i}$ corresponding to the former, and 
$\unicode[STIX]{x1D70E}_{i}=0$ to the latter. We now look for the conditions for instability. Equation (4.9) reads as
$$\begin{eqnarray}\int _{-\infty }^{\infty }[(\unicode[STIX]{x1D70E}_{r}^{2}-\unicode[STIX]{x1D70E}_{i}^{2}+\text{i}2\unicode[STIX]{x1D70E}_{r}\unicode[STIX]{x1D70E}_{i}){\mathcal{X}}+(\unicode[STIX]{x1D70E}_{r}+\text{i}\unicode[STIX]{x1D70E}_{i}){\mathcal{Y}}+{\mathcal{Z}}]\,\text{d}x=0.\end{eqnarray}$$Considering the imaginary part, we have
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \int _{-\infty }^{\infty }\unicode[STIX]{x1D70E}_{i}(2\unicode[STIX]{x1D70E}_{r}{\mathcal{X}}+{\mathcal{Y}})\,\text{d}x=0,\\ \displaystyle \int _{-\infty }^{\infty }\unicode[STIX]{x1D70E}_{i}(2\unicode[STIX]{x1D70E}_{r}[|D\hat{A}|^{2}+|k\hat{A}|^{2}]+[(1+\unicode[STIX]{x1D70F})[(D^{2}-k^{2})\hat{A}]^{2}+\bar{\unicode[STIX]{x1D70C}}_{x}|k\hat{A}|^{2}])\,\text{d}x=0.\end{array}\right\}\end{eqnarray}$$ If we consider 
$\unicode[STIX]{x1D70E}_{i}\neq 0$ (oscillatory modes), the terms in parentheses are strictly positive, with the exception of the gradient of the base-state density profile 
$\bar{\unicode[STIX]{x1D70C}}_{x}$, which could be positive or negative or both within the domain, and 
$\unicode[STIX]{x1D70E}_{r}$. For instability, we know that 
$\unicode[STIX]{x1D70E}_{r}>0$. This implies that there must be a region within the domain where 
$\bar{\unicode[STIX]{x1D70C}}_{x}<0$, which is a sufficient condition for instability.
So we have our first necessary criterion for instability with oscillatory modes: there must be a region within the domain where 
$\bar{\unicode[STIX]{x1D70C}}_{x}<0$. On the contrary, if 
$\bar{\unicode[STIX]{x1D70C}}_{x}$ is positive throughout the domain, then the system is uniquely unstable to stationary modes such that 
$\unicode[STIX]{x1D70E}_{i}=0$.
We now look for the conditions for the flow to be unstable when 
$\bar{\unicode[STIX]{x1D70C}}_{x}$ is positive throughout the domain. Using the expression for 
$\bar{\unicode[STIX]{x1D70C}}$ (2.6) and differentiating with respect to 
$x$, we have the expression for 
$\bar{\unicode[STIX]{x1D70C}}_{x}$ given by
$$\begin{eqnarray}\bar{\unicode[STIX]{x1D70C}}_{x}=\frac{-\exp [-x^{2}/4t]}{2\sqrt{\unicode[STIX]{x03C0}t}}+\frac{\exp [-x^{2}/4t\unicode[STIX]{x1D70F}]R}{2\sqrt{\unicode[STIX]{x03C0}t\unicode[STIX]{x1D70F}}}.\end{eqnarray}$$The density profiles are monotonically increasing downwards when
$$\begin{eqnarray}1\leqslant \unicode[STIX]{x1D70F}\leqslant R^{2},\end{eqnarray}$$ and 
$\bar{\unicode[STIX]{x1D70C}}_{x}$ is positive throughout the domain (the shaded region in figure 2). For these parameter settings, the system is unstable to purely real eigenvalues (
$\unicode[STIX]{x1D70E}_{i}=0$). Using 
$\unicode[STIX]{x1D70E}_{i}=0$ in (4.13), we have
$$\begin{eqnarray}\int _{-\infty }^{\infty }(\unicode[STIX]{x1D70E}_{r}^{2}{\mathcal{X}}+\unicode[STIX]{x1D70E}_{r}{\mathcal{Y}}+{\mathcal{Z}})\,\text{d}x=0,\end{eqnarray}$$ and substituting the expressions for 
${\mathcal{X}}$, 
${\mathcal{Y}}$ and 
${\mathcal{Z}}$ gives
$$\begin{eqnarray}\displaystyle & & \displaystyle \int _{-\infty }^{\infty }(\unicode[STIX]{x1D70E}_{r}^{2}[|D\hat{A}|^{2}+|k\hat{A}|^{2}]+\unicode[STIX]{x1D70E}_{r}[(1+\unicode[STIX]{x1D70F})[(D^{2}-k^{2})\hat{A}]^{2}+\bar{\unicode[STIX]{x1D70C}}_{x}|k\hat{A}|^{2}]\nonumber\\ \displaystyle & & \displaystyle \quad +\,\unicode[STIX]{x1D70F}[|D(D^{2}-k^{2})\hat{A}|^{2}+|k(D^{2}-k^{2})\hat{A}|^{2}]-k^{2}\unicode[STIX]{x1D701}(x)[|D\hat{A}|^{2}+|k\hat{A}|^{2}])\,\text{d}x=0.\qquad\end{eqnarray}$$ The terms corresponding to 
${\mathcal{X}}$ and 
${\mathcal{Y}}$ are strictly positive as 
$\bar{\unicode[STIX]{x1D70C}}_{x}$ is positive throughout the domain. So we now turn our attention to 
$\unicode[STIX]{x1D701}(x)$ as the remaining terms are positive. We have
$$\begin{eqnarray}\unicode[STIX]{x1D701}(x)=\bar{A}_{x}(1-\unicode[STIX]{x1D70F})-\bar{\unicode[STIX]{x1D70C}}_{x},\end{eqnarray}$$ where 
$\bar{A}_{x}$ is given by
$$\begin{eqnarray}\bar{A}_{x}=\frac{-\exp [-x^{2}/4t]}{2\sqrt{\unicode[STIX]{x03C0}t}}.\end{eqnarray}$$Using the above, we have
$$\begin{eqnarray}\unicode[STIX]{x1D701}(x)=\frac{\exp [-x^{2}/4t]\unicode[STIX]{x1D70F}}{2\sqrt{\unicode[STIX]{x03C0}t}}-\frac{\exp [-x^{2}/4t\unicode[STIX]{x1D70F}]R}{2\sqrt{\unicode[STIX]{x03C0}t}\sqrt{\unicode[STIX]{x1D70F}}}.\end{eqnarray}$$The second derivative of the base-state density profile is given by
$$\begin{eqnarray}\bar{\unicode[STIX]{x1D70C}}_{xx}=\frac{\exp [-x^{2}/4t]x}{4\sqrt{\unicode[STIX]{x03C0}}t^{3/2}}-\frac{\exp [-x^{2}/4t\unicode[STIX]{x1D70F}]Rx}{4\sqrt{\unicode[STIX]{x03C0}}t^{3/2}\unicode[STIX]{x1D70F}^{3/2}}.\end{eqnarray}$$ Using this, the equation for 
$\unicode[STIX]{x1D701}(x)$ simplifies as
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D701}(x)\frac{x}{2\unicode[STIX]{x1D70F}t}=\bar{\unicode[STIX]{x1D70C}}_{xx}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D701}(x)=\frac{2\unicode[STIX]{x1D70F}t}{x}\bar{\unicode[STIX]{x1D70C}}_{xx}. & \displaystyle\end{eqnarray}$$ If 
$\unicode[STIX]{x1D701}(x)$ is negative throughout the domain, from (4.18) 
$\unicode[STIX]{x1D70E}_{r}<0$ and the system is linearly stable. For instability, 
$\unicode[STIX]{x1D701}(x)$ must be greater than zero within some region of the domain. Thus we have that if 
$\bar{\unicode[STIX]{x1D70C}}_{xx}>0$ for 
$x<0$ (above the initial interface) and if 
$\bar{\unicode[STIX]{x1D70C}}_{xx}<0$ for 
$x>0$ (below the initial interface), then 
$\unicode[STIX]{x1D70E}_{r}<0$ and the system is linearly stable. The instability develops on either side of the interface, which implies that above the interface (
$x<0$) there must be region where 
$\bar{\unicode[STIX]{x1D70C}}_{xx}<0$ and, correspondingly, below the interface (
$x>0$) there must be a region where 
$\bar{\unicode[STIX]{x1D70C}}_{xx}>0$. Also, it is important to note that, at 
$t=0$, 
$\unicode[STIX]{x1D701}(x)=0$, and from this it follows that 
$\unicode[STIX]{x1D70E}_{r}<0$. So the flows where 
$\bar{\unicode[STIX]{x1D70C}}_{x}>0$ throughout the domain are initially linearly stable with the instability developing in time.
The density profile has an inflection point at the interface corresponding to 
$x=0$, as we can see from the expression (4.22) that 
$\bar{\unicode[STIX]{x1D70C}}_{xx}=0$. Instability, however, requires that 
$\bar{\unicode[STIX]{x1D70C}}_{xx}<0$ above the interface (
$x<0$) and 
$\bar{\unicode[STIX]{x1D70C}}_{xx}>0$ below the interface (
$x>0$). Since 
$\bar{\unicode[STIX]{x1D70C}}_{x}>0$ throughout the domain, above and below the interface, there has to be an additional point such that 
$\bar{\unicode[STIX]{x1D70C}}_{xx}=0$ on either side of the domain such that 
$\bar{\unicode[STIX]{x1D70C}}_{xx}$ goes from positive to negative above the initial interface, and from negative to positive below the interface. Otherwise it follows from (4.18) that 
$\unicode[STIX]{x1D70E}_{r}<0$. Thus the boundary of stability when 
$\bar{\unicode[STIX]{x1D70C}}_{x}>0$ throughout the domain can be obtained by using 
$\bar{\unicode[STIX]{x1D70C}}_{xx}=0$ (for 
$x\neq 0$) and is given by
$$\begin{eqnarray}R>\unicode[STIX]{x1D70F}^{3/2}\exp \left[\left(\frac{x^{2}}{4t}\right)\left(\frac{1}{\unicode[STIX]{x1D70F}}-1\right)\right],\end{eqnarray}$$ where all the eigenvalues in this regime are strictly real. The above function has a minimum value at 
$x=0$, which gives the limiting boundary for instability beyond which the flows are linearly stable:
$$\begin{eqnarray}R=\unicode[STIX]{x1D70F}^{3/2}.\end{eqnarray}$$The fact that all the eigenvalues are strictly real in such regimes has been numerically observed by Trevelyan et al. (Reference Trevelyan, Almarcha and De Wit2011, p. 52), which follows as a direct consequence from the discussion presented here. A similar asymptotic limit had been shown by Huppert & Manins (Reference Huppert and Manins1973) by considering linear base states, and is consistent with the numerical stability analysis carried out by Trevelyan et al. (Reference Trevelyan, Almarcha and De Wit2011). The result presented here when the density profiles are monotonically increasing downwards is qualitatively similar to the Rayleigh criterion observed in the case of shear flows (Lord Rayleigh Reference Rayleigh1880; Taylor Reference Taylor1915).
Figure 3. (a–d) Typical base-state density profiles 
$\bar{\unicode[STIX]{x1D70C}}(x,t)$ at 
$t>0$ for different combinations of 
$(R,\unicode[STIX]{x1D70F})$; (e–h) the corresponding first derivatives 
$\bar{\unicode[STIX]{x1D70C}}_{x}(x,t)$ (dotted line) and second derivatives 
$\bar{\unicode[STIX]{x1D70C}}_{xx}(x,t)$ (continuous line). Parameter settings: (a,e) 
$\unicode[STIX]{x1D70F}=8.0$, 
$R=0.50$; (b,f) 
$\unicode[STIX]{x1D70F}=0.80$, 
$R=0.25$; (c,g) 
$\unicode[STIX]{x1D70F}=0.25$, 
$R=0.75$; and (d,h) 
$\unicode[STIX]{x1D70F}=0.10$, 
$R=1.25$.
We shall now briefly overview the various base-state density profile configurations in the (
$R,\unicode[STIX]{x1D70F}$) plane to summarise the results. Figure 3(a–c) shows some typical base-state density profiles for 
$R<1$ (RT instability) whereas figure 3(d) corresponds to a DLC scenario. The density profiles have a region of 
$x$ where 
$\bar{\unicode[STIX]{x1D70C}}_{x}<0$, which can be seen in figure 3(e,f) and they are linearly unstable to stationary or oscillatory disturbances. It is also interesting to note that the regions corresponding to 
$\bar{\unicode[STIX]{x1D70C}}_{x}<0$ include the initial interface 
$x=0$ or are close to it for the scenarios corresponding to RT instabilities, which is also where the eigenfunctions are localised (Trevelyan et al. Reference Trevelyan, Almarcha and De Wit2011). The second derivatives 
$\bar{\unicode[STIX]{x1D70C}}_{xx}$ may have both positive and negative regions on either side of the interface or could be either positive or negative, as can be seen in figure 3(e,f).
Base-state profiles that are monotonically increasing are shown in figure 4(a–d), with their corresponding 
$\bar{\unicode[STIX]{x1D70C}}_{x}$ and 
$\bar{\unicode[STIX]{x1D70C}}_{xx}$ in figure 4(e,f). In all these cases, 
$\bar{\unicode[STIX]{x1D70C}}_{x}>0$ throughout the domain. In figure 4(e,f) we can see that for 
$x<0$ there is a region where 
$\bar{\unicode[STIX]{x1D70C}}_{xx}<0$ and a region where 
$\bar{\unicode[STIX]{x1D70C}}_{xx}>0$ for 
$x>0$, or in other words 
$\bar{\unicode[STIX]{x1D70C}}_{xx}=0$ at some location on either side of the interface. In such regimes, maximum growth rates are real and instability arises as stationary modes. A case corresponding to the asymptotic stability limit is shown in figure 4(g) beyond which 
$\bar{\unicode[STIX]{x1D70C}}_{xx}>0$ above the initial interface and 
$\bar{\unicode[STIX]{x1D70C}}_{xx}<0$ below (see figure 4h).
Figure 4. (a–d) Base-state density profiles that are monotonically increasing; (e–h) corresponding 
$\bar{\unicode[STIX]{x1D70C}}_{x}(x,t)$ (dotted line) and 
$\bar{\unicode[STIX]{x1D70C}}_{xx}(x,t)$ (continuous line). Parameter settings: 
$\unicode[STIX]{x1D70F}=4$ and (a,e) 
$R=2.0$, (b,f) 
$R=4.0$, (c,g) 
$R=8.0$, and (d,h) 
$R=15.0$.
4.1 Single-species RT instability
The single-species RT instability has recently received renewed attention due to its potential role in geological carbon dioxide storage (Huppert & Neufeld Reference Huppert and Neufeld2014; Slim Reference Slim2014; Sardina et al. Reference Sardina, Brandt, Boffetta and Mazzino2018; De Paoli et al. Reference De Paoli, Giurgiu, Zonta and Soldati2019a). Here we shall briefly discuss the special case of a single-species system that amounts to considering 
$\unicode[STIX]{x1D70F}=1$. We have from (3.14d)
$$\begin{eqnarray}\left[\frac{\unicode[STIX]{x1D70E}}{k^{2}}{\mathcal{L}}-\frac{1}{k^{2}}{\mathcal{L}}^{2}-\bar{\unicode[STIX]{x1D70C}}_{x}\right]\hat{\unicode[STIX]{x1D713}}=0.\end{eqnarray}$$The above can be rewritten as
$$\begin{eqnarray}-\unicode[STIX]{x1D70E}(D^{2}-k^{2})\hat{\unicode[STIX]{x1D713}}+(D^{2}-k^{2})^{2}\hat{\unicode[STIX]{x1D713}}=-\bar{\unicode[STIX]{x1D70C}}_{x}k^{2}\hat{\unicode[STIX]{x1D713}}\end{eqnarray}$$ because 
${\mathcal{L}}=(D^{2}-k^{2})$. We have
$$\begin{eqnarray}(D^{2}-k^{2})(D^{2}-k^{2}-\unicode[STIX]{x1D70E})\hat{\unicode[STIX]{x1D713}}=-\bar{\unicode[STIX]{x1D70C}}_{x}k^{2}\hat{\unicode[STIX]{x1D713}}.\end{eqnarray}$$ Letting 
$(D^{2}-k^{2})\hat{\unicode[STIX]{x1D713}}={\mathcal{G}}$ we rewrite the above equation as
$$\begin{eqnarray}(D^{2}-k^{2}-\unicode[STIX]{x1D70E}){\mathcal{G}}=-\bar{\unicode[STIX]{x1D70C}}_{x}k^{2}\hat{\unicode[STIX]{x1D713}}.\end{eqnarray}$$ Multiplying the above equation by 
$\hat{\unicode[STIX]{x1D713}}^{\ast }$ (the complex conjugate of 
$\hat{\unicode[STIX]{x1D713}}$) and integrating over the range of 
$x$ (the vertical extent of the domain) we get
$$\begin{eqnarray}\int _{-\infty }^{\infty }\hat{\unicode[STIX]{x1D713}}^{\ast }(D^{2}-k^{2}-\unicode[STIX]{x1D70E}){\mathcal{G}}\,\text{d}x=-\int _{-\infty }^{\infty }\hat{\unicode[STIX]{x1D713}}^{\ast }\bar{\unicode[STIX]{x1D70C}}_{x}k^{2}\hat{\unicode[STIX]{x1D713}}\,\text{d}x.\end{eqnarray}$$Considering the left-hand side of (4.31) we have
$$\begin{eqnarray}\int _{-\infty }^{\infty }\hat{\unicode[STIX]{x1D713}}^{\ast }(D^{2}-k^{2}-\unicode[STIX]{x1D70E}){\mathcal{G}}\,\text{d}x=\int _{-\infty }^{\infty }\hat{\unicode[STIX]{x1D713}}^{\ast }(D^{2}-k^{2}){\mathcal{G}}\,\text{d}x-\int _{-\infty }^{\infty }\unicode[STIX]{x1D70E}\hat{\unicode[STIX]{x1D713}}^{\ast }{\mathcal{G}}\,\text{d}x.\end{eqnarray}$$ As 
${\mathcal{G}}=(D^{2}-k^{2})\hat{\unicode[STIX]{x1D713}}$, equation (4.32) simplifies to
$$\begin{eqnarray}\displaystyle & & \displaystyle \int _{-\infty }^{\infty }|{\mathcal{G}}|^{2}\,\text{d}x-\int _{-\infty }^{\infty }\unicode[STIX]{x1D70E}\hat{\unicode[STIX]{x1D713}}^{\ast }(D^{2}-k^{2})\hat{\unicode[STIX]{x1D713}}\,\text{d}x\nonumber\\ \displaystyle & & \displaystyle \quad =\int _{-\infty }^{\infty }|{\mathcal{G}}|^{2}\,\text{d}x-\int _{-\infty }^{\infty }\unicode[STIX]{x1D70E}\hat{\unicode[STIX]{x1D713}}^{\ast }D^{2}\hat{\unicode[STIX]{x1D713}}\,\text{d}x+\int _{-\infty }^{\infty }\unicode[STIX]{x1D70E}\hat{\unicode[STIX]{x1D713}}^{\ast }\hat{\unicode[STIX]{x1D713}}k^{2}\,\text{d}x.\end{eqnarray}$$Using this in (4.31) gives
$$\begin{eqnarray}\int _{-\infty }^{\infty }[|{\mathcal{G}}|^{2}+\unicode[STIX]{x1D70E}(|D\hat{\unicode[STIX]{x1D713}}|^{2}+k^{2}|\hat{\unicode[STIX]{x1D713}}|^{2})]\,\text{d}x+k^{2}\int _{-\infty }^{\infty }\bar{\unicode[STIX]{x1D70C}}_{x}|\hat{\unicode[STIX]{x1D713}}|^{2}\,\text{d}x=0.\end{eqnarray}$$The real and imaginary parts of this equation must vanish separately; and the vanishing of the imaginary part gives
$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{i}\int _{-\infty }^{\infty }[|D\hat{\unicode[STIX]{x1D713}}|^{2}+k^{2}|\hat{\unicode[STIX]{x1D713}}|^{2}]\,\text{d}x=0,\end{eqnarray}$$ where 
$\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D70E}_{r}+\text{i}\unicode[STIX]{x1D70E}_{i}$. But the quantity inside the brackets is positive definite. Hence
$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{i}=0.\end{eqnarray}$$ This establishes that 
$\unicode[STIX]{x1D70E}$ is real. So we have
$$\begin{eqnarray}\displaystyle & \displaystyle \int _{-\infty }^{\infty }[|{\mathcal{G}}|^{2}+\unicode[STIX]{x1D70E}_{r}(|D\hat{\unicode[STIX]{x1D713}}|^{2}+k^{2}|\hat{\unicode[STIX]{x1D713}}|^{2})]\,\text{d}x+k^{2}\int _{-\infty }^{\infty }\bar{\unicode[STIX]{x1D70C}}_{x}|\hat{\unicode[STIX]{x1D713}}|^{2}\,\text{d}x=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}_{r}=\frac{\displaystyle -k^{2}\int _{-\infty }^{\infty }\bar{\unicode[STIX]{x1D70C}}_{x}|\hat{\unicode[STIX]{x1D713}}|^{2}\,\text{d}x-|{\mathcal{G}}|^{2}\,\text{d}x}{\displaystyle \int _{-\infty }^{\infty }(|D\hat{\unicode[STIX]{x1D713}}|^{2}+k^{2}|\hat{\unicode[STIX]{x1D713}}|^{2})\,\text{d}x}. & \displaystyle\end{eqnarray}$$ From the above, we can see that, if 
$\bar{\unicode[STIX]{x1D70C}}_{x}$ is positive throughout the domain, then 
$\unicode[STIX]{x1D70E}_{r}<0$ and the flow is linearly stable. Only when there is a region in the base flow profile where 
$\bar{\unicode[STIX]{x1D70C}}_{x}$ is negative, can we have 
$\unicode[STIX]{x1D70E}_{r}>0$, which is a necessary and sufficient condition for instability in single-species RT flows. Also, the eigenvalues are real and instability arises as non-oscillatory modes. In these flows, the onset of instability is marked by a transition from the background state to another steady state with the principle of exchange of stabilities being valid, which holds if all non-decaying disturbances are non-oscillatory in time (Chandrasekhar Reference Chandrasekhar1961; Davis Reference Davis1969). It is also interesting to note that Rayleigh–Bénard convection and Taylor–Couette flow are two other well-known problems in which 
$\unicode[STIX]{x1D70E}$ is real (Chandrasekhar Reference Chandrasekhar1961; Davis Reference Davis1969).
4.2 Immobile species B: 
$\unicode[STIX]{x1D70F}=0$
Figure 5. Plots of (a,b) 
$\bar{A}_{x}$ (dotted line) and 
$\bar{B}_{x}$ (continuous line) at 
$t>0$ for 
$R=0.50$. Parameter settings: (a) 
$\unicode[STIX]{x1D70F}=0.10$, and (b) 
$\unicode[STIX]{x1D70F}=0.001$.
We discuss the case 
$\unicode[STIX]{x1D70F}=0$, which physically corresponds to an immobile species B. It also provides a good approximation to the case when A diffuses much faster than B. For 
$\unicode[STIX]{x1D70F}=0$ from (4.14) we have
$$\begin{eqnarray}\int _{-\infty }^{\infty }\unicode[STIX]{x1D70E}_{i}(2\unicode[STIX]{x1D70E}_{r}[|D\hat{A}|^{2}+|k\hat{A}|^{2}]+[(D^{2}-k^{2})\hat{A}]^{2}+\bar{\unicode[STIX]{x1D70C}}_{x}|k\hat{A}|^{2})\,\text{d}x=0.\end{eqnarray}$$ As noted earlier, if 
$\bar{\unicode[STIX]{x1D70C}}_{x}$ is positive throughout the domain, then the system could be unstable only to purely real eigenvalues with 
$\unicode[STIX]{x1D70E}_{i}=0$. Otherwise, they could be complex, and oscillatory instabilities may be present. From (4.15) we can obtain the limit of 
$R$ above which 
$\bar{\unicode[STIX]{x1D70C}}_{x}$ is strictly positive, which is given by
$$\begin{eqnarray}R_{c}=\exp \left[-\frac{x^{2}}{4t}\left(1-\frac{1}{\unicode[STIX]{x1D70F}}\right)\right]\sqrt{\unicode[STIX]{x1D70F}}.\end{eqnarray}$$ The gradient at 
$x=0$ is given by
$$\begin{eqnarray}\bar{\unicode[STIX]{x1D70C}}_{x}|_{x=0}=\frac{1}{2\sqrt{\unicode[STIX]{x03C0}t}}\left(\frac{R}{\sqrt{\unicode[STIX]{x1D70F}}}-1\right),\end{eqnarray}$$ which indicates that, as 
$\unicode[STIX]{x1D70F}\rightarrow 0$, the function goes to 
$\infty$ at 
$x=0$. When 
$\unicode[STIX]{x1D70F}\rightarrow 0$, 
$\bar{B}_{x}=\exp [-x^{2}/4t\unicode[STIX]{x1D70F}]/(2\sqrt{\unicode[STIX]{x03C0}t\unicode[STIX]{x1D70F}})$, and can be treated as a Dirac delta function 
$\unicode[STIX]{x1D6FF}(x)$ as
$$\begin{eqnarray}\bar{B}_{x}=\bar{B}_{x}[\unicode[STIX]{x1D6FF}(x)\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FF},0}+(1-\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FF},0})],\end{eqnarray}$$ where 
$\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FF},0}$ is the Kronecker delta. We have 
$\bar{B}_{x}\unicode[STIX]{x1D6FF}(x)$ when 
$\unicode[STIX]{x1D70F}=0$, and 
$\bar{B}_{x}$ for 
$\unicode[STIX]{x1D70F}\neq 0$. So (4.39) reads as
$$\begin{eqnarray}\int _{-\infty }^{\infty }\unicode[STIX]{x1D70E}_{i}(2\unicode[STIX]{x1D70E}_{r}[|D\hat{A}|^{2}+|k\hat{A}|^{2}]+[(D^{2}-k^{2})\hat{A}]^{2}+[\bar{A}_{x}+R\bar{B}_{x}[\unicode[STIX]{x1D6FF}(x)\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FF},0}+(1-\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FF},0})]]|k\hat{A}|^{2})\,\text{d}x=0.\end{eqnarray}$$ For 
$\unicode[STIX]{x1D70F}=0$ we have
$$\begin{eqnarray}\int _{-\infty }^{\infty }\unicode[STIX]{x1D70E}_{i}(2\unicode[STIX]{x1D70E}_{r}[|D\hat{A}|^{2}+|k\hat{A}|^{2}]+[(D^{2}-k^{2})\hat{A}]^{2}+[\bar{A}_{x}+R\bar{B}_{x}\unicode[STIX]{x1D6FF}(x)]|k\hat{A}|^{2})\,\text{d}x=0.\end{eqnarray}$$ As 
$\bar{A}_{x}<0$ throughout the domain (4.20), the system could be unstable to both stationary or oscillatory disturbances. Some representative profiles of 
$\bar{A}_{x}$ and 
$\bar{B}_{x}$ illustrating 
$\unicode[STIX]{x1D70F}\rightarrow 0$ are shown in figure 5.
4.3 Discussion
In this section, we briefly discuss some of the implications of the results obtained in this study in the context of porous media flows. For the scenarios corresponding to a single-species system (
$\unicode[STIX]{x1D70F}=1$), and in the double-diffusive regime where density profiles are monotonically increasing downwards (
$1\leqslant \unicode[STIX]{x1D70F}\leqslant R^{2}$), the instability arises as stationary modes. Whereas in the remaining 
$(R,\unicode[STIX]{x1D70F})$ parameter space bounded by the curve 
$R=\unicode[STIX]{x1D70F}^{3/2}$, both stationary and oscillatory modes may destabilise the system. Oscillatory instabilities generally occur when the density gradients are altered due to both energy and mass transfer across a fluid layer provided both the mechanisms work in opposition: one of the gradients must induce a stabilising effect and the other must encourage instability. For instance, the presence of such oscillatory modes can influence the solidification process of a binary mixture where a porous mushy layer is formed via double-diffusive effects (Chen, Lu & Yang Reference Chen, Lu and Yang1994; Anderson & Worster Reference Anderson and Worster1996). The fact that the instability arises as stationary modes in the double-diffusive regime where the base-state density profiles are monotonically increasing downwards such that the density gradient is positive throughout the domain may in turn be used in such solidification processes where oscillatory instabilities are non-desirable.
The density profiles in porous media systems may also be altered due to changes in temperature and composition as a result of chemical reactions (De Wit Reference De Wit2020). In the present study, as the stability boundary and the nature of instabilities were determined based on the characteristics of the density profile, the obtained results may be used in the context of buoyancy-driven instabilities in reactive systems. Such chemohydrodynamic flows may be used to tune the spatiotemporal distribution of the species, thereby providing control over pattern formation in these systems.
5 Concluding remarks
The interface separating two miscible solutions containing different species in porous media and Hele-Shaw cells can destabilise in the gravity field due to an unstable density stratification, or due to differential diffusion of the species. The goal of this study was to obtain the conditions for buoyancy-driven instabilities to arise from the base-state density profiles using analytical methods. If there is an interval in the spatial domain where the first derivative of the base-state density profile is negative (
$\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D70C}}/\unicode[STIX]{x2202}x<0$), the flows are linearly unstable to stationary or oscillatory disturbances. This is a sufficient condition for instability in these flows. When the base-state density profiles are strictly monotonically increasing downwards (
$1\leqslant \unicode[STIX]{x1D70F}\leqslant R^{2}$, where 
$R$ is the buoyancy ratio and 
$\unicode[STIX]{x1D70F}$ is the diffusion coefficient ratio), with the first derivative of the density profile positive throughout the domain, instability develops if and only if the second derivative changes its sign on either side of the initial interface, or alternatively there is a point where the second derivative of the base-state density profile is zero on either side of the initial interface. In such regimes the system is unstable to uniquely stationary disturbances (purely real growth rates). This is a necessary and sufficient condition for instability when the first derivative of the base-state density profile is positive throughout the domain. A necessary condition for oscillatory instabilities in these flows is that there is an interval in the spatial domain where the first derivative of the base-state density profile is negative. The asymptotic neutral stability curve (
$R=\unicode[STIX]{x1D70F}^{3/2}$) obtained in this analysis is consistent with earlier studies where linear constant base-state profiles were used, and with numerical linear stability analyses. The fact that single-species RT instabilities arise as stationary modes also follows from the study presented in this paper. The analysis could be extended to another member of the buoyancy-driven instabilities in porous media, namely the viscous fingering instability, which is being currently investigated by the author.
Acknowledgements
The author would like to thank Dr A. De Wit and Dr B. Knaepen for stimulating discussions at the beginning of this work. A special word of gratitude goes to Dr A. Ranjbar for insightful feedback on some of the tools used in this study. The author gratefully acknowledges the constructive feedback from the anonymous referees during the reviewing process, which greatly helped in improving the manuscript.
Declaration of interests
The author reports no conflict of interest.
Appendix A. Reduction of the integrals $\mathbb{I}_{2}$ and $\mathbb{I}_{3}$
The expression for 
$\mathbb{I}_{2}$ reads as
$$\begin{eqnarray}\displaystyle \mathbb{I}_{2} & = & \displaystyle \int _{-\infty }^{\infty }-\unicode[STIX]{x1D70E}\hat{A}^{\ast }[(\unicode[STIX]{x1D70F}+1){\mathcal{L}}^{2}+\bar{\unicode[STIX]{x1D70C}}_{x}k^{2}]\hat{A}\,\text{d}x\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \int _{-\infty }^{\infty }-\unicode[STIX]{x1D70E}\hat{A}^{\ast }(\unicode[STIX]{x1D70F}+1)(D^{2}-k^{2})^{2}\hat{A}\,\text{d}x-\int _{-\infty }^{\infty }\unicode[STIX]{x1D70E}\hat{A}^{\ast }\bar{\unicode[STIX]{x1D70C}}_{x}k^{2}\hat{A}\,\text{d}x.\end{eqnarray}$$The first of the two integrals on the right-hand side of (A 2) can be simplified as
$$\begin{eqnarray}\displaystyle & & \displaystyle \int _{-\infty }^{\infty }-\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D70F}+1)\hat{A}^{\ast }(D^{2}-k^{2})^{2}\hat{A}\,\text{d}x\nonumber\\ \displaystyle & & \displaystyle \quad =\int _{-\infty }^{\infty }-\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D70F}+1)\hat{A}^{\ast }(D^{2}-k^{2})(D^{2}-k^{2})\hat{A}\,\text{d}x\nonumber\\ \displaystyle & & \displaystyle \quad =\int _{-\infty }^{\infty }-\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D70F}+1)\hat{A}^{\ast }D^{2}[(D^{2}-k^{2})\hat{A}]\,\text{d}x+\int _{-\infty }^{\infty }\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D70F}+1)k^{2}\hat{A}^{\ast }(D^{2}-k^{2})\hat{A}\,\text{d}x\nonumber\\ \displaystyle & & \displaystyle \quad =-\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D70F}+1)\left[\int _{-\infty }^{\infty }\hat{A}^{\ast }D^{2}[(D^{2}-k^{2})\hat{A}]\,\text{d}x-k^{2}\int _{-\infty }^{\infty }\hat{A}^{\ast }(D^{2}-k^{2})\hat{A}\,\text{d}x\right].\end{eqnarray}$$After two successive integrations by parts, the first of the two integrals on the right-hand side of the foregoing equation becomes
$$\begin{eqnarray}\displaystyle \int _{-\infty }^{\infty }\hat{A}^{\ast }D^{2}[(D^{2}-k^{2})\hat{A}]\,\text{d}x & = & \displaystyle [\hat{A}^{\ast }D[(D^{2}-k^{2})\hat{A}]]_{-\infty }^{\infty }-\int _{-\infty }^{\infty }D\hat{A}^{\ast }D[(D^{2}-k^{2})\hat{A}]\,\text{d}x\nonumber\\ \displaystyle & = & \displaystyle -\left([D\hat{A}^{\ast }[(D^{2}-k^{2})\hat{A}]]_{-\infty }^{\infty }-\int _{-\infty }^{\infty }D^{2}\hat{A}^{\ast }[(D^{2}-k^{2})\hat{A}]\,\text{d}x\right)\nonumber\\ \displaystyle & = & \displaystyle \int _{-\infty }^{\infty }D^{2}\hat{A}^{\ast }[(D^{2}-k^{2})\hat{A}]\,\text{d}x,\end{eqnarray}$$with the integrated part vanishing on account of the boundary conditions (3.11). Using the above expression in (A 3) we get
$$\begin{eqnarray}\displaystyle & & \displaystyle \int _{-\infty }^{\infty }-\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D70F}+1)\hat{A}^{\ast }(D^{2}-k^{2})^{2}\hat{A}\,\text{d}x\nonumber\\ \displaystyle & & \displaystyle \quad =-\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D70F}+1)\left[\int _{-\infty }^{\infty }D^{2}\hat{A}^{\ast }[(D^{2}-k^{2})\hat{A}]\,\text{d}x-\cdots -k^{2}\int _{-\infty }^{\infty }\hat{A}^{\ast }(D^{2}-k^{2})\hat{A}\,\text{d}x\right]\nonumber\\ \displaystyle & & \displaystyle \quad =-\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D70F}+1)\int _{-\infty }^{\infty }[(D^{2}-k^{2})\hat{A}]^{2}\,\text{d}x.\end{eqnarray}$$ So expression 
$\mathbb{I}_{2}$ reads as
$$\begin{eqnarray}\mathbb{I}_{2}=\int _{-\infty }^{\infty }-\unicode[STIX]{x1D70E}[(\unicode[STIX]{x1D70F}+1)[(D^{2}-k^{2})\hat{A}]^{2}+\bar{\unicode[STIX]{x1D70C}}_{x}|k\hat{A}|^{2}]\,\text{d}x.\end{eqnarray}$$ Considering next 
$\mathbb{I}_{3}$ (4.4) we have
$$\begin{eqnarray}\displaystyle \mathbb{I}_{3} & = & \displaystyle \int _{-\infty }^{\infty }\hat{A}^{\ast }(\unicode[STIX]{x1D70F}{\mathcal{L}}^{3}\hat{A}-k^{2}\unicode[STIX]{x1D701}(x){\mathcal{L}}\hat{A})\,\text{d}x\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \int _{-\infty }^{\infty }\hat{A}^{\ast }(\unicode[STIX]{x1D70F}[D^{2}-k^{2}]^{2}[D^{2}-k^{2}]\hat{A}-k^{2}\unicode[STIX]{x1D701}(x)(D^{2}-k^{2})\hat{A})\,\text{d}x\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \int _{-\infty }^{\infty }\hat{A}^{\ast }\unicode[STIX]{x1D70F}(D^{2}-k^{2})^{2}(D^{2}-k^{2})\hat{A}\,\text{d}x-\int _{-\infty }^{\infty }\hat{A}^{\ast }k^{2}(D^{2}-k^{2})\unicode[STIX]{x1D701}(x)\hat{A}\,\text{d}x.\end{eqnarray}$$ Let the two integrals in the above equation on the right-hand side be 
$\mathbb{I}_{3a}$ and 
$\mathbb{I}_{3b}$, respectively. From the first integral we have
$$\begin{eqnarray}\displaystyle \mathbb{I}_{3a} & = & \displaystyle \int _{-\infty }^{\infty }\unicode[STIX]{x1D70F}\hat{A}^{\ast }(D^{2}-k^{2})^{2}(D^{2}-k^{2})\hat{A}\,\text{d}x\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \int _{-\infty }^{\infty }\unicode[STIX]{x1D70F}\hat{A}^{\ast }[(D^{2}-k^{2})^{2}]D^{2}\hat{A}\,\text{d}x-\int _{-\infty }^{\infty }\unicode[STIX]{x1D70F}k^{2}\hat{A}^{\ast }(D^{2}-k^{2})^{2}\hat{A}\,\text{d}x.\end{eqnarray}$$Simplifying the first integral on the right-hand side by successive integrations by parts we have
$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \int _{-\infty }^{\infty }\unicode[STIX]{x1D70F}\hat{A}^{\ast }[(D^{2}-k^{2})]^{2}D^{2}\hat{A}\,\text{d}x\nonumber\\ \displaystyle & & \displaystyle \quad =\int _{-\infty }^{\infty }\unicode[STIX]{x1D70F}\hat{A}^{\ast }D^{2}(D^{2}-k^{2})^{2}\hat{A}\,\text{d}x\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle \quad =\displaystyle [\unicode[STIX]{x1D70F}\hat{A}^{\ast }D(D^{2}-k^{2})^{2}\hat{A}]_{-\infty }^{\infty }-\int _{-\infty }^{\infty }\unicode[STIX]{x1D70F}D\hat{A}^{\ast }D(D^{2}-k^{2})^{2}\hat{A}\,\text{d}x\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle \quad =\displaystyle -\left([\unicode[STIX]{x1D70F}D\hat{A}^{\ast }(D^{2}-k^{2})^{2}\hat{A}]_{-\infty }^{\infty }-\int _{-\infty }^{\infty }\unicode[STIX]{x1D70F}D^{2}\hat{A}^{\ast }(D^{2}-k^{2})^{2}\hat{A}\,\text{d}x\right)\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle \quad =\displaystyle \int _{-\infty }^{\infty }\unicode[STIX]{x1D70F}D^{2}\hat{A}^{\ast }(D^{2}-k^{2})^{2}\hat{A}\,\text{d}x.\end{eqnarray}$$ The expression 
$\mathbb{I}_{3a}$ reads as
$$\begin{eqnarray}\displaystyle \mathbb{I}_{3a} & = & \displaystyle \int _{-\infty }^{\infty }\unicode[STIX]{x1D70F}D^{2}\hat{A}^{\ast }(D^{2}-k^{2})^{2}\hat{A}\,\text{d}x-\int _{-\infty }^{\infty }\unicode[STIX]{x1D70F}k^{2}\hat{A}^{\ast }(D^{2}-k^{2})^{2}\hat{A}\,\text{d}x\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \int _{-\infty }^{\infty }\unicode[STIX]{x1D70F}(D^{2}-k^{2})\hat{A}^{\ast }(D^{2}-k^{2})^{2}\hat{A}\,\text{d}x.\end{eqnarray}$$ Let 
${\hat{W}}=(D^{2}-k^{2})\hat{A}$. The complex conjugate of 
${\hat{W}}$ is given by 
${\hat{W}}^{\ast }=(D^{2}-k^{2})\hat{A}^{\ast }$. The expression 
$\mathbb{I}_{3a}$ can be rewritten and simplified as
$$\begin{eqnarray}\displaystyle \mathbb{I}_{3a} & = & \displaystyle \int _{-\infty }^{\infty }\unicode[STIX]{x1D70F}{\hat{W}}^{\ast }(D^{2}-k^{2}){\hat{W}}\,\text{d}x\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \int _{-\infty }^{\infty }\unicode[STIX]{x1D70F}{\hat{W}}^{\ast }D^{2}{\hat{W}}\,\text{d}x-\int _{-\infty }^{\infty }\unicode[STIX]{x1D70F}k^{2}{\hat{W}}^{\ast }{\hat{W}}\,\text{d}x\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \int _{-\infty }^{\infty }-\unicode[STIX]{x1D70F}|D{\hat{W}}|^{2}\,\text{d}x-\int _{-\infty }^{\infty }\unicode[STIX]{x1D70F}|k{\hat{W}}|^{2}\,\text{d}x\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \int _{-\infty }^{\infty }-\unicode[STIX]{x1D70F}[|D{\hat{W}}|^{2}+|k{\hat{W}}|^{2}]\,\text{d}x\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \int _{-\infty }^{\infty }-\unicode[STIX]{x1D70F}[|D(D^{2}-k^{2})\hat{A}|^{2}+|k(D^{2}-k^{2})\hat{A}|^{2}]\,\text{d}x.\end{eqnarray}$$ We now turn our attention to 
$\mathbb{I}_{3b}$, which can be transformed in the manner
$$\begin{eqnarray}\displaystyle \mathbb{I}_{3b} & = & \displaystyle -\int _{-\infty }^{\infty }\hat{A}^{\ast }k^{2}(D^{2}-k^{2})\unicode[STIX]{x1D701}(x)\hat{A}\,\text{d}x\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle -\int _{-\infty }^{\infty }k^{2}\unicode[STIX]{x1D701}(x)\hat{A}^{\ast }D^{2}\hat{A}\,\text{d}x+\int _{-\infty }^{\infty }k^{4}\unicode[STIX]{x1D701}(x)\hat{A}^{\ast }\hat{A}\,\text{d}x\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \int _{-\infty }^{\infty }k^{2}\unicode[STIX]{x1D701}(x)|D\hat{A}|^{2}\,\text{d}x+\int _{-\infty }^{\infty }\unicode[STIX]{x1D701}(x)k^{4}|\hat{A}|^{2}\,\text{d}x\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \int _{-\infty }^{\infty }k^{2}\unicode[STIX]{x1D701}(x)[|D\hat{A}|^{2}+|k\hat{A}|^{2}]\,\text{d}x.\end{eqnarray}$$ So 
$\mathbb{I}_{3}$ is given by
$$\begin{eqnarray}\mathbb{I}_{3}=\int _{-\infty }^{\infty }-\unicode[STIX]{x1D70F}[|D(D^{2}-k^{2})\hat{A}|^{2}+|k(D^{2}-k^{2})\hat{A}|^{2}]\,\text{d}x+\int _{-\infty }^{\infty }k^{2}\unicode[STIX]{x1D701}(x)[|D\hat{A}|^{2}+|k\hat{A}|^{2}]\,\text{d}x.\end{eqnarray}$$
