Stability and asymptotic analysis for instationary gas transport via relative energy estimates

In this paper, we investigate the stability of solutions to the nonlinear system (1)–(3) with respect to perturbations of the initial and boundary values as well as the model parameters $\varepsilon$ and $\gamma$. Our analysis is done under the following structural assumptions:

• (A1)

  The pressure potential $P:\mathbb{R}_{+}\to\mathbb{R}$ is smooth and strictly convex and there exist positive constants $\underaccent{\bar}{\rho},\bar{\rho},\bar{w}$ and $\bar{\varepsilon}$ such that

  $\displaystyle\rho P^{\prime\prime}(\rho)\geq 4\bar{\varepsilon}^{2}|\bar{w}|^{2}\qquad\forall\underaccent{\bar}{\rho}\leq\rho\leq\bar{\rho}.$

  (7)

  Moreover, $0<\underaccent{\bar}{a}\leq a\leq\bar{a}$ and $|gz|\leq\bar{g}\bar{z}$ for appropriate constants $\underaccent{\bar}{a},\bar{a}$, and $\bar{g}\bar{z}$.

• (A2)

  Sufficiently smooth solutions $(\rho,w)$ and $(\hat{\rho},\hat{w})$ to the system (1)–(5) exist for parameters $0\leq\varepsilon,\hat{\varepsilon}\leq\bar{\varepsilon}$ and $0<\underaccent{\bar}{\gamma}\leq\gamma,\hat{\gamma}\leq\bar{\gamma}$ with $\underaccent{\bar}{\gamma},\bar{\gamma}$ constant, which satisfy

  $\displaystyle\underaccent{\bar}{\rho}\leq\rho,\hat{\rho}\leq\bar{\rho}\qquad\text{and}\qquad-\bar{w}\leq w,\hat{w}\leq\bar{w}.$

  (8)

Conditions (A1) and (A2) imply uniform bounds for $P$ and its derivatives and ensure that the flow is subsonic; we refer to [6] for details. Under these conditions, we will show the following stability result: Let $(\rho,w)$ and $(\hat{\rho},\hat{w})$ be sufficiently regular solutions of (1)–(3) with parameters $\varepsilon,\gamma$ and $\hat{\varepsilon},\hat{\gamma}$, and boundary values $\hat{h}_{\partial}$ and $h_{\partial}$ as described in (6). Then

	$\displaystyle\|\rho(\tau)-\hat{\rho}(\tau)\|_{L^{2}}^{2}$	$\displaystyle+\varepsilon^{2}\|w(\tau)-\hat{w}(\tau)\|^{2}_{L^{2}}+\int_{0}^{\tau}\|w(s)-\hat{w}(s)\|^{3}_{L^{3}}ds$
		$\displaystyle\leq\hat{C}e^{\hat{c}\tau}\Big{(}\|\rho(0)-\hat{\rho}(0)\|_{L^{2}}^{2}+\|\varepsilon w(0)-\varepsilon\hat{w}(0)\|_{L^{2}}^{2}$
		$\displaystyle\qquad\qquad\qquad+|\gamma-\hat{\gamma}|^{3/2}+|\varepsilon^{2}-\hat{\varepsilon}^{2}|+\int_{0}^{\tau}|h_{\partial}(s)-\hat{h}_{\partial}(s)|\,ds\Big{)},$

with constants $\hat{c},\hat{C}$ depending only on the bounds in assumptions (A1)–(A2) and on bounds for time derivatives of $\hat{\rho}$ and $\hat{w}$. Note that the energy (4) and dissipation functional (5), and as a consequence also the co-state variables (3) depend explicitly on the parameters $\varepsilon$ and $\gamma$, so their definition changes for the perturbed problem. A precise statement of our results and of the regularity assumption on the solutions is given in Section 4, where we also discuss various generalizations including perturbations in other model parameters and the choice of different boundary conditions. Let us mention two immediate consequences of the above stability estimate, namely

