Numerical spectral analysis of standing waves in quantum hydrodynamics with viscosity

1.   Introduction

In this paper, we consider two quantum hydrodynamics (QHD) systems for which we investigate the spectrum of the linearized operator about standing wave profiles. The system with a linear viscosity term reads

where $\rho \geq 0$ is the density, $m = \rho u$ is the momentum, where $u$ is the fluid velocity, $p(\rho) = \rho^{\gamma}$ with $\gamma \geq 1$ is the pressure, $t \geq 0$ and $x \in \mathbb{R}$. Moreover, $\mu, k > 0$ are viscosity and dispersive coefficients, respectively. The function $(\sqrt{\rho})_{xx}/\sqrt{\rho}$ is known as the Bohm potential, providing the model with a nonlinear third order dispersive term. These systems are used for instance in the theory of superfluidity ^[15] or to model semiconductor devices. In the case of nonlinear viscosity (see (3.1)) the term $\mu m_{xx}$ is replaced with $\mu \rho (m_x/\rho)_x$. The latter describes the interaction between the superfluid and a normal fluid, or it can be interpreted as describing the interaction of the fluid with a background. Systems (1.1) and (3.1) can be viewed as mean-field limits of the de Broglie-Bohm pilot wave theory ^[4,5,6] with added viscosity terms of superfluid type.

Standing waves for (1.1) are solutions of the form

such that

The first studies of models with dispersive terms were ^[12,22]; see also ^[11,13,20]. Moreover, existence of weak solutions for QHD systems has been considered in ^[1,2,3]. The spectral stability of traveling wave profiles for the $p$-system has been discussed in ^[14]. Concerning viscous QHD models, existence of traveling wave profiles for the system (1.1) was shown in ^[16] provided that the viscosity is sufficiently strong. The result was improved in ^[24] where global existence of dispersive shocks was proved without any restriction on the viscosity and dispersive coefficients and for arbitrary shock amplitude. A similar result about a related QHD system with nonlinear viscosity (see (3.1)) was obtained in ^[19]. Moreover, stability of traveling waves was established numerically in ^[17,18] in the case of non-monotone shocks, and analytically in ^[8,9] for small-amplitude profiles. These studies were complemented in ^[21] with a proof of local well-posedness and nonlinear decay of perturbations of subsonic constant states for (1.1).

The aim of this paper is to study the spectrum of the linearization around standing waves for (1.1) and (3.2). Existence theory for such solutions was established in ^[24]. There, it was shown that the standing waves admited by the system (1.1) are pulses, i.e., $P^+ = P^-$ and $J^+ = J^-$. Here below we recall the existence result [24, Theorem 4.1]. Let us denote by $c_s(\rho) = \sqrt{\gamma \rho^{\gamma-1}}$ the speed of sound and by $U^+ = J^+/P^+$ the velocity at the end state.

Theorem 1.1 (existence of standing waves). There exists a non-constant standing wave solution of (1.1) with $\lim_{x \rightarrow \pm \infty} [P(x), J(x)] = [P^+, J^+],$ if and only if $0 < |U^+| < c_s(P^+)$.

This result shows, in particular, that there exist no standing waves with supersonic or sonic end states.

The spectrum of the linearization around standing wave profiles for systems (1.1) and (3.2) consists of two parts: the essential spectrum and the point spectrum. It can be proved that the essential spectrum of such profiles is always stable (see Theorem 2.2 below). In the present paper, we carry out a numerical study of the point spectrum by using the Evans function method. We show that there exists an unstable real simple eigenvalue for both QHD systems with linear and nonlinear viscosity. Thereby, we provide numerical evidence for spectral instability of standing waves. This result is in contrast with the analytical and numerical results about spectral stability of traveling waves in ^[8,9,17,18]. A possible cause of the instability is the presence of a supersonic region along the standing wave profiles (see Section 4 and Conjecture 4.1 below).

The paper is organized as follows. In Section 2 we consider the QHD system with linear viscosity (1.1). We introduce the equation solved by standing wave profiles and discuss their numerical computation. Then, we derive the linearization around such profiles, describe its essential spectrum, and recast the system in integrated variables. Section 2.6 contains the numerical result about point spectrum instability. In Section 3 we discuss the QHD system with nonlinear viscosity (see Eq (3.2) below) following the steps for the linear viscosity case. Finally, we conclude the paper with a discussion of the numerical results.

Remark 1.2. The linearizations in Sections 2.3 and 3.2 can be obtained by setting $s = 0$ in the corresponding linearizations in ^[16,18] which hold also for $s = 0$. We provide them here for completeness.

