Numerical spectral analysis of standing waves in quantum hydrodynamics with viscosity

We study the spectrum of the linearization around standing wave profiles for two quantum hydrodynamics systems with linear and nonlinear viscosity. The essential spectrum for such profiles is stable; we investigate the point spectrum using an Evans function technique. For both systems we show numerically that there exists a real unstable eigenvalue, thus providing numerical evidence for spectral instability.


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 ρ ≥ 0 is the density, m = ρu is the momentum, where u is the fluid velocity, p(ρ) = ρ γ with γ ≥ 1 is the pressure, t ≥ 0 and x ∈ R.Moreover, µ, k > 0 are viscosity and dispersive coefficients, respectively.The function ( √ ρ) xx / √ ρ 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 µm xx is replaced with µρ(m x /ρ) 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 ρ(t, x) = P (x), m(t, x) = J(x), such that lim The first studies of models with dispersive terms were [22] and [12]; see also [11], [20] and [13].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] and [18] in the case of non-monotone shocks, and analytically in [8] and [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 (Theorem 4.1 in [24]).Let us denote by c s (ρ) = γρ γ−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 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 equation (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] and [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 equation (2.2) vanishes.Susbstituting equation (2.3) in (2.2) and proceeding as in [16, Section 2], we conclude that the density profile P solves the ODE where 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 steadystates 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 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 The solid line is the density profile, while the dashed line corresponds to sonic states.In the region where P (x) is below the dashed line the flow is supersonic.
Figure 1).This type of reduced density pulse is known as dark soliton (see, for instance [7]).

2.3.
Linearization.Let us denote by (ρ, m) the deviation from (P, J).We obtain the full linearized operator around the profile where .
The associated eigenvalue problem reads (2.9) 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 The associated eigenvalue problem is Let us rewrite it as a first order system , and the limit matrix M + is given by (2.10) The characteristic equation of M + is det(νId − M + ) = 0.It reads ) 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 Σ j = {λ = λ j (ξ) ∈ C : ξ ∈ R}, j = 1, 2, where λ j (ξ) are determined by the roots of (2.12).Let us now define stability of essential spectrum.
Definition 2.1.The essential spectrum of L is stable, if it is contained in the closed left half-plane, that is The essential spectum of 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.(i) 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 λ is to the right of the curves Σ j .2.5.System in integrated variables.Following [16] (see 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 ρ and m in terms of ρ and m and integrating from −∞ to x we get ) where Let us rewrite the system (2.13)-(2.14)as V ′ = M (x, λ) V , where V = [ρ, m, û1 , û2 ] ⊤ , ρ′ = û1 , û′ 1 = û2 , and The limiting matrix of M (x, λ) as x → ±∞ is (2.10).Since the profile (P, J) is a pulse, the limits of M (x, λ) at ±∞ coincide and are equal to the matrix M + given by (2.10).We assume that M + is hyperbolic.This is always true for λ to the right of the curves Σ j .Denote by ν − 1 , ν − 2 and by ν + 1 , ν + 2 the eigenvalues of M + with positive and negative real parts, respectively, and indicate with v ± i the corresponding (normalized) eigenvectors.Let Then, the Evans function can be defined by To numerically compute the Evans function, we use the compound matrix method (see [14] and [17], [18]).This method is used in order to get a stable numerical procedure, in spite of the fact that the system Y ′ = M (x, λ)Y is numerically stiff.Specifically, the compound matrix B(x, λ) is given by: be the piecewise linear interpolant of (x 1 , ζ 1 ), ..., (x N , ζ N ).We obtain the matrix B(x, λ) using ζ(x).Similarly we integrate the equation φ ′ = (B(x, λ) − µ + )φ on [0, L 1 ] backwards, where this time µ + is the stable eigenvalue of B at +∞ with minimal (negative) real part.Then, the coefficients µ ± 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 φ ± = (φ ± 1 , . . ., φ ± 6 ) as follows: .
To compute the initial conditions we integrate the reduced Kato ODE where P ± denotes the spectral projection of B± = lim x→+∞ B(x, λ) corresponding to µ ± .This choice of initial condition provides an analytic Evans function E(λ).
To numerically integrate (2.16) we use the algorithm from [25], that is |r 1 ± | = 1 eigenvector as before (referring to maximal/minimal decay/growth rate of B± ) and for k > 0, Then, provided E(λ) does not vanish on a closed contour Γ which does not intersect the essential spectrum of 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 λ = 0 and vertical segment on the imaginary axis.Furthermore, we do not evaluate the Evans function at 0, but evaluate it up to ±i10 −6 .Along the contour we integrate the Kato ODE using 5 • 10 4 points.Then, we compute E(λ) with the solver ode45 in Matlab with relative tolerance 10 −6 and we set L 1 = 40.Finally, we apply the symmetry of (2.17) The Evans function E(λ) is plotted in Figure 2. Contrary to the case of traveling waves studied in [17] and [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 ν would be a root of E(λ) = 0 with nonzero imaginary part, then by (2.17) ν 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(λ) on the real line.Equation (2.17) implies that E(λ) ∈ R for λ ∈ R. We initialize the Kato ODE at λ = 20 and integrate it along the real line until λ = 10 −6 using 2 • 10 4 points.The result is plotted in Figure 3.The Evans function has a real root λ 0 ≈ 0.0496.
The instability possibly is related to the presence of a supersonic region along the profile as will be discussed in Section 4.

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 (ρ, u) variables as where the enthalpy h(ρ) is [10].We are searching for standing wave profiles of the form Integrating equation (3.4) up to ±∞ we get where A = R + U + .Substituting (3.6) into (3.5) and integrating we obtain where 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 ≈ 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.

Linearization.
Here we perform a linearization of (3.2) around a standing wave profile.Denoting by (ρ, 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 ∞ is given by The eigenvalue problem associated to L ∞ is Let us rewrite it as a first-order system where Setting ν = iξ, ξ ∈ 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.Expressing ρ and u in terms of ρ and û in (3.10) and integrating from −∞ to x we obtain the system in integrated variables: ) 2 with Let us rewrite (3.12)-(3.13)as a first-order system V ′ = M (x, λ) V , where V = [ρ, û, û1 , û2 ] ⊤ , with ρ′ = û1 , and û′ 1 = û2 , and 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 φ ′ = (B(x, λ) − µ − )φ on [−L 1 , 0] and φ ′ = (B(x, λ) − µ + )φ on [0, L 1 ] backwards, where L 1 = 40 as before, however here B(x, λ) is the compound matrix of M (x, λ) defined in (3.14) and µ ± is the stable (unstable) eigenvalue of B at +∞ 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(λ) with the solver ode15s in Matlab with relative tolerances 10 −8 and 10 −10 , the latter being used if E(λ) is evaluated wtih λ close to 0. The Evans function E(λ) 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(λ) is small in absolute value for λ ≈ 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(λ) on the real line.We initialize the Kato ODE at λ = 20 and integrate it along the real line until λ = 10 −6 using 4 • 10 5 points.The result is depicted in Figure 8.The Evans function E(λ) has a real root λ 0 ≈ 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] (see Remark 4.7).

Discussion
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] and [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] and [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 |U (x)| ≤ c s (R(x)), x ∈ R. If this inequality holds, then taking limits as x → ±∞ we obtain |U ± | ≤ c s (R ± ), 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.

2. 6 .
Numerical evaluation of the Evans function.To define the Evans function, let us consider the equation Y ′ = M (x, λ)Y with M (x, λ) defined in (2.15).

Figure 2 .Figure 3 .
Figure 2. The image of a semi-circular contour with radius 20 through the Evans function E(λ) for the QHD system with linear viscosity.The origin is marked in red.

Figure 5 .
Figure 5.The image of a semi-annular region with radii 20 and 10 −6 through the Evans function E(λ) for the QHD system with nonlinear viscosity.The origin is marked in red.

Figure 6 .
Figure 6.Same as Figure 5 with enlargement of a region around the origin.

Figure 7 .Figure 8 .
Figure 7. Same as Figures 5 and 6 with further enlargement of a region around the origin.