• •

  uniqueness of regular subsonic bounded state solutions for specified initial and boundary values, and their stability with respect to perturbations in these problem data as well as the model parameters;

• •

  convergence of solutions to those of the parabolic limit problem which results by formally setting $\varepsilon=0$ in equations (1)–(6).

Existence and uniqueness of solutions to the parabolic limit problem can be proven rigorously by variational arguments; see [1, 23] or [25]. This parabolic problem also serves as the basis for simulation codes utilized in the gas network community [2, 22], and the stability estimate above allows us to obtain a quantitative justification for the use of this model. Due to the energy-based modelling, the results can be generalized almost verbatim to gas networks by utilizing appropriate coupling conditions; details will be discussed in Section 5.

The proof of our main result is based on the observation that (1)–(3) can be written as an abstract port-Hamiltonian system [27]

$\displaystyle\mathcal{C}\partial_{\tau}{\boldsymbol{u}}+(\mathcal{J}+\mathcal{R}({\boldsymbol{u}})){\boldsymbol{z}}({\boldsymbol{u}})$

$\displaystyle=\mathcal{B}_{\partial}{\boldsymbol{z}}({\boldsymbol{u}}),$

(9)

with ${\boldsymbol{u}}=(\rho,w)$ and ${\boldsymbol{z}}({\boldsymbol{u}})=(h,m)$ denoting the state and co-state variables, which are linked via the energy functional $\mathcal{H}({\boldsymbol{u}})=\mathcal{H}(\rho,w)$; moreover, $\mathcal{C}$, $\mathcal{J}$, $\mathcal{R}({\boldsymbol{u}})$, and $\mathcal{B}_{\partial}$, are appropriate operators, the last one incorporating the boundary conditions; see Section 2.

Any smooth function $\widehat{\boldsymbol{u}}=(\hat{\rho},\hat{w})$, e.g., a solution of (1)–(3) with perturbed model parameters, may be considered as a solution of the system

$\displaystyle\mathcal{C}\partial_{\tau}\widehat{\boldsymbol{u}}+(\mathcal{J}+\mathcal{R}(\widehat{\boldsymbol{u}})){\boldsymbol{z}}(\widehat{\boldsymbol{u}})$

$\displaystyle=\mathcal{B}_{\partial}{\boldsymbol{z}}(\widehat{\boldsymbol{u}})+\widehat{\boldsymbol{r}},$

(10)

with the same operators $\mathcal{C}$, $\mathcal{J}$, $\mathcal{R}(\cdot)$, $\mathcal{B}_{\partial}$, and the same energy functional $\mathcal{H}(\cdot)$ and state to co-state mapping ${\boldsymbol{z}}(\cdot)$, up to some perturbation $\widehat{\boldsymbol{r}}$. The coincidence of the underlying Hamiltonian structure will allow us to estimate the difference between ${\boldsymbol{u}}$ and $\widehat{\boldsymbol{u}}$ in terms of the perturbations $\widehat{\boldsymbol{r}}$ by means of relative energy estimates. For the stability analysis on the abstract level, we require some general conditions that will be verified in detail for the gas transport problem under investigation using the assumptions (A1)–(A2).

The use of an abstract problem structure greatly simplifies the analysis and allows to generalize our results in various directions. We briefly discuss the incorporation of perturbations in other model parameters and the extension to gas networks.

Relative entropy or energy techniques have been used intensively for the existence, stability, and discretization error analysis of time dependent partial differential equations. We refer to [16] for a recent summary of corresponding results for parabolic evolution problems. In the present paper, we are interested in hyperbolic problems, where the use of relative entropy arguments goes back to the seminal works of DiPerna [7] and Dafermos [5]; also see [6] for an introduction to the field. Typical aspects that are addressed are: convergence to steady states, stable dependence of solutions on initial data and parameters, and asymptotic limits. Examples for the latter include the low Mach limit of Euler and Navier-Stokes equations, which are investigated in [10] for instance. Long time convergence of solutions to damped Euler equations to Barenblatt solutions were investigated by Huang and coworkers in a series of papers [11].