2.   Quantum hydrodynamics with linear viscosity

2.1.   Profile equation

In this section we shall derive the equation satisfied by a standing wave profile for system (1.1). The Bohm potential can be rewritten in conservative form

Substituting the standing wave profiles $P$ and $J$ we get the following system of ODEs

From (2.1) it follows that the momentum is conserved along the profile

where $A = J^+$. Since $J$ is constant we have $J'' = 0$, hence the term coming from the viscosity in Eq (2.2) vanishes. Susbstituting Eq (2.3) in Eq (2.2) and proceeding as in [16, Section 2], we conclude that the density profile $P$ solves the ODE

where $f(P) = P^{\gamma} - B + \frac{A^2}{P},$ with $B = \left (\frac{m^2}{\rho} + \rho^{\gamma} \right)^+.$

2.2.   Numerical computation of standing wave profiles

We shall compute numerically standing waves whose existence is guaranteed by Theorem 1.1. They correspond to homoclinic loops for the profile ODE (2.4). Since the steady-state for this dynamical system to which the profile converges is a saddle, the method to compute the profile from ^[17] cannot be used because it relies on one of the steady-states being stable. Instead, the profile will be computed by solving a nonlinear boundary value problem.

To compute standing wave profiles we choose $P^+ > 0$ and $U^+$ that satisfy the condition of Theorem 1.1. Then, we choose a sufficiently large domain size $L_1 > 0$ and integrate numerically (2.4) on the domain $[-L_1, L_1]$ subject to the boundary conditions $P(-L_1) = P(L_1) = P^+,$ with the boundary value solver bvp4c in Matlab. As an initial guess we use the function

for $-L_1 \leq x \leq L_1$, where $c > 0$ is a parameter to be chosen so that the solver converges. From the qualitative analysis in [24, Section 4] we know that the profile has a single minimum, hence the initial guess has an appropriate shape. Using this initial guess the solver will not converge to the trivial constant solution.

Because of the fact that the momentum is conserved along the profile we have $u(x) = J^+/P(x)$. Hence, the flow is supersonic in the region where the following inequality holds

where the right-hand side is the density corresponding to sonic states. For our computations we use the parameter values

We have $|U^+| = 0.8 < 1.22 \approx c_s(P^+)$, therefore the condition of Theorem 1.1 is satisfied. Moreover, we set $L_1 = 20$ and $c = 0.8$. For parameters (2.7) the right hand side of (2.6) is $0.71$ and the flow is supersonic in the interval $|x| < 1.93$ (see Figure 1). This type of reduced density pulse is known as dark soliton (see, for instance ^[7]).

2.3.   Linearization

Let us denote by $(\rho, m)$ the deviation from $(P, J)$. We obtain the full linearized operator around the profile

where

The associated eigenvalue problem reads

2.4.   Essential spectrum

Since the profile is a pulse, the essential spectrum of (2.8) is obtained by considering the asymptotic operator at the end state

where

The associated eigenvalue problem is

Let us rewrite it as a first order system

where $V = [\rho, m, u_1, u_2]^{\top}$, with $\rho' = u_1$, $u_1' = u_2$, and the limit matrix $M^+$ is given by

The characteristic equation of $M^+$ is $\det(\nu Id - M^+) = 0$. It reads

Setting $\nu = i \xi$, $\xi \in \mathbb{R}$ in (2.11) and dividing by $2/k^2$ we obtain the dispersion relation

Since the profile $(P, J)$ is a pulse, the essential spectrum is given by the union of the two curves $\Sigma_j = \{ \lambda = \lambda_j(\xi) \in \mathbb{C}: \xi \in \mathbb{R} \}$, $j = 1, 2$, where $\lambda_j(\xi)$ are determined by the roots of (2.12). Let us now define stability of essential spectrum.

Definition 2.1. The essential spectrum of $\mathcal{L}$ is stable, if it is contained in the closed left half-plane, that is

The essential spectum of $\mathcal{L}$ is unstable, if it intersects the unstable half-plane:

Lemma 3 in ^[16] and Lemma 3.2 in ^[21] imply the following description of the stability of the essential spectrum.

Theorem 2.2 (Stability of essential spectrum).

(ⅰ) If $|U^+| \leq c_s(P^+)$, then the essential spectrum of $\mathcal{L}$ is stable. Moreover, it holds that $\mathrm{Re}\, \lambda_{1, 2} < 0$, provided $\xi \neq 0$.

(ⅱ) If $|U^+| > c_s(P^+)$, then the essential spectum of $\mathcal{L}$ is unstable.

