Dispersive shocks in Quantum Hydrodynamics with viscosity
Corrado Lattanzio Pierangelo Marcati Delyan Zhelyazov



In this paper we study existence and stability of shock profiles for  a 1-D compressible Euler system in the context of Quantum Hydrodynamic models. The dispersive term is originated by the quantum effects described through the Bohm potential; moreover we introduce a  (linear) viscosity to analyze its interplay with the former while proving existence, monotonicity and stability of travelling waves connecting a Lax shock for the underlying Euler system.
The existence of monotone profiles  is proved  for sufficiently small shocks;  while the case of large shocks leads to  the (global) existence for an oscillatory profile, where dispersion plays a significant role. 
The spectral analysis of the linearized problem about a profile is also provided. In particular, we derive a sufficient condition for the stability of the essential spectrum and we estimate the maximum modulus of the eigenvalues in the unstable plane, using a careful analysis of the Evans function.



A Multi-Band Catalog of 10978 Star Clusters, Associations and Candidates in the Milky Way
    Eliade F. Lima
    December 30, 2023
=========================================================================================






§ INTRODUCTION

The aim of this paper is to study how dissipation interacts with dispersion in terms of existence and stability of traveling wave profiles, or dispersive shocks. We consider Quantum Hydrodynamics with a linear viscosity term

    {[                                             ρ_t+m_x=0,; m_t+(m^2/ρ+ρ^γ)_x=ϵμ m_xx+ϵ^2 k^2 ρ((√(ρ))_xx/√(ρ))_x, ].

where ρ=ρ(t,x)>0, m=m(t,x), γ≥ 1, 0<ϵ≪1, μ>0, k>0. Here ϵμ and ϵ^2 k^2 are the viscosity and dispersive coefficients, respectively. Moreover ρ^γ is the pressure and we consider t≥0 and x ∈ℝ. The  dispersive term is due to the Bohm potential and the resulting system is used for instance in superfluidity or to model semiconductor devices. 

The first studies concerning dispersive terms can be found
by <cit.>; see also <cit.>, and  <cit.> (and the reference therein) for a quite complete analysis via the Whitham modulation theory.
The first attempt to analyze the spectral theory of the linearized operator around dispersive shocks has been discussed in  <cit.> regarding the case of  p-system with real viscosity and linear capillarity, while the mathematical theory  for  the 
Quantum Hydrodynamics   can be found in <cit.>.
In the present paper, we study in particular the effects of the combination of dispersion and dissipation effects in terms of existence and stability of profiles for such models. In addition, we are able to discuss  the spectrum of the linearized operator also  in the case of non monotone shocks for our model in Eulerian coordinates; we also underline here that   a more detailed numerical description of the behavior of the Evans function close to the zero eigenvalue is presented in the companion paper <cit.>.

 We first  present a local result, concerning the study of profiles for sufficiently small shocks for the underlying Euler system, where both viscosity and capillarity terms are neglected. The existence result is proved by means of  a bifurcation argument, where the bifurcation parameter is the difference between the end states of the shock. After this quite standard result, we focus on our main interest, namely the existence and stability of profiles for large shocks, showing in particular the combined effect of dissipation (coming from the viscosity term) and dispersion (coming from the capillarity term).
 Even if our result will require the viscosity coefficient to be “dominant”, our results include cases when the  dispersion plays a significant role, giving rise to the existence of oscillatory profiles. This result is achieved by showing that the dynamical system solved by the profile possesses an invariant region if the dissipation is sufficiently large.

Then we focus on the study of stability properties of such profiles, also taking advantage of the numerical tests described in full details in the companion paper <cit.>.
We start by investigating the spectrum of the linearized system about a profile: we analyze the essential spectrum and show that it is stable for subsonic or sonic end states. Moreover we give a bound for the real parts of the eigenvalues, by using an energy estimate, and exclude the presence of eigenvalues with non–negative real part for  |λ| sufficiently big. This last result is not obtained only via an asymptotic result as |λ|→ +∞, but rather we give an explicit estimate for the constant  which bounds from above the modulus of possible eigenvalues. This last more “quantitative” result is  fundamental to explicitly localize the region of the unstable half plane where eigenvalues may lie and then to numerically analyze such region, giving numerical evidence of spectral stability; see <cit.> for further details.

The paper is organized as follows. In the Section <ref> we recall basic facts on the underlying Euler system to then derive the second order equation satisfied by a traveling wave profile for  <ref>. In Section <ref> we present the (local and global) existence results for the profile. Finally, the   last section is devoted to the study of the linearized system about the profile and to present the results about spectral stability. 



§ PROFILE EQUATION
 

Let us start by recalling some well known facts  concerning the Euler system


    {[          ρ_t+m_x=0,; m_t+(m^2/ρ+ρ^γ)_x=0 ].

which will be used later for the analysis of the Quantum Hydrodynamics system.
The flow velocity is denoted by u=m/ρ. Let U=(ρ,m). The eigenvalues of the Jacobian (characteristic speeds) are 

    λ_1(U)=u-c_s(ρ),λ_2(U)=u+c_s(ρ),

where

    c_s(ρ)=√(γρ^γ-1)=√(d (ρ^γ)/dρ)

denotes the sound speed. A shock wave with end states ρ^± and m^± and shock speed s satisfies the Rankine-Hugoniot conditions:

    m^+-m^-=s(ρ^+-ρ^-),
    (m^2/ρ+ρ^γ)^+-(m^2/ρ+ρ^γ)^-=s(m^+-m^-).

A k-shock satisfies the Lax entropy condition:

    λ_k(U^+)<s<λ_k(U^-).



Let us consider a traveling wave profile

    ρ=P(x-s t/ϵ)m=J(x-s t/ϵ)

As customary,  the speed  s ∈ℝ  of the travelling wave and its limiting end states 

    lim_y→±∞P(y)=P^± and lim_y→±∞J(y)=J^±

satisfy the Rankine–Hugoniot conditions (<ref>)-(<ref>).

We rewrite the Bohm potential in   conservative form

    ρ((√(ρ))_xx/√(ρ))_x=1/2(ρ (lnρ)_xx)_x

After substituting the profiles P and J in the system (<ref>) and multiplying by ϵ we obtain

    - s P'+J'   =0,
    
    -s J' + (J^2/P+P^γ)'   =μ J”+k^2/2(P(ln P)”)'

where ' denotes d/dy and P=P(y), J=J(y). 







Integrating equation (<ref>), we get 





    J(y)-sP(y)=J^–sP^-

We can also integrate (<ref>) from y to +∞ to get




    J(y)-sP(y)=J^+-sP^+

So we obtain 

    J(y)=sP(y)-A,

 with

    A=s P^±-J^±,

as follows from the Rankine-Hugoniot condition (<ref>).

Substituting the expression for J(y) into equation (<ref>) and integrating we get

    ∫_-∞^y( -s(sP(x)-A)' + ((sP(x)-A)^2/P(x)+P(x)^γ)' )dx 
    
     = ∫_-∞^y (μ sP(x)”+k^2/2(P(x)(ln P(x))”)' )dx.

We can also integrate from y to +∞. We obtain the planar ODE

    P”=2/k^2 f(P) - 2 s μ/k^2 P' + P'^2/P,

where

    f(P)   =-s(s P-A)+(sP-A)^2/P+P^γ-B
       =P^γ-(As+B)+A^2/P.

Here the constant B is given by

    B=-s J^±+(J^2/P+P^γ)^±.

Let P'=Q, we get

    P'=Q   =f_1,
    
    Q'=2 f(P)/k^2-2μ s/k^2 Q+Q^2/P   =f_2.

The constants A,B in f(P) can be expressed in terms of P^±:

    f(P)   = P^γ +P^-P^+/P(P^+)^γ-(P^-)^γ/P^+-P^-
        -(P^+)^γ+1-(P^-)^γ+1/P^+-P^-.

We have

    f”(P)=γ (γ-1)P^γ-2+2 P^- P^+/P^3(P^+)^γ-(P^-)^γ/P^+-P^-,

so f”(P)>0 as a sum of a non-negative and a positive term.

The function f(P) has two zeros P^±. The system (<ref>)-(<ref>) has two equilibria [P^-,0] and [P^+,0]. 

Suppose P^+<P^- and s>0. The jacobian evaluated at the equilibria is

    M =  [             0             1; 2 f'(P^±)/k^2    -2 μ s/k^2 ].

Since f”(P)>0, we have f'(P^+)<0 and f'(P^-)>0.

The eigenvalues of M at P^+ are

    λ_1,2(P^+)=-s μ±√(2 k^2 f'(P^+) + s^2 μ^2)/k^2.

Since f'(P^+)<0 we have either 2k^2f'(P^+)+s^2μ^2<0 or 2k^2f'(P^+)+s^2μ^2 ≥ 0 and √(2k^2f'(P^+)+s^2 μ^2)<s μ. Hence λ_1(P^+), λ_2(P^+) <0. If 2k^2f'(P^+)+s^2μ^2<0 the eigenvalues have nonzero imaginary parts.

The eigenvalues of M at P^- are

    λ_1,2(P^-)=-s μ±√(2 k^2 f'(P^-) + s^2 μ^2)/k^2.

At P^- we have f'(P^-)>0, so √(2 k^2 f'(P^-)  + s^2 μ^2)> s μ. Therefore λ_1(P^-)>0 and λ_2(P^-)<0.


§ EXISTENCE OF SHOCK PROFILES

Here we prove existence of shock wave profiles for the system (<ref>), connecting the end states P^± and J^±, which satisfy the Rankine-Hugoniot conditions (<ref>)-(<ref>). We consider both existence of sufficiently small shocks, and possibly oscillatory profiles for large shocks, where the dispersion plays a significant role. The ratio μ/k controls how oscillatory the shocks are.



 §.§ Local existence of profiles

In this section we consider local existence of profiles using a bifurcation theory argument about the variable P^+ -P^- =h. In particular, proving  a transcritical bifurcation at h=0, we obtain existence of profiles for small shocks without  restrictions on the coefficients μ and k.

We start by rewriting f in (<ref>) in the variables u=P-P^- and P^+ -P^- =h as follows:

    f̃(u,h)   =(P^-+u)^γ +P^-(P^-+h)/P^-+u(P^-+h)^γ-(P^-)^γ/h
        -(P^-+h)^γ+1-(P^-)^γ+1/h.

Let u'=v. The system (<ref>)-(<ref>) becomes

    u'=v=f̃_1(u,v,h),
    
    v'=2/k^2f̃(u,h)-2 s μ/k^2v+v^2/P^-+u=f̃_2(u,v,h).

In the following lemma we transform the system (<ref>)-(<ref>) to a normal form, which is given by  two scalar decoupled equations, for which we can prove the existence of a heteroclinic connection.

For any P^->0 there is an ϵ>0, such that if s>0 (s<0), |P^+-P^-|<ϵ and P^+<P^- (P^-<P^+), there exists a heteroclinic for the system (<ref>)-(<ref>), connecting [P^-,0] to [P^+,0].


Suppose h< 0, s>0. Let A(u,0,h)= D(f̃_1,f̃_2)/D(u,v)|_(v=0) be the jacobian evaluated for v=0,  h< 0, and any u, that is

    A(u,0,h)=  [                 0                 1; 2/k^2∂f̃(u,h)/∂ u        -2 μ s/k^2 ].

Then, system (<ref>)–(<ref>) admits two equilibria [0,0] and [h,0], corresponding to [P^-,0] and [P^+,0].
The calculation in the end of Section <ref> shows that A(0,0,h) has two eigenvalues with negative real parts, and A(h,0,h) has a positive and a negative eigenvalue. 
Therefore, the equilibrium [0,0] is a saddle and [h,0] is stable.

Now, in order to study this system close to the bifurcation value h=0, 
we  extend f̃ and its derivatives by continuity:

    f̃(0,0)   =0,∂f̃(0,0)/∂ u=0,∂f̃(0,0)/∂ h=0,
    ∂^2 f̃(0,0)/∂ u ∂ h   =-γ(γ+1)/2(P^-)^γ-2<0,
    ∂^2 f̃(0,0)/∂ u^2   =γ(γ+1)(P^-)^γ-2>0.

As a consequence, we get

    A=A(0,0,0)=  [          0          1;          0 -2 μ s/k^2 ],

and clearly its eigenvalues   are λ_1=0, λ_2 = -2 s μ /(k^2).
Let K be the matrix of column eigenvectors of A, namely

    K=  [          1 -k^2/2 μ s;          0          1 ].

 We change the variables (u,v)=K(w_1,w_2), that is we transform the original system according to the eigenbasis at [0,0,0]:

    [ w_1; w_2 ]'=
    [          0          0;          0 -2 s μ/k^2 ][ w_1; w_2 ]
    +
    [ k^2/2 s μ g;           g ],

where

    g=g(w_1-k^2/2 s μw_2, w_2,h),
    
    g(u,v,h)=2/k^2f̃(u,h)+v^2/P^-+u.

In view of (<ref>), we have

    g(0,0,0)=0,∂_u g(0,0,0)=0,∂_v g(0,0,0)=0,

so that (<ref>) is of the form of a linear part and a perturbation.

There is a center manifold w_2=ψ(w_1,h), which satisfies the tangency conditions:

    ψ(0,0)=0,∂ψ(0,0)/∂ w_1=0,∂ψ(0,0)/∂ h=0.

We perform center manifold reduction:

    w_1'=ζ(w_1,h),


where the function ζ(w_1,h) is given by the following expression:

    ζ(w_1,h): =k^2/2 s μ g(w_1-k^2/2 s μψ(w_1,h),ψ(w_1,h),h).

Then we have

    ζ(0,0)   =k^2/2 s μ g(0,0,0)=0,∂ζ(0,0)/∂ w_1=0,
    ∂ζ(0,0)/∂ h   =0, a=∂^2 ζ(0,0)/∂ h ∂ w_1,
    2b=∂^2 ζ(0,0)/∂ w_1^2.

We get the normal form of a transcritical bifurcation

    w_1' = a h w_1 + b w_1^2,

where

    a =-γ(γ+1)/2 s μ (P^-)^γ-2,
    b= -a.

Indeed, Theorem 5.4 from <cit.>, p. 159 implies that the system (<ref>)-(<ref>) is locally topologically equivalent to the scalar nonlinear ODE (<ref>) for w_1 augmented with the linear equation w_2'=sgn(λ_2) w_2 = - w_2.
The nondegenerancy condition a,b≠ 0 is satisfied, so the parameter h unfolds the bifurcation with normal form given in (<ref>). Since λ_2<0, the center manifold is attracting. We have a<0, b>0, so the trivial equilibrium w_1=0 is unstable and the nontrivial equilibrium h<0 is stable for (<ref>).
Correspondingly, (<ref>) has the unstable (trivial) equilibrium w_1=0, a stable (nontrivial) equilibrium w_1^0(h) and, as a consequence, 
there exists a heteroclinic connecting w_1 to w_1^0(h). The homeomorphism preserves the number of eigenvalues with positive (negative) real parts of the equilibria. The decoupled system with first equation (<ref>) has a trivial equilibrium [0,0] with one positive and one negative eigenvalue, and a non-trivial equilibrium [w_1^0(h),0] with two negative eigenvalues. So the equilibrium [0,0] is mapped to [P^-,0] and  [w_1^0(h),0] is mapped to [P^+,0]. The heteroclinic, connecting [0,0] to [w_1^0(h),0], corresponds to a heteroclinic, connecting [P^-,0] to [P^+,0].

Now let s<0. We have λ_2>0, so the center manifold is not attracting. We have a>0 and b<0. Let P^+>P^- (that is h>0). The equilibrium w_1=0 is unstable, and w_1^0(h)=h>0 is stable. The decoupled system has a heteroclinic, connecting [0,0] to [w_1^0(h),0]. The equilibrium [0,0] has two positive eigenvalues, and it is  mapped to [P^-,0]. The equilibrium [w_1^0(h),0] has one positive and one negative eigenvalue, and it is mapped to [P^+,0]. The heteroclinic of the decoupled system corresponds to a heteroclinic, connecting [P^-,0] to [P^+,0].

The heteroclinic orbit for P constructed  in Lemma <ref> gives a traveling wave profile for our model, by using J(y) from  equation (<ref>).
The following corollary phrases this existence result of profiles in terms of Rankine-Hugoniot  and   Lax entropy conditions of end states, as well as super- or sub-sonic conditions.


For any P^->0, there is an ϵ>0, such that if P^+>0, |P^+-P^-|<ϵ, and the end states P^±, J^± and the speed s satisfy the Lax condition for a 2-shock  with supersonic right state and u^+>c_s(P^+) or a 1-shock with a subsonic left state, then there is a traveling wave profile.


Let P^±,J^±,s satisfy the Rankine-Hugoniot conditions. We are going to express the end states J^± in terms of P^± and s, using the Rankine-Hugoniot conditions. In this way we obtain an equivalent expression.

Expressing J^+ from (<ref>) we get

    J^+=J^-+s(P^+-P^-).

Substituting J^+ into (<ref>) and dividing by the quadratic coefficient

    P^–P^+/P^-P^+≠ 0,

we get the quadratic equation

    (J^-)^2-2sP^- J^-+P^-P^+((P^+)^γ-(P^-)^γ-P^-s^2 +((P^-)^2s^2)/(P^+))/P^–P^+=0.

The quadratic equation (<ref>)  has two solutions J^-_1,2=s P^--A_1,2, where

    A_1,2=∓√(P^+ P^-)√((P^+)^γ-(P^-)^γ/P^+-P^-)

Substituting these solutions in equation (<ref>) yields two solutions J^+_1,2=s P^+-A_1,2.



Suppose the shock satisfies the Lax condition

    λ_2(U^+)<s<λ_2(U^-).

Then λ_2(U^+)<s, so u^++c_s(P^+)<s. Suppose J^-=s P^--A_1, with A_1<0. This implies J^+=s P^+-A_1. Then using u^+=J^+/P^+, we get -A_1/P^++c_s(P^+)<0. Since c_s(P^+)>0, we get a contradiction. So J^-=s P^--A_2, where A_2>0, which implies J^+=s P^+-A_2.

From the Lax condition we get λ_2(U^+)<λ_2(U^-), hence u^++c_s(P^+)<u^-+c_s(P^-). Using u^±=J^±/P^±, we get

    A_2(P^+-P^-)/P^+ P^-<c_s(P^-)-c_s(P^+).

Now, suppose P^+>P^-. Since the sound speed c_s(P) is nondecreasing, we get c_s(P^-)-c_s(P^+) ≤ 0. On the other hand (A_2(P^+-P^-))/(P^+ P^-)>0. So we get a contradiction. Hence P^+<P^-.

Suppose we have a 2-shock with a supersonic right state and u^+>c_s(P^+). Then the Lax condition implies s>0. Therefore, we have the condition of Lemma <ref> for a local existence of a profile.



Suppose we have a 1-shock:

    λ_1(U^+)<s<λ_1(U^-).

In this case s<λ_1(U^-). Suppose J^-=s P^--A_2, with A_2>0. This implies J^+=s P^+-A_2. Then using u^-=J^-/P^- we get A_2/P^-<-c_s(P^-). Since c_s(P^-)>0, this is a contradiction. Hence J^-=s P^--A_1, where A_1<0, which implies J^+=s P^+-A_1.

We have from the Lax condition λ_1(U^+)<λ_1(U^-), so u^+-u^-<c_s(P^+)-c_s(P^-). Using u^±=J^±/P^±, we get

    A_1(P^+-P^-)/P^+ P^-<c_s(P^+)-c_s(P^-).

Suppose P^+<P^-. The sound speed c_s(P) is nondecreasing, and

    A_1(P^+-P^-)/P^+ P^->0,

so we get a contradiction. Hence we have P^+>P^-.

Suppose we have a 1-shock with a subsonic left state. Then |u^-|<c_s(P^-), and in particular u^-<c_s(P^-). The Lax condition implies s<0. Hence, we have the second condition of Lemma <ref> for a local existence of a shock profile.




 §.§ Global existence of profiles

This section concerns the proof of existence of the profiles in the case of large amplitude shocks. In contrast to the case of local existence, here we shall need conditions between  viscosity and dispersion coefficients, and in particular the latter needs to be sufficiently strong, see 
  Lemma <ref>.
  
We start by rescaling  the velocity Q in terms of k as Q̃=k^2 P'/2. Then the system can be rewritten as follows: 

    P'   =2/k^2Q̃=:f̃_1,
    Q̃'   =f(P)-2 s μ/k^2Q̃+2/k^2Q̃^2/P=:f̃_2.

The crucial observation is that the   reduced (indeed conservative) system 

    P'=2/k^2Q̃,
    Q̃'=f(P)

admits a   homoclinic loop, which confines the heteroclinic connection  we are looking for. 
 More specifically, the homoclinic loop turns out to be the boundary of an invariant region for the full system.
 Finally, to exclude closed trajectories  inside such invariant region, we shall use the Poincaré–Bendixon criterion. 
 
 
 We start with the case s>0 and P^+<P^-.


This condition is implied by the condition for a 2-shock with a supersonic right state, and u^+>c_s(P^+).

We have

    P^2 f'(P)=γ P^γ+1-A^2=0

has a unique solution P_0=(A^2/γ)^1/γ+1. Since f'(P)=γ P^γ-1-P^-2A^2, for P ∈ (0,P_0), P^γ+1<A^2/γ, so f'(P)<0. Moreover f'(P)>0 for P>P_0. We are going to consider choices of parameters, for which f(P)=0 has one or two positive solutions. If there are two solutions, they will be the limiting values of the traveling wave profile P^+<P^-. We have P^+< P_0<P^-. Also we have 
f”(P)>0.


Let us now consider  the related  reduced system (<ref>)-(<ref>), or its second order counterpart: 

    k^2/2P”=f(P),

where we truncate the last two terms in (<ref>). 






Let

    F(P)=∫ f(z) dz = P^γ+1/γ+1-(As+B)P+A^2 ln P,

The system (<ref>) has energy

    H(P,Q̃)=F(P)-Q̃^2/k^2-F(P^-).

Since F(P)-F(P^-)=∫_P^-^P f(z)dz and f(P)<0 for P^+<P<P^-, we have ∫_P^-^P f(z)dz>0 for P^+<P<P^-. Also f(P)>0 for 0<P<P^+, and lim_P ↓ 0F(P)=-∞, so there is a point P^⋆∈ (0,P^+), such that F(P)-F(P^-)>0 for P^⋆<P<P^-, and F(P)-F(P^-)<0 for 0<P<P^⋆. So we have H(P,0)<0 for 0<P<P^⋆, and H(P,0)>0 for P^⋆<P<P^-.
The system (<ref>) has a homoclinic loop starting at P^-, which corresponds to H(P,Q̃)=0. Moreover the energy levels, contained inside the homoclinic loop, are compact and correspond to H(P,Q̃)=C>0.

Let

    F_1(P)=F(P)-(s μ P/k)^2-F(P^-),

be the restriction of H on the line Q̃=s μ P.

Now, let s<0 and P^-<P^+. The function f(P) has two roots P^±. We are going to consider the reverse parameter ỹ=-y and let ξ, η correspond to P, -Q̃, respectively.  We are going to prove existence of a heteroclinic [ξ,η], connecting [P^+,0] to [P^-,0], which corresponds to a heteroclinic, connecting [P^-,0] to [P^+,0]. 
We rewrite the system (<ref>)-(<ref>) in terms of ỹ:

    ξ'   =2/k^2η,
    η'   =f(ξ)-2 s̃μ/k^2η+2/k^2η^2/ξ,

where s̃=-s, s̃>0.
In the new variables the reduced system <ref> becomes

    ξ'   =2/k^2η,
    η'   =f(ξ),

with '=d/dỹ.

The system (<ref>)-(<ref>) has energy

    H̃(ξ,η)=F(ξ)-η^2/k^2-F(P^+).

The system (<ref>)-(<ref>) has a homoclinic loop, starting at [P^+,0], contained in the set H̃(ξ,η)=0.

Let

    F̃_1(ξ)=F(ξ)-(s μξ/k)^2-F(P^+).

Let s>0 and P^+<P^-. If F_1(P)≤ 0 for 0<P<P^+, then there is a traveling wave profile for the system (<ref>)-(<ref>), connecting [P^-,0] to [P^+,0]. If in addition

    s μ/k<√(-2 f'(P^+)),

then the traveling wave profile is non-monotone.

Let s<0 and P^-<P^+. If F̃_1(ξ)≤ 0 for 0<ξ<P^-, then there is a traveling wave profile for the system (<ref>)-(<ref>), connecting [P^-,0] to [P^+,0]. If in addition

    -s μ/k<√(-2 f'(P^-)),

then the traveling wave profile is non-monotone.


We are going to show that the homoclinic loop of (<ref>)-(<ref>) is confining, the orbit, contained in the unstable manifold of the saddle [P^-,0] is inside it, and its ω-limit set is [P^+,0] (see Figure <ref>).

Let [P(y),Q̃(y)] be a solution of (<ref>)-(<ref>) and ℋ(y)=H(P(y),Q̃(y)). Then

    ℋ' = ∂ H/∂ PP'+∂ H/∂Q̃Q̃'=4 Q̃^2/k^4(s μ - Q̃/P).

If ℋ' ≥ 0 for any point of the homoclinic loop, then it is confining, that is  a trajectory that starts in the homoclinic loop will stay inside it for all y ≥ 0. Also s μ P - Q̃≥ 0 implies ℋ' ≥ 0 . We require that the homoclinic loop is contained in the region s μ P ≥Q̃. We would like to prove that the trajectory, contained in the unstable manifold of the saddle [P^-,0] will converge to [P^+,0] as y → +∞.


We can show that the ω-limit set of the trajectory, contained in the unstable manifold of the right equilibrium is the left equilibrium, if the left equilibrium is stable and we can exclude loops. For this we can use the Poincaré-Bendixon criterion, which states if f ∈ C^1(E), where E is a simply connected region in ℝ^2, and if there exists a function B ∈ C^1(E), such that the divergence of the vector field B f, ∇· (B f) is not identically zero and does not change sign in E, then the planar system x'=f(x) has no closed orbits, lying entirely in E (see <cit.>, p.265, Theorem 2). Under the conditions of this criterion, there are no separatrix cycles or graphics of x'=f(x), lying entirely in E (see <cit.>, p.265).

The divergence of f̃/P is 

    ∂/∂ P(f̃_1/P)+∂/∂ Q(f̃_2/P)=2/k^2 P^2(Q̃-s μ P),

therefore it does not change sign and is not identically zero in the same region, where H' ≥ 0.

Now let us consider the jacobians at the equlibria. Let a=2/k^2>0, b=2 s μ/k^2>0. Let

    J_1=  [       0       a; f'(P^±)       0 ],
         J_2=  [       0       a; f'(P^±)      -b ].

Here J_1 is the linearization of (<ref>) at P^±, Q̃=0 and J_2 is the linearization of (<ref>) at P^±, Q̃=0. The eigenvalues of J_2 are (-b ±√(b^2+4 a f'(P^±)))/2. At P^+, since f'(P^+)<0, we have either b^2+4 a f'(P^+)<0, or b^2+4 a f'(P^+)≥ 0 and √(b^2+4 a f'(P^+))<b. So in any case λ_1, λ_2 <0. If b^2+4 a f'(P^+)<0 the eigenvalues have nonzero imaginary parts.

Now we will show that the vector, which is tangent to the unstable manifold of the saddle [P^-,0] is directed inside the homoclinic loop. At P^- f'(P^-)>0, so √(b^2+4 a f'(P^-))>b. Therefore we have a positive and a negative eigenvalue - a saddle. The eigenvector, corresponding to (-b + √(b^2+4 a f'(P^-)))/2 is 

    v_2=  -[ (b + √(b^2+4 a f'(P^-)))/2 f'(P^-);                                  1 ]

The eigenvector, corresponding to the positive eigenvalue of J_1 is 

    ṽ_2= -[ √(a/f'(P^-));            1 ]

If ṽ_2,1>v_2,1, then the eigenvector v_2, which is tangent to the unstable subspace is pointing inside the homoclinic loop H(P,Q̃)=0. This is true, because 2 b^2+2b√(b^2+4af'(P^-))>0, which implies (b+√(b^2+4af'(P^-)))^2>a/f' (2 f'(P^-))^2 and ṽ_2,1>v_2,1 follows by taking a square root.



Now we would like to consider the sufficient condition for the existence of an oscillatory profile. From (<ref>) the jacobian at the equilibrium [P^+,0] has imaginary eigenvalues if and only if (<ref>) holds.

Let us express Q̃ as a function of P from (<ref>) and take the positive branch:

    Q̃=g(P)=k√(F(P)-F(P^-)).

Its derivative is

    g'(P)=k/2f(P)/√(F(P)-F(P^-)).

So g'(P) vanishes if and only if f(P)=0. The positive branch of the homoclinic has a maximum at P=P^+. So if the line Q̃=s μ P does not intersect the homoclinic loop in the interval (0,P^+], it will not intersect it also in the interval (0,P^-). If for all P in the interval (0,P^+] the expression (<ref>), which the restricion of H on the line Q̃=s μ P is nonpositive, then the line Q̃=s μ P does not intersect the homoclinic loop (since H(P,Q̃)>0 inside the loop).

Now, let use consider the case s<0.
The system (<ref>)-(<ref>) has the same form as (<ref>)-(<ref>) and the above proof applies to it. So it has a heteroclinic, connecting [P^+,0] to [P^-,0], which corresponds to a heteroclinic, connecting [P^-,0] to [P^+,0] in the parameter y.

We can show numerically that Lemma <ref> applies to the parameters

    A=1, B=1.1, s=1,γ=3/2,μ=1, k=√(2).

 These parameters correspond to a non-monotone profile (see <cit.>).
 The minimum value of μ for a given k, for which we can guarantee a heteroclinic is  the value, for which (s μ P)^2=k^2(F(P)-F(P^-)) has a unique solution. This is the maximum μ, for which this equation has a solution and then the tangency condition is

    k f(P)/2√(F(P)-F(P^-))=s μ.


 For γ=1, the condition of Lemma <ref> can be verified analytically. The derivative of F_1(P) is 

    F_1'(P)=f(P)-2(s μ/k)^2 P.

We have F_1”(P)=f'(P)-2(s μ/k)^2<0 for 0<P<P^+, so F_1'(P) is monotonically decreasing for 0<P<P^+. Also lim_P ↓ 0 F_1'(P)=+∞, F_1'(P^+)=-2(s μ/k)^2<0, hence F_1'(P) has one zero in the interval (0,P^+).

If γ∈ℕ, then P(f(P)-2(s μ/k)^2P) is a polynomial.
For γ = 1, let

    D=(As+B)^2k^4-4A^2k^2(k^2-2s^2μ^2).
 
The roots of P(f(P)-2(s μ/k)^2P)=0 are

    (As+B)k^2±√(D)/2(k^2-2s^2μ^2)

Consider the zero of F_1'(P) in the interval (0,P^+). At this zero F_1 will have a maximum. The condition of Lemma <ref> will be verified, if F_1(P)≤ 0 at the zero. This condition in conjunction with (<ref>) guarantees that there is an oscaillatory profle. 	We can write similar formulas for γ = 2,3.


 The vector, tangent to the unstable manifold of the saddle [P^-,0] is directed inside the homoclimic loop, correspoding to H(P,Q̃)=0. So for any trajectory [P(y),Q̃(y)], contained in the unstable manifold there will be y_0, such that H(P(y_0),Q̃(y_0))> ϵ_0>0. The forward trajectory will be contained in the energy level H(P,Q̃)≥ϵ_0, which contains the only steady-state [P^+,0].

For fixed values of γ (γ∈{1,3/2,5/3}) there are non-empty intervals for μ, for which a oscillatory heteroclinic exists.



§ LINEARIZATION AND STABILITY RESULTS



Using the change of variables τ=t/ϵ, y=(x-st)/ϵ, u(x,t)=u̅(τ,y), we get the full linearized operator around the profile for (<ref>):

    L
    [ ρ̃; J̃ ]
    =[                                                s ρ̃_y - J̃_y; s J̃_y +(J^2/P^2ρ̃)_y-(2 J/PJ̃)_y-γ (P^γ-1ρ̃)_y+μJ̃_yy+L_Vρ̃ ],

where

    L_V ρ̃ = k^2/2ρ̃_yyy-2 k^2 ( (√(P))_y(ρ̃/√(P))_y)_y,

with associated eigenvalue problem given by

    λ[ ρ̃; J̃ ] = L[ ρ̃; J̃ ].


A related constant coefficient linear operator is clearly obtained linearizing our original system about a constant state, thus obtaining the same operator, but for s=0.
 Denote

    α=J̅^2/P̅^2-γP̅^γ-1β= -2 J̅/P̅.

Then the operator, corresponding to the linearization around the constant steady–state (P̅, J̅) is

    L_c
    [ ρ̃; J̃ ]
    =[                    - J̃'; αρ̃'+βJ̃'+μJ̃”+k^2/2ρ̃”' ],

where '=d/dy.
The asymptotic operators  at ±∞ for (<ref>)  are given by

    L_±∞[ ρ̃; J̃ ]
    =[                  s ρ̃' - J̃'; α^±ρ̃'+β^±J̃'+μJ̃”+k^2/2ρ̃”' ]

where 

    α^±=(J^±)^2/(P^±)^2-γ (P^±)^γ-1β^±=s-2 (J^±)/P^±.

Finally, wee may rewrite the equation

    λ[ ρ̃; J̃ ] = L_±∞[ ρ̃; J̃ ].

as a first order system V'=M^±V, with V=[ρ̃,J̃,u_1,u_2]^T.



 §.§ Essential spectrum

As it is well known, the spectrum of L consists of two   parts: the essential spectrum and the point spectrum; we start here by investigating the former.

To this end, we recall that the dispersion relation can be found from det(i ξ Id-M)=0:

    λ^2+ξ(μξ - i (s+β^±))λ+ξ^2(-α^±+k^2 ξ^2/2-s(β^±+iμξ))=0.

To simplify notation, in what follows, we are going to drop the superscript of α and β.


The essential spectrum is not always stable.
The characteristic equation of M is det(ν Id - M)=0, that is

    ν^4+2 s μ/k^2ν^3+2/k^2(α+s β - λμ)ν^2-2s+β/k^2λν + 2 λ^2/k^2=0.


The condition for a subsonic steady-state is |u|<c_s( P), , which becomes
α<0 after squaring. The condition α^±≤ 0 corresponds to end states that are either subsonic or sonic.


If μ^2 ≠ 2 k^2 (generically) as long as λ is to the right of the curve λ(ξ), solving (<ref>), we have 2 roots with positive real parts and 2 roots with negative real parts of (<ref>) that is we have consistent splitting. Moreover if α^±≤ 0, then the bound for the essential spectrum is in the closed left half-plane. Also if ξ≠ 0, then λ_1,2<0.


Suppose λ∈ℝ, λ≫ 1. The Discriminant of (<ref>) is

    Δ = (512/k^10μ^4-2048/k^8μ^2+2048/k^6)λ^6+𝒪(λ^5).

Also 

    D = 64 a^4 (ŝ-q̂^2/4)=(-64/k^4μ^2+128/k^2)λ^2+𝒪(λ),

where a=1 is the fourth order coefficient of (<ref>), ŝ and q̂ are the zeroth and second order coefficients of the depressed quartic equation, associated to (<ref>).

If μ^2<2 k^2, then for sufficiently large λ, D>0 and if μ^2 ≠ 2 k^2, then for λ≫ 1, Δ>0, since the leading order terms in the expansion in λ have the appropriate signs. Therefore (<ref>) has two pairs of complex conjugated roots, that are not real. In particular the roots are simple.

Since

    P=q̂=-16 μ/k^2λ + 𝒪(1),

in all cases P<0, and if μ^2>2 k^2, then D<0, hence the four roots are real and distinct in the regime of real large positive λ.



Now we are going to apply the Descartes' rule of signs in the case of real roots (λ≫ 1). Suppose first that s+β<0. Then the number of sign differences between consecutive coefficients is 2, hence there are at most 2 positive roots. If we substitute ν→ -ν, then we have again two sign changes, so there are at most 2 negative roots. Suppose that the number of positive roots is less than 2. Then it has to be 0, because it is less than the upper bound by an even number. Then there must be 4 negative roots. This is a contradiction with the upper bound. Hence the polynomial has 2 negative and 2 positive roots.

Suppose s+β>0. Then the number of sign chages is 2. In the case ν→-ν the number of sign changes
is again 2. Hence similarly as before we get that the polynomial has 2 positive and 2 negative roots.

Now, suppose s+β=0. Then there are 2 sign changes. This is the case also for ν→ -ν. Hence again the polynomial has 2 positive and 2 negative roots.



Now consider the case μ^2<2 k^2 for real λ≫ 1. Suppose (<ref>) does not have a purely imaginary root. The Routh stability criterion is a necessary and sufficient condition for the existence of roots only in the left half-plane. If this is the case, all the coefficients of (<ref>) must be positive. However the second order coefficient is negative. Hence there are roots in the open right half-plane. If we substitute ν→ -ν the second order coefficient is again negative, therefore (<ref>) must have roots also in the left half-plane. Hence there are two complex conjugate roots in the left half-plane and two complex conjugated roots in the right half-plane.

Suppose we have a polynomial ν^4 + a_1 ν^3 + a_2 ν^2 + a_3 ν + a_4 with all a_i ≠ 0, a_i - real. If this polynomial has a purely imaginary root, then a_3/a_1>0 and a_3^2+a_1^2 a_4 = a_1 a_2 a_3.

Let s + β>0. Then all a_i ≠ 0, a_1>0 and a_3<0, hence a_3/a_1<0, hence (<ref>) does not have a purely imaginary root.

Now, let s+β<0. In this case a_i ≠ 0 and a_3/a_1>0. However in the regime λ≫ 1, which we are considering:

    a_3^2+a_1^2 a_4 - a_1 a_2 a_3 = 4/k^6(k^2(s+β)^2+2 s^2μ^2-2 s(s+β)μ^2)λ^2+𝒪(λ).

The coefficient of λ^2 is sum of three positive terms, therefore it is positive. Hence a_3^2+a_1^2 a_4 ≠ a_1 a_2 a_3. Therefore also in this case (<ref>) does not have a purely imaginary root.

Now, suppose s+β=0. It follows from the factorization (ν^2+p)(ν^2+qν+r), with p ≥ 0 that ν^4+a_1 ν^3 + a_2 ν^2 + a_4=0, with a_i ≠ 0, a_i real, does not have a purely imaginary root. If it would have a purely imaginary root, then this would imply p=0 or q=0. But we would have a_4=0 in the first case and a_1=0 in the second. We have a_1=2 s μ/k^2 ≠ 0 and a_4=2 λ^2/k^2 ≠ 0, so this is not the case.

It follows that (<ref>) does not have a purely imaginary root.



Now we are going to derive a sufficient condition for the stability of the essential spectrum. It applies to constant steady-state and also to the stability of the asymptotic steady-states as y →±∞. Note that the constant steady-state is not always stable.

The roots of the dispersion relation (<ref>) are

    λ_1,2=1/2(-μξ^2 + i (s + β)ξ±√(D))

where

    D = (ξ (μξ -i (s+β)))^2-4 ξ^2 (-α + k^2 ξ^2/2 -s (β + i μξ ))=p+i q,

with

    p    = -s^2 ξ^2 +4 αξ^2 + 2 s βξ^2-β^2 ξ^2 -2 k^2 ξ^4 + μ^2 ξ^4,
    
    q    = 2 s μξ^3 - 2 βμξ^3.

The condition

    -μξ^2 + |√(D)| ≤ 0

 guarantees that λ_1,2 are in the closed left half-plane. Since μξ^2 ≥ 0, it is equivalent to
 
    (√(D))^2 ≤μ^2 ξ^4.

 We have
 
    √(D) = √(2)/2√(√(p^2+q^2)+p)

The condition (<ref>) is equivalent to
 
    √(p^2+q^2)≤ 2 μ^2 ξ^4 - p.

 Let
 
    A = 2 μ^2 ξ^4 - p = ξ^2 (s-β)^2-4 αξ^2 + 2 k^2 ξ^4 + μ^2 ξ^4.

 If α≤ 0, then A ≥ 0 for all ξ∈ℝ as a sum of four non-negative terms. Also, if ξ≠ 0, A>0. If A ≥ 0, then (<ref>) is equivalent to A^2 - p^2 - q^2 ≥ 0. Let B = A^2 - p^2 - q^2. Then (<ref>) is equivalent to B ≥ 0. We have
 
    B/8 μ^2=-2 αξ^6 + k^2 ξ^8.

 Therefore if α≤ 0, A ≥ 0 and B ≥ 0. Also if ξ≠ 0, B>0. Hence (<ref>) holds.




 §.§ Point spectrum for a profile

In contrast with the situation described in the proposition above, the localization of the point spectrum of the linearized operator along a profile is more involved (note that our operatori is not self adjoint), and in particular the energy estimate in Lemma <ref> is proved only for (λ) sufficiently big.
Thus, to locate the point spectrum in that case, an efficient method is to locate   the zeros of the Evans function,  the latter being exactly the eigenvalues of the operator under consideration. The argument needed requires a careful analysis of the behavior of such function  for large |λ|, which gives a quantitative version of  the asymptotic results of <cit.>, excluding the presence of eigenvalues for |λ|> C with an explicit bound for the constant (see Lemma <ref>). To give the definition of the Evans function we shall use later on, we first rewrite  here below our problem in itegrated variables.


  §.§.§ System in integrated variables

For the analysis of the eigenvalue problem (<ref>) we will need the Evans function and to this end it is also important to re-express the above linearized system in terms of integrated variables, 
because  this transformation removes the zero eigenvalue (always present, being its eigenfunction given by the derivative of the profile), without further modifications of the spectrum; see, for instance <cit.>.
To this end, consider

    ρ̂(x)=∫_-∞^x ρ̃(y)dy,Ĵ(x)=∫_-∞^x J̃(y) dy,
 
Integrating the equation (<ref>) it follows that for λ≠ 0 the integrated variables ρ̃ and J̃ decay exponentially as |x|→ +∞. Expressing ρ̃ and J̃ in terms of ρ̂ and Ĵ, and integrating (<ref>) from -∞ to x we get the system in integrated variables:

    λρ̂   = s ρ̂'-Ĵ',
    λĴ   = f_1 ρ̂' + f_2 Ĵ'+μĴ”+k^2/2ρ̂”'-2 k^2 (√(P))'(ρ̂'/√(P))',

with

    f_1(x)    = J(x)^2/P(x)^2-γ P(x)^γ-1,
    
    f_2(x)    = s - 2 J(x)/P(x).

We can rewrite (<ref>)-(<ref>) as V'=M̂(x,λ)V, where V=[ρ̂,Ĵ,u_1,u_2]^T and 

    M̂(x,λ) = [                                         0                                         0                                         1                                         0;                                        -λ                                         0                                         s                                         0;                                         0                                         0                                         0                                         1;                               2 λ f_2/k^2                                   2 λ/k^2 2 λμ/k^2-2 f_1/k^2-2 s f_2/k^2-(P')^2/P^2                          2 P'/P-2 s μ/k^2 ].

The limit of M̂(x,λ) as x →±∞ is given by

    M^± = [                       0                       0                       1                       0;                      -λ                       0                       s                       0;                       0                       0                       0                       1;              2 β^±λ/k^2                 2 λ/k^2 2/k^2(μλ - s β^± - α^±)              -2 s μ/k^2 ],






  §.§.§ The Evans function

To define  the Evans function, let us consider the equation Y'=M̂(y,λ)Y, where M̂(y,λ) is defined in (<ref>).
As it is manifest, its limits at ±∞ are given by  the matrices M^±, defined by (<ref>), corresponding to limit states P^±, and we assume  these matrices are hyperbolic. This is always true if we are to the right of the bound for the essential spectrum. In addition, we assume that M^- has k unstable eigenvalues ν^-_1,…,ν^-_k (i.e. (ν^-_i)>0), and M^+ has n-k stable eigenvalues ν^+_1,...,ν^+_n-k (i.e. (ν^+_i)<0), and denote the corresponding (normalized) eigenvectors by v^±_i. In our case n=4 and k=2. Let Y^-_i be a solution of Y'=M(y,λ)Y, satisfying exp(ν^- y)Y^-(y) tends to v^-_i as y → -∞ and exp(ν^+ y)Y^+(y) tends to v^+_i as y → +∞. Then,  the Evans function can be defined by

    E(λ) = det(Y^-_1(0), .., Y^-_k(0),  Y^+_1(0), ... ,Y^+_n-k(0)).

As a consequence, a point λ∈ℂ is in the point spectrum of L if and only if E(λ)=0.


  §.§.§ Estimate for (λ) for variable coefficients

Now we consider the eigenvalue equation (<ref>) with variable coefficients. 
We are going to derive a bound of the real part of the eigenvalues (λ) using an energy estimate. 


For each eigenvalue of (<ref>) with variable profile we have
 (λ) ≤max(C_1,C_2, 2 C_3/k^2), where

    C_1 = 1/2 ϵ_1+M_2/ϵ_2+M_4/ϵ_3,
    C_2 = ϵ_2/2+M_3,
    C_3 = M_1/ϵ_2+M_5/ϵ_3,

with ϵ_2>0, ϵ_1/2+ϵ_3/2<μ, and explicit constants M_i.


We multiplying the first equation of (<ref>) by ρ̅̃̅ and the second equation of (<ref>) by J̅̃̅.
Let ρ̃=ρ̃_̃r̃+i ρ̃_̃ĩ and J̃=J̃_̃r̃+i J̃_̃ĩ. 
We have

    -α s (ρ̃'ρ̃)   =-α s(ρ̃_̃r̃' ρ̃_̃r̃+ρ̃_̃ĩ' ρ̃_̃ĩ)=-α s/2(|ρ̃|^2)',
    β(J̃'J̃)    = β/2(|J̃|^2)'

By integration by parts we obtain 

    ∫μ(J̃_̃r̃”J̃_̃r̃+J̃_̃ĩ”J̃_̃ĩ)dy = - μ∫ ((J̃_̃r̃')^2+(J̃_̃ĩ')^2)dy=-μ∫ |J̃'|^2 dy.

We have 

    (ρ̃”'J)=ρ̃_̃r̃”'J̃_̃r̃+ρ̃_̃ĩ”'J̃_̃ĩ,

hence by integration by parts

    ∫k^2/2(ρ̃”'J̃) dy=-k^2/2∫ (J̃_̃r̃' ρ̃_̃r̃”+J̃_̃ĩ' ρ̃_̃ĩ”)dy

Substituting J̃_̃r̃' and J̃_̃ĩ' from the first equation of (<ref>) yields

    -k^2/2∫ (J̃_̃r̃' ρ̃_̃r̃”+J̃_̃ĩ' ρ̃_̃ĩ”)dy 
       =-k^2/2∫( (s ρ̃_̃r̃'- (λ) ρ̃_̃r̃
       +(λ) ρ̃_̃ĩ)ρ̃_̃r̃”+(s ρ̃_̃ĩ'- (λ) ρ̃_̃ĩ-(λ) ρ̃_̃r̃)ρ̃_̃ĩ”)dy
       =-k^2/2∫ (t_1+t_2+t_3)dy,

where

    t_1    = s(ρ̃_̃r̃”ρ̃_̃r̃'+ρ̃_̃ĩ”ρ̃_̃ĩ')=s/2((ρ̃_̃r̃')^2+(ρ̃_̃ĩ')^2)',
    
    t_2    = (λ) (ρ̃_̃ĩρ̃_̃r̃”-ρ̃_̃ĩ”ρ̃_̃r̃),
    
    t_3    = -(λ) (ρ̃_̃r̃ρ̃_̃r̃”+ρ̃_̃ĩρ̃_̃ĩ”).

Therefore

    -k^2/2∫ t_1 dy = 0.

Using integration by parts we get

    -k^2/2∫ t_2 dy = -k^2 (λ)/2∫ (-ρ̃_̃r̃'ρ̃_̃ĩ'+ρ̃_̃ĩ'ρ̃_̃r̃')dy=0,

and again by integration by parts

    -k^2/2∫ t_3 dy = -(λ) k^2/2∫( (ρ̃_̃r̃')^2+(ρ̃_̃ĩ')^2) dy=-(λ) k^2/2∫ |ρ̃'|^2 dy,

which yields 

    ∫k^2/2(ρ̃”'J̃) dy = -(λ) k^2/2∫ |ρ̃'|^2 dy.

Taking the real part of s  ρ̃' ρ̃ and integrating we get

    ∫(s ρ̃' ρ̃)dy=0

by (<ref>). Moreover

    ∫(-J̃' ρ̃)dy ≤1/2 ϵ_1∫ |ρ̃|^2 dy + ϵ_1/2∫ |J̃'̃|^2 dy

by Young inequality. From the second equation of (<ref>) we have

    ∫(s J̃' J̃)dy=0

by (<ref>). Also

    ∫(((J^2/P-γ P^γ-1) ρ̃)' J̃) dy 
       ≤∫(ϵ_2/2 |J̃|^2 +1/2 ϵ_2(|(f_1 ρ̃_̃r̃)'|^2 + |(f_1 ρ̃_̃ĩ)'|^2) )dy 
       ≤∫(ϵ_2/2 |J̃|^2 +1/ϵ_2(|f_2|^2 |ρ̃|^2+(|f_1|^2 |ρ̃'̃|^2))dy 
       ≤ϵ_2/2∫ |J̃|^2 dy + M_2/ϵ_2∫ |ρ̃|^2 dy+M_1/ϵ_2∫ |ρ̃'|^2 dy,

where f_1 = J^2/P-γ P^γ-1, f_2 = f_1', M_1,2 = sup_y ∈ℝf_1,2 (y)^2. Further

    ∫((-2 J/PJ̃)' J̃)dy   =-∫(2 J/P)(J̃_̃r̃J̃_̃r̃'+J̃_̃ĩJ̃_̃ĩ')dy 
    
    =∫(J/P) (|J̃|^2)'dy    =-∫(J/P)'|J̃|^2dy ≤ M_3 ∫ |J̃|^2 dy,

by integration by parts, where f_3 = (J/P)' and M_3 = sup_y ∈ℝ|f_3 (y) |. Moreover

    ∫(-2k^2 ((√(P))'(ρ̃/√(P))' )' J̃)dy
        =- ∫ 2k^2 (√(P))' ((ρ̃_̃r̃/√(P))' J̃_̃r̃' +(ρ̃_̃ĩ/√(P))' J̃_̃ĩ' ) dy 
       ≤∫(|2 k^2 (√(P))' (ρ̃_̃r̃/√(P))' J̃_̃r̃'| +|2 k^2 (√(P))' (ρ̃_̃ĩ/√(P))' J̃_̃ĩ'| ) dy 
       =∫(|f_4 ρ̃_̃r̃ + f_5 ρ̃_̃r̃'||J_r'|+|f_4 ρ̃_̃ĩ + f_5 ρ̃_̃ĩ'||J̃_̃ĩ'| )
       ≤ϵ_3/2∫ |J̃'|^2 dy + 1/2 ϵ_3∫(|f_4 ρ̃_̃r̃ + f_5 ρ̃_̃r̃'|^2+|f_4 ρ̃_̃ĩ + f_5 ρ̃_̃ĩ'|^2 )dy  
       ≤ϵ_3/2∫ |J̃'|^2 dy + 1/ϵ_3∫ f_4^2 |ρ̃|^2 dy + 1/ϵ_3∫ f_5^2 |ρ̃'|^2 dy 
       ≤ϵ_3/2∫ |J̃'|^2 dy + M_4/ϵ_3∫ |ρ̃|^2 dy+M_5/ϵ_3∫ |ρ̃'|^2 dy,

where f_4 = k^2 P' (√(P))'/P^3/2, f_5 = 2 k^2 (√(P))'/√(P), M_4,5 = sup_y ∈ℝf_4,5 (y)^2. Using (<ref>) and (<ref>) and collecting the inequalities we get

    ((λ)-C_1)∫ |ρ̃|^2 dy + ((λ)-C_2)∫ |J̃|^2 dy 
       + (k^2/2(λ)-C_3)∫ |ρ̃'|^2 dy ≤ (-μ + ϵ_1/2+ϵ_3/2)∫ |J̃'|^2 dy
 
Hence if ϵ_1/2+ϵ_3/2 < μ and (λ)>max(C_1,C_2, 2 C_3/k^2) the left hand side of (<ref>) is positive and the right hand side is negative, which is a contradiction. Therefore (λ) ≤max(C_1,C_2, 2 C_3/k^2).









  §.§.§ The Evans function for large |λ|

Now we are going to show that we can consider a scalar equation to locate the eigenvalues.

Let us consider the eigenvalue problem with linear dispersion:

    λρ̃   = s ρ̃'-J̃' ,
    λJ̃   = (f_1(x) ρ̃)' + (f_2(x) J̃)'+μJ̃” + k^2/2ρ̃”'

with

    f_1(x)    = J(x)^2/P(x)^2-γ P(x)^γ-1, 
    
    f_2(x)    = s - 2 J(x)/P(x)

Consider also

    ρ̂^(4)+2 s μ/k^2ρ̂^”' + 2/k^2(f_1 + f_2 s - λμ)ρ̂^”+2/k^2(f_1'+s f_2'-λ f_2 - s λ)ρ̂' 
    
    +2/k^2(-λ f_2'+λ^2)ρ̂=0

and the eigenvalue problem with nonlinear dispersion:

    λρ̃   = s ρ̃'-J̃',
    λJ̃   = (f_1(x) ρ̃)' + (f_2(x) J̃)'+μJ̃” + f_3(x) ρ̃ + f_4(x) ρ̃' 
       + f_5(x) ρ̃”+ k^2/2ρ̃”',

which is equivalent to (<ref>),
where f_1(x) and f_2(x) are given by (<ref>), and

    f_3(x)    = k^2 (P'(x)P”(x)/P(x)^2 - (P'(x)/P(x))^3 ), 
    
    f_4(x)    = k^2 (3/2(P'(x)/P(x))^2-P”(x)/P(x)), 
    
    f_5(x)    = -k^2 P'(x)/P(x),

and the scalar equation

    ρ̂^(4) + 2/k^2(s μ + f_5)ρ̂^”'+2/k^2(f_1 + f_2 s + f_4 -λμ) ρ̂^”
       +2/k^2(f_1'+f_2' s-λ f_2+f_3-s λ) ρ̂^'
    +2/k^2(λ f_2'+λ^2)ρ̂=0.

In the proof of the following lemma we rescale the variable y.

If λ≠ 0, the eigenvalue equations (<ref>)-(<ref>),  (<ref>)-(<ref>) are  equivalent to (<ref>) and  (<ref>), respectively. In particular, if λ≠ 0 is not an eigenvalue of (<ref>)((<ref>)), it is also not an eigenvalue of (<ref>)-(<ref>)((<ref>)-(<ref>)). The Evans function for (<ref>)((<ref>))  does not vanish for (λ)≥ 0 and |λ| large enough.


Suppose that λ≠ 0. Integrating (<ref>)-(<ref>) we get

    ∫ρ̃ dx=0,∫J̃ dx = 0.

We shall use the integrated variable

    ρ̂(x)=∫_-∞^x ρ̃(y)dy.

We have that ρ̃(x) decays exponentially as |x|→ +∞. In particular 
|ρ̃(x)
| ≤ C_1 exp(-C_2 x). So

    |ρ̂(x)|   =|∫_x^∞ρ̃(y)dy| ≤∫_x^∞ |ρ̃(y)|dy ≤ C_1 ∫_x^∞exp(-C_2 y) dy 
       =C_1/C_2exp(-C_2 x).

Therefore ρ̂(x) decays exponentially as x → +∞. Similarly ρ̂(x) decays exponentially as x → -∞.



Integrating equation (<ref>) from -∞ to x yields

    J̃ = s ρ̂'-λρ̂.

After expressing ρ̃ in terms of ρ̂, substituting J̃ from (<ref>) and dividing by k^2/2 we get (<ref>). If λ≠ 0  is not an eigenvalue of (<ref>), it is also not an eigenvalue of (<ref>)-(<ref>).



Now, we make a change of variable

    x = y/|λ|^1/2.

Dividing (<ref>) by |λ|^2 yields

    d^4 ρ̂/dy^4+2 s μ/k^2 |λ|^1/2d^3 ρ̂/dy^3
    + 2/k^2(f_1(y/|λ|^1/2) + f_2(y/|λ|^1/2) s/|λ| - λ/|λ|μ)d^2 ρ̂/dy^2
       +2/k^2 |λ|^3/2(f_1'(y/|λ|^1/2)+s f_2'(y/|λ|^1/2)-λ f_2(y/|λ|^1/2) - s λ)d ρ̂/dy
       +2/k^2(-λ/|λ|^2 f_2'(y/|λ|^1/2)+λ^2/|λ|^2)ρ̂=0

Denote λ̃ = λ/|λ|. Taking lim_|λ|→ +∞ in (<ref>) yields:

    d^4 ρ̂/dy^4 - 2 μ/k^2λ̃d^2 ρ̂/dy^2 +2/k^2λ̃^2 ρ̂=0.

Rewriting (<ref>) as a first-order system gives:

    d/dy[   ρ; u_1; u_2; u_3 ]
     = [         0         1         0         0;         0         0         1         0;         0         0         0         1;  -2/k^2λ̃         0 2 μ/k^2λ̃         0 ][   ρ; u_1; u_2; u_3 ].

The characteristic equation of (<ref>) is

    z^4 - 2 μ/k^2λ̃ z^2 + 2 λ̃^2/k^2=0.

Let (λ) ≥ 0. If 

    D = 4/k^2(μ^2/k^2-2)λ̃^2 ≠ 0,

that is μ^2/k^2≠ 2, then (<ref>) has 4 distinct roots, with (z_1),(z_2)<0 and 
(z_3),(z_4)>0.

The condition D ≠ 0 holds, when the dispersion and dissipation  terms do not exactly balance.

We make the change of variable w=z^2. Then (<ref>) becomes:

    w^2 - 2 μ/k^2λ̃ w + 2 λ̃^2/k^2=0

If (λ) ≥ 0, then (λ̃) ≥ 0. If D ≠ 0, since 2 λ̃^2/k^2≠ 0, the equation (<ref>) has two distinct nonzero roots. Suppose μ^2/k^2 < 2, that is the dispersion is dominating. The roots are:

    w_1,2 = μλ̃/k^2±i/k√(2-μ^2/k^2)λ̃.

Moreover |w_j| = √(2)/k, j=1,2. Let λ = exp(i θ), and θ∈ [0, π/2] ∪ [3 π/2, 2 π[. We have

    Arg(μ/k^2 + i/k√(2-μ^2/k^2))∈ ]0, π/2[.

If θ∈ [0, π/2], then Arg(w_1) ∈ ]0, π[. If θ∈ [3 π/2, 2 π[, then (w_1) > 0. Also

    Arg(μ/k^2 - i/k√(2-μ^2/k^2))∈ ]3 π/2, 2 π[.
 
If θ∈ [0, π/2], then (w_2) > 0. If θ∈ [3 π/2, 2 π[, then Arg(w_1) ∈ ]π, 2 π[. So w_1,2 are not negative. Hence z^2 = w_k has one solution with positive and one with negative real part for k=1,2. Therefore z_j are not purely imaginary. Also z_j are distinct, because distinct nonzero numbers cannot have equal square roots.

The equation (<ref>) has consant coefficients, so its Evans function can be computed.
Let z_k  be a simple eigenvalue of the matrix

    [         0         1         0         0;         0         0         1         0;         0         0         0         1;  -2/k^2λ̃         0 2 μ/k^2λ̃         0 ].

The associated eigenvector is v_k = [1,z_k,z_k^2,z_k^3]^T. The Evans function for (<ref>) is 

    Ẽ(λ) = ([v_1, v_2,v_3,v_4]) = ∏_j<k(z_j-z_k) 
       =(z_1-z_2)(z_1-z_3)(z_1-z_4)(z_2-z_3)(z_2-z_4)(z_3-z_4) ≠ 0.

We have that Ẽ(λ) ≠ 0, since the eigevalues z_k are distinct. The coefficients of (<ref>) and (<ref>) are uniformly in y close to each other, their Evans functions are uniformly close in λ as a consequence of Theorem 3.1, <cit.>. Therefore the Evans function for (<ref>) never vanishes for (λ)≥ 0, and |λ|>C, where C is some constant. So for any eigenvalue with (λ)≥ 0, we have |λ| ≤ C.



If the viscosity μ=0, choosing λ̃=i yields purely imaginary roots of (<ref>). Therefore the matrix (<ref>) is not hyperbolic.

We rewrite (<ref>) as (<ref>)
Expressing ρ̃ in terms of ρ̂, which decays exponentially as |x| → + ∞, integrating the first equation of (<ref>) and substituting J̃ in the second equation of (<ref>) we get (<ref>).
Making the change of variable (<ref>) in (<ref>) yields

    d^4 ρ̂/dy^4+2/k^2s μ + f_5/|λ|^1/2d^3 ρ̂/dy^3+2/k^2(f_1 + f_2 s + f_4/|λ| -μλ/|λ|)d^2 ρ̂/dy^2
       +2/k^2f_1'+f_2' s-λ f_2+f_3-s λ/|λ|^3/2d ρ̂/dy
    +2/k^2(λ f_2'/|λ|^2 +λ^2/|λ|^2)ρ̂=0,

where all the functions f_k and their derivatives are evaluated at y/|λ|^1/2. Taking limit as |λ| → +∞ in (<ref>) we obtain (<ref>). Therefore the same conclusion follows. 


  §.§.§ Estimate for the maximum of |λ|

We will estimate the constant C, such that each eigenvalue of (<ref>) with non-negative real part has absolute value less than C. In preparation to stating Lemma <ref>, let

    θ_1,2   = arg(μ/2± i ( 1 - μ^2/4)^1/2), 
    α   = (exp(i(θ+θ_1)/2))θ=π/2, 
    
    z_1,3   = ∓exp(i(θ+θ_1)/2),z_2,4 = ∓exp(i(θ+θ_2)/2).

The distance between the roots |z_j-z_k| does not depend on θ and we can compute it e.g. for θ = 0.

Also

    m_1(λ)    = |λ|^-1sup_x≥ 0|f_2'|, 
    
    m_2(λ)    = |λ|^-3/2sup_x≥ 0|f_1'+s f_2'+f_3| + |λ|^-1/2sup_x≥ 0|f_2+s|, 
    
    m_3(λ)    = |λ|^-1sup_x≥ 0|f_1+s f_2 + f_4|, 
    
    m_4(λ)    = |λ|^-1/2sup_x ≥ 0|s μ + f_5|.

Moreover

    r_j,k   = p_k/g_j, p_k = ∑_l=1^4 m_l(λ)|z_l|^l-1=∑_l=1^4 m_l(λ), 
    
    g_1    = |z_1-z_2||z_1-z_3||z_1-z_4|, g_2 = |z_1-z_2||z_2-z_3||z_2-z_4|, 
    
    g_3    = |z_1-z_3||z_2-z_3||z_3-z_4|, g_4 = |z_1-z_4||z_2-z_4||z_3-z_4|, 
    
    R_+   = [r_j,k].

The matrix R_- has the suprema in the definition of m_k(λ) taken for x ≤ 0.

    δ_± = ‖ R_±‖_F,ϵ_± = 4 α^-1δ_±.

In the following lemma we decompose the system into a constant coefficients part, which depends only on the direction λ/|λ|, and a perturbation, which becomes small for large |λ|. We use exponential dichotomies to estimate the difference between the Evans functions of the constant coefficient system and the perturbed system. We suppose for simplicity that k=√(2).

Suppose μ < 2. Let √(2)√(ϵ_+^2+ϵ_-^2)<1, and δ_±<α/4 for some |λ|=C. Then for all λ with (λ)≥ 0 and |λ|≥ C, the Evans function E(λ) for (<ref>) has no zeroes.


Consider the system

    du/dx=A(λ̃)u

where λ̃ = λ/|λ|. The matrix A(λ̃) does not depend on x and has simple eigenvalues z_j, j=1,...,4, with (z_1), (z_2) < 0 and (z_3), (z_4) > 0. We may fix λ̃. The system (<ref>) has an exponential dichotomy (see <cit.>, Chapter 4) on R_0^+. That is there are constants K, α and projection P such that

    X(t)PX^-1(s)_2 ≤ K e^-α (t-s), t ≥ s    ≥ 0,
    X(t)(Id-P)X^-1(s)_2 ≤ K e^-α (s-t), s ≥ t    ≥ 0,

where X(t) is the fundamental solutions matrix for (<ref>) with X(0)=Id. We consider the scalar product ⟨ x,y ⟩= x ·y̅. The vector norm is |x| = √(⟨ x, x ⟩). The 2-norm is A_2 = sup_|x|=1|A x|. Since the matrix A(λ̃) has constant coefficients, the constants K and α can be computed.

Now, consider the perturbed system

    du/dx=A(λ̃)u + B(x,λ) u.

We have

    lim_| λ | → +∞sup_x ≥ 0B(x,λ)_2 = 0.

More precisely, sup_x ≥ 0B(x,λ)_2 = 𝒪(|λ|^-1/2). Let δ = sup_x ≥ 0B(x,λ). If δ < α/(4 K^2), then the perturbed system (<ref>) also has an exponential dichotomy with projection Q. Moreover (see <cit.>, Chapter 4, Prop. 1)

    P-Q_2 ≤ 4 α^-1 K^3 δ.

Let P and Q be projections onto the subspaces M and N. There exist unique orthogonal projections P̃ and Q̃ onto M and N. Moreover (see <cit.> p.58, Theorem 6.35)

    P̃-Q̃_2 ≤ P - Q _2.

We have M =  N = 2. Denote ϵ = 4 α^-1 K^3 δ. Then P̃-Q̃_2 ≤ϵ. Let u_k, k=1,...,4 be the eigenvectors of A(λ̃). We suppose they are normalized, that is |u_k| = 1. We have u_k ∈ M, that is P̃ u_k = u_k. Denote h_k = u_k - Q̃ u_k. Then

    |h_k| =|u_k - Q̃ u_k | = | P̃ u_k - Q̃ u_k | ≤P̃ - Q̃_2 |u_k | ≤ϵ.

Let v_k = Q̃ u_k. We have ⟨ v_k, h_k ⟩ = ⟨Q̃ u_k, u_k - Q̃ u_k ⟩ = 0. Hence |u_k|^2 = |v_k|^2 + |h_k|^2. Therefore |v_k|^2 ≥ |u_k|^2 - ϵ^2 = 1-ϵ^2. There is an ϵ_0, such that ∀ϵ < ϵ_0, v_k ≠ 0. Also |⟨ v_1,  v_2 ⟩| ≤ |⟨ u_1, u_2 ⟩| + 2 ϵ + ϵ^2 and |v_1||v_2| ≥ 1- ϵ^2. So if 1-ϵ^2 > |⟨ u_1, u_2 ⟩| + 2 ϵ + ϵ^2, then |⟨ v_1,  v_2 ⟩| < |v_1| |v_2|. There is an ϵ_0, such that if ϵ < ϵ_0, this inequality holds. We showed that {v_1,v_2} is a basis of N, and N = span({v_1,v_2}).

We do the same computation for ℝ_0^-, and obtain the vectors v_3 and v_4.

Let E(λ) be the Evans function for (<ref>). Then E(λ) = ([u_1,u_2,u_3,u_4]). Let E_p(λ) be the Evans function for (<ref>). Then E_p(λ) = ([v_1,v_2,v_3,v_4]). Denote H = [h_1,...,h_4], U = [u_1,...,u_4]. Let A_F = √(∑_j,k |a_j,k|^2). The inequality A_2 ≤A_F holds. 
Suppose U is the identity matrix. Let q be an eigenvalue of U+H. The Bauer-Fike theorem implies that

    |1-q|≤H_2.

Since H_2 ≤H_F, H_F<1 shows that 0 cannot be an eigenvalue of U+H, hence det(U+H)≠ 0. Therefore E_p(λ)≠ 0.

Let s_j, j=1,...,4 be the eigenvectors of A(λ̃) and S = [s_1,s_2,s_3,s_4]. Let D(λ̃) = diag(z_1,z_2,z_3,z_4) = S^-1 A(λ̃)S. We make the change of variable u = S v. Then (<ref>) becomes

    dv/dx = D(λ̃)v + S^-1B(x,λ)S v.

Since D(λ̃) is diagonal its eigenvectors are u_j = e_j, the standard basis vectors, for j=1,...,4. Hence U = Id, (U) = U_2=U^-1_2=1, and κ = 1. Moreover P = diag(1,1,0,0) and K = 1.

We suppose for simplicity that k=√(2), although all the calculations can be done for arbitrary value of k.
We have 0<θ_1<π/2, -π/2 < θ_2 < 0 and θ_2 = - θ_1. Let λ̃ = exp(i θ). Then

    -π/4+θ_j/2≤θ + θ_j/2≤π/4+θ_j/2,j=1,2.

The roots of (<ref>) are z_k. We have α = min_k (z_k). From (<ref>) it follows that we can take α = (exp(i(θ+θ_1)/2)) for θ=π/2. Since α depends on θ, we can get a non-ciruclar region where the eigenvalues are contained by computing α for different θ.

We can directly compute S^-1B(x,λ)S. Since B(x,λ) has all elements except the last row equal to zero, that is b_4,k, we obtain 

    S^-1B(x,λ)S_2 ≤S^-1B(x,λ)S_F ≤R_F

and for the matrix  R from (<ref>),
where m_k(λ) are some upper bounds for b_4,k(λ), |b_4,k(x,λ)|≤ m_k(λ).
For equation (<ref>) we can choose m_k(λ) from (<ref>).
Note that the matrix R(λ) from (<ref>) has monotonically decreasing with |λ| elements. We get two matrices R_± corresponding to ℝ_0^+ and ℝ_0^- respectively, and two values for δ_± and ϵ_±. Moreover H_F ≤√(2)√(ϵ_+^2 + ϵ_-^2). 

Alternatively (see <cit.>, page 770, Corollary 2.14)

    |E_p(λ)-E(λ)| = |(U) - (U+H)| 
       ≤((κH_2/U_2+1)^n-1) |(U)| ≤((κH_F/U_2+1)^n-1) |(U)|,

where κ =  U _2  U^-1_2 and n=4.
To show E_p(λ) ≠ 0, it suffices that |E(λ)-E_p(λ)| <|E(λ)|. If E_p(λ) ≠ 0, then {v_1,v_2} (and {v_3,v_4}) are linearly independent sets, because every subset of a linearly independent set is linearly independent.

Since A_2 ≤√(n)A_∞,where A_∞ = max_j (∑_k |a_j,k|), we can estimate

    S^-1B(x,λ)S_2 ≤√(n)S^-1B(x,λ)S_∞
       ≤√(n)S^-1_∞B(x,λ)_∞S_∞.

Moreover S_∞ = 4. Moreover S^-1_∞ can be estimated by (see <cit.>)

    S^-1_∞≤max_1≤ k ≤ n∏_j=1,
    j≠ k^n 1+|z_j|/|z_k-z_j|.


For the set of parameters (<ref>), to which Lemma <ref> applies, we can show numerically using Lemma <ref> that there are no eigenvalues with (λ)≥ 0 and |λ| ≥ 1.9 · 10^4. By Evans function computations we can check numerically that there are no eigenvalues with (λ)≥ 0 and |λ| < 1.9 · 10^4 (see <cit.>). This is numerical evidence of point spectrum stability.

10

AM1  P. Antonelli, P. Marcati, 
On the finite energy weak solutions to a system in Quantum Fluid Dynamics , Comm. Math. Phys. 287, 657-686 (2009)

AM2 P. Antonelli, P. Marcati,  The Quantum Hydrodynamics system in two space dimensions, Arch. Ration. Mech. Anal. 203 , 499-527  (2012)

AMtf P. Antonelli, P. Marcati, Finite Energy Global Solutions to a Two-Fluid Model Arising in Superfluidity, Bull. Inst. Math. Acad. Sin. 10,349–373 (2015)

AMDCDS P. Antonelli, P. Marcati, Quantum hydrodynamics with nonlinear interactions, Discrete Contin. Dyn. Syst. Ser. S 9, 1–13 (2016)

AS P. Antonelli, S. Spirito, Global existence of finite energy weak solutions of quantum Navier-Stokes equations, Arch. Ration. Mech. Anal. 225, 1161-1199 (2017)

Coppel
W. A. Coppel, Dichotomies in Stability Theory, Springer-Verlag Berlin Heidelberg, 1978

Michele1
F. Di Michele, P. Marcati, B. Rubino, Steady states and interface transmission conditions for heterogeneous quantum-classical 1-D hydrodynamic model of semiconductor devices. Physica D: Nonlinear Phenomena, 243(1), pp. 1-13, 2013 

Michele
F. Di Michele, P. Marcati, B. Rubino, Stationary solution for transient quantum hydrodynamics with bohmenian-type boundary conditions, Computational and Applied Mathematics, 36(1), pp. 459-479, 2017

DFMD. Donatelli, E. Feireisl, P. Marcati, Well/ill posedness for the Euler- Korteweg-Poisson system and related problems, Comm. Partial Differential Equations, 40, 1314-1335 (2015)

DM1D. Donatelli, P. Marcati, Quasineutral limit, dispersion and oscillations for Korteweg type fluids, SIAM J. Math. Anal. 47, 2265-2282  (2015)

DM D. Donatelli, P. Marcati,  Low Mach number limit for the quantum hydrodynamics system, Res. Math. Sci. 3, 3-13  (2016)

GLT J. Giesselmann, C. Lattanzio, and  A.E. Tzavaras,  Relative Energy for the Korteweg Theory and Related Hamiltonian Flows in Gas Dynamics, Arch. Ration. Mech. Anal. 223 , 1427-1484  (2017)





Gurevich1A. V. Gurevich and A. P. Meshcherkin. Expanding self-similar discontinuities and shock waves in dispersive hydrodynamics, Sov. Phys. JETP, 60(4):732-740, 1984.

Gurevich
A. V. Gurevich and L. P. Pitaevskii, Nonstationary structure of a collisionless shock wave, Sov. Phys. JETP, 38:291-297 (1974)

Hoefer1M. A. Hoefer, Shock Waves in Dispersive Eulerian Fluids, J. Nonlinear Sci., Volume 24, Issue 3, pp 525-577, (2014)

HoeferM. A. Hoefer, M. J. Ablowitz, I. Coddington, E. A. Cornell, P. Engels, and V. Schweikhard, Dispersive and classical shock waves in Bose-Einstein condensates and gas dynamics, Phys. Rev. A, 74:023623 (2006)

HumpherysJ. Humpherys, On the shock wave spectrum for isentropic gas dynamics with capillarity, J. Differential Equations, 246(7):2938-2957 (2009)



Kuznetsov 
Y. Kuznetsov, Elements of applied bifurcation theory, second edition, Springer-Verlag New York, 1998.




Kato
T. Kato, Perturbation Theory for Linear Operators, second edition, Springer-Verlag Berlin Heidelberg, 1995

LMZ18_num
C. Lattanzio, P. Marcati, D. Zhelyazov, Dispersive shocks and spectral analysis for linearized Quantum Hydrodynamics, arXiv preprint arXiv:1904.09885  (2019).

NovS. Novikov, S. V. Manakov, L. P. Pitaevskii, and V. E. Zakharov, Theory of Solitons, Consultants Bureau, New York, 1984.

Perko
L. Perko, Differential Equations and Dynamical Systems, third edition, Springer-Verlag New York, 2001

SagdeevS.R. Z. Sagdeev,  Kollektivnye protsessy i udarnye volny v razrezhennol plazme (Collective processes and shock waves in a tenuous plasma), in: Voprosy teorii plazmy (Problems of Plasma Theory), Vol. 5, Atomizdat, 1964.

Sandstede
B. Sandstede, Stability of Travelling Waves, Handbook of Dynamical Systems II, Elsevier (2002) 983-1055

Zak V. E. Zakharov, Stability Of Periodic Waves Of Finite Amplitude On The Surface Of A Deep Fluid, Zhurnal Prildadnoi Mekhaniki i Tekhnicheskoi Fiziki, 9(2), 86-94 (1968)