One particular aspect that we want to address in the present study are parabolic limits of quasilinear hyperbolic equations; see [21, 15, 17, 18] for some exemplary results in this direction. The latter reference as well as [4, 12] strongly rely on the underlying dissipative Hamiltonian structure of many equations in fluid mechanics, which will also play a major role in our analysis below. The previous papers, however, use formulations in conservative variables, which allows to deal with one solution being a weak solution in the classical sense of hyperbolic conservation laws [6]. Parabolic limits of hyperbolic $2\times 2$ systems have also been studied in [20, 29] using compensated compactness arguments [13, 14]. This has the advantage that no additional regularity of solutions to the parabolic limit problem is required, but no quantitative information about the speed of convergence is obtained. Spectral estimates for the linear part of the problem are used in [8] to derive convergence in the parabolic limit. Let us note that most of the above works consider only linear friction laws and unbounded domains or periodic boundary conditions.

In contrast to work mentioned previously, our study is based on a formulation in primitive variables, which requires both solutions to have some minimal smoothness. On the other hand, this formulation allows us to incorporate boundary conditions more naturally and to extend our results to networks in a straight-forward manner using appropriate coupling conditions at network junctions [24, 9] that guarantee energy conservation or dissipation at network junctions. Similar formulations for compressible flow were also considered in the context of port-Hamiltonian systems; see [27, 26] for the models and [3, 19] for corresponding discretization strategies. Other systems, that fit into the general framework that we develop in this paper include the Euler-Korteweg system, the system of quantum hydrodynamics and the Euler-Poisson equations. An overview about corresponding results can be found in [12].

The remainder of the manuscript is organized as follows: Section 2 is concerned with the perturbation analysis for the abstract system (9). In Section 3, we then briefly review the derivation of the model equations (1)–(3) starting from the usual balance equations of gas dynamics, and we show that they fit into the abstract form (9). In Section 4, we verify the assumptions required for our abstract analysis for the gas transport problem under investigation, which allows us to state and prove our main results. Their generalization to gas networks is discussed in Section 5.

Let $\mathbb{W}\subset\mathbb{V}$ be real Hilbert spaces, with $\mathbb{V}^{\prime}\subset\mathbb{W}^{\prime}$ denoting the corresponding dual spaces. We use $\langle\cdot,\cdot\rangle$ to denote the duality product on $\mathbb{V}^{\prime}\times\mathbb{V}$ and $\mathbb{W}^{\prime}\times\mathbb{W}$. Furthermore, let $\mathcal{H}:\mathbb{D}\subset\mathbb{V}\to\mathbb{R}$ be a smooth and strictly convex energy functional, with $\mathbb{D}\subset\mathbb{V}$ denoting some appropriate, e.g., convex and closed, subset. We consider abstract evolution problems of the form (11)–(12), with operators satisfying the following conditions.

From conditions (13)–(14), one can immediately see that

$\displaystyle\langle{\boldsymbol{u}},{\boldsymbol{v}}\rangle_{\mathcal{C}}:=\langle\mathcal{C}{\boldsymbol{v}},{\boldsymbol{u}}\rangle$

(18)

defines a scalar product on $\mathbb{V}$ and the associated norm $\|{\boldsymbol{v}}\|_{\mathcal{C}}^{2}=\langle{\boldsymbol{v}},{\boldsymbol{v}}\rangle_{\mathcal{C}}=\langle\mathcal{C}{\boldsymbol{v}},{\boldsymbol{v}}\rangle$ is equivalent to the standard norm $\|\cdot\|_{\mathbb{V}}$. The expression ${\boldsymbol{z}}({\boldsymbol{u}})=\mathcal{C}^{-1}\mathcal{H}^{\prime}({\boldsymbol{u}})=\operatorname{grad}_{\mathcal{C}}\mathcal{H}({\boldsymbol{u}})$ then denotes the gradient of the functional $\mathcal{H}$ at ${\boldsymbol{u}}$ with respect to this scalar product. We further introduce the symbol