That is, the stability of the essential spectrum is related to the end state $(P^+, J^+)$ being subsonic or sonic. Moreover, the existence condition in Theorem 1.1 is also related to stability of the essential spectrum. Indeeed, a standing wave profile exists if and only if the end state is subsonic with nonzero velocity. Hence, the essential spectrum of standing waves is always stable. Furthermore, Lemma 3 in ^[16] implies that we have consistent splitting, that is the limit matrix $M^+$ has 2 eigenvalues with positive real parts and 2 eigenvalues with negative real parts provided $\lambda$ is to the right of the curves $\Sigma_j$.

2.5.   System in integrated variables

Following [16, Section 4.2.1], we shall recast the eigenvalue problem (2.9) in terms of integrated variables

This transformation removes the zero eigenvalue without further modifications of the spectrum. Expressing $\rho$ and $m$ in terms of $\hat{\rho}$ and $\hat{m}$ and integrating from $-\infty$ to $x$ we get

where

Let us rewrite the system (2.13)-(2.14) as $\hat{V}' = \hat{M}(x, \lambda) \hat{V}$, where $\hat{V} = [ \hat{\rho}, \hat{m}, \hat{u}_1, \hat{u}_2]^{\top}$, $\hat{\rho}' = \hat{u}_1$, $\hat{u}_1' = \hat{u}_2$, and

The limiting matrix of $\hat{M}(x, \lambda)$ as $x \rightarrow \pm \infty$ is (2.10).

2.6.   Numerical evaluation of the Evans function

To define the Evans function, let us consider the equation $Y' = \hat{M}(x, \lambda)Y$ with $\hat{M}(x, \lambda)$ defined in (2.15). Since the profile $(P, J)$ is a pulse, the limits of $\hat{M}(x, \lambda)$ at $\pm \infty$ coincide and are equal to the matrix $M^+$ given by (2.10). We assume that $M^+$ is hyperbolic. This is always true for $\lambda$ to the right of the curves $\Sigma_j$. Denote by $\nu^-_1, \nu^-_2$ and by $\nu^+_1, \nu^+_2$ the eigenvalues of $M^+$ with positive and negative real parts, respectively, and indicate with $v_i^{\pm}$ the corresponding (normalized) eigenvectors. Let $Y^-_i$ be a solution of $Y' = \hat{M}(x, \lambda)Y$, satisfying $\exp(\nu^-_i x)Y^-(x)$ tends to $v^-_i$ as $x \rightarrow - \infty$ and $\exp(\nu^+_i x)Y^+(x)$ tends to $v^+_i$ as $x \rightarrow +\infty$. Then, the Evans function can be defined by

Consequently, $\lambda$ is in the point spectrum of $\mathcal{L}$ if and only if $E(\lambda) = 0$.

To numerically compute the Evans function, we use the compound matrix method (see ^[14,17,18]). This method is used in order to get a stable numerical procedure, in spite of the fact that the system $Y' = \hat{M}(x, \lambda)Y$ is numerically stiff. Specifically, the compound matrix $B(x, \lambda)$ is given by:

where $m_{jk}$ are the entries of $\hat{M}(x, \lambda)$ defined in (2.15). We integrate the equation $\phi' = (B(x, \lambda)-\mu^-)\phi$ numerically on a sufficiently large interval $[-L_1, 0]$, where $\mu^-$ is the unstable eigenvalue of $B$ at $-\infty$ with maximal (positive) real part. Denote the profile $(P(x), J(x))$ by $\zeta(x)$. Given a numerical approximation $(\zeta_k)_{k = 1}^N$ of $\zeta(x)$ at points $(x_k)_{k = 1}^N$ with $-L_1 = x_1 < x_2 < ... < x_N = L_1$, let $\tilde{\zeta}(x)$ be the piecewise linear interpolant of $(x_1, \zeta_1), ..., (x_N, \zeta_N)$. We obtain the matrix $B(x, \lambda)$ using $\tilde{\zeta}(x)$. Similarly we integrate the equation $\phi' = (B(x, \lambda)-\mu^+)\phi$ on $[0, L_1]$ backwards, where this time $\mu^+$ is the stable eigenvalue of $B$ at $+\infty$ with minimal (negative) real part. Then, the coefficients $\mu^{\pm}$ compensate for the growth/decay at infinity. Finally, the Evans function can be constructed by means of linear combination of the components of the two solutions $\phi^\pm = (\phi^\pm_1, \dots, \phi^\pm_6)$ as follows:

To compute the initial conditions we integrate the reduced Kato ODE