$\displaystyle\mathcal{G}({\boldsymbol{u}})=\mathcal{C}^{-1}\mathcal{H}^{\prime\prime}({\boldsymbol{u}})$

(19)

for the Hessian operator $\mathcal{G}({\boldsymbol{u}}):\mathbb{V}\to\mathbb{V}^{\prime}$, and note that

$\displaystyle\langle\mathcal{G}({\boldsymbol{u}}){\boldsymbol{v}},{\boldsymbol{w}}\rangle_{\mathcal{C}}=\langle\mathcal{H}^{\prime\prime}({\boldsymbol{u}}){\boldsymbol{v}},{\boldsymbol{w}}\rangle=\langle\mathcal{H}^{\prime\prime}({\boldsymbol{u}}){\boldsymbol{w}},{\boldsymbol{v}}\rangle=\langle\mathcal{G}({\boldsymbol{u}}){\boldsymbol{w}},{\boldsymbol{v}}\rangle_{\mathcal{C}},$

(20)

i.e., the Hessian is symmetric with respect to the scalar product induced by $\mathcal{C}$.

As a direct consequence of the underlying port-Hamiltonian structure and our assumptions on the operators, we obtain the following power balance relation.

We now study the stability of solutions to (11)–(12) with respect to perturbations. To this end, let $\widehat{\boldsymbol{u}}$ denote a classical solution of

	$\displaystyle\mathcal{C}\partial_{\tau}\widehat{\boldsymbol{u}}+[\mathcal{J}+\mathcal{R}(\widehat{\boldsymbol{u}})]{\boldsymbol{z}}(\widehat{\boldsymbol{u}})$	$\displaystyle=\mathcal{B}_{\partial}{\boldsymbol{z}}({\boldsymbol{u}})+\widehat{\boldsymbol{r}},$		(22)
	$\displaystyle z(\widehat{\boldsymbol{u}})$	$\displaystyle=\mathcal{C}^{-1}\mathcal{H}^{\prime}(\widehat{\boldsymbol{u}}),$		(23)

with appropriate perturbation described by the residual functional $\widehat{\boldsymbol{r}}\in C^{0}([0,T];\mathbb{W}^{\prime})$. As a measure for the difference of ${\boldsymbol{u}}$ and $\widehat{\boldsymbol{u}}$, we utilize the relative energy [6], defined by

$\displaystyle\mathcal{H}({\boldsymbol{u}}|\widehat{\boldsymbol{u}}):=\mathcal{H}({\boldsymbol{u}})-\mathcal{H}(\widehat{\boldsymbol{u}})-\langle\mathcal{H}^{\prime}(\widehat{\boldsymbol{u}}),{\boldsymbol{u}}-\widehat{\boldsymbol{u}}\rangle.$

(24)

Using the particular problem structure and some elementary computations, we can prove the following basic identity for the temporal change of the relative entropy.

We now derive quantitative estimates for the difference of the solutions ${\boldsymbol{u}}$ and $\widehat{\boldsymbol{u}}$ to (11)–(12) and (22)–(23) with respect to perturbations in the right hand side and the initial and boundary values. To do so, we make some abstract assumptions that will later be verified for the gas transport problem under consideration.

Together with the relative energy identity stated in the previous lemma, these abstract conditions immediately lead to the following stability estimate.

We may multiply (1)–(2) by appropriate test functions $q$, $r$, integrate over the spatial domain, and use integration-by-parts in the second equation, to obtain

	$\displaystyle(a\partial_{\tau}\rho,q)+(\partial_{x}m,q)$	$\displaystyle=0,$		(25)
	$\displaystyle(\varepsilon^{2}\partial_{\tau}w,r)-(h,\partial_{x}r)+(\gamma\frac{|w|}{a\rho}m,r)$	$\displaystyle=-hr\big{|}_{0}^{\ell},$		(26)

where $(a,b):=\int_{0}^{\ell}a(x)b(x)dx$ is used to abbreviate the $L^{2}$-scalar product. Note that boundary conditions (6) for $h$ could be incorporated naturally in the last term of (26). The two variational identities (25)–(26) hold for all time $t>0$ of relevance and for all smooth test functions $q$, $r$, independent of time, and they are satisfied, in particular, by all smooth solutions of (1)–(3). In compact notation, the system (25)–(26) can be stated as

$\displaystyle\langle\mathcal{C}\partial_{\tau}{\boldsymbol{u}},{\boldsymbol{w}}\rangle+\langle\mathcal{J}{\boldsymbol{z}}({\boldsymbol{u}}),{\boldsymbol{w}}\rangle+\langle\mathcal{R}({\boldsymbol{u}}){\boldsymbol{z}}({\boldsymbol{u}}),{\boldsymbol{w}}\rangle$

$\displaystyle=\langle\mathcal{B}_{\partial}{\boldsymbol{z}},{\boldsymbol{w}}\rangle,$

(27)

with state variable ${\boldsymbol{u}}=(\rho,w)$, co-state variable ${\boldsymbol{z}}({\boldsymbol{u}})=(h,m)$ defined by (3), time independent test function ${\boldsymbol{v}}=(q,r)$, and operators $\mathcal{C}$, $\mathcal{J}$, $\mathcal{R}({\boldsymbol{u}})$, and $\mathcal{B}_{\partial}$ given by

	$\displaystyle\langle\mathcal{C}{\boldsymbol{u}},{\boldsymbol{v}}\rangle$	$\displaystyle=(a\rho,q)+(\varepsilon^{2}w,r),$	$\displaystyle\langle\mathcal{R}({\boldsymbol{u}}){\boldsymbol{z}},{\boldsymbol{v}}\rangle$	$\displaystyle=(\gamma\frac{|w|}{a\rho}m,r)$
	$\displaystyle\langle\mathcal{J}{\boldsymbol{z}},{\boldsymbol{v}}\rangle$	$\displaystyle=(\partial_{x}m,q)-(h,\partial_{x}r),$	$\displaystyle\langle\mathcal{B}{\boldsymbol{z}}_{\boldsymbol{\partial}},{\boldsymbol{v}}\rangle$	$\displaystyle=-hr\big{|}_{0}^{\ell}.$

From these variational characterizations, it is not difficult to see that $\mathcal{J}$ is anti-symmetric and that $\mathcal{C}$ and $\mathcal{R}({\boldsymbol{u}})$ are symmetric and at least positive semi-definite, if the parameters $a$, $\varepsilon^{2}$, $\gamma$ and the density $\rho$ are positive. The operator $\mathcal{B}_{\partial}$ associated with the boundary terms does not have a particular property, except being supported only at the boundary.

A quick inspection of the above derivations shows that the operators $\mathcal{C}$, $\mathcal{R}({\boldsymbol{u}})$, and $\mathcal{J}$, and $\mathcal{B}_{\partial}$, can be formally identified with

$\displaystyle\mathcal{C}=\begin{pmatrix}a&0\\ 0&\varepsilon^{2}\end{pmatrix},\qquad\mathcal{R}({\boldsymbol{u}})=\begin{pmatrix}0&0\\ 0&\frac{\gamma|w|}{a\rho}\end{pmatrix},\qquad\text{and}\qquad\mathcal{J}-\mathcal{B}_{\partial}=\begin{pmatrix}0&\partial_{x}\\ \partial_{x}&0\end{pmatrix}.$