where $\mathcal{P}_{\pm}$ denotes the spectral projection of $\hat{B}_{\pm} = \lim_{x \rightarrow + \infty} \hat{B}(x, \lambda)$ corresponding to $\mu^{\pm}$. This choice of initial condition provides an analytic Evans function $E(\lambda)$. To numerically integrate (2.16) we use the algorithm from ^[25], that is $|r_{\pm}^1| = 1$ eigenvector as before (referring to maximal/minimal decay/growth rate of $\hat{B}^{\pm}$) and for $k > 0$,

Then, provided $E(\lambda)$ does not vanish on a closed contour $\Gamma$ which does not intersect the essential spectrum of $\mathcal{L}$, we can use the winding number

to count the number of zeros (taking into account multiplicity) inside the contour (see ^[23]).

For our calculations we use the set of parameters (2.7). Recall that the associated profile is depicted in Figure 1. We compute the Evans function on a semi-circular contour with radius 20, center at $\lambda = 0$ and vertical segment on the imaginary axis. Furthermore, we do not evaluate the Evans function at 0, but evaluate it up to $\pm i 10^{-6}$. Along the contour we integrate the Kato ODE using $5 \cdot 10^4$ points. Then, we compute $E(\lambda)$ with the solver ode45 in Matlab with relative tolerance $10^{-6}$ and we set $L_1 = 40$. Finally, we apply the symmetry of

The Evans function $E(\lambda)$ is plotted in Figure 2. Contrary to the case of traveling waves studied in ^[17,18], the winding number of the Evans function is (approximately) 1. This is a numerical evidence for the presence of an unstable real simple eigenvalue in the interval $(0, 20)$. Indeed, if $\nu$ would be a root of $E(\lambda) = 0$ with nonzero imaginary part, then by (2.17) $\overline{\nu}$ would also be a root, and hence the winding number would be greater than 1. Therefore, there exists a real root with multiplicity 1.

Let us compute $E(\lambda)$ on the real line. Equation (2.17) implies that $E(\lambda) \in \mathbb{R}$ for $\lambda \in \mathbb{R}$. We initialize the Kato ODE at $\lambda = 20$ and integrate it along the real line until $\lambda = 10^{-6}$ using $2 \cdot 10^4$ points. The result is plotted in Figure 3. The Evans function has a real root $\lambda_0 \approx 0.0496$.

The instability possibly is related to the presence of a supersonic region along the profile as will be discussed in Section 4.

3.   Quantum hydrodynamics with nonlinear viscosity

In this section we consider the following quantum hydrodynamics system with nolinear viscosity

Let us rewrite system (3.1) in terms of $(\rho, u)$ variables as

where the enthalpy $h(\rho)$ is

see ^[10]. We are searching for standing wave profiles of the form

such that

Substituting (3.3) into (3.2) we obtain

Integrating Eq (3.4) up to $\pm \infty$ we get

where $A = R^+ U^+$. Substituting (3.6) into (3.5) and integrating we obtain

where

with

Theorem 6.1 in ^[24] provides a characterization for the existence of homoclinic loops to (3.7). These homoclinic loops correspond to standing wave profiles for (3.2).

3.1.   Numerical computation of profiles

We use the method from Section 2.2, that is we compute the profile numerically using a boundary value solver. As an initial guess we use a function of the form (2.5) with $c = 0.5$,

We compute the profile for parameters

Since $|U^+| = 0.9 < 1.22 \approx c_s(R+)$, the condition for existence of profile of Theorem 6.1 in ^[24] is verified; see Figure 4. The flow along the profile is supersonic in the region where the following inequality holds:

For the parameter values (3.8) the right hand side of (3.9) is 0.78 and the flow is supersonic for $|x| < 2.29$.

3.2.   Linearization

Here we perform a linearization of (3.2) around a standing wave profile. Denoting by $(\rho, u)$ the deviation from the profile $(R, U)$, we obtain the following full linearized operator:

where

The associated eigenvalue problem reads

The asymptotic operator at $\infty$ is given by

The eigenvalue problem associated to $\mathcal{L}_{\infty}$ is

Let us rewrite it as a first-order system

where

3.3.   Essential spectrum

The characteristic equation $\det (\nu Id - M^+) = 0$ is

Setting $\nu = i \xi$, $\xi \in \mathbb{R}$, in (3.11) and dividing by $2/k^2$, we obtain the dispersion relation:

The stability condition for the essential spectrum is identical to Theorem 2.2. Therefore, standing wave profiles for the system (3.2) always have stable essential spectrum.

3.4.   System in integrated variables

Let us now consider the integrated variables