The latter follows rigorously by reversing the order of arguments in the derivations of the weak formulation. The energy of the system is here given by

$\displaystyle\mathcal{H}({\boldsymbol{u}})=\int_{0}^{\ell}a(\varepsilon^{2}\rho\frac{w^{2}}{2}+P(\rho)+\rho gz)dx,$

and by elementary computations, we obtain the formulas

$\displaystyle{\boldsymbol{z}}({\boldsymbol{u}})=\begin{pmatrix}\varepsilon^{2}\frac{w^{2}}{2}+P^{\prime}(\rho)\\ a\rho w\end{pmatrix}\qquad\text{and}\qquad\mathcal{G}({\boldsymbol{u}})=\begin{pmatrix}P^{\prime\prime}(\rho)&\varepsilon^{2}w\\ aw&a\rho\end{pmatrix}$

for the gradient ${\boldsymbol{z}}({\boldsymbol{u}})=\mathcal{C}^{-1}\mathcal{H}^{\prime}({\boldsymbol{u}})$ and Hessian $\mathcal{G}({\boldsymbol{u}})=\mathcal{C}^{-1}\mathcal{H}^{\prime\prime}({\boldsymbol{u}})$ of the energy functional.

As a next step, we briefly discuss the choice of suitable function spaces for the gas transport problem under consideration. We define

$\displaystyle\mathbb{V}:=L^{2}(0,\ell)\times L^{2}(0,\ell)\qquad\text{and}\qquad\mathbb{W}:=H^{1}(0,\ell)\times H^{1}(0,\ell),$

(28)

where $H^{1}(0,\ell)$ denotes the standard Sobolev space on the interval $(0,\ell)$. We further introduce the set of admissible states

$\displaystyle\mathbb{D}=\{(\rho,w)\in\mathbb{V}:(A1)-(A2)\text{ are satisfied}\}.$

Let us recall that classical solutions of (1)–(3) satisfy ${\boldsymbol{u}}=(\rho,w)\in C^{1}([0,T];\mathbb{V})\cap C^{0}([0,T];\mathbb{D})$ and ${\boldsymbol{z}}({\boldsymbol{u}})=(h,m)\in C^{0}([0,T];\mathbb{W})$. Then $\mathcal{C}$ and $\mathcal{R}({\boldsymbol{u}})$ with ${\boldsymbol{u}}\in\mathbb{D}$ can be understood as self-adjoint and positive semi-definite bounded linear operators mapping from $\mathbb{V}$ or $\mathbb{W}$ to the dual spaces $\mathbb{V}^{\prime}$ or $\mathbb{W}^{\prime}$. For $\varepsilon$, $a$ uniformly positive and bounded, $\mathcal{C}$ induces a norm

$\displaystyle\|{\boldsymbol{u}}\|_{\mathcal{C}}^{2}=\|\sqrt{a}\rho\|^{2}_{L^{2}(0,\ell)}+\|\varepsilon w\|^{2}_{L^{2}(0,\ell)},$

(29)

which is equivalent to the standard norm on $\mathbb{V}=L^{2}(0,\ell)\times L^{2}(0,\ell)$. Adopting the previous notation, we write $\langle\cdot,\cdot\rangle$ for the duality products on $\mathbb{V}^{\prime}\times\mathbb{V}$ and $\mathbb{W}^{\prime}\times\mathbb{W}$, respectively, and use $(a,b)=\int_{0}^{\ell}a(x)b(x)\,dx$ to denote the scalar product of $L^{2}(0,\ell)$. From the variational definition of the operator $\mathcal{J}$, one can see that

$\displaystyle\langle\mathcal{J}{\boldsymbol{w}},\widetilde{\boldsymbol{w}}\rangle:=$

$\displaystyle(\partial_{x}m,\tilde{h})-(h,\partial_{x}\tilde{m})=-\langle\mathcal{J}\widetilde{\boldsymbol{w}},{\boldsymbol{w}}\rangle,$