Expressing $\rho$ and $u$ in terms of $\hat{\rho}$ and $\hat{u}$ in (3.10) and integrating from $-\infty$ to $x$ we obtain the system in integrated variables:

with

Let us rewrite (3.12)-(3.13) as a first-order system

where $\hat{V} = [ \hat{\rho}, \hat{u}, \hat{u}_1, \hat{u}_2]^{\top}$, with $\hat{\rho}' = \hat{u}_1$, and $\hat{u}_1' = \hat{u}_2$, and

with

3.5.   Numerical evaluation of the Evans function

In this section we compute the Evans function for the QHD system with nonlinear viscosity (3.2). We use the method from Section 2.6. We integrate the ODE $\phi' = (B(x, \lambda)-\mu^-)\phi$ on $[-L_1, 0]$ and $\phi' = (B(x, \lambda)-\mu^+)\phi$ on $[0, L_1]$ backwards, where $L_1 = 40$ as before, however here $B(x, \lambda)$ is the compound matrix of $\hat{M}(x, \lambda)$ defined in (3.14) and $\mu^{\pm}$ is the stable (unstable) eigenvalue of $B$ at $+\infty$ with minimal (maximal) real part.

We evaluate the Evans function on the contour surrounding a semi-annular region with radii 20 and $10^{-6}$ with center at the origin and vertical segment on the imaginary axis. We integrate the reduced Kato ODE using $10^6$ points, then we compute $E(\lambda)$ with the solver ode15s in Matlab with relative tolerances $10^{-8}$ and $10^{-10}$, the latter being used if $E(\lambda)$ is evaluated wtih $\lambda$ close to 0. The Evans function $E(\lambda)$ is ploted in Figure 5, and in Figures 6 and 7 with enlargement of a region around the origin. Notice that the Evans function $E(\lambda)$ is small in absolute value for $\lambda \approx 0$ along the contour, hence a higher accuracy is required for its numerical computation than the one used for the QHD system with linear viscosity (cf. Section 2.6). The winding number of the Evans function is (approximately) 1. As before, this is numerical evidence for the presence of a simple real eigenvalue in the interval $(0, 20)$.

Now, let us evaluate $E(\lambda)$ on the real line. We initialize the Kato ODE at $\lambda = 20$ and integrate it along the real line until $\lambda = 10^{-6}$ using $4 \cdot 10^5$ points. The result is depicted in Figure 8. The Evans function $E(\lambda)$ has a real root $\lambda_0 \approx 0.0026$. Notice that both the viscosity coefficient and the unstable eigenvalue are smaller in the nonlinear viscosity case, compared to their counterparts in case with the linear viscosity (cf. the parameter set (2.7)), even though increasing the viscosity tends to stabilize the system. Therefore, our numerics suggests that the nonlinear viscosity has stabilizing effect, corroborating the analytical result in [9, Remark 4.7].

4.   Discussion and conclusions

It was shown numerically in Sections 2.6 and 3.5 that the standing wave profiles considered for the QHD systems (1.1) and (3.2) with linear and nonlinear viscosity are spectrally unstable. The essential spectrum is stable, however there is a simple real unstable eigenvalue. This numerical result is in contrast with the numerical stability results for traveling wave profiles in ^[17,18] where it was shown that non-monotone traveling wave profiles are spectrally stable. We notice that in both cases there exists a supersonic region, where the velocity along the profile is larger than the speed of sound (see Figures 1 and 4). On the other hand, one can show numerically that the velocity along the profiles is subsonic everywhere for the cases considered in ^[17,18]. Therefore, we make the following:

Conjecture 4.1. A standing or traveling wave profile for the QHD systems with linear or nonlinear viscosity (1.1) and (3.2) is spectrally stable if and only if the velocity is subsonic or sonic along the profile.

In other words, the conjecture states that a profile $(R, U)$ is spectrally stable if and only if

If this inequality holds, then taking limits as $x \rightarrow \pm \infty$ we obtain $|U^{\pm}| \leq c_s(R^{\pm})$, that is the end states are subsonic or sonic, in accordance with the necessary and sufficient condition in Theorem 2.2 for stability of the essential spectrum. Let us now compare the numerical results for the QHD systems with linear and nonlinear viscosity. The unstable eigenvalue is larger for the former system even though the viscosity coefficient is smaller for the latter (cf. the parameter sets (2.7) and (3.8)). Therefore, our numerical results suggest that the nonlinear viscosity has stabilizing effect.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research work was conducted in Gran Sasso Science Institute, University of L'Aquila and Universidad Nacional Autónoma de México.

Conflict of interest

The author declares no conflict of interest.