(30)

for all ${\boldsymbol{w}}=(h,m)$, $\widetilde{\boldsymbol{w}}=(\tilde{h},\tilde{m})\in\mathbb{W}$, i.e., $\mathcal{J}$ is skew-symmetric on $\mathbb{W}$. The formula $\langle\mathcal{B}_{\partial}{\boldsymbol{z}}({\boldsymbol{u}}),{\boldsymbol{v}}\rangle=-hr\big{|}_{0}^{\ell}$ with ${\boldsymbol{v}}=(q,r)$ finally shows that the boundary operator $\mathcal{B}_{\partial}$ acts on the co-state variables. Note that the required boundary values are well-defined for functions ${\boldsymbol{z}}({\boldsymbol{u}})=(h,m)$ and ${\boldsymbol{v}}=(q,r)\in\mathbb{W}=H^{1}(0,\ell)\times H^{1}(0,\ell)$.

In summary, we thus have shown that (1)–(3) can be interpreted as an abstract port-Hamiltonian system (11)–(12), and under conditions (A1)–(A2) also Assumption 1 is valid.

Let $(\hat{\rho},\hat{w})$ denote a solution to the perturbed equations

	$\displaystyle a\partial_{\tau}\hat{\rho}+\partial_{x}(a\hat{\rho}\hat{w})$	$\displaystyle=0,$		(31)
	$\displaystyle\hat{\varepsilon}^{2}\partial_{\tau}\hat{w}+\partial_{x}(P^{\prime}(\hat{\rho})+\hat{\varepsilon}^{2}\frac{\hat{w}^{2}}{2}+gz)+\hat{\gamma}\frac{|\hat{w}|}{a\hat{\rho}}\hat{m}$	$\displaystyle=0,$		(32)

which are again assumed to hold for all $0<x<\ell$ and time $t>0$. The terms carrying spatial derivatives are the corresponding co-state variables

$\displaystyle\hat{h}(\hat{\rho},\hat{v})=\hat{\varepsilon}^{2}\frac{\hat{w}^{2}}{2}+P^{\prime}(\hat{\rho})+gz\qquad\text{and}\qquad\hat{m}(\hat{\rho},\hat{v})=a\hat{\rho}\hat{w},$

(33)

which are again directly related to the derivatives of the corresponding energy functional $\hat{\mathcal{H}}(\hat{\rho},\hat{w})=\int_{0}^{\ell}a(\hat{\varepsilon}^{2}\frac{\hat{w}^{2}}{2\hat{\rho}}+P(\hat{\rho})+\hat{\rho}gz)\,dx$. By the some elementary manipulations and the same arguments as employed above, this system can again be written in the abstract form

	$\displaystyle\mathcal{C}\partial_{\tau}\widehat{\boldsymbol{u}}+(\mathcal{J}+\mathcal{R}(\widehat{\boldsymbol{u}})){\boldsymbol{z}}(\widehat{\boldsymbol{u}})$	$\displaystyle=\mathcal{B}_{\partial}{\boldsymbol{z}}(\widehat{\boldsymbol{u}})+\widehat{\boldsymbol{r}}$		(34)
	$\displaystyle{\boldsymbol{z}}(\widehat{\boldsymbol{u}})$	$\displaystyle=\mathcal{C}^{-1}\mathcal{H}^{\prime}(\widehat{\boldsymbol{u}}),$		(35)

with the functionals $\mathcal{H}(\cdot)$ and ${\boldsymbol{z}}(\cdot)$, and the operators $\mathcal{C}$, $\mathcal{J}$, and $\mathcal{R}(\cdot)$ denoting those for the unperturbed problem with parameters $\varepsilon$ and $\gamma$, and residual $\widehat{\boldsymbol{r}}=(\widehat{\boldsymbol{r}}_{1},\widehat{\boldsymbol{r}}_{2})$ given by

$\displaystyle\widehat{\boldsymbol{r}}_{1}=0\qquad\text{and}\qquad\widehat{\boldsymbol{r}}_{2}=(\varepsilon^{2}-\hat{\varepsilon}^{2})(\partial_{\tau}\hat{w}+\tfrac{1}{2}\partial_{x}|\hat{w}|^{2})+(\gamma-\hat{\gamma})|\hat{w}|\hat{w}.$

(36)

In the following section, we verify the abstract assumptions required for our stability analysis, without explicitly taking into account the special form of $\widehat{\boldsymbol{r}}_{1}$ and $\widehat{\boldsymbol{r}}_{2}$.

We always assume in the following that assumptions (A1)–(A2) are valid. Constants arising in the estimates may depend on the bounds of these conditions. Since we later consider the case $\varepsilon\to 0$, we will make explicit the dependence of the constants on this scaling parameter.

From Taylor’s formula, we know that

	$\displaystyle f({\boldsymbol{u}}|\hat{\boldsymbol{u}})$	$\displaystyle=f({\boldsymbol{u}})-f(\widehat{\boldsymbol{u}})-f^{\prime}(\widehat{\boldsymbol{u}})({\boldsymbol{u}}-\widehat{\boldsymbol{u}})$
		$\displaystyle=\frac{1}{2}\int_{0}^{1}\langle(1-s)f^{\prime\prime}(\widehat{\boldsymbol{u}}+s({\boldsymbol{u}}-\widehat{\boldsymbol{u}}))({\boldsymbol{u}}-\widehat{\boldsymbol{u}}),{\boldsymbol{u}}-\widehat{\boldsymbol{u}})\rangle ds.$

Now let $f({\boldsymbol{u}})=a(\varepsilon^{2}\rho\frac{|w|^{2}}{2}+P(\rho))$ denote the integrand of the energy functional defined above, and note that its Hessian is given by

$\displaystyle f^{\prime\prime}(\rho,w)=a\begin{pmatrix}P^{\prime\prime}(\rho)&\varepsilon^{2}w\\ \varepsilon^{2}w&\varepsilon^{2}\rho\end{pmatrix}.$

By multiplying with ${\boldsymbol{v}}=(x,y)$ from the left and right, one can see that

	$\displaystyle\frac{1}{a}\langle f^{\prime\prime}(\rho,w){\boldsymbol{v}},{\boldsymbol{v}}\rangle$	$\displaystyle=P^{\prime\prime}(\rho)x^{2}+2\varepsilon^{2}wxy+\varepsilon^{2}\rho y^{2}$
		$\displaystyle=P^{\prime\prime}(\rho)x^{2}+2\varepsilon wx(\varepsilon y)+\rho(\varepsilon y)^{2}.$

The second term in the second line can be estimated by

$\displaystyle|2\varepsilon wx(\varepsilon y)|$

$\displaystyle\leq 2\frac{\varepsilon^{2}w^{2}}{\rho}x^{2}+\frac{1}{2}\rho(\varepsilon y)^{2}.$

From condition (7), one can see that $2\frac{\varepsilon^{2}w^{2}}{\rho}\leq P^{\prime\prime}(\rho)/2$, which leads to the lower bound

	$\displaystyle\langle f^{\prime\prime}(\rho,w){\boldsymbol{v}},{\boldsymbol{v}}\rangle\geq\frac{a}{2}\left(P^{\prime\prime}(\rho)x^{2}+\rho|\varepsilon y|^{2}\right),$
and the corresponding upper bound
	$\displaystyle\langle f^{\prime\prime}(\rho,w){\boldsymbol{v}},{\boldsymbol{v}}\rangle\leq\frac{3a}{2}\left(P^{\prime\prime}(\rho)x^{2}+\rho|\varepsilon y|^{2}\right).$

Using the uniform bounds for $\rho$, $P^{\prime\prime}(\rho)$, and $a$, we arrive at the following result.